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,  and  examine the atomization of non-Newtonian fluids, while    and  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 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  implemented a modified approach for wall modeling in pipes and channels and Habla et al.  created a solver for viscoelastic two-phase flows. Moreover, in  the ISPH algorithm was extended by an enhanced interface treatment procedure, while in  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 .
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 OpenFOAM , 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.
2.1. 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 compressibleInterFoam. 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 . 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 corresponds to a cell completely filled with liquid and 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 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 velocity vector and the numerical implementation can be found in . 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 . Considering the fluid temperature T and the thermophysical behavior of liquids and gases , 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 (  ), and common properties for water ( ,  ) 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 and the gravitational acceleration is defined as stated in Equation (7).
The surface integral with the constant surface tension denotes the force acting at the liquid-gas interface , where determines the unit normal to the interface and twice the mean curvature of the interface. The Dirac delta in three dimensions is expressed by , containing , a point on the surface and , the point where the equation is calculated .
Finally, the energy equation applied in compressibleInterFoam  is stated in Equation (8).
The kinetic energy is expressed by and the thermal diffusivity of the mixture by . The specific isochoric heat capacity of the liquid and the gas phase is represented by and , respectively, where values of 4182 J∙kg−1∙K−1 and 1007 J∙kg−1∙K−1  have been used for each simulation shown in this paper.
2.2. 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 .
One of the first developed non-Newtonian model is the well-known two parameter equation called Power Law, published by Reiner in 1926 . Beside the strain rate, the fluid specific flow consistency index and the flow behavior index are used to calculate the viscosity field.
In the very same year, Herschel and Bulkley propagated their widely used minima function model , considering a minimal viscosity and a threshold strain stress combined with the Power Law equation.
However, the model presented by Casson in 1959  consists only of a threshold strain stress and one additional fluid specific parameter .
In contrast, the Cross Power Law model developed in 1965  contains a minimum and a maximum viscosity and , and the two rheological parameters and .
The Bird-Carreau equation published in 1972 is also expressed by a minimum and a maximum value as well as by two parameters and .
2.3. 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) as written in Listing 2. In order to pass the strain rate to each model specific viscosity calculation function, e.g. CrosspowerlawTransortI.H, the transferred parameters had to be extended as shown in Listing 3. Furthermore, the specific viscosity computation is implemented for each model. To avoid calculation errors, a strain rate of 0 s−1 is replaced by a value of 1E-10 s−1. Related to Equation (13), Listing 4 exemplarily shows the implementation of the Cross Power Law model in OpenFOAM. To ensure a pleasant input option for the characteristic fluid properties, the standard structure of the software is maintained so that the model-specific parameters can be entered in the respective file in the constant folder, e.g. thermophysicalProperties.water. For this purpose, appropriate variables need to be added to the respective header files and initialized in the c files. Listing 5 illustrates the initialization for the four new parameters related to the
Listing 1. OpenFOAM: Accessing the current velocity field.
Listing 2. OpenFOAM: Computation of the strain rate field.
Listing 3. OpenFOAM: Extension of the passed variables by the strain rate components.
Listing 4. OpenFOAM: Computation of the Cross Power Law transport model.
Listing 5. OpenFOAM: Initialization of the new fluid parameters in CrosspowerlawTransport.C.
Cross Power Law. Finally, the storage of all relevant fields in respect to the defined output conditions enables a comprehensive analysis of the simulation results. As an example, the appropriate source code for the strain rate is provided in Listing 6.
2.4. Numerical Test Case
To verify the implementations, pressure driven test simulations for each viscosity model have been carried out. The used computational domain is a three-dimensional horizontal backward facing step as depicted in Figure 1. The mesh consists of 246,800 hexahedral elements with a maximum edge length of 0.85 mm. In combination with a smooth transition to a minimum grid spacing of 0.2 mm at all walls and close to the step area with no hanging nodes, an acceptable numerical accuracy is achieved. Since the focus within this work is on the validation of the novel implementations, a grid study, as well as further optimization of the numerical solution properties has been neglected.
At the initial time t = 0 the whole geometry is filled with the respective liquid, and air is entering the domain with a constant viscosity of 1.8E−5 Pa∙s . The examined non-Newtonian liquids have been chosen from literature in relation to the implemented models. Their rheological properties are listed in 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
Listing 6. OpenFOAM: Storage of the strain rate field in heRhoThermo.C.
Figure 1. Computational domain: Three-dimentional backward facing step.
Table 1. Rheological flow parameters of water and the examined non-Newtonian liquids.
turbulence parameters for the applied k-ω-SST model developed by Menter . 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 (p: 0.2; T, and ω: 0.3; U: 0.7) and 4 alphaSubCycles are defined to improve the stability of the multiphase computation. For the calculation of the time derivatives,
Table 2. Settings of boundary conditions.
a first order accurate implicit Euler scheme was utilized. The convective terms in the volume fraction equation were processed with a bounded limited linear differencing scheme. Also, in view of stability, the convective terms of the turbulence equations were approximated with a diffusive first-order upwind differencing scheme. Apart from the terms in the temperature equation (limited linear differencing scheme), all other convective terms were discretized with an unbounded second-order linear upwind differencing scheme.
The simulations have been computed in parallel using 12 processors (Intel Xeon Platinum 8160 and Intel Xeon E5-2650, respectively). Depending on the non-Newtonian liquid, the CPU time to reach the instances in time shown in Figure 2 varied between 30 h and 35 h, worth mentioning that the simulation with water took just 24 h.
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 2 illustrates the cross section distribution of the volume fraction when the inflowing air has passed the step. Due to the highly varying properties of the examined mixtures, a very different air penetration of the cavity can be observed. Therefore, two different instances in time (0.15 s and 0.23 s) have been chosen to point out the specific vortex formation at the step.
The two higher viscous liquids, utilized to verify the Herschel-Bulkley and the Power Law model, move slower and tend to prevent the formation of vortices right after the step. These liquids also adhere more strongly to the wall, which is particularly noticeable close to the inlet area at the upper boundary. This results from the higher inertia of the liquids, which prevent a faster inflow of the compressed
Figure 2. Phase distribution of the examined mixtures. Transport model of the liquid: a) Const, b) Bird-Carreau, c) Cross Power Law, d) Casson, e) Herschel-Bulkley, f) Power Law.
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
Figure 3. Viscosity distribution of the examined mixtures. Transport model of the liquid: a) Const, b) Bird-Carreau, c) Cross Power Law, d) Casson, e) Herschel-Bulkley, f) Power Law.
Figure 4. Phase distribution for certain instances in time. Transport model of the liquid: Cross Power Law.
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. Finally, 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.
4. 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
Figure 5. Viscosity distribution for certain instances in time. Transport model of the liquid: Cross Power Law.
Figure 6. Chosen sample points for the validation of the Cross Power Law model.
Figure 7. Comparison of theoretical and simulated viscosity quantities. Transport models: Power-Law, Casson, Bird-Carreau.
Figure 8. Comparison of theoretical and simulated viscosity quantities. Transport models: Herschel-Bulkley, Cross Power Law.
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.  or  ), but none of them have yet been implemented in OpenFOAM. This should be part of a future work.
The authors gratefully acknowledge the computing resources granted by RWTH Aachen University on the supercomputer RWTH Compute Cluster, the computing time at the Center for Computational Sciences and Simulation of the University of Duisburg-Essen on the supercomputer magnitUDE, DFG grants INST 20876/209-1 FUGG, INST 20876/243-1 FUGG, the support by the Open Access Publication Fund of the University of Duisburg-Essen and the financial support provided by the Federal Ministry for Economic Affairs and Energy on the basis of a decision by the German Bundestag. The presented results are part of the project MoNNitor in the Central Innovation Programme for SMEs (ZIM).
 Shin, J.H., Lee, I., Kim, H. and Koo, J. (2012) Spray Atomization and Structure of Supersonic Liquid Jet with Various Viscosities of Non-Newtonian Fluids. Open Journal of Fluid Dynamics, 2, 297-304.
 Mansour, M., Kawahara, A. and Sadatomi, M. (2015) Experimental Investigation of Gas-Non-Newtonian Liquid Two-Phase Flows from T-Junction Mixer in Rectangular Microchannel. International Journal of Multiphase Flow, 72, 263-274.
 Santoso, A., Goto, D., Takehira, T., Aslam, A., Kawahara, A. and Sadatomi, M. (2016) Non-Newtonian Two-Phase Flow Characteristics across Sudden Expansion in Horizontal Rectangular Minichannel. World Journal of Mechanics, 6, 257-272.
 Gkotsis, P.K., Evgenidis, S.P. and Karapantsios, T.D. (2019) Influence of Newtonian and Non-Newtonian Fluid Behaviour on Void Fraction and Bubble Size for a Gas-Liquid Flow of Sub-Millimeter Bubbles at Low Void Fractions. Experimental Thermal and Fluid Science, 109, Article ID: 109912.
 Sepulveda, J., Montillet, A., Della Valle, D., Loisel, C. and Riaublanc, A. (2020) Deformation of Gas-Liquid Interfaces in a Non-Newtonian Fluid at High Throughputs inside a Microfluidic Device and Effect of an Expansion on Bubble Breakup Mechanisms. Chemical Engineering Science, 213, Article ID: 115377.
 Habla, F., Marschall, H., Hinrichsen, O., Dietsche, L., Jasak, H. and Favero, J.L. (2011) Numerical Simulation of Viscoelastic Two-Phase Flows Using OpenFOAM. Chemical Engineering Science, 66, 5489-5496.
 Zainali, A., Tofighi, N., Shadloo, M. and Yildiz, M. (2013) Numerical Investigation of Newtonian and Non-Newtonian Multiphase Flows Using ISPH Method. Computer Methods in Applied Mechanics and Engineering, 254, 99-113.
 Amani, A., Balcazar, N., Naseri, A. and Rigola, J. (2020) A Numerical Approach for Non-Newtonian Two-Phase Flows Using a Conservative Level-Set Method. Chemical Engineering Journal, 385, Article ID: 123896.
 Weller, H.G., Tabor, G., Jasak, H. and Fureby, C. (1998) A Tensorial Approach to Computational Continuum Mechanics Using Object-Oriented Techniques. Computers in Physics, 12, 620-631.
 Miller, S.T., Jasak, H., Boger, D.A., Paterson, E.G. and Nedungadi, A. (2013) A Pressure-Based, Compressible, Two-Phase Flow Finite Volume Method for Underwater Explosions. Computers & Fluids, 87, 132-143.
 Suponitsky, V., Froese, A. and Barsky, S. (2014) Richtmyer-Meshkov Instability of a Liquid-Gas Interface Driven by a Cylindrical Imploding Pressure Wave. Computers & Fluids, 89, 1-19.
 Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S. and Jan, Y.-J. (2001) A Front-Tracking Method for the Computations of Multiphase Flow. Journal of Computational Physics, 169, 708-759.
 Koch, M., Lechner, C., Reuter, F., Kohler, K., Mettin, R. and Lauterborn, W. (2016) Numerical Modeling of Laser Generated Cavitation Bubbles with the Finite Volume and Volume of Fluid Method, Using OpenFOAM. Computers & Fluids, 126, 71-90.
 Cho, Y. and Kensey, K.R. (1991) Effects of the Non-Newtonian Viscosity of Blood on Flows in a Diseased Arterial Vessel. Part 1: Steady Flows. Biorheology, 28, 241-262.
 Khalili Garakani, A., Mostoufi, N., Sadeghi, F., Hosseinzadeh, M., Fatourechi, H., Sarrafzadeh, M. and Mehrnia, M.R. (2011) Comparison between Different Models for Rheological Characterization of Activated Sludge. Journal of Environmental Health Science and Engineering, 8, 255-264.
 Koocheki, A., Ghandi, A., Razavi, S.M.A., Mortazavi, S.A. and Vasiljevic, T. (2009) The Rheological Properties of Ketchup as a Function of Different Hydrocolloids and Temperature. International Journal of Food Science and Technology, 44, 596-602.
 Gavrilov, A.A. and Rudyak V. (2015) Reynolds-Averaged Modeling of Turbulent Flows of Power-Law Fluids. Journal of Non-Newtonian Fluid Mechanics, 227, 45-55.
 Mehta, D., Thota Radhakrishnan, A.K., van Lier, J. and Clemens, F. (2018) A Wall Boundary Condition for the Simulation of a Turbulent Non-Newtonian Domestic Slurry in Pipes. Water, 10, 124.