Determination of sea water intake height caused by transient waves using CFD model
Abstract
Computational Fluid Dynamics is a mathematical tool, used to simulate fluid flow problems. The simulation results, are then used to consider kinematic and thermodynamics of flow particles within a specified geometry using Euler and Navier-Stoks solvers to simulate the most complex geometries in two and three dimensions. Sea water intake is one of complex geometry. In order to understand flow behavior incorporating control volume definition, it is essential to study all the forces (both internal and external).The transient analysis of the seawater intake in case of a power failure of Seawater pumps has been performed in this model. In this study, we have used Finite Volume Method to solve Navier-Stokes equations along with standard K-ε turbulent equations as well as Volume of Fluid (VOF) equations which are governing the free surface fluid motion.
Determination of sea water intake height caused by transient waves using CFD model
Mehdi Assadi Niazi
1- MS of water structures Engineering, Ardabil Regional Water Company.
Abstract
Computational Fluid Dynamics is a mathematical tool, used to simulate fluid flow problems. The simulation results, are then used to consider kinematic and thermodynamics of flow particles within a specified geometry using Euler and Navier-Stoks solvers to simulate the most complex geometries in two and three dimensions. Sea water intake is one of complex geometry. In order to understand flow behavior incorporating control volume definition, it is essential to study all the forces (both internal and external).The transient analysis of the seawater intake in case of a power failure of Seawater pumps has been performed in this model. In this study, we have used Finite Volume Method to solve Navier-Stokes equations along with standard K-ε turbulent equations as well as Volume of Fluid (VOF) equations which are governing the free surface fluid motion.
Keywords: Transient, CFD(computational fluid dynamics)Sea water,intake, Navier-Stokes flow.
1. INTRODUCTION
Computational Fluid Dynamics, which referred to more often as CFD, is a mathematical tool, used to simulate a fluid flow problems. The simulation results, are then used to consider kinematic and thermodynamics of flow particles within a specified geometry. This would give useful hints to optimize the existing design and relevant facilities. As it is obvious, all flow fields can be described by appropriate partial differential equations which are referred as government equations. In fact, CFD is an approach to solve these equations numerically and consequently obtain flow variables. The number of unknown variables e.g. velocity, pressure, temperature … are required to completely define a flow fields variables case by case. Nevertheless, it can be shown mathematically that a unique numerical solution exists for most cases even for more complicated flows such as turbulent and free surface motion flows. However, the number of complicated problems, which can be calculated and evaluated by CFD, has been increasing by developing more advanced numerical techniques. The main matter which limits the engineering application of CFD is difficulty in describing complex geometries. In recent years, many researchers have been conducted to overcome this shortcoming. Developing and using of unstructured grids is the most appropriate one. This gives more ability to Euler and Navier-Stoks (flow governing equations) solvers to simulate the most complex geometries in two and three dimensions.
generally, the performance of generating unstructured grids is so high and it does not require the user skills. Figure 1 shows the complex geometry of Bandar Abbas sea water intake facility. For such a complex geometry, the flow field can be claimed to be known whenever velocity and some static quantities are determined in every point of the model at every time. The number of static quantities, used for particular problems, depends on the nature of flow. Incompressible flows may be described by only one quantity such as pressure. However, for compressible flows one additional quantity e.g. density or temperature is required to define the state of fluid in every point of flow domain. This is not the case for the current project where the fluid (Sea Water) is definitely incompressible.
In order to understand flow behavior incorporating control volume definition, it is essential to study all the forces (both internal and external) that can be applied by flow on it. Conservative laws will describe the interaction between internal and external forces with fluid and flow characteristic. In fact, conservative governing equations, define that all flow properties such as energy and mass will remain constant under the action of internal and boundary forces.
The flux of these properties has two parts. First part is due to the bulk fluid motion, which is known as “convection” terms. The second part or diffusion term is due to statistical motion of molecules that exists even for a quiescent fluid. It has been well known that these conservative laws expressed by differential equations can fully describe incompressible flow of Newtonian fluids whether in laminar or turbulent regimes.
As we mentioned earlier, it is difficult to use structured grids for complex geometry of considered basin. Usually, this results to highly stretched grids, which affects the accuracy of solutions. Hence, unstructured grids are used here to have enough flexibility to generate appropriate grid cells.
Figure 1. 3D Seawater Intake Geometry
2. General Used Multiphase Model
The Volume of Fluid (VOF) model can model two or more immiscible fluids by solving a single set of momentum equations and tracking the volume fraction of each of the fluids throughout the domain. Typical applications include the prediction of jet breakup, the motion of large bubbles in a liquid, the motion of liquid after a dam break, and the steady or transient tracking of any liquid-gas interface such as Seawater intake.
The VOF formulation is generally used to compute a time-dependent solution, but for problems in which you are concerned only with a steady-state solution, it is possible to perform a steady-state calculation. A steady-state VOF calculation is sensible only when your solution is independent of the initial conditions and there are distinct in flow boundaries for the individual phases. For example, since the shape of the free surface inside a rotating cup depends on the initial level of the fluid, such a problem must be solved using the time-dependent formulation. On the other hand, the flow of water in a channel with a region of air on top and a separate air inlet can be solved with the steady-state formulation.
The VOF formulation relies on the fact that two or more fluids (or phases) are not interpenetrating. For each additional phase that you add to your model, a variable is introduced: the volume fraction of the phase in the computational cell. In each control volume, the volume fractions of all phases sum to unity. The fields for all variables and properties are shared by the phases and represent volume-averaged values, as long as the volume fraction of each of the phases is known at each location. Thus the variables and properties in any given cell are either purely representative of one of the phases, or representative of a mixture of the phases, depending upon the volume fraction values. In other words, if the qth fluid’s volume fraction in the cell is denoted as then the following three conditions are possible:
: The cell is empty (of the qth fluid)
: The cell is full (of the qth fluid)
: The cell contains the interface between the qth fluid and one or more other fluids
Based on the local value of 
, the appropriate properties and variables will be assigned to each control volume within the domain.
The tracking of the interface(s) between the phases is accomplished by the solution of a continuity equation for the volume fraction of one (or more) of the phases. For the qth phase, this equation has the following form:
(1)
where 
is the mass transfer from phase q to phase p and 
is the mass transfer from phase p to phase q. By default, the source term on the right-hand side of equation (1), 
, is zero, but you can specify a constant or user-defined mass source for each phase.
The volume fraction equation will not be solved for the primary phase; the primary-phase volume fraction will be computed based on the following constraint:
(2)
The properties appearing in the transport equations are determined by the presence of the component phases in each control volume. In a two-phase system, for example, if the phases are represented by the subscripts 1 and 2, and if the volume fraction of the second of these is being tracked, the density in each cell is given by
(3)
In general, for an n-phase system, the volume-fraction-averaged density takes on the following form:
(4)
All other properties (e.g., viscosity) are computed in this manner.
A single momentum equation is solved throughout the domain, and the resulting velocity field is shared among the phases. The momentum equation, shown below, is dependent on the volume fractions of all phases through the properties ρ and µ.
(5)
One limitation of the shared-fields approximation is that in cases where large velocity differences exist between the phases, the accuracy of the velocities computed near the interface can be adversely affected.
3 . Interpolation Near the Interface
Control-volume formulation requires that convection and diffusion fluxes through the control volume faces be computed and balanced with source terms within the control volume itself. There are four schemes for the calculation of face fluxes for the VOF model: geometric reconstruction, donor-acceptor, Euler explicit, and implicit.
In the geometric reconstruction and donor-acceptor schemes, We applies a special interpolation treatment to the cells that lie near the interface between two phases. Figure (2) shows an actual interface shape along with the interfaces assumed during computation by these two methods.
Figure 2. Interface calculations
The Euler explicit scheme and the implicit scheme treat these cells with the same interpolation as the cells that are completely filled with one phase or the other (i.e., using the standard upwind, second-order, QUICK, or modified HRIC scheme), rather than applying a special treatment. For more information, you can refer to CFD text books.
4. The Euler Explicit Scheme
In the Euler explicit approach, standard finite-difference interpolation schemes are applied to the volume fraction values that were computed at the previous time step.
(6)
Where n+1= index for new (current) time step
n = index for previous time step
= face value of the qth volume fraction, computed from the first or second-order upwind, QUICK, or modified HRIC scheme
V= volume of cell
= volume flux through the face, based on normal velocity
This formulation does not require iterative solution of the transport equation during each time step, as is needed for the implicit scheme.
5. The Implicit Scheme
In the implicit interpolation method, standard finite-difference interpolation schemes, including the modified HRIC scheme, are used to obtain the face fluxes for all cells, including those near the interface.
(7)
Since this equation requires the volume fraction values at the current time step (rather than at the previous step, as for the Euler explicit scheme), a standard scalar transport equation is solved iteratively for each of the secondary-phase volume fractions at each time step. The implicit scheme can be used for both time-dependent and steady-state calculations.
6. Solution method
In order to solve the conservation equations, a Finite-Volume Method is used to discredits them on each CV (control volume). For each CV, one algebraic equation is obtained; each of the equations involves the unknown from the CV-centre and from all neighboring CVs with which the current CV has common faces. Since the equations are non-linear, they have to be linearised in order to solve by an iterative solution method. The equations are also coupled but are solved in a segregated manner, i.e. for each variable in turn, whereby other variables are treated as known by using the best update available from previous iteration. In the course of obtaining an algebraic equation for each CV, three levels of approximation are applied:
- The integrals over surface, volume, and time are evaluated using midpoint-rule approximations, which use the value of the integrand at the centre of the integration domain.
- Since the variable values are calculated at CV centers only, values at cell-faces centers required for the evaluation of integrals have to be obtained by interpolation; here the simplest second-order approximation, namely linear distribution, is chosen with upstream value.
- In order to evaluate stress terms (forces) at CV-faces, numerical differentiation is needed to compute derivatives of the Cartesian velocity components with respect to Cartesian coordinates. This is done using Gauss method. The time derivative at the current time level is computed by differentiating a parabola fit through values at the current time level.
All of the above mentioned approximations are nominally of second order.
The mass conservation equation is transformed into a pressure-correction equation following the
Well-known SIMPLE-algorithm for collocated arrangements of variables.
In this case of turbulent flow, the standard k&ε with wall function is used. The additional transport equations introduced by the model for the turbulent kinetic energy k and its dissipation rate ε have a form similar to that of the momentum equations and can be discredited and solved using the same technique. With k and ε the eddy viscosity 
can be determined.
Computation of three-dimensional flows with free surfaces, especially when they are unsteady, requires a lot of both memory and computing time.
7. Case study characteristics
The intake basin is a 84 m long concrete building consisting of a screening and a pump station with forebays. The pumping station will receive flow from three 1500 m long inlet culverts. The inlet culverts discharge directly to a optimized tapered inlet forebay feeding six drum screens. Downstream of the drum screens, a tapered forebay is located conveying the flow to the pump sumps. The general arrangement of the intake basin are:
- Inlet basin length: 34.4m
- Filters Diameter : 14.5
- Outlet basin length: 28.6 m
- Pump inlet house length: 6.5 m
In order to eliminate the vortices, diffuser blocks and a vortex suppression boom and concrete triangle blooks are planned to be installed in the pump forebay in attempt to disrupt the low level jets leaving the drum screens reducing the mass rotation of flow. Various diffuser locations were investigated and the best arrangement is proposed as described in the previous reports.
To evaluate transient behavior, free surface flow field at the maximum level situation were considered. In this case, the Sea water level is about z=13.4 m.
Please note that z=0 is equal to Elev.= -8.50 m CD in this report.
8. Transient conditions for pump trip.
For the transient calculation of the flow, it is assumed that all operating pumps (without standby pumps) are working under maximum flow conditions when the power failure occurs. After power failure, the flow rate will decrease to zero within a time period, ∆t of less than t = 2.0 s. The effect of this time discrepancy is negligible considering a sloshing period of 820 seconds. The time ∆t = 2.0 s was therefore taken as the stopping time for all pumps. Please note, as the pumps are axial pumps, in reality the flow will not completely decrease to zero within such a short time. Therefore the above made assumption represents the most unfavorable case. The time depending flow rate is shown in the Figure 3.
Figure 3. Flow rate as a function of time
9. Boundary and Initial Conditions
There are three types of boundaries in fluids:
- ‘Natural’ Boundaries, consisting of impermeable walls
- Inlet and outlet Flow boundaries
- At the walls, the no-slip boundary condition applies, i.e. the velocity of the fluid is equal zero.
The inlet boundary conditions can be specified at that portion of the boundary where the fluid enters the solution domain from sea where the hydraulic pressure distribution or fixed pressure on the surface of sea water is known. The inlet-type boundary conditions can be specified even where the fluid leaves the computational domain. This case was also applied. The intake to the pumps is treated as an inlet of the solution domain. The velocities were adjusted in such a way that the flow rates agree with the given flow rates. At the inlet of the culverts, where the water usually comes into the solution domain a pressure boundary condition is defined. The pressure is set equal to the hydrostatic pressure. Steady state flow condition in optimum geometry is set for initial condition of transient analysis. In a first numerical simulation all six drum screens had been modeled. In this case the number of cells was about 3 millions. Fig. (4) shows a 3D and 2D horizontal cut at z = 1 m through the calculation domain respectively. The flow in the drums is extremely complicated therefore the calculation time for each time step rises quickly. In order to obtain sufficient results within an acceptable time frame, the drum region was simplified and the six drum screens were replaced by six culverts. Porous media baffles were introduced in order to simulate in these culverts the same pressure loss as it occurs in the drum screens perpendicular to the flow.
|
Figure 4. horizontal cut at z = 1 m through the calculation domain-3D calculation domain
10. Transient Mathematical model
In this study, we have used Finite Volume Method to solve Navier-Stokes equations along with standard K-ε turbulent equations as well as Volume of Fluid (VOF) equations which are governing the free surface fluid motion. The flow is fully turbulent and there are two different situations where turbulence is needed to be modeled. The first situation is where the turbulence structure is to be studied. The second is when the effect of turbulence on a system is to be modeled. While studying the structure of turbulent flows, it is important to capture (simulate) all of the scales present or model them correctly. DNS or LES handle the simulation option while RANS like models (such as K-ε) handle the modeling option. In modeling, the model constants should be tuned in a case by case manner to ensure that the model behaves correctly in abstracting all of the scales present. Since the study of the structure of turbulence is not the concern of this work, none of simulation or fine tuned models are necessary. While studying the impact of turbulence on a system, if turbulence fine structure affects the system behavior then simulation or fine tuned modeling should be considered. Otherwise (which is the case in this work) a model based on a collection of modeling and experimental data is the only thing necessary. Such models like Standard k-e or k-w or their variants are long available in the literature. In this work the Realizable k-e model (T.-H. Shih, W. W. Liou, A. Shabbir, Z. Yang, and J. Zhu.A New k-epsilon Eddy-Viscosity Model for High Reynolds Number Turbulent Flows Models Development and Validation. It should be noted again that these constants are found by model developers based on the best fit to experimental or direct numerical simulation data and their use is a standard industrial practice. The only time it is necessary to tune the model constants is when the detail of the modeling (not its overall impact) is of importance. In this work turbulence affects the analysis via its diffusion and dispersion behaviors and not via its fine detailed structure. It should be emphasized that whatever method used to compute the effects of turbulence as long as the structure of it is not of concern will not affect the system modeling accuracy. Pressure coupling is done by SIMPLE algorithm. As discussed before, the flow is fully turbulent and standard k-epsilon model is used to evaluate turbulence variables. For outlet boundaries in the pump suction sections, constant mass flow rate are applied as boundary condition. At the free surface boundary, the outlet vent condition is applied. The performance of the basin will be evaluated for standard operation of pump station. The Fluent software was selected for these study on basis of its effectiveness in computation of such phenomena contain turbulent flow, swirling free surface flow, severe eddying and complex flow distribution.
11. TIME STEP
After the power failure occurs, the flow rate decreases within a period of t to zero. In the beginning of the transient simulation the time step was set to 0.006 s. But, after some iterations, the time step have been evaluated dynamically according Courant number.
12. WATER LEVEL
After the pumps stop due to a power failure, the water level starts to rise slowly. After about 817 s, the averaged water level reaches its highest value in the basin top surface of about 14.8 m. Fig. (5) illustrates the averaged water level as a function of time. The construction height of the screen and the pump forebay is 15.2 m (Elev. 6.70 CD) which is about 0.4 m above the highest water level. In the attached movie, you will find the time dependent rate of free surface level increasing. In Figures (6) and (7) water level has been demonstrated in the time of zero and after 820 sec.
Figures 5. Water level as a function of time
Figure 6. Water level after 820 s Figure 7. Water level in the time of zero
13. FLOW CONDITIONS
As the velocity of the flow through the screen drums is relatively small, the water level in the pump forebay nearly follows the water level of the screen forebay and reaches its maximum value with a small time shift. The water rises with a vertical velocity of 0.006 m/s and the water level remains plane. Figure (8) shows the inlet flow rate and volume of fluid in the water surface.
|
Figure 8. inlet flow rate- of fluid in the water surface. (Z=0 is equal to Elev. – 8.50 CD).
14. Conclusion
A transient analysis of the flow behavior in the sea water intake in case of a power failure of the SW pumps has been performed in order to determine the increase of the water level. The transient analysis was performed assuming the high Sea level condition (Z=13.4/Elev.= +4.90 CD) and a new culvert, as these conditions represent the most conservative case. After the pumps stop due to a power failure the water level in the two forebays will starts to rise slowly with a vertical velocity of 0.006 m/s. After about 817 s the water level reaches its highest value in the basin of about Z=14.8m (Elev. = +6.30 CD ). The maximum resulting water level is not critical compared to the height of the top covering the basin. The construction height is almost sufficient for all culvert conditions. However, client can decide to increase the height of basin.
15. References
1. Bergant,A & Tijsseling.A; (2001) ” Parameters Affecting Water Hammer Wave Attenuation, Shape and Timing”,IAHR 10th, , Trondheim, Norway, Calibration and Validation work.
2. H.K. Versteeg, W. Malalasekera,(2007), “An introduction to computational fluid dynamics, the finite volume methods” Pearson Educationv,USA.
3. D. L. Youngs,(2004),” Time-Dependent Multi-Material Flow with Large Fluid Distortion.
4.In K. W. Morton and M. J. Baines, editors, “(1982)Numerical Methods for Fluid Dynamics. Academic Press.
5. Streeter, V.L., Wylie, E.B.(1979),” Fluid Mechanics, MacGraw-Hill, New York.
6.Streeter, V.L and Benjamin. (1960), “Hydraulic Transient” McGraw-Hill company, New York.
گروه: فنی و مهندسی کلمات کلیدی:
