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 .
Allen-Cahn (AC) and Cahn-Hilliard (CH) equations are solved in a square Ω = [ 0 , 50 ] × [ 0 , 50 ] \Omega=[0,50]\times[0,50] Ω = [ 0 , 50 ] × [ 0 , 50 ] using a partitioned algorithm (but the Allen-Cahn equations are tightly coupled).
( C H ) ∂ ρ ∂ t = ∇ ⋅ D ∇ ( ∂ f ∂ ρ − κ ρ ∇ 2 ρ ) ( A C ) ∂ η i ∂ t = − 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}
( C H ) ∂ t ∂ ρ ( A C ) ∂ t ∂ η i = ∇ ⋅ D ∇ ( ∂ ρ ∂ f − κ ρ ∇ 2 ρ ) = − L ( ∂ η i ∂ f − κ η ∇ 2 η i ) , i = 1 , 2
The mobility coefficient D D D in the Cahn-Hilliard equations is given by:
D = D v o l p ( ρ ) + D v a p [ 1 − p ( ρ ) ] + D s u r f ρ ( 1 − ρ ) + D g b ∑ i ∑ i ≠ m η i η m p ( ρ ) = ρ 3 ( 6 ρ 2 − 15 ρ + 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}
D p ( ρ ) = D vol p ( ρ ) + D vap [ 1 − p ( ρ ) ] + D surf ρ ( 1 − ρ ) + D gb i ∑ i = m ∑ η i η m = ρ 3 ( 6 ρ 2 − 15 ρ + 10 )
The free energy density is defined by:
F = ∫ V [ f ( ρ , η 1 … p ) + κ ρ 2 ( ∇ ρ ) 2 + ∑ i κ η 2 ( ∇ η i ) 2 ] d v
\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}
F = ∫ V [ f ( ρ , η 1 … p ) + 2 κ ρ ( ∇ ρ ) 2 + i ∑ 2 κ η ( ∇ η i ) 2 ] d v
where
f ( ρ , η i … p ) = 16 ρ 2 ( 1 − ρ ) 2 + [ ρ 2 + 6 ( 1 − ρ ) ∑ i η i 2 − 4 ( 2 − ρ ) ∑ i η i 3 + 3 ( ∑ i η i 2 ) 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}
f ( ρ , η i … p ) = 16 ρ 2 ( 1 − ρ ) 2 + ⎣ ⎡ ρ 2 + 6 ( 1 − ρ ) i ∑ η i 2 − 4 ( 2 − ρ ) i ∑ η i 3 + 3 ( i ∑ η i 2 ) 2 ⎦ ⎤
Initial condition
The initial condition for the Cahn-Hilliard equation is given by:
ρ ( x , y ) = { 1 , ( x − 25 ) 2 + ( y − 20 ) 2 < 100 , 1 , ( x − 25 ) 2 + ( y − 35 ) 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}
ρ ( x , y ) = ⎩ ⎨ ⎧ 1 , 1 , 0 , ( x − 25 ) 2 + ( y − 20 ) 2 < 100 , ( x − 25 ) 2 + ( y − 35 ) 2 < 25 , otherwise .
The initial conditions for the Allen-Cahn equations are given by:
η 1 ( x , y ) = { 1 , ( x − 25 ) 2 + ( y − 20 ) 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}
η 1 ( x , y ) = { 1 , 0 , ( x − 25 ) 2 + ( y − 20 ) 2 < 100 , otherwise .
η 2 ( x , y ) = { 1 , ( x − 25 ) 2 + ( y − 35 ) 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}
η 2 ( x , y ) = { 1 , 0 , ( x − 25 ) 2 + ( y − 35 ) 2 < 25 , otherwise .
Parameters used for the test
Description
Symbol
Value
mobility coefficient (AC)
L L L
10.0 10.0 10.0
energy gradient coefficient (AC)
κ η \kappa_\eta κ η
5 5 5
energy gradient coefficient (CH)
κ ρ \kappa_\rho κ ρ
2 2 2
bulk diffusivity in the lattice (CH)
D v o l D_{\mathrm{vol}} D vol
0.04 0.04 0.04
diffusivity of the vapor phase (CH)
D v a p D_{\mathrm{vap}} D vap
0.002 0.002 0.002
surface diffusivity (CH)
D s u r f D_{\mathrm{surf}} D surf
16 16 16
grain boundary diffusivity (CH)
D g b D_{\mathrm{gb}} D gb
1.6 1.6 1.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] t ∈ [ 0 , 2 ] with a time-step δ t = 1 0 − 4 \delta t=10^{-4} δ t = 1 0 − 4 .
Spatial discretization for convergence analysis: uniform grid with N = 200 N={200} N = 200 nodes in each spatial direction, with Q 1 \mathcal{Q}_1 Q 1 finite elements
Newton solver: relative tolerance 1 0 − 10 10^{-10} 1 0 − 10 , absolute tolerance 1 0 − 14 10^{-14} 1 0 − 14
Results
Figure 1 shows the evolution of conserved and non conserved phase-fields at time step 100 100 100 , 5000 5000 5000 , 12500 12500 12500 , and 20000 20000 20000 (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 (see Figure 4.12 in the reference).
Figure 1: evolution of conserved and non conserved phase-fields at time step 100, 5000, 12500, and 20000 (from top to bottom).
Figure 2: time evolution of two unequal sized particles.