Application of Onsager Machlup integral in solving dynamic equations in non-equilibrium systems

TL;DR

This paper introduces an improved variational method based on the Onsager-Machlup integral for solving kinetic equations in non-equilibrium systems, overcoming previous limitations and enabling steady state determination.

Contribution

The paper presents a novel approach using the Onsager-Machlup integral that allows for better approximate solutions and steady state analysis in non-equilibrium kinetic equations.

Findings

Successfully applied to diffusion, capillary, and free boundary problems.

Enables determination of steady states via variational calculus.

Improves accuracy over previous small-time incremental methods.

Abstract

In 1931, Onsager proposed a variational principle which has become the base of many kinetic equations for non-equilibrium systems. We have been showing that this principle is useful in obtaining approximate solutions for the kinetic equations, but our previous method has a weakness that it can be justified, strictly speaking, only for small incremental time. Here we propose an improved method which does not have this drawback. The new method utilizes the integral proposed by Onsager and Machlup in 1953, and can tell us which of the approximate solutions is the best solution without knowing the exact solution. The new method has an advantage that it allows us to determine the steady state in non-equilibrium system by a variational calculus. We demonstrate this using three examples, (a) simple diffusion problem, (b) capillary problem in a tube with corners, and (c) free boundary problem in liquid coating, for which the kinetic equations are written in second or fourth order partial differential equations.

Equations

$$\frac{d{x}_i}{dt}=-\sum _j{\mu }_{ij}\frac{\partial A}{\partial x_j},$$

$$\displaystyle\sum_{ij}\mu_{ij}\dot{x}_{i}\dot{x}_{j}\geq 0\quad\mbox{for arbitrary $\dot{x}_{i}$},$$

$${\mu }_{ij}={\mu }_{ji}.$$

$$-\frac{\partial A}{\partial x_i}-\sum _j{\zeta }_{ij}\frac{d{x}_j}{dt}=0.$$

$$\mathcal{R}(\dot{x};x)=\frac{1}{2}\sum _{i,j}{\zeta }_{ij}{\dot{x}}_i{\dot{x}}_j+\sum _i\frac{\partial A}{\partial x_i}{\dot{x}}_i.$$

$$\Phi =\frac{1}{2}\sum _{i,j}{\zeta }_{ij}{\dot{x}}_i{\dot{x}}_j,$$

$$\dot{A}=\sum _i\frac{\partial A}{\partial x_j}{\dot{x}}_i.$$

$${\dot{x}}_i=\sum _p\frac{\partial x_i}{\partial {\alpha }_p}{\dot{\alpha }}_p,$$

$$\mathcal{R}=\sum _{pq}\left[\sum _{ij}{\zeta }_{ij}\frac{\partial x_i}{\partial {\alpha }_p}\frac{\partial x_j}{\partial {\alpha }_q}\right]{\dot{\alpha }}_p{\dot{\alpha }}_q+\sum _p\left[\sum _i\frac{\partial A}{\partial x_i}\frac{\partial x_i}{\partial {\alpha }_p}\right]{\dot{\alpha }}_p.$$

$$\tilde{\mathcal{R}}(\dot{x};x)=\mathcal{R}(\dot{x};x)-{\mathcal{R}}_{\min }(x),$$

$${\dot{x}}_i^{\ast }(x)=-\sum _j{\mu }_{ij}\frac{\partial A}{\partial x_j}.$$

$$\tilde{\mathcal{R}}(\dot{x};x)=\frac{1}{2}\sum _{ij}{\zeta }_{ij}({\dot{x}}_i-{\dot{x}}_i^{\ast })({\dot{x}}_j-{\dot{x}}_j^{\ast }).$$

$$\mathcal{O}[x(t)]=\int_{0}^{t}dt^{\prime}\,\tilde{\mathcal{R}}(\dot{x}(t^{\prime});x(t^{\prime}))=\frac{1}{2}\int_{0}^{t}dt^{\prime}\sum_{ij}\zeta_{ij}(x(t^{\prime}))\big{[}\dot{x}_{i}(t^{\prime})-\dot{x}_{i}^{*}(t^{\prime})\big{]}\big{[}\dot{x}_{j}(t^{\prime})-\dot{x}_{j}^{*}(t^{\prime})\big{]}.$$

$${\dot{x}}_i(t)=\frac{d{x}_i(t)}{dt},{\dot{x}}_i^{\ast }(t)=-{\sum _j({\zeta }^{-1}{)}_{ij}\frac{\partial A(x)}{\partial x_j}|}_{x=x(t)}.$$

$$\tilde{\Phi }=\frac{1}{2}\sum _{ab}{\tilde{\zeta }}_{ab}{\dot{y}}_a{\dot{y}}_b$$

$${\dot{x}}_i=\sum _{a=1}^{f^{\prime }}{c}_{ia}(x){\dot{y}}_a,$$

$${\mathcal{R}}^{(y)}(\dot{y};x)=\frac{1}{2}\sum _{a,b}{\tilde{\zeta }}_{ab}{\dot{y}}_a{\dot{y}}_b+\sum _{i,a}\frac{\partial A}{\partial x_i}{c}_{ia}{\dot{y}}_a.$$

$${\tilde{\mathcal{R}}}^{(y)}(\dot{y};x)=\frac{1}{2}\sum _{a,b}{\tilde{\zeta }}_{ab}({\dot{y}}_a-{\dot{y}}_a^{\ast })({\dot{y}}_b-{\dot{y}}_b^{\ast })$$

$$\mathcal{O}[y(t)]=\frac{1}{2}{\int }_0^{t}d{t}^{\prime }\sum _{a,b}{\tilde{\zeta }}_{ab}({\dot{y}}_a-{\dot{y}}_a^{\ast })({\dot{y}}_b-{\dot{y}}_b^{\ast })$$

$$\dot{n}=-\nabla \cdot \mathbf{j}$$

$$A=\int d\mathbf{r}f(n(\mathbf{r}))$$

$$\dot{A}=\int d\mathbf{r}\frac{\partial f}{\partial n}\dot{n}=-\int d\mathbf{r}\frac{\partial f}{\partial n}\nabla \cdot \mathbf{j}=\int d\mathbf{r}\mathbf{j}\cdot \nabla \left[\frac{\partial f}{\partial n}\right]=\int d\mathbf{r}(\mathbf{j}\cdot \nabla n)\frac{{\partial }^{2}f}{\partial n^2}$$

$$\Phi =\frac{1}{2}\int d\mathbf{r}n\zeta {\mathbf{v}}_p^2=\frac{1}{2}\int d\mathbf{r}\frac{1}{n}\zeta {\mathbf{j}}^2$$

$$\mathcal{R}=\frac{1}{2}\int d\mathbf{r}\frac{1}{n}\zeta {\mathbf{j}}^2+\int d\mathbf{r}(\mathbf{j}\cdot \nabla n)\frac{{\partial }^{2}f}{\partial n^2}$$

$${\mathbf{j}}^{\ast }=-\frac{n}{\zeta }\frac{{\partial }^{2}f}{\partial n^2}\nabla n.$$

$$\frac{\partial n}{\partial t}=\nabla \cdot \left[D(n)\nabla n\right]$$

$$\dot{n}=-\frac{\partial j}{\partial x}.$$

$$A=k_{B}T\int dxn(x)\ln n(x),$$

$$\mathcal{R}=\int dx\left[\frac{1}{2}\frac{1}{n}\zeta j^2+k_{B}T\frac{j}{n}\frac{\partial n}{\partial x}\right].$$

$$j^{\ast }=-D\frac{\partial n}{\partial x}$$

$$\frac{\partial n}{\partial t}=D\frac{{\partial }^{2}n}{\partial x^2}.$$

Full text

I Introduction

In 1931, Onsager published two papers on the dynamics of non-equilibrium systems Onsager (1931a, b). He had shown that if the rate of change of the system is written as a linear function of thermodynamic forces, the coefficients must be symmetric and positive definite. The symmetry of the coefficients is called reciprocal relation, and has become the base of classical irreversible thermodynamics de Groot and Mazur (1984); Kjelstrup and Bedeaux (2008); Öttinger (2005); Beris and Edwards (1994).

In the same papers, Onsager proposed a variational principle which is a direct consequence of the reciprocal relation, and states the essence of his theory (the kinetic equation plus the reciprocal relation) in a compact form. Many kinetic equations that have been used to describe the non-equilibrium phenomena, such as Stokes equation in hydrodynamics, Smoluchowskii equation in particle diffusion, Nernst-Planck equation in electro-kinetics, Ericksen-Leslie equation in nematic liquid crystals, and etc, can be derived from this principle Doi (2011, 2013). The variational principle has also been shown to be useful to derive time evolution equations for complex systems, where two or more irreversible processes are coupled, or where the process is constrained by geometry, say on curved surfaces or lines.

We have been showing that the variational principle is useful not only for deriving the time evolution equations, but also useful for obtaining the solutions of the kinetic equations Doi (2015); Xu et al. (2016); Zhou and Doi (2018). A way of doing this is to assume certain forms for the solutions which involve some time dependent parameters, and determine the parameters by the variational principle. We have demonstrated the utility of this method for many examples Di et al. (2016); Man and Doi (2016, 2017); Zhou and Doi (2018) in fluid and soft matter systems. The method was also applied to the study of solid-state dewetting Jiang et al. (2019).

On the other hand, our previous method has a weakness that it is only valid, strictly speaking, to predict the state in near future, i.e., to predict the state at time $t+\Delta t$ from the knowledge of the state at time $t$ for infinitesimally small $\Delta t$. Although we can use this method repeatedly to predict the state at times $t+2\Delta t$, $t+3\Delta t$, $\cdots$, the error may increase with time, and we cannot tell which solution is the best among all possible evolution paths.

Here we propose a new variational method which does not have this weakness. In the new method, we do not use the original variational principle proposed by Onsager in 1931, but use the integral proposed by Onsager and Machlup in 1953 Onsager and Machlup (1953). Compared with the previous method, the calculation of the new method is more cumbersome, but it has an advantage that we can directly determine which kinetic paths is the best among various possible kinetic paths. Also, it allows us to determine the steady state directly without solving the time evolution equations.

In this paper, we first explain the general framework of our new variational method (Sec. II). Next we demonstrate this method for a few examples, i.e., by solving simple diffusion equation (Sec. III), and by solving two hydrodynamic problems, the liquid wetting in a tube (Sec. IV), and a free boundary problem associated with coating (Sec. V).

II Onsager principle

II.1 Onsager principle in its original form

First we explain the original variational principle of Onsager Onsager (1931a, b). Since our target here is the flow and diffusion in soft matter, we shall limit our discussion to isothermal systems where temperature is assumed to be constant.

We consider a non-equilibrium system which is characterized by a set of state variables $x=(x_{1},x_{2},\cdots,x_{f})$. Our objective is to calculate the time dependence of $x(t)$ for a given initial state $x(t=0)$. (In this paper, we do not consider fluctuations, but fluctuations can be included in the present framework by considering the distribution function of $x$, see Ref. Doi (2011).) The time evolution of $x$ is obtained by solving Onsager’s kinetic equations

[TABLE]

where $A(x)$ is the free energy of the system and $\mu_{ij}(x)$ are kinetic coefficients, both are functions of $x$. Onsager showed that if the kinetic equation is written in the form of Eq. (1), $\mu_{ij}(x)$ must be positive definite and symmetric.

[TABLE]

Equation (3) is called Onsager’s reciprocal relation.

The kinetic equation (1) can be written in the form of force balance equation. Let $\zeta_{ij}$ be the inverse matrix of $\mu_{ij}$ ($\sum_{k}\zeta_{ik}\mu_{kj}=\delta_{ij}$). Equation (1) is then written as

[TABLE]

The first term represents the thermodynamic force that drives the system to the state of minimum free energy, and the second term represents the frictional force which resists against this change.

The force balance equation (4) can be cast in a form of minimum principle. Consider the following quadratic function of the rate of change of the state $\dot{x}=(\dot{x}_{1},\dot{x}_{2},\cdots,\dot{x}_{f})$:

[TABLE]

The kinetic equation (4) is equivalent to the condition that $\mathcal{R}$ is minimum with respect to $\dot{x}$, $\partial\mathcal{R}/\partial\dot{x}_{i}=0$, i.e., the time evolution of non-equilibrium system is determined by the minimum principle of $\mathcal{R}(\dot{x};x)$. This is called Onsager’s variational principle, or simply the Onsager principle.

The function $\mathcal{R}(\dot{x};x)$ is called Rayleighian. The first term is called the dissipation function

[TABLE]

and the second term is called the free energy change rate

[TABLE]

The dissipation function is a quadratic function of $\dot{x}$, while the free energy change rate is a linear function of $\dot{x}$.

II.2 Onsager principle as a tool of approximation

Onsager’s variational principle can be used to solve problems by a variational method. To do this, we assume that $x_{i}(t)$ ($i=1,2,...$) is written as a certain function of some parameter set $\alpha=(\alpha_{1},\alpha_{2},\cdots)$, i.e., $x_{i}(t)$ is written as $x_{i}(t)=x_{i}(\alpha_{1}(t),\alpha_{2}(t),\cdots)$. Then $\dot{x}_{i}$ is written as

[TABLE]

and the Rayleighian is written as

[TABLE]

The time evolution of the parameter $\alpha(t)$ is determined by minimizing $\mathcal{R}$ with respect to $\dot{\alpha}$. This method is useful when we have an idea for the kinetic path and can write down the functions $x_{i}(\alpha)$. It has been applied to many problems in soft matter Di et al. (2016); Man and Doi (2016); Xu et al. (2016); Man and Doi (2017); Zhou et al. (2017a, b); Zhou and Doi (2018).

II.3 Onsager principle in a modified form

The above variational principle determines the state at $t+\Delta t$ from the knowledge of the state at time $t$. It says that among all possible states allowed for the system to be at time $t+\Delta t$, the state chosen by nature is given by the state which minimizes $\mathcal{R}((x(t+\Delta t)-x(t))/\Delta t;x(t))$. This principle can be used to conduct approximate calculation. Suppose we have many candidates for the state at time $t+\Delta t$, then we can tell which state is the best candidate: the state which gives the smallest value of the Rayleighian is the best.

This variational principle is a local principle: it can predict the state in near future, but cannot predict the state in far future. Suppose we have two kinetic paths, starting from the same initial state $x(0)$ and ending at different states $x(t)$ and $x^{\prime}(t)$, then we cannot tell which path is better: we can tell which path is better if $t$ is close to zero, but cannot tell if $t$ is far away from zero.

This problem can be resolved if we use a slightly modified definition for Rayleighian

[TABLE]

where $\mathcal{R}_{\rm min}(x)$ is the minimum value of $\mathcal{R}(\dot{x};x)$ in the space of $\dot{x}$. The minimum is given by the velocity

[TABLE]

Note that $\dot{x}_{i}^{*}$ is the actual velocity of the system at state $x$, and depends only on the state variable $x$.

With the use of Eq. (11), it is easy to show that $\tilde{\mathcal{R}}(\dot{x};x)$ is written as

[TABLE]

Since $\zeta_{ij}$ is positive definite, $\tilde{\mathcal{R}}(\dot{x};x)$ is larger than or equal to zero.

Now consider the following integral, which is a functional of certain kinetic path $x(t)$

[TABLE]

The integral is positive definite and is equal to zero only when $x(t)$ is equal to the actual kinetic path. Hence the variational principle can be stated simply that nature chooses the path which minimizes the functional $\mathcal{O}[x(t)]$.

The integral of Eq. (13) was first introduced by Onsager and Machlup in their discussion on the fluctuation of kinetic paths described by linear Langevin equation Onsager and Machlup (1953). They have shown that the probability of finding the kinetic path $x(t)$ is proportional to \exp\left[-\mathcal{O}[x(t)]/{\color[rgb]{0,0,0}2}k_{B}T\right]. It is easy to show that their theory can be extended to the general case where both $A$ and $\zeta_{ij}$ are functions of $x$. Indeed in some literature Qian et al. (2006), this form has been referred to as Onsager’s variational principle. We shall call this variational principle Onsager Machlup principle, and call the integral of Eq. (13) Onsager Machlup integral.

II.4 Variational calculus using Onsager Machlup integral

The Onsager Machlup variational principle can be used to obtain the best guess for the kinetic path of the system. We consider certain kinetic path which involves a parameter set $\alpha=(\alpha_{1},\alpha_{2},\cdots)$. The best guess for the actual path is the path which gives the smallest value of the Onsager Machlup integral. In the following sections, we shall show actual calculation of this method. Here we show a few tips that will help such calculations.

The Rayleighian that appears in the Onsager Machlup integral has two velocities, $\dot{x}_{i}$ and $\dot{x}_{i}^{*}$. For given kinetic path $x(t)$, they are calculated by

[TABLE]

Here $\dot{x}_{i}$ is defined by the time derivative of the state variables $x(t)$, while $\dot{x}_{i}^{*}$ is defined by the thermodynamic force at state $x$. The Onsager Machlup integral becomes minimum when $\dot{x}_{i}$ is equal to $\dot{x}_{i}^{*}$.

At first sight, such variational principle may look useless since to calculate $\mathcal{O}[x(t)]$, we need to know $\dot{x}_{i}^{*}$, but this is the quantity difficult to calculate; if we can calculate \dot{x}_{{\color[rgb]{0,0,0}i}}^{*}, we can directly use it to solve the time evolution equation, $dx_{i}/dt=\dot{x}_{i}^{*}$ and there is no need to use such variational principle.

In fact, the variational principle is useful since we have a freedom to choose the kinetic path $x(t)$. We need not calculate $\dot{x}_{i}^{*}$ for a general state; we need to calculate $\dot{x}_{i}^{*}$ only for the state we have chosen. By proper choice of the kinetic path, $\dot{x}_{i}^{*}$ may be calculated.

For example, consider the problem of deformation of a droplet in emulsions under flow Doi (1983). It is difficult to calculate the shape change of the droplet for general shapes, but it can be calculated for special cases where the droplet takes simple shapes such as ellipsoids. Therefore, if the droplet keeps a form which can be approximated by ellipsoids, the variational principle can be used to describe the deformation of the droplet. We shall see this more clearly in the examples given later.

II.5 Note on the dissipation function

So far we have been discussing a simple case where the dissipation function is given explicitly as a function of $\dot{x}$, as in Eq. (6). In many situations, however, the dissipation function is not given in this form. Quite often, the dissipation function is written as a quadratic function of other velocity variables $\dot{y}_{a}$ $(a=1,2,...f^{\prime})$ in an extended space ($f^{\prime}>f$)

[TABLE]

subject to constraints that $\dot{x}_{i}$ and $\dot{y}_{a}$ are linearly related to each other,

[TABLE]

where the coefficients $c_{ia}$ are some functions of $x$. Example of such situation is given in the next section.

The dissipation function $\Phi$ in Eq. (6) is the minimum of $\tilde{\Phi}$ subject to the constraint of Eq. (16). Hence one can state the Onsager principle in the extended space spanned by $\dot{y}_{a}(a=1,2,...f^{\prime})$ Doi (2011), where the Rayleighian is given by

[TABLE]

It is easy to show that the modified Rayleighian in the extended space is given by

[TABLE]

where $\dot{y}_{a}^{*}$ is the velocity which minimizes $\mathcal{R}^{(y)}$. The Onsager Machlup integral is also given by

[TABLE]

Equations (18) and (19) includes two velocities, $\dot{y}_{a}$ and $\dot{y}_{a}^{*}$. The former is determined by the assumed kinetic path $x(\alpha(t))$, while the latter is determined by the thermodynamic forces acting at state $x(\alpha(t))$. If the kinetic equations are solved rigorously, they agree with each other, and $\tilde{\mathcal{R}}^{(y)}$ is equal to zero. In our approximate treatment, we determine the parameter $\alpha$ so that $\tilde{\mathcal{R}}^{(y)}$ becomes minimum. This is the general strategy that we will use in the subsequent examples.

III Diffusion

III.1 Diffusion equation

As the first example, we consider the diffusion of particles. The state variable in this problem is $n(\bm{r})$, the number density of particles at position $\bm{r}$, and we consider the time evolution of this function $n(\bm{r};t)$. The “velocity” of the state variable is $\dot{n}=\partial n/\partial t$, but it is not possible to write down the dissipation function in terms of $\dot{n}$. This is because in diffusion, $\dot{n}$ must satisfy the conservation equation

[TABLE]

where $\bm{j}(\bm{r})$ be the flux of the particles. Equation (20) represents the kinetic constraints that $n$ can change only by diffusion, and that there is no process (such as chemical reaction) which creates or annihilates the particles. The variable $\bm{j}$, introduced to represent such constraints, is an example of the variables $\dot{y}_{a}$ given in the previous section. The Rayleighian for the diffusion problem can be constructed using $\bm{j}$.

The free energy of the system is written as

[TABLE]

where $f(n)$ is the free energy density of the particles. $\dot{A}$ is then calculated by Eqs. (20) and (21) as

[TABLE]

The dissipation function $\Phi$ is proportional to the square of the particle velocity $\bm{v}_{p}=\bm{j}/n$, and can be written as

[TABLE]

where $\zeta$ is the friction constant of one particle.

From Eqs. (22) and (23), the Rayleighian is given by

[TABLE]

The minimum of this is given by

[TABLE]

Therefore the time evolution equation for $n$ is given by

[TABLE]

where $D=(n/\zeta)\partial^{2}f/\partial n^{2}$. Equation (26) is the standard diffusion equation Doi (2013).

III.2 Onsager principle as an approximation tool

We now demonstrate how to use Onsager principle to obtain approximate solutions. For simplicity, we consider the diffusion in 1-dimension system and in dilute solution. The conservation equation (20) is written as

[TABLE]

The free energy of the system is given by

[TABLE]

and the Rayleighian is given by [see Eq. (24)]

[TABLE]

This becomes minimum when $j$ is equal to

[TABLE]

where $D=k_{B}T/\zeta$. Equation (30) and the conservation equation (27) give the 1-dimension diffusion equation

[TABLE]

The exact solution of this equation for the initial condition $n(x,0)=N_{0}\delta(x)$ is given by

[TABLE]

To use the variational principle, we assume that $n(x,t)$ is approximated by the following polynomial function

[TABLE]

where $m$ is a positive integer, $a(t)$ is a parameter characterizing the particle spreading, and $N$ is a constant which has been determined by the condition $\int_{-\infty}^{\infty}dxn(x,t)=N_{0}$.

We shall derive the kinetic equations for $a(t)$ using Onsager principle. From Eqs. (28) and (33), the free energy is calculated as

[TABLE]

where $z=|x|/a$. Hence the free energy change rate is calculated as

[TABLE]

On the other hand, the flux $j(x)$ is calculated by integrating Eq. (27) with respect to $x$,

[TABLE]

The dissipation function is then calculated as

[TABLE]

The Rayleighian is the summation of Eqs. (35) and (37)

[TABLE]

Minimizing this with respect to $\dot{a}$, we have

[TABLE]

The solution of this equation for the initial condition $a(0)=0$ is

[TABLE]

It is seen that the spreading of the particle concentration follows the scaling $a\sim\sqrt{t}$ for all $m$.

Figure 1 shows the comparison of the approximate solution $n(x,t)$ calculated by Eqs. (33) and (39) and the exact solution (32). From this figure, we can see that the curve of $m=1$ seems to be the best solution (the solution closest to the exact solution). However, if we do not know the exact solution, we cannot tell which solution is the best. This question can be answered if we use the Onsager Machlup integral.

III.3 Onsager Machlup integral

We now analyze the problem using the modified Rayleighian. In the present problem, the modified Rayleighian is obtained by Eqs. (18) and (23) as

[TABLE]

where $j$ is the flux determined by the assumed time dependence of $n(x;t)$, i.e., the flux given by Eq. (36)

[TABLE]

and $j^{*}$ is the exact flux at state $n(x,t)$, which is given by

[TABLE]

The modified Rayleighian is then calculated as

[TABLE]

The evolution is obtained by $\partial\tilde{\mathcal{R}}/\partial\dot{a}=0$, which leads to the same equation as Eq. (39).

Substituting Eq. (39) into the modified Rayleighian (44)

[TABLE]

The integration of the last term diverges as $z\rightarrow 1$. We performed the integration by setting the upper bound to $1-\epsilon$. For $m=1,2,3$, the results are

[TABLE]

If we compare the front factor of the divergent terms, the trial function $m=1$ gives the smallest value. Also for a typical small value of $\epsilon=10^{-3}$, the terms in the brackets are evaluated to be 0.3233, 0.5902, 0.6695 for $m=1,2,3$, respectively. In both cases, the trial function of $m=1$ gives the smallest value, and we can conclude that the function of $m=1$ is the best. This is consistent with the comparison shown in Fig. 1.

IV Finger flow in a square tube

As the second example, we consider the liquid wetting in a square tube (see Fig. 2). The liquid is contained in the left reservoir and connected to a tube having a square cross section. The end of the tube is closed, so the meniscus of the bulk liquid cannot move, but the liquid can advance along the corner of the tube forming a “finger” to reduce the surface energy. The tube is assumed to be placed horizontally, so gravity does not play a role in this problem.

The wetting dynamics of such situation was first studied by Dong and Chatzis Dong and Chatzis (1995) using a variational method Mayer et al. (1983). (See also Ref. Wei et al. (2011) for the application.) Here we use the same model as theirs, and solve the problem using Onsager principle. The governing equation of this problem turns out to be the same as the non-linear diffusion equation, and we shall show that the equation can be solved in good approximation by using the Onsager Machlup integral.

IV.1 Dynamical equation

We take $x$ coordinate along the tube axis. We define a dimensionless quantity called saturation $s(x;t)$ , which is the fraction of the liquid area in the cross section at position $x$ and time $t$. We focus on the motion of the finger part, and ignore the motion of the bulk part. We set the bottom of the finger at $x=0$, and assume, as in the previous works Dong and Chatzis (1995); Yu et al. (2018), that $s(x,t)$ takes a constant value $s^{*}$ at $x=0$

The time evolution equation for $s(x,t)$ can be derived from the Onsager principle. Since the gravity is ignorable, the free energy of the system is given by the surface energy of the liquid. The surface area of the liquid in the region between $x$ and $x+dx$ is proportional to $a\sqrt{s}dx$. Hence the free energy of the system can be written as

[TABLE]

where $\alpha_{1}$ is a certain numerical constant ($\alpha_{1}=8\sqrt{1-\pi/4}$, see Ref. Yu et al. (2018)).

The saturation $s(x;t)$ must satisfy the conservation equation

[TABLE]

The dissipation function is a quadratic function of $j(x)$. By using the lubrication approximation, one can show that the dissipation function is written as

[TABLE]

where $\eta$ is the liquid viscosity and $\alpha_{2}$ is another numerical constant. The hydrodynamic calculation Ransohoff and Radke (1988); Yu et al. (2018) indicates that $\alpha_{2}\simeq 90$.

The Rayleighian is given by Eq. (46) and Eq. (48), and the same procedure as in the previous section gives the following expression for the fluid flux

[TABLE]

where $D$ is defined by

[TABLE]

which has the same dimension as the diffusion constant.

The time evolution equation is given by

[TABLE]

This is the same equation as that derived by Dong et al.Dong and Chatzis (1995).

The boundary conditions and the initial conditions for $s(x,t)$ are

[TABLE]

Dong et al. solved Eq. (51) under these conditions, and showed that the liquid front $h(t)$ (see Fig. 2) increases with time as

[TABLE]

In the following we shall obtain $h(t)$ using Onsager Machlup integral.

IV.2 Onsager Machlup integral

We assume the profile can be written in the form

[TABLE]

This function is chosen to satisfy the boundary condition $s(0,t)=s^{*}$ and the condition $s(h(t);t)=0$. The function includes two parameters, $h(t)$ and $n$, which we shall determine using the Onsager principle.

The flux at position $x$ is given by integrating Eq. (47) and using the condition that $j$ becomes 0 at $x=h(t)$

[TABLE]

where $z=x/h(t)$, and the function $F(z)$ is given by

[TABLE]

On the other hand, the flux $j^{*}$ is given by Eq. (49)

[TABLE]

The modified Rayleighian is given by

[TABLE]

Using Eqs. (56) and (58), we can get

[TABLE]

The minimum of this equation is given by

[TABLE]

with

[TABLE]

The tip of the finger advances with the Lucas-Washburn scaling $h\sim t^{1/2}$. Different profile predicts a different front factor $H_{n}$ for the tip position.

Figure 3 shows the comparison of the approximate solution $s(x,t)$ calculated by Eq. (61) and the numerical solutions to Eq. (51). We see that the saturation profile of $n=2$ is the closest to the numerical solution [Fig. 3(a)]. We can also compare the front factor $H_{n}$ in the tip position of the finger. They are given by Eq. (62) and shown in Fig. 3(b). Numerical solution to the partial differential equation (51) with boundary condition (52) and initial condition (53) gives $H\simeq 0.1281$ Dong and Chatzis (1995). Again, the graph indicates that the best estimate for $h(t)$ is obtained for $n=2$.

We can confirm this observation by calculating the minimum value of the modified Rayleighian (60)

[TABLE]

with

[TABLE]

Fig. 4 shows the value of $C_{n}$, which is proportional to the Onsager Machlup integral, as a function of $n$. It is seen that $C_{n}$ takes minimum at integer number $n=2$, in agreement with the results shown in Fig. 3. If we relax and allow $n$ to be a real number, we can see that $C_{n}$ has a minimum around $n\simeq 1.68$. The saturation profile for $n=1.68$ is also shown in Fig. 3(a) in black, which resembles the numerical solution even better than $n=2$. This also demonstrates that the best approximation can be obtained by evaluating the Onsager Machlup integral.

V Film Coating

As the last example, we consider the problem of liquid coating on a solid substrate (see Fig. 5). A substrate is moving out from a reservoir of viscous liquid with velocity $U$ through a gap $h_{0}$, and is coated with the liquid. The liquid in the reservoir is kept at a pressure $P_{\rm in}$ higher than the atmospheric pressure. (Here we consider the situation that $P_{\rm in}$ is negative, so the liquid is sucked in to the reservoir.) We also assume that the gap $h_{0}$ is much smaller than the capillary length and ignore the effect of gravity.

We focus on the steady state of the system, and ask what is the thickness $h_{f}$ of the liquid film formed on the substrate in the steady state. This question is slightly different from the previous ones. In the previous examples, we discussed the time evolution of non-equilibrium systems, but here we discuss the steady state of non-equilibrium systems. The purpose of this example is to show that we can calculate the steady state directly, without solving the time evolution equation, using Onsager Machlup integral.

V.1 Time evolution equation for the film profile

First we consider the time evolution of the film profile. We take $x$ coordinate parallel to the substrate, the origin of which is located at the exit of the reservoir. Let $h(x,t)$ be the thickness of the liquid film at point $x$ and time $t$. The free energy of the system is given by the sum of the surface energy and the potential energy due the reservoir pressure $P_{\rm in}$

[TABLE]

where $\gamma$ is the surface tension, and we have assumed $(\partial h/\partial x)^{2}\ll 1$.

Let $v(x,t)$ be the depth averaged velocity of the fluid. The conservation condition of the fluid is written as

[TABLE]

$h(x,t)$ should satisfy the boundary condition

[TABLE]

The free energy change rate $\dot{A}$ is calculated from Eqs. (65) and (66)

[TABLE]

where we have used that $\dot{h}$ is equal to 0 at $x=0$. On the other hand, the energy dissipation function is calculated by the lubrication approximation Di et al. (2018) as

[TABLE]

where $\eta$ is the viscosity of the fluid. The velocity which minimizes the Rayleighian $\Phi+\dot{A}$ is given by

[TABLE]

The minimization of the Rayleighian also gives the following boundary condition:

[TABLE]

Substituting the velocity (70) into the conservation law (66), we have the evolution equation for the film thickness

[TABLE]

where $\text{Ca}=U/U^{*}$ is the capillary number and $U^{*}=\gamma/\eta$ is the capillary velocity.

V.2 Variational calculus for steady state using Onsager Machlup integral

We now use the variational principle to seek an approximate solution for the steady state of Eq. (72). The variational principle says that among all possible kinetic paths allowed for $h(x;t)$, nature chooses the path which minimizes the Onsager Machlup integral (19).

In the present problem, the the modified Rayleighian is given by

[TABLE]

where we have again used the formula of Eq. (18). In Eq. (73), $v$ is the velocity determined by the conservation equation (66), and $v^{*}$ is the velocity given by Eq. (70).

In the steady state, $v$ and $v^{*}$ become independent of time. Therefore the minimization principle of the Onsager Machlup integral is equivalent to the minimization of the modified Rayleighian $\tilde{\mathcal{R}}$.

To seek the minimum of Eq. (73), we take the same strategy as in the previous examples; we consider certain form for $h(x)$ which includes certain parameters, and seek the minimum in the parameter space. A simple choice for the liquid profile $h(x)$ is

[TABLE]

This function is chosen to satisfy the condition $h(0)=h_{0}$ and $h(\infty)=h_{f}$. Equation (74) includes 4 parameters $a,b,x_{f}$ and $h_{f}$. The number of parameters is reduced to 1 since $h(x)$ has to satisfy the boundary condition (71), and the continuity condition for $h(x)$ and $h^{\prime}(x)$ at $x=x_{f}$. However, Eq. (74) has a problem that the second order differential $h^{\prime\prime}(x)$ is discontinuous at $x=x_{f}$, and therefore $v^{*}$ has a singularity of delta function type $\delta(x-x_{f})$ which makes the modified Rayleighian $\tilde{\mathcal{R}}[h(x)]$ diverge. Therefore, we considered the following form

[TABLE]

This function includes 6 parameters, $a,b,x_{1},h_{1},h_{f}$ and ,$\kappa$, but the number of independent parameters can be reduced to 2 if we use the boundary condition (71), and the 3 continuity conditions for $h$, $h^{\prime}$ and $h^{\prime\prime}$ at $x=x_{1}$. We have chosen $h_{1}$ and $h_{f}$ as independent parameters. The other parameters are expressed by $h_{1}$ and $h_{f}$ as follows.

[TABLE]

At steady state, the flux $hv$ is constant and is given by $h_{f}U$. Hence the velocity $v(x)$ is given by

[TABLE]

On the other hand, the velocity $v^{*}(x,t)$ is given by Eq. (70).

Using Eqs. (76) and (70), we calculated the modified Rayleighian for the profile (75), and obtained the following result

[TABLE]

Using dimensionless parameters $\tilde{h}_{0}=ah_{0}/\text{Ca}^{2/3}$, $\tilde{h}_{1}=ah_{1}/\text{Ca}^{2/3}$, $\tilde{h}_{f}=ah_{f}/\text{Ca}^{2/3}$, the above equation can be written as

[TABLE]

Minimization of Eq. (78) with respect to $\tilde{h}_{1}$ and $\tilde{h}_{f}$ gives the steady state film thickness $h_{f}$. Equation (78) indicates that $\tilde{h}_{f}$ depends on $\tilde{h}_{0}$, but does not depend on Ca. Figure 6 shows the contour maps of $\tilde{\mathcal{R}}$ for various values of $\tilde{h}_{0}$, $\tilde{h}_{0}=10,100,1000$. The $x$-axis represents $\tilde{h}_{f}$ and the $y$-axis represents $\tilde{h}_{1}-\tilde{h}_{f}$.

Figure 6 shows that the minimum position is almost independent of of $\tilde{h}_{0}$, indicating that the film thickness is given by $\tilde{h}_{f}(\approx 1.07)$ and is almost independent of $\tilde{h}_{0}$， which gives us the celebrated scaling law

[TABLE]

This can be compared with result of the asymptotic solution for the steady state for Eq. (72) Bretherton (1961); Carvalho and Kheshgi (2000), which gives the same scaling law as Eq. (79), but the numerical coefficient is different, equal to 1.34. Figure 7 shows the comparison between the asymptotic analysis and the variational calculus. Our calculation indicates that the effect of gap distance $h_{0}$ is small.

VI Conclusion

In this paper, we have shown a new usage of Onsager Machlup variational principle. It is based on the minimum principle that nature chooses the kinetic path which makes the Onsager Machlup integral $\mathcal{O}[x(t)]$ minimum. We have demonstrated that this principle can be used to get approximate solutions for the evolution equations for linear and non-linear diffusion equations. These approximate solutions in general involve one or a few time-dependent parameters, such as $a(t)$ in the diffusion problem (33) and $h(t)$ in the finger flow problem (55). These parameters are optimized step by step in time, determined by the kinetic equations derived from Onsager principle. On the other hand, there are also one or a few parameters that remain constant in time, such as the power $m$ in Eq. (33) and $n$ in Eq. (55). These time-independent parameters represent the functional form of the approximate solutions. They cannot be optimized locally in time, yet they can still be optimized globally over a period of time through the minimization of Onsager Machlup integral. We have shown that the Onsager Machlup integral can be used to obtain the steady state without solving the evolution equation in our coating example.

The present method may be regarded as a special kind of least square minimization method, but there is a significant difference. In the existing least square minimization methods, the choice of the evaluation function (the function which tells the goodness of the solution) is arbitrary. If the evolutions equations is given by a set of equations, one can put any weight for each equation to get the evaluation function, and the final result depends on the choice of the weight. What we are proposing in this paper is that the Onsager Machlup integral is the most natural choice for the evaluation function since it represents the fundamental quantity which governs the physics in the problem.

Many extensions are possible based on the Onsager Machlup integral. The method may be used to find new numerical scheme to solve the kinetic equations. It can also be used to obtain oscillatory solutions in space and time for non-linear partial differential equations. More examples will be shown in future which demonstrate the powerfullness of Onsager’s variational principle.

Acknowledgements.

M.D. thanks Prof. Hans Christian Öttinger of ETH for his illuminating discussion on the Onsager principle which initiated this work. This work was supported by the National Natural Science Foundation of China (NSFC) through the Grant No. 21504004 and 21774004. M.D. acknowledges the financial support of the Chinese Central Government in the Thousand Talents Program.