Energy based surface correction

[kaipartmann]

Add the energy based surface correction algorithm for bond-based, ordinary state-based and the continuum-kinematics-inspired peridynamics formulation, as described in [1].

We are asked to use surface corrections for a revision of a paper, so this is a high priority now.

[1] Le, Q.V., Bobaru, F. Surface corrections for peridynamic models in elasticity and fracture. Comput Mech 61, 499–518 (2018). https://doi.org/10.1007/s00466-017-1469-1

[kaipartmann]

Already thought about this in #64. Will think about volume correction later.

[kaipartmann]

@ShiWeiHuH, the surfcorr#74 branch now contains the code with completed surface correction.

I added an abstract type AbstractCorrection which allows to implement own correction types and add these to the BondSystem. In order to do that, you can look at the file https://github.com/kaipartmann/Peridynamics.jl/blob/surfcorr%2374/src/discretization/bond_system_corrections.jl Maybe I find the time in the next days to add your recent code in https://github.com/kaipartmann/Peridynamics.jl/issues/32#issuecomment-2063321965 with this new feature of the BondSystem.

I checked that everything is type-stable and fast with Cthulhu. To achieve that, users now just have to add a type parameter to the material, if they want to opt in to use surface correction:

BBMaterial{EnergySurfaceCorrection}() # or
OSBMaterial{EnergySurfaceCorrection}()

[kaipartmann]

I want to validate the results first, before adding them to the main branch.

[ShiWeiHuH]

@kaipartmann Thank you very much for the new code. It works great and successfully corrected the surface effects of the forces. Next, I will try the surface correction for thermal effects according to the new code. I will continue to work on coupling issues and conduct tests. Thanks again for your help.

[kaipartmann]

Validation of the surface correction factors

@ShiWeiHuH I validated the energy based surface correction factors with a simple example:

function xwave(mat, lx, lyz, npyz, T, vmax, path)
    Δx = lyz / npyz
    pos, vol = uniform_box(lx, lyz, lyz, Δx)
    ids = sortperm(pos[1,:])
    body = Body(mat, pos[:, ids], vol[ids])
    failure_permit!(body, false)
    material!(body, horizon=3.015Δx, rho=7850.0, E=210e9, epsilon_c=0.01)
    point_set!(x -> x < -lx / 2 + 1.2Δx, body, :left)
    velocity_bc!(t -> t < T ? vmax * sin(2π / T * t) : 0.0, body, :left, :x)
    vv = VelocityVerlet(steps=2000)
    job = Job(body, vv; path=path)
    return job
end

Parameters:

lx, lyz, npyz = 0.2, 0.002, 5
T, vmax = 1.0e-5, 2.0

In this calculation, a displacement wave travels through the bar. The wave speed is then compared to the analytically calculated value.

Bond based no correction:

mat = BBMaterial() # defaults to BBMaterial{NoCorrection}()

--- WAVE SPEED SUMMARY ---
Theoretical wave speed c_0 = √(E/ρ):
  c_0 =  5172.1941530 m/s
Calculated wave speed in simulation:
  c_w =  4405.4743475 m/s
  Δc  =  -766.7198055 m/s (-14.824 %)

Bond based with correction:

mat = BBMaterial{EnergySurfaceCorrection}()

--- WAVE SPEED SUMMARY ---
Theoretical wave speed c_0 = √(E/ρ):
  c_0 =  5172.1941530 m/s
Calculated wave speed in simulation:
  c_w =  5312.4271502 m/s
  Δc  =   140.2329972 m/s (2.711 %)

Ordinary state-based without correction:

Note, that for the OSBMaterial, the keyword nu=0.25 was added to the material!(body, ...) call.

mat = OSBMaterial() # defaults to OSBMaterial{NoCorrection}()

--- WAVE SPEED SUMMARY ---
Theoretical wave speed c_0 = √(E/ρ):
  c_0 =  5172.1941530 m/s
Calculated wave speed in simulation:
  c_w =  5630.0295537 m/s
  Δc  =   457.8354006 m/s (8.852 %)

Ordinary state-based with correction:

mat = OSBMaterial{EnergySurfaceCorrection}()

--- WAVE SPEED SUMMARY ---
Theoretical wave speed c_0 = √(E/ρ):
  c_0 =  5172.1941530 m/s
Calculated wave speed in simulation:
  c_w =  5362.4009210 m/s
  Δc  =   190.2067680 m/s (3.677 %)

Additionally, the correspondence formulation that has no surface effect:

mat = NOSBMaterial()

--- WAVE SPEED SUMMARY ---
Theoretical wave speed c_0 = √(E/ρ):
  c_0 =  5172.1941530 m/s
Calculated wave speed in simulation:
  c_w =  5166.3805166 m/s
  Δc  =    -5.8136364 m/s (-0.112 %)