Download Computational Hydraulics Process, Lecture Notes- Physics - 1 and more Study notes Physics in PDF only on Docsity!
9. THE CFD PROCESS SPRING 2012
9.1 Introduction 9.2 The computational mesh 9.3 Boundary conditions 9.4 Flow visualisation Examples
9.1 Introduction
9.1.1 Stages of a CFD Analysis
A complete CFD analysis consists of:
- pre-processing ;
- solving ;
- post-processing.
This course has focused on the “solving” process, but this is of little use without pre- processing and post-processing facilities. Commercial CFD vendors supplement their flow solvers with grid-generation and flow-visualisation tools, as well as graphical user interfaces (GUIs) to simplify the setting-up of a CFD analysis.
Pre-Processing The pre-processing stage consists of:
- determining the equations to be solved;
- specifying the boundary conditions ;
- generating a computational mesh or grid.
It depends upon:
- the desired outputs of the simulation (e.g. force coefficients, heat transfer, ...);
- the capabilities of the solver.
Solving In commercial CFD packages the solver is often operated as a “black box”. Nevertheless, user intervention is necessary – to set under-relaxation factors and input parameters, for example – whilst an understanding of discretisation methods and internal data structures is required in order to supply mesh data in an appropriate form and to analyse output.
Post-Processing The raw output of the solver is a huge set of numbers corresponding to the values of each field variable ( u , v , w , p , …) at each point of the mesh. This must be reduced to some meaningful subset and/or manipulated further to obtain the desired predictive quantities. For example, a subset of surface pressures and cell-face areas is required to compute a drag coefficient.
Commercial packages routinely provide:
- plotting tools to visualise the flow;
- analysis tools to extract and manipulate data.
9.1.2 Commercial CFD
The table below lists some of the more popular commercial CFD packages.
Developer/distributor Code(s) Web address (liable to change!) CD adapco STAR CD STAR CCM+
http://www.cd-adapco.com/
Ansys FLUENT CFX
http://www.ansys.com/Products/
Flow Science FLOW3D http://www.flow3d.com/ EXA PowerFLOW http://www.exa.com/ CHAM PHOENIX http://www.cham.co.uk/
A popular non-commercial CFD solver is OpenFOAM (http://www.openfoam.com/)
An excellent web portal for all things CFD is http://www.cfd-online.com/.
9.2 The Computational Mesh
9.2.1 Mesh Structure
The purpose of the mesh generator is to decompose the flow domain into control volumes.
The primary outputs are:
- cell vertices ;
- connectivity information.
Precisely where the nodes are relative to the vertices depends on whether the solver uses, for example, cell-centred or cell-vertex storage. Further complexity is introduced if a staggered velocity grid is employed.
The shapes of control volumes depend on the capabilities of the solver. Structured-grid codes use quadrilaterals in 2-d and hexahedra in 3-d flows. Unstructured-grid solvers often use triangles (2-d) or tetrahedra (3-d), but newer codes can use arbitrary polyhedra.
tetrahedron hexahedron
cell-centred storage (^) cell-vertex storage
u p
v
staggered velocity mesh
Volumes
If r ≡( x , y , z )is the position vector, then = 3 ∂
z
z y
y x
x r. Hence, integrating over
an arbitrary control volume and using the divergence theorem gives for the volume of a cell:
⌡
= ⌠^ •
∂ V
V r d A 3
If the cell has plane faces this can be evaluated as
faces
V r f A f 3
where r f is any convenient position vector on a face and A f is the face vector area, since, for any other vector r in that face,
142 4 434 0
=
r • A f = r f • A f + r − r f • A f
The last term vanishes because r – r f is perpendicular to A f for any point on a plane face.
If the cell faces are not planar, then the volume depends on how these faces are broken down into triangles. Typically the area of a particular face can be regarded as an assemblage of triangular elements connecting the vertices of that face with a central reference point formed, e.g.
vertices
r (^) f (^) N r i
Examples for the commonest shapes follow.
Tetrahedra The volume of a tetrahedron formed from side vectors s 1 , s 2 , s 3 (taken in a right-handed sense) is
6 1 2 3 V =^1 s • s ∧ s
Hexahedra Volumes of arbitrary polyhedral cells are:
faces
V (^) 3 r f A f
For hexahedra, on each face the reference point and vector area are: r (^) f = 41 ( r 1 + r 2 + r 3 + r 4 ) A (^) f = 21 ( r 1 − r 3 )∧( r 2 − r 4 )
To obtain an outward -directed face vector, the points with position vectors r 1 , r 2 , r 3 , r 4 should be in clockwise order when viewed along the outward normal (right-hand screw).
This is a generalisation to arbitrary hexahedra of the result for Cartesian control volumes.
r f
s 3 s 1
s 2
2-d Cases
In 2-d cases one can formally consider cells to be of unit depth.
The “volume” of the cell is then numerically equal to its planar area.
The side “area” vectors are most easily obtained from their Cartesian projections :
d
d d x
y A
To obtain an outward -directed d A , the edge increments (d x ,d y ,0) must be taken anticlockwise around the cell.
Volume-Averaged Derivatives
By applying the divergence theorem to φ e x , where e x is the unit vector in the x direction, the volume-averaged derivative of a scalar field φ is
⌡
= ⌠^ φ ⌡
∂φ
∂
∂φ ∂ V
x V
A
V
V
x V x
d
d
or
faces V f A f
Thus, for arbitrarily-shaped cells, average derivatives may be derived from the values of φ on the cell faces, together with the components of the face area vectors and the cell volume.
e.g. for hexahedra:
V w w e e s s n n b b t^ t z
y
x = φ A +φ A +φ A +φ A +φ A +φ A ∂φ ∂
∂φ ∂
∂φ ∂
The values on cell faces, φ w , φ e etc. can be obtained by interpolation from the nodes on either side.
φe
φn
φw
φs A^ e
A n
A w
A s
s
d A
d dy
dx
V
9.2.4 Fitting Complex Boundaries With Structured Grids
Blocking Out Cells
The range of flows which can be computed in a rectangular domain is rather limited. Nevertheless, a number of significant bluff- body flows can be computed using a single- block Cartesian mesh by the process of blocking out cells. Solid-surface boundary conditions are applied to cell faces abutting the blocked-out region, whilst values of velocity and other flow variables are forced to zero by a modification of the source term for those cells. In the notation of Section 4, where the scalar-transport equations for a single cell are discretised as
a P φ P −∑ aF φ F = bP + sP φ P
the source terms are simply re-set: bP → 0 , sP →−
where is a large number (e.g. 10^30 ). Rearranging for φ P , this ensures that the computational variable φ P is effectively forced to zero in these cells. However, the computer still stores and carries out operations for these points, so that it is essentially performing a lot of redundant work. An alternative approach is to fit several structured grid blocks around the bluff body. Multi-block grids will be discussed further below.
Volume-of-Fluid Approach
The numerical simplicity and solver efficiency associated with Cartesian grids mean that some practitioners still attempt to retain this grid geometry even for complex curved boundaries: both for solid walls and free surfaces. In this volume-of-fluid approach the fraction f of the cell filled with fluid is stored: 0 outside the fluid domain, 1 within the interior of the fluid and 0 < f < 1 for cells which are cut by a boundary. At solid boundaries f is determined by the surface contour. At free surfaces, f emerges as part of the solution procedure.
mesh for 2-d rib with blocked-out cells
f=
f=
0 Body-Fitted Grids
The majority of general-purpose codes employ body- fitted (curvilinear) grids. The mesh lines/coordinate surfaces are distorted so as to fit snugly around curved boundaries. Accuracy in turbulent-flow calculations demands a high density of grid cells close to solid surfaces and the use of body-fitted meshes means that the grid need only be refined in the direction normal to the surface, with consequent saving of computer resources.
However, the use of body-fitted grids has important consequences.
- It is necessary to store detailed geometric components for each control volume; for example, in our research code STREAM we need to store ( x , y , z ) components of the cell-face-area vector for “east”, “north” and “top” faces of each cell, plus the volume of the cell itself – a total of 10 arrays.
- Unless the mesh happens to be orthogonal, the diffusive flux through the east face (say) is no longer given exactly by
A A n (^) PE
= φ E^ −φ P ∂
∂φ
because the discretised derivative of φ normal to the face involves cross-derivative terms parallel to the cell face and nodes other than P and E. The extra off-diagonal diffusion terms are typically transferred to the source term.
- A similar alignment problem occurs with the advection terms because, in general, interpolated values of all three velocity components are needed to evaluate the mass flux through a single face. This necessitates approximations in the pressure- correction equation that can slow down the solution algorithm.
Multi-Block Structured Grids
In multi-block structured grids the domain is decomposed into a small number of regions, in each of which the mesh is structured (i.e. cells can be indexed by ( i , j , k )).
A common arrangement (and that assumed by our own code STREAM), is that grid lines match at the interface between two blocks, so that there are cell vertices that are common to two blocks. However, some solvers do allow overlapping blocks ( chimera grids ) or block boundaries where cell vertices do not align. Interpolation is then needed at the boundaries of blocks.
2 3 4
1 5 2-d rib with multi-block mesh
P
E^ η^ =const.
const. ξ (^) =
u
v
vectors; e.g.
( , , ) A (^ ^ )
x y z V
where V is volume and )
(^
A is the area vector normal to face = constant (in the direction of increasing ). Curvilinear coordinates can be taken as the node indices, in which case = 1.
Non-orthogonal meshes are more computationally-intensive, but their additional geometric flexibility makes them desirable in general-purpose solvers.
9.2.6 Disposition of Grid Cells
Only two points are necessary to resolve a straight-line profile, but the more curved the profile the more points are needed to resolve it accurately. In general, more points are needed in rapidly-changing regions of the flow, such as:
- solid boundaries;
- separation, reattachment and impingement points;
- flow discontinuities; e.g. shocks, hydraulic jumps.
Simulations must demonstrate grid-independence , i.e. that a finer-resolution grid would not significantly modify the solution. This generally requires a sequence of calculations on successively finer grids.
Boundary conditions for turbulence modelling do impose some limitations on the cell size near walls. Low-Reynolds-number models (resolving the near-wall viscous sublayer) generally require that y +^ < 1 for the near-wall node, whilst high-Reynolds-number models relying on wall functions require (in principle) the near-wall node to lie in the log-law region, ideally 30 < y +^ < 150.
9.2.7 Multiple Grids
Multiple grids – combining cells so that there are 1, 1/2, 1/4, 1/8, ... times the number of cells in some basic fine grid – are used deliberately in so-called multigrid methods. These calculate the solution on alternately coarser and finer grids. The idea is that the solution is obtained quickly on the coarsest grid where the number of cells is small and changes propagate rapidly across the domain. The solution is then refined locally on the finest grid to obtain the most accurate solution.
If two levels of grid are used then a process known as Richardson extrapolation may be used to both estimate the error and, possibly, refine the solution. If the basic advection-diffusion discretisation is known to be of order n and the exact (but unknown) solution of some property φ is denoted φ, then one would expect the error φ – φ^ to be proportional to n ,
where is the mesh spacing. Thus, for solutions φ^ and φ 2 ^ respectively on two grids with
mesh spacing and 2 ,
n
n
C
C
2 * (^2 )
φ =φ +
φ =φ +
Two equations give us two equations for two unknowns, φ*^ and C , which we can solve to get a better estimate of the solution:
*^2 2
φ − φ φ = n
n
and the error in the fine-grid solution:
2 ^
φ − φ C n = (^) n
Unstructured-grid methods offer the possibility of local grid refinement : that is, adding more cells in regions where the error is estimated to be high.
9.3 Boundary Conditions
There are a number of common boundary types.
INLET
- Inflow transported variables specified on the boundary, either by a predefined profile or (my preference) by doing an initial 1-d, fully-developed-flow calculation.
- Stagnation (or reservoir ) total pressure and total temperature (in compressible flow) or total head (in incompressible flow) fixed; usual inflow condition for compressible flow.
OUTLET
- Outflow zero normal gradient (∂φ/∂ n = 0) for all variables.
- Pressure as for outflow, except fixed value of pressure; usual outlet condition in compressible flow if the exit is subsonic.
- Radiation (or convection ) prevent wave-like motions from reflecting at outflow boundaries by solving a simplified first-order wave equation with outward-directed wave velocity.
WALL
- Non-slip wall the default case for solid boundaries (zero velocity relative to the wall; wall stress computed by viscous-stress or wall-function expressions).
- Slip wall only the velocity component normal to the wall vanishes; used if it is not necessary to resolve a thin boundary layer on an unimportant wall boundary.
OTHER
- symmetry plane ∂φ/∂ n = 0, except for the velocity component normal to the boundary, which is zero; used where there is a geometric plane of symmetry, but also as a far-field boundary condition, because it ensures that there is no flow through, nor viscous stresses on, the boundary.
- periodic used in repeating flow; e.g. rotating machinery, regular arrays.
Line Contour Plots
A contour line (“isoline”) is a line along which some property is constant. The equivalent in 3 dimensions is an “isosurface”. Any field variable may be contoured. In contrast to line graphs, contour plots give a global view of the flow field, but are less useful for reading off precise numerical values. If the domain is linearly scaled then detail occurring in small regions is often obscured.
The actual numerical values of the isolines are sometimes less important than their overall disposition. If contour intervals are the same then clustering of lines indicates rapid changes in flow quantities. This is particularly useful in locating shocks and discontinuities.
Shaded Contour Plots
WARNING: colour plots are expensive to print or photocopy! This proviso apart, colour is an excellent medium for conveying information and good for on-screen and presentational analysis of data. Simple packages flood the region between isolines with a fixed colour for that interval. The most advanced packages allow a pixel-by-pixel gradation of colour between values specified at the cell vertices, together with lighting and other special effects such as translucency.
Grey-scale shading is an option if plots are to be reproduced in black and white.
Vector Plots
Vector plots display vector quantities (usually velocity; occasionally stress) with arrows whose orientation indicates direction and whose size (and sometimes colour) indicates magnitude. They are a popular and informative means of illustrating the flow field in two dimensions, although if grid densities are high then either interpolation to a uniform grid or a reduced set of output positions are necessary to prevent the number density of arrows blackening the plot. There can sometimes be
problems when selecting scales for the arrows when large velocity differences are present, especially in important areas of recirculation where the mean flow speed is low. In three dimensions, vector plots can be deceptive because of the angle from which they are viewed.
Streamline Plots
Streamlines are parallel to the mean velocity vector. They can generally be obtained by integration:
u
x
d t
d
(the only option in three dimensions).
In 2-d incompressible flow a more accurate method is to contour the streamfunction. This function is defined (and calculated) by fixing arbitrarily at one point and then setting the change in between two points as the volume flow rate (per unit depth) across any curve joining them:
⌡
− =⌠^ •
2
1
2 1 u n d s
(The sense used here is clockwise about the start point, although the opposite sign convention is equally valid). is well-defined because the flow rate across any path connecting two points must be the same or mass would accumulate between them. Contouring produces streamlines because a curve of constant is one across which there is no flow.
For a short path making orthogonal increments d x and d y then d = u d y − v d x and hence the velocity components are related to by
y
u ∂
x
v ∂
Computationally, is stored at cell vertices. This is convenient because the flow rates across the cell edges are already stored as part of the calculation. In the Cartesian cell shown: 2 = 1 − vs^ y 3 = 2 + ue^ x
If isolines are at equal increments in , then clustering of lines corresponds to high velocities, whilst regions where the streamlines are further apart signify low velocities. However, as with vector plots, this has the effect of making it difficult to visualise the flow pattern in low- velocity regions such as recirculation zones and a smaller increment in is needed here.
Particle paths can also be traced along solid surfaces using the wall stress vector. This often reveals important features connected with separation/reattachment/impingement on 3-d surfaces, which defines the basic flow topology.
vs
ue
ψ 3
ψ 1 ψ 2
1
2
dy
dx
v
u
Examples
Q1. (a) Explain what is meant, in the context of computational meshes, by: (i) structured; (ii) multi-block structured; (iii) unstructured; (iv) chimera.
(b) Define the following terms when applied to structured meshes: (i) Cartesian (ii) curvilinear (iii) orthogonal
Q2.
Suppose that your CFD solver requires multi-block, structured, curvilinear grids, where grid nodes must match at block interfaces. Sketch suitable grid topologies that might be used to compute the following flow configurations. Use symmetry where appropriate and “body- fitting” meshes, not “blocked-out” regions.
(a) Harbour configuration (treat as 2-d)
(b) Circular to square-section pipe expansion
(c) Outfall (treat as 2-d)
Q3. (Computational Hydraulics Exam, May 2011 – part) (a) For wind-loading calculations a CFD calculation of airflow is to be performed about the building complex shown in the figure. Sketch a suitable computational domain indicating the specific types of boundary condition that are applied at each boundary of the fluid domain. For each boundary type summarise the mathematical conditions that are applied to each flow variable (velocity, pressure, turbulent scalar).
(b) Define the drag coefficient for a building and explain how it would be calculated from the raw data obtained in a CFD simulation.
Q4.
(a) Two adjacent cells in a 2-dimensional Cartesian mesh are shown below, along with the cell dimensions and some of the velocity components (in m s–1) normal to cell faces. The value of the stream function at the bottom left corner is (^) A = 0. Find the value of the stream function at the other vertices B – F. (You may use either sign convention for the stream function).
A
D F
C
2 1
3
5
0.3 m 0.2 m
0.1 m
E
B
3
(b) Sketch the pattern of streamlines across the two cells in part (a).
Q5. (Computational Hydraulics Exam, June 2009 – part) The figure below shows a quadrilateral cell, together with the coordinates of its vertices and the velocity components on each face. If the value of the stream function at the bottom left corner is = 0, find: (a) the values of at the other vertices; (b) the unknown velocity component vn.
Q6. (Computational Hydraulics Exam, May 2010 – part) A numerical scheme known to be second-order accurate is used to calculate a steady-state solution on two regular Cartesian meshes A and B where the finer mesh A has half the grid spacing of mesh B. The values of the solution φ at a particular point are found to be 0. using mesh A and 0.78 using mesh B. Use Richardson extrapolation to: (a) estimate an improved value of the solution at this point; (b) estimate the error at this point using the mesh-A solution.
velocity face u v w 3 – e 3 1 s 0 2 n 1 vn
(4,5)
(1,1) (3,1)
(1,3)
s
e
n
w
Q11. (MSc Exam, May 2008) In a 2-dimensional unstructured mesh, one cell has the form of a pentagon. The coordinates of the vertices are as shown in the figure, whilst the average values of a scalar φ on edges 1 – 5 are: φ 1 = –7, φ 2 = 8, φ 3 = –2, φ 4 = 5, φ 5 = 0
Find: (a) the area of the pentagon;
(b) the cell-averaged derivatives ∂ x
∂φ and ∂ y
∂φ .
Q12. (*** Advanced ***) (MSc Exam, June 2009) (a) Give a physical definition of the stream function in a 2-d incompressible flow, and show how it is related to the velocity components.
(b) Show that if the flow is irrotational then satisfies Laplace’s equation. If the flow is not irrotational how is ∇^2 related to the vorticity?
The figure below and the accompanying table show a 2-d quadrilateral cell, with coordinates at the vertices and velocity components on the cell edges. The flow is incompressible. All quantities are non-dimensional.
(c) Assuming that the stream function at the bottom left (SW) corner is 0, calculate the stream function at the other vertices.
(d) Find the y -velocity component vn , whose value is not given in the table.
(e) By calculating (^) ⌡
⌠ (^) u • d s for the quadrilateral cell estimate the cell vorticity.
velocity face u v w 3 2 e 3 – s 0 2 n 1 vn
e
n
w s
x
y
(3,1)
Q13. (MSc Exam, June 2010 – part) (a) The vertices of a particular tetrahedral cell in a finite-volume mesh are (0,0,0), (4,0,0), (1,4,0), (1,2,4) Find: (i) the outward vector area of each face; (ii) the volume of the cell.
(b) For the cell in part (a) a scalar φ has value 6 on the largest face, 2 on the smallest face
and 3 on the other two faces. Find the cell-averaged derivatives ∂ x
∂φ , ∂ y
∂φ , ∂ z
∂φ .
(c) For the cell in part (a) determine whether the point (1,2,1) lies inside, outside or on the boundary of the cell. (You must justify your answer.)
Q14. (*** Advanced ***) (MSc Exam, May 2011 – part) A quadrilateral cell in a 2-d finite-volume mesh is shown in the figure below. The coordinates of the vertices are marked in the figure and the velocities on the edges are given in the adjacent table.
(a) Find the area A of the quadrilateral.
(b) Find the cell-averaged derivatives x
u ∂
y
u ∂
x
v ∂
y
v ∂
(c) By calculating the line integral from the velocity and geometric data, confirm that both the cell-averaged velocity derivatives and the discrete form of Stokes’ law,
⌡
= ⌠^ •
∂ A z (^) A u^ d s
give the same value for the cell-averaged vorticity, (^) z. (Note: this is, in fact, true in general)
(0,-1)
(2.5,2)
(-1,3)
(-2,0)
e
n
w
s
x
y
edge u v e 10 12 n 13 6 w 9 - s 1 0
x y
z
O
A(4,0,0)
B(1,4,0)
C(1,2,4)
