2.3 Vertical and Horizontal Axis Wind Rotors
3.2.2 RNG k- ε Model
The RNG-based k-ε turbulence model is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group” (RNG) methods. The analytical derivation results in a model with constants different from those in the standard k-ε model. The additional terms and functions in the transport equations for k and ε are also different.
Transport Equations for the RNG k-ε Model
The turbulence kinetic energy, k, and its rate of dissipation, ε, are obtained from the following transport equations:
24 The term of Gk represents the production of turbulence kinetic energy. Using the exact equation for the transport of k, this term may be defined as 𝐺� =
−𝜌𝑢��𝑢�� �����
�. The quantities αk and αε are the inverse effective Prandtl numbers for k and ε, respectively. Sk and Sε are user-defined source terms.
Modeling the Turbulent Viscosity
The scale elimination procedure in RNG theory produces a differential equation for turbulent viscosity:
Eq. (3-7) is integrated to obtain an accurate description of how the effective turbulent transport varies with the effective Reynolds number (or eddy scale), allowing the model to be better handled in the low-Reynolds-number and near-wall flows.
In the high-Reynolds-number limit, Eq. (3-7) gives
𝜇� = 𝜌𝐶���� (3-8) with Cμ= 0.0845, derived using RNG theory.
25
RNG Swirl Modification
Turbulence, in general, is affected by rotation or swirl in the mean flow. The RNG model in Fluent provides an option to account for the effects of swirl or rotation by modifying the turbulent viscosity appropriately. The modification takes the following functional form:
𝜇� = 𝜇�� 𝑓(𝛼�, 𝛺,��) (3-9) where 𝜇�� is the value of turbulent viscosity calculated without the swirl modification using either Eq. (3-7) or Eq. (3-8). Ω is a characteristic swirl number evaluated within Fluent, and αs is a swirl constant that assumes different values depending on whether the flow is swirl-dominated or only mildly swirling. This swirl modification always takes effect for axisymmetric, swirling flows and three-dimensional flows when the RNG model is selected. For mild swirling flows (the default in Fluent), αs is set to 0.07. For strong swirling flows, however, a higher value of αs can be used.
Calculating the Inverse Effective Prandtl Numbers
The inverse effective Prandtl numbers, k and ε, are computed using the following formula derived analytically by the RNG theory:
�����.����
The main difference between the RNG and standard k-ε models lies in the additional term in the ε equation given by
26
𝑅�����������(���/�� �)��� (3-11) where η ≡ Sk / ε, η0 = 4.38, β = 0.012.
The effects of this term in the RNG ε equation can be seen more clearly by rearranging Eq. (3-6). Using Eq. (3-11), the third and fourth terms on the right-hand side of Eq. (3-6) can be merged, and the resultant ε equation can be rewritten as becomes larger than C2ε. In the logarithmic layer, for instance, it can be shown that 𝜂 ≈ 3.0 gives 𝐶��∗ ≈ 2.0, which is close in magnitude to the value of 𝐶��∗ in the standard k-ε model (1.92). As a result, for weakly to moderately strained flows, the RNG model tends to give results largely comparable to the standard k-ε model.
In regions of large strain rate (η > η0), however, the R term makes a negative contribution, making the value of 𝐶��∗ less than C2ε. In comparison with the standard k-ε model, the smaller destruction of ε augments ε, reducing k and, eventually, the effective viscosity. As a result, in rapidly strained flows, the RNG model yields a lower turbulent viscosity than the standard k-ε model.
Thus, the RNG model is more responsive to the effects of rapid strain and streamlining curvature than the standard k-ε model, which explains the superior performance of the RNG model for certain classes of flows.
27
Model Constant
The model constants C1ε and C2ε in Eq. (3-6) have values derived analytically by the RNG theory. These values, used by default in Fluent, are
𝐶�� = 1.42; 𝐶�� = 1.68 3.2.3 Standard Wall Functions
The standard wall functions in Fluent are based on the proposal of Launder and Spalding (1974), and have been most widely used for industrial flows.
Momentum
The law-of-the-wall for mean velocity yields
𝑈∗ = ��𝑙𝑛(𝐸𝑦∗) (3-14) where
𝑈∗ ≡ ������ ��� ��� ��
�� (3-15) 𝑦∗ ≡���� �� ���� �� �� (3-16) In which
k = von Karman constant (= 0.487) E = empirical constant (= 9.793)
UP = mean velocity of the fluid at point P KP = turbulent kinetic energy at point P yP = distance from point P to the wall μ = dynamic viscosity of the fluid
28
3.3 Boundary Conditions
In the domain of interest mentioned above, boundary conditions are described at the rotating wind rotor, inlet surfaces, outlet surfaces, side surfaces (atmosphere), and walls (curtain).
1. Rotation boundary condition
According to Eq. (2-3), the angular speed ω (rad/s) of the wind rotor is expressed as
𝜔 = ����� (3-23) where D is the outer wind rotor diameter, ω the angular wind rotor speed, v1
the wind velocity, and λ the tip-speed ratio (TSR).
2. The inlet boundary condition The inlet boundary conditions are:
𝑢 = 𝑢��
𝑣 = 0 𝑤 = 0
where u, v and w represent the velocity components in X, Y and Z directions, respectively.
3. The outflow boundary condition
Outflow boundary conditions are applied at downstream flow exit, where the details of the local velocity and pressure are not known in advance. It is set by Fluent (2010) internally that the mass conservation is definitely maintained.
29
4. The symmetrical boundary condition
In the atmospheric case, the free surface boundary conditions, where the local velocity gradient approximate zero, are applied for side surfaces, provided that the distances are far enough from the center line of the domain. Via a series of numerical tests, the distance between the free surface and center line is chosen five times of the rotor diameter as Akwa et al. [21] do.
5. Wall boundary condition
The wall boundary conditions satisfy the no-slip condition that are u, v, w = 0.
3.4 Introduction to Fluent Software
Fluent is a state-of-the-art computer program for modeling fluid flow and heat transfer in complex geometries. It provides complete mesh flexibility, including the ability to solve the flow problems using unstructured meshes that can be generated about complex geometries with relative ease. Supported mesh types include 2-D triangular/quadrilateral, 3-D tetrahedral/hexahedral/pyramid, and mixed (hybrid) meshes. Fluent also allows refining or coarsening grid based on the flow solution.
Fluent is written in the C computer language and makes full use of the flexibility and power offered by the language. Consequently, true dynamic memory allocation, efficient data structures, and flexible solver control are all possible. In addition, Fluent uses a client/server architecture, which allows it to run as separate simultaneous processes on client desktop workstations and powerful computational servers. This architecture allows for efficient execution, interactive control, and complete flexibility between different types of machines
30
or operating systems.
All functions required to compute a solution and display the results are accessible in Fluent through an interactive, menu-driven interface.
3.5 Numerical Method
This study employs the computational fluid dynamics software Fluent to analyze the flow fields around rotating Savonius wind rotors. The finite volume iteration and SIMPLE algorithm are put in use to solve the governing equations of a transient flow field. And the corresponding grid movement is also solved by using sliding mesh method.
Fluent uses Segregated Solver method to solve the governing integral equations for the conservation of mass and momentum, and (when appropriate) for energy and other scalars such as turbulence and chemical species. In case a control-volume-based technique is used that consists of:
1. Division of domain into discrete control volumes using a computational grid.
2. Integration of the governing equations on the individual control volumes to construct algebraic equations for the discrete dependent variables such as velocities, pressure, temperature, and conserved scalars.
3. Linearization of the discretized equations and solutions of the resultant linear equation system yield updated values of the dependent variables.
3.5.1 Segregated Solution Method
Using this approach, the governing equations are solved sequentially (i.e., segregated from one another). Because the governing equations are non-linear (and coupled), several iterations of the solution loop must be performed before a
31
converged solution is obtained. Each time of iteration consists of the steps illustrated in Fig. 3.5 and outlined below:
1. Fluid properties are updated, based on the current solution. (If the calculation has just begun, the fluid properties will be updated based on the initialized solution.)
2. The u, v, and w momentum equations are each solved in turn using current values for pressure and face mass fluxes, in order to update the velocity field.
3. Since the velocities obtained in Step 2 may not satisfy the continuity equation locally, a Poisson-type equation for the pressure correction is derived from the continuity equation and the linearized momentum equations. This pressure correction equation is then solved to obtain the necessary corrections to the pressure and velocity fields and the face mass fluxes that continuity is satisfied.
4. Where appropriate equations for scalars such as turbulence, energy, species, and radiation are solved using the previously updated values of the other variables.
5. When inter-phase coupling is to be included, the source terms in the appropriate continuous phase equations may be updated with a discrete phase trajectory calculation.
6. A check for convergence of the equation set is made.
These steps are continued until the convergence criteria are met.
3.5.2 Linearization: Implicit
In the segregated solution method the discrete, non-linear governing equations are linearized to produce a system of equations for the dependent variables in every computational cell. The resultant linear system is then solved to
32
yield an updated flow-field solution.
The manner in which the governing equations are linearized may take an implicit form with respect to the dependent variable (or set of variables) of interest. The implicit form is described in the following:
Implicit
For a given variable, the unknown value in each cell is computed using a relation that includes both existing and unknown values from neighboring cells.
Therefore each unknown will appear in more than one equation in the system, and these equations must be solved simultaneously to give the unknown quantities.
In the segregated solution method each discrete governing equation is linearized implicitly with respect to that equation's dependent variable. This will result in a system of linear equations with one equation for each cell in the domain.
Because there is only one equation per cell, this is sometimes called a scalar system of equations. A point implicit (Gauss-Seidel) linear equation solver is used in conjunction with an algebraic multi-grid (AMG) method to solve the resultant scalar system of equations for the dependent variable in each cell. For example, the x-momentum equation is linearized to produce a system of equations in which u-velocity is the unknown. Simultaneous solution of this equation system (using the scalar AMG solver) yields an updated u-velocity field.
In summary, the segregated approach solves for a single variable field (e.g., p) by considering all cells at the same time. It then solves for the next variable field by again considering all cells at the same time, and so on. There is no explicit option for the segregated solver.
33
3.5.3 Discretization
Fluent uses a control-volume-based technique to convert the governing equations to algebraic equations that can be solved numerically. This control volume technique consists of integrating the governing equations about each control volume, yielding discrete equations that conserve each quantity on a control volume basis.
Discretization of the governing equations can be illustrated most easily by considering the steady-state conservation equation for transport of a scalar quantity 𝜙. This is demonstrated by the following equation written in integral form for an arbitrary control volume V as follows:
∮ 𝜌𝜙𝑣⃗ ∙ 𝑑𝐴⃗ = ∮ 𝛤�𝛻𝜙 ∙ 𝑑𝐴⃗ + ∮ 𝑆� �𝑑𝑉 (3-24)
Eq. (3-24) is applied to each control volume, or cell, in the computational domain. The two-dimension, triangular cell shown in Fig. 3.6 is an example of such a control volume. Discretization of Eq. (3-24) on a given cell yields
∑�������𝜌�𝑣����⃗𝜙� � ∙ 𝐴����⃗ =� ∑�������𝛤� (𝛻𝜙)�∙ 𝐴����⃗ + 𝑆� �𝑉 (3-25) where
𝑁����� = number of faces enclosing cell
34
𝜙� = value of 𝜙 convected through face f 𝜌�𝑣����⃗ ∙ 𝐴� ����⃗ = mass flux through the face �
𝐴����⃗ = area of face f �
(𝛻𝜙)� = magnitude of 𝛻𝜙 normal to face f V = cell volume
The equations solved by Fluent take the same general form as the one given above and apply readily to multi-dimension, unstructured meshes composed of arbitrary polyhedral.
By default, Fluent stores discrete values of the scalar 𝜙 at the cell center (c0 and c1 in Fig. 3.6). However, face values 𝜙� are required for the convection terms in Eq. (3-25) and must be interpolated from the cell center values. This is accomplished using an upwind scheme.
First-Order Upwind Scheme
When first-order accuracy is desired, quantities at cell faces are determined by assuming that the cell-center values of any field variable represent a cell-average value and hold throughout the entire cell; the face quantities are identical to the cell quantities. Thus when first-order upwind is selected, the face value 𝜙� is set equal to the cell-center value of 𝜙 in the upstream cell.
3.5.4 Simple Algorithm
The SIMPLE algorithm uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field.
If the momentum equation is solved with a guessed pressure field p*, the
35
resulting face flux 𝐽�∗, computed from 𝐽� = 𝐽�� + 𝑑�(𝑝�� − 𝑝��) (where pc0 and pc1 are the pressures within the two cells on either side of the face, and 𝐽�� contains the influence of velocities in these cell. The term df is a function of 𝑎�, the average of the momentum equation 𝑎� coefficients for the cells on either side of face f.)
𝐽�∗ = 𝐽�∗� + 𝑑�(𝑝��∗ − 𝑝��∗ ) (3-26) does not satisfy the continuity equation. Consequently, a correction 𝐽�� is added to the face flux 𝐽�∗ so that the corrected face flux, 𝐽�
𝐽� = 𝐽�∗ + 𝐽�� (3-27) satisfies the continuity equation. The SIMPLE algorithm postulates that 𝐽�� be written as
𝐽�� = 𝑑�(𝑝��� + 𝑝��� ) (3-28) where 𝑝� is the cell pressure correction.
The SIMPLE algorithm substitutes the flux correction equations, Eq. (3-27) and (3-28), into the discrete continuity equation (∑�������𝐽�𝐴� = 0) to obtain a discrete equation for the pressure correction 𝑝� in the cell:
𝑎�𝑝� = ∑ 𝑎�� ��𝑝��� + 𝑏 (3-29) where the source term b is the net flow rate into the cell:
𝑏 = ∑�������𝐽�∗𝐴� (3-30) The pressure-correction equation, Eq. (3-29), may be solved using the algebraic multigrid (AMG) method. Once a solution is obtained, the cell pressure
36
and the face flux are used correctly.
𝑝 = 𝑝∗+ 𝛼�𝑝� (3-31) 𝐽� = 𝐽�∗ + 𝑑�(𝑝��� − 𝑝��� ) (3-32) Here 𝛼� is the under-relaxation factor for pressure. The corrected face flux 𝐽� satisfies the discrete continuity equation identically during each time of iteration.
3.5.5 Sliding Mesh
The sliding mesh model allows adjacent grids to slide relative to one another.
In doing so, the grid faces do not need to be aligned on the grid interface. This setup requires a means of computing the flux across the two non-conformal interface zones of each grid interface.
To compute the interface flux, the intersection between the interface zones is determined at each new time step. The resulting intersection produces one interior zone (a zone with fluid cells on both sides) and one or more periodic zones. If the problem is not periodic, the intersection produces one interior zone and a pair of wall zones (which will be empty if the two interface zones completely intersect), as shown in Fig. 3.7. The resultant interior zone corresponds to where the two interface zones overlap; the resultant periodic zone corresponds to where they do not. The number of faces in these intersection zones will vary as the interface zones move relative to one another. Principally, fluxes across the grid interface are computed using the faces resulting from the intersection of the two interface zones (rather than from the interface zone faces themselves).
In the example shown in Fig. 3.8, the interface zones are composed of faces
37
A-B and B-C, and faces D-E and E-F. The intersection of these zones produces the faces a-d, d-b, b-e, etc. Faces produced in the region where the two cell zones overlap (d-b, b-e, and e-c) are grouped to form an interior zone, while the remaining faces (a-d and c-f) are paired up to form a periodic zone. To compute the flux across the interface into cell IV, for example, face D-E is ignored and faces d-b and b-e are used instead, bringing information into cell IV from cells I and III, respectively.
3.6 Computational Procedure of Simulation
The complete operating procedure by using Fluent package software is carried out through the following processes sequentially.
3.6.1 Model Geometry
Before Fluent calculations, it is necessary to build a model. This study used the pre-processor software Gambit to build the geometry of the model. Divide the geometry into finite volumes in order to generate grids conveniently. The details of geometry information can be referred to Section 3.1.
3.6.2 Grid Generation
After building the geometry, the model has to use the pre-processor Gambit to generate grids as shown in Fig. 3.9. This step defines the different grid sizes in different volumes. The smaller grid size for the small volume will increase the accuracy of the simulation, but it also produce larger grid number which cause calculation difficulty. To consider the appropriate grid size for grid generation is important. The grid generation usually reduced the calculation cost under
38 parameters need to be confirmed. Select the solver formulation and choose the basic equations (e.g., laminar, turbulent, inviscid, chemical species, or heat transfer models, etc.) to solve the problem. Identify additional models needed such as fans, heat exchangers, porous media, etc. Specify material properties and the boundary conditions. Adjust the solution control parameters. Give an initialized value for iterate the flow field model. Finish those foregoing steps, start to calculate a solution and examine the results. If necessary, refine the grid or consider revisions to the numerical or physical model.
3.6.4 Grid-independence Test
The grid-independence test should be taken in advance to have a trade-off between the guarantee of acceptable accuracy and an affordable calculation cost.
As described in Section 3.1, there are two types of rectangular calculation domains: a single Savonius wind rotor and the parallel matrix system with four Savonius wind rotors. The grid-independence tests in 3-D simulation of this study can be separated into z-axis direction and x-y plane direction in one single Savonius wind rotor with the boundary conditions of wind velocity 7 m/s and TSR 0.8.
39
In the case of one single Savonius wind rotor, grid-independence test of z-axis direction is carried out. Grid numbers of 60, 80, 100, 120 and 140 in z-axis are tested and the simulation results are shown in Fig. 3.10(a) and Table 3.3.
Because the changing rate of Cp from grid number 100 to 120 is small enough and remains almost the same value while the grid number increases, the grid number of 100 is selected.
Table 3.3 Grid-independence Tests in z-axis Grid Number (z-axis) Cp Changing Rate (%)
60 0.173
80 0.178 0.025%
100 0.181 0.015%
120 0.180 -0.005%
140 0.181 0.005%
The grid numbers of 4,750, 4,840, 5,024, 5,232, 5,452, 5,898 and 6,492 in x-y plane are tested in the case of one single Savonius wind rotor. The results are shown in Fig. 3.10 (b) and Table 3.4. Similar to the reason as described in the case in z-axis, grid number of 5,452 is chosen.
40
Table 3.4 Grid-independence Tests in x-y Plane
Grid Number (x-y Plane) Cp Changing Rate (%)
4,750 0.175
4,840 0.176 0.0011%
5,024 0.178 0.0011%
5,232 0.180 0.00096%
5,452 0.181 0.00045%
5,898 0.180 -0.00022%
6,492 0.181 0.00017%
Following the procedure for the grid-independence tests of one single Savonius wind rotor in z-axis and x-y plane, the grid number of the parallel matrix system with four Savonius wind rotors domain is set accordingly, and the set grid numbers of three domains are listed in Table 3.5.
Table 3.5 Grid Numbers of Two Domains
Grid Number (3-D)
One Single Rotor 657,400
Parallel Matrix System without Curtain 1,827,520 Parallel Matrix System with Curtain 2,980,600
41
Fig. 3.1 Schematics of Savonius wind rotor geometry in experimental study
Fig. 3.2 The domain of a single Savonius wind rotor
42
Fig. 3.3 The domain of four two-bladed Savonius wind rotors without curtain in parallel matrix system
Fig. 3.4 The domain of four two-bladed Savonius wind rotors with curtain in parallel matrix system
43
Fig. 3.5 Overview of the segregated solution method
Fig. 3.6 Control volume used to illustrate discretization of a scalar transport equation
44
Fig. 3.7 Zones created by non-periodic interface intersection
Fig. 3.8 Two-dimensional grid interface
45
Fig. 3.9 User interface of Gambit
40 60 80 100 120 140
0.172 0.174 0.176 0.178 0.180 0.182
Cp
Grid numbers in z-axis
(a)
46
6600 6400 6200 6000 5800 5600 5400 5200 5000 4800 0.175
0.176 0.177 0.178 0.179 0.180 0.181 0.182
Cp
Grid numbers in x-y plane
(b)
Fig. 3.10 Grid-independent test: (a) z-axis; (b) x-y plane
47
CHAPTER 4
EXPERIMENTAL APPARATUSES
4.1 Experiment Layout
The experimental layout is shown in Fig. 4.1. The system can be separated
The experimental layout is shown in Fig. 4.1. The system can be separated