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.

1 Like