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

高普陽(yáng)

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

高普陽(yáng)

（長(zhǎng)安大學(xué) 理學(xué)院，陜西 西安 710064）

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

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

0　引言

非牛頓流體充填過(guò)程中涉及兩種不同類型流體以及液體之間的運(yùn)動(dòng)界面，在三維情形下進(jìn)行數(shù)值模擬，一直是關(guān)注的熱點(diǎn)和難點(diǎn)［1-7］。數(shù)值模擬的關(guān)鍵是精確求解由兩種流體構(gòu)成流場(chǎng)信息、準(zhǔn)確捕捉運(yùn)動(dòng)界面以及保證流體質(zhì)量守恒性。

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

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

1　數(shù)學(xué)模型

非牛頓流體充填過(guò)程屬于非牛頓-牛頓兩相流動(dòng)問(wèn)題，其數(shù)學(xué)模型主要包括流場(chǎng)控制方程和界面演化方程。

1.1　水平集方程

用守恒型水平集方法追蹤兩相流動(dòng)的自由界面。用水平集函數(shù)（）的零等值面表示兩流體間的運(yùn)動(dòng)界面，流場(chǎng)中某個(gè)點(diǎn)到界面的距離用函數(shù)（）的絕對(duì)值表示：

（）在空氣區(qū)域內(nèi)取正值，在液體區(qū)域內(nèi)取負(fù)值。

在速度場(chǎng)影響下，自由界面的形態(tài)和位置發(fā)生變化，運(yùn)動(dòng)界面的水平集函數(shù)為

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

其中，（）為光滑的Heaviside函數(shù)，

1.2　Navier-Stokes方程

非牛頓流體充填過(guò)程屬于非牛頓-牛頓兩相流動(dòng)問(wèn)題，流動(dòng)過(guò)程涉及兩種不同特性的流體。為方便計(jì)算，用統(tǒng)一形式：

其中，=（，，）表示三維速度矢量，為壓力。和分別為流動(dòng)區(qū)域內(nèi)的密度和黏度，

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

2　數(shù)值算法

2.1　統(tǒng)一形式Navier-Stokes方程的求解

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

由中間速度*計(jì)算下一時(shí)刻的壓力，半離散格式為

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

在進(jìn)行空間離散前，需對(duì)求解區(qū)域Ω進(jìn)行三角剖分，得到Ω，用Ω（=1，2，…，）表示三角網(wǎng)格上的小三角形。為保證計(jì)算過(guò)程穩(wěn)定，需分別對(duì)速度和壓力采用二次和一次多項(xiàng)式插值。定義2個(gè)有限元空間：

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

2.2　水平集方程的求解

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

定義間斷有限元空間

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

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

2.3　算法流程

算法流程如圖1所示。

圖1　算法流程

3　數(shù)值算例

3.1　三維圓球剪切流動(dòng)

主要研究三維情形下的圓球剪切流動(dòng)［18-19］，以驗(yàn)證數(shù)值算法在處理大變形自由界面問(wèn)題時(shí)的穩(wěn)定性、準(zhǔn)確性和質(zhì)量守恒性。假設(shè)初始狀態(tài)下流場(chǎng)中圓球的球心坐標(biāo)為（0.35，0.35，0.35），半徑為0.15，如圖2所示。流場(chǎng)內(nèi)的速度為

圖3　當(dāng)t =3.0時(shí)不同網(wǎng)格的自由界面形態(tài)

圖4　不同時(shí)刻的自由界面形態(tài)

圖5　網(wǎng)格3上相對(duì)質(zhì)量誤差隨時(shí)間的變化

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

圖6　三維型腔示意

圖7　計(jì)算網(wǎng)格3示意

圖8　當(dāng)t =1.9 s時(shí)不同網(wǎng)格的熔體前沿界面形態(tài)

圖9給出了不同時(shí)刻非牛頓流體前沿界面形態(tài)的實(shí)驗(yàn)結(jié)果、文獻(xiàn)［9］的數(shù)值結(jié)果和本文的數(shù)值結(jié)果?？芍跏茧A段液體前沿界面呈圓弧狀，隨著時(shí)間的推進(jìn)，液體前沿界面的弧度逐漸變小，本文方法的數(shù)值結(jié)果和實(shí)驗(yàn)結(jié)果［9］吻合較好。圖9中，（a）1.1 s，（b）1.9 s，（c）2.7 s，（d）3.5 s，（e）4.3 s分別對(duì)應(yīng)液體區(qū)域所占型腔體積20%，40%，60%，80%和100%時(shí)的情形。網(wǎng)格1、網(wǎng)格2和網(wǎng)格3計(jì)算得到的熔體質(zhì)量的最大相對(duì)誤差分別為2.20%，1.00%和0.67%。數(shù)值結(jié)果表明，耦合算法能較穩(wěn)定、準(zhǔn)確地求解三維非牛頓流體的充填過(guò)程，并保證非牛頓流體具有較好的質(zhì)量守恒性。隨著充填過(guò)程的進(jìn)行，型腔內(nèi)熔體的總質(zhì)量逐漸增加，圖10給出了網(wǎng)格3熔體質(zhì)量數(shù)值結(jié)果和精確結(jié)果之間的相對(duì)誤差隨時(shí)間的變化。圖11進(jìn)一步展示了不同時(shí)刻前沿界面的三維視圖。

圖9　不同時(shí)刻的實(shí)驗(yàn)結(jié)果（左）、文獻(xiàn)［9］的數(shù)值結(jié)果（中）和本文的數(shù)值結(jié)果（右）

圖10　網(wǎng)格3熔體質(zhì)量相對(duì)誤差隨時(shí)間的變化

圖11　不同時(shí)刻前沿界面的三維視圖

4　結(jié)論

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

［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］ 周文，歐陽(yáng)潔，楊斌鑫，等. 三維非等溫非牛頓流體充模過(guò)程的建模與模擬［J］. 化工學(xué)報(bào)， 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.

國(guó)家自然科學(xué)基金資助項(xiàng)目（11901051，11971075）；陜西省自然科學(xué)基礎(chǔ)研究計(jì)劃青年項(xiàng)目（2020JQ-338）；長(zhǎng)安大學(xué)中央高?；究蒲袠I(yè)務(wù)費(fèi)項(xiàng)目（300102122107）；陜西省科學(xué)技術(shù)協(xié)會(huì)青年人才托舉計(jì)劃項(xiàng)目（20220504）.

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