Wrong max displacement | Linear Elasticity

General
1

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)

2

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
Screenshot 2023-11-11 at 10.44.46 PM2206×1204 305 KB

The maximum of the displacement could be found like this:

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:

The MAX displacement is:  0.0028940563971543896

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

1 Like
3

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

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

Which using my values would be

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

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 .
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:

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:

The MAX displacement from FEM is:  0.00028106529009377603
The MAX displacement from theory is:  0.00028571428571428574
2 Likes
5

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))

6

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