Method Chaiwoot Boonyasiriwat September 3, 2020 Introduction to - - PowerPoint PPT Presentation

Finite-Difference Time-Domain (FDTD) Method

Chaiwoot Boonyasiriwat

September 3, 2020

2

▪ FDTD is a numerical method for solving time-domain wave equations using the finite difference method. ▪ Explicit FDTD schemes are easy to implement but only suitable for problem domains that can be discretized as structured grids, e.g., rectangle, circle, cube, cylinder, and sphere. ▪ In addition, explicit FDTD schemes are conditionally stable, i.e., a stability condition must be satisfied. ▪ Typically, the second-order centered FD formula is used to approximate the time derivative while the spatial derivatives can be approximated using arbitrary-order FD formulas.

Introduction to FDTD

Approximating the derivatives in the wave equation by the second-order FD approximations yields the finite-difference equation where

Explicit FDTD in 1D

Rearranging the finite difference equation yields the explicit FDTD scheme: where is called the Courant number.

Explicit FDTD in 1D

Consider the wave equation in 3D Applying the Fourier transform in space and time to the wave equation yields the dispersion relation Phase speed is frequency independent. Hence, the wave equation governs nondispersive waves.

Dispersion Relation

Cohen (2002, p. 25-26)

Let’s find a dispersion relation for the semi-discrete wave equation Inserting the plane wave solution of the form into the wave equation and rearranging yields the dispersion relation

Dispersion Relation

Cohen (2002, p. 65-66)

Similarly, for the fully discrete wave equation inserting the plane wave solution of the form and rearranging yields the dispersion relation

Dispersion Relation

Cohen (2002, p. 66)

In a homogeneous medium the wave speed is given by Using an FD approximation, the numerical speed can be computed by Let’s define the numerical dispersion coefficient q as When there is no dispersion error, q = 1.

Numerical Dispersion

Numerical dispersion coefficient can be computed using Phase velocity: Group velocity:

Numerical Dispersion

Cohen (2002, p. 101-102)

▪ “Numerical dispersion produces parasitic waves since the numerical velocity is frequency dependent.” ▪ “These parasitic waves can leave the physical wave and produce some ringing features in the waveform.”

Concept of Numerical Dispersion

Cohen (2002, p. 102-103)

Semi-discrete wave equation using the second-order scheme: using a fourth-order scheme: Note the leading error term in each case.

Order of Numerical Dispersion

Cohen (2002, p. 104)

Fully-discrete wave equation using second-order in time and second-order in space: and fourth-order in space:

Order of Numerical Dispersion

Cohen (2002, p. 105)

Dispersion curve for second-order scheme

Dispersion Curve

K=1/n (n = #points/wavelength) Dispersion Coefficient

Cohen (2002, p. 114)

Order 2 14 Orders 2, 4, 6, 8, 10, 12, 14

▪ To avoid a serious numerical dispersion issue, we need to determine the number of grid points needed for sampling a wavelength. ▪ This leads to a numerical dispersion condition as shown, e.g., in the previous slide. ▪ This condition is used to determine the grid spacing in terms of medium property, e.g. velocity c, and wave frequency f: where n is the number of points per wavelength

Numerical Dispersion Condition

“Accuracy of approximation depends on the direction of wave propagation.”

Isotropy Curve

Cohen (2002, p. 118-119)

Numerical dispersion versus propagation angle

16

▪ Von Neumann stability analysis is the most widely used method for analyzing the stability of a finite difference scheme for time-dependent problems, e.g., wave equations, diffusion equation. ▪ In wave propagation problems, plane wave solutions are used for stability analysis. ▪ However, the method can only be used for linear PDEs. ▪ For nonlinear problems, Lyapunov stability analysis has been widely used.

Linear Stability Analysis

17

Consider an acoustic wave propagating in a 1D homogeneous medium governed by Assuming a plane wave solution given by

• ne can obtain the dispersion relation

Acoustic Plane Wave

18

Consider the case when  is complex The plane wave solution becomes When I > 0, we obtain an evanescent wave. When I < 0, wave amplitude increases exponentially with time. This can happen when a numerical solution is unstable.

Instability

19

Dispersion relation for the second-order scheme is We then obtain This leads to the CFL condition (Courant et al., 1928)

Stability Analysis

20

Consider the first-order coupled acoustic wave equations where p is pressure, u particle velocity, K bulk modulus, and  mass density. Let’s use the second-order FD In both space and time on the staggered grid to solve the acoustic wave equation

1D Acoustic Wave Equations

21

We then obtain the finite difference equations

1D Acoustic Wave Equations

22

Assuming plane wave solutions we can obtain where

Stability Analysis: Acoustics

23

▪ The coefficient matrix is called the amplification matrix whose eigenvalues determine the stability of the FD scheme. ▪ If , wave amplitudes are amplified every time step leading to a blow up and the scheme is unstable. ▪ Otherwise the scheme is stable as long as a stability condition is satisfied. ▪ In this case, the characteristic equation is ▪ The eigenvalues are

Stability Analysis: Acoustics

24

▪ If , “the eigenvalues are real and of magnitude greater than 1, so giving instability (Woolfson and Pert, 1999, p. 163). ▪ If , the eigenvalues are complex. ▪ If , both eigenvalues are equal to -1. ▪ In the last two cases, the scheme is stable. ▪ The stability conditions are satisfied when the CFL condition is satisfied ▪ This is in agreement with the stability analysis of the second-order wave equation using the second-order FDTD scheme.

Stability Analysis: Acoustics

25

▪ Recall the plane wave solution ▪ At the next time step the wave field will become ▪ Phase shift due to wave propagation is ▪ Recall the dispersion relation of the numerical solution ▪ Using the identity , the dispersion relation becomes

Phase Shift

26

▪ Consequently, we obtain the numerical phase error ▪ When C = 1, and there is no phase error, i.e., no numerical dispersion – the FDTD solution is equal to the exact solution.

Phase Shift (continued)

27

In practice, there usually exists an energy source and the governing equation will have an additional source term FD implementation becomes for i = 2,n-1

Source Term

where is is the index of source position.

28

▪ In many situations we want to model wave propagation in unbounded domains. ▪ It is impossible to construct a physical model representing an unbounded domain. ▪ Only part of the domain can be represented. ▪ The region of interest (ROI) is only used in a computational study. ▪ Reflection from ROI boundaries is undesirable and must be reduced as much as possible to avoid the interference of unwanted reflected waves with other waves.

Wave in Unbounded Domains

29

Consider the 1D wave equation where L is the two-way propagation operator which can be decomposed into a concatenation

• f two one-way (paraxial) operators:

Clayton-Engquist ABC

30

Applying each one-way operator to a wave field yields one-way wave equations which govern waves that propagate only in one direction, i.e., left or right.

One-Way Wave Equations

Left-going wave Right-going wave

31

The paraxial operators can be used as an absorbing boundary condition to avoid boundary reflection as follows. This method can perfectly absorbs only normally incident waves.

Clayton-Engquist ABC

Allow left-going wave only Allow right-going wave only

32

▪ Clayton-Engquist (1977) proposed to use paraxial wave equations as absorbing boundary conditions for 2D acoustic and elastic wave propagation. ▪ Consider the dispersion relation in 2D ▪ Rearranging the last equation yields

Paraxial Approximation

33

Padé or rational approximation is usually used to expand the square root term as a rational function where the derivatives of the rational function are the same as the original function up to derivative of order m+n.

Padé Approximation

34

This is equivalent to where is the truncation error. Suppose is an analytic function and can be expanded as a Maclaurin series

Padé Approximation

35

We then have This leads to a system of n+m+1 linear equations whose solution is the coefficients of the rational function.

Padé Approximation

36

▪ Padé approximation has been widely used in computational science research because it usually provides a more accurate function approximation than Taylor’s expansion of the same order. ▪ When Taylor series is divergent, Padé series is usually convergent.

Padé Approximation

37

Using rational approximation, we can obtain various approximations

Paraxial Approximation

0 degree 15 degree 45 degree

Reference: Clayton and Engquist (1977, 1980)

38

Paraxial wave equations corresponding to the dispersion relations are

Paraxial Wave Equations

0 degree 15 degree 45 degree

Reference: Clayton and Engquist (1977)

39

Paraxial Approximation

Image Source: Clayton and Engquist (1977)

40

▪ Applying these paraxial wave equations at top and bottom boundaries will reduce spurious reflections from the boundaries. ▪ The boundary conditions for left and right boundaries can be obtained by switching the locations of x and z in the previous equations. ▪ Higher-order ABCs are too complicated and

• nly the lower-order (2 and 4) ABCs are used.

▪ Many ABCs have been proposed later including that of Keys (1985)

Clayton-Engquist ABC

41

Derive the paraxial dispersion relations

Exercise

where are locations of the left and right boundaries, respectively.

42

▪ Sponge ABL is an artificial layer attached to the physical domain at the boundaries. ▪ Wave fields in the sponge zones are gradually damped to reduce boundary reflection.

Cerjan Sponge ABL

ABL 1D Physical Domain ABL

43

▪ ABL proposed by Cerjan et al. (1985) can be used for modeling in both time-space and time-wavenumber (Fourier) domains. ▪ They argued that previously proposed ABCs cannot be used for modeling in the Fourier domain. ▪ Later Clayton and Engquist (1980) extended their method to the Fourier domain. ▪ Sponge ABL is simpler but more expensive than ABC.

Cerjan Sponge ABL

44

▪ Perfectly matched layer (PML) proposed by Bérenger in 1994 for modeling EM wave propagation is the most widely used method for simulation of wave propagation in unbounded domains. ▪ PML is also an artificial layer to absorb wave energy.

Bérenger PML

PML 1D Physical Domain PML

45

▪ Bérenger’s PML for EM is an artificial medium whose impedance is equal to that of free space, i.e., ▪ In 2D TE case, EM waves in the PML satisfy

Bérenger PML

 = electrical conductivity * = magnetic conductivity

46

▪ To match the impedance of free space, Bérenger’s PML must satisfy the condition ▪ An additional complication in Bérenger work is that the magnetic field must be split into two components ▪ The resulting method is then called split-field PML (SPML) which is conditionally stable.

Bérenger PML

Reference: Abarbanel and Gottlieb (1998)

47

▪ Many efforts have been done to improve the stability of PML leading to nonsplit PML formulations (See Zhou (2003)). ▪ Modern PML formulations can be directly derived by using coordinate stretching (Chew and Liu, 1996) or analytic continuation (Johnson, 2010), namely where vanishes in the physical region.

Nonsplit PML

48

▪ Inserting the stretched coordinate into a plane wave yields an evanescent wave ▪ Using the stretched coordinate, the partial derivative becomes

Coordinate Stretching

49

▪ Let’s apply PML to the acoustic wave equations ▪ The first step is to temporal Fourier transform to the equations and we obtain ▪ Then apply the stretched coordinate to the spatial derivatives.

PML for Acoustic Wave

50

▪ We then obtain ▪ Taking inverse Fourier transform yields the PML formulation of acoustic wave equations

PML for Acoustic Wave

51

▪ In EM problems, the damping function has the physical meaning of conductivity. ▪ In acoustic and elastic cases, there is no direct analogy to any physical quantity but it could be related to viscosity of the medium. ▪ Bérenger (1994, p. 191) used where is the PML thickness.

PML Damping Function

52

Collino and Tsogka (2001, p. 301) used where and R is reflection coefficient.

PML Damping Function

53

▪ Many researchers have combined various methods to develop hybrid methods. ▪ Recently, Liu and Sen (2010) proposed a hybrid method to combine Clayton-Engquist ABC with an absorbing layer. ▪ In the absorbing layer, wave fields are computed by the two-way wave equation and one-way wave equation ▪ Total wave field is a linear combination

Hybrid Methods

▪ Abarbanel, S., and D. Gottlieb, 1988, On the construction and analysis of absorbing layers in CEM, Applied Numerical Mathematics, 27, no. 4, 331-340. ▪ Bérenger, J. P., 1994, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 114, 185-200. ▪ Cerjan, C., D. Kosloff, R. Kosloff, M. Reshef, 1985, A nonreflecting boundary condition for discrete acoustic and elastic wave equations, Geophysics, 50, no. 4, 705-708. ▪ Clayton, R., and B. Engquist, 1977, Absorbing boundary conditions for acoustic and elastic wave equation, Bulletin of the Seismological Society of America, 67, no. 6, 1529-1540. ▪ Clayton, R., and B. Engquist, 1980, Absorbing boundary conditions for wave-equation migration, Geophysics, 45, no. 5, 895-904. ▪ Cohen, G. C., 2002, Higher-Order Numerical Methods for Transient Wave Equations, Springer. ▪ Collino, F., and C. Tsogka, 2001, Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics, 66,

• no. 1, 294-307.

▪ Courant, R., K. Friedrichs, and H. Lewy, 1928, On the partial differential equations of mathematical physics, Physik. Math. Ann., 100, 32-74. ▪ Johnson, S. G., 2010, Notes on perfectly matched layers, Online MIT Course Notes. ▪ Woolfson, M. M., and G. J. Pert, 1999, An Introduction to Computer Simulation, Oxford University Press.

References