Merge pull request #74 from MultiSimOLab/fix-thermo

Fix thermo-mechanical coupling at constitutive model definition

Author: miguelmaso

src/PhysicalModels/ThermoElectroMechanicalModels.jl

@@ -117,7 +117,7 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro
   end
 
   function (obj::ThermoElectroMech_Bonet)(Λ::Float64=1.0; kwargs...)
-    @unpack Cv,θr, α, κ, γv, γd = obj.thermo
+    @unpack Cv, θr, α, κ, γv, γd = obj.thermo
     em = ElectroMechModel(obj.electro, obj.mechano)
     Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = em(;kwargs...)
     gd(θ) = 1/(γd+1) * ((θ/θr)^(γd+1) -1)
@@ -135,18 +135,18 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro
     ∂ηR∂F(F) = ∂ηR∂J(F)*H(F)
     ∂2ηR∂FF(F) = ×ᵢ⁴(∂ηR∂J(F) * F)
 
-    Ψ(F, E, θ, X...) = Ψem(F, E, X...)*(1.0+gd(θ))+gv(θ)*ηR(F)
+    Ψ(F, E, θ, X...) = Ψem(F, E, X...)*(1.0+gd(θ)) - θr*gv(θ)*ηR(F)
 
-    ∂Ψ_∂F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂F(F, E, X...) + gv(θ)*∂ηR∂F(F)
+    ∂Ψ_∂F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂F(F, E, X...) - θr*gv(θ)*∂ηR∂F(F)
     ∂Ψ_∂E(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂E(F, E, X...)
-    ∂Ψ_∂θ(F, E, θ, X...)  =  ∂gd(θ) *Ψem(F, E, X...) + ∂gv(θ)*ηR(F)
+    ∂Ψ_∂θ(F, E, θ, X...)  =  ∂gd(θ) *Ψem(F, E, X...) - θr*∂gv(θ)*ηR(F)
 
-    ∂2Ψ_∂2F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂FF(F, E, X...) + gv(θ)*∂2ηR∂FF(F)
+    ∂2Ψ_∂2F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂FF(F, E, X...) - θr*gv(θ)*∂2ηR∂FF(F)
     ∂2Ψ_∂2E(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂EE(F, E, X...)
-    ∂2Ψ_∂2θ(F, E, θ, X...)  =  ∂∂gd(θ) *Ψem(F, E, X...) + ∂∂gv(θ)*ηR(F)
+    ∂2Ψ_∂2θ(F, E, θ, X...)  =  ∂∂gd(θ) *Ψem(F, E, X...) - θr*∂∂gv(θ)*ηR(F)
 
     ∂ΨEF(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂EF(F, E, X...)
-    ∂ΨFθ(F, E, θ, X...) =  ∂gd(θ) *∂Ψem∂F(F, E, X...) + ∂gv(θ)*∂ηR∂F(F)
+    ∂ΨFθ(F, E, θ, X...) =  ∂gd(θ) *∂Ψem∂F(F, E, X...) - θr*∂gv(θ)*∂ηR∂F(F)
     ∂ΨEθ(F, E, θ, X...) =  ∂gd(θ) *∂Ψem∂E(F, E, X...)
 
     return (Ψ, ∂Ψ_∂F, ∂Ψ_∂E, ∂Ψ_∂θ, ∂2Ψ_∂2F, ∂2Ψ_∂2E, ∂2Ψ_∂2θ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ)

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -462,9 +462,11 @@ end
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
   ∇φ = VectorValue(1.0, 2.0, 3.0)
   θt = 3.4 - 1.0
+  θr = 293.0
+  Cv = 17.385
   modelMR = MooneyRivlin3D(λ=0.0, μ1=0.5, μ2=0.5)
   modelID = IdealDielectric(ε=1.0)
-  modelT = ThermalModel(Cv=17.385, θr=293.0, α=0.00156331, γv=2.0, γd=2.0)
+  modelT = ThermalModel(Cv=Cv, θr=θr, α=0.00156331, γv=2.0, γd=2.0)
 
   modelTEM = ThermoElectroMech_Bonet(modelT, modelID, modelMR)
   Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ = modelTEM()
@@ -497,6 +499,10 @@ end
   @test isapprox(∂ΨFθ(F(∇u), E(∇φ), θt), ∂2Ψ_∂2Fθ(F(∇u), E(∇φ), θt); rtol=1e-14)
   @test isapprox(∂ΨEF(F(∇u), E(∇φ), θt), ∂2Ψ_∂EF(F(∇u), E(∇φ), θt); rtol=1e-14)
 
+  F0 = I3
+  E0 = VectorValue(0.,0.,0.)
+  cv(F,E,θ,x...) = -θ*∂2Ψ∂2θ(F,E,θ,x...)
+  @test isapprox(Cv, cv(F0, E0, θr); rtol=1e-14)
 end