Hi,
I am solving a finite viscoelastic problem. The second PK stress tensor for the material is given as
where C1e and C1oe are material stiffness parameters, C is right Cauchy-Green deformation tensor. C can be decomposed into an elastic, Ce, and an inelastic part, Ci. So in the equation above Ci is the inelastic part of C. It follows a time-dependent equation as follows
where 
is the viscosity of the material and (.)^D is the deviatoric part of the tensor.
For small strain, the solution of this model should match with that of the standard linear solid (SLS) model. Here is my code to solve this problem. When the loading rate is very high during an uniaxial loading in the x direction, i.e., stretching is done in T=0.1 sec, Fenics solution matches with the analytical solution of SLS model. However, for a slower loading rate, say stretching in T =1 sec, Fenics solution does not match with analytical solution.
I believe, there is something wrong in the way I update my variables. Could anyone please help me to identify the error.
Here is my code
from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import math
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
def geometry_3d():
#mesh
mesh = UnitCubeMesh(6, 6, 6)
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
x0 = AutoSubDomain(lambda x: near(x[0], 0))
x1 = AutoSubDomain(lambda x: near(x[0], 1))
y0 = AutoSubDomain(lambda x: near(x[1], 0))
z0 = AutoSubDomain(lambda x: near(x[2], 0))
x0.mark(boundary_parts, 1)
y0.mark(boundary_parts, 2)
z0.mark(boundary_parts, 3)
x1.mark(boundary_parts, 4)
return boundary_parts
# Kinematics
def CalculatePKStress(u,pressure,Ci,ce,coe,eta):
I = Identity(V.mesh().geometry().dim()) ### Identity tensor
F = I + grad(u) ### Deformation gradient
J = det(F)
C = F.T*F ### Deformation tensor
T = 2*ce*I+2*coe * inv(Ci)-pressure * inv(C) ### second PK stress tensor
return T, (J-1)
# slope of Ci, inelastic strain tensor
def CalcGradCi(u_old,Ci_old,coe,eta):
I = Identity(V.mesh().geometry().dim())
F_old = I + grad(u_old)
C_old = F_old.T*F_old
Toe_old= 2*coe*inv(Ci_old)
slope =(2/eta)*C_old*dev(Toe_old)*Ci_old
return slope
# Create mesh and define function space ============================================
facet_function = geometry_3d()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())
# Limit quadrature degree
dx = dx(degree=4)
# Time stepping parameters
Tend = 1 # simulation time
Nsteps = 10 # number of steps
time = np.linspace(0, Tend, Nsteps+1)
dt = time[1]-time[0]
# material parameters
ce, coe, eta = 2, 2, 13.34 #3D model parameters
E1, E2, eta_sls = 12, 12, 20 # SLS model parameters
tau_sigma, tau_eps = (eta_sls/E1)*(1+E1/E2), eta_sls/E2 # SLS model time constant
eps_dot = 0.01/Tend # applied strain rate to stretch the cube in x direction
# Create function space
element_v = VectorElement("P", mesh.ufl_cell(), 1)
element_s = FiniteElement("P", mesh.ufl_cell(), 1)
mixed_element = MixedElement([element_v, element_s])
V = FunctionSpace(mesh, mixed_element)
# Define test and trial functions
dup = TrialFunction(V)
_u, _p = TestFunctions(V)
_u_p = Function(V)
u, p = split(_u_p)
_u_p_old = Function(V)
u_old, p_old = split(_u_p_old)
# define Ci_old
WF = TensorFunctionSpace(mesh, "DG", degree=0)
I = Identity(WF.mesh().geometry().dim())
initial_value = project(I, WF)
Ci_old = interpolate(initial_value, WF)
# Create tensor function spaces for stress and stretch output
defGrad = Function(WF, name='F')
PK_stress = Function(WF, name='PK1/2')
# Define Dirichlet boundary
bc0 = DirichletBC(V.sub(0).sub(0), Constant(0.), facet_function, 1)
bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facet_function, 2)
bc2 = DirichletBC(V.sub(0).sub(2), Constant(0.), facet_function, 3)
tDirBC = Expression(('eps_d*time_'),eps_d = eps_dot, time_=0.0 , degree=0)
bc3 = DirichletBC(V.sub(0).sub(0), tDirBC, facet_function, 4)
bcs = [bc0,bc1,bc2,bc3]
# # Total potential energy
I = Identity(V.mesh().geometry().dim())
pkstrs, hydpress = CalculatePKStress(u,p,Ci_old,ce,coe,eta)
F_cur = I + grad(u)
F1 = inner(dot(F_cur,pkstrs), grad(_u))*dx # weak form
F2 = hydpress*_p*dx # weak form for pressue term
F = F1+F2
J = derivative(F, _u_p,dup)
# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(F, _u_p, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
# Save solution in VTK format
file_results = XDMFFile("./Results/TestUniaxialLoading/Uniaxial.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
stretch_vec = np.zeros((Nsteps+1,1))
stress_vec = np.zeros((Nsteps+1,1))
stress_ana = np.zeros((Nsteps+1,1))
for i in range(len(time)):
tt = time[i]
print(i, 'time: ', tt)
# update time-dependent variable
tDirBC.time_ = tt
# Ci calculation
slope = CalcGradCi(u_old,Ci_old,coe,eta)
Ci = Ci_old+Constant(dt)*slope
Ci_proj = project(Ci, WF)
Ci_old.vector()[:] = Ci_proj.vector()
# solve and save disp
solver.solve()
# update old variables
_u_p_old.assign(_u_p)
# Extract solution components
u_print, p_print = _u_p.split()
# get stretch at a point for plotting
point = (0.5,0.5,0.5) # pick a point for plotting stretch and stress
DF = I + grad(u_print)
defGrad.assign(project(DF, WF))
stretch_vec[i] = defGrad(point)[0]
# get stress at a point for plotting
PK_s,tempPress = CalculatePKStress(u_print,p_print,Ci,ce,coe,eta)
PK_stress.assign(project(PK_s, WF))
stress_vec[i] = PK_stress(point)[0]
# get analytical solution for SLS model for uniaxial stretching
for i in range(len(time)):
stress_ana[i] = E1*eps_dot*(tau_sigma-tau_eps)*(1-np.exp(-time[i]/tau_eps))+E1*(stretch_vec[i]-1 )
# plot solutions to compare Fenics and analytical solution
f = plt.figure(figsize=(12,6))
plt.plot(stretch_vec, stress_vec,'r-')
plt.plot(stretch_vec, stress_ana,'k.')
This code produces results like this.
where the red line is the Fenics solution and the black dots are the analytical solution.
