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: