Amr_guaily--directionally Adaptive Least Squares Finite Element Method

DIRECTIONALLY ADAPTIVE LEAST SQUARES FINITE ELEMENT METHOD FOR THE COMPRESSIBLE EULER EQUATIONS By

Amr Gamal Mohammad Guaily B.Sc. in Aerospace Engineering, 2002 A Thesis Submitted to the Faculty of Engineering, Cairo University in Partial Fulfillment of the Requirements for the Degree of Master of Science in Engineering Mechanics

FACULTY OF ENGINEERING, CAIRO UNIVERSITY GIZA, EGYPT 2006

ACKNOWLEDGMENTS All gratitude is due to Allah almighty. The author wishes to express his gratitude to all those who provided help in various ways at the different stages of this work. I wish to express my deepest and sincere gratitude and appreciation to my main supervisor, prof. Dr. A. A. Megahed, professor of engineering mechanics, for consideration, for suggesting the problem, and for his sincere guidance during this work. I also wish to express my gratefulness to my supervisor Prof. Dr. M. M. Abd-El-Rahmann, professor of aerodynamics, for his continuous support and guidance during this work. Assist. Prof. M. W. El-Mallah has played an important role in this work through his valuable discussions that have been very useful in overcoming technical difficulties encountered during the work, so I wish to express my deepest and sincere gratitude to him. My family has always played an important role in my studies, so I would like to express my deep gratitude and appreciation to my parents and my elder brothers for their continuous support.

i

ABSTRACT The least-squares finite element method is used to solve the compressible Euler equations in both 2-D Cartesian and axisymmetric forms.

Since the

method is naturally diffusive, no explicit artificial viscosity is added to the formulation. The inherent artificial viscosity, however, is usually large and hence does not allow sharp resolution of discontinuities unless extremely fine grids are used.

To remedy this problem, while retaining the advantages of the least

squares method, a moving-node grid adaptation technique is used. The outstanding feature of the adaptive method is its sensitivity to directional features like shock waves, leading to the automatic construction of adapted grids where the element edge(s) are strongly aligned with such flow phenomena. Using well-known transonic and supersonic test cases, it is demonstrated that by coupling the least squares method with a robust adaptive method, shocks can be captured with high resolution despite using relatively coarse grids. A paper extracted from the thesis was accepted to be presented at the IASTED international conference on modeling and simulation (MS 2006), which will be held May 24 to May 26, 2006, at Montreal, Canada. Paper title “Enhanced Adaptive Finite Element Method for The Cartesian and Axisymmetric Inviscid Compressible Flows”

ii

Table of Contents ACKNOWLEDGMENTS

i

ABSTRACT

ii

Table of Contents

iii

List of Figures

v

Nomenclature

x

Chapter 1

1

Introduction and Literature Review

1.1

Introduction

1

1.2

The Finite Element Method

1

1.3

Advantages of The Finite Element Method

2

1.4

Approaches of the FEM Formulation

3

1.4.1

Direct Approach

3

1.4.2

Variational Approach

3

1.4.3

Weighted residual Approach

3

1.5

Literature Review of the Compressible Euler FEM Schemes 3

1.6

Literature Review of the Adaptation Techniques

5

1.7

Current Work

8

Chapter 2

Least Squares FEM for Inviscid Compressible Flows

10

2.1

Introduction

10

2.2

Least Squares Formulation

10

2.3

Interpretation of the Inherent Viscosity

14

2.4

Finite Element Approximation

15

2.5

Boundary Conditions

17

2.6

Solution Method

19

iii

2.7

Numerical Results

2.7.1

20

Planer Problems

20

2.7.1.1 Shock-Reflection Problem (SRP)

20

2.7.1.2

Supersonic Channel Flow (SCP)

24

2.7.1.3

Airfoil-Spoiler Configuration

33

2.7.2

Axisymmetric Problems

39

2.7.2.1

Nozzle Flow (Shock Free)

39

2.7.2.2

Nozzle Flow (With Shock)

45

2.7.2.3

Jet Flow (Nearly Fully Expanded)

49

2.7.2.4

Jet Flow (Fully Expanded)

53

2.7.2.5

Under Expanded Jet Flow (UEJF)

54

Chapter 3

Directionally Adaptive Technique for FEM

60

3.1

Introduction

60

3.2

Mathematical Analysis

60

3.2.1

Edge-Based Error Estimate

60

3.2.2

Moving-Node Scheme

64

3.2.3

Grid Smoothening

66

3.2.4

The Grid Adaptation Procedure

67

3.3

Numerical Results

68

3.3.1

Analytical Test Case

68

3.3.2

Shock Reflection Problem (SRP)

71

3.3.3

Supersonic Channel Problem (SCP)

94

Chapter 4

Summary and Conclusions

100

4.1

Thesis Summary

100

4.2

Conclusions

101

4.3

Recommendations for future work

103

References

104

iv

List of Figures Figure 1-1

Original mesh to be adapted

6

Figure 1-2

r-method used in adaptation

6

Figure 1-3

h-method used in adaptation

6

Figure 2-1

Element shape in global and local coordinates

16

Figure 2-2

Discontinuous angle at wall nodes

18

Figure 2-3

Computational domain (SRP)

20

Figure 2-4

Pressure contours for different time steps, (SRP)

21

Figure 2-5

Pressure distribution at y= 0.05 (SRP)

22

Figure 2-6

Density contours for different time steps, (SRP)

22

Figure 2-7

Mach number contours for different time steps, (SRP)

23

Figure 2-8

Convergence history of the flow solver (SRP)

23

Figure 2-9

Computational domain (SCP)

24

Figure 2-10

Computational domain and grid (SCP)

26

Figure 2-11

Density contours, ∆t=0.15 (SCP)

26

Figure 2-12

Density contours, ∆t=0.1 (SCP)

27

Figure 2-13

Density contours, ∆t=0.05 (SCP)

27

Figure 2-14

Pressure contours, ∆t=0.15 (SCP)

28

Figure 2-15

Pressure contours, ∆t=0.1 (SCP)

28

Figure 2-16

Pressure contours, ∆t=0.05 (SCP)

29

Figure 2-17

Mach number contours, ∆t=0.15 (SCP)

29

Figure 2-18

Mach number contours, ∆t=0.1 (SCP)

30

Figure 2-19

Mach number contours, ∆t=0.05 (SCP)

30

Figure 2-20

Mach number distribution for the upper wall (SCP)

31

Figure 2-21

Mach number distribution for the lower wall (SCP)

31

Figure 2-22

Convergence history of the flow solver (SCP)

32

Figure 2-23

Velocity vector plot (SCP)

32

Figure 2-24

Zoom on the lower boundary (SCP)

33

Figure 2-25

Oil drop flow visualization

34

Figure 2-26

Airfoil-spoiler details

35

v

Figure 2-27

Grid used in airfoil-spoiler configuration

35

Figure 2-28

Zoom on the grid near the airfoil

36

Figure 2-29

Pressure coefficient contours (airfoil-spoiler problem)

36

Figure 2-30

Coefficient of pressure on the upper surface

37

Figure 2-31

Coefficient of pressure on the lower surface

37

Figure 2-32

Growth of wake in (airfoil-spoiler problem)

38

Figure 2-33

Physical nozzle and the used grid (shock free nozzle flow)

40

Figure 2-34

Mach number contours (shock free nozzle flow)

40

Figure 2-35

Pressure contours (shock free nozzle flow)

41

Figure 2-36

Radial velocity contours (shock free nozzle flow)

41

Figure 2-37

Pressure distribution (shock free nozzle flow)

42

Figure 2-38

Mach number distribution (shock free nozzle flow)

43

Figure 2-39

Temperature distribution (shock free nozzle flow)

43

Figure 2-40

Density distribution (shock free nozzle flow)

44

Figure 2-41

Velocity vector plot (shock free nozzle flow)

44

Figure 2-42

Zoom on the lower boundary (shock free nozzle flow)

45

Figure 2-43

Mach number contours (nozzle flow with shock)

46

Figure 2-44

Mach number distribution (Nozzle Flow With Shock)

46

Figure 2-45

Pressure contours (Nozzle Flow With Shock)

47

Figure 2-46

Pressure distribution (Nozzle Flow With Shock)

47

Figure 2-47

Convergence history (Nozzle Flow With Shock)

48

Figure 2-48

Velocity vector plot (Nozzle Flow With Shock)

48

Figure 2-49

Supersonic free jet main regions

50

Figure 2-50

Axisymmetric jet flow boundaries

50

Figure 2-51

Computational domain and center line Mach number

distribution (jet flow problem)

51

Figure 2-52

Axial Mach number distributions (jet flow problem)

52

Figure 2-53

Mach number contours (jet flow problem)

52

Figure 2-54

Mach Number Contours (Jet Flow Full Expanded)

53

Figure 2-55

Radial velocity Contours (Jet Flow Full Expanded)

54

Figure 2-56

Computational domain (UEJF)

55

Figure 2-57

Contour plot for the Mach number (UEJF)

56

vi

Figure 2-58

Axial Mach number distribution (UEJF)

56

Figure 2-59

Axial pressure distribution (UEJF)

57

Figure 2-60

Axial density distribution (UEJF)

57

Figure 2-61

Axial distribution of radial-velocity (UEJF)

58

Figure 2-62

Radial distribution of axial Mach number (UEJF)

58

Figure 2-63

Radial distribution of pressure (UEJF)

59

Figure 2-64

Radial distribution of radial-velocity (UEJF)

59

Figure 3-1

Transformation of a unite circle to an ellipse by S

63

Figure 3-2

Spring analogy for a patch of elements.

64

Figure 3-3

A node on the grid and its surroundings

66

Figure 3-4

Initial mesh (analytical test case).

69

Figure 3-5

Isocontours for the initial mesh (analytical test case).

69

Figure 3-6

Adapted mesh (analytical test case).

69

Figure 3-7

Adapted isocontours (analytical test case).

70

Figure 3-8

Magnification of grid (analytical test case).

70

Figure 3-9

The computation domain (SRP).

71

Figure 3-10

The initial and 1st adapted pressure (SRP)

72

nd

rd

Figure 3-11

The 2 and 3 adapted pressure (SRP)

72

Figure 3-12

The initial grid and the and the pressure contours (SRP)

74

Figure 3-13

The 1st adapted grid and the pressure contours (SRP)

74

Figure 3-14

The 2nd adapted grid and the pressure contours (SRP)

75

Figure 3-15

The 3rd adapted grid and the pressure contours (SRP)

75

Figure 3-16

The 4th adapted grid and the pressure contours (SRP)

76

th

Figure 3-17

The 5 adapted grid and the pressure contours (SRP)

76

Figure 3-18

Zoom on the grid in the incident shock region SRP

77

Figure 3-19

Zoom on the grid in the reflected shock region (SRP)

77

Figure 3-20

Pressure distribution at y=0.5 (SRP)

78

Figure 3-21

Convergence history (SRP)

78

Figure 3-22

The initial grid and the Mach number contours (SRP)

79

st

Figure 3-23

The 1 adapted grid and Mach number SRP

79

Figure 3-24

The 2nd adapted grid and Mach number (SRP)

80

Figure 3-25

The 3rd adapted grid and Mach number (SRP)

80

vii

Figure 3-26

The 4th adapted grid and Mach number (SRP)

81

Figure 3-27

The 5th adapted grid and Mach number (SRP)

81

Figure 3-28

The initial grid and pressure (density adapted) SRP

82

Figure 3-29

The 1st adapted grid and pressure (density adapted) SRP

82

Figure 3-30

the 2nd adapted grid and pressure (density adapted) SRP

83

Figure 3-31

The 3rd adapted grid and pressure (density adapted) SRP

83

th

Figure 3-32

The 4 adapted grid and pressure (density adapted) SRP

84

Figure 3-33

The 5th adapted grid and pressure (density adapted) SRP

84

Figure 3-34

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 1st adp. SRP Figure 3-35

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 2nd adp. SRP Figure 3-36

88

Comparison of pressure adapted (upper) and density adapted

(lower) grid for the 2nd adp. SRP Figure 3-41

88

Comparison of pressure adapted (upper) and density adapted

(lower) grid for the 3rd adp. SRP Figure 3-42

89

Comparison of pressure adapted (upper) and density adapted

(lower) grid for the 4th adp. SRP Figure 3-43

87

Comparison of pressure adapted (upper) and density adapted

(lower) grid for the 1st adp. SRP Figure 3-40

86

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 5th adp. SRP Figure 3-39

86

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 4th adp. SRP Figure 3-38

85

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 3rd adp. SRP Figure 3-37

85

89

Comparison of pressure adapted (upper) and density adapted

(lower) grid for the 5th adp. SRP

90

Figure 3-44

Comparison of pressure distribution at y=0.5 for SRP

Figure 3-45

Comparison of pressure adapted (upper) and density adapted

(lower) Mach number contours for the 1st adp. SRP

viii

90 91

Figure 3-46

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 2nd adp. SRP Figure 3-47

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 3rd adp. SRP Figure 3-48

92

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 4th adp. SRP Figure 3-49

91

92

Comparison of pressure adapted (upper) and density adapted

(lower) pressure contours for the 5th adp. SRP

93

Figure 3-50

Comparison of pressure contours (SRP)

93

Figure 3-51

The initial grid and the corresponding pressure contours

95

Figure 3-52

The 1st adapted grid and the corresponding pressure cont. 95

Figure 3-53

The 2nd adapted grid and the corresponding pressure cont. 96

Figure 3-54

The 3rd adapted grid and the corresponding pressure cont. 96

Figure 3-55

The 4th adapted grid and the corresponding pressure cont. 97

Figure 3-56

The 5th adapted grid and the corresponding pressure cont. 97

Figure 3-57

Evolution of density contours during adaptation

98

Figure 3-58

Magnification of grid in reflected shock region

99

Figure 3-59

Convergence history

99

ix

Nomenclature A x , A y , Aaxis

Jacobian matrices

E

Nodal error

g

Adaptation parameter

H

Hessian matrix

h

The element edge length

I

Functional to be minimized

k

Spring constant

L2

Residual error

P

Potential energy

Q

Vector of unknowns

ℜ

Residual vector (scalar)

T

Transformation matrix

t

Time

W

Weight function

x, y

Cartesian coordinates

Ω

Solution domain

φ

Basis function

δ

Dirac delta function

ζ

Local element coordinate

Γ

Contours of the elements

Λ

Diagonal matrix of the eigen-values

ω

Relaxation parameter.

x

Subscripts i, j

Nodal indices

∞

Free stream value

dim

Dimensional value

t, n

Normal and tangential coordinates

Superscripts T

Transposes

n

Time level

e

Element

nd

Number of nodes per element

m

Level of iteration

xi

Chapter 1

Introduction and Literature Review

1.1 Introduction This research aims at building a least squares finite element Euler solver for the transonic regime combined with an error estimate and a directionally adaptive grid algorithm, allowing the use of anisotropic (stretched) elements. Adaptive finite element methods place more fine scale elements where more resolution is needed. Isotropic or shape regular adaptive methods use only elements with bounded aspect ratio (stretched elements are avoided). Anisotropic adaptive methods fit high aspect ratio elements (highly stretched elements) along the regions of rapid variation of the solution for situations like shocks or boundary layers. Anisotropic adaptive methods give a bigger saving in terms of computational cost (number of elements and degrees of freedom) than the isotropic ones if stretched elements are placed appropriately [1]. The current interest in the area of high-speed flows has increased the need for advanced computational fluid dynamics (CFD) codes, which have become the primary tools for the prediction of aero-thermal loads. Such flows are characterized by regions with steep directional gradients of flow variables, embedded in regions where the flow variables vary more smoothly. One approach [2] for improving the solution accuracy of such problems is to apply grid adaptation techniques.

1.2 The Finite Element Method It was originally introduced by civil engineers: Harold Martin at the university of Washington (also a Structural Analysis Consultant to the Boeing Company), and H. Argyris working at Stuttgart and at Imperial Collage.

1

The essence of the method was certainly available for a long time, but the above two names made it a practical proposition. In addition both names were among the first to suggest and actually use the method for fluid dynamics problems. The FEM is essentially a numerical method for solving Partial Differential Equations (PDE's) [4]. In general, if an analytical (exact) solution cannot be found for the PDE’s governing a given continuum problem, then a numerical solution must be attempted at a discrete number of points. The continuum field in fluid mechanics is hence subdivided into a grid at which the dependent variable(s) is (are) to be calculated. This dependent variable(s) is (are) assumed to behave locally in a given way (linearly, quadratically, according to a spline, etc.). In the finite difference method (FDM) the PDE is written in discrete form involving the grid points. For the FEM the PDE is first recast in an integral form over the entire domain and its residual minimized by several means. Both the FDM and the FEM result in a large system of simultaneous linear algebraic equations. The dependent variable is solved for at all grid points simultaneously through direct or iterative methods.

1.3 Advantages of The Finite Element Method The FEM approximates complicated geometrical boundaries easily. Almost invariably the FDM starts by regularizing the domain, i.e. mapping complex regions into regular (rectangular) regions.

2

The FEM accounts for boundary conditions in an easy, straightforward manner. Unlike the FDM, the FEM has no numerical boundary conditions. The no penetration boundary condition is implemented by a Dirichlet condition even with a curved boundary.

1.4 Approaches of the FEM Formulation Given below are a different approaches of the FEM formulation. For more details see ref. [4] 1.4.1

Direct Approach This approach is applicable only to simple problems governed by an alge-

braic relationship or a simple first order ordinary differential equation (ODE) for 1-D problems. 1.4.2

Variational Approach This approach is applicable to physical problems governed by extremiza-

tion (minimum or maximum) laws. While this is mostly the case in structural analysis, it is not the general case for fluid mechanics where not many variational principles exist. 1.4.3

Weighted residual Approach This approach is applicable to problems where a variational principle does

not exist, for nonlinear problems and for unsteady state problems. The residual of the PDE is minimized, weighted by a certain weight that depends on the particular residual approach.

1.5 Literature Review of the Compressible Euler FEM Schemes So far, many different approaches have been adopted in developing numerical schemes to solve the compressible Euler equations. Using the idea of upwind in the finite difference method, Brooks and Hughes [5] introduced the Streamline Upwind Petrov-Galerkin (SUPG) method, in which the weight function is modi-

3

fied by adding a perturbation to the standard Galerkin test function. The added perturbation creates an upwind effect by weighting more heavily the upstream nodes. Another approach was proposed by Baruzzi et al.[6], where the Laplacian of dependent variables was added to the continuity and momentum equations. The amount of artificial viscosity was then controlled by a single parameter as the coefficient of the Laplacians. They later extended this first-order artificial viscosity method to second-order. Another method is based on the least squares weighed residual method. The method has very good stability properties due to its minimization nature, and has been applied for the solution of a variety of problems. As one of the earliest efforts in this field one can mention the technique presented by Polk and Lynn [7] for the solution of unsteady gas dynamic equations, with elements that are constructed in both space and time. Fletcher [8] used the least-squares method to solve the Euler equations for the subcritical compressible flows. The special feature of his method was to represent groups of variables rather than single variables. Application of the least squares method to a governing equation of the form: L(Q)=f leads to the favorable result of a symmetric and positive definite coefficient matrix, if L is a first-order differential operator. Taghaddosi et al. [2] applied this approach combined with adaptive grid algorithm to the Euler equations. Also Moussaoui [9] applied this approach to solve both the compressible and incompressible flows using one formulation. If L is a higher order operator, this property is completely lost during the integration by parts and moreover elements with higher order continuity requirements, must be employed. Lynn and Arya [10] proposed to break down the higher order system to its first order counterpart as a way of eliminating this disadvantage. Pontaza et al. [11] used this approach to solve both the Euler and Navier-Stokes equations for the compressible regime.

Bolton and Thatcher [12] used the

LSFEM to solve the Navier-Stokes equations in the form of stress and stream functions

4

1.6 Literature Review of the Adaptation Techniques For a given mesh as shown in Figure 1-1, grid adaptation methods are mainly composed of an adaptive strategy such as mesh movement (r-method) as in Figure 1-2; mesh refinement/coarsening (h-method); as in Figure 1-3, higherorder interpolation (p-method). Even when the same error estimate is used to assess the accuracy of the solution, the resulting adapted grid depends strongly on the selected adaptation strategy. Classical methods such as grid refinement produce isotropic meshes in which the length scales in all directions are essentially the same. These methods are optimal only for those flow regions possessing nearly equal gradients in all spatial directions. As a result, directional flow features such as shocks, contact discontinuities and boundary layers are not necessarily adapted efficiently and the number of elements needed to represent them may increase disproportionally with each isotropic refinement. An alternative approach would be to seek solutions on anisotropic meshes where more resolution is introduced along those directions with rapidly changing flow variables. This idea was introduced by Peraire et al. [13], who used an adaptive remeshing procedure that incorporated directional stretching for the solution of the 2D Euler’s equations on triangular grids. Another method may be used as proposed by Fortan et al. [14], who used a metric as a measure of error, coupled to an h-r strategy, to achieve directionally adapted grids with high aspect ratios.

5

Figure 1-1

Original mesh to be adapted

Figure 1-2

r-method used in adaptation

Figure 1-3

h-method used in adaptation

6

The above approaches have primarily been used on unstructured meshes. This trend is mainly driven by the intrinsic ability of triangular elements in 2D and tetrahedral elements in 3D to deal with arbitrarily complex geometries. Unstructured adaptation algorithms can yield highly stretched grids as well as locally refined/coarsened meshes. In contrast, most refinement techniques for structured grids avoid propagating the refinement to the boundaries by allowing sides to have hanging nodes. Despite these advantages for unstructured meshes, structured grids of quadrilateral elements in 2D and hexahedral elements in 3D are still used with great success in CFD. One reason is that the integration of the governing equations on a structured grid requires less CPU time than on an unstructured one for the same number of nodes. Structured grids are also more suitable for turbulence modeling, particularly near solid walls where normals to the wall may be needed. Furthermore, a certain degree of grid anisotropy may also be introduced for structured grids through an improved moving-node scheme. A moving-node technique was originally introduced by Gnoffo [15], and generalized by Nakahashi and Deiwert [16] in the context of finite volume methods (FVM). All these schemes are based on spring analogy where the grid is viewed as a network of springs whose stiffness constants represent a measure of error [15]. The grid vertices are displaced until the equilibrium state of the spring forces is reached. Such techniques are characterized by their low cost and the conservation of nodal connectivity, but often can stall or diverge and tolerate only a limited range of nodal movement. A directionally adaptive FEM using an edge-based error estimate on quadrilateral grids was proposed by Ait-Ali-Yahia et al.[18]. In this method the use of an appropriate error estimate, combined with the vector nature of spring forces (i.e. their magnitude and direction), permits one to design a convergent adaptive procedure capable of achieving wider nodal movement and a high degree of grid anisotropy. The error of the numerical solution is evaluated using a bound avail-

7

able from finite element interpolation theory. The Hessian of a selected solution variable is computed and then modified to produce a positive definite matrix allowing one to define a measure of error, namely a Riemannian metric. The edgebased error estimate is thus expressed as the length of the edge of the elements in this Riemannian metric. The construction of an anisotropic mesh may thus be interpreted as being a uniform mesh in the defined metric. This metric introduces and controls the magnitude as well as the direction of the grid anisotropy. A mesh movement scheme is then applied as the adaptive strategy, which in contrast to the spring analogy technique used by Nakahashi and Deiwert [16], has no constraint on grid orthogonality. This leads to a simple and efficient nodal redistribution algorithm offering a greater range of point displacements. In this method, the optimal grid for a fixed number of nodes is thus defined as one in which the error is equidistributed over the edges.

1.7 Current Work The least squares finite element method is used to solve the compressible Euler's Equations for both the Cartesian and axisymmetric flows in the nonconservative form. The quality of the numerical results indicates the remarkable performance of the used technique. This is quite evident from the final results despite using coarse mesh. Also the robustness of the technique allows one to use large time step to reach the steady state solution in few iterations, which saves computational time considerably. But a disadvantage existed; namely the high value of the inherent artificial viscosity in the method which forbidden the sharp resolution of discontinuities, this disadvantages has been remedied. To remedy the disadvantage mentioned above in the least-squares finite element method, an adaptive algorithm with directional features, using an edge based error estimate on quadrilateral meshes, is used. The error of the numerical solution is measured by its second derivatives and the resulting Hessian tensor is used to define a Riemannian metric. A mesh movement strategy with no or-

8

thogonality constraints is used to equidistribute the lengths of the edges of the elements of the edges of the elements in the defined metric. The adaptive procedure has been proven to be effective on analytical test case. The used technique enforces the solid wall boundary condition since it enforces the boundary condition as a Dirichlet type at each node using a transformation matrix for each node. This results in a faster convergence than other methods and it can be used for both steady and unsteady flows unlike other techniques. The flow solver, combined with the proposed grid adaptation method is then validated on a variety of problems. The quality of the numerical results indicates the remarkable performance of the adaptive method, and demonstrates its superiority to many existing techniques. This is quite evident in the final adapted grids. Overall, the combination of the least-squares method with the directionally adaptive method has produced promising results, despite using an artificial viscosity mechanism that has no free parameters, on relatively coarse grids. The rest of the thesis is structured as follows Chapter 2 deals with the least squares finite element formulation followed by a variety of internal and external flow problems to validate the method Chapter 3 deals with the formulation of the adaptation scheme followed by a variety of flow problems to validate the method. Chapter 4 concludes the results of preceding chapters and gives recommendations for future work.

9

Chapter 2

Least Squares FEM for Inviscid Compressible Flows

2.1 Introduction In this chapter the least-squares finite element method is introduced in details. The least squares finite element method is easy to implement, is naturally diffusive, contains no free parameter(s), is stable thus allowing equal-order interpolation of all variables, and more importantly, results in a symmetric positive definite coefficient matrix. Application of other finite element methods such as Galerkin or the SUPG leads to non-symmetric matrices. The present work, therefore, is an attempt to investigate the least squares finite element method for compressible flows with shocks.

2.2 Least Squares Formulation The Euler equations for both Cartesian and axisymmetric flows can be written in the nondimensional form as a first order system in terms of primitive variables [2] as:

∂Q ∂Q ∂Q + Ax + Ay + α AaxisQ = 0 ∂t ∂x ∂y where

Q T = ( ρ ,u ,v , p ) is the vector of unknowns.

10

(2-1)

u ρ 0  0 u 0 Ax =  0 0 u 0 γ p 0 

v 0  Ay =  0  0 

Aaxis

0 v

ρ

0

v

0

0 γp

v  y  0 =  0   0 

0  1  ρ  0  u 

    ρ v 

0 0 1

0   0 0 0   0 0 0   0 0 γv  y 0 0

with

α =1

for axisymmetric flow

α =0

for planar flow

in which ρ is the density, (u , v ) the velocity components, p the pressure, and

γ is the specific heat ratio. The non-dimensionlization is as follow, ρ = ρdim / ρ∞, u = udim / C∞, v = vdim / C∞, p = pdim / (ρ∞ C∞2), t = tdim C∞ / L. where (dim stands for dimensional) ρ

Density.

u

Longitudinal velocity.

v

Lateral velocity.

p

Pressure.

t

Time.

L

Characteristic length.

11

ρ∞

Free stream density.

C∞

Free stream speed of sound.

The system given in equation (2-1) may be linearized using Newton method, by setting:

Q n +1 = Q n + ∆Q

neglecting the higher order terms, and discratizing the unsteady term using backward difference as:

∂ Q ∂ t

∆ Q ∆ t

≈

so, equation (2-1) can be rewritten as:

L ∆Q n +1 = −f

(2-2)

where

L = A

f

f

= (A

ax is

∂ + A ∂x

n x

n x

∂Q ∂x

 ρ nv =  y 

n

∂ 1 + ( I + A ∂y ∆t

n y

+ A

n y

∂Q ∂y

n

0

0

12

n

+ α f

γ p nv y

ax is

n

  

T

n

)

)

∂ρ ∂x ∂u ∂x ∂v ∂x ∂p ∂x

v  ∂u ∂v ( ∂x + ∂y + α y )  −1 ∂p   ρ 2 ∂x  A = −1 ∂p   ρ 2 ∂y   0  

∂ρ ρ +α ∂y y ∂u ∂y ∂v ∂y ∂p γp +α ∂y y

0 0 0

γ(

∂u ∂v + ∂x ∂y

Defining the residual vector as

ℜ

=

L ∆ Q

n + 1

+ f

the least squares functional is

I (∆Q n +1 ) =

1 (ℜ)T (ℜ)d Ω ∫∫ Ω 2

(2-3)

Now we can introduce the finite element approximation.

ne

∆Q n +1 ≈ ∆Q hn +1 = ∑ N j ∆Q jn +1

(2-4)

j =1

where ne

is the number of nodes per element.

N j = N jI

is the element shape function.

and the weighting function is

W =

∂ℜ = LN ∂ (∆Q n +1 )

so the least squares weak form is

13

         v  + α ) y 

∫∫

Ω

(LN )T (L ∆Q n +1 + f )d Ω = 0

(2-5)

Using equation (2-4) in equation (2-5) results the linear algebraic equations:

[ K ] {∆ Q } =

−

{R }

(2-6)

where

K

e ij

=

∫∫ ( L N

i

)T ( L N

j

)d Ω

e

Ωe

rie = ∫∫ (L N i )T fd Ω

(2-7)

e

Ω

are evaluated using Gauss-Legendre quadrature.

2.3 Interpretation of the Inherent Viscosity For the Euler’s equation (2-1) the numerical viscosity inherent in the least squares formulation can be demonstrated as follows [2]. The right-hand side vector of equation (2-6) as defined in equation (2-7) is

ri e =

∫∫ ( L N

i

)T fd Ω

Ωe

  1 ∂N i ∂N i = − ∫∫  A xn + A yn +( I + A n )N i  ∂x ∂y ∆t  Ωe  n  n ∂Q n  n ∂Q ( ) + + α A A f  x y ax is  d Ω ∂x ∂y  

14

T

T

n   ∂N i ∂N i   n ∂Q n n ∂Q e A A = − ∫∫  A xn + A yn +  d Ω + ...   x y x y x y ∂ ∂ ∂ ∂    Ωe 

After integrating this equation by parts (Green’s formula), the first term on the right-hand side yields the following term T

n   n ∂ ∂   n ∂Q n n n ∂Q + + A A A A dΩ     x y x y ∫∫e  ∂ x ∂ ∂ ∂ y x y    Ω

which represents an inherent numerical viscosity of the least squares method. The effect of the time step on the amount of artificial viscosity is examined for each problem separately by numerical experiment to determine the appropriate time step.

2.4 Finite Element Approximation In the finite element discritization, both variables and problem geometry are approximated using shape functions. If the same shape functions are used for approximation of geometry and variables, the resulting element is called “isoparametric”. In this work iso-parametric quadrilateral bilinear elements will be used. This type of elements relies on the bilinear local shape function given by

N i (ξ , η ) = (1 + ξξ i )(1 + ηη i ) / 4

(2-8)

(ξi , ηi) are the local coordinate of the corner node “i” of any quadrilateral element, while (ξ , η) are the element local axes as illustrated in Figure 2-1. The integer index “i” varies between 1 and 4 according to the corresponding node.

15

Figure 2-1

Element shape in global and local coordinates

The problem variables are approximated by

∆ Q

4

= Σ ∆ Q iN i=1

(2-9)

i

Where ∆ Q i is the value of ∆ Q , at node “i”. Transformation between local coordinates and global Cartesian coordinates is governed by

4

x = Σ xi N i i =1

(2-10)

4

y = Σ yi N i i =1

Partial derivatives w.r.t. x and y are related to partial derivatives w.r.t. ξ and η by

16

 ∂Ni   ∂x  ∂ξ   ∂ξ  = N ∂ i    ∂x  ∂η   ∂η

∂y   ∂Ni   ∂Ni    ∂ξ   ∂x    ∂N  = [J ] ∂∂Nx  ∂y   i   i ∂η   ∂y   ∂y 

(2-11)

We note that the Jacobian matrix J is a function of ξ and η only. By inversion of J, derivatives w.r.t. x and y, become functions in ξ and η. Area integrals are transformed according to:

dx dy = dA = det(J) dξ dη.

(2-12)

The last relation equation (2-12) is very important since all calculus is done in terms of the local coordinates. Integration over each element is done numerically using Gauss-Legendre quadrature of the 5th order accuracy. This method of integration involves evaluation of the integrand at fixed points, and multiplication of the resulting value by a corresponding weight, followed by summation ([4], [17]).

2.5 Boundary Conditions In general we have four types of boundary conditions 1- Dirichlet type in which Q "the dependent variable" is specified. 2- Neumann type in which

∂Q is specified, n is the direction normal to the ∂n

boundary. 3- Cauchy type in which part of the boundary is Dirichlet type and the rest it Neumann type.

17

4- Robbins

type

in

which

c 1Q + c 2

∂Q = c 3 is ∂n

specified,

where,

c1 , c 2 , and c 3 are constants

For Euler’s equations, the number of boundary conditions to be imposed on the domain is determined by the theory of characteristics. In the finite element method, the numerical boundary conditions are naturally imposed which is considered one of the most important features of the finite element method. Therefore no special treatments are required. For our first order system we have Dirichlet type of boundary conditions except for the non-Dirichlet boundary condition at a solid wall U.n=0 where n is the outward unit normal.

Figure 2-2

Discontinuous angle at wall nodes

An interesting procedure is the one known as coordinate rotation method. In this formulation ([9],[19]), the velocity components at the solid wall are rotated from the global Cartesian coordinate to a local coordinate system which permits the tangency condition to be imposed as a Dirichlet boundary condition as follow. Let the normal direction at a solid wall point A be given by ( cosθ , sinθ ), and denote by un and ut the normal and tangential components of the velocity U, respectively. Now, we have the following relation:

18

 u   cos θ  = v   sin θ

sin θ   u n    − cos θ   ut 

For nodes on solid wall boundary, a new nodal vector of unknowns is defined by replacing the Cartesian components (u, v) by the normal and tangential components in the global matrix [ K ] . The resulting matrix is nonsymmetrical. To restore the symmetry we pre-multiply the corresponding rows by the matrix T defined as: 0 1  0 cos θ T=   0 sin θ  0 0

0 sin θ − cos θ 0

0  0 0  1

Due to the discrete representation of the computational domain, the angle at a typical wall node A (Figure 2-2) can have two different values. It is necessary to assign a unique angle for such nodes. To accomplish this, a weighted average of the angles from the two adjacent wall edges can be used:

θA =

L 2θ1 + L1θ 2 L1 + L 2

2.6 Solution Method The resulting system of linear algebraic equations (2.6) is solved as follows: we first compute a U DUT factorization of the coefficient matrix [ K ] . D is a block diagonal matrix with blocks of order 1 or 2, and U is a matrix composed of the product of a permutation matrix and a unit upper triangular matrix. The solution of the linear system is then found. The previous method is implemented using a built in subroutine in the used language (Compaq Visual FORTRAN 6.5).

19

2.7 Numerical Results 2.7.1

Planer Problems

2.7.1.1 Shock-Reflection Problem (SRP)

We consider compressible, inviscid flow along a wall where the reflection of an oblique shock occurs. The computational domain is a 4.1 X 1.0 rectangle, which is discretized uniformly into (60 X 20) bi-linear elements. We consider flow at an inlet Mach number of Min = 2.9 and specify the following boundary conditions at the inlet (ρ, u, v, p) = (1.0, 2.9, 0.0, 0.7143). To simulate the oblique shock we specify the following at the top boundary :( ρ, u, v, p) = (1.7, 2.6193, -0.5063, 1.5282) ([20], [21]). The lower boundary is impermeable and the exit boundary is left free (supersonic exit)

Figure 2-3

Computational domain (SRP)

The effect of the inherent viscosity that is controlled by the time step, on the shock resolution is shown in Figure 2-4, which shows the pressure contours, when the time step changes from 0.15 to 0.05 (from top to bottom ∆ t=0.15, 0.1, 0.05). As expected oscillations near the shock start growing after reducing the artificial viscosity (reducing the time step) and the shock becomes more smeared by increasing it (increasing the time step), as indicated in Figure 2-5. Figure 2-6

20

shows the density contours while Figure 2-7 shows the Mach number contours where the time step has changed from 0.15 to 0.05. This shows the effect of the time step on the solution as the shock becomes more refined for smaller time step. Figure 2-8 shows the convergence history using different time steps. As expected as the time step decreased the solution reaches steady state slower.

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1

0.5

0 1

0.5

0

Figure 2-4

Pressure contours for different time steps, (SRP)

21

Figure 2-5

Pressure distribution at y= 0.05 (SRP)

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1

0.5

0 1

0.5

0

Figure 2-6

Density contours for different time steps, (SRP)

22

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1

0.5

0 1

0.5

0

Figure 2-7

Figure 2-8

Mach number contours for different time steps, (SRP)

Convergence history of the flow solver (SRP)

23

2.7.1.2 Supersonic Channel Flow (SCP)

A more challenging problem for this algorithm exists in the analysis of supersonic flow past the parabolic arc bump (Figure 2-9) of height to length ratioα = 4%. This problem creates an interesting shock interaction behind the bump. The shock from the leading edge of the bump is reflected down off the ceiling of the tunnel and intersects with the shock formed at the trailing edge of the bump (see Figure 2-9 for details).

Figure 2-9

Computational domain (SCP)

The curved lower wall for the circular arc is given by

0   y =  R 2 − x 2 +b  0 

−1.5 < x < −0.5 −0.5 < x < 0.5 0.5 < x < 1.5

where

b= R=

(α

2

− 0.25 ) 2α

0.25 + b 2

24

The boundary conditions at the inlet are:

ρ =1.0, u=4 uinf y (1-y) , v=0.0 and at the exit, the pressure is specified as p = 1

γ

On the walls the flow tangency condition is imposed and the exit boundary is left free (supersonic exit). The grid consists of 64 X 16 uniformly distributed bi-linear rectangular elements with 16 elements in the y-direction, 22 elements on the bump, and 21 elements in each side of it as shown in Figure 2-10. The effect of the inherent viscosity, controlled by the time step, on the shock resolution is shown in the figures starting from Figure 2-11 to Figure 2-13, which show the density contours, calculated with time steps ranging from 0.15 to 0.05. Figure 2-14 to Figure 2-16 show the pressure contours, which assure the important role of the time step on shock resolution. The leading- and trailing edge shocks as well as the interaction of the trailing-edge shock with the reflected shock are qualitatively well captured. The Mach number contours are plotted in the figures starting from Figure 2-17 to Figure 2-19 to assure the ability of the scheme in capturing such a complex shock structure. Figure 2-20 and Figure 2-21 show the Mach number distribution of the upper and lower walls respectively compared to Eidelman et al. [22] which shows good agreement. As expected oscillations near the shock start growing after reducing the artificial viscosity and the shock becomes more smeared by increasing it.

25

Figure 2-22 shows the convergence history using different time steps. As expected as the time step decreased the solution reaches steady state slower. The quality of the used method in applying the solid wall boundary condition is demonstrated in the figures starting from Figure 2-23 to Figure 2-24, which show the velocity vector plot and a zoom on the lower wall. The tangency condition is applied very well as seen in these figures.

Figure 2-10

Computational domain and grid (SCP)

1

0.9 0.8 0.7 0.6

0.5 0.4 0.3 0.2 0.1 0 -1.5

-1

Figure 2-11

-0.5

0

0.5

1

Density contours, ∆t=0.15 (SCP)

26

1.5

1 0.9 0.8 0.7 0.6

0.5 0.4 0.3 0.2 0.1 0 -1.5

-1

Figure 2-12

-0.5

0

0.5

1

1.5

Density contours, ∆t=0.1 (SCP)

1 0.9 0.8 0.7 0.6

0.5 0.4 0.3

0.2 0.1 0 -1.5

-1

Figure 2-13

-0.5

0

0.5

1

Density contours, ∆t=0.05 (SCP)

27

1.5

1 0.9

0.8 0.7 0.6

0.5

0.4 0.3 0.2

0.1 0 -1.5

-1

Figure 2-14

-0.5

0

0.5

1

1.5

Pressure contours, ∆t=0.15 (SCP)

1 0.9 0.8 0.7 0.6

0.5 0.4 0.3 0.2 0.1 0 -1.5

-1

Figure 2-15

-0.5

0

0.5

1

Pressure contours, ∆t=0.1 (SCP)

28

1.5

1 0.9

0.8 0.7

0.6

0.5

0.4

0.3

0.2 0.1

0 -1.5

-1

Figure 2-16

-0.5

0

0.5

1

1.5

Pressure contours, ∆t=0.05 (SCP)

1 0.9

0.8 0.7

0.6

0.5

0.4

0.3

0.2 0.1

0 -1.5

-1

Figure 2-17

-0.5

0

0.5

1

Mach number contours, ∆t=0.15 (SCP)

29

1.5

1

0.9 0.8 0.7

0.6

0.5 0.4

0.3 0.2

0.1

0 -1.5

-1

Figure 2-18

-0.5

0

0.5

1

1.5

Mach number contours, ∆t=0.1 (SCP)

1 0.9

0.8 0.7 0.6

0.5

0.4 0.3 0.2 0.1

0 -1.5

-1

Figure 2-19

-0.5

0

0.5

1

Mach number contours, ∆t=0.05 (SCP)

30

1.5

Figure 2-20

Mach number distribution for the upper wall (SCP)

Figure 2-21

Mach number distribution for the lower wall (SCP)

31

Figure 2-22

Convergence history of the flow solver (SCP)

1 0.9

0.8 0.7

0.6

0.5

0.4

0.3 0.2

0.1 0 -1.5

-1

-0.5

Figure 2-23

0

0.5

Velocity vector plot (SCP)

32

1

1.5

Figure 2-24

Zoom on the lower boundary (SCP)

2.7.1.3 Airfoil-Spoiler Configuration

Spoilers are widely used in aircraft as lateral control devices, as speed breaks, and as lift dampers during landing. Despite their wide usage, very little theoretical information exists. Almost all design work is done by comparison with experimental results followed by wind tunnel tests of trail models. The flow field past an airfoil with spoiler is complex. In fact the flow separated from the upper airfoil surface, due to the adverse pressure gradient generated by the presence of the spoiler, then reattaches to the spoiler surface at small angles of attack [23]. A recirculating bubble is formed upstream of the spoiler hinge. The flow separates again from the spoiler tip and converts into wake as a free shear layer. The flow on the lower airfoil surface also leaves the trailing edge and converts into wake. These shedding vortices make the wake highly turbulent and oscillatory. This unsteady wake affects the effectiveness of the spoiler. Extensive experimental investigations on steady spoiler characteristics have been taken by several authors ([24], [25], [26],[27]). On the theoretical side, Abdelrahman [24] developed a numerical scheme based on the steady state incompressible Navier-Stokes equations. Choi et al. [28] studied the transient response of an airfoil to a rapidly deploying spoiler using turbulent compressible Navier-Stokes equations with a turbulence model.

33

Flow Description The flow was examined using both water table and oil streaks on a flowsplitter plate mounted on a wind tunnel model [23]. These examinations indicate that the wake has two main characteristics. First it has an unsteady nature due to a shed vortex street. Second, despite this periodic behavior there is a region of wake that remains nearly constant in shape and closes a short distance downstream of the trailing edge. Figure 2-25 represents the flow on the splitter plate as investigated by Wentz and Ostowari [25] at zero angle of attack, 40o splitter deflection angle, 2.2x106 Reynolds number, and 0.13 Mach number. It clearly shows two basic regions: first, an outer essentially potential flow region, and second, a near wake region behind the spoiler. Also Figure 2-25 shows that the near wake may be subdivided into two more parts. The part upstream of the wing trailing edge shows very little fluid motion and thus may be considered at constant pressure. The part downstream of the trailing edge is characterized by a pair of vortices; the upper one rotating clockwise (for left to right flow) and smaller one below it rotating counterclockwise

Figure 2-25

Oil drop flow visualization

34

Figure 2-26

Airfoil-spoiler details

Figure 2-26 describes the details of the used configuration, while Figure 2-27 shows the O-type grid of [65x20] elements used to solve the problem of the airfoil-spoiler configuration, followed be a zoom on the grid near the airfoil surface in Figure 2-28, where the parameters in Figure 2-26 are as follows for NACA 23012: d = 0.75 c, t= 0.01 c, δ = 30o, and L=0.15 c. The inlet Mach number used was M= 0.2.

Figure 2-27

Grid used in airfoil-spoiler configuration

35

Figure 2-28

Zoom on the grid near the airfoil

Figure 2-29 shows the pressure coefficient isocontours, while Figure 2-30 and Figure 2-31 show the distribution of the pressure coefficient over the upper and lower surfaces compared to the experimental results of [24]. The coefficient of pressure suddenly decreases to a negative value due to the presence of the spoiler.

Figure 2-29

Pressure coefficient contours (airfoil-spoiler problem)

36

Exp.

Num.

Cp

2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 0

0.2

0.4

0.6

0.8

1

X

Figure 2-30

Coefficient of pressure on the upper surface

Cp

Exp.

Num.

1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 0

0.2

0.4

0.6

0.8

X

Figure 2-31

Coefficient of pressure on the lower surface

37

1

Figure 2-32

Growth of wake in (airfoil-spoiler problem)

Figure 2-32 shows the growth of the wake behind the spoiler which in qualitative agreement with the experimental results of [25].

38

2.7.2

Axisymmetric Problems

2.7.2.1 Nozzle Flow (Shock Free)

Nozzle flow will be considered here as it is one of the most important test cases used to show code robustness. In particularly it tests the efficiency of the flow tangency condition, which is considered one of the most difficult aspects when talking about Euler's equations. Nozzle Geometry  ( 0.1398 + 0.0347 tanh ( −9x + 2.5 ) )2  A (x ) =   0.1398 + 0.0347 tanh ( 9 ( x − 1) + 2.5 )

(

0.0 ≤ x ≤ 0.5

)

2

0.5 ≤ x ≤ 1.0

where A is the nozzle cross-sectional area and the radius, r=(A / Π)0.5. The computational domain is discretized uniformly into [60 X 20] bi-linear elements as shown in Figure 2-33. We consider flow at an inlet Mach number of Min = 0.224 and specify the following boundary conditions at the inlet "subsonic inlet" (ρ, u, v) = (1.0, 0.224, 0.0). The lower and upper boundaries are impermeable and the exit boundary is left free (supersonic exit). Figure 2-34 shows the Mach number contours through the nozzle. Despite solving the whole domain, the full symmetry is evident from the Mach number contours as well as the pressure contours shown Figure 2-35. Figure 2-36 shows the contours of the radial velocity, which assures the capability of the used technique to capture the symmetry of the flow.

39

0.2

0.15

0.1

r

0.05

0

-0.05

-0.1

-0.15

-0.2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

Figure 2-33

Physical nozzle and the used grid (shock free nozzle flow)

0.15

0.1

0.05

0

-0.05

-0.1

-0.15 0

0.1

Figure 2-34

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Mach number contours (shock free nozzle flow)

40

1

0.15

0.1

0.05

0

-0.05

-0.1

-0.15 0

0.1

0.2

Figure 2-35

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Pressure contours (shock free nozzle flow)

0.15

0.1

0.05

0

-0.05

-0.1

-0.15 0

0.1

Figure 2-36

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Radial velocity contours (shock free nozzle flow)

41

1

The flow accelerates (Mach number increases) as it goes through the nozzle as shown in Figure 2-37 which presents the pressure distribution of the nozzle center line as well as the nozzle wall compared to the exact one dimensional solution [29]. Figure 2-38 shows the Mach number distribution for both the nozzle wall and the centerline compared to the one-dimensional solution. The wall effect is evident in the distribution of the wall Mach number as it starts to decrease at an axial location of 0.7 due to the change in the area (acts like a compression corner). Figure 2-39 shows the distribution of the temperature across the nozzle in which again the two dimensional effect is clear. The density distribution is shown in Figure 2-40.

Figure 2-37

Pressure distribution (shock free nozzle flow)

42

Figure 2-38

Mach number distribution (shock free nozzle flow)

Figure 2-39

Temperature distribution (shock free nozzle flow)

43

Figure 2-40

Density distribution (shock free nozzle flow)

To show how far the flow tangency condition is imposed, Figure 2-41 shows the velocity vector plot of the flow inside the nozzle. Figure 2-42 is a zoom on the lower nozzle wall, which assures the robustness of the used method in enforcing the solid wall boundary condition.

0.15

0.1

r

0.05

0

-0.05

-0.1

-0.15 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x

Figure 2-41

Velocity vector plot (shock free nozzle flow)

44

1

Figure 2-42

Zoom on the lower boundary (shock free nozzle flow)

2.7.2.2 Nozzle Flow (With Shock)

The previous nozzle shape is solved but with a predetermined back pressure to generate a normal shock inside the divergent part. So now we have a different exit boundary condition, which is a subsonic exit. According to the theory of characteristics we should define only one dependent variable (pressure) at that boundary. The used backpressure is = 0.5. Figure 2-43 shows the Mach number contours of the resulting solution. In this figure the two dimensional shock is evident and as expected the shock takes the shape of a cap. Figure 2-44 shows the Mach number distribution for the centerline and for the wall compared to the exact 1-D solution as illustrated in [29]. As seen the shock is well captured despite using coarse grid [60X20].

45

0.15

0.1

0.05

0

-0.05

-0.1

-0.15 0

0.1

Figure 2-43

Figure 2-44

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Mach number contours (nozzle flow with shock)

Mach number distribution (Nozzle Flow With Shock)

46

Figure 2-45 shows the pressure contours in which again the shock is evident Figure 2-46 shows the pressure distribution, which assures the robustness of the used method despite using coarse grid. Figure 2-47 shows the convergence history of the flow solver. In this figure we started the solution by large time step to reach the steady state solution faster then the time step was reduced to refine the solution. In this figure, the quadratic convergence feature is evident which assures the robustness of the flow solver. Figure 2-48 shows a velocity vector plot, which again assures the shock existence.

0.15

0.1

0.05

0

-0.05

-0.1

-0.15 0

0.1

Figure 2-45

Figure 2-46

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Pressure contours (Nozzle Flow With Shock)

Pressure distribution (Nozzle Flow With Shock)

47

1

Figure 2-47

Convergence history (Nozzle Flow With Shock)

Velocity Vector Plot 0.4

0.3

0.2

0.1

0

-0.1

-0.2

-0.3

-0.4

0

0.1

Figure 2-48

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Velocity vector plot (Nozzle Flow With Shock)

48

1

2.7.2.3 Jet Flow (Nearly Fully Expanded)

The current test case is that of a supersonic jet flow. An experimental investigation of flow of a moderate Reynolds-number and Mach-number axisymmetric, near isentropically expanded, cold jet has been performed by Troutt and McLaughlin [30], with the following jet conditions Mach number (Mj)

2.1

Reynolds number (Re)

70,000

Jet / Nozzle Radius (rj)

0.005 mm

Total Temperature (To)

294

Static Pressure (P)

5060 N/m2

Center line Jet Temperature (Tj)

156.4 K

Center line Sonic Speed (Cj)

250 m/s

Center line Jet Speed (Uj)

525 m/s

Center line Jet Density (ρj)

0.1114 Kg/ m3

Before we go further we have to clear some definitions related to the mean flow field shown in Figure 2-49. Potential Core Region: It is the region where the axial velocity equals the jet velocity and no radial velocity exists. The Shear Mixing Region: It is the region where the axial velocity is less than half of its centerline value. To well capture the free jet phenomenon the solution domain is used as shown in Figure 2-50

49

Figure 2-51 (top) shows the domain of solution used in the current study which extended 70Rj in the axial direction and 20Rj in the radial direction using clustering in the normal direction near the nozzle end to well capture the rapid variations in the mixing layer region. Bi-linear elements are used with 60*26 elements in the axial and normal directions respectively.

Figure 2-49

Figure 2-50

Supersonic free jet main regions

Axisymmetric jet flow boundaries

50

Figure 2-51

Computational domain and center line Mach number distribution (jet flow problem)

The centerline Mach number distribution for both the experiment and the current work is shown in Figure 2-51 (bottom). The numerical results are in good agreement with the experiment despite using coarse grid. The figure shows the potential core length of the jet to be 10 radii while that of the experiment between 8 and 10 radii. The sonic point in the jet is reached between 22 and 24 radii while that of the experiment is between 18 and 20 radii. Figure 2-52 shows the axial Mach number distributions for different radial positions in which the jet decelerates and the outer stagnant air is accelerated. Figure 2-53 shows the Mach number contours. In this figure we can see the generic features of the supersonic jet flow; e.g. the potential core region. The symmetry of the flow is depicted in the figure, which assures the quality of the used method despite using coarse grid.

51

Figure 2-52

Axial Mach number distributions (jet flow problem)

Figure 2-53

Mach number contours (jet flow problem)

52

2.7.2.4 Jet Flow (Fully Expanded)

The current test case is that of a supersonic full expanded jet flow with Mj = 2.1. The importance of this test case is that we know the theoretical solution. To well capture the free jet phenomenon the solution domain shown in Figure 2-50 is used again. Figure 2-51 (top) shows the domain of solution used in the current study which extended 70Rj in the axial direction and 20Rj in the radial direction using clustering in the normal direction near the nozzle end to well capture the rapid variations in the mixing layer region. Bi-linear elements are used with 40*20 elements in the axial and normal directions respectively. Figure 2-54 shows the contour plot for the Mach number, which supports the expected theoretical solution, i.e. the expected cylindrical contours. Figure 2-55 shows the contour plot for the radial velocity. The symmetry of the solution is evident (the full domain was solved not half) which assures the quality of the used technique.

Figure 2-54

Mach Number Contours (Jet Flow Full Expanded)

53

Figure 2-55

Radial velocity Contours (Jet Flow Full Expanded)

2.7.2.5 Under Expanded Jet Flow (UEJF)

The current test case is that of a supersonic under expanded jet flow with Mj = 2.1 and pressure ratio pj/pamb. = 1.7. The importance of this test case is that tests the capability of the scheme to capture the well-known cell-structure. To well capture the free jet phenomenon the solution domain shown in Figure 2-56 is used with a [40 X 36] bi-linear elements which extended 50Rj in the axial direction and 20Rj in the radial direction using clustering in the normal direction near the nozzle end to well capture the rapid variations in the mixing layer region. Figure 2-57 shows the contour plot for the Mach number, which supports the expected solution, i.e. the expected cell structure (since the jet flow has a higher pressure value than the ambient, an expansion fan is created at the nozzle exit. This expansion fan is reflected when it reaches the centerline which acts like a wall. The reflected fan turns the slip line inwards while being reflected again and when it reaches the centerline again, it is reflected again to meet the slip line again and turns it outward, thus forming the first cell and so on) in con-

54

clusion the expansion fan and the slip line affects each other to make the slip line takes the shape of a cell. The symmetry of the solution is evident (the full domain was solved not only one half) which assures the quality of the used technique. Figure 2-58 shows the axial Mach number distribution at different radial locations. As expected the Mach number increases till the pressure reaches the ambient pressure as shown in Figure 2-59. Figure 2-60 shows the axial distribution of the density. Figure 2-61 shows the axial velocity distribution at different radial location,. Figures starting from Figure 2-62 to Figure 2-64show the radial variations of the Mach number, pressure, and the radial velocity respectively, which again assures the capability of the technique to capture the symmetry despite solving the whole domain.

Figure 2-56

Computational domain (UEJF)

55

Figure 2-57

Figure 2-58

Contour plot for the Mach number (UEJF)

Axial Mach number distribution (UEJF)

56

Figure 2-59

Axial pressure distribution (UEJF)

Figure 2-60

Axial density distribution (UEJF)

57

Figure 2-61

Figure 2-62

Axial distribution of radial-velocity (UEJF)

Radial distribution of axial Mach number (UEJF)

58

Figure 2-63

Figure 2-64

Radial distribution of pressure (UEJF)

Radial distribution of radial-velocity (UEJF)

59

Chapter 3

Directionally Adaptive Technique for FEM

3.1 Introduction The current interest in the area of high-speed flows has increased the need for advanced computational fluid dynamics (CFD) codes, which have become the primary tools for the prediction of aero-thermal loads. Such flows are characterized by regions with steep directional gradients of flow variables, embedded in regions where the flow variables vary more smoothly. One approach for improving the solution accuracy of such problems is to apply grid adaptation techniques.

3.2 Mathematical Analysis The mathematical analysis is divided into two sections as follow 3.2.1

Edge-Based Error Estimate

Consider a 1D element in which the solution variable g is approximated by gh with linear interpolation [18]. A local error Ee defined over an element e can be estimated as the difference between a quadratic interpolation gq and the actual linear one provided that the error is zero at the nodes and maximum at the middle of the element

Ee =

ζ 2

(h e − ζ )

2 d gq

dx

(3-1)

2 e

where

ζ is the local element coordinate and he the element length.

60

A measure of the overall error in the element is then considered to be the root-mean-square value of Ee

1

RMS

Ee

2  h e E e2 2 gq 1 2 d dζ  = = ∫ he 2  he  120 dx 0 

(3-2) e

If we define an optimal mesh as one for which the error is equidistributed over elements, the following should hold for each element [18]:

h

2 e

2 d gq

dx

=c

2

(3-3)

e

where c is a positive constant. The second derivative in equation (3-2) is based on gq, which is the solution being sought and hence not available. So it is here approximated by the second derivative of the numerical solution, i.e., d 2 g h / dx 2 . The above methodology can be extended to 2D based on the fact that each edge of a 2D element can be considered as a 1D element ([18], [31], [32], [33]). So the second derivative in equation (3-2) may be replaced with the Hessian matrix as follows

 ∂2 g h  ∂x 2  H = 2 ∂ g h   ∂y ∂x

∂2 g h   ∂x ∂y  ∂2 g h   ∂y 2 

(3-4)

61

The second derivative in equation (3-4) will vanish since we use a bi-linear element. However, a mass lumping (center of mass) [18], can be applied to recover an estimate of the second derivative. This yields the expression g

h , ij I

=

∫A I g

h , ij

N I dA

(3-5)

∫A I N I d A

where AI represents the elements sharing node I. After integration of equation (35) by parts, the nodal values of the Hessian reduce to

g h ,ij = I

where Γ

I

∫Γ I g h ,i N I n j d Γ − ∫A I g h ,i N I , j dA ∫A I N I dA

(3-6)

represents contours of the elements sharing node I.

The Hessian matrix given by equation (3-4) may be diagonalized as follows:

H = R (α )ΛR T (α )

(3-7)

where Λ is the diagonal matrix of the eigen-values of H and R is the matrix of the eigenvectors. The transformation Λ is a scaling in the direction of the axes and R is a rotation with angle α that the eigenvector corresponding to the smallest eigen-value makes with the x-axis. In order to obtain a symmetric, positive definite matrix, the Hessian is modified by taking the absolute value of its eigenvalues. This results in

H = R (α ) Λ R T ( α ) = S ( α ) S T ( α )

(3-8)

where S (α ) = R (α ) Λ . The transformation S of a unite circle would be an ellipse, rotated through an angle α , whose semi-major axis is the reciprocal of the

62

square root of the smallest eigenvalue, and semi-minor axis is the reciprocal of the square root of the largest eigenvalue [18], as shown in Figure 3-1. Therefore one can obtain a directionally stretched grid by mapping a uniform mesh using the transformation S. However, in the current approach a mesh with edges of equal length is sought in the transformed plane ST, where the length of a curve B is given by

1

d (B ) =

∫ 0

 s ` ( l ) T H ( l ) s ` ( l )  d l

(3-9)

and s (l ) is a parametric representation of the curve B. Since H is a function of the space coordinates, equation (3-9) defines a Riemannian metric. The modified Hessian is computed and stored on a background mesh and thus the value of H at any position of the domain interpolated during the adaptive process on this mesh. The edge-based error estimate can then be numerically evaluated from equation (3-9) for each edge of the element.

Figure 3-1

Transformation of a unite circle to an ellipse by S

63

3.2.2

Moving-Node Scheme

The adaptive strategy modifies the grid under the guidance of the error estimate to improve the quality of the numerical solution. Thus the use of an appropriate adaptive scheme is crucial for achieving the desired directionally adapted mesh. The used strategy relies on a node-moving scheme, also called nodal redistribution. As illustrated in Figure 3-2, the mesh may be viewed as a network of springs [15] whose stiffness constants represent the edge-based error estimate.

Figure 3-2

Spring analogy for a patch of elements.

The positions of the grid vertices may then be interpreted as the solution of an energy minimization problem. This yields for each vertex I

m in P xI

I

= m in ∑ ( x I − x J ) 2 k IJ xI

(3-10)

J

where PI denotes the potential energy of the four springs sharing a node I and kIJ are the associated four stiffness constants. These constants may be specified as

64

k IJ =

d (x I − x J ) xI −xJ

(3-11)

where x I − x J indicates the Euclidian norm

d (x I − x J ) is the length of the edge [x I − x J ] in the Riemannian metric defined by equation (3-9). After simplification, equation (3-10) reduces to the following system describing the equilibrium state of a spring network:

∑ (x Im +1 − x Jm +1)k

m IJ

J

=0

(3-12)

By lagging xJ and kIJ at the previous iteration m, equation (3-12) becomes

∆x I =

−∑ (x mI − x Jm )k m IJ J

∑k J

m

(3-13)

IJ

and the position of the vertex I is updated according to the expression

m +1

xI

= x Im + ω∆x I

(3-14)

where ω is a relaxation parameter (vector). The convergence of this scheme can be enhanced by using Gauss-Seidel algorithm with the latest values of xJ and kIJ when using equation (3-13). The iteration process in equation (3-14) can be applied to all nodes in the domain in order to adapt the mesh to the solution. Boundary nodes can also move in the same way as internal nodes but they are then re-projected onto the boundary to maintain the geometric integrity of the domain. The moving-node scheme is applied to grid points in a sweeping manner. The reason is to allow checking the quality of each newly oriented element during the mesh movement and thus avoiding formation of elements with a nega-

65

tive or nearly zeroing Jacobian. The adaptive method uses the solution of one of the scalar variables to adapt the mesh. Then goes back to the flow solver with the adapted mesh. Each mesh adaptation followed by the flow solver is called one adaptive cycle.

3.2.3

Grid Smoothening

The previous scheme does not guarantee the smoothness of the resulting grid, i.e. the resulting grid may contain elements with angles greater than 170o or less than 10o, which may cause ill-posedness of the resulting global matrix. So we propose to use additional diagonal elements to act as semi-torsional springs, or to add a grid smoothening step, which guarantees the smoothness of the adapted grid. Which is performed after each iteration of the adaptation scheme by simply equating the gradient of the adapted grid lines at each point (see Figure 3-3) as in equation (3-15) for the y direction and the same applies for the x direction.

Figure 3-3

A node on the grid and its surroundings

yi,j = ((xi,j-xi-1,j)yi+1,j+(xi+1,j-xi,j)yi-1,j)/ (xi+1,j-xi-1,j)

66

(3-15)

3.2.4

The Grid Adaptation Procedure

The grid adaptation procedure may be summarized in the following steps: -----------------------------------------------------------------------------------------Read a background mesh and the corresponding solution Compute H on the background mesh Current mesh is initialized by an initial mesh guess (optional) Move the nodes of the current mesh as follows DO m=1,MAXITER DO inod=1,NNODE DO iedge=1,NEDGE Determine H by interpolating on the background mesh (optional) Compute springs constants by numerical integration of (3.11) ENDDO Find new position of inod Move inod to its new position Check quality of elements sharing node I ENDDO Grid smoothness (optional) If (MAXDISP .<=. tolerance) exit the m loop ENDDO ----------------------------------------------------------------------------------------where MAXITER: maximum number of iterations. NNODE: total number of nodes. NEDGE: number of edges surrounding the current adapted node. inod: counter over the NNODE. iedge: counter over the NEDGE. MAXDISP: maximun allowable displacment. tolerance: allowable tolerance.

67

3.3 Numerical Results 3.3.1

Analytical Test Case

It is important to examine the behavior of the adaptation algorithm on an analytical function with strong gradient to represent its quality and robustness. A function g in the form

(

g ( x,y ) = tan −1 1000 ( x 4 y 4 − 0.25 )

)

(3-16)

is used in the solution domain [0, 2] x [0, 1]. The initial coarse mesh [30*20] is shown in Figure 3-4 and the corresponding isocontours of g are shown in Figure 3-5. After 60 iterations ( ω =0.3) of the mesh movement scheme the adapted mesh shown in Figure 3-6 is obtained. As illustrated in Figure 3-7 this mesh permits a better representation of the function g (the contour lines became smoother) using the same number of nodes. The magnification of the mesh in the discontinuity region presented in Figure 3-8 shows that the quadrilateral elements are strongly re-oriented in the direction of the discontinuity with a very high aspect ratio. Also we can see that the re-oriented element has a very small length in the direction normal to the shock to represent well the strong gradient normal to the shock, and a very high length in the shock direction to reduce the used number of elements, which raises the quality of the used technique.

68

Figure 3-4

Figure 3-5

Initial mesh (analytical test case).

Isocontours for the initial mesh (analytical test case).

Figure 3-6

Adapted mesh (analytical test case).

69

Figure 3-7

Adapted isocontours (analytical test case).

Figure 3-8

Magnification of grid (analytical test case).

70

3.3.2

Shock Reflection Problem (SRP)

The second test case is the shock reflection problem discussed previously in Chapter 2 (section 2.7.1.1). Figure 3-9 shows the computational domain again. This example tests certain features of the algorithm including the resolution of a system of two oblique shocks and their proper angles. This test case will be solved with and without the proposed smoothening step. Figure 3-10 shows the initial grid (a) and the corresponding pressure contours (b). Also in this figure the 1st adapted grid is shown (c) as well as the corresponding pressure contours (d). Figure 3-11 shows the results for the 2nd and 3rd adaptation cycles. In these figures one can see that the resulting grids are not adapted well to the required solution. Therefore in the following figures we will use the proposed smoothening step.

Figure 3-9

The computation domain (SRP).

71

1

(a)

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1

(b)

0.5

0 1

(c)

0.5

0 1

(d)

0.5

0

Figure 3-10

The initial and 1st adapted pressure (SRP)

1

(a)

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1

(b)

0.5

0 1

(c) 0.5 0 1

(d)

0.5

0

Figure 3-11

The 2nd and 3rd adapted pressure (SRP)

72

The initial coarse mesh composed of [61 x 21] nodes and the corresponding isocontours are shown in Figure 3-12 . Using ∆t=0.1, the results are adapted through five cycles using the pressure or the density as the adaptation parameter as specified in the title of each figure. The artificial viscosity is reduced beginning at the third cycle by reducing the time step to 0.05, since its amount was very high for the size of the grid near the shock. Figure 3-13 through Figure 3-17 show the pressure contours and the grids after each adaptation. The improvement in the shock resolution after adaptation is quite evident. Figure 3-18 and Figure 3-19 show how the elements are re-oriented to be aligned with the shock, creating very high aspect ratio elements, which assures the robustness of the technique. The pressure distribution at (y=0.5) is shown in Figure 3-20. The adapted solution captures the shock more sharply. The convergence history is shown in Figure 3-21. Each jump in the Figure corresponds to an adaptation cycle. The quadratic convergence of the Newton linearization is quite evident. The Mach number contours for the initial solution and for the adapted solution are shown in Figure 3-22 through Figure 3-27, which again assures the importance of the grid adaptation to sharply capture shocks despite using coarse mesh.

73

1 0.8

y

0.6 0.4 0.2 0

0

0.5

1

1.5

2 x

2.5

3

3.5

4

0

0.5

1

1.5

2 x

2.5

3

3.5

4

1 0.8

y

0.6 0.4 0.2 0

Figure 3-12

The initial grid and the and the pressure contours (SRP)

1 0.8

y

0.6 0.4 0.2 0

0

0.5

1

1.5

2 x

2.5

3

3.5

4

0

0.5

1

1.5

2 x

2.5

3

3.5

4

1 0.8

y

0.6 0.4 0.2 0

Figure 3-13

The 1st adapted grid and the pressure contours (SRP)

74

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-14

The 2nd adapted grid and the pressure contours (SRP)

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-15

The 3rd adapted grid and the pressure contours (SRP)

75

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-16

The 4th adapted grid and the pressure contours (SRP)

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-17

The 5th adapted grid and the pressure contours (SRP)

76

0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.4

Figure 3-18

0.5

0.6

0.7

0.8

0.9

1

Zoom on the grid in the incident shock region SRP

0.56

0.54

0.52

0.5

0.48

0.46

0.44 2.6

Figure 3-19

2.8

3

3.2

3.4

3.6

3.8

Zoom on the grid in the reflected shock region (SRP)

77

Figure 3-20

Pressure distribution at y=0.5 (SRP)

Figure 3-21

Convergence history (SRP)

78

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-22

The initial grid and the Mach number contours (SRP)

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-23

The 1st adapted grid and Mach number SRP

79

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-24

The 2nd adapted grid and Mach number (SRP)

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-25

The 3rd adapted grid and Mach number (SRP)

80

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-26

The 4th adapted grid and Mach number (SRP)

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-27

The 5th adapted grid and Mach number (SRP)

The problem of reflected shock is re-solved using the density as the adaptation parameter in order to know how far the effect of the adaptation parameter is. See Figure 3-28 through Figure 3-33 for the evolution of the grid and the solution during adaptation, which again supports the importance of the adaptation.

81

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-28

The initial grid and pressure (density adapted) SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-29

The 1st adapted grid and pressure (density adapted) SRP

82

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

the 2nd adapted grid and pressure (density adapted) SRP

Figure 3-30

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-31

The 3rd adapted grid and pressure (density adapted) SRP

83

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

The 4th adapted grid and pressure (density adapted) SRP

Figure 3-32

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-33

The 5th adapted grid and pressure (density adapted) SRP

Comparison of pressure contours is shown in Figure 3-34 through Figure 3-38 using the pressure as the adaptation parameter as well as the density as the adaptation parameter. From these figures we can notice that the density-adapted results are less oscillatory than those of the pressure adapted ones during adaptation cycles, but the final result is biased toward the pressure adapted ones

84

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-34

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 1st adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-35

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 2nd adp. SRP

85

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-36

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 3rd adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-37

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 4th adp. SRP

86

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-38

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 5th adp. SRP

Figure 3-39 through Figure 3-43 also show the resulting grids using the pressure and the density as the adaptation parameter. From these figures it can be seen that the density-adapted grids are more refined and oriented to the shock direction. The difference between the adapted results using density and pressure as an adaptation parameter is evident in Figure 3-44. To complete the comparison, Figure 3-45 through Figure 3-49 show the Mach number contours for the pressure-adapted case (top) and the density adapted case (bottom). From these figures one can conclude that the density as an adaptation parameter is recommended since it produces less oscillations during adaptation cycles.

87

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-39

Comparison of pressure adapted (upper) and density adapted (lower) grid for the 1st adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-40

Comparison of pressure adapted (upper) and density adapted (lower) grid for the 2nd adp. SRP

88

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-41

Comparison of pressure adapted (upper) and density adapted (lower) grid for the 3rd adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-42

Comparison of pressure adapted (upper) and density adapted (lower) grid for the 4th adp. SRP

89

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-43

Comparison of pressure adapted (upper) and density adapted (lower) grid for the 5th adp. SRP

3.5 Exact press adp. dens.adp no adp

3

2.5

2

1.5

1

0.5

0

Figure 3-44

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Comparison of pressure distribution at y=0.5 for SRP

90

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-45

Comparison of pressure adapted (upper) and density adapted (lower) Mach number contours for the 1st adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-46

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 2nd adp. SRP

91

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-47

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 3rd adp. SRP

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-48

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 4th adp. SRP

92

1 0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

1 0.8 0.6 0.4 0.2 0

Figure 3-49

Comparison of pressure adapted (upper) and density adapted (lower) pressure contours for the 5th adp. SRP

Figure 3-50 shows how important is the role of the adaptation. This figure shows the pressure contours for reference [34] (middle) and the pressure contours for the current work (top). The importance of this figure lies in that reference [34] uses a [121 X 81] grid i.e. (9801) nodes, while in the current work the used grid is [61 X 21]; i.e. (1281) nodes only, which saves computational time considerably. Also shown in this figure (bottom) are the pressure contours for [11].

Figure 3-50

Comparison of pressure contours (SRP)

93

3.3.3

Supersonic Channel Problem (SCP)

A more challenging problem for this algorithm exists in the analysis of supersonic flow past the parabolic arc bump discussed previously in Chapter 2 (section 2.7.1.2). The mesh used for this problem consists of 1024 elements and 1088 nodes evenly distributed in a domain which has a bump height of 4%. This problem creates an interesting shock interaction behind the bump. The shock from the leading edge of the bump is reflected down off the ceiling of the tunnel and crosses, then combines, with the shock formed at the trailing edge of the bump. The initial coarse mesh composed of [65 x 17] nodes and the corresponding isocontours are shown in Figure 3-51. Using ∆t=0.1, the results are adapted through five cycles. The artificial viscosity is reduced beginning at the second cycle by reducing the time step to 0.05, since its amount was very high for the size of the grid near the shock. Figures starting from Figure 3-52 to Figure 3-56 show the pressure contours and the grid after each adaptation. The improvement in the shock resolution after adaptation is quite evident. The density contours, for the initial solution and for the adapted solution, are shown in Figure 3-57, which again assures the importance of the grid adaptation to sharply capture shocks despite using coarse mesh. Figure 3-58 shows a magnification of the grid at the zone of shock reflection at the upper wall, which demonstrates the robustness of the algorithm in reorienting the elements in the shock direction with very high efficiency. Finally Figure 3-59 shows the convergence history of the flow solver, each jump corresponds to an adaptation cycle, the quadratic convergence is evident in the figure also, which assures the flow solver and adaptation algorithm robustness.

94

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-51

The initial grid and the corresponding pressure contours

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-52

The 1st adapted grid and the corresponding pressure cont.

95

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-53

The 2nd adapted grid and the corresponding pressure cont.

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-54

The 3rd adapted grid and the corresponding pressure cont.

96

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-55

The 4th adapted grid and the corresponding pressure cont.

1 0.8 0.6 0.4 0.2 0 -1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

1 0.8 0.6 0.4 0.2 0 -1.5

Figure 3-56

The 5th adapted grid and the corresponding pressure cont.

97

1

initial solution dt = 0.1

0.5

0 -1.5 1

2nd adp. dt = 0.05

4th adp. dt = 0.05

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

-1

-0.5

0

0.5

1

1.5

0.5

0 -1.5 1

3rd adp. dt = 0.05

-1

0.5

0 -1.5 1 0.5

0 -1.5

Figure 3-57

Evolution of density contours during adaptation

98

0.95

0.9

0.85

0.8

0.75

0.7

0.65

0.2

Figure 3-58

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Magnification of grid in reflected shock region

Figure 3-59

Convergence history

99

Chapter 4

Summary and Conclusions

4.1 Thesis Summary In chapter 1, the aim of the thesis has been introduced in section

1.1‘building a least squares finite element Euler solver for the transonic regime combined with an error estimate and a directionally adaptive grid algorithm’, and then the finite element method has been introduced in section 1.2. In section 1.3 the finite element advantages has been introduced followed by the different approaches used in the finite element formulation in section 1.4. An extensive literature review for the compressible Euler’s equation has been presented in section 1.5. Then an extensive literature review for the adaptation techniques has been presented in section 1.6. Finally, an overview of the current work has been presented. In chapter 2, the least squares finite element method for the compressible

Euler’s equation has been presented. First, the Least squares finite element advantages have been introduced in section 2.1. Then, the least squares formulation has been presented in details in section 2.2. An interpretation of the inherent artificial viscosity has been presented in section 2.3. The finite element approximation has been presented and the problem details as shape function and local coordinates have been reviewed in section 2.4. The boundary condition, one of the most important aspects in inviscid flows has been presented in section 2.5 using a very strong technique in applying the no-slip condition. In section 2.6 the solution method has been presented.

100

To show how robust and stable is the used technique, in section 2.7, a variety of a well-known problems concerning the planar and Axisymmetric flows have been solved, analyzed, and compared to exact or experimental results showing the superiority of the used technique over the existing ones specifically in that it uses a coarse mesh to capture discontinuities. In chapter 3, a directionally adaptive methodology for the finite element

method has been presented in details. An introduction has been introduced in section 3.1. In section 3.2, the mathematical analysis has been divided into two sections. First, in section 3.2.1, the error estimation has been presented using the second derivative of a given flow variable “Hessian matrix” compared to first derivative used by most other adaptive methods. Second, in section 3.2.2, a moving node scheme based on spring analogy was presented in details followed by an adaptation procedure summarized in few steps. To validate the used scheme, an analytical test case has been used demonstrating the capability of the used scheme in presenting the contours of a strong-gradient function despite using coarse mesh. Two numerical test cases have been used to show the capabilities of the used scheme. The first one was the shock reflection problem, which tests the capability of the scheme in capturing a system of two oblique shocks at the correct angle. Then, a complicated shock structure problem has been solved to show the capability of the scheme in re-orienting the grid lines with shocks despite such a complex shock structure.

4.2 Conclusions The least squares finite element method has been used to solve the compressible Euler's Equations for both the Cartesian and Axisymmetric flows in the non-conservative form. The quality of the numerical results indicates the remarkable performance of the used technique. This is quite evident from the final results especially for the nozzle flow as well as the jet flow problem despite using coarse mesh. Also the robustness of the technique allows one to use large time

101

step to reach the steady state solution in few iterations, as done in the shock reflection problem (∆t=0.15), which saves computational time considerably. But a disadvantage has been encountered; namely, the large inherent artificial viscosity in the method, which has prevented the sharp resolution of discontinuities, this disadvantage, has been remedied. To remedy this disadvantage in the leastsquares finite element method, an adaptive algorithm with directional features, using an edge based error estimate on quadrilateral meshes, has been used. The error of the numerical solution is measured by its second derivatives and the resulting Hessian tensor is used to define a Riemannian metric. A mesh movement strategy with no orthogonality constraints is used to equidistribute the lengths of the edges of the elements in the defined metric. The adaptive procedure has been proven to be effective for an analytical test case. Dealing with the boundary conditions, the used technique in treating the solid wall boundary is very robust since it enforces the boundary condition as a Dirichlet type at each node using a transformation matrix for each node this results in a faster convergence than other methods and it can be used for both steady and unsteady flows unlike other techniques. The flow solver, combined with the proposed grid adaptation method is then validated on shock-reflection problem, capturing the system of two shocks with high resolution and the correct angle despite using coarse mesh. The methodology is also tested on the supersonic channel flow problem. The quality of the numerical results indicates the remarkable performance of the adaptive method, and demonstrates its superiority to many existing techniques. This is quite evident from the final adapted grids especially for the shock-reflection and supersonic channel problems where elements are strongly aligned with the shocks. Overall, the combination of the least-squares method with the directionally adaptive method has produced promising results, despite relatively coarse grids using an artificial viscosity mechanism, which has no free parameters.

102

4.3 Recommendations for future work The following topics should be addressed in subsequent studies Concerning the least-squares method, using the conservative variables to overcome the problem of the oscillations after the leading- and trailing edge shocks. Concerning the adaptive scheme, modifying the adaptive methodology, by adding torsional springs to account for element's angles deformations to prevent elements having zero or negative Jacobian.

103

References [1]. J. C. .A. Villegas, “Anisotropic Adaptive Refinement Algorithms for Finite Element Methods”, Ph.D. thesis, Graduate School of Arts and Science, Department of Mathematics, New York University, September 2000. [2]. F. Taghaddosi, W.G. Habashi, G. Guevremont, and D. Ait-Ali-Yahia “An Adaptive Least-Squares Method for the Compressible Euler Equations”, Int. J. Num. Meth. Fluids, 31 (1999) 1121-1139. [3]. D. Ait-Ali-Yahia and W.G. Habashi “Finite Element Adaptive Method for Hypersonic Thermochemical Nonequilibrium Flows”, AIAA J., 35 (1997) 1294-1302. [4]. W.G. Habashi 'The FEM in Fluid Mechanics' Coarse Notes, Concordia University, Canada, ENGR 626. [5]. A.L. Brooks and T.I.R. Hughes, “Streamline Upwind Petrov-Galerkin formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations”, Comp. Meth. Appl. Mech. Eng. 32 (1982) 199-259. [6]. G.S. Baruzzi, W. G. Habashi and M.M. Hafez, “Finite Element Solutions of the Euler Equations for Transonic External Flows”, AIAA J. 29 (1991) 1886-1893. [7]. J.F. Polk and P.P. Lynn, “A Least Squares Finite Element Approach to Unsteady Gas Dynamics.”,' J. Comp. Phy. 12 (1978) 3-10.

104

[8]. C.A.J. Fletcher, “A Primitive Variable Finite Element Formulation for Inviscid Compressible Flows”,' J. Comp. Phy. 33 (1979) 301-312. [9]. F. Moussaoui “A Unified Approach for Inviscid Compressible and Nearly Incompressible Flow by Least-Squares Finite Element Method”, Applied Numerical Mathematics 44 (2003) 183-199. [10]. P.P. Lynn and S.K. Arya, “Use of the Least Squares Criterion in the Finite Element Formulation”, Int. J. Num. Meth. Eng. 6 (1973) 75-88. [11]. J.P. Pontaza, Xu Diao, J.N. Reddy, and K.S. Surana ”Least-Squares Finite Element Models of Two-Dimensional Compressible Flows”, Finite Element in Analysis and Design 40 (2004) 629-644. [12]. P. Bolton and R.W. Thatcher,”A Least-Squares Finite Element Method for the Navier–Stokes Equations”, J. of Comput. Phys., 213 (2006) 174– 183. [13]. J. Peraire, M. Vahdati, K. Morgan and O. C. Zienkiewicz, “Adaptive Remeshing for Compressible Flow Computations”, J. Comp. Phy. 72 (1987), 449-466. [14]. M. Fortan, M. –G. Vallet, D. Poirier and W. G. Habashi, “Error Estimation and Directionally Adaptive Meshing”, AIAA Conf. Paper 94-2221, 1994. [15]. P. A. Gnoffo, “A Finite-Volume Adaptive Algorithm Applied to Planetary Flow Fields”, AIAA J., 21 (1983) 1249-1254. [16]. K. Nakahashi and G. S. Deiwert, “Self-Adaptive Grid Method with Application to Airfoil Flow”, AIAA J. 25 (1987) 513-520.

105

[17]. G. F. Carey and J. T. Oden, “Finite Elements Computational Aspects”, Vol. III, Prentice-Hall, Inc. new Jersey, 1984. [18]. D. Ait-Ali-Yahia, W.G. Habashi, A Tam, M.-G. Vallet, and M Fortan ”A Directionally Adaptive Methodology Using An Edge-Based Error Estimate on Quadrilateral Grids”, Int. J. Num. Meth. Fluids, 23 (1996) 673-690. [19]. J.N. Reddy, “An Introduction to the Finite Element Method”, McGrawHill, New York, 1993. [20]. G.j. Le Beau, T.e. Tezduyar “Finite Element Computation Of Compressible Flows With The SUPG Formulation”, Advances in Finite Element Analysis in Fluid Dynamics. ASME 123 (1991) 21-27. [21]. N.E. Elkadri E, A. Soulaaimani, C. Deschenes, “A Finite Element Formulation of Compressible Flows Using Various Sets of Independent Variables”, Comput. Methods Appl. Mech. Engrg. 181 (2000) 161-189. [22]. S. Eidelman, P. Colella, and R. P. Shreeve,”Application of the Godunov Method and its Second-Order Extension to Cascade Flow Modeling”, AIAA J. 22 (1984) 1609-1615. [23]. N. J. Pfeiffer, and G. W Zumwalt, “Computation Model for Low Speed Flows Past Airfoils with Spoilers”, AIAA J. 20 (1982), 376-382. [24]. M. M. Abdelrahman “Numerical and Experimental Investigations of the Flow Field Past Airfoil-Spoiler Configuration” Proceedings of the Picast conference”, Tainan, 6-9 Dec. (1993), 225-232. [25]. W. H. Jr Wentz. and C. Ostowari, “Effect of Design Variables on Spoiler Control Effectiveness, Hinge Moments and Wake Turbulence”, AIAA Conf. Paper 81-0072, Jan. 1981.

106

[26]. M. D.Mack, H. C. Seetharam, Kuhn, W.G., and Bright, J. T. “Aerodynamics Spoiler Control Devices”, AIAA Conf. Paper 79-1873, 1979. [27]. C. S. Lee, and S. Bodapati, “Experimental Investigation of the Flow Field of an Airfoil with Spoiler”, AIAA J., 25 (1987), 1411-1416. [28]. S. W. Choi, K. S. Chang, H. Ok “Parametric Study of A Deploying Spoiler with Two-Equations Turbulence Models”, AIAA 2001-0864. [29]. J. D. Anderson, Jr., “Modern Compressible Flow”, McGraw-Hill, New York, 1982. [30]. T.R. Troutt and D.K. McLaughlin, “Experiments on the Flow and Acoustic Properties of a Moderate-Reynolds-Number Supersonic Jet”, J. Fluid Mech. 116 (1982) 123-156. [31]. F. S.-Gulick, C. Y. Lepage, and W.G. Habashi “Anisotropic 3-D Mesh Adaptation for Turbulent Flows”, AIAA Conf. Paper 2004-2533,2004. [32]. G. Xia, D. Li, and C. L. Merkle “Anisotropic Grid Adaptation on Unstructured Meshes”, AIAA Conf. Paper 2001-0443,2001 [33]. C. Y. Lepage, F. S.-Gulick, and W. G. Habashi, “Anisotropic 3-D Mesh Adaptation on Unstructured Hybrid Meshes”, AIAA Conf. Paper 20020859,2002. [34]. Z.-C. Zhang and S.-T. Yu, “Shock Capturing without Riemannian Solver- A Modified Space-Time CE/SE Method for Conservation Laws”AIAA Paper 99-0904,1999.

107