Derivation of Two-Fluid Model Based on Onsager Principle

Abstract

Using the Onsager variational principle, we study the dynamic coupling between the stress and the composition in a polymer solution. In the original derivation of the two-fluid model of Doi and Onuki the polymer stress was introduced a priori; therefore, a constitutive equation is required to close the equations. Based on our previous study of viscoelastic fluids with homogeneous composition, we start with a dumbbell model for the polymer, and derive all dynamic equations using the Onsager variational principle.

1. Introduction

In the studies of flow for polymeric liquids [1,2], the inclusion of polymers introduces two new variables in the system, namely the polymer concentration and the polymer conformation, which are absent in the flow of simple liquids. It is well recognized that the time evolution of the microscopic state variable, i.e., the local conformation of the polymer chain, is critically important in governing the dynamics of the polymer solution [3,4]. The total stress of a polymer solution therefore has two contributions, one from the polymer and another from the solvent. A prescribed constitutive equation is required to relate the polymer stress to the local flow conditions. On the other hand, in the standard treatment for the flow of polymeric liquids [1], the polymer concentration is assumed to be uniform in space. Therefore, the polymer concentration appears as a parameter in the framework, and there is no time evolution equation for the polymer concentration.

Experimentally, it has been shown that the polymer concentration can become non-uniform when the velocity gradient is not uniform [5,6,7]. Theoretically, a phenomenological two-fluid model has been developed [8,9,10], which incorporates the coupling between polymer stress and polymer diffusion in the continuum framework. A simple Hookean dumbbell model is used for the polymer chain, and there are a few studies based on different strategies [11,12,13,14].

In this manuscript, we shall re-derive the two-fluid model based on the Onsager principle [15]. In Section 2, we present a general derivation including all viscous coupling in the dissipation. In Section 3, we repeat Doi–Onuki’s original derivation from Ref. [10] for reference. In Section 4, we start with a dumbbell model for the polymer, and derive the time evolution equations using the Onsager variational principle. We conclude with a summary in Section 5.

2. Onsager Principle

First proposed by Onsager in his celebrated papers on the reciprocal relation [16,17], the Onsager principle is a variational principle to systematically derive the time evolution equations for out-of-equilibrium systems. The first step is to identify a set of state variables, $x=(x_1,x_2,\cdots )$, which characterizes the non-equilibrium state of the system under study. Then, the time evolution of the system is determined by the condition that the following quadratic function of $\dot{x}=({\dot{x}}_1,{\dot{x}}_2,\cdots )$ to be minimized with respect to $\dot{x}$,

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

(1)

Here, we use the dot for the partial time derivative, $\dot{x}=\partial x/\partial t$.

Equation (1) defines the Rayleighian of the system. It consists of two parts: one is the time derivative of the free energy $\dot{A}\left(x\right)={\sum }_i(\partial A/\partial x_i){\dot{x}}_i$. The other part is called the dissipation function $\mathsf{\Phi }=\frac{1}{2}{\sum }_{i,j}{\zeta }_{ij}{\dot{x}}_i{\dot{x}}_j$, where $2\mathsf{\Phi }$ represents the energy dissipated in the system per unit time when the state variables are changing at rate $\dot{x}$. The coefficient ${\zeta }_{ij}$ is called the friction coefficient, which is generally a function of state variables $x_i$. The dissipation function must be a quadratic function of ${\dot{x}}_i$. The minimum condition of the Rayleighian $\partial \mathcal{R}/\partial {\dot{x}}_i=0$ determines the time evolution of the state variables:

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

(2)

Equation (2) is an analogue to the force balance equation, where the left-hand side is the thermodynamic driving force and the right-hand side is the friction force. The reciprocal relation ${\zeta }_{ij}={\zeta }_{ji}$ is required in this derivation.

The Onsager principle is particularly useful for soft matter systems when inertia is not important. Many time evolution equations used in soft matter, such as the Stokes equation, Fick’s diffusion equation, Nernst–Planck equation, Cahn–Hilliard equation, Ericksen–Leslie equation, etc., can be derived based on the Onsager principle [15,18,19]. In a previous work [20], we have shown that the continuum mechanical equation for viscoelastic fluids can also be derived from the Onsager principle. Here, we use the same framework to derive the time evolution equations of a two-fluid model for polymer solutions, by taking into consideration the coupling between stress and diffusion.

2.1. State Variables

We first need to identify the state variables that characterize the non-equilibrium state of flowing polymer solutions. We choose the state variables as follows:

• Volume fraction of the polymer $\varphi$.
  The corresponding “velocity” variable is ${\mathbf{v}}^{\left(p\right)}$, the polymer velocity. The polymer volume fraction and the polymer velocity are related by the conservation law

  $$\dot{\varphi }=-\mathsf{\nabla }·\left(\varphi {\mathbf{v}}^{\left(p\right)}\right)=-{\nabla }_{\alpha }\left(\varphi v_{\alpha }^{\left(p\right)}\right).$$

  (3)
  Here, $\dot{\varphi }=\partial \varphi /\partial t$ and ${\nabla }_{\alpha }=\partial /\partial x_{\alpha }$.
• Conformation tensor $\mathit{c}$.
  $\mathit{c}$ is a non-dimensional tensor to characterize the microscopic state of the polymer chain. The $\mathit{c}$-tensor is equal to unit tensor $\mathit{I}$ when the polymer is at equilibrium, and deviates from $\mathit{I}$ when the polymer is deformed. Later, we will introduce the dumbbell model, which presents the polymer chain as a dumbbell consisting of two beads at positions ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$. These two beads are connected by an elastic spring that has an end-to-end vector $\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2$ and a spring constant k. The conformation of the dumbbell is then specified by the $\mathit{c}$-tensor defined by $\mathit{c}=\frac{k}{k_{B}T}\langle \mathbf{r}\mathbf{r}\rangle$.
  The corresponding “velocity” variable is the material time derivative of $\mathit{c}$ defined by

  $$D_t^{\left(p\right)}\mathit{c}=\frac{\partial \mathit{c}}{\partial t}+{\mathbf{v}}^{\left(p\right)}·\mathsf{\nabla }\mathit{c},\mathrm{or}{D}_t^{\left(p\right)}{c}_{\alpha \beta }=\frac{\partial c_{\alpha \beta }}{\partial t}+v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }{c}_{\alpha \beta }.$$

  (4)
  Notice that, here, we are using the polymer velocity ${\mathbf{v}}^{\left(p\right)}$ to define the material time derivative, not the medium velocity $\mathbf{v}$, which will be introduced next.
• In order to discuss the phenomena of the diffusion or migration of polymers, we need to introduce another “velocity” variable representing the velocity of the surroundings. This can be represented by the solvent velocity ${\mathbf{v}}^{\left(s\right)}$, or the medium velocity (volume-average velocity) defined by

  $$\mathbf{v}=\varphi {\mathbf{v}}^{\left(p\right)}+(1-\varphi ){\mathbf{v}}^{\left(s\right)}.$$

  (5)
  Here, we will use $\mathbf{v}$ following rheology convention.

The local flow condition is characterized by the velocity gradient tensor $\mathsf{\nabla }\mathbf{v}$

$${\nabla }_{\alpha }{v}_{\beta }=\frac{\partial v_{\beta }}{\partial x_{\alpha }},$$

(6)

and the related rate-of-strain tensor

$$\dot{\mathbf{\gamma }}={(\mathsf{\nabla }\mathbf{v})}^t+\mathsf{\nabla }\mathbf{v}=\left(\frac{\partial v_{\alpha }}{\partial x_{\beta }}+\frac{\partial v_{\beta }}{\partial x_{\alpha }}\right).$$

(7)

2.2. Free Energy

The general form of the free energy of a polymer solution can be written as

$$A=\int \mathrm{d}\mathbf{r}a(\varphi ,\mathit{c}).$$

(8)

where $a(\varphi ,\mathit{c})$ is the free energy density. We assume that $a(\varphi ,\mathit{c})$ has the following form

$$a(\varphi ,\mathit{c})=f\left(\varphi \right)+\varphi g\left(\mathit{c}\right).$$

(9)

The first term $f\left(\varphi \right)$ is the free energy density of polymer solutions at equilibrium. This term includes the entropic term $\varphi ln\varphi$, the interaction term of Flory–Huggins form $\chi \varphi (1-\varphi )$, and the interfacial energy, which depends on the concentration gradient ${|\nabla \varphi |}^2$. The second term includes $g\left(\mathit{c}\right)$, which represents the elastic energy of deformed polymer chains and is a function of the conformation tensor $\mathit{c}$. The second term is proportional to the polymer volume fraction $\varphi$.

The change rate of the free energy is given by

$$\begin{array}{ccc}\hfill \dot{A}& =& \int \mathrm{d}\mathbf{r}\left[\frac{\partial a}{\partial \varphi }\dot{\varphi }+\frac{\partial a}{\partial c_{\alpha \beta }}{\dot{c}}_{\alpha \beta }\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[-\frac{\partial a}{\partial \varphi }{\nabla }_{\gamma }\left(\varphi v_{\gamma }^{\left(p\right)}\right)+\frac{\partial a}{\partial c_{\alpha \beta }}{\dot{c}}_{\alpha \beta }\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[\varphi v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\left(\frac{\partial a}{\partial \varphi }\right)+\frac{\partial a}{\partial c_{\alpha \beta }}\left(D_t^{\left(p\right)}{c}_{\alpha \beta }-v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }{c}_{\alpha \beta }\right)\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\left(\varphi \frac{\partial a}{\partial \varphi }\right)-v_{\gamma }^{\left(p\right)}\frac{\partial a}{\partial \varphi }{\nabla }_{\gamma }\varphi -v_{\gamma }^{\left(p\right)}\frac{\partial a}{\partial c_{\alpha \beta }}{\nabla }_{\gamma }{c}_{\alpha \beta }+\frac{\partial a}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\left(\varphi \frac{\partial a}{\partial \varphi }\right)-v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }a+\frac{\partial a}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\left(\varphi \frac{\partial a}{\partial \varphi }-a\right)+\frac{\partial a}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right]\hfill \\ & =& \int \mathrm{d}\mathbf{r}\left[v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\Pi +\frac{\partial a}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right].\hfill \end{array}$$

(10)

where $\Pi$ stands for the osmotic pressure and is defined by

$$\Pi \equiv \varphi \frac{\partial a}{\partial \varphi }-a.$$

(11)

Notice that the elastic energy term $g\left(\mathit{c}\right)$ has no contribution to the osmotic pressure

$$\Pi =\varphi \frac{\partial (f+\varphi g)}{\partial \varphi }-(f+\varphi g)=\varphi \frac{\partial f}{\partial \varphi }-f.$$

(12)

2.3. Dissipation Function

The dissipation function also includes several terms. The first one accounts for the relative motion of the center-of-mass of the polymer with respect to the medium velocity

$${\mathsf{\Phi }}_{pv}=\frac{1}{2}\int \mathrm{d}\mathbf{r}\xi {({\mathbf{v}}^{\left(p\right)}-\mathbf{v})}^2,$$

(13)

where the friction coefficient $\xi =\xi \left(\varphi \right)$ is, in general, concentration-dependent. Notice that, by use of Equation (5), the integrand of Equation (13) is written as $\xi {(1-\varphi )}^2{({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})}^2$. Therefore, ${\mathsf{\Phi }}_{pv}$ can be understood as the dissipation due to the relative motion between polymer and solvent.

The second term represents the coupling between $D_t^{\left(p\right)}\mathit{c}$ and $\mathsf{\nabla }\mathbf{v}$. We write it in a very general form:

$${\mathsf{\Phi }}_{cv}=\frac{1}{2}\int \mathrm{d}\mathbf{r}\left\{{\xi }_{\alpha \beta \mu \nu }^{\left(cc\right)}(D_t^{\left(p\right)}{c}_{\alpha \beta })(D_t^{\left(p\right)}{c}_{\mu \nu })+2{\xi }_{\alpha \beta \mu \nu }^{\left(cv\right)}(D_t^{\left(p\right)}{c}_{\alpha \beta })\left({\nabla }_{\mu }{v}_{\nu }\right)+{\xi }_{\alpha \beta \mu \nu }^{\left(vv\right)}\left({\nabla }_{\beta }{v}_{\alpha }\right)\left({\nabla }_{\mu }{v}_{\nu }\right)\right\}$$

(14)

The first term is the inter-coupling of $D_t^{\left(p\right)}\mathit{c}$, the second term is the cross-coupling between $D_t^{\left(p\right)}\mathit{c}$ and $\mathsf{\nabla }\mathbf{v}$, and the last term is related to the solvent viscosity. Since $D_t^{\left(p\right)}\mathit{c}$ and $\mathsf{\nabla }\mathbf{v}$ are tensors of rank 2, the frictional coefficients $\zeta$ are tensors of rank 4. These coefficients must be positive-definite to ensure that the dissipation function is non-negative.

2.4. Time Evolution Equations

From the change rate of the free energy (11) and the dissipation functions (13) and (14), the Rayleighian can be written as

$$\begin{array}{ccc}\hfill \mathcal{R}& =& \int \mathrm{d}\mathbf{r}\left[v_{\alpha }^{\left(p\right)}{\nabla }_{\alpha }\Pi +\frac{\partial a}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right]+\frac{1}{2}\int \mathrm{d}\mathbf{r}\xi {(v_{\alpha }^{\left(p\right)}-v_{\alpha })}^2\hfill \\ & & +\frac{1}{2}\int \mathrm{d}\mathbf{r}\left\{{\xi }_{\alpha \beta \mu \nu }^{\left(cc\right)}(D_t^{\left(p\right)}{c}_{\alpha \beta })(D_t^{\left(p\right)}{c}_{\mu \nu })+2{\xi }_{\alpha \beta \mu \nu }^{\left(cv\right)}(D_t^{\left(p\right)}{c}_{\alpha \beta })\left({\nabla }_{\mu }{v}_{\nu }\right)+{\xi }_{\alpha \beta \mu \nu }^{\left(vv\right)}\left({\nabla }_{\beta }{v}_{\alpha }\right)\left({\nabla }_{\mu }{v}_{\nu }\right)\right\}\hfill \\ \hfill & & -\int \mathrm{d}\mathbf{r}p({\nabla }_{\alpha }{v}_{\alpha }),\hfill \end{array}$$

(15)

where the last term accounts for the incompressibility condition $\mathsf{\nabla }·\mathbf{v}=0$.

By variational calculation with respect to the three “velocity” variables, we obtain the following set of equations:

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta v_{\alpha }^{\left(p\right)}}=0& \Rightarrow & \xi (v_{\alpha }^{\left(p\right)}-v_{\alpha })+{\nabla }_{\alpha }\Pi =0\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta v_{\alpha }}=0& \Rightarrow & -\xi (v_{\alpha }^{\left(p\right)}-v_{\alpha })-{\nabla }_{\beta }\left({\xi }_{\mu \nu \alpha \beta }^{\left(cv\right)}(D_t^{\left(p\right)}{c}_{\mu \nu })\right)-{\nabla }_{\beta }\left({\xi }_{\mu \nu \alpha \beta }^{\left(vv\right)}\left({\nabla }_{\mu }{v}_{\nu }\right)\right)+{\nabla }_{\alpha }p=0\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta D_t^{\left(p\right)}{c}_{\alpha \beta }}=0& \Rightarrow & {\xi }_{\alpha \beta \mu \nu }^{\left(cc\right)}(D_t^{\left(p\right)}{c}_{\mu \nu })+{\xi }_{\alpha \beta \mu \nu }^{\left(cv\right)}\left({\nabla }_{\mu }{v}_{\nu }\right)+\frac{\partial a}{\partial c_{\alpha \beta }}=0\hfill \end{array}$$

(18)

Equation (16) and the conservation Equation (3) lead to the time evolution equation for $\varphi$. Equations (16) and (17) give the following force balance equation:

$${\nabla }_{\beta }\left(-{\xi }_{\mu \nu \alpha \beta }^{\left(cv\right)}(D_t^{\left(p\right)}{c}_{\mu \nu })-{\xi }_{\mu \nu \alpha \beta }^{\left(vv\right)}\left({\nabla }_{\mu }{v}_{\nu }\right)+(\Pi +p){\delta }_{\alpha \beta }\right)=0,$$

(19)

from which we can identify the stress tensor as the terms in the large brackets. Equation (18) gives the constitutive equation

$$D_t^{\left(p\right)}{c}_{\mu \nu }=-\left({\nabla }_{\gamma }{v}_{ϵ}\right){\xi }_{\alpha \beta \gamma ϵ}^{\left(cv\right)}{\xi }_{\alpha \beta \mu \nu }^{{\left(cc\right)}^{-1}}-\frac{\partial a}{\partial c_{\alpha \beta }}{\xi }_{\alpha \beta \mu \nu }^{{\left(cc\right)}^{-1}}.$$

(20)

These results are quite general, but their usefulness is limited because many phenomenological parameters are introduced in the model. It is not clear how to assign values to these parameters for a practical polymer model. In the following, we shall study two existing theories from the viewpoint of this general formulation.

3. Doi–Onuki Derivation

3.1. Polymer Solution

Doi and Onuki [10] and Mavrantzas and Beris [12] proposed a two-fluid model for polymer solutions. Their theories are not equivalent to each other, but are based on the same idea that the elastic stress created in the polymer will contribute to the motion of the polymer relative to the solvent if the stress is not uniform. Here, we will discuss the two-fluid model focusing on the Doi–Onuki theory. They started with a free energy change rate of the form

$$\dot{F}={\dot{F}}_{\mathrm{mix}}+{\dot{F}}_{\mathrm{el}}.$$

(21)

The first term ${\dot{F}}_{\mathrm{mix}}$ is the change rate of the mixing free energy

$$\begin{array}{ccc}\hfill F_{\mathrm{mix}}& =& \int \mathrm{d}\mathbf{r}f(\varphi ),\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\dot{F}}_{\mathrm{mix}}& =& \int \mathrm{d}\mathbf{r}\frac{\partial f}{\partial \varphi }\dot{\varphi }=-\int \mathrm{d}\mathbf{r}\frac{\partial f}{\partial \varphi }\mathsf{\nabla }·(\varphi {\mathbf{v}}^{\left(p\right)})=\int \mathrm{d}\mathbf{r}{\mathbf{v}}^{\left(p\right)}·\mathsf{\nabla }\Pi \hfill \end{array}$$

(23)

where $\Pi$ is the osmotic pressure given by Equation (11). The second term ${\dot{F}}_{\mathrm{el}}$ is the change rate of the elastic free energy (Equation (2.26) of Ref. [10]) and defines the stress exerted on the polymer ${\mathbf{\sigma }}^{\left(p\right)}$ (called the “network stress” ${\mathbf{\sigma }}^{\left(n\right)}$ in Ref. [10]).

$${\dot{F}}_{\mathrm{el}}=\int \mathrm{d}\mathbf{r}{\mathbf{\sigma }}^{\left(p\right)}:(\mathsf{\nabla }{\mathbf{v}}^{\left(p\right)}).$$

(24)

Some comments on this term are in order:

• This derivation is different from the standard derivation based on the Onsager principle. We normally start with a free energy as a function of the state variables and then calculate the change rate by performing the time derivative.
• Here, the polymer stress ${\mathbf{\sigma }}^{\left(p\right)}$ is input by hand; therefore, we still need a constitutive equation to relate the stress to the state variables.
• Here, the polymer velocity is used. In Ref. [10], it was noted that “it is the deformation of the polymer which causes the change of the free energy”. It turns out that it is very important to specify which velocity is coupled to the polymer stress.

The dissipation function is given by

$$\mathsf{\Phi }=\frac{1}{2}\int \mathrm{d}\mathbf{r}\zeta {({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})}^2,$$

(25)

due to the relative motion of the polymer chain to the solvent.

The Rayleighian is

$$\mathcal{R}=\int \mathrm{d}\mathbf{r}\left[\frac{1}{2}\zeta {({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})}^2+{\mathbf{v}}^{\left(p\right)}·\mathsf{\nabla }\Pi +{\mathbf{\sigma }}^{\left(p\right)}\mathsf{:}\mathsf{\nabla }{\mathbf{v}}^{\left(p\right)}-p\mathsf{\nabla }·(\varphi {\mathbf{v}}^{\left(p\right)}+(1-\varphi ){\mathbf{v}}^{\left(s\right)})\right]$$

(26)

which gives

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta {\mathbf{v}}^{\left(p\right)}}=0& \Rightarrow & \zeta ({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})-\mathsf{\nabla }·{\mathbf{\sigma }}^{\left(p\right)}+\mathsf{\nabla }\Pi +\varphi \mathsf{\nabla }p=0,\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta {\mathbf{v}}^{\left(s\right)}}=0& \Rightarrow & \zeta ({\mathbf{v}}^{\left(s\right)}-{\mathbf{v}}^{\left(p\right)})+(1-\varphi )\mathsf{\nabla }p=0.\hfill \end{array}$$

(28)

These are Equations (2.28) and (2.29) in Ref. [10]. Here, we have only two time evolution equations. To complete the formulation, a constitutive equation needs to be specified. This is slightly different from the previous derivation, where the constitutive Equation (20) is derived from the condition $\delta \mathcal{R}/\delta D_t^{\left(p\right)}{c}_{\alpha \beta }=0$.

3.2. Polymer Blends

For binary mixtures of the polymer melt, Doi and Onuki suggested to use the “tube velocity” (or friction-averaged velocity; see also Tanaka’s works [21,22]) in Equation (24)

$${\mathbf{v}}_{\mathrm{T}}=\frac{1}{{\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}}}({\zeta }_{\mathrm{L}}{\mathbf{v}}_{\mathrm{L}}+{\zeta }_{\mathrm{S}}{\mathbf{v}}_{\mathrm{S}})$$

(29)

where L and S stand for the long and short polymers.

The Rayleighian is written as

$$\begin{array}{ccc}\hfill \mathcal{R}& =& \int \mathrm{d}\mathbf{r}[\frac{1}{2}\zeta {({\mathbf{v}}_{\mathrm{L}}-{\mathbf{v}}_{\mathrm{S}})}^2-\mu \mathsf{\nabla }·\left({\varphi }_{\mathrm{L}}{\mathbf{v}}_{\mathrm{L}}\right)+{\mathbf{\sigma }}^{\left(n\right)}:\mathsf{\nabla }{({\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}})}^{-1}({\zeta }_{\mathrm{L}}{\mathbf{v}}_{\mathrm{L}}+{\zeta }_{\mathrm{S}}{\mathbf{v}}_{\mathrm{S}})\hfill \\ & & -p\mathsf{\nabla }·({\varphi }_{\mathrm{L}}{\mathbf{v}}_{\mathrm{L}}+(1-{\varphi }_{\mathrm{L}}){\mathbf{v}}_{\mathrm{S}})]\hfill \end{array}$$

(30)

which gives

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta {\mathbf{v}}_{\mathrm{L}}}=0& \Rightarrow & \zeta ({\mathbf{v}}_{\mathrm{L}}-{\mathbf{v}}_{\mathrm{S}})-\frac{{\zeta }_{\mathrm{L}}}{{\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}}}\mathsf{\nabla }·{\mathbf{\sigma }}^{\left(n\right)}+{\varphi }_{\mathrm{L}}\mathsf{\nabla }\mu +{\varphi }_{\mathrm{L}}\mathsf{\nabla }p=0,\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill \frac{\delta \mathcal{R}}{\delta {\mathbf{v}}_{\mathrm{S}}}=0& \Rightarrow & \zeta ({\mathbf{v}}_{\mathrm{S}}-{\mathbf{v}}_{\mathrm{L}})-\frac{{\zeta }_{\mathrm{S}}}{{\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}}}\mathsf{\nabla }·{\mathbf{\sigma }}^{\left(n\right)}+(1-{\varphi }_{\mathrm{L}})\mathsf{\nabla }p=0.\hfill \end{array}$$

(32)

These are Equations (4.3) and (4.4) of Ref. [10].

4. Dumbbell Model

A popular model of viscoelastic fluid is the Oldroyd-B model, also known as the dumbbell model [1,2]. It has been shown [3] that this model can be derived from an energetic principle similar to the Onsager principle. In our previous work [20], we derived a Rayleighian that gives the Oldroyd-B model, and showed that such a formulation is useful to obtain analytical solutions for certain problems. Here, we extend the framework to the two-fluid model, and derive a set of equations that accounts for the coupling of the flow and the diffusion.

A polymer chain is modeled as a dumbbell consisting of two beads connected by an elastic spring that has an end-to-end vector $\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2$. The conformation of the dumbbell is specified by the $\mathit{c}$-tensor defined by $\mathit{c}=\frac{k}{k_{B}T}\langle \mathbf{r}\mathbf{r}\rangle$.

4.1. Free Energy

For dilute solutions, the free energy function is given by

$$\begin{array}{ccc}\hfill a(\varphi ,\mathit{c})& =& \frac{k_{B}T}{v}\varphi ln\varphi +\varphi g\left(\mathit{c}\right),\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill g\left(\mathit{c}\right)& =& \frac{1}{2}\frac{k_{B}T}{v}\left[\mathrm{Tr}\left(\mathit{c}\right)-lndet\left(\mathit{c}\right)\right],\hfill \end{array}$$

(34)

where v is the volume of one single dumbbell. The change rate of the free energy is given by Equation (11)

$$\dot{A}=\int \mathrm{d}\mathbf{r}\left[v_{\gamma }^{\left(p\right)}{\nabla }_{\gamma }\Pi +\varphi \frac{\partial g}{\partial c_{\alpha \beta }}{D}_t^{\left(p\right)}{c}_{\alpha \beta }\right].$$

(35)

The variations of $\dot{A}$ with respect to the velocity variables are

$$\begin{array}{ccc}\hfill \frac{\delta \dot{A}}{\delta {\mathbf{v}}^{\left(p\right)}}& =& \mathsf{\nabla }\Pi ,\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill \frac{\delta \dot{A}}{\delta {\mathbf{v}}^{\left(s\right)}}& =& 0,\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \frac{\delta \dot{A}}{\delta D_t^{\left(p\right)}\mathit{c}}& =& \varphi \frac{\partial g}{\partial \mathit{c}}=\frac{1}{2}\frac{k_{B}T}{v}\varphi (\mathit{I}-{\mathit{c}}^{-1}).\hfill \end{array}$$

(38)

4.2. Dissipation Function

The dissipation function related to the $\mathit{c}$-tensor is given in Ref. [20], with $\dot{\mathit{c}}$ replaced by $D_t^{\left(p\right)}\mathit{c}$,

$${\mathsf{\Phi }}_c=\frac{1}{4}\frac{k_{B}T}{v}\int \mathrm{d}\mathbf{r}\tau \varphi \mathrm{Tr}\left[{\mathit{c}}^{-1}·(D_t^{\left(p\right)}{\mathit{c}}^t-\mathbf{\kappa }·\mathit{c}-\mathit{c}·{\mathbf{\kappa }}^t)·(D_t^{\left(p\right)}\mathit{c}-\mathbf{\kappa }·\mathit{c}-\mathit{c}·{\mathbf{\kappa }}^t)\right],$$

(39)

where $\mathrm{Tr}$ denotes the trace operation.

Here arises the main question: which velocity should we use in the velocity gradient $\mathbf{\kappa }$? We have a few options here:

• Using the polymer velocity, $\mathbf{\kappa }={\mathbf{\kappa }}^{\left(p\right)}={(\mathsf{\nabla }{\mathbf{v}}^{\left(p\right)})}^t$. This will give the same results of Doi–Onuki.
• Using the solvent velocity, $\mathbf{\kappa }={\mathbf{\kappa }}^{\left(s\right)}={(\mathsf{\nabla }{\mathbf{v}}^{\left(s\right)})}^t$. This will lead to a different dynamics.
• Using a combination of the polymer and solvent velocities, $\mathbf{\kappa }={(\mathsf{\nabla }\mathbf{V})}^t$, with

  $$\mathbf{V}={\alpha }_1{\mathbf{v}}^{\left(p\right)}+{\alpha }_2{\mathbf{v}}^{\left(s\right)},{\alpha }_1+{\alpha }_2=1.$$

  (40)
  This includes a special case of volume-average velocity $\mathbf{V}=\mathbf{v}$, with ${\alpha }_1=\varphi$ and ${\alpha }_2=1-\varphi$.

We will continue the derivation using the combination velocity $\mathbf{V}$. In general, one should expect ${\alpha }_i$ to be a function of the concentration $\varphi$, so it will be position-dependent.

$$\mathbf{\kappa }={(\mathsf{\nabla }\mathbf{V})}^t={\left[\mathsf{\nabla }({\alpha }_1{\mathbf{v}}^{\left(p\right)}+{\alpha }_2{\mathbf{v}}^{\left(s\right)})\right]}^t.$$

(41)

The total dissipation function is given by

$$\mathsf{\Phi }={\mathsf{\Phi }}_c+\frac{1}{2}\int \mathrm{d}\mathbf{r}\zeta {({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})}^2.$$

(42)

The variations with respect to the dynamic variables are

$$\begin{array}{ccc}\hfill \frac{\delta \mathsf{\Phi }}{\delta {\mathbf{v}}^{\left(p\right)}}& =& \zeta ({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})+\mathsf{\nabla }·\left[{\alpha }_1\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)\right],\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill \frac{\delta \mathsf{\Phi }}{\delta {\mathbf{v}}^{\left(s\right)}}& =& \zeta ({\mathbf{v}}^{\left(s\right)}-{\mathbf{v}}^{\left(p\right)})+\mathsf{\nabla }·\left[{\alpha }_2\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)\right],\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill \frac{\delta \mathsf{\Phi }}{\delta D_t^{\left(p\right)}\mathit{c}}& =& \frac{1}{2}\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)·{\mathit{c}}^{-1}.\hfill \end{array}$$

(45)

4.3. Time Evolution Equations

The Rayleighian is written as

$$\mathcal{R}=\dot{A}+\mathsf{\Phi }-\int \mathrm{d}\mathbf{r}p\mathsf{\nabla }·(\varphi {\mathbf{v}}^{\left(p\right)}+(1-\varphi ){\mathbf{v}}^{\left(s\right)}).$$

(46)

The variation of the Rayleighian with respect to $D_t^{\left(p\right)}\mathit{c}$ gives

$$D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})=-\frac{1}{\tau }(\mathit{c}-\mathit{I}).$$

(47)

This is the constitutive equation of the Oldroyd-B fluid model, with a small modification in which the velocity gradient is given by $\mathbf{V}$.

The variation with respect to ${\mathbf{v}}^{\left(p\right)}$ gives

$$\zeta ({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})+\mathsf{\nabla }·\left[{\alpha }_1\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)\right]+\mathsf{\nabla }\Pi +\varphi \mathsf{\nabla }p=0.$$

(48)

The variation with respect to ${\mathbf{v}}^{\left(s\right)}$ gives

$$\zeta ({\mathbf{v}}^{\left(s\right)}-{\mathbf{v}}^{\left(p\right)})+\mathsf{\nabla }·\left[{\alpha }_2\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)\right]+(1-\varphi )\mathsf{\nabla }p=0.$$

(49)

Equations (48) and (49) give

$$\mathsf{\nabla }·\left[{\mathbf{\sigma }}^{\left(p\right)}-(\Pi +p)\mathit{I}\right]=0,$$

(50)

where ${\mathbf{\sigma }}^{\left(p\right)}$ is defined by

$${\mathbf{\sigma }}^{\left(p\right)}=-\frac{k_{B}T}{v}\tau \varphi \left(D_t^{\left(p\right)}\mathit{c}-{(\mathsf{\nabla }\mathbf{V})}^t·\mathit{c}-\mathit{c}·(\mathsf{\nabla }\mathbf{V})\right)=\frac{k_{B}T}{v}\varphi (\mathit{c}-\mathit{I}).$$

(51)

Equation (50) indicates that the tensor ${\mathbf{\sigma }}^{\left(p\right)}-(\Pi +p)\mathit{I}$ is the total stress tensor and ${\mathbf{\sigma }}^{\left(p\right)}$ is the polymer contribution to the stress tensor. This definition does not depend on the choice of ${\alpha }_i$. For a homogeneous solution, the polymer number density is given by $n_p=\varphi /v$, and then Equation (51) becomes the standard form ${\mathbf{\sigma }}^{\left(p\right)}=n_p{k}_{B}T(\mathit{c}-\mathit{I})=G(\mathit{c}-\mathit{I})$ with the shear modulus $G=n_p{k}_{B}T$.

Using the expression (51), Equations (48) and (49) can then be written as

$$\begin{array}{ccc}\hfill & & \zeta ({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})-\mathsf{\nabla }·({\alpha }_1{\mathbf{\sigma }}^{\left(p\right)})+\mathsf{\nabla }\Pi +\varphi \mathsf{\nabla }p=0,\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill & & \zeta ({\mathbf{v}}^{\left(s\right)}-{\mathbf{v}}^{\left(p\right)})-\mathsf{\nabla }·({\alpha }_2{\mathbf{\sigma }}^{\left(p\right)})+(1-\varphi )\mathsf{\nabla }p=0.\hfill \end{array}$$

(53)

With the setting ${\alpha }_1=1$ and ${\alpha }_2=0$, we recover the Doi–Onuki results for the polymer solutions (27) and (28). With the setting ${\alpha }_1={\zeta }_{\mathrm{L}}/({\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}})$ and ${\alpha }_2={\zeta }_{\mathrm{S}}/({\zeta }_{\mathrm{L}}+{\zeta }_{\mathrm{S}})$, we recover the Doi–Onuki results for the polymer blends (31) and (32).

Combining the above two equations, we can obtain

$$\zeta ({\mathbf{v}}^{\left(p\right)}-{\mathbf{v}}^{\left(s\right)})=-(1-\varphi )\mathsf{\nabla }\Pi +({\alpha }_1-\varphi )\mathsf{\nabla }·{\mathbf{\sigma }}^{\left(p\right)}+(\mathsf{\nabla }{\alpha }_1)·{\mathbf{\sigma }}^{\left(p\right)},$$

(54)

where the last term only appears if ${\alpha }_1$ is position-dependent.

The above equation shows that the relative motion of the polymers with respect to the solvent has two origins: one is the gradient of the osmotic pressure $\mathsf{\nabla }\Pi$, which corresponds to the usual diffusion due to the concentration gradient. The other one is the gradient of the polymer contribution to the stress tensor $\mathsf{\nabla }·{\mathbf{\sigma }}^{\left(p\right)}$. This is the essence of the two-fluid model: the polymer contribution to the stress tensor should induce the polymer diffusion, i.e., the stress and the diffusion are coupled. Using the conservation equation, we can see that the time derivative of the volume fraction has contributions from ${\nabla }^2\Pi$ and $\mathsf{\nabla }\mathsf{\nabla }$:${\mathbf{\sigma }}^{\left(p\right)}$. These two terms are consistent with previous two-fluid models [14]. Existing models have used ${\alpha }_i$ that are independent of the position; therefore, the last term in Equation (54) vanishes.

The magnitude of the stress contribution to the diffusion depends on the choice of ${\alpha }_i$. For ${\alpha }_1=1$ and ${\alpha }_2=0$, the stress contribution is on the order of $O(1-\varphi )$, the same order of the contribution from the osmotic pressure. This is the original Doi–Onuki result [10]. For ${\alpha }_1=0$ and ${\alpha }_2=1$, the stress contribution is reduced to order $O\left(\varphi \right)$, which is small in comparison to the osmotic pressure term. For ${\alpha }_1=\varphi$ and ${\alpha }_2=1-\varphi$, the stress contribution vanishes.

For homogeneous solutions, the polymer concentration is uniform in space; therefore, the polymer velocity ${\mathbf{v}}^{\left(p\right)}$ and the solvent velocity ${\mathbf{v}}^{\left(s\right)}$ are the same, ${\mathbf{v}}^{\left(p\right)}={\mathbf{v}}^{\left(s\right)}$. Different choices in the velocity gradient $\kappa$ lead to the same time evolution equations.

5. Conclusions

We have used the Onsager variational principle to derive the time evolution equations for polymer solutions, taking into consideration the coupling between the stress and the composition. The strength of the current framework is that we start with a microscopic model for the polymer chains, and then the constitutive equation is a natural outcome from the variational calculation. The exact form of the time evolution equations will depend on the choice of dissipation function (see Equation (40)), which then determines the strength of the stress–diffusion coupling. The choice of ${\alpha }_1$ and ${\alpha }_2$ in Equation (40) should be based on the system considered. However, the derived stress tensor (51) is independent of this choice and is determined by the specific polymer model.