Merge pull request #73 from MultiSimOLab/thermo-visco

Thermo-electro-mechanical constitutive model

Author: miguelmaso

src/ComputationalModels/CartesianTags.jl

@@ -2,45 +2,83 @@
 "Shortcuts for the tags of cartesian discrete models."
 module CartesianTags
 
-  "Tags indicating points, edges and faces at X0."
-  const faceX0 = [1,3,5,7,13,15,17,19,25]
+"Tags indicating points, edges and faces at plane X0."
+const faceX0 = [1, 3, 5, 7, 13, 15, 17, 19, 25]
 
-  "Tags indicating points, edges and faces at X1."
-  const faceX1 = [2,4,6,8,14,16,18,20,26]
+"Tags indicating points, edges and faces at plane X1."
+const faceX1 = [2, 4, 6, 8, 14, 16, 18, 20, 26]
 
-  "Tags indicating points, edges and faces at Y0."
-  const faceY0 = [1,2,5,6,9,11,17,18,23]
+"Tags indicating points, edges and faces at plane Y0."
+const faceY0 = [1, 2, 5, 6, 9, 11, 17, 18, 23]
 
-  "Tags indicating points, edges and faces at Y1."
-  const faceY1 = [3,4,7,8,10,12,19,20,24]
+"Tags indicating points, edges and faces at plane Y1."
+const faceY1 = [3, 4, 7, 8, 10, 12, 19, 20, 24]
 
-  "Tags indicating points, edges and faces at Z0."
-  const faceZ0 = [1,2,3,4,9,10,13,14,21]
+"Tags indicating points, edges and faces at plane Z0."
+const faceZ0 = [1, 2, 3, 4, 9, 10, 13, 14, 21]
 
-  "Tags indicating points, edges and faces at Z1."
-  const faceZ1 = [5,6,7,8,11,12,15,16,22]
+"Tags indicating points, edges and faces at plane Z1."
+const faceZ1 = [5, 6, 7, 8, 11, 12, 15, 16, 22]
 
-  "Tag indicating the point at corner X0, Y0, Z0."
-  const corner000 = [1]
 
-  "Tag indicating the point at corner X1, Y0, Z0."
-  const corner100 = [2]
+"Tag indicating the edge at  X, Y0, Z0."
+const edgeX00 = [9]
 
-  "Tag indicating the point at corner X0, Y1, Z0."
-  const corner010 = [3]
+"Tag indicating the edge at  X, Y1, Z0."
+const edgeX10 = [10]
 
-  "Tag indicating the point at corner X1, Y1, Z0."
-  const corner110 = [4]
+"Tag indicating the edge at  X, Y0, Z1."
+const edgeX01 = [11]
 
-  "Tag indicating the point at corner X0, Y0, Z1."
-  const corner001 = [5]
+"Tag indicating the edge at  X, Y1, Z1."
+const edgeX11 = [12]
 
-  "Tag indicating the point at corner X1, Y0, Z1."
-  const corner101 = [6]
+"Tag indicating the edge at  X0, Y, Z0."
+const edge0Y0 = [13]
 
-  "Tag indicating the point at corner X0, Y1, Z1."
-  const corner011 = [7]
+"Tag indicating the edge at  X1, Y, Z0."
+const edge1Y0 = [14]
 
-  "Tag indicating the point at corner X1, Y1, Z1."
-  const corner111 = [8]
+"Tag indicating the edge at  X0, Y, Z1."
+const edge0Y1 = [15]
+
+"Tag indicating the edge at  X1, Y, Z1."
+const edge1Y1 = [16]
+
+"Tag indicating the edge at  X0, Y0, Z."
+const edge00Z = [17]
+
+"Tag indicating the edge at  X1, Y0, Z."
+const edge10Z = [18]
+
+"Tag indicating the edge at  X0, Y1, Z."
+const edge01Z = [19]
+
+"Tag indicating the edge at  X1, Y1, Z."
+const edge11Z = [20]
+
+
+"Tag indicating the point at corner X0, Y0, Z0."
+const corner000 = [1]
+
+"Tag indicating the point at corner X1, Y0, Z0."
+const corner100 = [2]
+
+"Tag indicating the point at corner X0, Y1, Z0."
+const corner010 = [3]
+
+"Tag indicating the point at corner X1, Y1, Z0."
+const corner110 = [4]
+
+"Tag indicating the point at corner X0, Y0, Z1."
+const corner001 = [5]
+
+"Tag indicating the point at corner X1, Y0, Z1."
+const corner101 = [6]
+
+"Tag indicating the point at corner X0, Y1, Z1."
+const corner011 = [7]
+
+"Tag indicating the point at corner X1, Y1, Z1."
+const corner111 = [8]
 end

src/Exports.jl

@@ -85,7 +85,8 @@ end
 @publish PhysicalModels EvolutiveKinematics
 @publish PhysicalModels get_Kinematics
 @publish PhysicalModels getIsoInvariants
-  
+
+@publish PhysicalModels Dissipation
 @publish PhysicalModels initializeStateVariables
 @publish PhysicalModels updateStateVariables!
 

src/PhysicalModels/ThermoElectroMechanicalModels.jl

@@ -107,7 +107,7 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro
   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
@@ -116,41 +116,40 @@ struct ThermoElectroMech_Bonet{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoElectro
     new{T,E,M}(thermo, electro, mechano)
   end
 
-  function (obj::ThermoElectroMech_Bonet)(Λ::Float64=1.0)
+  function (obj::ThermoElectroMech_Bonet)(Λ::Float64=1.0; kwargs...)
     @unpack Cv,θr, α, κ, γv, γd = obj.thermo
-    Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE = _getCoupling(obj.electro, obj.mechano, Λ)
-    gd(δθ)=1/(γd+1) * (((δθ+θr)/θr)^(γd+1) -1)
-    ∂gd(δθ)=(δθ+θr)^γd / θr^(γd+1)
-    ∂∂gd(δθ)= γd*(δθ+θr)^(γd-1) / θr^(γd+1)
-    gv(δθ)=1/(γv+1) * (((δθ+θr)/θr)^(γv+1) -1)
-    ∂gv(δθ)=(δθ+θr)^γv / θr^(γv+1)
-    ∂∂gv(δθ)=γv*(δθ+θr)^(γv-1) / θr^(γv+1)
+    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)
+    ∂gd(θ) = θ^γd / θr^(γd+1)
+    ∂∂gd(θ) = γd*θ^(γd-1) / θr^(γd+1)
+    gv(θ) = 1/(γv+1) * ((θ/θr)^(γv+1) -1)
+    ∂gv(θ) = θ^γv / θr^(γv+1)
+    ∂∂gv(θ) = γv*θ^(γv-1) / θr^(γv+1)
 
     J(F) = det(F)
     H(F) = det(F) * inv(F)'
 
-    η(F) =α*(J(F) - 1.0)+Cv/γv
-    ∂η∂J(F)=α
-    ∂η∂F(F)=∂η∂J(F)*H(F)
-    ∂2η∂FF(F)=×ᵢ⁴(∂η∂J(F) * F)
-
-    Ψ(F,E,δθ) = Ψem(F,E)*(1.0+gd(δθ))+gv(δθ)*η(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)
 
-    ∂Ψ_∂F(F, E, δθ)  =   (1.0+gd(δθ)) *∂Ψem∂F(F, E) + gv(δθ)*∂η∂F(F)
-    ∂Ψ_∂E(F, E, δθ)  =   (1.0+gd(δθ)) *∂Ψem∂E(F, E)
-    ∂Ψ_∂δθ(F, E, δθ) =   ∂gd(δθ) *Ψem(F, E) + ∂gv(δθ)*η(F)
+    Ψ(F, E, θ, X...) = Ψem(F, E, X...)*(1.0+gd(θ))+gv(θ)*ηR(F)
 
-    ∂2Ψ_∂2F(F, E, δθ) =  (1.0+gd(δθ)) *∂Ψem∂FF(F, E) + gv(δθ)*∂2η∂FF(F)
-    ∂2Ψ_∂2E(F, E, δθ) =  (1.0+gd(δθ)) *∂Ψem∂EE(F, E)
-    ∂2Ψ_∂2δθ(F, E, δθ) =  ∂∂gd(δθ) *Ψem(F, E) + ∂∂gv(δθ)*η(F)
+    ∂Ψ_∂F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂F(F, E, X...) + 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)
 
-    ∂ΨEF(F, E, δθ) =  (1.0+gd(δθ)) *∂Ψem∂EF(F, E)
-    ∂ΨFδθ(F, E, δθ) =  ∂gd(δθ) *∂Ψem∂F(F, E) + ∂gv(δθ)*∂η∂F(F)
-    ∂ΨEδθ(F, E, δθ) =  ∂gd(δθ) *∂Ψem∂E(F, E)
+    ∂2Ψ_∂2F(F, E, θ, X...)  =  (1.0+gd(θ)) *∂Ψem∂FF(F, E, X...) + 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)
 
-    entropy(F, E, δθ) = -∂Ψ_∂δθ(F, E, δθ)
+    ∂Ψ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)
+    ∂ΨEθ(F, E, θ, X...) =  ∂gd(θ) *∂Ψem∂E(F, E, X...)
 
-    return (Ψ, ∂Ψ_∂F, ∂Ψ_∂E, ∂Ψ_∂δθ, ∂2Ψ_∂2F, ∂2Ψ_∂2E, ∂2Ψ_∂2δθ, ∂ΨEF, ∂ΨFδθ, ∂ΨEδθ, entropy)
+    return (Ψ, ∂Ψ_∂F, ∂Ψ_∂E, ∂Ψ_∂θ, ∂2Ψ_∂2F, ∂2Ψ_∂2E, ∂2Ψ_∂2θ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ)
   end
 end
 
@@ -159,7 +158,7 @@ function Dissipation(obj::ThermoElectroMech_Bonet, Δt)
   Dvis = Dissipation(obj.mechano, Δt)
   gd(δθ) = 1/(γd+1) * (((δθ+θr)/θr)^(γd+1) -1)
   ∂gd(δθ) = (δθ+θr)^γd / θr^(γd+1)
-  D(F, E, δθ, A...) = (1 + gd(δθ)) * Dvis(F, A...)
-  ∂D∂θ(F, E, δθ, A...) = ∂gd(δθ) * Dvis(F, A...)
+  D(F, E, δθ, X...) = (1 + gd(δθ)) * Dvis(F, X...)
+  ∂D∂θ(F, E, δθ, X...) = ∂gd(δθ) * Dvis(F, X...)
   return(D, ∂D∂θ)
 end

src/WeakForms/WeakForms.jl

@@ -85,9 +85,16 @@ function residual(physicalmodel::ThermoElectroMechano, ::Type{Electro}, kine::NT
     return -1.0*∫(∇(vφ)' ⋅ (∂Ψφ ∘ (F∘(∇(u)'), E∘(∇(φ)), θ)))dΩ
 end
 
-function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0)
-    κ=physicalmodel.thermo.κ
-    return ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ
+function residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0, Δt=0.0, vars...)
+    κ = physicalmodel.thermo.κ
+    D, ∂D = Dissipation(Δt)
+    return ∫(κ * ∇(θ) ⋅ ∇(vθ) -D(u, φ, θ, vars...))dΩ
+end
+
+function transient_residual(physicalmodel::ThermoElectroMechano, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), (un, φn, θn), vθ, dΩ, Λ=1.0, Δt=0.0, vars...)
+    DΨ = physicalmodel(Δt=Δt)
+    η = -DΨ[4]
+    return ∫((1/Δt)*(θ*η(F,E,θ,vars...) - θn*η(Fn,En,θn,vars...) + η(F,E,θ)*(θ - θn)))dΩ
 end
 
 
@@ -410,4 +417,13 @@ function jacobian(physicalmodel::ThermoElectroMech_PINNs, kine::NTuple{2,Kinemat
 end
 
 
+
+
+function transient_residual(physicalmodel::ThermoElectroMech_Bonet, ::Type{Thermo}, kine::NTuple{3,KinematicModel}, (u, φ, θ), vθ, dΩ, Λ=1.0)
+    κ = physicalmodel.thermo.κ
+
+    return ∫(κ * ∇(θ) ⋅ ∇(vθ))dΩ
+end
+
+
 end

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -467,7 +467,7 @@ end
   modelT = ThermalModel(Cv=17.385, θr=293.0, α=0.00156331, γv=2.0, γd=2.0)
 
   modelTEM = ThermoElectroMech_Bonet(modelT, modelID, modelMR)
-  Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ, η = modelTEM()
+  Ψ, ∂Ψu, ∂ΨE, ∂Ψθ, ∂ΨFF, ∂ΨEE, ∂2Ψθθ, ∂ΨEF, ∂ΨFθ, ∂ΨEθ = modelTEM()
 
   K = Kinematics(Mechano, Solid)
   F, _, _ = get_Kinematics(K)