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

Finite-Difference Time-Domain (FDTD) Method Chaiwoot Boonyasiriwat September 3, 2020

Introduction to FDTD ▪ 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. 2

Explicit FDTD in 1D 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.

Dispersion Relation 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. Cohen (2002, p. 25-26)

Dispersion Relation 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 Cohen (2002, p. 65-66)

Dispersion Relation Similarly, for the fully discrete wave equation inserting the plane wave solution of the form and rearranging yields the dispersion relation Cohen (2002, p. 66)

Numerical Dispersion 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: Cohen (2002, p. 101-102)

Concept of Numerical Dispersion ▪ “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.” Cohen (2002, p. 102-103)

Order of Numerical Dispersion Semi-discrete wave equation using the second-order scheme: using a fourth-order scheme: Note the leading error term in each case. Cohen (2002, p. 104)

Order of Numerical Dispersion Fully-discrete wave equation using second-order in time and second-order in space: and fourth-order in space: Cohen (2002, p. 105)

Dispersion Curve Dispersion curve for second-order scheme Orders 2, 4, 6, 8, 10, 12, 14 Dispersion Coefficient 14 Order 2 K=1/n (n = #points/wavelength) Cohen (2002, p. 114)

Numerical Dispersion Condition ▪ 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

Isotropy Curve “Accuracy of approximation depends on the direction of wave propagation.” Numerical dispersion versus propagation angle Cohen (2002, p. 118-119)

Linear Stability Analysis ▪ 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. 16

Acoustic Plane Wave Consider an acoustic wave propagating in a 1D homogeneous medium governed by Assuming a plane wave solution given by one can obtain the dispersion relation 17

Instability 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. 18

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

1D Acoustic Wave Equations 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 20

1D Acoustic Wave Equations We then obtain the finite difference equations 21

Stability Analysis: Acoustics Assuming plane wave solutions we can obtain where 22

Stability Analysis: Acoustics ▪ 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 23

Stability Analysis: Acoustics ▪ 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. 24

Phase Shift ▪ 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 25

Phase Shift (continued) ▪ 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. 26

Source Term 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 where is is the index of source position. 27

Wave in Unbounded Domains ▪ 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. 28

Clayton-Engquist ABC Consider the 1D wave equation where L is the two-way propagation operator which can be decomposed into a concatenation of two one-way (paraxial) operators: 29

One-Way Wave Equations Applying each one-way operator to a wave field yields one-way wave equations Left-going wave Right-going wave which govern waves that propagate only in one direction, i.e., left or right. 30

Clayton-Engquist ABC The paraxial operators can be used as an absorbing boundary condition to avoid boundary reflection as follows. Allow left-going Allow right-going wave only wave only This method can perfectly absorbs only normally incident waves. 31

Paraxial Approximation ▪ 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 32

Padé Approximation 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 . 33

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

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

Padé Approximation ▪ 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. 36

Paraxial Approximation Using rational approximation, we can obtain various approximations 0 degree 15 degree 45 degree 37 Reference: Clayton and Engquist (1977, 1980)

Paraxial Wave Equations Paraxial wave equations corresponding to the dispersion relations are 0 degree 15 degree 45 degree 38 Reference: Clayton and Engquist (1977)