Provably stable overset grid methods for computational aeroacoustics

Journal of Sound and Vibration 330 (2011) 4161–4179

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Provably stable overset grid methods for computational aeroacoustics$ Daniel J. Bodony a,, George Zagaris b, Adam Reichert c, Qi Zhang d a

Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, United States Computational Science and Engineering, United States Department of Computer Science, United States d Department of Aerospace Engineering, United States b c

a r t i c l e i n f o

abstract

Article history: Accepted 7 February 2011 The peer review of this article was organised by the Guest Editor Available online 25 February 2011

The simulation of sound generating flows in complex geometries requires accurate numerical methods that are non-dissipative and stable, and well-posed boundary conditions. A structured mesh approach is often desired for a higher-order discretization that better uses the provided grids, but at the expense of complex geometry capabilities relative to techniques for unstructured grids. One solution is to use an overset mesh-based discretization where locally structured meshes are globally assembled in an unstructured manner. This article discusses recent advancements in overset methods, also called Chimera methods, concerning boundary conditions, parallel methods for overset grid management, and stable and accurate interpolation between the grids. Several examples are given, some of which include moving grids. & 2011 Elsevier Ltd. All rights reserved.

1. Introduction Early studies of sound generation and propagation fall into two different classes: (1) fundamental studies of geometrically simple flows, and (2) engineering studies of geometrically complex flows. While there are likely many reasons for this division, the two most apparent concern the state of numerical methods and boundary conditions at the time of the simulations. In the jet noise community, the two-dimensional shear layer study of Colonius et al. [1] and the heated, supersonic three-dimensional jet of Estivales and Gamet [2] are examples of the first and second category, respectively. Continued demands for reduced acoustic radiation from existing engineering systems motivate bridging the geometric complexity gap between the two categories, and the series of increasingly complex computational aeroacoustic benchmarks [3–6] echo this trend along with the need to ensure that new algorithms are well tested and verified. It follows, then, that simulations are focusing on models with increased geometrical complexity but are still achievable on current computational resources. Noise prediction from commercial aircraft landing gear [7] and from axisymmetric jet nozzles with chevrons [8] are two pacing examples. These examples also illustrate the points of view being taken regarding how one should approach complex geometries. In the first example an existing code using an unstructured discretization was modified to reduce the numerical dissipation inherent in upwind-based schemes while still retaining stable solutions. In the second example a multi-block structured code used neighboring overlapping blocks and data were passed between them using interpolation; the underlying numerical scheme was non-dissipative with filtering for regularization. Unstructured methods are also seeing significant improvement [9–12].

$

This article is an expanded, peer-reviewed version of an article previously published in Procedia Engineering, 6C, 2010, pp. 234–243.

 Corresponding author. Tel.: +1 217 244 3844.

E-mail address: [email protected] (D.J. Bodony). 0022-460X/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2011.02.010

4162

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

Background Grid

Rotor Grid (feature grid)

Stator Grid (feature grid)

Fig. 1. NASA Glenn Source Diagnostic Test Fan (a) and its component grids taken at the 75 percent radius (b). Note that the interblade passages are not included in (b).

Methods based on structured grids are more easily able to use higher-order discretizations which reduce the number of grid points required to resolve a flow feature of interest and, hence, the required number of points necessary in an unsteady three-dimensional simulation by an order of magnitude or more compared to second-order methods [13]. Single block implementations, where only one grid is used, have limited capability to model complex geometries even when grid blanking and curvilinear coordinates are used. Multi-block implementations, where more than one grid is used, have improved geometry modeling capabilities at the expense of determining the communication between the blocks. Traditional multiblock codes require that the intergrid interfaces match. Overset-based multi-block codes, on the other hand, do not require interface matching and are capable of more complex geometries. They also offer advantages for moving grid problems that are not possible with traditional multi-block codes in that the grids may be in relative motion without restriction. The overset method, also termed the Chimera1 method, can be considered a locally structured—globally unstructured approach where multiple, independent grids are assembled and data are passed between them. The earliest known numerical application was for the solution of elliptic partial differential equations [14], though early fluid simulations using the approach followed soon thereafter [15,16]. This article describes recent developments in improving numerical algorithms for overset methods. 1.1. Issues relevant to complex geometries on overset structured meshes Several issues arise when using overset structured meshes. As with other high-order methods, including those with single grids, boundary conditions are crucial and remain an active area of research. Further, overset codes present unique challenges to ensure that their behavior is predictable and accurate. Because the grids need not match one has to determine how the individual grids overlap. One must allow data from one grid to be passed to another in a meaningful and stable manner for general nonlinear problems. If the grids are in relative motion then the inter-grid mappings are functions of time and the communication links must be updated as part of the solution. Thus there are the following issues that affect codes based on overset meshes: (1) well-posed boundary conditions and their accuracy, (2) management of overset grids, and (3) stable and accurate exchange of data between grids. 1 The term ‘Chimera’ is derived from the Greek mythological creature of the same name described as ‘a thing of immortal make, not human, lion-fronted and snake behind, a goat in the middle’ by Homer in the Illiad.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4163

The latter two topics are discussed below while a detailed presentation concerning the first item can be found in Ref. [17]. Examples are given throughout the text. The application of interest motivating this work is the Source Diagnostic Test (SDT) Fan at NASA Glenn shown in Fig. 1(a) with the objective of developing a numerical methodology suitable for broadband fan noise investigations that was more general than sliding-mesh or chorochronic methods [18] . Specifically, the objective was to develop a numerical approach capable of handling the non-periodic fluctuations of broadband fan noise which arise, in part, from the turbulent wake of the rotors and stators interacting with the downstream blade rows. The test case was simplified to a twodimensional ‘representative problem’ using the SDT geometry at the 75 percent radius for a single stage, shown in Fig. 1(b), at approach conditions (7809 RPM at 62 percent power). (Further details of the simulations are given below.) The interblade passages were neglected to further simplify the geometry, but the rotor moves at the proper rate and its attached grid moves with it while the stator and background grid are static. 2. Management of overset grids Determination of the grid overlap and establishing the inter-grid communication links is a unique aspect of Chimerabased methods and has been the focus of a number of studies, primarily for static grids, and several grid assembly codes currently are available. The flow solvers Overflow [19] and Beggar [20] have internal overset grid assembly capability; however, this capability is not easily accessible to third party flow-solvers. PEGASUS [21,22] has been incorporated within a loosely integrated simulation framework and applied to moving grid simulations [23,24], and information can be shared between the fluid solver and PEGASUS, either through the filesystem or in-core. For both structured and unstructured grids SUGGAR and DiRTlib have been employed to solve a wide range of moving grid problems where the grid motion is known a priori [25], or prescribed by a 6DOF library as demonstrated in Ref. [26]. SUGGAR++, the successor of SUGGAR currently being developed, utilizes a direct hole-cutting method [27] and improves the automation, robustness and efficiency of SUGGAR [28]. Ogen, the grid generator for the Overture package, has recently added moving grid capability [29]. The overset grid approach has inherent coarse-level parallelism since the domain of interest, O, is defined by a S composite grid system, SðOÞ ¼ N1 i ¼ 0 ðGi Þ, consisting of N independent grids, Gi , with N Z 2. For moving grid simulations, the grid assembly process is iteratively invoked by the flow solver at each timestep to transfer the solution from one grid to another. Moreover, in a parallel setting, each grid, Gi 2 SðOÞ, is distributed to Mi processors by first partitioning into Mi smaller sub-grids, g(i,j), where i 2 f0,N1g denotes the global grid ID and j 2 f0,Mi 1g denotes the partition ID. While tools like SUGGAR and PEGASUS provide parallel overset grid assembly capability, their use within a moving body simulation framework requires a merge-and-repartition step at each time-step, where the sub-grids g(i,j) are merged to reconstruct the original grid Gi prior to determining the overset parameters. Algorithms for handling an already distributed composite grid system are needed to alleviate this overhead. A unique aspect of this work is the development of an algorithm for distributed overset grid assembly [30] capable of handling arbitrarily partitioned and distributed grids. Our algorithms are in library-form which enables an in-core and tightly coupled integration with the flow-solver. Generally, the grid assembly operation consists of four main steps: 1. Establish inter-grid communication: determines which grids need to pass information. 2. Hole-cutting: determines that part of the background which is inside a solid object and should not be used. 3. Donor–receiver pair search: determines the relationship between nodes needing data on one grid and the computational cell on another grid which will provide them. 4. Interpolation: computes the data from the donor cell and passes it to the receiving node. These steps are shown in Fig. 2 using the motivating example of the SDT fan. Each step is described in more detail below. 2.1. Establish inter-grid communication The first and most important step in the algorithm is to establish the inter-grid communication information. This information identifies the overlapping grid regions within a given overlapping grid system S N ¼ fG0 ,G1 ,G2 , . . . ,GN1 g distributed to a number of processors. Each overlapping grid region, Rij , is approximated by an axis-aligned bounding box computed as the bounding box intersection of two colliding grids, Rij ¼ BðGi Þ \ BðGj Þ. Rij is used in subsequent steps of the grid assembly process to minimize the search space. The primary steps needed to establish inter-grid communication are outlined below: 1. Each process, Pi, computes the bounding boxes of each grid it ‘owns.’ The bounding boxes for the rotor and stator blades are sketched in the upper-left pane of Fig. 2. 2. Using collective communication, all of the grid bounding boxes are distributed to all processes. Each process stores the grid bounding boxes in a matrix B where each row number, i, corresponds to process ID and each column, j, corresponds to the bounding box of the jth grid owned by process Pi.

4164

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

Fig. 2. Schematic of steps involved in overset grid assembly.

3. Each process Pi then performs a pair-wise bounding box collision of all the bounding boxes it owns with all other bounding boxes. That is, 8Bði,jÞ collide with all B(i,k) where B(i,j) and B(i,k) correspond to bounding boxes of two different grids. 4. If B(i,j) and B(i,k) collide, then compute the intersection Bði,jÞ \ Bði,kÞ and store it in the remote-connectivity list or local-connectivity list accordingly. The union of all collisions yields Rij .

2.2. Hole-cutting With the potential inter-grid communication maps Rij known, the next step is to remove those points that are outside the domain of interest. To do this, each grid Gi in the composite grid system is classified as either a background grid or a feature grid during the problem setup stage (see Figs. 1 and 3). Feature grids consist of solid boundaries, denoted by S, or are specially placed refinement grids. All other grids are called background grids. The cells of a background grid that are within, or intersect with, the solid boundaries of a feature grid are marked as inactive and must be removed from the calculation. This process is formally known as hole-cutting. A suitable definition of the cutting surface S is required to satisfy geometric queries. The present implementation employs a discrete, explicit, representation of S, that can be automatically extracted from the feature grids based on boundary condition information. Given the inter-grid connectivity Rij , hole points are identified. In addition, the mesh points surrounding the hole points that are within a user-supplied distance are identified to receive data from the background grid(s). Similarly, the surrounding fringe points of the feature grid are marked as receivers as well, thus forming a bidirectional exchange of the solution between the overlapping grids. To facilitate marking hole points and fringe points an integer-valued array, historically known as ‘IBLANK’ and given the symbol ib in the text to follow, is used at each point in the domain. Traditionally, the IBLANK array consists of three possible states, {0, 1, N}, whose values imply that the current point is inactive (blanked out), active, or receives its information from grid N, respectively. For example, a simple usage of ib might be qq ¼ minfjibj,1gRðqÞ qt for the system of equations q_ ¼ RðqÞ. When ib=1 the equations are updated as described by R(q), but for any other value ib r0, q_ ¼ 0 and another method is needed to update to the solution. In order to provide more fine-grained control, and to support periodic domains and partitioned grids, the traditional IBLANK definition has been extended as described in [31]. In this context, the main steps for the hole-cutting procedure are as follows. First, the values of ib are initialized to include periodic or interior boundaries and to identify the user-supplied fringe for interpolation (Fig. 3(a)). Those nodes in the fringe have an IBLANK value of ib =  N, where N is the grid number which will provide the required interpolation data. Second, the intersection region, Rij , of the feature grid and background grid is used to identify the portion of the

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

Background grid

4165

y-periodic boundary fringe imprint of Σ hole

y-periodic boundary

Fig. 3. Hole-cutting and donor–receiver search results using the 2-D rotor–stator 3-grid configuration. (a) The background grid after the hole-cutting and donor–receiver pair search operations are complete. (b) Sample spatial decomposition employed for the donor–receiver search on the rotor grid and the identified donor points.

background grid that is to be cut by the feature grid(s) to minimize the search space used during hole-cutting (Fig. 3(a)). Third, the cutting surface S on the feature grid(s) is determined whose interior(s) define a volume in which the background grid is to be turned inactive. Usually the interior volume corresponds to the interior of solid object, such as the rotor and stator blades of the SDT fan, which should be ignored for the current timestep. The surface S is then imprinted on the background grid, as shown in Fig. 3(a). Finally, a flood/fill operation [27] is used to mark all background nodes which are in the interior of S with an IBLANK value of ib= 0. These nodes will not be used by the flow solver in the current timestep. For periodic domains or grids that are partitioned to several processors the imprinted surface may not be watertight so that the flood/fill operation employed in the current implementation has been modified from Ref. [27] accordingly. 2.3. Donor–receiver pair search After all hole points are excluded, the next step is to identify the donor–receiver pairs within the associated overlapping region Rij . The process begins by extracting the candidate receiver points, R, in grid partition g(i,k) and the candidate donor cells, D, in grid partition g(j,m). Further, if the overlapping grids do not reside on the same process, the candidate receiver

4166

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

points are sent to the process that owns the donor grid. For every receiver point r 2 R the list of donors D is searched to identify the cell, d, that contains r. This is achieved using a Newton–Raphson iteration to solve for the local coordinates Ur ¼ ðx, Z, zÞT of r(x,y,z) with respect to the donor cell d as þ1 m 1 Um ¼ Um F ðUm r r ½ J  r Þ,

(1)

where m is the Newton–Raphson iteration counter and J is the jacobian matrix of the donor cell d [30]. A donor–receiver pair is then formed if 0e r Ur r1 þ e,8Ur 2 fx, Z, zg, where e is a user-supplied tolerance, typically 1  10  10. For all donor–receiver pairs formed, the ib value of the receiver point r is updated to indicate that a donor cell d has been found. Typically, this step entails communication if the receiver grid resides on a remote process. Fig. 3(b) shows the donating points on the rotor feature grid which will be used to pass data to the corresponding background grid receiving point in the fringe shown in Fig. 3(a). An example r–d pair is shown in the lower-right pane of Fig. 2 for a stator receiving point and a background grid donating cell.

2.4. Interpolation The last step in the grid assembly process is determining the interpolation weights for transferring the solution from the donor cell nodal data d to the corresponding receiver point r. (The flow solver used in this work is node-based and hence the focus of the paper is on node-based data instead.) Lagrangian interpolation is used for the solution transfer, hence, the data of the donor cell d is sufficient for forming the interpolation stencil. Given Ur , the interpolation weights are computed by Wi ¼ Ni ðUr Þ where Ni are the Lagrange shape functions for the cell d and i 2 ½1,k where k are the total number of nodes for the given donor cell d. Next, the value at the point is interpolated by multiplying the solution at the nodes of the donor cell, fi, as finterp ¼

k X

Wi fi :

(2)

i¼1

Finally, the solution associated with the receiver points r is updated with the interpolated solution from the corresponding donor grid. Further details on inter-grid data transfer and, specifically, what data are transferred is given in Section 3 below. The order of accuracy is controlled through choice of the interpolation shape functions. Other interpolation schemes which can be written in the form of Eq. (2) are possible.

2.5. Performance improvements for searching The approach outlined in Sections 2.1–2.3 is general and robust, but needs further refinement for dynamic problems where the grids are in relative motion, as is the case for the SDT fan of interest here. In effect Sections 2.1–2.3 outline a useful approach for static overset mesh problems suitable for parallel platforms. In a dynamic problem, where the grids are in relative motion, special care is taken to optimize all of the search queries since the steps in Sections 2.1–2.3 must be performed at least once per timestep. To this end accelerated searching techniques are employed which are based on spatial decomposition. The main idea is to minimize the search space of the algorithm by further subdividing the space of Rij such that for a given point query only a subset of the cells is visited. This is achieved using spatial data-structures such as an oct-tree [32,33], used in an earlier implementation [30], or a virtual grid. A virtual grid, shown in Fig. 4, is essentially a 2-D or 3-D hash table which divides the space of Rij into several buckets, S bi, viz, VðRij Þ ¼ ki ¼ 1 bi . It is virtual in the sense that a grid is not actually constructed. Instead, it is defined by (1) its bounds, Cmin = (xmin, ymin, zmin) and Cmax = (xmax, ymax, zmax), which are provided by Rij , and (2) the sampling step-size D ¼ ðdx , dy , dz Þ defined as a factor of the largest edge-length within region Rij .

Fig. 4. Example virtual grid. (a) Topology of the virtual grid in 2-D. (b) Virtual grid in 3-D enclosing of an object.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

The dimensions of VðRij Þ, denoted by ðDx , Dy , Dz Þ, can be computed as follows:       x x y y z z Dx ¼ max min , Dy ¼ max min , Dz ¼ max min :

dx

dy

dz

4167

(3)

Given this implicit definition and a query point q =(x,y,z) the corresponding bin bH is determined by first identifying the virtual grid coordinates, (i,j,k), given by i¼



xxmin

dx



,

j¼

  yymin ,

dy

k¼

  zzmin

dz

,

(4)

and the hashing function, Hði,j,kÞ, Hði,j,kÞ ¼ ðk1ÞDx Dy þ ðj1ÞDx þ i:

(5)

Time (s)

The cells of within region Rij are first mapped to the corresponding buckets using Eqs. (4) and (5) and each query q is mapped to the corresponding bucket bH in constant time. Then, within each bucket bH , the abbreviated list of potential

110 100 90 80 70 60 50 40 30 20 10 0

No cache

0.25% cache

0

1

2

3

4

0.15% cache

5 6 7 Iterations

8

9

10

11

Fig. 5. Sample results from the wing-store separation moving grid configuration. (a) Initial state of the wing-store configuration. (b) Wing-store after two iterations. (c) Graph of execution times per iteration with and without cache.

4168

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

cells is searched for the one containing q using the natural coordinates Fr described earlier. In doing so the spatial search uses only the subset of cells in the intersection region Rij . Once the subset of search cells bH is found they must still be examined. However, by observation it is easy to see that consecutive donor-cells are spatially close together. Thus, spatial caching is employed to pre-fetch the donor-cells for the upcoming queries and exploit locality to improve the performance. Once a donor-cell is found, its neighbors are cached into a search-pool of candidate donor cells. Consecutive queries then first look at the search-pool for a donor-cell and only resort to the nominal spatial search if a donor-cell is not found in the cache. The present cache implementation employs a least-recently used (LRU) donor-cell insertion/eviction policy implemented using a doubly linked-list. New cache entries are always placed in the beginning of the cache. Moreover, if there is a cache hit, it is moved to the beginning of the cache. In doing so, the most recently used cells and their neighbors are always in the beginning of the cache to exploit locality while the least-recently used cells are at the end. Using a doubly linked-list as the underlying data-structure enables the use of pointers to implement in constant time all cache insertion and eviction operations. However, searching the cache is currently an O(N) operation, where N is the number of cells in the cache, and care must be taken in choosing the cache size. Future work is focused in improving the cache-access-time by employing more frequent cache-purging as well as reducing the computational complexity to O(log(N)). The performance benefits from caching are illustrated using a three-grid system shown in Fig. 5, where the object undergoes vertical motion in the negative-z, downward, direction. Fig. 5(c) shows the execution times to find the donorcells of the corresponding points on the outer-boundary of the store grid for 10 iterations on a single processor. As it is illustrated, when no caching is employed the execution time is longer and increases for the first five iterations as the

4

y [in]

2 0 -2 -4 -5

0

5 x [in]

10

15

-5

0

5 x [in]

10

15

4

y [in]

2 0 -2 -4

2.5

y [in]

2 1.5 1 0.5 0 -1

0

1

2

3

4

x [in] Fig. 6. Magnitude of z-vorticity for the 2-D linear cascade model of the fan-rotor system, (a) shortly after start-up and (b) after the trailing edge of the rotor reaches the stator. Observe the lack of wake distortion (c) as it passes through the grid overlap regions.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4169

objects moves through an increasingly dense region of the grid prior to reaching a coarser region. In contrast, when caching is employed, the performance is drastically improved, even for the first iteration, and for consecutive iterations also. It is expected that logarithmic access to the cache, as well as a more sophisticated cache eviction policy, will yield additional improvements. 2.6. 2-D linear cascade fan-rotor system Returning to the motivating application, consider a 2-D linear cascade model corresponding to a 75 percent span cut of the 22-in SDT fan tested in the 90  150 acoustic wind tunnel facility at NASA Glenn. The system consists of three grids as was depicted in Fig. 1(b). The simulation proceeds with the rotor grid moving upwards along the y-direction. An approach condition of (7809 RPM 62 percent power) was used, yielding a rotor linear speed of local Mach number M= 0.56. To account for the y-periodic motion of the rotor through the implied infinite linear cascade array of the background mesh, a fourth grid was added to the system. The fourth grid is an exact copy of the rotor grid and moves at the same velocity. The main purpose of the fourth grid is to facilitate the periodic wrapping of the rotor grid moving passed the background grid. The grids were further partitioned and distributed to 64 processors for this computation. A sample flow field is visualized in Fig. 6 using the magnitude of the z-vorticity joz j ¼ jqv=qxqu=qyj to highlight the boundary layers and wakes stemming from the blades. Since the objective was on the overset mesh methodology, no attempt was made to quantitatively compare the simulation data with existing data. Note the grid-to-grid exchange from the boundary layers shedding from the rotor and stator grids. The rotor wake is approximately one chord away from the stator’s leading edge. As the rotor moves, the rotor wake is transferred to the background grid and the trailing vortices propagate and impact the stator. There is no distortion of the vorticity visible in the transfer, as shown in Fig. 6(c). The total time spent for grid-to-grid data transfer was approximately 0.060 s for a single iteration, with 35 percent of that time performing the hole-cutting and 57 percent of the time performing the donor–receiver search. For reference the wall time of the inviscid fluxes was 0.020 s while the viscous fluxes were 0.025 s. Thus using the moving grid capability cost roughly one extra evaluation of the fluxes, which is a manageable expense for all but the largest calculations possible on today’s resources. More importantly, the algorithms discussed above are not focused on turbomachinery problems, as are the chorochronic and some sliding mesh implementations, but are fully general to other moving grid problems. 3. Provably stable exchange of data between overset grids In the preceding section it was discussed how to efficiently transfer data between grids in an overset mesh framework, with particular attention paid to parallel efficiency. Now attention is turned to what data should be transferred. Previous work on multi-block solution transfer focused on ensuring conservation across the blocks by transforming the overlapped interface into a patched interface (see, for example, Wang [34]). Most Chimera-based codes, however, interpolate the conservative flow variables and ‘over-write’ the receiving points with the donated values (using the terminology of the previous section), much like an injection-based boundary condition. It is known, however, that injection boundary conditions can lead to instabilities, especially when using high-order non-dissipative numerical methods [35]. How to stably couple overset grids has received little attention, though Ref. [35] suggests that it can be done. In the traditional multi-block setting, where grids share faces but do not overlap, recent work has shown how to stably pass data [36]; when one of the grids is unstructured one may use the method of Ref. [37]. To this end we present one provably stable methodology for transferring data between overlapping grids using a class of existing finite difference operators. Proving stability for classes of operators is challenging and a few open issues remain, as we highlight below, but the effort is beneficial in that one is not restricted to implement a very specific scheme to yield the stated stability. Instead, a family of schemes is available. The essence of the method is given in this section for a simple hyperbolic equation; the more general problem of a hyperbolic system of equations is given in Section 3.2 with application to the linearized Euler equations. It is important to note that the methods to be presented do not use artificial dissipation, of any kind, for stabilization. Instead they utilize a provably stable method for passing data into and out of computational domains recently developed for hyperbolic and incompletely parabolic equations [38–40]. This approach is different than existing methods which rely on filtering [41] or on artificial dissipation combined with least squares interpolation [42] and ties more completely the numerical discretization of the spatial derivatives and of the interpolation. 3.1. Provably stable overset method for the advection equation To introduce the necessary ideas for proving an overset method is stable, we will first consider the simple, linear advection equation ut þ ux ¼ 0,

uðx,0Þ ¼ gðxÞ,

uða,tÞ ¼ f ðtÞ

(6)

on the one-dimensional domain O ¼ ða,bÞ. Eq. (6) has been nondimensionalized to have unit speed. The continuous problem is well-posed [43]. Before proceeding, a review of notation and necessary concepts is presented. Our interest is to find how to pass information from two grids which overlap in a region (such as in Fig. 7(a) and in the SDT example

4170

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

VL

x=0

hL

VR x=b

x=a x = –δ

hR

1.5 t=0

u (x, t) [-]

1

t=0.5

t=1

0 x [-]

0.5

0.5

0

−0.5 −1

−0.5

1

Fig. 7. One dimensional advection test problem showing (a) notation and setup and (b) solution at times t = 0, 0.5, 1. At t = 0.5 the solution is shown in the left (–&–) and right ð2J2Þ domains.

of Section 2) to yield a convergent, consistent, and stable numerical method. To do this we will consider how the numerical error, e, defined as the difference between the current numerical solution, v, and the exact solution u(x,t) evaluated at the grid points x ¼ fxi gN i ¼ 1, eðtÞ ¼ vðtÞuðx,tÞ, JeJ22

(7) JeJ2H

T

T

¼ e e we will be using matrix weighted norms ¼ e He, where behaves as t-1. In addition to the usual ‘2 norm H is a symmetric, positive definite matrix. For composite grid systems, it is beneficial to refer to approximate solutions and errors corresponding to the left grid and sets of values corresponding to the right grid (see Fig. 8(a)). Let vL and eL be the approximate solution and error on the left uniform grid with the vectors vR and eR defined similarly for the right, such that " # " # vL eL v¼ , e¼ : vR eR Our discussion proceeds by introducing the particular finite difference and boundary condition operators to be used before proving stability of the advection equation on two domains.

3.1.1. Summation-by-parts finite difference operators To discretize the derivative operators on the domain, h ¼ ðbaÞ=ðN1Þ, xi ¼ a þ ði1Þh,

i ¼ 1, . . . ,N,

(8)

we use a finite difference approximation which has a particular property, termed ‘summation-by-parts,’ or SBP, as follows. The pair of matrices {P,Q} is a pair of summation-by-parts (SBP) matrices of approximation order p if (a) P  1Qv is an order hp approximation to q=qx, (b) P is a symmetric positive definite matrix and (c) Q þ Q T ¼ diagð1,0,0, . . . ,0,1Þ. These conditions, as presented in Ref. [44], ensure that a discrete version of the ‘integration by parts’ property holds; namely, /P 1 Q x,ySP ¼ xN yN x1 y1 /x,P 1 Q ySP , where /x,ySP ¼ y Px is the usual P-weighted inner product in ‘2 . T

(9)

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4171

0.8 Energy ‘increase’ due to overlap double counting

2

||v||2 [-]

0.7

0.6

0.5

0.4 0

0.5

1

1.5

1

1.5

x [-] 1.4

1

2

||v||H [-]

1.2

0.8 0.6

0.4 0

0.5 x [-]

~ Fig. 8. Energy of solution, v ¼ ½vL vR T , in Fig. 7 in the (a) 2-norm JvJ2H and in the (b) energy stable H-norm JvJ2H~ . Dashed line corresponds to the energy at t =0 in the corresponding norm.

An example pair of second-order SBP matrices is 2 6 6 6 6 P ¼ h6 6 6 4

1 2

1 & 1 1 2

3

2

7 7 7 7 7, 7 7 5

6 1 62 6 6 Q ¼6 6 6 4

 12

3

1 2

0 &

1 2

&  12

& 0  12

7 7 7 7 7: 7 17 25 1 2

Higher-order SBP matrices may be found in Strand [44]. SBP operators based on the dispersion-relation-preserving approach of Tam and Webb [45] are also available [46]. There are many ways to construct pairs of SBP matrices {P,Q}. One can construct pairs for which P1 Q is up to a locally eighth-order approximation to qx [35,38,44,47]. If P is diagonal, then the operator P 1 Q is termed explicit, otherwise it is called implicit. One can show that assuming a certain structure for Q precludes P from being proportional to the identity [48]. We can extend SBP matrices to higher dimensions using the Kronecker product, denoted by . Consider a uniform discretization of a subset O  R2 , ðax ,bx Þ  ðay ,by Þ, with nx nodes in the x-direction and ny nodes in the y-direction. Let {Px, Qx} be a pair of nx  nx SBP matrices of approximation order p. Let {Py,Qy} be a pair of ny  ny SBP matrices of approximation order q. Then, define the following: H ¼ Px  Py , Gx ¼ Qx  Py , Gy ¼ Px  Qy :

4172

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

3.1.2. Simultaneous approximation term boundary condition The simultaneous approximation term (SAT) methodology is a process using SBP finite difference operators and enforcing boundary conditions such that the resulting method is energy stable [38]. The power of the SAT methodology is to reduce the problem of finding energy stable methods to enforcing that a small number of penalty parameters, say t, are set to the correct values. SAT boundary conditions have been shown to be effective in aeroacoustic problems on single domains by Bodony [17]. As an example, consider the advection problem (6) on the interval domain (a,b). Let {PQ} be a pair of SBP matrices. One possible SAT formulation of this problem using these matrices is, following Ref. [38], v_ ¼ P1 Q vtP1 sðv1 ðtÞf ðtÞÞ, vð0Þ ¼ gðxÞ, s ¼ ½1,0, . . . ,0,0T :

(10)

Observe that the boundary condition at the left point x =a is enforced weakly through the penalty term tP wðv1 f Þ. The method does not require v1, the value that approximates u(a,t), to be equal to f(t). Strictly enforcing boundary conditions can result in an unstable system [38]. Note also that the penalty term contains the penalty parameter t, a value that does not correspond to any quantity in the original PDE (6). The purpose of the penalty parameter is to enforce energy stability and convergence by being ‘large enough.’ For example, as shown in Ref. [38], when t Z 1=2 then method (10) is energy stable in the P-norm and converges to the true solution u at fixed time tf with order p, 1

Jeðtf ÞJP ¼ Oðhp Þ: Moreover, if t 4 12, then v also satisfies the stronger estimate rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi JvðtÞJP r Jvð0ÞJP þ

t

t2

max ðjf ðt 0 ÞjÞ: 2t1 t0 2½0,t

(11)

There are examples of SAT methods for several families of equations, including the compressible viscous Navier–Stokes equations [39,40,49,50]. For example, in [49], the authors construct an SAT method for the advection–diffusion equation with two penalty parameters. Energy stability requires that one penalty parameter be set to the diffusion coefficient, and the other to the diffusion coefficient multiplied by negative one. 3.1.3. Provably stable overset methodology for the advection equation We now return to stably coupling overset grids using the class of SBP non-dissipative finite difference operators. To simplify notation we adopt the following abbreviations: D ¼ P 1 Q þ tP1 ssT , c ¼ tP 1 s,

~ ¼ P 1 Q tP1 rrT , D

(12)

c~ ¼ tP1 r,

(13)

T

where r ¼ ½0, . . . ,0,1 . Using these, the ODE (10) is written v_ ¼ Dv þ f ðtÞc. Let {PL,QL} be a pair of NL  NL SBP matrices of approximation order p. Let tL be the left penalty parameter. Using PL, QL, and tL , let DL and cL be the derivative operator and penalty vector defined by (12) and (13) when P=PL, Q=QL and t ¼ tL . The right operators are defined similarly. Let I be an 1  NL matrix describing the interpolation from the left grid to the point x ¼ d on the right grid (as in Fig. 7(a)) such that IuðxL Þ ¼ uðd,tÞ þ einterp , where einterp is the interpolation error. The stable overset method is given by " # " # cL 0NL ,NR DL v_ ¼ Dv þf ðtÞd where D ¼ and d ¼ 0NR cR I DR

(14)

subject to the initial condition v(0) =g(x). This system will be Lyapunov stable if t1 Z 12 and t2 Z 12. To see this note that the system matrix D is block lower diagonal and hence its eigenvalue spectrum lies in the left-half-plane since the individual left and right domains are energy stable following Ref. [38]. Proving Lyapunov stability is not sufficient to show that the numerical method yields a convergent method because one must show how the eigenvalues of the system change as the grids are refined. It is shown in the Appendix that the method is also convergent as NL ,NR -1 and thus comprises a stable, consistent, and convergent numerical method. 3.1.4. Numerical example of provably stable overset methods for the advection equation Application of the above algorithm to the advection equation is shown in Fig. 7(b) where the initial condition uðx,0Þ ¼ expfðx þ 0:5Þ2 g convects to the right at unit speed. The solution at time t =0.5 shows that for a range of time the disturbance exists on two grids, leading to the increase in the 2-norm shown in Fig. 8(a), while the energy of the disturbance itself has not changed, as shown by the constant amplitude of the disturbance in Fig. 7(b) even after passing

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4173

~ is needed (discussed below) to remove the overlap-induced energy increase in a through the interface. A different norm H manner consistent with the underlying numerical discretization. 3.2. Provably stable overset methodology for systems of hyperbolic equations The methodology can be constructed for systems of equations in a manner similar to the advection equation, following a simple change of variables. To show this let a k  k, real symmetric matrix A be diagonalized by the orthogonal matrix T, A ¼ T LA T T

(15)

with LA ¼ diagðl1 , l2 , . . . , lk Þ where l1 Z l2 Z    Z lk are the necessarily real eigenvalues of A. Using A we can then define A þ ¼ TIðLA ÞLA T T and A ¼ AA þ and where I(A), the indicator of A, is a k  k matrix defined by ( 1 fAgi,j 4 0, fIðAÞgi,j ¼ 0 fAgi,j r 0: Using these definitions a general hyperbolic system for symmetric A is ut ¼ Aux , uðx,0Þ ¼ gðxÞ, A uða,tÞ ¼ fa ðtÞ, A þ uðb,tÞ ¼ fb ðtÞ:

(16)

The continuous system is well-posed [51]. The hyperbolic system (16) can be solved by decomposing it into characteristic variables and using the scalar method of (14). Let T be any orthogonal matrix that diagonalizes A, T T AT ¼ LA , and let wðx,tÞ ¼ T T uðx,tÞ be the system’s corresponding characteristic variables. Then it follows that w solves wt ¼ LA wx , wðx,0Þ ¼ T T gðxÞ, T L A wða,tÞ ¼ T fa ðtÞ,

LAþ wðb,tÞ ¼ T T fb ðtÞ, which is a system of decoupled advection equations. Applying the diagonal method of (14) to this new system and then transforming back gives the desired result that when tL 4 12 and tR 4 12, the following converges to the solution of (16) and is energy stable: ~  A þ Þ þ ðD  A ÞÞv þðd  fa ðtÞÞ þ ðd~  fb ðtÞÞ, v_ ¼ ððD vð0Þ ¼ gðxÞ:

(17)

To prove this claim we apply a similarity transformation to the system matrix to yield a decoupled set of advection equations. Because similarity transformations preserve eigenvalues, all eigenvalues of (17) also lie in the left-half-plane. 3.2.1. Numerical example of provably stable overset methods for systems of hyperbolic equations To demonstrate the above method (17) we use it to solve the nondimensional linearized Euler equations in two dimensions in the form of the First CAA Benchmark, Category 3 [3], qr q q ðMy r þ vÞ ¼ 0, þ ðMx r þuÞ þ qx qy qt

(18)

qu q q þ ðMx u þ pÞ þ ðMy uÞ ¼ 0, qt qx qy

(19)

qv q q þ ðMx vÞ þ ðMy v þ pÞ ¼ 0, qt qx qy

(20)

qp q q þ ðMx p þ uÞ þ ðMy p þ vÞ ¼ 0 qt qx qy

(21)

on a rectangular domain with initial conditions



rðx,y,0Þ ¼ pðx,y,0Þ ¼ exp ðln2Þ

x2 þ y2 9



,

(22)

4174

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

uðx,y,0Þ ¼ vðx,y,0Þ ¼ 0:

(23)

2 All quantities have been nondimensionalized by ambient quantities: pressure by r1 c1 , velocity by c1 , and a lengthscale ‘1 defined such that the wave travels one unit distance in unit time at the speed of sound. The exact solution is given in Ref. [3].

y

x

uL

uR

Left grid Right grid 1

y [-]

0.8

0.6

0.4

0.2

0 −1

−0.5

0 x [-]

0.5

1

0 −1

−0.5

0 x [-]

0.5

1

1

y [-]

0.8

0.6

0.4

0.2

Fig. 9. Overset mesh solution of the CAA Benchmark 1 Problem showing (a) grid arrangement, (b) initial condition (t = 0), and (c) pressure at time t= 0.5.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4175

||eL||HL,||eR||HR [-]

10-3

10-4 2 1 10-5 10-1.8

10-1.7 hx,L [-]

10-1.6

||eL||HL,||eR||HR [-]

10-4

10-5

10-6 4 1 10-7 10-1.8

10-1.7 hx,L [-]

10-1.6

Fig. 10. Numerical example of an acoustic pulse propagating in zero mean flow. (a) Overset decomposition with (b) convergence rate using second-order interpolation and (c) fourth-order interpolation. Legend: –&–, left grid; –n–, right grid.

Shown in Fig. 9(a) is the overset decomposition of the square domain into a ‘left grid’ and a ‘right grid.’ The numerical discretization uses the traditional fourth-order Runge–Kutta at constant timestep coupled to the fourth-order explicit SBP scheme [44]. The interpolation uses Lagrangian basis functions of either second or fourth order to transfer the characteristic variables across the overlap. An example of the solution with Mx = 0.5 and My =0 using second-order interpolation is shown in Fig. 9(b,c) where it is seen that the cylindrical pulse passes through the interface without distortion, even with a relatively low-order method. Some error is visible, however, near the pulse center in Fig. 9(c) which is controllable by increasing the scheme order and/or by increasing the mesh density. To perform the grid resolution study, the bounds of the individual domains were fixed as were the ratios of stepsizes to the x-direction spacing on the left grid, hx,L. As shown in Fig. 10(a) and (b) the order of accuracy in the entire domain, as measured by JeL JHL and JeR JHR , is dictated by the error of the interpolation scheme. Both the convergence rate and the absolute error at fixed grid resolution are strong functions of interpolation. When the location of the initial pulse is placed solely within the left grid we observe the convergence rates shown in Fig. 11(a) and (b). Notice that the left grid has a convergence rate close to the maximum achievable while the right grid convergence rate follows that of the interpolation error. The reason for this may be found by examining the problem in terms of the characteristic variables. For when the disturbance initiates in the left grid, as it propagates it passes into the right grid via the interpolation and so picks up the interpolation error. The right grid, on the other hand, is initially at uniform conditions and so the left-traveling characteristics have zero amplitude and are thus not affected by the interpolation. 3.3. Lyapunov versus ‘energy’ stability and viscosity Using these results one has a Lyapunov stable method wherein the stability is defined for the Cauchy problem as limt-1 uðtÞ ¼ 0, which is a weaker form of stability than ‘energy stable’ where dJwJ2H =dt r 0 8t where JwJ2H ¼ wT Hw for

4176

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

||eL||HL,||eR||HR [-]

10-1

10-2 2

1 3.5

1

10-3 10-1.8

10-1.7 hx,L [-]

10-1.6

10-1.7 hx,L [-]

10-1.6

||eL||HL,||eR||HR [-]

10-1

10-2 3.5

1

10-3 1

4.4

10-4 10-1.8

Fig. 11. Numerical example of an acoustic pulse propagating in zero mean flow, but with initial condition located solely within the left grid. (a) convergence rate using second-order interpolation and (b) fourth-order interpolation. Legend: –&–, left grid; –n–, right grid.

norm H. Moreover, Lyapunov stability does not guarantee convergence since the eigenvalues of the system may move towards the imaginary axis as the grid is refined, and some effort is required to show convergence independently of stability, such as given in the Appendix. Once the method is shown to be Lyapunov stable and convergent, however, the Lyapunov Theorem tells us that there exists at least one norm H for which dJwJ2H =dt r0 8t, and one can use the continuous Lyapunov equation to find the H. Proving that overset schemes are energy stable directly is challenging but avoids having to prove the method is Lyapunov stable and convergent since energy stability also implies convergence. Part of the difficulty stems from the overlap itself where the solution may exist on multiple grids simultaneously. An example is given in Fig. 8(b) where the two domain advection problem is solved with a Gaussian pulse initial condition which advects to the right at unit speed. As the pulse straddles the overlap, the energy in the norm H= diag(PL PR) increases for a brief period such that in this norm ~ in which the energy the scheme is not energy stable. By the Lyapunov theorem it is known there exists at least one norm H is strictly decreasing; computing one of these norms gives Fig. 8(c). The utility of the aforementioned theoretical results is muted somewhat by the observation that they depended on the system matrices between block triangular or, in other words, information passes from one domain to another along characteristic paths. When viscosity is added, even at the ‘simple’ level of the advection–diffusion equation ut þ ux ¼ euxx the proofs do not hold because although the individual domains are energy stable [38] the system is not necessarily so as the system matrix is block full. In on-going work we have developed an alternative approach which utilizes families of generalized SBP operators which modify the third SBP property in Section 3.1.1 to Q þQ T ¼ block_diagðBu ,0, . . . ,0,B‘ Þ where Bu and B‘ are negative definite and positive definite blocks, respectively, to overcome this difficulty.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4177

4. Conclusion In this article two aspects of numerical procedures were examined in the view of codes based on overset methods: (1) overset mesh management, and (2) stable interpolation. These items are unique to overset methods but their utility is based on the degree to which stable boundary conditions can be applied, such as discussed in Ref. [17]. For robust and accurate predictions of sound generating flows to be possible, methodologies for all three items must be present, and this article summarizes recent work towards that aim. In particular, details are given which demonstrate that application of modern computer science methodologies to overset grid management can open up new classes of complex geometry problems amenable to highly accurate structured-grid methodologies with moving components on massively parallel platforms. Moreover we have shown how to accurately and stably transfer information between overset grids in a manner consistent with a family of finite difference operators for problems involving strictly hyperbolic systems of equations. Further development of the two technologies presented in Sections 2 and 3 is ongoing and are focused on further performance improvement for the overset grid management and incorporating diffusion into the provably stable interpolation methodology.

Acknowledgments The first author (D.J.B.) is indebted to his students and collaborators at the University of Illinois for their contributions, especially Mr. M. Campbell (Computational Science and Engineering), Professor J. Freund (Aerospace Engineering and Mechanical Science and Engineering) and Professor M. Heath (Computer Science). In addition, D.J.B. acknowledges financial support from NASA Glenn (SBIR, Dr. E. Envia, monitor), NASA NRA (Supersonic fixed wing, Drs. J. Bridges and J. Debonis, monitors; Subsonic fixed wing, Dr. L. Hultgren, monitor), and from the Department of Energy. Appendix A. Proof of convergence of the overset method for the advection equation The proof assumes the following. First, let tL Z 12 and tR 4 12 as required for single domain energy stability [38]. Second, the pair of SBP matrices {PL, QL} approximate qx to order p, JuðxÞx PL1 QL uðxÞJPL rC1 hpL : and the pair of SBP matrices {PR, QR} approximate qx to order q, JuðxÞx PR1 QR uðxÞJPR r C2 hqR : The interpolation I has the property that (a) it is accurate to order r, juðd,tÞIuðxL Þj rC3 hrL and that (b) there is a constant C such that for vectors x 2 R

NL

(A.1)

the interpolation is bounded,

xT IT Ix o CxT x:

(A.2)

Finally, it is assumed that there is a constant m associated with the P-norm relating PL and RL ¼ hL INL ,

mJxJ2RL r JxJ2PL : With these assumptions the proof proceeds by first observing that, due to the lower block triangular structure of the system matrix, we can use the energy stability of the left domain to show that vL converges to u(xL). Then we note that the approximate solution on the right domain vR is the sum of three parts: (a) a component corresponding to exact boundary conditions and (b) two components corresponding to errors. Using bound (11), we show that (a) converges to u(xR) and (b) goes to zero. The rest of this section concerns the details of the proof and will use the following lemma: Lemma A.1. Let eL be the cumulative error associated with the left discrete domain. Then sffiffiffiffiffiffiffiffi C max Je J : max jIeL jr t2½0,tf  mhL t2½0,tf  L PL Proof. The proof follows from (A.2), qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi jIeL j ¼ eTL IT IeL r CeTL eL ¼ We now return to main proof.

sffiffiffiffiffi sffiffiffiffiffiffiffiffi C C Je J : JeL JRL r hL mhL L PL

&

4178

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

Proof. To show convergence, consider the left and right discrete domains separately. The left approximate solution vL is independent of the right approximate solution vR. Thus, because the single domain is energy stable it satisfies the bound JeL ðtf ÞJPL r tf C1 hpL : Now, decompose IvL into several pieces IvL ¼ IðuðxL Þ þ eL Þ ¼ uðaR ,tÞ þ einterp þIeL , where einterp is the interpolation error. Consider three subproblems ( _ ¼ DR w þ uðaR ,tÞcR , w R

wð0Þ ¼ g^ ð0Þ, (

(A.3)

:

y_ ¼ DR y þ ðeinterp ÞcR , yð0Þ ¼ 0NR ,

(

z_ ¼ DR z þ ðIeL ÞcR , vð0Þ ¼ 0NR :

By linearity, vR = w+ y+z. The error eR is equal to the error in w plus the solutions of the second and third ODEs, eR ¼ vR uðxR Þ ¼ w þy þ zuðxR Þ  ew þy þ z: Single domain energy stability [38] results thus bounds Jew JPR , JyJPR , and JzJPR : Jew ðtf ÞJPR rtf C2 hqR , sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jyðtf ÞJPR r

tf

pffiffiffiffi t2 max ðje jÞ ¼ S tf ðC3 hrL Þ, ð2t1Þ t2½0,tf  interp

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

sffiffiffiffiffiffiffiffi

Ct f t2 max ðjIeL jÞ rS Jzðtf ÞJPR r tf max Je J ð2t1Þ t2½0,tf  mhL t2½0,tf  L PL sffiffiffiffiffiffiffiffi sffiffiffiffi ! pffiffiffiffi Ct f C p1=2 p C th rS C t h ¼ S tf : mhL 1 f L m 1f L

Putting this together yields the desired result pffiffiffiffi JeR JPR ¼ Jew þ y þzJPR rJew JPR þJyJPR þ JzJPR rtf C3 hqR þ S tf C3 hrL þ

sffiffiffiffi C

! p1=2

C1 tf hL

m

,

which shows that the error on the right domain goes to zero as both the left and right domains are refined, and is thus convergent. & References [1] T. Colonius, S.K. Lele, P. Moin, Sound generation in a mixing layer, Journal of Fluid Mechanics 330 (1997) 375–409. [2] J.L. Estivalezes, L. Gamet, From jet flow computations to far-field noise prediction, FED-Vol. 238, 1996 (Fluids Engineering Division Conference, vol. 3, 1996). [3] J.C. Hardin, J.R. Ristorcelli, C.K.W. Tam (Eds.), ICASE/LaRC Workshop on Benchmark Problems in Computational Aeroacoustics (CAA), No. NASA CP-3300, 1995. [4] C.K.W. Tam, J.C. Hardin (Eds.), Second Computational Aeroacoustics (CAA) Workshop on Benchmark Problems, No. NASA CP-3352, 1997. [5] M. Dahl (Ed.), Third Computational Aeroacoustics (CAA) Workshop on Benchmark Problems, No. NASA CP-209790, 2000. [6] M. Dahl (Ed.), Fourth Computational Aeroacoustics (CAA) Workshop on Benchmark Problems, No. NASA CP-212954, 2004. [7] R.B. Langtry, P.R. Spalart, DES investigation of a baffle device for reducing landing-gear cavity noise, AIAA Paper 2008-0013, Presented at the 46th Aerospace Sciences Meeting and Exhibit, 2008. [8] A. Uzun, M.Y. Hussaini, High-fidelity numerical simulation of a chevron nozzle jet flow, AIAA Paper 2009-3194, Presented at the 15th AIAA/CEAS Aeroacoustics Conference and Exhibit, 2009. [9] K. Mahesh, G. Constantineescu, P. Moin, A numerical method for large-eddy simulation of compressible flows, Journal of Computational Physics 197 (1) (2004) 215–240. [10] J. Nordstrom, K. Forsberg, C. Adamsson, P. Eliasson, Finite volume methods, unstructured meshes and strict stability for hyperbolic problems, Applied Numerical Mathematics 45 (2003) 453–473. [11] Y. Hou, K. Mahesh, A robust, colocated, implicit algorithm for direct numerical simulation of compressible, turbulent flows, Journal of Computational Physics 205 (1) (2005) 205–221. [12] S. Mendez, M. Shoebyi, A. Sharma, F.E. Ham, S.K. Lele, P. Moin, Large-eddy simulations of perfectly-expanded supersonic jets: quality assessment and validation, AIAA Paper 2010-0271, Presented at the 48th Aerospace Sciences Meeting and Exhibit, 2010. [13] M. Wang, J.B. Freund, S.K. Lele, Computational prediction of flow-generated sound, Annual Review of Fluid Mechanics 38 (2006) 483–512. [14] E.A. Volkov, The method of composite meshes for finite and infinite regions with piecewise smooth boundary, Proceedings of the Steklov Institute of Mathematics 96 (1968) 145–185. [15] R.J. Magnus, H. Yoshihara, Inviscid transonic flow over airfoils, AIAA Paper 1970-0047, 1970.

D.J. Bodony et al. / Journal of Sound and Vibration 330 (2011) 4161–4179

4179

[16] J.A. Benek, J.L. Steger, F.C. Dougherty, A flexible grid embedding technique with applications to the euler equations, Proceedings of the 6th AIAA Computational Fluid Dynamics Conference, Danvers, MA, USA, 1983, pp. 373–382. [17] D.J. Bodony, Accuracy of the simultaneous-approximation-term boundary condition for time-dependent problems, Journal of Scientific Computing 43 (1) (2010) 118–133. [18] G.A. Gerolymos, G.J. Michon, J. Neubauer, Analysis and application of chorochronic periodicity in turbomachinery rotor/stator computations, Propulsion and Power 18 (6) (2002) 1139–1152. [19] R.L. Meakin, A.M. Wissink, Unsteady aerodynamic simulation of static and moving bodies using scalable computers, AIAA Paper 99-3302, Presented at the 14th Computational Fluid Dynamics Conference, Norfolk, VA, June 28–July 1 1999. [20] D. Belk, R. Maple, Automated assembly of structured grids for moving body problems, AIAA Paper 1995-1680, Presented at the 12th AIAA Computational Fluid Dynamics Conference, 1995. [21] N. Suhs, R. Tramel, PEGSUS 4.0 User’s Manual, Arnold Engineering Development Center Report AEDC-TR 91-8, Arnold Air Force Station, 1991. [22] N. Suhs, S. Rogers, W. Dietz, PEGASUS 5: an automated pre-processor for overset-grid CFD, American Institute of Aeronautics and Astronautics 41 (6) (2003) 1037–1045. [23] R.L. Meakin, N. Suhs, Unsteady aerodynamic simulaton of multiple bodies in relative motion, AIAA Paper 1989-1996, Proceedings of the 9th AIAA Computational Fluid Dynamics Conference, Buffalo, NY, 1989. [24] W.L. Sickles, A.G. Denny, R.H. Nichols, Time-accurate CFD predictions for the JDAM separation from an F-18C aircraft, AIAA Paper 2000-796, Proceedings of the 38th AIAA Aerospace Science and Exhibit, Reno, NV, 2000. [25] R.W. Noack, D.A. Boger, Improvements to suggar and dirtlib for overset store separation simulations, AIAA Paper 2009-340, Presented at the 47th AIAA Aerospace Science and Exhibit, Orlando, FL, 2009. [26] R.T. Biedron, J.L. Thomas, Recent enhancements to the fun3d flow solver for moving-mesh applications, AIAA Paper 2009-1360, Presented at the 47th AIAA Aerospace Sciences Meeting, Orlando, FL, 2009. [27] R.W. Noack, A direct cut approach for overset hole cutting, AIAA Paper 2007-3835, Presented at the 18th AIAA Computational Fluid Dynamics Conference, Miami, FL, 2007. [28] R.W. Noack, D.A. Boger, R.F. Kunz, P.M. Carrica, Suggar++: an improved general overset grid assembly capability, AIAA Paper 2009-3992, Presented at the 47th AIAA Aerospace Science and Exhibit, Orlando, FL, January 2009. [29] W.D. Henshaw, Overture: an object-oriented framework for overlapping grid applications, AIAA Paper 2002-3189, Presented at the 32nd AIAA Conference on Applied Aerodynamics, St Louis, MO, 2002. [30] G. Zagaris, D.J. Bodony, M. Brandyberry, M.T. Campbell, E.G. Shaffer, J.B. Freund, A collision detection approach to chimera grid assembly for high fidelity simulations of turbofan noise, AIAA Paper 2010-0836, Presented at the 2010 AIAA Aerospace Sciences Meeting, Orlando, FL, 2010. [31] G. Zagaris, M.T. Campbell, D.J. Bodony, E. Shaffer, M.D. Brandyberry, A toolkit for parallel overset grid assembly targeting large-scale moving body aerodynamic simulations, Presented at the 19th International Meshing Roundtable, October 2010. [32] H. Samet, The Design and Analysis of Spatial Data Structures, Addison-Wesley, 1990. [33] H. Samet, Applications of Spatial Data Structures: Computer Graphics, Image Processing, and GIS, Addison-Wesley, 1995. [34] Z.J. Wang, A fully conservative interface algorithm for overlapped grids, Journal of Computational Physics 122 (1995) 96–106. [35] B. Gustafsson, H. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, Wiley-Interscience, 1995. ¨ ¨ [36] J. Nordstrom, J. Gong, E.V. der Weide, M. Svard, A stable and conservative high order multi-block method for the compressible Navier–Stokes equations, Journal of Computational Physics 228 (2009) 9020–9035. ¨ ¨ [37] J. Nordstrom, F. Ham, M. Shoebyi, E.V. der Weide, M. Svard, K. Mattsson, G. Iaccarino, J. Gong, A hybrid method for unsteady inviscid fluid flow, Computers and Fluids 38 (2009) 875–882. [38] M.H. Carpenter, D. Gottlieb, S. Abarbenel, Time-stable boundary conditions for finite difference schemes involving hyperbolic systems: methodology and application for high-order compact schemes, Journal of Computational Physics 111 (1994) 220–236. ¨ ¨ [39] M. Svard, M.H. Carpenter, J. Nordstrom, A stable high-order finite difference scheme for the compressible Navier–Stokes equations, far-field boundary conditions, Journal of Computational Physics 225 (2007) 1020–1038. ¨ ¨ [40] M. Svard, J. Nordstrom, A stable high-order finite difference scheme for the compressible Navier–Stokes equations: no-slip wall boundary conditions, Journal of Computational Physics 227 (2008) 4805–4824. [41] S.E. Sherer, J.N. Scott, High-order compact finite-difference methods on general overset grids, Journal of Computational Physics 210 (2) (2005) 459–496. [42] C.K.W. Tam, F.Q. Hu, An optimized multi-dimensional interpolation scheme for computational aeroacoustics applications using overset grids, Presented at the 10th AIAA/CEAS Aeroacoustics Conference, AIAA Paper 2004-2812, 2004. [43] P. DuChateau, D. Zachman, Applied Partial Differential Equations, Dover Publications, Inc., New York, NY, USA, 2002. [44] B. Strand, Summation by parts for finite difference approximations for d/dx, Journal of Computational Physics 110 (1994) 47–67. [45] C.K.W. Tam, J.C. Webb, Dispersion-relation-preserving finite difference schemes for computational acoustics, Journal of Computational Physics 107 (1993) 262–281. [46] S. Johansson, High order finite difference operators with the summation by parts property based on DRP schemes, Technical Report 2004-050-NC, Uppsala, 2004. [47] B. Gustafsson, High Order Difference Methods for Time dependent PDE, Springer-Verlag, 2008. [48] H. Kreiss, G. Scherer, Finite element and finite difference methods for hyperbolic partial differential equations, in: C. De Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equations: Proceedings of a Symposium Conducted by the Mathematics Research Center, the University of Wisconsin-Madison, April 1–3, 1974. ¨ ¨ [49] K. Mattsson, M. Svard, J. Nordstrom, Stable and accurate artificial dissipation, Journal of Scientific Computing 21 (1) (2004) 57–79. ¨ [50] K. Mattsson, M. Svard, M. Shoebyi, Stable and accurate schemes for the compressible Navier–Stokes equations, Journal of Computational Physics 227 (2008) 2293–2316. ¨ ¨ [51] J. Nordstom, M. Svard, Well-posed boundary conditions for the Navier–Stokes equations, SIAM Journal of Numerical Analysis 43 (3) (2005) 1231–1255.