Skip to content

Commit d20c76b

Browse files
ShiWeiHuHShiWeiHuH
committed
new module: add temperature-dependent thermal properties (kp_T, cv_T, aph_T) with tests and examples
1 parent 92320fe commit d20c76b

File tree

4 files changed

+281
-2
lines changed

4 files changed

+281
-2
lines changed

‎src/HyperFSI.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ export Post2D, Post3D
2323
#export glob, write_sprata_bc_file_2d, read_bc_from_sparta, write_sparta_files, writedlm, readdlm
2424

2525
# New Material models
26-
export BBTMaterial, BBTMMaterial
26+
export BBTMaterial, BBTMMaterial, BBTTMaterial
2727

2828

2929
# New Discretization
@@ -45,6 +45,7 @@ include("structure/time_solvers/dual_timesteps.jl")
4545
include("structure/physics/bond_based_thermal_diffusion.jl")
4646
include("structure/physics/bond_based_thermomechanics.jl")
4747
include("structure/physics/bond_based_dualstep_thermomechanics.jl")
48+
include("structure/physics/bond_based_thermal_diffusion_temp_dependent.jl")
4849

4950
include("fluid/FlowTimesolver.jl")
5051
include("fluid/Evolution.jl")

‎src/structure/physics/bond_based_thermal_diffusion.jl‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -160,7 +160,7 @@ function Peridynamics.Peridynamics.allowed_material_kwargs(::BBTMaterial)
160160
return (thermal_kwargs())
161161
end
162162

163-
# Thermal
163+
#= Thermal
164164
@inline function update_temperature!(b::Peridynamics.AbstractBodyChunk, Δt::Float64)
165165
for point_id in Peridynamics.each_point_idx(b)
166166
param = Peridynamics.get_params(b, point_id)
@@ -169,6 +169,21 @@ end
169169
end
170170
return nothing
171171
end
172+
=#
173+
174+
@inline function update_temperature!(b::Peridynamics.AbstractBodyChunk, Δt::Float64)
175+
for point_id in Peridynamics.each_point_idx(b)
176+
param = Peridynamics.get_params(b, point_id)
177+
if isdefined(param, :cv_T)
178+
cv = param.cv_T(param.cv, b.storage.temperature[1, point_id])
179+
else
180+
cv = param.cv
181+
end
182+
k = Δt / (param.rho * cv)
183+
_update_temperature!(b.storage.temperature, b.storage.pflux, b.storage.hsource, k, point_id)
184+
end
185+
return nothing
186+
end
172187

173188
@inline function _update_temperature!(temperature, pflux, hsource, k, i)
174189
temperature[1, i] += (pflux[1, i] + hsource[1, i]) * k

‎src/structure/physics/bond_based_thermal_diffusion_temp_dependent.jl‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
"""
2+
Bond-Based Thermal Diffusion Module with Temperature-Dependent Properties
3+
4+
This module extends the standard bond-based thermal diffusion model
5+
by incorporating temperature-dependent thermal properties, including:
6+
7+
Thermal conductivity in PD (kp)
8+
Specific heat capacity (cv)
9+
Thermal expansion coefficient (aph).
10+
11+
Customizable Functions:
12+
kp_T(kc, kp, T): 3-parameter function (classical Thermal conductivity, PD Thermal conductivity , temperature T).
13+
cv_T(cv, T): 2-parameter function (specific heat cv, temperature T).
14+
aph_T(aph, T): 2-parameter function (expansion coefficient aph, temperature T).
15+
16+
Default Behavior:
17+
If no function is assigned (e.g., kp_T is undefined), the property is treated as temperature-independent.
18+
"""
19+
20+
@inline thermal_temp_kwargs() = (:E, :nu, :G, :K, :lambda, :mu, :rho, :horizon,
21+
:Gc, :epsilon_c, :kc, :aph, :cv, :rft, :h, :hσ, :hϵ, :tem∞, :thick, :kp_T, :cv_T, :aph_T)
22+
23+
struct BBTTMaterial{Correction,DM} <: Peridynamics.AbstractBondSystemMaterial{Correction}
24+
dmgmodel::DM
25+
function BBTTMaterial{C}(dmgmodel::DM) where {C,DM} #BB_Temperature_dependent_Theermal
26+
new{C,DM}(dmgmodel)
27+
end
28+
end
29+
30+
function BBTTMaterial{C}(; dmgmodel::Peridynamics.AbstractDamageModel=CriticalStretch()) where {C}
31+
return BBTTMaterial{C}(dmgmodel)
32+
end
33+
34+
BBTTMaterial(; kwargs...) = BBTTMaterial{NoCorrection}(; kwargs...)
35+
36+
struct BBTTPointParameters <: Peridynamics.AbstractPointParameters
37+
δ::Float64
38+
rho::Float64
39+
E::Float64
40+
nu::Float64
41+
G::Float64
42+
K::Float64
43+
λ::Float64
44+
μ::Float64
45+
Gc::Float64
46+
εc::Float64
47+
bc::Float64
48+
kc::Float64 # thermal conductivity at room temperature
49+
kp::Float64 # microconductivity at room temperature
50+
aph::Float64 # thermal expansion at room temperature
51+
cv::Float64 # specific heat capacity at room temperature
52+
rft::Float64 # Reference temperature at room temperature
53+
h::Float64 # convective heat transfer coefficient,
54+
hσ::Float64 # Stefan-Boltzman constant,
55+
hϵ::Float64 # emissivity
56+
tem∞::Float64 # temperature of the surrounding medium
57+
kp_T::Function # function describing kp variation with temperature f(kp,T)
58+
cv_T::Function # function describing kp variation with temperature f(kp,T)
59+
aph_T::Function # function describing kp variation with temperature f(kp,T)
60+
end
61+
62+
function const_kp(kc, kp, T)
63+
return kp
64+
end
65+
66+
function const_cv(cv, T)
67+
return cv
68+
end
69+
70+
function const_aph(aph, T)
71+
return aph
72+
end
73+
74+
function BBTTPointParameters(mat::BBTTMaterial, p::Dict{Symbol,Any})
75+
par = Peridynamics.get_given_elastic_params(p)
76+
(; E, nu, G, K, λ, μ) = par
77+
78+
if haskey(p, :thick)
79+
p[:nu] = 1/3
80+
else
81+
p[:nu] = 1/4
82+
end
83+
84+
(; δ, rho, E, nu, G, K, λ, μ) = Peridynamics.get_required_point_parameters(mat, p)
85+
Gc, εc = Peridynamics.get_frac_params(mat.dmgmodel, p, δ, K)
86+
kc, aph, cv, rft, h, hσ, hϵ, tem∞ = get_thermal_params(p, δ)
87+
88+
if haskey(p, :thick) #2D
89+
thick = float(p[:thick])
90+
bc = 9 * E / (π * thick * δ^3) # bond constant
91+
kp = 6 * kc / (π * thick * δ^3) # microcndicitvity constant
92+
else #3D
93+
bc = 12 * E / (π * δ^4)
94+
kp = 6 * kc / (π * δ^4)
95+
end
96+
97+
kp_T = get(p, :kp_T, const_kp)
98+
cv_T = get(p, :cv_T, const_cv)
99+
aph_T = get(p, :aph_T, const_aph)
100+
101+
return BBTTPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc, kc, kp, aph, cv, rft, h, hσ, hϵ, tem∞, kp_T, cv_T, aph_T)
102+
end
103+
104+
@Peridynamics.params BBTTMaterial BBTTPointParameters
105+
106+
@Peridynamics.storage BBTTMaterial struct BBTTStorage <: Peridynamics.AbstractStorage
107+
@lthfield position::Matrix{Float64}
108+
@pointfield displacement::Matrix{Float64}
109+
@pointfield velocity::Matrix{Float64}
110+
@pointfield velocity_half::Matrix{Float64}
111+
@pointfield velocity_half_old::Matrix{Float64}
112+
@pointfield acceleration::Matrix{Float64}
113+
@pointfield b_int::Matrix{Float64}
114+
@pointfield b_int_old::Matrix{Float64}
115+
@pointfield b_ext::Matrix{Float64}
116+
@pointfield density_matrix::Matrix{Float64}
117+
@pointfield damage::Vector{Float64}
118+
bond_stretch::Vector{Float64}
119+
bond_active::Vector{Bool}
120+
@pointfield n_active_bonds::Vector{Int}
121+
@lthfield temperature::Matrix{Float64}
122+
@pointfield pflux::Matrix{Float64}
123+
@pointfield hsource::Matrix{Float64}
124+
end
125+
126+
function Peridynamics.Peridynamics.init_field(::BBTTMaterial, ::Peridynamics.AbstractTimeSolver, system::Peridynamics.BondSystem,
127+
::Val{:bond_stretch})
128+
return zeros(Peridynamics.get_n_bonds(system))
129+
end
130+
131+
function Peridynamics.Peridynamics.allowed_material_kwargs(::BBTTMaterial)
132+
return (thermal_temp_kwargs())
133+
end
134+
135+
function pflux_point!(storage::BBTTStorage, system::Peridynamics.BondSystem,
136+
::BBTTMaterial, param::BBTTPointParameters, i::Int, mbd_t::Vector{Float64})
137+
138+
for bond_id in system.bond_ids[i]
139+
bond = system.bonds[bond_id]
140+
j, L = bond.neighbor, bond.length
141+
142+
mof_th = mbd_t[bond_id]
143+
144+
Δtem = storage.temperature[1, j] - storage.temperature[1, i]
145+
146+
kp_tm = 0.5*(param.kp_T(param.kc, param.kp, storage.temperature[1, i]) + param.kp_T(param.kc, param.kp, storage.temperature[1, j]))
147+
148+
storage.pflux[1, i] += storage.bond_active[bond_id] * Δtem / L * kp_tm * system.volume[j] * mof_th
149+
end
150+
return nothing
151+
end
152+
153+
function pflux_point!(storage::BBTTStorage, system::Peridynamics.BondSystem,
154+
::BBTTMaterial, paramhandler::Peridynamics.ParameterHandler, i::Int, mbd_t::Vector{Float64})
155+
156+
params_i = Peridynamics.get_params(paramhandler, i)
157+
for bond_id in system.bond_ids[i]
158+
bond = system.bonds[bond_id]
159+
j, L = bond.neighbor, bond.length
160+
161+
mof_th = mbd_t[bond_id]
162+
163+
Δtem = storage.temperature[1, j] - storage.temperature[1, i]
164+
165+
params_j = Peridynamics.get_params(paramhandler, j)
166+
167+
kp_tm = 0.5*(params_j.kp_T(params_j.kc, params_j.kp, storage.temperature[1, j]) + params_i.kp_T(params_i.kc, params_i.kp, storage.temperature[1, i]))
168+
169+
storage.pflux[1, i] += storage.bond_active[bond_id] * kp_tm * Δtem / L * system.volume[j] * mof_th
170+
end
171+
return nothing
172+
end

‎test/structure/physics/test_bond_based_thermal_diffusion_temp_dependent.jl‎

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
@testitem "BBTTMaterial_const_homogeneous" begin
2+
using HyperFSI: BBTTMaterial
3+
4+
pos = [ -0.25 0.25 -0.25 0.25
5+
-0.25 -0.25 0.25 0.25
6+
0.0 0.0 0.0 0.0]
7+
vol = [0.125, 0.125, 0.125, 0.125]
8+
height = 0.5
9+
body = Body(BBTTMaterial(), pos, vol)
10+
Peridynamics.material!(body; horizon=2, E=2.0e11, rho=1.0, epsilon_c=1e-2, kc=1.0, aph=1.0e-6, cv=1.0, rft=273,
11+
h=25.0, hσ=5.73e-8, hϵ=0.9, tem∞=0.0, thick=height)
12+
13+
@test body.mat isa BBTTMaterial
14+
@test body.point_params[1].kp_T == HyperFSI.const_kp
15+
@test body.point_params[1].cv_T == HyperFSI.const_cv
16+
@test body.point_params[1].aph_T == HyperFSI.const_aph
17+
@test HyperFSI.const_kp(1.0, 2.0, 300.0) == 2.0
18+
@test HyperFSI.const_kp(1.0, 2.0, 1000.0) == HyperFSI.const_kp(1.0, 2.0, 300.0)
19+
@test HyperFSI.const_cv(1.0, 300.0) == 1.0
20+
@test HyperFSI.const_cv(1.0, 1000.0) == HyperFSI.const_cv(1.0, 300.0)
21+
@test HyperFSI.const_aph(1.0e-6, 300.0) == 1.0e-6
22+
@test HyperFSI.const_aph(1.0e-6, 1000.0) == HyperFSI.const_aph(1.0e-6, 300.0)
23+
end
24+
25+
@testitem "BBTTMaterial_k(T)_homogeneous" begin
26+
using HyperFSI: BBTTMaterial
27+
28+
standard_thermal_params = Dict{Symbol, Float64}(
29+
:k => 1.0,
30+
:cv => 1.0,
31+
:rho => 1.0,
32+
:δ => 2.0
33+
)
34+
pos = [ -0.25 0.25 -0.25 0.25 -0.25 0.25 -0.25 0.25
35+
-0.25 -0.25 0.25 0.25 -0.25 -0.25 0.25 0.25
36+
-0.25 -0.25 -0.25 -0.25 0.25 0.25 0.25 0.25]
37+
vol = [0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
38+
39+
function k_linear(kc, kp, T)
40+
kc_act = (1.0 + (T - 300)/300) # linear increase of k with T
41+
kp_act = kc_act * kp / kc
42+
return kp_act
43+
end
44+
45+
body = Body(BBTTMaterial(), pos, vol)
46+
material!(body; horizon=2, E=2.0e11, rho=1.0, epsilon_c=1e-2, kc=1.0, aph=1.0e-6, cv=1.0, rft=273,
47+
h=25.0, hσ=5.73e-8, hϵ=0.9, tem∞=0.0, kp_T=k_linear)
48+
49+
@test body.point_params[1].kp_T == k_linear
50+
@test body.point_params[1].cv_T == HyperFSI.const_cv
51+
@test body.point_params[1].aph_T == HyperFSI.const_aph
52+
@test body.point_params[1].kc == standard_thermal_params[:k]
53+
@test body.point_params[1].kp == 6*standard_thermal_params[:k] / (π * standard_thermal_params[:δ]^4)
54+
@test k_linear(1.0, body.point_params[1].kp, 300.0) == body.point_params[1].kp
55+
@test k_linear(1.0, body.point_params[1].kp, 600.0) == 2*body.point_params[1].kp
56+
end
57+
58+
@testitem "BBTTMaterial_k(T)_cv(T)_homogeneous" begin
59+
using HyperFSI: BBTTMaterial
60+
standard_thermal_params = Dict{Symbol, Float64}(
61+
:k => 1.0,
62+
:cv => 1.0,
63+
:rho => 1.0,
64+
:δ => 2.0
65+
)
66+
pos = [ -0.25 0.25 -0.25 0.25 -0.25 0.25 -0.25 0.25
67+
-0.25 -0.25 0.25 0.25 -0.25 -0.25 0.25 0.25
68+
-0.25 -0.25 -0.25 -0.25 0.25 0.25 0.25 0.25]
69+
vol = [0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
70+
71+
function k_linear(kc, kp, T)
72+
kc_act = (1.0 + (T - 300)/300) # linear increase of k with T
73+
kp_act = kc_act * kp / kc
74+
return kp_act
75+
end
76+
77+
cv_linear(cv, T) = (1.0 + (T - 300)/300) # linear increase of cv with T
78+
body = Body(BBTTMaterial(), pos, vol)
79+
material!(body; horizon=2, E=2.0e11, rho=1.0, epsilon_c=1e-2, kc=1.0, aph=1.0e-6, cv=1.0, rft=273,
80+
h=25.0, hσ=5.73e-8, hϵ=0.9, tem∞=0.0, kp_T=k_linear, cv_T=cv_linear)
81+
82+
@test body.point_params[1].kp_T == k_linear
83+
@test body.point_params[1].cv_T == cv_linear
84+
@test body.point_params[1].aph_T == HyperFSI.const_aph
85+
@test body.point_params[1].kc == standard_thermal_params[:k]
86+
@test body.point_params[1].kp == 6*standard_thermal_params[:k] / (π * standard_thermal_params[:δ]^4)
87+
@test k_linear(1.0, body.point_params[1].kp, 300.0) == body.point_params[1].kp
88+
@test k_linear(1.0, body.point_params[1].kp, 600.0) == 2*body.point_params[1].kp
89+
@test cv_linear(1.0, 300.0) == body.point_params[1].cv
90+
@test cv_linear(1.0, 600.0) == 2*body.point_params[1].cv
91+
end

0 commit comments

Comments
 (0)