聚合物流體數(shù)值模擬的多層蒙特卡羅方法

后翠紅,蘇 進

(西安工程大學(xué) 理學(xué)院,陜西 西安 710048)

0 引 言

為了預(yù)測加工過程[1]和微流動控制[2]等領(lǐng)域內(nèi)聚合物的復(fù)雜流動規(guī)律,數(shù)值仿真技術(shù)逐漸成為一種重要研究手段。1993年,Laso和?ttinger[3-4]基于高分子的微觀模型,提出非牛頓流體的計算方法，即有限元和隨機模擬方法求解聚合物流體的流動問題。該方法需要追蹤大量分子軌跡,得到的應(yīng)力可能不光滑且成本昂貴,所以無法應(yīng)用到工業(yè)中。1997年,Hulsen等[5]提出了描述聚合物應(yīng)力的Brown構(gòu)型場(Brownian configuration field, BCF)模型。其基本思想是在宏觀上用有限元方法計算速度場,在微觀上用布朗分子模擬計算大量聚合物分子的構(gòu)型,再通過蒙特卡羅系綜平均方法得到聚合物分子對應(yīng)力的貢獻。該方法克服了追蹤分子軌跡問題,得到的應(yīng)力通常是空間光滑。BCF方法逐步成為研究聚合物流體流變性質(zhì)的有力工具。索哲[6]利用BCF方法研究了簡單流場中的應(yīng)力參數(shù);代向艷等[7]通過珠-簧鏈模型建立了基于BCF的有限體積方法;張僡鳳等[8]建立了耦合有限體積法的BCF方法的并行算法;Somali等[9]基于有限元方法的BCF研究了Hooke模型和FENE模型在不同條件下的線性與非線性穩(wěn)定性。BCF方法中的一個重要環(huán)節(jié)是利用隨機模擬方法，計算大量聚合物分子的構(gòu)型,并通過其蒙特卡羅方法得到聚合物分子對應(yīng)力的貢獻[5,10]。然而,在低Wi數(shù)問題中，由于隨機項占優(yōu)會帶來隨機誤差大和計算成本高的不足。Bonvin等[11]早期提出了一類基于控制變量思想的蒙特卡羅方法方差縮減技術(shù)使得隨機誤差得到部分改善，但并沒有綜合考慮減小隨機誤差并優(yōu)化計算成本。

2008年,針對傳統(tǒng)蒙特卡羅方法的方差縮減技術(shù),Giles[12]提出了多層蒙特卡羅(multilevel Monte Carlo,MLMC)的基本思想,即保持精度不變采用較大時間步長的計算結(jié)果以降低計算量,最終縮減成本。在金融領(lǐng)域中,2009年,Giles[13]證明了MLMC方法很容易拓展到幾個基本量的加權(quán)平均籃子期權(quán);2010年,Korn[14]將MLMC方法應(yīng)用到了金融和保險的理論研究中;2014年,孫健蘭[15]驗證了MLMC方法計算期權(quán)Greeks問題的有效性;2016年,宋斌等[16]將MLMC方法應(yīng)用于巴黎期權(quán)的定價計算中,擴展了巴黎期權(quán)的數(shù)值算法適用范圍,并提高了巴黎期權(quán)定價的計算精度;2017年,Kazeem[17]采用MLMC方法研究了路徑依賴期權(quán)的定價和期權(quán)價格敏感度問題。MLMC方法在金融領(lǐng)域的應(yīng)用越來越多,但缺乏在聚合物流動應(yīng)力計算中的深入應(yīng)用研究。聚合物應(yīng)力的BCF方法屬于微觀模型,與傳統(tǒng)宏觀模型相比,在微流體控制等鄰域的應(yīng)用更為廣泛。然而傳統(tǒng)BCF方法存在隨機誤差大及計算成本高的缺點。

本文為克服傳統(tǒng)方法的不足,在聚合物流動應(yīng)力的計算中,引入了一種基于Euler格式隨機Brown軌道的多層蒙特卡羅(multilevel Monte Carlo, MLMC)方法。首先通過Hooke模型驗證程序的正確性和計算結(jié)果的可靠性,然后研究了MLMC方法在計算聚合物應(yīng)力方面的計算效率,最后討論了Wi數(shù)對該方法計算成本的影響。

1 Brown分子模型

采用Hooke啞鈴模型和FENE啞鈴模型描述聚合物大分子鏈,該模型是將聚合物分子鏈簡化為一根無質(zhì)量的彈簧連接的2個質(zhì)量為m的剛性小珠,用Q表示2個小珠間的構(gòu)型向量,則Hooke模型和FENE模型的彈簧力分別為

F(Q)=HQ

(1)

(2)

式中:H表示彈簧的彈性常數(shù);b表示彈簧的最大拉伸量的平方。

構(gòu)型向量場的演化方程可表示為

(3)

式中：κ11、κ12、κ22為速度應(yīng)變參數(shù)。

(4)

式中：Wi=λU/L為Weissenberg數(shù)，表示聚合物大分子的松弛時間和特征時間的比值。式(4)右端第一項表示速度梯度引起的變化，第二項表示彈性恢復(fù)構(gòu)成的漂移分量，第三項表示Brown熱運動引起的擴散分量。

設(shè)Q(t)=(p(t),q(t)),W(t)=(w(t),v(t)),其中p(t)和q(t)分別表示在水平和垂直方向上的構(gòu)型分量,w(t)和v(t)表示2個相互獨立的Wiener過程。Brown構(gòu)型場Hooke模型演化方程分量形式為

(5)

(6)

FENE模型演化方程分量形式為

確定構(gòu)型場,則聚合物偏應(yīng)力張量為

(9)

式中：?表示向量的并矢運算；E[Q?F(Q)]表示Q?F(Q)的期望。

2 分子模型的MLMC方法

2.1 隨機微分方程離散的Euler格式

為了得到方程(4)的數(shù)值格式,需要將時間區(qū)間[0,T]離散。選取自然數(shù)M以確定時間步長Δt=T/M,得到離散時間點0=t0

(10)

2.2 MLMC方法

為了降低計算量并提高計算效率,考慮具有不同時間步長Δt=M-lT,l=1,2,…,L,M∈[2,8]的MLMC路徑模擬。對于隨機變量Q的函數(shù)g(Q),給出第l層的逼近序列{gl(Q)},其最細層的期望可表示為最粗層的期望加預(yù)期修正的和,即

(11)

用無偏估計估計E[gL(Q)],即

(12)

式中：Nl表示第l層的樣本數(shù)。以V(·)表示隨機變量的方差，則標準理論給出的均方誤差(MSE)為

(13)

由式(13)知

(14)

對于誤差ε,有

ε=E[g(Q)]-λML=E[g(Q)-gL(Q)]+

(gL(Q)-λML)=εT+εS

(15)

式中:εT表示因離散而造成的離散誤差；εS表示求N個樣本期望的統(tǒng)計誤差。

要使得MLMC方法均方誤差達到最小值ε,則

(16)

式中：V[gL(Q)]表示隨機變量函數(shù)gL(Q)的方差。采用不同時間步長對軌道離散，從而控制離散誤差εT,得到

(17)

式中：ρm是誤差密度函數(shù)。

平均步數(shù)M*與誤差密度和離散誤差的關(guān)系為

(18)

式中：m=1,2,…,M。如果

(19)

則對該步長進行細分;如果

(20)

則所有步長的誤差都比平均步長的誤差小,即保證整個離散誤差在εT內(nèi)。

通過控制樣本數(shù)控制統(tǒng)計誤差εS,當滿足V[gL(Q)]≤εS/2時即滿足設(shè)定的精度。對于第l層的步數(shù)Ml,根據(jù)以上時間步長的構(gòu)造,有

(21)

對于最高層L,存在常數(shù)a,有

(22)

通過拉格朗日方法求樣本最優(yōu)解

(23)

由式(23)得到最優(yōu)樣本數(shù)為

(24)

的唯一整數(shù)。

故由式(12)知,

(25)

在|E[gL(Q)-g(Q)]|≤C12-αl,Vl≤C22-βl,Cl≤C32γl的特定情況下,隨著l→∞,則總計算成本C為

(26)

式中,α,β,γ,C1,C2,C3均為正常數(shù)。

在Hooke模型中,剪切應(yīng)力τxy和第一法向應(yīng)力差τxx-τyy對應(yīng)的g(Q)分別為

(27)

(28)

因此相應(yīng)的MLMC計算公式分別為

(29)

(30)

在FENE模型中，剪切應(yīng)力和第一法向應(yīng)力差對應(yīng)的g(Q)分別為

(31)

(32)

故對應(yīng)的MLMC計算公式分別為

(33)

(34)

2.3 MLMC算法基本步驟

步驟1 各參數(shù)初始化。

步驟2 若層數(shù)l<1,輸出Ns=N0;若層數(shù)l≥1,對附加樣本數(shù)dNl(樣本數(shù)與收斂測試樣本數(shù)的差)進行評估。

步驟3 若附加樣本數(shù)dNl的和小于0,輸出NL；若附加樣本數(shù)dNl的和大于0,則計算每層的附加樣本數(shù)。

步驟4 根據(jù)每層的附加樣本數(shù)計算并更新每層的樣本Nl,l=0,1,…,L。

步驟5 根據(jù)每層的樣本Nl計算每層的方差Vl,l=0,1,…,L。

步驟6 根據(jù)每層的樣本Nl計算每層的成本Cl,l=0,1,…,L。

步驟7 根據(jù)每層的方差Vl和成本Cl計算并更新每層最優(yōu)樣本數(shù)Ns。

步驟8 對剩余誤差進行弱收斂性檢驗。若收斂,結(jié)束程序;若不收斂,設(shè)置L:=L+1,返回步驟2重新計算并更新Ns。

步驟9 若得到的Ns大于收斂樣本N,則該層繼續(xù)模擬Ns-N條軌道,直到滿足設(shè)定精度。

3 數(shù)值結(jié)果

3.1 正確性檢驗

為了驗證程序的正確性和計算結(jié)果的可靠性,用Hooke模型計算穩(wěn)態(tài)剪切流場中的剪切應(yīng)力,并與文獻[7]中解析解進行對比。Hooke模型穩(wěn)態(tài)Poiseuille流的解析解為

數(shù)值模擬中選擇的參數(shù)為:初始樣本數(shù)N0=20,精度ε=0.005,收斂測試樣本數(shù)N=1 000,Wi=1; 速度梯度轉(zhuǎn)置的4個分量為κ11=0,κ12=-0.5,κ21=0,κ22=0; 時間步長Δt=2-l(l=1,2，…,L)。圖1給出了Hooke模型在剪切流場中剪切應(yīng)力隨時間變化圖。

由圖1可知，當t=29時,應(yīng)力的演化狀態(tài)已經(jīng)達到穩(wěn)態(tài)。由圖1還可以發(fā)現(xiàn)，MLMC數(shù)值方法得到的數(shù)值解與解析解在穩(wěn)態(tài)時的結(jié)果一致吻合,說明了本文方法計算聚合物偏應(yīng)力的可靠性。圖2進一步給出了Hooke模型剪切流場中剪切應(yīng)力的均方剩余誤差隨時間演化圖，其中剪切應(yīng)力的均方剩余誤差Err(τxy)為

(35)

圖2顯示,隨時間增加,均方剩余誤差逐漸減小并趨于穩(wěn)定。

3.2 有效性分析

為了討論MLMC方法在簡單流場中聚合物流體應(yīng)力計算的效率及適應(yīng)性,分別計算了Hooke和FENE模型在剪切流場和拉伸流場的剪切應(yīng)力及第一法向應(yīng)力差的計算成本,并與標準蒙特卡羅(standard Monte Carlo, StdMC)方法相比較。

3.2.1 Hooke啞鈴模型 在Hooke模型中,設(shè)定精度ε=0.05,時間t=15,時間步長Δt=2-l,l=1,2，…,L,樣本數(shù)N=9 000,初始樣本數(shù)N0=20。計算剪切應(yīng)力和第一法向應(yīng)力差時,選擇Wi數(shù)及速度梯度4個分量κ11、κ12、κ21、κ22的參數(shù)如表1所示。

圖3給出了剪切流場中Hooke模型的剪切應(yīng)力計算成本(C)在不同精度(ξ)下的對比結(jié)果。由圖3可知,在不同精度下,MLMC方法比StdMC方法節(jié)約計算成本。為了進一步說明MLMC方法的計算效率,給出了具有不同精度的剪切應(yīng)力和第一法向應(yīng)力差的總成本降低率φC(%),即

(36)

經(jīng)計算，精度為0.05時,MLMC方法的計算成本比StdMC方法降低了大約68.6%。

圖4給出了拉伸流場中Hooke模型的第一法向應(yīng)力差計算成本在不同精度下的對比結(jié)果。圖4表明,在不同精度下,MLMC方法比StdMC方法節(jié)約計算成本,在同一精度下計算成本降低了大約99.4%。通過以上分析,發(fā)現(xiàn)在簡單流場中,當Wi=1時(低Wi數(shù))MLMC方法均比StdMC方法更節(jié)約計算成本。

3.2.2 FENE啞鈴模型 在FENE模型中,設(shè)定精度ε=0.005,時間t=2,時間步長Δt=2-l,l=0,1,…,L,收斂樣本N=9 000。計算剪切應(yīng)力和第一法向應(yīng)力差時,選擇Wi數(shù)、速度梯度4個分量κ11、κ12、κ21、κ22,初始樣本N0及最大拉伸長度b等 參數(shù)，如表2所示。

表 2 FENE模型的參數(shù)Tab.2 Parameters of FENE model

圖5給出了剪切流場中FENE模型的剪切應(yīng)力計算成本在不同精度下的對比結(jié)果。圖5表明,在不同精度下,MLMC方法比StdMC方法更節(jié)約計算成本。精度為0.01、0.04、0.05時,計算成本分別降低了約99.3%、92.6%、72.4%。圖6給出了拉伸流場中FENE模型的第一法向應(yīng)力差計算成本在不同精度下的對比結(jié)果。

圖6表明,在不同精度下,MLMC方法比StdMC方法節(jié)約了計算成本。精度高于0.04時,MLMC方法比StdMC方法節(jié)約計算成本約37.3%；精度低于0.04 時,MLMC方法比StdMC方法節(jié)約計算成本約15.2%。通過以上分析,發(fā)現(xiàn)在簡單流場中,當Wi=1/2時(低Wi數(shù))，MLMC方法比StdMC方法更節(jié)約計算成本。

圖7給出了剪切流場中，Hooke模型剪切應(yīng)力計算成本隨Wi數(shù)的變化關(guān)系。如圖7所示,Wi=2時,StdMC方法與MLMC方法的計算成本基本相近,而且MLMC方法的計算成本達到最小值；當Wi<2時,二者的計算成本均呈現(xiàn)隨Wi數(shù)的增大而減少的趨勢,但MLMC方法的計算成本比StdMC方法低；當Wi>2時,計算成本都隨Wi的增大而增加,MLMC方法的計算成本反而高于StdMC方法。圖8是剪切流場中，F(xiàn)ENE模型剪切應(yīng)力計算成本隨Wi數(shù)的變化圖。由圖8可知，數(shù)值結(jié)果與Hooke模型相類似：當Wi<1時,MLMC方法比StdMC方法更節(jié)約計算成本；然而,當Wi>1時,MLMC方法并不能減小計算成本。

3.2.3 參數(shù)分析Wi數(shù)是聚合物流體黏彈效應(yīng)的重要參數(shù),討論Wi對MLMC方法計算成本的影響。在Hooke模型中,除t=2,κ22=-0.25外,其他參數(shù)不變。使用FENE模型計算剪切流場的剪切應(yīng)力時，設(shè)定精度ε=4×10-4,除κ22=-0.7外其他參數(shù)不變。

綜上所述,隨著Wi數(shù)增加,Vl(Δτxy)在低層開始不發(fā)生變化,增加層數(shù)對于MLMC方法而言起不到方差縮減的效果,反而增加計算量。從物理學(xué)角度看,Wi數(shù)增加會引起聚合物彈性的增大,隨機項作用減弱,利用多層思想增加計算層數(shù)起不到縮減方差的效果。但是,當Wi數(shù)較低時,布朗力比彈性作用占優(yōu),隨機作用引起的誤差對層數(shù)比較敏感,采用MLMC方法可以降低方差,提高計算效率。因此,MLMC方法在解決低Wi數(shù)聚合物流動高精度應(yīng)力計算問題時具有更明顯的優(yōu)勢。

4 結(jié) 語

給出了一種基于Euler格式隨機Brown軌道計算聚合物流體應(yīng)力的多層蒙特卡羅(MLMC)方法。數(shù)值模擬結(jié)果表明,在低Wi數(shù)聚合物應(yīng)力計算中,與傳統(tǒng)蒙特卡羅(StdMC)方法相比，MLMC方法可以有效降低隨機誤差,節(jié)約計算成本。因此,該方法尤其適合解決低Wi數(shù)聚合物流動中的應(yīng)力高精度計算問題。在后續(xù)的研究中,可以進一步考慮將MLMC方法與流場的介觀計算方法相結(jié)合,研究微流動或微生物環(huán)境中大分子的流變學(xué)行為。