Hello, I am working with a 3d matrix-fibre viscoelasticity model with 2 internal variables. The following equation is the free energy function:
\Psi\left(\overline{I}_1,\overline{I}_e,\overline{I}_4,\overline{I}_{4e},J\right) = c\left(\overline{I}_1-3\right) + \frac{1}{2}\kappa \left(J-1\right)^2 + \frac{k_1}{2k_2}\left(\exp\left[k_2\left(\overline{I}_4-1\right)^2\right]-1\right) \\ \hspace{100pt} + \beta_1 c (\overline{I}_e-3) + \beta_2 \frac{k_1}{2k_2}\left(\exp\left[k_2\left(\overline{I}_{4e}-1\right)^2\right]-1\right)
, where \beta_1, \beta_2 are the viscous parameters, and the internal variables are \overline{I}_{4e} and C^m_v (viscous-matrix part of right cauchy-green), from \overline{I}_{4}=\overline{I}_{4v}\overline{I}_{4e} and \overline{I}_e=\overline{C}:{C^m_v}^{-1}, respectively.
The evolution equations for these are:
\dot{\overline{{C^m_v}^{-1}}}=-\frac{1}{\tau_1}{C^m_v}^{-1}+\frac{3C^{-1}}{\tau_1\mathbf{I}:(C^m_vC^{-1})} \hspace{20pt} \text{and} \hspace{20pt} \dot{\overline{{I_{4e}}}}=-\frac{1}{\tau_2}\frac{\overline{I}_{4e}\overline{\Psi}_{4e}}{\overline{\Psi}_{4e4e}\overline{I}_{4e}+\overline{I}_{4e}}
Backward Euler was used to update these variables by means of an algorithm at each iteration.
My question is, what is the best way to define these internal variables? So far I did it like this:
# Tensorial Space
ET2 = TensorFunctionSpace(mesh, 'DG', 0, shape=(3,3))
# Internal Variables
Cmv = Function(ET2, name = "inverse matrix Cv") # (Cᵐᵥ)⁻¹
I4e_ = inner(C_,A0) # Ī₄ₑ
# Variational problem
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)
Then to update them I used project and reshape with numpy to implement the algorithm, I used them as array because in both cases I had to use newton rhapson at a certain point using scipy, and because I didn’t know how to do the same using fenics. For example, I must find [\overline{I}_{4e}]^{n+1} by solving this equation by using NR:
[\overline{I}_{4e}]^{n+1}+\frac{\frac{\Delta t}{\tau_2}[\overline{I}_{4e}]^{n+1}([\overline{I}_{4e}]^{n+1}-1)}{[\overline{I}_{4e}]^{n+1}[1+2k_2([\overline{I}_{4e}]^{n+1}-1)^2]+([\overline{I}_{4e}]^{n+1}-1)}=[\overline{I}_{4e}]^n
Comparing with a numerical example where they give the dissipation curves I could check that the values of my internal variables are well calculated using numpy and scipy (see this figure). However, my final results are not the best (Relaxation graph). There are two problems, the elastic part through the volumetric response is very high, and the rate of decay due to viscosity is very small.
Is there a way to define this as in the linear viscoelasticity fenics tutorial?, like this:
Ve = VectorElement("CG", mesh.ufl_cell(), 1)
Qe = TensorElement("DG", mesh.ufl_cell(), 0) # Internal variable viscous strain
W = FunctionSpace(mesh, MixedElement([Ve, Qe]))
But in my case if I defined both in that way I would not know how to deal with them within the algorithm I have to apply and the Newton-Rhapson procedure that goes with it.
Thanks for any help!
