Implementation of Non-Newtonian Fluid Properties for Compressible Multiphase Flows in OpenFOAM

The paper presents the implementation of non-Newtonian fluid properties for compressible multiphase solver in the open source framework OpenFOAM. The transport models for Power Law, Cross Power Law, Casson, Bird-Carreau and Herschel-Bulkley fluids were included in the thermophysical model library. Appropriate non-Newtonian liquids have been chosen from literature, and pressure driven test simulations are carried out. Therefore, the solver compressibleInterFoam is used to compute air-liquid mixture flows over a backward facing step. A validation of the novel models has been performed by means of a sample-based comparison of the strain rate viscosity relation. The theoretical rheological properties of the selected liquids agree well with the results of the simulated data.


Introduction
Non-Newtonian liquids are part of almost all areas of our daily life, such as toothpaste, ketchup, concrete, lubrication oils, polymer melts or blood, to name just a few. In terms of process engineering and transport processes, flows of non-Newtonian fluids frequently occur in combination with compressible media (e.g. polyethylene foam and protein foam). Therefore, this field is of great interest for scientific research. For instance, [1] and [2] examine the atomization of non-Newtonian fluids, while [3] [4] [5] and [6] investigate the flow characteristics of gas bubbles in non-Newtonian liquids under certain conditions. With the increase in computational power, the use of numerical simulations How to cite this paper: Westermaier, S. and Kowalczyk, W. (2020) Implementation of Non-Newtonian Fluid Properties for Compressible Multiphase Flows in Open-FOAM. Open Journal of Fluid Dynamics, to optimize these complex flow patterns becomes more and more attractive. Concerning licensing costs and scalability, a computation with open source software is highly desirable. Thus, different numerical models and techniques for non-Newtonian multiphase flows were developed and implemented in open source software already. In this context, Sawko [7] implemented a modified approach for wall modeling in pipes and channels and Habla et al. [8] created a solver for viscoelastic two-phase flows. Moreover, in [9] the ISPH algorithm was extended by an enhanced interface treatment procedure, while in [10] a novel stress formulation was added. Numerical instabilities concerning the high Weissenberg number problem or stability issues caused by surface tensions were some of the objectives of the method proposed in [11].
However, the majority of the published studies focus on incompressible cases and a lack of open source applications for compressible issues have to be noted. Even within the widespread computational fluid dynamics framework Open-FOAM [12], no non-Newtonian models for compressible multiphase flows are available up to now. Therefore, the purpose of the current study is the implementation and validation of five common non-Newtonian models for compressible multiphase solver in OpenFOAM 5.x.

Compressible Multiphase Solver CompressibleInterFoam
The implementation in OpenFOAM is preceded by an extension of the thermophysical model library, which is used inter alia by the solver compressibleInter-Foam. On this note, a short explanation of the numerical application is given.
CompressibleInterFoam is a solver for two compressible, immiscible, non-isothermal phases (liquids and/or gases), whereby the interface is captured by the volume of fluid approach [13]. Consistently the Navier-Stokes equations (see Equations (6)-(9)) are solved for one fluid, where the considered fluid parameters are related to the phase distribution within the cell. The volume fraction of the liquid is represented by α , accordingly 1 α = corresponds to a cell completely filled with liquid and 0 α = to a cell full of gas. Based on the assumption of a homogeneous mixture, the cell specific density and viscosity ρ and µ are calculated by To track the interface of the two compressible phases, a transport equation for the volume fraction related to the fluid velocity vector U is used.
The third term on the left-hand side of Equation (3) was created to sharpen the interface and avoid numerical diffusion. Details regarding the relative veloc-ity vector rel U and the numerical implementation can be found in [14]. The right-hand side takes into account the pressure p and its influence on the densities of both phases in relation to their specific compressibility ψ [15]. Considering the fluid temperature T and the thermophysical behavior of liquids and gases [16], the equations of states can be obtained.
In terms of comparability, all presented studies have been carried out with general properties related to air and water. Hence, the gas constant for air ( [18]) have been used. Based on the calculated mixture viscosity (see Equation (1)), the total mass continuity equation can be written as Furthermore, the single momentum equation, including the dynamic viscosity µ , the unit tensor I and the gravitational acceleration g is defined as stated in Equation (7).
The surface integral with the constant surface tension σ denotes the force acting at the liquid-gas interface [19], where N determines the unit normal to the interface and κ twice the mean curvature of the interface. The Dirac delta in three dimensions is expressed by ( ) δ − x x , containing x , a point on the surface and x , the point where the equation is calculated [20].
Finally, the energy equation applied in compressibleInterFoam [21] is stated in Equation (8).
The kinetic energy is expressed by  [17] have been used for each simulation shown in this paper. Open Journal of Fluid Dynamics

Non-Newtonian Models
Many different strain rate based viscosity models are existing for non-Newtonian fluids. Therefore, the five most common ones have been implemented within this work and are presented in the following section, whereas the strain rate γ is calculated in the same way as stated in the strainRateFunction of OpenFOAM [22].
One of the first developed non-Newtonian model is the well-known two parameter equation called Power Law, published by Reiner in 1926 [23]. Beside the strain rate, the fluid specific flow consistency index 2 k and the flow behavior index n are used to calculate the viscosity field.
In the very same year, Herschel and Bulkley propagated their widely used minima function model [24], considering a minimal viscosity 0 µ and a threshold strain stress 0 τ combined with the Power Law equation.
However, the model presented by Casson in 1959 [25] consists only of a threshold strain stress and one additional fluid specific parameter m .
The Bird-Carreau equation published in 1972 is also expressed by a minimum and a maximum value as well as by two parameters 3 k and n [27].

Implementation in OpenFOAM
The above stated non-Newtonian transport models were implemented in the thermophysical model library of OpenFOAM 5.x. To this end, the major implementation steps are briefly explained. First, the appropriate file and folder structure was created for each model and the header files are included in psiThermos.C and rhoThermos.C, respectively.
Up to now, only temperature and pressure-based viscosity models were available in the thermophysical model library of OpenFOAM. Hence, the current velocity vector has to be accessed within the viscosity calculation loop according to Listing 1. Subsequently, a scalar field for the strain rate can be implemented in both heRhoThermo.C and hePsiThermo.C and derived in relation to Equation (9)     turbulence parameters for the applied k-ω-SST model developed by Menter [28].
A residual controlled Pimple algorithm with 50 outer and 4 inner corrector loops has been used. The first time step was determined as 1E−7 s and an automatic modification in relation to a maximum Courant number of 0.2 has been enabled. The tolerance criteria were given by 1E−6 and the residual control criteria related to the initial value were implemented with 1E−5, both with a relative tolerance of 0. Furthermore, under relaxation factors for the non-final steps

Numerical Test Case
To verify the implementations, pressure driven test simulations for each vis-  Table 1. In the interests of comparability, a simulation with water, using the already implemented const transport model, has also been conducted.
All other properties of the numerical setup remain the same for all simulations. Concerning this, Table 2 shows the boundary conditions including the Open Journal of Fluid Dynamics

Results
As expected, the numerical results of the performed test cases reveal that the usage of liquids with different viscosity characteristics leads to significant changes in the behavior of the compressible two-phase flow. Figure   air phase and the development of larger velocity gradients. The corresponding visualization of the mixture-dependent dynamic viscosity values is depicted in Figure 3. Different scales have been used to ensure better visibility of the distribution. The results show a physically consistent viscosity variation related to the local velocity gradient for the non-Newtonian liquids and a constant viscosity for water. As an example for the temporal development of the investigated inflow processes, Figure 4 and Figure 5 show the phase and viscosity distribution for certain instances in time of the simulation performed with the Cross Power Law model. The general objective of the presented work was the implementation and validation of non-Newtonian models for compressible multiphase solvers. Hence, a quantitative comparison of the theoretical and the simulated viscosity values was carried out for each model, based on the time steps depicted in Figure 2. Initially after the simulation, the calculated velocity field was used to calculate the local strain rates with the post-processing utility of OpenFOAM. In general, the strain rate field computed during the simulation (see Listing 2) could also be applied. In order to detect possible errors in the recent implementation, the verification of the models was performed with results determined by the existing source code of the software. Subsequently, the calculated strain rates and the already computed viscosity field were extracted. To achieve comparability with the theoretical values, only cells completely filled with liquid were taken into account. The extracted cells were sorted according to their viscosity value and twenty Open Journal of Fluid Dynamics  samples were selected in such a way, that a wide range of the present viscosity spectrum is represented. Figure 6 exemplarily shows the chosen cells for the validation of the Cross Power Law transport model, illustrated by yellow dots. Fi-nally, these simulation data sets are plotted against the theoretical quantities of the rheological fluid properties in respect to the applied strain rate. Figure 7 and Figure 8 show a perfect fitting for all non-Newtonian models. Thus, the validation confirms, that the implementation and the workspace management were successful and no numerical influences interfere with the new transport models.

Conclusions and Outlook
The main objective of the presented work was the allocation of non-Newtonian models for compressible multiphase solver in OpenFOAM. Therefore, the thermophysical model library has been extended by new transport models using the equations for Power Law, Cross Power Law, Casson, Bird-Carreau and Herschel-Bulkley fluids. Considering the structure of the software and the ease of use, an appropriate input option for the necessary fluid parameters has been realized.
Pressure driven test simulations have been carried out for each non-Newtonian model using the solver compressibleInterFoam. Themultiphase flows over a backward facing step reveal a physical flow behavior of the gas liquid mixtures. The increase of computation time compared to a Newtonian mixture was within reasonable limits (about 25% to 45%). Furthermore, the sample-based comparisons of the simulated strain rate viscosity relations show an exact match with regard to the theoretical flow properties of the selected liquids. Thus, the presented results    verify the correct implementation of the rheological transport models in the source code of OpenFOAM 5.x. As part of an open source software without licence costs, these models can be used for large-scale parallel computations of compressible multiphase flows with non-Newtonian fluids. However, as a limiting factor of this work has to be mentioned that the usage of common turbulence models and wall functions for the computation of non-Newtonian liquids may lead to accuracy issues. Specific approaches for some applications have already been developed (e.g. [33] or [34]), but none of them have yet been implemented in OpenFOAM. This should be part of a future work.