Some Recent Advances in Energetic Variational Approaches.

In this paper, we summarize some recent advances related to the energetic variational approach (EnVarA), a general variational framework of building thermodynamically consistent models for complex fluids, by some examples. Particular focus will be placed on how to model systems involving chemo-mechanical couplings and non-isothermal effects.

1. Introduction

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 the elasticity of deformable particles, interaction between charged ions, and bulk elasticity endowed by polymer molecules [1,2]. Due to their strong nonlinear and non-equilibrium nature, building thermodynamically consistent models for complex fluids has been an interesting and challenging problem. The difficulty arises from complicated coupling and competition of different electro-chemo-mechanic mechanisms, such as long-range interaction and thermal fluctuation, in different spatio-temporal scales.

Motivated by the non-equilibrium thermodynamics, especially the seminal works of Rayleigh [3] and Onsager [4,5], the energetic variational approach (EnVarA) has proven to be a powerful tool in studying numerous complex fluids systems in physics, chemistry, and biochemistry [6,7], including liquid crystals [8,9], viscoelastic fluids [10], multiple-phase flows [11,12], ionic solutions [13], etc. The idea of EnVarA is to describe a complex system by an energy–dissipation law where is the sum of the kinetic energy and the Helmholtz free energy , is the rate of energy dissipation. The energy–dissipation law (1) can be obtained by combining the first and second laws of thermodynamics for an isothermal and mechanically isolated system [14]. Indeed, the first law of thermodynamics states that where is the kinetic energy, is the internal energy, is the rate at which heat absorbed from the environment, and is the rate of the external work done by the environment. To analyze heat, one introduces the second law of thermodynamics where T is the absolute temperature, S is the entropy, and is the rate of entropy production. By subtracting the first law (2) and the second law (3), we obtain an energy–dissipation law (1) for an isothermal () and mechanically isolated () system with being the Helmholtz free energy. So, in this case, the rate of energy dissipation equals the rate of entropy production.

Starting with an energy–dissipation law (1), EnVarA derives the dynamics of the systems through the least action principle (LAP) and the maximum dissipation principle (MDP). The LAP, which states the equation of motion for a Hamiltonian system can be derived by taking variation of the action functional with respect to the trajectory , gives a unique procedure to derive the conservative force in the system. The MDP, variation of the dissipation potential with respect to , i.e., the velocity, derives the dissipation force in the system. In turn, the force balance condition leads to the underlying evolution equation of the system According to the Onsager theory [4,5], the rate of entropy production is quadratic in terms of in the linear response regime, and the dissipation potential in this case.

There are other forms of variational principles, such as the general equation for non-equilibrium reversible–irreversible coupling (GENERIC) [1,15,16], Onsager’s variational principle [17,18,19,20], and conservation–dissipation formalism [21,22], that have also been helpful in studying complex fluids. Although these variational principles are equivalent to EnVarA in most cases, these approaches are based on the principle of virtual work (PVW), in which the variation is often taken with respect to the state variables directly. The state variables are defined in Eulerian coordinates in many cases, and can be viewed as generalized coordinates in these systems [17].

In contract to other variable principles, the classical EnVarA formulation is rooted in the continuum mechanics [23,24,25] and a Lagrangian formulation of an underlying system. Continuum mechanics is a generalization of Newtonian particle mechanics. In continuum mechanics, it is assumed that there are infinite many particles form a continuum body. The motion of these particles is described by a flow map , where is the Lagrangian coordinate and is the Eulerian coordinate. More importantly, in the context of continuum mechanics, all dynamics of the employed variables are determined by the flow map and its derivatives (the velocity and the deformation tensor) through their kinematics. The energy–dissipation law, together with the kinematics of the employed variables, describes all the physics and assumptions in the system. It is important to realize that systems derived by the LAP can also be derived by the PVW, but not vice versa. In fact, the discrepancy between the LAP the PVW is of great interest in the mathematical theory of weak solutions and the theory of singularities [8,11].

The purpose of this paper is to provide a brief overview of some recent developments of the energetic variational approach, especially for the systems involving chemo-mechanical couplings and non-isothermal effects, through examples. We refer interested readers to [6,7] for more details of the classical energetic variational approach. The rest of the paper is organized as follows. In Section 2, we briefly review the EnVarA in continuum mechanics, and show its applications in modeling generalized diffusion, a dilute polymeric fluids, and the kinetic Fokker–Planck equation. The EnVarA formulation for reaction kinetics and its applications in modeling chemo-mechanical systems, such as reaction–diffusion systems, the dynamical boundary conditions, the Boltzmann equation and the reaction/active fluids, are given in Section 3. In Section 4, we discuss the extension of the EnVarA to non-isothermal systems. Some applications of the EnVarA in developing structure-preserving numerical schemes and thermodynamically consistent coarse-grained model are reviewed in Section 5.

2. EnVarA in Continuum Mechanics

As mentioned in the introduction, the classical energetic variational approaches are variational principles for continuum mechanics [7], and the variable in (4) should be understood as the flow map from a reference domain to a physical domain . Here is the Lagrangian coordinate and is the Eulerian coordinate. For a fixed , describes a trajectory of a particle (or a material point) labeled by , while for a fixed t, is a diffeomorphism from to (see Figure 1 for an illustration).

Figure 1

An illustration of the flow map.

An important feature of a continuum mechanical system is that the evolution of physical variables, such as the density function, are determined by the evolution of the flow map through kinematics. All the kinematic transport information of these variables is carried by the deformation tensor . For a given flow map , the deformation tensor is defined by A direct computation shows that satisfies a transport equation [7] Without ambiguity, we will not distinguish and throughout this paper. Due to the conservation of mass, , the kinematics of a density function can be written as in Lagrangian coordinates, where is the initial density. The kinematics (6) is equivalent to the continuity equation in Eulerian coordinates. For a scalar variable that is purely transported by the flow map, the kinematics is given by in Eulerian coordinates. For a nematic liquid crystal that consists of rod–like molecule, the kinematics of the nematic order parameter can be written as

2.1. Generalized Diffusion

One of the simplest classes of mechanical processes is generalized diffusion. Generalized diffusion is concerned with a conserved quantity satisfying the kinematics where is the average velocity. In the framework of EnVarA, a generalized diffusion can be described an energy–dissipation law where is the free energy density, is the friction coefficient. Due to the kinematics (7), the free energy can be reformulated as a functional of in Lagrangian coordinates. A direct computation shows that where is the test function satisfying with being the outer normal of in Eulerian coordinates (Here we will not distinguish and without ambiguity). Pushing forward to Eulerian coordinates, we have which indicates that where is the chemical potential. In the notion of the principle of virtual work [24], one can obtain (9) by using the relation

For the dissipation part, since , it is easy to compute that . As a consequence, we have the force balance equation Combining the force balance Equation (10) with the kinematics (7), one can obtain a generalized diffusion equation It is worth mentioning that the above derivation is rather formal. A certain analysis is needed to show the existence of the flow map . We refer interested readers to [26,27] for some related discussions.

Many classical models can be viewed as generalized diffusions with different forms of the free energy and the dissipation. For example, the porous medium equation (PME) can be obtained by taking and , the Cahn–Hilliard equations can be obtained by taking and [28], and the Poisson–Nernst–Planck (PNP) equations for ion transport can be obtained from the energy–dissipation law along with the Poisson equation [13].

Another interesting example of generalized diffusions is the nonlinear Fokker–Planck equation associated with a random process It is well known that different definitions of stochastic integrals can lead to different diffusion equations. For the stochastic process (13), Itô calculus leads to a diffusion equation which corresponds to the energy–dissipation law if [29]. Meanwhile, the Stratonovich integral yields which corresponds to the energy–dissipation law if . In both cases, the condition , known as the fluctuation–dissipation theorem in statistical physics [30], is crucial for the existence of the energy–dissipation law, as well as the existence of an equilibrium state [7].

In the case that , generalized diffusions can be viewed as Wasserstein gradient flows in space of all probability densities having finite second moments [31]. Formally, the Wasserstein gradient flow can be defined as a continuous time limit () of the semi-discrete scheme, known as the JKO scheme, where and is the Wasserstein distance between and . The Wasserstein distance between two probability densities and can be computed through a Benamou–Brenier formulation [32] where is the admissible set.

The Wasserstein gradient flows are Eulerian descriptions to generalized diffusions [33]. Other choices of dissipation can define other metrics in the space of probability measures [34,35]. Within (15), the JKO scheme can be roughly approximated as under a constant velocity assumption, where . The scheme (17) can be viewed as an implicit Euler discretization to the flow map Equation (10) in . Mathematically, it is interesting to establish the equivalence or discrepancy between (14) and (17).

2.2. Micro-Macro Model for Polymeric Fluids

One successful application of EnVarA is building thermodynamically consistent micro-macro models for many complex fluids [36,37]. Compared with the macroscopic continuum mechanics approach, micro–macro models couple the macroscopic hydrodynamic equations with a microscopic kinetic theory, which describes the origin of the macroscopic stress tensor [38,39].

In the simplest micro–macro models of complex fluids, the polymer molecules are modeled as an elastic dumbbell consisting of two “beads” joined by a one-dimensional spring [36,37]. The microscopic configuration of an elastic dumbbell is described by an end-to-end vector between the two beads, . Let be the number density distribution function of finding a molecule with end-to-end vector at position at time t. To build a micro–macro model by EnVarA, in addition to the macroscopic flow map at the physical space, one also needs to introduce a flow map at the configurational space, denoted by , where are Lagrangian coordinates in the configurational space. For a given , the microscopic velocity can be defined as Due to the conservation of mass, the density distribution function satisfies which leads to the kinematics in Eulerian coordinates.

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, T 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 and respectively, where is the elastic constant and 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, due to the macroscopic flow, which indicates From the energy–dissipation law (21), one can derive the dynamics of the system by performing EnVarA in both micro- and macro scales. First, we look at the dynamics at the macroscopic scale. Due to the “separation of scale” [7], the second term in the dissipation (21) 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 the variation of with respect to , we obtain in Eulerian coordinates. Indeed, consider a perturbation where satisfying . Then Pushing forward to Eulerian coordinates, we have which leads to (23). For the dissipation part, the MDP, i.e., taking the variation of with respect to , leads to where p is the Lagrangian multiplier for the incompressible condition . Hence, the macroscopic force balance results in the momentum equation where is the induced stress from the configuration space, representing the microscopic contributions to the macroscopic level.

On the microscopic scale, similar to the generalized diffusions discussed in the last subsection, by taking variations with respect to and , we obtain Combining with Equation (20), we obtain the equation on the microscopic scale: In summary, the final micro–macro system reads as follows: subject to a suitable boundary condition.

2.3. Kinetic Fokker–Planck Equation

Using a similar micro-macro approach, one can also formulate the kinetic Fokker–Planck equation into an EnVarA formulation. The kinetic Fokker–Planck equation, given by describes the evolution of the probability density function of a system of particles that satisfies a Langevin dynamics where is the standard Wiener process and is the potential function. Due to the fluctuation–dissipation theorem (FDT), , which ensures that the system admits an energy–dissipation law and can reach an equilibrium state [40].

The total energy of the system can be defined by Then the kinetic Fokker–Planck equation can be written as two parts [41] where conserves the total energy , while dissipates . Indeed, a direct computation shows that Since where the last equality follows the fact that and . Hence, conserves the total energy. In the meantime, Here, we need to obtain the quadratic form. is known as the fluctuation–dissipation theorem [30]. This calculation also reveals the connection between the fluctuation–dissipation theorem and linear response theorem.

Similar to the micro–macro approach for polymeric fluids introduced in a previous subsection, one can treat as a microscopic variable, and as a macroscopic variable; the micro–macro coupling is imposed through . Due to the conservation property, the probability density satisfies the kinematics where can be decomposed as Here is a conserve velocity. In other words, the equation conserves the total energy, i.e., if . Notice that which indicates that . In other word, we can obtain by applying the least action principle at the macroscopic scale.

The dissipative velocity can be obtained through the energy–dissipation law A standard variational procedure results in Here plays a role of the internal energy in the macroscopic scale. The EnVarA formulation to the kinetic Fokker–Planck equation may open a new door for both theoretical and numerical study of these types of equations.

3. EnVarA for Chemical Reactions

Chemical reactions play important roles in many physical, chemical, and biological processes [42]. The fundamental “law” for chemical reaction kinetics is the law of mass action (LMA), which is discovered by Waage and Guldberg, and H. van’t Hoff independently in the 19th century. A reversible chemical reaction system with N species and M reactions can often be represented by where and are constant coefficients. Let be concentrations of all species. Then satisfies the reaction kinetics where is the reaction rate for the l-th chemical reaction, and is the stoichiometric coefficients. From (38), it is noticed that In turn, one can define N-rank linearly independent conserved quantities for the reaction network.

The classical LMA states that the reaction rate is directly proportional to the product of the reactant concentrations, i.e., in which and are the forward and backward reaction constants for the l-th reaction. The LMA was designed for the ideal gas of a perfect gas of non-interacting point particles without charge or diameter. To generalize the LMA to more complicated systems, it is important to establish a thermodynamics basis of this empirical law.

Since 1950s, there has been a large amount of work aiming to build an Onsager-type variational theory for systems involving chemical reactions [1,43,44,45,46,47,48,49,50,51,52].The key idea is to build on analogies between continuum mechanics and reaction kinetics [1,52,53]. Inspired by these prior works, in a recent paper [54], the authors extended the EnVarA formulation to a system that involves chemical reactions. By introducing the reaction trajectory , which accounts for the “number” of forward chemical reactions that has occurred by time t, the reaction kinetics can be reformulated in terms of . The reaction trajectory , which is known as the internal state variable in [43], is analogous to the flow map in mechanical systems [55]. The relation between species concentration and the reaction trajectory is given by where is the initial concentration, and is the stoichiometric matrix with . The reaction rate , defined as , is the reaction velocity [52]. Within the reaction trajectory, we can describe the chemical kinetics of a reaction network by an energy–dissipation law in terms of and : where is the rate of energy dissipation due to the chemical reaction.

Unlike mechanical processes, chemical reactions are usually far from equilibrium. The linear response assumption may be not valid unless at the last stage of chemical reactions [1,56]. As a consequence, is not quadratic in terms of in general. For a general nonlinear dissipation the reaction rate can be derived as [54,57]: which is the “force balance” equation for the chemical part [54,57]. It is often assumed that . Equation (43) specifies the reaction rate of the l-th chemical reaction. In this formulation, the free energy determines the chemical equilibrium, while the dissipation functional determines the reaction rate. It is worth mentioning that the variational principle (43) can be applied to other systems with non-quadratic dissipations, such as the systems studied in [58,59].

The classical LMA can be obtained by taking where is an equilibrium of the reaction kinetics system that satisfies the detailed balance condition, i.e., , and . Similar to the cases in [58,59], near a chemical equilibrium, is very small, and the dissipation of the LMA can be well approximated as a classical Onsager quadratic form of .

As an illustration, we consider a single chemical reaction in which . According to the previous discussion, the variational approach produces where From (46), we can obtain the law of mass action where .

The energetic variational formulation of chemical reactions opens a new door to model a general chemo-mechanical system in a unified variational way. One of the simplest chemo-mechanical systems is reaction–diffusion type system. For a reaction–diffusion system, the concentration of each species satisfies the kinematics: where is the average velocity of each species due to its own diffusion, represents various reaction trajectories involved in the system, with being the stoichiometric matrix as defined earlier. Then the reaction–diffusion equation can be modeled through the energy–dissipation law [54]: We can employ EnVarA to obtain equations for the reaction and diffusion parts, respectively, i.e., to obtain the “force balance equation” of the chemical and mechanical subsystems, which leads to By taking , we have a reaction–diffusion system One special reaction–diffusion system is a one-species system, in which the concentration or density satisfies the kinematics where is the average velocity that describes the diffusion of the species, and is the reaction trajectory that describes the birth and death process. We can model this type system by an energy–dissipation law These systems are related to the unbalanced optimal transport in machine learning [60,61,62,63,64]. Different systems can be derived by choosing the free energy and the dissipation differently, including the Kantorovich–Fisher–Rao gradient flows [65].

3.1. Dynamical Boundary Condition

For a better description of short-range interactions between the material and the boundary, PDE systems with dynamical boundary conditions have drawn lots of attention recently [66,67]. These type of systems have a wide application in studying moving contact lines [68], batteries [69], the integrate-and-fire model for neuron networks [70], and sticky Brownian motion [71], etc.

Consider a bounded domain , and denote the densities in as and the density on as . In the models with dynamical boundary conditions, the conservation of total mass indicates that Hence, and satisfy the kinematics in Eulerian coordinates where is the outer normal of and is the reaction trajectory for the chemical reaction that represents the density exchange between the bulk and surface.

In general, a system with a dynamical boundary condition can be modeled through an energy–dissipation law where and are free energies in the bulk and surface, respectively. The standard variational procedures lead to the force balance equations for the mechanical and chemical parts where is the affinity of the bulk–surface reaction. Different choices of the free energy and the dissipation lead to different systems.

A typical example is the Cahn–Hilliard equation with dynamic boundary condition, which can be derived by taking Here the linear response assumption is used for the bulk–surface exchange. The final equation can be written as [67], subject to a suitable initial condition. In the limit , the model is reduced to the model proposed in [66], while limit corresponds to the model in [72].

3.2. Boltzmann Equation

Another application of the EnVarA formulation for reaction kinetics is to reformulate the Boltzmann equation of ideal gas into a variational form, which was also explored in some previous works [73,74]. The key point is to view the collisions between particles as generalized chemical reactions. From a historical perspective, Maxwell and Boltzmann used an analogue of the law of mass action for collisions and discovered the principle of detailed balance [75].

Consider a Boltzmann equation where is the density for particles at the point , having velocity at time t, is the collision term that can be written as Here , and is the transition probability for two particles with velocity and before a collision to have velocities and after the collision.

The collisions can be viewed as generalized chemical reactions. Since collisions conserve mass, momentum, and kinetic energy, satisfies Moreover, due to the conservation of momentum and kinetic energy one can assume that only depends on and , where is a parameter that determines and , i.e., For given and , a fixed defines a “chemical reaction” with , where is the collision kernel. At the Maxwellians, vanishes identically, so the system satisfies the detailed balance condition.

Like the kinetic Fokker–Planck equation, the kinematics of can be written as the sum of the conservative and dissipative parts, i.e., where and . The total energy of the system can be defined as For the conservative part, For the dissipative part, a direct computation shows that Here, and are used.

We can reformulate the energy–dissipation law of the Boltzmann Equation (68) in terms of reaction trajectories , by imposing the kinematics The dissipation can be written in terms of and , i.e., Notice that By a variational procedure, we can obtain which leads to It is worth mentioning that is the affinity of the chemical reaction (64). Indeed, the chemical potential associated with can be computed as where plays the role of the internal energy. Due to the conservation of the kinetic energy, we have , which indicates that the affinity of the collision equals to .

3.3. Reactive Fluids

Reactive and active fluids, which are systems involving coupling and competition between chemical reactions and various types of mechanical processes, have drawn lots of attention recently [76,77,78,79,80,81].

A simple example of reactive complex fluids is wormlike micellar (WLM) solutions, also known as “living polymer”. These are polymeric fluids that consist of long, cylindrical aggregates of self-assembled surfactants that can break and reform reversibly [82]. During the last couple of decades, a number of mathematical models have been proposed for wormlike micellar solutions [82,83,84,85,86]. The reptation–reaction model, proposed by Cates [82], is one of the first models that accounts for the reversible breaking and reforming of micellar chains. The Cates’ model assumes that a chain can break with a fixed probability per unit time per unit length anywhere along its length. Let be the number density of chains of length L, and the governing equation of can be formulated as where and are rate constants for the breakage and recombination, respectively. It can be shown that the steady state () is given by where is a constant that satisfies

Inspired by Cates’ seminal work, a two-species model for wormlike micellar solutions was proposed in [86], known as the Vasquez–Cook–McKinley (VCM) model. Although the VCM model was derived from a highly simplified discrete version of Cates’ model [87], it can capture the key rheological properties of wormlike micellar solutions [86,88,89]. However, due to the assumption that the breakage rate depends on the velocity gradient explicitly, the VCM model may not be thermodynamically consistent [87]. Using the generalized bracket approach [1], a thermodynamically consistent revision to the VCM model was proposed in [87,90], known as the GCB (Germann–Cook–Beris) model. They also extended such a model to a three species cases. Using the GENERIC framework [15,16], Grmela et al. also extended the VCM model into a three-species thermodynamically consistent model [84].

In a recent work [91], the authors developed a thermodynamically consistent two-species micro–macro model for wormlike micellar solutions, which incorporates a breakage and combination process of polymer chains into the classical micro–macro dumbbell model of polymeric fluids in the general framework of EnVarA. The model assumes that there exists only two species in the system. A molecule of species A can break into two molecules of species B, and two molecules of species B can reform species A. A polymer molecule of both species can be modeled as an elastic dumbbell consisting of two “beads” joined by a one-dimensional spring [37]. The microscopic configuration 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 t for species .

In general, the breakage and combination processes can be regarded as chemical reactions where and are end-to-end vectors of species B, and is an end-to-end of species A (see Figure 2a for illustration). We can denote the forward and backward reaction rates of (75) by and , respectively. The kinematics of and can be written as where and are effective macroscopic and microscopic velocities. Different models can be obtained by choosing and differently. In the simplest case, one can take which corresponds to the case that an A molecule at position with end-to-end vector can only break into two B molecules with the same end-to-end vector, and the combination process can only happen between two B molecules at the same position with the same end-to-end vector, as illustrated in Figure 2b, with . Within this assumption, one can have a detailed balance condition for each and , and the kinematics can reduce to where is the reaction trajectory for the breakage and combination for given and .

Figure 2

Schematic diagram of breakage and combination processes in wormlike micellar solutions, in which different species are indicated by different colors. (a) General reaction mechanism (75); (b) The reaction mechanism (76) considered in this paper ().

Within the above kinematics assumptions, the overall system can be modeled by an energy–dissipation law where is the number density distributions of each species, is the spring potential associated with each species, is the macroscopic velocity satisfying the incompressible condition , the constant is the density of the macroscopic flow, is the average microscopic velocity of each species in the configuration space, is the constant that represents the ratio between the kinetic energy and the elastic energy, and is a constant related to the relaxation time of each species. By taking , the final micro–macro system can be derived as (see [91] for details) where and the stress tensor is given by with being the number density of each species. The breakage and combination process indeed created an active stress in the macroscopic momentum equation. The global existence of the class solution near the global equilibrium of this model is given in [92]. The underlying variational structure (78) plays a crucial role in the theoretical analysis.

4. EnVarA for Non-Isothermal Systems

The above EnVarA framework works for isothermal cases, in which the whole system can be well-described by a single energy-dissipation law (1). However, in non-isothermal cases, one cannot obtain a single energy–dissipation law from the first and second laws of thermodynamics. Motivated by the classical treatments of temperature in rational mechanics [14,43], some recent work [93,94,95,96] extended EnVarA to non-isothermal cases.

To model a non-isothermal system, the extended EnVarA approach starts with a suitable form of the non-isothermal Helmholtz free energy density , where is the absolute temperature and represents all the mechanical variables in the system, such as the density and the deformation tensor , which are determined by the flow map . The non-isothermal Helmholtz free energy density needs to be concave with respect to the temperature , i.e., [14].

According to the basic thermodynamical relations [14,97], the entropy density is defined by while the internal energy density is defined by Since , s is a monotonically increasing function of , which indicates that can be represented as a function of and s, i.e., As a consequence, the internal energy density e can be reformulated as a function of and s by using (84). Without ambiguity, we denote . According to the chain rule, it is straightforward to show and As in isothermal systems, it is crucial to specify the kinematics of the mechanical variable as well as the kinematics of the temperature in non-isothermal systems. Different kinematic relations lead to different dynamics. To illustrate the idea, in the following, we only consider a simple case, with the mechanical variable being a density function only, which satisfies the kinematics . For the temperature , we assume that i.e., the temperature is purely transported along the trajectory. Moreover, we focus on cases where the kinetic energy throughout this section.

Although the whole system is no longer determined by the energy–dissipation law (1), we can still define the total energy as well as the dissipation potential , and derive the force balance Equation (4) for the mechanical part of the system by using the LAP and MDP. To simplify the calculations, we assume that throughout this section, although other forms of the dissipation can be handled in a similar manner. By applying the LAP and the MDP, we can obtain the force balance equation for the mechanical part, i.e., equation for the flow map . The computation procedure is almost the same as that in Section 2.1. Indeed, a direct computation shows that Pushing forward to the Eulerian coordinates, we have which indicates that . Then the force balance equation for the mechanical part can be written as To determine the equation for temperature or entropy s, we need to use both the first and second laws of thermodynamics. The first law of the thermodynamics can be written as where V is an arbitrary control volume, stands for the specific work done to the system, and is the absorbed heat flux. According to (85) and (86), it is easy to compute that where the kinematic assumption on the entropy s is used. Here, is the entropy flux defined by the Clausius–Duham relation, and is the density of the entropy production. Comparing (90) with (89), we obtain that Recall (88). We have On the other hand, by the Fourier law, we have where . Hence, the density of the entropy production is which is consistent with the second law of thermodynamics.

The final system can be written as For ideal gas models, it is often assumed that Then In this case, the non-isothermal system (96) can be simplified as where the last equation follows (90).

One can derive a different system by imposing a different kinematics assumption on the temperature and the entropy. For instance, we can assume that the temperature is independent of the flow map during the purely mechanics processes, which leads to a simple kinematics . The force balance Equation (88) becomes . In this case, the kinematics of the entropy becomes Again, is the density of the entropy production, and is the entropy flux defined by Clausius–Duhem relation. Following the same procedure, one can show that the entropy production is same as in (95). We refer the interested reader to [98] for a detailed derivation.

The formulations (91) and (100) are the differential version of the entropy balance equation [52,99] introduced by Prigogine in [100] as a modern formulation of the second law of thermodynamics (See 3.4 in [52]). For mechanical processes, is the entropy in the control volume V, is the entropy change due to the exchange of matter and energy with being the entropy flux [52], while is the entropy change due to irreversible process, i.e., the entropy production. The Clausius–Duhem relation specifies .

For systems involving chemical reactions, since the system is not determined by the flow map alone, it is natural to assume that the temperature also depends on the reaction trajectory [43]. The remaining derivation is almost same as for the mechanical system.

5. Coarse-Graining and Numerical Realization

In this section, we briefly review some applications of the EnVarA in developing structure-preserving numerical schemes [28,57,101] and thermodynamically consistent coarse-grained models [91].

One of the fundamental questions in developing numerical schemes is how to preserve the properties of the original systems as much as possible. Different mathematical representations of physical principles can lead to different structure-preserving methods [102]. For dissipative systems modeled by the EnVarA, the most important structures are the underlying variational structures and the kinematic constraints of physical variables, such as positivity and the conservation property of a mass density, and the conservation of current for a charged system.

In [28], the authors introduced a numerical framework called the discrete energetic variational approach to preserve the energetic variational structure in the semi-discrete level. The approach first seeks a finite dimensional realization to a continuous energy–dissipation law, given by where is the discrete state variable, is the discrete total energy that is the sum of the discrete kinetic energy and the discrete free energy , and is the discrete energy dissipation. The discrete energy–dissipation (102) law can be obtained by introducing a suitable spatial discretization to the original system by Eulerian [103], Lagrangian [28,101], or particle approximations [104]. From the discrete energy–dissipation law (102), one can obtain an ODE system of , given by by the energetic variational approach. Here, is the discrete action functional. In the case that and is a quadratic function of , the ODE system of can be written as where is a matrix.

The discrete energetic variational approach follows the idea of “discretize-then-variation”, which has been used to study both Hamiltonian and dissipative systems [18,105,106,107,108,109,110]. An advantage of this Ritz-type variational approach is that the resulting system inherits the variational structure from the continuous level, which enables us to apply these variational temporal discretizations [111,112,113]. For example, an implicit Euler discretization to Equation (104) could be reformulated as where is the admissible set of and is the numerical mobility that is independent with . The scheme is also known as the minimizing movement scheme [114,115]. Although in general, the optimization problem (105) is non-convex in a non-convex admissible set , one can use some line search-based optimization method to update such that . The existence of a minimizer of that decreases the total energy and the convergence of the numerical scheme at the discrete level can be proved under suitable conditions on [101,112,116]. A smaller can be chosen such that the optimization problem can be solved by some standard optimization method.

Similarly, the BDF2 scheme or a modified Crank–Nicolson scheme can also be reformulated as a minimization problem due to the variational structure in the semi-discrete level [113,117]. Other temporal discretization techniques, such as convex splitting [118,119,120,121], invariant energy quadratization (IEQ) [122,123], scalar auxiliary variable (SAV) [124,125], the discrete variational derivative method [105,126], and the strong stability-preserving (SSP) time discretizations [127] can also be applied to the semi-discrete system. In some applications, if the goal is to compute equilibrium states and the phase diagram, one can use some optimization method, such as the L-BFGS, to minimize the discrete free energy directly. Using the approach, we studied various confined liquid crystal systems [103,128,129,130].

The idea of “approximation-then-variation” can also be used to develop a dynamic coarse-grained model from a detailed multiscale model [18]. Similar to numerical approximations, the aim of coarse-graining is to reduce the degrees of freedom of the original system but preserve the properties of the original systems as much as possible. Closure approximation is one typical approach to develop coarse-grained models. The idea of closure approximation in to use the evolution equation of moments to represent the evolution of the probability density in the original Fokker–Planck type equations.

In [91], the authors studied maximum closure approximations to the micro–macro model of a wormlike micellar solution (79) by using “closure-then-variation”. The idea is to apply the closure approximation to the energy dissipation law first. Let and be the second moment of and , and be the number density of species . By applying the energetic variational approach in the coarse-grained level, we can obtain a thermodynamically consistent closure model for both mechanical and chemical parts of the system. In more detail, within the assumption that and , the free energy, obtained by the maximum entropy closure, can be reformulated in terms of number density and , and the conformation tensor of two species and , given by We impose the kinematics for the number density to account for the macroscopic breakage and reforming procedure: where is the macroscopic reaction trajectory.

The dissipation of the moment closure system consists of three parts: the viscosity of the macroscopic flow, the evolution of the conformation tensors, and the reaction on the number density, which can be formulated as where is the kinematic transport of the conformation tensor, and are mobility matrices. is the dissipation for breakage and reforming at the macroscopic scale. A typical choice of is The choice of determines the macroscopic reaction rate in the closure system. By taking with , and applying the energetic variational approach, we can obtain the final moment closure system where , , and . One can view (110) as a dynamics restricted in the sub–manifold formed by quasi–equilibrium states, in which and . It is clear that the fluctuations on the number density and will create an active stress tensor in the macroscopic momentum equation. The resulting closure system is similar to the VCM [86] and GCB models [87]. Numerical simulations in [91] show that the momentum closure model can capture the key rheological features of wormlike micellar solutions as the VCM and GCB models.

6. Conclusions

In this paper, we briefly review some recent advances related to the energetic variational approaches, by examples. The general framework of EnVarA provides a thermodynamically consistent description to a complicated chemo-mechanical system involving couplings and competitions, such as the energetics vs. kinematics, the macroscopic hydrodynamics vs. micro-structures, the reversible vs. irreversible, the mechanics vs. chemistry, and the deterministic vs. stochastic. It is a powerful tool to study multiscale and multi-physics problem arising in physics, chemistry, and biology. The current review only sketches the underlying variational structure of several complex fluid models with chemo-mechanical coupling and non-isothermal effects; all these models require lots of investigation, both theoretically and numerically.