Bond-based fatigue damage model

[Lumos1110]

Hello, do you have an idea for a fatigue model?
It would be great to make this board

[kaipartmann]

Hello @Lumos1110,
Thank you for suggesting this idea; I have already thought about this in the past.
Currently, the main focus of the package development is the v1.0.0 release scheduled for summer.
The entire package was newly written from scratch for this release with a completely new API and MPI support for extensive parallel simulations.
After the new API and features mature, I think this could be a feature worth considering. This model could then be an option and implemented: https://doi.org/10.1098/rsta.2021.0217.

I will mark this as a possible feature. Please feel free to ping me for updates.

[kaipartmann]

Fast draft for the implementation in https://doi.org/10.1098/rsta.2021.0217:

"""
TODO
"""
struct BBFMaterial <: AbstractMaterial end

struct BBFPointParameters <: AbstractPointParameters
    δ::Float64
    rho::Float64
    E::Float64
    nu::Float64
    G::Float64
    K::Float64
    λ::Float64
    μ::Float64
    Gc::Float64
    εc::Float64
    bc::Float64
end

@inline point_param_type(::BBFMaterial) = BBFPointParameters

@inline function allowed_material_kwargs(::BBFMaterial)
    return DEFAULT_POINT_KWARGS
end

function get_point_params(::BBFMaterial, p::Dict{Symbol,Any})
    δ = get_horizon(p)
    rho = get_density(p)
    if haskey(p, :nu)
        msg = "keyword `nu` is not allowed for BBMaterial!\n"
        msg *= "Bond-based peridynamics has a limitation on the possion ratio.\n"
        msg *= "Therefore, when using BBMaterial, `nu` is hardcoded as 1/4.\n"
        throw(ArgumentError(msg))
    else
        p[:nu] = 0.25
    end
    E, nu, G, K, λ, μ = get_elastic_params(p)
    Gc, εc = get_frac_params(p, δ, K)
    bc = 18 * K / (π * δ^4) # bond constant
    return BBFPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
end

@inline discretization_type(::BBFMaterial) = BondDiscretization

@inline function init_discretization(body::Body{BBFMaterial}, args...)
    return init_bond_discretization(body, args...)
end

struct BBFVerletStorage <: AbstractStorage
    position::Matrix{Float64}
    displacement::Matrix{Float64}
    velocity::Matrix{Float64}
    velocity_half::Matrix{Float64}
    acceleration::Matrix{Float64}
    b_int::Matrix{Float64}
    b_ext::Matrix{Float64}
    damage::Vector{Float64}
    bond_active::Vector{Bool}
    n_active_bonds::Vector{Int}
end

const BBFStorage = Union{BBFVerletStorage}

@inline function storage_type(::BBFMaterial, ::VelocityVerlet)
    return BBFVerletStorage
end

function init_storage(::Body{BBFMaterial}, ::VelocityVerlet, bd::BondDiscretization,
                      ch::ChunkHandler)
    n_loc_points = length(ch.loc_points)
    position = copy(bd.position)
    displacement = zeros(3, n_loc_points)
    velocity = zeros(3, n_loc_points)
    velocity_half = zeros(3, n_loc_points)
    acceleration = zeros(3, n_loc_points)
    b_int = zeros(3, n_loc_points)
    b_ext = zeros(3, n_loc_points)
    damage = zeros(n_loc_points)
    bond_active = ones(Bool, length(bd.bonds))
    n_active_bonds = copy(bd.n_neighbors)
    return BBFVerletStorage(position, displacement, velocity, velocity_half, acceleration,
                            b_int, b_ext, damage, bond_active, n_active_bonds)
end

@inline reads_from_halo(::Type{BBFMaterial}) = (:position,)
@inline writes_to_halo(::Type{BBFMaterial}) = ()

function force_density_point!(s::BBFStorage, bd::BondDiscretization, ::BBFMaterial,
                              param::BBFPointParameters, i::Int)
    for bond_id in each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        # current bond length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)

        # bond strain
        ε = (l - L) / L

        # calculate D, eq. (2.6), here:
        D = ...

        # apply the failure mechanism with D,
        # maybe something like this:
        if ε > param.εc && bond.fail_permit
            s.damage[i] = D
        end

        # update of force density
        temp = (1 - D) * param.bc * ε / l * bd.volume[j]
        s.b_int[1, i] += temp * Δxijx
        s.b_int[2, i] += temp * Δxijy
        s.b_int[3, i] += temp * Δxijz
    end
    return nothing
end

# you would have to disable the damage calculation in the timestepping scheme, because it
# would override the calculated D in the s.damage field:
function solve_timestep!(dh::ThreadsDataHandler{B{BBFMaterial,P,D,S}},
                         options::ExportOptions, Δt::Float64, Δt½::Float64,
                         n::Int) where {B<:AbstractBodyChunk,P,D,S}
    t = n * Δt
    @threads :static for chunk_id in eachindex(dh.chunks)
        chunk = dh.chunks[chunk_id]
        update_vel_half!(chunk, Δt½)
        apply_bcs!(chunk, t)
        update_disp_and_pos!(chunk, Δt)
    end
    @threads :static for chunk_id in eachindex(dh.chunks)
        exchange_read_fields!(dh, chunk_id)
        calc_force_density!(dh.chunks[chunk_id])
    end
    @threads :static for chunk_id in eachindex(dh.chunks)
        exchange_write_fields!(dh, chunk_id)
        chunk = dh.chunks[chunk_id]
        # calc_damage!(chunk)
        update_acc_and_vel!(chunk, Δt½)
        export_results(dh, options, chunk_id, n, t)
    end
    return nothing
end

[Lumos1110]

Okay, I'm looking forward to the new version.
May I ask what is your e-mail address?

[kaipartmann]

kai.partmann@uni-siegen.de

[Lumos1110]

Ok，i know.
Thank you!

[Lumos1110]

Did you get my e-mail, please? Maybe the title is too random to be considered spam :D

[kaipartmann]

Did you get my e-mail, please? Maybe the title is too random to be considered spam :D

You are right; it was filtered as spam, so I did not see it - just replied :D