A two species micro-macro model of wormlike micellar solutions and its maximum entropy closure approximations: An energetic variational approach
aa r X i v : . [ phy s i c s . f l u - dyn ] J a n A two species micro-macro model of wormlike micellar solutions and itsmaximum entropy closure approximations: An energetic variationalapproach
Yiwei Wang a , Teng-Fei Zhang b , Chun Liu a a Department of Applied Mathematics, Illinois Institute of Technology, Chicago, IL 60616, USA b School of Mathematics and Physics, China University of Geosciences, Wuhan,430074, China
Abstract
Wormlike micelles are self-assemblies of polymer chains that can break and reform reversibly. In thispaper, we derive a thermodynamically consistent two species micro-macro model of wormlike micellarsolutions by employing an energetic variational approach. The model incorporates a breaking and re-forming process of polymer chains into the classical micro-macro dumbbell model of polymeric fluids in aunified variational framework. We also study different maximum entropy closure approximations to thenew model by “variation-then-closure” and “closure-then-variation” approaches. By imposing a properdissipation in the coarse-grained level, the closure model, obtained by “closure-then-variation”, preservesthe thermodynamical structure of both mechanical and chemical parts of the original system. Severalnumerical examples show that the closure model can capture the key rheological features of wormlikemicellar solutions in shear flows.
1. Introduction
Wormlike micelles, also known as “living polymers”, are long, cylindrical aggregates of self-assembledsurfactants that can break and reform reversibly [8]. There are substantial interests in studying thewormlike micelles for fundamental research and industrial applications [10, 11, 63]. In particular, it hasbeen observed that many wormlike micellar solutions exhibit shear banding, where the material splitsinto layers with different viscosities when undergoing strong shearing deformation [50]. Theoretically,the shearing banding is thought to arise from a non-monotone rheological constitutive curves of the shearstress versus the applied shear rate for steady homogeneous flow [22, 48]. Understanding this unusualrheological behavior of wormlike micelles has been a focus of many theoretical and experimental studies[59, 23].During last a couple of decades, a number of mathematical models have been proposed for wormlikemicellar solutions [1, 8, 27, 49, 59]. Many theoretical models, such as Johnson-Segalman model [49] andRolie-Poly model [1], are one species models, which didn’t reflect the “living” nature of wormlike micelles.To account for the reversible breaking and reforming of micellar chains, Cates proposed a reptation-reaction model, in which the reaction kinetics is introduced to account for the reversible breaking andreforming process [8, 9]. Inspired by Cates’ seminal work, a two species, scission-reforming model forwormlike micellar solutions is proposed in [59], known as the the Vasquez-Cook-McKinley (VCM) model.Although the VCM model was derived from a highly simplified discrete version of Cates’ model [23], itcan capture the key rheological properties of wormlike micellar solutions [59, 56, 67]. As pointed outin [24], the VCM model is thermodynamically inconsistent, due to the assumption that the break ratedepends on the velocity gradient explicitly. The VCM model was later revisited into a thermodynamicallyconsistent form [23, 22], which is known as the Germann–Cook–Beris (GCB) model, by using generalizedbracket approach [5]. They also extended such model to a three species cases. Using the GENERICframework [28, 55], Grmela et al. also extend the VCM model into a multiple species model [27].From a modeling perspective, both the VCM and GCB models are phenomenological macroscopicmodels, which cannot provide detailed information on polymer structures. For many complex fluids,
Email addresses: [email protected] (Yiwei Wang), [email protected] (Teng-Fei Zhang), [email protected] (ChunLiu)
Preprint submitted to Elsevier wo-scale macro-micro models, which couple the evolution of the microscopic probability distributionfunction of polymeric molecules, with the macroscopic flow, have been proven successful in describingtheir dynamics [37, 38, 41]. In these models, the micro-macro interaction is coupled through a transport ofthe microscopic Fokker-Planck equation and the induced elastic stress tensor in the macroscopic equation.The competition between the kinetic energy and the multiscale elastic energies give arise to differenthydrodynamical the rheological properties. The goal of this paper is to extend such a micro-macroapproach to model wormlike micellar solution by incorporating the microscopic breaking and reformingreaction kinetics. Following the setting in the VCM model [59], we represent the wormlike micelles bytwo species of dumbbells respectively. Instead of constructing some empirical constitutive equation, weemploy an energetic variational approach to derive the governing equation from an prescribed energy-dissipation law.We also study maximum entropy closure approximations to the new micro-macro model. We adoptboth “closure-then-variation” and “variation-then-closure” approaches. The first approach, which hasbeen widely used in literature [32, 60], apply the maximum entropy closure at the PDE level, whilethe later approach first reformulates an energy-dissipation law at the course grained level and derivesthe closure system by a variation procedure [17]. Due to the presence of reaction kinetics, these twoapproaches are not equivalent. Although the “closure-then-variation” approach can obtain model with aenergy-dissipation property, the model fails to produce a non-monotone rheological constitutive curves ofthe shear stress versus the applied shear rate for steady homogeneous flow. In contrast, by formulatingthe dissipation part in the coarse-grained level properly, the “closure-then-variation” approach can resultin a model that preserves the thermodynamical structures of the mechanical and chemical parts of theoriginal system. The resulting closure system takes the same form as the VCM [59] and GCB models[23]. Numerical simulations show that the momentum closure model can capture the key rheologicalfeatures of wormlike micellar solutions as the VCM [59] and GCB [23] models.The paper is organized as follows. In section 2, we formally derive the micro-macro model for wormlikemicellar solutions by employing an energetic variational approach. A detailed investigation of maximumentropy closure approximations to the micro-macro model is presented in section 3. In section 4, weshow the closure model, obtained by “closure-then-variation” can capture the key rheological features ofwormlike micellar solutions in planar shear flows.
2. Energetic variational formation of the new micro-macro model
In this section, we employ an energetic variational approach to derive a thermodynamically consistenttwo-species micro-macro model for wormlike micellar solutions.
Originated from pioneering work of Onsager [51, 52] and Rayleigh [58], the energetic variationalapproach (EnVarA) provides a general framework to derive the dynamics of a nonequilibrium thermo-dynamic system from a prescribed energy-dissipation law through two distinct variational processes:the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP) [42, 25]. The energy-dissipation law, which comes from the first and second laws of thermodynamics [20, 25], can be formulatedas dd t E total ( t ) = −△ , (2.1)for an isothermal closed system. Here E total is the total energy, which is the sum of the Helmholtz freeenergy F and the kinetic energy K ; △ is the rate of energy dissipation, which is equal to the entropyproduction in this case. The LAP states that the dynamics of a Hamiltonian system is determinedas a critical point of the action functional A ( x ) = R T ( K − F )d t with respect to x (the trajectory inLagrangian coordinates for mechanical systems) [2, 25], i.e., δ A = Z T Z Ω( t ) ( f inertial − f conv ) · δ x d x d t. (2.2)2n the meantime, for a dissipative system, the dissipative force can be determined by minimizing thedissipation functional D = △ with respect to the “rate” x t in the linear response regime [15], i.e., δ D = Z Ω( t ) f diss · δ x t d x . (2.3)This principle is known as Onsager’s MDP [51, 52]. Thus, according to force balance (Newton’s secondlaw, in which the inertial force plays role of ma ), we have δAδ x = δ D δ x t (2.4)in Eulerian coordinates, which is the dynamics of the system. In the framework of EnVarA, the dynamicsof the system is totally determined by the choice of the energy-dissipation law and the kinematic relation,which shifts the main task of modeling to the construction of an energy-dissipation law. The EnVarAframework has been proved to be a powerful tool to build up thermodynamically consistent mathematicalmodels for many complicated system, especially those in complex fluids [42, 25].Complex fluids are fluids with complicated rheological phenomena, arising from the interaction be-tween the microscopic elastic properties and the macroscopic flow motions [41, 42]. A central problemin modeling complex fluids is to construct a constitutive relation, which links the stress tensor τ andthe velocity field ∇ u [36]. Unlike a Newtonian fluid, there is no simple linear relation τ = µ ˙ γ , where˙ γ = ( ∇ u + ∇ u T ) is the strain rate and µ is the viscosity, for complex fluids. Instead of constructingan empirical constitutive equation that often takes the form of ∂ t τ + ( u · ∇ ) τ = f ( τ , ∇ u ) , (2.5)the EnVarA framework derives the constitutive relation from the giving energy-dissipation law throughthe variation procedure. Hence, the multiscale coupling and competition among multiphysics can bedealt with systematically.As an illustration, we first give a formal derivation of a one-species incompressible micro-macro modelof dilute polymeric fluid by employing the EnVarA. A more detailed description to this model can befound in [41, 25]. In this model, it is assumed that the polymeric fluid consists of beads joined by springs,and a molecular configuration is described by an end-to-end vector q ∈ R d [7, 18]. At the microscopiclevel, the system is described by a Fokker-Planck equation of the number distribution function ψ ( x , q , t )with a drift term depending on the macroscopic velocity u . While the macroscopic motion of the fluidis described by a Navier-Stoke equation with an elastic stress depending on the ψ ( x , q , t ).To derive the dynamics of the system by the EnVarA, we need to introduce Lagrangian descriptionsin both microscopic and macroscopic scales. In the macroscopic domain Ω, we define the flow map x ( X , t ) : Ω → Ω, where X are Lagrangian coordinates and x are Eulerian coordinates. For fixed X , x ( X , t ) is the trajectory of a particle labeled by X , while for fixed t , x ( X , t ) is an orientation-preservingdiffeomorphism between the initial domain to the current domain. For a given flow map x ( X , t ), we candefine the associated velocity u ( x ( X , t ) , t ) = dd t x ( X , t ) , (2.6)and the deformation tensor e F ( x ( X , t ) , t ) = F ( X , t ) = ∇ X x ( X , t ) . (2.7)Without ambiguity, we will not distinguish F and e F in the following. It is easy to verify that F ( x , t )satisfies the transport equation [41] F t + u · ∇ F = ∇ u F , (2.8)in Eulerian coordinates, where u · ∇ F stands for u k ∂ k F ij . The deformation tensor F carries all thekinematic information of the microstructures, patterns, and configurations in complex fluids [39]. Similarto the macroscopic flow map x ( X , t ), we can also introduce the microscopic flow map q ( X , Q , t ), where X and Q are Lagrangian coordinates in physical and configuration spaces respectively. The correspondingmicroscopic velocity V is defined as V ( x ( X , t ) , q ( X , Q , t ) , t ) = dd t q ( X , Q , t ) . (2.9)3ue to the conservation of mass, the number density distribution function ψ ( x , q , t ) satisfies the followingkinematics ∂ t ψ + ∇ x · ( u ψ ) + ∇ q · ( V ψ ) = 0 , (2.10)where u and V are effective velocities in the macroscopic domain the the microscopic configurationspace respectively. In the framework of EnVarA, the micro-macro system can be modeled through theenergy-dissipation lawdd t Z Ω (cid:20) ρ | u | + λ Z ψ (ln ψ −
1) + ψU d q (cid:21) d x = − Z Ω (cid:20) η |∇ u | + Z R d λξ ψ | V − ∇ uq | d q (cid:21) d x , (2.11)where K = R Ω 12 ρ | u | d x is the kinetic energy, λ is a constant that represents the ratio between thekinetic energy and the elastic energy, U ( q ) is the spring potential energy, and ξ is the constant thatis related to the polymer relaxation time. We assume that u satisfies the incompressible condition ∇ · u = 0. The second-term in the dissipation accounts for the relative friction of microscopic particle tothe macroscopic flow, where ( ∇ u ) q is velocity induced by the macroscopic flow due to the Cauchy-Bornrule [41]. The Cauchy-Born rule states that the movement in configuration space follows the flow on themacroscopic level, i.e., q = F Q without the microscopic evolution, where Q are Lagrangian coordinatesin the configuration space. A direct computation shows that e V = dd t ( F Q ) = ∇ u F Q = ∇ uq , (2.12)which is the microscopic velocity induced by the macroscopic flow.Now we are ready to perform the energetic variational approach in both microscopic and macroscopicscales. It is import to keep the “separation of scales” in mind when applying the LAP and the MDPin both scales. On the microscopic scale, since x ( X , t ) is treated being independent from q ( Q , X , t ), astandard energetic variational approach results in1 ξ ψ ( V − ∇ uq ) = − ψ ∇ q (ln ψ + U ( q )) , (2.13)where the right-hand side is obtained by the LAP, taking the variation of −F q = − R ψ ln ψ + ψU d q withrespect to q , and the left-hand side is obtained by the MDP, taking the variation of D q = ξ R ψ | V −∇ uq | d q with respect to V [25, 44]. Combining (2.13) with the kinematics (2.10), we have ψ t + ∇ · ( u ψ ) + ∇ q · ( ∇ uq ψ ) = ξ (∆ q ψ + ∇ q · ( ψ ∇ q U )) . (2.14)On the macroscopic scale, due to the “separation of scales”, we treat q ( X , Q , t ) as being independentfrom x ( X , t ). The Cauchy-Born rule is taken into account by the dissipation term | V − ( ∇ u ) q | . Theaction functional is defined by A ( x ) = Z T Z Ω (cid:20) ρ | u | − λ Z R ( ψ (ln ψ −
1) +
U ψ )d q (cid:21) d x d t, (2.15)By the LAP, i.e., taking variation of A ( x ) with respect to x , we obtain δ A δ x = − ρ x tt = − ρ ( u t + u · ∇ u ) . (2.16)Meanwhile, for the dissipation part, the MDP results in δ D δ x t = − η ∆ u + λξ ∇ · Z ψ ( V − ∇ uq ) ⊗ q d q . (2.17)Notice λξ ∇ · Z ψ ( V − ∇ uq ) ⊗ q d q = λ ∇ · Z ( −∇ q ψ ⊗ q − ∇ q U ⊗ q ψ ) d q = − λ (cid:18)Z ∇ q U ⊗ q ψ d q − n I (cid:19) , (2.18)4here the first equality is obtained by using (2.13). Thus, the force balance condition leads to themacroscopic momentum equation: ρ ( u t + u · ∇ u ) + ∇ p = η ∆ u + ∇ · τ , (2.19)where p is the Lagrangian multiplier for the incompressible condition ∇ · u , and τ is the induced elasticstress tensor given by τ = λ (cid:18)Z R ψ ∇ q U ⊗ q d q − n I (cid:19) . (2.20)The form of the induced elastic stress tensor τ is exactly the Kramers’ expression of the polymeric stress[38], which reflects the microscopic contribution to the macroscopic flow. Remark 2.1.
In the above derivation, the induced elastic stress tensor is derived from the dissipation partof the energy-dissipation law. Alternatively, one can derive the equivalent induced elastic stress tensorfrom the conservative part [41, 33]. Due to the Cauchy-Born rule, we can assume that the configurationspace follows the flow in the macroscopic scale, i.e., q = F Q and V = ∇ uq . Thus, the macroscopicaction functional can be defined by A ( x ) = Z T Z Ω ρ | x t | − λ Z ψ (ln ψ −
1) + U ( F Q ) ψ d Q d X (2.21) and the macroscopic dissipation is simply D = R η |∇ u | d x (the second term in the dissipation vanishessince the Cauchy-Born rule is used). By taking variation of A ( x ) with respect to x , we have [41] δ A δ x = − ρ x tt + λ ∇ · ( Z ψ ∇ q U ⊗ q d q ) . (2.22) Meanwhile, δ D δ x t = − η ∆ u . Hence, we end up with the same macroscopic equation with τ given by τ = λ Z R ψ ∇ q U ⊗ q d q , (2.23) which is equivalent to the (2.20) in the incompressible case since ∇ · ( − n I ) will contribute to the pressureand can be dropped [38]. The classic energetic variational approach, which is indeed based on classical mechanics, cannot beapplied to systems involving chemical reactions. Since 1950’s, many papers tried to developed a Onsagertype variational theory for chemical reactions by building analogies between Newtonian mechanics andchemical reactions [30, 53, 6, 4, 5, 47]. In a recent work [62], the authors proposed an energetic variationalformulation to chemical reactions by using the reaction trajectory R as the state variable, which enableus to couple the chemical reactions with other mechanical mechanisms in a unified energetic variationalframework. The reaction trajectory, also known as the extent of reaction or degree of advancement, wasoriginally introduced by De Donder [13, 14]. For a general reversible chemical reaction system containing N species { X , X , . . . X N } and M reactions can be represented by α l X + α l X + . . . α lN X N −− ⇀↽ −− β l X + β l X + . . . β lN X N , l = 1 , . . . , M. (2.24)We can define a reaction trajectory R ∈ R M , where R l accounts for the “number” of forward chemicalreactions that has occurred by time t . The relation between species concentration c ∈ R N + and thereaction trajectory R is given by c = c + σR , (2.25)where c is the initial concentration, and σ ∈ R N × M is the stoichiometric matrix with σ il = β li − α li .One can view (2.25) as the kinematics of the chemical reaction system [43]. With the kinematics (2.25),one can reformulate the free energy F , which is a functional of c , in terms of the reaction trajectory R [62, 43], Moreover, notice that δ F δR l = N X i =1 σ il δ F δc i = N X i =1 σ il µ i , (2.26)5hich is exactly the affinity of l − the chemical reaction, as defined by De Donder [14]. The affinityplays a role of the “force” that drives the chemical reaction, which vanishes at the chemical equilibrium[35]. It is worth mentioning that R is the internal state variable in [12], which is analogous to the flowmap x ( X , t ) in mechanical systems [53]. The reaction rate r is defined as ∂ t R , which can be viewedas the reaction velocity [35]. Similar to a mechanical system, the reaction rate can be obtained from aprescribed energy-dissipation law in terms of R and ∂ t R :dd t F [ R ] = −D chem [ R , ∂ t R ] , (2.27)where D chem [ R , ∂ t R ] is the rate of energy dissipation due to the chemical reaction procedure. Since thelinear response assumption for chemical system may be not valid unless at the last stage of chemicalreactions [5, 15], D chem is not quadratic in terms of ∂ t R in general. For a general nonlinear dissipation D chem [ R , ∂ t R ] = ( Γ ( R , ∂ t R ) , ∂ t R ) = M X l =1 Γ l ( R , ∂ t R ) ∂ t R l ≥ , the reaction rate can be derived as [62, 43]:Γ l ( R , ∂ t R ) = − δ F δR l , (2.28)which is the “force balance” equation for the chemical part [62, 43]. It is often assumed that Γ l ( R , ∂ t R ) =Γ l ( R l , ∂ t R l ). So equation (2.28) specify the reaction rate of l − the chemical reaction. In this formulation,the choice of the free energy determines the chemical equilibrium, while the choice of the dissipationfunctional D chem [ R , ∂ t R ] determines the reaction rate. Now we are ready to derive a thermodynamically consistent two-species micro-macro model for worm-like micellar solutions. Following the setting of the VCM model[59], we consider there exist only twospecies in the system. A molecule of species A can break into two molecules of species B , and twomolecules of species B can reform species A . At a microscopic level, both of polymer molecules are mod-eled as elastic dumbbells as in classical models of dilute polymeric fluids [7, 18]. We denote the numberdensity distribution of finding each molecule with end-to-end vector q at position x by ψ A ( x , q , t ) and ψ B ( x , q , t ) respectively. Then, the number density of species α is defined by n α ( x , t ) = Z ψ α d q . (2.29)Due to the breakage and reforming of polymer chains, in the microscopic scale, ψ A and ψ B satisfies thekinematics ( ∂ t ψ A + ∇ · ( u A ψ A ) + ∇ q · ( V A ψ A ) = − ∂ t R∂ t ψ B + ∇ · ( u B ψ B ) + ∇ q · ( V B ψ B ) = 2 ∂ t R, (2.30)where u α is the effective macroscopic velocity of each species, V α is the effective microscopic velocity,and R is the reaction trajectory for the breakage and reforming procedure. Throughout this paper,we disregard the diffusive effects of A and B , and assume u A = u B = u , which is the velocity of themacroscopic fluids. Moreover, we assume that the velocity field u satisfies the incompressible condition ∇ · u = 0. Similar to the one species micro-macro model, the total energy of the system can be writtenas E total = Z Ω h ρ | u | + λ Z ψ A (ln ψ A −
1) + ψ A U A + ψ B (ln ψ B −
1) + ψ B U B d q i d x , (2.31)where λ is the ratio between the macroscopic kinetic energy and microscopic elastic energy, and U α isthe potential energy associated with each species. The choice of U α also determines the microscopicequilibrium of the breakage and reforming process. Indeed, since the affinity vanishes at the chemicalequilibrium [14, 35], i.e., 2 µ B − µ A = 0 , (2.32)6here µ A and µ B are chemical potentials defined by µ A = δ F δψ A = ln ψ A + U A , µ B = δ F δψ B = ln ψ B + U B . (2.33)Hence, we have ln ψ ∞ A + U A = 2 ln ψ ∞ B + 2 U B or ln (cid:18) ψ ∞ A ( ψ ∞ B ) (cid:19) = 2 U B − U A . (2.34)at the chemical equilibrium. We can define the microscopic equilibrium constant K eq = ψ ∞ A ( ψ ∞ B ) = exp (2 U B − U A ) (2.35)In general, K eq depends on q unless 2 U B − U A is a constant.The dissipation part is also similar to the one-species model (2.11), which takes the form △ = − Z Ω h η |∇ u | + λ Z ψ A ξ A | V A − ∇ uq | + ψ B ξ B | V B − ∇ uq | + ∂ t R Γ( R, ∂ t R )d q i d x , (2.36)where ξ α is a constant related to the relaxation time of each species, and ∂ t R Γ( R, ∂ t R ) ≥ R, ∂ t R ) determinedifferent reaction kinetics. A typical choice of Γ( R, ∂ t R ) isΓ( R, ∂ t R ) = ln (cid:18) ∂ t Rη ( R ) + 1 (cid:19) , (2.37)Recall (2.28), we have Γ( R, ∂ t R ) = µ A − µ B , (2.38)which implies that ∂ t R = η ( R ) (cid:18) exp ( − ( U B − U A )) ψ A ψ B − (cid:19) (2.39)If η ( R ) = k ( q ) ψ B , (2.39) can be further simplified as ∂ t R = k ( q ) ψ A − k ( q ) ψ B , (2.40)where k ( q ) = k ( q ) K eq ( q ) , which can be viewed as law of mass action at the microscopic level. Remark 2.2.
From a modeling perspective, Eq. (2.40) corresponds to the case that an A molecule atposition x with end-to-end vector q can only break into two B molecules with same end-to-end vector,and the reforming process can only happen between two B molecules at the same position x with thesame end-to-end vector. We should emphasize that the assumption here is only for the mathematicalsimplicity, which may not fully reflect the complicated physical scenario. In the original VCM model, theauthor assumes that ∂ t R = k ψ A − k ψ B ∗ ψ B , (2.41) where ψ B ∗ ψ B = Z ψ B ( x , q − ˜ q .t ) ψ ( x , ˜ q , t )d˜ q (2.42) The advantage of the assumption in the VCM model is that the system will satisfies the law of massaction for n A and n B in the macroscopic scale, that is ( n A + ∇ · ( n A u A ) = − k n A + k n B n B + ∇ · ( n B u B ) = 2 k n A − k n B . (2.43) So it would be easy to get a closed macroscopic equation under the assumption (2.41). However, it seemsto difficult to obtain a variational structure for the breakage and reforming mechanism (2.41).
The derivation for the dynamics in the mechanical part of the two-species model is almost same to7hat in the one-species case. In the microscopic scale, a standard EnVarA leads to ψ α ∇ q µ α = − ξ α ψ α ( V α − ( ∇ u ) q ) (2.44)that is ψ α V α = − ξ α ( ψ α ∇ q U α + ∇ q ψ α ) + ( ∇ u ) q ) ψ α (2.45)Hence, the microscopic equation is given by ( ∂ t ψ A + u · ∇ ψ A + ∇ q · ( ∇ uq ψ A ) − ξ A ∇ q · ( ∇ q ψ A + ∇ q U A ψ A ) = − ∂ t R,∂ t ψ B + u · ∇ ψ B + ∇ q · ( ∇ uq ψ B ) − ξ B ∇ q · ( ∇ q ψ B + ∇ q U B ψ B ) = 2 ∂ t R, (2.46)where ∂ t R is defined in (2.40). On the macroscopic scale, similar to the one species case, by an energeticvariational approach, we can obtain ρ ( u t + ( u · ∇ ) u ) + ∇ p = η ∆ u + λ ∇ · τ (2.47)where τ is the induced stress from the microscopic configurations τ = Z ψ A ∇ q µ A ⊗ q d q + Z ψ B ∇ q µ B ⊗ q d q = Z ( ∇ q U A ⊗ q ψ A + ∇ q U B ⊗ q ψ B ) d q − ( n A + n B )I (2.48)Hence, the final macro-micro system is given by ρ ( ∂ t u + ( u · ∇ ) u ) + ∇ p = η ∆ u + λ ∇ · τ ∇ · u = 0 ∂ t ψ A + u · ∇ ψ A + ∇ q · ( ∇ uq ψ A ) − ξ A ∇ q · ( ∇ q ψ A + ∇ q U A ψ A ) = − ∂ t R∂ t ψ B + u · ∇ ψ B + ∇ q · ( ∇ uq ψ B ) − ξ B ∇ q · ( ∇ q ψ B + ∇ q U B ψ B ) = 2 ∂ t R (2.49)where ∂ t R = k ( q ) ψ A − k ( q ) ψ B , (2.50)and τ is the stress tensor given by (2.48). According to the previous derivation, it is easy to show thatthe system satisfies the following energy-dissipation law:dd t Z (cid:20) ρ | u | + λ Z ψ A (ln ψ A − U A ) + ψ B (ln ψ B − U B )d q (cid:21) d x = − Z (cid:20) η |∇ u | + λξ A Z ψ A |∇ q (ln φ A + U A ) | d q + λξ B Z ψ B |∇ q (ln ψ B + U B ) | d q + λ Z ( k ( q ) ψ A − k ( q ) ψ B ) ln (cid:18) k ( q ) ψ A k ( q ) ψ B (cid:19) d q (cid:21) d x , (2.51)
3. Moment closure models
The micro-macro model (2.49) provides a thermodynamically consistent multi-scale description towormlike micellar solutions. However, it might be difficult to study this model directly, as the microscopicequation (2.46) is high dimensional. Notice that the macroscopic stress tensor only involves the zerothand second moments of the number distribution functions of two species, it is a natural idea to derivea coarse-grained macroscopic equation from the original micro-macro model through moment closure.Moment closure has been a powerful tool to obtain coarse-grained macroscopic constitutive equationsfrom more detailed micro-macro models for complex fluids [16, 19, 21, 54, 61]. One challenge in momentclosure is to preserve the thermodynamic structures, i.e., the coarse-grained system should satisfy aenergy-dissipation law analogous to the energy-dissipation law of the original system [54, 32]. Thepresence of the chemical reaction imposes additional difficulties for closure approximations.8hroughout this section, we assume the potential energy U α to be U A = 12 H A | q | + σ A , U B = 12 H B | q | + σ B , (3.1)where σ A and σ B are constants related to the equilibrium of the breakage and reforming procedure, H A and H B are Hookean spring constants associated with species A and B . Moreover, we assume that H A = 2 H B , (3.2)then K eq = exp(2 σ B − σ A ) is a constant, which enables us to have a model with both k and k beingconstants. Same assumption is used in the GCB model [23]. Other types of potential energies can beconsidered but will result in more complicated closure systems. With the assumption H A = 2 H B , theglobal equilibrium distribution of the system is given by ψ ∞ A = n ∞ A (2 πH − A ) d/ exp (cid:18) − H A q T q (cid:19) , ψ ∞ B = n ∞ B (2 πH − B ) d/ exp (cid:18) − H B q T q (cid:19) , (3.3)where n ∞ A and n ∞ B are number densities at the global equilibrium. Correspondingly, the second momentsat the global equilibrium are given by A eq = Z q ⊗ q ψ ∞ A d x = n ∞ A H A I , B eq = Z q ⊗ q ψ ∞ B d x = n ∞ B H B I . (3.4)Let K macro eq = n ∞ A ( n ∞ B ) be the macroscopic equilibrium constant and a direct computation shows that K eq = e σ B − σ A = ψ ∞ A ( ψ ∞ B ) = 2 d π d/ H d/ B K macro eq , (3.5)which reveals the connection between K eq and K macro eq . Remark 3.1.
In the original VCM model [59], it is assumed that H B = 3 kTl N B , H A = 3 kTl N A = 3 kT l N B , (3.6) where N α is the number of Kuhn steps of length l in the segment of species α . So H B = 2 H A . Thisassumption is consistent with the case that two spring is connected in series. While the assumption H A =2 H B corresponds to the case that two spring is connected in parallel. The reason we take H A = 2 H B isonly for the mathematical simplicity.3.1. Maximum entropy closures Maximum entropy closures, also known as quasi-equilibrium approximations [26, 32, 60] have beensuccessfully used to derive effective macroscopic equations from the micro-macro multi-scale models forpolymeric fluids, including nonlinear dumbbell models [60, 32] and liquid crystal polymers [3, 34, 29, 65].For nonlinear dumbbell models with FENE potential, it has been shown that maximum entropy closurecan capture the hysteretic behavior and maintain the energy-dissipation property [32, 60].The idea of the maximum entropy closure is to maximize the “relative entropy” subjected to moments[26, 32, 54, 60]. For our system, we can approximate ψ α ( α = A, B ) based on its zeroth moment n α andsecond moment M α by solving the constrained optimization problem ψ ∗ α = argmin A Z R d ψ ln ψ + ψU α ( q )d q , (3.7)where A = { ψ : R d → R , ψ ≥ | Z ψ d q = n α , Z ( q ⊗ q ) ψ d q = M α } . (3.8) Proposition 3.1.
For the Hookean potential U α = H α | q | + σ α , the minimization problem (3.7) hasa unique minimizer ψ ∗ α in the class A for given n α > and a symmetric positive-definite matrix M α . oreover, ψ ∗ α is given by ψ ∗ α ( q ) = n α (2 π ) d/ (det f M α ) / exp( − q T f M − α q ) , where f M α = M α /n α . We call ψ ∗ α is the quasi-equilibrium state associated with n α and M α .Proof. The solution to the constrained optimization problem (3.7) is given by δδψ n Z ψ ln ψ + U α ( q ) ψ d q + λ (cid:20)Z ψ d q − n α (cid:21) + X ij λ ij (cid:20)Z q i q j ψ d q − ( M α ) ij (cid:21)o = 0 , (3.9)where λ and λ ij are Lagrangian multipliers. From (3.9), one can obtain that ψ ∗ α = C exp( − H α | q | − σ α ) exp − λ − X ij λ ij q i q j , (3.10)where C > R ψ ∗ α d x = n α , ψ ∗ α can be written as ψ ∗ α = n α Z ( λ ij ) exp − H α | q | − X ij λ ij q i q j . (3.11)where Z ( λ ij ) = R exp( − H α | q | ) exp (cid:16) − P ij λ ij q i q j (cid:17) d q is the normalizing constant. Since ψ ∗ α /n α isthe multivariate normal distribution N (0 , Σ ) with the covariance matrix given by Σ = ( H α I + 2 Λ ) − , (3.12)which is uniquely determined by its second moment, i.e., Σ = ( H α I + 2 Λ ) − = M α /n α [31].Thus, for given n A > n B >
0, positive-definite matrices A and B , we can define the uniquequasi-equilibrium states ψ ∗ A = n A (2 π ) d/ (det e A ) / exp (cid:18) − q T e A − q (cid:19) .ψ ∗ B = n B (2 π ) d/ (det e B ) / exp (cid:18) − q T e B − q (cid:19) , (3.13)where e A = A /n A and e B = B /n B are conformation tensors [23]. We call the manifold formed by allquasi-equilibrium distributions as quasi-equilibrium manifold , denoted by M ∗ = ( ψ ∗ = n (2 π ) d/ (det f M ) / exp (cid:18) − q T f M − q (cid:19) | n > , f M symmetric, positive-definite ) (3.14)For ψ ∗ α ∈ M ∗ , its second moment M α = n α f M α depends on its zeroth moment n α . We can apply the maximum entropy closure to the micro-macro model (2.49) directly. Since k and k are constants, by integrating (2.30) over q , we have ( ∂ t n A + ∇ · ( n A u ) = − k n A + k R ψ B d q ∂ t n B + ∇ · ( n B u ) = 2 k n A − k R ψ B d q . (3.15)Meanwhile, multiplying both side of (2.30) by q ⊗ q arrives at ( ∂ t A + ( u · ∇ ) A − ( ∇ u ) A − A ( ∇ u ) T = ξ A (2 n A I − H A A ) − k A + k R q ⊗ q ψ B d q ∂ t B + ( u · ∇ ) B − ( ∇ u ) B − B ( ∇ u ) T = ξ B (2 n B I − H B B ) + 2 k A − k R q ⊗ q ψ B d q . (3.16)10herefore, for the Hookean spring potential, the moment closure is needed only due to the nonlinearreaction term in the microscopic scale. With the maximum entropy closure (3.13), these two terms canbe computed out explicitly. Indeed, notice that Z ( ψ ∗ B ) d q = Z R d n B Z B exp( − q T e B − q )d q , (3.17)by letting q = √ e q , we have Z ( ψ ∗ B ) d q = 2 − d/ Z R d n B Z B exp( − e q T e B − e q )d e q = n d/ B d π d/ (det B ) / n B . (3.18)Hence, Z k ψ ∗ A − k ( ψ ∗ B ) d q = k n A − e k ( B ) n B , (3.19)where e k ( B ) is given by e k ( B ) = n d/ B d ( π ) d/ (det( B )) / k . (3.20)Interestingly, in this case, the maximum entropy closure gives us the law of mass action on numberdensities A k −− ⇀↽ −− e k k Z ( ψ ∗ B ) ( q ⊗ q )d q = ( √ − d Z R d n B Z B exp( − e q T e B − e q ) ˜ q ⊗ ˜ q d ˜ q = 12 e k ( B ) n B B . (3.21)Therefore, applying the maximum entropy approximation to (2.49), we can obtain a moment closuresystem ρ ( ∂ t u + ( u · ∇ ) u ) + ∇ p = η ∆ u + λ ∇ · ( H A A + H B B − ( n A + n B ) I ) ∇ · u = 0 ∂ t n A + ∇ · ( n A u ) = − k n A + e k ( B ) n B ,∂ t n B + ∇ · ( n B u ) = 2 k n A − e k ( B ) n B ,∂ t A + ( u · ∇ ) A − ( ∇ u ) A − A ( ∇ u ) T = 2 ξ A ( n A I − H A A ) − k A + e k ( B ) n B B ∂ t B + ( u · ∇ ) B − ( ∇ u ) B − B ( ∇ u ) T = 2 ξ B ( n B I − H B B ) + 2 k A − e k ( B ) n B B . (3.22)where e k ( B ) = n d/ B d ( π ) d/ (det( B )) / k . This is the model obtained by the “variation-then-closure”, i.e., applying the maximum entropy closureat the PDE level. One can prove that the closure system (3.22) possesses an energy-dissipation law. Toshow this, we first look at the case with u = 0. Proposition 3.2.
In absence of the flow field u = 0 , given n A > , n B > and symmetric, positive-definite matrices A and B , the closure system (3.22) satisfies the energy-dissipation law dd t F CL ( n A , n B , A , B ) = −△ CL ≤ , (3.23)11 here F CL ( n A , n B , A , B ) is the coarse-grained free energy given by F CL ( n A , n B , A , B ) = Z n A (cid:18) ln (cid:18) n A n ∞ A (cid:19) − (cid:19) + n B (cid:18) ln (cid:18) n B n ∞ B (cid:19) − (cid:19) − n A (cid:18) H A A n A (cid:19) + 12 tr( H A A − n A I ) − n B (cid:18) H B B n B (cid:19) + 12 tr( H B B − n B I )) d x , (3.24) and △ CL is the rate of energy dissipation, given by △ CL = Z ξ A tr (cid:0) ( H A I − n A A − ) A (cid:1) + ξ B tr (cid:0) ( H B I − n B B − ) B (cid:1) + ( k n A − e k ( B ) n B ) ln (cid:18) n A n ∞ A (cid:19) − (cid:18) n B n ∞ B (cid:19) + ln det( H B B /n B ) p det( H A A /n A ) ! + tr (cid:18) ( k A − e k ( B ) n B B )( n B B − − n A A − ) (cid:19) d x . (3.25) In particular, under the condition that n A > , n B > , and A and B are symmetric positive-definite, △ CL ≥ . Remark 3.2.
The coarse-grained free energy F CL ( n A , n B , A , B ) is same to the macroscopic free energygiven in [23]. The free energy contains two part: the Oldroyd-B type elastic energy associated with species A and B [31, 66], and the Lyapunov function of the chemical reaction A k −−− ⇀↽ −−− e k on number densitywith k n ∞ A = e k eq ( n ∞ B ) and e k eq2 = H d/ / (2 d π d/ ) .Proof. We first show that we have the identity (3.23) if n A , n B , A and B satisfy equation (3.22) with u = 0. Indeed, for F CL ( n A , n B , A , B ), a direct computation leads todd t F CL = dd t Z n A (ln ( n A /n ∞ A ) −
1) + n B (ln ( n B /n ∞ B ) − − n A ln det ( H A A /n A ) / H A A − n A I ) / − n B ln det ( H B B /n B ) / H B B − n B I )) / x . = Z (ln n A − ln n ∞ A − ln det( H A A ) / d ln( n A ) / ∂ t n A + (ln n B − ln n ∞ B − ln det( H B B ) / d ln( n B ) / ∂ t n A + tr(( H A I − n A A − ) ∂ t A )) / H B I − n B B − ) ∂ t B )) / x . (3.26)Substituting (3.22) into (3.26), and rearranging term, we havedd t F CL = − Z ξ A tr (cid:0) ( H A I − n A A − ) A (cid:1) + ξ B tr (cid:0) ( H B I − n B B − ) B (cid:1) + ( k n A − e k ( B ) n B ) ln (cid:18) n A n ∞ A (cid:19) − (cid:18) n B n ∞ B (cid:19) + ln det H B ( B /n B ) p det( H A A /n A ) ! + tr (cid:16) ( k A − e k ( B ) n B B / n B B − − n A A − / (cid:17) d x . To prove △ CL ≥
0, we first define the quasi-equilibrium state ψ ∗ A and ψ ∗ B for given n A > n B > A and B . The existence and uniqueness of ψ ∗ A and ψ ∗ B havebeen shown in proposition 3.1. Notice that for ψ ∗ α ( q ) = n α p (2 π ) d det( M ) exp( − q T M − α q ) . ψ ∞ α = n ∞ α q (2 π ) d H − dα exp( − H α q T q ) , (3.27)12e have Z ψ ∗ α ln (cid:18) ψ ∗ α ψ ∞ α (cid:19) d q = Z ψ ∗ α ln n α n ∞ α + ln 1 p det( H α M α ) + 12 (cid:0) − q T M − α q + H α q T q (cid:1)! d q = n α ln n α n ∞ α − n α H α M α ) + tr( H α n α M α − n α I )and Z ψ ∗ α (cid:12)(cid:12)(cid:12)(cid:12) ∇ q (cid:18) ln ψ ∗ α ψ ∞ α (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) d q = Z ψ α (cid:12)(cid:12)(cid:12)(cid:12) ∇ q ( − q T M − α q + 12 H α q T q ) (cid:12)(cid:12)(cid:12)(cid:12) d q = tr( − M − α + H α I ) n α M α (3.28)Moreover, by using the fact that k ψ ∞ A = k ( ψ ∞ B ) , we have Z ( k ψ ∗ A − k ( ψ ∗ B ) ) (cid:18) ln ψ ∗ A ψ ∞ A − ψ ∗ B ψ ∞ B (cid:19) d q = Z ( k ψ ∗ A − k ( ψ ∗ B ) ) " ln n A n ∞ A + ln 1 q det( H A e A ) + 12 (cid:16) − q T e A − q + H A q T q (cid:17) − n B n ∞ B + ln 1 q det( H B e B ) + 12 (cid:16) − q T e B − q + H B q T q ) (cid:17) = ( k n A − e k ( B ) n B ) ln (cid:18) n A n ∞ A (cid:19) − (cid:18) n B n ∞ B (cid:19) + ln det H B e B q det( H A e A + 12 tr (cid:18) ( − e A − + H A I + 2 e B − − H B I )( k A − k ( B ) n B e B ) (cid:19) , (3.29)where the last equality follows (3.19) and (3.21). Using H A = 2 H B and combining the above calculations((3.2), (3.28) and (3.29)), we can show the (3.23) is exactly same todd t Z Z ψ ∗ A (cid:18) ln (cid:18) ψ ∗ A ψ ∞ A (cid:19) − (cid:19) + ψ ∗ B (cid:18) ln (cid:18) ψ ∗ B ψ ∞ B (cid:19) − (cid:19) d q d x = − Z Z ξ A ψ A (cid:12)(cid:12)(cid:12)(cid:12) ∇ q (cid:18) ln (cid:18) ψ ∗ A ψ ∞ A (cid:19)(cid:19)(cid:12)(cid:12)(cid:12)(cid:12) + ξ B ψ B (cid:12)(cid:12)(cid:12)(cid:12) ∇ q (cid:18) ln (cid:18) ψ ∗ B ψ ∞ B (cid:19)(cid:19)(cid:12)(cid:12)(cid:12)(cid:12) + ( k ψ ∗ A − k ( ψ ∗ B ) ) ln (cid:18) k ψ ∗ A k ( ψ ∗ B ) (cid:19) d q d x , (3.30)which is obtained by replacing ψ α by ψ ∗ α in the original micro-macro energy-dissipation law (2.51). It isclear that the right-hand side of (3.30) is nonnegative, i.e., △ CL ≥ t (cid:18)Z ρ | u | d x + F CL ( n A , n B , A , B (cid:19) = − (cid:18)Z η |∇ u | d x + △ CL (cid:19) ≤ n A > n B > A and B . However, it is not straightforwardto derive the equation (3.22) from the the energy-dissipation law (3.31). Moreover, due to presence of thereaction procedure, the dynamics (3.22) no longer lies on the quasi-equilibrium manifold M ∗ . Indeed,the maximum entropy closure only use the information of the free energy part of the original system, itis unclear whether it is suitable for the dissipation part. As discussed in the next section, the closuremodel (3.22) fails to procedure a non-monotonic curves of the shear stress versus the applied shear ratefor steady homogeneous flow. Such a closure approximation may only valid when the elastic part reachesits equilibrium much faster then the reaction part in the original system, i.e., the solution will move to M ∗ rapidly [26]. Unfortunately, in a high shear rate region, the macroscopic flow prevents the elastic13art to reach its equilibrium. To obtain a thermodynamically consistent macroscopic model that suitable for the high shear rateregions, in this section, we consider a different closure approximation procedure, known as closure-then-variation . The idea is to apply the closure approximation to the energy dissipation law first, and derivethe closure system by applying the energetic variational approach in the coarse-grained level. Thisapproach is similar to the Onsager principle based dynamic coarse graining method proposed in [17].By imposing a proper dissipation mechanism on the quasi-equilibrium manifold M ∗ , we can have athermodynamically consistent closure model for both mechanical and chemical part of the system.On the quasi-equilibrium manifold M ∗ , we have A = n A e A and B = n B e B . So the free energy F CL ( n A , n B , A , B ) for the closure system, defined in (3.24), can be reformulated in terms of numberdensity n A and n B , and the conformation tensor of two species e A and e B , given by e F CL ( n A , n B , e A , e B ) = Z n A (cid:18) ln (cid:18) n A n ∞ A (cid:19) − (cid:19) + n B (cid:18) ln (cid:18) n B n ∞ B (cid:19) − (cid:19) + n A (cid:16) − ln det (cid:16) H A e A (cid:17) + tr (cid:16) H A e A − I (cid:17)(cid:17) + n B (cid:16) − ln det (cid:16) H B e B (cid:17) + tr (cid:16) H B e B − I (cid:17)(cid:17) (3.32)We can impose the kinematics for the number density to account for the macroscopic breakage andreforming procedure: ( ∂ t n A + ∇ · ( n A u ) = − ∂ t R n ∂ t n B + ∇ · ( n B u ) = 2 ∂ t R n . (3.33)where R n is the macroscopic reaction trajectory.The dissipation of the macroscopic moment closure system on M ∗ consists of three parts: the viscosityof the macroscopic flow, the evolution of the conformation tensors and the reaction on the number density,which can be formulated as e △ ∗ = Z η |∇ u | + tr M A d e A d t ! + tr M B d e B d t ! d x + e D chem ( R n , ∂ t R n ) (3.34)where d • d t = ∂ t • +( u · ∇ ) • − ( ∇ u ) • − • ( ∇ u ) T is the kinematic transport of the conformation tensor [40], M A ( n A , e A ) and M B ( n B , e B ) are mobility matrices. e D chem ( R n , ∂ t R n ), defined by e D chem ( R n , ∂ t R n ) = ∂ t R n ln ( η n ( R n ) ∂ t R n + 1) . (3.35)is the dissipation for breakage and reforming process at the macroscopic scale. The choice of η n ( R n )determines the macroscopic reaction rate in the closure system. One can view (3.34) as a projection ofthe original dissipation on the quasi-equilibrium manifold M ∗ . We then apply the energetic variationalapproach to obtain the dynamics on M ∗ , i.e, the moment closure system. The chemical reaction on the number density : By performing energetic variational approach withrespect to R n and ∂ t R n , we obtainln ( η n ( R n ) ∂ t R n + 1) = − δ e F CL δR n = µ nA − µ nB . (3.36)For the closure energy e F CL [ n A , n B , e A , e B ], we can compute the corresponding chemical potential ofnumber density n A and n B as µ nA = ln n A − ln n ∞ A −
12 ln (cid:16) det (cid:16) H A e A (cid:17)(cid:17) + 12 tr (cid:16) H A e A − I (cid:17) ,µ nB = ln n B − ln n ∞ B −
12 ln (cid:16) det (cid:16) H B e B (cid:17)(cid:17) + 12 tr (cid:16) H B e B − I (cid:17) , (3.37)which is same to the generalized chemical potential defined in [23].14t the chemical equilibrium for given e A and e B , µ nA = 2 µ nB , we have K neq eq = n neq A ( n neq B ) = n ∞ A exp( − tr( τ A /n A )) r det (cid:16) H A e A (cid:17) ( n ∞ B ) exp( − tr( τ B /n B )) det (cid:16) H B e B (cid:17) , (3.38)where τ A = H A A − n A I , τ B = H B B − n B I (3.39)is the induced stress tensor associated with species A and B respectively. Following [23], we take1 /η n ( R n ) = e k exp(tr( τ B ) /n B ) / det( H B e B ) n B , (3.40)which gives k neq1 = k eq exp( tr( τ A /n A )) q det( H A e A ) , k neq2 = e k eq exp(tr( τ B /n B )det( H B e B ) . (3.41)Thus, the number densities satisfy ∂ t n A + ∇ · ( n A u ) = − k neq1 n A + k neq2 n B ,∂ t n B + ∇ · ( n B u ) = 2 k neq1 n A − k neq2 n B , (3.42)The resulting non-equilibrium reaction rates of number densities is exact same to those in the GCBmodel [22, 23]. Gradient flows with convection on conformation tensors : The evolution of conformation tensorcan be obtained by performing energetic variational approach in terms of A ( B ) and d A d t ( d B d t ) [25, 66],which result in M A ( ∂ t e A + ( u · ∇ ) e A − ( ∇ u ) e A − e A ( ∇ u ) T ) = − δ F CL δ A M B ( ∂ t e B + ( u · ∇ ) e B − ( ∇ u ) e B − e B ( ∇ u ) T ) = − δ F CL δ B . (3.43)By taking M A = n A e A − / ξ A and M B = n B e B − / ξ B , we have ∂ t e A + ( u · ∇ ) e A − ( ∇ u ) e A − e A ( ∇ u ) T = ξ A (2 I − H A e A ) ∂ t e B + ( u · ∇ ) e B − ( ∇ u ) e B − e B ( ∇ u ) T = ξ B (2 I − H B e B ) . (3.44) Macroscopic flow equation:
Now we compute the macroscopic flow equation by performing theenergetic variational approach with respect to the flow map x ( X , t ). When writing the macroscopicforce balance, we should assume that the number densities and the conformation tensors to be purelytransported with flow. Under the incompressible condition (det F = 1), we have the kinematics [40] e A = F e A F T , e B = F e B F T , n A = n A , n B = n B , (3.45)and the action functional for the moment closure system can be given by e A [ x ] = Z T Z ρ | x t | − λ (cid:20) n A (cid:16) H A F e A F T (cid:17) + n B (cid:16) H B F e B F T (cid:17)(cid:21) d X d x (3.46)after dropping all the constant terms. A direct computation results in δ e A δ x = − ρ ( u t + u · ∇ u ) + λ ∇ · ( H A A + H B B )) . (3.47)The only dissipation term for the macroscopic flow is the viscosity part D η = R η |∇ u | d x [25], so the15issipative can be computed as δ D η δ x t = − η ∆ u . The final macroscopic force balance can be written as ρ ( u t + u · ∇ u ) + ∇ ˜ p = η ∆ u + λ ∇ · ( H A A + H B B )) , (3.48)where ˜ p is a Lagrangian multiplier for the incompressible condition. Equation (3.48) is equivalent to ρ ( u t + u · ∇ u ) + ∇ p = η ∆ u + λ ∇ · ( H A A + H B B − ( n A + n B ) I ) , (3.49)due to the incompressible condition.Finally, we get the the closure system ρ ( ∂ t u + ( u · ∇ ) u ) + ∇ p = η ∆ u + λ ∇ · ( H A A + H B B − ( n A + n B ) I ) ∇ · u = 0 ∂ t n A + ∇ · ( n A u ) = − k neq1 n A + k neq2 n B ,∂ t n B + ∇ · ( n B u ) = 2 k neq1 n A − k neq2 n B ,∂ t e A + ( u · ∇ ) e A − ( ∇ u ) e A − e A ( ∇ u ) T = 2 ξ A ( I − H A e A ) ∂ t e B + ( u · ∇ ) e B − ( ∇ u ) e B − e B ( ∇ u ) T = 2 ξ B ( I − H B e B ) , (3.50)where k neq1 and k neq2 are defined in (3.41) . One can view (3.50) as a dynamics restricted in the quasi-equilibrium manifold M ∗ . Recall that A = n A e A and B = n B e B on M ∗ . Combining (3.44) with (3.42),we have the second moment equations ∂ t A + ( u · ∇ ) A − ( ∇ u ) A − A ( ∇ u ) T = 2 ξ A ( n A I − H A A ) − k neq1 A + k neq2 n B e A ∂ t B + ( u · ∇ ) B − ( ∇ u ) B − B ( ∇ u ) T = 2 ξ B ( n B I − H B B ) + 2 k neq1 n A e B − k neq2 n B B (3.51) Remark 3.3.
We notice that the reaction terms in the celebrate VCM [59] and GCB models [23] takesa different form. As mentioned in remark 2.2, the VCM model assume the microscopic reaction takesthe form k ψ A − k ψ B ∗ ψ B from (3.51), which leads to the term k A A − k n B B in the second momentequation. The GCB model also take such a form as a starting point. To obtain the same form of reactionterms, one need further assume e A = e B , then − k neq1 A + k neq2 n B e A ≈ − k neq1 A + 12 k neq2 n B B , k neq1 n A e B − k neq2 n B B ≈ k neq1 A − k neq2 n B B . (3.52) The assumption (3.52) is reasonable, since for given number densities n A and n B , we have e A eq = H A I , e B eq = H B I , which implies that e A eq = e B eq at the local equilibrium. So e A ≈ e B is valid at leastnear the local equilibrium. Under the approximation (3.52), we can reach a closure model ρ ( ∂ t u + ( u · ∇ ) u ) + ∇ p = η ∆ u + λ ∇ · ( H A A + H B B − ( n A + n B ) I ) ∇ · u = 0 ∂ t n A + ∇ · ( n A u ) = − k neq1 n A + k neq2 n B ,∂ t n B + ∇ · ( n B u ) = 2 k neq1 n A − k neq2 n B ,∂ t A + ( u · ∇ ) A − ( ∇ u ) A − A ( ∇ u ) T = 2 ξ A ( n A I − H A A ) − k neq1 A + k neq2 n B B ∂ t B + ( u · ∇ ) B − ( ∇ u ) B − B ( ∇ u ) T = 2 ξ B ( n B I − H B B ) + 4 k neq1 A − k neq2 n B B , (3.53) which has the same form of the VCM and GCB models. Although the dynamics (3.53) no longer lies onthe quasi-equilibrium, it can produce more reasonable shear-stress curve at in a high shearing rate region.Compare with (3.50), (3.53) can force (cid:12)(cid:12) e A − e B (cid:12)(cid:12) to be small due to the approximation (3.52). We willcompare these two models in details in the future work. Remark 3.4.
In the above derivation, we assume the second moments can be written as the multiplicationof number density and the conformation tensor, i.e. A = n A e A and B = n B e B . Such a decomposition isvalid on the submanifold formed by the quasi-equilibrium states, but may not true in general. A differentmoment closure system can be obtained if one treat the number densities and the second moments to be ndependent. Then the free energy of the closure system (3.24) can be written as F CL ( n A , n B , A , B ) = Z n A (cid:18) ln (cid:18) n A n ∞ A (cid:19) −
12 det H A A − (cid:19) + n B (cid:18) ln (cid:18) n B n ∞ B (cid:19) − −
12 det( H B B ) (cid:19) + d n A ln n A − n A + n B ln n B − n B ) + 12 tr( H A A ) + 12 tr( H B B ) , (3.54) which implies that µ nA = ln n A − ln n ∞ A −
12 det( H A A /n A ) , µ nB = ln n B − ln n ∞ B −
12 det( H B B /n B ) . (3.55) We can simply modify the reaction rate in (3.22) by k neq1 = k eq / det( H A A /n A ) , k neq2 = k eq / det( H B B /n B ) (3.56) to obtain another closure model. We’ll explore this in the future. Remark 3.5.
It is worth mentioning that the derivation in this section can be viewed as a pure macro-scopic approach to model wormlike micellar solutions in the framework of EnVarA starting with the freeenergy F [ n A , n B , e A , e B ] and the dissipation given by (3.32) and (3.34) respectively. As a pure macro-scopic approach, it is not necessary to assume H A = 2 H B . If we assume H B = 2 H A , and adopt theapproximation A /n A ≈ B /n B , then we have − k neq1 A + k neq2 n B e A ≈ − k neq1 A + 2 k neq2 n B B , k neq1 n A e B − k neq2 n B B ≈ k neq1 A − k neq2 n B B , (3.57) which is exactly same to those in the GCB model [23, 22].
4. Numerics
In this section, we discuss the prediction of the above moment closure models through a few toyexamples. Detailed numerical studies for the original micro-macro model and the closure models will becarried out in future work.
First we consider a steady homogeneous shear flow with the velocity field given by u = ( κy, , where κ is the constant shear rate. Moreover, we assume that the number density of each species isspatial homogeneous. So the original PDE system is reduced to an ODE system of n α , A and B . Wesolve the ODE system by the standard explicit Euler scheme. We take the initial condition as n A = 1 , n B = 2 . , A = n A H A I B = n B H B I . (4.1)For each κ , we compute the number density of each species, and the induced shear stress τ = H A A + H B B . We compare the predictions for three models (3.22), (3.50) and (3.53). The terminal criterionfor the numerical calculation is T = 2 or n A < − . Fig. 4.1 showed the calculated shear stress andthe number densities of two species as a function of κ for three models ( k eq = 0 . e k eq = 0 . ξ A = 0 . ξ A /ξ B = 1 e − H A = 2, H B = 1). At small shear rates, all three models can produce similar results, dueto the fact that e A and e B are close to their equilibria. The closure model (3.22), obtained by applying themaximum entropy closure to the equation directly, fails to obtain a non-monotonic shear-stress curve.The main reason might be the fact that the break rate k is independent with the shear rate in thismodel, which can not lead to a pronounced breakage of species A . The predictions of model (3.50) and(3.53) are also different in the high shear rate region. The model (3.50) leads to a rapidly breakageof species A (Fig. 4.1 (c) - (d)), which does not seem to match previous experimental and simulation17 -2 -1 -4 -2 shear rate s hea r s t r e ss (a) -2 -1 shear rate nu m be r den s i t y (b) -2 -1 -4 -3 -2 -1 shear rate s hea r s t r e ss (c) -2 -1 shear rate nu m be r den s i t y (d) -2 -1 -4 -3 -2 -1 shear rate s hea r s t r e ss (e) -2 -1 shear rate nu m be r den s i t y (f) Figure 4.1:
Calculated shear stress and species number density as a function of κ ( k eq = 0 . e k eq = 0 . ξ A = 0 . ξ A /ξ B = 10 − , H A = 2, H B = 1). (a)-(b) Model (3.22); (c)-(d) Model (3.50); (e)-(f) Model (3.53). results [23]. The curves produced by the model 3.53, shown in Fig. 4.1(e)-(f), is consistent with theresults by the VCM and GCB models qualitatively [23]. As mentioned earlier, the approximation (3.52)can be viewed as an implicit regularization term such that | e A − e B | to be small, which prevent e A tobe too large. This simple numerical test show the importance of choosing a proper dissipation in thecourse-grained level in order to capture the non-equilibrium rheological properties of wormlike micellarsolutions. A detailed comparison of different closure models will be made in future work. In this subsection, we investigate the transient behavior of the model in a planar shear flow for theclosure model (3.53). Let u = ( u ( y ) ,
0) and u ( y ) satisfies ( u t = η∂ yy u + λ∂ y ( H A A + H B B ) ,u ( l ) = κ ( t ) , u (0) = 0 , (4.2)18e take κ ( t ) = γ tanh( at ), where a is a parameter control how the wall velocity approaches steady-state[67, 22]. Other parameters are set as: l = 0 . H A = 2, H B = 1, η = 1, λ = 1, k = 1 and e k eq2 = 6 . γ = 50 through this subsection. S hea r s t r e ss (a) S hea r S t r e ss (b) N u m be r den s i t y (c) Figure 4.2:
Transient behavior of the closure model (3.53) in a planar Couette flow: (a) Calculated shear stress at movingwall for different ramp-up rate a . (b) The wall shear stress as a function of t , (c) Species number density at the wall. Fig. 4.2(a) shows the transient response of the wall shear stress tensor for different ramp-up rates a . In all three cases, the shear stress will reach its maximum during the ramp-up process. Differentramp-up rates do not significantly affect the steady-state. Fig. 4.2 (b) shows temporal evolution of thetotal stress at the moving surface for κ ( t ) = 50 tanh(5 t ), the individual contributions of species of A and B are represented by dashed and dash-dotted lines. The number densities of species A and B are plottedin Fig. 4.2 (c). The above results are qualitatively agree with rheological characteristics predicted bythe GCB model in circular Taylor-Couette flow (see Fig. 4 and Fig. 5 in [22]).
5. Summary
In this paper, inspired the celebrated VCM type models [59, 23], we derive a thermodynamicallyconsistent two-species micro-macro model of wormlike micellar solutions by employing an energetic vari-ational approach. Our model incorporates a breaking and reforming process of polymer chains into aclassical micro-macro dumbbell model for polymeric fluids in unified variational framework. The ener-getic variational formulation for the micro-macro model opens a new door for both numerical studiesand theoretical analysis [46]. New structure-preserving numerical methods can be developed based onthe energetic variational formulation by using the approach in [43, 44, 45]. The modeling approachalso provides a framework to integrate other mechanism, and can be adopted to other chemomechanicalsystems beyond the wormlike micellar solutions, such as active soft matter systems [57, 64].We also study the maximum entropy closure approximation to the micro-macro model of wormlikemicellar solutions. In particular, we construct closure approximations by both “variation-then-closure”and “closure-then-variation” approaches. We show that these two approaches result in different closuremodels due to presence of the chemical reaction. In principle, maximum entropy closure only uses theinformation from the free energy part of the original system [26]. Applying the closure approximationon the PDE level cannot guarantee the thermodynamical consistency. By a “closure-then-variation”approach, we can restrict the dynamics on the coarse-grained manifold by choosing the dissipationproperly. As a consequence, the closure system preserves the thermodynamical structures of the original19ystem for both chemical and mechanical parts. Several numerical examples show that the closuremodel, obtained by “closure-then-variation” can capture the key rheological features of wormlike micellarsolution. A detailed numerical study for our models will be carried out in future work.
Acknowledgement
Y. Wang and C. Liu are partially supported by the National Science Foundation (USA) NSF DMS-1759536 and the United States-Israel Binational Science Foundation (BSF)
References [1] Adams, J., Fielding, S.M., Olmsted, P.D., 2011. Transient shear banding in entangled polymers: A study using therolie-poly model. J. Rheol. 55, 1007–1032.[2] Arnol’d, V.I., 2013. Mathematical methods of classical mechanics. volume 60. Springer Science & Business Media.[3] Ball, J.M., Majumdar, A., 2010. Nematic liquid crystals: from Maier-Saupe to a continuum theory. Molecular crystalsand liquid crystals 525, 1–11.[4] Bataille, J., Edelen, D., Kestin, J., 1978. Nonequilibrium thermodynamics of the nonlinear equations of chemicalkinetics. J. Non-Equilib. Thermodyn. 3, 153–168.[5] Beris, A.N., Edwards, B.J., Edwards, B.J., 1994. Thermodynamics of flowing systems: with internal microstructure.36, Oxford University Press on Demand.[6] Biot, M.A., 1977. Variational-lagrangian irreversible thermodynamics of initially-stressed solids with thermomoleculardiffusion and chemical reactions. J. Mech. Phys. Solids 25, 289–307.[7] Bird, R.B., Armstrong, R.C., Hassager, O., 1987. Dynamics of polymeric liquids. Vol. 1: Fluid mechanics. John Wileyand Sons Inc.[8] Cates, M., 1987. Reptation of living polymers: dynamics of entangled polymers in the presence of reversible chain-scission reactions. Macromolecules 20, 2289–2296.[9] Cates, M., 1990. Nonlinear viscoelasticity of wormlike micelles (and other reversibly breakable polymers). J. Phys.Chem. 94, 371–375.[10] Cates, M., Candau, S., 1990. Statics and dynamics of worm-like surfactant micelles. J. Phys.: Condens. Matter 2,6869.[11] Cates, M.E., Fielding, S.M., 2006. Rheology of giant micelles. Advances in Physics 55, 799–879.[12] Coleman, B.D., Gurtin, M.E., 1967. Thermodynamics with internal state variables. The Journal of Chemical Physics47, 597–613.[13] De Donder, T., 1927. L’affinit´e. M´emoires de la Classe des sciences. Acad´emie royale de Belgique. Collection in 8 9,1–94.[14] De Donder, T., 1936. Thermodynamic theory of affinity. volume 1. Stanford university press.[15] De Groot, S.R., Mazur, P., 2013. Non-equilibrium thermodynamics. Courier Corporation.[16] Doi, M., 1981. Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropicand liquid crystalline phases. Journal of Polymer Science: Polymer Physics Edition 19, 229–243.[17] Doi, M., 2016. A principle in dynamic coarse graining–Onsager principle and its applications. The European PhysicalJournal Special Topics 225, 1411–1421.[18] Doi, M., Edwards, S.F., Edwards, S.F., 1988. The theory of polymer dynamics. volume 73. oxford university press.[19] Du, Q., Liu, C., Yu, P., 2005. Fene dumbbell model and its several linear and nonlinear closure approximations.Multiscale Modeling & Simulation 4, 709–731.[20] Ericksen, J.L., 1998. Introduction to the Thermodynamics of Solids. volume 131 of
Applied Mathematical Sciences .Springer, New York.[21] Feng, J., Chaubal, C., Leal, L., 1998. Closure approximations for the Doi theory: Which to use in simulating complexflows of liquid-crystalline polymers? J. Rheol. 42, 1095–1119.[22] Germann, N., Cook, L., Beris, A., 2014. Investigation of the inhomogeneous shear flow of a wormlike micellar solutionusing a thermodynamically consistent model. J. Non-Newtonian Fluid Mech. 207, 21–31.[23] Germann, N., Cook, L., Beris, A.N., 2013. Nonequilibrium thermodynamic modeling of the structure and rheology ofconcentrated wormlike micellar solutions. J. Non-Newtonian Fluid Mech. 196, 51–57.[24] Germann, N., Kate Gurnon, A., Zhou, L., Pamela Cook, L., Beris, A.N., Wagner, N.J., 2016. Validation of constitutivemodeling of shear banding, threadlike wormlike micellar fluids. J. Rheol. 60, 983–999.[25] Giga, M.H., Kirshtein, A., Liu, C., 2017. Variational modeling and complex fluids. Handbook of mathematical analysisin mechanics of viscous fluids , 1–41.[26] Gorban, A.N., Karlin, I.V., Ilg, P., ¨Ottinger, H.C., 2001. Corrections and enhancements of quasi-equilibrium states.J Non-Newtonian Fluid Mech 96, 203–219.[27] Grmela, M., Chinesta, F., Ammar, A., 2010. Mesoscopic tube model of fluids composed of worm-like micelles. Rheol.Acta 49, 495–506.
28] Grmela, M., ¨Ottinger, H.C., 1997. Dynamics and thermodynamics of complex fluids. i. development of a generalformalism. Physical Review E 56, 6620.[29] Han, J., Luo, Y., Wang, W., Zhang, P., Zhang, Z., 2015. From microscopic theory to macroscopic theory: a systematicstudy on modeling for liquid crystals. Arch. Ration. Mech. Anal. 215, 741–809.[30] Horn, F., Jackson, R., 1972. General mass action kinetics. Archive for rational mechanics and analysis 47, 81–116.[31] Hu, D., Lelievre, T., et al., 2007. New entropy estimates for the oldroyd-b model and related models. Commun. Math.Sci. 5, 909–916.[32] Hyon, Y., Carrillo, J.A., Du, Q., Liu, C., 2008. A maximum entropy principle based closure method for macro-micromodels of polymeric materials. Kinetic and Related Models 1, 171–184.[33] Hyon, Y., Liu, C., et al., 2010. Energetic variational approach in complex fluids: maximum dissipation principle.Discrete & Continuous Dynamical Systems-A 26, 1291.[34] Ilg, P., Karlin, I.V., Kr¨oger, M., ¨Ottinger, H.C., 2003. Canonical distribution functions in polymer dynamics.(ii).liquid-crystalline polymers. Physica A: Statistical Mechanics and its Applications 319, 134–150.[35] Kondepudi, D., Prigogine, I., 2014. Modern thermodynamics: from heat engines to dissipative structures. John Wiley& Sons.[36] Le Bris, C., Lelievre, T., 2009. Multiscale modelling of complex fluids: a mathematical initiation, in: Multiscalemodeling and simulation in science. Springer, pp. 49–137.[37] Le Bris, C., Lelievre, T., 2012. Micro-macro models for viscoelastic fluids: modelling, mathematics and numerics.Science China Mathematics 55, 353–384.[38] Li, T., Zhang, P., 2007. Mathematical analysis of multi-scale models of complex fluids. Commun. Math. Sci. 5, 1–51.[39] Lin, F., 2012. Some analytical issues for elastic complex fluids. Commun. Pure Appl. Math. 65, 893–919.[40] Lin, F.H., Liu, C., Zhang, P., 2005. On hydrodynamics of viscoelastic fluids. Communications on Pure and AppliedMathematics 58, 1437–1471.[41] Lin, F.H., Liu, C., Zhang, P., 2007. On a micro-macro model for polymeric fluids near equilibrium. Commun. PureAppl. Math. 60, 838–866.[42] Liu, C., 2009. An introduction of elastic complex fluids: an energetic variational approach, in: Multi-Scale Phenomenain Complex Fluids: Modeling, Analysis and Numerical Simulation. World Scientific, pp. 286–337.[43] Liu, C., Wang, C., Wang, Y., 2020. A structure-preserving, operator splitting scheme for reaction-diffusion equationsinvolving the law of mass action. arXiv preprint arXiv:2010.16320 .[44] Liu, C., Wang, Y., 2020a. On Lagrangian schemes for porous medium type generalized diffusion equations: a discreteenergetic variational approach. J. Comput. Phys. , 109566.[45] Liu, C., Wang, Y., 2020b. A variational Lagrangian scheme for a phase-field model: A discrete energetic variationalapproach. SIAM J. Sci. Comput. 42, B1541–B1569. doi: .[46] Liu, C., Wang, Y., Zhang, T.F., 2021. On a two-species micro-macro model for wormlike micellar solutions: dynamicstability analysis. submitted .[47] Mielke, A., 2011. A gradient structure for reaction–diffusion systems and for energy-drift-diffusion systems. Nonlin-earity 24, 1329.[48] Mohammadigoushki, H., Dalili, A., Zhou, L., Cook, P., 2019. Transient evolution of flow profiles in a shear bandingwormlike micellar solution: experimental results and a comparison with the VCM model. Soft Matter 15, 5483–5494.[49] Olmsted, P., Radulescu, O., Lu, C.Y., 2000. Johnson–segalman model with a diffusion term in cylindrical couetteflow. J. Rheol. 44, 257–275.[50] Olmsted, P.D., 2008. Perspectives on shear banding in complex fluids. Rheol. Acta 47, 283–300.[51] Onsager, L., 1931a. Reciprocal relations in irreversible processes. I. Physical review 37, 405.[52] Onsager, L., 1931b. Reciprocal relations in irreversible processes. II. Physical review 38, 2265.[53] Oster, G.F., Perelson, A.S., 1974. Chemical reaction dynamics. Archive for rational mechanics and analysis 55,230–274.[54] ¨Ottinger, H.C., 2009. On the stupendous beauty of closure. J Rheol 53, 1285–1304.[55] ¨Ottinger, H.C., Grmela, M., 1997. Dynamics and thermodynamics of complex fluids. II. illustrations of a generalformalism. Physical Review E 56, 6633.[56] Pipe, C., Kim, N., Vasquez, P., Cook, L., McKinley, G., 2010. Wormlike micellar solutions: II. comparison betweenexperimental data and scission model predictions. J. Rheol. 54, 881–913.[57] Prost, J., J¨ulicher, F., Joanny, J.F., 2015. Active gel physics. Nature physics 11, 111–117.[58] Strutt, J., 1871. Some general theorems relating to vibrations. Proceedings of the London Mathematical Society 1,357–368.[59] Vasquez, P.A., McKinley, G.H., Cook, L.P., 2007. A network scission model for wormlike micellar solutions: I. modelformulation and viscometric flow predictions. J. Non-Newtonian Fluid Mech. 144, 122–139.[60] Wang, H., Li, K., Zhang, P., 2008. Crucial properties of the moment closure model FENE-QE. J. Non-NewtonianFluid Mech. 150, 80–92.[61] Wang, Q., 1997. Comparative studies on closure approximations in flows of liquid crystal polymers: I. elongationalflows. Journal of Non-Newtonian Fluid Mechanics 72, 141–162.[62] Wang, Y., Liu, C., Liu, P., Eisenberg, B., 2020. Field theory of reaction-diffusion: Law of mass action with an energeticvariational approach. Physical Review E 102, 062147.[63] Yang, J., 2002. Viscoelastic wormlike micelles and their applications. Current opinion in colloid & interface science 7,276–281.[64] Yang, X., Li, J., Forest, M.G., Wang, Q., 2016. Hydrodynamic theories for flows of active liquid crystals and thegeneralized Onsager principle. Entropy 18, 202.[65] Yu, H., Ji, G., Zhang, P., 2010. A nonhomogeneous kinetic model of liquid crystal polymers and its thermodynamicclosure approximation. Communications in Computational Physics 7, 383.[66] Zhou, J., Doi, M., 2018. Dynamics of viscoelastic filaments based on Onsager principle. Physical Review Fluids 3,084004.[67] Zhou, L., Cook, L.P., McKinley, G.H., 2012. Multiple shear-banding transitions for a model of wormlike micellar olutions. SIAM Journal on Applied Mathematics 72, 1192–1212.olutions. SIAM Journal on Applied Mathematics 72, 1192–1212.