From 5b355cae6cbd1cd873414401a70747c7ff8fa0dd Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:27:51 +0300 Subject: [PATCH] chore(test): delete plasticity new-API experiment Remove unmaintained plasticity API scaffold tests. - Drop `test/materials/test_plasticity_new_api.jl`. --- test/materials/test_plasticity_new_api.jl | 747 ---------------------- 1 file changed, 747 deletions(-) delete mode 100644 test/materials/test_plasticity_new_api.jl diff --git a/test/materials/test_plasticity_new_api.jl b/test/materials/test_plasticity_new_api.jl deleted file mode 100644 index ced58ce..0000000 --- a/test/materials/test_plasticity_new_api.jl +++ /dev/null @@ -1,747 +0,0 @@ -""" -# Plasticity - NEW API (Test-Driven Development) - -**What:** Shows how plasticity models SHOULD work with the NEW API - -**Why:** -- **Permanent deformation** - Irreversible (metals, soils) -- **Yield criterion** - von Mises, Tresca, Drucker-Prager -- **Hardening** - Isotropic, kinematic, mixed -- **Rate-independence** - Path-independent (classical plasticity) -- **History-dependent** - Internal state variables - -**NEW API Concepts:** -1. **Plastic material types** - J2Plasticity, DruckerPrager -2. **Yield function** - f(σ, α) ≤ 0 (elastic domain) -3. **Flow rule** - Plastic strain rate direction -4. **Hardening laws** - Isotropic (expanding yield surface), kinematic (translation) -5. **Return mapping** - Radial return, closest point projection - -**Test Problems:** - -## Test 1: J2 Plasticity (von Mises) -- Yield: f = √(3J₂) - σ_y(ε_p) -- Isotropic hardening -- Validates return mapping algorithm - -## Test 2: Kinematic Hardening -- Backstress α (yield surface translates) -- Armstrong-Frederick model -- Validates ratcheting behavior - -## Test 3: Perfect Plasticity -- No hardening: σ_y = constant -- Validates elastic-perfectly plastic -- Tests limit load - -## Test 4: Cyclic Loading (Bauschinger Effect) -- Load → Unload → Reverse load -- Validates kinematic hardening -- Tests hysteresis loop - -**Expected Behavior (when implemented):** -✅ Yield criterion correctly evaluated -✅ Elastic-plastic split accurate -✅ Return mapping converges -✅ Hardening modulus computed correctly -✅ Consistent tangent for Newton -✅ Path-independence validated - -**Status:** 🚧 VISIONARY TEST - Implementation in progress -""" - -using Test -using JuliaFEM -using Tensors -using LinearAlgebra -using Statistics - -@testset "Plasticity - NEW API (TDD)" begin - - # ============================================================================= - # J2 PLASTICITY (VON MISES) - # ============================================================================= - - @testset "J2 Plasticity - Isotropic Hardening (Visionary)" begin - @test_skip begin # Skip until implemented - - # Material parameters - E = 200e3 # Young's modulus (MPa) - ν = 0.3 # Poisson's ratio - σ_y0 = 250.0 # Initial yield stress (MPa) - H = 2000.0 # Hardening modulus (MPa) - - # NEW: J2 plasticity material - material = J2Plasticity( - E=E, - ν=ν, - yield_stress=σ_y0, - hardening=IsotropicHardening(H=H), - hardening_law=:linear # or :exponential, :voce - ) - - # Strain history (uniaxial tension) - ε_max = 0.005 # 0.5% total strain - n_steps = 100 - ε_history = range(0, ε_max, length=n_steps) - - # Strain tensor (uniaxial) - σ_history = [] - ε_p_history = [] - - # Internal state - state = PlasticState( - ε_p=zero(SymmetricTensor{2,3}), # Plastic strain - ε_p_eq=0.0, # Equivalent plastic strain - α=zero(SymmetricTensor{2,3}) # Backstress (if kinematic) - ) - - for ε in ε_history - # Strain tensor (uniaxial tension in x) - ε_total = SymmetricTensor{2,3}(( - ε, 0.0, 0.0, - 0.0, -ν * ε, 0.0, - 0.0, 0.0, -ν * ε - )) - - # Compute stress (with return mapping) - σ, state_new = compute_stress(material, ε_total, state) - - push!(σ_history, σ[1, 1]) # Axial stress - push!(ε_p_history, state_new.ε_p_eq) - - state = state_new - end - - # Validate elastic region - ε_elastic = σ_y0 / E - elastic_indices = findall(ε_history .<= ε_elastic) - - for i in elastic_indices - # Elastic: σ = E ε - @test isapprox(σ_history[i], E * ε_history[i], rtol=0.01) - @test ε_p_history[i] == 0.0 - end - - # Validate plastic region - plastic_indices = findall(ε_history .> ε_elastic) - - for i in plastic_indices - # Plastic: σ_y(ε_p) = σ_y0 + H ε_p - ε_p = ε_p_history[i] - σ_y_current = σ_y0 + H * ε_p - - # Stress should be at yield - @test isapprox(σ_history[i], σ_y_current, rtol=0.01) - end - - end - end - - # ============================================================================= - # RETURN MAPPING ALGORITHM - # ============================================================================= - - @testset "Radial Return Mapping (Visionary)" begin - @test_skip begin - - material = J2Plasticity( - E=200e3, - ν=0.3, - yield_stress=250.0, - hardening=IsotropicHardening(H=2000.0) - ) - - # Trial elastic step (exceed yield) - ε_trial = SymmetricTensor{2,3}(( - 0.003, 0.001, 0.0, - 0.001, 0.002, 0.0, - 0.0, 0.0, 0.0 - )) - - state = PlasticState( - ε_p=zero(SymmetricTensor{2,3}), - ε_p_eq=0.0, - α=zero(SymmetricTensor{2,3}) - ) - - # Elastic predictor - σ_trial = elastic_stress(material, ε_trial - state.ε_p) - - # Yield function - s_trial = dev(σ_trial) # Deviatoric stress - q_trial = sqrt(1.5 * dcontract(s_trial, s_trial)) # von Mises stress - - f_trial = q_trial - material.yield_stress - - if f_trial > 0 - # Plastic: Return mapping required - σ, state_new = return_mapping(material, σ_trial, state) - - # Validate yield criterion satisfied - s = dev(σ) - q = sqrt(1.5 * dcontract(s, s)) - σ_y_current = material.yield_stress + material.H * state_new.ε_p_eq - - @test isapprox(q, σ_y_current, atol=1e-6) - - # Validate plastic strain increased - @test state_new.ε_p_eq > state.ε_p_eq - - else - # Elastic: No return mapping - @test f_trial <= 0 - end - - end - end - - # ============================================================================= - # KINEMATIC HARDENING (ARMSTRONG-FREDERICK) - # ============================================================================= - - @testset "Kinematic Hardening (Visionary)" begin - @test_skip begin - - # Material with kinematic hardening - material = J2Plasticity( - E=200e3, - ν=0.3, - yield_stress=250.0, - hardening=KinematicHardening( - C=5000.0, # Kinematic hardening modulus - γ=50.0 # Armstrong-Frederick parameter - ), - mixed_hardening=false - ) - - # Cyclic loading: tension → compression - ε_max = 0.005 - n_cycles = 3 - - ε_history = [] - σ_history = [] - - state = PlasticState( - ε_p=zero(SymmetricTensor{2,3}), - ε_p_eq=0.0, - α=zero(SymmetricTensor{2,3}) # Backstress - ) - - for cycle in 1:n_cycles - # Tension - for ε in range(0, ε_max, length=50) - ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) - σ, state = compute_stress(material, ε_tensor, state) - - push!(ε_history, ε) - push!(σ_history, σ[1, 1]) - end - - # Compression - for ε in range(ε_max, -ε_max, length=100) - ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) - σ, state = compute_stress(material, ε_tensor, state) - - push!(ε_history, ε) - push!(σ_history, σ[1, 1]) - end - - # Back to tension - for ε in range(-ε_max, 0, length=50) - ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) - σ, state = compute_stress(material, ε_tensor, state) - - push!(ε_history, ε) - push!(σ_history, σ[1, 1]) - end - end - - # Validate Bauschinger effect - # Yield stress in compression < initial yield - σ_y_compression = minimum(σ_history[ε_history.<0]) - @test abs(σ_y_compression) < material.yield_stress - - # Validate hysteresis loop closes - # (For stabilized cycle) - @test length(ε_history) > 0 - - end - end - - # ============================================================================= - # PERFECT PLASTICITY (NO HARDENING) - # ============================================================================= - - @testset "Perfect Plasticity (Visionary)" begin - @test_skip begin - - # No hardening: H = 0 - material = J2Plasticity( - E=200e3, - ν=0.3, - yield_stress=250.0, - hardening=NoHardening() # H = 0 - ) - - # Large strain (well into plastic) - ε_max = 0.01 # 1% strain - ε_history = range(0, ε_max, length=100) - - σ_history = [] - state = PlasticState( - ε_p=zero(SymmetricTensor{2,3}), - ε_p_eq=0.0, - α=zero(SymmetricTensor{2,3}) - ) - - for ε in ε_history - ε_tensor = SymmetricTensor{2,3}((ε, 0.0, 0.0, 0.0, 0.0, 0.0)) - σ, state = compute_stress(material, ε_tensor, state) - push!(σ_history, σ[1, 1]) - end - - # After yield, stress should be constant - ε_yield = material.yield_stress / material.E - plastic_indices = findall(ε_history .> ε_yield) - - σ_plastic = σ_history[plastic_indices] - - # All plastic stresses ≈ σ_y (no hardening!) - @test all(isapprox.(σ_plastic, material.yield_stress, rtol=0.01)) - - end - end - - # ============================================================================= - # CONSISTENT TANGENT (FOR NEWTON) - # ============================================================================= - - @testset "Consistent Tangent (Visionary)" begin - @test_skip begin - - material = J2Plasticity( - E=200e3, - ν=0.3, - yield_stress=250.0, - hardening=IsotropicHardening(H=2000.0) - ) - - # Strain state (plastic) - ε = SymmetricTensor{2,3}((0.003, 0.001, 0.0, 0.001, 0.002, 0.0)) - state = PlasticState( - ε_p=SymmetricTensor{2,3}((0.001, 0.0, 0.0, 0.0, 0.0, 0.0)), - ε_p_eq=0.001, - α=zero(SymmetricTensor{2,3}) - ) - - # Compute stress and tangent - σ, state_new, C_ep = compute_stress_tangent(material, ε, state) - - # Validate tangent via finite difference - δε = 1e-8 - - for i in 1:6 # Voigt notation - ε_pert = ε + δε * basis_symmetric_tensor(i) - σ_pert, _ = compute_stress(material, ε_pert, state) - - dσ_numerical = (σ_pert - σ) / δε - dσ_tangent = C_ep ⊡ basis_symmetric_tensor(i) - - @test isapprox(dσ_numerical, dσ_tangent, rtol=0.01) - end - - # Validate symmetry (major) - for i in 1:6, j in 1:6 - @test isapprox(C_ep[i, j], C_ep[j, i], atol=1e-10) - end - - end - end - - # ============================================================================= - # MULTI-AXIAL LOADING - # ============================================================================= - - @testset "Multi-Axial Loading (Visionary)" begin - @test_skip begin - - material = J2Plasticity( - E=200e3, - ν=0.3, - yield_stress=250.0, - hardening=IsotropicHardening(H=2000.0) - ) - - # Combined tension + shear - ε_axial_max = 0.003 - ε_shear_max = 0.002 - - n_steps = 100 - σ_history = [] - - state = PlasticState( - ε_p=zero(SymmetricTensor{2,3}), - ε_p_eq=0.0, - α=zero(SymmetricTensor{2,3}) - ) - - for i in 1:n_steps - # Proportional loading - ε_axial = ε_axial_max * i / n_steps - ε_shear = ε_shear_max * i / n_steps - - ε = SymmetricTensor{2,3}(( - ε_axial, ε_shear, 0.0, - ε_shear, 0.0, 0.0, - 0.0, 0.0, 0.0 - )) - - σ, state = compute_stress(material, ε, state) - push!(σ_history, σ) - end - - # Validate von Mises yield criterion - for σ in σ_history - s = dev(σ) - q = sqrt(1.5 * dcontract(s, s)) - σ_y_current = material.yield_stress + material.H * state.ε_p_eq - - # Should be at or below yield - @test q <= σ_y_current + 1e-6 - end - - end - end - - # ============================================================================= - # PSEUDO-CODE: PLASTICITY INTEGRATION - # ============================================================================= - - @testset "Plasticity Integration Pattern (Visionary)" begin - # Pseudo-code showing return mapping - - println("\n" * "="^70) - println("PLASTICITY INTEGRATION (RETURN MAPPING)") - println("="^70) - - integration_pseudo = """ - # Return mapping algorithm (radial return for J2) - - function compute_stress_plastic(material, ε_total, state_old) - # 1. Elastic predictor - ε_elastic_trial = ε_total - state_old.ε_p - σ_trial = C_elastic ⊡ ε_elastic_trial - - # 2. Check yield - s_trial = dev(σ_trial) # Deviatoric - q_trial = sqrt(1.5 * s_trial : s_trial) # von Mises - - σ_y = material.σ_y0 + H * state_old.ε_p_eq - f_trial = q_trial - σ_y - - if f_trial <= 0 - # Elastic: Accept trial state - return σ_trial, state_old - end - - # 3. Plastic corrector (return mapping) - # Solve for Δλ (plastic multiplier) - # f = q - σ_y(ε_p + Δλ) = 0 - - # Newton iteration - Δλ = 0.0 - for iter in 1:max_iter - σ_y_current = material.σ_y0 + H * (state_old.ε_p_eq + Δλ) - q_current = q_trial - 3*G*Δλ # G = shear modulus - - f = q_current - σ_y_current - - if abs(f) < tol - break - end - - # Derivative: df/dΔλ - df_dΔλ = -3*G - H - - # Update - Δλ -= f / df_dΔλ - end - - # 4. Update stress and state - n = s_trial / norm(s_trial) # Flow direction - - σ = σ_trial - 2*G*Δλ * n - ε_p_new = state_old.ε_p + Δλ * n - ε_p_eq_new = state_old.ε_p_eq + Δλ - - state_new = PlasticState(ε_p_new, ε_p_eq_new, state_old.α) - - return σ, state_new - end - """ - - println(integration_pseudo) - println("="^70) - println("✓ Elastic predictor: Assume elastic step") - println("✓ Check yield: f(σ_trial) ≤ 0?") - println("✓ Return mapping: Project back to yield surface") - println("✓ Newton iteration: Solve for plastic multiplier Δλ") - println("✓ Update state: ε_p, ε_p_eq, α") - println("="^70) - end - - # ============================================================================= - # KEY ARCHITECTURAL INSIGHTS - # ============================================================================= - - println("\n" * "="^70) - println("PLASTICITY ARCHITECTURE INSIGHTS (NEW API)") - println("="^70) - println("✓ J2 plasticity: von Mises yield, isotropic/kinematic hardening") - println("✓ Yield function: f(σ, α) = √(3J₂) - σ_y(ε_p)") - println("✓ Return mapping: Radial return, closest point projection") - println("✓ Consistent tangent: C_ep for Newton quadratic convergence") - println("✓ Internal state: ε_p, ε_p_eq, α (per integration point!)") - println("✓ Isotropic hardening: Yield surface expands") - println("✓ Kinematic hardening: Yield surface translates (Bauschinger)") - println("✓ Path-independent: Same final state for same strain path") - println("✓ Works with Newton-Krylov (tangent from return mapping)") - println("="^70) - -end - -""" -# IMPLEMENTATION NOTES - -## J2 Plasticity (von Mises) - -### Yield Function - -**Definition:** -f(σ, ε_p) = √(3J₂) - σ_y(ε_p) - -where: -- J₂ = (1/2) s:s (second deviatoric invariant) -- s = σ - (1/3)tr(σ)I (deviatoric stress) -- σ_y(ε_p) = yield stress (function of plastic strain) - -**Equivalent form:** -f = q - σ_y - -where q = √(3J₂) = von Mises stress. - -**Elastic domain:** f ≤ 0 - -**Yield surface:** f = 0 - -### Flow Rule - -**Associative plasticity:** Plastic strain rate direction = yield gradient - -ε̇_p = λ̇ ∂f/∂σ = λ̇ (3/2) s/q = λ̇ n - -where: -- λ̇ = plastic multiplier (rate) -- n = (3/2) s/q = flow direction (unit deviatoric) - -**Properties:** -- Incompressible: tr(ε̇_p) = 0 (volume preserving) -- Radial: ε̇_p ∝ s (proportional to deviatoric stress) - -### Hardening Laws - -**Isotropic (linear):** -σ_y(ε_p) = σ_y0 + H ε_p_eq - -where: -- σ_y0 = initial yield stress -- H = hardening modulus -- ε_p_eq = ∫ √(2/3 ε̇_p:ε̇_p) dt = equivalent plastic strain - -**Isotropic (exponential/Voce):** -σ_y(ε_p) = σ_∞ - (σ_∞ - σ_y0) exp(-b ε_p_eq) - -Saturates to σ_∞. - -**Kinematic (Armstrong-Frederick):** -α̇ = C ε̇_p - γ α λ̇ - -where: -- α = backstress (2nd order tensor) -- C = kinematic hardening modulus -- γ = recall parameter - -**Modified yield:** -f = √(3/2 (s-α):(s-α)) - σ_y - -### Return Mapping Algorithm - -**Problem:** Given ε_{n+1}, find σ_{n+1} and state_{n+1}. - -**Elastic predictor:** -``` -ε_e_trial = ε_{n+1} - ε_p_n -σ_trial = C_elastic : ε_e_trial -``` - -**Check yield:** -``` -f_trial = q_trial - σ_y(ε_p_eq_n) -``` - -**If f_trial ≤ 0:** Elastic, return (σ_trial, state_n) - -**If f_trial > 0:** Plastic, solve for Δλ: - -**Consistency condition:** -f(σ_{n+1}, ε_p_eq_{n+1}) = 0 - -**Discretized flow rule:** -ε_p_{n+1} = ε_p_n + Δλ n - -**Stress update:** -σ_{n+1} = σ_trial - 2G Δλ n - -where G = shear modulus. - -**Yield condition:** -q_{n+1} = q_trial - 3G Δλ = σ_y(ε_p_eq_n + Δλ) - -**Solve for Δλ (Newton):** -```julia -function return_mapping(σ_trial, state, material) - s_trial = dev(σ_trial) - q_trial = sqrt(1.5 * dcontract(s_trial, s_trial)) - n = s_trial / norm(s_trial) - - # Initial guess - Δλ = 0.0 - ε_p_eq_old = state.ε_p_eq - G = material.E / (2*(1 + material.ν)) - H = material.H - - for iter in 1:max_iter - # Current yield stress - σ_y = material.σ_y0 + H * (ε_p_eq_old + Δλ) - - # Residual - f = q_trial - 3*G*Δλ - σ_y - - if abs(f) < tol - break - end - - # Derivative - df_dΔλ = -3*G - H - - # Newton update - Δλ -= f / df_dΔλ - end - - # Update stress - σ = σ_trial - 2*G*Δλ * n - - # Update state - ε_p_new = state.ε_p + Δλ * n - ε_p_eq_new = ε_p_eq_old + Δλ - - return σ, PlasticState(ε_p_new, ε_p_eq_new, state.α) -end -``` - -### Consistent Tangent - -**For Newton convergence:** Need C_ep = dσ/dε (algorithmic tangent). - -**Elastic:** -C_ep = C_elastic - -**Plastic:** More complex! - -C_ep = C_elastic - (2G)² / (3G + H) * (n ⊗ n) - -where ⊗ = outer product. - -**Derivation:** Chain rule through return mapping. - -**Properties:** -- Symmetric (major symmetry) -- Positive-definite (for H > 0) -- Converges to C_elastic as Δλ → 0 - -## Internal State Storage - -**Per integration point:** -```julia -struct PlasticState{dim} - ε_p::SymmetricTensor{2,dim} # Plastic strain - ε_p_eq::Float64 # Equivalent plastic strain - α::SymmetricTensor{2,dim} # Backstress (kinematic) -end -``` - -**Element-level:** -```julia -struct PlasticElement - topology::AbstractTopology - basis::AbstractBasis - nodes::NTuple{N,Int} - state::Vector{PlasticState} # One per integration point! -end -``` - -**Key:** State is HISTORY-DEPENDENT, must be stored! - -## Nodal Assembly (Plasticity) - -```julia -function tangent_matvec_plastic!(w, v, u_current, material, elements, states) - Threads.@threads for node_i in 1:n_nodes - w_local = zero(Vec{3}) - - for elem in node_to_elements[node_i] - for (ip_idx, ip) in enumerate(integration_points(elem)) - # Current state at this integration point - state = states[elem][ip_idx] - - # Strain - ε = compute_strain(elem, ip, u_current) - - # Consistent tangent (elastic or plastic) - σ, state_new, C_ep = compute_stress_tangent(material, ε, state) - - for node_j in elem.nodes - # Tangent block - K_t_ij = compute_plastic_tangent_block(elem, node_i, node_j, C_ep, ip) - - v_j = Vec{3}(v[3*(node_j-1)+1:3*node_j]) - w_local += K_t_ij ⊡ v_j - end - - # Update state (for next iteration) - states[elem][ip_idx] = state_new - end - end - - w[3*(node_i-1)+1:3*node_i] = w_local - end -end -``` - -**Key:** State updated during tangent computation! - -## Next Steps - -1. Implement `J2Plasticity` material type -2. Implement `PlasticState` struct -3. Implement `return_mapping` algorithm -4. Implement `compute_stress_plastic` -5. Implement `consistent_tangent_plastic` -6. Implement hardening laws (isotropic, kinematic) -7. Implement state storage (per integration point) -8. Validate against analytical solutions -9. Validate against experimental data -10. Performance benchmarks - -"""