Transparent boundary conditions for the elastic Transparent boundary - - PowerPoint PPT Presentation

HYP 2006, Lyon, July 17-21 1

Transparent boundary conditions for the elastic Transparent boundary conditions for the elastic waves in anisotropic media waves in anisotropic media

Ivan L. Sofronov, Nickolai A. Zaitsev Keldysh Institute of Applied Mathematics RAS, Moscow

Outline Outline

• Problem formulation
• Generation of LRBC operator, main steps
• Numerical tests
• Conclusions

HYP 2006, Lyon, July 17-21 2

Motivation Motivation

Generation of low-reflecting boundary conditions on open boundaries for anisotropic elastodynamics is a challenging problem. The PML method is not stable in this case, see

• Beaches E., Fauqueux S., Joly P., Stability of perfectly matched layers group

velocities and anisotropic waves, JCP, 188, 2003, 399 – 433;

• D. Appelö and G. Kreiss, A New Absorbing Layer for Elastic Waves, JCP, 215,

2006, 642 – 660.

HYP 2006, Lyon, July 17-21 3

Low Low-

• reflecting boundary conditions (LRBC)

reflecting boundary conditions (LRBC)

LRBC on “open boundaries” are needed for modeling of wave propagation using a bounded computational domain. Typical setup: Remark: In the domain of interest the governing equations and geometry can be much more complex than outside C. Generation of LRBC on G is an additional problem which is set up by considering auxiliary external Initial Boundary Value Problems (IBVPs) outside C. Domain of interest Computational domain C

Γ

HYP 2006, Lyon, July 17-21 4

Here is density; is the Cartesian velocities; are (constant) elastic coefficients of the Hook’s law written in the matrix form, . In polar coordinates for the velocity vector the system reads:

Anisotropic elasticity, 2D orthotropic media Anisotropic elasticity, 2D orthotropic media

We consider 2D elastodynamic equations:

2 2 2 2 11 22 12 1 2 2 2 2

f f f f f f A A A A A A f t r r r θ θ θ ∂ ∂ ∂ ∂ ∂ ∂ = + + + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ( ) ,

r

f v vθ = ( , ), ( , ) θ θ

i j i

A r A r

( ) ( )

2 2 2 2 11 33 33 12 1 1 1 2 2 2 2 1 2 1 2 2 2 2 2 33 22 33 12 2 2 2 1 2 2 2 1 2 1 2

, . v v v v c c c c t x x x x v v v v c c c c t x x x x ρ ρ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ ∂

where are 2x2 matrices

ρ ( )

1 2

, v v

{ }

nm

c

n nm m

c σ ε =

HYP 2006, Lyon, July 17-21 5

where is Dirac’s delta function; is a basis on G , i.e. (in 2D it is constructed, e.g., by using )

Generation of LRBC operator, main steps Generation of LRBC operator, main steps

Governing equation in 2D space:

{sin ,cos } θ θ k k

Main steps: Stage 1: consider a set of auxiliary external IBVPs outside C (set wrt “m”):

2

0 in / | (1) | ) ( ( )

m m t m t m t m

L R C t

= Γ

 − =  =   =  δ ϕ θ ϕ θ E E E E Γ − =

tt

f Lf { ( )} ϕ θ

∞ = m m

( , ) ( ) ( ) θ ϕ θ = ∑

m m m

f t c t

( ) t δ

HYP 2006, Lyon, July 17-21 6

Generation of LRBC operator (cont.) Generation of LRBC operator (cont.)

Stage 2: make Laplace transform and pass to elliptic BVPs (parameterized by ):

C

2 2

ˆ ˆ 0 in / ˆ | ( ) (2) ˆ |

m m m m m r

s L R C ϕ θ

Γ →∞

 − =   =   =   E E E E Γ

Stage 3: solve (numerically) the problems and evaluate

Thus we obtain the Dirichlet-to-Neumann maps

ˆ ( , ) on θ ∂ Γ ∂

m r

n E s ˆ ( ) ( )[ ] ( , ), ϕ θ ψ θ θ

Γ

∂  ≡ =   ∂   a

m m m

s r r R n E

HYP 2006, Lyon, July 17-21 7

Thus we obtain the Poincare-Steklov operator in space of Fourier coefficients:

• r, in matrix form:

Generation of LRBC operator (cont.) Generation of LRBC operator (cont.)

Stage 4: form matrix of the Poincare-Steklov operator:

we take arbitrary data on G

ˆ ˆ ˆ ( ) ( ) ( ) = ∑

n m m n n

d s P s c s ˆ ˆ ˆ ( ) ( ) ( ), = s s s d P c ˆ ˆ ( , ) ( ) ( ) θ ϕ θ = ∑

m m m

f s c s

{ }

1

ˆ ˆ ˆ , , ...

T

c c = c

and write the representation of its normal derivative on G ˆ ˆ( , ) ( ) ( ) θ ϕ θ   ∂ =   ∂  

∑

n n n

d s f s n

ˆ ˆ ˆ ( , ) ( ) ( )[ ] ( ) ( ) ( ) θ ψ θ ϕ θ ∂ = = ∂

∑ ∑ ∑

m m m m n n m m n

f s c s s c s P s n

HYP 2006, Lyon, July 17-21 8

Generation of LRBC operator (cont.) Generation of LRBC operator (cont.)

First we represent matrix by sum of three matrices to take into account asymptotic at

Stage 5: make inverse Laplace transform for the P-S operator.

ˆ ( ) s P : → ∞ s

1 1

ˆ ˆ ˆ ( ) ( ); , are , ( ) (1) = + + = s s s consts s

• P

P P K P P K

Then we calculate rational approximations to each entry in such that all poles have negative real parts, i.e.

ˆ ( ) s K

1

ˆ ˆ ( ) ( ) , Re( )

m n

L m m m m n l n n n l m l n l

K s K s s

=

≈ ≡ ≤ < −

∑

% α β δ β

Finally the inverse Laplace transform of gives:

1

( ) ( ) ( ) ( ) ( ) ∂ = + + ∗ ∂ % t t t t t t c d P P c K c ˆ ˆ ˆ ( ) ( ) ( ) = s s s d P c

HYP 2006, Lyon, July 17-21 9

Stage 6: compose azimuth modes:

Denote by the operator of Fourier decomposition for i.e. Consequently,

Remark: the explicit form of kernels is

that permits to treat convolutions by stable recurrent formulas.

Generation of LRBC operator (cont.) Generation of LRBC operator (cont.)

1

( ) exp( ), Re( ) α β β δ

=

= ≤ <

∑

%

m n

L m m m m n n l n l n l l

t K t

{ }

1

( ) ∂ ∂ − + + ∗ = ∂ ∂ % f f

• f

t f t n

1 1 1

Q PQ Q P Q Q K Q ( ) % m

n

K t ( , ) ( ) ( ) θ ϕ θ = ∑

m m m

f t c t Q

{ }

: ( , ) ( ) . θ →

m

f t c t Q

The LRBC in the physical space reads:

{ } { }

1 :

( ) , ( ) ,

−

∂ → ∂

m m

f c t d t f n Q

HYP 2006, Lyon, July 17-21 10

Governing equations are implemented in the polar system of coordinates. We consider the task in a circle:

• At

2 R =

we prescribe Dirichlet data (pulse) to initiate elastic waves;

• G with radius

10 RG =

is the external boundary where we put LRBC (with 36 azimuth harmonics).

R0

Γ

RG (LRBC)

Setup of the test problem Setup of the test problem

HYP 2006, Lyon, July 17-21 11

Test problem: reference solution Test problem: reference solution

RG (LRBC)

The verification of LRBC is made by comparing with the reference solution of a second problem having 8x bigger external radius (extended domain). Comparison of two solutions in C-norm is made at r < RG.

RExtended

HYP 2006, Lyon, July 17-21 12

Test calculations Test calculations

Parameters of anisotropic media are taken from anisotropic medium, case IV:

( ) ( )

2 2 2 2 11 33 33 12 1 1 1 2 2 2 2 1 2 1 2 2 2 2 2 33 22 33 12 2 2 2 1 2 2 2 1 2 1 2

, . v v v v c c c c t x x x x v v v v c c c c t x x x x ρ ρ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ ∂

Beaches E., Fauqueux S., Joly P., Stability of perfectly matched layers group velocities and anisotropic waves, JCP, 188, 2003, 399 – 433:

HYP 2006, Lyon, July 17-21 13

Anisotropic case Anisotropic case-

• IV, t=0.52

IV, t=0.52

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 14

Anisotropic case Anisotropic case-

• IV, t=1.92

IV, t=1.92

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 15

Anisotropic case Anisotropic case-

• IV, t=2.62

IV, t=2.62

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 16

Anisotropic case Anisotropic case-

• IV, t=3.32

IV, t=3.32

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 17

Anisotropic case Anisotropic case-

• IV, t=4.01

IV, t=4.01

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 18

Anisotropic case Anisotropic case-

• IV, t=4.71

IV, t=4.71

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 19

Anisotropic case Anisotropic case-

• IV, t=5.41

IV, t=5.41

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 20

Anisotropic case Anisotropic case-

• IV, t=6.11

IV, t=6.11

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 21

Anisotropic case Anisotropic case-

• IV, t=6.81

IV, t=6.81

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 22

Anisotropic case Anisotropic case-

• IV, t=7.50

IV, t=7.50

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 23

Anisotropic case Anisotropic case-

• IV, t=12.0

IV, t=12.0

LRBC Reference Ur Utheta

HYP 2006, Lyon, July 17-21 24

Difference between LRBC and Reference solution Difference between LRBC and Reference solution

t=2.6 t=7.5 t=5.4 t=12.7

HYP 2006, Lyon, July 17-21 25

LRBC stability at large times LRBC stability at large times

HYP 2006, Lyon, July 17-21 26

• A novel approach to generate numerical LRBCs for anisotropic time-domain

wave propagation problems is proposed and implemented in 2D

• LRBCs are efficient in both isotropic and anisotropic media, including the case

where PML approach fails (case IV)

• LRBC operator does not depend on the meshing inside the computational

domain

• LRBC efficiency (amplitude of reflected waves) can be automatically controlled

during generation of the operator

• It is possible to generate in advance a library of LRBC matrices for given media

parameters and geometry of computational domain

• The algorithm is highly parallelized

Conclusions Conclusions

HYP 2006, Lyon, July 17-21 27

• The approach is still computationally expensive and must be revised and

improved to enhance the performance

• Real*16 accuracy
• 3D problems

Further research Further research

HYP 2006, Lyon, July 17-21 28

Stage 2 [Laplace image]: we take a finite interval

and choose set of knots which are representative enough to consider discrete counterparts of kernels for rational approximation (Stage 5). In particular, the Chebyshev’s nodes are used: Main features of the algorithm for discretized system:

Stages 1–6: we take discrete azimuth basis

LRBC operator, numerical aspects LRBC operator, numerical aspects

( )

{sin ,cos }, 0,1,..., ; 20 ; θ θ = = M M k k k

max

[0, ] S

max

{ } [0, ]

j

s S ∈

ˆ ( )

m n j

P s ~ 100 J

max

0.5 1 cos , 1,..., 2

j

S j s j J J π −     = − =        

1

ˆ ( )

m m m n n n

P s P s P − −

HYP 2006, Lyon, July 17-21 29

Stages 3-4 [P-S operator]: we discretise governing equations by pseudospectral derivatives

in azimuth direction (uniform grid) and finite differences in radial direction (exponential grid) and make massive calculations to evaluate each with a guarantied accuracy, (controlled by mesh convergence). Number of separate tasks is Thus we form the matrix

LRBC operator, numerical aspects (cont.) LRBC operator, numerical aspects (cont.)

ˆ ( )

m n j

P s

10

~ 10

PS

ε

−

ˆ ( ).

j

s P 4*( 1)* ; + M J

C

2 2

ˆ ˆ 0 in / ˆ | ( ) ˆ |

m m m m m r

s L R C ϕ θ

Γ →∞

 − =   =   =   E E E E Γ

HYP 2006, Lyon, July 17-21 30

Stage 5 [inverse Laplace transform]:

Constant matrices are estimated from the rational approximations

• n the interval (

is calculated by Chebyshev-Pade algorithm), and the reminder matrix is formed: Rational approximations are calculated for each : (indices at are omitted); L=8 . The minimization is made by an optimization algorithm.

1

, P P

1

ˆ ˆ ( ): ( ) = − −

j j j

s s s K P P P

2 , 1

{ }

ˆ ( ) min , Re( )

l l

L l j l j l l

s j

K s s

α β

α β δ β

=

− → ≤ < −

∑

ˆ , , , K L α β α β

m n

ˆ ( )

m n j

K s

LRBC operator, numerical aspects (cont.) LRBC operator, numerical aspects (cont.)

ˆ ( ) ( )

m m n j n j

R s P s ≈

max

[0, ] S

( )

m n

R s

HYP 2006, Lyon, July 17-21 31

As a result we obtain the rational functions satisfying They are explicitly inverted from the Laplace space:

Stage 6 [discrete NRBC]: introducing matrices of the discrete

Fourier transform for vector-functions in sin-cos basis, we write

• ut our discrete BC:

{ }

1

( )

M M M M M M

f f

• f

t f t n ∂ ∂ − + + ∗ = ∂ ∂

1 1 1

Q PQ Q P Q Q K Q % ,

M M

• 1

Q Q ( ) ,

r

f v vθ = ˆ ( )

m n

K s %

1

ˆ ( ) ( ) exp( ), Re( ) α β β δ

=

≡ ≤ <

∑

% % a

m n

L m m m m m n n n l n l n l l

K s K t t

8

ˆ ˆ | ( ) ( ) | ( ~ 10 )

R m m n j n j R

K s K s ε ε

−

− < %

LRBC operator, numerical aspects (cont.) LRBC operator, numerical aspects (cont.)

Remark: in isotropic case

1 1

, , ( ) is diagonal = = % p p t P I P I K

HYP 2006, Lyon, July 17-21 32

LRBC matrix, residual LRBC matrix, residual

ˆ ˆ max | ( ) ( ) | − %

j

m m n j n j s

K s K s { }

1

( ) ∂ ∂ − + + ∗ = ∂ ∂ %

M M M M M M

f f

• f

f t n t

1 1 1

Q P Q Q P Q Q Q K

HYP 2006, Lyon, July 17-21 33

LRBC matrix, amplitude LRBC matrix, amplitude max |

( ) | % m

n t

K t { }

1

( ) ∂ ∂ − + + ∗ = ∂ ∂ %

M M M M M M

f f

• f

f t n t

1 1 1

Q P Q Q P Q Q Q K

HYP 2006, Lyon, July 17-21 34

LRBC matrix, amplitude LRBC matrix, amplitude

1

| |

m n

P

{ }

1

( ) ∂ ∂ − + + ∗ = ∂ ∂ %

M M M M M M

f f

• f

t f t n

1 1 1

Q Q Q P Q Q K Q P

HYP 2006, Lyon, July 17-21 35

LRBC matrix, entry LRBC matrix, entry

(14,1) (18,3) = → = n m ˆ ˆ ( ) ( ) − %

m m n n

K s K s ˆ ( )

m n

K s ˆ poles of ( ) % m

n

K s ( ) % m

n

K t