Skip to content

Commit 463985e

Browse files
committed
Added the dual-horizon bond-based formulation
1 parent 33fa0da commit 463985e

File tree

4 files changed

+260
-7
lines changed

4 files changed

+260
-7
lines changed

src/Peridynamics.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,8 @@ end
99
import LibGit2, Dates
1010

1111
# Material models
12-
export BBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial, CKIMaterial
12+
export BBMaterial, DHBBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial,
13+
CKIMaterial
1314

1415
# CMaterial related types
1516
export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling
@@ -75,6 +76,7 @@ abstract type AbstractCorrection end
7576
abstract type AbstractStorage end
7677
abstract type AbstractCondition end
7778
abstract type AbstractBondSystemMaterial{Correction} <: AbstractMaterial end
79+
abstract type AbstractBondBasedMaterial{CM} <: AbstractBondSystemMaterial{CM} end
7880
abstract type AbstractCorrespondenceMaterial{CM,ZEM} <: AbstractBondSystemMaterial{ZEM} end
7981
abstract type AbstractBondAssociatedSystemMaterial <: AbstractMaterial end
8082
abstract type AbstractConstitutiveModel end
@@ -132,6 +134,7 @@ include("time_solvers/velocity_verlet.jl")
132134
include("time_solvers/dynamic_relaxation.jl")
133135

134136
include("physics/bond_based.jl")
137+
include("physics/dh_bond_based.jl")
135138
include("physics/continuum_kinematics_inspired.jl")
136139
include("physics/ordinary_state_based.jl")
137140
include("physics/constitutive_models.jl")

src/physics/bond_based.jl

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
7070
- `damage::Vector{Float64}`: Damage of each point.
7171
- `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
7272
"""
73-
struct BBMaterial{Correction,DM} <: AbstractBondSystemMaterial{Correction}
73+
struct BBMaterial{Correction,DM} <: AbstractBondBasedMaterial{Correction}
7474
dmgmodel::DM
7575
function BBMaterial{C}(dmgmodel::DM) where {C,DM}
7676
new{C,DM}(dmgmodel)
@@ -83,6 +83,14 @@ end
8383
BBMaterial(; kwargs...) = BBMaterial{NoCorrection}(; kwargs...)
8484

8585
function StandardPointParameters(mat::BBMaterial, p::Dict{Symbol,Any})
86+
(; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters_bb(mat, p)
87+
(; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
88+
bc = 18 * K /* δ^4) # bond constant
89+
return StandardPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
90+
end
91+
92+
function get_required_point_parameters_bb(mat::AbstractBondBasedMaterial,
93+
p::Dict{Symbol,Any})
8694
par = get_given_elastic_params(p)
8795
(; E, nu, G, K, λ, μ) = par
8896
if isfinite(nu) && !isapprox(nu, 0.25)
@@ -100,9 +108,7 @@ function StandardPointParameters(mat::BBMaterial, p::Dict{Symbol,Any})
100108
msg *= "Please define either only one or two fitting elastic parameters!\n"
101109
throw(ArgumentError(msg))
102110
end
103-
(; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
104-
bc = 18 * K /* δ^4) # bond constant
105-
return StandardPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
111+
return (; δ, rho, E, nu, G, K, λ, μ)
106112
end
107113

108114
@params BBMaterial StandardPointParameters
@@ -130,8 +136,8 @@ function init_field(::BBMaterial, ::AbstractTimeSolver, system::BondSystem,
130136
end
131137

132138
# Customized calc_failure to save the bond stretch ε for force density calculation
133-
function calc_failure!(storage::BBStorage, system::BondSystem,
134-
::BBMaterial, ::CriticalStretch,
139+
function calc_failure!(storage::AbstractStorage, system::BondSystem,
140+
::AbstractBondBasedMaterial, ::CriticalStretch,
135141
paramsetup::AbstractParameterSetup, i)
136142
(; εc) = get_params(paramsetup, i)
137143
(; position, n_active_bonds, bond_active, bond_stretch) = storage

src/physics/dh_bond_based.jl

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
"""
2+
DHBBMaterial()
3+
DHBBMaterial{Correction}()
4+
5+
A material type used to assign the material of a [`Body`](@ref) with the dual-horizon
6+
bond-based formulation of peridynamics.
7+
8+
# Keywords
9+
- `dmgmodel::AbstractDamageModel`: Damage model defining the fracture behavior.
10+
(default: `CriticalStretch()`)
11+
12+
Possible correction methods are:
13+
- [`NoCorrection`](@ref): No correction is applied. (default)
14+
- [`EnergySurfaceCorrection`](@ref): The energy based surface correction method of
15+
Le and Bobaru (2018) is applied.
16+
17+
# Examples
18+
19+
```julia-repl
20+
julia> mat = DHBBMaterial()
21+
DHBBMaterial{NoCorrection}()
22+
23+
julia> mat = DHBBMaterial{EnergySurfaceCorrection}()
24+
DHBBMaterial{EnergySurfaceCorrection}()
25+
```
26+
---
27+
28+
```julia
29+
DHBBMaterial{Correction}
30+
```
31+
32+
Material type for the dual-horizon bond-based peridynamics formulation.
33+
34+
# Type Parameters
35+
- `Correction`: A correction algorithm type. See the constructor docs for more informations.
36+
- `DM`: A damage model type.
37+
38+
# Allowed material parameters
39+
When using [`material!`](@ref) on a [`Body`](@ref) with `DHBBMaterial`, then the following
40+
parameters are allowed:
41+
Material parameters:
42+
- `horizon::Float64`: Radius of point interactions.
43+
- `rho::Float64`: Density.
44+
Elastic parameters:
45+
- `E::Float64`: Young's modulus.
46+
- `G::Float64`: Shear modulus.
47+
- `K::Float64`: Bulk modulus.
48+
- `lambda::Float64`: 1st Lamé parameter.
49+
- `mu::Float64`: 2nd Lamé parameter.
50+
Fracture parameters:
51+
- `Gc::Float64`: Critical energy release rate.
52+
- `epsilon_c::Float64`: Critical strain.
53+
54+
!!! note "Poisson's ratio and bond-based peridynamics"
55+
In bond-based peridynamics, the Poisson's ratio is limited to 1/4 for 3D simulations.
56+
Therefore, only one additional elastic parameter is required.
57+
Optionally, the specification of a second keyword is allowed, if the parameter
58+
combination results in `nu = 1/4`.
59+
60+
# Allowed export fields
61+
When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
62+
`DHBBMaterial`, the following fields are allowed:
63+
- `position::Matrix{Float64}`: Position of each point.
64+
- `displacement::Matrix{Float64}`: Displacement of each point.
65+
- `velocity::Matrix{Float64}`: Velocity of each point.
66+
- `velocity_half::Matrix{Float64}`: Velocity parameter for Verlet time solver.
67+
- `acceleration::Matrix{Float64}`: Acceleration of each point.
68+
- `b_int::Matrix{Float64}`: Internal force density of each point.
69+
- `b_ext::Matrix{Float64}`: External force density of each point.
70+
- `damage::Vector{Float64}`: Damage of each point.
71+
- `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
72+
"""
73+
struct DHBBMaterial{Correction,DM} <: AbstractBondBasedMaterial{Correction}
74+
dmgmodel::DM
75+
function DHBBMaterial{C}(dmgmodel::DM) where {C,DM}
76+
new{C,DM}(dmgmodel)
77+
end
78+
end
79+
80+
function DHBBMaterial{C}(; dmgmodel::AbstractDamageModel=CriticalStretch()) where {C}
81+
return DHBBMaterial{C}(dmgmodel)
82+
end
83+
DHBBMaterial(; kwargs...) = DHBBMaterial{NoCorrection}(; kwargs...)
84+
85+
function StandardPointParameters(mat::DHBBMaterial{C,D}, p::Dict{Symbol,Any}) where {C,D}
86+
(; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters_bb(mat, p)
87+
(; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
88+
bc = 0.5 * 18 * K /* δ^4) # half of the normal bond constant
89+
return StandardPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
90+
end
91+
92+
@params DHBBMaterial StandardPointParameters
93+
94+
@storage DHBBMaterial struct DHBBStorage <: AbstractStorage
95+
@lthfield position::Matrix{Float64}
96+
@pointfield displacement::Matrix{Float64}
97+
@pointfield velocity::Matrix{Float64}
98+
@pointfield velocity_half::Matrix{Float64}
99+
@pointfield velocity_half_old::Matrix{Float64}
100+
@pointfield acceleration::Matrix{Float64}
101+
@htlfield b_int::Matrix{Float64}
102+
@pointfield b_int_old::Matrix{Float64}
103+
@pointfield b_ext::Matrix{Float64}
104+
@pointfield density_matrix::Matrix{Float64}
105+
@pointfield damage::Vector{Float64}
106+
bond_stretch::Vector{Float64}
107+
bond_active::Vector{Bool}
108+
@pointfield n_active_bonds::Vector{Int}
109+
end
110+
111+
function init_field(::DHBBMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int})
112+
return zeros(3, get_n_points(system))
113+
end
114+
115+
function init_field(::DHBBMaterial, ::AbstractTimeSolver, system::BondSystem,
116+
::Val{:bond_stretch})
117+
return zeros(get_n_bonds(system))
118+
end
119+
120+
function force_density_point!(storage::DHBBStorage, system::BondSystem, ::DHBBMaterial,
121+
params::StandardPointParameters, t, Δt, i)
122+
(; position, bond_stretch, bond_active, b_int) = storage
123+
(; bonds, correction, volume) = system
124+
for bond_id in each_bond_idx(system, i)
125+
bond = bonds[bond_id]
126+
j = bond.neighbor
127+
Δxij = get_vector_diff(position, i, j)
128+
ε = bond_stretch[bond_id]
129+
ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
130+
b = ω * params.bc * ε .* Δxij
131+
update_add_vector!(b_int, i, b * volume[j])
132+
update_add_vector!(b_int, j, -b * volume[i])
133+
end
134+
return nothing
135+
end
136+
137+
function force_density_point!(storage::DHBBStorage, system::BondSystem, ::DHBBMaterial,
138+
paramhandler::ParameterHandler, t, Δt, i)
139+
(; position, bond_stretch, bond_active, b_int) = storage
140+
(; bonds, correction, volume) = system
141+
params_i = get_params(paramhandler, i)
142+
for bond_id in each_bond_idx(system, i)
143+
bond = bonds[bond_id]
144+
j = bond.neighbor
145+
Δxij = get_vector_diff(position, i, j)
146+
ε = bond_stretch[bond_id]
147+
params_j = get_params(paramhandler, j)
148+
ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
149+
b = ω * (params_i.bc + params_j.bc) / 2 * ε .* Δxij
150+
update_add_vector!(b_int, i, b * volume[j])
151+
update_add_vector!(b_int, j, -b * volume[i])
152+
end
153+
return nothing
154+
end
Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
@testitem "Internal force density dual-horizon bond based" begin
2+
ref_position = [0.0 1.0; 0.0 0.0; 0.0 0.0]
3+
volume = [1.0, 1.0]
4+
δ = 1.5
5+
E = 210e9
6+
body = Body(DHBBMaterial(), ref_position, volume)
7+
material!(body; horizon=δ, rho=1, E)
8+
no_failure!(body)
9+
10+
dh = Peridynamics.threads_data_handler(body, VelocityVerlet(steps=1), 1)
11+
chunk = dh.chunks[1]
12+
(; position, b_int) = chunk.storage
13+
(; bonds) = chunk.system
14+
15+
@test position == ref_position
16+
@test b_int == zeros(3, 2)
17+
18+
# Boundary Condition:
19+
# Point 2 with v_z = 1 m/s with Δt = 0.0015 s
20+
position[1, 2] = 1.0015
21+
22+
Peridynamics.calc_force_density!(chunk, 0, 0)
23+
24+
b12 = 18 * E / (3 * (1 - 2 * 0.25)) /* δ^4) * 1.0015 * 0.0015/1.0015 * 1.0
25+
@test b_int [b12 -b12; 0.0 0.0; 0.0 0.0]
26+
end
27+
28+
@testitem "Internal force density dual-horizon bond based with surface correction" begin
29+
ref_position = [0.0 1.0; 0.0 0.0; 0.0 0.0]
30+
volume = [1.0, 1.0]
31+
δ = 1.5
32+
E = 210e9
33+
body = Body(DHBBMaterial{EnergySurfaceCorrection}(), ref_position, volume)
34+
material!(body; horizon=δ, rho=1, E)
35+
no_failure!(body)
36+
37+
ts = VelocityVerlet(steps=1)
38+
dh = Peridynamics.threads_data_handler(body, ts, 1)
39+
Peridynamics.initialize!(dh, ts)
40+
chunk = dh.chunks[1]
41+
(; position, b_int) = chunk.storage
42+
(; bonds) = chunk.system
43+
44+
@test position == ref_position
45+
@test b_int == zeros(3, 2)
46+
47+
# Boundary Condition:
48+
# Point 2 with v_z = 1 m/s with Δt = 0.0015 s
49+
position[1, 2] = 1.0015
50+
51+
Peridynamics.calc_force_density!(chunk, 0, 0)
52+
53+
sc = 2 * 3.1808625617603665 # double the normal surface correction factor
54+
@test all(x -> x sc, chunk.system.correction.scfactor)
55+
56+
b12 = sc * 18 * E / (3 * (1 - 2 * 0.25)) /* δ^4) * 1.0015 * 0.0015/1.0015 * 1.0
57+
@test b_int [b12 -b12; 0.0 0.0; 0.0 0.0]
58+
end
59+
60+
@testitem "Internal force density dual-horizon bond based interface" begin
61+
ref_position = [0.0 1.0; 0.0 0.0; 0.0 0.0]
62+
volume = [1.0, 1.0]
63+
δ = 1.5
64+
Ea = 105e9
65+
Eb = 2 * Ea
66+
E = (Ea + Eb) / 2
67+
body = Body(DHBBMaterial(), ref_position, volume)
68+
point_set!(body, :a, [1])
69+
point_set!(body, :b, [2])
70+
material!(body, :a, horizon=δ, rho=1, E=Ea, Gc=1.0)
71+
material!(body, :b, horizon=δ, rho=1, E=Eb, Gc=1.0)
72+
no_failure!(body)
73+
74+
dh = Peridynamics.threads_data_handler(body, VelocityVerlet(steps=1), 1)
75+
chunk = dh.chunks[1]
76+
(; position, b_int) = chunk.storage
77+
(; bonds) = chunk.system
78+
79+
@test position == ref_position
80+
@test b_int == zeros(3, 2)
81+
82+
# Boundary Condition:
83+
# Point 2 with v_z = 1 m/s with Δt = 0.0015 s
84+
position[1, 2] = 1.0015
85+
86+
Peridynamics.calc_force_density!(chunk, 0, 0)
87+
88+
b12 = 18 * E / (3 * (1 - 2 * 0.25)) /* δ^4) * 1.0015 * 0.0015/1.0015 * 1.0
89+
@test b_int [b12 -b12; 0.0 0.0; 0.0 0.0]
90+
end

0 commit comments

Comments
 (0)