examples.reactiveWetting.liquidVapor1D

Solve a single-component, liquid-vapor, van der Waals system.

This example solves a single-component, liquid-vapor, van der Waals system as described by Wheeler et al. [5]. The free energy for this system takes the form,

(1)f=eρ2m2+RTm(lnρmv¯ρ)

where ρ is the density. This free energy supports a two phase equilibrium with densities given by ρl and ρv in the liquid and vapor phases, respectively. The densities are determined by solving the following system of equations,

(2)P(ρl)=P(ρv)

and

(3)μ(ρl)=μ(ρv)

where μ is the chemical potential,

(4)μ=fρ

and P is the pressure,

(5)P=ρμf

One choice of thermodynamic parameters that yields a relatively physical two phase system is

>>> molarWeight = 0.118
>>> ee = -0.455971
>>> gasConstant = 8.314
>>> temperature = 650.
>>> vbar = 1.3e-05

with equilibrium density values of

>>> liquidDensity = 7354.3402662299995
>>> vaporDensity = 82.855803327810008

The equilibrium densities are verified by substitution into Eqs. (2) and (3). Firstly, Eqs. (1), (4) and (5) are defined as python functions,

>>> from fipy import CellVariable, Grid1D, TransientTerm, VanLeerConvectionTerm, DiffusionTerm, ImplicitSourceTerm, ConvectionTerm, CentralDifferenceConvectionTerm, Viewer
>>> from fipy.tools import numerix
>>> def f(rho):
...     return ee * rho**2 / molarWeight**2 + gasConstant * temperature * rho / molarWeight * \
...            numerix.log(rho / (molarWeight - vbar * rho))
>>> def mu(rho):
...     return 2 * ee * rho / molarWeight**2 + gasConstant * temperature / molarWeight * \
...            (numerix.log(rho / (molarWeight - vbar * rho)) + molarWeight / (molarWeight - vbar * rho))
>>> def P(rho):
...     return rho * mu(rho) - f(rho)

The equilibrium densities values are verified with

>>> print(numerix.allclose(mu(liquidDensity), mu(vaporDensity)))
True

and

>>> print(numerix.allclose(P(liquidDensity), P(vaporDensity)))
True

In order to derive governing equations, the free energy functional is defined.

F=[f+ϵT2(jρ)2]dV

Using standard dissipation laws, we write the governing equations for mass and momentum conservation,

(6)ρt+j(ρuj)=0

and

(7)(ρui)t+j(ρuiuj)=j(ν[jui+iuj])ρiμNC

where the non-classical potential, μNC, is given by,

(8)μNC=δFδρ=μϵTj2ρ

As usual, to proceed, we define a mesh

>>> Lx = 1e-6
>>> nx = 100
>>> dx = Lx / nx
>>> mesh = Grid1D(nx=nx, dx=dx)

and the independent variables.

>>> density = CellVariable(mesh=mesh, hasOld=True, name=r'$\rho$')
>>> velocity = CellVariable(mesh=mesh, hasOld=True, name=r'$u$')
>>> densityPrevious = density.copy()
>>> velocityPrevious = velocity.copy()

The system of equations is solved in a fully coupled manner using a block matrix. Defining μNC as an independent variable makes it easier to script the equations without using higher order terms.

>>> potentialNC = CellVariable(mesh=mesh, name=r'$\mu^{NC}$')
>>> epsilon = 1e-16
>>> freeEnergy = (f(density) + epsilon * temperature / 2 * density.grad.mag**2).cellVolumeAverage

In order to solve the equations numerically, an interpolation method is used to prevent the velocity and density fields decoupling. The following velocity correction equation (expressed in discretized form) prevents decoupling from occurring,

(9)ui,fc=Afdfdf(ρiμNCfρfi,fμNC)

where Af is the face area, df is the distance between the adjacent cell centers and af is the momentum conservation equation’s matrix diagonal. The overbar refers to an averaged value between the two adjacent cells to the face. The notation i,f refers to a derivative evaluated directly at the face (not averaged). The variable uic is used to modify the velocity used in Eq. (6) such that,

(10)ρt+j(ρ[uj+uic])=0

Equation (10) becomes

>>> matrixDiagonal = CellVariable(mesh=mesh, name=r'$a_f$', value=1e+20, hasOld=True)
>>> correctionCoeff = mesh._faceAreas * mesh._cellDistances / matrixDiagonal.faceValue
>>> massEqn = TransientTerm(var=density) \
...           + VanLeerConvectionTerm(coeff=velocity.faceValue + correctionCoeff \
...                                         * (density * potentialNC.grad).faceValue, \
...                                   var=density) \
...           - DiffusionTerm(coeff=correctionCoeff * density.faceValue**2, var=potentialNC)

where the first term on the LHS of Eq. (9) is calculated in an explicit manner in the VanLeerConvectionTerm and the second term is calculated implicitly as a DiffusionTerm with μNC as the independent variable.

In order to write Eq. (7) as a FiPy expression, the last term is rewritten such that,

ρiμNC=i(ρμNC)μNCiρ

which results in

>>> viscosity = 1e-3
>>> ConvectionTerm = CentralDifferenceConvectionTerm
>>> momentumEqn = TransientTerm(coeff=density, var=velocity) \
...               + ConvectionTerm(coeff=[[1]] * density.faceValue * velocity.faceValue, var=velocity) \
...               == DiffusionTerm(coeff=2 * viscosity, var=velocity) \
...               - ConvectionTerm(coeff=density.faceValue * [[1]], var=potentialNC) \
...               + ImplicitSourceTerm(coeff=density.grad[0], var=potentialNC)

The only required boundary condition eliminates flow in or out of the domain.

>>> velocity.constrain(0, mesh.exteriorFaces)

As previously stated, the μNC variable will be solved implicitly. To do this the Eq. (8) is linearized in ρ such that

(11)μNC=μ+(μρ)(ρρ)ϵTj2ρ

The superscript denotes the current held value. In FiPy, μρ is written as,

>>> potentialDerivative = 2 * ee / molarWeight**2 + gasConstant * temperature * molarWeight / density / (molarWeight - vbar * density)**2

and μ is simply,

>>> potential = mu(density)

Eq. (11) can be scripted as

>>> potentialNCEqn = ImplicitSourceTerm(coeff=1, var=potentialNC) \
...                  == potential \
...                  + ImplicitSourceTerm(coeff=potentialDerivative, var=density) \
...                  - potentialDerivative * density \
...                  - DiffusionTerm(coeff=epsilon * temperature, var=density)

Due to a quirk in FiPy, the gradient of μNC needs to be constrained on the boundary. This is because ConvectionTerm’s will automatically assume a zero flux, which is not what we need in this case.

>>> potentialNC.faceGrad.constrain(value=[0], where=mesh.exteriorFaces)

All three equations are defined and an are combined together with

>>> coupledEqn = massEqn & momentumEqn & potentialNCEqn

The system will be solved as a phase separation problem with an initial density close to the average density, but with some small amplitude noise. Under these circumstances, the final condition should be two separate phases of roughly equal volume. The initial condition for the density is defined by

>>> numerix.random.seed(2011)
>>> density[:] = (liquidDensity + vaporDensity) / 2 * \
...    (1  + 0.01 * (2 * numerix.random.random(mesh.numberOfCells) - 1))

Viewers are also defined.

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewers = Viewer(density), Viewer(velocity), Viewer(potentialNC)
...     for viewer in viewers:
...         viewer.plot()
...     input('Arrange viewers, then press <return> to proceed...')
...     for viewer in viewers:
...         viewer.plot()

The following section defines the required control parameters. The cfl parameter limits the size of the time step so that dt = cfl * dx / max(velocity).

>>> cfl = 0.1
>>> tolerance = 1e-1
>>> dt = 1e-14
>>> timestep = 0
>>> relaxation = 0.5
>>> if __name__ == '__main__':
...     totalSteps = 1e10
... else:
...     totalSteps = 10

In the following time stepping scheme a time step is recalculated if the residual increases between sweeps or the required tolerance is not attained within 20 sweeps. The major quirk in this scheme is the requirement of updating the matrixDiagonal using the entire coupled matrix. This could be achieved more elegantly by calling cacheMatrix() only on the necessary part of the equation. This currently doesn’t work properly in FiPy.

>>> while timestep < totalSteps:
...
...     sweep = 0
...     dt *= 1.1
...     residual = 1.
...     initialResidual = None
...
...     density.updateOld()
...     velocity.updateOld()
...     matrixDiagonal.updateOld()
...
...     while residual > tolerance:
...
...         densityPrevious[:] = density
...         velocityPrevious[:] = velocity
...         previousResidual = residual
...
...         dt = min(dt, dx / max(abs(velocity)) * cfl)
...
...         coupledEqn.cacheMatrix()
...         residual = coupledEqn.sweep(dt=dt)
...
...         if initialResidual is None:
...             initialResidual = residual
...
...         residual = residual / initialResidual
...
...         if residual > previousResidual * 1.1 or sweep > 20:
...             density[:] = density.old
...             velocity[:] = velocity.old
...             matrixDiagonal[:] = matrixDiagonal.old
...             dt = dt / 10.
...             if __name__ == '__main__':
...                 print('Recalculate the time step')
...             timestep -= 1
...             break
...         else:
...             matrixDiagonal[:] = coupledEqn.matrix.takeDiagonal()[mesh.numberOfCells:2 * mesh.numberOfCells]
...             density[:] = relaxation * density + (1 - relaxation) * densityPrevious
...             velocity[:] = relaxation * velocity + (1 - relaxation) * velocityPrevious
...
...         sweep += 1
...
...     if __name__ == '__main__' and timestep % 10 == 0:
...         print('timestep: %e / %e, dt: %1.5e, free energy: %1.5e' % (timestep, totalSteps, dt, freeEnergy))
...         for viewer in viewers:
...             viewer.plot()
...
...     timestep += 1
>>> from fipy import input
>>> if __name__ == '__main__':
...     input('finished')
>>> print(freeEnergy < 1.5e9)
True
Last updated on Feb 14, 2025. Created using Sphinx 7.1.2.