SB·Saswata BhattacharyyaIIT HYDERABAD
MS5033 · WORKED EXAMPLE · JULIA + NUMPY

Elasticity reshapes phase separation.

Use the elastic energy kernel from my teaching notebook in a phase-field calculation.

The physical idea.

This is coherent spinodal decomposition. A composition fluctuation produces the dilatational eigenstrain \(\boldsymbol\epsilon^0(c)=\epsilon_c(c-\bar c)\mathbf I\). Under homogeneous linear elasticity, mechanical equilibrium can be solved in Fourier space. Eliminating the displacement field gives the elastic-energy kernel \(B(\mathbf n)\), which adds an elastic penalty to the chemical and interfacial terms in the Cahn–Hilliard equation.

\[\widehat c^{\,n+1}=\frac{\widehat c^{\,n}-\Delta t\,M k^2\mathcal F[f^{\prime}(c^n)]}{1+\Delta t\,M k^2\!\left(\kappa k^2+B(\mathbf n)\right)}\]

This example is based on the elasticity section of my teaching notebook, PhaseField_Examples_Spectral.ipynb. I wrote that section to explain how elastic strain energy affects phase separation. Here, you can compare the composition patterns with and without elastic effects.

LIVE 256 × 256 SYSTEM · NATIVE-PIXEL DISPLAY

Change the physics.

Changing a value restarts the calculation from the same deterministic fluctuation. This makes comparisons meaningful.

How the elastic constants are constructed and scaled
\[\begin{aligned}\epsilon_c&=\frac{\text{misfit (\%)}}{100},\\ C_{44}^*&=\mu^*\frac{2A_Z}{1+A_Z},\\ C_{11}^*&=\mu^*\!\left[\frac{2(2+A_Z)}{1+A_Z}-\frac{1-4\nu}{1-2\nu}\right],\\ C_{12}^*&=\mu^*\!\left[\frac{2A_Z}{1+A_Z}-\frac{1-4\nu}{1-2\nu}\right].\end{aligned}\]

These are the relations used in PhaseField_Examples_Spectral.ipynb. With the notebook defaults \(\epsilon_c=1\%\), \(\mu^*=400\), \(\nu=1/3\), and \(A_Z=3\), they give \(C_{11}^*=1400\), \(C_{12}^*=1000\), and \(C_{44}^*=600\). Here \(\mu^*\) and \(C_{ij}^*\) are scaled, dimensionless energy coefficients, not values in Pa or GPa. Their scale is set relative to the chemical free-energy coefficient \(A=1\). Also, \(A_Z=2C_{44}^*/(C_{11}^*-C_{12}^*)\). Values outside the displayed stable range restore all defaults.

Time 0.0Mean \(c\) 0.50000Range 0.490 to 0.510

The solver and canvas are both 256 × 256. One simulation point maps to one displayed pixel. Blue: \(c=0\). Green: \(c=0.5\). Red: \(c=1\). The domain is periodic.

Model and numerical choices.

The interactive calculation uses the values in my elastochemical teaching notebook: \(A=1\), \(\kappa=1\), \(M=1\), \(\epsilon_c=0.01\) (1%), \(\mu^*=400\), \(\nu=1/3\), \(A_Z=3\), \(\Delta x=1\), \(\Delta t=0.2\), and \(t_{\max}=1000\). It evolves and displays a periodic 256 × 256 system at one canvas pixel per simulation point. The downloadable scripts below are an earlier simplified teaching adaptation with normalized elastic constants; their numbers should not be read as Pa or GPa. The real discrete elastic kernel is symmetrized across conjugate Fourier partners on Nyquist lines. The mean-stress-free zero mode has \(B(\mathbf 0)=0\). This is a dimensionless 2D constitutive example, not a full 3D plane-strain eigenstrain model.

Derive the kernel before coding it.

\[\begin{aligned}Q_{11}&=C_{11}n_1^2+C_{44}n_2^2\\Q_{22}&=C_{44}n_1^2+C_{11}n_2^2\\Q_{12}&=(C_{12}+C_{44})n_1n_2\\B(\mathbf n)&=\mathbf E:\mathbf C:\mathbf E-(\boldsymbol\sigma^0\mathbf n)\cdot\mathbf Q^{-1}\cdot(\boldsymbol\sigma^0\mathbf n)\end{aligned}\]

Here σ⁰ = C:E. The diagonal shear coefficient is C₄₄. The supplied notebook cell used C₁₂ in these positions; this adaptation corrects that contraction. For isotropic stiffness, this kernel is independent of direction. Positive elastic energy penalizes composition fluctuations.

To isolate elastic effects, compare chemical and elastic runs with the same 256 × 256 grid, initial condition, time and color scale. The default chemical gallery uses a different grid and is not a controlled side-by-side experiment.

Questions explored in this example.

Set C44 = (C11−C12)/2 and check that B becomes direction independent. Set εc = 0 to recover the chemical solver on the same 256 × 256 grid and initial field. Change εc and inspect suppression of fluctuations.

Bitter–Crum & Eshelby.

Understand the assumptions behind shape-independent elastic energy, uniform strain inside ellipsoidal inclusions, and the anisotropic patterns in the worked example.

Explore the theorems and their limits →
FOURIER MICROELASTICITY

Where the elastic kernel and driving forces come from.

Follow the derivation of the Khachaturyan elastic-energy kernel, then obtain the composition and order-parameter derivatives, including the real-space chain rule.

Read the derivations →

When matrix and precipitate have different stiffness.

Go beyond the homogeneous kernel with two explained equilibrium solvers and the corrected classroom Elastic.ipynb.

Explore direct iteration and equivalent eigenstrain →

My Python calculation, explained in steps.

Both versions implement the same model. Python uses NumPy; Julia uses FFTW for spectral examples. Scripts write CSV snapshots and diagnostics into the output folder.

Download the Python notebook ↓ · Read the code in steps → Download Julia version ↓

Run the notebook cells in order for the explained calculation. The commands and complete script below are an optional command-line route; the script is linked at the end of the notebook lesson.

python -m pip install numpy
python phase_field.py elastic results
Read the complete NumPy code
"""Phase-field example based on my teaching notebook PhaseField_Examples_Spectral.ipynb.
Purpose: explain spinodal decomposition and the effect of elastic strain energy.
Run: python phase_field.py [chemical|elastic] [output_directory]
Requires NumPy. CSV rows index x, columns index y; all units dimensionless.
Periodic 2D homogeneous elasticity; zero mean stress. Not a calibrated alloy.
The explicit bulk term means time-step convergence must be checked.
"""
import sys
from pathlib import Path
import numpy as np


def kernel(kx, ky, c11=4., c12=2., c44=3., eigenstrain=.2):
    """B = E:C:E - (sigma0.n).Q^-1.(sigma0.n), E=eigenstrain*I.
    Q_ik=C_ijkl*n_j*n_l; Q11=C11*nx²+C44*ny² (not C12).
    This is a 2D constitutive model, not a reduction with epsilon_zz eigenstrain.
    """
    k2 = kx*kx + ky*ky
    radius = np.sqrt(np.where(k2 == 0, 1., k2))
    nx, ny = kx/radius, ky/radius
    q11 = c11*nx*nx + c44*ny*ny
    q22 = c44*nx*nx + c11*ny*ny
    q12 = (c12+c44)*nx*ny
    det = q11*q22-q12*q12
    det = np.where(k2 == 0, 1., det)
    stress = (c11+c12)*eigenstrain
    b = 2*(c11+c12)*eigenstrain**2-stress**2*(q22*nx*nx-2*q12*nx*ny+q11*ny*ny)/det
    return np.where(k2 == 0, 0., b)


def run(mode='chemical', output='phase_output', n=None, dt=.2, steps=2000):
    out=Path(output);out.mkdir(parents=True,exist_ok=True)
    if n is None: n = 256 if mode == 'elastic' else 64
    if mode not in ('chemical','elastic'): raise ValueError('Choose chemical or elastic')
    k=2*np.pi*np.fft.fftfreq(n);kx,ky=np.meshgrid(k,k,indexing='ij');k2=kx*kx+ky*ky
    b=kernel(kx,ky) if mode == 'elastic' else np.zeros_like(k2)
    # Real-field discrete operator: symmetrize conjugate partners on even grids.
    # Only Nyquist lines can differ from the continuum directional kernel.
    partner=(-np.arange(n)) % n
    b=.5*(b+b[np.ix_(partner,partner)])
    i,j=np.meshgrid(np.arange(n),np.arange(n),indexing='ij')
    # Deterministic broadband perturbation, identical across Julia and NumPy.
    noise=np.sin((i+1)*12.9898+(j+1)*78.233)*43758.5453
    noise=.02*(noise-np.floor(noise)-.5);noise-=noise.mean();c=.5+noise
    initial_mean=c.mean();records=[]
    denominator=1+dt*k2*(k2+b)
    def record(step):
        h=np.fft.fft2(c)
        energy=np.mean(c*c*(1-c)**2)+.5*np.sum((k2+b)*np.abs(h)**2)/n**4
        records.append([step*dt,c.mean(),energy,c.min(),c.max()])
        np.savetxt(out/f'field_{step:04d}.csv',c,delimiter=',')
    record(0)
    for step in range(1,steps+1):
        h=(np.fft.fft2(c)-dt*k2*np.fft.fft2(2*c*(1-c)*(1-2*c)))/denominator
        c=np.fft.ifft2(h).real
        if step%25 == 0 or step==steps: record(step)
    assert np.isfinite(c).all()
    assert abs(c.mean()-initial_mean)<1e-10
    np.savetxt(out/'diagnostics.csv',records,delimiter=',',header='time,mean,energy_density,min,max',comments='')
    print(mode,'mean drift',abs(c.mean()-initial_mean),'energy',records[0][2],records[-1][2])
    return c,np.array(records)

if __name__=='__main__': run(*(sys.argv[1:3]))

About these calculations.

I use these examples to explain numerical methods. I provide the parameters, boundary conditions and numerical checks with the results. These are teaching calculations, not predictions calibrated for a particular material.

← All worked examples

How I develop my codes.

I start with a small calculation, check its result and then add the next part. These lessons follow the same sequence, with explanations and checks at each stage. Open the worked Python, Julia and C lessons →