PDDO-MultiPhysics Examples Phase field and microstructure

Glass-ceramic crystallization

Crystallization of a lithium disilicate glass-ceramic during a heating ramp: 41 grain fields and a composition field, with nucleation handled by Python hooks.

crystal_growth.yaml 262,144 nodes imex-euler, 123 min

1The Problem

As a Li2O-SiO2 glass is heated from 550 to 850 °C, four crystalline phases nucleate and grow in sequence, each in its own temperature window: lithium metasilicate, cristobalite, lithium disilicate and β-quartz. Nucleation is stochastic, so whether a phase appears also depends on the random draws of the run. Each grain is an Allen-Cahn field with its phase's mobility, anisotropy and driving force; the composition evolves by Cahn-Hilliard diffusion coupled to all of them.

$$\frac{\partial\phi_g}{\partial t} = L_0 L_g\,a_g(\nabla\phi_g)\Big[\kappa\nabla^2\phi_g - W'(\phi_g) + \Delta G_g\,p'(\phi_g) - 2\beta\,\phi_g\sum_{k\ne g}\phi_k^2\Big], \quad g = 1\ldots 41$$
$$\frac{\partial c}{\partial t} = M\,\nabla^2\frac{\partial f}{\partial c} - M\kappa_c\,\nabla^4 c$$

2The Problem File

The 41 grain fields are one indexed field, phi: {type: scalar, index: g}, and the equations, relations and phase parameters use the index (m_loc[g], ph_ceq: {index: p, values: [...]}). The parts that are not PDEs, such as stochastic nucleation, temperature windows, interface noise and phase bookkeeping, live in a hooks file, crystal_growth_hooks.py, whose functions the engine calls before and after each step.

The complete problem file, examples/crystal_growth/crystal_growth.yaml 
# Glass-ceramic crystallization during an in-situ XRD heating ramp (Li2O-SiO2 glass):
# multigrain Allen–Cahn (one field per grain slot) + Cahn–Hilliard composition.
# Same model as crystal_growth_pddo/crystal_growth_insitu_xrd_3d_pddo.py (parabolic driving
# force); the non-PDE parts (nucleation, temperature-dependent growth windows, interface
# noise, clipping, mass correction, phase fractions) are in crystal_growth_hooks.py.
#
#   dt(phi_g) = L0 Lm_g a_g(grad phi_g) [kappa lap(phi_g) - W'(phi_g) + dG_g p'(phi_g)
#                                        - 2 beta phi_g (sum_k phi_k^2 - phi_g^2)]
#   dt(c)     = M lap(df/dc) - M kappa_c bih(c),   df/dc = 2 A_c (c sum phi - sum phi c_eq)
"k">name: crystal_growth
"k">description: Glass-ceramic crystallization, in-situ XRD heating ramp 550 -> 850 C
"k">domain:
  "k">size: [64, 64, 64]
  "k">n: [64, 64, 64]
  "k">periodic: [true, true, true]
"k">indices:
  "k">g: 41            # grain slots: 3 Li2SiO3 + 8 cristobalite + 24 LS2 + 6 beta-quartz
  "k">p: 4             # crystalline phases
"k">parameters:
  # --- phase-field model
  "k">kappa: 2.0
  "k">W: 1.0
  "k">L0: 1.4
  "k">beta: 3.0
  "k">A_c: 6.0
  "k">M_c: 1.0
  "k">kappa_c: 1.0
  "k">c_init: 0.46
  "k">c_noise: 0.02
  # --- heating schedule (deg C; t_end = 1800 steps * dt)
  "k">T_start: 550.0
  "k">T_end: 850.0
  # --- phases (index p): Li2SiO3 (meta), cristobalite, LS2 (Li2Si2O5), beta-quartz
  "k">ph_ceq:    {index: p, values: [0.67, 0.05, 0.50, 0.25]}
  "k">ph_delta:  {index: p, values: [0.40, 0.45, 0.95, 0.30]}
  "k">ph_cubic:  {index: p, values: [0, 1, 0, 0]}          # anisotropy: 1 cubic, 0 axial
  "k">ph_onset:  {index: p, values: [640, 660, 690, 780]}
  "k">ph_peak:   {index: p, values: [670, 685, 730, 800]}
  "k">ph_tmax:   {index: p, values: [750, 720, 850, 980]}
  "k">ph_diss:   {index: p, values: [0.50, 0.70, 0.00, 0.00]}
  "k">ph_Lm:     {index: p, values: [0.22, 0.90, 1.70, 1.00]}
  "k">ph_cap:    {index: p, values: [3, 8, 24, 6]}
  "k">ph_twin:   {index: p, values: [0, 1, 0, 0]}          # twin-forming phase (cristobalite)
  # --- twinning (reduced Clayton-Knap, at the maximum of the twin-forming phase; hooks)
  "k">tw_on: 1
  "k">tw_A: 1.0
  "k">tw_kpar: 0.03
  "k">tw_kperp: 3.0
  "k">tw_L: 1.0
  "k">tw_dt: 0.05
  "k">tw_steps: 200
  "k">tw_period: 4.5
  # --- nucleation and noise (used by the hooks)
  "k">nuc_interval: 25
  "k">nuc_attempts: 24
  "k">nuc_prob0: 0.5
  "k">nuc_barrier: 0.15
  "k">glass_thresh: 0.25
  "k">r_seed: 2.0
  "k">min_nuc_dist: 6.0
  "k">noise_amp: 0.008
  "k">rng_seed: 11
  # --- per-grain values, set by the hooks when a grain nucleates (runtime parameters)
  "k">T:      {runtime: true, value: 550.0}
  "k">Lm:     {index: g, runtime: true, value: 0.0}
  "k">ceq:    {index: g, runtime: true, value: 0.0}
  "k">gate:   {index: g, runtime: true, value: 0.0}        # growth window at the current T
  "k">d_ax:   {index: g, runtime: true, value: 0.0}
  "k">d_cub:  {index: g, runtime: true, value: 0.0}
  "k">R:      {index: g, runtime: true, shape: matrix, value: identity}   # grain orientation
"k">fields:
  "k">phi:   {type: scalar, index: g}
  "k">c: scalar
  "k">phase: {type: scalar, external: true}     # phase id + 1 where crystalline (output only)
  "k">grain: {type: scalar, external: true}     # owning grain + 1 where crystalline
"k">relations:
  "k">sumsq: {expr: "sum(phi[k]**2, k in g)", cache: true}
  "k">dfdc: "2*A_c*(c*sum(phi[k], k in g) - sum(phi[k]*ceq[k], k in g))"
  # anisotropy in the grain frame: m = R^T grad(phi) / |grad(phi)|
  "k">m_loc[g]: "dot(transpose(R[g]), grad(phi[g]))"
  "k">g2[g]: "dot(grad(phi[g]), grad(phi[g])) + 1e-12"
  "k">aniso[g]: >-
    max(1 + d_ax[g]*(3*m_loc[g][2]**2/g2[g] - 1)
          + d_cub[g]*((m_loc[g][0]**2/g2[g])**2 + (m_loc[g][1]**2/g2[g])**2
                      + (m_loc[g][2]**2/g2[g])**2 - 0.6), 0.05)
  "k">dG[g]: "gate[g]*max(1 - A_c*(c - ceq[g])**2, 0)"
"k">equations:
  "k">phi: >-
    dt(phi[g]) = L0*Lm[g]*aniso[g]*(kappa*lap(phi[g])
                 - 2*W*phi[g]*(1 - phi[g])*(1 - 2*phi[g])
                 + dG[g]*30*phi[g]**2*(1 - phi[g])**2
                 - 2*beta*phi[g]*(sumsq - phi[g]**2))
  "k">c: "dt(c) = M_c*lap(dfdc) - imp(M_c*kappa_c*bih(c))"
"k">initial_conditions:
  "k">phi: "0"
  "k">c: "c_init + c_noise*randn()"
"k">seed: 11
"k">time:
  "k">end: 36.0              # 1800 steps
  "k">dt: 0.02
  "k">scheme: imex-euler
"k">output:
  "k">every: 360
  "k">fields: [c, phase, grain]
"k">monitors:
  "k">T_C: "mean(T)"
  "k">mean_c: "mean(c)"
  "k">solid: "mean(sumsq)"
"k">hooks: crystal_growth_hooks.py
The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
Figure 1. The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
Indexed equations: one Allen-Cahn equation per grain slot g, and the composition.
Figure 2. Indexed equations: one Allen-Cahn equation per grain slot g, and the composition.
Phase parameters as indexed arrays: equilibrium composition, onset and peak temperatures, mobility.
Figure 3. Phase parameters as indexed arrays: equilibrium composition, onset and peak temperatures, mobility.
The Runtime Hooks node: the Python file with nucleation and the phase bookkeeping.
Figure 4. The Runtime Hooks node: the Python file with nucleation and the phase bookkeeping.

3Results

Nodes
64×64×64 = 262,144
Periodic
xyz
PDDO order
4, horizon 3.015 Δx
Time scheme
imex-euler
Steps
1,800 × Δt = 0.02
Newton iterations
1,800
Factorizations
0
Run time
123 min

Lithium metasilicate nucleates first, near 670 °C, and reaches 1.6 % of the volume at the stored times; lithium disilicate follows near 730 °C and grows at its expense, so by the end of the ramp the metasilicate has dissolved again. β-quartz appears near 790 °C. At 850 °C 87 % of the nodes are crystalline: 65 % lithium disilicate and 22 % β-quartz. Cristobalite does not nucleate in this run. The mean composition stays at 0.460 (the hooks correct it globally). The hooks also write these fractions to phase_fractions.csv at every output step.

Phase map at the end of the ramp: each colour is one crystalline phase, the glass between them.
Figure 5. Phase map at the end of the ramp: each colour is one crystalline phase, the glass between them.
Composition on three slices.
Figure 6. Composition on three slices.
The monitor solid = mean(Σφg2) against time: crystallization sets in after t ≈ 15 (about 700 °C).
Figure 7. The monitor solid = mean(Σφg2) against time: crystallization sets in after t ≈ 15 (about 700 °C).
The heating ramp.
Figure 8. The heating ramp.

4Run It Yourself

The problem file ships with PDDO-MultiPhysics. Open it in the graphical interface and press Compute, or run it from the command line:

python -m pddo_mp gui examples/crystal_growth/crystal_growth.yaml
python -m pddo_mp run examples/crystal_growth/crystal_growth.yaml

Results go to results/crystal_growth/ next to the problem file: a ParaView time series (.pvd), the monitors as CSV, and the final state. Any entry can be changed on the command line, for example --set time.end=18 for a run to half the end time. The user manual describes every part of the problem file and of the interface.

Have a coupled problem of your own?

Send us the physics and the geometry — we will write the problem file with you, or build the model you need.

Talk to an engineer About PDDO-MultiPhysics