Thermal coupling with Bonet formulation

src/Exports.jl

@@ -44,6 +44,7 @@ end
 @publish PhysicalModels TransverseIsotropy3D
 @publish PhysicalModels TransverseIsotropy2D
 @publish PhysicalModels ThermalModel
+@publish PhysicalModels ThermalModel3rdLaw
 @publish PhysicalModels IdealDielectric
 @publish PhysicalModels Magnetic
 @publish PhysicalModels IdealMagnetic
@@ -57,7 +58,6 @@ end
 @publish PhysicalModels FlexoElectroModel
 @publish PhysicalModels ThermoElectroMech_Govindjee
 @publish PhysicalModels ThermoElectroMech_PINNs
-@publish PhysicalModels ThermoElectroMech_Bonet
 @publish PhysicalModels MagnetoMechModel
 @publish PhysicalModels ARAP2D
 @publish PhysicalModels ARAP2D_regularized

src/PhysicalModels/PhysicalModels.jl

@@ -45,14 +45,14 @@ export IdealMagnetic2D
 export HardMagnetic
 export HardMagnetic2D
 export ThermalModel
+export ThermalModel3rdLaw
 export ElectroMechModel
 export ThermoElectroMechModel
 export ThermoMechModel
 export ThermoMech_EntropicPolyconvex
 export FlexoElectroModel
 export ThermoElectroMech_Govindjee
 export ThermoElectroMech_PINNs
-export ThermoElectroMech_Bonet
 export MagnetoMechModel
 export GeneralizedMaxwell
 export ViscousIncompressible

src/PhysicalModels/ThermalModels.jl

@@ -20,5 +20,31 @@ struct ThermalModel <: Thermo
     ∂Ψθθ(δθ) = -obj.Cv / (δθ + obj.θr)
     return (Ψ, ∂Ψθ, ∂Ψθθ)
   end
+end
+
+
+struct ThermalModel3rdLaw <: Thermo
+  cv0::Float64
+  θr::Float64
+  α::Float64
+  κ::Float64
+  γv::Float64
+  γd::Float64
+  function ThermalModel3rdLaw(; cv0::Float64, θr::Float64, α::Float64, κ::Float64, γv::Float64, γd::Float64)
+    new(cv0, θr, α, κ, γv, γd)
+  end
+end
 
+function (obj::ThermalModel3rdLaw)()
+  @unpack cv0, θr, α, κ, γv, γd = obj
+  g(θ,θr,γ) = 1/(γ+1) * ((θ/θr)^(γ+1) -1)
+  ∂g(θ,θr,γ) = θ^γ / θr^(γ+1)
+  ∂∂g(θ,θr,γ) = γ*θ^(γ-1) / θr^(γ+1)
+  gd(θ) = g(θ,θr,γd)
+  ∂gd(θ) = ∂g(θ,θr,γd)
+  ∂∂gd(θ) = ∂∂g(θ,θr,γd)
+  gv(θ) = g(θ,θr,γv)
+  ∂gv(θ) = ∂g(θ,θr,γv)
+  ∂∂gv(θ) = ∂∂g(θ,θr,γv)
+  return (gv, ∂gv, ∂∂gv, gd, ∂gd, ∂∂gd)
 end

src/PhysicalModels/ThermoElectroMechanicalModels.jl

@@ -7,13 +7,19 @@ function update_state!(obj::TEM, state, F, E, θ, args...) where {TEM<:ThermoEle
   update_state!(obj.mechano, state, F, args...)
 end
 
+function update_time_step!(obj::TEM, Δt::Float64) where {TEM<:ThermoElectroMechano}
+  update_time_step!(obj.thermo,  Δt)
+  update_time_step!(obj.electro, Δt)
+  update_time_step!(obj.mechano, Δt)
+end
+
 struct ThermoElectroMechModel{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectroMechano
   thermo::T
   electro::E
   mechano::M
   fθ::Function
   dfdθ::Function
-  
+
   function ThermoElectroMechModel(thermo::T, electro::E, mechano::M; fθ::Function, dfdθ::Function) where {T<:Thermo,E<:Electro,M<:Mechano}
     new{T,E,M}(thermo, electro, mechano, fθ, dfdθ)
   end
@@ -22,6 +28,10 @@ struct ThermoElectroMechModel{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectroM
     new{T,E,M}(thermo, electro, mechano, fθ, dfdθ)
   end
 
+  function ThermoElectroMechModel(thermo::ThermalModel3rdLaw, electro::E, mechano::M) where {E<:Electro, M<:Mechano}
+    new{ThermalModel3rdLaw,E,M}(thermo, electro, mechano)
+  end
+
   function (obj::ThermoElectroMechModel)(Λ::Float64=1.0)
     Ψt, ∂Ψt_θ, ∂Ψt_θθ = obj.thermo(Λ)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
@@ -111,72 +121,33 @@ struct ThermoElectroMech_Govindjee{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoEle
 end
 
 
-struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectroMechano
-  thermo::T
-  electro::E
-  mechano::M
-
-  function ThermoElectroMech_Bonet(thermo::T, electro::E, mechano::M) where {T<:Thermo,E<:Electro,M<:Mechano}
-    new{T,E,M}(thermo, electro, mechano)
-  end
-
-  function ThermoElectroMech_Bonet(; thermo::T, electro::E, mechano::M) where {T<:Thermo,E<:Electro,M<:Mechano}
-    new{T,E,M}(thermo, electro, mechano)
-  end
-end
-
-g(θ,θr,γ) = 1/(γ+1) * ((θ/θr)^(γ+1) -1)
-∂g(θ,θr,γ) = θ^γ / θr^(γ+1)
-∂∂g(θ,θr,γ) = γ*θ^(γ-1) / θr^(γ+1)
-
-function (obj::ThermoElectroMech_Bonet)(Λ::Float64=1.0)
-  @unpack Cv, θr, α, κ, γv, γd = obj.thermo
+function (obj::ThermoElectroMechModel{ThermalModel3rdLaw,<:Electro,<:Mechano})(Λ::Float64=0.0)
+  @unpack cv0, θr, α, κ, γv, γd = obj.thermo
   em = ElectroMechModel(obj.electro, obj.mechano)
   Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = em()
-
-  gd(θ) = g(θ,θr,γd)
-  ∂gd(θ) = ∂g(θ,θr,γd)
-  ∂∂gd(θ) = ∂∂g(θ,θr,γd)
-  gv(θ) = g(θ,θr,γv)
-  ∂gv(θ) = ∂g(θ,θr,γv)
-  ∂∂gv(θ) = ∂∂g(θ,θr,γv)
-
-  J(F) = det(F)
-  H(F) = cof(F)
-
-  ηR(F) = α*(J(F) - 1.0)+Cv/γv
-  ∂ηR∂J(F) = α
-  ∂ηR∂F(F) = ∂ηR∂J(F)*H(F)
-  ∂2ηR∂FF(F) = ×ᵢ⁴(∂ηR∂J(F) * F)
+  gv, ∂gv, ∂∂gv, gd, ∂gd, ∂∂gd = obj.thermo()
+  ηR, ∂ηR∂F, ∂∂ηR∂FF = _getCoupling(obj.thermo, obj.mechano)
 
   Ψ(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...) - θ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...) - θr*∂gv(θ)*ηR(F)
 
-  ∂∂Ψ∂FF(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂FF(F, E, X...) - θr*gv(θ)*∂2ηR∂FF(F)
+  ∂∂Ψ∂FF(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂FF(F, E, X...) - θr*gv(θ)*∂∂ηR∂FF(F)
   ∂∂Ψ∂EE(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂EE(F, E, X...)
   ∂∂Ψ∂θθ(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...) - θr*∂gv(θ)*∂ηR∂F(F)
   ∂∂Ψ∂Eθ(F, E, θ, X...) =  ∂gd(θ) *∂Ψem∂E(F, E, X...)
-
   return (Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂θ, ∂∂Ψ∂FF, ∂∂Ψ∂EE, ∂∂Ψ∂θθ, ∂∂Ψ∂EF, ∂∂Ψ∂Fθ, ∂∂Ψ∂Eθ)
 end
 
-function update_time_step!(obj::ThermoElectroMech_Bonet, Δt::Float64)
-  update_time_step!(obj.thermo, Δt)
-  update_time_step!(obj.electro, Δt)
-  update_time_step!(obj.mechano, Δt)
-end
-
-function Dissipation(obj::ThermoElectroMech_Bonet)
-  @unpack Cv, θr, α, κ, γv, γd = obj.thermo
+function Dissipation(obj::ThermoElectroMechModel{ThermalModel3rdLaw,<:Electro,<:Mechano})
+  @unpack cv0, θr, α, κ, γv, γd = obj.thermo
+  gv, ∂gv, ∂∂gv, gd, ∂gd, ∂∂gd = obj.thermo()
   Dvis = Dissipation(obj.mechano)
-  gd(θ) = g(θ,θr,γd)
-  ∂gd(θ) = ∂g(θ,θr,γd)
   D(F, E, θ, X...) = (1 + gd(θ)) * Dvis(F, X...)
   ∂D∂θ(F, E, θ, X...) = ∂gd(θ) * Dvis(F, X...)
   return(D, ∂D∂θ)

src/PhysicalModels/ThermoMechanicalModels.jl

@@ -1,4 +1,21 @@
 
+# ===================
+# Common functions
+# ===================
+
+function initialize_state(obj::TM, points::Measure) where {TM<:ThermoMechano}
+  initialize_state(obj.mechano, points)
+end
+
+function update_state!(obj::TM, state, F, θ, args...) where {TM<:ThermoMechano}
+  update_state!(obj.mechano, state, F, args...)
+end
+
+function update_time_step!(obj::TM, Δt::Float64) where {TM<:ThermoMechano}
+  update_time_step!(obj.thermo,  Δt)
+  update_time_step!(obj.mechano, Δt)
+end
+
 # ===================
 # MultiPhysicalModel models
 # ===================
@@ -17,6 +34,10 @@ struct ThermoMechModel{T<:Thermo,M<:Mechano} <: ThermoMechano
     new{T,M}(thermo, mechano, fθ, dfdθ)
   end
 
+  function ThermoMechModel(thermo::ThermalModel3rdLaw, mechano::M) where {M<:Mechano}
+    new{ThermalModel3rdLaw,M}(thermo,mechano)
+  end
+
   function (obj::ThermoMechModel)(Λ::Float64=1.0)
     Ψt, ∂Ψt_θ, ∂Ψt_θθ = obj.thermo(Λ)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
@@ -35,6 +56,41 @@ struct ThermoMechModel{T<:Thermo,M<:Mechano} <: ThermoMechano
 end
 
 
+function (obj::ThermoMechModel{ThermalModel3rdLaw,<:Mechano})(Λ::Float64=0.0)
+  @unpack cv0, θr, α, κ, γv, γd = obj.thermo
+  gv, ∂gv, ∂∂gv, gd, ∂gd, ∂∂gd = obj.thermo()
+  Ψm, ∂Ψm∂F, ∂∂Ψm∂FF = obj.mechano()
+  ηR, ∂ηR∂F, ∂∂ηR∂FF = _getCoupling(obj.thermo, obj.mechano)
+  Ψ(F, θ, X...)       =  Ψm(F, X...)*(1.0+gd(θ)) - θr*gv(θ)*ηR(F)
+  ∂Ψ∂F(F, θ, X...)    =  (1.0+gd(θ)) *∂Ψm∂F(F, X...) - θr*gv(θ)*∂ηR∂F(F)
+  ∂Ψ∂θ(F, θ, X...)    =  ∂gd(θ) *Ψm(F, X...) - θr*∂gv(θ)*ηR(F)
+  ∂∂Ψ∂FF(F, θ, X...)  =  (1.0+gd(θ)) *∂∂Ψm∂FF(F, E, X...) - θr*gv(θ)*∂∂ηR∂FF(F)
+  ∂∂Ψ∂θθ(F, θ, X...)  =  ∂∂gd(θ) *Ψm(F, X...) - θr*∂∂gv(θ)*ηR(F)
+  ∂∂Ψ∂Fθ(F, θ, X...)  =  ∂gd(θ) *∂Ψm∂F(F, X...) - θr*∂gv(θ)*∂ηR∂F(F)
+  return (Ψ, ∂Ψ∂F, ∂Ψ∂θ, ∂∂Ψ∂FF, ∂∂Ψ∂θθ, ∂∂Ψ∂Fθ)
+end
+
+function _getCoupling(thermo::ThermalModel3rdLaw, ::Mechano)
+  @unpack cv0, θr, α, κ, γv, γd = thermo
+  J(F) = det(F)
+  H(F) = cof(F)
+  ηR(F) = α*(J(F) - 1.0) + cv0/γv
+  ∂ηR∂J(F) = α
+  ∂ηR∂F(F) = ∂ηR∂J(F)*H(F)
+  ∂∂ηR∂FF(F) = ×ᵢ⁴(∂ηR∂J(F) * F)
+  return (ηR, ∂ηR∂F, ∂∂ηR∂FF)
+end
+
+function Dissipation(obj::ThermoMechModel{ThermalModel3rdLaw,<:Mechano})
+  @unpack cv0, θr, α, κ, γv, γd = obj.thermo
+  gv, ∂gv, ∂∂gv, gd, ∂gd, ∂∂gd = obj.thermo()
+  Dvis = Dissipation(obj.mechano)
+  D(F, θ, X...) = (1 + gd(θ)) * Dvis(F, X...)
+  ∂D∂θ(F, θ, X...) = ∂gd(θ) * Dvis(F, X...)
+  return(D, ∂D∂θ)
+end
+
+
 struct ThermoMech_EntropicPolyconvex{T<:Thermo,M<:Mechano} <: ThermoMechano
   thermo::T
   mechano::M

src/PhysicalModels/ViscousModels.jl

@@ -165,9 +165,18 @@ function return_mapping_algorithm!(obj::ViscousIncompressible,
   res, ∂res = JacobianReturnMapping(γα, Ce, Se(Ce), Se_trial, ∂Se∂Ce(Ce), C, λα)
   maxiter = 20
   tol = 1e-6
-  for _ in 1:maxiter
+  for i in 1:maxiter
     #---------- Update -----------#
-    Δu = -∂res \ res[:]
+    local Δu
+    try
+      Δu = -∂res \ res[:]
+    catch e
+      if e isa LinearAlgebra.SingularException
+        error("Singular jacobian in return mapping algorithm (singular value at pos $(e.info), iteration $i)")
+      else
+        rethrow()
+      end
+    end
     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 ---------#

src/WeakForms/WeakForms.jl

@@ -412,7 +412,7 @@ function jacobian(physicalmodel::ThermoElectroMech_PINNs, kine::NTuple{2,Kinemat
   jacobian(physicalmodel, ThermoElectro, kine,(u, φ, θ), dθ, vφ, dΩ, Λ)
 end
 
-function transient_residual(physicalmodel::ThermoElectroMech_Bonet, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0)
+function transient_residual(physicalmodel::ThermoElectroMechModel, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0)
   κ = physicalmodel.thermo.κ
   return ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ
 end

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -5,7 +5,7 @@ using StaticArrays
 using Test
 using HyperFEM.PhysicalModels
 using HyperFEM.TensorAlgebra
-
+using HyperFEM.IO
 
 
 import Base: -
@@ -470,12 +470,11 @@ end
   ∇φ = VectorValue(1.0, 2.0, 3.0)
   θt = 3.4 - 1.0
   θr = 293.0
-  Cv = 17.385
+  cv0 = 17.385
   modelMR = MooneyRivlin3D(λ=0.0, μ1=0.5, μ2=0.5)
   modelID = IdealDielectric(ε=1.0)
-  modelT = ThermalModel(Cv=Cv, θr=θr, α=0.00156331, γv=2.0, γd=2.0)
-
-  modelTEM = ThermoElectroMech_Bonet(modelT, modelID, modelMR)
+  modelT = ThermalModel3rdLaw(cv0=cv0, θr=θr, α=0.00156331, κ=1.0, γv=2.0, γd=2.0)
+  modelTEM = ThermoElectroMechModel(modelT, modelID, modelMR)
   Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ = modelTEM()
 
   K = Kinematics(Mechano, Solid)
@@ -509,7 +508,8 @@ end
   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)
+  @test isapprox(cv0, cv(F0, E0, θr); rtol=1e-14)
+  # @test isapprox(0, Ψ(F0, E0, 0); atol=1e-14)
 end