Kinematics and strain
Large displacement shell behavior begins with deformation measures that remain valid when geometry changes the stiffness.
F = ∂x / ∂X, E = 1/2 (FTF - I)
Section 2 / finite element mechanics
Non-linear finite element analysis for three-dimensional shell structures.
FeToolKit is positioned around rigorous shell and stress analysis mathematics, with practical workflow tools for engineers who need traceable answers quickly.
Infinitesimal strain theory
When displacement gradients are small compared with unity, the Lagrangian and Eulerian finite strain tensors collapse to a single linear measure, the infinitesimal (Cauchy, small, or linear) strain tensor. Second-order terms are dropped, so material and spatial descriptions become indistinguishable. This is the strain measure behind linear elastic stress analysis and most first-order margin checks. Reference: Infinitesimal strain theory (Wikipedia).
Small-gradient condition: ||∇u|| ≪ 1
E = 1/2 (∇Xu + (∇Xu)T + (∇Xu)T∇Xu) ≈ 1/2 (∇Xu + (∇Xu)T)
e = 1/2 (∇xu + (∇xu)T - ∇xu (∇xu)T) ≈ 1/2 (∇xu + (∇xu)T)
E ≈ e ≈ ε = 1/2 (∇u + (∇u)T), εij = 1/2 (ui,j + uj,i)
εxx = ∂ux/∂x, εxy = 1/2 (∂ux/∂y + ∂uy/∂x), εxz = 1/2 (∂ux/∂z + ∂uz/∂x), (and cyclic)
With F = ∇u + I: ε = 1/2 (FT + F) - I
Generalized: E(m) = 1/(2m) (U2m - I) ≈ ε, e(m) = 1/(2m) (V2m - I) ≈ ε
E = 1/2 (∇Xu + (∇Xu)T + (∇Xu)T∇Xu) ≈ 1/2 (∇Xu + (∇Xu)T)
e = 1/2 (∇xu + (∇xu)T - ∇xu (∇xu)T) ≈ 1/2 (∇xu + (∇xu)T)
E ≈ e ≈ ε = 1/2 (∇u + (∇u)T), εij = 1/2 (ui,j + uj,i)
εxx = ∂ux/∂x, εxy = 1/2 (∂ux/∂y + ∂uy/∂x), εxz = 1/2 (∂ux/∂z + ∂uz/∂x), (and cyclic)
With F = ∇u + I: ε = 1/2 (FT + F) - I
Generalized: E(m) = 1/(2m) (U2m - I) ≈ ε, e(m) = 1/(2m) (V2m - I) ≈ ε
Principle of virtual work
Internal and external work terms are assembled into a residual and driven to equilibrium.
δW = ∫Ω δεTσ dΩ - ∫Ω δuTb dΩ - ∫Γ δuTt dΓ = 0
Shell resultants
Three-dimensional stress states reduce to membrane and bending quantities that map directly to aerospace checks.
N = ∫-t/2t/2 σ dz, M = ∫-t/2t/2 zσ dz
Stress recovery and margins
Field solutions become design decisions through equivalent stress, ratios, factors, and margins of safety.
σvm = sqrt(σx2 - σxσy + σy2 + 3τxy2), MS = σallow/σapplied - 1
Airy stress function
For plane problems with no body force, a single scalar potential satisfies equilibrium identically; compatibility then reduces to the biharmonic equation. The same potential can be written in polar coordinates for holes, discs, and notches, and a body-force potential adds a term when gravity or inertia loads are present.
σx = ∂2φ/∂y2, σy = ∂2φ/∂x2, τxy = -∂2φ/∂x∂y
∇4φ = ∂4φ/∂x4 + 2∂4φ/∂x2∂y2 + ∂4φ/∂y4 = 0
Polar: σr = (1/r) ∂φ/∂r + (1/r2) ∂2φ/∂θ2, σθ = ∂2φ/∂r2, τrθ = -∂/∂r [(1/r) ∂φ/∂θ]
With body-force potential V (b = -∇V): σx = ∂2φ/∂y2 + V, σy = ∂2φ/∂x2 + V, ∇4φ = -(1 - ν) ∇2V (plane stress)
∇4φ = ∂4φ/∂x4 + 2∂4φ/∂x2∂y2 + ∂4φ/∂y4 = 0
Polar: σr = (1/r) ∂φ/∂r + (1/r2) ∂2φ/∂θ2, σθ = ∂2φ/∂r2, τrθ = -∂/∂r [(1/r) ∂φ/∂θ]
With body-force potential V (b = -∇V): σx = ∂2φ/∂y2 + V, σy = ∂2φ/∂x2 + V, ∇4φ = -(1 - ν) ∇2V (plane stress)
Cauchy's stress quadric
The traction on any plane through a point follows from the stress tensor and the plane normal (Cauchy's formula). Plotting the reciprocal square root of the normal stress along every direction gives a quadric surface whose axes are the principal stress directions, and whose shape shows at a glance whether the point is in tension, compression, or saddle-type shear.
ti = σji nj, σn = niσijnj, τn = sqrt(|t|2 - σn2)
Quadric: σij xi xj = ±k2, radius vector r = k / sqrt(|σn|) along n
Principal axes: σ1x2 + σ2y2 + σ3z2 = ±k2, det(σ - λI) = 0
λ3 - I1λ2 + I2λ - I3 = 0, I1 = tr σ, I2 = 1/2 [(tr σ)2 - tr(σ2)], I3 = det σ
Property: the traction on the plane normal to OP is parallel to the normal of the quadric at P
Quadric: σij xi xj = ±k2, radius vector r = k / sqrt(|σn|) along n
Principal axes: σ1x2 + σ2y2 + σ3z2 = ±k2, det(σ - λI) = 0
λ3 - I1λ2 + I2λ - I3 = 0, I1 = tr σ, I2 = 1/2 [(tr σ)2 - tr(σ2)], I3 = det σ
Property: the traction on the plane normal to OP is parallel to the normal of the quadric at P
Finite element theory
The displacement field is approximated element by element with shape functions. Strain follows by differentiating the shape functions, and the Galerkin weak form turns equilibrium into a sparse stiffness system. Isoparametric mapping and Gauss quadrature make arbitrary element shapes practical, and Newton iteration extends the same machinery to non-linear problems.
Weak form: ∫Ω δεTσ dΩ = ∫Ω δuTb dΩ + ∫Γt δuTt dΓ
Interpolation: u(ξ) ≈ N(ξ) de, ε = B de, B = ∂N (strain-displacement operator), σ = D ε
Isoparametric: x(ξ) = Σ Ni(ξ) xi, J = ∂x/∂ξ, dΩ = det J dξ
ke = ∫ BT D B det J dξ ≈ Σg wg BT(ξg) D B(ξg) det J(ξg)
Assembly: K = Ae ke, K u = f, fe = ∫ NTb dΩ + ∫ NTt dΓ
Non-linear: KT(ui) Δu = fext - fint(ui), ui+1 = ui + Δu until ||residual|| < tol
Interpolation: u(ξ) ≈ N(ξ) de, ε = B de, B = ∂N (strain-displacement operator), σ = D ε
Isoparametric: x(ξ) = Σ Ni(ξ) xi, J = ∂x/∂ξ, dΩ = det J dξ
ke = ∫ BT D B det J dξ ≈ Σg wg BT(ξg) D B(ξg) det J(ξg)
Assembly: K = Ae ke, K u = f, fe = ∫ NTb dΩ + ∫ NTt dΓ
Non-linear: KT(ui) Δu = fext - fint(ui), ui+1 = ui + Δu until ||residual|| < tol
Bauschinger effect
Prior plastic strain in one direction lowers the yield stress on reversal. Kinematic hardening captures this by translating the yield surface with a back stress.
f = |σ - α| - σy = 0, dα = (2/3) C dεp - γ α dp
Ideal kinematic case: reverse yield σr = σmax - 2σy
Ideal kinematic case: reverse yield σr = σmax - 2σy
Elastic stress analysis
Linear elasticity couples equilibrium, strain-displacement, and Hooke's law; the finite element method discretizes this into a stiffness system.
σij,j + bi = 0, εij = 1/2 (ui,j + uj,i), σij = λδijεkk + 2μεij
K u = f, K = ∫Ω BT D B dΩ
K u = f, K = ∫Ω BT D B dΩ
Buckling analysis
Classical critical loads for columns and plates, and the linear eigenvalue form used in finite element buckling solutions.
Pcr = π2EI/(KL)2, σcr = kπ2E / [12(1 - ν2)] · (t/b)2
(K + λKg) φ = 0, Pcr = λ1Pref
(K + λKg) φ = 0, Pcr = λ1Pref
Stability analysis
Equilibrium is stable when the second variation of total potential energy is positive. Limit points and bifurcations occur where the tangent stiffness becomes singular.
δ2Π = 1/2 δuT KT δu > 0, KT = K0 + Kσ + Ku
Critical point: det KT = 0 (arc-length: ΔuTΔu + ψ2Δλ2 = Δs2)
Critical point: det KT = 0 (arc-length: ΔuTΔu + ψ2Δλ2 = Δs2)
Random vibration analysis
Broadband loading is described by a power spectral density. Response PSD follows from the transfer function, and Miles' equation gives the single-mode RMS shortcut.
Gy(f) = |H(f)|2 Gx(f), σrms = sqrt(∫ Gy(f) df)
Miles: grms = sqrt[(π/2) fn Q · ASD(fn)], design level = 3 × rms
Miles: grms = sqrt[(π/2) fn Q · ASD(fn)], design level = 3 × rms
Impact analysis
Energy balance gives a quick estimate of impact force and equivalent static load factor; explicit time integration resolves the full transient within a stability limit.
1/2 m v2 = 1/2 k δmax2 → Fmax = v sqrt(k m), n = 1 + sqrt(1 + 2h/δst)
M ü + C u̇ + fint = fext, Δt ≤ L/c, c = sqrt(E/ρ)
M ü + C u̇ + fint = fext, Δt ≤ L/c, c = sqrt(E/ρ)
Goodman fatigue analysis
Mean stress reduces the allowable alternating stress. The Goodman line bounds the safe region, and Miner's rule accumulates damage over a load spectrum.
σa/Se + σm/Ftu = 1/n, MS = 1 / (σa/Se + σm/Ftu) - 1
σa = σf' (2N)b, D = Σ ni/Ni ≤ 1
σa = σf' (2N)b, D = Σ ni/Ni ≤ 1
Fracture mechanics
The stress intensity factor characterizes the crack-tip field. Fracture occurs when it reaches toughness, and subcritical growth follows the Paris law.
K = Y σ sqrt(πa), K ≥ KIc → fracture, ac = (1/π)(KIc/Yσ)2
da/dN = C (ΔK)m, G = K2/E', J = ∫Γ (W dy - T · ∂u/∂x ds)
da/dN = C (ΔK)m, G = K2/E', J = ∫Γ (W dy - T · ∂u/∂x ds)
Elastic-plastic bending to rupture
Beam bending beyond first yield: the plastic zone spreads inward from the outer fibers until a full plastic hinge forms, and rupture is governed by the ultimate bending (modulus of rupture) stress or the outer-fiber failure strain.
ε = σ/E + 0.002 (σ/σ0.2)n (Ramberg-Osgood), κ = σy / (E yy)
Rectangular section: My = σybh2/6, M = σyb (h2/4 - yy2/3), Mp = σybh2/4 = 1.5 My
Rupture: Fb = kbFtu, Mu = Fb I/c, MS = Mu/Mapplied - 1, εouter = κh/2 ≤ εu
Rectangular section: My = σybh2/6, M = σyb (h2/4 - yy2/3), Mp = σybh2/4 = 1.5 My
Rupture: Fb = kbFtu, Mu = Fb I/c, MS = Mu/Mapplied - 1, εouter = κh/2 ≤ εu
Buckling of plates and shells
Thin plates buckle elastically and then carry more load in the post-buckled state, so plate buckling is a stiffness limit rather than an immediate failure. Cylinders and spheres are imperfection-sensitive: test buckling loads fall well below the classical value, so a knockdown factor is applied. The classical results below are the starting point for panel, skin, and fuselage checks.
Plate equation: D∇4w + Nx ∂2w/∂x2 = 0, D = Et3 / [12(1 - ν2)]
Compression: σcr = kc π2E / [12(1 - ν2)] · (t/b)2, shear: τcr = ks π2E / [12(1 - ν2)] · (t/b)2
Post-buckling effective width (von Kármán): be ≈ 1.9 t sqrt(E/σedge)
Cylinder, axial: σcl = E t / [R sqrt(3(1 - ν2))] ≈ 0.605 E t/R, design: σcr = γ σcl, γ < 1 (NASA SP-8007)
Sphere, external pressure: pcl = 2E / sqrt(3(1 - ν2)) · (t/R)2, design: pcr = γ pcl
Non-linear shell buckling: det KT(u, λ) = 0 (limit point), (K0 + λKg)φ = 0 (bifurcation)
Compression: σcr = kc π2E / [12(1 - ν2)] · (t/b)2, shear: τcr = ks π2E / [12(1 - ν2)] · (t/b)2
Post-buckling effective width (von Kármán): be ≈ 1.9 t sqrt(E/σedge)
Cylinder, axial: σcl = E t / [R sqrt(3(1 - ν2))] ≈ 0.605 E t/R, design: σcr = γ σcl, γ < 1 (NASA SP-8007)
Sphere, external pressure: pcl = 2E / sqrt(3(1 - ν2)) · (t/R)2, design: pcr = γ pcl
Non-linear shell buckling: det KT(u, λ) = 0 (limit point), (K0 + λKg)φ = 0 (bifurcation)
Critical resolved shear stress
Slip begins in a crystal when the shear stress resolved on a slip system reaches a critical value. That value is not a constant: it is the sum of an athermal part set by obstacles that heat cannot help overcome, and a thermally activated part that falls with temperature and rises with strain rate. Solute atoms, vacancies, and irradiation defects add to it as point-defect obstacles.
Schmid's law: τR = σ cosφ cosλ = τc, σy = M τc (Taylor factor M for polycrystals)
τc(T, γ̇) = τath + τ*(T, γ̇)
Thermal activation: γ̇ = γ̇0 exp[-ΔG(τ*) / kT], ΔG = ΔG0 [1 - (τ*/τ0)p]q
τ* = τ0 {1 - [ (kT/ΔG0) ln(γ̇0/γ̇) ]1/q}1/p (τ* → 0 above the athermal temperature)
Strain-rate sensitivity: m = ∂ ln σ / ∂ ln ε̇, activation volume V* = kT ∂ ln γ̇ / ∂τ
Solid solution (Fleischer): Δτss ∝ G εs3/2 c1/2, Labusch: Δτss ∝ G εs4/3 c2/3
Forest hardening: τ = τ0 + α G b sqrt(ρ), irradiation / dispersed barriers: Δτ = α G b sqrt(N d), grain size: σy = σ0 + ky d-1/2
τc(T, γ̇) = τath + τ*(T, γ̇)
Thermal activation: γ̇ = γ̇0 exp[-ΔG(τ*) / kT], ΔG = ΔG0 [1 - (τ*/τ0)p]q
τ* = τ0 {1 - [ (kT/ΔG0) ln(γ̇0/γ̇) ]1/q}1/p (τ* → 0 above the athermal temperature)
Strain-rate sensitivity: m = ∂ ln σ / ∂ ln ε̇, activation volume V* = kT ∂ ln γ̇ / ∂τ
Solid solution (Fleischer): Δτss ∝ G εs3/2 c1/2, Labusch: Δτss ∝ G εs4/3 c2/3
Forest hardening: τ = τ0 + α G b sqrt(ρ), irradiation / dispersed barriers: Δτ = α G b sqrt(N d), grain size: σy = σ0 + ky d-1/2
Phase transitions and states of matter
A phase is stable where its Gibbs free energy is lowest. At a first-order transition the free energies of two phases cross, so entropy and volume jump and latent heat is exchanged. Continuous transitions show no jump but a discontinuity in the response functions. For structural alloys the practical consequences are solidification, solid-state precipitation, martensitic transformation, and the loss of strength near melting.
Gibbs free energy: G = H - TS, dG = -S dT + V dP, equilibrium: Gα(T, P) = Gβ(T, P)
First order: S = -(∂G/∂T)P and V = (∂G/∂P)T are discontinuous, latent heat L = T ΔS
Clausius-Clapeyron: dP/dT = L / (T Δv) = Δh / (T Δv)
Gibbs phase rule: F = C - P + 2 (degrees of freedom, components, phases)
Second order: CP = -T ∂2G/∂T2, κT and α discontinuous (Ehrenfest)
Landau expansion: G(η) = G0 + a (T - Tc) η2 + b η4 + c η6, order parameter η
Fluid states: (P + a/v2)(v - b) = RT (van der Waals), critical point where (∂P/∂v)T = (∂2P/∂v2)T = 0
Lever rule for a two-phase field: fα = (Cβ - C0) / (Cβ - Cα)
First order: S = -(∂G/∂T)P and V = (∂G/∂P)T are discontinuous, latent heat L = T ΔS
Clausius-Clapeyron: dP/dT = L / (T Δv) = Δh / (T Δv)
Gibbs phase rule: F = C - P + 2 (degrees of freedom, components, phases)
Second order: CP = -T ∂2G/∂T2, κT and α discontinuous (Ehrenfest)
Landau expansion: G(η) = G0 + a (T - Tc) η2 + b η4 + c η6, order parameter η
Fluid states: (P + a/v2)(v - b) = RT (van der Waals), critical point where (∂P/∂v)T = (∂2P/∂v2)T = 0
Lever rule for a two-phase field: fα = (Cβ - C0) / (Cβ - Cα)
Elasto-plastic strain
Total strain is split into a recoverable elastic part and a permanent plastic part. A yield function marks the elastic limit, a flow rule gives the direction of plastic straining, and hardening moves or grows the yield surface. In a finite element solver this is applied as an elastic trial step followed by a return to the yield surface, and the tangent stiffness is made consistent with that update.
Decomposition: ε = εe + εp, σ = D (ε - εp)
Yield: f(σ, κ) = σvm - σy(ε̄p) ≤ 0, σvm = sqrt(3/2 s:s), s = σ - 1/3 tr(σ) I
Flow rule (Prandtl-Reuss): dεp = dλ ∂f/∂σ = dλ (3/2) s / σvm
Kuhn-Tucker: dλ ≥ 0, f ≤ 0, dλ f = 0, consistency: df = 0
Return mapping: σn+1 = σtrial - 3μ Δλ n, n = (3/2) strial / σvmtrial
Tangent: Dep = D - (D a aT D) / (h + aT D a), a = ∂f/∂σ, h = hardening modulus
Uniaxial: ε = σ/E + 0.002 (σ/σ0.2)n, true strain εt = ln(1 + e), true stress σt = s (1 + e)
Notch estimate (Neuber): σ ε = (Kt S)2 / E
Yield: f(σ, κ) = σvm - σy(ε̄p) ≤ 0, σvm = sqrt(3/2 s:s), s = σ - 1/3 tr(σ) I
Flow rule (Prandtl-Reuss): dεp = dλ ∂f/∂σ = dλ (3/2) s / σvm
Kuhn-Tucker: dλ ≥ 0, f ≤ 0, dλ f = 0, consistency: df = 0
Return mapping: σn+1 = σtrial - 3μ Δλ n, n = (3/2) strial / σvmtrial
Tangent: Dep = D - (D a aT D) / (h + aT D a), a = ∂f/∂σ, h = hardening modulus
Uniaxial: ε = σ/E + 0.002 (σ/σ0.2)n, true strain εt = ln(1 + e), true stress σt = s (1 + e)
Notch estimate (Neuber): σ ε = (Kt S)2 / E
Non-linear fracture mechanics
When the plastic zone at a crack tip is no longer small, linear elastic K loses meaning and the J-integral or crack-tip opening displacement takes over. J is path-independent for a non-linear elastic body, which makes it a good proxy for elastic-plastic behavior up to the point of significant crack growth. Tearing resistance is then described by the J-R curve and the tearing modulus, and structural integrity codes combine fracture and plastic collapse in a failure assessment diagram.
J-integral: J = ∫Γ (W dy - T · ∂u/∂x ds), W = ∫ σij dεij, J = -∂Π/∂a
Linear elastic limit: J = G = K2/E', E' = E (plane stress), E/(1 - ν2) (plane strain)
HRR field: σij = σ0 [ J / (α ε0 σ0 In r) ]1/(n+1) σ̃ij(θ, n)
Opening: δ = J / (m σy), EPRI estimate: J = Je + Jp, Jp = α ε0 σ0 c (1 - a/W) h1(a/W, n) (P/P0)n+1
Stable tearing: J ≥ JR(Δa), Tapp = (E/σ02) dJ/da ≤ Tmat
Failure assessment diagram: Kr = KI/Kmat ≤ f(Lr), Lr = σref/σy, f(Lr) = (1 - 0.14 Lr2) [0.3 + 0.7 exp(-0.65 Lr6)] (R6 Option 1)
Linear elastic limit: J = G = K2/E', E' = E (plane stress), E/(1 - ν2) (plane strain)
HRR field: σij = σ0 [ J / (α ε0 σ0 In r) ]1/(n+1) σ̃ij(θ, n)
Opening: δ = J / (m σy), EPRI estimate: J = Je + Jp, Jp = α ε0 σ0 c (1 - a/W) h1(a/W, n) (P/P0)n+1
Stable tearing: J ≥ JR(Δa), Tapp = (E/σ02) dJ/da ≤ Tmat
Failure assessment diagram: Kr = KI/Kmat ≤ f(Lr), Lr = σref/σy, f(Lr) = (1 - 0.14 Lr2) [0.3 + 0.7 exp(-0.65 Lr6)] (R6 Option 1)
Cumulative damage
Real components see variable-amplitude loading, so life is estimated by adding up the damage from each load level. Miner's rule is linear and ignores sequence, which is unconservative when high loads come first for notched parts. Non-linear rules and crack-growth integration with load-interaction models capture sequence effects, and an equivalent constant-amplitude stress gives a quick comparison between spectra.
Miner's rule: D = Σi ni / Ni, failure at D = 1 (test scatter typically 0.3 to 3)
S-N life (Basquin): σa = σf' (2Nf)b → Ni = 1/2 (σa,i/σf')1/b
Equivalent stress: σeq = [ Σ ni σa,im / Σ ni ]1/m, m = -1/b
Non-linear damage (Marco-Starkey): D = Σ (ni/Ni)xi
Continuum damage (Chaboche): dD/dN = [1 - (1 - D)β+1]α [ σa / (M (1 - D)) ]β
Crack-growth form: N = ∫a0ac da / [C (ΔK)m], Wheeler retardation: da/dN = ΦR C (ΔK)m, ΦR = [ rp,i / (aOL + rp,OL - ai) ]p
Damage tolerance margin: life factor = Ncalc / Nrequired, Nrequired includes a scatter factor
S-N life (Basquin): σa = σf' (2Nf)b → Ni = 1/2 (σa,i/σf')1/b
Equivalent stress: σeq = [ Σ ni σa,im / Σ ni ]1/m, m = -1/b
Non-linear damage (Marco-Starkey): D = Σ (ni/Ni)xi
Continuum damage (Chaboche): dD/dN = [1 - (1 - D)β+1]α [ σa / (M (1 - D)) ]β
Crack-growth form: N = ∫a0ac da / [C (ΔK)m], Wheeler retardation: da/dN = ΦR C (ΔK)m, ΦR = [ rp,i / (aOL + rp,OL - ai) ]p
Damage tolerance margin: life factor = Ncalc / Nrequired, Nrequired includes a scatter factor
Pagoda-roof rainflow fatigue counting
Rainflow counting reduces an irregular stress or strain history to a set of closed hysteresis cycles. In the pagoda-roof picture the history is turned on its side and imagined as a stack of roofs with rain flowing down from each peak and valley. A flow stops when it meets an earlier flow or runs past a reversal of greater magnitude, and each flow path gives a half-cycle whose range is its vertical extent. Matched half-cycles form full cycles, and what remains at the end is a residue of half-cycles.
Cycle parameters: Δσ = σmax - σmin, σa = Δσ/2, σm = (σmax + σmin)/2, R = σmin/σmax
Pagoda rules: start a flow at each reversal; stop when it meets a reversal of equal or greater magnitude on the opposite side, or meets a flow from above, or reaches the end of the history
Three-point (ASTM E1049) form: X = |s3 - s2|, Y = |s2 - s1|, if X ≥ Y, count a cycle of range Y and remove s1, s2
Result: cycle matrix n(Δσj, σm,k), residue counted as half-cycles (or closed by repeating the history)
Mean stress correction: Goodman σa,0 = σa / (1 - σm/Ftu), Morrow σa = (σf' - σm)(2N)b, SWT: σmax εa E = (σf')2 (2N)2b + σf' εf' E (2N)b+c
Damage: D = Σj,k njk / Njk, Njk from the corrected S-N curve
Pagoda rules: start a flow at each reversal; stop when it meets a reversal of equal or greater magnitude on the opposite side, or meets a flow from above, or reaches the end of the history
Three-point (ASTM E1049) form: X = |s3 - s2|, Y = |s2 - s1|, if X ≥ Y, count a cycle of range Y and remove s1, s2
Result: cycle matrix n(Δσj, σm,k), residue counted as half-cycles (or closed by repeating the history)
Mean stress correction: Goodman σa,0 = σa / (1 - σm/Ftu), Morrow σa = (σf' - σm)(2N)b, SWT: σmax εa E = (σf')2 (2N)2b + σf' εf' E (2N)b+c
Damage: D = Σj,k njk / Njk, Njk from the corrected S-N curve
Finite element formulations
The same continuum problem can be discretized in several ways. The displacement formulation is the workhorse: it minimizes total potential energy and gives a stiffness matrix directly. Mixed and hybrid formulations treat stress or pressure as extra unknowns, which avoids locking in nearly incompressible or thin-section problems at the cost of a larger, indefinite system. Convergence is judged with the patch test and by the order of the error as the mesh is refined.
Displacement: Π(u) = 1/2 uTK u - uTf, δΠ = 0 → K u = f
Hellinger-Reissner: ΠHR(u, σ) = ∫Ω [ σ:ε(u) - 1/2 σ:C-1:σ - b·u ] dΩ - ∫Γt t·u dΓ
Mixed system: [ -A G ; GT 0 ] {σ; u} = {0; f}, stability requires the inf-sup (LBB) condition
Dynamics: M ü + C u̇ + K u = f(t), modes: (K - ω2M) φ = 0, M = ∫ ρ NTN dΩ
Error for polynomial degree p: ||u - uh||E = O(hp), ||u - uh||L2 = O(hp+1)
Patch test: a constant-strain state must be reproduced exactly on any patch of elements
Hellinger-Reissner: ΠHR(u, σ) = ∫Ω [ σ:ε(u) - 1/2 σ:C-1:σ - b·u ] dΩ - ∫Γt t·u dΓ
Mixed system: [ -A G ; GT 0 ] {σ; u} = {0; f}, stability requires the inf-sup (LBB) condition
Dynamics: M ü + C u̇ + K u = f(t), modes: (K - ω2M) φ = 0, M = ∫ ρ NTN dΩ
Error for polynomial degree p: ||u - uh||E = O(hp), ||u - uh||L2 = O(hp+1)
Patch test: a constant-strain state must be reproduced exactly on any patch of elements
Gauss and Lobatto nodes
Element matrices are integrated numerically. Gauss-Legendre points are the interior roots of a Legendre polynomial and give the highest accuracy per point. Gauss-Lobatto points include the element ends, which is convenient for continuity and for spectral elements, where placing nodes at the Lobatto points makes the mass matrix diagonal. The choice of rule controls both accuracy and whether elements lock or develop zero-energy modes.
Gauss-Legendre (n points): ∫-11 f dξ ≈ Σi wi f(ξi), Pn(ξi) = 0, wi = 2 / [(1 - ξi2) (Pn'(ξi))2], exact to degree 2n - 1
n = 2: ξ = ±1/sqrt(3), w = 1 n = 3: ξ = 0, ±sqrt(3/5), w = 8/9, 5/9
Gauss-Lobatto (n points): (1 - ξi2) Pn-1'(ξi) = 0, wi = 2 / [n(n - 1) Pn-1(ξi)2], exact to degree 2n - 3
n = 3: ξ = -1, 0, 1, w = 1/3, 4/3, 1/3 (Simpson)
Tensor product: ∫ f dΩ ≈ ΣiΣjΣk wiwjwk f(ξi, ηj, ζk) det J
Spectral elements: Lagrange shape functions on the Lobatto nodes, Ni(ξj) = δij, lumped Mii = ρ wi det Ji
n = 2: ξ = ±1/sqrt(3), w = 1 n = 3: ξ = 0, ±sqrt(3/5), w = 8/9, 5/9
Gauss-Lobatto (n points): (1 - ξi2) Pn-1'(ξi) = 0, wi = 2 / [n(n - 1) Pn-1(ξi)2], exact to degree 2n - 3
n = 3: ξ = -1, 0, 1, w = 1/3, 4/3, 1/3 (Simpson)
Tensor product: ∫ f dΩ ≈ ΣiΣjΣk wiwjwk f(ξi, ηj, ζk) det J
Spectral elements: Lagrange shape functions on the Lobatto nodes, Ni(ξj) = δij, lumped Mii = ρ wi det Ji
Pressurized tube and plastic bending
A pressurized tube carries hoop and axial stress from the pressure, with a steep gradient through the wall once the tube is thick. First yield starts at the bore. Beyond it the plastic zone spreads outward until the whole wall yields. Bending adds a plastic hinge, and the axial and hoop stresses already present reduce the bending moment the section can carry.
Thin wall: σθ = pR/t, σz = pR/(2t) (closed ends)
Lamé: σr = A - B/r2, σθ = A + B/r2, A = (piri2 - poro2)/(ro2 - ri2), B = (pi - po) ri2ro2/(ro2 - ri2)
First yield (Tresca, bore): py = σy (ro2 - ri2) / (2ro2), full plastic: pc = σy ln(ro/ri) (von Mises: 2/sqrt(3) times this)
Bending: My = σy π(D4 - d4) / (32D), Mp = σy (D3 - d3)/6, shape factor Mp/My ≈ 1.27 (thin tube)
Axial load, thin tube: Mp,red = Mp cos(πN / 2Np), Np = σy A
Pressure interaction: σz2 - σzσθ + σθ2 = σy2 reduces the axial fiber stress available for bending
Lamé: σr = A - B/r2, σθ = A + B/r2, A = (piri2 - poro2)/(ro2 - ri2), B = (pi - po) ri2ro2/(ro2 - ri2)
First yield (Tresca, bore): py = σy (ro2 - ri2) / (2ro2), full plastic: pc = σy ln(ro/ri) (von Mises: 2/sqrt(3) times this)
Bending: My = σy π(D4 - d4) / (32D), Mp = σy (D3 - d3)/6, shape factor Mp/My ≈ 1.27 (thin tube)
Axial load, thin tube: Mp,red = Mp cos(πN / 2Np), Np = σy A
Pressure interaction: σz2 - σzσθ + σθ2 = σy2 reduces the axial fiber stress available for bending
Kirchhoff stress and Kirchhoff plate laws
In large-deformation analysis the Kirchhoff stress is the Cauchy stress weighted by the volume change, and it is the stress that pairs naturally with the rate of deformation in elasto-plastic updates. For thin plates the Kirchhoff-Love hypotheses say that normals stay straight, unstretched, and perpendicular to the mid-surface. They give the plate equation, the moment-curvature laws, and the Kirchhoff effective shear at a free edge.
Kirchhoff stress: τ = J σ, J = det F, P = τ F-T, S = F-1 τ F-T, τ:d = S:É (work conjugate)
Kirchhoff-Love kinematics: u = u0 - z ∂w/∂x, v = v0 - z ∂w/∂y, ε = ε0 + zκ, κ = { -w,xx; -w,yy; -2w,xy }
Moments: Mx = -D (w,xx + ν w,yy), My = -D (w,yy + ν w,xx), Mxy = -D (1 - ν) w,xy, D = Et3/[12(1 - ν2)]
Equilibrium: Mx,xx + 2Mxy,xy + My,yy = -q → D ∇4w = q
Shear: Qx = -D ∂(∇2w)/∂x, free edge x = a: Mx = 0, Vx = Qx + ∂Mxy/∂y = 0 (Kirchhoff condition)
Corner force of a simply supported plate: R = 2 Mxy
Kirchhoff-Love kinematics: u = u0 - z ∂w/∂x, v = v0 - z ∂w/∂y, ε = ε0 + zκ, κ = { -w,xx; -w,yy; -2w,xy }
Moments: Mx = -D (w,xx + ν w,yy), My = -D (w,yy + ν w,xx), Mxy = -D (1 - ν) w,xy, D = Et3/[12(1 - ν2)]
Equilibrium: Mx,xx + 2Mxy,xy + My,yy = -q → D ∇4w = q
Shear: Qx = -D ∂(∇2w)/∂x, free edge x = a: Mx = 0, Vx = Qx + ∂Mxy/∂y = 0 (Kirchhoff condition)
Corner force of a simply supported plate: R = 2 Mxy
Isoparametric flexure
In an isoparametric bending element, geometry, deflection, and rotations are interpolated with the same shape functions. Mindlin-Reissner kinematics keep rotation independent of slope, which lets the element carry transverse shear and use simple C0 shape functions. Thin plates then cause shear locking, which is cured by reduced or selective integration, assumed shear strain (MITC), or stabilization against hourglass modes.
Mapping: x = Σ Ni(ξ, η) xi, w = Σ Ni wi, θ = Σ Ni θi, J = [ ∂x/∂ξ ∂y/∂ξ ; ∂x/∂η ∂y/∂η ]
Gradients: {∂N/∂x; ∂N/∂y} = J-1 {∂N/∂ξ; ∂N/∂η}, dΩ = det J dξ dη
Mindlin strains: κ = Lb θ, γ = ∇w - θ
ke = ∫ BbT Db Bb det J dΩ + ∫ BsT Ds Bs det J dΩ
Db = [Et3/12(1 - ν2)] [1 ν 0; ν 1 0; 0 0 (1 - ν)/2], Ds = k G t I, k = 5/6
Locking: shear stiffness ∝ t, bending stiffness ∝ t3, so as t → 0 the shear term over-constrains; Q4 remedy: 2×2 Gauss for bending, 1×1 for shear
Beam analog (Timoshenko): Φ = 12EI / (k G A L2), shear-flexible stiffness reduced by 1/(1 + Φ)
Gradients: {∂N/∂x; ∂N/∂y} = J-1 {∂N/∂ξ; ∂N/∂η}, dΩ = det J dξ dη
Mindlin strains: κ = Lb θ, γ = ∇w - θ
ke = ∫ BbT Db Bb det J dΩ + ∫ BsT Ds Bs det J dΩ
Db = [Et3/12(1 - ν2)] [1 ν 0; ν 1 0; 0 0 (1 - ν)/2], Ds = k G t I, k = 5/6
Locking: shear stiffness ∝ t, bending stiffness ∝ t3, so as t → 0 the shear term over-constrains; Q4 remedy: 2×2 Gauss for bending, 1×1 for shear
Beam analog (Timoshenko): Φ = 12EI / (k G A L2), shear-flexible stiffness reduced by 1/(1 + Φ)
Double and Gaussian curvature analysis
The Gaussian curvature of a surface is the product of its two principal curvatures. It decides how a shell carries load. A cylinder or cone has zero Gaussian curvature and can be unrolled flat, so it bends easily. A dome has positive curvature and a saddle has negative curvature, and neither can change shape without stretching the middle surface. That is why doubly curved shells are stiff, why they carry load mainly in membrane action, and why they are imperfection-sensitive when buckling.
Principal curvatures: det(II - κ I) = 0 → κ1, κ2, K = κ1κ2 = (LN - M2) / (EG - F2), H = (κ1 + κ2)/2
K > 0 synclastic (dome), K < 0 anticlastic (saddle), K = 0 developable (cylinder, cone)
Theorema Egregium: K depends only on the first fundamental form, so bending a K ≠ 0 surface requires stretching it
Gauss-Bonnet: ∫A K dA + ∮ kg ds = 2πχ
Membrane equilibrium (Laplace): N1/R1 + N2/R2 = p, sphere: N = pR/2, cylinder: Nθ = pR, Nx = pR/2
Shallow shell (Donnell-Mushtari-Vlasov): D∇4w - ∇k2F = q, (1/Eh) ∇4F + ∇k2w = 0, ∇k2 = ky ∂2/∂x2 + kx ∂2/∂y2
K > 0 synclastic (dome), K < 0 anticlastic (saddle), K = 0 developable (cylinder, cone)
Theorema Egregium: K depends only on the first fundamental form, so bending a K ≠ 0 surface requires stretching it
Gauss-Bonnet: ∫A K dA + ∮ kg ds = 2πχ
Membrane equilibrium (Laplace): N1/R1 + N2/R2 = p, sphere: N = pR/2, cylinder: Nθ = pR, Nx = pR/2
Shallow shell (Donnell-Mushtari-Vlasov): D∇4w - ∇k2F = q, (1/Eh) ∇4F + ∇k2w = 0, ∇k2 = ky ∂2/∂x2 + kx ∂2/∂y2
Super-element finite element analysis
A super-element replaces a large part of the model with a smaller equivalent defined at its boundary. Interior degrees of freedom are condensed out, so the rest of the model sees only the interface. Static condensation is exact for static loads. Component mode synthesis adds interior vibration modes so that dynamic behavior is also captured. This is how large airframe and engine models are split between teams and assembled into a global model.
Partition: [ Kaa Kao ; Koa Koo ] {ua; uo} = {fa; fo} (a = retained boundary, o = omitted interior)
Static condensation (Guyan): K* = Kaa - Kao Koo-1 Koa, f* = fa - Kao Koo-1 fo
Recovery: uo = Koo-1 (fo - Koa ua), reduced mass: M* = TT M T, T = [ I ; -Koo-1Koa ]
Craig-Bampton: {ua; uo} = [ I 0 ; Ψ Φk ] {ua; q}, Ψ = -Koo-1Koa (constraint modes)
Fixed-interface modes: (Koo - ω2Moo) Φk = 0, truncated to the lowest k modes
Reduced system: KCB = [ K* 0 ; 0 Ω2 ], MCB = [ M̂aa M̂aq ; M̂qa I ]
Assembly: Kglobal = Σs LsT Ks Ls + Kresidual, interface compatibility ua(1) = ua(2)
Static condensation (Guyan): K* = Kaa - Kao Koo-1 Koa, f* = fa - Kao Koo-1 fo
Recovery: uo = Koo-1 (fo - Koa ua), reduced mass: M* = TT M T, T = [ I ; -Koo-1Koa ]
Craig-Bampton: {ua; uo} = [ I 0 ; Ψ Φk ] {ua; q}, Ψ = -Koo-1Koa (constraint modes)
Fixed-interface modes: (Koo - ω2Moo) Φk = 0, truncated to the lowest k modes
Reduced system: KCB = [ K* 0 ; 0 Ω2 ], MCB = [ M̂aa M̂aq ; M̂qa I ]
Assembly: Kglobal = Σs LsT Ks Ls + Kresidual, interface compatibility ua(1) = ua(2)
Shock and impact
Shock loads are short and broadband, so the structural response is judged with stress-wave and energy methods rather than a static load. A shock response spectrum gives the peak response of a family of single-degree-of-freedom oscillators to the same input, and is the usual way of specifying equipment shock environments. Drop and impact events are bounded by momentum and energy balances, and impact on stiff bodies sends stress waves whose size depends on the material impedances.
Base shock on an oscillator: z¨ + 2ζωn ż + ωn2 z = -y¨(t), SRS(fn) = maxt |absolute acceleration|
Half-sine pulse of duration T: peak amplification ≈ 1.77 near fnT ≈ 0.8 (undamped), tends to the pulse peak at high fn
Drop: v = sqrt(2 g h), impulse: Favg Δt = m Δv, energy: 1/2 m v2 = 1/2 k δmax2, n = 1 + sqrt(1 + 2h/δst)
Stress wave: c = sqrt(E/ρ), σ = ρ c v, impedance Z = ρ c
Interface: σt/σi = 2Z2 / (Z1 + Z2), σr/σi = (Z2 - Z1) / (Z1 + Z2)
Two bodies: e = (v2' - v1') / (v1 - v2), fully plastic: v' = m1v1 / (m1 + m2), energy lost = 1/2 [m1m2/(m1 + m2)] v2
Strain-rate effect (Johnson-Cook): σ = (A + Bεn) [1 + C ln(ε̇/ε̇0)] [1 - T*m]
Explicit integration: Δt ≤ Lmin/c
Half-sine pulse of duration T: peak amplification ≈ 1.77 near fnT ≈ 0.8 (undamped), tends to the pulse peak at high fn
Drop: v = sqrt(2 g h), impulse: Favg Δt = m Δv, energy: 1/2 m v2 = 1/2 k δmax2, n = 1 + sqrt(1 + 2h/δst)
Stress wave: c = sqrt(E/ρ), σ = ρ c v, impedance Z = ρ c
Interface: σt/σi = 2Z2 / (Z1 + Z2), σr/σi = (Z2 - Z1) / (Z1 + Z2)
Two bodies: e = (v2' - v1') / (v1 - v2), fully plastic: v' = m1v1 / (m1 + m2), energy lost = 1/2 [m1m2/(m1 + m2)] v2
Strain-rate effect (Johnson-Cook): σ = (A + Bεn) [1 + C ln(ε̇/ε̇0)] [1 - T*m]
Explicit integration: Δt ≤ Lmin/c
Hertzian contact
Hertz theory gives the contact area, pressure, and subsurface stress when two smooth elastic bodies are pressed together. The contact grows more slowly than the load, the pressure peaks at the center of the contact, and the maximum shear stress lies below the surface. This is the basis for bearing, gear, roller, and pin-on-plate checks, and for estimating the force and duration of an elastic impact.
Effective modulus and radius: 1/E* = (1 - ν12)/E1 + (1 - ν22)/E2, 1/R = 1/R1 + 1/R2
Sphere: a = (3FR / 4E*)1/3, δ = a2/R, F = (4/3) E* sqrt(R) δ3/2
Pressure: p(r) = p0 sqrt(1 - r2/a2), p0 = 3F / (2πa2) = (6FE*2 / π3R2)1/3
Maximum shear (ν = 0.3): τmax ≈ 0.31 p0 at depth z ≈ 0.48 a, first yield at p0 ≈ 1.6 σy
Line contact (per unit length F'): b = sqrt(4F'R / πE*), p0 = 2F' / (πb) = sqrt(F'E* / πR), τmax ≈ 0.30 p0 at z ≈ 0.78 b
Elastic impact: 1/2 m v2 = (8/15) E* sqrt(R) δmax5/2 → δmax = [15 m v2 / (16 E* sqrt(R))]2/5, Fmax = (4/3) E* sqrt(R) δmax3/2, tc ≈ 2.94 δmax / v
Sphere: a = (3FR / 4E*)1/3, δ = a2/R, F = (4/3) E* sqrt(R) δ3/2
Pressure: p(r) = p0 sqrt(1 - r2/a2), p0 = 3F / (2πa2) = (6FE*2 / π3R2)1/3
Maximum shear (ν = 0.3): τmax ≈ 0.31 p0 at depth z ≈ 0.48 a, first yield at p0 ≈ 1.6 σy
Line contact (per unit length F'): b = sqrt(4F'R / πE*), p0 = 2F' / (πb) = sqrt(F'E* / πR), τmax ≈ 0.30 p0 at z ≈ 0.78 b
Elastic impact: 1/2 m v2 = (8/15) E* sqrt(R) δmax5/2 → δmax = [15 m v2 / (16 E* sqrt(R))]2/5, Fmax = (4/3) E* sqrt(R) δmax3/2, tc ≈ 2.94 δmax / v
1Idealize geometry, material, boundary conditions, and load cases.
2Assemble element stiffness, stress resultants, and global residuals.
3Iterate with tangent stiffness until force and displacement norms converge.
4Recover stresses, compute margins, and document assumptions.