Wave propagation across material interface

Based on the wave propagation tutorial, this tutorial features a displacement wave crossing a material interface, as investigated in [PDW24]. This means that a bar containing two sections with different material properties is regarded.

First import the Peridynamics.jl package:

using Peridynamics

Then the geometric parameters are set, which are the same as in the wave propagation tutorial. These are the length lx, width and height lyz of the cuboid as well as the number of points in the width npyz, which determines the point spacing Δx.

lx = 0.2
lyz = 0.002
npyz = 4
Δx = lyz / npyz

0.0005

With these parameters we now create a body, here using the non-ordinary state-based correspondence formulation.

pos, vol = uniform_box(lx, lyz, lyz, Δx)
body = Body(NOSBMaterial(), pos, vol)

6400-point Body{NOSBMaterial}:
  1 point set(s):
    6400-point set `all_points`

Again, failure is not allowed in the whole body.

failure_permit!(body, false)

Then the material parameters for one half of the body are assigned to the whole body first.

material!(body, horizon=3.015Δx, rho=7850.0, E=210e9, nu=0.25, epsilon_c=0.01)

Now a point set containing the other half of all points is created.

point_set!(x -> x < 0, body, :set1)

The parameters for this point set are then overwritten with their new parameters.

material!(body, :set1, horizon=3.015Δx, rho=7850.0, E=105e9, nu=0.25, epsilon_c=0.01)

Except for the Young's modulus, these are the same in both sections:

To employ the boundary conditions creating a displacement wave, the point set :left is created:

point_set!(x -> x < -lx / 2 + 1.2Δx, body, :left)

As in the wave propagation tutorial, the applied velocity boundary condition is

\[ {v}_x (t) = \begin{cases} v_\mathrm{max} \cdot \sin(2\pi \cdot \frac{t}{T}) \qquad & \forall \; 0 \leq t \leq T \\ 0 &\text{else.} \end{cases}\]

T, vmax = 1.0e-5, 2.0
velocity_bc!(t -> t < T ? vmax * sin(2π / T * t) : 0.0, body, :left, :x)

The Velocity Verlet algorithm is used as time integration method and 2000 time steps are calculated:

vv = VelocityVerlet(steps=2000)

VelocityVerlet:
  n_steps        2000
  safety_factor  0.7

Finally the job is defined and submitted.

job = Job(body, vv; path="results/xwave_interface")
submit(job)

The following video shows the results for the displacement in x-direction $u_x$ and the first entry of the stress tensor, i.e. the normal stress in x-direction $\sigma_{11}$.