Advice on solving a more complex 2D linear elasticity problem

1

I’m very new to Python and the FEniCS project and I have managed to follow some of the tutorials to solve a 2D cantilever beam with a constant UDL. I then represented the mesh and displacement data through Paraview.

I wish to develop this further by adding multiple materials along the cantilever (eg. E = 150GPa from 0 to 2m, E = 200GPa from 2 to 3m) and adding either multiple point loads along the beam or representing these as a force distribution which will be in the form of a polynomial. The data I wish to retrieve is still just the displacement along the beam.

I would really appreciate any advice on how I should start to develop this.

Best regards.

2

Welcome,

this is a very good starting point for what you want to accomplish (constant parameter values over different subdomains):
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/materials.html

The force distribution is most easily represented as an Expression for the traction vector and incorporated directly into the boundary term in your weak form. That is, for the weak form of balance of linear momentum

\int_\Omega \sigma:\operatorname{grad} v\,\text d x - \int_{\Gamma_\mathrm{N}} t\cdot v\,\text ds=0

you would prescribe the traction t on a part of your boundary.

This tutorial
https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html
might get you started in the right direction (there’s also a part on different material parameters on subdomains)

3

Hi klunkean,

Thanks for your help, I really appreciate it.

I’ve managed to set different material properties over various subdomains and validated the results for the simple uniformly distributed load. So thanks for that!

I’ve tried a few different approaches to adding a force distribution by prescribing the traction but I have been unable to get it to solve and have hit a bit of a wall. I’ll try to add the part of my code that details this.

# Create finite element function space, V
V = VectorFunctionSpace(mesh, 'CG', 1)

# Define force distribution as traction expression
el = V.ufl_element()
t = UserExpression(('2e-11*x[1]**3-5e-7*x[1]**2+0.0009*x[1]+35.996'), element = el, degree = 2)

# Define boundary conditions

def root_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] < tol

bc = DirichletBC(V, (0, 0), root_boundary)

def t_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] > tol

bc1 = DirichletBC(V, t, t_boundary)

bcs = [bc, bc1]

# Mark boundaries
boundaries = MeshFunction('size_t', mesh, 2)
boundaries.set_all(0)
AutoSubDomain(root_boundary).mark(boundaries, 1)
AutoSubDomain(t_boundary).mark(boundaries, 2)
ds = ds[boundaries]

# Define stress and strain
## Strain tensor
def eps(v):
    return sym(grad(v))

## Stress tensor 
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) 

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() 
v = TestFunction(V)
f = Constant((0, 0))
F = inner(sigma(u), eps(v))*dx - dot(t,v)*ds(2)
a, L = lhs(F), rhs(F)

# Compute solution 
uh = Function(V)
solve(a == L, uh, bcs)

Apologies if this isn’t in the best format to read.

I get the error:

*** Error: Unable to evaluate expression.
*** Reason: Missing eval() function (must be overloaded).
*** Where: This error was encountered inside Expression.cpp.
*** Process: 0

Which I have been unable to find help for when looking it up prior to posting this.

Again, thanks for your help and I hope this question isn’t too tedious to answer!

4

It should work if you just write

t = Expression('2e-11*x[1]**3-5e-7*x[1]**2+0.0009*x[1]+35.996', element = el)

Can you also please format your code using

```python    
code goes here
```

Then we can copy it and run it directly.

5

Great, that definitely moved it along, thank you.

And I see, I was trying to figure that one out.

Sorry to keep pestering but I am now getting the error:

error: invalid type argument of unary '’ (have ‘int’)
values[0] = 2e-11x[0]**3-5e-7x[0]**2+0.0009x[0]+35.996;
^
error: invalid type argument of unary '’ (have ‘int’)
values[0] = 2e-11x[0]**3-5e-7x[0]**2+0.0009x[0]+35.996;
^
Do you know what could be the issue?

Thanks again.

6

Oh of course. I didn’t think too much about the string in your Expression.

That actually is C++ code so you need to use functions defined in the C++ standard libraries (like the cmath library). You need to use pow(x[0],2) instead of the python x[0]**2. ** has a different meaning in C++.

7

Ah, that sorted that issue. Again, great to know.

Sorry if it feels like your basically solving my whole code. There seems to be multiple issues.

This the error I’m now getting:

—> bc1 = DirichletBC(V, t, t_boundary)

*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Expecting a vector-valued boundary value but given function is scalar.
*** Where: This error was encountered inside DirichletBC.cpp.
*** Process: 0

I’ve tried replacing the line causing the issue with bc1 = DirichletBC(V, (0, t), t_boundary) but then receive the error:

TypeError: Expression.float returned non-float (type NotImplementedType)

Apologies for so many questions, you’ve helped a lot already so no worries if this is becoming too tedious for you to answer.

Cheers, again.

8

Oh okay, so there we found another problem in your code that I missed when first skimming through. Mainly due to the wrong formatting. Please consider using the edit function to properly format your code in the first post.

What you are trying to do is prescribing Dirichlet boundary conditions and Neumann boundary conditions on the same surface (t_boundary with ID 2). That’s not even physically sound. I’d suggest reading up a bit more on the boundary conditions in the second tutorial link I posted above.

If you want to prescribe forces, or to be more precise, tractions, you just leave the Neumann term dot(t,v)*ds(2) in your form and define a proper, that is vector valued, expression for t. For traction in y-Direction you could for example do
t = Expression(('0.0', '5.0'), element = el)

The specific error you get arises from trying to prescribe a scalar valued expression as a Dirichlet value for a vector valued function space. But since you don’t even want to prescribe Dirichlet BC on your traction boundary we don’t need to care about that error right now. Just remove your bc1 from your code.

1 Like
9

Okay, so I’ve started from scratch to produce a simpler model with a triangular distributed load from 50N/mm to 0 acting in the positive y-direction on a 2d cantilever beam (1000*20mm). This should be easier to understand how to create a sound model.

from dolfin import *
from mshr import *
from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt

# Material properties
L = 1000                        # Length mm
H = 20                          # Length mm 
nu = 0.3                        # Poisson's ratio
rho = 1.5e-9                    # Density kg/mm3   
E = 1e9                         # Young's modulus Pa

# Lame's parameters
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Generate simple beam (1000*20)
Beam = Rectangle(Point(0, 0), Point(1000, 20))
mesh = generate_mesh(Beam, 10)

# Create finite element function space, V
V = VectorFunctionSpace(mesh, 'CG', 1)

# Define root boundary condition
def root_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), root_boundary)

# Strain tensor
def eps(v):
    return sym(grad(v))

# Stress tensor 
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) 

# Define force distribution
el = V.ufl_element()
f = Expression(('0.0', '(-0.05*x[0]+50)'), element = el) # Equation for triangular distributed load with root N/mm = 50

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() 
v = TestFunction(V)
t = Constant((0, 0))
a = inner(sigma(u), eps(v))*dx
L = dot(f, v)*dx + dot(t, v)*ds

# Compute solution 
u = Function(V)
solve(a == L, u, bc)

plot(u, title='Displacement', mode='displacement')

# Compute magnitude of displacement
V = FunctionSpace(mesh, 'CG', 1)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
     u_magnitude.vector().get_local().min(),
     u_magnitude.vector().get_local().max())

Euler = 50*1000**4/(30*E*((20**4)/12))
print(Euler)

plot(u, title='Displacement', mode='displacement')

# Test expression is providing correct distribution
Vf = VectorFunctionSpace(mesh, 'CG', 1)
f = interpolate(f, Vf)
values = f.vector()
plt.figure()
plt.plot(values)

Upon comparing it to beam theory calculations, the FEniCS model is way off. I’ve also noticed that changing the mesh resolution impacts the deflection also which must be connected to the expression. I’m guessing that the more elements I add, the more force it applies. I added a little test at the end of the script to display the distribution and also noticed that at a certain resolution it flips.

Could you give me some pointers in what’s incorrect and hopefully I’ll be able to gain some understanding and get a working model.

Thanks again.

Edit:
I just tested a theory using mesh.num_cells() and divided the expression by this number (Basically getting force per cell) which resulted in a much more reasonable displacement. What’s your take on this?

10

Hi,
three problems:

  • your are using P1 elements with a mesh consisting of only one element along the beam height, this is very bad for bending. Use either a much refined mesh or P2 elements, or even both
  • your comparison with beam theory is wrong. The beam theory deflection is w=qL^4/(30EI) with I=bh^3/12 with b the beam width in the out-of-plane direction. In the above formula, q is a load per unit length (in N/m) where in your calculation you use a body force f=50 N/mm^3. In the comparison you should therefore put for the load per unit length q=f\cdot bh so that the deflection now becomes w=fL^4/(30h^2/12).
  • finally, here you are doing a 2D calculation in plane strain condition whereas beam theory assumes plane-stress conditions. You must therefore change the effective lmbda Lamé coefficient for plane stress conditions (see https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html) or use a zero Poisson ratio.
2 Likes
11

Perfect, this resulted in a difference of 0.0011mm in deflection between beam theory and FEniCS. Thank you.

Just out of curiosity, do you know why adding more elements appears to flip the distribution? I’ll add the final code and you can see that either changing the resolution of the mesh form 10 to 50 flips the distribution plot as well as changing the ‘Vf’ VectorFunctionSpace (in the final part of the script) to P2 elements. This doesn’t seem to change the deflection so doesn’t appear to be an issue so far.

from dolfin import *
from mshr import *
from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt

# Material properties
L = 1000                        # Length mm
H = 20                          # Length mm 
nu = 0                          # Poisson's ratio = 0 for plane-stress conditions
rho = 1.5e-9                    # Density kg/mm3   
E = 1e9                         # Young's modulus Pa

# Lame's parameters
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Generate simple beam (1000*20)
Beam = Rectangle(Point(0, 0), Point(1000, 20))
mesh = generate_mesh(Beam, 10)

# Create finite element function space, V
V = VectorFunctionSpace(mesh, 'CG', 2) # P2 Elements?

# Define root boundary condition
def root_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), root_boundary)

# Strain tensor
def eps(v):
    return sym(grad(v))

# Stress tensor 
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) 

# Number of cells in mesh
nc = mesh.num_cells()

# Define force distribution
el = V.ufl_element()
f = Expression(('0.0', '(-0.05*x[0]+50)'), element = el) # Equation for triangular distributed load with root N/mm = 50

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() 
v = TestFunction(V)
t = Constant((0, 0))
a = inner(sigma(u), eps(v))*dx
L = dot(f, v)*dx + dot(t, v)*ds

# Compute solution 
u = Function(V)
solve(a == L, u, bc)

plot(u, title='Displacement', mode='displacement')

# Compute magnitude of displacement
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
     u_magnitude.vector().get_local().min(),
     u_magnitude.vector().get_local().max())

#Euler = (50*20*20)*1000**4/(30*E*((20**4)/12))
Euler = (50*1000)/(30*(20**2)/12)
print(Euler)

plot(u, title='Displacement', mode='displacement')

# Test expression is providing correct distribution
Vf = VectorFunctionSpace(mesh, 'CG', 1)
f = interpolate(f, Vf)
values = f.vector()
plt.figure()
plt.plot(values)

Thanks again.

13

The vector of degrees of freedom is absolutely not ordered in any meaningful way so that the plot of values makes no sense and changes when refining the mesh or changing the interpolation degree since the dof ordering changes each time.

1 Like
16

Hi again,

I’ve still been unable to completely solve this problem. This is my current script.

from dolfin import *
from mshr import *
from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt

# Material properties
Le = 10000                      # Length mm
H = 50                          # Length mm 
nu = 0.3                        # Poisson's ratio = 0 for plane-stress conditions
model = 'plane_stress'

# Generate spar
sparinner = Rectangle(Point(0, 0), Point(5000, 50))
sparouter = Rectangle(Point(5000, 0), Point(10000, 50))

spar = sparinner + sparouter
mesh = generate_mesh(spar, 10)

# Create finite element function space, V (P2 Elements)
V = VectorFunctionSpace(mesh, 'CG', 2) 

# Create two materials of different Young's modulus throughout the beam using expressions
E = [100e9, 100e9]
mu_list = []
for x in E:
    mu_list.append(x/2/(1+nu))

lmbda_list = []
for x in E:
    lmbda_list.append(x*nu/(1+nu)/(1-2*nu))

lmbdaplane_list = []
for lmbda_value, mu_value in zip(lmbda_list, mu_list):
    lmbdaplane_list.append(2*mu_value*lmbda_value/(lmbda_value+2*mu_value))

if model == 'plane_stress':
    lmbda_list = lmbdaplane_list

tol = 1E-14
mu0 = mu_list[0]
mu1 = mu_list[1]
mu = Expression('x[0] <= 5000 + tol ? mu0 : mu1' , degree = 0, tol=tol, mu0=mu0, mu1=mu1)

tol = 1E-14
lmbda0 = lmbda_list[0]
lmbda1 = lmbda_list[1]
lmbda = Expression('x[0] <= 5000 + tol ? lmbda0 : lmbda1' , degree = 0, tol=tol, lmbda0=lmbda0, lmbda1=lmbda1)

# Define root boundary condition
def root_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), root_boundary)

# Strain tensor
def eps(v):
    return sym(grad(v))

# Stress tensor 
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) 

# Define force distribution
H2 = H*2
el = V.ufl_element()
t = Expression(('0.0', '((-5*x[0]+50e3))/H2'), H2=H2, element = el)   # Equation for triangular distributed load with root kN/mm = 50
                                                                      # Body force = 0.02kN/mm^3
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() 
v = TestFunction(V)
f = Constant((0, 0))
a = inner(sigma(u), eps(v))*dx
L = dot(f, v)*dx + dot(t, v)*ds

# Compute solution 
u = Function(V)
solve(a == L, u, bc)

plot(u, title='Displacement', mode='displacement')

# Compute magnitude of displacement
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
     u_magnitude.vector().get_local().min(),
     u_magnitude.vector().get_local().max())

# Estimate max displacement for comparison
E = 100e9
Euler = (50e3)*Le**4/(30*E*((H**4)/12))
print(Euler)

# Test expression is providing correct distribution
Vf = VectorFunctionSpace(mesh, 'CG', 1)
f = interpolate(t, Vf)

values = f.vector()
plt.figure()
plt.plot(values)

My main issue now is that I have to divide the traction expression by H*2 to retrieve the correct displacement. Any ideas why? This is a cause for concern when changing the geometry in these lines:

sparinner = Rectangle(Point(0, 0), Point(5000, 50))
sparouter = Rectangle(Point(5000, 0), Point(10000, 50))

As I’m not sure the same trick of dividing by 2*50 will still work for when spar inner and outer have different heights or when I have more complex geometry. Hence, it’d be great to get some understanding here.

If there are any other issues with the code i’d really appreciate extra advice also.

Note, I have kept the material properties the same for simplicity at this current state.

Thanks again.

17

First, you still assume a 3D beam model with out-of-plane width equal to H, this is why you must divide your traction by H. Just assume a traction per unit length along the out-of-plane direction and use Euler = (50e3)*Le**4/(30*E*((H**3)/12)) instead.
Second, now you changed your loading to a distributed loading on the exterior boundary but your ds measure corresponds to the total boundary. Hence, loading is applied on the left part (not important since displacement is fixed), on the right (not important either since loading is 0 on the right), and on both the top and bottom boundaries. This is why you need an additional factor 2 since loading is applied twice.
If you just want to apply it on the top boundary for instance do:

def top(x, on_boundary):
    return near(x[1], H) and on_boundary
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
AutoSubDomain(top).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)

the top facets are now marked as 1 and replace the ds measure by ds(1) in your linear form definition

18

Okay, that all seems to work. Thank you.

So if I was to have a more complex geometry of varying H, this expression would no longer be correct for the full beam as the out of plane width would also change and hence, couldn’t simply divide by H. How would I adjust his model for a varying H?

For example,

sparinner = Rectangle(Point(0, 0), Point(5000, 50))
sparouter = Rectangle(Point(5000, 0), Point(10000, 30))

So H0 = 50 and H1 = 30. I’m assuming there must be a way to incorporate this into the traction expression.

And I’ve already adjusted the boundary to,

# Define top boundary for traction
def top(x, on_boundary):
    return between(x[1], (25,50)) and on_boundary
facets = MeshFunction('size_t', mesh, 1)
facets.set_all(0)
AutoSubDomain(top).mark(facets, 1)
ds = Measure('ds', subdomain_data = facets)

t = Expression(('0.0', 'x[0] <= 5000 + tol ? ((-5*x[0]+50e3))/H0 : ((-5*x[0]+50e3))/H1'), H0=H0, H1=H1, tol=tol, element = el)

Does this seem correct?

Full code for 2 materials and 2 widths

from dolfin import *
from mshr import *
from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt

# Material properties
Le = 10000                      # Length mm
H = 50                          # Length mm 
nu = 0.3                        # Poisson's ratio = 0 for plane-stress conditions
model = 'plane_stress'

# Generate spar
sparinner = Rectangle(Point(0, 0), Point(5000, 50))
sparouter = Rectangle(Point(5000, 10), Point(10000, 40))
spar = sparinner + sparouter
mesh = generate_mesh(spar, 10)

# Create finite element function space, V (P2 Elements)
V = VectorFunctionSpace(mesh, 'CG', 2) 

# Create two materials of different Young's modulus throughout the beam using expressions
E = [100e9, 50e9]
mu_list = []
for x in E:
    mu_list.append(x/2/(1+nu))

lmbda_list = []
for x in E:
    lmbda_list.append(x*nu/(1+nu)/(1-2*nu))

lmbdaplane_list = []
for lmbda_value, mu_value in zip(lmbda_list, mu_list):
    lmbdaplane_list.append(2*mu_value*lmbda_value/(lmbda_value+2*mu_value))

if model == 'plane_stress':
    lmbda_list = lmbdaplane_list

tol = 1E-14
mu0 = mu_list[0]
mu1 = mu_list[1]
mu = Expression('x[0] <= 5000 + tol ? mu0 : mu1' , degree = 0, tol=tol, mu0=mu0, mu1=mu1)

tol = 1E-14
lmbda0 = lmbda_list[0]
lmbda1 = lmbda_list[1]
lmbda = Expression('x[0] <= 5000 + tol ? lmbda0 : lmbda1' , degree = 0, tol=tol, lmbda0=lmbda0, lmbda1=lmbda1)

# Define root boundary condition
def root_boundary(x, on_boundary):
    tol = 1e-14
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), root_boundary)

# Define top boundary for traction
def top(x, on_boundary):
    return between(x[1], (25,50)) and on_boundary
facets = MeshFunction('size_t', mesh, 1)
facets.set_all(0)
AutoSubDomain(top).mark(facets, 1)
ds = Measure('ds', subdomain_data = facets)

# Strain tensor
def eps(v):
    return sym(grad(v))

# Stress tensor 
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) 

# Define force distribution
el = V.ufl_element()
H0  = 50
H1 = 30 
t = Expression(('0.0', 'x[0] <= 5000 + tol ? ((-5*x[0]+50e3))/H0 : ((-5*x[0]+50e3))/H1'), H0=H0, H1=H1, tol=tol, element = el)   # Equation for triangular distributed load with root kN/mm = 50
                                                                   # Divide by H as traction assumes 3D model
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() 
v = TestFunction(V)
f = Constant((0, 0))
a = inner(sigma(u), eps(v))*dx
L = dot(f, v)*dx + dot(t, v)*ds(1)

# Compute solution 
u = Function(V)
solve(a == L, u, bc)

plot(u, title='Displacement', mode='displacement')

# Compute magnitude of displacement
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
     u_magnitude.vector().get_local().min(),
     u_magnitude.vector().get_local().max())

I’m now considering creating a model with several (say 10) point loads along the beam to represent the force distribution rather than have the expression. Do you have any guidance on where to start with this?

I’ve attempted to use,

# Apply point load at x = L
A, b = assemble_system(a, L, bc)
ps = PointSource(V, Point(10000, 50), 100e3)
ps.apply(b)

to apply a point load of 100kN at the end of the beam. But the resultant deflection is incorrect (for some reason out by a factor of exactly 50). And upon investigation, the point load also needs to be divided by H, does this also assume 3D?