中 華 大 學
碩 士 論 文
題目：Direct Numerical Simulation of a Turbulent Flow
over a Wavy Channel
(直接數值模擬紊流流經波形流道之流場分析)
系 所 別：機械與航太工程研究所 學號姓名：M09408003 周 中 祺
指導教授：蔡 永 培 博 士
中華民國 九十六 年 七 月
中文摘要
波形流道流場是自然界中相當重要的一種流動現象，它在工程與科 學領域上有相當重要的應用與重要的地位。亦可應用於熱傳導方面，波 形壁有助於表面積的增加，可利加快散熱速度。再者亦可發展溝槽薄 膜，法國空中巴士A320試驗機將其貼覆於飛機70%表面達到節省1~2%
之燃料之功效。因此針對該特性進行數值模擬研究，在數值方法的選擇 上，以加權型基本不震盪算則(Weighted Essentially Non-Oscillatory)直接 求解奈維爾-史托克 (Navier-Stokes)方程式。而程式撰寫技巧上採用 Lower-upper Symmetric successive over relaxation ( LU-SSOR method ) 以增強數值之穩定性及加速程式之收斂，並與Ts a i [ 1 ] , K i m[ 2 ] 及 C h e r u k a t e t a l . [ 3 ] 等 文 獻 進行模擬值比較。於前階段結果發現 WENO方法解析與比較的文章有相當不錯的精確度，均在波形壁的下坡 處準確的解析出流體分離及再重合中所產生的迴流區。進而在將振幅加 倍進行模擬得到流體於波形壁下波處的分離點不僅會提早發生，迴流區 的範圍也會加大的結果。故本研究以此為基準進而針對波形壁流動現象 進行速度場；幾何外型二大部分進行系統探討。
中文關鍵詞：波形壁；超音速；加權型基本不震盪算則。
Abstract
The flows over a wavy wall are very important phenomena in nature. They play quite important roles in engineering and scientific applications. Since wavy wall can increase boundary’s surface and speed up radiation, therefore wavy wall flows can be applied in the area of heat conduction. Moreover a concave-convex thin film coating was developed and applied on the external surface of an Airbus A320 test airplane.
The experiment results indicate that a surface coverage of 70% could achieve a fuel savings of 1-2%. In the direct simulation methods of flows over a wavy wall, Weighted Essentially Non-Oscillatory by the compressible flow domain was selected to solve Navier-Stokes equation. And the LU-SSOR method was applied to increase numerical stability and converging of the formula. A good agreement was found between simulation analysis results based on the WENO method and the experimental values generated by Tsai[1] and Kim[2] and Cherukat et al.[3]. The simulation analysis has accurately predicted the fluid separation and superposition of the return flow in the wavy wall downhill areas. As the amplitude doubles, it was predicted that the detachment point of the fluid can occur earlier near the profile wall and the region of return flow can increase. This research then carries on the velocity field take this as the datum in view of the wavy wall phenomenon of flow, geometry outlook. The wavy wall applies major part three carries on the system discussion.
Keyword：Wavy wall, Supersonic, WENO scheme.
Acknowledgements
The work reported in this master's degree describes my research activities at Chung Hua University during the period from September 2005 to July 2007. I now wish to thank those who provided me with help and encouragement during the present study.
First and foremost, I wish to express my sincere gratitude to my advisor Professor Dr.
Y. P. Tsai for this valuable advice and encouragement during the course of this study.
Without his patience and careful guidance this work would not be accomplished. I am indebted to the members of my oral defence committee, Dr. H. S. Huang and Dr. T. S.
Cheng for their constructive criticism. Besides, my sincere thanks also go to the scholarships supported by the CFD lab partner, James p. Oyan, C. Y. Lin et al., and Litter Tiger, Y. R. Lin.
Finally, and most importantly, I would like to dedicate this work to my dear Phoebe for her support and understanding. I wish to thank my parents for all their encouragement.
They gave much more than they expected, but they contributed more than they can imagine.
Contents
中文摘要 · · · 1
Abstract · · · 2
Acknowledgements · · · 3
Contents · · · 4
Figure contents · · · 5
1. Introduction · · · 8
2. Numerical method · · · 11
2.1 Flow configuration · · · 11 2.2 Governing equation · · · 1 2 2.3 Spatial Discretization· · · 1 3 2.4 Temporal Discretization · · · 1 4 2.5 Computational Grid · · · 1 6 2.6 Initial condition and boundary condition · · · 1 7 2.6.1 Initial condition · · · 1 7 2.6.2 Boundary condition · · · 1 7 2.6.3 Boundary condition on top of base and below · · · 1 7 2.6.4 Inflow and outflow boundary condition · · · 1 8 3. Results · · · 1 9 3.1 The flows over a wavy wall analyze · · · 1 9 3.1.1. Vortex change· · · 1 9 3.1.2. Pressure change· · · 2 0 3.2 Wavy wall system · · · 2 1 3.2.1 The subsonic and supersonic flow over a solid wavy wall · · · 2 1 3.2.2 The to compare of subsonic and supersonic flow over a solid wavy wall · · · 2 3 3.2.3 By physical space comparison · · · 2 4 4. Conclusions · · · 2 6 References · · · 2 7 Figure · · · 2 8
Figure contents
Fig.2-1 Flow over a wavy wall of a computational domain · · · 11 Fig.2-2 A grid arrangement uses grid number of 251×41· · · 1 6 Fig.2-3 Computational Grid resolutions · · · 1 7 Fig.A The vortex in flows over a wavy wall · · · 1 9 Fig.B Pressure change of flows over a bilateral wavy wall · · · 2 0 Fig.3-1 Computational Grid · · · 2 8 Fig.3-2 Streamwise velocity contour at T=2.3msec · · · 2 8 Fig.3-3 Streamwise velocity contour · · · 2 9 Fig.3-4 Spanwise velocity contour for M=0.6 at T=2.3msec · · · 2 9 Fig.3-5 Comparison with streamwise velocity contour for M=0.6 at T=2.3 msec. · · · 3 0 Fig.3-6 Comparison with spenwise velocity contour for M=0.6 at T=2.3 msec · · · 3 0 Fig.3-7 Streamwise velocity and spanwise velocity contour · · · 3 1 Fig.3-8 Comparison with crest and trough r.m.s · · · 3 2 Fig.3-9 Comparison with pressure contour for M=0.6 at T=2.3 msec · · · 3 3 Fig.3-10 Supersonic single wavy wall time T=2.3msec in the density distribute
contours · · · 3 4 Fig.3-11 The point streamline for much number M=1.5, A/λ=0.1· · · 3 5 Fig.3-12 Reynolds shear stress (normalized by )x/λ=0.0 · · · 3 5 Fig.3-13 The point density distribute for much number M=1.5, A/λ=0.1 · · · 3 6 Fig.3-14 The streamline for much number M=1.5, A/λ=0.1 · · · 3 6 Fig.3-15 The point vector for much number M=1.5, A/λ=0.1 · · · 3 7 Fig.3-16 Mean velocity profiles at crest and trough streamline locations · · · 3 8 Fig.3-17 Reynolds shear stress (normalized by
2
Ub
)x/λ=0.0 · · · 3 9 Fig.3-18 Reynolds shear stress (normalized by
2
Vb
)x/λ=0.0 · · · 4 0 Fig.3-19 Reynolds shear stress (normalized by
2
Vb
)x/λ=0.0 · · · 4 1 Fig.3-20 The point contours for subsonic time T=1.5msec in the mean streamwise
and normal velocity · · · 4 2 Fig.3-21 The point contours for subsonic time T=2.3 msec in the mean streamwise
and normal velocity · · · 4 2 Fig.3-22 The point contours for supersonic time T=2.3 msec in the mean streamwise
and normal velocity · · · 4 2 Fig.3-23 Shows subsonic instantaneous velocity chart for time T=2.3msec in the
spanwise direction · · · 4 2 Fig.3-24 Shows supersonic instantaneous velocity chart for time T=2.3msec in the
2
Ub
cross streamwise direction · · · 4 3 Fig.3-25 The point streamline for much number M=0.6, time T=2.3msec, amplitude
A=2.54mm · · · 4 3 Fig.3-26 The point streamline for much number M=1.5, time T=2.3msec, amplitude
A=2.54mm · · · 4 4 Fig.3-27 Pressure curve from the studies for much number M=0.6 and 1.5 time
T=2.3msec· · · 4 4 Fig.3-28. Comparison the turbulence intensities also agree reasonably well with
subsonic and supersonic · · · 4 5 Fig.3-29 Computational domain of flow over a single wavy wall · · · 4 6 Fig.3-30 Computational domain of flow over a bilateral wavy wall · · · 4 6 Fig.3-31 The computational grid of single wavy wall oscillation amplitude
A=2.54(mm) · · · 4 6 Fig.3-32 The computational grid of single wavy wall oscillation amplitude
A=5.08(mm)· · · 4 7 Fig. 3-33 The computational grid of bilateral wavy wall oscillation amplitude
A=5.08(mm)· · · 4 7 Fig. 3-34 The computational grid of flat oscillation amplitude A=0(mm) · · · 4 7 Fig. 3-35 The point one single wavelength for amplitude A=5.08, T=0.5msec
streamwise velocity· · · 4 8 Fig. 3-36 The point one single wavelength for amplitude A=5.08, T=1.5msec
streamwise velocity· · · 4 8 Fig. 3-37 The point one single wavelength for amplitude A=5.08, T=1.5msec
spanwise velocity· · · 4 9 Fig. 3-38 The point streamline for much number M=0.6, time T=1.5msec, amplitude
A=2.54mm · · · 4 9 Fig. 3-39 The point streamline for much number M=0.6, time T=2.3msec, amplitude
A=2.54mm · · · 5 0 Fig. 3-40 The point streamline for much number M=0.6, time T=1.5msec, amplitude
A=5.08mm · · · 5 0 Fig. 3-41 Shows the speed contours from Cherukat et al.[3] for varying x/λ from 0.0
to 1.0 · · · 5 1 Fig.3-42 Shows the speed contours for much number M=0.6, time T=2.3msec, x/λ= 0.05
· · · 5 1
Fig. 3-44 Point contours for supersonic single wavy wall time T=2.3msec in the
density distribute contours· · · 5 2 Fig. 3-45 Point contours for supersonic bilateral wavy wall time T=2.3msec in the
density distribute contours· · · 5 2 Fig. 3-46 Point contours for supersonic flat boundary time T=2.3msec in the density
distribute contours · · · 5 2 Fig. 3-47 Point contours for supersonic single wavy wall time T=2.3msec in the mean
streamwise and normal velocity · · · 5 2 Fig. 3-48 Point contours for supersonic bilateral wavy wall time T=2.3msec in the
mean streamwise and normal velocity· · · 5 3 Fig. 3-49 Point contours for supersonic flat boundary time T=2.3msec in the mean
streamwise and normal velocity · · · 5 3 Fig. 3-50 Shows supersonic single wavy wall instantaneous velocity chart for time
T=2.3msec in the cross streamwise direction· · · 5 3 Fig. 3-51 Shows supersonic bilateral wavy wall instantaneous velocity chart for time
T=2.3msec in the cross streamwise direction· · · 5 3 Fig. 3-52 Shows supersonic flat boundary instantaneous velocity chart for time
T=2.3msec in the cross streamwise direction· · · 5 4 Fig. 3-53 The point streamline for single wavy wall, much number M=1.5, time
T=2.3msec, amplitude A=5.08mm · · · 5 4 Fig. 3-54 The point streamline for bilateral wavy wall · · · 5 4 Fig. 3-55 The point streamline for flat boundary · · · 5 4 Fig. 3-56 Pressure curve from the studies for much number M=1.5, T=2.3msec· · · 5 5
1. Introduction
Recent researches on wavy channel concentrate on the areas of revelation of fluxion structure and the control and reduction of anti-obstruction mechanism. The research was aimed to theorize this phenomenon and use it in the engineering applications to reduce obstruction. However the effects of fluid structure change and the wall surface pressure gradient on the change of wall surface obstruction in the near wall area are very complex. The reduction in obstruction from the wavy channel is apparently due to the improvement of the partial turbulent flow structure in the near wall areas, which has reduced the friction in the fluids and between the fluid the solid wall. Therefore in order to further the revelation of wavy channel in obstruction reduction, it is believed that the near wall flowing should be characterized.
Experimental research of flow over the wavy wall was performed as early as 1932.
The theory of the surface wall variation has been developed. Zilker and Hanratty[4], Zilker et al.[5] continued to conduct studies including the experiments of gauging wall surface pressure, measuring mean velocity of stream and the wall surface shearing stress and analyzing the mobile characteristics. In 1985 Patel[6] has improved simulation for several kinds of near wall with low Reynolds number kmodel.
Patel and Chen[7] have used double-decked model to solve turbulent complex stalled flow. Since this model saves the grids, therefore it can effectively save the computation space and time and enhances the computation feasibility. Patel et al.[8] calculation of the wavy wall turbulent flow was conducted on two kinds of wave height ratios in 1991. The simulated streamlines picture in the near wall vortex area and the cross section velocity distribution were in a good agreement with the experimental results.
Comparison of friction coefficient curves with pressure coefficients was also made between the wavy wall surface and flat surface. The effects of the profile on flows in the near wall area were analyzed. Ferrira and Lopes [9] carried on the multi-group wind tunnel experiments on the unimodular flows of sinusoidal wavy wall flow field.
The low Reynolds number k-εmodel with the control volumetric method was used to compute wavy wall near zone stalled flow for various wave height ratios (A/). Montalbano and McCready [10] used the wave stability theory and added small perturbation quantity to Orr-Sommerfeld equation to develop the relation between the laminar wall surface pressure and the shear stresses.
Airiau and Giovannini [11] utilized statistical simulation to obtain the stream function and vortex charts at different time for the sinusoidal wavy wall. Also the average stream function and the average wall surface pressure coefficient curves were calculated using the time average method. With these tools vorticity and variation of
Bontozoglou[12] the dimensionless Navier-Stokes equations are solved in the whole range of the laminar flow regime. Numerical predictions are compared with available experimental data for very low Reynolds numbers. The emphasis in the discussion of results is given in the presentation of free surface profiles, streamlines, velocity, and pressure distributions along the free surface and the wall. The interaction of the dimensionless numbers of the flow is studied, criteria for flow reversal are established, and a resonance phenomenon at high Reynolds numbers is investigated. Boersma[13].
The evolution in space and time of particles are released in this flow will be examined.
It will be shown that small waves on the channel bottom can generate large longitudinal vortices similar to Langmuir vortices that are observed in flows with waves at the free-surface. The simulation results show that the concentration of the particles is maximal on the downstream side of the wave crest. Nakagawa et al.[14], Measurements of turbulence with laser Doppler velocimetry (LDV) are compared for turbulent flow over a flat surface and a surface with sinusoidal waves of small wavelength. The wavy boundary was highly rough in that the flow separated. The Reynolds number based on the half-height of the channel and the bulk velocity was 46,000.
The wavelength was 5 mm and the height to wavelength ratio was 0.1. The root-mean-squares of the velocity fluctuations are approximately equal if normalized with the friction velocity. This can be explained as a consequence of the approximate equality of the correlation coefficients of the Reynolds shear stress. Calculations with a direct numerical simulation (DNS) are used to show that the fluid interacts with the wall in quite different ways for flat and wavy surfaces. They show similarity in that large quadrant 2 events in the outer flow, for both cases, are associated with plumes that emerge from the wall region and extend over large distances. Measurements of skewness of the streamwise and wall-normal velocity fluctuations and quadrant analyses of the Reynolds shear stresses are qualitatively similar for flat and wavy surfaces. However, the skewness magnitudes and the ratio of the quadrant 2 to quadrant 4 contributions are larger for the wavy surface. Thus, there is evidence that turbulent structures are universal in the outer flow and for quantitative differences in the statistics that reflect differences in the way in which the fluid interacts with the wall. Zilker et al.[15] measurements of the shear-stress variation along and the velocity profiles above a solid wavy wall bounding a turbulent flow are presented for waves with height-to-length ratios of 2a/λ=0.0312 and 0.05. These are compared with previous measurements of the wall shear stress reported by Thorsness (1975) and by Morrisroe(1970)for2a/λ = 0.012. The investigation covered a range of conditions from those for which a linear behaviour is observed to those for which a separated flow is just being initiated.
In the direct simulation methods of flows over a wavy wall, Weighted Essentially Non-Oscillatory by the compressible flow domain was selected to solve Navier-Stokes equation. And the LU-SSOR method was applied to increase numerical stability and converging of the formula. A good agreement was found between simulation analysis results based on the WENO method and the experimental values generated by Tsai[1]
“Direct Numerical Simulation of a Fully Developed Turbulent Flow over a Wavy Wall”and Kim[2] “Numerical investigation of confined supersonic mixing flow”and Cherukat et al.[3] “DirectNumericalSimulation ofaFullyDeveloped Turbulent Flow over a Wavy Wall” The simulation analysis has accurately predicted the fluid separation and superposition of the return flow in the wavy wall downhill areas. As the amplitude doubles, it was predicted that the detachment point of the fluid can occur earlier near the profile wall and the region of return flow can increase. This research then carries on the velocity field take this as the datum in view of the wavy wall phenomenon of flow, geometry outlook. The wavy wall applies major part three carries on the system discussion.
2. Numerical method
2.1 Flow configuration
The parameters of the simulation are the same as used by Tsai[1] and Cherukat et al.[3]. In their numerical, the amplitude and the wavelength of the wavy wall were 2.54 mm and 50.8 mm, respectively. The distance between the mean location of the wavy surface and the flat wall was 50.8 mm. A schematic diagram of the three-dimensional computational domain is shown in Figure 2-1. It consists of a channel which is unbounded in both the streamwise (x) and spanwise (z) directions. The lower wall has
) (4
N waves with sinusoidal shape and a mean position at the y = 0 plane (y is the vertical direction). The flat wall is located at y = h. The location of the wavy wall, y ,
is given by cos(2 )
a x
y where a is the amplitude of the wave and is the wavelength. The mean flow in the streamwise direction is pressure driven. In the present study, wavelength ¸ and amplitude a were set equal to h and 0:05h to match the parameters of Hudson’s (1993) measurements.The flow is assumed to be homogeneous in the spanwise direction, justifying the use of periodic boundary conditions. The flow is also assumed to be periodic in the streamwise direction. Thus, the computational box size in the streamwise (xN) and spanwise directions( x) should be large enough to include the largest length scale of the turbulent structures. The extents of the computational domain were, respectively, chosen to be 4h in the streamwise direction and 2h in the spanwise direction.
F i g . 2 - 1 F l o w o v e r a w a v y w a l l o f a c o m p u t a t i o n a l d o m a i n .
2.2 Governing equation
2D Navier-Stokes differential equation：
E E F F 0 Q
t x y
( 2 - 1 )
where Q is a conservation variable in equation：
Q u v e
( 2 - 2 a )
where is density; u and v are quantities of velocity; e is unit volume energy in the above equation. E and F are Non-stickiness fluxes︰
2
u u P
E uv
e P u
( 2 - 2 b )
2
v F uv
v P e P v
( 2 - 2 c )
where P is pressure which is related to conservation variable as follows:
1
2 2
/ 2P r eu v ( 2 - 3 )
E and F are stickiness fluxes defined as:
0
xx e yx
x
E M
R
( 2 - 4 a )
0
xy e yy
y
F M
R
( 2 - 4 b )
Where the shear stressand heat flux are:
2xx ux vy ux
( 2 - 5 a )
2yy ux vy vy
( 2 - 5 b )
xy yx ux vy
( 2 - 5 c )
1x xx xy
r
T u v
P x
( 2 - 6 a )
1y yx yy
r
T u v
P y
( 2 - 6 b )
where 2
3 and viscosity coefficient is：
1 c 1.5
T c T
( 2 - 7 )
where T is a non-dimensional temperature; c isSutherlan’sconstant110.4
T ; r p t
P c k
, generally pr is set at 0.72.
2.3 Spatial Discretization
The framework for constructing finite-volume WENO scheme is to calculate the fluxes to be used at the cell boundary in the methods. The numerical flux at the cell boundary xi1/2 ,Fi1/2, for example is defined as a montone function of left and right extrapolated values QiL1/2,
Q
iR1/2：) Q , Q F(
= (t)) F(Q
=
Fi+1/2 i+1/2 Li1/2 Ri1/2 ( 2 - 8 )
These extrapolated values are obtained from cell averages by means of a high order polynomialreconstruction such asgiven in Shi,Hu and Shu’spapers.Forexample,for a scalar function q(x) the fifth order WENO reconstruction defines the left extrapolated value qiL1/2 as
2 2 1 1 0 0 L
2 / 1
q
i
( 2 - 9 )where vk is the extrapolated value obtained from cell averages in the kth stencil
) 2 k i , 1 k i , k i (
S
k
( 2 - 1 0 )) 2 5 6(
1
1 2
0 qi qiqi
, ( 2 - 11 a )
) 2 5 6(
1
1 1
1 qiqi qi
, ( 2 - 11 b )
) 11 7
2 6( 1
1 2
2 qi qi qi
, ( 2 - 11 c )
and are nonlinear WENO weights given byk
2 1 0
0
0
( 2 - 1 2 a )
2 1 0
1
1
( 2 - 1 2 b )
2 1 0
2
2
, ( 2 - 1 2 c )
2 0 6
) 10
( 10
3 IPSO
( 2 - 1 3 a )
2 1 1 6
) 10
( 5
3
IPS
( 2 - 1 3 b )
2 2 2 6
) 10
( 10
1
IPS
( 2 - 1 3 c )
The smoothness indicators IPSk are
2 2 1 2
2 1
0 (3 4 )
4 ) 1 2
12( 13
qi qi qi qi qi qi
IPS , ( 2 - 1 4 a )
2 1 1 2
1 1
1 ( )
4 ) 1 2
12( 13
qi qi qi qi qi
IPS , ( 2 - 1 4 b )
2 1 2
2 1 2
2 ( 4 )
4 ) 1 2
12( 13
i i i
i i
i q q q q q
q
IPS ( 2 - 1 5 c )
Complete and detailed description about WENO scheme for solving 2D Navier-Stokes equation can be found in the papers by Casper, Atkins, Hu, Shu, Liu, Osher. and Chan.
2.4 Temporal Discretization
In the time integral, this study uses implicit expression time integration of LU-SSOR by Yoon & Jameson (1987) to the two orders precision. Also only the diagonal matrix of the concealed mode of the inverse matrix was calculated to save the computing time. The vectorization process is quite succinct and easy to understand.
Therefore the expression may apply only in the incompressible flow field.
Firstly, making linearization processing to the n+1 time step flux,
1 2 1 2 1 2 1 2
Q o Q B F F
Q o Q A E E
Q o Q B F F
Q o Q A E E
n v n v n v
n v n v n
v
n n n
n n n
( 2 - 1 6 )
where Aand B
are Jocobian matrixes of non- viscosity value fluxes E
and F
; and Av
andBv
are Jocobian matrixes of viscosity value fluxesEv ,Fv
. Q
is the micro variable of conservation variable in the time step. Using the sign of characteristic values, non- viscosity flux Jacobian matrix can be further decomposed into:
1
1
i i i i i i i ii
A A R R R R
A
( 2 - 1 7 ) where diagonal matrix i is composed by non-negative characteristic values of the
and matrixi i is composed by the non-positive characteristic values of .i
Using the implicit Euler scheme for space-time dicretization, Equation (3.1) becomes：
} ] S )
~F F~ [(
] S ) F~ F~ [(
] S ) E~ E~ [(
] S ) E~ E~ t [(
) Q Q ( V
1 n
j 2, i 1 v 1
n j 2, i 1 v
1 n
j 2, i 1 v 1
n j 2, i 1 v n
j , i 1 n
j , i j , i
( 2 - 1 8 )
where n is the time target.
Substituting the linearized flux from equation (3.22) and decomposed Jacobian matrix into eq.(3.24) and neglecting the second and higher order items, a diagonal line implicit expression can be obtained.
RHS S
F F S
F F
S E E S
E E
Q S
B B Q
S B B
Q S
B B Q
S B B
Q S
A A Q
S A A
Q S
A A Q
S A A Q
t I V
j v i j
v i
j v i
j v i
n j j i
v i j
j i v i
j j i
v i j
j i v i
j j i v i
j j i v i
j j i v i
j j i v i
j i j i
} ]
~) [(~ ]
~) {[( ~
} ]
~ ) [(~ ]
~ ) {[( ~
} ]
) [(
] ) [(
] ) [(
] ) [(
] ) [(
] ) [(
] ) [(
] ) {[(
2 , 1 2
, 1
2, , 1
2 1
, 2 , 1 1
, 2 , 1
1 , 2 , 1 ,
2 , 1
, , 2 1 ,
, 1 2 1
, , 1 2 1 , ,
2 1 ,
,
( 2 - 1 9 )
where I is the unit matrix. Using viscosity flux of the Jacobian matrix in the implicit scheme can increase the converging speed.
2.5 Computational Grid
The flow field contains several regions, with large gradients, which require careful resolution. Grid clustering is required in the vertical direction near the walls to resolve the boundary layer. The separated shear layer, which is located over the separation bubble, also needs special attention. Because of these requirements, a fine resolution should be maintained up to the middle of the channel in order to resolve the separated shear layer properly. The (x, y) plane was discretized into 10291 elements (see Figure 2-1). The number of grid points in the x direction was 251, in the y direction was 41. To reduce the computational cost, the governing equations were first integrated in time on a very coarse grid until the flow reached a statistically steady state. Then the flow field was interpolated onto a finer grid to get a restart flow field for the simulation with higher resolution. The following results were obtained:
(a) A streamwise resolution of at least 251 elements is required; higher resolution in this direction does not cause significant changes.
(b) Significant changes were found if less than 21 elements were used in the vertical direction. Higher resolutions than this improved the first-order statistics only slightly.
(c) As the resolutions in the streamwise and vertical directions improve, the shear stress along the upper wall increases.
Comparisons of streamlines, , of the mean flow fields are presented in Figure 2-2.
The boundary of the separated region corresponds to the ª = 0 contour. The size of the separation bubble is seen to change only slightly as the resolution increases.
F i g . 2 - 2 A g r i d a r r a n g e m e n t u s e s g r i d n u m b e r o f 2 5 1 × 4 1 .
Fig. 2-3 Comparison of the streamlines () with different resolutions. (a) 12541 elements; (b) 25141 elements; (c) 50141 elements.
2.6 Initial condition and boundary condition 2.6.1 Initial condition
The inflow speed is at Mach number 0.6, and the Reynolds number is 3.65518x105.
The pressure is in the static atmospheric condition and the density is 1.2019 kg/m3.
2.6.2 Boundary condition
Boundary value of the flow field is assumed to be the initial flow field exit condition and kept constant.
2.6.3 Boundary condition on top of base and below
y 0 u y 0
( 2 . 2 0 )
y 0 e
0 v
2.6.4 Inflow and outflow boundary condition:
The boundary condition is subsonic of sound and expressed as follows:
2 2 2 1
) 2(
1
1 p u v C
C u
( 2 . 2 1 ) where C1 and C2 are constants determined by the initial condition. The suffix 0 and e indicate exterior value of outflow and extrapolated internal value respectively.
2 1 0
2 1 2 1 0 0
) (
P) (
) (
e te t
e e e
P e P e
v v
M u
P P P P
( 3 . 8 )
3. Results
3.1 The flows over a wavy wall analyze 3.1.1. Vortex change
At boundary curved process, has three important vortex in this flow, see the Fig.
A.
F i g . A . T h e v o r t e x i n f l o w s o v e r a w a v y w a l l . Vortex A:
The wavy wall in the curving process, the curved angle increases gradually, because above in the space the right half is a partial high pressure region, causes the fluid backflow and produces the turbulence to the nearness inward flow boundary.
Vortex B:
This vortex and the vortex production way differs from. First in a (a) place production small vortex behind, after but only then forms the main vortex, but the small vortex finally to merge into the vortex to vanish gradually.
Vortex C:
This vortex production is from the underneath boundary production. The vortex size increases along with the boundary bending increases. When in the future will revise the wavy wall surface, the vortex could not immediately destroy, but will be along with the mainstream toward the downstream merge, and gradually dissipation.
3.1.2. Pressure change
Fig. B Pressure change of flows over a bilateral wavy wall.
Fig. B shows, has the very tremendous pressure gradient in front of wavy wall downhill a place B. Mainly in x and the y direction, has its different function respectively.
X direction:
When boundary unceasing curved, unceasing changes small in the flow field crook cross-sectional area, this speed of flow presents the increase in this the phenomenon.
Therefore in enters the neck of a bottle place, in front of the flow field can have x direction high-pressured gradient to use to accelerate the fluid.
curving section to have the high-pressured force gradient in the Y direction, but causes the fluid curve.
Flow over a downhill process.
Flow over a downhill process time, a (c) place spatial gradual expansion. But in order to enable the fluid to fill because of expands the space which are many, in the boundary has the obvious partial high pressure region, enables the fluid to reduce in x direction speed. Therefore in expansion time, here wall surface receives the large external pressure.
Flow over a above process.
Flow over a above process time, a (d) place the space changes small, causes the fluid upward to push. Now, the partial high pressure falls on (d) place, differs from with the boundary curving expansion.
Because the boundary migration speed to reduce along with the bending increases, therefore needs to produce the forward pressure gradient in the y direction, causes the fluid to enlarge gradually in the y direction speed, fluid which promotes, has not merged in the mainstream, but is used for to fill (e) place because of the space which contracts appears.
3.2 Wavy wall system
3.2.1 The subsonic and supersonic flow over a solid wavy wall
The main objective of this study is to investigate the characteristics of the wavy wall in subsonic flow field using Mach number 0.6 and compare results with numeric simulation made by Cherukat et al.[1]. For comparison purpose, this study uses the same flow field contour as in the two-dimensional flow field simulation of wavy wall.
Therefore the streamline direction is 230.2mm and cross streamline direction is 50.8 mm. The top of flow field board is a smooth wall surface and the bottom is a sinusoidal profile wall which has wavy oscillation amplitude(A) of 2.54mm and wave length (λ+) of 50.8 mm. Since the wavy wall and the flat surface are the critical areas, grids at both sides of the board are densified. Fig. 3-2 show the return flow region produced in the wavy wall downhill quarter for three different grid numbers at the same density and Mach number. The 251x41 grid number was seen to provide precise analysis for the return flow region. However the 125x41 grid number was found insufficient and unable to reveal the return flow region. Grid number of 501x41 has exhibited similar result to that of 251x41. Fig. 3-3 show instantaneous velocity distribution maps for time duration of 1.5 and 2.3msecs. At the same Mach number the wavy wall downhill return flow region increases with time. Fig. 3-4 shows instantaneous velocity chart for time T=2.3msec in the cross streamline direction. In order to demonstrate the
formation of return flow region in the wavy wall area, single wavelength to span stream-line direction to blink time uniform velocity line flow linear velocity, M=0.6 A=2.54mm T=2.3msec Comparing these results with reference in Fig. 3-5. We has a good agreement is observed. Then the perpendicular velocity line. Comparing these results with reference in Fig. 3-6. We has a good agreement is observed too. See the Fig. 3-7, 3-8 the streamwise velocity and spanwise velocity contour for M=0.6 at T=2.3msec, comparison with Hudson’s (1993) measurement and open circles is Cherukat et al., and the filled circles is present, Hudson at x/λ=0.0 (crest) and Cherukat et al. is open square, present is filled square. The turbulence intensities also agree reasonably wellwith Hudson’smeasurements.Theagreement is better for urms since vurms can be measured more accurately. In the recirculating zone, the mean velocity and the velocity fluctuations are very small and of the same order. This introduces errors in measurements and is one of the main reasons for differences in computed and measured turbulence intensities, see the Fig.3-9, 3-10. The Fig.3-11 shows pressure contour for M=0.6 at T=2.3msec. The Fig.3-12 shows Reynolds shear stress comparison with Cherukat et al. Comparison with Cherukat et al and present. This contour shows the low-pressure and highest speed occur at the peaks of wavy wall surface. The supersonic flow main objective of this study is to investigate the characteristics of the wavy wall in subsonic flow field using Mach number 1.5 and compare results with numeric simulation made by Cherukat et al.[3]. For comparison purpose, this study uses the same flow field contour as in the two-dimensional flow field simulation of wavy wall. Therefore the streamline direction is 230.2mm and cross streamline direction is 50.8mm. The top of flow field board is a smooth wall surface and the bottom is a sinusoidal profile wall which has wavy oscillation amplitude (A) of 2.54mm and wave length (λ+)of50.8 mmsee the Fig.3-1. Since the wavy wall and the flat surface are the critical areas, grids at both sides of the board are densified. See the Fig.13. Comparison with Kim(2001)’s experiment result about the point contours for supersonic single wavy wall time T=2.3msec in the density distribute contours. And Fig.3-14 the point streamline for much number M=1.5, A/λ=0.1. Comparison with Kim(2001) experiment result The flow over a wavy wall vector for much number M=1.5, A/λ=0.1 see the Fig.3-15. Mean velocity profiles at two streamwise locations (streamwise velocity and spanwise velocity) in the crest and thetrough see the Fig.3-16.
The Reynolds shear stress (see the Fig.3-17 through 19) are u'2 and v'2 can be expressed based on the time average method as follows:
2
2 ' 0T(u U) dt
u T
( 3 . 1 )2
2 ' 0T(v V) dt
v T
( 3 . 2 )Where T is the statistical time; for this statistical time, U and V are the time average speed in the streamline and cross streamline directions respectively. In this study, the simulation uses the following ensemble average relations:
2
2
' 1
( )
N n n
u U
u N
( 3 . 3 )
2
2
' 1
( )
N n n
v V
v N
( 3 . 4 )Where N is folding sequence, the statistical time T N t.
3.2.2 The to compare of subsonic and supersonic flow over a solid wavy wall The present paper research main purpose is observes the wavy wall in subsonic and supersonic of sound flow field characteristic using Mach M=0.6 and M=1.5, and with Tsai[1] and Cherukat et al.[3]. In the two-dimensional wavy wall in the flow field value simulation, this article uses with Tsai[1] and Cherukat et al.[3]. The same flow field contour is advantageous for the comparison, therefore the streamwise is 203.2 (mm), spenwise 50.8 (mm). In flow field boarding up for the smooth wall surface, gets down the board is the Sinusoidal wavy wall, but under board wavy wall oscillation amplitude A=2.54(mm), the wave length λ=50.8(mm) the computation territory length picks with the comparison article same 4 wave lengths, the grid number is 251x41, meaning is grid of number the streamline direction is 251, the cross streamline direction grid number is 41. In the grid selection, by the 125x41 test discovery, the grid number too few causes to be unable first to simulate in the flow field the wavy wall downhill return flow region, then by the high grid number test, the result and grid of number the 251x41 and does not have the too big difference, in order to save the time, therefore uses grid of number the 251x41 to carry on the simulation. Fig.3-20,21 show instantaneous velocity distribution maps for time duration of 1.5 and 2.3 msecs. At the same Mach number (M=0.6) the wavy wall downhill return flow region increases with time. Fig.3-22 shows Mach number (M=1.5) instantaneous velocity distribution maps for time duration of 2.3 msecs. At the same time duration of 2.3 msecs the wavy wall downhill vortex region increases of subsonic. Fig.3-23, 3-24 shows subsonic and supersonic instantaneous velocity chart for time T=2.3msec in the cross streamline direction. Fig.3-25 and Fig.3-26 show the comparison of streamline patterns which clearly show the streamline direction, the formation of vortex and the relation between magnitude of subsonic and supersonic flow. Apparently the formation of vortex occurs earlier and the size increases of supersonic flow. Pressure curve from the studies for
much number M=0.6 and 1.5 time T=2.3msec. Comparison with Tsai[1] M=0.6 and Cherukat[3] and present M=1.5. In Fig.3-27 is from the studies of Tsai[1]. Lowest pressure and highest speed occur at the peaks of wavy wall surface. Fig.3-28 shows the agreement is better for urms since vurms can be measured more accurately. In the recirculating zone, the mean velocity and the velocity fluctuations are very small and of the same order. This introduces errors in measurements and is one of the main reasons for differences in computed and measured turbulence intensities.
3.2.3 By physical space comparison
The same flow field contour is advantageous for the comparison, therefore the streamline direction is 203.2mm, cross streamline direction 50.8mm. In flow field boarding up for the smooth wall surface, gets down the board is the Sinusoidal wavy wall, but under board wavy wall oscillation amplitude A=5.08mm and A=2.54mm, the wave length λ=50.8mm. By physical space see the Fig3-29, Fig3-30. This research aims at the amplitude of two kinds see the Fig3-31, Fig3-32. And single wavy wall, bilateral wavy wall and flat, three different computational grid see the Fig.3-31 through Fig3-34. The computation territory length picks with the comparison article same 4 wave lengths, the grid number is 251x41, meaning is grid of number the streamline direction is 251, the cross streamline direction grid number is 41. In the grid selection, by the 125x41 test discovery, the grid number too few causes to be unable first to simulate in the flow field the wavy wall downhill return flow region, then by the high grid number test, the result and grid of number the 251x41 and does not have the too big difference, in order to save the time, therefore uses grid of number the 251x41 to carry on the simulation. The grid arrangement uses grid number of 251×41. The turbulent flow occurs mostly in the flow field nearby the wavy wall surface.
Fig.3-35 and Fig.3-36 show instantaneous velocity distribution graphs at the same Mach number of 0.6 and the same oscillation amplitude of 5.08mm for time of 0.5 and 1.5msec. Comparing with the results for the oscillation amplitude of 2.54 mm, obviously the fluid detachment point in the wavy wall downhill area has occurred earlier. The time for formation of the return flow region was found to be in reverse relation to the magnitude of oscillation amplitude, i.e., the higher magnitude of oscillation amplitude, the quicker to form a bigger return flow region in the wavy wall downhill area. Fig.3-37 is the instantaneous vertical velocity graph when the return flow region forms. Fig.3-38 through Fig.3-40 show the comparison of streamline patterns which clearly show the streamline direction, the formation of return flow region and the relation between magnitude of oscillation amplitude size, and time duration. Fig.3-41 shows the speed contours from Cherukat et al.[3] for varying x/λ
this study which also reveals the formation of return flow region in the wavy wall downhill area.
Speed contour in Fig.3-43 is under the same Mach number M=0.6 and the time duration but the oscillation amplitude doubles (A=5.08mm). See the Fig.3-44 through Fig.3-46, shows three different density contour. The by way of density contour our to be clear about in the flow field the density change different has a difference along with the boundary, especially of bilateral wavy wall. The supersonic single wavy wall, bilateral wavy wall and flat boundary, three different instantaneous uniform velocity graphs of for Fig.3-14 in the streamline direction are shown in Fig.3-47 through Fig.3-49. And Fig.3-50 through Fig.3-52 is the instant uniform velocity graph for Fig.3-13 in the difference along with the boundary of cross streamline direction.
Fig.3-53 through Fig.3-55 show the comparison of streamline patterns which clearly show the streamline direction, the formation of vortex and the relation between magnitude of subsonic and supersonic flow. Apparently the formation of vortex occurs earlier and the size increases of supersonic flow. See the Fig.3-56, pressure curve from the studies for much number M=1.5 time T=2.3msec. Comparison with single and bilateral wavy channel wall and flat boundary. The lowest pressure and highest speed occur at the peaks of wavy wall surface and the boundary range of variation is bigger regarding the pressure influence is also bigger.
4. Conclusions
This study successfully simulates the flow field characteristics of the two-dimensional profile wall using the WENO and LU-SSOR schemes for space-time discretization and a grid of 251x41 with both sides densified. The local turbulence, when scaled with this velocity, is less influenced. The simulation agrees with the picture presented by Cherukat[3] and Hudson et al. (1995). Turbulence near a wavy surface is very different from that near a flat wall, in that production is not associated with the floworiented vortices described by Brooke and Hanratty (1993). Rather, it is associated with a shear layer which is formed by flow separation and extends over the full wavelength. Flow near the wavy surface shows very large spanwise fluctuations in the region between the trough and the next crest where reattachment occurs. Also numerical method used in this study has been validated and the results for supersonic;
subsonic and three different physical space comparison and supersonic mixing flow were found very similar to those simulated by Tsai[1], Kim[2] and Cherukat et al[3].
Observations were also made on increase of two times in oscillation amplitude. Further studies on comparing the two methods are recommended. As oscillation amplitude of the wavy wall increase for two times, the fluid separation and the contact phenomenon can occur earlier.
References
[1] Yeong-Pei Tsai and Chung-Chyi Chou 2007. “DirectNumericalSimulation ofa Fully DevelopedTurbulentFlow overaWavy Wall”Chung Hua Journal of Science and Engineering, Vol.5,No.2,pp.9-14.
[2] S. K. Kim 2001. “An Experimental Study of Flow in a Wavy Channel by PIV”The 6THAslan Symposium on Visualization.
[3] P. Cherukat, Y. Na, T.J. Hanratty, and J.B. McLaughlin.1998. “Direct Numerical Simulation of a Fully Developed Turbulent Flow over a Wavy Wall”Theoretical and Computational Fluid Dynamics 11,109–134.
[4] P.Zilker, J.Hanratty. 1979. “Influence of the amplitude of a solid wavy wall on a turbulence flow. Part 2”separated flows, J. Fluid Mech 90, 257-271.
[5] P.Zilker, W.Cook, and J. Hanratty. 1977. “Influence of the amplitude of a solid wavy wall on a turbulent flow. Part 1”Non-separated flows 82 , 29-51.
[6] V. C. Patel, W. Rodi and G. Scheuerer. 1985. “Turbulence models for near-wall and low Reynolds number flows: a review”AIAA J. 23. 1308-1319.
[7] V. C. Patel and H. C. Chen. 1988. “Near-wall turbulence models for complex flows including separation”AIAA J. 26. 641-648.
[8] V. C. Patel, J. T. Chon and J. Y. Yoon. 1991 “Turbulent flow in a channel with a wavy wall”, J. of Fluids Eng. 113, 579-586.
[9] A. D. Ferreira and A. M. G. Lopes etc. 1995. “Experimental and numerical simulation of flow around two-dimensional hills”, J. of Wind Eng. and Industrial Aerodynamics. (54/55), 173-181.
[10] E. D. Montalbano and M. J. McCready. 1998. “Laminar channel flow over long and moderate waves”, J. Comp. Phys. 6. 1-23.
[11] C. Airiau and A. Giovannini. 1999. Vorticity evolution on a separated wavy wall flow. Third International Workshop on Vortex Flows and Related Numerical Methods. 7. 1-11.
[13]Daniel P, Zilker, Gerald W. Cook and Thomas J. Hanrattya1 2006. “Influenceofthe amplitude of a solid wavy wall on a turbulent flow. Part 1. Non-separated flows” Cambridge University Press 12Apr.
[12] N. A. Malamataris1 and V. Bontozoglou. 1999 “Computer Aided Analysis of ViscousFilm Flow along an Inclined Wavy Wall” JournalofComputationalPhysics 154, 372–392
[13] Bendiks Jan Boersma 2000. “Particledistributionsin theflow overawavy wall” Center for Turbulence Research Proceedings of the Summer Program 109.
[14] S. Nakagawa, Y. Na, T.J. Hanratty. 2003. “Influence of a wavy boundary on turbulenceI.”Highly rough surfaceExperiments in Fluids. 35. 422–436
[15] Harten, A.,(1989), “ENO Schemes with Subcell Resolution. Journal of Computational Physics”, 83. 148-184.
[16] Liu, et al.(1994), “Weighted Essentially Non-oscillatory Schemes. Journal of Computational Physics”115. 200-212.
[17] Jiang, G.S., and Shu, C. W. 1996. “Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics”, 126. 202-228.
Figure
2 5 141 element in (x,y)-plane
Fig. 3-1. Element distribution in the present study (10291 elements).
Fig. 3-3 Streamwise velocity contour for M=0.6 at. (A) T=1.5 msec (B)T=2.3 msec.
Fig. 3-4 Spanwise velocity contour for M=0.6 at T=2.3msec.
Fig. 3-5 Comparison with streamwise velocity contour for M=0.6 at T=2.3 msec.
(A)Present; (B)Cherukat et al.
(B)
(A)
(B) (A)
Fig. 3-7 Comparison with Hudson’s (1993) measurement;─, present;●, Hudson at x/λ=0.0 (crest); ----, present;■. (a) Streamwise velocity; (b) Spanwise velocity.
(A)
(B)
Fig. 3-8 Comparison with Cherukat et al.;○, present;●, Cherukat et al. at x/λ=0.0 (crest);□present;■. (a) Streamwise velocity; (b) Spanwise velocity.
(B) (A)
Fig. 3-9 Comparison with Hudson’s (1993) measurement;─, present;●. (a) x/λ= 0:0 (crest); (b) x/λ= 0:5 (trough).
(A)
(B)
Fig.3-10 Comparison with Cherukat et al.;○, present;●. (a) x/λ= 0:0 (crest);(b)x/λ=
0:5 (trough).
(A)
(B)
Fig. 3-11 Comparison with pressure contour for M=0.6 at T=2.3 msec. Cherukat et al.;
─, present;---.
Fig. 3-12 Reynolds shear stress (normalized by U2b ) x/λ=0.0.
Fig. 3-13 Point contours for supersonic single wavy wall time T=2.3msec in the density distribute contours, (a) Present; (b) Kim(2001). experiment result.
Fig.3-14 The point streamline for much number M=1.5, A/λ=0.1; (a) Present；(b) (B)
(A)
(A)
(B)
Fig.3-15 The point vector for much number M=1.5, A/λ=0.1;
(a) Present；(b) Kim(2001). experiment result.
(A)
(B)
Fig.3-16 Mean velocity profiles at two streamwise locations (the crest and the trough).
(B) (A)
Fig.3-17 Reynolds shear stress (normalized by U2b)x/λ=0.0 (A)step 1;(B) step 2;
(C)step 3; (D)step4.
(A) (B)
(C) (D)
Fig.3-18 Reynolds shear stress (normalized by Vb2)x/λ=0.0 (A)step 1;(B) step 2;
(C)step 3; (D)step4.
(A) (B)
(C) (D)
Fig. 3-19 Reynolds shear stress (normalized by Vb2)x/λ=0.0(A)step 1;(B) step 2;
(C)step 3; (D)step4.
(A) (B)
(C) (D)
Fig. 3-20 The point contours for subsonic time T=1.5msec in the mean streamwise and normal velocity.
Fig. 3-21 The point contours for subsonic time T=2.3msec in the mean streamwise and normal velocity.
Fig.3-22 The point contours for supersonic time T=2.3msec in the mean streamwise and normal velocity.
Fig.3-24 Shows supersonic instantaneous velocity chart for time T = 2 . 3 m s e c i n the cross streamwise direction.
Fig.3-25 The point streamline for much number M=0.6, time T=2.3msec, amplitude A=2.54mm.
Fig.3-26 The point streamline for much number M=1.5, time T=2.3msec, amplitude A=2.54mm.
Fig.3-27 Pressure curve from the studies for much number M=0.6 and 1.5 time T=2.3msec. Comparison with Tsai[1] M=0.6;---, Cherukat[3] M=0.6;─,
Fig.3-28 Comparison the turbulence intensities also agree reasonably well with Subsonic;─, Supersonic;－－.(a) crest (b) trough.
(B) (A)
Fig.3-29 Computational domain of flow over a single wavy wall.
Fig.3-30 Computational domain of flow over a bilateral wavy wall.
Fig.3-32 The computational grid of single wavy wall oscillation amplitude A=5.08(mm), the wave length λ=50.8(mm).
Fig3-33 The computational grid of bilateral wavy wall oscillation amplitude A=5.08(mm), the wave length λ=50.8(mm).
Fig3-34 The computational grid of flat oscillation amplitude A=0(mm), the wave length λ=50.8(mm).
X/CL
Y/CL
1 1.25 1.5 1.75
-0.1 0 0.1 0.2 0.3 0.4
Fig.3-35 The point one single wavelength was used in the simulation contours for amplitude A=5.08, time T=0.5msec in the streamwise and normal velocity.
X/CL
Y/CL
1 1.25 1.5 1.75
-0.1 0 0.1 0.2 0.3 0.4
Fig.3-36 The point one single wavelength was used in the simulation contours for amplitude A=5.08, time T=1.5msec in the streamwise and normal velocity.
X/CL
Y/CL
1 1.25 1.5 1.75
-0.1 0 0.1 0.2 0.3 0.4
Fig.3-37 The point one single wavelength was used in the simulation contours for amplitude A=5.08, time T=1.5msec in the spanwise and normal velocity.
X/CL
Y/CL
1 1.25 1.5 1.75
0 0.1 0.2 0.3 0.4
Fig.3-38 The point streamline for much number M=0.6, time T=1.5msec, amplitude A=2.54mm.
