Finite element solution of the Navier—Stokes equations

ON THE RESOLUTION OF THE NAVIER-STOKES EQUATIONS BY THE FINITE ELEMENT METHOD USING A SUPG STABILIZATION TECHNIQUE. Application to some wastewater treatment problems *P. Vellando, **J. Puertas, *I. Colominas ETS. de Ingenieros de Caminos, Canales y Puertos. *Dpto. de Métodos Matemáticos y de Representación. **Dpto. de Tecnología de la Construcción. Universidad de La Coruña. Campus de Elviña. 15071 La Coruña, Spain. Tel: (34) 981 167000, fax: (34) 981 167170. e-mail: vellando@iccp.udc.es ABSTRACT In this paper an analysis of the incompressible flow has been carried out, from the very definition of the governing equations, up to the resolution of some practical problems, passing through the comprehensive study of the stabilized finite element techniques used in their resolution. As a consequence of this analysis, a code based upon a realistic interpretation of the forces has been written, which allows for the modelling of the open channel flow, with optimum results in the resolution of some benchmark and real flow problems related with the wastewater industry. Key words: Viscous incompressible flow, FEM, Navier Stokes, SUPG. INTRODUCTION The Finite Element Method was first developed in the fifties by Turner and Clough so as to solve some structural problems of the aeronautical industry [Turner 56]. The application of the Finite Element Method to the flow problems requires some modifications with respect to the formulation used for the structural stress analysis problems that were its first application. Some of these modifications have been borrowed from the finite difference or finite volume approaches, and many others have been specifically developed for finite elements. In the early seventies we find many works regarding not only the mere existence and consistency of these flow problems [Ladyzhenskaya 69], [Babuska 71], [Brezzi 74], but also many works that give a finite element solution to the Navier-Stokes equations [Baker 71], [Oden 72], [Fortin 72], [Taylor 73], [Shen 76], [Zienkiewicz 76]. Since then, the Finite Element Method is a powerful tool for the resolution of the Navier Stokes equations, which will be used in this work to solve the incompressible flow. The material we are going to deal with, when solving the flow, is of a fluid nature, and therefore it has not a fixed shape, which is instead a function of time. In addition to Newton’s second law, that rules any dynamic problem, an equation that ensures for the conservation of mass should be verified. Moreover, the Navier-Stokes equations are a set of differential equations with respect to both space and time in which both the pressure and the velocity are the unknowns. As a consequence, the finite element formulation used for the conventional structural analysis cannot be applied straightforwardly. When applying the finite element analysis to the problems of the rigid body, the weighted residual method can be exclusively applied to the Newton second law, which 1 for statics clearly turns out to be the equilibrium equation. On the contrary, when dealing with fluids, the shape is not any more conserved, and apart from stating the equilibrium of momentum, we have to ensure for the continuity of mass. Consequently, we have two equations to be verified at the same time, and the finite element formulation should also account for the verification of both. The only set of unknowns in the conventional structural analysis is that of the displacements, as a consequence, the system obtained thanks to the application of the Finite Element Method, gives the displacements in the structure depending on the stiffness matrix (that features the structure), and the load vector. In the flow problems, we are headed towards the socalled mixed Finite Element Methods, in which both the velocity and pressure set of unknowns have to be treated simultaneously. Depending on how these two sets of equations and unknowns are tackled, several different approaches have been developed. Some authors [Kim 88], [Choi 94], agree to classify the different approaches to solve the viscous incompressible flow by the Finite Element Method, into three different categories, these are the mixed (or velocitypressure integrated), the penalty and the segregated methods. The mixed method is the most intuitive of these approaches, and is based upon carrying out a similar analysis for the continuity equation to that used for the momentum equation, keeping both velocity and pressure as the unknowns up to the end of the problem, [Baker 71] [Oden 72], [Zienkiewicz 76]. This apparently straightforward way of dealing with our equations is not as simple as it appears to be, and it may be the reason of the obtaining of a meaningless solution when used in connection with a faulty basic element [Babuska 71], [Taylor 73], [Brezzi 74]. Besides a big expense in the storing memory, the socalled mixed formulation, leads to some consistency problems in the obtaining of the solution when a wrong choice in the basic functions has been made. As a consequence, the penalty and segregated methods have been developed to overcome these difficulties. In this work, an algorithm of each type will be developed, programmed, verified, employed and discussed. The 2D Navier-Stokes equations assume a flow that takes place on a two-dimensional plane, and it is therefore laminar in that sense. A fully 3D resolution of the flow, apart from leading us to high rate memory requirements, would involve some problems in the modelling of the free water surface. The Shallow Water formulation has been considered as a way of including the third dimension in the calculations, being able to give a meaningful solution for flows in which the depth is small compared to the horizontal dimension. The integration in depth of the 3D Navier-Stokes formulation causes the dependence of the continuity equation with respect to depth, and consequently the appearance of some quasi-non-linear terms that depend on both the velocity and the depth. These equations are solved thanks to a newly developed iterative algorithm, which will be solved on a mixed formulation basis to be described later in the text. The use of a Galerkin formulation, that takes weighting functions equal to trial functions when solving the Navier-Stokes equations, may lead to some problems of instability in the flow solution by the Finite Element Method. To avoid this difficulty, some so-called stabilization procedures have been released since the MAFELAP conference in 1975 [Zienkiewicz 76]. The stiffness matrix resulting from structural problems solved by the Finite Element Method is symmetric, instead the ‘stiffness’ matrix obtained for fluids is non-symmetric and the use of symmetric weighting functions may lead to some instability problems. The faster the flow turns, the more non-symmetric the coefficient matrix becomes. In practice this is featured by the appearance of some spurious node-to-node oscillations also known as ‘wiggles’. One 2 way of avoiding these oscillations is to carry out a refinement in the mesh, such that convection no longer dominates on an element level, but this refinement turns to be a memory resources sink. This point will be avoided in this work by the use of a stabilization technique of the SUPG type, for all the algorithms considered in it. The SUPG (Streamline/Upwinding Petrov-Galerkin) technique, first developed by Brookes [Brookes 82], succeeds in eliminating the spurious velocity field, without carrying out a severe refinement in the mesh, by considering weighting functions that differ from trial functions in an upwinding term. This method was first released for the transport equation, and its generalisation to the Navier-Stokes equations brings an additional problem; that is the appearance of an excessive diffusion normal to the flow. The SUPG method eliminates this spurious crosswind diffusion by considering an artificial diffusion that acts only in the direction of the flow. All the particulars regarded in this introduction and some others, will be further considered in the following chapters. GOVERNING EQUATIONS As in any other dynamic problem, the equation we are going to refer to, is the Newton second law, which gives the variation in the momentum as the summation of the acting forces on the volume of integration. To this condition we should add another one, due to the fact that we are dealing with a shape-changing matter in which we have to ensure the conservation of mass. Both equations make up the Navier-Stokes equations. These relations are named after their discoverer, the French civil engineer Claude-Louis Navier (1785-1836), who in 1821 formulated the equations that rule the incompressible flow. The Navier-Stokes equations also bear the name of the Irish mathematician George Gabriel Stokes (1819-1903), who not knowing the previous discoveries made by Navier, Poisson and Saint-Venant, re-obtained the Navier-Stokes equations for slightly different assumptions, and published these works in 1845. The Irish mathematician gives also his name to the simplified version of the Navier-Stokes equations, in which the convective terms are dropped. Using the indicial notation in which the indices after commas will stand for derivatives with respect to the variables specified in the index, we can express the Navier-Stokes equations as: 1 ui ,t + u j ui , j = − p,i + νui , jj + f i ρ in ui ,i = 0 together with the boundary and initial conditions: σ ij n j ui ]Γ1 = bi ui (x j ,0 ) = ui 0 (x j ) ] Γ2 Ω (1) = ti (2) where u i is the velocity, p is the pressure, f i is the body force per unit mass, ρ is density, ν is the cinematic viscosity, Γ1 and Γ2 are two non overlapping subsets of the piecewise smooth domain boundary Γ , b i is the velocity vector prescribed in Γ1 , t i is the traction vector prescribed on Γ2 , σ ij is the stress along the boundary Γ2 , and n j is the outward unit vector normal to Γ2 . The 2D Navier-Stokes equations will be used in this thesis to solve many benchmark problems of the related literature with very good results, as will became clear later in the 3 text. The 2D or laminar (in the sense of planar) Navier-Stokes equations do not take into account the third dimension in space, and provide with the velocities and pressures of a theoretical planar flow. Nevertheless, for many real flow problems, the third dimension in space is very important and the 3D Navier-Stokes equations should be considered. The three-dimensional Navier-Stokes equations result in a very large-dimensioned system of equations that involves very high computational costs. Moreover the 3D schemes present a great difficulty in the treatment of the free surface. For flows in which the horizontal dimension is small compared to depth, the Shallow Water formulation can be employed as a simplification of the 3D Navier-Stokes equations [Weiyan 92]. The Shallow Water equations are a simplification of the Navier-Stokes equations, which can be used when the main direction of the flow is the horizontal one and the distribution of the horizontal velocity along the vertical direction can be assumed as uniform. These equations assume that the vertical acceleration of the fluid is negligible and that a hydrostatic distribution of the pressure can be adopted. The Shallow Water equations are obtained by integrating the 3D Navier-Stokes equations in depth, and give a meaningful solution for flows in which the horizontal dimension is small compared with the depth. When a 2D Navier-Stokes equation is used, no attention is paid to the third dimension in space, and the results are based upon a 2D approach to the flow problem. Therefore, the continuity equation is only held on a 2D basis. So as to get some information about the variations in depth along the flow, either a 3D NavierStokes equation or the Shallow Water equations (if the flow can be regarded as shallow), should be used. The Shallow Water equations are solved in this work for that purpose. Integrating in depth the 3D Navier-Stokes equations, it is obtained: h,t + hui ,i + ui h,i = 0 (3) ui ,t + u j ui , j = − gh,i + νui , jj + g (S 0i − S fi ) where h is the depth, g is the gravity force and S 0i y S fi are the geometric and friction slopes. Let us make some comments on the evaluation of the friction slope. All the flows found in civil engineering practice can be featured by the Reynolds number (UL/ν, where U and L are the characteristic velocity and length of the flow and ν is the kinematic viscosity that depends on the fluid nature). For small Reynolds numbers, the flow can be regarded as laminar, and the streamlines are parallel to each other. As the Reynolds number is increased, a chaotic, random and intrinsically unsteady type of motion appears. If these turbulent effects are to be solved by using the Navier-Stokes equations, a very refined mesh would be required to capture the eddies taking place on a wide range of length scales, and a special attention should be devoted to the unsteady resolution of the turbulent phenomena, that take place at a very high frequency [Versteeg 95]. The mesh refinement and the time step required for this purpose are not yet computationally affordable and a turbulence model should be implemented in order to evaluate these turbulent eddies. Most of these turbulence models are based upon decomposing the involved variables into a mean value (within a time increment) and a fluctuating term that depends on time. As a consequence of this approach, a term that evaluates the turbulent losses as a function of a so-called eddy viscosity ν t , is obtained. To evaluate this eddy or turbulent viscosity, a specific turbulence model such as the k-ε model should be introduced. Making use of these turbulence models, the turbulent viscosity is calculated for each time step and position, allowing for the capturing of these eddies [Rodi 93]. Some other flow models evaluate this eddy viscosity as a constant within the flow domain, such as the RMA2 flow model 4 developed by the Brigham University, which is one of the most commonly used programs to evaluate the flow in channels. Another approach to the turbulent problems would be to evaluate the friction slope on a Manning basis. The integration in depth of the 3D Navier-Stokes equations allows for the empirical evaluation of the energy losses taking place in flows that can be regarded as shallow. The Manning formula evaluates empirically the overall energy losses taking place in the fluid flow, including those related with the turbulent effects. This formulation does not capture the turbulent eddies taking place within the fluid flow but takes into account the overall turbulent energy losses, moreover keeps the full NavierStokes formulation, being ready to incorporate a k-ε model. By doing so, the frictional slope can be taken as: n 2 ui u 2j (4) S fi = Rh4 / 3 where R h is the hydraulic radius and n is the Manning coefficient. FINITE ELEMENT FORMULATION Mixed laminar formulation This formulation is directly obtained from the application of the weighted residuals method on both the dynamic and continuity equations. In the Galerkin method, the weighting functions are taken equal to the shape functions, in terms of which the unknowns are interpolated. We are going to introduce an SUPG stabilization of the algorithm by considering an additive term pi in the weighting functions. Applying the weighted residuals method in this way on both the dynamic and continuity NavierStokes equations, we obtained: 1 h h h h h h h h h h ∫ wi u i,t + u j u i, j − f i dΩ + ν ∫ wi, j u i, j dΩ − ∫ wi,i pdΩ − ∫ t i wi dΓ2 + Ωh ( ) ρ Ωh Ωh   1 h + ∑ ∫ pih  uih,t + u hj ui,jh −νui,jj + p ,ih − f i h  dΩ = 0 ; ρ e Ωe   ∫q u h Ωh h i ,i Γ2 dΩ = 0 (5) where wi = wi + pi are the SUPG weighting functions for the dynamic equation and q are the weighting functions for the continuity equation. The h superscript stands for the discretization being carried out in our formulation, which will be made in terms of a Q1P0 basic element: uih = ∑ uij N i j ph = ∑ p j χ j 1 4 j =1 j =1 (6) This formulation should verify in addition to the governing equations, some consistency conditions, the most restrictive of which is the LBB (Ladyzhenskaya-Babuska-Brezzi), or divergence-stability condition [Babuska 71], [Brezzi 74]. A wrong election in this basic element may fail to verify this condition and prevent the whole algorithm from converging. The Q1P0 basic element, even not being strictly div-stability stable, has shown not to produce the well known checkerboard pressure mode for the flow 5 problems to be considered in this work as will be shown later, providing with stable and efficient solutions. Once the elementary matrices are evaluated and assembled, the integral equation (5) can be expressed in matrix notation as: ∂υ Mυ + Cυ (u, v )υ + νAυ υ − Bp = f ∂t (7) BT υ = 0 where M υ is the mass matrix, Cυ (u, v ) is the convective matrix, A υ is the viscous matrix, B is the pressure matrix, f is the external forces vector, p is the pressure vector, u is the velocity vector in the x direction, v is the velocity vector in the y direction and υ is the velocity vector. For the derivatives with respect to time included in the first term of the first equation in (7), a finite difference approach will be used in order to transform our partial differential system into an algebraic one. This term can be evaluated on a backward differencing scheme as:  υ n − υ n −1  ∂υ  (8) Mυ ≈ Mυ    t ∆ ∂t   is the unknowns where υ is the unknowns vector at the present iteration and υ vector obtained in the former iteration. Therefore, the second term in expression (8) can be brought to the right hand side of the equality. The convective term Cυ (u, v )υ that appears in this formulation, is not the product of a coefficient matrix times a vector of unknowns, but a non-linear velocity-dependent function. This term should be eliminated in order to transform the resulting system into a linear system of equations. The numerical scheme to be used for this transformation is going to be a so-called successive approximation method, that can be mathematically expressed as: n −1 n (9) ∫ wi u j ui, j dΩ ≈ ∫ wi u j ui, j dΩ n −1 n Ωh Ωh and has shown to provide with good results in the in the resolution of the Navier-Stokes equations [Gartling 74]. This method converges linearly to the solution in opposition to some others, like the Newton-Raphson method, which do it in a quadratic way in the vicinity of the solution. Still the necessity of an appropriate initial guess in the NewtonRaphson method may prevent the solutions from converging [Jamet 73], and a continuation technique, is often required. The mixed formulation so-implemented is still quite expensive in terms of storing memory requirements, with the associated coefficient matrix of the resulting system being 2M+N dimensional, where M and N are the number of the velocity and pressure unknowns respectively. So as to avoid this problem we are going to introduce the penalty and the segregated formulations. Penalty laminar formulation The main difficulty found when obtaining a numerical solution for the Navier-Stokes equations is that apart from verifying the dynamic constitutive equation, the solutions must satisfy in addition the incompressibility condition. The mixed finite element formulation leads to a system in which both velocity and pressure are taken as unknowns. Nonetheless, besides the problems entailed in the election of the basic elements in order to allow for the div-stability condition to be held, mixed methods 6 result in a large-dimensioned system. Therefore, not only a larger dimension has to be handled with its corresponding increased memory requirements, but also the stiffness matrix is found to be radically different to the narrow-band type of matrix, which is preferred for the direct resolution of the system of equations. To overcome these shortcomings, the penalty formulation is presented here. The penalty formulation provides with the possibility of imposing the incompressibility constraint without solving an auxiliary pressure equation, by replacing the continuity equation with the expression: (10) ui ,i = −εp where the so-called penalty parameter ε is a number close to zero. This equation is incorporated into the dynamic equation, and therefore a system that depends on both velocity and pressure is transformed into a velocity-dependant single equation, that converges to the fully incompressible problem as ε approaches zero [Hughes 79], [Heinrich 81], [Sohn 90], [Hannani 95]. Once we have applied the weighted residuals method, the following integral dynamic equation is obtained: 1 h h h h h h h h h h ∫ wi ui ,t + u j ui , j − f i dΩ +ν ∫ wi , j ui , j dΩ + ∫ ui ,i wi ,i dΩ − ∫ ti wi dΓ2 + Ωh ( ) Ωh ε ( ) Ωh Γ2 1   h (11) + ∑ ∫ pih  uih,t + u hj ui,jh −νui,jj − uih,i ,i − f i  dΩ = 0 ε   e Ωe which is solved for pressure. Once the velocity field has been obtained, the pressure field can be calculated as a post-processing result, by using the formula: 1 (12) p h = − uih,i ε The solution to equation (11) will approximate that of the initial problem as ε tends to zero, provided that the penalty consistency condition is verified. If not, the use of the penalty formulation could lead to the obtaining of a non-singular coefficient matrix associated to the penalty term. As ε tends to zero, this term may dominate the system of equations, therefore the whole problem could be over-constrained, and the only possible solution could be the trivial one. When carrying out an exact integration of the penalty term, ‘locking’ occurs and the only possible solution is the trivial one. This is a problem totally analogous to the one obtained when a linear velocity and a constant pressure is employed when using a mixed formulation. The discrete formulation in (11) would not be consistent and the algorithm would not achieve convergence [Hughes 79]. This problem can be avoided by making a so-called selective reduced integration of the elementary matrices involved in the resolution of the problem. A reduced numerical integration consists in using a quadrature rule that is not exact for the polynomials considered. The use of a one point Gauss quadrature rule for the integration of the quadratic functions in the penalty term, transforms the associated ‘penalty’ matrix into a rank deficient matrix and consequently ‘unlocks’ the obtaining of a non-trivial solution. For more details on this topic you can refer to [Carey 84]. Once the elementary matrices are evaluated and assembled, the integral equation in (11) can be expressed in matrix notation as: Mυ 1 ∂υ + Cυ (u, v )υ + νA υ υ + B ε υ = f ∂t ε 7 (13) where M υ is the mass matrix, Cυ (u, v ) is the convective matrix, A υ is the viscous matrix, B ε is the penalty matrix, u is the velocity vector in the x direction, v is the velocity vector in the y direction, f is external forces vector and υ is the velocity vector. The non-linearities and the derivatives with respect to time are solved by using a successive approximation method and a backward differencing scheme, in the same way as we proceeded for the mixed formulation. Segregated laminar formulation To overcome the drawbacks arising from the resolution of the integrated velocitypressure and penalty formulations of the viscous flow, the so-called segregated methods are developed in order to reduce the memory requirements when solving the NavierStokes equations. Some of the most commonly used of these segregated methods, that obtain the flow variables in a sequential way, are the fractional step method [Donea 82], [Ramaswamy 92], [Choi 97], and those based upon a SIMPLE algorithm [Benim 86], [Rice 86], [Choi 94], [du Toit 98]. An algorithm based upon the SIMPLE method, first released for finite volumes, is described in this section. The penalty method succees in solving the Navier-Stokes Equations with great memory savings due to the smaller number of equations to be solved, producing meaningful and stable solutions thanks to the use of the so-called reduced integration as seen in the previous section. Anyhow, the accurateness of the method depends on the election of the parameter ε . For very small values of ε , the weight of the penalty term in the stiffness matrix happens to cancel the amount of information contributed by the viscous term, which is very small in comparison. This information is consequently truncated and dropped from the equations. The penalty parameter should consequently be chosen depending on the word length of the computer. On the other hand, if the selected penalty parameter is too large, this choice may spoil the whole procedure, as ε is wanted to tend to zero so as to allow for convergence. Consequently, the choice of ε is not a trivial task, and a wrong selection in the parameter may lead to a meaningless solution. Moreover, the penalty formulation achieves a great reduction in the storing requirements, compared to the mixed formulation (the 2N+M equations in the mixed formulation are reduced to a 2N dimensioned system in the penalty formulation). Still, the stiffness matrix is far from being a narrow band type of matrix despite the renumbering of the nodes. Many of these shortcomings are not present in the so-called segregated methods, which are broadly employed by many authors so as to solve the Navier-Stokes equations. Following the success of the Finite Volumes Method [Patankar 80], several authors adopted the formulation in the SIMPLE and SIMPLEST methods to the finite element approach [Benim 86], [Rice 86], [Shaw 91], [Haroutunian 93], [Ferzinger 96]. These segregated finite element schemes give solution to the problem of the viscous incompressible flow, by employing a procedure in which the velocity and pressure unknowns are not obtained simultaneously but in a sequential way. The segregated formulations calculate velocities and pressures in an alternative iterative sequence, requiring much less storing needs than the mixed methods. Moreover, these algorithms not only achieve a greater reduction in the number of equations compared to the penalty method (in this formulation the dimension of the system to be solved is equal to the number of nodes), but also allow for the production of narrow band stiffness matrices, when a proper renumbering of the nodes has been carried out. The segregated method also avoids the use of the sometimes inconvenient penalty parameter. 8 Another gain of these segregated algorithms is that a mixed-order interpolation can be used [Schneider 78], [Rice 86]. As has already been said, the mixed and penalty methods require a velocity approximation different from that of the pressure. The easier-to-implement discretization of the domain in terms of the same basic functions for both velocity and pressure, leads to oscillation-free solutions, and the tendency to produce the checkerboard pressure distribution is therefore eliminated. The use of the weighted residuals method on the dynamic equation gives the dynamic system: h h h h h h h h h h ∫ wi u j ui , j − f i dΩ +ν ∫ wi , j ui , j dΩ + ∫ wi p,i dΩ − ∫ ti wi dΓ2 + Ω ( ) ( Ω Ω ) Γ2 h + ∑ ∫ pih u hj ui,jh −νui,jj + p,hi − f i dΩ = 0 e Ωe (14) The matrix version of the dynamic system can be written as: ∂N j p j dΩ C(u, v )u + νAu = Gu = f x − ∫ wi Ω ∂x ∂N j (15) C(u, v )v + νAv = Gv = f y − ∫ wi p j dΩ Ω ∂y where the right-hand side member of these equations is considered as a given value and will be taken as zero for the first iteration. We can rewrite now equation (14) as: ∂N j g ii ui = −∑ g ij u j + f xi − ∫ wi p j dΩ Ω ∂x j ≠i ∂N j (16) g ii vi = −∑ g ij v j + f yi − ∫ wi p j dΩ Ω ∂y j ≠i where g ij is the coefficient matrix of the dynamic system. Now we can express the velocities as: ∂N j  1    − ∑ g ij u j + f xi − ∫ wi Ω p d ui = j  Ω g ii  j ≠i ∂x  ∂N j  1   − ∑ g ij v j + f yi − ∫ wi (17) p j dΩ  vi =  Ω g ii  j ≠i ∂y  A set of ‘intermediate’ velocities, the pseudo-velocities u~i y v~i , are defined as follows:   1  1   − ∑ gij u j + f xi   − ∑ gij v j + f yi  u~i = v~i = (18)     gii  i ≠ j gii  i ≠ j   The velocities can be now approximated in terms of the pseudo-velocities by using the expressions: ∂N j ∂N j (19) pj ; ui ≈ u~i − K ip vi ≈ v~i − K ip pj ∂x ∂y Where the velocity-pressure coupling coefficients K ip are given by the equation: 1 (20) K ip = wi dΩ g ii ∫Ω Using the weighted residuals method in the continuity equation we obtain: h h h h h (21) ∫ wi , j u j dΩ − ∫ wi u j n j dΓ2 = 0 Ωh Γ2 or in expanded form: 9 ∂N j ∂N j ∂wi ∂wi p p N K p + N K p j dΩ = k k j k k ∫Ωh ∂x ∂x ∂y ∂y ∂wi ~ ∂wi (22) =∫ N ju j + N j v~j dΩ − ∫ wi (N j u j n x + N j v j n y ) dΓ2 Ω h ∂x Γh ∂y Once the velocity system has been calculated, the pressure system is obtained and solved, the velocities are then corrected by using the following expression: ∂N j ∂N j 1 ~ − 1 w p d Ω p j dΩ u i = u~i − w = v v (22) i j i i i ∂x ∂y g ii ∫Ω g ii ∫Ω to ensure continuity. The iterative process is based upon assuming a zero pressure field as a first guess in the resolution of the dynamic equation, providing the velocity field as the output. Once the pseudo-velocities and the pressure-velocity coupling coefficients have been calculated, the continuity system is assembled and solved, and thus the values for the pressure field are obtained. Finally, the velocities are updated, making use of the newly determined pressure field, and with both the new velocities and pressures the dynamic equations are reassembled, solved and the same procedure is repeated until convergence is achieved. When using a segregated algorithm, the use of uncoupled velocity and pressure fields may lead to the divergence of the whole process. To avoid this problem, an underrelaxation of the unknowns can be introduced so as to guarantee the convergence of the process. The linear relaxation formulae to be used for this purpose is: (23) φ n = φ n −1 + α φ n − φ n −1 where φ n and φ n−1 are the values of the unknowns (either velocity or pressure) at the present and former iterations. The momentum equations are also under-relaxed making use of an inertial relaxation factor r i defined as: ( ) ri = ∫ wi dΩ (24) Ω with r i being added to the terms in the diagonal of the dynamic coefficient matrix as follows: ∂N (g ii + ri )uin + ∑ g ij u nj = f xi − ∫Ω wi j p j dΩ + ri uin−1 ∂x j ≠i ∂N (g ii + ri )vin + ∑ g ij v nj = f yi − ∫Ω wi j p j dΩ + ri vin−1 (25) ∂y j ≠i with uin−1 and vin−1 being the values of the velocities obtained in the previous iteration. The use of this formulation, leads to a N-dimensioned narrow band coefficient matrix, and consequently to further memory savings in the resolution of the Navier-Stokes equations. Shallow Water formulation For the resolution of the Shallow Water equations we will use now a Finite Element mixed approach that will share the same advantages and shortcomings of that of the laminar mixed formulation explained before. If we apply the weighted residuals method on both the dynamic and continuity Shallow Water equations, we obtain: h h h h h h h h h h h h ∫ wi ui,t + u j ui, j − g S0i − S fi dΩ + ν ∫ wi, jui, j dΩ − g ∫ wi,i h dΩ − ∫ h ti wi dΓ2 Ωh ( ( )) Ωh 10 Ωh Γ2 +∑ ∫ p (u e Ωe h i h i ,t ( )) + u hjuih, j − νuih, jj + gh,hi − g S0hi − S hfi dΩ = 0 ∫ q (h h Ωh h ,t ) + h h uih,i + uih h,hi dΩ = 0 (26) But now we have a new difficulty that did not appear in the numerical approach to the 2D Navier-Stokes equations presented previously: we have the depth itself and the gradient of depth being included as part of the continuity equation. In fact, the inclusion of the depth and the gradient of depth in the continuity equation, allows for the verification of the conservation of mass in a pseudo-3D basis and not on a 2D laminar sense, as a consequence of having carried out an integration in depth of the NavierStokes equations. As a consequence, some pseudo-non-linearities show up in the continuity equation, which should be considered in addition to the non-linearities resulting from the convective quadratic term. The Shallow Water equations will be integrated in order to cope with this problem. Let us introduce the following approach: we are going to assume that the depth values in the continuity equation are going to be constant all over the domain for the first iteration, and equal to the outflow given depth. In the following iterations carried out in order to solve the convection, the depths and gradients of depth in the continuity equation will be evaluated from the results of the former iteration, and this evaluation will be carried out in terms of a finite difference approach. Since a non-equal order interpolation of the unknowns must be used to achieve converge, the velocities and the depths are calculated on a different mesh. The depths to be re-fed in the continuity equation for the second and the following iterations are going to be evaluated on the velocity mesh points. Recall that the basic element used in this formulation is the Q1P0 basic element, or in other words, the velocity is interpolated in terms of bilinear continuous functions with respect to a four-nodded basic element, and the depth is interpolated in terms of constant discontinuous functions within the basic element. The depth at a velocity node hi* will be taken as the mean value of the depths for the former iteration in the surrounding basic elements, and the gradients of depth on the velocity mesh h * , will be evaluated from the star depths hi* on a finite difference basis, i.e.: hi* = 1 n k ∑ hi n k =1 hi*, x = hi*, y 1 n ∑ n k =1 1 n = ∑ n k =1 h*1i1 ∑x 3 j =1 * − hi* hkij ∑y 3 j =1 − xi h − hi* kij * kij kij − yi h*1i2 k=1 k=4 h*i (27) k=3 h*1i3 k=2 where n takes a value of 1, 2, 3 or 4 depending on the velocity node being a convex corner, a side, a concave corner or an inside node, and h i j is the constant depth in the * are the star depths on the velocity nodes in the basic surrounding elements. Where hkij element k, that shares a common node (i) with hi* , xkij and y kij are the x and y coordinates of the nodes in the basic element k, that shares a common node (i) with hi* , xi and yi are the co-ordinates of the node where the gradient of depth is being evaluated and n is defined in the same way as for the star depths. The contribution to the 11 derivative with respect to x by nodes with the same abscise is ignored, and so is the contribution to the derivative with respect to y by nodes with the same ordinate, in order to avoid a division by zero. After each iteration for convection has been solved, the star depths and star gradients of the depth are calculated and re-fed into the continuity equation. The use of this algorithm developed by the authors in the resolution of the Shallow Water equations achieves very good numerical results as will be seen in the numerical examples. The general procedure for the obtaining of the steady system of differential equations could be written in its matrix form as: Cυ (u, v )υ + νA υ υ − Bh = f (28) D h * υ + E h * υ = 0 ( ) ( ) ( ) where C υ (u, v) is the convective matrix, A υ is the viscous matrix, B is the depth matrix, f is external forces vector, D(h*) is the star depth matrix, E h * is the star gradient of depth matrix, f is the external forces vector, h is the depth vector and υ is the velocity vector. SUPG stabilization of the algorithms The above formulations include a SUPG term in the dynamic weighting functions. This SUPG contribution will affect all the terms in the dynamic equation and will be based upon an optimal rule function for the obtaining of the multi-dimensional diffusion coefficient. The SUPG formulation will be the following: k uˆ hj wih, j with (29) pih = wi = wi + pi uh where the multi-dimensional definition of the diffusion coefficient k is given by: ξ uξh hξ + η uηh hη u 2 uˆ i = i , u = ui ui k = u 2 where   1  1  ξ =  cothα ξ −  η =  cothαη −  αξ  αη    u hh u hh αξ = ξ ξ αη = η η 2ν 2ν uξh = eξi ueih uηh = eηi ueih η (30) ξ 1 hη eη eξ ξ -15 x2 -10 -5 5 -1 x1 hξ Figure 1. Characteristic basic-element lengths and unit vectors. Optimal rule for the approximation of ξ y η . 12 10 15 α where hξ , hη and eξi , eηi are the characteristic basic-element lengths and unit vectors in the direction of the local axes ξ and η (see figure 1). The parameters α ξ and αη are the directional Reynolds numbers of the basic element, u eih is the velocity in the interior of the element and ν is the kinematic viscosity of the fluid. Different versions of the streamline upwind formulation have been used by other authors and can be found in [Franca 92], [Sampaio 91], [Hill 95], [Hannani 95], [Cruchaga 96], and [Choi 97]. For the present work, the stabilization technique will be based upon the streamline upwind Petrov-Galerkin weighting functions as defined in (29, 30). These weighting functions will be applied on the formulation as specified in formulas (5, 11, 14, 26), with very good results as will be seen in the numerical examples shown in the following section. NUMERICAL RESULTS Once we have presented the algorithms to be used in this work, we are going to proceed to validate them by using some benchmark problems with known solution. Validation of the laminar algorithms, the Cavity flow and the Backward-facing step benchmark problems The Cavity flow benchmark problem The driven cavity flow is a classical test used by many authors to check the quality of the methodology employed in the resolution of the 2D Navier-Stokes equations. This benchmark problem is based upon the analysis of the flow in a square cavity with prescribed horizontal velocity in the upper side and solid boundaries in the lateral and bottom sides. This is a challenging problem due to the presence of several re-circulating regions in which the solution changes rapidly, and because of the pressure singularities that show up in the upper corners. This benchmark test will be used, to validate the algorithms developed in this thesis by its comparison with the results obtained by other authors. These results to compare with, will be those of Ghia [Ghia 82] obtained by employing a mesh of 129x129 nodes; Hannani [Hannani 95] with non-uniform meshes of 32x32, 45x45 and 80x80 Q1P0 basic elements; and the results by Kondo [Kondo 91], making use of a 40x40 element mesh of four-node, non regular basic elements. All of them can be considered as reference results, specially those of Ghia, that are commonly employed to check the validity of the algorithms by most of the authors in the related bibliography. The boundary conditions used for this problem have been of the Dirichlet type in all the boundaries. A unitary horizontal velocity heading towards the right hand side has been prescribed for the top side (including the upper corners), and the no-slip condition has been considered for the rest of the boundary. The discretized domain used in the calculations has been a 1681-node non-regular mesh with 1600 Q1P0 elements (see figure 2). 13 400 350 300 Y 250 200 150 100 50 0 0 100 200 300 400 X Figure 2. Cavity Flow 41x41 non-regular mesh The results for all the formulations considered have been the same and can be seen in figures 3-6. This benchmark problem has been solved by making used of the mixed, penalty and segregated algorithms and the Reynolds numbers used have been 100, 1000, 5000 and 10000. The Reynolds number has been defined as Re = U·L / ν, where U is the velocity in the upper side, L is the length of the side of the square domain, and ν is the kinematic viscosity. The value of Reynolds = 10000 is considered as a limit for the steady Cavity Flow calculations, since has been shown through detailed numerical experiments [Shen 76] that above this bound, the stationary solution ceases to be stable. 400.00 400 400 350 350 300 300 250 250 200 200 Y Y 350.00 150 150 100 100 50 50 vel 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 300.00 250.00 200.00 150.00 100.00 50.00 0 0 100 200 300 400 0 500 0 100 200 300 400 500 X X 0.00 0.00 50.00 100.00 150.00 200.00 250.00 300.00 350.00 400.00 250.00 300.00 350.00 400.00 250.00 300.00 350.00 400.00 Figure 3. Cavity flow. Velocity field, streamlines, and pressure field Re=100 400.00 400 400 350 350 300 300 250 250 200 200 Y Y 350.00 150 150 100 100 50 50 vel 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 300.00 250.00 200.00 150.00 100.00 50.00 0 0 100 200 300 400 0 500 X 0 100 200 300 400 500 X 0.00 0.00 50.00 100.00 150.00 200.00 Figure 4. Cavity flow. Velocity field, streamlines, and pressure field Re=1000 400.00 400 400 350 350 300 300 250 250 200 Y Y 350.00 vel 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 200 150 150 100 100 50 50 300.00 250.00 200.00 150.00 100.00 50.00 0 0 100 200 300 400 500 0 0 100 200 300 X X 400 500 0.00 0.00 50.00 100.00 150.00 200.00 Figure 5. Cavity flow. Velocity field, streamlines, and pressure field Re=5000 14 400.00 400 400 350 350 300 300 250 250 Y Y 350.00 200 Level 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 200 150 150 100 100 50 50 vel 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 300.00 250.00 200.00 150.00 100.00 50.00 0 0 50 100 150 200 250 300 350 400 450 0 500 0 50 100 150 200 250 300 350 400 450 500 0.00 0.00 X X 50.00 100.00 150.00 200.00 250.00 300.00 350.00 400.00 Figure 6. Cavity flow. Velocity field, streamlines, and pressure field Re=10000 The most commonly used comparison results for this benchmark problem are the horizontal velocities along a vertical central line, which have been plotted in figure 7. Re=100 Re = 1000 Re=100 1.2 1.2 1.2 1.0 0.8 1.0 Present 41x41 Ghia 129x129 0.6 1 0.4 0.2 0.8 0.8 0.0 -0.5 0 0.5 1 1.5 0.6 Ghia 129x129 Present 41x41 0.4 0.4 Hannani 45x45 0.2 0.2 Present 41x41 Ghia 129x129 0.6 0.0 -0.5 Hannani 80x80 0 0 0.5 1 -0.5 1.5 0 0.5 1 1.5 Figure 7. Horizontal velocities along a central vertical line for Re=100, and 1000. Re = 10000 Re = 5000 1.2 1.2 1 1 0.8 0.8 Ghia 129x129 Present 41x41 0.6 0.4 0.4 0.2 0.2 0 -0.5 Ghia 129x129 Present 41x41 Kondo 40x40 0.6 Hannani 31x31 0 0.5 1 0 -0.5 1.5 0 0.5 1 1.5 Figure 8. Horizontal velocities along a central vertical line for Re= 5000, and 10000. The horizontal velocity profiles along a central vertical line adjust to the reference values of [Ghia 82], with a much finer mesh and are also substantially better than those of [Hannani 95] and [Kondo 91], for a mesh with a similar refinement and even a finer one. No substantial differences are observed among the results of the three formulations used for the velocity results nor for the pressure field results, which are also in good agreement with the benchmark solutions of the problem obtained by those authors. The good results obtained in the velocity profiles have made useless the employment of a finer mesh, which would necessitate a much longer CPU time. The calculation times are shorter for the mixed algorithm and of increasing magnitude for the penalty and segregated method. For the penalty solution, the introduction of the penalty parameter makes the linear system of equations more difficult to be solve, since the penalty parameter tends to zero. This computational time can be reduced, anyhow, by the use of a properly weighted penalty parameter. For the segregated resolution of the flow, a direct solver has been used in the calculations, with a definitively greater computational cost, and the convergence process is consequently slowed down. The algorithms implemented have proved to give very accurate results even for a less refined mesh, showing that the upwind weighting implemented in the numerical scheme 15 is a powerful tool to solve some flow problems without using very refined meshes, and with no wiggles in the so-obtained solution. The Backward-facing step benchmark problem The laminar Backward Facing Step benchmark problem is presented next, as one of the most commonly used benchmark problems in the literature, in order to validate the algorithms that give solution to the Navier-Stokes equations. The backward step is based upon a simple geometry where flow separation and reattachment occur. Experimental data for this problem can be found in [Armaly 83], who also solved this problem numerically by using a control-volume-based Finite Difference Method. The problem of the backward step flow will be solved in this section by using the penalty algorithm and its results will be compared with those of Armaly, which are generally used as verification data. The geometry and boundary conditions considered for this benchmark problem have been those used in [Armaly 83]. An expansion ratio of 1:1.94 has been considered for the widening of the channel, which has a total length of 50 so as to allow for the vortices to take place. The inlet boundary has been located at 3.5 step heights upstream of the expansion corner. The domain has been split into 2850 Q1P0 basic non-regular elements with 3021 nodes. The mesh is coarser at the outlet and more refined at the lefthand side of the channel, so as to allow for a better accuracy in the regions where the primary vortices occur. A bias parameter of 0.5 has been used for this purpose along the x-axis, therefore the width of the basic elements at the inlet is one half of that of the elements at the outlet, and the height of the basic elements is uniform within the whole domain. The mesh can be seen in figure 9, where a magnifying factor of two has been used for the y-axis. A parabolic horizontal velocity profile has been imposed at the inlet with a maximum velocity of 1, and the velocity is equal to zero at the boundaries. The lateral sides have been considered as solid boundaries and the no-slip condition has been imposed on them. Finally, a zero traction condition has been imposed at the outlet. 4 3 2 1 0 0 5 10 15 20 25 30 35 40 45 50 X Figure 9. Backward Facing Step. Mesh The flow has been obtained for a Reynolds number between 100 and 1200. The Reynolds number has been defined as Re = u·D / ν, where u is the average inlet velocity, D is the hydraulic diameter and the kinematic viscosity ν has been altered so as to make the Reynolds number vary. As foretold by the experimental results in [Armaly 83], there exists a single re-circulation zone at the expansion corner up to a Reynolds number of about 450, beyond which a second vortex shows up at the top boundary, and gets bigger as the Reynolds number is increased. 2 Y 1 0 0 5 10 15 X Figure 10. Flow in a Backward Facing Step. Streamlines for Reynolds = 400 16 3 Y 2 1 0 0 2 4 6 8 10 12 14 16 18 20 X Figure 11. Flow in a Backward Facing Step. Streamlines for Reynolds = 500 8 7 6 Y 5 4 3 2 1 0 0 5 10 15 20 25 30 35 40 45 50 X Figure 12. Flow in a Backward Facing Step. Streamlines for Reynolds = 1200 The size of the reattachment zones si versus the Reynolds number is compared with the experimental and numerical results by Armaly; these results can be seen in figures 10, 11 and 12. The reattachment locations of the vortices are defined as follows; s1 is the reattachment location of the primary vortex, s2 is the separation location of the secondary top boundary vortex and s3 is the reattachment location of the secondary vortex. All of them have been measured from the expansion corner, as depicted in figure 13. s3 s2 s1 Figure 13. Flow over a Backward Facing Step. Sketch of the recirculation lengths As seen in figures 14, 15 and 16, the computed results obtained in the present work compare more favourably with experimental data than the numerical results from Armaly. Although the present results are totally analogous to the experimental data in [Armaly 83] for s3 and for all the Reynolds numbers considered, when taking about s2 and specially s1, the experimental data differ from the calculated results beyond a Reynolds number of about 400. This difference between measured and calculated values is not only shown in the numerical results by Armaly, but also in the results by [Kim 88] and [Kwack 85] among many others. The differences in these values are due to the fact that the 3D effect becomes very important as the Reynolds number is increased. As pointed out by Armaly, these effects became predominant beyond a Reynolds number of 1300. Reattachment length s1 20 18 16 14 12 10 8 6 4 2 0 s1 Armaly exp s1 Armaly cal Present 0 200 400 600 800 1000 1200 1400 Figure 14. Reattachment length s1 versus Reynolds number for the Backward Facing Step 17 Reattachment length s2 18 16 14 12 s2 Armaly exp 10 s2 Armaly cal 8 Present 6 4 2 0 0 200 400 600 800 1000 1200 1400 Figure 15. Reattachment length s2 versus Reynolds number for the Backward Facing Step Reattachment length s3 30 25 20 s3 Armaly exp s3 Armaly cal 15 Present 10 5 0 0 200 400 600 800 1000 1200 1400 Figure 16. Reattachment length s3 versus Reynolds number for the Backward Facing Step Verification of the Shallow Water Algorithm The Shallow Water algorithm has been verified by making use of the numerical example of a channel in which the width is sharply doubled in the flow direction. The flow in the widening channel has been obtained by using the mixed Shallow Water equations with a kinematic viscosity equal to 10-6 m2/s. The channel has a length of 100 m, and spreads from 10 meters of width at the inlet, up to 20 m at a distance of 13 m from the inlet. A hydrostatic pressure of 1 m has been imposed at the outlet, and a roughness Manning coefficient of 0.01 m–1/3s has been considered throughout the channel length. The domain has been split into a 1171-node, non-regular mesh with 1090 Q1P0 basic elements (see figure 17). Frame 004  22 Jun 2000  ITERAC= 2Visc= 1.000000000000000E-006 50 Y 40 30 20 10 0 0 25 50 75 100 X Figure 17. Widening canal. Finite Elements Mesh An inflow normal uniform velocity of 3 cm/s has been imposed along the inlet and an adverse slope against the main flow direction of magnitude 10-2 has been considered all over the domain. For these flow conditions, a parallel flow is achieved at the expansion corner. As we are moving within a Froude number much smaller than the unity in all these examples, the flow may be described as subcritical with Fr<<1, where the Froude number is defined as: 18 Fr = v gh (31) where v stands for the one-dimensional velocity, g is the gravity acceleration and h is the depth. The general equation of the one dimensional gradually varied flow (see for instance [Chadwick 86]) can be written as: dh S 0 − S f = dx 1 − Fr 2 (32) where dh/dx is the variation in depth along the length of the channel, and So and Sf are the geometric and friction slopes respectively. For the present conditions, the formula 32 can be simplified into: ∆x 1 ≈ ∆h S 0 (33) The x-component of the velocity is plotted along the domain for both the twodimensional Navier-Stokes and the Shallow-Water algorithms, so as to compare them. As can be seen in plot 18, the Shallow Water algorithm behaves following (33). For the Navier-Stokes formulation the discharge at the inflow is 0.440 m3/s, and at the outlet the discharge is 0.286 m3/s. The continuity equation is not verified as this algorithm can only be applied to the simplification of a 2D flow. On the contrary, when the ShallowWater equation is used both discharges at the inlet and at the outlet are equal to 0.587 m3/s. The mass is therefore conserved along the channel. As can be seen from these examples, the 2D laminar Navier-Stokes equations can provide a good evaluation of the shallow flow in which no major changes in the depth are taking place, but when an important variation in the depth occurs, the Shallow Water formulation is a more appropriate formulation that allows for the conservation of mass throughout the domain. Frame 001  25 Jun 2000  ITERAC= Frame 001  25 Jun 2000  ITERAC= 6Visc= 1.000000000000000E-006 6Visc= 1.000000000000000E-006 Z Z Y = 3 cm/s, S = -10-2). Velocity module and depths Figure 18. Flow in a widening canal (vel 0 X Laminar Navier-Stokes (lower red surfaces) versus Shallow Water(upper blue surfaces) Y X 2 1.75 0.04 1.5 VELX 1 0.02 0.5 20 0.01 20 0.75 0 0 0.25 40 40 X 0 20 60 80 10 100 0 X 60 0 20 80 Y 10 100 19 h 1.25 0.03 0 Y APPLICATION TO SOME WASTEWATER TREATMENT PROBLEMS Once the code has been checked on some well-known benchmark problems with optimum results, it has been used to solve some real flow problems related with the civil engineering technology and in particular with the wastewater treatment industry. Some authors [Espert 96] have used the potential flow equations to evaluate the flow in clarifiers and other wastewater treatment basins. When we use these simplifications, we can obtain an approximation of the flow for slow creeping conditions, but only the resolution of the all-term-including Navier-Stokes equations will allow us to detect the real streamlines and the vortices that show up even for very slow water flows. Let us now use the previously explained formulations in the resolution of some wastewater treatment basins. Flow in clarification basins The flow of water in a rectangular and circular conventional clarification basins has been considered. Clarification has two main applications in the water treatment processes. Its most usual aim is to reduce the solids load after coagulation and flocculation have taken place. Its second application is the removal of heavy settleable solids from a turbid source to lessen the solids load in water. The aim of a good clarification basin design is the obtaining of a sufficiently stable flow, so as to achieve a better sedimentation. There is a large number of non-conventional devices for high rate clarification, such as tube or plate settlers, dissolved air flotation clarifiers, sludge blanket or slurry recirculation clarifiers. The choice of one of those depends on the features of the inflow water, the outflow water requirements, time, space and budget availability to carry out the purification of the water. The flow in a prototype of plate settlers with biological treatment, being developed in the School of Civil Engineering of La Coruña (ETSICCPC), has been also regarded. The description of the flow may be a powerful tool to attain an optimum shape in the designing of these structures, in order to make the most of the plant resources. Flow in a rectangular and circular conventional clarifiers The rectangular and circular basins are the most commonly used clarification devices. In spite of their simplicity, they have achieved excellent results with scant maintenance costs. These basins were originally designed with the capacity to store sludge for several months and were periodically taken out of service for manual cleaning. Today, most of the clarification basins include a continuous cleaning mechanical equipment, such as dragging chains that plow the sludge along the basin floor to hoppers. Nevertheless, these mobile devices for cleaning and other purposes do not have an important influence in the streamline distribution, and can be ignored when the flow is calculated (for further details on clarification basins you can refer to [Metcalf 95]). We are going to evaluate the flow features in the vertical section of a rectangular and circular clarifiers, with the following dimensions and design parameters: Clarifier Rectangular Circular Dimensions (m) Depth (m) Slope (%) Nodes&Elements Detention time (h) Surf. Load. Rate (m/h) W: 9, L: 24 3.3 1.2 1052, 949 3 1 D: 17.5 3.65 0.8 817x2, 756x2 3 1 20 Walkway Influe inta Influent intake Foam launder Baffle plate Foam sweeper Foam sweepers Overflow launder Baffle plate Foam launder Baffle slab Overflow launder Effluent pipe Sluddge scrappers Sluge hoppers Influent pipe Sluge scraper Sludge removal pipe Sludge hopper Sluge withdrawal Figure 19. Rectangular and circular clarifier with bottom sludge scraper. Sketch When working with clarifying basins, one of the criteria to be used in their definition is that of achieving a maximum head loss at the inlet, so as not to disturb the slow flow of the water mass being treated. Therefore, we should avoid turbulence by placing some kind of energy dissipating structure in the faster zone, that is the inlet (see figure 19). One of these maze-looking dissipating structures has been considered for the inlet of our clarifiers, being placed in the left-hand side for the rectangular one and in the centre of the basin for the circular, where the respective inlets are. For the outlets (in the righthand side for the rectangular one and in the perimeter for the circular one), a conventional overflow launder has been disposed. The domain in which the flow takes place has been split into 949 Q1P0 basic elements with 1052 nodes for the rectangular and 817x2 nodes and 756x2 elements for the circular, which has benefited from its symmetry property. For the working parameters chosen, a velocity of 1 cm/s has been imposed at the inlet in both of them. The no-slip condition has been imposed at the bottom and lateral sides, and the spillway has been left free with a zero traction boundary condition. For the topside, the vertical velocity has been fixed as zero and the horizontal velocity has been left free. The results for the velocity and pressure fields can be seen in figures 20 and 21. Frame 001  29 Jun 2000  ITERAC Frame 001  01 Aug 2000  ITERAC 700 600 800 700 500 Y 600 Y 500 400 400 300 300 200 200 100 0 100 0 1000 2000 X 0 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 X Figure 20. Flow in a rectangular and circular clarifying basin. Streamlines 300.00 200.00 200.00 100.00 200.00 600.00 1000.00 1400.00 1800.00 2200.00 100.00 200.00 300.00 400.00 500.00 Figure 21. Flow in a rectangular and circular clarifying basin. Pressures 21 600.00 700.00 800.00 Flow in a lamellar ‘LUPA’ clarifier The behaviour of the water flow in a lamellar settler being developed in the ETSICCPC is also regarded. In the plate settler, the settling area is increased by the disposal of several lamella plates along the basin. By doing so, the settled solids in the plates are dropped and removed, and therefore the settling rate is improved. Figure 22 shows the cross-section of a standard plate settler. The environmental and sanitary engineering area of the ETSICCPC is developing the design of a lamellar settler equipped with a permeable bio-film that carries out a biological treatment of the water. As part of this researching work, a prototype of the ‘LUPA’ clarifier (see Picture 1) has been constructed and the flow along the model has been evaluated. Discharge flumes Lamella plates Flocculation chamber Sludge hopper Figure 22. Conventional plate settler. Cross section Picture 1. The ‘LUPA’ prototype The code has been used to evaluate the flow in the interlamellar vertical section, with dimensions 80x30x10 cm3, with a water inflow at the bottom, and a free overflow at the top. The lamellas are placed at a 50.1 degree angle with respect to the horizontal. The flow has been calculated with the same assumptions to those considered for the horizontal flow clarifiers, that is, with no-slip condition on solid boundaries and vertical velocity equal to zero on the free surface. The definition of the boundary conditions is completed with a Dirichlet relationship at the inlet and a Newman equality with zero traction at the outlet. The longitudinal section of the clarifier has been divided into 1760 Q1P0 basic elements and 1916 nodes. The velocity at the inflow is set parallel to the walls. The flow has been calculated for a discharge of 15 l/day and 1500 l/day corresponding with inflow velocities of 5·10-3 and 5·10-1 cm/sg. Frame 001  30 Jun 2000  ITERAC= 1Visc= 1.000000000000000E-002 50 50 40 40 30 30 Y Y Frame 003  03 Jul 2000  ITERAC= 20 20 10 10 0 0 10 20 30 40 50 60 0 70 X 0 1Visc= 1.000000000000000E-002 10 20 30 40 X Figure 23. Flow in the lupa prototype, q = 15 and 1500 l/day. Streamlines 22 50 60 From the plots, it can observe how two primary vortices show up at both sides of the inlet for very slow flows (15 l/day). The bigger one is generated at the right-hand side, which is the lower and larger side, and a smaller one can be seen in the left-hand side. As the discharge is increased, not only the size of the primary vortices is increased, in particular for the right-hand side one, but a secondary vortex showing up, can also be seen. This secondary vortex not only alters the course of the flow, but also allows for an increase in the contact time of the water with the upper lamella, and consequently with the bio-film on it, improving the behaviour of the clarifier. Flow in a maze flocculator The flow along a maze chamber, often used in the flocculation processes, has been calculated. Flocculation is defined as the agglomeration of small particles and colloids to form settleable or filterable particles. A separate flocculation process, where chemical aids are added to water, is very often included in the treatment train to enhance contact of destabilized particles and to build dense floc particles of optimum size. The hydraulic flocculators, in opposition to the mechanical ones, allow for the formation of the flocs without the help of any mechanical device. This type of flocculation is simple and effective, especially for relatively constant flows. This sort of chambers is also used in chlorination processes. Chlorination forms part of the chemical disinfection treatments that are carried out on supply water in order to achieve its purification and transformation into drinkable water. The aim of this winding design is to achieve a slow and steady flow over a long distance to allow for the flocs to form. In chlorine disinfection processes, this slowness enables water to maintain contact with the chemical reagent over a long period of time (see [Metcalf 95] for further details on maze flocculators). The velocities involved are quite slow, and a laminar flow is expected, however, small vortices can show up and the Stokes evaluation of the flow could not detect them. For this reason, a convective-term-including formulation is required. A rectangular chamber, in which water is re-circulated along a winding path, often constitutes this kind of basins, and for this particular case will be modelled as a prismatic tank with dimensions 8m wide and 10 m long, in which a twisting channel is inscribed, split into 10 straight segments. The design parameters chosen for the chlorination tank are the following: - Tank dimensions: 8x10x2 m3 Channel width: 1m Channel length: 80 m Horizontal velocity: 6.6 cm/s Contact time: 20.2 minutes A 2091-node regular mesh with 2000 Q1P0 basic elements has been chosen to model the tank. The mixed Shallow Water algorithm has been used with a Manning coefficient of 0.012 m-1/3/s. A Dirichlet boundary condition has been prescribed at the inlet, where a parabolic velocity of 6.6 cm/s has been settled at the six lower left-hand-side nodes. At the outlet, the velocity on the six lower right-hand-side nodes has been considered as an unknown, and a hydrostatic pressure boundary condition of 2 m depth has been prescribed. A slope of 10-3 has been considered falling rightward all over the domain. A viscosity of 1.0·10-6 m2/s has been used for the wastewater. 23 Frame 002  21 Jun 2000  ITERAC 800 700 600 Y 500 400 300 200 100 0 0 250 500 750 1000 X Figure 24. Flow in a maze flocculator. Mesh As can be seen in figure 25, the results for the velocity field show the appearance of small re-circulation zones besides the corners. This re-circulation is the responsible for both the appearance of sediments in this region and also is the cause of a certain energy loss, as can be seen in the pressure plot (figure 25). These effects, if unwanted, could be removed by decreasing the velocity of the flow or the re-shaping of the channel. Frame 002  21 Jun 2000  ITERAC Frame 001  22 Jun 2000  ITERAC Z 900 X Y 800 700 2.00004 500 2.00003 Y 600 2.00002 400 2.00001 300 0 2 100 200 0 200 300 100 400 200 500 300 100 600 400 Y 0 0 250 500 750 700 500 800 600 1000 700 X X 800 900 1000 Figure 25. Maze flocculator. Streamlines and depths (in m) CONCLUSIONS In this work, an exhaustive analysis of the incompressible flow has been carried out, from the very definition of the governing equations, up to the resolution of some practical problems, passing through the comprehensive study of the numerical techniques used in their resolution. As a direct consequence of this study, a code has been written based upon this analysis, which allows for a modelling of the incompressible flow based upon a realistic interpretation of the forces taking place within the flow, and giving optimum results. The three different approaches: mixed, penalty and segregated, have been implemented and their results have been checked and verified by the comparison of the three of them among themselves and also against some reference results. As a consequence, several conclusions have been reached. The first is that, as expected, the results obtained by the three of them in the resolution of some benchmark problems have been identical, in a comparison study that had not been carried out prior to this work. The different approaches result in a different computer efficiency that depends not only on the algorithm employed, but also on the numerical solver used to obtain the solution to the resulting system of equations. Nonetheless the algorithm used does not affect the 24 accuracy of the solutions when an adequate selection of the numerical parameters has been carried out. The second conclusion is that all the results compare very favourably with the reference numerical and empirical results by other authors. As a consequence, the code not only enables a comparison study of the available finite element numerical techniques for the resolution of the Navier-Stokes equations, but also, as proved by the examples provided, contributes to a better and faster approach to these problems. A newly developed algorithm for the resolution of the Shallow Water equation, making use of the finite difference tools within the finite element frame, has been implemented with optimum results. The evaluation of this friction term is based upon on a Manning type formula, which makes use of the empirically determined Manning roughness coefficient. This term accounts not only for the energy losses that take place because of the friction with the wetted perimeter, but also for the overall turbulent losses that take place over the whole domain of integration. The turbulent eddies taking place within the flow are not detected, but the turbulent energy losses are taken into account thanks to this empirically determined formula, which provides a meaningful solution for practical flows. Some of the most commonly used hydrodynamic models for the flow calculations (such as the RMA2 model developed by the Brigham University, which is broadly used world-wide), incorporate a turbulence model featured by a constant eddy viscosity, which is not hydraulically speaking well justified. In contrast, the Shallow Water algorithm developed by the author includes an empirically determined turbulent losses term but also keeps the Navier-Stokes formulation of the problem, being ready to incorporate a k-ε turbulent model that has been developed within the research group and provides an eddy viscosity that varies in time and space. The accuracy of the numerical solutions so-obtained has been checked by using some reference benchmark numerical and empirical solutions with great success, and once the program has been validated, it has been used in the resolution of some wastewater treatment flow problems. The so defined code creates an optimum frame for the evaluation of the flow in some wastewater treatment basins, which is an essential tool in the designing of the wastewater treatment plants for the optimisation of their behaviour. The evaluation of the pressure and velocity of the flow in these basins provides very useful information about the flow properties. The data about the streamlines and velocity field distribution allows us to know where the main recirculation regions are taking place. This information will be priceless for the purpose of obtaining the geometrical parameters of the basins in order to achieve a better performance for the treatment plant. The obtaining of this optimum geometry will allow for a further recirculation, if the energy losses are required; or will enable its avoidance if unwanted, modifying in this way the detention times within the basin. The velocity and pressure fields also provide invaluable information about the distribution of the discharge among the outlets, which again can be redefined in order to improve the behaviour of the plant. Thanks to the information obtained by this numerical evaluation of the flow, the water treatment basins and channels can consequently be designed to fit the requirements of the processes being carried out. REFERENCES  Armaly B.F., Durst F., Pereira J.C.F. and Schönung B., Experimental and theoretical investigation of backward-facing step flow. J. Fluid Mech. 127 (1983), pp 473-496. 25                   Babuska I., Error bounds for finite element method, Numerische mathematic 16 (1971), pp 322-333. Baker A.J., A finite element computational theory for the mechanics and thermodynamics of a viscous compressible multi-species fluid. Bell Aerospace Research Rept. 9500-920200 (1971). Benim A. C. and Zinser W., A segregated formulation of Navier-Stokes equations with finite elements. Comput. Methods Appl. Mech. Engrg. 57 (1986) 223-237. Brezzi F., On the existence, uniqueness and approximation of saddle point problems arising from Lagrange multipliers. Rev. Française Automatique Informatique Reserche Operationnelle, Ser. Rouge Anal. Numér., 8 (R2) (1974), pp 129-151. Brooks A.N. and Hughes J.R., Streamline Upwind / Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg 32 (1982), pp 199259. Carey G. and Oden J., Finite Elements. Prentice-Hall (1984). Chadwick A. and Morfett J., Hydraulics in Civil Engineering. Allen & Unwin (1986). Choi H.G. and Yoo J.Y., Streamline upwind scheme for the segregated formulation of the Navier-Stokes equation. Numerical Heat Transfer 25-Part B (1994), pp 145161. Choi H. G., Choi H. and Yoo J.Y., A fractional four-step finite element formulation of the unsteady incompressible Navier-Stokes equations using SUPG and linear equal-order element methods. Comput. Methods Appl. Mech. Engrg. 143 (1997), pp 333-348. Cruchaga M. A. and Oñate E., A finite element formulation for incompressible flow problems using a generalized streamline operator. Comput. Methods Appl. Mech. Engrg. 143 (1997), pp 49-67. Donea J., Giuliani S., and Laval H., Finite Element solution of the unsteady Navier-Stokes equation by a fractional step method. Comput. Methods. Appl. Mech. Engrg. 30 (1982), pp 53-73. Espert V., García M., Sancho H. and López A., Modelo matemático bidimensional para el estudio de flujo de agua a través de un decantador rectangular con lamelas. Ingeniería del Agua 3 (1996) pp 16. Ferziger J.H. and Peric M., Computational methods for fluid dynamics. Springer (1996). Fortin M., Analysis of the convergence of mixed finite element method. RAIRO Anal. Numer. 11 (1977), pp 341-354. Franca L. P. and Frey S. L., Stabilized finite element methods: II The incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg 99 (1992), pp 209-233. Gartling D.K., Finite element analysis of viscous incompressible flow. Ph.D. dissertation. University of Texas at Austin, Austin (1974). Ghia U., Ghia K. N. and Shin C.T., High Re solutions for incompressible flow using the Navier-Stokes equation and the multigrid method. J. Comput. Phys. 48 (1982), pp 387-411. Hannani S.K., Stanislas M. and Dupont P., Incompressible Navier-Stokes computations with SUPG and GLS formulations- A comparison study. Comput. Methods Appl. Mech. Engrg 124 (1995), pp 153-170. 26                    Haroutunian V., Engelman M. S. and Hasbani I., Segregated finite element algorithms for the numerical solution of large scale incompressible flow problems. Int. J. Numer. Meth. Fluids 17 (1993), pp 323-348. Heinrich J.C., Huyakorn P.S., Mitchell A.R. and Zienkiewicz O.C., An upwind finite element scheme for the two dimensional convective transport equation. Int. J. Numer. Methods Engrg 11 (1977), pp 131-143. Heinrich J.C., and Marshall R.S., Viscous incompressible flow by a penalty function Finite Element Method. Comput. Fluids 9 (1981), pp 73-83. Hughes T.J.R., Liu W.K. and Brooks A., Review of Finite Element analysis of incompressible viscous flow by the penalty function formulation. J. Comput. Phys. 30 (1979), pp 1-60. Jamet P., and Raviart P.A., Numerical solution of the stationary Navier-Stokes equation by Finite Element Methods. Computing methods in applied sciences and engineering (Ed. By Glowinski and J.L. Lions). Springer, Berlin (1973), pp 193223. Kim S.W., A fine grid Finite Element computation of two-dimensional high Reynolds numbers flows. Computers and Fluids 16-4 (1988), pp 429-444. Kondo N., Third-Order upwind finite element formulations for incompressible viscous flow problems. Comput. Methods Appl. Mech. Engrg. 93 (1991), pp 169187. Kwak D., and Chang J.L.C., A three dimensional incompressible Navier-Stokes flow solver. Part I. –INS3D Code. CFD Workshop. University of Tennessee Space Institute, Tullahoma, TN (1985). Ladyzhenskaya O., The mathematical theory of viscous incompressible flow. Gordon & Breach (1969). Laval H. and Quartapelle L., A fractional-step Taylor-Galerkin method for unsteady incompressible flows. Int. J. Numer. Meth. Fluids 11 (1990), pp 501-513. Metcalf & Eddy INC., Ingeniería de aguas residuales, tratamiento, vertido y reutilización. Mc Graw Hill (1995). Oden, J.T., Finite Elements of non-linear continua. McGraw Hill, New York (1972). Patankar S.V., Numerical heat transfer and fluid flow. McGraw-Hill (1980). Ramaswamy B., and Jue T.C., Some recent trends and developments in Finite Element analysis of incompressibe thermal flow. Int. J. Numer. Methods Engrg. 35 (1992), pp 671-707. Rice J.G. and Schnipke R. J., An equal-order velocity pressure formulation that does not exhibit spurious velocity modes. Comput. Methods Appl. Mech. Engrg. 58 (1986), pp 135-149. Rodi W., Turbulence models and their application in hydraulics: A state of the art rewiew. Balkema, Rotterdam (1993). Sampaio P. A. B. de, A Petrov-Galerkin formulation for the incompressible NavierStokes equations using equal order interpolation for velocity and pressure. Int. J. Numer. Meth. Engrg 31 (1991), pp 1135-1149. Schneider G.E. Raithby G.D., and M.M. Yovanovich, Finite Element solution procedures for solving the incompressible Navier-Stokes equations using equal order variable interpolation., Numer. Heat Trans. I (1978), pp 433-451. Shaw C.T., Using a segregated finite element scheme to solve the incompressible Navier-Stokes equations. Int. J. Numer. Meth. Fluids 12 (1991), pp 81-92. 27            Shen S. F.and Habashi W., Local linearization of the Finite Element methods and its application to compressible flow. Int. J. Numer. Methods Engrg. 10 (1976), pp 565. Sohn J.L., and Heinrich J.C. A Poisson equation formulation for pressure calculations in penalty Finite Element models for viscous incompressible flows. Int. J. Numer. Meth. Eng. 30 (1990), pp 349-361. Taylor C. and Hood P., A numerical solution of the Navier-Stokes equation using FEM technique . Compt. &Fluids I (1973), pp 73-100. Temam R., Navier-Stokes equations. North Holland, Amsterdam (1977). Toit C.G. du, Finite element solution of the Navier-Stokes equations for incompressible flow using a segregated algorithm. Comput. Methods Appl. Mech. Engrg. 151 (1998), pp 131-141. Turner, M., Clough R., Martin H., and Topp L., Stiffness and deflection analysis of complex structures, J. Aero Sci. 23/ 9 (1956), pp.805-823. Vellando, P.(Pablo Rodríguez-Vellando Fernández-Carvajal). On the resolution of the Navier-Stokes equations by the Finite Element Method using a SUPG stabilization technique. Application to some wastewater treatment problems. Dotoral thesis. Universidad de La Coruña (Spain), Marzo, 2001. Versteeg H.K. and Malalasekera W., An introduction to computational fluid Dynamics. Harlow, UK (1995). Weiyan T., Shallow water hydrodinamics. Elsevier (1992). Zienkiewicz O.C., Constrained variational principles and penalty function methods in Finite Element analysis, in Lecture Notes in Mathematics, Conference of the Numerical Solution of Differential Equations, ed. G. A. Watson. Springer-Verlag, Berlin (1974), pp 207-214. Zienkiewicz O.C., R.H. Gallagher and P. Hood, Newtonian and non-Newtonian viscous incompressible flow. Temperature induced flows. Finite element solution, J.R. Whiteman, ed. The Mathematics of Finite Elements and Applications II (MAFELAP 1975). Academic Press, London (1976). NOTATION Aυ K ip m M Mυ n Viscous coefficient matrix B Gradient of pressure matrix Penalty matrix Bε bi Prescribed velocity Cυ (u, v ) Convective coefficient matrix eξi , eηi Unit vectors in the ξ and η directions f fi g h h h* h,*j hξ , hη k ij k ij N Ni p p p q Q1P0 r Re Rh S0 Sf t Body force vector Body forces components Gravity force Depth vector Depth or size of the grid Star depth Characteristic lengths along ξ and η Star gradient of depths Diffusion Artificial diffusion 28 Pressure-velocity coupling coefficient Mass Number of pressure nodes Unsteady coefficient matrix Manning roughness coef. or outward unit vector normal to the interface Number of velocity nodes Velocity shape function Pressure vector Pressure SUPG weighting term Weighting function Bilinear velocity-constant pressure BE Inertial relaxation factor Reynolds number Hydraulic radius Geometric slope Friction slope Time ti u u u~i v v v~i Velocity along the y direction Velocity vector along the y direction Pseudo-velocity in the y direction αp Weighting functions Pressure relaxation parameter wi η µ Traction vector Velocity along the x direction Velocity vector in the x direction Pseudo-velocities in the x direction ν ξ ρ σ ij τ ij υ χi Velocity relaxation parameter αu α ξ ,αη Reynolds numbers along ξ and η Γ ε Ω ∂Ω Boundary of the Ω domain Penalty parameter 29 Local spatial variable Dynamic viscosity of the fluid Kinematic viscosity of the fluid Local spatial variable Density of the fluid Stress along the boundary Stress Velocity vector Pressure shape function Domain of integration Domain boundary