- Email: [email protected]
Software Applications
CHAPTER 12
Software Applications Contents 12.1 12.2 12.3 12.4 12.5
Programs for Stability Analysis Structured 1-D Grid Generator Structured 2-D Grid Generators Structured to Unstructured Grid Converter Quasi 1-D Euler Solver 12.5.1 Example run 12.6 Structured 2-D Euler/Navier-Stokes Solver 12.6.1 Example run 12.7 Unstructured 2-D Euler/Navier-Stokes Solver 12.7.1 Example run 12.8 Parallelization References
397 397 398 399 399 400 400 403 405 405 406 407
The software package contains the source codes and executables of several 1-D and 2-D flow solvers and grid generators. Provided are also input datasets and grids for various 1-D and 2-D examples of flow cases. Furthermore, there are two programs for the von Neumann stability analysis of explicit and implicit time-stepping schemes. Finally, you will find examples for code parallelization using different approaches. The aim of the software is to demonstrate how to translate the theoretical principles of the computational fluid dynamics, which were presented in the previous chapters, into a computer code. The programs should be conceived as a basis for further experimentation and enhancements. The source codes are provided under the terms of the GNU General Public License, which means that they may be freely copied and redistributed free of charge. See the file LICENSE.txt in one of the directories for more details. Source codes of the flow solvers and grid generators are written in Fortran-77 and in Fortran-90. The source of the 2-D unstructured flow solver is also provided in modern C++. The programs do not contain any references to external libraries. The source codes are kept as simple as possible, but still flexible enough. No extensive attempts were made to optimize the software for speed or for memory consumption. With a few exceptions, the source codes are organized in such a way that there is just one subroutine or function per file. Sources of the 2-D flow solvers are fully documented using the Doxygen1 tool. All grid, solution, and convergence files are stored in a plain ASCII format. The convergence and the solution files are written out in a form suitable for visualization 1
http://www.doxygen.org.
Computational Fluid Dynamics: Principles and Applications
Copyright © 2015 Elsevier Ltd. All rights reserved.
395
396
Computational Fluid Dynamics: Principles and Applications
with the program Vis2D.2 The format of the convergence and solution files is always the same—there is one column for each variable. The names of the variables are given in the header. In the case of the solution files generated by the unstructured flow solver, the indexes of the nodes of each element are stored after the coordinates and the flow quantities. The first column represents node 1, the second one node 2, etc. The file format is described in detail in the documentation of Vis2D. Because of its simple form, the format can be adapted easily for other visualization tools if necessary. The software package is organized into three highest-level directories: C++, Fortran, and Parallelization. These contain the following subdirectories: • C++ • unstructured_2d—solution of 2-D Euler/Navier-Stokes equations on unstructured triangular grids. • Fortran • analysis—von Neumann stability analysis of 1-D model equations • grid_1d—1-D grid generation (Laval nozzle) • grid_struct_2d—2-D structured grid generation for external and internal flows • grid_unstr_2d—conversion of 2-D structured into unstructured triangular grids • structured_1d—solution of quasi 1-D Euler equations (flow trough Laval nozzle); demonstration of the multigrid method • structured_2d—solution of 2-D Euler/Navier-Stokes equations on structured grids • unstructured_2d—solution of 2-D Euler/Navier-Stokes equations on unstructured triangular grids. • Parallelization • agents—Asynchronous Agents Library • cuda—CUDA • mpi—MPI • openmp—OpenMP • serial—serial version. The contents of the above subdirectories are further explained in the sections below. In general, the subdirectory of each simulation code contains a README.txt file. It describes the particular files, how to compile and to run the application, and how the results are stored. Additionally, the README.txt file explains the meaning of the main variables and of the input parameters. In the case of the 2-D flow solvers, the subdirectory doc contains the documentation in HTML format (index.html). The subdirectory of each flow solver also has the folder run, where you can find input files, grids, solutions, and graphics for the various demonstration cases. 2
Available from CFD Consulting & Analysis (http://www.cfd-ca.de).
Software Applications
The last remark concerns the convergence history. The convergence of all flow solvers provided here is measured as the 2-norm of the difference of the density variable from two consecutive time steps, i.e., N ρ2 = (ρIn+1 − ρIn )2 . (12.1) I=1
The summation is carried out over all N control volumes. It is convenient to normalize the convergence measure in Eq. (12.1) with its value from the first iteration. The programs store the normalized convergence history directly in logarithmic scale.
12.1 PROGRAMS FOR STABILITY ANALYSIS The directory analysis contains two programs for the von Neumann stability analysis (Section 10.3) of linear 1-D model equations. The first program computes the Fourier symbol and the magnitude of the amplification factor for the explicit multistage time-stepping scheme (Section 6.1.1) and the hybrid scheme (Section 6.1.2). The source code is provided in the subdirectory mstage. The second program analyses the damping properties of a general implicit scheme (Section 6.2). It is contained in the subdirectory implicit. Both programs can deal with either the 1-D convection or the 1-D convection-diffusion model equation (Sections 10.3.2 and 10.3.3). The spatial discretization can be conducted either by the central scheme with artificial dissipation, by the 1st-order upwind, or by the 2nd-order upwind scheme, respectively. The same schemes can be also applied to the implicit operator. In the case of the explicit scheme, the effect of the (central or upwind) implicit residual smoothing (Section 9.3) can be investigated as well.
12.2 STRUCTURED 1-D GRID GENERATOR The directory grid_1d contains a program for the generation of 1-D structured grids. Each grid node is associated with a certain area. The area distribution along the x-axis corresponds to the Laval nozzle. The area of the nozzle is evaluated using the relation 1 πx A(x) = 1 + (A1 − 1) 1 + cos for 0 ≤ x ≤ xthr 2 xthr (12.2)
1 π(x − xthr ) A(x) = 1 + (A2 − 1) 1 − cos for xthr < x ≤ 1, 2 1 − xthr where xthr denotes the location of the throat. The length of the nozzle is supposed to be one. Furthermore, in Eq. (12.2) A1 represents the inlet area and A2 the outlet area
397
398
Computational Fluid Dynamics: Principles and Applications
A2
A1
Ai+1/2
Flow 1
2
i
ib2
imax
Control volume
Figure 12.1 Grid and control volume for the 1-D Euler solver; points i = 1 and i = imax are dummy points.
(see the sketch in Fig. 12.1). The area of the throat results from Eq. (12.2) to A∗ (x = xthr ) = 1. The grid generated by this program serves as an input to the quasi 1-D Euler solver from Section 12.5.
12.3 STRUCTURED 2-D GRID GENERATORS The directory grid_struct_2d contains three different programs for the generation of 2-D structured grids for external and internal flows. The first program, which source code is provided in the subdirectory cgrid, generates a C-type grid (see Section 11.1.1) around an airfoil. An example can be seen in Fig. 11.6. The initial grid is generated algebraically by using the linear TFI method Eq. (11.5). Afterward, elliptic PDEs (Section 11.1.3) are employed to produce boundary-orthogonal grid with specified wall spacing. The airfoil contour is approximated by a Bézier spline (see Refs. [1, 2] for an introduction to Bézier splines). The second program is intended for the generation of an algebraic grid inside a 2-D channel with circular bump. The grid points are clustered around the front and the rear point of the bump. The source code can be found in the subdirectory channel. The next program, which is located in the subdirectory hgrid, generates an H-type grid (Section 11.1.1) for a cascade. An example of grid created by the program is presented in Fig. 11.8. The linear TFI technique (Eq. (11.5)) for algebraic grid generation is employed in this case. The contour of the blade is again described by a Bézier spline, which interpolates the given points. The programs cgrid and hgrid rely on the library provided in the subdirectory srccom. This library contains routines for spline interpolation, linear TFI, elliptic grid generation, and for grid stretching.
Software Applications
12.4 STRUCTURED TO UNSTRUCTURED GRID CONVERTER The directory grid_unstr_2d contains a program for the conversion of 2-D structured grids into unstructured triangular grids. The program divides each quadrilateral cell of the structured grid into two triangles by connecting diagonal nodes with the shortest distance. The triangulated grids can be used as an input for the unstructured flow solver in Section 12.7.
12.5 QUASI 1-D EULER SOLVER The directory structured_1d contains a program for the solution of quasi 1-D Euler equations. The equations govern the inviscid flow in a nozzle or in a tube. They can be written in conservative, differential form as [3] ∂W ∂ Fc + = Q. (12.3) ∂t ∂x The vectors of conservative variables, convective fluxes, and the source term read ⎡ ⎡ ⎤ ⎤ ⎡ ⎤ ρuA ρA 0 = ⎣ ρuA ⎦, Fc = ⎣(ρu2 + p)A⎦, Q = ⎣p dA/dx⎦, W (12.4) ρEA ρHuA 0 where A denotes the nozzle area. The total enthalpy H is given by the formula (2.12), and the pressure results according to Eq. (2.29) from u2 p = (γ − 1) ρ E − . (12.5) 2 The governing equations (12.3) are discretized on a structured grid using the dual control-volume methodology (see Section 4.2.3). A sketch of the grid and of the control volume is displayed in Fig. 12.1. Three different schemes are employed for the spatial discretization: • Central scheme with scalar artificial dissipation (Section 4.3.1) • CUSP flux-vector splitting scheme (Section 4.3.2) • Roe’s flux-difference splitting scheme (Section 4.3.3). The discretized equations (12.3) are advanced in time using the explicit multistage scheme (Section 6.1.1 or 6.1.2). The program uses local time-stepping (Section 9.1) and the central implicit residual smoothing technique (Section 9.3.1) for convergence acceleration. The source code can be found in the subdirectory src. The input files and the results are provided in the subdirectory run. A second version of the flow solver employs the multigrid method from Section 9.4 for convergence acceleration. The spatial and temporal discretization schemes are otherwise the same as above. The multigrid scheme uses the FAS formulation for the
399
400
Computational Fluid Dynamics: Principles and Applications
solution on coarse grids and FMG to initialize the flow on the finest grid. The source code is contained in the subdirectory srcmg, exemplary input files can be found in the subdirectory runmg. The boundary conditions at the inlet and the outlet plane are implemented in characteristic variables as presented in Section 8.4. The concept of dummy points, described in Section 8.1, is utilized for this purpose (see Fig. 12.1).
12.5.1 Example run As a test example, let us consider flow through a Laval nozzle with the inlet area A1 = 1.5, outlet area A2 = 2.5, and with the following boundary conditions: • Inlet total pressure pt,in = 1.0 × 105 Pa • Inlet total temperature Tt,in = 288.0 K • Outlet static pressure pout = 7.0 × 104 Pa The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and multigrid with five grid levels: σ = 4.5, = 0.8, limiter coefficient K = 1.5, entropy correction coefficient δ = 0.05 · c. The grid for the example can be generated using the program described in Section 12.2. It is also provided as laval.grd in the subdirectory runmg. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named input_r. Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to runmg and typing: % eul1dmg < input_r
This will produce the convergence history displayed in Fig. 12.2. Figure 12.3 shows the resulting Mach-number distribution along the nozzle length.
12.6 STRUCTURED 2-D EULER/NAVIER-STOKES SOLVER The directory structured_2d contains a program for the solution of 2-D Euler and Navier-Stokes equations on structured body-fitted grids. The spatial discretization is based on the cell-centered finite-volume approach described in Section 4.2.1. It employs the central discretization scheme with scalar artificial dissipation (Section 4.3.1), as well as Roe’s upwind scheme presented in Section 4.3.3. The governing equations are integrated in time using an explicit multistage scheme (Section 6.1.1 or 6.1.2), accelerated by local time-stepping (Section 9.1) and the central implicit residual smoothing (Section 9.3.1). Gradients of the velocity components and of the temperature are computed using the dual control-volume approach described in Section 4.4.1.
Software Applications
340
1
0
320
–1
300
–2
280
–3
260
–4
240
Mass flow (kg/s)
log (residual)
Convergence history Mass flow
220
–5 0
50
100
150
200
250
300
Iteration
Figure 12.2 Convergence history of the density residual and of the mass flow.
2.00
1.75
Mach number
1.50
1.25
1.00
0.75
0.50
0.25 0.0
0.1
0.2
0.3
0.4
0.5 x (m)
Figure 12.3 Mach-number distribution over nozzle length.
0.6
0.7
0.8
0.9
1.0
401
402
Computational Fluid Dynamics: Principles and Applications
The program also offers the possibility to simulate incompressible flows at very low Mach numbers (see Section 9.5). In order to keep the code easy to understand, no turbulence model is implemented. Source code of the flow solver is provided in the subdirectory src. The input files, grids, and the results can be found in the subdirectory run. The subdirectory doc contains the documentation (files, modules, functions, variables) in HTML format. The program utilizes the concept of dummy cells (Section 8.1) for the treatment of the boundary conditions. Two layers of dummy cells are employed. A sketch of the grid in the computational domain is displayed in Fig. 12.4 (cf. Fig. 8.1). The program can deal only with single-block grids. However, it is very flexible with respect to the specification and the type of the boundary conditions. The following eight boundary types are provided: • Coordinate cut (for C- or O-type grids) • Far-field (optionally with vortex correction; sub- and supersonic) • Inflow (subsonic only) • Outflow (sub- and supersonic) • Fluid injection (or extraction) • Line periodic • Solid wall (viscous or inviscid) • Symmetry (along x = const. or y = const. line). The implementation of the above boundary conditions closely follows the discussion in Chapter 8. The four boundaries in the computational space can be divided into an arbitrary number of segments. Each of the segments can be associated with a different boundary J 3 JM JL J2
2
4
2 1 0 0
1
2
1
I2
IL
IM
I
Figure 12.4 Grid in the computational space for the structured 2-D Euler/Navier-Stokes solver. There are two layers of dummy cells at each boundary (dashed line). Numbers in circles denote the sides of the computational domain.
Software Applications
i + 1, j + 1
i, j+1
I,J
SI(I,J)
I+1, J
VOL(I,J) i, j
i +1, j SJ(I,J)
Figure 12.5 Control volume and face vectors for the structured 2-D Euler/ Navier-Stokes solver. Indexes (I, J) denote the cell, indexes (i, j) represent the grid node.
condition. This approach is quite similar to the description of block interfaces presented in Section 8.9. In fact, the program could be relatively easily extended to multiblock grids. Definitions of the segments are stored separately from the grid in the topology file (with the file extension top). The definition of the face vectors SI and SJ (i.e., n · S) employed in the solver can be seen in Fig. 12.5. The face vectors are associated with the left and the bottom face of the control volume VOL(I,J) (corresponds to I,J ). They point outwards of the control volume. The face vectors of the remaining two sides of the control volume are obtained as -SI(I+1,J) and -SJ(I,J+1). In this way, we have to store only two face vectors for each control volume.
12.6.1 Example run As a test example, let us compute transonic flow past the symmetric NACA 0012 airfoil with the following boundary conditions along the far-field: • Mach number M∞ = 0.8 • Angle of attack α = 1.25◦ • Static pressure p∞ = 1.0 × 105 Pa • Static temperature T∞ = 288.0 K The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and explicit time-stepping: σ = 5.0, = 1.2, limiter coefficient K = 30.0, entropy correction coefficient δ = 0.05 · c. The grid for the example can be generated using the program described in Section 12.3. It is also provided as n0012.grd in the subdirectory run. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named n0012_input_r.
403
Computational Fluid Dynamics: Principles and Applications
0
0.40 0.35
-1 Convergence history Lift coefficient
-2
0.25 0.20
-3
Lift coefficient
0.30 log(residual)
404
0.15 -4 0.10 -5 0
500
1000
1500
2000
0.05 2500
Iteration
Figure 12.6 Convergence history of the density residual and of the lift coefficient.
Figure 12.7 Mach-number distribution around the airfoil.
Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to run and typing: % Struct2D n0012_input_r
This will produce the convergence history displayed in Fig. 12.6. Figure 12.7 shows the resulting Mach-number contours around the airfoil.
Software Applications
12.7 UNSTRUCTURED 2-D EULER/NAVIER-STOKES SOLVER The directories unstructured_2d and C++/unstructured_2d contain a program for the solution of 2-D Euler and Navier-Stokes equations on unstructured triangular grids. The median-dual cell-vertex scheme described in Section 5.2.2 is utilized for the spatial discretization. The computation of the convective and viscous fluxes, as well as of the gradients employs an edge-based data structure. The program also offers the possibility to simulate incompressible flows at very low Mach numbers (see Section 9.5). Source code of the flow solver is provided in the subdirectory src. The input files, grids, and the results can be found in the subdirectory run. The subdirectory doc contains the documentation (files, modules, functions, variables) in HTML format. Convective fluxes are evaluated according to Roe’s flux-difference splitting scheme (Sections 4.3.3 and 5.3.2). The solution at the faces of the control volume is obtained by piecewise linear reconstruction of the left and right states as given by Eq. (5.42). The gradients of the flow variables are computed with the Green-Gauss approach (Eqs. (5.49) and (5.50)). The reconstructed solution is limited using Venkatakrishnan’s limiter described in Section 5.3.5 (Eq. (5.67)). The gradients at the faces of the control volume, which are required for the evaluation of the viscous fluxes, are obtained by the methodology presented in Section 5.4.2. The discretized governing equations are integrated in time using the explicit multistage scheme (Section 6.1.1 or 6.1.2), enhanced by local time-stepping (Section 9.1) and the central implicit residual smoothing (Section 9.3.2). The same boundary conditions as described previously for the structured flow solver (except the injection and the coordinate cut) are included. The data structure is such that each boundary face could have a different boundary condition. The implementation of the far-field, inlet and outlet boundary conditions employs the concept of dummy nodes (with one layer). All topological data related to the boundaries are stored together with the grid coordinates and the node indexes in one common file (with the extension ugr).
12.7.1 Example run As a test example, let us compute transonic flow through the VKI-1 turbine cascade with the following boundary conditions [4]: • Inlet total pressure pt,in = 1.0 × 105 Pa • Inlet total temperature Tt,in = 300.0 K • Inlet flow angle αinl = 30◦ • Outlet static pressure pout = 5.283 × 104 Pa The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and explicit time-stepping: σ = 5.5, = 0.4, limiter coefficient K = 1.0, entropy correction coefficient δ = 0.05 · c.
405
Computational Fluid Dynamics: Principles and Applications
0
70 65
–1
55
–2
Convergence history Mass flow
50 –3
45
Mass flow (kg/s)
60
log (residual)
406
40 –4 35 –5 0
500
1000
1500
2000
30 2500
Iteration
Figure 12.8 Convergence history of the density residual and of the mass flow.
The grid for the example can be generated using the programs described in Sections 12.3 and 12.4. An unstructured grid for the case generated by the advancing-front method (see Section 11.2.2) is provided as vki1.ugr in the subdirectory run. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named vki1_input. Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to run and typing: % Unstruct2D vki1_input
This will produce the convergence history displayed in Fig. 12.8. Figure 12.9 shows the resulting pressure contours inside the cascade.
12.8 PARALLELIZATION The directory Parallelization contains examples of code parallelized using Asynchronous Agents Library, CUDA, MPI, and OpenMP (Section 9.6). The program to be parallelized solves the 2-D Laplace equation φxx + φyy = 0
(12.6)
Software Applications
Figure 12.9 Pressure contours inside the cascade.
on a structured, rectangular grid using the Jacobi iteration. The boundary values are kept fixed during the iterations. The examples are written in C (CUDA and MPI) and in C++. In the case of MPI, compilation is controlled by a Makefile. In the other cases, project files for Visual Studio (VS 2012) are provided. Note that all examples, except for OpenMP, require the presence of corresponding libraries and compilation tools (see references in Section 9.6). The OpenMP program does not need any special provisions, any recent C++ compiler is sufficient.
REFERENCES [1] Farin GE. Curves and surfaces for computer aided geometric design—a practical guide. 3rd ed. Boston: Academic Press; 1993. [2] Hoschek J, Lasser D. Fundamentals of computer aided geometric design. Wellesley, MA: A.K. Peters Ltd; 1993. [3] Shapiro AH. The dynamics and thermodynamics of compressible fluid flow. New York: Ronald Press; 1953. [4] Kiock R, Lehthaus F, Baines NC, Sieverding CH. The transonic flow through a plane turbine cascade as measured in four European wind tunnels. ASME J Eng Gas Turbines Power 1986;108:277-84.
407
Software Applications Contents 12.1 12.2 12.3 12.4 12.5
Programs for Stability Analysis Structured 1-D Grid Generator Structured 2-D Grid Generators Structured to Unstructured Grid Converter Quasi 1-D Euler Solver 12.5.1 Example run 12.6 Structured 2-D Euler/Navier-Stokes Solver 12.6.1 Example run 12.7 Unstructured 2-D Euler/Navier-Stokes Solver 12.7.1 Example run 12.8 Parallelization References
397 397 398 399 399 400 400 403 405 405 406 407
The software package contains the source codes and executables of several 1-D and 2-D flow solvers and grid generators. Provided are also input datasets and grids for various 1-D and 2-D examples of flow cases. Furthermore, there are two programs for the von Neumann stability analysis of explicit and implicit time-stepping schemes. Finally, you will find examples for code parallelization using different approaches. The aim of the software is to demonstrate how to translate the theoretical principles of the computational fluid dynamics, which were presented in the previous chapters, into a computer code. The programs should be conceived as a basis for further experimentation and enhancements. The source codes are provided under the terms of the GNU General Public License, which means that they may be freely copied and redistributed free of charge. See the file LICENSE.txt in one of the directories for more details. Source codes of the flow solvers and grid generators are written in Fortran-77 and in Fortran-90. The source of the 2-D unstructured flow solver is also provided in modern C++. The programs do not contain any references to external libraries. The source codes are kept as simple as possible, but still flexible enough. No extensive attempts were made to optimize the software for speed or for memory consumption. With a few exceptions, the source codes are organized in such a way that there is just one subroutine or function per file. Sources of the 2-D flow solvers are fully documented using the Doxygen1 tool. All grid, solution, and convergence files are stored in a plain ASCII format. The convergence and the solution files are written out in a form suitable for visualization 1
http://www.doxygen.org.
Computational Fluid Dynamics: Principles and Applications
Copyright © 2015 Elsevier Ltd. All rights reserved.
395
396
Computational Fluid Dynamics: Principles and Applications
with the program Vis2D.2 The format of the convergence and solution files is always the same—there is one column for each variable. The names of the variables are given in the header. In the case of the solution files generated by the unstructured flow solver, the indexes of the nodes of each element are stored after the coordinates and the flow quantities. The first column represents node 1, the second one node 2, etc. The file format is described in detail in the documentation of Vis2D. Because of its simple form, the format can be adapted easily for other visualization tools if necessary. The software package is organized into three highest-level directories: C++, Fortran, and Parallelization. These contain the following subdirectories: • C++ • unstructured_2d—solution of 2-D Euler/Navier-Stokes equations on unstructured triangular grids. • Fortran • analysis—von Neumann stability analysis of 1-D model equations • grid_1d—1-D grid generation (Laval nozzle) • grid_struct_2d—2-D structured grid generation for external and internal flows • grid_unstr_2d—conversion of 2-D structured into unstructured triangular grids • structured_1d—solution of quasi 1-D Euler equations (flow trough Laval nozzle); demonstration of the multigrid method • structured_2d—solution of 2-D Euler/Navier-Stokes equations on structured grids • unstructured_2d—solution of 2-D Euler/Navier-Stokes equations on unstructured triangular grids. • Parallelization • agents—Asynchronous Agents Library • cuda—CUDA • mpi—MPI • openmp—OpenMP • serial—serial version. The contents of the above subdirectories are further explained in the sections below. In general, the subdirectory of each simulation code contains a README.txt file. It describes the particular files, how to compile and to run the application, and how the results are stored. Additionally, the README.txt file explains the meaning of the main variables and of the input parameters. In the case of the 2-D flow solvers, the subdirectory doc contains the documentation in HTML format (index.html). The subdirectory of each flow solver also has the folder run, where you can find input files, grids, solutions, and graphics for the various demonstration cases. 2
Available from CFD Consulting & Analysis (http://www.cfd-ca.de).
Software Applications
The last remark concerns the convergence history. The convergence of all flow solvers provided here is measured as the 2-norm of the difference of the density variable from two consecutive time steps, i.e., N ρ2 = (ρIn+1 − ρIn )2 . (12.1) I=1
The summation is carried out over all N control volumes. It is convenient to normalize the convergence measure in Eq. (12.1) with its value from the first iteration. The programs store the normalized convergence history directly in logarithmic scale.
12.1 PROGRAMS FOR STABILITY ANALYSIS The directory analysis contains two programs for the von Neumann stability analysis (Section 10.3) of linear 1-D model equations. The first program computes the Fourier symbol and the magnitude of the amplification factor for the explicit multistage time-stepping scheme (Section 6.1.1) and the hybrid scheme (Section 6.1.2). The source code is provided in the subdirectory mstage. The second program analyses the damping properties of a general implicit scheme (Section 6.2). It is contained in the subdirectory implicit. Both programs can deal with either the 1-D convection or the 1-D convection-diffusion model equation (Sections 10.3.2 and 10.3.3). The spatial discretization can be conducted either by the central scheme with artificial dissipation, by the 1st-order upwind, or by the 2nd-order upwind scheme, respectively. The same schemes can be also applied to the implicit operator. In the case of the explicit scheme, the effect of the (central or upwind) implicit residual smoothing (Section 9.3) can be investigated as well.
12.2 STRUCTURED 1-D GRID GENERATOR The directory grid_1d contains a program for the generation of 1-D structured grids. Each grid node is associated with a certain area. The area distribution along the x-axis corresponds to the Laval nozzle. The area of the nozzle is evaluated using the relation 1 πx A(x) = 1 + (A1 − 1) 1 + cos for 0 ≤ x ≤ xthr 2 xthr (12.2)
1 π(x − xthr ) A(x) = 1 + (A2 − 1) 1 − cos for xthr < x ≤ 1, 2 1 − xthr where xthr denotes the location of the throat. The length of the nozzle is supposed to be one. Furthermore, in Eq. (12.2) A1 represents the inlet area and A2 the outlet area
397
398
Computational Fluid Dynamics: Principles and Applications
A2
A1
Ai+1/2
Flow 1
2
i
ib2
imax
Control volume
Figure 12.1 Grid and control volume for the 1-D Euler solver; points i = 1 and i = imax are dummy points.
(see the sketch in Fig. 12.1). The area of the throat results from Eq. (12.2) to A∗ (x = xthr ) = 1. The grid generated by this program serves as an input to the quasi 1-D Euler solver from Section 12.5.
12.3 STRUCTURED 2-D GRID GENERATORS The directory grid_struct_2d contains three different programs for the generation of 2-D structured grids for external and internal flows. The first program, which source code is provided in the subdirectory cgrid, generates a C-type grid (see Section 11.1.1) around an airfoil. An example can be seen in Fig. 11.6. The initial grid is generated algebraically by using the linear TFI method Eq. (11.5). Afterward, elliptic PDEs (Section 11.1.3) are employed to produce boundary-orthogonal grid with specified wall spacing. The airfoil contour is approximated by a Bézier spline (see Refs. [1, 2] for an introduction to Bézier splines). The second program is intended for the generation of an algebraic grid inside a 2-D channel with circular bump. The grid points are clustered around the front and the rear point of the bump. The source code can be found in the subdirectory channel. The next program, which is located in the subdirectory hgrid, generates an H-type grid (Section 11.1.1) for a cascade. An example of grid created by the program is presented in Fig. 11.8. The linear TFI technique (Eq. (11.5)) for algebraic grid generation is employed in this case. The contour of the blade is again described by a Bézier spline, which interpolates the given points. The programs cgrid and hgrid rely on the library provided in the subdirectory srccom. This library contains routines for spline interpolation, linear TFI, elliptic grid generation, and for grid stretching.
Software Applications
12.4 STRUCTURED TO UNSTRUCTURED GRID CONVERTER The directory grid_unstr_2d contains a program for the conversion of 2-D structured grids into unstructured triangular grids. The program divides each quadrilateral cell of the structured grid into two triangles by connecting diagonal nodes with the shortest distance. The triangulated grids can be used as an input for the unstructured flow solver in Section 12.7.
12.5 QUASI 1-D EULER SOLVER The directory structured_1d contains a program for the solution of quasi 1-D Euler equations. The equations govern the inviscid flow in a nozzle or in a tube. They can be written in conservative, differential form as [3] ∂W ∂ Fc + = Q. (12.3) ∂t ∂x The vectors of conservative variables, convective fluxes, and the source term read ⎡ ⎡ ⎤ ⎤ ⎡ ⎤ ρuA ρA 0 = ⎣ ρuA ⎦, Fc = ⎣(ρu2 + p)A⎦, Q = ⎣p dA/dx⎦, W (12.4) ρEA ρHuA 0 where A denotes the nozzle area. The total enthalpy H is given by the formula (2.12), and the pressure results according to Eq. (2.29) from u2 p = (γ − 1) ρ E − . (12.5) 2 The governing equations (12.3) are discretized on a structured grid using the dual control-volume methodology (see Section 4.2.3). A sketch of the grid and of the control volume is displayed in Fig. 12.1. Three different schemes are employed for the spatial discretization: • Central scheme with scalar artificial dissipation (Section 4.3.1) • CUSP flux-vector splitting scheme (Section 4.3.2) • Roe’s flux-difference splitting scheme (Section 4.3.3). The discretized equations (12.3) are advanced in time using the explicit multistage scheme (Section 6.1.1 or 6.1.2). The program uses local time-stepping (Section 9.1) and the central implicit residual smoothing technique (Section 9.3.1) for convergence acceleration. The source code can be found in the subdirectory src. The input files and the results are provided in the subdirectory run. A second version of the flow solver employs the multigrid method from Section 9.4 for convergence acceleration. The spatial and temporal discretization schemes are otherwise the same as above. The multigrid scheme uses the FAS formulation for the
399
400
Computational Fluid Dynamics: Principles and Applications
solution on coarse grids and FMG to initialize the flow on the finest grid. The source code is contained in the subdirectory srcmg, exemplary input files can be found in the subdirectory runmg. The boundary conditions at the inlet and the outlet plane are implemented in characteristic variables as presented in Section 8.4. The concept of dummy points, described in Section 8.1, is utilized for this purpose (see Fig. 12.1).
12.5.1 Example run As a test example, let us consider flow through a Laval nozzle with the inlet area A1 = 1.5, outlet area A2 = 2.5, and with the following boundary conditions: • Inlet total pressure pt,in = 1.0 × 105 Pa • Inlet total temperature Tt,in = 288.0 K • Outlet static pressure pout = 7.0 × 104 Pa The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and multigrid with five grid levels: σ = 4.5, = 0.8, limiter coefficient K = 1.5, entropy correction coefficient δ = 0.05 · c. The grid for the example can be generated using the program described in Section 12.2. It is also provided as laval.grd in the subdirectory runmg. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named input_r. Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to runmg and typing: % eul1dmg < input_r
This will produce the convergence history displayed in Fig. 12.2. Figure 12.3 shows the resulting Mach-number distribution along the nozzle length.
12.6 STRUCTURED 2-D EULER/NAVIER-STOKES SOLVER The directory structured_2d contains a program for the solution of 2-D Euler and Navier-Stokes equations on structured body-fitted grids. The spatial discretization is based on the cell-centered finite-volume approach described in Section 4.2.1. It employs the central discretization scheme with scalar artificial dissipation (Section 4.3.1), as well as Roe’s upwind scheme presented in Section 4.3.3. The governing equations are integrated in time using an explicit multistage scheme (Section 6.1.1 or 6.1.2), accelerated by local time-stepping (Section 9.1) and the central implicit residual smoothing (Section 9.3.1). Gradients of the velocity components and of the temperature are computed using the dual control-volume approach described in Section 4.4.1.
Software Applications
340
1
0
320
–1
300
–2
280
–3
260
–4
240
Mass flow (kg/s)
log (residual)
Convergence history Mass flow
220
–5 0
50
100
150
200
250
300
Iteration
Figure 12.2 Convergence history of the density residual and of the mass flow.
2.00
1.75
Mach number
1.50
1.25
1.00
0.75
0.50
0.25 0.0
0.1
0.2
0.3
0.4
0.5 x (m)
Figure 12.3 Mach-number distribution over nozzle length.
0.6
0.7
0.8
0.9
1.0
401
402
Computational Fluid Dynamics: Principles and Applications
The program also offers the possibility to simulate incompressible flows at very low Mach numbers (see Section 9.5). In order to keep the code easy to understand, no turbulence model is implemented. Source code of the flow solver is provided in the subdirectory src. The input files, grids, and the results can be found in the subdirectory run. The subdirectory doc contains the documentation (files, modules, functions, variables) in HTML format. The program utilizes the concept of dummy cells (Section 8.1) for the treatment of the boundary conditions. Two layers of dummy cells are employed. A sketch of the grid in the computational domain is displayed in Fig. 12.4 (cf. Fig. 8.1). The program can deal only with single-block grids. However, it is very flexible with respect to the specification and the type of the boundary conditions. The following eight boundary types are provided: • Coordinate cut (for C- or O-type grids) • Far-field (optionally with vortex correction; sub- and supersonic) • Inflow (subsonic only) • Outflow (sub- and supersonic) • Fluid injection (or extraction) • Line periodic • Solid wall (viscous or inviscid) • Symmetry (along x = const. or y = const. line). The implementation of the above boundary conditions closely follows the discussion in Chapter 8. The four boundaries in the computational space can be divided into an arbitrary number of segments. Each of the segments can be associated with a different boundary J 3 JM JL J2
2
4
2 1 0 0
1
2
1
I2
IL
IM
I
Figure 12.4 Grid in the computational space for the structured 2-D Euler/Navier-Stokes solver. There are two layers of dummy cells at each boundary (dashed line). Numbers in circles denote the sides of the computational domain.
Software Applications
i + 1, j + 1
i, j+1
I,J
SI(I,J)
I+1, J
VOL(I,J) i, j
i +1, j SJ(I,J)
Figure 12.5 Control volume and face vectors for the structured 2-D Euler/ Navier-Stokes solver. Indexes (I, J) denote the cell, indexes (i, j) represent the grid node.
condition. This approach is quite similar to the description of block interfaces presented in Section 8.9. In fact, the program could be relatively easily extended to multiblock grids. Definitions of the segments are stored separately from the grid in the topology file (with the file extension top). The definition of the face vectors SI and SJ (i.e., n · S) employed in the solver can be seen in Fig. 12.5. The face vectors are associated with the left and the bottom face of the control volume VOL(I,J) (corresponds to I,J ). They point outwards of the control volume. The face vectors of the remaining two sides of the control volume are obtained as -SI(I+1,J) and -SJ(I,J+1). In this way, we have to store only two face vectors for each control volume.
12.6.1 Example run As a test example, let us compute transonic flow past the symmetric NACA 0012 airfoil with the following boundary conditions along the far-field: • Mach number M∞ = 0.8 • Angle of attack α = 1.25◦ • Static pressure p∞ = 1.0 × 105 Pa • Static temperature T∞ = 288.0 K The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and explicit time-stepping: σ = 5.0, = 1.2, limiter coefficient K = 30.0, entropy correction coefficient δ = 0.05 · c. The grid for the example can be generated using the program described in Section 12.3. It is also provided as n0012.grd in the subdirectory run. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named n0012_input_r.
403
Computational Fluid Dynamics: Principles and Applications
0
0.40 0.35
-1 Convergence history Lift coefficient
-2
0.25 0.20
-3
Lift coefficient
0.30 log(residual)
404
0.15 -4 0.10 -5 0
500
1000
1500
2000
0.05 2500
Iteration
Figure 12.6 Convergence history of the density residual and of the lift coefficient.
Figure 12.7 Mach-number distribution around the airfoil.
Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to run and typing: % Struct2D n0012_input_r
This will produce the convergence history displayed in Fig. 12.6. Figure 12.7 shows the resulting Mach-number contours around the airfoil.
Software Applications
12.7 UNSTRUCTURED 2-D EULER/NAVIER-STOKES SOLVER The directories unstructured_2d and C++/unstructured_2d contain a program for the solution of 2-D Euler and Navier-Stokes equations on unstructured triangular grids. The median-dual cell-vertex scheme described in Section 5.2.2 is utilized for the spatial discretization. The computation of the convective and viscous fluxes, as well as of the gradients employs an edge-based data structure. The program also offers the possibility to simulate incompressible flows at very low Mach numbers (see Section 9.5). Source code of the flow solver is provided in the subdirectory src. The input files, grids, and the results can be found in the subdirectory run. The subdirectory doc contains the documentation (files, modules, functions, variables) in HTML format. Convective fluxes are evaluated according to Roe’s flux-difference splitting scheme (Sections 4.3.3 and 5.3.2). The solution at the faces of the control volume is obtained by piecewise linear reconstruction of the left and right states as given by Eq. (5.42). The gradients of the flow variables are computed with the Green-Gauss approach (Eqs. (5.49) and (5.50)). The reconstructed solution is limited using Venkatakrishnan’s limiter described in Section 5.3.5 (Eq. (5.67)). The gradients at the faces of the control volume, which are required for the evaluation of the viscous fluxes, are obtained by the methodology presented in Section 5.4.2. The discretized governing equations are integrated in time using the explicit multistage scheme (Section 6.1.1 or 6.1.2), enhanced by local time-stepping (Section 9.1) and the central implicit residual smoothing (Section 9.3.2). The same boundary conditions as described previously for the structured flow solver (except the injection and the coordinate cut) are included. The data structure is such that each boundary face could have a different boundary condition. The implementation of the far-field, inlet and outlet boundary conditions employs the concept of dummy nodes (with one layer). All topological data related to the boundaries are stored together with the grid coordinates and the node indexes in one common file (with the extension ugr).
12.7.1 Example run As a test example, let us compute transonic flow through the VKI-1 turbine cascade with the following boundary conditions [4]: • Inlet total pressure pt,in = 1.0 × 105 Pa • Inlet total temperature Tt,in = 300.0 K • Inlet flow angle αinl = 30◦ • Outlet static pressure pout = 5.283 × 104 Pa The simulation shall be carried out using Roe’s upwind scheme for spatial discretization and explicit time-stepping: σ = 5.5, = 0.4, limiter coefficient K = 1.0, entropy correction coefficient δ = 0.05 · c.
405
Computational Fluid Dynamics: Principles and Applications
0
70 65
–1
55
–2
Convergence history Mass flow
50 –3
45
Mass flow (kg/s)
60
log (residual)
406
40 –4 35 –5 0
500
1000
1500
2000
30 2500
Iteration
Figure 12.8 Convergence history of the density residual and of the mass flow.
The grid for the example can be generated using the programs described in Sections 12.3 and 12.4. An unstructured grid for the case generated by the advancing-front method (see Section 11.2.2) is provided as vki1.ugr in the subdirectory run. The flow solver further requires the input of various parameters in order to set up the boundary conditions, the numerical schemes, as well as the output of convergence history and of plot files. These data are provided in the user input file named vki1_input. Assuming the source code of the solver is already compiled (see README.txt for details), the simulation can be started by opening a terminal window (Command Prompt on Windows), changing to run and typing: % Unstruct2D vki1_input
This will produce the convergence history displayed in Fig. 12.8. Figure 12.9 shows the resulting pressure contours inside the cascade.
12.8 PARALLELIZATION The directory Parallelization contains examples of code parallelized using Asynchronous Agents Library, CUDA, MPI, and OpenMP (Section 9.6). The program to be parallelized solves the 2-D Laplace equation φxx + φyy = 0
(12.6)
Software Applications
Figure 12.9 Pressure contours inside the cascade.
on a structured, rectangular grid using the Jacobi iteration. The boundary values are kept fixed during the iterations. The examples are written in C (CUDA and MPI) and in C++. In the case of MPI, compilation is controlled by a Makefile. In the other cases, project files for Visual Studio (VS 2012) are provided. Note that all examples, except for OpenMP, require the presence of corresponding libraries and compilation tools (see references in Section 9.6). The OpenMP program does not need any special provisions, any recent C++ compiler is sufficient.
REFERENCES [1] Farin GE. Curves and surfaces for computer aided geometric design—a practical guide. 3rd ed. Boston: Academic Press; 1993. [2] Hoschek J, Lasser D. Fundamentals of computer aided geometric design. Wellesley, MA: A.K. Peters Ltd; 1993. [3] Shapiro AH. The dynamics and thermodynamics of compressible fluid flow. New York: Ronald Press; 1953. [4] Kiock R, Lehthaus F, Baines NC, Sieverding CH. The transonic flow through a plane turbine cascade as measured in four European wind tunnels. ASME J Eng Gas Turbines Power 1986;108:277-84.
407