Development of an Algebraic Model of Empirical Parameterization of Near Wakes around a Vehicle

In this paper, we describe and evaluate an algebraic model that has been adopted in a diagnostic wind flow model. Our numerical model is based on the Röckle’s wind modelling approach and we intend to reproduce the steady-state flow patterns of recirculation vortices that are generated in the near-wake region behind a vehicle. The evaluation of the practicality of our proposed model is performed by comparing the wind tunnel experiments of the flow around a vehicle conducted by the Loughborough University. We also compare the numerical results of our model with the CFD model. The proposed model reproduces the flow patterns behind a vehicle and it has significant advantages, such as low numerical costs. We expect that further improvements in the algebraic model when considering the vehicle’s shape will improve its practicality for the numerical analysis of flow fields around a vehicle.


Introduction
Global warming has been one of the major concerns of the world, which needs to be addressed urgently. Alternative fuel vehicles (AFVs), such as hydrogen fuel cells, compressed natural gas, and liquefied petroleum gas vehicles, are expected to be promising vehicles that could serve as counter-measures for reducing the emission of greenhouse gasses. However, in the maritime safety field, significant concerns regarding the fire risks of a pure car carrier or car ferry during marine transportation are expected to increase, owing to the spread of AFVs. Notably, in the vast vehicle space of ships, ventilation efficiency may still be insufficient, owing to the ventilation by hull structures, such as girders, transverse webs, and densely loaded vehicles. Hence, in the case of fuel leakage, flammable gases may accumulate in the vehicle spaces. In order to prevent these undesired situations, it is necessary to conduct ventilation analysis using numerical simulations, when considering the effect of the various obstacles and locations of ventilation equipment in the preliminary design step.
The computational fluid dynamics (CFD) model is well known as one of the most popular approaches to ventilation analysis that provides detailed and precise numerical data. However, the vehicle space of a car carrier is generally vast. Therefore, the CFD model generally requires long computational times, although recent computers have become cheaper and more powerful. Hence, this model is unsuitable in a case study that deals with various leakage scenarios. Accordingly, it is crucial to develop a numerical model with considerable balance in terms of accuracy and computational costs.
As a practical alternative to the CFD, we focus on a mass-consistent wind field (MAS-CON) model, which is a diagnostic wind field model used for wind-flow estimation in applied meteorology (Dickerson, 1978;Sherman, 1978) [1,2]. In the MASCON model of applied meteorology, the atmospheric wind field data are assimilated by correcting initial values, which are estimated by the actual data at the monitoring posts, in order to satisfy the continuity equation. Because the MASCON model solely comprises the continuity equation, it has the advantage of providing numerical results significantly faster than the CFD model, which requires highly intensive iterations to satisfy its continuity and momentum equations simultaneously. However, in the case where the MASCON model is applied for the flow around the obstacle, no initial velocity data are available, especially around the obstacle. Therefore, this model requires operations equivalent to solving the momentum equation in the CFD model; otherwise, the MASCON model must solely provide the solutions of potential flow. To compensate for the lack of the momentum equation, Röckle (1990) [3] developed algebraic models that are expressed by the geometric dimensions of the obstacle. Röckle's model is an empirical approach that is parameterised by the previous experimental data, which enables the estimation of initial values around the obstacle. The diagnostic microscale wind-field (DMW) code [4] provided by Verein Deutscher Ingenieure (VDI) adopted these models and reproduced urban wind flows considering the presence of numerous buildings as obstacles. Recently, a quick urban and industrial complex (QUIC) model [5] was developed by U.S. Los Alamos National Laboratory (LANL), which contains various algebraic models, such as the revision of Röckle's models. This model enables the quick evaluation of the atmospheric dispersion of toxic chemicals in urban areas.
Regarding the parameterisation approach for the flow around a vehicle, an algebraic model of the far-wake that was proposed by Eskridge and Hunt (1979) [6] is well known. Eskridge and Hunt aimed to apply the model to predict the exhaust pollution behavior from vehicles. Accordingly, they developed the algebraic model for estimating the velocity distribution in the far-wake region. Recently, Carpentieri et al. (2012) [7] developed the CAR-BUILD model that roughly estimates the velocity fields and diffusion parameters of the standard Gaussian plume in the far-wake region. They pointed out that algebraic models that are applicable to reproducing the velocity and concentrations in the near-wake region should be developed to further improve prediction accuracy.
In this study, we intend to adopt the MASCON model for the estimation of ventilation flow-fields in closed spaces, such as the vehicle space of a car carrier. Accordingly, an algebraic model of the near-wake of a vehicle has been developed as the first step. The near-wake model has been implemented into the MASCON model that was coded in the Fortran 90 programming language. The reproducibility of the proposed model is verified by comparing numerical simulations of the proposed model with the wind tunnel experiments using 25% scaled models that were performed by Loughborough University [8]. Furthermore, the CFD simulations using an open-source code have been performed, and the practicality of the proposed model is discussed in terms of its prediction accuracy and computation times.

Computational Procedure of the MASCON Model
The computational procedure of the MASCON model (Dickerson, 1978) [1] consists of two steps, namely, prediction and correction steps. In the prediction step, the initial velocity field of the mainstream is estimated by a uniform flow in accordance with the inflow boundary conditions, whereas, in the vicinity of the obstacle, initial velocity is estimated by an algebraic equation that is expressed as a function of the geometrical dimensions of an obstacle. In the correction step, the initial velocity field, which does not satisfy the mass-conservation law, is corrected by a variational method with a constraint on the continuity equation. Here, at any control volume, V, the functional, F, which consists of the total sum of the square of the difference between initial and true values, imposed a constraint on the continuity equation that is multiplied by the Lagrangian multiplier λ and it is defined as where u, v, and w, and u 0 , v 0 , and w 0 represent the velocity components for the x, y, and z directions, respectively, and their initial values. α h and α v denote the weight coefficients of correction in the horizontal and vertical directions, respectively. In this study, the weight coefficients for the horizontal direction are assumed to be the same for the x and y directions. This problem is classified as an optimisation problem with a constraint. Therefore, the first variation of F, δF, is considered to be a stationary value, such that the following Euler-Lagrange equations are obtained: The following boundary conditions are also obtained: where U = (u, v, w) t and n represent the velocity and normal vectors of each boundary of the computational domains. Finally, the above Euler-Lagrange equations substitute in the continuity equation, the following Poisson equation of λ, which corrects the divergence of the initial velocity field is obtained: In the case Dirichlet's velocity condition is given as in a solid (U 0 = 0) and fixed velocity inlet/outlet boundaries, the final velocity must be kept constant as an initial value, such that the following boundary condition for λ is set, for example, in the x-direction: However, in the case where the Neumann boundary condition for velocity is given, such as in free inflow/outflow boundaries, the following condition is imposed for λ outside the boundary: Therefore, the true values of velocity are determined by substituting λ in Equation (2) in order to determine λ by solving Equation (4) under the given boundary conditions expressed as Equations (5) and (6). Most of the computational cost of the MASCON model is spent by solving Poisson's equation of λ to obtain divergence-free flows. In our code, we adopted the staggered-grid arrangement, such that the velocity components are stored at cell faces, and λ is defined at the cell center. Then, two conventional algorithms of the Poisson solver, the biconjugate gradient stabilized (BiCGStab) [9] and successive overrelaxation (SOR) methods [10,11] are implemented and parallelized by OpenMP, which is a popular programming API of the shared-memory type parallelization. Here, we repeatedly state that the solution obtained by the MASCON model is corrected under the constraints with the initial assumptions, and these assumptions are regarded as the best estimates that are closest to the true solution.
Therefore, the method that provides the initial assumption is vital, because the quality of the initial estimate determines the accuracy of the numerical solution of the MAS-CON model. Figure 1 presents the typical flow pattern around a wall-mounted cube. In the QUIC model, the reattachment length L R from behind the cube is determined by the following equation that is derived by Fackrell and Pearce [12]:

Algebraic Models of Empirical Parameterization
where L, W and H denote the length, width, and height of the cube, respectively. The nearwake region of the leeward cube is defined by an ellipsoid shape, as in the following region: where x represents the distance from the obstacle in the streamwise direction. The initial velocity component of x-direction is given by the following equation: where the reference velocity U re f (H) was set to the inflow velocity at the height of the windward cube. d l is the distance from the back of the cube to the arbitrary x position in the near-wake region, while d NN is the distance between the boundary of the near-wake region and behind the cube. Oka (2014) [13], one of the co-authors of this paper, reported Röckle's model and its revisions implemented in the QUIC, which were applied to the flow around the surfacemounted cube in the channel of the parallel flat plates. The obtained results agreed well with the experiment conducted by Hussein and Martinuzzi (1996) [14]. Accordingly, Oka concluded that the algebraic approach with the MASCON model is a significantly costefficient computational method, and it is suitable for the fast simulation of the steady state of the flow field.
In contrast to the wall-mounted cube, a different flow pattern appears in the near wake region of a vehicle. Figure 2 presents the schematic diagram of the flow pattern of a vehicle. Two streams, one flowing over the ceiling (upper stream) and the other under the vehicle (lower stream), are instantaneously detached from the back of the vehicle and merged at a slightly downstream location. At the inner region of the merged streams, two recirculation vortices, rotating in opposite directions to each other, are formed. Conventional models can only apply a single recirculation vortex. Therefore, no model can reproduce the two recirculation vortices that are generated by the two upper and downstream streams behind the vehicle. In this study, we have developed an algebraic model of the two recirculation vortices that are applicable to the near-wake region based on the revised Röckle model. When considering the elliptical shape of the near-wake region behind the vehicle in the x-z section and the two recirculation vortices formed in the upper and lower halves of the ellipse, Equation (8) is modified for the near-wake region, as follows: where H f and H v represent the floor height and distance between the ground and top of the vehicle's rear, respectively. Similarly, for the geometric dimension H shown in Equation (7), which determines the lengthL R of the near-wake region, considering the shape of the rear of the vehicle generating the recirculation vortices, is modified in the following equation: The initial velocity component was obtained in accordance with Equation (9), similar to the conventional Röckle model. In addition, we tentatively implemented the far-wake model of the QUIC model that was reported by Singh et al. (2006) [15], which is a modified version of Röckle's model. Because the virtual origin of the velocity deficit equation was not explicitly defined in the paper by Carpentieri et al. (2012) [7], the correct results have not been obtained in our trial implementation. In the QUIC model, the far-wake region is defined at the back end of the near-wake region and within the region that is given by the following equation: Regarding the initial estimate of the velocity, d N denotes the distance between the boundary of the far-wake region and the back of the cube, and it is similar to the near-wake model, and it is defined, as follows:

Computational Fluid Dynamics Model (CFD)
In this study, the OpenFOAM is adopted as the CFD model. OpenFOAM, which is the CFD toolbox written in C++ language (available at https://github.com/OpenFOAM accessed on 21 November 2020), is one of the most popular open-source codes in the world. The OpenFOAM model is parallelized by an efficient distributed-memory parallel programming via message passing interface (MPI). The OpenFOAM version 7.0 (The OpenFOAM Foundation, London, UK), which is the stable and previous version, is adopted for the incompressible steady flow simulation.

Test Cases of Wind Tunnel Tests with 25% Scaled Model of the DrivAer Model
The DrivAer model is a generic vehicle model developed by the Institute of Aerodynamics and Fluid Mechanics at the Technische Universität München [16] for the aerodynamic investigations of passenger vehicles, as shown in Figure 3. The dimensions of the full-scale model of the DrivAer model are 4613 mm in length, 1820 mm in width, and 1418 mm in height. In this study, wind tunnel tests with a 25% scaled model of the DrivAer model, which were conducted by the Loughborough University and reported by Varney et al. [8], are adopted as the test cases for the numerical simulation. Hence, the dimensions of the scaled vehicle model are 1153.25 mm in length, 455 mm in width, and 354.5 mm in height. The wind tunnel facility consists of a test section with dimensions of 3600 mm × 1920 mm × 1320 mm (length × width × height). The representative velocity is set to 40 m/s. Here, the kinematic viscosity of air is 1.512 × 10 −5 m 2 /s, and the total length of the vehicle is assumed to be a representative length, whereas the Reynolds number (Re l ) is calculated as 3.05 × 10 6 . Our algebraic model is applicable when the hypothesis of the Reynolds number independence is satisfied. The Reynolds number independence appears between 11,000 and 99,000 in the case of the surface-mounted cube described in the EPA Guideline [17]. In the blunt-shaped Bus case, Solmaz et al. [18] investigated the Reynolds number independence on the drag coefficient experimentally and via CFD, and they reported that it appeared approximately more than 55,000. Strangfeld et al. [19] reported that the results of a wind tunnel experiment with the 25% scaled fastback model of the DrivAer model indicated a Reynolds number independence in the range of Re l ≥ 2.0 × 10 6 . Because our algebraic model intends to reproduce the flow pattern and not for the aerodynamic study, the applicable range of Re l in this study may be wider than that of the aerodynamics case. However, we expect that 10 5 is the minimum value of the applicable range of Reynolds number.

Computational Conditions
Computational domains are set similar to that of the wind tunnel experiment's test section, except for its length, in order to obtain sufficient numerical data in the wake region, as shown in Figure 4. The validation cases of the DrivAer model provided by Wolf Dynamics (available at http://www.wolfdynamics.com/tutorials.html?id=152 accessed on 21 November 2020) are adopted as a basis, and they are modified to be suitable for our calculation. A number of cells are set to approximately 3.8 × 10 6 cells for the half model with a cut through the symmetry plane and center of the y-direction in the CFD simulation, which is 960 mm in width as a half of the wind tunnel.
In the CFD simulation, the base mesh that consists of 0.2 m of regular grids is generated using blockMesh application and refined to five levels of fineness using the snappy-HexMesh application. Figure 5 shows the schematic diagram of the mesh configuration. The finest mesh, denoted as Lv 5 in Figure 5   The experimental data of inflow velocity were reported by [20]. Subsequently, to assume the boundary condition of inlet velocity, we adopt the Blasius boundary layer profile following a 1/7th power-law with a boundary layer thickness (δ) of 60 mm, as illustrated in Figure 6. Tables 1 and 2 summarise the boundary conditions in the CFD simulation.  In the CFD simulation, the pressure-velocity coupling is performed by using the SIMPLE algorithm implemented in simpleFoam. The k-omega SST turbulence model [21] is adopted for the Reynolds Averaged Navier Stokes (RANS) simulations. The average, minimum and maximum values of the y + are 20.917, 1.0346 and 266.46 in the estateback model, respectively. Although the k-Omega SST turbulence model should be used with y + < 2, the near-wall treatment developed by Menter et al. [22] enables the computed wall shear-stress varies by less than 2%, and all of the solutions follow the logarithmic profile in the wide range of the y + values from 0.2 to 100. Accordingly, we consider that the mesh configuration of the CFD simulations generally satisfies the above recommendation for practical applications. Table 3 summarises the solver settings of the CFD simulation. The pressure equation is solved by the GAMG solver with the Gauss-Seidel smoother and the tolerance and relative tolerance are set to 1 × 10 −6 and 0.01, respectively. On the other hand, the PBiCGStab solver with the DILU preconditioner is utilised to solve the rest equations, such as U and turbulence equations. The tolerance and relative tolerance are set to 1 × 10 −8 and 0.001, respectively. In addition, because we assume a vehicle is in parking condition, the motion of wheels are changed to the stationary state from the original rotational configuration of the Wolf dynamics. bounded Gauss upwind div(phi,omega) bounded Gauss upwind div((nuEff*dev(grad(U).T()))) Gauss linear Laplacian Gauss linear limited 0.5

Results and Discussion
In this section, we present the results and comparisons of the simulations that were obtained from two models and experiments, and we then discuss their practicality in terms of reproducibility and computational cost.
First, we investigate the optimal mesh resolution of our proposed method on the reproducibility of recirculation vortices. Figure 7 shows the comparison between the initial and corrected velocity fields of the proposed model. In the figure on the left (Figure 7), the initial velocity vector is given in the opposite direction to the mainstream around the recirculation region. After the correction by the MASCON model, as shown in the right figure, it can be observed that the velocity vectors in the recirculation region are varied in a swirling manner. Figure 8 presents the results of the simulations that were obtained from the proposed method, with varying mesh sizes. Streamlines of recirculation vortices are clearly illustrated in the near-wake region along with the increment of the mesh resolutions. In the case of dx = dy = dz = 0.025 m, the recirculation vortices are unclear. Therefore, we determine that the minimum resolutions of the mesh size has been 0.0166 m, while considering the numerical costs, as described in Table 4.    [23] pointed out that the RANS analysis for the DrivAer model using k-omega SST turbulence model experienced fluctuations due to a pseudo-unsteadiness in the flow. We consider that the pseudo-unsteadiness appears in our simulations. However, to avoid choosing the result of an arbitrary step, we reasonably calculate the numerical solutions by averaging over the additional 1000 steps from the end of the simulation in the same way as past study conducted by Heft et al. (2012) [16]. Figure 9 presents the comparisons of the streamlines and velocity profiles in the nearwake region between the experiment and two models. It can be observed that there are two recirculation vortices in the near-wake region for the experimental data. In the CFD result, the lower vortex appears successfully, whereas the upper vortex is unclear. On the other hand, although our proposed method successfully reproduces the two recirculation vortices, the streamline trajectories are slightly different from the experiments. To elucidate this result, we compare the vertical profile of the velocity at the four positions of the x-direction. The magnitude of the velocity deficit, which is expressed as the difference between the velocity profile of the approaching flow and actual velocity, is reproduced well by both CFD and the proposed method, but the shape of the distribution is not optimally reproduced. Therefore, it is inferred that the proposed method can generally reproduce the velocity profile of the experiment in the near-wake region.  Figure 10 presents comparisons of the streamlines and velocity profiles in the nearwake region between the experiment and two numerical models in the case of the fastback model. In the results that were obtained from the experiment, two recirculation vortices also appear in the near-wake region. In contrast, two numerical models successfully reproduce two recirculation vortices that are similar to the experiment. However, the length of the near-wake region is over-estimated in the downwind direction of the two models. Although the trajectory of the near-wake region, which is the blue colored region presented in Figure 10a-c, extends towards the ground in the experiment, the near-wake region shows a horizontal extension in the two numerical models. These differences are also observed in the vertical profile of the local velocity, as illustrated in Figure 10d. These differences may be induced, owing to the alterations in the direction of the near-wake flow by the shape of the vehicle's rear. The proposed algebraic model solely considers the vehicle's geometric dimensions and it does not consider the effect of the vehicle's actual rear shape. Hence, we conclude that the proposed algebraic model requires further improvement in the inclusion of the effects of the vehicle's rear shape.  Figure 11 presents the comparisons of the streamlines and velocity profiles in the near-wake region between the experiment and two numerical models. In the experimental result, two recirculation vortices are formed in a similar manner to the previously discussed two vehicle models. Regarding the result of the CFD model, although the trend in the lower recirculation vortex is slightly different from that of the experiment, the simulation generally reproduces the experiment. On the other hand, the formation of two recirculation vortices in the near-wake region is also observed in the results that were obtained from the proposed model. In addition, in the near-wake region, the differences in velocity deficit in the vertical velocity profiles appear to be similar to those discussed for the fastback model. Therefore, we infer that modifications to our algebraic model, as concluded in the Fastback model subsection, will also improve the reproducibility of the flow of recirculation vortices in the near-wake region for the Notchback model.

Practicality
Finally, we discuss the practicality of the two models in terms of the accuracy and computational costs.

Accuracy
Regarding accuracy, 143 points of data are sampled in the near-wake region, 11 from 0.95985 m to 1.251 m in the x-direction, and 13 positions from 0.01591 m to 0.34586 m in the z-direction. Local velocity data that are normalised by the representative velocity, U ∞ = 40 m/s, are compared between the experiment and the two models, as presented in Figure 12. If the plots are located closer to a diagonal line, it indicates that the model reproduces exactly in the experiment. It can be seen that the CFD model is better than the proposed method in terms of accuracy, with the overall plots being closed to diagonal lines in all three vehicle models. On the other hand, in the proposed model, the flow pattern of the estateback model can be reproduced better than the other two vehicle models. Figure 13 plots each sampled position colored in accordance with the difference in values between the experiment and the two numerical models. Because our algebraic model does not accurately consider the shape of the recirculation vortices, the error of the proposed method is particularly distributed at the boundary between the near-wake region and main flow, as indicated by the deeper red plots. To statistically and quantitatively evaluate these results, the root mean square error (RMSE) of the difference between the experimental data and the two numerical models is calculated and summarised in Table 5. For each of the vehicle models, the CFD model shows less error from the experimental values than from the proposed model. However, calculations using the proposed model verify that the estateback model indicates the least error and most reproducible results among the three models.   Figure (a-c) illustrate the estateback, fastback, and notchback models obtained from the results of the CFD model. In addition, Figure (d-f) depict the estateback, fastback, and notchback models that were obtained from the proposed model. Regarding computational costs, Table 5 summarises the elapsed times for the computation. In the case of the CFD model, residues of velocity, pressure, turbulence energy, and omega in the iteration have not been reduced at approximately after 1000 steps, as illustrated in Figure 14, such that the CFD calculations have not reached the steady state, despite the elapses of 100,000. The CFD model has taken approximately 300 times longer than the proposed method to obtain the results at 5000 steps. Therefore, we have concluded that the proposed method can be more practical than the CFD model in terms of computational costs.

Conclusions
In this study, we have developed an algebraic model for reproducing the flow pattern in the near-wake region of a vehicle, and evaluated it by comparing with the experiment. We have verified its practicality in terms of accuracy when compared with the experiment using a 25% scaled wind tunnel. In addition, we have confirmed and compared the computational costs of the CFD and our algebraic model with the MASCON model. Consequently, the conclusions drawn from this study are summarised, as follows: • In general, our algebraic model that was implemented in the MASCON model successfully reproduces two recirculation vortices in the near-wake region of the vehicle for the three vehicle shapes. • The proposed method reproduces the flow pattern of the estateback model better than the other two vehicle models. However, the CFD model reproduces the flow pattern of the experiments better than our proposed method in the fastback and notchback models. • In the CFD model, the residues are almost constant from approximately 1000 iterations, and they have never reduced below the criterion for determining the steady-state solution. Therefore, we have regarded the solution by averaging over the additional 1000 steps from the end of the simulation. If we assume that the time elapsed to reach 5000 steps is the computational times in the CFD model, then it takes approximately 300 times longer than the proposed method. • The accuracy of reproducibility for the fastback and notchback while using the proposed method may be improved by revising the algebraic model to account for the rear shape of the vehicle. • We conclude that the proposed method is practical in terms of the computational costs. Therefore, the algebraic model combined with the MASCON model would be a promising model that alternates the CFD model, by further improvement of the reproducibility.
Finally, we will continue to improve the accuracy of the algebraic model in the nearwake region. Accordingly, we will develop a model that can be applied to the far-wake region in order to more accurately reproduce the flow pattern around the vehicle.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

OpenMP
Open