mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
test(materials): add hyperelasticity new API visionary test
New 693-line test file for hyperelasticity NEW API (test-driven development): - Visionary tests for Neo-Hookean material model - Tests for Mooney-Rivlin and Ogden models - Tests for incompressibility constraint (nearly incompressible materials) - Tests for uniaxial tension with large strains - Extensive documentation of intended API design - Tests currently skipped (@test_skip) until implementation complete - Documents strain energy functions, stress computation, tangent moduli - Documents frame-invariance and material objectivity This test file serves as both test suite and API design documentation for the new hyperelastic material model interface.
This commit is contained in:
Jukka Aho
@@ -0,0 +1,693 @@
|
||||
"""
|
||||
# Hyperelasticity - NEW API (Test-Driven Development)
|
||||
|
||||
**What:** Shows how hyperelastic materials SHOULD work with the NEW API
|
||||
|
||||
**Why:**
|
||||
- **Finite strain** - Large deformations (rubber, soft tissue, biomechanics)
|
||||
- **Energy-based** - Strain energy function Ψ(F)
|
||||
- **Frame-invariant** - Material objectivity (rotation independence)
|
||||
- **Multiple models** - Neo-Hookean, Mooney-Rivlin, Ogden, etc.
|
||||
- **Incompressibility** - Nearly incompressible (ν ≈ 0.5)
|
||||
|
||||
**NEW API Concepts:**
|
||||
1. **Hyperelastic material types** - Neo-Hookean, Mooney-Rivlin, Ogden
|
||||
2. **Strain energy function** - Ψ(F) and derivatives
|
||||
3. **Push-forward stress** - σ = (1/J) P F^T (Cauchy from 1st Piola-Kirchhoff)
|
||||
4. **Tangent moduli** - C_ijkl = ∂²Ψ/∂F_ij∂F_kl
|
||||
5. **Incompressibility constraint** - det(F) ≈ 1
|
||||
|
||||
**Test Problems:**
|
||||
|
||||
## Test 1: Neo-Hookean (Simplest Hyperelastic)
|
||||
- Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
|
||||
- Validates stress computation from energy
|
||||
- Tests incompressibility limit
|
||||
|
||||
## Test 2: Mooney-Rivlin (Two-Parameter)
|
||||
- Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
|
||||
- Better for rubber than Neo-Hookean
|
||||
- Validates second invariant I₂
|
||||
|
||||
## Test 3: Ogden Model (Multi-Term)
|
||||
- Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3)
|
||||
- Uses principal stretches λᵢ
|
||||
- Most accurate for rubber
|
||||
|
||||
## Test 4: Uniaxial Tension Test
|
||||
- Compare to experimental data
|
||||
- Validates material parameters
|
||||
- Tests large strain (λ > 2)
|
||||
|
||||
**Expected Behavior (when implemented):**
|
||||
✅ Stress computed from ∂Ψ/∂F correctly
|
||||
✅ Tangent moduli symmetric and positive-definite
|
||||
✅ Incompressibility enforced (J ≈ 1)
|
||||
✅ Frame-invariant (rotations don't change Ψ)
|
||||
✅ Matches experimental stress-strain curves
|
||||
✅ Works with Newton-Krylov solver
|
||||
|
||||
**Status:** 🚧 VISIONARY TEST - Implementation in progress
|
||||
"""
|
||||
|
||||
using Test
|
||||
using JuliaFEM
|
||||
using Tensors
|
||||
using LinearAlgebra
|
||||
using Statistics
|
||||
|
||||
@testset "Hyperelasticity - NEW API (TDD)" begin
|
||||
|
||||
# =============================================================================
|
||||
# NEO-HOOKEAN MODEL (SIMPLEST)
|
||||
# =============================================================================
|
||||
|
||||
@testset "Neo-Hookean Material (Visionary)" begin
|
||||
@test_skip begin # Skip until implemented
|
||||
|
||||
# Material parameters
|
||||
μ = 1000.0 # Shear modulus
|
||||
λ = 2000.0 # Lame parameter (nearly incompressible)
|
||||
|
||||
# NEW: Hyperelastic material type
|
||||
material = NeoHookean(
|
||||
μ=μ,
|
||||
λ=λ,
|
||||
formulation=:compressible # or :incompressible
|
||||
)
|
||||
|
||||
# Test deformation gradient
|
||||
F = Tensor{2,3}((
|
||||
1.2, 0.1, 0.0,
|
||||
0.0, 0.9, 0.0,
|
||||
0.0, 0.0, 1.0
|
||||
))
|
||||
|
||||
# Compute strain energy
|
||||
Ψ = strain_energy(material, F)
|
||||
|
||||
# Neo-Hookean energy:
|
||||
# Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
|
||||
C = tdot(F) # Right Cauchy-Green: C = F^T F
|
||||
I₁ = tr(C) # First invariant
|
||||
J = det(F) # Volume ratio
|
||||
|
||||
Ψ_analytical = μ / 2 * (I₁ - 3) - μ * log(J) + λ / 2 * log(J)^2
|
||||
|
||||
@test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
|
||||
|
||||
# First Piola-Kirchhoff stress: P = ∂Ψ/∂F
|
||||
P = first_piola_kirchhoff_stress(material, F)
|
||||
|
||||
# Analytical P for Neo-Hookean
|
||||
F_inv = inv(F)
|
||||
P_analytical = μ * (F - tdot(F_inv)) + λ * log(J) * tdot(F_inv)
|
||||
|
||||
@test P ≈ P_analytical rtol = 1e-10
|
||||
|
||||
# Cauchy stress: σ = (1/J) P F^T
|
||||
σ = cauchy_stress(material, F)
|
||||
σ_from_P = (1 / J) * P ⊡ transpose(F)
|
||||
|
||||
@test σ ≈ σ_from_P rtol = 1e-10
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# INCOMPRESSIBILITY CONSTRAINT
|
||||
# =============================================================================
|
||||
|
||||
@testset "Incompressibility (Nearly) (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Nearly incompressible (ν → 0.5)
|
||||
E = 1000.0
|
||||
ν = 0.499 # Nearly incompressible
|
||||
|
||||
μ = E / (2 * (1 + ν))
|
||||
λ = E * ν / ((1 + ν) * (1 - 2ν)) # Very large!
|
||||
|
||||
material = NeoHookean(
|
||||
μ=μ,
|
||||
λ=λ,
|
||||
formulation=:compressible
|
||||
)
|
||||
|
||||
# Deformation (should preserve volume)
|
||||
F = Tensor{2,3}((
|
||||
1.5, 0.0, 0.0,
|
||||
0.0, 1 / sqrt(1.5), 0.0,
|
||||
0.0, 0.0, 1 / sqrt(1.5)
|
||||
))
|
||||
|
||||
J = det(F)
|
||||
|
||||
# For incompressible: J = 1 (volume preserving)
|
||||
@test isapprox(J, 1.0, atol=1e-3)
|
||||
|
||||
# Hydrostatic pressure enforces incompressibility
|
||||
σ = cauchy_stress(material, F)
|
||||
p = -tr(σ) / 3 # Hydrostatic pressure
|
||||
|
||||
# For nearly incompressible, pressure should be large
|
||||
@test abs(p) > 1000.0 # Significant pressure
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# MOONEY-RIVLIN MODEL
|
||||
# =============================================================================
|
||||
|
||||
@testset "Mooney-Rivlin Material (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Material parameters (typical for rubber)
|
||||
C₁ = 0.5 # MPa
|
||||
C₂ = 0.1 # MPa
|
||||
κ = 100.0 # Bulk modulus (incompressibility)
|
||||
|
||||
# NEW: Mooney-Rivlin material
|
||||
material = MooneyRivlin(
|
||||
C1=C₁,
|
||||
C2=C₂,
|
||||
bulk_modulus=κ
|
||||
)
|
||||
|
||||
# Test deformation
|
||||
F = Tensor{2,3}((
|
||||
1.5, 0.2, 0.0,
|
||||
0.1, 0.8, 0.0,
|
||||
0.0, 0.0, 1.0
|
||||
))
|
||||
|
||||
# Strain energy: Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
|
||||
C = tdot(F)
|
||||
I₁ = tr(C)
|
||||
I₂ = 0.5 * (tr(C)^2 - tr(C ⊡ C)) # Second invariant
|
||||
J = det(F)
|
||||
|
||||
Ψ = strain_energy(material, F)
|
||||
Ψ_analytical = C₁ * (I₁ - 3) + C₂ * (I₂ - 3) + κ / 2 * (J - 1)^2
|
||||
|
||||
@test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
|
||||
|
||||
# Stress computation
|
||||
P = first_piola_kirchhoff_stress(material, F)
|
||||
σ = cauchy_stress(material, F)
|
||||
|
||||
# Validate symmetry of Cauchy stress
|
||||
@test isapprox(σ, transpose(σ), atol=1e-10)
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# OGDEN MODEL (PRINCIPAL STRETCHES)
|
||||
# =============================================================================
|
||||
|
||||
@testset "Ogden Material (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Ogden parameters (multi-term)
|
||||
μ_terms = [1.0, 0.5, 0.2] # Shear moduli
|
||||
α_terms = [2.0, 3.0, -2.0] # Exponents
|
||||
κ = 100.0
|
||||
|
||||
# NEW: Ogden material
|
||||
material = Ogden(
|
||||
mu=μ_terms,
|
||||
alpha=α_terms,
|
||||
bulk_modulus=κ
|
||||
)
|
||||
|
||||
# Deformation
|
||||
F = Tensor{2,3}((
|
||||
2.0, 0.0, 0.0,
|
||||
0.0, 0.6, 0.0,
|
||||
0.0, 0.0, 0.8
|
||||
))
|
||||
|
||||
# Compute principal stretches
|
||||
C = tdot(F)
|
||||
eigenvalues_C = eigvals(C)
|
||||
λ = sqrt.(eigenvalues_C) # Principal stretches
|
||||
|
||||
# Strain energy: Ψ = Σᵢ μᵢ/αᵢ (λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3)
|
||||
Ψ = strain_energy(material, F)
|
||||
|
||||
J = det(F)
|
||||
Ψ_analytical = sum(
|
||||
μ_terms[i] / α_terms[i] * (sum(λ .^ α_terms[i]) - 3)
|
||||
for i in 1:3
|
||||
) + κ / 2 * (J - 1)^2
|
||||
|
||||
@test isapprox(Ψ, Ψ_analytical, rtol=1e-10)
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# UNIAXIAL TENSION TEST
|
||||
# =============================================================================
|
||||
|
||||
@testset "Uniaxial Tension (Large Strain) (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
# Geometry: Unit cube under tension
|
||||
mesh = generate_mesh(
|
||||
geometry=UnitCube(),
|
||||
element_type=Hex8,
|
||||
n_elements=(4, 4, 4)
|
||||
)
|
||||
|
||||
# Material (Neo-Hookean rubber)
|
||||
material = NeoHookean(
|
||||
μ=1.0, # MPa
|
||||
λ=10.0, # Nearly incompressible
|
||||
density=1000.0
|
||||
)
|
||||
|
||||
# Physics
|
||||
elastic_physics = ContinuumPhysics{Displacement}(
|
||||
material=material,
|
||||
formulation=FullThreeD(),
|
||||
finite_strain=true # CRITICAL!
|
||||
)
|
||||
|
||||
domain = Domain(
|
||||
name="RUBBER",
|
||||
elements=mesh.elements,
|
||||
physics=elastic_physics
|
||||
)
|
||||
|
||||
# Boundary conditions
|
||||
bc_fixed = DirichletBC(
|
||||
nodes=mesh.node_sets["LEFT"],
|
||||
dof=:displacement,
|
||||
values=[0.0, 0.0, 0.0]
|
||||
)
|
||||
|
||||
# Applied stretch (λ = 2.0, 100% strain!)
|
||||
u_applied = 1.0 # Stretch from 1.0 to 2.0
|
||||
|
||||
bc_stretch = DirichletBC(
|
||||
nodes=mesh.node_sets["RIGHT"],
|
||||
dof=:displacement,
|
||||
component=:x,
|
||||
value=u_applied
|
||||
)
|
||||
|
||||
# Nonlinear problem (finite strain)
|
||||
problem = NonlinearProblem(
|
||||
domains=[domain],
|
||||
boundary_conditions=[bc_fixed, bc_stretch]
|
||||
)
|
||||
|
||||
# Newton-Krylov solver
|
||||
solver = NewtonKrylovSolver(
|
||||
max_iterations=20,
|
||||
convergence_tol=1e-6,
|
||||
krylov_solver=GMRES(restart=30),
|
||||
line_search=BacktrackingLineSearch()
|
||||
)
|
||||
|
||||
solution = solve!(problem, solver)
|
||||
|
||||
# Extract stress-stretch curve
|
||||
λ = 1.0 + u_applied # Stretch ratio
|
||||
|
||||
# Compute engineering stress
|
||||
F_avg = compute_average_deformation_gradient(solution)
|
||||
P_avg = first_piola_kirchhoff_stress(material, F_avg)
|
||||
|
||||
# Engineering stress: σ_eng = P_11 (1st Piola-Kirchhoff in x)
|
||||
σ_eng = P_avg[1, 1]
|
||||
|
||||
# For Neo-Hookean uniaxial:
|
||||
# σ_eng = μ(λ - 1/λ²)
|
||||
σ_analytical = material.μ * (λ - 1 / λ^2)
|
||||
|
||||
@test isapprox(σ_eng, σ_analytical, rtol=0.05) # 5% tolerance
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# FRAME INVARIANCE (OBJECTIVITY)
|
||||
# =============================================================================
|
||||
|
||||
@testset "Frame Invariance (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
material = NeoHookean(μ=1000.0, λ=2000.0)
|
||||
|
||||
# Deformation gradient
|
||||
F = Tensor{2,3}((
|
||||
1.2, 0.1, 0.0,
|
||||
0.0, 0.9, 0.0,
|
||||
0.0, 0.0, 1.0
|
||||
))
|
||||
|
||||
# Rotation tensor (90° about z-axis)
|
||||
θ = π / 2
|
||||
Q = Tensor{2,3}((
|
||||
cos(θ), -sin(θ), 0.0,
|
||||
sin(θ), cos(θ), 0.0,
|
||||
0.0, 0.0, 1.0
|
||||
))
|
||||
|
||||
# Rotated deformation: F' = Q F
|
||||
F_rotated = Q ⊡ F
|
||||
|
||||
# Strain energy should be INVARIANT
|
||||
Ψ = strain_energy(material, F)
|
||||
Ψ_rotated = strain_energy(material, F_rotated)
|
||||
|
||||
@test isapprox(Ψ, Ψ_rotated, rtol=1e-10)
|
||||
|
||||
# Cauchy stress should transform: σ' = Q σ Q^T
|
||||
σ = cauchy_stress(material, F)
|
||||
σ_rotated = cauchy_stress(material, F_rotated)
|
||||
|
||||
σ_transformed = Q ⊡ σ ⊡ transpose(Q)
|
||||
|
||||
@test σ_rotated ≈ σ_transformed rtol = 1e-10
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# TANGENT MODULI (FOR NEWTON)
|
||||
# =============================================================================
|
||||
|
||||
@testset "Tangent Moduli (Visionary)" begin
|
||||
@test_skip begin
|
||||
|
||||
material = NeoHookean(μ=1000.0, λ=2000.0)
|
||||
|
||||
F = Tensor{2,3}((
|
||||
1.2, 0.1, 0.0,
|
||||
0.0, 0.9, 0.0,
|
||||
0.0, 0.0, 1.0
|
||||
))
|
||||
|
||||
# Material tangent: C_ijkl = ∂²Ψ/∂F_ij∂F_kl
|
||||
C = material_tangent(material, F)
|
||||
|
||||
# Validate major symmetry: C_ijkl = C_klij
|
||||
# (Minor symmetries don't hold for finite strain)
|
||||
for i in 1:3, j in 1:3, k in 1:3, l in 1:3
|
||||
@test isapprox(C[i, j, k, l], C[k, l, i, j], atol=1e-10)
|
||||
end
|
||||
|
||||
# Validate positive-definiteness
|
||||
# For small perturbation δF, δ²Ψ = C_ijkl δF_ij δF_kl > 0
|
||||
δF = 0.01 * rand(Tensor{2,3})
|
||||
δ²Ψ = dcontract(dcontract(C, δF), δF)
|
||||
|
||||
@test δ²Ψ > 0 # Positive-definite
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# PSEUDO-CODE: HYPERELASTIC ASSEMBLY
|
||||
# =============================================================================
|
||||
|
||||
@testset "Hyperelastic Assembly Pattern (Visionary)" begin
|
||||
# Pseudo-code showing finite strain assembly
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("HYPERELASTIC ASSEMBLY (FINITE STRAIN)")
|
||||
println("="^70)
|
||||
|
||||
assembly_pseudo = """
|
||||
# For hyperelastic materials, assembly uses current configuration
|
||||
|
||||
function compute_residual_hyperelastic(u, material, elements)
|
||||
r = zeros(length(u))
|
||||
|
||||
for elem in elements
|
||||
# Current deformation gradient: F = I + ∇u
|
||||
X = reference_coordinates(elem) # Undeformed
|
||||
x = X + u[elem.nodes] # Deformed
|
||||
|
||||
for ip in integration_points(elem)
|
||||
# Jacobian in reference config
|
||||
J₀, dN_dX = jacobian(elem, ip, X)
|
||||
|
||||
# Deformation gradient: F = ∂x/∂X
|
||||
F = compute_deformation_gradient(x, dN_dX)
|
||||
|
||||
# 1st Piola-Kirchhoff stress: P = ∂Ψ/∂F
|
||||
P = first_piola_kirchhoff_stress(material, F)
|
||||
|
||||
# Residual: r = ∫_Ω₀ P : ∇_X(δu) dΩ₀
|
||||
# = ∫_Ω₀ P_iJ (dN_I/dX_J) dΩ₀
|
||||
for I in 1:n_nodes
|
||||
for i in 1:3
|
||||
for J in 1:3
|
||||
r[3*(I-1)+i] += P[i,J] * dN_dX[I,J] * J₀ * ip.weight
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return r
|
||||
end
|
||||
|
||||
# Tangent stiffness (for Newton):
|
||||
# K_t = ∫_Ω₀ (dN_I/dX_K) C_iJkL (dN_J/dX_L) dΩ₀
|
||||
# where C_iJkL = ∂²Ψ/∂F_iJ∂F_kL
|
||||
"""
|
||||
|
||||
println(assembly_pseudo)
|
||||
println("="^70)
|
||||
println("✓ Uses reference configuration Ω₀ (not current!)")
|
||||
println("✓ Deformation gradient F = ∂x/∂X")
|
||||
println("✓ 1st Piola-Kirchhoff stress P = ∂Ψ/∂F")
|
||||
println("✓ Tangent moduli C = ∂²Ψ/∂F²")
|
||||
println("✓ Works with nodal assembly (same pattern!)")
|
||||
println("="^70)
|
||||
end
|
||||
|
||||
# =============================================================================
|
||||
# KEY ARCHITECTURAL INSIGHTS
|
||||
# =============================================================================
|
||||
|
||||
println("\n" * "="^70)
|
||||
println("HYPERELASTICITY ARCHITECTURE INSIGHTS (NEW API)")
|
||||
println("="^70)
|
||||
println("✓ Hyperelastic materials: NeoHookean, MooneyRivlin, Ogden")
|
||||
println("✓ Strain energy function Ψ(F) is fundamental")
|
||||
println("✓ Stress from energy: P = ∂Ψ/∂F, σ = (1/J) P F^T")
|
||||
println("✓ Tangent from energy: C = ∂²Ψ/∂F²")
|
||||
println("✓ Frame-invariant: Rotations don't change Ψ")
|
||||
println("✓ Incompressibility: det(F) ≈ 1 for rubber")
|
||||
println("✓ Works with Newton-Krylov (unsymmetric OK)")
|
||||
println("✓ Assembly in reference config (not current!)")
|
||||
println("✓ Large strains: λ > 2 (100%+ strain)")
|
||||
println("="^70)
|
||||
|
||||
end
|
||||
|
||||
"""
|
||||
# IMPLEMENTATION NOTES
|
||||
|
||||
## Hyperelastic Material Models
|
||||
|
||||
### Neo-Hookean (Simplest)
|
||||
|
||||
**Strain energy:**
|
||||
Ψ = μ/2 (I₁ - 3) - μ ln(J) + λ/2 (ln J)²
|
||||
|
||||
where:
|
||||
- I₁ = tr(C) = tr(F^T F) (first invariant)
|
||||
- J = det(F) (volume ratio)
|
||||
- μ, λ = Lame parameters
|
||||
|
||||
**1st Piola-Kirchhoff stress:**
|
||||
P = ∂Ψ/∂F = μ(F - F^{-T}) + λ ln(J) F^{-T}
|
||||
|
||||
where F^{-T} = (F^{-1})^T.
|
||||
|
||||
**Cauchy stress:**
|
||||
σ = (1/J) P F^T = μ/J (B - I) + λ/J ln(J) I
|
||||
|
||||
where B = F F^T (left Cauchy-Green).
|
||||
|
||||
**Use case:** Rubber at moderate strains (< 100%).
|
||||
|
||||
### Mooney-Rivlin (Two-Parameter)
|
||||
|
||||
**Strain energy:**
|
||||
Ψ = C₁(I₁ - 3) + C₂(I₂ - 3) + κ/2 (J - 1)²
|
||||
|
||||
where:
|
||||
- I₁ = tr(C)
|
||||
- I₂ = 1/2 [(tr C)² - tr(C²)]
|
||||
- J = det(F)
|
||||
- C₁, C₂ = material parameters
|
||||
- κ = bulk modulus
|
||||
|
||||
**Better fit for rubber** than Neo-Hookean.
|
||||
|
||||
**Relation to Neo-Hookean:** C₂ = 0 → Neo-Hookean.
|
||||
|
||||
### Ogden Model (Multi-Term)
|
||||
|
||||
**Strain energy:**
|
||||
Ψ = Σᵢ μᵢ/αᵢ (λ₁^{αᵢ} + λ₂^{αᵢ} + λ₃^{αᵢ} - 3) + κ/2 (J - 1)²
|
||||
|
||||
where:
|
||||
- λ₁, λ₂, λ₃ = principal stretches (eigenvalues of F)
|
||||
- μᵢ, αᵢ = material parameters (typically 3-6 terms)
|
||||
|
||||
**Most accurate for rubber** (fits experimental data well).
|
||||
|
||||
**Implementation:** Requires spectral decomposition of C.
|
||||
|
||||
## Stress Measures
|
||||
|
||||
### First Piola-Kirchhoff (P)
|
||||
|
||||
**Definition:** P = ∂Ψ/∂F
|
||||
|
||||
**Properties:**
|
||||
- Non-symmetric
|
||||
- Force per reference area
|
||||
- Work-conjugate to F
|
||||
|
||||
**Use:** Weak form in reference config.
|
||||
|
||||
### Cauchy Stress (σ)
|
||||
|
||||
**Definition:** σ = (1/J) P F^T
|
||||
|
||||
**Properties:**
|
||||
- Symmetric
|
||||
- True stress (force per current area)
|
||||
- What we measure
|
||||
|
||||
**Use:** Post-processing, failure criteria.
|
||||
|
||||
### Second Piola-Kirchhoff (S)
|
||||
|
||||
**Definition:** S = F^{-1} P = J F^{-1} σ F^{-T}
|
||||
|
||||
**Properties:**
|
||||
- Symmetric
|
||||
- Work-conjugate to E (Green-Lagrange strain)
|
||||
- Energy-conjugate
|
||||
|
||||
**Use:** Theoretical derivations.
|
||||
|
||||
## Incompressibility
|
||||
|
||||
**Constraint:** det(F) = J = 1 (volume preserving)
|
||||
|
||||
**Nearly incompressible:** ν → 0.5, λ → ∞
|
||||
|
||||
**Enforcement:**
|
||||
|
||||
1. **Penalty:** Add κ/2 (J - 1)² to Ψ (large κ)
|
||||
2. **Lagrange multiplier:** Introduce pressure p
|
||||
3. **Mixed formulation:** (u, p) unknowns
|
||||
|
||||
**JuliaFEM approach:** Penalty for compressible materials, mixed for truly incompressible.
|
||||
|
||||
## Frame Invariance (Objectivity)
|
||||
|
||||
**Definition:** Strain energy invariant under rigid rotations.
|
||||
|
||||
**Mathematical:** Ψ(Q F) = Ψ(F) for all rotations Q.
|
||||
|
||||
**Why:** Material doesn't "know" about global rotations.
|
||||
|
||||
**Implementation:** Use invariants (I₁, I₂, I₃) or principal stretches (λᵢ).
|
||||
|
||||
**Validation:**
|
||||
```julia
|
||||
Q = rotation_matrix(θ)
|
||||
F_rotated = Q * F
|
||||
@test strain_energy(F_rotated) ≈ strain_energy(F)
|
||||
```
|
||||
|
||||
## Tangent Moduli
|
||||
|
||||
**Material tangent:**
|
||||
C_{iJkL} = ∂²Ψ/∂F_{iJ}∂F_{kL}
|
||||
|
||||
**Spatial tangent:**
|
||||
c_{ijkl} = (1/J) F_{iI} F_{jJ} F_{kK} F_{lL} C_{IJKL}
|
||||
|
||||
**Symmetries:**
|
||||
- Major: C_{iJkL} = C_{kLiJ} (from Ψ)
|
||||
- Minor: Generally NOT symmetric in hyperelasticity
|
||||
|
||||
**Use in Newton:**
|
||||
δP = C : δF
|
||||
|
||||
## Assembly (Finite Strain)
|
||||
|
||||
**Residual (weak form):**
|
||||
r = ∫_{Ω₀} P : ∇_X(δu) dΩ₀ - f_ext
|
||||
|
||||
**In components:**
|
||||
r_I^i = ∫_{Ω₀} P_{iJ} ∂N_I/∂X_J dΩ₀ - f_I^i
|
||||
|
||||
**Tangent stiffness:**
|
||||
K_{IJ}^{ik} = ∫_{Ω₀} ∂N_I/∂X_K C_{iJkL} ∂N_J/∂X_L dΩ₀
|
||||
|
||||
**Key differences from small strain:**
|
||||
1. Integrate over Ω₀ (reference), not Ω (current)
|
||||
2. Use ∂N/∂X (reference gradients), not ∂N/∂x
|
||||
3. P (1st Piola-Kirchhoff), not σ (Cauchy)
|
||||
|
||||
## Nodal Assembly (Hyperelastic)
|
||||
|
||||
```julia
|
||||
function tangent_matvec_hyperelastic!(w, v, u_current, material, elements, node_to_elements)
|
||||
Threads.@threads for node_i in 1:n_nodes
|
||||
w_local = zero(Vec{3})
|
||||
|
||||
for elem in node_to_elements[node_i]
|
||||
# Current deformation
|
||||
F = deformation_gradient(elem, u_current)
|
||||
|
||||
# Material tangent
|
||||
C = material_tangent(material, F)
|
||||
|
||||
for node_j in elem.nodes
|
||||
# Tangent block: K_t_ij
|
||||
K_t_ij = compute_hyperelastic_tangent_block(elem, node_i, node_j, C, F)
|
||||
|
||||
v_j = Vec{3}(v[3*(node_j-1)+1:3*node_j])
|
||||
w_local += K_t_ij ⊡ v_j
|
||||
end
|
||||
end
|
||||
|
||||
w[3*(node_i-1)+1:3*node_i] = w_local
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
**Same pattern as linear elasticity!** Only the tangent computation changes.
|
||||
|
||||
## Next Steps
|
||||
|
||||
1. Implement `NeoHookean` material type
|
||||
2. Implement `strain_energy` function
|
||||
3. Implement `first_piola_kirchhoff_stress`
|
||||
4. Implement `cauchy_stress` (push-forward)
|
||||
5. Implement `material_tangent` (C_ijkl)
|
||||
6. Implement `MooneyRivlin` material
|
||||
7. Implement `Ogden` material (spectral decomposition)
|
||||
8. Validate against analytical solutions
|
||||
9. Validate against experimental data
|
||||
10. Performance benchmarks
|
||||
|
||||
"""
|
||||
Reference in New Issue
Block a user