Hybrid Phase Field Model for Molten Solidification

Contained in this repository is a simulation for molten solidification. The goal is to simulate a droplet landing on a substrate. It is trying to capture phase transitions and how different phases react with eachother. The model makes a few key assumptions over a regular fluid dymanics simulation:

• The density of air, liquid and solid are all the same
• The force term due to gravity is ignored for air

All of the code is written in C++ for version 17.

Directory Layout

├── README.md                    # Project overview and build instructions
├── include
|   ├── shared_structs.h         # Includes shared definitions, structs, and Scalar Field and Velocity Field class implementations (class function implementations are here too...)
│   ├── *.h                      # Other header files
├── input.txt                    # txt file for inputting simulation parameters
└── src
    ├── bc_functions.cpp         # Implementation of boundary condition functions (Neumman and No-Slip)
    ├── discrete_functions.cpp   # Implementation of discretised functions
    ├── initialisation.cpp       # Implementation of initilisations functions
    ├── main.cpp                 # Programme Entry Point
    ├── output_functions.cpp     # Functions for outputting data
    ├── update_functions.cpp     # Functions for updating fields
    └── update_velocity.cpp      # Functions for updating velocity fields

Equations

The simulation works by having three key areas:

• Air area ($\Omega_a$)
• Liquid area ($\Omega_l$)
• Solid area ($\Omega_s$)

These are described by two different fields:

• $\psi$ describes the part of the simulation that is either $\Omega_a$ or $\Omega_l \cup \Omega_s$
• $\varphi$ describes the part of the simulation that is either $\Omega_l$ or $\Omega_s$

Phase Field Equations

The values the fields take are as follows:

$$\psi = \begin{cases} 0, & \text{for air} \\1, & \text{for droplet} \end{cases}$$ $$\varphi = \begin{cases} 0, & \text{for solid} \\1, & \text{for liquid} \end{cases}$$

Psi-Field

The $\psi$ field is described by a standard Cahn-Hilliard equation:

$$D_t \psi = M_{\psi} \nabla^2 \mu_{\psi}$$

Where

• $D_t$ is the material derivative defined as by $D_t = \frac{\partial}{\partial t} + \vec{v} \cdot \nabla$
• $M_{\psi}$ is the diffusion coefficient
• And $\mu_{\psi}$ is the chemical potential of the system

$\mu_{\psi}$ is given by

$$\mu_{\psi} = a_{\psi} g'(\psi) - k_{\psi} \nabla^2 \psi$$

Where

• The $a_{\psi}$ is a constant describing the free energy associated with the bulk energy
• $k_{\psi}$ is a constant describing the free energy assocuiated with boundaries
• $g'(\psi)$ is the derivative of function the $g(\psi)$ with minima at 0 and 1 given by $g(\psi) = \psi^2(1-\psi)^2$

Phi-Field

To describe the $\varphi$ field, an Allen-Cahn phase field equation is used, where the evolution of $\varphi$ is given by

$$D_t \varphi = - M_{\varphi} (a_{\varphi} g'(\varphi) - k_{\varphi} \nabla^2 \varphi + \frac{\rho_l L}{T_M}p'(\varphi)(T_M - T))$$

To enfore the boundary condition, a factor of \psi is added to the above equation:

$$D_t (\psi\varphi) = - M_{\varphi} (a_{\varphi}\psi g'(\varphi) - k_{\varphi} \nabla \cdot (\psi \nabla \varphi) + \frac{\rho_l L}{T_M}\psi p'(\varphi)(T_M - T))$$

Where

• $a_{\varphi}$ and $k_{\varphi}$ are analogous to the ones given in the $psi$ field equation
• $L$ is the latent heat
• $T$ is the temprature
• $T_{M}$ is the melting point of the substance
• $p'(\varphi)$ is the derivative of an interpolation function that goes between 0 and 1

Specifically $p'(\varphi)$ is described by

$$p'(\varphi) = 30\varphi^2(\varphi - 1)^2$$

With the special condition

$$p'(\varphi) = \begin{cases} 0, & \varphi < 0\\\ 0, & \varphi > 1\\\ 1, & \varphi = 1 \text{ and } T > T_M\\\ 1, & \varphi = 1 \text{ and } T < T_M\\\ p'(\varphi), & \text{else} \end{cases}$$

Heat Field and Property Fields

The heat diffusion equation is given by

$$D_t T = \nabla \cdot (D \nabla T)$$

Heat diffusion ($D$), density ($\rho$) and viscosity ($\eta$) are given by the general property field:

$$\beta(\psi, \psi\varphi) = \beta_g + (\beta_s - \beta_g)\psi + (\beta_l - \beta_s)\psi\varphi$$

For material properties (density, viscosity and diffusion) the above equation is used.

• $\beta_g$ is the value of that property in pure air
• $\beta_s$ is the value in pure solid
• $\beta_l$ is the value in pure liquid

Momentum Equation

The momentum equation is given by

$$\frac{d (\rho \textbf{u})}{dt} + \rho(\textbf{u} \cdot \nabla)\textbf{u} = - \nabla P + \nabla \cdot [\eta (\nabla \mathbf{u} + \nabla \mathbf{u}^T)] + \psi \rho \mathbf{g} + \mathbf{f_s}$$

with the surface term being given by

$$\mathbf{f_s} = \varphi \mu_{\psi} \nabla \psi$$

Implementation and Discretisation

The field equations are discretised using a simple explicit euler step method in a 3D grid.

For velocity, a projection step method is used.

Discrete Operators

The spacial derivatives are implemented using a second order central finite difference scheme.

$$\nabla^2 \phi_{i,j,k} \approx \frac{1}{h^2} \left( \phi_{i+1,j,k} + \phi_{i-1,j,k} + \phi_{i,j+1,k} + \phi_{i,j-1,k} + \phi_{i,j,k+1} + \phi_{i,j,k-1} - 6\phi_{i,j,k} \right)$$ $$\frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1} - \phi_{i-1}}{2\Delta x}$$ $$\frac{\partial^2 \phi}{\partial x^2} \approx \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{\Delta x^2}$$

Phase Fields and Heat Equation

The heat equation and psi are the easiest to solve, just being

$$T^{n+1} = T^{n} + dt (\widetilde{\nabla} \cdot (D^{n} \widetilde{\nabla}T^{n}) + \textbf{u} \cdot \widetilde{\nabla} T^{n})$$

For $T$, a Dirichlet boundary condition is used.

$$\psi^{n+1} = \psi^{n} + dt(M_{\psi} \widetilde{\nabla}^2 \mu_{\psi}^{n})$$

Solving for $\varphi$ is a bit more complicated. PhiPsi can be thought of as its own seperate field. First an interim value is calculated.

$$(\psi\varphi)^* = (\psi\varphi)^{n} + dt(- M_{\varphi} (a_{\varphi}\psi g(\varphi) - k_{\varphi} \widetilde{\nabla} \cdot (\psi \widetilde{\nabla} \varphi) + \frac{L}{T_M}\psi p'(\varphi)(T_M - T)))$$

Then $\varphi$ is soled for from $(\psi\varphi)^*$

$$\varphi^{n+1} = \frac{(\psi\varphi)^*}{\psi^{(n+1)}}$$

After calculating $\varphi^{n+1}$, it is clamped between the values 0 and 1.

Then we get $(\psi\varphi)^{n+1}$ from

$$(\psi\varphi)^{n+1} = \psi^{n+1}\varphi^{n+1}$$

For $\psi$ and $\varphi$, Neumman boundary conditions are used.

When property fields are updated ($D$, $\rho$, etc), after the $\psi$ and $\varphi$ fields have been, Dirilicht Boundary conditions are used for them also.

Momentum Equation

Velocity is defined at the centre of each grid point. At each time step, the velocity field is solved for using a projection step method:

$$\frac{d (\rho \textbf{u}^*)}{dt} + (\textbf{u}^{n} \cdot \widetilde{\nabla})\textbf{u}^{n} = + \widetilde{\nabla} \cdot [\eta (\widetilde{\nabla} \mathbf{u}^{n} + \widetilde{\nabla} (\mathbf{u}^{n})^T)] + \psi \rho \mathbf{g} + \widetilde{\mathbf{f_s}}$$

$$\widetilde{\nabla}^2 P' = \frac{\rho}{dt} \widetilde{\nabla} \cdot \mathbf{u}^*$$

$$\mathbf{u}^{n+1} = \mathbf{u}^* - \frac{dt}{\rho}\widetilde{\nabla} P'$$

The No-Slip boundary condition is applied to the new velocity field such at, at the edge of the box, the velocity in the ghost cells have will be the negative of the velocity bordering them in the box.

Field and Configuration Data Structures

Fields are stored in classes. They are defined in shared_structs.h. VectorField3D is used for vector fields (this is only used for the velocity field), while ScalarField3D is used for scalar fields (every other field). They both use flattened arrays for increased performance. A side effect of this, is that in order for index accessing to work (e.g (i, j, k)), W > L > H. All of their class methods have been implemented within the header file, for now. Also, within shared_structs.h, there is a cfg struct that stores the initial values of constants from the input file, as well as definitions for selecting only inner cells (non-ghost cells) and all cells. Additionally, there is a box struct which defines the dimensions of the box.

Within initialisation.h, there are four important structs that store field data:

• PsiFieldData stores all the field data (constants, previous fields, chemical potentials, fluxes etc.) for $\psi$
• PsiPhiFieldData is similar but stores the data for both $\varphi$ and $\varphi\psi$
• TFieldData stores the field data for both the T field and the diffusion field along with the diffusion value constants
• vFieldData stores the velocity field, the pressure field, the viscosity field and the density fields

Each of them has a constructor method which requires at least the cfg struct and the box struct to generate themselves.

PsiPhiFieldData needs psi_field_n to generate, which is a member variable of PsiFieldData. TFieldData and vFieldData both require psi_field_n and psi_phi_field_n, the latter being a member variable of PsiPhiFieldData respectively. Because of this, PsiFieldData must be generated first, followed by PsiPhiFieldData and then the other two can be generated.

The point of this is to keep the code less cluttered.

Using the Simulation

Input File

The simulation has an input.txt file. This will give the values of the constants within the simulation:

Parameter(s)	Description
W, L, H	Box sizes in x, y, z. Padded by +2 internally
R, x, y, z	Initial droplet radius and position (in grid units)
T_droplet	Initial temperature of the droplet
T_init	Initial temperature of rest of the box and walls
T_star	Melting Temperature of solid phase
M, a, k	Parameters for ψ and φ phase field equations
Da, Dl, Ds	Diffusivity values for air, liquid, and solid
Dg	Diffusivity for the bottom plate (not implemented
g	Acceleration due to gravity
rho_a, rho_s, rho_l	Densities (must be equal currently)
visc_a, visc_s, visc_l	Viscosities (can be different)
dx	Cell width
dt	Timestep size
steps	Total number of steps
save_step	Interval to save output

Pressure Solver

The way the Pressure Solver uses a Gauss-Seidel Red-Black Successive over Relaxation (SOR) method.

It is important to note that alot of values for the pressure solver have been hard coded in.

These include:

• "tolerance" the minimum difference between $P_{n+1}$ and $P_n$ for the solver to decide it has converged to a suffiecent P field. This has been set to a default value of 0.1. Note though this is the difference across both of the entire fields
• "max_tries" is the number of attempts that solver will try to converge $P_n$ to. If $|P_{n+1} - P_{n}|$ ($\Delta P_n$) does not go below tolerance in less steps than max tries, then the solver will give up and just take the field it is on when max_tries is reached. In the code this is set to 1,000,000
• "omega" is the SOR factor. Its value has just sort of been guessed as 1.7

The Pressure field has neumann boundary conditions. Because of this, it is important to fix the value of P at some point in the field. Otherwise, the solver will never converge. The point fixed in the code is $P(3,3,3) = 0$. This should not affect too much as the momentum equation only cares about pressure gradients, though for trying to gain information about the pressure, it is problematic.

Output folder

Currently, the sim will export fields into their respective folder within the main output folder:

output/
  Ds/                               Diffusion Field Values
  PhiPsis/                          ψφ Field Values
  Phis/                             φ Field Values
  Ps/                               Pressure Field Values
  Psis/                             ψ Field Values
  Ts/                               Temperature Field Values
  vs/                               Magnitude of velocity Field
  vxs/                              Value of Velocity x-compomenent Field
  vys/                              Value of Velocity y-compomenent Field
  vzs/                              Value of Velocity z-compomenent Field
  intput.txt                        Input file

The sim will also copy the initial input.txt file into the output folder.

Output files are binaries and are name in the following way:

Take some field, say $\psi$, it will be exported as

Psi_000000.bin

on the 0th save

When the number of save_steps has been reached again, the next file will be called

Psi_000001.bin

So the number in the file name simply corresponds to (no. of steps at time of saving)/(save_step).

Binaries can be pretty easily read into Python using numpy, for example:

field_location = "output/Psis/Psi_000000.bin"

# W_padded, L_padded, H_padded = W + 2, L + 2, H + 2 to account for ghost cells
field = np.fromfile(field_location, dtype=np.float64).reshape((W_padded, L_padded, H_padded)

Running the Sim

To run the simulation it is also important you have the eigen library installed.

An example compilation is as follows:

export OMP_NUM_THREADS=4
g++-15 -std=c++17 -O3 -march=native -fopenmp -I./include -I/opt/homebrew/include/eigen3 -o simulation src/*.cpp
./simulation

It is important for speed to set the number of threads for OpenMP otherwise your computer will probably default to a single thread.

When the simulation is running, at the save_step interval, the sim will print out diagnostics information:

--------------------------------
On step: 0
OpenMP max threads: 4
Total psi: 2109
Total phi: 29791
Total phi psi: 2109
Total D is 3190
Total T is 2109
Total P is 0
Total v is 0
random div v is 0

This provides useful information while the sim is running. Total X just gives you the sum of that entire field which is quite useful for debugging. random div v is just computing the divergence of the v field at (13,13,13). This is arbitrary and can be changed. The point of this is to monitor that the divergence condition is being met.

At the beginning and ocassionally during the simulation running, after the step 0 diagnostics print, you will see the following:

on try 1000: abs_sum_difference: 3.32962
on try 2000: abs_sum_difference: 2.80918
on try 3000: abs_sum_difference: 2.37007
on try 4000: abs_sum_difference: 1.99957
on try 5000: abs_sum_difference: 1.68697
on try 6000: abs_sum_difference: 1.42321
on try 7000: abs_sum_difference: 1.20067
on try 8000: abs_sum_difference: 1.01291
on try 9000: abs_sum_difference: 0.854479
on try 10000: abs_sum_difference: 0.720809

This is showing the progress the pressure solver is making. It is simply the value $\Delta P_n$. Normally at the start the pressure solver will take longer. This is because the Pressure field has intialised to zero across all points. It should not normally print after this point unless it has gotten stuck.

Bibliography and Credits

The work for this code was completed for the Shendruk Research Group and was funded by the University of Edinburgh, School of Physics and Astronomy Career Development Scholarship.

The simulation gets its system of equations from the following paper:

Huang, Z., Lin, G., & Ardekani, A. M. (2022). A consistent and conservative phase-field model for thermo-gas-liquid-solid flows including liquid-solid phase change. Journal of Computational Physics, 449, 110795. https://doi.org/10.1016/j.jcp.2021.110795

With some modifications.