Hi,
I’m not sure if this is still a problem but I think that - using Cardanos formula - you can make an Ogden material model working in Fenics. Below, I just copied some code that I implemented for exactly this purpose: it calculates eigenvalues and projection tensors for 2D- and 3D tensors and allows you to get the Ogden material model working.
###########################
# load required libraries #
###########################
import dolfin as df # dolfin module
import ufl # ufl module
import numpy as np # numpy
import copy as cp # copy for deepcopies
####################################################
# Method for analytical calculation of eigenvalues #
####################################################
def eigenValues(T, tol=df.DOLFIN_EPS):
'''
Method to calculate the (real) eigenvalues of a second order tensor T using
its principal invariants.
'''
dim=T.geometric_dimension() # get geometric dimension of tensor T to decide which subroutine to use for eigenvalue computations
# choose projection tensor routine based on dimensionality of problem
if dim == 2: # 2D-problem
return _eigVal2D(T,tol)
else: # 3D-problem
return _eigVal3D(T,tol)
########################################################################
# Method for the analytical calculation of eigenvalues for 2D-Problems #
########################################################################
def _eigVal2D(T, tol):
'''
Analytically calculate eigenvalues for a two-dimensional tensor with a
characteristic polynomial equation of the form
lambda^2 - I1*lambda + I3 = 0 .
Since the characteristic polynomial is a reduced quadratic equation, the
eigenvalues can be determined using the p-q-formula.
NOTE:
The method implemented here, implicitly assumes that the polynomial has
only real roots, since imaginary ones should not occur in this use case.
In order to ensure eigenvalues with algebraic multiplicity of 1, the idea
of numerical perturbations is adopted from "Computation of isotropic tensor
functions" by C. Miehe (1993). Since direct comparisons with conditionals
have proven to be very slow, not the eigenvalues but the coefficients
occuring during the calculation of them are perturbated to get distinct
values.
'''
# determine perturbation from tolerance
pert = 2*tol
# get required invariants
I1 = df.tr(T) # trace of tensor
I3 = df.det(T) # determinant of tensor
# determine discriminant
# -> If the discriminant is 0, eigenvalues with a algebraic multiplicity
# -> of 2 could occur. Here, this is prevented by "artificially creating"
# -> eigenvalues with an algebraic multiplicity of 1 by adding a numerical
# -> perturbation to discriminants that are close to 0. This is especially
# -> necessary since a zero discriminant will cause problems during the
# -> automatic differentiation procedure
D = I1**2-4*I3 # preliminary value for discriminant
D=ufl.conditional(ufl.lt(D,tol),abs(D)+pert,D) # add numerical perturbation to discriminant, if close to zero; ensure positiveness of D
# calculate polynomial roots
# -> Here, the polynomial solution formula is only used to calculate the
# -> "bigger" root, since - if I1 and sqrt(D) are of similar magnitude - a
# -> loss of significance could lead to erroneous results. For the second
# -> root, Vietas formula is used.
lambda1 = 0.5*(I1+ufl.sqrt(D))
lambda2 = I3/lambda1
# return polynomial roots (eigenvalues)
return lambda1, lambda2
########################################################################
# Method for the analytical calculation of eigenvalues for 3D-Problems #
########################################################################
def _eigVal3D(T, tol):
'''
Analytically calculate eigenvalues for a three-dimensional tensor T with a
characteristic polynomial equation of the form
lambda^3 - I1*lambda^2 + I2*lambda - I3 = 0 .
Since the characteristic polynomial is in its normal form , the eigenvalues
can be determined using Cardanos formula. This algorithm is based on:
"Efficient numerical diagonalization of hermitian 3 × 3 matrices" by
J. Kopp (equations: 21-34, with coefficients: c2=-I1, c1=I2, c0=-I3).
NOTE:
The method implemented here, implicitly assumes that the polynomial has
only real roots, since imaginary ones should not occur in this use case.
In order to ensure eigenvalues with algebraic multiplicity of 1, the idea
of numerical perturbations is adopted from "Computation of isotropic tensor
functions" by C. Miehe (1993). Since direct comparisons with conditionals
have proven to be very slow, not the eigenvalues but the coefficients
occuring during the calculation of them are perturbated to get distinct
values.
'''
# determine perturbation from tolerance
pert = 2*tol
# get required invariants
I1 = df.tr(T) # trace of tensor
I2 = 0.5*(df.tr(T)**2-df.inner(T,T)) # 2nd invariant of tensor
I3 = df.det(T) # determinant of tensor
# determine terms p and q according to the paper
# -> Follow the argumentation within the paper, to see why p must be
# -> positive. Additionally ensure non-zero denominators to avoid problems
# -> during the automatic differentiation
p = I1**2 - 3*I2 # preliminary value for p
p = ufl.conditional(ufl.lt(p,tol),abs(p)+pert,p) # add numerical perturbation to p, if close to zero; ensure positiveness of p
q = 27/2*I3 + I1**3 - 9/2*I1*I2 # preliminary value for q
q = ufl.conditional(ufl.lt(abs(q),tol),q+ufl.sign(q)*pert,q) # add numerical perturbation (with sign) to value of q, if close to zero
# determine angle phi for calculation of roots
phiNom2 = 27*( 1/4*I2**2*(p-I2) + I3*(27/4*I3-q) ) # preliminary value for squared nominator of expression for angle phi
phiNom2 = ufl.conditional(ufl.lt(phiNom2,tol),abs(phiNom2)+pert,phiNom2) # add numerical perturbation to ensure non-zero nominator expression for angle phi
phi = 1/3*ufl.atan_2(ufl.sqrt(phiNom2),q) # calculate angle phi
# calculate polynomial roots
lambda1 = 1/3*(ufl.sqrt(p)*2*ufl.cos(phi)+I1)
lambda2 = 1/3*(-ufl.sqrt(p)*(ufl.cos(phi)+ufl.sqrt(3)*ufl.sin(phi))+I1)
lambda3 = 1/3*(-ufl.sqrt(p)*(ufl.cos(phi)-ufl.sqrt(3)*ufl.sin(phi))+I1)
# return polynomial roots (eigenvalues)
return lambda1, lambda2, lambda3
##############################################################
# Method for an analytical calculation of projection tensors #
##############################################################
def projectionTensors(T, *eigenValues):
'''
Method to calculate the projection tensors of a second order tensor T using
its already known eigenvalues.
'''
dim=T.geometric_dimension() # get geometric dimension of tensor T to decide which subroutine to use for projection tensor computations
# choose projection tensor routine based on dimensionality of problem
if dim == 2: # 2D-problem
return _projTen2D(T,*eigenValues)
else: # 3D-problem
return _projTen3D(T,*eigenValues)
##############################################################################
# Method for an analytical calculation of projection tensors for 2D-problems #
##############################################################################
def _projTen2D(T, lambda1, lambda2):
'''
Analytically calculate projection tensors for a two-dimensional tensor T
with two eigenvalues of algebraic multiplicity 1 (distinct eigenvalues).
Mathematical manipulations show that the projection tensors of a simplified
two-dimensional problem can be computed independent of the third eigenvalue
by just using the simplified version of the three-dimensional problem:
M_a = T-lambda_b*I/(lambda_a-lambda_b)
with a,b={1,2} and b!=a
'''
# get required quantities:
I = ufl.Identity(2)
# calculate individual projection tensors
M1 = (T-lambda2*I)/(lambda1-lambda2)
M2 = (T-lambda1*I)/(lambda2-lambda1)
# return projection tensors
return M1, M2
##############################################################################
# Method for an analytical calculation of projection tensors for 3D-problems #
##############################################################################
def _projTen3D(T, lambda1, lambda2, lambda3):
'''
Analytically calculate projection tensors for a three-dimensional tensor T
with 3 eigenvalues of algebraic multiplicity 1 (distinct eigenvalues).
Using some mathematical manipulations on the initial eigenvalue problem and
the properties of projection tensors, they can be expressed in terms of the
tensor T and its eigenvalues as "Frobenius-covariants"
M_a = prod_b (T-lambda_b*I)/(lambda_a-lambda_b)
for a,b={1,2,3} and b!=a.
'''
# get required quantities
I = ufl.Identity(3)
# calculate the individual projection tensors
M1 = (T-lambda2*I)*(T-lambda3*I)/((lambda1-lambda2)*(lambda1-lambda3))
M2 = (T-lambda3*I)*(T-lambda1*I)/((lambda2-lambda3)*(lambda2-lambda1))
M3 = (T-lambda1*I)*(T-lambda2*I)/((lambda3-lambda1)*(lambda3-lambda2))
# return projection tensors
return M1, M2, M3
Using the projection tensors for the evaluation of stresses does not converge as good as using the derivative of the energy - I think this is normal since there might be problems, when you have eigenvalues with an algebraic muliplicity of more than 1. In my case, a simple version (2D with projection tensors) of using the model would look like this:
# material parameters
nOgden = 2
nu = 0.49
mu = [-6.77e3, 1.88e3] # mu parameters for matrix
alpha = [-6.88, 2.0] # alpha parameters (exponents) for matrix
kappa = sum([mu[i]*alpha[i] for i in range(nOgden)])*(1+nu)/(3*(1-2*nu)) # bulk modulus for matrix
# Kinematics
dim = u.geometric_dimension()
I = Identity(dim) # identity tensor
F = I+grad(u) # deformation gradient
detF = det(F) # Jacobi determinant
C = F.T*F # right Cauchy-Green tensor
# spectral decomposition for Ogden material
lam21, lam22 = eigenValues(C,tol=1e-10) # eigenvalues of C
M1, M2 = projectionTensors(C, lam21, lam22) # projection tensors for C
lam1_iso = detF**(-1./3.)*sqrt(lam21) # first eigenvalue of F_iso
lam2_iso = detF**(-1./3.)*sqrt(lam22) # second eigenvalue of F_iso
lam3_iso = detF**(-1./3.) # third eigenvalue of F_iso (2D case)
# stress
g = kappa/2*(detF*detF - 1) # reoccuring function for volumetric part of stress
T1_dev = 1/lam21*sum([mu[i]*(2/3*lam1_iso**alpha[i]-1/3*(lam2_iso**alpha[i]+lam3_iso**alpha[i])) for i in range(nOgden)]) # contribution of first eigenvalue to deviatoric part of 2nd Piola-Kirchhoff stress
T2_dev = 1/lam22*sum([mu[i]*(2/3*lam2_iso**alpha[i]-1/3*(lam1_iso**alpha[i]+lam3_iso**alpha[i])) for i in range(nOgden)]) # contribution of second eigenvalue to deviatoric part of 2nd Piola-Kirchhoff stress
T_vol = 1/lam21**2*g*M1 + 1/lam22*g*M2 # volumetric part of 2nd Piola-Kirchhoff stress
T_dev = T1_dev*M1 + T2_dev*M2 # deviatoric part of 2nd-Piola-Kirchhoff stress
T = T_dev + T_vol # 2nd Piola-Kirchhoff stress in matrix
P = T*F.T # 1st Piola-Kirchhoff stress in matrix
I also tried it in 3D and at least for simple cases, it worked quite good.
Hope this helps…