Hi everyone,
In the below MWE, I have a three-layered annulus, with the radii being 3, 3.5, 4, and 4.5. I am defining the incompressibility factor to be \nu=0.49, where \nu=0.5 defines perfect incompressibility. I am minimizing an energy functional of the form \Pi[\mathbf{u}]=\sum_{i} \int_{\Omega_{i}}J_{g,i}W(\mathbf{F}_{e}), where i=1,2,3 stands for each of the layers, \Omega_{i} is the reference stress free configuration for each of the layers, \mathbf{F} = \mathbf{I}+\nabla\mathbf{u}, \mathbf{u} is the displacement vector. In the above equation, we are assuming that \mathbf{F}=\mathbf{F}_{e}\mathbf{F}_{g}, where \mathbf{F}_{g}=\text{diag}(1+t, 1+t, 1) for the innermost layer and Identities for the other two layers, t stands for time, \mathbf{F}_{e} registers the elastic deformation due to growth, and J_{g,i} is the determinant of the growth tensor \mathbf{F}_{g,i}. I am subjecting only the innermost layer to growth.
If the area of the innermost layer is A at time T=0, then after T=t, the area will be close to A(1+t)^2. But I am getting a substantial deficit in the area. If I start with A=10.210239232039262, then after T=0.05, the area should be close to A(1.05)^2=11.256788753323287 but from the simulation, I am getting 10.210429082293246, which is way off considering the very fine mesh I am working with.
Any suggestions would be extremely helpful. If you have any questions, please let me know. Below are the MWE, followed by a .geo script I am working with. If you think there is a mistake in the coding, please let me know. It will help me figure out whether it’s a coding error or there is some mathematical modelling required for compressible cases. I have a feeling that it’s the latter.
MWE:
msh = meshio.read("mymesh.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("dS", domain=mesh, subdomain_data=mf)
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)
#************************************************************************************
V_ima = VectorFunctionSpace(mesh, 'P', 1) # Changed to second order
W_ima_ref = FunctionSpace(mesh, 'P', 1)
T_ima = TensorFunctionSpace(mesh, 'P', 1)
initial_int_area = assemble(Constant(1.0)*dx_custom(5))
initial_med_area = assemble(Constant(1.0)*dx_custom(6))
initial_ad_area = assemble(Constant(1.0)*dx_custom(7))
print('initial intima area =', assemble(Constant(1.0)*dx_custom(5)))
print('initial media area =', assemble(Constant(1.0)*dx_custom(6)))
print('initial adventitia area =', assemble(Constant(1.0)*dx_custom(7)))
#***********************************************************************************************
ug1 = Function(V_ima)
vg1 = TestFunction(V_ima)
# Elasticity parameters
nu = 0.4999
mu1, mu2, mu3 = 27.9, 1.27, 7.56
ds = 0.001
NumStepg = 51
for j in range(0,NumStepg):
tj = j*ds
dg = mesh.geometry().dim()
Ig = Identity(dg) # Identity tensor
Fg = Ig + grad(ug1) # Deformation gradient
g_1 = Expression((('1.0+p1*t1','0.0','0.0'),
('0.0','1.0+p2*t1','0.0'),
('0.0','0.0','1.0')),degree=1,p1=1.0,p2=1.0,t1=tj)
ginv = Expression((('1.0/(1.0+p1*t1)','0.0','0.0'),
('0.0','1.0/(1.0+p2*t1)','0.0'),
('0.0','0.0','1.0')),degree=1,p1=1.0,p2=1.0,t1=tj)
g_i = project(det(g_1),W_ima_ref)
Ce = (Fg*ginv).T*(Fg*ginv) # Right Cauchy-Green tensor for the Intima
Fe = Fg*ginv
Mg = cofac(Fg)
# Invariants of deformation tensors
I1 = tr(Ce) # Ic for the Intima
Jg_e = det(Fg)
psig1 = (mu1/2)*(I1 - 3) + ((nu)*(mu1)/(1-2*nu))*(Jg_e-1)**2 - mu1*(ln(Jg_e))
psig2 = (mu2/2)*(I1 - 3) + ((nu)*(mu2)/(1-2*nu))*(Jg_e-1)**2 - mu2*(ln(Jg_e))
psig3 = (mu3/2)*(I1 - 3) + ((nu)*(mu3)/(1-2*nu))*(Jg_e-1)**2 - mu3*(ln(Jg_e))
Pi = (det(g_1))*psig1*dx_custom(5) + psig2*dx_custom(6) + psig3*dx_custom(7)
# Compute first variation of Pi (directional derivative about u in the direction of v)
Fg1 = derivative(Pi, ug1, vg1)
bcgs_ima = []
ug1.set_allow_extrapolation(True)
# Solving variational problem
print(f"Solving for IMA displacement field at time T={tj}", flush=True)
solve(Fg1 == 0, ug1, bcgs_ima, solver_parameters={"newton_solver": {"relative_tolerance": 1e-8,
"absolute_tolerance": 1e-8,"maximum_iterations": 200}})
ug1.set_allow_extrapolation(True)
##########################################################################################
ug_ima = Function(V_ima)
vg1 = TestFunction(V_ima)
ug1.set_allow_extrapolation(True)
ug_ima.assign(ug1)
ug1 = interpolate(ug_ima, V_ima)
ug1.set_allow_extrapolation(True)
print('end of the growth loop ',tj, flush=True)
print('T =', tj, flush=True)
ALE.move(mesh, ug1)
print('initial intima area =', initial_int_area)
print('final intima area =', assemble(Constant(1.0)*dx_custom(5)))
print('initial media area =', initial_med_area)
print('final media area =', assemble(Constant(1.0)*dx_custom(6)))
print('initial adventitia area =', initial_ad_area)
print('final adventitia area =', assemble(Constant(1.0)*dx_custom(7)))
.geo script (Global mesh size factor: 0.1, Frontal Delaunay):
Circle(1) = {0.0, 0.0, 0.0, 3, 0, 2*Pi};
//+
Circle(2) = {0.0, 0.0, 0.0, 3.5, 0, 2*Pi};
//+
Circle(3) = {0.0, 0.0, 0.0, 4.0, 0, 2*Pi};
//+
Circle(4) = {0.0, 0.0, 0.0, 4.5, 0, 2*Pi};
//+
Curve Loop(1) = {1};
//+
Curve Loop(2) = {2};
//+
Curve Loop(3) = {2};
//+
Curve Loop(4) = {1};
//+
Plane Surface(1) = {3, 4};
//+
Curve Loop(5) = {3};
//+
Curve Loop(6) = {2};
//+
Plane Surface(2) = {5, 6};
//+
Curve Loop(7) = {4};
//+
Curve Loop(8) = {3};
//+
Plane Surface(3) = {7, 8};
//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {3};
//+
Physical Curve("C4", 4) = {4};
//+
Physical Surface("S1", 5) = {1};
//+
Physical Surface("S2", 6) = {2};
//+
Physical Surface("S3", 7) = {3};