A Lattice Boltzmann Method and Asynchronous Model Coupling for Viscoelastic Fluids

The numerical algorithms of viscoelastic flows can appear a tremendous challenge as the Weissenberg number (Wi) enlarged sufficiently. In this study, we present a generalized technique of time-stably advancing based on the coupled lattice Boltzmann method, in order to improve the numerical stability of simulations at a high Wi number. The mathematical models of viscoelastic fluids include both the equation of the solvent and the Oldroyd-B constitutive equation of the polymer. In the two-dimensional (2D) channel flow, the coupled method shows good agreements between the corresponding exact results and the numerical results obtained by our method. In addition, as the Wi number increased, for the viscoelastic flows through contractions, we show that the prediction of our presented method can reproduce the same numerical results that were reported by previous studies. The main advantage of current method is that it can be applied to simulate the complex phenomena of the viscoelastic fluids.


Introduction
Over the last few decades, there has been substantial development in the numerical simulation of viscoelastic fluids, which are commonly encountered in many important industries, such as oil and chemical engineering.However, the numerical solution of the viscoelastic constitutive models is a tremendous challenge, due to breakdown in convergence of algorithms when the Weissenberg number (Wi, is a non-dimensional number characterizing the effect of elasticity on the flow) is increased [1,2].At large levels of elasticity, the stress field often contains stress singularities at sharp corners or stagnation points [3].Many researchers have already had some attempts to develop robust and stable numerical algorithm to solve viscoelastic fluid flows at moderately high values of Wi number, such as discontinuous Galerkin method [4], finite-volume methods [5][6][7], and the lattice Boltzmann method.Among the various numerical methods of the viscoelastic fluids in literatures, the lattice Boltzmann method (LBM) [8][9][10] has attracted much attention due to its kinetic nature for modeling the complex fluids.
Since 1990s, Qian [11] and Giraud [12] have been developed the lattice Boltzmann method to account for non-Newtonian behaviors, but they did not focus on the model with strong elastic effects.Later, several approaches have been proposed to incorporate constitutive equations in lattice Boltzmann simulations for viscoelastic fluids [13][14][15][16][17]. Recently, Papenkort et al. [18] investigated viscoelastic shear-thinning fluid using lattice Boltzmann method; Wang et al. [19] employed a multiple-relaxation-time lattice Boltzmann flux solver for non-Newtonian power-law fluid flows.Xie et al. [20] extended a lattice Boltzmann method scheme to multiphase viscoplastic fluids by applying the Herschel-Bulkley constitutive relationship.It is worth mentioning that a coupled lattice Boltzmann method have been devised to solve the equations of an Oldroyd-B viscoelastic fluid in the reference [21], which did not need to adding extra distribution functions for the stress tensor.The coupled systems of the viscoelastic fluid include the macro scale constitutive equation and the mesoscopic lattice Boltzmann equation.There is a more disparity in timescales between the macro solver and meso solver.To distinguish this different time scale and enhance the computational efficiency, the asynchronously coupled method [22,23] was applied for the coupling time in reference [21].That means, the data between macro constitutive equation and the lattice Boltzmann equation exchange at every time step.The constitutive solver may run at a slower pace than the meso-scale dynamics.
The asynchronous coupling method allows for the coupled solvers to exchange the data at every loop, although the meso and macro time steps are unequal.This special treatment is actually a numerical approximation to the true solution.For viscoelastic fluids at low Wi number, asynchronously coupled method may get steady result in simulations if the macro time step can guarantee the adequate resolutions of each solvers because of the small scale-separation.However, for viscoelastic fluids at high Wi number, the computational stability of asynchronously coupling is greatly affected by the polymer stress as the elasticity interaction becomes relatively large.Therefore, the asynchronously coupling method has the limitations, such as poor numerical stability, which were limited to solve the low Wi number problem.
The aim of the current paper is to develop a generalized technique for time advancing of the viscoelastic coupled lattice Boltzmann method, seeking to be immune to high Wi number problem.The basic idea of our new method is as follows.(i) Run the lattice Boltzmann solver using its own time step δt, while allowing for a flexible number of lattice Boltzmann solver steps N before the exchange of coupling variables; (ii) Run the constitutive solver using its own time step ∆t.(iii) Exchange data between the lattice Boltzmann and constitutive solvers.Note that, the number of flexible steps, N, provides an additional control of the coupling.This skill enables the applicability of the coupled lattice Boltzmann method at high Wi number, while minimizing the number of lattice Boltzmann and constitutive time steps needed.We refer to this method as a generalized asynchronously coupling method.This article is organized as follows: Section 2 describes the mathematical models that are used in the simulation of viscoelastic fluids.In Section 3, we give the coupled lattice Boltzmann scheme of an Oldroyd-B viscoelastic fluid and the details of the generalized asynchronously coupling method.In Section 4, we show the results of simulations for viscoelastic fluids at high Wi number using the generalized method.Finally, we give the conclusions of this paper and some future prospects.

The Descriptions of Mathematical Models
In the simulations, we use the Oldroyd-B constitutive model to describe the elasticity effect of viscoelastic fluids, in which the polymer molecules are modeled by the two beads connected by a spring [24].We treat the solvent of the viscoelastic model as an incompressible flow at low Reynolds number.The equations of the momentum and mass conservation are given by where u, ρ 0 , p, η s represent, respectively, the velocity, the density, pressure, and solvent viscosity.τ is the extra-stress tensor accounting for the viscoelastic effects of the stretching and the relaxation of the polymer macro-molecules diluted in the solvent.The constitutive equation of the Oldroyd-B model reads where the Dτ/Dt, the superscript "T" and d = 1 2 (∇u + ∇u T ) denote, respectively, the material derivative, transpose operation and the rate of strain tensor; λ p , η p are the relaxation time of polymer solutions and the dynamic viscosity of the polymer, respectively.In the simulations, we present all the data in non-dimensional forms, normalized in terms of the characteristic length L and the characteristic velocity U (characteristic time is L/U).Weissenberg (Wi) number is defined by Wi = λ p U/L to characterize the viscoelastic effect for shear flows.The Reynolds number, Re, is defined as Re = ρ 0 UL/η t , respectively.The parameter β = η s /η t is the ratio of the solvent viscosity to total viscosity, η t = η s + η p is the total viscosity.Finally, we get

The Coupled LBM of Viscoelastic Fluids
We will briefly describe the numerical implementations simulating the viscoelastic fluids in this section.The two-dimensional (2D) incompressible lattice Boltzmann method is utilized to solve the equations of the solvent.In order to avoid adding a new set of distribution functions, the constitutive equation is split into the relaxation-stretching and advection operators, where the advection operator of the stress tensor is directly solved by implementing the particle distribution functions of solvent.Specifically, to seek the coupled LBM that can be immune to high Wi number problem, we develop a generalized asynchronously coupled method of advancing the time step.All of these details will be described as follow.

Incompressible LBGK Model with the Extra Force
In this study, we use the 2D incompressible Lattice Boltzmann Bhatnagar-Gross-Krook (LBGK) model [25] based on the two-dimensional nine velocity (D2Q9) model to recover the Navier-Stokes equations.The particle velocity e α may be written as here, δx, δt are, respectively, the lattice grid spacing and the time step, c = δx δt represents the particle velocity.
The lattice Boltzmann equation (LBE) reads where f α (x, t) denotes the distribution function at the node of x at the time t, and the BGK collision model, Ω α ( f (x, t)), is given by where λ = λ 0 /δt represents the dimensionless relaxation time, λ 0 is relaxation time.The equilibrium distribution function (EDF), f (eq) α , is given by where u and p are the velocity and the fluid pressure, respectively.The weight coefficients are here c s = c/ √ 3 represents the acoustic speed, ρ 0 is a constant of density.The pressure and velocity of flow are deduced from the ensemble average of distribution functions In the simulation of viscoelastic fluids, the treatment of the elastic force may be crucial for the coupled lattice Boltzmann method.Here, we incorporate the extra forces F = ∇•τ into lattice Boltzmann equation, and treat the elastic force of polymers as an extra forced of the lattice Boltzmann BGK [26].The incompressible LBGK Equation ( 5) with the additional elastic forced term is updated by where the equilibrium distribution f (eq) α is updated by with u(x, t) = u(x, t) the right-hand-side term F α is given by where The flow velocity is finally updated by averaging the value of the after and fore collision The incompressible Navier-Stokes Equations ( 1) and (2) could be derived from lattice Boltzmann Equation (12) using the Chapman-Enskog expansion, with the kinematic viscosity η s = c 2 s (λ − 1/2)δt.Detailed mathematical deduce is given in the Appendix A.

Solving the Oldroyd-B Constitutive Model
In this section, we describe briefly the way to solve the Oldroyd-B Equation (3) using the particle distribution functions of flow fields.To check the information conveniently, we rewrite the constitutive equation as here, the advection operator and relaxation-stretching operator acted on the stress tensor τ are given, respectively, by Here, the physical variables of the Equations ( 19) and ( 20) are distributed on the same grid points as the velocity field.To ensure the numerical precision of simulation, the advection operator is calculated by The relaxation-stretching operators A rel-str (τ, t) are stretched through the velocity gradient tensor.The derivatives of velocity are computed by In the simulation, the time stepping follows as Here, the time step ∆t may normally be different from δt, and the velocity field is fixed during the time when the constitutive equation advance from n to n + 1 time step.Finally, combing with the Equations ( 21) and ( 23), we obtained the temporal advancing of the stress tensor, as follows The detail derivation of the numerical scheme is has been shown in the in reference [21].

The Temporal Marching of Coupled Solvers
The computing steps of the velocity and viscoelastic stress fields are arranged as follows: (I) Equation ( 24) are run firstly to simulate the constitutive equations, of whose results are utilized to calculate the elastic forces; (II) The velocity fields are simulated using Equation ( 14), and the obtained In this part, we mainly discuss the methods of time marching in coupled solvers.For viscoelastic fluid flows at low Re number, the coupled systems have a very huge difference of the time scales between the constitutive modes and lattice Boltzmann evolution.So, we must carefully handle with the time scale when simulating the real geometry device.One approach to time advancement in a coupled framework is to operate the solvers at the smallest scale that required by the lattice Boltzmann equation.However, the computational cost will normally be much larger than it needed.Such a type of method is displayed in Figure 1a and we named it as the fully coupled (FC) method.To improve the efficiency of simulations, the asynchronously coupled method was applied in previous work [21].There are two time steps being used in this coupled framework.The lattice Boltzmann equation runs by the suitable time step corresponding to its physical model, but the macro constitutive equations are solver with fewer numbers of time steps (with larger time-step sizes) than the mesoscale dynamics.The named asynchronously coupled methods mean that the data exchange mode between the lattice Boltzmann equation and macro constitutive equation at every time step.On a separate note, the numerical error that was caused by the scale-separated method may be acceptable if the time step of the macro constitutive is chosen properly.This method is shown schematically in Figure 1b and we refer to it as the asynchronously coupled (AC) method.The coupled algorithm using a homogeneous time step is given, as following Algorithm 1.
Algorithm 1.The coupled algorithm using a homogeneous time step.
For viscoelastic fluids at low Wi number, the asynchronous time marching of the coupled solvers could get steady result in simulations if the macro time step ∆t is chosen properly.However, as Wi number increases, the asynchronously coupled method possesses fairly poor numerical stability.To widen the application of the coupled lattice Boltzmann model, we propose a generalized asynchronously coupled method for time-stably advancing of the viscoelastic lattice Boltzmann method.In this method, the asynchronously coupled way adding an iterative procedure is performed, which may allow for the lattice Boltzmann solver to run with a flexible number of steps (N) before the exchange of coupling variables.The number of added steps, N, gives a further role to relax the coupling instances system.We refer to this method as the generalized asynchronously coupling (GAC) method (see Figure 1c), which can extend the applicability over a wide range of Wi number.The similar ideas have been used in literature [27] where the disparate time scales were applied to study dynamical systems.In brief, the procedure of the generalized asynchronously coupling method is given as following Algorithm 2.
Algorithm 1.The coupled algorithm using a homogeneous time step.

Require:
( , ) For viscoelastic fluids at low Wi number, the asynchronous time marching of the coupled solvers could get steady result in simulations if the macro time step t  is chosen properly.However, as Wi number increases, the asynchronously coupled method possesses fairly poor numerical stability.To widen the application of the coupled lattice Boltzmann model, we propose a generalized asynchronously coupled method for time-stably advancing of the viscoelastic lattice Boltzmann method.In this method, the asynchronously coupled way adding an iterative procedure is It's should be noted at that the GAC scheme we presented in this paper is equivalent to the AC scheme (used in [21]) when N = 1.The flexibility that is provided through the assignment of N enables the best of these methods to be combined in an adaptable manner.In this paper the assignment of N is made as straightforwardly as possible: to minimise the number of macro time steps performed.As such, it is set as large as the restriction on macro time step will allow.This inherent adaptability gives the GAC scheme the range of applicability of the AC method, while maintaining the macro-time-step efficiency as much as possible.

Numerical Results and Discussions
The generalized asynchronously coupled method is first validated by a planar channel flow.Subsequently, to prove reliability of the coupled solvers, the 4:1 contraction problem is also implemented in the simulation.We will give the detail numerical results as follow.

2D Channel Flows
In order to verify the numerical accuracy of coupled solvers, we choose a 2D channel shaped by two parallel panels over distance of about L with H = 10L long (see Figure 2).The characteristic length L and the characteristic velocity U are the height and maximum velocity of at the inlet, respectively.
The velocity boundary conditions paralleled to the x-axis wall are, respectively, The boundary conditions of stress paralleled to the x-axis wall could be obtained by velocity values.The Neumann boundary conditions at the exit are, respectively, used as follows The stress tensor values are initially imposed to zero everywhere, and the fully-developed conditions are used at the entry, which is taken by then the maximum velocity at the inlet is U = 0.1.The Oldroyd-B constitutive equation has an exact solution when the velocity field has reached the steady state, which is written by The boundary conditions of stress paralleled to the x-axis wall could be obtained by velocity values.The Neumann boundary conditions at the exit are, respectively, used as follows The stress tensor values are initially imposed to zero everywhere, and the fully-developed conditions are used at the entry, which is taken by then the maximum velocity at the inlet is U = 0.1.The Oldroyd-B constitutive equation has an exact solution when the velocity field has reached the steady state, which is written by Re ( 1 ) So, we also can study the numerical precision of the presented method by comparing with the exact solution.The relative errors of numerical solutions for stress based on the L2-norm are given by  So, we also can study the numerical precision of the presented method by comparing with the exact solution.The relative errors of numerical solutions for stress based on the L 2 -norm are given by where M is the total number of nodes, τ numerical xy , τ numerical xx denote the steady-state numerical solutions and τ exact xy , τ exact xx are the exact solutions.Generally speaking, at low Wi number, asynchronously coupled method can get a stability result in simulations only if the δt and ∆t can guarantee the adequate resolutions of each solvers because of the small scale-separation.The value of N is related with the ratio of time-step sizes (∆t/δt), but is usually not linear with respect to it.In practical applications, N is found by the trial and error approach.For instance, we initially set N 0 = 2 h , h = lg(∆t/δt) , present round to ceiling.We can take N = N 0 , if the iterative error of was decreased when N = N 0 , other else we use N 0 = 2 h+1 to test once again, and so on.We first run simulations using δt ≈ 1 × 10 −5 and ∆t ≈ 1 × 10 −2 with the fined mesh of 400 × 40, and N = 5 for generalized asynchronously coupled method, at lower Wi = 1, and get the results in the Figure 3, comparing AC scheme with GAC scheme.We find that despite both AC scheme and GAC scheme can get a stability result in simulations, the convergence speed of GAC scheme is fast than the AC scheme.
Appl.Sci.2018, 8, x FOR PEER REVIEW 10 of 17 where M is the total number of nodes, numerical Generally speaking, at low Wi number, asynchronously coupled method can get a stability result in simulations only if the t  and t  can guarantee the adequate resolutions of each solvers because of the small scale-separation.The value of N is related with the ratio of time-step sizes , but is usually not linear with respect to it.In practical applications, N is found by the trial and error approach.For instance, we initially set ceiling.We can take N = N0, if the iterative error of was decreased when N = N0, other else we use N0 = 2 h+1 to test once again, and so on.We first run simulations using with the fined mesh of 400 × 40, and N = 5 for generalized asynchronously coupled method, at lower Wi = 1, and get the results in the Figure 3, comparing AC scheme with GAC scheme.We find that despite both AC scheme and GAC scheme can get a stability result in simulations, the convergence speed of GAC scheme is fast than the AC scheme.To show the validity of the generalized asynchronously coupled method, we then take the dimensionless units U as 0.1, L as 1, Wi = 15.0,parameter  as 0.5, and  = 1, Re number as 1.0, then Wi 10   , and compare the numerical results of our method with the analytical solution and asynchronously coupled scheme.We also set for coupled using the fined mesh of 400 × 40, and N = 5 for generalized asynchronously coupled method.Figure 4 shows that the relative error of x y  and x x  simulated by our present scheme is extremely small at Wi = 15, whereas the relative error of obtained by asynchronously scheme is clearly exceeded more than 100%.In other words, the asynchronously coupled method is considerably subject to numerical error at Wi = 15.Nevertheless, when comparing with the standard asynchronously coupled methods, the generalized method persists to give stable solutions up to Wi = 15.
Overall, seeing from the analysis of relative error, we found that the present method could be immune to problem of Wi = 15.Moreover, we plot the velocity components u, as obtained by numerical method at To show the validity of the generalized asynchronously coupled method, we then take the dimensionless units U as 0.1, L as 1, Wi = 15.0,parameter β as 0.5, and ρ = 1, Re number as 1.0, then λ p = LWi/U = 10, η t = ρUL/Re = 1, and compare the numerical results of our method with the analytical solution and asynchronously coupled scheme.We also set δt ≈ 1 × 10 −5 and ∆t ≈ 1 × 10 −2 for coupled using the fined mesh of 400 × 40, and N = 5 for generalized asynchronously coupled method.Figure 4 shows that the relative error of τ xy and τ xx simulated by our present scheme is extremely small at Wi = 15, whereas the relative error of obtained by asynchronously scheme is clearly exceeded more than 100%.In other words, the asynchronously coupled method is considerably subject to numerical error at Wi = 15.Nevertheless, when comparing with the standard asynchronously coupled methods, the generalized method persists to give stable solutions up to Wi = 15.Overall, seeing from the analysis of relative error, we found that the present method could be immune to problem of Wi = 15.Moreover, we plot the velocity components u, as obtained by numerical method at centre point (5, 0.5) on time t in Figure 5, comparing with the results of the asynchronously coupled method at β = 0.5 and Re = 1.0.From Figure 5, we found that, before eventually reaching a steady state, the velocity profiles overshoot the terminal velocity due to the action of the elasticity.However, the result of the asynchronously coupled method is divergent as the time increasing because of its poor stability.centre point (5, 0.5) on time t in Figure 5, comparing with the results of the asynchronously coupled method at 0.5   and Re = 1.0.From Figure 5, we found that, before eventually reaching a steady state, the velocity profiles overshoot the terminal velocity due to the action of the elasticity.However, the result of the asynchronously coupled method is divergent as the time increasing because of its poor stability.  and Re = 1.0, respectively.

Contraction Flows
Viscoelastic fluids through contractions are important in the rapidly developed research field of micro fluidics, such as the polymer processing.As we know, the contraction flow has singularity at the reentrant corner, which makes calculations at high Weissenberg numbers challenging.The planar 4:1 contraction was seen as a benchmark [28] model in 1987 because of the geometrical simplicity and the known numerical difficulty.In this section, to further test the validity of the method presented in this paper, we simulate the Oldroyd-B fluids through a 4:1contraction geometry at Wi = 5.
Figure 6 shows the flow domain of which the ratio between the downstream and upstream height is 1 to 4. The characteristic length L and the characteristic velocity U are the half height and average velocity at the outlet, respectively.A parabolic Poiseuille flow is set at inflow by 2 3 ( 16) 128 while the fully-developed conditions are given at the outlet, so the average velocity U at the out let is 1.No-slip boundary conditions are given along the stationary walls.Symmetry boundary  centre point (5, 0.5) on time t in Figure 5, comparing with the results of the asynchronously coupled method at 0.5   and Re = 1.0.From Figure 5, we found that, before eventually reaching a steady state, the velocity profiles overshoot the terminal velocity due to the action of the elasticity.However, the result of the asynchronously coupled method is divergent as the time increasing because of its poor stability.

Contraction Flows
Viscoelastic fluids through contractions are important in the rapidly developed research field of micro fluidics, such as the polymer processing.As we know, the contraction flow has singularity at the reentrant corner, which makes calculations at high Weissenberg numbers challenging.The planar 4:1 contraction was seen as a benchmark [28] model in 1987 because of the geometrical simplicity and the known numerical difficulty.In this section, to further test the validity of the method presented in this paper, we simulate the Oldroyd-B fluids through a 4:1contraction geometry at Wi = 5.
Figure 6 shows the flow domain of which the ratio between the downstream and upstream height is 1 to 4. The characteristic length L and the characteristic velocity U are the half height and average velocity at the outlet, respectively.A parabolic Poiseuille flow is set at inflow by 2 3 ( 16) 128 while the fully-developed conditions are given at the outlet, so the average velocity U at the out let is 1.No-slip boundary conditions are given along the stationary walls.Symmetry boundary

Contraction Flows
Viscoelastic fluids through contractions are important in the rapidly developed research field of micro fluidics, such as the polymer processing.As we know, the contraction flow has singularity at the reentrant corner, which makes calculations at high Weissenberg numbers challenging.The planar 4:1 contraction was seen as a benchmark [28] model in 1987 because of the geometrical simplicity and the known numerical difficulty.In this section, to further test the validity of the method presented in this paper, we simulate the Oldroyd-B fluids through a 4:1contraction geometry at Wi = 5.
Figure 6 shows the flow domain of which the ratio between the downstream and upstream height is 1 to 4. The characteristic length L and the characteristic velocity U are the half height and average velocity at the outlet, respectively.A parabolic Poiseuille flow is set at inflow by while the fully-developed conditions are given at the outlet, so the average velocity U at the out let is 1.
No-slip boundary conditions are given along the stationary walls.Symmetry boundary conditions are imposed along the center-line.We take the dimensionless units U as 1, L as 1, parameter β as 1/9, and Re number as 1.0.conditions are imposed along the center-line.We take the dimensionless units U as 1, L as 1, parameter  as 1 9 , and Re number as 1.0.To verify the mesh convergence, we carried out the calculations of a steady flow with 1.0  using three Cartesian meshes, M1, M2 and M3.The total number of cells and the grid size of each mesh are reported in Table 1. Figure 7 displays the temporal evolution of the local velocity at the intersection of the symmetry line and the contraction plane of the channel, for the different mesh.From Figure 7, we can see that the velocity tends to be meshing independent as the number of the nodes increased.To juggle the computational precision and efficiency, the numerical results are got in the next parts by utilizing the dense mesh of M3.For the generalized asynchronously coupled method, we set    To verify the mesh convergence, we carried out the calculations of a steady flow with Wi = 1.0 using three Cartesian meshes, M1, M2 and M3.The total number of cells and the grid size of each mesh are reported in Table 1. Figure 7 displays the temporal evolution of the local velocity at the intersection of the symmetry line and the contraction plane of the channel, for the different mesh.From Figure 7, we can see that the velocity tends to be meshing independent as the number of the nodes increased.To juggle the computational precision and efficiency, the numerical results are got in the next parts by utilizing the dense mesh of M3.For the generalized asynchronously coupled method, we set δt ≈ 1.6 × 10 −5 , ∆t ≈ 1 × 10 −3 and N = 4. Figure 8 displays the numerical results of velocity components along the horizontal center lines of the contraction with Wi = 1.0, comparing between the results of GAC and AC method, respectively.We can see that the difference of those two schemes is not obvious at low Wi number flows.That study of low Wi number flows is a benchmark example when considering it as a preview to the high Wi number flows in the following part.To verify the mesh convergence, we carried out the calculations of a steady flow with 1.0 Wi  using three Cartesian meshes, M1, M2 and M3.The total number of cells and the grid size of each mesh are reported in Table 1. Figure 7 displays the temporal evolution of the local velocity at the intersection of the symmetry line and the contraction plane of the channel, for the different mesh.From Figure 7, we can see that the velocity tends to be meshing independent as the number of the nodes increased.To juggle the computational precision and efficiency, the numerical results are got in the next parts by utilizing the dense mesh of M3.For the generalized asynchronously coupled method, we set        To obtain some complex viscoelastic behaviors, we consider the dynamics of viscoelastic fluids at Wi = 5.0 using the GAC in this part.Figure 9 displays a set of instantaneous streamlines of flows with Wi = 5.0 on the finest mesh M3.Concerning the dynamics of viscoelastic fluids, we can find a series of different flow regimes, while increasing the Wi number to 5.0.From Figure 9, we observe that the elastic lip vortex grows in size as time goes by, eventually running up to the corner vortex region, and mixing with each other in a fairly complex dynamic process.The dynamical transitions are signed by a complex interaction between pulsating lip and corner vortices, which tend to approach and separate.In the literatures, through a full 4:1 planar contraction, Comminal et al. [5] and Afonso et al. [29], by carrying out numerical simulations of the Oldroyd-B model, also predicted the similar coalescence of lip and salient corner vortices into one big vortex with increasing elasticity.When comparing quantitatively with the previous works, we find that the prediction of our presented method reproduces the numerical results reported in the previous studies.However, under the same condition, the procedure using the AC scheme does not obtain the converged solution as the Wi number increases to 5.0.Overall, thorough the stringent test case, we have found that the generalized coupled lattice Boltzmann method could be fit to simulate the Oldroyd-B fluids when the Wi number is increased.To obtain some complex viscoelastic behaviors, we consider the dynamics of viscoelastic fluids at Wi = 5.0 using the GAC in this part.Figure 9 displays a set of instantaneous streamlines of flows with Wi = 5.0 on the finest mesh M3.Concerning the dynamics of viscoelastic fluids, we can find a series of different flow regimes, while increasing the Wi number to 5.0.From Figure 9, we observe that the elastic lip vortex grows in size as time goes by, eventually running up to the corner vortex region, and mixing with each other in a fairly complex dynamic process.The dynamical transitions are signed by a complex interaction between pulsating lip and corner vortices, which tend to approach and separate.In the literatures, through a full 4:1 planar contraction, Comminal et al. [5] and Afonso et al. [29], by carrying out numerical simulations of the Oldroyd-B model, also predicted the similar coalescence of lip and salient corner vortices into one big vortex with increasing elasticity.When comparing quantitatively with the previous works, we find that the prediction of our presented method reproduces the numerical results reported in the previous studies.However, under the same condition, the procedure using the AC scheme does not obtain the converged solution as the Wi number increases to 5.0.Overall, thorough the stringent test case, we have found that the generalized coupled lattice Boltzmann method could be fit to simulate the Oldroyd-B fluids when the Wi number is increased.

Conclusions
In conclusions, we have employed a coupled lattice Boltzmann method, with a generalized asynchronously coupling scheme, to simulate the Oldroyd-B fluids when the Weissenberg number is increased.This coupled method includes two types of equations, one is Lattice Boltzmann

Conclusions
In conclusions, we have employed a coupled lattice Boltzmann method, with a generalized asynchronously coupling scheme, to simulate the Oldroyd-B fluids when the Weissenberg number is increased.This coupled method includes two types of equations, one is Lattice Boltzmann equation of the solvent and the other is constitutive equation of the stress tensor.The two equations are solved using two different time steps-one is of macroscopic scale and the other is of mesoscopic scale.For viscoelastic fluids at a higher Wi number, the computational stability of asynchronously coupling is greatly affected by the polymer stress as the elasticity interaction becomes relatively large.Therefore, we have presented a generalized technique for time-stably advancing of the coupled lattice Boltzmann method, which can improve the numerical stability of the simulation at a higher Wi number.The main idea of the new method is that one can add an adaptable number of micro solver time steps before the transmission of data between the coupling variables.
We have shown that the numerical results, as calculated by the new method, are in good agreement with analytical results of the two-dimensional channel flow.In addition, for viscoelastic fluids through contractions at Wi = 5, we have found that the prediction of our presented method can reproduce the same numerical results that were obtained by previous literatures.The main advantage of our presented method is that it could be use to simulate the Oldroyd-B fluids when the Wi number is increased.Furthermore, our new method could also be easily implemented, which keeps the traditional advantageous of LBM.Hence, the coupled lattice Boltzmann method that was described in this paper could be seen as an alternative numerical technique for simulating the viscoelastic fluid.
to calculate the related velocity gradient; (III) The velocity and elastic stress fields are updated iteratively at each step.

Figure 2 .
Figure 2. Channel flow set up parameters.

Figure 2 .
Figure 2. Channel flow set up parameters.

Figure 3 .
Figure 3.Comparison of the relative L2-error of x y  and x x  for y at position x = 5 at Wi = 1.0 with

Figure 3 .
Figure 3.Comparison of the relative L 2 -error of τ xy and τ xx for y at position x = 5 at Wi = 1.0 with β = 0.5 and Re = 1.0, respectively.

Figure 4 .
Figure 4. Comparison of the relative L2-error of x y  and x x  for y at position x = 5 at Wi = 15.0

Figure 5 .
Figure 5.The time development of the velocity components u central point (5, 0.5) for Wi = 15.0 (right) with .5

Figure 4 .
Figure 4. Comparison of the relative L 2 -error of τ xy and τ xx for y at position x = 5 at Wi = 15.0 with β = 0.5 and Re = 1.0, respectively.

Figure 4 .
Figure 4. Comparison of the relative L2-error of x y  and x x  for y at position x = 5 at Wi = 15.0

Figure 5 .
Figure 5.The time development of the velocity components u at central point (5, 0.5) for Wi = 15.0 (right) with 0 .5   and Re = 1.0, respectively.
and N = 4. Figure8displays the numerical results of velocity components along the horizontal center lines of the contraction with Wi = 1.0, comparing between the results of GAC and AC method, respectively.We can see that the difference of those two schemes is not obvious at low Wi number flows.That study of low Wi number flows is a benchmark example when considering it as a preview to the high Wi number flows in the following part.

Figure 7 .
Figure 7. Mesh dependence of the velocity components through the horizontal at Wi = 1.0 using mesh of M1, M2, and M3.

Figure 8 .
Figure 8.The numerical solutions of velocity components through the horizontal centerlines with Wi = 1.0, comparing between the results of GAC and the AC, respectively.
Appl.Sci.2018, 8, x FOR PEER REVIEW 12 of 17conditions are imposed along the center-line.We take the dimensionless units U as 1, L as 1, parameter  as 1 9 , and Re number as 1.0.
and N = 4. Figure8displays the numerical results of velocity components along the horizontal center lines of the contraction with Wi = 1.0, comparing between the results of GAC and AC method, respectively.We can see that the difference of those two schemes is not obvious at low Wi number flows.That study of low Wi number flows is a benchmark example when considering it as a preview to the high Wi number flows in the following part.

Figure 7 .
Figure 7. Mesh dependence of the velocity components through the horizontal at Wi = 1.0 using mesh of M1, M2, and M3.

Figure 8 .
Figure 8.The numerical solutions of velocity components through the horizontal centerlines with Wi = 1.0, comparing between the results of GAC and the AC, respectively.

Figure 7 .
Figure 7. Mesh dependence of the velocity components through the horizontal at Wi = 1.0 using mesh of M1, M2, and M3.

Figure 7 .
Figure 7. Mesh dependence of the velocity components through the horizontal at Wi = 1.0 using mesh of M1, M2, and M3.

Figure 8 .
Figure 8.The numerical solutions of velocity components through the horizontal centerlines with Wi = 1.0, comparing between the results of GAC and the AC, respectively.

Figure 8 .
Figure 8.The numerical solutions of velocity components through the horizontal centerlines with Wi = 1.0, comparing between the results of GAC and the AC, respectively.

Table 1 .
Mesh characteristics used in this study.

Table 1 .
Mesh characteristics used in this study.
Mesh Total Cell Δ x (x/L)