- Email: [email protected]
A higher-order velocity post-processing for the Helmholtz equation
Finite Elements in Analysis and Design 41 (2005) 459 – 492 www.elsevier.com/locate/finel
A higher-order velocity post-processing for the Helmholtz equation Rivânia H. Paulinoa , Fernando Alves Rochinhab,∗ a Department of Mechanical Engineering, Federal University of Espirito Santo, UFES, Vitória, ES, Brazil b Department of Mechanical Engineering, Federal University of Rio de Janeiro, UFRJ, Ilha do Fundao, Cx. Postal 68503, Rio
de Janeiro, RJ 21945-970, Brazil Received 29 October 2003; accepted 20 May 2004 Available online 25 December 2004
Abstract A macroelement post-processing recovery technique of the velocity field based on residuals of the equilibrium equation and irrotational conditions for the Helmholtz equation is proposed. The derivatives are recovered by solving local variational problems, that are based on the residuals aforementioned and involve superconvergent derivatives of the acoustical pressure field on special points. Improved accuracy for derivative finite element approximation is obtained with low computational cost and easy implementation. A number of comprehensive simulations is presented in order to show the main features of the proposed methodology. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Acoustics; Post-processing; Flux recovery
1. Introduction In acoustic problems, the steady-state sound response in a non-viscous medium is governed by the Helmholtz equation [1]. One major problem in numerically solving this equation by means of the finite element method (FEM) is linked to a potential loss of ellipticity with an increasing excitation frequency. This degradation of the standard Galerkin FEM approximation is specially observable in the convergence behavior of the acoustic velocity, a field directly related to the gradient of the primal variable. ∗ Corresponding author. Tel.: +55 21 25628315.
E-mail address: [email protected] (F.A. Rochinha). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.05.005
460
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Classical Galerkin FEM provides good phase and amplitude accuracy as long as the mesh is fine enough with respect to the wave number in the propagation region. However, “fine enough” is often too expensive for adequate resolution, even for moderate wave numbers. Hence, one of the main concerns in acoustic finite element analysis is the adequacy of the finite element mesh. Acousticians often use the, so-called, “rule of the thumb” [2] which prescribes a relation between the minimal number of elements and the wave number. This rule only takes into account the capacity of the interpolation given by a specific mesh and is not based on stability aspects. Therefore, as it is well documented in the literature [1–4], one can obtain very bad results even when the finite element mesh obeys this rule. Those approximation problems have been addressed in a number of publications leading to several solutions like the ones contained in [1,5–9]. The numerical phenomena was analyzed in [2] where the following error estimate for the pressure gradient is presented, when linear elements are used. eh C1 kh + C2 k 3 h2 ,
kh < 1
(1)
in which constants C1 and C2 are independent of wave number k and element size h. The first term on the right-hand side of Eq. (1) represents the interpolation error and the second one represents the pollution error. In order to keep the local error constant, maintaining kh constant is enough. That is the “rule of the thumb” which corresponds to taking a certain number of elements per wavelength. As one can see from Eq. (1), this is not sufficient to maintain the pollution error under control. Thus, even following the heuristic rule, the pollution error increases with k. In some applications, ranging from very specific situations like the absorption analysis in silencers [10] to more general scenarios involving inverse problems [11,12] or error estimators [1,5,13,14], the velocity field which is directly linked to the pressure gradient is of main interest. As previously mentioned, the differentiation of pressure computed by standard finite element approximation leads to poor results. Thus, alternative formulations have been developed aimed at improving the velocity field. Here, a macroelement post-processing approach [15,16] is developed to recover the gradient of the pressure with improved accuracy. The macroelement post-processing is constructed over subregions of the finite element mesh and is based on the superconvergence of the derivatives of the pressure finite element solution at special points combined with a residual form of a balance equation and the irrotationality condition. Numerical analysis and error estimates are presented which demonstrate that higher asymptotic convergence rates are obtained with the proposed approach. An outline of the paper follows. A brief review of the basic problem leading to the Helmholtz equation is presented in Section 2. The macroelement post-processing is introduced and analyzed in Section 3. Numerical experiments confirming the predicted rates of convergence are reported in Section 4. Finally, Section 5 contains some concluding remarks. 2. Problem statement and finite element formulation 2.1. Continuous problem The expression acoustic problem used here refers to the sound propagation in a non-viscous fluid in equilibrium. All non-linear effects are neglected and only a time-harmonic response is considered leading to the Helmholtz equation, in which the fluctuation of the sound pressure is seek. Although it is a very
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
461
well-known modelling, the basic equations will be briefly reviewed below in order to help the developing of the recovery strategy which will be proposed later on. The modelling relies on three equations, namely: continuity, Euler and state. They involve the velocity field V=V(x, t), the pressure P =P (x, t) and the density ¯ = ¯ (x, t). The first two equations are presented right below in the same order they were mentioned in the beginning of this paragraph j¯
+ div(¯V) = 0,
jt jV
1 + (V div(∇V)) + ∇P = 0, jt ¯
(2) (3)
where div( ) = j( )/jx + j( )/jy is the divergence operator and ∇( ) = (j( )/jx)ex + (j( )/jy)ey is the gradient with ex and ey as the unitary vectors x and y directions. Assuming that the sound propagation will produce only slight perturbations on the flow, the following fluctuations of the basic fields are introduced: = ¯ − 0 , v = V and p = P − p0 , where the subindex 0 refers to the uniform state. With this in mind, the non-linear terms contained in the above equations are linearized giving rise to j jt 0
+ 0 div v = 0, jv jt
+ ∇p = 0.
(4) (5)
The above linearized equations are complemented by the linear state equation p = c2 , where c stands for the sound velocity in the fluid. As it was mentioned earlier, only harmonic solutions will be treated here. Therefore, all fields can be recast as f (x, t) = fˆ(x)eiwt , where fˆ(x) is a stationary function and w is the circular frequency. Thus, Eq. (5) assumes the following form: vˆ = −
1 ∇ p. ˆ iw0
(6)
Eq. (6) will be referred to as the constitutive equation due to its resemblance to the constitutive relations in structural dynamics. The combination of Eq. (4) and the state equation leads to iw pˆ + div vˆ = 0. c 2 0
(7)
Moreover, directly from (6) one achieves the irrotationality condition rot(ˆv) = 0,
(8)
where rot stands for the rotational operator. The last three equations form the basis for the proposed velocity recovery method. For the remainder of this analysis, the symbol ˆ. will be omitted because all the variables appearing till the end are the stationary value of the fields. Therefore, there is no need to distinguish in the adopted notation.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
462
Finally, the Helmholtz equation can be retrieved by combining Eqs. (6) and (5) leading to the following mathematical problem Find p : → C such that ∇ 2 p + k 2 p = 0,
(9)
where wave number k is defined by k=
(10)
c
and ∇ 2 stands for the Laplacian operator and C is the set of complex numbers. In order to have a complete description of the wave propagation phenomena, it is necessary to identify the appropriate boundary conditions, which, in general terms, are set below. Only internal acoustic will be addressed here, so boundary conditions associated to infinite domains like the absorbing ones will not be included. • Prescribed pressure on D (Dirichlet boundary condition) p = p˜
on D .
(11)
• Prescribed acoustic velocity on N (Neumann boundary condition) vn = v˜ n
on N .
(12)
• Fixed impedance in the boundary R (Robin boundary condition) ∇p · n = −ikc0 An p
on R
(13)
where An is the admittance coefficient, the inverse of impedance in the boundary. 2.2. Variational formulation and finite element discretization In order to solve the interior acoustic problem introduced above by means of the FEM, the following variational formulation, based on the primal variable p, is stated. Problem P. Find p ∈ V such that 2 ∇p.∇q d − k pq d − ik 0 c ∇p · nq d = 0 ∀q ∈ V , (14)
N ∪ R
where . stands for the ordinary scalar product in the complex plane and the function space of admissible pressures V is defined as V = {p ∈ H1 () | q = 0 on D }, where H1 () denotes the Sobolev space of complex-valued functions 2 jf 1 2 H () = f : → C, |f (x)| dx < ∞, (x) dx < ∞ . jx H1 () is a subset of L2 () the space of square integrable complex functions.
(15)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
463
The Galerkin finite element approximation of Eq. (14) is obtained by considering a finite dimensional space Vh ⊂ V , and then seeking an approximate solution to (14). Let Shs be a quadrilateral or triangular C0 Lagrangian finite element subspace of H1 () of degree s. The discrete version of Problem P in Vh = Shs ∩ V reads Problem Ph : Find ph ∈ Vh , such that 2 ∇ph .∇qh d − k ph qh d − ik 0 c ∇ph · nqh d = 0 ∀qh ∈ Vh , (16)
N ∪R
where h represents the mesh parameter. The above discrete formulation leads to the matrix equation (K + ik cC − k 2 M)q = −ik cb,
(17)
where K, M are often referred to as rigidity and mass matrices and q is the vector of Np nodal values of the discrete pressure ph . Besides, C and b are associated to boundary terms. This terminology is motivated by the similarity of (17) to the structural dynamic equations. The components of these matrices and vectors are given by Kij = ∇ j · ∇ i d, Mij = j i d, Cij =
R
(∇ j · n)i d,
bi =
N
(v · n)i d
with j (j = 1, . . . , Np ) denoting the standard finite element basis. The convergence analysis of the above finite element problem is addressed in [2] where the emphasis is placed on the pollution effect. That numerical misbehavior reflects the lack of ellipticity of the finite element formulation in the presence of high wave numbers. In many finite element codes an heuristic strategy of choosing the mesh’s size as a function of k, often referred to as the “rule of the thumb”, is adopted in order to keep track of the high oscillatory behavior of p. The following estimate, presented in [2], for a one-dimensional prototype situation, demonstrates that even with the adoption of this rule the pollution effect strongly disturbs the quality of the solution |ph − p|1 C1 kh + C2 k 3 h2 , |p|1
kh < 1 (“ rule of the thumb” ),
(18)
where |.|1 is the H1 semi-norm and C1 and C2 are positive numbers independent of the mesh. The second term of the above bound dominates the error in the pre-asymptotic range of the finite element approximation, which implies that for reasonable meshes, on an engineering perspective, the standard analysis might entail unreliable results.
3. Macroelement post-processing technique Usually, after solving the primal problem introduced in the previous section, the velocity field is obtained by using (6). This direct approach leads to lower-order accurate derivatives when compared to
464
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
the primal variable. Nevertheless, the accuracy of the derivatives is often of main interest in engineering problems. This motivates the development of the so-called recovery strategies, in which a more accurate approximation of the derivative variable is sought. Several approaches are found in literature [17,18] and it is worth mentioning that some of the most important ones have already been tested on the Helmholtz equation with relative success [1,5,19], in the context of finite element adaptivity. The macroelement post-processing technique proposed in [15] will be extended to the acoustic problem introduced earlier. The basic idea of the method is to find a better approximation for the derivative (∇ph ), and consequently for the acoustic velocity field, by solving a local variational problem on each macroelement, understood as the union of a set of neighboring elements with common edges. This problem involves least-squares residuals of the balance equation (7) and of the irrotationality condition (8) and relation (6) at special points, which are likely to present superconvergence [15,19]. The inclusion of the residuals of the two equations is twofold, it enhances the stability of the macroelement problem and includes explicitly in the formulation conditions to be obeyed by the actual velocity field. This last aspect insures that, in a least-squares sense, the approximated velocity reproduces some important features presented by its continuous counterpart. In order to establish the proposed method, one can suppose that the domain is decomposed into macroelements not necessarily disjoints and define MP rh ⊂ L2 () and MQrh ⊂ L2 (), respectively, triangular and quadrilateral lagrangian finite element spaces of piecewise polynomial of degree r on each element and class C 0 in each macroelement but discontinuous on the macroelement boundaries. Considering Zh = MP rh × MP rh × MP rh or Zh = MQrh × MQrh × MQrh , the following residual form for the problem of retrieving a better approximation of the three-dimensional acoustic velocity field is proposed Problem ME h : Given ph ∈ Vh solution of Problem Ph , find vh∗ ∈ Zh , such that iw i ∗ 2 bh (vh , uh ) = ∀uh ∈ Zh , ∇ph , uh − 1 h ph , div uh (19) w 0 0 c 2 h where the sesquilinear form bh (., .) is defined as bh (vh∗ , uh ) = (vh∗ , uh )h + 1 h2 (div vh∗ , div uh ) + 2 h2 (rot vh∗ , rot uh ), where vh∗ is the post-processed velocity. A compact notation was adopted, in which (., .) represents the L2 () product and (., .)h denotes that this term is evaluated by a suitable integration rule, taking into account only specific points where superconvergence of the derivatives are supposed to occur. Taking, for instance, macroelements composed by four bilinear quadrilaterals in a 2 × 2 arrange, as depicted in Fig. 1(a), and evaluating (·, ·)h with one Gauss point per element, one has a configuration similar to the biquadratic element. The scalars 1 and 2 are positive real parameters that might be adjusted in order to have better convergence performance. Whenever N = 0, the velocity at this part of the boundary can be strongly enforced, and Problem ME h is slightly modified, incorporating the boundary condition in the discrete space Zh . Concerning the computational implementation, the macroelement post-processing leads to a number of similar and independent linear problems, each one corresponding to a macroelement. The degrees of freedom associate to each node are the real and imaginary parts of the velocity components. Therefore, for two-dimensional problems, using a macroelement formed by four elements, which will be adopted in some numerical examples of this article, the total number of unknowns is 36.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
(a)
465
(b)
Fig. 1. Two different macroelement configurations composed of bilinear quadrilaterals. (•), nodal points; (+), superconvergent points; (a) 2 × 2 arrange; (b) 2 × 1 arrange.
This computational implementation is rather simple, requiring only a standard finite element code in order to obtain the matrix problem to be solved at each macroelement. It is important to point out that parallel computation could be explored as the proposed method leads to independent problems covering all the domain. In such case, groups formed by a number of macroelements could be assigned to different processors. A reasonable alternative to a post-processing scheme in order to obtain accurate flux approximations is provided by mixed formulations [15]. Unfortunately so far, in the authors knowledge, there are no successful mixed methods for the Helmholtz equation encompassing general boundary conditions reported in the literature. 3.1. Convergence analysis Constructing a well-posed Problem ME h depends on the number of elements in the macroelement, on the choice of the finite element space Zh and on the number of superconvergence points adopted to evaluate the constitutive Eq. (6). In the present section those topics are investigated leading to the post-processing’s convergence analysis. It is also important to mention that no regularity issue need to be addressed, though the presence in the formulation of the differential operators div and rot, as the postprocessed velocity, over each macroelement, is a finite element function of, at least, class C 0 , therefore belonging locally to H1 . In general, convergence of a numerical approximation is assessed by combining two main issues: consistency and stability. As consistency is guaranteed by the construction of the formulation, stability remains the only aspect to be addressed in order to analyze the post-processing’s convergence. A stability check can be performed by means of a simple numerical test involving an eigenvalue problem, as proposed in [15]. In the present context, this methodology was applied to several different configurations of macroelements, that are summarized in Table 1 for the situation described in the first numerical example of the next section.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
466
Table 1 Stability study of macroelement configurations Arrangement
Elements
min
max
Stability
(a) 1 × 1
Q4 Q9 T3 T6 Q4 Q9 T3 T6 Q4 Q9 T3 T6
— 0.4110 — 0.1448 0.2649 0.7516 0.3291 0.2777 0.7500 0.9009 0.8320 0.4857
— 140.01 — 213.680 26.725 139.430 40.093 178.220 26.721 137.64 35.928 164.40
No Yes No Yes Yes Yes Yes Yes Yes Yes Yes Yes
(b) 1 × 1 (c) 2 × 1 (d) 1 × 2 (e) 2 × 2 (f) 2 × 4
The numerical results which will be presented later demonstrate that the proposed post-processing method has good accuracy and robustness. It is also remarkable that when compared to the direct calculation of the velocity field, a better performance is achieved. Stability of the proposed formulation is proved with the aid of the following associated generalized eigenvalue problem j
j
bh (h , uh ) = j (h , uh )
∀ uh ∈ Zh |ME (j = 1, . . . , N )
(20)
j
where j and h are, respectively, an eigenvalue and its correspondent eigenvector and Zh |ME corresponds to the restriction of space Zh to a single macroelement. Besides, N stands for the number of degrees of freedom associated to each macroelement. The j ’s in (20), which are all real and non-negative due to the symmetry of the involved operators, can be easily computed for different macroelement configurations by using standard eigenvalue algorithms, within a finite element code. Table 1 reports the minimum and the maximum eigenvalues for several different arrangements. The post-processed velocity vh∗ can be expressed as a linear combination of the eigenvectors yielding N N ∗ ∗ ∗ i ∗ i (21) vi h , v i h , bh (vh , vh ) = bh i=1
i=1
vi∗
are the components of the flux variable in the basis formed by the eigenvectors. From the where bilinearity of the bh (., .), considering (20) and taking profit of the orthogonality among eigenvectors, one has bh (vh∗ , vh∗ ) =
N i=1
i (vi∗ )2 (i , i )
(22)
which in turn implies ∗ 2 bh (vh∗ , vh∗ ) Min j vh .
(23)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
467
If Min j , the lowest eigenvalue, is strictly positive, inequality (23) implies that bh is a positive definite operator, which entails the stability of the formulation in a local sense. The stability on the whole domain follows by summing (23) over all macroelements. The coercivity constant given by the aforementioned eigenvalue depends on the aspect ratio of the elements contained on the macroelement but not on the mesh’s parameter h. This is assessed computing the matrices related to the eigenvalue problem using integrals performed with the aid of on a reference macroelement of unity dimensions. In order to derive error estimates for the proposed methodology, the following norm, equivalent to the canonical one in Zh , is introduced |uh | = bh (uh , uh )1/2 .
(24)
The stability results stated above are now used to prove the post-processing’s convergence, leading to a priori error estimates. Those results are presented in Lemma 2, which relies on Lemma 1. Lemma 1. Let p be the solution of the Problem P, v given by (6) and ph and vh∗ be their respective approximations by solving problems Ph and ME h . Then, considering macroelements leading to stable formulations, there exists a positive number C independent of h such that v − vh∗ C(v − uh + |v − uh | +
1 |p − ph |h + khp − ph ) k
∀ uh ∈ Zh ,
(25)
where | · |h represents the semi-norm endowed by the special integration (., .)h . A slight abuse of notation was adopted as |||.||| is only proved to be a norm in Zh . Therefore the term |||v − uh |||, as v − uh ∈ / Zh , must be understood as a compact notation for bh (v − uh , v − uh ). Proof. The variational consistency referred to above, which states that the true solution satisfies the discrete problem, is expressed as iw i 2 bh (v, uh ) = ∀ uh ∈ Zh ∇p, uh − 1 h p, div uh (26) w0 0 c2 h which along with the stability estimate (23) leads to vh∗ − uh 2 Cbh (vh∗ − uh , vh∗ − uh ) C(bh (v − uh , vh∗ − uh ) + bh (vh∗ − v, vh∗ − uh )) i ∗ ∗ C bh (v − uh , vh − uh ) + ∇(ph − p), vh − uh w h iw −1 h2 (ph − p), div (vh∗ − uh ) c 2 1 C |||v − uh ||| ||vh − uh || + |ph − p|h ||vh − uh || w w + ˜ 1 h 2 ||ph − p|| ||vh − uh || , c
(27)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
468
where ˜ 1 stands for the former parameter 1 multiplied by a constant that comes from the inverse estimate used to obtain the third term of the last inequality. Hence, from (27) 1 k ∗ ˜ vh − uh C |||v − uh ||| + |p − ph |h + 1 h ||p − ph || (28) kc c which yields the desired result after using the triangular inequality. Now the convergence result is stated for the macroelement post-processing. Lemma 2. Under the same hypotheses of Lemma 1 and assuming that the continuous fields p and v are regular enough and the existence of a superconvergence behavior of the derivative expressed as |p − ph |h (C1 k r hr+ + C2 k 2r+1 h2r+ )|p|r
0 1
(29)
there exist positive numbers C1 and C2 , independent of h, such that v − vh∗ (C1 (kh)r h + C2 k(kh)2r+1 ) |p|r
kh < 1.
(30)
Proof (Error estimate). Considering the classical finite element interpolation results inf |||v − uh ||| Chr+1 |v|r+1 ,
uh ∈Zh
inf ||v − uh || Chr+1 |v|r+1 ,
uh ∈Zh
the approximation errors established for the one-dimensional case in [2] considering kh < 1 p − ph (C1 (kh)r + C2 k 2r+1 h2r )|p|r ,
(31)
the superconvergence hypothesis (29) and that for an oscillatory solution |p|r /||p|| = O(k r ), one has v − vh∗ (Ck r+1 hr+1 + C1 (k r−1 hr+ + (kh)2r+ ) + C2 ((kh)r+1 + k 2(r+1) h2r+1 ))|p|r .
(32)
Now, clustering the similar terms, the desired result is achieved. There are a few remarks to be made about the above lemma. First, it is worth calling attention to the two parcels of the final error estimate. Similarly to the analysis presented in [2], pre and asymptotic convergence behaviors are observed. The second parcel is a direct consequence of the pollution effect. The second point to be highlighted is about the error estimate for pressure field. It is only rigorously proved to the one-dimensional case, but a number of numerical simulations reported in the literature have demonstrated that a similar behavior is observed for two-dimensional and three-dimensional situations. Therefore, the obtained estimation for the post-processed velocity can be also valid to those cases. The superconvergence hypothesis used in Lemma 1 is also inspired in numerical observations not relying on any sort of theoretical developments. Anyway, this sort of behavior is not only confirmed by the numerical experiments in Section 4, but also in [5,19]. Estimate (32), in the asymptotic range, represents a higher convergence rate when compared to the velocity obtained by the direct differentiation of the finite element approximation of the pressure field.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
469
This is fairly confirmed by numerical examples presented in the next chapter. On the other hand, the pre-asymptotic range is a result of the pollution contained in ph . This error on the post-processed velocity would be controlled by a mesh, probably, not suitable for practical computations, in which parameter h = O(1/k 1.5 ) when linear finite elements are used. Note that this mesh parameter also represents an improvement. Indeed, an optimal error control for the H1 norm of ph would demand a finer mesh with h = O(1/k 2 ) [2].
4. Numerical examples In this section, the post-processing technique performance is illustrated by analyzing a number of typical situations comprising different domains and boundary conditions. In some of them, the exact solution is known which allows to confirm the rate of convergence predicted in the analysis presented before. The numerical experiments include linear, bilinear, quadratic and biquadratic finite elements. For all examples, the following fluid properties are assumed: c = 340 m s−1 and = 1.225 kg m−3 . Three performance indices, based on L2 norms, are used in the foregoing analysis: • eh = v − vh the difference between the exact and finite element solutions, • e∗ = v − vh∗ the difference between the exact and recovered solutions, • eh∗ = vh∗ − vh the difference between the recovered and finite element solutions, where vh = (−1/iw 0 )∇ph is an approximation of the velocity field resulting from the pressure finite element solution gradient. The first two indices are, indeed, the errors associated to approximating the velocity field by, respectively, the finite element solution gradient and the proposed recovery. The last one can be interpreted as an estimation of the former error by using the recovered velocity. Several numerical tests, concerning all the examples reported in this section, were performed in order to evaluate the sensitivity of the obtained results to different choices of the parameters 1 and 2 introduced in (3). Those tests revealed that results were not significantly affected by the values of the parameters. Thus, they were fixed equal to one in all simulations presented below. 4.1. Sound wave propagation in a planar duct The first situation treated here concerns sound wave propagation in a planar duct, which is schematically depicted in Fig. 2. The tube occupies a rectangular domain which dimensions are given by length L=1.0 m and width H = 0.2 m. The excitation surface is at tube’s left end (x = 0).
v0
V.n=0
L
Fig. 2. Duct geometry, boundary conditions and typical mesh.
x
470
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Two distinct problems, representing typical applications, are analyzed. In both of them, the propagated wave is essentially one-dimensional due the combination of symmetry with homogeneity of boundary conditions and by limiting the range of excitation by the cut off frequency. This one-dimensional setting is fairly interesting as it provides appropriate conditions, like the possibility of having exact solutions, to assess the features of the post-processing. In the first problem, a vertical uniform noise source is placed at the tube’s entry (x =0) and the interface of the exit (x = L) has a non-frequency-dependent fixed impedance. This is modeled by the following boundary conditions: p(0, y) = P 0,
vx (L, y) = vn (L, y) =
P0 . 0 c
(33)
Therefore, the x-component of the exact velocity is vx (x, y) =
P0 [cos(kx) − i sin(kx)] 0 c
(34)
which will be used to check the approximations furnished by the proposed recovery technique. In the second problem the excitation is modeled by a given normal velocity, a situation that typically arises in fluid–solid interaction, and the right end surface is considered to have infinity impedance. The boundary conditions are then given by v(0, y) = V0
v(L, y) = vn = 0
(35)
and the x-component of the exact velocity is thus obtained as vx (x, y) =
V0 sin[k(L − x)]. sin(kL)
(36)
The two problems involving the duct are solved with uniform meshes of 10, 14, 18, 22, 26 elements in the x direction and 2 in y. Linear and quadratic quadrilaterals and triangular elements are used, composing, respectively, arrangements of 2 × 2 and 2 × 4 macroelements. The first issue to be addressed now is the superconvergence of the pressure gradient approximation. It is investigated in a significant number of numerical tests which involve different elements, boundary conditions and excitation frequencies. The numerical evidences of the superconvergence occurrence are summarized in Figs. 5–8. They contain the results obtained for k = 5(f = 54 Hz) and k = 12(f = 649 Hz) and adopting linear and quadratic triangular and quadrilateral lagrangian elements. Those graphics, in which the notation Ti and Qj refer to triangles and quadrilaterals, present eh , the difference between the exact solution and ∇ph , using full integration (fi) and reduced integration(ri). The latter is related to special points where the superconvergence is likely to occur. As it can be observed, the expected behavior was obtained for quadrilaterals, in which the index introduced in (29) attained values between 0.7 and 1.0. The same was not true for triangles, although the obtained rates of convergence seemed to be slightly better for the reduced integration scheme. By no means it implies that superconvergence does not occur for those elements, although different sets of points were tried leading to the same behavior. Now, the convergence rates of the macroelement post-processing are checked. They confirm the a priori error estimates developed before. The results concerning the duct’s second problem are presented in Figs. 9–12. In those graphics, there are two performance indices being plotted, namely e∗ and eh . The former
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
471
y x
Fig. 3. Non-regular mesh, 26 × 26 quadrilateral elements.
0.2 0.2 P=Po 0.2
0.6
1.2
0.6
Fig. 4. Expansion chamber geometry, boundary conditions and typical mesh.
represents the error committed by approximating the velocity field using the post-processed solution. The second one aims to stress the improvement on the results by comparing the recovered results with the direct differentiation of the finite element solution ph . It is important to highlight that this was still observed for high values of the wave number, as k = 20 (f = 1082 Hz): Figs. 13 and 14. As it was stated previously, there is a pre-asymptotic behavior for high wave numbers which can be clearly seen in Fig. 13. Some numerical tests were carried out using non-uniform meshes leading to irregular elements, as shown in Fig. 3. The obtained post-processed velocity fields did not present high sensitivity to this nonregularity, although there was an expectation of non-optimal convergence rates due to a potential loss of superconvergence. Besides, the wave propagation becomes no longer parallel to the grid, which tends to increase the pollution effect.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -3.2 -3.4 -3.6 1.0
-3.8
log||e||
1.0
-4 1.1
-4.2 -4.4 2.0
-4.6 -4.8 0.7
fi=4 (Q4) ri=1(Q4) fi=3 (T3) ri=1(T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 5. 1D duct: superconvergence study at f = 54 Hz, linear elements.
fi=4 (Q4) ri=1 (Q4) fi=3 (T3) ri=1 (T3)
-2.8
-3
-3.2
1.2
log||e||
472
-3.4 1.2 1.6
-3.6 1.9
-3.8 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
-log (h)
Fig. 6. 1D duct: superconvergence study at f = 649 Hz, linear elements.
1.5
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-5 2.0
log||e||
-5.5 2.0
-6
-6.5 3.0
fi=9 (Q9) ri=4 (Q9) fi=6 (T6) ri=3 (T6)
-7
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 7. 1D duct: superconvergence study at f = 54 Hz, quadratic elements.
-4 -4.2 2.0
-4.4 -4.6 2.0
log||e||
-4.8 -5 -5.2
3.4
-5.4 fi=9 (Q9) ri=4 (Q9) fi=6 (T6) ri=3 (T6)
-5.6 -5.8 -6 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
-log (h)
Fig. 8. 1D duct: superconvergence study at f = 649 Hz quadratic elements.
1.5
473
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -4 -4.2 -4.4 1.0
-4.6
log (e)
1.0
-4.8 -5 -5.2 eh (Q4) e* (Q4) eh (T3) e* (T3)
-5.4
1.8 1.2
-5.6 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 9. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, linear elements.
-3.3 -3.4 -3.5 -3.6 -3.7
log (e)
474
1.3
-3.8 1.6 1.3
-3.9 -4 -4.1 -4.2
0.7
1.6
eh (Q4) e* (Q4) eh (T3) e* (T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 10. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, linear elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-5.5
log (e)
-6
2.0
-6.5
-7 2.3
-7.5
-8 0.7
eh (Q9) e* (Q9) eh (T6) e* (T6) 0.8
3.2
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 11. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, quadratic elements.
-4.4 -4.6 -4.8 -5
2.0
log (e)
-5.2 -5.4 -5.6 3.0
-5.8 -6
eh (Q9) e* (Q9) eh (T6)
-6.2
e* (T6)
-6.4 0.7
0.8
0.9
3.8
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 12. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, quadratic elements.
475
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2
log (e)
-2.5
-3 3.0
-3.5
-4 0.7
3.0
eh (Q4) e* (Q4) eh (T3) e* (T3) 0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 13. 1D duct: convergence study for macroelement post-processing at f = 1082 Hz, linear elements.
-4.2 -4.4 -4.6 -4.8
log (e)
476
2.0
-5 -5.2 3.0
-5.4 -5.6 -5.8 0.7
3.7
eh (Q9) e* (Q9) eh (T6) e* (T6) 0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 14. 1D duct: convergence study for macroelement post-processing at f = 1082 Hz, quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
477
Those tests were performed on 10 × 10, 14 × 14, 18 × 18, 22 × 22 and 26 × 26 meshes. To confirm once more the robustness of the proposed method, several tests were accomplished for the two proposed examples, but only results related to the duct problem involving a unitary impedance condition in the left end are graphically depicted, as similar conclusions could be obtained from all the tested examples. The 26 × 26 mesh is illustrated in Fig. 3, emphasizing the distortion of the elements. Figs. 19–22 report the performance of the post-processing for non-regular meshes. The error and its convergence rate are comparable to the ones obtained by adopting regular meshes, which are depicted in Figs. 15–18. The quality of an approximate solution can also be assessed by means of its derivatives. In Figs. 23 and 24, the divergence of the recovered velocity is compared to the exact solution and to the error of the divergence of the gradient of the finite element approximation as well. The better performance, for the analyzed frequency range, of the proposed methodology, when quadrilaterals were used, is observed once more. As it has been mentioned, there are some drawbacks in applying the finite element method to the Helmholtz equation. This can be illustrated by analyzing the dispersive characteristics of the obtained numerical wave. A valuable tool for it, often utilized by the modal analysis community, is the frequency response function (FRF). In general terms, it relates the acoustical response, here taken as e∗ and eh , to the excitation frequencies.
-3.2 -3.4 -3.6 1.0
-3.8
log (e)
1.0
-4 -4.2
1.7
-4.4
eh (Q4) e* (Q4) eh (T3)
-4.6 -4.8 -5 0.7
2.0
e* (T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 15. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, regular meshes: linear elements.
478
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-2.8
-3
log (e)
-3.2 0.9
-3.4 1.2
-3.6
-3.8
-4 0.7
1.96
eh (Q4) e* (Q4) eh (T3) e* (T3)
0.8
1.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 16. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, regular meshes: linear elements.
-5 2.0
log (e)
-5.5
-6 2.4
eh (Q9) e* (Q9) eh (T6) e* (T6)
-6.5
3.2
-7 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 17. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, regular meshes: quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
479
-4 -4.2 -4.4 2.0
log (e)
-4.6 -4.8 -5 -5.2
2.5
-5.4
eh (Q9) e* (Q9) eh (T6)
-5.6
e* (T6)
-5.8 -6 0.7
3.4
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 18. 1D duct: convergence study for macroelement post-processing at f = 6494 Hz, regular meshes: quadratic elements.
-2.8 eh (Q4) e* (Q4)
-3
1.0
log (e)
-3.2
-3.4
-3.6 1.97 -3.8
-4
-4.2 0.8
1
1.2
-log (h)
Fig. 19. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, non-regular meshes: linear elements.
480
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.4 -2.5
eh (Q4) e* (Q4)
-2.6 1.2
log (e)
-2.7 1.8 -2.8 -2.9 -3 -3.1 -3.2 0.8
1
1.2
-log (h)
Fig. 20. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, non-regular meshes: linear elements.
-4 eh (Q9) e* (Q9) 2.0
-4.5
log (e)
-5
-5.5
3.1
-6
-6.5 0.8
1
1.2
-log (h)
Fig. 21. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, non-regular meshes: quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
481
-3.2 eh (Q4) e* (Q4)
-3.4 -3.6
1.97
-3.8
log (e)
-4 -4.2
3.4
-4.4 -4.6 -4.8 -5 -5.2 0.8
1
1.2
-log (h)
Fig. 22. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, non-regular meshes: quadratic elements.
-2
log (div e)
-2.5
-3 div eh (Q4) div e* (Q4) div eh (T3) div e* (T3)
1.5
-3.5 1
1.05
1.1
1.15
1.2
1.25
1.3
1.35
1.4
1.45
-log (h)
Fig. 23. 1D duct: convergence study for macroelement post-processed divergence at f = 271 Hz, linear elements.
482
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.1 -2.2 -2.3
log (div e)
-2.4 -2.5 -2.6 1.2
-2.7 div eh (Q4) div e* (Q4) div eh (T3) div e* (T3)
-2.8 -2.9 -3 1
1.05
1.1
1.15
1.2
1.25
1.3
1.35
1.4
1.45
-log (h)
Fig. 24. 1D duct: convergence study for macroelement post-processed divergence at f = 271 Hz, linear elements.
Those FRFs corresponding to meshes of 10 and 26 linear elements in the x direction are presented in Figs. 25–28. The following conclusions could be drawn. • For low wave numbers, e∗ and eh present a good and similar behavior, what was already expected from the analysis developed in the previous section. • For high wave numbers, the post-processing presents a deterioration in its performance, but it still represents a better approximation of the acoustic velocity than the gradient of the finite element solution. • All FRFs manifest the same tendency of growth deviation between approximations and exact value when the wave number tends to magnify. This phenomenon is in accordance with the increasing influence of the pollution mentioned before. 4.2. Two-dimensional cavity Once again, two situations are treated involving the same domain form: a two-dimensional cavity. In the first problem, the cavity occupies the region (0 x 10, 0 y 10) and normal velocity is set to be zero at boundaries, therefore the acoustic mode corresponding to natural frequency of 49.02 Hz, is given by the closed form x y x y 1
v(x, y) = sin cos ex + cos sin ey , (37) i0 w 5 5 5 5 5 5 where ex and ey stand for the ordinary vectorial basis in the plane.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 100 eh (T3) eh (Q4) e* (T3) e* (Q4)
90 80 70
error%
60 50 40 30 20 10 0 0
2
4
6
8
10 k
12
14
16
18
20
Fig. 25. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.1 m.
100 eh (T3) eh (Q4) e* (T3) e* (Q4)
90 80 70
error%
60 50 40 30 20 10 0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 26. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.38 m.
483
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
45 eh (T6) eh (Q9) e* (T6) e* (Q9)
40 35
error%
30 25 20 15 10 5 0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 27. FRFs relating the global error eh and e∗ to the wave number k, quadratic elements: h = 0.1 m.
3
eh (T6) eh (Q9) e* (T6) e* (Q9)
2.5
2
error%
484
1.5
1
0.5
0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 28. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.38 m.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
485
The same situation was analyzed in [19] where a finite element error estimator for the acoustic problem was proposed. The potentiality of using the recovery technique proposed here as the basis of an error estimator is now briefly addressed. This estimator is simply eh and its performance is assessed by comparing it to eh∗ . The indices will be referred to as estimated error and exact error and their distribution along the domain is presented in Figs. 29 and 30. In order to highlight the good quality of the results, the difference between the exact and estimated error is plotted in Fig. 31. The second problem consists in a unity cavity [7] with D comprising all edges of the square and N = 0. The imposed Dirichlet boundary condition is given as p(x, y) = exp[ik(x cos + y sin )],
(x, y) ∈ D .
(38)
Therefore the exact solution to this problem is a plane wave propagating at angle to the cartesian x-axis, and is given by p(x, y) = exp[ik(x cos + y sin )],
(x, y) ∈
(39)
and v(x, y) = p(x, y)ik exp[− sin ex + cos ey ],
(x, y) ∈ .
(40)
The solution to this problem was sought on uniform meshes of bilinear finite elements. Convergence rates are shown in Fig. 32 for a single =45◦ value and frequencies of f =649 and 1082 Hz for several meshes.
Fig. 29. 2D cavity: exact error (eh ) distribution.
486
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Fig. 30. 2D cavity: estimated error (eh∗ ).
Fig. 31. 2D cavity: difference between exact and estimated error distributions.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
487
-2
-2.5
log (e)
-3 1.5 -3.5
-4 eh (f=649 Hz) e* (f=649 Hz) eh (f=1082 Hz)
-4.5
3.5
e* (f=1082 Hz) -5 1
1.2
1.4
1.6
1.8
2
-log (h)
Fig. 32. 2D cavity: convergence study for = 45◦ .
Two performance indices are plotted, namely e∗ and eh . Please note that convergence rates were visibly better for the post-processing and the approximate solution didn’t present any earnings, independent of the wave number. A similar behavior is presented in Fig. 33, in which and k were varied and the mesh was fastened with 40 × 40 elements. It is worth to observe the deterioration on the post-processing performance with the increasing of k due to the pollution effect. That happens for all searched wave propagation direction , although for high wave numbers as k = 20, the proposed recovery still presents good solutions for = 45◦ .
4.3. Expansion chamber This is a more challenging situation, as it deals with a two dimensional real like problem. It consists in an expansion chamber with a perforated outlet pipe, schematically shown in Fig. 4. The excitation stems from a noise source with p(0, y) = 2.0 at the left-hand side. At the right-hand side, the interface between the interior fluid and the exterior is modeled by a constant impedance. A comparison of the imaginary and real components in the middle section’s chamber between vh and v∗ for k =10 (f =541 Hz) is presented in Figs. 34–37. From them it is possible to figure out the improvement on using the proposed post-processing for obtaining the acoustic velocity field. Those results show that the recovered velocity was able to capture significative variations due to abrupt changes in the domain geometry.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.5 -3 -3.5
eh (theta = 22.5) e* (theta = 22.5) eh (theta = 45) e* (theta = 45) eh (theta = 75) e* (theta = 75)
log (e)
-4 -4.5 -5 -5.5 -6 -6.5 -7 1
2
3
4
5
6
7
8
9
10
K
Fig. 33. 2D cavity: post-processing’s numerical performance for several ’s.
x 10-3
1.5
v vh v*
1
0.5 real velocity (ur)
488
0
-0.5
-1
-1.5 0
0.5
1
1.5
2
x
Fig. 34. Expansion chamber: horizontal real velocity component along the center line.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 0.01 v vh v*
0.008
imaginary velocity (ui)
0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 0
0.5
1
1.5
2
x
Fig. 35. Expansion chamber: horizontal imaginary velocity component along the center line.
x 10-4 6 v vh v*
4 2
real velocity (vr)
0 -2 -4 -6 -8 -10
0
0.5
1
1.5
2
x
Fig. 36. Expansion chamber: vertical real velocity component along the center line.
489
490
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
2
x 10-3
1
imaginary velocity (vi)
0 -1 -2 -3 -4 v vh v*
-5 -6 -7 0
0.5
1
x
1.5
2
Fig. 37. Expansion chamber: vertical imaginary velocity component along the center line.
5. Concluding remarks The macroelement post-processing recovery technique introduced in [15] was extended here for the steady-state acoustic problem governed by the Helmholtz equation. The obtained analytical and numerical results have demonstrated the ability of this method for handling the inherent numerical difficulties linked to the finite element approximation of this equation. The proposed method achieves good accuracy for the velocity acoustic field at low computational cost. The convergence of the method is proved leading to a priori error estimates, which are confirmed by several numerical examples usually explored in the literature. The computation implementation is very simple, only requiring a standard finite element code and is amenable for parallel computing. The velocity post-processing has been successfully used in noise source estimation [12], enhancing the numerical performance of the iterative algorithm used to solve the resulting inverse problem, and also as a local error estimator [14] to be explored in adaptive finite element analysis. At this point, it is important to mention that the post-processing could be regarded in a different perspective of the one adopted here. It could be undertaken not based on a standard finite element framework. This means that either vh∗ could be approximated by using different basis functions, such as harmonic functions [6,20], or by using as input data (the pressure field) not a finite element solution of
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
491
Helmholtz equation. In this last case, a typical application would be obtaining a smooth representation of the velocity field from a cloud of pressure experimental values.
Acknowledgements We gratefully acknowledge the partial support from the Conselho Nacional de Desenvolvimento e Pesquisa (CNPQ).
References [1] Ph. Bouilarrd, F. Ihlenburg, Error estimation and adaptivity for the finite element method in acoustics: 2D and 3D applications, Comput. Methods Appl. Mech. Eng. 176 (1999) 147–163. [2] F. Ihlenburg, Finite element analysis of acoustic scattering, Appl. Math. Sci. 132 (1998). [3] A. Deraemaeker, I. Babuska, P. Bouillard, Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Meth. Eng. 46 (1999) 417–499. [4] I. Harari, T.J.R. Hughes, Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains, Comput. Methods Appl. Mech. Eng. 97 (1992) 103–124. [5] R. Baus˜ys, N.E. Winberg, Adaptive finite element strategy for acoustics problems, J. Sound Vibr. 226 (5) (1999) 905– 922. [6] C. Farhat, I. Harari, L.P. Franca, The discontinuous enrichment method, Comput. Meth. Appl. Mech. 190 (2001) 6455– 6479. [7] A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Meth. Eng. 49 (2000) 399–419. [8] S. Suleau, Ph. Bouillard, One-dimensional dispersion analysis for the element-free Galerkin method for the Helmholtz equation, Int. J. Numer. Meth. Eng. 47 (2000) 1169–1188. [9] L.L. Thompson, P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Meth. Eng. 38 (1995) 371–397. [10] K.S. Peat, K.L. Raithi, A finite element analysis of the convected acoustic wave motion in dissipative silencers, J. Sound Vibr. 184 (3) (1995) 529–545. [11] T. Lin, D.L. Russel, A superconvergent method for approximating bending moment of elastic beams with hysteris damping, Appl. Numer. Math. 38 (2001) 145–165. [12] R.H.P. Romero, F.A. Rochinha, A post-processing technique applied to error estimation of FEM approximations of Helmholtz equation, Proceedings of the Seventh US National Congress on Computational Mechanics, vol. 107, New Mexico, USA, 2003. [13] R.H.P. Romero, A post-processing scheme to the Helmholtz equation and applications, Doctor in Science Thesis, Federal University of Rio de Janeiro, Brazil, 2003 (in Portuguese). [14] R.H.P. Romero, F.A. Rochinha, Estimation of noise sources using an enhanced finite element formulation, Proceedings of the Seventh US National Congress on Computational Mechanics, vol. 612, New Mexico, USA, 2003. [15] A.F.D. Loula, F.A. Rochinha, M.A. Murad, Higher-order gradient post-processings for second-order elliptic problems, Comput. Methods Appl. Mech. Eng. 128 (1995) 361–381. [16] A.F.D. Loula, M.A. Murad, F.A. Rochinha, 3D gradient post-processings schemes for second order elliptic problems, in: Proceedings of the IV World Conference on Computational Mechanics, Buenos Aires, Argentina, 1998. [17] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Numer. Meth. Eng. 33 (1992) 1331–1364.
492
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
[18] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 2: error estimates and adaptivity, Int. J. Numer. Meth. Eng. 33 (1992) 1365–1382. [19] Ph. Bouilarrd, J.F. Allard, G. Warzée, Superconvergent patch recovery technique for the finite element method in acoustics, Comm. Num. Methods Eng. 12 (1996) 581–594. [20] C. Farhat, I. Harari, U. Hetmaniuk, A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime, Comput. Meth. Appl. Mech. 192 (2003) 1389–1419.
A higher-order velocity post-processing for the Helmholtz equation Rivânia H. Paulinoa , Fernando Alves Rochinhab,∗ a Department of Mechanical Engineering, Federal University of Espirito Santo, UFES, Vitória, ES, Brazil b Department of Mechanical Engineering, Federal University of Rio de Janeiro, UFRJ, Ilha do Fundao, Cx. Postal 68503, Rio
de Janeiro, RJ 21945-970, Brazil Received 29 October 2003; accepted 20 May 2004 Available online 25 December 2004
Abstract A macroelement post-processing recovery technique of the velocity field based on residuals of the equilibrium equation and irrotational conditions for the Helmholtz equation is proposed. The derivatives are recovered by solving local variational problems, that are based on the residuals aforementioned and involve superconvergent derivatives of the acoustical pressure field on special points. Improved accuracy for derivative finite element approximation is obtained with low computational cost and easy implementation. A number of comprehensive simulations is presented in order to show the main features of the proposed methodology. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Acoustics; Post-processing; Flux recovery
1. Introduction In acoustic problems, the steady-state sound response in a non-viscous medium is governed by the Helmholtz equation [1]. One major problem in numerically solving this equation by means of the finite element method (FEM) is linked to a potential loss of ellipticity with an increasing excitation frequency. This degradation of the standard Galerkin FEM approximation is specially observable in the convergence behavior of the acoustic velocity, a field directly related to the gradient of the primal variable. ∗ Corresponding author. Tel.: +55 21 25628315.
E-mail address: [email protected] (F.A. Rochinha). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.05.005
460
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Classical Galerkin FEM provides good phase and amplitude accuracy as long as the mesh is fine enough with respect to the wave number in the propagation region. However, “fine enough” is often too expensive for adequate resolution, even for moderate wave numbers. Hence, one of the main concerns in acoustic finite element analysis is the adequacy of the finite element mesh. Acousticians often use the, so-called, “rule of the thumb” [2] which prescribes a relation between the minimal number of elements and the wave number. This rule only takes into account the capacity of the interpolation given by a specific mesh and is not based on stability aspects. Therefore, as it is well documented in the literature [1–4], one can obtain very bad results even when the finite element mesh obeys this rule. Those approximation problems have been addressed in a number of publications leading to several solutions like the ones contained in [1,5–9]. The numerical phenomena was analyzed in [2] where the following error estimate for the pressure gradient is presented, when linear elements are used. eh C1 kh + C2 k 3 h2 ,
kh < 1
(1)
in which constants C1 and C2 are independent of wave number k and element size h. The first term on the right-hand side of Eq. (1) represents the interpolation error and the second one represents the pollution error. In order to keep the local error constant, maintaining kh constant is enough. That is the “rule of the thumb” which corresponds to taking a certain number of elements per wavelength. As one can see from Eq. (1), this is not sufficient to maintain the pollution error under control. Thus, even following the heuristic rule, the pollution error increases with k. In some applications, ranging from very specific situations like the absorption analysis in silencers [10] to more general scenarios involving inverse problems [11,12] or error estimators [1,5,13,14], the velocity field which is directly linked to the pressure gradient is of main interest. As previously mentioned, the differentiation of pressure computed by standard finite element approximation leads to poor results. Thus, alternative formulations have been developed aimed at improving the velocity field. Here, a macroelement post-processing approach [15,16] is developed to recover the gradient of the pressure with improved accuracy. The macroelement post-processing is constructed over subregions of the finite element mesh and is based on the superconvergence of the derivatives of the pressure finite element solution at special points combined with a residual form of a balance equation and the irrotationality condition. Numerical analysis and error estimates are presented which demonstrate that higher asymptotic convergence rates are obtained with the proposed approach. An outline of the paper follows. A brief review of the basic problem leading to the Helmholtz equation is presented in Section 2. The macroelement post-processing is introduced and analyzed in Section 3. Numerical experiments confirming the predicted rates of convergence are reported in Section 4. Finally, Section 5 contains some concluding remarks. 2. Problem statement and finite element formulation 2.1. Continuous problem The expression acoustic problem used here refers to the sound propagation in a non-viscous fluid in equilibrium. All non-linear effects are neglected and only a time-harmonic response is considered leading to the Helmholtz equation, in which the fluctuation of the sound pressure is seek. Although it is a very
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
461
well-known modelling, the basic equations will be briefly reviewed below in order to help the developing of the recovery strategy which will be proposed later on. The modelling relies on three equations, namely: continuity, Euler and state. They involve the velocity field V=V(x, t), the pressure P =P (x, t) and the density ¯ = ¯ (x, t). The first two equations are presented right below in the same order they were mentioned in the beginning of this paragraph j¯
+ div(¯V) = 0,
jt jV
1 + (V div(∇V)) + ∇P = 0, jt ¯
(2) (3)
where div( ) = j( )/jx + j( )/jy is the divergence operator and ∇( ) = (j( )/jx)ex + (j( )/jy)ey is the gradient with ex and ey as the unitary vectors x and y directions. Assuming that the sound propagation will produce only slight perturbations on the flow, the following fluctuations of the basic fields are introduced: = ¯ − 0 , v = V and p = P − p0 , where the subindex 0 refers to the uniform state. With this in mind, the non-linear terms contained in the above equations are linearized giving rise to j jt 0
+ 0 div v = 0, jv jt
+ ∇p = 0.
(4) (5)
The above linearized equations are complemented by the linear state equation p = c2 , where c stands for the sound velocity in the fluid. As it was mentioned earlier, only harmonic solutions will be treated here. Therefore, all fields can be recast as f (x, t) = fˆ(x)eiwt , where fˆ(x) is a stationary function and w is the circular frequency. Thus, Eq. (5) assumes the following form: vˆ = −
1 ∇ p. ˆ iw0
(6)
Eq. (6) will be referred to as the constitutive equation due to its resemblance to the constitutive relations in structural dynamics. The combination of Eq. (4) and the state equation leads to iw pˆ + div vˆ = 0. c 2 0
(7)
Moreover, directly from (6) one achieves the irrotationality condition rot(ˆv) = 0,
(8)
where rot stands for the rotational operator. The last three equations form the basis for the proposed velocity recovery method. For the remainder of this analysis, the symbol ˆ. will be omitted because all the variables appearing till the end are the stationary value of the fields. Therefore, there is no need to distinguish in the adopted notation.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
462
Finally, the Helmholtz equation can be retrieved by combining Eqs. (6) and (5) leading to the following mathematical problem Find p : → C such that ∇ 2 p + k 2 p = 0,
(9)
where wave number k is defined by k=
(10)
c
and ∇ 2 stands for the Laplacian operator and C is the set of complex numbers. In order to have a complete description of the wave propagation phenomena, it is necessary to identify the appropriate boundary conditions, which, in general terms, are set below. Only internal acoustic will be addressed here, so boundary conditions associated to infinite domains like the absorbing ones will not be included. • Prescribed pressure on D (Dirichlet boundary condition) p = p˜
on D .
(11)
• Prescribed acoustic velocity on N (Neumann boundary condition) vn = v˜ n
on N .
(12)
• Fixed impedance in the boundary R (Robin boundary condition) ∇p · n = −ikc0 An p
on R
(13)
where An is the admittance coefficient, the inverse of impedance in the boundary. 2.2. Variational formulation and finite element discretization In order to solve the interior acoustic problem introduced above by means of the FEM, the following variational formulation, based on the primal variable p, is stated. Problem P. Find p ∈ V such that 2 ∇p.∇q d − k pq d − ik 0 c ∇p · nq d = 0 ∀q ∈ V , (14)
N ∪ R
where . stands for the ordinary scalar product in the complex plane and the function space of admissible pressures V is defined as V = {p ∈ H1 () | q = 0 on D }, where H1 () denotes the Sobolev space of complex-valued functions 2 jf 1 2 H () = f : → C, |f (x)| dx < ∞, (x) dx < ∞ . jx H1 () is a subset of L2 () the space of square integrable complex functions.
(15)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
463
The Galerkin finite element approximation of Eq. (14) is obtained by considering a finite dimensional space Vh ⊂ V , and then seeking an approximate solution to (14). Let Shs be a quadrilateral or triangular C0 Lagrangian finite element subspace of H1 () of degree s. The discrete version of Problem P in Vh = Shs ∩ V reads Problem Ph : Find ph ∈ Vh , such that 2 ∇ph .∇qh d − k ph qh d − ik 0 c ∇ph · nqh d = 0 ∀qh ∈ Vh , (16)
N ∪R
where h represents the mesh parameter. The above discrete formulation leads to the matrix equation (K + ik cC − k 2 M)q = −ik cb,
(17)
where K, M are often referred to as rigidity and mass matrices and q is the vector of Np nodal values of the discrete pressure ph . Besides, C and b are associated to boundary terms. This terminology is motivated by the similarity of (17) to the structural dynamic equations. The components of these matrices and vectors are given by Kij = ∇ j · ∇ i d, Mij = j i d, Cij =
R
(∇ j · n)i d,
bi =
N
(v · n)i d
with j (j = 1, . . . , Np ) denoting the standard finite element basis. The convergence analysis of the above finite element problem is addressed in [2] where the emphasis is placed on the pollution effect. That numerical misbehavior reflects the lack of ellipticity of the finite element formulation in the presence of high wave numbers. In many finite element codes an heuristic strategy of choosing the mesh’s size as a function of k, often referred to as the “rule of the thumb”, is adopted in order to keep track of the high oscillatory behavior of p. The following estimate, presented in [2], for a one-dimensional prototype situation, demonstrates that even with the adoption of this rule the pollution effect strongly disturbs the quality of the solution |ph − p|1 C1 kh + C2 k 3 h2 , |p|1
kh < 1 (“ rule of the thumb” ),
(18)
where |.|1 is the H1 semi-norm and C1 and C2 are positive numbers independent of the mesh. The second term of the above bound dominates the error in the pre-asymptotic range of the finite element approximation, which implies that for reasonable meshes, on an engineering perspective, the standard analysis might entail unreliable results.
3. Macroelement post-processing technique Usually, after solving the primal problem introduced in the previous section, the velocity field is obtained by using (6). This direct approach leads to lower-order accurate derivatives when compared to
464
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
the primal variable. Nevertheless, the accuracy of the derivatives is often of main interest in engineering problems. This motivates the development of the so-called recovery strategies, in which a more accurate approximation of the derivative variable is sought. Several approaches are found in literature [17,18] and it is worth mentioning that some of the most important ones have already been tested on the Helmholtz equation with relative success [1,5,19], in the context of finite element adaptivity. The macroelement post-processing technique proposed in [15] will be extended to the acoustic problem introduced earlier. The basic idea of the method is to find a better approximation for the derivative (∇ph ), and consequently for the acoustic velocity field, by solving a local variational problem on each macroelement, understood as the union of a set of neighboring elements with common edges. This problem involves least-squares residuals of the balance equation (7) and of the irrotationality condition (8) and relation (6) at special points, which are likely to present superconvergence [15,19]. The inclusion of the residuals of the two equations is twofold, it enhances the stability of the macroelement problem and includes explicitly in the formulation conditions to be obeyed by the actual velocity field. This last aspect insures that, in a least-squares sense, the approximated velocity reproduces some important features presented by its continuous counterpart. In order to establish the proposed method, one can suppose that the domain is decomposed into macroelements not necessarily disjoints and define MP rh ⊂ L2 () and MQrh ⊂ L2 (), respectively, triangular and quadrilateral lagrangian finite element spaces of piecewise polynomial of degree r on each element and class C 0 in each macroelement but discontinuous on the macroelement boundaries. Considering Zh = MP rh × MP rh × MP rh or Zh = MQrh × MQrh × MQrh , the following residual form for the problem of retrieving a better approximation of the three-dimensional acoustic velocity field is proposed Problem ME h : Given ph ∈ Vh solution of Problem Ph , find vh∗ ∈ Zh , such that iw i ∗ 2 bh (vh , uh ) = ∀uh ∈ Zh , ∇ph , uh − 1 h ph , div uh (19) w 0 0 c 2 h where the sesquilinear form bh (., .) is defined as bh (vh∗ , uh ) = (vh∗ , uh )h + 1 h2 (div vh∗ , div uh ) + 2 h2 (rot vh∗ , rot uh ), where vh∗ is the post-processed velocity. A compact notation was adopted, in which (., .) represents the L2 () product and (., .)h denotes that this term is evaluated by a suitable integration rule, taking into account only specific points where superconvergence of the derivatives are supposed to occur. Taking, for instance, macroelements composed by four bilinear quadrilaterals in a 2 × 2 arrange, as depicted in Fig. 1(a), and evaluating (·, ·)h with one Gauss point per element, one has a configuration similar to the biquadratic element. The scalars 1 and 2 are positive real parameters that might be adjusted in order to have better convergence performance. Whenever N = 0, the velocity at this part of the boundary can be strongly enforced, and Problem ME h is slightly modified, incorporating the boundary condition in the discrete space Zh . Concerning the computational implementation, the macroelement post-processing leads to a number of similar and independent linear problems, each one corresponding to a macroelement. The degrees of freedom associate to each node are the real and imaginary parts of the velocity components. Therefore, for two-dimensional problems, using a macroelement formed by four elements, which will be adopted in some numerical examples of this article, the total number of unknowns is 36.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
(a)
465
(b)
Fig. 1. Two different macroelement configurations composed of bilinear quadrilaterals. (•), nodal points; (+), superconvergent points; (a) 2 × 2 arrange; (b) 2 × 1 arrange.
This computational implementation is rather simple, requiring only a standard finite element code in order to obtain the matrix problem to be solved at each macroelement. It is important to point out that parallel computation could be explored as the proposed method leads to independent problems covering all the domain. In such case, groups formed by a number of macroelements could be assigned to different processors. A reasonable alternative to a post-processing scheme in order to obtain accurate flux approximations is provided by mixed formulations [15]. Unfortunately so far, in the authors knowledge, there are no successful mixed methods for the Helmholtz equation encompassing general boundary conditions reported in the literature. 3.1. Convergence analysis Constructing a well-posed Problem ME h depends on the number of elements in the macroelement, on the choice of the finite element space Zh and on the number of superconvergence points adopted to evaluate the constitutive Eq. (6). In the present section those topics are investigated leading to the post-processing’s convergence analysis. It is also important to mention that no regularity issue need to be addressed, though the presence in the formulation of the differential operators div and rot, as the postprocessed velocity, over each macroelement, is a finite element function of, at least, class C 0 , therefore belonging locally to H1 . In general, convergence of a numerical approximation is assessed by combining two main issues: consistency and stability. As consistency is guaranteed by the construction of the formulation, stability remains the only aspect to be addressed in order to analyze the post-processing’s convergence. A stability check can be performed by means of a simple numerical test involving an eigenvalue problem, as proposed in [15]. In the present context, this methodology was applied to several different configurations of macroelements, that are summarized in Table 1 for the situation described in the first numerical example of the next section.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
466
Table 1 Stability study of macroelement configurations Arrangement
Elements
min
max
Stability
(a) 1 × 1
Q4 Q9 T3 T6 Q4 Q9 T3 T6 Q4 Q9 T3 T6
— 0.4110 — 0.1448 0.2649 0.7516 0.3291 0.2777 0.7500 0.9009 0.8320 0.4857
— 140.01 — 213.680 26.725 139.430 40.093 178.220 26.721 137.64 35.928 164.40
No Yes No Yes Yes Yes Yes Yes Yes Yes Yes Yes
(b) 1 × 1 (c) 2 × 1 (d) 1 × 2 (e) 2 × 2 (f) 2 × 4
The numerical results which will be presented later demonstrate that the proposed post-processing method has good accuracy and robustness. It is also remarkable that when compared to the direct calculation of the velocity field, a better performance is achieved. Stability of the proposed formulation is proved with the aid of the following associated generalized eigenvalue problem j
j
bh (h , uh ) = j (h , uh )
∀ uh ∈ Zh |ME (j = 1, . . . , N )
(20)
j
where j and h are, respectively, an eigenvalue and its correspondent eigenvector and Zh |ME corresponds to the restriction of space Zh to a single macroelement. Besides, N stands for the number of degrees of freedom associated to each macroelement. The j ’s in (20), which are all real and non-negative due to the symmetry of the involved operators, can be easily computed for different macroelement configurations by using standard eigenvalue algorithms, within a finite element code. Table 1 reports the minimum and the maximum eigenvalues for several different arrangements. The post-processed velocity vh∗ can be expressed as a linear combination of the eigenvectors yielding N N ∗ ∗ ∗ i ∗ i (21) vi h , v i h , bh (vh , vh ) = bh i=1
i=1
vi∗
are the components of the flux variable in the basis formed by the eigenvectors. From the where bilinearity of the bh (., .), considering (20) and taking profit of the orthogonality among eigenvectors, one has bh (vh∗ , vh∗ ) =
N i=1
i (vi∗ )2 (i , i )
(22)
which in turn implies ∗ 2 bh (vh∗ , vh∗ ) Min j vh .
(23)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
467
If Min j , the lowest eigenvalue, is strictly positive, inequality (23) implies that bh is a positive definite operator, which entails the stability of the formulation in a local sense. The stability on the whole domain follows by summing (23) over all macroelements. The coercivity constant given by the aforementioned eigenvalue depends on the aspect ratio of the elements contained on the macroelement but not on the mesh’s parameter h. This is assessed computing the matrices related to the eigenvalue problem using integrals performed with the aid of on a reference macroelement of unity dimensions. In order to derive error estimates for the proposed methodology, the following norm, equivalent to the canonical one in Zh , is introduced |uh | = bh (uh , uh )1/2 .
(24)
The stability results stated above are now used to prove the post-processing’s convergence, leading to a priori error estimates. Those results are presented in Lemma 2, which relies on Lemma 1. Lemma 1. Let p be the solution of the Problem P, v given by (6) and ph and vh∗ be their respective approximations by solving problems Ph and ME h . Then, considering macroelements leading to stable formulations, there exists a positive number C independent of h such that v − vh∗ C(v − uh + |v − uh | +
1 |p − ph |h + khp − ph ) k
∀ uh ∈ Zh ,
(25)
where | · |h represents the semi-norm endowed by the special integration (., .)h . A slight abuse of notation was adopted as |||.||| is only proved to be a norm in Zh . Therefore the term |||v − uh |||, as v − uh ∈ / Zh , must be understood as a compact notation for bh (v − uh , v − uh ). Proof. The variational consistency referred to above, which states that the true solution satisfies the discrete problem, is expressed as iw i 2 bh (v, uh ) = ∀ uh ∈ Zh ∇p, uh − 1 h p, div uh (26) w0 0 c2 h which along with the stability estimate (23) leads to vh∗ − uh 2 Cbh (vh∗ − uh , vh∗ − uh ) C(bh (v − uh , vh∗ − uh ) + bh (vh∗ − v, vh∗ − uh )) i ∗ ∗ C bh (v − uh , vh − uh ) + ∇(ph − p), vh − uh w h iw −1 h2 (ph − p), div (vh∗ − uh ) c 2 1 C |||v − uh ||| ||vh − uh || + |ph − p|h ||vh − uh || w w + ˜ 1 h 2 ||ph − p|| ||vh − uh || , c
(27)
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
468
where ˜ 1 stands for the former parameter 1 multiplied by a constant that comes from the inverse estimate used to obtain the third term of the last inequality. Hence, from (27) 1 k ∗ ˜ vh − uh C |||v − uh ||| + |p − ph |h + 1 h ||p − ph || (28) kc c which yields the desired result after using the triangular inequality. Now the convergence result is stated for the macroelement post-processing. Lemma 2. Under the same hypotheses of Lemma 1 and assuming that the continuous fields p and v are regular enough and the existence of a superconvergence behavior of the derivative expressed as |p − ph |h (C1 k r hr+ + C2 k 2r+1 h2r+ )|p|r
0 1
(29)
there exist positive numbers C1 and C2 , independent of h, such that v − vh∗ (C1 (kh)r h + C2 k(kh)2r+1 ) |p|r
kh < 1.
(30)
Proof (Error estimate). Considering the classical finite element interpolation results inf |||v − uh ||| Chr+1 |v|r+1 ,
uh ∈Zh
inf ||v − uh || Chr+1 |v|r+1 ,
uh ∈Zh
the approximation errors established for the one-dimensional case in [2] considering kh < 1 p − ph (C1 (kh)r + C2 k 2r+1 h2r )|p|r ,
(31)
the superconvergence hypothesis (29) and that for an oscillatory solution |p|r /||p|| = O(k r ), one has v − vh∗ (Ck r+1 hr+1 + C1 (k r−1 hr+ + (kh)2r+ ) + C2 ((kh)r+1 + k 2(r+1) h2r+1 ))|p|r .
(32)
Now, clustering the similar terms, the desired result is achieved. There are a few remarks to be made about the above lemma. First, it is worth calling attention to the two parcels of the final error estimate. Similarly to the analysis presented in [2], pre and asymptotic convergence behaviors are observed. The second parcel is a direct consequence of the pollution effect. The second point to be highlighted is about the error estimate for pressure field. It is only rigorously proved to the one-dimensional case, but a number of numerical simulations reported in the literature have demonstrated that a similar behavior is observed for two-dimensional and three-dimensional situations. Therefore, the obtained estimation for the post-processed velocity can be also valid to those cases. The superconvergence hypothesis used in Lemma 1 is also inspired in numerical observations not relying on any sort of theoretical developments. Anyway, this sort of behavior is not only confirmed by the numerical experiments in Section 4, but also in [5,19]. Estimate (32), in the asymptotic range, represents a higher convergence rate when compared to the velocity obtained by the direct differentiation of the finite element approximation of the pressure field.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
469
This is fairly confirmed by numerical examples presented in the next chapter. On the other hand, the pre-asymptotic range is a result of the pollution contained in ph . This error on the post-processed velocity would be controlled by a mesh, probably, not suitable for practical computations, in which parameter h = O(1/k 1.5 ) when linear finite elements are used. Note that this mesh parameter also represents an improvement. Indeed, an optimal error control for the H1 norm of ph would demand a finer mesh with h = O(1/k 2 ) [2].
4. Numerical examples In this section, the post-processing technique performance is illustrated by analyzing a number of typical situations comprising different domains and boundary conditions. In some of them, the exact solution is known which allows to confirm the rate of convergence predicted in the analysis presented before. The numerical experiments include linear, bilinear, quadratic and biquadratic finite elements. For all examples, the following fluid properties are assumed: c = 340 m s−1 and = 1.225 kg m−3 . Three performance indices, based on L2 norms, are used in the foregoing analysis: • eh = v − vh the difference between the exact and finite element solutions, • e∗ = v − vh∗ the difference between the exact and recovered solutions, • eh∗ = vh∗ − vh the difference between the recovered and finite element solutions, where vh = (−1/iw 0 )∇ph is an approximation of the velocity field resulting from the pressure finite element solution gradient. The first two indices are, indeed, the errors associated to approximating the velocity field by, respectively, the finite element solution gradient and the proposed recovery. The last one can be interpreted as an estimation of the former error by using the recovered velocity. Several numerical tests, concerning all the examples reported in this section, were performed in order to evaluate the sensitivity of the obtained results to different choices of the parameters 1 and 2 introduced in (3). Those tests revealed that results were not significantly affected by the values of the parameters. Thus, they were fixed equal to one in all simulations presented below. 4.1. Sound wave propagation in a planar duct The first situation treated here concerns sound wave propagation in a planar duct, which is schematically depicted in Fig. 2. The tube occupies a rectangular domain which dimensions are given by length L=1.0 m and width H = 0.2 m. The excitation surface is at tube’s left end (x = 0).
v0
V.n=0
L
Fig. 2. Duct geometry, boundary conditions and typical mesh.
x
470
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Two distinct problems, representing typical applications, are analyzed. In both of them, the propagated wave is essentially one-dimensional due the combination of symmetry with homogeneity of boundary conditions and by limiting the range of excitation by the cut off frequency. This one-dimensional setting is fairly interesting as it provides appropriate conditions, like the possibility of having exact solutions, to assess the features of the post-processing. In the first problem, a vertical uniform noise source is placed at the tube’s entry (x =0) and the interface of the exit (x = L) has a non-frequency-dependent fixed impedance. This is modeled by the following boundary conditions: p(0, y) = P 0,
vx (L, y) = vn (L, y) =
P0 . 0 c
(33)
Therefore, the x-component of the exact velocity is vx (x, y) =
P0 [cos(kx) − i sin(kx)] 0 c
(34)
which will be used to check the approximations furnished by the proposed recovery technique. In the second problem the excitation is modeled by a given normal velocity, a situation that typically arises in fluid–solid interaction, and the right end surface is considered to have infinity impedance. The boundary conditions are then given by v(0, y) = V0
v(L, y) = vn = 0
(35)
and the x-component of the exact velocity is thus obtained as vx (x, y) =
V0 sin[k(L − x)]. sin(kL)
(36)
The two problems involving the duct are solved with uniform meshes of 10, 14, 18, 22, 26 elements in the x direction and 2 in y. Linear and quadratic quadrilaterals and triangular elements are used, composing, respectively, arrangements of 2 × 2 and 2 × 4 macroelements. The first issue to be addressed now is the superconvergence of the pressure gradient approximation. It is investigated in a significant number of numerical tests which involve different elements, boundary conditions and excitation frequencies. The numerical evidences of the superconvergence occurrence are summarized in Figs. 5–8. They contain the results obtained for k = 5(f = 54 Hz) and k = 12(f = 649 Hz) and adopting linear and quadratic triangular and quadrilateral lagrangian elements. Those graphics, in which the notation Ti and Qj refer to triangles and quadrilaterals, present eh , the difference between the exact solution and ∇ph , using full integration (fi) and reduced integration(ri). The latter is related to special points where the superconvergence is likely to occur. As it can be observed, the expected behavior was obtained for quadrilaterals, in which the index introduced in (29) attained values between 0.7 and 1.0. The same was not true for triangles, although the obtained rates of convergence seemed to be slightly better for the reduced integration scheme. By no means it implies that superconvergence does not occur for those elements, although different sets of points were tried leading to the same behavior. Now, the convergence rates of the macroelement post-processing are checked. They confirm the a priori error estimates developed before. The results concerning the duct’s second problem are presented in Figs. 9–12. In those graphics, there are two performance indices being plotted, namely e∗ and eh . The former
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
471
y x
Fig. 3. Non-regular mesh, 26 × 26 quadrilateral elements.
0.2 0.2 P=Po 0.2
0.6
1.2
0.6
Fig. 4. Expansion chamber geometry, boundary conditions and typical mesh.
represents the error committed by approximating the velocity field using the post-processed solution. The second one aims to stress the improvement on the results by comparing the recovered results with the direct differentiation of the finite element solution ph . It is important to highlight that this was still observed for high values of the wave number, as k = 20 (f = 1082 Hz): Figs. 13 and 14. As it was stated previously, there is a pre-asymptotic behavior for high wave numbers which can be clearly seen in Fig. 13. Some numerical tests were carried out using non-uniform meshes leading to irregular elements, as shown in Fig. 3. The obtained post-processed velocity fields did not present high sensitivity to this nonregularity, although there was an expectation of non-optimal convergence rates due to a potential loss of superconvergence. Besides, the wave propagation becomes no longer parallel to the grid, which tends to increase the pollution effect.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -3.2 -3.4 -3.6 1.0
-3.8
log||e||
1.0
-4 1.1
-4.2 -4.4 2.0
-4.6 -4.8 0.7
fi=4 (Q4) ri=1(Q4) fi=3 (T3) ri=1(T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 5. 1D duct: superconvergence study at f = 54 Hz, linear elements.
fi=4 (Q4) ri=1 (Q4) fi=3 (T3) ri=1 (T3)
-2.8
-3
-3.2
1.2
log||e||
472
-3.4 1.2 1.6
-3.6 1.9
-3.8 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
-log (h)
Fig. 6. 1D duct: superconvergence study at f = 649 Hz, linear elements.
1.5
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-5 2.0
log||e||
-5.5 2.0
-6
-6.5 3.0
fi=9 (Q9) ri=4 (Q9) fi=6 (T6) ri=3 (T6)
-7
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 7. 1D duct: superconvergence study at f = 54 Hz, quadratic elements.
-4 -4.2 2.0
-4.4 -4.6 2.0
log||e||
-4.8 -5 -5.2
3.4
-5.4 fi=9 (Q9) ri=4 (Q9) fi=6 (T6) ri=3 (T6)
-5.6 -5.8 -6 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
-log (h)
Fig. 8. 1D duct: superconvergence study at f = 649 Hz quadratic elements.
1.5
473
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -4 -4.2 -4.4 1.0
-4.6
log (e)
1.0
-4.8 -5 -5.2 eh (Q4) e* (Q4) eh (T3) e* (T3)
-5.4
1.8 1.2
-5.6 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 9. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, linear elements.
-3.3 -3.4 -3.5 -3.6 -3.7
log (e)
474
1.3
-3.8 1.6 1.3
-3.9 -4 -4.1 -4.2
0.7
1.6
eh (Q4) e* (Q4) eh (T3) e* (T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 10. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, linear elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-5.5
log (e)
-6
2.0
-6.5
-7 2.3
-7.5
-8 0.7
eh (Q9) e* (Q9) eh (T6) e* (T6) 0.8
3.2
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 11. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, quadratic elements.
-4.4 -4.6 -4.8 -5
2.0
log (e)
-5.2 -5.4 -5.6 3.0
-5.8 -6
eh (Q9) e* (Q9) eh (T6)
-6.2
e* (T6)
-6.4 0.7
0.8
0.9
3.8
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 12. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, quadratic elements.
475
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2
log (e)
-2.5
-3 3.0
-3.5
-4 0.7
3.0
eh (Q4) e* (Q4) eh (T3) e* (T3) 0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 13. 1D duct: convergence study for macroelement post-processing at f = 1082 Hz, linear elements.
-4.2 -4.4 -4.6 -4.8
log (e)
476
2.0
-5 -5.2 3.0
-5.4 -5.6 -5.8 0.7
3.7
eh (Q9) e* (Q9) eh (T6) e* (T6) 0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 14. 1D duct: convergence study for macroelement post-processing at f = 1082 Hz, quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
477
Those tests were performed on 10 × 10, 14 × 14, 18 × 18, 22 × 22 and 26 × 26 meshes. To confirm once more the robustness of the proposed method, several tests were accomplished for the two proposed examples, but only results related to the duct problem involving a unitary impedance condition in the left end are graphically depicted, as similar conclusions could be obtained from all the tested examples. The 26 × 26 mesh is illustrated in Fig. 3, emphasizing the distortion of the elements. Figs. 19–22 report the performance of the post-processing for non-regular meshes. The error and its convergence rate are comparable to the ones obtained by adopting regular meshes, which are depicted in Figs. 15–18. The quality of an approximate solution can also be assessed by means of its derivatives. In Figs. 23 and 24, the divergence of the recovered velocity is compared to the exact solution and to the error of the divergence of the gradient of the finite element approximation as well. The better performance, for the analyzed frequency range, of the proposed methodology, when quadrilaterals were used, is observed once more. As it has been mentioned, there are some drawbacks in applying the finite element method to the Helmholtz equation. This can be illustrated by analyzing the dispersive characteristics of the obtained numerical wave. A valuable tool for it, often utilized by the modal analysis community, is the frequency response function (FRF). In general terms, it relates the acoustical response, here taken as e∗ and eh , to the excitation frequencies.
-3.2 -3.4 -3.6 1.0
-3.8
log (e)
1.0
-4 -4.2
1.7
-4.4
eh (Q4) e* (Q4) eh (T3)
-4.6 -4.8 -5 0.7
2.0
e* (T3)
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 15. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, regular meshes: linear elements.
478
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
-2.8
-3
log (e)
-3.2 0.9
-3.4 1.2
-3.6
-3.8
-4 0.7
1.96
eh (Q4) e* (Q4) eh (T3) e* (T3)
0.8
1.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 16. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, regular meshes: linear elements.
-5 2.0
log (e)
-5.5
-6 2.4
eh (Q9) e* (Q9) eh (T6) e* (T6)
-6.5
3.2
-7 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 17. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, regular meshes: quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
479
-4 -4.2 -4.4 2.0
log (e)
-4.6 -4.8 -5 -5.2
2.5
-5.4
eh (Q9) e* (Q9) eh (T6)
-5.6
e* (T6)
-5.8 -6 0.7
3.4
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
-log (h)
Fig. 18. 1D duct: convergence study for macroelement post-processing at f = 6494 Hz, regular meshes: quadratic elements.
-2.8 eh (Q4) e* (Q4)
-3
1.0
log (e)
-3.2
-3.4
-3.6 1.97 -3.8
-4
-4.2 0.8
1
1.2
-log (h)
Fig. 19. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, non-regular meshes: linear elements.
480
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.4 -2.5
eh (Q4) e* (Q4)
-2.6 1.2
log (e)
-2.7 1.8 -2.8 -2.9 -3 -3.1 -3.2 0.8
1
1.2
-log (h)
Fig. 20. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, non-regular meshes: linear elements.
-4 eh (Q9) e* (Q9) 2.0
-4.5
log (e)
-5
-5.5
3.1
-6
-6.5 0.8
1
1.2
-log (h)
Fig. 21. 1D duct: convergence study for macroelement post-processing at f = 54 Hz, non-regular meshes: quadratic elements.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
481
-3.2 eh (Q4) e* (Q4)
-3.4 -3.6
1.97
-3.8
log (e)
-4 -4.2
3.4
-4.4 -4.6 -4.8 -5 -5.2 0.8
1
1.2
-log (h)
Fig. 22. 1D duct: convergence study for macroelement post-processing at f = 649 Hz, non-regular meshes: quadratic elements.
-2
log (div e)
-2.5
-3 div eh (Q4) div e* (Q4) div eh (T3) div e* (T3)
1.5
-3.5 1
1.05
1.1
1.15
1.2
1.25
1.3
1.35
1.4
1.45
-log (h)
Fig. 23. 1D duct: convergence study for macroelement post-processed divergence at f = 271 Hz, linear elements.
482
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.1 -2.2 -2.3
log (div e)
-2.4 -2.5 -2.6 1.2
-2.7 div eh (Q4) div e* (Q4) div eh (T3) div e* (T3)
-2.8 -2.9 -3 1
1.05
1.1
1.15
1.2
1.25
1.3
1.35
1.4
1.45
-log (h)
Fig. 24. 1D duct: convergence study for macroelement post-processed divergence at f = 271 Hz, linear elements.
Those FRFs corresponding to meshes of 10 and 26 linear elements in the x direction are presented in Figs. 25–28. The following conclusions could be drawn. • For low wave numbers, e∗ and eh present a good and similar behavior, what was already expected from the analysis developed in the previous section. • For high wave numbers, the post-processing presents a deterioration in its performance, but it still represents a better approximation of the acoustic velocity than the gradient of the finite element solution. • All FRFs manifest the same tendency of growth deviation between approximations and exact value when the wave number tends to magnify. This phenomenon is in accordance with the increasing influence of the pollution mentioned before. 4.2. Two-dimensional cavity Once again, two situations are treated involving the same domain form: a two-dimensional cavity. In the first problem, the cavity occupies the region (0 x 10, 0 y 10) and normal velocity is set to be zero at boundaries, therefore the acoustic mode corresponding to natural frequency of 49.02 Hz, is given by the closed form x y x y 1
v(x, y) = sin cos ex + cos sin ey , (37) i0 w 5 5 5 5 5 5 where ex and ey stand for the ordinary vectorial basis in the plane.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 100 eh (T3) eh (Q4) e* (T3) e* (Q4)
90 80 70
error%
60 50 40 30 20 10 0 0
2
4
6
8
10 k
12
14
16
18
20
Fig. 25. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.1 m.
100 eh (T3) eh (Q4) e* (T3) e* (Q4)
90 80 70
error%
60 50 40 30 20 10 0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 26. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.38 m.
483
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
45 eh (T6) eh (Q9) e* (T6) e* (Q9)
40 35
error%
30 25 20 15 10 5 0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 27. FRFs relating the global error eh and e∗ to the wave number k, quadratic elements: h = 0.1 m.
3
eh (T6) eh (Q9) e* (T6) e* (Q9)
2.5
2
error%
484
1.5
1
0.5
0 0
2
4
6
8
10
12
14
16
18
20
k
Fig. 28. FRFs relating the global error eh and e∗ to the wave number k, linear elements: h = 0.38 m.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
485
The same situation was analyzed in [19] where a finite element error estimator for the acoustic problem was proposed. The potentiality of using the recovery technique proposed here as the basis of an error estimator is now briefly addressed. This estimator is simply eh and its performance is assessed by comparing it to eh∗ . The indices will be referred to as estimated error and exact error and their distribution along the domain is presented in Figs. 29 and 30. In order to highlight the good quality of the results, the difference between the exact and estimated error is plotted in Fig. 31. The second problem consists in a unity cavity [7] with D comprising all edges of the square and N = 0. The imposed Dirichlet boundary condition is given as p(x, y) = exp[ik(x cos + y sin )],
(x, y) ∈ D .
(38)
Therefore the exact solution to this problem is a plane wave propagating at angle to the cartesian x-axis, and is given by p(x, y) = exp[ik(x cos + y sin )],
(x, y) ∈
(39)
and v(x, y) = p(x, y)ik exp[− sin ex + cos ey ],
(x, y) ∈ .
(40)
The solution to this problem was sought on uniform meshes of bilinear finite elements. Convergence rates are shown in Fig. 32 for a single =45◦ value and frequencies of f =649 and 1082 Hz for several meshes.
Fig. 29. 2D cavity: exact error (eh ) distribution.
486
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
Fig. 30. 2D cavity: estimated error (eh∗ ).
Fig. 31. 2D cavity: difference between exact and estimated error distributions.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
487
-2
-2.5
log (e)
-3 1.5 -3.5
-4 eh (f=649 Hz) e* (f=649 Hz) eh (f=1082 Hz)
-4.5
3.5
e* (f=1082 Hz) -5 1
1.2
1.4
1.6
1.8
2
-log (h)
Fig. 32. 2D cavity: convergence study for = 45◦ .
Two performance indices are plotted, namely e∗ and eh . Please note that convergence rates were visibly better for the post-processing and the approximate solution didn’t present any earnings, independent of the wave number. A similar behavior is presented in Fig. 33, in which and k were varied and the mesh was fastened with 40 × 40 elements. It is worth to observe the deterioration on the post-processing performance with the increasing of k due to the pollution effect. That happens for all searched wave propagation direction , although for high wave numbers as k = 20, the proposed recovery still presents good solutions for = 45◦ .
4.3. Expansion chamber This is a more challenging situation, as it deals with a two dimensional real like problem. It consists in an expansion chamber with a perforated outlet pipe, schematically shown in Fig. 4. The excitation stems from a noise source with p(0, y) = 2.0 at the left-hand side. At the right-hand side, the interface between the interior fluid and the exterior is modeled by a constant impedance. A comparison of the imaginary and real components in the middle section’s chamber between vh and v∗ for k =10 (f =541 Hz) is presented in Figs. 34–37. From them it is possible to figure out the improvement on using the proposed post-processing for obtaining the acoustic velocity field. Those results show that the recovered velocity was able to capture significative variations due to abrupt changes in the domain geometry.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 -2.5 -3 -3.5
eh (theta = 22.5) e* (theta = 22.5) eh (theta = 45) e* (theta = 45) eh (theta = 75) e* (theta = 75)
log (e)
-4 -4.5 -5 -5.5 -6 -6.5 -7 1
2
3
4
5
6
7
8
9
10
K
Fig. 33. 2D cavity: post-processing’s numerical performance for several ’s.
x 10-3
1.5
v vh v*
1
0.5 real velocity (ur)
488
0
-0.5
-1
-1.5 0
0.5
1
1.5
2
x
Fig. 34. Expansion chamber: horizontal real velocity component along the center line.
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492 0.01 v vh v*
0.008
imaginary velocity (ui)
0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 0
0.5
1
1.5
2
x
Fig. 35. Expansion chamber: horizontal imaginary velocity component along the center line.
x 10-4 6 v vh v*
4 2
real velocity (vr)
0 -2 -4 -6 -8 -10
0
0.5
1
1.5
2
x
Fig. 36. Expansion chamber: vertical real velocity component along the center line.
489
490
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
2
x 10-3
1
imaginary velocity (vi)
0 -1 -2 -3 -4 v vh v*
-5 -6 -7 0
0.5
1
x
1.5
2
Fig. 37. Expansion chamber: vertical imaginary velocity component along the center line.
5. Concluding remarks The macroelement post-processing recovery technique introduced in [15] was extended here for the steady-state acoustic problem governed by the Helmholtz equation. The obtained analytical and numerical results have demonstrated the ability of this method for handling the inherent numerical difficulties linked to the finite element approximation of this equation. The proposed method achieves good accuracy for the velocity acoustic field at low computational cost. The convergence of the method is proved leading to a priori error estimates, which are confirmed by several numerical examples usually explored in the literature. The computation implementation is very simple, only requiring a standard finite element code and is amenable for parallel computing. The velocity post-processing has been successfully used in noise source estimation [12], enhancing the numerical performance of the iterative algorithm used to solve the resulting inverse problem, and also as a local error estimator [14] to be explored in adaptive finite element analysis. At this point, it is important to mention that the post-processing could be regarded in a different perspective of the one adopted here. It could be undertaken not based on a standard finite element framework. This means that either vh∗ could be approximated by using different basis functions, such as harmonic functions [6,20], or by using as input data (the pressure field) not a finite element solution of
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
491
Helmholtz equation. In this last case, a typical application would be obtaining a smooth representation of the velocity field from a cloud of pressure experimental values.
Acknowledgements We gratefully acknowledge the partial support from the Conselho Nacional de Desenvolvimento e Pesquisa (CNPQ).
References [1] Ph. Bouilarrd, F. Ihlenburg, Error estimation and adaptivity for the finite element method in acoustics: 2D and 3D applications, Comput. Methods Appl. Mech. Eng. 176 (1999) 147–163. [2] F. Ihlenburg, Finite element analysis of acoustic scattering, Appl. Math. Sci. 132 (1998). [3] A. Deraemaeker, I. Babuska, P. Bouillard, Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Meth. Eng. 46 (1999) 417–499. [4] I. Harari, T.J.R. Hughes, Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains, Comput. Methods Appl. Mech. Eng. 97 (1992) 103–124. [5] R. Baus˜ys, N.E. Winberg, Adaptive finite element strategy for acoustics problems, J. Sound Vibr. 226 (5) (1999) 905– 922. [6] C. Farhat, I. Harari, L.P. Franca, The discontinuous enrichment method, Comput. Meth. Appl. Mech. 190 (2001) 6455– 6479. [7] A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Meth. Eng. 49 (2000) 399–419. [8] S. Suleau, Ph. Bouillard, One-dimensional dispersion analysis for the element-free Galerkin method for the Helmholtz equation, Int. J. Numer. Meth. Eng. 47 (2000) 1169–1188. [9] L.L. Thompson, P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Meth. Eng. 38 (1995) 371–397. [10] K.S. Peat, K.L. Raithi, A finite element analysis of the convected acoustic wave motion in dissipative silencers, J. Sound Vibr. 184 (3) (1995) 529–545. [11] T. Lin, D.L. Russel, A superconvergent method for approximating bending moment of elastic beams with hysteris damping, Appl. Numer. Math. 38 (2001) 145–165. [12] R.H.P. Romero, F.A. Rochinha, A post-processing technique applied to error estimation of FEM approximations of Helmholtz equation, Proceedings of the Seventh US National Congress on Computational Mechanics, vol. 107, New Mexico, USA, 2003. [13] R.H.P. Romero, A post-processing scheme to the Helmholtz equation and applications, Doctor in Science Thesis, Federal University of Rio de Janeiro, Brazil, 2003 (in Portuguese). [14] R.H.P. Romero, F.A. Rochinha, Estimation of noise sources using an enhanced finite element formulation, Proceedings of the Seventh US National Congress on Computational Mechanics, vol. 612, New Mexico, USA, 2003. [15] A.F.D. Loula, F.A. Rochinha, M.A. Murad, Higher-order gradient post-processings for second-order elliptic problems, Comput. Methods Appl. Mech. Eng. 128 (1995) 361–381. [16] A.F.D. Loula, M.A. Murad, F.A. Rochinha, 3D gradient post-processings schemes for second order elliptic problems, in: Proceedings of the IV World Conference on Computational Mechanics, Buenos Aires, Argentina, 1998. [17] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Numer. Meth. Eng. 33 (1992) 1331–1364.
492
R.H. Paulino, F.A. Rochinha / Finite Elements in Analysis and Design 41 (2005) 459 – 492
[18] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 2: error estimates and adaptivity, Int. J. Numer. Meth. Eng. 33 (1992) 1365–1382. [19] Ph. Bouilarrd, J.F. Allard, G. Warzée, Superconvergent patch recovery technique for the finite element method in acoustics, Comm. Num. Methods Eng. 12 (1996) 581–594. [20] C. Farhat, I. Harari, U. Hetmaniuk, A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime, Comput. Meth. Appl. Mech. 192 (2003) 1389–1419.