[Phase Field] Coupled Vectorized Order Parameters

variational formulation
1

Hey, hoping someone smarter than me can help me answer the following problem I am experiencing in a project I am working.

At present, I am working to solve a phase-field (PF) model (based on the Allen-Cahn/Time-Dependent Ginzburg-Landau equations) for ferroelectric domain structures (see e.g. here or here ). However, I think this question can be helpful to others who have vectorized order parameters as most of the PF models I have seen here are for scalar order parameters (crack propagation, Cahn-Hilliard, etc.)

from fenics import *
import itertools

def T4th(v1111, v1122, v1212):
    return as_tensor([[[[v1111, 0.0, 0.0], [0.0, v1122, 0.0], [0.0, 0.0, v1122]],
                       [[0.0, v1212, 0.0], [0.0, 0.0, 0.0],     [0.0, 0.0, 0.0]],
                       [[0.0, 0.0, v1212], [0.0, 0.0, 0.0],     [0.0, 0.0, 0.0]]],
                      [[[0.0, 0.0, 0.0],     [v1212, 0.0, 0.0], [0.0, 0.0, 0.0]],
                       [[v1122, 0.0, 0.0], [0.0, v1111, 0.0], [0.0, 0.0, v1122]],
                       [[0.0, 0.0, 0.0],     [0.0, 0.0, v1212], [0.0, 0.0, 0.0]]],
                      [[[0.0, 0.0, 0.0],     [0.0, 0.0, 0.0],     [v1212, 0.0, 0.0]],
                       [[0.0, 0.0, 0.0],     [0.0, 0.0, 0.0],     [0.0, v1212, 0.0]],
                       [[v1122, 0.0, 0.0], [0.0, v1122, 0.0], [0.0, 0.0, v1111]]]])


G1111 = 0.6*0.98*1E-10; G1122 = 0.0/G1111; G1212 = 0.6*0.98*1E-10/G1111
G1111 = 1
Gijkl = T4th(G1111, G1122, G1212)

mesh = BoxMesh(Point(0, 0, 0), Point(10e-9,10e-9,10e-9), 8,8,8)
V = VectorFunctionSpace(mesh, 'P', 1,dim=3)
P, P0 = Function(V), Function(V)
P0.vector()[:] = 0.1 * (0.5 - np.random.rand(*P0.vector()[:].shape))
dt = 1e-2
dP, v = TrialFunction(V), TestFunction(V)

a0 = abs(3.8*(25-479)*1e5)
p0=0.757
a1,a11,a12,a111,a112,a123 = 3.8*(25-479)*1e5/a0, -7.3e7 * p0**2 / a0, 7.5e8 * p0**2 / a0, 2.6e8 * p0**4 / a0, 6.1e8 * p0**4 / a0, -3.7e9 * p0**4 / a0

Landau = 0
for i in range(3):
  Landau += P[i]**2 * a1 + P[i]**4 * a11 + P[i]**6 * a111
  for j in range(i,3):
    Landau += P[i]**2 * P[j]**2 * a12 + (P[i]**4 * P[j]**2 + P[i]**2 * P[j]**4) * a112
  
Landau += P[0]**2 * P[1]**2 * P[2]**2 * a123

Grad = 0
for i,j,k,l in itertools.product(range(3),range(3),range(3),range(3)):
  Grad += 0.5 * Gijkl[i,j,k,l] * P[i].dx(j) * P[k].dx(l)

F = (Landau + Grad)
dFdP = diff(F,P)

At this point, is where I get stumped. I have both P and P0 defined (P represents the n+1 timestep I am solving for and I want to do some implicitly, and P0 represents the n timestep I would currently be at). The variable dFdP represents the driving force (\frac{\delta F}{\delta P_i}) in the equation below.

\frac{\partial P_i}{\partial t} = -\frac{\delta F}{\delta P_i}

Which using backward Euler notation can be written as

P_i^{n+1} - P_i^{n} = -\Delta t \frac{\delta F}{\delta P_i^{n+1}}

I have tried to solve this a few different ways such as

fun = P*v*dx - P0*v*dx + dFdP*v*dx
solve(fun==0,P)

But the answer is not correct (note I have a spectral solver for these equations already but need an FEM solver for more complex boundary conditions).

Any help on the proper way to set up the problem to be solved would be very much appreciated. Open to any questions along the way as well.

Cheers,
PF

1 Like
2

Found my issue. All of my variables where renormalized to my length scale, but I hadn’t properly changed simulation box length (i.e. Point(10e-9,10e-9,10e-9) should be Point(1,1,1)).

After this, the problem solves correctly.

Thanks and hopes this helps someone.