wave propagation

In this example, we demonstrate an example of wave propagation in a 2D solid.

Essentially, it is a 2D dynamics problem. The CP4 element is used to model the solid. Various time integration methods will be used to compare their performance regarding numerical energy dispersion.

The model script can be downloaded here.

This model contains 16641 nodes and 16384 elements. The memory usage is about 1.2 GB.

Due to small step size and need to export visualization data, the full analysis takes around 10 minutes to complete on an average PC platform (6 physical cores @ 3 GHz).

Model

A square solid of size $3200\times3200$ is used. The structured mesh can be generated by using whatever mesh generator available. It is not difficult to generate an array of squares using scripting languages such as Python or Matlab. Here, Gmsh is used.

The left boundary is constrained along the horizontal direction. The bottom boundary is constrained along the vertical direction. The top and right boundaries are free. A force load is applied to the top right corner.

Material

Whether plane stress or plane strain assumption is adopted is not the focus of this example, we simply use a plane stress element with a unit thickness. The material is defined as follows.

material Elastic2D 1 2e4 .0 1e-4 1

The elastic modulus is $20000$ while the density is $0.0001$.

Visualisation

For visualisation, we define a Visualisation recorder. We record von Mises stress to represent the propagation of stress field.

hdf5recorder 1 Visualisation MISES width 5

IBC

The boundaries can be extracted by generating node groups.

generatebyrule nodegroup 1 1 1. 0. # left
generatebyrule nodegroup 2 2 1. 0. # bottom

Then BCs can be applied via groupmultiplierbc.

groupmultiplierbc 1 1 1
groupmultiplierbc 2 1 2

The point load is defined using a Ricker wavelet.

var r:=(x-.1)*20*pi;
(1-2*r^2)*exp(-r*r)

We define an expression to use the above Ricker wavelet as the load.

expression SimpleScalar 1 x ricker
amplitude custom 1 1
cload 1 1 -100 2 3
cload 2 1 -100 1 3

Time Integration

We use both implicit and explicit time integration methods. The step size is unity.

Implicit

The implicit time integration methods are the default. If no integrator is defined, a default Newmark integrator will be used.

An implicit integrator shall be used with an implicit step.

step ImplicitDynamic 1 1
# or just
# step Dynamic 1 1

Explicit

Most explicit methods use acceleration as the primary variable, the equations of motion are often expressed as a function of acceleration. This differs from implicit methods that often use displacement as the primary variable. In order to adopt such a difference, one needs to define an ExplicitDynamic step, similar to the setting in ABAQUS.

step ExplicitDynamic 1 1

Please be aware that most displacement-based constraints cannot be used in explicit analysis.

Other Settings

Since we are using a linear elastic material with the CP4 elements, the elemental stiffness is symmetric. As there are no other non-linear constraints defined in the model, the global stiffness/mass matrix is also symmetric. It is possible to then turn on symmetric banded storage to save memory space.

Also, since the system is linear, the global stiffness/mass matrix does not change once assembled. It is possible to indicate the solver to skip iterations.

set symm_mat 1
set band_mat 1
set linear_system 1

Results

For a not-so-rigorous comparison, different spectral radii are used for different methods, mainly for the purpose of showcasing different methods.

Also, the chosen model parameters are quite arbitrary. Sufficiently accurate results often require an accurate estimation of the highest frequency of the model, which governs the time step size.

Implicit

The Bathe two-step method appears to have the best numerical dispersion performance among implicit methods.

The GSSSS optimal scheme and OALTS perform fine if the spectral radius is chosen properly.

The second-order, unconditionally stable Newmark method has significant high-frequency noise. This explains why it is mainly used for structural dynamics in which the low-frequency response is of interest.

The chosen time step size is $0.005$ seconds for all cases. The following radii are used:

• Implicit Bathe: 0.6

• Implicit GSSSS Optimal: 0.6

• Implicit OALTS: 0.6

• Implicit Newmark: 1.0

implicit methods

Explicit

The explicit methods show better numerical dispersion.

The Tchamwa method is first-order accurate, and does not require corrector.

The Noh-Bathe two-step explicit method, as discussed in the original reference, shows superior performance. However, it requires a corrector step, which requires an additional element-wise computation for each substep. This increases the computational cost.

The WAT2 is a newly developed framework.

The chosen time step size is $0.0005$ seconds. The following radii are used:

• Explicit Bathe: 0.6

• Explicit Tchamwa: 0.6

• Explicit Tchamwa: 0.3

• Explicit WAT2: ~0.36

explicit methods