Computational Engineering and V&V manual

A. Laser beam welding of AISI SS 304

In this case study, macroscale modeling including laser dynamics, melt pool dynamics, phase change, solidification, Marangoni effects, heat, and mass transfer, and surface tension effects are captured [1].

This validation study compares the depth of penetration and weld width reported in experiments [2] for the laser beam welding process of AISI SS 304 with the numerical results.  

Different sets of laser power and welding speed parameters are used and the consequent melt pool dimensions are reported in experiments summarized in Table 2A.1. The welding operation is performed on a 1.6 mm thick austenitic stainless steel (AISI 304 SS). Argon gas is used to avoid atmospheric contamination. The length and breadth of the steel are 60 mm and 20 mm, respectively. 

Table 2A.1. Laser parameters used for the simulation

Figure 2A.1. Comparison between the experimental (left) and simulated (right) results for AISI 304 SS for laser power of 1250 W and welding speed of 750 mm/min

Figure 2A.2. Simulated results for AISI 304 SS for (a) laser power of 1000 W and welding speed of 1000 mm/min and (b) laser power of 1250 W and welding speed of 1250 mm/min

Figure 2A.1 illustrates the excellent level of similarity between the experimental image of the melt pool and the numerically generated melt pool. Deep penetration is observed in experiments and numerical estimates where the heat is immediately transported to the material volume. Figure 2A.2 shows the melt pool profile for different values of laser power and welding speed. 

Figure 2A.3. Histogram plot showing the comparison of depth of penetration between experimental, AM PravaH, and another CFD software for the four cases of laser parameters

Figure 2A.4. Histogram plot showing the comparison of weld width between experimental, AM PravaH, and another CFD software for the four cases of laser parameters

Figures 2A.3 and 2A.4 show the histogram comparing AM PravaH and another CFD software with the experimental data. AM PravaH results are closer to the experimental data for both the depth of penetration and the weld width.  Table 2A.2 shows that the error percentages for depth of 

penetration and weld width are under 4% and 1% respectively for AM PravaH - an excellent validation.

Table 2A.2. Comparison of the depth of penetration and the weld width results with experiment for different sets of laser parameters

The numerically obtained melt pool widths and depths of penetration are in good quantitative and qualitative agreement with the experiment.

B. Melt pool evolution during laser powder bed fusion of Ti-6Al-4V

This validation study compares the melt pool width and depth reported in experiments for the laser powder bed fusion (LPBF) process for Ti-6Al-4V with the numerical results [3].

The laser power is kept at 300 watts, and the scanning speed is at 1750 mm/s. The laser power along with the scanning speed, hatch spacing, and layer thickness are mentioned in Table 2B.1. Figure 2A.1 depicts the computational domain considered for the simulation where the Ti-6Al-4V powder particles are distributed randomly over the metal substrate.

Table 2B.1. process parameters for LPBF modeling

Figure 2B.1. Schematic showing the computational domain with random packing of Ti-6Al-4V powder bed

(a)

(b)

Figure 2B.2. (a) Isometric view of melt pool progression and temperature profile map during deposition of LPBF, (b) longitudinal section of melt pool at Laser power: 300 W, scanning speed: 1750 mm/s

Figures 2B.2 (a) and (b) show the melt pool progression for the LPBF model for both isometric and longitudinal view. As the laser moves on top of the powdered bed, it is seen that the powder melts. Figure 2B.2 (b) illustrates the effect of Marangoni surface tension forces and proper evaluation of the melt pool. Comparing the experimental and simulation results in Table 2B.2, the agreement between the values is excellent where the overall error percentage is less than 4%.

Table 2B.2. Molten pool width and depth comparison for experiment and simulation

The numerically obtained melt pool widths are in excellent quantitative agreement with the experiments.

C. Melt pool evolution during laser powder bed fusion of Inconel 718

This validation study compares the melt pool width and depth reported in experiments for the laser powder bed fusion (LPBF) process for In718 with the numerical results [4].

The laser power is kept at 300 watts, and the scanning speed is at 1500 mm/s. The diameter of the particles is kept at 30 μm and the layer thickness of the powder bed is kept at 30 μm as mentioned in Table 2C.1. Figure 2C.1 depicts the computational domain considered for the simulation where the In718 powder particles are distributed randomly over the metal substrate.

Table 2C.1. Process parameters for LPBF modeling

Figure 2C.1. Schematic showing the computational domain with random packing of Ti-6Al-4V powder bed.

(a)

(b)

Figure 2C.2. (a) Isometric view of melt pool progression and temperature profile map during deposition of LPBF, (b) longitudinal section of melt pool at Laser power: 200 W, scanning speed: 1500 mm/s

Figures 2C.2 (a) and (b) show the melt pool progression for the LPBF model for both isometric and longitudinal views. As the laser moves on top of the powdered bed, it is seen that the powder melts. Figure 2C.2 (b) illustrates the effect of Marangoni surface tension forces and proper evaluation of the melt pool. Comparing the experimental and simulation results in Table 2C.2, the agreement between the values is excellent where the overall error percentage of width and depth are 1.8% and 2.7% respectively.

Table 2C.2. Molten pool width and depth comparison for experiment and simulation

Due to stringent rules posed by the government for reducing carbon emissions, the automotive industry is facing a shift from internal combustion engine-driven vehicles to electric vehicles. There are many challenges that need to be overcome to match the performance of electric vehicles to that of conventional vehicles. The major challenge lies in designing an efficient Battery thermal management system. This is applicable to hybrid as well as fully electric vehicles from two-wheelers to cars and heavy multi-utility vehicles. 

A. Natural convection in a sealed battery pack

1. Introduction

In this case, the transient thermal response of a 15-cell, 48 V, lithium-ion battery pack for an unmanned ground vehicle (UGV) is simulated. Results were then compared and validated against experimental results [15]. Heat generation rates and specific heat capacity of a single cell were experimentally measured [5] and used as input to the thermal model. A heat generation load was applied to each battery, and natural convection film boundary conditions were applied to the exterior of the enclosure. The buoyancy-driven natural convection inside the enclosure was modeled along with the radiation heat transfer between internal components. The maximum temperature of the batteries reached 65.6 ⁰C after 630 s of usage at a simulated peak power draw of 3600 W or roughly 85 A. 400 -170-66-90

The goal of this analysis was to determine the maximum temperatures the batteries would experience under a prescribed high power loading of 240 W per battery or about 85 A. 

Figure 3A.1. CAD model top section view in .stl format 

Figure 3A.2. CAD model side section view in .stl format

2. Simulation details

The material properties of the batteries, aluminum, plastic, and air at 20 ⁰C are gathered from the literature [5]. The viscosity of the air in the model was 1.85105 kg/m-s, and the air thermal expansion coefficient was 0.00343 1/K. Heat generation in the batteries was modeled as uniformly distributed. The heat generation rate calculated theoretically using the Entropy method through experimentation in the report was used as input for the volumetric heat source in the simulation [5]. A constant value of 31.97 W was applied individually as an absolute heat source (Q) to all the cells in the battery pack.

A no-slip wall condition was imposed at the inner enclosure walls, battery walls, and plastic walls where the enclosed air comes in contact with the solid surface. At the start of the simulation, all the components of the system were initialized to atmospheric pressure and a temperature of 25 ⁰C. The external free surfaces of the enclosure were given surface film boundary conditions to simulate free convection on the vertical and horizontal surfaces. The sink temperatures for the film boundary conditions were set to an ambient temperature of 25 ⁰C. To be conservative, a smaller-than-predicted, uniform convection coefficient boundary condition of 2 W/m2 K was applied to the exterior surfaces of the pack.

A transient analysis was required in the present case since the simulation only reached a steady state long after the battery pack would be depleted. Experiments showed that the batteries were depleted after 630 s at the 5 P rate. A variable adjustable time step with a maximum Courant number of 0.95 and a maximum Diffusion number of 10 was used to run a 630 s flow time study. This sufficiently small Courant number ensured the accuracy of the results and was based on the characteristic time of the modeled system.

The Boussinesq approximation was used by assigning the modeled air an operating density of 1.205 kg/m3. The Boussinesq approximation neglects any changes in the density of the fluid, except when the density terms appear multiplied by the gravity term. This leads to the conclusion that inertia differences are negligible while gravity is large enough to create a buoyancy-driven flow.

3. Results and Discussion

Figure 3A.3. Temperature back view (oC) after 630 s at 5P discharge rate [Left - ANSYS, Right - OpenFOAM]

Figure 3A.4. Temperature top section view (oC) after 630s at 5P discharge rate [Left - ANSYS, Right - OpenFOAM]

Figures 3A.3 and 3A.4 illustrate that the aluminum enclosure stayed at a relatively uniform and near ambient temperature throughout the battery discharge. This result is significant as it illustrates that there is a lack of heat flow to the enclosure which would have transferred thermal energy away from the batteries. This results in a high-temperature increase of the battery pack inside the aluminum enclosure.

Figure 3A.5. The vertical plane of velocity vectors showing temperature (oC) for air in the enclosure after 630 s at 5P discharge rate [Left - ANSYS, Right - OpenFOAM]

Figure 3A.5 shows an image of the buoyancy-driven flow in the pack, plotted as vectors. For the sake of clarity, the batteries are removed from this image. A plume of hot air is observed collecting at the top of the enclosure and relatively colder air at the bottom. This figure also shows that the air rises along the battery walls and flows downward along the relatively cool inside walls of the enclosure. Thus a convection current is set up inside the enclosure. Hot air collected from the battery goes to the top. Then it comes in contact with the cool enclosure walls and again flows down to the bottom of the enclosure. Thus air keeps circling in the space between the battery and the enclosure.

4. Validation and Conclusion

The CFD model was validated by comparing the surface temperature of the battery at the center of the pack to the experimentally measured surface temperature of the battery performing the same discharge profile while insulated by polyurethane. Figure 3A.6 shows that the surface temperature in the model followed the same trend and had a slightly smaller magnitude compared to the experimentally measured temperature. The temperature rise at the end of the simulation was completely conformal till 7 minutes. The results matched closely as expected since the modeled battery pack has near insulating conditions due to the close proximity of other batteries and the weak buoyancy-driven flow predicted by the low Rayleigh number.

Figure 3A.6. Middle battery surface temperature (oC) comparison between simulation, and experiment

Figure 3A.7 shows that simulated and experimental are almost conformal till 7 minutes of simulation time. After 7 minutes the predicted CFD values are lower than the experimental values. Percentage mean error was calculated for all the data points. The following Table 3A.1 shows the percentage mean error value comparison between simulated, and experimental. 

Figure 3A7. Middle battery surface temperature (oC) CFD results error between simulation and experiment

Table 3A.1. CFD results in error analysis

Simulated results were found to be 98.27% accurate when compared with the experimental values.

B. Phase change material cooling in a single cell

1. Introduction

Li-ion energy storage is exponentially increasing because of its high energy density, capacity, high cycle life, and low self-discharge rate. However, when the battery is operational it generates heat at the cell level, and the temperature can increase beyond 40 ⁰C. This increase in temperature affects the battery performance and life adversely.

Figure 3B.1. Heat transfer mechanism in the phase change material battery cell cooling [6]

Phase change material cooling is an efficient, light weighted, compact, low-cost, uniform temperature distribution, energy-intensive way of battery cooling. 

Figure 3B.2. Li-ion cell PCM module (a) schematic (b) CAD model (top view), (c) schematic of the cell core (d) isometric view of cell with PCM module 

Figure 3B.3. The meshing of the simulated model (a) top view (b) isometric view

2. Simulation Details

The CFD model of phase change material cooling is validated by comparing the average temperature of the PCM module. This includes the cell as well as the filled PCM domain [7].  Figure 3B.4 shows the temperature profile of the PCM module at 2000 s. The maximum temperature of the cell reaches up to 312.7 K at 2000 s and starts to saturate at the same temperature as shown in the plot (Figure 3B.7).

Figure 3B.4. Temperature profile at 2000 s (a) top view (b) isometric view

3. Results and Discussion

Figure 3B.5 shows the fraction of PCM that is in direct contact with the cell and starts melting after experiencing the melting point of the PCM i.e. 311.15 K. A constant heat was provided to the system corresponding to the current rate of 2C. It was observed from previous studies that the latent heat of fusion and thermal conductivity plays a vital role in dissipating the heat and in turn improving the performance of a battery. 

Figure 3B.5. The PCM melted fraction surrounded by the cell 

Figure 3B.6 represents the temperature variation only in the PCM material. The maximum temperature experienced by the PCM is 312.3 K, which is lower than the cell temperature.  Yet, it is efficiently absorbing the heat from the cell. The region that experiences a temperature greater than 311.15 K starts melting and the fraction is shown in Figures 3B.5 and 3B.6 (b).

Figure 3B.6. Computational results for the PCM domain (a) Temperature profile (top view) (b) PCM melted fraction (top view)

Figure 3B.7 shows the graphical representation of the average temperature inside the domain including the cell and the PCM that starts saturating after 1800s.

Figure 3B.7. The average temperature of the PCM module with time  

4. Validation and conclusion

Figure 3B.8 shows a comparison plot of numerical analysis performed with respect to the experimental data. The simulated result shows a good agreement with the experimental data [16]. The error analysis was performed and found the simulated result showed promising results with an error of only 0.13 %, shown in Figure 3B.9.

Figure 3B.8. Validation plot of phase change material cooling for a single battery cell. Inhouse Paanduv results were compared with experimental data.

Figure 3B.9. Error analysis for simulated results while comparing it with the experimental data obtained from Chaudhari et. al. [7]

C. Phase change material cooling in a Prismatic Cell

1. Introduction 

The electrification of transportation is crucial for addressing the challenges of greenhouse gas emissions and energy scarcity, aligning with the goals of the Paris Agreement that took effect in November 2016. Lithium-ion batteries (LIBs) are widely used in electric vehicles (EVs) and hybrid electric vehicles (HEVs) due to their high specific energy, long cycle life, and low self-discharge rate, making them ideal for energy storage in these vehicles.

Despite the rapid advancement of electric vehicles (EVs), challenges remain regarding vehicle lifespan and dynamic performance. One critical factor influencing the Li-ion battery's safety, performance, and lifespan—the primary energy source for EVs—is its operating temperature.

Elevated temperatures, resulting from continuous heat buildup during cycling, can accelerate the depletion of active battery materials and increase the risk of thermal runaway. The acceptable operating temperature for LIBs is typically -20 °C to 60 °C.

Efforts are underway to develop more efficient and advanced Battery Thermal Management (BTM) systems to address these challenges. These systems can be categorized based on the cooling medium used, including air cooling, liquid cooling, PCM cooling, and combinations thereof.

PCM cooling, which utilizes latent heat during the phase change process, is an innovative solution for BTM systems.

2. Simulation Details

The validation of the CFD model for the phase change process in the prismatic cell involves comparing the maximum temperature of the battery module, focusing solely on the battery module itself [EV1]. 

Fig.3.C.1 CAD model of the Prismatic Battery cell

Following the validation of the module with experimental results, three different cases with varying PCM configurations are tested.

Fig. 3.C.2. Schematic showing different configurations of the PCM used in the simulations (case 1: PCM on the top, case 2: PCM on the sides case, 3: PCM on top and the sides)

• In this study, a prismatic LiFeO4 battery is taken as the sample for the test. The nominal cell capacity is 12 A h and the dimensions are (height X width X depth) 90 mm X 70 mm X 27 mm. Fig. 3.C.1.

• The convective heat transfer coefficient for cell and PCMP with the value of 5 Wm-2K-1 and 2.25 Wm-2K-1 under the condition of natural convection. The initial and ambient temperatures are set as 298.15 K.

• Heat generation rates for cell core, positive and negative terminals during 5C discharge rate are 12.70e4 W/m3, 26.52e4 W/m3, and 16.39e4 W/m3, respectively.
  
  

To study the impact of PCMP configuration on thermal performance, three cases are developed in this section, as illustrated in Fig. 3.C.2. These cases have identical PCMP volumes, with dimensions a, b, and c set at 2 mm, 5.18 mm, and 1.4 mm, respectively.

Fig. 3.C.3. Meshing of the computational domain

3. Results and Discussions:

To validate the developed battery model, the surface temperature evolution of the battery is monitored during a 5C discharge. The simulation involves a high discharge rate of 5C, lasting 720 seconds.

Fig. 3.C.4. Illustrates the comparison between the maximum surface temperature between the experimental results and the simulation results.

Fig. 3.C.4.  Comparison between experimental and simulated results of battery cell 

Fig. 3.C.5. illustrates the evolution of T max under various PCMP configurations. After the 5C discharge, the T max of the battery without PCMP reaches 336.35 K. In contrast, for Case 1, Case 2, and Case 3, the T max values are 326.91 K, 330.46 K, and 326.16 K, respectively. This represents reductions of 9.79 K, 6.25 K, and 10.54 K, respectively.

Fig. 3.C.5. Temperature plots of different PCM configurations at the end of 720 seconds. 

A. Single-phase reaction kinetics

1. Introduction

Water splitting reaction assisted by solar light is a widely explored pathway to generate hydrogen. However, multiphase reactions such as this one are tough to model as include liquid and gas both in the reaction and their transport. These reactions usually take place at room temperature at atmospheric pressure in presence of the catalyst either in slurry form or thin film coated at a substrate. The gas and bubble nucleation starts at the catalyst site. In this case, CFD modeling is performed and an appropriate multiphase solver is validated providing the zero-order reaction rate for the H2 generation reaction [8, 9]. Since the solid phase is not considered the mass transfer effects are neglected while modeling the system. The bubble dynamics studied in the system resulted from the multiphase reaction where the evolution of H2 and O2 is happening by the given reaction rate. 

Water splitting reaction kinetics

Reaction kinetics is governed by Arrhenius law as below

Numerical validations of water splitting

a.    Reaction kinetics only

The verification below shows the comparison of the numerically calculated time duration of the water decay with that of the analytically calculated value (t= 5.5x106 s) for a pure (Navier-Stokes ignored) water splitting reaction. 

The exact thermophysical properties of water, hydrogen, and oxygen such as molecular weights, viscosity, and density are used.  The reaction is modeled at 300 K and atmospheric pressure.

Figure 4A.1. Comparison of H2O decay obtained from simulation vs analytical solution for a pure reaction kinetics system

The results (Figure 4A.1) are in excellent agreement.

b    Reaction kinetics in a single phase

The verification below shows the comparison of the numerically calculated time duration of the water decay with that of the analytically calculated value (t= 5.5106 s) for water splitting reaction in a single phase (Navier-Stokes included). 

The exact thermophysical properties of water, hydrogen, and oxygen such as molecular weights, viscosity, and density are used.  The reaction is modeled at 300 K and atmospheric pressure.

Figure 4A.2. Comparison of H2O decay obtained from simulation vs analytical solution for a reaction kinetics system in a single phase

The results (Figure 4A.2) are in excellent agreement.

c.     Numerical validation of water splitting reactions kinetics in a multiphase system with phase change 

In this section, backed with the verifications, results obtained from the modeling of a full-fledged water-splitting reaction in a multiphase (liquid-gas) system are reported. The splitting reaction happens inside the liquid phase containing water as a reactant and the products (hydrogen and oxygen) are released into the gaseous phase. Therefore, phase change modeling is also included.

The verification below shows the comparison of the numerically calculated k=0.7 10-5 kmolm3- s of water decay with that of the analytically calculated value of 110-5 kmolm3- s for a water splitting reaction in a multiphase system with phase change from liquid to gas. The exact thermophysical properties of water, hydrogen, and oxygen such as molecular weights, viscosity, and density are used. The reaction is modeled at 300 K and atmospheric pressure. The initial state of the system is shown in Figure 4A.3.

Figure 4A.3. Initial water-air system (1 m x 2 m)

The results (Figure 4A.4) are in good agreement.

Figure 4A.4. Comparison of rate constant obtained from simulation vs analytical solution for a reaction kinetics system in a multiphase and phase-changing setup

The hydrogen generation can vary based on the nature of the flow (laminar or turbulent) and the dynamicity (stagnant or disturbed). 

Different combinations of these conditions cause varying flow patterns and change the extent of bubble formation (Figure 4A.5). Therefore, three variations of the same problem are modeled - laminar, turbulent, and disturbed system - with a fine mesh.

Figure 4A.5. Comparison of water splitting reacting flow dynamics for the laminar, turbulent, and disturbed systems

1. The laminar system: No turbulence is used.

2. The turbulent system: The k-𝜺 turbulence model is used with a 5% turbulent intensity.

3. The disturbed system: The bottom wall is disturbed in ∓x direction with a random velocity magnitude of under 0.5 m/s. 

In all the scenarios (Figure 4A.6), a fairly steady rate of generation of hydrogen is observed with complex mixing, decaying water, and varying velocity profiles. 

Figure 4A.6. Hydrogen production for the laminar, turbulent, and disturbed systems respectively.

The amount of hydrogen becomes steady inside the control volume after some time as the rate of production of hydrogen becomes almost equal to its escape from the control volume.

B. Flow through Kenics static mixer 

1. Introduction

The flow and mixing behavior of two miscible liquids is simulated in a Kenics static mixer, in Reynolds number 160 flow regime. The performance of the Kenics mixer is numerically evaluated and validated with verified experimental literature [10]. The pressure drop ratio (Z-factor) and coefficient of variation (CoV) features have been used to evaluate the mixing performance in Kenics static mixer. Numerous empirical formulas are already available in previous literature for Z-factor calculation [10]. The model is firstly validated based on experimental data measured for the pressure drop ratio and the coefficient of variation. CFD results are consistent with measured data and those obtained by available correlations in the literature. 

Figure 4B.1. CAD model (Kenics Mixer) isometric view in .stl format

Figure 4B.2. CAD model (Kenics Mixer) side view in .stl format

Figure 4B.3. Isometric view of the meshed domain 

Figure 4B.4. Side view of the meshed domain 

2. Simulation Details

The general assumptions of the model used are as follows: the flow is single-phase, three-dimensional, steady state, incompressible, isothermal, and without any reaction (without mass production or mass consumption). 

Re=ρVDμ 

where V stands for the mean velocity, D for the diameter of the tube including the mixing elements and mu stands for the dynamic viscosity of the fluid. 

In the numerical modeling of a static mixer, a no-slip boundary condition is used on the tube wall and element surfaces. The atmospheric pressure, i.e., zero gauge pressure, is utilized at the outlet of the static mixers. Brine 25 % (wt.) is injected into the tube co-currently to the main flow. The tracer (brine) velocity chosen at each Reynolds number (based on the inlet condition of the main liquid) corresponds to a ratio of 0.1 between the tracer flow rate and the main flow rate. Velocity at water and brine inlets is calculated using equation (1) and then given as input in the model.

The accuracy of the numerical model is an essential aspect to be quantified prior to using it for predicting new situations. To assess the validity of the CFD results, the simulated pressure drop ratio and the coefficient of variation are compared with the experimental data and available correlations. The Z-factor is defined as the ratio of the static mixer pressure drop, ΔP, to the pressure drop in the empty tube ΔPo. In this study, the pressure drop is calculated as the absolute difference between the area-weighted average pressure at the planes located at 1.5D upstream and downstream of the mixing elements. This is an important parameter for static mixers and provides a measure of the energy required for the static mixer in the pipeline. 

Moreover, the coefficient of variation is calculated at the outlet face of the static mixer. alpha.brine (volume fraction of brine) quantity was integrated on the outlet face. Then the point data of alpha.water (volume fraction of water) on each mesh cell of the outlet face was exported to an excel file. The mean was calculated for all the data points using the mean value function in excel. After that, the standard deviation was calculated by using the function inbuilt in excel. The formula of covariance is the standard deviation divided by the mean. Thus the ratio of these two quantities gives us the value of covariance which is the CoV of brine at the outlet. CoV with a value of 0 means full mixing. CoV with a value of 1 means no mixing. The value of CoV lies between 0 and 1. The value of CoV should be less than 0.05 for the static mixer to pass as an experimental design. 

Many reports also compare the values computed for Z-factor in the range Re=160 and for CoV in the range Re=160 with the data measured in the Kenics mixer [11] and in a static mixer with baffle shape elements [12], respectively. From the literature, it is deduced that an increase in the Reynolds number leads to an increase in pressure drop, which in turn leads to an increase in the Z-factor.

3. Results and Discussion

The empirical expressions that are commonly used for estimating the Z-factor in Kenics static mixers are listed in Table 4B.1. There is a considerable discrepancy in the results predicted by the different correlations for the Z-factor.

Figure 4B.5. Cut section view of domain showing pressure contour 

This discrepancy is explained by the differences in the dimensions of the static mixers, such as diameter and length, and because the impact of some geometrical features has been ignored in the correlations. Unlike the variability predicted by the existing correlations for the Kenics static mixer, the numerical results are consistent. The value of ΔP obtained from the simulation was 0.565 and the value of  ΔPo obtained using the formula was 0.074. The ratio of ΔP to ΔPo gives us the numerical value of the Z-factor as 7.66.

Figure 4B.6. Cut section view of domain showing velocity contour 

Table 4B.1. Correlation values for Z-factor of Kenics static mixer. Image is taken from the given reference [10] 

A good agreement between the results is also attained for CoV. Overall, This shows that the current model is able to predict the hydrodynamics and mass transfer of brine into water flowing through the various static mixers at considered Reynolds numbers with reasonable accuracy. As a result, the current CFD model is suitable for studying the other static mixers.

Figure  4B.7. Section view of domain showing brine mass fraction contour 

4. Validation and Conclusion

Table 4B.2. Z-factor and CoV calculated by experiments, simulation, and error analysis

Figure 4B.8. CFD results error in Z-factor compared to correlation 

                        Figure 4B.9. CFD results error in CoV compared to experimental 

Table 4B.2 shows the comparison between the experimental and numerical solutions obtained through simulated results. The Z-factor and CoV for the Kenics static mixer have been validated. From Figure 4B.9, it is clear that an increase in the Reynolds number leads to a reduction in CoV.

A three-dimensional CFD simulation was performed for a static mixer, namely, Kenics static mixer, for Reynolds numbers 160. There are a number of correlations in the open literature to predict the pressure drop ratio for well-known static mixers, such as the Kenics one. The correlations have been mainly developed for laminar and turbulent flow regimes. However, most of these correlations do not consider the geometric parameters, such as the blade thickness and the number of elements, and they are therefore under great uncertainty. To validate the numerical model, experimental data available for a static mixer with a baffle shape element and a correlation proposed for the Kenics mixer have been taken as benchmarks. The pressure drop ratio and the coefficient of variation predicted by the CFD model have been compared with experimental data reported in the open literature, and a good agreement (5.27% and 4.86% for Z-factor and CoV, respectively) has been found. Moreover, some available correlations for the Z-factor have been evaluated for a Kenics mixer, and the best fit was obtained for the one listed in Table 4B.1 (Shah and Kale et. al.) for Re = 160.

The velocity field, pressure drop ratio (Z-factor), CoV, and concentration distribution have been determined at various conditions in Kenics static mixer.

Finally, the following remarks may be drawn from this study: 

The current solver can be used to simulate the incompressible mixing of two different fluids with reasonable accuracy in evaluating Z-factor and CoV values for any static mixer. 

C. Transport of gas in a microfluidic channel

1. Introduction

The simplistic case of the transport of gas in different media is modeled to mimic the human gut system. Figure 4C.1 (a) shows a simulation domain that mimics the human gut system indicating the flow of different species. Figure 4C.1 (b) shows the computational domain for assessing the transport of O2 in liquid and air and how the diffusion occurs due to the difference in solubility of O2 in different media. 

Figure 4C.1. (a) Computational domain of a human gut system, (b) computational domain for O2 transport in water and air media

Henry’s law 

Henry's Law is a gas law that states that the amount of gas that is dissolved in a liquid is directly proportional to the partial pressure of that gas. Transport of O2 is carried out from water medium to air, therefore, Henry’s law is important to consider.

KO2 is Henry’s coefficient of O2 = 3.1810-2 (dimensionless) 

This is the ratio of the concentration of the species in the liquid

to the corresponding concentration in the gas.

Along with the momentum equation (Navier-Stokes equation), species

transport equations are solved with Henry’s law of gas solubility, 

which is also taken into account to deduce the actual concentrations of the species. 

2. Numerical analysis of O2 transport 

The numerical simulation includes multiphase modeling where the inlet has a ratio of H2O: O2 mass fraction as 0.9:0.1. O2 is dissolved in water initially, once it comes in contact with air in the domain, the diffusion of O2 starts in gaseous form in the air. The water and dissolved O2 enter the domain with 0.03 mm/s velocity. The partial pressure of O2 is calculated using the given formulation, which is gained through the simulated parameters. 

Figure 4C.2. The partial pressure of O2 in the domain (a) at 0.1 s (b) 2 s and  (c) 5 s

The plot shown below Figure 4C.3 shows the total integrated partial pressure of O2 in the domain at different times. A clear increase in the partial pressure can be observed with an increase in time. This indicates 

that the amount of O2 is continuously increasing as the fluid enters the domain.

Figure 4C.3. The total integrated partial pressure of O2 in the domain with respect to time 

The exact thermophysical properties of water, air, and O2 are taken into account for the numerical simulations. The properties include molecular weight, density, and viscosity The transport is modeled at 300 K and at atmospheric pressure. 

Comparison with Kim et. al. 

The trend of O2 concentration at a particular position in the microfluidic channel is computed and plotted. This was compared with the literature by Kim et. al. (2013) [13]. A good agreement in the trend was seen and shown in Figures 4C.4 (a) and (b). 

Figure 4C.4. (a) O2 mass fraction at varying y position at fixed x=0.003 m and time = 5 s and (b) O2 concentration in the domain of water and PDMS at varying y directions and at x = 0.003 m from the channel entrance and different time 5, 10, 30 and 70 s

D. Mixing in a Split and Recombine Mixture

1. Introduction 

A mixer is one of the basic building blocks in microfluidics, along with components such as  pumps and valves, it is a critical component in several microfluidic devices. For example,  the mixing of reactants is essential for providing homogeneous reaction environments for  chemical and biological reactions. The efficiency of many devices, such as biosensors,  depends on mixing. Microfluidic mixers are thus integral components essential for the  proper functioning of microfluidic devices for a wide range of applications. In this case study, we have performed the computational analysis of the mixing of two micro-fluids through multiphase CFD simulation using a passive mixture called split and recombine mixture. Furthermore, we have calculated the pressure drop obtained in a single mixing step. The results (like the quality of lamination pattern, and pressure drop) obtained are then compared with the results from various CFD commercial software. This case study not only validates the findings but also provides a check on the quality of the results obtained from free source software [14]. 

Mixing is an essential step in many microfluidics processes thus it becomes important to develop methods for more efficient mixing of micro-fluids. This can be done by reducing the mixing time (tmix) which depends on the diffusivity (D) and length scale over which diffusion must act in order to homogenize the concentration, known as the striation length  (lst). The mixing time is then given by the following equation:

(Diffusivity of microfluidics is usually around 10-09 m2/s to 10-11 m2/s)

From the above equation, it can be clearly concluded that to decrease the mixing time we need to decrease the striation length as it will be more convenient than increasing the diffusivity. This can be done in two ways:

• Stretching and folding of the fluid

• By breakup and rejoining of the fluid

2. Methodology

1. The mixing in our split and recombine mixture is modeled through multiphase CFD simulation. 

2. Several Multi-processing methods have been employed for faster numerical calculations.

The end motivation of this validation report is to provide a detailed comparison of the result obtained with commercial CFD software. This is accomplished by generating the volume fraction’s contour of the two fluids on four selected cross-sections of the mixture (A-A, B-B, C-C, D-D) as shown below (Figure 4D.3). 

These contours were obtained as 2D PNG Images which were then imported into a MATLAB script file to calculate the final required comparison-based results.

The MATLAB code first resizes all the input images into a binary image of 256×256 dimensions. These scaled images are then checked for similarity by calculating the percentage of pixels exactly matching in two sets of images that were supposed to be compared.

Geometry and meshing

The geometry development and mesh generation are done successfully. Figures 4D.1 (a, b) and 4D.2 show the model and meshing of the split and recombine mixer, respectively. 

Figure 4D.1. Model of the split and recombine mixer: computational domain of microfluidic channel (a) front view (b) isometric view 

Figure 4D.2. Meshing of the computational domain of microfluidic channel

Excluding 2 inlets, and 1 outlet all other external faces are specified as walls.

Figure 4D.3. Representation of cross sections of the split and recombine mixer (A-A, B-B, C-C, D-D) in the model 

The pure structured mesh is used for the simulation, which is composed of square channels with a side length of d = 500 𝜇m. The structured grid of the mixer consists of 336,000 cells with a quadratic base area with a side length of 25 𝜇m. Thus containing 20 cells on one edge.

3. Results and Discussion

1. Evaluation of pressure drop produced in the mixture as measured in different analyses by respective CFD tools

2. Visual Comparison of results from each aforementioned CFD tool

3. Tabular view of similarity data obtained between the results generated from various commercial CFD software with OpenFOAM

Figure 4D.4. Two microf-fluids with the same inlet velocity enter from the left and are shown in respective colors

Evaluation of pressure drop 

The correctness of the pressure drop across a fluidic structure at a given flow rate, i.e. the fluidic resistance, is one of the central requirements for consistent simulation results. Usually, for a fully developed flow, there exists a linear pressure drop across the length of the channel, (according to Hagen-Poiseuille’s equation given below for circular pipes)

But, for a complex geometry like ours, the pressure drop is hardly predictable because a  fully developed Poiseuille flow is not established and additional pressure losses in constrictions, entrance losses, and turns have to be considered.

Table  4D.1. Pressure drop over the split-and-recombine structure

Visual comparison of results

The qualitative lamination patterns of the mixer were analyzed using a scalar marker with a  numerical diffusion coefficient (ND)of 2 × 10−20 m2 s-1 which is too small to create real physical diffusion as mentioned before (Figure 4D.5).

• Since ND is too small to create any physical diffusion so all the blurriness obtained at the interference of the fluids represents the accuracy of the respective tools as ideally, the interference surface should remain undistorted.

• The scalar field provides also a check of the correctness of the calculated flow field since the streamlines of the flow determine the lamination pattern.

Thus, for the considered problem, identical lamination patterns provided by different simulation tools prove the consistency of the underlying flow fields.

Figure 4D.5. Qualitative comparison of volume fraction of different software with respect to OpenFOAM 

The interphase lines obtained from the simulations are relatively sharp as compared to results from other tools.  This statement can be contoured by the fact that the mesh used in this commercial tools simulation contained two times more cells than the one which was used in our simulation.

The initial lamination (A–A) is (practically)  identical for all codes. After the first twist  (cross-section (B–B)) slight differences are visible. The cross-sections (C–C) and  (D–D) show a corresponding behavior with a slight variation from our simulated results due to the presence of small black corners.

ACE+ and CFX show no significant deviation, while Flow-3D shows a  residual amount of white liquid (𝜙 = 0) in the corners as a result of the used methodology. Fluent shows a slightly different interface region. Overall all tools provide a satisfactory and very similar result. 

Shown below is Table 4D.2 describing the similarity between the various results obtained from our simulation with all the aforementioned CFD tools.

Figure 4D.6. Plot comparing the similarity % of OpenFOAM with other CFD software

Table 4D.2. Comparison of pressure drop calculated from difference simulation software

We can conclude that the results obtained from our simulations resemble most of the results from CFX.

E. The motion of an agitator 

1. Introduction

Impellers, stirrers, or agitators are commonly used in industry for enhancing circulation between multiple fluids, chemicals, species, or phases. Consequently, better designs can be made to minimize dead zones and improve circulation via mechanical means. The same concept is applicable to the bio-processing sector where the enhanced circulation of biomass inside a reactor is important to increase the yields of the gaseous/ liquid products. While some enhancement can be achieved using experimental and empirical means, a thorough computational fluid dynamics (CFD) study can result in even better-optimized designs for circulation. The rotor  inside a tank is simulated and validated against the experimental data and also compared with the commercial software tool, Ansys [15]. 

2. Methodology 

Geometry 

The sample agitator used in this study is a standard Pitched Turbine commonly deployed in industry containing four blades (image attached below, Figure 4E.1). 

Meshing 

The total number of cells in the domain is 568000. The number of cells in the reference work [15] is 600000 cells. 

Unstructured meshing- Unstructured meshing is important to accurately capture the CAD geometry file. An effective mesh is produced with the best mesh efficiency. 

Dynamic meshing- The motion, in this case, rotational motion; of a solid body is captured using a dynamic meshing approach with AMI (Arbitrary Mesh Interface) to accurately capture the agitator motion. 

3. Results and Discussion 

The agitator rotation results inside a tank with water-only phase [15] are validated against the experimental data. Keeping the same number of cells, simulated results are 5.6% more accurate. 

Figure 4E.1. Schematic of the problem setup, velocity contours, unstructured meshing, and the CAD of pitched blade turbine raytraced in ParaView 

4. Results 

The results of generated turbulent kinetic energy are studied. The term K/Utip² is the ratio of turbulent kinetic energy to agitator tip velocity squared. The term r/T denotes the ratio of radial coordinate to tank diameter. Table 4E.1 and the plot in Figure 4E.2 below show the comparison of turbulent kinetic energy values obtained from the simulation and compared with the experimental values and simulation results produced in the reference paper [15]. 

Table 4E.1. Parameter and experimental error comparison

Figure 4E.2. Comparison of simulated results with PIV (Experimental) data

Multiphase interface tracking

The multiphase flows concerning N phases, the volume of fluid (VOF) model developed by Paanduv proves to be an efficient method for interface tracking between N phases. All interpolation schemes, used for flux calculation in advection, operate with the existing MULES-based solvers. They are inherently robust with bounded, conservative numerics that can operate 2nd-order in time and efficiently with large time steps. Numerous simulations were performed and some of them were validated with that of the experimental results and it was found that VOF method solvers gave closer results to that of the experiments thus making it an important method for interface tracking for the multiphase flows. 

Turbulence Modeling
Many turbulence models are present in the literature but the three main models are k-epsilon, k-omega, and k-omega SST. Based on numerous simulations performed by us, k-omega SST is found to be more accurate as compared to the other models. The shear stress transport (SST) formulation combines the best of two worlds. The use of a k-ω formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sublayer, hence the SST k-ω model can be used as a Low-Re turbulence model without any extra damping functions. The SST formulation also switches to a k-ε behavior in the free-stream and thereby avoids the common k-ω problem that the model is too sensitive to the inlet free-stream turbulence properties. 

Multiphase, drag, dynamic motion, and turbulence solvers can be successfully used to perform agitator case studies. As per the validation for a single phase, a good agreement with the experimental data makes the solver applicable for turbomachinery applications. 

F. Effect of surface skewness on biofilm 

1. Introduction
It is well recognized that biofilms are of great importance to human health and welfare, including their potential negative effects on public health (such as increased pathogen resistance to antibiotics) and economics (such as pipeline corrosion and increased friction). Biofilms are also highly beneficial in many engineered systems, such as bioreactors for industrial processes and wastewater treatment. The effects of surface properties on microbial adhesion and biofilm development are therefore of great interest in a variety of environments, and the goals of surface selection/design can be to reduce or enhance biofilm growth and microbial activity [16].

2. Methodology
Reynolds numbers were calculated considering flow between 2 coaxial cylinders (Taylor-Couette flow) using equation (1), where Ω was the angular velocity of the inner cylinder (rad/s), d was the distance between the inner and outer cylinders (m), ro was the radius of the outer cylinder, and ν was the kinematic viscosity of water (m2/s). For comparison with CFD model predictions, shear stress was calculated for a simplified system consisting of two parallel, flat, sliding plates (planar Couette flow) using equation (2),

where = dynamic viscosity 

We simulated velocity gradients, within the negative and positive skew surface features. The meshing is prepared such that the complex geometries and wall shear stress estimations are accurately captured. We have used second-order numerical schemes for the subsequent analysis. Each attachment surface was modeled as a flat top surface moving at 0.83 m/s at a distance of 0.4 cm from a stationary flat bottom surface (the attachment surfaces), which matched the rotational speed and the gap between the inner and outer cylinders in the experiment. The dynamic viscosity of water was set to 8.9x10-4 Pa-s. The flat surface control was used to verify the simulated results with the analytical solution.

The 3-dimensional mesh was effectively simplified to 2 dimensions to reduce computational times by modeling a single cell in the Z-direction, resulting in a control volume containing 7200 cells. Viscous forces were set to solve for a laminar steady state. A no-slip condition was specified at the solid-liquid interfaces. Periodic boundary conditions were set for the X-axis extents to remove any boundary effects: at the minimum and maximum extents of the X-axis (X ±), flow outputs from the X+ boundary were inputs to the X- boundary. The one-cell Z ± boundaries were dictated by symmetry, and the Y ± boundaries were set by the solid surfaces and the no-slip condition. Initial conditions included a stationary fluid and a pressure of 0 Pa.

3. Results and Discussion

CFD modeling was used to better understand how the different topographies may have affected biofilm attachment and growth (Figure 4F.1). As described above, flow conditions were simulated in a simplified system consisting of 2 flat parallel plates. The focus of these analyses was on local shear stresses, as these provide an indication of velocity gradients, which in turn drive mass transport (both for transport of cells to the surfaces during initial adhesion and transport of metabolites to and from attached biofilms).

The CFD simulations predicted the same uniform shear stress on the moving flat surface (0.18 Pa) as did the analytical solution for concentric cylinders (Equation (2)), which provided some confirmation of the CFD model. It is therefore hypothesized that the predicted greater area with high shear values on the negative skew surface relative to the positive skew surface was indicative of increased mass transport to the surface, which would explain the greater rates of biofilm attachment, growth, and activity on the negative skew surface relative to the other surfaces. 

Figure 4F1. Plan view schematic of an annular bioreactor with an inner cylinder (radius ri) rotating with annular speed (𝛀) inside of the stationary cylinder (radius ro), with 3 different surfaces attached to the inner cylinder shown in detail to the right.

The CFD results, therefore, provide a possible explanation for the higher rates of biofilm development and activity found on the negative skewness surfaces than on the positive skewness and flat surfaces. In contrast, the biomass on the negative skewness surface was concentrated along the rounded peaks, obscuring the attachment surface and appearing denser than the biofilm on the positive skewness surface. The troughs were clearly visible, indicating little or no growth, corresponding to the areas with the predicted lowest shear values (Figure 4F.2 (c)). These observations were consistent with the hypothesis suggested above that the higher shear stress along the rounded peaks of the negative skewness surfaces was preferred attachment sites and/or they provided zones of increased mass transfer to accelerate biofilm growth. In contrast, the high shear zones predicted on the positive skew surface (near the sharp peaks, Figure 4F.2 (a) and (c)) were not obviously associated with increased attached biomass. It is possible that the very high shear rates over a small area near the positive skew peaks (the maximum shear of 8.5 Pa was over twice the maximum shear of 4.1 Pa on the negative shear surface) negatively affected biomass adhesion and biofilm development because of the shear forces, which may have overwhelmed any benefits associated with increased mass transfer rates.

A. The rise of the bubble in the channel 

1. Introduction

The rise of bubbles is a critical phenomenon in many industries such as metallurgical, nuclear, microfluidic, wastewater, etc. The problem statement here presents a validation case showing the comparison of simulated terminal bubble shape and rising velocity with the experimental data for four different mesh elements [17]. 

2. CAD and Meshing

Figure 5A.1. (a) CAD with the bubble, (b) Mesh with the bubble

3. Simulation Detail

Table 5A.1. Properties of fluid used in the simulation

Both liquid and bubble are assumed to be stationary at the initial state. The mesh used is uniform and structured with the number of cells 3950, 15400, 34350, and 60800 respectively.

4. Results and Discussion

Figure 5A2. Comparison of the bubble shape with four different mesh cells for ρ = 1384 kgm⁻³, ν = 0.000151 m²s⁻¹

Figure 5A.3. (a) The plot of bubble rising velocity with time for various mesh cells, (b) Clock time taken by traditional VOF method for different mesh cells

Figure 5A.2 shows the bubble profile for four mesh grids viz. 3950, 15400, 34350, and 60800. The bubble shape is almost the same as the experiments [17] with a slight change, the tail or skirt portion of the bubble is slightly elongated. As seen in Figure 5A.3 (b), the rising velocity is compared with the experimental data (black color) using four mesh grids viz. 3950, 15400, 34350, and 60800 respectively. The plot clearly shows that a finer mesh grid gives closer simulation results when compared with the experiment [17]. This clearly proves that to get accurate results, mesh refinement is important. The values provided in Table 5A.2 show an increase in error in rising velocity with reduced mesh cells. However, increasing the mesh cells will certainly increase the computational clock time thus reducing the performance as shown in Table 5A.3. 

Table 5A.2. Error in rising velocity percentages for different mesh cells

Table 5A.3. Clock time in minutes for different mesh cells

B. The flow of the water from a vertical step 

1. Introduction

The vertical water drop is commonly used in both natural and artificial channels such as irrigation systems, weir flow, spillways, and sustainable hydraulic structures. The simulation performed here shows the validation and verification of the simulated results with that of the experimental results performed by Rajaratnam et. al. [18]. Basically, the pool dimensions and the velocity profiles are compared.

Figure 5B.1. Schematic sketch of the drop

2. CAD and Meshing

             (a)                                         (b)   

Figure 5B2. (a) CAD with the bubble, (b) Mesh with the bubble

3. Simulation Detail

The density and kinematic viscosity of water is taken as 1000 kg/m³ and 1e-6 m²/s respectively. The velocity of water flow from the inlet is 0.84 m/s which resembles the discharge rate used in the experiments.

4. Results and Discussion

                                         (a)                                                                                     (b)

                                  (a)                                                                                                (b)

Figure 5B.3. Water flowing through the water flow step at a time (a) 5 s, (b) 10 s, (c) 15 s, and (d) 20 s

As observed in Figures 5B.4 (a) and (b), the plot of inlet velocity varying along with the inlet distance and edge velocity varying along the edge distance are in good agreement with the experimental results. The overall error percentage for the velocity comes to around 4.5%. The overall error of pool dimensions which includes the pool length and height is 5.6%.

                            (a)                                                                                           (b)

Figure 5B.4. The plot of (a) inlet velocity vs distance above the inlet, (b) edge velocity vs distance above the edge

Table 5B.1. Error percentages for base velocity and pool dimensions

C. Sloshing in a cuboidal container

1. Introduction

Sloshing is a phenomenon commonly found in partially filled liquid storage vessels, such as liquid cargo tanks and fuel tanks, in motion. Violent sloshing may result in serious structural damage to the liquid tank or even overturn the liquid cargo ship. Thus, a reliable prediction of the sloshing is crucial for the design and deployment of such structures. For this purpose, theoretical analyses, physical experiments, and numerical simulations are commonly utilized.

To study the sloshing phenomenon inside a cuboidal container due to a sinusoidal motion. The visualization results and pressure plots at probe locations are compared and validated with experimental data [19].

2. CAD and Meshing

The dimensions of the container are 600 mm × 300 mm × 650 mm, with the size of the mesh kept at 8 mm uniform along all three directions. The total mesh cell count is 230850.

Figure 5C.1. Schematic of the 3D domain with locations of pressure probes

                              (a)                                                                                                  (b)

                                                                                                                      Figure 5C.2. (a) CAD and, (b) Mesh of the simulation domain

3. Simulation Detail

An external sinusoidal motion x = -A sin(ωt) is applied to move the domain. The amplitude and frequency of the sinusoidal motion are 7 mm and 4.749 rad/s. The initial height of the water is kept at 90 mm.

4. Results and Discussion

Figure 5C.3 shows the sloshing behavior comparison between the experiments and simulations at time instant  (a) & (b) t = 8.55 s, (c) & (d) t = 8.7 s, (e) & (f) t = 8.92 s and (g) & (h) t = 9.17 s respectively. The comparison simulation images are in good mutual agreement with the experimental setup.

The pressure plot at probe location P1 as shown in Figure 5C.4, exhibits a sinusoidal-type motion when water comes into and fro contact with the probe.

                   (a)                                                                                         (b)

                  (a)                                                                                                 (b)

                (a)                                                                                        (b)

Figure 5C.3. Visual comparison of simulation with experiments at (a) & (b) t=8.55 s, (c) & (d) t=8.7 s, (e) & (f) t=8.92 s and (g) & (h) t=9.17 s

D. Transport of sediment and scouring near a bridge pier

1. Introduction

Scouring of sediments is found to be a major problem in bridge piers, abutments, embankments in rivers, sea, etc. Basically, the simulation performed is in the case of sediment scour of a bridge pier. The sediment scour depth for a streamlined bridge pier is validated with the experimental results [20].

2. CAD and Meshing

Figure 5D.1. (a) CAD and, (b) Mesh of the simulation domain

3. Simulation Detail

The velocity of the water flow is kept at 0.3 m/s. Figure 5D.2 represents the simulation domain containing sediment (red color) and water (blue color). The flow of water at a velocity of 0.3 m/s is from left to right.

Figure 5D.1. (a) CAD and, (b) Mesh of the simulation domain

3. Simulation Detail

The velocity of the water flow is kept at 0.3 m/s. Figure 5D.2 represents the simulation domain containing sediment (red color) and water (blue color). The flow of water at a velocity of 0.3 m/s is from left to right.

Figure 5D.3. Experimental image of scouring in a streamlined bridge pier at velocity 0.3 m/s

Figure 5D.4. Simulation image showing the sediment transport velocity along with the scouring happening close to the bridge pier

From experimental and simulation results shown in Figures 5D.3 and 5D.4 clearly illustrate the scouring phenomenon occurring close to the streamline bridge pier structure. The sediment transport velocity is maximum at the sideways of the piers which means maximum scouring occurs sideways. Also at the front of the pier scouring depth is visible. 

Table 5D.1. Comparison of scouring depth for both experiment and simulation setup

E. Resistance forces in a Duisburg Test Case (DTC) hull

1. Introduction

In this case study, we have calculated the total resistance force which is a sum of pressure and viscous forces on the DTC hull, and compared it with experimental results [21].

2. CAD and Meshing

                              (a)                                                                                         (b)

Figure 5E.1. (a) Side view and (b) Isometric view of the DTC hull

3. Simulation Detail

The velocity of the water which hits the hull is kept at 1.668 m/s and k-omega turbulence modeling is used to capture the turbulence phenomenon.

4. Result and Discussion

                                 (a)                                                                                              (b)

Figure 5E.2. (a) DTC hull over water (b) Velocity behavior of water close to the DTC hull

Figure 5E.2 (b) shows the velocity of the water surrounding the DTC hull. The velocity is reduced when water is close to the DTC hull. Table 5E.1 clearly shows the comparison between the experimental and simulation results for the resistance force and the values are in very good agreement.

Table 5E.1. Comparison of scouring depth for both experiment and simulation setup

F. Flow analysis in a pipe (Hagen-Poiseuille flow)

1. Introduction

A simple flow through a pipe is modeled in this case study. A simple Hagen-Poiseuille validation for the laminar and turbulent case condition for a cylindrical pipe flow is accomplished.

2. CAD and Meshing

                             (a)                                                                                                   (b)

Figure 5F.1. (a) CAD of the pipe, (b) Mesh of the pipe

3. Simulation Detail

To achieve that we take a 2 m long pipe with a diameter of 0.05 m. The Reynolds number is 1000 which falls in the laminar flow regime and for the turbulent regime, the Reynolds number is taken to be 200000.

The simulation performed is henceforth validated with Hagen Poiseuille’s equations. Equations (2) and (3) depict the velocity profile and pressure drop equations respectively.

                     (a)                                                                                                       (b)

Figure 5F.2. (a) Radial Velocity map and (b) radial velocity plot at an axial distance of 1.9 m of the pipe

Figure 5F.2 (a) shows the radial velocity profile is maximum at the center of the pipe and decreases as it moves toward the pipe wall. Figure 5F.2 (b) shows the traditional parabolic profile of laminar flow in a pipe with a maximum velocity of 0.2 m/s. The simulated result (red line) is validated with the analytical result (blue line) of Equation (2).

                                   (a)                                                                                             (b)

A. Supersonic flow around of a bullet 

1. Introduction

Firearms industry is one of the areas of application in CFD engineering for external aerodynamic flow over bullets and ballistics. The purpose of this case is to present the simulation of a Berger Bullet (USA) 0.308 caliber 155-grain VLD, which is subjected to supersonic (greater than 1.2 Mach) speeds. This information will be useful to understand the stability of the warhead during its flight trajectory (external ballistics). Calculated drag data is compared with experimental data to validate our case setup. Once we understand and validate this case we can try simulating CFD code by replacing the bullet body with a different body and thus, this tutorial demonstration can prove useful in a wide variety of external aerodynamic problems.

The bullet geometry is symmetric in a cylindrical coordinate system. Axisymmetric wedge problems require less computation time compared to symmetric problems. It is suggested that the wedge angle should be (< 5 degrees). We use this default wedge angle value of 5 degrees to solve our wedge problem. If we try to reduce the wedge angle below 5 degrees we will end up having very small cells near the axis of the wedge which may lead to high non-orthogonality, skewness, and aspect ratio, in short, bad cells. If we increase the wedge angle above 5 degrees we are spending extra time and resources by just solving larger size cells in the wedge direction. 5-degree wedge fluid domain gives us accurate results that are conformal with the experimental values as showcased further in this case study.

2. Simulation details

Geometry- We received the dimensions of the Berger bullet and .stl file is constructed [22]. The wake region length is ~(twice) the total length of the bullet. Dimensions are in meters. The geometry is symmetric about the Y and Z axis. That means the X-axis lies exactly along the center axis of this cylindrically symmetric geometry.

Figure 6A.1. CAD of the bullet

Meshing- We generated a hexahedral mesh with boundary layers to validate our CFD problem. 

Figure 6A.2. (a) CAD of the bullet and the domain around the bullet (b) CAD of the bullet and background meshing of the external domain

Figure 6A.3. Meshing of the entire computational domain (a) & (b) isometric and front view (c) & (d) mesh quality at the corners 

In Figure 6A.4 (b), it is clearly visible that the smallest edge of the bullet has been divided into 10 parts and this is the criteria for fine mesh in this case. For coarse and medium mesh this smallest edge has been divided into 2 and 5 parts respectively. 

Figure 6A.4. Snippets of the meshing (a) isometric view of the external domain excluding the bullet (b) & (c) closer look at the corners  (d) 2D view of meshing of the computational domain 

3. Results and Discussion

The wave drag trend in the supersonic flow of the bullet is captured well in the simulations and the fine mesh results are reasonably accurate (~ 3.5%) in drag prediction. Calculations with available data interpolation suggest that super fine mesh results can provide (~2%) accuracy and is possible to achieve using our approach.

Figure 6A.5. Pressure distribution around the bullet                                                          Figure 6A.5. Velocity profile around the bullet

Fig 6A.6. Temperature profile around the bullet                                              Fig 6A.7. Drag coefficient variation considering simulation time

The above results show values for Mach no. 2.5. It can be observed that the accuracy of drag coefficient prediction is increasing as we increase the refinement of the mesh. Accurate simulations up to ~96.5% can be simulated using a fine mesh with boundary layers. The accuracy of the solution can be improved by creating super fine mesh by increasing the number of divisions for the smallest edge of the bullet to 15 or 20 along with better computational resources.

Table 6A.1. Details of the parameters, error deduced from the simulation for three mesh count 

B. Flow around a drone propeller

1. Introduction

The objective of this validation is to compare the experimental [23] and numerical results of the coefficient of the thrust of a drone propeller. 

2. CAD and Meshing

                         (a)                                                                                               (b)

Figure 6B.1. (a) CAD of the APC slow flyer 10”x7”, (b) Mesh with the Arbitrary Mesh Interface (AMI)

3. Simulation Detail

The simulation used AMI (Arbitrary Mesh Interface) for the dynamic meshing which rotates the propeller, kept inside the AMI. The total mesh cell count was 4,12,278. In the following scenario, we are posed with a problem consisting of initial boundary conditions of a range of free-stream velocities and constant turbulence intensity (0.1%). Outlet conditions were set as inlet-outlet type with uniform pressure and velocity.

Table 6B.1. Initial free stream velocity and rotational velocity

4. Results and Discussion

Table 6B.2. Error percentages of coefficient of thrust 

Table 6B.2 and Figure 6B.2 shows the error percentages for three different cases with three different free stream velocities. The error percentage of the coefficient of thrust was calculated. Thus simulated results showed less percentage of errors. The thrust coefficient decreases with an increase in velocity as observed in Figure 6B.3. 

    Fig 6B.2. Histogram showing the error for the simulated results                     Fig 6B.3. Plot showing the coefficient of thrust vs velocity 

                                                       Figure 6B.4. (a) Pressure and (b) Velocity visualization in the propeller

From the above results, we can conclude that the simulated results showed good accuracy and the variation of thrust vs free stream velocity can be understood very well.

C. Aerodynamics around a blunt body

1. Introduction

The forebody of missiles, rockets, reentry capsules, etc., generally adopt blunt shapes with a motive to withstand the severe heating experienced during high-speed flight regimes and as well to optimize the availability of volume to accommodate larger payloads. At high speeds (supersonic/hypersonic), the formation of a detached shock wave leads to higher aerodynamic drag and heating, leading to either loss in performance or erosion/ablation of the surface. The use of gas jets, energy deposition techniques, breathing noses, spike-tipped noses, etc., are some of the methods which are adopted to suppress these problems. Due to its simplicity, flow over the spiked body has been studied in detail in the last few decades. The flow field around a blunt body and its modification due to the presence of a spike having either a sharp or blunt or flat aero-disk is schematically shown in Figure 6C.2. The presence of the spike leads to the formation of a spike shock, separated zone, and separation shock due to the adverse pressure gradient, recirculation zone, shear layer, etc., which will depend on the shape of the tip of the spike. There exists the possibility of the formation of separation shock which will impinge on the main body resulting in reattachment shock. The modification of the flow field over the main body due to the presence of a spike leads to a reduction in drag and heating. The main reason for the reduction is due to the alteration of strong shock on the main body with a spike shock. A major part of the main body will experience lesser pressure due to the presence of a recirculation zone and hence a reduction in drag. The flow field between the tip of the spike and the main body will strongly depend on incoming freestream flow conditions, the shape of the tip of the spike, the length, and diameter of the stem connecting the spike tip to the body, etc. The reduction of drag/heating will depend upon the shape at the tip of the spike. [24]

It is established that the use of a spike does reduce drag, and the amount of reduction depends on the shape and size of the spike. Furthermore, the formation of a recirculation zone in between the spike and body might even lead to unsteadiness in the flow, and hence its knowledge is of importance.

2. CAD Details

The geometry adopted in the present study is a blunt hemispherical body having a base diameter (D) of 15 mm and an overall length of 1.5D. The aerospike used has a basic stem diameter of 2 mm (0.133D). Flat head aero-disk (Aerospike) having a diameter of 0.33D with a flare angle of 130 degrees. Spike length was kept at the range of L/D = 1.

Figure 6C.1. Blunt body CAD model actual and .stl file with domain 

Figure 6C.2. Aerospike CAD model actual and .stl file with domain

3. Mesh Details

To capture the unsteady flow field, only 2D laminar computation has been made. This was done to reduce the computation time and as well the expected uncertainty in the suitability of a turbulence model to capture an unsteady flow field. For unsteady flow computation, the grid used has to be different from the grid used for steady computation. Very fine grids near the surface had to be used, The overall size of the grid is 20,500 and 22,500 elements respectively. 

Figure 6C.3. Meshed model blunt body                                                                             Figure 6C.4. Meshed model blunt body with aerospike 

SnappyHexMesh tool was used in parallel mode for the grid generation. The meshing was carried out on 4 cores of an 8 GB system. The mesh generated by snappyHexMesh is extruded to create an axisymmetric wedge domain.

4. Simulation Details

All the experiments were performed using the blowdown type Supersonic Wind Tunnel at Birla Institute of Technology, Mesra, Ranchi, India having a test section size of 50 mm X 100 mm. The present series of experiments have been made at a fixed Mach number of 2.0, with

settling chamber pressure of about 3.2x105 N/m2, which was measured using a pressure transducer (Make Sensym, Model ASCX150DN), and Reynolds number (Re) of 3.5 x 105 based on the base diameter of the model.

The 2D simulation was performed, which uses a finite volume approach to solve compressible Reynolds-Averaged Navier-Stokes equations. For the unsteady state solution, axisymmetric computations have been made adopting an explicit coupled solver and using the "k-𝝐" turbulence model. The use of this turbulence model has been arrived at, after making necessary grid sensitivity tests, convergence history, and obtaining a good comparison with the experimental and numerical results reported in the literature.

The supersonic freestream boundary condition at the inlet was specified. inletOutlet boundary condition was imposed on the outlet to avoid backflow. No slip wall boundary conditions with suitable near-wall treatment for turbulent flows were enforced on the hemispherical body and spike. The overall grid, computational domain, and boundary conditions in addition to the convergence history for drag were also monitored during the entire solution period. Results were analyzed only when it was ascertained that the residuals had converged to the order of 10-5.

Computations were continued with very small time steps (Δt= 10-8s) to capture the unsteadiness of flow. During simulation, specific points were monitored, and pressure, temperature, and drag histories at those points were recorded for further processing of data. 

Table 6C.1. Important parameters used in the simulation

5. Results and Discussion

These results indicate the occurrence of different types of flow features with and without attachment of aerospike. The pressure contours for both cases are shown in Figure 6C.8 along with some representative values of pressures. Pressure is observed to be maximum on the main body in the vicinity of the reattachment location for all the cases. For sharp spikes, the values are higher than the blunt and aerospike of the same length. [24]

All the results presented till now (density, velocity, pressure contours, etc.) clearly indicate the existence of a variety of separated zones on the main body due to the presence of a spike. Therefore, the possibility of unsteadiness in the flow field could also exist due to these complex flows. Very limited information is available about the unsteadiness of the flow field around the hemispherical main body with a spike attached to it. Only laminar simulations of the flow have been made to capture the unsteady flow field, with a reason to reduce the computation time and avoid the uncertainty of the turbulence model.

Figure 6C5. Schematic of the flow field around a blunt body without and with spike [24]

Figure 6C.6. Comparison of the velocity flow field around a blunt body without and with a spike with experimental images [24]

Figure 6C.7. Comparison of the temperature flow field around a blunt body without and with spike

Figures 6C.6, 6C.7, and 6C.8 show velocity, temperature, and pressure distribution, respectively. For velocity, a recirculation low-velocity zone can be seen forming in between the tip of the spike and the blunt body. In temperature contours, it is seen that the temperature is distributed over a region with an aerospike attached in comparison to a blunt body. Pressure contours show that the high-pressure region is dispersed when an aerospike is attached compared to the high-pressure region forming on the nose of a blunt body without an aerospike. This protects the blunt body from experiencing direct shock due to supersonic speed.

Figure 6C.8. Comparison of the pressure flow field around a blunt body without and with spike

6. Validation and Conclusion

Experiments and computations have been made to obtain the flow field on a main hemispherical body at a supersonic Mach number of 2, in the presence of a sharp spike of different lengths and blunt and aerospike for L/D=1. The drag changes considerably with the change in spike shape and as well with a change in spike length for the sharp spike [24].

The computed results compare well with corresponding experiments. The maximum drag reduction of around 46.80 % is recorded for aerospike through experiments and 46.22 % through simulations, which indicates the importance of spike tip shape on the overall flow field. Drag reduction up to 45.40 % is obtained with the change in spike length to L/D=1.50 for a sharp spike showing decreasing drag with the increasing length of the spike [24]. 

Table 6C.2. Validation and comparison results for drag coefficient value

A. Flow around a DrivAer

1. Introduction

Nowadays with the help of computational fluid dynamics (CFD) and high-performance computing (HPC) technology, vehicle aerodynamics engineers are reducing the wind tunnel experiments and accelerating the vehicle design cycle. A lot of the investigations in automotive aerodynamics are still based on strongly simplified generic bodies such as the Ahmed Body [25]. These simple car models like Ahmed Body or SAE body make it easy to relate the observed phenomena to specific areas and thus help to understand basic flow structures. At the same time, more complex flow phenomena, e.g. at the underbody, wheels/wheelhouses, and around the rearview mirrors, etc., cannot be reproduced due to the oversimplification of these geometries. On the other hand, it is usually not feasible to investigate these phenomena on a specific production vehicle, as, due to its short life span and restricted access, typically little validation data is available. Recognizing the need for a model combining the strengths of both approaches, various more or less generic models, such as the VW reference car and the MIRA reference car, have been proposed. However, while these reference cars mark a step in the right direction, these models are still too generic to completely understand the complex phenomena occurring in realistic vehicles.

To close this gap, the Institute of Aerodynamics and Fluid Mechanics of the Technische Universitat Munchen (TUM), in cooperation with two major car companies, the Audi AG and the BMW Group, has proposed a new realistic generic car model called “DrivAer Model”. The body is based on two typical medium-class vehicles (Audi A4 and BMW3 series) and includes three interchangeable tops and two different underbody geometries to allow for a high universality. To encourage the use of the DrivAer model in independent research projects, the TUM research group is open to sharing the geometry, and a comprehensive experimental database is published in different papers.

The aim of the current work is to use the DrivAer body model and its experimental results to validate our CFD solver. 

2. CAD Details

This current study is focused on one DrivAer body namely Fastback and smooth underbody type as shown in Figure 7A.1. Hence, this model is simulated at ground conditions without ground effect (WoGS). Different parts of the Fastback model like mirrors and wheels are also considered for the simulation study as shown in Figure 7A.1. The individual drag contribution of these components is also calculated along with the total car drag. The vehicle velocity considered for this numerical study is 40 m/s, the Reynolds number is 4.87⨯106, and the turbulence model used is k-w-SST.

To capture the flow around the car model the computational domain is created around the car which is often called a numerical wind tunnel as shown in Figure 7A.2. This computational domain is like a physical wind tunnel test section. The length of the computational domain is 28 m, its width is 10 m and its height is 7 m. While the length of the physical wind tunnel test section is 4.8 m, the width is 2.4 m and the height is 1.8 m. The blockage ratio used in the computational domain is 0.57% while in the physical wind tunnel, it was 6.17%.

Figure 7A1. CAD model isometric view in .stl format 

Figure 7A.2. CAD model with domain isometric view

Figure 7A.2. CAD model with domain isometric view

3. Mesh Details

To capture flow physics around the car more accurately the mesh size kept near the car is fine and becomes coarse when going away from the car as shown in Figure 7A.4. To achieve this variable mesh size distribution various refinement boxes are created with different dimensions. To capture the boundary layer around the car, multiple layers of very fine and fine element sizes are kept around the full car and along the road which is shown in Figure 7A.4. The y+ obtained is 30 which is supposed to be good for incompressible simulations.

Figure 7A.3. Meshed model isometric view using snappyHexMesh

Figure 7A.4. A meshed model with domain isometric view using snappyHexMesh

The hexahedral meshing was generated on 4 cores of an 8 GB system. The mesh generated on smooth underbody models is about 2.08 million volume cells.

4. Simulation Details

Boundary conditions are implemented on the computational domain, car body, and road. Velocity was specified at the inlet and pressure was specified at the outlet. The road was given a translation speed of 40 m/s in the direction of flow in simulations without ground simulations. Car wheels are considered stationary without ground simulation. All other car parts are defined as stationary walls.

Table 7A.1.  Important parameters used in the simulation

The model considered here is with wheels and mirrors and was simulated for 500 iterations. The convergence achieved was 5 decades fall for flow variables P and U and 6 decades fall for turbulence quantities omega and k. Convergence of drag force is also monitored and achieved up to 15 counts in the last 100 iterations.

5. Results and Discussion

The total drag values obtained for the fastback model with a smooth underbody without ground effect are given in Table 7A.2. Table 7A.2 shows the drag values obtained for the fastback model with a smooth underbody considering without ground effect. The percentage error found in the simulated drag values is in the range of 1.18% as compared to the experimental drag values. The simulated drag value is taken as the average drag value of the last 100 iterations.

Pressure and velocity contours have also been plotted as shown in Figures 7A.5 and 7A.6, respectively.

        Figure 7A.5. Pressure contour displayed on car                                           Figure 7A.6. Section plane of velocity contour

Back vortices and side vortices which are responsible for drag generation are also captured as shown in Figures 7A.7 and 7A.8, respectively.

Figure 7A.7. Section plane of velocity showing back vortices

Figure 7A.8. Velocity streamlines showing side vortices

6. Validation and Conclusion

This current study is focused on one model of the DrivAer body namely FastBack and smooth underbody. Hence, one configuration is simulated without ground effect conditions. The coefficients of drag (Cd) values are within the 1.18% error band as compared to the experimental values published by the TUM. 

Table 7A.2. Validation results and error percentage