Non linear constitutive law

General
1

Hi all,

I am trying to model a non linear material for a simple axial force problem where the constitutive law is defined as,

\sigma = E_0*\epsilon + E_1*\epsilon^3

I have tried to work it out a Non linear poisson equation but, coefficient of \del u is longer A(u) but it is a function of \del u.

Can anyone please provide a starting point how I can model this problem?

Summary

This text will be hidden

2

Hi,
you should be more precise about the shape of \sigma, \epsilon ? are they scalar quantities ? are you solving in 1D ?

For general nonlinear elasticity, I recommend that you look at the hyperelasticity tutorial
https://jsdokken.com/dolfinx-tutorial/chapter2/hyperelasticity.html

It is not required to start with a potential \psi, you could also directly implement the UFL expression for P as a function of grad(u). The main point is that you need to define a non-linear 1-form, in your case something like:

eps = ufl.grad(u)
sig = E_0*eps + E_1*ufl.power(eps, 3)
F = ufl.inner(ufl.grad(v), sig)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)

but again you need to make things more explicit regarding the shapes of eps, sig, the meaning of eps**3 in this case, etc…

1 Like
3

@bleyerj Hi, thankyou for the reply.

Sorry for not being clear. Yes, I am solving a 1d axial force problem.
To get started with modelling non linear constitutive laws, I took cube of the strain

I have actually tried it solving this way but from energy formulation following Defining non-linear constitutive laws

from dolfin import *
from dolfin.cpp.mesh import *
from mshr import *
#import matplotlib.pyplot as plt
import numpy as np

L = 40.
H = 10.

res = 100

base = Rectangle(Point(0,0), Point(L,H))

mesh = generate_mesh(base, res)
x = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

class Left_face(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0)
     
class Right_face(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],L)

left_face = Left_face()
right_face = Right_face()

u = Function(V, name='Displacement')

eps = sym(grad(u))

E_0 = 2e5
E_1 = E_0/100

sigma = (E_0*eps + (E_1*(eps**3)))

bc1 = DirichletBC(V, Constant((0,0)), sub_domains, 1)
#bc = DirichletBC(V, Constant((0.,0.)), left_face)

f1 = Constant((100,0)) 
du = TrialFunction(V)
u_ = TestFunction(V)
l = dot(f1,u)*ds(2) 

# Energy:
energy = (0.5*inner(sigma,eps) - l)*dx



# Solve with homogeneous Dirichlet boundary conditions:
solve(derivative(energy,u)==0, u, bc1)

I was getting an error “Cannot take the power of a non-scalar expression .”

I guess I was not using ufl properly. I will do the appropriate changes and try the same again

4

The cube of a strain means something for a scalar strain only. For a tensorial quantity in 2D or 3D, you can define a cube such as:
dot(eps, dot(eps, eps))
but I don’t think this makes a realistic stress/strain relation