![]() |
deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
|
This program was contributed by Umair Hussain <husain.umair2010@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
![]() |
This code solves the solidification problem based in the famous work by Ryo Kobayashi (1993) [1]. The model is based on the Allen-Cahn [2] phase field equation coupled with the transient heat equation Though we have covered only the isotropic directional solidification from the paper in the results, the same code can be modified and used for other types of solidification problems. Let us quickly go through the governing equations and the boundary conditions solved in this problem.
\begin{align*} \tau \frac{\partial p}{\partial t} = \nabla \cdot( \left(\epsilon^2\right)\nabla p) + p(1-p)(p-\frac{1}{2}+m) \label{heateq} \\ \frac{\partial T}{\partial t}=\nabla^2T+K\frac{\partial p}{\partial t} \end{align*}
where \(m(T) = \frac{a}{\pi}\tan^{-1}(\gamma(T_e-T))\)
The problem is subjected to the boundary conditions:
\begin{align*} p(0,y,t)= 1 \\ T(0,y,t)= T_\gamma -\Delta T \end{align*}
and the initial conditions:
\begin{align*} p(x,y,0)= 0 \\ T(x,y,0)= T_\gamma -\Delta T \end{align*}
Here, \(\Delta T\) is the degree of undercooling.
Using this code, we have reproduced one of the study from Kobayashi's work regarding the dendritic behaviour during directional solidification. The latent heat parameter 'K' in the equation determines the amount of heat released as the phase solidifies. If this value is high enough, we would observe an unstable interface between the solid and liquid phase, which would lead to the formation of dendrites as shown in these images. To assist this growth we need to add a random perturbation term ' \(a \chi p (1-p)\)' to the dynamic term of the phase field equation.
[1] Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3–4), 410–423. https://doi.org/10.1016/0167-2789(93)90120-P
[2] Allen, S. M., & Cahn, J. W. (1979). A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica, 27(6), 1085–1095. https://doi.org/10.1016/0001-6160(79)90196-2
For Parallel Computation
Initial values class
Prescribing p=1 at the left face (this will be maintained in the subsequent iterations when zero BC is applied in the Newton-Raphson iterations)
To apply the boundary values only to the solution vector without the Jacobian Matrix and RHS Vector
Separating each variable as a scalar to easily call the respective shape functions
To copy values and gradients of solution from previous iteration Old Newton iteration
Old time step iteration
Copying old solution values
Old solution values
Shape Functions
Shape Functions
Adding random noise at the interface
Assembling Jacobian matrix
Finding f1 and f2 at previous iteration for rhs vector
Adding random noise at the interface
Assembling RHS vector
Applying zero BC
Reading mesh
using only one process to output the result
Initialise the solution
Applying Boundary Conditions at t=0
Plotting initial solution
Time steps begin here:
Newton-Raphson iterations begin here:
Saving parallel vectors as non-parallel ones
Initialise the delta solution as zero
Assemble Jacobian and Residual
Solving to get delta solution
Checking for convergence
pcout << "Nothing wrong till here!!!!!!" << std::endl;
Transfer the converged solution to the old_solution vector to plot output
output the solution at only specific number of time steps
Using a direct parallel solver
Updating the solution by adding the delta solution