 |
EPFL-SCR No 8
|
 |
| |
EPFL-SCR No 8
|
|
Mark L. Sawley, Olivier Byrde and David Cobut, EPFL - Laboratoire de Mécanique des Fluides
| Surface pressure and streamribbons computed for the flow around a solar racing car designed by the Biel School of Engineering for the across-Australia World Solar Challenge |
La simulation numérique joue un rôle toujours plus important pour la résolution des problèmes d'écoulements incompressibles. Cette étude montre comment la combinaison d'un solveur de type Navier-Stokes et d'un ordinateur massivement parallèle permet de simuler des écoulements incompressibles dans des situations tri-dimensionnelles complexes. L'écoulement d'un fluide fortement visqueux dans un mélangeur statique, ainsi que l'écoulement de l'air autour d'une voiture solaire, tous deux calculés sur un Cray T3D, illustrent l'éventail des possibilités d'un tel système.
Computational Fluid Dynamics plays an increasing role in the design process for a wide variety of applications involving incompressible flows. The present study is concerned with the use of a Navier-Stokes solver combined with the large computational capacities of a parallel computer system to perform complex three-dimensional incompressible flow simulations. Examples are presented of numerical simulations of the flow in an in-line static mixer and around a solar racing car, which have been undertaken on a Cray T3D system.
In addition to the traditional empirical and experimental techniques, numerical methods are currently employed - both to understand complex flow phenomena and in the design phase - for an increasing number of flow applications. For a wide range of application areas, the flow can be considered to be incompressible, necessitating the use of appropriate physical modelling and numerical methods to define and solve the governing set of partial differential equations. At the LMF, fundamental studies of turbulence and stability are being undertaken, in conjunction with the numerical simulation of complex flows in pumps, water turbines and chemical mixers, as well as for automotive aerodynamics. Of particular current importance is the use of high-performance parallel computer systems to perform complex flow simulations. This paper presents results obtained in studies from two of the above-mentioned incompressible flow application areas. These studies have been undertaken within the framework the Cray Research - EPFL Parallel Application Technology Program, using the 256-processor Cray T3D system installed at the EPFL.
The first application involves the simulation of three-dimensional laminar flow in an in-line static mixer. The tracking of particle paths through the mixer is shown to provide detailed insights into the basic flow phenomena, and thus clarify the mixing procedure. To obtain sufficient resolution of the complex flow behaviour, as well as qualitative and quantitative measures of the mixing properties via the use of a large number of particle trajectories, both the numerical simulation and the particle tracking have been undertaken on the Cray T3D.
The simulation of the three-dimensional incompressible flow around a solar racing car is also presented. This study has been undertaken to assess the aerodynamic behaviour of a prototype designed for the 1996 World Solar Challenge. Good agreement has been obtained between the computed values of the forces and moments on the car surface and those measured in windtunnel experiments.
Numerical simulations were performed using the multi-block flow solver sagarmatha [1], developed within the institute for the numerical simulation of steady / unsteady, inviscid / laminar / turbulent, incompressible flows. This code solves the Euler or Reynolds-averaged Navier-Stokes equations on 3D block-structured computational meshes. The numerical method employs a cell-centered finite volume discretization with an artificial compressibility method to couple the pressure and velocity fields. A high-order upwind spatial discretization scheme based on the approximate Riemann solver of Roe is employed for the advection terms, while the diffusion terms are discretized using a central approximation. The time integration of the unsteady Navier-Stokes equations is performed using an implicit two-stage Runge-Kutta scheme. At each time step, a non-linear system resembling the equations for stationary flow is solved using an ADI method.
The numerical simulation of complex flows requires significant computational resources, both in performance and memory. The present study has therefore employed a high-performance parallel computer system, the
256-processor Cray T3D system. Parallelization of the sagarmatha code is achieved by resolving the flow equations in all blocks (sub-domains) in parallel by assigning one block to each processor [2]. Communication between processors, necessary to exchange data at the edges of neighbouring blocks, is performed using message passing (PVM). A preliminary study [3] has shown that this programming model provides the desired combination of performance and portability. The communication overhead has been measured to account for less than 5% of the total resolution time. Load balancing is determined by the amount of computation and communication to be performed by each processor, and is optimized by choosing blocks with an approximately equal number of mesh cells.
Static mixers are in-line mixing devices which consist of mixing elements inserted in a length of pipe. Mixers of this type are used in continuous operation, with the energy for mixing being derived from the pressure loss incurred as the process fluids flow through the mixing elements. A variety of element designs are available from various manufacturers, with the number of elements required for a particular application being dependent on the difficulty of the mixing task (more elements are necessary for difficult tasks). The range of applications for static mixers is very wide, and includes both laminar and turbulent processes [4].
The optimization of chemical mixers is traditionally performed by trial and error, with much depending on previous experience and wide safety margins. Some of the inherent difficulties of this approach can be overcome by numerical flow simulation. The ability to simulate numerically the flow through a mixer can contribute significantly to the understanding of the mixing process and provide for better, faster and cheaper design optimization.
The geometry considered in the present study corresponds to the Kenics KM static mixer manufactured by Chemineer Inc. [5]. This mixer design has been employed in the process industry since the mid-1960s, mainly for the in-line blending of liquids under laminar flow conditions. The Kenics mixer is comprised of a series of mixing elements aligned at 90°, each element consisting of a short helix of one and a half pipe diameters in length. Each helix has a twist of 180° with right-hand and left-hand elements being arranged alternately in the pipe, as illustrated in Figure 1. It is generally proposed that flow division and azimuthal flow reversal, together with the induced radial motion, all contribute to the performance of the Kenics mixer.
| Figure 1 - Perspective view of a 4-element Kenics KM static mixer |
In the present study, numerical simulations have been performed for Kenics mixers having up to 12 elements, for a range of Reynolds number (which is a measure of the ratio of inertia forces to viscous forces), 
. The results presented in this paper will focus on flow in a 6-element mixer with Re = 25 and 100.
The symmetry of the Kenics mixer described above allows for the computation of the flow in only one half of the pipe. The 3D meshes for each of the mixing elements and inlet and outlet pipes are obtained by rotation and replication of a basic 2D mesh along the z axis [ [6]. The complete mesh for the whole mixer is constructed by juxtaposition of these component meshes. Multi-block decomposition is very natural for this mesh. The partition used for the present study involved the subdivision of each of the different pipes and mixing elements into eight blocks.
Three meshes have been considered in this study. The coarse mesh contains 9,216 cells per element, the medium mesh 47,040 and the fine mesh 188,160. For a 6-element mixer (i.e. 6 mixing elements plus inlet and outlet pipes), this corresponds to a total number of cells of 73,728, 376,320 and 1,505,280, respectively (for half the mixer cross-section). It is noted that the coarse mesh is of similar resolution to that currently used by industry. Since the number of blocks per element is 8, a total of 64 processors is used to compute the flow in a 6-element mixer. More processors could be used by further subdivision of blocks; due to the high level of code scalability [2] this will provide higher performance. All blocks have exactly the same number of cells to ensure almost-optimal load balancing between processors.
The computational time to obtain a converged flow solution was found to be dependent both on the mesh and on the Reynolds number, significantly more iterations being required for convergence of low Re flow solutions [6]. For the medium mesh and Re = 25, 4 hours were required on 64 processors of the Cray T3D, while for Re = 100, the flow solution was obtained in 2 hours.
To determine the efficiency of a chemical mixer, it is necessary to establish means by which the extent of the fluid mixing can be gauged both qualitatively and quantitatively. For the computation of a multi-species flow in a static mixer, it is possible to determine the extent to which an inhomogeneous species concentration at the inlet is distributed with the passage of the fluid through the mixer. Numerically, this approach has been found to be inappropriate due to the domination of numerical diffusion in the resolution of the species transport equations compared to the physical diffusion.
A more accurate determination of the mixing phenomena can be obtained by calculating the trajectories of fluid particles in the flowfield of the mixer. Such an approach, which can be applied to a single species fluid, has been adopted for the present study. To obtain a global evaluation of the mixing, it is necessary to study the trajectories of a large number of particles; the particle trajectories were therefore calculated in parallel on the Cray T3D.
Two different approaches were envisaged for the parallel computation of the trajectories. The first involves maintaining the block distribution of the flowfields across the processors and calculating the section of each trajectory by the processor whose memory contains the corresponding flowfield section. Alternatively, by replicating the velocity flowfield on each processor, complete trajectories for different particles can be calculated by each processor. The second approach has been adopted here, since it enables a large number of particle trajectories to be computed in an "embarrassingly parallel" manner, with the replication of the velocity flowfield requiring generally less than 1% of the post-processing computation time. Using this approach up to 1 million particle trajectories could be computed in a reasonable time on the Cray T3D. Table 1 indicates representative times for the calculation of the particle trajectories (neglecting the additional time required to undertake statistical analysis and for I/O). Note that due to symmetry, distributing the particles initially in only half of the mixer allowed the trajectories of twice as many particles to be determined.
| Number of particles | Cray Y-MP time | Cray T3D time |
| 4,160 | 30 min. | 8 min. (16 proc.) |
| 16,012 | 120 min. | 15 min. (32 proc.) |
| 65,792 | -- | 30 min. (64 proc.) |
| 262,656 | -- | 60 min. (128 proc.) |
| 1,049,600 | -- | 120 min. (256 proc.) |
Table 1 - Computational times required for different numbers of particle trajectories
Computations have been performed on the Cray T3D for each of the configurations and meshes mentioned above. Results presented here were obtained for a 6-element mixer using the medium mesh, these results being almost identical to those computed using the fine mesh. On the contrary, differences observed with the results obtained using the coarse mesh indicate that this mesh may not provide sufficient resolution of the flow within the mixer.
In order to obtain more information about the fluid motion in the mixer, the flowfields at different axial positions have been studied. At the inlet of the mixer, uniform flow is imposed; an approximately parabolic profile of the streamwise velocity develops by the beginning of the first mixing element. Further downstream, it is observed that in each half-domain section, the fluid rotates in the opposite sense to the helical twist of the mixing element, leading to the creation of a vortical structure. For certain Reynolds numbers (e.g. Re = 100), a secondary vortex develops on the suction side of each element, which aids in the mixing procedure [6].
 |
Figure 2 - Computed particle locations at different axial positions along the mixer (from top left to bottom right: at inlet, near the end of each of the six mixing elements, and at the outlet) for flow with Re = 25 |
Further insight into the mixing process can be gained by tracing the passage of different fluid particles through the mixer. Figures 2 and 3 show, for flows with Re = 25 and 100, respectively, the distribution of two initially semi-circular disks of particles (distinguished by colours red and blue) at different distances along the mixer. In these plots, the position of 65,792 particles of each colour are shown at each axial location. These plots illustrate the randomization of the particle positions via the combined effect of flow division and reversal.
|
Figure 3 - Computed particle locations at different axial positions along the mixer (from top left to bottom right: at inlet, near the end of each of the six mixing elements, and at the outlet) for flow with Re = 100 |
Comparing Figures 2 and 3, it can be qualitatively observed that the mixing is more efficient for the Re = 100 flow. For this flow case, each mixing element divides by two the size of the structures observed, with the fluid wrapping around the main and secondary vortices. The presence of flow division between elements both increases the effective interface length and enhances the action of flow reversal. For Re = 25, however, a number of relatively large regions containing particles of only one colour are present at the mixer outlet.
In order to quantify the mixing process, we define the global structure radius as the value for which 99% of the particles have at a closer distance a particle of a different colour. In Figure 4 is plotted, as a function of the number of particles trajectories, the global structure radius calculated for Re = 25 and 100. Also shown for comparison is the global structure radius calculated for a random distribution of particles. It can be observed that for the random distribution of particles, as expected, the global structure radius decreases with the number of particle considered. However, for each of the flow simulations, a limiting value is reached. Figure 4 indicates that a large number of particle trajectories must be considered in order to obtain the limiting value of this quantitative measure of mixing, particularly for efficient mixing for which only small-scale flow structures persist.
|
Figure 4 - Calculated global structure radius at the mixer outlet as a function of the number of particle trajectories analyzed |
The motion of the particles in the Kenics mixer presented above is characteristic of chaotic behaviour. The computed exponential reduction of the size of the flow structures with the number of mixer elements is consistent with the presence of chaotic advection (see, e.g. [7], [8]). At the small spatial scales, molecular diffusion (which has been neglected in the present analysis) will cause randomization of the particles, thus completing the mixing process. Work is presently being undertaken to quantify the mixing procedure in terms of chaos measures, such as Lyaponov exponents. This should enhance a detailed study of the efficiency of the Kenics mixer as a function of physical and operational parameters, which will be beneficial in the design optimization procedure.
The World Solar Challenge (WSC) is the world's premier solar racing car event, traversing 3000 km across Australia from Darwin to Adelaide [9], [10]]. This race has been held every three years since its inception in 1987, and attracts increasing attention, particularly from automotive manufacturers. The 1996 WSC includes over 50 entries from 14 different countries. The WSC provides a showcase for solar science and a challenge in creative vehicle design, energy efficiency and innovation in electrical power systems.
Substantial progress has been achieved in solar racing car performance over the past decade. During the 1993 WSC, the winning entry from Honda Motor Corporation completed the race with an average speed of approximately 85 km/h, attaining top speeds in excess of 100 km/h. To achieve and improve on such performances, external aerodynamics is of paramount importance, an optimal design providing minimum drag and maximum car stability. To date, empirical techniques together with limited windtunnel testing have generally been used for the design of the shape of solar racing cars. Many different solar racing car shapes have been used in the past (see references in [9]), with a corresponding wide range of aerodynamic properties. This variation indicates that more research and development is required in the path towards optimal design and hence improved performance.
While possessing various specific features, solar racing cars also exhibit many of the features typical of automotive aerodynamics. These include complex 3D flowfields with turbulent boundary layers and wakes, flow separation and ground effects. Numerical simulation of such flows present a number of challenges, since the basic physical phenomena and their modelling are not always well understood. In addition, complete numerical simulations are generally computationally intensive, requiring substantial computer power and memory.
The present study investigates the use of numerical simulation on the Cray T3D to aid in the design optimization of a solar racing car. Comparison of numerical predictions with measured values of the aerodynamic forces and moments are undertaken for a particular car design.
Within the framework of the present study, different solar racing car designs have been examined both experimentally and numerically. This series of designs has been provided by the Biel School of Engineering, who have gained 3rd, 1st and 2nd places, respectively, in the first three World Solar Challenges [11]. The present study was motivated by the aerodynamic design optimization for the 1996 WSC. Figure 5 shows the 1993 WSC entry "Spirit of Biel-Bienne III" [11], and one of the prototypes for the 1996 WSC (denoted here as WSC96-P2) for which numerical simulations have been undertaken.
Each of the racing car designs are built around the solar collection panel, which consists of an approximately flat surface, the area of which is limited to 8 m2 by the World Solar Challenge race regulations. The vehicles considered have three wheels, however the 1993 design has two front wheels, while the WSC96-P2 prototype has two rear wheels. All wheels are enclosed in casing to minimize drag. Contrary to earlier designs, the WSC96-P2 prototype has the cockpit in a forward position above the front wheel.
 |
Figure 5 - The "Spirit of Biel-Bienne III" [11] (top) and windtunnel model of the WSC96-P2 prototype |
 |
Windtunnel tests were undertaken for each of the solar racing car designs using 1:8 scale models. Experiments were conducted in a closed return windtunnel at the LMF having an open test section of 50x70 cm, and a flow speed of 18.75 m/s. The model was positioned to be parallel to the surface of a stationary flat ground plate, with the lower surface of the wheels located just above the plate surface. Of interest in the present study are the aerodynamic forces and moments on the car (as defined in Figure 6). Measurements were made for a range of angles of sideslip, corresponding to that encountered due to lateral winds in the 1993 WSC [12]. The windtunnel results indicate that changes in the basic design shape of the solar racing car results in significant modifications to its aerodynamic properties, with the desired decrease in drag being obtained at the cost of decreased aerodynamic stability.
 |
Figure 6 - Definition of the aerodynamic force and moment components (top) and sideslip angle (bottom). |
 |
Numerical flow simulations have been undertaken corresponding to the windtunnel conditions for flow around the WSC96-P2 prototype, enabling validation of the computed flow results. As a first step in the numerical simulation procedure, inviscid flow calculations were performed. Such simulations alleviate the complex mesh generation problem, avoid the uncertainty of turbulence modelling, and reduce computational costs.
To construct a suitable computational mesh, the surface of the windtunnel model of the WSC96-P2 prototype was first digitized using a mechanical 3D 3-axis scanner. The data points obtained were treated using a CAD system to construct a structured surface mesh. The main difficulty with this geometry is the coupling of the high curvature areas with the remaining essentially flat regions. Special attention was therefore paid to the meshing of the leading edge, cockpit and wheel casing surfaces.
Due to the lateral symmetry of the model, a computational mesh was generated for only half the flow domain; this mesh being reflected in the symmetry plane to perform simulations of flow with a sideslip angle. For the simulation of flow with zero sideslip angle, the assumption of flowfield symmetry enabled the computation to be performed using only half the flow domain.
To facilitate the construction of the volume mesh, the computational domain was divided into seven blocks. Classical transfinite interpolation was used to generate the mesh within each of the blocks, with mesh cells clustered close to the car surface. The resulting multi-block computational mesh was comprised of a total of 362,280 mesh points (for half the flow domain) [12].
In order to improve the load balance on the parallel computer system, the blocks containing the largest number of mesh cells were subdivided, resulting in a total of 16 approximately equal size blocks, as shown in Figure 7.
 |
Figure 7 - The number of mesh points per block (top), and the computation and synchronization / communication times measured for each of the blocks (bottom) |
 |
The 3D computation of inviscid flow with zero sideslip angle using the 16-block mesh required 1.5 hours on 16 processors of the Cray T3D to obtain a reduction of the residual of four orders of magnitude. The measured computation times for the individual blocks presented in Figure 7 show that adequate load balancing could be achieved by block subdivision. Nevertheless, the synchronization overhead could be further reduced by the use of a greater number of processors. It should be noted, however, that on the Cray T3D the number of processors allocated must be a power of two.
Similar performance results were obtained for flow with non-zero sideslip angle using a 32-block mesh (constructed by reflection) and 32 processors of the Cray T3D. However, due to the different number of iterations required for convergence, the resolution time differed from just under one hour for 5O sideslip to almost two hours for 10O sideslip.
An initial exploratory study concentrated on the general inviscid features of the flow with zero sideslip angle around the WSC96-P2 model assuming windtunnel conditions. Figure 8 shows the computed pressure on the surface of the car, and the velocity streamlines calculated just above the upper surface. A strong pressure peak is observed on the leading edge surface, with depressions in the region around the cockpit and wheel casings. The pressure is seen to be approximately constant over the majority of the upper surface. However, the presence of the wheels results in a greater variation of the pressure on the lower surface, with the presence of a region of low pressure between the rear wheels.
 |
Figure 8 - Computed surface pressure contours and streamlines. |
An examination of the computed streamlines presented in Figure 8 indicates that along the lateral edges the streamlines are deflected slightly from the upper to the lower surface, resulting in the formation of weak trailing edge vortices. The possible presence of flow separation near the trailing edge of the cockpit and subsequent reattachment are not captured by inviscid flow simulations. However, preliminary 2D flow studies [13] have predicted that flow separation occurs. Assuming the flow to be laminar results in a large flow recirculation zone, leading to unsteady flow with the convection of vortices along the upper surface of the car. If the flow is assumed to be turbulent (as is the case in the windtunnel), separation is delayed and reattachment occurs downstream.
Numerical simulations have also been performed for sideslip angles of 5°, 10° and 15°. A comparison of the forces and moments computed from the surface flowfields with those measured in the windtunnel is presented in Figure 9. It should be emphasized that the present numerical study has considered only inviscid flow simulations. In the assumed absence of global regions of flow separation, these simulations should therefore yield a good approximation of the pressure contribution to the forces and moments. However, the contribution due to skin friction is not simulated numerically.
 |
Figure 9 - Comparison of numerical simulations (x)
and windtunnel measurements ( )
for the three components of the forces (left) and moments (right) as a
function of sideslip angle.. |
 |
|
 |
|
 |
|
 |
|
 |
The results presented in Figure 9 show that generally good agreement is obtained for the dependence of the force and moment components on the sideslip angle. Such agreement suggests that the pressure contribution dominates the forces and moments exerted on the car. Nevertheless, some quantitative discrepancies can be observed. These can be explained by the following considerations:
The present study has shown that modifying the car shape can result in significantly different aerodynamic properties. Further optimization of the shape of the solar racing car should take advantage of numerical techniques to supplement empirical and experimental methods. Future numerical studies, concentrating on viscous flow simulations and enhanced physical modelling (e.g. by the use of appropriate advanced turbulence models), will provide a greater degree of applicability of the numerical simulations. In particular, the simulation of "real world" influences on the full-scale car - such as tyre friction on the road surface, relative motion of the car and road, rotating wheels, and change of the centre of gravity due to the motor, pilot and equipment - should present a valuable contribution to the design process.
The present study has shown that the numerical simulation of three-dimensional incompressible flows can necessitate a considerable computational effort, not only for the calculation of the flowfields, but also - as is the case for the chemical mixer - in the post-processing phase. Such simulations, which are often not possible or practical on conventional vector supercomputer systems or workstations, benefit greatly from the increased capability of a high-performance parallel computer system, such as the Cray T3D. The possibility of increased spatial and temporal resolution, more adequate physical and numerical modelling, enhanced post-processing and reduced simulation time and costs will all contribute to the integration of numerical flow simulation into the design process in a number of industrial application areas.
Discussions with R. LaRoche (Cray Research Inc.) and A. Bakker (Chemineer Inc.) concerning the flow simulations in the chemical mixer are gratefully acknowledged. The study of the solar racing car was motivated by an interaction with the Biel School of Engineering; G. Insom is particularly thanked for providing the geometry and windtunnel models as well as a large measure of enthusiasm. P. Barbey and V. Colmar are acknowledged for their contributions via their undergraduate projects. J.-D. Reymond generated the computational meshes, and Th. Ursenbacher and T.-V. Truong provided the experimental data and numerous comments on their interpretation. Financial support for this study was partially provided by the Cray Research - EPFL Parallel Application Technology Program.
refer to contents |
©EPFL-SCR # 8 - 1996 |
| your comments |