Tutorial Kalthoff-Winkler

[kaipartmann]

Add a tutorial for the Kalthoff-Winkler experiment.

[oldninja]

I would like to help with this one. What would be the process to add the tutorial? Should I just follow https://github.com/kaipartmann/Peridynamics.jl/blob/main/docs/src/literate/tutorial_tension_dynfrac.jl?

[oldninja]

I started putting something together, but it is not working right now. I tried to update to the latest version, 0.3.0-DEV, but unfortunately, I could not get it to work, so I am using the previous version of 0.2.0 format.

# Newton, meter, Pascal
using Pkg
using Peridynamics
len = 200.0E-3
wid = 100.0E-3
thi =   9.0E-3
Δx  =   1.0E-3
δ   = 4.015Δx
a   =  50.0E-3
pc = PointCloud(len, wid, thi, Δx)
mat= BondBasedMaterial(; horizon=4.015Δx, E=191.0e+9, rho=8000.0, Gc=22120.4)
set_a = findall(p -> -a/2+δ ≤ p[1] ≤ -a/2+δ && 0 ≤ p[2] < wid/2, eachcol(pc.position))
set_b = findall(p -> -a/2-δ ≤ p[1] ≤ -a/2-δ && 0 ≤ p[2] < wid/2, eachcol(pc.position))
precrack1 = PreCrack(set_a, set_b)
set_c = findall(p ->  a/2+δ ≤ p[1] ≤  a/2+δ && 0 ≤ p[2] < wid/2, eachcol(pc.position))
set_d = findall(p ->  a/2-δ ≤ p[1] ≤  a/2-δ && 0 ≤ p[2] < wid/2, eachcol(pc.position))
precrack2 = PreCrack(set_c, set_d)
set_top  = findall(p -> -a/2+2.015δ ≤ p[1] ≤  a/2-2.015δ && p[2] > wid/2-2.015Δx, eachcol(pc.position))
set_bot  = findall(p -> -δ ≤ p[1] ≤ δ && p[2] < -wid/2+2.015Δx, eachcol(pc.position))
bc_top    = VelocityBC(t ->  -32., set_top, 2)
bc_bot_x  = VelocityBC(t ->  0.0, set_bot, 1)
bc_bot_y  = VelocityBC(t ->  0.0, set_bot, 2)
vv = VelocityVerlet(2000, 0.95)
jobname = "KW"
path = joinpath("results", jobname)
if ispath(path)
    rm(path; recursive=true) # Delete existing results
end
!ispath(path) && mkpath(path) # create the path if it does not exist
es = ExportSettings(path, 100)
job = PDSingleBodyAnalysis(name=jobname, pc=pc, mat=mat, precracks=[precrack1,precrack2], bcs=[bc_top,bc_bot_x,bc_bot_y], td=vv, es=es)
submit(job)

[kaipartmann]

Great! I really appreciate your help!

The tutorial you mentioned can be a good starting point, so using Literate would be preferred. Just ensure that the submit(job) command will be shown as a plain code block that does not run during documentation generation.

That code looks great! Version v0.3.0 will probably be ready next week, so I may try to convert it to the current state of the main branch.

[kaipartmann]

I just tried your code and the model looks good. However, I get NaN values for every field after the first time step.

Converting to v0.3.0-DEV solved it:

# Newton, meter, Pascal
using Peridynamics
len = 200.0E-3
wid = 100.0E-3
thi =   9.0E-3
Δx  =   1.0E-3
δ   = 4.015Δx
a   =  50.0E-3
pos, vol = uniform_box(len, wid, thi, Δx)
body = Body(BBMaterial(), pos, vol)
material!(body; horizon=4.015Δx, E=191.0e+9, rho=8000.0, Gc=22120.4)
point_set!(p -> -a/2+δ ≤ p[1] ≤ -a/2+δ && 0 ≤ p[2] < wid/2, body, :set_a)
point_set!(p -> -a/2-δ ≤ p[1] ≤ -a/2-δ && 0 ≤ p[2] < wid/2, body, :set_b)
precrack!(body, :set_a, :set_b)
point_set!(p ->  a/2+δ ≤ p[1] ≤  a/2+δ && 0 ≤ p[2] < wid/2, body, :set_c)
point_set!(p ->  a/2-δ ≤ p[1] ≤  a/2-δ && 0 ≤ p[2] < wid/2, body, :set_d)
precrack!(body, :set_c, :set_d)
point_set!(p -> -a/2+2.015δ ≤ p[1] ≤  a/2-2.015δ && p[2] > wid/2-2.015Δx, body, :set_top)
point_set!(p -> -δ ≤ p[1] ≤ δ && p[2] < -wid/2+2.015Δx, body, :set_bot)
velocity_bc!(t -> -32., body, :set_top, :y)
velocity_bc!(t -> 0.0, body, :set_bot, :x)
velocity_bc!(t -> 0.0, body, :set_bot, :y)
vv = VelocityVerlet(steps=2000, safety_factor=0.95)
path = joinpath("results", "KW")
if ispath(path)
    rm(path; recursive=true) # Delete existing results
end
!ispath(path) && mkpath(path) # create the path if it does not exist
job = Job(body, vv; path=path, freq=100)
submit(job)

                                                               _
           _____         _     _                             _(_)_
          | ___ \       (_)   | |                           (_) (_)
          | |_/ /__ _ __ _  __| |_   _ _ __   __ _ _ __ ___  _  ___ ___
          |  __/ _ \ '__| |/ _` | | | | '_ \ / _` | '_ ` _ \| |/ __/ __|
          | | |  __/ |  | | (_| | |_| | | | | (_| | | | | | | | (__\__ \
          \_|  \___|_|  |_|\__,_|\__, |_| |_|\__,_|_| |_| |_|_|\___|___/
                                  __/ |
                                 |___/   Copyright (c) 2024 Kai Partmann

MULTITHREADING SIMULATION WITH 6 THREADS
BODY
  POINT CLOUD
    number of points ...................................................... 180000
    min, max values x-direction .................................. -0.0995, 0.0995
    min, max values y-direction .................................. -0.0495, 0.0495
    min, max values z-direction .................................... -0.004, 0.004
  POINT SETS
    number of points in set `set_bot` ........................................ 144
    number of points in set `all_points` .................................. 180000
    number of points in set `set_top` ........................................ 612
    number of points in set `set_c` ............................................ 0
    number of points in set `set_d` ............................................ 0
    number of points in set `set_a` ............................................ 0
    number of points in set `set_b` ............................................ 0
  CONDITIONS
    number of BC's ............................................................. 3
  MATERIAL
    material type ....................................... BBMaterial{NoCorrection}
    horizon ............................................................. 0.004015
    density ................................................................. 8000
    Young's modulus ..................................................... 1.91e+11
    Poisson's ratio ......................................................... 0.25
    shear modulus ....................................................... 7.64e+10
    bulk modulus .................................................... 1.273333e+11
DATA HANDLER CREATION COMPLETED ✓
BOND SYSTEM
  number of bonds ................................................... 3.778458e+07
VELOCITY VERLET TIME SOLVER
  number of time steps ...................................................... 2000
  time step size .................................................... 2.326834e-07
  time step safety factor ................................................... 0.95
  simulation time ................................................... 0.0004653667
TIME INTEGRATION LOOP 100%|████████████████████████████████████████| Time: 0:01:09

[kaipartmann]

There seems to be an issue with the point sets for the predefined cracks, as they are empty.

Also if you have any problems using v0.3.0-DEV please contact me, I would be very happy to help you!

[kaipartmann]

I think this resolved the issue with the point sets. I also added the usage of the energy based surface correction, which does not really have a large overhead. Also

!ispath(path) && mkpath(path) # create the path if it does not exist

is no longer needed, as the path gets created during simulation :)

# Newton, meter, Pascal
using Peridynamics

function kw(len, wid, thi, Δx, δ, a, path)
    pos, vol = uniform_box(len, wid, thi, Δx)
    body = Body(BBMaterial{EnergySurfaceCorrection}(), pos, vol)
    material!(body; horizon=4.015Δx, E=191.0e+9, rho=8000.0, Gc=22120.4)
    point_set!(p -> -a/2-δ ≤ p[1] ≤ -a/2 && 0 ≤ p[2] < wid/2, body, :set_a)
    point_set!(p -> -a/2 ≤ p[1] ≤ -a/2+δ && 0 ≤ p[2] < wid/2, body, :set_b)
    precrack!(body, :set_a, :set_b)
    point_set!(p ->  a/2-δ ≤ p[1] ≤  a/2 && 0 ≤ p[2] < wid/2, body, :set_c)
    point_set!(p ->  a/2 ≤ p[1] ≤  a/2+δ && 0 ≤ p[2] < wid/2, body, :set_d)
    precrack!(body, :set_c, :set_d)
    point_set!(p -> -a/2+2.015δ ≤ p[1] ≤  a/2-2.015δ && p[2] > wid/2-2.015Δx, body, :set_top)
    point_set!(p -> -δ ≤ p[1] ≤ δ && p[2] < -wid/2+2.015Δx, body, :set_bot)
    velocity_bc!(t -> -32., body, :set_top, :y)
    velocity_bc!(t -> 0.0, body, :set_bot, :x)
    velocity_bc!(t -> 0.0, body, :set_bot, :y)
    vv = VelocityVerlet(steps=2000, safety_factor=0.95)
    job = Job(body, vv; path=path, freq=100)
    @mpitime submit(job)
    return nothing
end

len = 200.0E-3
wid = 100.0E-3
thi =   9.0E-3
Δx  =   1.0E-3
δ   = 4.015Δx
a   =  50.0E-3
path = joinpath("results", "KW")
ispath(path) && rm(path; recursive=true) # Delete existing results
kw(len, wid, thi, Δx, δ, a, path)

[kaipartmann]

@oldninja Please feel free to open a PR regarding the tutorial at any time.

Also, I would like to know more about the issues you ran into when trying v0.3.0-DEV, as they could potentially be bugs that are unknown to me.