Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D

ARTICLE IN PRESS

JID: AMC

[m3Gsc;June 8, 2015;11:7]

Applied Mathematics and Computation 000 (2015) 1–10

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D Kathrin Eisenschmidt, Moritz Ertl, Hassan Gomaa, Corine Kieffer-Roth∗, Christian Meister, Philipp Rauschenberger, Martin Reitzle, Karin Schlottke, Bernhard Weigand Institute of Aerospace Thermodynamics, Pfaffenwaldring 31, Stuttgart 70569, Germany

a r t i c l e

i n f o

Article history: Available online xxx Keywords: DNS Droplet Spray Phase change Multi-component Non-Newtonian

a b s t r a c t A brief overview of the main features of the code Free Surface 3D (FS3D) is presented here. FS3D is a Direct Numerical Simulation (DNS) code based on a Volume of Fluid (VOF) method to solve the incompressible Navier–Stokes equations. The energy equation is included and multiphase flows can be considered, in the sense that evaporation and ice growth are already implemented. Furthermore, the code has been expanded to consider multi-component systems and non-Newtonian fluids. Complex phenomena with strong computational effort can be examined, since the code works on massive parallel architectures. © 2015 Published by Elsevier Inc.

1. Introduction The numerical simulation of multiphase flows is important to describe many processes in nature and engineering. Applications are as diverse as weather forecast, administering of medication or spreading of substances in agriculture. Relevant technical applications are, amongst others, fuel injection in form of sprays in combustion chambers as well as water injection in compressors of gas turbines. The program FS3D is able to solve the incompressible Navier–Stokes equations as well as the energy equation with temperature dependent thermo-physical properties by means of Direct Numerical Simulation. It is continuously optimized and expanded with new features and has been in use for more than fifteen years. FS3D includes procedures to account for sophisticated processes, like the simulation of phase transition or turbulent inflow conditions. Applying DNS, and thus resolving even the smallest temporal and spatial scales, FS3D does not need to apply any turbulence modeling. To face the resulting huge computational load, FS3D is parallelized using OpenMP as well as MPI. Many phenomena have been investigated until now. To begin with, basic drop and bubble dynamic processes can be considered. For instance, the impact of a droplet onto a thin film (“splash”) can be simulated for a better understanding of the instabilities at the border of the resulting liquid crown. In air conditioning systems, the prerequisites for the occurrence of secondary droplets (“ splashing”) can be determined in case of the impact of an incident droplet on a droplet adhering to a surface. The oscillation behavior of an evaporating droplet with and without oncoming flow has been investigated [1]. Moreover, even much more complex phenomena like the primary jet breakup can be simulated [2].

∗

Corresponding author. Tel.: +49 71168562413; fax: +49 71168562317. E-mail address: [email protected], [email protected] (C. Kieffer-Roth).

http://dx.doi.org/10.1016/j.amc.2015.05.095 0096-3003/© 2015 Published by Elsevier Inc.

Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC 2

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

The present survey begins with the description of the basic fluid mechanics in the most fundamental case of two adjoining, immiscible, incompressible and Newtonian fluids. Subsequently, the energy equation is included to allow the liquid/solid and the liquid/gas phase transition. Jump conditions at phase boundaries have to be taken into account and the important numerical schemes are explained. Then, the extension to multi-component systems is presented. Finally, non-Newtonian fluids with a shear rate dependent viscosity are considered. 2. Fundamental method 2.1. Numerical approach The code Free Surface 3D (FS3D) allows for Direct Numerical Simulations of multiphase flows including sharp interfaces between immiscible phases. The flow field is computed by solving the conservation equations in a one-field formulation, where the different phases and fluids are regarded as a single fluid with variable thermo-physical properties that discontinuously change across the separating interface. The conservation equations for mass and momentum read accordingly

ρt + ∇ · [ρ u] = 0,

(1)

[ρ u]t + ∇ · [ρ uu] =

∇ · [S − Ip] + ρ g + fγ ,

(2)

with u denoting the velocity vector, ρ the density, p the static pressure and g the gravitational acceleration. For Newtonian fluids, the shear stress tensor S is given by

S=μ



 ∇ u + ( ∇ u )T .

(3)

The treatment of non-Newtonian fluids is discussed separately in Section 5. fγ represents the influence of surface tension in the vicinity of the interface, and is handled as a volume force that acts on every cell containing an interface [3]. The finite volume method is used to discretize the Navier–Stokes equations for incompressible flow in space. According to the Marker and Cell (MAC) method [4], velocities are stored on the cell faces, while scalar variables are stored at the cell centers. The continuity and momentum equations are solved simultaneously to provide the velocities, with a set of one-dimensional transport equations for each spatial direction. A second-order upwind scheme is applied to discretize the convective terms, while a second-order accurate central difference discretization scheme is used for the diffusive terms. The fluxes are computed by Godunov type schemes and a Total Variation Diminishing (TVD) limiter [5] is used to restrict the spatial derivatives to avoid the generation of new velocity maxima. To identify different phases, additional indicator variables fi are used according to the classical Volume of Fluid (VOF) method [6], where fi is defined as

fi (x, t ) =

⎧ ⎨1

in the disperse phase,

]0, 1[ for interfacial cells,

⎩

0

(4)

in the continuous phase.

Here, i is equal to 1 in the liquid, i = 2 in the vapor and i = 3 in the solid phase. For general expressions, the notation f is used in the following. Consistent with the one-field formulation, the local variables, e.g., the density ρ and the dynamic viscosity μ, are defined by the local f value and the respective values for, e.g., the pure gaseous (index g) and liquid (index l) phases

ρ (x, t ) = ρl f (x, t ) + ρg [1 − f (x, t )],

(5)

μ(x, t ) = μl f (x, t ) + μg [1 − f (x, t )].

(6)

The volume fraction f is advected by solving the transport equation:

ft + ∇ · [ f u] = DV.

(7)

The dummy variable DV is equal to zero if no phase change takes place. It represents a source term in case of solidification, and a source term and a diffusion term in case of evaporation, since the vapor phase can diffuse into the gas phase. For the computation of the VOF fluxes, FS3D makes use of the piecewise linear interface reconstruction (PLIC) method [7,8], which reconstructs a plane separating the fluid phase and the gaseous phase in interface cells. This plane is orthogonal to the ˆ γ , which is defined by the negative gradient of the volume fraction f. Therefore, the fluxes can be calculated local normal vector n with first-order accuracy. Second-order accuracy for the interface reconstruction step is achieved by prolongating/extending the interface obtained with the first-order procedure and minimizing the distance to the first-order interfaces in the neighbor cells, thus obtaining a new normal. In the class of second-order advection methods for approximating the volume fluxes, FS3D makes use of two methods described in Section 2.2. The surface tension is computed in terms of either the conservative continuous surface stress (CSS) model suggested by Lafaurie et al. [9], or the continuum surface force model (CSF) by Brackbill et al. [3], or the recently introduced balanced force approach by Popinet [10], allowing a significant reduction of parasitic currents as presented in Section 2.2. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

3

As a result of volume conservation in incompressible flow, a Poisson equation needs to be solved for the pressure

∇·

 1  ∇ ·u ∇p = . ρ( f ) t

(8)

Here, a multigrid solver is applied to the resulting set of linear equations. A Red-Black Gauss–Seidel algorithm is used as a smoother and the whole algorithm can be run as a V- or W-cycle scheme [8,11]. Two time integration schemes are implemented: a first-order explicit Euler scheme and a second-order Runge–Kutta scheme. The maximum time step sizes are restricted by the Courant–Friedrichs–Lewy (CFL) condition [12]. Usually, time step limitations arising from the viscous and surface tension terms are overcome by corresponding implicit treatments. To be able to perform simulations with very high spatial resolutions, it is indispensable to make use of supercomputers like the Cray-XC40 at the High Performance Computing Center Stuttgart (HLRS). Therefore, FS3D is fully parallelized using MPI as well as OpenMP, which allows to use a spatial resolution of more than two billion computational cells per simulation. In Table 3, T the speedup S = T1 for the droplet collision test case with a resolution of 768 × 128 × 128 grid cells is shown. Here, TN (resp. T1 ) N is the wall-clock time of a computation on N (resp. 1) cores. The pressure correction causes up to 90% of the total computational effort, most of it contributed by the smoothing routines. Hence, especially the Red-Black Gauss–Seidel based multigrid solver is highly optimized to work in parallel environments. Therefore, a routine agglomerates the whole domain at a certain level of coarsening and gathers all the restricted solutions on the root process. The process then computes all the coarser levels while all other processes are idling, totally avoiding any MPI communication due to the exchange of boundary conditions after every smoothing cycle. This, in turn, leads to a much better overall performance by reducing the total MPI communications by a factor of four. 2.2. Improvements 2.2.1. Mass advection: numerical methods and test cases In the original advection method, three one-dimensional non-conservative transport equations are considered successively [8]. For that reason, interfaces have to be reconstructed three times. Due to permutation of the equation sequence [13] and a particular choice for the divergence correction, the method is of second-order accuracy in space and time. Fluid is transported parallel to the velocities, which are perpendicular to the cell faces. From a geometric point of view, the volume fluxes are cuboids with strong overlaps. In the alternative formulation, fluid is advected in one step with a more realistic approximation for the volume fluxes [14]. A general, six-faced polyhedron is used for the volume fluxes. Its lateral faces are determined by velocities which are interpolated onto the middle of the edges. Consequently, the fluxes for adjacent faces do not overlap and repeated advection of fluid is reduced. The position of the backface is adjusted in such a way that the volume of the polyhedron is equal to the volumetric flux resulting from the normal velocity components stored at the center of the cell faces. This ensures volume conservation. Numerical errors for both methods have been calculated for deformation-free test cases (pure translation and rotation) as well as for two deformation test cases from literature: the vortex-in-a-box¥ case of Rider and Kothe [7] combined with a laminar pipe flow as in the work of Liovic et al. [15] and the deformation flow field test case of LeVeque [16]. In fact, the new method is more accurate, while the increase in computational effort and memory demand is reasonable. Robustness was achieved by accounting for special cases. 2.2.2. Curvature calculation based on height functions: numerical approach and application Instead of considering the interface as a sharp discontinuity, the interface normal is calculated in the CSF model from the gradient of the volume fraction. This leads to the following expression for the surface tension, considered to be of continuous nature:

˜fγ = σ κ∇ f.

(9)

To guarantee balanced force discretization, according to Francois et al. [17], surface tension has to be calculated at the centers of the cell faces. The evaluation of the gradient of the volume fraction is based on the direct neighbors of the considered cell face. The use of a larger stencil for this evaluation would be one of the causes for parasitic currents, and is not consistent with the discretization of the Poisson equation, where a local coupling of the pressure values is carried out [18]. Moreover, an accurate curvature calculation is needed. It is based on a height function approach coupled with a local paraboloid fitting. In this procedure, inspired by the article of Popinet [10], an adaptive stencil, orientated preferably in the direction of the biggest component of the interface normal, has to be found in order to calculate the mean curvature. In the CSF formulation, the advantages of the height function approach compared to a curvature calculation for an implicit ∇f ) can be shown in Fig. 1. In case of a static droplet in equilibrium, spurious currents could be reduced by function κ = ∇ · ( |∇ f|

two orders of magnitude [18].

Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

(a) Previous CSF model |vmax |

= 33.09 cm s

(b) New CSF model |vmax | = 0.026 cms

Fig. 1. Velocity field around a droplet with radius R = 10−3 m at t = 0.03 s , in an xy plane through the center of the droplet (Cubic computational domain: edge length of 4 × 10−3 m , 643 grid cells in each direction) [18].

3. Phase change models 3.1. Energy equation A couple of simplifications and assumptions are made which on the one hand are necessary to be able to solve the problem and on the other hand are due to negligible physical effects: 1. 2. 3. 4. 5. 6.

Thermal radiation is neglected. Heat production due to viscous dissipation is neglected. The thermodynamic properties are far from the critical point. There are no chemical reactions. Local thermodynamic equilibrium at the interface is assumed. The densities of both liquid water and ice are considered to be equal.

Keeping these simplifications in mind and assuming small Eckert numbers Ec, the energy equation in its temperature form [19] can be written as

∂ (ρ c p T ) + ∇ · (ρ c p T u ) = ∇ · (k∇ T ) + q˙  , ∂t

(10)

q˙ 

is a source term, which will be explained in greater detail in Sections 3.2.2 and 3.3.3. where k is the heat conductivity and In order to solve the energy equation (10), a multiscalar approach is used. This means that separate temperature fields are introduced for each phase. The temperature fields are only defined within the respective phase and are coupled at the boundary interfaces via jump conditions [20]. One advantage of this approach is that fluxes over cell boundaries can be calculated separately and without averages. 3.2. Jump conditions at the phase boundary 3.2.1. Freezing/melting For freezing processes, phase change takes place at the interface with temperature Tγ . This temperature can be obtained from the Gibbs–Thomson equation and depends on the local orientation and curvature κ of the interface as well as the specific heat capacities of the solid (cs ) and the liquid (cp, l ) phase [21]

Tγ = Tm







 c p,l − cs Tm 1

1− σ κ + σζ ζ κ1 + σθ θ κ2 + ρL L

Tγ ln

Tγ Tm





+ Tm − Tγ





.

(11)

Here, Tm is the equilibrium melting temperature for a planar surface, L is the fusion enthalpy and σ is the surface energy. ζ ˆ re f and the vector normal to the interface n ˆ γ in the plane of the and θ describe the relative angle between a reference normal n principal curvatures κ 1 and κ 2 , respectively, as can be seen in Fig. 2. The jump condition for the energy conservation, which is also known as the Stefan condition, then reads as follows:





m˙  L − c p,l − cs



Tm − Tγ



ˆ γ + ks ∇ Ts · n ˆγ , = −kl ∇ Tl · n

(12)

ˆ γ is the unit vector normal to the interface. This equation expresses the fact that where m˙  is the area specific mass source and n the released fusion enthalpy equals the sum of the heat flow from the interface into the ice and into the water. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

[m3Gsc;June 8, 2015;11:7] 5

z

nˆ ref θ nˆ γ

ζ

y

x Fig. 2. Relationship between the interface and the reference normal.

3.2.2. Evaporation/condensation When looking at evaporation processes, the jump condition for the energy conservation is not used. Instead, the occurring evaporation enthalpy is taken care of at the phase interface through the source term q˙  = hv m˙  , where hv is the latent heat of evaporation and m˙  the volumetric mass source of vapor. Due to the phase change, the velocity field can no longer be considered as divergence-free. This is a result of the velocity jump at the phase interface, which is caused by the formation (evaporation) or disappearance (condensation) of volume. A common approach to account for the divergence of the velocity field is

∇ · u = −m˙ 

1 1 − , ρv ρl

(13)

a volume balance, as used by various authors such as [22]. However, as FS3D uses the VOF-method, u is a mass averaged velocity field and therefore cannot be considered as divergence-free for interface cells. It is crucial to understand that there is a difference between a volume weighted velocity and a mass weighted velocity, which additionally contains the densities of the phases (for more details see [23]). These two velocities are identical only in case of fully wetted or completely empty cell faces and only then volume is conserved. However, the common case consists of partly wetted cell faces. Therefore, the volume weighted and the mass weighted velocities differ due to the large density ratio (e.g. water–air: ρ f /ρ gp ≈ 1000). As a consequence, the volume balance of Eq (13) cannot be used to calculate the divergence of u. The problem is solved by using the velocities of the interface and the gaseous phase (consisting of vapor and inert gas), u

and ugp , respectively, to calculate a virtual mass averaged velocity field u∗ . This velocity field in turn is used as the source term in the actual continuity equation

∇ · u = ∇ · u∗ .

(14)

For a description of how to obtain the velocities u and ugp , refer to Section 3.3.4. 3.3. Numerical approach Within FS3D, the convective transport, the diffusion and the phase change are treated in separate steps. In the following sections, the convective transport is neglected and the focus lies on the phase change coupled with its influence on the temperature fields. A prerequisite for the simulation of solidification problems is the handling of three different phases including a solid one. 3.3.1. Rigid body motion in the Eulerian framework of FS3D In FS3D, rigid bodies are advected with the same transport equation as a fluid in the Eulerian framework. This allows to add volume changes during a phase change on the corresponding VOF variable. The implemented method is based mainly on the work of Patankar [24]. The whole computational domain (including the particle domain) is treated as a fluid in the first step. The rigid body velocity, which is calculated from the known divergence-free velocity field as ur = U + ω × r, is then imposed within the zone occupied by the solid. When the VOF variable representing the solid is now advected with velocities fulfilling the rigidity constraint, the particle is truly rigid. This is in contrast to other approaches, where the particle is treated as a Lagrangian one (e.g. [25]). The technique implemented in FS3D is described and validated in [26]. Several different rigid bodies can be identified and a collision model allows their interaction. 3.3.2. Multiphase VOF The treatment of fluid flow problems with more than two phases is rarely researched. The implemented approach is based on the VOF method with Piecewise Linear Interface Calculation (PLIC), which is less costly than e.g. the Moment of Fluid (MOF) method [27], and is suitable for the purpose. The volume fraction of ice, which is the third immiscible phase, is represented by the VOF variable f3 according to Eq. (4). Its advection is described by the transport Eq. (7). The reconstruction via PLIC is difficult if three phases exist in one cell. The interfaces must not overlap and the correct volume fractions have to be included. For the moment, parallel interfaces are forced [28]. A contact angle model will allow any arrangement in the future. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

3.3.3. Freezing/melting The Stefan condition  (12) can be used to derive an expression for the growth velocity Vγ by simply substituting the mass flux per unit area m˙  = ρ Vγ . The heat fluxes in the solid q˙ s and the liquid q˙ l are determined through the Fourier heat conduction law q˙  = −k∇ T, where the gradient of the temperature fields is calculated through second-order one-sided finite differences [29]. Because of the iterative algorithm, the heat fluxes are closely linked to the mass source term m˙  and the released enthalpy of fusion at the interface. Therefore, they can directly be used as the source term of the energy Eq. (10). The mass source m˙  is then calculated using the ˆ γ . The difference of the new position and the PLIC surfaces: The PLIC surface is translated with Vγ t along the normal vector n initial position of the interface can then be used to determine the mass source. If the cell is filled up in the current time step tn → tn+1 , a polyhedron of the initial PLIC surface and its translated counterpart is built and cut with the cell sides. The parts of the polyhedron which protrude into an adjacent cell are added to that cell if it has been empty beforehand. This allows for a natural way of freezing over cell boundaries. Anisotropic growth can be achieved using the full Gibbs–Thomson Eq. (11), if an orientation-dependent surface energy is used [30]. It needs to be mentioned that as of now, three-dimensional particles are only allowed to grow spherically. Anisotropy is currently limited to two-dimensional problems. Because the timestep for an explicit method is proportional to the smallest grid size h cubed [31], an implicit algorithm is used for the treatment of the diffusive heat transport. For an implicit treatment, an additional CFL-condition is given by

αγ =

Vγ t ≤ 1. h

(15)

As a criterion for convergence of this implicit algorithm, an averaged residual of the growth velocities is calculated, which has to be smaller than a certain threshold. An under-relaxation of the heat fluxes and the growth velocity can be used to stabilize the algorithm and improve its convergence behavior. Because the phase change does not only change the mass in every cell but also the cell energy, the temperature in all boundary cells needs to be updated after every inner iteration step [29]. 3.3.4. Evaporation/condensation The volumetric mass source of vapor m˙  is obtained via the area specific vapor mass source m˙  , which is calculated based on the equation

m˙  =

Dvp ρgp ∇ Xd · nˆ γ , 1 − Xd

(16)

following the derivation in [32]. Here, Dvp is the binary diffusion coefficient of vapor inside the gas, ρ gp is the density of the gaseous phase and Xd represents the vapor mass fraction. A detailed explanation of how to obtain the local density of the gaseous phase, the vapor mass fraction and its gradient can be found in [32] and [23]. Once Eq. (16) is solved, the volumetric mass source of vapor can be calculated using the local interface density a

a =

SPLIC , V

(17)

with SPLIC being the exact area of the reconstructed interface and V the cell volume. It is important to note that a needs to be determined based on the PLIC reconstruction instead of using the gradient of the f1 -field. To avoid undershoots of the f1 -variable, the mass source is limited to m˙  = f1 ρ1 /t in cells containing almost no liquid. The increase in volume due to evaporation necessitates for the consideration of two different transport velocities: the gaseous phase gets advected with ugp , the liquid phase with u . These two velocities only differ from each other in cells containing a part of the interface and hence are considered in cells containing the interface only. The iterative algorithm used to compute the two velocities is the following. First, preliminary velocities ugp and u are calculated using a volume flux balance based on Eq. (13) across the non-wetted portion of the cell faces and the mass weighted velocity uρ of the last time step. Second, a preliminary volume weighted velocity uvol is obtained through ugp and u . Its divergence is then compared to the desired divergence produced by the evaporation (see Eq. (13)) and an error E(∇ ·u ) is determined for each cell containing the interface. This divergence error is then distributed vol

to the spatial directions according to the local unit normal vector of the surface, which allows the computation of the final volume weighted velocity uvol . Now ugp and u can be obtained using again the same volume flux balance as before, but now in combination with uvol instead of uρ . Finally, a virtual mass weighted velocity field u∗ is computed by

u∗ = u fA ρl + ugp (1 − fA )ρgp .

(18)

Its divergence is employed as the source term for the continuity equation as given in Eq. (14). It should be mentioned, that a non-iterative method was proposed later [33]. 3.4. Application One example for a phase transition is the growth of a small ice crystal within a supercooled droplet. For this example, the previously mentioned anisotropy of the surface energy is taken into account, thus resulting in dendritic growth. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

128x128

196x196

7

256x256

Fig. 3. Dependence of the orientation of the grid on the 2D crystal growth. Blue crystal is turned −60◦ with respect to the black reference case [21,34]. Table 1 Material properties of liquid water and ice at 263.1 K.

Water Ice

T [K]

[kg/m3 ]

ρ

L [J/kg]

σ

λ

[N/m]

[W/(m K)]

cp [J/(kg K)]

c [J/(kg K)]

263.15

998.2

333600

0.028

0.5386 2.3190

4284 –

– 2030

The material properties for the setup are listed in Table 1. The evolution of an initially round, two-dimensional ice crystal for different resolutions was examined in [21,34] and is shown in Fig. 3. It can be deduced that a resolution of approximately six to eight cells of the tip of the dendrite is sufficient to suppress any grid dependency. It is, however, challenging to determine the size of the smallest structures a priori when the resolution of the setup is set. The method described above will therefore not be able to accurately capture structures which are smaller than the resolution defined before. 4. Multi-component treatment 4.1. Numerical approach With FS3D, it is also possible to study problems involving multiple components in the liquid and vapor phase. To account for the different components within the phases, and to allow a seamless integration into the already implemented differentiation between the liquid, vapor and gaseous phases, a two-scalar approach for each species is chosen. This means, that the volume fraction of the mixed liquid, composed of n different species, is still represented by the VOF variable f, and the volume fraction of the mixed vapor phase is still represented by the VOF variable f2 . In addition, each species i is identified by its own volume fractions  i and Yi , depicting the volume amount of species i in the liquid and gaseous phase of each control volume, respectively. The implemented method itself is thereby subject of the following assumptions: • The whole system is still considered to be incompressible. • Due to a lack of an in-phase surface reconstruction scheme, which could represent liquid-liquid interfaces, the species are represented as an ideal mixture (linear mixing behavior) within each control volume. While FS3D is currently only used to solve advection dominated multi-component problems, the implemented method is especially advantageous for subsequent studies of problems influenced by diffusion processes or even heat and mass transfer between the liquid and the vapor phases. In this case, for each species i, the corresponding volume fractions  i and Yi are transported by:

   m˙   ∇ ρl Di j ∇ i − i , [ f i ]t + ∇ [ f i u ] = ρ ρl,i l j=1 n 

(19)

j =i

[ f2 ϒi ]t + ∇ [ f2 ϒi ugp ] =

   m˙  ϒi ∇ ρgp Di j ∇ + i . ρgp ρv,i j=1

n 

(20)

j =i

Here, u and ugp are the transport velocities of the interface and gaseous phase within each interfacial cells. The mixed densities are represented by ρ l and ρ gp for the liquid and gaseous phases, respectively. Furthermore, Dij denote the different diffusion coefficients and m˙  are the mass sources of the different species due to phase change processes. Not shown here, but i also important in this context, is the inclusion of the mass source terms into the volume balances of f and f2 as well as into the energy equation. Currently, the local fluid properties in each control volume are defined by the volume and species fractions f, f2 ,  i and Yi , following the one-field-formulation and a linear mixing behavior. As an example, the density field is defined by the densities of Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

(a) t 1

(b) t 2

(c) t 3

Fig. 4. Multi-component droplet/wall film interaction.

the mixed liquid (index l) and vapor (index v) phase

ρl =

n 

ρl,i i ρv =

i=1

n 

ρv,i ϒi

(21)

i=1

resulting in the local densities, defined by

ρ = ρl f + ρv f2 + ρgp (1 − f − f2 ).

(22)

4.2. Application A very interesting application of this multi-component treatment is the impact of a droplet consisting of species A onto a wall film of species B, which is depicted in Fig. 4. In this application case one would, amongst others, especially be interested in: • The species composition of the resulting secondary droplets. • The mass amount of species A remaining in the wall film. • The spreading of species A within the wall film. All three can differ significantly, depending on the properties of the different species. 5. Non-Newtonian fluids In many technical and pharmaceutical applications as well as in nature we can observe that many fluids have to be described by using non-Newtonian properties. FS3D has been employed for non-Newtonian droplet oscillations and collisions in the past [35,36]. A more complex example of an application of non-Newtonian liquid are disintegrating jets used in spray drying of polymeric solutions to generate solid particles with defined properties. For fluids with non-Newtonian properties, the viscosity μ is not constant. FS3D is able to simulate fluids that have a viscosity which changes with the shear rate γ˙ . To calculate the shear dependent viscosity we use the Carreau–Yasuda model [37]: )  (n−1 μ(γ˙ ) − μ∞  a a = 1 + (τ γ˙ ) . μ0 − μ∞

(23)

The model uses five constants: the viscosity at zero shear μ0 , the viscosity at very high shear rates μ∞ and three model parameters τ , a and n. The values for these parameters are obtained from experimental data. Two other models for shear dependent viscosities, the Zhou model and the Cross model, are also available in FS3D [38]. The viscosity is used to calculate the shear stress tensor S of the momentum conservation equation (2). The shear rate is calculated from the velocity field at the last time step, as follows:



 0.5 γ˙ = 0.5 ∇ u + (∇ u )T 2 .

(24)

Currently, the viscosity models in FS3D are only able to model shear dependent viscosities. Therefore, the simulated flows are not allowed to have visco–elastic properties. There are two dimensionless groups which are used to check for visco–elastic behavior. The elasto-capillary number and the Deborah number, defined by:

El =

λσ λ , De =  , η0 l ρ l3

(25)

σ

where λ is the characteristic relaxation time of the fluid, η0 is the dynamic viscosity and l is a characteristic length. The boundaries for a solely shear dependent fluid are De < <0.35 and El < <2.35. The injection of a liquid jet of an aqueous polyvinylpyrrolidone (PVP K30) solution with 20 wt.% from a dual nozzle velocity profile into an air atmosphere at 30 bar is simulated as an example. The material properties of the jet and the atmosphere are given in Table 2. The parameters for the Carreau–Yasuda model are μ∞ = 0.0254 mPa s, a = 22.045, n = −0.4886 and τ = 7.14. The Deborah number for PVP K30 20 wt.% is De = 8.17 × 10−8 and the elasto-capilary number is El = 1.27 × 10−8 , putting the used solution in the region of a shear dependent fluid. The results of the simulation are shown in Fig. 5. The high shear introduced into the jet by the inlet velocity profile is clearly visible in the upper section. The corresponding reduction in viscosity for that region can be observed in the lower section. The lower viscosity allows for an easier breakup of the jet and an improved atomization. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

JID: AMC

ARTICLE IN PRESS K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

[m3Gsc;June 8, 2015;11:7] 9

Fig. 5. PVP K30 20 wt.% jet injected into air at ambient pressure p = 30 bar. 2D slice of the inflow plane displaying the shear rate and the velocity vectors in the upper section and the viscosity μ/μ0 in the lower section. Table 2 Material properties for the simulation.

PVP K30 20 wt.% Air at 30 bar

Density ρ [kg/m3 ]

Viscosity μ0 [Pa s]

Surface tension σ [N/m]

1049 38.1

5.72 × 10−5 1.79 × 10−8

64.373 × 10−3

Table 3 Speedup with MPI (droplet collision, 768 × 128 × 128 cells). Number of cores Speedup S

1 1

8 3

64 27

512 96

6. Conclusion In this paper, it has been seen that with the multiphase code FS3D for incompressible flows, a wide range of physical phenomena can be simulated. Among these, we investigate the freezing and evaporation phase transitions. Multi-component systems and fluids with Newtonian as well as non-Newtonian properties can be considered. Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095

JID: AMC 10

ARTICLE IN PRESS

[m3Gsc;June 8, 2015;11:7]

K. Eisenschmidt et al. / Applied Mathematics and Computation 000 (2015) 1–10

The major findings have been • • • •

Phase changes can be reproduced precisely. The behavior of complex rigid body systems is captured. Multicomponent systems are described properly. An appropriate modeling of the viscosity for non-Newtonian fluids is disposable.

Acknowledgments Part of the simulations were performed on the national supercomputers at the High Performance Computing Center Stuttgart (HLRS). The authors would like to thank the German Research Foundation (DFG) for financial support of the project within the SFB-TRR 75 and within the Cluster of Excellence in Simulation Technology at the University of Stuttgart. References [1] J. Schlottke, E. Dulger, B. Weigand, A VOF-based 3D numerical investigation of evaporating, deformed droplets, Prog. Comput. Fluid Dyn. 9 (6–7) (2009) 426–435. [2] W. Sander, B. Weigand, Direct numerical simulation and analysis of instability enhancing parameters in liquid sheets at moderate Reynolds numbers, Phys. Fluids 20 (5) (2008) 053301. [3] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface-tension, J. Comput. Phys. 100 (2) (1992) 335–354. [4] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182–2189. [5] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys. 32 (1) (1979) 101–136. [6] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201–225, doi:10.1016/00219991(81)90145-5. [7] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (2) (1998) 112–152, doi:10.1006/jcph.1998.5906. [8] M. Rieber, Numerische Modellierung der Dynamik freier Grenzflächen in Zweiphasenströmungen, Dissertation, Universitt Stuttgart, Germany (2004). [9] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with surfer, J. Comput. Phys. 113 (1) (1994) 134–147. [10] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (16) (2009) 5838–5866. [11] P. Wesseling, An Introduction to Multigrid Methods, John Wiley & Sons, 1992. [12] C. Hirsch, Numerical Computation of Internal & External Flows, John Wiley & Sons, Ltd., 2007. [13] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (3) (1968) 506–517. [14] J. Hernández, J. López, P. Gómez, C. Zanzi, F. Faura, A new volume of fluid method in three dimensions—part I: multidimensional advection method with face-matched flux polyhedra, Int. J. Numer. Methods Fluids 58 (8) (2008) 897–921. [15] P. Liovic, M. Rudman, J.-L. Liow, D. Lakehal, D. Kothe, A 3D unsplit-advection volume tracking algorithm with planarity-preserving interface reconstruction, Comput. Fluids 35 (2006) 1011–1032. [16] R.J. LeVeque, High-resolution conservative algorithms for advection in incompressible flow, SIAM 33 (1996) 627–665. [17] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (1) (2006) 141–173. [18] M. Boger, J. Schlottke, C.-D. Munz, B. Weigand, Reduction of parasitic currents in the DNS VOF code FS3D, in: 12th Workshop on Two-Phase Flow Predictions, March 22–25, 2010, Halle (Saale), Germany, 2010. [19] H.D. Baehr, K. Stephan, Heat and Mass Transfer, 2nd ed., Springer Verlag, 2006. [20] J. Schlottke, P. Rauschenberger, B. Weigand, C. Ma, D. Bothe, Volume of fluid direct numerical simulation of heat and mass transfer using sharp temperature and concentration fields, in: ILASS – Europe 2011, 24th European Conference on Liquid Atomization and Spray Systems, Estoril, Portugal, 2011 http://www.ilass.uci.edu/ . [21] P. Rauschenberger, Ein Verfahren für die direkte numerische Simulation des Gefriervorgangs von unterkühltem Wasser, Dissertation, Universität Stuttgart (2014). [22] G. Son, V.K. Dhir, Numerical simulation of film boiling near critical pressures with a level set method, J. Heat Transfer-Trans. ASME 120 (1) (1998) 183–192. [23] J. Schlottke, B. Weigand, Direct numerical simulation of evaporating droplets, J. Comput. Phys. 227 (10) (2008) 5215–5237. [24] N.A. Patankar, A formulation for fast computations of rigid particulate flows, Center for Turbulence Research, Ann. Res. Briefs (2001) 185–196. [25] N. Sharma, N.A. Patankar, A fast computation technique for the direct numerical simulation of rigid particulate flows, J. Comput. Phys. 205 (2005) 439–457. [26] P. Rauschenberger, J. Schlottke, K. Eisenschmidt, B. Weigand, Direct numerical simulation of multiphase flow with rigid body motion in an Eulerian framework, in: ILASS – Europe 2011, 24th European Conference on Liquid Atomization and Spray Systems, Estoril, Portugal, 2011 http://www.ilass.uci.edu/ . [27] V. Dyadechko, M. Shashkov, Moment-of-fluid interface reconstruction, Los Alamos National Laboratory, Los Alamos, USA, 2005, pp. 1–39. [28] K. Eisenschmidt, P. Rauschenberger, B. Weigand, Freezing processes of supercooled droplets in clouds: a multiple vof variable approach, in: ILASS – Europe 2011, 24th European Conference on Liquid Atomization and Spray Systems, Estoril, Portugal, 2011 http://www.ilass.uci.edu/ . ´ Ž. Tukovic, ´ I.V. Roisman, B. Weigand, C. Tropea, Comparative assessment [29] P. Rauschenberger, A. Criscione, K. Eisenschmidt, D. Kintea, S. Jakirlic, of volume-of-fluid and level-set methods by relevance to dendritic ice growth in supercooled water, J. Comput. Fluids 79 (0) (2013) 44–52, doi:10.1016/j.compfluid.2013.03.010. [30] M. Muschol, D. Liu, H.Z. Cummins, Surface-tension-anisotropy measurements of succinonitrile and pivalic acid: comparison with microscopic solvability theory, Phys. Rev. A 46 (2) (1992) 1038–1050, doi:10.1103/PhysRevA.46.1038. [31] T.Y. Hou, J.S. Lowengrub, M.J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114 (1994) 312–338. [32] A.F. Mills, Mass Transfer, Prentice Hall, 2001. [33] C. Ma, D. Bothe, Numerical modeling of thermocapillary two-phase flows with evaporation using a two-scalar approach for heat transfer, J. Comput. Phys. 233 (2013) 552–573. [34] B. Weigand, C. Tropea, A. Birkefeld, The Collaborative Research Center SFB-TRR 75: droplet dynamics under extreme boundary conditions, in: ILASS – Europe 2011, 26th European Conference on Liquid Atomization and Spray Systems, Bremen, Germany, 2014 http://www.ilass.uci.edu/ . [35] M. Ertl, N. Roth, G. Brenn, H. Gomaa, B. Weigand, Simulations and experiments on shape oscillations of newtonian and non-newtonian liquid droplets, in: ILASS – Europe 2013, 25th European Conference on Liquid Atomization and Spray Systems, Chania Crete, Greece, 2013, p. 7. [36] M. Motzigemba, N. Roth, D. Bothe, H.-J. Warnecke, J. Prüss, K. Wielage, B. Weigand, The effect of non-newtonian flow behaviour on binary droplet collisions: VOF-simulation and experimental analysis, in: ILASS – Europe 2011, 26th European Conference on Liquid Atomization and Spray Systems, Zaragoza, Spain, 2002 http://www.ilass.uci.edu/ . [37] R.I. Tanner, Engineering Rheology, second ed., Oxford Engineering Science Series, Oxford University Press, 2002. [38] M.M. Cross, Rheology of non-newtonian fluids: a new flow equation for pseudoplastic systems, J. Colloid Sci. 20 (5) (1965) 417–437.

Please cite this article as: K. Eisenschmidt et al., Direct numerical simulations for multiphase flows: An overview of the multiphase code FS3D, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.05.095