The free energy of this model contains a term that is derived from an equation of state. The equation of state and its parametrization determines the density of the liquid and gas phase.
Here we use the Carnahan-Starling EOS, with the parametrization from the paper
Next, we use a function that determines the gas and liquid density, using the Maxwell construction rule. This function uses an iterative procedure, that terminates when the equal-area condition is fulfilled with a certain tolerance.
Now the free energy depends only on ρ and φ. This transformed form is later used to derive expressions for the chemical potential, pressure tensor and force computations.
We define one function that takes an expression with derivative objects in it, substitutes the spatial derivatives with finite differences using the strategy defined in the `fd_discretization` function and compiles a kernel from it.
%% Cell type:code id: tags:
``` python
defmake_kernel(assignments):
# assignments may be using the symbols ρ and φ
# these is substituted with the access to the corresponding fields here
field_substitutions={
ρ:ρ_field.center,
φ:φ_field.center
}
processed_assignments=[]
forainassignments:
new_rhs=a.rhs.subs(field_substitutions)
# ∂∂f representing the laplacian of f is replaced by the explicit carteisan form
# ∂_0 ∂_0 f + ∂_1 ∂_1 f (example for 2D)
# otherwise the discretization would not do the correct thing
new_rhs=replace_generic_laplacian(new_rhs)
# Next the "∂" objects are replaced using finite differences
In the next cell the kernel to compute the chemical potential is created. First an analytic expression for μ is obtained using the free energy, which is then passed to the discretization function above to create a kernel from it. We only have to store the chemical potential of the φ coordinate explicitly, which enters the Cahn-Hilliard lattice Boltzmann for φ.
For the pressure tensor a trick for enhancing numerical stability is used: the bulk component is not stored directly in the pressure tensor field, but the related quantity called `pbs` is stored in a separate field.
$ pbs = \sqrt{|ρ c_s^2 - p_{bulk} |} $
The force is then calculated as $ \nabla \cdot P_{if} + 2 (\nabla pbs) pbs$
In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression.
#### Lattice Boltzmann schemes for time evolution of ρ and φ
- ρ is handled by a normal LB method (compressible, entropic equilibrium)
- stream and collide are splitted into separate kernels
- macroscopic values are computed after the stream, but inside the stream kernel
- velocity field stores the velocity which was not corrected for the forces yet
- the φ collision kernel corrects the velocity itself, because then the updated forces are used for the correction. When u is computed, the updated forces are not computed yet
- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system
- exact difference method is used to couple the force into the ρ-LBM
%% Cell type:markdown id: tags:
The following cell handles the clipping of the order parameters:
%% Cell type:code id: tags:
``` python
ifclipping:
defclip(ac,symbol,min_value,max_value):
"""Function to clip the value of a symbol which is on one of lhs of the assignments
