Examples¶

exs_spring¶

Purpose:

Show the basic steps in a finite element calculation.

Description:

The general procedure in linear finite element calculations is carried out for a simple structure. The steps are:

• Define the model

• Generate element matrices

• Assemble element matrices into the global system of equations

• Solve the global system of equations

• Evaluate element forces

Consider the system of three linear elastic springs, and the corresponding finite element model. The system of springs is fixed in its ends and loaded by a single load \(F=100 N\) and a stiffness \(k=1500 N/m\)

Example:

The computation is initialized by importing CALFEM and NumPy, then defining the topology matrix edof containing the degrees of freedom for each element:

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu

edof = np.array([
    [1, 2],                 # element 1 between nodes 1 and 2
    [2, 3],                 # element 2 between nodes 2 and 3
    [2, 3]                  # element 3 between nodes 2 and 3
])

The global stiffness matrix K (3×3) of zeros:

K = np.zeros((3, 3))        # empty global stiffness matrix
f = np.zeros((3, 1))        # empty global load vector

And the load vector f (3×1) with the load \(F=100\) at DOF 2:

f[1] = 100.0                # (N), f[1] corresponds to dof 2

Element stiffness matrices are generated by the function cfc.spring1e. The element property ep for the springs contains the spring stiffnesses \(k\) and \(2k\) respectively, where \(k=1500\):

CodeOutput

k = 1500.                   # (N/m)
ep1 = 2.*k                  # spring stiffness, element 1
ep2 = k                     # spring stiffness, element 2
ep3 = 2.*k                  # spring stiffness, element 3
Ke1 = cfc.spring1e(ep1)     # element stiffness matrix, element 1
Ke2 = cfc.spring1e(ep2)     # element stiffness matrix, element 2
Ke3 = cfc.spring1e(ep3)     # element stiffness matrix, element 3

cfu.disp_h2("Ke1")
cfu.disp_array(Ke1, tablefmt='plain')
cfu.disp_h2("Ke2")
cfu.disp_array(Ke2, tablefmt='plain')
cfu.disp_h2("Ke3")
cfu.disp_array(Ke3, tablefmt='plain')

The element stiffness matrices are assembled into the global stiffness matrix K according to the topology:

CodeOutput

cfc.assem(edof[0, :], K, Ke1)   # assemble element stiffness matrix 1
cfc.assem(edof[1, :], K, Ke2)   # assemble element stiffness matrix 2
cfc.assem(edof[2, :], K, Ke3)   # assemble element stiffness matrix 3

cfu.disp_h2("Stiffness matrix K (N/m):")
cfu.disp_array(K, tablefmt='plain')

The global system of equations is solved considering the boundary conditions. DOFs 1 and 3 are fixed (value 0):

CodeOutput

bc_dofs = np.array([1, 3])       # dofs with prescribed displacment (0)
bc_vals = np.array([0.0, 0.0])   # prescribed displacments

a, r = cfc.solveq(K, f, bc_dofs, bc_vals)

cfu.disp_h2("Displacements a (m):")
cfu.disp_array(a, tablefmt="plain")

cfu.disp_h2("Reaction forces r (N):")
cfu.disp_array(r, tablefmt="plain")

Element forces are evaluated from the element displacements. These are obtained from the global displacements a using the function cfc.extract_ed and the spring forces are evaluated using the function cfc.spring1s:

CodeOutput

ed1 = cfc.extract_ed(edof[0, :], a) # Nodal displacements, element 1
ed2 = cfc.extract_ed(edof[1, :], a) # Nodal displacements, element 2
ed3 = cfc.extract_ed(edof[2, :], a) # Nodal displacements, element 3

es1 = cfc.spring1s(ep1, ed1)        # Element force, element 1
es2 = cfc.spring1s(ep2, ed2)        # Element force, element 2
es3 = cfc.spring1s(ep3, ed3)        # Element force, element 3

cfu.disp_h2("Element forces (N):")
print("N1 = "+str(es1))
print("N2 = "+str(es2))
print("N3 = "+str(es3))

exs_flw_temp1¶

Purpose:

Analysis of one-dimensional heat flow.

Description:

Consider a wall built up of concrete and thermal insulation. The outdoor temperature is \(-17°\text{C}\) and the temperature inside is \(20°\text{C}\). At the inside of the thermal insulation there is a heat source yielding \(10\) W/m².

The wall is subdivided into five elements and the one-dimensional spring (analogy) element spring1e is used. Equivalent spring stiffnesses are \(k_i=\lambda A/L\) for thermal conductivity and \(k_i=A/R\) for thermal surface resistance. Corresponding spring stiffnesses per m² of the wall are:

\(k_1 =\)	\(1/0.04\)	\(=\)	\(25.0\)	W/K
\(k_2 =\)	\(1.7/0.070\)	\(=\)	\(24.3\)	W/K
\(k_3 =\)	\(0.040/0.100\)	\(=\)	\(0.4\)	W/K
\(k_4 =\)	\(1.7/0.100\)	\(=\)	\(17.0\)	W/K
\(k_5 =\)	\(1/0.13\)	\(=\)	\(7.7\)	W/K

Example:

The computation is initialized by importing CALFEM and NumPy. A global system matrix K and a heat flow vector f are defined. The heat source inside the wall is considered by setting \(f_4=10\). The element matrices Ke are computed using cfc.spring1e, and the function cfc.assem assembles the global stiffness matrix. The system of equations is solved using cfc.solveq with boundary conditions. The prescribed temperatures are \(T_1=-17°\text{C}\) and \(T_6=20°\text{C}\):

CodeOutput

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu

# ----- Topology -------------------------------------------------

edof = np.array([
    [1, 2],
    [2, 3],
    [3, 4],
    [4, 5],
    [5, 6]
])

# ----- Stiffness matrix K and load vector f ---------------------

K = np.array(np.zeros((6, 6)))
f = np.array(np.zeros((6, 1)))
f[3] = 10.0

# ----- Element stiffness and element load matrices  -------------

ep1 = 25
ep2 = 24.3
ep3 = 0.4
ep4 = 17
ep5 = 7.7

Ke1 = cfc.spring1e(ep1)
Ke2 = cfc.spring1e(ep2)
Ke3 = cfc.spring1e(ep3)
Ke4 = cfc.spring1e(ep4)
Ke5 = cfc.spring1e(ep5)

# ----- Assemble Ke into K ---------------------------------------

cfc.assem(edof[0, :], K, Ke1)
cfc.assem(edof[1, :], K, Ke2)
cfc.assem(edof[2, :], K, Ke3)
cfc.assem(edof[3, :], K, Ke4)
cfc.assem(edof[4, :], K, Ke5)

# ----- Solve the system of equations ----------------------------

bc = np.array([1, 6])
bcVal = np.array([-17, 20])
a, r = cfc.solveq(K, f, bc, bcVal)

cfu.disp_h2("Temperatures a:")
cfu.disp_array(a, tablefmt="plain")

cfu.disp_h2("Reaction flows r:")
cfu.disp_array(r, tablefmt="plain")

The temperature values \(a_i\) at the node points are given in the vector a and the boundary heat flows in the vector r.

After solving the system of equations, the heat flow through each element is computed using cfc.extract_ed and cfc.spring1s:

CodeOutput

ed1 = cfc.extract_ed(edof[0, :], a)
ed2 = cfc.extract_ed(edof[1, :], a)
ed3 = cfc.extract_ed(edof[2, :], a)
ed4 = cfc.extract_ed(edof[3, :], a)
ed5 = cfc.extract_ed(edof[4, :], a)

q1 = cfc.spring1s(ep1, ed1)
q2 = cfc.spring1s(ep2, ed2)
q3 = cfc.spring1s(ep3, ed3)
q4 = cfc.spring1s(ep4, ed4)
q5 = cfc.spring1s(ep5, ed5)

cfu.disp_h2("Element flows:")

print("q1 = ")
print(q1)
print("q2 = ")
print(q2)
print("q3 = ")
print(q3)
print("q4 = ")
print(q4)
print("q5 = ")
print(q5)

The heat flow through the wall is \(q=14.0\) W/m² in the part of the wall to the left of the heat source (elements 1-3), and \(q=4.0\) W/m² in the part to the right of the heat source (elements 4-5). This demonstrates how the internal heat source affects the heat flow distribution through the wall.

exs_bar2¶

Purpose:

Analysis of a plane truss.

Description:

Consider a plane truss consisting of three bars with the properties \(E=200\) GPa, \(A_1=6.0 \times 10^{-4}\) m², \(A_2=3.0 \times 10^{-4}\) m² and \(A_3=10.0 \times 10^{-4}\) m², and loaded by a single force \(P=80\) kN. The corresponding finite element model consists of three elements and eight degrees of freedom.

Example:

The computation is initialized by importing CALFEM and NumPy and defining the element topology matrix containing only the degrees of freedom for each element:

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu
import calfem.vis_mpl as cfv

edof = np.array([
    [1, 2, 5, 6],
    [5, 6, 7, 8],
    [3, 4, 5, 6]
])

The global stiffness matrix K and the load vector f are initialized. The load of 80 kN is applied in the negative y-direction at DOF 6:

K = np.array(np.zeros((8, 8)))
f = np.array(np.zeros((8, 1)))

f[5] = -80e3

The element property vectors ep1, ep2 and ep3 are defined by

E = 2.0e11          # (N/m^2)
A1 = 6.0e-4         # (m^2)
A2 = 3.0e-4         # (m^2)
A3 = 10.0e-4        # (m^2)
ep1 = [E, A1]
ep2 = [E, A2]
ep3 = [E, A3]

and the element coordinate vectors ex1, ex2, ex3, ey1, ey2 and ey3 by

ex1 = np.array([0.0, 1.6])
ex2 = np.array([1.6, 1.6])
ex3 = np.array([0.0, 1.6])

ey1 = np.array([0.0, 0.0])
ey2 = np.array([0.0, 1.2])
ey3 = np.array([1.2, 0.0])

The element stiffness matrices Ke1, Ke2 and Ke3 are computed using bar2e:

CodeOutput

Ke1 = cfc.bar2e(ex1, ey1, ep1)
Ke2 = cfc.bar2e(ex2, ey2, ep2)
Ke3 = cfc.bar2e(ex3, ey3, ep3)

Ke1 = cfc.bar2e(ex1, ey1, ep1)
Ke2 = cfc.bar2e(ex2, ey2, ep2)
Ke3 = cfc.bar2e(ex3, ey3, ep3)

cfu.disp_h2("Ke1")
cfu.disp_array(Ke1, tablefmt='plain')
cfu.disp_h2("Ke2")
cfu.disp_array(Ke2, tablefmt='plain')
cfu.disp_h2("Ke3")
cfu.disp_array(Ke3, tablefmt='plain')

Based on the topology information, the global stiffness matrix can be generated by assembling the element stiffness matrices.

CodeOutput

K = cfc.assem(edof[0, :], K, Ke1)
K = cfc.assem(edof[1, :], K, Ke2)
K = cfc.assem(edof[2, :], K, Ke3)

cfu.disp_h2("Stiffness matrix K (N/m):")
cfu.disp_array(K, tablefmt='plain')

The system of equations \(\boldsymbol{Ka}=\boldsymbol{f}\) is solved by specifying the boundary conditions and using solveq():

CodeOutput

bc_dofs = np.array([1, 2, 3, 4, 7, 8])
bc_vals = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
a, r = cfc.solveq(K, f, bc_dofs, bc_vals)

cfu.disp_h2("Displacements a (m):")
cfu.disp_array(a)

cfu.disp_h2("Reaction forces r (N):")
cfu.disp_array(r)

The vertical displacement at the point of loading is 1.15 mm. The element forces are calculated by extracting element displacements and computing stresses:

CodeOutput

ed1 = cfc.extract_ed(edof[0, :], a)
ed2 = cfc.extract_ed(edof[1, :], a)
ed3 = cfc.extract_ed(edof[2, :], a)

es1 = cfc.bar2s(ex1, ey1, ep1, ed1)
es2 = cfc.bar2s(ex2, ey2, ep2, ed2)
es3 = cfc.bar2s(ex3, ey3, ep3, ed3)

cfu.disp_h2("Element normal forces (N):")

cfu.disp("N1")
cfu.disp_array(es1)
cfu.disp("N2")
cfu.disp_array(es2)
cfu.disp("N3")
cfu.disp_array(es3)

The normal forces are \(N_1=-29.84\) kN (compression), \(N_2=57.62\) kN (tension) and \(N_3=37.31\) kN (tension).

Visualization of the results can be performed using CALFEM’s visualization functions:

plotpar1 = [2, 1, 0]
sfac = cfv.scalfact2(ex3, ey3, ed1, 0.1)
print("sfac =", sfac)

cfv.figure(1)
cfv.eldraw2(ex1, ey1, plotpar1)
cfv.eldraw2(ex2, ey2, plotpar1)
cfv.eldraw2(ex3, ey3, plotpar1)

plotpar2 = [1, 2, 1]
cfv.eldisp2(ex1, ey1, ed1, plotpar2, sfac)
cfv.eldisp2(ex2, ey2, ed2, plotpar2, sfac)
cfv.eldisp2(ex3, ey3, ed3, plotpar2, sfac)
cfv.axis([-0.4, 2.0, -0.4, 1.4])
plotpar3 = 2
cfv.scalgraph2(sfac, [1e-3, 0, -0.3], plotpar3)
cfv.title("Displacements")

# ----- Draw normal force diagram --------------------------------

plotpar = [2, 1]
sfac = cfv.scalfact2(ex1, ey1, es2[:, 0], 0.1)
cfv.figure(2)
cfv.secforce2(ex1, ey1, es1[:, 0], plotpar, sfac)
cfv.secforce2(ex2, ey2, es2[:, 0], plotpar, sfac)
cfv.secforce2(ex3, ey3, es3[:, 0], plotpar, sfac)
cfv.axis([-0.4, 2.0, -0.4, 1.4])
cfv.scalgraph2(sfac, [5e4, 0, -0.3], plotpar3)
cfv.title("Normal force")

cfv.show_and_wait()

Figures:

exs_bar2_l¶

Purpose:

Analysis of a plane truss.

Description:

Consider a plane truss, loaded by a single force \(P=0.5\) MN.

The corresponding finite element model consists of ten elements and twelve degrees of freedom.

Material properties:

• Cross-sectional area: \(A=25.0 \times 10^{-4}\) m²

• Young’s modulus: \(E=2.10 \times 10^{5}\) MPa

Example:

The computation is initialized by importing CALFEM and NumPy. The element topology matrix contains only the degrees of freedom for each element:

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu

edof = np.array([
    [1, 2, 5, 6],
    [3, 4, 7, 8],
    [5, 6, 9, 10],
    [7, 8, 11, 12],
    [7, 8, 5, 6],
    [11, 12, 9, 10],
    [3, 4, 5, 6],
    [7, 8, 9, 10],
    [1, 2, 7, 8],
    [5, 6, 11, 12],
])

The global stiffness matrix K and load vector f are initialized. The load \(P=0.5\) MN is divided into x and y components:

K = np.zeros([12, 12])
f = np.zeros([12, 1])
f[10] = 0.5e6 * np.sin(np.pi / 6)
f[11] = -0.5e6 * np.cos(np.pi / 6)

The material and geometric properties are defined, along with element coordinate matrices:

A = 25.0e-4
E = 2.1e11
ep = [E, A]

ex = np.array([
    [0.0, 2.0],
    [0.0, 2.0],
    [2.0, 4.0],
    [2.0, 4.0],
    [2.0, 2.0],
    [4.0, 4.0],
    [0.0, 2.0],
    [2.0, 4.0],
    [0.0, 2.0],
    [2.0, 4.0],
])

ey = np.array([
    [2.0, 2.0],
    [0.0, 0.0],
    [2.0, 2.0],
    [0.0, 0.0],
    [0.0, 2.0],
    [0.0, 2.0],
    [0.0, 2.0],
    [0.0, 2.0],
    [2.0, 0.0],
    [2.0, 0.0],
])

The element stiffness matrices are computed and assembled in a loop:

for elx, ely, eltopo in zip(ex, ey, edof):
    Ke = cfc.bar2e(elx, ely, ep)
    cfc.assem(eltopo, K, Ke)

The system of equations is solved by specifying boundary conditions and using solveq():

CodeOutput

bc_dofs = np.array([1, 2, 3, 4])
bc_vals = np.array([0.0, 0.0, 0.0, 0.0])
a, r = cfc.solveq(K, f, bc_dofs, bc_vals)

cfu.disp_h2("Displacements a:")
cfu.disp_array(a, tablefmt='plain')

cfu.disp_h2("Reaction forces r:")
cfu.disp_array(r, tablefmt='plain')

The displacement at the point of loading is \(-1.7 \times 10^{-3}\) m in the x-direction and \(-11.3 \times 10^{-3}\) m in the y-direction. At the upper support the horizontal force is \(-0.866\) MN and the vertical \(0.240\) MN. At the lower support the forces are \(0.616\) MN and \(0.193\) MN, respectively.

Normal forces are evaluated from element displacements. These are obtained from the global displacements using extract_ed() and the forces are calculated using bar2s():

CodeOutput

ed = cfc.extract_ed(edof, a)
N = np.zeros([edof.shape[0]])

cfu.disp_h2("Element forces:")

i = 0
for elx, ely, eld in zip(ex, ey, ed):
    es = cfc.bar2s(elx, ely, ep, eld)
    N[i] = es[0][0]
    print("N%d = %g" % (i + 1, N[i]))
    i += 1

The largest normal force \(N=0.626\) MN is obtained in element 1 and is equivalent to a normal stress \(\sigma=250\) MPa.

To reduce the quantity of input data, the element coordinate matrices ex and ey can alternatively be created from a global coordinate matrix coord and a global topology matrix dof using coord_extract():

coord = np.array([
    [0, 2],
    [0, 0],
    [2, 2],
    [2, 0],
    [4, 2],
    [4, 0]
])

dof = np.array([
    [1, 2],
    [3, 4],
    [5, 6],
    [7, 8],
    [9, 10],
    [11, 12]
])

ex, ey = cfc.coord_extract(edof, coord, dof, 2)

exs_beam1¶

Purpose:

Analysis of a simply supported beam.

Description:

Consider a beam with the length \(9.0\) m. The beam is simply supported and loaded by a point load \(P=10000\) N applied at a point \(3.0\) m from the left support. The corresponding computational model has six degrees of freedom and consists of two beam elements with four degrees of freedom. The beam has Young’s modulus \(E=210\) GPa and moment of inertia \(I=2510 \times 10^{-8}\) m⁴.

Example:

The computation is initialized by importing CALFEM and NumPy. The element topology matrix contains only the degrees of freedom for each element:

import numpy as np
import matplotlib.pyplot as plt
import calfem.core as cfc
import calfem.utils as cfu
import calfem.vis_mpl as cfv

edof = np.array([
    [1, 2, 3, 4],
    [3, 4, 5, 6]
])

The global stiffness matrix K and load vector f are initialized. The point load P = 10000 N is applied at DOF 3:

K = np.array(np.zeros((6, 6)))
f = np.array(np.zeros((6, 1)))
f[2] = -10e3

The material and geometric properties are defined, along with element coordinates and stiffness matrices:

CodeOutput

E = 210e9
I = 2510e-8

ep = np.array([E, I])
ex1 = np.array([0, 3])
ex2 = np.array([3, 9])
eq1 = np.array([0])
eq2 = np.array([0])

Ke1 = cfc.beam1e(ex1, ep)
Ke2 = cfc.beam1e(ex2, ep)

cfu.disp_h2("Ke1")
cfu.disp_array(Ke1, tablefmt='plain')
cfu.disp_h2("Ke2")
cfu.disp_array(Ke2, tablefmt='plain')

Based on the topology information, the global stiffness matrix can be generated by assembling the element stiffness matrices.

cfc.assem(edof[0, :], K, Ke1)
cfc.assem(edof[1, :], K, Ke2)

The system of equations is solved by defining boundary conditions and using solveq():

CodeOutput

bc_dofs = np.array([1, 5])
bc_vals = np.array([0.0, 0.0])
a, r = cfc.solveq(K, f, bc_dofs, bc_vals)

cfu.disp_array(a, ["a"])
cfu.disp_array(r, ["r"])

The section forces and element displacements are calculated from global displacements:

CodeOutput

ed = cfc.extract_ed(edof, a)

es1, ed1, ec1 = cfc.beam1s(ex1, ep, ed[0, :], eq1, nep=4)
es2, ed2, ec2 = cfc.beam1s(ex2, ep, ed[1, :], eq2, nep=7)

cfu.disp_h2("es1")
cfu.disp_array(es1, ["V1", "M1"])
cfu.disp_h2("ed1")
cfu.disp_array(ed1, ["v1"])
cfu.disp_h2("ec1")
cfu.disp_array(ec1, ["x1"])
cfu.disp_h2("es2")
cfu.disp_array(es2, ["V2", "M2"])
cfu.disp_h2("ed2")
cfu.disp_array(ed2, ["v2"])
cfu.disp_h2("ec2")
cfu.disp_array(ec2, ["x2"])

Results:

The solution shows the expected behavior for a simply supported beam with a point load:

• Maximum deflection: 22.8 mm at the loading point (DOF 3)

• Support reactions: 6667 N at left support, 3333 N at right support

• Maximum shear force: 6667 N (in element 1)

• Maximum moment: 20000 Nm (at the loading point)

Visualization of the results using matplotlib:

cfv.figure(1)
plt.plot([0, 9], [0, 0], color=(0.8, 0.8, 0.8))
plt.plot(
    np.concatenate(([0], ec1[:, 0], 3 + ec2[:, 0], [9]), 0),
    np.concatenate(([0], ed1[:, 0], ed2[:, 0], [0]), 0),
    color=(0.0, 0.0, 0.0),
)
cfv.title("displacements")

# ----- Draw shear force diagram----------------------------------

cfv.figure(2)
plt.plot([0, 9], [0, 0], color=(0.8, 0.8, 0.8))
plt.plot(
    np.concatenate(([0], ec1[:, 0], 3 + ec2[:, 0], [9]), 0),
    -np.concatenate(([0], es1[:, 0], es2[:, 0], [0]), 0) / 1000,
    color=(0.0, 0.0, 0.0),
)
cfv.title("shear force")

# ----- Draw moment diagram----------------------------------

cfv.figure(3)
plt.plot([0, 9], [0, 0], color=(0.8, 0.8, 0.8))
plt.plot(
    np.concatenate(([0], ec1[:, 0], 3 + ec2[:, 0], [9]), 0),
    -np.concatenate(([0], es1[:, 1], es2[:, 1], [0]), 0) / 1000,
    color=(0.0, 0.0, 0.0),
)
cfv.title("bending moment")

cfv.show_and_wait()

exs_beam2¶

Purpose:

Analysis of a plane frame.

Description:

A frame consists of one horizontal and two vertical beams according to the figure.

Material and geometric properties:

\(E\)	\(=\)	\(200\) GPa
\(A_1\)	\(=\)	\(2.0 \times 10^{-3}\) m²
\(I_1\)	\(=\)	\(1.6 \times 10^{-5}\) m⁴
\(A_2\)	\(=\)	\(6.0 \times 10^{-3}\) m²
\(I_2\)	\(=\)	\(5.4 \times 10^{-5}\) m⁴
\(P\)	\(=\)	\(2.0\) kN
\(q_0\)	\(=\)	\(10.0\) kN/m

The corresponding finite element model consists of three beam elements and twelve degrees of freedom.

Example:

The computation is initialized by importing CALFEM and NumPy. The element topology matrix contains only the degrees of freedom for each element:

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu
import calfem.vis_mpl as cfv

edof = np.array([
    [4, 5, 6, 1, 2, 3],
    [7, 8, 9, 10, 11, 12],
    [4, 5, 6, 7, 8, 9]
])

The global stiffness matrix K and load vector f are initialized. The point load P = 2000 N is applied at DOF 4:

CodeOutput

K = np.array(np.zeros((12, 12)))
f = np.array(np.zeros((12, 1)))
f[3] = 2.0e3

The material and geometric properties are defined for each element type:

E = 200.0e9
A1 = 2.0e-3
A2 = 6.0e-3
I1 = 1.6e-5
I2 = 5.4e-5

ep1 = [E, A1, I1]
ep3 = [E, A2, I2]

ex1 = np.array([0, 0])
ex2 = np.array([6, 6])
ex3 = np.array([0, 6])

ey1 = np.array([4, 0])
ey2 = np.array([4, 0])
ey3 = np.array([4, 4])

eq1 = np.array([0, 0])
eq2 = np.array([0, 0])
eq3 = np.array([0, -10e3])

The element stiffness matrices and load vectors are computed and assembled:

Ke1 = cfc.beam2e(ex1, ey1, ep1)
Ke2 = cfc.beam2e(ex2, ey2, ep1)
Ke3, fe3 = cfc.beam2e(ex3, ey3, ep3, eq3)

K = cfc.assem(edof[0, :], K, Ke1)
K = cfc.assem(edof[1, :], K, Ke2)
K, f = cfc.assem(edof[2, :], K, Ke3, f, fe3)

The system of equations is solved by specifying boundary conditions and using solveq():

CodeOutput

bc_dofs = np.array([1, 2, 3, 10, 11])
bc_vals = np.array([0.0, 0.0, 0.0, 0.0, 0.0])

a, r = cfc.solveq(K, f, bc_dofs, bc_vals)

cfu.disp_array(a, ["a"])
cfu.disp_array(r, ["r"])

The element displacements are extracted and section forces are computed along each element:

CodeOutput

ed = cfc.extract_ed(edof, a)

es1, edi1, ec1 = cfc.beam2s(ex1, ey1, ep1, ed[0, :], eq1, nep=21)
es2, edi2, ec2 = cfc.beam2s(ex2, ey2, ep1, ed[1, :], eq2, nep=21)
es3, edi3, ec3 = cfc.beam2s(ex3, ey3, ep3, ed[2, :], eq3, nep=21)

cfu.disp_h2("es1")
cfu.disp_array(es1, ["N", "Vy", "Mz"])
cfu.disp_h2("edi1")
cfu.disp_array(edi1, ["u1", "v1"])
cfu.disp_h2("es2")
cfu.disp_array(es2, ["N", "Vy", "Mz"])
cfu.disp_h2("edi2")
cfu.disp_array(edi2, ["u1", "v1"])
cfu.disp_h2("es3")
cfu.disp_array(es3, ["N", "Vy", "Mz"])
cfu.disp_h2("edi3")
cfu.disp_array(edi3, ["u1", "v1"])

Visualization of the frame analysis results using CALFEM’s visualization functions:

plotpar = [2, 1, 0]
sfac = cfv.scalfact2(ex3, ey3, edi3, 0.1)
print("sfac=")
print(sfac)

cfv.figure(1)
cfv.eldraw2(ex1, ey1, plotpar)
cfv.eldraw2(ex2, ey2, plotpar)
cfv.eldraw2(ex3, ey3, plotpar)

plotpar = [1, 2, 1]
cfv.dispbeam2(ex1, ey1, edi1, plotpar, sfac)
cfv.dispbeam2(ex2, ey2, edi2, plotpar, sfac)
cfv.dispbeam2(ex3, ey3, edi3, plotpar, sfac)
cfv.axis([-1.5, 7.5, -0.5, 5.5])
plotpar1 = 2
cfv.scalgraph2(sfac, [1e-2, 0.5, 0], plotpar1)
cfv.title("Displacements")

# ----- Draw normal force diagram --------------------------------

plotpar = [2, 1]
sfac = cfv.scalfact2(ex1, ey1, es1[:, 0], 0.2)
cfv.figure(2)
cfv.secforce2(ex1, ey1, es1[:, 0], plotpar, sfac)
cfv.secforce2(ex2, ey2, es2[:, 0], plotpar, sfac)
cfv.secforce2(ex3, ey3, es3[:, 0], plotpar, sfac)
cfv.axis([-1.5, 7.5, -0.5, 5.5])
cfv.scalgraph2(sfac, [3e4, 1.5, 0], plotpar1)
cfv.title("Normal force")

# ----- Draw shear force diagram ---------------------------------

plotpar = [2, 1]
sfac = cfv.scalfact2(ex3, ey3, es3[:, 1], 0.2)
cfv.figure(3)
cfv.secforce2(ex1, ey1, es1[:, 1], plotpar, sfac)
cfv.secforce2(ex2, ey2, es2[:, 1], plotpar, sfac)
cfv.secforce2(ex3, ey3, es3[:, 1], plotpar, sfac)
cfv.axis([-1.5, 7.5, -0.5, 5.5])
cfv.scalgraph2(sfac, [3e4, 0.5, 0], plotpar1)
cfv.title("Shear force")

# ----- Draw moment diagram --------------------------------------

plotpar = [2, 1]
sfac = cfv.scalfact2(ex3, ey3, es3[:, 2], 0.2)
print("sfac=")
print(sfac)

cfv.figure(4)
cfv.secforce2(ex1, ey1, es1[:, 2], plotpar, sfac)
cfv.secforce2(ex2, ey2, es2[:, 2], plotpar, sfac)
cfv.secforce2(ex3, ey3, es3[:, 2], plotpar, sfac)
cfv.axis([-1.5, 7.5, -0.5, 5.5])
cfv.scalgraph2(sfac, [3e4, 0.5, 0], plotpar1)
cfv.title("Moment")

cfv.show_and_wait()

exs_beambar2¶

Purpose:

Analysis of a combined beam and bar structure.

Description:

Consider a structure consisting of a beam with \(A_1=4.0 \times 10^{-3}\) m² and \(I_1=5.4 \times 10^{-5}\) m⁴ supported by two bars with \(A_2=1.0 \times 10^{-3}\) m². The beam as well as the bars have \(E=200\) GPa. The structure is loaded by a distributed load \(q=10\) kN/m. The corresponding finite element model consists of three beam elements and two bar elements and has 14 degrees of freedom.

Example:

The computation is initialized by importing CALFEM and NumPy. The topology matrices are defined separately for beam and bar elements, containing only the degrees of freedom:

import numpy as np
import calfem.core as cfc
import calfem.utils as cfu
import calfem.vis_mpl as cfv
edof1 = np.array([
    [1,  2,  3,  4,  5,  6],
    [4,  5,  6,  7,  8,  9],
    [7,  8,  9, 10, 11, 12]
])

edof2 = np.array([
    [13, 14,  4,  5],
    [13, 14,  7,  8]
])

K = np.array(np.zeros((14, 14)))
f = np.array(np.zeros((14, 1)))

The material and geometric properties are defined for beam and bar elements:

E = 200.e9
A1 = 4.e-3
I1 = 5.4e-5
A2 = 1.e-3

ep1 = [E, A1, I1]
ep4 = [E, A2]

eq1 = [0, 0]
eq2 = [0, -10e+3]

ex1 = np.array([0, 2])
ex2 = np.array([2, 4])
ex3 = np.array([4, 6])
ex4 = np.array([0, 2])
ex5 = np.array([0, 4])
ey1 = np.array([2, 2])
ey2 = np.array([2, 2])
ey3 = np.array([2, 2])
ey4 = np.array([0, 2])
ey5 = np.array([0, 2])

The element stiffness matrices are computed using beam2e() for beam elements and bar2e() for bar elements. Element load vectors from distributed loads are also computed:

Ke1 = cfc.beam2e(ex1, ey1, ep1)
Ke2, fe2 = cfc.beam2e(ex2, ey2, ep1, eq2)
Ke3, fe3 = cfc.beam2e(ex3, ey3, ep1, eq2)

Ke4 = cfc.bar2e(ex4, ey4, ep4)
Ke5 = cfc.bar2e(ex5, ey5, ep4)

The global stiffness matrix and load vector are assembled using assem():

K = cfc.assem(edof1[0, :], K, Ke1)
K, f = cfc.assem(edof1[1, :], K, Ke2, f, fe2)
K, f = cfc.assem(edof1[2, :], K, Ke3, f, fe3)
K = cfc.assem(edof2[0, :], K, Ke4)
K = cfc.assem(edof2[1, :], K, Ke5)

The system of equations is solved by specifying boundary conditions and using solveq(). The vertical displacement at the end of the beam is 13.0 mm:

CodeOutput

bc_dofs = np.array([1, 2, 3, 13, 14])
bc_vals = np.array([0.0, 0.0, 0.0, 0.0, 0.0])

a, r = cfc.solveq(K, f, bc_dofs, bc_vals)
cfu.disp_h2("Displacements a:")
cfu.disp_array(a)

cfu.disp_h2("Reaction forces r:")
cfu.disp_array(r)

Maximum vertical displacement: 0.009 m = 9.5 mm

The section forces are calculated using beam2s() and bar2s() from element displacements. This yields normal forces of -35.4 kN and -152.5 kN in the bars and maximum moment of 20.0 kNm in the beam:

CodeOutput

ed1 = cfc.extract_ed(edof1, a)
ed2 = cfc.extract_ed(edof2, a)

es1, _, _ = cfc.beam2s(ex1, ey1, ep1, ed1[0, :], eq1, nep=11)
es2, _, _ = cfc.beam2s(ex2, ey2, ep1, ed1[1, :], eq2, nep=11)
es3, _, _ = cfc.beam2s(ex3, ey3, ep1, ed1[2, :], eq2, nep=11)
es4 = cfc.bar2s(ex4, ey4, ep4, ed2[0, :])
es5 = cfc.bar2s(ex5, ey5, ep4, ed2[1, :])

cfu.disp_h2("es1 = ")
cfu.disp_array(es1, headers=["N", "Q", "M"])
cfu.disp_h2("es2 = ")
cfu.disp_array(es2, headers=["N", "Q", "M"])
cfu.disp_h2("es3 = ")
cfu.disp_array(es3, headers=["N", "Q", "M"])
cfu.disp_h2("es4 = ")
cfu.disp_array(es4, headers=["N"])
cfu.disp_h2("es5 = ")
cfu.disp_array(es5, headers=["N"])