I hereby declare that the work presented in this thesis is solely my own work and ..... A finite volume discretisation is used, where the domain being studied is divided into a ..... The steady transport equation is obtained by dropping the transient term ...... register. Щd memory. (a). (b). (c). Figure 6.3: Memory hierarchy for three ...
A Parallel Navier–Stokes Solver for Natural Convection and Free Surface Flow A Thesis submitted in partial fulfilment of the requirements for the Degree of Doctor of Philosophy by S.E. Norris
Department of Mechanical Engineering University of Sydney September 2000
Abstract A parallel numerical method has been implemented for solving the Navier–Stokes equations on Cartesian and non-orthogonal meshes. To ensure the accuracy of the code first, second and third order differencing schemes, with and without flux-limiters, have been implemented and tested. The most computationally expensive task in the code is the solution of linear equations, and a number of linear solvers have been tested to determine the most efficient. Krylov space, incomplete factorisation, and other iterative and direct solvers from the literature have been implemented, and have been compared with a novel black-box multigrid linear solver that has been developed both as a solver and as a preconditioner for the Krylov space methods. To further reduce execution time the code was parallelised, after a series of experiments comparing the suitability of different parallelisation techniques and computer architectures for the Navier–Stokes solver. The code has been applied to the solution of two classes of problem. Two natural convection flows were studied, with an initial study of two dimensional Rayleigh–B´enard convection being followed by a study of a transient three dimensional flow, in both cases the results being compared with experiment. The second class of problems modelled were free surface flows. A two dimensional free surface driven cavity, and a two dimensional flume flow were modelled, the latter being compared with analytic theory. Finally a three dimensional ship flow was modelled, with the flow about a Wigley hull being simulated for a range of Reynolds and Froude numbers.
i
Declaration I hereby declare that the work presented in this thesis is solely my own work and that to the best of my knowledge the work is original except where otherwise indicated by reference to other authors. No part of this work has been submitted for any other degree or diploma.
Stuart Norris September 2000
ii
Acknowledgements Firstly, I would like to thank my supervisor, Steve Armfield, without whom no thesis would have eventuated. His support has always been gratefully appreciated, and his knowledge and scholarship has proved invaluable. The Department of Mechanical Engineering at the University of Sydney, for their financial support with an APA and travel scholarships. The School of Mathematics at the University of New South Wales, and especially Bill McKee, for their support in the first year of my PhD. Mark Sagar, Paul Charette, Peter Hunter, and the PTM/University of Auckland collaboration for providing lots of interesting and rewarding distractions. Kwok, Sandy, Chen and Dugald for many interesting discussions. Peter Jackson and Gordon Mallinson who are partly to blame for my interest in Fluid Mechanics and Numerical Modelling. Saskia and May for being good flatmates. Brian Maher and Neville Turner, for some practical work on Port Jackson. Long may Starlight Express vanquish the opposition. My Parents, for their never ending support. All of which is nothing compared to that due to Fiona.
iii
Contents Abstract
i
Declaration
ii
Acknowledgements
iii
List of Symbols
viii
Commonly Used Acronyms
ix
1
Introduction
1
1.1
1
2
An Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Differencing Schemes
7
2.1
Equations and Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2
The Steady Transport Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2.1
Numerical Stability Issues . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.2.2
First Order Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
2.2.3
Second Order Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2.2.4
Third Order Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.2.5
Summary of the Advective Discretisation Schemes . . . . . . . . . . . . . .
26
2.3
Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
2.4
The Transient Transport Equation . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
2.4.1
The Forward Euler Differencing Scheme . . . . . . . . . . . . . . . . . . .
30
2.4.2
The Backward Euler Differencing Scheme . . . . . . . . . . . . . . . . . .
31
2.4.3
Crank–Nicolson or Centred Differencing . . . . . . . . . . . . . . . . . . .
33
2.4.4
Adams–Bashforth Differencing . . . . . . . . . . . . . . . . . . . . . . . .
34
iv
CONTENTS
2.4.5 2.5
2.6 3
Summary of the Temporal Discretisation Schemes . . . . . . . . . . . . . .
35
A Comparison of the Discretisation Methods . . . . . . . . . . . . . . . . . . . . .
36
2.5.1
The Advecting Witches Hat Problem . . . . . . . . . . . . . . . . . . . . .
36
2.5.2
The Smith–Hutton Problem . . . . . . . . . . . . . . . . . . . . . . . . . .
47
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
Linear Solvers
52
3.1
A Description of the Linear Solvers . . . . . . . . . . . . . . . . . . . . . . . . . .
52
3.1.1
Linear Equations Resulting from a Finite Volume Discretisation of a PDE . .
53
3.1.2
Direct Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
3.1.3
Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
3.1.4
Simple Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
3.1.5
Incomplete Factorisation Methods . . . . . . . . . . . . . . . . . . . . . . .
61
3.1.6
Krylov Space Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
3.1.7
Multigrid Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
A Comparison of the Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
3.2.1
The Solver Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
3.2.2
Convergence of the Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
3.2.3
Scaling of the Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
3.2.4
Memory Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
91
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
92
3.2
3.3 4
v
Solution of the Navier–Stokes Equations
93
4.1
The Navier–Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
93
4.2
Solution of the Navier–Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . .
95
4.2.1
The SIMPLE Velocity–Pressure Coupling Scheme . . . . . . . . . . . . . .
95
4.2.2
Rhie–Chow Velocity Interpolation . . . . . . . . . . . . . . . . . . . . . . .
99
4.2.3
A Fractional-Step Method for Transient Flow . . . . . . . . . . . . . . . . . 100
4.3
Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.4
Benchmark Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.5
4.4.1
The Driven Cavity Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.4.2
The Natural Convection Flow . . . . . . . . . . . . . . . . . . . . . . . . . 109
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
CONTENTS
5
Non-Orthogonal Meshes
Non-Orthogonal Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.2
Non-Orthogonal Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.3
Non-Orthogonal Advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.4
Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.5
Navier–Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.7
5.5.1
The Pressure Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.5.2
Velocity Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.5.3
Pressure Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
Examples of Non-Orthogonal Problems . . . . . . . . . . . . . . . . . . . . . . . . 133 5.6.1
Conduction in a Skew Domain . . . . . . . . . . . . . . . . . . . . . . . . . 133
5.6.2
Driven Cavity Flow with a Distorted Mesh . . . . . . . . . . . . . . . . . . 138
5.6.3
Driven Cavity Flow in a Skew Cavity . . . . . . . . . . . . . . . . . . . . . 143
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
Parallelisation and Architecture
147
6.1
The Finite Volume Solver Test Code . . . . . . . . . . . . . . . . . . . . . . . . . . 147
6.2
Architecture and Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
6.3
7
115
5.1
5.6
6
vi
6.2.1
Architecture and Program Speed . . . . . . . . . . . . . . . . . . . . . . . . 149
6.2.2
Implementation and Program Speed . . . . . . . . . . . . . . . . . . . . . . 153
6.2.3
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
Parallelisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 6.3.1
Parallelising the Test Code . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
6.3.2
Parallel Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
6.3.3
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
6.4
Parallelisation of CFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
6.5
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
Natural Convection 7.1
171
Rayleigh–B´enard Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 7.1.1
Using CFD to Measure the Critical Rayleigh Number . . . . . . . . . . . . . 173
7.1.2
Using CFD to Model Wavelength Selection . . . . . . . . . . . . . . . . . . 178
7.1.3
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
CONTENTS
7.2
Three Dimensional Transient Convective Flow . . . . . . . . . . . . . . . . . . . . . 181 7.2.1
7.3 8
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
Viscous Free Surface Flow 8.1
195
The Modelling of Flow About Ships . . . . . . . . . . . . . . . . . . . . . . . . . . 195 8.1.1
A Review of CFD Schemes Described in the Literature . . . . . . . . . . . . 196
8.2
A Numerical Method to Model Free Surface Flow . . . . . . . . . . . . . . . . . . . 198
8.3
Examples of Free Surface Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200
8.4 9
vii
8.3.1
The Free Surface Driven Cavity . . . . . . . . . . . . . . . . . . . . . . . . 201
8.3.2
Two-Dimensional Channel Flow . . . . . . . . . . . . . . . . . . . . . . . . 202
8.3.3
Free Surface Ship Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215
Conclusions
216
List of Symbols
!
# $
%'& ,
.
-
/10 / 6 7
8'9 80
* @ A C E
" : H
3I J
K N
M
G
L
Courant number constant pressure specific heat diffusion flux Fourier number Froude number acceleration due to gravity specific enthalpy thermal conductivity mass flux mass fraction of species " length number of equations number of processors Nusselt number () +* order multigrid prolongation operator pressure 2 43 P´eclet number 5" Prandtl number multigrid restriction operator residual Rayleigh number :)+*; ?" Reynolds number 2 5 temperature velocity B axis component of velocity D axis component of velocity F axis component of velocity
G thermal diffusivity coefficient of volumetric expansion for temperature error diffusivity wavelength dynamic viscosity J G kinematic viscosity volume streamfunction density non-dimensional time dimensionless temperature *; >+* viii
Commonly Used Acronyms ADI BiCG BiCGSTAB CFD CG CGS CMSSL FLOPS FOU GMRES HPF IC ICCG ILU K&R LDL LU MAC MG MPI MSI MSOU PDE PVM QMR QUICK RBSOR SAXPY SIMPLE SD SIP SOR SSOR SOU ULTRA VOF
Alternating Direction Implicit (linear solver) Bi-Conjugate Gradient (linear solver) Bi-Conjugate Gradient Stabilised (linear solver) Computational Fluid Dynamics Conjugate Gradient (linear solver) Conjugate Gradient Squared (linear solver) Connection Machine System Scientific Library (numerics library) Floating point Operations Per Second First Order Upwind (discretisation scheme) General Minimalised Residual (linear solver) High Performance Fortran Incomplete Cholesky factorisation (linear solver) Incomplete Cholesky–Conjugate Gradient (linear solver) Incomplete Lower–Upper factorisation (linear solver) Kernighan and Ritchie (the original dialect of C) Lower Diagonal Lower O factorisation (linear solver) Lower Upper factorisation (linear solver) Marker And Cell Multi-Grid (linear solver) Message Passing Interface (parallelisation library) Modified Strongly Implicit procedure (linear solver) Monotonic Second Order Upwind (flux-limited discretisation scheme) Partial Differential Equation Parallel Virtual Machine (parallelisation library) Quasi-Minimalised Residual (linear solver) Quadratic Upwind Interpolation for Convective Kinematics (discretisation scheme) Red Black Successive Over Relaxation (linear solver) A Scalar PRQ BTSUD operation Semi-Implicit Method for Pressure-Linked Equations Steepest Descent (linear solver) Strongly Implicit Procedure, Stone’s method (linear solver) Successive Over Relaxation (linear solver) Symmetric Successive Over Relaxation (linear solver) Second Order Upwind (discretisation scheme) Universal Limiter for Tight Resolution and Accuracy (flux-limiter) Volume Of Fluid
ix
Chapter 1
Introduction The modelling of fluid flows is of great interest to Engineers and Scientists alike, with many engineering problems and issues of scientific interest depending upon complex flow phenomena. Due to the complexity of the flow problems there are few analytic models of fluid flows, but the advent of digital computers has stimulated the development of numerical methods for the modelling of flow. This thesis describes the development of a computer code to model the Navier–Stokes equations, the equations of motion for a viscous incompressible fluid. A finite volume primitive variable formulation is used, with particular care being paid to the accuracy and the computational speed of the method. The code has been tested at each stage of it’s development, with the developed code being used to model two gravity dominated flows; natural convection in an enclosure, and free surface flow.
1.1
An Outline of the Thesis
This thesis discusses the development of a computer code for modelling fluid mechanics problems. The code developed is sufficiently general to allow the modelling of a wide class of problems, and care is taken in the development to ensure the accuracy and efficiency of the numerical methods used. After the description and testing of the code, it is applied to the modelling of two problems in fluid mechanics. Chapter 2 describes the finite volume discretisation of the transport equation, an equation for the conservation of some property (such as enthalpy or linear momentum) in the presence of diffusion and fluid flow. A finite volume discretisation is used, where the domain being studied is divided into a number of discrete cells, on each of which a conservation equation is imposed. The technique ensures global conservation, an advantage over finite difference schemes, whilst being simpler to derive and implement than finite element schemes. The initial discretisation is derived on a Cartesian mesh, where boundaries of cells are aligned with the Cartesian axes. The diffusive terms in the transport equation are modelled using a central difference approximation, while a number of different discretisations are implemented and tested for the advection and transient terms. For the advection terms, which quantify transport due to motion of the fluid, discretisations based upon first, second and third order approximations are described, along with two advanced techniques that use flux limiters to prevent the occurrence of physically unrealistic oscillations in the solution. For the time discretisation first and second order methods are described and implemented. To compare the accuracy of the different discretisation schemes they were used to model two standard
1
CHAPTER 1. INTRODUCTION
2
Exact
1
Central
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
-0.2
-0.2
Figure 1.1: The final scalar distribution for the transient advection test problem. The exact solution (left), and the solution calculated using a central difference approximation (right). test problems, the first being the transient advection of a scalar blob in the absence of diffusion shown in Figure 1.1[168], whilst the second problem was the steady advection and diffusion of a scalar field shown in Figure 1.2[157]. The implemented discretisation schemes were tested for a range of mesh sizes and time steps, with the calculated solutions compared with exact solutions. The finite volume discretisations used in this study generate large sparse sets of linear equations, and in Chapter 3 methods available for the solution of these linear systems are discussed. The two main classes of linear solver described are direct and iterative methods. Direct methods solve a system of linear equations in a set number of steps, with examples of such solvers being dense LU decomposition, and tridiagonal and block tridiagonal methods that take advantage of the sparse nature of the linear systems generated by finite volume discretisations. In contrast to the direct solvers, iterative methods solve a system in an undetermined number of operations, with the solution operation being repeatedly applied to an initial estimate of the solution, and with the error in the estimated solution being measured. When the error is reduced below some arbitrary tolerance the equations are considered solved and the process terminates. Whilst the number of operations is potentially unlimited, the methods are typically faster and require less memory than direct methods.
2.5 2
1 1.5 1 0.5
0.5
0 1 0.8 0.6
-1
0 -1
-0.5
0
0.5
1
-0.5
0.4 0
0.2
0.5 1 0
Figure 1.2: The Smith–Hutton steady advection–diffusion problem. The streamfunction for the flow (left), and the solution for the scalar distribution (right). Following the taxonomy of Barrett et al[7], the iterative methods described in this thesis are classified into four categories; the simple iterative methods, such as Jacobian Iteration and Successive Over Relaxation, incomplete factorisation methods, such as Incomplete LU factorisation, Krylov space meth-
CHAPTER 1. INTRODUCTION
3
2D/Dirichlet: Convergance of Iterative Solvers
2D/Dirichlet: Solution Time for Solvers
1
100 Jacobi MSI MG-Jacobi MG-SIP CG CG-Jacobi CG-MSI CG-MG-SIP BiCGSTAB BiCGSTAB-Jacobi BiCGSTAB-MSI BiCGSTAB-MG-SIP GMRES-MSI GMRES-MG-SIP
Residual
0.01
0.001
10
1
Solution Time
0.1
0.1 LDL Decomposition Block Tridiagonal/LDL Jacobi MSI MG-Jacobi MG-SIP CG-MSI CG-MG-SIP BiCGSTAB-MSI BiCGSTAB-MG-SIP GMRES-MSI GMRES-MG-SIP
0.01
0.0001 0.001
1e-05
0.0001
1e-06
1e-05 0
5
10
15
10
100
Time
1000 10000 Number of Equations
100000
1e+06
Figure 1.3: The rate of convergence, and the scaling of the solution time with the number of equations, for linear methods solving a two dimensional discretisation of the Laplace equation. ods (sometimes referred to as Conjugate Gradient methods), and multigrid (or multilevel) methods. Solvers belonging to one class may employ methods from another class to speed their convergence, and so Krylov space methods commonly use incomplete factorisation and simple iterative methods as preconditioners, whilst multigrid schemes use incomplete factorisation and simple iterative methods as smoothers. In Chapter 3 these direct and iterative solvers are described and implemented for the solution of equations arising from the finite volume discretisation of the two- and three-dimensional transport equation. The methods are compared in their ability to solve two- and three-dimensional test problems, being rated in terms of the speed of solution and the memory usage, with the scaling of the methods as the number of equations increases being examined. The convergence of a selection of iterative solvers is shown in Figure 1.3 along with the variation of the solution time as the number of equations is increased.
Figure 1.4: The streamfunction for the Re driven cavity (left) and the Ra natural convection (right) benchmark test problems.
, Pr
The motion of a viscous fluid is described by the Navier–Stokes equations. Using the finite volume discretisation techniques developed in Chapter 2 two algorithms for the solution of the Navier–Stokes equations are described in Chapter 4. A primitive variable form of the equations is used, with the flow being directly solved for the velocity and the pressure field , along with any ancillary scalar fields such as temperature . The first solution technique couples the velocity and pressure fields using the
CHAPTER 1. INTRODUCTION
4
Figure 1.5: The calculated streamfunction for the Re
skew driven cavity problem.
the SIMPLE algorithm[123, 122], as modified by Rhie–Chow[139] to enable a steady state solution to be obtained on a collocated or non-staggered mesh. A second coupling scheme, a fractional step method[4, 5], is derived for calculating transient flows, although it may be used in the limit as to calculate steady flow, whilst the steady SIMPLE formulation may used in a Crank–Nicolson or backwards Euler time stepping manner to calculate transient flows.
To test the solver the steady state SIMPLE coupling scheme has been applied to the solution of two benchmark flow problems, a driven cavity flow[50] and natural convection in a square cavity[30], the streamfunctions for the two flows being shown in Figure 1.4. The solutions have been calculated using the different discretisation schemes developed in Chapter 2, and the effect of mesh refinement on solution accuracy is examined. The solution techniques described in Chapters 2 and 4 were developed for Cartesian meshes, where the mesh is aligned with the coordinate axes. The finite volume methods are extended in Chapter 5, following the methods of Peri´c and Jones[72, 129, 43], such that they may be used on grids that are ). The finite neither planar (with surfaces remaining flat) nor Orthogonal (mesh angles being at volume discretisations are rederived and implemented for a general non-orthogonal mesh, and the SIMPLE coupling scheme is extended for use with these non-orthogonal discretisations. The code is tested upon three test problems, with the solutions being compared with published benchmarks (such as the skew driven cavity flow shown in Figure 1.5), and with solutions calculated upon Cartesian grids.
Speed for DEC Alpha 500/266 - HPF (forall)
Speed for Thinking Machines CM-5 - CMF (cshift)
400
2500 1 cpu 2 cpu 4 cpu 6 cpu 8 cpu 10 cpu
350
1 cpu 16 cpu 32 cpu 64 cpu 128 cpu 2000
250
Speed (MFLOPS)
Speed (MFLOPS)
300
200
150
1500
1000
100 500 50
0
0 0
20
40
60
80 Grid Size
100
120
140
160
0
50
100
150
200
250
Grid Size
Figure 1.6: The speed (in MFLOPS) of a three dimensional conjugate gradient solver, plotted as a function of array size. Speeds were calculated running on a cluster of networked DEC Alpha workstations (left), and on a Thinking Machines CM-5 parallel computer (right). The efficiency of the solution methods are studied in Chapter 6, with the calculation speed on parallel and vector supercomputers being compared with that on Intel PC’s and DEC Alpha RISC worksta-
CHAPTER 1. INTRODUCTION
5
Streamfunction
Temperature
Figure 1.7: Streamfunction (top) and isotherms (bottom) for supercritical Rayleigh–B´enard convection, for Ra = 3500. tions. To examine issues in code implimentation a test code solving a three dimensional conduction problem using a conjugate gradient solver was written using C, FORTRAN 77 and Fortran 90. These codes were used to compare the speed of executables generated by compilers for the different languages, and to compare the executable speed of different code structures. This code was ported to vector and parallel supercomputers, and modified to run in parallel on a cluster of networked workstations. Parallelisation was implemented using High Performance Fortran, a data parallel language, and using two message passing libraries, PVM and MPI. The codes were compared in terms of both performance and code complexity, with performance being measured for a range of array sizes to examine the scaling of performance, and to understand the effect of memory hierarchy on program performance. Examples of the speed of the code when run on two parallel platforms, a DEC Alpha workstation cluster and a Thinking Machines CM-5, are shown in Figure 1.7, plotted as a function of mesh size. After studying the performance of the test codes, the developed Navier–Stokes solver and an existing CFD code were parallelised using High Performance Fortran and run on a variety of parallel computers to measure the performance gains that may be attained in a real fluid modelling code.
Figure 1.8: Experimentally observed structures in an intrusion flow viewed from above (left), from Sch¨opf and Stiller[153], and the isotherms for a numerical simulation of the flow (right), viewed from the front of the cavity. After these tests, the codes were used to model two fluids problems. First, a study of transient internal natural convection modelled in two and three dimensions is discussed in Chapter 7. An initial study modelled two dimensional Rayleigh–B´enard convection (shown in Figure 1.7), with the results being
CHAPTER 1. INTRODUCTION
6
Free Surface Elevation at Fr 0.3 Re 10000
Figure 1.9: Free surface for the flow around a Wigley ship’s hull, for Fr
and Re
.
compared with experimental results and the predictions of the analytic theory of Rayleigh[137]. The code was shown to accurately predict the onset of natural convection, although the relationship between wavelength and Rayleigh number differed between the numerical calculation and experiment. The CFD codes were then applied to a transient natural convection flow in a side wall heated cavity, similar to the benchmark natural convection flow shown in Figure 1.4. Experimental observations of the onset of convection[153] revealed three dimensional Rayleigh–B´enard like structures in the thermal intrusion at the onset of convection, and the code was used to examine the nature of these structures. Finally, in Chapter 8 the application of the code to the modelling of three free surface flow problems is discussed. The method used for capturing the free surface geometry was developed and then tested on two simple two-dimensional flow problems, a free surface driven cavity and a two-dimensional approximation of a flume flow. The code was then applied to the modelling of the three dimensional flow about a ship’s hull. A realistic wave wake pattern for flow past a Wigley ship’s hull is then obtained using the code. (shown in Figures 1.9 and 1.10).
Free Surface Elevation at Fr 0.3 Re 100000
Figure 1.10: Contours of free surface elevation for the flow around a Wigley ship’s hull, for Fr and Re .
!
Chapter 2
Finite Volume Differencing Schemes This chapter discusses the basic techniques for the numerical solution of Partial Differential Equations (PDEs) using Finite Volume approximations. With analytic methods the solution to a PDE is found for all locations within the domain of interest. However, with finite volume or finite difference methods the solution is found only for a set of discrete points within the space of the domain. For finite difference methods the domain has a number of points placed within it, and the equations are solved at each point. In contrast, for finite volume methods the domain is divided up into a number of control volumes, with the value at the centre of the control volume being held to be representative for the value over the entire control volume. By integrating the original PDE over the control volume the equation is cast into a form that ensures conservation. The derivatives at the faces of the volume are approximated by finite difference equations, and a system of sparse linear equations is generated that can be solved using a standard linear method. The process of transforming a PDE into a set of linear equations is termed discretisation. In this chapter we will only consider the discretisation of a generic transport equation upon a Cartesian mesh (ie: one aligned with the , and axis of the domains coordinate system). The notation used for the discretisation is described in Section 2.1, with the steady transport equation being discretised in Section 2.2, and the transient transport equation being discretised in Section 2.4. The accuracy of the different schemes is tested in Section 2.5 for both transient and steady flow problems. The methods used to solve the linear systems generated by the finite volume discretisations are given in Chapter 3, whilst the methods used to solve the Navier–Stokes equations will be described in Chapter 4. Further finite volume discretisations that are not restricted to Cartesian meshes and can be used on a general curvilinear mesh are derived in Chapter 5.
2.1
Equations and Notation
The transient transport equation for the conservation of a specific1 scalar in a fluid undergoing advection and diffusion can be written as
(2.1)
where is the diffusion coefficient, any source term for the scalar per unit volume, the density of the fluid and the fluid’s velocity field. The first term in Equation (2.1) is the time derivative measuring the rate of change of with respect to time, and the second term is the advective component representing the transport of by the ambient velocity field. On the right hand side of the equation 1 ie:
one defined per unit mass
7
CHAPTER 2. DIFFERENCING SCHEMES
8
the first term is the transport due to diffusion, whilst the final term is the contribution due to sources of within the field. For a one dimensional domain Equation (2.1) can be written
!" # " $ %
&
(' )
(2.2)
whilst in two and three dimensions it may be written as
& ,' ' )
(2.3)
" $ # " * +- , & & & (' ' ' ) respectively, with $ , * and - being the , and components of velocity.
(2.4)
" " $ +*
&
and
For the steady transport equation the time derivative is set to zero, and so the leading term may be dropped from Equations (2.1) to (2.4). For a finite volume discretisation the domain over which the transport equation is to be solved is divided up into a number of discrete volumes, as in Figure 2.1, with each volume having a representative value located at it’s centre. For a one dimensional domain the cells are numbered from . to / , with the cell at point 0 having neighbours 0#12. and 03 4. , the interfaces with those neighbouring volumes being 0 1675 and 0# 875 . For a two dimensional domain the cells are numbered 09 ;:< , and for a three dimensional domain they are numbered 0= >:+ =? .
@
A
A
A
5 @
H H
A
B C 5
A
@
@
B
@
A
BED 5
A
@
@
@
@ @
F @
F#I @ A @
A
A
A
A
A
A
A
@
@
A
A
A
A
A
A
A
@
A
A
A
A
A
A
@
@
A
A
A
A
A J MH L B"K A
A
A
@
@
A
A
A
A
A
A
A
@
@
A
D 5 H @ C 5
5
@ 5
@
A
A @
A
A @
A
B C 5
A
A
A @ B
@
A
BED 5
A @
A
A @
@
@
@ F G )
Figure 2.1: The discretisation of a 1D (top) and 2D (bottom) domain into Cartesian finite volumes. For convenience, the commonly used compass notation originating from Imperial College[43, 122] will be used when referring to the neighbours of a cell. The subscript N signifies the cell upon which the equation is being discretised, whilst its immediate neighbours in the axis are labelled O and P (for East and West), in the axis / and (for North and South), and in the axis Q and R (for Top and Bottom). A capital letter signifies a value at the centre of the neighbouring cell, whilst a
CHAPTER 2. DIFFERENCING SCHEMES
9
Grid location
Compass notation
09 ;:+ =? 0S1T.U >:+ =? 0V .U >:+ =? 09 ;:1T.U =? 09 ;: .U =? 09 ;:+ 9? 1T. 09 ;:+ 9?W . 0S1 75 >:U 9? 0V 75 >:U 9? 0= >:W1 75 =? 0= >: 75 =? 09 ;:+ 9? 1 75 09 ;:+ 9?W X75
P W E S N B T w e s n b t
Table 2.1: Conversion between the mesh index and the compass notation. lower case letter signifies a value at the interface between the two cells. Cells further away from the N cell are named using the chain of locations that link them back to the central point N . For instance the cell that is to the North of the East neighbour is called /YO , whilst the cell to the West of the West neighbour is PZP . A diagram of these neighbouring nodes which make up the computational molecule at a point is given in Figure 2.2, whilst a table converting between this compass notation and the cell indexing previously referred to is given in Table 2.1.
a9b [ \[
F
[]
[c
ed ab
\[
][
^[
\_[ \
\[
[
^[
\_[ \
][
FF [ \ [ [ ^ F F F [ \ [ ] [ ^ ^`[ ^ [c\ [c [c^ [ cMc
Figure 2.2: The computational molecule at point N illustrating the compass notation used.
^[
`^ [ ^ [ F
[] \[ [ f h f f f [c
[
g ^[ i k ij
ed ab
for one-, two- and three-dimensional meshes,
This chapter will only describe discretisations on Cartesian meshes, with the complications of nonorthogonal meshes being left to Chapter 5. The geometry for a cell in a one-dimensional mesh is given in Figure 2.3 with the geometry for a cell in a two-dimensional Cartesian mesh being given in Figure 2.4. The cell has width l and height lm , whilst the distance between the centre of the cell and that of it’s immediate neighbour on the right (or east) side is n o , and that with it’s left (or west) neighbour is n p . For a one-dimensional mesh the area of the cell faces are taken to be unitary, ie: q o rq p s. , whilst the volume of the cell is lt68l . With the two-dimensional Cartesian mesh the cell face areas are q o uq p vlm and qxwy8q{z|vl , whilst the volume is lt}8l lm . For a threedimensional Cartesian mesh the cell face areas are q o rq p ~lmlm Mq{wrqz~l|@?
(3.32)
This suggests an iterative method defined by,
c Q e >@? >@c ? V f $ $ T @> ?
(3.33)
where the terms on the right hand side of the equation are all from the previous iteration. This algorithm is given in Figure 3.6.
U Q Õi¢D£¤ andc Ä ¥ Ä ¢D£¤ while Ó c c T V V Ó Ä -È Ó È Î %×Ö % Q set
8
Figure 3.6: The Jacobi method.
The Jacobi method requires words of storage over that required for the storage of the linear system and its solution. Note that there is no intrinsic order in which to update each equation, so the method is trivial to parallelise.
CHAPTER 3. LINEAR SOLVERS
61
Successive Over Relaxation (SOR) If the equations in (3.32) are solved in sequence so that the updated values of they are available, then we obtain the Gauss–Seidel algorithm,
are used as soon as
ÛÜ >@c ? >@c J c Q @ > ? @ > ? ? V T @> ? $ $ÙØÚ T @> ? + / 2 (
(3.34)
assuming a case where the equations are arranged in increasing order from west to east, south to north, bottom to top. If the iterative update is under- or over-relaxed then we arrive at the Successive Over Relaxation scheme,
>@c ? c Ö @ > ? >@? T T $ Ý > @ ? > @ ? $ ÚØ + / 2 (
>@c ?J V ÛÜ
c $
T
IQ
K cV T)Ö $ "
(3.35)
depending on the previous The equations now contain a data dependency, with the new value of equations updated values. This limits the ability to parallelise or vectorise the algorithm. The dependency can be removed by changing the order in which the equations are solved (such as is done with RBSOR) but this in turn will affect the rate of convergence, typically for the worse.
U Q Õ¡¢D£¤ and Ä ¥ Ä ¢D£¤ while Q 6Z "]"^" 8 for Ó c Þ c « k T j V c T j U c V m ħ !È Ó ÈnÎ V %×Ö Ó % Q set
Figure 3.7: The SOR algorithm.
Ö
The choice of the value of is crucial to the convergence of the SOR method. A theorem due to Kahan[74] shows that SOR fails to converge if is outside the interval (0,2). If the relaxation term then SOR reduces to the Gauss-Seidel method.
Ö Q
Ö
Q 6Z ^" "]" 8 Q [ "]"^" 8 T Q
Symmetric SOR and Red-Black SOR utilise the same algorithm but with changes made to the order of the operations. With Symmetric SOR every sweep through the equations in the order is followed by a sweep in the reverse order, . With Red-Black SOR the operations are done in two passes, the first pass operating on the odd elements of the array, , followed by the second pass which operates on the even elements, .
ß8& 8 T Q ^" "]" Q
-Z à "^"]" 8
The advantages of Red-Black SOR is that dependencies between adjacent array values are removed, enabling the method to be vectorised or parallelised. However, as is seen in Figures (3.22) to (3.25), this is at the cost of a greatly reduced rate of convergence. The SOR methods require no storage above that of the original system.
3.1.5
Incomplete Factorisation Methods
FaG
A second class of iterative methods is developed by using an approximate decomposition of the array into it’s factors. Instead of doing a complete factorisation only a restricted number of offdiagonal elements of the and matrices are allowed to take on a non-zero value, typically only
F
G
CHAPTER 3. LINEAR SOLVERS
62
the diagonals that are non-zero in the matrix and one or two further bands. For each iteration a forward and backward substitution process using this incomplete factorisation is applied to the residual formulation of the system of equations. These methods are efficient in their own right, but also have value as preconditioners for the Krylov space solvers, and as smoothers for the multigrid schemes. The most commonly encountered implementations are Incomplete Cholesky factorisation (commonly used as a preconditioner for the Conjugate Gradient method) and the Strongly Implicit method of Stone[163].
Incomplete Cholesky Factorisation (IC)
á Å
FHF
The simplest incomplete factorisation scheme is the Incomplete Cholesky method, where a symmetric matrix is factored into a system, with infill only in the locations where .
U Q Õ¡¢D£¤ c Ä ¥ Ä ¢D£¤ whilec Ç V È c Èoand Ä Pâ Ç V CZ T Î V Q "^"] 6"
