Why are there nonzero PK1 and Cauchy stress in The Y-direction when applying growth only in the X-direction?

General
1

In the below MWE, I have a unit square. I am applying an exponential growth only in the X-direction with the help of the growth tensor Fg = as_tensor([[exp(r*T),0.0],[0.0,1.0]]), where r and T are the rate of growth and time, respectively. I am implementing a hyperelastic model to solve \eta\frac{\partial \mathbf{u}}{\partial t}=S \nabla.\mathbf{T}, where \eta is the friction coefficient arising from interaction of the growing square with a substrate it is immersed in, \mathbf{u} is a displacement vector, \mathbf{T} is the first Piola-Kirchoff stress tensor, and S is a constant that has dimension of area.

I solve the variational form subject to the condition that the Y-displacement is 0 for all time steps (otherwise, it generates a Y-displacement and I want the square to grow only in the X-direction).

When I calculate the Cauchy and the PK1 stress, I find that the stress in the Y-direction is nonzero. I think it should be zero.

Here is the MWE:

from ufl import *
from ufl.classes import *
from dolfin import *     
import numpy as np


mesh = UnitSquareMesh(8,8)

ds_custom = Measure("ds", domain=mesh)   # Boundary measure                                                           
dx_custom = Measure("dx", domain=mesh)   # Surface measure


FS = FunctionSpace(mesh, 'P', 1)
VS = VectorFunctionSpace(mesh, 'P', 1) 
TS = TensorFunctionSpace(mesh, 'P', 1)

u = Function(VS) 
u_n = Function(VS) 
v = TestFunction(VS)

NumStepg = 5
Time = 2  # Units: mins
dt = Time/NumStepg
j = 0
eta = 5.0 # Friction coefficient 
r = 0.02  # Rate of growth
mu1 = 500.0/(1.0*(0.1**3)/12.0)  # Flexural rigidity  
nu1 = 0.499  # poisson ratio

csa3 = 0.01


vtk1 = File("displacement.pvd")
vtk2 = File("PK1stress.pvd")
vtk3 = File("CauchyStress.pvd")
vtk4 = File("Fe.pvd")
vtk5 = File("dW_dFe.pvd")

for T in np.linspace(0,Time,NumStepg): 

 d = mesh.geometry().dim()
 I = Identity(d)             # Identity tensor
 F = I + grad(u)             # Deformation gradient
 Fg = as_tensor([[exp(r*T),0.0],[0.0,1.0]])       # Growth tensor only in the X-direction
 Fg_inv = as_tensor([[exp(-r*T),0.0],[0.0,1.0]])  # Inverse of the growth tensor
 
 Fe = variable(F*Fg_inv)
 Ce = variable(Fe.T*Fe)
 Je  = det(Fe)
 
 # Invariants of deformation tensors
 I1 = tr(Ce)  

 W = ((mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-(2*nu1)))*(Je-1)**2 - mu1*(ln(Je)))
 
 TT = ((det(Fg)*diff(W,Fe)*(Fg_inv.T)))  #   1st PK Stress for compressible morphoelasticity
 
 PE = (eta)*dot(u,v)*dx_custom + (csa3)*(dt)*inner(TT,grad(v))*dx_custom - (eta)*dot(u_n, v)*dx_custom
 
 all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
 bcs2 = DirichletBC(VS.sub(1), Constant(0.0), all_domains, 1)              # making Y- component of displacement zero
 bcs = [bcs2]
 
  
 # Solve variational problem
 solve(PE == 0, u , bcs,
      solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,   
      "absolute_tolerance": 5e-10,"maximum_iterations": 50}})
 
 print(f"Number of iterations: {j}", flush=True)
 print(f"Time = {T} mins", flush=True)
 
 j += 1
 
 u_n.assign(u)
 
 cauchystress = (2/Je)*((mu1/2)*(Fe*Fe.T - I) + (nu1*mu1*(Je-1)*Je*I)/(1-2*nu1))   # Cartesian Cauchy Stress for compressible morphoelasticity
 
 PK1 = project(TT, TS)
 disp = project(u, VS)
 CS = project(cauchystress, TS)
 dwdfe = project(diff(W,Fe), TS)
 fe = project(Fe, TS)
 
 PK1.rename('PK1','PK1')
 disp.rename('disp','disp')
 CS.rename('CS','CS')
 dwdfe.rename('dwdfe','dwdfe')
 fe.rename('fe','fe')
 
 vtk1<<(disp, T)
 vtk2<<(PK1, T)
 vtk3<<(CS, T)
 vtk4<<(fe, T)
 vtk5<<(dwdfe, T)

Furthermore, when I calculate \frac{\partial W}{\partial \mathbf{F}_e} manually, it gives me zero values in the (1,1) position when \mathbf{F}_e has the (1,1) component 1 for all time points. Since W is a scalar function and \mathbf{F}_e has the value 1 in the (1,1) position for all time points, \frac{\partial W}{\partial \mathbf{F}_e} must have 0 at the (1,1) position since W can be rewritten as a function of the first (0,0) component of \mathbf{F}_e only.

Please let me know if there is anything wrong with the code or my understanding.
Thanks in advance!

2

I did not look at the details but if you need to constrain the displacement to zero in the Y direction to avoid growth in this direction then it is natural to have a non-zero stress in this direction as a reaction to prevent growth

3

Hi @bleyerj ,

I understand what you mentioned. The problem becomes more interesting when the same mesh is created in Gmsh and imported to FEniCS. In the below MWE, I am importing a square mesh created in Gmsh (script given after MWE) and solving the same problem with restricting displacement in the Y-direction but not in the Z- direction. In Paraview, we can see that the displacement is always 0 in the Z-direction (without even restricting it) but shows a stress in the Z-direction! I cannot figure out why such is the case.

MWE:

from dolfin import *     
import meshio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from ufl import cofac, sqrt
from scipy import linalg as la
from ufl import *
from ufl.classes import *
import pandas as pd

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=True):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})    
    return out_mesh
    
######################################################################################################################

msh = meshio.read("square.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)

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

FS = FunctionSpace(mesh, 'P', 1)
VS = VectorFunctionSpace(mesh, 'P', 1) 
TS = TensorFunctionSpace(mesh, 'P', 1)

u = Function(VS) 
u_n = Function(VS) 
v = TestFunction(VS)

NumStepg = 5
Time = 2  # Units: mins
dt = Time/NumStepg
j = 0
eta = 5.0 # Friction coefficient 
r = 0.02  # Rate of growth
mu1 = 500.0/(1.0*(0.1**3)/12.0)  # Flexural rigidity  
nu1 = 0.499  # poisson ratio

csa3 = 0.01


vtk1 = File("displacement.pvd")
vtk2 = File("PK1stress.pvd")
vtk3 = File("CauchyStress.pvd")
vtk4 = File("Fe.pvd")
vtk5 = File("dW_dFe.pvd")

for T in np.linspace(0,Time,NumStepg): 

 d = mesh.geometry().dim()
 I = Identity(d)             # Identity tensor
 F = I + grad(u)             # Deformation gradient
 Fg = as_tensor([[exp(r*T),0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])        # Growth tensor only in the X-direction
 Fg_inv =as_tensor([[exp(-r*T),0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])   # Inverse of the growth tensor
 
 Fe = variable(F*Fg_inv)
 Ce = variable(Fe.T*Fe)
 
 Je  = det(Fe)
 
 # Invariants of deformation tensors
 I1 = tr(Ce)  

 W = ((mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-(2*nu1)))*(Je-1)**2 - mu1*(ln(Je)))
 
 TT = ((det(Fg)*diff(W,Fe)*(Fg_inv.T)))  #   1st PK Stress for compressible morphoelasticity
 
 PE = (eta)*dot(u,v)*dx_custom(5) + (csa3)*(dt)*inner(TT,grad(v))*dx_custom(5) - (eta)*dot(u_n, v)*dx_custom(5)
 
 all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
 bcs2 = DirichletBC(VS.sub(1), Constant(0.0), all_domains, 1)              # making Y- component of displacement zero
 bcs = [bcs2]
 
  
 # Solve variational problem
 solve(PE == 0, u , bcs,
      solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,   
      "absolute_tolerance": 5e-10,"maximum_iterations": 50}})
 
 print(f"Number of iterations: {j}", flush=True)
 print(f"Time = {T} mins", flush=True)
 
 j += 1
 
 u_n.assign(u)
 
 cauchystress = (2/Je)*((mu1/2)*(Fe*Fe.T - I) + (nu1*mu1*(Je-1)*Je*I)/(1-2*nu1))   # Cartesian Cauchy Stress for compressible morphoelasticity
 
 PK1 = project(TT, TS)
 disp = project(u, VS)
 CS = project(cauchystress, TS)
 dwdfe = project(diff(W,Fe), TS)
 fe = project(Fe, TS)
 
 PK1.rename('PK1','PK1')
 disp.rename('disp','disp')
 CS.rename('CS','CS')
 dwdfe.rename('dwdfe','dwdfe')
 fe.rename('fe','fe')
 
 vtk1<<(disp, T)
 vtk2<<(PK1, T)
 vtk3<<(CS, T)
 vtk4<<(fe, T)
 vtk5<<(dwdfe, T)

Gmsh script:

SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Line(1) = {4, 3};
//+
Line(2) = {3, 2};
//+
Line(3) = {2, 1};
//+
Line(4) = {1, 4};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("top", 3) = {1};
//+
Physical Curve("bottom", 4) = {3};
//+
Physical Curve("left", 1) = {4};
//+
Physical Curve("right", 2) = {2};
//+
Physical Surface("body", 5) = {1};