From 99784c40b94eff44961afa1a139df958ae1e40ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 30 Sep 2025 17:54:42 +0200 Subject: [PATCH 1/2] fixes in const models --- src/PhysicalModels/KinematicModels.jl | 2 +- src/PhysicalModels/PhysicalModels.jl | 16 ++++++----- src/PhysicalModels/ViscousModels.jl | 38 ++++++++++++++++++++++----- 3 files changed, 42 insertions(+), 14 deletions(-) diff --git a/src/PhysicalModels/KinematicModels.jl b/src/PhysicalModels/KinematicModels.jl index 2b1cd74..b89a783 100644 --- a/src/PhysicalModels/KinematicModels.jl +++ b/src/PhysicalModels/KinematicModels.jl @@ -18,7 +18,7 @@ struct Kinematics{A,B} <: KinematicModel function Kinematics(::Type{Visco}; F::Function=(∇u) -> one(∇u) + ∇u) C(F) = F' * F - Ce(C,Uvα⁻¹) = Uvα⁻¹ * C * Uvα⁻¹ + Ce(C,Uvα⁻¹) = Uvα⁻¹' * C * Uvα⁻¹ metrics = (F, C, Ce) A = typeof(metrics) new{A,Visco}(metrics) diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index b9e5284..62dbc67 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/src/PhysicalModels/PhysicalModels.jl @@ -443,8 +443,10 @@ end struct ComposedElasticModel <: Elasto Model1::Elasto Model2::Elasto + Kinematic function ComposedElasticModel(model1::Elasto,model2::Elasto) - new(model1,model2) + @assert model1.Kinematic === model2.Kinematic + new(model1,model2,model1.Kinematic) end function (obj::ComposedElasticModel)(Λ::Float64=1.0) DΨ1 = obj.Model1(Λ) @@ -547,8 +549,8 @@ struct VolumetricEnergy{A} <: Elasto Ψ(F) = (λ / 2.0) * (J(F) - 1)^2 ∂Ψ_∂J(F) = λ * (J(F) - 1) ∂Ψ2_∂J2(F) = λ - ∂Ψu(F) = ∂Ψ_∂J(F) * H(F) - ∂Ψuu(F) =∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F) + ∂Ψu(F) = ∂Ψ_∂J(F) * H(F) + ∂Ψuu(F) = ∂Ψ2_∂J2(F) * (H(F) ⊗ H(F)) + ×ᵢ⁴(∂Ψ_∂J(F) * F) return (Ψ, ∂Ψu, ∂Ψuu) end end @@ -1102,16 +1104,16 @@ struct IncompressibleNeoHookean3D{A} <: Elasto function (obj::IncompressibleNeoHookean3D)(::KinematicDescription{:SecondPiola}, Λ::Float64=1.0) Ψ(C) = obj.μ / 2 * tr(C) * det(C)^(-1/3) S(C) = begin - J = det(C) + detC = det(C) invC = inv(C) - obj.μ * J^(-1/3) * I3 - obj.μ / 3 * tr(C) * J^(-1/3) * invC + obj.μ * detC^(-1/3) * I3 - obj.μ / 3 * tr(C) * detC^(-1/3) * invC end ∂S∂C(C) = begin - J = det(C) + detC = det(C) trC = tr(C) invC = inv(C) IinvC = I3 ⊗ invC - 1/3 * obj.μ * J^(-1/3) * (4/3*trC*invC⊗invC -(IinvC+IinvC') -trC/J*×ᵢ⁴(C)) + 1/3 * obj.μ * detC^(-1/3) * (4/3*trC*invC⊗invC -(IinvC+IinvC') -trC/detC*×ᵢ⁴(C)) end return (Ψ, S, ∂S∂C) end diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl index 93682fc..467d1c7 100644 --- a/src/PhysicalModels/ViscousModels.jl +++ b/src/PhysicalModels/ViscousModels.jl @@ -15,9 +15,10 @@ struct ViscousIncompressible{T} <: Visco end function (obj::ViscousIncompressible)(Λ::Float64=1.0; Δt::Float64) Ψe, Se, ∂Se∂Ce = obj.ShortTerm(KinematicDescription{:SecondPiola}()) + Ψ(F, Fn, state) = Energy(obj, Δt, Ψe, Se, ∂Se∂Ce, F, Fn, state) ∂Ψ∂F(F, Fn, state) = Piola(obj, Δt, Se, ∂Se∂Ce, F, Fn, state) ∂Ψ∂F∂F(F, Fn, state) = Tangent(obj, Δt, Se, ∂Se∂Ce, F, Fn, state) - return Ψe, ∂Ψ∂F, ∂Ψ∂F∂F + return Ψ, ∂Ψ∂F, ∂Ψ∂F∂F end end @@ -117,7 +118,6 @@ function return_mapping_algorithm!(obj::ViscousIncompressible, Δt::Float64, for _ in 1:maxiter #---------- Update -----------# Δu = -∂res \ res[:] - # Ce += reshape(Δu[1:end-1], 3, 3) Ce += TensorValue{3,3}(Tuple(Δu[1:end-1])) # TODO: Check reconstruction of TensorValue. ERROR: MethodError: no method matching (TensorValue{3, 3})(::Vector{Float64}) λα += Δu[end] #---- Residual and jacobian ---------# @@ -156,7 +156,7 @@ function JacobianReturnMapping(γα, Ce, Se, Se_trial, ∂Se∂Ce, F, λα) ∂res1_∂Ce = ∂Se∂Ce - (1-γα) * λα * ×ᵢ⁴(Ce) ∂res1_∂λα = -(1-γα) * Ge ∂res2_∂Ce = Ge - res = [get_array(res1)[:]; res2] #TODO: Check the TensorValue interface + res = [get_array(res1)[:]; res2] ∂res = MMatrix{10,10}(zeros(10, 10)) ∂res[1:9, 1:9] = get_array(∂res1_∂Ce) ∂res[1:9, 10] = get_array(∂res1_∂λα)[:] @@ -180,7 +180,7 @@ Viscous 1st Piola-Kirchhoff stress - `Pα::SMatrix` """ function ViscousPiola(Se::Function, Ce::TensorValue, invUv::TensorValue, F::TensorValue) - Sα = invUv * Se(Ce) * invUv + Sα = invUv' * Se(Ce) * invUv F * Sα end @@ -232,7 +232,7 @@ end """ -Tangent operator of Ce for at fixed Uv +Tangent operator of Ce at fixed Uv """ function ∂Ce_∂C_Uvfixed(invUv) invUv ⊗₁₃²⁴ invUv @@ -298,7 +298,7 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64, # Sv:(D(δC_{Uvfixed})[ΔC]) #------------------------------------------ Sv = invUv_Se * invUv - C3 = Sv ⊗₁₃²⁴ I3 + C3 = I3 ⊗₁₃²⁴ Sv #------------------------------------------ # Total Contribution #------------------------------------------ @@ -307,6 +307,32 @@ function ViscousTangentOperator(obj::ViscousIncompressible, Δt::Float64, end +function Energy(obj::ViscousIncompressible, Δt::Float64, + Ψe::Function, Se_::Function, ∂Se∂Ce_::Function, + F::TensorValue, Fn::TensorValue, stateVars::VectorValue) + state_vars = get_array(stateVars) + Uvn = TensorValue{3,3}(Tuple(state_vars[1:9])) # TODO: Update tensor slicing until next Gridap version has been released + λαn = state_vars[10] + #------------------------------------------ + # Get kinematics + #------------------------------------------ + invUvn = inv(Uvn) + _, C_, Ce_ = get_Kinematics(obj.Kinematic) + C = C_(F) + Cn = C_(Fn) + Ceᵗʳ = Ce_(C, invUvn) + Cen = Ce_(Cn, invUvn) + #------------------------------------------ + # Return mapping algorithm + #------------------------------------------ + Ce, _ = return_mapping_algorithm!(obj, Δt, Se_, ∂Se∂Ce_, F, Ceᵗʳ, Cen, λαn) + #------------------------------------------ + # Elastic energy + #------------------------------------------ + Ψe(Ce) +end + + """ First Piola-Kirchhoff for the incompressible case From 9c11246358d9f9caa7ef108a1ee61dacab6b7280 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3?= Date: Tue, 30 Sep 2025 17:54:57 +0200 Subject: [PATCH 2/2] added viscous tests --- .../ViscousModelsTests.jl | 167 ++++++++++++++++++ test/TestConstitutiveModels/runtests.jl | 6 + 2 files changed, 173 insertions(+) create mode 100644 test/TestConstitutiveModels/ViscousModelsTests.jl diff --git a/test/TestConstitutiveModels/ViscousModelsTests.jl b/test/TestConstitutiveModels/ViscousModelsTests.jl new file mode 100644 index 0000000..cb390ed --- /dev/null +++ b/test/TestConstitutiveModels/ViscousModelsTests.jl @@ -0,0 +1,167 @@ + +using Gridap.TensorValues +using Gridap.Arrays +using HyperFEM.TensorAlgebra +using HyperFEM.PhysicalModels +using StaticArrays +using Test +import Base:≈ +≈(a::MultiValue, b::MultiValue; kwargs...) = ≈(get_array(a), get_array(b); kwargs...) + + +μ = 1.367e4 # Pa +N = 7.860e5 # - +λ = μ*100 # Pa +μ1 = 3.153e5 # Pa +τ1 = 10.72 # s +μ2 = 5.639e5 # Pa +τ2 = 0.82 # s +μ3 = 1.981e5 # Pa +τ3 = 498.8 # s + + +function isochoric_F(F) + J = det(F) + @assert J > 0 "Non-physical deformation, got det(F) < 0 (det(F) = $J)" + J^(-1/3) * F +end + + +function numerical_piola(Ψ, F, ϵ=1e-6) + P = MMatrix{3,3}(zeros(Float64,9)) + for i in 1:9 + Fp = mutable(F) + Fm = mutable(F) + Fp[i] += ϵ + Fm[i] -= ϵ + Ψp = Ψ(TensorValue(Fp)) + Ψm = Ψ(TensorValue(Fm)) + P[i] = (Ψp - Ψm) / 2ϵ + end + return TensorValue(P) +end + + +function numerical_tangent(Ψ, F, ϵ=1e-5) + H = MMatrix{9,9}(zeros(Float64,81)) + for j in 1:9 + ej = zeros(9); ej[j] = ϵ + for i in 1:9 + ei = zeros(9); ei[i] = ϵ + if i == j + Ψp = Ψ(F + TensorValue(ei...)) + Ψ0 = Ψ(F) + Ψm = Ψ(F - TensorValue(ei...)) + H[i,j] = (Ψp - 2Ψ0 + Ψm) / ϵ^2 + else + Ψpp = Ψ(F + TensorValue(( ei+ej)...)) + Ψpm = Ψ(F + TensorValue(( ei-ej)...)) + Ψmp = Ψ(F + TensorValue((-ei+ej)...)) + Ψmm = Ψ(F + TensorValue((-ei-ej)...)) + H[i,j] = (Ψpp - Ψpm - Ψmp + Ψmm) / 4ϵ^2 + end + end + end + return TensorValue(H) +end + + +function richardson_expansion(func, x, ϵ) + (4.0func(x,ϵ) - func(x,2ϵ)) / 3.0 +end + + +function test_viscous_derivatives_numerical(model; rtolP=1e-12, rtolH=1e-12) + Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2) + F = TensorValue(1.:9...) * 1e-3 + I3 + Fn = TensorValue(1.:9...) * 5e-4 + I3 + Uvn = isochoric_F(TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3) + λvn = 1e-3 + Avn = VectorValue(Uvn.data..., λvn) + piola = richardson_expansion((F, ϵ) -> numerical_piola(F -> Ψ(F, Fn, Avn), F, ϵ), F, 1e-5) + tangent = richardson_expansion((F, ϵ) -> numerical_tangent(F -> Ψ(F, Fn, Avn), F, ϵ), F, 1e-4) + @test isapprox(∂Ψu(F, Fn, Avn), piola, rtol=rtolP) + @test isapprox(∂Ψuu(F, Fn, Avn), tangent, rtol=rtolH) +end + + +function test_elastic_derivatives_numerical(model; rtolP=1e-12, rtolH=1e-12) + Ψ, ∂Ψu, ∂Ψuu = model() + F = TensorValue(1.:9...) * 1e-3 + I3 + piola = richardson_expansion((F,ϵ) -> numerical_piola(Ψ,F,ϵ), F, 1e-5) + tangent = richardson_expansion((F,ϵ) -> numerical_tangent(Ψ,F,ϵ), F, 1e-4) + @test isapprox(∂Ψu(F), piola, rtol=rtolP) + @test isapprox(∂Ψuu(F), tangent, rtol=rtolH) +end + + +struct EmptyElastic <: Elasto + Kinematic::KinematicModel + function EmptyElastic() + new(Kinematics(Elasto)) + end + function (::EmptyElastic)(Λ=0.0) + Ψ(F) = 0.0 + ∂Ψu(F) = TensorValue(zeros(9)...) + ∂Ψuu(F) = TensorValue(zeros(81)...) + return (Ψ, ∂Ψu, ∂Ψuu) + end +end + + +@testset "VolumetricEnergy" begin + hyper_elastic_model = VolumetricEnergy(λ=λ) + test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-10, rtolH=1e-9) +end; + +@testset "EightChain" begin + hyper_elastic_model = EightChain(μ=μ, N=N) + test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-3, rtolH=1e-2) +end; + +@testset "EightChain+VolumetricEnergy" begin + hyper_elastic_model = EightChain(μ=μ, N=N) + VolumetricEnergy(λ=λ) + test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-5, rtolH=1e-4) +end; + +@testset "NeoHookean3D" begin + hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ) + test_elastic_derivatives_numerical(hyper_elastic_model, rtolP=1e-10, rtolH=1e-9) +end; + +@testset "ViscousIncompressible" begin + visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1) + test_viscous_derivatives_numerical(visco, rtolP=1e-3, rtolH=1e-3) +end + +@testset "ViscousIncompressible2" begin + visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1.0), 10.0) + Ψ, ∂Ψu, ∂Ψuu = visco(Δt = 0.1) + F = 1e-2*TensorValue(1,2,3,4,5,6,7,8,9) + I3 + Fn = 5e-3*TensorValue(1,2,3,4,5,6,7,8,9) + I3 + Fvn = 2e-2*TensorValue(1.0,2.0,3.0,4.0,5.0,8.7,6.5,4.3,6.5) + I3 + Cvn = Fvn'*Fvn + Uvn = sqrt(Cvn) + Avn = VectorValue(Uvn.data...,0.0) + @test isapprox(norm(∂Ψu(F, Fn, Avn)), 0.20303772905627682, rtol=1e-10) + @test isapprox(norm(∂Ψuu(F, Fn, Avn)), 4.847586088299776, rtol=1e-10) +end + +@testset "GeneralizedMaxwell EightChain 0-branch" begin + hyper_elastic_model = EightChain(μ=μ, N=N) + VolumetricEnergy(λ=λ) + cons_model = GeneralizedMaxwell(hyper_elastic_model) + test_viscous_derivatives_numerical(cons_model, rtolP=1e-5, rtolH=1e-4) +end; + +@testset "GeneralizedMaxwell NeoHookean 0-branch" begin + hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ) + cons_model = GeneralizedMaxwell(hyper_elastic_model) + test_viscous_derivatives_numerical(cons_model, rtolP=1e-10, rtolH=1e-9) +end; + +@testset "GeneralizedMaxwell NeoHookean 1-branch" begin + hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ) + branch1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1) + cons_model = GeneralizedMaxwell(hyper_elastic_model, branch1) + test_viscous_derivatives_numerical(cons_model, rtolP=1e-3, rtolH=1e-2) +end; diff --git a/test/TestConstitutiveModels/runtests.jl b/test/TestConstitutiveModels/runtests.jl index 95e7bd6..ff18c32 100644 --- a/test/TestConstitutiveModels/runtests.jl +++ b/test/TestConstitutiveModels/runtests.jl @@ -1,7 +1,13 @@ +using HyperFEM + @testset "ConstitutiveModels" verbose = true begin @time begin include("PhysicalModels.jl") end + @time begin + include("ViscousModelsTests.jl") + end + end \ No newline at end of file