Phase-Field (Diffuse Interface)

Presentation of the model

The Phase Field model is a diffuse-interface model dedicated to the direct numerical simulation of incompressible two-phase flows of immiscible fluids. In diffuse-interface models, the interfaces separating the phases of the system are not modelled as surfaces of discontinuity but as volumetric transition zones across which all the physical quantities vary continuously. These models are based on the introduction of an order parameter characteristic of the phases. In the model presented here, this order parameter is representative of the volume fraction of one of the phases, denoted \(\varphi\). In order to take into account the capillary phenomena, important in direct numerical simulation, the energy of the system is assumed to depend on \(\varphi\) but also on \(\nabla \varphi\), in the following form:

\[E = W\left(\varphi\right)+\frac{\alpha}{2}\left(\nabla \varphi\right)^2 \]

where \(\alpha\) is called the internal capillarity coefficient. The function \(W\left(\varphi\right)\) is a double-well function whose simplest form is:

\[W\left(\varphi\right)=\beta\varphi^2\left(1-\varphi\right)^2 \]

where \(\beta\) is a constant characteristic of the model. The coefficients \(\alpha\) and \(\beta\) are characteristic of the interfacial tension \(\sigma\) and of the thickness \(h\) of the interfaces. Indeed, the following relations hold:

\[\sigma=\frac{\sqrt{2\alpha\beta}}{6} \]

\[h=4\sqrt{\frac{\alpha}{2\beta}} \]

The interfacial tension \(\sigma\) is a characteristic of the fluids present, and the thickness \(h\) is generally chosen equal to about 4 or 5 times the mesh size.

By application of the second principle of thermodynamics, it can be shown that the equations of motion are as follows:

\[\frac{d\varphi}{dt}=\nabla\cdot\left[\kappa\left(\varphi\right)\nabla\left(\tilde{\mu}+\frac{d\rho}{d\varphi}\frac{u^2}{2}\right)\right] \]

\[\rho\left(\varphi\right)\frac{du}{dt}=-\nabla P + \tilde{\mu}\nabla\varphi+\rho\left(\varphi\right)g+\nabla\cdot\tau \]

\[\nabla\cdot u=0 \]

where

\[\tilde{\mu}=\frac{dW}{d\varphi}-\nabla\cdot\left(\alpha\nabla\varphi\right) \]

is the generalised chemical potential and

\[\rho\left(\varphi\right)=\rho_{1}+\varphi\left(\rho_{2}-\rho_{1}\right) \]

is the density, \(\rho_{1}\) and \(\rho_{2}\) representing the densities of the phases corresponding to \(\varphi=0\) and \(\varphi=1\) respectively. In the previous system, \(d/dt\) represents the convective derivative, \(P\) is the pressure and \(\tau\) the viscous stress tensor.

The evolution equation of \(\varphi\) is called the Cahn-Hilliard equation. The coefficient \(\kappa\) appearing in this equation is called the mobility and may depend on \(\varphi\) quadratically:

\[\kappa\left(\varphi\right)=a\kappa_{0}\varphi\left(1-\varphi\right) \]

where \(a\) and \(\kappa_{0}\) are constants. Note that the term \(\tilde{\mu}\nabla\varphi\) appearing in the momentum balance equation represents the capillary forces.

Change of variable

For reasons of programming simplicity, the Cahn-Hilliard equation is treated as a convection-diffusion equation of concentration. The right-hand side of the Cahn-Hilliard equation is treated as a source term, as is the \(\tilde{\mu}\nabla\varphi\) term of the momentum balance equation. For reasons of symmetry, the "concentration", i.e. the unknown of the Cahn-Hilliard equation, is defined by:

\[c=\varphi-\frac{1}{2} \]

\(c\) then varies from \(-1/2\) to \(+1/2\) across an interface. The closure relations \(W\left(\varphi\right)\), \(\rho\left(\varphi\right)\) and \(\kappa\left(\varphi\right)\) must therefore be modified accordingly:

\[W\left(\varphi\right)=\beta\left(c-\frac{1}{2}\right)^2\left(c+\frac{1}{2}\right)^2 \]

\[\rho\left(\varphi\right)=\frac{\rho_{1}+\rho_{2}}{2}+c\left(\rho_{2}-\rho_{1}\right) \]

\[\kappa\left(\varphi\right)=a\kappa_{0}\left(\frac{1}{2}+c\right)\left(\frac{1}{2}-c\right) \]

Boussinesq approximation

The general model can be degraded by using the Boussinesq approximation, in which the density is assumed constant except in the gravity term, where it is linearised. This approximation is the one classically used in TrioCFD, and the variables used for the numerical resolution of the resulting system of equations are those classically used in TrioCFD. In particular, the "pressure" used for the coupled resolution of the momentum balance and zero-divergence-of-velocity equations is in fact \((P/\rho_{m}-gz)\), where \(\rho_{m}=(\rho_{1}+\rho_{2})/2\) is the mean density and \(z\) the vertical coordinate (along the direction of gravity). Given the form of the function \(\rho(c)\), the source term linked to gravity has the following expression:

\[\frac{\rho_{2}-\rho_{1}}{\rho_{m}}\, c\, g \]

This form makes it possible to define the value of the parameter beta_co (i.e. \((\rho_{2}-\rho_{1})/\rho_{m}\)).

Forms of the momentum balance

Different forms of the momentum balance are possible, and four are currently available. Although these forms are equivalent at the continuous level, they may potentially give different results because of the spatial discretization and numerical resolution. These forms modify the definition of the "pressure" used for the resolution of the Poisson equation during the resolution of the momentum balance (see the numerical resolution of the momentum balance below). The available forms are the following:

Form 1

\[-\nabla P + \tilde{\mu}\nabla c = -\nabla\left(P - \tilde{\mu}c\right) - c\nabla\tilde{\mu} \]

Form 2

\[-\nabla P + \tilde{\mu}\nabla c = -\nabla\left(P+W+c\nabla\cdot\left(\alpha\nabla c\right)\right)+c\nabla\left[\nabla\cdot\left(\alpha\nabla c\right)\right] \]

Form 3

\[-\nabla P + \tilde{\mu}\nabla c = -\nabla\left(P+W-\alpha\frac{\left(\nabla c\right)^2}{2}+c\nabla\cdot\left(\alpha\nabla c\right)\right)+c\nabla\left[\nabla\cdot\left(\alpha\nabla c\right)\right]-\nabla\left(\alpha\frac{\left(\nabla c\right)^2}{2}\right) \]

Form 4

\[-\nabla P + \tilde{\mu}\nabla c = -\nabla\left(P+W\right)-\nabla\cdot\left(\alpha\nabla c\right)\nabla c \]

The first form is generally preferred to the others, which are more historical. This form has the numerical advantage of having a zero source term at equilibrium, the latter being in particular characterised by the condition \(\tilde{\mu}=\text{cste}\). On the other hand, the "pressure" \((P-\tilde{\mu}c)\) is constant across a curved interface at equilibrium, which does not correspond to the physicists' usual pressure; the physicists' usual pressure is in fact \(P\).

Numerical resolution

The system of equations is solved numerically in two main steps:

1. resolution of the Cahn-Hilliard equation, which makes it possible to determine \(\varphi^{n+1}\) and \(\tilde{\mu}^{n+1}\);
2. resolution of the momentum balance equation coupled to the zero-divergence constraint by a projection method, which makes it possible to determine \(u^{n+1}\) and \(P^{n+1}\).

The resolution of the Cahn-Hilliard equation can be done either explicitly or implicitly. The implicit resolution is performed either by a fixed-point method or by a Newton-Krylov method and essentially makes it possible to increase the stability time step (sometimes significantly). The time-discretization schemes of the Cahn-Hilliard and momentum balance equations may be different.

Cahn-Hilliard equation

The time advance of the \(c\) field is performed in two main sub-steps:

1. determination of \(c^{n+1/2}\) taking into account only the right-hand side of the Cahn-Hilliard equation; this step can be solved either explicitly or implicitly.
2. determination of \(c^{n+1}\) taking into account the convection term; it is for this second step that the time scheme (Euler, Runge-Kutta, etc.) can be chosen.

The explicit resolution of the first sub-step solves the following equation:

\[\frac{c^{n+1/2}-c^n}{\Delta t}=\nabla\cdot\left[\kappa\left(c^n\right)\nabla\left(\tilde{\mu}^n+\frac{d\rho}{d\varphi}\frac{\left(u^n\right)^2}{2}\right)\right] \]

The implicit resolution deserves more explanation because it is based on the use of a semi-implicit approximation of \(\tilde{\mu}\). For any field \(\phi\), we denote:

\[\phi^{n+1/4}=\theta\phi^{n+1/2}+\left(1-\theta\right)\phi^{n},\quad \text{with } \theta \in \left[0;1\right] \]

For reasons of stability and accuracy, the value of \(\theta\) is set to 0.6.

\[\frac{c^{n+1/2}-c^n}{\Delta t}=\nabla\cdot\left[\kappa\left(c^n\right)\nabla\left(\tilde{\mu}^{n+1/4}+\frac{d\rho}{d\varphi}\frac{\left(u^n\right)^2}{2}\right)\right] \]

with

\[\tilde{\mu}^{n+1/4}=W\left(c^{1/4}\right)+\nabla\cdot\left(\alpha\nabla c^{n+1/4}\right) \]

This system is nonlinear (because of the nonlinearity of \(W(c)\)) and is solved iteratively either by a fixed-point method or by a Newton-Krylov method (also called nonlinear GMRES method). Note that, for reasons of simplification of the nonlinear system, the mobility is always determined at time step \(n\) and not at time step \((n+1)\).

Momentum balance equation

The momentum balance equation, coupled to the zero-divergence constraint of the velocity, is solved by a projection method whose main steps are recalled here. To simplify, only an explicit Euler scheme is considered, and the momentum balance equation is put in the following form:

\[\rho\frac{\partial u}{\partial t}=-\nabla P + F \]

where the vector \(F\) represents all the forces other than that of pressure. This equation is discretized as follows:

\[\rho^n \frac{u^{n+1}-u^n}{\Delta t}=-\nabla P^{n+1}+F \]

where the time step at which \(F\) is evaluated is not specified because it may depend on the problems treated; however, the dependencies of \(F\) on \(u\) must be explicit. After division by \(\rho^n\), application of the divergence operator and accounting for the zero-divergence condition, the following Poisson equation on the pressure is obtained:

\[\nabla\cdot\left(\frac{1}{\rho^n}\nabla P^{n+1}\right)=\nabla\cdot\left(\frac{F}{\rho^n}\right) \]

The pressure matrix associated with the discretization of this equation depends on the field \(\rho^n\) and must therefore be assembled at each time step. This is the reason why the methods related to the assemblers had to be modified.

Generalization of the closure laws (Navier-Stokes / Cahn-Hilliard binary model)

Since version v1.8.3, the Phase Field baltik has had changes in order to generalize the thermodynamic landscape, the density law and so on of this binary Cahn-Hilliard model coupled with the (isothermal) Navier-Stokes equations.

Chemical potential function \(\frac{dW}{dc}\)

The chemical potential function \(\frac{dW}{dc}\) is associated with the homogeneous contribution \(W\) (so-called "thermodynamic landscape") function of \(c=\varphi - \frac{1}{2}\) to the energy density \(E\).

TrioCFD v1.8.2: the potential has the following prescribed form:

\[\frac{dW}{dc} = \beta \left( 4 c \left(c-\frac{1}{2}\right) \left(c+\frac{1}{2}\right)\right) \]

Modification since v1.8.3: the potential is written as:

\[\frac{dW}{dc} = \beta \times \begin{cases} 4 c \left(c-\frac{1}{2}\right) \left(c+\frac{1}{2}\right) & \text{if } \texttt{Op}_W = \texttt{defaut} \\ \texttt{f}_W & \text{if } \texttt{Op}_W = \texttt{fonction}~\texttt{f}_W \end{cases} \]

where \(\texttt{f}_W\) is a univariate function (as read by lire_f of the Table class of TRUST) whose variable \(c\) should be referred to as var, e.g. (4.*val*(val+0.5)*(val-0.5)).

Syntax evolution — input file excerpt in v1.8.2:

Convection_diffusion_Phase_field {
  ...
  sources {
    Source_Con_Phase_field {
      ...
      beta <beta>
      ...
    }
  }
  ...
}

Input file excerpt since v1.8.3:

Convection_diffusion_Phase_field {
  ...
  sources {
    Source_Con_Phase_field {
      ...
      beta <beta>
      potentiel_chimique { Op_W }
      ...
    }
  }
  ...
}

Kinetic parameter ("mobility") \(\kappa\)

The kinetic parameter (so-called "mobility") \(\kappa\) is either a constant or a function of \(c\) in the Cahn-Hilliard equation.

TrioCFD v1.8.2: the mobility is written as:

\[\kappa = \kappa_0 \times \begin{cases} 1 & \text{if } \texttt{Op}_K = \texttt{non} \\ a \left(\frac{1}{2}+c\right) \left(\frac{1}{2}-c\right) & \text{if } \texttt{Op}_K = \texttt{oui} \end{cases} \]

Modification since v1.8.3: the mobility is written as:

\[\kappa = \kappa_0 \times \begin{cases} 1 & \text{if } \texttt{Op}_K = \texttt{non} \\ a \left(\frac{1}{2}+c\right) \left(\frac{1}{2}-c\right) & \text{if } \texttt{Op}_K = \texttt{defaut} \\ a\, \texttt{f}_K & \text{if } \texttt{Op}_K = \texttt{fonction}~\texttt{f}_K \end{cases} \]

where \(\texttt{f}_K\) is a univariate function (as read by lire_f of the Table class of TRUST) whose variable \(c\) should be referred to as var.

Syntax evolution — input file excerpt in v1.8.2:

Convection_diffusion_Phase_field {
  ...
  sources {
    Source_Con_Phase_field {
      ...
      kappa <kappa_0>
      kappa_variable Op_K
      multiplicateur_de_kappa <a>
      ...
    }
  }
  ...
}

Input file excerpt since v1.8.3 (kappa_variable { Op_K }):

Convection_diffusion_Phase_field {
  ...
  sources {
    Source_Con_Phase_field {
      ...
      kappa <kappa_0>
      kappa_variable { Op_K }
      multiplicateur_de_kappa <a>
      ...
    }
  }
  ...
}

Mass density \(\rho\)

The mass density \(\rho\) is described as a function of \(c\).

TrioCFD v1.8.2: the mass density has the following prescribed form:

\[\rho = \frac{\rho_1+\rho_2}{2}+c\left(\rho_2-\rho_1\right) \]

Note

With the concentration convention used here, \(c=\varphi-\frac{1}{2}\in[-\frac{1}{2},\frac{1}{2}]\), this is the same density law as in the Change of variable section above: it returns \(\rho_1\) at \(c=-\frac{1}{2}\) and \(\rho_2\) at \(c=+\frac{1}{2}\), and matches the implementation ( \(\frac{d\rho}{dc}=\rho_2-\rho_1\)).

Modification since v1.8.3: different mass density formulations are available depending on whether the Boussinesq approximation (input parameter \(\texttt{Op}_{\text{Boussinesq}}\)) is made ( \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\)) or not ( \(\texttt{Op}_{\text{Boussinesq}} = \texttt{non}\)):

• if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\), it can be written in terms of the fluid properties as \(\rho = \rho_0 \left( 1 + \beta_c c\right)\);
• otherwise, for both values of \(\texttt{Op}_{\text{Boussinesq}}\), it can take the prescribed form \(\rho = \frac{\rho_1+\rho_2}{2}+c\left(\rho_2-\rho_1\right)\);
• or be of the general form \(\rho = \texttt{F}_\rho\), where \(\texttt{F}_\rho\) is a TRUST input field giving \(\frac{d\rho}{dc}\).

Syntax evolution — if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\) (not available in v1.8.2); since v1.8.3:

Lire fluide
{
  ...
  rho Champ_Uniforme 1 <rho_0>
  beta_co Champ_Uniforme 1 <beta_c>
  ...
}
...
Navier_Stokes_phase_field {
  ...
  approximation_de_boussinesq oui { }
  ...
}

Whatever the value of \(\texttt{Op}_{\text{Boussinesq}}\) (both v1.8.2 and since v1.8.3):

Convection_diffusion_Phase_field {
  ...
  rho_1 <rho_1>
  rho_2 <rho_2>
  ...
}

If \(\rho\) has a general form (not available in v1.8.2); since v1.8.3:

Navier_Stokes_phase_field {
  ...
  approximation_de_boussinesq Op_Boussinesq {
    rho_fonc_c F_rho F_rho'
  }
  ...
}

Dynamic viscosity \(\mu\) as a function of \(c\)

The dynamic viscosity \(\mu\) is either constant or described as a function of \(c\). The dynamic viscosity (and actually the kinematic viscosity \(\nu\)) is treated differently depending on the input parameters \(\texttt{Op}_{\text{Boussinesq}}\) (Boussinesq approximation) and \(\texttt{Op}_{\text{constant }\mu}\) (constant dynamic viscosity approximation), which can both be set to oui (activated) or non (deactivated):

Navier_Stokes_phase_field {
  ...
  approximation_de_boussinesq Op_Boussinesq
  viscosite_dynamique_constante Op_constant_mu
  ...
}

TrioCFD v1.8.2:

• if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\) or \(\texttt{Op}_{\text{constant }\mu} = \texttt{oui}\), the dynamic viscosity is constant and given by \(\mu = \mu_0\);
• if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{non}\) and \(\texttt{Op}_{\text{constant }\mu} = \texttt{non}\), the dynamic viscosity has the prescribed form \(\mu = \frac{\mu_1+\mu_2}{2}+c\left(\mu_2-\mu_1\right)\).

Note

Mirroring the density law (and the Change of variable convention \(c=\varphi-\frac{1}{2}\in[-\frac{1}{2},\frac{1}{2}]\)), the linear interpolation \(\mu(\varphi)=\mu_1+\varphi(\mu_2-\mu_1)\) gives \(\mu=\frac{\mu_1+\mu_2}{2}+c(\mu_2-\mu_1)\), returning \(\mu_1\) at \(c=-\frac{1}{2}\) and \(\mu_2\) at \(c=+\frac{1}{2}\). This matches the implementation, which builds mu = 0.5*(mu1+mu2) + c*(mu2-mu1) ( \(\frac{d\mu}{dc}=\mu_2-\mu_1\)). The \(\frac{\mu_2-\mu_1}{2}\) form in the model report is a transcription error (the same factor-2 error as in the density law) and has been corrected here.

Modification since v1.8.3:

• if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\) or \(\texttt{Op}_{\text{constant }\mu} = \texttt{oui}\), the dynamic viscosity is constant and given by \(\mu = \mu_0\);
• if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{non}\) and \(\texttt{Op}_{\text{constant }\mu} = \texttt{non}\), the dynamic viscosity can take the prescribed form \(\mu = \frac{\mu_1+\mu_2}{2}+c\left(\mu_2-\mu_1\right)\);
• or be of the general form \(\mu = \texttt{F}_\mu\), where \(\texttt{F}_\mu\) is a TRUST input field.

Syntax evolution — if \(\texttt{Op}_{\text{Boussinesq}} = \texttt{oui}\) or \(\texttt{Op}_{\text{constant }\mu} = \texttt{oui}\) (both v1.8.2 and since v1.8.3):

Lire fluide
{
  ...
  mu Champ_Uniforme 1 <mu_0>
  ...
}

If \(\texttt{Op}_{\text{Boussinesq}} = \texttt{non}\) and \(\texttt{Op}_{\text{constant }\mu} = \texttt{non}\) — v1.8.2:

Convection_diffusion_Phase_field {
  ...
  mu_1 <mu_1>
  mu_2 <mu_2>
  ...
}

Since v1.8.3:

Navier_Stokes_phase_field {
  ...
  viscosite_dynamique_constante non {
    mu_1 <mu_1>
    mu_2 <mu_2>
  }
  ...
}

If \(\mu\) has a general form (not available in v1.8.2); since v1.8.3:

Navier_Stokes_phase_field {
  ...
  viscosite_dynamique_constante non {
    rho_fonc_c F_mu
  }
  ...
}