Skip to content

Commit 0d0f657

Browse files
kaipartmannkaipartmann
authored
Merge pull request #239 from kaipartmann/energy_surface_correction
Surface correction with model-specific strain energy density
2 parents 4462dd9 + d16269e commit 0d0f657

File tree

4 files changed

+209
-112
lines changed

4 files changed

+209
-112
lines changed

‎src/discretization/bond_system_corrections.jl‎

Lines changed: 44 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -75,81 +75,66 @@ function initialize!(dh::AbstractMPIBodyDataHandler{BondSystem{EnergySurfaceCorr
7575
end
7676

7777
function calc_mfactor!(chunk::AbstractBodyChunk{BondSystem{EnergySurfaceCorrection}})
78-
(; system, storage, paramsetup) = chunk
78+
(; mat, system, storage, paramsetup) = chunk
7979
(; mfactor) = system.correction
80-
ka, a = 1.001, 0.001
80+
λ = 1.001 # deformation stretch factor
8181
for d in 1:3
82-
defposition = copy(system.position)
83-
defposition[d, :] .*= ka
82+
# reset to undeformed positions
83+
storage.position .= system.position
84+
# apply uniform deformation only in dimension d
85+
@views storage.position[d, :] .= λ .* system.position[d, :]
86+
# needed to calculate initial bond lengths
87+
calc_force_density!(storage, system, mat, paramsetup, 0.0, 0.0) # dummy t, Δt
88+
# calculate the mfactor for dimension d
8489
for i in each_point_idx(chunk)
85-
stendens, E_mean, nu_mean = stendens_point(system, paramsetup, storage,
86-
defposition, i)
87-
mfactor[d, i] = analytical_stendens(E_mean, nu_mean, a) / stendens
90+
lame = get_averaged_lame_parameters(system, storage, paramsetup, i)
91+
Ψ_theory = stendens_uniext_small_strain(lame, λ)
92+
strain_energy_density_point!(storage, system, mat, paramsetup, i)
93+
Ψ_pd = storage.strain_energy_density[i]
94+
mfactor[d, i] = Ψ_theory / Ψ_pd
8895
end
96+
# this is needed since the position is modified during the calculation
97+
storage.position .= system.position
98+
# also the artificial strain energy density and b_int should be reset to zero
99+
storage.strain_energy_density .= 0.0
100+
storage.b_int .= 0.0
89101
end
90102
return nothing
91103
end
92104

93-
function stendens_point(system::BondSystem{C}, params::AbstractPointParameters,
94-
storage::AbstractStorage, defposition::AbstractMatrix{Float64},
95-
i::Int) where {C<:EnergySurfaceCorrection}
96-
stendens = 0.0
97-
for bond_id in each_bond_idx(system, i)
98-
bond = system.bonds[bond_id]
99-
j, L = bond.neighbor, bond.length
100-
Δxij = get_vector_diff(defposition, i, j)
101-
l = norm(Δxij)
102-
ε = (l - L) / L
103-
failure = storage.bond_active[bond_id]
104-
stendens += failure * 0.25 * params.bc * ε * ε * L * system.volume[j]
105-
end
106-
return stendens, params.E, params.nu
105+
function get_averaged_lame_parameters(::BondSystem, ::AbstractStorage,
106+
params::AbstractPointParameters, i)
107+
return SVector{2,Float64}(params.λ, params.μ)
107108
end
108109

109-
function stendens_point(system::BondSystem{C}, paramhandler::AbstractParameterHandler,
110-
storage::AbstractStorage, defposition::AbstractMatrix{Float64},
111-
i::Int) where {C<:EnergySurfaceCorrection}
112-
stendens = 0.0
113-
params_i = get_params(paramhandler, i)
114-
E_mean, nu_mean = 0.0, 0.0
110+
function get_averaged_lame_parameters(system::BondSystem, storage::AbstractStorage,
111+
paramsetup::AbstractParameterHandler, i)
112+
(; bonds) = system
113+
(; bond_active) = storage
114+
lame = zero(SVector{2,Float64}) # lame[1] = λ, lame[2] = μ
115+
params_i = get_params(paramsetup, i)
116+
n_active_bonds = 0
115117
for bond_id in each_bond_idx(system, i)
116-
bond = system.bonds[bond_id]
117-
j, L = bond.neighbor, bond.length
118-
Δxij = get_vector_diff(defposition, i, j)
119-
l = norm(Δxij)
120-
ε = (l - L) / L
121-
params_j = get_params(paramhandler, j)
122-
bc_mean = (params_i.bc + params_j.bc) / 2
123-
failure = storage.bond_active[bond_id]
124-
stendens += failure * 0.25 * bc_mean * ε * ε * L * system.volume[j]
125-
E_mean += failure * (params_i.E + params_j.E) / 2
126-
nu_mean += failure * (params_i.nu + params_j.nu) / 2
118+
if bond_active[bond_id]
119+
bond = bonds[bond_id]
120+
j = bond.neighbor
121+
params_j = get_params(paramsetup, j)
122+
λ = (params_i.λ + params_j.λ) / 2
123+
μ = (params_i.μ + params_j.μ) / 2
124+
lame += SVector{2,Float64}(λ, μ)
125+
n_active_bonds += 1
126+
end
127127
end
128-
E_mean /= storage.n_active_bonds[i]
129-
nu_mean /= storage.n_active_bonds[i]
130-
return stendens, E_mean, nu_mean
128+
lame /= n_active_bonds
129+
return lame
131130
end
132131

133-
@inline function analytical_stendens(E, nu, a)
134-
return E / (2 * (1 + nu)) * (nu / (1 - 2 * nu) + 1) * a^2
132+
@inline function stendens_uniext_small_strain(lame, λ)
133+
return (0.5 * lame[1] + lame[2]) * (λ - 1)^2
135134
end
136135

137-
function calc_stendens!(stendens, defposition, chunk)
138-
system = chunk.system
139-
for i in each_point_idx(chunk)
140-
params = get_params(chunk, i)
141-
for bond_id in each_bond_idx(system, i)
142-
bond = system.bonds[bond_id]
143-
j, L = bond.neighbor, bond.length
144-
Δxijx = defposition[1, j] - defposition[1, i]
145-
Δxijy = defposition[2, j] - defposition[2, i]
146-
Δxijz = defposition[3, j] - defposition[3, i]
147-
l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
148-
ε = (l - L) / L
149-
stendens[i] += 0.25 * params.bc * ε * ε * L * system.volume[j]
150-
end
151-
end
152-
return nothing
136+
@inline function stendens_uniext_finite_strain(lame, λ)
137+
return (lame[1] + 2 * lame[2]) / 8 * (λ^2 - 1)^2
153138
end
154139

155140
@inline get_mfactor(chunk::AbstractBodyChunk) = chunk.system.correction.mfactor

‎test/discretization/test_bond_system_corrections.jl‎

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
@testitem "EnergySurfaceCorrection BBMaterial" begin
2+
Δx = 0.1
3+
pos, vol = uniform_box(1, 1, 1, Δx)
4+
horizon = 3.01 * Δx
5+
rho = 8000
6+
E = 210e9
7+
nu = 0.25
8+
mat = BBMaterial{EnergySurfaceCorrection}()
9+
body = Body(mat, pos, vol)
10+
material!(body; horizon, rho, E, nu)
11+
ts = VelocityVerlet(steps=1)
12+
13+
dh = Peridynamics.threads_data_handler(body, ts, 1)
14+
Peridynamics.initialize!(dh, ts)
15+
chunk = dh.chunks[1]
16+
(; system) = chunk
17+
(; correction) = system
18+
(; mfactor, scfactor) = correction
19+
20+
@test minimum(mfactor[1, :]) ≈ 1.0 atol=0.08
21+
@test maximum(mfactor[1, :]) ≈ 3.7 atol=1.0
22+
@test minimum(mfactor[2, :]) ≈ 1.0 atol=0.08
23+
@test maximum(mfactor[2, :]) ≈ 3.7 atol=1.0
24+
@test minimum(mfactor[3, :]) ≈ 1.0 atol=0.08
25+
@test maximum(mfactor[3, :]) ≈ 3.7 atol=1.0
26+
@test minimum(scfactor) ≈ 1.0 atol=0.08
27+
@test maximum(scfactor) ≈ 3.5 atol=1.0
28+
end
29+
30+
@testitem "EnergySurfaceCorrection OSBMaterial" begin
31+
Δx = 0.1
32+
pos, vol = uniform_box(1, 1, 1, Δx)
33+
horizon = 3.01 * Δx
34+
rho = 8000
35+
E = 210e9
36+
nu = 0.25
37+
mat = OSBMaterial{EnergySurfaceCorrection}()
38+
body = Body(mat, pos, vol)
39+
material!(body; horizon, rho, E, nu)
40+
ts = VelocityVerlet(steps=1)
41+
42+
dh = Peridynamics.threads_data_handler(body, ts, 1)
43+
Peridynamics.initialize!(dh, ts)
44+
chunk = dh.chunks[1]
45+
(; system) = chunk
46+
(; correction) = system
47+
(; mfactor, scfactor) = correction
48+
49+
@test minimum(mfactor[1, :]) ≈ 0.75 atol=0.2
50+
@test maximum(mfactor[1, :]) ≈ 1.5 atol=0.3
51+
@test minimum(mfactor[2, :]) ≈ 0.75 atol=0.2
52+
@test maximum(mfactor[2, :]) ≈ 1.5 atol=0.3
53+
@test minimum(mfactor[3, :]) ≈ 0.75 atol=0.2
54+
@test maximum(mfactor[3, :]) ≈ 1.5 atol=0.3
55+
@test minimum(scfactor) ≈ 0.75 atol=0.2
56+
@test maximum(scfactor) ≈ 1.5 atol=0.3
57+
end
58+
59+
@testitem "EnergySurfaceCorrection DHBBMaterial" begin
60+
Δx = 0.1
61+
pos, vol = uniform_box(1, 1, 1, Δx)
62+
horizon = 3.01 * Δx
63+
rho = 8000
64+
E = 210e9
65+
nu = 0.25
66+
mat = DHBBMaterial{EnergySurfaceCorrection}()
67+
body = Body(mat, pos, vol)
68+
material!(body; horizon, rho, E, nu)
69+
ts = VelocityVerlet(steps=1)
70+
71+
dh = Peridynamics.threads_data_handler(body, ts, 1)
72+
Peridynamics.initialize!(dh, ts)
73+
chunk = dh.chunks[1]
74+
(; system) = chunk
75+
(; correction) = system
76+
(; mfactor, scfactor) = correction
77+
78+
@test minimum(mfactor[1, :]) ≈ 1.0 atol=0.08
79+
@test maximum(mfactor[1, :]) ≈ 3.7 atol=1.0
80+
@test minimum(mfactor[2, :]) ≈ 1.0 atol=0.08
81+
@test maximum(mfactor[2, :]) ≈ 3.7 atol=1.0
82+
@test minimum(mfactor[3, :]) ≈ 1.0 atol=0.08
83+
@test maximum(mfactor[3, :]) ≈ 3.7 atol=1.0
84+
@test minimum(scfactor) ≈ 1.0 atol=0.08
85+
@test maximum(scfactor) ≈ 3.5 atol=1.0
86+
end
87+
88+
@testitem "Analytical strain energy density functions" begin
89+
test_cases = [
90+
(210e9, 0.25, 1.01),
91+
(200e9, 0.30, 1.05),
92+
(150e9, 0.20, 1.10),
93+
(300e9, 0.35, 1.02),
94+
]
95+
for (E, nu, λ) in test_cases
96+
λ_lame = E * nu / ((1 + nu) * (1 - 2 * nu))
97+
μ_lame = E / (2 * (1 + nu))
98+
lame = [λ_lame, μ_lame]
99+
Ψ_small = 1/2 * λ_lame * (λ - 1)^2 + μ_lame * (λ - 1)^2
100+
Ψ_finite = 1/8 * λ_lame * (λ^2 - 1)^2 + 1/4 * μ_lame * (λ^2 - 1)^2
101+
@test Peridynamics.stendens_uniext_small_strain(lame, λ) ≈ Ψ_small
102+
@test Peridynamics.stendens_uniext_finite_strain(lame, λ) ≈ Ψ_finite
103+
end
104+
end
105+
106+
@testitem "Get averaged Lamé parameters" begin
107+
pos = [0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0
108+
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
109+
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
110+
vol = [1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]
111+
horizon = 1.01
112+
rho = 8000
113+
E = 105e9
114+
nu = 0.25
115+
mat = BBMaterial{EnergySurfaceCorrection}()
116+
body = Body(mat, pos, vol)
117+
point_set!(x -> x < 4.5, body, :left)
118+
point_set!(x -> x > 4.5, body, :right)
119+
material!(body, :left; horizon, rho, E, nu)
120+
material!(body, :right; horizon, rho, E=2E, nu)
121+
ts = VelocityVerlet(steps=1)
122+
123+
dh = Peridynamics.threads_data_handler(body, ts, 1)
124+
Peridynamics.initialize!(dh, ts)
125+
chunk = dh.chunks[1]
126+
(; system, storage, paramsetup) = chunk
127+
128+
params1 = Peridynamics.get_params(paramsetup, 1)
129+
params4 = Peridynamics.get_params(paramsetup, 4)
130+
params5 = Peridynamics.get_params(paramsetup, 5)
131+
params6 = Peridynamics.get_params(paramsetup, 6)
132+
params7 = Peridynamics.get_params(paramsetup, 7)
133+
params10 = Peridynamics.get_params(paramsetup, 10)
134+
135+
# point #1
136+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 1)
137+
@test lame[1] ≈ params1.λ
138+
@test lame[2] ≈ params1.μ
139+
140+
# point #4
141+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 4)
142+
@test lame[1] ≈ params4.λ
143+
@test lame[2] ≈ params4.μ
144+
145+
# point #5
146+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 5)
147+
@test lame[1] ≈ (params4.λ + 2 * params5.λ + params6.λ) / 4
148+
@test lame[2] ≈ (params4.μ + 2 * params5.μ + params6.μ) / 4
149+
150+
# point #6
151+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 6)
152+
@test lame[1] ≈ (params5.λ + 2 * params6.λ + params7.λ) / 4
153+
@test lame[2] ≈ (params5.μ + 2 * params6.μ + params7.μ) / 4
154+
155+
# point #7
156+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 7)
157+
@test lame[1] ≈ params7.λ
158+
@test lame[2] ≈ params7.μ
159+
160+
# point #10
161+
lame = Peridynamics.get_averaged_lame_parameters(system, storage, paramsetup, 10)
162+
@test lame[1] ≈ params10.λ
163+
@test lame[2] ≈ params10.μ
164+
end

‎test/integration/b_int_dh_bond_based.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ end
5050

5151
Peridynamics.calc_force_density!(chunk, 0, 0)
5252

53-
sc = 2 * 3.1808625617603665 # double the normal surface correction factor
53+
sc = 3.1808625617603665 # double the normal surface correction factor
5454
@test all(x -> x ≈ sc, chunk.system.correction.scfactor)
5555

5656
b12 = sc * 18 * E / (3 * (1 - 2 * 0.25)) / (π * δ^4) * 1.0015 * 0.0015/1.0015 * 1.0

‎test/integration/energy_surface_correction.jl‎

Lines changed: 0 additions & 52 deletions
This file was deleted.

0 commit comments

Comments
 (0)