Stochastic Brownian dynamics is an extremely powerful way to simulate the polymer dynamics in solutions and melts. Mathematically, these models are described by stochastic differential equations. The most challenging problems are the Monte Carlo algorithm, which simulates the motion of a large number of model particles and hence requires an enormous amount of computer time. It is therefore necessary to develop an efficient numerical method in operational emergency response applications. In this paper, we give an improved multilevel Monte Carlo (improved-MLMC) method based on equilibrium control variables at each level to calculate the propagation of polymers. The improved-MLMC method can be shown to result in asymptotically optimal random errors and reduce total cost when compared to the standard Monte Carlo and MLMC methods. Finally, the effect of the Wi number (dimensionless parameter) on the total cost of the presented MLMC method is also discussed in detail.

In recent years, with the development of the processing industry, numerical simulation technology has gradually become an important way to study the polymer flow. Traditionally, a continuum approach to the mathematical modeling of polymeric liquids has been adopted in which the stress on a macroscopic fluid element is related to the deformation by a constitutive equation. Although those approaches are sufficient to provide a qualitative description of important flow phenomena in some situations, they are incapable of providing quantitative agreement with experimental measurements. In the 1990s, Öttinger1 proposed a very flexible and robust stochastic method to determine the polymeric contribution to the extra-stress rather than solving closed form constitutive equations. Mathematically, a stochastic differential equation for chain dynamics is solved via Brownian dynamics (BD) to obtain the polymeric stress at the mesoscopic level. The approach combines a macroscopic description of the kinematics with a coarse-grained microscopic description of the evolution of the stress. Laso and Öttinger2 termed this hybrid method CONNFFESSIT (Calculation of Non-Newtonian Flows: Finite Element and Stochastic Simulation Technique). Later, the BCF (Brownian Configuration Field) techniques3 have also been addressed to overcome the drawbacks of the CONNFFESSIT. Specifically, the BCF approach relies on an ensemble of continuous fields of spatially correlated macromolecules as opposed to individual polymer molecules to describe the polymer dynamics under flow. However, the stochastic simulations methods, as it was originally described in Refs. 2 and 3, are too expensive to allow for industrial applications because of the slow convergence rate of stochastic methods.

The computational bottleneck is the Monte Carlo algorithm, which simulates the motion of a large number of model dumbbells to obtain the polymer extra-stress. Efficient Monte Carlo algorithms for integrating the trajectories of the individual polymer are crucial for fast and precise forecasts. Hence, some variance reductions have been used to reduce the statistical error in a stochastic simulation without increasing the number of trajectories that need to be simulated.4,5 Hulsen et al.,3 Öttinger et al.,4 and Halin et al.6 also discussed the variance reduction technique based on control variables. The idea was to use strongly correlated local ensembles of dumbbells, and to subtract equilibrium ensembles of dumbbells. In order to obtain a more efficient method, Bonvin and Picasso7 extended these ideas to nonequilibrium ensembles of dumbbells. Chauvière8 also derived an efficient method for solving the stochastic differential, which describes the evolution of the configuration fields. A more efficient version of this method that is competitive with macroscopic computations was developed by Chauvière and Lozinski.9 

The multilevel Monte Carlo (MLMC) method is a recently established technique for efficient computation of an observable’s statistics by approximate sampling in the case when generation of samples of different accuracy is possible. The idea is based on the observation that coarse sample approximations can be used as control varieties for more accurate sample approximations. As a variance reduction technique of the Monte Carlo estimator, the MLMC method10–12 has been used successfully in a wide range of applications from mathematical finance13 to uncertainty quantification in subsurface flows.14 Anderson and Higham15 showed how to extend the multilevel Monte Carlo approach to the continuous time Markov chain. Pauli and Arbenz16 investigated a fault tolerant MLMC method; the proposed optimization procedure can react on experienced faults. Giles et al.17 constructed and analyzed the multilevel Monte Carlo method for the approximation of distribution functions and densities of univariate random variables. Ullmann and Papaioannou18 introduced a multilevel estimator that is motivated by the need to estimate the probability of rare events in engineering systems with random inputs. Most importantly, the MLMC method reduces the computational cost by at least one order of magnitude compared to the standard Monte Carlo method, greatly lowering the computational complexity needed to compute expected values of functions of the state of the system to a specified accuracy. In 2016, Giles et al.19 gave Brownian dynamics simulations of the finitely extensible nonlinear elastic (FENE) model with equilibrium fields using the multilevel Monte Carlo method incorporated within non-nested adaptive time steps. Belomestny and Nagapetyan20 discussed the possibility of using the multilevel Monte Carlo (MLMC) approach for weak approximation schemes. The aim of the present article is the further extension of the MLMC methodology for estimation of the polymer extra-stress combining with other variance reduction techniques, called an improved-MLMC method. In our work, to explore the Hooke and FENE model simulations in more detail, we will also give some results for molecules interacting with simple velocity fields.

The outline of this paper is as follows. In Sec. II, we introduce the stochastic model for simulating the chain dynamics in the frames of Hooke and FENE dumbbells for simple velocity fields. Section III introduces the proper treatment of variance reduction techniques using correlated local ensembles of dumbbells in each level context of the MLMC method and contains a review of the MLMC algorithm with particular focus on the integrating of shear stress and first normal stress difference. In Sec. IV, we present numerical results for two different flows corresponding to shearing action and biaxial extension. This includes comparisons between the StdMC and MLMC methods so as to identify the most efficient combination. The effect of the Wi number on the results and MLMC method are also discussed in this section. We finally conclude and outline ideas for further work in Sec. V.

The Hooke and FENE dumbbell models are used to describe polymer macro-molecular chains in this article. The spatial configuration of a dumbbell is given in terms of a connector vector G, linking pairs of beads in dumbbells. Starting from the linear (Hooke) force law, the spring force is given by

FG=HG.
(1)

The property of finite extensibility can be incorporated by introducing a finitely extensible nonlinear elastic (FENE) spring force,

FG=H1G2/bG,
(2)

where H is the spring constant and b is the square of the dimensionless segmental maximum extensibility.

The evolution equation21 for the configuration vector field is given by

dGx,t=κx,tGx,t+FGζdt+2KBZζdWt,
(3)

where ζ is the friction coefficient of the dumbbells’ beads in the solvent, KB is the Boltzmann constant, Z is the absolute temperature, F(G) is the spring force, x is the vertical abscissa, t is the current time, W(t) is a multi-dimensional Wiener process, and the components of dW(t) are independent processes, i.e., Gaussian. κ(t) is the velocity gradient transpose, where the simple shear flow and extensional flow fields are, respectively, given by

κS=0κ120 0,κE=κ11 00κ22.

If we scale the length of the connector vector G by KBZ/H, then the evolution equation (3) for the dimensionless configuration field becomes

dGx,t=κx,tGx,tFG2Widt+1WidWt,
(4)

where Wi = λU/L is the Weissenberg number and λ is the relaxation time of the bead–spring chain. The first term on the right-hand side of Eq. (4) is caused by the action of the flow field; the second term is contributed by elastic recovery; the third term is the thermal motion of Brown.

For convenience, the dumbbells are characterized by their horizontal and vertical elongations, p(y,t) and q(x,t), respectively; w(t), v(t) are two independent standard Wiener processes. Then, the component equations (4) of the Hooke model are:

dpx,t=κ11px,t+κ12qx,tpx,t2Widt+1Widwt,
(5)
dqx,t=κ21px,t+κ22qx,tqx,t2Widt+1Widvt.
(6)

Similarly, the component form of the FENE model is

dpx,t=κ11px,t+κ12qx,t12Wipx,t1p2x,t+q2x,t/bdt+1Widwt,
(7)
dqx,t=κ21px,t+κ22qx,t12Wiqx,t1p2x,t+q2x,t/bdt+1Widvt.
(8)

The Kramers expression for the polymeric contribution to the extra-stress tensor becomes

τp=1WiI+GFG,
(9)

where ⊗ is the tensor product of two vectors, n is the number of polymer molecules per unit volume, I is the second order identity tensor, and GFG denotes the GF(G) mathematical expectation.

We now present the original discretization of Eq. (4) in the frame of Ref. 2; one can see Refs. 21 and 22 for details about the discretization of stochastic differential equations. The configuration vector G is discretized into 0 = t0 < t1 < < tm < < tM = T in time [0, T]. Then, the Euler discretization of Eq. (4) at tm is

Gtm+1=Gtm+κx,tmGtmFG2WiΔtm+1WiΔWtm,
(10)

where Δtm=tmtm1,ΔWtm=WtmWtm1=TMli=tm1+1tmIi,m=1,,M1.

Given a sequence of approximations Gl of a stochastic variable G, we consider Monte Carlo path simulations for a sequence of approximations glG of the configuration vector function g(G) with different time steps Δtl = MlT, M ∈ [2, 8], l = 0, 1, …, L. According to Giles and Michael,12 a trivial telescoping sum gives

gLG=g0G+l=0LglGgl1G
(11)

expressing the expected value on the finest level as the expected value on the coarsest level of approximation plus a sum of expected corrections. Given that the accuracy of Eq. (11) is ε, independent samples Nl are used at each level, and therefore, we can use the following unbiased estimator:

N01n=1N0g00,nG+l=1LNl1n=1Nlgll,nGgl1l,nG,
(12)

where g−1(G) ≡ 0.

Then, a standard piece of theory gives the mean square error (MSE) as

MSE=gLGEgG2=l=0LNl1VglGgl1G+gLGgG2.
(13)

Known from Eq. (13),

gGgLG=l=0LglG=l=0Ln=1Nlgll,nGgl1l,nGNl=l=0LΔgll,nGNl:=¯λML.
(14)

Then, the error ε of g(G) is

ε=gG¯λML=gGgLG+gLG¯λML=εT+εS,
(15)

where εT is the discrete error caused by time dispersion and εS is the statistical error caused by the expectation of Nl samples.

Firstly, to control the time discretization error εT, the different time step on the level l Δtl = MlT, thus we get

gGgLGl=1LφlΔtl2,
(16)

where φl measures the density of the global error. The relationship between the mean number of time steps f¯, φl, and εT is

φlΔtl2=εTf¯l=0,1,2,,L
(17)

if

φlΔtl2εTf¯
(18)

the computation of time step is fulfilled. Otherwise, the time step is refined by splitting it into two equal parts if

φlΔtl2>εTf¯.
(19)

Second, to control the sample statistical error εS, the number of Ns is the optimal samples, which satisfies V[gl(G) − gl−1(G)] ≤ εT/2. On level l, the number of steps is fl and we get

fl+fl1=OεS1+1f¯OεS1=OεS1.
(20)

The optimal sample Ns can be obtained by Lagrange Optimization method in Ref. 13.

Ns=2ε2Vl/Cll=0LVlCl.
(21)

Where Cl and Vl are the calculation cost and calculation variance of gl(G) − gl−1(G) on the l level, respectively. Consequently, the configuration vector function is

gG=l=0LglGl=0Ln=1Nlgll,nGgl1l,nG/Ns.
(22)
  • In the particular case in which gLGgGC12αl, VlC22βl, ClC32γl, the total cost to achieve the ε2 MSE accuracy is

C=Oε2,    β>γOε2log ε12,β=γOε2βγ/α,β<γ,
(23)

where α, β, γ, C1, C2, C3 are positive constants.

  • (b)

    The values of g(G) of shear stress τxy and first normal stress difference τxxτyy in the Hooke model are, respectively,

gSG=1Wipq,gNG=1Wip2q2.
(24)

The corresponding MLMC method’s calculation formulas are

τxy=1Wil=0Ls=0Nspql,sNs,τxxτyy=1Wil=0Ls=0Nspl,s2-ql,s2Ns.
(25)

The values of g(G) corresponding to shear stress and first normal stress difference in the FENE model are, respectively,

gSG=(b+4)pqWi1p2+q2/b,gNG=(b+4)p2q2Wi1p2+q2/b.
(26)

Therefore, the corresponding MLMC method’s calculation formulas are

τxy=(b+4)Wil=0Ls=0Nspql,sNs1pl,s2+ql,s2b,
(27)
τxxτyy=(b+4)Wil=0Ls=0Nspl,s2-ql,s2Ns1pl,s2+ql,s2b.
(28)

In order to optimize the random error or reduce the calculation amount of given error on the basis of the MLMC method, an equilibrium control variable23G¯ is introduced based on the stochastic differential equation (4). G¯ has essentially the same fluctuations, and the average remains unchanged as the random variable, where G¯ satisfies the stochastic differential equation

dG¯x,t=FG¯2Widt+1WidWt.
(29)

gG¯ has essentially the same fluctuations, and the average remains unchanged as the configuration vector function g(G). gG¯ and g(G) are independent Gaussian random variables. Then, the average of the quantity g(G) remains unchanged while the fluctuations gGgG¯ are reduced. The expectance of g(G) is split into two parts,

gG=gG¯+gGgG¯.
(30)

The evolution equation (29) of the Brownian configuration field Hooke model in the component form p¯x,t and q¯x,t is

dp¯x,t=p¯x,t2Widt+1Widwt,
(31)
dq¯x,t=q¯x,t2Widt+1Widvt.
(32)

The evolution equation (29) of the Brownian configuration field FENE model in the component form p¯x,t and q¯x,t is

dpx,t=12Wipx,t1p2x,t+q2x,t/bdt+1Widwt,
(33)
dqx,t=12Wiqx,t1p2x,t+q2x,t/bdt+1Widvt.
(34)

The Euler discretization of Eq. (28) is

G¯tm+1=G¯tmFG¯2WiΔtm+1WiΔWtm.
(35)

Therefore, the values of g(G) of shear stress and first normal stress difference in the Hooke model are, respectively,

gSG=1Wipqp¯q¯,gNG=1Wip2p¯2q2q¯2.
(36)

The corresponding shear stress and first normal stress difference are, respectively,

τxy=1Wil=0Ls=0Nspqp¯q¯l,sNs,τxxτyy=1Wil=0Ls=0Nspl,s2ql,s2p¯l,s2q¯l,s2Ns.
(37)

The FENE model is similar to the Hooke model, given as follows:

gSG=(b+4)pqWi1p2+q2/b(b+4)p¯q¯Wi1p¯2+q¯2/b,
(38)
gNG=(b+4)p2q2Wi1p2+q2/b(b+4)p¯2q¯2Wi1p¯2+q¯2/b.
(39)

The corresponding shear stress and first normal stress difference are, respectively,

τxy=(b+4)Wil=0Ls=0Nspql,sNs1pl,s2+ql,s2bp¯q¯l,sN¯s1p¯l,s2+q¯l,s2b,
(40)
τxxτyy=(b+4)Wil=0Ls=0Nspl,s2ql,s2Ns1pl,s2+ql,s2bp¯l,s2q¯l,s2N¯s1p¯l,s2+q¯l,s2b.
(41)

In this section, in order to further optimize the random error and reduce the total cost, an improved-MLMC method is proposed by combining the MLMC method with the equilibrium control variable method. Our main idea is to reduce the variance of the estimator glG at each level l by using the equilibrium control variable method. When the control variable(equilibrium variable) is subtracted from the variable of the configuration, then the average of the quantity of the configuration vector function remains unchanged. The first term on the right-hand side of Eq. (29) is zero from the definition of the equilibrium configuration, whereas the second one is computed using Brownian dynamics simulations. Thus, glG satisfies

glG=glGglG¯.
(42)

Therefore, the sum of stress at each level is

gG=l=0LglG=l=0LglGglG¯=l=0Ln=1Nlgll,nGgl1l,nG/Nsn=1Nlg¯ll,nGg¯l1l,nG/N¯s.
(43)

A detailed description of the improved-MLMC method is given in Algorithm 1.

ALGORITHM 1

Improved-MLMC method.

Input: M, N, L, N0, ε, Wi, κ11, κ12, κ21, κ22, T, b 
Output: τxy, τxxτyy, C, Vl 
Data initialization 
If l == 0 
while t < T do 
t = tf (fine level time) 
ΔW=ΔtIn 
Gtm+1=Gtm+κx,tGtmFG2WiΔtm+1WiΔWtm 
G¯tm+1=G¯tmFG2WiΔtm+1WiΔWtm 
G1=G,G¯ 
Δt = min(Δt(G1, l), Ttf
tf = min(tf + Δt, T
end if 
else 
while t < T do 
t = tc(coarse level time) 
ΔW=ΔtIn 
Gtm+1=Gtm+κx,tGtmFG2WiΔtm+1WiΔWtm 
G¯tm+1=G¯tmFG2WiΔtm+1WiΔWtm 
G1=G,G¯Δt = min(Δt(G1, l − 1), Ttc
tc = min(tc + Δt, T
end 
Gl=GlG¯l, 
If l > 0 
Gl1=Gl1G¯l1 
End 
glG=glGglG¯ 
Evaluate the additional sample number dNl of gl(G) and glG¯ 
at each level l, ifl=0LdNl<0, output Nl = NL 
Otherwise 
Compute and update the variance and mean Vl, mll = 0, 1, …, 
L at each level 
Estimate α, β, γ 
Compute and update the cost Cl, l = 0, 1, …, L at each level 
Compute and update the optimal number of samples Ns at each 
level 
Weak convergence test 
If not converged, set L: = L + 1 and initialize target NL
end 
gG=l=0LglGglG¯ 
Input: M, N, L, N0, ε, Wi, κ11, κ12, κ21, κ22, T, b 
Output: τxy, τxxτyy, C, Vl 
Data initialization 
If l == 0 
while t < T do 
t = tf (fine level time) 
ΔW=ΔtIn 
Gtm+1=Gtm+κx,tGtmFG2WiΔtm+1WiΔWtm 
G¯tm+1=G¯tmFG2WiΔtm+1WiΔWtm 
G1=G,G¯ 
Δt = min(Δt(G1, l), Ttf
tf = min(tf + Δt, T
end if 
else 
while t < T do 
t = tc(coarse level time) 
ΔW=ΔtIn 
Gtm+1=Gtm+κx,tGtmFG2WiΔtm+1WiΔWtm 
G¯tm+1=G¯tmFG2WiΔtm+1WiΔWtm 
G1=G,G¯Δt = min(Δt(G1, l − 1), Ttc
tc = min(tc + Δt, T
end 
Gl=GlG¯l, 
If l > 0 
Gl1=Gl1G¯l1 
End 
glG=glGglG¯ 
Evaluate the additional sample number dNl of gl(G) and glG¯ 
at each level l, ifl=0LdNl<0, output Nl = NL 
Otherwise 
Compute and update the variance and mean Vl, mll = 0, 1, …, 
L at each level 
Estimate α, β, γ 
Compute and update the cost Cl, l = 0, 1, …, L at each level 
Compute and update the optimal number of samples Ns at each 
level 
Weak convergence test 
If not converged, set L: = L + 1 and initialize target NL
end 
gG=l=0LglGglG¯ 

To demonstrate the fidelity and computational efficiency of our algorithm, the results of benchmark flow problems, namely, Hooke model in a simple shear flow field, are presented first. The implementation of our algorithms are realized in MATLAB, and the details code can be referred to a recent review paper,12 shear stress and first normal stress difference are computed in steady state, and comparing them with the analytical solution. The steady state analytical solution of the Hooke model is

τxy=κ12,τxxτyy=2Wiκ122.

In the numerical simulations, the values of the physical parameters are set, respectively, as follows: initial sample N0 = 20, accuracy ε = 0.005, convergence test sample N = 1000, Wi = 1, time step Δt = 2l, and l = 0, 1, …, L. The four components of the velocity gradient transpose are, respectively, κ11 = 0, κ12 = −0.5, κ21 = 0, and κ22 = 0. Figure 1 shows the comparisons between the numerical result and the steady analytical solution of stress evolution with time by the Hooke model, presenting both shear stress and first normal stress difference.

FIG. 1.

Comparison between the numerical result and the steady analytical solution of stress evolution with time by the Hooke model. (a) Shear stress. (b) First normal stress difference.

FIG. 1.

Comparison between the numerical result and the steady analytical solution of stress evolution with time by the Hooke model. (a) Shear stress. (b) First normal stress difference.

Close modal

As evidenced in Fig. 1, the numerical solution from the improved-MLMC method matches the analytical solution, thus demonstrating the stress fields from the presented MLMC method have been solved accurately.

Figures 2 and 3 further show the influence of the accuracy, ε, on the evolution of the residual errors of shear stress and first normal stress difference. The residual errors Err(τxy) and Err(τxxτyy) of shear stress and first normal stress difference are, respectively,

Errτxy=pqpqLp¯q¯p¯q¯L=l=L+1pqlpql1p¯q¯l+p¯q¯l1,
Errτxxτyy=ppqqppqqLp¯p¯q¯q¯p¯p¯q¯q¯L=L=L+1ppqqlppqql1p¯p¯q¯q¯l+p¯p¯q¯q¯l1.

As shown in Fig. 3, the residual errors decrease asymptotically with increasing accuracy. Thus, the numerical solution is closer to the analytical solution.

FIG. 2.

Evolution of the residual error of shear stress for the Hooke model.

FIG. 2.

Evolution of the residual error of shear stress for the Hooke model.

Close modal
FIG. 3.

Evolution of the residual error of the first normal shear stress flow difference for the Hooke model.

FIG. 3.

Evolution of the residual error of the first normal shear stress flow difference for the Hooke model.

Close modal

In this section, in order to check the efficiency of our algorithm, we have studied the total cost of shear stress and first normal stress difference in both shear flow and extensional flow fields, with both Hooke and FENE models, and compared it with StdMC and MLMC methods.

As the examples of the applications of the improved-MLMC method developed in this paper to a shear flow field situation, we have tested the stress performance of the two models, namely, the Hooke and FENE models. The values of the physical parameters are κ11 = 0, κ12 = −0.5, κ21 = 0, κ22 = 0, T = 1, N0 = 20, Δt = 2l, and l = 0, 1, …, L. The dimensionless finite extensibility parameter in the FENE model is b = 50. Other parameter selections are shown in Table I.

TABLE I.

Parameter selection of two models in the shear flow field.

εN
Hooke model 0.005 1000 
FENE model 0.0005 9000 
εN
Hooke model 0.005 1000 
FENE model 0.0005 9000 

Figures 4(a) and 4(b) depict the comparison of total cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in Hooke and FENE models at Wi = 1.5 and 2, respectively. From the figure, one can clearly observe that the total cost (Definition is similar to Ref. 13.) of the presented MLMC method is almost lower than that of the MLMC and StdMC methods at any given accuracy. In addition, in order to further illustrate the computational efficiency of the improved-MLMC method, the total cost reduction rate ϕ(%) of shear stress with varying accuracy is calculated in Tables II and III. Specifically, ϕ(%) is, respectively, evaluated as follows:

ϕ%MLMC=CMLMCCimprovedMLMCCMLMC×100%,
(44)
ϕ%StdMC=CStdMCCimprovedMLMCCStdMC×100%.
(45)
FIG. 4.

Comparison of computational cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in the shear flow field: (a) shear stress in the Hooke model; (b) shear stress in the FENE model.

FIG. 4.

Comparison of computational cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in the shear flow field: (a) shear stress in the Hooke model; (b) shear stress in the FENE model.

Close modal
TABLE II.

Total cost reduction rate of shear stress in the shear flow field in the Hooke model.

ε0.0010.0020.0030.0040.005
ϕ(%)MLMC 27.72 27.02 10.96 37.23 29.41 
ϕ(%)StdMC 99.85 99.94 99.79 99.95 99.83 
ε0.0010.0020.0030.0040.005
ϕ(%)MLMC 27.72 27.02 10.96 37.23 29.41 
ϕ(%)StdMC 99.85 99.94 99.79 99.95 99.83 
TABLE III.

Total cost reduction rate of shear stress in the shear flow field in the FENE model.

ε0.00010.00020.00030.00040.0005
ϕ(%)MLMC 42.53 42.65 30.40 30.59 18.01 
ϕ(%)StdMC 67.10 67.10 60.08 60.07 52.40 
ε0.00010.00020.00030.00040.0005
ϕ(%)MLMC 42.53 42.65 30.40 30.59 18.01 
ϕ(%)StdMC 67.10 67.10 60.08 60.07 52.40 

However, it should be noted that the shear stress based on the Hooke model performs better than that based on the FENE model. Tables II and III illustrate that the total cost reduction of the Hooke model is more than that of the FENE model, and the differences are more pronounced when the accuracy is decreased.

We shall focus on the total cost of the extensional flow field. Two different models have been examined, namely, Hooke and FENE models. The values of physical parameters when calculating first normal stress differences of the FENE model are κ12 = 2, κ22 = −1. The values of the physical parameters in the Hooke model are κ11 = 0.5, κ12 = 0, κ21 = 0, and κ22 = −0.5. Other parameters are consistent with those in the shear flow field.

Figures 5(a) and 5(b) present the comparison of total cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in Hooke and FENE models at Wi = 1.5 and 2, respectively. Clearly, over the range of accuracies considered, the total costs of the presented MLMC method are also much lower as compared to those of other methods. Furthermore, we still illustrate the efficiency of this method by calculating the total cost reduction rate. From Tables IV and V, obviously, it is seen that this method works well, since the total cost reduction rate is reduced. Generally speaking, it is demonstrated that the improved-MLMC method is the preferred choice for the Hooke model, while it is more suited for the problem of high accuracy in the FENE model.

FIG. 5.

Comparison of computational cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in the extensional flow field: first normal stress difference in the (a) Hooke model and (b) FENE model.

FIG. 5.

Comparison of computational cost with varying accuracy among StdMC, MLMC, and improved-MLMC methods in the extensional flow field: first normal stress difference in the (a) Hooke model and (b) FENE model.

Close modal
TABLE IV.

Total cost reduction rate of first normal stress differences in the extensional flow field in the Hooke model.

ε0.0010.0020.0030.0040.005
ϕ(%)MLMC 16.68 14.27 14.90 17.12 16.59 
ϕ(%)StdMC 70.89 70.22 70.47 71.00 70.57 
ε0.0010.0020.0030.0040.005
ϕ(%)MLMC 16.68 14.27 14.90 17.12 16.59 
ϕ(%)StdMC 70.89 70.22 70.47 71.00 70.57 
TABLE V.

Total cost reduction rate of first normal stress differences in the extensional flow field in the FENE model.

ε0.00010.00030.00050.00070.001
ϕ(%)MLMC 25.75 49.16 68.37 24.13 59.63 
ϕ(%)StdMC 99.98 99.40 99.60 99.20 99.41 
ε0.00010.00030.00050.00070.001
ϕ(%)MLMC 25.75 49.16 68.37 24.13 59.63 
ϕ(%)StdMC 99.98 99.40 99.60 99.20 99.41 

It is well known that the Wi number is an important parameter to describe the effectiveness and adaptability of polymer fluids. This section discusses effects of the Wi number on the total cost of the improved-MLMC method. In the two models, the values of the physical parameters are ε = 0.005, N = 9000, N0 = 2, and b = 40. Other parameters are consistent with those in the shear flow field and extensional flow field.

Figure 6 shows the total cost in shear flow and extensional flow fields of the Hooke model as a function of the Wi number using StdMC, MLMC, and improved-MLMC methods, respectively. As can be seen from Fig. 6, we find the total costs of shear stress and first normal stress difference both decrease when the Wi number is increased for all the three types of methods, StdMC, MLMC, and improved-MLMC methods. As we expected, the total cost of the presented MLMC method is also much lower as compared to that of StdMC and MLMC methods. At this time, we note that the total cost reduction rate decreases with the increase in the Wi number, that is, the computational efficiency decreases.

FIG. 6.

Total cost(C) in the Hooke model as a function of the Wi number: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

FIG. 6.

Total cost(C) in the Hooke model as a function of the Wi number: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

Close modal

We have followed an analogous research method with the FENE model. The results of total cost in shear flow and extensional flow fields of the FENE model as a function of the Wi number are presented in Fig. 7 using StdMC, MLMC, and improved-MLMC methods, respectively. As seen from Fig. 7, the total cost in two kinds of stresses in different flow fields decreases with the increase in the Wi number; the improvement would be less dramatic at higher Wi numbers as the total cost reduction rate required would decrease, and it may be impossible to compute the coarse level total cost exactly.

FIG. 7.

Total cost(C) in the FENE dumbbell model as a function of the Wi number: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

FIG. 7.

Total cost(C) in the FENE dumbbell model as a function of the Wi number: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

Close modal

In order to illustrate the relationship between the total cost of the improved-MLMC method and the Wi number, we compute Vlg(G)) with the parameter set: Wi = 0.5, 2, and 5. As evidenced in Fig. 8(a) of the Hooke model, increasing the Wi number from 0.5 to 5 leads to the decrease in the level where variance reduction is not obvious. In Fig. 8(b) of the Hooke model, when l < 6, Vlg(G)) decreases with the number of levels, and in case of l = 6, the Wi number is larger and Vlg(G)) is smaller. In general, at larger Wi numbers, Vlg(G)) decreases faster and reaches the minimum. At this time, increasing the number of levels can no longer reduce the variance.

FIG. 8.

Variance Vlg(G)) as a function of the level in the Hooke model: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

FIG. 8.

Variance Vlg(G)) as a function of the level in the Hooke model: (a) shear stress in the shear flow field; (b) first normal stress difference in the extensional flow field.

Close modal

In this paper, an improved multilevel Monte Carlo (improved-MLMC) method is introduced to solve the problems of high total cost in polymer flows with the low Wi number for the traditional BCF method. First, the numerical solutions have been verified by comparing them with the analytical solution from the Hooke model. The empirical study demonstrates that this method works well. Second, the total costs from the improved-MLMC method have been validated by comparing them with the numerical solution from StdMC and MLMC methods. This can be explained based on the fact that the efficiency of the presented method is optimized not only by the MLMC method but also by equilibrium control variables at each level. Finally, we discussed the effect of the Wi number on the total cost. The dominant cost lies on the coarsest level for low Wi numbers, but with the increase in the Wi number, the main cost lies on the finest level. Therefore, the presented MLMC method is more efficient in low Wi number problems. Based on these results, we conclude that for the simple model, one new improved-MLMC method is proposed, which is much superior to commonly used StdMC and MLMC methods, especially in regard to the total cost of low Wi number and high accuracy.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the National Natural Science Foundation of China (Grant No. 11601411) and Scientific Research Program funded by Shaanxi Provincial Education Department (Program No. 20JY025).

1.
H. C.
Öttinger
,
J. Polym. Int.
43
,
290
(
1997
).
2.
M.
Laso
and
H. C.
Öttinger
,
J. Non-Newtonian Fluid Mech.
47
,
1
(
1993
).
3.
M. A.
Hulsen
,
A. P. G.
van Heel
, and
B. H. A. A.
van den Brule
,
J. Non-Newtonian Fluid Mech.
70
,
79
(
1997
).
4.
H. C.
Öttinger
,
B. H. A. A.
van den Brule
, and
M. A.
Hulsen
,
J. Non-Newtonian Fluid Mech.
70
,
255
(
1997
).
5.
M.
Melchior
and
H. C.
Öttinger
,
J. Chem. Phys.
103
,
9506
(
1995
).
6.
P.
Halin
,
G.
Lielens
,
R.
Keunings
, and
V.
Legat
,
J. Non-Newtonian Fluid Mech.
79
,
387
(
1998
).
7.
J.
Bonvin
and
M.
Picasso
,
J. Non-Newtonian Fluid Mech.
84
,
191
(
1999
).
8.
C.
Chauvière
,
SIAM J. Sci. Comput.
23
,
2123
(
2002
).
9.
C.
Chauvière
and
A.
Lozinski
,
SIAM J. Sci. Comput.
24
,
1823
(
2002
).
10.
S.
Heinrich
,
Lect. Notes Comput. Sci.
2179
,
58
(
2001
).
11.
A.
Brandt
and
V.
Ilyin
,
J. Mol. Liq.
105
,
245
(
2003
).
12.
M. B.
Giles
and
B.
Michael
,
Acta Numer.
24
,
259
(
2015
).
13.
M.
Giles
and
L.
Szpruch
,
Recent Developments in Computational Finance
, edited by
T.
Gerstner
and
P.
Kloeden
(
World Scientific
,
2013
), p.
3
47
.
14.
K.
Prashant
,
L.
Peiyao
,
J. G.
Francisco
, and
W. O.
Cornelis
,
J. Comput. Phys.
371
,
382
(
2018
).
15.
D. F.
Anderson
and
D. J.
Higham
,
SIAM Multiscale Modell. Simul.
10
,
146
(
2012
).
16.
S.
Pauli
and
P.
Arbenz
,
Comput. Math. Appl.
70
,
2638
(
2015
).
17.
M. B.
Giles
,
T.
Nagapetyan
, and
K.
Ritter
,
SIAM J. Uncertainty Quantif.
3
,
267
(
2015
).
18.
E.
Ullmann
and
I.
Papaioannou
,
SIAM J. Uncertainty Quantif.
3
,
922
(
2015
).
19.
M. B.
Giles
,
C.
Lester
, and
J.
Whittle
,
Monte Carlo and Quasi-Monte Carlo Methods
edited by
R.
Cools
and
D.
Nuyens
(
Springer
,
2016
), pp.
303
314
.
20.
D.
Belomestny
and
T.
Nagapetyan
,
Bernoulli
23
,
927
(
2017
).
21.
H. C.
Öttinger
,
Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing Simulation Algorithms
(
Springer
,
1996
), p.
1158
.
22.
P. E.
Kloeden
and
E.
Platen
,
Numerical Solution of Stochastic Differential Equations
(
Springer
,
1992
), pp.
305
307
.
23.
M.
Melchior
and
H. C.
Öttinger
,
J. Chem. Phys.
105
,
3316
(
1996
).