Modified integration rules for reducing dispersion error in ... - CiteSeerX
Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287 www.elsevier.com/locate/cma
Modified integration rules for reducing dispersion error in finite element methods Murthy N. Guddati *, Bin Yue Department of Civil Engineering, North Carolina State University, Campus Box 7908, Raleigh, NC 27695-7908, USA Received 10 July 2003; accepted 17 September 2003
Abstract This paper describes a simple but effective technique for reducing dispersion errors in finite element solutions of timeharmonic wave propagation problems. The method involves a simple shift of the integration points to locations away from conventional Gauss or Gauss–Lobatto integration points. For bilinear rectangular elements, such a shift results in fourth-order accuracy with respect to dispersion error (error in wavelength), as opposed to the second-order accuracy resulting from conventional integration. Numerical experiments involving distorted meshes indicate that the method has superior performance in comparison with other dispersion reducing finite elements. 2003 Elsevier B.V. All rights reserved. Keywords: Wave propagation; Numerical dispersion; Numerical integration; Helmholtz equation
1. Introduction Finite element methods, when applied to time-harmonic wave propagation problems, incur solution errors that can be classified into error in the amplitude and error in the wavelength. In large-scale problems where the domain size is much larger than the wavelength, the error in the wavelength, known as the dispersion error, tends to accumulate and produce large phase errors, resulting in completely erroneous results [1]. On the other hand, amplitude error is relatively benign in that it does not increase with growing domain size. This paper focuses on the reduction of dispersion error. The dispersion error could be reduced either by using a finer mesh, or by using higher order finite elements. Both of these options come with large increase in computational cost, making some of the largescale simulations prohibitively expensive. It is thus desirable to reduce the dispersion for lower order finite elements. It can be quickly realized that dispersion error in one-dimensional problems can be completely eliminated by resorting to analytical solutions [1]. For higher dimensions, the dispersion error cannot be completely eliminated [2], but can be reduced. Over the past two decades, several researchers have proposed
*
Corresponding author. Tel.: +1-919-515-7699; fax: +1-919-515-7908. E-mail address: [email protected] (M.N. Guddati).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.09.010
276
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
various methods for reducing the dispersion in higher dimensions, especially for low-order finite elements. The following is a brief survey. The peculiar aspect of dispersion error in higher dimensions is its anisotropy, i.e., the dispersion error depends on the angle between the wave propagation direction and mesh orientation. Mullen and Belytchko [3] have analyzed this property, and compared the performance of various finite element as well as standard finite difference techniques. Marfurt [4] reduced the dispersion error by taking generalized average of consistent and lumped mass matrices, but the anisotropy still remains. Later, Krenk [5] proposed selective scaling of stiffness and mass matrices to reduce the dispersion error. The implication of his method for distorted meshes is not clear. More recently, Harari and Hughes [6] have used Galerkin/least squares (GLS) method [7] to reduce dispersion on square meshes. The resulting dispersion error, although reduced, is still significant. Thompson and Pinsky [8] modified this method to further reduce the dispersion, but their method still suffers from anisotropy. Babuska and Ihlenburg [9] used the concepts of generalized finite element method (GFEM) and finite difference stencils to derive a method that reduces the dispersion by several orders of magnitude, but only for square meshes. Oberai and Pinsky [10] used the concept of edge residuals to derive the same scheme in the context of finite elements. However, when extended to rectangular or unstructured meshes, their scheme becomes heuristic in nature and loses its accuracy properties. Babuska and Melenk [11] used a partition of unity method to reduce the dispersion, but for the method to be practical, one needs to know the behavior of the solution a priori, which is often not the case. The method of residual free bubbles, developed by Franca et al. [12], is a robust method that can improve the accuracy by effective enrichment of shape functions, but such enrichment demands significant increase in computational cost, especially for non-uniform meshes. Similar comment applies to the sub-grid modeling approach proposed by Cipolla [13]. Weighted-averaging FEM proposed by Min et al. [14] significantly reduces dispersion, but only on square meshes and at significant increase in computational cost. A fundamental limitation of most of the above methods is that they are applicable only for propagating waves in homogeneous media discretized by regular square meshes. Many important and practical problems of wave propagation involve complex geometries and material heterogeneities, and cannot be discretized using square meshes. Although dispersion reducing methods for triangular meshes may offer added flexibility [15], the current versions are limited to uniform structured meshes. With the aim of reducing dispersion on unstructured meshes, we recently proposed a dispersion-reducing framework named the Local Mesh-Dependent Augmented Galerkin (L-MAG) methods [16]. L-MAG methods provide a systematic procedure to reduce the dispersion error on rectangular meshes. The methods are more flexible than the other dispersion reducing methods in that they are accurate for meshes containing rectangular elements of varying size. However, although they perform better than other dispersion reducing methods, L-MAG methods lose their accuracy properties when used for distorted meshes. In this paper, we develop an effective dispersion reducing technique that is significantly simpler than all the existing dispersion reducing methods. The method is based on the simple observation that by shifting the Gauss quadrature to some unconventional locations would reduce the anisotropic dispersion significantly. In fact, for rectangular elements, the proposed method makes the numerical wavelength fourth-order accurate, as opposed to the second-order accuracy obtained from conventional integration rules. Needless to say, the proposed method can be directly applied to distorted meshes. The remainder of the paper focuses on the derivation of the method, its application to structured and unstructured meshes and the comparison of its performance with some existing dispersion reducing finite element methods. The following is the outline of the paper. Section 2 contains discussion of the variational boundary value problem and the Galerkin finite element discretization of the problem. The idea of shifting the integration points is discussed in Section 3. The locations of the integration points that reduce the dispersion error are obtained in Section 4. Section 5 focuses on evaluating the performance of proposed method with the help of numerical examples. The paper is concluded in Section 6 with some closing remarks.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
277
2. Finite element formulation of time-harmonic wave propagation problems Helmholtz equation, also called the reduced wave equation, is used to analyze scalar wave propagation encountered in linear acoustic and anti-plane shear problems. The corresponding boundary value problem takes the form: find u : X ! C, such that r2 u k 2 u ¼ f
ð1Þ
in X
with boundary conditions u¼g
ð2Þ
on Cg ;
ru n þ bðuÞ ¼ F
ð3Þ
on Ch :
In the above, u is the complex-valued field variable (acoustic pressure, or anti-plane displacement), k is the wave number associated with the excitation frequency (k ¼ x=c), g is the specified field variable on the Dirichlet boundary Cg , whereas a Robin type boundary condition is specified on the remainder of the boundary Ch . Here, b is an operator ranging from zero (for pure Neumann boundary condition) to the DtN map in the context of absorbing boundaries. Using classical variational calculus, the above differential (strong) form of the boundary value problem can be converted to the following variational (weak) form: find u : X ! C, with u ¼ g on Cg , such that ðrv; ruÞX k 2 ðv; uÞX þ ðv; bðuÞÞCh ¼ ðv; f ÞX þ ðv; F ÞCh
8v; v ¼ 0 on Cg
The notation ð; ÞX is used to define inner product, over the domain X, as follows: Z ðf ; gÞX ¼ f g dX;
ð4Þ
ð5Þ
X
where the f is the complex conjugate of f . The notation ð; ÞCh represents inner product taken over the boundary Ch . For numerical solution of the above variational boundary value problem, (Bubnov) Galerkin FEM can be used, where the Galerkin approximation of the variational boundary value problem takes the form: find uh : X ! C, with uh ¼ gh on Cg , such that ðrvh ; ruh ÞX k 2 ðvh ; uh ÞX þ ðvh ; bðuh ÞÞCh ¼ ðvh ; f ÞX þ ðvh ; F ÞCh
8vh ; vh ¼ 0 on Cg :
ð6Þ
In the above, uh and vh are approximations of u and v, respectively, and gh is the consistent approximation of the Dirichlet data. In the context of conventional finite element method, uh and vh take the form uh ðxÞ ¼ NðxÞU and
vh ðxÞ ¼ NðxÞV;
ð7Þ
where NðxÞ is the shape function (row) vector, while U is the discretized field variable (degree of freedom vector) and V is its variation. By substituting (7) in (6), we obtain the algebraic equation SU ¼ F;
ð8Þ
where S is the coefficient matrix given by S ¼ K k 2 M þ ðN; BðNÞÞCh
ð9Þ
and F is the force vector F ¼ ðN; f ÞX þ ðN; F ÞCh :
ð10Þ
In (9), note that B is the discretized version of the operator b in (3) and K ¼ ðrN; rNÞX ;
M ¼ ðN; NÞX
ð11Þ
278
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
are stiffness and mass matrices, respectively. These system matrices are obtained by assembling the element matrices that are given by Z Z T ðrNÞ rN dX; Me ¼ NT N dX; ð12Þ Ke ¼ Xe
Xe
where Xe is the domain associated with the element.
3. Generalized integration rules for evaluation of stiffness and mass matrices Since finite elements can be distorted, it is customary to map the finite elements from a square parent element (f1 < n < þ1g f1 < g < þ1g), and transform the above integrals into the natural coordinate system as follows: Z þ1 Z þ1 Z þ1 Z þ1 T Ke ¼ ðrNÞ rN detðJÞ dn dg; Me ¼ NT N detðJÞ dn dg: ð13Þ 1
1
1
1
The above integrals are traditionally evaluated using Gauss quadrature Z þ1 Z þ1 ng X / dn dg ¼ Wi /ðni ; gi Þ; 1
1
ð14Þ
i¼1
where ng is the number of Gauss points, Wi are the weights and ðni ; gi Þ are the local coordinates of Gauss pffiffiffi points. For pffiffiffi the most commonly used bilinear quadrilateral element, a 2 · 2 Gauss rule (ni ¼ 1= 3; gi ¼ 1= 3; Wi ¼ 1) is used to evaluate the stiffness matrix. Such an integration rule is considered desirable, as it represents exact integration for square and rectangular elements. On the other hand, the mass matrix is pffiffiffi pffiffiffi evaluated either by 2 · 2 Gauss rule (ni ¼ 1= 3; gi ¼ 1= 3; Wi ¼ 1) resulting in the consistent mass matrix, or Gauss–Lobatto rule (ni ¼ 1; gi ¼ 1; Wi ¼ 1) resulting in lumped mass matrix. Lumped mass matrix is preferred for explicit computations in time-domain, while there is no special preference for frequency domain solutions. Unlike consistent integration, lumped mass evaluation incurs errors in numerical integration. However, since the integration errors are of the same (second) order as the solution error, lumped mass is considered an acceptable alternative to consistent mass matrices. These observations indicate that any integration rule that has second-order accuracy would be a viable alternative to the conventional integration schemes. Here we introduce such a generalized 2 · 2 integration rule as given below Ke ¼
4 X
T
Wi K ðrNðnKi ; gKi ÞÞ rNðnKi ; gKi Þ detðJÞ
ð15Þ
i¼1
and Me ¼
4 X
M M M Wi M NðnM i ; gi ÞNðni ; gi Þ detðJÞ:
ð16Þ
i¼1
Considering the symmetry of the element, we propose to use symmetric, yet general, locations of integration M points with nKi ¼ a and gKi ¼ a for the stiffness terms, and nM i ¼ b and gi ¼ b for the mass terms, K M with the weights Wi ¼ Wi ¼ 1. As will be shown in the next section, the above generalization gives us the flexibility to reduce the dispersion error through appropriate choice of the location of integration points.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
279
4. Dispersion reducing integration rule In this section, by taking advantage of the flexibility offered by the generalized integration rule, we try to fix the integration points so that the dispersion error is minimized for rectangular meshes. The coefficient matrix (S ¼ K k 2 M) of a four-node rectangular finite element of size Dx Dy, when computed using Eqs. (15) and (16), is given by 2 3 S0 Sx Sxy Sy 6 Sx S0 Sy Sxy 7 7 ð17Þ Se ¼ 6 4 Sxy Sy S0 Sx 5; Sy Sxy Sx S0 where ðDx2 þ Dy 2 Þða2 þ 1Þ DxDy 2 ð1 þ b2 Þ k 2 ; 4DxDy 16 ðDx2 Dy 2 Þ a2 ðDx2 þ Dy 2 Þ DxDy ð1 b4 Þk 2 ; Sx ¼ 4DxDy 16 ðDy 2 Dx2 Þ a2 ðDx2 þ Dy 2 Þ DxDy ð1 b4 Þk 2 ; Sy ¼ 4DxDy 16 ðDx2 þ Dy 2 Þða2 1Þ DxDy 2 ð1 b2 Þ k 2 : Sxy ¼ 4DxDy 16
S0 ¼
ð18Þ
In order to obtain the equations associated with the central node (x; y) of the finite element patch in Fig. 1, we need to assemble all four elements connected to the node. Such an assembly will result in the following form of the equation: 4S0 uhx;y þ 2Sx ðuhxDx;y þ uhxþDx;y Þ þ 2Sy ðuhx;yDy þ uhx;yþDy Þ þ Sxy ðuhxDx;yDy þ uhxþDx;yDy þ uhxþDx;yþDy þ uhxDx;yþDy Þ ¼ 0:
ð19Þ
This equation can also be obtained by an equivalent finite difference method using a nine-node stencil shown in the Fig. 1. Noting that (19) is applicable for all nodes of the discretized system, it can be shown that the solution for the difference equations takes the form of a plane wave uh ¼ eik
h ðx cos hþy cos hÞ
ð20Þ
:
4
∆x
3 ∆y
(x
i +1
(x
, y j −1 )
i −1
j +1
i
, yj )
i −1
, y j− 1 )
(x , y ) i
j −1
(x
i +1
(x , y ) i
(x
1
(x , y )
, y j +1 )
(x
, yj )
(x
, y j−1)
i +1
2
j
i +1
Fig. 1. Nine-node four-element rectangular finite element patch for dispersion analysis (nine-node stencil in case of finite differences).
280
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
The only difference between the approximate solution and the exact solution (u ¼ eikðx cos hþy sin hÞ ) is that the wave numbers, k h and k, are different. In fact, the relative dispersion error can be quantified as the relative error in wavelengths k kh k h k ¼ ð21Þ erelative ¼ dispersion k k h : In the above equations, the approximate wave number k h can be obtained by substituting the approximate waveform in (20) into the difference equation (19). Such substitution, after some simplification, results in ½S0 þ Sx Cx þ Sy Cy þ Sxy Cx Cy uhx;y ¼ 0;
ð22Þ
where Cx ¼ cosðk h Dx cos hÞ;
Cy ¼ cosðk h Dy sin hÞ:
ð23Þ
The expression in the brackets in (22) must be zero for the non-trivial solution (ux;y 6¼ 0). Thus the discrete dispersion relationship relating the exact and approximate wave numbers is given by S0 þ Sx Cx þ Sy Cy þ Sxy Cx Cy ¼ 0: Using (18), the dispersion relation can be written as 4 ðDx2 þ Dy 2 Þða2 þ Cx Cy a2 Cx a2 Cy a2 þ 1 Cx Cy Þ þ ðDx2 Dy 2 ÞðCx Cy Þ h i Dx2 Dy 2 k 2 ð1 þ b2 Þ2 þ Cx Cy ð1 b2 Þ2 þ ðCx þ Cy Þð1 b4 Þ ¼ 0
ð24Þ
ð25Þ
and the relation between the exact wave number k and the approximated wave number k h is obtained as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðDx2 þ Dy 2 Þða2 þ Cx Cy a2 Cx a2 Cy a2 þ 1 Cx Cy Þ þ ðDx2 Dy 2 ÞðCx Cy Þ k¼ : ð26Þ 2 2 DxDy ð1 þ b2 Þ þ Cx Cy ð1 b2 Þ þ ðCx þ Cy Þð1 b4 Þ In order to obtain the relative dispersion error, we take Taylor expansion of k in terms of k h . After some simplifications, the dispersion error can be written as k kh k2 2 2 2 2 2 2 2 2 2 2 k h ¼ 24 ð3a 2ÞðDx þ Dy Þ cos h sin h þ ð2 3b ÞðDx cos h þ Dy sin hÞ þ Oðk 4 ðDx þ DyÞ4 Þ:
ð27Þ
Using the above equation, we can obtain the estimates pffiffiffiof the dispersion errors for classical integration rules. When consistent integration is used (a ¼ b ¼ 1= 3), the dispersion error becomes k kh k2 2 4 4 4 2 4 ð28Þ k h ¼ 24 Dx cos h þ Dy sin h þ O k ðDx þ DyÞ : pffiffiffi When lumped mass is used (a ¼ 1= 3; b ¼ 1) we obtain an analogous result k kh k2 ¼ ðDx2 þ Dy 2 Þ cos2 h sin2 h þ ðDx2 cos2 h þ Dy 2 sin2 hÞ þ O k 4 ðDx þ DyÞ4 : ð29Þ k h 24 From the above two estimates, it is clear that conventional integration rules result in second-order accuracy of the numerical wave number. theffi other hand, On pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi ffi we observe that a simple shift of the integration points to unconventional location of 2=3; 2=3 , i.e.,
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
2 a2 ¼ b2 ¼ ; 3
281
ð30Þ
eliminates the second-order term in the error (27). It can be shown that, for this case, the dispersion error is given by k kh k 4 4 6 6 6 4 6 ð31Þ k h ¼ 480 Dx cos h þ Dy sin h þ Oðk ðDx þ DyÞ Þ: Note that the resulting dispersion error is similar to that of L-MAG2 [16], but the proposed method is much simpler. pffiffiffiffiffiffiffiffi The choice of b ¼ 2=3 for mass integration is not completely unconventional. For rectangular elements, it can be shown to be equivalent to averaging lumped and consistent mass matrices (see e.g. [4]). pffiffiffiffiffiffiffi ffi However, the choice of a ¼ 2=3 is unconventional and is the key to reducing the anisotropy in numerical dispersion. Incorporating the proposed modified integration rule into finite element software is trivial and imposes no additional computational cost. The proposed modification is completely local (element-level). Furthermore, unlike other dispersion reducing methods, the implementation of this method is independent of the element distortion. These observations indicate that the method is not only effective for uniform rectangular meshes, but could also work well for non-uniform distorted meshes. Further it is also straightforward to apply the proposed method to three-dimensional problems.
5. Numerical examples The performance of the proposed modified integration rule is tested using three examples: simulation of plane wave, simulation of radiation from a point source, and simulation of asymmetric radiation from a circular cavity. All the problems are solved using the proposed method, as well as some existing methods, namely Galerkin method with consistent integration (referred to as Galerkin in the figures), modified Galerkin least squares method [8] (GLSm), generalized finite element [9] or residual-based method [10] (Residual), and the second- and sixth-order local mesh-dependant Galerkin (L-MAG) methods proposed recently by the authors [16] (L-MAG2 and L-MAG6). 5.1. A plane wave problem We consider a plane wave u ¼ eikðx cos hþy sin hÞ with k ¼ 75 propagating in a square domain. Two types of meshes (uniform and distorted rectangular elements, Fig. 2) are used to discretize the domain. Robin boundary conditions are used for the analyses. Since the dispersion error accumulates whereas amplitude error does not, several simulations with increasing domain sizes are performed to differentiate the two types of errors. The associated errors are measured in terms of the relative L2 error given by e¼
ku uh kL2 : kukL2
ð32Þ
5.1.1. Uniform rectangular mesh Five domain sizes, 0.2 · 0.2, 0.5 · 0.5, 1 · 1, 1.5 · 1.5 and 2 · 2 are considered. The domains are discretized using 0.01 · 0.005 rectangular elements implying that kDx ¼ 0:75 and kDy ¼ 0:375 for each case. The simulation is performed for nine different propagation angles ranging from 0 to 90. For the domain size of 1 · 1, Fig. 3 shows comparison of errors resulting from various methods. Errors for other domain sizes have
282
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
Fig. 2. Meshes used in numerical examples: (a) uniform rectangle element and (b) distorted rectangular element.
10 Galerkin
GLSm
Residual
L-MAG2
L-MAG6
Proposed
Relative error
1
0.1
0.01
0.001 0
11.25
22.5
33.75
45
56.25
67.5
78.75
90
Propagation direction (degrees)
Fig. 3. Error for varying angles of propagation (with uniform rectangular elements).
2 Galerkin GLSm Relative error
1.5
Residual L-MAG2 L-MAG6
1
Proposed
0.5
0 0.2
0.5
0.8
1.1 Domain size
1.4
1.7
2
Fig. 4. Error for varying domain sizes (with uniform rectangular elements).
similar variation. It can be clearly seen that L-MAG and proposed methods have the best performance, with their errors much less than the errors from Galerkin, GLSm and Residual/GFEM methods. In order to illustrate the effect of the reduction in dispersion error, the error for various domain sizes is plotted with propagation angle fixed at 0. Fig. 4 shows the increasing error as the domain size is increased,
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
283
which is essentially due to the accumulative nature of the dispersion error. It can be clearly seen that the error increases significantly for Galerkin, GLSm and residual/GFEM methods, while the error growth is minimal for L-MAG methods as well as the proposed method. This clearly illustrates that the proposed method, as well as L-MAG methods, can be used to solve large-scale wave propagation problems with rectangular meshes. 5.1.2. Distorted rectangular mesh The distorted mesh is made up of rectangular patches consisting of four distorted elements as shown in Fig. 5. Except for the element distortion, all the parameters are identical to the case of uniform rectangular meshes. The errors from various methods with different propagation directions for the unit domain are shown in Fig. 6. Errors for other domain sizes have similar variation. The increasing error as the domain size is increased is plotted in Fig. 7. It is clear that the L-MAG methods have smaller errors than other existing methods. However, the striking feature is that the proposed method even outperforms the L-MAG methods, indicating that the proposed method is best suited for meshes when there is distortion. 5.2. Radiation from point source The time harmonic radiation from a point source in infinite domain is considered. The computational domain is a unit square with unit pulse at the lower left corner. The wave number is chosen as 75. The
Fig. 5. Layout of the distorted mesh (h ¼ 0:005).
10 Galerkin
GLSm
Residual
L-MAG2
L-MAG6
Proposed
Relative error
1
0.1
0.01
0.001 0
11.25
22.5 33.75 45 56.25 67.5 Propagation direction (degrees)
78.75
90
Fig. 6. Error for varying angles of propagation (with distorted rectangular elements).
284
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287 2 Galerkin GLSm Relative error
1.5
Residual L-MAG2
1
L-MAG6 Proposed
0.5
0 0.2
0.5
0.8
1.1 Domain size
1.4
1.7
2
Fig. 7. Error for varying domain sizes (with distorted rectangular elements).
Real part of displacement
0.2 0.15
Galerkin
GLSm
Residual
L-MAG6
Proposed
Exact
L-MAG2
0.1 0.05 0 0.90
0.92
0.94
0.96
0.98
1.00
-0.05 -0.1 Distance from center (r)
Fig. 8. Radiation from point source (with uniform rectangular elements): response along the radius.
meshes used in the previous example are used for this example, i.e. elements have aspect ratio of 2, with khx ¼ 0:75, khy ¼ 0:375 (see Fig. 2). A PML type absorbing boundary condition is used to simulate wave radiation into the exterior [17]. The results from various methods are compared with the exact solution, which is similar to the GreenÕs function of the two-dimensional HelmholtzÕs operator [18], and is given by ð1Þ
uðrÞ ¼ iH0 ðkrÞ;
ð33Þ
ð1Þ
where H0 is the zeroth-order Hankel function of the first kind. The real parts of the response along h ¼ 0 with r ranging from 0.9 to 1.0 are shown in Figs. 8 and 9 for uniform and distorted rectangular meshes, respectively. It can be seen that, consistent to the previous observations, L-MAG methods and proposed method have superior accuracy properties for rectangular meshes, while the proposed method outperforms all other methods for distorted meshes. 5.3. Radiation from an infinitely long circular cylinder This example considers the radiation from an infinitely long, circular cylinder (R ¼ 1) into an unbounded fluid domain. It reduces to a two-dimensional problem when the prescribed pressure on the wet surface does not vary along the axis of the cylinder. The wave number for this problem is chosen to be 4p. The Dirichlet boundary condition on the surface of the cylinder (r ¼ R) is chosen as u ¼ cos h.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
285
Fig. 9. Radiation from point source (with distorted rectangular elements): response along the radius.
5
4
3
2
1
0
0
1
2
3
4
5
Fig. 10. Radiation from an infinitely long circular cylinder: finite element mesh.
The unbounded domain is truncated using a 10 · 10 square around the cavity, with PML absorbers applied at the computational boundary. Furthermore, due to symmetry, only one-quarter of the problem is analyzed. The mesh is shown in Fig. 10. Sixty elements are used along the radial direction and 100 in circumferential direction. The mesh consists of elements with aspect ratio as high as 6 (see inset A), kh as high as 1.9 (see inset B), indicating only 3–4 elements per wavelength. Furthermore, most elements are distorted. Various methods are used to analyze this problem and the obtained solution is compared with the exact solution given by ð1Þ
uðr; hÞ ¼ ð1Þ
H1 ðkrÞ ð1Þ H1 ðkRÞ
cos h;
ð34Þ
where H1 is the first-order cylindrical Hankel function of the first kind. The real parts of results along h ¼ 0 with r ranging from 4.5 to 5.0 are shown in Fig. 11. It is observed that the dispersion errors from LMAG and proposed methods are much smaller than other methods, while the results from L-MAG6 and the proposed method are the closest to the exact solution.
286
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287 0.8
Real part of displacement
0.6
Galerkin
GLSm
Residual
L-MAG6
Proposed
Exact
L-MAG2
0.4 0.2 0 4.5
4.6
4.7
4.8
4.9
5.0
-0.2 -0.4 -0.6 Distance from center (r)
Fig. 11. Radiation from source point: response along the radius.
6. Concluding remarks In this paper, we developed a simple but effective method for reducing numerical dispersion in finite element solutions of time-harmonic wave pffiffiffiffiffiffiffipropagation ffi pffiffiffiffiffiffiffiffi problems. By simply shifting the integration points to the unconventional location of ð 2=3; 2=3Þ, we were able to obtain fourth-order accuracy in the wave number, as opposed to the second-order accuracy resulting from conventional methods. In spite of its simplicity, the proposed method appears to outperform more involved dispersion reducing methods, especially for rectangular and distorted elements. Needless to say, the implementation of the proposed method into existing finite element programs is a trivial issue and there is no additional computational cost associated with the method. It should be noted that when only square elements are used, the GFEM/residual methods, as well as L-MAG6, would have the best performance, resulting in sixth-order accuracy in the wavelength. When the mesh contains distorted and rectangular elements, GFEM/residual methods quickly lose their accuracy. L-MAG6 method maintains its accuracy for rectangular meshes, but deteriorates when applied to distorted meshes. On the other hand, the proposed method retains accuracy even for distorted meshes, while being almost as accurate as L-MAG6 for rectangular meshes. Based on these observations, we advocate the use of the proposed modified integration rule whenever the elements are (even slightly) distorted, L-MAG 6 when rectangular, L-MAG6/GFEM/residual methods, when the mesh contains only square elements. The proposed method, although developed for two-dimensional problems, is directly extendible to threedimensional problems. It may be possible to further increase the accuracy of the proposed method by combining the proposed idea with the ideas of L-MAG methods. It may also be possible to devise modified integration rules for higher order finite elements, as well as for finite elements for more complex waves such as elastic and flexural waves. These extensions are subjects of ongoing investigations.
Acknowledgements This material is based upon work supported by the National Science Foundation under Grant No. 0100188. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
287
References [1] F. Ihlenburg, in: Finite Element Analysis of Acoustic Scattering, Applied Mathematical Sciences, vol. 132, Springer-Verlag, New York, 1998. [2] I.M. Babuska, S.A. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, S.I.A.M. J. Numer. Anal. 34 (1997) 2392–2423. [3] R. Mullen, T. Belytschko, Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation, Int. J. Numer. Methods Engrg. 18 (1982) 11–29. [4] K.J. Marfurt, Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation, Geophysics 49 (1984) 533–549. [5] S. Krenk, Optimal formulation of simple finite elements, in: Variational Methods in Engineering, Proceedings of the 2nd International Conference, Springer, Berlin, 1985. [6] I. Harari, T.J.R. Hughes, Galerkin/least squares finite element method for the reduced wave equation with non-reflecting boundary conditions, Comput. Methods Appl. Mech. Engrg. 98 (1992) 411–454. [7] T.J.R. Hughes, Multiscale phenomena: GreenÕs functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles, and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. [8] L.L. Thompson, P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38 (1995) 371–397. [9] I.M. Babuska, F. Ihlenburg, A generalized finite element method for solving the Helmholtz equation in two-dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128 (1995) 325–359. [10] A.A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 49 (2000) 399–419. [11] I.M. Babuska, J.M. Melenk, The partition of unity finite element method, Int. J. Numer. Methods Engrg. 40 (1997) 727–758. [12] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residue-free bubbles for the Helmholtz equation, Int. J. Numer. Methods Engrg. 40 (1997) 4003–4009. [13] J.L. Cipolla, Subgrid modeling in a Galerkin method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 177 (1999) 35–49. [14] D.J. Min, H.S. Yoo, C. Shin, H.J. Hyun, J.H. Suh, Weighted-averaging finite-element method for scalar wave equation in the frequency domain, J. Seism. Explor. 11 (2002) 197–222. [15] I. Harari, C.L. Nogueira, Reducing dispersion of linear triangular elements for the Helmholtz equation, J. Engrg. Mech. ASCE 128 (2002) 351–358. [16] M.N. Guddati, B. Yue, Local mesh-dependent augmented Galerkin methods for reducing dispersion error for Helmholtz equation, Int. J. Numer. Methods Engrg. (submitted). [17] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Computat. Phys. 114 (1994) 185–200. [18] P.M. Morse, K.U. Ingard, Theoretical Acoustics, Princeton University Press, Princeton, NJ, 1986.
Modified integration rules for reducing dispersion error in finite element methods Murthy N. Guddati *, Bin Yue Department of Civil Engineering, North Carolina State University, Campus Box 7908, Raleigh, NC 27695-7908, USA Received 10 July 2003; accepted 17 September 2003
Abstract This paper describes a simple but effective technique for reducing dispersion errors in finite element solutions of timeharmonic wave propagation problems. The method involves a simple shift of the integration points to locations away from conventional Gauss or Gauss–Lobatto integration points. For bilinear rectangular elements, such a shift results in fourth-order accuracy with respect to dispersion error (error in wavelength), as opposed to the second-order accuracy resulting from conventional integration. Numerical experiments involving distorted meshes indicate that the method has superior performance in comparison with other dispersion reducing finite elements. 2003 Elsevier B.V. All rights reserved. Keywords: Wave propagation; Numerical dispersion; Numerical integration; Helmholtz equation
1. Introduction Finite element methods, when applied to time-harmonic wave propagation problems, incur solution errors that can be classified into error in the amplitude and error in the wavelength. In large-scale problems where the domain size is much larger than the wavelength, the error in the wavelength, known as the dispersion error, tends to accumulate and produce large phase errors, resulting in completely erroneous results [1]. On the other hand, amplitude error is relatively benign in that it does not increase with growing domain size. This paper focuses on the reduction of dispersion error. The dispersion error could be reduced either by using a finer mesh, or by using higher order finite elements. Both of these options come with large increase in computational cost, making some of the largescale simulations prohibitively expensive. It is thus desirable to reduce the dispersion for lower order finite elements. It can be quickly realized that dispersion error in one-dimensional problems can be completely eliminated by resorting to analytical solutions [1]. For higher dimensions, the dispersion error cannot be completely eliminated [2], but can be reduced. Over the past two decades, several researchers have proposed
*
Corresponding author. Tel.: +1-919-515-7699; fax: +1-919-515-7908. E-mail address: [email protected] (M.N. Guddati).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.09.010
276
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
various methods for reducing the dispersion in higher dimensions, especially for low-order finite elements. The following is a brief survey. The peculiar aspect of dispersion error in higher dimensions is its anisotropy, i.e., the dispersion error depends on the angle between the wave propagation direction and mesh orientation. Mullen and Belytchko [3] have analyzed this property, and compared the performance of various finite element as well as standard finite difference techniques. Marfurt [4] reduced the dispersion error by taking generalized average of consistent and lumped mass matrices, but the anisotropy still remains. Later, Krenk [5] proposed selective scaling of stiffness and mass matrices to reduce the dispersion error. The implication of his method for distorted meshes is not clear. More recently, Harari and Hughes [6] have used Galerkin/least squares (GLS) method [7] to reduce dispersion on square meshes. The resulting dispersion error, although reduced, is still significant. Thompson and Pinsky [8] modified this method to further reduce the dispersion, but their method still suffers from anisotropy. Babuska and Ihlenburg [9] used the concepts of generalized finite element method (GFEM) and finite difference stencils to derive a method that reduces the dispersion by several orders of magnitude, but only for square meshes. Oberai and Pinsky [10] used the concept of edge residuals to derive the same scheme in the context of finite elements. However, when extended to rectangular or unstructured meshes, their scheme becomes heuristic in nature and loses its accuracy properties. Babuska and Melenk [11] used a partition of unity method to reduce the dispersion, but for the method to be practical, one needs to know the behavior of the solution a priori, which is often not the case. The method of residual free bubbles, developed by Franca et al. [12], is a robust method that can improve the accuracy by effective enrichment of shape functions, but such enrichment demands significant increase in computational cost, especially for non-uniform meshes. Similar comment applies to the sub-grid modeling approach proposed by Cipolla [13]. Weighted-averaging FEM proposed by Min et al. [14] significantly reduces dispersion, but only on square meshes and at significant increase in computational cost. A fundamental limitation of most of the above methods is that they are applicable only for propagating waves in homogeneous media discretized by regular square meshes. Many important and practical problems of wave propagation involve complex geometries and material heterogeneities, and cannot be discretized using square meshes. Although dispersion reducing methods for triangular meshes may offer added flexibility [15], the current versions are limited to uniform structured meshes. With the aim of reducing dispersion on unstructured meshes, we recently proposed a dispersion-reducing framework named the Local Mesh-Dependent Augmented Galerkin (L-MAG) methods [16]. L-MAG methods provide a systematic procedure to reduce the dispersion error on rectangular meshes. The methods are more flexible than the other dispersion reducing methods in that they are accurate for meshes containing rectangular elements of varying size. However, although they perform better than other dispersion reducing methods, L-MAG methods lose their accuracy properties when used for distorted meshes. In this paper, we develop an effective dispersion reducing technique that is significantly simpler than all the existing dispersion reducing methods. The method is based on the simple observation that by shifting the Gauss quadrature to some unconventional locations would reduce the anisotropic dispersion significantly. In fact, for rectangular elements, the proposed method makes the numerical wavelength fourth-order accurate, as opposed to the second-order accuracy obtained from conventional integration rules. Needless to say, the proposed method can be directly applied to distorted meshes. The remainder of the paper focuses on the derivation of the method, its application to structured and unstructured meshes and the comparison of its performance with some existing dispersion reducing finite element methods. The following is the outline of the paper. Section 2 contains discussion of the variational boundary value problem and the Galerkin finite element discretization of the problem. The idea of shifting the integration points is discussed in Section 3. The locations of the integration points that reduce the dispersion error are obtained in Section 4. Section 5 focuses on evaluating the performance of proposed method with the help of numerical examples. The paper is concluded in Section 6 with some closing remarks.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
277
2. Finite element formulation of time-harmonic wave propagation problems Helmholtz equation, also called the reduced wave equation, is used to analyze scalar wave propagation encountered in linear acoustic and anti-plane shear problems. The corresponding boundary value problem takes the form: find u : X ! C, such that r2 u k 2 u ¼ f
ð1Þ
in X
with boundary conditions u¼g
ð2Þ
on Cg ;
ru n þ bðuÞ ¼ F
ð3Þ
on Ch :
In the above, u is the complex-valued field variable (acoustic pressure, or anti-plane displacement), k is the wave number associated with the excitation frequency (k ¼ x=c), g is the specified field variable on the Dirichlet boundary Cg , whereas a Robin type boundary condition is specified on the remainder of the boundary Ch . Here, b is an operator ranging from zero (for pure Neumann boundary condition) to the DtN map in the context of absorbing boundaries. Using classical variational calculus, the above differential (strong) form of the boundary value problem can be converted to the following variational (weak) form: find u : X ! C, with u ¼ g on Cg , such that ðrv; ruÞX k 2 ðv; uÞX þ ðv; bðuÞÞCh ¼ ðv; f ÞX þ ðv; F ÞCh
8v; v ¼ 0 on Cg
The notation ð; ÞX is used to define inner product, over the domain X, as follows: Z ðf ; gÞX ¼ f g dX;
ð4Þ
ð5Þ
X
where the f is the complex conjugate of f . The notation ð; ÞCh represents inner product taken over the boundary Ch . For numerical solution of the above variational boundary value problem, (Bubnov) Galerkin FEM can be used, where the Galerkin approximation of the variational boundary value problem takes the form: find uh : X ! C, with uh ¼ gh on Cg , such that ðrvh ; ruh ÞX k 2 ðvh ; uh ÞX þ ðvh ; bðuh ÞÞCh ¼ ðvh ; f ÞX þ ðvh ; F ÞCh
8vh ; vh ¼ 0 on Cg :
ð6Þ
In the above, uh and vh are approximations of u and v, respectively, and gh is the consistent approximation of the Dirichlet data. In the context of conventional finite element method, uh and vh take the form uh ðxÞ ¼ NðxÞU and
vh ðxÞ ¼ NðxÞV;
ð7Þ
where NðxÞ is the shape function (row) vector, while U is the discretized field variable (degree of freedom vector) and V is its variation. By substituting (7) in (6), we obtain the algebraic equation SU ¼ F;
ð8Þ
where S is the coefficient matrix given by S ¼ K k 2 M þ ðN; BðNÞÞCh
ð9Þ
and F is the force vector F ¼ ðN; f ÞX þ ðN; F ÞCh :
ð10Þ
In (9), note that B is the discretized version of the operator b in (3) and K ¼ ðrN; rNÞX ;
M ¼ ðN; NÞX
ð11Þ
278
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
are stiffness and mass matrices, respectively. These system matrices are obtained by assembling the element matrices that are given by Z Z T ðrNÞ rN dX; Me ¼ NT N dX; ð12Þ Ke ¼ Xe
Xe
where Xe is the domain associated with the element.
3. Generalized integration rules for evaluation of stiffness and mass matrices Since finite elements can be distorted, it is customary to map the finite elements from a square parent element (f1 < n < þ1g f1 < g < þ1g), and transform the above integrals into the natural coordinate system as follows: Z þ1 Z þ1 Z þ1 Z þ1 T Ke ¼ ðrNÞ rN detðJÞ dn dg; Me ¼ NT N detðJÞ dn dg: ð13Þ 1
1
1
1
The above integrals are traditionally evaluated using Gauss quadrature Z þ1 Z þ1 ng X / dn dg ¼ Wi /ðni ; gi Þ; 1
1
ð14Þ
i¼1
where ng is the number of Gauss points, Wi are the weights and ðni ; gi Þ are the local coordinates of Gauss pffiffiffi points. For pffiffiffi the most commonly used bilinear quadrilateral element, a 2 · 2 Gauss rule (ni ¼ 1= 3; gi ¼ 1= 3; Wi ¼ 1) is used to evaluate the stiffness matrix. Such an integration rule is considered desirable, as it represents exact integration for square and rectangular elements. On the other hand, the mass matrix is pffiffiffi pffiffiffi evaluated either by 2 · 2 Gauss rule (ni ¼ 1= 3; gi ¼ 1= 3; Wi ¼ 1) resulting in the consistent mass matrix, or Gauss–Lobatto rule (ni ¼ 1; gi ¼ 1; Wi ¼ 1) resulting in lumped mass matrix. Lumped mass matrix is preferred for explicit computations in time-domain, while there is no special preference for frequency domain solutions. Unlike consistent integration, lumped mass evaluation incurs errors in numerical integration. However, since the integration errors are of the same (second) order as the solution error, lumped mass is considered an acceptable alternative to consistent mass matrices. These observations indicate that any integration rule that has second-order accuracy would be a viable alternative to the conventional integration schemes. Here we introduce such a generalized 2 · 2 integration rule as given below Ke ¼
4 X
T
Wi K ðrNðnKi ; gKi ÞÞ rNðnKi ; gKi Þ detðJÞ
ð15Þ
i¼1
and Me ¼
4 X
M M M Wi M NðnM i ; gi ÞNðni ; gi Þ detðJÞ:
ð16Þ
i¼1
Considering the symmetry of the element, we propose to use symmetric, yet general, locations of integration M points with nKi ¼ a and gKi ¼ a for the stiffness terms, and nM i ¼ b and gi ¼ b for the mass terms, K M with the weights Wi ¼ Wi ¼ 1. As will be shown in the next section, the above generalization gives us the flexibility to reduce the dispersion error through appropriate choice of the location of integration points.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
279
4. Dispersion reducing integration rule In this section, by taking advantage of the flexibility offered by the generalized integration rule, we try to fix the integration points so that the dispersion error is minimized for rectangular meshes. The coefficient matrix (S ¼ K k 2 M) of a four-node rectangular finite element of size Dx Dy, when computed using Eqs. (15) and (16), is given by 2 3 S0 Sx Sxy Sy 6 Sx S0 Sy Sxy 7 7 ð17Þ Se ¼ 6 4 Sxy Sy S0 Sx 5; Sy Sxy Sx S0 where ðDx2 þ Dy 2 Þða2 þ 1Þ DxDy 2 ð1 þ b2 Þ k 2 ; 4DxDy 16 ðDx2 Dy 2 Þ a2 ðDx2 þ Dy 2 Þ DxDy ð1 b4 Þk 2 ; Sx ¼ 4DxDy 16 ðDy 2 Dx2 Þ a2 ðDx2 þ Dy 2 Þ DxDy ð1 b4 Þk 2 ; Sy ¼ 4DxDy 16 ðDx2 þ Dy 2 Þða2 1Þ DxDy 2 ð1 b2 Þ k 2 : Sxy ¼ 4DxDy 16
S0 ¼
ð18Þ
In order to obtain the equations associated with the central node (x; y) of the finite element patch in Fig. 1, we need to assemble all four elements connected to the node. Such an assembly will result in the following form of the equation: 4S0 uhx;y þ 2Sx ðuhxDx;y þ uhxþDx;y Þ þ 2Sy ðuhx;yDy þ uhx;yþDy Þ þ Sxy ðuhxDx;yDy þ uhxþDx;yDy þ uhxþDx;yþDy þ uhxDx;yþDy Þ ¼ 0:
ð19Þ
This equation can also be obtained by an equivalent finite difference method using a nine-node stencil shown in the Fig. 1. Noting that (19) is applicable for all nodes of the discretized system, it can be shown that the solution for the difference equations takes the form of a plane wave uh ¼ eik
h ðx cos hþy cos hÞ
ð20Þ
:
4
∆x
3 ∆y
(x
i +1
(x
, y j −1 )
i −1
j +1
i
, yj )
i −1
, y j− 1 )
(x , y ) i
j −1
(x
i +1
(x , y ) i
(x
1
(x , y )
, y j +1 )
(x
, yj )
(x
, y j−1)
i +1
2
j
i +1
Fig. 1. Nine-node four-element rectangular finite element patch for dispersion analysis (nine-node stencil in case of finite differences).
280
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
The only difference between the approximate solution and the exact solution (u ¼ eikðx cos hþy sin hÞ ) is that the wave numbers, k h and k, are different. In fact, the relative dispersion error can be quantified as the relative error in wavelengths k kh k h k ¼ ð21Þ erelative ¼ dispersion k k h : In the above equations, the approximate wave number k h can be obtained by substituting the approximate waveform in (20) into the difference equation (19). Such substitution, after some simplification, results in ½S0 þ Sx Cx þ Sy Cy þ Sxy Cx Cy uhx;y ¼ 0;
ð22Þ
where Cx ¼ cosðk h Dx cos hÞ;
Cy ¼ cosðk h Dy sin hÞ:
ð23Þ
The expression in the brackets in (22) must be zero for the non-trivial solution (ux;y 6¼ 0). Thus the discrete dispersion relationship relating the exact and approximate wave numbers is given by S0 þ Sx Cx þ Sy Cy þ Sxy Cx Cy ¼ 0: Using (18), the dispersion relation can be written as 4 ðDx2 þ Dy 2 Þða2 þ Cx Cy a2 Cx a2 Cy a2 þ 1 Cx Cy Þ þ ðDx2 Dy 2 ÞðCx Cy Þ h i Dx2 Dy 2 k 2 ð1 þ b2 Þ2 þ Cx Cy ð1 b2 Þ2 þ ðCx þ Cy Þð1 b4 Þ ¼ 0
ð24Þ
ð25Þ
and the relation between the exact wave number k and the approximated wave number k h is obtained as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðDx2 þ Dy 2 Þða2 þ Cx Cy a2 Cx a2 Cy a2 þ 1 Cx Cy Þ þ ðDx2 Dy 2 ÞðCx Cy Þ k¼ : ð26Þ 2 2 DxDy ð1 þ b2 Þ þ Cx Cy ð1 b2 Þ þ ðCx þ Cy Þð1 b4 Þ In order to obtain the relative dispersion error, we take Taylor expansion of k in terms of k h . After some simplifications, the dispersion error can be written as k kh k2 2 2 2 2 2 2 2 2 2 2 k h ¼ 24 ð3a 2ÞðDx þ Dy Þ cos h sin h þ ð2 3b ÞðDx cos h þ Dy sin hÞ þ Oðk 4 ðDx þ DyÞ4 Þ:
ð27Þ
Using the above equation, we can obtain the estimates pffiffiffiof the dispersion errors for classical integration rules. When consistent integration is used (a ¼ b ¼ 1= 3), the dispersion error becomes k kh k2 2 4 4 4 2 4 ð28Þ k h ¼ 24 Dx cos h þ Dy sin h þ O k ðDx þ DyÞ : pffiffiffi When lumped mass is used (a ¼ 1= 3; b ¼ 1) we obtain an analogous result k kh k2 ¼ ðDx2 þ Dy 2 Þ cos2 h sin2 h þ ðDx2 cos2 h þ Dy 2 sin2 hÞ þ O k 4 ðDx þ DyÞ4 : ð29Þ k h 24 From the above two estimates, it is clear that conventional integration rules result in second-order accuracy of the numerical wave number. theffi other hand, On pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi ffi we observe that a simple shift of the integration points to unconventional location of 2=3; 2=3 , i.e.,
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
2 a2 ¼ b2 ¼ ; 3
281
ð30Þ
eliminates the second-order term in the error (27). It can be shown that, for this case, the dispersion error is given by k kh k 4 4 6 6 6 4 6 ð31Þ k h ¼ 480 Dx cos h þ Dy sin h þ Oðk ðDx þ DyÞ Þ: Note that the resulting dispersion error is similar to that of L-MAG2 [16], but the proposed method is much simpler. pffiffiffiffiffiffiffiffi The choice of b ¼ 2=3 for mass integration is not completely unconventional. For rectangular elements, it can be shown to be equivalent to averaging lumped and consistent mass matrices (see e.g. [4]). pffiffiffiffiffiffiffi ffi However, the choice of a ¼ 2=3 is unconventional and is the key to reducing the anisotropy in numerical dispersion. Incorporating the proposed modified integration rule into finite element software is trivial and imposes no additional computational cost. The proposed modification is completely local (element-level). Furthermore, unlike other dispersion reducing methods, the implementation of this method is independent of the element distortion. These observations indicate that the method is not only effective for uniform rectangular meshes, but could also work well for non-uniform distorted meshes. Further it is also straightforward to apply the proposed method to three-dimensional problems.
5. Numerical examples The performance of the proposed modified integration rule is tested using three examples: simulation of plane wave, simulation of radiation from a point source, and simulation of asymmetric radiation from a circular cavity. All the problems are solved using the proposed method, as well as some existing methods, namely Galerkin method with consistent integration (referred to as Galerkin in the figures), modified Galerkin least squares method [8] (GLSm), generalized finite element [9] or residual-based method [10] (Residual), and the second- and sixth-order local mesh-dependant Galerkin (L-MAG) methods proposed recently by the authors [16] (L-MAG2 and L-MAG6). 5.1. A plane wave problem We consider a plane wave u ¼ eikðx cos hþy sin hÞ with k ¼ 75 propagating in a square domain. Two types of meshes (uniform and distorted rectangular elements, Fig. 2) are used to discretize the domain. Robin boundary conditions are used for the analyses. Since the dispersion error accumulates whereas amplitude error does not, several simulations with increasing domain sizes are performed to differentiate the two types of errors. The associated errors are measured in terms of the relative L2 error given by e¼
ku uh kL2 : kukL2
ð32Þ
5.1.1. Uniform rectangular mesh Five domain sizes, 0.2 · 0.2, 0.5 · 0.5, 1 · 1, 1.5 · 1.5 and 2 · 2 are considered. The domains are discretized using 0.01 · 0.005 rectangular elements implying that kDx ¼ 0:75 and kDy ¼ 0:375 for each case. The simulation is performed for nine different propagation angles ranging from 0 to 90. For the domain size of 1 · 1, Fig. 3 shows comparison of errors resulting from various methods. Errors for other domain sizes have
282
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
Fig. 2. Meshes used in numerical examples: (a) uniform rectangle element and (b) distorted rectangular element.
10 Galerkin
GLSm
Residual
L-MAG2
L-MAG6
Proposed
Relative error
1
0.1
0.01
0.001 0
11.25
22.5
33.75
45
56.25
67.5
78.75
90
Propagation direction (degrees)
Fig. 3. Error for varying angles of propagation (with uniform rectangular elements).
2 Galerkin GLSm Relative error
1.5
Residual L-MAG2 L-MAG6
1
Proposed
0.5
0 0.2
0.5
0.8
1.1 Domain size
1.4
1.7
2
Fig. 4. Error for varying domain sizes (with uniform rectangular elements).
similar variation. It can be clearly seen that L-MAG and proposed methods have the best performance, with their errors much less than the errors from Galerkin, GLSm and Residual/GFEM methods. In order to illustrate the effect of the reduction in dispersion error, the error for various domain sizes is plotted with propagation angle fixed at 0. Fig. 4 shows the increasing error as the domain size is increased,
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
283
which is essentially due to the accumulative nature of the dispersion error. It can be clearly seen that the error increases significantly for Galerkin, GLSm and residual/GFEM methods, while the error growth is minimal for L-MAG methods as well as the proposed method. This clearly illustrates that the proposed method, as well as L-MAG methods, can be used to solve large-scale wave propagation problems with rectangular meshes. 5.1.2. Distorted rectangular mesh The distorted mesh is made up of rectangular patches consisting of four distorted elements as shown in Fig. 5. Except for the element distortion, all the parameters are identical to the case of uniform rectangular meshes. The errors from various methods with different propagation directions for the unit domain are shown in Fig. 6. Errors for other domain sizes have similar variation. The increasing error as the domain size is increased is plotted in Fig. 7. It is clear that the L-MAG methods have smaller errors than other existing methods. However, the striking feature is that the proposed method even outperforms the L-MAG methods, indicating that the proposed method is best suited for meshes when there is distortion. 5.2. Radiation from point source The time harmonic radiation from a point source in infinite domain is considered. The computational domain is a unit square with unit pulse at the lower left corner. The wave number is chosen as 75. The
Fig. 5. Layout of the distorted mesh (h ¼ 0:005).
10 Galerkin
GLSm
Residual
L-MAG2
L-MAG6
Proposed
Relative error
1
0.1
0.01
0.001 0
11.25
22.5 33.75 45 56.25 67.5 Propagation direction (degrees)
78.75
90
Fig. 6. Error for varying angles of propagation (with distorted rectangular elements).
284
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287 2 Galerkin GLSm Relative error
1.5
Residual L-MAG2
1
L-MAG6 Proposed
0.5
0 0.2
0.5
0.8
1.1 Domain size
1.4
1.7
2
Fig. 7. Error for varying domain sizes (with distorted rectangular elements).
Real part of displacement
0.2 0.15
Galerkin
GLSm
Residual
L-MAG6
Proposed
Exact
L-MAG2
0.1 0.05 0 0.90
0.92
0.94
0.96
0.98
1.00
-0.05 -0.1 Distance from center (r)
Fig. 8. Radiation from point source (with uniform rectangular elements): response along the radius.
meshes used in the previous example are used for this example, i.e. elements have aspect ratio of 2, with khx ¼ 0:75, khy ¼ 0:375 (see Fig. 2). A PML type absorbing boundary condition is used to simulate wave radiation into the exterior [17]. The results from various methods are compared with the exact solution, which is similar to the GreenÕs function of the two-dimensional HelmholtzÕs operator [18], and is given by ð1Þ
uðrÞ ¼ iH0 ðkrÞ;
ð33Þ
ð1Þ
where H0 is the zeroth-order Hankel function of the first kind. The real parts of the response along h ¼ 0 with r ranging from 0.9 to 1.0 are shown in Figs. 8 and 9 for uniform and distorted rectangular meshes, respectively. It can be seen that, consistent to the previous observations, L-MAG methods and proposed method have superior accuracy properties for rectangular meshes, while the proposed method outperforms all other methods for distorted meshes. 5.3. Radiation from an infinitely long circular cylinder This example considers the radiation from an infinitely long, circular cylinder (R ¼ 1) into an unbounded fluid domain. It reduces to a two-dimensional problem when the prescribed pressure on the wet surface does not vary along the axis of the cylinder. The wave number for this problem is chosen to be 4p. The Dirichlet boundary condition on the surface of the cylinder (r ¼ R) is chosen as u ¼ cos h.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
285
Fig. 9. Radiation from point source (with distorted rectangular elements): response along the radius.
5
4
3
2
1
0
0
1
2
3
4
5
Fig. 10. Radiation from an infinitely long circular cylinder: finite element mesh.
The unbounded domain is truncated using a 10 · 10 square around the cavity, with PML absorbers applied at the computational boundary. Furthermore, due to symmetry, only one-quarter of the problem is analyzed. The mesh is shown in Fig. 10. Sixty elements are used along the radial direction and 100 in circumferential direction. The mesh consists of elements with aspect ratio as high as 6 (see inset A), kh as high as 1.9 (see inset B), indicating only 3–4 elements per wavelength. Furthermore, most elements are distorted. Various methods are used to analyze this problem and the obtained solution is compared with the exact solution given by ð1Þ
uðr; hÞ ¼ ð1Þ
H1 ðkrÞ ð1Þ H1 ðkRÞ
cos h;
ð34Þ
where H1 is the first-order cylindrical Hankel function of the first kind. The real parts of results along h ¼ 0 with r ranging from 4.5 to 5.0 are shown in Fig. 11. It is observed that the dispersion errors from LMAG and proposed methods are much smaller than other methods, while the results from L-MAG6 and the proposed method are the closest to the exact solution.
286
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287 0.8
Real part of displacement
0.6
Galerkin
GLSm
Residual
L-MAG6
Proposed
Exact
L-MAG2
0.4 0.2 0 4.5
4.6
4.7
4.8
4.9
5.0
-0.2 -0.4 -0.6 Distance from center (r)
Fig. 11. Radiation from source point: response along the radius.
6. Concluding remarks In this paper, we developed a simple but effective method for reducing numerical dispersion in finite element solutions of time-harmonic wave pffiffiffiffiffiffiffipropagation ffi pffiffiffiffiffiffiffiffi problems. By simply shifting the integration points to the unconventional location of ð 2=3; 2=3Þ, we were able to obtain fourth-order accuracy in the wave number, as opposed to the second-order accuracy resulting from conventional methods. In spite of its simplicity, the proposed method appears to outperform more involved dispersion reducing methods, especially for rectangular and distorted elements. Needless to say, the implementation of the proposed method into existing finite element programs is a trivial issue and there is no additional computational cost associated with the method. It should be noted that when only square elements are used, the GFEM/residual methods, as well as L-MAG6, would have the best performance, resulting in sixth-order accuracy in the wavelength. When the mesh contains distorted and rectangular elements, GFEM/residual methods quickly lose their accuracy. L-MAG6 method maintains its accuracy for rectangular meshes, but deteriorates when applied to distorted meshes. On the other hand, the proposed method retains accuracy even for distorted meshes, while being almost as accurate as L-MAG6 for rectangular meshes. Based on these observations, we advocate the use of the proposed modified integration rule whenever the elements are (even slightly) distorted, L-MAG 6 when rectangular, L-MAG6/GFEM/residual methods, when the mesh contains only square elements. The proposed method, although developed for two-dimensional problems, is directly extendible to threedimensional problems. It may be possible to further increase the accuracy of the proposed method by combining the proposed idea with the ideas of L-MAG methods. It may also be possible to devise modified integration rules for higher order finite elements, as well as for finite elements for more complex waves such as elastic and flexural waves. These extensions are subjects of ongoing investigations.
Acknowledgements This material is based upon work supported by the National Science Foundation under Grant No. 0100188. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.
M.N. Guddati, B. Yue / Comput. Methods Appl. Mech. Engrg. 193 (2004) 275–287
287
References [1] F. Ihlenburg, in: Finite Element Analysis of Acoustic Scattering, Applied Mathematical Sciences, vol. 132, Springer-Verlag, New York, 1998. [2] I.M. Babuska, S.A. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, S.I.A.M. J. Numer. Anal. 34 (1997) 2392–2423. [3] R. Mullen, T. Belytschko, Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation, Int. J. Numer. Methods Engrg. 18 (1982) 11–29. [4] K.J. Marfurt, Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation, Geophysics 49 (1984) 533–549. [5] S. Krenk, Optimal formulation of simple finite elements, in: Variational Methods in Engineering, Proceedings of the 2nd International Conference, Springer, Berlin, 1985. [6] I. Harari, T.J.R. Hughes, Galerkin/least squares finite element method for the reduced wave equation with non-reflecting boundary conditions, Comput. Methods Appl. Mech. Engrg. 98 (1992) 411–454. [7] T.J.R. Hughes, Multiscale phenomena: GreenÕs functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles, and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. [8] L.L. Thompson, P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38 (1995) 371–397. [9] I.M. Babuska, F. Ihlenburg, A generalized finite element method for solving the Helmholtz equation in two-dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128 (1995) 325–359. [10] A.A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 49 (2000) 399–419. [11] I.M. Babuska, J.M. Melenk, The partition of unity finite element method, Int. J. Numer. Methods Engrg. 40 (1997) 727–758. [12] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residue-free bubbles for the Helmholtz equation, Int. J. Numer. Methods Engrg. 40 (1997) 4003–4009. [13] J.L. Cipolla, Subgrid modeling in a Galerkin method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 177 (1999) 35–49. [14] D.J. Min, H.S. Yoo, C. Shin, H.J. Hyun, J.H. Suh, Weighted-averaging finite-element method for scalar wave equation in the frequency domain, J. Seism. Explor. 11 (2002) 197–222. [15] I. Harari, C.L. Nogueira, Reducing dispersion of linear triangular elements for the Helmholtz equation, J. Engrg. Mech. ASCE 128 (2002) 351–358. [16] M.N. Guddati, B. Yue, Local mesh-dependent augmented Galerkin methods for reducing dispersion error for Helmholtz equation, Int. J. Numer. Methods Engrg. (submitted). [17] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Computat. Phys. 114 (1994) 185–200. [18] P.M. Morse, K.U. Ingard, Theoretical Acoustics, Princeton University Press, Princeton, NJ, 1986.
