Theory & Formulation

This page covers the finite element theory implemented in Cauchy: the weak form of 2D elasticity, isoparametric shape functions, Gauss quadrature, element stiffness assembly, and solver algorithms.

Weak Form of 2D Elasticity

The strong form seeks displacement field u(x) satisfying equilibrium div(sigma) + b = 0 with boundary conditions. The weak (variational) form multiplies by a test function v and integrates by parts:

$$ \int_\Omega \sigma : \varepsilon(v) \, d\Omega = \int_\Omega b \cdot v \, d\Omega + \int_{\Gamma_t} t \cdot v \, d\Gamma $$

where sigma is the stress tensor, epsilon is the strain tensor, b is body force, and t is traction on the Neumann boundary.

Constitutive Law (D-Matrix)

For linear isotropic elasticity, the stress-strain relationship sigma = D * epsilon is:

Plane Stress (thin plates)

$$ D = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} $$

Plane Strain (thick sections)

$$ D = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix} $$

Q4 Shape Functions

The 4-node bilinear quad uses natural coordinates (xi, eta) in [-1, +1]:

$$ N_i(\xi, \eta) = \frac{1}{4}(1 + \xi_i \xi)(1 + \eta_i \eta) $$

Node ordering: 0=(-1,-1), 1=(+1,-1), 2=(+1,+1), 3=(-1,+1) (CCW). Evaluated using 2x2 Gauss quadrature.

Q8 Shape Functions (Serendipity)

The 8-node serendipity element uses quadratic shape functions with 4 corner nodes and 4 midside nodes. This eliminates the internal node of a Q9 Lagrange element while retaining quadratic completeness:

Corners: $$ N_i = \frac{1}{4}(1 + \xi_i \xi)(1 + \eta_i \eta)(\xi_i \xi + \eta_i \eta - 1) $$

Horizontal midside (xi=0): $$ N_i = \frac{1}{2}(1 - \xi^2)(1 + \eta_i \eta) $$

Vertical midside (eta=0): $$ N_i = \frac{1}{2}(1 + \xi_i \xi)(1 - \eta^2) $$

Node ordering (CCW): corners 0-3, horizontal midside 4-7, vertical midside 8-11. Evaluated using 3x3 Gauss quadrature with points at +/-sqrt(3/5) and 0, weights 5/9, 8/9, 5/9.

Q4 vs Q8: Element Comparison

All three major commercial FEA codes (Abaqus, ANSYS, Nastran) recommend quadratic elements as the primary choice for general structural analysis. Q8 is the default recommendation in Abaqus (CPS8R), ANSYS (PLANE183), and Nastran (CQUAD8).

Property	Q4 (Bilinear)	Q8 (Serendipity)
Nodes	4 (corners only)	8 (4 corners + 4 midside)
Shape functions	Bilinear: 1, xi, eta, xi*eta	Biquadratic: +xi^2, eta^2, xi^2eta, xieta^2
Integration	2x2 Gauss (4 pts)	3x3 Gauss (9 pts)
Stiffness matrix	8x8	16x16
Convergence rate	O(h^2) energy norm	O(h^3) energy norm
Shear locking	Yes (full integration)	No (quadratic displacement field)
Curved boundaries	No (straight edges)	Yes (midside nodes)
Distortion sensitivity	Low	Moderate (midside position matters)
Best for	Large deformation, explicit dynamics	General analysis, stress concentrations
Abaqus equivalent	CPS4, CPS4R	CPS8, CPS8R
ANSYS equivalent	PLANE182	PLANE183

Shear Locking in Q4 Elements

Shear locking is an artificial stiffening that occurs when fully-integrated lower-order elements model bending-dominated problems. In pure bending, the shear strain should be zero, but Q4's bilinear shape functions introduce parasitic shear strains that make the element appear artificially stiff. This is why Q4 tip deflections can be off by 50%+ on coarse meshes, while Q8 converges to the correct solution with far fewer elements.

Jacobian & B-Matrix

The Jacobian maps natural to physical coordinates:

$$ J = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} = \sum_{i=1}^{4} \begin{bmatrix} \frac{\partial N_i}{\partial \xi} x_i & \frac{\partial N_i}{\partial \xi} y_i \\ \frac{\partial N_i}{\partial \eta} x_i & \frac{\partial N_i}{\partial \eta} y_i \end{bmatrix} $$

The strain-displacement matrix B relates nodal displacements to strain:

$$ B = \begin{bmatrix} \frac{\partial N_1}{\partial x} & 0 & \cdots & \frac{\partial N_4}{\partial x} & 0 \\ 0 & \frac{\partial N_1}{\partial y} & \cdots & 0 & \frac{\partial N_4}{\partial y} \\ \frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} & \cdots & \frac{\partial N_4}{\partial y} & \frac{\partial N_4}{\partial x} \end{bmatrix} $$

Element Stiffness

$$ K_e = \int_{-1}^{1} \int_{-1}^{1} B^T D B \, \det(J) \, t \, d\xi \, d\eta $$

Evaluated using 2x2 Gauss quadrature with points at +/-1/sqrt(3) and unit weights.

Penalty Method for Dirichlet BCs

Instead of eliminating rows/columns, a large penalty term alpha is added to the diagonal of constrained DOFs:

$$ K_{ii} \leftarrow K_{ii} + \alpha, \quad f_i \leftarrow f_i + \alpha \cdot u_i^{prescribed} $$

where alpha = max(1e4 * K_max, 1e8). This enforces the constraint approximately, with error proportional to 1/alpha.

Von Mises Stress

$$ \sigma_{vm} = \sqrt{\sigma_{xx}^2 - \sigma_{xx}\sigma_{yy} + \sigma_{yy}^2 + 3\sigma_{xy}^2} $$

Strain-Displacement Relations

The strain tensor is related to the displacement field through the strain-displacement matrix B:

$$ \varepsilon = \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy} \end{bmatrix} = \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{bmatrix} = B \cdot d_e $$

where u and v are displacement components in x and y directions, and d_e is the element nodal displacement vector.

Principal Stresses

Principal stresses are the maximum and minimum normal stresses at a point, occurring on planes where shear stress is zero:

$$ \sigma_{1,2} = \frac{\sigma_{xx} + \sigma_{yy}}{2} \pm \sqrt{\left(\frac{\sigma_{xx} - \sigma_{yy}}{2}\right)^2 + \sigma_{xy}^2} $$

The principal angle (orientation of principal planes) is:

$$ \tan(2\theta_p) = \frac{2\sigma_{xy}}{\sigma_{xx} - \sigma_{yy}} $$

Strain Energy

The strain energy stored in the structure is:

$$ U = \frac{1}{2} \int_\Omega \sigma : \varepsilon \, d\Omega = \frac{1}{2} d^T K d $$

For linear elastic problems, the strain energy equals the work done by external forces:

$$ U = W = \frac{1}{2} f^T d $$

Mesh Convergence (GCI Method)

The Grid Convergence Index (GCI) quantifies discretization error in mesh refinement studies:

$$ GCI = \frac{F_s}{r^p - 1} \left| \frac{f_2 - f_1}{f_1} \right| $$

where F_s is the safety factor (3.0 for two meshes), r is the refinement ratio (h_coarse/h_fine), p is the observed order of convergence, and f_1, f_2 are the fine and coarse mesh solutions.

Richardson extrapolation provides an estimate of the exact solution:

$$ f_{exact} \approx f_1 + \frac{f_1 - f_2}{r^p - 1} $$

Gauss Quadrature

Numerical integration uses Gauss quadrature. For 2x2 integration of Q4 elements:

$$ \int_{-1}^{1} \int_{-1}^{1} f(\xi, \eta) \, d\xi \, d\eta \approx \sum_{i=1}^{2} \sum_{j=1}^{2} w_i w_j f(\xi_i, \eta_j) $$

Integration points are at +/-1/sqrt(3) with unit weights (w_i = w_j = 1).

Assembly Process

Element stiffness matrices are assembled into the global system using connectivity matrices:

$$ K_{global} = \sum_{e=1}^{n_{elem}} L_e^T K_e L_e $$

where L_e is the Boolean connectivity matrix mapping element DOFs to global DOFs. The assembly uses COO format (append-only) and converts to CSR for efficient matrix-vector products.

Boundary Conditions

Dirichlet (Essential)

Prescribed displacements on boundary Gamma_u:

$$ u_i = \bar{u}_i \quad \text{on } \Gamma_u $$

Neumann (Natural)

Prescribed tractions on boundary Gamma_t:

$$ \sigma \cdot n = t \quad \text{on } \Gamma_t $$

In Cauchy, Dirichlet BCs are enforced using the penalty method, which adds a large stiffness to constrained DOFs without restructuring the matrix.

Conjugate Gradient Solver

For large sparse systems, CG solves Kx=f iteratively:

Initialize: r_0 = f - K*x_0, p_0 = r_0
For k = 0, 1, 2, ...: alpha_k = (r_k . r_k) / (p_k . K*p_k)
x_{k+1} = x_k + alpha_k * p_k
r_{k+1} = r_k - alpha_k * K*p_k
beta_k = (r_{k+1} . r_{k+1}) / (r_k . r_k)
p_{k+1} = r_{k+1} + beta_k * p_k

Cauchy uses Jacobi (diagonal) preconditioning and terminates when the residual drops below 1e-10.

ZZ Error Estimator (Zienkiewicz-Zhu)

The ZZ error estimator uses superconvergent patch recovery (SPR) to compute a smoothed stress field sigma* from the FE solution sigma_h. The element-wise error indicator is:

$$ \eta_e = \sqrt{\int_{\Omega_e} |\sigma^* - \sigma_h|^2 \, d\Omega} $$

SPR Recovery: For each node, average the Gauss-point stresses from surrounding elements using inverse-distance weighting. This produces a superconvergent stress field that converges at a higher rate than the FE solution.

Marking Strategy: "Largest first" -- sort elements by error indicator descending, mark the top elements until the sum of marked errors exceeds theta * total_error (theta = 0.5 default).

Adaptive Loop: Solve -> recover -> estimate -> mark -> refine -> repeat. The error estimator drives mesh refinement toward regions of high stress gradient, achieving similar accuracy with 10-14x fewer nodes than uniform refinement.

Red-Green h-Refinement

Marked elements undergo geometric subdivision:

Red refinement: Split a quad into 4 equal sub-quads (add 5 nodes: center + 4 edge midpoints)
Green refinement: Transition zone -- split into 2 triangles or quads to maintain conforming edges

After refinement, elements inside geometric cutouts (e.g., the circular hole) are removed, and the mesh is re-conformed.