ABSTRACT
The present study analyzes the 2D scalar wave propagation problems in a time domain BEM through large and complex geophysical structures. To achieve this, a time domain BEM formulation is extended to multi-layered geophysical structures. Using appropriate temporal variations for the field variables the time-related kernels are obtained explicitly. The BEM solutions presented for layered media are generated using synthetic seismograms and are seen to be stable. The seismograms obtained in terms of the BEM and FDM are compared.
PDF Abstract XML References Citation
How to cite this article
DOI: 10.3923/jas.2006.1703.1711
URL: https://scialert.net/abstract/?doi=jas.2006.1703.1711
INTRODUCTION
In recent decades, very great advances in computer science are closely linked to the growth of research in wave modelling. The wave modelling is a standard method for generating the wave response for an earth model. The use of wave equations to model wave response of complex geophysical materials with interfaces is a topic that has been intensively studied by seismologists. The challenge of this work is to build on satisfactory numerical models that are capable of solving the equations accurately and also cope with layered media.
A number of techniques are used for such modelling. Finite difference methods (FDM) (Reynolds, 1978; Demir, 1998; Phadke et al., 2000; Koparan, 2005), finite element methods (FEM) (Di et al., 2000), hybrid methods (Lie and Yu, 2002) are generally seen among these techniques in use. Well-known advantages of the Boundary Element Method (BEM) to tackle many problems in various disciplines (Beskos, 1997) lead us to use the BEM to solve the wave equation in large and complex geometries.
The general time domain BEM formulation for the scalar wave equation was established by Mansur (1983). Since then many BEM formulations have been employed to analyse various problems in the time domain (Israil and Banerjee, 1990; Wang and Takemiya, 1992; Carrer and Mansur, 1994; 2002; Mansur et al., 1998; Yu et al., 2000; Abreu et al., 2003; Marrero and Dominguez, 2003; Yokoi, 2003; Karabalis, 2004; Rizos and Zhou, 2006).
Acceptable BEM solutions depend upon numerical stability, as is the case for all numerical methods. This requirement led to the following important research:
• | Treatment of singular integrals, |
• | Temporal and spatial variation of the field variables, |
• | Use of time or frequency domain approaches, |
• | Consideration of direct or indirect BEM types. |
There has been a debate whether the equation of the BEM should employ the fundamental solutions in explicit or implicit form. Israil and Banerjee (1990) used the explicit method whilst Mansur (1983) and Dominguez (1993) adopted the implicit method. Carrer and Mansur (1999) also adopted the explicit approach for studying internal components in 2D-elastodynamic problems. The advantages and disadvantages of the earlier works were examined by Birgisson and Crouch (1998) in the context of the 2D-elastodynamic problems.
The present study applies the best (Gallego and Dominguez, 1996) and commonest scheme that is constant and linear time variations for the flux and potential, respectively. For the space approximation of the field variables constant elements are used. In the present study, the time-related kernel functions are considered explicitly.
Physically, the study focuses on the determination of the wave behaviour in layered media. Despite increased volume of research dedicated to the analysis of wave propagation by using the time domain BEM; consideration of layered, large and complex media in elasticity problems producing seismograms is not so common. Some of the examples are Ahmad and Banerjee, (1988) who use a frequency domain approach, Cheung et al. (1993) who use cylindrical co-ordinates for half-space geometries and Birgisson and Crouch (1998).
Consideration of the wave propagation problems in piecewise homogeneous media through the BEM is a profound problem due to its structural complexity. The equation is solved for various geophysical structures in the present study. The synthetic seismograms generated from the BEM solution of the wave equation for layered media are presented and a comparison of the BEM and FDM results gives qualitative agreement. To this end, some computer codes were produced for each subdomain as subprograms. The present study determines the wave behaviour through multi-layered media by substructuring and also by ensuring that the dynamical interfacial conditions are satisfied across the common interfaces.
BOUNDARY INTEGRAL EQUATION
The 2D scalar wave equation corresponding to a homogeneous isotropic elastic body Ω enclosed by the boundary B is given by Morse et al. (1953), as follows:
(1) |
In this equation φ, f and are functions of position and time and represent potential, body source and acceleration respectively, whilst c is the speed of wave propagation. In the above φ,ii and are the second derivatives of the potential φ with respect to the direction xi and the time t, respectively. Under consideration of zero body source and zero initial conditions, the scalar reciprocal equation in terms of two dynamical states may be given as:
(2) |
where, αi takes the values 1, 0.5 or 0 depending on whether the load point yi is in the domain, on the smooth boundary or outside the domain, respectively. Here the upper limit t+ is used to avoid ending the integration at the peak of the Dirac delta function (Morse and Feshbach, 1953). In the above equation:
(3) |
(4) |
where r is the distance between the load point yi and the field point x. In the above φ and φs represent the actual and fundamental solution states of the scalar potential, respectively, whilst H stands for the Heaviside function. Since the radiation conditions (Eringen and Suhubi, 1975) are automatically satisfied, the boundary integral Eq. 2 remains valid for semi-infinite (or infinite) domains, too. It may, therefore, be seen that there is no need to discretize an external boundary when it is of infinite extent. This makes the BEM advantageous over the domain techniques.
BOUNDARY ELEMENT EQUATION
In the study the boundary is approximated by straight-line elements. Temporal integration of the kernels of Eq. 2 is used to obtain a numerical solution of the partial differential equation by applying the boundary integral equations. The kernels explicitly obtained here for Eq. 2 agree with the kernels implicitly given by Dominguez (1993). The flux and potential are interpolated by constant and linear time variations, respectively. The choice of time interpolation functions is not arbitrary even if it is so in principle. However, when the piecewise linear time interpolation function is used for the flux as well as for the potential the solution process is prone to become unstable (Cole et al., 1978). Discussions have been going on the selection of the time functions and descriptions of the time functions used here can be found in, for example, Gallego and Dominguez (1996), Mansur et al. (1998), Sari (2000).
In this analysis time is divided into n equal intervals, t = nΔt. To obtain the approximate solution of the boundary integral equation, suppose that:
(5) |
where ηm(τ) and μm(τ) are temporal interpolation functions. Moreover φm and qm indicate the potential and flux, respectively, at time tm = mΔt at point x. The interpolation functions used by us explicitly are:
(6) |
(7) |
The discretized form of Eq. 2 is:
(8) |
where
(9) |
and:
(10) |
Using Eq. 7 and 9 the effect of the load at the field point x at time t = tn can be found as follows:
(11) |
where
During the evaluation of the time integral, three cases are categorised for expression (11) depending on the position of the wave front of the fundamental solution.
Similarly consideration of the integral (10) by using (6) gives the flux kernel explicitly:
(12) |
where
Use of the causality property classifies the expression (12) in four different cases depending on the position of the wave front of the fundamental solution.
It is worth noting that the flux kernel (12) has been obtained by considering the fundamental solution explicitly with the same time function ηm unlike Mansurs approach. In other words, the above flux kernel has been found without using further transformations (Mansur, 1983) on the Dirac delta function in the fundamental solution in Eq. 2. In addition, it is important to realise that there is no singularity appearing at the wave front in the potential and flux kernels in the above.
DISCRETIZATION
The boundary of the domain is discretized into a number of elements. Over each element, the co-ordinates
are expressed by means of their nodal values by using linear elements whilst the field variables are represented by constant elements. Then the geometry of element j can be represented by the co-ordinates:
(13) |
where
(14) |
The nodal values of the field variables on the boundary are approximated using the spatial interpolation function Ψj for the node j so that (5) becomes:
(15) |
where, φmj and qmj denote the potential and its normal derivative at node j for time tm = mΔt time whilst Ψj is the spatial interpolation function for the field variables. For the points x(ζ) inside an element, Ψj = Ψj(ζ). When the boundary nodal quantities are constant over the element in approximation (15) Ψj = 1. Using the spatial variations for the boundary node yi with a set of discrete elements Bj, j = 1, 2, , N on the boundary B, Eq. 8 can be written as:
(16) |
In Eq. 16 n shows the final time, t = nΔt, whilst φni denotes the unknown potential at the load point yi, at time step n.
INTEGRATION
The basic idea here is to solve Eq. 2 numerically by discretizing boundary values spatially as well as temporally. The boundary B is discretized to integrate spatially the kernels Unm and Qnm over the all boundary elements. Regular integrals are evaluated using a standard Gaussian quadrature. The integration to be evaluated is expressed by means of the homogeneous co-ordinate -1≤ζ≤1 along the elements. To evaluate the integrals the differential is expressed, in xy-plane, as:
(17) |
where, , J , is the Jacobian of the transformation. With the spatial discretization Eq. 16 takes the following form for the 2D scalar wave problems,
(18) |
The boundary is divided into N elements and the field variables are assumed to be constant over each element and equal to the value at the mid-element node.
A singularity exists if and only if the load and field points coincide at the first time step. Once the load point becomes also a field point the integrals are treated analytically (Dominguez, 1993).
After the regular and singular integrals have been evaluated, for each element j, one can write,
(19) |
(20) |
Use of the value enables Eq. 20 to be redefined as follows:
(21) |
Now Eq. 16 can be written as:
(22) |
where, N represents the total number of boundary elements. Taking all boundary elements, Eq. 22 can be rewritten in a more compact form:
(23) |
where, Gnm and Hnm are square matrices which are calculated by spatial integration for each element and φm and qm are the column vectors of boundary nodal quantities. At time t, there are as many unknowns as the number of equations in the matrix Eq. 23. If the boundary quantities φm and qm are known for the time m = 1, 2, , n-1, then for each time step n, the solution can be found. The detail of the solution procedure used here can be found in Dominguez (1993).
LAYERED MEDIA
The numerical implementation has mostly been considered for homogeneous materials. Geophysical structures are virtually infinitely variable, namely, they are inhomogeneous. In this study, two or more different homogeneous media will be considered. Figure 1 depicts a domain consisting of two subdomains which are Ω1 and Ω2. Applying the procedure for each subdomain separately gives, from Eq. 23,
(24) |
for subdomain Ω1.
Here, for time tm = mΔt, and indicate the nodal potentials and fluxes at the boundary B1, whilst and stand for the nodal potentials a nd fluxes at the interface when it belongs to B1. It may be seen that the corresponding coefficient matrices and are obtained from the interface when it belongs to the boundary B1. The matrices and are derived from the boundary B1, excluding the interface. Similarly
(25) |
for subdomain Ω2.
Where two subdomains have a common interface, it is necessary that the common nodal points on the boundaries B1 and B2 are at the same potential and that there is a dynamical equilibrium of the flux across the corresponding interface elements.
Fig. 1: | Two subdomains of the medium |
These conditions can be written as:
(26) |
(27) |
Note that for a subdomain the integral equation is constructed in terms of the interface variables of the subdomain as well as the boundary variables. At the interface both field variables are unknown, so that the number of unknowns and knowns are not equal. However considering Eq. 26 and 27 corresponding to the interface, the number of unknowns and knowns becomes equal. Thus, the discrete equation system is constructed for the piecewise homogeneous medium and simultaneously solved for all subdomains. By considering the interface equations it follows that the algebraic equation system can be written in a compact form as,
(28) |
where, Dn and Xn, at time step n, are the system matrix which consists of the matrix coefficients of both subdomains and the unknown vector which also contains the interface unknowns. The Dn matrix can be set up in two ways. In the first method, the Eq. 22 and 27 are added directly to the system matrix obtained by using the subdomains. In the second approach, the interface variables are eliminated. Here the first technique is used to consider piecewise homogeneous media. For the establishment of system matrix the first one is more flexible than the latter. On the other hand, in the latter, the dimension of the system matrix is smaller than the first.
In the coefficient matrix (28) interface contributions can be scaled. The argument in the above can be repeated for three or more subdomains. So for each time step the boundary and interface unknowns can be calculated from Eq. 28.
In terms of the boundary variables computed, for time t = Δt, 2Δt,..., nΔt the potential value at internal points selected can be found. From the boundary node x and internal point yi with the time difference n-m, internal potentials can be easily computed numerically at time step n.
NUMERICAL RESULTS
The behaviour of the scalar wave through 2D rock structures is examined. In the following examples, wave propagation in layered media due to excitations within a geophysical structure is modelled. The boundary conditions are homogeneous and inhomogeneous for external and internal boundaries, respectively.
Synthetic seismograms are computed for three different earth models. The seismograms correspond to the models have time (in seconds) on the horizontal axis. All displacement values received at the selected receiver points are on the vertical axis. The vertical axis shows the length of topside of the physical model in the seismograms and its unit is meters (m).
Although study has been going on the choice of the time variations, in this study the time variations used was the scheme given by Gallego and Dominguez (1996). The size of the time step affects the stability properties of the BEM system. A convenient factor for measuring stability is β = cΔt/Δx with time step size Δt and the length of element Δx, where 0.5≤β≤1.5, although for smaller or larger time values it is possible to obtain stable results. In the present study, the BEM results obtained are in this quoted interval. As was pointed out (Dominguez and Gallego, 1991) for β≥1 good results are obtained for straight-line constant elements.
TWO-LAYERED MEDIUM
This example concerns the study of solution of the wave Eq. (1) for a two-layered medium, with differing acoustic properties, of infinite extent at the sides and the bottom. The geometry of the problem is shown in Fig. 2 and is considered for comparison with the FDM results.
Fig. 2: | Physical model showing the geometry used to generate seismograms in the scalar BEM and FDM programs for a two-layered medium |
Fig. 3: | Synthetic seismogram generated from the BEM solution of the 2D scalar wave equation for the two-layered medium |
Fig. 4: | Synthetic seismogram generated from the FDM solution of the 2D scalar wave equation for the two-layered medium with transparent boundary conditions |
The velocity of the upper medium is 1500 m s-1 for seawater while the lower medium has a velocity of 2130 m s-1 for shale, as given by Reynolds (1978). As can be seen from Fig. 2 the receivers are accommodated in a horizontal line at 10 m below the top boundary. The Neumann boundary condition, is specified for the topside.
The number of receivers used is 120. The distance between the adjacent receiver points is equal. To obtain this synthetic seismogram the boundary elements used are uniform with 4 m element length and the total elapsed time is 0.7 s with 0.004 s the time increment. The source position is (180, 445). In the present study, the source is taken as an internal boundary, which is a small square. The Dirichlet boundary condition is prescribed for the sides of the small square with the length being 4 m. The results obtained from the BEM (Fig. 3) are in qualitative agreement on comparison to the FDM results (Fig. 4) obtained using the transparent boundary conditions to simulate a medium of infinite extent on the sides. Reflections have been following primary wave originates from the top boundary. It can also be noticed that the interface waves are visible. The first interface wave strikes the interface and then reflects. As for the second interface wave first hits the topside and then follows the same path as the first.
TWO-LAYERED MEDIUM (SALT DOME)
The salt dome geometry consists of two different media with the dome rising from the lower medium into the upper medium as can be seen from the physical geometry (Fig. 5). Figure 5 shows the problem by setting to the wave speeds 2500 m s-1 and 4500 m s-1. The source is located at (420,1790) and its size is 8 m.
Fig. 5: | The geometry of the salt dome model used to generate seismograms in the scalar BEM and FDM programs |
Fig. 6: | Synthetic seismogram generated from the BEM solution of the 2D scalar wave equation for the salt dome model |
Fig. 7: | Synthetic seismogram generated from the FDM solution of the 2D scalar wave equation for the salt dome model with transparent boundary conditions |
The element length, time increment and number of receivers used are 16 m, 0.01 s and 249, respectively. The Neumann boundary conditions are used for the topside of the salt dome geometry.
Figure 6 indicates the results obtained by using the BEM for the salt dome model. Figure 7 shows the FDM solutions obtained for the same model. Good agreement between the results can be seen. The energy of the incident wave is divided in two parts at the interface: going upward and going downward. Two interface waves from the dome can be seen from Fig. 6 and 7 and the last part of the first wave is followed by the second wave reflected from the dome. This wave is slightly clearer in the FDM results.
THREE-LAYERED MEDIUM
Three-layered medium is considered to model heterogeneous rocks. Three flat layers are used as is shown in Fig. 8. The velocities are given from the top layer to the bottom layer respectively: 1500 m s-1 of sea-water, 2440 m s-1 of marls and 4000 m s-1 of granite. The position of the source is (240,445) and its size is 5 m. The element size used is 5 m in this problem and the results obtained for the time increment of 0.0025 s is stable.
The results are depicted in a seismogram representation of the BEM solution. Figure 9 shows the results obtained by the BEM incorporating the interfaces of the medium. Two interface waves for each one can be easily seen from Fig. 9.
The BEM results in Fig. 9 and FDM results in Fig. 10 (Koparan, 2005) are in qualitative agreement except the last two weak waves caused by the third interface in the FDM. Relatively increase of the source size used in the BEM can redress this problem. In other words, according to our experience in the BEM, increase in the internal source size will be easily provide us to see the last interface waves as in the FDM results. However, this time the reflected waves from the sides of the source come out. Therefore, to clearly be able diagnose the waves and not to cause unnecessary waves, a suitable choice of the internal source size is important in multi-layered media. Arithmetic length can also cause this problem.
Fig. 8: | The geometry of a three-layered medium used to generate seismograms in the scalar BEM and FDM programs |
Fig. 9: | Synthetic seismogram generated from the BEM solution of the 2D scalar wave equation for the three-layered medium |
Fig. 10: | Synthetic seismogram generated from the FDM solution of the 2D scalar wave equation for the three-layered medium with transparent boundary conditions |
Yet, length arithmetic used is not the same in this problem. So those, in the FDM results, numerical errors clearly come out. As will be pointed out later on, different discretizations of the two methods also pay an important role.
These layered media problems also show the versatility of the BEM.
All the results presented were depicted in MATHCAD 2001 environment except Fig. 10 obtained with the use of MATLAB 6.5.
CONCLUSIONS AND FUTURE WORK
In the present study the time-domain 2D scalar wave BEM has been extended to multi-layered media structures by adding the interfacial equations. Since, in nature, geomaterials are heterogeneous requiring much intricacy through dealing with the BEM, it is believed that the layered media results presented here are important.
Although the existence of the interfaces makes the solution of the problem more prone to instability, the present study showed that the BEM is capable of treating layered media problems. It has been recognized that the method needs to be used carefully, especially for piecewise homogeneous media, because of structural intricacy and because of the possibility of the emergence of the stability problems.
It is also concluded that to increase the numerical stability of the results for the time domain BEM, an increase in the number of elements is suggested. However it should not be forgotten that if the elements size is taken to be very small then the desired stability might not be obtained. Despite the obtained qualitative agreement in the two methods, there are some quantitative differences between the two results in general. The most striking reasons of the differences can be given by:
• | The two methods are based on different discretization. |
• | For the two techniques the source used is different. |
Moreover, comparison between the two methods clearly states the power of the BEM over the FDM.
The presented results give confidence for future work where consideration may be given to the problem of 2D elastic wave propagation which are similarly based but more involved.
REFERENCES
- Abreu, A.I. and J.A.M. Carrer and W.J. Mansur, 2003. Scalar wave propagation in 2D: A BEM formulation based on the operational quadrature method. Eng. Anal. Boundary Elem., 27: 101-105.
Direct Link - Ahmad, S. and P.K. Banerjee, 1988. Time domain transient elastodynamic analysis of 3D solids by BEM. Int. J. Num. Meth. Eng., 26: 1709-1728.
Direct Link - Beskos, D.E., 1997. Boundary element methods in dynamic analysis: Part II (1986-1996). Applied Mech. Rev., 50: 149-197.
CrossRefDirect Link - Carrer, J.A.M. and W.J. Mansur, 1994. Space derivatives in the time domain BEM analysis for the scalar wave equation. Eng. Anal. Boundary Elem., 13: 67-74.
Direct Link - Carrer, J.A.M. and W.J. Mansur, 2002. Time dependent fundamental solution generated by a not impulsive source in the boundary element method analysis of the 2D scalar wave equation. Commun. Num. Methods Eng., 18: 277-285.
Direct Link - Cheung, Y.K., L.G. Tham and Z.X. Lei, 1993. Wave propagation in layered media by time domain BEM. Earthquake Eng. Struct. Dyn., 22: 225-244.
Direct Link - Cole, D.M., D.D. Kosloff and J.B. Minster, 1978. A numerical boundary integral equation method for elastodynamics I. Bull. Seismo. Soc. Am., 68: 1331-1357.
Direct Link - Di, Q., L. Zhu and M. Wang, 2000. 2-D finite element modeling for seismic wave response in media with sand bodies. Phys. Earth Planetary Interiors, 120: 245-254.
Direct Link - Gallego, R. and J. Dominguez, 1996. Hypersingular BEM for transient elastodynamics. Int. J. Num. Meth. Eng., 39: 1681-1705.
Direct Link - Israil, A.S.M. and P.K. Banerjee, 1990. Advanced development of time-domain BEM for 2D scalar wave propagation. Int. J. Num. Meth. Eng., 29: 1003-1020.
Direct Link - Karabalis, D.L., 2004. Non-singular time domain BEM with applications to 3D inertial soil-structure interaction. Soil Dyn. Earthquake Eng., 24: 281-293.
Direct Link - Lie, S.T. and G.Y. Yu, 2002. Stability improvement to BEM/FEM coupling scheme for 2D scalar wave problems. Adv. Eng. Software, 33: 17-26.
Direct Link - Mansur, W.J., J.A.M. Carrer and E.F.N. Siqueira, 1998. Time discontinuous linear traction approximation in time-domain BEM scalar wave propagation analysis. Int. J. Num. Meth. Eng., 42: 667-683.
Direct Link - Marrero, M. and J. Dominguez, 2003. Numerical behaviour of the time domain BEM for three-dimensional transient elastodynamic problems. Eng. Anal. Boundary Elem., 27: 39-48.
Direct Link - Phadke, S., D. Bhardwaj and S.K. Dey, 2000. An explicit predictor-corrector solver with application to seismic wave modelling. Comput. Geosci., 26: 1053-1058.
Direct Link - Rizos, C.D. and S. Zhou, 2006. An advanced direct time domain BEM for 3-D wave propagation in acoustic wave. J. Sound Vibrat., 293: 196-212.
Direct Link - Yokoi, T., 2003. The higher order born approximation applied to improve the solution seismic response of a three-dimensional canyon computed by the indirect boundary eleament method. Phys. Earth Planetary Interiors, 137: 97-106.
Direct Link - Wang, C. and H. Takemiya, 1992. Analytical elements of time domain BEM for two-dimensional scalar wave problems. Int. J. Num. Meth. Eng., 33: 1737-1754.
Direct Link - Yu, G., W.J. Mansur, J.A.M. Carrer and L. Gong, 2000. Stability of Galerkin and collocation time domain boundary element methods as applied to the scalar wave equation. Comput. Struct., 74: 495-506.
Direct Link