Jingxin Zhang, Dongfang Liang, Hua Liu

(1. MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China 2. Key Laboratory of Estuarine & Coastal Engineering, Ministry of Transport, Shanghai 201201, China 3. School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

?

A Detached Eddy Simulation Model forFree Surface Flows

Jingxin Zhang1,2,3,*, Dongfang Liang1,3, Hua Liu1,3

(1. MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China 2. Key Laboratory of Estuarine & Coastal Engineering, Ministry of Transport, Shanghai 201201, China 3. School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

Geophysical flows in oceanic shelves, estuaries, and rivers are often studied with the shallow water equations under either hydrostatic or hydrodynamic pressure assumptions. The Reynolds-Averaged Navier Stokes (RANS) equations are widely used for simulating the complex turbulent flows because of their high efficiency. The Directed Numerical Simulation (DNS) and Large Eddy Simulation (LES) are commomly unaffordable under the constraint of the present computer power. In practical applications, one kind of hybrid models, i.e. Detatched-Eddy Simulation (DES), is a good compromise because of the lower computational cost for simulating near-wall flows with RANS. In the paper, a DES model is proposed based on a fully hydrodynamic pressure model instead of hydrostatic model. The numerical scheme is based on the Finite Volume Method (FVM) on unstructured grids in the horizontal plane, and the σ coordinate in the vertical direction to fix the free surface and the uneven bottom. The in-house codes are paralleled using OpenMP. The proposed model is shown to be particularly effective in the prediction of small-scale vortical structures.

Detached-eddy simulation, Coherent structure, Free surface flow

0 Introduction

Prediction of shallow flows with free surface is important for practical applications in coastal and civil engineering for the design of nautical constructions, management of waterways and prevention of natural disasters such as floods and accidental spills. Numerical modeling is one of the primary tools used in engineering practice for the prediction of shallow flows. The Reynolds-Averaged Navier Stokes (RANS)

equation model is widely adopted for turbulent flows because of its affordable computational cost, especially for flows with complex geometry with high Reynolds number. It is efficient for RANS to predict the time-averaged flow properties, however, some turbulent characteristics, i.e. coherent structures or vortical structures can’t be exactly predicted. The small vortical structures are not resolved in these averaged models. Compared to RANS, the Large-Eddy Simulation (LES) is more advanced for the investigation of small-scale vortical structures, however, its computational cost increases sharply with high Reynolds number and geometric complexity. The coupling of the Reynolds-Averaged Navier Stokes (RANS) equation model and the Large Eddy Simulation (LES) is a natural strategy in a wide range of complex flow simulations with higher Reynolds number. As one of the hybrid RANS and LES models, Spalart et al[1] proposed the Detached-Eddy Simulation (DES) by means of modifying the turbulence length in the primary transport equation of the eddy viscosity devised by Spalart et al[2]. By now, there are altogether three progressive models in the DES family, i.e. DES97, DDES and IDDES [3-4]. Spalart [5] made a comprehensive review about the developments, applications and problems of DES. The DES, DDES or IDDES, has been successfully applied in aerodynamics, but there is few applications in hydrodynamics. For environmental hydrodynamics, the RANS model is still the preferential numerical tool for turbulence prediction.

For large-scale shallow water flows, the hydrostatic model is widely used for its efficiency and accuracy, however, the non-hydrostatic pressure is vital in some cases, where the vertical acceleration cannot be neglected compared to the gravity. To improve the predictions of shallow flows, the model needs to account for a non-hydrostatic pressure. Casulli [6-7] reported on a model with fully hydrodynamic pressure terms in which the pressure is represented as the sum of the hydrostatic and non-hydrostatic constituents. Jankowski [8] presented the details of a predictor-corrector method for calculating the pressure field. The predictor step is used to calculate the hydrostatic pressure, while the non-hydrostatic pressure is calculated in the following corrector step. Li [9] simulated water waves by directly solving Navier-Stokes equations using a similar decomposition method and adopting a fractional method to implement the numerical model. Kocyigit et al[10] and Chen et al[11] solved the three-dimensional non-hydrostatic equations in the Cartesian coordinate system. Fringer et al[12] presented a non-hydrostatic model for ocean flows based on unstructured grids with the finite volume method (FVM), which was implemented in a highly-efficient parallel computational code. Because of these improvements, non-hydrostatic modeling becomes increasingly feasible for flows in the natural environments. For large-scale shallow water flow simulation, the grid size is usually too coarse to depict local geometry, such as small dunes. In other words, the uneven bottom is actually smoothed out in the numerical simulation. However, the grid size must be fine enough in DES simulations, and some small geometry should be accurately depicted. As a result, the local hydrostatic assumption is no longer valid, and it is important to extend the hydrostatic model to the fully hydrodynamic model.

With apurpose to investigate small turbulent flow structures in geophysical flows, an in-house code was first established based on the hydrostatic pressure, which had been widely used in hydraulic engineering, and then the non-hydrostatic pressure was introduced to update the model for flows with complex boundaries. The DES strategy has been implemented in the current version, and the model can be easily switched from RANS to local LES. So far, it is not practical to carry out the DES simulations of free-surface flows in a natural river with decades of meters width and kilometers length, but it is feasible for the simulations in laboratory scales to investigate turbulent flow mechanisms.

1 Numerical Model

1.1 Governing Equations

Following the work of Casulli [6], the total pressure is split into the hydrostatic componentph=ρg(ζ-z)andnonhydrostaticcomponentpn. The three dimensional governing equations in a rotating frame can be written in a conversation formwheregisgravitationalacceleration,andf=2ωsinφistheCoriolisforceterm,whereφisthelatitude,andωistheEarth’sangularvelocity.ζis the free surface elevation. When the nonhydrostatic componentpnisignored,theverticalmomentumequationisalsoneglected,andthesystemdegeneratestothecommonshallowwaterequationswithonlytwohorizontalmomentumequationsandonecontinuityequation.

(1)

(2)

(3)

(4)

Inordertofitthefreesurfaceandtheunevenbottomboundary,theverticalcoordinatezistransformedtotheσcoordinate [13], and the transformed equations are written as:

(5)

(6)

(7)

(8)

(9)

1.2 Turbulence Model

(10)

Theeddyviscosityυtisgivenby

(11)

where

υis the molecular viscosity.

The functionfwiscalculatedas

(12)

[21]simulatedtheopenchannelflowwithfreesurfacebyalevel-setmethod,inwhichtheairmotionwasconsidered.Theresultsrevealthattheeddyviscositynearthefreesurfaceisaboutseveraltimesthemolecularviscosity.

1.3 Overwhelm the Model-Stress Depletion (MSD) by Zonal-Detached-Eddy Simulation (ZDES)

The MSD is one severe problem with the RANS/LES hybrid model, which is pronounced for the DES when the switch from RANS to LES takes place in the boundary layer, so Spalart [3] proposed a new version of DES, i.e. DDES, to delay the switch from RANS to LES. Although the MSD issue is partially overcome, the “LES-content” is not sufficient in the boundary layer with the model running in RANS mode. Shur et al.

[4] developed a new model, named IDDES, for the simulation with an ambiguous grid scale, which can solve the MSD problem by modifying the turbulent length scale, initial condition and inflow condition. In fact, the MSD issue comes from the insufficient velocity fluctuation when the simulation is switched from RANS to LES. Generating fluctuations by numerical method is one feasible way to introduce much more “LES-content” at the interface between RANS and LES. Keating et al.

[22] presented a dynamic stochastic forcing method, which significantly speeds up the transition resulting in more accurate predictions of the velocity fluctuation. Compared to numerical fluctuation generator, the Zonal-DES approach is well adopted to handle separated flows in which strong instabilities develop rapidly, thus overwhelming the turbulence inherited from upstream boundary layers. Breuer et al.

fv1=1,fv2=0,fw=1

(13)

Consequently,inthefirststepofmodelrunningthefunctionsfv1,fv2andfwweresetinaccordancetoEq.(13)intheLESzone.

1.4 Numerical Scheme

1.4.1 Coordinate Transformation

J=xξyη-xηyξ

(14)

whereξdirectsfromthecontrolcellcentertotheneighboringcellcenteracrossthecommonface,andηisfromonevertextotheotheronealongthefacepointinginananticlockwisedirectionaroundthecontrolcellcenter(Fig.1).

Fig.1 Local coordinate on the control cell face.

1.4.2 Predictor Step for Flow Under Hydrostatic Pressure

The semi-implicit method is introduced by a parameterθ[27].Theimportanceofthevalueofθhasbeenshowninpreviousresearches[10, 11, 27].Thevalueofθistakentobe0.5inthefollowingcomputations.

Firstly,thecontinuityequationisdiscretizedinto:

(15)

andthediscretizedmomentumequationsare:

(16)

(17)

(18)

whereFisanoperatorthatincludestheexplicitdiscretizationoftheconvectiveterms,thehorizontalviscoustermsandtheCoriolisforceterm.

(19)

andthediscretizedmomentumequationsaregivenby

(20)

(21)

(22)

Further,discretizingtheaboveequationsgives:

(23)

(24)

(25)

Hence,thediscretizedequationscanbewritteninacompactMatrixnotation:

(26)

(27)

(28)

(29)

Thedetailsofthematricesareasfollows:

Aiis a series of tri-diagonal matrices.

(30)

(31)

(32)

(33)

RewritingthegoverningEq. (26)into:

(34)

(35)

IntegratingEq.(35)overthehorizontalcontrolvolumebasedonGauss’stheorem:

(36)

wherethesymbolfindicates the cell face, ΔSiis the horizontal area of the celli, Δlisis the length of theSthface,NSis the total number of faces, and (cosαis, sinαis) is thex- andy- projection of the face’s normal unit vector. The derivative with respect to the Cartesian coordinate (x,y) is then transformed to that with respect to (ξ,η), and then the Eq. (36) is deduced as

(37)

Eq. (37)iswritteninacompactform:

(38)

where

1.4.3 Corrector Step for Flow Under Nonhydrostatic Pressure

The fully discretized equations at time stepn+1, following the resolved hydrostatic fluid flow, can be re-written as:

(39)

(40)

(41)

SubtractingEqs. (20), (21)and(22)from(39), (40)and(41),onecandeducethefollowingdiscretizedequationsinthecorrectorstepforflowsinducedbythenonhydrostaticpressurepn.

(42)

(43)

(44)

Although the predicted flow field satisfies the depth-averaged continuity, it does not meet the local continuity condition. Therefore, a nonhydrostatic pressure must enforce the corrected flow to satisfy the continuity condition. In the transformedσframe,thecontinuityEq. (5)doesnotexplicitlyincludetheverticalvelocityqz,andanotherformwhichcontainingqzcanbededuced:

(45)

Someresearchers[10]adoptedacontinuityequationintheCartesiancoordinate,whichwasthesameastheaboveequationbutonlypreservedthefirstthreeterms,andtheσcoordinatewasreplacedbyzcoordinate, i.e. ?z=D?σ.ItwasnotveryclearinKocyigit’spaperhowtodealwiththeresidualtermsintroducedbythecoordinatetransformationbetweentheCartesianandσcoordinates.IntheCartesiancoordinate,thevariablesatcontrolpointsmustbeinterpolatedateverytimestepbecauseofthetransformationbetweenthedifferentcoordinates,whichisasimilartechniquepresentedbyStansby[28]orHuangetal.

[29].Ifitonlyreplaces?z=D?σandneglectsthelastterminEq.(45),somepracticaltestsverifiedthattheaccuracyisnotgood,especiallyforcaseswithsharpchangeofbottomelevations.

(46)

(47)

where

Atthesolidwallandbottom,azeronormalgradientisimposed,andatthefreesurface,azerononhydrostaticpressureisspecified.ThematrixsystemiscalculatedbythepreconditioningBi-CGSTABmethod.

Atlast,thenewwatersurfacecanbeupdatedafterthecorrectorstep.Casulli[7]presentedonemethodforupdatingthewatersurfaceafterthecorrectorprocedure,andZhangetal.

[30]implementedanotherthirdstepsimilartothefirststeptoupdatethenewsurfaceelevationsandvelocities.Inthefollowingcasestudies,itdoesnotrevealtheobviousimprovementwithathirdupdating,andthusithasbeenturnedoffinthecodes.

1.4.4TVDschemefortheadvectionterms

Fortheadvectionterms,theTVDschemeisusedtocalculatethecellfaceflux.Anarbitraryvariable〈φ〉fatthecellfaceiscalculatedbymeansofasecondorderTVDscheme:

(48)

wheresymbolfindicatesthecellface,φDandφCarerespectivelythecell-centeredvaluesinthedownwindandupwindnodesofthefacef.Usingafluxlimiterfunction,ψ(rf),whichissimplyalinearfunctionofrf,onecanconstructhighorderTVDschemes.TheratiorfiscalculatedbymethodextendedbyDarwishetal.

[31].ThefollowingsecondorderTVDschemeshavebeenimplementedinthepresentmodel

SUPERBEElimiter

MINMODlimiter

OSHERlimiter

MUSCLlimiter

VanLEERlimiter

SWEBYlimiter

QUICKlimiter

MC limiter

For the CC scheme, the problem induced by the pressure-velocity coupling on the collocated grid is avoided by means of an interpolation scheme suggested by Rhie and Chow [32].

2 Applications

Two test cases were explored in this computational study to validate the numerical model and to investigate the coherent vortical structures of turbulent water flows with free surface. The first test case is about flows over a series of dunes, and the second one is focused on turbulent flows at the channel confluences with a 90° junction.

2.1 Open Channel Flow Over Series of Dunes

The single dune geometry is the same as that in the experiments of Balachandar et al.

[33], in which the measurements were made at the 17th dune of a train of 22 dunes mounted in a hydraulic flume. Fig.2 shows the single dune geometry on the central vertical plane. A series of five dunes are designed in the study case, and the simulation domain is then extended upstream and downstream in order to yield a fully developed turbulent flow. The horizontal grids with a finest 5mm×5mm resolution cover the dune zone, and then are coarsened to 20cm×20cm when reaching the upstream and downstream open boundaries. The design of the computational grids ensures the LES model for flow over the dunes and the RANS model beyond the dunes (see Fig.3). One idea is to investigate the influence of upstream velocity fluctuation on thevortical structure simulated in the following dunes. The flow passing one dune contains fluctuations, which can be considered as a physical method to generate velocity fluctuations. For natural river flows, the curvilinear riverside and uneven bottom provides the driving force to generate the velocity fluctuation. Based on this idea, the fluctuation introduced at the inlet is not vital only if there is enough space before the focusing zone for turbulent flows to develop.

Fig.2 Dune geometry and measurement locations(mm).

Fig.3 Schematic of computation domain division.

The Reynolds number, based on the water depthLand the free-surface velocityU0attheinlet,isabout5.7×104,andtheFroudenumber,whichisbasedonthesamevariables,isabout0.44.Inthevertical,thefirstgridpointtothebottomismaintainedatabout1.0z+, and a stretching ratio between grid cells in the log layer is about 1.15. The 5mm grid size is selected by means of a few tests to ensure that the mode is switched from RANS to LES at tens ofz+, which is important for DES application [22]. For this open channel flow, the horizontal grid size in the LES zone is about 1/24L, which has been calibrated to be fine enough for LES [16]. In order to improve the simulation efficiency, the simulation was firstly run with the RANS mode with hydrostatic pressure assumption until computation converged, and then the non-hydrostatic mode was turned on. Later on, the mode was switched from RANS to DES, and the computation lasted about 15 large-eddy turnover time (L/uτ), in whichuτis the friction velocity.

2.2 Time-Averaged Mean Flow

The results simulated by RANS are ultimately steady in the present case. However, the flow simulated by DES is unsteady because of the velocity fluctuation successfully simulated, and the flow signals during 15 large-eddy turnover time were time-averaged to obtain the mean flow. In experiments, altogether six vertical lines were laid to measure the velocity.

Three groups of streamwise velocities at six locations were extracted from the first, third and fifth dune zone, respectively, with the relative distances to the crest of the dunebeing the same as those in experiments, i.e.x/h=2, 4, 5, 6, 12, and 18. Fig.4 shows the comparison of streamwise velocity simulated by DES with the experimental data, and in each plot, the three lines are respectively for the data obtained at the first, third and fifth dunes. At each location, the calculated velocity at the fifth dune matches the measurements with the highest accuracy. For the grid design in the study case, the inlet was covered by RANS, and no velocity fluctuation was input. As a result, insufficient “LES-content” was predicted over the first dune, but the velocity fluctuation inspired locally can be transported into the following dunes, which acts like a velocity fluctuation generator.

The mean Reynolds shear stress, 〈-u′w′〉iscalculatedbytime-averagingthesignalofinstan-taneousReynoldsshearstress-u′w′calculatedfromtheresolvedturbulentscale.Fig.5showstheprofilesof〈-u′w′〉alongdifferentverticallinesatthefirst,thirdandfifthdune,togetherwiththemeasureddata.Thepeakoccurredataboutz/L=0 is clearly predicted numerically. The location of the predicted peak is well in agreement with the experimental data, although the magnitude is considerably smaller especially over the first dune. The improvement of the agreement between measurements and calculations from the first to fifth dune verifies the velocity fluctuation from upstream evidently affects the turbulence simulation over the downstream dunes. Only the results by ZDES are presented here because of the superiority of ZDES to DES. The improvement of ZDES lies in the much more “LES-content” inspired at the interface between RANS and LES, indicating that much more turbulent fluctuation can be directly simulated. Compared to the experimental data, the calculated

Fig.4 Comparisons of stream velocity at different sites.

Fig.5 Comparisons of predicted mean Reynolds stress with experimental data.

Reynolds stress over the third dune is the most accurate among the three selected dunes, and there are a pronounced overshoot over the fifth dune, especially at the locationx/h=2. The mean flow velocity results suggest that the accuracy will be gradually improved from the first dune to the fifth dune because the turbulent flow gradually develops downstream, which, however, is not supported by the predicted mean Reynolds stress. One fact should be highlighted that the unsteady turbulent flow has not converged to a steady statistic status until the fifth dune in the present study, so it cannot be concluded that the accuracy will decrease with extending dunes. In fact, the experimental data were collected at the 17th dune to ensure the turbulent flow fully developing.

2.3 Instantaneous Flow

It is expected that small-scale vortical structures are predicted by the DES model. A positive isovalue of criterionQis used to show the turbulent structures, which defines as the vortex tubes in the

regions where the second invariant of the velocity gradient tensor is positive, i.e.Q=(ΩijΩij-SijSij)/2>0, in whichSijandΩijrepresentthesymmetricandasymmetriccomponentsofV, respectively. Fig.6 shows a snapshot of the instantaneous isosurface ofQ=5whichclearlyshowstheevolutionofvorticalstructuresfromthefirsttothefifthdune.The“LES-content”isgraduallyinspireddownstream.Commonly,turbulencefluctuationmustbeinputattheupstreamforLEStoinvestigateflowoverdunes.Althoughonlythemeanflowboundaryconditionwasspecifiedattheinlet,thelocaldunecaninspireturbulence,andcanbetakenasakindofupstreamboundaryconditionforthenextdunes,whichcanbeconsideredasphysicalmethodtogenerateturbulenceotherthannumericalmethod.TheauthorsdidthesimilarsimulationonlyturningontheRANSmodel,andtherewerelittlediscrepancyoftheflowstructuresamongthefivedunes.However,boththeinstantaneousvorticalstructureandthemeanvelocitydistribution(Figs.4and5)revealexplicitdifferenceamongthefivedunes.

Fig.6 Snapshot of the criterion Q isosurfaces contoured by vorticity.

(a) 1st dune

(b) 3rd dune

(c) 5th duneFig.7 Instantaneous flows at different cross-sections (x/λ=0.2, with λ the dune length).

(a) 1st dune

(b) 3rd dune

(c) 5th duneFig.8 Instantaneous velocity vectors (u′,w′) in the central x-z plane.

ThedifferencebetweenRANSandDESagainillustratestheimportanceoftheinjectingsufficient“LES-content”attheupstreamboundary.

Selectingthecross-sectionsatarelativelocation0.2λin each dune, and the instantaneous velocity and isolines of vorticity simulated are shown in Fig.7. These plots again emphasize that the upstream turbulence fluctuation affects the “LES-content” simulation.

Subtracting the mean velocity from the instantaneous flow, Fig.8 reveals the instantaneous fluctuation velocity (u′,w′) in the vertical central plane (y=0). Abundant turbulent coherent structures are predicted in the fifth dune rather than that in the first one, which further contributes to the illustration of the above view points. Drawing an analogy between the velocity fluctuation generated by numerical method and by the frontal four dunes, the fluctuation injected into the fifth dune zone can be considered as a result of a physical method instead of numerical method. Fortunately, in natural rivers, there are sufficient disturbing forces induced by uneven bottom and curvilinear river bank, which can generate velocity fluctuation for the flow simulation in the studied area by LES, provided that there is enough upstream length for turbulence to fully develop.

2.4 Flow at a 90° Open- Channel Junction

The junction of two open channels is a common occurrence in many natural and hydraulic structures, especially in waste water treatment facilities, and fish passage conveyance structures. There are complex turbulent flow structures in the junction flow (Fig.9). Several pieces of work on CFD of the flows at open channel junctions were reported, and most of them were based on RANS. LES is generally superior to RANS in the prediction of complex turbulent flows, but the cost of a LES simulation is much greater than that of a RANS. The DES simulation was carried out to satisfy a balance between the computational effort and the prediction accuracy, especially, in the flow separation zone.

Fig.9 Schematic flow structure (after Weber et al., 2001).

Fig.10 View of the horizontal local meshes.

The RANS simulations resolve only the mean motion, and a turbulence model is used to account for the effects of the turbulence on the mean flows, so certain vortical structures cannot be obtained. The superiority of DES over RANS is a more complete description of the turbulence characteristics, especially, the small vortical structures. For the open-channel junction flow, two streams combine at the junction, and then form a mixing layer downstream. The mixing process depends on the ratio of the discharge and momentum between the main and side channels. Fig.11 shows the vortical structures by means of a criterionQ, which is contoured by vorticity magnitude. The left panel reveals abundant “LES-content”. There are two separation points at the confluence for the side flow, one occurs at the left corner of the side channel, and another one occurs at the right corner. The separation flow originated from the left corner forms a circulation in the main channel with a reattached point downstream the channel. The vortical structures have a complex spatial and temporal evolution. The separation flow from the right corner mixes with the main flow and forms stream wise oriented vortical cells. In RANS simulation, the most of small turbulent structures are wiped off as a result of the time-averaging, which is evidently shown in the right panel as a comparison with DES simulation. The free surface positions shown in Fig.12 also present fluctuations, which can be simulated by DES rather than by RANS. The instantaneous upwelling and downdraft are predicted by DES simulation, however, the RANS simulation shows a much smoother free surface in a quasi-steady status.

Fig.11 Snapshot of the criterion Qisosurface contoured by vorticity.

Fig.12 View of free surface.

The eddies are the basic dynamic cells to advect and diffuse mass and momentum. The eddy intensity is quantified by vorticity. Fig. 13 reveals the vorticity distribution in two different horizontal planes, one near the surface and the other near the bottom. In the near surface plane, the shear line, defined as the centre line of a mixing layer formed by the main and side channel flow, outlines a zone with circulation flow justly downstream of the confluence. The contracted flow is firstly formed and then develops in a distance downstream. When passing the reattachment point, the two streams gradually mix well until an even distribution of the vorticity across whole channel width. The vorticity distribution near the bottom also reveals the similar evolution of the eddies (right panel in Fig.13). The evolution of junction flow is further illustrated in Fig.14, in which the velocity fields in different cross-sections along the main stream wise are plotted with the simultaneous vorticity contour. The velocity field clearly depicts the evolution of eddies downstream the channel.

Fig.13 Vorticity in different horizontal planes.

The large eddies induced by junction flows control the mixing layer evolution formed by two streams, which is important to investigate waste water draining, local bed erosion and so on. The RANS simulation is based on the time-averaging technique, and the small vortical structures are commonly under-predicted. The DES simulation is undoubtfully an advanced method for investigating the turbulence mechanisms with a superiority in prediction of turbulent coherent structures.

Fig.14 Instantaneous velocity and vorticity at different sections.

3 Conclusions and Future Work

A hybrid RANS/LES model, DES, was developed to simulate the open channel hydrodynamics. The numerical model was extended from the hydrostatic model to the fully hydrodynamic model. The numerical strategy is based on the FVM with unstructured grids to fit irregular boundaries, which is common in geophysical free surface flows. The current version of the code is implemented using OpenMP, and it is expected to be further extended by combining MPI and OpenMP techniques to take advantage of the more and more advanced hardware.

In the framework of DES, the velocity fluctuation, i.e. the “LES-content” introduced upstream, affects the prediction of small-scale vortical structures, which is verified by the test cases. However, the inflow condition can be compensated by physical disturbing in the upstream part of the domain. The flow passing a dune generates velocity fluctuation, which can be considered as a physical fluctuation generator instead of a numerical generator. As a velocity fluctuation generator at the interface, ZDES speeds up the transition from RANS to LES, and more “LES-content” can be produced. ZDES improves the accuracy of turbulent flow simulation, in terms of not only the mean flow, but also the instantaneous flow, especially, the turbulent vortical structures.

The natural geophysical flow with free surface contains a large range of temporal and spatial length scales, and it is unaffordable to carry out LES, especially for cases with complex boundaries and high Reynolds number. It is a feasible technique to do LES only in the zone of interest, and RANS simulation elsewhere, which is the idea of DES. As mentioned in the study case, the insufficient fluctuation from the upstream is one severe problem, but the fluctuation can be gradually generated by extending the upstream computational domain. Fortunately, the irregular bank and uneven bottom act as physical generator of turbulent fluctuations in natural geophysical flows.

Commonly, the advantage of DES over RANS is to predict the small-scale vortical structures rather than the mean flow properties. The simulation of junction flow, which is always used as hydraulic engineering, reveals that abundant “LES-content” can be well predicted. It is expected to take advantage of the DES to investigate much more turbulence characteristics for natural flows with free surface, and to contribute to research of some environmental fluid mechanics issues.

[1]Spalart P R, Jou W H, Strelets M, et al. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[M]//Liu C, Liu Z, Advances in DNS/LES. Greyden Press, 1997.

[2]Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows[J]. La Recherche Aérospatiale, 1994, 1: 5-21.

[3]Spalart P R, Deck S, Shur M L, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J]. Theor. Comput. Fluid Dyn., 2006, 20: 181-195.

[4]Shur K L, Spalart P R, Strelets M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capa-bilities[J]. International Journal of Heat and Fluid Flow, 2008, 29: 1638-1649.

[5]Spalart P R. Detached-Eddy simulation[J]. Annu. Rev. Fluid Mech., 2009, 41: 181-202.

[6]Casulli V. A semi-implicit finite difference method for non-hydrostatic, free-surface flows[J]. International Journal for Numerical Methods in Fluids 1999, 30:425-440.

[7]Casulli V, Zanolli P. Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems[J]. Mathematical and computer modeling, 2002; 36: 1131-1149.

[8]Jankowski J A. A non-hydrostatic model for free surface flows[D]. Germany, University of Hannover, Ph.D. Dissertation, 1999.

[9]Li B, Fleming C A. Three-dimensional model of Navier-Stokes equations for water waves[J]. Journal of Waterway, Port, Coastal and ocean Engineering, 2001, 127(1): 16-25.

[10]Kocyigit M B, Falconer R A, Lin B. Three-dimensional numerical modeling of free surface flows with non-hydrostatic pressure[J]. International Journal for Numerical Methods in Fluids, 2002; 40: 1145-1162.

[11]Chen X J. A fully hydrodynamic model for three-dimensional, free-surface flows[J]. International Journal for Numerical Methods in Fluids, 2003, 42: 929-952.

[12]Fringer O B, Gerritsen M, Street R L. An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator[J]. Ocean Modeling, 2006, 14: 139-173.

[13]Phillips N A. A coordinate system having some special advantages for numerical forecasting[J]. Journal of Meteorology, 1957, 14: 184-185.

[14]Spalart P R. Strategies for turbulence modelling and simulations[J]. International Journal of Heat Fluid Flow, 2000, 21: 252-263.

[15]Shur M, Spalart P R, Strelets M, et al.Detached-eddy simulation of an airfoil at high angle ofattack[C]//Rodi W, Laurence D. Fourth International Symposium on Engineering Turbulence Modelling and Experiments.Corsica, 1999 (Elsevier: New York, 1999).

[16]Hinterberger C, Frohlich J, Rodi W. Three-dimensional and depth-averaged Large-Eddy Simulations of some shallow water flows[J]. Journal of Hydraulic engineering, 2007, 133, 8: 857-872.

[17]Pope S B. Turbulent flows[M]. Cambridge University Press, 2000.

[18]Spalart P R, Rumsey C L. Effective inflow conditions for turbulence models in aerodynamic calculations[J]. AIAA Journal, 2007, 45(10): 2544-2553.

[19]Aupoix B, Spalart P R. Extensions of the Spalart-Allmaras turbulence model to account for wall roughness[J]. International Journal of Heat and Fluid Flow, 2003, 24: 454-462.

[20]Eca L, Hoekstra M, Hay A, et al. A manufactured solution for a two-dimensional steady wall-bounded incompressible turbulent flow[J]. International Journal of Computational Fluid Dynamics, 2007, 21(3-4): 175-188.

[21]Yue W, Lin C L, Patel V C. Large eddy simulation of turbulent open-channel flow with free surface simulated by level set method[J]. Physics of Fluids, 2005, 17: 1-12.

[22]Keating A, Piomelli U. A dynamic stochastic forcing method as

a wall-layer model for large-eddy simulation[J]. Journal of Turbulence, 2006, 7(12): 1-24.

[23] Breuer M, Jovicˇic＇ N, Mazaev K. Comparison of DES, RANS and LES for the separated flow around a flat plate at high incidence[J]. International Journal for Numerical Method in Fluids, 2003, 41: 357-388.

[24]Deck S. Numerical simulation of transonic buffet over a supercritical airfoil[J]. AIAA J., 2005, 43(7): 1556-1566.

[25]Deck S. Zonal-detached eddy simulation of the flow around a high-lift configuration[J]. AIAA J., 2005, 43(11): 2372-2384.

[26]Deck S, Weiss P, Pamiès M, et al. Zonal detached eddy simulation of a spatially developing plate turbulent boundary layer[J]. Computers & Fluids, 2011, 48: 1-15.

[27]Casulli V, Cattani E. Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow[J]. Computers and Mathematics with Applications, 1994, 27(4): 99-112.

[28]Stansby P K. Semi-implicit finite volume shallow-water flow and solute transport solver with κ-ε turbulence model[J]. International Journal for Numerical Methods in Fluids, 1997, 25: 285-313.

[29]Huang W, Spaulding M. Reducing horizontal diffusion errors in σ -coordinate coastal ocean models with a second order Lagrangian-interpolation finite-difference scheme[J]. Ocean Engineering, 2002, 29: 495-512.

[30]Zhang J X, Liu H, Xue L P. A vertical 2-D mathematical model for hydrodynamic flows with free surface inσcoordinate[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(1): 82-90.

[31]Darwish M S, Moukalled F. TVD schemes for unstructured grids[J]. International Journal of heat and Mass Transfer, 2003, 46: 599-611.

[32]Rhie C M, Chow W L. Numerical study of the turbulent flow past an airfoil with trailing edge separation[J]. AIAA J, 1983, 21: 1525-1532.

[33]Balachandar R, Polatel C, Hyun B -S, et al. LDV, PIV and LES investigation of flow over a fixed dune[M]//Gyr A, Kinzelbach W. Sedientation and Sediment Transport, Procedings of the symposium held in Monte Verita, Switzerland, Kluwer Academic, Dordrecht, 2002: 171-178.

[34]Weber L J, Schumate E D, Mawer N. Experiments on flow at a 90° open-channel junction[J]. Journal of Hydraulic Engineering, 2001, 127(5): 340-350.

0258-1825(2016)02-0239-13

帶自由表面水流的分離渦模擬

張景新1,2,3,*, 梁東方1,3, 劉 樺1,3

(1. 上海交通大學 水動力學教育部重點實驗室, 上海 200240; 2. 上海河口海岸科學研究中心 河口海岸交通行業(yè)重點實驗室, 上海 201201; 3. 上海交通大學 船舶海洋與建筑工程學院, 上海 200240)

目前水動力學問題的數(shù)值求解仍然廣泛采用RANS模型,特別對于大尺度的地表水流的數(shù)值模擬，該類模型計算效率較高，但其難以給出湍流場的小尺度渦結(jié)構(gòu)。深入研究復雜地形條件下的水流運動，需要發(fā)展諸如LES等的高級數(shù)值模型。天然地表水流，地形、邊界復雜，雷諾數(shù)通常較高，高級數(shù)值模擬受限于計算能力的限制，還難以工程應(yīng)用。RANS和LES混合模型是目前具有工程應(yīng)用的一種混合數(shù)值模型，其已經(jīng)在CFD領(lǐng)域獲得了長足發(fā)展。然而對于地表水流運動的數(shù)值模擬，相關(guān)工作還較少。本文將分離渦模型(DES)，即一種RANS和LES的混合模型，應(yīng)用于帶自由表面的地表水流運動，建立了一套數(shù)值仿真模型。模型基于有限體積法，水平面內(nèi)采用非結(jié)構(gòu)計算網(wǎng)格，垂向為結(jié)構(gòu)化網(wǎng)格，對流項離散格式采用二階TVD格式，并行基于OpenMP語言庫。算例表明DES模型有助于揭示復雜地形條件下帶自由表面水流的大渦擬序結(jié)構(gòu)。

分離渦;擬序結(jié)構(gòu);自由表面流動

O352

A doi: 10.7638/kqdlxxb-2016.0006

*Associate Professor, Department of Engineering Mechanics; zhangjingxin@sjtu.edu.cn

format: Zhang J X, Liang D F, Liu H. A Detached Eddy Simulation model for free surface flows[J]. Acta Aerodynamica Sinica, 2016, 34(2): 239-251.

10.7638/kqdlxxb-2016.0006. 張景新,梁東方,劉樺. 帶自由表面水流的分離渦模擬(英文)[J]. 空氣動力學學報, 2016, 34(2): 239-251.

Received： 2015-12-15; Revised：2016-02-10

Jointly supported by the National Natural Science Foundation (No.11572196), the National Basic Research Program of China (No. 2014CB046200), and the Non-Profit Industry Financial Program of MWR (201401027).