Time Dependent Wave Propagation Modeling Using Finite Difference Scheme of 2D Wave Equation Based on Absorbing and Reflecting Boundaries

M. Joni Alam^1*, W. Z. Loskor², M. Mohiuddin¹

¹Department of Mathematics, Comilla University, Cumilla, Bangladesh.
²Department of Science and Humanities, Bangladesh Army International University of Science and Technology, Cumilla, Bangladesh.

DOI: 10.4236/jamp.2021.99148                

Abstract

Boundary procedure is an important phenomenon in numerical simulation. To reduce or eliminate the spurious reflections significantly which is occurred in boundary is a challenging and vital approach. The appropriate artificial numerical boundaries can be applied to eliminate the effect of unnecessary spurious reflections in case of the numerical simulations of wave propagation phenomena problems. Typically, to reduce the artificial reflections, the absorbing boundary conditions are necessary. In this paper, we overview and investigate the appropriate typical absorbing boundary conditions and analyzed the boundary effect of two dimensional wave equation numerically. Reflections over the wide-ranging incident angles are complicated to eliminate, but the absorbing boundary conditions that we have applied are computationally cost efficient, easy to apply and able to reduce reflections significantly. For numerical solution, finite difference method is applied to develop numerical scheme using 2D wave equation. Using the developed numerical scheme, we obtain the numerical solution of the governing equation as an initial boundary value problem and realize the qualitative behavior of the solution in infinite space. The finite difference numerical scheme has been investigated by developing MATLAB programming language code. Numerical results have been discussed and analyzed with presenting different qualitative behavior of the numerical scheme. The accuracy and efficiency of the numerical scheme has been illustrated. The stability analysis was discussed and verified stability condition. Using the numerical scheme and absorbing boundary conditions, the boundary effects and absorption of spurious reflection of boundary have been demonstrated.

Keywords

Absorbing Boundary Conditions, 2D Wave Equation, Finite Difference Method, Numerical Solution

1. Introduction

The solution of wave propagation phenomena in spatially discrete method for the numerical simulation always requires the elimination of the spurious events. The unnecessary spurious events have been produced by the boundaries of the numerical grids. In numerical simulation, a finite region of space is covered by the numerical mesh. As a result, these excessive boundary effects events arise in the system. The boundary effects have been appeared usually as reflection. These unwanted events are always inessential to the real physical events. As a result, the elimination of these events is desired in the study.

In recent years, study says that the explicit finite difference schemes are being used widely in finding the approximate solution for wave propagation problems. The problems of wave propagation phenomena are typically solved for an infinite medium. The solution also can only be solved at finite number of points because of introducing the finite core of computation domain. Besides, the unwanted boundary effect caused by Dirichlet or Neumann conditions cannot be neglected in obtaining finite difference approximation; we should use boundary conditions which reduce edge reflections. Alpert, Greengard, and Hagstrom (2002) described a scalar wave equation with a new and efficient non-reflecting boundary condition and compared the performance with different existing method [1]. The technique described here appears particularly advantageous in case of only long time integrations and the primary limitation of the technique is the propagation medium that must be uniform at the boundary. Appelö, Hagstrom, and Kreiss (2006) presented a model with perfectly matched layer boundary condition with the limitation of some special cases [2]. Benito, Urena, Gavete, Salete, and Muelas (2013) pointed out perfectly matched layer absorbing boundary condition in unbounded domain and shown the avoided reflection over the grid point from edges [3]. Clayton and Enquist (1977) described the elimination of boundary effects [4]. The findings remain some major difficulties in cooperating with corner points of the rectangular domain. Costen and Bérenger (2012) has shown the implementation of absorbing boundary condition in corner regions of computational domain [5]. Gächter and Grote (2003) derived an exact non-reflecting boundary condition which is investigated for three space dimensions [6]. Givoli and Neta (2003) incorporated a finite difference scheme and proposed a non-reflecting boundary scheme for time-dependent wave problems in unbounded domains [7]. The boundary effects along the boundary discussed here caused by normal and shear stress modules. Grote and Keller (2000) demonstrated finite difference method and showed high improvement in accuracy about the implementation of absorbing boundary conditions over standard methods [8]. Hagstrom, Mar-Or, and Givoli (2008) introduced artificial absorbing boundary condition to two dimensional wave guides [9]. Hall and Wang (2009) discussed both reflecting and absorbing boundary conditions in two dimensional models of heterogeneous, anisotropic and fractured media [10]. The artificial boundary is able to reduce the edges reflection but no significance approach is analyzed. Lysmer and Kuhlemeyer (1969) proposed a possible solution concentrating on viscous damping situation [11]. The proposed method can describe the energy propagation in case of outward direction only. However, this method almost reduces reflected compressional waves in boundary. It also works over a long range of incident angles. But, the method is not able to eliminate reflected shear waves completely. Mur (1981) presented absorbing boundary condition for maximum elimination of boundary effect [12]. We overviewed the Mur’s absorbing boundary condition as more efficient and accurate in implementing boundary effect. Reynolds (1978) developed a finite difference models that shows unwanted reflections are produced from the edges of boundary [13]. This happens due to use of boundary conditions modeled by Dirichlet or Neumann. The model was used for generating synthetic seismograms. Reynold in his study shows boundary conditions and the reduction of the edge reflection in a significant manner as well. Smith (1974) proposed a method which was claimed to work perfectly for all incident angles [14]. In this method, the simulation is experimented twice for each absorbing boundary. The first implementation was with Dirichlet boundary conditions and next once with Neumann boundary conditions. Subsequently, the reflections produced by two boundary conditions are in opposite sign. Therefore, sum of such two opposite cases terminate the reflections. Tsynkov (2003) constructed non-local artificial boundary conditions in investigation the numerical simulation of time-dependent acoustic waves for unbounded unobstructed space [15]. Burgos and Santos (2016) presented the formulation of finite elements based on Deslauriers-Dubuc interpolating scaling functions for introducing wave propagation modeling. Results obtained in 2-D are compared with the standard Finite Difference Method for validation [16].

In this study, a systematic method is presented through applying the set of absorbing boundary conditions. These conditions are noticed to be able to eliminate the unnecessary boundary events. The approached method is applied to the two dimensional wave equation in solving initial boundary value problems. The applied absorbing boundary conditions would be able to reduce or absorb the reflecting effect by the boundaries thorough the outward-moving energy.

2. Model and Method

In this study, two dimensional wave equation is considered to solve and analysis the boundary effects in the space-time domain $\Omega \times \left(0,T\right]$, where $\Omega =\left(0,L_x\right)\times \left(0,L_y\right)$ is a rectangular spatial domain. The complete initial-boundary value problem is defined by

$u_{tt}=c^2\left(u_{xx}+u_{yy}\right)+f\left(x,y,t\right),\left(x,y\right)\in \Omega ,t\in \left(0,T\right]$ (1)

$u\left(x,y,0\right)=I\left(x,y\right),\left(x,y\right)\in \Omega$ (2)

$u_t\left(x,y,0\right)=V\left(x,y\right),\left(x,y\right)\in \Omega$ (3)

$u\left(0,y,t\right)=0$ and $u\left(Lx,y,t\right)=0$ (4)

$u\left(x,0,t\right)=0$ and $u\left(x,Ly,t\right)=0$ (5)

where the four sides of the rectangle $\Omega =\left[0,L_x\right]\times \left[0,L_y\right]$ specifies as $x=0,x=L_x$ and $y=0,y=L_y$. The initial condition (2) represents the initial shape of the incidence waves while another initial condition (3) specifies initial velocity of waves. Furthermore, in solving PDEs requires boundary conditions. The Equation (4) and (5) represents fixed boundary in the rectangular domain. Here, the fixed boundaries reflect the incidence waves at the boundary. Later, the investigation is carried out to eliminate the reflection from the boundary by introducing the absorbing boundary conditions.

Waves in one-dimension move along the line. The reflection is occurred here by some particular change or disturbance in the line and consequently waves move backwards. Contrary, when we go through the higher dimensions to know the reflection effects and to know the wave disturbance occurred initially in some localized region and spreads it out is far from obvious. In this study, we initiate absorbing boundary conditions and shows how waves is being absorbed rather than reflected in the boundaries of the rectangular region. The phenomena can be accounted to be investigated like falling a pebble or other solid substance into still water create an outward forwarding circle of ripples. The boundary effects of the produced ripples are minimized in the study by introducing explicit finite difference method.

3. Finite Difference Discretization for 2D Wave Equation

The finite differences approximation is an extensively used method in solving partial differential equations. Considering, d space dimension and constant wave velocity , the wave equation can be written as

$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u\text{for}x\in \Omega \subset {ℝ}^d,t\in \left(0,T\right]$

Consider 2D wave equation with some free region

$\frac{{\partial }^{2}u}{\partial t^2}=c^2\frac{{\partial }^{2}u}{\partial x^2}+c^2\frac{{\partial }^{2}u}{\partial y^2}+f$

which can be discretized as

$\frac{u_{i,j}^{n+1}-2u_{i,j}^n+u_{i,j}^{n-1}}{\Delta t^2}=c^2\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+c^2\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}+f_i^n$

For constant mesh spacing let $\Delta x=\Delta y$. The numerical scheme as follows

$u_{i,j}^{n+1}=2u_{i,j}^n-u_{i,j}^{n-1}+{\left[c\frac{\Delta t}{\Delta x}\right]}^2\left(u_{i+1,j}^n+u_{i-1,j}^n-4u_{i,j}^n+u_{i,j+1}^n+u_{i,j-1}^n\right)+\Delta t^2{f}_i^n$

$u_{i,j}^{n+1}=2u_{i,j}^n-u_{i,j}^{n-1}+{\left[CFL\right]}^2\left(u_{i+1,j}^n+u_{i-1,j}^n-4u_{i,j}^n+u_{i,j+1}^n+u_{i,j-1}^n\right)+\Delta t^2{f}_i^n$

$u_{i,j}^{n+1}=2\left(1-2{\left[CFL\right]}^2\right)u_{i,j}^n-u_{i,j}^{n-1}+{\left[CFL\right]}^2\left(u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n\right)+\Delta t^2{f}_i^n$

4. Numerical Solution

In general, it is very difficult in most case to obtain the exact solution of partial differential equation. In some cases, though analytical solution can be obtained but that is in very complex form. Therefore, one needs to obtain the approximate solution of PDE. Finite difference method is one of the standard methods that have been used to find the approximate solution in solving partial differential equations. In this section, the graphical representation of wave propagation in two dimensions using the developed numerical scheme is presented for different time-step selection and investigated the boundary effects.

The time-dependent numerical simulation for two dimension wave equation can be described by the above numerical results. Here, the four numerical results demonstrated by using the developed finite difference scheme. We consider the constant propagation speed $c=1$ for each case. The temporal domain was considered as [0, 10]. The propagation of waves created from a source towards the boundary is shown here for 1st sec in Figure 1(a), 2nd sec in Figure 1(b), 5th sec in Figure 1(c) and finally at last 10th sec in Figure 1(d) respectively.

5. Stability Analysis of Numerical Scheme

The stability of numerical scheme depends on CFL condition. The CFL stands for the developed numerical scheme as $\frac{c\Delta t}{\Delta x}$. Optimization of the plane sampling intervals $\left(\Delta x,\Delta y\right)$ depends on some criteria while we select time sampling interval that also depends on some criteria as well. The numerical scheme shows

Figure 1. Numerical solution of 2-dimensional wave equation for different time-step.

instable results when these criteria are exceeded. The attained CFL values which are more different than the values are supposed to be obtained then the numerical solution are not structured to be a validated model.

Stability depends on the term $\left(1-2{\left[CFL\right]}^2\right)$ which comes from the developed finite difference scheme in this study. The value of $\left(1-2{\left[CFL\right]}^2\right)$ should be considered as less than or equal to zero otherwise stable solution cannot be obtained so far. As a result, $CFL\le \pm 1/\sqrt{2}$ is found from $\left(1-2{\left[CFL\right]}^2\right)\le 0$. Because the CFL value cannot be taken as negative value in the implementation, $\frac{c\Delta t}{\Delta x}\le 1/\sqrt{2}$ must be carried out to obtain the stable solution.

Here, Figure 2(a) shows the rate of propagation of waves towards the boundaries with different angle of incidence. The CFL condition does not exceed the limit criteria as we choose the CFL as 0.5 which is less than $1/\sqrt{2}$. Contrary, Figure 2(b) shows no visible wave is being propagated here as CFL condition exceeds here. In this case, we select CFL value as 0.8 which is greater than $1/\sqrt{2}$. So, the stability criteria $\frac{c\Delta t}{\Delta x}\le 1/\sqrt{2}$ is failed which results the instable graph. The investigation satisfied the stability criteria properly.

6. Reflection of Waves in Boundaries

The wave can be undergoing by reflection or be absorbed by the boundary. In the wave propagation model described in Equation (1) we consider the end points as Dirichlet boundary conditions given by

$u\left(0,y,t\right)=0,u\left(Lx,y,t\right)=0$ and $u\left(x,0,t\right)=0,u\left(x,Ly,t\right)=0$

The following numerical results are experimented by using the developed

(a) (b)

Figure 2. Stability test of finite difference scheme.

numerical scheme of 2D wave equation with the above mentioned Dirichlet boundary conditions.

Here, the Dirichlet boundary condition is applied and the reflections of waves in boundaries have been noticed. We investigate the effect of spurious reflection for different time step selection. We carried out the investigation for 10 sec. We noticed that, at 7 sec the waves hit the boundary and reflected backward which is shown in Figure 3(a). The reflection is growing significantly with time. Figure 3(b) shows the reflection at 8sec while Figure 3(c) gives the results of reflection at 9 sec. Finally, we experienced the maximum reflection at the boundary at 10sec in Figure 3(d).

7. Absorbing Boundary Condition for 2D Wave Equation

In this section we consider a finite two-dimensional plane. We reform the unbounded boundary restraint by introducing the following artificial boundary. In this computation, we restrict the computation to the square like $\Omega =\left[-L,L\right]\times \left[-L,L\right]$. Under our assumption, there is no other source term is available outside the computational domain Ω. Also we assume that no initial perturbations are present there. The following absorbing boundary conditions ensure the significant reduction of spurious reflection of waves. Even though the boundary effects are reduced here in case of the propagation occurred in infinitely many directions. The sustainable solution of boundary effect for different angle of incidence can be described by the following boundary condition. In this article, we focus on absorbing boundary conditions known as nonreflecting boundary conditions. The boundary effects are computed in both space and time.

Figure 3. Numerical results of 2D wave equation for Dirichlet boundary.

For absorption across the boundary, the absorbing boundary conditions are

${\frac{\partial u}{\partial x}|}_{x=0}=c{\frac{\partial u}{\partial t}|}_{x=0}$, ${\frac{\partial u}{\partial x}|}_{x=Lx}=-c{\frac{\partial u}{\partial t}|}_{x=Lx}$

${\frac{\partial u}{\partial y}|}_{y=0}=c{\frac{\partial u}{\partial t}|}_{y=0}$, ${\frac{\partial u}{\partial y}|}_{y=Ly}=-c{\frac{\partial u}{\partial t}|}_{y=Ly}$

which can be discretized by

$u_{1,j}^{n+1}=u_{2,j}^n+\frac{CFL-1}{CFL+1}\left[u_{2,j}^{n+1}-u_{1,j}^n\right]$

$u_{n_x,j}^{n+1}=u_{n_x-1,j}^n+\frac{CFL-1}{CFL+1}\left[u_{n_x-1,j}^{n+1}-u_{n_x,j}^n\right]$

$u_{i,1}^{n+1}=u_{i,2}^n+\frac{CFL-1}{CFL+1}\left[u_{i,2}^{n+1}-u_{i,1}^n\right]$

$u_{i,n_y}^{n+1}=u_{i,n_y-1}^n+\frac{CFL-1}{CFL+1}\left[u_{i,n_y-1}^{n+1}-u_{i,n_y}^n\right]$

The discretized form of the described absorbing condition is implemented in finite computational domain. The result is illustrated using MATLAB and investigated how the reflection is reduced or eliminated. The boundary conditions presented here are more stable and it is said to be computationally efficient as well. The numerical experiments and boundary effect of absorbing boundary condition is demonstrated in this study graphically.

Figure 4. Results of absorbing boundary condition for different time steps.

Here, we have presented the effect of non-reflecting and absorbing boundary condition for different time steps selection. After creating the waves from a source, the results of last four seconds have been shown here in this article. We see that, the waves are not reflected in each case as the propagation goes on. In Figure 4(a), we see that at 7 sec of the computation, the waves are being absorbed by the boundaries. Then, gradually in Figures 4(b)-(d), we have chosen the time selection for 8 sec, 9 sec and 10 sec respectively and found that all waves are being absorbed by the boundaries. The spurious reflections are being eliminated significantly here is noticed.

8. Conclusion

In numerical simulation for seismic wave fields, absorbing boundary condition plays a significant role in reducing boundary reflection. It is often essential to introduce artificial boundaries to limit the area of computation. Elimination of reflection effects by the boundaries may lead the controlling of different devastating disaster. The absorbing boundary condition and its effect are investigated in this paper. The applied boundary condition not only guarantees stable difference approximation, but also minimizes the spurious reflections that occur at the boundaries.

Acknowledgements

This work is supported by Comilla University Research Fund, Cumilla-3506, Bangladesh.

Conflicts of Interest

The authors confirm that there is no conflict of interest to declare for this publication.

References

[1]	Alpert, B., Greengard, L. and Hagstrom, T. (2002) Nonreflecting Boundary Condition for the Time-Dependent Wave Equation. Journal of Computational Physics, 180, 270-296. https://doi.org/10.1006/jcph.2002.7093
[2]	Appelö, D., Hagstrom, T. and Kreiss, G. (2006) Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-Posedness, and Stability. SIAM Journal of Applied Mathematics, 67, 1-23. https://doi.org/10.1137/050639107
[3]	Benito, J.J., Urena, F., Gavete, L., Salete, E. and Muelas, A. (2013) A GFDM with PML for Seismic Wave Equations in Heterogeneous Media. Journal of Computational and Applied Mathematics, 252, 40-51. https://doi.org/10.1016/j.cam.2012.08.007
[4]	Clayton, R.W. and Enquist, B. (1977) Absorbing Boundary Conditions for Acoustic and Elastic Wave Equations. Bulletin of the Seismological Society of America, 67, 1529-1539.
[5]	Costen, F. and Bérenger, J.-P. (2010) Implementation of the Huygens Absorbing Boundary Condition in Corner Regions. IEEE Transactions on Electromagnetic Compatibility, 54, 367-374. https://doi.org/10.1109/TEMC.2012.2186302
[6]	Gächter, G.K. and Grote, M.J. (2003) Dirichlet-to-Neumann Map for Three Dimensional Elastic Waves. Wave Motion, 37, 293-311. https://doi.org/10.1016/S0165-2125(02)00091-4
[7]	Givoli, D. and Neta, B. (2003) High-Order Non-Reflecting Boundary Scheme for Time-Dependent Waves. Journal of Computational Physics, 186, 24-46. https://doi.org/10.1016/S0021-9991(03)00005-6
[8]	Grote, M.J. and Keller, J.B. (2000) Nonreflecting Boundary Conditions for Elastic Waves. SIAM Journal of Applied Mathematics, 60, 803-819. https://doi.org/10.1137/S0036139998344222
[9]	Hagstrom, T., Mar-Or, A. and Givoli, D. (2008) High-Order Local Absorbing Conditions for the Wave Equation: Extensions and Improvements. Journal of Computational Physics, 227, 3322-3357. https://doi.org/10.1016/j.jcp.2007.11.040
[10]	Hall, F. and Wang, Y. (2009) Elastic Wave Modeling by an Integrated Finite Difference Method. Geophysical Journal International, 177, 104-114. https://doi.org/10.1111/j.1365-246X.2008.04065.x
[11]	Lysmer, J. and Kuhlemeyer, R.L. (1969) Finite Dynamic Model for Infinite Media. Journal of Engineering Mechanics Division, 95, 859-878. https://doi.org/10.1061/JMCEA3.0001144
[12]	Mur, G. (1981) Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic Field Equations. IEEE Transactions on Electromagnetic Compatibility, 23, 377-382. https://doi.org/10.1109/TEMC.1981.303970
[13]	Reynolds, A.C. (1978) Boundary Conditions for the Numerical Solution of Wave Propagation Problems. Geophysics, 43, 1099-1110. https://doi.org/10.1190/1.1440881
[14]	Smith, W.D. (1974) A Non-Reflecting Plane Boundary for Wave Propagation Problems. Journal of Computational Physics, 15, 492-503. https://doi.org/10.1016/0021-9991(74)90075-8
[15]	Tsynkov, S.V. (2003) Artificial Boundary Conditions for the Numerical Simulation of Unsteady Acoustic Waves. Journal of Computational Physics, 189, 626-650. https://doi.org/10.1016/S0021-9991(03)00249-3
[16]	Burgos, R.B. and Santos, M.A.C. (2016) Finite Elements Based on Deslauriers-Dubuc Wavelets for Wave Propagation Problems. Applied Mathematics, 7, 1490-1497. https://doi.org/10.4236/am.2016.714128