李路程,許曉陽

(1-西北工業(yè)大學應用數(shù)學系，西安 710129;2-陜西理工學院數(shù)學與計算機科學學院，漢中 723001)

1 引言

采用Brown構形場(Brownian configuration fields,BCF)方法[1,2]模擬粘彈性流體問題時，需要耦合求解宏觀尺度的守恒方程和介觀分子尺度的構形方程.該方法繞開了采用本構方程封閉求解流場控制方程的難點，因而可以摒棄本構方程建立中潛在的不正確性.同時，該方法從介觀分子尺度求解聚合物分子構形，故可以得到聚合物的分子信息.

守恒方程和構形方程耦合求解的常見方法有有限元法(finite element method,FEM)[1,2]和有限體積法(finite volume method,FVM)[3-5].FEM能很好地適應不規(guī)則區(qū)域，而在對流項的處理及對守恒方程的原始變量法求解方面實施起來比較繁瑣.雖然FVM的格式精度不高，但是它具有占用內存少、穩(wěn)定性強和對流項處理簡單的優(yōu)勢；同時該方法求解的原守恒方程在每個控制體上都守恒，且物理意義明確.因此，將FVM與BCF方法相結合，是獲得聚合物分子信息的有效途徑.目前，基于BCF的有限體積方法(finite volume method based on Brownian configuration fields,FVMBCF)鮮有文獻研究.代向艷等[6,7]運用FVMBCF方法模擬了Couette流和Poiseuille流，張僡鳳等[8]研究了FVMBCF方法的并行算法，并用之模擬了管道流.文獻[6–8]中流動問題的求解區(qū)域簡單，流動特性單一.因此，本文擬應用FVMBCF方法模擬具有復雜流動特性的平板收縮流動問題，以考證該方法模擬復雜流動問題的可靠性.

2 守恒方程

對于不可壓縮的等溫粘彈性流動，在忽略體積力的情況下，無量綱形式的流場守恒方程[1]為連續(xù)方程

動量方程

其中u表示流體速度，Reynolds數(shù)Re表示流體慣性力與粘性力的比值，p表示壓力，τp表示聚合物應力，粘度比β表示溶劑粘度與流體總粘度的比值.

粘彈性流動問題的流場控制方程包含守恒方程和封閉的應力本構方程.粘彈性流體本構關系的復雜性和機理的不確定性，增加了本構方程建立的難度，BCF方法恰恰繞開了這個難題.

3 BCF方法

BCF方法通過計算大量聚合物分子的構形，根據(jù)其系綜平均得到聚合物應力.再將介觀尺度獲得的應力與宏觀尺度的守恒方程結合，從而求解粘彈性流場的控制方程.BCF方法求解得到的構形在空間上光滑，確保計算所得的應力光滑[3]；該方法還可以避免CONNFFESSIT(calculation of non-Newtonian fluid:finite element and stochastic simulation techniques)方法[9]中需要追蹤大量粒子軌道的問題.

3.1 分子模型

常用的聚合物分子模型有圖1所示的珠-簧(鏈)模型和珠-桿(鏈)模型.為了表征聚合物分子的粘彈特性，本文選用有限拉伸非線性彈性(finite extensible nonlinear elastic,FENE)珠-簧模型.

圖1:聚合物分子模型示意圖

FENE珠-簧模型的彈簧力[1]

其中Q表示彈簧的構形向量，|Q|是彈簧的長度，Q0是彈簧的最大拉伸量，H是彈簧的彈性系數(shù).

3.2 分子構形演化方程

推導分子構形演化方程，首先需要分析稀溶液中聚合物分子模型中珠子的受力情況.根據(jù)Newton第二定律，加速度為零時，珠子i所受的合力公式

其中和分別代表珠子i所受的彈力、粘滯力和Brown隨機力.

彈力

由式(3)給出；粘滯力Brown隨機力

計算中常近似為

其中ri表示珠子i的位置向量，是i號珠子的速度，表示珠子i處溶劑的速度，ζ為小球與溶劑的摩擦系數(shù)，kB為Boltzmann常數(shù)，T表示絕對溫度，δ表示單位張量，是Markov過程，服從分布

整理方程(4)，得到

上面兩式相減，令Q=r2?r1，根據(jù)隨體導數(shù)公式得到

(7)式兩端同時除以又有得到

對(8)式進行無量綱處理，得到

其中Weissenberg數(shù)，L為流場的特征長度，U為流場的特征速度，Q表示分子模型彈簧的構形向量，κ表示速度梯度的轉置，F(xiàn)表示分子模型中彈簧的彈力.

求解方程(9)得到分子構形向量Q后，則由Kramers表達式[3]計算聚合物應力

其中nkBT為求解參數(shù)，I表示二階單位矩陣，〈Λ〉表示構形空間上的系綜平均，運算符?表示兩個向量的外積，其表達式為

3.3 施加方差縮減技術的BCF方法

1)設置構形場：在任一計算節(jié)點處，設置Nf個編號為j(j=1,K,Nf)的珠-簧模型；在所有節(jié)點處編號為j(j=1,K,Nf)的珠-簧模型都具有相同的初始構形，并經歷相同的隨機過程；

2)引入構形控制變量：初始時刻流體處于靜止狀態(tài)，分子構形服從Gauss隨機分布N(0,1)；流動過程中，分子構形的演化由方程(9)來描述；流動趨于穩(wěn)定時，分子構形最終服從均勻隨機分布U(0,1).這表明當前狀態(tài)應力的變化與平衡狀態(tài)應力的變化相關，因此，可以選取初始狀態(tài)服從均勻隨機分布，且流動過程中與Q經歷相同隨機過程的構形作為Q的控制變量[9].針對FENE模型，滿足方程

根據(jù)式(10)及其恒等變形，有

設根據(jù)式(12)和外積定義，應力各分量為

至此，給出了施加方差縮減技術的BCF方法計算聚合物應力各分量的表達式.根據(jù)統(tǒng)計學知識，控制變量的引入，不會改變統(tǒng)計量

的數(shù)學期望；同時，由于控制變量法的統(tǒng)計學性質，它們三個統(tǒng)計量的隨機誤差將大大地降低.因此，控制變量法的統(tǒng)計學性質使得構形控制變量的引入，能夠有效地減少應力的隨機誤差.

4 數(shù)值計算方法

基于Brown構形場方法求解粘彈性流動問題時，首先通過求解聚合物分子的構形方程(9)和控制變量方程(11)，再由式(13)計算聚合物應力，然后將介觀尺度求解的應力與宏觀尺度的守恒方程(1),(2)結合，從而得到描述粘彈性流動問題的宏觀物理量.

4.1 同位網格FVM

本文采用同位網格FVM[5]離散守恒方程(1),(2)，其中對流項的離散采用一階迎風格式，該格式簡單且易于實施.計算中，所有變量置于圖2所示的同一套網格.圖2中實線交點為計算物理量的存儲點，虛線為控制體界面，節(jié)點位于控制體的中心.

我撕到最后一頁時，聽到有人在樓下跟我媽講話，聽我媽那反應就知道來人不是扒鍋街的人，她平時的嗓門震得門框都會抖三抖，哪像那會，明明想笑卻拼命壓著喉嚨，像只有痰的老母雞。

圖2:同位網格示意圖

將連續(xù)方程(1)在圖2所示的P控制體上進行積分，其離散格式

二維動量守恒方程(2)的通量表達形式

其中

時間項的離散格式精度為一階；用控制容積積分法對對流-擴散項進行離散

精度為一階；其中

源項中應力采用顯式差分格式離散，壓力采用線性插值計算，精度均為一階；從而有限體積法部分的整體格式精度為一階.BCF方法是隨機方法，難以進行精度分析.設維數(shù)d，流動時間步FLAG，收斂判定條件數(shù)flag，構形數(shù)目Nf，計算點個數(shù)M，本文算法的計算量約為

4.2 FVMBCF方法的算法流程

基于并行程序的FVMBCF方法模擬粘彈性流動的具體算法流程如下：

步驟1準備階段設用P個處理機進行計算，在根進程上給定初始的速度u0，在每個進程上給定初始的分子構形和構形控制變量，初始構形和構形控制變量在空間上分別服從Gauss隨機分布和均勻隨機分布，它們的數(shù)值可由相應的分布函數(shù)N(0,1)和U(0,1)獨立取樣得到；

步驟2計算速度和壓力將初始或上一時間步的聚合物應力代入方程(1)和(2)，并用FVM求解，即可得到新的速度和壓力；

步驟3計算構形和構形控制變量用新的速度u計算速度梯度κ，并將速度梯度κ廣播到每一個進程中去.當已知速度梯度κ及上一時間步的分子構形和構形控制變量時，根據(jù)方程(9)和(11)，可計算得到當前時間步的分子構形和構形控制變量

步驟4計算聚合物應力在每個進程上，將步驟3中計算得到的當前時間步的分子構形和構形控制變量代入式(13)，并歸約到0進程上求和，計算聚合物應力

步驟5若收斂條件滿足，計算結束；否則返回步驟2.

5 數(shù)值算例

5.1 Couette流

為了驗證施加方差縮減技術的FVMBCF方法的有效性，本文首先基于FENE珠-簧模型，模擬聚合物稀溶液的平板Couette流.

在該流動問題中，聚合物流體被限制在兩個距離為L=1.0的平行平板間.當t＜0時，流體和兩個平板靜止；當t=0，下平板開始以恒定的水平速度u=1.0沿x軸正方向移動，壁面邊界滿足無滑移邊界條件，初始狀態(tài)為平衡態(tài).為了與文獻中的結果進行對比，本文采用文獻[11]中的模型參數(shù)，即

圖3(a)刻畫了不同時刻、不同y值處的速度.圖3(b)描述了四個不同位置y=0.2,y=0.4,y=0.6,y=0.8處的速度演化.由圖3(b)可知，越靠近下平板的位置，速度過沖現(xiàn)象發(fā)生的越早.圖3(c)描述了圖3(b)中的四個位置處的應力演化，應力曲線光滑.模擬結果與文獻[11]的結果吻合，從而驗證了施加方差縮減技術的FVMBCF方法的有效性.

圖3:平板Couette流速度和應力分布

5.2 平板收縮流

4:1平板收縮流能夠反映流體經旋轉、剪切和拉伸變形后的流動特性.收縮口附近存在的應力奇點，使流動問題的數(shù)值解在收縮口附近存在大梯度變化.因此，對4:1粘彈性收縮流動問題的模擬，更能檢驗本文數(shù)值方法的有效性.

同時，由于應力奇點的存在及FVMBCF方法在計算過程中不可避免地產生隨機誤差，導致在收縮口附近應力出現(xiàn)振蕩.所以本節(jié)先給出施加方差縮減技術的FVMBCF方法的數(shù)值模擬結果，然后在5.3節(jié)給出關于施加方差縮減技術必要性的討論.

由于收縮流的對稱特性，僅選取圖4所示y≥0的計算區(qū)域.將計算區(qū)域劃分成四個部分，采用疏密不同的拼接網格.四個區(qū)域的網格分別取為：42×12、42×77、32×12和22×12.計算中固壁面邊界滿足無滑移邊界條件.無量綱參數(shù)選取

收斂判定條件閾值5×10?5.入口速度為

圖4:收縮流動示意圖

圖5給出了不同We下收縮口附近的流線分布.由圖5可知，在收縮口上游和下游處流動充分發(fā)展；在收縮口附近形成強拉伸區(qū)域；在上游拐角處，流動以旋轉為主，形成角渦，且隨著We的增大，渦流區(qū)的范圍減小，渦心位置上升，角渦變小.

圖5:不同We下收縮口附近的的流線分布

圖6給出了We=1時收縮口附近的應力等值線分布.由圖6可以看出，在收縮口附近，應力等值線光滑且分布密集，說明在收縮口處應力變化劇烈.圖5與圖6的數(shù)值結果與文獻[12,13]相一致.

圖6:We=1.0時收縮口附近的應力等值線分布

圖7描述了流動穩(wěn)定時收縮口附近y=1水平層面的速度、剪切應力及第一法向應力差的分布情況.由圖7可知，隨著We的增大，速度過沖現(xiàn)象更加明顯.另外，剪切應力和第一法向應力差過沖現(xiàn)象加劇，收縮口附近的應力梯度很大，各應力分量均在收縮口附近取得最大值，且峰值隨著We的增大而增加.這說明在收縮口附近，剪切拉伸作用隨We的增大而加強.

圖7:收縮口附近y=1的速度u、剪切應力τxy、第一法向應力差分布

5.3 方差縮減技術

圖8描述了未施加及施加方差縮減技術情形下收縮口附近應力等值線的分布情況(左列是未施加方差縮減技術的結果，右列是施加方差縮減技術的結果).由圖8可知，未施加方差縮減技術時，收縮口附近的應力等值線出現(xiàn)振蕩現(xiàn)象；而施加方差縮減技術后，振蕩現(xiàn)象消失，應力各分量等值線均變得光滑.這說明方差縮減技術能夠很好地抑制收縮口附近的應力振蕩.

圖9給出了y=0.52位置處應力的數(shù)值結果，圖中小方框內曲線是對收縮口附近(圖中圓形區(qū)域)應力分布情形的放大處理.由圖9可知施加方差縮減技術后應力的數(shù)值曲線較未施加的應力數(shù)值曲線光滑，振蕩現(xiàn)象消失；但應力各分量的數(shù)值減小.這說明了FVMBCF方法施加方差縮減技術時具有能夠確保獲得光滑應力解的優(yōu)點，但是同時存在著消弱應力解數(shù)值精度的缺陷.

圖8:未施加及施加方差縮減技術情形下收縮口附近應力等值線的分布情況

圖9:y=0.52位置處應力的數(shù)值結果

施加方差縮減技術時，由于控制變量參與計算，使得構形向量的總計算量將近增加一倍.但是模擬結果顯示：同樣基于四個并行進程，在選取相同收斂閾值的條件下，施加方差縮減技術的數(shù)值模擬耗時(57285.5s)比未施加方差縮減技術的數(shù)值耗時(56237.4s)僅多1.8%.這表明方差縮減技術能夠在計算時間微量增加的條件下，很好地抑制應力的振蕩現(xiàn)象，其穩(wěn)定性能良好.

6 結束語

本文給出了施加方差縮減技術的FVMBCF方法及其求解流程，模擬了Couette流和4:1收縮流的粘彈性流動問題.數(shù)值結果表明：

1)該方法可以準確地模擬包含應力奇異點的復雜流動問題，計算結果可靠，算法易于實現(xiàn)；

2)該方法具有穩(wěn)定性強和對流項處理簡單的優(yōu)點，且應力數(shù)值解具有光滑特性；

3)方差縮減技術能夠很好地抑制收縮口附近的應力振蕩現(xiàn)象.

參考文獻：

[1]Abedijaberi A,Khomami B.Sedimentation of a sphere in a viscoelastic fluid:a multiscale simulation approach[J].Journal of Fluid Mechanics,2012,694(1):78-99

[2]Abedijaberi A,Khomami B.Continuum and multi-scale simulation of mixed kinematics polymeric flows with stagnation points:closure approximation and the high Weissenberg number problem[J].Journal of Non-Newtonian Fluid Mechanics,2011,166(11):533-545

[3]Owens R G,Phillips T N.Computational Rheology[M].London:Imperial College Press,2002

[4]陶文銓.數(shù)值傳熱學(第2版)[M].西安：西安交通大學出版社，2001 Tao W Q.Numerical Heat Transfer(2nd Edition)[M].Xi’an:Xi’an Jiaotong University Press,2001

[5]宋道云.同位網格有限體積法及其在粘彈性液體收縮流數(shù)值模擬中的應用研究[D].上海：華東理工大學，2002 Song D Y.The algorithm of collocated-grid finite volume method and its application of numerical simulation in a contraction flow for viscoelastic fluids[D].Shanghai:East China University of Science and Technology,2002

[6]代向艷，歐陽潔，張小華，等.基于Brown構形場的有限體積法在聚合物流動模擬中的應用[J].工程數(shù)學學報，2011,28(5):642-654 Dai X Y,Ouyang J,Zhang X H,et al.Simulation of the polymeric flows using finite volume method based on Brownian con figuration fields[J].Chinese Journal of Engineering Mathematics,2011,28(5):642-654

[7]代向艷，歐陽潔.平板Couette流場中聚合物流變性質及分子構象的數(shù)值模擬[J].高分子材料科學與工程，2011,27(3):157-161 Dai X Y,Ouyang J.Rheological properties and molecular con figuration of polymer solution in planar Couette flow[J].Polymer Materials Science and Engineering,2011,27(3):157-161

[8]張僡鳳，歐陽潔，代向艷.耦合有限體積法的Brown構型場并行算法研究[J].計算物理學報，2012,29(1):17-24 Zhang H F,Ouyang J,Dai X Y.Paralle algorithm for Brownian con figuration fields with finite volume method[J].Chinese Journal of Computational Physics,2012,29(1):17-24

[9]Bonvin J,Picasso M.Variance reduction methods for CONNFFESSIT-like simulations[J].Journal of Non-Newtonian Fluid Mechanics,1999,84(2-3):191-215

[10]ttinger H C,van den Brule B,Hulsen M A.Brownian con figuration fields and variance reduced CONNFFESSIT[J].Journal of Non-Newtonian Fluid Mechanics,1997,70(3):255-261

[11]Tran-Canh D,Tran-Cong T.Element-free simulation of dilute polymeric flows using Brownian con figuration fields[J].Korea-Australia Rheology Journal,2004,16(1):1-15

[12]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(2-3):83-106

[13]Edussuriya S S,Williams A J,Bailey C.A cell-centred finite volume method for modelling viscoelastic flow[J].Journal of Non-Newtonian Fluid Mechanics,2004,117(1):47-61