ArityMismatch when implementing GreenLagrangeStrain / Hyperelasticity

1

Hello,

i am trying to implement the Green-Lagrange-Strain into my 3D-Problem.
As long as i use the InfinitesimalStrain and the PK2_lin stress tensor everything works out quite well, but when i try to use the GreenLagrangeStrain with PK2, i get the following error message:

ArityMismatch: Adding expressions with non-matching form arguments () vs (‘v_0’,).

# import modules
from fenics import *
from mshr import *
import numpy as np

# define geometry & parameters
thickness = 50e-6                               #thickness in m
diameter = 1.5e-3                               #diameter in m
radius = diameter/2
E = Constant(330e6)                             #Young's Modulus in Pa
nu = Constant(0.25)                             #Poisson's ratio
mu = Constant(E/(2.0*(1.0+nu)))
lmbda = Constant(E*nu/((1.0+nu)*(1.0-2.0*nu)))

# mesh geometry
No_Segments = int((diameter/thickness)+1)
membrane_domain = Cylinder(Point(0, 0, 0), Point(0, 0, thickness), radius, radius, 100) 
mesh = generate_mesh(membrane_domain, No_Segments)

#defining subdomains
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0.0) and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], thickness) and on_boundary

class Outer(SubDomain):
    def inside(self, x, on_boundary):
        return near(sqrt(x[0]*x[0]+x[1]*x[1]), radius, 4e-6) and on_boundary
    
#define function space for displacement, velocity and acceleration
V = VectorFunctionSpace(mesh, "CG", 1)

#test and trial functions
du = TrialFunction(V) 
u_ = TestFunction(V)
u = Function(V, name="Displacement")

# define subdomains
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_subdomains.set_all(0) #mark all with '0'
outer = Outer()
outer.mark(boundary_subdomains, 1)
top = Top()
top.mark(boundary_subdomains, 2)
bottom = Bottom()
bottom.mark(boundary_subdomains, 3)

#define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)

#set up boundary condition on outer region (dirichlet)
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, outer)
bcs = [bc]

file = XDMFFile("example/results.xdmf")
file.parameters["functions_share_mesh"] = True
file.parameters["flush_output"] = True
file.parameters["rewrite_function_mesh"] = False

def InfinitesimalStrain(u):
    return variable(0.5*(grad(u) + grad(u).T))

def PK2_lin(u):
    epsilon = InfinitesimalStrain(u)
    psi = lmbda/2*(tr(epsilon)**2) + mu*tr(epsilon*epsilon)
    S = variable(diff(psi, epsilon))
    return S

def GreenLagrangeStrain(u):
    d = len(u)
    I = Identity(d)
    F = I+grad(u)
    C = F.T * F
    E = 0.5*(C - I)
    return E
    
def PK2_GLS(u):
    E = GreenLagrangeStrain(u)
    psi = lmbda/2*(tr(E)**2) + mu*tr(E*E) #St. Venant Kirchhoff material
    #E = variable(E)
    S = variable(diff(psi, E))
    return S

#form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

p0 = 200000 #load amplitude in Pa
p = Constant((0, 0, p0))

#elastic stiffness form
def k(u, u_):
    return inner(PK2_lin(u), InfinitesimalStrain(u_))*dx 
    #return inner(PK2_GLS(u), GreenLagrangeStrain(u_))*dx 

#work of external forces
def Wext(u_):
    return dot(u_, p)*dss(3)

res = k(u, u_) - Wext(u_)

J = derivative(res, u)
problem = NonlinearVariationalProblem(res, u, bcs=bcs, J=J)
NLsolver = NonlinearVariationalSolver(problem)
NLsolver.parameters['newton_solver']['relative_tolerance'] = 1e-6
NLsolver.parameters['newton_solver']['linear_solver'] = 'mumps'

NLsolver.solve()
file.write(u, 0)

The used measures of stress and strain are defined in the elastic stiffness form at the end section of the code.

Do i have a general mistake in my formulation?
Any help to resolve this would be great :slight_smile:
Thanks in advance!

Joey

2

Please format your code using 3x` formatting, as the current formatting does not respect indentation. It would also be great if you could remove all code that is not related to the error message, and make this a minimal example.

3

ok, now it´s a minimum example :sweat_smile:

# import modules
from fenics import *
from mshr import *
import numpy as np

E = Constant(330e6)                             #Young's Modulus in Pa
nu = Constant(0.25)                             #Poisson's ratio
mu = Constant(E/(2.0*(1.0+nu)))
lmbda = Constant(E*nu/((1.0+nu)*(1.0-2.0*nu)))

mesh = UnitCubeMesh(20, 20, 20)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0.0) and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 1.0) and on_boundary
    
V = VectorFunctionSpace(mesh, "CG", 1)
du = TrialFunction(V) 
u_ = TestFunction(V)
u = Function(V, name="Displacement")

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_subdomains.set_all(0) #mark all with '0'
top = Top()
top.mark(boundary_subdomains, 2)
bottom = Bottom()
bottom.mark(boundary_subdomains, 3)
dss = ds(subdomain_data=boundary_subdomains)
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, top)
bcs = [bc]

file = XDMFFile("example/results.xdmf")
file.parameters["functions_share_mesh"] = True
file.parameters["flush_output"] = True
file.parameters["rewrite_function_mesh"] = False

def InfinitesimalStrain(u):
    return variable(0.5*(grad(u) + grad(u).T))

def PK2_lin(u):
    epsilon = InfinitesimalStrain(u)
    psi = lmbda/2*(tr(epsilon)**2) + mu*tr(epsilon*epsilon)
    S = variable(diff(psi, epsilon))
    return S

def GreenLagrangeStrain(u):
    d = len(u)
    I = Identity(d)
    F = I+grad(u)
    C = F.T * F
    E = 0.5*(C - I)
    return E
    
def PK2_GLS(u):
    E = GreenLagrangeStrain(u)
    psi = lmbda/2*(tr(E)**2) + mu*tr(E*E) #St. Venant Kirchhoff material
    #E = variable(E)
    S = variable(diff(psi, E))
    return S

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
p = Constant((0, 0, 200000))

#elastic stiffness form
def k(u, u_):
    return inner(PK2_lin(u), InfinitesimalStrain(u_))*dx 
    #return inner(PK2_GLS(u), GreenLagrangeStrain(u_))*dx 

#work of external forces
def Wext(u_):
    return dot(u_, p)*dss(3)

res = k(u, u_) - Wext(u_)

J = derivative(res, u)
problem = NonlinearVariationalProblem(res, u, bcs=bcs, J=J)
NLsolver = NonlinearVariationalSolver(problem)
NLsolver.parameters['newton_solver']['relative_tolerance'] = 1e-6
NLsolver.parameters['newton_solver']['linear_solver'] = 'mumps'

NLsolver.solve()
file.write(u, 0)

thanks for any help :slight_smile:

4

Hi, first psi must be defined in terms of the variable not in terms of E prior to declaring it as a variable:

def PK2_GLS(u):
    E = variable(GreenLagrangeStrain(u))
    psi = lmbda / 2 * (tr(E) ** 2) + mu * tr(E * E)  # St. Venant Kirchhoff material
    S = diff(psi, E)
    return S

second your residual definition is wrong, the form related to internal forces is:
\int_\Omega PK2(u) : \delta E(\widehat{u}) \text{dx}
where \delta E is the first-variation of the Green Lagrange strain in direction \widehat{u}. So something like:

# elastic stiffness form
def k(u, u_):
    dE = derivative(GreenLagrangeStrain(u), u, u_)
    return inner(PK2_GLS(u), dE) * dx
1 Like
5

Thank you a lot @bleyerj :grin: