Computation of aeroacoustic sources for a Gulfstream G550 nose landing gear model using adaptive FEM

Computers and Fluids 124 (2016) 136–146

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Computation of aeroacoustic sources for a Gulfstream G550 nose landing gear model using adaptive FEM Rodrigo Vilela de Abreu a,∗, Niclas Jansson b, Johan Hoffman b a b

Computational Technology Laboratory and Linné FLOW Centre, CSC, KTH, Stockholm, SE-10044 Sweden Computational Technology Laboratory, CSC, KTH, Stockholm, SE-10044 Sweden

a r t i c l e

i n f o

Article history: Received 20 May 2014 Revised 29 July 2015 Accepted 26 October 2015 Available online 31 October 2015 Keywords: Landing gear noise Computational fluid dynamics Computational aeroacoustics Adaptive finite element methods Turbulence CAA CFD FEM

a b s t r a c t This work presents a direct comparison of unsteady, turbulent flow simulations with measurements performed using a Gulfstream G550 nose landing gear model. The experimental campaign, which was carried out by researchers from the NASA Langley Research Center, provided a series of detailed, well documented wind-tunnel measurements for comparison and validation of computational fluid dynamics (CFD) and computational aeroacoustics (CAA) methodologies. Several computational efforts were collected and presented at the Benchmark for Airframe Noise Computation workshops, BANC-I and II. For our simulations, we used a General Galerkin finite element method (G2), where no explicit subgrid model is used, and where the computational mesh is adaptively refined with respect to a posteriori estimates of the error in a quantity of interest, here the source term in Lighthill’s equation. The mesh is fully unstructured and the solution is time-resolved, which are key ingredients for solving problems of industrial relevance in the field of aeroacoustics. Moreover, we choose to model the boundary layers on the landing gear geometry with a free-slip condition for the velocity, which we previously observed to produce good results for external flows at high Reynolds numbers, and which considerably reduces the amount of cells required in the mesh. The comparisons presented here are an attempt to quantify the accuracy of our models, methods and assumptions; to that end, several results containing both time-averaged and unsteady flow quantities, always side by side with corresponding experimental values, are reported. The main finding is that we are able to simulate a complex, unsteady flow problem using a parameter-free methodology developed for high Reynolds numbers, external aerodynamics and aeroacoustics applications. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Since 2010, our research group has contributed with incompressible flow computations to the Benchmark for Airframe Noise Computation workshops, BANC-I and II, where our contributions clearly differ from other workshop contributions in what regards the computational method. We use a General Galerkin (G2) finite element (FE) method and since our first simulation with the Rudimentary Landing Gear (RLG) geometry by Spalart et al. at Boeing [1,2], we have emphasized that adaptive mesh generation, combined with unstructured meshes and simple wall-layer models, is a strong set of tools for computing the flow around complex geometries at high Reynolds numbers (Re). In this way, we (1). eliminate the need for a boundarylayer mesh, since we allow slip velocities which are related to the wall shear stresses by a “friction coefficient”, as in earlier Large Eddy

∗

Corresponding author. Tel.: +4687906267. E-mail address: [email protected] (R. Vilela de Abreu).

http://dx.doi.org/10.1016/j.compfluid.2015.10.017 0045-7930/© 2015 Elsevier Ltd. All rights reserved.

Simulation (LES) wall-layer models [3]; (2). eliminate the need for ad hoc mesh generation, which is the standard technique employed by the majority of the CFD/CAA community, including most of the participants of the BANC workshops [4,5]; (3). are able to simulate very complex geometries with ease, without the need to use overset or block meshes to describe the geometry. We believe our methodology is therefore highly suitable for several industrial applications, including CAA, where complex geometries and intricate unsteady flows are common. At the BANC-II workshop following the 18th AIAA/CEAS Aeroacoustics Conference (33rd AIAA Aeroacoustics Conference) in Colorado Springs, 2012, most of the participants in Category 4 (the flow past a partially-dressed cavity-closed nose landing gear) used unstructured meshes (except the Exa Corporation group [4], who used the Lattice–Boltzmann method and, therefore, a blockstructured Cartesian mesh). Furthermore, at the workshop, Vatsa et al. [6] introduced mesh adaptivity to their CFD/CAA methodology for computing the flow around the landing gear. We consider this an indication that the methodology we have developed since our early

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

137

Fig. 1. Overview of the initial mesh.

publications [7–10], and in particular since the BANC-I workshop, is a natural direction to pursuit for future CFD/CAA simulation methodologies. However, one fundamental difference between the mesh adaptivity strategy used by Vatsa et al. and the one we develop is that we base our mesh refinement on carefully derived a posteriori error estimates, which relate the residuals of the equations being solved, and the solution of an adjoint problem, to the discretization error in a functional of the computed solution. There are several aspects that we would like to highlight with this publication. The nose landing gear is a highly complex geometry, for which there are only a few CFD frameworks that currently could produce accurate unsteady flow simulations, mainly due to the limitations imposed by mesh size requirements. The main objective of this article is to compare our simulation results with the experiments performed in the NASA Langleys closed-wall Basic Aerodynamic Research Tunnel (BART) [11]. Such a comparison was previously reported for the RLG geometry [1], the difference now being, besides a significant increase in geometry complexity, that several other measurement techniques were employed, which allows for a more quantitative comparison of the flow field. These techniques include, among others, PIV measurements and unsteady pressure coefficients. Moreover, we would like to present to the CFD/CAA community our work on mesh adaptivity methods based on a posteriori error estimation. To the best of our knowledge, the application of this technique to unsteady flow CFD/CAA is highly uncommon [12]. Finally, since this is a subject surrounded by controversy, and a one highly prone to be met with skepticism, we would like to emphasize that we obtain good agreement with experiments without resolving any turbulent boundary layers on the surface of the landing gear. This has been a crucial assumption in our previous work and the one that has enabled us to tackle highly complex problems with relative ease. In the following sections, we describe the mesh adaptivity strategy and present some of the details of the computational mesh (Section 2); briefly present our method and formulate the a posteriori error estimate that forms its core (Section 3); give an overview of the computational resources and tools used in the simulations (Section 4); and, finally, compare our results with experiments (Section 5). 2. Adaptive mesh generation 2.1. The algorithm A simple description of a G2 adaptive algorithm reads: 1. For the mesh T n with cells K: compute the primal problem and the (linearized) adjoint problem.  2. If K∈Tn EK < T OL then stop, else:

3. Mark the 10% of cells with highest error indicator, EK , for refinement. 4. Generate the refined mesh Tn+1 , and goto 1. Here, the primal problem is to solve the Navier–Stokes equations (NSE) and EK is the error indicator for each cell K, which we describe more precisely in Section 3. For now, it suffices to say that EK is a function of the residual of the NSE and the solution of a linearized adjoint problem to these equations. The formulation of the adjoint problem includes the definition of a target functional for the refinement, which enters the adjoint equations as a boundary condition or a volume source term, and which should be chosen according to the problem we are solving. In other words, although the algorithm is fully automated, and will produce “converged solutions” as it is finished, one needs to ask the right question in order to obtain the desired answer from the algorithm. In our previous publications [7–10,13], we mostly solved aerodynamic problems where the aerodynamic drag constituted a good choice of functional. For CAA applications, however, we introduced a new functional based on Lighthill’s analogy [14], which produces different mesh refinement patterns than the drag functional, and which we believe is more suitable for CAA [15]. Apart from a suitable choice of data for the adjoint problem, the only other input required from the user is an initial discretization of the geometry, T0 . Since our FE method is designed for tetrahedral meshes that do not require any special treatment of the near-wall layers (no need for boundary-layer mesh), the initial mesh can be easily created with any standard mesh generation tool. 2.2. Overview of the mesh and features of the adaptive algorithm The initial mesh is shown in Fig. 1. The surface mesh clearly illustrates how we usually create the initial mesh. On the one hand, the surface of the landing gear is accurately represented in its finest details. On the other hand, by examining the volume mesh, we clearly see that the mesh is dramatically coarser close to the wind tunnel walls, away from the landing gear, with cell sizes that are as large as the wheel diameter. This way, we leave the decision of where to refine the volume mesh to the adaptive algorithm, which will then produce a fine mesh based on the choice of output functional. Figs. 2 and Fig. 3 show the final mesh, which was obtained after 6 adaptive iterations. In Fig. 2 we see three cuts of the final volume mesh at different Z-coordinates, and it is clear that the refinement pattern is not the same at the different cut plane positions. Moreover, the three geometrical shapes (manually drawn) illustrate the approximate regions where the refinement is more concentrated. Interesting to notice is that the blue shape was intentionally drawn unsymmetrically, since there are more cells refined on the port side of the geometry, extending up to the wind tunnel wall, than on the

138

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

Fig. 2. Overview of the final mesh.

∂2T

Fig. 3. Snapshot of the refinement target, ∂ x ∂ixj , on the final mesh. i j

starboard side (the flow is from left to right). This is also a feature of the mean flow as captured by the PIV measurements: although a less careful observation of the landing gear would suggest symmetry about its longitudinal plane, the geometry is not really symmetric, which results in an asymmetric flow (see Section 5.2). This illustrates one of the main advantages of adaptivity, which is the

ability of “revealing” features of the flow that cannot be predicted a priori. Another striking example of how mesh adaptivity works may be seen in Fig. 3, where we superimposed the visualization of the refinement target (Lighthill’s source term, ∂ 2 Tij /∂ xi ∂ xj , see Section 3.4) and the volume mesh at two different cut planes. In the top figure, we see

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

traces of Lighthill’s source term stretching several wheel diameters downstream of the landing gear. There is certainly not much sound being generated by the quadrupoles all the way to the right of the figure; however, we cannot be sure a priori of how much a coarser mesh in the wake would affect the numerical solution upstream, close to the landing gear. Moreover, an even less predictable phenomenon is whether the quadrupoles in the wake near the landing gear produce any significant amount of noise or not. What the adaptive refinement indicates is that there is a considerable amount of quadrupole type fluctuations in the wake, in particular up to 4 wheel diameters, 4D, downstream of the landing gear. The role of quadrupoles in landing gear noise is a difficult problem still under investigation [16,17] and we would certainly like to study that in future stages of our work, in particular when we are able to propagate the aeroacoustic sources calculated with G2 to the far field by solving Lighthill’s equation. It is a customary “best practice” in CFD/CAA to ensure that the computed solution is independent of the mesh size. In other words, a “mesh study”, where the same problem is solved on different meshes, should accompany any simulation in order to guarantee that possible errors in the solution do not come from the discretization. In practice, this is very difficult to achieve, especially for time-dependent problems, where the cost of each simulation is prohibitively high. What is unique with the adaptive algorithm in G2 is that no “mesh study” is needed: once the stop criterion is met (in “step 2” of the algorithm above), the method should guarantee that further mesh refinement would not change the solution with respect to the output functional. As we describe in the following section, the a posteriori error estimates in G2 are based on the computation of a functional of the solution, which typically involve averaging over time and space of the solution of the NSE. In Section 5.1, we analyze some results that indicate how the solution varies with the discretization for the current problem.

3.1. The cG(1)cG(1) method

∇ · u = 0,

(2) where U¯ n = 12 (U n + U n−1 ) is piecewise constant in time over In , with the stabilizing term

SDnδ (U¯ n , P n ; v, q) ≡ and

(v, w) =

(δ1 (U¯ n · ∇U¯ n + ∇ Pn − f ), U¯ n · ∇v + ∇ q) + (δ2 ∇ · U¯ n , ∇ · v),

 3 

K

(3)

v · w dx,

(i j (v), i j (w)),

i, j=1

n−1 |2 h−2 )−1/2 and with the stabilization parameters δ1 = κ1 (k−2 n + |U δ2 = κ2 h|U n−1 |, where κ 1 and κ 2 are positive constants of unit size.

3.2. A posteriori error estimate for cG(1)cG(1) The a posteriori error estimate is based on the following theorem (for a detailed proof, see Chapter 30 in Hoffman and Johnson [18]): Theorem 1. If Uˆ = (U, P ) solves (2), uˆ = (u, p) is a weak NSE solution, and ϕˆ = (ϕ , θ ) solves an associated adjoint problem with data M( · ), then we have the following a posteriori error estimate for the target functional M(Uˆ ) with respect to the reference functional M(uˆ):

|M(uˆ) − M(Uˆ )| ≤ +

  In K∈T n

 N   n=1

In

|R1 (Uˆ )|K · ω1 dt

K∈Tn

|R2 (U )|K ω2 dt +

  In K∈T n

 ˆ ϕ) |SDδ (U; ˆ K |dt n

=:



EK

K∈Tn

with

As the basic model for incompressible Newtonian fluid flow, we consider the NSE with constant kinematic viscosity ν > 0, enclosed in  ⊂ R3 , with boundary  , over a time interval I = (0, T ]:

u(x, 0) = u0 (x),

n ¯n ¯n ¯n ((U n − U n−1 )k−1 n + U · ∇ U , v) + (2ν(U ), (v)) − (P , ∇ · v) n n n n + (∇ · U¯ , q) + SDδ (U¯ , P ; v, q) = ( f, v), ∀vˆ = (v, q) ∈ V0 × W,

((v), (w)) =

The simulation methodology uses a finite element method with piecewise linear approximation in space and time, with a numerical stabilization in the form of a weighted least squares method based on the residual. The NSE are discretized directly, without applying any filter. Thus, the method does not produce a Large Eddy Simulation (LES) filtered solution, but is instead an approximation of a weak solution, satisfying the weak form of the NSE. For this method, we can derive a posteriori error estimates of quantities of interest with respect to a weak solution, which form the basis for our adaptive mesh refinement algorithm. The a posteriori error estimates are based on the solution of an associated adjoint problem, similar to an optimal control problem. We refer to this method as General Galerkin (G2) or Direct Finite Element Simulation (DFS). It is a time-resolved method, being related to an implicit LES, as no explicit subgrid model is used. G2 is validated for a number of standard benchmark problems in the literature [7–10]. In the following subsections, we describe the method and present the modifications proposed for CAA.

u˙ + (u · ∇)u + ∇ p − 2ν∇ · (u) = f,

tensor, with the strain rate tensor i j (u) = 1/2(∂ ui /∂ x j + ∂ u j /∂ xi ), and δ ij the Kronecker delta function. The relative importance of viscous and inertial effects in the flow is determined by the Reynolds number Re = UL/ν, where U and L are characteristic velocity and length scales. The cG(1)cG(1) method is based on the continuous Galerkin method cG(1) in space and time. With cG(1) in time, the trial functions are continuous, piecewise linear and the test functions piecewise constant. cG(1) in space corresponds to both test functions and trial functions being continuous, piecewise linear. Let 0 = t0 < t1 < · · · < tN = T be a sequence of discrete time steps, with associated time intervals In = (tn−1 , tn ) of length kn = tn − tn−1 , and let W ⊂ H1 () be a finite element space consisting of continuous, piecewise linear functions on a tetrahedral mesh T = {K } of mesh size h(x), with Ww the functions v ∈ W satisfying the Dirichlet boundary condition v| = w. We seek Uˆ = (U, P ), continuous piecewise linear in space and time, and the cG(1)cG(1) method for NSE with homogeneous Dirichlet (no slip) boundary conditions reads: for n = 1, . . . , N, find (Un , Pn ) ≡ (U(tn ), P(tn )), with Un ∈ V0 ≡ [W0 ]3 and Pn ∈ W, such that:

K∈T

3. The General Galerkin method (G2) and its a posteriori error estimate

139

(x, t ) ∈  × I, (x, t ) ∈  × I, x ∈ ,

(1)

with u(x, t) the velocity vector, p(x, t) the pressure, u0 (x) initial data and f(x, t) a body force. Moreover, σi j = −νi j (u) + pδi j is the stress

R1 (Uˆ ) = U˙ + (U · ∇)U + ∇ P − 2ν∇ · (u) − f, R2 (U ) =

∇ · U,

(4)

where SDnδ (·; ·)K

is a local version of the stabilization form (3), and the stability weights are given by

ω1 = C1 hK |∇ϕ|K , ω2 = C2 hK |∇θ |K , where hK is the diameter of element K in the mesh Tk , and C1, 2 represent interpolation constants. Moreover, |w|K ≡ ( w1 K , w2 K ,

140

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

w3 K ), with w K = (w, w)1/2 , and the dot denotes the scalar prodK uct in

R3 .

For simplicity, it is assumed that the time derivatives of the dual variables φˆ = (φ , θ ) can be bounded by their spatial derivatives. Given Theorem 1, we are able to formulate the adaptive algorithm more precisely, as follows. The error indicator, EK , is a function of the residual of the NSE and of the gradients of the solution of a linearized adjoint problem to these equations (a detailed formulation of the adjoint problem is given in Chapter 14 in Hoffman and Johnson [18]). Thus, on a given mesh, we must first solve the NSE to compute the residuals, R1 (Uˆ ) and R2 (U), and then a linearized adjoint problem to compute the weights multiplying the residuals, ω1 and ω1 . With that  information, we are able to compute K∈Tn EK and check it against the given stop criterion. This procedure of solving the forward and backward problems for the NSE is closely related to an optimization loop and can be understood as the problem of finding the “optimal mesh” for a given geometry and boundary conditions, i.e., the mesh with the least possible degrees of freedom for computing the output functional, M(uˆ), with a given degree of accuracy. 3.3. Free-slip velocity as a model of turbulent boundary layers In the experimental campaign at the NASA Langley Research Center, boundary layers were tripped to turbulent on the shock strut, and they were also expected to be either turbulent or transitional over the remaining components of the landing gear model due to surface roughness [11]. Thus, for the present simulations, we used a free-slip boundary condition for the velocity, which should be a good model for landing gear simulations at flight conditions, where Re is typically much larger than in a wind tunnel experiment. This choice of boundary condition is motivated by our previous work on external, high Re turbulent flow [1,19,20]. It is worth noting here that, as opposed to the circular cylinder case in [19], most of the flow separation from the landing gear surface is triggered by sharp edges in the geometry (e.g., see Fig. 6). This diminishes the role of boundary layers in the prediction of separation for the present geometry. 3.4. Lighthill’s analogy as refinement target In [15], we introduced a new refinement strategy especially derived for CAA applications. The new strategy is based on Lighthill’s analogy, proposed by Lighthill in his seminal paper from 1952 [14]. In that paper, Lighthill defined what today is known as the Lighthill tensor:

T ≈ ρ uuT .

(5)

Eq. (5) shows the simplified expression for the tensor, where it is assumed that the viscous stresses are negligible in the production of sound, and also that no heat is generated, nor advected, within the domain, which is true in our simulations. Lighthill derived a wave equation directly from the NSE without any approximations [14], where the second derivative of the tensor defined in Eq. (5) appears as a source term. We, therefore, introduced a new functional to be used to drive the adaptive algorithm, which is basically an averaging of the second derivative of the Lighthill tensor in the whole domain:

M(uˆ) =

1 |I|

  I



∂ 2 Ti j /∂ xi ∂ x j dx dt.

(6)

Eq. 6 is a volume integral of the source term of Lighthill’s wave equation in its differential form over the entire domain, and should not be interpreted as an integral solution of Lighthill’s equation using Green’s functions techniques (as Lighthill himself presented in his paper from 1952, and that Curle further developed to include surface

terms in 1955). The goal here is to quantify the error in all acoustic sources with Eq. (6), even those acting on solid walls. One could argue that, for low Mach number flows, it would be less expensive to include only dipole sources in the refinement strategy, e.g., by limiting the integral in (6) to a smaller domain near solid walls. However, this is not as general as including the entire volume – there are certainly examples where it is not known a priori whether dipoles or quadrupoles are the dominating sound source, or even other examples, where neither of them dominate, but both are relevant. We used this definition of M(uˆ) for computing the error indicators described in Theorem 1. 4. Computational tools and resources The G2 method is implemented in the Unicorn solver, which uses the open source software framework FEniCS-HPC [21–24]. This framework is designed for the automated solution of partial differential equations on massively parallel architectures using FE methods and provides, among other features, automated evaluation of variational forms given their high-level description in mathematical notation, adjoint-based adaptive error control, implicit turbulence modeling by use of a stabilized FE method, and strong linear parallel scaling up to thousands of cores. For this work, we used the computer Lindgren at PDC – Center for High Performance Computing at KTH in Stockholm, Sweden. Lindgren is a Cray XE6 system with 1516 nodes, each containing a dual AMD Opteron 12-core “Magny-Cours” (2.1 GHz) processor (36,384 cores in total), and 32 GB DDR3 of memory per node. The interconnect technology is Cray Gemini, with a 3D Torus topology. It is possible to divide our computational work in two stages. First, the adaptive mesh refinement is performed and a series of primal and adjoint problems are solved, with the mesh being refined as the algorithm advances (see Section 2); then, the finest mesh is used for sampling data for postprocessing where the simulation should be long enough to provide reliable statistical quantities. For the adaptive refinement of the landing gear mesh, we needed ca. 900,000 core hours on a varying number of processors; for the finest mesh, we used ca. 600,000 core hours on 2304 cores for sampling. If we break down the above figures to investigate the amount of core hours required for each adaptive iteration, we obtain (in thousands of core hours, from coarsest to finest mesh): 27, 41, 76, 120, 197 and 4391 . Thus, the amount of CPU resources required for each iteration grows exponentially as the mesh is refined (with an exponent of ca. 1.4). For instance, if the next mesh in the adaptive algorithm would contain ca. 35 million cells, it would require approximately 690,000 core hours for a complete iteration. A normalized parameter that should indicate the general performance of the code is the CPU-time required to run one time-step per cell. We computed this parameter for all meshes in our computations, obtaining (in microseconds and in ascending order of number of cells): 78, 87, 94, 101, 110 and 123. These numbers show that the CPU-time required per time-step per cell also increases as the mesh is refined. In order to compare the efficiency of our code with another competing code, we computed the same normalized parameter in the simulations by Vatsa et al. [25]. For those computations, the CPUtime required to run one time-step per cell is 216 microseconds. This is slower than our simulations by a factor of almost 2; however, one should have in mind that Vatsa et al. performed compressible simulations, which are, in general, more expensive than incompressible ones (due to larger number of variables and equations being computed, difficulties in time-stepping for low Mach number flows, etc). 1 This is the time required for solving primal and adjoint problems only, the sampling window was not included here.

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

199

Lighthill source term [dB]

198

197

196

195

194

193

192 0

0.5

1

1.5

number of cells

2

2.5 7

x 10

∂2T

Fig. 4. L2-norm of the refinement target, ∂ x ∂ixj , as a function of the number of mesh i j cells in the mesh.

During our final simulation, where statistics were sampled, a total of 1.4 TB of data was generated at a sample rate of 32 kHz, adding extra difficulty to the visualization and postprocessing work. Finally, the postprocessing of data for the PIV comparisons was performed with Saaz [26]. 5. Results The simulations in this publication were setup to reproduce similar conditions as in the experimental campaign at NASA [11], with an inlet velocity of 56.6 m/s. For the finest mesh, a time step of 6.4997 · 10−7 s was used and the total simulation time was 0.28 s, where 0.08 s was considered as start-up phase and was discarded. In the following paragraphs, we present detailed comparisons of the computational results with the measurements performed in the BART tunnel at the NASA Langley Research Center. But first, we discuss mesh convergence, which is an important result of our adaptive computations. 5.1. Mesh convergence As described in Section 3.4, the mesh refinement algorithm used the source term of Lighthill’s wave equation, ∂ 2 Tij /∂ xi ∂ xj , as the refinement target. Fig. 4 shows the L2 -norm of ∂ 2 Tij /∂ xi ∂ xj as a function of the number of cells in the mesh. In order to increase the physical meaning of the quantities in the figure, we chose to convert ∂ 2 Tij /∂ xi ∂ xj to decibels; to achieve that, we first multiplied 2 to obtain a quantity that has units of density, ∂ 2 Tij /∂ xi ∂ xj by D2 /U∞ where D is the wheel diameter and U∞ is the inlet flow speed; then, we used an analogy to the relationship between acoustic quantities valid for isentropic flow, p = ρ c0 , where ρ is the acoustic density, p is the acoustic pressure, and c0 is the speed of sound, to obtain a quantity that has units of pressure; the conversion to decibels is then straightforward. As expected, the L2 -norm of ∂ 2 Tij /∂ xi ∂ xj increases with the number of cells, which means that more of the energetic turbulent scales are resolved as the mesh is refined. Despite the fact that we expect that this functional, which is the second derivative of a non-linear turbulent quantity, will continue to increase with further refinements, Fig. 4 indicates that the growth becomes milder as the number of cells increases, which is a desirable behavior. Due to limitations of both computational resources and time, we decided to stop the refinement algorithm before clear convergence was achieved.

141

In future work, we will investigate convergence with respect to the solution of Lighthill’s wave equation (or any other acoustic analogy equation) and quantify the effect of mesh refinement on the sound pressure levels at the far field. This way, one would be sure that the aerodynamic sound emitted by the configuration, which not only depends on the intensity of the sources, but also on the phases and retarded times of the different dipole and quadrupole sources inside the wake, is being correctly predicted. A more qualitative assessment of convergence can be performed by observing contours of ∂ 2 Tij /∂ xi ∂ xj on different meshes. Fig. 5 shows snapshots of the solution in a plane passing through the center of the domain for the initial mesh and for the mesh after 2, 4 and 6 refinements. The difference is clear between the 2 finest and the 2 coarsest meshes in the figure. Besides an increase of intensity of Lighthill’s source term in a wider part of the wake, there is a clear improvement in the resolution of the shear layer produced by the wake of the wheels, close to the top wall of the wind tunnel. However, a visual comparison of the results on the 2 finest meshes (in the bottom half of Fig. 5) does not present any major discrepancies apart from an intense red spot close to the fuselage at the bottom of the wind tunnel of the finest mesh, which could be an artifact caused by the visualization of the instantaneous flow field. 5.2. Comparison with PIV measurements Fig. 6 brings an overview comparison with PIV plane measurements, where time-averaged contours of streamwise velocity, 2D turbulent kinetic energy, TKE2 , and Z-vorticity are shown; the quantities were plotted on two distinct planes, one passing through the side of the starboard wheel and the other through the wake of the landing gear door. There is good qualitative agreement between experiments and simulation. One remarkable feature of the comparison is that there is a clear asymmetry between port and starboard sides when observing the door plane, both in the streamwise velocity and the TKE fields. This is not an intuitive result; however, as discussed in Section 2.2, a closer look at the geometry reveals that the landing gear is indeed asymmetrical with respect to its longitudinal plane. Fig. 2 shows that the adaptive mesh refinement automatically captures this behavior, which would be a somewhat intricate task had we generated and refined the mesh manually. Another interesting feature of the comparison is that we seemingly get an overshoot of vorticity on the side of the starboard wheel; this is probably due to under resolution of the PIV measurement mesh as compared to the resolution of the computational mesh, and not an “overshoot”: the wheels constitute a strong source of unsteadiness for this configuration and we expect the mesh refinement to capture this behavior, producing a fine mesh in that area. Fig. 7 shows the mesh size projected on a plane passing through the center of the wheel in a region close to where the PIV measurement was performed. The PIV resolution is 1.3 mm [11] and the contours in Fig. 7 show that the mesh is indeed smaller than 1.3 mm over several cells coinciding with the region of strongest vorticity. 5.3. Steady pressure probes Several static pressure probes were placed on the landing gear geometry, both on the surface of the starboard wheel and on the surface of the landing gear door. Fig. 8 shows the time-averaged pressure coefficient, Cp , distribution on the starboard wheel plotted against the position angle with respect to the free stream. Noteworthy here is that, despite the fact that we do not resolve any 2 The kinetic energy on a measurement plane, i.e., T KE = 1/2((u )2 + (v )2 ), where u and v are the fluctuations of the velocity components on the measurement planes in the streamwise and spanwise directions, respectively.

142

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

∂2T

Fig. 5. Snapshots of ∂ x ∂ixj on different meshes; from top to bottom: initial mesh, and mesh after 2, 4 and 6 refinements. i j

Fig. 6. PIV, location of measurement planes; time-averaged fields of streamwise velocity, 2D turbulent kinetic energy and vorticity. In all figures: left, simulation, right, measurements [11].

boundary layer on the entire geometry, we still get a good match with the experimental static pressure distribution on the wheel. Moreover, we do not expect that the minor differences encountered in this comparison will have a great impact upon the overall sound pressure level at the far field and we consider this result satisfactory, given the ease and low cost that our boundary layer model provides us with during both the meshing process and the computations, respectively. Fig. 9 shows the Cp distribution on the surface of the door. The sensors were placed in rows that were numbered from 1 to 9 as to

the left in the figure. Here, the agreement is less satisfactory than on the wheel surface and some of the sensors are off by a factor 2. The result is, nonetheless, still quantitatively good for most of the rows, especially rows 4 through 8, and the general trend is well captured on the other rows (except for row 1). 5.4. Unsteady pressure probes A better indication of the near field sound sources is obtained by examining unsteady pressure signals. The simulation results were

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

143

Fig. 7. Mesh size distribution near the region where PIV measurements were performed. The contour shows that the computational mesh resolution is finer than the PIV resolution of 1.3 mm [11] on several cells.

sampled during 0.2 s (after a start-up phase of 0.08 s) at a frequency of 32 kHz, which results in a maximum frequency resolution of 16 kHz due to Nyquist’s theorem. For the calculation of the Power Spectral Density (PSD) we divided the time series into 3 hamming windows with 50% overlap. Fig. 10 shows the PSD of simulation and measurement sensors and the comparison is of varying quality. The worst result was obtained at channel 3, where the simulation fails to capture the pressure levels obtained in BART. Moreover, at channel 12 the simulation overpredicts the level for frequencies up to 1, 000 Hz. Another marked difference is the shape of the curves for frequencies above 2, 000 Hz: for frequencies higher than this, the simulation curves decay much faster than the measured ones, which indicates lack of mesh resolution near these sensors. However, at channels 5, 7 and 13, the computed signals decay in a similar fashion as the measured ones up to 10,000 Hz; even at channel 15, the decay is not as accentuated for higher frequencies as in the remaining channels. To our understanding, this could be a consequence of the automated mesh refinement strategy used in our computational methodology. For some regions, where more sound is being generated, or at least where the magnitude of Lighthill’s source term is more sensitive to the flow features, we expect to obtain a finer mesh, and consequently, a better resolution of the turbulent fluctuations. This could explain why the turbulent signals are better resolved in the frequency domains for some channels, while for other ones, only the lower frequencies signal are resolved. Another way of quantifying the sound sources is to compute the root mean square of the pressure, prms , by integrating the PSD. The results for each channel are shown in Table 1, where we compare the Sound Pressure Level (SPL) of measurements and simulation. For some channels, the level is reasonably well captured, namely for channels 4, 5, 8 and 15: they all lie within 2 dB from the measured level. The error in the remaining channels could depend on several factors, e.g., the separation from the geometry upstream some of the

Fig. 9. Static pressure probes, door; locations, right, and static pressure coefficient, left. Table 1 Comparison of SPL for each channel (numbered as in Fig. 10). channel #

SPLexp (dB)

SPLsim (dB)

SPLexp − SPLsim (dB)

3 4 5 7 8 10 12 13 15

144 137 143 138 141 142 137 136 145

136 137 141 144 141 138 143 139 147

8 0 2 −6 0 4 −6 −3 −2

channels could be difficult to predict and the exact spot for turbulence impingement is not correctly captured in the computation (as in the case of channel 3, where the simulation result undershoots the measured value by an order of magnitude). Another plausible explanation is the short sampling time of simulation signals when compared to experiments, which deteriorates the capacity of the simulation to predict the exact location of flow features that have longer time scales. In order to address and understand some of these differences, we performed a “postprocessing artifice”, in which we explored one of the biggest advantages of simulations over experiments, which is the prompt availability of the complete flow field in each point of the mesh over the entire simulation time. Based on that and on the idea that we do not expect computational solutions to have a pointwise agreement with measurements, at least not for the highly complex, unsteady and non-linear phenomena of high Re separation and turbulence impingement, we collected unsteady pressure signals in a small

Fig. 8. Static pressure probes, starboard wheel; locations, right (several probes were omitted for simplicity), and static pressure coefficient, left.

144

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

channel #3

channel #4

2

2

−2

10

2

PSD [Pa /Hz]

2 PSD [Pa /Hz]

0

10

0

10

−2

10

3

10

10

4

10

exp sim 2

10

3

10

10

4

10

channel #7

channel #8

channel #10

2

−2

2

2 0

10

−2

3

10

4

10

exp sim 2

10

3

10

10

4

2

10

3

10

10

freq [Hz]

freq [Hz]

channel #12

channel #13

channel #15

2

2

−2

10

2

PSD [Pa /Hz]

2 PSD [Pa /Hz]

2

0

10

0

10

−2

3

10

freq [Hz]

4

10

0

10

−2

10

exp sim

4

10

2

10

2

exp sim

freq [Hz]

10

10

0

10

−2

10

exp sim

4

10

10

PSD [Pa /Hz]

2 PSD [Pa /Hz]

2

0

10

10

10

freq [Hz]

10

2

3

10

freq [Hz]

2

10

exp sim 2

freq [Hz]

10

10

0

10

−2

10

exp sim 2

PSD [Pa /Hz]

2

10

2

PSD [Pa /Hz]

10

PSD [Pa /Hz]

channel #5

10

exp sim 2

10

3

10

4

10

exp sim 2

10

freq [Hz]

Fig. 10. Unsteady pressure probes, locations (top) and power spectral density (bottom).

3

10

freq [Hz]

4

10

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

145

2

10

0

2

PSD [Pa /Hz]

10

−2

10

−4

10

channel #3 sim channel #4 sim channel #3 exp channel #4 exp

−6

10

2

3

10

4

10

10

freq [Hz]

(b) PSD, superimposed. (a) Sampling areas.

2

10

0

2

PSD [Pa /Hz]

10

−2

10

−4

10

channel #3 sim channel #4 sim channel #3 exp channel #4 exp

−6

10

2

3

10

10

4

10

freq [Hz]

(c) PSD, average. Fig. 11. Unsteady pressure for channels 3 and 4. The blue and red curves were generated by computing the PSD for ca. 700 points inside an area close to channels 3 and 4, respectively.

Fig. 12. Volume rendering of the λ2 vorticity criterion.

area, rather than a point, near channels 3 and 4 (see Fig. 11(a)). After performing a PSD for several points inside those areas (ca. 700 points for each channel), we plotted the results in Fig. 11(b) and (c). We were able to recover the asymmetry shown in the experimental results for channels 3 and 4, channel 3 being the one subjected to a stronger un-

steadiness in the BART measurements. Here we show that, although a pointwise agreement is extremely difficult to achieve for such flows, a general quantitative agreement is still possible and, most likely, integral quantities such as drag and lift forces, and overall sound pressure levels are computable to a desired tolerance.

146

R. Vilela de Abreu et al. / Computers and Fluids 124 (2016) 136–146

The question that remains to be answered is how these local differences would influence the sound prediction at the far field, which is a subject for future investigations. 5.5. Turbulence visualization Fig. 12 shows a volume rendering of the λ2 vorticity criterion on the finest mesh. The objective of this visualization is to show the amount of resolved turbulence structures in the wake and their extent to several wheel diameters downstream of the landing gear, especially near the fuselage geometry. This result differs significantly from the results of most of the contributors to the BANC-II workshop [4,5], due to construction of meshes that resolve the flow structures in the vicinity of the landing gear, but not so strongly in the turbulent wake away from the geometry. As mentioned earlier, whether this has an impact on the far field sound prediction or not is a subject for future work. 6. Conclusions This paper presented simulations performed using the stabilized FE method G2 – which uses fully unstructured meshes, simple walllayer models, and adaptive mesh refinement to solve flow problems in complex geometries –, compared with wind tunnel measurements performed at the NASA Langley Research Center. The main goals were to elucidate some of the principal characteristics of G2, highlight the importance of using adaptive mesh generation for predicting the flow around complex geometries, as well as to validate the simulation results by direct comparisons with measurements. To that end, we described in detail the main features of the adaptive mesh generation and presented results that indicate a trend of mesh convergence for the chosen refinement target, which for this work was the source term of Lighthill’s equation, ∂ 2 Tij /∂ xi ∂ xj . Moreover, the comparisons between simulation results and experiments – which include PIV measurements, and steady and unsteady pressure probes – good quality in our predictions, and build a solid ground for further investigations of flow past geometries of similar or greater complexity. Moreover, we investigated the presence of quadrupole sources with the help of the solution computed on several meshes, which is a natural output of the refinement algorithm. No hard conclusions could be drawn, except that, initially, i.e., on coarser meshes, the refinement is concentrated on the dipole sources, whereas for finer meshes, there is a shift of behavior and the refinement seems to account for the quadrupole sources in the wake too. The next step in the pursuit of a complete adaptive framework for CAA applications is to solve the problem of the propagation of sound sources to the far field, using our time-resolved solutions as the source term in an acoustic analogy. As a final concluding remark, we would like to point out that the G2 method, when applied to high Re flows, is a parameter free model and, hence, constitutes a powerful computational tool. Acknowledgments The authors would like to acknowledge the financial support from the Swedish Foundation for Strategic Research, the European Research Council - ERC Starting Grant “UNICON”, proposal number 202984. - ERC Proof of Concept Grant “ADAPTIVE”, proposal number 665630 and the Swedish Research Council. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC – Center for High-Performance Computing.

The initial volume mesh was generated with ANSA from Beta-CAE Systems S. A., who generously provided an academic license for this project. Finally, we are very grateful to Mehdi R. Khorrami from the NASA Langley Research Center for organizing Category 4 of the BANCworkshop series and for generously providing the experimental data. References [1] Vilela de Abreu R, Jansson N, Hoffman J. Adaptive computation of aeroacoustic sources for a rudimentary landing gear. Int J Numer Meth Fluids 2014;74(6):406– 21. doi:10.1002/fld.3856. [2] Spalart PR, Mejia K. Analysis of experimental and numerical studies of the rudimentary landing gear. In: Proceedings for 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, Florida; 2011. [3] Schumann U. Subgrid-scale model for finite difference simulation of turbulent flows in plane channels and annuli. J Comput Phys 1975;18:376–404. [4] Casalino D. Pdcc-nlg. In: Proceedings for the Second Benchmark problems for Airframe Noise Computations (BANC-II), Colorado Springs, CO, 7-8; June 2012. [5] Ueno Y. Results of ddes with layered grid (using khis modified “cflow”). In: Proceedings for the Second Benchmark problems for Airframe Noise Computations (BANC-II), Colorado Springs, CO, 7-8; June 2012. [6] Vatsa VN, Khorrami M, Lockard D. Application of adaptive gridding in fun3d for simulation of flow over a nose landing gear. In: Proceedings for the Second Benchmark problems for Airframe Noise Computations (BANC-II), Colorado Springs, CO, 7-8; June 2012. [7] Hoffman J. Computation of mean drag for bluff body problems using adaptive dns/les. SIAM J Sci Comput 2005;27(1):184–207. [8] Hoffman J. Adaptive simulation of the turbulent flow past a sphere. J Fluid Mech 2006;568:77–88. [9] Hoffman J, Johnson C. A new approach to computational turbulence modeling. Comput Methods Appl Mech Engrg 2006;195:2865–80. [10] Hoffman J. Efficient computation of mean drag for the subcritical flow past a circular cylinder using general Galerkin g2. Int J Numer Meth Fluids 2009;59(11):1241–58. [11] Neuhart DH, Khorrami MR, Choudhari MM. Aerodynamics of a gulfstream g550 nose landing gear model. In: Proceedings for the 15th AIAA/CEAS Aeroacoustics Conference (30th AIAA Aeroacoustics Conference), Miami, FL; May 2009. [12] Fidkowski KJ, Darmofal DL. Review of output-based error estimation and mesh adaptation in computational fluid dynamics. AIAA Journal 2011;49(4):673–94. doi:10.2514/1.J050073. [13] Jansson N, Hoffman J, Nazarov M. Adaptive simulation of turbulent flow past a full car model. Proc SC’11 State of the practice reports 2011. [14] Lighthill MJ. On sound generated aerodynamically. Proc R Soc Lond A 1952;211:564–87. [15] Vilela de Abreu R, Jansson N, Hoffman J. Adaptive computation of aeroacoustic sources for a rudimentary landing gear using Lighthill’s analogy. In: Proceedings for the 17th AIAA/CEAS Aeroacoustics Conference (32nd AIAA Aeroacoustics Conference), Portland, OR; June 2011. [16] Spalart PR, Shur ML, Strelets MK, Travin AK. Initial noise predictions for rudimentary landing gear. J Sound and Vib 2011;330(17):4180–95. doi:10.1016/j.jsv.2011. 03.012. http://www.sciencedirect.com/science/article/pii/S0022460X11001921 [17] Spalart PR, Shur ML, Strelets MK, Travin AK. Sensitivity of landing-gear noise predictions by large-eddy simulation to numerics and resolution. In: Proceedings for the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Nashville, Tennessee; January 2012. [18] Hoffman J, Johnson C. Computational Turbulent Incompressible Flow. Applied Mathematics: Body and Soul, vol. 4. Springer; 2007. [19] Hoffman J, Jansson N. A Computational Study of Turbulent Flow Separation for a Circular Cylinder Using Skin Friction Boundary Conditions. Ercoftac, Vol. 16. Springer; 2010. [20] Hoffman J, Johnson C. Resolution of d’alembert’s paradox. J Math Fluid Mech December 2008;10.Published online first at www.springerlink.com [21] Hoffman J, Jansson J, Vilela de Abreu R, Degirmenci NC, Jansson N, Müller K, et al. Unicorn: parallel adaptive finite element simulation of turbulent flow and fluidstructure interaction for deforming domains and complex geometry. Comput Fluids 2012. [22] Hoffman J, Jansson J, Nazarov M, Jansson N. Unicorn: a unified continuum mechanics solver. In: Automated Scientific Computing. Springer 2011. [23] Hoffman J, Jansson J, Jansson N, Johnson C, Vilela de Abreu R. Turbulent flow and fluid-structure interaction. In: Automated Solution of Differential Equations by the Finite Element Method; ch. 28. Springer 2012. [24] FEniCS. Fenics project 2003. fenicsproject.org [25] Vatsa VN, Lockard DP, Khorrami MR. Fun3d solutions for nose landing gear. In: Proceedings for the First Benchmark problems for Airframe Noise Computations (BANC-I), Stockholm, Sweden; June 2010. [26] King A, Arobone E, Baden SB, Sarkar S. The saaz framework for turbulent flow queries. In: Proceedings for the 7th IEEE International Conference on e-Science, Stockholm, Sweden; December 2011.