Ogden model for finite elasticity 2D in variational formulation

variational formulation
1

Hi
I’m trying to build a 2D Ogden problem in fenics, but I obtain a constant zero displacement field as solution. Thanks for the attention and for the answers.

from dolfin import *
from mshr import*
import matplotlib.pyplot as plt
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

#Geometry
lpx=50
lpy=10
n_mesh=30
domain=Rectangle(Point(0,0),Point(lpx,lpy))
mesh=generate_mesh(domain, n_mesh)

U=VectorFunctionSpace(mesh,"P",2,2)
tollb=1E-13
def boundary_left(x): return abs(x[0]) < tollb
def boundary_right(x):return abs(x[0]-lpx) < tollb

# Material's parameter
alfa_1=1
mu_1=10
alfa_2=2
mu_2=10
D_1=10
u=Function(U)
Trial_u=TrialFunction(U)
Test_u=TestFunction(U)

F=grad(u)+Identity(2)
C=F.T*F
lambda_1=(tr(C)+ sqrt(tr(C)**2-4*det(C)))/2
lambda_2=(tr(C)- sqrt(tr(C)**2-4*det(C)))/2
t1=(mu_1/alfa_1)*((lambda_1**alfa_1)+(lambda_2**alfa_1)-2)
t2=(mu_2/alfa_2)*((lambda_1**alfa_2)+(lambda_2**alfa_2)-2)
t3=(1/D_1)*((det(C)-1)**2)
Energy=(t1+t2+t3)*dx

#Energy=Ogd_Energy(u)*dx
var1=derivative(Energy,u,Test_u)
var2=derivative(var1,u,Trial_u)
bc_left=DirichletBC(U, Constant((0.0,0.0)), boundary_left)
bc_right=DirichletBC(U,Constant((50.0,0.0)),boundary_right)
bc=[bc_left,bc_right]
problem = NonlinearVariationalProblem(var1,u,bc,J=var2)
solver = NonlinearVariationalSolver(problem)

plt.figure(0)
c = plot(u.sub(0), mode='color')
plt.colorbar(c)
2 Likes
2

The immediate problem is that you need to call the solve method of solver after creating it:

solver.solve()

This results in a nan residual on the first iteration, but that is a common issue when using derivative on functionals involving sqrt, where the argument to sqrt is zero (which is the case here when \mathbf{C} is identity on most of the domain in the initial guess, since, in 2D, \operatorname{tr}\mathbf{I} = 2 and \operatorname{det}\mathbf{I} = 1). (The derivative function differentiates symbolically and \frac{d}{dx}\sqrt{x} = \frac{1}{2\sqrt{x}}, so you divide by zero if evaluating the symbolic derivative at x = 0.) One way to get around this is to regularize the square root slightly, e.g.,

reg_sqrt = lambda x : sqrt(x + DOLFIN_EPS)
lambda_1=(tr(C)+ reg_sqrt(tr(C)**2-4*det(C)))/2
lambda_2=(tr(C)- reg_sqrt(tr(C)**2-4*det(C)))/2
4 Likes
3

Just to add to @kamensky 's response above, it may also help to apply the displacement gradually. Though it shouldn’t be a problem in this specific case, but in general it is a good idea to apply it gradually to allow Newton’s method to converge to the proper solution (also helpful for the case when you have a non-convex strain energy density).

As an aside you can play with how loose a regularization works for your problem. For instance, for this problem sqrt(DOLFIN_EPS + ..) doesn’t work for me but sqrt(2. * DOLFIN_EPS + ...) works.

2 Likes
4

thank you very much for the answers!!