三維非牛頓流體充填過程的有限元-間斷有限元數值模擬研究

高普陽

三維非牛頓流體充填過程的有限元-間斷有限元數值模擬研究

高普陽

（長安大學 理學院，陜西 西安 710064）

針對三維非牛頓流體充填問題，建立了有限元-間斷有限元耦合算法。對于兩相Navier-Stokes方程，基于壓力增量修正格式分三步求解，分別采用二次和一次拉格朗日插值多項式求解速度和壓力，以確保計算過程穩定。采用守恒型水平集（level set）方法追蹤運動界面，并依據間斷有限元方法求解水平集和重新初始化方程。以三維圓球剪切流動及非牛頓流體三維平板型腔充填過程為例，并與已有文獻的數值和實驗結果進行比較，以驗證數值算法的穩定性、準確性以及流體的質量守恒性。

有限元；間斷有限元；水平集；非牛頓充填過程

0　引言

非牛頓流體充填過程中涉及兩種不同類型流體以及液體之間的運動界面，在三維情形下進行數值模擬，一直是關注的熱點和難點［1-7］。數值模擬的關鍵是精確求解由兩種流體構成流場信息、準確捕捉運動界面以及保證流體質量守恒性。

TOMé等［8］采用有限差分法求解流場控制方程、用VOF法追蹤運動界面，并基于CFD Freeflow3D程序數值模擬了三維容器內非牛頓流體的充填過程，研究了入口速度對充填過程中自由界面形態及裹氣現象的影響，并將數值結果與實驗結果進行了對比。MUKRAS等［9］用ANSYS-CFX數值模擬了三維型腔內高密度聚乙烯的充填過程，分析了不同注射速度條件下的充填模式。LIU等［10-11］數值模擬了三維型腔內黏性非牛頓流體的共注射成型過程，對一些特有的流動現象進行了研究。ZHUANG等［12］基于有限體積法，用水平集方法追蹤運動界面，數值模擬了三維型腔的充填過程，并分析了溫度、黏彈性應力等對流動性態的影響。HE等［13］用光滑粒子動力學方法（SPH）模擬了冪律型非牛頓流體的充填過程，并分析了流動過程中前沿界面的形態以及纖維的取向。BORZENKO等［14］用SIMPLE算法求解流場控制方程，用VOF方法追蹤運動界面，研究了三維矩形腔的非牛頓流體充填過程，分析了數和數等參數對自由界面的影響機理。目前尚鮮見相關用有限元-間斷有限元方法模擬三維充填過程研究的報道，且充填過程大多未考慮流體的質量守恒性。

本文構建了三維充填過程的統一計算模型，并采用有限元-間斷有限元耦合算法對模型進行數值模擬。首先，采用壓力增量修正算法對統一計算模型進行分裂；然后，基于有限元法對其進行數值求解，對速度和壓力分別取二次和一次拉格朗日插值基函數，以保證計算過程穩定；最后，為準確追蹤運動界面并保證流體較好的質量守恒性，用守恒型水平集方法追蹤運動界面，并用間斷有限元法數值求解水平集及重新初始化方程。

1　數學模型

非牛頓流體充填過程屬于非牛頓-牛頓兩相流動問題，其數學模型主要包括流場控制方程和界面演化方程。

1.1　水平集方程

用守恒型水平集方法追蹤兩相流動的自由界面。用水平集函數（）的零等值面表示兩流體間的運動界面，流場中某個點到界面的距離用函數（）的絕對值表示：

（）在空氣區域內取正值，在液體區域內取負值。

在速度場影響下，自由界面的形態和位置發生變化，運動界面的水平集函數為

為確保界面附近水平集函數的符號距離特性和流體的質量守恒性，引入了修正型水平集，重新初始化方程［15］：

其中，（）為光滑的Heaviside函數，

1.2　Navier-Stokes方程

非牛頓流體充填過程屬于非牛頓-牛頓兩相流動問題，流動過程涉及兩種不同特性的流體。為方便計算，用統一形式：

其中，=（，，）表示三維速度矢量，為壓力。和分別為流動區域內的密度和黏度，

不同于牛頓流體，非牛頓流體的黏度受形變速率影響。本文用冪律型本構方程［16］描述黏度與形變速率張量之間的非線性關系：

2　數值算法

2.1　統一形式Navier-Stokes方程的求解

用基于有限元的壓力增量修正格式［17］對統一形式的Navier-Stokes方程進行數值求解。首先，引入中間速度*，對式（8）左端的時間項進行空間離散，顯式處理壓力項，對非線性對流項做線性化處理，式（8）右端的速度選用中間時刻的值：

由中間速度*計算下一時刻的壓力，半離散格式為

對式（10）和式（11）兩邊取散度，由速度+1的不可壓縮條件，可得

在進行空間離散前，需對求解區域Ω進行三角剖分，得到Ω，用Ω（=1，2，…，）表示三角網格上的小三角形。為保證計算過程穩定，需分別對速度和壓力采用二次和一次多項式插值。定義2個有限元空間：

式（10）～式（12）的空間離散格式分別為

2.2　水平集方程的求解

在得到速度場之后，需要進一步求解水平集方程，以更新運動界面。水平集方程為雙曲型方程，本文采用間斷有限元方法進行數值求解。首先，將水平集方程（式（1））在時間上進行隱式離散，得到半離散方程

定義間斷有限元空間

因此，式（16）的全離散格式為

同理，水平集重新初始化方程的全離散格式為

2.3　算法流程

算法流程如圖1所示。

圖1　算法流程

3　數值算例

3.1　三維圓球剪切流動

主要研究三維情形下的圓球剪切流動［18-19］，以驗證數值算法在處理大變形自由界面問題時的穩定性、準確性和質量守恒性。假設初始狀態下流場中圓球的球心坐標為（0.35，0.35，0.35），半徑為0.15，如圖2所示。流場內的速度為

圖3　當t =3.0時不同網格的自由界面形態

圖4　不同時刻的自由界面形態

圖5　網格3上相對質量誤差隨時間的變化

3.2　非牛頓流體三維平板型腔充填過程

圖6　三維型腔示意

圖7　計算網格3示意

圖8　當t =1.9 s時不同網格的熔體前沿界面形態

圖9給出了不同時刻非牛頓流體前沿界面形態的實驗結果、文獻［9］的數值結果和本文的數值結果。可知，初始階段液體前沿界面呈圓弧狀，隨著時間的推進，液體前沿界面的弧度逐漸變小，本文方法的數值結果和實驗結果［9］吻合較好。圖9中，（a）1.1 s，（b）1.9 s，（c）2.7 s，（d）3.5 s，（e）4.3 s分別對應液體區域所占型腔體積20%，40%，60%，80%和100%時的情形。網格1、網格2和網格3計算得到的熔體質量的最大相對誤差分別為2.20%，1.00%和0.67%。數值結果表明，耦合算法能較穩定、準確地求解三維非牛頓流體的充填過程，并保證非牛頓流體具有較好的質量守恒性。隨著充填過程的進行，型腔內熔體的總質量逐漸增加，圖10給出了網格3熔體質量數值結果和精確結果之間的相對誤差隨時間的變化。圖11進一步展示了不同時刻前沿界面的三維視圖。

圖9　不同時刻的實驗結果（左）、文獻［9］的數值結果（中）和本文的數值結果（右）

圖10　網格3熔體質量相對誤差隨時間的變化

圖11　不同時刻前沿界面的三維視圖

4　結論

提出了三維情形下的有限元-間斷有限元耦合算法，并基于非牛頓-牛頓兩相流計算模型，數值模擬了三維非牛頓流體的充填過程。在三維圓球剪切流動的數值算例中，圓球形狀演化趨勢和文獻［18-19］的描述一致，當自由界面發生最大變形時亦未出現任何不穩定現象，在充填過程中，網格3熔體質量最大相對誤差為1.5%。在三維矩形型腔非牛頓流體的充填過程中，前沿界面非常穩定，在初始階段，前沿界面呈圓弧狀，隨著時間的進展，逐漸變平緩，前沿界面的演化趨勢與實驗結果及文獻［9］的結果相吻合，在充填過程中，網格3非牛頓流體質量的最大相對誤差為0.67%。數值結果表明，本文的有限元-間斷有限元耦合算法能穩定、有效、準確地模擬三維非牛頓流體的充填過程，并保證了流體具有較好的質量守恒性。

［1］ XU X Y， OUYANG J， YANG B X， et al. SPH simulations of three-dimensional non-Newtonian free surface flows［J］. Computer Methods in Applied Mechanics and Engineering， 2013， 256： 101-116. DOI：10.1016/j.cma.2012.12.017

［2］ FIGUEIREDO R A， OISHI C M， CUMINATO J A， et al. Three-dimensional transient complex free surface flows： Numerical simulation of XPP fluid［J］. Journal of Non-Newtonian Fluid Mechanics， 2013， 195： 88-98. DOI：10.1016/j.jnnfm.2013.01.004

［3］ ZHANG C H， CHEN S G， SUN Q C， et al. Free-surface simulations of Newtonian and non-Newtonian fluids with the lattice Boltzmann method［J］. Acta Geologica Sinica（English Edition）， 2016， 90（3）： 999-1010. DOI：10.1111/1755-6724.12740

［4］ 周文，歐陽潔，楊斌鑫，等. 三維非等溫非牛頓流體充模過程的建模與模擬［J］. 化工學報， 2011， 62（3）： 618-627.

ZHOU W， OUYANG J， YANG B X， et al. Modeling and simulation of 3D mold filling process for non-isothermal non-Newtonian fluid［J］. CIESC Journal， 2011， 62（3）： 618-627.

［5］ BAUM M， ANDERS D. A numerical simulation study of mold filling in the injection molding process［J］. Computer Methods in Materials Science， 2021， 21（1）： 25-34. DOI：10.1002/fld.166

［6］ LI Q， QU F C. A level set based immersed boundary method for simulation of non-isothermal viscoelastic melt filling process［J］. Chinese Journal of Chemical Engineering， 2021， 32： 119-133. DOI：10.1016/j.cjche.2020.09.057

［7］ XU X Y， TIAN L Y， PENG S， et al. Development of SPH for simulation of non-isothermal viscoelastic free surface flows with application to injection molding［J］. Applied Mathematical Modelling， 2022， 104： 782-805. DOI：10.1016/j.apm.2021.12.015

［8］ TOMé M F， CASTELO A，NóBREGA J M， et al. Numerical and experimental investigations of three-dimensional container filling with Newtonian viscous fluids［J］. Computers & Fluids， 2014， 90： 172-185. DOI：10.1016/j.compfluid.2013.11.015

［9］ MUKRAS S M S， AL-MUFADI F A. Simulation of HDPE mold filling in the injection molding process with comparison to experiments［J］. Arabian Journal for Science and Engineering， 2016， 41（5）： 1847-1856. DOI：10.1007/s13369-015-1970-9

［10］LIU Q S， OUYANG J， ZHOU W， et al. Numerical simulation of the sequential coinjection molding process based on level set method［J］. Polymer Engineering & Science， 2015， 55（8）： 1707-1719. DOI：10.1002/pen.24009

［11］LIU Q S， OUYANG J， LI W M， et al. Three-dimensional modelling and simulation of sequential co-injection moulding with application to a toothbrush handle［J］. The Canadian Journal of Chemical Engineering， 2016， 94（2）： 382-390. DOI：10.1002/cjce.22394

［12］ZHUANG X， OUYANG J， LI W M， et al. Three-dimensional simulations of non-isothermal transient flow and flow-induced stresses during the viscoelastic fluid filling process［J］. International Journal of Heat and Mass Transfer， 2017， 104： 374-391. DOI：10. 1016/j.ijheatmasstransfer.2016.08.064

［13］HE L P， LU G， CHEN D C， et al. Three-dimensional smoothed particle hydrodynamics simulation for injection molding flow of short fiber-reinforced polymer composites［J］. Modelling and Simulation in Materials Science and Engineering， 2017， 25（5）： 055007.

［14］BORZENKO E I， HEGAJ E I. Three-dimensional simulation of a tank filling with a viscous fluid using the VOF method［J］. Journal of Siberian Federal University （Mathematics & Physics）， 2020， 13（6）： 670-677. DOI：10.17516/1997-1397-2020-13-6-670-677

［15］GU Z H， WEN H L， YU C H， et al. Interface-preserving level set method for simulating dam-break flows［J］. Journal of Computational Physics， 2018， 374： 249-280. DOI：10.1016/j.jcp.2018.07.057

［16］AGUIRRE A， CASTILLO E， CRUCHAGA M， et al. Stationary and time-dependent numerical approximation of the lid-driven cavity problem for power-law fluid flows at high Reynolds numbers using a stabilized finite element formulation of the VMS type［J］. Journal of Non-Newtonian Fluid Mechanics， 2018， 257： 22-43. DOI：10.1016/j.jnnfm.2018.03.014

［17］GUERMOND J L， MINEV P， SHEN J. An overview of projection methods for incompressible flows［J］. Computer Methods in Applied Mechanics and Engineering， 2006， 195（44-47）： 6011-6045. DOI：10.1016/j.cma.2005.10.010

［18］GUERMOND J L， de LUNA M Q， THOMPSON T. An conservative anti-diffusion technique for the level set method［J］. Journal of Computational and Applied Mathematics， 2017， 321： 448-468. DOI：10.1016/j.cam.2017.02.016

［19］LUO K， SHAO C X， YANG Y， et al. A mass conserving level set method for detailed numerical simulation of liquid atomization［J］. Journal of Computational Physics， 2015， 298： 495-519. DOI：10.1016/j.jcp.2015.06.009

The numerical investigation of non-Newtonian fluid filling process via finite element and discontinuous Galerkin method

GAO Puyang

（，，710064，）

In this paper, we develop a coupled finite element and discontinuous Galerkin method in three dimension and study the non-Newtonian fluid filling process. To solve two phase Navier-Stokes equations, we employ the incremental pressure correction scheme to accomplish it in three steps. In order to guarantee the computational stability, we take the second order and first order interpolation polynomials for the velocity and pressure, respectively. In addition, the conservative Level Set method is employed to capture the moving interface. The discontinuous Galerkin method is used to solve the Level Set and its re-initialization equations. We take the three dimensional vortex shearing problem and the three dimensional non-Newtonian fluid filling process to verify the proposed approach, compare the result with the numerical results and existing experimental data to illustrate the stability, accuracy and the mass conservation property of the coupled scheme.

finite element; discontinuous Galerkin; level set; non-Newtonian filling process

O 242.1

A

1008?9497（2023）01?049?07

2021?12?13.

國家自然科學基金資助項目（11901051，11971075）；陜西省自然科學基礎研究計劃青年項目（2020JQ-338）；長安大學中央高校基本科研業務費項目（300102122107）；陜西省科學技術協會青年人才托舉計劃項目（20220504）.

高普陽（1991—），ORCID：https：//orcid.ord/0000-0001-8620-1783，男，博士，講師，主要從事非牛頓流動問題的數值算法研究，E-mail： gaopuyang@chd.edu.cn.