Skip to content

Example 1: sintering of two unequal sized particles

Files

Statement of the problem

This test is a Finite Element version of the sintering test of two unequal sized particles from 1.

Allen-Cahn (AC) and Cahn-Hilliard (CH) equations are solved in a square Ω=[0,50]×[0,50]\Omega=[0,50]\times[0,50] using a partitioned algorithm (but the Allen-Cahn equations are tightly coupled).

(CH)ρt=D(fρκρ2ρ)(AC)ηit=L(fηiκη2ηi),i=1,2 \begin{align} (CH)\quad \frac{\partial \rho}{\partial t} &= \nabla \cdot D \nabla \left( \frac{\partial f}{\partial \rho} - \kappa_{\rho}\nabla^{2}\rho \right) \\[5pt] (AC)\quad \frac{\partial \eta_i}{\partial t} &= -L\left( \frac{\partial f}{\partial \eta_i} - \kappa_{\eta}\nabla^{2}\eta_i \right), \quad i=1,2 \end{align}

The mobility coefficient DD in the Cahn-Hilliard equations is given by:

D=Dvolp(ρ)+Dvap[1p(ρ)]+Dsurfρ(1ρ)+Dgbiimηiηmp(ρ)=ρ3(6ρ215ρ+10) \begin{align} D &= D_{\mathrm{vol}}p(\rho) + D_{\mathrm{vap}}\left[1-p(\rho)\right] + D_{\mathrm{surf}}\rho(1-\rho) + D_{\mathrm{gb}} \sum_i \sum_{i\neq m}\eta_i\eta_m \\[5pt] p(\rho)&=\rho^3(6\rho^2-15\rho +10) \end{align}

The free energy density is defined by:

F=V[f(ρ,η1p)+κρ2(ρ)2+iκη2(ηi)2]dv \begin{align} F &= \int_V \left[ f(\rho,\eta_{1\ldots p}) + \frac{\kappa_{\rho}}{2}(\nabla \rho)^2 + \sum_i \frac{\kappa_{\eta}}{2} (\nabla \eta_i)^2 \right] \, dv \end{align}

where

f(ρ,ηip)=16ρ2(1ρ)2+[ρ2+6(1ρ)iηi24(2ρ)iηi3+3(iηi2)2] \begin{align} f(\rho,\eta_{i\ldots p}) &= 16\rho^{2}(1-\rho)^{2} + \left[ \rho^{2} + 6(1-\rho)\sum_i \eta_i^{2} - 4(2-\rho)\sum_i \eta_i^{3} + 3\left(\sum_i \eta_i^{2}\right)^{2} \right] \end{align}

Initial condition

The initial condition for the Cahn-Hilliard equation is given by:

ρ(x,y)={1,(x25)2+(y20)2<100,1,(x25)2+(y35)2<25,0,otherwise. \begin{align} \rho(x,y)= \begin{cases} 1, & (x-25)^2+(y-20)^2 < 100, \\[4pt] 1, & (x-25)^2+(y-35)^2 < 25, \\[4pt] 0, & \text{otherwise}. \end{cases} \end{align}

The initial conditions for the Allen-Cahn equations are given by:

η1(x,y)={1,(x25)2+(y20)2<100,0,otherwise. \begin{align} \eta_1(x,y)= \begin{cases} 1, & (x-25)^2+(y-20)^2 < 100, \\[4pt] 0, & \text{otherwise}. \end{cases} \end{align}

η2(x,y)={1,(x25)2+(y35)2<25,0,otherwise. \begin{align} \eta_2(x,y)= \begin{cases} 1, & (x-25)^2+(y-35)^2 < 25, \\[4pt] 0, & \text{otherwise}. \end{cases} \end{align}

Parameters used for the test

Description Symbol Value
mobility coefficient (AC) LL 10.010.0
energy gradient coefficient (AC) κη\kappa_\eta 55
energy gradient coefficient (CH) κρ\kappa_\rho 22
bulk diffusivity in the lattice (CH) DvolD_{\mathrm{vol}} 0.040.04
diffusivity of the vapor phase (CH) DvapD_{\mathrm{vap}} 0.0020.002
surface diffusivity (CH) DsurfD_{\mathrm{surf}} 1616
grain boundary diffusivity (CH) DgbD_{\mathrm{gb}} 1.61.6

Boundary conditions

Neumann boundary conditions are prescribed on boundary of the domain.

Numerical scheme

  • Time integration: Euler Implicit over the interval t[0,2]t\in[0,2] with a time-step δt=104\delta t=10^{-4}.
  • Spatial discretization for convergence analysis: uniform grid with N=200N={200} nodes in each spatial direction, with Q1\mathcal{Q}_1 finite elements
  • Newton solver: relative tolerance 101010^{-10}, absolute tolerance 101410^{-14}

Results

Figure 1 shows the evolution of conserved and non conserved phase-fields at time step 100100, 50005000, 1250012500, and 2000020000 (from top to bottom). Figure 2 shows the time evolution of two unequal sized particles. The results are in good agreement with those presented in 1 (see Figure 4.12 in the reference).

sintering sintering sintering sintering
Figure 1: evolution of conserved and non conserved phase-fields at time step 100, 5000, 12500, and 20000 (from top to bottom).
sintering
Figure 2: time evolution of two unequal sized particles.

  1. S Bulent Biner and others. Programming phase-field modeling. Springer, 2017.