Hi All,
Hope you are safe. I am solving a coupled phase field and elasticity equations. The weak form of these equations are given below:
where u is the displacement, and eta is the order parameter. The problem needs to be solved in nanoscale. The code for solving this problem is:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
L = 40. # in nm
H = 40.
Nx = 300
Ny = 300
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
# GL constants
thta = Constant(pi/6)
gamma_0 = Constant(0.1295)
AA = Constant(1.4025) # in GPa
Dt = 1 # in ps
t = 0
tMax = 1000
lmbda = Constant(24.0) # GPa
mu = Constant(19.4) # GPa
kappa = Constant(0.0878) # nJ/m
s = as_vector(( cos(thta), sin(thta)))
m = as_vector((-sin(thta), cos(thta)))
SM = outer(s,m)
P = sym(SM)
BoundaryU = Expression(('gam*x[1] + 0.001*t'), gam = 0.15, t = 0, degree = 1)
# Function Spaces
V_u = VectorElement('CG', mesh.ufl_cell(), 1) # displacement finite element
V_eta = FiniteElement('CG', mesh.ufl_cell(), 1) # damage finite element
V = FunctionSpace(mesh, MixedElement([V_u, V_eta]))
# Mechanical Boundary Conditions
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 40.0)
top = CompiledSubDomain("near(x[1], side) && on_boundary", side = 40.0)
bot= CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.0)
bc = []
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
top.mark(facets, 1)
bot.mark(facets, 2)
left.mark(facets, 3)
right.mark(facets, 4)
ds = Measure("ds", domain = mesh, subdomain_data=facets)
bc.append(DirichletBC(V.sub(0).sub(1), BoundaryU, facets, 1))
bc.append(DirichletBC(V.sub(0).sub(1), BoundaryU, facets, 2))
bc.append(DirichletBC(V.sub(0).sub(0), Constant(0.0), facets, 1))
bc.append(DirichletBC(V.sub(0).sub(0), Constant(0.0), facets, 2))
bc.append(DirichletBC(V.sub(1), Constant(0.0), bot))
bc.append(DirichletBC(V.sub(1), Constant(0.0), top))
##########################################################
U_ = TestFunction(V)
(u_, eta_) = split(U_)
dU = Function(V)
(du, deta) = split(dU)
Uold = Function(V)
(uold, eta_old) = split(Uold)
eps = sym(grad(du))
eps_pl = gamma_0*(deta**2*(3-2*deta))*sym(SM)
eps_el = eps - eps_pl
sigma = lmbda*tr(eps)*Identity(2) + 2.0*mu*eps_el
df0_detta1 = 2*AA*deta - 6*AA*deta**2 + 4*AA*deta**3
df0_detta2 = -2*inner(grad(du), sym(SM))
df0_detta3 = (3*deta**2 - 2*deta**3)*gamma_0
df0_detta4 = 6*deta*(1-deta)
df0_detta = df0_detta1 + mu*gamma_0*(df0_detta3+df0_detta2)*df0_detta4
dF = (df0_detta*eta_ + 2.0*kappa*dot(grad(deta), grad(eta_)))*dx
mech_form = inner(sigma, grad(u_))*dx
res = dF + mech_form
U_init = Expression(('0.', '0.', '(x[0] - 20)*(x[0] - 20) + (x[1] - 20)*(x[1] - 20) < 3*3 - tol ? 1.0 : 0'), tol = 1e-3, degree = 2)
Uold.assign(project(U_init,V))
dU.assign(Uold)
etaFile = File("results5/eta.pvd")
uFile = File("results5/u.pvd")
sigma_file = File("results5/sig.pvd")
while t < tMax:
uFile << dU.split()[0]
etaFile << dU.split()[1]
BoundaryU.t = t
solve(res==0, dU, bc, solver_parameters = {"newton_solver":{"linear_solver": \
"gmres", "preconditioner": "ilu" , "relative_tolerance": 1e-6, "absolute_tolerance": 1e-14} }, \
form_compiler_parameters = {"cpp_optimize": True, \
"quadrature_degree": 2})
sigma_project = project(sigma, TensorFunctionSpace(mesh, 'P', 1))
sigma_file << sigma_project
Uold.assign(dU)
t += Dt
It should be noted that a circular region of initial radius 3nm embedded in a rectangular domain as an initial condition for the order parameter. For the boundary condition, a simple tension with Neuman (free) conditions on traction and order parameter along the lateral edges of the domain was considered as:
The deformation is applied in small increments of magnitude 0.001. Unfortunately, the problem is not converged, and I actually do not know which part I am doing mistakes. I really do appreciate any help regarding this issue.
Bests,
Ben