二維均勻流及剪流流經正方阻塊之數值分析

國

立

交

通

大

學

機械工程學系

碩

士

論

文

二維均勻流及剪流流經正方阻塊之數值分析

Numerical Analysis of Two-Dimensional Uniform and

Shear Flows past a Square Cylinder

研 究 生：柯 伯 翰

指導教授：崔 燕 勇 教授

Numerical Analysis of Two-Dimensional Uniform and

Shear Flows past a Square Cylinder

二維均勻流及剪流流經正方阻塊之數值分析

Student: Ko, Po-Han Advisor: Dr. Tsui, Yeng-Yung 學生: 柯伯翰 指導教授: 崔燕勇 博士

國 立 交 通 大 學 機械工程學系

碩士論文

A Thesis

submitted to Department of Mechanical Engineering College of Engineering

National Chiao Tung University in partial fulfillment of the requirements

for the Degree of Master of Science in

Mechanical Engineering July 2006

Acknowledgments

I would like to thank Prof. Tsui, Yeng-Yung for his recommendations, patience, and so much guidance.

I would like to thank the individuals in the CFD Laboratory.

Thanks to one of the senior students Hu, Yu-Chang for his helping and teaching when I was really confused.

Thanks to Clark Lin and Yang, Ching-Shiang for their friendships. I do improve myself through all discussions with them.

Thanks to Flyduck Cheng, Eatrol Lu, and Just Jhu for also their friendships. Their so many suggestions helped me to move further.

I would like to thank my parents, papa Kenns, mom Zung, and my brother Kurt for their support and encouragement.

Numerical Analysis of Two-Dimensional Uniform and

Shear Flows past a Square Cylinder

Student: Ko, Po-Han Advisor: Dr. Tsui, Yeng-Yung

Department of Mechanical Engineering National Chiao Tung University

Abstract

The two-dimensional unsteady flow field over a square cylinder has been numerically predicted by finite volume method (FVM) with a non-iterative method PISO, standing for Pressure-Implicit with Splitting of Operators, solving implicitly discretised time-dependent fluid equations. The uniform and the shear free stream were considered. The shear flow was studied with two different velocity ratio 3:1 and 3:2 as well as two different distances 6 and 10 in unit length between the square cylinder and the inlet boundary of the computational domain while Reynolds number was in the range from 50 to 500. The higher velocity ratio of free stream brings the lower Strouhal frequency. The negative vortices dissipate more rapidly then the positive ones do even thought the negative background dominates. The direction of lifting force is pointing to the faster free stream side of the shear flow cases.

二維均勻流及剪流流經正方阻塊之數值分析

研究生 柯伯翰 指導教授 崔燕勇 博士 國立交通大學機械工程學系 摘要 本論文主要探討二維非穩態的流場中，流體流經一方形阻塊，數值方 法使用有限體積法(FVM)來離散，使用隱式非疊代的 PISO 法，即 Pressure-Implicit with Splitting of Operators，來求解暫態的流體方程 式。當中分別探討了均勻流及剪力流兩種不同的入口條件；於剪力流 的條件下分別以 3：1 及 3：2 兩種速度比，以及從入口邊界到方形阻 塊之間以 6 和 10 不同的單位長度加以比較， 探討的雷諾數範圍從 50 到 500。結果發現到，在入口速度較高的速度比之下，會產生較低的 Strouhal 頻率；而在整個流場當中都是負的渦度所主導下卻發現，負 的渦度耗散還比正的渦度耗散還要迅速；在剪力流的案例中，舉升力 是指向速度較快的那一方。

Contents

Abstract………i List of Figures……….v Nomenclature………...viii

1 Introduction

………...………1 1.1 Bluff Body……… ..……..2 1.2 Imposed Frequency……….…..5 1.3 Shear Flow………6

1.4 Problem under Consideration………...………….7

2 Mathematical Model

……….9 2.1 Physical Model……….……….……9 2.2 Governing Equation………..9 2.3 Dimensionless Group……….……….10 2.4 Boundary Conditions……….………..………10

3 Discretization

………..……..11 3.1 Introduction………...11

3.2 Finite Volume Method...…...11

4 Numerical Method for Unsteady Flow Calculation

…………..15

4.1 Introduction………...15

4.2 Details of the PISO Algorithm...…...16

4.3 Implementation of Boundary Conditions………...……….20

4.4 Solution Procedure of PISO………21

5 Validation of Numerical Code

……..………...22

5.1 Grid Test………..22

5.2 Data Comparison……….22

5.3 Time Step Size………24

6 Results and Discussion

………...25

6.1 Relation of St and Re with Shear Free Stream………25

6.3 Lift Coefficient………27

6.4 Drag Coefficient………..29

6.5 Instantaneous Vorticity Contours………30

6.5.1 Uniform Flow……….30

6.5.2 Shear Flow………..31

6.6 Instantaneous Streamlines………...32

6.6.1 Uniform Flow……….32

6.6.2 Shear Flow………..33

6.7 Time Evolution of Flow Fields………...33

6.8 Elongation of Computational Domain………34

7 Conclusions

………35

List of Figures

Figure 2-1. Illustration of the physical model……….…..39 Figure 3-1. Illustration of the primary cell P and the neighbor cell nb with a

face f in between………...39 Figure 3-2. Illustration of the primary cell P and the neighbor cell nb with a

face f in between………...39 Figure 4-1. Block diagram of solution procedure of PISO algorithm……...40 Figure 5-1. The computational mesh and geometry………..40 Figure 5-2. Variation of Strouhal number with Reynolds number (Uniform

free stream cases of cell numbers 198x97 and 232x123)……….…………..41 Figure 5-3. Variation of Strouhal number with Reynolds number (Uniform

free stream)………...41 Figure 5-4. Variation of transverse velocity [at point (x, y)= (1.23, 0)]

spectra with Reynolds number (Uniform free stream)……….42 Figure 5-5. Variation of phase space of velocity components at (x=1.23, y=0)

with Reynolds number (Uniform free stream)………..44 Figure 5-6. Variation of mean drag and lift coefficients with Reynolds

number (Uniform free stream)………..45 Figure 5-7. Variation of Strouhal number with Reynolds number (Uniform

free stream cases of time step sizes 0.01 and 0.005)…………46 Figure 6-1. Variation of Strouhal number with Reynolds number (Uniform

and shear free streams)……… 46 Figure 6-2. Variation of phase space of velocity components at (x=1.23, y=0)

with Reynolds number (Velocity ratio 3:1 and L1=6) ………..………47 Figure 6-3. Variation of phase space of velocity components at (x=1.23, y=0)

with Reynolds number (Velocity ratio 3:2 and L1=6) ………...…49 Figure 6-4. Variation of phase space of velocity components at (x=1.23, y=0)

with Reynolds number (Velocity ratio 3:1 and L1=10)………..51 Figure 6-5. Variation of phase space of velocity components at (x=1.23, y=0)

with Reynolds number (Velocity ratio 3:2 and L1=10)………..53

Figure 6-6. Variation of transverse velocity [at point (x, y)= (1.23, 0)] spectra with Reynolds number (Velocity ratio 3:1 and L1=6)………....55 Figure 6-7. Variation of transverse velocity [at point (x, y)= (1.23, 0)]

spectra with Reynolds number (Velocity ratio 3:2 and L1=6)………....57 Figure 6-8. Variation of transverse velocity [at point (x, y)= (1.23, 0)]

spectra with Reynolds number (Velocity ratio 3:1 and L1=10)………..59 Figure 6-9. Variation of transverse velocity [at point (x, y)= (1.23, 0)]

spectra with Reynolds number (Velocity ratio 3:2 and L1=10)………..61 Figure 6-10. Variation of mean lift coefficient with Reynolds number (Shear free streams)……….63 Figure 6-11. Variation of mean drag coefficient with Reynolds number

(Uniform and shear free streams)……….63 Figure 6-12. Variation of instantaneous vorticity contours [Broken line (ωmin,

max

ω , Δω)≡(-20.22, 0, 0.8088); Solid line (ωmin, ωmax, Δω)≡(0,

20.22, 0.8088)] with Reynolds number (Uniform free stream)………...64 Figure 6-13. Variation of instantaneous vorticity contours [Broken line (ωmin,

max

ω , Δω )≡(-20.22, -0.608, 0.862); Solid line (ωmin, ωmax, ω

Δ )≡(0.218, 14.26, 0.862)] with Reynolds number (Velocity ratio 3:1 and L1 =6) ………...66 Figure 6-14. Variation of instantaneous vorticity contours [Broken line (ωmin,

max

ω , Δω )≡(-20.22, -0.608, 0.862); Solid line (ωmin, ωmax, ω

Δ )≡(0.218, 14.26, 0.862)] with Reynolds number (Velocity ratio 3:2 and L1=6)……… ……...68 Figure 6-15. Variation of instantaneous vorticity contours [Broken line (ωmin,

max

ω , Δω )≡(-20.22, -0.608, 0.862); Solid line (ω_min, ω_max,

ω

Δ )≡(0.218, 14.26, 0.862)] with Reynolds number (Velocity ratio 3:1 and L1=10)……….70 Figure 6-16. Variation of instantaneous vorticity contours [Broken line (ω_min,

max

ω , Δω )≡(-20.22, -0.608, 0.862); Solid line (ω_min, ω_max,

ω

Δ )≡(0.218, 14.26, 0.862)] with Reynolds number (Velocity ratio 3:2 and L1=10)…… ………...72 Figure 6-17. Variation of instantaneous streamlines with Reynolds number (Uniform free stream)………...74

Figure 6-18. Variation of instantaneous streamlines with Reynolds number (Velocity ratio 3:1 and L1=6)………..………….76 Figure 6-19. Variation of instantaneous streamlines with Reynolds number

(Velocity ratio 3:2 and L1=6)………...78 Figure 6-20. Variation of instantaneous streamlines with Reynolds number

(Velocity ratio 3:1 and L1=10)………..….. 80 Figure 6-21. Variation of instantaneous streamlines with Reynolds number

(Velocity ratio 3:2 and L1=10)……….82 Figure 6-22. Time evolution of instantaneous flow within a period T

(Uniform free stream; Re=100)……….84 Figure 6-23. Time evolution of instantaneous flow within a period T

(Uniform free stream; Re=200)……….86 Figure 6-24. Time evolution of instantaneous flow within a period T

(Uniform free stream; Re=300)……….88 Figure 6-25. Time evolution of instantaneous flow within a period T

(Uniform free stream; Re=400)……….90 Figure 6-26. Time evolution of instantaneous flow within a period T (Shear

free stream; Velocity ratio 3:1; L1=6; Re=100)………92 Figure 6-27. Time evolution of instantaneous flow within a period T (Shear

free stream; Velocity ratio 3:1; L1=6; Re=200)………94 Figure 6-28. Time evolution of instantaneous flow within a period T (Shear

free stream; Velocity ratio 3:1; L1=6; Re=300)………96 Figure 6-29. Time evolution of instantaneous flow within a period T (Shear

free stream; Velocity ratio 3:1; L1=6; Re=400)………98 Figure 6-30. Variation of Strouhal number with Reynolds number (Uniform

free stream cases of computational domain 10x24 and 10x32)………..100 Figure 6-31. Variation of mean drag coefficient with Reynolds number

(Uniform free stream cases of computational domain 10x24 and 10x32)………..100

Nomenclature

A Coefficients of algebraic equation

D

C Drag Coefficient

L

C Lift Coefficient

D

C Mean Drag Coefficient

L

C Mean Lift Coefficient

Cr Courant number u_cΔt Δx

p

f Weighting factor

s

f Strouhal frequency

H Height of computational domain

L Side length of square cylinder

•

f

m Mass flow rate cross face

P Pressure Re Reynolds number VLν S Source term Sv Surface vector St Strouhal number f_sL V t Time c u Convective velocity

u,v Cartesian components of velocities

V Average value of inlet velocity

x,y Components of Cartesian coordinate

β Blockage ratioL H

Pnb

δv Vector pointing from center of primary cell to the center of neighboring one

ν Kinetic viscosity φ Property Γ Diffusivity ∀ Δ Volume ρ Fluid density μ Viscosity coefficient

SIMPLE Semi-Implicit Method for Pressure-Linked Equations PISO Pressure-Implicit with Splitting of Operators

∇ Gradient

n New time level

Chapter 1

Introduction

The new technologies are dominating the world and people are so depending on those technologies to have the lives easier and enjoyable. Once the wind breezes, poets and artists are about to remember it in any kinds of beautiful ways. It happens but it was the only thing happened. The wind comes from nature and nature is not just beautiful it is amazingly helpful. Although sails started to be utilized to move boats by the ancient Egyptians thousands years ago, the applications were broadly presented since fifteenth century. Until 1903, Wright brothers introduced a craft bringing people to the sky. The interactions between fluid and structures were quite doing the jobs. Some interesting phenomena were well described by a German scientist Ludwig Prandtl with the discovery of boundary layer occurred at the fluid-structure interface, 1904. After the giant breakthrough, researches and studies were following the trend so fast to move even further. The airfoil by far keeps playing the most important role for sure but the fluid-structure interaction is not just about the airfoil which is a big topic but still partial. The structure suffered from the blowing fluid does not have to be streamlined in shape. It may be of any kind of funny shapes being seen anywhere. A special vibrating motion would happen in the wake of the unstreamlined bodies, i.e. bluff bodies, and somehow this motion would be vibrating the bodies as well. Break step march is also the similar way by soldiers to induce an imposed frequency which happens to create the resonance on the bridge, and then to do the damage. Due to the existence of wake oscillation, engineers have to really pay attention on the fluid-structure interaction of

unstreamlined bodies instead of the airfoil industries only. Modern skyscrapers are also part of the studies. Because of the boundary layer effect, the wind at the higher part of the buildings is undoubtedly blowing so fast, and so creating the stronger oscillation and some unexpected harmful instabilities to the structure. It does not have to wait for earthquakes to come; they would be easily destroyed just by the beautiful breeze. Bluff body studies are so much mentioned with the safety-maintenances which are so about our living lives.

1.1 Bluff Body

Bluff, unstreamlined bodies have been working on so many places. People had been unconsciously ignoring the most interesting part for thousands years until the wire was vibrating so hard with the Aeolian tone coming out. Strouhal at 1878 started the experiment of the interaction between bluff bodies and the blowing wind. Even the wind went smoothly and uniformly, the existence of fluctuating forces which was exciting Rayleigh so much at 1896 was attributed to the vortex shedding from the leeward surface of the circular cylinder body into the wake with the alternating motion periodically to the downstream [1] once Reynolds number went beyond a certain number otherwise a pair of symmetric vortices stayed right behind the body steadily [2]. Vortex shedding phenomena were caused by the separated boundary layers from the upper and lower surfaces of the bluff body [3]. Lots of scientists and engineers were getting the researches begun after the fascination had been widely spread. Strouhal number was, hence, theoretically defined by Roshko [4], Bearman [5], and Griffin [6] for scaling the vortex street. Due to the computer had not been sufficiently well-performance yet, people were so much hitting the road from the experimental parts. The fluctuating vortex street was experimentally found to be

complex and varied along the increase of Reynolds number from 40 to 10,000 while the phenomena changed from stable to transition, and then to the irregular movement [4]. Furthermore, Hama [7] visually found out the three-dimensionality in the wake of a circular cylinder. The decreasing span-wise frequency was observed with the increasing Reynolds number great than 150. Despite the wake instability interested many researchers devoting their sophisticated experimental skills, the Strouhal-Reynolds number relation from different labs somehow deviated from each other by about 20 percent still. Until 1988, Williamson and Roshko proposed that the three-dimensional effect introduced by the so-called oblique shedding mode happened to be the critical reason of those disagreements for years. Although they all brought us the different experimental data, the contributions were still made. The discontinuities of Strouhal frequency along the increase of Reynolds number were reported by Trittin [8]. The causes of these discontinuous variations were explained to be attributed to the unskilled flow control which did not maintain the flow uniformly [9]. The forced cylinder vibration was told to be one of the possible reasons by observing the existence of quasi-periodic state in the laminar wake and the transformation to chaotic state [10]. The cylinder vibration triggering the discontinuities was still supported [11]. The quasi-periodic state was actually confirmed but the flow non-uniformity and the forced vibration were on the contrary examined not to be exact the main reasons causing the discontinuities [12]. Williamson and Roshko [12] indicated an intrinsic three dimensional wake instability which potentially caused the discontinuities of Strouhal frequency as the Reynolds number went beyond 178 and 260 even under a controlled parallel shedding, i.e. two-dimensional flow, of circular cylinder. This critical Reynolds number 178 was said to be the onset point of the

transition of wake from two-dimensional to three-dimensional. Although the wire was the research topic, the bluff body does not have to be always a circular cylinder. The shape of interest was thought to be sharp-corners or to be something not that smooth. The mechanism responsible for vortex shedding of rectangular cylinder was sort of different from the one of circular cylinder. The development of Vortex Street was triggered by the separation from the leading edge of a square cylinder instead of from the leeward side of circular one while Reynolds number was sufficiently high. The data of square and of other rectangular cylinders reported in detail and systematical were achieved by Okajima [13] although the width-to-height ratio of a rectangular cylinder had already been considered to be the contributing factor by Nakaguchi et al. since 1968 [14]. The vortex street had been a popular topic with Reynolds number staying high, whereas the finding of the onset critical Reynolds number for the vortices turning from symmetric to asymmetric and from steady to unsteady was the topic as well. Works had never been too easy even the numerical analysis came to be the one doing the prediction of the phenomena of interests. The onset Reynolds number turned out to be more than just one certain number due to the blockage effect. A bigger blockage ratio brought a higher value of onset Reynolds number [15, 16] as well as the higher Strouhal frequency and drag coefficient [17]. Two critical values of Reynolds number are mentioned so far. One is the onset point for flow turning from steady to unsteady and the other one is the onset point of three-dimensionality. The 2D simulation would over predict the drag and lift coefficient once the flow goes beyond the critical Reynolds number of three-dimensionality [18]. Many researchers were quite focusing on the second onset point with the interest of three-dimensionality and they defined

Once the flow features turn into three-dimensional, span-wise instabilities appear with different wavelengths which classify the particular features with modes A and B [12]. Mode A was characterized by longer span-wise wavelength by about two to four diameters of cylinder with larger vortex dislocations while second onset Reynolds number was achieved. Mode B happened with shorter wavelength by about one diameter instead while Reynolds number was about 230 whose value was higher than the one of second onset point. Although these measurements or predictions of mode classification were done with circular cylinder flow, those of rectangular or sharp-corners cylinders possessed the similarities as well [19]. Even a new Mode S with span-wise wavelength of 2.8 diameter between those of Mode A and B was observed in the wake of square cylinder while Reynolds number equaled 200 [20]. More complex phenomena could be found once people were doing the studies of the Karman vortex street in the wake of the rectangular or square cylinder. The nonlinear dynamic system of the wake was truly interesting. The frequency spectrum was studied from first onset Reynolds number to what induced the chaos. This transition process was somehow indicating the existence of three-dimensionality.

1.2 Imposed Frequency

Once the unfixed body is with the incident flow coming with an imposed frequency, the body as well as the wake may follow the imposed frequency to oscillate, the resonance is usually termed lock-on which could be also induced by self-vibration bluff body with certain range of the imposed frequency either from the incident flow or bluff body itself [21, 22, 23, 24, 25]. The self-vibrating bluff body could move in various ways for researchers to observe the lock-on effect. Armstrong et al. [26] and Barbi et al. [27] examined the vortex resonance

on the body oscillating in-line with the coming uniform flow, and the results were quite closed to the earlier experiments done by Griffin and Ramberg [28]. These experiments were also showing that the certain and appropriate imposed frequency being dominating the lock-on effect had to be somewhere in the range closing to twice of the natural or Strouhal frequency of the wake. The span width of the range was proportional to the magnitude or amplitude of the imposed oscillation. The result of cross-flow oscillation also described the similar shape of the lock-on regime measured by Koopmann [22]; the difference was that the cross-flow oscillation was locking the body and the wake in the lifting direction instead of in the dragging one. Moreover, the oscillation frequency in the lifting direction was naturally one half of the one in the dragging direction of the body [21], the lock-on regime of the cross-flow oscillation was, hence, locating near the Strouhal frequency rather than twice of it. Vortex resonance or lock-on could be also observed with small rotational oscillation of bluff body. Although it was rarely studied, Tokumaru and Dimotakis [29] and Filler et al. [30] analogically found the phenomena of amplitude-dependent range of lock-on frequencies.

1.3 Shear Flow

Researches of uniform flow have been widely studied and well-developed experimentally and numerically although there may be some unrevealed phenomena possibly remained. Again back to the study of bridges, the structures so attached to the ground are obviously immersed in the boundary layer which possesses the shear flow inducing something even more complicated and practical. The most investigations were done with an uniform shear flow of which the velocity profile varied linearly along the transverse direction. The

identical assumption of the numerical study was even surprisingly leading two opposite results and the mean lift and the mean drag forces on the cylinder were quite an issue never allowed to happen by scientists and engineers. The mean lift force was found to be likely to point to either the slower free-stream velocity site [31, 32, 33, 34] or the faster site [35, 36]. Despite the diverse shear rates applied by those authors might be the reason [37], the conflict kept on being discussed because the discrepancies were not just simply coming from the forces exerted to the cylinder. The existence of vortex shedding was even argued [33, 38]. The non-uniform flow is by far more complicated but unfortunately it is practical and easily met. The unnecessary vibration makes the noise, and devastatingly brings the failure of bridges or even the very tall skyscrapers in the future.

1.4 Problem under Consideration

The present study is focusing on the flow over a square cylinder. The transition process of the wake instabilities will be numerically predicted by two-dimensional simulation based on finite volume discretization with PISO algorithm. The Strouhal-Reynolds relation will be made and be compared with the data from some of the remarkable studies in the past decades to have the certification made. The shear flow is also managed to be studied by controlling the velocity profile of the encountering flow with certain velocity jump. Instead of applying the inlet flow with linear variation, the velocity profiles at the upper half and at the lower half of inlet boundary are both uniform but different in speed, due to which the shear flow would be naturally grown and more real. The flow phenomena are going to be studied within a controlled range of moderate Reynolds number from 50 to 500 where the flow remains unsteady with alternating vortices shedding to the downstream, and then achieves the chaotic

state before which the periodic, the quasi-periodic state, the change of dominating frequency, and so on could be numerically observed and discussed.

Chapter 2

Mathematical Model

2.1 Physical Model

A fixed square cylinder confronted by a constant free stream has been studied. The velocity profile of the free stream is divided into two parts with different values of velocity to have the shear flow effect studied, which is depicted in Fig. 2-1. The model is described by Cartesian coordinate system with x-axis being aligned with the direction of free stream and the origin located at the center of the cylinder. The average velocity of the free stream, V , and the side length of the cylinder, L, are the characteristic velocity and length in the study, respectively. The blockage ratio, β =L H, is 0.1 which means that the vertical width of the computational domain, H, is defined as 10 in unit length. The default L1 is 6, whereas it will be changed to 10 for further studies.

2.2 Governing Equations

The study is considered to be two-dimensional with the assumption of a laminar Newtonian working fluid under unsteady and incompressible conditions without body force. The heat transferring phenomenon does not affect the fluid field being studied under the isothermal condition.

Continuity equation 0 = ∂ ∂ i i x u (2.1)

Navier-Stokes equation (Momentum equations)

(

)

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ i j j i j i j i j i x u x u x x P x u u t u ρ _μ ρ (2.2)

designated by ρ, μ, P, u, respectively.

2.3 Dimensionless Groups

The dimensionless quantities are introduced as,

L tV t∗ = L x x∗ = , L y y∗ = , V u u∗ = , V v v∗ = , ₂ V P P ρ = ∗

,where both L and V represent the characteristic length and characteristic velocity of the domain of interest, respectively.

The governing equations can be written in dimensionless form as

0 = ∂ ∂ ∗ ∗ i i x u (2.3)

( )

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ i j j i j i j i j i x u x u x x P x u u t u Re 1 (2.4)

,where the dimensionless number Reynolds number is defined as

ν

VL =

Re .

2.4 Boundary Conditions

Inlet: the velocity profile is divided into two equal parts with identical uniform

shapes but different values in velocity, which is depicted in Fig. 2-1. Due to the average value of inlet velocity is always expected to be the characteristic velocity, the velocity values of half upper and half lower part are set to be

α

+ = 1

u and u= 1−α , respectively, with v=0.

Outlet: convective boundary condition =0 ∂ ∂ + ∂ ∂ x u t c φ φ [42],

where φ represents the transported property and uc is the convective velocity.

Cylinder surface: Non-slip boundary condition (u =0, v=0) is applied at the wall surfaces. Free Boundaries: =0 ∂ ∂ y u ,v=0

Chapter 3

Discretization

3.1 Introduction

Finite volume approach with unstructured-grid arrangement is employed to discretize the two-dimensional unsteady incompressible Navier-Stokes equations with implicit differencing in time, which is allowing the bigger time step for unsteady problem without losing stability, and the second order central differencing in space.

In the present study, the Navier-Stokes equations can be simplified to a general form as

( )

V

(

)

P t +∇⋅ =∇⋅ Γ∇ −∇ ∂ ∂φ _{ρ φ φ} ρ v (3.1) where φ designates each velocity component.

The equation is also presented with volume integral form as

( )

_∫∫∫

(

)

_∫∫∫

∫∫∫

∫∫∫

∀ ∀ ∀ ∀ ∀ ∇ − ∀ ∇ Γ ⋅ ∇ = ∀ ⋅ ∇ + ∀ ∂ ∂ Pd d d V d t ρ φ φ φ ρ v (3.2)

integrated over a control volume.

3.2 Finite Volume Method

The numerical methods commonly used for fluid dynamics are finite difference, finite volume, finite element, and boundary element methods. Before applying one method for the work, the form of grids is chosen to be unstructured to well model the extremely complex geometries. The troubles of coordinate transform for finite difference are, hence, avoided. Finite volume method is applied in this study because it is more fluid dynamics than finite element method is and has less mathematical treatments than boundary element method

does. The smoothness assumption is avoided due to finite volume is based on the integral form of equations. Despite the physical phenomena are not always coming smoothly, the solution with high gradient can be well performed and the conservation property is preserved as well. Dealing with the integral form of the transport equation, the terms of the equation are known as unsteady and convection term at the left hand side while diffusion and pressure term at the right hand side of the equation. The further deduced control surface approach by Gauss divergence theorem simplifies each term other than unsteady term as

∫∫

∫∫

∫∫

∫∫∫

∀+ ⋅ = Γ∇ ⋅ − ∂ ∂ ∀ S S S S Pd S d S d V d t v v v v_{φ φ} ρ φ ρ (3.3) Unsteady term

The volume integral of the unsteady term

∫∫∫

∀ ∀ ∂ ∂ d t ρφ

can be discretized on the center of each control volume with first-order approximation as

(

n o

)

t φ φ ρ − Δ ∀ Δ (3.4) where the superscripts n and o of propertyφ are denoting the new and old time

levels, respectively. The old time level term would be put into the source term.

Convection term

The surface integral of this term

∫∫

⋅ S S d Vvφ v ρ

was approximated on the entire surface of one control volume (cell) as

(

)

∑

∑

= ⋅ f n f f f f f C f V S F ρ v v φ (3.5) with subscript f denoting that the properties on the surfaces of a control volume

linear combination of upwind and central difference schemes and can be written as

[ ]

[

(

_UD

)

]

o f CD f n UD f C f F F F F = + γ − (3.6) in which the second term at the right hand side of the equation is moved into the source term, which is explicit. In the present study, the central difference scheme is mainly applied via γ =1

The first-order upwind difference approximation is

nb f P f UD f m m F φ ⎟φ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ =max • ,0 min • ,0 (3.7) as well as the second-order central difference approximation is

(

)

[

p P p nb

]

f CD f m f f F = − φ + φ • 1 (3.8) where fpdenotes the weighting factor. The subscripts P and nb denote the

properties on primary and neighboring control volumes (cells), respectively.

Diffusion term

The surface integral of this term

∫∫

Γ∇ ⋅

S

S dv

φ

was approximated over the entire control volume surface as

(

)

∑

∑

= ∇ ⋅ f f n f f f D f S F μ φ v (3.9) whereSf d

(

Sf d

)

v v v v _{= + −}

. dvis the vector pointing in the direction from primary volume (cell) center to the neighboring volume (cell) center (Fig. 3-2). The length dv was considered to be the factor affecting the diffusion dominancy and numerical stability. Hence, the over-relaxed approach of dv was introduced. [39]

, 2 Pnb f Pnb f S S d δ δ v v v v v ⋅ = (3.10)

The diffusion flux approximated by the chosen scheme over-relaxed approach is rewritten as

(

)

(

S d

)

S S F _nbn _Pn _f o_{f f} f Pnb f f D f v v v v v − ⋅ ∇ + − ⋅ = φ φ μ φ δ μ 2 (3.11) in which the second term at the right hand side of the equation is also moved into the source term, which is explicit.

Pressure term

The surface integral is

∫∫∫

∫∫

∀ ∀ ∇ = Pd S d P S r (3.12) which is approximated as ∀ Δ ∇ =

∑

o f f o f S P P v (3.13)

Arrangement of the difference transport equation

The volume integral form of the transport equation was approximated to a linear algebraic equation also known as

∀ Δ ∇ − + =

∑

o nb n nb nb n P P A S P A φ φ (3.14) in which t A A nb nb P Δ ∀ Δ + =

∑

ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− + ⋅ = max • ,0 2 f f Pnb f f nb m S S A v v v δ μ

(

)

[

]

(

)

{

}

o f f o f f o UD f CD f t d S F F S γ μ φ ρ φ Δ ∀ Δ + − ∇ + − − =

∑

v v

where the source term S includes the second terms of unsteady, convection, and diffusion terms.

Chapter 4

Numerical Method for Unsteady Flow Calculation

4.1 Introduction

The governing equation was approached implicitly by finite volume method but the calculation efficiency is expected to be promoted by ignoring the conventional iteration procedure. A non-iterative method named PISO which stands for Pressure-Implicit with Splitting of Operators proposed by Issa 1984 [40] to numerically analyze the implicitly discretised time-dependent fluid equations is being applied. The tedious iteration is removed with the advantaged large time step of implicit differencing remained. The process of solution was split into a series of predictor and few correctors of momentum equations at each time step with an accurate approximation well achieved. During the process, the discretised algebraic momentum equations would be implicitly solved alone to temporally predict the values of velocity at first, of which this moment therefore is termed the predictor step. The accurate solution would never be approached without continuity equation being involved. The mass conservation law is then fulfilled at the following few corrector steps by solving the algebraic pressure correction equations derived from the algebraic momentum equations combined with mass conservation law. The concept of the successive correctors of splitting operators is similar to an iterative method SIMPLE [41]. SIMPLE has to conventionally go all over the iteration steps to have the solution approached at a single time step but PISO saves the time from the unnecessary repeats. This breakthrough of PISO is attributed to the few more correctors at each time step while SIMPLE has only one corrector proceeded at each iteration step. Even the first corrector satisfies continuity equation, mass

conservation has not fully obtained yet [40]. PISO method possesses good temporal accuracy and good stability for large time step which additionally profits the approaches of steady-state problems as well due to the fast and efficient calculation speed.

4.2 Details of the PISO Algorithm

Predictor step

The predictor step is to implicitly solve the algebraic momentum equation deduced by finite volume the task of previous chapter using the prevailing pressure field.

(

−∇ Δ∀

)

+ =

∑

∗ ∗ o P nb nb nb P PV A V S P A v v (3.14) which is thus solved to yield the Vv∗field but the mass conservation law has not satisfied yet. The corrector steps are then taking care of the mass conservation law of the flow field by updating the corresponding pressure.

First corrector step

The new velocities and the corresponding new pressure assumed to be obtained from the very first corrector step are denoted with superscript∗∗ and∗, respectively. The momentum equation is taken as

(

−∇ Δ∀

)

+ = ∗ ∗ ∗ ∗

∑

P nb nb nb P PV A V S P A v v (4.1) The first velocity corrector is certainly the difference between these two algebraic momentum equations and deduced as

f f P f P A V _⎟⎟ ∇ ′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ′ v (4.2) whereV′≡V∗∗ −V∗, o P P P′≡ ∗ −

f f f f V S m&′ =ρ v′ ⋅ v _{f f} f P f P S A v ⋅ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ρ (4.3)

Over-relaxed approach is employed to letSf d

(

Sf d

)

v v v

v _{= + −}

for better numerical diffusion control.

(

)

P

(

S d

)

A P S S A m f f f P f Pnb f f Pnb f f P f f v v v v v v & _⎟⎟ ∇ ′ ⋅ − ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ′ ∇ ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ′ δ ρ δ ρ 2 (4.4)

Then, replacing the term Pf δPnb

v ⋅ ′

∇ in the bracket byPnb′ −PP′, and summing up the

mass flux correction equations all over the control surfaces to completely achieve the properties on one control volume with one equation goes to

(

S d

)

P A m P A P A _{f f} f P f f f f nb P nb P P P v v ₋ ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ′ + ′ = ′

∑

∑

• ρ (4.5) where

∑

= f P nb P P A A f Pnb f f P f P nb S S A A _{v v} v ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = δ ρ 2

∑

∑

∑

_{′ =} ∗∗ ₋ ∗ f f f f f f m m

m& & &

Mass conservation is directly achieved by assuming

∑

∗∗ =0

f f

m& over each control volume. The first pressure correction equation is deduced as

1 2 1 1 P P f nb P nb P P P P A P S S A ′ =

∑

′ + + (4.6) where

∑

∗ − = f f P m S11 &

(

)

∑

_⎟⎟ ∇ ′⋅ − ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = f f f f P f P P S d A S12 ρ v v

More corrector steps

The new velocities and the corresponding new pressure are denoted with superscript∗∗∗ and∗∗, respectively. The momentum equation is taken as

(

−∇ Δ∀

)

+ = ∗∗ ∗∗ ∗ ∗ ∗

∑

P nb nb nb P PV A V S P A v v (4.7) The difference between the momentum equation right above and the previous one with velocity denoted with∗∗is deduced as

f f P f P f nb nb P P A A V A V _⎟⎟ ∇ ′′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ = ′′

∑

v v (4.8) where V ′′≡V∗∗∗ −V∗∗,V′≡V∗∗ −V∗,P′′≡P∗∗ −P∗

Again, the mass flow rate corrector, hence, is

f f f P f P f nb nb f f P S A A V A m v v & ⋅ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ′′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ = ′′ ρ

∑

(4.9)

Over-relaxed approach and ∇Pf′ ⋅δPnb = Pnb′ −PP′

v

are applied. By following the same procedure of first corrector, the present corrector is

(

)

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⋅ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ − − ⋅ ′′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ′′ + ′′ = ′′

∑

∑

∑

_f f P f nb nb f f f f P f f f f nb P nb P P P S A V A d S P A m P A P A v v v v & ρ ρ (4.10) where

∑

= f P nb P P A A f Pnb f f P f P nb S S A A v v v ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = δ ρ 2

∑

∑

∑

_{′′ =} ∗∗∗₋ ∗∗ f f f f f f m m

m& & &

volume while

∑

∗∗

f f

m& =0 has been done from the assumption of the first corrector. The present corrector, therefore, is

2 2 2 1 P P f nb P nb P P P P A P S S A ′ =

∑

′ + + (4.11) where f f P f nb nb f P S A V A S v v ⋅ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ − = ρ

∑

2 1

(

)

∑

_⎟⎟ ∇ ′′⋅ − ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = f f f f P f P P S d A S22 ρ v v

Although more corrector steps are needed for the conservation law, two corrector steps are quite sufficient to have the accuracy of solution within temporal truncation error [40].

The source term of each algebraic pressure correction equation, i.e. corrector, has been divided into two parts. The first one is accounting for the orthogonality of structured girds, whereas the second one is for the non-orthogonality of unstructured grids. Namely, the first one is sufficiently needed if the gird arrangement is regular and rectangular. On the other hand, the second one is also required if the grid arrangement is irregular and distorted. Each pressure correction equation, i.e. corrector, hence, is solved by two steps with the first and the second source term separately. The correcting calculation goes through the first step once and then the second step of which the further repeats are necessary if the grids are quite distorted. Due to the grid arrangement of present study is structured, surface vectorSf

v

parallels the vectorδPnb

v

, i.e. Sf

v

equalsdv, only the first step is gone through.

4.3 Implementation of Boundary Conditions

The inlet boundary condition locating at the upstream of the square cylinder is set to be a shear flow. The inlet is equally divided into an upper part and a lower part. These two parts possess the uniform flow with different velocities from each other to create a jump in between. At the outlet one locating at the downstream of the cylinder, the convective boundary condition

0 = ∂ ∂ + ∂ ∂ x u t c φ φ

[42] a wave equation is applied. The implicit method is used as

0 = Δ − + Δ − x u t n C n B c o B n B φ φ φ φ (4.12)

The boundary and the neighboring inner cell properties are denoted with the subscripts B and C, respectively. The convective velocity uc is set to be the

local flow velocity at the outlet boundary. The equation is also shown as

Cr Cr Cn o B n B + + = 1 φ φ φ (4.13) in which x t u Cr c Δ Δ

= is defined as Courant number.

At the upper and lower site boundaries of the entire computational domain, the boundary velocities in the x-direction are set to be identical with the x-velocities at the immediate neighboring cells inside the domain along y-direction ( =0

∂ ∂ y u

) and v=0 which are both utilized to have the free boundary condition applied. At the surface of the square cylinder, the non-slip boundary condition on wall is applied.

4.4 Solution Procedure of PISO

Step 1: read the velocities and pressure of the flow field from the old time level. Step 2: solve the momentum equation (3.14) to get V∗.

Step 3: compute P′ by solving the first pressure correction equation (4.6) to update velocities and pressure to get V∗∗and P∗.

Step 4: compute P ′′ by solving the second pressure correction equations (4.11) to further update velocities and pressure to get V∗∗∗and P∗∗.

Step 5: if the required time step is achieved, then stop the calculation and output

the data otherwise proceed to the next time step and repeat all over the way from step 1 to step 4.

Chapter 5

Validation of Numerical Code

5.1 Grid Test

The computational domain is discretized into a number of rectangular cells with different sizes. The grid distribution is not uniform throughout the domain in which the denser part is located close to the square cylinder. It is depicted in Fig. 5-1. The flow pattern is expected to be complex at the near wake of the cylinder.

The grid test has been done with 198×97 and 232×123 cells in x and y directions, respectively. The relations between Strouhal number and Reynolds number obtained by these two cases are depicted in Fig. 5-2. The case with 232×123 cells presents a smoother result and is applied for the present study.

5.2 Data Comparison

The data of the present numerical prediction with uniform free stream has been validated by making the comparison of the relations between Strouhal number and Reynolds number with those reported by Okajima [13], Davis et al. [17], Franke et al.[44], Sohankar et al.[19], and Saha et al.[43], which is depicted in Fig. 5-3. The Strouhal frequency is detected from the transverse velocity of a monitor point located at (x=1.23, y=0) in the near wake of the cylinder. The detected signals are identical if x<5 [45]. The range of Reynolds number for validation is from 50 to 500. During these predictions, some interesting phenomena of the Strouhal frequency were observed from the power spectrum analyses. The relation of St-Re shows a big valley from Re=175~375. With starting from the deep bottom of the valley, i.e. Re=250, and then going all

over the way up to the hill along the increase of Reynolds number, the period doubling which means that the second peak of lower frequency emerges shows up at the spectra (Fig. 5-4) while the starting point was found to be 218 and 250 by Saha et al. [45] and Sohankar et al. [19], respectively, although there is also the third peak which is still small. The frequency of this additional peak is somehow higher than the half value of Strouhal frequency from Re=250, but in the range of around Re=375~425, the frequency is exactly the half value of Strouhal frequency. This 0.5 frequency-locking mode was also reported in the range Re=325~375 by Saha et al. [45]. The clear spectrum are shown in Re=375~425, the range of half frequency locking, but before and after which the noises are quiet strong as Re>300. About Re>450, fS

2 1

possesses the highest peak at the spectra and dominates. This feature was also found at Re=500 by Sohankar et al. [19]. As Re=500, the chaotic state is determined by the space diagram (Fig. 5-5).

The mean drag and mean lift coefficient of the forces acting on the square cylinder are also calculated. These mean coefficients are got from the time average values over a certain range of time. The drag and lift coefficients, CD

and CL, are defined as FD V L

2

2 ρ and FL V L

2

2 ρ , respectively.

The comparison of the mean drag coefficient has been also made with the data of Davis et al. [17], Franke et al.[44], and Saha et al.[43], depicted in Fig. 5-6. The value of the mean lift coefficient is starting to deviate from zero as Re>325. On the other hand, C_L =−0.04 was found at Re=400 by Sohankar et al. [19].

5.3 Time Step Size

The test of step size marching in time is also done by making the comparison with the data reported by Saha et al. [43]. The time step size is chosen to be 0.01 or even the finer one 0.005. Applying the data by Saha et al. for the validation is because the blockage ratio was 0.1 and the square cylinder was 6 in unit length away from the inlet boundary of the computational domain, which are identically applied for the uniform free stream case of present study. The comparison is depicted in Fig. 5-7 which shows the relations of Reynolds number and Strouhal number obtained from the computations with time step size 0.01 and 0.005 as well as the data by Saha et al.. The present computations with different step sizes are quite consistent with each other in the whole range of various Reynolds numbers under consideration. The time step size 0.01 is, therefore, chosen for the present study to shorten the calculation time.

Chapter 6

Results and Discussion

Shear flow has been studied by changing the velocity ratio of free stream and the distance between the square cylinder and the inlet boundary of the computational domain, which provides different shear rate, the velocity gradient confronted by the square cylinder. The parameter α is considered to be 0.5 and 0.2 to have the velocity ratio of free stream 3:1 and 3:2, respectively, whereas the parameter L1 is 6 and 10 to have different distances for shear flow to grow, see Fig. 2-1, with Reynolds number from 50 to 500 in the present study. The shear rate effect can be studied with wide range of Reynolds number.

6.1 Relation of St and Re with Shear Free Stream

The relation of St-Re is depicted in Fig. 6-1 with Strouhal number St,

V L

f_s , and Reynolds number Re, VLν . The similarities between uniform and shear flow show up as Re<100 but the discrepancies appear as Re goes high. The St-Re relations of shear flows even with different velocity ratios of free stream and L1s are showing the St-drop about Re>150, which might indicates the transition process to three-dimensionality. Around Re=300, the chaotic states are observed by the phase space analyzes (Fig. 6-2 ~6-5), before which the interesting feature of Strouhal frequency is found with different shear rates. The higher velocity ratio generally brings the lower Strouhal frequency as L1 remains the same, which was also reported by Saha et al. [46] and Kang [47]. On the contrary, the longer L1 possessing the lower shear rate of the more developed shear flow brings the lower Strouhal frequency, whose discrepancy may be affected by the different transverse spans of shear layer naturally grown.

As Reynolds number goes beyond around 325, the shear flow frequency spectrum, see Fig. 6.6 ~6.9, become so confusing. The main Strouhal frequencies are quite hard to be defined for the shear flow cases besides the case of velocity ratio 3:1 and L1=6. Hence, Instead of defining the Strouhal frequencies, the dominating frequencies are applied for Fig. 6.1 as Re>325.

6.2 Period Doubling

The signals in the near wake of square cylinder were detected by the transverse component of velocity at point (x=1.23, y=0) and were analyzed by Power spectrum analysis. The signals detected from the flows with shear free stream are also repeating some of the phenomena observable in the uniform free stream case, especially the cases of shear rate 3:2 (Fig. 6-7 & 6-9). Once the velocity ratio goes higher to 3:1, the process of period-doubling somehow disappears. The periodic state lasts to higher Reynolds number as well as the disorder vibrating motion of the wake, hence, comes along with lower Reynolds number. The period-doubling process becomes rarer even with the longer L1. The most interesting phenomena happen in the case of velocity ratio 3:2 and L1=10. As Reynolds number goes beyond around 375, the period-doubling is so obvious no matter how big the Reynolds number is. The dominating frequencies are all 0.10376 Hz.

6.3 Lift Coefficient

The coefficient of the lifting force acting on the square cylinder has been numerically predicted. As the free stream is uniform, the mean lift coefficient should be zero theoretically unless Reynolds number exceeds a certain number, and it was well studied without too many discrepancies between the data from researchers. It is unfortunately that as the free stream is shear flow, the direction of mean lifting force acting on the cylinder becomes an issue. Two opposite directions were observed and predicted with similar physical assumptions. Moreover, the direction of the mean lift could be even reversed at certain Reynolds number, which was reported by Kurose and Komori [48]. Hence, the mean lift coefficient under a shear flow condition is an important topic for the present study.

The mean lift coefficients from cases in the present study were obtained by the average values over a period of time in which the wake flow was alternatively shedding in a stationary periodic state. Despite the chaotic state might be reached with higher Reynolds number, the time-averaged values of the mean lift coefficients were still calculated by the identical time-span even the periodic motion was not observed.

In the present study, the direction of the mean lifting force is pointing toward the higher velocity part, i.e. pointing up, all over the whole span of Reynolds number from 50 to 500. Known from Fig. 6-10, the mean lift coefficient does not always go up along the increase of Reynolds number. After a certain value of Reynolds number, the mean lift goes down. The good agreement between each other is achieved by the cases of shear rate 3:2. Even the parameter L1 is not the same, the both Lift-Re relations of L1=6 and L1=10

are quite coinciding with each other before Reynolds number reaches 325. As Re>325, the discrepancies and the fluctuating motions are observed and considered be triggered by the chaotic state. The mean lift coefficient dramatically rises once the velocity ratio is increased to 3:1. Not only does the discrepancy happen between the 3:1 cases of L1=6 and L1=10, but the mean lift coefficient of L1=10 also keeps going up until the chaotic state is approached, which is quite diverse from the other three shear flow cases in the present study. The turning points of the mean lift-Re relation is at about Re=225 for the cases of velocity ratio 3:2 while it is at about Re=175 for those of velocity ratio 3:1. The fluctuating motion of mean lift-Re relation is also observed as the flow of velocity ratio 3:1 reaches the chaotic state about Re=300. Especially the one with L1=10, the range for the values of mean lift coefficient to go up and down is almost a unit order. The one with L1=6 fluctuates but as the value of Reynolds number goes beyond 350, the decrease of mean lift coefficient is present until about Re=475.

The explanation of this phenomenon may be attributed to the theory proposed by Bernoulli. The distinguishing feature of present study is that although the shear flow was applied, the shear flow has not been fully developed yet as well as the shear layer did not cover the whole square cylinder in transverse direction. The shear flow effect seems to dominate less. Therefore, the upward lift force was caused because the lower pressure field was formed by the faster free stream on top of the cylinder.

6.4 Drag Coefficient

The mean drag force acting on the square cylinder is studied and presented by the relation of the mean drag coefficient and Reynolds number, which is depicted in Fig. 6-11. The mean values of them were also calculated from the average values over an identical time-span with the one applied for the mean values of lift coefficient. The study of mean drag coefficient presents neither an issue nor the big discrepancy, but the power of it may be related to the velocity ratio. The possible damages of a bridge or of other structures encountering the shear flow can be imagined by this investigation. The mean drag coefficients of the cases with higher velocity ratio are obviously much higher than those with lower one whereas the parameter L1 seems to affect less in this study. The good agreement of mean drag-Re relations with different values of L1 shows up as the velocity ratio is remained the same. The dissimilarities present after the chaotic state was achieved. The comparison is also made with the case of uniform free stream. As the velocity ratio is staying 3:2, the drag forces are just a little bit higher than those of the case with velocity ratio 3:3, i.e. uniform free stream. The coincidence between 3:2 and 3:3 of velocity ratio with L1=6 is even approached when the values of Reynolds number are around 400. As the velocity ratio goes up to 3:1, the drag forces become much stronger.

6.5 Instantaneous Vorticity Contours

The instantaneous flow structure certainly changes with the increase of Reynolds number. During this change, the movement of the flow is going from period to quasi-period, and then to the chaotic state, whose corresponding flow patterns of uniform and shear free stream are all depicted. In order to tell the differences between the flow patterns with different values of Reynolds number precisely, the instantaneous flow fields have to be caught under the same condition. The lift coefficient is applied to be the criterion for the present study. All the flow structures being displayed with different Reynolds number are caught as the lift coefficients reach the corresponding mean values. The vortices alternatively shedding from the square cylinder to the farther downstream are also depicted.

6.5.1 Uniform Flow

The instantaneous vorticity contours are depicted in Fig. 6-12 with different values of Reynolds number. Due to the uniform free stream, the strengths of the positive and negative vortices shedding from the bottom and the top sides of the square cylinder, respectively, are quite compatible. At Re=50, the vortices are elongated to about 9 in unit length and they become much more rounded in shape at about Re>75. The vortex farther away from the cylinder to the downstream is more rounded and also dissipated by the mechanism of momentum diffusion. With the increase of Reynolds number, the transverse moving range of vortices is becoming wider. Before Reynolds number reaches 300, the vortices in the wake flow are well-behaved leading the succeeding one to the downstream step by step. One positive vortex leads another negative one and vise versa. The distance between two vortices is well maintained as well. As

Reynolds number goes beyond about 300, the vortices start to move sort of arbitrarily. The distances are no more maintained. The positive vortex and the negative one are so against each other. As Reynolds number moves further to about 450, the vortices are even being piled up and clustered together in any possible direction.

6.5.2 Shear Flow

All the shear flow fields in the present study possess the negative background vorticity because the higher value of free stream velocity was arranged to be located at the upper half part of the inlet boundary. The negative vortices would be made at where the shear layer grows. The strengthened clockwise vortices, i.e. negative, as well as the weakened counterclockwise ones, i.e. positive, are hence shown in the shear flow cases (Fig. 6-13 ~6-16). The elongated positive vortices are also observed due to the positive vortices are shed from the square cylinder’s bottom side where the velocity of the free stream is slower. As the positive vortices are shedding to the downstream, they are then confronting the faster free stream from the top and elongated. On the contrary, the negative vortices shedding from the top of the square cylinder are even more rounded because of the lower free stream from the bottom. With increasing the velocity ratio, the phenomena just mentioned are more obvious. The positive vortices are longer and thinner while the negative vortices are more rounded and bigger. On the other hand, the distance between the vortex and the succeeding vortex becomes longer with the higher velocity ratio.

As Reynolds number stays low, the counterclockwise and clockwise vortices are extending to the farther downstream continuously, especially the counterclockwise ones with solid lines. Due to the convection dominates much less at Re=50; the clockwise vortices are not rounded yet until Re=75. Along the

increase of Reynolds number, the vortices break into pieces and occupy bigger range in the transverse direction by the wider shear layer in which the velocity gradient is lower from the bottom to the top. Similarly, the farther downstream engages the vortices with the wider span in transverse. The dissipation of vortices at the downstream of the wake is also observed. The clockwise vortices dissipate more rapidly than the counterclockwise ones do. The higher velocity ratio, 3:1, elongates the counterclockwise vortices in the near wake, but whose succeeding counterclockwise vortices about to grow at around 10 unit length away from the cylinder are easy to be torn apart.

6.6 Instantaneous Streamlines

The instantaneous flow patterns of each case with the increase of Reynolds number are described by streamlines depicted in Fig. 6-17 ~6.21. These instantaneous flow fields are also caught as the instantaneous lift coefficients reach the corresponding mean values. These instantaneous streamlines of these cases are corresponding to those illustrations of the instantaneous vorticity contours mentioned in chapter 6.5.

6.6.1 Uniform Flow

As the free stream is uniform, one eddy right behind the square cylinder can be observed in Fig. 17. As the Reynolds number stays low, the curvatures of streamlines at the downstream of the wake flow are small and smooth. As the Reynolds number goes higher, the wake fluctuates more with bigger curvature. The streamlines are also illustrating the phenomena of alternative motion in the wake of a square cylinder by showing the wave-like motion. As the Reynolds number goes higher, the distances between waves are getting closer until around Re=250 beyond which the wave-like motion starts to go arbitrarily and influence

the flow field wider in the transverse direction.

6.6.2 Shear Flow

Fig. 6.18 ~6.21 are illustrating the flow patterns of the flow fields with shear free stream. Two different velocity ratios and distances between the cylinder and inlet boundary are also considered. The shear free stream with higher velocity at the top side brings another eddy right under the square cylinder and even some smaller eddies in the farther wake. The eddy adhere to the bottom of the cylinder may grow up as the velocity ratio goes high, while on the contrary, it may turn smaller as the distance between the cylinder and the inlet boundary elongates. The shear rate of the present cases rises while the velocity ratio goes up as well as the inlet-cylinder distance is shortened, hence, this eddy grows. The biggest eddy formed beneath the square cylinder can be observed in the cases of velocity ratio 3:1 and L1=6.

6.7 Time Evolution of Flow Fields

Fig. 6-22 ~6-25 are showing the time evolution of an uniform free stream passing over the square cylinder with Reynolds number 100, 200, 300, and 400. The instantaneous flow fields evolve through a period obtained from the corresponding Strouhal frequency. The strength of vorticity dominates the behaviors of eddy behind the cylinder. As the negative vorticity occupies bigger area immediate behind the cylinder than the positive one does, the dominating eddy shown by streamlines goes clockwise. Once the area occupied by the positive vorticity catches up the one occupied by the negative vorticuty, the counterclockwise eddy grows while on the contrary the clockwise one diminishes and vice versa. The positive and the negative vorticities are alternatively shedding from the cylinder and regularly dominating the flow field

in turn, whose phenomena are being describing by the instantaneous streamlines. Moreover, the weaker positive vorticity at the farther downstream rolls up streamlines while the negative one rolls down them. As Re=300, the counterclockwise eddy at 5T/6 even breaks the rule to grow up competing against the clockwise one. The time-evolving instantaneous flow fields with shear free stream are also studied. In order to see the big difference to the uniform free stream one, the cases are selected with velocity ratio 3:1 and L1=6 (Fig. 6-26 ~6-29). Due to the elongated positive vorticity and the rounded negative one, the streamlines in the wake show the smoother concaves and the rugged convexes. Furthermore, although the positive and the negative vorticities dominate in turn, the positive vorticity is always eager to grow before the negative vorticty diminishes, because of which, the counterclockwise eddy appears with the still big clockwise one, and hence is pushed beneath the cylinder. As the clockwise eddy becomes smaller, the counterclockwise one goes behind the cylinder. As Reynolds number is 200 and 300, the negative vorticity is even strong at the farther downstream which can be told by the clockwise eddy of the streamlines. As Re=400, both the positive and the negative vorticities are grown strong and so effecting the farther wake dramatically and chaotically.

6.8 Elongation of Computational Domain

The computational domain is elongated from 24 to 32 in x-direction to have the studied more certified. The Strouhal frequency and the mean drag coefficient are again obtained with longer computational domain to make the comparison with the default assumption (Fig. 6-30 & 6-31). The good agreement is shown.