study of cantilever beam

Binder

Consider a cantilever beam using the Finite Element Method library GetFEM++.

Creating a Model

Now we are ready, import the library.

[1]:
import getfem as gf
import numpy as np
import numpy.testing as npt
import pandas as pd
import pyvista as pv
from IPython.display import Markdown

from pyvirtualdisplay import Display
import os
display = Display(visible=0, size=(1280, 1024))
display.start()
[1]:
<pyvirtualdisplay.display.Display at 0x7ff96dbb9d90>

The review cases are as follows. GetFEM++ uses FEM_PRODUCT and IM_PRODUCT to two-dimensional the finite element method and the integration method, respectively. For quadratic elements, the Gaussian integration point is 3. IM_GAUSS1D(K) represents the integration point of \(K/2+1\) points. The element uses plane strain elements. Set up these meshes, the finite element method and the integration method.

[2]:
cases = [
    "case11",
    "case12",
    "case13",
    "case14",
    "case21",
    "case22",
    "case23",
    "case24",
    "case31",
    "case32",
    "case33",
    "case34",
    "case41",
    "case42",
    "case43",
    "case44",
]
[3]:
xs = [
    4,
    4,
    4,
    16,
    4,
    4,
    4,
    16,
    4,
    4,
    4,
    16,
    4,
    4,
    4,
    16,
]
ys = [
    1,
    2,
    4,
    8,
    1,
    2,
    4,
    8,
    1,
    2,
    4,
    8,
    1,
    2,
    4,
    8,
]
[4]:
fem_names = [
    "FEM_PK(1, 2)",
    "FEM_PK(1, 2)",
    "FEM_PK(1, 2)",
    "FEM_PK(1, 2)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK(1, 1)",
    "FEM_PK_WITH_CUBIC_BUBBLE(1, 1)",
    "FEM_PK_WITH_CUBIC_BUBBLE(1, 1)",
    "FEM_PK_WITH_CUBIC_BUBBLE(1, 1)",
    "FEM_PK_WITH_CUBIC_BUBBLE(1, 1)",
]
[5]:
methods = [
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(2)",
    "IM_GAUSS1D(2)",
    "IM_GAUSS1D(2)",
    "IM_GAUSS1D(2)",
    "IM_GAUSS1D(0)",
    "IM_GAUSS1D(0)",
    "IM_GAUSS1D(0)",
    "IM_GAUSS1D(0)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
    "IM_GAUSS1D(4)",
]
[6]:
pd.options.display.float_format = "{:.2f}".format
data = []
columns = ["Case Name", "Mesh", "Finite Element Method", "Integration Method"]
for case, x, y, fem_name, method in zip(cases, xs, ys, fem_names, methods):
    data.append([case, str(x) + "x" + str(y), fem_name, method])
df = pd.DataFrame(data=data, columns=columns)
Markdown(df.to_markdown())
[6]:
  Case Name Mesh Finite Element Method Integration Method
0 case11 4x1 FEM_PK(1, 2) IM_GAUSS1D(4)
1 case12 4x2 FEM_PK(1, 2) IM_GAUSS1D(4)
2 case13 4x4 FEM_PK(1, 2) IM_GAUSS1D(4)
3 case14 16x8 FEM_PK(1, 2) IM_GAUSS1D(4)
4 case21 4x1 FEM_PK(1, 1) IM_GAUSS1D(2)
5 case22 4x2 FEM_PK(1, 1) IM_GAUSS1D(2)
6 case23 4x4 FEM_PK(1, 1) IM_GAUSS1D(2)
7 case24 16x8 FEM_PK(1, 1) IM_GAUSS1D(2)
8 case31 4x1 FEM_PK(1, 1) IM_GAUSS1D(0)
9 case32 4x2 FEM_PK(1, 1) IM_GAUSS1D(0)
10 case33 4x4 FEM_PK(1, 1) IM_GAUSS1D(0)
11 case34 16x8 FEM_PK(1, 1) IM_GAUSS1D(0)
12 case41 4x1 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4)
13 case42 4x2 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4)
14 case43 4x4 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4)
15 case44 16x8 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4)

Mesh

The overall size of the model is L = 10 mm in length, h = 1 mm in height, and b = 1 mm in depth. In general, a slender ratio of 1: 10 is considered a beam element.

[7]:
L = 10.0
b = 1.0
h = 1.0
meshs = []
for case, x, y in zip(cases, xs, ys):
    X = np.arange(x + 1) * L / x
    Y = np.arange(y + 1) * h / y
    mesh = gf.Mesh("cartesian", X, Y)
    meshs.append(mesh)
    mesh.export_to_vtk("mesh_" + case + ".vtk", "ascii")

Outputs an image of each mesh.

[8]:
p = pv.Plotter(shape=(4, 4))

for i in range(4):
    for j in range(4):
        p.subplot(i, j)
        mesh = pv.read("mesh_" + cases[i * 4 + j] + ".vtk")
        p.add_text(cases[i * 4 + j], font_size=10)
        p.add_mesh(mesh, color="tan", show_edges=True)

p.show(cpos="xy")
_images/cantilever_13_0.png
[8]:
[(5.0, 0.5, 19.41486882686712),
 (5.0, 0.5, 0.0),
 (0.0, 1.0, 0.0)]

Region

Sets the area on the left side of the mesh where the Dirichlet condition is set. The right side sets the area for setting the Neumann condition.

[9]:
TOP_BOUND = 1
RIGHT_BOUND = 2
LEFT_BOUND = 3
BOTTOM_BOUND = 4

for mesh in meshs:
    fb1 = mesh.outer_faces_with_direction([0.0, 1.0], 0.01)
    fb2 = mesh.outer_faces_with_direction([1.0, 0.0], 0.01)
    fb3 = mesh.outer_faces_with_direction([-1.0, 0.0], 0.01)
    fb4 = mesh.outer_faces_with_direction([0.0, -1.0], 0.01)
    mesh.set_region(TOP_BOUND, fb1)
    mesh.set_region(RIGHT_BOUND, fb2)
    mesh.set_region(LEFT_BOUND, fb3)
    mesh.set_region(BOTTOM_BOUND, fb4)

Finite Element Method

Create a MeshFem object and associate the mesh with the finite element method.

[10]:
fems = []
for fem_name in fem_names:
    fems.append(gf.Fem("FEM_PRODUCT(" + fem_name + "," + fem_name + ")"))
[11]:
mfus = []
for mesh, fem in zip(meshs, fems):
    mfu = gf.MeshFem(mesh, 2)
    mfu.set_fem(fem)
    mfus.append(mfu)

Integral method

Associate the integration method with the mesh.

[12]:
ims = []
for method in methods:
    ims.append(gf.Integ("IM_PRODUCT(" + method + ", " + method + ")"))
[13]:
mims = []
for mesh, im in zip(meshs, ims):
    mim = gf.MeshIm(mesh, im)
    mims.append(mim)

Variable

Define the model object and set the variable “u”.

[14]:
mds = []
for mfu in mfus:
    md = gf.Model("real")
    md.add_fem_variable("u", mfu)
    mds.append(md)

Properties

Define properties as constants for the model object. The Young’s modulus of steel is \(E = 205000\times 10^6N/m^2\). Also set Poisson’s ratio \(\nu = 0.0\) to ignore the Poisson effect.

[15]:
E = 10000  # N/mm2
Nu = 0.0

for md in mds:
    md.add_initialized_data("E", E)
    md.add_initialized_data("Nu", Nu)

Plane Strain Element

Defines the plane strain element for variable ‘u’.

[16]:
for md, mim in zip(mds, mims):
    md.add_isotropic_linearized_elasticity_brick_pstrain(mim, "u", "E", "Nu")

Boundary Conditions

Set the Dirichlet condition for the region on the left side.

[17]:
for (md, mim, mfu, fem) in zip(mds, mims, mfus, fems):
    if fem.is_lagrange():
        md.add_Dirichlet_condition_with_simplification("u", LEFT_BOUND)
    else:
        md.add_Dirichlet_condition_with_multipliers(mim, "u", mfu, LEFT_BOUND)

Set the Neumann boundary condition on the right side.

[18]:
F = 1.0  # N/mm2
for (md, mfu, mim) in zip(mds, mfus, mims):
    md.add_initialized_data("F", [0, F / (b * h)])
    md.add_source_term_brick(mim, "u", "F", RIGHT_BOUND)

Solve

Solve the simultaneous equations of the model object to find the value of the variable ‘u’.

[19]:
for md in mds:
    md.solve()

The constraint on the left end has a displacement of 0.0.

[20]:
for md, mfu, case in zip(mds, mfus, cases):
    u = md.variable("u")
    dof = mfu.basic_dof_on_region(LEFT_BOUND)
    npt.assert_almost_equal(abs(np.max(u[dof])), 0.0)

Review results

Output and visualize the results of each case to vtk files.

[21]:
for md, mfu, case in zip(mds, mfus, cases):
    u = md.variable("u")
    mfu.export_to_vtk("u_" + case + ".vtk", "ascii", mfu, u, "u")

Computation of Theoretical Solutions

Calculates the deflection at each coordinate for comparison with the theoretical solution. The theoretical solution for a cantilever beam subjected to a concentrated load is as follows.

\(w\left(x\right)=\dfrac{FL^3}{3EI}\)

[22]:
I = b * h ** 3 / 12
w = F * L ** 3 / (3 * E * I)
w
[22]:
0.4

deformation diagram of case11

[23]:
p = pv.Plotter(shape=(4, 4))

for i in range(4):
    for j in range(4):
        p.subplot(i, j)
        mesh = pv.read("u_" + cases[i * 4 + j] + ".vtk")
        p.add_text(cases[i * 4 + j], font_size=10)
        p.add_mesh(mesh.warp_by_vector("u"), color="tan", show_edges=True)

p.show(cpos="xy")
_images/cantilever_41_0.png
[23]:
[(5.015016078948975, 0.7012025117874146, 19.565021099812608),
 (5.015016078948975, 0.7012025117874146, 0.0),
 (0.0, 1.0, 0.0)]

Comparison with Theoretical Solutions

[24]:
b = 1.0
I = b * h ** 3.0 / 12.0
dmax = 1.0 / 3.0 * (F * L ** 3) / (E * I)
dmax
[24]:
0.4

Calculates the ratio of deformation to theoretical solution for each case.

[25]:
data = []

columns = [
    "case name",
    "mesh",
    "finite element method",
    "integration method",
    "ratio with theory",
]

for case, x, y, fem_name, method, md, mfu in zip(
    cases, xs, ys, fem_names, methods, mds, mfus
):
    u = md.variable("u")
    dof = mfu.basic_dof_on_region(RIGHT_BOUND)
    ratio = "{:.3f}".format(max(u[dof] / dmax))
    data.append([case, str(x) + "x" + str(y), fem_name, method, ratio])
df = pd.DataFrame(data=data, columns=columns)
Markdown(df.to_markdown())
[25]:
  case name mesh finite element method integration method ratio with theory
0 case11 4x1 FEM_PK(1, 2) IM_GAUSS1D(4) 0.999
1 case12 4x2 FEM_PK(1, 2) IM_GAUSS1D(4) 1
2 case13 4x4 FEM_PK(1, 2) IM_GAUSS1D(4) 1
3 case14 16x8 FEM_PK(1, 2) IM_GAUSS1D(4) 1.006
4 case21 4x1 FEM_PK(1, 1) IM_GAUSS1D(2) 0.244
5 case22 4x2 FEM_PK(1, 1) IM_GAUSS1D(2) 0.244
6 case23 4x4 FEM_PK(1, 1) IM_GAUSS1D(2) 0.244
7 case24 16x8 FEM_PK(1, 1) IM_GAUSS1D(2) 0.841
8 case31 4x1 FEM_PK(1, 1) IM_GAUSS1D(0) 6.52355e+10
9 case32 4x2 FEM_PK(1, 1) IM_GAUSS1D(0) 1.317
10 case33 4x4 FEM_PK(1, 1) IM_GAUSS1D(0) 1.056
11 case34 16x8 FEM_PK(1, 1) IM_GAUSS1D(0) 1.021
12 case41 4x1 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4) 0.999
13 case42 4x2 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4) 1
14 case43 4x4 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4) 1
15 case44 16x8 FEM_PK_WITH_CUBIC_BUBBLE(1, 1) IM_GAUSS1D(4) 1.006
[ ]: