Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0011602/00001
## Material Information- Title:
- A Finite Volume, Cartesian Grid Method for Computational Aeroacoustics
- Creator:
- POPESCU, MIHAELA (
*Author, Primary*) - Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Acoustics ( jstor )
Aeroacoustics ( jstor ) Boundary conditions ( jstor ) Cartesianism ( jstor ) Error rates ( jstor ) Pistons ( jstor ) Sound ( jstor ) Sound waves ( jstor ) Spacetime ( jstor ) Wave equations ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Mihaela Popescu. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 7/30/2007
- Resource Identifier:
- 74493182 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads:
popescu_m ( .pdf )
popescu_m_Page_062.txt popescu_m_Page_017.txt popescu_m_Page_002.txt popescu_m_Page_043.txt popescu_m_Page_094.txt popescu_m_Page_078.txt popescu_m_Page_127.txt popescu_m_Page_146.txt popescu_m_Page_038.txt popescu_m_Page_145.txt popescu_m_Page_110.txt popescu_m_Page_134.txt popescu_m_Page_114.txt popescu_m_Page_130.txt popescu_m_Page_103.txt popescu_m_Page_142.txt popescu_m_Page_023.txt popescu_m_Page_039.txt popescu_m_Page_086.txt popescu_m_Page_029.txt popescu_m_Page_047.txt popescu_m_Page_032.txt popescu_m_Page_107.txt popescu_m_Page_042.txt popescu_m_Page_074.txt popescu_m_Page_030.txt popescu_m_Page_085.txt popescu_m_Page_125.txt popescu_m_Page_102.txt popescu_m_Page_020.txt popescu_m_Page_061.txt popescu_m_Page_034.txt popescu_m_Page_079.txt popescu_m_Page_058.txt popescu_m_Page_026.txt popescu_m_Page_093.txt popescu_m_Page_052.txt popescu_m_Page_108.txt popescu_m_Page_087.txt popescu_m_Page_013.txt popescu_m_Page_083.txt popescu_m_Page_099.txt popescu_m_Page_076.txt popescu_m_Page_081.txt popescu_m_Page_049.txt popescu_m_Page_112.txt popescu_m_Page_008.txt popescu_m_Page_025.txt popescu_m_Page_141.txt popescu_m_Page_090.txt popescu_m_Page_091.txt popescu_m_Page_096.txt popescu_m_Page_139.txt popescu_m_Page_133.txt popescu_m_Page_003.txt popescu_m_Page_124.txt popescu_m_Page_022.txt popescu_m_Page_121.txt popescu_m_Page_101.txt popescu_m_Page_006.txt popescu_m_Page_095.txt popescu_m_Page_147.txt popescu_m_Page_009.txt popescu_m_Page_054.txt popescu_m_Page_123.txt popescu_m_Page_069.txt popescu_m_Page_016.txt popescu_m_Page_119.txt popescu_m_Page_109.txt popescu_m_Page_118.txt popescu_m_Page_024.txt popescu_m_Page_073.txt popescu_m_Page_068.txt popescu_m_Page_010.txt popescu_m_Page_082.txt popescu_m_Page_129.txt popescu_m_Page_048.txt popescu_m_Page_007.txt popescu_m_Page_070.txt popescu_m_Page_105.txt popescu_m_Page_075.txt popescu_m_Page_046.txt popescu_m_Page_132.txt popescu_m_Page_036.txt popescu_m_Page_092.txt popescu_m_Page_113.txt popescu_m_Page_018.txt popescu_m_Page_106.txt popescu_m_Page_011.txt popescu_m_Page_072.txt popescu_m_Page_064.txt popescu_m_Page_135.txt popescu_m_Page_028.txt popescu_m_Page_137.txt popescu_m_Page_122.txt popescu_m_Page_089.txt popescu_m_Page_019.txt popescu_m_Page_012.txt popescu_m_Page_104.txt popescu_m_Page_098.txt popescu_m_Page_067.txt popescu_m_Page_037.txt popescu_m_Page_051.txt popescu_m_Page_117.txt popescu_m_Page_050.txt popescu_m_Page_027.txt popescu_m_Page_005.txt popescu_m_Page_063.txt popescu_m_Page_041.txt popescu_m_Page_035.txt popescu_m_Page_056.txt popescu_m_Page_120.txt popescu_m_Page_015.txt popescu_m_Page_084.txt popescu_m_Page_065.txt popescu_m_Page_148.txt popescu_m_Page_071.txt popescu_m_Page_144.txt popescu_m_Page_136.txt popescu_m_Page_033.txt popescu_m_Page_055.txt popescu_m_Page_066.txt popescu_m_Page_040.txt popescu_m_Page_014.txt popescu_m_Page_021.txt popescu_m_Page_111.txt popescu_m_Page_128.txt popescu_m_Page_080.txt popescu_m_Page_143.txt popescu_m_Page_001.txt popescu_m_Page_115.txt popescu_m_Page_031.txt popescu_m_Page_057.txt popescu_m_Page_140.txt popescu_m_Page_116.txt popescu_m_Page_059.txt popescu_m_Page_045.txt popescu_m_Page_077.txt popescu_m_Page_004.txt popescu_m_pdf.txt popescu_m_Page_060.txt popescu_m_Page_100.txt popescu_m_Page_044.txt popescu_m_Page_053.txt popescu_m_Page_088.txt popescu_m_Page_126.txt popescu_m_Page_131.txt popescu_m_Page_097.txt popescu_m_Page_138.txt |

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 Nan-Chu 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 DISPERSION-RELATION-PRESERVING AND SPACE- TIME CE/SE SCHEMES FOR WAVE EQUATIONS..............................30 Introduction ....................... ......................................... ................... 30 The Dispersion-Relation Preservation (DRP) Scheme................... ....... .........32 D iscretization in Sp ace ......................... .. ............................................ .. .. ..32 Tim e D iscretization ....................................... ...... ............ .. ............... 37 The Space-Time 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 Space-Time 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 DISPERSION-RELATION- 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 volume-based DRP scheme (DRP-fv) ............................................66 Boundary treatment of the DRP scheme............. ..... .................69 O P C S ch em e ..................... .. .. ...... ...... ............................................ 7 0 Finite-difference-based optimized prefactored compact (OPC-fd) scheme.70 Finite volume-based OPC scheme (OPC-fv) ............................................71 The boundary treatment of the OPC scheme .............................................71 Time Discretization The Low Dispersion and Dissipation Runge-Kutta (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: One-Dimensional Linear Wave Equation .................................82 Test problem 2: One-Dimensional Nonlinear Wave Equation............................85 Test problem 3: One-Dimensional Nonlinear Burgers Equation ......................89 Test problem 4: Two-Dimensional Acoustic Scattering Problem.....................91 Sum m ary and C onclu sions .............................................................. .....................94 4 A FINITE VOLUME-BASED HIGH ORDER CARTESIAN CUT-CELL METHOD FOR COMPUTATIONAL AEROACOUSTICS.........................98 Intro du action ...................................... ................................................ 9 8 Cut-Cell 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 Space-Time CE/SE Scheme...........................................121 Finite-Volume Treatment of Dispersion-Relation-Preserving and Optimized Prefactored Com pact Schem es ....................... ... ........ ...... ............... .... 122 Cartesian Grid, Cut-Cell 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 2-1 The stencil coefficient for N = 3 ..................................................... ............. 36 3-1 A summary of proposed CAA algorithms............................................ ..........96 3-2 The computational cost for DRP and OPC schemes............................ ............97 4-1 Published cut-cell approach for different problems .............. ...... ....................118 LIST OF FIGURES Figure p 1-1 Sound propagation aw ay from a source ........................................ .....................3 1-2 Moving surface Ffowcs Williams and Hawkings equation ...............................15 1-3 Schematic diagram showing the (2n + 1) point stencil on a nonuniform g rid .............................................................................. 2 4 1-4 Cartesian boundary treatment of curved wall surface; b) detail around boundary ..................................... .................. ............... ........... 27 2-1 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 2-2 Dispersive characteristics of DRP scheme........................ ..............36 2-3 Scheme of the solution elements (SEs) and conservation elements (CEs) ..............43 2-4 Comparison between analytical and numerical solutions Effect of E on the accuracy of space-tim e a- scheme ........................................ ............... 52 2-5 The dependence of the error on E for the space-time a-c scheme at t = 200: a) bAx = 3; b) b/Ax = 6; c) b/x = 20............................. .. ............. 53 2-6 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 2-7 Effect of von the accuracy of space time a-c scheme: b/Ax = 3, t = 200 ..............56 2-8 The behavior of the error in function of the wavelength: comparison between DRP and space-tim e a-E schemes ....................................................... 57 2-9 Accumulation of the error in time for short wave b/Ax = 3 ..............................59 2-10 Accumulation of the error in time for intermediate wave b/Ax = 6.......................60 2-11 Accumulation of the error in time for long wave b/Ax = 20..............................61 3-1 Grid points cluster for one-dimensional problem ......................................... 67 3-2 Grid notation for two-dimensional problem ................................. ............... 69 3-3 Four-six-stage optimized Runge-Kutta of order four scheme: a) dissipation error; b) phase error ........................................ ........ ............... 74 3-4 Dispersive characteristics of the schemes ..... ......... ....................................... 80 3-5 Phase speed error on a logarithmic scale ...................................... ............... 80 3-6 Errors with respect to the time step size under a fixed space Ax, at t = 50 linear w ave equation ...................................... ............... .... ....... 84 3-7 Errors under a fixed CFL = 0.5, at t = 50 linear wave equation: .......................85 3-8 Errors with respect to the space step size under a fixed CFL = 0.5, at t = 5; nonlinear wave equation .......................... ................... ............ ... 87 3-9 DRP-fd solution nonlinear wave equation; t = 5; CFL = 0.5 ................................87 3-10 DRP-fv solution nonlinear wave equation; t = 3; CFL = 0.5 ................................88 3-11 OPC-fd solution nonlinear wave equation; t= 5; CFL = 0.5................................88 3-12 OPC-fv solution nonlinear wave equation; t= 5; CFL = 0.5................................88 3-13 Error as a function of Pe nonlinear Burgers equation ................ ...................90 3-14 Numerical solution obtained by DRP schemes nonlinear Burgers q u a tio n ................................................... ..................... ................ 9 1 3-15 Numerical solution obtained by OPC schemes nonlinear Burgers equ atio n ............................................................................ 9 1 3-16 Instantaneous pressure contours at time t = 7; two-dimensional acoustic scattering problem .................. ........................................... ...... 93 3-17 The pressure history at point A, B and C ...................................... ............... 93 4-1 Illustration of the interfacial cells and cut-and-absorption procedures.................. 100 4-2 M modified cut cell approach for CAA..................................... ............... 101 4-3 Radiation from a baffled piston (test problem (ii)) ............................................ 105 4-4 Radiation from a baffled piston: Beam patterns |D(O); ka = 2 ...........................108 4-5 Radiation from a baffled piston: Beam patterns |D(O)|; ka = 7.5 ........................108 4-6 Radiation from a baffled piston: Beam patterns |D(O)|; ka = 12.5 ......................108 4-7 Radiation from a baffled piston: Real of piston radiation impedance....................109 4-8 Radiation from a baffled piston: Numerical solution on axis (ka = 2) ................10 4-9 Radiation from a baffled piston: Comparison between analytical and com puted solutions on axis (ka = 2) ............... ............................... .............. 111 4-10 Radiation from a baffled piston: Contour plot of pressure (ka = 2)....................111 4-11 Radiation from a baffled piston: Comparison between analytical and computed solutions on axis (ka = 7.5) ............. ..... ...................... .......... 112 4-12 Radiation from a baffled piston: Contour plot of pressure (ka = 7.5) ...................112 4-13 Reflection of the pulse on an oblique wall (test problem (ii)): general description ........... .................................................. ............. ...... ... 113 4-14 Reflection of the pulse on an oblique wall (test problem (ii)): cell arou n d b ou n d ary ..................... .. ........................................................ .... 1 14 4-15 Reflection of the pulse on an oblique wall: history of pressure for different angle of w all ................................................. ... ...... .. ........ .... 115 4-16 Reflection of the pulse on an oblique wall: a = 630 ............. ......... ... ............115 4-17 Wave generated by a baffled piston and reflects on an oblique wall: General description of the domain ................... ... ................. .............. 116 4-18 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 dispersion-relation-preservation (DRP) scheme, proposed by Tam and Webb, and the space-time a-c method, developed by Chang. Between them, it seems that for long waves, errors grow slower with the space-time a-c 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 dispersion-relation-preserving (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, cut-cell method is combined with the high-order finite-volume 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 near-acoustic field of a supersonic jet (at about 10 jet-diameters 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 Courant-Friederich-Lewy (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 1-1. 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 = 10-3 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 10-4 and 10-5 meters and the wavelength is about one-third of a meter, so the displacement amplitude is only about 10-4 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 large-scale; 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 non-uniform medium. Specific examples include as: radiation from a duct opening or engine-inlet due to a specified source or a specified duct-mode disturbance, radiation from a specified duct-mode disturbance, and radiation from a specified source across a finite barrier/sound wall with an absorbing treatment. Sound propagation in a specified non-uniform 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 mean-flow 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 small-signal 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 high-speed jets), and sonic boom propagation through atmospheric turbulence, sound propagation in complex fluids and multi-phase systems (such as in lithotripsy), and internal flows of thermo-acoustic 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'= P-Po 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,[(p-po)- 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 flow-acoustics 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 cross-flow 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 cross-flow is simply u.n =u.n (1.25) Ve Real flow field ii V,(t) Fluid at rest S(t) Surface with rigid-body motion Figure 1-2. 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 cold-air jet flow (which he called self-noise, i.e., noise generated by turbulent-turbulent 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 10-2 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 y-1 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 = 10-2, 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 near-acoustic-field. 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 Navier-Stokes 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 time-independent 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 vorticity-dominated flow, an expansion about the time-dependent, 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 non-trivial 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 near-field fluctuation, and because they must propagate with little attenuation over long distances. In practice this has dictated the use of high-order-accurate 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, higher-order schemes are more suitable for CAA than lower-order schemes since, overall, the former are less dissipative. That is why higher-order 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 Dispersion-Relation-Preserving (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 6-8 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 DRP-fv), extending the concept embodied in the original, finite difference-based DRP scheme (which we call DRP-fd). Similarly, for the OPC scheme, we extend the basic concepts of the original, finite difference-based OPC (OPC-fd) scheme, in a finite volume formulation, named OPC-fv. Our ultimate goal is to extend the finite volume version of suitable schemes into a cut- cell type of Cartesian-grid 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 grid-optimized dispersion-relation-preserving 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 Xn-1 Xn (=x) Figure 1-3. 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, o-is the half-width 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 grid-optimization process, high-order finite difference equations can be solved with curvilinear grids with a guarantee of local and thus resultant global dispersion-relation-preserving properties. Hence, the approach can be used with success for a body-fitted 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 N-1 f(xo+7Ax)= ZS f(x), x xo-jAx (1.44) J=0 where S, (= 0, 1, 2, ..., N-1) 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 N-1 N-1 L= e"'Y-ISJe- dy+ iSJ-1 (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 1-4 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 1-4. 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 dispersion-relation-preservation (DRP) scheme, proposed by Tam and coworkers (Tam and Webb, 1993; Tam et al., 1993), and the space-time a-s method, developed by Chang (1995). The purpose is to examine the characteristics of existing schemes capable of handling acoustic problems. The space- time a-c 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, cut-cell method is extended along with the OPC-based 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 Dispersion-Relation-Preserving (DRP) scheme, and the space-time a-c 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 Runge-Kutta 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, cut-cell 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 DISPERSION-RELATION-PRESERVING AND SPACE-TIME 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, non-oscillatory, 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 dispersion-relation - preservation (DRP) scheme based on an optimized high-order finite difference concept, proposed by Tam (Tam and Webb, 1993; Tam et al., 1993), and the space-time 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 higher-order finite difference (DRP) method (Tam and Webb 1993) and the space-time 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 space-time CE/SE method views the flux calculations in a joined space-time 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 Dispersion-Relation 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) * Fourth-order approximation is U -l+2 +8u+ 8u1 +-2 + O(Ax4) (2.4) Sxi 12Ax * Sixth-order 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 = a-e-N 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 anti-symmetric, 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= k-2 a sin (k.j) dk (2.15) where k = cAx Hence, with the anti-symmetric 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(N-1)), 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 6th-order, 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 7-point 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 2-1). In the following, we will only consider the DRP scheme with a symmetric stencil and in particular the 7-point 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 2-1. 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 2-1. 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 2-2. 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 4-level 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 w-plane parallel to the real- co -axis above all poles and singularities. Also, the shifting theorem for Laplace transformation is f (t+ A)=e e-Af (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 four-level 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(A-A]2 d (At) (2.24) dEi i.e., d= 0. dbo In Eq.(2.24) a is the weight and ; is the non-dimensionalized 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.10-4 (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 Fourier-Laplace 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 2-1 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 7-point in space and 4-point 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 +a-u = -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 Space-Time 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 space-time 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. a-Ip 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 space-time 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 2-3a, 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) j-3/2 j-1 j-1/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 n-1/2 n-1 (j,n) Ax/2 (j,n) 1/2..... n-1/2) / SSE(j -1/2,n-1/2) c SE(j-1/2,n+1/2) (j,n) A --- n-1/2) c SE(j +1/2,n-1/2) c SE(j +1/2,n +1/2) I . ............ f ) ... ...... ...... I ......- .I ------SE ( j, n ----.- (j,n) (j,n)L- ---- Figure 2-3. 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 first-order 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 2-3c 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(j-1/2,n-1/2). To specify that a part of surface is in a certain SE in Figure 2-3c, that part is around by a certain type of line. Figure 2-3d, Figure 2-3e, and Figure 2-3f 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 g-d 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 2-3c) 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,n-1/2)+Q q(j+1/2,n-1/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 space-time 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 non-present. This is the second reason, in addition to previous observation, why the present scheme can be highly dispersive. a-e 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: (1-v2)(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 a-s 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 a-c 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 a-c scheme as: q(j,n)=Mq(j- 1/2,n-1/2)+M q(j+1/2,n-1/2) (2.53) where M =()+v = ( ; -M = () (-2 (2.54) E\-1 2-1+vj \1- 2c-l+v ) To assess the leading truncation errors in the a-c 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 -+--+ -- + --E--E 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 a-c 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 a-c 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 a-c 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 long-wavelength. 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 Space-Time Schemes To assess the individual and relative merits of the DRP and space-time 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 7-point 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) S-3 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 2-1). da For the space-time a-E 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 (2-1-+)) (u)" J+ ( 4 J+ 2 2 Before evaluating the DRP and space-time a-c method, we further address the influence of the parameters: E and v for the space-time a-c 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 2-4 and Figure 2-5). 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 2-4a and Figure 2-4b 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 2-5 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 space-time 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 space-time a-c scheme. The DRP scheme prefers smaller v (or At), indicating that it is more expensive to compute than the space-time a-c 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 10E-01 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 2-4. Comparison between analytical and numerical solutions Effect of s on the accuracy of space-time a-s 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 v--0.99 0.25 0.5 0.75 1 v-0.99 \ 0.25 0.5 0.75 1 S error v--0.5 5.5E-02 - 4.5E-02 3.5E-02 2.5E-02 1.5E-02 0.25 0.5 0.75 error v-0.5 2.5E-02 2.0E-02 1.5E-02 - 1 57E-03 1.0E-02- 1 55E-03 1 53E-03 5.0E-03 1 025 05 075 0.25 0.5 0.75 1 S error V=.1 8.0E-02 6.0E-02 - 4.0E-02 2.0E-02 , 0.25 0.5 0.75 1 S error 1.2E-01 1.0E-01 8.0E-02 6.0E-02 4.0E-02 2.0E-02 v0O.1 2 12E-03 2 10E-03 S208E-03 2 06E-03 / S204E-03 I 02 04 0.25 0.5 0.75 1 S c) Figure 2-5. The dependence of the error on E for the space-time a-s scheme at t = 200: a) b/Ax = 3; b) b/x = 6; c) b/Ax = 20 error error 2.60E-03 2.20E-03 1.80E-03 error v=O.1 error 4.95E-04 4.85E-04 4.75E-04 4.65E-04 error 4.09E-05 4.08E-05 4.07E-05 4.06E-05 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 2-6 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 space-time a-c 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 2-6 and Figure 2-7). 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 16E-03 t=200 16E-03- b/Ax=3. 1 5E-03 - 1 5E-03 - 1 4E-03 1 3E-03 1 3E-03 al) error 7 7E-03 7 6E-03 -1 j 10.0 01 02 DRP scheme t=2000 b/Ax=3. 01 02 Space time CE SE scheme 3 3E-0 23E-0 1 3E-0 2 5E-0 2 2 b/Ax=3. t=200 2 -v= 3 102 10 100 error 5 0E-04 4 0E-04 3 0E-04 DRP scheme t=200 b/Ax=6. error DRP scheme 1 OE-03 t=2000 90E-04 b/Ax=6. 8 0E-04[- 1 7 OE-04 6 0E-04 10.1 01 02 error 50E-0 40E-0 30E-0 20E-0 1 OE-0 error 90E-0 80E-0 7 0E-0 6 0E-0 50E-0 Space -time CE SE scheme 1 9E-02 1 4E-02 b/Ax=6. 9 5E-03 t=200 V=S 5 0E-03 0 5 OE-04 . 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 9E-03 14E-03 b/Ax=20. 9 5E-04 t=200 V=S 5 OE-04 5 0E-05 102 10I 100 V DRP scheme DRP scheme error 1 0E-03 t=2000 1 E-03- t=2000 b/A x=3. 9 OE-04 b/Ax=6. 80E-04 7 OE-04 - 0.1 6 0E-04- 01 02 0.4 error 9 0E-04 - 8 0E-04 - 7 0E-04- 6 0E-04 - 5 0E-04 01 02 DRP scheme t=2000 b/Ax=20. 0.4 1/ S0.1 01 02 01 02 b2) Figure 2-6. 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 a-E scheme, at two different time instants t = 200 and t = 2000. In all cases At = v error 7 7E-03 7 6E-03 7 6E-03 7 6E-03 7 5E-03 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 2-7. Effect of von the accuracy of space time a-s scheme: b/Ax = 3, t = 200 Figure 2-8 depicts the dependence of error for b/Ax between 3 and 300. The graphs show that the space-time a-s 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 space-time 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 10-2 v=0.02 10-2 v==O0.1 t=200 t=200 g10-3 410-3 10-4 10-4 1 \ 0.01 10-5 --- I 10-5 101 b/Ax 102 101 b/Ax 102 DRP scheme Space Time a -s scheme 10-3^ v=0.1 10-3- t=200 v=s=0.99 t=200 10-4 10-4 : 0.1-------- \ 1 05 1 10-6- 10-6 \ 101 b/Ax 102 101 b/Ax 102 Figure 2-8. The behavior of the error in function of the wavelength: comparison between DRP and space-time a-E schemes Short Wave: b/Ax = 3 Figure 2-9 summarizes various aspects of error accumulation in time, along with selected solution profiles for both schemes. It is apparent that the space time a-s 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 space-time a-E scheme, at the beginning, the rate of accumulation is slower. Then, around t (103), the dissipation becomes more substantial. Around t _(5-103), 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 2-10, for this case, the DRP scheme has lower error than the space-time a-E scheme. The space-time a-c 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 2-11 presents the solutions for long wave computations. Both schemes perform well over the interval of time. But the space-time a-c 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 solution-t=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 6E-02 1 8E-02 1 OE-02 2 0E-03 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 2-9. Accumulation of the a-E scheme; b)DRP sct error in time for short wave b/Ax 3; a) space-time error 3 0E-02 2 5E-02 2 0E-02 1 5E-02 1 OE-02 1 error 2 5E-02 2 0E-02 1 5E-02 1 OE-02 5 0E-03 error 5.0E-03- 3.0E-03 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.OE-03 - 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 2-10. Accumulation of the error in time for intermediate wave b/Ax = 6; a) space-time a-E 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.5E-03 3.5E-03 2.5E-03 1.5E-03 5.0E-04 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 2-11. Accumulation of the error in time for long wave b/Ax = 20, a) space-time a-E scheme; b) DRP scheme error 5.8E-03 5.6E-03 5.4E-03 5.2E-03 5. OE-03 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 multi-step method, which requires more boundary conditions and initial data, while the space-time a-c scheme is a one-step method. Combined with the fact that the DRP scheme performs better with smaller v, it can be more expensive to compute than for space-time a-c 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 1-D linear wave equation. Obviously, more issues will arise when multi-dimensional 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 well-defined framework, the merits of the two schemes in a well-defined context. In this sense, the present study can be viewed to establish an upper bound of the performance characteristics of the DRP and space-time CE/SE methods. CHAPTER 3 FINITE VOLUME TREATMENT OF DISPERSION-RELATION-PRESERVING AND OPTIMIZED PREFACTORED COMPACT SCHEMES FOR WAVE PROPAGAION In computational aero-acoustics (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 lower-order schemes since, overall, the former are less dissipative. That is why higher-order spatial discretization schemes have gained considerable interest in computational acoustics (Hixon, 1997; Kim et al., 1997; Lin an Chin, 1995). Table 3-1 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 Dispersion-Relation-Preserving (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 6-8Ax (defined as 6-8 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 difference-based DRP scheme (which we call DRP-fd) to a finite volume formulation (which we call DRP-fv). Similarly, for the OPC- scheme, we extend the basic concepts of the original, finite difference-based OPC (OPC- fd) scheme, to a finite volume formulation, called OPC-fv. Our overall goal is to develop the finite volume version of DRP and OPC schemes into a cut-cell type of Cartesian-grid 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 aero-acoustic 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 one-dimensional 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 Low-Dispersion Runge-Kutta (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 volume-based DRP scheme (DRP-fv) To incorporate the DRP-fd concept into a finite volume framework, let us consider a one-dimensional linear wave equation: +c = 0 (3.2) at ax To derive the discredized equation, we employ the grid point cluster shown in Figure 3-1. 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 one-dimensional 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. i-2 i-1(E) w i e i+l(W) i+2 i+3 --Ax Figure 3-1.Grid points cluster for one-dimensional 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 1-D 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), = a1-2 + a A21 + 4 40- + 5O+2 + 6+3 (3.7) (&), = 810 3 +M-2 + 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 DRP-fd. 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 above-described 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 3-2): di7 1 a F (3.12) at Sabcd where the resulting DRP-fv scheme is S= -[c,(u Ay, +ue Aye +u";Ay, +u c2,(u A, +u' Ax +u" Ax + u ^'Ax)] uj = AZ-2,j + P2- 1, + 3Uz,, + /4Uz+1,j +1/5+2,j + /6Uz+3,j (3.14) uuj = A u-3 +- 2U-2,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,j-3 + 2U,j-2 + 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 3-2. Grid notation for two-dimensional 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 Finite-difference-based optimized prefactored compact (OPC-fd) 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 wide-band performance, especially for the short wave components. Finite volume-based OPC scheme (OPC-fv) 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 OPC-fd 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 OPC-fd scheme with the follow boundary stencil: D, =,, D=- eu (3.28) AX = AX j=N-3 D 1 iN D eN+- = -D SN+-j 1(3.29) J=1 j=N-3 where the coefficients s, and e, are determined by matching the Taylor series of the forward and backward compact interior stencils to third-order accuracy. The boundary treatment in case of OPC-fv approach is similar to that of OPC-fd, 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 Runge-Kutta (LDDRK) Method Hu et al.(1996) consider time integration using the Runge-Kutta algorithm of the differential equation au = F (u) (3.34) at where the operator F is a function of u. An explicit p-stage 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= p-j+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 Runge-Kutta 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 3-3. Four-six-stage optimized Runge-Kutta of order four scheme: a) dissipation error; b) phase error In this work we use a two-step alternating scheme: in odd steps we use four stages and in the even steps we use six stages. The scheme is a fourth-order accurate scheme in time for a linear problem and second-order 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. * four-stage 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) * six-stage 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 * DRP-fv. F = -c( ul u )/Ax (3.44) where e = u1-2 + 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. * OPC-fd F -cDf+DDF) (3.47) where D1 and Df are obtained from the following system of equations: 71 +PDzf= [b(ul +-u)+d(uz 1-u)], 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. * OPC-fv 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 DRP-fd 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 DRP-fv 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) DRP-fv requires a greater number of operations than DRP-fd: six multiplications and eleven additions to compute the first derivatives at a given point. To see the computational cost of the OPC-fd 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). OPC-fv 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,-, u-u] (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 3-2 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 3-4. 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 3-5, which shows the phase speed error, abs( x -1 as a function of wave number on a log-arithmetic 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 I|rl -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 Fourier-Laplace 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- -1-1-,- 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 3-4. 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 3-5. Phase speed error on a logarithmic scale Figure 3-4 that- -- DR <1.8 for DRP scheme ''A 2.1 for OPC scheme (3- 00 13__- -- -__ __ _____1____3 10-1 -E -----i i -i --=== === = = < ''--------- ------ ---r---- S------- -y^--/- ----- -----r -i 10-4 3- that a Ax < 1.8 for DRP scheme S_<----- 2.1 for O scheme (3.69) 1Figure3-4 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 3-3a 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 4-6 LDRRK if ||r|-l| < 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 ||r|-l1 < 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 one-dimensional 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 one-dimensional 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 one-dimensional 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: One-Dimensional 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 4-6 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 3-6 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 DRP-fv and OPC-fv schemes have essentially the same behavior as the corresponding finite difference approach; hence, we only present comparison for DRP-fv and OPC-fv 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 3-7, 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 3-7b). 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 10-2 DRP fv N 10-3 -a-DRP -fd -- OPC fd -- DRP fv OPC fv -e- OPC fd B | -B -*--OPC f a Q 10-4 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 3-6. Errors with respect to the time step size under a fixed space Ax, at t= 50 - linear wave equation -2 "" U.UU4 DRP 10 _--_- DRP-fd -e-- OPC-fPC 10- OPC-fd 10 0.002 S10- 0 O 10_4 O S0.001 0 w 1 10-6 1 10-7 10-2 10-1 100 1000 2000 AX t a) b) Figure 3-7. 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: One-Dimensional 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 1-D 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: * DRP-fd 3 2 F =-0.5 ak (u+k) (3.80) k=-3 * DRP-fv I= -0.5((U)2 (U )2) (3.81) where ue and u, are as defined before * OPC-fd F =-0.25(D +DF) (3.82) where D1B and DF is backward and forward derivative of u2 in place of u * OPC-fv 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). |

Full Text |

PAGE 1 A FINITE VOLUME, CARTESIAN GRID METHOD FOR COMPUTATIONAL AEROACOUSTICS By MIHAELA POPESCU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 PAGE 2 Copyright 2005 by Mihaela Popescu PAGE 3 To my daughter Ileana Klein, my parents, my brother and my sister for their love and support PAGE 4 iv ACKNOWLEDGMENTS I would like to express my deepest gratit ude to my advisor Dr. Wei Shyy for his subtle yet effective methods of encouraging all of his student s. I will always value my good fortune to have been one of them. Se ldom does one encounter individuals with intellectual caliber, scie ntific temperament and a spirit of humanity that he embodies. I am equally grateful for his enduring en thusiasm and boundless patience during the process of preparing me to be a researcher and contributor to the global scientific community. Without his enthusiasm and commitme nt to excellence, my research in this field/area could not have been accomplished. I would like also like to express my deep est 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 committ ee members Dr. Lou Cattafesta, Dr. Nam Ho Kim and Dr. Jacob Nan-Chu Chung for th eir 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 gr oup 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 PAGE 5 v close and sincere relationship between me mbers 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 Cara bet, Cosmin Carabet and Lascu Luminita for finding the most extraordinary and creativ e ways of being there for me whenever I needed them. I would like to thank to my family memb ers: my parents, Aurica and Conatantin Popescu; my brothers and sisters, Pavel, Teofil, Marinela, Daniel Popescu, and Tabita Chirita and their family, for their sustaini ng support and love throughout years. They were the undisclosed source of energy and st rength that I relied on during challenging times Finally, I would like to express my deep est love and appreciation to my daughter Ileana Klein, who demonstrated perseverance and resilien ce in face of unfavorable prognosis and the challenging circumstance of her birth. PAGE 6 vi TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES.............................................................................................................ix LIST OF FIGURES.............................................................................................................x ABSTRACT.....................................................................................................................xi ii CHAPTER 1 COMPUTATIONAL ASPECTS IN AEROACOUSTICS..........................................1 Introduction................................................................................................................... 1 Definition of Sound......................................................................................................3 The Characteristics of Sound.................................................................................3 Viscous Effect in Sound Wave..............................................................................5 Classification of Aeroacoustic Problems......................................................................6 Linear Problems in Aeroacoustics.........................................................................6 Nonlinear Problems in Aeroacoustics...................................................................9 Lightills 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 Computational Issues..................................................................................................21 The Numerical Approach to Re duce Dissipation and Dispersion.......................22 Complex Geometry.............................................................................................24 Scope.......................................................................................................................... .27 2 ASSESMENT OF DISPERSION-R ELATION-PRESERVING AND SPACETIME CE/SE SCHEMES FOR WAVE EQUATIONS..............................................30 Introduction.................................................................................................................30 The Dispersion-Relation Preservation (DRP) Scheme...............................................32 Discretization in Space........................................................................................32 PAGE 7 vii Time Discretization.............................................................................................37 The Space-Time Conservation Element and Solution Element Method....................41 a Scheme.........................................................................................................41 aScheme..........................................................................................................46 Numerical Assessment of the DRP and Space-Time Schemes..................................49 Short Wave: b / x = 3..........................................................................................57 Intermediate Wave: b / x = 6..............................................................................58 Long Wave: b / x = 20.........................................................................................58 Summary and Conclusions.........................................................................................62 3 FINITE VOLUME TREATMEN T OF DISPERSION-RELATIONPRESERVING AND OPTIMIZED PREFACTORED COMPACT SCHEMES FOR WAVE PROPAGAION.....................................................................................63 Numerical Schemes....................................................................................................65 DRP Scheme........................................................................................................66 Finite volume-based DRP scheme (DRP-fv)...............................................66 Boundary treatment of the DRP scheme......................................................69 OPC Scheme........................................................................................................70 Finite-difference-based optimized pr efactored compact (OPC-fd) scheme.70 Finite volume-based OPC scheme (OPC-fv)...............................................71 The boundary treatment of the OPC scheme...............................................71 Time Discretization The Low Disp ersion and Dissipation Runge-Kutta (LDDRK) Method............................................................................................72 Analytical Assessment of DRP and OPC Schemes....................................................77 Operation Count..................................................................................................77 Dispersion Characteristics...................................................................................78 Stability................................................................................................................79 Computational Assessment of the DRP and OPC Schemes.......................................82 Test problem 1: One-Dimensional Linear Wave Equation.................................82 Test problem 2: One-Dimensional Nonlinear Wave Equation............................85 Test problem 3: One-Dimensiona l Nonlinear Burgers Equation........................89 Test problem 4: Two-Dimensiona l Acoustic Scattering Problem.......................91 Summary and Conclusions.........................................................................................94 4 A FINITE VOLUME-BASED HIGH ORDER CARTESIAN CUT-CELL METHOD FOR COMPUTATIONAL AEROACOUSTICS.....................................98 Introduction.................................................................................................................98 Cut-Cell Procedure.....................................................................................................99 Test Cases.................................................................................................................104 Radiation from a Baffled Piston........................................................................104 General description....................................................................................104 Directive factor D .......................................................................................107 Pressure on the face of the piston...............................................................109 Low frequency ( ka = 2)..............................................................................110 High frequency ( ka = 7.5)..........................................................................111 PAGE 8 viii Reflection of a Pulse on an Oblique Wall.........................................................112 Wave Generated by a Baffled Piston and Reflected on an Oblique Wall.........115 Conclusion................................................................................................................117 5 SUMMARY AND FUTURE WORK......................................................................119 Assessment of DRP and Space-Time CE/SE Scheme..............................................121 Finite-Volume Treatment of Dispersi on-Relation-Preserving and Optimized Prefactored Compact Schemes............................................................................122 Cartesian Grid, Cut-Cell Approach for Complex Boundary Treatment...................123 Future Work..............................................................................................................124 LIST OF REFERENCES.................................................................................................125 BIOGRAPHICAL SKETCH...........................................................................................134 PAGE 9 ix LIST OF TABLES Table page 2-1 The stencil coefficient for N = 3...............................................................................36 3-1 A summary of proposed CAA algorithms................................................................96 3-2 The computational cost for DRP and OPC schemes................................................97 4-1 Published cut-cell approach for different problems...............................................118 PAGE 10 x LIST OF FIGURES Figure page 1-1 Sound propagation away from a source.....................................................................3 1-2 Moving surface Ffowcs Williams and Hawkings equation...................................15 1-3 Schematic diagram showing the (2 n + 1) point stencil on a nonuniform grid........................................................................................................................... 24 1-4 Cartesian boundary treatment of curv ed wall surface; b) detail around boundary...................................................................................................................27 2-1 x versus x for the optimized DRP 4th order scheme, 7 point stencil, standard 6th order central scheme and 4th order central scheme................35 2-2 Dispersive characteristics of DRP scheme...............................................................36 2-3 Scheme of the solution elements ( SE s) and conservation elements ( CEs )..............43 2-4 Comparison between analytical a nd numerical solutions Effect of on the accuracy of space-time ascheme....................................................................52 2-5 The dependence of the error on for the space-time ascheme at t = 200: a) b / x = 3; b) b / x = 6; c) b / x = 20.........................................................53 2-6 The dependence of the error as function of for short ( b / x = 3), intermediate ( b / x = 6) and long ( b / x = 20) waves................................................55 2-7 Effect of on the accuracy of space time ascheme: b / x = 3, t = 200...............56 2-8 The behavior of the error in f unction of the wavelength: comparison between DRP and space-time a-schemes..............................................................57 2-9 Accumulation of the error in time for short wave b / x = 3...................................59 2-10 Accumulation of the error in time for intermediate wave b / x = 6.......................60 2-11 Accumulation of the erro r in time for long wave b / x = 20..................................61 3-1 Grid points cluster for one-dimensional problem....................................................67 PAGE 11 xi 3-2 Grid notation for two-dimensional problem.............................................................69 3-3 Four-sixstage optimized Runge -Kutta of order four scheme: a) dissipation error; b) phase error...........................................................................74 3-4 Dispersive characteristics of the schemes................................................................80 3-5 Phase speed error on a logarithmic scale.................................................................80 3-6 Errors with respect to th e time step size under a fixed space x at t = 50 linear wave equation..............................................................................................84 3-7 Errors under a fixed CFL = 0.5, at t = 50 linear wave equation:...........................85 3-8 Errors with respect to the space step size under a fixed CFL = 0.5, at t = 5; nonlinear wave equation.................................................................................87 3-9 DRPfd solution nonlinear wave equation; t = 5; CFL = 0.5................................87 3-10 DRPfv solution nonlinear wave equation; t = 3; CFL = 0.5................................88 3-11 OPC-fd solution nonlinear wave equation; t = 5; CFL = 0.5................................88 3-12 OPC-fv solution nonlinear wave equation; t = 5; CFL = 0.5................................88 3-13 Error as a function of Pe nonlinear Burgers equation............................................90 3-14 Numerical solution obtained by DRP schemes nonlinear Burgers quation......................................................................................................................91 3-15 Numerical solution obtained by OPC schemes nonlinear Burgers equation....................................................................................................................91 3-16 Instantaneous pressu re contours at time t = 7; two-dimensional acoustic scattering problem....................................................................................................93 3-17 The pressure history at point A, B and C.................................................................93 4-1 Illustration of the interfacial cel ls and cut-and-absorption procedures..................100 4-2 Modified cut cell approach for CAA..................................................................101 4-3 Radiation from a baffled piston (test problem (ii))...............................................105 4-4 Radiation from a baffled piston: Beam patterns | D ( )|; ka = 2.............................108 4-5 Radiation from a baffled piston: Beam patterns | D ( )|; ka = 7.5..........................108 4-6 Radiation from a baffled piston: Beam patterns | D ( )|; ka = 12.5........................108 PAGE 12 xii 4-7 Radiation from a baffled piston: Real of piston radiation impedance....................109 4-8 Radiation from a baffled pi ston: Numerical solution on axis ( ka = 2) ..................110 4-9 Radiation from a baffled piston: Comparison between analytical and computed solutions on axis ( ka = 2).......................................................................111 4-10 Radiation from a baffl ed piston: Contour pl ot of pressure ( ka = 2)......................111 4-11 Radiation from a baffled piston: Comparison between analytical and computed solutions on axis (ka = 7.5)...................................................................112 4-12 Radiation from a ba ffled piston: Contour plot of pressure ( ka = 7.5)...................112 4-13 Reflection of the pulse on an obliq ue wall (test problem (ii)): general description..............................................................................................................113 4-14 Reflection of the pulse on an obli que wall (test problem (ii)): cell around boundary.....................................................................................................114 4-15 Reflection of the pulse on an obli que wall: history of pressure for different angle of wall............................................................................................115 4-16 Reflection of the pulse on an oblique wall: = 630...............................................115 4-17 Wave generated by a ba ffled piston and reflec ts on an oblique wall: General description of the domain.........................................................................116 4-18 Wave generated by a baffled piston and reflects on an oblique wall: = 630; a) t = 9; b) t = 14;.....................................................................................116 PAGE 13 xiii Abstract of Dissertation Pres ented 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: Mechanic al 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 computa tions, it is essential that the numerical techniques for acous tic 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 fi rst investigate two promising numerical methods for treating convective transport: the disper sion-relation-preservation (DRP) scheme, proposed by Tam and Webb, and the space-time amethod, developed by Chang. Between them, it seems that for long waves, errors grow slower with the space-time ascheme, while for short waves, often cri tical for acoustics co mputations, errors accumulate slower with the DRP scheme. Based on these findings, two optimized PAGE 14 xiv numerical schemes, the dispersion-relationpreserving (DRP) scheme and the optimized prefactored compact (OPC) scheme, origina lly 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, cut-cel l method is combined with the high-order finite-volume schemes to offer additional cap abilities 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 PAGE 15 1 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 extern ally applied forces or motions of classical acoustics. Thus, the sounds generated by vibr ating violin strings and loudspeakers fall into the category of classical acoustics, wher eas the sound generated by turbulent flow or the unsteady aerodynamic forces on propellers fa lls into the domain of aeroacoustics. A main feature of aeroacoustics is the la rge 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 near-acoustic field of a s upersonic jet (at about 10 jet-diameters away) the acoustic disturbance amplitudes are about th ree 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 ca ses, the efficiency can be much smaller; thus the amplitude of disturbance is much less. Another issue in acoustics is the differen ce between the scales of the unsteady flow and sound. This phenomenon is evid ent in situations in which flow speed is much less than the sound speed that characterizes the me dium. This is because the time scale of the flow and the sound must match. In low Mach number flow (M 1) this will give an acoustic wavelength that is M-1 times longer than the flow le ngth scale. In this case, a direct numerical simulation of the dynamic fl ow and generation of s ound field will not be possible. PAGE 16 2 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 in space, as well as a frequency = 2 f in time. These are coupled by the relation f = c where c is the speed of sound in the medi um (assumed quiescent). Thus, to determine an acoustic wave, one must reso lve both its wavelength and its frequency. Numerical spatial and temporal approximati on meets the same problem: the dissipation must be reduced in space, like that in tim e. 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 dissi pation 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. Howe ver, it is very diffi cult to propagate an accurate wave in this kind of grid. The dissipation and disper sion characteristics depend upon the Courant-Friederich-Lewy (CFL) number (= c t / x in one dimension), which can be interpreted as the distance an acousti c 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 refl ected such that it starts propagating back in the other direction. Anis otropy can appear because of the different sized grid spacing along different directions An approach to solving the problem was suggested by Goodrich (1995), who recomme nded approximating the solution of the PAGE 17 3 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 co mputational 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 electroma gnetic waves, as they do not require a medium). Sound is a wave that moves para llel to the directio n of the propagation Figure 1-1. 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 w ith a uniform pressure p0 and density 0. When this is perturbed by a sound wave the pressure at position x and time t changes to p0 + p ( x t ), the density source PAGE 18 4 changes to 0 + ( x t ), and the fluid particle moves with a velocity vxt The ratios 0'/ pp and 0'/ 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 f 20 kHz (1.1) The range described by (1.1) is called th e audible range. The maximum sensitivity of the ear is within the freque ncy range of 2 kHz to 3 kHz (policemans whistle tone) The acoustic pressure of intense rock et engine noise can be as mu ch as 9 orders of magnitude greater than the pressure of the weakest sound detectable by hum an 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 inte nsity level. Although the leve ls are unitless, they are expressed in decibels (dB). SPL = 20 log10( p/pref) (1.2) On this scale a fluctuation of one atmo sphere 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 / p0 = 10-3. At the threshold of pain, flui d 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 wh ich the sound wave is traveling. The displacement amplitude of the particles is thus between 10-4 and 10-5 PAGE 19 5 meters and the wavelength is about one-third of a meter, so the di splacement amplitude is only about 10-4 of the wavelength sound (Dow ling and Ffows Williams, 1983). In the case where sounds can be approximated by small perturbations, several effects are noteworthy. First, there is no inte raction 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 behavi or of the sound numeri cally, we have to use a scheme that adds a low leve l of numerical dissipation. Viscous Effect in Sound Wave The effect of Reynolds number on C AA is ambiguous. Although all sound is ultimately dissipated into heat by viscosity, acoustics is generally considered to be inviscid fluid phenomena. If the viscous te rm is considered in the standard linear analysis, one finds its value to be 0Re2/ c (1.3) where c is the speed of sound, is the wavelength, 0 is the density and is the coefficient of viscosity. The values in air of these parameters for most practical interests are (Blackstock, 2000) PAGE 20 6 0 5331.6 0.33 1.293 1.710 c (1.4) Based on these values, Re is approximately 4.8 107 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 be ing essentially a weak motion of an inviscid fluid. The effect of viscosity can be taken in to consideration after the sound travels about 20/ wavelengths (i.e., after 1.5 x 108 m). Thus, dissipative loss becomes important for high frequency sound propa gation over long distances. On the other hand, if one considers the ge neration 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 sour ce of sound in flows is the acceleration of the vorticity (Powell, 1964), whic h is only present in viscous flow. Even there, a curious independence is observed that is apparently du e to the fact that the large-scale; 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 PAGE 21 7 scattered into another type of disturbance). A linearized solution of the unsteady Euler equations subject to appropriate boundary cond ition can provide useful prediction of the noise generation. The problems of linear inte raction 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 bounda ry value problems of linear acoustics: sound propagation in a uniform media in pr esence of reflecting surfaces, barriers, absorbing walls, and ducts, and also propagatio n and scattering of sound in a prescribed non-uniform medium. Specific examples incl ude as: radiation from a duct opening or engine-inlet due to a specifi ed source or a specified duct-mode disturbance, radiation from a specified duct-mode disturbance, and radiation from a specified source across a finite barrier/sound wall with an absorbing treatment. Sound propagation in a specified non-unifo rm time dependent medium including refraction/scattering in steady and unstea dy vortical flows, sound propagation in nonuniform ducts including the interaction with geometrical changes, linear impedance and mean-flow variations are considered as linear problems of scattering. The equation of a linear acoustics wa ve is deduced for a homogeneous fluid characterized by = 0, P = p0, and 0' uuu The sound waves minutely disturb the status of the quiescent fluid 00 2 000 0, 0',' Ppppc uuuc (1.5) PAGE 22 8 To obtain the wave equation, we will star t with the Euler equa tion and the equation of state for an isentropic process 0 0t tuu uuuP (1.6) and 0pc (1.7) Substituting Eq.(1.5) into (1.6), the following is obtained 0 000 0t ttuuu uuuuuup (1.8) The first order terms are small because th e perturbation is very small; hence, second order terms or higher can be neglecte d. The underlined terms are second order or higher. The result is 0 0 2 00 0t tu up pc (1.9) To reduce the set to a single equation in u the equation of state Eq.(1.7) first is used to eliminate p from the second expression in Eq. (1.9) 2 000tuc (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 22 00ttucu (1.11) PAGE 23 9 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 small-signal 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 s ound beam (acoustics near field of high-speed jets), and sonic boom propagation through atmo spheric turbulence, sound propagation in complex fluids and multi-phase systems (such as in lithotripsy), and internal flows of thermo-acoustic cooling system. Also included in this category are problems of scattering of nonlinear disturbances into s ound, such as airframe noise. The basic equations that describe the nonlin ear behavior are the same that describe the general motion of a viscous, heat conduc ting fluid: mass conservation, momentum conservation, entropy balance, and thermodynamic state. The mass conservation or continuity equation is 0 D u Dt (1.12) where is the mass density, u is the fluid velocity vector, and D/Dt = / t + u is the total, or material, time derivative. The momentum equation may be written as 21 3BDu Puu Dt (1.13) where P is the pressure, is the shear viscosity, and B is the bulk viscosity. Shear viscosity accounts for diffusi on of momentum between adj acent fluid elements having different velocity. Bulk visc osity provides an approximation description, valid at low PAGE 24 10 frequency, of nonechilibrum deviation between actual local pressure and thermodynamic pressure. Equation of state is given by , PPT or PPs (1.14) where T is the absolute temperature, and s is the specific entropy. A commonly used equation of state of a perfect gas is 0 00exp/vP ssc P (1.15) where P0, 0, and s0 are reference values of the pressure, density, a nd specific entropy. = cP/ cv is the ratio of the specific heats at heat pressure ( cp) and constant volume ( cv). In case of arbitrary fluid the equa tion of state is obt ained by expanding of Eq.(1.14) in a Taylor series about (0, s0). A more general description take implies relaxation, like vibr ational relation of diatomic molecules (as in air) and chem ical relation in s eawater (Hamilton and Blackstock, 1997). The former occurs when the energy associated with molecular vibration fails to keep in step with mol ecular translation energy associated with the fluctuating temperature in gas. Lightills acoustic analogy The study of flow that generates acoustic waves began with Gu tins theory (1948) of propeller noise, which was developed in 1937. However, the th eory could not be considered until 1952, when Lighthill (1952, 1954) introduced his acoustic analogy to PAGE 25 11 deal with the problem of calcu lating acoustic radiation from a relatively small region of turbulent flow embedded in an infinite homogeneous fluid. It is well known that a gas (not monoato mic) allows three distinct fundamental modes of oscillation: the vortical and entr opy 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 oscilla tions are small, and ii) the base flow is the uniform medium at rest or is in un iform 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 (195 2, 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 c0 and density 0 is constant. Density and pressure fluctuations are defined by 0 0' p pp (1.16) where the subscript is used to denote cons tant reference values that are usually taken to be corresponding properties at large distance from the flow. Lighthills basic idea was to reformulate th e general equations of gas dynamics in order to derive a wave equation suit able to describe sound propagation. The continuity and momentum equation are 0j j ij i ji jiju tx u p uu txxx (1.17) PAGE 26 12 where ij is the component of the viscous stress tensor. For a Stokestian gas it can be expressed in terms of the velocity gradient 2 3j ik ijij jiku uu x xx (1.18) where is the viscosity of the fluid. Multiplying the continuity equation by ui, adding the result to the momentum equati on, and combining terms yields iijijij juuup tx (1.19) which upon adding and subtracting the term c0 2/ xi, becomes 2 0ij i iiT u c txx (1.20) where 2 000 ijijijijTuuppc (1.21) By differentiating the continuity equation with respect with t, and subtracting the divergence of Eq.(1.20), Light hills equation is obtained 2 2 22 0 2' 'ij ijT c txx (1.22) Equation (1.22) has the same form as th e wave equation that governs the acoustic field produced by a quadrupole source 2Tij/ xi xj in a nonmoving medium (Goldstein, 1976). Hence, there is an exact analo gy 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 Tij) in a fictitious (nonmoving) acoustic medium with sound speed c0. This source will vanish in the region outside certain types PAGE 27 13 of turbulent flows such that Eq. (1.22) re duces to the homogeneous wave equation in such regions. In the case where the acoustic s generated are not due to a jet with high temperature and the flow is isentropic, the second term of Lighthills turbulence stress in Eq.(1.21) can be neglected. For a low Mach number, the third term in Eq.(1.21) can be negl ected. It can be neglected because viscosity a nd heat conduction cause a slow damping due to conversion of acoustic energy into heat and have a si gnificant effect only over large distances. We have therefore shown that Tij is approximately equal to uiuj inside the flow and is equal to zero outside this region. Assuming that the density fluctuation is negligible within the moving fluid, then Tij 0uiuj (1.23) The Lighthill approximation has a great advant age 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. Lighthills 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 Li ghthill is relevant on ly 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 c onversion of mechanical energy into acoustic energy is only one way, and this is the r eason why acoustics can be inferred from an incompressible description of the flow. PAGE 28 14 Ffowcs Williams Hawkings equation Let us consider a finite volume of sp ace containing a distur bed flow and rigid bodies in an arbitrary motion with the surround ing fluid being at rest. The bodies and the flow generate a sound. In this case, it is pos sible to replace both th e flow and the surface by an equivalent acoustic source, assuming that the whole medium is perfectly at rest. The key assumption is that no flow-acous tics coupling occurs, i.e., the acoustic field does not affect the fl ow from which the sound orig inates. Consequently, this approach is not valid when some resonant conditions induce an acoustic feedback on the flow. To represent the real me dium with the flow and the ob stacle in a convenient way, Ffowcs Williams and Hawkings (1969), a nd Ffowcs Williams (1969,1992) defined an equivalent medium where the rigid bodies are replaced by mathema tical surfaces. In order to preserve the kinematics of the fl ow and the boundary condition of no cross-flow on the surface, a discontinuity must be im posed at the surface location by introducing mass and momentum sources. The mass and momentum e quations are written as 0''jsi ji iijijij jjf uuf txx f uuuf txx (1.24) In Equation (1.24), 0 is the mean density, usi is the velocity field of a point on the surface, is the Dirac delta function, ij (= ij (P P0)ij) is the viscous stress tensor (P being the static pressure with mean value P0) and f( x ,t) = 0 defines the kinematics of the surfaces. If the normal un it vector on the surface is n then the boundary condition of no cross-flow is simply PAGE 29 15 sunun (1.25) Figure 1-2. Moving surface Ffowcs Williams and Hawkings equation Following the procedure used to obtain Lighthills equation, we can derive the equation for density variation = 0 2 22 2 00 22'' 'ij ijsi iijijiT f f cfuf txxxxxtx (1.26) The values of and Tij 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 ri gid 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 2Tij / xi xj in the outer region of the surface due to the flow a surface distribution / xi (ij (f) f/ xj ) due to the interaction of the flow with the moving bodies a surface distribution /t ( 0 usi(f) f/ xi ) due to the kinematics of the bodies Like in the Lighthill analogy, the anal ogy 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. Vi(t) Fluid at rest S(t) Surface with rigid-body motion Ve Real flow field s u n PAGE 30 16 The merits of Lighthill and Ffowcs Williams Hawkings analogies The Lighthill analogy marked an important milestone in the study of flowgenerated acoustics. His theory is able to ex plain a number of importa nt characteristics of acoustic radiation from a jet. The soluti on 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 understand ing aeroacoustic propagation. For example, his theory gave a mathematical representation of sound and pseudosound, an approximation of the sound emission from subs onic cold-air jet flow (which he called self-noise, i.e., noise generate d by turbulent-turbulent interac tion, and shear noise, i.e., noise generated by turbulent mean shear), for the first time. Lighthills 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 aircra ft engine fan, a compressor, or a turbine. The solution of the second analogy is again solved with the Greens function. In this, it is difficult to obtain a general solution, because the analyti cal 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 usef ul 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. PAGE 31 17 Nonlinear problem: Shock wave formation An example where a linear approximation is not suitable is subsoni c flow in a pipe, because waves do not attenuate fast with propag ation. 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 10-2 in the wave propagation can have a significant effect (Hirschberg and Rienstra, 1994). The most sp ectacular nonlinear effect is the formation of shock waves as a result of the steepe ning 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 ma ss and momentum conservation: 0dp ucu txc (1.27) which implies that dp Ju c (1.28) is invariant along the characteristic path C in the (x, t) plane which are described by the equation Cdx uc dt (1.29) When a C+ wave propagates into a uniform region, the Cwave emanating from the uniform region will all carry the same information: J= J0 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 (J+ + J0) (1.30) PAGE 32 18 and 01 2dp JJ c (1.31) In other words u and dp/c are constant along C+, the thermodynamic variable. Hence, speed of sound can be consid ered constant, so that (u + c) of C+ is constant. For an ideal gas there is the relation 2 1dpc c (1.32) 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 ts: 1 0 sduc t dx (1.33) which for an ideal gas with constant can be expressed in te rms of the space derivative of pressure at x = 0: 0 01sp t p c x (1.34) Because the variation of the speed of sound is negligible, the distance of appearance of the shock wave is given by 0 0 01ssp xct p x (1.35) For harmonic waves 0 max dp kp dx (1.36) PAGE 33 19 0 0 1sp kx p (1.37) Hence, for an amplitude 0 / p p= 10-2, the shock is expect ed appear within 10 acoustic waves (Goldstein 1976). Computational Techniques for Aeroacoustics Direct Numerical Simulation (DNS) Direct simulation methods aim to co mpute both unsteady flow and the sound generated by it. These methods must use a domain that includes th e noise producing flow region and at least a part of the near-acousticfield. The computational mesh needs to be selected so that both th e flow and its sound can be well repr esented. The first issue in this approach is that the computati onal cost of such direct com putations is large, hence only simple flow configurations can be studied using this direct method. Direct simulation of the acoustic field solves the compressible Navier-Stokes equations (or Euler equations in those cases where viscosit y is not important) without further approximation. These equations gov ern 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 fr om 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 nu mber is high due to the range of scales present in the flow field. The characteristic frequency of a source th at radiates sound is given by the Strouhal number, St = fL/U O(1), where L and U are the characteristic length and velocity. This implies that the wavelength of sound produced is L/M, PAGE 34 20 where M is the Mach number of the flow. On th e other hand, the dissipation of turbulent energy takes place at the Kolmogorov length scale = L/Re3/4 (Tennekes and Lumley 1972), where Re is the Reynolds number of the flow. T hus, the ratio of th e wavelength of sound to the size of the energy dissipating eddies, /Re3/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 ut ilized 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 vi scous and acoustic phenomena appears to be preferable for high Reynolds number flows. Perturbation Technique A perturbation technique consists of splitti ng up the flow field calculation into two parts. First, a time-independent viscous background flow is calculated, and then a perturbation equation (that describes the s ound) about this background flow is derived, and viscous action on the pert urbation is neglected. In this approach, an initial disturbance is introduced upstream which cause s an instability and causes the wave to grow resulting in the radiation of sound. Ta m and Morris (1980) developed this approach analytically for shear layers, and Tam and Bu rton (1984) developed it for subsonic jets. Hardin and Pope (1995) developed a s lightly different a pproach to address vorticity-dominated flow, an expans ion about the time-dependent, viscous, incompressible, subsonic flow. If the density = 0 is constant, then the continuity and momentum equations become a set of four equations for the three incompressible velocity components and incomp ressible pressure. One solves this set of equations by PAGE 35 21 using a grid and numerical scheme designed fo r the viscous aspects of the flow for the time dependent viscous incompressible flow. Us ing 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 (Coloni us and Lele, 2004). The numerically difficulty of this approach stems from the fact the full LEE se t admit non-trivial inst ability wave solution of the homogeneous equation. Propagation of linear acoustic waves thr ough 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 co mputational 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 a pproach 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 near-field fluctuation, and because they mu st propagate with little attenuation over long distances. In practice th is has dictated the use of highorder-accurate numerical methods. PAGE 36 22 The Numerical Approach to Reduce Dissipation and Dispersion In computational aeroacoustics (CAA), accu rate prediction of the generation of sound is demanding due to the requirement of preserving the shape and frequency of wave propagation and generation. Furthermor e, the numerical schemes need to handle multiple scales, including long and short wa ves, 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) th at in order to conduct satisfactory CAA, numerical schemes should induce minimal disp ersion and dissipation errors. In general, higher-order schemes are more suitable for C AA than lower-order schemes since, overall, the former are less dissipative. That is w hy higher-order spatial discretization schemes have gained considerable interest in com putational acoustics (Hixon, 1997; Kim et al., 1997; Lin and Chin, 1995). For longer wavelengths, the formal order of accuracy is suffici ent to indicate the performance of a scheme. However, for shorte r waves relative to the grid size, it is known that the leading truncation error term s 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 as sociated with a particular range of wave numbers has been used over the years by ma ny researchers, e.g., Hu et al. (1996), Stanescu and Habashi (1998), Shu (1994), Na nce 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 Dispersion-Relation-Preserving (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 th e coefficients to satisfactorily resolve short PAGE 37 23 waves with respect to the computational grid, namely, waves with wavelengths of x (defined as 6-8 points per wave or PPW) or shorter. It maximizes the accuracy by matching the wave number and frequency char acteristics 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 stenc il 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 (199 8) 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 th e 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 ph ysics, 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 irre gular geometry and moving bounda ries. In this work, we investigate a finite volume formulation (whi ch we call DRP-fv), extending the concept embodied in the original, finite differencebased DRP scheme (which we call DRPfd). PAGE 38 24 Similarly, for the OPC scheme, we extend the basic concepts of the original, finite difference-based OPC (OPC-fd) scheme, in a finite volume formulation, named OPC-fv. Our ultimate goal is to extend the finite volume version of suitable schemes into a cutcell type of Cartesian-grid computational t echnique 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 intere st, 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, Che ong and Lee (2001), proposed a grid-optimized dispersion-relation-preserving scheme (GODRP). They considered the approximation of the derivative by 1n ji jnu auxxx x x (1.38) where /2nn x xxn (1.39) ()/(,1,...,1,)ii x xxxinnnn (1.40) x0x1x2x-1x-2 xn-1xnx-nx-n+1( =x ) x0x1x2x-1x-2 xn-1xnx-nx-n+1( =x ) Figure 1-3. Schematic diagram showing the (2n + 1) point stencil on a nonuniform grid PAGE 39 25 The wavenumber of the scheme over a nonuniform grid is obtained by using a Fourier transformation: in ixx j jni ae x (1.41) The scheme is derived by requiring that i) th e 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 2 0 2 0expln2r ie real e imagExxdx x x Sgncxdx (1.42) The terms er and ei denote the upper limits of the in tegration intervals of the real and imaginary parts, respectively. The term is the weighting factor, is the half-width 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 th e accuracy of the finite difference method for curvilinear mesh es from the wave number point of view. Through the grid-optimization process, high-o rder finite difference equations can be solved with curvilinear grids with a guarantee of local and thus resultant global dispersion-relation-preserving pr operties. Hence, the approach can be used with success for a body-fitted 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+1 point stencil GODRP spatial disc retization, the scheme implies a PAGE 40 26 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 x 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 f(x) = ei[x + ()] (1.43) The total effect on the function f(x) will be weighted by the amplitude A(). The general formula for extrapolation is given by 1 00 0,N jjj j f xxSfxxxjx (1.44) where Sj (j = 0, 1, 2, N-1) are the stencil coefficients Their values are obtained by imposing the following: i) the difference betwee n 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 constraine d optimization problem is handled by the Lagrange multiplier 2 11 00 01k NN iyijy jj jjLeSedyS (1.45) 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 1-4 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 PAGE 41 27 and G2E is perpendicular to AB. The value of th e 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 instabi lity when the wavenumber is high. Hence, we do not expect that the approach gives accurate results for short waves or high wavenumber. G1G4G3G2G5 A B 1237 6 5 4 1 2 3 7 6 5 4 G1G4G3G2G5 A B 1237 6 5 4 1 2 3 7 6 5 4 G2A B E 1 2 12 G2A B E 1 2 12 a) b) Figure 1-4. Cartesian boundary treatment of curved wall surface; b) detail around boundary Scope The present thesis has three main contribu tions. First we investigate two numerical methods for treating convective transport: the dispersion-relati on-preservation (DRP) scheme, proposed by Tam and coworkers (T am and Webb, 1993; Tam et al., 1993), and the space-time amethod, developed by Chang (1995). The purpose is to examine the characteristics of existing schemes capable of handling acoustic problems. The space- PAGE 42 28 time amethod directly controls the level of dispersion and dissipation via a free parameter, 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,, are different between them. Different performance characteristics are observed betw een 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 develo ped by Ashcroft and Zhang (2003), are extended from their original finite differen ce 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. Lin ear 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 linea r and nonlinear wave equations with different wavelengths as test cases, the performance of these approaches is studied. Finally, the Cartesian grid, cut-cell method is extended along with the OPC-based finite volume scheme so that this high order scheme can treat curved geometry associated with practical acoustic applications. The a pproach 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 us ing flow properties at the cell surface interpolated from cell nodes while maintaining desirable accuracy level. PAGE 43 29 The rest of the thesis is structured as follows. Chapter 2 presents the investigation of the Dispersion-Relation-Preservi ng (DRP) scheme, and the space-time amethod. The characteristics of these two schemes are emphasized us ing a simple wave equation with the initial disturbance in the form of the Gaussian function Chapter 3 presents principal characterist ics and introduces two space discretization schemes: the DRP scheme and the OPC sche me. A low dispersion and low dissipation Runge-Kutta proposed by Hu and coworkers (1996) is employed for the time stepping procedure, and combined with the DRP a nd 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 a nd nonlinear test problems are presented to evaluate the merits of these schemes. Chapter 4 presents the principal characte ristics for the Cartesian grid, cut-cell approach. Additionally, we present the proposed adjustment for the acoustic approach. Finally, several test problems are presented to demonstrate th e performance of the present approach. Finally, in Chapter 5, summary and future perspectives are presented. PAGE 44 30 CHAPTER 2 ASSESMENT OF DISPERSION-RELATI ON-PRESERVING AND SPACE-TIME CE/SE SCHEMES FOR WAVE EQUATIONS Introduction A number of numerical schemes based on va rious concepts have been proposed to treat wave propagation and c onvective transport, including the concepts of upwinding, flux splitting, total variation diminishing, monotonicity, non-oscillatory, higher order differencing, and Riemann solvers. For some representative works, see van Leer (1979), Roe (1981), Harten (1983, 1989), Osher a nd Chakravarthy (1983), Shyy (Shyy, 1983; Shyy et al., 1997), Leonard (1988), Shu a nd Osher (1988), Hirs ch (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 dispersion-relation preservation (DRP) scheme based on an optimi zed high-order finite difference concept, proposed by Tam (Tam and Webb, 1993; Tam et al., 1993), and the space-time 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: 0uu c tx where a > 0 is a constant (2.1) PAGE 45 31 The exact solution of the wave equation for the initial value problem with initial data (,0)() (,)() uxUx uxtUxct (2.2) 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 soluti on without creating wiggles. In other words, the goal is to offer satisfactory resolution, especially for the high wave number components, by balancing numerical dissipati on and dispersion. Here we use the term balancing because of the finite resolution possessed by numerical computation. The two methods optimized higher-order finite difference (DRP) method (Tam and Webb 1993) and the space-time conservation element and solution element (CE/SE) (Chang, 1995) address the aforementioned issu es from different angles. The philosophy of the DRP method is to maximize the accuracy by matching the wave number and frequency characteristics between the analytic al and the numerical ope rators in the range of resolvable scales. The space-time CE/S E method views the flux calculations in a joined space-time conservation element, taki ng into account the unifi ed 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 appro aches 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 PAGE 46 32 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 Webbs paper (1993) and in Changs paper (1995) will be summarized to help discuss the related analysis and numerical asse ssment conducted in the present work The Dispersion-Relation Preservation (DRP) Scheme In the following, we summarize the met hod by Tam and Webb (1993) to highlight the salient features of DRP scheme. Discretization in Space As illustrations, x is approximated by the central difference schemes of various orders considering uniform grids with spacing to x we give three examples. The second order central di fference approximation is 2 11() 2ll luu u Ox xx (2.3) Fourth-order approximation is 4 211288 () 12llll luuuu u Ox xx (2.4) Sixth-order approximation is 6 321123945459 () 60llllll luuuuuu u Ox xx (2.5) The construction of the scheme (Tam and Webb, 1993) 1N j jNu x auxjx xx (2.6) PAGE 47 33 is based on two goals: (i) the behavior of th e numerical solution in the resolvable wave number range closely matches that of the ex act solution, and (ii) the formal order of accuracy of scheme spanning 2 N + 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 1 2ix f fxedx (2.7) ()()ix f xfed (2.8) as well as the derivative and shift theorems ~if xfx (2.9) ~ifxef (2.10) With these tools, the relationship between the differential and discrete operators (given in Eq.(2.6)) can be approximated as 1N ijx j jNifaefif x (2.11) where N ijx j jNi ae x (2.12) We see that is a periodic function of x with a period of 2 Obviously the goal is to ensure that is as close to as possible. To accomplish this goal the error is minimized over a certain wavenumber range, x [ the numerical dispersion is reduced by narrowing the range of optimization (Tam and Webb, 1993) PAGE 48 34 2() Exxdx (2.13) It is noted that is real, and hence the coefficients aj must be anti-symmetric, i.e., a0 = 0 and a-j = 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 2 12sinN j jEkakjdk (2.15) where k = x Hence, with the anti-symme tric condition, one can us e the method of the leastsquares fitting to minimize E over a wave number range, namely: 0,lE a l = 1, 2, N (2.16) Furthermore, to ensure that the scheme is accurate to O ( x2( N -1)), additional conditions can be made by expand ing the right side of Eq.(2. 6) as a Taylor series and then equating terms. In this way we have N -1 independent equations with N unknowns: ( aj)j=1,N. Note that aj = a-j, j = 1,, N and a0 = 0. Tam (Tam and Webb, 1993) s hows the relationship between x and x for the 4th order, DRP scheme and, the aforementioned 6th-order, 4th -order and 2th -order central difference schemes over the interval 0 to For x up to c x the individual curves are nearly the same as the straight line x x Thus, the finite difference scheme can provide reasonable approximation for wave number so x becomes less than c x If we wish to resolve a short wave with PAGE 49 35 straightforward central differen ce approximations using a fixed size mesh, we need to use a scheme with a large stencil. On the other hand, the 7-point DRP scheme has the widest favorable range of x Tam and Webb (1993) demonstrate that DRP can be effective in improving the performance of a given stencil w ithin certain wave number (see Figure 2-1). In the following, we will only c onsider the DRP scheme with a symmetric stencil and in particular the 7-point formula. 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 x DRP 6th order 4th order Figure 2-1. x versus x 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 importa nt to evaluate the group velocity of a finite difference scheme. The gr oup velocity is characterized by /dd which should be almost one in order to reproduce exact resu lt. A way to reduce dispersion is to adjust the range of the wave number in the optimization process. b shows / dd curves of the DRP scheme, for different ranges of the optimization parameter Upon examining the corresponding /dd curves, x = 1.1 gives the best overall fit. In this case based on the criterion 1.0.01 d d the optimized scheme has a bandwidth 15% wider than the standard 6th order central difference scheme. PAGE 50 36 For long wave the important range of x is small, and hence any value of less than 1.1 is reasonable. For short waves, on the other hand, disper sion can be noticeable with any choice of Table 2-1. The stencil coefficient for N = 3 a1 a2 a3 /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 0 0.2 0.4 0.6 0.8 1 0.99 0.992 0.994 0.996 0.998 1 1.002 1.004 1.006 1.008 1.01 x standard-2th order standard-4th order standard-6th order DRP scheme=1.1 f(x)=1 () () dx dx 0 0.2 0.4 0.6 0.8 1 0.99 0.992 0.994 0.996 0.998 1 1.002 1.004 1.006 1.008 1.01 x standard-2th order standard-4th order standard-6th order DRP scheme=1.1 f(x)=1 () () dx dx a) 0 0.2 0.4 0.6 0.8 1 0.99 0.994 0.998 1.002 1.006 1.01 x =0.7 =0.9 =1.1 eta= /2 f(x)=1 () () dx dx 0 0.2 0.4 0.6 0.8 1 0.99 0.994 0.998 1.002 1.006 1.01 x =0.7 =0.9 =1.1 eta= /2 f(x)=1 () () dx dx b) Figure 2-2. Dispersive charac teristics of DRP scheme: a) /dxdx versus x for optimized and standard schemes. .b) /dxdx versus x for the 4th order optimized scheme for different wave number range of optimization: range is from to PAGE 51 37 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 t To advance to the next level, we use the following 4-level finite difference approximation: () 3 (1)() 0 nj nn j ju uutb t (2.17) 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: b0, b1, b2, b3. To determine these constants and create a scheme of O ( t3) accuracy, the terms in Eq.(2.17) are expanded in a Taylor series to match exactly up to order t2. This leaves one free parameter, b0. The relationship between ot her coefficients and b0 are (Tam and Webb, 1993; Tam et al., 1993): 102030531623 3;3; 12312 bbbbbb (2.18) One can utilize the Laplace transform to determine b0. First, the Laplace transform and its inverse of a function f ( t ) are related by 01 ; 2itit f ftedtftfed (2.19) where is a line in the upper half -plane parallel to the real-axis above all poles and singularities. Also, the shifting th eorem for Laplace transformation is ~i f tef (2.20) To apply the Laplace transfor m, we need to generalize the equation to one with a continuous variable, namely, PAGE 52 38 3 0() ()()j jutjt uttuttb t (2.21) On applying the shifting theorem to four-level scheme presented in Eq.(2.21), we find that ~ 3 0itijt j ju ueutbe t (2.22) Hence, 3 01it ijt j jie tbe (2.23) To optimize the time stepping scheme, the error is optimized 2 2 1Re()1Im Ettttdt (2.24) i.e., 1 00. dE db In Eq.(2.24) is the weight and is the non-dimensionalized frequency range needed for the numerical scheme to match th e exact solution. Substituting Eq.(2.18) and (2.23) into Eq.(2.24) E1 becomes a function of b0 alone, For and the scheme becomes (Tam and Webb, 1993) 01232.30256;2.49100;1.57434;0.38589 bbbb (2.25) Eq.(2.23) indicates that the relations hip between t and t is not one to one. This means that spurious solutions appear. Th e stability will be established in function of the real solution. It can be written in the form: PAGE 53 39 432 32100 ii bzbzbzbz tt (2.26) where z = e i t and the values of bj ( j = 0, 1, 2, 3) are given by Eq.(2.26). To maintain satisfactory temporal resolution while being stable, the imaginary part of the solution ( t ) should be negative but close to zero. The interval 00.41t satisfies these expectations. Furthermore, Re()Re() tt in this range. To ensure that the damping numerical is minimized, Tam and Webb (1993) suggest the condition 4Im()0.11810 t (i.e., 00.19t is adopted). This condition guarantees numerical stability and negligible damping. To compute stability of the scheme we take into consideration Fourier-Laplace transformation of the wave equation (Eq.(2.1)) 1 2initialiuciuu (2.27) where characterize the PDE For the long wave we can approximate wavenumber of the scheme with wavenumber of the PDE (2.28) which leads to initialucuku (2.29) Hence *ckk (2.30) The condition of the numerical stability is that amplification factor for time discretization is less than 1, and hence *t <0.41. It is also noted from Figure 2-1 that 1.8x (2.31) PAGE 54 40 hold true. By introducing Eq. (2.31) into (2.30) and upon multiplying by t it is found that 1.8M1c tkkt x (2.32) where M is mach number. Therefore to ensure nu merical stability it is sufficient by Eq. (2.32) to restrict t to be less than tmax, where tmax is given by max0.41 1.8M1 x t kkc (2.33) Therefore, for t < tmax the schemes are numerically stable. Consequently, the schemes yield the following criteria for numerical stability: 0.21CFL (2.34) Based on Eqs.(2.1), (2.6), and (2.17) one can obtain the final form of the DRP scheme with 7-point in space and 4-point in time: 33 1 03nnnj lljklk jkt uucbau x (2.35) The leading truncation error of the scheme (given in Eq.(2.35)) can be evaluated using the Taylor series expansion, yielding: 3 3 33 0 454 33 435 031 64 ()3 2 4!105!txjxxxx j xxxxxk jk jkucucjbux uxOxc cjjbak (2.36) Replacing in Eq.(2.36) by the numerical va lues of the various coefficients, the scheme becomes: 334451.63181.2230.012814()tx xxxxxxxxxuau auxaauxOx (2.37) PAGE 55 41 The following can be summarized in regard to the present DRP scheme: the scheme is forth order in space, and third order in time < 0.21 assures stability and reasonable accuracy the first term on the RHS of Eq.(2.37) is dissipative. The accuracy bound, < 0.21, indicates that the coefficient of the l eading dissipative term is small. This observation suggests that, dependi ng on the relative magnitude of and x dispersive patterns may be dominant in numerical solutions. The Space-Time Conservation Element and Solution Element Method The tenet in this method is to treat local and global fl ux conservation in a unified space and time domain. To meet this requ irement, Chang (1995) introduces solution elements, which are subdomains in the space -time coordinates. Within each solution element, any flux vector is then approximate d in terms of some simple smooth functions. In the last step, the computa tional domain is divided into conservation elements within which flux conservation is enforced. Note th at a solution element generally is not the same as a conservation element. We su mmarize in the following the key concepts adopted in Changs approach. a Scheme Consider Eq.(2.1), and define F1 = -u F2 = au x1 = x , and x2 = t Applying Green theo rem we obtain: 21 1122 12SFF dxdyFdxFdx xx (2.38) where 1122 SFdxFdx is contour integral on closed curve S. Following a vector notation, Eq .(2.38) can be written as 1122 SFdxFdx = SVgds (2.39) PAGE 56 42 where 12(,) gFF and 12, dsdxdx (dx1, dx2) is a differential vector associated with a point (x1, x2) on the closed curve S. Because Eq.(2.1) is valid anywhere in the definition domain, we obtain 0SVgds (2.40) where (i) S(V) is the boundary of an arbitrary space-time region V in E2 (a twodimensional Euclidean space) with x1 = x and x2 = t; (ii) SVgds is contour integral on closed curve S(V); (iii) (F2, -F1) is a current density in E2. Note that gds is the spacetime flux of (F2, -F1) leaving the region V through the surface element. As shown in Figure 2-3a, is a set of mesh points (j, n) that is adapted to discretize a physical domain, where 2,ijn with i = 0, 1, 2, 3, To facilitate the construction of the present scheme a solution element (SE) associated with (j, n) is illustrated in b. For any (x, y) SE(j, n), u(x, t), is approximated by u(x, t; j, n): *(,;,)()()()()nnnn jxjjtjuxtjnuuxxutt (2.41) with n ju, n x ju, and n t ju are constants in SE(j,n), ,n j x t are coordinates of the mesh index ( j,n ), and ** *(,;,)(,;,) (,;,),,nn n jjxt jjuxtjnuxtjn uxtjnuuu xt (2.42) PAGE 57 43 Figure 2-3. Scheme of th e solution elements ( SE s) and conservation elements ( CEs ): a) The index positions of SE s and CE s; b) SE (j,n); c) CE-(j,n) ; d) CE+(j,n) ; e) CE+(j-1/2, n+1/2) ; f) CE_(j+1/2, +1/2) ; n+1 n+1/2 n n-1/2 n-1 j-3/2 j-1 j1/2 j j+ 1/2 j+1 j+3/2 a ) /2 x /2 x /2 t /2t (j,n) b ) c) d ) (-12,-12) SEjn (j,n) A B (j-1/2, n-1/2) (,) SEjn (12,-12) SEjn (j,n) A (j+1/2, n-1/2) (j,n) (j,n) f) e ) (,) SEjn (12,12) SEjn (-12,12) SEjn PAGE 58 44 Eq.(2.41) has the form of the first-order Ta ylor series expansion. Furthermore, if one takes into account that nn tx jjucu, because u(x, t; j, n) satisfies wave equation, then Eq.(2.41) becomes *(,;,)()[()()]nnn jxjjuxtjnuuxxctt (2.43) In each SE(j, n), (,)((,),(,)) gxtuxtcuxt is approximated by ***(,;,,)((,;,),(,;,))gxtjnuxtjncuxtjn (2.44) In order to develop appropriate approxima tion for the flux, one divides the physical domain into nonoverlapping rectangular elemen ts, referred to as conservation elements (CEs). Specifically, CE receive sign - or + in func tion of the slope of the line that connects the vertex from of a CE. If the slope is negative, CE receives the positive sign; otherwise the sign is negative In Figure 2-3c the line that unifies the vertices founded in is positive, and CE index is (j,n). The final notation is CE-(j, n). The surface of CE belongs to two different SE, SE(j, n) and SE(j-1/2,n-1/2). To specify that a part of surfa ce is in a certain SE in Figure 2-3c, that part is around by a certain type of line. Figure 2-3d, Figure 2-3e, and Figure 2-3f illustrate three other corresponding cases relating the conservation element to the solution element. The approximation of Eq.(2.40) is ,,0SCEjnFjngds (2.45) for all (j, n). ,SCEjngds is contour integral over closed path S(CE(j,n)) and represents the total flux l eaving the boundary of any conser vation element is zero. The PAGE 59 45 flux at any interface separating two neighboring CEs is calculated using the information from a single SE. As an example, the interface AC de picted in Figure 2-3c) and Figure 23d) 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: 12 2212 1 12 2 2 1221 4 ,11nn nn xxjj jjFjnuuuu x x (2.46) where, again, ct x is CFL number. With the aid of Eqs.(2.45) and (2.46), n ju and n x ju can be solved in terms of 12 12n ju and 12 12n x ju if (1 2) 0, i.e., ,12,1212,12 qjnQqjnQqjn (2.47) where (,) /4n j n x j u qjn xu for all (j, n) (2.48) 2 211 11 11 ; 22 11 11QQ (2.49) We have employed the Taylor series expansion to eval uate the lead ing truncation errors of Eq.(2.46), yielding th e following two expressions for F+ and Frespectively: PAGE 60 46 2 24 2 24() 48161648 () 48161648txxxx txxxxcccc ucuuxOx cccc ucuuxOx (2.50) From our derivation indicat ed by Eqs.(2.49), and (2.50) we can establish the following observation for the space-time CE/SE scheme: The method is second order in space and time; The method is stable if 2 < 1 The dispersion strength of the scheme increases when (or time step) is reduced because appears in the denominator in th e leading truncation error term in Eq.(2.50). Hence a small value of may not be desirable. The third order term O( x3) is zero in Eq(2.50), i ndicating that not only the 2u/x2 term but also the 4u/x4 terms are non-present. Th is is the second reason, in addition to previous observation, why the present scheme can be highly dispersive. aScheme In order to help to reduce the dispersi ve 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: 2 21 (,) 4n x jx Fjndu (2.51) where is an independent parameter of the numerical variables, and 1/21/2 1/21/2 1 1/21/2 2 1/21/2/nnn nn xxxjj jjjduuuuux (2.52) Because the magnitude of the added terms in the scheme is controlled by the numerical dissipation is controlled by However, because F( j n ) 0 if 0, strictly specking CE+( j,n ) and CE-( j,n ) are no longer conservation elements in the ascheme. On the other hand, although the net flux entering the interface separating CE+( j,n ) and PAGE 61 47 CE-( j,n ) is not zero, but the sum of F+( j n ) and F-( j n ) is zero. Hence the total flux leaving CE(j,n) vanishes. As a result, CE(j,n ) is a conservation element in the ascheme. 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 ascheme as: ,12,1212,12 qjnMqjnMqjn (2.53) where 2 2 11 2211 11 ; 121 121 MM (2.54) To assess the leading truncation errors in the ascheme a modified equation is derived by substituting Taylor series expansion into Eq.(46) for 1 2 1 2n x ju 1 2 1 2n ju The resulting equations are: 2 2 3 2 2 3124412332 662 124412332 662txxxx xxxx txxxx xxxxccccccx ucuu ccx uHOT ccccccx ucuu ccx uHOT (2.55) From Eq.(2.55), we can see that the ascheme introduce dissipative term uxxxx and reduce the strength of the dispersive uxxx. 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 which is different from the DRP scheme discussed previously. As PAGE 62 48 already mentioned, in the ascheme, both dispersive and dissipative aspects are affected by In Refs. (Chang, 1995; Chang et al., 1999), insi ght is offered in regard to the choice of and based on the von Neumann stability analysis of the ascheme, without detailed information of the tr uncation error. The analysis reported in these references suggests that: (i) should be large; (ii) shouldnt be close to 0 or 1; (iii) the spurious solution can be effectively suppressed for = 0.5 in the case of long-wavelength. Reference (Loh et al., 1998) suggests that for ae roacoustic computation, it is essential to choose a large CFL number and a small Further insight can be gained based on the present truncation error analysis. If we relate and ( CFL number) by = + 1, then, in Eqs. (2.55): Coefficient of uxxx in first equation of the system Eq. (2.55) is equal to: 2 11(1) 34c (2.56) Coefficient of uxxx in second equation of system Eq. (2.55) is equal to: 2 161 1 34c (2.57) The above equations indicate that if is close to 1, and and are close to each other, i.e., if 1 is small, then the leading dispersion error shown in Eq.(2.55) is small. On the other hand, to maintain adequate numer ical dissipation, it is helpful to let vary in the same manner as These observations as well as the analysis in Chang (1995) indicate that the value of and should be coordinated. Further evaluation will be made based on numerical computations, to help establish a more explicit guideline. PAGE 63 49 Numerical Assessment of the DRP and Space-Time Schemes To assess the individual and relative mer its of the DRP and space-time schemes the following simple test problem is adopted. 0. uu tx (2.58) The initial condition imposed at time t = 0 are: 2 2expln2 2ln2expln2xx u b xx u bb (2.59) which is a Gaussian profile. The exact solution is: 2 2expln2 2ln2expln2xxt U b x txt U bb (2.60) In this study we evaluate the performance of the schemes in short, intermediate, and long waves relative to the grid sp acing, which is assured by the value b / x. In all cases x =1, and c = 1, so that t = For the optimized DRP scheme we adopt the one with 4th level in time and 7-point in space, and 4th order formal accuracy in space and 3rd order formal accuracy in time: 3 1 0n nnj lljl juubF (2.61) where t x (because c = 1) PAGE 64 50 3 3kk ljlj jFau (2.62) The values of the of bj are given in Eq.(24), and the value aj are as chosen to minimize difference between d d and 1 (see Table 2-1). For the space-time ascheme, the following formulas result: 11 2 22 11 22 11 2 22 11 22 11 22 11 22 11 22 11 221 1(1)() 24 1(1)() 4 1 ()1(21)() 24 1(21)() 4nn n jx jj nn x jj nn n xjx jj nn x jjx uuu x uu x uuu x uu (2.63) Before evaluating the DRP and space-time amethod, we furt her address the influence of the parameters: and for the space-time ascheme, and for the DRP scheme. An effort is made in this study to constr uct 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 / x = 3), and are preferable to be close to each other (see Figure 2-4 and Figure 2-5) For long and intermediate waves there are virtually no need to introduce mu ch numerical dissipation, hence can be reduced close to zero. On the other hand, mismatched and can significantly reduce the accuracy of the scheme. To demonstrate this fact, Figure 2-4a and Figure 2-4b show two solutions for PAGE 65 51 the intermediate wave, all with = 0.5, respectively for = 0.9. Changing from 0.5 to 0.99 causes significant numerical di ssipation. On the other hand, = 0 can be a very acceptable choice in for long wave. Figure 25 compiles long and intermediate wave solution with different value of and Taking into consideration the previous observations, it is decided that = = 0.99 is a good choice, and is used in the results presented for the space-time 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 coor dinated choice is beneficial. For the present DRP scheme a value of less than 0.21 guarantees numerical stability and negligible numerical da mping. If we decrease the value of much further, then, as indicated in Eq.(2.37), the numerical damping, as indicated by the leading dissipation term uxxxx, may be less than adequate due to the small value of For these reasons we decided to use = 0.1 for longer waves and = 0.2 for short wave ( b / x = 3). To summarize, the DRP scheme is highe r formal order than the space-time ascheme. The DRP scheme prefers smaller (or t ), indicating that it is more expensive to compute than the space-time ascheme. To evaluate the solution accuracy, we define an error vector as: 1,...,T NEEE (2.64) U(xi) is the exact solution at the point xi, and ui is the numerical solution at the point xi. The error norm is adopted as norm one PAGE 66 52 1 N i iE E N (2.65) where (),1iiiEUxuiN (2.66) which will be used to help to measure the order of accuracy in actual computations. x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.10E-01 analyticalsolution x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.10E+00 analyticalsolution x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.50E+00 analyticalsolution x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.70E+00 analyticalsolution x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.90E+00 analyticalsolution x u(x,t) 160 180 200 220 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution=0.99E+00 analyticalsolution a) x u(x,t) 20 30 40 50 60 0 0.2 0.4 0.6 0.8 1numericalsolution=0.50E+00 analyticalsolution x u(x,t) 20 30 40 50 60 0 0.2 0.4 0.6 0.8 1numericalsolution=0.99E+00 analyticalsolution b) Figure 2-4. Comparison between analyti cal and numerical solutions Effect of on the accuracy of space-time ascheme: a) = 0.9, b/ x = 3, t = 200; b) b/ x = 6, t = 200, = 0.5 PAGE 67 53 0.25 0.5 0.75 1 1.80E-03 2.20E-03 2.60E-03 error 0.25 0.5 0.75 1 2.5E-02 3.5E-02 4.5E-02 error 0.25 0.5 0.75 1 4.0E-02 6.0E-02 error a) 0.25 0.5 0.75 1 2.0E-02 4.0E-02 6.0E-02 8.0E-02 error 0.25 0.5 0.75 1 1.5E-02 2.5E-02 3.5E-02 4.5E-02 5.5E-02 error 0.25 0.5 0.75 1 4.65E-04 4.75E-04 4.85E-04 4.95E-04 error b) 0.25 0.5 0.75 1 2.0E-02 4.0E-02 6.0E-02 8.0E-02 1.0E-01 1.2E-01 error 0.25 0.5 0.75 1 5.0E-03 1.0E-02 1.5E-02 2.0E-02 2.5E-02 error 0.25 0.5 0.75 1 4.06E-05 4.07E-05 4.08E-05 4.09E-05 error e error 0. 2 5 0. 5 0.7 5 1.53E-03 1.55E-03 1.57E-03 e error 0. 2 0.4 2.04E-03 2.06E-03 2.08E-03 2.10E-03 2.12E-03 c) Figure 2-5. The depend ence of the error on for the space-time ascheme at t = 200: a) b/ x = 3; b) b/ x = 6; c) b/ x = 20 PAGE 68 54 In the following we will present the results based on three categories: long wave (b/ x = 20) intermediate wave (b/ x = 6) short wave (b/ x = 3) They are defined according to the ratio b/ x, where b is a parameter that characterizes the wavelength. Figure 2-6 offer an overview of the two schemes for long, intermediate and short wave computation at two different time instants: t/ t = 200 and 2000. It is noted that since x is fixed in each plot, varying is the same as varying t. 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 : (iii) a slower increase with respect to 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 can cause insufficient numerical damping. Hence, it is advisable to use a somewhat larger value of but still within the stability bound of The spacetime ascheme has exhibited unusual e rror variations as a function of (or t): the error changes very little with and then suddenly decreases rapidly for close to 1 (see Figure 2-6 and Figure 2-7). This trend holds for all resolution. This type of behavior is unusual at first glance, bu t can be explained from the viewpoint of truncation error analysis. Specifically, with = both uxxx and uxxxx terms are small as approach 1, making the scheme sharply improves the apparent order of accuracy. Next, we present the error as a function of spatial resolution (b/ x). PAGE 69 55 0.1 0.2 1.3E-03 1.3E-03 1.4E-03 1.4E-03 1.5E-03 1.5E-03 1.6E-03 1.6E-03DRPschemet=200 b/x=3.error1 2 0.01 1 0.1 0.2 1.0E-04 2.0E-04 3.0E-04 4.0E-04 5.0E-04DRPschemet=200 b/x=20.error 1 0.90.4 1 0.1 0.2 1.0E-04 2.0E-04 3.0E-04 4.0E-04 5.0E-040.9 1DRPschemet=200 b/x=6.error0.3 1 a1) 0.1 0.2 7.5E-03 7.6E-03 7.6E-03 7.6E-03 7.7E-03DRPschemet=2000 b/x=3.error 0.1 0.2 6.0E-04 7.0E-04 8.0E-04 9.0E-04 1.0E-030.4 1DRPschemet=2000 b/x=6.error0.1 1 0.1 0.2 5.0E-04 6.0E-04 7.0E-04 8.0E-04 9.0E-04DRPschemet=2000 b/x=20.error 1 0.40.1 1 a2) error 10-2 10-1 100 5.0E-05 5.0E-04 9.5E-04 1.4E-03 1.9E-03Space-timeCESEschemeb/x=20. t=200 error 10-2 10-1 100 5.0E-04 5.0E-03 9.5E-03 1.4E-02 1.9E-02Space-timeCESEschemeb/x=6. t=200 error 10-2 10-1 100 2.5E-03 1.3E-02 2.3E-02 3.3E-02b/x=3. t=200 Space-timeCESEscheme b1) 0.1 0.2 7.5E-03 7.6E-03 7.6E-03 7.6E-03 7.7E-03DRPschemet=2000 b/x=3.error 0.1 0.2 6.0E-04 7.0E-04 8.0E-04 9.0E-04 1.0E-030.4 1DRPschemet=2000 b/x=6.error0.1 1 0.1 0.2 5.0E-04 6.0E-04 7.0E-04 8.0E-04 9.0E-04DRPschemet=2000 b/x=20.error 1 0.40.1 1 b2) Figure 2-6. The dependence of the error as function of for short (b/ x = 3), intermediate (b/ x = 6) and long (b/ x = 20) waves for: a) DRP scheme; b) Space time ascheme, at two different time instants t = 200 and t = 2000. In all cases t = PAGE 70 56 x u(x,t) 180 190 200 210 220 0 0.2 0.4 0.6 0.8 1numericalsolution=0.99E+00 analyticalsolution x u(x,t) 180 190 200 210 220 0 0.2 0.4 0.6 0.8 1numericalsolution=0.90E+00 analyticalsolution Figure 2-7. Effect of on the accuracy of space time ascheme: b/ x = 3, t = 200 Figure 2-8 depicts the dependence of error for b/ x between 3 and 300. The graphs show that the space-time amethod is second order accuracy in space. In contrast with the DRP scheme does not exhibit consistent tr ends 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 independe nt of space (the slope is close to 0.1). Comparing the performance of these two sche mes, we can deduce that the DRP scheme gives a better solution for b/x less than 10 (for a short and intermediate wave). For long wave, the space-time CE/SE scheme gives a be tter order of accuracy ; nevertheless both schemes are satisfactory in practical terms. Finally, we compare these two methods with respect to their performance with suffi cient accumulation of the time steps. PAGE 71 57 b/x error 101 102 10-5 10-4 10-3 10-2 1Space-Timeschemea 2 t=200 b/x error 101 102 10-5 10-4 10-3 10-2DRPscheme1 4 0.01 1 t=200 b/x error 101 102 10-6 10-5 10-4 10-3 2 1Space-Timeschemea t=200 b/x error 101 102 10-6 10-5 10-4 10-3DRPscheme4 1 0.11 t=200 Figure 2-8. The behavior of the error in f unction of the wavelength: comparison between DRP and space-time a-schemes Short Wave: b /x = 3 Figure 2-9 summarizes various aspects of error accumulation in time, along with selected solution profiles for both schemes. It is apparent that the space time ascheme introduces both dissipation and dispersion errors, while the DRP scheme shows mainly dispersion errors with less dissipation. In addition, the two schemes exhibit diffe rent levels and growth rates. For the space-time ascheme, at the beginning, the rate of accumulation is slower. Then, around t (103), the dissipation becomes more substantial. Around t (5103), both dissipation PAGE 72 58 and dispersion characterize the solution, but the overall error gr owth 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 (2104), dispersion becomes highly visible while dissipation is also present at a modest level. Intermediate Wave: b /x = 6 As shown in Figure 2-10, for this case, the DRP scheme has lower error than the space-time ascheme. The space-time ascheme exhibits similar growth trends in both cases, intermediate and short wave. Again, both dispersion and dissipation errors are noticeable. The DRP scheme behaves differently betw een short and intermediate waves. The error growth rate at the begi nning is low, and then becomes higher. Throughout the entire computation, only minor dispersion error appears. Long Wave: b /x = 20 Figure 2-11 presents the solutions fo r long wave computations. Both schemes perform well over the interval of time. But the space-time ascheme exhibits noticeably slower rate of error accumulation. Neither dispersion nor dissipation errors are significant. PAGE 73 59 x u(x,t) 19980 20000 20020 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.20E+05 analyticalsolutiont=0.20E+05 D x u(x,t) 880 890 900 910 920 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.90E+03 analyticalsolutiont=0.90E+03 B x u(x,t) 5260 5280 5300 5320 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.53E+04 analyticalsolutiont=0.53E+04 C x u(x,t) 80 90 100 110 120 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.10E+03 analyticalsolutiont=0.10E+03 A t 103 104 1.0E-02 1.5E-02 2.0E-02 2.5E-02 3.0E-02b/x=3Space-TimeschemeA B C D 1 1erroraa) x u(x,t) 1475 1500 1525 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.150E+0 4 analyticalsolutiont=0.150E+04 B x u(x,t) 80 100 120 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.100E+03 analyticalsolutiont=0.100E+03 A x u(x,t) 19950 20000 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.200E+05 analyticalsolutiont=0.200E+05 C t 102 103 104 2.0E-03 1.0E-02 1.8E-02 2.6E-021 0.8 1 0.65 A B CDRPschemeb/x=3. error b) Figure 2-9. Accumulation of the error in time for short wave b/ x = 3; a) space-time ascheme; b)DRP scheme PAGE 74 60 x u(x,t) 20000 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.20E+05 analyticalsolutiont=0.20E+05 C x u(x,t) 80 100 120 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.10E+03 analyticalsolutiont=0.10E+03 A x u(x,t) 3575 3600 3625 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.36E+04 analyticalsolutiont=0.36E+04 B t 103 104 5.0E-03 1.0E-02 1.5E-02 2.0E-02 2.5E-02b/x=6Space-TimeschemeA B C 1 1erroraa) x u(x,t) 80 90 100 110 120 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.100E+03 analyticalsolutiont=0.100E+03 A x u(x,t) 19980 20000 20020 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.200E+05 analyticalsolutiont=0.200E+05 B t 102 103 104 1.0E-03 3.0E-03 5.0E-03b/x=6. DRPschemeA B 1 1 error b) Figure 2-10. Accumulation of the error in time for intermediate wave b/ x = 6; a) space-time a-scheme; b) DRP scheme PAGE 75 61 x u(x,t) 50 100 150 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.10E+0 analyticalsolutiont=0.10E+03 A x u(x,t) 19950 20000 20050 0 0.2 0.4 0.6 0.8 1numericalsolution-t=0.20E+05 analyticalsolutiont=0.20E+05 B t 103 104 5.0E-03 5.2E-03 5.4E-03 5.6E-03 5.8E-03b/x=20 Space-TimeschemeA Berror 1 0.12aa) x u(x,t) 50 75 100 125 150 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.100E+03 analyticalsolutiont=0.100E+03 A x u(x,t) 19950 19980 20010 20040 -0.2 0 0.2 0.4 0.6 0.8 1numericalsolutiont=0.200E+05 analyticalsolutiont=0.200E+05 B t 102 103 104 5.0E-04 1.5E-03 2.5E-03 3.5E-03 4.5E-03b/x=20. DRPschemeA B 1 0.9error b) Figure 2-11. Accumulation of the error in time for long wave b/x = 20, a) space-time ascheme; b) DRP scheme PAGE 76 62 Summary and Conclusions It should be noted that he DRP scheme is a multi-step method, which requires more boundary conditions and initial data, while the space-time ascheme is a one-step method. Combined with the fact that the DRP scheme performs better with smaller it can be more expensive to compute than for space-time ascheme. 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 optimi zation procedures, is employed. The present study is restricted to the i nvestigation of a simple 1-D linear wave equation. Obviously, more issues will arise when multi-dimensional geometry, nonlinearity of the physics, and coupling of the dependent variables need to be considered. Nonuniform grid also creates va rying CFL number even with a constant convection speed. Nevertheless, it seems useful to consider, in a well-defined framework, the merits of the two schemes in a well-define d context. In this sense, the present study can be viewed to establish an upper bound of the performance characteristics of the DRP and space-time CE/SE methods. PAGE 77 63 CHAPTER 3 FINITE VOLUME TREATMENT OF DI SPERSION-RELATION-PRESERVING AND OPTIMIZED PREFACTORE D COMPACT SCHEMES FOR WAVE PROPAGAION In computational aero-acoustics (CAA) accurate prediction of generation of sound is demanding due the requirement of pres erving the shape and frequency of wave propagation and generation. Furthermore, th e 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 recogn ized (Hardin and Hussa ini, 1992; Tam and Webb, 1993; Tam et al., 1993) that in orde r to conduct satisfactory CAA, numerical schemes should induce minimal dispersion and dissipation errors. In general, higherorder schemes would be more suitable for CAA than the lower-order schemes since, overall, the former are less dissipative. That is why higher-order spatial discretization schemes have gained considerable interest in computational acoustics (Hixon, 1997; Kim et al., 1997; Lin an Chin, 1995). Table 3-1 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 shorte r 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 w ith a particular range of wave numbers has been used over the years by many researchers, e.g., Refs.(Hu et al., 1996; Stanescu and PAGE 78 64 Habashi, 1998; Nance et al., 1997; Wang an d Sankar, 1999; Cheong and Lee, 2001; Wang and Chen, 2001; Ashcroft and Zhang, 2003) A successful approach is the Dispersion-Relation-Preserving (DRP) fin ite 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 6-8 x (defined as 6-8 points per wave or PPW) or shorter. It maximizes th e accuracy by matching the wave number and frequency characteristics between the analytical and the numer ical operators in the range of resolvable scales. Recently, Ashcroft an d 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 th e biased operators such that their dispersion characteristics match those of the original central comp act 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 back ward operators are equal to those at the corresponding wavenumber of the origin al 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 appr oach (Udaykumar et al., 1999; Yang et al., 1999; Udaykumar et al., 1999), which ensures that fluxes estimated from different sides PAGE 79 65 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 shoc k and turbulence aspects of the aeroacoustic computations. Furthermore, a finite volume fo rmulation 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 difference-based DRP sc heme (which we call DRPfd) to a finite volume fo rmulation (which we call DRPfv). Similarly, for the OPCscheme, we extend the basic concepts of th e original, finite difference-based OPC (OPCfd) scheme, to a finite volume formulation, called OPC-fv. Our overall goal is to develop the finite volume version of DRP and OPC schemes into a cut-cell type of Cartesian-grid computational technique that we have de veloped earlier for moving and complex boundary computations (Yang et al., 1999; Udaykumar et al., 1999; Ye et al., 1999) to treat aero-acoustic problems need ed for engineering practices. We present the finite volume formula tion of both DRP and OPC schemes, and assess both fd and fv versions of the DR P and OPC schemes, using well defined test problems to facilitate systematic evaluation s. Both linear and non linear wave equations with different wavelengths and viscous effects are utilized for direct comparisons. In the following, we first summarize the essenc e of the individual schemes, including derivations, then present assessment of the test cases. Numerical Schemes We use the following one-dimensional wave equation to facilitate the development and presentation of the concept and numerical procedures: 0uu c tx (3.1) PAGE 80 66 The equation contains time and space deriva tive. In our work the space derivative term is treated with DRP or OPC scheme and the time derivative by a Low Dissipation and Low-Dispersion Runge-Kutta (L DDRK) 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 volume-based DRP scheme (DRP-fv) To incorporate the DRP-fd concept into a finite volume framework, let us consider a one-dimensional linear wave equation: 0c tx (3.2) To derive the discredized equation, we employ the grid point cluster shown in Figure 3-1. 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 cont rol volume. For the one-dimensional problem under consideration, we assume a unit thickness in the y and z directions; thus, we obtain 0w ew edxcAA t (3.3) PAGE 81 67 where (A)e and (A)w are the flux across the east and west face, respectively. Figure 3-1.Grid points cluste r for one-dimensional problem Hence, the discretized wave equation (3.2) can be written as 0ewxcAA t (3.4) where is the averaged value of over a control volume. Taking into account Eq. (3.4) we describe the general form of the approximation of x in 1-D using the control volume concept: ` 1ewAA xx (3.5) The general form of the DRP scheme is: 3 31k klk k la xx (3.6) where x is the space grid, and coefficients aj 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 e in the neighborhood: 12213415263 iiiiii e (3.7) 13223145162 iiiiii w (3.8) i e i+1 ( W ) i+2 i+3 w i-1 ( E ) i-2 x (x)w ( x)e PAGE 82 68 Taking into consideration Eqs. (3.5), (3.6), (3.7) and, (3.8) we obtain the values of the i, i = 1,,6 by imposing that the value of at the same locations has the same values as that of the DRP-fd. 3 3k klkew ka (3.9) Hence, the values of coefficients s are 163 2523 34123a aa aaa (3.10) To illustrate the above-described concept, we consider the following equation: 120uuu cc txy (3.11) If we integrate Eq. (3.11) on the surface we have (see Figure 3-2): ,1ij abcdu F tS (3.12) where the resulting DRP-fv scheme is ,1,,,, 2,,,,() ()senw ijijsijeijnijw senw ijsijeijnijwFcuyuyuyuy cuxuxuxux (3.13) ,12,21,3,41,52,63,e ijijijijijijijuuuuuuu (3.14) ,1322,31,4,51,62,w ijiijijijijijuuuuuuu (3.15) ,1,22,13,4,15,26,3n ijijijijijijijuuuuuuu (3.16) ,1,32,23,14,5,16,2s ijijijijijijijuuuuuuu (3.17) PAGE 83 69 121 2abcdsseennww VS sseennwwuu dvcudycudxScuyuyuyuy tt cuxuxuxux (3.18) d c e s w n P b a S W E N d c e s w n P b a S W E N Figure 3-2. Grid notation for two-dimensiona l 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 an d 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 re quires 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 eq ual 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. PAGE 84 70 OPC Scheme Finite-difference-based optimized prefactored compact (OPC-fd) scheme To derive the factorized compact sche me Ashcroft and Zhang (2003) define forward and backward operators F iD and B iD, such that 1 2BF ii iu DD x (3.19) The generic stencil for the forward and backward derivative operators are then defined as: 121121FFFFFFFFF iiiiiiiDDaubucudueu x (3.20) 121121BBBBBBBBB iiiiiiiDDaubucudueu x (3.21) 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 st encil are equal in ma gnitude 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 erro r (weighted deviation) as: 0 r E xxWxdx (3.22) where W( x) 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 th e weighting function. The coefficients are obtained by imposing that, within a given asympt otic order, the error is minimal. In space PAGE 85 71 discretization, one sacrifices formal or der of accuracy in favor of wide-band performance, especially for the short wave components. Finite volume-based OPC scheme (OPC-fv) Taking into account Eq.(3.3) that describes the approxima tion 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 th e function for points at the center of the cell face, namely, e and w assumes similar forms : 0.5() 0.5()FeBe e FwBw wuuu uuu (3.23) where uFe, uBe, uFw and uBw are determined from: 11 FeFe iiiiuubudu (3.24) 11 FwFw iiiiuubudu (3.25) 11 BeBe iiiiuubudu (3.26) 11 BwBw iiiiuubudu (3.27) where the coefficients are the same as those in the OPC-fd scheme: = F= B, = F= B, b = bF= -dB, d= dF = -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 OPC-fd sc heme with the follow boundary stencil: 4 1 1311 ,N BB jjNjj jjNDsuDeu xx (3.28) and PAGE 86 72 4 111 1311 ,N FF NjjNNjj jjNDeuDsu xx (3.29) where the coefficients sj and ej are determined by matching the Taylor series of the forward and backward compact interior stencils to third-order accuracy. The boundary treatment in case of OPC-fv approach is similar to that of OPCfd, but the boundary stencil is computed on the face: Di = (uA)i e (uA)i w (3.30) 33 1 11 33 111 11 ww iiNiNi ii ee iiNiNi iiuauuru uauuru (3.31) where the value of the coefficients are: 1 11 2122 1 3123 32 1 F B N BF NN B F NNNae as assaee asss aeee (3.32) 1 11 2212 1 3123 32 1 B F N BF NN F B NNNre rs reerss rsss reee (3.33) Time Discretization Th e Low Dispersion and Dissipa tion Runge-Kutta (LDDRK) Method Hu et al.(1996) consider time integrati on using the Runge-Kutta algorithm of the differential equation u Fu t (3.34) where the operator F is a function of u. An explicit p-stage algorithm advances the solution of Equation (3.34) from the nth to the (n + 1)th iteration as PAGE 87 73 0 0 (1) 1 () () 1... 1,..., ...n i i i ni i p nuu KtFu KtFu uubKip uu (3.35) where bp = 1, u( p ), 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+1 can be written on short like 1 1 1jjn p p nnj l j j lpju uubt t (3.36) 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 1 11n p j j n ju rit u (3.37) The exact amplification factor is *iti eree (3.38) The numerical amplification factor r in (3.37) is viewed as an approximation of the exact factor. The order of the optimized Ru nge-Kutta scheme is indicated by the leading PAGE 88 74 coefficient in (3.37) that matche s the Taylor series expansion of e-i. For instance, the third order algorithm is obtained by setting j = 1/j! for j = 1, 2 and 3. To compare the numerical and exact solutions we take into consideration the ratio: i er re r (3.39) where |r| represents the dissipation rate (obv iously, the correct value should be 1), and represents the phase error (or dispersive erro r) where the correct value should be zero. Hu et al (1996) obtained coefficients of th e low dispersion and dissipation RungeKutta (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 dissipati on errors are minimized. In other words the following integral is minimized: 2 0 11minp j i j jcied (3.40) and iii) the amplification factor of the scheme is less than 0ne within the given stability limit t |r|-1 0 1 2 3 -0.04 -0.02 0 0.02 0.04 L RLLDDRK=1.64 RLDDDR=2.52LDDRK 0 1 2 3 -0.04 -0.02 0 0.02 0.04 t LDDRKLLDDRK=1.85L a) b) Figure 3-3. Four-sixstage optimized Runge-K utta of order four sc heme: a) dissipation error; b) phase error PAGE 89 75 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 fourth-order accurate scheme in time for a linear problem and second-order 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. four-stage 1 21 32 43 4 11 4 1 3 1 2n n n n nnKtFu KtFuK KtFuK KtFuK uuK (3.41) six-stage (1) 1 (2) 2 (3) 3 (4) 4 (5) 5 (6) 6 10.17667 0.38904 1 4 1 3 1 2n n n n n n nnKtFu KtFuK KtFuK KtFuK KtFuK KtFuK uuK (3.42) PAGE 90 76 In the follow we will give an implementa tion example of LDDRK scheme when we use OPC and DRP scheme for space discreti zation. Base on Eq.(3.1), the value of F in point l is defined as following: DRP fd: 3 3 lili ic Fau x (3.43) DRP-fv : Fl = c ( ul e ul w )/ x (3.44) where 12213415263 e llllllluuuuuuu (3.45) 13223145162 w llllllluuuuuuu (3.46) In the linear case, the fv and fd schemes are equivalent. OPC-fd 2BF lllc FDD (3.47) where Dl B and Dl F are obtained from the following system of equations: 1111 ,1,...,FF iiiiiiDDbuuduuiN x (3.48) 1111 ,1,...,BB iiiiiiDDduubuuiN x (3.49) where N represent the number of grid points in space. OPC-fv Fl = c ( ul e ul w )/ x (3.50) where ul e = 0.5( ul Be + ul Fe), and ul w = 0.5( ul Bw + ul Fw) (3.51) The value of ul Be, ul Fe, ul Bw, ul Fw are obtained by solving the follow system of equations 111,...,FeFe iiiiuubuduiN (3.52) 111,...,FwFw iiiiuubuduiN (3.53) 111,...,BeBe iiiiuubuduiN (3.54) 111,...,BwBw iiiiuubuduiN (3.55) PAGE 91 77 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 3332221111iiiiiiiFaxxaxxaxx x (3.56) This scheme requires a total of three multiplications and five additions to evaluate the first derivatives in a certain point. In cas e of DRP-fv the most efficient form of the computations scheme is 132221311eiiiiiiuuuuuuu x (3.57) 132221311wiiiiiiuuuuuuu x (3.58) DRPfv requires a greater number of opera tions than DRP-fd: 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 11111 222FF iiiiiiDbuuduuD x (3.59) 11111 222BB iiiiiiDbuuduuD x (3.60) PAGE 92 78 where the relation between the coefficients of the forward and backward stencils have been substituted to highlight the equivalent terms in the tw o stencils. The operation count is then four multiplications and five add itions per point (Ashcroft and Zhang, 2003). OPC-fv can be written in the form: 11 22FFeFw iiiDuu x (3.61) 11 22BBeBw iiiDuu x (3.62) where 11 111 1FeFe iiii F FwFw iiiiubuduu ubuduu (3.63) 11 111 1BeBe iiii BwBe iiiiubuduu ubuduu (3.64) In this case the operation count is eleven additions and six multiplications per point. So we can see also, in Tabl e 3-2 the finite volume a pproach is computationally more expensive. Dispersion Characteristics The characteristics of the OP C and DRP schemes, in the finite difference form over the interval 0 to are shown in Figure 3-4. 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 x < 1.30 for the DRP scheme, and x < 1.84 for the OPC scheme. The disp ersive characteristics of these PAGE 93 79 schemes can be more clearly seen in Figur e 3-5, which shows the phase speed error, 1 dx abs dx as a function of wave number on a l og-arithmetic scale. We see that the DRP scheme has a somewhat larger e rror than the OPC scheme until around 3 /4. The error is maintained to be within 2% for x less than 0.85 for the DRP scheme, and less than 1.53 for the OPC scheme. Overall, the OP C scheme yields sli ghtly less dispersion error than the DRP scheme. The dispersive characteristics of LDDR K are obtained by studying the value of | r | and i.e., dissipation rate and dispersion erro r (see Eq.(3.39)), resp ectively. In Figure 33 we can see conditions of stability: | r | < 1 for t 2.52 .To obtain an accurate solution the dispersive characteristics (| r | and ) should be close to the exact solution (| r | close to one and close to zero). Hu et al. (1996) considered time accurate criterion || r | -1| 0.001 (i.e t 1.64), and 0.001 (i.e., t 1.85). These two conditions are satisfied if t 1.64. Stability The Fourier-Laplace transformation of the wave equation (Eq.(3.1)) is 1 2initialiuciuu (3.65) where characterize the PDE For the long wave we can approximate wavenumber of the scheme with wavenumber of the PDE (3.66) which leads to initialucuku (3.67) PAGE 94 80 Hence *ckk (3.68) 0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5 Ideal DRP OPC Figure 3-4. Dispersive charac teristics of the schemes:. x versus x 0 0.5 1 1.5 2 2.5 3 3.5 10-4 10-3 10-2 10-1 100 101 x DRP OPC Figure 3-5. Phase speed error on a logarithmic scale The condition of the numerical stability is that amplification factor for time discretization is le ss than 1, and hence t 2.52 (see Figure 3-3a). It is also noted from Figure 3-4 that 1.8for DRPscheme 2.1for OPCscheme x x (3.69) Phase speed erro r PAGE 95 81 hold true. By introducing Eq. (3.69) into (3.68) and upon multiplying by t it is found that 1.8M1for DRP scheme 2.1M1for OPC scheme c tkkt x c tkkt x (3.70) where M is mach number. From Figure 3-3a it is clear that the condition of stability is satisfied if | t | is less than 2.52. Therefore to ensu re numerical stability it is sufficient by Eq. (3.70) to restrict t to be less than tmax, where tmax is given by max max2.52 for DRP scheme 1.8M1 2.52 for OPC scheme 2.1M1 x t kkc x t kkc (3.71) Therefore, for t < tmax the schemes are numerically stable. Consequently, the schemes yield the following criteria for numerical stability: 1.4for DRP scheme 1.2for OPC scheme CFL CFL (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 soluti on. In the previous analysis we have established that the solution is time accurate for 4-6 LDRRK if ||r|-1| 0.001 and | | 0.001. But this limit is not fixed, but depe nds 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 ||r|-1| 0.02 and | | 0.02, or t 2.0. Hence, in this case the condition of being both accu rate and stable is CFL 1.1 PAGE 96 82 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 one-dimensional wave equation with constant speed. The purpose of this test is to check the accuracy, stability, di ssipation and dispersi on of the scheme. The second test problem is one-dimensional nonlinear wave equati on with no viscous dissipation. The purpose of this test case is to i) check the influen ce of singularities on the performance of the scheme, and ii) anal yze dispersion properties when waves are coupled. In the third test problem, we c onsider the one-dimensional viscous Burgers equation, which contains unsteady, nonlinear convection and viscous terms. In this case we pay attention to the influence of the vi scosity 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: One-Dimens ional Linear Wave Equation To assess the behavior of the DRP and OPC schemes the following simple test problem is studied first. 0. uu c tx (3.73) 2 0expln2 xx u r ; at t = 0 (3.74) PAGE 97 83 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 2 0expln2 x xct U r (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 / x. For time discretization, we previously presented the detailed formulas for the 4-6 LDDRK, see equations (3.41), and (3.42). Tam (Tam and Webb,1993; Tam et al., 1993) show that x is related to x and in function of x they divided the wave spectrum into two categories; i) the long waves (waves for which x in this case x is less than xc, ii) the short waves (waves for which is not close to ). This difference between l ong and short waves is totally dependent upon the grid space. Hence, by insp ecting 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 / x = 20) intermediate wave ( r / x = 6) short wave ( r / x = 3) The categories are defined according to the ratio r / x where r is a parameter that characterizes the wavelength of this problem. Th is test problem is linear; hence we do not expect differences between finite di fference and finite volume approach. PAGE 98 84 In regard to the time step selection, the CFL number () limit is similar for all schemes. We can see from Fi gure 3-6 that the critical CF L number of both schemes is close to 1.1. From the study of the error in time for lin ear equations with constant convection speed it is clear that the DRP-fv and OPC-fv schemes have essentially the same behavior as the corresponding finite difference approach; he nce, we only present comparison for DRPfv and OPC-fv 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 confirme d in Figure 3-7, where the CFL number is maintained at 0.5. For the long time scale solution, the accumula tion of error for both DRP and OPC schemes is very close (as seen in Fi gure 3-7b). 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. t ERROR 0.1 0.2 0.3 0.4 10-4 10-3 10-2 DRP-fd DRP-fv OPC-fd OPC-fv t ERROR 0.05 0.1 10-5 10-4 10-3 DRP-fd DRP-fv OPC-fd OPC-fv a) r / x = 3 (short wave); b) r / x = 10 long wave Figure 3-6. Errors with respect to the time step size under a fixed space x at t = 50 linear wave equation PAGE 99 85 x ERROR 10-2 10-1 100 10-7 10-6 10-5 10-4 10-3 10-2 DRP-fv DRP-fd OPC-fv OPC-fd 4 1 t ERROR 1000 2000 0.001 0.002 0.003 0.004 1 1 DRP OPC a) b) Figure 3-7. Errors under a fixed CFL = 0.5, at t = 50 linear wave equa tion: a) error with respect to the space size; b) accu mulation of the error in time Test problem 2: One-Dimensio nal Nonlinear Wave Equation The finite volume and finite difference schemes are equivalent for a linear equation. The difference between them appear s for the nonlinear convective equation. To observe the merits and similarities of DRP, a nd OPC schemes, we restrict ourselves to the 1-D case. In this test, a nonlinear wave equati on with a different speed is solved 0 uu u tx (3.76) This equation is solved in the conservative form 20.50 u u tx (3.77) To better understand the effect of high gr adients and discontinuities, we chose the following initial conditions 00 (,0) 10 x ux x (3.78) The solution for this problem can be written PAGE 100 86 00 (,)0 1 x x uxtxt t x t (3.79) In this case, for both DRP and OPC scheme s, 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: DRP-fd 2 3 30.5ikik kFau (3.80) DRP-fv 220.5ew iiiFuu (3.81) where ue and uw are as defined before OPC-fd 0.25BF iiiFDD (3.82) where Di B and Di F is backward and forward derivative of u2 in place of u OPC-fv 2 20.25iew iiFuu (3.83) where ue and uw are defined by (3.52) (3.55) The similarities and differences for all three categories (short waves [ x / U = 1.0], intermediate waves [ x / U = 0.25], and long waves [ x / U = 0.06]) are first presented. It should be noted again that the short, interm ediate and long waves are defined based on the numerical resolution. Here, U is defined as the jump ( umax umin); in our case U = 1, hence in the following we discuss only the effect of the grid space step ( x ). PAGE 101 87 ERROR 10-2 10-1 100 10-3 10-2 DRP-fv DRP-fd OPC-fv OPC-fdx2 1 Figure 3-8. Errors with resp ect to the space step size unde r a fixed CFL = 0.5, at t = 5; nonlinear wave equation The evolution of the error as a function of grid spacing ( x ) is similar for both DRP and OPC schemes; the difference between the fi nite volume and finite difference versions are far greater, as shown in Figure 3-8. In th e case of finite volume, error decreases with decreasing grid space (see Figure 3-10 and Figur e 3-12). For finite difference, a totally different behavior is seen. The error not only does not decrease when grid spacing decreases, but in fact increases, as seen in Figure 3-9 and Figure 3-11). X U -4 -2 0 2 4 6 8 10 -0.2 0 0.2 0.4 0.6 0.8 1 DRP-fd Exact X=1. X U -4 -2 0 2 4 6 8 10 -0.5 -0.25 0 0.25 0.5 0.75 1 DRP-fd Exact X=0.25 X U -4 -2 0 2 4 6 8 10 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 DRP-fd Exact X=0.06 a) x = 1 b) x = 0.25 c) x = 0.06 Figure 3-9. DRPfd solution nonlinear wave equation; t = 5; CFL = 0.5 PAGE 102 88 X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 DRP-fv Exact X=1. X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 DRP-fv Exact X=0.25 X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 DRP-fv Exact X=0.06 a) x = 1 b) x = 0.25 c) x = 0.06 Figure 3-10. DRPfv solution nonlinear wave equation; t = 3; CFL = 0.5 For short waves, all solutions show substantial errors, but the finite difference schemes perform noticeably worse. In the case of intermediate or long waves, the finite volume schemes exhibit satisfactory or bette r performance than the finite difference schemes. X U -4 -2 0 2 4 6 8 10 -0.25 0 0.25 0.5 0.75 1 OPC-fd Exact X=1. X U -4 -2 0 2 4 6 8 10 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 OPC-fd Exact X=0.25 X U -4 -2 0 2 4 6 8 10 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 OPC-fd Exact X=0.06 a ) x = 1 b) x = 0.25 c) x = 0.06 Figure 3-11. OPC-fd solution nonlinear wave equation; t = 5; CFL = 0.5 X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 OPC-fv Exact X=1. X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 OPC-fv Exact X=0.25 X U -4 -2 0 2 4 6 8 10 0 0.25 0.5 0.75 1 OPC-fv Exact X=0.06 a) x = 1 c) x = 0.25 d) x = 0.06 Figure 3-12. OPC-fv solution nonlinear wave equation; t = 5; CFL = 0.5 PAGE 103 89 Test problem 3: One-Dimensional Nonlinear Burgers Equation In this test the solution for the onedimensional nonlinear Burgers equation is evaluated. 2 2uuu u txx (3.84) The numerical solution will approach e quation (3.84) in conservative form; 22 20.5 uuu txx (3.85) The initial condition is: 0,01tanh 2 x x ux (3.86) In this case the exact solution is: 0,1tanh 2 x xt uxt (3.87) The scheme described earlier for inviscid Bu rgers equation can also be applied to the current equation. This is accomplishe d by simply adding a second-order centraldifference expression for the viscous term uxx In other words Fi is replaced by Hi 11(2)/iiiiiHFuuux (3.88) Because of the viscosity that characterizes the scheme in this case, it is expected that the solution of both approaches would be stable and similar. Hence this term will have a large influence over the value of the error. In our discussion, we will distinguish the following three categories of results: short wave ( x/ = 10) intermediate wave ( x / = 3 ) PAGE 104 90 long wave ( x/ = 1) In this case the numerical performance is affected by two parameters: the CFL number and the Peclet number (Pe = U x /). First we compare the solution of all four schemes as function of the Peclet number (Pe) under constant CFL number (0.2). The valu e of CFL number is fixed at 0.2, because the critical value for all schemes is much lower in the present case than for the linear case. The behavior of the error is simila r among DRP-fv, DRP-fd, OPC-fd and OPC-fv: the error increases with incr easing Peclet number, until a certain value beyond which the schemes can no longer perform satis factorily (see Figure 3-13). Pe ERROR 100 101 10-4 10-3 10-2 10-1 DRP-fv DRP-fd OPC-fv OPC-fd* * * * Figure 3-13. Error as a function of Pe nonlinear Burgers equation; x = 0.25; CFL = 0.2; t = 20 For the four schemes (DRP-fv, DRP-fd, and OPC-fv, OPC-fd), the solution and error are very similar for all categories of wa ve, as shown in Figure 3-13, Figure 3-14 and Figure 3-15. For long waves the solution is re produced with high accura cy with all four schemes, but the finite volume approach pres ents a slightly higher accuracy than the finite difference schemes. The error for the inte rmediate wave is nearly the same with all four approaches. PAGE 105 91 X U 16 18 20 22 0 0.5 1 1.5 2 2.5 DRP-fv DRP-fd Exact X U 16 18 20 22 0 0.5 1 1.5 2 DRP-fv DRP-fd Exact X U 16 18 20 22 0 0.5 1 1.5 2 DRP-fv DRP-fd Exact a) Pe = 10 b) Pe = 3 c) Pe = 1 Figure 3-14. Numerical solution obtained by DRP schemes nonlinear Burgers quation; x = 0.25; CFL = 0.2; t = 20 X U 16 18 20 22 0 0.5 1 1.5 2 2.5 DRP-fv DRP-fd Exact X U 16 18 20 22 0 0.5 1 1.5 2 OPC-fv OPC-fd Exact X U 16 18 20 22 0 0.5 1 1.5 2 OPC-fv OPC-fd Exact a) Pe = 10 b) Pe = 3 c) Pe = 1 Figure 3-15. Numerical solution obtained by OPC schemes nonlinear Burgers equation; x = 0.25; CFL = 0.2; t = 20 Test problem 4: Two-Dimensiona l Acoustic Scattering Problem To check the accuracy of the finite volume schemes in multi-dimensional situations, we consider a test problem from the Second CA A Workshop (Tam and Webb,1993; Tam et al., 1993): the two-dimens ional acoustic scattering problem. The physical problem is to find the sound field generated by a propeller scattered off by the fuselage of an aircraft. The pressure loading on the fuselage is an input to the interior noise problem. The fuselage is idealized by a circular cylinder and the noise source (propeller) as a line source so that the co mputational problem is two-dimensional. The cylinder has a radius of R = 0.5 and is located at the center of the domain. The linearized Euler equations in polar coordinates are: PAGE 106 92 000 11 000 0r rru p up trrr uuu p (3.89) At time t = 0, the initial conditions are: ur = u = 0 (3.90) 2 24 (,,0)expln2 0.04 xy pxy (3.91) The test problem asks for the unsteady pressure time history at three points A( r = 5, = 900), B( r =5, = 1350) and C( r = 5, = 1800), over the interval t = 5 10. The numerical computations were performed over the domain: R [0.5, 10.5] and [0, 2 ]. For this problem three kinds of the boundary conditions are needed: Wall condition on the wall of the cylinder at R = 0.5 Periodic condition along both azimuthally boundaries at = 0 and = 2 Outfield boundary condition, along of th e far field boundary, is the acoustic radiation of Bayliss a nd Turkel (1982)]. The wall condition is based on the wall condition of Tam and Dong (1994). This requires that: 0rdvdp dtdr (3.92) This condition is satisfied by imposing the pres sure derivatives on the wall to be zero, and vr = 0 on the wall. PAGE 107 93 X Y -10 -5 0 5 10 -10 -5 0 5 10 X Y -10 -5 0 5 10 -10 -5 0 5 10 a) DRP-fv b)OPC-fv Figure 3-16. Instantaneous pr essure contours at time t = 7 two-dimensional acoustic scattering problem x Y -10 0 10 -10 -5 0 5 10A B C time P 6 8 10 -0.02 0 0.02 0.04 0.06 OPC-fv DRP-fv Exact a) position of the testing points b) A: R = 5 = 900 time p 6 8 10 -0.02 0 0.02 0.04 0.06 OPC-fv DRP-fv Exact time P 6 8 10 -0.02 -0.01 0 0.01 0.02 0.03 0.04 OPC-fv DRP-fv Exact c) = 1350 d) = 900 Figure 3-17. The pressure hi story at point A, B and C two-dimensional acoustic scattering problem: finite volume approach PAGE 108 94 For this calculation, a uniformly spaced grid of 101 radial points and 153 azimuthally points was used, with a time st ep of CFL = 0.5. Figure 3-16 shows an instantaneous pressure at t = 7. In this figure, the ac oustic pulse is reflected by the cylinder and reaches the outer boundary. We ca n see that two transients are shown: the first and larger transient trav els directly from the source; the second and smaller transient is reflected from the cylinder. Both sche mes reproduce both transients with acceptable accuracy. Figure 3-17 compares the solution gi ven by the fourth order schemes: DRP-fv and OPC-fv. Between the two schemes the OPC-fv scheme performs better. Summary and Conclusions The DRP and OPC schemes, originally propos ed in the finite di fference form, have been assessed. To better handle nonlineari ty and geometric complexities, the finite volume versions of both schemes have also been developed. Linear and nonlinear wave equations, with and without visc ous dissipation, have been adopt ed as the test problems. For the linear wave equation with cons tant convection speed, the numerical stability bound posed by the CFL number is comparable between the DRP and OPC schemes. Both OPC and DRP produce solutions of a comparable order of accuracy, but the magnitude of the error of the OPC scheme is lower. For the nonlinear wave equation, the fin ite volume schemes can produce noticeably better solutions and can handle the discontinu ity or large gradients more satisfactorily. However, as expected, all schemes have di fficulties when there is insufficient mesh resolution, as reflected in some of the short wave cases. In conclusion, the finite volume versi on of both DRP and OPC schemes improve the capabilities of the original version of th e finite difference formulas in regard to nonlinearity and high gradients. They can e nhance performance of the original DRP and PAGE 109 95 OPC schemes for many wave propagation problems enc ountered in engineering applications. PAGE 110 96 Table 3-1. A summary of proposed CAA algorithms Scheme The philosophy of the scheme Applications DRP (Tam and Webb, 1993; Tam et al., 1993) In this scheme a central difference is employed to approximate first derivative. The coefficients are optimized to minimize a particular type of error Wave propagation LDDRK (Hu et al., 1996 ; Stanescu et al., 1998) Traditionally, the coefficients of the RungeKutta scheme are optimized to minimize the dissipation and propagation waves. Wave propagation problem LDFV (Nance et al, 1997; Wang et al., 1999) Scheme minimizes the numerical dispersion errors that arise in modeling convection phenomena, while k eeping dissipation errors small: special high-order polynomials that interpolate the properties at the cell centers to the left and right sides of cell faces. A low pass filter has been implemented Shock noise prediction GODPR (Cheong and Lee, 2001) It is derived, base d on optimization that gives finite difference equations locally the same dispersion relation as the original partial differential equations on the grid points in the nonuniform Cartesian or curvilinear mesh Acoustic radiation from an oscillating circular cylinder in a wall Scattering of acoustic pulse from a cylinder OWENO (Wang and Chen, 2001) The idea is to optimize WENO in wave number schemes, following the practice of DRP scheme to achieve high resolution for short wave. But in the same time it retains the advantage of WENO scheme in that discontinuity are captured without extra numerical damping. Simulation of the shock/broadband acoustic wave CE/SE (Chang, 1995; Chang et al., 1999) The method is developed by imposing that: (i) space and time to be unified and treated as a single entity; (ii) local and global flux conservation in space and time to be enforced; Flow involving shock; acoustic wave FDo, RKo (Bogey and Bailly, 2004) Optimized schemes are obtained by similar approach as DRP (space discretization), respectively LDDRK (time discretization). The difference consists that: i) error is minimized taking into account logarithm of the wavenumber; ii) the error is minimized on an interval that starts from ln( /16). The stability and accuracy increase for these schemes a)convective wave equation b)subsonic flows past rectangular open cavities; c) circular jet PAGE 111 97 Table 3-2. The computational cost for DRP and OPC schemes Scheme Number of operation DRP-fd 8 DRP-fv 17 OPC-fd 9 OPC-fv 17 PAGE 112 98 CHAPTER 4 A FINITE VOLUME-BASED HIGH OR DER CARTESIAN CUT-CELL METHOD FOR COMPUTATIONAL AEROACOUSTICS A finite volume-based high-order scheme with optimized dispersion and dissipation characteristics in cooperation with the Cart esian cut-cell technique is presented for aeroacoustic computations involving geometri c complexities and non linearities. The field equation is solved based on an optimized pr efactored compact finite volume (OPC-fv) scheme. The cut-cell approach handles the boundary shape by sub-dividing the computational cells in accordance with the lo cal geometric characteri stics and facilitates the use of numerical procedures with a desirable level of accu racy. The resulting technique is assessed by severa l test problems that demonstr ate satisfactory performance. Introduction In computational aero-acoustics (CAA), one needs to resolve dispersion and dissipation characteristics in order to preser ve form and amplitude of the wave (Hardin and Hussaini (1992), Tam and coworkers (T am and Webb, 1993; Tam et al., 1993)). Furthermore, special care needs to be exer cised in treating the boundaries to prevent the creation of spurious waves while still accounti ng for wave reflection and/or scattering. In order to meet these expectations, there are sp ecific issues and chal lenges associated with the employment of a structured multiblock grid, an overlapping grid, or an unstructured grid ((Delfs, 2001), (Shyy et al., 2001), (Henshaw, 2004), (Sherer, 2004), (Dumbser, 2003), (Basel and Grnewald, 2003 )), other alternatives should also be pursued. In this work, we present an approach utilizing the Cartesian cut-cell approach based on the PAGE 113 99 finite-volume based scheme aimed at op timizing the dispersion and dissipation treatment. The Cartesian, cut-cell approach has been developed extensively for moving and complex boundary computations ( see Table 4-1). It uses a background Cartesian grid for the majority of the flow domain while creating sub-cells in th e boundary regions to meet the geometric requirements. Special algorithms for complex geometries are constructed using flow properties and appr opriate interpolation and extrapolation procedures. The low dispersion and dissipation scheme is based on an extension of the finite difference-based, optimized prefactored comp act (OPC) scheme originally developed by Ashcroft and Zhang (2003) that we call OPC-fd. In th e previous chapter we present the extension of this method using the finite vol ume concept that we call OPC-fv, creating a scheme to better satisfy the nonlin earity and conservation laws. The present technique (which consists of the optimized disper sion and dissipation characteristics and the cut-cell technique) offers desirable capabilities in both interior and boundary regions. Several test problems will be presented to highlight the performance characteristics of the proposed approach. Cut-Cell Procedure The cut-cell method rearrang es the computational cells in the vicinity of the interface via sub-division to conform to the specified boundary shape. Depending on the intersection between the grid line and the interface, the subdivided, or cut, cells can remain independent or can be merged into a neighboring cell in a gi ven direction, e.g., a direction approximately normal to the solid face. Accordingly, the interfacial cells are reorganized along with their neighboring cells to form new cells with triangular, trapezoidal, or pentagonal shapes. PAGE 114 100 The flux across the cell boundaries can be approximated by 1m kk k S f nfn (4.1) where S f n is contour integral on the closed path S The flux on the face is computed based on the multi-dimensional interpolation method (Ye et al., 1999; Ye, 2001). Since the OPC scheme considered here is fourth order, it is desirable to preserve the same order of accuracy around the boundary. Figure 4-1. Illustration of the interfacial cells and cut-and-absorption procedures The Cartesian cut cell mesh approach follows the subsequent steps: Locate the intersection of the boundary using a Cartesian mesh. Construct the background Cartesian mesh: Th e cells are flagged as solid cells, flow cells, or boundary cells. The boundary cells are those that eith er intersect the boundary or have a face in common with the boundary. Determine the geometric characteristics of the boundary cells su ch as cell volume, the direction normal to the boundary, and other information. Merge cells as necessary. A minimum accepta ble cell area Smin is specified, and when the area of a cut cell is smaller than this value, it is merged into a neighboring cell. Determine the new characteristics of the merged cells. interface interface domain 1 domain 2 domain 2 domain 1 PAGE 115 101 A B CE D 29 28 27 26 25 19 22 21 20 23 31 33 32 37 35 34 13 16 15 14 17 7 10 9 8 11 1 4 3 2 5 6 12 18 24 30 36 A B CE D 29 28 27 26 25 19 22 21 20 23 31 33 32 37 35 34 13 16 15 14 17 7 10 9 8 11 1 4 3 2 5 6 12 18 24 30 36 BD AE C fswfe fw G (xi,yj) BD AE C fswfe fw G (xi,yj) a) b) Figure 4-2. Modified cut cell approach for CAA: a) Cartesian cut cell approach; b) Detail around of the cut cell For the trapezoid ACDE from Figure 4-2, the finite volume approach can be used to approximate the wave equation: 0ABCDEuEF dv txy (4.2) where dv is a volume element. Appling Stokes s theorem to Eq.(4.2), we have uc/ t SABCDE + ABCDE E dyFdx (4.3) The value of the integral can be approximated by ABCDEACCDDEEA E dyFdxEdyEdyFdxEdyFdx (4.4) where ABCDEEdyFdx is contour integral on closed path ABCDE. The function F and E can be represented ge nerically by the function f The fluxes on the faces AC and DE can be approximate by PAGE 116 102 ACABBC f dyfdyfdy (4.5) An approximation of the value of the flux is given by: wABswBC AC f dyfyyfyy (4.6) The value of the flux at point w is given by the specific formula for the boundary cell. The value of fsw is approximated using a polynomial that is fourth order in the x and y direction: 44 00ij swij ij f bxy (4.7) where the coefficients bij are unknown. This interpolation has fourth-order accuracy in the evaluation of the flux on the cut cells. In this case, the value of the coefficients is obtained using the values of f in 25 grid points. An example is given in Figure 4-2 where the value of the function is approx imated using the following 25 points. To solve bij we use the following system of equations by expressing the function f at 25 locations: 443443 11111111 11 443443 22 22222222 443443 2525 25252525252525251 1 1 xyxyxyxy f b f b xyxyxyxy f b xyxyxyxy (4.8) The coefficients bij from Eq.(4.7) become the coefficients b1, b2, b25 in Eq. (4.8) These coefficients can now be expressed in terms of values of f at the twenty five points by inverting Eq.(4.8), i.e., 25 1,1,...25nnjj jbafn (4.9) PAGE 117 103 where anj are the elements of the invers e of the matrix in Eq.(4.8). After bn is obtained, the value of f at the center of BC is expressed in the form 443443 123232425swswswswswswswswsw f bxybxybxybxbyb (4.10) and using Eq.(4.9), the value of fsw can be rewritten as 25 1 s wjj j f f (4.11) Note that bi i = 1, 2, 25 are coefficients that de pend only on the mesh, the location, and the orientation of the boundary. Therefore, with a fixed geometry, these coefficients can be computed once at th e beginning of the solution procedure. Now we turn to the calculation of the flux on the immersed face CD of the cell ( i, j ). To compute the flux on a solid boundary, we use the reflection boundary condition. In Figure 4-2b, ( xi, yj) is the mass center of the boundary cell. We introduce a point across the boundary that is symmetrically oppos ite of the center of mass of face CD and note it with G The variable values in the ghost point G are 2Gij Gijijpp uuunn (4.12) where n is the normal vector of the solid boundary, the index ij indicates the value of variable at point ( xi, yj), and the index G indicates values at the ghost point G The value of the flux on face CD is approximate by 2 2Gij CD Gij CD p p p uu u (4.13) PAGE 118 104 Following simple algebra, we can observe that Eq.(4.13) assures that the velocity is zero in the normal direction, and the variation of pressure in the normal direction around the boundary is zero. Test Cases Three test cases have been selected to assess the performance of the present approach: i) radiation from a baffled pist on in an open domain; ii) reflection of an acoustic pulse on an oblique wall; iii) a wave generated by a baffled piston and reflected on an oblique wall. The numerical methods used in the test pr oblem are: OPC-fv for discretization in space, and LDDRK for discretization in time. Radiation from a Baffled Piston A piston in a large baffle is a good st arting approximation fo r investigating the radiation of sound from a boxed loudspeaker. The physical problem is to find the sound field generated by a piston that is baffled. The problem is solved in a two-dimensional planar domain. This problem is chosen to evaluate the field equation solver in the absence of the rigid wall. General description We will use a system of coordinate centered in origin of the piston as shown in Figure 4-3. Dimensionless variables with resp ect to the following scales are to be used Length scale = diamet er of piston, 2 a Velocity scale = speed of sound, c Time scale = 2 a/c Density scale = undisturbed density, 0 Pressure scale = 0c2 PAGE 119 105 Frequency = c /2 a The linearized Euler equations are used: 0 0. p u t u p t (4.14) The initial conditions are: ,0. ,0. ,0. uxy vxy pxy (4.15) y 2a y 2a x y x y x y a) b) Figure 4-3. Radiation from a ba ffled piston (test problem (ii)): a) general description; b) boundary condition For this problem, two boundary conditions are used. The boundary condition on the wall with piston is 0cos,0 ,0, 0pVtxpiston uxt otherwise (4.16) where is the frequency of the piston, and V0 is the amplitude of the displacement. The solution is obtained using linear equation; hence the value of V0 influences only the PAGE 120 106 amplitude of the solution and not its behavior From this reason we take the value of V0 equal to one. The outflow boundary condition is based on the acoustic radiation condition of Tam and Web (1994): 0. cossin0. 2 u p t pppp txyr (4.17) where is the angular coordinate of the boundary point, and r is distance from the origin (center of the piston ) to the boundary point. An analytical solution exists for this problem (Williams, 1999; Morse and Ingard, 1968): 22 0 22',0, ..' 2p suxt c p xytddx (4.18) where the x values are points on the source, 2 = ( x-x )2 + y2 is distance from a point of the piston, and pu is time derivative of the di splacement of the piston. The behavior of the piston is presente d in terms of the Helmholtz number, ka where k is wavenumber, and a is the radius of the piston a ka c (4.19) In the following we will pres ent two cases: low frequency ( ka = 2) and high frequency ( ka = 7.5). PAGE 121 107 Directive factor D The directional charac teristic of a source is descri bed by the amplit ude directivity factor D defined as the pressure at any angle to the pressure on the angle of maximum pressure. ,/2 pr D pr (4.20) where the pressure is com puted at any arbitrary time t It is clear that the radiation is strongest on the y axis, that is why we took maximum pressure at = /2. We note that r and are spherical coordinate. Because wavenumber k is in the same direction as r then the same spherical angles describe both of them. Thus in spherical coordinate we have: cos cos sin sinx ykk xr kk yr (21) It is clear that the strongest radiation is on the y axis. In case that 1ka the directive function has nulls, and between the nulls are secondary radiation maxima, of monotonically decreasing prominence. The num ber of nulls and secondary maxima is determined by the size of ka, where a = radius of the piston, and ka = 2 a/ In other words the number of the lobs increases with the value of ka. In our calculation we compare the analytical and numeri cal values of beam pattern for ka = 2 and ka = 7.5. As shown in Figure 4-2 and Figure 4-3, the numer ical solution recovers with high accuracy the analytical beam pattern, which proves th e capacity of the appr oach to recover the solution. PAGE 122 108 Figure 4-3 and Figure 4-4 represent, respectively, the beam pattern for ka = 7.5 and ka = 12.5. In these figures we can see that the number of lobs in creases in the same time with Helmholtz number (ka). Figure 4-4. Radiati on from a baffled pist on: Beam patterns |D( )|; ka = 2 Figure 4-5. Radiati on from a baffled pist on: Beam patterns |D( )|; ka = 7.5 Figure 4-6. Radiati on from a baffled pist on: Beam patterns |D( )|; ka = 12.5 PAGE 123 109 Pressure on the face of the piston In the following the radiance impedance of the piston is calcul ated. The radiation impedance of the piston is given as: av p av p Z u (4.22) where the average values of th e pressure and velocity on th e face of the piston are given by: /2 /2 /2 /2(,) (,)a av a a av a p pxydx uuxydx (4.23) The function R1, the resistance is defined by: 10Re/p R Zc (4.24) To contrast the behavior of resistance is in the doma in of low frequency (ka 1) and high frequency (ka 1), we can see in Figure 4-7 that the resistance is close to zero at low frequency, and close to one at high frequency. Ov erall, the numerical solution closely follows the analytical solution. Figure 4-7. Radiation from a baffled piston: Real of piston radiation impedance PAGE 124 110 Low frequency ( ka = 2) The computation is done on a uniformly spaced grid x = y = 0.05 and a time step t based on CFL = 0.2, where CFL = t c x The challenge of this problem is to solve the discontinuity at time t and at the points yct when maxcty where c is the speed of sound. In this test case, the variables are non-dimensionalized. Hence, c = 1, and the point of discontinuity will be at y = t, for t < ymax (see Figure 4-8). Figure 4-9 compares the solu tion given by the OPC-fv scheme with the analytical solution demonstrating that there is good agreement between them. Figure 4-8. Radiation from a baffled piston: Numerical solution on Oy axis (ka = 2), t = 7.3 < ymax Figure 4-10 shows contour plots of pressure at t = 20 and t = 40. The piston radiation starts as a collimated plane-wa ve beam that moves from the face to r radius of piston, beyond which the beam propagates sphe rically. This behavior is in accordance with Blackstocks study (2000). PAGE 125 111 Figure 4-9. Radiation from a baffled piston: Comparis on between analytical and computed solutions on axis (ka = 2) a) b) Figure 4-10. Radiation from a baffled piston: Contour plot of pressure (ka = 2) numerical solution at: a) t = 20.0; b) t = 40.0 High frequency ( ka = 7.5) The computation is done on a uniformly spaced grid x = y = 0.02 and a time step t based on CFL = 0.2, where CFL = t c x In this case we increase the resolution to better capture the wave characteristics. Fi gure 4-11 compares the solution given by the OPC-fv scheme with the anal ytical solution dem onstrating that there is good agreement between them. Figure 4-12 shows c ontour plots of pressure at t = 20, highlighting the entire domain as well as the region close to the piston. PAGE 126 112 Figure 4-11. Radiation from a baffled piston: Comparis on between analytical and computed solutions on axis (ka = 7.5) a) b) Figure 4-12. Radiation from a baffled piston: Contour plot of pressure (ka = 7.5) numerical solution at t = 20: a) entire domain; b) reduced domain around piston Reflection of a Pulse on an Oblique Wall A simple example that enables us to check the performance of the cut-cell approach is the reflected sound on an oblique wall. Even though this problem can be solved by placing the wall parallel to the grid line, we purposely arrange the wall so it is at an oblique angle to the grid line. This offers a direct evaluation of the cut-cell technique. PAGE 127 113 In this example, the sound hits the wall and reflects. The problem is characterized by the linearized Euler equations, Eq. (4.14). The initial condition is 22 00(,)expln(2) xxyy pxy bb (4.25) where b = 1./6. The values x0 and y0 are chosen so that the distance from the point (x0, y0) to the oblique solid wall is equal to 1.5. Figure 4-13. Reflection of th e pulse on an oblique wall (t est problem (ii)): general description The outflow boundary condition is the same as that used in the previous test problem. The solid wall boundary conditions for pressure and velocity are given by reflecting boundary condition. The computation is done on a Cartesia n grid that is characterized by x = y = 0.05 and CFL = 0.5. Figure 4-12 i llustrates the solid boundary intersecting a Cartesian mesh. The boundary cell can remain independent, like in ABCK, or it can merge with a neighboring cell. For example, CID merges with KCIJ. The computation of the parameter in the cut cell CDIJK can be written for pressure fr om wave Eq.(4.14) in the form: PAGE 128 114 1 0.CDCDJKJKDIDIIJIJCKCKCDCD CDIJKp udyudyvdxvdxvdxvdx tS (4.26) where the values of u and v represent the values of the function in the middle of the segment. The value of the parameter is approximated using i) the OPC scheme on faces IJ, JK and KC, ii) fourth order polynomial on faces DI and iii) refl ecting boundary condition on face DF (boundary). Figure 4-14. Reflection of the pulse on an oblique wall (tes t problem (ii)): cell around boundary To evaluate the performance of the cu t-cell approach, Fi gure 4-15 shows the pressure history from time = 0 4 and for three wall angles: 900, 810, and 630. As apparent from the figure, there is a good ag reement between these solutions. Figure 4-16 shows the pressure contours at three time instants by placing the solid wall at 630 inclinations. PAGE 129 115 Figure 4-15. Reflection of the pulse on an obliq ue wall: history of pressure for different angle of wall a) b) c) Figure 4-16. Reflection of th e pulse on an oblique wall: = 630 a) t = 0.8; b) t = 1.6; c) t = 3.2 Wave Generated by a Baffled Piston and Reflected on an Oblique Wall This problem is based on the combined ch aracteristics of the previous two cases. The wave generated by a piston is reflected by an oblique wall. The linearized Euler equations (Eq. (4.14)) is used in this test case, and the initial condition is the same as that of the baffled piston, Eq. (4.16). The rectangular domain over which we do the computation is: ,6,60,15 xy The bottom of the domain is a piston mounted on a plane rigid PAGE 130 116 baffle; hence, the boundary conditions are give n by Eq.(4.16). The velocity and pressure on the solid and open boundaries are identical to those used in the previous case. Figure 4-17. Wave generated by a baffled piston and reflects on an oblique wall: General description of the domain The piston presents the fo llowing characteristics: V0 = 1 and = 4 (ka = 2). The solution is obtained using the uniform grid x = y = 0.05, and with a time step of CFL = 0.5. Figure 4-18 highlights the pressu re contours at different time instants. a) b) Figure 4-18. Wave generated by a baffled piston and reflects on an oblique wall: = 630; a) t = 9; b) t = 14; The challenges of this problem are to recover: PAGE 131 117 The wave generate by the piston The reflection wave The combination of the two waves: wave generated by the piston and reflected on the oblique wall. Figure 4-18 shows the outcome of the com putation, exhibiting a series of local maxims and nulls. Conclusion A method based on OPC-fv scheme aime d at optimizing the dispersion and dissipation properties, and the cut-cell technique aimed at ha ndling geometric variations is presented. The approach is motivated by the need for handling acoustic problems with nonlinearities (using finite vol ume technique) and a complex ge ometry (using the cut-cell technique). Selected test cases have been used to demonstrate the performance of the method. The present approach can offer accu rate and versatile treatment to some important and challenging as pects of acoustic problems. PAGE 132 118 Table 4-1. Published cut-cell a pproach for different problems Authors Test Cases Objective De Zeeuw et al, 1999 ; Transonic NACA Airfoil Subsonic Three-Element Airfoil Supersonic Double Ellipse Supersonic Channel Flow A method for adaptive refinement of the steady equation Udaykumar et al, 1997 Deformation of viscous droplets in axisymetric Extensional Stakesian flow Deformation of Droplets in extensional flows with Inertia Effects Deformation of droplets in constricted tubes The ELAFINT algorithm is developed and applied to compute flows with solid-fluid and fluid-fluid interface Yang et al, 1999; Causon et al, 1999; Ingram et al, 2003 Open-Ended shock tube Fifteen-degree wedge flow at Mach 2 Muzzle brake flowfields Muzzle Sabot/Projectile flowfields Shallow water flow The flow is an upwind scheme of the Godunov-type based on MUSCL reconstruction and a suitable approximate Riemann solver Udaykumar et al, 1999 Inviscid flow around circular cylinder Track solid-liquid boundaries on a fixed underlying grid. The interface is treated as a discontinuity conditions and explicitly tracked Ye et al, 1999 Two-dimensional stokes flow past a circular cylinder Flow past a circular cylinder immersed in a freestream Flow past a circular cylinder in a channel Application to complex geometries: (i) flow past a random array of cylinders (ii) flow past a cascade of airfoils A Cartesian grid method has been developed for simulating two-dimensional unsteady, viscous, incompressible flows with complex immersed boundaries Lahur et al, 2000 Moving piston Moving cylinder Treat moving body problem D.Calhoun et al, 2000 Advection and diffusion of a plane wave in a channel Advection and diffusion in an annulus Advection through a field of irregular objects A fully conservative, highresolution, finite volume algorithm for advectiondiffusion equations in irregular geometries Verstappen et al., 2000 Flow over a circular cylinder To introduce a novel cut-cell Cartesian grid method that preserve the spectral properties of convection and diffusion Ye et al, 2001 Bubble dynamics with large liquid-to vapor density ratio and phase change Treating sharp discontinuity interface for bubble dynamics and phase change PAGE 133 119 CHAPTER 5 SUMMARY AND FUTURE WORK Predicting noise radiation associated with unsteady flows is the central theme in aeroacoustics. Since the generation of the s ound by the flow typica lly involves complex physical mechanisms, approximate models are often employed. Colonius and Lele (2004) have reviewed the modeling as well as the numerical aspects based on recent progress. Although it is not possible to se parate the flow from sound, oftentimes, we can separate the flow prediction from noise prediction, based on pressure, fre quency and length scale separations. In these regards, CAA presents a series of challenges. In particular, it is important to extract small quantities relative to other phenomena in the course of computations. For example, the pressure fluctuations in a turbulent flow are on order of 10-2 of the mean pressure (Blake, 1986). However, loud noise can be produced by pressure fluctuation at the level of 10-5 of the mean pressure (Hardin, 1996). In order to accurately capture the acoustic field, numerical techniques need to yi eld low artificial dissipation. In addition, the numerical dispersion related to frequency and phase constitutes another concern. For this reason, many tests of goodness of CAA so lutions have been developed. In light of these issues, we have inves tigated schemes aimed at offe ring refined performance in numerical dissipation and disp ersion, including the disp ersion-relation-preservation (DRP) scheme, proposed by Tam and Webb (1993), the space-time amethod, developed by Chang (1995), and the optimi zed prefactored compact (OPC) scheme, proposed by Ascroft and Zhang (2003). The space-time amethod directly controls the PAGE 134 120 level of dispersion and dissi pation via a free parameter, the DRP scheme minimizes the error by matching the characteristics of the wave, and the OPC scheme aims at employing a compact stencil to achieve optimized outcome. Insight into the dispersive and dissipative aspects in each scheme is gained from analyzing the truncation error. Another challenge in CAA consists in the reproduction of the nonlinear propagation, including nonlinear steepening an d decay, viscous effects in intense sound beams (e.g., acoustic near field of high-sp eed jets), sonic boom propagation, hydrofoil and trailing-edge noise at low speed (G oldstein, 1976; Wang and Moin, 2000; Manoha and Herraero, 2002; Takeda et al, 2003; Howe 1978, 1999, 2000, 2001; Casper and Farassat, 2004; Roger and Mor eau, 2005) and high speed jet noi se ( Davies et al, 1963; Hong and Mingde, 1999; Shur et al, 2003; Fisher et al, 1998; Morris et al, 1997; Bogey et al, 2003). In this regard, we have develope d finite volume treatments of the DRP and OPC schemes. These schemes, originally ba sed on the finite difference, attempt to optimize the coefficients for better resolution of short waves with respect to the computational grid while maintaining pre-dete rmined formal orders of accuracy. In the present study, finite volume formulations of both schemes are presented to better handle the nonlinearity and complex geometry enc ountered in many engin eering applications. Linear and nonlinear wave equations, with and without viscous dissipation, have been adopted as the test problems. Highlighting the principal characteri stics of the schemes and utilizing linear and nonlinea r wave equations with differe nt wavelengths as the test cases, the performance of thes e approaches is documented. For nonlinear wave equations, the finite volume version of both DRP a nd OPC schemes offer substantially better solutions in regions of high gradient or discontinuity. PAGE 135 121 Boundary treatment is another critical issue. Instead of employing boundary conforming grid techniques, we propose a Cartesian grid, cut-cell approach, aimed at aero-acoustics computations involving geomet ric complexities and nonlinearities. The cut-cell approach handles the boundary shape by sub-dividing the computational cells in accordance with the local geometric characteri stics and facilitates the use of numerical procedures with desirable accuracy. The resu lting technique is assessed by several test problems that demonstrate satisfactory performance. In the following we summarize the outcome of the present research, as presented in the previous chapters Assessment of DRP and Space-Time CE/SE Scheme First we investigate two num erical methods for treating convective transport: the dispersion-relation-preserving (DRP) by Ta m and Webb (1993), and the unified spacetime amethod, developed by Chang (1995). Ba sed on the investig ation of the DRP CE/SE scheme we offer the following summary: The two methods exhibit different perfor mance with regard to the CFL number, Both truncation error analysis and numerical testing indicate that better solutions can be obtained for the DRP scheme method if is close to 0.2 for short wave, and close to 0.1 for intermediate and long wave. For the space-time amethod it is preferable if is less than but close to 1. Both schemes are dispersive for short waves. The space-time amethod directly controls the level of di spersion and dissipation, via the free parameter, The DRP scheme, on the other hand, minimizes the error by adjusting the scheme to match the characteristics of the wave. Formally, based on the truncation error analys is, the DRP scheme is fourth order in space and third order in ti me, while the space-time ascheme is second order. Evidence based on the truncation error analys is and the test pr oblems indicate that in order to reduce numerical dispersion a nd maintain satisfactory resolution, for short wave ( e.g., b/ x = 3), and are preferable to be close to each other. For long and intermediate waves there is virtua lly no need to introduce much numerical PAGE 136 122 dissipation and hence can be reduced close to zero. On the other hand, mismatched and can substantially worsen the pe rformance of the scheme. Hence to achieve the best performance of th e CE/SE scheme, one chooses the largest possible and then a matching The DRP scheme exhibits mainly dispersive errors, while the space-time ascheme exhibits both disper sive and dissipative errors. It is advisable to use space-time ascheme for long wave computation because its error grows slower. For short and intermediate waves, the DRP scheme yields errors with lower level and slower accumulation rate. DRP does not offer an accurate solution for a short wave, x greater than cx (in our case greater than 1.1), because the wave length is not adequately to condition of the method. Hence the DRP scheme offers a good guidance in regard to temporal and spatial sizes. It should also be noted that he DRP sche me is a multi-step method, which requires more boundary conditions and initia l data, while the space-time ascheme is a one-step method. Combined with the fact that the DRP scheme performs better with smaller ( = CFL), it can be more expensive to compute than for space-time ascheme. Finite-Volume Treatment of Dispersi on-Relation-Preserving and Optimized Prefactored Compact Schemes The DRP and OPC schemes, originally propos ed in the finite di fference form, have been further developed. To better handle nonlinearity and geometric complexities, the finite volume version of both schemes has al so been constructed. Linear and nonlinear wave equations, with and wit hout viscous dissipation, have been adopted as the test problems. For the linear wave equation, the nu merical stability bound posed by the CFL number is comparable between the DRP a nd OPC schemes. Both OPC and DRP schemes PAGE 137 123 produce solutions of comparable accuracy, but the magnitude of the error of the OPC scheme is somewhat lower. For the nonlinear wave equation, the fin ite volume schemes can produce noticeably better solutions and can handle the discontinu ity or large gradients more satisfactorily. However, as expected, all schemes have di fficulties when there is insufficient mesh resolution, as reflected in some of the short wave cases. In conclusion, the finite volume versi on of both DRP and OPC schemes improve the capabilities of the original version of th e finite difference formulas in regard to nonlinearity and high gradients. They can e nhance performance of the original DRP and OPC schemes for many wave propagation problems enc ountered in engineering applications. Cartesian Grid, Cut-Cell Approach for Complex Boundary Treatment A method based on high order, finite-vol ume schemes aimed at optimizing the dispersion and dissipation proper ties, and the Cartesian grid, cut-cell technique aimed at handling geometric complexities is presented. The approach is motivated by the need for handling acoustic problems for practical pr oblems involving realistic geometries. The finite volume-based Optimized Prefact ored Compact scheme and the Cartesian cut-cell approach are combined to offer 4th order accuracy and geometric flexibility. The computational overhead of the cut-cell a pproach is modest because the following information needs to be computed only once, unless, of course, if the geometry is time dependent: Data communication between cells affected by the boundary treatment Calculation of area and ot her geometric information Interpolation procedures required for th e flux computation in the boundary region. PAGE 138 124 Based on the evaluation of the test cases investigated, we conclude that the present approach can be effective in treating aero acoustics problems with irregular geometry. Future Work Even though the radiant pressure fluctuation is less than 10-4 of the ambient pressure, there are situations when the non linearity behavior can not be neglected. Examples include intense sound beam resulting from the near field high-speed after burner flows, sonic boom, and scattering of high intensity noise caused by the blade and vortex interaction. For example, the interactio n of unsteady disturbanc es with leading and trailing edges of fan and compressor blades has motivated numerous studies (Goldstein, 1976; Wang and Moin, 2000; Ma noha and Herraero, 2002; Ta keda et al., 2003; Howe 1978, 1999, 2000, 2001; Casper and Farassat, 2004; Roger and Moreau, 2005). In this problem, the edge is usually a source of hi gh frequency sound associated with smallerscale boundary layer turbulence. A general description of high fre quency, nonlinear acoustics involves nonequilibrium effects, which arises when e.g., the acoustic time scale becomes comparable to that of the vibrational relaxation process of polyatomic gases (Hamilton and Blackstock, 1997). Another important example of nonlinear be havior is landing gear noise. This problem is recognized as one of the major components of airframe noise for commercial aircraft. The computational investigations should take into consideration the following aspects: the flow conditions with varying Mach number from 0.18 to 0.24 (Guo, 2005), the noise radiation involves complex fl ows around a complex geometry. To be computationally efficient, adaptive grid re finement techniques ar e attractive. Recent efforts reported by Singh et al. (2005) are directly applicable. PAGE 139 125 LIST OF REFERENCES Ashcroft, G. and Zhang, X., 2003, O ptimized Prefactored Compact Scheme, J. Comput. Phys., Vol.190, pp.459-477 Basel, B., and Grnewald, M., 2003,High Orde r Unstructured Finite Difference Method in Aeroacoustics, Presented at CAA: From Acoustic Sources Modeling to FarField Radiated Noise Prediction, EUROMECH Colloquim no. 449, Chamonix, France Bayliss, A., Turkel, E. and Manthey, J ., 1982, Far Field Boundary Conditions for Compressible Flows, J. Comput. Phys., Vol.48, pp.182-199 Blacke, W.K., 1986, Mechanics of Flow-Induc ed Sound and Vvibration, Vols.I and II, Academic Press, New York Blackstock, D.T., 2000, Fundamentals of Physical Acoustic, A Wiley-Interscience publication John Wiley & Sons inc, USA Bogey, C., Bailly, C. and Juve, D., 2003, N oise Investigation of a High Subsonic, Moderate Reynolds Number Jet Using a Cmpressible Large Eddy Simulation, Theor Comput Fluid Dyn, Vol.16, pp.273-297 Bogey, C. and Bailly, C.,2004, A Family of Low Dispersive and Low Dissipative Explicit Schemes for Flow and Noise Computations, J. Comput. Phys. Vol.194, pp.194-214 Calhoun, D. and Le Veque, R.J., 2000, A Cartesian Grid Finite -Volume Method for the Advection-Diffusion Equation in Irregular Geometries, J. Comput. Phys., Vol.157, pp.143-180 Casper, J. and Farassat, F, 2004, Broadband Tr ailing Edge Noise Predictions in the Time Domain, J. Sound Vib, Vol.271, pp.159-176 Chang, S.C., 1995, The method of SpaceTime Conservation Element and Solution Element A New Approach for Solving th e Navier Stokes and Euler Equations, J. Comp. Phys., Vol.119, pp.295-324 Chang, S.C., Wang, X.Y. and Chow, C. Y., 1999, The Space-Time Conservation Element and Solution Element Methods: a New High-Resolution and Genuinely Multidimensional Paradigm for Solving Conservation Laws, J. Comput. Phys., Vol.156, pp.89-136 PAGE 140 126 Cheong, C. and Lee, S., 2001, Grid-Optimi zed Dispersion-Relatio n-Preserving Schemes on General Geometries for Co mputational Aeroacoustics, J. Comput. Phys., Vol.174, pp.248-276 Chu, B.T. and Kovasznay, L.S.G., 1958, N on-Linear Interaction in Viscous HeatConducting Compressible gas, J. Fluid Mech., Vol.3 (5), pp.494-514 Colonius, T. and Lele, S.K., 2004, Computa tional Aeroacoustics: Progress on Nonlinear Problems of Sound Generation, Progress in Aerospace Sciences, Vol.40, pp.345416 Davies, P., Fisher, M.J. and Baratt, M.J., 1963, The Characteristics of Turbulence in the Mixing Region of a Round Jet, J. Comput. Phys., Vol.15, pp.337-367 Delfs, J.W., 2001, An Overlapping Grid Technique for High Resolution CAA Schemes for Complex Geometries, AIAA Paper 2001-2199 Dowling, A.P. and Ffows Williams, J.E., 1983, Sound and Source of Sound, Ellis Horwood Limited, Chicester, England, UK Dumbser, M., 2003, ADER Discontinuous Ga lerkin Schemes for Aeroacoustics, CAA: From Acoustic Sources Modeling to Fa r-Field Radiated Noise Prediction, EUROMECH Colloquim no. 449 Efraimsson, G. and Kreiss, G., 1996, A No te on the Effect of Artificial Viscosity on Solution of Conservation Laws, Applied Numerical Mathematics, Vol.21, pp.155173 Efraimsson, G., 1997, A 2D Analysis of th e Influence of Artific ial Viscosity Terms on Solutions of the Euler Equation, J. Comput. Phys., Vol.138, pp.103-120 Ffowcs-Williams, J.E., 1969, Sound Generation by Turbulence and Surfaces Arbitrary Motion, Phil. Trans. Roy. Soc. A264 Ffowcs-Williams, J.E., 1992, Modern Methods in Analytical Acoustics, D.G.Crighton et al. Springer Verlag Ffowcs-Williams, J.E., and Hawkings, D.L., 1969, Theory Relating to the Noise of Rotating Machinery, J. Sound and Vib. Vol.10(1) Fisher, M.J., Preston, G.A. and Bryce, W.D., 1998, A M odelling of the Noise from Simple Coaxial Jets, Part I: With Unheated Primary Flow, J. Sound Vibr., Vol.209(3), pp.385-403 Gloerfelt, X., Bailly, C., Juv, D., 2003, Direct Computation of the Noise Radiated by A Subsonic Cavity Flow and Application of Integral Methods, J. Sound and Vibration, Vol.266, pp.119-146 PAGE 141 127 Goldstein, M.E., 1976, Aeroacoustics, McGraw-Hill International Book Company, New York Goodrich, J.W., 1995, An Approach to th e Development of Numerical Algorithm for First Order Linear Hyperbolic System in Multiple Space Dimensions: The Constant Coefficient Case, NASA TM 106928 Guo, Y., 2005, A Statistical Model For Landing Fear Noise Prediction, J.Sound and Vibration, Vol.282, pp.6187 Gutin, L., 1948, On the Sound Field of a Rotating Propoller, NACA TM 1195 Hamilton, M.F. and Blackstock, D.T., 1997, Nonlinear Acoustic, Academic Press Limited, London Hardin, J.C., 1996, Introduction to Computational Aero-Acoustic, Applied AeroAcoustic Lecture, Lecture Series 1996-2004, von Karman Institute for Fluid Dynamics, Chauss de Waterloo, 72, Rhode Saint Gense, Belgium Hardin, J. and Hussaini, M.Y., 1992, Computational Aeroacoustics, Springer-Verlag, New York,Berlin Hardin, J.C and Pope, D.S., 1995, Sound Ge neration by Flow over a Two-Dimensional Cavity, AIAA Journal, Vol.33, No.3, pp.407-412 Harten, A., 1989, ENO Scheme s with Subcell Resolution, J. Comput. Phys., Vol.83, pp.148-184 Harten, A., 1983, High Resolution Sche mes for Hyperbolic Conservation law, J. Comput. Phys., Vol.49, pp.357-393 Hataue, I., 2003, On Analogy a nd Dissimilarity of Dependen ce of Stability on Several Parameters in Flow Simulation, J. Computational and Applied Mathematics, Vol.159, pp.45-53 Henshaw, B., 2004, A high-Order Accurate Solver for Maxwells Equations on Overlaping Grids, 7th Symposium on Overset Composite Grids and Solution Technology, Huntington Beach, Ca Hirschberg, A., 1994, Applied Aeroacoustics, Presented at Applied Aero-Acoustic Lecture, Lecture Series 1994-2004, von Karman Institute for Fluid Dynamics, Chauss de Waterloo, 72, Rhode Saint Gense, Belgium Hirschberg, A. and Rienstra, S.W ., 1994, Elements of Aeroacoustics, Applied AeroAcoustic Lecture, Lecture Series 1994-2004, von Karman Institute for Fluid Dynamics, Chauss de Waterloo, 72, Rhode Saint Gense, Belgium PAGE 142 128 Hirsh, C., 1990, Numerical Computation of Internal and External Flows, Vol.II: Computational Methods for Inviscid and Viscous Flows, Wiley Series in Numerical Methods in Engineering, Massachusettes Hixon, R., 1997, Evaluation of High-Accura cy MacCormack Type Scheme Using Benchmark Problems, NASA Contractor Report 202324, ICOMP-97-03 Hixon, R. and Turkel, E., 1998, High-Accu racy Compact MacCormack-Type Schemes for Computational Aeroacoustics, NASA CR 1998 208672 Hong, Y. and Mingde, S.U., 1999, Applicati on and Comparison of Two SGS Models in Large Eddy Simulation of Fr ee Turbulent Jet Flow, Communication in Nonlinear Science & Numerical Simulation, Vol.4(1) Howe, M.S., 1978, A Review of the Theory of Trailing Edge Noise, J. Sound Vib., Vol.62, pp.437-465 Howe, M.S., 1999, Trailing Edge Noise at Low Mach Numbers, J. Sound Vib., Vol.225(2), pp.211-238 Howe, M.S., 2000, Trailing Edge Noise at Lo w Mach Numbers, Part 2: Attached and Separated Edge Flows, J. Sound Vib., Vol.234(5), pp.761-775 Howe, M.S., 2001, On the Hydroacoustics of Trailing Edge Noise with a Detached Flap, J. Sound Vib., Vol.239(4), pp.801-817 Hu, F.Q., Hussaini, M.Y. and Manthey, J.L., 1996 Low Dissipation and Dispersion Runge-Kutta for Computational Acoustics, J.Comput. Phys., Vol.124, pp.177-191 Hu, Z.W., Morfey, C.L. and Sandham, N. D., 2002, Aeroacoustics of Wall-Bounded Turbulent Flows, AIAA J., Vol.40(3), pp.465 -473 Ingram, D.M., Causon, D.M., and Mingham, C. G., 2003, Development in Cartesian Cut Cell Methods, Mathematics and Comput. In Simulation, Vol.61, pp.561-572 Kim,. J.W., and Lee, D.J., 1996, Optimized Compact Finite Difference Schemes with Maximum Resolution, AIAA J., Vol.34, No. 5, pp.887-893 Kim, C., Roe, P.L. and Thomas, J.P., 1997, Accurate Schemes for Advection and Aeroacoustics, Technical Paper 97-2091, AIAA Press, Washington DC Lele, S.K., 1997, Computationa l Aeroacoustics: A Review, AIAA paper 97-0018 Leonard., B.P., 1988, Simple High-Accur acy Resolution Program for Convective Modeling of Discontinuities, Int. J. Numer. Methods in Fluids, Vol.8, pp.12911318 PAGE 143 129 Lighthill, M.J., 1952, On Sound Generated Aerodynamically. I. General Theory, Proc. Roy. Soc. (London), 211A, 1107, pp.564-587 Lighthill, M.J., 1954, On Sound Generated Aerodynamically. II. Turbulence as a Source of Sound, Proc. Roy. Soc. (London), 22A, 1148, pp.1-32 Lin, S.Y. and Chin, Y.S., 1995, Comparis on of Higher Resolution Euler Schemes for Aeroacoustics Computation, AIAA J., Vol.33, pp.237 Loh C.Y., Hultgren, L.S., Chang, S.C. and Jorgenson, P.C.E., 2000, Noise Computation of a Shock-Containing Supersonic Axis ymmetric Jet by the CE/SE Method, AIAA Paper 2000-0475 Loh, C.Y., Hultgren, L.S. and Chang, S. C., 1998, Computing Waves in Compressible Flow Using the Space-Time Conservation Element and Solution Element Method, AIAA Paper 98-0369 Loh, C.Y., Hultgren, L.S., Chang, S.C. a nd Jorgenson, P.C.E., 1999, Vortex Dynamics Simulation in Aeroacoustics by the Sp aceTime Conservation Element and Solution Element Method, AIAA Paper 99-0359 Loh, C.Y., Hultgren, L.S., Wang, X.Y., Ch ang, S.C. and Jorgenson, P.C.E., 2000, Aeroacoustics Computation for Nearly Fully Expanded Supersonic Jets Using the CE/SE Method, NASA/TM 210225 Liou, M.S. and Steffen, C.J.,1993, A New Flux Splitting Scheme, J.Comput. Phys., Vol.107, pp.23-39 Margot, G. and Olsson, P., 1996, Designing an Efficient Soluti on Strategy for Fluid Flows 1. A Stable High Order Finite Difference Scheme and Sharp Shock Resolution for the Euler Equations, J. Comput. Phys., Vol.129, pp.245-262 Manoha, E. and Herraero, C., 2002, Numer ical Prediction of Airfoil Aerodynamic Noise, AIAA Paper 2002-2573 Morris, P.J., Long, L.N., Bangalore, A. and Wang, Q., 1997, A Parallel ThreeDimensional Computational Aeroacoustic s Method Using Nonlinear Disturbance Equations, J.Comput. Phys., Vol.133, pp.56-74 Morse, P.M., and Incard, K.U., 1968, Theoretical Acoustic, McGraw-Hill Book Company, New York Nance, D.V., Viswanathan, K. and Sankar, L.N., 1997, Low-Dispersion Finite Volume Scheme for Aeroacoustic Application, AIAA Journal, Vol.35, No.2, pp.255-262 Oran, E.S. and Boris, J., 2002, Numerical Simulation of Reactive Flow, Cambridge University Press, DC/New York, second edition PAGE 144 130 Osher, S. and Chakravarthy, S., 1983, U pwind Schemes and Boundary Conditions with Applications to Euler Equati ons in General Geometries, J. Comput. Phys., Vol.50, pp.447-481 Popescu, M. and Shyy, W., 2002 Assessmen t of Dispersion-Rela tion-Preserving and Space-Time CE/SE Schemes for Wave Equations, Numerical Heat Transfer, Vol.42, No 2, pp.93-118 Popescu, M., Shyy, W., and Garbey, M ., 2004, A Study of Dispersion-RelationPreserving and Optimized Prefactored Compact Schemes for Wave Equation, AIAA paper 2004-0519, also accepted for publication in J. Comput. Phys. Powell, A., 1964, Theory of Vortex Sound, JASA Vol.33(1), pp.177-195 Roe, P.L., 1981, Numerical Algorith ms for the Linear Wave Equation, Technical Report 81047, Royal Aircraft Establishment, Bedford, UK Roger, M., 1995, Applied Aeroacoustics Me thods for the Acoustic Design of Shrouded Tail Rotor, AIAA Aeroacoustics Conference, Mnschen Roger, M. and Moreau, S., (accepted Oct ober 2004), Back-Scattering Correction and Further Extensions of Amiets Traili ng-Edge Noise Model. Part 1: Theory J.Sound Vibr. Sherer, S., 2004, High-Order OversetGrid Acitivities in AFRL /VAAC, Air Forces Research Laboratory, Presented at 7th Symposium on Overset Composite Grids and Solution Technology, Huntington Beach, California Shu, C. and Oscher, S., 1988, Efficient Impl ementation of Essentially Non-Oscillatory Shock-Capturing Schemes II, J. Comput. Phys., Vol.83, pp.32-78 Shur, M.L., Spalart, P.R., Strelets, M.Kh. and Travin, A.K., 2003, Towards of Noise from Jet Engine, Int.J.Heat and Fluid Flow, Vol.24, pp.551-561 Shyy, W., 1985 A Study of Finite Differe nce Approximations to Steady-State, Convection-Dominated Flow Problems, J. Comput. Phys., Vol.57, pp.415-438 Shyy, W., 1994, revised printing 1997, Computational Modeling For Fluid Flow and Interfacial Transport, Elsevier, Amsterdam, The Netherlands Shyy, W., Francois, M. and Udaykumar, H. S., 2001, Cartesian and Curvilinear Grid Method for Multi-Domain, Moving Boundary Domain Thirteenth International Conference on Domain Decomposition Methods, CIMNE, Barcelona, M. Debit et al. editors Shyy, W., Francois, M., Udaykumar, H.S ., NDri, N, and Tran-Son-Tay, R., 2001, Moving boundaries in microscale biofluid dynamics, Applied Mechanics Reviews, Vol.54, pp.405-453 PAGE 145 131 Shyy, W., Thakur, S.S., Ouyang, H ., Liu, J. and Blosch, E., 1977, Computational Techniques for Complex Transport Phenomena, Cambridge University Press, New York Singh, R.K., NDri, N.N., Uzgoren, E., Shyy, W., and Garbey, M., 2005 ThreeDimensional Adaptive, Cartesian Grid Method for Multiphase Flow Computations, AIAA Paper 2005-1389, 43rd Aerospace Sciences Meeting & Exhibit Stanescu, D. and Habashi, W.G., 1998 N-Storage Low Dissipa tion and Dispersion Runge-Kutta for Computational Acoustics, J. Comput. Phys., Vol.143, 674-681 Sun, M. and Takayama, K., 1999, C onservative Smoothing on an Adaptive Quadrilateral Grid, J. Comput. Phys., Vol.150, pp.143-180 Tam, C. K. W. and Dong, Z., 1994, Wa ll Boundary Condition for High-Order Finite Difference Schemes in Computational Aeroacoustics , Theor. Comput. Fluid Dynam. Vol.6, pp.303 Tam, C. K. W. and Hardin, J. C., 1997 (Eds), Second Computational Aeroacoustics Workshop on Benchmark Problems, NASA CP-3352 Tam, C.K.W. and Burton, D.E., 1984, S ound Generation by Instability Waves of Supersonic Flows. Part2. Axisymetric Jets, J. Fluid. Mech., Vol.138 pp.273-295 Tam, C.K.W. and Morris, P.J., 1980, The Radiation of Sound by Instability Waves of a Compressible Plane Turbulent Shear Layer, J. Fluid. Mech., Vol.98(2), pp.349381 Tam, C.K.W., and Kurbatskii, K.A., 2000, A Wavenumber Based Extrapolation and Interpolation Method for Use in Conjunc tion with High-Order Finite Difference Schemes, J. Comput. Phys., Vol.157, pp.588-617 Tam, C.K.W., and Webb, J.C,., 1994, Radiation Boundary Condition and Anisotropy Correction for Finite Difference Solu tions of the Helmholtz Equation, J. Comput. Phys., Vol.113, pp.122-133 Tam, C.K.W., and Webb, J.C., 1993, Dispe rsion-Relation-Preservi ng Finite Difference Schemes for Computational Acoustics, J. Comput. Phys., Vol.107, pp.262-281 Tam, C.K.W., Webb, J.C., and Dong, Z., 1993, A study of the Short Wave Components in Computational Acoustics, J. Comput. Acoust., Vol.1, pp.1-30 Takeda, K., Zhang, X. and Nelson, P.A., 2003, Computational Aeroacoustic Simulation of Leading-Edge Slat Flow, J.Sound and Vibration, Vol.270, pp.559-572 Tang, L. and Baeder, J., 1999, Uniformly Accurate Finite Difference Schemes for pRefinement, SIAM J. Sci. Comput, Vol.20, No.3, pp.1115-1131 PAGE 146 132 Tannehill, C., Anderson, D.A. and Pletcher, R.H., 1997, Computational Fluid Mechanics and Heat Transfer, Taylor & Francis Group, Philadelphia, PA Tennekes, H. and Lumley, L., 1972 ,A First Course in Turbulence, The MIT Press, Cambridge, Massachusetts, and London, England Thakur, S., Shyy, W. and Liou, M.-S., 1996 a, "Convection Treatment and Pressure Splitting for Sequential Solution Procedur es. Part II: Pressure-Based Algorithm, Numerical Heat Transfer; Part B, Vol.29 pp.29-42 Thakur, S.S., Shyy, W., and Liou, M.S., 1996b, Convection Treatment and Pressure Splitting for Sequential Solution Procedures. Part I: Theory and OneDimensional Test Cases, Numerical Heat Transfer, Part B, Vol.29, pp.1-27 Toro, E.F. and Billet, S.J., 1996 Centred TVD Schemes for Hyperbolic Conservation Laws. Technical Report MMU-9603, Department of Mathematics and Physics, Manchester Metropoli tan University, UK Udaykumar, H.S., Kan, H.C., Shyy, W. and Tran-Son-Tay, R ., 1997 Multiphase Dynamics in Arbitrary Geometri es on Fixed Cartesian Grids J. Comp. Phys., Vol.137, pp.366-405 Udaykumar, H.S., Mittal, R. and Shyy, W ., 1999 Computational of Solid-Liquid Phase Fronts in the Sharp Interface Limit on Fixed Grids, J. Comput. Phys., Vol.153, pp.535-574 van Leer, B., 1979, Towards the Ultimate Conservation Difference Scheme V. A Second Order Sequel to Godunovs Method, J. Comput. Phys., Vol.32, pp.101136 Verstappen, R.W.C.P and Veldman, A.E.P., 2000 Numerical Computation of a Viscous Flow around a Circular Cylin der on a Cartesian Grid, European Congress Methods in Applied Scie nces and Engineering, Barcelona, 11-14 September Vichnevetsky, R., 1987, Wave propagation and Reflection in Irregular Grids for Hyprebolic Equations, Applied Numerical Mathematics, Vol.3, pp.133-166 Wang, G. and Sankar, L.N., 1999, Prediction of Rotorcraft Noise with a Low-Dispersion Finite Volume Scheme, AIAA-99-0480 Wang, X.Y., Chang, S.C. and Jorgenson, P.C.E., 2000, Accuracy Study of the SpaceTime CE/SE Method for Computational Ae roacoustics Problems Involving Shock Waves, AIAA Paper No.2000-0474, 38th Aerospace Sciences Meeting & Exhibit, Reno, NV, January 10-13 Wang, Z.J. and Chen, R.F., 2001, Optimi zed Weighted Essent ially Nonoscillatory Schemes for Linear Waves with Discontinuity, J. Comput. Phys., Vol.174, pp.381-404 PAGE 147 133 Wang, M. and Moin, P., 2000, Computationa of Trailing-Edge Flow and Noise Using Large-Eddy Simulation, AIAA J., Vol.38, pp.2201-2209 Williams, E.G., 1999, Fourier Acoustics Sound Radia tion and Nearfield Acoustical Holography, Academic Press, Dan Diego Yang, G. and Ingram, D.M., 1999 Cartesi an Cut-Cell Method for Axisymmetric Separating Body Flows, AIAA Journal, Vol.37, No.8, pp.905-911 Ye, T., 2001, Direct Numerical Simulation of a Translating Vapor Bubble With Phase Change, Ph.D. Dissertation, University of Fl orida, Department of Mechanical Engineering, Gainesville FL Ye, T., Mittal, R., Udaykumar, H.S. and Shyy, W., 1999 An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries, J. Comput. Phys., Vol.156, pp.209-240 Ye, T., Shyy, W. and Chung,J.N., 2001, A Fixed-Grid, Sharp-Interface Method for Bubble Dynamics and Phase Change, J. Comput. Phys. Vol.174, pp.781-815 Yu, S.T., Chang, S.C., 1997, Treatments of Stiff Source Terms in Conservation Laws by the Methods of SpaceTime Conser vation Element and Solution Element, AIAA Paper No.97-0435, 35th Aerospace Sciences Mee ting & Exhibit, Reno, NV, January PAGE 148 134 BIOGRAPHICAL SKETCH Mihaela Popescu was born in Romania, in a small village, Berca, known for rare natural phenomenamud volcanoe s. Mihaela received two b achelor degrees from the University of Bucharest, Romania. Her firs t BS was in mathema tics, specializing in mathematics applied in mechanics. Her second BS was in bi ology, specializing in ecology. While pursuing her degree in biol ogy, Mihaela worked with modeling of Chironomid community based on ecosystem char acteristics. In adition, for over a year, Mihaela worked for the Nuclear Energy R eactor Institute in Bucharest, Romania. Mihaelas other qualification also entail working as a teacher for two years. In 1999, she enrolled in the University of Florida to pursue a Doctor of Philosophy degree in Aerospace Engineering under the guidance of Dr. Wei Shyy. During her Ph.D. studies she was the reci pient of the Alumni Fellowship at the University of Florida. Mud volcano cones in the Berca Anticline Depression of the Curvature Subcarpathians. [Photograph Credit: Dr Dan Balteanu, Romanian Academy] |