The above figure shows the numbering and sign convention used for the six end displacements and forces for a general 2D frame member. These are shown relative to the local coordinate system of the member, with origin at the left, or j-end, and x-axis coincident with the axis of the member. The y-axis is positive upwards. It is sometimes necessary to distinguish member ends; the term j-end is used to refer to the left end or origin, and k-end is the other end (note that the CAL manual uses the i- and j- end for the same purpose).
Of course $f_0$ and $f_3$ are axial forces, $f_1$ and $f_4$ are shear forces, and $f_2$ and $f_5$ are bending moments.
$\newcommand{\mat}[1]{\left[\begin{matrix}#1\end{matrix}\right]}$
The local member stiffness matrix, $\mat{K_l}$, expresses the end forces as a function of the end displacements, thus: $$ \mat{f_0\\ f_1\\ f_2\\ f_3\\ f_4\\ f_5} = \mat{K_l} \mat{\delta_0\\ \delta_1\\ \delta_2\\ \delta_3\\ \delta_4\\ \delta_5} $$
In [1]:
from sympy import *
init_printing(use_latex='mathjax')
from sympy.matrices import Matrix
We will first determine the inverse relationship for end moments -- the rotation flexibility matrix, $\mat{F}$, that expresses the end rotations as a function of the end moments, as shown in the following figure.
$$ \mat{\delta_2\\ \delta_5} = \mat{F} \mat{f_2\\f_5} $$Note that at each end, there are two other forces acting - axial and shear: $f_0$, $f_1$ and $f_3$, $f_4$. These are not shown, for simplicity of figure.
The method of virtual work is probably the simplest and most direct way to do develop the rotation flexibility matrix.
The following figure shows unit moments placed individually at each end (at location #s 2 and 5). Below each segment is the corresponnding bending moment diagram, $m_i$:
The figure also shows the resulting rotational displacements $\alpha_{ij}$ for each segment; these are the coefficients of the flexibility matrix:
$$ \mat{F} = \mat{\alpha_{22} & \alpha_{52}\\ \alpha_{25} & \alpha_{55}} $$The coefficients of $\mat{F}$, $\alpha_{ij}$, are the end displacements due to unit values of the end forces. Specifically $\alpha_{ij}$ is the rotation at $i$ due to a unit moment at $j$ and is calculated using the method of virtual work (integrating products of bending moments):
$$ \mat{F} = \mat{\int_0^L\frac{m_2 m_2}{EI} dx & \int_0^L\frac{m_2 m_5}{EI} dx \\ \int_0^L\frac{m_5 m_2}{EI} dx & \int_0^L\frac{m_5 m_5}{EI} dx } $$
In [2]:
L,E,I,A,x = symbols('L E I A x')
m2 = x/L - 1 # unit moment at DOF 2
m5 = x/L # unit moment at DOF 5
EI = E*I
F = Matrix([[m2*m2, m2*m5],
[m5*m2, m5*m5]]).integrate((x,0,L))/EI
F
Out[2]:
In [3]:
M = F**-1
M
Out[3]:
In [4]:
# transform end moments to end thrust,shear,moment
Tf = Matrix([[0,0],[1/L,1/L],[1,0],[0,0],[-1/L,-1/L],[0,1]])
Tf
Out[4]:
In [5]:
f2,f5 = symbols('f2 f5') # confirm that that is OK
Tf*Matrix([[f2],
[f5]])
Out[5]:
Column $j$ of the member stiffness matrix $\mat{K_l}$ consists of all six end forces that are consistent with a unit value of displacement $j$ and zero values of the other five displacements.
Therefore column 2 consists of the six end forces that occur when $\delta_2 = 1$ and all other displacements are 0. That can be computed:
$$ \mat{T_f} \mat{M} \mat{1\\0} $$
In [6]:
Kl = zeros(6,6) # build the member stiffness matrix one column at a time
Kl[:,2] = Tf * M * Matrix([[1],[0]]) # column 2
Kl[:,2]
Out[6]:
Column 5 is similar, but with $\delta_5 = 1$.
In [7]:
Kl[:,5] = Tf * M * Matrix([[0],[1]]) # column 5
Kl[:,5]
Out[7]:
The forces due to end translations are determined by mapping the unit translation to end rotations of the elastic curve from the chord.
That is, for the case of $\delta_1 = 1$, (and notabley $\delta_2 = \delta_5 = 0$), then the end moments will be consistent with end rotations of $1/L$. Column 1 will be calculated:
$$ \mat{T_f} \mat{M} \mat{1/L\\1/L} $$
In [8]:
Kl[:,1] = Tf * M * Matrix([[1/L],[1/L]]) # column 1
Kl[:,1]
Out[8]:
The end forces due to a unit value of $\delta_4$ are the same as those for $\delta_1$, but with reversed signs.
In [9]:
Kl[:,4] = Tf * M * Matrix([[-1/L],[-1/L]]) # column 4
Kl[:,4]
Out[9]:
Now the axial forces and displacements, which are de-coupled from shears and moments:
In [10]:
ac = Matrix([E*A/L,0,0,-E*A/L,0,0])
ac
Out[10]:
In [11]:
Kl[:,0] = ac # column 0
Kl[:,3] = -ac # column 3
In [12]:
Kl
Out[12]:
And, for ease of copying into other software, here it is again.
In [13]:
print(Kl)
In [ ]: