Abstract

The boundary element method (BEM) is a very effective numerical tool which has been widely applied in engineering problems. Wave equation is a very important equation in applied mathematics with many applications such as wave propagation analysis, acoustics, dynamics, health monitoring and etc. This paper presents to solve the nonlinear 2-D wave equation defined over a rectangular spatial domain with appropriate initial and boundary conditions. Numerical solutions of the governing equations are obtained by using the dual reciprocity boundary element method (DRBEM). Two-dimension wave equation is a time-domain problem, with three independent variables . At the first the Laplace transform is used to reduce by one the number of independent variables (in the present work ), then Salzer's method which is an effective numerical Laplace transform inversion algorithm is used to recover the solution of the original equation at time domain. The present method has been successfully applied to 2-D wave equation with satisfactory accuracy.

1. Introduction

Partial differential equations (PDEs) are governing equations of many physical phenomena such as fluid dynamics, acoustics, electrostatics, electrodynamics, fluid flow, elasticity, health monitoring, or elasticity. The boundary element method is a well-known powerful technique to solve PDEs. The main idea of BEM is to convert the original PDEs to an equivalent boundary integral equation by using Green’s theorem and a fundamental solution. This method has many advantages such as high accuracy due to the semi-analytical nature and use of integrals, increase the efficiency with dimension reduction of problems that resulting a considerable reduction in computer time and cost, high-precision modeling in infinite domain problems and thin shell-like structures/materials, solving large-scale problems with complicated geometries and etc ^{1, 2, 3, 4}. During the last decade many engineering problems have been solved using BEM, such as three dimensional static fracture and fatigue crack growth ⁵, numerical simulation of three-dimensional double-diffusive natural convection in porous media ⁶, solving acoustic problems in open-boundary domains ⁷, solving heat conduction problems ⁸, stability analysis of shallow tunnels subjected to eccentric loads ⁹, modeling of fluid flow in periodic cell with porous cylinder ¹⁰, hydrofoil analysis ¹¹, steady state convection diffusion radiation problems ¹².

Wave equation is a second-order linear hyperbolic partial differential equation that describes the propagation of waves, such as sound or water waves. It is a very important equation with many applicable to the engineering problems. For this reason it has attracted a significant amount of attention in applied researches, such as earthquake engineering, ocean wave propagation, Tsunami prediction, traffic induced vibrations, soil dynamics, acoustics, electromagnetic, fluid dynamics, nondestructive testing, ultrasonic imaging, machinery induced vibrations, safety of nuclear facilities and etc.

BEM is especially suitable for modeling of wave propagation in bounded and unbounded domains. The dual reciprocity method (DRM), which is introduced for the first time by Brebbia and Nardini ¹³ is one of the most efficient methods to eliminate the domain integrals. The main idea of DRM is to expand the inhomogeneous term in terms of its values at the nodes which lie in domain and boundary. The inhomogeneous term is approximated by interpolation in terms of some well-known functions , called radial basis functions (RBFs). By coupling BEM and DRM methods the PDE problem will be reduced to a boundary formulation and there is no any domain integration in the boundary integral equation. This is the main idea behind of DRBEM. Recently, Ghassemi et al. developed to solve Helmholtz equation and its particular solutions using DRBEM employed to the acoustic problem ^{14, 15, 16}.

The main aim of this paper is to solve 2-D wave equation using DRBEM. The 2-D wave is a time-domain PDE problem. In general there are three types of time-domain numerical solution strategies to solve wave equation: 1) Space-time integral equations, 2) Laplace transform method, and 3) Time-stepping methods. In the present study Laplace transform is used. The main idea of Laplace transform is that after applying the Laplace transform to the PDE, the number of independent variables is reduced by one but the equation involves instead a complex parameter . The resulting (elliptic) equation will be solved by any suitable numerical method giving over a spatial grid, for several specified values of the parameter for a given time point. The solution in the time domain will be recovered by a numerical Laplace inversion algorithm.

This paper is organized as follows: in the sections 2 and 3, the Laplace transform and its numerical inverse is introduced, respectively. The governing equation and boundary conditions are described in section 4. The details of Dual Reciprocity Boundary Element Method and how to use it to solve PDEs are summary presented in section 5. The section 6 includes numerical implementation of DRBEM and its results. Finally in section 7, the conclusions drawn from the present work are given.

2. The Laplace Transform

The Laplace transform can be useful in solving ordinary and partial deferential equations. In this present work it is used to eliminate the time-dependency and to produce a subsidiary equation which is then solved in complex arithmetic by DRBEM. Let be a function of . The Laplace transform of is defined as follows:

(1)

When the Laplace transform exists, it is denoted by .

Consider the th order linear PDE in just two independent spatial variables , and time , with coefficients independentof time, defined for over the finite region , with sufficient initial and boundary conditions to give a unique solution to the problem. This function represented as follows ¹⁷:

(2)

where belongs to a certain class of linear operators, functions and are assumed to be continuous.

The Laplace transform of the th order equation with respect to gives:

(3)

where is the Laplace transform of with respect to , is the Laplace transform of and containsall the initial condition information. Eq. (3) is a th order equationin only and which has to be solved for over the domain subject to the transformed boundary conditions for particular values of the complex variable (thus is itself complex). The boundary conditions for Eq. (3) are obtained by taking the Laplace transform with respect to of the boundary conditions specified for Eq. (2).

3. Numerical Inversion of the Laplace Transforms

Analytic inversion of the Laplace transform is defined as a contour integration in the complex plane. The Bromwich contour is commonly chosen. The inversion formula for the Laplace transform is as follows:

(4)

where is a suitable complex number, is a real nonnegative parameter and must be analytic in the right-half plane of the complex plane.

There are several approaches available for the numerical inversion of Laplace transforms ^{18, 19, 20, 21}. In this work the Salzer's method is used for the numerical in version of the Laplace transform. This method is based on the numerical evaluation of the Bromwich integral. Salzer's method is in fact Gaussian quadrature of the Bromwich integral and gives the following approximation to ¹⁶:

(5)

where and are the abscissae and weights of the Gaussian quadrature formula of order . The valuesof and are tabulated for various values of by ^{22, 23, 24}. In the following how to use of Eq. (5) is explained to recover the solution in the time domain.

It is possible to compute from the Eq. (3) using any suitable numerical approaches with its transformed boundary conditions for a fixed value of the complex parameter . The approximate solution is now obtained from the values of . The value of can be identified with , whichis the solution of Eq. (3) with . For each value of , Eq. (3) must be solved for times. Therefore using this process it can be recover the solution of Eq. (2) in the time domain.

4. Governing Equation and Boundary Conditions

Wave equation is a time-space dependent problem. General 2-D wave equation is as follows:

(6)

for which is the unknown solution, and denotes the Laplacian and is a constant speed of the wave propagation, which we assume to be constant in the whole environment. To get the unique solve of Eq. (6) the behavior of the solution on the boundaries of the environment, initial conditions and other specifications of solution needs to be specified. We adopt the unit square as the spatial solution domain with 100 elements per each side and 100 interior points, with initial condition , and on the boundaries. The analytical solution of Eq. (6) is as flows:

(7)

Taking the Laplace transform with respect to converts Eq. (6) into the subsidiary equation, as flows:

(8)

The Eq. (8) has to be solved for .

The subsidiary equation will be solved using DRBEM. The next section describes how to solve this equation.

5. The Use of Dual Reciprocity Boundary Element Method

This section describes how to solve subsidiary Eq. (8) using DRBEM. The general form of Eq. (8) is as follows:

(9)

where includes the coefficients of all the dependent terms and represents the remaining terms.

In this paper to solve the Eq. (9), modified Bessel function is used as the weighting functions and only is considered as source term. The source term is then transferred to the boundary by using dual reciprocity principles. So the interpolation approximation is needed only for source term ²⁵. A brief description of the DRM is given next.

Consider a set of collocation points arbitrarily located in the domain the function may be approximated by the following expression:

(10)

where is sum of boundary nodes and interior (DRM) nodes as seen in Figure 1, subscript denotes any interpolation node include boundary and DRM nodes, and are interpolation coefficients and radial basis function (RBF), respectively. Many functions such as polynomials can be selected as functions. In this paper the following function is selected as (for 2D):

(11)

where . The inverse of matrix , denoted by , can be used to determined vector . The values at the nodal points may be approximated by the following expression:

(12)

So with having values, the values of at any points (other than nodal points) can be calculated using Eq. (10). Therefore the interpolation procedure is completed. It is necessary to note that the interpolation improves with an increase in the number of nodal points.

In this section, the discretization of Eq. (9) is described. The boundary is discretized into constant elements. For more details refer to ^{24, 25}. Based on Eq. (10), the right hand side (RHS) of Eq. (9) can be written as follows:

(13)

for 2D, is as follows:

(14)

The parameter is easily determind by .

Now the Eq. (9) can be written as follows:

(15)

For the source point located at point , let be the weighting function. The weighted residual formulation of Eq. (15) is then given by:

(16)

where is the fundamental solution of , as follows:

(17)

where is the zeroth order Bessel functions of the second kinds.

Now it is possible to take terms on the LHS and RHS to the boundary. Using the Gauss-Green theorem, the RHS of Eq. (16) can be expressed in terms of boundary integral as follows:

(18)

The term is the normal derivateve of and can be expressed as:

(19)

The resulting discretized equation is:

(20)

where

(21)

(22)

(23)

in the above formulas, and denote source point and integration point, respectively. The integrals of Eq. (21) and Eq. (22) can be evaluated by numerical and analytical methods. In this paper numerical gauss quadrature formula is used.

Using the Gauss-Green theorem, the LHS of Eq. (16) can be expressed in terms of boundary integral as follows:

(24)

where , and are obtained by Eq. (21), Eq. (22), and Eq. (23), respectively. The and are boundary condition at source point which denote the values of the boundary quantity and its normalderivative. The detailes are prsented in ²⁶.

By combining the LHS (Eq. (24)) with RHS (Eq.(20)), the problem will be reduced to a boundary only formulation. So it is a discretize version of Eq. (9) using the DRBEM method.

6. Numerical implementation and results

In the present work a simple square domain with 100 elements per each side and 100 interior (DRM) points, is considered. The body is discretized into constant equal elements. The geometry of the all nodes on the domain is shown in the Figure 2. A computer program in Matlab software is developed to all implementations. The results for and are given in Table 1, when the node point is located on the middle of each element. For each point in Table 1, two comparisons are made; firstly with the results from the analytical solution and secondly with those calculated from the DRBEM. By considering Table 1, it can be seen that the analytic solution and DRBEM results are in close agreement almost everywhere. The results are compared in terms of relative error, defined by:

(25)

where and are the approximate and exact values of problem under consideration.

Figure 3 and Table 1 show a comparison of mean relative error for various t (0 to 0.9) and P (1 to 6) values and total mean relative error for all times between 0 to 0.9, respectively. Based on the results of Table 2, in general it can be said that by increasing value, the general accuracy will be increased. In this study the best and the worst results are obtained with and , respectively. Figure 4 shows the numerical solution of the two-dimensional wave equation for various values and . In general the obtained results show that DRBEM is able to solve wave equation with satisfactory results.

7. Conclusions

In this paper we considered the DRBEM solution of 2D wave equation which is a fundamental partial differential equation in many engineering problems. The method involves taking the Laplace transform of the equation, using DRBEM to solve the transformed equation and a numerical Laplace transform inversion algorithm to obtain the solution at specified points. The Salzer's method was used to obtain the numerical Laplace inverse. The use of this method is easy and can give good result for not too large. The experiments show that the results obtained using the DRBEM method have a good relationship with those obtained by the analytical solution for suitable large values. It is possible to use other numerical Laplace inverse algorithms such as Gaver-Stehfest method, Schapery’s method, Talbot method, and Fourier series method, to solve the subsidiary equation. These methods may be providing better results than Salzer's method. It should be noted that the robust methods require complex arithmetic evaluations, but usually provide more reliable results.

Acknowledgements

All computations presented here have been performed on the parallel machines of the high performance computing research center (HPCRC) of Amirkabir University of Technology; their supports are gratefully acknowledged.

References