Computational Fluid Dynamics Analysis Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Computational Fluid Dynamics (CFD) is categorized as a subdivision of fluid mechanics which deals with the computational solution of fluid flows. Computational Fluid Dynamics comprise the governing differential equations and applicability of these governing equations depends on the nature of the flow. These governing equations basically consist of Continuity Equation, Conservation of Momentum also known as Navier Stokes Equation and Conservation of Energy. Fluent deals with inviscid and viscous flow both but in the present study we are dealing with the viscous flow, so mainly Continuity Equation and Navier Stokes Equation will be the area of concern. Some of the prominent advantages of using CFD methods include the accuracy and reliability of the results and lower cost of application of CFD as compared to the expensive experimental methods. CFD uses computational software which offers a user friendly platform that enables user to simulate any flow with various sets of test conditions.

CFD works on a principle of discretization where a flow domain is discretized in very small units called cells. This unit cell structure is known as mesh or grid. Several

discretization schemes are available in Fluent and choices can be made on the basis of the needs of the end result. These cells are used for the analysis of the flow problem. Fluent gives the properties of the fluid at every single node. Spatial discretization schemes available in Fluent are Least Squares Cell based, Green Gauss Node and Cell based which is to be chosen according to the flow pattern. A pre-processing is required before proceeding to the post-processing. ICEM CFD has been used as a pre-processor and Fluent as a post-processor. In addition to this, the boundary condition of the problem is one of the decisive factors that plays vital role in determination of the accuracy of the simulation process therefore, it is the matter of great importance to select an appropriate boundary condition so as to achieve a desired result.

This chapter provides a detailed explanation of various parameters used to set up the simulation of a fluid flow around a two dimensional rotating Vertical Axis Wind Turbine using NACA 0018 airfoils. This chapter deals with the types of viscous model, boundary conditions, discretization schemes, time step calculations, reference values and solver used. Grid is generated using ICEM CFD which is considered to be more advanced meshing tool than GAMBIT and flow is analyzed using FLUENT. Numerical values of forces obtained from the simulation are validated with the research article by Guerri et al (2007). by plotting a graph between force and angle of rotation. Another validation is done by comparing a graph between Cp and λ obtained by CFD computations with experimental results by M. C. Claessens (2006).

2.2 Fundamental Equations

As it has already been mentioned above that Computational Fluid Dynamics comprises the governing differential equations and applicability of these governing equations depends on the nature of the flow. These equations have their mathematical representations which can be employed individually or in a group depending on the need of the desired output. Three basic principles which govern the characteristics of the flow of any fluid are conservation of mass, momentum and energy. In present case, we are dealing with the equation of continuity with the application of K-Omega models.

The Continuity Equation or conservation of mass is given by (White, 2005):-

= 0 (2.1)

Navier Stokes Equation for incompressible flow is given by (White, 2005):-


2.3 Coordinate System of NACA 4 digit series

NACA 4 digit airfoil family is defined by a series of numbers where each number has its own significance for instance in NACA 0018 the first digit represents maximum camber as percentage of chord and it is denoted by 'm'. Second digit refers to distance of maximum camber from the leading edge in tenth of percentage denoted by 'p'. Rest of the two digits designates to maximum thickness of airfoil as a percentage of chord denoted by 't'. The coordinates of NACA 4 digit airfoil family is given by following equation (Source: Wikipedia):-



c is the chord length of the airfoil

x is the position along chord from 0 to c

y is the half thickness at a given value of x

t is the maximum thickness as a fraction of chord

The position of the coordinates on the upper curve of the airfoil (Xu, Yu) and lower curve (XL, YL) of the airfoil is given by

Xu = XL = X

Similarly, the position of the coordinates on lower curve of the airfoil (XL, YL) is given by

Yu = +Y

YL = -Y

2.4 Geometry Creation

ICEM CFD has been used as a pre-processor for grid generation. Symmetrical airfoil NACA 0018 has been used for the present study as shown in figure 2.1. Rotor is equipped with three airfoils each with the chord length of 0.107 m. located at 1200 from one another. Diameter of the rotor is set as 2 m. and Far field is located at 168 chords from the center of the rotor. Formatted point data for airfoil come from different sources. For this case NACA ASCII 4 digit series has been used. Two curves with 100 node points on each of them is drawn using create/modify tool going from leading edge to trailing edge. All the curves and surfaces are assigned separately with different part names, this helps setting up the boundary conditions.

2.5 Grid Generation

Two separate zones are created, rotor being rotary and square far-field being stationary. The far-field mesh which is of hexahedral type is less dense as compared to hexahedral mesh in the rotary zone. 2D planar blocking is created around the 1200 section of rotor as depicted in figure 2.2. Edge-curves are associated, this helps get a nice fit of mesh around the edges. Of several types of grids available in ICEM CFD like H-grid, O-grid, C-grid and Y or quarter O-grid, current meshing uses quarter O-grid along with C-grid around the airfoil. O-grid allows a uniform orientation of mesh around the geometry and C-grid captures the geometry of the airfoil. Block is then split to capture the airfoil geometry. Blades of the vertical axis wind turbine are made rotating by making a pair of opposite nodes periodic. This is done by setting up base, axis and angle in global mesh parameters and then selecting periodic vertices using the edit block tool. Hexahedral meshing is used considering the fact that hexahedral mesh provides more uniform and smooth meshing over tetrahedral or quad. Edge meshing parameters are used to get a more uniform mesh distribution. Mesh around the airfoil needs to be dense enough so as to have a smooth gradient change in fluid. This 1200 section is then copied and rotated to get a complete 3600 VAWT as shown in figure 2.5. Fluent does not accept structured mesh pattern therefore, the structured hexahedral mesh is converted into unstructured mesh before exporting it to Fluent. Figure 2.6 shows the complete set up of VAWT mesh.

2.6 Turbulence model

The end result of the flow problem primarily depends on the Reynolds number. Working with low Reynolds number is comparatively complex as it requires more precision and accuracy to deal with. Computational Fluid Dynamics offers a gamut of flow models which can be used individually as per the requirement of the end result. Various turbulent modeling and simulation techniques like Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), Detached Eddy Simulation model (DES), Reynolds stress model (RSM), K-Epsilon model, K-Omega model, Spalart-Allmaras model are available and each one of them can be effectively used in particular areas of applications.

According to Fluent manual:- "It is an unfortunate fact that no single turbulence model is universally accepted as being superior for all classes of problems. The choice of turbulence model will depend on considerations such as the physics encompassed in the flow, the established practice for a specific class of problem, the level of accuracy required, the available computational resources, and the amount of time available for the simulation. To make the most appropriate choice of model for a particular set of application, one needs to understand the capabilities and limitations of the various options".

In the present study, wall bounded turbulent flows around the vertical axis wind turbine has been modeled using SST K-omega model. K-omega model is categorized as standard and shear stress transport (SST). SST model for K-omega differs from standard model in a context that SST provides a change in a gradual manner from standard K-omega model in the inner region, to high Reynolds number flow with K-epsilon model in outer region. Currently, our area of consideration is to determine the forces acting on each of the three rotating airfoils and to obtain an optimum value of tip speed ratio which gives the maximum power output when wind passes the turbine at a speed of 10 m/s. In the present study Reynolds number is set as 1.086*106 for a rotor diameter of 2 m.

2.7 Boundary Conditions

Boundary condition in Fluent defines the flow parameters at the boundaries of the flow domain. The end result depends on the boundary condition to a great extent. There are various boundary types available in Fluent like pressure, velocity and mass flow inlet, pressure outlet, pressure far-field AND outflow. Stationary and moving wall, axis are also the types of boundaries. In our current research, periodic boundary condition is applied to set the airfoils rotating.

In order to use velocity inlet as a boundary type the magnitude and direction of the velocity must be known. The possible pairs of boundary types at the inlet and exit are:-

Pressure inlet - Pressure outlet

Mass flow inlet - Pressure outlet

Velocity inlet - Pressure outlet or Outflow

For our study a boundary pair of velocity inlet and outflow has been used. Outflow boundary condition is generally suitable for the simulation of airfoil related problems. Airfoils are considered to be a moving wall in a moving fluid zone.

2.8 Problem set up in Fluent

The Rotating hexahedral mesh of the rotor is merged with the stationary mesh of the far-field. Mesh is then checked into ICEM CFD for all possible errors like unconnected vertices, periodicity, uncovered faces and single elements. ANSYS is used as a common structural solver and Fluent_V6 as an output solver to write a mesh file which could be read into Fluent. Two dimensional double precision with parallel solver is used. Mesh is checked again in Fluent for any negative volumes and skewness. Fluent also allows scaling the size of the working domain and at the same time user can set the units in SI, CGS and other format. In order to apply periodicity Turbo-outer and Far-inner both the parts are changed from wall boundary types to interface boundary types. Now mesh interface is created by selecting both the interface zones which allows us to set the rotational periodic boundary condition. Moving mesh technique is applied for this simulation where rotor is set to rotate at 380 RPM and the far-field remains stationary. Various types of pressure velocity coupling schemes are available in Fluent and their selection depends on various factors. The present study involves the application of PISO scheme. Among several special discretization schemes available in Fluent Green-gauss node based gradient with Presto pressure and second order upwind scheme is applied. Simulation begins with first order upwind scheme and then continues with second order after the first convergence is achieved, just to avoid instability in flow. A convergence criterion for the solution is set at 10-6.

2.9 Time step calculations

Unsteady simulation involves time dependent calculations. Time step is calculated using speed of the rotor Ω (rpm).

Rotational speed of the rotor Ω = 380 rpm = 380/60 revolution/sec

= 6.33 rps

Or, Rotor makes 6.33 revolutions in one second

Or, it takes 0.1579 seconds to make 1 revolution (or 360o)

Therefore, time taken to rotate 360o is 0.1579 seconds

Time step size is given as 1.052*10-4 seconds

Therefore number of time steps required for one revolution

= 0.1579/1.052*10-4 = 1500 steps

Maximum iterations/time step is assumed to be 80.

2.10 Reference values

Chord length = 0.107 m

Reference Length = Radius of the Rotor = 1 m.

Area = Rotor diameter* span = 2*1 = 2 m2 for 2D Span = 1

Depth = Span = 2.64 m

Enthalpy = 0 jule/kg

Pressure = 0 atm

Density = 1.225 kg/m3

Temperature = 288.16 k

Reynolds number = 1.086e06

Viscosity = 1.8421e-05 kg-s/m

Turbulent Kinetic energy = 1.5 m2/s2

Turbulent dissipation rate = 1386 s-1

2.11 Calculation of Torque imposed by horizontal and vertical forces acting on Airfoils

Forces acting on each airfoil are sum of the two components that is pressure force and viscous force. In order to calculate torque produced by each airfoil, this force is further split into horizontal and vertical direction. Component of these two forces do not always contribute to the torque due to rotary motion of the turbine. The component of force for which axis of rotation of turbine lies in its direction, produces no torque. Forces acting on each airfoil are set to be reported in FLUENT in horizontal and vertical direction. Torque produced by each airfoil is then calculated by following equations:-

T1 = - Fx1.Rcosθ - Fy1.R.sinθ

T2 = -Fx2.Rcos(θ+120) - Fy2.R.sin(θ+120)

T3 = -Fx3.Rcos(θ+240) - Fy3.R.sin(θ+240)

Total Torque T = T1 + T2 + T3 (Total Torque 'T' is calculated for every 6 degrees of rotation)

Average Torque Σ T/n   (0- 360 degrees) where n = 360/6 = 60

2.12 Grid independent study

Computational results obtained by CFD simulation must be grid independent. The results should not vary with the number of cells in mesh. Therefore grid independency is one of the important parameters to check the accuracy of the solution. Simulation is run for a cell size of 65000 and 140,000 and then component of the forces are plotted as a function of angle of attack. Results are found to be independent of the number of cells with negligible difference. Figure 2.7 and Figure 2.9 shows the grid independent solution obtained by simulations. These graphs also serve the purpose of validation as figure 2.7 and 2.9 follows the same trend as shown in figure 2.8 and 2.10 by Guerri et. al (2007).

2.13 Validation of coefficient of performance of VAWT obtained by 2D CFD simulation

Simulation is set to run at several tip speed ratios ranging from 1 to 6 at Reynolds number of 106 and then a graph is plotted between Cp and λ. Results are found to be in a good agreement with the experimental result by M. C. Claessens (2006) as shown in figure 2.11. Maximum Cp = 0.34 is obtained at λ= 3.8.

2.14 Von Karman Vortex Street

Figure 2.13(b) shows a beautiful phenomenon in fluid dynamics. When fluid flows over a cylindrical body it creates periodic pattern of swirling and bubbly flows which is called Von Karman Vortex Street. This wake shedding is the result of flow separation.

Figure 2.1 Schematic view of the geometry of rotor (120o) with NACA0018 airfoil.

Figure 2.2 Depiction of (a) Blocking scheme with the application of quarter O-grid (b) Periodic vertices to make a moving mesh

Figure 2.3 Schematic view of the Hexahedral meshing of rotor with NACA0018.

Figure 2.4 Closer view of the O-type grid around NACA0018 airfoil.

Figure 2.5 360o view of rotor (unstructured hexahedral mesh) with three airfoils.

Figure 2.6 schematic views of the stationary far-field and rotor (unstructured hexahedral mesh).

Figure 2.7: Grid Independent Result for cell size of 65000 and 140000 (Horizontal component of the blade force at λ =1.88, V =10 m/s)

Figure 2.7 The blade force component Fx at 180 rpm (Guerri et al. 2007).

Figure 2.9: Grid Independent Result for cell size of 65000 and 140000 (Vertical component of the blade force at λ =1.88, V =10 m/s)

Figure 2.8 The blade force component Fy at 180 rpm (Guerri et al. 2007).

Figure 2.11: Validation of Cp of VAWT as a function of λ with Experimental result

by M. C. Claessens (2006) at Re = 200000

Figure 2.12: Depiction of Leading and Trailing edge Vortex formation

(a) Wake leaving from LE and TE (b) Von Karman Vortex Street

Figure 2.13: Formation of Swirling vortices at λ = 2