Space-time FEM for elastic wave propagation — no time stepping

Specific Domains Numerics
1

I have been experimenting with a simple space-time finite element formulation for 1D elastodynamics, and I thought the example might be interesting for the Numerics category.

The problem is an elastic bar immediately after impact with a rigid wall. The impact itself is not modeled; the calculation starts from the post-impact initial state and follows the subsequent stress-wave propagation, reflection and release.

Instead of discretizing space first and then integrating the semi-discrete system in time, I introduce the scaled coordinate

y = ct,

where

c = \sqrt{\frac{E}{\rho}}

is the longitudinal wave speed.

The complete problem can then be treated as an ordinary 2D finite-element problem in the (x,y) plane.

Using particle velocity v and normalized stress

s = \frac{\sigma}{\rho c},

the first-order equations become

v_{,y} - s_{,x} = 0,
s_{,y} - v_{,x} = 0.

I used a least-squares formulation based on the residuals

r_1 = v_{,y} - s_{,x}, \qquad r_2 = s_{,y} - v_{,x}.

Minimizing

J(v,s) = \frac12 \int_Q \left( r_1^2 + r_2^2 \right)\,dQ

leads to a symmetric coupled formulation that can be discretized directly with standard continuous Lagrange finite elements.

In LowLevelFEM.jl, the four bilinear forms are written as

Kvv = ∫(Grad(V) ⋅ I ⋅ Grad(V))
Kvs = ∫(Grad(V) ⋅ C ⋅ Grad(S))

Ksv = ∫(Grad(S) ⋅ C ⋅ Grad(V))
Kss = ∫(Grad(S) ⋅ I ⋅ Grad(S))

with

I = [1.0 0.0;
     0.0 1.0]

C = [ 0.0 -1.0;
     -1.0  0.0]

and the complete multifield matrix is simply

K = SystemMatrix([
    Kvv  Kvs
    Ksv  Kss
])

The entire evolution over the chosen space-time domain is then obtained with one solve:

v, s = solveField(K, f, support=bc)

There is no time-stepping loop and no separate transient integration algorithm.

The numerical solution shows the two expected characteristic fronts in the (x,ct) plane:

  • the compression wave traveling from the constrained end to the free end,
  • and the reflected release wave traveling back toward the wall.

The solution also reproduces the classical 1D result. During the compressed phase,

\sigma_c = -\rho c v_0,

and the release wave returns to the wall at

t = \frac{2L}{c}.

What I find particularly interesting in this example is how directly the mathematical formulation maps to executable Julia code.

The complete notebook, including the derivation and plots, is here:

https://github.com/perebalazs/LowLevelFEM.jl/blob/main/examples/space-time-FEM.ipynb

I would be very interested in comments from people who have worked with space-time FEM, least-squares FEM, or related formulations for hyperbolic problems.

2

Hey, I’m not familiar with all the details, but since you mentioned impact and elastic wave propagation in 1D, have a look at

Carlos Yoong, Anders Thorin, Mathias Legrand. The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints. The ASME 2015 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Aug 2015, Boston, United States. ⟨10.1115/DETC2015-46919⟩. ⟨hal-01194922⟩

3

Thanks, this is interesting — the physical benchmark is indeed very similar.
The numerical approach is quite different, though: their WFEM is a time-marching method with \Delta t=\Delta x/c and special treatment of unilateral contact, while my example treats (x,ct) as a single 2D finite-element domain and solves the complete space-time problem in one global least-squares solve.
My example was mainly intended as a compact demonstration of the weak-form/multifield capabilities of LowLevelFEM rather than as a study of space-time methods themselves.

4

Yes, It’s more for the benchmark, these problems are notoriously sensitive to plenty of small effects, numerical dissipation… that can quickly become problematic or not depending on the application. So it was more if you were looking for some comparisons, there you go.
And the paper is just well explained and beautifully written is believe.

5

Time-space FE is great in 1D. The complexity explodes in higher dimensions.

6

True — but for a compact 1D demonstration, I found the formulation too elegant not to try. :slightly_smiling_face: