ML052220086

E-mail Gucwa, Entergy Nuclear Northeast, to Ennis, NRR, VY EPU Supplement 30

ML052220086
Person / Time
Site:	Vermont Yankee
Issue date:	08/02/2005
From:	Gucwa L Entergy Nuclear Northeast
To:	Richard Ennis Office of Nuclear Reactor Regulation
References
Download: ML052220086 (50)

Text

I-Rick Ennis - -BVY-0-5--0-5--2--E-xhi-b-it-s--S-R-X-B--A-18,--,-l.---E,-M"-E--B--B,--l--3--8--l-,2-,a-n-d Paq IRick Ennis - BVYO5-052 Exhibits SRXB-A-18-1, EMEB-B-138-1, 2, and 3 Paei 11 From: "Gucwa, Len" <LGUCW90@entergy.com>

To: "Rick Ennis" <RXE@nrc.gov>

Date: 8/2/05 2:23PM

Subject:

BVY 05-052 Exhibits SRXB-A-18-1, EMEB-B-138-1, 2, and 3

<<BVY 05-072 Ex.SRXB-A-18-1.pdf>> <<BVY 05-072 Ex.EMEB-B-138-1.pdf>> <<BVY 05-072 Ex.EMEB-B-138-2.pdf>> <<BVY 05-072 Ex.EMEB-B-138-3.pdf>>

That's everything!

Len T. Gucwa, P.E.

VY Licensing Igucw90@entergy.com 802/451-3193 CC: <Douglas.Rosinski~pillsburylaw.com>

CC:

Ic:\temp\GW)00001 .TMP Page I c:\temp\GWIOOOOI .TMP Page 11 Mail Envelope Properties (42EFB9F0.A8A: 23: 2698)

Subject:

BVY 05-052 Exhibits SRXB-A-18-1, EMEB-B-138-1, 2, and 3 Creation Date: 8/2/05 2:14PM From: "Gucwa, Len" <LGUCW90(entergy.com>

Created By: LGUCW90@entergy.com Recipients nrc.gov owf4po.OWFNDO RXE (Rick Ennis) pillsburylaw.com Douglas.Rosinski CC Post Office Route ovvf4_po.OWFNDO nrc.gov pillsburylaw.com Files Size Date & Time MESSAGE 232 08/02/05 02:14PM TEXT.htm 2431 BVY 05-072 Ex.SRXB-A-18-1.pdf 600742 BVY 05-072 Ex.EMEB-B-138-1.pdf 832848 BVY 05-072 Ex.EMEB-B-138-2.pdf 972271 BVY 05-072 Ex.EMEB-B-138-3.pdf 1256329 Mime.822 5017391 Options Expiration Date: None Priority: Standard Reply Requested: No Return Notification: None Concealed

Subject:

No Security: Standard

BVY 05-072 Docket No. 50-271 Exhibit SRXB-A-18-1 Vermont Yankee Nuclear Power Station Proposed Technical Specification Change No. 263 - Supplement No. 30 Extended Power Uprate Response to Request for Additional Information "RS-001, BWR Template SE" Total number of pages in this Exhibit (excludina this cover sheet) Is 8.

(2)draft GDC-9, insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage

.throughoutits design lifetime; 2.8.5 Accident and Transient Analyses 2.8.5.1 Decrease in Feedwater Temperature, Increase in Feedwater Flow, Increase in Steam Flow, and Inadvertent Opening of a Main Steam Relief or Safety Valve Reaulatorv Evaluation Excessive heat removal causes a decrease In moderator temperature which Increases core reactivity and can lead to a power level increase and a decrease in shutdown margin. Any unplanned power level Increase may result In fuel damage or excessive reactor system pressure. Reactor protection and safety systems are actuated to mitigate the transient. The NRC staff's review covered (1) postulated Initial core and reactor conditions, (2) methods of thermal and hydraulic analyses, (3) the sequence of events, (4) assumed reactions of reactor system components, (5)functional and operational characteristics of the reactor protection system, (6) operator actions, and (7)the results of the transient analyses. The NRC's acceptance criteria are based on (1)draft GDC-6, insofar as It requires that the reactor core be designed to function throughout its design lifetime without exceeding acceptable fuel damage limits, ( draft GDC-14 and 15, Insofar as they require that the core protection system be designed to act automatically to prevent or suppress conditions that could result in exceeding acceptable fuel damage limits and that protection systems be prgvided for sensing accident situations and Initiating the operation of necessary ESFs; and (k) draft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained In SRP Section 15.1.1-4 and other guidance provided in Matrix 8 of RS-001.

Technical Evaluation

[insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements In the regulatory evaluation and (2)provide a clear link to the conclusions reached by the NRC staff, as documented In the conclusion section.)

Conclusion The NRC staff has reviewed the licensee's analyses of the excess heat removal events a described above and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The'NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that theA sand the RCPB pressure limits will not be exceeded as a result of these events. Basedon this, the NRC sta q concludes that the plant will continue to meet the requirements of draft GDC6,14, 15, 27, and 28 following implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the events stated.

INSERT 8 FOR SECTION 32 . BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

(2) draft GDC-9, insofar as it required that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2..5.2 Decrease in Heat Removal by the Secondary System 2.8.5.2 .1 Loss of External Load; Turbine Trip; Loss of Condenser Vacuum; Closure of Main Steam Isolation Valve; and Steam Pressure Regulator Fallure (Closed)

ReQulatorv Evaluation A number of initiating events may result in unplanned decreases in heat removal by the secondary system. These events result in a sudden reduction In steam flow and, consequently, result In pressurization events. Reactor protection and safety systems are actuated to mitigate the transient. The NRC staffs review covered the sequence of events, the analytical models used for analyses, the values of parameters used In the analytical models, and the results of the transient analyses. The NRC's acceptance criteria are based on (1) draft GDC-6, Insofar as it requires that the reactor core be designed to functon throughout Its design lifetime without exceeding acceptable tuel damage lirnlt and (urdraft GDC-27 and 28, insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained in SRP Section 15.2.1-5 and other guidance provided in Matrix 8 of RS-001.

.*Technical Evaluation

[Insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements in the regulatory evaluation and (2) provide a clear link to the conclusions reached by the NRC staff, as documented in the conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the'decrease in heat removal events described above and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that theF~g)and the RCPB pressure limits will not be exceeded as a result of these events. Baseeon this, the NRC staff q) concludes that the plant will continue to meet the requirements of draft GDC-6,27, and 28 following implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the events stated.

INSERT8 FOR SECTION 3.2-BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

(2) draft GDC-9, insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2.8.5.2.2 Loss of Nonemergency AC Power to the Station Auxiliaries Regulatory Evaluation The loss of nonemergency ac power is assumed to result in the loss of all power to the station auxiliaries and the simultaneous tripping of all reactor coolant circulation pumps. This causes a flow coastdown as well as a decrease In heat removal by the secondary system, a turbine trip, an increase In pressure and temperature of the coolant, and a reactor trip. Reactor protection and safety systems are actuated to mitigate the transient. The NRC staffs review covered (1) the sequence of events, (2) the analytical model used for analyses, (3) the values of parameters used in the analytical model, and (4) the results of the transient analyses. The NR.C's acceptance criteria are based on (1) draft GDC-6, insofar as It requires that the reactor ed to function throughout Its design lifetime without exceeding acceptable fuel damage limits, nd (tpdraft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained In SRP Section 15.2.6 and otherguidance provided in Matrix 8 of RS-001.

Technical Evaluation

[Insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements Inthe regulatory evaluation and

• (2) provide a clear link to the conclusions reached by the NRC staff, as documented In the conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the loss of nonemergency ac power to station auxiliaries event and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee bggmderonstrated that the reactor protection and safety systems will continue to ensure that the bnd the RCPB pressure limits will not be exceeded as a result of this event. Based on thi the NR staf concludes that the plant will continue to meet the requirements of draft GDC-6,, and 28 following Implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the loss of nonemergency ac power to station auxiliaries event INSERT 8 FOR SECTION 3.2- BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

(2)draft GDC-9, insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2.8.5.2.3 Loss of Normal Feedwater Flow Regulatorv Evaluation A loss of normal feedwater flow could occur from pump failures3 valve malfunctions, or a LOOP.

Loss of feedwater flow results Inan increase Inreactor coolant temperature and pressure which eventually requires a reactor trip to prevent fuel damage. Decay heat must be transferred from fuel following a loss of normal feedwater flow. Reactor protection and safety systems are actuated to provide this function and mitigate other aspects of the transient. The NRC staff's review covered (1)the sequence of events, (2)the analytical model used for analyses, (3)the values of parameters used Inthe analytical model, and (4)the results of the transient analyses.

The NRC's acceptance criteria are based on (1)draft GDC-6, Insofar as it requires that the reactor core be designed to function throughout its design lifetime without exceeding acceptable fuel damage limits, and td(raft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable

.fuel damage limits. Specific review criteria are contained in SRP Section 15.2.7 and other guidance provided InMatrix 8 of RS-001.

Technical Evaluation

[nsert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements In the regulatory evaluation and (2) provide a clear link to the conclusions reached by the NRC staff, as documented Inthe conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the loss of normal feedwater flow event and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that th4X7tF9 and the RCPB pressure limits will not be exceeded as a result of the loss of normatfeUEwater flow. Based on this, the NRC staff concludes that the plant will continue to meet the requirements of draft GDC6, 7, and 28 e following implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the loss of normal feedwater flow event.

INSERT 8 FOR SECTION 3.2 - BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

(2)draft GDC-9, Insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2.8.5.3 Decrease In Reactor Coolant System Flow 2.8.5.3.1 Loss of Forced Reactor Coolant Flow Regulatory Evaluation A decrease in reactor coolant flow occurring while the plant Is at power could result In a degradation of core heat transfer. An increase in fuel temperature and accompanying fuel damage could then result iA Ilare exceeded during the transient. Reactor protection and safety systems are actuated to mitigate the transient. The NRC staffs review covered (1) the postulated initial core and reactor conditions, (2)the methods of thermal and hydraulic analyses, (3) the sequence of events, (4) assumed reactions of reactor systems components, (5)the functional and operational characteristics of the reactor protection system, (6)operator actions, and (7)the results of the transient analyses. The NRC's acceptance criteria are based on (1) draft GDC-6, insofar as it requires that the reactor core be designed to function throughout its design lifetime without exceeding acceptable fuel damage limitswand (kdraft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained In SRP Section 15.3.1-2 and other guidance provided In Matrix 8 of RS-001.

Technical Evaluation

[Insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements Inthe regulatory evaluation and (2) provide a clear link to the conclusions reached by the NRC staff, as documented In the conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the decrease in reactor coolant flow event and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models.

The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that the(A L)and the RCPB pressure limits will not be exceeded as a result of this event. Based on this, the NRC staff concludes that the plant will continue to meet the requirements of draft GDC-6P7, and 28 following implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the decrease in reactor coolant flow event.

INSERT 8 FORSECTION 32 -BWRTEMPLATE SAFETY EVALUATION DECEMBER 2003

I (2)draft GDC-9, insofar as It requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2;8.5.5 Inadvertent Operation of ECCS or Malfunction that Increases Reactor Coolant Inventory Regulatory Evaluation Equipment malfunctions, operator errors, and abnormal occurrences could cause unplanned Increases in reactor coolant Inventory. Depending Qn the temperature of the Injected water and the response of the automatic control systems, a power level increase may result and, without adequate controls, could lead to fuel damage or ovdrpressurization of the RCS. Alternatively, a power level decrease and depressurization may result. Reactor protection and safety systems are actuated to mitigate these events. The NRC staffs review covered (1)the sequence of events, (2) the analytical model used for analyses, (3)the values of parameters used In the analytical model, and (4)the results of the transie'nt analyses. The NRC's acceptance criteria are based on (1) draft GDC6, insofar as it requires hat the reactor core bi designed t9 function throughout its design Ii etimewithpot exce eing acceptable fuel damage limitsand (d)draft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained in SRP Section 15.5.1-2 and other guidance provided in Matrix 8 of RS-001.

Technical Evaluation

[Insert technical evaluation. The technical evaluation should (1)clearly explain why the proposed changes satisfy each of the requirements In the regulatory evaluation and (2)provide a clear link to the conclusions reached by the NRC staff, as documented in the conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the inadvertent operation of ECCS or malfunction that increases reactor coolant Inventory and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that the AFDgbnd the RCPB pressure limits will not be exceeded as a result of this event.

Based on t is,the NRC staff concludes that the plant will continue to meet the requirements of draft GDC are, the NRC staff finds the proposed EPU acceptable with respect to the inadvertent operation of ECCS or malfunction that increases reactor coolant inventory.

INSERT8 FOR SECTION 3.2- BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

(2) draft GDC-9, Insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an exceedingly low probability of gross rupture or significant leakage throughout its design lifetime; 2.8.5.4.3 Startup of a Recirculation Loop at an Incorrect Temperature and Flow Controller Malfunction Causing an Increase in Core Flow Rate Regulatorv Evaluation A startup of an Inactive loop transient may result Ineither an Increased core flow or the introduction of cooler water into the sore. This event causes an Increase In core reactivity due to decreased moderator temperature and gore void fraction. The NRC staffs review covered (1) the sequence of events, (2)the analytical model, (3)the values of parameters used In the analytical model, and (4)the results of the transient analyses. The NRC's acceptance criteria are based on (1) draft GDC-6, Insofar as it requires that the reactor core be designed to function throughout its design lifetime Without eceeding acceplda el (raft GDC-14 and 15, insofar as they require that the pore protection systems be designed to act automatically to prevent or suppress conditions that could result In exceeding acceptable fuel damage limits and that protection systems be provided for sensing accident situations and Initiating the operation of necessary ESFs; (pdraft GPC-32, Insofar as It requires that limits, which Include considerable margin, be placed on the maximum reactivity worth of control rods or elements and on rates at which reactivity can be Increased to ensure that the potential effects of a sudden or large change of reactivity cannot (a)rupture the reactor coolant pressure boundary or (b) disrupt the core, its support structures, or other vessel internals sufficiently to Impair the effectiveness of emergency core cooling; and (X)<draft GDC-27 and 28, Insofar as they require that at least two reactivity control systems be provided and be capable of making and holding tMe core subcritical from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained InSRP Section 15.4.4-5 and other guidance provided In Matrix 8 of RS-001.

Technical Evaluation

[insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements In the regulatory evaluation and (2) provide a clear link to the conclusions reached by the NRC staff, as documented in the conclusion section.]

Conclusion The NRC staff has reviewed the lcensee's analyses of the increase Incore flow event and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that theT Ds)and the RCPB pressure limits will not be exceeded as a result of this event. Based on this fe NRC staff concludes that the planwl continue to meet the requirements of draft GDC ti4, 15,27,28, and 2 following implementation of the proposed EPU. Therefore, the NRC staff finds the proposed EPU acceptable with respect to the increase in core flow event.

INSERT 8 FOR SECTION 3.2 -1BWRTEMPLATE SAFETY EVALUATION DECEMBER2003

1 (2) draft GDC-9, insofar as it requires that the reactor coolant pressure boundary shall be designed and constructed so as to have an

-exceedingly low probability of gross rupture or significant leakage throughout its design lifetime:

I .

.8.5.6 Decrease In Reactor Coolant Inventory 2.8.5.6.1 Inadvertent Opening of a Pressure Relief Valve Reciulatorv Evaluation The Inadvertent opening of a pressure relief valve results in a reactor coolant Inventory decrease and a decrease in RCS pressure. The pressure relief valve discharges into the suppression pool. Normally there Is no reactor trip. The pressure regulator senses the RCS pressure decrease and partially closes the turbine control valves (TCVs) to stabilize the reactor at a lower pressure. The reactor power settles out at nearly the Initial power level. The coolant inventory is maintained by the feedwater control system using water from the condensate storage tank via the condenser hotwell. The NRC staff's review covered (1) the sequence of events, (2) the analytical model used for analyses, (3) the values of parameters used in the analytical model, and (4) the results of the transient analyses. The NRC's acceptance criteria are based on (1) dft GDC6, insofa as t uires that the re be desi ned to function throughout its design etime without exceeding acceptable fuel damage limitsand raft GDC-27 and 28, insofar as they require that at least two reactivity control systems be provided and be capable of making and holding the core subcritcal from any hot standby or hot operating condition sufficiently fast to prevent exceeding acceptable fuel damage limits. Specific review criteria are contained in SRP Section 15.6.1 and other guidance provided in Matrix 8 of RS-001.

Technical Evaluation

[Insert technical evaluation. The technical evaluation should (1) clearly explain why the proposed changes satisfy each of the requirements In the regulatory evaluation and (2) provide a clear link to the conclusions reached by the NRC staff, as documented In the conclusion section.]

Conclusion The NRC staff has reviewed the licensee's analyses of the inadvertent opening of a pressure relief valve event and concludes that the licensee's analyses have adequately accounted for operation of the plant at the proposed power level and were performed using acceptable analytical models. The NRC staff further concludes that the licensee has demonstrated that the reactor protection and safety systems will continue to ensure that theDL)and the RCPB L.

pressure limits will not be exceeded as a result of this event. Based on tHishe NRC staff concludes that the plant will continue to meet the requirements. of draft GD and 28 following implementation of the proposed EPU. Therefore, the'NRC staff finds the proposed EPU acceptable with respect to the Inadvertent opening of a pressure relief valve event.

INSERT 8 FOR SECTION 32 - BWR TEMPLATE SAFETY EVALUATION DECEMBER 2003

BVY 05-072 Docket No. 50-271 Exhibit EMEB-B-138-1 Vermont Yankee Nuclear Power Station Proposed Technical Specification Change No. 263 - Supplement No. 30 Extended Power Uprate Response to Request for Additional Information Prediction of Unsteady Loading on A Circular Cylinder in High Reynolds Number Flows using Eddy Simulation.

l Total number of pages in this Exhibit (excludina this cover sheet) is 9.

Proceedings of OMAE 2005:

24 th International Conference on Offshore Mechanics and Arctic Engineering June 12 - 16, 2005, Halkidiki, Greece OMAE 2005-67044 PREDICTION OF UNSTEADY LOADING ON A CIRCULAR CYLINDER IN HIGH REYNOLDS NUMBER FLOWS USING LARGE EDDY SIMULATION Sung-Eun Kim L. Srinivasa Mohan Fluent Incorporated, Lebanon, N.H. Fluent India 03766, U.S.A. Pune, India ABSTRACT of large-scale turbulent structure around a circular cylinder, which in turn affects unsteady loading on the cylinder.

Large eddy simulations were carried out for the flow Apparently, turbulence modeling is an issue here.

around a hydrodynamically smooth, fixed circular cylinder at There are largely three approaches being explored by CFD two Reynolds numbers, one at a subcritical Reynolds number practitioners to modeling high-Re turbulent flows around (Re = 1.4 x 105) and the other at a supercritical Reynolds circular cylinders and bluff bodies in general. One approach number (Re= 1.0 x 106). The computations were made using a employs unsteady Reynolds-averaged Navier-Stokes (URANS) parallelized finite-volume Navier-Stokes solver based on a equations supplemented by turbulence models as the governing multidimensional linear reconstruction scheme that allows use equations. URANS-based approach, despite its low of unstructured meshes. Central differencing was used for computational cost mainly due to less demanding mesh discretization of both convection and diffusion terms. Time- resolution requirement, is fundamentally ill equipped to capture advancement scheme, based on an implicit, non-iterative large-scale turbulent structure present in the flows. Results fractional-step method, was adopted in conjunction with a obtained using URANS computations typically underpredict three-level, backward second-order temporal discretization. the amplitudes of fluctuating forces. Obviously, this Subgrid-scale turbulent viscosity was modeled by a dynamic shortcoming has a serious negative implication to accurate Smagorinsky model adapted to arbitrary unstructured meshes prediction of VIV.

with the aid of a test-filter applicable to arbitrary unstructured Fundamentally more viable than URANS-based approach meshes. The present LES results closely reproduced the flow for bluff-body flows is large eddy simulation (LES). LES is features observed in experiments at both Reynolds numbers. designed to directly resolve large eddies, with the effects of The time-averaged mean drag coefficient, root-mean-square unresolved smaller eddies, namely subgrid-scale turbulence, on force coefficients and the frequency content of fluctuating the resolved large eddies accounted for by turbulence models.

forces (vortex-shedding frequency) are predicted with a However, LES is computationally very expensive, often commendable accuracy. prohibitively, especially when thin turbulent boundary layer is a predominant feature of the flow in question to be accurately INTRODUCTION resolved in a LES. Resolving near-wall turbulence with infinitesimal length and time scales requires extremely fine Unsteady loading on a circular cylinder caused by the flow mesh and small time-step size. For that matter, LES is more is a precursor to its vortex-induced vibration (VIV). The major feasible for free flows such as jets, mixing layer, and wakes, difficulty of computing flows around circular cylinders and massively separated flows.

originate from the fact that circular cylinder flows of practical There are very few LES studies published in the literature interest are high Reynolds number (Re) flows, often involving a of flows around circular cylinders at high Reynolds numbers.

laminar-to-turbulent transition in various regions of the flows Breuer [I] was the one of the very few who tackled the case of such as boundary layer, separated shear-layer, or near-wake. a high subcritical Reynolds number, Re = 1.4 x 105, at which The location of transition and the extent of laminar or turbulent the experiment of Cantwell and Coles [2] was conducted.

flow regime, in real flows, depend on such factors as Reynolds Using a multi-block, structured-mesh-based finite-volume number, freestream turbulence, surface roughness, span- solver with an explicit time-marching scheme, Breuer was able diameter ratio (LD), and blockage, among others. The state of to predict the global parameters of the flow, and the mean flow the flow in those regions dictates the formiation and evolution and turbulence with a commendable accuracy. More recently, I Copyright © 2005 by ASME

Catalano et al. 13] attempted LES for three Reynolds numbers triangular, tetrahedral, pyramidal, prismatic, and hybrid in critical to supercritical Reynolds regime (Re = 0.5 x 106, 1.0 meshes. The solution gradients at cell centers that are needed x 106, 2.0 x 106). Their predictions of the global flow to compute convective and diffusive fluxes are determined parameters were in reasonable agreement with the experimental applying the Green-Gauss theorem [7]. Diffusive fluxes are data, although the skin-friction and the Reynolds number- discretized using central differencing. Discretization of dependency of the mean drag coefficient were poorly predicted. convective fluxes requires great caution in LES. Upwind-A more recent trend in turbulence modeling of bluff-body biased schemes such as second-order upwind scheme have been flows is to employ what may be called "hybrid" approaches in widely used for RANS computations. Unfortunately, numerical which one attempts to adjust turbulence models to local mesh diffusion introduced by upwind schemes, which might be resolution in one way or another. In one such approach acceptable in RANS computations for high-Re flows where referred to, in the literature, as detached eddy simulation eddy-viscosity is larger than molecular viscosity by orders of (DES), either a RANS-based or a subgrid-scale (SGS) magnitude, can overwhelm physical diffusion that is typically turbulence model is invoked depending on whether or not the much smaller in LES; subgrid-scale turbulent viscosity is filter-length (local mesh size) is larger than the local integral smaller than RANS-based eddy viscosity by orders of length-scale of turbulence. In DES, turbulence models magnitude. For this reason, a central differencing scheme [13]

essentially reduce to RANS models in near-wall region or when was employed for the discretization of convection terms in this the local mesh size is too coarse to explicitly resolve energy- study.

containing eddies. One fundamental criticism about DES comes down to the lingering question of how one can possibly TIM E-ADVANCEIM ENT SCIIEMIES reconcile a RANS turbulence model and a subgrid-scale turbulence model, two very different models in terms of their An implicit fractional-step method (FSM) [8] in spectral content, at the common interface. It should also be combination with a second-order accurate, three-level realized that, as a consequence of falling back to a RANS backward-differencing scheme for time-discretization was model in near-wall region, the fidelity of DES would be solely employed to advance the solution in time. In this algorithm, the determined by that of the embedded RANS model. Quite a few momentum equations are decoupled from the continuity studies employing DES have appeared recently on circular equation using an approximate factorization of the coupled cylinder flows. Among others, Travin et al.[4], Vatsa and Navier-Stokes equations. For incompressible flows, the FSM Singer[5], and more recently Pandya et al. [6] all employed preserves the formal second-order temporal accuracy without DES based on an eddy-viscosity transport model to study the having to perform, at each time-step, costly outer iterations to flow around a smooth circular cylinder at subcritical and couple velocity and pressure. The FSM thus provides a highly supercritical Reynolds numbers. The results of these DES are efficient algorithm for CPU-intensive transient computations mixed, insofar as some of the global flow parameters such as like LES.

the r.m.s. force coefficients and the mean flow in the near-wake Consider a semi-discrete forn of the Navier-Stokes were predicted poorly. equations in "pressure-correction" (0.41 = p"' - p") form, The objective of this study is to assess the feasibility of LES for high Reynolds number flows around a fixed, smooth A G](up' )=Or) circular cylinder. To that end, two Reynolds numbers, one in subcritical (Re = 1.4 x 105) and the other in supercritical (Re =

1.0 X 106) regime were deliberately chosen so that they straddle

[ D 0j 47) t (I) the critical Reynolds number (- 3 x 105). The main concern is where u"'} and r are the velocity vector and the momentum how well LES can reproduce the known characteristics of the source vector, respectively, and the integer n the time level, flow at the two Reynolds numbers and the changes in the global with n+J refers to the current time level. A is the coefficient flow parameters such as drag coefficient, r.m.s. force matrix defined in terms of the discrete convective and diffusive coefficients, and vortex-shedding frequency as the Reynolds operators and the time-advancement scheme chosen, and G and number increases. D the discrete gradient and divergence operators, respectively.

The results of the LES will be compared whenever The coupled system of equations in Equation (I) is possible, with the experimental data and other LES and DES extremely difficult to solve as it stands, since the matrix A has results reported in the literature. to be inverted for every iteration. In the FSM, the original coupled equations in Equation (I) are approximated by NUMERICAL METHOD The computations were carried out using the segregated A 10 1 == r) [~o(A,2)] (2) solver in FLUENT, a general-purpose CFD software. The ID A0DG][iO I] to (r'nj 2 details of the numerical method are described in the following.

Factorizing equation (2), we have a series of split operations SPATIAL DISCRETIZATION SCIIEMIES that consist of FLUENT employs a cell-centered finite-volume method A6 =-Gp' +r" based on a multidimensional linear reconstruction scheme that AtDG4V = Dfi (3) permits use of computational elements (cells) with arbitrary polyhedral topology, including quadrilateral, hexahedral, u = u -AtG47"'

2 Copyright C 2005 by ASME

On a per-iteration basis, the series of operations in (5,, _ =-C.Y _ (5)

Equation (3) closely resemble the SIMPLEC scheme. The difference is that, in the iterative SIMPLEC scheme, the series of operations in Equation (3) are repeated in a loop until the We determine the model constant, Cx, using the dynamic Smagorinsky model (DSM) originally proposed by Germano et solutions converge, while the FSM needs only one sweep. To preserve second-order accuracy with the FSM, however, sub- al. [12]. To implement the dynamic procedure for the present iterations are needed for the set of three momentum equations finite-volume solver requires a special test-filter applicable to and individual scalar equations to account for the nonlinearity arbitrary unstructured meshes. The dynamic procedure employed in and coupling among the individual equations and high-order in the present study is "local" in the sense that it does not require source terms. Yet, the non-iterative FSM is takes much less the existence of any statistically homogeneous directions, being CPU time than the iterative SIMPLEC scheme. applicable to arbitrary complex three-dimensional flows. The The system of discretized governing equations are solved details of the implementation of the DSM in the framework of using a point-implicit, Gauss-Seidel relaxation along with an unstructured mesh based finite-volume solver can be found in algebraic multigrid (AMG) method to accelerate solution Kim [13]. The DSM has been validated for a number of wall-convergence. The solver and the subgrid-scale turbulence bounded Ilows such as a square cylinder and a sphere.

model are fully parallelized.

As a validation for the spatial discretization schemes and the transient algorithm, we computed laminar flow around a SOLUTION DOMAIN AND COMPUTATIONAL MESH circular cylinder at the Reynolds number of 100 with both non-iterative FSM and iterative SIMPLEC scheme. At this Our choices of the solution domain and the computational Reynolds number, the flow exhibits a periodic vortex-shedding. mesh were guided by the findings in Breuer's LES study [1].

A C-type structured mesh with 48,000 cells was used for the What one should keep in mind in determining the size of the solution domain is that the spanwise extent of the domain computation. A dimensionless time-step size of At = 0.04 (non- should be larger than the spanwise correlation length of dimensionalized by D/U& D being the cylinder diameter, U0 the turbulence. Fortunaiely, the spanwise correlation length is freestream velocity) was used, with which one period of the known to decrease as the Reynolds number increases. We took vortex-shedding was resolved with approximately 150 time-advantage of this fact, having decided to use a spanwise length steps. The predicted mean drag coefficient (ED), r.m.s. lift of 2.OD for both the subcritical and the critical Re numbers.

coefficient (ci), and Strouhal number (St) are summarized in The top and bottom boundaries of the domain were placed at Table 1. The FSM gives practically the same predictions as the 10.5D from the cylinder axis. Thus, the blockage ratio of our iterative SIMPLEC scheme, and the global parameters numerical model (HID, where H is the height of the domain) is predicted by both methods are seen to closely match the other approximately 4.8 %. The upstream inlet and the downstream predictions and the experimental data. exit boundaries are located at 8.5D upstream and 20.5D downstream of the cylinder axis, respectively.

A locally refined, hexahedral mesh with a total of 6.8 Table 1. Summary or the global flow parameters predicted million cells was used for the computations for both Reynolds for laminar flow (Re = 100) - Norberg 110] compiled numbers. The overall mesh resolution is comparable to case numerical results whose references can be found in ref. 110] -3BI" in Breuer [1]. The near-wall mesh resolution is such that the distance from the cylinder surface is 104D at the wall-

_ _ __ _ _ _ C E.__ C i. St adjacent cells, which translates to y+(=ury/l well under 1.0 for Present - FSM 1.336 0.28 0.165 the subcritical flow, and below 6.0 for the supercritical flow.

Present - SIMPLEC 1.345 0.28 0.166 For the subcritical flow, the boundary layer was resolved with 10 - 18 cells in the laminar region.

Park et al. 19] 1.33 0.33 0.165 Norberg I10] - 0.23/0.29 - BOUNDARY CONDITIONS AND OTHER DETAILS OF Williamson [ ] - - 0.164 COMPUTATION On the cylinder wall, we employ a law-of-the-wall that SUBGRID-SCALE TURBULENCE MODELING invokes proper wall-laws depending on the mesh resolution, namely, ye at the wall adjacent cells. The top and bottom wall For incompressible flows, the filtered Navier-Stokes boundaries were treated as a slip (zero-stress) wall. A equations can be WTitten as periodicity condition was applied on the pair of lateral boundaries in the span-wise direction. At the upstream inlet boundary, the freestream condition was specified. At the au~aiu=_I a ha-ri+ _ a IVawi 4 downstream exit boundary, we extrapolated the solution

a. axi pax, ax, axeJ axJ variables from the adjacent interior cells in a mass-conserving manner.

where r., = uu, - jF1, is the subgrid-scale turbulent stress. In The time-step size (dz) used for the computations is 0.005 this study, the subgrid-scale (SGS) stress is modeled using in a dimensionless unit for both Reynolds numbers. The time-isotropic eddy-viscosity as step size was determined based on the estimate of the 3 Copyright C 2005 by ASME

characteristic length and time scales of the smallest resolved eddies, r=1u', where t was taken as 0.05D, and u as Table 2. Summary of the global now pararneters predicted by 0.214 From these rough estimates and with the dimensionless the present LES for Re = 1.4 x 10, compared wilti other time-step size of 0.005, one turnover time of the smallest numerical results (Breuer [1]; DES-TS - Tram'in et td.14]; DES-LS

- Travin et al.[4]; and experimental data (CC - Caniwell &

resolved eddies will be resolved with 50 time steps. Colesl2], WVA - West & Apelt 114], SB: Szepessy & Bearmin [I15],

To shorten the initial transient period of the solution and to NO - Norberg 110], ZDM Zdravkovich [16])

quickly attain a statistically stationary state, we took, as the initial condition, a partially converged steady RANS solution lCD Ca -C l S i with pseudo-random fluctuating velocity components superimposed on the mean velocity field. LES (Present) 1.13 0.59 1.20 0.67 0.205 The computations were carried out on a 16-CPU lnx86 LES (Breuer) 1.45 - 1.76 0.34 0.204 cluster with AMD Opteron processors and Infiniband DES-TS 0.59 0.06 0.67 1.2 0.31 interconnects. With the non-iterative FSM, the computation DES-LS 1.08 0.29 1.04 1.1 0.21 took approximately 12 CPU-seconds per a time-step, which Exp.(CC) 1.24 - 1.21 0.5 0.18 translates into being able to complete roughly 7,200 time-steps Exp.(WA) 1.3 0.58 - - -

in 24 hours2.777778e-4 days <br />0.00667 hours <br />3.968254e-5 weeks <br />9.132e-6 months <br />.

Exp.(SB) 1.35 0.5 1.05 -

Exp.(NO) - 0.52 - 0.18/0.195 RESULTS Exp. (ZD) 1.15 0.5/0.6 1.05/1.2 - 0.18/0.21 Subcritical Reynolds number (Re = 1.4 x105 )

Figure I depicts the flow structure around the cylinder. The present LES gave the mean CD around 1.13, which Although the near-wake flow appears chaotic, the figure clearly falls within the scatter of the experimental data (1.1 - 1.35) for shows the existence of a large-scale motion, namely alternate smooth cylinders. Interestingly, our prediction came out vortex-shedding known to occur in the subcritical Re number substantially lower than Breuer's prediction despite the range. The color map of the resolved turbulent kinetic energy comparable mesh resolution, the same formal order of spatial in Figure 1 also indicates that the boundary layer remains accuracy (second-order) and the same dynamic SGS turbulence laminar up to the separation point, and transition occurs in the model used in both computations.

separated shear layer.

Figure 1. The flow structure behind the circular cylinder at 400 450 500 Re = 1.4 x ]05 - iso-surface of the second invariant of Non-dimensional time, tg (= tUJ/D) velocity deformation tensor, colored by the resolved turbulent kinetic energy Figure 2. Time histories of the drag and lift coefficients predicted by the present LES for Re = 1.4 x 105 The global flow parameters predicted by the present LES are summarized in Table 2 along with the LES prediction of The DES-TS result [4] is shown to severely underpredict Breuer [1] and the DES predictions of Travin et al. [4], and the the mean CD, better matching the data obtained with rough experimental data. As for Breuer's LES result, among the cylinders 1171, which is consistent with the DES predictions by cases involving different mesh resolution and SGS turbulence others [5, 6]. It has been found that DES, run in a normal models, we picked the result of case "Bl" whose mesh and mode with finite freestream turbulence, essentially yields a SGS turbulence model (dynamic Smagorinsky mode) closely tripped boundary layer, leading to a premature transition to match ours. In selecting the experimental data for the turbulence and an early drag crisis even at a subcritical comparison, we gave more weight to the ones measured on Reynolds number. The DES-LS result [4] based on the so-smooth cylinder surface, and with sufficiently large span-wise called "tripless" approach, in which laminar flow is essentially lengths (LID > 5) and smaller blockage-ratio (HID < 0. 1). enforced in the boundary layer, better predicts the mean CD.

4 Copyright © 2005 by ASME

10 I06 02 1* .I I I I I! I1 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 x/D St Figure 3. Power spectral density of CL time-history for Re Figure 5. Time-averaged axial velocity profile in the wake

1.4 x10 predicted by the present LES for Re = 1.4 X 11)5 The r.m.s. lift coefficient predicted by the present LES is Figure 5 depicts the predicted time-averaged axial velocity within the scatter of the experimental data, although it is closer profile along the wake centerline at the mid-span. The overall to the upper bound. The CL history Figure 2 clearly shows agreement between the prediction and the measurement is quite that, despite modulation of the amplitude, there is a distinct good. The rapid relaxation of the mean axial velocity in x/D

frequency of vortex-shedding, which is consistent with what 1.0 - 3.0 shown in the measurement is an indication of has been known for smooth cylinders in the subcritical regime. vigorous mixing of momentum taking place in the near-wake The Strouhal number corresponding to the main shedding region. That our prediction closely captures the recovery of the frequency was found to be 0.205 as shown in Figure 3. mean axial velocity verifies that the dynamics of the large-scale structure in the near-wake is well predicted in the present LES.

2.0 The length of the time-averaged recirculation bubble (see also Table 2 for L,) from our LES was found to be slightly larger than what Cantwell and Coles [2] reported. Yet the present

-Present LES prediction is much closer to the data than the DES 1.0 o Exp. (Cantwell & Coles.1983) predictions. Considering that hot-wire measurements become

• Exp. (Sarioglu & Yavuz. 2002) increasingly difficult and less reliable near the recirculation region, our LES prediction is considered quite good.

i3 0.0 -

.1.0 S 5.0 _U Lower surface

-20

.0 30 60 90 120 150 180 0 (deg.)

Figure 4. Mean static pressure (Cp) distribution on the cylinder surface for Re = 1.4 X lO1 The time-averaged pressure distribution on the cylinder

-5.0 surface is shown in Figure 4, along with the experimental data 12, 18] obtained at the same Reynolds number. The negative . I . . I . I . . I . . . .

peak predicted by the present LES is larger than the measured 0 30 60 90 120 150 ISO 0 (deg.)

ones. However, the pressure level in the separated region - the plateau after 9= 90 degrees, and the mean base pressure (CPB) Figure 6. Scaled instantaneous skin-friction coefficient are all accurately predicted by the present LES. distribution on the cylinder surface at Re = 1.4 x 105 5s Copyright © 2005 by ASME

peaks at 3.0 around 0= 50 degrees, which is closely reproduced 0.00 i IE . I i . . I . . I . .

by the present LES.

O.OW~t l--Upper surrace 0.004e l_% Lower surfacel Supercritical Reynolds number Case (Re = 1.0 x10')

0.003 0.002 .

0.001, I . I . .

0 30 60 90 120 ISO 180 0 kdeg.)

Figure 7. Resolved turbulent kinetic energy distribution at 0.0001D above the cylinder surface, predicted the present LES at Re = 1.4 x 105 Figure 8. The flow structure behind the circular cylinder at Re = 1.0 X 106 - an iso-surface of the second invariant of velocity deformation tensor, colored by the resolved As mentioned earlier, the experimental evidence [19] turbulent kinetic energy indicates that, in the subcritical regime, the boundary layer remains laminar until it finally separates from the cylinder surface. Therefore, it would be of much interest to see how Figure 8 portrays an impression of the instantaneous flow closely LES can reproduce the physics associated with the structure. Compared to the flow structure for the subcritical separation and the transition. Figures 6 and 7 depict the flow shown in Figure 1, the wake is much narrower and more circumferential variations of two quantities along the cylinder chaotic without any trace of a large-scale alternate vortex surface that shed some light on how the LES predicts the shedding, suggesting a delayed turbulent separation. This is separation and transition. typical of the flow in the supercritical regime. The color map Figure 6 shows a plot of the instantaneous skin-friction on the iso-suraface - representing resolved turbulent kinetic energy - also indicates that the laminar boundary layer is coefficient - rescaled to match the non-dimensionalization sustained farther downstream than the subcritical flow.

adopted by Achenbach [19] - against the circumferential angle.

First, it is noticeable that the skin-friction distribution exhibits a Table 3. Summary of the global quantities predicted by LES for considerable top-bottom asymmetry. This is a clearly the Re = 1.0 x 10', compared with other LES result and the (upstream) influence of the large-scale, alternate vortex- experimental data; CA - Catalano el al.13]; SZ - Szechenyi 120];

shedding. What is most interesting in the figure is the S11 -ShIli etal. [211; Z) - Zdravkovich 116])

appearance of small transient separation bubbles on both upper and lower surface starting as early as. 70 - 75 degrees. CD E. CL -E. St Interestingly, this range largely coincides with the range of LES (present) 0.27 0.12 0.28 -

separation angles reported by many others 12, 19] deduced from LES (CA) 0.32 - 0.32 0.35 the inflection point of mean Cp curve; 77 degree by Cantwell Exp.(SZ) 0.25 and Coles [2], and 78 degrees by Achenbach [19] at Re = 1.0 x Exp.(SH) 0.24 0.33 -

105. Although not shown here, the separation angle determined Exp. (ZD) 0.2/0.4 0.1/0.15 0.2/0.34 0.18/0.5 based on the time-averaged velocity field obtained from the present LES was found to be much larger than the ones determined from C, distribution, reaching around 98 degree, The global flow parameters predicted by the present LES which was also found by Breuer [I]. Based on the are summarized in Table 3. Our prediction of the mean drag instantaneous skin-friction distribution, we surmise that a coefficient (E.) closely matches the data of Szechenyi [20] and transient laminar separation occurs much earlier than the mean Shih et al. 121] interpolated at the present Reynolds number.

separation angle. We consider their data more reliable than others in view of the The distribution of the resolved turbulent kinetic energy sufficiently large spanwise lengths used (LID = 9.3 and 8.0, depicted in Figure 7 supports the possibility of an intermittent respectively), relatively low blockage ratio (7% and I 1%,

laminar separation around 70 - 75 degrees, insofar as the respectively), surface smoothness, and relatively low turbulent kinetic energy is still quite low there. The complete turbulence intensity of the wind tunnels - 0.08 % for Shih et transition occurs a little downstream as indicated by the rapid al.'s data [21]. The r.m.s. lift coefficient prediction by the increase in the resolved turbulent kinetic energy between 75 present LES falls within the scatter of the experimental data and 85 degrees. The magnitude of the skin-friction coefficient compiled by Zdravkovich [16]. The base pressure is also seen predicted by the present LES is also in good agreement with to be in the range of the experimental data.

Achenbach's data [191 measured at Re = 1.0 x ]0 5.

Achenbach's data show that the scaled skin-friction coefficient 6 Copyright C 2005 by ASME

As before, we plotted the circumferential variations of the instantaneous skin-fiiction and the resolved turbulent kinetic energy in Figures I I and 12. The top-bottom asymmetry of the skin-friction distribution is seen to be much less pronounced than for the subcritical flow. Figure 12 suggests that the boundary layer sustains being laminar down to 90 degrees. The separation angle of the timc-averaged velocity field from the 04 . . .i. r \ . . . . A A I. present LES wvas found to be around 108 degrees, which is U A close to 106 degrees estimated by Shih et al.[21] using their experimental data on a smooth cylinder at approximately the same Reynolds number. All this taken together, our LES result suggests that the boundary layer undergoes transition starting as

-0.2- early as 90 degrees, before it finally separates - "turbulently" -

V.

at around 108 degrees.

-0.4 420 . T .

440 460 480 Non-dimensional lime. 1* (= tU01D)

Figure 9. Time histories of drag and lift coefficients for the circular cylinder at Re = 1.0 X lo6 l---Upper surfacel 5.0 [ ower surfaceI The time histories of drag and lift coefficients are shown in Figure 9. Consistent with the flow structure portrayed in Figure 8, the CL-history at this supercritical Reynolds number is more random than that for the subcritical flow. Indeed, no single frequency of vortex-shedding could be identified at this Reynolds number, as shown in the plot of the power spectral density of CL in Figure 10. This finding is consistent with the others' experimental observations, for instance, Shih et al. [21]

who found that coherent vortex shedding disappeared on a smooth cylinder beyond Re = 4 x 105. 0 30 60 90 120 150 180 (deg.)

• I l I . I -. I i l . I. I I

• I Figure 11. Scaled instantaneous skin-friction cocrficient distribution on the cylinder surf ice at Re =1.0 x 106 10-

-3 U

0 IM iVv Ho,1

. .. 2 . . . . . . .

0a-'

..' . . . I.'. .' ....

'u 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Strouhal Number Figure 10. Power spectral density of the lift coefficient (Ct) for the critical Reynolds number (Re = 1.0 x 106) 0 (de".)

For a smooth cylinder, the experimental data collected by Figure 12. Resolved turbulent kinetic energy distribution at Zdravkovich [16] and Achenbach's data [19] indicate that, 0.OOOID above the cylinder surface, predicted the present above the Reynolds number of 1.5 x 106 referred to as "TrBIA" LES at Re = 1.0 x 1O0 in [16], transition clearly precedes the boundary ]ayer separation. The exact location of the transition depends on such factors as freestream turbulence level and surface The skin-friction prediction, including the location of the roughness. At Re = 1.0 x 106 falling in the range referred to as peak and the magnitude, is also in good agreement with the "TrBL3", however, both sets of data seem to suggest that experimental data of Achenbach [19] measured at Re = 8.5 x transition occurs near the separation. 10 and Re = 3.6 x 106. At these two Reynolds numbers, 7 Copyright C)2005 by ASME

Achenbach's data show that the scaled skin-friction University in Turkey for kindly allowing me to use his distributions have the maximum values in the range of 3.5 - 4.0 experimental data.

occurring around 0= 60 - 65 degrees. As shown in Figure 11, our LES prediction closely matches Achenbach's data. This is REFERENCE in a sharp contrast with the LES prediction by Catalano et al.

13], which significantly overpredicted the overall level of the [I] Breuer, M., 2000, "A Challenging Test Case for Large skin-friction on the front half of the cylinder. Contrary to what Eddy Simulation: High Reynolds Number Circular is suggested by the experimental evidence, their LES seem to Cylinder Flow," Int'l J. Beat and Fluid Flow, 21, pp. 648 -

yield too early a transition, giving a skin-friction level typical 654.

of a turbulent boundary layer from immediately downstream of [2] Cantwell, B. and Coles, D., 1983, "An Experimental Study the front stagnation point. As they pointed out [3], one most of Entrainment and Transport in the Turbulent Near-Wake likely cause for the discrepancy is the substantially coarser of a Circular Cylinder," J. FluidMech., 136, pp. 321 - 374.

mesh (ye around 50) used in their LES predictions. [3] Catalano, P., Wang, M., Iccarino, G. and Moin, P., 2003, "Numerical Simulation of the Flow Around a Circular

SUMMARY

AND CONCLUSIONS Cylinder at High Reynolds Numbers," Ini'l J. Hlear cnd Fluid Flow, 24, pp. 463 - 469.

Turbulent flow around a smooth fixed circular cylinder [4] Travin, A., Shur, M., Strelets. M. and Spalart, P., 1999, was computed using large eddy simulation (LES) at two "Detached-Eddy Simulation Past a Circular Cylinder,"

Reynolds numbers (Re = 1.4 x 105, 1.0 x 106). The Flow, Turbulence and Combustion, 63, pp. 293 - 313.

computations were carried out using a second-order accurate, [5] Vatsa, V. N. and Singer, B. A., 2003, "Evaluation of a parallelized finite-volume Navier-Stokes solver capable of Second-Order Accurate Navier-Stokes Code for Detached handling arbitrary unstructured meshes. An implicit, non- Eddy Simulation Past a Cricular Cylinder," AIAA Paper iterative fractional-step method was employed in favor of its 2003-4085.

high efficiency in LES for incompressible flows. [6] Pandya, M. J., Frink, N. T., Abdol-Hamid, K. S., and The LES results for both Reynolds numbers predicted, Chung, J. J., 2004, "Recent Enhancements to USM3D with a good accuracy, the global flow parameters such as mean Unstructured Flow Solver for Unsteady Flows," AIAA drag coefficient, r.m.s. lift coefficient, and the Strouhal number. Paper 2004-5201.

Furthermore, the salient features of the subcritical and the [7] Kim, S. E., Makarov, B., and Caraeni, D., 2004, "Multi-supercritical flows experimentally observed and measured on a Dimensional Reconstruction Scheme for Unstructured smooth cylinder, such as the separation angle, mean velocity in Meshes," AIAA Paper 2004-2548.

the wake, length of the recirculation bubble, and transition [8] Kim, S. E. and Makarov, B., 2005, "An Implicit location, are largely reproduced by the LES. Despite the use of Fractional-Step Method for Efficient Transient Simulation a large number of computational elements (6.8 million cells), of Incompressible Flows," To be presented at 17 1h AIAA the solution turnaround time was quite reasonable. Computational Fluid Dynamics Conference, June 6 - 9, Toronto, Ontario.

To summarize what has been achieved through this study: [9] Park, J., Kwon, K., and Choi, H., 1998, "Numerical Solutions of Flow Past a Circular Cylinder at Reynolds

• The spatial discretization method (finite-volume on Numbers up to 160," KSME lnt. J., 12, pp. 1200- 1205.

unstructured meshes) and the solution algorithm [I0]Norberg, C., 2003, "Fluctuating Lift on a Circular (implicit fractional-step method) employed in this study Cylinder: Review and New Measurements," J. Fluids and have been shown to be sufficiently robust and accurate Structures, 17, pp. 57 - 96.

to be used in LES for high Reynolds number flows. I I]Williamson, C.H.K., 1989, "Oblique and Parallel Modes of

• The dynamic Smagorinsky model adapted to Vortex Shedding in the Wake of a Circular Cylinder at unstructured meshes using a novel test-filter [13] has Low Reynolds Numbers," J. Fluid Mech., 206, pp. 579 -

been shown to work robustly and accurately for 627.

complex, high Reynolds number turbulent flows. [12]Germano, M/, Piomelli, U., Moin, P., and Cabot, W. H.,

• The present LES capability has been shown to predict 1991, "Dynamic Subgrid Scale Eddy Viscosity Model,"

the salient features of turbulent flow around a smooth Physics of Fluids A, 3, 19, pp. 1760 - 1765.

circular cylinder at both subcritical and supercritical [13]Kim, S. E., 2004, "Large Eddy Simulation Using regimes with a good accuracy and reasonable Unstructured Mesh and Dynamic Subgrid-Scale computational cost. Turbulence Models," AIAA Paper 2004-2548.

[14]West, G. S. and Apelt, C. J., 1993, "Measurements of It is planned to continue this work to include higher fluctuating Pressures and Forces on a Circular Cylinder in Reynolds numbers. the Reynolds Number Range 104 to 2.5 x 10'," J. Fluids and Structures, 7, pp. 227 - 244.

[15] Szepessy, S. and Bearman, P. W., 1992, "Aspect Ratio and End Plate Effects on Vortex Shedding from a Circular ACKNOWLEDGEMENT Cylinder, J. Fluid Mech., 234, pp. 191 - 217.

16]Zdravkovich, M. M., 1997, "Flow Around Circular The authors acknowledge that FLUENT and GAMBIT software were used for this study. The authors are also Cylinders, Fundamentals, Vol. 1, Oxford Univ. Press thankful to Dr. Mustafa Sarioglu at Karadeniz Technical (Chapter 6).

8 Copyright © 2005 by ASME

17]Roshko, A., 1961, *Experiments on the Flow Past a Circular Cylinder at a Very High Reynolds Number," .1.

Fluid Mech., 10, No. 3, pp. 345 - 356.

[18]Sarioglu, M. and Yavuz, T., 2002, *Subcritical Flowv Around Bluff Bodies," AIAA J. Vol. 40, No. 7, July.

1[9]Achenbach, E., 1968, "Distribution of Local Pressure and Skin-Friction Around a Circular Cylinder in Cross-Flow up to Re = 5 x I06,,, J. FluidMech., 34, Part4, pp. 625 - 639.

[20]Szechenyi, E., 1975, "Supercritical Reynolds Number Simulation for Two-Dimcnsional Flow Over Circular Cylinders," J. Fluid.Mech., 70, pp. 529 - 542.

[21]Shih, W. C. L., Wang, C., Coles, D., and Roshko, A., 1993, "Experiments on Flow Past Rough Circular Cylinders at Large Reynolds Numbers," J. Wind Eng. And InduIs.

Aerodynamics, 49, pp. 351 - 368.

9 Copyright C 2005 by ASME

BVY 05-072 Docket No. 50-271 Exhibit EMEB-B-138-2 Vermont Yankee Nuclear Power Station Proposed Technical Specification Change No. 263 - Supplement No. 30 Extended Power Uprate Response to Request for Additional Information Large Eddy Simulation of Confined Swirling Coaxial Jets.

Total number of pages in this Exhibit I (excludina this cover sheet) is 10. I

Proceedings of ASME-FED:

2 nd Symposium on Measurements and Modeling of Large-Scale Turbulent Structures June 19 - 23, 2005 FEDSM 2005-77085 LARGE EDDY SIMULATION OF CONFINED SWIRLING COAXIAL JETS Sung-Eun Kim, Fluent Inc., Xuelei Zhu, Fluent Inc., Lebanon, Stefano Orsino, Fluent Inc.,

Lebanon, NH 03766, U.S.A. NH 03766, U.S.A. Lebanon, NH 03766, U.S.A.

ABSTRACT concentration were much less satisfactory. The major discrepancy was found in the core region, where the mean Large eddy simulations (LES) have been carried out for species concentration was grossly underpredicted by the RANS confined swirling coaxial jets discharged into a suddenly computations. Evidently, the RANS model employed gives expanded pipe, which was studied experimentally by Roback rise to an excessive mixing. Brankovic et al. [3] investigated and Johnson [1, 2]. The computations were made using a the sensitivity of their RANS predictions to such parameters as parallelized finite-volume-based Navier-Stokes solver that is the turbulent Schmidt number and the inlet turbulence level.

second-order accurate in time and space, and permits use of However, the influence of these parameters was negligibly unstructured meshes. The computational domain starts from an small. The poor prediction of the species mixing has been inlet placed upstream of the swirl generator, which makes the attributed to the inability of RANS-based turbulence models to inlet boundary condition easy to specify. Subgrid-scale accurately represent the mixing by large-scale turbulent turbulent stresses were modeled using a dynamic Smagorinsky structure (or large eddies).

model applicable to complex three-dimensional flows without Playing a major role in mixing, large eddies are also any statistically homogeneous directions. Subgrid-scale responsible for undesirable yet somewhat subtler phenomena turbulent scalar flux is modeled using a constant Schmidt such as combustion instability and acoustic excitation inside number in conjunction with the dynamically computed subgrid- combustors. To address these issues properly, one has to turn scale turbulent viscosity. The LES predictions were found to to high-level turbulent simulation like LES. Akselvoll and closely reproduce the salient features of the flow and the Moin [5] used LES to compute the flow and the species species concentration downstream of the expansion. One of the concentration in non-swirling coaxialjets. Pierce and Moin [6]

conclusions was that a good resolution of the mean flow and attempted LES for the same confined coaxial swirling jet as is turbulence in the upstream region is crucial in accurately considered here. These studies demonstrated that LES can predicting the mixing downstream of the expansion. indeed predict the flow and the species concentration in confined coaxial jets with a commendable accuracy. LES also allows one to directly compute r.m.s. velocity fluctuations and INTRODUCTION r.m.s. species fluctuation which are important quantities in the context of turbulent combustion.

Confined swirling co-axial jets discharging into a suddenly The present study concerns a LES computation for the expanded pipe [1, 2] have been studied numerically by several swirling coaxial jets using a finite-volume solver that permits others using the Reynolds-Averaged Navier-Stokes (RANS) use of unstructured meshes. Unlike the study of Pierce and equations [3, 4] and large eddy simulation [6]. The flow in Moin [6], however, our computational domain starts from an question consists of a non-swirling jet in the center, and an inlet boundary placed upstream of the swirl generator. Despite outer annual jet with a swirl imparted by a 8-vane, variable- the larger solution domain and the implied increase in the angle swirl generator, with a swirl number of approximately computational resource, having the inlet boundary upstream of 0.41. the swirl generator makes it straightforward to specify the inlet Earlier numerical studies using RANS-based turbulence boundary conditions. Meshing the computational domain models [3, 4] have shown that the mean velocity field in the including the swirl generator, which could become a difficult mixing region - the shear layer between the jets and the inner and time-consuming task with structured meshes, is made recirculation region - was predicted with a reasonably good easier by the flexibility offered by unstructured mesh allowed accuracy. However, the predictions of the mean species by the present finite-volume solver.

I Copyright © 2005 by ASME

This paper is, in many aspects, a progress report about an ongoing study whose ultimate goal is to find an optimal Au=-Gp +r" strategy based on the technique of LES for modeling turbulent AtDGtp"' = DO (3) flow and species mixing in coaxial combustors with or without swirl. u..I = u - ArG,7,+'

The paper is organized as follows. First, we briefly discuss the numerical method and the subgrid-scale turbulence models On a per-iteration basis, the series of operations in for stresses and scalar flux. Special emphasis is laid on a Equation (3) closely resembles SIMPLEC scheme. The transient algorithm whose efficiency significantly benefits the difference is that, in the iterative SIMPLEC scheme, the present study. This is followed by the details of the LES operations in Equation (3) are repeated in an outer loop until computations regarding the choices of the solution domain, the all the solution variables converge, whereas the FSM needs mesh (resolution), boundary conditions, time-step size, and the only one sweep. To preserve second-order accuracy with the overall solution strategy. The results are then presented. FSM, however, sub-iterations are needed for the set of three momentum equations and individual scalar equations to account for the nonlinearity in and coupling among the NUMERICAL METHOD individual equations and high-order source terms. Yet, non-iterative FSM is takes much less CPU time than iterative The computation were carried out using the segregated SIMPLEC scheme.

solver in FLUENT, a general-purpose CFD software. The system of discretized governing equations are solved FLUENT employs a cell-centered finite-volume method based using a point-implicit, Gauss-Seidel relaxation along with on a multi-dimensional linear reconstruction scheme that algebraic multi-grid (AMG) method to accelerate solution permits use of computational elements (cells) with arbitrary convergence. The N-S solver and the SGS turbulence model polyhedral cell topology including quadrilateral, hexahedral, are fully parallelized.

triangular, tetrahedral, pyramidal, prismatic, and hybrid meshes. The solution gradients at cell centers that are needed to compute convective and diffusive fluxes are obtained by SUBGRID-SCALE TURBULENCE MODELING applying Green-Gauss theorem [9]. Diffusive and convective fluxes are discretized using central differencing [8]. For incompressible flows, the filtered Navier-Stokes An implicit fractional-step method (FSM) [10] in equations can be written as conjunction with a second-order accurate, three-level backward- differencing scheme for time-discretization was employed to advance the solution in time. In this algorithm, the momentum equations are decoupled from the continuity awi ,, 2a =_I ap -1IlY+ a Iv2iL (4) equation applying an approximate factorization of the coupled at arir waxy , x ax.4I ax, I Navier-Stokes equations. For incompressible flows, the FSM preserves the formal second-order temporal accuracy without where r,, =uuj- -ii is the subgrid-scale turbulent stress. In having to perform, at each time-step, costly outer iterations to this study, the subgrid-scale stress is modeled using isotropic couple velocity-field and pressure. The FSM thus provides a eddy-viscosity as highly efficient algorithm for CPU-intensive transient computations like LES.

Consider the semi-discrete form of the Navier-Stokes rig- TOr* -Cvl (5) equations in "pressure-correction" (" p"+I -p") form:

WVe determined the model constant, Cn, using the dynamic Smagorinsky model (DSM) originally proposed by Germano et al. [7, 8].

[D 0] 6p- ) 0) The filtered transport equation for a passive scalar is given where u"nt and r are the velocity vector and the momentum by source vector, respectively, and the integer n is the time level.

A is the coefficient matrix defined in terms of the discrete ap+ai7H= a ( aa qj (6) convective and diffusive operators and the time-advancement scheme chosen, and G and D are the discrete gradient and at ax, ax, ax, divergence operators, respectively. The coupled system of equations in Equation (1) is extremely difficult to solve as it is, where q1 is the subgrid-scale turbulent flux of the species since the matrix A has to be inverted for every iteration. In the concentration (V), and a the molecular diffusivity. The subgrid-fractional-step method, the original coupled equations in scale flux is modeled using Equation (I) are approximated by EDA 0DG][I A ]( u) =[ (r+ [°(At )]

(2)

Sc, ax, (7)

ID A:DGJ0 J

] 6p-+ )+L~

Factorizing equation (2), we have a series of split operations as 2 Copyright © 2005 by ASME

where Sc, is the subgrid-scale turbulent Schmidt number. In this study, we used a constant value of 0.9.

To implement the dynamic procedure for the present finite-volume solver requires a test-filter applicable to arbitrary unstructured meshes. The test-filter finally adopted for this study is a top-hat filter involving a volume comprising the cell in the center plus the neighboring cells that share the faces with the center cell. To make the dynamic procedure tractable, an approximation was made that is tantamount to a non-varying Cv over the test-filter volume. Thus, the dynamic procedure employed in this study is "approximately local" in the sense that, despite the ad hoc assumption, it does not require an existence of any statistically homogeneous directions, being applicable to complex three-dimensional flows. To avoid numerical instability Figure 2. A sectional view of the computational mesh (medium-sized mesh with 2.7 million elements) used in the likely to be caused by a large fluctuation of the model constant, we smoothed the model constant by applying the test-filter on it, computations.

and also clipped it so that the effective viscosity remains positive.

The additional details of the implementation of the DSM in the The computational meshes were built using Gambit. We framework of an unstructured mesh based finite-volume solver used three progressively refined, unstructured hexahedral can be found in the reference [8]. The DSM has been validated meshes for the computations. The three meshes have 1.1 for a number of wall-bounded flows such as fully-developed million (coarse), 2.7 million (medium), and 4.8 million (fine) channel flows and flows around bluff bodies such as a square elements, respectively. A sectional view of the medium-size cylinder and a sphere. mesh is shown in Figure 2. The resolution of the medium mesh is such that, in the mixing region, it can resolve the integral DETAILS OF COMPUTATIONS length-scale estimated around 20 mm based on the experimental observation [1, 2], which is equivalent to D16, with Solution domain, swirl vane geometry and meshes approximately 20 elements. It can resolve the smallest observable eddies of a size around 6 mmn, which were observed A partial view of the computational domain is depicted in in the experiment [ I, 2], with 6 - 7 elements.

Figure 1, along with the coordinate system. The inlet boundary It is worth mentioning here that the mesh resolution in and is placed at L.OD upstream of the swirl generator, where D is around the inner and the annular tubes and the swirl vanes is the downstream pipe diameter (D = 122 nmn). The exit too coarse to resolve the energy-containing eddies originating boundary is at 15D downstream of the pipe expansion. The from and transported in that upstream region. We recognize computational domain thus has the swirl generator in it. The that poor resolution of the turbulent structure in the upstream entire inlet tube, 8 swirl vanes, and the downstream pipe were region will negatively impact the prediction in the downstream, modeled without using any periodic boundaries. particularly in the shear layers of the mixing region, inasmuch as the energy-carrying eddies from upstream feed the turbulent structure being developed in the downstream mixing region.

pipe expansion Furthermore, none of the three meshes are fine enough to accurately resolve the viscous sublayer on the pipe wall downstream of the expansion. This was deemed justifiable in light of a relatively passive role played by the wall downstream swuirvane of the expansion in the mixing occurring near the central core.

. - '7:, .', - ' Boundary conditions, time-step size, and solution strategy I I I - - _.. t -,

i-_ - -4:.

.: On the annular jet inlet, a uniform axial velocity of 1.667 ni/s was specified according to the data given by Roback and Johnson [1]. For the velocity boundary condition at the inner inlet XLZ jet inlet, we used 0.797 ni/s, a value derived from the volumetric flow-rate (6.2 gallon per minute) given in Ref. [I],

Figure 1. A partial view of the cor npulational domain instead of the inlet velocity mentioned in the same report (0.52 mis). The inlet velocity of 0.52 mis quoted in the report does The swirl generator consists of 8 identical vanes mounted not match the given volume flow-rate, and is inconceivably too on the hub with an equal spacing in the azimuthal direction. low considering that the measured axial velocity immediately The blade-section has a camber an(d a vane-angle that change downstream of the expansion is around 0.8 ni/s.

with the radius. The blade geoi netry was taken from the On wall boundaries, we employed a blended law-of-the-design data given in the Ref. [1]. iA NURB surface was built wall that invokes proper wall-laws depending on the local mesh based on the digitized surface geome try. resolution, namely ye at the wall-adjacent cells. The downstream exit boundary was treated as a "pressure-outlet" boundary offered in FLUENT. In essence, the solution 3 Copyright © 2005 by ASME

variables are extrapolated in a mass-conserving manner on this The LES was started using a steady RANS solution as the boundary. initial guess. To accelerate the solution to a statistically The influence of the exit boundary condition was stationary state, we superimposed a pseudo-random fluctuating discussed at length by Pierce and Moin [6]. They found that velocity-field on the mean velocity field taken from the RANS the usual convective outflow boundary condition applied on an solution. Before the statistics are collected, the LES exit boundary with a cross-section of the downstream pipe computation was run until the initial transients in the solution yielded a central recirculation zone that is far smaller than what are completely washed out, which typically took 2 - 3 flow-through times; one flow-through time is taken to be LlUb, where L is the length of the downstream pipe, and Ub the axial mean bulk velocity.

RESULTS Figure 4. Instantaneous velocity vectors in the r-z plane

-x Figure 3. Instantaneous velocity vectors at four axial locations - top-left, z = 25 nim; top-right, z= 51 mm; bottom-left, z=102 mm, boitom-rigltt,z=203mnm was observed in the experiment. They were able to bring their LES prediction much closer to the experimental data by putting a second expansion before the exit boundary, which was conceived based on an argument that it better represents the settling chamber used in the experiment.

We also investigated the effect of the second expansion, comparing the result to what was obtained without he Figure 5. Time-averaged velocity vectors in the r-z plane expansion. However, we could not notice any appreciable difference between the two. The insignificantly small difference found in our computations is perhaps due to the way To obtain stable statistics of the solution, the transient the solution variables are extrapolated at "pressure-outlet" computations were continued for a sufficiently long period of boundary that could be different from their treatment of the time, typically for more than 7 - 8 flow-through times, until convective boundary. the time-averaged velocity and species concentration fields The time-step size used in the present LES was determined largely recover an axisymmetry.

based on the estimated characteristic time-scale of the smallest Regarding the mesh-dependency of the solutions, it was eddies to be resolved in the LES. We took v = 1.0 m/s and t found that the LES results from the medium mesh (2.7 million

= 0.05D (D = 0.122 m) as the characteristic velocity-scale and cells) and the fine mesh (4.8 million cells) showed little length-scale of the smallest eddies. These estimates give an difference, while the coarse mesh result deviates a little farther eddy-turnover time of r = /v = 0.006 second. It was finally from other two. Unless stated otherwise, the results presented decided to use a time-step size of At = 0.0002 second, which in this paper are for the medium mesh.

will resolves one eddy-turnover time (i) of the smallest resolved eddies roughly in 30 time steps.

4 Copyright C)2005 by ASME

the penetration length of the center jet with high species concentration (colored pink in the figure in Figure 7) was found to be around 50 mm, which is closely reproduced by our LES results as shown in Figure 7.

Figure 6. Instantaneous velocity vector plots at two axial locations - left, z = 25 mm; right, z= 51 imm Overall flow structure Figures 3 and 4 show the instantaneous velocity vectors projected on four crossflow (r-0) planes and a r-z plane. The vector plots portray turbulent eddies with widely varying length-scales throughout the mixing region. Small eddies are Figure 8. Vortical flow structure in the mixing region shown to form in the shear layers between the inner jet and the visualized by the !so-surface of the second-invariant of the outer annular jet, growing in size in the downstream direction. velocity deformation tensor, colored by velocity magnitude.

At z = 25 nim, small eddies are confined near the shear layer between the inner and the outer annual jets. Yet, larger eddies Figure 8 gives an overall impression of the turbulent are also seen to have formed in the annular recirculation region vortical structure in the mixing region. This figure indicates behind the step. One can visually tell from the figure that the that the flow in this coaxial jets are well mixed.

smallest structure at z = 25 nim resolved in the LES is roughly D/20 (D = 122 mm), which is close to the size of the smallest Velocity field eddy, 6 mm, observed by Roback and Johnson [I]. The size of the largest eddy at z = 25 mm, which is about two-thirds of the Figure 9 shows the mean axial velocity profile along the step-height as shown in the figure, also closely matches the centerline of the pipe. The LES prediction quite closely experimentally observed value, 20 imm, quoted in ref. [1]. reproduces the general trend such as the rapid drop of the mean Figure 5 depicts the mean velocity vectors on a r-z plane. axial velocity immediately downstream of the expansion and the gradual recovery further downstream. However, the recirculation zone predicted by the present LES appears to be shifted slightly downstream compared to what is indicated by the measurement.

1.5 LES (present) 1.0 o Measured (Roback and Johnson, 1983)1 0.5 0 0

0. 00 Figure 7. Contours of the instantaneous species concentration in r-z plane Figure 6 and 7 depict the contours of the instantaneous Figure 9. Mean axial velocity along the centerline of the species concentration on two crossflow planes and a r-z plane. pipe.

The contours on the crossflow planes give an idea of the length-scale of the turbulent structure. The observations from these figures are consistent with those from the instantaneous velocity vector plots discussed earlier. In the experiment (I],

5 Copyright © 2005 by ASME

z=5mm z=5 mm 8

.9U C C

'U

E z=25mm

>1

>1 C

Ma U U9 U4 z=51 mm i 2:

.C V

Ca Ut a

Ca

'4 r/rO r/rO Figure 10. Mean axial and azimuthal velocity profiles at two axial locations (z = 25 mm, z = 51 mm)

Figure 10 shows the radial profiles of the mean axial and obtained by averaging the radial profiles taken at four azimuthal velocity components at three axial locations (z = 5 azimuthal locations (0= 00, 90°, 180°, 270°). The predictions mm, 25 mm, 51 mm). The profiles shown in the figure were are seen to capture the overall trends and the peak values of the 6 Copyright C)2005 by ASME

velocity components. However, it is clearly noticeable that the used in this study, including the fine mesh (4.8 million cells),

locations of the peaks predicted by the LES are generally are fine enough to accurately resolve the energy-containing shifted toward the centerline. This implies that the outer eddies generated in the upstream region.

annular jet expands or opens up less than it does in reality. The r.m.s. axial velocity fluctuation shown in Figure 11 This also means that the annular recirculation zone predicted by seems to support our reasoning. The r.m.s. axial velocity the present LES is longer (in the axial direction) than in the fluctuation at the centerline (r = 0) at z = 25 mm shown in the experiment. The underprediction of the jet angle is also top-figure, and the one at z = 5 mm (not shown here) are correlated with the downstream shift of the central recirculation severely underpredicted in the LES computation. As shown in zone discussed earlier. the bottom-figure for z = 51 mm, further downstream of the jet exit, the r.m.s. axial velocity fluctuation catches up with the z= 25 mm data, as the turbulent eddies generated in the shear layer become full-fledged.

The impact of the incoming turbulent eddies on the mixing in the downstream as suggested above and the cost implication of using an extremely fine mesh in the upstream region beg a question; what would be the best practice that can be adopted to obtain a sufficiently accurate prediction of the flow and species mixing occurring in a coaxial jet combustor with LES? We will ponder a little upon this question at the very end.

Species concentration Figure 12 shows the profile of the mean species concentration along the pipe centerline. The LES prediction closely reproduces the trend - the plunge of the mean species concentration occurring near z = 50 mm. One noticeable i discrepancy between the prediction and the measurement is in the length of the inviscid core for the species concentration. The LES yields an inviscid core for the species concentration almost down to z = 35 mm, whereas the experimental data indicates that the mean species concentration starts being diffused away almost immediately after the expansion. We think that this is again, in a major part, due to the poor resolution of the incoming turbulent E eddies which would take part in "tearing" the inner jet.

Z-

.9 a1 C1 rlr.

Figure 11. r.ms. axial velocity fluctuation at two cross-nlow planes (z = 25 mm, z = 51 m)

We surmise that, among others things, the most likely 50 100 150 culprit for this discrepancy is the lack of mesh resolution in the Axial distance from the expansion (mm) region upstream of the expansion, in and around the inner and the outer annular pipes, and the swirl vanes, As mentioned Figure 12. Mean species concentration along the centerline earlier, turbulent eddies coming from the upstream region feed of the pipe the shear layers developed downstream of the expansion, enhancing the mixing of the momentum and the species concentration in the inner and the annular jets with the The radial profiles of the mean species concentration at two surrounding flow, which will lead to an increase in the jet axial locations are shown in Figure 13. The predictions are in angle. In our LES computations, these energy-feeding eddies good agreement with the measurements at both locations. The are almost missing, since none of the computational meshes mean species concentration profile at z = 25 mm predicted by the LES computation shows a sign of being "under-diffused", insofar 7 Copyright ©2005 by ASME

as it has a fuller profile than the measured one. The prediction at r.m.s. value at z = 25 nmn, located in the shear layer between the z = 51 mm is somewhat lower than the measurement near the inner and the outer annular jets, is also considerably centerline. However, the overall agreement at this location is underpredicted. Apparently, the LES underpredicts the remarkably good in view of the steep change of the mean species entrainment of the ambient fluid occurring at this location. At z concentration near z = 51 mm, as can be seen in Figure 12. It = 51 mm, the prediction comes much closer to the should be noted, in passing, that the RANS predictions reported measurement, as the r.m.s. velocity fluctuation did at the same in the literature 13] severely under-predicted the mean species axial location.

concentration at this axial location.

z =25mm z =25mm z=51 mm z =51 mm 0

(fJ 0

U C

r/ro rlro Figure 13. lean species concentration at two axial Figure 14. r.m.s. species concentration at two axial locations locations (z = 25 mm, 51 mm) (z = 25mm, 51 mm)

Equally important - probably even more important than the The same remarks given regarding what might have caused mean species concentration in the context of modeling the discrepancy in the r.m.s. velocity fluctuation (Figure 11) turbulent combustion, is the fluctuation of species largely apply to the results for the r.m.s. fluctuating species concentration. The usual RANS-based turbulence models concentration.

cannot directly predict the r.m.s. fluctuating species concentration, unless the transport equation for the variance of fluctuating species concentration is explicitly solved. One PRELMINARY RESULTS WITH LOCALLY REFINED obvious benefit from LES is that one can directly compute it. MESHES Figure 14 depicts the r.m.s. fluctuating species concentration at the two crossflow planes. At z = 25 mm, as in the case of the The fact that the fine mesh (4.8 million cells) offers a r.m.s. axial velocity fluctuation, the LES result grossly meager improvement over the medium mesh (2.7 million cells) underpredicts the r.m.s. value in the core region. The peak in terms of the mesh resolution (spacing) warrants additional 8 Copyright ©2005 by ASME

computations with substantially finer mesh, for instance, by halving the grid spacing in all three (axial, radial, and CONCLUSION azimuthal) directions. However, doubling the number of elements in the three directions results in an eight-fold increase The flow and the species transport in confined swirling in the total cell counts, which becomes unwieldy. coaxial jets were computed using LES. An unstructured mesh-based finite-volume solver was employed for the computations.

A highly efficient time-advancement scheme was used in z=51 mm conjunction with second-order accurate temporal and spatial

'5..

discretization schemes. A dynamic Smagorinsky model

_ I Measured (Roback & Johnson [1]) adapted to general three-dimensional flows was employed as 2.0 - LES (M.7Mcell mesh) the subgrid-scale turbulence model.

-- LES (2.7M cell & locally refined)

For LES, it is evidently a bold attempt and a costly 1.51 proposition to include the upstream components of a combustor like the swirl generator and the upstream tubes in the 1.0 computational domain. Nevertheless, the present LES computations closely reproduce the salient features of the flow 0.5 and the species concentration in the mixing region. For an accurate prediction of the mixing in the downstream (e.g.,

combustion chamber), a good resolution of the mean flow and 0.0 the turbulence in and around the inner and the outer annular

' a 0~ B1 X tubes, and the swirl vanes turned out to be more important than

-0.5 originally thought. The numerical evidence found from this I , I , I I oo I 0.2 0.4 0.6 0.8 study indicates that the three globally refined meshes used in rTrN the present study, despite the largest cell counts reaching up to 4.8 million cells, still cannot provide a sufficient resolution of z=Sl mm the upstream region. We believe that the lack of mesh resolution is responsible for the discrepancy between the predictions and the measurements, most notably the overall shift of the velocity peaks toward the centerline.

Finally, a preliminary result was presented which was obtained with a new mesh locally refined in the upstream part of the domain. The significant improvement from this locally adapted mesh supports our conclusion, and at the same time, provides us with an avenue to improving the accuracy of LES prediction for the subject flow.

ACKNOWLEDGEMENT The authors acknowledge that FLUENT and GAMBIT software were used for this study. The LES-cluster in Fluent r/rO Inc. was used for the computations.

Figure 15. Mean axial and azimuthal velocity profiles at two axial locations - predicted with a locally refined mesh REFERENCES To keep the mesh size under a tractable limit, the medium [I] Roback, R. and Johnson, B. V., "Mass and Momentum mesh was locally refined only in the domain upstream of the Turbulent Transport Experiments with Confined Swirling expansion. Furthermore, only the cells within a specified Coaxial Jets," NASA CR-168252, 1983.

distance from the wall, roughly 0.07D in this study, were 12] Johnson, B. V. and Roback, R., "Mass and Momentum refined to further save the cell counts. The local refinement Turbulent Transport Experiments with Confined Swirling resulted in a total of 6.5 million cells. The computation has Coaxial Jets - Part 1," AIAA-84-1380, Presented at been carried out on this mesh. The results are shown in Figure AIAA/SAE/ASME 2 0 th Joint Ppropulsion Conference, 15, being compared with the medium mesh results, for the Cincinnati, OH, June 11 - 13, 1984.

radial profiles of the mean axial and azimuthal velocity [3] Brankovic, A., Ryder, Jr. R. C., and Syed, S. A., 1998, components at z = 51 mm. Clearly, the predictions with the "Mixing and Combustion Modeling for Gas Turbine adapted mesh are significantly improved over the predictions Combustors Using Unstructured CFD Technology," AIAA with the medium mesh. The overall shift of the velocity peaks Paper 98-3854, Presented at 341h AIAA/ASME/SAE/ASEE is now far smaller than what we saw earlier, and the predictions Joint Propulsion Conference & Exhibit, July 13 - 15, 1998, capture the profiles of the mean velocity components much Cleveland, OH.

more closely.

9 Copyright C 2005 by ASME

(4] Lin, C. A., "Modeling a Confined Swirling Coaxial Jet,"

Annual Research Brief, 1998, Center for Turbulence Research., Stanford University, 1998.

[5] Akselvoll, K. and Moin, P., 1996, "Large-Eddy-Simulation of Turbulent Confined Coannualr Jets," J. Fluid inech.,

315, pp.387-411

[6] Pierce, C. D. and Moin, P., 1998, "Large Eddy Simulation of a Confined Jet with Swirl and Heat Release," AIAA Paper 98-2892.

[7] Germano, M/, Piomelli, U., Moin, P., and Cabot, W. H.,

1991, "Dynamic Subgrid Scale Eddy Viscosity Model,"

Physics of Fluids A, 3,19, pp. 1760- 1765.

[8] Kim, S. E., 2004, "Large Eddy Simulation Using Unstructured Mesh and Dynamic Subgrid-Scale Turbulence Models," AIAA Paper 2004-2548.

[9] Kim, S. E., Makarov, B., and Caraeni, D., 2004, "Multi-Dimensional Reconstruction Scheme for Unstructured Meshes," AIAA Paper 2004-2548.

[10]Kim, S. E. and Makarov, B., 2005, "An Implicit Fractional-Step Method for Efficient Transient Simulation of Incompressible Flows," To be presented at 17 h AIAA Computational Fluid Dynamics Conference, June 6 - 9, Toronto, Ontario.

10 Copyright © 2005 by ASME

BVY 05-072 Docket No. 50-271 Exhibit EMEB-B-138-3 Vermont Yankee Nuclear Power Station Proposed Technical Specification Change No. 263 - Supplement No. 30 Extended Power Uprate Response to Request for Additional Information Large Eddy Simulation Using Unstructured Meshes and Dynamic Subgrid-Scale Turbulence Models.

Total number of pages in this Exhibit l I (excludino this cover sheet) is 17. I

34th AIAA Fluid Dynamics Ccrterence end Exhibit. June 28 - July 1 Pcrtland Oregon Large Eddy Simulation Using Unstructured Meshes and Dynamic Subgrid-Scale Turbulence Models Sung-Eun Kim*

Fluent Inc. Lebanon, New Ilampshire, 03766, U.S.A.

This paper concerns development of a large eddy simulation (LES) capability based on a finite-volume solver that permits use or unstructured meshes. The solver employs a cell-centered scheme along with a multi-dimensional linear reconstruction. Convection and diffusion terms are discretized using a second-order central-differencing scheme. A three-level second-order scheme is used for temporal discretization. Subgrid-scale turbulent stresses are closed using dynamic Smagorinsky model and dynamic turbulent kinetic energy transport model. A test-filter was designed for the dynamic procedure which is applicable to unstructured meshes of arbitrary cell-topology. The dynamic procedure also avails itself to three-dimensional flows without any statistically homogenous directions. Wall boundary conditions are imposed using a wall-function approach that applies appropriate wall-laws depending on near-wall mesh resolution. The LES capability is validated for a wide range of wall-bounded flows. We present here the results for a fully-developed channel flow and two bluff-body flows. The predictions are in good agreement with direct numerical simulation (DNS) results and the experimental daa.

1. Introduction v

W JE encounter many industrial applications of computational fluid dynamics (CFD) where the flows are dominated by unsteady, large-scalc coherent structures. Those large-scalc structures impact, to a great extent, various aspects of the flows such as energy consumption, safety, comfort, and noise. The ramification of whether or not one can harness the large-scale structures is therefore quite significant. Attempts to numerically predict such flows using unsteady Reynolds-Averaged Navier-Stokes (URANS) equations have been met with limited success. Some of the better RANS models seem to be capable of capturing the "largest" scale occurring often in the form of alternate vortex-shedding. However, the remaining coherent structures are left largely unresolved.

Large eddy simulation (LES) is fundamentally suited to the task of predicting coherent structures. The major obstacle in using LES for practical high Reynolds-number (Re) flows, from the practitioners' standpoint, is its high cost incurred by an unwieldy number of computational elements and painfully long solution time. Yet today's ever-increasing computing power is rapidly making LES feasible. Another difficulty often encountered when attempting LES for industrial applications comes from complex geometry. In CFD, meshing for industrial applications involving complex geometry by itself can become a grand challenge. It is hugely time-consuming or often impossible to generate high-quality structured meshes for complex configurations, which has led industrial CFD practitioners to opt for unstructured meshes. Although unstuctured meshes are routinely used today in RANS computations for industrial applications, attempts to conduct LES with unstructured meshes are just starting to appear in the literature." 5 As yet the efficacy of unstructured meshes for LES for practical high-Re flows has not been fully established. Among the issues yet to be addressed are numerical accuracy, stability, and subgrid-scale turbulence modeling.

This paper is concerned with evaluating a LES capability developed using a finite-volume solver based on second-order numerics. Permitting use of unstructured meshes, the resulting LES capability can easily handle complex ge-ometnes encountered in industrial applications. In addition, it lends itselt to local mesh adaptation that can be utilized

• Principal Engineer, Fluent Inc., Lebanon N.H., Member AIAA.

Copyright 2004 by Fluent Inc.. Published by the American Institute of Aeronautics and Astronautics, Inc. with permission.

I of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

to efficiently allocate computational cells, substantially reducing the computational cost. Adequacy of second-order spatial discretization for LES has often been questioned, and there are some misgivings about using it for LES. It will be shown in this paper, however, that the second-order central differencing scheme adopted in the present work yields a commendable accuracy for LES. For subgrid-scale (SGS) turbulence closure, we implemented two dynamic SGS vis-cosity models in the framework of unstructured meshes, namely, the dynamic Smagorinsky model originally proposed by Germano et aO., 6 and Lilly 7 and the dynamic turbulence kinetic energy transport model of Kim and Menon. 8' 9 For the dynamic procedure, a test-filter readily applicable to unstructured meshes was designed. The resulting dynamic SGS models can be used for three-dimensional flows without any statistically homogeneous directions.

The paper is organized as follows. WVe start with a brief description of the two dynamic subgrid-scale turbulence models and the details of their implementations. This will be followed by an overview of the numerical methods and algorithms adopted in this work. Finally, validations will be presented for a selected number of wall-bounded flows ranging from a fully-developed channel flow to a couple of bluff-body flows including one around a square-cylinder with salient edges and another past a smooth sphere.

II. Filtered Navier-Stokes Equation and Subgrid-Scale Turbulence Modeling A. Implicit filter with finite-volume discretization The governing equations for LES are generally obtained by filtering the Navier-Stokes equations in either Fourier (wave-number) or physical space. In the present work, the filtering operation (denoted by an overbar) is defined as a spatially convoluted integral of the variable in question as t(x) = fl t(y)G(x,y)dy (I) where D is the computational domain, y E D, and G is the filter function.

With the cell-centered finite-volume discretization and the linear reconstruction scheme employed in this work, the discrete solution variable at a cell-center (co) can be written as 4P(co) = V fl(y)dy, y E v (2) where V is the volume of a computational cell.

The definition in Equation (2) of a discrete solution variable at ccl center can be interpreted as a filtering operation

¢(x)- (co) = Vl y Ee (3)

The implied filter function, G(x,y), is then a top-hat filter G(xy) = /V for Ix- y E V (4) 0 otherwise Using the filtering operation in Equation (3), the filtered Navier-Stokes equations for incompressible flows (as-sumed for brevity) can be written as a

a.

_i+

axj

_ I o atij Paxi axj axj HODj

((aa5) ax; = °(6) where ¶ij is the subgrid-scale stress defined by;

,ij _7-j i-j (7) which is unknown and needs a closure.

Thus, we used in this work grid and finite-volume discretization as an implicit filter.

2 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

B. Subgrid-scale (SGS) stress models The SGS models based on the concept of isotropic eddy-viscosity compute the SGS turbulent stress from Tij tkkttoy = -2v,3ii (8) where v, is the SGS eddy-viscosity, and Sij the resolved rate-of-strain tensor defined by 2 \axj ax,,

The task of SGS turbulence modeling is to express SGS viscosity, v,, as a functional of known quantities. In the present work, we employed two dynamic SGS eddy-viscosity models. They are described below.

1. Dynanzic Smagorinsky Model (DSM)

The underpinning of DSM is the algebraic eddy-viscosity model originally proposed by Smagorinsky.10 In the original Smagorinsky's model, SGS eddy-viscosity is computed from

= c.V v, JSJ (9) where Cv is a model constant (Cv = 0.1 0.2), S , / the modulus of rate-of-strain tensor of the resolved scales, and a _ V'- 3.

The subgrid-scale stress can then be written as T

ij-- k = -2CvA ,13 S3j (10)

Despite its simplicity, this model has several shortcomings. The most problematic one, from a practical standpoint, has to do with the model constant, Cv. There is no single value of the constant which is universally applicable to a wide range of flows. Another serious drawback is that the SGS viscosity model in Equation (9) with a constant value of Cv is not applicable to transitional flows where the flow in question is laminar either locally or intermittently, since Equation (9) always gives a finite SGS viscosity even in laminar region as long as there is velocity gradient.

Germano et aL6 and subsequently Lilly7 conceived a procedure that resolves these problems. In this procedure, Cv is dynamically computed during LES, on-the-fly, using the information provided by the smaller scales of the resolved (known) fields. To separate the smaller scales from the resolved field, the dynamic procedure needs a so-called "test-filter" having a width (A) that is larger than the grid-filter width (A). Denoting the test-filtered variables by a tilde and putting the "grid-filtered" Navier-Stokes equations through the test-filter, we obtain "test-filtered" Navier-Stokes equations as a:i. ai,-, I ap azi1 a a..,

. a pa. av a aXJ

-0 (12) ax, where Tj is now a "subtest-scale" stress defined by Trj - - (13)

The underpinning of the dynamic modeling is a similarity concept that T.j, the subtest-scale stress, can be written as a functional of the larger resolved scales in a manner similar to 'rj,6 which leads to 3 of t7 American Institute of Aeronautics and Astronautics Paper 2004-2548

- Bij (14) 2 Tij-TTkk = -2ca Isis (14) where a is the test-filter width. Note that the same model coefficient, Cv, is used in the expressions for both -Tjand Tlj.

Equation (14) alone is not helpful at all in determining Cv, because Tij is not known in LES. The breakthrough came from the realization that 'rij and Tij are related to each other by6 Th ij = uw-Fhii-

= Itiij - Th1tj= Lij (15)

The stress, Lj, which might be called subtest-scale Leonard stress, can be interpreted as the stress associated with the smaller resolved scales between the test-filter scale (a) and the grid-filter scale (A). Since Lij can be directly computed from the resolved scales in LES, the identity given by Equation (15) can be used to determine Cv. Thus, we have Lii 3iLk = aXjCv - ,}cv (16) where atij = -2A2 1S1S 11 (17)

= -2A 2ISISij (18)

One predicament that makes it difficult to determine Cv from Equation (16) is the fact that Cv in the second term on the right-hand side of the equation is under the test-filtering operator. As it stands, Equation (16) is an integral equation for Cv as discussed at length by Ghosal et aL. This difficulty can be avoided by taking out Cv from the test-Tilter operation as LU - 3 Lok CV (aij - F) (I19)

We followed this rather ad hoc approach in this work despite its mathematical inconsistency, which amounts to as-suming that Cv remains constant in the computational cells associated with the test-filter.

Since there are more equations in Equation (19) than the unknowns, the model constant Cv is obtained by seeking for Cv which minimizes the error norm defined by E = (Lij- °J Lkk - CvMij) (20) where

-aij- ij = -2 ( S2sl A-a2 1 jj)

Taking DEfDCv and setting it zero, we obtain Cv = _Mj_ (21)

Cv determined in this way varies with time and space. In fact, it varies in a wide range, often taking either large negative or positive value. Although negative Cv and consequently negative eddy-viscosity is often interpreted as representing "back-scatter" (flow of energy from smaller to larger scales), too large a negative eddy-viscosity causes numerical instability, ultimately leading to divergence of numerical solutions. The usual remedy for this numerical difficulty is to average Cv in statistically homogeneous directions. Obviously, this workaround can be exploited only when there are such homogeneous directions, which is a rarity in practical applications. Even if there are any statisti-cally homogeneous directions, the averaging becomes extremely cumbersome with unstructured meshes. In the present 4 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

work, we simply "condition" the Cv computed by Equation (21) by test-filtering it. This simple volume-weighted aver-aging better preserves the locality of the model constant. To ameliorate the potential numerical difficulty caused by Cv being negative for an extended period of time, we evaluated two alternative approaches. In one approach, we simply clip Cv at 0 when it becomes negative. This option therefore rules out any chance for the model to mimic backscatter.

In another, effective viscosity (V+ v,) instead of Cv is clipped at zero, permitting small negative SGS viscosity to happen. As yet we do not have any conclusive evidence that supports superiority of one approach to another in terms of prediction accuracy, except an indication that computations with the first approach appear numerically more robust.

For this reason, we used the first approach in the computations presented in this paper.

As mentioned earlier, the dynamic procedure described above requires a test-filter. The most important criterion the test-filter should satisfy specifically for this work is that it should be applicable to unstructured meshes of arbitrary cell topology without incurring unduly high cost. The test-filter finally adopted is a top-hat filter that involves a volume comprising the cell itself plus the neighboring cells that share the cell faces with the center cell. Thus, the test-filter operation amounts to a volume-weighted averaging of the variable in question, which is easily implementable and takes advantage of the data structure of the underlying finitc-volumc solver. With hexahedral meshes, the ratio of the test-filter to the gird-filter scale (Z/Z) is approximately 2.1 (= 91/3). The ratio for tetrahedra is smaller, being around 1.7(= 51/3).

2. Localized Dynanic Kinetic Energy Mlodel (LDKEM)

The dynamic Smagorinsky's model described so far is essentially an algebraic model in which subgrid-scale stresses are modeled using the resolved velocity field. A more elaborate SGS stress model would be the one which is directly based on SGS turbulence and can be used to parametrize SGS stresses. The most widely used ones among others are what can be called "one-equation" models in which SGS turbulent kinetic energy, k-8g. = (Tt2 - Uq)/2, is explicitly computed by solving its transport equation. 8 . 9. 11.12 In the present.work, a localized dynamic k5g-equation model of Kim and Menon8 .9 was chosen in favor of its overall cfficacy.13 .14 In the LDKEM, the subgrid-scale eddy viscosity, V,, is computed from V, = CkA (22)

Consequently, SGS stress is can be written as

- 2ij

- ksg, =-2CAkQ, ASij (23) k585 is obtained by solving its transport equation8 ' 9 SRS _____ _

-,r2 C + UV +V ) (24) ax ax, A ax aXJ The only difference between the original formulation and the present one is that the contribution from the molecular diffusion of kss, is included in this work.

As shown above, there are three model parameters appearing in these equations, namely, Ck, CE, and Ok. which need to be specified. In the current implementation of the LDKEM, the first two are determined from the dynamic procedures to be described in the following, whereas crk is simply set to a constant value of 1.0. The underpinning of the dynamic procedure employed in the LDKEM is the hypothesis corroborated by the experimcntal cvidcnccl5. 16 that there is a strong correlation between the subgrid scale stress, tij,and the subtest-scale Leonard's stress, Lij. In place of parameterizing Tij and utilizing the Germano's identity as was done in the DSM, the LDKEM models Lij directly as Lj - 3!Lkk = -2CkAk,'/,,Sij (25) 5 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

where k,,, is the resolved turbulent kinetic energy associated with the scales between the test-filter (Z) and the grid-filter (A). It can be directly computed from kes Itkl M lk Ilk (26)

The model parameter, Ck, can then be determined from Equation (25) by minimizing the error norm as in the DSM.

Consequently we have CA = LjjMjj (27) with Mij defined by Mij = -2Ake! Suj (28)

The model parameter, CE, of the dissipation term in Equation (24) is also determined by a dynamic procedure, whose underpinning is the hypothesis that the dissipation-rate of keen (e) can be expressed in the same functional form as the dissipation-rate of kz,.

e = CckEST (29)

The dissipation-rate of kel can also be computed from e = (v +v,) aaxJi ax}-, axi ax (30)

From Equation (29) and Equation (30), CE is given in a closed form C k(t/,)

kest axi at,-.;;at

(.Axj axi axi axi a-;) (31)

The DKEM has several desirable attributes that the DSM lacks. First, as a consequence of parameterizing Lij directly, the dynamic procedure in the LDKEM does not involve any test-filter operation on the model parameter, Ck.

Thus, unlike the DSM model, Ck is a genuine, local quantity free from any mathematical inconsistency. Secondly, Ck in the LDKEM behaves numerically more benignly than Cv in the DSM, having much less fluctuation. As a small premium, one can even save a small amount of computational effort with the LDKEM, inasmuch as the test-filtering on the SGS stress done on in the DSM is not needed. Lastly and probably most important, the LDKEM enioys the benefits of a high-order turbulence model. Adopting SGS turbulent kinetic energy to parametrize the SGS stress renders the LDKIEM better suited to non-equilibrium flows. By accounting the kg,,-budget more rigorously, backscatter of kinetic energy is allowed in the LDKEM on a much sounder physical basis than in the DSM.

III. Boundary conditions Wall boundaries are the most crucial and yet difficult ones to handle in LES. In LES resolving all the way down to viscous sublayer, no-slip would be clearly the choice for wall boundary condition for the resolved velocity field.

However, the cost of such "wall-resolving" LES is prohibitively high for practical flows involving high-Re flows. A practical alternative is to use the law-of-the-wall bridging the wall and the first grid point (cell center) off the wall. The simplest implementation of the wall-law can done using

_+ = + fory+ <y+ (32) lnE,)for y+>y where E = 9.793, ic = 0.419, y+ _ 77,y/v, u+- 7/Ii, y is the "cross-over"at which the two wall laws intersect. This demarcation of the entire inner layer into the two distinct layers is apparently a simplification which is at odds with the 6of17 American Institute of Aeronautics and Astronautics Paper 2004-2548

presence of a buffer region in reality between the viscous sublayer and the log-layer. However, this is much better than exclusively relying on no-slip condition. One can mimic the presence of the buffer layer by blending the linear and the logarithmic laws using an adequate blending function. The blending has some merits, inasmuch as, besides giving a smooth transition between the two layers (numerically more stable) and representing the mean velocity profile in the buffer layer more accurately, the blending also allows the respective laws in the two regions to be independently modified or extended to take into account other effects such as pressure gradient, surface roughness, and transpiration.

In the present work, we employed the blending function suggested by Kader.17 In the finite-volume discretization adopted in this work, the blended wall-law is employed to compute wall-shear, essentially diffusion flux at wall. To that end, the wall-law is applied to the parallel components of the resolved velocity at wall-adjacent cells to compute the friction-velocity (En')in an iterative manner.

The LDKEM needs a wall boundary condition for ks 1 . To that end, either a Dirichlet-type boundary condition (k 8g, = 0) or a Neumann boundary condition (ak3gl/an = 0, in view of k.g, - 2) is conceivable. In the present work, we simply set the diffusion flux of kegs at walls to zero, which is essestially equivalent to the Neumann boundary condition.

IV. Flow solver The present work was carried out using FLUENT, a general-purpose CFD code. FLUENT employs a cell-ccntcred finite-volume method based on a multi-dimensional linear reconstruction scheme, which permits use of computational elements (cells) with arbitrary polyhedral topology, including quadrilateral, hexahedral, triangular, tetrahedral, pyra-midal, prismatic, and hybrid meshes. There are several choices of the solver algorithms in FLUENT including coupled explicit, coupled implicit, and segregated solvers. For the computations presented in this paper, we used the segregated solver exclusively.

In the segregated solver, the governing equations are solvedie- UeniEal .7Several difterent solution algorithms are offered including SIMPLE, SIMPLEC, PISO, and fractional-step method (FSM). The temporal discretization in the segregated solver employs a fully-implicit, three-level second-order scheme. Time-accurate solutions can be obtained using either iterative time-advancement (ITA) scheme or non-iterative time-advancement (NITA) scheme. The NITA scheme greatly saves CPU time, since the costly outer iterations arc not needed. Unless stated otherwise, we used the fractional-step method in conjunction with the NITA scheme in this study.

Accurate spatial discretization is crucial in LES. The spatial discretization schemes employed in this work are based on a multi-dimensional linear reconstruction scheme. 18-20 Diffusive fluxes are discretized using central differ-encing. Discretization of convective fluxes requires caution in LES. Upwind-biased schemes such as second-order upwind, QUICK, and third-order MUSCL schemes have been most widely used for RANS computations. Unfortu-nately, numerical diffusion introduced by upwind schemes, which might be acceptable in RANS computations for high Reynolds-number flows, is detrimental to LES. This is because, in LES, numerical diffusion, however small it is, can easily overwhelm physical diffusion. This is the case even with high-order upwind schemes. For this reason, for LES, central-differencing schemes have been preferred for their meritoriously low - or zero in ideal conditions -

numerical diffusion. Thus we added a second-ordcr central diffecrncing (CD) scheme for discrctization of convective terms specifically for LES.2 1 Unfortunately, any pure CD schemes are susceptible to producing unphysical oscillations in the solution fields, which becomes especially pervasive in high Peclet-numbersituations - low diffusivity and coarse mesh which is almost a norm in LES for industrial applications. The usual remedy is to add a modicum of numerical dissipation, either explicitly or implicitly, to suppress the oscillations at the price of sacrificing spatial accuracy. In our implementation of CD, however, no numerical dissipation was explicitly added. And, unless stated otherwise, the CD scheme was used in this work. To back up the CD scheme in case it fails, we also developed what may be called a bounded central differencing (BCD) scheme that essentially detects in the solution fields any wiggles with a wave length of 2Ar or less (%< 2&r) and suppress them by switching to upwind schemes of varying orders depending on the severity of the wiggles, while retaining the CD elsewhere. It should be emphasized that the BCD scheme significantly differs from the often-employed hybrid schemes blending central differencing and upwind schemes with a fixed ratio.

The BCD scheme is reserved for industrial applications involving high-Reynolds number flows and less-than-ideal 7 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

meshes.

The discretized algebraic equations are solved using a point-wise Gauss-Seidel iterative algorithm. An algebraic multi-grid (AMG) method is employed to accelerate solution convergence. The solver is fully parallelized, which is crucial in LES for industrial applications.

V. Validations The LES capability described so far has been validated for a wide range of wall-bounded flows from simple to complex ones. In this paper, we present the results of a fully-developed channel flow and two bluff-body flows. The channel flow case is a fundamentally important case whose subtlety offers an opportunity to critically evaluate various aspects of SGS turbulence modeling such as the dynamic procedure used to determine the model constants.

The bluff-body flows include the one around a cylinder with square cross-section at a moderately high Reynolds number and the one around a sphere at two Reynolds numbers. Deliberately chosen, both involve large-scale, co-herent structures around the bodies and in the wake, representing typical bluff-body flows encountered in industrial applications.

A. Fully-developed channel flow at Rez = 180 A fully-developed channel flow was computed for the Reynolds number of Ret = 180 (ReH = 3,300) using the two dynamic models. The computational domain is a box of the size [27tH x 2H x itH] in the axial, normal, and spanwise directions, respectively. The computational domain is bounded by two walls on the bottom and the top of the channel, two pairs of periodic boundaries in the axial and the spanwise directions. The computations were carried out using two hexahedral grids; a coarse grid with 36 x 36 x 36 cells and a globally refined mesh with 72 x 72 x 72 cells. The resolutions of the meshes are such that, with the coarser mesh, y+ value at the wall-adjacent cells is approximately 0.6, and the cell size is Ay+ = 27 near the channel centerplane. The channel walls are treated effectively as no-slip boundaries due to sufficiently low y+ values at the wall-neighboring cells. On the pair of periodic boundaries in the axial direction, a pressure-drop across the pair derived using the given wall-shear (new = put) was specified, with the flow-rate determined as a part of the solution. The time-step size of At+ = 0.3 was used, where A+ = At u/V. The CD scheme was used for the discretization of convection terms.

The mean axial velocity (U+) and the three r.m.s.velocity components (it'+, V+, W+) predicted using the two dy-namic models are shown in Figure 1 on the following page along with the DNS results of Kim et aL.2 2 The predictions with the coarse mesh are seen to ovcrpredict U+ by about 8- 12% near the center of the channel. The peak in the us+

profile is also overpredicted with the coarse mesh, whereas the peaks in the profiles of V+ and w1+ are underpredicted.

Our results with the coarse mesh show largely the same trends as found by others who employed grids of somewhat finer resolutions (32 x 64 x 32 mesh, 23 65 x 65 x 65 meshl), yet closely matching their predictions despite the coarser mesh employed in this work. However, our predictions of the r.m.s. velocity components near the channel center is relatively poor. We surmise that the much larger grid spacing near the channel center (Ay+ = 27) is responsible for that.

The predictions improve greatly with the fine mesh. Both dynamic models reproduce the DNS results remarkably well. The mean axial velocity and the r.m.s. fluctuating velocity components are accurately predicted throughout the entire range of y+. Particularly noteworthy is the excellent agreement with the DNS data in terms of the peaks values and their locations of the r.m.s.fluctuating velocity components. Overall, the results obtained with the fine mesh compare favorably with other results mentioned earlier. 1 23 . For instance, our predictions are substantially closer to the DNS data than the results of Haworth and JansenI who computed the same channel flow using LES on a 65 x 65 x 65 node mesh using the Lagrangian dynamic Smagorinsky's model.

The present results are promising, inasmuch as they demonstrate that the second-order CD scheme in conjunction with the dynamic models is able to accurately predict this fundamentally important wall-bounded flow carrying an intricate near-wall physics. Regarding the impact of SGS modeling, we did not find any significant difference between the results from the two dynamic models. This should not come as a big surprise, however, since the fully-developed channel flow is near equilibrium in the mean.

8 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

v (a) Mean velocirv (U+) vredictions (b) r.m.s. velocity (us+I Dredictions t,

3:

y y (c) r.m.s. velocity (/+) predictions (d) r.m.s. velocity (R'+) predictions Figure 1. The results of LES using two dyanaic SGS turbulence models for Reg = 180 9of17 American Institute of Aeronautics and Astronautics Paper 2004-2548

-4.

B. Flow around a square cylinder 24 The flow past a square cylinder measured by Lyn et aL was considered. The Reynolds number based on the freestream velocity (Uo) and the width of the cylinder (H) is 22,000. The subject flow is featured by a massive flow separation accompanied by unsteady large-scale structures of widely varying length scales. As such, it aptly represents 25 turbulent flows around bluff bodies with sharp edges. Some others have also tackled this flow using LES.13.

3 The domain size and the mesh resolution were chosen in reference to the earlier studies by others. 25 A more comprehensive study using different mesh resolutions and domain sizes (spanwise in particular) is deterred Tor a luture study. Our objective here is to evaluate the efficacy of the present LES capability by comparing the predictions with other results based on meshes with comparable resolution. The computational domain is bounded by an upstream inlet boundary, top and bottom boundaries located 7.0 H from the center of the cylinder, and an exit boundary at 20 H from thc cylinder axis. Frecstream conditions were specified at the inlct. The top and the bottom boundaries were treated as symmetry planes (frictionless walls). The exit boundary was modeled as a pressure boundary where the solution variables are extrapolated in a mass-conserving manner. A pair of periodic boundaries separated with a span of 3.0H were used in the spanwise direction. The computational domain is filled with a hexahedral mesh with 660,000 elements. We took advantage of our unstructured mesh capability, embedding a block of locally refined mesh around the cylinder to better resolve the near-wall and wake regions as depicted in Figure 2 on the next page. The averaged wall-distance at the wall-adjacent cells is 0.012 H. The time-step size (Al) used for the present computations is 0.02 3

time unit (H/Uo) which is comparable to that used by others.1 The CD scheme was used for the discretization of convective terms. The statistics were obtained during the LES for a sufficiently long period of time, typically for more than several scores of time units.

Figure 2 on the following page shows the profiles of the mean axial velocity (U/Uo) and the r.m.s.velocity fluctua-tions (i' and a') along the centerplane (y = 0) in the wake predicted using the two dynamic models. The time-averaged axial velocity distributions show that the length of the recirculation bubble behind the cylinder is predicted by the two dynamic models to be around L, = 0.9 which agrees remarkably well with the cxperimental value (L, -- 0.9). The negative peak and the recovery of the mean velocity in the near-wake are also captured very closely. However, the predictions start to deviate from the measurement for x1H > 2.0, reaching 0.8 Uo asymptotically in the far-wake. As 24 shown in the figure, this value is considerably larger than what the measurement indicates (0.62Uo). Interestingly, t 3 25 have grossly overpredicted the recovery of the axial mean velocity. Our others who computed the same flow . also LES predictions of the mean axial velocity in the far-wake with both dynamic models were found to be largely com-3 parable to the prediction by Sohankar et aL.1 based on their dynamic one-equation model denoted by "OEDSMA' in their paper. However, the present LES predictions reproduce the mean velocity profile in the recirculation bubble more accurately than others. Particularly noteworthy is that our DSM yields somehow a much better prediction than the DSM model used by Sohankar et al. in terms of the recirculation bubble size and the asymptotic value of the axial mean velocity in the far-wakc, which begs a question of what could possibly contribute to this sizable difference.

One possible cause is the effectively finer mesh used in the present computations which was made possible by the embedded region of fine mesh around the cylinder. It is also quite likely that the differences in the details of the DSM implementation is responsible. In this regard, it should be noted that Sohankar et al. average Cv in the spanwise (homogeneous) direction, whereas the DSM used here does not.

The r.m.s.fluctuating velocity components are also predicted with a reasonable accuracy by the present LES. The peak values are appreciably underpredicted. However, the locations of the peaks are closely captured. Another obser-vation worthy of mentioning is that v' is relatively poorly predicted by both dynamic models.

Table I on page 12 summarizes other global parameters predicted by the present computations, along with the results predicted by others.13 .25 Our LES predictions of the mean drag coefficient, Strouhal number, and r.m.s.lift coefficient well match the measurements and the predictions by others.

10 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

1111E11111!3!

U1111 Tl ll 1 II I

(. wi (a) Mesh with an embedded line mesh (b) Axial mean velocity 02 0.4 -

LES with DSMI (prexent) 0.1 / l --- LES with LDICEM (prcscnt)_

Lyn et al. (1995) 0.2 -

1.B a

°.O.0I I 5.0 10.0 .0 2.0 4n x/H XJH (c) r~m.s. axial fluctuating velocity (d) r~m.s. verticvi1 fluctuating velocity Figure 2. Mesh used for the flow around a square cylinder and the mean axial velocity and normal stress distributions along the wake centerline II of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

Table 1. Summary of the LES predictions of the global quantities for the square-cylinder case (ReH = 22,000)

Methods L, St CD C'1 DSM (present) 0.9 0.133 2.19 1.19 DSM (Sohankar etaL.1 3 ) :0.6 0.126 2.03 1.23 DSM (Fureby et aL25 ) 0.83 0.132 2.0 1.34 LDKEM (present) 0.9 0.131 2.14 1.17 OEDSM (Sohankar etaL.3 ) 0.6 0.132 2.32 1.54 LDKEM (Fureby et al25 ) 0.74 0.130 2.10 1.32 Measured (Lyn et al.24) m0.9 0.13 st: 2.1 1.2 C. Flow around a sphere

1. DirectsitnuationforReD= 300 Before tackling the turbulent flow cases, laminar flow at ReD = 300 was computed on a hybrid unstructured mesh using direct numerical simulation. At this Reynolds number, the flow exhibits a weak unsteadiness leading to oscillations in drag and lift forces. The hybrid mesh has a total of 860,000 cells, consisting of prismatic cells in the near-wall region grown from the surface triangles on the body surface and tetrahedral cells filling the rest of the solution domain. The time-step size of 0.04 D/U was used. Several others have computed this flow to validate their numerics. 26-28 Table 2. Summary or the prediction ror the laminar flow over a sphere (Ren = 300)

Methods St CD Present 0.133 0.667, Tomboulides et aL' 0.136 0.671 Johnson nndPatel26 0 1327 .lO.65m Measured 2 9 0.15 -0.16 0.629 The results are summarized in Table 2 along with others' predictions. 2 62 8 Tomboulides et aL 27 and Johnson and Patel 2 6 used, respectively, a high-order spectral element method and a second-order upwind finite difference scheme on high-quality structured meshes. As shown, our predictions agree well with others' results, which is remarkable considering that a hybrid unstructured mesh was used in this work in conjunction with the second-order discretization scheme.

2. Iurbulentjtows WeV considered two Reynolds numbers (ReD = 1.0 x 104,1.14 x 106), one being in sub-critical and the other in super-critical regime. The subcritical flow case has been numerically studied by several others,2 8' 30 while the supereritical case was studied experimentally by Achenbach. 31 In this work, a hybrid unstructured mesh was deliberately used for both Re cases. The hybrid mesh has in total 2.46 million cells, consisting of 0.6 million prismatic cells in the near-wall region and tetrahedral cells filling the rest of the solution domain, with a large fraction of the total cell counts clustered in the near-wake region (see Figure 3).

The mesh quality is not exceptionally high and yet quite reasonable, except the rapid expansion of cell size around the 12 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

(a) Hybrid mesh consisting of prismatic and tetrahedral elements (b) Enlarged view of the near-wall mesh Figure 3. Ilybrid unstructured mesh used for the flow around a sphere cell-clustered region. The average distance from the wall at the wall-adjacent cells is around 1.1 x 10-3 D. For the lower-Re case, the near-wall mesh is sufficiently fine to resolve the boundary layer, with they+ values at the wall cells below y+ = 1.0 for most part of the wall. For the higher-Re case, however, the near-wall mesh is far from being fine enough to accurately resolve the boundary layer which is much thinner than the lower-Re case, and the y i values at the wall cells increase by almost two orders of magnitude. Thus, the wall adjacent cells are most likely to penetrate the fully turbulent region (log-layer) on a significant portion of the wall, especially near 0 = 90° where the skin-friction reaches a maximum. The mesh being not ideal, the higher-Re case offers a good opportunity to assess the wall-function based approach adopted in this work. Partial views of the mesh are shown in Figure 3.

Attempts to use the pure CD scheme have not been successful for this case. Numerical oscillations were observed sporadically in a few spots rather remote from the body where the cell size increases rapidly, being accompanied by abnormally large velocity magnitude. Although the oscillations were not catastrophic and affected the global quantities very little, the subsequent computations were carried out using the BCD scheme discussed earlier. A time-step size of 0.02 D/U was used for both Reynolds numbers. The data were collected for more than hundreds of time units.

Figure 4 shows the time-histories of CD for the lower-Re case recorded during the LES using the two dynamic models.

Table 3. Summary or the LES prediction for turbulent flow past a sphere (RCD = 10,000)

Methods CD St 4S Of LES with DSM (present) 0.438 0.182 86- 87 86- 88 LES with LDKEM (present) 0.433 0.185 86- 87 86- 88 DES (Pclacz et aL)30 0.430 - _ _

DES (Constantinescu et aL.)2 8 0.397 0.200 84 - 87 93 - 108 LES (Constantinescu et al.)2 8 0.393 0.195 84 - 86 86 - 88 Measured 31 .3 2 0.40 0.185 - 0.19 Correlation 33 0.46 -

13 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

The global quantities predicted for the lower Re case are summarized in Table 3 on the page before along with others' results. The predicted CD values (0.438 and 0.433 for the DSM and LDKEM, respectively) are in fair agreement with the often-quoted experimental value of 0.4 measured in 1920's 3 4 and other predictions. Constantinescu et aL28 predicted CD at around 0.4 using LES and detached eddy simulation (DES) on a structured hexahedral mesh having 450,000 nodes. Our predictions are closer to the value obtained by Pelaez et al.30 (CD FL0.43) who carried out a DES on an unstructured mesh with 770,000 nodes. The Strouhal numbers predicted by the two dynamic models came out very close to each other, matching the measured one quite closely in view of the scatter in the experimental data.

The data of Achenbach3 ' and Kim and Durbin35 favor lower Strouhal number around St = 0.15, whereas Sakamoto's data 32 suggests a higher value between 0.18 and 0.19. The locations of flow separation (45) predicted by the two dynamic models are nearly identical and were found in 86° - 870. The locations of laminar-to-turbulent transition (0,), which were obtained in this work by reading off the angle beyond which v, or kegs increases rapidly, were found in 860 88° for both dynamic models, which compare well with the LES predictions by Constantinescu et aL28 0.5

_ __ _ __ __ _ _ _ _ _ _ _ _ _ _ _ _ l ~

lR0 _ _ _ _1-0.50 I 10 LU LU.

m 2Experimental _ , LoESwith DSM _X 0.5 jJami10 ac j , .&1l]l . 1.fit I _ II l I i @I fl IM II 11 1 I 1et J 045 - SA

,.l!

e10 10"_ __ .- as P__ ___ .=

10 . Il = :: lii- - miiL

-,_ = :i 0 40 - 0IM _- " _ A 21

.2 _ 1.ii

]lll I EaS I lt 0.3~0J, I .

to 200 300 400 500 600 102 io leo t Uu/D ReD (a) Time histories of CD for Re 0 = 10,000 (b) Comparison of the predicted time-averaged CD values (ReD =

300,10,000,1.14 x 106) and the experimental mean CD-curve for a range of Reynolds number Figure 4. 'itme histories of drag coefficient (CD) and and the predicted mean drag coefticients for the sphere With the higher-Re case, the flow in reality has already undergone the "drag crisis", and the flow structure has changed drastically from those of subcritical regime. The change in the flow structure can be seen from Figure 5 and Figure 6. Depicted in these figures are the iso-contours of the second-invariant of the deformation tensor, (niQjfl -

SjSij)/2. Both figures aptly portray the hairpin-like vortical structures in the wake observed in experiments. The wake for the higher-Re case (Figure 6) is much narrower than in Figure 5, which is the consequence of the delayed onset of flow separation. The CD values predicted by the DSM and the LDKEM are 0.139 and 0.142, respectively. These values are fairly close to the range of values (CD = 0.12 0.14) measured by Achenbach. 3 '

For the higher-Re case, the predictions of the locations of the separation and the onset of transition were much less satisfactory. At Re = 1.14 x 106, the experimental results 31 show that the transition occurs near 970 98°well before the separation occurring near 1200. The present predictions failed to reproduce this experimental finding. The present results exhibit too early a separation at around 1000 and a delayed transition. This discrepancy is most likely due to the use of too coarse a mesh in this work to resolve the very thin boundary layer for the high-Re case.

The drag coefficients predicted in this study for the three Reynolds numbers are plotted in Figure 4 along with the mean experimental curve.

14 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

Figure 5. Vortical structure in the near-wake Ot the sphere for ReD = 10,000 - visualized using the iso-contour or the second-invariant Or the velocity deformation tensor Figure 6. Vortical structure In the near-wake of the sphere for ReD = 1.14 x 106 - visualized using the Iso-contour or the second-Invariant of the velocity deformation tensor 15 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

VI. Summary and Conclusion A large eddy simulation capability based on a finite-volume solver has been developed and validated for a num-ber of wall-bounded flows. The finite-volume solver employs second-order numerics and permits use of unstructured meshes, thus being able to easily handle industrial applications involving complex geometry. Turbulence closure for subgrid-sale stresses is effected using two dynamic subgrid-scale viscosity models, namely, the dynamic Smagorin-sky model (DSM) and the localized dynamic k-equation model (LDKEM). These two dynamic models allow one to compute arbitrary three-dimensional flows without any statistically homogeneous directions. The validations demon-strated that the present LES capability is capable of predicting the wall-bounded flows of varying complexity with a commendable accuracy, having a potential to provide a practical tool for high-level simulation of turbulent flows encountered in industrial applications.

VII. Acknowledgments The author gratefully acknowledges use of Fluent Inc.'s software and thanks the members of the development group at Fluent Inc. Special thanks go to Sunil Vijay Unaune at the Fluent India who contributed the mesh used for the computation of the flow past a sphere.

References

'llaworth. D. C. and Jansen. K., "Large-FAdy Simulation on Unstructured Deforming Meshes: Towards reciprocating IC Engines:' Cotnputers Fluids. Vol. 29. 2000. pp. 493-524.

2 K. Mahesh. G. Constantinescu. P. M.."Large Eddy Simulation of GasTurbine Combustors, AnnualReserach Briefs, CenterforTurbulence Research, NASA Ames/Stanford University, 2000. pp. 219-229.

3 Urbin, G. and Knight, D., 'Large Eddy Simulation of a Supersonic Boundary Layer Using an Unstructured Grid:' AIAA Journal. Vol. 39, No. 7, July 2001, pp. 927-935.

4 K. Jansen. T. NM.and Reba. R., "Finite-Element Based Large-Eddy Simulation of the Near-Nozzle Region of a Compressible Round Jet,"

Aiaa paper 2002-2358, 2002.

"S. Benhamadouche. K. Mi.and Constantinescu, G.. 'Colocated Finite-Volume Schemes for Large Eddy Simulation on Unstructured Grids,"

Proceedings of the Summer Program. Center for Turbulence Research. NASA Aames/Stanford University. 2002, pp. 143-154.

6 M. Germano, U. Piomelli. P. M. and Cabot. W. 11.. "Dynamic Subgrid Scale Eddy Viscosity Model." Physics of Fluids A. Vol. 3. 1991.

pp. 1760-1765.

7 Lilly. D. K., 'A Proposed Modification of the Gennano Subgrid Scale Closure Method." Physics of Fluids A, Vol. 4. No. 3, 1992. pp. 633-635.

8Kim. W.-W. and Menon. S. 'A New Dynamic One-Equation Sub-Grid Scale Model for Large Eddy Simulations." Aiaa paper 95-0356.

1995.

9 Kim, W.-W. and Menon, S., "Application of the Localized Dynamic Subgrid-Scale Model to Turbulent Wall-Bounded Flows, Aiaa paper 97-0210, 1997.

tfSinagorinsky, J., "General Circulation Experiments with the Primitive Equations, part I: The Basic Experiment," Monthly Wecather Reiview.

Vol. 91, 1963. pp.99-165.

1IS. Ghosal, T. Lund. P. M. and Akselvoll. K., "A Dynamic Localization Model for Large Eddy Simulation of Turbulent Flows," Journalof Fluid Mechanics, Vol. 286, 1995, pp. 229-255.

12 Davidson. L.. 'Large Eddy'Simulation: A Dynamic One-Equatio Subgrid Model for Three Dimensional Recirculating Flow'"Proceedings of the 11th International Symposium Turbulent ShearFlows, Vol. 3,1997, pp. 26.1-26.6.

13 A. Sohankar. L. D. and Norberg. C., "Large Eddy Simulation of Flow Past a Square Cylinder. Comparison of Different Subgrid Scale Models," Journal of Fluids Engineering, Vol. 122. No. 1.2000. pp. 376-404.

t4 Krajnovic. S. and Davidson. L., "Large Eddy Simulation of the Flow Around a Bluff Body." AIAA Journal, Vol. 40, No. 5, May 2002, pp. 927-935.

15S. Liu, C. M. and Katz, J., "On the Properties of Similarity SUbgrod-Scale Models as Deduced from Measurements in a Turbulent Jet."

Journal of Fluid Mechanics, Vol. 275, 1994, pp.83-119.

t6S. Liu. C. MI. and Katz, J., Experimental Study of Similarity Subgrid-Scale Models of Turbulence in the Far-Field of a Jet, Direct and Large Eddy Simulation, Kluwer, 1994.

17 Kader, B. A., "Temperature and Concentration Profiles in Fully Turbulent Boundary Layers." tnt. J. Heat Mass Transfer. Vol. 24, No. 9.

1981. pp. 1541-1544.

"tMathur.S. R. and Munhy.J. Y.,"A Prcssure-Based Method for Unstructured Nieshes." NumericallHeatTransfer. Vol.31.1997. pp. 195-215.

16 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548

' 9 S.-E. Kim. S. R. Niathur. J. Y. M. and Choudhury, D., "A Reynolds-Averaged Navier-Stokes Solver Using Unstructured Mesh-Based Finite-Volume Scheme" Aiaa paper 98-0231. 1998.

"S.-E. Kim, B. K. and Caraeni, D.. "A Nlultidimensional Linear Reconstruction Scheme for Arbitrary Unstructured Grids. Aiaa paper 2003-3990. Jan. 2003.

21 L. Davidson. D. Cokljat, J. F. M. A. L C M. NV.R.. 'LESFOIL: Large Eddy Simulation of Flow Around a High Lift Airfoil," Notes on Numnerical Fluid Mechanics and Mudidisciplinary Design - Volune 83. Springer. 2002.

22J. Kim, P.M. and Moser. R.. "Turbulent Statistics in Fully Developed Channel Flow at Low Reynolds Number" Journal of Fluid Mechanics.

Vol. 177. 1987. pp. 133-166.

2D. Choi. D. Prasad. M. W. and Pierce, C., "Evaluation of an Industrial CFD code for LES Applications," Proceedings of the 2000 Sunmner Prograin. Center for Turbulence Research, NASA Aames/Stanford University, 2000, pp.221-228.

24 D. A. Lyn, S. Einav, W. R. and Park, J. H., "A Laser Doppler Velocifimetry Study of Ensemble Averaged Characteristics of the Turbulent Near-Wake of a Square Cylinder,' Journal of Fluid Mechanics. Vol. 304. 1995. pp. 285-319.

2'C. Fureby, G. Tabor. H. G. W. and Gosman, A. D., "Large Eddy Simulation of the Flow Around a Square Prism:' AMA4 Journal. Vol. 38, No. 3, May 2000. pp. 442-452.

26 johnson. T. A. and Patel, V. C.. "Flow Past a Sphere up to a Reynolds Number of 300,- Journal of Fluid Mechanics. Vol. 378, 1999.

pp. 19-70.

2A. G. Tomboulides, S. A. Orszag. G. A. K., "Direct and Large-Eddy Simulations of Axisymmetric Wakes," Aiaa paper 93-0546, Jan. 1993.

28G. Constantinescu. M. C. and Squires. K.. "Turbulence Modeling Applied to Flow over a Sphere.' AIAA Journal, Vol. 41, No. 9, 2003, pp. 1733-1742.

2')Roos, F.W. and Williams, W. %V.. "Some Experimental Results on Sphere and Disk Drag:'AIAA Journal, Vol. 9, No.2. 1971, pp. 285-291.

30J. Pelaez. D. J. M. and Kandil. O.. "Unsteady Analysis of Separated Aerodynamic Flows Using an Unstructured NMultigrid Algorithm." Aiaa paper 201TU6Jia-n200V.

3 tAchenbach. E., 'Experiments on the Flow past Spheres at Very High Reynolds Numbers." Journal of Fluid Mechanics. Vol. 54. No. 3.1972.

pp. 565-575.

3-Sakamoto. H. and Haniu. H.. "A Study on Vortex Shedding from Spheres in an Unifomi Flow,' Journal of Fluids Engineering, Vol. 112, 1990. pp.386-393.

33 White, F. NI.. Viscous Fluid Flow. NMcGraw-Hill. 1974.

34Schlichting. H.. Boundary-Layer Theory, McGraw-Hill, 1979.

35 Kim, K. J. and Durbin, P. A.. "Observation of the Frequencies in a Sphere Wake and Drag Increase by Acoustic Excitation." Physics of Fluids. Vol. 31. No. 11. 19S8, pp. 3260-3265.

17 of 17 American Institute of Aeronautics and Astronautics Paper 2004-2548