The subject of multiphase flows encompasses many processes in nature and a broad range of engineering applications, such as weather forecasting, fuel injection, sprays, and spreading of substances in agriculture. To investigate these processes the Institute of Aerospace Thermodynamics (ITLR) uses the direct numerical simulation (DNS) in-house code Free Surface 3D (FS3D). The code is continuously optimized and expanded with new features and has been in use for more than 20 years.
The program FS3D was specially developed to compute the incompressible Navier-Stokes equations as well as the energy equation, with free surfaces. Complex phenomena demanding strong computational effort can be simulated because the code works on massive parallel architectures. Due to DNS, and thus resolving the smallest temporal and spatial scales, no turbulence modeling is needed. In the last years a vast number of investigations were performed with FS3D: for instance, phase transitions like freezing and evaporation, basic drop and bubble dynamics processes, droplet impacts on a thin film (“splashing”), and primary jet breakup, as well as spray simulations, studies involving multiple components, wave breaking processes, and many more.
The flow field is computed by solving the conservation equations of mass, momentum, and energy in a one-field formulation on a Cartesian grid using finite volumes. The different fluids and phases are treated as a single fluid with variable thermophysical properties that change across the interface. Based on the used Volume-of-Fluid (VOF) method additional indicator variables are used to identify different phases. The VOF variables 𝒇i are defined as and represent the different phases liquid (i=1), vapour (i=2) and solid (i=3). To ensure a successful advection of the VOF variable, a sharp interface, as well as its exact position, is required. This is done using the piecewise linear interface reconstruction (PLIC) method, which reconstructs a plane on a geometrical basis and, therefore, can determine the liquid and gaseous fluxes across the cell faces. The advection can be achieved with second-order accuracy by using two different methods . For the computation of the surface tension several models are implemented in FS3D; for instance, the conservative continuous surface stress model (CSS), the continuum surface force model (CSF) or a balanced force approach (CSFb), which allows a significant reduction of parasitic currents. Due to the volume conservation in incompressible flow Poisson’s equation of pressure needs to be solved, which is achieved by using a multigrid solver. In order to perform simulations with high spatial resolutions, FS3D is fully parallelized using MPI and OpenMP. This makes it possible to perform simulations with more than a billion cells on the supercomputer Cray-XC40 at HLRS. Some applications of FS3D and results are presented in the following.
Supercooled water droplets exist in liquid form at temperatures below the freezing point. They are present in atmospheric clouds at high altitude and are important for phenomena like rain, snow, and hail. The understanding of the freezing process, its parametrization, and the link to a macrophysical system such as a whole cloud is essential for the development of meteorological models.
The diameter of a typical supercooled droplet, as it exists in clouds, is on the order of 100 μm whereas the ice nucleus is in the nanometer range. This large difference in the scales requires a fine resolution of the computational grid. To capture the complex anisotropic structures that develop as the supercooled droplet solidifies, an anisotropic surface energy density is considered at the solid-liquid boundary using the Gibbs-Thomson equation. The energy equation is solved implicitly in a two-field formulation in order to remove the severe timestep constraints of solidification processes. The density of both ice and water are considered equal. This is a reasonable assumption and greatly simplifies the problem at hand. A typical setup consists of a computational grid with 512 × 512 × 512 cells where the initial nucleus is resolved by roughly 20 cells. A visualization of a hexagonally growing ice particle embedded in a supercooled water droplet is shown in Fig. 1.
Evaporation of supercooled water droplets
Not only freezing processes but also the evaporation of supercooled water droplets need to be understood for the improvement of meteorological models. In the presented study the evaporation rate, depending on the relative humidity of the ambient air, is in the focus of numerical investigations with FS3D.
Several simulations of levitated supercooled water droplets are performed at different constant ambient temperatures and varying relative humidities Φ, with one example shown in Fig. 2. The evaporation rate β is determined and compared to experimental measurements . The setup consists of an inflow boundary on the left side, an outflow boundary on the right side, and free slip conditions on all lateral boundaries. The grid resolution is 512 × 256 × 256 cells and the diameter of the spherical droplet is resolved by approximately 26 cells.
The resulting dependency of the evaporation rate on the relative humidity is depicted in Fig. 3, for an ambient temperature of T∞=268,15 K. The numerical results agree very well with experimental data. This shows that FS3D is capable of simulating the evaporation of supercooled water droplets and therefore can help to improve models for weather forecast. For example, future numerical simulations of the evaporation of several supercooled water droplets and their interaction could be investigated, a goal that is currently not feasible experimentally.
Non-newtonian jet break up
Liquid jet break up is a process in which a fluid stream is injected into a surrounding medium and disintegrates into many smaller droplets. It appears in many technical applications; for instance, fuel injection in combustion gas turbines, water jets for firefighting, spray painting, spray drying, or ink jet printing. In some of these cases an additional level of complexity is introduced if the injected liquids are non-Newtonian; i.e., they have a shear dependent viscosity. Due to the complex physical processes, which happen on very small scales in space and time, it is hard to capture jet break up by experimental methods in great detail. For this reason it is a major subject for numerical investigations, and therefore, for investigations with FS3D.
We are simulating the injection of aqueous solutions of the polymer Praestol into ambient air. The shear-thinning behavior is incorporated by using the Carreau-Yasuda model. The largest simulations are done on a 2304 × 768 × 768 grid, using over 1.3 billion cells, where the cells in the main jet region have an edge length of 4∙10-5 m . The simulated real time is in the order of 10 ms.
We investigate the influence of different destabilizing parameters on the jet (see Fig. 4), such as the Reynolds number, the velocity profile at the nozzle or the concentration of the injected solutions (and therefore the severity of the non-Newtonian properties). We analyze the influence of these parameters on the jet break up behavior, quantified by the liquid surface area, the surface waves disturbing the jet surface and the droplet size distribution . We then investigate the three-dimensional simulation data, such as the velocity field or the internal viscosity distribution, in detail to explain the differences in jet behavior (see Fig. 5).
The interaction between an airflow and a water surface influences many environmental processes. This is particularly important for the formation and amplification of hurricanes. Water waves, wave breaking processes, and entrained water droplets play a crucial role in the momentum, energy, and mass transfer in the atmospheric boundary layer.
In order to simulate a wind wave from scratch a quiescent water layer with a flat surface and an air layer with a constant velocity field are initialized. The computational domain, corresponding to one wavelength of λ=5 cm, has a resolution of 512 × 256 × 1024 cells. Every simulation is performed on the Cray-XC40 at HLRS with at least several thousand processors. Due to transition, the air interacts with the water surface and a wind wave develops, shown in Fig. 6. In the first step the occurring parasitic capillary waves on the frontside of the wind wave are evaluated. Wave steepnesses and the different wave lengths of all parasitic capillary waves offer detailed insights into energy dissipation mechanisms, which could not be gained from experiments. In a second step the wind is enhanced by applying a wind stress boundary condition at the top of the computational domain. This leads to the growth of the wave amplitude and finally to wave breaking. Not only phenomenological comparison of this process with experiments, but also information about temporal evolution of the wave energy, structures in the water layer, or dynamics of vortices are remarkable results of these simulations. For future investigations of wind waves and, for example, droplet entrainment from the water surface higher velocities, higher resolutions, and therefore, higher computational power will be needed. Such simulations requiring more than one billion cells makes the use of supercomputers indispensable.
If a liquid droplet impacts on a thin wall film, the resulting phenomena can be very complex. Impact velocity, droplet size and wall film thickness have a large influence on the shape and morphology of the observed crown. If the conditions are such that secondary droplets are ejected, this phenomenon is called splashing.
The splashing process is highly unsteady and its appearance is dominated by occurring instabilities that have a wide range of different scales. However, only a limited amount of properties are accessible through experiments. For example, thickness of the crown wall and velocity profiles are difficult to obtain experimentally.
Currently, we are able to perform simulations with up to one billion cells. A rendering of an exemplary simulation is shown in Fig. 7. In order to capture splashing processes on the smallest scale, a very high resolution is required. Therefore, often only a quarter of the physical domain is simulated by applying symmetry boundary conditions.
When the droplet and the wall film consist of two different liquids, additional phenomena occur that cannot be explained anymore with single-component splashing theories. One reason for this is that not only the properties of the liquids themselves but also their ratio matters.
Due to this, a multi-component module is implemented in FS3D, which captures the concentration distribution of each component within the liquid phase. This makes it possible to evaluate, for example, composition of the secondary droplets. One technical application for which this is important is the interaction of fuel droplets with the lubricating oil film on the cylinder in a diesel engine. This interaction occurs during the regeneration of the particle filter and leads to both a dilution of the engine oil wall film and to higher pollutant emissions. Here, a better understanding of two-component splashing dynamics can be a great advantage in order to minimize both engine emissions and lubrication losses.
The FS3D team gratefully acknowledges support by the High Performance Computing Center Stuttgart over all the years. In addition we kindly acknowledge the financial support by the Deutsche Forschungsgemeinschaft (DFG) in the projects SFB-TRR75, WE2549/35-1, and SimTech.
 Eisenschmidt, K., Ertl, M., Gomaa, H., Kieffer-Roth, C., Meister, C., Rauschenberger, P., Reitzle, M., Schlottke, K., Weigand, B.:
Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation, 272, pp. 508-517, 2016.
 Ertl, M., Weigand, B.:
Analysis methods for direct numerical simulations of primary breakup of shear-thinning liquid jets. Atomization and Sprays 27(4), 303–317, 2017.
 Reitzle, M., Kieffer-Roth, C., Garcke, H., Weigand, B.:
A volume-of-fluid method for three-dimensional hexagonal solidification processes, J. Comput. Phys. 339: 356-369, 2017.
 Ruberto, S., Reutzsch, J., Roth, N., Weigand, B.:
A systematic experimental study on the evaporation rate of supercooled water droplets at subzero temperatures and varying relative humidity, Exp Fluids, 58:55, 2017.