Multilevel Monte Carlo method for the Brownian configuration field of polymer fluids Open Access

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¹ 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² termed this hybrid method CONNFFESSIT (Calculation of Non-Newtonian Flows: Finite Element and Stochastic Simulation Technique). Later, the BCF (Brownian Configuration Field) techniques³ 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.,³ Öttinger et al.,⁴ and Halin et al.⁶ 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⁷ extended these ideas to nonequilibrium ensembles of dumbbells. Chauvière⁸ 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.⁹ 

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¹³ to uncertainty quantification in subsurface flows.¹⁴ Anderson and Higham¹⁵ showed how to extend the multilevel Monte Carlo approach to the continuous time Markov chain. Pauli and Arbenz¹⁶ investigated a fault tolerant MLMC method; the proposed optimization procedure can react on experienced faults. Giles et al.¹⁷ constructed and analyzed the multilevel Monte Carlo method for the approximation of distribution functions and densities of univariate random variables. Ullmann and Papaioannou¹⁸ 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.¹⁹ 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²⁰ 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

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²¹ 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(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

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⊗FG$ denotes the G ⊗ F(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 = t₀ < t₁ < $⋯$ < t_m < $⋯$ < t_M = T in time [0, T]. Then, the Euler discretization of Eq. (4) at t_m is

where $Δtm=tm−tm−1,ΔWtm=Wtm−Wtm−1=TMl∑i=tm−1+1tmIi,m=1,…,M−1$⁠.

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 Δt_l = M^−lT, M ∈ [2, 8], l = 0, 1, …, L. According to Giles and Michael,¹² 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₋₁(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^lT, thus we get

where φ_l measures the density of the global error. The relationship between the mean number of time steps $f¯$⁠, φ_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−gG≤C12−αl$⁠, V_l ≤ C₂2^−βl, C_l ≤ C₃2^γl, the total cost to achieve the ε² MSE accuracy is

where α, β, γ, C₁, C₂, C₃ 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,

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

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²³ $G¯$ 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

$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 $gG−gG¯$ are reduced. The expectance of g(G) is split into two parts,

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

The evolution equation (29) of the Brownian configuration field FENE model in the component form $p¯x,t$ and $q¯x,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,

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.

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,¹² 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₀ = 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, κ₁₁ = 0, κ₁₂ = −0.5, κ₂₁ = 0, and κ₂₂ = 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.

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 κ₁₁ = 0, κ₁₂ = −0.5, κ₂₁ = 0, κ₂₂ = 0, T = 1, N₀ = 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.

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:

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 κ₁₂ = 2, κ₂₂ = −1. The values of the physical parameters in the Hooke model are κ₁₁ = 0.5, κ₁₂ = 0, κ₂₁ = 0, and κ₂₂ = −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.

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₀ = 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(G)) decreases with the number of levels, and in case of l = 6, the Wi number is larger and V_l(Δg(G)) is smaller. In general, at larger Wi numbers, V_l(Δg(G)) decreases faster and reaches the minimum. At this time, increasing the number of levels can no longer reduce the variance.

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