Hello, I am implementing the following matrix-fiber model (Liu et al., 2019):
,
where Ie and I4e contain the internal variables. To validate the code I am comparing with a numerical example, from the same paper, that deals with a relaxation simulation on a biaxially stretched plate. The plate is stretched in two directions for 1 second until it reaches a stretch of 1.1, then the displacement remains constant.
My problem is that after I have exhausted all attempts, the results of the Cauchy stress are different, using the analytical expression and also using diff(), see Figure 1.
The internal variables are well updated in each iteration since the paper delivers the dissipation curves of the Clasius-Duhem inequality of the example I am using to validate, see Figure 2.
Because of the above, I think there is something in my code that is not correct. It is worth mentioning that I also tried to do it with mixed u-p formulation, but it reached almost the same results. The following is part of my code:
## Function Space
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
un = Function(V, name = "Displacement")
v = TestFunction(V)
du = TrialFunction(V)
# Tensorial Space
ET2 = TensorFunctionSpace(mesh, 'DG', 0, shape=(3,3))
# Scalar Space
ET0 = FunctionSpace(mesh, "CG", 1)
# Internal Variable
Cmv = Function(ET2, name = "inverse matrix Cv")
#** The 2nd internal variable I treat it as a scalar
## Dirichlet Boundary Conditions
bc1 = DirichletBC(V.sub(0), Constant(0.), facets, 1) # x = 0 >> UX = 0
bc2 = DirichletBC(V.sub(2), Constant(0.), facets, 3) # z = 0 >> UZ = 0
bc3 = DirichletBC(V.sub(1), Constant(0.), facets, 5) # y = 0 >> UY = 0
tDirBC = Expression(('eps_d*time_'),eps_d = 0.0, time_=0.0 , degree=0)
bc4 = DirichletBC(V.sub(0), tDirBC, facets, 2) # x = 1 >> UX = disp
bc5 = DirichletBC(V.sub(2), tDirBC, facets, 4) # z = 1 >> UZ = disp
bcs = [bc1,bc2,bc3,bc4,bc5]
## Neumann Boundary Conditions
B = Constant((0.0, 0.0, 0.0)) # Body force per unit volume
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
## Fiber direction
theta = np.radians(0.0)
a0 = as_vector([cos(theta),sin(theta),0.0]) # Directional vector
A0 = outer(a0,a0) # Structural 2nd order tensor
## Initialize variables
I4_ = Constant(1.0) # [I₄ᶦˢᵒ]⁰ = 1 (Coefficient)
I4v = Constant(1.0) # [I₄ᵥ]⁰ = 1 (Coefficient)
I4v_ = Constant(1.0) # [I₄ᵥᶦˢᵒ]⁰ = 1 (Coefficient)
# Internal Variables
Cmv.assign(project(I, ET2)) # [(Cᵐᵥ)⁻¹]⁰ = I (Coefficient)
I4e_ = Constant(1) # [I₄ₑᶦˢᵒ]⁰ = 1 (Coefficient)
## Energy Functional
def FreeEnergy(u,Cmv,I4e_):
"""
[Input] u : [u]ⁿ (displacement field)
Cmv : [(Cᵐᵥ)⁻¹]ⁿ (matrix)
I4e_ : [I₄ₑᶦˢᵒ]ⁿ (fiber)
[Output] Ψₙ₊₁ : Free Energy function
"""
## Kinematics
(...)
# Invariants
I1 = tr(C)
I4 = inner(C,A0)
# Isochoric Invariants
I1_ = J**(-2/3)*I1
Ie_ = inner(C_,Cmv)
I4_ = inner(C_,A0)
## Free Energy
PSImeq = Constant(c)*(I1_ - Constant(3))
PSImvs = Constant(beta1)*Constant(c)*(Ie_ - Constant(3))
PSIfeq = Constant(k1/(2*k2))*(exp(Constant(k2)*(I4_ - 1)**2) - Constant(1))
PSIfvs = Constant(beta2*k1/(2*k2))*(exp(Constant(k2)*(I4e_ - 1)**2) - Constant(1))
PSIvol = Constant(0.5*k)*(J - 1)**2
return PSImeq + PSIfeq + PSIvol + PSImvs + PSIfvs
Pi = FreeEnergy(u,Cmv,I4e_)*dx - dot(B, u)*dx - dot(T, u)*ds(4)
DPi = derivative(Pi, u, v)
Jac = derivative(DPi, u, du)
problem = NonlinearVariationalProblem(DPi, u, bcs=bcs, J=Jac)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
for i, d in enumerate(Time):
if StretchLimit < 1.1:
tDirBC.time_ = d
else:
tDirBC.time_ = Time[i-j]
j += 1
## Solve Non-linear variational problem
# Input: [u]ⁿ, [(Cᵐᵥ)⁻¹]ⁿ, [I₄ₑᶦˢᵒ]ⁿ # Ouput: [u]ⁿ⁺¹
solver.solve()
un.assign(u)
## Update Kinematics
(...)
## Update Internal Variables
Cmv0 = Cmv # [(Cᵐᵥ)⁻¹]ⁿ (Coefficient)
I4v0 = I4v # [I₄ᵥ]ⁿ (Coefficient)
Cmv, I4v, I4e_, alpha = IntVar(Cn,Cn_,Cmv0,I4v0,pt,dt) # [(Cᵐᵥ)⁻¹]ⁿ⁺¹, [I₄ᵥ]ⁿ⁺¹, [I₄ₑᶦˢᵒ]ⁿ⁺¹, α
## Update Strain
Farray = np.reshape(evaluate_function(local_project(Fn, ET2),pt),[3,3]) # [F]ⁿ⁺¹ (Array)
R, U = polar(Farray) # Rotation, Stretch Tensor
H = np.log(U) # [H]ⁿ⁺¹ Hencky Strain
StretchLimit = U[0][0] # Relaxation stretch point [Liu et al. 2019]
## Update Stress
# 2nd Piola-Kirchhoff
Sniso_eq, Sniso_vis, Snani_eq, Snani_vis, hydro = Stress(un,Cmv,I4_,I4e_)
Sn = Sniso_eq + Sniso_vis + Snani_eq + Snani_vis + hydro
Seq = Sniso_eq + Snani_eq + hydro
# Cauchy
Cauchy = inv(Jn)*Fn*Sn*Fn.T
CauchyEq = inv(Jn)*Fn*Seq*Fn.T
I appreciate any help you can give me, I am really stuck on this.
