An improved axisymmetric convected boundary element method formulation with uniform ﬂow

– This paper presents an improved form of the convected Boundary Element Method (BEM) for axisymmetric problems in a subsonic uniform ﬂow. The proposed formulation based on the axisymmetric Convected Helmholtz Equation (CHE) and its fundamental solution that describes the sound radiation from a monopole source. The variables in the new axisymmetric boundary integral formulation can be expressed explicitly in terms of the acoustic pressure and its particular normal derivative. Also, the constant coeﬃcients are expressed only in terms of the axisymmetric convected Green’s function and its convected normal derivative. The particular and convected derivatives reduce the ﬂow eﬀects of the normal and the ﬂow direction derivatives incorporated in the conventional convected boundary integral formulas. The advanced form of the axisymmetric boundary integral representation with ﬂow is a similar form of the axisymmetric boundary element method without ﬂow. Precisely, the two new operators signiﬁcantly reduce the computational burden of the classical BEM and then becomes the CPU time of BEM without ﬂow. The formula is veriﬁed comparing to both analytical and Finite Element Methods (FEM) of an axisymmetric inﬁnite rigid duct in a subsonic uniform ﬂow.


Introduction
The numerical studies of some acoustic problems with fluid flow require only the special methods, such as the computational techniques of the Finite Element Method (FEM) and boundary element method (BEM) for solving any axisymmetric medium in the science and engineering. The finite element method is currently the most widely used numerical method for solving the interior medium with mean flow and that the boundary element method is used for solving the exterior medium with a constant flow [1].
The finite element method has been implemented for solving the Helmholtz equation for the acoustic wave diffraction problems without flow [2]. It has been implemented with a displacement formulation for the modelling and simulation of the fluid-structure interaction problems in which the coupling conditions are formulated [3]. A numerical technique of the FEM for the modeling of the acoustic pressure inside an axisymmetric lined duct in a Corresponding author: bassem.barhoumi@enit.rnu.tn the presence of a constant flow was developed and presented in Reference [4]. This method is based on the resolution of the convected Helmholtz equation in the modal representation which gives good results compared with the analytical ones. The FEM applied to determine the time-harmonic acoustic field in a 2D infinite waveguide with walls covered of an absorbing material, and in the presence of a mean flow assumed uniform far from the source [5]. The boundary conditions of the rigid and lined ducts are represented by transparent boundary conditions suitable for use in a finite element scheme which based on the Dirichlet-to-Neumann (DtN) boundary condition [5]. For particular acoustic domains [6], has determined the acoustic propagation in a transformed domain by the Prandtl-Glauert transformation.
It is known that the finite element method has its limitations in modeling infinite domains. The unbounded medium contains the non-reflection condition which is represented by the Sommerfeld radiation condition in a medium at rest [7] or in a subsonic uniform flow [6]. Using the boundary element method (BEM), which requires a discretization of only the generator of the acoustic domain and that this Sommerfeld condition is automatically Article published by EDP Sciences fulfilled by the fundamental solution. This makes the BEM superior to the FEM in many cases, especially for the unbounded problems, such as the wave propagation in the infinite lined ducts.
The boundary element method can be formulated in the transformed acoustic medium [8] and in the original acoustic medium [9]. In order to take advantage of both the FEM and the BEM, a coupled FEM/BEM approach has been proposed [10]. A direct collocation boundary element formulation for the calculation of acoustic propagation in a subsonic uniform mean flow has been presented in Reference [11] in which the model is directed towards the calculation of the propagation inside and radiation from axisymmetric lined ducts. An advanced boundary element/fast Fourier transform (BE/FFT) methodology for solving axisymmetric acoustic wave scattering and radiation problems with non-axisymmetric boundary conditions has been developed in Reference [12].
However, the main disadvantages of using boundary integral formulations in an original acoustic medium are the same as in the non-flow case [13]. In addition, the convection effects in the presence of flow substantially increase these difficulties and makes significantly more complicated the conventional integral representations [10,11]. These axisymmetric formulations are derived from the three-dimensional boundary integral formulation developed in Reference [9] in which the convected terms due to the Green's function and its derivatives such as the normal derivative and the derivative in the flow direction.
Also, the authors in Reference [14] have been used a reformulation of the three-dimensional convected boundary element to solve the axisymmetric acoustic problems. The reformulation based on the non-standard normal derivative similar to the non-flow normal derivative of the Green's kernel. The non-standard derivative leads to a global form of the integral formula contains another supplementary convective terms in which make difficulty the numerical implementation of the boundary element method. The acoustic variable of the classical and reformulation axisymmetric boundary integral representations in References [10,11,14] in terms of the acoustic pressure as well as its normal and tangential derivatives.
The aim of this paper is to contribute to reduce the complexity of convected boundary integral formulas for axisymmetric domains in a subsonic uniform flow. Thus, a new reformulation of the boundary element method for the axisymmetric mediums (BEMA) which based on the axisymmetric convected Helmholtz equation and its fundamental solution is presented. This new axisymmetric boundary element method requires only two operators which substantially reduce the flow effects of the normal and the flow direction derivatives of the axisymmetric convected Green's function incorporated in the classical convected boundary integral formulations. Through the use of these new convected terms, the new axisymmetric BEM formula has a similar form of the axisymmetric BEM without flow. The acoustic variables in the new axisymmetric boundary integral equation can be expressed explicitly in terms of the acoustic pressure and its particular normal derivative.
The numerical implementation of BEMA is developed for the acoustic propagation in an axisymmetric rigid duct and validated by comparison with FEMA method and the analytical solution in which the CPU times are presented.

An improved Axisymmetric BEM Formulation
Consider an acoustic medium Ω with a compressible fluid, inviscid, isentropic, which is characterized by the Mach flow vector M ∞ , of density ρ ∞ and speed of sound c ∞ in the axisymmetric physical space (o, r, z). In Figure 1, ∂Ω = Γ ∞ ∪ Γ is the generator limiting the axisymmetric medium Ω, such as Γ ∞ and Γ are the exterior and interior generators of the outgoing normal vector n q (n rq , n zq ) and the tangential vector t q (t rq , t zq ) at the point q, respectively. The generator ∂Ω not perturbed the streamlines and unchanged the direction of the flow.

Axisymmetric convected helmholtz equation
A combination of the Axisymmetric Linearized Euler Equations can be produced the Axisymmetric Convected Helmholtz Equation in a harmonic regime (exp(+iωt)) in which t is the time, ω is the angular frequency and i is the imaginary unit. The wave equation wearing on axisymmetric acoustic pressure p is decomposed into a non-flow part and a convection part [10] where Δ and ∇ are the axisymmetric Laplace and Gradient operators. k = ω/c ∞ is the original wavenumber. The fundamental solution of Equation (1) is the axisymmetric convected Green's function G k 0 = 2π 0 G k dβ, which is the integration of the three-dimensional convected Green function G k over the azimuthal angle β. The Green's kernel G k (m, q) is given for the points m and q by where m(r m , z m ) and q(r q , z q ) are the axisymmetric source-observer points, α is the Prandtl-Glauert factor, k * = k/α is the convected wavenumber, and that M * ∞ = M ∞ /α is the convected Mach vector. The distance R * β is the convected radius which depends on the azimuthal order β, it's given for the points m and q by in which r * (m, q) = mq 2 + (mq.M * ∞ ) 2 is the axisymmetric convected radius and mq is the physical distance between source-observer. Also, the axisymmetric convected Green's function satisfies the following inhomogeneous convected Helmholtz equation in which (·) q is the operator at point q and δ is the Dirac delta function such as

Axisymmetric convected Helmholtz integral equation
To obtain a boundary integral equation formulation for the acoustic problem, the Helmholtz Equation (1) is multiplied by the axisymmetric convected Green's function G k 0 (m, q) and its fundamental Equation (4) by p(q). Subtracting these two equations and integrating over the axisymmetric exterior medium Ω yields Using the axisymmetric second Green's formula, the axisymmetric Gradient property and the axisymmetric divergence theorem (Appendix A), the axisymmetric acoustic pressure Equation (5) takes the following form where ∇ q = σ n n q + σ t t q is the axisymmetric Gradient operator in which the couple (σ n = ∂/∂ nq , σ t = ∂/∂ tq ) designates the normal and tangential derivatives at point q, respectively. M ∞n = M ∞ .n q is the normal Mach number.
In order to simplify the boundary integral Equation (6), the normal derivative ∂/∂n q and the flow direction derivative (M ∞ .∇ q ) can be converted into a particular normal derivative d/dn q , which is given by the following relation The particular normal derivative of the axisymmetric convected Green function dG k 0 /dn q = 2π 0 dG k /dn q dβ is the integration of the particular normal derivative of the three-dimensional Green kernel dG k /dn q over the azimuthal angle β. Using the normal and the flow direction derivatives of the three-dimensional Green's function in References [10,11,14], one obtains the following derivative function where r nβ = mq axi .n q is the normal distance and that the vector mq axi of coordinate (r q − r m cos(β), z q − z m ).
The particular normal derivative of the convected Green's function Equation (8) is decomposed into two terms. The first term generalizes the expression of the normal derivative in a non-flow case and that the second term explains the convective influence. Substituting the particular normal derivatives Equations (7) and (8) in the integral formula Equation (6) yields The axisymmetric boundary integral formulation Equation (9) is decomposed into two integrals in which the acoustic variables are expressed only in terms of the acoustic field as well as its particular normal derivative which interpreting as a condition applied to the interior generator Γ . Thus, the axisymmetric BEM Equation (9) is similar to the boundary integral equation solution for solving the axisymmetric potential flows [15] and the nonflow case (the right first term in the Equation (6)) [16]. From the particular normal derivative operator Equation (7), we obtain the convected normal derivative operator

313-page 3
Substituting the particular normal derivative of the threedimensional Green's function Equation (8) into the convected normal derivative operator Equation (10), one obtains The integration over the circular generator Γ ∞ given by Equation (9), represents contributions waves which come in the interior region due to the source is very far from observer point; q → ∞. Thus, this integral is null. According to the three-dimensional Green's function at infinity [17], the axisymmetric Green's kernel G k 0 is proportional to 1/ R * β . Also, taking into account the term r nβ /R * β is equivalent to the tangential Prandtl-Glauert factor This leads to two convected radiation conditions which verified by the axisymmetric acoustic pressure and its Green's kernel lim r→∞ √ r dp dr (r)+i k where M ∞t = M ∞ .t q is the tangential Mach number. The first condition Equation (12) indicates the non-reflection of waves, and the second condition Equation (13) generalizes the classical Sommerfeld radiation condition in a uniform flow subsonic [7,8]. Substituting the non-reflection conditions into the boundary Equation (9), we obtain (14) in which the signs (+, −) designates the exterior and interior acoustic domains. When the point m is taken on the generator at point q, the Green's kernel and its convected normal derivative Equations (2) and (11) contain singular and regular parts. Thus, Equation (14) requires to study of behavior of the kernels near the singularity.

Singular integrals
A classical procedure can be used to isolate the singularity problems of the axisymmetric Green's function and its convected derivative. In fact, the singularities can be isolated by decomposing the axisymmetric convected Green's function G k 0 into a singular static part G 0 0 (k = 0) and regular part g k 0 depending on the wavenumber. Then, the axisymmetric Green's function takes the following form Thus, the convected normal derivative is given by Because the difference between the singular and regular behavior, the singular functions in Equations (15) and (16) are evaluated analytically and that the regular parts are evaluated numerically with standard Gauss-Legendre quadrature of high order. In fact, the static axisymmetric Green's function is given by [10,11] G 0 0 (m, q) = Using the normal and the flow direction derivatives of the static Green's kernel in Reference [10], one obtains the following static convected derivative where K(τ ) and E(τ ) are the complete elliptic integrals of the first and second kinds, the parameters τ and τ * are given by According to the complete elliptic integrals at r → 0 [17] , the static Green's function exhibits a logarithmic behavior near the singularity, as follows Thus, the static convected normal derivative Equation (18) near the singularity is characterized by Since r n0 /r * = O(r * ), then the singularity of the static convected function d k G k 0 /dn q is thus only apparent. Following the exclusion procedure [18,19] to isolate the singularity problems in the boundary formula Equation (14), which can be rewritten as where s ε is the disk of center m and of radius ε = r bounded by the circle Γ + ε ∪ Γ − ε and Γ = lim ε→0 (Γ ε ∪ Γ − ε ) for the exterior domain and that Γ = lim ε→0 (Γ ε ∪ Γ + ε ) for the interior domain (see Fig. 2).
Substituting Equations (20) and (21) in Equation (22), we get a new original form of the integral boundary representation with uniform flow wearing on the axisymmetric acoustic pressure at any point m where c ± (m) is the free term at any point m, which is given by in which θ ± (m) is the convected angle. When the point m in the exterior and interior domains, the convected angle θ ± (m) = 0 and that the coefficient c ± (m) = 1. If m is a regular point of the generator, θ ± (m) = π and that c ± (m) = 0.5.

Numerical implementation
The numerical implementation of the axisymmetric boundary integral Equation (23) for a body of arbitrary shape is obtained by discretizing the generator Γ of the body with N one-dimensional quadratic isoparametric elements, according to the following convergence criterion where n e is the number of finite elements for a wavelength λ and h e is the diameter of the finite element. Then, by using the quadratic shape functions, the acoustic pressure p and its normal derivative σ n of any point q on the generator Γ are assumed to be given in terms of the nodal coordinates. Thus, the tangential derivative in the local coordinates system along the generator Γ at point q is given by [11,20] with −1 ≤ ξ ≤ 1 is the local coordinate and J i is the Jacobien of reference-local transformation. Using the collocation method so that the point m coincides with the nodes q i of the discretized generator Γ , Equation (14) can be written in complex matrix form with Robin condition where p, p c and σ n are the nodal vectors containing the values of the axisymmetric acoustic pressure p, its Dirichlet acoustic pressure p c and its Neumann normal derivative σ n at each node, respectively, whereas [A] and [B n ] are the axisymmetric acoustic matrices with the elements of matrix [A] contain the convected angles and elementary integrals related to the axisymmetric convected Green's function G k 0 and its convected normal derivative d k G k 0 /dn q , while the elements of matrix [B n ] is composed only by the elementary integrals containing the integrands of axisymmetric Green's function and its convected derivative. [B c ] is the acoustic matrix containing the integrands G k 0 and d k G k 0 /dn q .

Axisymmetric finite element method
The axisymmetric Finite Element Method is based on the weak variational formulation, which can be obtained by multiplying the axisymmetric wave Equation (1) by a test function p * (the axisymmetric acoustic pressure p * satisfies the axisymmetric convected Helmholtz Eq. (1)). The axial and radial wavenumbers satisfies the convergence criterion of the discretization Equation (25) in which the number of triangles is 1652 elements (see Fig. 4).
The evaluation of the acoustic pressure at point m in an interior axisymmetric rigid duct Ω given by Equation (23) is numerically where Γ = Γ − ∪Γ R ∪Γ + is the generator of this duct, and taking into account that the free term at the corner points m = Γ − ∩ Γ R and m = Γ R ∩ Γ + is c − (m) = 3/4.

Numerical results
We apply BEMA, FEMA techniques and the analytical solution (S.Anal) associated to the axisymmetric acoustic pressure which propagated in the increased zdirection of axisymmetric excitation (p + n /A n ). Figures 5  and 6 give the cartographies of the real parts Re(P + n /A n ) for different values of the Mach vector by S.BEMA, S.FEMA and S.Anal methods. Figures 5 and 6 are presented the real parts of the axisymmetric acoustic pressure Re(P + n /A n ) by FEMA and BEMA which are in very good concord with analytical solution of relative error given by the L 2 norm. For FEMA, the relative error is less than 0.5% and for the axisymmetric boundary integral method, the error is less than 0.8%.
We observe that for the initial mode (0), the axisymmetric acoustic pressure is uniform along the radial and axial waves. Also, when the subsonic uniform flow is increased, the wavelength increases, which is explained by the offset of the wave in the positive z-direction (see Fig. 5).
In the other mode, the axisymmetric acoustic pressure along the radial and axial axis is non-uniform in the presence and absence of the flow. Indeed, in the upstream region, the acoustic field is radially amplified to be greater than that obtained in the downstream region in which we observe the Sommerfeld radiation condition Equation (12), and we notice also that the uniform flow is proportional to the wavelength (see Fig. 6).
In addition, the CPU times of the Numerical boundary element method, the finite element method and that the classical boundary element method in Reference [14] with mean flow are obtained by a machine 2.93 GHz using MATLAB-ACOUSTIC codes, which are given by the following Table 1.
A comparison between the CPU times of FE and BE methods in Table 1, show that the axisymmetric Convected boundary element formulation Equation (23) reduces 80% of the computation time of Finite element method Equation (32).
Compared to the classical axisymmetric boundary integral formulae, the proposed method reduces a CPU time equal to 60% of the classical axisymmetric BEM. Also, the relative error for the classical axisymmetric boundary integral equation is less than 2%.
We observe that when the number of the discretization increase in term of the numerical implementation of the axisymmetric boundary element method with mean flow, the gain time between BEM-FEM and BEM-Classical BEM increase.
For the BEM case, we observe that the numerical execution time for the non-flow case is the numerical execution time for the flow case. This is because the improved form of the axisymmetric boundary integral representation Equation (23) with the mean flow is similar to the boundary element method without flow.

Conclusion and prospects
In this work, an improved form of the axisymmetric Helmholtz integral equation and conventional finite element method for an axisymmetric medium in a subsonic uniform flow have been developed and presented. The axisymmetric convected boundary element method is derived from two new operators such as the convected and particular derivatives, which substantially reduces the complexity due to the presence of flow in the conventional formulations. The formulation is derived to be easy to implement as a numerical tool for computational codes of the axisymmetric acoustic medium in a mean flow.
Also, the improved boundary integral representation significantly reduces the computational burden of the convenctional axisymmetric finite element method. Precisely, for the rigid duct case, the axisymmetric BEM require three Matrices in terms of integrals over the generator and that axisymmetric FEM require six assembly matrices in terms of integrals over the medium and its generator. The proposed methodologies can be applied to simplified the modal boundary element method and its derivatives [11,23].
Finally, we think that the method could be extended for the analysis of physical phenomena in other noise pollution systems such as compressors, aircraft engines, and ventilation systems that contain different materials. It is hoped that this could have a significant impact in design optimization of flow duct systems in the industry sector.

Appendix A. Axisymmetric acoustic properties
The axisymmetric Green's identity formula as a function of the test function p * is given by Ω (Δ q p * (q) p (q) − Δ q p (q) p * (q)) dΩ q = ∂Ω p * (q) ∂p ∂n q (q) − p (q) ∂p * ∂n q (q) r q dΓ q , (A.1) 313-page 7  where ∂/∂ nq is the normal derivative and that Δ q is the axisymmetric Laplace operator in an acoustics problem r−z Δ q (·) = ∂ 2 ∂r 2 q (·) + 1 r q ∂ ∂r q (·) + ∂ 2 ∂z 2 q (·) , (A.2) Using the axisymmetric Green's identity formula Equation (A.1) for the axisymmetric Green's function 3) The axisymmetric Gradient property is given by And that the axisymmetric Divergence theorem is given by