ch2_theory.md

Chapter 2: Theory

This chapter presents the element formulation, global assembly, solution procedures, stress recovery, and the optimisation algorithm implemented in ShepherdLab LW FEA-OPT.

Element Formulation

Degrees of Freedom

Each node carries three degrees of freedom (DOFs):

$$ \mathbf{u}_i = \begin{bmatrix} u_{x,i} \ u_{y,i} \ \theta_i \end{bmatrix} $$

where u_x and u_y are translational displacements and θ is the in-plane rotation. A single element connecting nodes i and j therefore has a 6×1 displacement vector:

$$ \mathbf{u}^{(e)} = \begin{bmatrix} u_{x,i} & u_{y,i} & \theta_i & u_{x,j} & u_{y,j} & \theta_j \end{bmatrix}^\mathsf{T} $$

Local Stiffness Matrix

In the local coordinate system aligned with the element axis, the standard 2-D Euler–Bernoulli frame stiffness matrix is:

$$ \mathbf{k}^{(e)} = \begin{bmatrix} \frac{EA}{L} & 0 & 0 & -\frac{EA}{L} & 0 & 0 \\ 0 & \frac{12EI}{L^3} & \frac{6EI}{L^2} & 0 & -\frac{12EI}{L^3} & \frac{6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{4EI}{L} & 0 & -\frac{6EI}{L^2} & \frac{2EI}{L} \\ -\frac{EA}{L} & 0 & 0 & \frac{EA}{L} & 0 & 0 \\ 0 & -\frac{12EI}{L^3} & -\frac{6EI}{L^2} & 0 & \frac{12EI}{L^3} & -\frac{6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{2EI}{L} & 0 & -\frac{6EI}{L^2} & \frac{4EI}{L} \end{bmatrix} $$

where E is the elastic modulus, A = w·t is the cross-sectional area, I = w·t³/12 is the second moment of area for a rectangular section of width w and thickness t, and L is the element length.

All local stiffness matrices are built simultaneously as a (M, 6, 6) NumPy array using vectorised operations, avoiding per-element Python loops.

Coordinate Transformation

An element oriented at angle α (measured counter-clockwise from the global x-axis) is rotated into the global frame by:

$$ \mathbf{T}^{(e)} = \begin{bmatrix} \mathbf{R} & \mathbf{0} \\ \mathbf{0} & \mathbf{R} \end{bmatrix}, \qquad \mathbf{R} = \begin{bmatrix} \cos\alpha & \sin\alpha & 0 \\ -\sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

The global element stiffness matrix is then:

$$ \mathbf{K}^{(e)} = {\mathbf{T}^{(e)}}^\mathsf{T}, \mathbf{k}^{(e)}, \mathbf{T}^{(e)} $$

The batched matrix product for all elements is computed as a single (M, 6, 6) operation using NumPy's @ operator on 3-D arrays.

Rigid Elements

Rigid elements (e.g. attachment fixtures) are modelled using a stiffness penalty method: the elastic modulus is multiplied by a large factor β (default 10²), giving E_rigid = β · E. This is simpler than Lagrange multiplier approaches and integrates directly with the standard assembly procedure.

Cross-Section Properties

All elements use a rectangular cross-section:

$$ A = w,t, \qquad I = \frac{w,t^3}{12} $$

where w is the (constant) out-of-plane width and t is the thickness. During optimisation, only t varies; w remains fixed.

---

Global Assembly and Solution

Assembly

The global stiffness matrix is assembled from element contributions by the standard direct stiffness method:

$$ \mathbf{K} = \sum_{e=1}^{M} \mathbf{L}^{(e)\mathsf{T}}, \mathbf{K}^{(e)}, \mathbf{L}^{(e)} $$

where L⁽ᵉ⁾ is the Boolean localisation (scatter) matrix that maps the six element DOFs to the global DOF vector. In practice, the assembly is implemented by building vectorised row, column, and value arrays from all elements simultaneously and constructing a scipy.sparse.csr_matrix in a single call.

Partitioning

The global system K u = P is partitioned into free (subscript f) and fixed (subscript c) DOFs:

$$ \mathbf{K}_{ff},\mathbf{u}_f + \mathbf{K}_{fc},\mathbf{u}_c = \mathbf{P}_f $$

$$ \mathbf{K}_{cf},\mathbf{u}_f + \mathbf{K}_{cc},\mathbf{u}_c = \mathbf{R}_c $$

For prescribed zero displacements ($\mathbf{u}_c = \mathbf{0}$), the reduced system is:

$$ \mathbf{K}_{ff}, \mathbf{u}_f = \mathbf{P}_f $$

which is solved by scipy.sparse.linalg.spsolve.

---

Newton–Raphson Large-Deflection Analysis

For problems with large displacements, the equilibrium equation becomes nonlinear because the element geometry (lengths and orientations) changes with deformation. LW FEA-OPT uses an incremental load-stepping scheme with Newton–Raphson iterations at each step.

Load Stepping

The total load P is divided into N equal increments:

$$ \mathbf{P}_n = \frac{n}{N}, \mathbf{P}_{\text{ref}}, \qquad n = 1, 2, \ldots, N $$

Equilibrium and Residual

At load step n, the goal is to find displacements u such that the internal forces balance the external forces:

$$ \mathbf{R}(\mathbf{u}) = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}}(\mathbf{u}) = \mathbf{0} $$

The internal force vector is assembled element-by-element:

$$ \mathbf{f}_{\text{int}} = \sum_{e=1}^{M} {\mathbf{T}^{(e)}}^\mathsf{T}, \mathbf{k}^{(e)}, \mathbf{T}^{(e)}, \mathbf{u}^{(e)} $$

where T⁽ᵉ⁾ and k⁽ᵉ⁾ are evaluated at the current deformed configuration. The assembly uses np.einsum for batched matrix–vector products and np.add.at for scatter-addition into the global vector.

Iteration

At each Newton–Raphson iteration k:

1. Update element geometry (lengths, orientations) from the current displacement field.
2. Assemble the tangent stiffness matrix K_t = Σ_e T⁽ᵉ⁾ᵀ k⁽ᵉ⁾ T⁽ᵉ⁾.
3. Compute the residual R⁽ᵏ⁾ = f_ext − f_int⁽ᵏ⁾.
4. Solve for the displacement increment: K_{t,ff} Δu_f = R_f⁽ᵏ⁾.
5. Update: u_f⁽ᵏ⁺¹⁾ = u_f⁽ᵏ⁾ + Δu_f.

Convergence Criteria

Three criteria are checked simultaneously:

$$ \text{Force:} \quad \frac{\lVert\mathbf{R}\rVert}{\lVert\mathbf{f}_{\text{ext}}\rVert} < \varepsilon_f $$

$$ \text{Displacement:} \quad \frac{\lVert\Delta\mathbf{u}\rVert}{u_{\max}} < \varepsilon_u $$

$$ \text{Energy:} \quad |\mathbf{R}^\mathsf{T} \Delta\mathbf{u}| < \varepsilon_E $$

Convergence is declared when the force criterion is satisfied and either the displacement or energy criterion is met. On the first iteration, only the force criterion is used (the displacement increment has no meaningful predecessor).

The default tolerances are ε_f = 10⁻², ε_u = 10⁻⁶, and ε_E = 10⁻⁶.

---

Stress Recovery

After solving for displacements, member end forces are recovered in local coordinates:

$$ \mathbf{f}_{\text{local}}^{(e)} = \mathbf{k}^{(e)}, \mathbf{T}^{(e)}, \mathbf{u}^{(e)} $$

giving force components f_local = [N_i, V_i, M_i, N_j, V_j, M_j]ᵀ at the start (i) and end (j) nodes. This is computed for all elements simultaneously using np.einsum.

Stresses at each node of an element are:

$$ \sigma_{\text{axial}} = \frac{N}{A} $$

$$ \sigma_{\text{bending}} = \frac{M \cdot c}{I}, \quad c = t/2 $$

$$ \tau_{\text{shear}} = \frac{3V}{2A} $$

where the shear stress uses the rectangular-section maximum (at the neutral axis). The combined von Mises stress is:

$$ \sigma_{\text{vM}} = \sqrt{(\sigma_{\text{axial}} + \sigma_{\text{bending}})^2 + 3,\tau_{\text{shear}}^2} $$

The factor of safety (FOS) for element e is:

$$ \text{FOS}^{(e)} = \frac{\sigma_{\text{ult}}}{\max!\big(\sigma_{\text{vM},i}^{(e)},; \sigma_{\text{vM},j}^{(e)}\big)} $$

---

Stress-Ratio Thickness Optimisation

The optimiser implements a fully stressed design (FSD) approach, which iteratively adjusts element thicknesses so that every element operates near its allowable stress. The result is a more uniform stress distribution and, typically, a lighter structure.

Update Rule

At iteration k, the thickness of element e is updated as:

$$ t^{(k+1)}_e = t^{(k)}_e \left( \frac{\sigma^{(k)}_e}{\sigma_{\text{allow}}} \right)^{!\alpha} $$

where σ⁽ᵏ⁾_e is the maximum von Mises stress in element e at the current iteration, σ_allow = σ_ult / FOS_target, and α ∈ (0, 1] is a damping exponent that controls aggressiveness.

Choosing α: Smaller values of α produce more cautious updates and are less likely to oscillate, at the cost of slower convergence. A typical range is 0.3–0.5. For stiff structures, α = 0.4 works well; for flexible trusses, α = 0.3 may be more stable.

Move Limits

To prevent oscillation, the relative change per iteration is clamped:

$$ t^{(k+1)}_e \in \Big[, t^{(k)}_e (1 - \delta),;; t^{(k)}_e (1 + \delta) ,\Big] $$

where δ is the move limit (default 0.2, i.e. ±20% per iteration).

Absolute Bounds

After applying the move limit, thickness is clamped to user-specified bounds:

$$ t^{(k+1)}_e \in [, t_{\min},; t_{\max} ,] $$

Setting a small t_min (e.g. 0.2 mm) allows elements to effectively vanish, enabling topology-like optimisation where load paths emerge organically.

Convergence

Convergence is declared when the maximum relative thickness change across all optimisable elements falls below a tolerance ε_t:

$$ \max_e \frac{|t^{(k+1)}_e - t^{(k)}_e|}{t^{(k)}_e} < \varepsilon_t $$

Stress Uniformity

The quality of the optimised design is measured by the coefficient of variation (CoV) of the maximum element stresses:

$$ \text{CoV} = \frac{\text{std}(\sigma_e)}{\text{mean}(\sigma_e)} $$

A lower CoV indicates more uniform material utilisation.

---

Truss Subdivision

To enable topology-like optimisation of truss structures, each original member can be subdivided into n smaller frame elements using subdivide_truss. This creates n−1 intermediate nodes per member, uniformly spaced along the original member axis.

The fine mesh then allows the thickness optimiser to vary material along each member, not just between members. Combined with a small t_min, the optimiser can effectively remove material from low-stress regions, producing organic, topology-optimised shapes.

The member_map output tracks which original member each sub-element belongs to, which is useful for grouping, colouring, and interpreting results.