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.

## I. INTRODUCTION

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, Öttinger^{1} 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 Öttinger^{2} termed this hybrid method CONNFFESSIT (Calculation of Non-Newtonian Flows: Finite Element and Stochastic Simulation Technique). Later, the BCF (Brownian Configuration Field) techniques^{3} 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 Picasso^{7} extended these ideas to nonequilibrium ensembles of dumbbells. Chauvière^{8} 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 method^{10–12} has been used successfully in a wide range of applications from mathematical finance^{13} to uncertainty quantification in subsurface flows.^{14} Anderson and Higham^{15} showed how to extend the multilevel Monte Carlo approach to the continuous time Markov chain. Pauli and Arbenz^{16} 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 Papaioannou^{18} 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 Nagapetyan^{20} 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.

## II. THE STOCHASTIC MODEL OF POLYMERS IN SIMPLE FLOW FIELD

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

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

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

The evolution equation^{21} for the configuration vector field is given by

where *ζ* is the friction coefficient of the dumbbells’ beads in the solvent, *K*_{B} is the Boltzmann constant, *Z* is the absolute temperature, ** F**(

**) is the spring force,**

*G**x*is the vertical abscissa,

*t*is the current time,

**(**

*W**t*) is a multi-dimensional Wiener process, and the components of d

**(**

*W**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

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

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:

Similarly, the component form of the FENE model is

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

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 $G\u2297FG$ denotes the

**⊗**

*G***(**

*F***) mathematical expectation.**

*G*## III. IMPROVED-MLMC METHOD BASED ON EQUILIBRIUM CONTROL VARIABLES

### A. MLMC method

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 =

*t*

_{0}<

*t*

_{1}< $\cdots $ <

*t*

_{m}< $\cdots $ <

*t*

_{M}=

*T*in time [0,

*T*]. Then, the Euler discretization of Eq. (4) at

*t*

_{m}is

where $\Delta tm=tm\u2212tm\u22121,\Delta Wtm=Wtm\u2212Wtm\u22121=TMl\u2211i=tm\u22121+1tmIi,m=1,\u2026,M\u22121$.

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*(

**) with different time steps Δ**

*G**t*

_{l}=

*M*

^{−l}

*T*,

*M*∈ [2, 8],

*l*= 0, 1, …,

*L*. According to Giles and Michael,

^{12}a trivial telescoping sum gives

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 *N*_{l} are used at each level, and therefore, we can use the following unbiased estimator:

where *g*_{−1}(** G**) ≡ 0.

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

Known from Eq. (13),

Then, the error *ε* of *g*(** G**) is

where *ε*_{T} is the discrete error caused by time dispersion and *ε*_{S} is the statistical error caused by the expectation of *N*_{l} samples.

Firstly, to control the time discretization error *ε*_{T}, the different time step on the level *l* Δ*t*_{l} = *M*^{l}*T*, thus we get

where *φ*_{l} measures the density of the global error. The relationship between the mean number of time steps $f\xaf$, *φ*_{l}, and *ε*_{T} is

if

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

Second, to control the sample statistical error *ε*_{S}, the number of *N*_{s} is the optimal samples, which satisfies *V*[*g*_{l}(** G**) −

*g*

_{l−1}(

**)] ≤**

*G**ε*

_{T}/2. On level

*l*, the number of steps is

*f*

^{l}and we get

The optimal sample *N*_{s} can be obtained by Lagrange Optimization method in Ref. 13.

Where *C*_{l} and *V*_{l} are the calculation cost and calculation variance of *g*_{l}($G$) − *g*_{l−1}($G$) on the *l* level, respectively. Consequently, the configuration vector function is

In the particular case in which $gLG\u2212gG\u2264C12\u2212\alpha l$,

*V*_{l}≤*C*_{2}2^{−βl},*C*_{l}≤*C*_{3}2^{γl}, the total cost to achieve the*ε*^{2}MSE accuracy is

where *α*, *β*, *γ*, *C*_{1}, *C*_{2}, *C*_{3} are positive constants.

- (b)
The values of

*g*() of shear stress*G**τ*_{xy}and first normal stress difference*τ*_{xx}−*τ*_{yy}in the Hooke model are, respectively,

The corresponding MLMC method’s calculation formulas are

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

Therefore, the corresponding MLMC method’s calculation formulas are

### B. Variance reduction technique based on equilibrium control variable method

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 variable^{23} $G\xaf$ is introduced based on the stochastic differential equation (4). $G\xaf$ has essentially the same fluctuations, and the average remains unchanged as the random variable, where $G\xaf$ satisfies the stochastic differential equation

$gG\xaf$ has essentially the same fluctuations, and the average remains unchanged as the configuration vector function ** g**(

**). $gG\xaf$ and**

*G***(**

*g***) are independent Gaussian random variables. Then, the average of the quantity**

*G***(**

*g***) remains unchanged while the fluctuations $gG\u2212gG\xaf$ are reduced. The expectance of**

*G***(**

*g***) is split into two parts,**

*G*The evolution equation (29) of the Brownian configuration field Hooke model in the component form $p\xafx,t$ and $q\xafx,t$ is

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

The Euler discretization of Eq. (28) is

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

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

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

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

### C. Improved-MLMC method

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

Therefore, the sum of stress at each level is

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

Input: M, N, L, N_{0}, ε, W_{i}, κ_{11}, κ_{12}, κ_{21}, κ_{22}, T, b |

Output: τ_{xy}, τ_{xx} − τ_{yy}, C, V_{l} |

Data initialization |

If l == 0 |

while t < T do |

t = t^{f} (fine level time) |

$\Delta W=\Delta tIn$ |

$Gtm+1=Gtm+\kappa x,t\u22c5Gtm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G\xaftm+1=G\xaftm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G1=G,G\xaf$ |

Δt = min(Δt(1, Gl), T − t^{f}) |

t^{f} = min(t^{f} + Δt, T) |

end if |

else |

while t < T do |

t = t^{c}(coarse level time) |

$\Delta W=\Delta tIn$ |

$Gtm+1=Gtm+\kappa x,t\u22c5Gtm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G\xaftm+1=G\xaftm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G1=G,G\xaf$Δt = min(Δt(1, Gl − 1), T − t^{c}) |

t^{c} = min(t^{c} + Δt, T) |

end |

$Gl=Gl\u2212G\xafl$_{,} |

If l > 0 |

$Gl\u22121=Gl\u22121\u2212G\xafl\u22121$ |

End |

$glG=glG\u2212glG\xaf$ |

Evaluate the additional sample number dN_{l} of g_{l}() and $glG\xaf$ G |

at each level l, if$\u2211l=0LdNl<0$, output N_{l} = N_{L} |

Otherwise |

Compute and update the variance and mean V_{l}, m_{l}l = 0, 1, …, |

L at each level |

Estimate α, β, γ |

Compute and update the cost C_{l}, l = 0, 1, …, L at each level |

Compute and update the optimal number of samples N_{s} at each |

level |

Weak convergence test |

If not converged, set L: = L + 1 and initialize target N_{L}. |

end |

$gG=\u2211l=0LglG\u2212glG\xaf$ |

Input: M, N, L, N_{0}, ε, W_{i}, κ_{11}, κ_{12}, κ_{21}, κ_{22}, T, b |

Output: τ_{xy}, τ_{xx} − τ_{yy}, C, V_{l} |

Data initialization |

If l == 0 |

while t < T do |

t = t^{f} (fine level time) |

$\Delta W=\Delta tIn$ |

$Gtm+1=Gtm+\kappa x,t\u22c5Gtm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G\xaftm+1=G\xaftm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G1=G,G\xaf$ |

Δt = min(Δt(1, Gl), T − t^{f}) |

t^{f} = min(t^{f} + Δt, T) |

end if |

else |

while t < T do |

t = t^{c}(coarse level time) |

$\Delta W=\Delta tIn$ |

$Gtm+1=Gtm+\kappa x,t\u22c5Gtm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G\xaftm+1=G\xaftm\u2212FG2Wi\Delta tm+1Wi\Delta Wtm$ |

$G1=G,G\xaf$Δt = min(Δt(1, Gl − 1), T − t^{c}) |

t^{c} = min(t^{c} + Δt, T) |

end |

$Gl=Gl\u2212G\xafl$_{,} |

If l > 0 |

$Gl\u22121=Gl\u22121\u2212G\xafl\u22121$ |

End |

$glG=glG\u2212glG\xaf$ |

Evaluate the additional sample number dN_{l} of g_{l}() and $glG\xaf$ G |

at each level l, if$\u2211l=0LdNl<0$, output N_{l} = N_{L} |

Otherwise |

Compute and update the variance and mean V_{l}, m_{l}l = 0, 1, …, |

L at each level |

Estimate α, β, γ |

Compute and update the cost C_{l}, l = 0, 1, …, L at each level |

Compute and update the optimal number of samples N_{s} at each |

level |

Weak convergence test |

If not converged, set L: = L + 1 and initialize target N_{L}. |

end |

$gG=\u2211l=0LglG\u2212glG\xaf$ |

## IV. NUMERICAL RESULTS

### A. The benchmark test

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

In the numerical simulations, the values of the physical parameters are set, respectively, as follows: initial sample *N*_{0} = 20, accuracy *ε* = 0.005, convergence test sample *N* = 1000, *Wi* = 1, time step Δ*t* = 2^{−l}, 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.

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,

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

### B. The efficiency of improved-MLMC method in a simple flow field

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, *N*_{0} = 20, Δ*t* = 2^{−l}, 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.

. | $\epsilon $ . | N
. |
---|---|---|

Hooke model | 0.005 | 1000 |

FENE model | 0.0005 | 9000 |

. | $\epsilon $ . | 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:

$\epsilon $ . | 0.001 . | 0.002 . | 0.003 . | 0.004 . | 0.005 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 27.72 | 27.02 | 10.96 | 37.23 | 29.41 |

ϕ(%)_{StdMC} | 99.85 | 99.94 | 99.79 | 99.95 | 99.83 |

$\epsilon $ . | 0.001 . | 0.002 . | 0.003 . | 0.004 . | 0.005 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 27.72 | 27.02 | 10.96 | 37.23 | 29.41 |

ϕ(%)_{StdMC} | 99.85 | 99.94 | 99.79 | 99.95 | 99.83 |

$\epsilon $ . | 0.0001 . | 0.0002 . | 0.0003 . | 0.0004 . | 0.0005 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 42.53 | 42.65 | 30.40 | 30.59 | 18.01 |

ϕ(%)_{StdMC} | 67.10 | 67.10 | 60.08 | 60.07 | 52.40 |

$\epsilon $ . | 0.0001 . | 0.0002 . | 0.0003 . | 0.0004 . | 0.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.

$\epsilon $ . | 0.001 . | 0.002 . | 0.003 . | 0.004 . | 0.005 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 16.68 | 14.27 | 14.90 | 17.12 | 16.59 |

ϕ(%)_{StdMC} | 70.89 | 70.22 | 70.47 | 71.00 | 70.57 |

$\epsilon $ . | 0.001 . | 0.002 . | 0.003 . | 0.004 . | 0.005 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 16.68 | 14.27 | 14.90 | 17.12 | 16.59 |

ϕ(%)_{StdMC} | 70.89 | 70.22 | 70.47 | 71.00 | 70.57 |

$\epsilon $ . | 0.0001 . | 0.0003 . | 0.0005 . | 0.0007 . | 0.001 . |
---|---|---|---|---|---|

ϕ(%)_{MLMC} | 25.75 | 49.16 | 68.37 | 24.13 | 59.63 |

ϕ(%)_{StdMC} | 99.98 | 99.40 | 99.60 | 99.20 | 99.41 |

$\epsilon $ . | 0.0001 . | 0.0003 . | 0.0005 . | 0.0007 . | 0.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, *N*_{0} = 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.

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.

In order to illustrate the relationship between the total cost of the improved-MLMC method and the *Wi* number, we compute *V*_{l}(Δ*g*(** 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,

*V*

_{l}(Δ

*g*(

**)) decreases with the number of levels, and in case of**

*G**l*= 6, the

*Wi*number is larger and

*V*

_{l}(Δ

*g*(

**)) is smaller. In general, at larger**

*G**Wi*numbers,

*V*

_{l}(Δ

*g*(

**)) decreases faster and reaches the minimum. At this time, increasing the number of levels can no longer reduce the variance.**

*G*## V. CONCLUSIONS AND DISCUSSION

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.

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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).