- Email: [email protected]
Parallel Navier–Stokes simulations for high speed compressible flow past arbitrary geometries using FLASH
Computers & Fluids 110 (2015) 27–35
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Parallel Navier–Stokes simulations for high speed compressible flow past arbitrary geometries using FLASH Benzi John ⇑, David R. Emerson, Xiao-Jun Gu Scientific Computing Department, STFC Daresbury Laboratory, Warrington WA4 4AD, United Kingdom
a r t i c l e
i n f o
Article history: Received 15 November 2013 Received in revised form 4 December 2014 Accepted 7 December 2014 Available online 16 December 2014 Keywords: FLASH code Navier–Stokes Hypersonic flow Parallel
a b s t r a c t We report extensions to the FLASH code to enable high-speed compressible viscous flow simulation past arbitrary two- and three-dimensional stationary bodies. The body shape is embedded in a blockstructured Cartesian adaptive mesh refinement grid by implementing appropriate computer graphics algorithms. A high mesh refinement level is required for an accurate body shape representation which results in large grid sizes especially for three-dimensional simulations. Simulations are done in parallel on IBM Blue Gene/Q computing system on which the code performance has been assessed in both pure MPI and hybrid MPI-OpenMP modes. We also implement appropriate wall boundary conditions in FLASH to model viscous-wall effects. Navier–Stokes (NS) solutions for various two-dimensional test cases like a shock–boundary layer interaction problem as well as for hypersonic flow past blunted cone–cylinder–flare and double-cone geometries are shown. Three dimensional NS simulations of micro vortex generators employed in hypersonic flow control have also been carried out and the computed results have been found to be consistent with experimental results. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The FLASH code [1] developed at the Flash center at the University of Chicago is an open-source code with a wide user base. It was originally designed for application to astrophysics problems like super novae, galaxy clusters and stellar structure. More recently the code has also been used to study High Energy Density Physics (HEDP) problems like laser-driven fusion experiments. The FLASH software system [2,3] follows a modular structure consisting of several inter-operable modules like hydrodynamic, material property and nuclear physics solvers that can be combined to solve various problems in cosmology, high energy density physics, etc. The compressible hydrodynamics code has already been validated for typical benchmark problems like the Sod shock tube test case and the classic wind tunnel step problem [4] by solving the inviscid hydrodynamic (Euler) equations. However, it has not yet been applied to practical high-speed CFD applications mainly due to the fact that it relies on a block-structured adaptive mesh refinement scheme using Cartesian cells to generate the grid. The inherent Cartesian grid structure means that special schemes need to be devised to embed geometries of arbitrary shape in the flow domain. Also, FLASH originally being designed as an astrophysics ⇑ Corresponding author. E-mail address: [email protected] (B. John). http://dx.doi.org/10.1016/j.compfluid.2014.12.008 0045-7930/Ó 2014 Elsevier Ltd. All rights reserved.
code is not designed to simulate any viscous wall effects. It has not yet been applied for any compressible Navier–Stokes CFD simulations for flow past arbitrary body shapes, to the best of our knowledge. The FLASH code is modular and extensible which enables users to extend its functionality for their own applications. It has one of the best adaptive mesh compressible hydrodynamics solvers among various open source CFD codes and methodologies to extend its capability for hypersonic flow simulation past arbitrary geometries will be beneficial and of interest to the hypersonic flow community. FLASH is also a scalable parallel code currently featuring both MPI and hybrid MPI-OpenMP modes and stands in good stead in comparison with other codes and parallelization methods [5–8]. The latest version of FLASH (FLASH 4) [1] includes a strategy to incorporate stationary rigid bodies in a computational domain. The solid body is essentially treated as part of the fluid domain and a reflecting boundary condition is applied at the solid/fluid interface. The surface of the rigid body is represented by stair steps due to the regular Cartesian grid structure in FLASH. This scheme, however, can be done only for simple rectangular, spherical or any other shape that can be represented by an analytical expression, which implies that arbitrary complex geometries for real CFD applications cannot be handled currently. In this work, we extend this scheme to generate grids around arbitrary two-dimensional (2D) or three-dimensional (3D) geometric shape by implementing
28
B. John et al. / Computers & Fluids 110 (2015) 27–35
appropriate computer graphics algorithms in FLASH. We also invoke the material property-viscosity module in FLASH to solve the Navier–Stokes equations. Additionally, we implement the noslip wall boundary conditions and Sutherland’s law of viscosity to accurately model viscous-wall effects. A brief discussion of this implementation and preliminary flow results were reported by the authors in [9]. In this work, we elaborate on this and carry out additional Navier–Stokes simulations for several two-dimensional test cases as well as a three-dimensional viscous simulation of micro vortex generators employed in hypersonic flow control. 2. Numerical method The FLASH code has two compressible gas hydrodynamic solvers [1] based on the Finite Volume Method (FVM), built around different operator splitting methods, viz. directionally split and unsplit solvers. They solve the standard compressible Navier–Stokes equations for continuity, momentum and energy, with the pressure field determined from the equation of state [1]. The split solver is based on the piecewise parabolic method (PPM) [10], which is essentially a higher order version of the Godunov scheme. The directionally unsplit solver [11–13] is based on a Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) Hancock type second-order scheme. The unsplit hydro implementation can solve 1D, 2D and 3D problems with added capabilities of exploring various numerical implementations: different types of Riemann solvers; slope limiters; first, second, third and fifth order reconstruction methods as well as a strong shock/rarefaction detection algorithm. One of the notable features of the unsplit hydro scheme is that it particularly improves the preservation of flow symmetries as compared to the splitting formulation. Also, the scheme used in this unsplit algorithm can take a wide range of CFL stability limits for all three dimensions when compared to the directionally split algorithm [11–13]. Grid generation in FLASH is mainly based on a block-structured adaptive mesh refinement (AMR) scheme using PARAMESH [14]. In block-structured AMR, the fundamental data structure is a block of cells arranged in a logically Cartesian fashion, which implies that each cell can be specified using a block identifier (processor number and local block number) and a coordinate triple (i, j, k), where i, j and k refer to cell number in the x-, y-, and z-directions, respectively. PARAMESH handles the filling of guard cells that surrounds each block with information from other blocks or, at the boundaries of the physical domain from an external boundary routine. If the neighbor block has a different level of refinement, the data from the neighbor’s cells is adjusted by either prolongation (interpolation from a coarse to finer level of resolution) or restriction (averaging from a fine to a coarser level). PARAMESH also enforces flux conservation at jumps in refinement across block boundaries, as described by Berger and Colella [15]. In this work, we implement a computer graphics algorithm in FLASH to enable the generation of random complex three dimensional shapes that represents a rigid body. The algorithm is based on a point-in-polyhedron test using spherical polygons proposed by Carvalho and Cavalcanti [16] coded by John Burkardt. The algorithm determines if a given point is inside or outside a threedimensional polyhedron based on a method using spherical polygons. The user needs to define all faces of the 3D geometry by specifying the co-ordinate points of each face listed according to the orientation with respect to the outward normal at that face as input to the algorithm. Each face of the polyhedron is then projected onto a unit sphere and the resulting signed area of the spherical polygon thus formed, determines whether the point is inside or outside the polyhedron. A new variable called bdry_var is defined and is specified as an adaptive mesh refinement (AMR) variable for all the grid points. The grid points of the computational
domain computed to be outside the solid body represents the fluid domain, while those inside represents the solid domain. The values of bdry_var for all grid points within the fluid domain are assigned negative values, while those falling within the solid domain are assigned positive values. An illustration of the demarcation between fluid and solid regions in an AMR Cartesian grid is shown in Fig. 1. The fluid/solid interface represents the wall at which appropriate wall boundary conditions should be imposed. An advantage of this scheme is that mesh generation is fairly simple and extremely quick as there is no time consuming phase associated with surface mesh generation and the associated volume mesh. A limitation of this scheme is that, as the grid is nonbody-fitted, the body shape is represented by stair steps. A high level of refinement is needed to accurately represent a shape which can result in very large grid sizes for 3D simulations. This can be resolved to a great extent by resorting to high performance parallel computing. We have also implemented the no-slip boundary condition at the wall (fluid/solid interface) in FLASH to model viscous-wall effects [17]. Boundary conditions in FLASH need to be handled with the aid of guard cells in each coordinate direction which surrounds each block of local data. The data need to be carefully set such that fluxes at the boundary are physically correct. The ghost cell adjacent to the wall is denoted by g whereas the cell inside the computational domain is denoted by 1. The velocity components (u, v, w) and pressure p in the ghost cells are set as:
ug ¼ 2V w u1
v g ¼ v 1 wg ¼ w1
ð1Þ
pg ¼ þp1 where Vw is the wall velocity. For adiabatic wall boundary conditions, the density gradient in the normal direction vanishes and hence density, q can be specified as
qg ¼ þq1
ð2Þ
For an isothermal wall, density and energy e can be specified as
qg ¼ 2qw q1 eg ¼ 2ew e1
ð3Þ
The Sutherland’s law of viscosity [11] has been implemented in FLASH to model the variation of absolute viscosity with temperature.
Fig. 1. AMR Cartesian grid demarcating the solid domain (denoted in red color) and fluid domain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
B. John et al. / Computers & Fluids 110 (2015) 27–35
3. Computed results and discussion We have carried out CFD simulations for flow past several body shapes that are characteristic for flow in hypersonic flight applications. For all computations reported here, the directionally unsplit hydrodynamic solver in FLASH with a third order accurate data reconstruction scheme based on the PPM method has been employed. The Reimann solver chosen is the HLLC scheme by Toro [18] and the slope limiter employed is based on the Van-Leer scheme. For all the simulations, a very high AMR level corresponding to level = 6 is specified for the bdry_var to get a good representation of the body shape. The same refinement level is also specified for density and pressure (specified as pres_var and dens_ var in FLASH). Our investigations showed that this refinement level is more than adequate to produce grid independent results. 3.1. Laminar shock–boundary layer interaction problem Shock–boundary layer interaction (SBLI) phenomenon is one of the most extensively studied phenomena in compressible fluid dynamics and is typically used as a benchmark case for the validation of Navier–Stokes flow solvers. The computational domain extent for this case is (x, y) e [0, 2.4] [0, 1.164] following the works by Zhang et al. [19] and Jiang et al. [20]. The predicted flow contours which also shows the computational domain for this test case is shown in Fig. 2. The left and upper boundaries are supersonic inflow boundaries and are defined such that an oblique shock is generated with a shock wave angle of 32.6°. The right boundary is defined as a non-reflective supersonic outflow. The bottom boundary is defined as a combination of a symmetry boundary (x < 0.8) and a solid wall (x > 0.8) on which the shock impinges. The free stream left inflow Mach number is 2.0 and the Reynolds number based on the shock impingement length from the wall leading edge is 296,000. The impinging shock, reflecting shock and the boundary layer separation region can be clearly seen from the flow contours. A zoom up of the SBLI and flow separation region with velocity streamlines overlaid on it indicating the presence of a recirculation zone can be visualized from Fig. 2c. The non-dimensional wall pressure and the skin friction coefficient Cf along the wall are plotted and compared with experimental data [21] and other simulation results in the literature [19,20] in Fig. 3. Our results are in good agreement with experimental data and fully consistent with computational results obtained by other researchers for this test case [19–20,22]. 3.2. Hypersonic blunted-cone–cylinder–flare (HB-2) configuration Hypersonic flow past a blunted cone–cylinder–flare (HB-2) configuration is studied in this section. We carry out simulations at two flow conditions for which experimental setup details [23,24] are available in the literature. The free stream Mach number, density and temperature [25,26] for the first case correspond to Ma = 7.5, qa = 0.0091 kg/m3 and Ta = 138.9 K, while those for the second case are Ma = 5 and qa = 0.0051 kg/m3. The flow conditions considered are such that the laminar flow assumption is valid and real gas effects are minimal. A two-dimensional axi-symmetric simulation is considered with the left and right boundaries defined as supersonic inflow and outflow boundaries respectively. The wall is modeled as an isothermal no-slip surface. The cone axis is treated as an axi-symmetry boundary. The computed density on the full computational domain as well as the pressure contours near the blunted cone for the case with Ma = 7.5 is shown in Fig. 4. The strong bow shock upstream of the body as well as a weaker shock over the cylinder–cone
29
juncture that is characteristic for this type of flow is well captured from our computations. The predicted wall pressure distribution non-dimensionalized with respect to the computed stagnation point pressure P0 is compared with experimental data for the two cases and shown in Fig. 5. The computed pressure profile is slightly over predicted for a small region over the cylinder, but in general the predicted pressure distribution is in good agreement with the experimental results. 3.3. Hypersonic double-cone configuration The flow associated with a double cone geometry is very complex as it involves shock–shock and shock–-boundary layer interactions. Moreover, there is an inherent unsteadiness associated with the separation region near the juncture between the two cones. Consequently, the double-cone configuration is considered to be a very challenging case from a computational point of view. Here, we carry out a Navier–Stokes simulation for flow past the double cone geometry at a Mach number of Ma = 12.73 [25,27]. Flow conditions are such that the laminar flow assumption is valid. The free stream conditions correspond to low enthalpy flow conditions [27] and so real gas effects can be reasonably neglected. The free stream density and temperature for this case are qa = 5.9 104 kg/m3 and Ta = 46.1 K respectively. A two-dimensional axi-symmetric simulation is performed with appropriate boundary conditions. The left and right boundaries are defined as supersonic inflow and outflow boundaries respectively, and the wall is modeled as an isothermal no-slip surface. The computed velocity and temperature contours for the full geometry are shown in Fig. 6. The shock wave features and flow separation region at the cone–cone juncture can be visualized from the contour plots. The numerical schlieren calculated from the density gradient is shown in Fig. 7 and this shows the flow features in much more detail. The first cone generates an oblique shock and interacts strongly with the detached bow shock created by the second cone resulting in a transmitted shock impinging on the second cone. The boundary layer separation region near the cone–cone juncture generates its own shock and interacts with the bow shock and this causes a shift in the location of the interaction point and the flow separation region. Also, the schlieren indicates the presence of supersonic jets along the surface of the second cone. All these computed flow features are very much consistent with those expected for hypersonic double-cone flows. The computed wall pressure coefficient, Cp along the extent of the double-cone axis is compared with experimental data [27] in Fig. 8. The predicted pressure profile is non-smooth and in general seems to be shifted toward the right in comparison to the experimental profile. The discrepancy may be due to the complexity of this type of flow combined with the fact that the grid is non-body fitted. Nevertheless, in general the predicted flow features and wall Cp variations from our simulations are qualitatively sound. We also carried out these simulations with a 5th order WENO reconstruction scheme [1], but they produced highly unsteady results and the flow results failed to converge. The reason for this is currently unknown and will be investigated in the future. 3.4. 3D Micro vortex generator (MVG) simulations It is well known that shock–boundary layer interactions cause adverse pressure gradients and boundary layer separation leading to reduced propulsive efficiency in hypersonic vehicles. Micro-ramps are a recently developed passive flow control device which has been shown to suppress the adverse effects of SBLI. The majority of experimental and numerical studies on micro-ramps carried out previously [28–31] have been limited to supersonic flow conditions. Effects of micro-ramps in hypersonic flow control have
30
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 2. Flow results for the shock–boundary layer interaction case: (a) pressure (P/P1) contours, (b) density (q/q1) contours and (c) streamlines overlaid on velocity contours showing recirculation in the flow separation region.
Fig. 3. SBLI validation results: (a) pressure distribution and (b) skin friction distribution along the wall.
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 4. Computed flow contours of (a) density (kg/m3) and (b) pressure (Pa) for the Ma = 7.5 case.
31
distances before being lift-off further downstream. These vortices improve the health of the boundary layer by enhancing flow transport between high shear and low shear regions near the wall, thereby reducing the chance of a flow separation. Additionally, experimental results also show the presence of horseshoe vortices and a pair of secondary vortices. We have carried out 3D Navier– Stokes simulations of this experimental work [32] using the FLASH code. In the experiment, the free stream Mach number, pressure and temperature were Maa = 5, Pa = 1354.93 Pa and Ta = 62.79 K respectively. The computational domain showing the full experimental setup comprising of the shock generator and flat plate with micro ramp is shown in Fig. 9. A high level of AMR (level = 6) was required to accurately represent the body shape, resulting in an initial grid size of about 50 million Cartesian cells. We perform two separate simulations corresponding to the experiments, one with and another without the shock generator. The left inflow boundary is treated as a supersonic inflow boundary, whereas all other domain boundaries as treated as supersonic outflow boundaries. The wall boundaries are treated as no-slip adiabatic walls. First we consider the case without the shock generator. The computed flow streamlines overlaid on velocity contours in the vicinity of the MVG is shown in Fig. 10a. From the streamlines it can be observed that the flow from the top surface of the ramp moves toward the slant edges on both sides in a rotational manner and forms counter-rotating vortices before being eventually lift-off further downstream of the ramp. Downstream of the ramp, the streamlines pass through the vortices which is essentially in a low velocity region often characterized as a momentum deficit region in the literature [28–32]. From the density contours in Fig. 10b traces of a pair of secondary vortices formed immediately downstream of the ramp can be visualized. Density contours along a plane passing through the MVG center is plotted in Fig. 10c from which the leading edge and trailing edge shocks as well the flow wake past the ramp can be seen. All these trends in flow observed from the computed results are characteristic of flow past a MVG and very much consistent with experimental results. We have also performed CFD simulation of the micro-ramp configuration with a shock generator in the form of a wedge. This generates an oblique shock wave that impinges on the flat plate downstream of the MVG resulting in shock–boundary layer interaction. The computed numerical schlieren along a plane passing through the MVG center is shown in Fig. 11. The computed schlieren confirms the presence of the incident shock formed from the shock generator, the SBLI separation region, the subsequent re-attachment shock as well as the leading edge and trailing edge shock generated by the MVG. It is to be noted that we have assumed laminar flow conditions for our MVG simulations as turbulence modeling is beyond the current capability of FLASH. For the case with no shock generator, turbulent effects are not expected to significantly affect the global flow field due to the small size of the micro-ramp and the fact that compressible effects dominate viscous effects at hypersonic speeds. Turbulence modeling, however, needs to be taken into consideration for a more complete study of flow structure in the boundary layer past the MVG, especially for the case with the shock generator. Nevertheless, our simulations have captured all the flow features characteristic of flow past micro-ramps and demonstrated that three dimensional simulations for flow past arbitrary bodies, as in the case of a MVG, can be done effectively using FLASH with the aid of high performance computing.
Fig. 5. Comparison of wall pressure distribution for the case with (a) Ma = 7.5 and (b) Ma = 5.
3.5. Parallel performance of the code
been studied experimentally only very recently [32]. Experimental results show that MVG generates a pair of counter-rotating primary vortices which remain in the boundary layer for relatively long
Parallelization and adaptive mesh refinement (AMR) portion of FLASH is handled mainly by the PARAMESH package. The message passing paradigm is based on MPI and single program multiple
32
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 6. Computed flow contours of (a) velocity (m/s) contours and (b) temperature (K) plotted for the full double-cone geometry.
Fig. 8. Comparison of predicted wall pressure coefficient Cp with experiment. Fig. 7. Computed numerical Schlieren on top-half of the double-cone geometry.
data (SPMD) model. PARAMESH consists of a suite of subroutines which handle distribution of work to processors, guard cell filling, and flux conservation between blocks. FLASH uses various parallel input/output (IO) libraries to simplify and manage the output of the large amounts of data. In the present work, all I/O operations have been handled using parallel HDF5 [33] I/O libraries. The latest version of FLASH also features the capability to run multithreads using OpenMP within each MPI task. However, presently the OpenMP parallel regions do not cover all units in FLASH and is only focused to units which are exercised most often like the various hydrodynamic units. The level of OpenMP threading is expected to increase with future FLASH releases. In the present work, we have assessed the parallel performance of the code on an IBM Blue Gene/Q system called the Blue Joule using both pure MPI and hybrid MPI-OpenMP modes. Blue Joule
is the second fastest supercomputer in the United Kingdom and rated number 23 in the Top500 list of supercomputers published in June 2014. It consists of several racks, each containing 1024 nodes and each such node has 16-core, 64 bit, 1.60 GHz A2 PowerPC processors. Each node is also equipped with 16 GB of system memory, providing a total 96 TB of distributed memory across the system. To assess the scalability of the code in the pure MPI mode, we have considered a FLASH 3D Navier–Stokes micro vortex generator (MVG) simulation with three different block sizes, i.e. 83, 123 and 163 cells per block. The variation of wall-clock time taken with increasing cells per block as a function of number of cores is shown in Fig. 12 for up to 16,384 cores. The code performs consistently well with both increase in load and increase in number of MPI tasks. We have also carried out performance studies using hybrid MPI-OpenMP to assess the potential benefits. The OpenMP threading strategy is based on the fact that the solution updates
B. John et al. / Computers & Fluids 110 (2015) 27–35
33
Fig. 9. Experimental setup of MVG with shock generator.
Fig. 10. Flow in the vicinity of MVG (a) streamlines overlaid on velocity contours (b) density contours downstream of MVG and (c) density contours on an x–y plane passing through the MVG center.
on different block units are independent, enabling multiple threads to safely update the solution on different blocks at the same time. This is included in the FLASH application by adding threadBlockList = True to the FLASH setup line [1]. To assess the parallel performance using hybrid MPI-OpenMP mode, we have considered the
3D MVG simulation with 83 cells per block. The speedup result as a function of the number of OpenMP threads per MPI task is plotted in Fig. 13, for four different mixes of nodes and cores. The performance plot is purely based on the wall-clock time taken for the hydrodynamic part of the code and does not take into
34
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 11. Computed numerical Schlieren along a plane passing through the center of the micro-ramp.
Nevertheless significant savings in computational time can be potentially attained with comparatively lesser number of MPI tasks using the hybrid MPI-OpenMP mode. 4. Conclusions
Fig. 12. Scalability of FLASH code on BGQ.
consideration the time taken for other parts of the code like initialization and I/O which has not been threaded yet. Speedup is calculated as the ratio of the wall-clock time taken for the multithreaded case to the corresponding time taken for the case with one OpenMP thread. The best speedup among the cases considered is about 4 for the case with 1024 cores and 8 threads. One of the reasons for the reduced speedup is that hydrodynamic solvers in FLASH call Paramesh subroutines which have not been multithreaded yet to exchange guard cells and conserve fluxes. Higher number of MPI tasks leads to further loss in performance, which may be due to the increased synchronization time involved.
We have carried out extensions to the FLASH code to enable compressible Navier–Stokes simulations for flow past arbitrary body shapes. The solid body shape is represented in an AMR Cartesian grid by implementing appropriate computer graphics algorithms. A high level of refinement is specified to accurately represent the body shape, which is done with the aid of high performance computing. The advantage with this scheme is that mesh generation is fairly simple and extremely quick, the downside being the approximations at the wall surface due to the fact that the grid is non-body fitted. We also invoke the material property-viscosity module in FLASH to solve the Navier–Stokes equations and implement appropriate boundary conditions at the fluid–solid interface to accurately handle viscous-wall effects. We apply these modifications and carry out several 2D Navier–stokes test simulations to first check FLASH’s predictive capability. The computed results compare well with experimental data for the shock–boundary layer interaction problem as well as for flow past blunted-cone–cylinder flare geometry. For flow past the double cone case, minor discrepancies are noted between the experimental and simulation results. Nevertheless, all complex flow features associated with the double-cone geometry have been captured qualitatively. Finally we carry out 3D Navier–Stokes simulation of micro vortex generator used in hypersonic flow simulation. Our simulations have captured all the essential flow features like the presence of primary, secondary, horseshoe vortices and momentum deficit region downstream of the MVG. The flow simulation and code scalability results demonstrate that FLASH can be used as an effective tool to carry out high-speed compressible Navier–Stokes flow simulation past arbitrary body shapes. The impact of higher-order schemes like the fifth order WENO scheme on the accuracy of the flow results will be investigated in the future. Acknowledgements
Fig. 13. Speedup of the hydrodynamics module in hybrid MPI-OpenMP mode.
The authors would like to thank the Engineering and Physical Sciences Research Council (EPSRC) for its support of Collaborative Computational Project 12 (CCP12). The authors also thank Prof. Konstantinos Kontis and the group at the University of Manchester for providing experimental setup details of the MVG. They would also like to thank Dongwook Lee at Flash Center, University of Chicago for providing FLASH user support and useful comments regarding capabilities of the code.
B. John et al. / Computers & Fluids 110 (2015) 27–35
References [1] FLASH user guide-version 4.0; 2012. [2] Dubey A, Antypas K, Ganapathy MK, Reid LB, Riley KM, Sheeler D, et al. Extensible component based architecture for FLASH, a massively parallel, multiphysics simulation code. Parallel Comput 2009;35:512–22. [3] Fryxell B, Olson K, Ricker P, Timmes FX, Zingale M, Lamb DQ, et al. FLASH: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. Astrophys J Suppl 2000;131:273–334. [4] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73. [5] Gorobets A, Trias FX, Oliva A. A parallel MPI + OpenMP + OpenCL algorithm for hybrid supercomputations of incompressible flows. Comput Fluids 2013;88:764–72. [6] Ye Y, Li K. Entropic Lattice Boltzmann Method based high Reynolds number flow simulation using CUDA on GPU. Comput Fluids 2013;88:241–9. [7] Li K, Yang W, Li K. Performance analysis and optimization for SpMV on GPU using probabilistic modeling. IEEE Trans Parallel Distrib Syst 2014:99. http:// dx.doi.org/10.1109/TPDS.2014.230822. [8] Liang J, Li K. Accelerating solidification process simulation for large-sized system of liquid metal atoms using GPU with CUDA. J Comput Phys 2014;257:521–35. [9] John B, Emerson DR, Gu XJ. Parallel compressible viscous flow simulations using FLASH code: implementation for arbitrary three dimensional geometries. Procedia Eng 2013;61:52–6. [10] Woodward P, Colella P. The piecewise parabolic method (PPM) for gasdynamical simulations. J Comput Phys 1984;54:174–201. [11] Gardiner T, Stone J. An unsplit Godunov method for ideal MHD via constrained transport in three dimensions. J Comput Phys 2008;227:4123–41. [12] Lee D, Deane A. An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics. J Comput Phys 2009;228:952–75. [13] Lee D. A solution accurate, efficient and stable unsplit staggered mesh scheme for three dimensional magnetohydrodynamics. J Comput Phys 2013;243:269–92. [14] Olson K. PARAMESH: a parallel adaptive grid tool. In: Parallel computational fluid dynamics: theory and applications – proceedings of the parallel CFD 2005 conference. Elsevier; 2006. p. 341. [15] Berger MJ, Colella P. Local adaptive mesh refinement for shock hydrodynamics. J Comput Phys 1989;82:64–84.
35
[16] Carvalho P, Cavalcanti P. Point in polyhedron testing using spherical polygons. Graphics Gems V. Academic Press; 1995. [17] White FM. Viscous fluid flow. McGraw-Hill; 2005. [18] Toro EF. Riemann solvers and numerical methods for fluid dynamics, a practical introduction. Springer; 2009. [19] Zhang M, Yu J, Chang SC. Solving the Navier–Stokes equations by the CESE method. AIAA Paper 2004-0075; 2004. [20] Jiang C, Cui S, Feng X. Solving the Euler and Navier–Stokes equations by the AMR-CESE method. Comput Fluids 2012;54:105–17. [21] Hakkinen RJ, Greber I, Trilling L, Abarbanel SS. The interaction of an oblique shock wave with a laminar boundary layer. NASA Memo 2-18-59W; 1959. [22] Yao Y, Krishnan L, Sandham ND, Roberts G. The effect of Mach number on unstable disturbances in shock/boundary-layer interactions. Phys Fluids 2007;19:054104. [23] Gray J, Lindsay E. Force tests of standard hypervelocity ballistic models HB-1 and HB-2 at Mach 1.5 to 10, AEDC-TDR-64-246; 1964. [24] Gilchrist AR, Williams MJ. Pressure distributions and forces on AGARD models HB1 and HB2 at M equals 7.5. Tech. Rep., Aerodynamics Note 346, ARL; 1974. [25] Tissera S. Assessment of high-resolution methods in hypersonic real-gas flows. PhD Thesis. Cranfield University; 2010. [26] Birch T, Prince S, Ludlow DK, Qin N. The application of a parabolized Navier– Stokes solver to some hypersonic flow problems. AIAA paper 2001-1753; 2001. [27] Coblish JJ, Smith MS, Hand T, Candler GV, Nompelis I. Double-cone experiment and numerical analysis at AEDC hypervelocity wind tunnel no. 9. AIAA paper 2005-0902; 2005. [28] Babinsky H, Li Y, Ford CWP. Microramp control of supersonic oblique shockwave/boundary-layer interactions. AIAA J 2009;47:668–75. [29] Li Q, Liu C. LES for supersonic ramp control flow using MVG at M = 2.5 and Reh = 1440. AIAA Paper 2010-0592; 2010. [30] Holden H, Babinsky H. Effect of microvortex generators on separated normal shock/boundary layer interactions. J Aircr 2007;44:170–4. [31] Anderson BH, Tinapple J, Surber L. Optimal control of shock wave turbulent boundary layer interactions using micro-array actuation. AIAA Paper 20063197; 2006. [32] Mohd R, Hossein Z, Azam C, Konstantinos K. Micro-Ramps for hypersonic flow control. Micromachines 2012;3:364–78. [33] HDF5 user’s guide. Version 1.8.7; 2011.
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Parallel Navier–Stokes simulations for high speed compressible flow past arbitrary geometries using FLASH Benzi John ⇑, David R. Emerson, Xiao-Jun Gu Scientific Computing Department, STFC Daresbury Laboratory, Warrington WA4 4AD, United Kingdom
a r t i c l e
i n f o
Article history: Received 15 November 2013 Received in revised form 4 December 2014 Accepted 7 December 2014 Available online 16 December 2014 Keywords: FLASH code Navier–Stokes Hypersonic flow Parallel
a b s t r a c t We report extensions to the FLASH code to enable high-speed compressible viscous flow simulation past arbitrary two- and three-dimensional stationary bodies. The body shape is embedded in a blockstructured Cartesian adaptive mesh refinement grid by implementing appropriate computer graphics algorithms. A high mesh refinement level is required for an accurate body shape representation which results in large grid sizes especially for three-dimensional simulations. Simulations are done in parallel on IBM Blue Gene/Q computing system on which the code performance has been assessed in both pure MPI and hybrid MPI-OpenMP modes. We also implement appropriate wall boundary conditions in FLASH to model viscous-wall effects. Navier–Stokes (NS) solutions for various two-dimensional test cases like a shock–boundary layer interaction problem as well as for hypersonic flow past blunted cone–cylinder–flare and double-cone geometries are shown. Three dimensional NS simulations of micro vortex generators employed in hypersonic flow control have also been carried out and the computed results have been found to be consistent with experimental results. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The FLASH code [1] developed at the Flash center at the University of Chicago is an open-source code with a wide user base. It was originally designed for application to astrophysics problems like super novae, galaxy clusters and stellar structure. More recently the code has also been used to study High Energy Density Physics (HEDP) problems like laser-driven fusion experiments. The FLASH software system [2,3] follows a modular structure consisting of several inter-operable modules like hydrodynamic, material property and nuclear physics solvers that can be combined to solve various problems in cosmology, high energy density physics, etc. The compressible hydrodynamics code has already been validated for typical benchmark problems like the Sod shock tube test case and the classic wind tunnel step problem [4] by solving the inviscid hydrodynamic (Euler) equations. However, it has not yet been applied to practical high-speed CFD applications mainly due to the fact that it relies on a block-structured adaptive mesh refinement scheme using Cartesian cells to generate the grid. The inherent Cartesian grid structure means that special schemes need to be devised to embed geometries of arbitrary shape in the flow domain. Also, FLASH originally being designed as an astrophysics ⇑ Corresponding author. E-mail address: [email protected] (B. John). http://dx.doi.org/10.1016/j.compfluid.2014.12.008 0045-7930/Ó 2014 Elsevier Ltd. All rights reserved.
code is not designed to simulate any viscous wall effects. It has not yet been applied for any compressible Navier–Stokes CFD simulations for flow past arbitrary body shapes, to the best of our knowledge. The FLASH code is modular and extensible which enables users to extend its functionality for their own applications. It has one of the best adaptive mesh compressible hydrodynamics solvers among various open source CFD codes and methodologies to extend its capability for hypersonic flow simulation past arbitrary geometries will be beneficial and of interest to the hypersonic flow community. FLASH is also a scalable parallel code currently featuring both MPI and hybrid MPI-OpenMP modes and stands in good stead in comparison with other codes and parallelization methods [5–8]. The latest version of FLASH (FLASH 4) [1] includes a strategy to incorporate stationary rigid bodies in a computational domain. The solid body is essentially treated as part of the fluid domain and a reflecting boundary condition is applied at the solid/fluid interface. The surface of the rigid body is represented by stair steps due to the regular Cartesian grid structure in FLASH. This scheme, however, can be done only for simple rectangular, spherical or any other shape that can be represented by an analytical expression, which implies that arbitrary complex geometries for real CFD applications cannot be handled currently. In this work, we extend this scheme to generate grids around arbitrary two-dimensional (2D) or three-dimensional (3D) geometric shape by implementing
28
B. John et al. / Computers & Fluids 110 (2015) 27–35
appropriate computer graphics algorithms in FLASH. We also invoke the material property-viscosity module in FLASH to solve the Navier–Stokes equations. Additionally, we implement the noslip wall boundary conditions and Sutherland’s law of viscosity to accurately model viscous-wall effects. A brief discussion of this implementation and preliminary flow results were reported by the authors in [9]. In this work, we elaborate on this and carry out additional Navier–Stokes simulations for several two-dimensional test cases as well as a three-dimensional viscous simulation of micro vortex generators employed in hypersonic flow control. 2. Numerical method The FLASH code has two compressible gas hydrodynamic solvers [1] based on the Finite Volume Method (FVM), built around different operator splitting methods, viz. directionally split and unsplit solvers. They solve the standard compressible Navier–Stokes equations for continuity, momentum and energy, with the pressure field determined from the equation of state [1]. The split solver is based on the piecewise parabolic method (PPM) [10], which is essentially a higher order version of the Godunov scheme. The directionally unsplit solver [11–13] is based on a Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) Hancock type second-order scheme. The unsplit hydro implementation can solve 1D, 2D and 3D problems with added capabilities of exploring various numerical implementations: different types of Riemann solvers; slope limiters; first, second, third and fifth order reconstruction methods as well as a strong shock/rarefaction detection algorithm. One of the notable features of the unsplit hydro scheme is that it particularly improves the preservation of flow symmetries as compared to the splitting formulation. Also, the scheme used in this unsplit algorithm can take a wide range of CFL stability limits for all three dimensions when compared to the directionally split algorithm [11–13]. Grid generation in FLASH is mainly based on a block-structured adaptive mesh refinement (AMR) scheme using PARAMESH [14]. In block-structured AMR, the fundamental data structure is a block of cells arranged in a logically Cartesian fashion, which implies that each cell can be specified using a block identifier (processor number and local block number) and a coordinate triple (i, j, k), where i, j and k refer to cell number in the x-, y-, and z-directions, respectively. PARAMESH handles the filling of guard cells that surrounds each block with information from other blocks or, at the boundaries of the physical domain from an external boundary routine. If the neighbor block has a different level of refinement, the data from the neighbor’s cells is adjusted by either prolongation (interpolation from a coarse to finer level of resolution) or restriction (averaging from a fine to a coarser level). PARAMESH also enforces flux conservation at jumps in refinement across block boundaries, as described by Berger and Colella [15]. In this work, we implement a computer graphics algorithm in FLASH to enable the generation of random complex three dimensional shapes that represents a rigid body. The algorithm is based on a point-in-polyhedron test using spherical polygons proposed by Carvalho and Cavalcanti [16] coded by John Burkardt. The algorithm determines if a given point is inside or outside a threedimensional polyhedron based on a method using spherical polygons. The user needs to define all faces of the 3D geometry by specifying the co-ordinate points of each face listed according to the orientation with respect to the outward normal at that face as input to the algorithm. Each face of the polyhedron is then projected onto a unit sphere and the resulting signed area of the spherical polygon thus formed, determines whether the point is inside or outside the polyhedron. A new variable called bdry_var is defined and is specified as an adaptive mesh refinement (AMR) variable for all the grid points. The grid points of the computational
domain computed to be outside the solid body represents the fluid domain, while those inside represents the solid domain. The values of bdry_var for all grid points within the fluid domain are assigned negative values, while those falling within the solid domain are assigned positive values. An illustration of the demarcation between fluid and solid regions in an AMR Cartesian grid is shown in Fig. 1. The fluid/solid interface represents the wall at which appropriate wall boundary conditions should be imposed. An advantage of this scheme is that mesh generation is fairly simple and extremely quick as there is no time consuming phase associated with surface mesh generation and the associated volume mesh. A limitation of this scheme is that, as the grid is nonbody-fitted, the body shape is represented by stair steps. A high level of refinement is needed to accurately represent a shape which can result in very large grid sizes for 3D simulations. This can be resolved to a great extent by resorting to high performance parallel computing. We have also implemented the no-slip boundary condition at the wall (fluid/solid interface) in FLASH to model viscous-wall effects [17]. Boundary conditions in FLASH need to be handled with the aid of guard cells in each coordinate direction which surrounds each block of local data. The data need to be carefully set such that fluxes at the boundary are physically correct. The ghost cell adjacent to the wall is denoted by g whereas the cell inside the computational domain is denoted by 1. The velocity components (u, v, w) and pressure p in the ghost cells are set as:
ug ¼ 2V w u1
v g ¼ v 1 wg ¼ w1
ð1Þ
pg ¼ þp1 where Vw is the wall velocity. For adiabatic wall boundary conditions, the density gradient in the normal direction vanishes and hence density, q can be specified as
qg ¼ þq1
ð2Þ
For an isothermal wall, density and energy e can be specified as
qg ¼ 2qw q1 eg ¼ 2ew e1
ð3Þ
The Sutherland’s law of viscosity [11] has been implemented in FLASH to model the variation of absolute viscosity with temperature.
Fig. 1. AMR Cartesian grid demarcating the solid domain (denoted in red color) and fluid domain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
B. John et al. / Computers & Fluids 110 (2015) 27–35
3. Computed results and discussion We have carried out CFD simulations for flow past several body shapes that are characteristic for flow in hypersonic flight applications. For all computations reported here, the directionally unsplit hydrodynamic solver in FLASH with a third order accurate data reconstruction scheme based on the PPM method has been employed. The Reimann solver chosen is the HLLC scheme by Toro [18] and the slope limiter employed is based on the Van-Leer scheme. For all the simulations, a very high AMR level corresponding to level = 6 is specified for the bdry_var to get a good representation of the body shape. The same refinement level is also specified for density and pressure (specified as pres_var and dens_ var in FLASH). Our investigations showed that this refinement level is more than adequate to produce grid independent results. 3.1. Laminar shock–boundary layer interaction problem Shock–boundary layer interaction (SBLI) phenomenon is one of the most extensively studied phenomena in compressible fluid dynamics and is typically used as a benchmark case for the validation of Navier–Stokes flow solvers. The computational domain extent for this case is (x, y) e [0, 2.4] [0, 1.164] following the works by Zhang et al. [19] and Jiang et al. [20]. The predicted flow contours which also shows the computational domain for this test case is shown in Fig. 2. The left and upper boundaries are supersonic inflow boundaries and are defined such that an oblique shock is generated with a shock wave angle of 32.6°. The right boundary is defined as a non-reflective supersonic outflow. The bottom boundary is defined as a combination of a symmetry boundary (x < 0.8) and a solid wall (x > 0.8) on which the shock impinges. The free stream left inflow Mach number is 2.0 and the Reynolds number based on the shock impingement length from the wall leading edge is 296,000. The impinging shock, reflecting shock and the boundary layer separation region can be clearly seen from the flow contours. A zoom up of the SBLI and flow separation region with velocity streamlines overlaid on it indicating the presence of a recirculation zone can be visualized from Fig. 2c. The non-dimensional wall pressure and the skin friction coefficient Cf along the wall are plotted and compared with experimental data [21] and other simulation results in the literature [19,20] in Fig. 3. Our results are in good agreement with experimental data and fully consistent with computational results obtained by other researchers for this test case [19–20,22]. 3.2. Hypersonic blunted-cone–cylinder–flare (HB-2) configuration Hypersonic flow past a blunted cone–cylinder–flare (HB-2) configuration is studied in this section. We carry out simulations at two flow conditions for which experimental setup details [23,24] are available in the literature. The free stream Mach number, density and temperature [25,26] for the first case correspond to Ma = 7.5, qa = 0.0091 kg/m3 and Ta = 138.9 K, while those for the second case are Ma = 5 and qa = 0.0051 kg/m3. The flow conditions considered are such that the laminar flow assumption is valid and real gas effects are minimal. A two-dimensional axi-symmetric simulation is considered with the left and right boundaries defined as supersonic inflow and outflow boundaries respectively. The wall is modeled as an isothermal no-slip surface. The cone axis is treated as an axi-symmetry boundary. The computed density on the full computational domain as well as the pressure contours near the blunted cone for the case with Ma = 7.5 is shown in Fig. 4. The strong bow shock upstream of the body as well as a weaker shock over the cylinder–cone
29
juncture that is characteristic for this type of flow is well captured from our computations. The predicted wall pressure distribution non-dimensionalized with respect to the computed stagnation point pressure P0 is compared with experimental data for the two cases and shown in Fig. 5. The computed pressure profile is slightly over predicted for a small region over the cylinder, but in general the predicted pressure distribution is in good agreement with the experimental results. 3.3. Hypersonic double-cone configuration The flow associated with a double cone geometry is very complex as it involves shock–shock and shock–-boundary layer interactions. Moreover, there is an inherent unsteadiness associated with the separation region near the juncture between the two cones. Consequently, the double-cone configuration is considered to be a very challenging case from a computational point of view. Here, we carry out a Navier–Stokes simulation for flow past the double cone geometry at a Mach number of Ma = 12.73 [25,27]. Flow conditions are such that the laminar flow assumption is valid. The free stream conditions correspond to low enthalpy flow conditions [27] and so real gas effects can be reasonably neglected. The free stream density and temperature for this case are qa = 5.9 104 kg/m3 and Ta = 46.1 K respectively. A two-dimensional axi-symmetric simulation is performed with appropriate boundary conditions. The left and right boundaries are defined as supersonic inflow and outflow boundaries respectively, and the wall is modeled as an isothermal no-slip surface. The computed velocity and temperature contours for the full geometry are shown in Fig. 6. The shock wave features and flow separation region at the cone–cone juncture can be visualized from the contour plots. The numerical schlieren calculated from the density gradient is shown in Fig. 7 and this shows the flow features in much more detail. The first cone generates an oblique shock and interacts strongly with the detached bow shock created by the second cone resulting in a transmitted shock impinging on the second cone. The boundary layer separation region near the cone–cone juncture generates its own shock and interacts with the bow shock and this causes a shift in the location of the interaction point and the flow separation region. Also, the schlieren indicates the presence of supersonic jets along the surface of the second cone. All these computed flow features are very much consistent with those expected for hypersonic double-cone flows. The computed wall pressure coefficient, Cp along the extent of the double-cone axis is compared with experimental data [27] in Fig. 8. The predicted pressure profile is non-smooth and in general seems to be shifted toward the right in comparison to the experimental profile. The discrepancy may be due to the complexity of this type of flow combined with the fact that the grid is non-body fitted. Nevertheless, in general the predicted flow features and wall Cp variations from our simulations are qualitatively sound. We also carried out these simulations with a 5th order WENO reconstruction scheme [1], but they produced highly unsteady results and the flow results failed to converge. The reason for this is currently unknown and will be investigated in the future. 3.4. 3D Micro vortex generator (MVG) simulations It is well known that shock–boundary layer interactions cause adverse pressure gradients and boundary layer separation leading to reduced propulsive efficiency in hypersonic vehicles. Micro-ramps are a recently developed passive flow control device which has been shown to suppress the adverse effects of SBLI. The majority of experimental and numerical studies on micro-ramps carried out previously [28–31] have been limited to supersonic flow conditions. Effects of micro-ramps in hypersonic flow control have
30
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 2. Flow results for the shock–boundary layer interaction case: (a) pressure (P/P1) contours, (b) density (q/q1) contours and (c) streamlines overlaid on velocity contours showing recirculation in the flow separation region.
Fig. 3. SBLI validation results: (a) pressure distribution and (b) skin friction distribution along the wall.
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 4. Computed flow contours of (a) density (kg/m3) and (b) pressure (Pa) for the Ma = 7.5 case.
31
distances before being lift-off further downstream. These vortices improve the health of the boundary layer by enhancing flow transport between high shear and low shear regions near the wall, thereby reducing the chance of a flow separation. Additionally, experimental results also show the presence of horseshoe vortices and a pair of secondary vortices. We have carried out 3D Navier– Stokes simulations of this experimental work [32] using the FLASH code. In the experiment, the free stream Mach number, pressure and temperature were Maa = 5, Pa = 1354.93 Pa and Ta = 62.79 K respectively. The computational domain showing the full experimental setup comprising of the shock generator and flat plate with micro ramp is shown in Fig. 9. A high level of AMR (level = 6) was required to accurately represent the body shape, resulting in an initial grid size of about 50 million Cartesian cells. We perform two separate simulations corresponding to the experiments, one with and another without the shock generator. The left inflow boundary is treated as a supersonic inflow boundary, whereas all other domain boundaries as treated as supersonic outflow boundaries. The wall boundaries are treated as no-slip adiabatic walls. First we consider the case without the shock generator. The computed flow streamlines overlaid on velocity contours in the vicinity of the MVG is shown in Fig. 10a. From the streamlines it can be observed that the flow from the top surface of the ramp moves toward the slant edges on both sides in a rotational manner and forms counter-rotating vortices before being eventually lift-off further downstream of the ramp. Downstream of the ramp, the streamlines pass through the vortices which is essentially in a low velocity region often characterized as a momentum deficit region in the literature [28–32]. From the density contours in Fig. 10b traces of a pair of secondary vortices formed immediately downstream of the ramp can be visualized. Density contours along a plane passing through the MVG center is plotted in Fig. 10c from which the leading edge and trailing edge shocks as well the flow wake past the ramp can be seen. All these trends in flow observed from the computed results are characteristic of flow past a MVG and very much consistent with experimental results. We have also performed CFD simulation of the micro-ramp configuration with a shock generator in the form of a wedge. This generates an oblique shock wave that impinges on the flat plate downstream of the MVG resulting in shock–boundary layer interaction. The computed numerical schlieren along a plane passing through the MVG center is shown in Fig. 11. The computed schlieren confirms the presence of the incident shock formed from the shock generator, the SBLI separation region, the subsequent re-attachment shock as well as the leading edge and trailing edge shock generated by the MVG. It is to be noted that we have assumed laminar flow conditions for our MVG simulations as turbulence modeling is beyond the current capability of FLASH. For the case with no shock generator, turbulent effects are not expected to significantly affect the global flow field due to the small size of the micro-ramp and the fact that compressible effects dominate viscous effects at hypersonic speeds. Turbulence modeling, however, needs to be taken into consideration for a more complete study of flow structure in the boundary layer past the MVG, especially for the case with the shock generator. Nevertheless, our simulations have captured all the flow features characteristic of flow past micro-ramps and demonstrated that three dimensional simulations for flow past arbitrary bodies, as in the case of a MVG, can be done effectively using FLASH with the aid of high performance computing.
Fig. 5. Comparison of wall pressure distribution for the case with (a) Ma = 7.5 and (b) Ma = 5.
3.5. Parallel performance of the code
been studied experimentally only very recently [32]. Experimental results show that MVG generates a pair of counter-rotating primary vortices which remain in the boundary layer for relatively long
Parallelization and adaptive mesh refinement (AMR) portion of FLASH is handled mainly by the PARAMESH package. The message passing paradigm is based on MPI and single program multiple
32
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 6. Computed flow contours of (a) velocity (m/s) contours and (b) temperature (K) plotted for the full double-cone geometry.
Fig. 8. Comparison of predicted wall pressure coefficient Cp with experiment. Fig. 7. Computed numerical Schlieren on top-half of the double-cone geometry.
data (SPMD) model. PARAMESH consists of a suite of subroutines which handle distribution of work to processors, guard cell filling, and flux conservation between blocks. FLASH uses various parallel input/output (IO) libraries to simplify and manage the output of the large amounts of data. In the present work, all I/O operations have been handled using parallel HDF5 [33] I/O libraries. The latest version of FLASH also features the capability to run multithreads using OpenMP within each MPI task. However, presently the OpenMP parallel regions do not cover all units in FLASH and is only focused to units which are exercised most often like the various hydrodynamic units. The level of OpenMP threading is expected to increase with future FLASH releases. In the present work, we have assessed the parallel performance of the code on an IBM Blue Gene/Q system called the Blue Joule using both pure MPI and hybrid MPI-OpenMP modes. Blue Joule
is the second fastest supercomputer in the United Kingdom and rated number 23 in the Top500 list of supercomputers published in June 2014. It consists of several racks, each containing 1024 nodes and each such node has 16-core, 64 bit, 1.60 GHz A2 PowerPC processors. Each node is also equipped with 16 GB of system memory, providing a total 96 TB of distributed memory across the system. To assess the scalability of the code in the pure MPI mode, we have considered a FLASH 3D Navier–Stokes micro vortex generator (MVG) simulation with three different block sizes, i.e. 83, 123 and 163 cells per block. The variation of wall-clock time taken with increasing cells per block as a function of number of cores is shown in Fig. 12 for up to 16,384 cores. The code performs consistently well with both increase in load and increase in number of MPI tasks. We have also carried out performance studies using hybrid MPI-OpenMP to assess the potential benefits. The OpenMP threading strategy is based on the fact that the solution updates
B. John et al. / Computers & Fluids 110 (2015) 27–35
33
Fig. 9. Experimental setup of MVG with shock generator.
Fig. 10. Flow in the vicinity of MVG (a) streamlines overlaid on velocity contours (b) density contours downstream of MVG and (c) density contours on an x–y plane passing through the MVG center.
on different block units are independent, enabling multiple threads to safely update the solution on different blocks at the same time. This is included in the FLASH application by adding threadBlockList = True to the FLASH setup line [1]. To assess the parallel performance using hybrid MPI-OpenMP mode, we have considered the
3D MVG simulation with 83 cells per block. The speedup result as a function of the number of OpenMP threads per MPI task is plotted in Fig. 13, for four different mixes of nodes and cores. The performance plot is purely based on the wall-clock time taken for the hydrodynamic part of the code and does not take into
34
B. John et al. / Computers & Fluids 110 (2015) 27–35
Fig. 11. Computed numerical Schlieren along a plane passing through the center of the micro-ramp.
Nevertheless significant savings in computational time can be potentially attained with comparatively lesser number of MPI tasks using the hybrid MPI-OpenMP mode. 4. Conclusions
Fig. 12. Scalability of FLASH code on BGQ.
consideration the time taken for other parts of the code like initialization and I/O which has not been threaded yet. Speedup is calculated as the ratio of the wall-clock time taken for the multithreaded case to the corresponding time taken for the case with one OpenMP thread. The best speedup among the cases considered is about 4 for the case with 1024 cores and 8 threads. One of the reasons for the reduced speedup is that hydrodynamic solvers in FLASH call Paramesh subroutines which have not been multithreaded yet to exchange guard cells and conserve fluxes. Higher number of MPI tasks leads to further loss in performance, which may be due to the increased synchronization time involved.
We have carried out extensions to the FLASH code to enable compressible Navier–Stokes simulations for flow past arbitrary body shapes. The solid body shape is represented in an AMR Cartesian grid by implementing appropriate computer graphics algorithms. A high level of refinement is specified to accurately represent the body shape, which is done with the aid of high performance computing. The advantage with this scheme is that mesh generation is fairly simple and extremely quick, the downside being the approximations at the wall surface due to the fact that the grid is non-body fitted. We also invoke the material property-viscosity module in FLASH to solve the Navier–Stokes equations and implement appropriate boundary conditions at the fluid–solid interface to accurately handle viscous-wall effects. We apply these modifications and carry out several 2D Navier–stokes test simulations to first check FLASH’s predictive capability. The computed results compare well with experimental data for the shock–boundary layer interaction problem as well as for flow past blunted-cone–cylinder flare geometry. For flow past the double cone case, minor discrepancies are noted between the experimental and simulation results. Nevertheless, all complex flow features associated with the double-cone geometry have been captured qualitatively. Finally we carry out 3D Navier–Stokes simulation of micro vortex generator used in hypersonic flow simulation. Our simulations have captured all the essential flow features like the presence of primary, secondary, horseshoe vortices and momentum deficit region downstream of the MVG. The flow simulation and code scalability results demonstrate that FLASH can be used as an effective tool to carry out high-speed compressible Navier–Stokes flow simulation past arbitrary body shapes. The impact of higher-order schemes like the fifth order WENO scheme on the accuracy of the flow results will be investigated in the future. Acknowledgements
Fig. 13. Speedup of the hydrodynamics module in hybrid MPI-OpenMP mode.
The authors would like to thank the Engineering and Physical Sciences Research Council (EPSRC) for its support of Collaborative Computational Project 12 (CCP12). The authors also thank Prof. Konstantinos Kontis and the group at the University of Manchester for providing experimental setup details of the MVG. They would also like to thank Dongwook Lee at Flash Center, University of Chicago for providing FLASH user support and useful comments regarding capabilities of the code.
B. John et al. / Computers & Fluids 110 (2015) 27–35
References [1] FLASH user guide-version 4.0; 2012. [2] Dubey A, Antypas K, Ganapathy MK, Reid LB, Riley KM, Sheeler D, et al. Extensible component based architecture for FLASH, a massively parallel, multiphysics simulation code. Parallel Comput 2009;35:512–22. [3] Fryxell B, Olson K, Ricker P, Timmes FX, Zingale M, Lamb DQ, et al. FLASH: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. Astrophys J Suppl 2000;131:273–334. [4] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73. [5] Gorobets A, Trias FX, Oliva A. A parallel MPI + OpenMP + OpenCL algorithm for hybrid supercomputations of incompressible flows. Comput Fluids 2013;88:764–72. [6] Ye Y, Li K. Entropic Lattice Boltzmann Method based high Reynolds number flow simulation using CUDA on GPU. Comput Fluids 2013;88:241–9. [7] Li K, Yang W, Li K. Performance analysis and optimization for SpMV on GPU using probabilistic modeling. IEEE Trans Parallel Distrib Syst 2014:99. http:// dx.doi.org/10.1109/TPDS.2014.230822. [8] Liang J, Li K. Accelerating solidification process simulation for large-sized system of liquid metal atoms using GPU with CUDA. J Comput Phys 2014;257:521–35. [9] John B, Emerson DR, Gu XJ. Parallel compressible viscous flow simulations using FLASH code: implementation for arbitrary three dimensional geometries. Procedia Eng 2013;61:52–6. [10] Woodward P, Colella P. The piecewise parabolic method (PPM) for gasdynamical simulations. J Comput Phys 1984;54:174–201. [11] Gardiner T, Stone J. An unsplit Godunov method for ideal MHD via constrained transport in three dimensions. J Comput Phys 2008;227:4123–41. [12] Lee D, Deane A. An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics. J Comput Phys 2009;228:952–75. [13] Lee D. A solution accurate, efficient and stable unsplit staggered mesh scheme for three dimensional magnetohydrodynamics. J Comput Phys 2013;243:269–92. [14] Olson K. PARAMESH: a parallel adaptive grid tool. In: Parallel computational fluid dynamics: theory and applications – proceedings of the parallel CFD 2005 conference. Elsevier; 2006. p. 341. [15] Berger MJ, Colella P. Local adaptive mesh refinement for shock hydrodynamics. J Comput Phys 1989;82:64–84.
35
[16] Carvalho P, Cavalcanti P. Point in polyhedron testing using spherical polygons. Graphics Gems V. Academic Press; 1995. [17] White FM. Viscous fluid flow. McGraw-Hill; 2005. [18] Toro EF. Riemann solvers and numerical methods for fluid dynamics, a practical introduction. Springer; 2009. [19] Zhang M, Yu J, Chang SC. Solving the Navier–Stokes equations by the CESE method. AIAA Paper 2004-0075; 2004. [20] Jiang C, Cui S, Feng X. Solving the Euler and Navier–Stokes equations by the AMR-CESE method. Comput Fluids 2012;54:105–17. [21] Hakkinen RJ, Greber I, Trilling L, Abarbanel SS. The interaction of an oblique shock wave with a laminar boundary layer. NASA Memo 2-18-59W; 1959. [22] Yao Y, Krishnan L, Sandham ND, Roberts G. The effect of Mach number on unstable disturbances in shock/boundary-layer interactions. Phys Fluids 2007;19:054104. [23] Gray J, Lindsay E. Force tests of standard hypervelocity ballistic models HB-1 and HB-2 at Mach 1.5 to 10, AEDC-TDR-64-246; 1964. [24] Gilchrist AR, Williams MJ. Pressure distributions and forces on AGARD models HB1 and HB2 at M equals 7.5. Tech. Rep., Aerodynamics Note 346, ARL; 1974. [25] Tissera S. Assessment of high-resolution methods in hypersonic real-gas flows. PhD Thesis. Cranfield University; 2010. [26] Birch T, Prince S, Ludlow DK, Qin N. The application of a parabolized Navier– Stokes solver to some hypersonic flow problems. AIAA paper 2001-1753; 2001. [27] Coblish JJ, Smith MS, Hand T, Candler GV, Nompelis I. Double-cone experiment and numerical analysis at AEDC hypervelocity wind tunnel no. 9. AIAA paper 2005-0902; 2005. [28] Babinsky H, Li Y, Ford CWP. Microramp control of supersonic oblique shockwave/boundary-layer interactions. AIAA J 2009;47:668–75. [29] Li Q, Liu C. LES for supersonic ramp control flow using MVG at M = 2.5 and Reh = 1440. AIAA Paper 2010-0592; 2010. [30] Holden H, Babinsky H. Effect of microvortex generators on separated normal shock/boundary layer interactions. J Aircr 2007;44:170–4. [31] Anderson BH, Tinapple J, Surber L. Optimal control of shock wave turbulent boundary layer interactions using micro-array actuation. AIAA Paper 20063197; 2006. [32] Mohd R, Hossein Z, Azam C, Konstantinos K. Micro-Ramps for hypersonic flow control. Micromachines 2012;3:364–78. [33] HDF5 user’s guide. Version 1.8.7; 2011.