Stabilization of the finite element method for the stationary Stokes problem for an incompressible fluid

Chapter 3 Applications and numerical examples This chapter presents an analysis of the results obtained solving the Stokes problem in both two and three dimensions. For each case two problems have been studied. The first one has an analytical solution in order to compare the accuracy of the numerical solution. The second one is the called cavity flow problem, a benchmark flow problem with determined boundary conditions. 3.1 Bidimensional case The application of the proposed finite element formulation in two dimensions implies certain type of elements, shape functions and Gauss points and weights. They are displayed in figures 3.1 to 3.4. In particular, the elements used in this work are quadrilateral and triangular as can be seen in figure 3.1. The node numeration for each type of element is shown at figure 3.2. The integrals are calculated using a Gauss quadrature. The Gauss points coordinates and their corresponding weights are displayed in figure 3.3. Finally, the shape functions are listed in figure 3.4. 17 Mixed elements Figure Name Velocity interpolation Pressure interpolation r r r r eQ1Q0 Linear Constant r r r r e e e e Q1Q1 Linear Linear r r r r rrr r r e e e e Q2Q1 Quadratic Linear r r r r rrr r r e e e e eee e e Q2Q2 Quadratic Quadratic  A A A r r r eP1P0 Linear Constant  A A A r r r e e e P1P1 Linear Linear  A A A r r r r e e e Mini Linear Linear  A A A rrr r r r e e e P2P1 Quadratic Linear  A A A rrr r r r eee e e e P2P2 Quadratic Quadratic rVelocity node ePressure node Figure 3.1: Types of mixed elements in 2 dimensions 18 Node numeration Type of element s s s s 1 2 34 Linear Quadrilateral s s s s s s s ss 1 2 34 5 6 7 8 9 Quadratic Quadrilateral       TTTTT T s s s 1 2 3 Linear Triangular       TTTTT T s s s s 1 2 3 4Linear Triangular (Mini)       TTTTT T s s ss s s 1 2 3 4 56 Quadratic Triangular Figure 3.2: Node numeration in 2D elements 19 Element type Approximation Gauss points Gauss weights Quadrilateral Linear (.577350,.577350) 1.000000 (.577350,-.577350) 1.000000 (-.577350,.577350) 1.000000 (-.577350,-.577350) 1.000000 Quadratic (-.774597,-.774597) -.308642 (0,-.774597) .493827 (.774597,-.774597) .308642 (-.774597,0) .493827 (0,0) .790123 (.774597,0) .493827 (-.774597,.774597) .308642 (0,.774597) .493827 (.774597,.774597) .308642 Triangular Linear (.500000,.500000) .166667 (0,.500000) .166667 (.500000,0) .166667 Linear (Mini) (.101286,.101286) .062970 (.101286,.797427) .062970 (.797427,.101286) .062970 (.470142,.470142) .066197 (.059716,.470142) .066197 (.470142,.059716) .066197 (.333333,.333333) .112500 Quadratic (.333333,.333333) -.281250 (.600000,.200000) .260417 (.200000,.600000) .260417 (.200000,.200000) .260417 Figure 3.3: Gauss Quadrature in 2 dimensions 20 Element type Approximation Shape functions Quadrilateral Constant N1=1 Linear N1=1 4(1 −ξ)(1 −η) N2=1 4(1 + ξ)(1 −η) N3=1 4(1 + ξ)(1 + η) N4=1 4(1 −ξ)(1 + η) Quadratic N1=1 4ξ(ξ−1)η(η−1) N2=1 4ξ(ξ+ 1)η(η−1) N3=1 4ξ(ξ+ 1)η(η+ 1) N4=1 4ξ(ξ−1)η(η+ 1) N5=1 2(1 −ξ2)η(η−1) N6=1 2ξ(ξ+ 1)(1 −η2) N7=1 2(1 −ξ2)η(η+ 1) N8=1 2ξ(ξ−1)(1 −η2) N9=(1 −ξ2)(1 −η2) Triangular Constant N1=1 Linear N1=ξ N2=η N3=1 −ξ−η Linear (Mini) N1=ξ N2=η N3=1 −ξ−η N4=27ξη(1 −ξ−η) Quadratic N1=ξ(2ξ−1) N2=η(2η−1) N3=(1 −2ξ−2η)(1 −ξ−η) N4=4ξη N5=4η(1 −ξ−η) N6=4ξ(1 −ξ−η) Figure 3.4: Shape functions in 2 dimensions 21 3.1.1 Problem with analytical solution There are several Stokes problems (2.14) that have an analytical solution. These problems allows the validation of the numerical methods used to approximate the solution. In order to obtain an analytical solution, the main issue is to find a velocity field such as ∇v= 0. A proposed velocity field that follows that condition in two dimensions can be written as, v1(x, y) = x2(1 −x)2(2y−6y2+ 4y3) v2(x, y) = −y2(1 −y)2(2x−6x2+ 4x3) In addition, the viscosity is taken as ν= 1. The pressure field can be chosen randomly because it is just needed its gradient. In this case, it has been chosen as, p(x, y) = x(1 −x) Then, according to equation (2.12), the body forces verify b1=−(24y−12)(x4−2x3+x2) + (4y3−6y2+ 2y)(12x2−12x+ 2) +1 −2x b2= (24x−12)(y4−2y3+y2) + (4x3−6x2+ 2x)(12y2−12y+ 2) The problem is solved in a square domain of unity sides, Ω = [0,1] ×[0,1]. This domain can be discretized in several ways depending on the finite elements used. In this work, the type of elements implemented are quadrilateral and triangular. Moreover, the discretization can be uniform or non-uniform distributed. Figure 3.5 shows an example of the different possible combinations implemented in this work. Figures 3.5(a) and 3.5(b) are a representation of 10x10 quadrilateral elements uniform and non-uniform distributed meshes. Similarly, figures 3.5(c) and 3.5(d) show 10x10 triangular uniform and non-uniform distributed meshes. It is important to note that for meshes of triangular elements, the elements located in the upper-right and lowerleft corners are modified in order to avoid the three nodes located on the boundary. 22 (a) Uniform quadrilateral (b) Non-uniform quadrilateral (c) Uniform triangular (d) Non-uniform triangular Figure 3.5: Type of 2D meshes 23 Comparison between stabilized and non-stabilized solution In this section it is presented a serie of graphics comparing the stabilized and non-stabilized solution of this problem with the analytical solution. In all cases, the stationary Stokes problem has been solved using a mesh with 20x20 elements with an uniform distribution. On one hand, figure 3.6 shows the velocity and pressure fields for the analytical and the numerical nonstabilized solution. Note that the numerical solution, for both the pressure and velocity fields, is unacceptable due to the instabilities appeared. On the other hand, figure 3.7 displays the comparison between the analytical and the numerical stabilized solution. In this case, smooth and stable velocity and pressure fields are obtained. Therefore, the numerical solution is acceptable. The solution of the problem using non-uniform distributed quadrilateral and triangular meshes leads to the same behavior as using uniform distributed meshes. In this particular case, the geometrical distribution of the elements does not affect the accuracy of the solution. 0 0.5 1 0 0.5 1 −5 0 5 x 1015 Numerical pressure field 0 0.5 1 0 0.5 1 0 0.2 0.4 Analytical pressure field (a) Pressure field 0 0.5 1 0 0.5 1Numerical velocity field 0 0.5 1 0 0.5 1Analytical velocity field (b) Velocity field Figure 3.6: Non-stabilized solution 24 0 0.5 1 0 0.5 1 0 0.2 0.4 Numerical pressure field 0 0.5 1 0 0.5 1 0 0.2 0.4 Analytical pressure field (a) Pressure field 0 0.5 1 0 0.5 1Numerical velocity field 0 0.5 1 0 0.5 1Analytical velocity field (b) Velocity field Figure 3.7: Stabilized solution Error depending on mesh size The calibration of the efficiency of this stabilization method is made measuring the error between the numerical approximation and the analytical solution. The error has been calculated as an absolute error: ||v|| =||vnum −vanalyt||∞(3.1) ||p|| =||pnum −panalyt||∞(3.2) Figure 3.8 shows the maximum absolute error value for both stabilized linear element approximations in an uniform distributed mesh. That is, the linear quadrilateral and the linear triangular. 25 3.2 Tridimensional case The numerical solution to the stationary Stokes problem in three dimensions follows the same patterns than in two dimensions. The type of elements used in this case are the generalization of the elements presented in two dimensions. That is, hexahedral and tetrahedral elements. The different kind of elements used for the tridimensional analysis are shown in figure 3.17. It is also important to define the node numeration in every type of element in order to implement the finite element method properly. The elements implemented in this work are: the bilinear hexahedral (8 nodes), two quadratic hexahedral (20 and 27 nodes) and a linear tetrahedral (4 nodes). The node numeration of every element type is presented in figure 3.18. The integrals are calculated using a Gauss quadrature as in the bidimensional case. The values of the Gauss points coordinates and their correspondent weights are shown in figure 3.19. Finally, the shape functions are presented in figure 3.20. 32 Mixed elements Figure Name Velocity interpolation Pressure interpolation      p p pp p p pp r r r r rr rr e e e e ee ee H1H1 Linear Linear      p p pp p p pp r r r r rr rr r r r r r r r r r r r r e e e e ee ee H2H1 (20 nodes) Quadratic Linear      p p pp p p pp r r r r rr rr r r r r r r r r r r r r e e e e ee ee e e e e e e e e e e e e H2H2 (20 nodes) Quadratic Quadratic      p p pp p p pp r r r r rr rr r r r r r r r r r r r r r r r r r r r e e e e ee ee H2H1 (27 nodes) Quadratic Linear      p p pp p p pp r r r r rr rr r r r r r r r r r r r r r r r r r r r e e e e ee ee e e e e e e e e e e e e e e e e e e eH2H2 (27 nodes) Quadratic Quadratic ppppppp HH   J J J r r r r e e e e T1T1 Linear Linear rVelocity node ePressure node Figure 3.17: Types of mixed elements in 3 dimensions 33 Node numeration Type of element       p p p p p pp p p p p p pppp r r r r r r r r 1 2 3 4 5 6 7 8 Linear Hexahedral       p p p p p pp p p p p p pppp r r r r r r r r r r r r r r r r r r r r 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1920 Quadratic Hexahedral (20 nodes)       p p p p p pp p p p p p pppp r r r r r r r r r r r r r r r r r r r r r r r r r r r 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1920 21 22 23 24 25 26 27 Quadratic Hexahedral (27 nodes) pppppppppp P P P P P P   @@@@@@ r r r r 1 2 3 4 Linear Tetrahedral Figure 3.18: Node numeration in 3D elements 34 Element type Approximation Gauss points Gauss weights Hexahedral Linear (-.577350,-.577350,-.577350) 1.000000 (.577350,-.577350,-.577350) 1.000000 (-.577350,.577350,-.577350) 1.000000 (.577350,.577350,-.577350) 1.000000 (-.577350,-.577350,.577350) 1.000000 (.577350,-.577350,.577350) 1.000000 (-.577350,.577350,.577350) 1.000000 (.577350,.577350,.577350) 1.000000 Quadratic (-.774597,-.774597,-.774597) .171468 (0,-.774597,-.774597) .274348 (.774597,-.774597,-.774597) .171468 (-.774597,0,-.774597) .274348 (0,0,-.774597) .438957 (.774597,0,-.774597) .274348 (-.774597,.774597,-.774597) .171468 (0,.774597,-.774597) .274348 (.774597,.774597,-.774597) .171468 (-.774597,-.774597,0) .274348 (0,-.774597,0) .438957 (.774597,-.774597,0) .274348 (-.774597,0,0) .438957 (0,0,0) .702332 (.774597,0,0) .438957 (-.774597,.774597,0) .274348 (0,.774597,0) .438957 (.774597,.774597,0) .274348 (-.774597,-.774597,.774597) .171468 (0,-.774597,.774597) .274348 (.774597,-.774597,.774597) .171468 (-.774597,0,.774597) .274348 (0,0,.774597) .438957 (.774597,0,.774597) .274348 (-.774597,.774597,.774597) .171468 (0,.774597,.774597) .274348 (.774597,.774597,.774597) .171468 Tetrahedral Linear (.585410,.138197,.138197) .015625 (.138197,.585410,.138197) .015625 (.138197,.138197,.585410) .015625 (.138197,.138197,.138197) .015625 Figure 3.19: Gauss Quadrature in 3 dimensions 35 Element type Approximation Shape functions Hexahedral Linear N1=1 8(1 −ξ)(1 −η)(1 −φ) N2=1 8(1 −ξ)(1 + η)(1 −φ) N3=1 8(1 + ξ)(1 + η)(1 −φ) N4=1 8(1 + ξ)(1 −η)(1 −φ) N5=1 8(1 −ξ)(1 −η)(1 + φ) N6=1 8(1 −ξ)(1 + η)(1 + φ) N7=1 8(1 + ξ)(1 + η)(1 + φ) N8=1 8(1 + ξ)(1 −η)(1 + φ) Quadratic (20 nodes) N1=1 8(1 −ξ)(1 −η)(1 −φ)(−ξ−η−φ−2) N2=1 8(1 −ξ)(1 + η)(1 −φ)(−ξ+η−φ−2) N3=1 8(1 + ξ)(1 + η)(1 −φ)(ξ+η−φ−2) N4=1 8(1 + ξ)(1 −η)(1 −φ)(ξ−η−φ−2) N5=1 8(1 −ξ)(1 −η)(1 + φ)(−ξ−η+φ−2) N6=1 8(1 −ξ)(1 + η)(1 + φ)(−ξ+η+φ−2) N7=1 8(1 + ξ)(1 + η)(1 + φ)(ξ+η+φ−2) N8=1 8(1 + ξ)(1 −η)(1 + φ)(ξ−η+φ−2) N9=1 4(1 −ξ2)(1 −η)(1 −φ) N10=1 4(1 −η2)(1 + ξ)(1 −φ) N11=1 4(1 −ξ2)(1 + η)(1 −φ) N12=1 4(1 −η2)(1 −ξ)(1 −φ) N13=1 4(1 −φ2)(1 −η)(1 −ξ) N14=1 4(1 −φ2)(1 −η)(1 + ξ) N15=1 4(1 −φ2)(1 + η)(1 + ξ) N16=1 4(1 −φ2)(1 + η)(1 + ξ) N17=1 4(1 −ξ2)(1 −η)(1 + φ) N18=1 4(1 −η2)(1 + ξ)(1 + φ) N19=1 4(1 −ξ2)(1 + η)(1 + φ) N20=1 4(1 −η2)(1 −ξ)(1 + φ) Quadratic (27 nodes) N1=1 8ξ(ξ−1)η(η−1)φ(φ−1) N2=1 8ξ(ξ+ 1)η(η−1)φ(φ−1) N3=1 8ξ(ξ+ 1)η(η+ 1)φ(φ−1) N4=1 8ξ(ξ−1)η(η+ 1)φ(φ−1) N5=1 8ξ(ξ−1)η(η−1)φ(φ+ 1) N6=1 8ξ(ξ+ 1)η(η−1)φ(φ+ 1) N7=1 8ξ(ξ+ 1)η(η+ 1)φ(φ+ 1) N8=1 8ξ(ξ−1)η(η+ 1)φ(φ+ 1) N9=1 4(1 −ξ2)η(η−1)φ(φ−1) N10=1 4(1 −η2)ξ(ξ+ 1)φ(φ−1) 36 Hexahedral Quadratic (27 nodes) N11=1 4(1 −ξ2)η(η+ 1)φ(φ−1) N12=1 4(1 −η2)ξ(ξ−1)φ(φ−1) N13=1 4(1 −ξ2)η(η−1)φ(φ+ 1) N14=1 4(1 −η2)ξ(ξ+ 1)φ(φ+ 1) N15=1 4(1 −ξ2)η(η+ 1)φ(φ+ 1) N16=1 4(1 −η2)ξ(ξ−1)φ(φ+ 1) N17=1 4ξ(ξ−1)η(η−1)(1 −φ2) N18=1 4ξ(ξ+ 1)η(η−1)(1 −φ2) N19=1 4ξ(ξ+ 1)η(η+ 1)(1 −φ2) N20=1 4ξ(ξ−1)η(η+ 1)(1 −φ2) N21=1 2(1 −ξ2)(1 −η2)φ(φ−1) N22=1 2(1 −ξ2)(1 −η2)φ(φ+ 1) N23=1 2(1 −ξ2)η(η−1)(1 −φ2) N24=1 2(1 −η2)ξ(ξ+ 1)(1 −φ2) N25=1 2(1 −ξ2)η(η+ 1)(1 −φ2) N26=1 2(1 −η2)ξ(ξ−1)(1 −φ2) N27=(1 −ξ2)(1 −η2)(1 −φ2) Tetrahedral Linear N1=1 N2=0 N3=0 N4=−1 Figure 3.20: Shape functions in 3 dimensions 37 3.2.1 Problem with analytical solution There are also several Stokes problems with analytical solution in three dimensions. The velocity field has also to verify the incompressibility condition, ∇v= 0. A suitable analytical velocity field can be found extending the analytical field presented for the bidimensional case to the three dimensional one. Hence, the analytical velocity field in three dimensions is chosen as, v1(x, y, z) = x2(1 −x)2(2y−6y2+ 4y3)(2z−6z2+ 4z3) v2(x, y, z) = y2(1 −y)2(2x−6x2+ 4x3)(2z−6z2+ 4z3) v3(x, y, z) = −2z2(1 −z)2(2x−6x2+ 4x3)(2y−6y2+ 4y3) The viscosity is taken as 1 whereas the pressure field has only to be continue because it is just needed its gradient, see (2.12). The pressure field chosen in this work is, p(x, y, z) = x2+y2+z2 Then, according to equation (2.12) the components of the body forces are, b1=−(x4−2x3+x2)(2y−6y2+ 4y3)(−12 + 24z) + (2z−6z2+ 4z3) (−12 + 24y)+ (2z−6z2+ 4z3)(2y−6y2+ 4y3)(2 −12x+ 12x2) + 2x b2=−(y4−2y3+y2)(2x−6x2+ 4x3)(−12 + 24z) + (2z−6z2+ 4z3) (−12 + 24x)+ (2z−6z2+ 4z3)(2x−6x2+ 4x3)(2 −12y+ 12y2) + 2y b3= 2(z4−2z3+z2)(2y−6y2+ 4y3)(−12 + 24x) + (2x−6x2+ 4x3) (−12 + 24y)+ 2(2x−6x2+ 4x3)(2y−6y2+ 4y3)(2 −12z+ 12z2) + 2z The computational domain in this case is a cube, Ω = [0,1] ×[0,1] ×[0,1]. Similar to the two dimensional example, there are several kinds of discretization for the three dimensional domain. Figure 3.21 shows four discretizations. Figures 3.21(a) and 3.21(b) show the meshes with hexahedral elements, uniform and non-uniform distributed, respectively. Similarly, figures 3.21(c) and 3.21(d) display two meshes composed by tetrahedral elements, with an uniform and a non-uniform distribution of the elements, respectively. In these cases, the mesh is composed by 5 elements in each dimension in order to see clearly the resulting mesh. Note that, as in the bidimensional case, the tetrahedral elements in the corners are modified to avoid elements with all nodes being belonging to the boundary. 38 (a) Uniform hexahedral (b) Non-uniform hexahedral (c) Uniform tetrahedral (d) Non-uniform tetrahedral Figure 3.21: Type of 3D meshes 39 Comparison between stabilized and non-stabilized solution In order to analyze the numerical solution obtained by the Galerkin Least- Squares method, figures from figure 3.22 to figure 3.25, show the computed velocity and pressure fields using a mesh composed by 10x10x10 linear hexahedral elements. The use of the stabilization method becomes critical in order to have an acceptable solution to the problem. In particular, figure 3.22(a) represents the horizontal components of the velocity field in the plane z = 0.3 of the solution obtained without implementing the stabilization method. Similarly, figure 3.22(b) shows the component vzof the velocity field in the same plane z = 0.3. These graphics clearly indicate that the obtained solution is unacceptable due to the instabilities appeared. On the contrary, figures 3.23(a) and 3.23(b) display the solution obtained implementing the stabilization method. In this case, the solution can be accepted in comparison to the analytical solution. The analysis of the pressure field is analog to the analysis of the velocity field. Figure 3.24 displays the pressure field in the plane z = 0.3 of the solution obtained without using the stabilization method. Similarly, figure 3.25 shows the pressure field in the plane z = 0.3 obtained with the implementation of the stabilization method. The differences between the behavior of the non-stabilized and the stabilized numerical solution can be clearly observed. The non-stabilized solution is unacceptable whereas the stabilized solution reflects the same behavior as the analytical solution. Both the non-stabilized and the stabilized solutions have been shown at plane z = 0.3 as an example. The behavior of the solutions in the other planes, from z = 0 to z = 1 is analog to the behavior in plane z = 0.3. Moreover, the solution obtained with meshes of tetrahedral elements follow the same pattern than the behavior shown in meshes with hexahedral elements. 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Numerical velocity field in z=3.000000e−001 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Analytical velocity field in z=3.000000e−001 (a) 00.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −0.5 0 0.5 1 Numerical velocity vz field in z=3.000000e−001 00.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −5 0 5 x 10−3 Analytical velocity vz field in z=3.000000e−001 (b) Figure 3.22: Velocity field of the non-stabilized solution 41 Viscosity ν= 0.1 C=1 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.7554 4.37·10−24.37·10−24.37·10−2 13x13x13 0.6539 3.31·10−23.31·10−23.31·10−2 15x15x15 0.5762 2.58·10−22.58·10−22.58·10−2 17x17x17 0.5149 2.06·10−22.06·10−22.06·10−2 19x19x19 0.4653 1.69·10−31.69·10−21.69·10−2 C=10 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.3471 1.15·10−21.15·10−21.16·10−2 13x13x13 0.2970 8.59·10−38.59·10−38.64·10−3 15x15x15 0.2595 6.58·10−36.58·10−36.61·10−3 17x17x17 0.2304 5.21·10−35.21·10−35.23·10−3 19x19x19 0.2071 4.23·10−34.23·10−34.24·10−3 C=100 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2401 1.61·10−31.61·10−31.70·10−3 13x13x13 0.2051 1.35·10−31.35·10−31.40·10−3 15x15x15 0.1790 1.06·10−31.06·10−31.09·10−3 17x17x17 0.1587 8.59·10−48.59·10−48.73·10−4 19x19x19 0.1426 7.14·10−47.14·10−47.23·10−4 C=1000 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2194 3.40·10−43.40·10−43.59·10−4 13x13x13 0.1877 2.75·10−42.75·10−43.04·10−4 15x15x15 0.1639 2.20·10−42.20·10−42.54·10−4 17x17x17 0.1455 1.93·10−41.93·10−42.17·10−4 19x19x19 0.1308 1.73·10−41.73·10−41.95·10−4 C=10000 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2164 3.59·10−43.59·10−43.74·10−4 13x13x13 0.1851 2.83·10−42.83·10−43.01·10−4 15x15x15 0.1617 2.29·10−42.29·10−42.36·10−4 17x17x17 0.1435 1.97·10−41.97·10−42.08·10−4 19x19x19 0.1290 1.76·10−41.76·10−41.87·10−4 Figure 3.28: Errors for viscosity ν= 0.1 48 Viscosity ν= 0.01 C=1 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.7554 4.38·10−14.38·10−14.38·10−1 13x13x13 0.6539 3.30·10−13.30·10−13.30·10−1 15x15x15 0.5762 2.57·10−12.57·10−12.57·10−1 17x17x17 0.5149 2.05·10−12.05·10−12.05·10−1 19x19x19 0.4653 1.68·10−11.68·10−11.68·10−1 C=10 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.3471 1.16·10−11.16·10−11.16·10−1 13x13x13 0.2970 8.51·10−28.51·10−28.51·10−2 15x15x15 0.2595 6.49·10−26.49·10−26.50·10−2 17x17x17 0.2304 5.12·10−25.12·10−25.12·10−2 19x19x19 0.2071 4.15·10−24.15·10−24.15·10−2 C=100 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2401 1.73·10−21.73·10−21.74·10−2 13x13x13 0.2051 1.28·10−21.28·10−21.29·10−2 15x15x15 0.1790 9.79·10−39.79·10−39.82·10−3 17x17x17 0.1587 7.73·10−37.73·10−37.75·10−3 19x19x19 0.1426 6.27·10−36.27·10−36.28·10−3 C=1000 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2194 1.78·10−31.78·10−31.91·10−3 13x13x13 0.1877 1.48·10−31.48·10−31.55·10−3 15x15x15 0.1639 1.16·10−31.16·10−31.20·10−3 17x17x17 0.1455 9.37·10−49.37·10−49.58·10−4 19x19x19 0.1308 7.77·10−47.77·10−47.89·10−4 C=10000 Elements Error P Error Vx Error Vy Error Vz 11x11x11 0.2164 3.72·10−43.72·10−43.75·10−4 13x13x13 0.1851 3.01·10−43.01·10−43.13·10−4 15x15x15 0.1617 2.44·10−42.44·10−42.62·10−4 17x17x17 0.1435 2.08·10−42.08·10−42.22·10−4 19x19x19 0.1290 1.87·10−41.87·10−41.98·10−4 Figure 3.29: Errors for viscosity ν= 0.01 49 100101102103104 10−4 10−3 10−2 10−1 100 Constant Error Vx/Vy error vs Constant Viscosity 10 Viscosity 1 Viscosity 0.1 Viscosity 0.01 (a) 100101102103104 10−4 10−3 10−2 10−1 100 Constant Error Vz error vs Constant Viscosity 10 Viscosity 1 Viscosity 0.1 Viscosity 0.01 (b) 100101102103104 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Constant Error Pressure error vs Constant Viscosity 10 Viscosity 1 Viscosity 0.1 Viscosity 0.01 (c) Figure 3.30: Error vs. Constant depending on the fluid viscosity. 50 3.2.2 Cavity flow problem The definition of the Cavity flow problem in three dimensions can be easily derived from the two dimensional case. The domain of the problem is a cubic cavity with unity sides, Ω = [0,1] ×[0,1] ×[0,1]. The upper plane of the cavity moves at unit speed while the other ones are fixed. These boundary conditions can be seen in figure 3.31:       p p p p p pp p p p p p pppp - - - - - - - - - vx= 1 vy= 0 vz= 0 Figure 3.31: Boundary conditions for the Cavity flow problem in 3D Figures from 3.32 to 3.34 show the numerical solution to the cavity flow problem in three dimensions in a mesh of 20x20x20 linear hexahedral elements with a non-uniform distribution in various planes z. Figures 3.32(a) to 3.34(a) represent the horizontal components of the velocity field. Figures 3.32(b) to 3.34(b) show the z-component of the velocity field. Note that the numerical solution obtained reproduces the boundary conditions of the cavity flow problem, see figures 3.34(a) and 3.34(b). Figures 3.32(c) to 3.34(c) display the pressure field. The behavior of the obtained solution is analog to the two dimensional case. The use of linear elements without the stabilization method leads to an unstable solution. In contrary, the implementation of the Galerkin Least-Squares method enables the use of linear elements obtaining an stable numerical solution. 51 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Numerical horizontal velocity field in the plane z=1.781169e−001 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 Numerical Vz component of the velocity field in the plane z=1.781169e−001 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Numerical pressure field in the plane z=1.781169e−001 (c) Figure 3.32: Stabilized solution of the Cavity flow problem in z= 1.78 ·10−1 52 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Numerical horizontal velocity field in the plane z=8.218831e−001 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 Numerical Vz component of the velocity field in the plane z=8.218831e−001 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −10 −5 0 5 10 Numerical pressure field in the plane z=8.218831e−001 (c) Figure 3.33: Stabilized solution of the Cavity flow problem in z= 8.22 ·10−1 53 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Numerical velocity field in z=1 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 Numerical velocity vz field in z=1 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 −300 −200 −100 0 100 200 300 Numerical pressure field in z=1 (c) Figure 3.34: Stabilized solution of the Cavity flow problem in z= 1.00 54