聚合物流体数值模拟的多层蒙特卡罗方法

后翠红,苏 进

(西安工程大学 理学院,陕西 西安 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方法与流场的介观计算方法相结合,研究微流动或微生物环境中大分子的流变学行为。