thermo-mechanical coupling

[vavrines]

Hi @kaipartmann

We're working on the thermo-mechanical coupling model. A script from scratch is here (https://github.com/vavrines/KitPD.jl/blob/main/couple/script.jl). We're just about to clean and refactor the codes. The idea is to reuse the structs and methods here as much as possible. I'd like to know if you're interested in implementing this functionality. If things go well, we should be able to contribute codes, but your assistance will likely be needed.

[kaipartmann]

Hello @vavrines, Thank you for reaching out!

I am very interested in implementing this into Peridynamics.jl. Your code is very well structured and comfortable to read, so implementing it using your example should be feasible. Do you have literature describing your thermo-mechanical coupling model?

Would you also be interested in us collaborating on a publication about this topic? We have good HPC access and can run large-scale simulations once the model works with the package.

[vavrines]

@kaipartmann Thanks for the feedback! Basically we follow Chapter 13 in the monograph Peridynamic Theory and Its Applications. I'd love to work together with you. We'll first try to implement the thermal model. I'll let you know our process and there may be some beginner's questions in between. Our long-term plan is to do fluid-solid interaction by coupling CFD and PD codes.

[kaipartmann]

@vavrines, this sounds like a good plan! I recently talked with the developers of TrixiParticles.jl. Maybe their repo is also worth looking for you.

Using your example, I hacked something together. I am unsure if it works correctly - the code is probably full of bugs/errors. The structure of Peridynamics.jl is currently not suited for coupled problems, so I had to override many functions, which resulted in a lot of complicated internal code. It also only works with the current state of the main branch and is likely to break.

I am still in the process of reworking the internal structure of the package. I will consider how a coupled model could be easily implemented without changing many internal functions for the next version.

However, you can probably try this out for your model development. If you have questions, please do not hesitate to contact me.

Example simulation: heating of a bar

The best thing is that you start from an empty new project and add the current main branch of the package:

shell> mkdir pd_thermal_coupling

shell> cd pd_thermal_coupling

shell> julia

julia> ]

(@v1.10) pkg> activate .

(pd_thermal_coupling) pkg> add Peridynamics#main

Then you create the following file:

bb_thermal_coupling_model.jl

```julia using Peridynamics using Base.Threads struct BBTMaterial <: Peridynamics.AbstractMaterial end struct BBTPointParameters <: Peridynamics.AbstractPointParameters δ::Float64 rho::Float64 E::Float64 nu::Float64 G::Float64 K::Float64 λ::Float64 μ::Float64 Gc::Float64 εc::Float64 bc::Float64 kc::Float64 # thermal conductivity kp::Float64 # microconductivity aph::Float64 # thermal expansion cv::Float64 # specific heat capacity end @inline Peridynamics.point_param_type(::BBTMaterial) = BBTPointParameters @inline function Peridynamics.allowed_material_kwargs(::BBTMaterial) return (Peridynamics.DEFAULT_POINT_KWARGS..., :kc, :aph, :cv) end function Peridynamics.get_point_params(::BBTMaterial, p::Dict{Symbol,Any}) δ = Peridynamics.get_horizon(p) rho = Peridynamics.get_density(p) if haskey(p, :nu) msg = "keyword `nu` is not allowed for BBTMaterial!\n" msg *= "Bond-based peridynamics has a limitation on the possion ratio.\n" msg *= "Therefore, when using BBTMaterial, `nu` is hardcoded as 1/4.\n" throw(ArgumentError(msg)) else p[:nu] = 0.25 end E, nu, G, K, λ, μ = Peridynamics.get_elastic_params(p) Gc, εc = Peridynamics.get_frac_params(p, δ, K) bc = 18 * K / (π * δ^4) # bond constant kc, kp, aph, cv = get_thermal_params(p, δ) return BBTPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc, kc, kp, aph, cv) end function get_thermal_params(p::Dict{Symbol,Any}, δ) haskey(p, :kc) || throw(UndefKeywordError(:kc)) haskey(p, :aph) || throw(UndefKeywordError(:aph)) haskey(p, :cv) || throw(UndefKeywordError(:cv)) kc::Float64 = float(p[:kc]) kc ≤ 0 && throw(ArgumentError("`kc` should be larger than zero!\n")) aph::Float64 = float(p[:aph]) aph ≤ 0 && throw(ArgumentError("`aph` should be larger than zero!\n")) cv::Float64 = float(p[:cv]) cv ≤ 0 && throw(ArgumentError("`cv` should be larger than zero!\n")) kp = 6 * kc / (π * δ^4) return kc, kp, aph, cv end @inline Peridynamics.discretization_type(::BBTMaterial) = Peridynamics.BondDiscretization @inline function Peridynamics.init_discretization(body::Body{BBTMaterial}, args...) return Peridynamics.init_bond_discretization(body, args...) end struct BBTVerletStorage <: Peridynamics.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} temperature::Vector{Float64} pflux::Vector{Float64} end const BBTStorage = Union{BBTVerletStorage} @inline Peridynamics.storage_type(::BBTMaterial, ::VelocityVerlet) = BBTVerletStorage function Peridynamics.init_storage(::Body{BBTMaterial}, ::VelocityVerlet, bd::Peridynamics.BondDiscretization, ch::Peridynamics.ChunkHandler) n_loc_points = length(ch.loc_points) n_points_with_halo = length(ch.point_ids) position = copy(bd.position) displacement = zeros(3, n_loc_points) velocity = zeros(3, n_points_with_halo) 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) temperature = zeros(n_points_with_halo) pflux = zeros(n_loc_points) return BBTVerletStorage(position, displacement, velocity, velocity_half, acceleration, b_int, b_ext, damage, bond_active, n_active_bonds, temperature, pflux) end Peridynamics.storage_field(s::BBTStorage, ::Val{:temperature}) = getfield(s, :temperature) Peridynamics.storage_field(s::BBTStorage, ::Val{:pflux}) = getfield(s, :pflux) @inline function Peridynamics.get_halo_read_fields(s::BBTStorage) return (s.position, s.velocity, s.temperature) end @inline Peridynamics.get_halo_write_fields(::BBTStorage) = () const BBTBodyChunk = Union{Peridynamics.BodyChunk{BBTMaterial,P,D,S}, Peridynamics.MultiParamBodyChunk{BBTMaterial,P,D,S}} where {P,D,S} function calc_pflux!(b::BBTBodyChunk) b.store.pflux .= 0.0 for point_id in Peridynamics.each_point_idx(b) pflux_point!(b.store, b.dscr, b.mat, get_param(b, point_id), point_id) end return nothing end function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial, param::BBTPointParameters, i::Int) for bond_id in Peridynamics.each_bond_idx(bd, i) bond = bd.bonds[bond_id] j, L = bond.neighbor, 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] Δvijx = s.velocity[1, j] - s.velocity[1, i] Δvijy = s.velocity[2, j] - s.velocity[2, i] Δvijz = s.velocity[3, j] - s.velocity[3, i] l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz) ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l Δtem = s.temperature[j] - s.temperature[i] s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L - ev * 0.5 * param.bc * param.aph) * bd.volume[j] end return nothing end @inline function update_temperature!(b::Peridynamics.AbstractBodyChunk, Δt::Float64) for point_id in Peridynamics.each_point_idx(b) param = get_param(b, point_id) k = Δt / (param.rho * param.cv) _update_temperature!(b.store.temperature, b.store.pflux, k, point_id) end return nothing end @inline function _update_temperature!(temperature, pflux, k, i) temperature[i] += pflux[i] * k return nothing end function Peridynamics.force_density_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial, param::BBTPointParameters, i::Int) for bond_id in Peridynamics.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 εm = (l - L) / L εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph ε = εm + εt # failure mechanism if ε > param.εc && bond.fail_permit s.bond_active[bond_id] = false end s.n_active_bonds[i] += s.bond_active[bond_id] # update of force density temp = s.bond_active[bond_id] * 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 have to override the time stepping schemes function Peridynamics.solve_timestep!(dh::Peridynamics.ThreadsDataHandler{BBTBodyChunk}, options::Peridynamics.ExportOptions, Δt::Float64, Δt½::Float64, n::Int) t = n * Δt @threads :static for chunk_id in eachindex(dh.chunks) chunk = dh.chunks[chunk_id] Peridynamics.update_vel_half!(chunk, Δt½) Peridynamics.apply_bcs!(chunk, t) Peridynamics.update_disp_and_pos!(chunk, Δt) end @threads :static for chunk_id in eachindex(dh.chunks) Peridynamics.exchange_read_fields!(dh, chunk_id) chunk = dh.chunks[chunk_id] calc_pflux!(chunk) update_temperature!(chunk, Δt) end @threads :static for chunk_id in eachindex(dh.chunks) Peridynamics.exchange_read_fields!(dh, chunk_id) Peridynamics.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 function Peridynamics.solve_timestep!(dh::Peridynamics.MPIDataHandler{BBTBodyChunk}, options::Peridynamics.ExportOptions, Δt::Float64, Δt½::Float64, n::Int) t = n * Δt chunk = dh.chunk Peridynamics.@timeit_debug Peridynamics.TO "update_vel_half!" Peridynamics.update_vel_half!(chunk, Δt½) Peridynamics.@timeit_debug Peridynamics.TO "apply_bcs!" Peridynamics.apply_bcs!(chunk, t) Peridynamics.@timeit_debug Peridynamics.TO "update_disp_and_pos!" Peridynamics.update_disp_and_pos!(chunk, Δt) Peridynamics.@timeit_debug Peridynamics.TO "exchange_read_fields!" Peridynamics.exchange_read_fields!(dh) Peridynamics.@timeit_debug Peridynamics.TO "calc_pflux!" calc_pflux!(chunk) Peridynamics.@timeit_debug Peridynamics.TO "update_temperature!" update_temperature!(chunk, Δt) Peridynamics.@timeit_debug Peridynamics.TO "exchange_read_fields!" Peridynamics.exchange_read_fields!(dh) Peridynamics.@timeit_debug Peridynamics.TO "calc_force_density!" Peridynamics.calc_force_density!(chunk) Peridynamics.@timeit_debug Peridynamics.TO "exchange_write_fields!" Peridynamics.exchange_write_fields!(dh) Peridynamics.@timeit_debug Peridynamics.TO "calc_damage!" Peridynamics.calc_damage!(chunk) Peridynamics.@timeit_debug Peridynamics.TO "update_acc_and_vel!" Peridynamics.update_acc_and_vel!(chunk, Δt½) Peridynamics.@timeit_debug Peridynamics.TO "export_results" Peridynamics.export_results(dh, options, n, t) return nothing end function temperature_bc!(f::F, b::Body, name::Symbol) where {F<:Function} Peridynamics.check_if_set_is_defined(b.point_sets, name) type = Peridynamics.check_condition_function(f) if type === :sdbc sdbc = Peridynamics.SingleDimBC(f, :temperature, name, 0x1) Peridynamics._condition!(b.single_dim_bcs, sdbc) elseif type === :pdsdbc pdsdbc = Peridynamics.PosDepSingleDimBC(f, :temperature, name, 0x1) Peridynamics._condition!(b.posdep_single_dim_bcs, pdsdbc) end return nothing end function temperature_ic!(b::Body, name::Symbol, value::Real) Peridynamics.check_if_set_is_defined(b.point_sets, name) sdic = Peridynamics.SingleDimIC(convert(Float64, value), :temperature, name, 0x1) Peridynamics._condition!(b.single_dim_ics, sdic) return nothing end @inline function Peridynamics.apply_sdbc!(field::Vector{Float64}, value::Float64, ::UInt8, point_ids::Vector{Int}) @simd for i in point_ids @inbounds field[i] = value end return nothing end @inline function Peridynamics.apply_pdsdbc!(field::Vector{Float64}, pos::Matrix{Float64}, bc::Peridynamics.PosDepSingleDimBC{F}, point_ids::Vector{Int}, t::Float64) where {F} @simd for i in point_ids value = bc(SVector{3}(pos[1,i], pos[2,i], pos[3,i]), t) if !isnan(value) field[i] = value end end return nothing end function Peridynamics.apply_ic!(b::BBTBodyChunk, ic::Peridynamics.SingleDimIC) for point_id in b.psets[ic.point_set] _apply_ic!(Peridynamics.get_storage_field(b.store, ic.field), ic.value, ic.dim, point_id) end return nothing end function _apply_ic!(field::Matrix{Float64}, value::Float64, dim::UInt8, point_id::Int) setindex!(field, value, dim, point_id) return nothing end function _apply_ic!(field::Vector{Float64}, value::Float64, ::UInt8, point_id::Int) setindex!(field, value, point_id) return nothing end ```

Then you can run a simulation with the new API of the package:

include("bb_thermal_coupling_model.jl")

l, Δx, δ = 1.0, 1/100, 3.015/100
pos, vol = uniform_box(l, 0.1l, 0.1l, Δx)
body = Body(BBTMaterial(), pos, vol)
# add material parameters as in your code
material!(body; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7, kc=51.9, aph=11.5e-6, cv=472)
point_set!(x -> x > 0.0, body, :heating_area)
point_set!(x -> x < -l/2+3Δx, body, :fixed_area)
# temperature boundary condition: increase temperature by 1e5 K/s
temperature_bc!(t -> 1e5 * t, body, :heating_area)
velocity_bc!(t -> 0.0, body, :fixed_area, :x)
vv = VelocityVerlet(steps=2000)
job = Job(body, vv; path=joinpath(@__DIR__, "thermal_coupling"))
@time submit(job);

How it works

I mainly defined a bond-based thermal material and added new internal methods for this material type.

struct BBTMaterial <: Peridynamics.AbstractMaterial end

The pflux (thermal loop) is calculated in the pflux_point! function:

function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial,
                      param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, 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]
        Δvijx = s.velocity[1, j] - s.velocity[1, i]
        Δvijy = s.velocity[2, j] - s.velocity[2, i]
        Δvijz = s.velocity[3, j] - s.velocity[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
        ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l
        Δtem = s.temperature[j] - s.temperature[i]
        s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L -
                      ev * 0.5 * param.bc * param.aph) * bd.volume[j]
    end
    return nothing
end

Then, the force density calculation of the plain BBMaterial is extended to account for temperature changes:

function Peridynamics.force_density_point!(s::BBTStorage,
                                           bd::Peridynamics.BondDiscretization,
                                           ::BBTMaterial, param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.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
        εm = (l - L) / L
        εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph
        ε = εm + εt

        # failure mechanism
        if ε > param.εc && bond.fail_permit
            s.bond_active[bond_id] = false
        end
        s.n_active_bonds[i] += s.bond_active[bond_id]

        # update of force density
        temp = s.bond_active[bond_id] * 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

The other methods are mostly internal functions.

[kaipartmann]

@vavrines I will introduce changes to the main branch in a few minutes that will break the provided code. Therefore, I created a branch thermal_coupling that works as described with the example.

So in the local environment just do

(pd_thermal_coupling) pkg> rm Peridynamics

(pd_thermal_coupling) pkg> add Peridynamics#thermal_coupling

[ShiWeiHuH]

Thank you for your assistance. I am a member of vavrines's team and am currently working on this project. I am benefiting greatly from reviewing your code.

[ShiWeiHuH]

@vavrines I will introduce changes to the main branch in a few minutes that will break the provided code. Therefore, I created a branch thermal_coupling that works as described with the example.

So in the local environment just do

(pd_thermal_coupling) pkg> rm Peridynamics

(pd_thermal_coupling) pkg> add Peridynamics#thermal_coupling

My apologies for the interruption, may I ask you a few questions?

1. I encountered a bug when trying out the “bb_thermal_coupling_model.jl” script. It prompted an UndefVarError: discretization_type not defined at line 63 of bb_thermal_coupling_model.jl.
2. I would like to know how the surface and volume corrections are implemented in the code, specifically in which function.
3. During the thermal-structural coupling calculation, due to the presence of virtual layers, the domains for force and heat calculations should be different. Previously, we found that if this is not taken into account, there are significant calculation deviations, especially near the boundaries. In other words, during heat calculation, the displacement virtual layer needs to be excluded, and during force calculation, the temperature virtual layer needs to be excluded. Which specific function(s) need(s) to be modified for this setting? Thank you again for your assistance.

[kaipartmann]

Hello @ShiWeiHuH! I am happy that the example benefits you! :)

1. I encountered a bug when trying out the “bb_thermal_coupling_model.jl” script. It prompted an UndefVarError: discretization_type not defined at line 63 of bb_thermal_coupling_model.jl.

It seems that you're encountering an issue with the branch on GitHub. Please try in the package mode

pkg> status
[4dc47793] Peridynamics v0.2.0 `https://github.com/kaipartmann/Peridynamics.jl.git#thermal_coupling`

and check if it reads thermal_coupling at the end of the Peridynamics.jl entry.

2. I would like to know how the surface and volume corrections are implemented in the code, specifically in which function.

We currently do not have volume or surface corrections implemented. If you need that, I can add it. I tried to implement it before but did not finish to debug the code.

3. During the thermal-structural coupling calculation, due to the presence of virtual layers, the domains for force and heat calculations should be different. Previously, we found that if this is not taken into account, there are significant calculation deviations, especially near the boundaries. In other words, during heat calculation, the displacement virtual layer needs to be excluded, and during force calculation, the temperature virtual layer needs to be excluded. Which specific function(s) need(s) to be modified for this setting?

This is very interesting, and I did not know that. I am very new to thermal models with peridynamics, so I am unsure if I correctly understand the actual consequences of this regarding the code. So, is it a problem if a point has some mechanical and thermal boundary conditions simultaneously? Can you highlight your implementation of these effects in the script you provided for me?

[ShiWeiHuH]

Hi, @kaipartmann.Thank you very much for your response. As for question 1: I'm very sorry, the bug was caused by my own operational mistake. It has been resolved now. Thank you for your guidance. As for question 2: With a sufficient number of PD points, surface corrections and volume corrections can also yield good calculation results. I'm also trying to make some modifications to enhance my understanding of your code. Your code is very inspiring, and I'm learning from it.

3. Regarding the issue of thermal-mechanical coupling, I primarily addressed it in the init.jl file, where I also implemented surface corrections and volume corrections. I'll mark the relevant parts for you. Once again, thank you for your assistance.

[ShiWeiHuH]

Hi, @kaipartmann. We have updated the init.jl file. In the new init.jl, compared to the original version, the main content remains unchanged, but detailed comments have been added. We have annotated the codes for handling the thermal domain and the mechanical domain separately, while also marking the codes for surface correction and volume correction. During thermal coupling calculations, it's not always the case that the mechanical domain and the thermal domain are different. If there are displacement or velocity boundary conditions, as well as temperature boundary conditions, it is necessary to add a fictitious layer. Only when the fictitious layers for thermal and deformation solutions are located at different positions, should the scenario of different domains for motion and heat transfer calculations be considered. There are also cases where both temperature and displacement or velocity boundary conditions are set in the same fictitious layer, such as in thermal pulling. On the side receiving the pulling force, it's possible to simultaneously set temperature and displacement boundary conditions. However, for the fixed end on the other side, there might not be temperature boundary conditions, leading to a scenario where the thermal and mechanical domains are different.

[kaipartmann]

Hi @ShiWeiHuH, Thank you for your updates in the init.jl file. I think I can now better understand what the fictitious layers do. I hope to get some time in the next few days to implement these into an updated version of the bb_thermal_coupling_model.jl script. I am unsure how long it will take, but I will ping you if I encounter problems or have any questions.

[ShiWeiHuH]

Hi @ShiWeiHuH, Thank you for your updates in the init.jl file. I think I can now better understand what the fictitious layers do. I hope to get some time in the next few days to implement these into an updated version of the bb_thermal_coupling_model.jl script. I am unsure how long it will take, but I will ping you if I encounter problems or have any questions. Hello, @kaipartmann, I want to express my gratitude once again for your assistance. Over this period, I've been working on implementing surface and volume corrections based on your provided code, particularly focusing on the 'bond_discretization.jl' section. As a result, I've developed a 'modify.jl' file dedicated to this task. The basic functionality has been implemented. By executing the following three lines of code, we can derive the necessary correction coefficients. Here is the modified code :

include("bb_thermal_coupling_model.jl")
include("modify.jl")
Peridynamics.get_git_info() = "*"
l, Δx, δ = 1.0, 1/100, 3.015/100
pos, vol = uniform_box(l+6*Δx, 0.1l, 0.1l, Δx)
body = Body(BBTMaterial(), pos, vol)
# add material parameters as in your code
material!(body; horizon=3.015*Δx, E=2.1e5, rho=8e-6, Gc=2.7, kc=51.9, aph=11.5e-6, cv=472)
point_set!(x -> x < -l/2, body, :tem_bc)
point_set!(x -> x >  l/2, body, :fix_bc)
# temperature boundary condition: increase temperature by 1e5 K/s
temperature_bc!(t -> 1e5 * t, body, :tem_bc)
velocity_bc!(t -> 0.0, body, :fix_bc, :x)
mp = mpm_init(Dict(:dx => Δx, :horizon => 3.015*Δx, :emod => 2.1e5, :th_con => 51.9)) ### creat parameters for modify
sbd = mbd_init(body, mp) ### creat the index and volume modify coefficient
mm = undate_mbd(sbd, body, mp) ### creat surface coefficient
vv = VelocityVerlet(steps=2000)
job = Job(body, vv; path=joinpath(@__DIR__, "thermal_coupling"), fields=(:displacement, :damage))
@time submit(job)

This part is newly added：

mp = mpm_init(Dict(:dx => Δx, :horizon => 3.015*Δx, :emod => 2.1e5, :th_con => 51.9)) ### creat parameters for modify
sbd = mbd_init(body, mp) ### creat the index and volume modify coefficient
mm = undate_mbd(sbd, body, mp) ### creat surface coefficient

The modify.jl

using Peridynamics
using NearestNeighbors

mutable struct Mbond neighbor::Int length::Float64 vol_m::Float64 ### voleme modify coefficient surf_th::Float64 ### surface modify coefficient of thermal zone surf_mh::Float64 ### surface modify coefficient of deformation zone end

struct MbondDiscretization <: Peridynamics.AbstractDiscretization bonds::Vector{Mbond} n_neighbors::Vector{Int} bond_ids::Vector{UnitRange{Int}} volume::Vector{Float64} end

struct Mparameters dx::Float64 δ::Float64 kc::Float64 # thermal conductivity kp::Float64 # microconductivity E::Float64 bc::Float64 end

function mpm_init(p::Dict{Symbol,Float64}) dx = p[:dx] δ = p[:horizon] kc = p[:th_con] kp = 6 kc / (π δ^4) E = p[:emod] bc = 12 E / (π δ^4) mpm = Mparameters(dx, δ, kc, kp, E, bc) return mpm end

function mbd_init(body::Body, mpm::Mparameters) bonds, n_neighbors = mfind_bonds(body, mpm) bond_ids = mfind_bond_ids(n_neighbors) volume = body.volume mbd = MbondDiscretization(bonds, n_neighbors, bond_ids,volume) return mbd end

function mfind_bonds(body::Body, mpm::Mparameters) balltree = BallTree(body.position) bonds = Vector{Mbond}() sizehint!(bonds, body.n_points * 100) n_neighbors = zeros(Int, body.n_points) for i in 1:body.n_points n_neighbors[i] = mfind_bonds!(bonds, balltree, body.position, mpm, i) end return bonds, n_neighbors end

function mfind_bonds!(bonds::Vector{Mbond}, balltree::BallTree, position::Matrix{Float64}, mpm::Mparameters, i::Int)

neigh_idxs_with_i = inrange(balltree, view(position, :, i), mpm.δ)
neigh_idxs = filter(x -> x != i, neigh_idxs_with_i)
for j in neigh_idxs
    L = sqrt((position[1, j] - position[1, i])^2 +
             (position[2, j] - position[2, i])^2 +
             (position[3, j] - position[3, i])^2)

             if L < (mpm.δ - 0.5 * mpm.dx) 
                vol_m = 1.0
            elseif L < (mpm.δ + 0.5 * mpm.dx)  
                vol_m  = (mpm.δ + 0.5 * mpm.dx - L) / mpm.dx
            else
                vol_m = 0.0
            end

            surf_th = 1.0
            surf_mh = 1.0
    push!(bonds, Mbond(j, L, vol_m, surf_th, surf_mh))
end
return length(neigh_idxs)

end

function mfind_bond_ids(n_neighbors::Vector{Int}) bond_ids = fill(0:0, length(n_neighbors)) bonds_start, bonds_end = 1, 0 for i in eachindex(n_neighbors) bonds_end = bonds_start + n_neighbors[i] - 1 bond_ids[i] = bonds_start:bonds_end bonds_start += n_neighbors[i] end return bond_ids end

function undate_mbd(mbd::MbondDiscretization, bo::Body, mpm::Mparameters)

tem = 0.001.*[sum(bo.position[:, i]) for i in 1:bo.n_points]
disx = copy(bo.position)
disx[1, :] .= 1.001 .* disx[1, :]
disy = copy(bo.position)
disy[2, :] .= 1.001 .* disy[2, :]
disz = copy(bo.position)
disz[3, :] .= 1.001 .* disz[3, :]

point_th = zeros(Float64, bo.n_points)
point_mx = zeros(Float64, bo.n_points)
point_my = zeros(Float64, bo.n_points)
point_mz = zeros(Float64, bo.n_points)

for i in 1:bo.n_points
    pd_th = 0
    pd_mechx = 0
    pd_mechy = 0
    pd_mechz = 0

    for bond_id in meach_bond_idx(mbd, i)
        bond = mbd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        if i in bo.point_sets[:fix_bc] || j in bo.point_sets[:fix_bc]
            mfth = 0
        else
            mfth = 1
        end

        if i in bo.point_sets[:tem_bc] || j in bo.point_sets[:tem_bc]
            mfmech = 0
        else
            mfmech = 1
        end

        nlengthx = sqrt((disx[1,j] - disx[1,i])^2 + (disx[2,j] - disx[2,i])^2 + (disx[3,j] - disx[3,i])^2)
        nlengthy = sqrt((disy[1,j] - disy[1,i])^2 + (disy[2,j] - disy[2,i])^2 + (disy[3,j] - disy[3,i])^2)
        nlengthz = sqrt((disz[1,j] - disz[1,i])^2 + (disz[2,j] - disz[2,i])^2 + (disz[3,j] - disz[3,i])^2)

        pd_th += 0.25 * mpm.kp * (tem[j] - tem[i])^2 / L * mbd.volume[j] * bond.vol_m * mfth
        pd_mechx += 0.25 * mpm.bc * ((nlengthx-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
        pd_mechy += 0.25 * mpm.bc * ((nlengthy-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
        pd_mechz += 0.25 * mpm.bc * ((nlengthz-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
    end 

    if pd_th != 0
        point_th[i,1] = 1e-6 *  mpm.kc * 1.5 / pd_th
    else
        point_th[i,1] = 0
    end

    if pd_mechx != 0
        point_mx[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechx
    else
        point_mx[i,1] = 0
    end

    if pd_mechy != 0
        point_my[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechy
    else
        point_my[i,1] = 0
    end

    if pd_mechz != 0
        point_mz[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechz
    else
        point_mz[i,1] = 0
    end
end

for i in 1:bo.n_points
    for bond_id in meach_bond_idx(mbd, i)
        bond = mbd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        bond.surf_th = 0.5 * (point_th[i]+point_th[j])

        if abs(bo.position[3,j] - bo.position[3,i]) < 1e-10
            if abs(bo.position[2,j] - bo.position[2,i]) < 1e-10
                θ = 0
            elseif abs(bo.position[1,j] - bo.position[1,i]) < 1e-10
                θ = π/2
            else
                θ = atan(abs((bo.position[2,j] - bo.position[2,i])/(bo.position[1,j] - bo.position[1,i])))
            end
            ϕ = π/2
            sx = (point_mx[i] + point_mx[j]) / 2
            sy = (point_my[i] + point_my[j]) / 2
            sz = (point_mz[i] + point_mz[j]) / 2
            bond.surf_mh = sqrt(1/((cos(θ)*sin(ϕ)/sx)^2 + (sin(θ)*sin(ϕ)/sy)^2 + (cos(ϕ)/sz)^2))
        elseif (abs(bo.position[2,j] - bo.position[2,i]) < 1e-10) && (abs(bo.position[1,j] - bo.position[1,i]) < 1e-10)
            sz = (point_mz[i] + point_mz[j]) / 2
            bond.surf_mh = sz
        else
            θ = atan(abs((bo.position[2,j] - bo.position[2,i])/(bo.position[1,j] - bo.position[1,i])))
            ϕ = acos((bo.position[3,j] - bo.position[3,i])/L)
            sx = (point_mx[i] + point_mx[j]) / 2
            sy = (point_my[i] + point_my[j]) / 2
            sz = (point_mz[i] + point_mz[j]) / 2
            bond.surf_mh = sqrt(1/((cos(θ)*sin(ϕ)/sx)^2 + (sin(θ)*sin(ϕ)/sy)^2 + (cos(ϕ)/sz)^2))
        end
    end  
end  
return mbd        

end

@inline meach_bond_idx(mbd::MbondDiscretization, point_id::Int) = mbd.bond_ids[point_id]

My next step is to pass the correction coefficients into 'bb_thermal_coupling_model.jl'. 
Additionally, you mentioned that the time stepping schemes need to be rewritten. What are the main points to pay attention to in this part? I will try to proceed with it next.

[kaipartmann]

Hi @ShiWeiHuH,

Thank you very much for providing your code!

First, I'm sorry about the bug with the get_git_info() function. It has already been solved in the main branch.

I ran your code, and it should work fine. You are correct in bridging the parallel computation features with your newly defined MBondDiscretization, because they cause problems when calculating the surface correction factors.

I encountered the same issue and am working on it (#73). I'm now prioritizing the implementation of the surface correction (see #74), and your assistance is greatly appreciated.

I already tried to do something similar to what you did, but my code still has some bugs. You can take a look at the surface_correction branch in the file: https://github.com/kaipartmann/Peridynamics.jl/blob/surface_correction/src/physics/bond_based_corrected.jl These lines look very similar to your code: https://github.com/kaipartmann/Peridynamics.jl/blob/1cb431789f94db74e4bd545b3e2284334054656c/src/physics/bond_based_corrected.jl#L205-L243

You can proceed by adding your MBondDiscretization as a discretization method with

@inline function Peridynamics.init_discretization(body::Body{BBTMaterial}, args...)
    return mbd_init(body, ...)
end

But your mbd_init function needs to return a ChunkHandler and has to be similar to: https://github.com/kaipartmann/Peridynamics.jl/blob/a71541fc6ef511b92031305f95f42c922dfd7aee/src/discretizations/bond_discretization.jl#L15-L28

However, please note that in the main branch, this process is already more straightforward when using the @storage macro. It might be easier to proceed by copying some of the code from there.

If I am ready with #73, I can provide a new script of the current state that includes your modify.jl with the thermal computation.

[vavrines]

Hi @kaipartmann ,

Thanks for the tips! Besides of #73 and #74, please let us know when the main branch reaches a relatively stable time. It might not be a bad idea to start pulling things before we go too deep in the obsolete codes.

[ShiWeiHuH]

@kaipartmann Thank you for your reply. After reviewing this portion of surface_correction, It effectively modifies surface effects. 💯 🥇 👍 This approach reminds me of the method used in the 《Peridynamic Theory and Its Applications》. Volume and surface corrections are computed as coefficients for bond force calculation in each iteration. Even though the code in this book calculates correction coefficients with each iteration when computing forces, the calculations still rely on the coordinates of the initial PD points. Correction coefficients are a physical property that reflects the initial geometric characteristics, and once the geometric structure is generated, these coefficients remain fixed and do not change. Some subsequent articles have suggested that parameters such as correction coefficients should be precomputed, and then directly called in subsequent iterations, which can significantly reduce computational overhead. Inspired by the part of your code related to bond indexing and length handling, I attempted to precompute correction coefficients and then call them in each iteration. Regarding the issue of different thermal dynamics in different regions, it can be addressed during surface correction. For example, in thermal dynamics computations, if any PD point of a bond lies within the virtual layer of displacement constraints, the correction coefficient is set to 0, thus eliminating the influence of the virtual layer. I will proceed with experimentation as you suggested. Thanks again for your assistance. 👍

[kaipartmann]

@ShiWeiHuH thank you for taking the time to review the code. Your assistance is greatly appreciated!

I wanted to do what you suggested but needed #73 solved for that. A halo exchange in the initialization phase is needed to calculate the bond correction factors, and I previously designed the exchange functions so that they would only work during the time integration. Therefore, I took this quick and dirty approach to calculate the factors during each time step because then it would be possible to exchange the values to the halo points.

I will try to add the surface correction as part of the BondSystem. If I am ready, you can probably copy the code in bond_system.jl and create a ThermalBondSystem that handles thermal correction factors.

Unfortunately, the surface correction still has some bugs; I will explain this closer in #74 and ping you there.

[ShiWeiHuH]

@kaipartmann Thank you. I understand what you mean now. My current approach is somewhat makeshift. Considering that the calculation of correction factors only occurs once at the beginning, I bypassed the issue of parallelization and simply calculated the correction factors separately and left them as is. It would be great if we could handle them uniformly within your parallel framework. Today, I managed to integrate the computed correction factors into the coupling program via function calls. Similarly, it's not the most elegant solution, haha. Below is the code I modified. Take a look when you have time. Lastly, looking forward to your new code. Many thanks, as always. I made just minor adjustments to bb_thermal_coupling_model.jl. I just added a function call in the heat flux and force calculations to incorporate the correction factors calculated in th.jl and modify.jl into the computations. (th.jl and modify.jl are both the same as yesterday, no changes.） thermal flux

function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial,
                      param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

############## introduced datas into the calculation #####newly added
        mbond = get_init_mbd(mm, bond_id) 
        volm = mbond.vol_m
        mof_th = mbond.surf_th
##############

        Δ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]
        Δvijx = s.velocity[1, j] - s.velocity[1, i]
        Δvijy = s.velocity[2, j] - s.velocity[2, i]
        Δvijz = s.velocity[3, j] - s.velocity[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
        ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l
        Δtem = s.temperature[j] - s.temperature[i]
        s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L -
                      ev * 0.5 * param.bc * param.aph) * bd.volume[j] * volm * mof_th   ### modifys coefficient
    end
    return nothing
end

force density

function Peridynamics.force_density_point!(s::BBTStorage,
                                           bd::Peridynamics.BondDiscretization,
                                           ::BBTMaterial, param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

######### introduced datas into the calculation#####newly added
        mbond = get_init_mbd(mm, bond_id)
        volm = mbond.vol_m
        mof_mh = mbond.surf_mh
#########

        # 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
        εm = (l - L) / L
        εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph
        ε = εm + εt

        # failure mechanism
        if ε > param.εc && bond.fail_permit
            s.bond_active[bond_id] = false
        end
        s.n_active_bonds[i] += s.bond_active[bond_id]

        # update of force density
        temp = s.bond_active[bond_id] * param.bc * ε / l * bd.volume[j] * volm * mof_mh  ### modifys coefficient
        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

Function

#####newly added
@inline get_init_mbd(modifys_coefficients::MbondDiscretization, point_id::Int64) = modifys_coefficients.bonds[point_id]

[kaipartmann]

@ShiWeiHuH @vavrines I have a question regarding your surface correction factors. In your code, you used:

0.6 * 1e-6 * mpm.E

How did you derive these analytical strain energy density? The surface correction works correctly with this value, but I would like to better understand the assumptions behind it.

[ShiWeiHuH]

Hi, @kaipartmann You can take a look at this simple derivation. In addition, I would like to ask you what is the command to execute multi-threading or multi-node parallelism. Also, I see that you have corrected the surface effect, and it seems that the volume has not been corrected. After all, the volume at the boundary of the near-field domain is not necessarily complete.

[kaipartmann]

Hi @ShiWeiHuH, Thank you very much for your detailed description! I added an issue (#94) because the OSBMaterial allows for the use of arbitrary Poisson's ratios.

On the parallelism:

The parallelism functionality is implemented at a very high level. With the PointDecomposition, the Body is chopped into multiple BodyChunks, where each chunk contains all relevant information for the simulation. Depending on the type of parallelism, either a ThreadsDataHandler or a MPIDataHandler is created that contain all the body chunks and all information for the simulation.

Then, everything is done on the body chunk level. For multithreading, each thread works on its chunk: https://github.com/kaipartmann/Peridynamics.jl/blob/9bb65ae8e467532ff53e34be67bd66e2434feb86/src/time_solvers/velocity_verlet.jl#L212-L233 For MPI, each process uses its chunk, and then the processes need to communicate to exchange the information of points at the chunk border (halo points): https://github.com/kaipartmann/Peridynamics.jl/blob/9bb65ae8e467532ff53e34be67bd66e2434feb86/src/time_solvers/velocity_verlet.jl#L250-L264

If you have further questions regarding the parallelism, please do not hesitate to contact me.

On volume correction:

I did not include the volume correction because it needs a uniformly distributed point cloud to be correct, which is not guaranteed with the current approach of our package. If someone uses an arbitrary mesh with Abaqus and converts that into a point cloud, the volume correction will introduce errors, and the point spacing is unknown.

However, from our observations, the correspondence formulation is also very precise with arbitrary meshes, and the volume error had only a very negligible effect on the simulation results, so we did not prioritize implementing it.

Do you have observations on how much the volume error influences the simulation results in your thermo-mechanical coupled simulations? If you use the correspondence formulation instead of the bond-based formulation, all the corrections could theoretically be dismissed.

However, I have some ideas about how to guarantee that a point cloud is uniformly distributed with a given point spacing. This would again break the current API but make it possible to use volume correction if a point cloud is guaranteed to be uniform and to ignore it otherwise. I will open an issue and ping you if my idea is clearer.

[ShiWeiHuH]

Hi, @kaipartmann It's been a while!

I noticed that PD has a new version with many great features. That's fantastic!

I have a couple of questions:

Currently, I'm using single-thread computation with one region. If I want to switch to using 64 threads, should I set the environment variable NUM=64 first? Will this automatically divide the region into 64 blocks?

Can computation be resumed if it shuts down unexpectedly?

[kaipartmann]

Hi @ShiWeiHuH, I am very happy you like the new features and changes we made in the last weeks!

Threads

To use multiple threads you can use a few different approaches, as also mentioned [here].(https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads)

Use Julia command line options

julia -t 64
julia --threads 64

Specify the number of threads in the VSCode settings If you use VSCode, then you can specify the number of threads for a Julia REPL in the settings.

Setting an environment variable

JULIA_NUM_THREADS=64

You can check how many threads you use with

Threads.nthreads()

The body is then divided in different blocks with distribute_equally https://github.com/kaipartmann/Peridynamics.jl/blob/bc81309d3057533e2c972a7b356d9f3be762a371/src/discretization/decomposition.jl#L18-L44

You can check the different regions with

body = ...
n_chunks = 64
point_decomp = Peridynamics.PointDecomposition(body, n_chunks)
point_decomp.decomp # Vector with n_chunks elements and a range of the point ids for each region

Resume simulations

For now, it is not possible to resume simulations. The current approach is to capture as many bugs as possible before submitting the simulation job. If a simulation stops unexpectedly, all exported data up to this point is preserved.

However, the main reason I have not considered it for now is performance. It would be necessary to export the current state of the data handler regularly, which would then take a significant fraction of the simulation time (I estimate up to 50% for some setups).

I made an issue https://github.com/kaipartmann/Peridynamics.jl/issues/144, so maybe there is a way to do it efficiently, and it will be included in a future version!