Finite Elements: 1D acoustic wave equation

Finite Elements: 1D acoustic wave equation

¾ Helmholtz (wave) equation (time-dependent)

¾ Regular grid

¾ Irregular grid

¾ Explicit time integration

¾ Implicit time integraton

¾ Numerical Examples

Scope: Understand the basic concept of the finite element method applied to the 1D acoustic wave equation.

Acoustic wave equation in 1D

How do we solve a time-dependent problem such as the acoustic wave equation?

where v is the wave speed.

using the same ideas as before we multiply this equation with

an arbitrary function and integrate over the whole domain, e.g. [0,1], and after partial integration

f u

v

t

u − Δ =

∂

dx f

dx u

v dx

u

t

∫ ∫

∫ ^{∂ − ∇ ∇ =}

0 1

0 2 1

0

2

ϕ ϕ ϕ

.. we now introduce an approximation for u using our previous basis functions...

Weak form of wave equation

) ( )

~ (

1

x t

c u

u

N

i

i

ϕ

∑

=

≈

together we obtain

) ( )

~ (

1 2 2

2

u u c t

x

N

i

i t

t

t

∑ ϕ

=

∂

=

∂

≈

∂

note that now our coefficients are time-dependent!

... and ...

∑ ∫ ∫

∑ ∫ ^⎥ _⎦ ^{⎤ + ⎢} _⎣ ^{⎡ ∇ ∇ ⎥} _⎦ ^{⎤ =}

⎢ ⎣

⎡ ∂

0 1

0 2

1

0 2

j j

i

i i

j i

i i

t

c ϕ ϕ dx v c ϕ ϕ dx f ϕ

which we can write as ...

Time extrapolation

∑ ∫ ∫

∑ ∫ ^⎥ _⎦ ^{⎤ + ⎢} _⎣ ^{⎡ ∇ ∇ ⎥} _⎦ ^{⎤ =}

⎢ ⎣

⎡ ∂

0 1

0 2

1

0 2

j j

i

i i

j i

i i

t

c ϕ ϕ dx v c ϕ ϕ dx f ϕ

... in Matrix form ...

g c

A v c

M

&& +

=

M A b

... remember the coefficients c correspond to the

actual values of u at the grid points for the right choice of basis functions ...

How can we solve this time-dependent problem?

stiffness matrix mass matrix

Time extrapolation

... let us use a finite-difference approximation for the time derivative ...

g c

A v c

M

&& +

=

... leading to the solution at time t_k+1 :

g c

A dt v

c c

M

c

⎟ +

=

⎠

⎜ ⎞

⎝

⎛

− +

2

1

1

2

[ ]

2 2

1

1

( )

( ) 2

+

=

−

+

−

k

M g v A c dt c c

c

we already know how to calculate the matrix A but how can we calculate matrix M?

Mass matrix

∑ ∫ ∫

∑ ∫ ^{⎥ =}

⎦

⎢ ⎤

⎣

⎡ ∇ ∇

⎥ +

⎦

⎢ ⎤

⎣

⎡ ∂

0 1

0 2

1

0 2

j j

i

i i

j i

i i

t

c ϕ ϕ dx v c ϕ ϕ dx f ϕ

... let’s recall the definition of our basis functions ...

∫

=

1

0

dx M

ϕ

ϕ

i i

i

i i

i x x x

elsewhere h x h for

x

x h

h for x

x = −

⎪⎪

⎪

⎩

⎪⎪

⎪

⎨

⎧

<

<

−

≤

<

− +

=

−

−

,~ 0

0 ~

~ 1

~ 0 1

~

~) (

1 1

ϕ

i=1 2 3 4 5 6 7 + + + + + + +

h

h

h

h

h

h

... let us calculate some element of M ...

Mass matrix – some elements

3 3

1

1

0

2

0

2

1 1

0

1

i i

h

i h

i i

i ii

h h

h dx dx x

h dx x

M

i i

+

=

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎛ −

⎟⎟ +

⎠

⎜⎜ ⎞

⎝

= ⎛

=

−

−

∫

∫

∫ ^{ϕ ϕ}

⎪⎪

⎪

⎩

⎪⎪

⎪

⎨

⎧

<

<

−

≤

<

− +

=

−

−

elsewhere h x h for

x

x h h for

x

x _i

i

i i

i

0

0 ~

~ 1

~ 0 1

~

~) (

1 1

ϕ

i=1 2 3 4 5 6 7 + + + + + + +

h

h

h

h

h

h

Diagonal elements: M_ii , i=2,n-1

h_i x_i j_i

h_i-1

Matrix assembly

% assemble matrix Mij

M=zeros(nx);

for i=2:nx-1, for j=2:nx-1,

if i==j,

M(i,j)=h(i-1)/3+h(i)/3;

elseif j==i+1 M(i,j)=h(i)/6;

elseif j==i-1 M(i,j)=h(i)/6;

else

M(i,j)=0;

end end end

% assemble matrix Aij

A=zeros(nx);

for i=2:nx-1, for j=2:nx-1,

if i==j,

A(i,j)=1/h(i-1)+1/h(i);

elseif i==j+1

A(i,j)=-1/h(i-1);

elseif i+1==j

A(i,j)=-1/h(i);

else

A(i,j)=0;

end end end

M_ij A_ij

Numerical example

Implicit time integration

... let us use an implicit finite-difference approximation for the time derivative ...

g c

A v c

M

&& +

=

... leading to the solution at time t_k+1 :

g c

A dt v

c c

M

c

⎟ +

=

⎠

⎜ ⎞

⎝

⎛ − +

− + +

1 2

2

1

1

2

[ ] ( ⁽

⁾ )

1 2 2

2

1 ⁻

2

+

=

+

+

−

k

M v dt A gdt M c c

c

How do the numerical solutions compare?

Summary

The time-dependent problem (wave equation) leads to the introduction of the mass matrix.

The numerical solution requires the inversion of a system matrix (it may be sparse).

Both explicit or implicit formulations of the time-dependent part are possible.

The time-dependent problem (wave equation) leads to the introduction of the mass matrix.

The numerical solution requires the inversion of a system matrix (it may be sparse).

Both explicit or implicit formulations of the time-dependent part are possible.