- Email: [email protected]
Direct numerical simulation of aerodynamic sounds by a compressible CFD scheme with node-by-node finite elements
Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910 www.elsevier.com/locate/cma
Direct numerical simulation of aerodynamic sounds by a compressible CFD scheme with node-by-node finite elements J. Tsuchida a, T. Fujisawa b, G. Yagawa a
c,*
Toyota Motor Corporation, 1 Toyota-cho, Toyota City, Aichi Prefecture 471-8571, Japan b Prometech Software Inc., 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan c Tokyo University, 1-32-19, Akatsutsumi, Setagaya-ku, Tokyo 156-0044, Japan
Received 5 November 2004; received in revised form 5 April 2005; accepted 26 May 2005
Abstract A direct numerical simulation (DNS) of aerodynamic sounds is performed by using a compressible CFD scheme with node-by-node finite elements. Although the computational simulations of aerodynamic sounds generated from three-dimensional objects with complex boundary shapes are required in various industrial applications, such simulations involve a lot of difficulties regarding instability in numerical schemes and enormous computational costs. In this paper, the present authors have combined the C-CUP (CIP-Combined Unified Procedure, CIP: Constrained Interpolation Profile) method and the CIVA (Cubic Interpolation with Volume/Area Coordinates) method with the Free Mesh Method (FMM) in order to achieve efficient simulation of aerodynamic sounds using unstructured computational grids of finite elements on parallel computers. As a numerical example, the edge tone is successfully computed using the present method with parallel computing with Hitachi SR8000 supercomputer. 2005 Elsevier B.V. All rights reserved. Keywords: Computational fluid dynamics; Finite element method; Edge tone; Sounding mechanism
1. Introduction Aerodynamic sound analyses are usually performed by employing some approximation models in order to save computational costs. The approximation models often used can be categorized into the following *
Corresponding author. Tel./fax: +81 3 3327 5680. E-mail address: [email protected] (G. Yagawa).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.05.039
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1897
two kinds: One is the hybrid method [1], which combines incompressible fluid analysis and acoustic analogy. This method is the most practical approach for calculation of aerodynamic sounds from industrial point of view. Another is the acoustic/viscous splitting method [2,3], which computes average flows and a gap from average flows separately. The average flows are computed by incompressible fluid analysis and a gap from average flows is computed by solving disturbance equations. However, the progress of computers in recent years has made it possible to compute aerodynamic sound by solving the compressible Navier–Stokes equations directly without introducing any approximation models, i.e. DNS (Direct Numerical Simulation of aerodynamic sound). The DNS has a fundamental merit to capture the interference between average flows and sound waves, especially when the resonance occurs in the fluid field. Nevertheless, the computational costs required in the DNS is very large, so that the DNS reported so far was limited to two-dimensional objects with simple boundary shapes. The DNS of threedimensional object with complex boundary shapes is still difficult to solve due to enormous computational costs. In order to overcome this situation, this study aims at the DNS of aerodynamic sound generated from three-dimensional objects, using unstructured computational grids of finite elements on parallel computers. We employ unstructured computational grids to concentrate computational power on the critical parts and then solve the problem efficiently even if the objects have complex boundary shapes. Since the sound is a phenomenon of propagation of compressible waves, both generation and propagation of aerodynamic sound can be captured by compressible fluid analysis from theoretical point of view. However, conventional density-based numerical schemes for compressible flows cause so-called stiffness problem under low Mach number and numerical simulations confront serious instability. In order to overcome this problem, we employ the C-CUP (CIP-Combined Unified Procedure, CIP: Constrained Interpolation Profile) method [4], which can treat both incompressible and compressible flows simultaneously. Therefore, it can perform compressible flow analysis under low Mach number. The C-CUP method is one of the splitting methods, which compute the governing equations by separating advection and non-advection terms. The method originally based on the finite difference method, and the advection terms are solved by the CIP (Constrained Interpolation Profile) method [5,6]. However, in the present study, we employ unstructured computational grids of finite elements, so that the CIVA (Cubic Interpolation with Volume/Area Coordinates) method [7] is also introduced. The CIVA method is one of the third-order interpolation methods for triangular or tetrahedral elements, which is a generalized CIP method using unstructured computational grids. In this study, the advection terms are solved by the CIVA method, making use of the node-based data structure of the finite elements in the FMM (Free Mesh Method) [8–10], which makes it possible to solve advection terms efficiently on large-scale finite element meshes.
2. Compressible flow simulation scheme based on finite elements 2.1. Outline of the numerical method Parallel computing methods for finite element analysis can be classified into element-based and nodebased ones. In the present study, we use the free mesh method (FMM) [8–10], which is a node-based finite-element method. In the FMM, a node-based data structure is employed to memorize sparse matrices and connectivity data between nodes (Fig. 1). The FMM performs both mesh generation and solution of equations in a seamless manner. Therefore, it is an efficient method for adaptive analysis on parallel computers. Although a fixed computational grid is employed in the present study, node-based data structures of the FMM are used in order to perform parallel computing effectively on distributed memory systems [11].
1898
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 1. Free mesh method.
On the other hand, the C-CUP method [4] is one of the splitting methods, which solves governing equation separating into the advection and the non-advection parts. The following is a summarized procedure of the C-CUP method. 1. 2. 3. 4.
Solve the advection parts using CIP method. Solve the pressure equation. Compute the velocity and the density of the next time step. Renew the spatial derivatives of primitive physical values. The governing equations are oq þ u rq ¼ qr u; ot ou rp þ u ru ¼ þ Qu ; ot q op þ u rp ¼ qC 2S u; ot
ð1Þ ð2Þ ð3Þ
where u is the velocity, p is the pressure, q is the density and CS is the local sound speed. Qu expresses the force, which includes gravity, buoyancy, elastic stress tensor and surface tension. oq þ u rq ¼ 0; ot ou þ u ru ¼ 0; ot op þ u rp ¼ 0. ot
ð4Þ ð5Þ ð6Þ
These equations are solved by the CIP method or other similar methods. In this study, we solve the equations using a new method based on the CIVA method and the node-based data structure of the FMM, which is described later.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1899
The following are the non-advection parts: oq ¼ qr u; ot ou rp ¼ þ Qu ; ot q op ¼ qC 2S u. ot
ð7Þ ð8Þ ð9Þ
By discretizing the time derivatives on the LHS of Eqs. (7)–(9) and estimating the spatial derivatives on the RHS from the values at the next time step, we obtain qnþ1 q ¼ q r unþ1 ; Dt unþ1 u rpnþ1 ¼ þ Qu ; Dt q pnþ1 p nþ1 ¼ q C 2 ; S u Dt
ð10Þ ð11Þ ð12Þ
where n is the time step and Dt is the increment of time. The superscript * indicates the value after solving the advection parts. By taking divergence of Eq. (11) and substituting it into Eq. (12), we obtain the following pressure equation: nþ1 rp pnþ1 p ru r . ð13Þ þ ¼ 2 q Dt q C 2 S Dt The velocity un+1 and the density qn+1 are respectively computed by the following equations: nþ1 rp unþ1 ¼ u Dt Q u ; q pnþ1 p . qnþ1 ¼ q þ C 2 S
ð14Þ ð15Þ
The artificial viscosity is needed in order to stabilize the solution involving shock waves, when inconservative governing equation is used. Therefore, in this study, the following von Neumann–Richtmyer artificial viscosity [12] is employed: 8 < a q C c þ 1 lr u lr u if r u < 0; S 2 ð16Þ qv ¼ : 0 otherwise; where a is a constant parameter, c is the specific heat ratio and l is the representative length of each element. In this study, the average distance between the central node and the associated satellite nodes is used as the representative length of elements (see Fig. 1). In order to compute aerodynamic sound precisely, non-physical numerical reflection of waves at the boundaries must be suppressed. For this purpose, the NSCBC (Navier–Stokes Characteristic Boundary Condition) [13] is employed and non-reflecting boundary conditions are set at outflow boundaries. Since flow simulation under low Mach number includes extremely different time scales between the advection and non-advection parts, the multi-time-step integration technique [14] is employed. In the multitime-step integration, the advection and non-advection parts are solved by different time increments in order to perform calculation efficiently.
1900
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Under consideration of artificial viscosity, the NSCBC and the multi-time-step integration technique, the numerical algorithm we propose in this study is summarized as follows: 1. Determine the boundary conditions using the NSCBC. The values of the boundary nodes are fixed to the values computed by the NSCBC (only at the first step): ðun ; pn ; qn ; run ; rpn ; rqn Þ. 2. Solve the advection parts: ðu ; p ; q ; ru ; rp ; rq Þ. 3. Solve the pressure equation using iterative method such as Bi-CGSTAB: ðu ; p ; q ; ru ; rp ; rq Þ. 4. Renew the velocity: rp u ¼ u Dt ; q ðu ; p ; q ; ru ; rp ; rq Þ. 5. Renew the density: qnþ1 ¼ q þ
ðp p Þ ; C 2 S
ðu ; p ; qnþ1 ; ru ; rp ; rq Þ. 6. Add the physical and artificial viscosity to the velocity obtained above: unþ1 ¼ u þ DtQu ; Qu
¼
qnþ1
osij oqv ; oxj oxi
1
oui ouj þ sij ¼ l oxj oxi
2 ouk ; dij l 3 oxk
ðunþ1 ; p ; qnþ1 ; ru ; rp ; rq Þ. 7. Add the artificial viscosity to the pressure obtained above: pnþ1 ¼ p þ DtQp ; Qp ¼ ðc 1Þqv ru ; ðunþ1 ; pnþ1 ; qnþ1 ; ru ; rp ; rq Þ. 8. Renew the spatial derivatives: ðunþ1 ; pnþ1 ; qnþ1 ; runþ1 ; rpnþ1 ; rqnþ1 Þ. 9. Determine the boundary conditions using the NSCBC. The values of the boundary nodes are fixed to the values computed by the NSCBC: ðunþ1 ; pnþ1 ; qnþ1 ; runþ1 ; rpnþ1 ; rqnþ1 Þ.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1901
10. Repeat Steps 3–9 several times, if the multi-time-step integration technique is employed. 11. Return to Step 2. In the above procedure, l is the viscosity coefficient and dij is the Kronecker delta. 2.2. Calculation of advection terms Regarding the calculation of advection terms, we use a new method combining the CIVA method [7] and the node-based data structure of the FMM. The advection equation is expressed in the hyperbolic equation as follows: of of þu ¼ 0. ot ox
ð17Þ
If the time increment Dt is small enough, it can be assumed that the velocity is constant during the time step. Then, the solution of the advection equation is given by f ðx; t þ DtÞ ¼ f ðx uDt; tÞ.
ð18Þ
Eq. (18) means that if we need the values of x at t + Dt, we can obtain the values by computing x uDt at t. The value at x uDt can be obtained by interpolating the nodal values of the element which contains the point x uDt. In the CIP method, both the physical values and their spatial derivatives are memorized at each node, and the physical values between two nodes are interpolated by the cubic polynomial, so that advection terms can be calculated with high accuracy. On the other hand, since we use unstructured computational grids of the finite elements, the CIVA method is used instead of the original CIP method. The CIVA method approximates the third order interpolation within triangle or tetrahedral elements. Interpolation in the CIVA method is expressed as f~ ðr1 ; r2 ; r3 ; r4 Þ ¼
4 X i¼1
fi ri þ d
4 X
bjk fr2j rk þ cðr1 r2 r3 þ r2 r3 r4 þ r1 r2 r4 þ r1 r3 r4 Þg
ð19Þ
j;k¼1ðj6¼kÞ
with bjk ¼ fj fk þ ðxk xj Þfjx þ ðy k y j Þfjy þ ðzk zj Þfjz
0 6 c 6 0:5;
d ¼ 0 or 1
where f~ is the interpolated value, f1, f2, f3 and f4 are the nodal values and r1, r2, r3 and r4 are the volume coordinates. Here, the parameter c cannot be determined uniquely from the nodal values, because both physical values and its spatial derivative are zero at each node of tetrahedral elements and the terms of ri, rj, rk in the RHS of Eq. (19) is also zero. Then, we employ constant curvature conditions to determine c. In this case, the constant parameter c = 0.5 could be used for all elements. If the parameter c less than 0.5 is used, the numerical viscosity is added. The parameter d is used for a switching between the first and the third order interpolations. Generally speaking, numerical oscillation occurs where physical values change suddenly. In order to stabilize numerical oscillation, the interpolation is switched from the third to the first order where physical values change suddenly. A criterion of the switching is given by the following equations: ( nþ1 1 f min 6 f~ 6 fmax ; d¼ ð20Þ 0 otherwise;
1902
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 2. Advection calculation using the CIVA method.
where fmin ¼ minðf1n ; f2n ; f3n ; f4n Þ; fmax ¼ maxðf1n ; f2n ; f3n ; f4n Þ; nþ1 f~ is the third order interpolated value.
Although Eq. (20) has no mathematical background, several experiments show this stabilization method works well, if the computational grid is fine enough. Fig. 2 shows the basic concept of solving advection terms using the CIVA method and the node-based data structure of the FMM. In the FMM, each node memorizes the list of satellite elements which surround the central node. First, an element which contains the point x uDt is searched. This element can be found by searching the satellite elements, in which all the volume coordinates have positive numbers. Then, the advection calculation is performed by the interpolation using Eq. (19) in the searched element. As long as the CFL condition is satisfied, one of satellite elements contains the point x uDt. Therefore, the node-based data structure of the FMM is efficient for searching the element which includes the point x uDt. Furthermore, since this procedure can be performed independently for each node, parallelization is easily implemented in the node-based finite-element calculation of the FMM. 2.3. Verification of accuracy 2.3.1. Advection terms Fig. 3 shows a test model for advection calculation. The model is a rectangular channel with the size 50 · 10 · 10, which is divided into tetrahedral finite elements. The number of nodes is 7928, and the number of elements is 38,842. Initially, the values of nodes which are included in a sphere with radius 3 drawn as a shaded object in the figure are set as 2, whereas those of the other nodes are set as 1. The profiles of nodal values are moved in the x-direction at the velocity u = 1.0. Fig. 4 shows the results of advection calculation of 1000 time steps at time increment Dt = 0.04. In the figure, nodal values in the initial condition and computational results using the first order interpolation, the third order interpolation and the coupling of first and third order interpolations are depicted at the cross section of Z = 5.0, respectively. The figure shows that first order interpolation causes serious numerical diffusion, whereas third order interpolation results in numerical oscillation (Fig. 4 (ii) and (iii)). It is shown that the coupling of the first and third order interpolations gives the best result among the three cases (Fig. 4 (iv)).
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1903
Fig. 3. Test model for advection calculation.
Fig. 4. Numerical results of advection calculation: (i) Initial condition; (ii) first order interpolation; (iii) third order interpolation; (iv) coupling of first and third interpolations.
2.3.2. Computation of incompressible flow In order to confirm the applicability of the present algorithm to incompressible flow simulation, the three-dimensional lid-driven cavity flow is calculated. The model is a cube with the size 1.0 · 1.0 · 1.0 and it is divided into a tetrahedral mesh. The number of nodes is 44,124 and the number of elements is 247,226. The boundary condition of u = 1.0 is set at the top surface and non-slip conditions are at the other surfaces. Specific heat ratio c is set as 1.0 · 10100, which means the flows are almost incompressible. Reynolds number Re is 400 and the time increment Dt is set as 1.0 · 103. Fig. 5 shows that the present results agree well with KuÕs one [15]. A little difference around Z = 0.0 may be caused by insufficiency of resolution of computational grids near the wall. 2.3.3. Computation of compressible flow In order to confirm the applicability of the present algorithm to compressible flow simulation, the SodÕs shock tube problem [16] is taken here. The model is a rectangular channel with the size 1.0 · 0.01 · 0.01. The number of nodes is 44,764 and the number of elements is 212,068. Fig. 6 shows the initial condition.
1904
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 5. Numerical result of lid-driven cavity problem.
p = 1.0 x = 0.0
ρ = 1.0
p = 0.1
ρ = 0.125 x = 1.0
Fig. 6. Initial conditions of the SodÕs shock tube problem.
Fig. 7. Numerical and theoretical results of the SodÕs shock tube problem.
The viscosity coefficient l is set as 0 and the specific heat ratio c is 1.4. The artificial viscosity given by Eq. (16) is added with coefficient a = 1.6. Slip conditions are given at side surfaces. Regarding both ends of the rectangular channel, Dirichlet boundary conditions of (u, p, q) = (0.0, 1.0, 1.0) at X = 0.0 and (u, p, q) = (0.0, 0.1, 0.125) at X = 1.0 are respectively given. The time increment Dt is 2.0 · 104. Fig. 7 is the comparison between the theoretical solution and the computational results after 500 time steps, which shows the good performance of the proposed method.
3. Computation of edge tone 3.1. Air-reed instruments and edge tone Air-reed instrument is a kind of aerophone such as pipe organs, flutes, recorders, bamboo flutes that do not contain mechanically oscillating parts nor sound by human lip like a brass instrument. These instru-
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1905
ments have wedge-shaped structure called edge against the direction of blown jet flow. The blown jet causes self-excited oscillation at the edge and sound is generated. The sounds generated in this way are called edge tone. In air-reed instruments, particular frequencies are selected and amplified by the resonance tube and finally compose the timbre of the musical instruments. The edge tone was precisely studied by Brown [17] from experimental point of view. Tezduyar et al. [18,19] succeeded in computing the edge tone using a finite element fluid analysis. The present authors performed numerical simulations of edge tone both in two-dimensional [20] and in three-dimensional [21] using a hybrid method of incompressible flow analysis and acoustic analogy in the past. In these studies, the CurleÕs equation [22] is used as the acoustic analogy, which is derived from the LighthillÕs theory [23]. By contrast, we conducted a direct numerical simulation of sound based on the compressible Navier–Stokes equations in the present study. 3.2. Brown’s experimental equation Brown [17] has performed precise experiments regarding the edge tone and obtained the following experimental equation: 1 f ¼ 0:466j ð100U 40Þ 0:07 with j ¼ 1:0; 2:3; 3:8; 5:4; ð21Þ 100l where f (Hz) is the frequency of the edge tone, U (m/s) is the velocity of the air jet and l (m) is the distance between the edge and the nozzle of air jet. The parameter j is the oscillation mode determined by U and l. If the distance from the edge is constant, j increases with flow velocity U. Fig. 8 shows the relationship between U, l and f given by Brown [17]. The step-like features in the figure represent the points at which the value of j changes. 3.3. Direct numerical simulation of edge tone A direct numerical simulation of the edge tone is performed here using the proposed method. Fig. 9 shows the analysis model for simulation of edge tone. The width of the nozzle and the shape of the edge
Fig. 8. Relationship between velocity of air jet, distance between nozzle and edge, and frequency of edge tone [17].
1906
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
are similar to those of BrownÕs experimental model [17]. The thickness of the model is 1 (mm), and the slip conditions are given at side surfaces. The model is divided into tetrahedral finite element mesh as shown in Fig. 10. The number of nodes is 236,675 and the number of elements is 1,128,418. The distance between the nozzle and the edge is 6 (mm) and the velocity of the blown jet is 5 (m/s). The width of nozzle and the shapes of edge are set to be the same as those of BrownÕs experiment. The properties of fluid are set as air, i.e. static sound pressure p0 = 1.01325 · 105 (Pa), density q0 = 1.184 (kg/m3), coefficient of viscosity l = 1.71 · 105 (Pa s), specific heat ratio c = 1.4, Prandtl number Pr = 0.72 and gas constant R = 8.31451.
Fig. 9. Analysis model for simulation of edge tone (unit: mm).
Fig. 10. Finite element mesh for simulation of the edge tone.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1907
In this calculation, the multi-time-step integration technique [14] is employed. The time increment of calculation of advection terms is set as Dt = 1.0 · 106 (s), whereas that of non-advection terms is set as Dt = 1.0 · 107 (s). Employed are 16 nodes (128 processors) of the Hitachi SR8000 supercomputer, and the computation of 30,000 time steps (0.03 (s) took approximately 72 h. Fig. 11 shows the distributions of velocity at Y = 0.5 (mm): The upper is the result at t = 0.020 (s) and the lower is the result at t = 0.021 (s). It is seen that the blown air jet oscillates at the upper and lower sides of the edge.
Fig. 11. Distributions of velocity at the cross-section Y = 0.5. Upper: t = 0.020 (s); lower: t = 0.021 (s).
Fig. 12. Distributions of pressure at the cross-section Y = 0.5 (mm). Left: t = 0.020 (s); right: t = 0.021 (s).
1908
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 13. Sound wave obtained at (X, Z) = (9, 15) (mm): 9 (mm) in the horizontal direction and 15 (mm) in the vertical direction from tip of the edge.
Fig. 14. Power spectrum of the sound wave.
Fig. 12 shows the distributions of pressure at Y = 0.5 (mm): The left is the result at t = 0.020 (s) whereas the right is the result at t = 0.021 (s). By comparing Fig. 12 with Fig. 11, it is seen that the up and down oscillation of the air jet is observed around the edge tip, which is considered to cause the alternate fluctuation of the pressure and then the sound source of the edge tone. Fig. 13 shows the sound wave observed at (X, Z) = (9, 15) (mm): 9 (mm) in the horizontal direction and 15 (mm) in the vertical direction from the tip of the edge. Since this study is a direct numerical simulation of aerodynamic sound, the profile of sound pressure is directly obtained by measuring the pressure of flow field at the observation point. The perturbation of pressure derived by the compressible Navier–Stokes equations is the sound itself. Fig. 14 shows the power spectrum of sound wave obtained by FFT analysis, which is employed after 0.01 (s), when the sound wave has stabilized. The figure shows that 335.69 (Hz) is the strongest frequency component. If U = 5.0, l = 0.006 and j = 1.0 (1st mode) are substituted into the BrownÕs experimental equation, the resulted frequency f becomes 342.26 (Hz) and it is concluded that the frequency obtained from the BrownÕs experimental equation is well coincident with the result of the present direct numerical simulation.
4. Conclusion Simulated is the aerodynamic sound using a node-based finite-element method combined with the CCUP method and the CIVA method on parallel computers. The node-based data structure of the FMM
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1909
is suitable for solving advection terms using the CIVA method. The three-dimensional lid-driven cavity problem and the SodÕs shock tube problem are calculated as verification of accuracy, and reasonable results are obtained. As a numerical example of DNS of aerodynamic sound, an edge tone is successfully simulated by the proposed method, agreed well with the BrownÕs experimental equation. The FMM is especially effective in the adaptive analysis on parallel computers, since the method can generate local finite elements robustly even on massively parallel processors. The obtained results in this study shows potentiality of DNS of aerodynamic sounds generated from arbitrary three-dimensional objects with complex shapes using unstructured computational grids.
Acknowledgements A part of this research was performed in the ÔFrontier Simulation Software for Industrial Science (FSIS)Õ project supported by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT).
References [1] J.C. Hardin, S.L. Lamkin, Aeroacoustic computation of cylinder wake flow, AIAA J. 22 (1984) 51–57. [2] W.Z. Shen, J.N. Sørensen, Comment on the aeroacoustic formulation of Hardin and Pope, AIAA J. 37 (1999) 141–143. [3] S.A. Slimon, M.C. Soteriou, D.W. Davis, Computational aeroacoustics simulations using the expansion about incompressible flow approach, AIAA J. 37 (1999) 409–416. [4] T. Yabe, P.Y. Wang, Unified numerical procedure for compressible and incompressible fluid, J. Phys. Soc. Jpn 60 (1991) 2105– 2108. [5] T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comput. Phys. Commun. 66 (1991) 219–232. [6] T. Yabe, T. Ishikawa, P.Y. Wang, T. Aoki, Y. Kadota, F. Ikeda, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comput. Phys. Commun. 66 (1991) 233–242. [7] N. Tanaka, The CIVA method for mesh-free approaches: improvement of the CIP method for n-simplex, Comput. Fluid Dyn. J. 8–1 (1999) 121–127. [8] G. Yagawa, T. Yamada, Free mesh method: a new meshless finite element method, Int. J. Comput. Mech. 18 (1996) 383–386. [9] G. Yagawa, T. Hosokawa, Application of the free mesh method to the 3-dimensional problem, Trans. Jpn. Soc. Mech. Eng. (A) 63–614 (1998) 2251–2256 (in Japanese). [10] G. Yagawa, Node-by-node parallel finite elements: a virtually meshless method, Int. J. Numer. Meth. Engrg 60 (2004) 69–102. [11] T. Fujisawa, M. Inaba, G. Yagawa, Parallel computing of high-speed compressible flows using a node-based finite-element method, Int. J. Numer. Meth. Engrg 58 (2003) 481–511. [12] J. Von Neumann, R.D. Richtmyer, A method for numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 232– 237. [13] T.J. Poinsot, S.K. Lele, Boundary conditions for direct simulations of compressible viscous flow, J. Comput. Phys. 101 (1992) 104–129. [14] M. Ida, An improved unified solver for compressible and incompressible fluids involving free surfaces: II. Multi-time-step integration and applications, Comput. Phys. Commun. 150 (2003) 300–322. [15] H.C. Ku, R.S. Hirsh, T.D. Taylor, A pseudospectral method for solution of the three-dimensional incompressible Navier–Stokes equations, J. Comput. Phys. 70 (1987) 439–462. [16] G.A. Sod, A survey of several finite difference methods for systems of nonlinear conservation laws, J. Comput. Phys. 27 (1978) 1– 31. [17] G.B. Brown, The vortex motion causing edge tones, Proc. Phys. Soc. London XLIX (1937) 493–507. [18] T.E. Tezduyar, J. Liou, Grouped element-by-element iteration schemes for incompressible flow computations, Comput. Phys. Commun. 53 (1989) 441–453. [19] J. Liou, T.E. Tezduyar, Iterative adaptive implicit–explicit methods for flow problems, Int. J. Numer. Meth. Fluids 11 (1990) 867– 880. [20] G. Yagawa, M. Shirazaki, Y. Hizaki, Analysis of sound generation mechanism of air-reed instruments with finite element method, in: Proc. the 13th CFD symposium, 47, 1999, (in Japanese).
1910
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
[21] J. Tsuchida, S. Ito, T. Fujisawa, G. Yagawa, Computational fluid dynamics on sounding mechanism of air-reed instruments, Trans. Jpn. Soc. Mechan. Eng. (B) 69–678 (2003) 32–37 (in Japanese). [22] N. Curle, The influence of solid boundaries upon aerodynamic sound, Proc. Roy. Soc. 231 (1955) 505–514. [23] M.J. Lighthill, On sound generated aerodynamically I. General theory, Proc. Roy. Soc. London A 211 (1952) 564–587.
Direct numerical simulation of aerodynamic sounds by a compressible CFD scheme with node-by-node finite elements J. Tsuchida a, T. Fujisawa b, G. Yagawa a
c,*
Toyota Motor Corporation, 1 Toyota-cho, Toyota City, Aichi Prefecture 471-8571, Japan b Prometech Software Inc., 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan c Tokyo University, 1-32-19, Akatsutsumi, Setagaya-ku, Tokyo 156-0044, Japan
Received 5 November 2004; received in revised form 5 April 2005; accepted 26 May 2005
Abstract A direct numerical simulation (DNS) of aerodynamic sounds is performed by using a compressible CFD scheme with node-by-node finite elements. Although the computational simulations of aerodynamic sounds generated from three-dimensional objects with complex boundary shapes are required in various industrial applications, such simulations involve a lot of difficulties regarding instability in numerical schemes and enormous computational costs. In this paper, the present authors have combined the C-CUP (CIP-Combined Unified Procedure, CIP: Constrained Interpolation Profile) method and the CIVA (Cubic Interpolation with Volume/Area Coordinates) method with the Free Mesh Method (FMM) in order to achieve efficient simulation of aerodynamic sounds using unstructured computational grids of finite elements on parallel computers. As a numerical example, the edge tone is successfully computed using the present method with parallel computing with Hitachi SR8000 supercomputer. 2005 Elsevier B.V. All rights reserved. Keywords: Computational fluid dynamics; Finite element method; Edge tone; Sounding mechanism
1. Introduction Aerodynamic sound analyses are usually performed by employing some approximation models in order to save computational costs. The approximation models often used can be categorized into the following *
Corresponding author. Tel./fax: +81 3 3327 5680. E-mail address: [email protected] (G. Yagawa).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.05.039
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1897
two kinds: One is the hybrid method [1], which combines incompressible fluid analysis and acoustic analogy. This method is the most practical approach for calculation of aerodynamic sounds from industrial point of view. Another is the acoustic/viscous splitting method [2,3], which computes average flows and a gap from average flows separately. The average flows are computed by incompressible fluid analysis and a gap from average flows is computed by solving disturbance equations. However, the progress of computers in recent years has made it possible to compute aerodynamic sound by solving the compressible Navier–Stokes equations directly without introducing any approximation models, i.e. DNS (Direct Numerical Simulation of aerodynamic sound). The DNS has a fundamental merit to capture the interference between average flows and sound waves, especially when the resonance occurs in the fluid field. Nevertheless, the computational costs required in the DNS is very large, so that the DNS reported so far was limited to two-dimensional objects with simple boundary shapes. The DNS of threedimensional object with complex boundary shapes is still difficult to solve due to enormous computational costs. In order to overcome this situation, this study aims at the DNS of aerodynamic sound generated from three-dimensional objects, using unstructured computational grids of finite elements on parallel computers. We employ unstructured computational grids to concentrate computational power on the critical parts and then solve the problem efficiently even if the objects have complex boundary shapes. Since the sound is a phenomenon of propagation of compressible waves, both generation and propagation of aerodynamic sound can be captured by compressible fluid analysis from theoretical point of view. However, conventional density-based numerical schemes for compressible flows cause so-called stiffness problem under low Mach number and numerical simulations confront serious instability. In order to overcome this problem, we employ the C-CUP (CIP-Combined Unified Procedure, CIP: Constrained Interpolation Profile) method [4], which can treat both incompressible and compressible flows simultaneously. Therefore, it can perform compressible flow analysis under low Mach number. The C-CUP method is one of the splitting methods, which compute the governing equations by separating advection and non-advection terms. The method originally based on the finite difference method, and the advection terms are solved by the CIP (Constrained Interpolation Profile) method [5,6]. However, in the present study, we employ unstructured computational grids of finite elements, so that the CIVA (Cubic Interpolation with Volume/Area Coordinates) method [7] is also introduced. The CIVA method is one of the third-order interpolation methods for triangular or tetrahedral elements, which is a generalized CIP method using unstructured computational grids. In this study, the advection terms are solved by the CIVA method, making use of the node-based data structure of the finite elements in the FMM (Free Mesh Method) [8–10], which makes it possible to solve advection terms efficiently on large-scale finite element meshes.
2. Compressible flow simulation scheme based on finite elements 2.1. Outline of the numerical method Parallel computing methods for finite element analysis can be classified into element-based and nodebased ones. In the present study, we use the free mesh method (FMM) [8–10], which is a node-based finite-element method. In the FMM, a node-based data structure is employed to memorize sparse matrices and connectivity data between nodes (Fig. 1). The FMM performs both mesh generation and solution of equations in a seamless manner. Therefore, it is an efficient method for adaptive analysis on parallel computers. Although a fixed computational grid is employed in the present study, node-based data structures of the FMM are used in order to perform parallel computing effectively on distributed memory systems [11].
1898
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 1. Free mesh method.
On the other hand, the C-CUP method [4] is one of the splitting methods, which solves governing equation separating into the advection and the non-advection parts. The following is a summarized procedure of the C-CUP method. 1. 2. 3. 4.
Solve the advection parts using CIP method. Solve the pressure equation. Compute the velocity and the density of the next time step. Renew the spatial derivatives of primitive physical values. The governing equations are oq þ u rq ¼ qr u; ot ou rp þ u ru ¼ þ Qu ; ot q op þ u rp ¼ qC 2S u; ot
ð1Þ ð2Þ ð3Þ
where u is the velocity, p is the pressure, q is the density and CS is the local sound speed. Qu expresses the force, which includes gravity, buoyancy, elastic stress tensor and surface tension. oq þ u rq ¼ 0; ot ou þ u ru ¼ 0; ot op þ u rp ¼ 0. ot
ð4Þ ð5Þ ð6Þ
These equations are solved by the CIP method or other similar methods. In this study, we solve the equations using a new method based on the CIVA method and the node-based data structure of the FMM, which is described later.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1899
The following are the non-advection parts: oq ¼ qr u; ot ou rp ¼ þ Qu ; ot q op ¼ qC 2S u. ot
ð7Þ ð8Þ ð9Þ
By discretizing the time derivatives on the LHS of Eqs. (7)–(9) and estimating the spatial derivatives on the RHS from the values at the next time step, we obtain qnþ1 q ¼ q r unþ1 ; Dt unþ1 u rpnþ1 ¼ þ Qu ; Dt q pnþ1 p nþ1 ¼ q C 2 ; S u Dt
ð10Þ ð11Þ ð12Þ
where n is the time step and Dt is the increment of time. The superscript * indicates the value after solving the advection parts. By taking divergence of Eq. (11) and substituting it into Eq. (12), we obtain the following pressure equation: nþ1 rp pnþ1 p ru r . ð13Þ þ ¼ 2 q Dt q C 2 S Dt The velocity un+1 and the density qn+1 are respectively computed by the following equations: nþ1 rp unþ1 ¼ u Dt Q u ; q pnþ1 p . qnþ1 ¼ q þ C 2 S
ð14Þ ð15Þ
The artificial viscosity is needed in order to stabilize the solution involving shock waves, when inconservative governing equation is used. Therefore, in this study, the following von Neumann–Richtmyer artificial viscosity [12] is employed: 8 < a q C c þ 1 lr u lr u if r u < 0; S 2 ð16Þ qv ¼ : 0 otherwise; where a is a constant parameter, c is the specific heat ratio and l is the representative length of each element. In this study, the average distance between the central node and the associated satellite nodes is used as the representative length of elements (see Fig. 1). In order to compute aerodynamic sound precisely, non-physical numerical reflection of waves at the boundaries must be suppressed. For this purpose, the NSCBC (Navier–Stokes Characteristic Boundary Condition) [13] is employed and non-reflecting boundary conditions are set at outflow boundaries. Since flow simulation under low Mach number includes extremely different time scales between the advection and non-advection parts, the multi-time-step integration technique [14] is employed. In the multitime-step integration, the advection and non-advection parts are solved by different time increments in order to perform calculation efficiently.
1900
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Under consideration of artificial viscosity, the NSCBC and the multi-time-step integration technique, the numerical algorithm we propose in this study is summarized as follows: 1. Determine the boundary conditions using the NSCBC. The values of the boundary nodes are fixed to the values computed by the NSCBC (only at the first step): ðun ; pn ; qn ; run ; rpn ; rqn Þ. 2. Solve the advection parts: ðu ; p ; q ; ru ; rp ; rq Þ. 3. Solve the pressure equation using iterative method such as Bi-CGSTAB: ðu ; p ; q ; ru ; rp ; rq Þ. 4. Renew the velocity: rp u ¼ u Dt ; q ðu ; p ; q ; ru ; rp ; rq Þ. 5. Renew the density: qnþ1 ¼ q þ
ðp p Þ ; C 2 S
ðu ; p ; qnþ1 ; ru ; rp ; rq Þ. 6. Add the physical and artificial viscosity to the velocity obtained above: unþ1 ¼ u þ DtQu ; Qu
¼
qnþ1
osij oqv ; oxj oxi
1
oui ouj þ sij ¼ l oxj oxi
2 ouk ; dij l 3 oxk
ðunþ1 ; p ; qnþ1 ; ru ; rp ; rq Þ. 7. Add the artificial viscosity to the pressure obtained above: pnþ1 ¼ p þ DtQp ; Qp ¼ ðc 1Þqv ru ; ðunþ1 ; pnþ1 ; qnþ1 ; ru ; rp ; rq Þ. 8. Renew the spatial derivatives: ðunþ1 ; pnþ1 ; qnþ1 ; runþ1 ; rpnþ1 ; rqnþ1 Þ. 9. Determine the boundary conditions using the NSCBC. The values of the boundary nodes are fixed to the values computed by the NSCBC: ðunþ1 ; pnþ1 ; qnþ1 ; runþ1 ; rpnþ1 ; rqnþ1 Þ.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1901
10. Repeat Steps 3–9 several times, if the multi-time-step integration technique is employed. 11. Return to Step 2. In the above procedure, l is the viscosity coefficient and dij is the Kronecker delta. 2.2. Calculation of advection terms Regarding the calculation of advection terms, we use a new method combining the CIVA method [7] and the node-based data structure of the FMM. The advection equation is expressed in the hyperbolic equation as follows: of of þu ¼ 0. ot ox
ð17Þ
If the time increment Dt is small enough, it can be assumed that the velocity is constant during the time step. Then, the solution of the advection equation is given by f ðx; t þ DtÞ ¼ f ðx uDt; tÞ.
ð18Þ
Eq. (18) means that if we need the values of x at t + Dt, we can obtain the values by computing x uDt at t. The value at x uDt can be obtained by interpolating the nodal values of the element which contains the point x uDt. In the CIP method, both the physical values and their spatial derivatives are memorized at each node, and the physical values between two nodes are interpolated by the cubic polynomial, so that advection terms can be calculated with high accuracy. On the other hand, since we use unstructured computational grids of the finite elements, the CIVA method is used instead of the original CIP method. The CIVA method approximates the third order interpolation within triangle or tetrahedral elements. Interpolation in the CIVA method is expressed as f~ ðr1 ; r2 ; r3 ; r4 Þ ¼
4 X i¼1
fi ri þ d
4 X
bjk fr2j rk þ cðr1 r2 r3 þ r2 r3 r4 þ r1 r2 r4 þ r1 r3 r4 Þg
ð19Þ
j;k¼1ðj6¼kÞ
with bjk ¼ fj fk þ ðxk xj Þfjx þ ðy k y j Þfjy þ ðzk zj Þfjz
0 6 c 6 0:5;
d ¼ 0 or 1
where f~ is the interpolated value, f1, f2, f3 and f4 are the nodal values and r1, r2, r3 and r4 are the volume coordinates. Here, the parameter c cannot be determined uniquely from the nodal values, because both physical values and its spatial derivative are zero at each node of tetrahedral elements and the terms of ri, rj, rk in the RHS of Eq. (19) is also zero. Then, we employ constant curvature conditions to determine c. In this case, the constant parameter c = 0.5 could be used for all elements. If the parameter c less than 0.5 is used, the numerical viscosity is added. The parameter d is used for a switching between the first and the third order interpolations. Generally speaking, numerical oscillation occurs where physical values change suddenly. In order to stabilize numerical oscillation, the interpolation is switched from the third to the first order where physical values change suddenly. A criterion of the switching is given by the following equations: ( nþ1 1 f min 6 f~ 6 fmax ; d¼ ð20Þ 0 otherwise;
1902
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 2. Advection calculation using the CIVA method.
where fmin ¼ minðf1n ; f2n ; f3n ; f4n Þ; fmax ¼ maxðf1n ; f2n ; f3n ; f4n Þ; nþ1 f~ is the third order interpolated value.
Although Eq. (20) has no mathematical background, several experiments show this stabilization method works well, if the computational grid is fine enough. Fig. 2 shows the basic concept of solving advection terms using the CIVA method and the node-based data structure of the FMM. In the FMM, each node memorizes the list of satellite elements which surround the central node. First, an element which contains the point x uDt is searched. This element can be found by searching the satellite elements, in which all the volume coordinates have positive numbers. Then, the advection calculation is performed by the interpolation using Eq. (19) in the searched element. As long as the CFL condition is satisfied, one of satellite elements contains the point x uDt. Therefore, the node-based data structure of the FMM is efficient for searching the element which includes the point x uDt. Furthermore, since this procedure can be performed independently for each node, parallelization is easily implemented in the node-based finite-element calculation of the FMM. 2.3. Verification of accuracy 2.3.1. Advection terms Fig. 3 shows a test model for advection calculation. The model is a rectangular channel with the size 50 · 10 · 10, which is divided into tetrahedral finite elements. The number of nodes is 7928, and the number of elements is 38,842. Initially, the values of nodes which are included in a sphere with radius 3 drawn as a shaded object in the figure are set as 2, whereas those of the other nodes are set as 1. The profiles of nodal values are moved in the x-direction at the velocity u = 1.0. Fig. 4 shows the results of advection calculation of 1000 time steps at time increment Dt = 0.04. In the figure, nodal values in the initial condition and computational results using the first order interpolation, the third order interpolation and the coupling of first and third order interpolations are depicted at the cross section of Z = 5.0, respectively. The figure shows that first order interpolation causes serious numerical diffusion, whereas third order interpolation results in numerical oscillation (Fig. 4 (ii) and (iii)). It is shown that the coupling of the first and third order interpolations gives the best result among the three cases (Fig. 4 (iv)).
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1903
Fig. 3. Test model for advection calculation.
Fig. 4. Numerical results of advection calculation: (i) Initial condition; (ii) first order interpolation; (iii) third order interpolation; (iv) coupling of first and third interpolations.
2.3.2. Computation of incompressible flow In order to confirm the applicability of the present algorithm to incompressible flow simulation, the three-dimensional lid-driven cavity flow is calculated. The model is a cube with the size 1.0 · 1.0 · 1.0 and it is divided into a tetrahedral mesh. The number of nodes is 44,124 and the number of elements is 247,226. The boundary condition of u = 1.0 is set at the top surface and non-slip conditions are at the other surfaces. Specific heat ratio c is set as 1.0 · 10100, which means the flows are almost incompressible. Reynolds number Re is 400 and the time increment Dt is set as 1.0 · 103. Fig. 5 shows that the present results agree well with KuÕs one [15]. A little difference around Z = 0.0 may be caused by insufficiency of resolution of computational grids near the wall. 2.3.3. Computation of compressible flow In order to confirm the applicability of the present algorithm to compressible flow simulation, the SodÕs shock tube problem [16] is taken here. The model is a rectangular channel with the size 1.0 · 0.01 · 0.01. The number of nodes is 44,764 and the number of elements is 212,068. Fig. 6 shows the initial condition.
1904
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 5. Numerical result of lid-driven cavity problem.
p = 1.0 x = 0.0
ρ = 1.0
p = 0.1
ρ = 0.125 x = 1.0
Fig. 6. Initial conditions of the SodÕs shock tube problem.
Fig. 7. Numerical and theoretical results of the SodÕs shock tube problem.
The viscosity coefficient l is set as 0 and the specific heat ratio c is 1.4. The artificial viscosity given by Eq. (16) is added with coefficient a = 1.6. Slip conditions are given at side surfaces. Regarding both ends of the rectangular channel, Dirichlet boundary conditions of (u, p, q) = (0.0, 1.0, 1.0) at X = 0.0 and (u, p, q) = (0.0, 0.1, 0.125) at X = 1.0 are respectively given. The time increment Dt is 2.0 · 104. Fig. 7 is the comparison between the theoretical solution and the computational results after 500 time steps, which shows the good performance of the proposed method.
3. Computation of edge tone 3.1. Air-reed instruments and edge tone Air-reed instrument is a kind of aerophone such as pipe organs, flutes, recorders, bamboo flutes that do not contain mechanically oscillating parts nor sound by human lip like a brass instrument. These instru-
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1905
ments have wedge-shaped structure called edge against the direction of blown jet flow. The blown jet causes self-excited oscillation at the edge and sound is generated. The sounds generated in this way are called edge tone. In air-reed instruments, particular frequencies are selected and amplified by the resonance tube and finally compose the timbre of the musical instruments. The edge tone was precisely studied by Brown [17] from experimental point of view. Tezduyar et al. [18,19] succeeded in computing the edge tone using a finite element fluid analysis. The present authors performed numerical simulations of edge tone both in two-dimensional [20] and in three-dimensional [21] using a hybrid method of incompressible flow analysis and acoustic analogy in the past. In these studies, the CurleÕs equation [22] is used as the acoustic analogy, which is derived from the LighthillÕs theory [23]. By contrast, we conducted a direct numerical simulation of sound based on the compressible Navier–Stokes equations in the present study. 3.2. Brown’s experimental equation Brown [17] has performed precise experiments regarding the edge tone and obtained the following experimental equation: 1 f ¼ 0:466j ð100U 40Þ 0:07 with j ¼ 1:0; 2:3; 3:8; 5:4; ð21Þ 100l where f (Hz) is the frequency of the edge tone, U (m/s) is the velocity of the air jet and l (m) is the distance between the edge and the nozzle of air jet. The parameter j is the oscillation mode determined by U and l. If the distance from the edge is constant, j increases with flow velocity U. Fig. 8 shows the relationship between U, l and f given by Brown [17]. The step-like features in the figure represent the points at which the value of j changes. 3.3. Direct numerical simulation of edge tone A direct numerical simulation of the edge tone is performed here using the proposed method. Fig. 9 shows the analysis model for simulation of edge tone. The width of the nozzle and the shape of the edge
Fig. 8. Relationship between velocity of air jet, distance between nozzle and edge, and frequency of edge tone [17].
1906
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
are similar to those of BrownÕs experimental model [17]. The thickness of the model is 1 (mm), and the slip conditions are given at side surfaces. The model is divided into tetrahedral finite element mesh as shown in Fig. 10. The number of nodes is 236,675 and the number of elements is 1,128,418. The distance between the nozzle and the edge is 6 (mm) and the velocity of the blown jet is 5 (m/s). The width of nozzle and the shapes of edge are set to be the same as those of BrownÕs experiment. The properties of fluid are set as air, i.e. static sound pressure p0 = 1.01325 · 105 (Pa), density q0 = 1.184 (kg/m3), coefficient of viscosity l = 1.71 · 105 (Pa s), specific heat ratio c = 1.4, Prandtl number Pr = 0.72 and gas constant R = 8.31451.
Fig. 9. Analysis model for simulation of edge tone (unit: mm).
Fig. 10. Finite element mesh for simulation of the edge tone.
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1907
In this calculation, the multi-time-step integration technique [14] is employed. The time increment of calculation of advection terms is set as Dt = 1.0 · 106 (s), whereas that of non-advection terms is set as Dt = 1.0 · 107 (s). Employed are 16 nodes (128 processors) of the Hitachi SR8000 supercomputer, and the computation of 30,000 time steps (0.03 (s) took approximately 72 h. Fig. 11 shows the distributions of velocity at Y = 0.5 (mm): The upper is the result at t = 0.020 (s) and the lower is the result at t = 0.021 (s). It is seen that the blown air jet oscillates at the upper and lower sides of the edge.
Fig. 11. Distributions of velocity at the cross-section Y = 0.5. Upper: t = 0.020 (s); lower: t = 0.021 (s).
Fig. 12. Distributions of pressure at the cross-section Y = 0.5 (mm). Left: t = 0.020 (s); right: t = 0.021 (s).
1908
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
Fig. 13. Sound wave obtained at (X, Z) = (9, 15) (mm): 9 (mm) in the horizontal direction and 15 (mm) in the vertical direction from tip of the edge.
Fig. 14. Power spectrum of the sound wave.
Fig. 12 shows the distributions of pressure at Y = 0.5 (mm): The left is the result at t = 0.020 (s) whereas the right is the result at t = 0.021 (s). By comparing Fig. 12 with Fig. 11, it is seen that the up and down oscillation of the air jet is observed around the edge tip, which is considered to cause the alternate fluctuation of the pressure and then the sound source of the edge tone. Fig. 13 shows the sound wave observed at (X, Z) = (9, 15) (mm): 9 (mm) in the horizontal direction and 15 (mm) in the vertical direction from the tip of the edge. Since this study is a direct numerical simulation of aerodynamic sound, the profile of sound pressure is directly obtained by measuring the pressure of flow field at the observation point. The perturbation of pressure derived by the compressible Navier–Stokes equations is the sound itself. Fig. 14 shows the power spectrum of sound wave obtained by FFT analysis, which is employed after 0.01 (s), when the sound wave has stabilized. The figure shows that 335.69 (Hz) is the strongest frequency component. If U = 5.0, l = 0.006 and j = 1.0 (1st mode) are substituted into the BrownÕs experimental equation, the resulted frequency f becomes 342.26 (Hz) and it is concluded that the frequency obtained from the BrownÕs experimental equation is well coincident with the result of the present direct numerical simulation.
4. Conclusion Simulated is the aerodynamic sound using a node-based finite-element method combined with the CCUP method and the CIVA method on parallel computers. The node-based data structure of the FMM
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
1909
is suitable for solving advection terms using the CIVA method. The three-dimensional lid-driven cavity problem and the SodÕs shock tube problem are calculated as verification of accuracy, and reasonable results are obtained. As a numerical example of DNS of aerodynamic sound, an edge tone is successfully simulated by the proposed method, agreed well with the BrownÕs experimental equation. The FMM is especially effective in the adaptive analysis on parallel computers, since the method can generate local finite elements robustly even on massively parallel processors. The obtained results in this study shows potentiality of DNS of aerodynamic sounds generated from arbitrary three-dimensional objects with complex shapes using unstructured computational grids.
Acknowledgements A part of this research was performed in the ÔFrontier Simulation Software for Industrial Science (FSIS)Õ project supported by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT).
References [1] J.C. Hardin, S.L. Lamkin, Aeroacoustic computation of cylinder wake flow, AIAA J. 22 (1984) 51–57. [2] W.Z. Shen, J.N. Sørensen, Comment on the aeroacoustic formulation of Hardin and Pope, AIAA J. 37 (1999) 141–143. [3] S.A. Slimon, M.C. Soteriou, D.W. Davis, Computational aeroacoustics simulations using the expansion about incompressible flow approach, AIAA J. 37 (1999) 409–416. [4] T. Yabe, P.Y. Wang, Unified numerical procedure for compressible and incompressible fluid, J. Phys. Soc. Jpn 60 (1991) 2105– 2108. [5] T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comput. Phys. Commun. 66 (1991) 219–232. [6] T. Yabe, T. Ishikawa, P.Y. Wang, T. Aoki, Y. Kadota, F. Ikeda, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comput. Phys. Commun. 66 (1991) 233–242. [7] N. Tanaka, The CIVA method for mesh-free approaches: improvement of the CIP method for n-simplex, Comput. Fluid Dyn. J. 8–1 (1999) 121–127. [8] G. Yagawa, T. Yamada, Free mesh method: a new meshless finite element method, Int. J. Comput. Mech. 18 (1996) 383–386. [9] G. Yagawa, T. Hosokawa, Application of the free mesh method to the 3-dimensional problem, Trans. Jpn. Soc. Mech. Eng. (A) 63–614 (1998) 2251–2256 (in Japanese). [10] G. Yagawa, Node-by-node parallel finite elements: a virtually meshless method, Int. J. Numer. Meth. Engrg 60 (2004) 69–102. [11] T. Fujisawa, M. Inaba, G. Yagawa, Parallel computing of high-speed compressible flows using a node-based finite-element method, Int. J. Numer. Meth. Engrg 58 (2003) 481–511. [12] J. Von Neumann, R.D. Richtmyer, A method for numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 232– 237. [13] T.J. Poinsot, S.K. Lele, Boundary conditions for direct simulations of compressible viscous flow, J. Comput. Phys. 101 (1992) 104–129. [14] M. Ida, An improved unified solver for compressible and incompressible fluids involving free surfaces: II. Multi-time-step integration and applications, Comput. Phys. Commun. 150 (2003) 300–322. [15] H.C. Ku, R.S. Hirsh, T.D. Taylor, A pseudospectral method for solution of the three-dimensional incompressible Navier–Stokes equations, J. Comput. Phys. 70 (1987) 439–462. [16] G.A. Sod, A survey of several finite difference methods for systems of nonlinear conservation laws, J. Comput. Phys. 27 (1978) 1– 31. [17] G.B. Brown, The vortex motion causing edge tones, Proc. Phys. Soc. London XLIX (1937) 493–507. [18] T.E. Tezduyar, J. Liou, Grouped element-by-element iteration schemes for incompressible flow computations, Comput. Phys. Commun. 53 (1989) 441–453. [19] J. Liou, T.E. Tezduyar, Iterative adaptive implicit–explicit methods for flow problems, Int. J. Numer. Meth. Fluids 11 (1990) 867– 880. [20] G. Yagawa, M. Shirazaki, Y. Hizaki, Analysis of sound generation mechanism of air-reed instruments with finite element method, in: Proc. the 13th CFD symposium, 47, 1999, (in Japanese).
1910
J. Tsuchida et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1896–1910
[21] J. Tsuchida, S. Ito, T. Fujisawa, G. Yagawa, Computational fluid dynamics on sounding mechanism of air-reed instruments, Trans. Jpn. Soc. Mechan. Eng. (B) 69–678 (2003) 32–37 (in Japanese). [22] N. Curle, The influence of solid boundaries upon aerodynamic sound, Proc. Roy. Soc. 231 (1955) 505–514. [23] M.J. Lighthill, On sound generated aerodynamically I. General theory, Proc. Roy. Soc. London A 211 (1952) 564–587.