A study of two dimensional separated flow by a combination of the finite element method and Navier-Stokes equations

KJELL HERFJORD

A STUDY OF TWO-DimENsioNAL

SEPARATED FLOW BY A COMBINATION

OF THE FINITE ELEMENT METHOD

AND NAVIER-STOKES EQUATIONS

NT

.

TECHNISCHE UNIVERSITHIT

Scheepshydrameehanica

Archief

Mekel7eg 2, 2628 CD Delft

Te1015-786873/Fax:781836

Q

Q0011

V-.1

I

ILI

.0, I

I

1111,

i

N

UNIVERS TETET TRONDHEIM

NoRGES TEKNISKE HOGSKOEJE NrEA-rapport,199&11i

DOKTOR INGENIORAVHANDLING 1996:2 INSTITUIT FOR MARIN HYDRODYNAMIKK TRONDHEIM

NTHTrykk 1996 TECHNISCHE UNIVERSITEIT Laboratorlum voor Scheepshydromechanica.

Archie

Mekelweg Z 2628 CD Delft kt 015-786673- Fax:0154 781838

A Study of Two-dimensional

Separated Flow by a Combination of

the Finite Element Method and

Navier-Stokes Equations

A Thesis Submitted for the Degree of Dr lug

Kjell Herfjord

Department of Marine Hydrodynamics

The Norwegian Institute of Technology

1995

Summary

A Navier-Stokes solver is implemented based on the finite element method. The

incom-pressible, laminar form of the equations are solved in two dimensions. The element method

is based on a standard Bubnov-Galerkin form in space and forward Euler in time. The inherent underdiffusivity of the Galerkin and the forward Euler methods is compensated

for by using a so-called balancing tensor viscosity. This is the two-dimensional form of the negative diffusivity deduced from the 1D advection-diffusion equation.

The element method is implemented in two dimensions by the use of triangular elements

and linear interpolation functions. The resulting algebraic equations are solved in three

steps according to the so-called velocity correction method. In two of the steps. the basic version is solved by simplifying the system of equations by applying "mass lumping'.

The validity of the basic solution method is verified by comparing with results from the

use of second order time integration (Runge-Kutta 2) and by solving the equation systems with consistent mass. The basic version is extended by the implementation of solution in

an accelerated coordinate system. The solution is referred to the earth-fixed coordinate system through a transformation term. The implementation is restricted to translatory (not rotating) coordinate system, and may be used to simulate e.g. a situation of lock-in.

In addition, different possibilities of implementing boundary conditions are tested.

The solution of the velocity/pressure field may be checked by a method of a posteriori error estimation which is implemented. Experience from using the method is reported.

Further, the conservation of momentum and kinetic energy is evaluated via two different ways of computing each quantity. The difference (residual or unbalance) between the two approaches expresses a measure for inaccuracy of the solution. The expressed inaccuracy may be used to correct the model of the problem, by refining the mesh and/or by adjusting the domain.

The method is compared with analytical solutions of the Navier-Stokes equations (Stokes' first problem and Couette flow) as well as a semi-analytical solution (Blasius' solution for

the boundary layer over a flat plate). The method is validated through comparison with

measurements and other computations. The validation tests have been relatively extensive, and contain tests on two-dimensional sections of both circular and square cross-sections. Both steady incident and planar oscillating flow have been used. Tests are performed with single cylinders as well as two in different relative positions (side-by-side and in tandem)

and in different distances to each other. The moving body option is used to simulate the

situation of lock-in and to assess the range of lock-in for different motion amplitudes. A demonstration of how the method can be applied in offshore engineering is given for a

case where a section of the columns of a semisubmersible is modelled. In this particular

icY

example, the effect on damping of the low frequency motions from superimposed wave

frequency relative velocities is shown. The influence on low frequency damping from wave frequency motions are shown to be strong. A further development of the methodology to handle full scale flow conditions is recommended.

Acknowledgements

This study has been carried out at the Department of Marine Hydrodynamics at the

Norwegian Institute of Technology (NTH), under the supervision of Professor Odd M.

Faltinsen. His continuous engagement and guidance is highly acknowledged. At Norsk Hydro, Dr. Finn Gunnar Nielsen has been supervising the work. The discussions with him

has been fruitful and ensured that the work has kept momentum.

Further. I am very grateful to Dr. Torbjorn Utnes for his valuable advices on the use of the finite element method to solve the Navier-Stokes equations at the outset of this

work. I also want to thank the fellow students and the staff at the Departments of Marine

Hydrodynamics and Marine Structures who contributed to make the process enjoyable

both socially and scientifically.

The study has received a grant from the University of Trondheim (UNIT). Financial sup-port has also been received from Norsk Hydro. The computing time on CRAY Y-MP has been supported by the Norwegian Supercomputing Committee (TRU).

Contents

Summary

Acknowledgements

4

Contents

Nomenclature

9 1

Introduction

13

1.1 Background and motivation 13

1.2 The organization of the thesis 17

2 The governing equations and solution methods

18

2.1 Introduction and discussion of methods 18

9.2 The Navier-Stokes equations, basic incompressible version 20 2.2.1 The introduction of accelerated coordinate system 22

2.3 The finite element method 23

2.3.1 Basis and background 93

2.3.2 Spatial discretization 94

2.3.3 Choice of basis functions 25

2.3.4 Time integration 26

2.3.5 Matrix equations 27

2.3.6 Solution technique 99

3 Convergence, consistency, stability

31

3.1 Introduction and statements on convergence and consistency 31

3.9 Stability 33

3.2.1 Introduction 33

:3.9.9 General consideration in terms of matrix expressions 33

3.2.3 Stability criteria based on the von Neumann method 36

3.2.4 The von Neumann method applied on the advection-diffusion equation 37

3.2.5 Modified stability criteria with balancing tensor viscosity 39

4 Boundary conditions for the Navier-Stokes equations

41

4.1 Introduction 41

4.2 Traction-free boundary condition 49

4.3 Variation of pressure boundary condition at outflow boundary 45 4.4 Variation of velocity boundary conditions at upper and lower boundaries 46

4.5 Conclusions 47 4 z . . . . . . . 2

'5

5

Consistent mass and higher order time integration

48

5.1 Introduction :1:i, 48' 5.2 Outline of methods 48 5.3 Numerical results . , 5.4 Discussion . .. .... , 49' .50

6,

Conservation of mass, momentum and energy'

52

6.1 Introduction 52.

6.2 Continuity . 52

6.3 Momentum 53

6.4 Energy 54

'6.4.1 Development of equilibrium for time rate of change of kinetic energy 54

8.4.2 Alternative approach 55

6.4.3 Analysis of difference between the two approaches 56

16.5 Numerical results . . . 57

6.5.1 Introduction . 57

6.5.2 Numerical methods. . 58

6.5.3 Result presentation . 60'

8.5.4 Discussion of results. 64

7 A posteriori error estimation

69

7.1 Introduction 69'

7.2 Basis for the estimator and the energy norm 69

7.3 Numerical results 79

.7.3.1

Discussion . 76

7.4 'Conclusions _it 77

8h Comparisons with, analytical solutions

80'

8.1 Introduction SO

8.2 Stokes' first problem

...

. 80

8.3 Flow formation in Couette motion . . . _{'1 V'} '82

8.4 The flat plate, Blasius' solution '84

8.5 Discussion 93

8.6 Conclusion..

Validation of the numerical method

96

9.1

Introduction ...

. 96

9.2 Lid driven cavity .... . 97

9.2.1 Introduction . . 97

9.2.2 Results . .... 97

9.3 Impulsively started flow . . r., 11 V,1.1. 100

9.3.1 Introduction

,

100

9.3.2 Results.

...

100

.9.4 The circular cylinder in uniform incident flow 104

7.3.1

95

6

9.4.1 Introduction 104

9.4.2 Numerical results 107

9.5 Circular cylinder in oscillatory ambient flow 116

9.5.1 Introduction 116

9.5.2 Method of analysis 116

9.5.3 Important parameters 117

9.5.4 Summary of observations from Sarpkaya (1986) 118

9.5.5 Numerical results 119

9.6 Quadratic section in steady incident flow 126

9.6.1 Introduction 126

9.6.2 Numerical results 127

9.6.3 Comparison of results 131

9.7 Quadratic section in planar oscillatory flow 134

9.7.1 Introduction 134

9.7.2 Numerical results 134

9.8 Steady incident flow past two cylinders in tandem 140

9.8.1 Introduction 140

9.8.2 Numerical results 141

9.9 Oscillating flow around two cylinders in tandem 142 9.10 Steady incident flow past two side-by-side cylinders 144

9.10.1 Introduction 144

9.10.2 Review of literature 144

9.10.3 Computations with present numerical method 148 9.10.4 Comparison between measured and computed results 161

9.11 The simulation of lock-in 165

9.11.1 Introduction 165

9.11.2 Numerical results and comparison with other investigations 169

10 Viscous damping of slow drift motions

190

10.1 Introduction 190

10.2 Computational procedure and results 190

10.3 Summary 198

11 Discussion and suggestions for further work

199

11.1 Suggestions for further work 200

12 Conclusions 202

References

204

A Summary of applied mesh models

213 . . . .

... .

. .

...

. . . .

. ...

. . . .

.. , ...

The

B.1 B.9

B.3 B.4

finite element method

Introduction A model problem

B.2.1 Discretization of the problem

B.9.2 Matrix equations Numerical example

The finite element method for the Navier-Stokes equations

B.4.1 Governing equations and preliminaries

B.4.9 A FEM of the Navier-Stokes equations

B.4.3 Time integration by the velocity correction method 8.4.4 Weighted residual formulation

B.4.5 Discretization of the spatial domain 8.4.6 Matrix formulation

8.4.7 Solution of the equations

B.4.8 Convergence, consistency and stability

B.4.9 SUPG (streamline upwind/Petrov-Galerkin) method B.4.10 Status of the method

215 215 215 216 217 218 221 221 221 222 299 225 226 229 999 230 231

C Element formulation of 1D advection diffusion equation

232

An introduction to the von Neumann method

236

The underdiffusivity of forward Euler.

238

F The flat plate

240

F.1 Introduction 240 F.2 Test conditions 241 F.3 Results 243 B D E . .

8

Nomenclature

Roman symbols:

A amplitude of flow particle motion, amplitude of motion

.A advection matrix in Na.vier-Stokes solution by FEM matrix in system of linear DE considered in Chapter 3

wall friction coefficient, Cf = A

cp pressure coefficient, cp = pfr

cell Courant number, C = in one dimension added mass coefficient of accelerating body

CD drag (in-line) force coefficient (also written Cd)

lift (transverse) force coefficient (also written Cl)'

Cm mass coefficient of cylinder in accelerated flow field (also written Cm)

AC matrix in system of linear DE considered in Chapter 3 LI chalacteristic length (e.g. diameter)

gradient matrix in Navier-Stokes solution by FEM,,

Dr = ](Dy,Dy) - divergense matrix

ie error

e 11,2 error in the energy. norm,

Eki

kinetic energy, Chapter 6

diagonalized matrix in system of linear DE considered in Chapter 3

(C diagonalized)

oscillation frequency of forced oscillation of cylinder in study of lock-in,,

'Section 9.11

natural frequency

ft

,Strouhal vortex shedding frequency for affixed cylinder

(to be distinguished from f a general vortex shedding

frequency,, which may be influenced by motion of the cylinder) vortex shedding frequency

acceleration of gravity

dimension transverse to the cylinder axis of a test section (wind or water tunnel)

imaginary unit, z = N/-1 unit vector in s-direction

unit vector in y-direction

KC Keulegan-Carpenter number, KC = LIT" = 2,4

If

subscript denoting "low frequency"

Lc length of cylinder exposed to flow in test section (wind or water tunnel)

mass matrix in Navier-Stokes solution by FEM

momentum M = (Mi. MOT in two dimensions. Chapter 6 outward pointing unit normal vector (on boundary F),

n

(nx,n)7. in two dimensions

vector of interpolation functions used in FEM dynamic pressure

Pe cell Reynolds number (or Peclet number),

Pe =

= in one dimension residual

Re Reynolds number, Re =

Ret Reynolds number based on length of plate

Re z Reynolds number based on the distance from the leading edge of a plate

s 112 energy norm of velocity/pressure field centre distance between two cylinders

tic nondimensional frequency of forced oscillation, S,

Si

Strouhal number, Si = 4

nondimensional vortex shedding frequency, S= 4

SBTV diffusion matrix in Navier-Stokes solution by FEM

Laplacian matrix in Navier-Stokes solution by FEM time variable

ti traction (or stress) in direction i on boundary

ii velocity vector, it = (u,v)T in two dimensions

10

general symbol for characteristic velocity

Uo wall velocity

Um amplitude of harmonic oscillating velocity

reduced velocity, Ur =

Uo. velocity at infinity (also written Uinf)

V, normal velocity of boundary

weighting function in finite element approximation

wf

subscript denoting "wave frequency" WD damping work from one oscillation period

space coordinate, x = (x , 07. in two dimensions

Greek symbols:

Stokes' parameter, 8 = 1-% = boundary of computational domain boundary layer thickness

Kronecker delta

Ax increment in space coordinate

Lt

increment in time variable rate-of-strain tensor

average (global) relative error

dimensionless coordinate for boundary layers, 77

0 angle in the range [0, 27r]

12 dynamic viscosity

ii kinematic viscosity amplification factor local relative error

density

total stress tensor

standard deviation

To(x) wall friction, To(x) = p (g)_y=0

if

=

tangent unit vector (on boundary F)

stress tensor

general unknown scalar function in DE phase angle stream function angular frequency vorticity computational domain

Abbreviations

BTV balancing tensor viscosity

CFD computational fluid dynamics

DE differential equation

DNS direct numerical simulation (of turbulence) FDM finite difference method

FEM finite element method FVM finite volume method

LES large eddy simulation

ODE ordinary differential equation PDE partial differential equation

SCPG streamline upwind:Petrov-Galerkin

VIC vortex-in-cell

1

Introduction

1.1

Background and motivation

Oil production platforms have traditionally been fixed to the sea bottom. Steel jackets and

concrete platforms as the Condeep are examples of platforms standing on the sea floor. At larger water depths, it is necessary to use floating production platforms. The floating

production units may grossly be divided into column stabilized semisubmersible types and

ship shaped structures. In most of the cases, the oil producing pipes (risers) are flexible and connected to subsea completed wells. In the case of tension leg platforms (described

below), the wells may be completed at the platform deck as well.

The platforms have to be kept stationary with limited oscillatory motions motions above

the oil production wells. There are different methods to keep the platforms stationary. The

traditional way is to use a catenary anchor line system. The anchor line is then hanging in a catenary from the platform to the anchor on the sea floor, typically 1 km away at a

water depth of around 300 m. An anchor system may consist of 8-16 anchors, with anchor

lines in a star formation around the platform. Another method to keep the platform in place, is to use vertical pre-tensioned "anchor lines", or tethers, from the platform to a fixed foundation on the bottom. This is called a tension leg platform (TLP). In addition. the position keeping system may be aided by thrusters in addition to a more traditional

system. A complete dynamic position system with thrusters only is also possible. Offshore platforms can hydrodynamically be categorized as small and large volume struc-tures (Faltinsen (1990)). Viscous effects and flow separation effects are particularly impor-tant for small volume structures. Examples on large volume structures are ships, TLP and large semisubmersibles. Jackets are small volume structures. Calculations of wave induced

motions and loads an large volume structures are mostly based on potential flow. The

incident wave slope is assumed small and a perturbation scheme with the wave slope as a small parameter is used. The lowest order terms are the linear, or first order hydrodynamic loads. This means that the hydrodynamic forces and resulting motions in regular incident waves are proportional to the incident wave amplitude. A steady state solution is assumed

so that the fluid motion and the body motions oscillate with the frequency of the wave excitation loads. Results in an irregular sea state is obtained by linear superposition. In

addition to the first order forces, second (and higher) order excitation forces are generated

by the waves. These forces act at both low frequencies and high frequencies (relative to the wave frequencies). In order to compute the responses to these second order forces, the viscosity in the water has to be included. An objective of this thesis is to be able to

compute viscous damping of slow drift motions due to second order forces. However, our computations are limited to laminar flow, which has relevance for model scale, but not for

full scale.

1.1. BACKGROUND AND MOTIVATION 13

Second order loads in both high and low frequency regimes may be important for offshore floating platforms. However, the second order forces are in general so small that they only

cause problems to the structure if they excite resonance oscillations. The sum frequency

loads (second order high frequency loads) may, for example. excite vertical resonant mo-tions of a tension leg platform (TLP), while the difference frequency loads may excite the

resonant horizontal motions of a TLP, a moored ship or platform. The second order low

frequency forces can also excite resonant heave, pitch and roll motions of a conventionally moored semisubmersible platform.

Damping is important for a resonant response. A main contribution to the damping of

the motions at the natural frequencies of an offshore platform is energy dissipation due to viscous effects. In order to compute these effects we have to include the boundary layer on the surface of a structure in the numerical model. The boundary layer cannot be described

without including viscosity in the computations. This is of paramount importance, since the boundary laver is controlling the formation of vortices in the fluid, which cause the

energy dissipation and damping of the low frequency motions. The vortex shedding from

both the structure itself and the anchor lines is important for the viscous damping. The methods used in practice to compute the viscous damping are often based on empirical

drag formulations and are not sufficiently accurate.

The background for this work is found in the statement just made. The ability to predict the damping of low frequency and high frequency resonant motions is not satisfactory. Floating platforms have been used for some time. and low frequency motions is not a new "invention'. However, several reasons made them less important in the past. The traditional drilling platforms have had a relatively small volume displacement, and the

operation has been on relatively shallow waters. The vertical (heave) motions of the drilling platforms was the major concern. This is dominated by first order effects. The positioning was taken care of by a mooring system dimensioned with large safety factors. The position

has thus been controlled within the range of tolerance, with the low frequency motions present, although not predicted properly. As drilling operations and production move to deeper waters there are factors that change the picture. The stiffness in the mooring system (both catenary and tension legs) will be reduced, which is one of the effects that

causes the resonant motions to be more pronounced, and in many cases to be dominating.

This is shown in Herfjord and Nielsen (1988) where results from both experiments and

computations for a deep draft floater were reported. The deep draft floater was positioned by catenary anchor lines, and the low frequency motions were dominating in all degress of freedom. The platform was optimized so that the first order motions were minimized. The second order responses became therefore relatively more important.

Further insight in methodology and accuracy of hydrodynamic computations was gained during the research program FPS2000 (1989-92). This research programme was sponsored by the Norwegian Research Council and several norwegian and foreign oil companies. The research programmme was covering most marine technology aspects of offshore floating

production. One of the sub-projects was named "Comparative Study of Computer Pro-grams". The aim of the project was to establish a state-of-the-art for the computation of

hydrodynamic forces and motion responses of moored structures in waves. Many universi-ties, research institutions and engineering companies were invited to supply computations for a deep draft floater and a turret moored ship in an operating and a survival wave con-dition. 23 institutions from practically all over the world delivered data to the comparison

study. Not all of the institutions delivered computations of responses to the second order

forces. However, the number was considered large enough to statistically evaluate the data

and conclude about the state-of-the-art. The prediction of the first order quantities by

the different institutions showed in general very small scatter. However, there were large

scatter in the calculations of the second order forces and low frequency motions. It was concluded that the damping was the main cause for the scatter in motion computations.

80 % of the scatter was attributed to the damping, while most of the rest, about 20 %, was

due to the excitation force. A negligible part was due to the added mass. The restoring

force stiffness was specified, and had no effect on the scatter. During this investigation, the

contributors had been interviewed and some human errors were discovered. These data were removed before the statistical analysis. The hydrodynamic damping is sometimes included by using too simple approximations in standard computations of low frequency motions. The lowest level of approximation is to specify the damping as a certain per-centage of the critical damping, sometimes only- as a constant value, independent of the

sea state. This is obviously too simple. A more common way is to use a force model for

the viscous damping, i.e. a Morison type of formulation. However, the comparative study showed different opinions about what drag coefficients should be used (see also Faltinsen 1994). In calculating the drag term in Morison's equation nearly all participants assumed that the low frequency motions were uncoupled from the linear wave frequency motions.

However, this may not be a good approximation. in particular for higher sea states (see

Aanesland et al. (1990) and Herfjord and Faltinsen (1994)).

The results from the referred study were reported in Herfjord and Nielsen (1991) and (1992). The first of these publications gave a summary of the results with statistics and a short summary of methods that were used. A more elaborate investigation of the data were reported in the second paper. This included the statisical analysis mentioned above. Some recommendations for computational practice were given as well.

Another project that emphasized the lack of proper methods for viscous damping prediction

was the work with qualification of a platform of semisubmersible type

for phase II of

the Troll field. Phase II is the development for exploiting the oil region of the field.

The platform is to be positioned with a catenary anchor line system at a water depth of 330 meters. It was revealed that an accurate computation of the low frequency part

of the motions was seriously hampered by the lack of reliable damping data. The

non-linear excitation (drift forces) could be reliable computed by three-dimensional diffraction

programs, often referred to as panel programs. The slowly varying excitation force was

computed by the socalled Newman's approximation, see Newman (1974).

The damping contributions could not be properly estimated. The damping is due to a

combination of a potential flow effects (wave drift damping) and viscous effects. The wave

Li. BACKGROUND AND MOTIVATION 15

drift damping could not be computed by standard diffraction programs. One program with this option was MULDIF (Zhao and Faltinsen (1989)), a program developed at the

Norwegian Institute of Technology (NTH) and Marintek.

The viscous damping is due to vortex shedding from the platform's structural parts as well as anchor lines and risers. The effect from the anchor lines is significant, while the effect from the risers is small. The large effect from the anchor lines is due to the relatively taut mean configuration, causing a large displacement of the configuration in its own plane due

to horizontal top displacement. The top displacement is amplified several times over a large portion of the line, producing a relatively large energy dissipation. The anchor line

damping was demonstrated by model scale experiments by Huse (1986). He used a Morison

equation formulation to theoretically estimate the damping. The effect of superimposing wave frequency motions was shown in Huse and Matsumoto (1989). From model tests

with a traditional semisubmersible drilling platform referred in Huse (1986), the anchorline damping may contribute with as much as 50 % of the total, while the study by Nielsen et al. (1994) for a large semisubmersible hull indicates that the effect from the hull becomes more and more important relative to anchor line damping when the volume displacement of the hull increases. The relative contributions of the total damping could typically be 10 % from anchor lines, and more than 50 % from the viscous hull damping.

An objective of the present study is to be able to predict slow drift viscous damping in model scale of offshore platforms consisting of slender structural members. The effect of ambient wave motion, current and body motion is included. The damping will be formulated by Morison's equation, and the drag coefficient will be rationally established in stead of using empirical coefficients. One can in that way rationally account for the effect of body shape, ambient flow and body motions. This will be done by solving the Navier-Stokes equations. Existing numerical methods are adopted. The purpose of this

thesis is to test carefully a Navier-Stokes solver for relevant structures and flow situations

and try to qualify the method for the described need. The work is limited to laminar flow, while the flow is turbulent in full scale. The ultimate tool for the described need is thus not given in this thesis. However, the study is a first step in developing a reliable

method for full scale. It should be noted that two-dimensional flow in the cross-sectional

plane of a slender structure is assumed. This means that three-dimensional effects due to

uncorrelated flow along the cylinder axis cannot be studied.

It is a major concern that the method should be time efficient as well as accurate. It is

also very important that one may have confidence in the results from the method when it

is used on problems and cases where there are no similar results to compare with. Thus

it is considered important to do correlation studies on many different examples and cases

in order to qualify the method on a broad basis. Two cases which are important in ocean

engineering are given special attention. Lock-in is an important load case for risers. For the platform itself an important effect in viscous slow drift damping is the effect of interaction from superimposed wave frequency motions on the low frequency oscillations.

1.2

The organization of the thesis

The thesis is divided into three major parts. There is one part covering the theoretical considerations, including the basic theory and developments of instruments for quality

assessment and verification as well as sections containing variation of the numerical

meth-ods. This part is contained in Chapters 2 to 8. There is another part containing tests

for validation of the numerical method implemented. The validation part is contained in

Chapter 9, where a wide variation of cases are covered. Concluding remarks are given for

each case. The third major part is an example of a practical application of the numerical

method. namely on evaluating the slow drift viscous damping of a floating platform. This is reported in Chapter 10. The effect of superimposed wave frequency motions on the low frequency resonant motions is investigated particularly. Lock-in is treated in the validation part of the thesis, since this case is treated extensively by other investigators and data for

comparison are available. The thesis ends with chapters on discussion and

recommenda-tions for further work as well as conclusions. The conclusions are given on an overall basis, since each case of validation has its own concluding remark where reported.

2

The governing equations and solution methods

2.1

Introduction and discussion of methods

In the present report, a method for solving the transient two-dimensional Navier-Stokes equations with the finite element method will be presented. Laminar flow is assumed. There exist other methods to compute the viscous flow. In the following, they will be

presented shortly including a brief outline of their advantages and disadvantages. The two-dimensional discrete vortex model (see for instance Sarplcaya and Shoaff, 1979) requires relatively small computational effort. In order to predict separation points on continuously curved surfaces it is necessary to combine the discrete vortex method with boundary layer calculations. However, the predictions of the separation points are difficult. Graham (1980) also used a method of discrete vortex tracking in otherwise inviscid fluid. No viscous diffu-sion was considered. The method was used to compute the forces on sharpedged cylinders in oscillatory flow at low KC-numbers, since the separation points are clearly defined at the edges. The computations where performed on isolated edges with different internal angles,

assuming that the shedding was independent of shedding from other edges. Results were

given for a flat plate and a square cylinder, edge to the flow, showing a good comparison

with experimental results. Already Gerrard (1967) employed a discrete vortex method to

compute the flow condition behind a circular cylinder. The computations were relying on

empirical quantities and judgement from the user of the program. Single vortices were released at a plane perpendicular to the flow at a position behind the cylinder. The release position (breadth of the flow) and strength of the vortices had to be predicided according

to the Reynolds number to be simulated. The referred work does not represent the final

breakthrough of the vortex method. However, it serves as a demonstration of the potential

of the method. A single vortex model for small KC-numbers was also used by Faltinsen and Sortland (1987). The vortex sheet model (Faltinsen and Pettersem 1987) follows the roll up of sheets of concentrated vorticity and is a further development of the discrete vortex model. The method is applicable to both two-dimensional and three-dimensional

structures and can include interaction between an arbitrary number of structural elements. This method has also problems in predicting separation from continuously curved surfaces.

An advantage of the method is that it can accurately handle vortex shedding from sharp

edges. It is in particular practical for small KC number flow around bodies with sharp corners. The effect of free surface waves can relatively easily be included (see Braathen

and Faltinsen, 1988 and Faltinsen, 1993).

The two-dimensional vortex-in-cell (VIC) method follows also discrete vortices in the flow' field, but no separate boundary layer calculation is needed. Parts of the problem are solved

in a fixed grid by the finite difference method. This method may thus be characterized as a mixed Eulerian-Lagrangian method. The vortex-in-cell method solves the

vorticity-stream function formulation of the Navier-Stokes equations. A so-called operator splitting

18 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

technique is used in the solution procedure. The diffusion and advection terms are solved

in two different steps in this method. This division of the problem constitute actually an approximation of the Navier-Stokes equations. The diffusion effect is simulated by a

random walk technique. A special treatment is needed of particles that are inside the body

after the random walk. The described technique is different from the velocity correction

method given later in this report, where the advection and diffusion terms are solved in the

same step. The process of generating a grid in the numerical solution can be performed in a couple of ways. Scolan and Faltinsen (1994) used a conformal mapping technique to transform the analysed sectional shape into a circular form in an auxilliary plane. A polar mesh is used to solve the Poisson equation for the stream function in the auxilliary

plane. Scolan and Faltinsen (1993) studied conformal mapping of two circular cylinders of

arbitrary distance and position relative to each other. Slaouti and Stansby (1992) used a technique with overlapping meshes in the physical plane. In principle a general

arrange-ment of cross-sections may be handled. Graham (1988) used a vortex particle method much similar to the vortex-in-cell method. The unsteady Navier-Stokes equations in the

vortic-ity/streamfunction form was solved. It is a mixed Eulerian-Lagrangian method, where parts of the problem were solved in a fixed grid. In this mesh, the streamfunction

equa-tion as well as the diffusion part of the momentum equaequa-tion were solved. The convecequa-tion

part were modelled by convecting the vorticity particles in the velocity field. A confor-mal transformation technique was used to transform the physical mesh to a more regular

mesh, either a rectangular domain or a polar mesh. As the diffusion was solved for in the

fixed mesh, it was necessary to reinterpolate the new vorticity field from the fixed mesh

back onto the moving point distribution. This makes Graham's method different from the method described in Scolan and Faltinsen (1994). who modelled the diffusion by a random walk of the vortex particles.

The finite difference method has been used for some time to solve the Navier-Stokes equa-tions. Actually, this method has followed the whole history of Navier- Stokes solvers and has up to now been dominating in commercial codes solving viscous flow problems. Lecointe

and Piquet (1985) have given a careful study of how to use the finite difference method

to solve the two-dimensional Navier-Stokes equations in the vorticityjstreamfunction for-mulation. Yeung and Ananthakrishnan (1992) solve the Navier-Stokes equations including free surface by the finite difference method. The problem is solved in a transformed

com-putational space, where a regular grid is obtained. The numerical procedure involves a fractional-step formulation which is exactly the same approach as used in the present work, see the following section, where the different names used on the method is

com-mented as well. However, the finite difference method is considered not so flexible as the

finite element method when we are talking about complicated shapes and flow problems depending on an unstructured grid. Based on the above considerations, it is decided to use the finite element method. In this report, only two-dimensional laminar flow will be

studied, but a generalization of the method to include three-dimensional flow and turbu-lence can be done. In principle, there is no inherent limitations in the employed method for doing this. However, problems both in connection with the physical resolution of the

lam-2.2. THE NAVIER-STOKES EQUATIONS, BASIC INCOMPRESSIBLE VERSION 19

mar flow, numerical problems may occur, in particular for conditions controlled by a wake

which for relatively low Reynolds numbers will be instable, and transition into turbulent

flow will take place. According to Wile (1960) Re --= 200 is the highest Reynolds number

for which the wake of a cylinder in a steady incident flow remains laminar. An example

of free surface flow solved by the finite element method is presented by Ramaswamy and Kawahara (1987). They used the standard Galerkin procedure to solve the Navier-Stokes

equations in a Lagrangian formulation. The technique implies that the fluid in the

inte-rior of a finite element always remains in that element and fluid boundaries always move

with the element boundaries. In an incompressible Lagrangian computation. the volume of each element must remain constant. The equations are solved by a splitting-up time discretization which is very similar to the present velocity-correction method. However, some iterations are employed each time step to satisfy the equations in an updated grid

according to the Lagrangian formulation.

The finite volume method (FVM) is often referred to as something between the finite

difference method (FDM) and the finite element method (FEM). The finite volume method consider an integral form of the equations, similar to FEM. However, without a weighting

function. By treatment of the integrals, the approximations of the derivatives end up as finite differences of the unknown variables at vertices in a mesh (just as in FDM), see

Hirsch (1988).

2.2

The Navier-Stokes equations, basic incompressible version

The numerical method solves the unsteady two-dimensional Navier-Stokes equations in

laminar flow together with the continuity equation and relevant boundary conditions. The Navier-Stokes equations are in its most wellknown version written as

au

at- + (u V)u

--1Vp + v72u + g (2.1)

where it = (u, v)T is the velocity vector, t is the time variable, p is the pressure, p is the

mass density of the fluid, P is the kinematic viscosity coefficient and g is the acceleration

of gravity. The socalled del-operator V is used, in two dimensions it is defined as V =

i + 43, where i and 3 are the unit vectors in x- and y-directions, respectively. The

gravity force will be ignored in the problems discussed in this report. The continuity equation for an incompressible fluid can be written as

V u

0 (2.2)

The computational domain will be called Q. On the boundaries we have to impose

bound-ary conditions of a form that render the problem to be well-posed and solvable. The

boundary conditions are imposed according to the flow case. The -classical" case of a steady incident current passed a cylindrical section may be used as an example. In that case it is standard to define a domain with two boundaries perpendicular to the incident

20 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

flow and two boundaries parallel to the incident flow. On the upstream boundary (inlet). the boundary condition is the current velocity. This is called a Dirichlet condition for the velocity. On the downstream boundary (outlet) there is no restriction for the velocity, in order to let the created vortices leave the domain freely. On the other hand, in order to control the pressure, a Dirichlet condition is specified for the pressure on this boundary.

This condition may be obtained from the gradient of the velocity, or alternatively using a

simpler condition of zero dynamic pressure on this boundary. On the boundaries parallel

to the inflow current it is standard to impose a mixed Dirichlet/Neumann condition on the components of the velocity. It is specified zero outflow, i.e. no "vertical" flow (v=0), which

is a Dirichlet condition for this component. In addition it is specified zero gradient for

the "horizontal" velocity (i.e. = 0, which is a Neumann condition for this component,

giving a symmetry of the problem. In order to model an infinite domain, it is necessary to place the boundaries a certain distance away from the body. A more elaborate discussion

an this will be given in a later chapter. On the body we impose a no slip condition. i.e.

that the local relative velocity between the body and the flow is zero. For computations

element by element, the domain of each element is called Cr .

The presented form of the Navier-Stokes equation is given with the so-called primitive variables (velocity/pressure formulation). It is possible to introduce vorticity and stream function as variables (vorticity/streamfunction formulation). When this formulation is used, the solution is restricted to two dimensions due to the use of the streamfunction.

When the equations are formulated in the pure incompressible version. the heat exchange may be described by an uncoupled energy equation. For compressible flow, the necessary coupling is considered through relations which are assumed known, e.g. an expression for

the density (p) as a function of the pressure and temperature.

In earlier implementations of incompressible flow solvers, the incompressibility constant was

imposed in an incomplete manner through a socalled penalty parameter (see e.g. Brooks

and Hughes (1982)). The pressure was then expressed as

p=

(2.3)

A is the penalty parameter, and is set to a large number. By this, the pressure (p) become finite, even if V is is very small (should be identical to zero in incompressible flow). Since

A is a predecided known number, the result is that one unknown is eliminated, and the continuity equation is imposed implicitely in the momentum equation by Eq. (2.3). It is emphasized that the present implementation is based on Equations (2.1) and (2.2).

Remark The equations are implemented in dimensional form. The default value for the

density is the sea-water value 1025 kg/m3. The kinematic viscosity (Li) is then expected to be an input of unit m2/s. Accordingly, the length and velocity are assumed to have dimensions (in) and (m1s), respectively. More or less as a standard, the characteristic

length dimension (e.g. the diameter D of a cylinder) is 1 measured in meters, and the

characteristic velocity (e.g. the inflow steady velocity (U) or the ambient oscillatingflow

2.2. THE NAVIER-STOKES EQUATIONS, BASIC INCOMPRESSIBLE VERSION 91

non-dimensional form by using a density equal to 1. The viscosity parameter is in that

case non-dimensional. being equal to I, = -`

2.2.1

The introduction of accelerated coordinate system

The simulation of a moving body in a flow field may be done in a couple of alternative ways. The closest to think of may be to move the body relative to the boundaries in the meshed domain. This would imply a modification of the mesh according to the motion. A much more efficient way to obtain the same effect is to introduce an accelerated coordinate system fixed to the body. This method is implying that the computational grid is following the body. The computations are performed in the moving grid. However, the solution may

be referred to the earth-fixed system. Let us refer to the global (earth-fixed) system by

the subscript g, and the local (moving) system by 1, then

dxi

x, = xi +

uB(t)dt.

= =us

0 di

where the x-coordinate is considered, and us is the body velocity in x-direction.

= u(x,,y,,t) = u(ri +

us(i)di,y( +

f vB(t)dt,t)

(2.5)

The spatial differenciation (x-direction) is expressed as follows:

att(xg, ,t) au(x, + fr; uB(t)dt, yr + vB(t)dt,i)

ax,

ax, ax

art

Showing that the spatial differentiation remain unchanged. The temporal differentiation

is expressed as follows:

au(xg, y9, t) au(xi + uB(t)dt. yi fr: vg(t)dt,t)

at

at

ou(si + fot us(i)dt.yi + vB(t)dt. I) agi On

at

au au au

=

-uB(t)--a-x-e -

vs(-The corresponding expression for the v-velocity is:

ay(x,,y,t)(9v

ay ay

at _ yE(t)5,- (2.8)

The above expressions give the necessary transformations between the earth-fixed and the moving coordinate systems. By this it is possible to express the Navier-Stokes equations

solved in a moving coordinate system, however, referring the solution to the earth-fixed

system as follows: (2.4) (2.6) (2.7) 0 + at

22 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

au

(us V)u + (at V)u =-

Vp + vV2u + g

0.9)

at

where us r- (us,v8)T is the velocity vector of the body relative to an earth-fixed coordinate system. it = (u,v)T means the flow velocity vector referred to the earth-fixed coordinate

system. means the rate of change with time for a fixed point in the accelerated coordinate system. It is emphasized that Eq. (2.9) is describing the case of a translatorv (non-rotating)

system. Further it is commented that the given equation describe the case of a moving body in infinite fluid domain. In the case of motion relative to another body, a wall or a

free surface, it will be necessary to modify the mesh.

2.3

The finite element method

'2c3.1

Basis and background

In the finite difference method, the equations are attacked by approximating the derivatives

directly by using the unknown variables defined at nodes in a grid. In the finite element

method it is performed a weighting integration over the domain, followed by partial

inte-gration which reduce the higher order derivatives by (at least) one. This makes the two methods to appear rather different, and to make the FEM more appealing, since lower order derivatives are easier to approximate accurately than those of higher order. The above brief introduction to the principles of the finite element method will be elaborated

in the following.

The basics of the finite element method is explained below, to a large extent following Zienkiewicz and Morgan (1983), by starting with a differential equation (DE) written in

general form.

.Cu + p = '0 (2.10)

In the outset it is assumed that is a linear differential operator and p is independent

of u. The unknown function it is defined on a domain II which may be one-, two- or

three-dimensional. The boundary of this domain is denoted T. On this boundary there are boundary conditions for the dependent variable (u). On one part called 1'9,it have known

values, g. We call this a Dirichlet boundary condition. On the rest of r, named Th, there

are derivative conditions on u, e.g. = h. This type of condition is named Neumann

boundary condition.

In the process of solving the posed problem, a weighted residual method is applied. First of all we need to construct socalled trial functions as candidate solutions to the problem. Next we define socalled weighting functions to be applied. The trial functions are defined so that they satisfy the Dirichlet condition on 1.19. The weighting functions

are defined to take the value zero at the same part of the boundary. The notations ü for

2.3. THE FINITE ELEMENT METHOD 93

solution ü into the original DE we will be left with a result different from zero, contrary to the case with the exact solution inserted. This residual we call Rn:

= + p (2.11)

The idea is now that by the use of the weighting function we may level out the residual

and make it zero in a weighted manner by integrating over the domain as follows:

fzoRndS1 = 0 (2.12)

This is by the present author considered as the basic idea behind the finite element method.

It is hard to think of trial- and weighting functions valid for the whole domain in any

problem. Consequently, we divide the domain into pieces and define functions piecewise

over the domain. The further development is elaborated in Appendix B. There the finite

element concept is described by solving a linear one-dimensional second order differential

equation. The solution of the Navier-Stokes equations by the finite element method is described in some detail. In the following, a more brief description of the present imple-mentation is given, only slightly modified compared to the presentation in Herfjord and

Faltinsen (1994).

2.3.2

Spatial discretization

The element method is treated extensively in many text books. In the present context the

books by Chung (1978), Zienkiewicz and Morgan (1983), Hughes (1987) and Zienkiewicz

and Taylor (1989, 1991) are found particularly useful. lathe following text some details

about the present implementation of the finite element method is given.

The spatial domain is divided into a large number of elements which may vary in size.

The element corners (common with the neighbouring elements), is called nodes (we may also choose to use internal nodes in addition). The unknowns are defined on these nodes. The element method involves using a socalled trial solution of an unknown function in the problem. The trial solution is made of a function of specified form which is interpolating the solution between the defined nodes. The interpolation function consists of basis functions existing in a number equal to the number of nodes. The basis functions are defined to have

the function value 1 at 'its own' node, and zero at all other nodes. The trial solution for the unknown function u (e.g. the velocity in x direction), may be written as below. As it

is an interpolated approximation, we call it "Ilh

repo

Uh = E

(2.13)

where Ni are the basis functions and npo= number of nodes. u i = 1,

. ,npo are

the function values at the nodes, and these will constitute the unknowns in the problem.

94 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

known values. These are the boundary conditions on the part of the boundary that have

prescribed values, socalled Dirichlet boundary conditions.

We now need to create equations in a number corresponding to the number of unknown

function values for each variable. There may be several variables (degrees of freedom) at each node. In a typical problem in fluid dynamics with incompressible flow, we solve for the velocity components and the pressure at each node. When we restrict ourselves to two dimensions, it means that we have three unknowns at each node.

In order to create the necessary number of equations we use the method of weighted resid-uals, which introduces still a new function, the weighting function. Wh. In many problems we prefer to create the weighting function by the basis functions already defined, and we may write

npo

h

EciAri

_(2.14)

i.1

We have described the function with npo basis functions. As mentioned above, some of the

npo function values are known, and do not need an equation in the system to be created. However, a standard procedure in 'real life' is to build up the total system with a number of npo equations. Then the system is treated according to the number of nodes with Dirichlet

boundary conditions. One technique involves eliminating rows and columns and leave a system of exactly the number of unknowns. In another technique, the rows and columns

corresponding to the known values are modified, leaving the original size of the system to

be treated by the solver. In the present implementation, the latter alternative is used (for

reference, see e.g. Dhatt and Touzot. 1984).

Using the same basis functions for the trial solution and the weighting function, is known as the Calerkin method. Fluid flow is actually one example where the weighting function sometimes is made of different basis functions (different from the trial solution). However, in the present implementation we use the Galerkin method.

When we insert the trial solution into the differential equations that we want to solve, we get a residual (or error), since it is not completely satisfying the equation. We seek

to nullify this residual in a weighted manner by multiplying with the weighting function and integrate over the domain we have defined. We want this operation to be valid for arbitrary coefficients c, thus we get as many equations as there are unknowns. In reality

the integrations are performed element by element, and finally an assembling is made into

the equation system that we want to solve. 2.3.3

Choice of basis functions

As mentioned earlier the spatial domain is divided into a large number of elements which

may vary in size. They have to be small in the areas where there are large gradients in

2.3. THE FINITE ELEMENT METHOD 25

the present implementation. we have chosen to use triangular elements. We use linear basis functions on these triangles. This means that between the nodes, the solution is

interpolated by a linear function.

Linear triangular elements are easy to implement and the computations are fast on each element. However, compared to a bilinear quadrilateral element, the number of elements

will be exactly doubled, so the gain in element by element computational effort is more or less eliminated by the larger number. However, when there is a need to employ extremely

unstructured grids in connection with complicated boundaries and automatic remeshing is used, it seems that triangular elements are preferred, see e.g. Peraire et al. (1990) and

Moller and Hansbo (1995).

When we insert the trial solution into the differential equation and perform the integrations

over the spatial domain, we eliminate the partial differentiation in space and end up with an algebraic equation when we deal with a stationary problem.

In the case of a time

dependent problem, we end up with ordinary differential equations in time. It is also possible to use a finite element approach in the time integration as well. This method is

described in e.g. (Zienkiewics and Taylor, 1991: Hughes, 1992; Mittal and Tezduyar, 1992; Johnson, 1992).

2.3.4

Time integration

We have chosen to use the forward Euler method to perform the time integration. This is actually the simplest possible method for performing a numerical integration. It is

a first order method that demands small time steps in order to be able to march the

solution forward in time in a stable manner and within acceptable accuracy limits. The

method is, however, in relatively common use (see Kovacs and Kawahara, 1991 and Gresho

et al. 1984). The method is discussed thoroughly in the latter publication. We use a procedure called the velocity correction method in the time integration. This is based on a split operator technique already proposed by Chorin (1968). The split operator method

is presented in Zienkiewicz and Taylor (1991) in a couple of variants. It is used by Donea

et al. (1982) who called it the fractional step method. If a turbulence model is used, the

velocity-pressure equations may be solved by use of the same technique.

The forward Euler method used on Navier-Stokes equations is written as follows un-" 1 vpn-1-1 vv2un (urli

v)un

(tin v )un

(2.15)

At

where the superscript n denotes time step number. One should note that we have left out

the gravity and kept the formulation with the accelerated coordinate system (see equation

(2.9)). The updated pressure (p'1) appears on the right hand side of Eq. (2.15). As will be shown shortly, the computation provide the updated velocity/pressure field at the end

of the computational procedure each time step. The updated velocity may now be written

-96 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

un+i un

_1vpn+i

pv2,10 (uri/3

v)un

(un 77)u)

_(2.16)

We see that the updated pressure term is on the right hand side of the equation. This is

of course still unknown. and in the first step of the velocity correction method we simply

forget about the pressure term and write an intermediate velocity as

U = un

At(vVun

. V)un + (un V)un)

The lacking part of Eq. (2.16) reads

un+i

_vpn+i

If this equation is substituted into Eq. (2.2), we get v2pn+1=

Pv.0

At

There is no term with un+' in Eq. 2.19, because this is the correct velocity at step n + 1, which satisfies the continuity equation. The equation above is a Poisson equation for the pressure. The last step is a velocity correction step, giving un+i by Eq. (2.18). Thus, the solution is performed in the three step given by equations (2.17), (2.19) and (2.18)

respectively.

2.3.5

Matrix equations

By approximating the time differentiation by the forward Euler scheme, we are left with the spatial differentiations. In the further treatment we use trial solutions and weighting function (as mentioned earlier) in order to obtain a weighting residual form giving the wanted matrix formulation. Trial solution based on linear interpolating basis functions is

inserted on each step. We write the functions as follows

"pa uh(x,

=

E ui(i)Ni(x,y)

(2.20) "pa

vh(x,y,t) = E vi(t)Ni(x,y)

(2.21) t=i "pa Ph(x,Y-t)=

EpicoNi(x,

y) (2.22) i=1

We will comment at this point that we use equal order basis functions for the velocities and the pressure. Others have used a basis function for the pressure of one order lower than for the velocities (Brooks and Hughes, 1982). This method imposes a stabilizing effect on the pressure. which otherwise tends to experience spurious oscillations. However, equal order interpolation is used by e.g. Tezduyar et al. (1992), who are using a stabilizing term on the

(2.17) (2.18) (2.19)

=

-y,

step 3

At

Mu"*"

_Dpn-1-1

Note that in these equations, u and p denote the vectors of unknowns at the nodes, (uT =

, unpo. v1... v,,,o), PT = (pl... pnpo)). The matrices involved in the equations are

given below as an element by element integration. The so-called mass matrix is given as

follows

M' =

f NNTdit

(2.26)

where N = [N1 N2 N3]T is the vector of the basis functions. As a triangular element with linear variation is chosen, the basis functions exist in a number of three per element. The

supersript (e) denotes that the integration is performed over each element. The global interactions are obtained by summing (assembling) the element contributions into the

system matrix.

In the diffusion matrix ,573nTv. we have included the balancing tensor viscosity:

(2.25)

where in and I.7n are the arithmetic mean of the velocities in x and y direction, respectively.

The velocities are updated each time step, thus the diffusion matrix is varying with time,

sn,v

[N,, N,y1 , -n - Un 2 At-77, 71 2 -n -n 7,' 2 + tf7nVn N T.L. N:Tv (2.27)

2.3. THE FINITE ELEMENT METHOD 27

pressure. We do not use any special treatment of the pressure in our implementation. By

using small time steps we have not experienced any unphysical behaviour in the pressure

field.

Matrix equations are obtained by introducing a weighting function in a manner attributed

to the russian mathematician Galerkin. The Galerkin method in combination with

ex-plicit integration (Euler) is known to introduce an effect of negative diffusivity (Gresho et a/.,1984; Brooks and Hughes, 1982). In order to counteract this negative viscosity, we add viscosity by using the so-called balancing tensor viscosity (Kovacs and Kawahara, 1991;

Gresho at al., 1984). This is a concept working on the diffusion term, and the modified diffusion matrix is given the subscript BTV. The matrix equations that evolve from the

technique, are given as follows

step 1

Mu = Mu n

At.5237-vu'l + At(unB.D Dy)un AtAnun + Ata-6,) (2.23)

step 2 Sa)11+1 =-1±DT + Gr(`74-1 (2.24) At

=

= (a11...

28 CHAPTER 2. THE GOVERNING EQUATIONS AND SOLUTION METHODS

and we denote it by n as superscript. The comma (,) followed by x as subscript means

differentiation by x.

The advective term of the equation is nonlinear, since it is a function of the velocity, and

is updated according to the solved velocities at each time step. The advection matrix is

written as follows:

A" =

L1V 7.0T N NT +

I

NeT N NT d1-1_,x (2.28)

1/`

The pressure term in step 2 has a matrix similar to the unmodified diffusion matrix,

however, without the viscosity coefficient. We prefer to call Sic, the Laplacian matrix, since is has its origin from a Laplacian operator.

S;

=

f

Nl\TTrdO + N,,,N,Tyclit (2.29)

Further DT = (D,, po is called the divergens matrix (from the divergens operator),

where

D,

=

f

N .N,Trdct, Dy =

f

N N.TyclO (2.30)

t

D is the gradient matrix, given as

D =

{Dx

Dy } (2.31)

The terms Gn) and Gn+1 are boundary terms evolving from an integration by parts of

Ku (p)

terms including second order differentiation in space.

2.3.6

Solution technique

In the solution of the system of equations, the procedure is simplified by applying a diagonalization of the mass matrix. Different diagonalization methods are discussed in

Zienkiewicz and Taylor (1989) [Appendix 8] and Zienkiewicz and Taylor (1991) [Chapter 9]. For a linear triangular element, all the alternatives gave equal results, the diagonal being equal to the sum along each row of the matrix, which may be written as follows:

mz,

E

f

NiNidS1 (2.32)

This technique is often referred to as mass lumping. an expression coming from structural

engineering, where the technique simply means that the mass of the structure is

concen-trated at nodes. This technique offer a large simplification of the solution, since the system

is decoupled, and there is one unknown per equation in the system. Further reference is made to Chapter 5. where computations also with consistent mass matrix are performed

2.3. THE FINITE ELEMENT METHOD 29

The mass lumping technique is used in step 1 and step 3. In step 2, where the pressure

is solved for, we do not simplify the matrix. However, we compress the storage by using the skyline storage technique, as described in i.e. Dhatt and Touzot (1984). The procedure

for the storing and the solution method (actually a Gaussian elimination method). is in

fact taken from the referred text book. A further reference in connection with the present implementation is Ren and Utnes (1993).

3

Convergence, consistency, stability

3.1

Introduction and statements on convergence and consistency

When a physical problem described by a differential equation is solved numerically, we want to know the -quality" of the solution, i.e. how far from the correct solution our

results are. And secondly we want to know how fast the results approach the correct value as we let the discretization size approach zero. These questions lead directly to the notion

of convergence and order of accuracy or rate of convergence.

In traditional numerical treatment of differential equations, as in solution by the finite

difference method, the question about convergence is treated indirectly, since convergence may be hard to prove directly. So, the question of convergence is treated indirectly by the

notions of consistency and stability. Consistency is about the ability of the numerical method to represent the original differential equation properly. The element formulation

is consistent if it represents the exact DE when the discretization size goes to zero. When a time-dependent problem is solved, stability means in addition that the time integration method is able to march the solution forward in time without unbounded growth of the

initial condition.

The indirect approach to convergence is ensured through Lax's equivalence theorem

(reproduced from Smith (1990)):

Given a properly posed linear initial-value problem and a linear finite-difference

ap-proximation to it that satisfies the consistency condition, stability is the necessary

and sufficient condition for convergence.

As seen, the theorem is stated for linear problems formulated in the finite difference method, but it is often referred to in connection with finite element analysis as well (see e..g Hughes

(1987) [Chapter 8]. Within the finite element method the problem is given an integral

formulation, resulting in a system of equations to be solved where the coefficients are results from integrations. The consistency of this equation to the original DE is not obvious. For the element method, the so-called Babuska-Brezzi condition is often referred to as a test of convergence. Babuska (1971. 1973) and Brezzi (1974) developed expressions for error

bounds for linear differential equations solved by the finite element method. The error

bounds were given in general terms, however, giving a rate of convergence depending on

the order of the approximation (i.e. interpolation functions). In Zienkiewicz etal. (1986) and Hughes et al. (1986) the Babuska-Brezzi condition is referred to in connection with the solution of a linear steady problem in fluid flow (containing only pressure and diffusion,

and no advection). In fluid flow the unknowns are velocities and pressure, and the FEM 30

3.1. INTRODUCTION AND STATEMENTS ON CONVERGENCE AND CONSISTENCY 31

formulation constitute a mixed formulation. Hughes et al. (1986) points out that

equal-order interpolations for velocities and pressure in connection with a pure Galerkin method

fail to satisfy the Babuska-Brezzi condition. The combination equal order interpolation with Galerkin is told mainly to have a destabilizing effect. causing oscillations in the

solution. In Zienkiewicz et al. (1986) the patch test developed by Irons is discussed and stated to be a simplified method for testing the satisfaction of the Babuska-Brezzi condition. The patch test is used to disqualify equal order interpolations on a mixed

formulation, (at least for a steady. linear problem), when standard Galerkin is used. Hughes et al. (1986) report a convergent solution using equal order interpolations on a stationary mixed problem by using the so-called SUPG (streamline-upwind/Petrov-Galerkin). SUPG is using a stabilizing attachment to all terms in the equations. Formal proof of convergence is carried out. Through the documentation of both consistency and stability, an expression

is obtained showing convergence and order of accuracy.

By this it is stated that the

Babuska-Brezzi condition is circumvented.

In Hughes et al. (1987) convergence analysis is performed on linear time-dependent multi-dimensional advective-diffusive systems. which may be called a linear version of the Navier-Stokes equations. The full non-linear case cannot be treated formally. The referred article

used a socalled space-time discretization of the problem by the element method. The er-ror estimate was then deduced as depending on a common space-time mesh parameter. Tezduyar et al. (1992) describe a similar method using SUPG and in addition a pressure

stabilizing attachment referred to as PSPG (pressure-stabilizing/Petrov-Galerkin).

An example of an error bound for an element method may be given with reference to

Johnson ( 1990). The proof is carried through for the case of a linear interpolation function over a triangular element. The result of the limiting error is given as follows:

492u - Uh < 2h2 I mar

19X aY

where

u=the unknown function

uh=a linear interpolant of u

h=mesh parameter (e.g. the largest side length of the triangle)

value over the triangular element

As an absolute value for the error it is of course of no use, since one does not know neither

the function nor the derivative of it. However, the equation shows that the error goes to zero as the mesh size goes to zero (convergence) and that the rate of convergence is h2. Which is a rather encouraging result, since only a linear interpolation function is used. This means that if the mesh size is halved, the error is reduced to one quarter. Another assumption behind this finding is that no side of the element is very small compared to

the others. If so, the rate of convergence will slow down.