# Wrong max displacement | Linear Elasticity

**URL:** <https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728>\
**Category:** General\
**Created:** [November 11, 2023, 9:28pm UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728 "2023-11-11T21:28:06Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![tarekhbb](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tarekhbb/32/5854_2.png) [@tarekhbb](https://fenicsproject.discourse.group/u/tarekhbb)\
**Post date:** [November 11, 2023, 9:28pm UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/1 "2023-11-11T21:28:07Z")

</div>

Hello, I’m trying to do a project for my class. Basically I’m trying to test this program to calculate the max displacement in a beam however the value it exports is different to the one that appears in Paraview and different to the one I calculated.  
Max displacements:  
Code: 0.0049421492043371365m  
Manual: 5.08e-4  
Paraview: 2.9e-3

Here’s the code:  
from fenics import \*  
import numpy as np

L = 10 #Longitud de viga en metros  
H = 0.5  
W = 0.3

malla = BoxMesh(Point(0, 0, 0), Point(L, H, W), 100, 10, 10)

V = VectorFunctionSpace(malla, ‘P’, 1)

def frontera(x, on\_boundary):  
return on\_boundary and near(x[0], 0)

cf = DirichletBC(V, Constant((0, 0, 0)), frontera)

E = 210e9 # Módulo de Young  
nu = 0.3 # Relación de Poisson  
mu = E / (2 \* (1 + nu))  
lambda\_ = E \* nu / ((1 + nu) \* (1 - 2 \* nu))

u = TrialFunction(V)  
d = u.geometric\_dimension()  
I = Identity(d)  
epsilon = lambda u: 0.5 \* (grad(u) + grad(u).T)  
sigma = lambda u: lambda\_ \* tr(epsilon(u)) \* I + 2 \* mu \* epsilon(u)

v = TestFunction(V)  
T = Constant((0, 1000, 0)) # Carga puntual  
a = inner(sigma(u), epsilon(v)) \* dx  
L = dot(T, v) \* ds

u = Function(V)  
solve(a == L, u, cf)  
vtkfile = File(‘desplazamiento.pvd’)  
vtkfile \<\< u

displacement\_values = u.compute\_vertex\_values(malla)

d = u.geometric\_dimension()

displacement\_magnitude = np.sqrt(sum(displacement\_values[i::d]\*\*2 for i in range(d)))

max\_displacement = np.max(displacement\_magnitude)  
print(“Desplazamiento máximo:”, max\_displacement)

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [November 12, 2023, 4:59am UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/2 "2023-11-12T04:59:45Z")

</div>

You are making mistake in calculation of the max displacement. The maximum of the displacement obviously happens at the free end of the beam and is equal to 0.0028940 as shown here:

 ![Screenshot 2023-11-11 at 10.44.46 PM](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/bd8f07f5ea2d5c631f0f85fd012a16321d9eabe0.png)

The maximum of the displacement could be found like this:

```auto
from fenics import *
import numpy as np

L = 10 #Longitud de viga en metros
H = 0.5
W = 0.3

malla = BoxMesh(Point(0, 0, 0), Point(L, H, W), 100, 10, 10)

V = VectorFunctionSpace(malla, 'P', 1)

def frontera(x, on_boundary):
    return on_boundary and near(x[0], 0)

cf = DirichletBC(V, Constant((0, 0, 0)), frontera)

E = 210e9 # Módulo de Young
nu = 0.3 # Relación de Poisson
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))

u = TrialFunction(V)
d = u.geometric_dimension()
I = Identity(d)
epsilon = lambda u: 0.5 * (grad(u) + grad(u).T)
sigma = lambda u: lambda_ * tr(epsilon(u)) * I + 2 * mu * epsilon(u)

v = TestFunction(V)
T = Constant((0, 1000, 0)) # Carga puntual
a = inner(sigma(u), epsilon(v)) * dx
L = dot(T, v) * ds

u = Function(V)
solve(a == L, u, cf)
vtkfile = File('desplazamiento.pvd')
vtkfile << u

#displacement_values = u.compute_vertex_values(malla)

#d = u.geometric_dimension()

#displacement_magnitude = np.sqrt(sum(displacement_values[i::d]**2 for i in range(d)))

#max_displacement = np.max(displacement_magnitude)
#print('Desplazamiento máximo:', max_displacement)

#print (u(10,0,0))

x = malla.coordinates()

All_Nodes_Displacements_XYZ = []
Final_Calculated_Displacement = []

for i in range(len(x)):
    All_Nodes_Displacements_XYZ.append(u(x[i][0],x[i][1],x[i][2]))

for j in range (len(All_Nodes_Displacements_XYZ)):
    Final_Calculated_Displacement.append(np.sqrt(pow(All_Nodes_Displacements_XYZ[j][0],2) +\
                                          pow(All_Nodes_Displacements_XYZ[j][1],2) +\
                                          pow(All_Nodes_Displacements_XYZ[j][2],2)))

print ('The MAX displacement is: ',np.max(Final_Calculated_Displacement))

```

The output of the above code is:

```auto
The MAX displacement is: 0.0028940563971543896

```

Which is exactly the same value you can visualize in Paraview.

---

<div class="post-metadata">

**Author:** ![tarekhbb](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tarekhbb/32/5854_2.png) [@tarekhbb](https://fenicsproject.discourse.group/u/tarekhbb)\
**Post date:** [November 12, 2023, 5:29am UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/3 "2023-11-12T05:29:07Z")

</div>

Omg you’re right! But still the max displacement doesn’t coincide with the max displacement calculated from:

```auto
max_displacement=(P*L^3)/(3*E*I)

```

Which using my values would be

```auto
max_displacement=(1000*10^3)/(3*210e9*3.125e-3) = 5.08e-4

```

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [November 13, 2023, 2:00am UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/4 "2023-11-13T02:00:37Z")

</div>

First of all if you are assuming a clamped beam under uniform distributed load, the maximum displacement is calculated as: \delta\_{max} = {wL^4}/{8EI}. You may want to check [Beam Deflection Formulas](https://www.scribd.com/document/398760596/Beam-Deflection-Formula).  
Note that the \omega should be multiplied by the area of the top face of the beam so its unit is consistent to N/m.  
That said, we have, \delta\_{max} = {w (b h) l^4}/{8E(\frac {1}{12}bh^3)}  
In addition, the linear part of the variational form (L) should be multiplied by dx.  
Here is how it works:

```auto
from fenics import *
import numpy as np

Length = 10 #Longitud de viga en metros
H = 0.5
W = 0.3

malla = BoxMesh(Point(0, 0, 0), Point(Length, H, W), 100, 10, 10)

V = VectorFunctionSpace(malla, 'Lagrange', degree=1)

def frontera(x, on_boundary):
    return on_boundary and near(x[0], 0)

cf = DirichletBC(V, Constant((0, 0, 0)), frontera)

E = 210e9 # Módulo de Young
nu = nu = Constant(0.3) # Relación de Poisson

mu = E / (2 * (1 + nu))
lambda_ = E*nu/(1-nu**2)
lmbda = 2*mu*lambda_/(lambda_+2*mu)

u = TrialFunction(V)
d = u.geometric_dimension()
I = Identity(d)
epsilon = lambda u: 0.5 * (grad(u) + grad(u).T)
sigma = lambda u: lambda_ * tr(epsilon(u)) * I + 2 * mu * epsilon(u)

v = TestFunction(V)
T = Constant((0, 1000, 0)) # Carga puntual
a = inner(sigma(u), epsilon(v)) * dx
L = inner(T, v) * dx

u = Function(V)
solve(a == L, u, cf)
vtkfile = File('desplazamiento.pvd')
vtkfile << u

x = malla.coordinates()

All_Nodes_Displacements_XYZ = []
Final_Calculated_Displacement = []

for i in range(len(x)):
    All_Nodes_Displacements_XYZ.append(u(x[i][0],x[i][1],x[i][2]))

for j in range (len(All_Nodes_Displacements_XYZ)):
    Final_Calculated_Displacement.append(np.sqrt(pow(All_Nodes_Displacements_XYZ[j][0],2) +\
                                          pow(All_Nodes_Displacements_XYZ[j][1],2) +\
                                          pow(All_Nodes_Displacements_XYZ[j][2],2)))

print ('The MAX displacement from FEM is: ', np.max(Final_Calculated_Displacement))

force = 1000 * H * W

print ("The MAX displacement from theory is: ", force * pow(Length,4) / (8 * E * W * pow(H,3)/12.))

```

The output of the above code is:

```auto
The MAX displacement from FEM is: 0.00028106529009377603
The MAX displacement from theory is: 0.00028571428571428574

```

---

<div class="post-metadata">

**Author:** ![tarekhbb](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tarekhbb/32/5854_2.png) [@tarekhbb](https://fenicsproject.discourse.group/u/tarekhbb)\
**Post date:** [November 14, 2023, 1:43am UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/5 "2023-11-14T01:43:19Z")

</div>

I’m working with a Point Load at the end of the beam, not a distributed load. I tried to point it out with the variable T = Constant((0, 1000, 0))

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [November 14, 2023, 2:56am UTC](https://fenicsproject.discourse.group/t/wrong-max-displacement-linear-elasticity/12728/6 "2023-11-14T02:56:10Z")

</div>

> I’m working with a Point Load at the end of the beam, not a distributed load. I tried to point it out with the variable T = Constant((0, 1000, 0))

In your initial code, you have never referred to a coordinate where the end load should be applied to. The T defines the vector of the distributed load. For implementation of a clamped beam under end load you may want to have a look at [this discussion](https://fenicsproject.discourse.group/t/how-can-point-loads-be-applied-to-the-end-of-the-cantilever-beam/8253/6)