Example: Bars

This example is from the CALFEM manual (exs2.py).

Purpose:

Analysis the displacement at point of loading on a plane truss

Description:

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

bars1

[1]:
import numpy as np
import calfem.core as cfc

Topology matrix Edof

Define element topology. Element number are not used in the Edof matrix in CALFEM for Python.

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

Stiffness matrix K and load vector f

[3]:
K = np.matrix(np.zeros((8,8)))
f = np.matrix(np.zeros((8,1)))

Element properties

[4]:
E = 2.0e11
A1 = 6.0e-4
A2 = 3.0e-4
A3 = 10.0e-4
ep1 = [E,A1]
ep2 = [E,A2]
ep3 = [E,A3]

Element coordinates

[5]:
ex1 = np.array([0., 1.6])
ex2 = np.array([1.6, 1.6])
ex3 = np.array([0., 1.6])

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

Element stiffness matrices

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

Assemble Ke into K

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

[7]:
cfc.assem(Edof[0,:],K,Ke1)
cfc.assem(Edof[1,:],K,Ke2)
cfc.assem(Edof[2,:],K,Ke3)
[7]:
matrix([[ 7.50e+07,  0.00e+00,  0.00e+00,  0.00e+00, -7.50e+07,
          0.00e+00,  0.00e+00,  0.00e+00],
        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
          0.00e+00,  0.00e+00,  0.00e+00],
        [ 0.00e+00,  0.00e+00,  6.40e+07, -4.80e+07, -6.40e+07,
          4.80e+07,  0.00e+00,  0.00e+00],
        [ 0.00e+00,  0.00e+00, -4.80e+07,  3.60e+07,  4.80e+07,
         -3.60e+07,  0.00e+00,  0.00e+00],
        [-7.50e+07,  0.00e+00, -6.40e+07,  4.80e+07,  1.39e+08,
         -4.80e+07,  0.00e+00,  0.00e+00],
        [ 0.00e+00,  0.00e+00,  4.80e+07, -3.60e+07, -4.80e+07,
          8.60e+07,  0.00e+00, -5.00e+07],
        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
          0.00e+00,  0.00e+00,  0.00e+00],
        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         -5.00e+07,  0.00e+00,  5.00e+07]])
[8]:
print("Stiffness matrix K:")
print(K)
Stiffness matrix K:
[[ 7.50e+07  0.00e+00  0.00e+00  0.00e+00 -7.50e+07  0.00e+00  0.00e+00
   0.00e+00]
 [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00
   0.00e+00]
 [ 0.00e+00  0.00e+00  6.40e+07 -4.80e+07 -6.40e+07  4.80e+07  0.00e+00
   0.00e+00]
 [ 0.00e+00  0.00e+00 -4.80e+07  3.60e+07  4.80e+07 -3.60e+07  0.00e+00
   0.00e+00]
 [-7.50e+07  0.00e+00 -6.40e+07  4.80e+07  1.39e+08 -4.80e+07  0.00e+00
   0.00e+00]
 [ 0.00e+00  0.00e+00  4.80e+07 -3.60e+07 -4.80e+07  8.60e+07  0.00e+00
  -5.00e+07]
 [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00
   0.00e+00]
 [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00 -5.00e+07  0.00e+00
   5.00e+07]]

Solve the system of equations

Considering the prescribed displacements in bc, the system of equations is solved using the function solveq, yielding displacements a and support forces r.

[9]:
bc = np.array([1,2,3,4,7,8])
f[5] = -80e3
a, r = cfc.solveq(K,f,bc)

print("Displacements a:")
print(a)

print("Reaction forces r:")
print(r)
Displacements a:
[[ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [-0.00039793]
 [-0.00115233]
 [ 0.        ]
 [ 0.        ]]
Reaction forces r:
[[ 29844.55958549]
 [     0.        ]
 [-29844.55958549]
 [ 22383.41968912]
 [     0.        ]
 [     0.        ]
 [     0.        ]
 [ 57616.58031088]]

Element forces

The vertical displacement at the point of loading is 1.15 mm. The negative values shows the direction is to the bottom, against the direction of node 6. The section forces es1, es2 and es3 are calculated using bar2s from element displacements ed1, ed2 and ed3 obtained using extract

[10]:
ed1 = cfc.extract_eldisp(Edof[0,:],a);
N1 = cfc.bar2s(ex1,ey1,ep1,ed1)
ed2 = cfc.extract_eldisp(Edof[1,:],a);
N2 = cfc.bar2s(ex2,ey2,ep2,ed2)
ed3 = cfc.extract_eldisp(Edof[2,:],a);
N3 = cfc.bar2s(ex3,ey3,ep3,ed3)
print("N1=",N1)
print("N2=",N2)
print("N3=",N3)
N1= -29844.559585492225
N2= 57616.580310880825
N3= 37305.69948186528