# Phase field for Sintering of Two Particles

Hi everyone,

I am trying to reproduce the results of Paper 1 or Paper 2 for sintering of two particles. For simplicity, M (in equation 6 or 23) is considered a scalar, and no rigid body motion is considered. the code seems to be running alright, however, I am not able to get the order parameters, phi’s, to look like those of Figures 3 and 2 in from the two papers respectively. I am trying to get phi’s output simultaneously as the microstructure evolves. I have attached my current “working” code, any pointers on what the issue might be will be appreciated. Thank you.
P.S. This is my first question on this group, so please pardon me if I flaunted any group guidelines on asking questions. Below is my code;

``````import matplotlib.pyplot as plt
import numpy as np

# Parameters
T = 2.0            # final time
num_steps = 100    # number of time steps
dt = T / num_steps # time step size
kappa_c = 1.0        # parameter in the Cahn-Hilliard equation
kappa_phi = 0.5
theta = 0.50
A = 16.0
B = 1.0
L = 10.0            # mobility for Allen-Cahn
M = 1.0            # mobility for Cahn-Hilliard
num_particles = 2  # Number of order parameters / particles
lambda_interface = 0.005  # Interface width for smooth transition

# Create mesh and define function space
nx = ny = 100
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 2)  # Using P2 elements for higher accuracy

# Create a mixed function space for multiple phi_i, c, and mu
elements = [V.ufl_element()] * (num_particles + 2)  # One element for c, one for mu, and one for each phi_i
W = FunctionSpace(mesh, MixedElement(elements))

# Define initial conditions for each phi_i and c with smooth transition
class InitialConditionPhi(UserExpression):
def __init__(self, center, radius, lambda_interface, **kwargs):
super().__init__(**kwargs)
self.center = center
self.lambda_interface = lambda_interface

def eval(self, values, x):
r = np.sqrt((x[0] - self.center[0])**2 + (x[1] - self.center[1])**2)
# Smooth transition using arctan function
values[0] = 0.5 * (1 - np.arctan((self.radius - r) / self.lambda_interface) / np.pi)

def value_shape(self):
return ()

# Initial conditions for concentration field c
class InitialConditionC(UserExpression):
def eval(self, values, x):
r1 = np.sqrt((x[0] - 0.3)**2 + (x[1] - 0.5)**2)
r2 = np.sqrt((x[0] - 0.7)**2 + (x[1] - 0.5)**2)
values[0] = 1.0 if r1 < 0.2 or r2 < 0.2 else 0.0

def value_shape(self):
return ()

# Set up initial conditions
phi_init = [InitialConditionPhi(center=(0.3, 0.5), radius=0.2, lambda_interface=lambda_interface, degree=2),
c_init = InitialConditionC(degree=2)

# Interpolate initial conditions into the mixed function space
u_n = Function(W)
for i in range(num_particles):
assign(u_n.sub(i), interpolate(phi_init[i], V))
assign(u_n.sub(num_particles), interpolate(c_init, V))
assign(u_n.sub(num_particles + 1), interpolate(c_init, V))  # Initial condition for mu can be derived from c

# Define boundary condition (no-flux boundary)
def boundary(x, on_boundary):
return on_boundary

bcs = [DirichletBC(W.sub(i), Constant(0.0), boundary) for i in range(num_particles + 2)]

# Define the variational problem
u = Function(W)
v = TestFunction(W)
phis = split(u)[:num_particles]  # All phi_i components
c = split(u)[num_particles]      # Concentration component
mu = split(u)[num_particles + 1]  # Auxiliary variable component
phis_n = split(u_n)[:num_particles]
c_n = split(u_n)[num_particles]
mu_n = split(u_n)[num_particles + 1]

###################################################################################################################################################

# NEW
# Step-by-step calculation
phi_cubed_sum = sum([phi**3 for phi in phis])  # Summation over i (phi_i^3)
phi_squared_sum = sum([phi**2 for phi in phis])  # Summation over i (phi_i^2)

f1 = A * (c**2) * (1 - c)**2
f2 = 6 * (1 - c) * phi_squared_sum
f3 = 4 * (2 - c) * phi_cubed_sum
f4 = 3 * (phi_squared_sum)**2

f_bulk = f1 + B * (c**2 + f2 - f3 + f4)  # Full expression

# Bulk free energy derivatives
dfdc = diff(f_bulk, variable(c))
df_dphi = [diff(f_bulk, variable(phi)) for phi in phis]

# Define mu_mid using the theta-method
mu_mid = (1 - theta) * mu_n + theta * mu  # mu_n & mu are the previous and current values of chemical potential

# Weak froms
F_c = (c - c_n) * v[num_particles] * dx + dt * M * (dot(grad(v[num_particles]), grad(mu_mid))) * dx
F_mu = mu * v[num_particles+1] * dx - dfdc * v[num_particles+1] * dx - (dot(grad(v[num_particles+1]), (kappa_c * grad(c)))) * dx
F_phis = [(phi - phi_n) * v[i] * dx + dt * L * (df_dphi[i] * v[i] * dx + kappa_phi * dot(grad(phi), grad(v[i])) * dx) for i, (phi, phi_n) in enumerate(zip(phis, phis_n))]

# Combine all weak forms
F = F_c + F_mu + sum(F_phis)

###################################################################################################################################################

# Define the Jacobian
J = derivative(F, u)

# Create VTK files for saving output
vtkfile_phi = [File(f'output/phi_{i}.pvd') for i in range(num_particles)]
vtkfile_c = File('output/c.pvd')
vtkfile_mu = File('output/mu.pvd')

# Time-stepping
t = 0
for n in range(num_steps):
t += dt
solve(F == 0, u, bcs, J=J)
u_n.assign(u)

# Extract solutions for each phi_i, c, and mu
phis = u.split()[:num_particles]
c = u.split()[num_particles]
mu = u.split()[num_particles + 1]

# Save to file (VTK)
for i, phi in enumerate(phis):
vtkfile_phi[i] << (phi, t)
vtkfile_c << (c, t)
vtkfile_mu << (mu, t)

# Plot solution at some time steps (optional)
if n % 10 == 0:
plt.figure()
plot(c)
plt.title(f'C at time step {n}')
plt.show()

for i, phi in enumerate(phis):
plt.figure()
plot(phi)
plt.title(f'Phi_{i} at time step {n}')
plt.show()

plt.figure()
plot(mu)
plt.title(f'Mu at time step {n}')
plt.show()`````````