Complex fluids comprise a large class of soft materials, such as polymeric solutions, liquid crystals, ionic solutions, and fiber suspensions. These are fluids with complicated rheological phenomena, arising from different “elastic” effects, such as elasticity of deformable particles, interaction between charged ions and bulk elasticity endowed by polymer molecules [liu2009introduction]. Modeling and numerical simulations of complex fluids have been interesting problems for a couple of decades [bird1987, larson1998, li2007mathematical, le2012micro].
Models for complex fluids can be classified as pure macroscopic models[Keunings1997, Lielens1998, Sizaire1999] and micro-macro models [bird1992transport, le2012micro]
. The macroscopic models utilize a closed-form constitutive equation for the stress tensorto supplement the conservation of mass and momentum [Keunings1997, Lielens1998, Sizaire1999]. The advantage of this approach is low computational cost, however, it is often difficult to construct a closed-form of the constitutive equation that captures the flow behaviors of complex fluids. The idea of micro–macro approaches is to couple the macroscopic conservation laws with a microscopic kinetic theory, which describes the origin of the macroscopic stress tensor [bird1992transport, le2012micro]. One of the simplest micro-macro models of complex fluids is the dumbbell model for dilute polymeric fluids given by [lin-liu-zhang]
In this model, the macroscopic motion of the fluid is described by a Navier–Stokes equation with an elastic stress depending on the microscopic configuration of polymer chains, where denotes the fluid velocity field, is the pressure, is the density of the fluid, and is a parameter related to the polymer density. On the microscopic level, a polymer chain is modeled by a elastic dumbbell with two beads connected by a single spring. The molecular configuration is described by an end-to-end vector of the dumbbell, see Fig. 1, and is the spring potential. The interactions among polymer chains are neglected due to the dilute assumption. In the case that ( is the elastic constant), the model is known as the Hookean dumbbell model. While the case that ( is the maximum dumbbell extension) is known as the FENE (Finite Extensible Nonlinear Elastic) dumbbell model. The microscopic dynamics is described by a Fokker–Planck equation of the number density distribution function with a drift term depending on the macroscopic velocity field , where is a constant related to the polymer relaxation time, is the Boltzmann constant and is the absolute temperature. Alternatively, the microscopic dynamics can also be described by a stochastic differential equation (SDE), or Langevin dynamics, given by [Bris2012]
is the standard multidimensional white noise.
Although micro-macro models give an elegant description to the origin of the macroscopic stress tensor, it is a long lasting challenge to simulate micro-macro models directly. During the past decades, various computational techniques have been developed for various micro-macro models.[weinan2009general, Halin1998, Hulsen1997, Laso1993, ottinger1996, Lozinski2011] There are mainly two approaches, the Langevin equation based stochastic simulation methods and direct simulation methods based on the Fokker–Planck equation [lozinski2011langevin]. The CONNFFESSIT (Calculation of Non-Newtonian Flow: Finite Elements and Stochastic Simulation Technique) algorithm , was the first Langevin equation based numerical method, which couples a finite element discretization to the macroscopic flow with a numerical solver for the microscopic SDE (2) [Laso1993, ottinger1996]. Along this direction, other stochastic approaches, the Lagrangian particle method (LPM) [Halin1998] and the Brownian configuration field (BCF) method [Hulsen1997]
were proposed in order to reduce the variance and the computational cost of the original CONNFFESSIT algorithm. Several extensions of these approaches and the corresponding numerical experiments have been extensively investigated in recent years[chauviere2002, Koppol2007, Jourdain2004, Boyaval2010, Bris2012, Griebel2014, XuSPH2014]. Stochastic approaches have been the dominant simulation methods of the micro-macro model. However, such stochastic approaches suffer from several shortcomings, including the high computational cost and the stochastic fluctuations. An alternative approach is the direct simulation of the equivalent Fokker–Planck equation in the configurational space. Examples include Galerkin spectral element technique[suen2002, Lozinski2003jcp, Chauviere2004, KnezevicFP, shen2012approximation] and the lattice Boltzmann technique[Ammar2010, Bergamasco2013]. As pointed out in Ref. Lozinski2011, the computational cost of Fokker-Planck-based methods increases rapidly for simulations in strong flows (with highly localized distribution function) or involving high dimensional configuration spaces. Besides, there are a huge literature devoted to developing macroscopic closure systems, as an approximation to the original macro-micro model [feng1998closure, du2005fene, hyon2008maximum]. Although the closure systems are easier to solve, they may not able to capture some flow behaviors, such as hysteresis, of the original micro-macro model. [Hyon2008].
During the past decades, there have been growing interests in solving the Fokker-Planck equation by a deterministic particle method [russo1990deterministic, degond1990deterministic, lacombe1999presentation, Carrilloblob, wangEVI]. Compared with the Langevin dynamics based stochastic particle methods, the deterministic particle methods are often computationally cheap and do not suffer from stochastic fluctuations. The main challenge of applying the deterministic particle method to the Fokker-Planck type equation is how to deal with the diffusion term. Various kernel regularization approaches have been developed to overcome this difficulty. [russo1990deterministic, degond1990deterministic, lacombe1999presentation, Carrilloblob]
The goal of this paper is to develop an efficient numerical method for the micro-macro model (1) by coupling a finite element discretization to the macroscopic flow with a deterministic particle discretization to the microscopic Fokker-Planck equation. As pointed out in Ref. wangEVI, applying the kernel regularization to the equation may not preserve variational structure at the particle level. To preserve the variational structure, we apply the particle approximation at the energy-dissipation law level and employ a discrete energetic variational approach [wang2020jcp, wangEVI] to derive a coarse-grained micro-macro model with a particle approximation first. A structure-preserving particle-FEM discretization is developed for the coarse-grained model. Various numerical experiments have been performed to validate the new scheme via several benchmark problems. Despite its simplicity, the deterministic particle method is robust, accurate and able to catch certain complex behaviors of polymeric fluids. The numerical results obtained by our scheme are in excellent agreement with those from the former work. Moreover, the deterministic particle discretization is shown to be more efficient than the stochastic particle approaches, where a large ensemble of realizations of the stochastic process are needed.
The rest of this paper is organized as follows. A formal derivation of the micro-macro model of dilute polymeric fluids by employing the energetic variational approach in given in section 2. Then we derive a coarse-grained micro-macro model with the particle approximation and construct a deterministic particle-FEM scheme in Section 3. Numerical simulations to validate our scheme are presented in Section 4.
2 Energetic variational approach to the micro-macro model
In this section, we give a formal derivation to micro-macro models of dilute polymeric fluids by employing the energetic variational approach (EnVarA).
Starting with a prescribed energy-dissipation law and a kinematic (transport) relation, the framework of EnVarA derives the dynamics of a non-equilibrium system through two distinct variational processes: the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP).[MDP2, onsager1931reciprocal2] This approach was developed from pioneering works of Rayleigh [Rayleigh1871] and Onsager [onsager1931reciprocal, onsager1931reciprocal2], and has been successfully applied to build up many mathematical models in physics, chemical engineering and biology.[wang2020field, Giga2017, liu2009introduction, wang2021variational] From a numerical perspective, the EnVarA formulation also provides a guide line to develop structure-preserving numerical schemes for different systems.[wang2020jcp, liu2020variational, liu2021structure]
For an isothermal and closed system, an energy-dissipation law is often given by
which comes from the first and second laws of thermodynamics [ericksen1992introduction, Giga2017]. Here is the total energy, which is the sum of the Helmholtz free energy and the kinetic energy ; stands for the rate of energy dissipation, which equals to the entropy production in this case. The LAP states that, the dynamics of the Hamiltonian part of the system can be determined by taking variation of the action functional with respect to (the trajectory in Lagrangian coordinates) [arnol2013mathematical, Giga2017], which derives the conservative force
where is the physical domain at time . For a dissipative system (), the dissipative force can be obtained by minimizing the Onsager dissipation functional with respect to the “rate” in the linear response regime, i.e.,
This principle is known as Onsager’s MDP [onsager1931reciprocal, onsager1931reciprocal2, Giga2017]. Consequently, the force balance condition (Newton’s second law, in which the inertial force plays a role of ) results in
which is the dynamics of the system.
Now we show how to derive the micro-macro model (1) formally under the EnVarA framework. As mentioned in the introduction, the polymer molecules are modeled as an elastic dumbbell consisting of two ”beads” joined by a one-dimensional spring at the microscopic scale. [li2007mathematical, lin-liu-zhang] The microscopic configuration between the two beads is described by an end-to-end vector . Let be the number density distribution function of finding a molecule with end-to-end vector at position at time . To derive the micro-macro model by EnVarA, one needs to introduce Lagrangian descriptions in both micro- and macro- scales. Let be the flow map at the physical space and be the flow map at the configurational space, where and are Lagrangian coordinates in physical and configurational space respectively. For given flow maps and , the corresponding macroscopic velocity and the microscopic velocity satisfy
Moreover, one can define the deformation tensor associated with the flow map by
Without ambiguity, in this paper, we will not distinguish and . Obviously, carries all the transport information of configurations in the system [linfh2012] and satisfies the transport equation in Eulerian coordinates [lin-liu-zhang]
Due to the conservation of mass, the density distribution function satisfies
which can be written as
in Eulerian coordinates. Here and are the macro- and microscopic velocities defined in (7).
In the framework of EnVarA, the micro-macro system can be modeled through an energy-dissipation law
where is the constant density of the fluid, is a constant that represents the polymer density, is the Boltzmann constant, is the absolute temperature, is the solvent viscosity, the constant is related to the polymer relaxation time, is the microscopic elastic potential of the polymer molecules. For Hookean and FENE models, the elastic potential is given by
respectively, where is the elastic constant, is the maximum dumbbell extension in FENE models. The second term of the dissipation accounts for the micro-macro coupling with being the macroscopic induced velocity. According to the Cauchy-Born rule, at the macroscopic scale, which indicates
Now we are ready to derive the dynamics of the system. First, we look at the dynamics at the macroscopic scale. Due to the ”separation of scale” [Giga2017], the second term on the right hand side of the dissipation (11) vanishes when deriving the macroscopic force balance. Since , the action functional can be written as
in Lagrangian coordinates, where is the initial number distribution function, and due to . By applying the LAP, i.e., taking variation of with respect to , we get
in Eulerian coordinates. Indeed, consider a perturbation
where satisfying . Then
Push forward to Eulerian coordinates, we have
which leads to (13).
For the dissipation part, the MDP, i.e., taking variation of with respect to , leads to
where is the Lagrangian multiplier for the incompressible condition . Hence, the macroscopic force balance results in the momentum equation
is the induced stress from the configuration space, representing the microscopic contributions to the macroscopic level. Here, denotes a tensor product and is a matrix for two vectors and .
On the microscopic scale, by taking variations with respect to and , we obtain
Then combining with Eq. (10), we get the equation on the microscopic scale:
And thus, the final coupled system reads as follows,
subject to a suitable boundary condition. For the sake of simplicity, in this paper, we consider Dirichlet boundary conditions for the velocity and non-flux boundary condition for with respect to .
3 Numerical Methods
In this section, we propose a numerical scheme to solve the micro-macro model by combining a finite element discretization for the macroscopic fluid dynamic equation [becker2008, chenrui2015, bao2021] with a deterministic particle method for the microscopic Fokker-Planck equation [wangEVI].
To preserve the variational structure after a particle discretization, we first derive a coarse-grained micro-macro system with a deterministic particle approximation in the configurational space by a discrete energetic variational approach [wang2020jcp]. The discrete energetic variational approach follows the idea of ”Approximation-then-Variation”, which first applies particle approximation to the continuous energy-dissipation law and derives a coarse-grained model by variational procedures. As an advantage, the derived equations of particles preserve the variational structure at the particle level [wangEVI].
For simplicity, we assume that
which means that the number density of polymer chains is spatially homogeneous. Thus, for fixed , can be approximated by
with a particle approximation at the microscopic level. Here is the number of particles at and time , is a set of representative particles at at time , is the weight of the corresponding particle satisfying . In the current work, we fix , i.e., all the particles are equally weighted.
can be viewed as representative particles that represent information of the number density distribution at . Since only needs to be computed at each time-step, the computational cost can be largely reduced. It is more like an Eulerian approach, rather than a Lagrangian approach.
Substitute the approximation (21) into the continuous energy-dissipation law (11), we can obtain a discrete energy-dissipation law in terms of and the macroscopic flow. However, from a computational point of view, the term can not be defined in a proper way. To deal with this issue, we introduce a kernel regularization, i.e., replacing by , where is a kernel function and
A typical choice of is the Gaussian kernel, given by
Here is the bandwidth which controls the inter-particle distances and is the dimension of the space. We take as a constant for simplicity. The values of will affect the numerical results. We’ll discuss the choices of in the next section.
Within the kernel regularization, the discrete energy can be written as
and the discrete dissipation is
where is the material derivative of .
The differential equations of can be derived via the discrete energetic variational approach [wang2020jcp], i.e., performing the EnVarA in terms of and , which leads to
By direct computation we get the equation for
Here we denote by for convenience. As an advantage of the “approximation-then-variation” approach, it can be noticed that Eq. (24) is a gradient flow with respect to in absence of the flow, i.e. at .
The variational procedure for the macroscopic flow is almost the same as that in the continuous case, shown in section 2. The final micro-macro system with particle approximation is given by
where satisfies (24).
One can view the macroscopic flow equation (25) along with the microscopic evolution equation (24) as a coarse-grained model for the original micro-macro model (19). To solve the system numerically, it is a natural idea to apply some decoupled schemes at each step. Precisely, we apply the following algorithm for the temporal discretization:
Step 1: Solve the equation by a finite element method to obtain updated values for the velocity and pressure by treating the viscoelastic stress field as a known term obtained from the last time step.
Step 2: Using the updated velocity field , solve the equation of .Update the values of the viscoelastic stress at each node, and project them into the finite element space of .
In the present work, we adopt the finite element method developed in Ref. becker2008 and Ref. chenrui2015 for Step 1. For more detail, let be the bounded computational domain, and be the triangulation of . consists of a set of simplexes and a set of nodal points . For a nonnegative integer , denote as the space of polynomial functions of degree less than or equal to on the simplex . We can construct finite-dimensional subspaces , and as follows,
Let the finite element spaces for the velocity and the stress tensor be and respectively, where is the dimension of space. Then the equation (25) can be solved by a standard velocity-correction projection method [Guermond2006]. The advantage of using the element for the velocity and pressure is that it avoids an inf-sup condition [mixedFEM].
Next we discuss how to solve equation (24) with a given flow field . The numerical difficulties arise from the fact that is a function of and due to the convection term . Many earlier numerical studies based on CONNFFESSIT algorithms either focus on the shear flows in which the convection term vanishes or ignore the convection term. [ottinger1996, Jourdainanalysis1] To deal with the convection term in stochastic methods, there are two types of methods have been developed. One is to introduce a spatial-temporal discretization to , as used in Brownian configuration field method [ottinger1997brownian, du2005fene, XuSPH2014]. Another way is to use a Lagrangian viewpoint to compute the convection term [Halin1998]. In the current study, we use the idea of the second approach, and use an operator splitting approach to solve (24). Initially, we assign ensemble of particles to each node (). We assume that is spatially homogeneous, and use the same ensemble of initial samples at all . Within the values , we solve the microscopic equation (24) by the following two steps:
Step 1: At each node , solve (24) without the convection term by
Step 2: To deal with the convention term, we view each node as a Lagrangian particle, and update it according to the Eulerian velocity field at each node
Hence, is an ensemble of samples at the new point . To obtain at
, we use some interpolation to get(at mesh with being the set of node) from (at mesh with being node) for each .
An advantage of the above update-and-projection approach is that it doesn’t require a spatial discretization on .
The operator splitting approach has been widely used in many previous Fokker-Planck based numerical approaches [helzel2006multiscale]. One important reason is that the system admits a variational structure without the convention terms. In addition, by treating the convention part separately, we treat the particles at each physical location independently, which largely saves the computational cost.
Since the first step in (26) admits a variational structure, the implicit Euler discretization can be reformulated as an optimization problem. In more detail, we define
In this case, the first step of (26) is the gradient of with respect to . Hence, we can solve the nonlinear system by solving the optimization problem
An advantage of this reformulation is that we can prove the existence of the . More precisely, we have the following result.
For any given , there exists at least one minimal solution of (28) that also satisfies
Moreover, when the velocity , then converges to a stationary point of as . In the case of , is reduced to
Let () be vectorized , namely, Denote and as and respectively. For given , we define
be the admissible set. Obviously, is non-empty and closed, since and is continuous. Moreover, it’s easy to prove that is bounded from below, since
And thus, is coercive and is a bounded set. Hence, admits a global minimizer in . Since is a global minimizer of , we have
which is equivalent to Eq. (29). For the series , we obtain
Summing from 0 to , we have
where is some constant independent of . Hence,
indicating the convergence of . Moreover, due to Eq. (26):
And thus, converges to a stationary point of as .
The full discretization scheme can be summarized as follows. Given the time step size , the initial conditions , , and , having computed for , , and for , we compute , , and by the following algorithm:
Step 1: Treat the stress tensor explicitly and solve the macroscopic flow equation
The equation (33) can be solved by the following two steps:
Step 1.1: Find , such that for any ,
Step 1.2: Find , such that for any ,
Then we obtain
Step 2: Within the values , we solve the microscopic equation (24) by the update-and-project approach. Within the ensemble of particles on each node , the updated values of the viscoelastic stress at each node, denoted as , can be obtained through the second equation of Eq. (25). And then project them into the finite element space of , i.e. . To this end, we choose the projection operator , such that, for each component of the stress with , , where denotes the nodal basis for .
In the previous numerical scheme, we estimate the macroscopic stress tensor by taking microscopic distribution function as the empirical measure for the finite number of particles. More advanced techniques can be applied to this stage to obtain a more accurate estimation to the stress tensor, such as the maximum-entropy based algorithm developed in Ref. arroyo2006local and Ref. rosolen2013adaptive that reconstructs basis functions from particles. We’ll explore this perspective in future works.
4 Numerical Experiments
In this section, we perform various numerical experiments to validate the proposed numerical scheme by studying various well-known benchmark problems for the micro-macro models [Laso1993, ottinger1996, Hyon2014]. For all the numerical experiments carried out in this section, we suppose that the flow is two-dimensional and the dumbbells lie in the plane of the flow, namely, the configuration vector is also two-dimensional. At each node, we use the same initial ensemble of
particles, sampled from the 2-dimensional standard normal distribution.
4.1 Hookean case: the accuracy test
We first consider the Hookean model in a simple shear flow [Bris2012, Laso1993], where the fluid is enclosed between two parallel planes of infinite length separated by a distance . At , the lower plane starts to move in the positive direction with a constant velocity , as shown in Fig. 2. We consider the no-slip boundary conditions at the walls. In this case, the velocity is in the -direction, and only depends on the -variable, namely, the velocity with . Obviously, the velocity field automatically satisfies the incompressible condition . Moreover, we assume that also only depends on , so .
Since the micro-macro model (19) with the Hookean potential is equivalent to a macroscopic viscoelastic model, the Oldroyd–B model. We can compare the simulation results for the micro-macro model and the corresponding Oldroyd–B model, as an accuracy test for the proposed numerical scheme. In this subsection, we take the parameters as follows,
Such parameters might not be physically reasonable, and we use them primarily to illustrate our numerical scheme. When the parameters are set as Eq. (34), the corresponding Oldroyd-B model can be written as follows [Jourdainanalysis1],
Here, denotes the component of the stress . In this work, these equations are solved sequentially in a segregated manner [Favero2010, Pimenta2017]. Namely, for , having computed and , we compute and as follows,
For the spatial discretization, we use the same finite element method mentioned in Section 3 to simulate the Oldroyd-B model. For both models, we set the number of elements and . The kernel bandwidth [liu2016stein], where med is the median of the pairwise distance between the particles .