Theoretical Basis#

PyCBA is based on the matrix stiffness method.

Note#

Until this section is completed, please also see the Matrix Stiffness Method lecture notes here for more information.

Elements#

The stiffness matrices for the elements used in PyCBA are given below.

Note that the releases relate only to internal hinges within the elements themselves, and not to the supports conditions. Note also, that in modelling joints, there must always be one element with a ‘fixed’ end at a joint; otherwise the resulting global stiffness matrix will be singular.

Type 1: Fixed-Fixed#

\[\begin{split}\bf{K} = \frac{EI}{L^3} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2 \\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}\end{split}\]

Type 2: Fixed-Pinned#

\[\begin{split}\bf{K} = \frac{3EI}{L^3} \begin{bmatrix} 1 & L & -1 & 0 \\ L & L^2 & -L & 0 \\ -1 & -L & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}\end{split}\]

Type 3: Pinned-Fixed#

\[\begin{split}\bf{K} = \frac{3EI}{L^3} \begin{bmatrix} 1 & 0 & -1 & L \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & -L \\ L & 0 & -L & L^2 \end{bmatrix}\end{split}\]

Type 4: Pinned-Pinned#

\[\begin{split}\bf{K} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}\end{split}\]

Displacements#

Displacement along an element are found using cumulative trapezoidal integration to approximate the integral. It is not accurate to use the beam element shape functions due to the presence of inter-nodal loading. Thus, the rotations along an element are found by integrating the curvatures:

\[\theta(x) = \theta_0 + \int_{x=0}^{x=L} \frac{M(x)}{EI} dx,\]

and similarly the displacements then found by integrating the rotations:

\[\delta(x) = \delta_0 + \int_{x=0}^{x=L} \theta(x) dx.\]

In both cases we approximate the integral with:

\[\int_0^L y(x) dx \approx \sum_{i=1}^{i=n} \frac{y_{i-1}+y_i}{2} \Delta x.\]

which is implemented using SciPy’s cumulative trapezoidal integration along the \(n\) integration points along the member.

Note that \(\theta_0\) and \(\delta_0\) are required: the rotation and translation at the \(i\)-node of the element. This is typically the output of the global system of equations; i.e., the nodal displacements. However, for elements with releases (types 2, 3, 4 above), we need to determine the starting rotation, \(\theta_0\). From the slope deflection equations we can find this starting rotation at the \(i\)-node:

\[\theta_i = \frac{\delta_j - \delta_i}{L} - \frac{L}{3EI}\left(-FEM_i + 0.5FEM_j + M_i - 0.5 M_j\right)\]

in which the \(FEM\) are the consistent nodal loading moments due to the inter-nodal loads, and the \(M\) moments are solved from the global stiffness equation.

Global Stiffness Matrix#

For an \(N\)-span beam there are \(N+1\) nodes. Each node has two degrees of freedom (DOFs): a vertical displacement \(v\) and a rotation \(\theta\). The global DOF vector is of length \(2(N+1)\), ordered node by node:

\[\mathbf{d} = \begin{bmatrix} v_0 & \theta_0 & v_1 & \theta_1 & \cdots & v_N & \theta_N \end{bmatrix}^T\]

The unrestricted global stiffness matrix \(\mathbf{K}\) of size \(2(N+1) \times 2(N+1)\) is assembled by the direct stiffness method. For span \(m\) (zero-indexed), the left node is global node \(m\) and the right node is global node \(m+1\), mapping to global DOF indices \([2m,\; 2m+1,\; 2m+2,\; 2m+3]\). The \(4\times 4\) element stiffness matrix \(\mathbf{k}^{(e)}\) (one of the four types above) is overlapped into \(\mathbf{K}\):

\[K_{[2m:2m+4],\,[2m:2m+4]} \mathrel{+}= \mathbf{k}^{(e)}\]

The overlap at interior nodes automatically enforces displacement compatibility between adjacent spans. The result is the banded, symmetric, positive semi-definite matrix characteristic of beam structures.

Spring Supports#

A spring support at DOF \(i\) with stiffness \(k_s > 0\) contributes a restoring force \(k_s u_i\) at that DOF. This is incorporated directly into the global stiffness matrix before applying boundary conditions:

\[K_{ii} \leftarrow K_{ii} + k_s\]

A spring DOF remains as a free unknown in the linear system — no row or column elimination is performed. The spring force is recovered after solving as:

\[F_s^{(i)} = k_s \, u_i\]

Boundary Conditions#

The governing system of equations is:

\[\mathbf{K}\,\mathbf{d} = \mathbf{f}\]

where \(\mathbf{f}\) is the equivalent nodal force vector assembled from the consistent nodal loads of each span (negated, as loads oppose the reactions they generate).

Boundary conditions are applied by the direct elimination method. A DOF \(i\) is eliminated whenever its displacement is known: either because an explicit value \(\bar{d}_i\) has been prescribed in the displacement vector \(\mathbf{D}\), or because the DOF is fully fixed (\(R_i = -1\)) and no override is given, in which case \(\bar{d}_i = 0\). For each such DOF the procedure is:

Transfer the constraint’s contribution to all other DOFs:

\[f_j \leftarrow f_j - K_{ji}\,\bar{d}_i \quad \forall\, j\]
Zero the \(i\)-th row and column:

\[K_{ij} = K_{ji} = 0 \quad \forall\, j\]
Set the diagonal to unity and the right-hand side to the prescribed value:

\[K_{ii} = 1, \qquad f_i = \bar{d}_i\]

This preserves symmetry and reproduces the prescribed displacement exactly in the solution. A support settlement is therefore modelled simply by providing a non-zero \(\bar{d}_i\) at a fixed DOF; a negative value corresponds to a downward settlement. Prescribed displacements may also be applied to free DOFs. Spring DOFs (\(R_i > 0\)) without a prescribed displacement are not eliminated; their stiffness is already on the diagonal and they remain as free unknowns.

Reaction Recovery#

After solving the restricted system for the nodal displacements \(\mathbf{d}\), support reactions are recovered using the unrestricted stiffness matrix \(\mathbf{K}_U\) (assembled before boundary conditions were applied) and the original nodal force vector \(\mathbf{f}_U\):

\[\mathbf{r}^* = \mathbf{K}_U\,\mathbf{d} - \mathbf{f}_U\]

The reaction at a fully-fixed DOF \(i\) (\(R_i = -1\)) is:

\[R_i = r^*_i\]

For spring DOFs the spring force is reported separately as \(F_s^{(i)} = -k_s\,u_i\) (upward positive), because the residual \(r^*_i\) at a spring DOF contains structural coupling contributions from neighbouring DOFs in addition to the spring force itself.