Problems with the evolution equation for finite-strain Viscoelasticity

Linear Algebra
1

Hello Everyone,

I’m trying to implement a finite-viscoelastic model but I have a problem with the evolution equation. I implemented the internal variable as a DG-Tensorelement like it was proposed in the tutorial for linear viscoelasticity by Jeremy Bleyer linear_viscoelasticity or done here https://fenicsproject.discourse.group/t/issues-with-a-viscoelastic-problem/3434/10.
But it fails to converge on the 2nd iteration.

I am using an Ogden-Model for the equilibrium and non-equilibrium strain. The evolution equation has the following form:
\dot{\boldsymbol{b}}_e = \boldsymbol{l} \cdot \boldsymbol{b}_e + \boldsymbol{b}_e \cdot \boldsymbol{l}^T -(2\mathcal{V}^{-1}:\tau_\text{NEQ})\cdot \boldsymbol{b}_e
with \boldsymbol{b}_e being the internal variable, \boldsymbol{l} the velocity gradient, :\tau_\text{NEQ} the non-equilibrium stress and \mathcal{V}^{-1} a 4th order Tensor.
In the model the non-equilibrium stress is a function of the internal variables \boldsymbol{b}_e and the equilibrium strain of the left-cauchy-Green \boldsymbol{b}.

implementation of the Elements:

P_disp = VectorElement("CG", mesh.ufl_cell(), 2)  # displacement
P_visc = TensorElement("DG", mesh.ufl_cell(), 0 , symmetry=True)  # element for the viscous left cauchy green
mixed_element = MixedElement([P_disp, P_visc])   
V = FunctionSpace(mesh, mixed_element)


### Define test and trial functions
du, dbe = TestFunctions(V)
w_trial = TrialFunction(V)  
w = Function(V, name="Variables at current step")
(u, be) = split(w)
w_old = Function(V, name="Variables at previous step")
(u_old, be_old) = split(w_old)

implementation of the evolution equation:

F = I + grad(u)     # deformation gradient, t=n+1
F_old = I + grad(u_old)     # deformation gradient, t=n
L = ((F - F_old)/dt)*inv(F) # velocity gradient

pred_step = L*be + be*L.T
tmp =  1/eta_d * dev(T_neq) + 2/(9*eta_v) * inner(T_neq, I)*I # inner(4th order tensor, T_NEQ) 
Lie_be = -dot(tmp, be) # zum Zeitpunkt t = n+1
evol = 1/dt*inner((be - be_old), dbe)*dx - inner(pred_step, dbe)*dx - inner( Lie_be, dbe)*dx

The total code is a bit lengthy, but if it helps I can also post it.

Thanks for your help in advance!

2

Hi Martin,

Finite viscoelasticity is a highly nontrivial problem. The simple solve syntactic sugar won’t work (I assume you would have figured that out already). It does, however, work for some simple meshes. A POC is here (FEniCS_Shrimali_Ghosh_Kumar_Lopez-Pamies/viscoFEM.py at main · bhaveshshrimali/FEniCS_Shrimali_Ghosh_Kumar_Lopez-Pamies · GitHub) in case if you are interested. The code isn’t super well documented, but I would be happy to help with any questions if you have.

For real problems, with more intricate meshes I would recommend designing your own solver using either of NewtonSolver or PETScSNESSolver classes or using petsc4py directly.

3 Likes
3

Hello Bhaveshshrimali,

thanks for your code and your kind offer!
the paper you referenced in your code is also very informative, I will look into it!

4

Hi,
A couple of more recent papers (see here and here) may also help you if want to look into the formulation, discretization and solvers in more detail than the paper referenced in the code

5

In your code you minimize the free energy and apply a jump condition

a_uv = derivative(freeEnergy(u, Cv), u, delu) * dx + qvals / \
       h_avg * inner(jump(u), jump(delu)) * dS

instead of using the classical momentum balance

a_uv = inner(stressPiola(u, Cv), grad(delu))*dx

is the first one more efficent/error prone?

6

Hi Martin,
The two are equivalent for a conforming discretization (such as one using Lagrange elements), however in the code we are using a Crouzeix-Raviart discretization for the displacement (which is continuous only at the facet midpoint), hence the stabilization term penalizing the jump.

I can elaborate a bit more on this when I’m back on my computer…