Text

E-mail Gucwa, Entergy Nuclear Northeast, to Ennis, NRR, VY EPU Supplement 30
	ML052220086
Person / Time
Site:	Vermont Yankee File:NorthStar Vermont Yankee icon.png
Issue date:	08/02/2005
From:	Gucwa L; Entergy Nuclear Northeast
To:	Richard Ennis; Office of Nuclear Reactor Regulation
References
	Download: ML052220086 (50)
	v • d • e

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 11 I Rick Ennis - BVYO5-052 Exhibits SRXB-A-18-1, EMEB-B-138-1, 2, and 3 Paei 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:

CC:

<Douglas.Rosinski~pillsburylaw.com>

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

Subject:

Creation Date:

From:

Created By:

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

LGUCW90@entergy.com Recipients nrc.gov owf4po.OWFNDO RXE (Rick Ennis) pillsburylaw.com Douglas.Rosinski CC Post Office ovvf4_po.OWFNDO Files Size MESSAGE 232 TEXT.htm 2431 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 Mime.822 5017391 Route nrc.gov pillsburylaw.com Date & Time 08/02/05 02:14PM 600742 832848 972271 1256329 Options Expiration Date:

Priority:

Reply Requested:

Return Notification:

Concealed

Subject:

Security:

None Standard No None No 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

.throughout its 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..5.2 2.8.5.2 (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; Decrease in Heat Removal by the Secondary System

.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 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 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 In an increase In reactor 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 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 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 In Matrix 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 In the 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 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 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 In either 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 e l (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 In SRP 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 In core 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

(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:

1 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 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.

L.

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:

24th 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 Fluent Incorporated, Lebanon, N.H.

03766, U.S.A.

ABSTRACT Large eddy simulations were carried out for the flow around a hydrodynamically smooth, fixed circular cylinder at two Reynolds numbers, one at a subcritical Reynolds number (Re = 1.4 x 105) and the other at a supercritical Reynolds number (Re= 1.0 x 106). The computations were made using a parallelized finite-volume Navier-Stokes solver based on a multidimensional linear reconstruction scheme that allows use of unstructured meshes. Central differencing was used for discretization of both convection and diffusion terms. Time-advancement scheme, based on an implicit, non-iterative fractional-step method, was adopted in conjunction with a three-level, backward second-order temporal discretization.

Subgrid-scale turbulent viscosity was modeled by a dynamic Smagorinsky model adapted to arbitrary unstructured meshes with the aid of a test-filter applicable to arbitrary unstructured meshes. The present LES results closely reproduced the flow features observed in experiments at both Reynolds numbers.

The time-averaged mean drag coefficient, root-mean-square force coefficients and the frequency content of fluctuating forces (vortex-shedding frequency) are predicted with a commendable accuracy.

INTRODUCTION Unsteady loading on a circular cylinder caused by the flow is a precursor to its vortex-induced vibration (VIV). The major difficulty of computing flows around circular cylinders originate from the fact that circular cylinder flows of practical interest are high Reynolds number (Re) flows, often involving a laminar-to-turbulent transition in various regions of the flows such as boundary layer, separated shear-layer, or near-wake.

The location of transition and the extent of laminar or turbulent flow regime, in real flows, depend on such factors as Reynolds number, freestream turbulence, surface roughness, span-diameter ratio (LD), and blockage, among others. The state of the flow in those regions dictates the formiation and evolution L. Srinivasa Mohan Fluent India Pune, India of large-scale turbulent structure around a circular cylinder, which in turn affects unsteady loading on the cylinder.

Apparently, turbulence modeling is an issue here.

There are largely three approaches being explored by CFD practitioners to modeling high-Re turbulent flows around circular cylinders and bluff bodies in general.

One approach employs unsteady Reynolds-averaged Navier-Stokes (URANS) equations supplemented by turbulence models as the governing equations.

URANS-based

approach, despite its low computational cost mainly due to less demanding mesh resolution requirement, is fundamentally ill equipped to capture large-scale turbulent structure present in the flows.

Results obtained using URANS computations typically underpredict the amplitudes of fluctuating forces.

Obviously, this shortcoming has a serious negative implication to accurate prediction of VIV.

Fundamentally more viable than URANS-based approach for bluff-body flows is large eddy simulation (LES). LES is designed to directly resolve large eddies, with the effects of unresolved smaller eddies, namely subgrid-scale turbulence, on the resolved large eddies accounted for by turbulence models.

However, LES is computationally very expensive, often prohibitively, especially when thin turbulent boundary layer is a predominant feature of the flow in question to be accurately resolved in a LES.

Resolving near-wall turbulence with infinitesimal length and time scales requires extremely fine mesh and small time-step size. For that matter, LES is more feasible for free flows such as jets, mixing layer, and wakes, and massively separated flows.

There are very few LES studies published in the literature of flows around circular cylinders at high Reynolds numbers.

Breuer [I] was the one of the very few who tackled the case of a high subcritical Reynolds number, Re = 1.4 x 105, at which the experiment of Cantwell and Coles [2] was conducted.

Using a multi-block, structured-mesh-based finite-volume solver with an explicit time-marching scheme, Breuer was able to predict the global parameters of the flow, and the mean flow and turbulence with a commendable accuracy. More recently, I

Catalano et al. 13] attempted LES for three Reynolds numbers in critical to supercritical Reynolds regime (Re = 0.5 x 106, 1.0 x 106, 2.0 x 106).

Their predictions of the global flow parameters were in reasonable agreement with the experimental data, although the skin-friction and the Reynolds number-dependency of the mean drag coefficient were poorly predicted.

A more recent trend in turbulence modeling of bluff-body flows is to employ what may be called "hybrid" approaches in which one attempts to adjust turbulence models to local mesh resolution in one way or another.

In one such approach referred to, in the literature, as detached eddy simulation (DES), either a RANS-based or a subgrid-scale (SGS) turbulence model is invoked depending on whether or not the filter-length (local mesh size) is larger than the local integral length-scale of turbulence.

In DES, turbulence models essentially reduce to RANS models in near-wall region or when the local mesh size is too coarse to explicitly resolve energy-containing eddies. One fundamental criticism about DES comes down to the lingering question of how one can possibly reconcile a RANS turbulence model and a subgrid-scale turbulence model, two very different models in terms of their spectral content, at the common interface.

It should also be realized that, as a consequence of falling back to a RANS model in near-wall region, the fidelity of DES would be solely determined by that of the embedded RANS model. Quite a few studies employing DES have appeared recently on circular cylinder flows. Among others, Travin et al.[4], Vatsa and Singer[5], and more recently Pandya et al. [6] all employed DES based on an eddy-viscosity transport model to study the flow around a smooth circular cylinder at subcritical and supercritical Reynolds numbers. The results of these DES are mixed, insofar as some of the global flow parameters such as the r.m.s. force coefficients and the mean flow in the near-wake were predicted poorly.

The objective of this study is to assess the feasibility of LES for high Reynolds number flows around a fixed, smooth 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 1 06) regime were deliberately chosen so that they straddle the critical Reynolds number (- 3 x 105). The main concern is how well LES can reproduce the known characteristics of the flow at the two Reynolds numbers and the changes in the global flow parameters such as drag coefficient, r.m.s.

force coefficients, and vortex-shedding frequency as the Reynolds number increases.

The results of the LES will be compared whenever possible, with the experimental data and other LES and DES results reported in the literature.

NUMERICAL METHOD The computations were carried out using the segregated solver in FLUENT, a general-purpose CFD software.

The details of the numerical method are described in the following.

SPATIAL DISCRETIZATION SCIIEMIES triangular, tetrahedral, pyramidal, prismatic, and hybrid meshes.

The solution gradients at cell centers that are needed to compute convective and diffusive fluxes are determined applying the Green-Gauss theorem [7]. Diffusive fluxes are discretized using central differencing.

Discretization of convective fluxes requires great caution in LES.

Upwind-biased schemes such as second-order upwind scheme have been widely used for RANS computations. Unfortunately, numerical diffusion introduced by upwind schemes, which might be acceptable in RANS computations for high-Re flows where eddy-viscosity is larger than molecular viscosity by orders of magnitude, can overwhelm physical diffusion that is typically much smaller in LES; subgrid-scale turbulent viscosity is smaller than RANS-based eddy viscosity by orders of magnitude. For this reason, a central differencing scheme [13]

was employed for the discretization of convection terms in this study.

TIM E-ADVANCEIM ENT SCIIEMIES An implicit fractional-step method (FSM)

[8]

in combination 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 equation using an approximate factorization of the coupled Navier-Stokes equations. For incompressible flows, the FSM preserves the formal second-order temporal accuracy without having to perform, at each time-step, costly outer iterations to couple velocity and pressure. The FSM thus provides a highly efficient algorithm for CPU-intensive transient computations like LES.

Consider a semi-discrete forn of the Navier-Stokes equations in "pressure-correction" (0.41 = p" ' - p") form,

[

A G](up' )=Or)

D 0 j

47) t (I) where u"'} and r are the velocity vector and the momentum source vector, respectively, and the integer n the time level, with n+J refers to the current time level. A is the coefficient matrix defined in terms of the discrete convective and diffusive operators and the time-advancement scheme chosen, and G and D the discrete gradient and divergence operators, respectively.

The coupled system of equations in Equation (I) is extremely difficult to solve as it stands, since the matrix A has to be inverted for every iteration.

In the FSM, the original coupled equations in Equation (I) are approximated by A

10 1

== r)

[~o(A,2)]

ID A0DG][iO I ]

(r'nj 2

to (2)

Factorizing equation (2), we have a series of split operations that consist of FLUENT employs a cell-centered finite-volume method based on a multidimensional linear reconstruction scheme that permits use of computational elements (cells) with arbitrary polyhedral

topology, including quadrilateral, hexahedral, A6 =-Gp' +r" AtDG4V = Dfi u

= u -AtG47"'

(3) 2 Copyright C 2005 by ASME

On a per-iteration basis, the series of operations in 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 solutions converge, while the FSM needs only one sweep. To preserve second-order accuracy with the 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 individual equations and high-order source terms. Yet, the non-iterative FSM is takes much less CPU time than the iterative SIMPLEC scheme.

The system of discretized governing equations are solved using a point-implicit, Gauss-Seidel relaxation along with an algebraic multigrid (AMG) method to accelerate solution convergence.

The solver and the subgrid-scale turbulence model are fully parallelized.

As a validation for the spatial discretization schemes and the transient algorithm, we computed laminar flow around a circular cylinder at the Reynolds number of 100 with both non-iterative FSM and iterative SIMPLEC scheme.

At this Reynolds number, the flow exhibits a periodic vortex-shedding.

A C-type structured mesh with 48,000 cells was used for the computation. A dimensionless time-step size of At = 0.04 (non-dimensionalized by D/U& D being the cylinder diameter, U0 the freestream velocity) was used, with which one period of the vortex-shedding was resolved with approximately 150 time-steps.

The predicted mean drag coefficient (ED), r.m.s. lift coefficient (ci), and Strouhal number (St) are summarized in Table 1. The FSM gives practically the same predictions as the iterative SIMPLEC scheme, and the global parameters predicted by both methods are seen to closely match the other predictions and the experimental data.

Table 1. Summary or the global flow parameters predicted for laminar flow (Re = 100) - Norberg 110] compiled numerical results whose references can be found in ref. 110]

(5,, _ =-C.Y _

(5)

We determine the model constant, Cx, using the dynamic Smagorinsky model (DSM) originally proposed by Germano et al. [12]. To implement the dynamic procedure for the present finite-volume solver requires a special test-filter applicable to arbitrary unstructured meshes. The dynamic procedure employed in the present study is "local" in the sense that it does not require the existence of any statistically homogeneous directions, being applicable to arbitrary complex three-dimensional flows. The details of the implementation of the DSM in the framework of unstructured mesh based finite-volume solver can be found in Kim [13]. The DSM has been validated for a number of wall-bounded Ilows such as a square cylinder and a sphere.

SOLUTION DOMAIN AND COMPUTATIONAL MESH Our choices of the solution domain and the computational mesh were guided by the findings in Breuer's LES study [1].

What one should keep in mind in determining the size of the solution domain is that the spanwise extent of the domain should be larger than the spanwise correlation length of turbulence.

Fortunaiely, the spanwise correlation length is known to decrease as the Reynolds number increases. We took advantage of this fact, having decided to use a spanwise length of 2.OD for both the subcritical and the critical Re numbers.

The top and bottom boundaries of the domain were placed at 10.5D from the cylinder axis. Thus, the blockage ratio of our numerical model (HID, where H is the height of the domain) is approximately 4.8 %. The upstream inlet and the downstream 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 million cells was used for the computations for both Reynolds numbers. The overall mesh resolution is comparable to case

-3BI" in Breuer [1]. The near-wall mesh resolution is such that the distance from the cylinder surface is 104D at the wall-adjacent cells, which translates to y+(=ury/l well under 1.0 for the subcritical flow, and below 6.0 for the supercritical flow.

For the subcritical flow, the boundary layer was resolved with 10 - 18 cells in the laminar region.

E.__

C C i.

St Present - FSM 1.336 0.28 0.165 Present - SIMPLEC 1.345 0.28 0.166 Park et al. 19]

1.33 0.33 0.165 Norberg I10]

0.23/0.29 Williamson [ ]

0.164 SUBGRID-SCALE TURBULENCE MODELING For incompressible flows, the filtered Navier-Stokes equations can be WTitten as au~aiu=_I a _

ha-ri+ a IVawi 4

a. axi

pax, ax, axeJ axJ where r., = uu, -

jF1, is the subgrid-scale turbulent stress. In this study, the subgrid-scale (SGS) stress is modeled using isotropic eddy-viscosity as BOUNDARY CONDITIONS AND OTHER DETAILS OF COMPUTATION On the cylinder wall, we employ a law-of-the-wall that invokes proper wall-laws depending on the mesh resolution, namely, ye at the wall adjacent cells. The top and bottom wall boundaries were treated as a slip (zero-stress) wall.

A 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 downstream exit boundary, we extrapolated the solution variables from the adjacent interior cells in a mass-conserving manner.

The time-step size (dz) used for the computations is 0.005 in a dimensionless unit for both Reynolds numbers. The time-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 0.214 From these rough estimates and with the dimensionless time-step size of 0.005, one turnover time of the smallest resolved eddies will be resolved with 50 time steps.

To shorten the initial transient period of the solution and to quickly attain a statistically stationary state, we took, as the initial condition, a partially converged steady RANS solution with pseudo-random fluctuating velocity components superimposed on the mean velocity field.

The computations were carried out on a 16-CPU lnx86 cluster with AMD Opteron processors and Infiniband interconnects. With the non-iterative FSM, the computation took approximately 12 CPU-seconds per a time-step, which translates into being able to complete roughly 7,200 time-steps in 24 hours2.777778e-4 days <br />0.00667 hours <br />3.968254e-5 weeks <br />9.132e-6 months <br />.

RESULTS Table 2. Summary of the global now pararneters predicted by the present LES for Re = 1.4 x 10, compared wilti other numerical results (Breuer [1]; DES-TS - Tram'in et td.14]; DES-LS

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

Colesl2], WVA - West & Apelt 114], SB: Szepessy & Bearmin [I15],

NO - Norberg 110], ZDM Zdravkovich [16])

lCD Ca

-C l

S i LES (Present) 1.13 0.59 1.20 0.67 0.205 LES (Breuer) 1.45 1.76 0.34 0.204 DES-TS 0.59 0.06 0.67 1.2 0.31 DES-LS 1.08 0.29 1.04 1.1 0.21 Exp.(CC) 1.24 1.21 0.5 0.18 Exp.(WA) 1.3 0.58 Exp.(SB) 1.35 0.5 1.05 Exp.(NO) 0.52 0.18/0.195 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.

Although the near-wake flow appears chaotic, the figure clearly shows the existence of a large-scale motion, namely alternate vortex-shedding known to occur in the subcritical Re number range. The color map of the resolved turbulent kinetic energy in Figure 1 also indicates that the boundary layer remains laminar up to the separation point, and transition occurs in the separated shear layer.

The present LES gave the mean CD around 1.13, which falls within the scatter of the experimental data (1.1 - 1.35) for smooth cylinders.

Interestingly, our prediction came out substantially lower than Breuer's prediction despite the comparable mesh resolution, the same formal order of spatial accuracy (second-order) and the same dynamic SGS turbulence model used in both computations.

Figure 1. The flow structure behind the circular cylinder at Re = 1.4 x ]05 -

iso-surface of the second invariant of velocity deformation tensor, colored by the resolved turbulent kinetic energy The global flow parameters predicted by the present LES are summarized in Table 2 along with the LES prediction of Breuer [1] and the DES predictions of Travin et al. [4], and the experimental data.

As for Breuer's LES result, among the cases involving different mesh resolution and SGS turbulence models, we picked the result of case "Bl" whose mesh and SGS turbulence model (dynamic Smagorinsky mode) closely match ours.

In selecting the experimental data for the comparison, we gave more weight to the ones measured on smooth cylinder surface, and with sufficiently large span-wise lengths (LID > 5) and smaller blockage-ratio (HID < 0. 1).

400 450 500 Non-dimensional time, tg (= tUJ/D)

Figure 2. Time histories of the drag and lift coefficients predicted by the present LES for Re = 1.4 x 105 The DES-TS result [4] is shown to severely underpredict the mean CD, better matching the data obtained with rough cylinders 1171, which is consistent with the DES predictions by others [5, 6].

It has been found that DES, run in a normal mode with finite freestream turbulence, essentially yields a tripped boundary layer, leading to a premature transition to turbulence and an early drag crisis even at a subcritical Reynolds number.

The DES-LS result [4] based on the so-called "tripless" approach, in which laminar flow is essentially enforced in the boundary layer, better predicts the mean CD.

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 St Figure 3. Power spectral density of CL time-history for Re

=1.4 x 10 The r.m.s. lift coefficient predicted by the present LES is within the scatter of the experimental data, although it is closer to the upper bound.

The CL history Figure 2 clearly shows that, despite modulation of the amplitude, there is a distinct frequency of vortex-shedding, which is consistent with what has been known for smooth cylinders in the subcritical regime.

The Strouhal number corresponding to the main shedding frequency was found to be 0.205 as shown in Figure 3.

2.0

-Present 1.0 o Exp. (Cantwell & Coles.1983)

Exp. (Sarioglu & Yavuz. 2002) i3 0.0 -

.1.0 S

-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 surface is shown in Figure 4, along with the experimental data 12, 18] obtained at the same Reynolds number. The negative peak predicted by the present LES is larger than the measured ones. However, the pressure level in the separated region - the plateau after 9= 90 degrees, and the mean base pressure (CPB) are all accurately predicted by the present LES.

x/D Figure 5. Time-averaged axial velocity profile in the wake predicted by the present LES for Re = 1.4 X 11)5 Figure 5 depicts the predicted time-averaged axial velocity profile along the wake centerline at the mid-span. The overall agreement between the prediction and the measurement is quite good. The rapid relaxation of the mean axial velocity in x/D =

1.0 - 3.0 shown in the measurement is an indication of vigorous mixing of momentum taking place in the near-wake region. That our prediction closely captures the recovery of the mean axial velocity verifies that the dynamics of the large-scale structure in the near-wake is well predicted in the present LES.

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 LES prediction is much closer to the data than the DES predictions. Considering that hot-wire measurements become increasingly difficult and less reliable near the recirculation region, our LES prediction is considered quite good.

5.0 _U Lower surface

-5.0 I

I I

I 0

30 60 90 120 150 ISO 0 (deg.)

Figure 6.

0.00 i

E.

I I

i I

I O.OW~t l--Upper surrace 0.004e l_% Lower surfacel 0.003 0.002 0.001, I

I peaks at 3.0 around 0= 50 degrees, which is closely reproduced by the present LES.

Supercritical Reynolds number Case (Re = 1.0 x10')

0 30 60 90 0 kdeg.)

120 ISO 180 Figure 7. Resolved turbulent kinetic energy distribution at 0.0001D above the cylinder surface, predicted the present LES at Re = 1.4 x 105 As mentioned earlier, the experimental evidence [19]

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 closely LES can reproduce the physics associated with the separation and the transition.

Figures 6 and 7 depict the circumferential variations of two quantities along the cylinder surface that shed some light on how the LES predicts the separation and transition.

Figure 6 shows a plot of the instantaneous skin-friction coefficient - rescaled to match the non-dimensionalization adopted by Achenbach [19] - against the circumferential angle.

First, it is noticeable that the skin-friction distribution exhibits a considerable top-bottom asymmetry.

This is a clearly the (upstream) influence of the large-scale, alternate vortex-shedding.

What is most interesting in the figure is the appearance of small transient separation bubbles on both upper and lower surface starting as early as. 70 - 75 degrees.

Interestingly, this range largely coincides with the range of separation angles reported by many others 12, 19] deduced from the inflection point of mean Cp curve; 77 degree by Cantwell and Coles [2], and 78 degrees by Achenbach [19] at Re = 1.0 x 105. Although not shown here, the separation angle determined 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, which was also found by Breuer [I].

Based on the instantaneous skin-friction distribution, we surmise that a transient laminar separation occurs much earlier than the mean separation angle.

The distribution of the resolved turbulent kinetic energy depicted in Figure 7 supports the possibility of an intermittent laminar separation around 70 - 75 degrees, insofar as the turbulent kinetic energy is still quite low there. The complete transition occurs a little downstream as indicated by the rapid increase in the resolved turbulent kinetic energy between 75 and 85 degrees. The magnitude of the skin-friction coefficient predicted by the present LES is also in good agreement with Achenbach's data

[191 measured at Re = 1.0 x

]0 5.

Achenbach's data show that the scaled skin-friction coefficient 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 turbulent kinetic energy Figure 8 portrays an impression of the instantaneous flow structure.

Compared to the flow structure for the subcritical flow shown in Figure 1, the wake is much narrower and more chaotic without any trace of a large-scale alternate vortex shedding, suggesting a delayed turbulent separation. This is typical of the flow in the supercritical regime. The color map on the iso-suraface - representing resolved turbulent kinetic energy - also indicates that the laminar boundary layer is sustained farther downstream than the subcritical flow.

Table 3. Summary of the global quantities predicted by LES for Re = 1.0 x 10', compared with other LES result and the experimental data; CA - Catalano el al.13]; SZ - Szechenyi 120];

S11 -ShIli etal. [211; Z) - Zdravkovich 116])

E.

CD CL

-E.

St LES (present) 0.27 0.12 0.28 LES (CA) 0.32 0.32 0.35 Exp.(SZ) 0.25 Exp.(SH) 0.24 0.33 Exp. (ZD) 0.2/0.4 0.1/0.15 0.2/0.34 0.18/0.5 The global flow parameters predicted by the present LES are summarized in Table 3.

Our prediction of the mean drag coefficient (E.) closely matches the data of Szechenyi [20] and Shih et al. 121] interpolated at the present Reynolds number.

We consider their data more reliable than others in view of the sufficiently large spanwise lengths used (LID = 9.3 and 8.0, respectively), relatively low blockage ratio (7% and I 1%,

respectively),

surface smoothness, and relatively low turbulence intensity of the wind tunnels - 0.08 % for Shih et al.'s data [21].

The r.m.s. lift coefficient prediction by the present LES falls within the scatter of the experimental data compiled by Zdravkovich [16]. The base pressure is also seen to be in the range of the experimental data.

6 Copyright C 2005 by ASME

04 r

\\

.i.

A A I.

U A

-0.2-V.

-0.4 T

420 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 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.

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 present LES wvas found to be around 108 degrees, which is 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 early as 90 degrees, before it finally separates - "turbulently" -

at around 108 degrees.

l---Upper surfacel 5.0

[

ower surfaceI 0

30 60 90 120 150 180 (deg.)

Figure 11.

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

I I

-. I l

i I.

I I

I

-3 U

0 IM iVv Ho,1 I.'....

0a-'

2

'u 0 0.2 0.4 0.6 0.8 1

1.2 Strouhal Number 1.4 1.6 1.8 2

Figure 10. Power spectral density of the lift coefficient (Ct) for the critical Reynolds number (Re = 1.0 x 106)

For a smooth cylinder, the experimental data collected by Zdravkovich [16] and Achenbach's data [19] indicate that, above the Reynolds number of 1.5 x 106 referred to as "TrBIA" 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 roughness. At Re = 1.0 x 106 falling in the range referred to as "TrBL3", however, both sets of data seem to suggest that transition occurs near the separation.

0 (de".)

Figure 12. Resolved turbulent kinetic energy distribution at 0.OOOID above the cylinder surface, predicted the present LES at Re = 1.0 x 1O0 The skin-friction prediction, including the location of the peak and the magnitude, is also in good agreement with the experimental data of Achenbach [19] measured at Re = 8.5 x 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 distributions have the maximum values in the range of 3.5 - 4.0 occurring around 0= 60 - 65 degrees. As shown in Figure 11, our LES prediction closely matches Achenbach's data. This is in a sharp contrast with the LES prediction by Catalano et al.

13], which significantly overpredicted the overall level of the skin-friction on the front half of the cylinder. Contrary to what is suggested by the experimental evidence, their LES seem to yield too early a transition, giving a skin-friction level typical of a turbulent boundary layer from immediately downstream of the front stagnation point.

As they pointed out [3], one most likely cause for the discrepancy is the substantially coarser mesh (ye around 50) used in their LES predictions.

SUMMARY

AND CONCLUSIONS Turbulent flow around a smooth fixed circular cylinder was computed using large eddy simulation (LES) at two Reynolds numbers (Re = 1.4 x 105, 1.0 x 106).

The computations were carried out using a second-order accurate, parallelized finite-volume Navier-Stokes solver capable of handling arbitrary unstructured meshes. An implicit, non-iterative fractional-step method was employed in favor of its high efficiency in LES for incompressible flows.

The LES results for both Reynolds numbers predicted, with a good accuracy, the global flow parameters such as mean drag coefficient, r.m.s. lift coefficient, and the Strouhal number.

Furthermore, the salient features of the subcritical and the supercritical flows experimentally observed and measured on a smooth cylinder, such as the separation angle, mean velocity in the wake, length of the recirculation bubble, and transition location, are largely reproduced by the LES. Despite the use of a large number of computational elements (6.8 million cells),

the solution turnaround time was quite reasonable.

To summarize what has been achieved through this study:

The spatial discretization method (finite-volume on unstructured meshes) and the solution algorithm (implicit fractional-step method) employed in this study have been shown to be sufficiently robust and accurate to be used in LES for high Reynolds number flows.

The dynamic Smagorinsky model adapted to unstructured meshes using a novel test-filter [13] has been shown to work robustly and accurately for complex, high Reynolds number turbulent flows.

The present LES capability has been shown to predict the salient features of turbulent flow around a smooth circular cylinder at both subcritical and supercritical regimes with a good accuracy and reasonable computational cost.

It is planned to continue this work to include higher Reynolds numbers.

ACKNOWLEDGEMENT The authors acknowledge that FLUENT and GAMBIT software were used for this study.

The authors are also thankful to Dr. Mustafa Sarioglu at Karadeniz Technical University in Turkey for kindly allowing me to use his experimental data.

REFERENCE

[I] Breuer, M., 2000, "A Challenging Test Case for Large Eddy Simulation: High Reynolds Number Circular Cylinder Flow," Int'l J. Beat and Fluid Flow, 21, pp. 648 -

654.

[2] Cantwell, B. and Coles, D., 1983, "An Experimental Study of Entrainment and Transport in the Turbulent Near-Wake of a Circular Cylinder," J. Fluid Mech., 136, pp. 321 - 374.

[3] Catalano, P., Wang, M., Iccarino, G. and Moin, P., 2003, "Numerical Simulation of the Flow Around a Circular Cylinder at High Reynolds Numbers," Ini'l J. Hlear cnd Fluid Flow, 24, pp. 463 - 469.

[4] Travin, A., Shur, M., Strelets. M. and Spalart, P., 1999, "Detached-Eddy Simulation Past a Circular Cylinder,"

Flow, Turbulence and Combustion, 63, pp. 293 - 313.

[5] Vatsa, V. N. and Singer, B. A., 2003, "Evaluation of a Second-Order Accurate Navier-Stokes Code for Detached Eddy Simulation Past a Cricular Cylinder," AIAA Paper 2003-4085.

[6] Pandya, M. J., Frink, N. T., Abdol-Hamid, K. S., and Chung, J. J., 2004, "Recent Enhancements to USM3D Unstructured Flow Solver for Unsteady Flows," AIAA Paper 2004-5201.

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

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

[9] Park, J., Kwon, K., and Choi, H., 1998, "Numerical Solutions of Flow Past a Circular Cylinder at Reynolds Numbers up to 160," KSME lnt. J., 12, pp. 1200- 1205.

[I0]Norberg, C., 2003, "Fluctuating Lift on a Circular Cylinder: Review and New Measurements," J. Fluids and Structures, 17, pp. 57 - 96.

I I]Williamson, C.H.K., 1989, "Oblique and Parallel Modes of Vortex Shedding in the Wake of a Circular Cylinder at Low Reynolds Numbers," J. Fluid Mech., 206, pp. 579 -

627.

[12]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.

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

[14]West, G. S. and Apelt, C. J., 1993, "Measurements of fluctuating Pressures and Forces on a Circular Cylinder in 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 Cylinder, J. Fluid Mech., 234, pp. 191 - 217.

16]Zdravkovich, M. M.,

1997, "Flow Around Circular Cylinders, Fundamentals, Vol. 1, Oxford Univ. Press (Chapter 6).

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. Fluid Mech., 34, Part 4, 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:

2nd 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.,

Lebanon, NH 03766, U.S.A.

Xuelei Zhu, Fluent Inc., Lebanon, NH 03766, U.S.A.

Stefano Orsino, Fluent Inc.,

Lebanon, NH 03766, U.S.A.

ABSTRACT Large eddy simulations (LES) have been carried out for confined swirling coaxial jets discharged into a suddenly expanded pipe, which was studied experimentally by Roback and Johnson [1, 2].

The computations were made using a parallelized finite-volume-based Navier-Stokes solver that is second-order accurate in time and space, and permits use of unstructured meshes. The computational domain starts from an inlet placed upstream of the swirl generator, which makes the inlet boundary condition easy to specify.

Subgrid-scale turbulent stresses were modeled using a dynamic Smagorinsky model applicable to complex three-dimensional flows without any statistically homogeneous directions.

Subgrid-scale turbulent scalar flux is modeled using a constant Schmidt number in conjunction with the dynamically computed subgrid-scale turbulent viscosity. The LES predictions were found to closely reproduce the salient features of the flow and the species concentration downstream of the expansion. One of the conclusions was that a good resolution of the mean flow and turbulence in the upstream region is crucial in accurately predicting the mixing downstream of the expansion.

INTRODUCTION Confined swirling co-axial jets discharging into a suddenly expanded pipe [1, 2] have been studied numerically by several others using the Reynolds-Averaged Navier-Stokes (RANS) equations [3, 4] and large eddy simulation [6].

The flow in question consists of a non-swirling jet in the center, and an outer annual jet with a swirl imparted by a 8-vane, variable-angle swirl generator, with a swirl number of approximately 0.41.

Earlier numerical studies using RANS-based turbulence models [3, 4] have shown that the mean velocity field in the mixing region - the shear layer between the jets and the inner recirculation region - was predicted with a reasonably good accuracy.

However, the predictions of the mean species concentration were much less satisfactory.

The major discrepancy was found in the core region, where the mean species concentration was grossly underpredicted by the RANS computations.

Evidently, the RANS model employed gives rise to an excessive mixing. Brankovic et al. [3] investigated the sensitivity of their RANS predictions to such parameters as the turbulent Schmidt number and the inlet turbulence level.

However, the influence of these parameters was negligibly small.

The poor prediction of the species mixing has been attributed to the inability of RANS-based turbulence models to accurately represent the mixing by large-scale turbulent structure (or large eddies).

Playing a major role in mixing, large eddies are also responsible for undesirable yet somewhat subtler phenomena such as combustion instability and acoustic excitation inside combustors. To address these issues properly, one has to turn to high-level turbulent simulation like LES.

Akselvoll and Moin [5] used LES to compute the flow and the species concentration in non-swirling coaxialjets. Pierce and Moin [6]

attempted LES for the same confined coaxial swirling jet as is considered here. These studies demonstrated that LES can 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 r.m.s. species fluctuation which are important quantities in the context of turbulent combustion.

The present study concerns a LES computation for the swirling coaxial jets using a finite-volume solver that permits use of unstructured meshes.

Unlike the study of Pierce and Moin [6], however, our computational domain starts from an inlet boundary placed upstream of the swirl generator. Despite the larger solution domain and the implied increase in the computational resource, having the inlet boundary upstream of the swirl generator makes it straightforward to specify the inlet boundary conditions.

Meshing the computational domain including the swirl generator, which could become a difficult and time-consuming task with structured meshes, is made easier by the flexibility offered by unstructured mesh allowed by the present finite-volume solver.

This paper is, in many aspects, a progress report about an ongoing study whose ultimate goal is to find an optimal strategy based on the technique of LES for modeling turbulent flow and species mixing in coaxial combustors with or without swirl.

The paper is organized as follows. First, we briefly discuss the numerical method and the subgrid-scale turbulence models for stresses and scalar flux.

Special emphasis is laid on a transient algorithm whose efficiency significantly benefits the present study.

This is followed by the details of the LES computations regarding the choices of the solution domain, mesh (resolution), boundary conditions, time-step size, and the overall solution strategy. The results are then presented.

NUMERICAL METHOD The computation were carried out using the segregated solver in

FLUENT, a

general-purpose CFD software.

FLUENT employs a cell-centered finite-volume method based on a multi-dimensional linear reconstruction scheme that permits use of computational elements (cells) with arbitrary polyhedral cell topology including quadrilateral, hexahedral, 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 applying Green-Gauss theorem [9]. Diffusive and convective fluxes are discretized using central differencing [8].

An implicit fractional-step method (FSM)

[10]

in 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 equation applying an approximate factorization of the coupled Navier-Stokes equations. For incompressible flows, the FSM preserves the formal second-order temporal accuracy without having to perform, at each time-step, costly outer iterations to couple velocity-field and pressure. The FSM thus provides a highly efficient algorithm for CPU-intensive transient computations like LES.

Consider the semi-discrete form of the Navier-Stokes equations in "pressure-correction" ("

p"+I -p") form:

[D 0] 6p- )

0) where u"nt and r are the velocity vector and the momentum source vector, respectively, and the integer n is the time level.

A is the coefficient matrix defined in terms of the discrete convective and diffusive operators and the time-advancement scheme chosen, and G and D are the discrete gradient and divergence operators, respectively.

The coupled system of equations in Equation (1) is extremely difficult to solve as it is, since the matrix A has to be inverted for every iteration. In the fractional-step method, the original coupled equations in Equation (I) are approximated by Au=-Gp +r" AtDGtp"' = DO (3) u..I = u - ArG,7,+'

On a per-iteration basis, the series of operations in Equation (3) closely resembles SIMPLEC scheme.

The difference is that, in the iterative SIMPLEC scheme, the operations in Equation (3) are repeated in an outer loop until the all the solution variables converge, whereas the FSM needs only one sweep. To preserve second-order accuracy with the 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 individual equations and high-order source terms. Yet, non-iterative FSM is takes much less CPU time than iterative SIMPLEC scheme.

The system of discretized governing equations are solved using a point-implicit, Gauss-Seidel relaxation along with algebraic multi-grid (AMG) method to accelerate solution convergence. The N-S solver and the SGS turbulence model are fully parallelized.

SUBGRID-SCALE TURBULENCE MODELING For incompressible

flows, the filtered Navier-Stokes equations can be written as awi,, 2a = _I ap -1 IlY+

a Iv 2iL at arir waxy x

ax.4I ax, I (4) where r,, =uuj- -ii is the subgrid-scale turbulent stress. In this study, the subgrid-scale stress is modeled using isotropic eddy-viscosity as rig-TOr* -Cvl (5)

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

The filtered transport equation for a passive scalar is given by ap+ai7H = a ( aa qj at ax, ax, ax (6) where q1 is the subgrid-scale turbulent flux of the species concentration (V), and a the molecular diffusivity. The subgrid-scale flux is modeled using EDA 0 DG] [I A ](

u)

=[ (r+

[°(At )]

ID A:DGJ0 J

] 6p-+

)+L~

(2)

Sc, ax, (7)

Factorizing equation (2), we have a series of split operations as 2

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 likely to be caused by a large fluctuation of the model constant, we smoothed the model constant by applying the test-filter on it, and also clipped it so that the effective viscosity remains positive.

The additional details of the implementation of the DSM in the framework of an unstructured mesh based finite-volume solver can be found in the reference [8]. The DSM has been validated for a number of wall-bounded flows such as fully-developed channel flows and flows around bluff bodies such as a square cylinder and a sphere.

DETAILS OF COMPUTATIONS Solution domain, swirl vane geometry and meshes A partial view of the computational domain is depicted in Figure 1, along with the coordinate system. The inlet boundary is placed at L.OD upstream of the swirl generator, where D is the downstream pipe diameter (D = 122 nmn). The exit boundary is at 15D downstream of the pipe expansion. The computational domain thus has the swirl generator in it.

The entire inlet tube, 8 swirl vanes, and the downstream pipe were modeled without using any periodic boundaries.

pipe expansion Figure 2. A sectional view of the computational mesh (medium-sized mesh with 2.7 million elements) used in the computations.

The computational meshes were built using Gambit. We used three progressively refined, unstructured hexahedral meshes for the computations.

The three meshes have 1.1 million (coarse), 2.7 million (medium), and 4.8 million (fine) elements, respectively. A sectional view of the medium-size mesh is shown in Figure 2. The resolution of the medium mesh is such that, in the mixing region, it can resolve the integral length-scale estimated around 20 mm based on the experimental observation

[1, 2],

which is equivalent to D16, with approximately 20 elements.

It can resolve the smallest observable eddies of a size around 6 mmn, which were observed in the experiment [ I, 2], with 6 - 7 elements.

It is worth mentioning here that the mesh resolution in and around the inner and the annular tubes and the swirl vanes is too coarse to resolve the energy-containing eddies originating from and transported in that upstream region. We recognize that poor resolution of the turbulent structure in the upstream region will negatively impact the prediction in the downstream, 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.

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 of the expansion in the mixing occurring near the central core.

swuir vane

'7:,

I I I - - _.. t -,

_ - i -

-4:.

.:

inlet XLZ Figure 1. A partial view of the cor The swirl generator consists of on the hub with an equal spacing The blade-section has a camber an(

with the radius.

The blade geoi design data given in the Ref. [1]. i based on the digitized surface geome Boundary conditions, time-step size, and solution strategy 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 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],

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 8 identical vanes mounted not match the given volume flow-rate, and is inconceivably too in the azimuthal direction.

low considering that the measured axial velocity immediately d a vane-angle that change downstream of the expansion is around 0.8 ni/s.

netry was taken from the On wall boundaries, we employed a blended law-of-the-A NURB surface was built wall that invokes proper wall-laws depending on the local mesh 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

variables are extrapolated in a mass-conserving manner on this boundary.

The influence of the exit boundary condition was discussed at length by Pierce and Moin [6]. They found that the usual convective outflow boundary condition applied on an exit boundary with a cross-section of the downstream pipe yielded a central recirculation zone that is far smaller than what The LES was started using a steady RANS solution as the initial guess.

To accelerate the solution to a statistically stationary state, we superimposed a pseudo-random fluctuating velocity-field on the mean velocity field taken from the RANS solution.

Before the statistics are collected, the LES computation was run until the initial transients in the solution 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 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 the solution variables are extrapolated at "pressure-outlet" boundary that could be different from their treatment of the convective boundary.

The time-step size used in the present LES was determined based on the estimated characteristic time-scale of the smallest eddies to be resolved in the LES. We took v = 1.0 m/s and t

= 0.05D (D = 0.122 m) as the characteristic velocity-scale and length-scale of the smallest eddies. These estimates give an eddy-turnover time of r = /v = 0.006 second. It was finally decided to use a time-step size of At = 0.0002 second, which will resolves one eddy-turnover time (i) of the smallest resolved eddies roughly in 30 time steps.

Figure 5. Time-averaged velocity vectors in the r-z plane To obtain stable statistics of the solution, the transient computations were continued for a sufficiently long period of time, typically for more than 7 - 8 flow-through times, until the time-averaged velocity and species concentration fields largely recover an axisymmetry.

Regarding the mesh-dependency of the solutions, it was found that the LES results from the medium mesh (2.7 million cells) and the fine mesh (4.8 million cells) showed little difference, while the coarse mesh result deviates a little farther from other two.

Unless stated otherwise, the results presented in this paper are for the medium mesh.

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 shown to form in the shear layers between the inner jet and the outer annular jet, growing in size in the downstream direction.

At z = 25 nim, small eddies are confined near the shear layer between the inner and the outer annual jets. Yet, larger eddies are also seen to have formed in the annular recirculation region behind the step. One can visually tell from the figure that the 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 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 step-height as shown in the figure, also closely matches the experimentally observed value, 20 imm, quoted in ref. [1].

Figure 5 depicts the mean velocity vectors on a r-z plane.

Figure 8. Vortical flow structure in the mixing region visualized by the !so-surface of the second-invariant of the velocity deformation tensor, colored by velocity magnitude.

Figure 8 gives an overall impression of the turbulent vortical structure in the mixing region. This figure indicates that the flow in this coaxial jets are well mixed.

Velocity field Figure 9 shows the mean axial velocity profile along the centerline of the pipe. The LES prediction quite closely reproduces the general trend such as the rapid drop of the mean 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 species concentration on two crossflow planes and a r-z plane.

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],

Figure 9. Mean axial velocity along the centerline of the pipe.

z=5mm z=5 mm

.9U C

'U

>1 CMa U9 U4 i

Ca aCa 8

C

E z=25mm

>1 U

z=51 mm 2:

.C V

Ut

'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 azimuthal velocity components at three axial locations (z = 5 mm, 25 mm, 51 mm). The profiles shown in the figure were obtained by averaging the radial profiles taken at four azimuthal locations (0= 00, 90°, 180°, 270°). The predictions 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 locations of the peaks predicted by the LES are generally shifted toward the centerline.

This implies that the outer annular jet expands or opens up less than it does in reality.

This also means that the annular recirculation zone predicted by the present LES is longer (in the axial direction) than in the experiment.

The underprediction of the jet angle is also correlated with the downstream shift of the central recirculation zone discussed earlier.

z= 25 mm used in this study, including the fine mesh (4.8 million cells),

are fine enough to accurately resolve the energy-containing eddies generated in the upstream region.

The r.m.s. axial velocity fluctuation shown in Figure 11 seems to support our reasoning.

The r.m.s. axial velocity fluctuation at the centerline (r = 0) at z = 25 mm shown in the top-figure, and the one at z = 5 mm (not shown here) are severely underpredicted in the LES computation. As shown in the bottom-figure for z = 51 mm, further downstream of the jet exit, the r.m.s. axial velocity fluctuation catches up with the 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 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 eddies which would take part in "tearing" the inner jet.

i E

Z-

.9 a1 C1 rlr.

Figure 11. r.ms. axial velocity fluctuation nlow planes (z = 25 mm, z = 51 m) at two cross-We surmise that, among others things, the most likely culprit for this discrepancy is the lack of mesh resolution in the region upstream of the expansion, in and around the inner and the outer annular pipes, and the swirl vanes, As mentioned earlier, turbulent eddies coming from the upstream region feed 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 surrounding flow, which will lead to an increase in the jet angle. In our LES computations, these energy-feeding eddies are almost missing, since none of the computational meshes 50 100 150 Axial distance from the expansion (mm)

Figure 12. Mean species concentration along the centerline of the pipe The radial profiles of the mean species concentration at two axial locations are shown in Figure 13.

The predictions are in good agreement with the measurements at both locations. The mean species concentration profile at z = 25 mm predicted by the LES computation shows a sign of being "under-diffused", insofar 7

as it has a fuller profile than the measured one. The prediction at z = 51 mm is somewhat lower than the measurement near the centerline. However, the overall agreement at this location is remarkably good in view of the steep change of the mean species concentration near z = 51 mm, as can be seen in Figure 12. It should be noted, in passing, that the RANS predictions reported in the literature 13] severely under-predicted the mean species concentration at this axial location.

r.m.s. value at z = 25 nmn, located in the shear layer between the inner and the outer annular jets, is also considerably underpredicted.

Apparently, the LES underpredicts the entrainment of the ambient fluid occurring at this location. At z

= 51 mm, the prediction comes much closer to the measurement, as the r.m.s. velocity fluctuation did at the same 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 locations (z = 25 mm, 51 mm) at two axial Equally important - probably even more important than the mean species concentration in the context of modeling turbulent combustion, is the fluctuation of species concentration.

The usual RANS-based turbulence models 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 obvious benefit from LES is that one can directly compute it.

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 r.m.s. axial velocity fluctuation, the LES result grossly underpredicts the r.m.s. value in the core region.

The peak Figure 14. r.m.s. species concentration at two axial locations (z = 25mm, 51 mm)

The same remarks given regarding what might have caused the discrepancy in the r.m.s. velocity fluctuation (Figure 11) largely apply to the results for the r.m.s. fluctuating species concentration.

PRELMINARY RESULTS WITH LOCALLY REFINED MESHES The fact that the fine mesh (4.8 million cells) offers a meager improvement over the medium mesh (2.7 million cells) in terms of the mesh resolution (spacing) warrants additional 8

computations with substantially finer mesh, for instance, by halving the grid spacing in all three (axial, radial, and azimuthal) directions.

However, doubling the number of elements in the three directions results in an eight-fold increase in the total cell counts, which becomes unwieldy.

z=51 mm

'5..

2.0 1.51 1.0 0.5 0.0

-0.5

_ I Measured (Roback & Johnson [1])

LES (M.7M cell mesh)

-- LES (2.7M cell & locally refined)

' a 0~

X B1 I

I I

I oo CONCLUSION The flow and the species transport in confined swirling 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 conjunction with second-order accurate temporal and spatial discretization schemes.

A dynamic Smagorinsky model adapted to general three-dimensional flows was employed as the subgrid-scale turbulence model.

For LES, it is evidently a bold attempt and a costly proposition to include the upstream components of a combustor like the swirl generator and the upstream tubes in the computational domain.

Nevertheless, the present LES computations closely reproduce the salient features of the flow 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 the turbulence in and around the inner and the outer annular tubes, and the swirl vanes turned out to be more important than originally thought. The numerical evidence found from this study indicates that the three globally refined meshes used in the present study, despite the largest cell counts reaching up to 4.8 million cells, still cannot provide a sufficient resolution of 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 Inc. was used for the computations.

I 0.2 0.4 0.6 0.8 rTrN z=Sl mm r/rO Figure 15. Mean axial and azimuthal velocity profiles at two axial locations - predicted with a locally refined mesh To keep the mesh size under a tractable limit, the medium mesh was locally refined only in the domain upstream of the expansion.

Furthermore, only the cells within a specified distance from the wall, roughly 0.07D in this study, were refined to further save the cell counts. The local refinement resulted in a total of 6.5 million cells. The computation has been carried out on this mesh. The results are shown in Figure 15, being compared with the medium mesh results, for the radial profiles of the mean axial and azimuthal velocity components at z = 51 mm.

Clearly, the predictions with the adapted mesh are significantly improved over the predictions with the medium mesh. The overall shift of the velocity peaks is now far smaller than what we saw earlier, and the predictions capture the profiles of the mean velocity components much more closely.

REFERENCES

[I] Roback, R. and Johnson, B. V., "Mass and Momentum Turbulent Transport Experiments with Confined Swirling Coaxial Jets," NASA CR-168252, 1983.

12] Johnson, B. V. and Roback, R., "Mass and Momentum Turbulent Transport Experiments with Confined Swirling Coaxial Jets -

Part 1," AIAA-84-1380, Presented at AIAA/SAE/ASME 2 0 th Joint Ppropulsion Conference, Cincinnati, OH, June 11 - 13, 1984.

[3] Brankovic, A., Ryder, Jr. R. C., and Syed, S. A., 1998, "Mixing and Combustion Modeling for Gas Turbine Combustors Using Unstructured CFD Technology," AIAA Paper 98-3854, Presented at 341h AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, July 13 - 15, 1998, Cleveland, OH.

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.

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 W

JE encounter many industrial applications of computational fluid dynamics (CFD) where the flows are dominated v

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.

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 Lilly7 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) = V l y

e E

(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

_i+ _

I o atij

((a a5)

a.

axj P axi axj axj HODj ax; = °(6) where ¶ij is the subgrid-scale stress defined by;

,ij _7-j -ij (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 v,

= c.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 p a.

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)

Tij-TTkk = -2ca 2 Isis (14) where a is the test-filter width. Note that the same model coefficient, Cv, is used in the expressions for both -Tj and 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 (

S2 sl 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 = Cc kEST (29)

The dissipation-rate of kel can also be computed from e = (v +v,)

a i -,

(30) axJ ax}

axi ax From Equation (29) and Equation (30), CE is given in a closed form kest axi axi axi axi C

k(t/,)

(.Axj at,-.;;at 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 solved ie-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 The flow past a square cylinder measured by Lyn et aL24 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 turbulent flows around bluff bodies with sharp edges. Some others have also tackled this flow using LES.13. 25 The domain size and the mesh resolution were chosen in reference to the earlier studies by others. 3 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 time unit (H/Uo) which is comparable to that used by others.13 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 shown in the figure, this value is considerably larger than what the measurement indicates (0.62Uo).24 Interestingly, others who computed the same flowt3.25 also have grossly overpredicted the recovery of the axial mean velocity. Our LES predictions of the mean axial velocity in the far-wake with both dynamic models were found to be largely com-parable to the prediction by Sohankar et aL.13 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

1111 E11111!3!

U1111 Tl ll 1

I I

I

(.

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

LES with DSMI (prexent) l

--- LES with LDICEM (prcscnt)_

0.1 /

a Lyn et al. (1995) 0.2 -

O.I I

1.B

°.0 5.0 10.0

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

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 13 27

.lO.65m Measured2 9 0.15 -0.16 0.629 The results are summarized in Table 2 along with others' predictions. 262 8 Tomboulides et aL27 and Johnson and Patel2 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 Measured31.3 2 0.40 0.185 - 0.19 Correlation33 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's3 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 FL 0.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 data32 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 l lR0 ~

_ _ _ _1-0.5 I

2Experimental 0.50 10 LU m

LU.

, LoESwith DSM

_X 0.5 jJami10 a c j,.& 1 l]l. 1. fit I

J 045 -

S A 0 40 -

0.3~ J, I

0to 200 300 400 500 600 t Uu/D

,.l!

e10 II l I i

@ fl I IM II 11 1 I 1et 10 10"_

as P _.=

. I l

= ::

lii-

- miiL

= :i 0IM 21 A

.2

]lll 1.ii I EaS I

lt 102 io leo ReD (a) Time histories of CD for Re0 = 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.

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

3Urbin, 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.

4K. 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.

6M. 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.

7Lilly. 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.

9Kim, 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.

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

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

13A. 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.

t4Krajnovic. 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.

17Kader, 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." NumericallHeat Transfer. Vol.31.1997. pp. 195-215.

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

'9S.-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.

24D. 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.

26johnson. 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.

3tAchenbach. 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.

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

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

35Kim, 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

ML052220086

Contents

Text

Subject:

Subject:

Subject:

SUMMARY

Navigation menu

ML052220086

Text

Subject:

Subject:

Subject:

SUMMARY

Navigation menu

Search