Transient Simulation Research on Automobile Aerodynamic Lift Based on LBM Method.

Yong, Zhang ; Zhengqi, Gu ; Shuichang, Liu 等

1. Introduction

In order to save fuel, reduction of aerodynamic drag has always been concerned by car-body design engineer, but little attention was paid to aerodynamic lift. However, in terms of automobile aerodynamic lift being proportional to the second-order of speed, the reduction of tire gripping capacity could be caused by the increasing speed, which would bring the vehicle handling and stability problems. The aerodynamic lift of transient changing characteristics is considered as the inducement for all kinds of traffic accidents under high speed or speeding in an instant, even known as the life-line of high-speed automobiles [1]. However, for many years the attention of automobile aerodynamic lift is devoted largely to steady state, while the transient characteristics research is rarely presented.

Based on the theory of the aircraft aerodynamic lift, the steady-state theory of automobile aerodynamic lift suggests that the air flowing velocity through the upper is faster than the bottom because the upper surface of car-body is longer than the bottom surface. According to the Bernoulli equation, the pressure on the upper surface is smaller than that on the bottom surface. It is believed that the pressure difference between the upper and the bottom surface brings about the aerodynamic lift [1]. But, in reality, carbody is a three-dimension geometry entity with protrusive objects and air flow surrounding the car is easy to be separated. Due to the phenomena of the ground effects, the aerodynamic lift of some car-body models is negative, which implies the 2D path steady-state theory is insufficient to explain the real situation.

Wind tunnel investigations on car-body aerodynamic lift have been performed by Joel A. Walte [2], Gerhard Wickern [3], B. Hetherington [4], Christoffer Land-strom [5], Ding Ning, et al [6]. These studies show that the car lift is very sensitive to the bottom airflow state, as well as the boundary layer suction rate and the model installation misalignment. In order to study the influence of the moving floor on aerodynamic force in wind tunnel investigations, S. Krajnovic [7] calculated the aerodynamic lift of the Ahmed body model for slant angle of 25 degrees using finite volume method. The calculation error was around 14% between the simulated results and the wind tunnel test data. Yoshihiro Okada et al. [8] calculated aerodynamic lift using the finite volume method to get unstructured tetrahedron mesh of the model, whose results suggested that the conventional aerodynamic evaluation method focusing on steady aerodynamic lift coefficient is insufficient to evaluate automobile straight-ahead stability at high speed. A new aerodynamic evaluation method for automobile stability was proposed to indicate the rear of the car playing the most important role in aerodynamic lift.

However, these studies mentioned above were restricted from the disadvantages of the steady state simulations in finite volume method. In general, the meshing depends on experience to unavoidably bring about the abnormal grid on the complicated surface of vehicle body, especially the fixed and equal thickness boundary layer mesh. It could be clearly seen in Fig.1, the mesh model was divided in ANSYS14.0-ICEM.Due to aerodynamic lift being very sensitive to car body boundary layer, it is thought that the precision of the simulation results for aerodynamic lift would be limited by the reasons mentioned above, while the simulation results of aerodynamic drag usually could meet the request of practice application, testified by the most of cases.

The air flow field around the car body is always changing and the transient characteristics of automobile aerodynamic lift at high speeds are rarely reported in the literature. Because of those, it is necessary to study the transient aerodynamic lift without referring to the mesh approach. Different from the finite volume method, the Lattice Boltzmann Method (LBM) divides the macroscopic fluid into a series of fluid particles, which would act as collision and migration in the discrete grid model node [9]. From the macroscopic, the particles are controlled by Navier-Stokes (N-S) equation using statistical method [10]. It is worth mentioning that making discrete space in LBM does not need to mesh the model, avoiding the calculated error caused by the distortion gird. Instead of that, the particle is considered as interpolation node to improve the calculation accuracy by controlling its density to realistically simulate the boundary layer.

LBM was originally developed from the Lattice Gas Automata (LGA) [11]. Recent research led to major improvements on this approach [12-15]. Some studies using commercial software such as Power Flow, XFLOW, etc., have indicated that the LBM method could offer a better description of situations and fluid flow mechanism [16, 17].

This paper was devoted to study transient automobile aerodynamic lift based on the LBM method. The D3Q27 model is implemented to obtain the discrete particles for international standard model. The adaptive wake refinement feature is employed to obtain a series of regular fluid particles for macroscopic fluid and the large eddy simulation (LES) combined with wale-viscosity model is carried on studying the transient characteristics of automobile aerodynamic lift. The simulate results are compared with the experimental results obtained from the HD-2 wind tunnel in Hunan University. The temporal and spatial characteristics of rear velocity field and vorticity field are contrasted to the PIV testing data as well.

2. Numerical methods

2.1. Lattice Boltzmann method

The standard LBM equation is defined as following [18]:

[mathematical expression not reproducible] (1)

where f is the distribution function, Q is collision operator, r is the position vector of one particle, [partial derivative] is the particle velocity along a certain direction.

Q is usually expressed by the BGK model, which can be written as following:

[mathematical expression not reproducible] (2)

here [f.sub.eq] is the local equilibrium distribution function, [tau] is the relaxation time spent on the particles collision to achieve equilibrium. It has been proved that, if the discrete speed set is symmetrical, equation (1) and the N-S equations can be unified by equations (3-4) like this [17, 18].

Density:

[mathematical expression not reproducible] (3)

Velocity:

[mathematical expression not reproducible] (4)

Internal energy:

[mathematical expression not reproducible] (5)

Here [C.sub.s] is the speed of sound. Because the space direction can be statistical, the D3Q27 discrete grid model, as the highest precision at present, is used for spatial discretization. When e is equal to 26, the model is shown in Fig. 2. In Eqs. (3), (4) and (5), it could be greater than or equal 0 and less than or equal 26.

Finally, the macroscopic viscosity is written as a function of the relaxation parameter as:

[mathematical expression not reproducible] (6)

For a positive viscosity, the relaxation time must be greater than 0.5 [18, 19].

2.2. Turbulence model

Because of the LBM constructing to describe the transient motion and clashing of the microscopic particle, the turbulence model of the LES fits well with this purpose.

In LES, large scale eddies can be solved directly, and sub-grid scale model was applied to solve small scale eddies. Various studies have shown that LES is an effective method. Taking consideration of local eddy-viscosity and near wall behaviours, the wall-adapting local eddy viscosity model is used.

The turbulent eddy viscosity is expressed as:

[mathematical expression not reproducible] (7)

[mathematical expression not reproducible] (8)

[mathematical expression not reproducible] (9)

[mathematical expression not reproducible] (10)

where k is Karman constant, d is distance from the wall and the constant [C.sub.[omega]] is 0.325. Uniform non-equilibrium wall function is used to simulate the boundary layer [20].

3. Simulation and conditions

3.1. The geometry model

In order to present a representative case, the inclined back MIRA model is used. MIRA model is the international standard model widely used in fundamental research of automotive aerodynamics. The length (L), width (W), and height (H) of the model body is 1389 mm, 542 mm and 474 mm. The real model as shown in Fig. 3 is located in the Hunan University HD-2 wind tunnel.

3.2. Spatial discretization

D3Q27 discrete grid model is used for spatial discretization based on the commercial software XFLOW 2013.Build9.0. In order to get a sound, result within a reasonable time, the spacing size of initial particle around far field is taken as 0.01 m, while the resolution near the body wall and trailing vortex is 0.003125 m. The lattice structure can be modified later by the solver in order to adapt to the flow patterns. The adaptive wake refinement feature is based on the module of the vortices field. To eliminate the effects of wall interference on the calculation process, the virtual wind tunnel should be ten times length, seven times width and five times height of the MIRA model. The distance from the entrance to the car is three times length. Finally, the length, width and height of the calculation domain size is 13890 mm, 3794 mm and 2370 mm, and the blocking ratio is 2.86%. The number of the initial particle element is 9 million. The initial discrete particle distribution is shown in Fig. 4.

3.3. Boundary condition

The boundary of the computational domain conditions is shown in Table 1.

4. Results and discussion

4.1. Drag and lift coefficient

The computer configurations include 64-bit Windows 7 system, Intel Xeon (RE) with 8 cores, dual CPU (2.80 GHz), and RAM of 48.0 GB. According to the computational domain length and inlet velocity boundary condition, the time of transient flow field analysis is 0.4 seconds, and time step [delta]t is taken as 0.001s. After taking around 149 h to solve this model in XFLOW 2013.Build9.0, the number of particle amounts to 50.82 million. The final discrete particle distribution is shown in Fig. 5. The curve of drag and lift coefficient of the MIRA model at all step times are shown in Figs. 6 and 7.

As shown in Fig. 6, the aerodynamic drag coefficient [c.sub.d] decreases rapidly from the maximum during the initial stage. After a while, the oscillation amplitude in a larger scope appears among the next 0.2 s. Then, oscillation starts to disappear and the aerodynamic drag coefficient [c.sub.d] quickly becomes to be stabilized with the amplitude staying from 0.278 to 0.286. The average cd is 0.283 between 0.25 s and 0.4 s, which is comparable to the experimental result 0.277 obtained from the HD-2wind tunnel of Hunan University. In addition, the IVK automotive model wind tunnel of Stuttgart University data is 0.278 [21], and the TJ-3 Boundary Layer Wind Tunnel of Tong Ji University is 0.274 [22]. The deviation of the simulated result based on the LBM between them is 0.006, 0.005 and 0.009 respectively. It is shown that this method is satisfactory.

However, the aerodynamic lift coefficient [c.sub.l] with time shown in Fig. 7 is obviously different from the aerodynamic drag coefficient [c.sub.d] time history. At the initial stage, [c.sub.l] increases suddenly from the least negative value to the maxi-mum positive value. Then, the irregular fluctuation appears. Starting from 0.25 s, the aerodynamic lift coefficient [c.sub.l] becomes stable and its amplitude stays in the range of -0.003 and 0.075. The average value [c.sub.l] is 0.057 between 0.25 s and 0.4 s. Due to the absence of aerodynamic lift experiment data for a unified standard car model, this simulated result is compared with the test data obtained from HD-2 wind tunnel. The deviation of the two results is around 0.009 when the tested result of the aerodynamic lift coefficient [c.sub.l] is 0.048. This error may be attributed to the measuring inaccuracy of aerodynamic lift in wind tunnel or data processing method. It could be believed that the error 0.009 is smaller, and acceptable.

4.2. Velocity field and vorticity field

In order to further explain the phenomenon of non-periodical oscillation for the aerodynamic lift coefficient [c.sub.l], the transient temporal and spatial characteristics of the auto body turbulent vorticity and rear velocity vector were studied. The car-body turbulence vorticity at 0.365 s is shown in the Fig. 8, and seven vertical sections are also intercepted behind the model. The section A is located at 1/4 width of the car, and the section B is at the end wall of the car, while the other section C, D, E, F, G are located parallel at 1/4, 1/2, 1, 5/4/, 9/4 times length of the car, respectively. Six sections instantaneous velocity vector diagrams at 0.365 s are showing in Fig. 9, and six transients (respectively at 0.340 s, 0.345 s, 0.350 s, 0.355 s, 0.360 s, 0.365 s) velocity vector diagrams of section A are shown in Fig. 10.

At the same time, to verify the reliability of the method presented in this paper, the particle image velocity (PIV) experiment about trailing vortex was implemented in the HD-2 wind tunnel, while the particles analysis system was MicroVecV3.3.2. ND-YAG double pulsed laser. Its single pulse energy is 500 MJ and 532nm, and it is green light, and the maximum working frequency is 10 Hz. The CCD digital camera with pixels of spatial resolution 4000 and 2672, and the collection rate 5 frames per second and the microPulse725 type synchronous controller with delay time precision 10ns are both shown in Fig. 11. System connection schematic diagram is shown in Fig. 12 and the HD-2 wind tunnel is shown in Fig. 13. The six continuous instantaneous automobile rear velocity vector images for the section A based on particle image velocimetry are shown in Fig. 14.

Combining the LBM method with the adaptive wake refinement feature and LES with wale-viscosity model, a detailed and explicit simulation of flow field around the model automobile is shown in Fig. 8.

According to the results shown in Fig. 9, the car rear flow field still shows a pair of inward twisting vortex system, such as two black circles (I and II) in each figures. But the boundary of the pair of vortex system becomes less sharp, irregular and asymmetric. Different scales of eddies crush and drag with each other appear in different positions. Judging by the evolving of motion vortex, six different positions are shown in Fig. 10, which has the same characteristics. The phenomenon of squeezing and dragging becomes more pronounced. At the time of 0.340 s, there is a large clockwise vortex (A) at the rear upper position of the luggage chamber, but with a counterclockwise vortex (B) in rear lower position. When the time passing away 0.005 s, the interactions between the vortex A and B in the original location result in the disorder air flow, and the originally large clockwise vortex A forms a new small vortex (A2) after crushing, which is then being pushed to the rear. At the next time, two new vortex cores appear in the rear. There are also two vortex patches with opposite revolving direction at the time of 0.355 s, but the position is further backwards. In the following 0.005 second, two new vortex cores (A and B) with opposite revolving direction also appear near the rear of the luggage chamber, but a large clockwise vortex (A2) is set on behind A. Then, at the last moment, because of crushing and dragging, the air flow is disorder again at the opposite of A2, and a mature turbulence vortex appears near the luggage chamber. In conclusion, the migration of the vorticity is clearly revealed in the six instantaneous velocity vector diagrams of section as shown in Fig. 10. The generation, migration, development and vanish of the vortex core evolves to be non-periodical but random shedding and oscillation.

Furthermore, the six continuous instantaneous images are collected once every 0.2 seconds based on PIV, which are shown in Fig. 14. There is a large clockwise vortex (A) in the rear upper position of the luggage chamber, and a small eddy (B) in rear lower position in the first picture. In the second picture, flow filed of A and B is extremely turbulent. In the third picture, two new vortex cores (A and B) also appear in the rear. In the fourth picture, it is confusing again, but after crushing, and it is being pushed to the rear, with a small vortex (B2) setting on behind B. In the fifth picture, it's similar with the first picture extremely. In the sixth image, air flow is disorder again.

Comparing the simulated result shown in Fig. 10 with the tested results shown in Fig. 14, this two results could not ensure synchronization in space and time, because the PIV is limited by the frame rate of CCD and LBM method is based on particle motion belonging to the Lagrange method. However, both of them show that the transient motion of the rear vortex is non-periodical turbulence. There is a high similarity in the description of the vortex transient characteristics between the simulated results and the experiment results obtained from the PIV wind tunnel test around the rear of the car. Both of them show that the generation, migration, development and vanish of the vortex core evolves to be non-periodical but random shedding and oscillation. This phenomenon is closer to real situation. It is concluded that this simulation method is applicable.

Combining the simulated and experimental results, it could be seen that the flow around car tail is a highly separating 3D turbulent flow. There is strong non-periodical pressure pulsation at different position. According to the existing research conclusions of automobile aerodynamic lift [22], it is shown that the automobile tail is a major area of turbulent kinetic energy dissipation. Because of the tail flow field of the car being significant effect on aerodynamic lift, it is concluded that the asymmetric, unstable oscillations vortex results in the unstable aerodynamic lift. This is consistent with the non-periodical change of the aerodynamic lift coefficient as discussed above.

5. Conclusions

According to the data analysis and discussion, following conclusion are made.

1. The LBM method combining with adaptive wake refinement feature and LES with wale-viscosity model brings about more elaborate information of flow field around automobiles. The simulation results show agreement with the experimental data from wind tunnel test.

2. Transient simulation using the method presented in this paper showed significant difference from the steady calculation based on the classic finite volume method. The rear flow field of the car has a pair of irregular and asymmetric inward twisting vortex system, which is close to real flow station. It is clearly shown that the instability of aerodynamic lift is caused by the asymmetric and unstable oscillations rear vortex.

Acknowledgments

This work is financially supported by the Natural Science Foundation of Hunan Province in China (2017JJ2074), the Innovative team of the National Financial Support of China (0420036017) and the Science Research Project of Hunan province ministry of education (16B072). All authors thank Dr. Max LAI and Ivan.Fan from Beijing SOYOTEC Technology co., LTD.

References

[1.] Hucho, W.H. 1998. Aerodynamics of Road Vehicles, Society of Automotive Engineers, Warrendale, printed in the United States America: 35.

[2.] Walter, J.A.; Canacci, V.; Rout, R.K.; Koester, W.; Williams, Jack; Walter, T.M.; Stephen, J.; Nagle, P.A. 2005. Uncertainty Analysis of Aerodynamic Coefficients in an Automotive Wind Tunnel, SAE: 01-0870. http://dx.doi.org/10.4271/2005-01-0870.

[3.] Wickern, G.; Wagner, A.; Zoerner, C. 2005. Induced Drag of Ground Vehicles and Its Interaction with Ground Simulation, SAE: 01-0872. http://dx.doi.org/ 10.4271/2005-01-0872.

[4.] Hetherington, B.; Sims-Williams, D.B. 2006. Support Strut Interference Effects on Passenger and Racing Car Wind Tunnel Models, SAE: 01-0777. http://dx.doi.org/10.4271/2006-01-0565.

[5.] Landstrom, C.; Walker, T.; Lofdahl, L. 2010. Effects of Ground Simulation on the Aerodynamic Coefficients of a Production Car in Yaw Conditions, SAE: 01-0755. http://dx.doi.org/10.4271/2010-01-0755.

[6.] Ding, N.; Yang, Z.; Li, Q. 2013. Numerical Simulation on the Effects of Moving Belt System in Wind Tunnel on Aerodynamic Lift Force, Automotive Engineering 35(2): 143-146. http://dx.doi.org/10.3969/j.issn.1000-680X.2013.02.009

[7.] Krajnovic, S.; Davidson, L. 2005. Influence of Floor Motions in Wind Tunnels on the Aerodynamics of Road Vehicles, Journal of Wind Engineering and Industrial Aerodynamics 93(9): 677-696. http://dx.doi.org/10.1016/j.jweia.2005.05.002.

[8.] Akiyama, Y.; Nouzawa, T.; Nakamura, T.; Okamoto, S.; Okada, Y. 2009. Flow Structures above the Trunk Deck of Sedan-Type Vehicles and Their Influence on High-Speed Vehicle Stability 1st Report: On-Road and Wind-Tunnel Studies on Unsteady Flow Characteristics that Stabilize Vehicle Behavior, SAE: 01-004. http://dx.doi.org/10.4271/2009-01-0004.

[9.] Chen, H.; Chen, S.; Matthaeus, W.H. 1992. Recovery of the Navier-Stokes Equations Using a Lattice-Gas Boltzmann Method, Phys Rev A 45(8): 5339-5342. http://dx.doi.org/10.2307/3185570.

[10.] Qian, Y.H.; Orszag, S.A. 1993. Lattice BGK model for the Navier-Stokes equation: Nonlinear deviation in compressible regimes, E. p. Letters 21(3): 255-259. http://dx.doi.org/10.1209/0295-5075/21/3/001.

[11.] Chen, S.; Doolen, G.D. 2003. Lattice Boltzmann Method for Fluid Flows, Annual Review of Fluid Mechanics 30(5): 329-364. http://dx.doi.org/10.1146/annurev-fluid-121108-145519.

[12.] Aidun, C.K.; Clausen, J.R. 2009. Lattice-Boltzmann Method for Complex Flows, Annual Review of Fluid Mechanics 42 (1): 439-472. http://dx.doi.org/10.1146/annurev-fluid-121108-145519.

[13.] Shan, X.; Chen, H. 2007. A General Multiple-relaxation-time Boltzmann Collision Model, International Journal of Modern Physics C 18(4): 635-643. http://dx.doi.org/10.1142/S0129183107010887.

[14.] Franck, G.; Nigro, N.; Storti, M.; Elia, J.D. 2009. Numerical Simulation of the Flow around the Ahmed Vehicle Model, Latin American Applied Research 39(4): 295-306. https://www.researchgate.net/publication/228610009

[15.] Chen, H.D.; Kandasamy, S.; Orszag, S.; Orszag, R.; Shock, S.; Yakhot, V. 2003. Extended Boltzmann Kinetic Equation for Turbulent Flows, Science 301(5633): 633-636. http://dx.doi.org/10.1126/science.1085048.

[16.] Exa Corporation. 2014. Power Flow User's Guide Version 4.0.

[17.] Next Limit Dynamics SL Corporation. 2013. XFlow User's Guide Version.

[18.] He, Y.L.; Wang, Y.; Li, Q. 2009. The Theory and Application of Lattice Boltzmann Method, Chain Science Press, Beijing.78 p. (in Chinese).

[19.] Premnath, K.N.; Banerjee, S. 2012. On the Three-dimensional Central Moment Lattice Boltzmann Method, J. Statistical Physics 143(4): 747-794. http://dx.doi.org/10.1007/s10955-011-0208-9.

[20.] Sewoong, J.; Phares, J.D.; Srinivasa, A.R. 2013. A Model for Tracking Inertial Particles in a Lattice Boltzmann Turbulent Flow Simulation, Int. J. Multiphase Flow 49(3): 1-7. http://dx.doi.org/10.1016/j.ijmulti-phaseflow.2012.09.001.

[21.] Hoffman, J.; Martindale, B.; Arnette, S.; Williams, J. 2003. Development of Lift and Drag Corrections for Open Jet Wind Tunnel Tests for an Extended Range of Vehicle Shapes, SAE: 01-0934. http://dx.doi.org/10.4271/2003-01-0934.

[22.] Pang, J.A.; Li, Z.X.; Yu. Z.P. 2002. Correction Methods for Automotive Model Tests in TJ-2 Wind Tunnel, Automotive Engineering 24(5): 371-375. http://dx.doi.org/10.3321/2002.05.001.

Zhang YONG (*,**), Gu ZHENGQI (*,**), Liu SHUICHANG (***)

(*) School of Mechanical Engineering, Hunan University of Technology, Zhu zhou, 412007, China, E-mail: zhangyong7051678@163. com

(**) State Key Laboratory of Advanced Design and Manufacture for Vehicle Body (Hunan University), Chang Sha, 410082, China

(***) School of Mechanical Engineering, Hunan University of Technology, No. 401, at Taishanxi Road, Tianyuan Zone, Zhu zhou City, 412007, China, E-mail: liushuichang100@163.com (corresponding author)

crossref http://dx.doi.org/10.5755/j01.mech.23.6.19847

Received September 30, 2016

Accepted December 07, 2017 Table 1 Boundary setting of the computational domain Boundary setting Remarks Inlet velocity Velocity [partial derivative]=35m/s, Turbulence intensity I=0.5% Outlet pressure Pressure p=101325 Pa, Turbulence intensity I=0.5% Ground Velocity boundary, [partial derivative]=35m/s Top Wall boundary The body surface Wall boundary with no slip