HYBRID SCHEME FOR COMPRESSIBLE TURBULENT FLOW AROUND CURVED SURFACE BODY

Xu Changyue,Ran Qian,Sun Jianhong

(College of Aerospace Engineering,NUAA,29 Yudao Street,210016,Nanjing,P.R.China)

INTRODUCTION

The prediction of compressible turbulent flow around a body with curved surface plays an important role in applications and fundamentals,especially theflow involving shock wave.To capture the discontinuity caused by shock wave,a dissipative shock-capturing scheme is applied,such as the total variation diminishing(TVD)schemes.For computing fine-scale feature,such as turbulence,away from shock wave,a non-dissipative numerical scheme must be used.These characteristics can prevent smearing of flow features and non-physical oscillations near the high gradient regions of flow variables.

A great effort has been made in the past decades to study high-order accurate shock-capturing scheme,e.g.,essentially nonoscillatory(ENO), weighted essentially nonoscillatory(WENO)and compact schemes[1-5].Recently,the high-order accurate hybrid schemes composed of central difference and high-order accurate shock-capturing sub-schemes have been investigated. A high-order accurate hybrid central-WENO scheme has been proposed in Ref.[6].The fifth order WENO scheme was divided into two parts,a central flux part and a numerical dissipation part,which were hybridized by means of a weighting function indicated the local smoothness of flow fields.Five flow problems were chosen for validation including interaction of a moving shock with a density wave,advection of an isentropic vortex,double Mach reflection of a strong shock,mixing-layer/shock interaction and weak-shock/vortex interaction.Shen and Yang[7]proposed a hybrid compact-WENO schemes,which were hybridized with compact difference and WENO scheme to draw lessons from the ENO scheme.Several typical cases,similar as those calculated by Kim and Kwon[6],were chosen for validating the proposed hybrid scheme.Large-eddy simulations of the Richtmyer-Meshkov instability with reshock have been performed in Ref.[8]using an improved version of the tuned centre-difference and WENO hybrid method.

Many researches on the high-order accurate hybrid scheme are applied to the compressible turbulent flow problems with simple geometries.However,many practical flows commonly take place on a curved surface. The finite volume method is widely used for compressible turbulent flow problems with curved surface and non-linear phenomena,such as compressible flow swirling flows injected into a coaxial dump chamber[9],transonic flows around an aerofoil[10]and a cylinder[11-12],etc.In this paper,a second-order accuracy Roe′s flux difference splitting scheme in the finite volume method is hybridized with the central scheme by means of a binary sensor function.Three typical cases arecalculated using detachededdy simulation(DES)technique,including the transonic turbulent flow over an axisymmetric bump,supersonic turbulent flow around an axisymmetric slender pointed body,and transonic turbulent flow past an aerofoil.

1 MATHEMATICAL FORMULATION AND NUMERICAL METHODS

1.1 Governing equations and boundary conditions

To investigate the compressible turbulent flow around a curved surface body,the 3-D Favre-averaged compressible Navier-Stokes equations in generalized coordinates are used.To nondimensionalize the equations,we use the freestream variables,including density d∞,temperature T∞ ,velocity U∞ and diameter of the body D for bump and slender body or cord length c for aerofoil as the characteristic scales.The governing equations in flux form can be written as

where Q=[d,d u,d v,d w,E]Tis the conservative flow variables with density d,velocity component ui,pressure p,and specific total energy E.J=?(a,Z,Y)/?(x,y,z)represents Jacobian of the transformation, (F,G,H)the inviscid flux terms,and(Fν,Gν,Hν)the viscous flux terms.

The diffusive fluxes are given by

where Re is the Reynolds number,_the molecular viscosity,Pr the Prandtl number and S ij the strain-rate tensor defined as Sij=(?ui/?xj+ ?uj/?xi)/2.Perfect gas relationship and Sutherland law for molecular viscosity coefficient_are employed.

In this study,the initial condition is set as the free-stream quantities.The far field boundary conditions are treated by local 1-D Riemann-invariants.No-slip and adiabatic conditions are applied to the body surface.

1.2 Turbulencemodeling

Three typical cases are calculated,including transonic turbulent flow over an axisymmetric bump,supersonic turbulent flow around an ax-isymmetric slender pointed body,and transonic turbulent flow past an aerofoil. As the high Reynolds number flows are considered,say about 106,the methodology used in this paper is DES introduced by Spalart[13].The model is derived from Spalart-Allmaras turbulence model[14],which is aone-equation model for theeddy viscos-solving a transport equation.The reader may refer to Ref.[14]for details on theconstants and the quantities involved.

The model is provided with a destruction term for the eddy viscosity that depends on the distance to the nearest solid wall d.This term adjusts the eddyscale with local deformation rateproducing an eddy viscosity givenSpalart[13]proposed to replace d to the closest wall withdefined by

whereΔrepresents a characteristic mesh length and is defined as the largest of grid spacing in all three directions,i.e., Δ=max(Δx,Δy,Δz),and the constant C DES is taken as 0.65 from a calibration of the model for isotropic turbulence[15].When d ＜ C DESΔ,themodel acts in a Reynolds averaged Navier-Stokes(RANS)mode and when d＞ CDESΔ the model acts in a Smagorinsky largeeddy simulation(LES)mode.

1.3 Numerical procedure

The governing equations are numerically solved by the finite-volume method.Both the convective and diffusive terms arediscretized with a second-order central schemes,and a fourth-order low artificial numerical dissipation is employed to prevent the numerical oscillations at high wavenumbers[11,16]. The temporal integration is performed using an implicit approximatefactorization method with sub-iterations to ensure the second-order accuracy.

綜上所述，建筑工程造價(jià)合理預(yù)算對(duì)建筑工程項(xiàng)目管理和質(zhì)量安全有重要作用，建筑企業(yè)應(yīng)該高度重視工程造價(jià)預(yù)算環(huán)節(jié)，分析在工作過程中出現(xiàn)超預(yù)算現(xiàn)象的原因，并制定行之有效的解決措施。培養(yǎng)專業(yè)的預(yù)算管理人員，分工明確，避免出現(xiàn)工程造價(jià)超預(yù)算現(xiàn)象，從而提高建筑企業(yè)的經(jīng)濟(jì)效益和市場(chǎng)中的競(jìng)爭(zhēng)力。

To capture the discontinuity caused by a shock wave,a second-order upwind scheme with Roe′s flux difference splitting is introduced into theinviscid flux.The spatial discretization is constructed explicitly to be shock capturing with the upwind scheme and to revert to a central stencil with low numerical dissipation in turbulent flow regions away from shock wave.The artificial dissipation is also turned off in the region in which the upwind scheme works.A binary sensor function H i+1/2 at cell interface i+ 1/2 is used for the detection of shock waves.H i+1/2 is determined by the pressure and density curvature criteria proposed in Ref.[8].wheredenotes the complement of l,Ti+p1/2and Td i+1/2represent the pressure and density relative curvatures at cell interface i+1/2.

The 3-D version of this detection is used in the simulations.Similar to the treatment[8],the value of c1 and c2 that proved to give the best results arechosen as 0.01.Based on this detection,Roe′s second-order upwind flux only operates at the cells in the vicinity of shock waves.

To illustrate the hybrid technique,the flux difference in the i direction(?a F)i can be written as

whereκ(4)andλi+1/2 are the coefficient and spectral radius of 4th Jameson′s artificial numerical dissipation[16],respectively. Δ(? )=(? )i+1-(? )iis the forward difference operator,▽ (? )=(? )i-(?)i-1 the backward difference operator,q=[d,u,v,w,p]Tthe non-conservative flow variable.q L and q R are theinterpolated values in the cell-interface position.And A conv denotes the inviscid ma-trix in Roe′s flux-difference splitting scheme[17].

where the diagonal matrixΛis the matrix of eigenvalues of A conv,T the matrix of right eigenvectors as columns,and T-1the matrix of left eigenvectors as rows.

When H i+1/2 is 1,thehybrid scheme Eq.(11)reverts to a Roe′s FDSscheme

However,when H i+1/2 equals zero,the hybrid scheme reverts to a central scheme with 4th Jameson′s artificial numerical dissipation

2 RESULTSAND DISCUSSION

2.1 Transonic turbulent flow over axisymmetric bump

Transonic flow over the bachalo-Johnson[18]axisymmetric bump is chosen to study the performance of the present hybrid scheme in predicting the compressible flow involving shock wave and turbulence. The free-stream Mach number is 0.875 and Reynolds number based on the pipediameter D is 2.66×106.Fig.1shows thegeometry and computational mesh of the bump.Here,the diameter of bump is 1.25D.O-H mesh is employed,and the grid number is 113× 81× 65 in the streamwise,normal,azimuthal direction,respectively,with the outer boundary 30D.Time step is 0.005D/a∞.

Fig.1 Geometry and computational mesh of bump

hybrid scheme.Here,the increment between the iso-lines is a constant value of 0.04.A local supersonic zone is formed above the bump surface,and a shock wave occurs at x/D≈0.65.Moreover,the shock wave is captured well,which indicates that the present hybrid scheme is capable of capturing the discontinuity caused by shock wave.

Fig.2 Time-averaged flow pattern by means of isolines of local Mach number

The measured and computational mean pressure coefficient distribution〈Cp〉 on the surface of the cylinder and bump are show n in Fig.3.Compared with experimental measurement[18],the upwind scheme overestimates the position of shockwave/turbulence interaction.The hybrid scheme gives a better agreement with the experimental data,which is closely related to the low-dissipation of the present hybrid scheme.The mean streamwise velocity〈u〉 distributions calculated by the two numerical schemes are compared with experimental data[17]at thevarious streamwise locations from x/D=0.563 to x/D=1.375,as shown in Fig.4.Upstream shock wave(x/D=0.563),the computed result obtained by the hybrid schemeis same as the upwind scheme.How-ever,downstream shock wave(x/D=0.563,0.875,1.000and 1.375),results obtained by the hybrid scheme are closer to the experimental data, which indicates that the present hybrid scheme can predict well the shock wave and turbulence over the bump.

2.2 Supersonic turbulent f low past axisymmetric slender pointed body

In this case,supersonic turbulent flow involving the shock wave around an axisymmetric

Fig.3 Distribution of wall pressure coefficient〈C p〉

Fig.4 Mean streamwisevelocity profiles obtained numerically and experimentally

slender pointed body is investigated. The body geometric equation can be written as

Fig.5 shows thegeometry and computational mesh of the axisymmetric slender pointed body.O-H-type grid is used with clustered distributions in the vicinity of the wall.The grid number is 113× 101× 81 in the streamwise,radial, azimuthal direction,respectively,with the outer boundary 30D,and the time step is 0.005D/a∞.According to the previous experiment[19], the free-stream Mach number Ma∞is chosen as 2.5,angle of attack T as 14o,and Reynolds number Re based on thediameter of thecylinder as 1.1×106.

Time-averaged flow pattern in the meridiansection calculated by the hybrid scheme is shown in Fig.6 using local Mach number iso-lines.Here,the increment between the iso-lines is a constant value of 0.04.We can found that an oblique

Fig.5 Geometry and computational mesh of axisymmetric slender pointed body

shock wave occurs at the windward surface,and some expanded waves exist on the leeward surface.Theazimuthal distributions of wall pressure coefficient〈C p〉 obtained numerically and experimentally are exhibited in Fig.7.Here,θ=0°is corresponding to the windward surface. At x/D=3.5,two calculated results agree well with the experimental data,as shown in Fig.7(a).Fig.7(b)shows that the mean separation position predicted by the upwind scheme is delayed to the leeward surface at x/D=5.5,however,the hybrid scheme can give an improved result.Moreover,〈C p〉 is overestimated by both the upwind and hybrid schemes on the leeward surfaceθ≈150°.Fig.7(c)shows that〈C p〉 is underestimated by the two numerical schemes on the leeward surface θ≈ 150°. Similar to 〈C p〉 distribution at x/D=5.5,the mean separation position predicted by theupwind scheme is delayed to the leeward surface at x/D=11.5,as shown in Fig.7(d).Furthermore,after we carefully examine the flow fields,it is revealed that the hybrid scheme can predict the main vortex bifurcating and the secondary vortex structures,which also indicatesthat the present hybrid scheme can give a reasonable result.

Fig.6 Time-averaged flow pattern in meridian-section using iso-lines of local Mach number

Fig.7 Azimuthal distributions of mean wall pressure coefficient〈C p〉 obtained numerically and experimentally

2.3 Transonic turbulent flow past aerofoil

Finally,another typical case,transonic turbulent flow past a RAE2822 aerofoil,is chosen to validate the present hybrid scheme.Experiment in Ref.[20]is considered.Here,the experimental data is obtained in the wind tunnel for the flow conditions:Ma∞=0.75,T=3.19°,Re=6.2× 106.

In order to compare the computational flow with experimental flow past the aerofoil in freeflight conditions,some corrections to the tunnel data are required.The wind tunnel correction proposed in Ref.[21]is used,i.e.,Ma∞=0.75,T=2.72o,Re=6.2×106.C-type mesh is used,and the grid number is 257×97,as shown in Fig.8.The far field is about 20 cord lengths from the aerofoil,and the time step is 0.005D/a∞.

Fig.9 shows time-averaged flow pattern calculated by the hybrid scheme using the iso-lines

Fig.8 Computational mesh of RAE2822 aerofoil

of local Mach number.Here,the increment between the iso-lines is a constant value of 0.04.Similar to the bump case,local supersonic zone and shock wave are formed over the aerofoil surface,suggesting that the present hybrid scheme can capture the shock wave.To show the performance of the present hybrid scheme in predicting the compressible turbulent flow involving shock wave quantitatively,comparisons of the calculated wall pressure coefficient〈Cp〉 and skin friction coefficient〈C f〉 distributions with the experimental measurement are given in Fig.10.The mean shock wave position is predicted downstream by the upwind scheme.However,result obtained by the hybrid scheme agrees well with the experimental data.The calculated〈C f〉 distributions obtained by the two numerical schemes give a good agreement with the experimental data upstream the shock wave.However,〈C f〉 is slightly underestimated by both upwind and hybrid schemes downstream the shock wave.For more validation,comparisons of thecomputational time-averaged lift coefficient 〈CL〉tand drag coefficient〈C D〉t with the experimental ones are exhibited in Table 1,which also shows a good agreement between the numerical results obtained by the hybrid scheme and the experimental measurements.

Fig.9 Time-averaged flow pattern using iso-lines of local Mach number

Fig.10 Comparison of present numerical results with experimental data

Table 1 Comparisons of present computational results with experimental data in Ref.[20]

3 CONCLUSION

A hybrid central-upwind scheme for the compressible turbulent flow around a curved surface body is proposed.Two sub-schemes,the central difference scheme and the Roe′s flux-difference splitting scheme,arehybridized by means of a binary sensor function.In order to examine the capability of the hybrid scheme in computing the compressible turbulent flow around a curved surface body,especially the flow involving shock wave,three typical cases are investigated using detached-eddy simulation technique,i.e.,the compressible flow including transonic and supersonic regime over an axisymmetric bump,axisymmetric slender pointed body and aerofoil.It is revealed that the numerical results obtained by the proposed hybrid scheme can give more reasonable agreements with the experiment than that obtained by the upwind scheme,which indicates that the present hybrid scheme can be applied to computing the compressible flow around a curved surface body involving shock wave and turbulence.

[1] Shu C W,Osher S.Efficient implementation of essentially non-oscillatory shock capturing schemes II[J].Journal of Computational Physics,1989,83(1):32-78.

[2] Jiang G S,Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics,1996,126(1):202-228.

[3] Balsara D S,Shu C W.Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[J].Journal of Computational Physics,2000,160(2):405-452.

[4] Tolstykh A I.High accuracy non-centered compact difference schemes for fluid dynamics applications[M].Singapore:World Scientific,1994.

[5] Lele S K.Compact finite difference schemes with spectral-like resolution[J].Journal of Computational Physics,1992,103(1):16-42.

[6] Kim D,Kwon J H.A high-order accurate hybrid scheme using a central flux scheme and a WENO scheme for compressible flowfield analysis[J].Journal of Computational Physics,2005,210(2):554-583.

[7] Shen Y,Yang G.Hybrid finite compact-WENO schemes for shock calculation [J]. International Journal for Numerical Methods in Fluids,2007,53(4):531-560.

[8] Hill D J,Pantano C,Pullin D I.Large-eddy simulation and multiscale modeling of a Richtmyer-Meshkov instability with reshock[J].Journal of Fluid Mechanics,2006,557(1):29-61.

[9] Lu X Y,Wang S W,Sung H G,et al.Large eddy simulations of turbulent swirling flows injected into a dump chamber[J].Journal of Fluid Mechanics,2005,527(1):171-195.

[10]Chen L W,Xu C Y,Lu X Y.Numerical investigation of the compressible flow past an aerofoil[J].Modern Physics Letters B,2009,23:233-236.

[11]Xu CY,Chen L W,Lu X Y.Large-eddy simulation of the compressible past a wavy cylinder[J].Journal of Fluid M echanics,2010,665(1):238-273.

[12]Xu C Y,Chen L W,Lu X Y.Effect of Mach number on transonic flow past a circular cylinder[J].Chinese Science Bulletin,2009,54(11):1886-1893.[13]Spalart P R,Jou W H,Strelets M,et al.Comments on the feasibility of LESfor wings and on a hybrid RANS/LES approach [C]//Advances in DNS/LES.Columbus,OH:Greyden Press,1997.

[14]Spalart P R,Allmaras S R.A one-equation turbulence model for aerodynamic flows[R].AIAA paper 92-0439,1992.

[15]Shur M L,Spalart P R.Detached-eddy simulation of an airfoil at high angle of attack[C]//the 4th International Symposium of Engineering Turbulence Modelling and Measurements.Amsterdam:Elsevier,1999:669-678.

[16]Swanson R C,Turkel E.On central-difference and upwind schemes [J]. Journal of Computational Physics,1990,101(2):292-306.

[17]Bachalo W D,Johnson D A.Transonic,turbulent boundary-layer separation generated on an axisymmetric flow model[J].AIAA Journal,1986,24(3):437-443.

[18]Xu Changyue.Large eddy simulation of the compressible flow past a circular cylinder and its flow control[D].Hefei:University of Science and Technology of China,2009.(in Chinese)

[19]Yang Yunjun,Zhou Weijiang.Numerical simulation of turbulence flows about a slender pointed body at high angle of attack[J].Acta Aerodynamica Sinica,2003,21(3):351-355.(in Chinese)

[20]Cook P H,Mcdonald M A,Firmin M C P.Aerofoil 2822-pressure distributions, boundary layer and wake measurement[R].AGARD AR 138,1979.

[21]Krist S,Biedron R,Rumsey C.CFL3D user′s manual(Version 5.0)[R].NASA/TM-1998-208444,1998.