UFDC Home  myUFDC Home  Help 



Full Text  
A FINITE VOLUME, CARTESIAN GRID METHOD FOR COMPUTATIONAL AEROACOUSTICS By MIHAELA POPESCU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Mihaela Popescu To my daughter Ileana Klein, my parents, my brother and my sister, for their love and support ACKNOWLEDGMENTS I would like to express my deepest gratitude to my advisor Dr. Wei Shyy for his subtle yet effective methods of encouraging all of his students. I will always value my good fortune to have been one of them. Seldom does one encounter individuals with intellectual caliber, scientific temperament and a spirit of humanity that he embodies. I am equally grateful for his enduring enthusiasm and boundless patience during the process of preparing me to be a researcher and contributor to the global scientific community. Without his enthusiasm and commitment to excellence, my research in this field/area could not have been accomplished. I would like also like to express my deepest appreciation for Dr. Mark Sheplak who assessed my strength and successfully directed my research toward acoustics. I have benefited substantially from his guidance, experience, knowledge, and philosophy both professionally and personally. Similarly, I thank sincerely my committee members Dr. Lou Cattafesta, Dr. Nam Ho Kim and Dr. Jacob NanChu Chung for their support, encouragement and sharing freely of their expertise whenever it was needed. Being part of the Computational Fluid Dynamics (CFD) group was a source of pride and honor. I have always felt fortune to be part of this group and am grateful especially for their genuine friendship and support. I benefited from both collaboration and the pleasant work environment offered by members and visitors of CFD group. They were a real team as well as family for me. Our advisor and friend, Dr. Wei Shyy, fostered close and sincere relationship between members and visitors which enhanced the experience for all of us. I would like to thank my cousin Catalina and Laly Chirita, my friend Rada Munteanu, my friends from home Florina Carabet, Cosmin Carabet and Lascu Luminita for finding the most extraordinary and creative ways of being there for me whenever I needed them. I would like to thank to my family members: my parents, Aurica and Conatantin Popescu; my brothers and sisters, Pavel, Teofil, Marinela, Daniel Popescu, and Tabita Chirita and their family, for their sustaining support and love throughout years. They were the undisclosed source of energy and strength that I relied on during challenging times Finally, I would like to express my deepest love and appreciation to my daughter Ileana Klein, who demonstrated perseverance and resilience in face of unfavorable prognosis and the challenging circumstance of her birth. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST O F TA B LE S ........ .... ... .......... ...... ........................ ........ ...... ....... ix LIST OF FIGURES ............................... ... ...... ... ................. .x A B S T R A C T .............................................. ..........................................x iii CHAPTER 1 COMPUTATIONAL ASPECTS IN AEROACOUSTICS .........................................1 In tro d u ctio n ...................................... .................................. ................. . Definition of Sound ...................... ..................................3 The C characteristics of Sound........................................................ ............... 3 V iscous Effect in Sound W ave.................................... ........................... ......... 5 Classification of Aeroacoustic Problems .......................................... ...............6 Linear Problems in Aeroacoustics..................... .............................. 6 N onlinear Problem s in Aeroacoustics ....................................... ............... 9 L ightill's acoustic analogy ........................................ ........ ............... 10 Ffowcs Williams Hawkings equation ......................................14 The merits of Lighthill and Ffowcs Williams Hawkings analogies.......... 16 Nonlinear problem: Shock wave formation ...........................................17 Computational Techniques for Aeroacoustics............... ..........................................19 Direct Numerical Simulation (DNS) ............ ............................................. 19 Perturbation Technique ............................................... ............................ 20 Linearized Euler Equation.................................................... 21 Com putational Issues ........ .................... ............................. ............... 21 The Numerical Approach to Reduce Dissipation and Dispersion.....................22 C om plex G eom etry ....................... .. ........................ .. ...... ........... 24 S c o p e ........................................................................... 2 7 2 ASSESSMENT OF DISPERSIONRELATIONPRESERVING AND SPACE TIME CE/SE SCHEMES FOR WAVE EQUATIONS..............................30 Introduction ....................... ......................................... ................... 30 The DispersionRelation Preservation (DRP) Scheme................... ....... .........32 D iscretization in Sp ace ......................... .. ............................................ .. .. ..32 Tim e D iscretization ....................................... ...... ............ .. ............... 37 The SpaceTime Conservation Element and Solution Element Method....................41 a S ch em e ....................................................... ...... 4 1 a E S ch em e .............................. ..................................................... ............ 4 6 Numerical Assessment of the DRP and SpaceTime Schemes ..............................49 Short W ave: b/ x = 3 ............................. ..... ........................................ 57 Interm ediate W ave: b/Ax = 6................................ ................................. 58 L ong W ave: b/ x = 20.............................................................. ............... 58 Sum m ary and C onclu sions .............................................................. .....................62 3 FINITE VOLUME TREATMENT OF DISPERSIONRELATION PRESERVING AND OPTIMIZED PREFACTORED COMPACT SCHEMES FO R W A V E PR O PA G A IO N ......................................................... .....................63 N um erical Schem es .......................... .............. ................. ..... ...... 65 D R P S ch em e ................................ ................................................ 6 6 Finite volumebased DRP scheme (DRPfv) ............................................66 Boundary treatment of the DRP scheme............. ..... .................69 O P C S ch em e ..................... .. .. ...... ...... ............................................ 7 0 Finitedifferencebased optimized prefactored compact (OPCfd) scheme.70 Finite volumebased OPC scheme (OPCfv) ............................................71 The boundary treatment of the OPC scheme .............................................71 Time Discretization The Low Dispersion and Dissipation RungeKutta (LDDRK) Method................................. ............... 72 Analytical Assessment of DRP and OPC Schemes..............................................77 O p eratio n C o u n t ........................................ .............................. .................... 7 7 Dispersion Characteristics ................... ......... ..... ............78 S ta b ility ................. ......... ......... .. ............................................ 7 9 Computational Assessment of the DRP and OPC Schemes ........... ...............82 Test problem 1: OneDimensional Linear Wave Equation .................................82 Test problem 2: OneDimensional Nonlinear Wave Equation............................85 Test problem 3: OneDimensional Nonlinear Burgers Equation ......................89 Test problem 4: TwoDimensional Acoustic Scattering Problem.....................91 Sum m ary and C onclu sions .............................................................. .....................94 4 A FINITE VOLUMEBASED HIGH ORDER CARTESIAN CUTCELL METHOD FOR COMPUTATIONAL AEROACOUSTICS.........................98 Intro du action ...................................... ................................................ 9 8 CutCell Procedure ............................... ... .. ....... .... ............... 99 T e st C a se s ................................................................... ...................................1 0 4 Radiation from a Baffled Piston............................................... .................. 104 G general description ............................................................................. 104 D irective factor D ......................................... .. ...... ... .. ........ .... 107 Pressure on the face of the piston............... ............................................. 109 Low frequency (ka = 2) ............... ....... ............ ........... .......... 110 High frequency (ka = 7.5) ............................................................. 111 Reflection of a Pulse on an Oblique W all ...................................................... 112 Wave Generated by a Baffled Piston and Reflected on an Oblique Wall .........115 C conclusion ........... ................................................... .................. 117 5 SUMM ARY AND FUTURE W ORK ................................................. ............... 119 Assessment of DRP and SpaceTime CE/SE Scheme...........................................121 FiniteVolume Treatment of DispersionRelationPreserving and Optimized Prefactored Com pact Schem es ....................... ... ........ ...... ............... .... 122 Cartesian Grid, CutCell Approach for Complex Boundary Treatment............... 123 Future W ork ............ .. ....... ................................... ......................... 124 LIST O F R EFER EN CE S ......... .................................................. ...... ............... 125 BIOGRAPHICAL SKETCH ...... ........ .. ................. .... ....................... 134 LIST OF TABLES Table p 21 The stencil coefficient for N = 3 ..................................................... ............. 36 31 A summary of proposed CAA algorithms............................................ ..........96 32 The computational cost for DRP and OPC schemes............................ ............97 41 Published cutcell approach for different problems .............. ...... ....................118 LIST OF FIGURES Figure p 11 Sound propagation aw ay from a source ........................................ .....................3 12 Moving surface Ffowcs Williams and Hawkings equation ...............................15 13 Schematic diagram showing the (2n + 1) point stencil on a nonuniform g rid .............................................................................. 2 4 14 Cartesian boundary treatment of curved wall surface; b) detail around boundary ..................................... .................. ............... ........... 27 21 a~x versus ahx for the optimized DRP 4th order scheme, 7 point stencil, standard 6th order central scheme and 4th order central scheme ..............35 22 Dispersive characteristics of DRP scheme........................ ..............36 23 Scheme of the solution elements (SEs) and conservation elements (CEs) ..............43 24 Comparison between analytical and numerical solutions Effect of E on the accuracy of spacetim e a scheme ........................................ ............... 52 25 The dependence of the error on E for the spacetime ac scheme at t = 200: a) bAx = 3; b) b/Ax = 6; c) b/x = 20............................. .. ............. 53 26 The dependence of the error as function of v for short (b/Ax = 3), intermediate (b/Ax = 6) and long (b/Ax = 20) waves.................................... 55 27 Effect of von the accuracy of space time ac scheme: b/Ax = 3, t = 200 ..............56 28 The behavior of the error in function of the wavelength: comparison between DRP and spacetim e aE schemes ....................................................... 57 29 Accumulation of the error in time for short wave b/Ax = 3 ..............................59 210 Accumulation of the error in time for intermediate wave b/Ax = 6.......................60 211 Accumulation of the error in time for long wave b/Ax = 20..............................61 31 Grid points cluster for onedimensional problem ......................................... 67 32 Grid notation for twodimensional problem ................................. ............... 69 33 Foursixstage optimized RungeKutta of order four scheme: a) dissipation error; b) phase error ........................................ ........ ............... 74 34 Dispersive characteristics of the schemes ..... ......... ....................................... 80 35 Phase speed error on a logarithmic scale ...................................... ............... 80 36 Errors with respect to the time step size under a fixed space Ax, at t = 50 linear w ave equation ...................................... ............... .... ....... 84 37 Errors under a fixed CFL = 0.5, at t = 50 linear wave equation: .......................85 38 Errors with respect to the space step size under a fixed CFL = 0.5, at t = 5; nonlinear wave equation .......................... ................... ............ ... 87 39 DRPfd solution nonlinear wave equation; t = 5; CFL = 0.5 ................................87 310 DRPfv solution nonlinear wave equation; t = 3; CFL = 0.5 ................................88 311 OPCfd solution nonlinear wave equation; t= 5; CFL = 0.5................................88 312 OPCfv solution nonlinear wave equation; t= 5; CFL = 0.5................................88 313 Error as a function of Pe nonlinear Burgers equation ................ ...................90 314 Numerical solution obtained by DRP schemes nonlinear Burgers q u a tio n ................................................... ..................... ................ 9 1 315 Numerical solution obtained by OPC schemes nonlinear Burgers equ atio n ............................................................................ 9 1 316 Instantaneous pressure contours at time t = 7; twodimensional acoustic scattering problem .................. ........................................... ...... 93 317 The pressure history at point A, B and C ...................................... ............... 93 41 Illustration of the interfacial cells and cutandabsorption procedures.................. 100 42 M modified cut cell approach for CAA..................................... ............... 101 43 Radiation from a baffled piston (test problem (ii)) ............................................ 105 44 Radiation from a baffled piston: Beam patterns D(O); ka = 2 ...........................108 45 Radiation from a baffled piston: Beam patterns D(O); ka = 7.5 ........................108 46 Radiation from a baffled piston: Beam patterns D(O); ka = 12.5 ......................108 47 Radiation from a baffled piston: Real of piston radiation impedance....................109 48 Radiation from a baffled piston: Numerical solution on axis (ka = 2) ................10 49 Radiation from a baffled piston: Comparison between analytical and com puted solutions on axis (ka = 2) ............... ............................... .............. 111 410 Radiation from a baffled piston: Contour plot of pressure (ka = 2)....................111 411 Radiation from a baffled piston: Comparison between analytical and computed solutions on axis (ka = 7.5) ............. ..... ...................... .......... 112 412 Radiation from a baffled piston: Contour plot of pressure (ka = 7.5) ...................112 413 Reflection of the pulse on an oblique wall (test problem (ii)): general description ........... .................................................. ............. ...... ... 113 414 Reflection of the pulse on an oblique wall (test problem (ii)): cell arou n d b ou n d ary ..................... .. ........................................................ .... 1 14 415 Reflection of the pulse on an oblique wall: history of pressure for different angle of w all ................................................. ... ...... .. ........ .... 115 416 Reflection of the pulse on an oblique wall: a = 630 ............. ......... ... ............115 417 Wave generated by a baffled piston and reflects on an oblique wall: General description of the domain ................... ... ................. .............. 116 418 Wave generated by a baffled piston and reflects on an oblique wall: a = 630; a) t= 9; b) t 14; ........... ................ ............... ......................... 116 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A FINITE VOLUME, CARTESIAN GRID METHOD FOR COMPUTATIONAL AEROACOUSTICS By Mihaela Popescu August 2005 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering Computational Aeroacoustics (CAA) combines the disciplines from both aeroacoustics and computational fluid dynamics and deals with the sound generation and propagation in association with the dynamics of the fluid flow, and its interaction with the geometry of the surrounding structures. To conduct such computations, it is essential that the numerical techniques for acoustic problems contain low dissipation and dispersion error for a wide range of length and time scales, can satisfy the nonlinear conservation laws, and are capable of dealing with geometric variations. In this dissertation, we first investigate two promising numerical methods for treating convective transport: the dispersionrelationpreservation (DRP) scheme, proposed by Tam and Webb, and the spacetime ac method, developed by Chang. Between them, it seems that for long waves, errors grow slower with the spacetime ac scheme, while for short waves, often critical for acoustics computations, errors accumulate slower with the DRP scheme. Based on these findings, two optimized numerical schemes, the dispersionrelationpreserving (DRP) scheme and the optimized prefactored compact (OPC) scheme, originally developed using the finite difference approach, are recast into the finite volume form so that nonlinear physics can be better handled. Finally, the Cartesian grid, cutcell method is combined with the highorder finitevolume schemes to offer additional capabilities of handling complex geometry. The resulting approach is assessed against several well identified test problems, demonstrating that it can offer accurate and effective treatment to some important and challenging aspects of acoustic problems CHAPTER 1 COMPUTATIONAL ASPECTS IN AEROACOUSTICS Introduction Aeroacoustics deals with the sound generated by aerodynamics forces or motions originating in a flow rather than by externally applied forces or motions of classical acoustics. Thus, the sounds generated by vibrating violin strings and loudspeakers fall into the category of classical acoustics, whereas the sound generated by turbulent flow or the unsteady aerodynamic forces on propellers falls into the domain of aeroacoustics. A main feature of aeroacoustics is the large difference between the energy levels of the unsteady flow fluctuation and the sound. This is true even for a very loud noise. For example, in the nearacoustic field of a supersonic jet (at about 10 jetdiameters away) the acoustic disturbance amplitudes are about three orders of magnitude smaller than flow disturbance (Lele, 1997). For supersonic jets one percent of the noise is from the mechanical power of the jet. In many other cases, the efficiency can be much smaller; thus the amplitude of disturbance is much less. Another issue in acoustics is the difference between the scales of the unsteady flow and sound. This phenomenon is evident in situations in which flow speed is much less than the sound speed that characterizes the medium. This is because the time scale of the flow and the sound must match. In low Mach number flow (M o 1) this will give an acoustic wavelength that is M1 times longer than the flow length scale. In this case, a direct numerical simulation of the dynamic flow and generation of sound field will not be possible. An important challenge in computational aeroacoustics (CAA) lies in the coupling of time and space that appears in the acoustic wave. An acoustic wave has a wavelength X in space, as well as a frequency o) = 27f in time. These are coupled by the relation Xf = c, where c is the speed of sound in the medium (assumed quiescent). Thus, to determine an acoustic wave, one must resolve both its wavelength and its frequency. Numerical spatial and temporal approximation meets the same problem: the dissipation must be reduced in space, like that in time. However, the CFD community has done much work to overcome these difficulties, and high accuracy numerical schemes have been developed. CAA schemes try to minimize dissipation and dispersion error in both space and time and still maintain a certain order of accuracy. A further problem lies in the effect of the grid on the solution. In CFD, grids are often stretched to provide high resolution in regions of high gradients with lesser resolution where gradients are less steep. However, it is very difficult to propagate an accurate wave in this kind of grid. The dissipation and dispersion characteristics depend upon the CourantFriederichLewy (CFL) number (= cAt/Ax in one dimension), which can be interpreted as the distance an acoustic wave travels in one time step. If CFL is changing because the grid has been stretched, strange phenomena can occur. Vichnevetsky (1987) showed that if a wave is propagated over an expanding grid, the wave could actually appear to change frequency and be reflected such that it starts propagating back in the other direction. Anisotropy can appear because of the different sized grid spacing along different directions. An approach to solving the problem was suggested by Goodrich (1995), who recommended approximating the solution of the governing partial differential equations rather than approximating the derivatives and solving the resulting numerical scheme. Boundary conditions for CAA are also problematic. Often one is interested in solving the problem in an unbounded space. But for such computation, boundary conditions have to be imposed; thus the computational domain is finite. The first problem is the selection of the boundary condition to be imposed. Clearly, the proper boundary condition depends upon what is outside the computational domain: e.g., source, reflection boundary and free space. The second problem is in the implementation of these conditions. Definition of Sound The Characteristics of Sound A wave is the movement of a disturbance or piece of information from one point to another in a medium (excluding electromagnetic waves, as they do not require a medium). Sound is a wave that moves parallel to the direction of the propagation source Figure 11. Sound propagation away from a source As a sound wave propagates, it disturbs the fluid from its mean state. These disturbances are almost always small. We will consider departures from a state in which the fluid is at rest with a uniform pressure po and density po. When this is perturbed by a sound wave the pressure at position x and time t changes to po + p' (2, t), the density changes to po + p'( t), and the fluid particle moves with a velocity i (x, t). The ratios p'/ p, and p'/ p, are usually much less than unity. To illustrate the magnitude of the disturbance in an acoustic wave we will consider pressure waves in air which can be detected by the human ear. The pressure waves are referred to as sound. For harmonic pressure fluctuation (a pure tone), the typical range of frequency that our ear can detected sound is 20 Hz The range described by (1.1) is called the audible range. The maximum sensitivity of the ear is within the frequency range of 2 kHz to 3 kHz (policeman's whistle tone) The acoustic pressure of intense rocket engine noise can be as much as 9 orders of magnitude greater than the pressure of the weakest sound detectable by human ear. A logarithm scale was necessary to be able to comprehend this large range: 9 decade in amplitude (18 decades in intensity) (Blackstock, 2000). The logarithm devised is called a level, e.g., sound pressure level and intensity level. Although the levels are unitless, they are expressed in decibels (dB). SPL = 20 loglo(p 'pf) (1.2) On this scale a fluctuation of one atmosphere in pressure corresponds to 194 dB. The threshold of pain is between 130 and 140 corresponding to a pressure variation of only one thousandth of an atmosphere, i.e. p'/po = 103 At the threshold of pain, fluid particles in a 1 kHz wave vibrate with a velocity of about 0.1 m/s, which is only about 1/104 of the speed at which the sound wave is traveling. The displacement amplitude of the particles is thus between 104 and 105 meters and the wavelength is about onethird of a meter, so the displacement amplitude is only about 104 of the wavelength sound (Dowling and Ffows Williams, 1983). In the case where sounds can be approximated by small perturbations, several effects are noteworthy. First, there is no interaction between different acoustic waves; the sound fields add linearly. Hence, for the study of sound that comes from more than one source, we can study the sound from each source separately, and then add the final solutions. Second, the flow variables satisfy the linearized equation of the fluid motion, so each flow variable is linearly related to every other variable. This leads to a great simplification in the mathematics. Knowledge of one variable, e.g., pressure, provides a basis for a simple evaluation of all the other variables such as the density or particle velocity. Third, if we wish to solve the behavior of the sound numerically, we have to use a scheme that adds a low level of numerical dissipation. Viscous Effect in Sound Wave The effect of Reynolds number on CAA is ambiguous. Although all sound is ultimately dissipated into heat by viscosity, acoustics is generally considered to be inviscid fluid phenomena. If the viscous term is considered in the standard linear analysis, one finds its value to be Re = 27tcAp 0/ (1.3) where c is the speed of sound, A is the wavelength, po is the density, and 7 is the coefficient of viscosity. The values in air of these parameters for most practical interests are (Blackstock, 2000) c= 331.6 A 0.33 p0 1.293 (1.4) r7 1.7. 105 Based on these values, Re is approximately 4.8 x107 at the most audible frequency. This high value of Reynolds number determines that viscous effects are negligible in a sound field, because the pressure represents a far greater stress field than that induced by viscosity. Hence, we will regard sound as being essentially a weak motion of an inviscid fluid. The effect of viscosity can be taken into consideration after the sound travels about 27iAlpo/r wavelengths (i.e., after 1.5x108 m). Thus, dissipative loss becomes important for high frequency sound propagation over long distances. On the other hand, if one considers the generation of sound by flows rather than propagation of sound through flows, then a significant Reynolds number effect can be observed. For example, the most common source of sound in flows is the acceleration of the vorticity (Powell, 1964), which is only present in viscous flow. Even there, a curious independence is observed that is apparently due to the fact that the largescale; hence efficient sound generation structures in the flow change little with Reynolds. Classification of Aeroacoustic Problems Diverse problems of aeroacoustics can be classified on the basis of the physical phenomena that are expected to play a central role. The main classification is between the linear and nonlinear problem. Linear Problems in Aeroacoustics A major strength offered by computational approaches is the generality to deal with linear interaction problems of aeroacoustics (where one type of physical disturbance is scattered into another type of disturbance). A linearized solution of the unsteady Euler equations subject to appropriate boundary condition can provide useful prediction of the noise generation. The problems of linear interaction are not limited to sound generation but also contain the "reciprocal" problem of the generation (or receptivity) of vertical disturbance due to the incident sound (Colonius and Lele, 2004). This category includes the classical boundary value problems of linear acoustics: sound propagation in a uniform media in presence of reflecting surfaces, barriers, absorbing walls, and ducts, and also propagation and scattering of sound in a prescribed nonuniform medium. Specific examples include as: radiation from a duct opening or engineinlet due to a specified source or a specified ductmode disturbance, radiation from a specified ductmode disturbance, and radiation from a specified source across a finite barrier/sound wall with an absorbing treatment. Sound propagation in a specified nonuniform time dependent medium including refraction/scattering in steady and unsteady vortical flows, sound propagation in non uniform ducts including the interaction with geometrical changes, linear impedance and meanflow variations are considered as linear problems of scattering. The equation of a linear acoustics wave is deduced for a homogeneous fluid characterized by p = po, P = po, and i = fi0 + u'. The sound waves minutely disturb the status of the quiescent fluid P = Po + S, Sp << po P P= o << Pc2 (1.5) S= 0+ u', i' < To obtain the wave equation, we will start with the Euler equation and the equation of state for an isentropic process p, + u Vp +pV u = 0 0 \(1.6) p [u,+ (u)u +VP=0 and p = Cp (1.7) Substituting Eq.(1.5) into (1.6), the following is obtained 3p, + u VSp + p0V u + 3pV u = 0 u + (U V u+,pu +p (U +V 0 (1.8) The first order terms are small because the perturbation is very small; hence, second order terms or higher can be neglected. The underlined terms are second order or higher. The result is 8p, + poV u = 0 Pout +Vp = 0 (1.9) p = c Sp To reduce the set to a single equation in ii, the equation of state Eq.(1.7) first is used to eliminate Vp from the second expression in Eq. (1.9) PoUt + cV V(p)=0 (1.10) If the time derivative of this equation is subtracted from the divergence of the first Eq. (1.9), the result is the classical wave equation u cV2U 0 (1.11) Nonlinear Problems in Aeroacoustics Nonlinear aeroacoustics is the field of acoustics devoted to very intense sound, specifically to waves with amplitude high enough that the smallsignal assumption is violated. Retaining the nonlinear term makes the conservation equation much harder to solve. This includes problems of nonlinear propagation, such as nonlinear steeping and decay, focusing, viscous effects in intense sound beam (acoustics near field of highspeed jets), and sonic boom propagation through atmospheric turbulence, sound propagation in complex fluids and multiphase systems (such as in lithotripsy), and internal flows of thermoacoustic cooling system. Also included in this category are problems of scattering of nonlinear disturbances into sound, such as airframe noise. The basic equations that describe the nonlinear behavior are the same that describe the general motion of a viscous, heat conducting fluid: mass conservation, momentum conservation, entropy balance, and thermodynamic state. The mass conservation or continuity equation is DP V i = 0 (1.12) Dt where p is the mass density, ii is the fluid velocity vector, and D/Dt = a8/t + i V is the total, or material, time derivative. The momentum equation may be written as +VP=/V2i++ /B +1 V(V ii) (1.13) Dt 3) where P is the pressure, /u is the shear viscosity, and / is the bulk viscosity. Shear viscosity accounts for diffusion of momentum between adjacent fluid elements having different velocity. Bulk viscosity provides an approximation description, valid at low frequency, of nonechilibmm deviation between actual local pressure and thermodynamic pressure. Equation of state is given by P = P(p,T) or (1.14) P =P(p,s) where Tis the absolute temperature, and s is the specific entropy. A commonly used equation of state of a perfect gas is P = P exp[(s s)/c ] (1.15) where Po, po, and so are reference values of the pressure, density, and specific entropy. y= cp/c, is the ratio of the specific heats at heat pressure (cp) and constant volume (c,). In case of arbitrary fluid the equation of state is obtained by expanding of Eq.(1.14) in a Taylor series about (po, so). A more general description take implies relaxation, like vibrational relation of diatomic molecules (as in air) and chemical relation in seawater (Hamilton and Blackstock, 1997). The former occurs when the energy associated with molecular vibration fails to keep in step with molecular translation energy associated with the fluctuating temperature in gas. Lightill's acoustic analogy The study of flow that generates acoustic waves began with Gutin's theory (1948) of propeller noise, which was developed in 1937. However, the theory could not be considered until 1952, when Lighthill (1952, 1954) introduced his acoustic analogy to deal with the problem of calculating acoustic radiation from a relatively small region of turbulent flow embedded in an infinite homogeneous fluid. It is well known that a gas (not monoatomic) allows three distinct fundamental modes of oscillation: the vortical and entropy modes, both of which are convective, and the acoustic mode, which is the solution of the wave equation (Chu and Kovasznay, 1958). The three modes may exist independently only if i) the oscillations are small, and ii) the base flow is the uniform medium at rest or is in uniform motion. If the flow is inhomogeneous, then they are coupled, and it is not easy to separate them. To overcome these limitations, an approximation is necessary. In this section, the acoustic approach introduced by Lighthill (1952, 1954) is presented. This approach is used to calculate acoustic radiation from relatively small regions of turbulent flow embedded in an infinite homogeneous fluid where the speed of the sound co and density po is constant. Density and pressure fluctuations are defined by P = P P (1.16) P'= PPo where the subscript '0' is used to denote constant reference values that are usually taken to be corresponding properties at large distance from the flow. Lighthill's basic idea was to reformulate the general equations of gas dynamics in order to derive a wave equation suitable to describe sound propagation. The continuity and momentum equation are ap a + pu = 0 at a8x Sp (1.17) au, + x 8a 8+ at 8 ax, ax where o, is the component of the viscous stress tensor. For a Stokestian gas it can be expressed in terms of the velocity gradient Uau+ Oiu 2 Iu o u = p 2 + 3 (1.18) [ ax x 3 ay where / is the viscosity of the fluid. Multiplying the continuity equation by u,, adding the result to the momentum equation, and combining terms yields Spu (puuj +3p jP ) (1.19) at axi which upon adding and subtracting the term co2ap x,, becomes apu p 2 u+2 (1.20) at 0a, x x, where Tj puu +3,[(ppo) c(p p)] C (1.21) By differentiating the continuity equation with respect with t, and subtracting the divergence of Eq.(1.20), Lighthill's equation is obtained p'c V2 p' (1.22) at2 0 axX Equation (1.22) has the same form as the wave equation that governs the acoustic field produced by a quadrupole source a2T, /cxOx in a nonmoving medium (Goldstein, 1976). Hence, there is an exact analogy between the density fluctuation that occurs in any real flow and the small amplitude density fluctuations that would result from a quadrupole source distribution (of strength T,,) in a fictitious nonmovingg) acoustic medium with sound speed co. This source will vanish in the region outside certain types of turbulent flows such that Eq. (1.22) reduces to the homogeneous wave equation in such regions. In the case where the acoustics generated are not due to a jet with high temperature and the flow is isentropic, the second term of Lighthill's turbulence stress in Eq.(1.21) can be neglected. For a low Mach number, the third term in Eq.(1.21) can be neglected. It can be neglected because viscosity and heat conduction cause a slow damping due to conversion of acoustic energy into heat and have a significant effect only over large distances. We have therefore shown that T, is approximately equal to pu,uj inside the flow and is equal to zero outside this region. Assuming that the density fluctuation is negligible within the moving fluid, then T, poulUj (1.23) The Lighthill approximation has a great advantage in that it is possible to solve the equation with a standard Green function. But, this approximation is only for isentropic, low Mach number flows. Lighthill's equation could be used to study the sound generated by unsteady flows where there are no solid boundaries (or more correctly, by flows where the effect of boundaries can be neglected). Another limitation of this theory is that the principle of sound generated aerodynamically as stated by Lighthill is relevant only when there is no back reaction of the acoustic waves in flow, such as at the trailing edge or in initial shear layer ( Hirschberg, 1994). In this case, the conversion of mechanical energy into acoustic energy is only one way, and this is the reason why acoustics can be inferred from an incompressible description of the flow. Ffowcs Williams Hawkings equation Let us consider a finite volume of space containing a disturbed flow and rigid bodies in an arbitrary motion with the surrounding fluid being at rest. The bodies and the flow generate a sound. In this case, it is possible to replace both the flow and the surface by an equivalent acoustic source, assuming that the whole medium is perfectly at rest. The key assumption is that no flowacoustics coupling occurs, i.e., the acoustic field does not affect the flow from which the sound originates. Consequently, this approach is not valid when some resonant conditions induce an acoustic feedback on the flow. To represent the real medium with the flow and the obstacle in a convenient way, Ffowcs Williams and Hawkings (1969), and Ffowcs Williams (1969,1992) defined an equivalent medium where the rigid bodies are replaced by mathematical surfaces. In order to preserve the kinematics of the flow and the boundary condition of no crossflow on the surface, a discontinuity must be imposed at the surface location by introducing mass and momentum sources. The mass and momentum equations are written as ap 0 0f + pu = pous,, (f) (1.24) a cl9j \ x In Equation (1.24), po is the mean density, us, is the velocity field of a point on the surface, 5is the Dirac delta function, o', (= (P Po)5) is the viscous stress tensor (P being the static pressure with mean value Po) andf(x,t) = 0 defines the kinematics of the surfaces. If the normal unit vector on the surface is ni, then the boundary condition of no crossflow is simply u.n =u.n (1.25) Ve Real flow field ii V,(t) Fluid at rest S(t) Surface with rigidbody motion Figure 12. Moving surface Ffowcs Williams and Hawkings equation Following the procedure used to obtain Lighthill's equation, we can derive the equation for density variation p' p po a2 2 a2 (1.26) 22 c U2 + "' a (f) +a Pou3(f) (1.26) at2 CO OX2 ,a + The values of p' and T, are zero inside the mathematical surface, because in these zones the fluid is considered to be at rest. Roger (1995) showed that in the presence of flow and rigid bodies, fluctuations in the fluid are exactly the same as those that would exist in an equivalent acoustic medium at rest, and they are forced by three source distribution: * a volume distribution a2T, /xx,axj in the outer region of the surface due to the flow * a surface distribution a/Ox, (o',j (f) f/ ,'\ ) due to the interaction of the flow with the moving bodies * a surface distribution a/at ( po us,,() f/ Ox, ) due to the kinematics of the bodies Like in the Lighthill analogy, the analogy of Ffowcs Williams and Hawkins is limited to flow where there is no flow acoustic coupling, or where the acoustic field does not affect the flow from which the sound originates. The merits of Lighthill and Ffowcs Williams Hawkings analogies The Lighthill analogy marked an important milestone in the study of flow generated acoustics. His theory is able to explain a number of important characteristics of acoustic radiation from a jet. The solution of the Lighthill analogy could be obtained analytically using the Green function. The approach was developed when computation capability was limited. That is why this analogy is a big step for understanding aeroacoustic propagation. For example, his theory gave a mathematical representation of sound and pseudosound, an approximation of the sound emission from subsonic coldair jet flow (which he called selfnoise, i.e., noise generated by turbulentturbulent interaction, and shear noise, i.e., noise generated by turbulent mean shear), for the first time. Lighthill's theory represents the basis for the Ffowcs Williams and Hawkings analogy, which accounts for sound generation along a boundary, as in a helicopter rotor, an airplane propeller, an aircraft engine fan, a compressor, or a turbine. The solution of the second analogy is again solved with the Green's function. In this, it is difficult to obtain a general solution, because the analytical solution is available only for certain boundary shapes. The theory of Ffowcs Williams and Hawkings appeared when the performance of the computers increased but was still substantially less than what a computer is capable of today. Hence, this analogy was very useful for studying sound emission from sources such as a thin strut in a turbulent flow, propeller noise (Gutin, 1948) and sound generation near a plate. The advent of powerful computers has made it easier to study more complex problems using these analogies since no analytical solutions are required. Nonlinear problem: Shock wave formation An example where a linear approximation is not suitable is subsonic flow in a pipe, because waves do not attenuate fast with propagation. In this case a wave can propagate as much as 102 wavelengths before friction has a significant effect. This implies that even a nonlinearity of order 102 in the wave propagation can have a significant effect (Hirschberg and Rienstra, 1994). The most spectacular nonlinear effect is the formation of shock waves as a result of the steepening of the compression side of a wave. The study of the equation of a shock wave is deduced for a homentropic fluid, namely, a fluid with constant entropy (Hamilton, 1998). The characteristic form of the equation is obtained from mass and momentum conservation: S+(u c) a u = 0 (1.27) at ax\ pc which implies that J = ui dp (1.28) pc is invariant along the characteristic path C in the (x, t) plane which are described by the equation dxj =u+c (1.29) CdtC+ When a C wave propagates into a uniform region, the C wave emanating from the uniform region will all carry the same information: J = Jo is constant. Typically for a simple wave, the characteristic C+ are simple lines in the (x, t) plane. Eq.(1.28) and Eq.(1.29) yield u = 1/2 ( + Jo) (1.30) and f = I(J J,) (1.31) pc 2 In other words u and fdp/pc are constant along C+, the thermodynamic variable. Hence, speed of sound can be considered constant, so that (u + c) of C+ is constant. For an ideal gas there is the relation dp 2c (1.32) pc y1 A shock wave appears at the point where two characteristics lines of the same family intersect. The intersection of two C occurs after traveling time t,: t = d(u+c) J (1.33) IL dx which for an ideal gas with constant y can be expressed in terms of the space derivative of pressure at x = 0: t, = P (1.34) Because the variation of the speed of sound is negligible, the distance of appearance of the shock wave is given by (') x, Cot = (1.35) For harmonic waves dp k, p (1.36) dx178 kx, = (1.37) (r2 +1)k Hence, for an amplitude p/p = 102, the shock is expected appear within 10 acoustic waves (Goldstein 1976). Computational Techniques for Aeroacoustics Direct Numerical Simulation (DNS) Direct simulation methods aim to compute both unsteady flow and the sound generated by it. These methods must use a domain that includes the noise producing flow region and at least a part of the nearacousticfield. The computational mesh needs to be selected so that both the flow and its sound can be well represented. The first issue in this approach is that the computational cost of such direct computations is large, hence only simple flow configurations can be studied using this direct method. Direct simulation of the acoustic field solves the compressible NavierStokes equations (or Euler equations in those cases where viscosity is not important) without further approximation. These equations govern the total flow field including the acoustics, so if one could solve them in a domain reaching out to the far field, then the one can obtain the acoustic radiation emerging from the fluctuating flow. Gloerfelt at al. (2003) did a comparison of this method with the acoustic analogy. Both methods gave solutions that agree with the experimental data. This approach is limited to flows where the viscosity is not important; it encounters fundamental difficulties when the Reynolds number is high due to the range of scales present in the flow field. The characteristic frequency of a source that radiates sound is given by the Strouhal number, St = fL/U 0(1), where L and U are the characteristic length and velocity. This implies that the wavelength of sound produced is A= L/M, where M is the Mach number of the flow. On the other hand, the dissipation of turbulent energy takes place at the Kolmogorov length scale r = L/Re3/4 (Tennekes and Lumley 1972), where Re is the Reynolds number of the flow. Thus, the ratio of the wavelength of sound to the size of the energy dissipating eddies, A/qlRe3/4/M, takes high values at low Mach numbers, even for moderate Reynolds numbers. Accordingly, the requirement for spatial and temporal resolutions will simply be beyond computational capabilities in the near future. Thus, although DNS has been utilized quite successfully in low Reynolds number (Gloerfelt et al., 2003; Hu et al., 2002) (Re < 1000), a less direct approach where separate grids can be designed for the viscous and acoustic phenomena appears to be preferable for high Reynolds number flows. Perturbation Technique A perturbation technique consists of splitting up the flow field calculation into two parts. First, a timeindependent viscous background flow is calculated, and then a perturbation equation (that describes the sound) about this background flow is derived, and viscous action on the perturbation is neglected. In this approach, an initial disturbance is introduced upstream which causes an instability and causes the wave to grow resulting in the radiation of sound. Tam and Morris (1980) developed this approach analytically for shear layers, and Tam and Burton (1984) developed it for subsonic jets. Hardin and Pope (1995) developed a slightly different approach to address vorticitydominated flow, an expansion about the timedependent, viscous, incompressible, subsonic flow. If the density p = po is constant, then the continuity and momentum equations become a set of four equations for the three incompressible velocity components and incompressible pressure. One solves this set of equations by using a grid and numerical scheme designed for the viscous aspects of the flow for the time dependent viscous incompressible flow. Using the solution from the previous time step, the compressible flow equation is perturbed about the time dependent incompressible flow. The differences between the compressible and incompressible flows are then assumed to be inviscid and isentropic. Linearized Euler Equation The linearized Euler equation (LEE) comes from the modified Euler equation. The principal issue appears when we deal with LEE with specific source term and LEE with projected source terms (Colonius and Lele, 2004). The numerically difficulty of this approach stems from the fact the full LEE set admit nontrivial instability wave solution of the homogeneous equation. Propagation of linear acoustic waves through a medium with known properties, its refraction and scattering due to the nonuniform of the medium or the base flow, and scattering from solid surfaces, and scattering from solid surfaces (with prescribed boundary conditions) are problem for which computational technique are well suited. Simulation of such phenomena can be based on linearized field equation subject to the physical boundary conditions; refraction and scattering effects are automatically obtained. However in such a direct approach the high frequency limit becomes computationally demanding. Computational Issues Numerical methods for problems of sound generation and propagation must overcome a host difficulties that arise because acoustic waves are very weak compared to nearfield fluctuation, and because they must propagate with little attenuation over long distances. In practice this has dictated the use of highorderaccurate numerical methods. The Numerical Approach to Reduce Dissipation and Dispersion In computational aeroacoustics (CAA), accurate prediction of the generation of sound is demanding due to the requirement of preserving the shape and frequency of wave propagation and generation. Furthermore, the numerical schemes need to handle multiple scales, including long and short waves, and nonlinear governing laws arising from sources such as turbulence, shocks, interaction between fluid flows and elastic structures, and complex geometries. It is well recognized (Hardin and Hussaini, 1992), Tam (Tam and Webb, 1993; Tam et al., 1993) that in order to conduct satisfactory CAA, numerical schemes should induce minimal dispersion and dissipation errors. In general, higherorder schemes are more suitable for CAA than lowerorder schemes since, overall, the former are less dissipative. That is why higherorder spatial discretization schemes have gained considerable interest in computational acoustics (Hixon, 1997; Kim et al., 1997; Lin and Chin, 1995). For longer wavelengths, the formal order of accuracy is sufficient to indicate the performance of a scheme. However, for shorter waves relative to the grid size, it is known that the leading truncation error terms are not good performance indicators (Shyy, 1985; Shyy, 1997). To handle broad band waves, the idea of optimizing the scheme coefficients by minimizing the truncation error associated with a particular range of wave numbers has been used over the years by many researchers, e.g., Hu et al. (1996), Stanescu and Habashi (1998), Shu (1994), Nance et al. (1997), Wang and Sankar(1999), Cheong and Lee (2001), Wang and Chen (2001), Ashcroft and Zang (2003), and Tang and Baeder (1999). A successful approach is the DispersionRelationPreserving (DRP) finite difference scheme proposed by Tam (Tam and Webb, 1993; Tam et al. 1993). The basic idea in the DRP scheme is to optimize the coefficients to satisfactorily resolve short waves with respect to the computational grid, namely, waves with wavelengths of Ax (defined as 68 points per wave or PPW) or shorter. It maximizes the accuracy by matching the wave number and frequency characteristics between the analytical and numerical operators in the range of resolvable scales. Recently, Ashcroft and Zhang (Ashcroft and Zang, 2003) have reported a strategy for developing optimized prefactored compact (OPC) schemes requiring smaller stencil support than DRP. The prefactorization strategy splits the central implicit schemes into forward and backward biased operators. Using Fourier analysis, they have shown that it is possible to select the coefficients of the biased operators such that their dispersion characteristics match those of the original central compact scheme. Hixon and Turkel (1998) proved that the prefactored scheme is equivalent to the initial compact scheme if the real components of the forward and backward operators are equal to those at the corresponding wave number of the original compact scheme, and the imaginary components of the forward and backward operators are equal in magnitude and opposite in sign. Both DRP and OPC schemes are originally designed based on the finite difference approach. In order to satisfy the governing laws of the fluid physics, it may be advantageous to adopt the finite volume approach (Udaykumar et al., 1997; Yang and Ingram, 1999; Udaykumar, 1999), which ensures that fluxes estimated from different sides of the same surface are identical, i.e., no spurious source/sink is generated due to numerical treatment. Furthermore, a finite volume formulation can offer an easier framework to handle the irregular geometry and moving boundaries. In this work, we investigate a finite volume formulation (which we call DRPfv), extending the concept embodied in the original, finite differencebased DRP scheme (which we call DRPfd). Similarly, for the OPC scheme, we extend the basic concepts of the original, finite differencebased OPC (OPCfd) scheme, in a finite volume formulation, named OPCfv. Our ultimate goal is to extend the finite volume version of suitable schemes into a cut cell type of Cartesiangrid computational technique that we have developed earlier for moving and complex boundary computations (Yang and Ingram, 1999; Udaykumar et al., 1999; Kim and Lee, 1996; Ye et al., 1999; Lahur et al., 2000; Calhoun, 2000). Complex Geometry To handle problems of practical interest, a CAA scheme needs to have the capability of handling irregular and curved geometries. It is challenging to develop methods exhibiting desirable order of accuracy and controlling dispersion and dissipation errors while being capable of handling complex geometries. In an attempt to address the need for flexible grid distributions, Cheong and Lee (2001), proposed a gridoptimized dispersionrelationpreserving scheme (GODRP). They considered the approximation of the derivative by I n a au (x+Ax, ) (1.38) 9x Axj > where A = (x x,)/2n (1.39) Ax =(x x)/Ax (i= n,n+l,..., n 1,n) (1.40) /9'1// h /I I / X_n Xn+1 X_ Xo X, X2 Xn1 Xn (=x) Figure 13. Schematic diagram showing the (2n + 1) point stencil on a nonuniform grid The wavenumber of the scheme over a nonuniform grid is obtained by using a Fourier transformation: = n a = eae 'AY(1.41) The scheme is derived by requiring that i) the real part of the scheme closely match the analytical solution in the chosen range of wavenumbers, and ii) the imaginary part of the scheme be as close to zero as possible. These requirements can be achieved by minimizing the integrated error E, defined as E= f ae, x a Ax d(aAx) 0 2 (1.42) +Ax + Sgn(c) exp ln2Ax d(aAx) The terms er and e, denote the upper limits of the integration intervals of the real and imaginary parts, respectively. The term A is the weighting factor, ois the halfwidth of a Gaussian function, and c is the speed of wave propagation in ut + cux = 0. GODRP schemes of curvilinear grids permit an assessment of the accuracy of the finite difference method for curvilinear meshes from the wave number point of view. Through the gridoptimization process, highorder finite difference equations can be solved with curvilinear grids with a guarantee of local and thus resultant global dispersionrelationpreserving properties. Hence, the approach can be used with success for a bodyfitted grid to study the generation of sound around a body with complex geometry. The coefficients are obtained based on local characteristics of the grid; in other words, for a 2n+l point stencil GODRP spatial discretization, the scheme implies a solution of 2n+1 equations for each point in the grid. This will lead to a considerable increase in the computational cost. Tam et al. (2000) proposed to solve the problem around a boundary using ghost points and an optimized extrapolation scheme around the complex boundary. They considered an extrapolation scheme that is of desirable accuracy over a wavenumber range k < aAx < k. The goal for the extrapolation is that the scheme works well for general functions. For this purpose, it is sufficient to consider waves with unit amplitudes over a desired band of wavenumbers fa(x)= eil c+a)] (1.43) The total effect on the functionf(x) will be weighted by the amplitude A(a). The general formula for extrapolation is given by N1 f(xo+7Ax)= ZS f(x), x xojAx (1.44) J=0 where S, (= 0, 1, 2, ..., N1) are the stencil coefficients. Their values are obtained by imposing the following: i) the difference between the left and right sides of Eq.(1.44) is zero when the single Fourier components of Eq.(1.43) is substituted in the formula; ii)the error is zero if the approximated function is zero. The constrained optimization problem is handled by the Lagrange multiplier k N1 N1 L= e"'YISJe dy+ iSJ1 (1.45) 0 J=0 \=0 The boundary curve is approximated by line segments joining the intersection points of the computation mesh and the boundary. For instance, the curved surface between A and B in Figure 14 is replaced by a straight line segment. G2 is a ghost point. A ghost value is assigned to G2 as a boundary condition. The enforcement point is at E and G2E is perpendicular to AB. The value of the derivatives of the function at points A and B are found via extrapolation from the points at (1', 2', 3', 4', 5', 6', 7') and (1, 2, 3, 4, 5, 6, 7), respectively. The method is accurate only for waves with a wavelength of 8 mesh spacings or longer. The method induces numerical instability when the wavenumber is high. Hence, we do not expect that the approach gives accurate results for short waves or high wavenumber. G5 G4 G3 G2 B G 1 '2 3 4 5 6 7 1 // .L* ^ ./ 1' 2' 3' 4' 5' 6' G2 B 1 2 A E 7' 2' a) b) Figure 14. Cartesian boundary treatment of curved wall surface; b) detail around boundary Scope The present thesis has three main contributions. First we investigate two numerical methods for treating convective transport: the dispersionrelationpreservation (DRP) scheme, proposed by Tam and coworkers (Tam and Webb, 1993; Tam et al., 1993), and the spacetime as method, developed by Chang (1995). The purpose is to examine the characteristics of existing schemes capable of handling acoustic problems. The space time ac method directly controls the level of dispersion and dissipation via a free parameter, E, while the DRP scheme minimizes the error by matching the characteristics of the wave. Insight into the dispersive and dissipative aspects in each scheme is gained from analyzing the truncation error. Even though both methods are explicit in time, the appropriate ranges of the CFL number, v, are different between them. Different performance characteristics are observed between the two schemes with regard to long and short waves. Second, the two low dispersion numerical schemes developed via optimization procedures in the wave number domain, namely, the DRP scheme and the optimized prefactored compact (OPC) scheme developed by Ashcroft and Zhang (2003), are extended from their original finite difference framework to the finite volume framework. The purpose is to extend the capability of these schemes to better handle nonlinearity and conservation laws of the fluid motion. Linear and nonlinear wave equations, with and without viscous dissipation, have been adopted as the test problems. By highlighting the principal characteristics of the schemes and utilizing linear and nonlinear wave equations with different wavelengths as test cases, the performance of these approaches is studied. Finally, the Cartesian grid, cutcell method is extended along with the OPCbased finite volume scheme so that this high order scheme can treat curved geometry associated with practical acoustic applications. The approach uses a background Cartesian grid for the majority of the flow domain with special treatment being applied to cells which are cut by solid bodies, thus retaining a boundary conforming capability. Surface integrals around complex geometries are computed using flow properties at the cell surface interpolated from cell nodes while maintaining desirable accuracy level. The rest of the thesis is structured as follows. Chapter 2 presents the investigation of the DispersionRelationPreserving (DRP) scheme, and the spacetime ac method. The characteristics of these two schemes are emphasized using a simple wave equation with the initial disturbance in the form of the Gaussian function Chapter 3 presents principal characteristics and introduces two space discretization schemes: the DRP scheme and the OPC scheme. A low dispersion and low dissipation RungeKutta proposed by Hu and coworkers (1996) is employed for the time stepping procedure, and combined with the DRP and OPC schemes. A study of the dispersive characteristics and stability is presented for these schemes. The boundary treatment is presented in this chapter. The DRP and OPC schemes are then extended to the finite volume approach. Four linear and nonlinear test problems are presented to evaluate the merits of these schemes. Chapter 4 presents the principal characteristics for the Cartesian grid, cutcell approach. Additionally, we present the proposed adjustment for the acoustic approach. Finally, several test problems are presented to demonstrate the performance of the present approach. Finally, in Chapter 5, summary and future perspectives are presented. CHAPTER 2 ASSESSMENT OF DISPERSIONRELATIONPRESERVING AND SPACETIME CE/SE SCHEMES FOR WAVE EQUATIONS Introduction A number of numerical schemes based on various concepts have been proposed to treat wave propagation and convective transport, including the concepts of upwinding, flux splitting, total variation diminishing, monotonicity, nonoscillatory, higher order differencing, and Riemann solvers. For some representative works, see van Leer (1979), Roe (1981), Harten (1983, 1989), Osher and Chakravarthy (1983), Shyy (Shyy, 1983; Shyy et al., 1997), Leonard (1988), Shu and Osher (1988), Hirsch (1990), Liou and Steffen (1993), Tam (Tam and Webb, 1993; Tam et al., 1993), Chang (Chang, 1995; Chang et al., 1999), Thakur et al. (Part I, 1996; Part II, 1996), Toro and Billet (1996), Yu and Chang (1997), Loh et al. (1998, 1999, 2000), Wang and Moin (2000), Oran and Boris(2001). In this chapter, we focus on two approaches: the dispersionrelation  preservation (DRP) scheme based on an optimized highorder finite difference concept, proposed by Tam (Tam and Webb, 1993; Tam et al., 1993), and the spacetime concept, proposed by Chang (Chang, 1995; Chang et al., 1999). These two schemes have derived based on interesting concepts, and have been applied to compute flow problems requiring balanced treatment for both dispersion and dissipation. The well studied wave equation, first order linear hyperbolic equation, is: 0u 0u +c= 0, where a > 0 is a constant (2.1) at x The exact solution of the wave equation for the initial value problem with initial data u(x, 0)= U(x) (2.2) u(x,t)= U(x Ct) ( A simple equation like this one has served useful purpose to guide the development in numerical techniques. The challenge for the numerical solution of the wave equation is to maintain the correct sharpness of the solution without creating wiggles. In other words, the goal is to offer satisfactory resolution, especially for the high wave number components, by balancing numerical dissipation and dispersion. Here we use the term "balancing" because of the finite resolution possessed by numerical computation. The two methods optimized higherorder finite difference (DRP) method (Tam and Webb 1993) and the spacetime conservation element and solution element (CE/SE) (Chang, 1995) address the aforementioned issues from different angles. The philosophy of the DRP method is to maximize the accuracy by matching the wave number and frequency characteristics between the analytical and the numerical operators in the range of resolvable scales. The spacetime CE/SE method views the flux calculations in a joined spacetime conservation element, taking into account the unified wave propagation in a solution element. The present work intends to complement the analyses offered in the original studies to gain further insight into issues related to dispersion, dissipation, and resolution. The main features of both approaches will be highlighted via truncation error analysis and numerical computation to evaluate order of accuracy, stability constraints, and the implication in spatial and temporal resolutions of each scheme. We contrast the performance of the two schemes by varying the ratio between the dominant wavelength and numerical mesh spacing from long ( 20 or higher) to short ( 3). In the following, we will first present the concepts and derivation of both schemes, and then conduct computational evaluations. Some of details presented in Tam and Webb's paper (1993) and in Chang's paper (1995) will be summarized to help discuss the related analysis and numerical assessment conducted in the present work The DispersionRelation Preservation (DRP) Scheme In the following, we summarize the method by Tam and Webb (1993) to highlight the salient features of DRP scheme. Discretization in Space As illustrations, is approximated by the central difference schemes of various orders considering uniform grids with spacing to Ax, we give three examples. * The second order central difference approximation is V 1 ' + O(2) (2.3) * Fourthorder approximation is U l+2 +8u+ 8u1 +2 + O(Ax4) (2.4) Sxi 12Ax * Sixthorder approximation is ( u1+ 9u, +2+45u,1 45u, +9u,2 3 O(6) (2.5) Ix), 60Ax The construction of the scheme (Tam and Webb, 1993) Ou 1N S(x) 1E a u(x+ jA) (2.6) C AX j= N is based on two goals: (i) the behavior of the numerical solution in the resolvable wave number range closely matches that of the exact solution, and (ii) the formal order of accuracy of scheme spanning 2N+ 1 nodes is 2(N 1). To facilitate the derivation, one can start with the Fourier Transform of a function  f(x) and its inverse f (a)= f (x) e adx (2.7) f(x) = f(c)e xda (2.8) as well as the derivative and shift theorems Ofx=i) (a) (2.9) f (x + A) = ef (a) (2.10) With these tools, the relationship between the differential and discrete operators (given in Eq.(2.6)) can be approximated as ia N[ a ej = iaf (2.11) i AX = aeN where N = a e'a (2.12) Axj= N We see that a is a periodic function of aAx with a period of 27i. Obviously the goal is to ensure that a is as close to ac as possible. To accomplish this goal the error is minimized over a certain wavenumber range, aAx [ 7; q] the numerical dispersion is reduced by narrowing the range of optimization (Tam and Webb, 1993) E j c Ax aAx2 d(aAx) (2.13) '7 It is noted that a is real, and hence the coefficients a, must be antisymmetric, i.e., ao = 0 and aj = aj. (2.14) On substituting Eq.(2.12) into Eq.(2.13) and taking Eq.Error! Reference source not found. into account, E can be written as E= k2 a sin (k.j) dk (2.15) where k = cAx Hence, with the antisymmetric condition, one can use the method of the least squares fitting to minimize E, over a wave number range, namely: OE = 0, = 1, 2,...N (2.16) da, Furthermore, to ensure that the scheme is accurate to O(A2(N1)), additional conditions can be made by expanding the right side of Eq.(2.6) as a Taylor series and then equating terms. In this way we have N independent equations with N unknowns: (aj)j=1,N. Note that aj = aj,j = 1,...,N, and ao = 0. Tam (Tam and Webb, 1993) shows the relationship between aAx and rAx for the 4th order, DRP scheme and, the aforementioned 6thorder, 4th order and 2th order central difference schemes over the interval 0 to ;r For ahx up to acAx the individual curves are nearly the same as the straight line aAx = aAx. Thus, the finite difference scheme can provide reasonable approximation for wave number so aAx becomes less than acAx. If we wish to resolve a short wave with straightforward central difference approximations using a fixed size mesh, we need to use a scheme with a large stencil. On the other hand, the 7point DRP scheme has the widest favorable range of hAx. Tam and Webb (1993) demonstrate that DRP can be effective in improving the performance of a given stencil within certain wave number (see Figure 21). In the following, we will only consider the DRP scheme with a symmetric stencil and in particular the 7point formula. 2.5   DRP  6th order 2 ...... 4th order  1.5  0.5  0 0.5 1 1.5 2 2.5 3 aAx Figure 21. aAx versus ahx for the optimized DRP 4th order scheme, 7 point stencil, standard 6th order central scheme and 4th order central scheme To follow wave propagation, it is important to evaluate the group velocity of a finite difference scheme. The group velocity is characterized by da da, which should be almost one in order to reproduce exact result. A way to reduce dispersion is to adjust the range of the wave number in the optimization process. b shows da I da curves of the DRP scheme, for different ranges of the optimization parameter q. Upon examining the corresponding da Ida curves, aAx =1.1 gives the best overall fit. In this case based on da the criterion 1.0 <.01, the optimized scheme has a bandwidth 15% wider than the da standard 6th order central difference scheme. For long wave the important range of aAx is small, and hence any value of 7 less than 1.1 is reasonable. For short waves, on the other hand, dispersion can be noticeable with any choice of r7. Table 21. The stencil coefficient for N = 3 1/ a1 a2 a3 7/2 0.7934 0.1848 0.0254 1.1 0.7688530 0.1650824 0.0204372 0.9 0.762145 0.1597162 0.0190957 0.85 0.7607435 0.1585945 0.0188153 0.7 0.7571267 0.1557013 0.0180920 0.6 0.7551720 0.1541376 0.0177011 a) 1 01 d(7Ax) d (aAx) 1006 1 002 b) Figure 22. Dispersive characteristics of DRP scheme: a) d(adx)/ d (aAx) versus aAx for optimized and standard schemes. .b)d (aAx)/d(aAx) versus aAx for the 4th order optimized scheme for different wave number range of optimization: range is from r to rq Time Discretization Following Tam and Webb (1993) suppose that u(t) is an unknown vector and the time axis is divided into a uniform grid with time step At. To advance to the next level, we use the following 4level finite difference approximation: n+) ~ u +Atb u (2.17) J=0 \at Similar to space discretization, the goal is to develop schemes which closely match the exact solution in the frequency domain, and exhibit desired formal orders of accuracy. In Eq.(2.17) there are four constants, namely: bo, bl, b2, b3. To determine these constants and create a scheme of O(At3) accuracy, the terms in Eq.(2.17) are expanded in a Taylor series to match exactly up to order At2. This leaves one free parameter, bo. The relationship between other coefficients and bo are (Tam and Webb, 1993; Tam et al., 1993): 53 16 23 b,= 3b, + ; b,= 3b, b3 = b + (2.18) 12 3 12 One can utilize the Laplace transform to determine bo. First, the Laplace transform and its inverse of a functionf(t) are related by S(0) = f(t)e'ct'dt; f(t)=if(o))e,tdcod (2.19) where F is a line in the upper half wplane parallel to the real co axis above all poles and singularities. Also, the shifting theorem for Laplace transformation is f (t+ A)=e eAf (2.20) To apply the Laplace transform, we need to generalize the equation to one with a continuous variable, namely, u(t + At)= u(t) + Atyb I( ja) (2.21) J=0 at On applying the shifting theorem to fourlevel scheme presented in Eq.(2.21), we find that ie at = u + At bJejAt u (2.22) Hence, ie'At 1) 0= 3 (2.23) Aty bjeyOAt J=0 To optimize the time stepping scheme, the error is optimized E, = o [Re(At At)]2+(L c)[Im(AA]2 d (At) (2.24) dEi i.e., d= 0. dbo In Eq.(2.24) a is the weight and ; is the nondimensionalized frequency range needed for the numerical scheme to match the exact solution. Substituting Eq.(2.18) and (2.23) into Eq.(2.24) El becomes a function of bo alone, For a= 0.36 and ,= 0.5 the scheme becomes (Tam and Webb, 1993) bo= 2.30256; b,= 2.49100; b2= 1.57434; b3= 0.38589 (2.25) Eq.(2.23) indicates that the relationship between &At and woAt is not one to one. This means that spurious solutions appear. The stability will be established in function of the real solution. It can be written in the form: b3z4 + bz3 +bz2 + b +z 0 (2.26) &At) gAt where z = e l"at and the values ofb, ( = 0, 1, 2, 3) are given by Eq.(2.26). To maintain satisfactory temporal resolution while being stable, the imaginary part of the solution (oAt) should be negative but close to zero. The interval 0 < At < 0.41 satisfies these expectations. Furthermore, Re(oAt) Re(&At) in this range. To ensure that the damping numerical is minimized, Tam and Webb (1993) suggest the condition Im(woAt) < 0.118.104 (i.e., 0 < At<0.19 is adopted). This condition guarantees numerical stability and negligible damping. To compute stability of the scheme we take into consideration FourierLaplace transformation of the wave equation (Eq.(2.1)) ii u = ciau + U tna (2.27) 2;r where a, co* characterize the PDE For the long wave we can approximate wavenumber of the scheme with wavenumber of the PDE a a (2.28) which leads to io*u = cau + kulftal (2.29) Hence o)* = c + kk (2.30) The condition of the numerical stability is that amplification factor for time discretization is less than 1, and hence )* At <0.41. It is also noted from Figure 21 that aAx <1.8 (2.31) hold true. By introducing Eq. (2.31) into (2.30) and upon multiplying by At it is found that o*At <1.8 c [kk.M+1].At (2.32) Ax where M is mach number. Therefore to ensure numerical stability it is sufficient by Eq. (2.32) to restrict At to be less than Atmax, where Atmax is given by 0.41 Ax Atax = (2.33) 1.8[kk.M+1] c Therefore, for At < Atmax the schemes are numerically stable. Consequently, the schemes yield the following criteria for numerical stability: CFL < 0.21 (2.34) Based on Eqs.(2.1), (2.6), and (2.17) one can obtain the final form of the DRP scheme with 7point in space and 4point in time: n+l n cAt 3 ur =u cA bJ Iaku; (2.35) Ax j=0 k= 3 The leading truncation error of the scheme (given in Eq.(2.35)) can be evaluated using the Taylor series expansion, yielding: u,+cu =c j'3b + u AX3 ( 4 (2.36) 2=0 + c y(j4 +2j3)b + + yakk U A)' + O(A)5) 4! 10 5! k I Replacing in Eq.(2.36) by the numerical values of the various coefficients, the scheme becomes: 1t +au = 1.6318a*v3 .Ax3 (1.223a./4 + a 0.012814) Ax4+O(Ax) (23 The following can be summarized in regard to the present DRP scheme: * the scheme is forth order in space, and third order in time * v< 0.21 assures stability and reasonable accuracy * the first term on the RHS of Eq.(2.37) is dissipative. The accuracy bound, v < 0.21, indicates that the coefficient of the leading dissipative term is small. This observation suggests that, depending on the relative magnitude of v and Ax, dispersive patterns may be dominant in numerical solutions. The SpaceTime Conservation Element and Solution Element Method The tenet in this method is to treat local and global flux conservation in a unified space and time domain. To meet this requirement, Chang (1995) introduces solution elements, which are subdomains in the spacetime coordinates. Within each solution element, any flux vector is then approximated in terms of some simple smooth functions. In the last step, the computational domain is divided into conservation elements within which flux conservation is enforced. Note that a solution element generally is not the same as a conservation element. We summarize in the following the key concepts adopted in Chang's approach. aIp Scheme Consider Eq.(2.1), and define F1 = u, F2 = au, x = x, and x2 = t Applying Green theorem we obtain: F2 O dxdy =J(Fdx, + Fdx2) (2.38) r F x, a x a2 S where f(F1dx1 +F2dx2) is contour integral on closed curve S. S Following a vector notation, Eq.(2.38) can be written as (F1,dx, +F2dx2)= I g.d (2.39) s S(V) where = (F,,F2) and ds = (dx1,dx2). (dxl, dx2) is a differential vector associated with a point (xi, x2) on the closed curve S. Because Eq.(2.1) is valid anywhere in the definition domain, we obtain J gcd= 0 (2.40) s(V) where (i) S(V) is the boundary of an arbitrary spacetime region V in E2 (a two dimensional Euclidean space) with xl = x and x2 = t; (ii) f g. d is contour integral on S(V) closed curve S(V); (iii) (F2, F1) is a current density in E2. Note that fI d is the space time flux of (F2, F1) leaving the region V through the surface element. As shown in Figure 23a, Q is a set of mesh points (, n) that is adapted to discretize a physical domain, where j, n = , with i = 0, 1, 2, 3, ... To facilitate the construction of the present scheme a solution element (SE) associated with (/, n) is illustrated in b. For any (x, y) ESE(j, n), u(x, t), is approximated by u*(x, t;j, n): u*(x, t; j, n) = u + (u)x x)+(u,)(t t) (2.41) with u", (u)", and (u,) are constants in SE(j,n), (x,t") are coordinates of the mesh index (j,n), and au *(x,t; j,n) P iu*(x,t; j,n) uI* (x,, t; j, n) = u=(ux) t = (u,) (2.42) 43 a) j3/2 j1 j1/2 j j+1/2 j+1 j+3/2  m t m a fl L _____ _____ L _____ _____ ______ At/2 Ax/2 n+1 n+1/2 n n1/2 n1 (j,n) Ax/2 (j,n) 1/2..... n1/2) / SSE(j 1/2,n1/2) c SE(j1/2,n+1/2) (j,n) A  n1/2) c SE(j +1/2,n1/2) c SE(j +1/2,n +1/2) I . ............ f ) ... ...... ...... I ...... .I SE ( j, n . (j,n) (j,n)L  Figure 23. Scheme of the solution elements (SEs) and conservation elements (CEs): a) The index positions of SEs and CEs; b) SE(j,n); c) CE (j,n); d) CE+(j,n); e) CE+ /12, n+ 12); f) CE 0i+ 12, + /2); Eq.(2.41) has the form of the firstorder Taylor series expansion. Furthermore, if one takes into account that (u,)" = c(uyx, because u*(x, t; j, n) satisfies wave equation, then Eq.(2.41) becomes u (x, t; j,n)= u + (u)"[(x x) c(t t)] (2.43) In each SE(j, n), g(x, t) = (u(x, t), cu(x, t)) is approximated by g* (x, t;, j, n) = (u (x, t; j, n), cu (x, t; j, n)) (2.44) In order to develop appropriate approximation for the flux, one divides the physical domain into nonoverlapping rectangular elements, referred to as conservation elements (CEs). Specifically, CE receive sign "" or "+" in function of the slope of the line that connects the vertex from Q of a CE. If the slope is negative, CE receives the positive sign; otherwise the sign is negative In Figure 23c the line that unifies the vertices founded in Q is positive, and CE index Q is (j,n). The final notation is CE_(j, n). The surface of CE belongs to two different SE, SE(j, n) and SE(j1/2,n1/2). To specify that a part of surface is in a certain SE in Figure 23c, that part is around by a certain type of line. Figure 23d, Figure 23e, and Figure 23f illustrate three other corresponding cases relating the conservation element to the solution element. The approximation of Eq.(2.40) is F+ (j, n) = J gd= 0 (2.45) S(CE+I ,I for all (/, n)eQ. f gd is contour integral over closed path S(CE+(j,n)) and S(CE+I 1 represents the total flux leaving the boundary of any conservation element is zero. The flux at any interface separating two neighboring CEs is calculated using the information from a single SE. As an example, the interface AC depicted in Figure 23c) and Figure 2 3d) is a subset of SE(j, n). Thus the flux at this interface is calculated using the information associated with SE(j, n). We integrate along the entire boundary S(CE+) of CE+ such that CE+ is on the left as we advance in the direction of integration (S(CE+) is traversed counterclockwise). With the above preliminaries, it follows from (39) that: (Ax F+(j,n) [(1 v )(u) (1 v)(ux) / ]+2(1v) ; u2 n) (2.46) cAt where, again, v = is CFL number. Ax Wn it t aan n12 With the aid of Eqs.(2.45) and (2.46), u1 and (uj) can be solved in terms of u"12 and (ux) J if(1 v2) 0, i.e., q(j,n)= Qq(j 1/2,n1/2)+Q q(j+1/2,n1/2) (2.47) where q(jn)= for all (f, n) e Q (2.48) ,n(AX/4)(ux)n QO Q = (2.49) 2 1 (1 v) + (1++v) We have employed the Taylor series expansion to evaluate the leading truncation errors of Eq.(2.46), yielding the following two expressions for F+ and F. respectively: ( cv cV c c 2 +CU = + Ax (Ax4 48 6 16 48v(2.50) (2.50) CV2 CV C C ut +cux = u  _ Ax2 +O(Ax4) S m 48 16 16 48v) From our derivation indicated by Eqs.(2.49), and (2.50) we can establish the following observation for the spacetime CE/SE scheme: * The method is second order in space and time; * The method is stable if v2 < 1 * The dispersion strength of the scheme increases when v (or time step) is reduced because v appears in the denominator in the leading truncation error term in Eq.(2.50). Hence a small value of v may not be desirable. * The third order term O(Ax3) is zero in Eq(2.50), indicating that not only the 2 u/Cx2 term but also the 4U/Cax4 terms are nonpresent. This is the second reason, in addition to previous observation, why the present scheme can be highly dispersive. ae Scheme In order to help to reduce the dispersive nature of the present scheme, Chang (1995) introduces an "artificial viscosity" into original CE/SE scheme. Instead of F+(j,n) = 0, the following modification is made: (1v2)(x)2 F,(j,n)= + 4 (dux)" (2.51) where E is an independent parameter of the numerical variables, and (dux) = ()[(u) 2 +(ux) (2 ( fl/Ax (2.52) Because the magnitude of the added terms in the scheme is controlled by E, the numerical dissipation is controlled by e. However, because F+(j,n)#0 if Es 0, strictly speaking CE+(j,n) and CE.(j,n) are no longer conservation elements in the as scheme. On the other hand, although the net flux entering the interface separating CE+(j,n) and CE.(j,n) is not zero, but the sum ofF+ (,n) and F.(,n) is zero. Hence the total flux leaving CE(j,n) vanishes. As a result, CE(j,n) is a conservation element in the ac scheme. These are seemingly subtle but relevant feature of the present scheme. Using Eqs.(2.47), (2.51), (2.52) one obtains the solution of the ac scheme as: q(j,n)=Mq(j 1/2,n1/2)+M q(j+1/2,n1/2) (2.53) where M =()+v = ( ; M = () (2 (2.54) E\1 21+vj \1 2cl+v ) To assess the leading truncation errors in the ac scheme a modified equation is derived by substituting Taylor series expansion into Eq.(46) for (ux) 1, u The resulting equations are: CV 2 C CV C C C A 2 U, + CuM = U ++  + EE 12 4 4 12v 3v 3) 2 C C C +u.s  +HOT u 6 6 2 (2.55) cv2 cv c c c c 2A x' ut+cux =u\ ++ +E  12 4 4 12v 3v 3) 2 +u.'  +HOT ( 6 6v 2 From Eq.(2.55), we can see that the ac scheme introduce dissipative term uxx and reduce the strength of the dispersive ux. The net outcome term is that the formal order of accuracy is maintained while the scheme becomes less dispersive. Again, it is undesirable to employ a small v, which is different from the DRP scheme discussed previously. As already mentioned, in the ac scheme, both dispersive and dissipative aspects are affected by E. In Refs. (Chang, 1995; Chang et al., 1999), insight is offered in regard to the choice of E and v, based on the von Neumann stability analysis of the ac scheme, without detailed information of the truncation error. The analysis reported in these references suggests that: (i)v should be large; (ii) E shouldn't be close to 0 or 1; (iii) the spurious solution can be effectively suppressed for e = 0.5 in the case of longwavelength. Reference (Loh et al., 1998) suggests that for aeroacoustic computation, it is essential to choose a large CFL number and a small E. Further insight can be gained based on the present truncation error analysis. If we relate sand v( CFL number) by e= v+ e1, then, in Eqs. (2.55): * Coefficient of uxxx in first equation of the system Eq. (2.55) is equal to: V+1L(V 41)2 _EJ (2.56) 3v 4 * Coefficient of uxxx in second equation of system Eq. (2.55) is equal to: vc (v 2 6v 1) (2.57) The above equations indicate that if v is close to 1, and e and v are close to each other, i.e., if E1 is small, then the leading dispersion error shown in Eq.(2.55) is small. On the other hand, to maintain adequate numerical dissipation, it is helpful to let E vary in the same manner as v. These observations as well as the analysis in Chang (1995) indicate that the value of v and E should be coordinated. Further evaluation will be made based on numerical computations, to help establish a more explicit guideline. Numerical Assessment of the DRP and SpaceTime Schemes To assess the individual and relative merits of the DRP and spacetime schemes the following simple test problem is adopted. + = 0. (2.58) at ax The initial condition imposed at time t = 0 are: u = exp In 2  J (2.59) u, = 21n2 exp ln2 which is a Gaussian profile. The exact solution is: U = exp In (2.60) U = 2In 2 exp In2 X Sb ) ( b In this study we evaluate the performance of the schemes in short, intermediate, and long waves relative to the grid spacing, which is assured by the value b/Ax. In all cases Ax =1, and c = 1, so that At v. For the optimized DRP scheme, we adopt the one with 4th level in time and 7point in space, and 4th order formal accuracy in space and 3rd order formal accuracy in time: 3 "n+l)= + v bj (2.61) J o At where v = (because c = 1) Ax 3 (k) a (k) (2.62) S3 The values of the of b, are given in Eq.(24), and the value a, are as chosen to doa minimize difference between and 1 (see Table 21). da For the spacetime aE scheme, the following formulas result: 2 nu 1 2 x 1 2 2 2 +(1 ) n2 (I_ V )2 ( ) 2 2 1 Ax 1n 2 2 4 2 (2.63) + E n n + _In2 +(1 2+(2 (21+)) (u)" J+ ( 4 J+ 2 2 Before evaluating the DRP and spacetime ac method, we further address the influence of the parameters: E and v for the spacetime ac scheme, and v for the DRP scheme. An effort is made in this study to construct a simple guidance appropriate for long as well as short waves. Evidences based on the test problems, aided by the truncation error analysis, indicate that in order to reduce numerical dispersion and to maintain satisfactory resolution, for short wave ( e.g., b/Ax = 3), v and E are preferable to be close to each other (see Figure 24 and Figure 25). For long and intermediate waves there are virtually no need to introduce much numerical dissipation, hence E can be reduced close to zero. On the other hand, mismatched E and v can significantly reduce the accuracy of the scheme. To demonstrate this fact, Figure 24a and Figure 24b show two solutions for the intermediate wave, all with v= 0.5, respectively for v= 0.9. Changing E from 0.5 to 0.99 causes significant numerical dissipation. On the other hand, e = 0 can be a very acceptable choice in for long wave. Figure 25 compiles long and intermediate wave solution with different value of v and E. Taking into consideration the previous observations, it is decided that v= e = 0.99 is a good choice, and is used in the results presented for the spacetime scheme. In more complicated computations involving coupled system with nonlinear effects, such a choice may not be stable. However, for a simple wave equation, this coordinated choice is beneficial. For the present DRP scheme a value of v less than 0.21 guarantees numerical stability and negligible numerical damping. If we decrease the value of v much further, then, as indicated in Eq.(2.37), the numerical damping, as indicated by the leading dissipation term uxx, may be less than adequate due to the small value of v. For these reasons we decided to use v= 0.1 for longer waves and v= 0.2 for short wave (b/Ax = 3). To summarize, the DRP scheme is higher formal order than the spacetime ac scheme. The DRP scheme prefers smaller v (or At), indicating that it is more expensive to compute than the spacetime ac scheme. To evaluate the solution accuracy, we define an error vector as: E= [E,...,EN (2.64) U(x,) is the exact solution at the point x,, and u, is the numerical solution at the point x,. The error norm is adopted as norm one N E =  N where E, =U(x,)u,, 1 which will be used to help to measure the order of accuracy in actual computations. numencal solution 8 =0 10E01 numerical solution =0 10E+00 numerical solution s =0 50E+0O r analytical solution 1 analytical solution r analytical solution 08 08 08 06 06 06 004 04 04 S0 0 02 02 02 ,j , , numerical solution 8=O 70E+00 analytical solution 08 06 S04 02 0 02 160 180 200 220 X) numerical solution S =0 50E+ analytical solution  / 7 \  \ 20 30 40 50 6C x numerical solution s=0 90E+00 analytical solution 04 xx 02 0 02 160 180 200 220 X 00 numerical solution s =0 99E+00 analytical solution 0.6 0.4 0.2 0 , 20 30 40 50 60 X (2.65) (2.66) x numerical solution S =0 99E+00 analytical solution 1 08 06 04 02 0  02 160 180 200 220 X b) Figure 24. Comparison between analytical and numerical solutions Effect of s on the accuracy of spacetime as scheme: a) v= 0.9, b/Ax = 3, t = 200; b) b/Ax = 6, t= 200, v= 0.5 1 0.8 0.6 0.4 0.2 C v=.99 \ \  0.25 0.5 0.75 1 S v0.99 0.25 0.5 0.75 1 v0.99 \ 0.25 0.5 0.75 1 S error v0.5 5.5E02  4.5E02 3.5E02 2.5E02 1.5E02 0.25 0.5 0.75 error v0.5 2.5E02 2.0E02 1.5E02  1 57E03 1.0E02 1 55E03 1 53E03 5.0E03 1 025 05 075 0.25 0.5 0.75 1 S error V=.1 8.0E02 6.0E02  4.0E02 2.0E02 , 0.25 0.5 0.75 1 S error 1.2E01 1.0E01 8.0E02 6.0E02 4.0E02 2.0E02 v0O.1 2 12E03 2 10E03 S208E03 2 06E03 / S204E03 I 02 04 0.25 0.5 0.75 1 S c) Figure 25. The dependence of the error on E for the spacetime as scheme at t = 200: a) b/Ax = 3; b) b/x = 6; c) b/Ax = 20 error error 2.60E03 2.20E03 1.80E03 error v=O.1 error 4.95E04 4.85E04 4.75E04 4.65E04 error 4.09E05 4.08E05 4.07E05 4.06E05 In the following we will present the results based on three categories: * long wave (b/Ax = 20) * intermediate wave (b/Ax = 6) * short wave (b/Ax = 3) They are defined according to the ratio b/Ax, where b is a parameter that characterizes the wavelength. Figure 26 offer an overview of the two schemes for long, intermediate and short wave computation at two different time instants: t/At = 200 and 2000. It is noted that since Ax is fixed in each plot, varying v is the same as varying At. First, we analyze errors for DRP scheme applied to both intermediate and long waves. The following are some of the main observations: (i) the characteristics are similar for short and long times; (ii) the error increases with v: (iii) a slower increase with respect to v, for little. For short wave, DRP scheme behaves differently as time progresses. As previously discussed in the context of truncation error analysis, Eq.(30), a small value of vcan cause insufficient numerical damping. Hence, it is advisable to use a somewhat larger value of v, but still within the stability bound of v< 0.21. The spacetime ac scheme has exhibited unusual error variations as a function of v (or At): the error changes very little with v, and then suddenly decreases rapidly for v close to 1 (see Figure 26 and Figure 27). This trend holds for all resolution. This type of behavior is unusual at first glance, but can be explained from the viewpoint of truncation error analysis. Specifically, with v= e, both ux and ux terms are small as v approach 1, making the scheme sharply improves the apparent order of accuracy. Next, we present the error as a function of spatial resolution (b/Ax). error DRP scheme 16E03 t=200 16E03 b/Ax=3. 1 5E03  1 5E03  1 4E03 1 3E03 1 3E03 al) error 7 7E03 7 6E03 1 j 10.0 01 02 DRP scheme t=2000 b/Ax=3. 01 02 Space time CE SE scheme 3 3E0 23E0 1 3E0 2 5E0 2 2 b/Ax=3. t=200 2 v= 3 102 10 100 error 5 0E04 4 0E04 3 0E04 DRP scheme t=200 b/Ax=6. error DRP scheme 1 OE03 t=2000 90E04 b/Ax=6. 8 0E04[ 1 7 OE04 6 0E04 10.1 01 02 error 50E0 40E0 30E0 20E0 1 OE0 error 90E0 80E0 7 0E0 6 0E0 50E0 Space time CE SE scheme 1 9E02 1 4E02 b/Ax=6. 9 5E03 t=200 V=S 5 0E03 0 5 OE04 . 102 10I 100 V DRP scheme t=200 b/Ax=20. / 4 0 1 02 V SDRP scheme 4 St=2000 0.4 4 b/Ax=20. 4 1 4  0.1 01 02 Space time CE SE scheme 1 9E03 14E03 b/Ax=20. 9 5E04 t=200 V=S 5 OE04 5 0E05 102 10I 100 V DRP scheme DRP scheme error 1 0E03 t=2000 1 E03 t=2000 b/A x=3. 9 OE04 b/Ax=6. 80E04 7 OE04  0.1 6 0E04 01 02 0.4 error 9 0E04  8 0E04  7 0E04 6 0E04  5 0E04 01 02 DRP scheme t=2000 b/Ax=20. 0.4 1/ S0.1 01 02 01 02 b2) Figure 26. The dependence of the error as function of v for short (b/Ax = 3), intermediate (b/Ax = 6) and long (b/Ax = 20) waves for: a) DRP scheme; b) Space time aE scheme, at two different time instants t = 200 and t = 2000. In all cases At = v error 7 7E03 7 6E03 7 6E03 7 6E03 7 5E03 r 4 4 4 E numerical solution S=0.90E+00 numerical solution S=0.99E+00 1  analytical solution 1 analytical solution 0.09 v 0.99 0.8 0.8 S l I Il I\ 0.6 \ 0.6 0.4 0.4 0.2 0.2 180 190 200 210 220 180 190 200 210 220 x X Figure 27. Effect of von the accuracy of space time as scheme: b/Ax = 3, t = 200 Figure 28 depicts the dependence of error for b/Ax between 3 and 300. The graphs show that the spacetime as method is second order accuracy in space. In contrast with the DRP scheme does not exhibit consistent trends for varying spatial resolutions: (i) for short and intermediate waves the order accuracy of this scheme is four; (ii) for long wave, the error tends to become almost independent of space (the slope is close to 0.1). Comparing the performance of these two schemes, we can deduce that the DRP scheme gives a better solution for b/Ax less than 10 (for a short and intermediate wave). For long wave, the spacetime CE/SE scheme gives a better order of accuracy; nevertheless both schemes are satisfactory in practical terms. Finally, we compare these two methods with respect to their performance with sufficient accumulation of the time steps. 57 DRP scheme Space Time a s scheme 102 v=0.02 102 v==O0.1 t=200 t=200 g103 4103 104 104 1 \ 0.01 105  I 105 101 b/Ax 102 101 b/Ax 102 DRP scheme Space Time a s scheme 103^ v=0.1 103 t=200 v=s=0.99 t=200 104 104 : 0.1 \ 1 05 1 106 106 \ 101 b/Ax 102 101 b/Ax 102 Figure 28. The behavior of the error in function of the wavelength: comparison between DRP and spacetime aE schemes Short Wave: b/Ax = 3 Figure 29 summarizes various aspects of error accumulation in time, along with selected solution profiles for both schemes. It is apparent that the space time as scheme introduces both dissipation and dispersion errors, while the DRP scheme shows mainly dispersion errors with less dissipation. In addition, the two schemes exhibit different levels and growth rates. For the spacetime aE scheme, at the beginning, the rate of accumulation is slower. Then, around t (103), the dissipation becomes more substantial. Around t _(5103), both dissipation and dispersion characterize the solution, but the overall error growth rate decreases slightly in comparison to the previous state. Initially, the DRP scheme produces lower error level and grows almost linearly with time. After t (103), the error growth rate slows down; the dispersion becomes substantial, but the solution is not smeared. Around t (2.104), dispersion becomes highly visible while dissipation is also present at a modest level. Intermediate Wave: b/Ax = 6 As shown in Figure 210, for this case, the DRP scheme has lower error than the spacetime aE scheme. The spacetime ac scheme exhibits similar growth trends in both cases, intermediate and short wave. Again, both dispersion and dissipation errors are noticeable. The DRP scheme behaves differently between short and intermediate waves. The error growth rate at the beginning is low, and then becomes higher. Throughout the entire computation, only minor dispersion error appears. Long Wave: b/Ax = 20 Figure 211 presents the solutions for long wave computations. Both schemes perform well over the interval of time. But the spacetime ac scheme exhibits noticeably slower rate of error accumulation. Neither dispersion nor dissipation errors are significant. Space Time a e scheme b/A x=3 v=s=0.99 A / D C D numerical solutiont=0 20E+05 analytical solution t=0 20E+05 08 06 Z 0 4  02 0 02 19980 20000 20020 x A numerical solution t=0 10E+03 S analytical solution t=0 10E+03 08 06 504 02 0 80 90 100 110 120 x error 2 6E02 1 8E02 1 OE02 2 0E03 numerical solution t=0 90E+03 S analytical solution t=0 90E+03 08 06 S04 02 80 880 890 900 910 920 DRP scheme b/Ax=3. v=0.2 r=0.85 1 0.65 C numerical solution t=0 53E+04 analytical solution t=0 53E+04 08 06 X04 02 0 5260 5280 x 5300 5320 numerical solution t=0 200E+05 1 analytical solution t=0 200E+05 A numerical solution t=0 100E+03 analytical solution t=0 100E+03 1475 x1500 1525 Figure 29. Accumulation of the aE scheme; b)DRP sct error in time for short wave b/Ax 3; a) spacetime error 3 0E02 2 5E02 2 0E02 1 5E02 1 OE02 1 error 2 5E02 2 0E02 1 5E02 1 OE02 5 0E03 error 5.0E03 3.0E03 Space Time a E scheme b/Ax=6 v= s0.99 B/ y c /A S numerical solution t=0 20E+05 Analytical solution t=O 20E+05 0o "0 I  i  I I I  I I 103 t 104 numerical solution t=0.10E+03 numerical solution t=0.36E+04 analytical solution t=0.10E+03 B analytical solution t=0.36E+04 08 // 06 504 02 0 30 100 120 3575 3600 3625 x x DRP scheme b/Ax=6. v=0.1 r]=0.6 1.OE03  A 103 t A numerical solution t=0 100E+03 B numerical solution t=0 200E+05 analytical solution t=0 100E+03 analytical solution t=0 200E+05 1F Y 1  Z0.4 5 80 90 100 110 x 20 0.2 19980 20000 20020 b) Figure 210. Accumulation of the error in time for intermediate wave b/Ax = 6; a) spacetime aE scheme; b) DRP scheme Space Time a s scheme B b/A x=20 v=s=0.99 1o 1 A r^, _______ 103 ( 104 ty A numerical solution t=0.10E+0 1 analytical solution t=0.10E+03 0.8 0.6 30.4 0.2 50 100 150 x error 4.5E03 3.5E03 2.5E03 1.5E03 5.0E04 1 0.8 0.6 40.4 0.2 C, 19950 20000 x Numerical solution t=0.100E+03 B analytical solution t=0.100E+03 :\ I5 I 10I 125I 1 50 75 100 125 150  numerical solution t=0.200E+05 analytical solution t=0.200E+05 0.8 0.6 0.4 0.2 0 0.21 0' ' 19950 19980 20010 20040 x b) Figure 211. Accumulation of the error in time for long wave b/Ax = 20, a) spacetime aE scheme; b) DRP scheme error 5.8E03 5.6E03 5.4E03 5.2E03 5. OE03 0.12 B numerical solution t=0.20E+05 analytical solution t=0.20E+05 20050 A 0.8 0.6 0.4 0.2 0 Summary and Conclusions It should be noted that he DRP scheme is a multistep method, which requires more boundary conditions and initial data, while the spacetime ac scheme is a onestep method. Combined with the fact that the DRP scheme performs better with smaller v, it can be more expensive to compute than for spacetime ac scheme. Tam and Webb (1993) show by means of the Laplace transform that the DRP scheme can be constructed to use the same number of initial data as the original PDE and the single time step method. For the boundary treatment, the DRP scheme requires additional points outside the computational domain. A combination of ghost points and backward difference operators, based on similar optimization procedures, is employed. The present study is restricted to the investigation of a simple 1D linear wave equation. Obviously, more issues will arise when multidimensional geometry, nonlinearity of the physics, and coupling of the dependent variables need to be considered. Nonuniform grid also creates varying CFL number even with a constant convection speed. Nevertheless, it seems useful to consider, in a welldefined framework, the merits of the two schemes in a welldefined context. In this sense, the present study can be viewed to establish an upper bound of the performance characteristics of the DRP and spacetime CE/SE methods. CHAPTER 3 FINITE VOLUME TREATMENT OF DISPERSIONRELATIONPRESERVING AND OPTIMIZED PREFACTORED COMPACT SCHEMES FOR WAVE PROPAGAION In computational aeroacoustics (CAA) accurate prediction of generation of sound is demanding due the requirement of preserving the shape and frequency of wave propagation and generation. Furthermore, the numerical schemes need to handle multiple scales, including long and short waves, and nonlinear governing laws arising from sources such as turbulence, shocks, interaction between fluid flows and elastic structures, and complex geometries. It is well recognized (Hardin and Hussaini, 1992; Tam and Webb, 1993; Tam et al., 1993) that in order to conduct satisfactory CAA, numerical schemes should induce minimal dispersion and dissipation errors. In general, higher order schemes would be more suitable for CAA than the lowerorder schemes since, overall, the former are less dissipative. That is why higherorder spatial discretization schemes have gained considerable interest in computational acoustics (Hixon, 1997; Kim et al., 1997; Lin an Chin, 1995). Table 31 summarizes several approaches proposed in the literature. For longer wavelengths, the formal order of accuracy is sufficient to indicate the performance of a scheme. However, for shorter waves relative to the grid size, it is known that the leading truncation error terms are not good indicators (Shyy, 1985; Shyy, 1997). To handle broad band waves, the idea of optimizing the scheme coefficients via minimizing the truncation error associated with a particular range of wave numbers has been used over the years by many researchers, e.g., Refs.(Hu et al., 1996; Stanescu and Habashi, 1998; Nance et al., 1997; Wang and Sankar, 1999; Cheong and Lee, 2001; Wang and Chen, 2001; Ashcroft and Zhang, 2003) A successful approach is the DispersionRelationPreserving (DRP) finite difference scheme proposed by Tam and coworkers (Tam and Webb, 1993; Tam et al., 1993). The basic idea in the DRP scheme is to optimize coefficients to satisfactorily resolve short waves with respect to the computational grid, namely, waves with wavelengths of 68Ax (defined as 68 points per wave or PPW) or shorter. It maximizes the accuracy by matching the wave number and frequency characteristics between the analytical and the numerical operators in the range of resolvable scales. Recently, Ashcroft and Zhang (2003) have reported a strategy for developing optimized prefactored compact (OPC) schemes, requiring smaller stencil support than DRP. The prefactorization strategy splits the central implicit schemes into forward and backward biased operators. Using Fourier analysis, they have shown that it is possible to select the coefficients of the biased operators such that their dispersion characteristics match those of the original central compact scheme. Hixon and Turkel (1998) proved that the "prefactored scheme is equivalent to the initial compact scheme if: i) the real components of forward and backward operators are equal to those at the corresponding wavenumber of the original compact scheme; ii) the imaginary components of the forward and backward operators are equal in magnitude and opposite in sign". Both DRP and OPC schemes are originally designed based on the finite difference approach. In order to satisfy the governing laws of the fluid physics, it can be advantageous to adopt the finite volume approach (Udaykumar et al., 1999; Yang et al., 1999; Udaykumar et al., 1999), which ensures that fluxes estimated from different sides of the same surface are identical, i.e., no spurious source/sink is generated due to numerical treatment. Such a requirement is particularly important when nonlinearity is involved, as is typically the case in shock and turbulence aspects of the aeroacoustic computations. Furthermore, a finite volume formulation can offer an easier framework to handle the irregular geometry and moving boundaries. In this work, we extend the concept embodied in the original, finite differencebased DRP scheme (which we call DRPfd) to a finite volume formulation (which we call DRPfv). Similarly, for the OPC scheme, we extend the basic concepts of the original, finite differencebased OPC (OPC fd) scheme, to a finite volume formulation, called OPCfv. Our overall goal is to develop the finite volume version of DRP and OPC schemes into a cutcell type of Cartesiangrid computational technique that we have developed earlier for moving and complex boundary computations (Yang et al., 1999; Udaykumar et al., 1999; Ye et al., 1999) to treat aeroacoustic problems needed for engineering practices. We present the finite volume formulation of both DRP and OPC schemes, and assess both fd and fv versions of the DRP and OPC schemes, using well defined test problems to facilitate systematic evaluations. Both linear and nonlinear wave equations with different wavelengths and viscous effects are utilized for direct comparisons. In the following, we first summarize the essence of the individual schemes, including derivations, then present assessment of the test cases. Numerical Schemes We use the following onedimensional wave equation to facilitate the development and presentation of the concept and numerical procedures: +c=0 (3.1) at ax The equation contains time and space derivative. In our work the space derivative term is treated with DRP or OPC scheme and the time derivative by a Low Dissipation and LowDispersion RungeKutta (LDDRK) scheme (Hu et al., 1996). The first scheme used for space derivative is DPR scheme. We present the finite volume version of DRP, and the boundary treatment of the DRP. The OPC scheme is the second method considered for the space derivative. The finite difference procedure of the OPC scheme is offered. This presentation is follow by the extension of this approach to a finite volume and the specific boundary treatment of OPC. LDDRK (Hu et al., 1996) scheme is used to approximate time derivative. The principal characteristics will be presented DRP Scheme Finite volumebased DRP scheme (DRPfv) To incorporate the DRPfd concept into a finite volume framework, let us consider a onedimensional linear wave equation: +c = 0 (3.2) at ax To derive the discredized equation, we employ the grid point cluster shown in Figure 31. We focus on the grid point i, which has the grid points i 1, and i+ 1 as its neighbors. The dashed lines define the control volume, and letters e and w denote east and west faces, respectively, of the control volume. For the onedimensional problem under consideration, we assume a unit thickness in they and z directions; thus, we obtain j dx +c(A), (A)w) =0 (3.3) e where (AO)e and (AO)w are the flux across the east and west face, respectively. i2 i1(E) w i e i+l(W) i+2 i+3 Ax Figure 31.Grid points cluster for onedimensional problem Hence, the discretized wave equation (3.2) can be written as Ax +c((A)) (AO)=) 0 (3.4) where q is the averaged value of q over a control volume. Taking into account Eq. (3.4) we describe the general form of the approximation of in 1D using the control volume concept: c ((OA) (OA) (3.5) ax Ax e The general form of the DRP scheme is: ak~+k (3.6) ax A, Ax 3 where Ax is the space grid, and coefficients a, are constant. The DRP scheme has a general form similar to the central difference approximation. Hence, one can adopt a central difference scheme to express Oe in the neighborhood: (0), = a12 + a A21 + 4 40 + 5O+2 + 6+3 (3.7) (&), = 810 3 +M2 + AO 1 +840 + P5i+1 6M +2 (3.8) Taking into consideration Eqs. (3.5), (3.6), (3.7) and, (3.8) we obtain the values of the f,, i = 1,...,6 by imposing that the value of q at the same locations has the same values as that of the DRPfd. Sak l+k =e w (3.9) k=3 Hence, the values of coefficients f s are S= A = a3 / =2 =5 a2 +a3 (3.10) A/3 =/4 =1a +a2 +a3 To illustrate the abovedescribed concept, we consider the following equation: au au au +c +c2 = 0 (3.11) at ax Oy If we integrate Eq. (3.11) on the surface we have (see Figure 32): di7 1 a F (3.12) at Sabcd where the resulting DRPfv scheme is S= [c,(u Ay, +ue Aye +u";Ay, +u c2,(u A, +u' Ax +u" Ax + u ^'Ax)] uj = AZ2,j + P2 1, + 3Uz,, + /4Uz+1,j +1/5+2,j + /6Uz+3,j (3.14) uuj = A u3 + 2U2,j +3U,1,j +4U,,j +5 +6U1, +2,j (3.15) j = P1U,, 2 + 2U,,J 1 +3U,,J +/4U ,J+1 + 5U,,j+2 + ,6U +3 (3.16) u~ = P1,j3 + 2U,j2 + 3,, 1 + 4U,,J + 5U,,J1 + 6U ,,J+2 (3.17) 69 d v + f (c1udy cudx) Sabcd + c1 (UAY + UeAYe + U + uw uAy)+ Sat s at (3.18) + c (uAx, UAx uAx u )Ax) + d N c W w P Je E .L L I I a S b Figure 32. Grid notation for twodimensional problem, where (i) P denotes the center of a cell, (ii) E, W, N, and S denote, respectively, the nodes corresponding to the east, west, north and south neighbors, (iii) e, w, n and s denote, respectively, the center of the east, west, north and south face of the cell, and (iv) a, b, c, and d denote, respectively, the corners of the cell Boundary treatment of the DRP scheme The current version of the DRP scheme requires seven grid points in space. Consequently, it is necessary to impose some supplementary condition for boundary treatments. In this regard, Tam and Webb (1993) devise ghost points. The minimum number of ghost points is equal to the number of boundary conditions. For example, for an inviscid flow the condition of no flux through the wall requires a minimum of one ghost value per boundary point on the wall. It is desirable to use a minimum number of ghost points to maintain simplicity in coding and structuring data. In this work we use only backward difference for grid points near the computational boundary and a ghost point is used only for wall boundary condition. OPC Scheme Finitedifferencebased optimized prefactored compact (OPCfd) scheme To derive the factorized compact scheme Ashcroft and Zhang (2003) define forward and backward operators Df and D,, such that : (D +DF) (3.19) The generic stencil for the forward and backward derivative operators are then defined as: 7F D,+1 +F .DF [aF u + bFu, I+ F .u, + dF . 1 + eu 2] (3.20) fiB D + yB D, = 1[aB u +bB + B .u, +d" u 1 +eB u, (3.21) Y1 +2 U+1 1 2 The coefficients of the scheme are chosen such that: i) the wavenumber of the scheme is close to the important wavenumber of the exact solution; ii) the imaginary components of the forward and backward stencil are equal in magnitude and opposite in sign, and the real components are equal and identical to original compact scheme; iii) the scheme preserves a certain order of accuracy. The authors (Ashcroft and Zhang, 2003) define the integrated error (weighted deviation) as: E = f (a'A Ax) W (aAx) daAx (3.22) where W(aAx) is a weighting function, and r is a factor to determine the optimization range (0< r < 1). The integrated error, defined in Eq.(3.22), is different from the one of Tam and Web (1993) in that it contains the weighting function. The coefficients are obtained by imposing that, within a given asymptotic order, the error is minimal. In space discretization, one sacrifices formal order of accuracy in favor of wideband performance, especially for the short wave components. Finite volumebased OPC scheme (OPCfv) Taking into account Eq.(3.3) that describes the approximation of the first derivative in the finite volume formulation, equations that describe the OPC scheme, Eq.(3.20) and Eq.(3.21), and the idea that the general form of approximation of the function for points at the center of the cell face, namely, e and w assumes similar forms : Fu =0.5(1Fe Be u (3.23) u = 0.5(uF + uB") where uF, u1e, uF and uB" are determined from: U1, + fue = bu,l du, (3.24) q"I + flF =bu du, (3.25) pu, + Be = bu du 1 (3.26) fuI," + rB = bu_ du, (3.27) where the coefficients are the same as those in the OPCfd scheme: 7 = rf= 7, fp= f= f/f, b = bF= d", d= d = bB. These relationships among forward and backward operators are obtained by Ashcroft and Zhang (2003). The boundary treatment of the OPC scheme Boundary Formulation of the OPC scheme employs a biased explicit stencil. Ashcroft and Zhang (2003) design OPCfd scheme with the follow boundary stencil: D, =,, D= eu (3.28) AX = AX j=N3 D 1 iN D eN+ = D SN+j 1(3.29) J=1 j=N3 where the coefficients s, and e, are determined by matching the Taylor series of the forward and backward compact interior stencils to thirdorder accuracy. The boundary treatment in case of OPCfv approach is similar to that of OPCfd, but the boundary stencil is computed on the face: D, = (uA)e (uA)"w (3.30) 3 3 u"I = Zau u = UN, 3 3 ,=1 (3.31) where the value of the coefficients are: B F B 2F  B tF a2 Sl S2 2a eN+eN1 (3.32) a3 S, S2 S3 +eN 1 e2 ~1 =eN F = S r eN eN 1 r2F = 1 S2 (3.33) r3NB = e+eN l eN2 3F = SI S2 S3 Time Discretization The Low Dispersion and Dissipation RungeKutta (LDDRK) Method Hu et al.(1996) consider time integration using the RungeKutta algorithm of the differential equation au = F (u) (3.34) at where the operator F is a function of u. An explicit pstage algorithm advances the solution of Equation (3.34) from the nth to the (n + 1)th iteration as u(O) nu K") =AtF(u() K') = AtFu 1) (3.35) u =u" +bK() i=l,...,p n +1 = (p) where * bp = 1, * u), where p indicates the stage in algorithm advances * un+1, where n indicates the number of iterations for time dependent computation The vale of the un+ can be written on short like n+l P P 0,i n u" "+ b1 i AtJ (3.36) j=1 i= pj+li 7j The resulting algorithm is obtained by optimizing the dispersion and dissipation properties. Assuming F(u) is linear and applying temporal Fourier transform to (3.36), the amplification factor is given by r= u = 1l+ (i At)J (3.37) ]=1 The exact amplification factor is r = eo^At = e (3.38) The numerical amplification factor r in (3.37) is viewed as an approximation of the exact factor. The order of the optimized RungeKutta scheme is indicated by the leading coefficient in (3.37) that matches the Taylor series expansion of e ". For instance, the third order algorithm is obtained by setting y, = 1/j! forj = 1, 2 and 3. To compare the numerical and exact solutions we take into consideration the ratio: r= r e (3.39) re where r represents the dissipation rate (obviously, the correct value should be 1), and 3 represents the phase error (or dispersive error) where the correct value should be zero. Hu et al (1996) obtained coefficients of the low dispersion and dissipation Runge Kutta (LDDRK) scheme by imposing that: i) the scheme has certain order of accuracy, ii) the error of the amplification factor of the scheme is minimized, which means that both dispersion and dissipation errors are minimized. In other words the following integral is minimized: 2 1+ Zcj (io do = min (3.40) j=1 and iii) the amplification factor of the scheme is less than One within the given stability limit 0.04  0.04 LDDRK LDDRK 0.02 LLDDRK =1.85 0.02 L o 0 RL    L i 0.02 0.02 LLDDRK =1.64 RLDDDR = 2.52 0.04 0.04 ______...., 0 1 2 3 o At o*At a) b) Figure 33. Foursixstage optimized RungeKutta of order four scheme: a) dissipation error; b) phase error In this work we use a twostep alternating scheme: in odd steps we use four stages and in the even steps we use six stages. The scheme is a fourthorder accurate scheme in time for a linear problem and secondorder accurate for a nonlinear problem. The advantage of the alternating schemes is that, when two steps are combined, the dispersion and the dissipation errors can be reduced and higher order of accuracy can be maintained. The specific procedure is given bellow. * fourstage K() = AtF(u") K(2) =AtF u" +IK() 4 K(3) = AtF u" + K(2) (3.41) K = AtF u + (3(3.41) Ku + 2 Un+1 n + K(4) * sixstage K(1) = AtF(u) K(2) =AtF(u +0.17667K(1)) K3) =AtF(un +0.38904K(2)) K(4) =AtFfu + K(3) (3.42) K5 =AtF (u K+(4) K(6)= AtF un +1 K(5 n2 K u+l = u + K(6) In the follow we will give an implementation example of LDDRK scheme when we use OPC and DRP scheme for space discretization. Base on Eq.(3.1), the value of F in point / is defined as following: * DRP fd: F c Y ul (3.43) AX =3 * DRPfv. F = c( ul u )/Ax (3.44) where e = u12 + A1 u +3 u/ + 841+1 + 5 u1+2 + 6 U1+3 (3.45) uw = /P ,l3 + ,1,2 + 23 /11 + A / 1 /5u 1+1 /6 u1+2 (3.46) In the linear case, the fv and fd schemes are equivalent. * OPCfd F cDf+DDF) (3.47) where D1 and Df are obtained from the following system of equations: 71 +PDzf= [b(ul +u)+d(uz 1u)], i=1,...,N (3.48) 8DB +r7DB,= [d(u u, u)+b(u ur,j)], i=1,...,N (3.49) where N represent the number of grid points in space. * OPCfv F= c(ulUlw )/Ax (3.50) where e B e u = 0.5(ul B + u ), and u = 0.5(ui" + uw) (3.51) The value of ule, ul B", ul are obtained by solving the follow system of equations 7u' + fue = bu+ du i = 1,..., N (3.52) 7uF + fuFw= bu dul, i 1,...,N (3.53) fuBe + ruB = bu du i = 1,..., N (3.54) fiuBw +rUBw =bu dul i= 1,..., N (3.55) Analytical Assessment of DRP and OPC Schemes Operation Count We will compare the cost between the alternative approaches only for the approximation of the first derivative, because we employ the same time stepping scheme for both scheme. The efficient form of general formula for the discretization in space of the DRPfd scheme is F 2[a3 (x,3 x,3 )+a2 (x+2 x,2)+a, (x,, x, )] (3.56) This scheme requires a total of three multiplications and five additions to evaluate the first derivatives in a certain point. In case of DRPfv the most efficient form of the computations scheme is u, = [A (u, 3+u 2)+, 2(u ~+ +u, 1)+ (u, +u,+)] (3.57) uw = [, (i(u,3 +u,+2)+fi(u 2+u 1)+ s(u, +u )] (3.58) DRPfv requires a greater number of operations than DRPfd: six multiplications and eleven additions to compute the first derivatives at a given point. To see the computational cost of the OPCfd scheme we adopt the most efficient form that is 1DF [b(u,, u,)+d(u,, u,)] D2,F 2 2fAx 2/ (3.59) DB = 1 b(u u 1,)+d(u u+)] r D 2 28Ax 2, 8 (3.60) where the relation between the coefficients of the forward and backward stencils have been substituted to highlight the equivalent terms in the two stencils. The operation count is then four multiplications and five additions per point (Ashcroft and Zhang, 2003). OPCfv can be written in the form: 1DF 1uFe uFw) (3.61) 2 2Ax ' D= (Ue uI)B (3.62) 2 2 2 Ax where uKFe 1 bu, du quFe] P F (3.63) uw = bu, du,, uu] (3.64) u! Bw u= 1 du r u e] In this case the operation count is eleven additions and six multiplications per point. So we can see also, in Table 32 the finite volume approach is computationally more expensive. Dispersion Characteristics The characteristics of the OPC and DRP schemes, in the finite difference form over the interval 0 to r, are shown in Figure 34. One can see that the difference between the effective wave number of the scheme and the real wave is maintained to be within 2% if cAx < 1.30 for the DRP scheme, and cAx < 1.84 for the OPC scheme. The dispersive characteristics of these schemes can be more clearly seen in Figure 35, which shows the phase speed error, abs( x 1 as a function of wave number on a logarithmetic scale. W e see that the daAx ) DRP scheme has a somewhat larger error than the OPC scheme until around 371/4. The error is maintained to be within 2% for aAx less than 0.85 for the DRP scheme, and less than 1.53 for the OPC scheme. Overall, the OPC scheme yields slightly less dispersion error than the DRP scheme. The dispersive characteristics of LDDRK are obtained by studying the value of Irl and 5, i.e., dissipation rate and dispersion error (see Eq.(3.39)), respectively. In Figure 3 3 we can see conditions of stability: Irl < 1 for o)*At < 2.52 .To obtain an accurate solution the dispersive characteristics (Irl and 6) should be close to the exact solution (Irl close to one and close to zero). Hu et al. (1996) considered time accurate criterion Irl 1 < 0.001 (i.e o)*At < 1.64), and 3< 0.001 (i.e., o)*At < 1.85). These two conditions are satisfied if o)*At < 1.64. Stability The FourierLaplace transformation of the wave equation (Eq.(3.1)) is i10*U = ciau + Umntal (3.65) 27z where a, co characterize the PDE For the long wave we can approximate wavenumber of the scheme with wavenumber of the PDE a Ua (3.66) which leads to io*u = cau + kuinial (3.67) 80 Hence o) =cc+kk (3.68) 3.5rLLIllllll71jn==========LLLLLllllj 3 5 1 i 1 1 1 1 1 1 1 1 1 i 1 i 1 1 1i 1i 1i 1i 1i 1 1 1 1 ii i Ideal 2 .5 .. Id l 1 11, 1^ 1* 1 1 1 2.5 t DRP  2 rT 2. ~ ~^~rv T~r ~^~ tz~^~ \~ ~ ~ ~~^ ~ ~ , r r r r r rr iL LL II JJi LLLLLLIIII 1.5 . 0.5 . 0 0 0.5 1 1.5 2 2.5 3 3.5 Figure 34. Dispersive characteristics of the schemes:. dAx versus aAx 10 S     100 : 10 3 I 10 === .. .. ====S == = 0 05 1 15 2 25 3 35 Figure 35. Phase speed error on a logarithmic scale Figure 34 that  DR <1.8 for DRP scheme ''A 2.1 for OPC scheme (3 00 13__  __ __ _____1____3 101 E i i i === === = = < ''  r S y^/  r i 104 3 that a Ax < 1.8 for DRP scheme S_< 2.1 for O scheme (3.69) 1Figure34 that hold true. By introducing Eq. (3.69) into (3.68) and upon multiplying by At it is found that a*At <1.8 c [kk.M+1].At forDRP scheme (3.70) *At< 2.1 [kk.M+1].At for OPC scheme Ax where M is mach number. From Figure 33a it is clear that the condition of stability is satisfied if Ico*At is less than 2.52. Therefore to ensure numerical stability it is sufficient by Eq. (3.70) to restrict At to be less than Atmax, where Atmax is given by 2.52 Ax Atmax 2.52 AY for DRP scheme max 1.8[kk.M +1] c (3.71) 2.52 Ax Atm 2.2 for OPC scheme x 2.1[kk.M+1] c Therefore, for At < Atmax the schemes are numerically stable. Consequently, the schemes yield the following criteria for numerical stability: CFL <1.4 for DRP scheme CFL <1.2 for OPC scheme (3.72) Although it is clear that CFL<1.4 is the stability condition for DRP scheme, this limit does not assure accuracy of the solution. In the previous analysis we have established that the solution is time accurate for 46 LDRRK if rl < 0.001 and 6 8 < 0.001. But this limit is not fixed, but depends on the scheme that is used for space discretization. For example, in the case of the DRP scheme, the solution is considered time accurate as long as rl1 < 0.02 and 6 1 < 0.02, or 0a*At <2.0. Hence, in this case the condition of being both accurate and stable is CFL < 1.1 The OPC scheme is less sensitive to the dispersive characteristics of the LDDRK scheme; hence CFL<1.2 is a condition of the stability and accuracy for the OPC scheme. This limit is in concordance with the stability analysis of Ashcroft and Zhang (2003). Computational Assessment of the DRP and OPC Schemes To investigate the behavior of the schemes, we will use four test problems. First, we consider a onedimensional wave equation with constant speed. The purpose of this test is to check the accuracy, stability, dissipation and dispersion of the scheme. The second test problem is onedimensional nonlinear wave equation with no viscous dissipation. The purpose of this test case is to i) check the influence of singularities on the performance of the scheme, and ii) analyze dispersion properties when waves are coupled. In the third test problem, we consider the onedimensional viscous Burgers equation, which contains unsteady, nonlinear convection and viscous terms. In this case we pay attention to the influence of the viscosity on the solution accuracy. The last test problem is a 2D acoustic scattering problem from the second CAA Workshop (Tam and Hardin, 1997). This problem tests the curved wall boundary and the capacity of the scheme to reproduce different wavelengths. Test problem 1: OneDimensional Linear Wave Equation To assess the behavior of the DRP and OPC schemes the following simple test problem is studied first. du du +c= 0. (3.73) at x u= exp ln2 X ; at t= 0 (3.74) which is a Gaussian profile. This is one of the test problems offered in the second CAA Workshop (Tam and Hardin 1997)] The exact solution is U= exp ln2 x ct (3.75) In this study we evaluate the performance of the schemes in short, intermediate, and long waves relative to the grid spacing, which is assured by the value r/Ax. For time discretization, we previously presented the detailed formulas for the 46 LDDRK, see equations (3.41), and (3.42). Tam (Tam and Webb,1993; Tam et al., 1993) show that dAx is related to aAx, and in function of aAx they divided the wave spectrum into two categories; i) the long waves (waves for which aAx, in this case aAx is less than aAxc, ii) the short waves (waves for which a is not close to a ). This difference between long and short waves is totally dependent upon the grid space. Hence, by inspecting the number of grid points on the wavelength, we can decide that we have a certain category of wave. In the following, we will present the results based on three categories * long wave (r/Ax = 20) * intermediate wave (r/Ax = 6) * short wave (r/Ax = 3) The categories are defined according to the ratio r/Ax, where r is a parameter that characterizes the wavelength of this problem. This test problem is linear; hence we do not expect differences between finite difference and finite volume approach. In regard to the time step selection, the CFL number (v) limit is similar for all schemes. We can see from Figure 36 that the critical CFL number of both schemes is close to 1.1. From the study of the error in time for linear equations with constant convection speed it is clear that the DRPfv and OPCfv schemes have essentially the same behavior as the corresponding finite difference approach; hence, we only present comparison for DRPfv and OPCfv schemes. The error decreases when the grid size in space decreases until a critical value is reached. For all schemes the errors have slopes consistent with the formal order of accuracy in space. This conclusion is confirmed in Figure 37, where the CFL number is maintained at 0.5. For the long time scale solution, the accumulation of error for both DRP and OPC schemes is very close (as seen in Figure 37b). Here, we consider: i) different grid space, so both schemes have almost the same initial error, and ii) the same CFL number (0.5). This behavior is expected because both schemes present the same discretization in time. a DRP fd 102 DRP fv N 103 aDRP fd  OPC fd  DRP fv OPC fv e OPC fd B  B *OPC f a Q 104 S 3 0 10  10 4 I I e I,,,,, S0.1 0.20.0.4 0.05 0.1 At At a) r/Ax = 3 (short wave); b) r/Ax = 10 long wave Figure 36. Errors with respect to the time step size under a fixed space Ax, at t= 50  linear wave equation 2 "" U.UU4 DRP 10 __ DRPfd e OPCfPC 10 OPCfd 10 0.002 S10 0 O 10_4 O S0.001 0 w 1 106 1 107 102 101 100 1000 2000 AX t a) b) Figure 37. Errors under a fixed CFL = 0.5, at t = 50 linear wave equation: a) error with respect to the space size; b) accumulation of the error in time Test problem 2: OneDimensional Nonlinear Wave Equation The finite volume and finite difference schemes are equivalent for a linear equation. The difference between them appears for the nonlinear convective equation. To observe the merits and similarities of DRP, and OPC schemes, we restrict ourselves to the 1D case. In this test, a nonlinear wave equation with a different speed is solved += 0 (3.76) at ax This equation is solved in the conservative form 0u (u2) +0.5 =0 (3.77) 8t 8x To better understand the effect of high gradients and discontinuities, we chose the following initial conditions 0 x0 (3.78) u(x, 0) = >0 (3.78) [1x>0 The solution for this problem can be written 0 x<0 u(x, t)= 0 < x < t (3.79) 1 x>t In this case, for both DRP and OPC schemes, the finite difference version behaves differently from the finite volume version. In the Eqs.(3.41) and (3.42) the function F takes the form: * DRPfd 3 2 F =0.5 ak (u+k) (3.80) k=3 * DRPfv I= 0.5((U)2 (U )2) (3.81) where ue and u, are as defined before * OPCfd F =0.25(D +DF) (3.82) where D1B and DF is backward and forward derivative of u2 in place of u * OPCfv S= 0.25 e)2 )2) (3.83) where ue and u, are defined by (3.52) (3.55) The similarities and differences for all three categories (short waves [Ax/U = 1.0], intermediate waves [Ax/U = 0.25], and long waves [Ax/U = 0.06]) are first presented. It should be noted again that the short, intermediate and long waves are defined based on the numerical resolution. Here, U is defined as the jump (max m,,,); in our case U= 1, hence in the following we discuss only the effect of the grid space step (Ax). 