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

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
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.

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

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

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

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

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")
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
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
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};