A Numerical Study of the Dynamics of Parachute Aided Descent

A Numerical Study of the Dynamics of Parachute Aided Descent

R. Darshankumar Department of Mechanical Engineering

Coimbatore Institute of Technology

Coimbatore, India

Abstract- Computational Fluid Dynamics (CFD) has been successfully applied to many engineering problems. This paper presents the influence of three different two-equation turbulence models on computing the drag coefficient for both hemispherical canopy and Ram air canopy at different velocities. The aim of this paper is to offer guidance to select appropriate turbulence models for this application. Using available different turbulence models from ANSYS FLUENT the drag coefficient of hemispherical canopy (HC) and Ram air canopy (RC) of different sizes were calculated in the range of velocity (2 m/s -50 m/s). These results were used to evaluate the terminal velocity and time taken to reach the ground by the paratrooper by solving the differential equations numerically using RK4 technique. Comparisons have been made for both the canopies based on the effects of shape and also the area of the canopies. It was observed that all three models predicted the drag coefficient to be a constant in the above mentioned velocity range and results also indicate that the canopys shape does not influence the terminal velocity.

Keywords: Turbulence Models, Parachutes Canopy, Drag Coefficient, Terminal Velocity.

1. INTRODUCTION

   Skydiving is an interesting sport and an amazing airborne technique used in military operations. Though named differently the physics involved is Newtons second law of motion which is mathematically expressed as

   = (1)

   Where F denotes the net force in the system, m denotes the mass of the object and a denotes the acceleration of fall. Parachutes are used to land safely from high altitude free fall. Basically there are two types of parachute in which one is the dome canopy and the other ram-air canopy model. In this project both these canopies are studied under no cross wind assumption. In the dome canopy the air is trapped in the canopys envelope which create a high pressure that retards movement in the direction opposite to the entering air flow whereas in the ram-air canopy model, the parafoil acts as a wing, allowing the jumper to fly towards a target.

2. PHYSICS INVOLVED

   Newtons Second Law: The acceleration of an object as produced by a net force is directly proportional to the magnitude of the net force, in the same direction as the net force, and inversely proportional to the mass of the object. As mentioned before it is mathematically expressed as (1).

   The force term (F) consists of both the paratroopers weight

   1. and the air resistance (FD) which opposes the fall. Skydiving consists of two phase namely,

      • Free Fall,

      • Parachute Aided Fall.

        Free fall is the period which starts from the point of jump and end when the parachute is deployed and the remaining period where the descent is using parachute is called Parachute aided fall. At the point of jump the paratroopers velocity would be zero and he falls the velocity will increase rapidly and at some point the velocity becomes independent of time. That velocity is known as Terminal velocity. This happens when the weight (W) and drag force (FD) becomes equal. It is known that the weight always remains constant and the drag force changes with velocity and when the velocity is such that the magnitude of drag force equals the weight then the net force (F) acting on the system would be zero.

        = 0 (2)

        = (3)

        Drag force (FD) takes into account the mass, shape of the object and it is mathematically expressed as

        = 0.5 2 (4)

        Where

        = (3),

        = (2),

        = ,

        = ().

        The paratrooper attains terminal velocity during both the

        phases if the fall is sufficiently large. Initially it is important to know the drag coefficients of both the canopies and the paratroopers body if he were to fall tummy facing the ground.

3. TURBULENCE MODELS

   Turbulence is an irregular motion which in general makes its appearance in fluids, gaseous or liquids, when they flow past solid surfaces or even when neighbouring streams of the same fluid flow past or over one another.

   Characteristics of Turbulent flow

   • Non-uniform,

   • Inertial dominant flow,

   • High Reynolds number,

   • Unsteady,

   • Presence of eddies,

   • Large length scale between eddies,

   • Large time scale between eddies,

   • Highly diffusive,

   • Highly dissipative,

   • Flow will be 3 dimensional.

   Turbulent flows are highly complex such that the physics of the problem changes for every moment. Therefore time averaging the Navier-Strokes (NS) equation is important in order to get time averaged results.

   The idea of Reynolds Time Averaging is to express the variable which is a function of space and time as the sum of mean and a fluctuating component.

   Standard , RNG and SST computes eddy viscosity as a function of turbulent kinetic energy (k) and either turbulence dissipation rate () or specific turbulence dissipation rate () respectively.

   Standard Model

   Standard is a two-equation model in which the solution of two separate transport equation allows the turbulent velocity and length scales to be determined. This model is valid only for fully turbulent flows [2, 3]. The turbulence kinetic energy, k, and its rate of dissipation, , are obtained from the following transport equations:

   () +

   () =

   [( + ) ]

   (, ) = (,) + (, ) (5)

   Where,

   (,) = (, ),

   + +

   Y + (10)

   () +

   () =

   [( + ) ]

   (, ) = (, ).

   Definition of time averaging of a property ,

   = lim 1

   (6)

   +1 ( + 3)

   2

   2

   2 + (11)

   Model computes as a function of &,

   0

   2

   Therefore by definition time average of a function independent of time,

   = (12)

   Model constants: = 0.09, 1 = 1.44, 2 = 1.92, =

   = lim

   1

   = ( lim

   1 ) =

   1, = 1.3

   0

   RNG Model

   RNG is the refined version standard .

   And time average of a fluctuating component is,

   The standard version is a high-Reynolds-number model but

   = lim

   1

   = 0

   the RNG version provides an option to also account for low-Reynolds-number effects [2]. This feature makes this

   0

   Time Averaging of NS equations give Continuity equation:

   + = 0 (7)

   model more accurate and reliable for a wider class of flows than the standard version. This model is often recommended for strongly strained turbulent flows [5] and this model also reduces the over production of turbulence around the stagnation points [7]. The transport equations

   for this model are as follows:

   :

   ( + + ) =

   p>

   2 2

   () + ( ) = [ ] +

   + ( + )

   2 2

   2 2

   + Y +

   (13)

   ( + ) (8)

   () +

   ( ) =

   [

   ]

   :

   ( + + ) =

   + (

   + )

   + (2 + 2)

   1

   3

   2 2

   2 R

   +

   (14)

   ( +

   (9)

   2

   )

   The formula for computing in both standard and

   After time averaging the momentum equations, few

   RNG are same for high-Reynolds-number.

   fluctuating terms survive. These fluctuating terms are

   =

   2

   (15)

   known as Reynolds stresses and are related to mean velocity gradients using Boussinesq hypothesis which introduces turbulence or eddy viscosity().

   Model constant: = 0.0845

   SST Model

   Shear Stress Transport combines

   standard and high Reynolds number version of

   =

   the model to capture the physics in the near wall region independent of far field and the physics in the farther portion of the boundary layer independent of near

   =

   0.5

   (23)

   wall [2]. The standard model is formulated into

   in order to achieve the above mentioned purpose. A blending function is used in SST which acts a sort of joint between the standard and the transformed

   and this blending function plays a very interesting part. It is designed that such that it activates the in the near wall region and the transformed far from the wall [2]. This feature ensures that the physics at both the regions are captured properly. The transport equations [6] for this model are as follows:

   Equation (20) gives the terminal velocity attained by the

   paratrooper under free fall. The terminal velocity can be the maximum velocity attained by the paratrooper. Both drag force and drag coefficient are unknown. So it is imperative to find the drag force either experimentally or using CFD tools. For obvious reasons CFD tools are chosen over experimentation to find the drag force on the paratrooper under free fall.

   During free fall the paratroopers body is modelled as a short cylinder. So the paratroopers height is taken as cylinder height (L) and the paratroopers chest

   () +

   ( ) =

   [( + ) ]

   measurement as cylinder diameter (D).

   According to Indian Army the minimum physical standards

   +

   Y

   +

   (16)

   are [4]

   = 166 ,

   () +

   ( ) =

   [( + ) ]

   = 77 ,

   = 48

   The typical measurements are given below

   + + + (17) The turbulent viscosity is computed as

   S.No.	Height, L (m)	Chest, D (m)	L/D	Weight (kg)
   1	1.66	0.77	2.155	62
   2	1.725	0.9271	1.86	69
   3	1.775	1.0026	1.77	73
   4	1.825	1.08	1.69	78
   5	1.88	1.155	1.62	82
   6	1.935	1.23	1.57	88
   7	1.98	1.31	1.51	91

   S.No.	Height, L (m)	Chest, D (m)	L/D	Weight (kg)
   1	1.66	0.77	2.155	62
   2	1.725	0.9271	1.86	69
   3	1.775	1.0026	1.77	73
   4	1.825	1.08	1.69	78
   5	1.88	1.155	1.62	82
   6	1.935	1.23	1.57	88
   7	1.98	1.31	1.51	91

   TABLE I. Typical Measurements of a Male

   = 1

   (18)

   [ 1 , 2 ]

   1

   Where S is the strain rate magnitude and

   = 1

   1 +(11)

   (19)

   ,1 .2

   Simulations were done for velocities 1 / & 100 /

   = 1

   1 +(11)

   (20)

   using the typical physical measurements and the difference between (1 ) (100 ) was about 0.01.

   ,1 .2

   TABLE II. Drag Coefficients for Different L/D values

   L/D	Velocity ()		Difference in %
   1.51	1	0.639	1.9
   	100	0.627
   1.57	1	0.649	1.5
   	100	0.639
   1.62	1	0.652	1.8
   	100	0.64
   1.69	1	0.653	1.8
   	100	0.641
   1.77	1	0.673	1.6
   	100	0.662
   1.86	1	0.682	1
   	100	0.675
   2.155	1	0.736	1.4
   	100	0.726

   L/D	Velocity ()		Difference in %
   1.51	1	0.639	1.9
   	100	0.627
   1.57	1	0.649	1.5
   	100	0.639
   1.62	1	0.652	1.8
   	100	0.64
   1.69	1	0.653	1.8
   	100	0.641
   1.77	1	0.673	1.6
   	100	0.662
   1.86	1	0.682	1
   	100	0.675
   2.155	1	0.736	1.4
   	100	0.726

   Where 1, 2are the blending functions, are given by

   1 = (4) (21)

   2 = (22) (22)

   Model constants: ,1

   = 1.176, ,1

   = 2, .2

   = 1, .2 =

   1.168, 1 = 0.31

4. FREE FALL

Initially after the jump the paratrooper falls at a varying rate starting from 0 / because the drag force and weight of the paratrooper wont be equal. And then the drag force increases and finally becomes equal to the weight of the paratrooper thereby the acceleration goes to zero. At this point the paratrooper attains terminal velocity (Net force = 0).

In all the cases the difference between

(1 ) (100 ) is about 0.01 and thus

can be concluded to be a constant based on L/D.

Cd vs L/D

ANSYS FLUENT (6.3.26). Standard was used in analysing the canopies for grid independency.

TABLE III. Properties of Air at 20oC [1]

0.76

0.72

Cd

Cd

0.68

0.64

0.6

Cd = 0.1331(L/D) + 0.4302

R² = 0.9816

Hemispherical Canopy (HC) Parachute diameter = 8.1048 m Projected area = 51.59 m2

1.5 1.6 1.7 1.8 1.9 2 2.1 2.2

L/D

Cd Linear (Cd)

Fig.1 Variation of with L/D

1. CFD SIMULATION

   The canopies have been designed and meshed in Gambit

   2.4.6 and analysed in ANSYS FLUENT 6.3.26. The hemispherical canopy is designed with a vent hole of diameter 0.305 m. Several hemispherical canopies are analysed but the vent hole dimension remains constant for all the canopies.

   The analysis was done for two different mesh count and they provide results which deviate by 1% which is negligible.

   Velocity m/s	Cd 4 lakh elements	Cd 8.9 lakh elements	Difference in %
   2	1.144	1.158	1.2
   10	1.141	1.155	1.2
   20	1.142	1.157	1.3
   30	1.146	1.163	1.5
   40	1.15	1.167	1.4
   50	1.152	1.171	1.6

   Velocity m/s	Cd 4 lakh elements	Cd 8.9 lakh elements	Difference in %
   2	1.144	1.158	1.2
   10	1.141	1.155	1.2
   20	1.142	1.157	1.3
   30	1.146	1.163	1.5
   40	1.15	1.167	1.4
   50	1.152	1.171	1.6

   TABLE IV. Comparison of Drag Coefficient Values for Two Different Mesh Count (HC)

   1.4

   1.2

   1

   Cd

   Cd

   0.8

   0.6

   0.4

   0.2

   0

   Grid Independency Test

   Fig.2 Gambit model for HC

   Fig.3 Gambit model for RC

2. GRID INDEPENDENCY TEST

The hemispherical canopy and ram air canopy were designed and meshed using Gambit 2.4.6 (3D modelling software) and analysis were done using

0.E+00 2.E+07 4.E+07

Re

Cd (E=4L) Cd (E=8.9L)

Fig.4 Grid independent result for HC

The mesh count used varies by a large number but the results are almost constant and it can be concluded that the drag coefficient for hemispherical canopy is independent of velocity in the range (2m/s -50 m/s).

Ram-air Canopy Span = 4.9 m Chord = 2.35 m Area = 11.59 m2

Projected span = 4.45 m Projected Chord = 1.94 m Projected area = 8.63 m2

TABLE V. Comparison of Drag Coefficient Values for Two Different Mesh Count (RC)

1.5

Cd

Cd

1

0.5

Cd vs Re

0

0.E+00 1.E+07 2.E+07 3.E+07

Grid independency test

Grid independency test

Cd – k-

1.2

1

0.8

0.6

0.4

0.2

0

0.E+00

1.2

1

0.8

0.6

0.4

0.2

0

0.E+00

Cd

Cd

Cd – SST k-

Re

Cd – RNG k-

1.E+07

Re

2.E+07

1.E+07

Re

2.E+07

Cd (E=4L) Cd (E=8L)

Cd (E=4L) Cd (E=8L)

Fig.5 Grid independent result for RC

Similar to hemispherical canopy analysis, using two mesh counts for ram-air canopy the drag coefficients were calculated and concluded to be a constant in the velocity range (2m/s 50 m/s).

VII.PREDICITION OF DIFFERENT TURBULENCE MODELS

RNG , SST turbulence models were used to analyse both the canopies using the same dimensions and properties of air and it was observed that both the models predicted drag forces which were almost same as the drag forces of standard model.

Hemispherical canopy

TABLE VI. Predictions of 3 different Two Equation Models (HC)

Fig.6 Predictions of 3 different 2 equation models for HC

Ram-air canopy

TABLE VII. Prediction of 3 different Two Equation Models (RC)

Velocity m/s Cd

k- RNG k- SST k-

2 1.045 1.042 1.042

10 1.042 1.04 1.04

20 1.042 1.039 1.039

30 1.041 1.039 1.038

Cd vs Re

1.5

Cd

Cd

1

0.5

0

0.E+00 1.E+07 2.E+07

Re

Cd – k- Cd – RNG k-

Cd – SST k-

Fig.7 Predictions of 3 different 2 equation models for RC

The drag coefficient for most of the object remains constant depending upon the mass and projected area when >

104 [8]. As mentioned before the drag coefficients predicted by all three models are same. It is concluded that

working with standard model will be sufficient for this application.

1. NUMERICAL TECHNIQUE

   Runge Kutta 4th order method is used to solve the differential equation involved in this problem to find the terminal velocity and time taken by the paratrooper to reach the ground are given below

   TABLE VIII. Dimensions of Canopies

   Canopy	Diameter m	Span m	Chord m	Area m2
   Dome	4.7	–	–	17.4
   Ram air	–	6.28	2.78	17.4

   Paratroopers dimensions (Cylinder dimensions) Length (L) = 1.935 m

   Diameter (D) = 1.23 m

   2 =

   (24)

   Mass (m) = 117 kg

   2

   2 = ()2 (25)

   Cd (Free Fall) = 0.64.

   Velocity m/s	Drag Coefficient, Cd
   	Hemispherical Canopy	Ram Air Canopy
   2	1.037	1.048
   10	1.035	1.046
   20	1.034	1.045
   30	1.034	1.045
   40	1.033	1.044
   50	1.033	1.044

   Velocity m/s	Drag Coefficient, Cd
   	Hemispherical Canopy	Ram Air Canopy
   2	1.037	1.048
   10	1.035	1.046
   20	1.034	1.045
   30	1.034	1.045
   40	1.033	1.044
   50	1.033	1.044

   TABLE IX. Drag coefficient of both canopies under parachute aided fall

   2

   2

   Let =

   2

   2 = ()2 (26)

   2

   Equation (23) is a second order PDE and therefore it is first converted into two first order ordinary differential equations.

   = () = () (27)

   = (28)

   {(0) = 0

   (0) = 0

   Table IX provides drag coefficients for both the canopies which has same projected area.

   Velocity vs Time

   Velocity vs Time

   40

   30

   20

   10

   0

   40

   30

   20

   10

   0

   Initially let = 0

   Velocity m/s

   Velocity m/s

   1(1) =

   1(2) = ( 2)

   2(1) = ( + 0.5 1(2))

   2

   2(2) = ( ( + 0.51(2)) )

   3(1) = ( + 0.5 2(2))

   3(2) = ( ( + 0.52(2))2)

   0

   0

   100

   100

   200

   200

   300

   300

   4(1) = ( + 3(2))

   2

   Time s

   HC Velocity (m/s) RC Velocity (m/s)

   Time s

   HC Velocity (m/s) RC Velocity (m/s)

   4(2) = ( ( + 3(2)) )

   = (1(1) + 22(1) + 23(1) + 4(1))6

   +1 = + (29)

   = (1(2) + 22(2) + 23(2) + 4(2))6

   +1 = + (30)

   "" is the time step and (29) and(30) can be used to find the distance travelled from the point of jump and values of velocity respectively according to the time step defined after every successful completion of the entire procedure. The above expressions were solved using C++ programming.

2. RESULTS AND DISCUSSION

   When the paratrooper jumps out of the plane his initial velocity would be zero and it increases rapidly until he reaches his terminal velocity in the free fall phase and then once 3/4th of the total distance is covered he deploys his parachute and due to the sudden enlargement in the projected area the velocity reduces quickly. Again after a certain distance the paratrooper attains the terminal velocity.

   Fig.8 Comparison between HC and RC of the entire fall

   Parachute Aided Fall

   Parachute Aided Fall

   40

   30

   20

   10

   0

   40

   30

   20

   10

   0

   105

   110

   115

   120

   105

   110

   115

   120

   Time s

   Time s

   HC Velocity m/s

   RC Velocity m/s

   HC Velocity m/s

   RC Velocity m/s

   Velocity m/s

   Velocity m/s

   Fig.9 Comparison of parachute aided fall for HC and RC

   The terminal velocity under parachute aided fall is 10.35 m/s for HC and 10.25 m/s for RC and also the time taken to reach the ground varies by just one second which is not significant so it can be concluded that the shape doesnt play any significant role in deciding the terminal velocity. The HC is subjected to 1603 Pa and RC is subjected to 1587 Pa for this particular set of values. The negligible difference between the pressures exerted on the canopies is a reason which explains why both the canopies have drag coefficients.

   From (23) it can be seen that the velocity varies with respect to projected area in the form

   = (31)

   The effect of area was studied for both HC and RC and it was observed that it varies as per (30). For a particular paratrooper the effect of the canopy area was studied and the results are as follows.

   Paratroopers dimensions (Cylinder dimensions) Length (L) = 1.98 m

   Diameter (D) = 1.31 m Mass (m) = 121 kg Cd (Free Fall) = 0.64.

   Hemispherical Canopy

   TABLE X. Diameters and corresponding Terminal velocities for HC

   Diameter m	Area m2	vt m/s
   5.3048	22.1	8.611
   6.3048	31.22	7.197
   7.3048	41.91	6.131
   8.1048	51.59	5.464
   9.3048	68	4.701

   Also the time taken to reach the ground (tg)by the paratrooper from different height of jumps (HOJ) has been studied and it was observed that as the HOJ increases for a particular parachute it varies linearly.

   tG vs Diamter

   tG vs Diamter

   750

   650

   550

   450

   350

   250

   750

   650

   550

   450

   350

   250

   5

   6 7

   8 9

   10

   5

   6 7

   8 9

   10

   Diameter m

   tG s,HOJ=5000 m tG s, HOJ=6000 m tG s, HOJ=7000 m tG s, HOJ=8000 m tG s, HOJ=9000 m

   tG s, HOJ=10000 m

   Diameter m

   tG s,HOJ=5000 m tG s, HOJ=6000 m tG s, HOJ=7000 m tG s, HOJ=8000 m tG s, HOJ=9000 m

   tG s, HOJ=10000 m

   tG s

   tG s

   Fig.11 Comparison of tg with respect to area for different of HOJ using HC

   9

   8

   7

   6

   5

   4

   Terminal velocity vs Area

   vt = 45.952A-0.54 R² = 0.9998

   9

   8

   7

   6

   5

   4

   Terminal velocity vs Area

   vt = 45.952A-0.54 R² = 0.9998

   20

   60

   20

   60

   vt m/s Power (vt m/s)

   vt m/s Power (vt m/s)

   40

   Area m2

   40

   Area m2

   vt m/s

   vt m/s

   Fig.10 Variation of terminal velocity with Area of HC

   Using best fit method the variation of terminal velocity with respect to area is

   = 45.9520.54 (32)

   Comparing (31) and (32) it is observed that = 45.952

   and = 0.54.

   Ram-Air Canopy

   TABLE XI. Area and corresponding Terminal velocities for RC

   Span m	Chord m	Area m2	Vt m/s
   3.92	1.71	6.73	16.501
   4.45	1.94	8.63	14.569
   4.99	2.19	10.95	12.931
   5.49	2.77	15.18	11.174
   6.28	2.78	17.44	10.423

   Here the terminal velocity varies with area in the form

   = 42.7960.5 (33)

   These results also agrees with (31) and by comparison between (31) and (33) shows

   = 42.796 and = 0.5.

   Terminal velocity vs Area

   17

   16

   15

   14

   13

   12

   11

   10

   Terminal velocity vs Area

   17

   16

   15

   14

   13

   12

   11

   10

3. CONCLUSION

   v = 42.796A-0.5

   v = 42.796A-0.5

   t

   t

   R² = 1

   6

   R² = 1

   6

   8 10

   8 10

   12

   12

   14

   14

   16 18

   16 18

   Area m2

   Vt m/s Power (Vt m/s)

   Area m2

   Vt m/s Power (Vt m/s)

   Velocity m/s

   Velocity m/s

   CFD has been applied to this particular problem and the drag forces for short cylinders of different L/D, hemispherical and ram-air canopies have been found. The drag coefficients remained almost constant and it is known that it remains constant for most geometry which is subjected to flows having > 104. And the terminal velocity and time taken to reach the ground also agrees with the common understanding. As the area increases the drag increases and thus results in lower terminal velocity and as the terminal velocity reduces the time to reach the ground will increase. And it can be concluded that any of the two equation turbulence model can be used for this application.

   tg s

   tg s

   Fig.12 Variation of terminal velocity with Area of HC

   tG vs Area

   500

   450

   400

   350

   300

   250

   200

   150

   6 8 10 12 14 16 18

   Area m2

   tG s,HOJ=5000 m tG s, HOJ=6000 m tG s, HOJ=7000 m tG s, HOJ=8000 m tG s, HOJ=9000 m tG s, HOJ=10000 m

   tG vs Area

   500

   450

   400

   350

   300

   250

   200

   150

   6 8 10 12 14 16 18

   Area m2

   tG s,HOJ=5000 m tG s, HOJ=6000 m tG s, HOJ=7000 m tG s, HOJ=8000 m tG s, HOJ=9000 m tG s, HOJ=10000 m

   Fig.11 Comparison of tg with respect to area for different of HOJ using RC

4. REFERENCE

1. C. P. Kothandaraman, S. Subramanyan, Heat and Mass Transfer Data Book, 7th edition, New Age International Publishers, 2010.

2. Fluent 6.3 Users Guide, Fluent Inc, 2006

3. H. K. Versteeg, W. Malalasekera, An introduction to Computational Fluid Dynamics, 1st edition, Longman Scientific & Technical, 1995.

4. R. Vaidya, R. Bhalwar, S. Bobdey, Anthropometric Parameters of Armed Forces Personnel, MJAFI 2009; 65: 313-318

5. O. Rouaud, M. Havet, Computation of the airflow in a pilot scale clean room using turbulence models, International Journal of Refrigeration 25(2002) 351-361

6. P. A. Costa Rocha, H. H. Barobosa Rocha, F. O. Moura Carneiro,

   M. E. Vieira da Silva, A.Valente Bueno, SST (shear stress transport) turbulence model calibration: A case study on a small scale horizontal axis wind turbine, Journal on Energy 65 (2014) 412- 418

7. Un Yuon Jeong, Hyun-Moo Koh, Hae Sung Lee, Finite element formulation for the analysis of turbulent wind flow passing bluff structures using the RNG model, Journal of Wind Engineering and Industrial Aerodynamics 90 (2002) 151-169

8. Yunus A. Cengel, John M. Cimbala, Fluid Mechanics: Fundamentals and Applications, 1st edition, Tata McGraw Hill, 2006.