A mixedfinite element methodfor the solution of the Navier-Stokes equations is presented. The scheme uses discontinuousfinite elements that allow the utilization of upwind derivatives. This makes the scheme more stable and accelerates the convergence. Moreover the boundary layer is better consideredfor high Reynolds numbers. The existence of a solution and the convergence of the se heme are established. Results inside a square cavity far Re = 10000 far the recirculating jlow, the stream function and their second derivatives, viscous stress tensor and velocity vector are obtained. Between Re = 1000 and Re = 2000 a third secondary vortex appears nea, the upper right comer

Figueroa Nieto, Jaime

Contribuciones Científicas y Tecnológicas

Contribuciones Científicas y Tecnológicas Nº 66 - Area Matemática VI, Mayo 1984. Mixed finite element solution for the navier-stokes equations* Jaime Figueroa Nieto•• SUMMARY: A mixedfinite element methodfor the solution of the Navier-Stokes equations is presented. The scheme uses discontinuousfinite elements that allow the utilization of upwind derivatives. This makes the scheme more stable and accelerates the convergence. Moreover the boundary /ayer is better consideredfor high Reyno/ds numbers. The existence of a solution and the convergence of the se heme are established. Results inside a square cavity far Re~ JOOOO far the recirculating jlow, the stream function and their second derivatives, viscous stress tensor and velocity vector are obtained. Between Re= 1000 and Re=2000 a third secondary vortex appears nea, the upper right comer. 1. INTRODUCCION The main difficulty of the numerical solution of the Navier-Stokes equations for an incompressible viscous fluid is to simulate the behaviour when viscosity becomes small , or equivalently, when the Reynolds numbers becomes large. In this case near the walls zones will appear where the gradient of the velocity is large and the viscous effects are importan!. One of the models proposed to describe separated flow corresponding to small viscosity is that dueto G. K. Batche-lor (1956). The model postulates at high Reynolds numbers the existence of a central core of constant vorticity sorroun-ded by a thin viscous !ayer. In studies done by Burgraff (1966) and receotly by Nallasamy & Krishna ( 1978) provi-de comparisons of the uniforrn vorticity model with nume-rical solutions of the Navier-Stokes equations. The square cavity with a steady sliding top wall is a geometrically simple case and for which the recirculating flow has no analytical solution, is a prototype of separated flow for which Batchelor's model would be applicable. From a computational point of view the cavity flow is an * Manuscrito revisado y aprobado en forma definitiva en sep-tiembre de 1982. ** Departamento de Matemática y Ciencia de la Computación, Facultad de Ciencia, Universidad de Santiago de Chile (USACH). easy model problem which allows comparisons among numerical methods. The method that will be presented was introduced by Hermann (1) in 1967 to solve the bending plate problem. Other forrnulations of the method were given by Hellan (2) in 1967 and Yisser (3) in 1969. In 1973 , Johnson (4) justified the method proving the existence of a solution and obtaining an estímate of the error. In 1976, Brezzi and Raviart (5) proved that the method is the best in the sense that the order of convergence is O(hk+ 1) with k the degree of polynomials emrloycct. The stationary Navier-Stokes cquatitln , in a boundcd domain í1 C IR" are: { -v1 ¡j + ¡j . 'v ¡j + 'vp = 1 in n ( l. 1) div u = O in í1 u = o on an where v is the kinematic viscosity; p the pressure; u = (u;)1c ;c n the velocity vector; 1 = (f;)1c ;c n the body externa( forces. A classical variational forrnulation of the problem is (Cf.(6)): (1.2) { Find fl a pair (u, p) in Y X L2(í1)/ IR satisfying: n + ++ . + v~ = I fooiu;oiv;dx +(u · v' u, v)- (p,d1v v) = (1, t); 't/ t E (Hó (Ü))" where Y = { t E (Hó (í1))": div t = O}. There is a solution (u, p) in y X L2(0)/IR of (1.2) (equiventely of ( 1 . 1)) if n is a bounded domain of IR", n ,s; 4, with a Lipschitz _continous boundary (Cf.(7)). The solution is unique if the viscosity is relatively large with respect the norrn of the body externa( forces. In two dimensions the Navier-Stokes equations can be stated in terrns of the stream function as follows: 25 Toe problems ( l. 2) and ( l. 1) are equivalen! . That is: - if (u, p) is the solution of ( 1.2) , the only function 'l'EHij(Ü) such that u= c-¡;l'I' = (il2'1' , - il 1'l') is the solution of (1 .3); - if 'I' is the solution of (1 .3) there exists pin L2(.0.)/IR such that the pair ( ti, p) = (c-¡;l'I' , p) is the solution of (1.2) . To interpret the problem (1 .3), the function q, satisfies { v.:\ 2'l' - a 1(.::\ '1'a2 'l') + ai(.::\ q, a 1 'l') = curl f in f1 aw (1.4) 'I' = - = Oon a.o. éln 2. FORMULA TION OF THE METHOD FOR THE STOKES PROBLEM The Stokes problem has the following fonnulation: (2. 1) { Find A in Hij(Ü) satisfying V .::\2A = curl 1 in n Let A in Ha(fl) be a solution of (2. 1 ). lf we set <T;j = ilfj A, J ,,,;; i, j ,,,;;2 then <T = (<T;j)1,,. ;, j,,.2 and A satisfy the equations 2 ;,~ 1 JO (<Tij - afjA) T;j dx = O 2 i.~1 JO <T;j afjµ dx - JO curlf dx = O Vµ E Ha(fl) We assume that !1 is a convex polygon . Let th be a triangulation of .O. by triangles K with sides of length less than h . Let T = (Tjj) 1,,. ;.j,,.2 in (L 2(0.))4. T12 = T21 so that T;/ K E H 1(K ), V K E th, then 2 i.j a ) T jj ªDA dx = - ¡ ¡ I 2 0 K~h i.j =1 J K aj T;j él;X dx + 26 .. .. where n = (n1, n2) is the unit normal vector and t = (n2, - n1) is the unit tangen! along the boundary aK. If we impose the condition denoted (MN): 2 . ¡ .. K1 K, _ ,.¡=1 T,/ K 1 n; nj -then a = (<T,j)1,,.;, j .. 2 and A satisfy the system { a(a, T) + b(T, A)= O ; VTEY._ (2.2) v b(<T, µ.) + f O µ. curlf dx = O Vµ.dl~O.) with 1 < s < 3/2 where the space 2 = { T = (T;j) t,,.i. j"'2 E (L2(!l))4, T1 2 = T21 , T;/K E H1(K) , VKEth,} , T satisfy the condition (MN) a(a, T) = . _l J <T;j T;j dx l.) ª 1 o b(a, µ ) = t. _'l <f aj Tj T;i il;µ dx Kt , ,.¡ • I K -f ªaµ. T¡j n; lj dS) aK t Brezzi and Raviart (5) proved that the method (2 .2) is optima!. 3. FORMULA TION OF THE METHOD FOR THE NA YIER-STOKES EQUATIONS We shall consider the fonnulation (1.2) . lt is easy obtain "centered" fonnulation of the method by introducing the tenn in the second equation of (2.2). But the problem is to simulate efficiently the behaviour of the fluid when the viscosity becomes small . For that Fortín (8) proposed to consider fonnulations called " de-centered''. Toe method is a generalization of that introdu-ced by Lesaint (9) for the approximation of hyperbolic equations of first degree by means of discontinuous finite elements. Toe difficulty is the discretization of the non linear tenn fo (u · 'vu)~ dx where U = CltiÜ, ~ = cut! µ. Since the functions A, µ. e HM'I'), 1 <s<3/2 are not sufficiently smooth. For that we introduce the space Mh defined as Mh = Hi(O) n {µ.eCº(O): µ./K e Pk , v'Ke,h}, with kelN . + -+ + + éJ). lf A, µ. e Mh and u = curl A then u · n = - "ar, that is u· ñ is continuous across the boundary aK, i.e. u· ñK = + + ' - u · nK, on K 1nK2 -f 0. This allows us to define a partition of the boundary ilK = ilK+ n ilK_ with ilK + the zone where the fluid exits, that is u· ñ> O and aK_ = aK -ilK+. By definition { ~K on ilK+ Uc = ~K· on aK_ n OK1 ' K 1 next to K We consider the next discretization of the nonlinear convective tenn by using upwind derivatives: 2 ;JI J u;uj(iljv;) dx K + l I · + ...... K•{h (u · v) ( u · n) ds '" denoted d(A, A, µ.) and introduced in the second equations of (2.2). We obtain the mixed fonnulation { Find a pair (<7, >,.) in Vh x Mh such that a(<Y, T) + b(T, A) = O ; v'TeVh v b(<7, µ.) + d(A, A, µ.) (3.1) + f n µ. curl f dx = O ; Vµ.eMh where Conceming the existence of a solution for the system (3 . 1) we have the next result (Cf. ( 10), ( 11)). Denoting Fh:Mh -+ Vh ; µ.h -+ ah where <Yh is the solution ofthe problem a(ah, Th) = - b(Th, µ.h) , VTh e Vh we have Theorem 3.1 lf i) d(µ.h , µ.h, µ.h) ~ O ; Vµ.h e Mh (3 .2) ii) 11 Ah 11 o.n ~ C II Fh(Ah) 11 o.n ; v'Ah E Mh (3.3) there exists at least one solution (CTh, Ah) in Vh x Mh for the system (3.1). Conceming the convergence ofthe scheme we have the following result (Cf. (10), (11)): Theorem 3.2 Suppose that (3.2) and (3.3) hold for ali h> O. Then there exista pair (CT, A) in (L2(0))3 x H'(íl) anda subsequence of the sequence (<Th, Ah) denoted (ah , Ah) where (ah, >..h) is a solution of (3. 1) so that ah-+ <7 weakly in (L2(0))3 Ah-+ A strongly in H1(0) . Theorem 3.3 The function A belongs to Hij (O) and is the solution ofthe Navier-Stokes equations in tenns of the stream function. Remark The error estímate for this type of scheme is an open problem. 4. THE METHOD OF SOLUTION If the boundary an is moving with a relative velocity g in H 112 ( ilil) then the conditions on the boundary are u · ñ = O and u· i = g. That is~ = O ~ = g . Let >. 0 1n a, ' an H2 (O) be such that >.0 = O and : ~, = g over ilil. Toen we write the problem Find >. e H2(0), >. - >.o in Hij(il) such that - v J 6 X61P dx + f 6A(il 1Ail21P - il2Xil 11P) dx + n n J 1P curl 7 dx = O; VIP e H2(0) 11 The boundary condition ~ = g is introduced in the first equation of (3.1): iln 1,J - 1 a(CT, T) + b(T, A) = Í J T;in;ni g ds "'' We shall use the iterative method 1) (ot >.~ is the initial value in Vh x Mh 2) (8h , lj,h) in Vh x Mh is the solution of a(9h, Th) + b(Th, lj,h) = Í J (Th);in;ni g ds; V Th E Vh 1,1• 1 ñU (4. 1) vb( eh' µ.h) + d( xg' lj,h' µ.h) + f µ.h curl f dx = O; Vµ.h E Mh 11 3) (CTg + 1, >,g+ 1) = (CTg, Ag) + "{ (9h - <Tg, lj,h - >,.g) with O < -y < 2. Remarks conceming the iterative method The existence theorem in the paragraph 3 may be adapted if we consider d(~ . >-., µ. ) with ~ fixed and the hypothese 27 STREAM LINES RE = 2000 1 d, FIG. 1 OEFINITION SKETCH TP PRIMARY VORTE X LL LOWER LEFT VORTEX LR LOWER R IGHT VORTEX UR UPPER RIGHTVORTEX h GRID SIZE Fig. l. Definition sketch. d(X.h, '1ih , '1ih) :S O, V IJih in Mh . Then there exista so lution for the linealized system (4.1) . The unicity is easy to see . Fortin and Thomasset (12) used d(>..~. >..~. µ.h ). The matrix is the same each iteration. The system (4.1 ) is solved by using the frontal method (Cf. (13)). Toe difficulty of the linearization (4 . 1) is the increase in the band with and that the matrix is different with each iteration. However the convergence is rapid. 5. PRACTICAL IMPLEMENTATION OF THE METHOD lt was demonstrated (5) that the degrees of freedom by e lement and for the space V h are Jk(~+ I) . Therefore , the total degrees of freedom by element is DL(k) = 3k(k + 1) 2 + (k + 1) (k+ 2) 2 Toe result on the degrees of freedom of a tensor-valued function 7 E Yh: "Let K be in 7h and let 7 = (7;j)1:s i. j:s 2 be a tensor-valued function such that k 7;i E PK - I and 712 = 721. Then 7 is uniquely determined by the values of 28 K: J 7;i q dx ; q E Pk - 2, 1 :S i, j :S 2 where K; , i = 1, 2, 3 denotes the sides of K · · was also established in (5). We refer to ( 11 ) for technical details . 6. NUMERICAL RESUL TS FOR THE RECIRCULATING FLOW INSIDE A SQUARE CA VITY Let n the unit square JO, 11 x JO , 11 . We assume that the wall y = 1 is moving with a relative velocity g = - 1 . We shall apply the method for k= 2. We use regular meshes that allow easy comparison with finite-difference methods . We approximate >,.h by continuous functions which are pieced together from second degree polynomials defined over the triangles and approximate CTlj by linear functions over triangles . 01 º' STOKES RE = 100 Fig. 2. Streamlines in the square cavity. STOKES RE= 100 _, RE = 5000 R = 10000 Fig. 3. Lines of constant vorticity in the square cavity. 29 We have fifteen degrees of freedom. six corresponding to Ah and nine to <Th. We calculate an approximation (>. , <1) of (Ah, <Th). We recall that Ah is an approximate value of A the stream function; <Tn approximatly equal to dfj A, 1 ~ i, j ~ 2. With those we may approximate 6 Ah the vorticity; ti = curi A = (d2A, - d1A) the ve locity vector; D;¡ , 1 s i , j ~ 2 the viscous s tre ss tensor whe r e D 1 1 2 CT~2. D 22 = - 2 CT~2. D12 = D 21 = CT~2 - CT~¡ . Referent the numerical performance of the method, for a grid of 200 elements, h = 1 /1 O and 1681 unknowns, the region needed was 650 K, the time 31 seconds, 7 iterations and 1050 output-input . The computations were carricd out on a IBM 370/158 . RE = 5000 Remarks about the results Let A denote the stream function and CT·· = a2 ' 1 < 1• J. IJ tJ /\, - • s 2. Thc velocity vector u= (u 1, u2) is equal to curi A. Then u, = d2A and u2 = - d1A. For that CT22 (respecti vely -<T11) may be interpreted as the fluid acceleration in the horizontal direction (vertical), that is CT22 ( - CT 11 ) is the variation of the horizontal componen! (vertical) of the velocity vector. Yelocities on the centerlines for Reynolds numbers up 10000 are known (Cf. ( 14). ( 12)). For instan-ce, for Re=O the acceleration is positive on x = .5 for Os y s. s. negative at ali other points. The graph of <Tn for Re= O confirme this result (Fig . 5). Similarly the accelera-tion is zero on y = .5 at x = .2 1 and x= .79, positive over RE = 10000 Fig. 4. lso- C111 lines in the square cavity. 30 (0,.21), (.79,1) negative on the remaining interval and symmetric. The graph of <1 11 (Fig. 4) for Re=O confirms this results. On the otherhand , cJ 1J\ is zeroon x= .5 then cJ 21 ( .5.y)=O forall y in (0, l ), that is <1 12 is zero on x= .5 a fact confirmed by the graph (Fig. 6). But the graph is not too exact at the centre of primary wortex for Re=O (where the velocity is zero) because it is hard to approximate; for Re= 400 the non intersection of lines may be interpreted by the viscous nature of the primary vortex. The viscous stress tensor The tensor can be stated in terms of the stream function as Figure 6 shows the iso- <1 12 lines in the square cavity and Figure 7 shows the iso- (<122- <111 ) lines again in the square cavity. The upper right vortex Between Re = 1000 and Re= 2000 (exactly Re = 1500 ac-cording to ( 15)) there appears a third secondary vortex near the upper upstream comer. It size increases with the Rey-nolds numbers. The dimensions are given in Table 1. Table I TP CENTRE ANO DIMENSIONS OF SECONDARY VORTICES Reynolds Vortex LL h TP Centre Numbers b1 d1 1/ 10 o .50 .75 1110 100 .40 .75 . 14 .18 1/ 10 400 .41 .65 .30 .35 1/10 1000 .43 .57 .27 .29 1/10 2000 .45 .55 . 15 . 17 1/ 10 5000 .46 .53 1/ 12 10000 .47 .51 (REYNOLDS) 10000 ~ 5000 X X 1000 •X • 400 " X X 100 •X 4( o L----M-~-------4------it'---- X o RE o 100 400 1000 5000 10000 .5 •NALLASAMY- KRISHNA h = 1/50 X THISSTUDY h = 1/10 NALLASAMY THISSTUOY . 21 .79 .21 .79 . 17 .75 .19 . 76 .19 . 73 .09 . 76 . 11 . 78 . 10 .89 .os .88 .10 .89 TABLE 2 POINTS OF VANISHING VERTICAL ACCELERATION Vortex LR Vortex UR b2 d2 b3 d3 .10 .10 .14 .10 .23 .13 .33 .24 .04 .07 .28 .27 .04 . 10 . 10 . 11 .06 .13 ·Y 1 .5 X o~------------+-'"" O 100 400 1000 5000 10000 • NALLASAMY - KRISHNA h = 1/50 X TH IS STUOY h = 1/10 RE NALLASAMV THISSTUOV o 53 .SJ 100 •• o •oo -- 35 1000 25 26 SOJO -- 11 10000 os 10 TABLE 3 POINTS OF VANISHING HORIZONTAL ACCELERATION (REYNOLDS) 31 -2 o. o o. .5 STOKES RE= 100 - 2. o. o. RE= 400 RE= 1000 - 35. - 2 - 2. - 20. o RE= 5000 RE= 10000 Fig. S. lso- a 22 lines in the square cavity. 32 STOKES RE= 100 RE= 400 RE= 1000 RE = 5000 RE 10000 Fig. 6. lso-a12 lines in the square cavity. 33 + + + NALLASAMY et al. 119771 h • 1/50 000 HUGHES et 11. 11978) FORTIN et al 119791 h • 1/ 12 THIS STUOY h • 1/ 10 - 1 RE= 1000 STOKES RE= 1000 RE = 10000 Fig. 7. lso- D12 lines in the square cavity. NALL ASAMY et al. (19711 h • 1/50 000 HUGHES et al (19781 FORTIN e t al (19791 h 1/ 12 THIS STUOY h • 1/1 2 + o - 1 RE = 10000 Fig. 8. Velocity profiles on the centreline of the cavity. Fig. 9. Velocity profiles on the centreline of the cavity. 34 7. CONCLUDING REMARKS The results for the stream function and the vonicity agree with others obtained by different methods. lt should be noted that the vorticity is not an unknown. The utilization of upwind derivatives is effecti ve and the frontal method for the linear system solution is very good. The results fo r the v iscous streess tensor seems reali stic. The method is convenient in time. prccis ion, s ize of cent ral memory. Four unknowns are calculated dircctl y. Threc additional oncs may be approximatcd. REFERENCES ( 1) L. H cRMANN. 1-'inite eleme/11 bendint: a1111l.vsis Jór ¡,lates. J. Eng. M cch. Di v. ASCE. EMS. 93 ( 1967). (2) K . HELLAN. Analysis of e/astic plates i11fl11x11re by II sim¡>li-fiedfinite e/eme/11 method. A cla Poly1echnica Scanel inavica. Ci 46. Trondheim ( 1967). (3) W . V1ss1cR. A refined 111ixed t_\pe plate be11di11g t' le111e111. A .L A.A. Joumal 7 ( 1969). (4) C. JoHNS<>N. 0 11 the com·ergence ofa 111ixed fi 11ite ele111e111 method for plate be11di1111 problems, Numcr. Mali, . 2 1 ( 1973). (5) F. BREZ/1 anel P.A . RAvtART. Mixedfinite ele111e111111e1ltods for 4'• arder elliptic eq1101io11s. i11: J . Miller (ed .). Topic in N11111erica/ A11alysis 11/ (Academic Press. 1977). p¡>. 33-56. (6) J.L. L tONS Q11e/q11es méthodes de réso/ution des problemes aux limites non linéaires, (Dunod, Paris, 1968). (7) V . GtRAULT and P.A. RAVIART. Fi11i1e Ele111e111 Approxima-lion of the Navier-Stokes Equations. (Leciures No1cs in M athematics 749. Springer-Verlag. Berlín . 1979). (8) M . FoRTIN. Résol111io111111mériq11e des éq11a1io11s de Na"ier-S toke s par des é léme111s finis de t_vpe mixte, IRIA-LABORIA Tech. Report 184 ( 1976). (9) P. L ESAINT. Sur la resolutio11 des sysremes hyperboliques d11 premier ordre par des méthodes d ' éléme111s fi11is, Thcsc . Univcrsité ele Paris VI. (1975). ( 10) C . JOHNSON. A mixedfinite eleme111111ethod for the Nm•ier-Stokes eq11a1ions. RAIRO 12 ( 1978). ( 1 1) J. F1cuEROA. La méthode d' éléments fin is de H errma1111-Johmo11 por les éq11a1io11s de Navier-Srokes srationnaires. These ele 3cme cycle . Uni versi1é de Paris V I ( 1979). ( 12) M . FORTIN and F. THOMASSIT Mixedfinite eleme/11 method., for incompressihlejlow problems. J. ofComputational Ph\'-sic 31 ( / 979). · ( / ]) P. JoLY. La méthode fro111a/e, Univcrsite Pari~ V I Tcch. Report 771008 LA 189 ( 1977). ( 14) M . NALLASAMY and K . KRISHNA PRASAO. 011 cavitvf/011• 111 high Reynolds 111u11ben. J. Fluid M echanics 79 N" 2 ( J 977). pp. 39 1-4 14 . ( 15) M .D . ÜLSON and S1t111-Yu T uANN, Newfi11ite e/eme111 re-sults for the square ca1·ity. Computcrs and Fluids. Vol. 7 (1978). pp. 123- 135 . 35 

Mixed finite element solution for the navier-stokes equations

http://www.revistas.usach.cl/ojs/index.php/contribuciones/article/download/4241/26003308/

Mixed finite element solution for the navier-stokes equations

Abstract

Similar works

Full text

Available Versions

Contribuciones Científicas y Tecnológicas