Skip to content

Commit ae07cff

Browse files
jmartfrutjmartfrut
committed
Adds thermo-electro-visco simulation
Implements a thermo-electro-visco simulation using HyperFEM and Gridap. It sets up the computational model, defines material properties, boundary conditions, and weak forms for thermo-electro-mechanical coupling. The simulation solves a staggered system of equations for electric potential, displacement, and temperature, and includes a post-processor to save simulation results.
1 parent a23b75d commit ae07cff

File tree

2 files changed

+251
-28
lines changed

2 files changed

+251
-28
lines changed

‎src/ComputationalModels/CartesianTags.jl‎

Lines changed: 65 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -2,45 +2,82 @@
22
"Shortcuts for the tags of cartesian discrete models."
33
module CartesianTags
44

5-
"Tags indicating points, edges and faces at X0."
6-
const faceX0 = [1,3,5,7,13,15,17,19,25]
5+
"Tags indicating points, edges and faces at X0."
6+
const faceX0 = [1, 3, 5, 7, 13, 15, 17, 19, 25]
77

8-
"Tags indicating points, edges and faces at X1."
9-
const faceX1 = [2,4,6,8,14,16,18,20,26]
8+
"Tags indicating points, edges and faces at X1."
9+
const faceX1 = [2, 4, 6, 8, 14, 16, 18, 20, 26]
1010

11-
"Tags indicating points, edges and faces at Y0."
12-
const faceY0 = [1,2,5,6,9,11,17,18,23]
11+
"Tags indicating points, edges and faces at Y0."
12+
const faceY0 = [1, 2, 5, 6, 9, 11, 17, 18, 23]
1313

14-
"Tags indicating points, edges and faces at Y1."
15-
const faceY1 = [3,4,7,8,10,12,19,20,24]
14+
"Tags indicating points, edges and faces at Y1."
15+
const faceY1 = [3, 4, 7, 8, 10, 12, 19, 20, 24]
1616

17-
"Tags indicating points, edges and faces at Z0."
18-
const faceZ0 = [1,2,3,4,9,10,13,14,21]
17+
"Tags indicating points, edges and faces at Z0."
18+
const faceZ0 = [1, 2, 3, 4, 9, 10, 13, 14, 21]
1919

20-
"Tags indicating points, edges and faces at Z1."
21-
const faceZ1 = [5,6,7,8,11,12,15,16,22]
20+
"Tags indicating points, edges and faces at Z1."
21+
const faceZ1 = [5, 6, 7, 8, 11, 12, 15, 16, 22]
2222

23-
"Tag indicating the point at corner X0, Y0, Z0."
24-
const corner000 = [1]
2523

26-
"Tag indicating the point at corner X1, Y0, Z0."
27-
const corner100 = [2]
24+
"Tag indicating the line at X, Y0, Z0."
25+
const edgeX00 = [9]
2826

29-
"Tag indicating the point at corner X0, Y1, Z0."
30-
const corner010 = [3]
27+
"Tag indicating the line at X, Y1, Z0."
28+
const edgeX10 = [10]
3129

32-
"Tag indicating the point at corner X1, Y1, Z0."
33-
const corner110 = [4]
30+
"Tag indicating the line at X, Y0, Z1."
31+
const edgeX01 = [11]
3432

35-
"Tag indicating the point at corner X0, Y0, Z1."
36-
const corner001 = [5]
33+
"Tag indicating the line at X, Y1, Z1."
34+
const edgeX11 = [12]
3735

38-
"Tag indicating the point at corner X1, Y0, Z1."
39-
const corner101 = [6]
36+
"Tag indicating the line at X0, Y, Z0."
37+
const edge0Y0 = [13]
4038

41-
"Tag indicating the point at corner X0, Y1, Z1."
42-
const corner011 = [7]
39+
"Tag indicating the line at X1, Y, Z0."
40+
const edge1Y0 = [14]
4341

44-
"Tag indicating the point at corner X1, Y1, Z1."
45-
const corner111 = [8]
42+
"Tag indicating the line at X0, Y, Z1."
43+
const edge0Y1 = [15]
44+
45+
"Tag indicating the line at X1, Y, Z1."
46+
const edge1Y1 = [16]
47+
48+
"Tag indicating the line at X0, Y0, Z."
49+
const edge00Z = [17]
50+
51+
"Tag indicating the line at X1, Y0, Z."
52+
const edge10Z = [18]
53+
54+
"Tag indicating the line at X0, Y1, Z."
55+
const edge01Z = [19]
56+
57+
"Tag indicating the line at X1, Y1, Z."
58+
const edge11Z = [20]
59+
60+
"Tag indicating the point at corner X0, Y0, Z0."
61+
const corner000 = [1]
62+
63+
"Tag indicating the point at corner X1, Y0, Z0."
64+
const corner100 = [2]
65+
66+
"Tag indicating the point at corner X0, Y1, Z0."
67+
const corner010 = [3]
68+
69+
"Tag indicating the point at corner X1, Y1, Z0."
70+
const corner110 = [4]
71+
72+
"Tag indicating the point at corner X0, Y0, Z1."
73+
const corner001 = [5]
74+
75+
"Tag indicating the point at corner X1, Y0, Z1."
76+
const corner101 = [6]
77+
78+
"Tag indicating the point at corner X0, Y1, Z1."
79+
const corner011 = [7]
80+
81+
"Tag indicating the point at corner X1, Y1, Z1."
82+
const corner111 = [8]
4683
end

‎thermo_electro_visco.jl‎

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
using HyperFEM
2+
using HyperFEM: jacobian, solve!
3+
using HyperFEM.ComputationalModels.PostMetrics
4+
using HyperFEM.ComputationalModels.CartesianTags
5+
using Gridap, GridapGmsh, GridapSolvers, DrWatson
6+
using GridapSolvers.NonlinearSolvers
7+
using Gridap.FESpaces
8+
9+
folder = joinpath(@__DIR__, "results", "ThermoViscoElectric")
10+
setupfolder(folder; remove=".vtu")
11+
12+
13+
long = 100 # m
14+
width = 100 # m
15+
thick = 1 # m
16+
domain = (0.0, long, 0.0, width, 0.0, thick)
17+
partition = (8, 2, 2)
18+
geometry = CartesianDiscreteModel(domain, partition)
19+
labels = get_face_labeling(geometry)
20+
21+
add_tag_from_tags!(labels, "fixed", [CartesianTags.edgeX00 ;CartesianTags.edgeX10 ;CartesianTags.edge0Y0 ;CartesianTags.edge1Y0;CartesianTags.corner000; CartesianTags.corner100 ; CartesianTags.corner010 ; CartesianTags.corner110 ])
22+
add_tag_from_tags!(labels, "bottom", CartesianTags.faceZ0)
23+
add_tag_from_tags!(labels, "top", CartesianTags.faceZ1)
24+
25+
26+
# Constitutive model
27+
hyper_elastic_model = NeoHookean3D(λ=10000.0, μ=2.0e5)
28+
viscous_branch_1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0.0, μ=4.0e5); τ=2.0)
29+
viscous_branch_2 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0.0, μ=2.0e4); τ=10.0)
30+
visco_elastic_model = GeneralizedMaxwell(hyper_elastic_model, viscous_branch_1, viscous_branch_2)
31+
elec_model = IdealDielectric(ε=3.5416e-11)
32+
therm_model = ThermalModel(Cv=100.0, θr=293.15, α=1.1631, κ=10.0)
33+
cons_model = ThermoElectroMech_Bonet(therm_model, elec_model, visco_elastic_model)
34+
ku = Kinematics(Mechano, Solid)
35+
ke = Kinematics(Electro, Solid)
36+
kt = Kinematics(Thermo, Solid)
37+
F, _... = get_Kinematics(ku)
38+
E = get_Kinematics(ke)
39+
40+
# Setup integration
41+
order = 1
42+
degree = 2 * order
43+
Ω = Triangulation(geometry)
44+
dΩ = Measure(Ω, degree)
45+
t_end = 1.0 # s
46+
Δt = 0.005 # s
47+
48+
# Dirichlet boundary conditions
49+
evolu(Λ) = 1.0
50+
dir_u_tags = ["fixed"]
51+
dir_u_values = [[0.0, 0.0, 0.0]]
52+
dir_u_timesteps = [evolu]
53+
dirichlet_u = DirichletBC(dir_u_tags, dir_u_values, dir_u_timesteps)
54+
55+
evolφ(Λ) = Λ
56+
dir_φ_tags = ["bottom", "top"]
57+
dir_φ_values = [0.0, 0.002]
58+
dir_φ_timesteps = [evolφ, evolφ]
59+
dirichlet_φ = DirichletBC(dir_φ_tags, dir_φ_values, dir_φ_timesteps)
60+
61+
dirichlet_θ = NothingBC()
62+
dirichlet_bc = MultiFieldBC([dirichlet_u, dirichlet_φ, dirichlet_θ])
63+
64+
# Finite Elements
65+
reffeu = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
66+
reffeφ = ReferenceFE(lagrangian, Float64, order)
67+
reffeθ = ReferenceFE(lagrangian, Float64, order)
68+
69+
# Test FE Spaces
70+
Vu = TestFESpace(geometry, reffeu, dirichlet_u, conformity=:H1)
71+
Vφ = TestFESpace(geometry, reffeφ, dirichlet_φ, conformity=:H1)
72+
Vθ = TestFESpace(geometry, reffeθ, dirichlet_θ, conformity=:H1)
73+
74+
# Trial FE Spaces and state variables
75+
Uu = TrialFESpace(Vu, dirichlet_u)
76+
Uφ = TrialFESpace(Vφ, dirichlet_φ)
77+
Uθ = TrialFESpace(Vθ, dirichlet_θ)
78+
uh⁺ = FEFunction(Uu, zero_free_values(Uu))
79+
φh⁺ = FEFunction(Uφ, zero_free_values(Uφ))
80+
θh⁺ = FEFunction(Uθ, zero_free_values(Uθ))
81+
82+
Uu⁻ = TrialFESpace(Vu, dirichlet_u)
83+
Uφ⁻ = TrialFESpace(Vφ, dirichlet_φ)
84+
Uθ⁻ = TrialFESpace(Vθ, dirichlet_θ)
85+
uh⁻ = FEFunction(Uu⁻, zero_free_values(Uu⁻))
86+
φh⁻ = FEFunction(Uφ⁻, zero_free_values(Uφ⁻))
87+
θh⁻ = FEFunction(Uθ⁻, zero_free_values(Uθ⁻))
88+
89+
Uη⁻ = TrialFESpace(Vθ)
90+
ηh⁻ = FEFunction(Uη⁻, zero_free_values(Uη⁻))
91+
92+
UD⁻ = TrialFESpace(Vθ)
93+
Dh⁻ = FEFunction(UD⁻, zero_free_values(UD⁻))
94+
95+
Fh = F∘∇(uh⁺)'
96+
Fh⁻ = F∘∇(uh⁻)'
97+
A = initializeStateVariables(cons_model, dΩ)
98+
99+
# =================================
100+
# Weak forms: residual and jacobian
101+
# =================================
102+
103+
Ψ, ∂Ψ∂F, ∂Ψ∂E, ∂Ψ∂θ, ∂∂Ψ∂FF, ∂∂Ψ∂EE, ∂∂Ψ∂θθ, ∂∂Ψ∂FE, ∂∂Ψ∂Fθ, ∂∂Ψ∂Eθ, _ = cons_model(Δt=Δt)
104+
D, ∂D∂θ = Dissipation(cons_model, Δt)
105+
η(x...) = -∂Ψ∂θ(x...)
106+
∂η∂θ(x...) = -∂∂Ψ∂θθ(x...)
107+
κ = cons_model.thermo.κ
108+
109+
Mechano_coupling(Λ) = uh⁻ + (uh⁺ - uh⁻) * Λ
110+
Electro_coupling(Λ) = φh⁻ + (φh⁺ - φh⁻) * Λ
111+
Thermo_coupling(Λ) = θh⁻ + (θh⁺ - θh⁻) * Λ
112+
113+
# Electro
114+
res_elec(Λ) = (φ, vφ) -> residual(cons_model, Electro, (ku, ke, kt), (Mechano_coupling(Λ), φ, Thermo_coupling(Λ)), vφ, dΩ, 0.0, Fh⁻, A...; Δt=Δt)
115+
jac_elec(Λ) = (φ, dφ, vφ) -> jacobian(cons_model, Electro, (ku, ke, kt), (Mechano_coupling(Λ), φ, Thermo_coupling(Λ)), dφ, vφ, dΩ, 0.0, Fh⁻, A...; Δt=Δt)
116+
117+
# Mechano
118+
res_mec(Λ) = (u, v) -> residual(cons_model, Mechano, (ku, ke, kt), (u, Electro_coupling(Λ), Thermo_coupling(Λ)), v, dΩ, 0.0, Fh⁻, A...; Δt=Δt)
119+
jac_mec(Λ) = (u, du, v) -> jacobian(cons_model, Mechano, (ku, ke, kt), (u, Electro_coupling(Λ), Thermo_coupling(Λ)), du, v, dΩ, 0.0, Fh⁻, A...; Δt=Δt)
120+
121+
# Thermo
122+
res_therm(Λ) = (θ, vθ) -> begin
123+
uhᵞ = Mechano_coupling(Λ)
124+
φhᵞ = Electro_coupling(Λ)
125+
∫( 1/Δt*(θ*η∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...) - θh⁻*ηh⁻)*vθ +
126+
-1/Δt*0.5*(η∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...) + ηh⁻)*(θ - θh⁻)*vθ +
127+
-0.5*(D∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...) + Dh⁻)*vθ +
128+
0.5*(κ * ∇(θ)·∇(vθ)) + 0.5*(κ * ∇(θh⁻)·∇(vθ))
129+
)dΩ
130+
end
131+
jac_therm(Λ) = (θ, dθ, vθ) -> begin
132+
uhᵞ = Mechano_coupling(Λ)
133+
φhᵞ = Electro_coupling(Λ)
134+
∫( 1/Δt*(η∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...) + θ*∂η∂θ∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...))*dθ*vθ +
135+
-1/Δt*0.5*(∂η∂θ∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...)*(θ - θh⁻) + (η∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...) + ηh⁻))*dθ*vθ +
136+
-0.5*∂D∂θ∘(F∘∇(uhᵞ), E∘∇(φhᵞ), θ, Fh⁻, A...)*dθ*vθ +
137+
0.5*κ*∇(dθ)·∇(vθ)
138+
)dΩ
139+
end
140+
141+
# nonlinear solver electro
142+
ls = LUSolver()
143+
nls = NewtonSolver(ls; maxiter=20, atol=1e-6, rtol=1e-6, verbose=true)
144+
comp_model_elec = StaticNonlinearModel(res_elec, jac_elec, Uφ, Vφ, dirichlet_φ; nls=nls, xh=φh⁺, xh⁻=φh⁻)
145+
146+
# nonlinear solver mechano
147+
comp_model_mec = StaticNonlinearModel(res_mec, jac_mec, Uu, Vu, dirichlet_u; nls=nls, xh=uh⁺, xh⁻=uh⁻)
148+
149+
# nonlinear solver thermo
150+
comp_model_mec = StaticNonlinearModel(res_therm, jac_therm, Uθ, Vθ, dirichlet_θ; nls=nls, xh=θh⁺, xh⁻=θh⁻)
151+
152+
# nonlinear staggered model
153+
comp_model = StaggeredModel((comp_model_elec, comp_model_mec), (φh⁺, uh⁺, θh⁺), (φh⁻, uh⁻, θh⁻))
154+
155+
count = Ref{Int}(0)
156+
157+
# Postprocessor to save results
158+
function driverpost(post)
159+
step = post.iter
160+
pvd = post.cachevtk[3]
161+
filePath = post.cachevtk[2]
162+
count[] += 1
163+
164+
push!(uz, component_LInf(uh⁺, :z, Ω))
165+
if post.cachevtk[1]
166+
pvd[step] = createvtk(Ω, filePath * "/TIME_$(count[])" * ".vtu", cellfields=["u" => uh⁺, "φ" => φh⁺])
167+
end
168+
updateStateVariables!(A, cons_model, Δt, F∘∇(uh⁺), F∘∇(uh⁻))
169+
end
170+
171+
post_model = PostProcessor(comp_model_mec, driverpost; is_vtk=true, filepath=folder)
172+
173+
args_elec = Dict(:stepping => (nsteps=1, maxbisec=1))
174+
args_mec = Dict(:stepping => (nsteps=1, maxbisec=1))
175+
args_therm = Dict(:stepping => (nsteps=1, maxbisec=1), :post=>post_model)
176+
args = (args_elec, args_mec, args_therm)
177+
178+
t = [0:Δt:t_end-Δt...]
179+
uz = Float64[]
180+
nsteps = t_end / Δt
181+
solve!(comp_model; stepping=(nsteps=nsteps, nsubsteps=1, maxbisec=1), kargsolve=args)
182+
183+
plot(t, uz)
184+
185+
@assert uz[end] ≈ 1.39059921e-5
186+
@assert norm(uz) ≈ 8.79961900e-5

0 commit comments

Comments
 (0)