粘彈流動(dòng)的間斷有限元模擬?

郭虹平,歐陽(yáng)潔,楊廣輝,周 文

(西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，西安 710129)

1 引言

近年來(lái)，隨著聚合物工業(yè)的飛速發(fā)展，粘彈流體流動(dòng)的數(shù)值研究備受關(guān)注．粘彈流體流動(dòng)的控制方程由質(zhì)量方程、動(dòng)量方程和本構(gòu)方程構(gòu)成，屬于橢圓-雙曲混合型偏微分方程組．?dāng)?shù)值模擬算法的不準(zhǔn)確可能會(huì)導(dǎo)致應(yīng)力張量正定性的改變，盡管此變化對(duì)橢圓型方程的求解影響不大，但對(duì)雙曲型方程的求解會(huì)帶來(lái)很大誤差[1]．因此，發(fā)展適合求解橢圓-雙曲混合型偏微分方程組的高精度方法很有必要．

目前模擬粘彈流動(dòng)的主要方法是有限差分法、有限體積法和有限元法．有限差分法求解規(guī)則區(qū)域的問(wèn)題時(shí)，算法實(shí)施容易，但對(duì)復(fù)雜區(qū)域的適應(yīng)性很差．有限體積法有很好的幾何靈活性，逼近及格式都是局部的，但不能在一般非結(jié)構(gòu)網(wǎng)格上取得高精度．傳統(tǒng)有限元法對(duì)不規(guī)則區(qū)域的適應(yīng)性較好、精度較高，但不適于處理間斷問(wèn)題，對(duì)網(wǎng)格加密或減疏處理時(shí)需考慮連續(xù)性的限制．間斷有限元(discontinuous Galerkin，簡(jiǎn)稱(chēng)DG)法兼顧有限體積法及有限元法的優(yōu)點(diǎn)，采用局部逼近，精度高且高度并行，易于處理非規(guī)則區(qū)域[2]．

Yurun[3]分別采用加入穩(wěn)定化方案的傳統(tǒng)有限元法和DG方法求解穩(wěn)態(tài)UCM本構(gòu)模型，驗(yàn)證了兩種方法均可有效求解光滑問(wèn)題．但對(duì)含應(yīng)力奇異點(diǎn)的平板收縮流模擬時(shí)，DG方法表現(xiàn)出更好的魯棒性．2005年，Kim等人[4]模擬了基于Oldroyd-B本構(gòu)模型描述的粘彈流動(dòng)，實(shí)現(xiàn)了質(zhì)量方程、動(dòng)量方程和本構(gòu)方程的耦合求解．其質(zhì)量方程和動(dòng)量方程的求解采用傳統(tǒng)有限元法，本構(gòu)方程的求解采用DG方法，且用并行算法提高了計(jì)算效率．但是，運(yùn)用非結(jié)構(gòu)網(wǎng)格上統(tǒng)一的DG方法對(duì)粘彈流動(dòng)進(jìn)行模擬研究卻鮮有報(bào)道．

鑒于Oldroyd-B本構(gòu)模型是檢驗(yàn)數(shù)值算法優(yōu)劣的最苛刻的粘彈本構(gòu)模型，而發(fā)展基于非結(jié)構(gòu)網(wǎng)格的數(shù)值方法是模擬復(fù)雜區(qū)域上實(shí)際問(wèn)題最有效的技術(shù)手段，本文基于非結(jié)構(gòu)網(wǎng)格，建立了Oldroyd-B粘彈流動(dòng)的統(tǒng)一DG求解框架．該框架包含采用IPDG(interior penalty DG)求解質(zhì)量方程和動(dòng)量方程，RKDG(Runge-Kutta DG)求解本構(gòu)方程這兩個(gè)核心．對(duì)平面Poiseuille流和4:1平板收縮流問(wèn)題的數(shù)值模擬，驗(yàn)證了本文方法的有效性．

2 粘彈流動(dòng)的控制方程

考慮由Oldroyd-B[5,6]本構(gòu)模型描述的瞬態(tài)、不可壓、等溫流動(dòng)問(wèn)題，其流場(chǎng)控制方程如下：

質(zhì)量方程

動(dòng)量方程

本構(gòu)方程

其中u為流體速度，p為壓力，τ為應(yīng)力張量；表示牛頓粘度ηs和總粘度η0的比；，ρ為密度，U為特征速度，L為特征長(zhǎng)度；，λ為特征松弛時(shí)間；為τ的上隨體導(dǎo)數(shù)，即

為便于構(gòu)造求解格式，化簡(jiǎn)(3)式，得

其中

3 本構(gòu)方程的RKDG方法求解

Cockburn和Shu提出了求解雙曲型方程的RKDG方法[7-9]．Hesthaven和Warburton[10]則通過(guò)巧妙合理地布置節(jié)點(diǎn)來(lái)構(gòu)造基函數(shù)，簡(jiǎn)化了DG方法在非結(jié)構(gòu)網(wǎng)格上的實(shí)施．本文將基于非結(jié)構(gòu)網(wǎng)格，用RKDG方法單獨(dú)求解Oldroyd-B本構(gòu)方程，以獲得新的時(shí)間層的應(yīng)力張量．

為構(gòu)造本構(gòu)方程的RKDG方法，結(jié)合(1)式和(5)式可得

令f(τ)=uτ，則不可壓縮流場(chǎng)中Oldroyd-B本構(gòu)方程的守恒形式可表示成

為求解(7)式，對(duì)區(qū)域?構(gòu)造三角剖分Th，其中剖分單元用?k表示，并定義有限元空間

Vh(?k)是單元?k上的基函數(shù)，其中np為基函數(shù)的個(gè)數(shù)．本文取Lagrange型基函數(shù)lki(x)，其系數(shù)便是相應(yīng)節(jié)點(diǎn)處物理量的值．

設(shè)單元?k上的逼近函數(shù)形式為

對(duì)(7)式兩邊乘以試驗(yàn)函數(shù)，且進(jìn)行一次分部積分，得DG方法的弱形式為

其中f?=uτ?為數(shù)值通量[11]，n為單元?k邊界上的外法向向量．為保證數(shù)值格式的穩(wěn)定性，τ?選用迎風(fēng)型數(shù)值通量，見(jiàn)圖1，即

圖1:二維迎風(fēng)通量示意圖

由(8)式知，在單元?k上

將(8),(10),(11)式代入(9)式，可得

定義如下矩陣：質(zhì)量矩陣

剛度矩陣

邊界矩陣

化簡(jiǎn)(12)式，可得如下常微分方程

其中

采用TVD Runge-Kutta方法對(duì)(13)式進(jìn)行顯式時(shí)間離散[11]即為RKDG方法，本文采用二階RK格式．

4 質(zhì)量方程和動(dòng)量方程的IPDG求解

4.1節(jié)和4.2節(jié)將分別給出求解質(zhì)量方程、動(dòng)量方程的時(shí)間分裂格式和空間離散格式(即IPDG方法)．其中動(dòng)量方程(2)中的應(yīng)力作為擬體力處理，且在動(dòng)量方程、質(zhì)量方程的求解過(guò)程中不斷地更新速度和壓力．

4.1 時(shí)間分裂格式

為求解動(dòng)量方程(2)，采用剛性穩(wěn)定的時(shí)間步長(zhǎng)法[12]，將第n層到第n+1層每個(gè)時(shí)間步的推進(jìn)分為三個(gè)階段．

步驟1非線(xiàn)性對(duì)流方程

步驟2壓力方程

步驟3粘性方程

其中n=0時(shí)，計(jì)算參數(shù)取為

以后的參數(shù)為

由分裂過(guò)程可知，第1步通過(guò)求解對(duì)流方程(14)式來(lái)計(jì)算預(yù)估速度；第2步引入中間速度完成壓力投影，通過(guò)求解壓力Poisson方程(16)式得到壓力；第3步通過(guò)(15)式更新中間速度，最后隱式求解粘性方程(17)式得到un+1．

4.2 空間離散

下面給出分裂方程(14),(16),(17)的空間離散格式．其中非線(xiàn)性對(duì)流方程(14)的DG方法求解中，選用合適的數(shù)值通量以提高數(shù)值解的穩(wěn)定性；橢圓方程(16)和(17)的DG方法求解中，選用單元計(jì)算模板較小的IPDG方法[13]．

為下文推導(dǎo)方便，定義為u的跳躍，為u的平均．選用Lax-Friedrichs通量對(duì)方程(14)的右端項(xiàng)進(jìn)行空間離散，得

橢圓型方程(16)和(17)分別為Poisson和Helmholtz方程，其統(tǒng)一形式為

對(duì)于Possion方程：；對(duì)于Helmholtz方程：

選用IPDG方法可得上述橢圓方程的離散結(jié)果為

其中??D和??N分別為Dirich let和Neumann邊界條件，Γ為Γh中的所有單元邊界的集合，κj為懲罰因子，與網(wǎng)格尺寸同階．

5 統(tǒng)一的DG算法

統(tǒng)一算法的宗旨是將屬于橢圓-雙曲混合型偏微分方程組的粘彈流動(dòng)控制方程采用相同的數(shù)值方法進(jìn)行統(tǒng)一求解，從而使方程之間的數(shù)據(jù)信息傳遞簡(jiǎn)單而且快捷，使其程序?qū)崿F(xiàn)統(tǒng)一化．本文建立的統(tǒng)一DG框架包含兩個(gè)核心：一是采用IPDG方法結(jié)合Asam s-Bashforth二階時(shí)間離散格式求解質(zhì)量方程和動(dòng)量方程；二是采用RKDG方法求解本構(gòu)方程．這三套方程的耦合即構(gòu)成了本文粘彈流動(dòng)模擬的統(tǒng)一DG方法．其算法流程如下：

步驟1初始化速度u0、壓力p0和應(yīng)力τ0；

步驟2對(duì)于i=0,…,N(t0→tN):

1)將應(yīng)力τi作為擬體力，引入中間速度,，運(yùn)用第4節(jié)的離散格式求解質(zhì)量方程和動(dòng)量方程，更新速度、壓力；

2)根據(jù)更新的速度ui+1，更新si+1；

3)運(yùn)用第3節(jié)的RKDG方法求解本構(gòu)方程，更新應(yīng)力場(chǎng)；

步驟3輸出數(shù)據(jù)．

6 數(shù)值算例

6.1 平面Poiseuille流

當(dāng)Re很小時(shí)，平面Poiseuille流具有解析解，可以定量檢驗(yàn)數(shù)值算法的穩(wěn)定性和精度．圖2給出了平面Poiseuille流的流動(dòng)示意．流場(chǎng)在出口處充分發(fā)展，固壁處采用無(wú)滑移邊界條件．入口速度設(shè)為u=0.4y(1?y),v=0．

圖2:平面Poiseuille流動(dòng)示意

Oldroyd-B模型的穩(wěn)態(tài)Poiseuille流動(dòng)的精確解為[14]：

圖3給出了不同We數(shù)下第一法向應(yīng)力和剪切應(yīng)力在x=4時(shí)的數(shù)值解與解析解比較．圖3(a)中各曲線(xiàn)從上到下依次表示Re=0.01,,We=1,We=2,We=3,We=4,We=5,We=6時(shí)第一法向應(yīng)力的數(shù)值解和解析解，圖3(b)為We=2時(shí)的剪切應(yīng)力模擬結(jié)果．可以看出，第一法向應(yīng)力和剪切應(yīng)力的數(shù)值解和解析解吻合很好，并且在We=6時(shí)仍可得到滿(mǎn)意結(jié)果．表1給出了應(yīng)力τxx在l2范數(shù)下的誤差分析．由表1可知該方法具較高的計(jì)算精度，且隨著網(wǎng)格的加密和插值次數(shù)的提高，模擬結(jié)果更加精確．

圖3: 平面Poiseuille流數(shù)值解與解析解的比較

表1: 應(yīng)力τxx在l2范數(shù)下的誤差分析

6.2 平板收縮流

4:1平板收縮流具有旋轉(zhuǎn)、拉伸和強(qiáng)剪切等復(fù)雜的流動(dòng)特征，在收縮口處存在應(yīng)力奇點(diǎn)，其數(shù)值解存在大梯度變化．對(duì)4:1粘彈收縮流動(dòng)問(wèn)題的模擬，可以進(jìn)一步驗(yàn)證算法對(duì)模擬非線(xiàn)性復(fù)雜粘彈流動(dòng)的有效性．

圖4給出了4:1平板收縮流的計(jì)算區(qū)域、邊界條件及網(wǎng)格剖分．計(jì)算時(shí)固壁處速度采用無(wú)滑移邊界條件，入口處設(shè)為

圖4:平板收縮流示意

為檢驗(yàn)本構(gòu)模型求解的正確性，首先考察忽略非線(xiàn)性項(xiàng)影響即(N(un)=0)的蠕動(dòng)流動(dòng)．圖5給出了Re=1時(shí)，在不同We數(shù)下收縮口附近蠕動(dòng)流的流線(xiàn)分布．從圖中可以看出隨著We數(shù)的增大，流場(chǎng)中的渦流區(qū)域逐漸減小，這與文獻(xiàn)[15]的數(shù)值結(jié)果一致．

圖5: 蠕動(dòng)流在不同We數(shù)下的流線(xiàn)比較

為真正模擬Oldroyd-B本構(gòu)模型的粘彈流體流動(dòng)，其次考慮具有非線(xiàn)性項(xiàng)影響的粘彈流動(dòng)．圖6至圖9給出了在Re=1,We=1時(shí)收縮口附近粘彈流動(dòng)的流線(xiàn)分布和應(yīng)力分布．從圖中可以看出，收縮口附近形成強(qiáng)拉伸區(qū)域，在上游的拐角處，流動(dòng)以旋轉(zhuǎn)為主，形成主渦，這與文獻(xiàn)[1,16]預(yù)測(cè)的流動(dòng)趨勢(shì)完全吻合．

圖6:收縮口附近的流線(xiàn)分布

圖7: 收縮口附近的應(yīng)力τxx分布

圖8: 收縮口附近應(yīng)力τxy分布

圖9: 收縮口附近的應(yīng)力τyy分布

7 結(jié)論

本文建立了一種求解Oldroyd-B粘彈流體流動(dòng)的數(shù)值方法，給出了粘彈流動(dòng)求解的統(tǒng)一DG框架．該方法運(yùn)用IPDG求解質(zhì)量方程和動(dòng)量方程，運(yùn)用RKDG求解本構(gòu)方程，成功實(shí)現(xiàn)了粘彈流動(dòng)流場(chǎng)控制方程的DG方法耦合求解．對(duì)平面Poiseuille流和4:1平板收縮流問(wèn)題的數(shù)值模擬結(jié)果表明：

1)本文方法在求解Oldroyd-B本構(gòu)模型時(shí)無(wú)需加入穩(wěn)定化方案，其實(shí)施比傳統(tǒng)有限元方法更為簡(jiǎn)便；

2)采用統(tǒng)一DG框架求解流場(chǎng)控制方程是模擬非牛頓流動(dòng)問(wèn)題的一種有效方法．該方法使得三套方程之間的數(shù)據(jù)信息傳遞簡(jiǎn)單而且快捷，其程序?qū)崿F(xiàn)易于統(tǒng)一化；

3)平面Poiseuille流模擬的計(jì)算結(jié)果與精確解吻合很好，從而說(shuō)明該方法具有很好的穩(wěn)定性和較高的計(jì)算精度．與文獻(xiàn)吻合的平板收縮流模擬結(jié)果進(jìn)一步說(shuō)明本文的統(tǒng)一DG方法可以在非結(jié)構(gòu)網(wǎng)格上模擬包含應(yīng)力奇異點(diǎn)的復(fù)雜粘彈流動(dòng)問(wèn)題；

4)本文的統(tǒng)一DG方法的實(shí)施并不局限于Oldroyd-B本構(gòu)模型，也可推廣到其它微分型本構(gòu)方程的求解．并且，該方法不依賴(lài)于問(wèn)題的維數(shù)，故也可方便地?cái)U(kuò)展到高維問(wèn)題的求解．

參考文獻(xiàn)：

[1]趙曉凱，歐陽(yáng)潔，李五明.平板收縮流的對(duì)數(shù)構(gòu)象模擬[J].工程數(shù)學(xué)學(xué)報(bào)，2012,29(5):703-714 Zhao X K,Ouyang J,Li W M.Numerical simulation of planar contraction flow using log-conformation tensor approach[J].Chinese Journal of Engineering Mathematics,2012,29(5):703-714

[2]Cheng J,Shu C W.High order schemes for CFD:a review[J].Chinese Journal of Computational Physics,2009,26(3):633-655

[3]Yurun F.A comparative study of the discontinuous Galerkin and continuous SUPG finite element methods for computation of viscoelastic flows[J].Computer Methods in Applied Mechanics and Engineering,1997,141(1):47-65

[4]Kim J M,Kim C,Kim J H,et al.High-resolution finite element simulation of 4:1 planar contraction flow of viscoelastic fluid[J].Journal of Non-Newtonian Fluid Mechanics,2005,129(1):23-37

[5]Bird R B,Armstrong R C,Hassager O.Dynamics of polymeric liquids[J].Fluid Mechanics,1987,95(1):310-349

[6]Na Y,Yoo J Y.A finite volume technique to simulate the flow of a viscoelastic fluid[J].Computational Mechanics,1991,8(1):43-55

[7]Cockburn B,Shu C W.The Runge-Kutta local projection p1-discontinuous Galerkin finite element method for scalar conservation laws[J].Journal of Numerical Analysis,1991,25(3):337-361

[8]Cockburn B,Shu C W.TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II:general framework[J].Mathematics of Computation,1989,52(186):411-435

[9]Cockburn B,Shu C W.The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems[J].Journal of Computational Physics,1998,141(2):199-224

[10]Hesthaven J S,Warburton T.Nodal Discontinuous Galerkin Methods:Algorithms,Analysis,and Applications[M].New York:Springer,2008

[11]Shu C W.Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[J].Journal of Scientific Computing,2001,16(3):173-261

[12]Karniadakis G E,Israeli M,Orszag S A.High-order splitting methods for the incompressible Navier-Stokes equations[J].Journal of Computational Physics,1991,97(2):414-443

[13]Arnold D N,Brezzi F,Cockburn B,et al.Unified analysis of discontinuous Galerkin methods for elliptic problems[J].SIAM Journal on Numerical Analysis,2002,39(5):1749-1779

[14]Baaijens FPT.Mixed finite element methods for viscoelastic flow analysis:a review[J].Journal of Non-Newtonian Fluid Mechanics,1998,79(2-3):361-385

[15]Phillips T N,Williams A J.Comparison of creeping and inertial flow of an Oldroyd-B fluid through planar and axisymmetric contractions[J].Journal of Non-Newtonian Fluid Mechanics,2002,108(1):25-47

[16]Aboubacar M,Webster M F.A cell-vertex finite volume/element method on triangles for abrupt contraction viscoelastic flows[J].Journal of Non-Newtonian Fluid Mechanics,2001,98(s2-3):83-106