Finite Element Solution Of The Two-dimensional .

Available athttp://pvamu.edu/aamAppl. Appl. Math.ISSN: 1932-9466Applications and AppliedMathematics:An International Journal(AAM)Vol. 13, Issue 1 (June 2018), pp. 535 - 565Finite Element Solution of the Two-dimensional IncompressibleNavier-Stokes Equations Using MATLAB1*Endalew Getnet Tsega and 2V.K. KatiyarDepartment of MathematicsIndian Institute of Technology RoorkeeUttarakhand, India1endalebdumath2016@gmail.com, 2vktmafma20@gmail.com*Corresponding AuthorReceived: July 5, 2017; Accepted: May 28, 2018AbstractThe Navier–Stokes equations are fundamental in fluid mechanics. The finite element methodhas become a popular method for the solution of the Navier-Stokes equations. In this paper,the Galerkin finite element method was used to solve the Navier-Stokes equations for twodimensional steady flow of Newtonian and incompressible fluid with no body forces usingMATLAB. The method was applied to the lid-driven cavity problem. The eight-nodedrectangular element was used for the formulation of element equations. The velocitycomponents were located at all of 8 nodes and the pressure variable is located at 4 corner ofthe element. From location of velocity components and pressure, it is obvious that thiselement consists of 16 unknowns for velocities and 4 unknowns for pressure. As a result, theunknown variables for velocities and pressure are 20 per each element. The quadraticinterpolation functions represent velocity components while bilinear interpolation functionrepresents pressure. Finite element codes were developed for implementation. The numericalresults were compared with benchmark results from the literature.Keywords: Navier–Stokes equations; finite element method; steady-state solution; eightnode rectangular element; MATLABMSC 2010 No.: 35Q30, 35Q35, 76D05, 76M10535

536Endalew Getnet Tsega et al.1. IntroductionThe Navier-Stokes equation is a set of nonlinear partial differential equations that describethe flow of fluids, which represent conservation of linear momentum. It is the cornerstone offluid mechanics as noted by Cengel et al. (2010). It is solved jointly with continuity equation.These equations cannot be solved exactly. So, approximations and simplifying assumptionsare commonly made to allow the equations to be solved approximately. Recently, high speedcomputers have been used to solve such equations by replacing with a set of algebraicequations using a variety of numerical techniques like finite difference, finite volume, andfinite element methods.Finite element method is the most powerful numerical technique for computational fluiddynamics which is readily applicable to domains of complex geometrical shape and providesa great freedom in the choice of numerical approximations. It reduces a partial differentialequation system to a system of algebraic equations that can be solved using traditional linearalgebra techniques. In finite element method, the domain of interest is subdivided into smallsubdomains called finite elements. Over each finite element, the unknown variable isapproximated by a linear combination of approximation functions called shape functionswhich are associated with the node of the element characterize the element. The piecewiseapproximations for elements are assembled together to obtain a global system to the wholedomain. One of the major advantages of the finite element method is that a general purposecomputer program can be developed easily to analyze various kinds of problems as noted byKwon et al. (1997). In particular, any complex shape of problem domain with prescribedboundary conditions can be handled with ease using the finite element method.Jiajan (2010) discussed the Galerkin finite element formulation of two dimensional unsteadyincompressible Navier-Stokes equations using the quadratic triangular element (6-nodes).The pressure variable was located at the corner nodes and the velocity components werelocated at all of the six nodes. Ghia et al. (1982) studied high Reynolds number solutions forincompressible flow using the Navier-Stokes equation and the multigrid method. Persson(2002) implemented a finite element based solver of the incompressible Navier-Stokesequations on unstructured two dimensional triangular meshes. He solved the lid-driven cavityflow problem for four different Reynold’s numbers: 100, 500, 1000 and 2000.Glaisner et al. (1987) discussed finite-element procedures for the Navier-Stokes equations inthe primitive variable formulation and the vorticity stream-function formulations. Steadystate solution of lid-driven cavity flow was obtained by the velocity-pressure formulationusing the nine-node rectangular element in their work. Taylor et al. (1981) and Smith et al.(2014) used eight-noded rectangular element mesh to solve two-dimensional incompressibleNavier-Stokes Equations with FORTRAN programming language. Rhaman et al. (2014)presented Galerkin finite element method to simulate the motion of fluid particles which satisfies theunsteady Navier-Stokes equations through a programming code developed in FreeFem .Simpson (2017) used nine noded rectangular elements with two degree of freedom on eachnode for finite element simulation of a coupled reaction-diffusion problem using MATLAB.Khennane (2013) developed MATLAB codes for 4-nodded and 8-noded quadrilateralelements for the linear elastic static analysis of a two dimensional problem using finiteelement method.

AAM: Intern. J., Vol. 13, Issue 1 (June 2018)5372. Governing EquationsFor steady Newtonian incompressible fluid with no body forces, the governing equations fortwo-dimensional flow are:Continuity equation: u v 0. x y(1)Navier-Stokes equation:x-component 2u 2u u u p v 2 2 , y x y x x u(2)y-component 2v 2v v v p v 2 2 . y y y x x u(3)where u and v are the x, y components of the velocity vector, p is static pressure, is densityand dynamic viscosity of the flowing fluid.Using L and V as a characteristics (reference) length and velocity respectively, we define thedimensionless variablesx* xyuvp, y* , u* , v* , and p* .LLVV V 2The governing equations, Equation(1), Equation (2),and Equation (3) can be written in theirdimensionless form (ignoring the astrix ‘*’) as: u v 0, x y(4) u u p 1 2 u 2 u ,u v x y x Re x 2 y 2 (5) v v p 1 2 v 2 v , v x y y Re x 2 y 2 (6)uwhere VL is the Reynolds number.Re

538Endalew Getnet Tsega et al.3. Formulation of Element EquationsTo describe the structure of finite element programming of a steady-state solution of theNavier–Stokes equations, let us consider a flow confined to a rectangular cavity driven by auniform horizontal velocity at the top. The velocities at the other three boundaries are set tozero. Eight-noded rectangular elements are used to model the flow. We use all the 8 nodes ofeach element to model velocities components u and v and the 4 corner nodes to model thepressure p. Meshes are numbered in x-direction. Freedoms numbered are in the order u p vas used by Smith et al. (2014) and Taylor et al. (1981).Figure 1. Lid-driven cavityLet us denote an element by . Shape functions for the rectangular elements are expressed interms of local coordinates and where 2(x-xc)/lx, 2(y-yc)/lyHere, (xc, yc)element, andlength in xu,vu,v,pu,vis the centroid of thelx, ly represent itsand y-direction.u,v,pu,vu,v,pu,vu,v,p

AAM: Intern. J., Vol. 13, Issue 1 (June 2018)539Figure 2. Eight-noded rectangular elementFigure 3. Eight-noded rectangular elements meshSuppose the nodes 1, 2, 3, 4, 5, 6, 7, 8 have coordinates (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1),(0, 1), (-1, 1), (-1, 0) in the local coordinate system. Then, the general form of shape functionsfor 4-noded bilinear rectangular element (considering corner nodes) using local coordinates isM a b c d .(7)The general form of shape functions for 8-noded quadratic rectangular element (consideringall nodes) using local coordinates isN a b c d 2 e f 2 g 2 h 2 .(8)Using Kronecker-delta property of shape functions, from Equation (7) and Equation (8) shapefunctions for 4-noded and 8-noded rectangular elements are

540Endalew Getnet Tsega et al.1(1 ),41M 2 (1 ),41M 3 (1 ),41M 4 (1 ).4M1 (9)N 1 14 (1 )(1 )(1 ),N 2 12 (1 2 )(1 ),N 3 14 (1 )(1 )(1 ),N 4 12 (1 )(1 2 ),(10)N 5 14 (1 )(1 )(1 ),N 6 12 (1 2 )(1 ),N 7 14 (1 )(1 )(1 ),N 8 12 (1 )(1 2 ).Thus, the quadratic interpolation functions are used for velocity components while bilinearinterpolation functions for pressure. As a result, the unknown variables for velocities andpressure are 20 per each element. Thus, the dependent variable u, v, and p are expressed as8u N i ui ,i 18v N i vi ,(11)i 14p M i pi .i 1where u i , vi and pi are velocity and pressure values at the nodes.Now, expressing Equation (4), Equation (5), and Equation (6) using these shape functions weget8 N jj 1 x8uj j 18 Nk 18k N juk j 1 y N j xvj 0 ,8u j N k vkk 1(12)8 j 1 N j y28 2N 1 8 N jj , u u jj2 Re j 1 x 2j 1 y 4u j l 1 M lpl x(13)

AAM: Intern. J., Vol. 13, Issue 1 (June 2018)8 Nk 18kuk N jj 18v j N k vk xk 1541 N j8 4v j yj 1l 1 M lpl y(14)8 N 1 8 N jj . v v jj2 Re j 1 y 2 yj 1 22Employing Galerkin weighted residual approach on Equation (13), we get 8NN k uk i k 1 N j8 xj 18u j N k vkk 18 N j yj 14uj l 1 M lpl x(15)28 2N 1 8 N jj dA 0 u u jj2 Re j 1 x 2j 1 y Using Gauss-Divergence Theorem, we have8 2 N8 2N 1 1 N ijj Nu Nu j dA ji 22 Re i Re xj 1 y j 1 x 8 N1j Nu j dSi Re j 1 n8 N j xj 1u j dA N i y8 j 1 N j u j dA y where is the boundary of the element , n n x , n y is the unit outward normal vector to N j N j N jnx n y is the directional derivative of N j in the direction n x ynormal to the boundary . Hence, we havethe element and8 NN k uk i k 1 N j8 xj 11 N i Re x8 8u j N i N k vkk 1 N j N iuj x yj 18 8 N j xj 14u j Ni l 1 M lpl x8 N 1j u j dA dA Ni u j dS 0. yRe j 1 n N jj 1i, j 1,2, ., 8.For Dirichilet boundary condition, we have8 NN k uk i k 18 N j xj 11 N i Re xi, j 1,2, ., 8.8 j 18u j N i N k vk N jk 1 N iuj x y8 j 18 j 1 N j x4u j Ni l 1 u j dA dA 0. y N j M lpl x(16)

542Endalew Getnet Tsega et al.Applying similar procedure for Equation (14), we get8 NN k uk i k 18 N j xj 11 N i Re x8 j 18v j N i N k vk N jk 1 N ivj x y8 j 18 N jj 1 x4v j Ni l 1 M lpl y v j dA dA 0. y N j(17)i, j 1,2, ., 8.Employing Galerkin weighted residual approach on Equation (12) using the weight functionsMl , we get8 N 8 N j jMu v j dA 0, j l j 1 xj 1 y or8 N8 N jjMu Mv dA 0. ljlj j 1 xj 1 y l 1, 2, 3, 4 .(18)Due to the nonlinearity, the set of algebraic equations that will be obtained here cannot besolved in a single shot, but an iterative solution is necessary. In such an iterative solutionnonlinear terms can be linearized in a number of different ways. The simplest possibility,which will be used in this paper, is known as Picard linearization, in which the nonlinearterms are replaced by u u u u v u v, x y x y v v v vu v u v. x y x yuwhere u and v are approximate values for the velocity components.We assume starting values u 01 , u 02 , . . ., u 08 and v 01 , v 02 , . . ., v 08 for the element and8u N k u0k ,k 18v N k v0 k .k 1The iteration process continues by replacing u 0 k and v0 k , k 1, 2, ., 8, by the average ofvelocity component values from the previous two iterations until tolerance is satisfied.From Equation (16), Equation (17) and Equation (18), we get a system of equations in matrixform as