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

后翠紅,蘇 進

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

0 引 言

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

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

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

1 Brown分子模型

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

F(Q)=HQ

(1)

(2)

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

構型向量場的演化方程可表示為

(3)

式中：κ11、κ12、κ22為速度應變參數。

(4)

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

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

(5)

(6)

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

確定構型場,則聚合物偏應力張量為

(9)

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

2 分子模型的MLMC方法

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

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

(10)

2.2 MLMC方法

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

(11)

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

(12)

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

(13)

由式(13)知

(14)

對于誤差ε,有

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

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

(15)

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

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

(16)

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

(17)

式中：ρm是誤差密度函數。

平均步數M*與誤差密度和離散誤差的關系為

(18)

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

(19)

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

(20)

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

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

(21)

對于最高層L,存在常數a,有

(22)

通過拉格朗日方法求樣本最優解

(23)

由式(23)得到最優樣本數為

(24)

的唯一整數。

故由式(12)知,

(25)

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

(26)

式中,α,β,γ,C1,C2,C3均為正常數。

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

(27)

(28)

因此相應的MLMC計算公式分別為

(29)

(30)

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

(31)

(32)

故對應的MLMC計算公式分別為

(33)

(34)

2.3 MLMC算法基本步驟

步驟1 各參數初始化。

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

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

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

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

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

步驟7 根據每層的方差Vl和成本Cl計算并更新每層最優樣本數Ns。

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

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

3 數值結果

3.1 正確性檢驗

為了驗證程序的正確性和計算結果的可靠性,用Hooke模型計算穩態剪切流場中的剪切應力,并與文獻[7]中解析解進行對比。Hooke模型穩態Poiseuille流的解析解為

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

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

(35)

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

3.2 有效性分析

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

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

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

(36)

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

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

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

表 2 FENE模型的參數Tab.2 Parameters of FENE model

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

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

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

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

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

4 結 語

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