PHOENICS provides a selection of representative models for the simulation of inelastic, time-independent, viscous non-Newtonian fluids.
Non-Newtonian viscous models belong to one of two categories, namely fluids which exhibit a yield stress, and those fluids which do not.
The following models are provided for fluids without a yield stress:
The following models are provided for fluids with a yield stress, also called viscoplastic fluids:
Such non-Newtonian fluids can either be pseudoplastic or dilitant, which means they are either shear-thinning or shear-thickening fluids. The former are fluids for which the apparent viscosity of the fluid decreases with increasing shear rate. Some examples are blood, milk, tomato ketchup and mayonnaise.
On the other hand, dilitant fluids are those for which the apparent viscosity increases with increasing shear rate. These fluids are less common than pseudoplastics, and are found only in a few applications, such as for example, china-clay suspensions, titanium dioxide, and mixtures of corn flour or sand in water.
As mentioned already, viscoplastic non-newtonian fluids have a yield stress, which must be exceeded to produce motion in the fluid. After this critical value of the yield stress is attained, these fluids behave like normal purely viscous fluids. The constitutive relations are distinguished by the presence of a yield stress, and this can cause difficulties in the numerical modelling of these fluids through their discontinuous character. Therefore, the Papanastasiou (1987) variant of these models has been provided for use in flows involving both yielded and unyielded flow.
The mathematical models of the constitutive equations for purely viscous fluids are founded on the assumption that the stress tensor is only a function of the rate of strain tensor, as described in the next section.
For both incompressible Newtonian fluids and non-Newtonian viscous fluids, the stress-strain relationship may be written in terms of the generalised Newtonian assumption ( see for example Singh et al (2017) ), as follows:
| τij = 2μaDij | (1.2.1) |
where τij is the shear stress tensor, μa is the coefficient of dynamic viscosity; and Dij is the strain rate tensor, defined by:
| Dij =½(∇u+∇uT) | (1.2.2) |
Equation (1.2.1) can be written as
| τij = μaγ | (1.2.3) |
where γ (in units of s-1) is the shear rate or second invariant of the strain rate tensor Dij, which is determined as:
| γ = √E | (1.2.4) |
where
| E = 2 DijDij | (1.2.5) |
is the generation function in PHOENICS, which is defined by the 3d stored variable GEN1. In the CFD solver, GEN1 is stored in the F-array location identified by the integer variable LGEN1.
In equation (1.2.3) μa is the apparent dynamic viscosity, that for non-Newtonian fluids, unlike the molecular viscosity of Newtonian fluids, depends not only the fluid's local pressure and temperature, but also on the shear rate γ.
In PHOENICS, it is the apparent kinematic viscosity ENUL, i.e. μ divided by the density, RHO1, which is specified, rather than the dynamic viscosity μ.
Temperature is an important factor in some applications because the apparent viscosity can vary significantly with temperature. For example, some food products thicken to high viscosity when cooled, whereas when heated, some, like honey, pour faster. The temperature dependence of the apparent viscosity is often of the following form:
| μa = μafT | (1.3.1) |
where fT is a temperature function, which has a value of unity in isothermal flow, but otherwise is evaluated from one of the four options defined in Section 4 below.
The following parameters can be set in the Q1 input file to activate whole-field storage and printout to the solution(PHIDA/PHI)file and the RESULT text file:
The formulae used for each of the non-Newtonian models coded in PHOENICS are given in Sections 2 and 3 below.
For a power-law fluid, the apparent dynamic viscosity is given by:
| μa = min( fT K γ(n-1), fT μ0 ) | (2.1.1) |
where K is the fluid consistency index; n is the power-law or flow-behaviour index; μ0 is the zero-shear viscosity limit; and fT is the temperature function.
The power-law model is activated by setting the following parameters in the Q1 input file:
Many purely viscous fluids encountered in polymer processing operations and thermal processing of liquid foods conform to the power-law model. Other examples of power-law fluids include rubber solutions, adhesives, polymer solutions or melts, and biological fluids.
For blood flow, Walburn and Schneck [1976] suggested K=14.67.10-3 Pa.sn and n=0.7755.
The power-law model (for pseudoplastics) predicts infinite (if μ0=0) and zero viscosities in the limits of very low and very high shear rates respectively.
The Sisko fluid model extends the power-law model to include the Newtonian component, and it is defind by:
| μa = [ μ∞ + K γ(n-1) ] fT | (2.2.1) |
where K is the consistency index, μ∞ is the infinite-shear-rate viscosity, and n is an exponent.
The Sisko power-law model may be activated in PHOENICS by introducing the following PIL commands in the Q1 file:
The Sisko model is considered the most appropriate model for lubricating oils and greases, and the model has also been successfully extended to the shear-thinning rheological behavior of concentrated non-Newtonian slurries.
The Cross model is a four-parameter model that covers the entire
shear rate range, and it is given by:
| μa =[ μ∞ + (μo - μ∞ )/( 1+ (λ γ )n ) ] fT | (2.3.1) |
where μo and μ∞ are the asymptotic values of the viscosity at very low and very high shear rates, λ is a time constant in units of s, and n is an exponent.It is clear that this model has finite, non-zero viscosity, at both zero and infinite shear rate limits.
The Cross model may be activated in PHOENICS by introducing the following PIL commands in the Q1 file:
The Cross model has been used extensively for representing the rheological behaviour of many shear-thinning fluids, such as for example foods, biofluids and polymers.
The Carreau model is a four-parameter model given by:
| μa = [ μ∞ + ( μo - μ∞ )( 1 + ( λ γ )2 )(n-1)/2 ] fT | (2.4.1) |
where μo is the zero-shear-rate viscosity, μ∞ is the infinite-shear-rate viscosity, λ is the time constant in s, and n is an exponent.
For (λγ)<<1, the model reduces to a Newtonian fluid with μa=μo. For (λγ)>>1, the model becomes a power-law fluid with μa = K.γ(n-1) .
The Carreau model is activated by setting the following parameters in the Q1 input file:
The Carreau-Yusada model is a five-parameter model given by:
| μa = [ μ∞ + ( μo - μ∞ )( 1 + ( λ γ )m )(n-1)/m ] fT | (2.5.1) |
where again μo is the zero-shear-rate viscosity, μ∞ is the infinite-shear-rate viscosity, λ is the time constant in sm, and m and n are exponents. It is evident that:
The Carreau-Yusada model is activated by setting the following parameters in the Q1 input file:
The Powell-Eyring model is a three-parameter model given by the following expression
| μa = [ μ∞ + ( μo - μ∞) sinh-1 ( λ γ )/( λ γ ) ] fT | (2.6.1) |
where μo is the zero-shear-rate viscosity in Pa.s, μ∞ is the infinite-shear-rate viscosity in Pa.s, and λ is the time constant in s.
The Newtonian model is recovered when μ∞=μo
The Powell-Eyring model is activated by setting the following parameters in the Q1 input file:
Examples of Powell-Eyring fluids are blood, polymer melts, pulp, and suspensions of solids in non-Newtonian liquids.
The Ellis model is a three-parameter model that describes shear-thinning non-yield-stress fluids. The model has seen extensive use for modelling the rheologies of a wide range of polymeric solutions. It can be used as a substitute for the power-law model and is appreciably better than the power-law in matching experimental measurements. Its distinctive feature is the low-shear Newtonian plateau without a high-shear plateau, and hence it may better match the observations in low and medium shear-rate regimes than in the high shear-rate regimes.
| μa = [ μ0/( 1 + ( τ/τ1/2 )α-1 ) ] fT | (2.7.1) |
where μo is the low-shear viscosity, τ is the shear stress, τ1/2 is the shear stress at which μa = μo/2 and α is a flow-behaviour index related to the power-law index by α= 1/n.
The Ellis model is activated by setting the following parameters in the Q1 input file:
For the Bingham-Plastic model, the apparent viscosity is given by:
| μa= [ K + fpτo/γ ] fT | (3.1.1) |
where K is the rigidity coefficient or plastic viscosity, τo is the yield stress, and fp is the Papanastasiou (1987) regularisation function, given by:
| fp = 1 - exp ( -mpγ ) | (3.1.2) |
where by default the stress-growth parameter, mp=1000. In the limit of mp = 0, a Newtonian fluid is recovered, whereas the limit of mp -> ∞ is equivalent to the discontinuous, standard Bingham-plastic model.
The standard Bingham-plastic model uses fp=1, and this is activated by setting the following parameters in the Q1 input file:
The Bingham-Papanastasiou model is selected by setting IENULA=-6, and then fp is computed from equation (3.1.2) above, with a default value of mp=1000s. This value of mp can be changed by setting:
As with the power-law option, the calculation of LGEN1 is activated automatically. When STORE(BTAU) appears in the Q1 file, the local values of τ as defined by equation (3.4) may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT. Similar to the power-law option, the model may also be activated for use in conjugate-heat-transfer applications by setting ENUL=GRND10 and using GRND5 in the "props" file to select the Bingham option for ENUL.
A Bingham fluid is a fluid for which the imposed stress must exceed a critical yield stress to initiate motion. Examples of fluids which behave as, or nearly as, Bingham plastics include water suspensions of clay, dilling mud, sewage sludge, toothpaste, mayonnaise, chocolate, mustard, some emulsions and thickened hydrocarbon greases, and slurries of uranium oxide in nuclear reactors.
The Herschel-Bulkley (HB) model is a three-parameter rheological model that describes the Newtonian and all main classes of the time-independent non-Newtonian fluids. It is given by the relation:
| μa = [ K γ(n-1) + fpτ0/γ ] fT | (3.2.1) |
where K is the consistency coefficient, γ is the shear rate, n is the flow behavior index, τo is the yield-stress above which the substance starts to flow, and fp is the Papanastasiou (1987) regularisation function, given by equation (3.1.2).
The HB model reduces to the power-law when τo=0, to the Bingham plastic when n = 1, and to the Newton’s law for viscous fluids when both these conditions are satisfied. There are six main classes to this model, as follows:
The standard HB model uses fp=1, and it is activated by setting the following parameters in the Q1 input file:
The HB-Papanastasiou model is selected by setting IENULA=-7, and then fp is computed from equation (3.1.2) above, with a default value of mp=1000s. This value of mp can be changed by setting:
The Casson model is a three-parameter model with a yield stress and two rheological parameters. A Casson fluid is a shear-thinning liquid which is assumed to have an infinite viscosity at zero rate of shear, a yield stress below which no flow occurs, and a zero viscosity at infinite shear rate.
| μa = [ ( fp√( τ0/γ )+√K )2] fT | (3.3.1) |
where K is the consistency index, τ0 is the yield stress and fp is the Papanastasiou regularisation function, which is given by:
| fp = 1 - exp ( -√(mpγ ) ) | (3.3.2) |
If the shear stress is smaller than the yield stress, the Casson fluid behaves like a solid, otherwise as a fluid.
This model has often been used for describing the steady-state behaviour of many foostuffs, such as for example yogurt, tomato puree, molten chocolate, etc. The flow behaviour of some particulate suspensions also closely approximate to this type of behaviour.
Biological materials, especially blood, can also be described very well by the Casson fluid model. For the description of blood flow, Charm et al (1964,1965) and Konozsy (2012) suggested the use of the Casson-Papanastasiou model with τ0 = 10.82 mPa and K = 3.1.10-3 Pa.s for mp >100.
The standard Casson model uses fp=1 and it is activated by setting the following parameters in the Q1 input file:
The Casson-Papanastasiou model is selected by setting IENULA=-8, and then fp is computed from equation (3.3.2) above, with a default value of mp=1000s. This value can be changed by setting:
There are four options to account for the dependence of the apparent viscosity on temperature, as follows:
| fT = exp{ β ( 1/T' - 1/T'r ) } | (4.1) |
| fT = exp{ β ( T' - T'r ) } | (4.2) |
| fT = { 1. + β ( T' - T'r ) }-n | (4.3) |
| fT = exp{ β/T' } | (4.4) |
where T'= (T + To), Tr'= (T + Tr), β is a coefficient, Tr is a reference temperature, and To is the PHOENICS reference temperature TEMP0.
Equation (4.3) is applicable to power-law fluids only, and n is the power-law exponent.
In equation (4.4), the coefficient β is often defined by Ea/R where Ea is the activation energy in J/mol and R is the universal gas constant, which is equal to 8.3147732 J/mol.
The following parameters may be set in the Q1 input file:
If IENULB is set to a negative value, then for the Cross, Carreau and Carreau-Yusada models, both the appararent viscosity μa, and the time constant λ, are multiplied by the temperature function, fT.
For highly-viscous fluids, slow convergence or other numerical difficulties may be encountered, such as for example oscillations or even divergence of the field variables. The first course of action might be to increase the relaxation on the moomentum equations and/or introduce linear relaxation of the apparent kinematic viscosity.
If the foregoing practices prove unsuccessful, then it might be beneficial to set ADDDIF=T. This setting modifies the pressure-correction equation in such a way as to accelerate convergence for highly-viscous flows.
As noted earlier, the discontinuity in apparent viscosity at the yield surface of viscoplastic fluids can present numerical difficulties, such as oscillations or even divergence. The simulation can be made more amenable to computation by adopting the regularisation of Papanastasiou (1987), which was described earlier for the Bingham-Plastic, Herschel Bulkley and Casson models.
The Bingham-plastic and Herschel-Bulkley models and their variants have the option to set LSG31=T, which linearises the apparent viscosity in terms of its old iterate values, as follows:
| τ = K γn + τo fp | (5.1) |
| γ* = τ/μa* | (5.2) |
| μa = K γ*(n-1) + τo fp/γ* | (5.3) |
This linearisation usually provides enough implicit relaxation to damp down any oscillations in the field values of μa, but it can lead to very slow convergence.
For the simulation of the turbulent flow of non-Newtonian fluids with wall functions, the user is advised that the use of standard wall functions in these flows is probably questionable, and whilst modified wall functions can be found in the literature ( see for example Sawko (2012)), possibly more accurate results are likely to be obtained via the use of a low-Reynolds-number turbulence model ( see for example Malin (1998) ) or from an enhanced wall-function treatment.
Several Q1's may be found in the PHOENICS input library which demonstrate and validate laminar and turbulent flow of viscous non-Newtonian fluids, as listed below.
Options Library - Non-Newtonian cases at \phoenics\d_earth\d_opt\d_non-newtonian\inplibAll the non-Newtonian expressions described above can be set via entries in the PROPS file. All seven real constants ENULA - ENULG and both integer constants IENULA and IENULB must be specified when the kinematic viscosity flag is set to GRND4.
For example, the lines
<53 > Non-Newtonian Fluid
Non-Newtonian fluid
53 850. GRND4 2100.0 0.15 6.8E-5
23. 0.5 3. 4. 5. 0.010372 240. 4 1
would declare that fluid number 53 had the name 'Non_Newtonian fluid, with a density of 850 kg/m3 and used the Carreau-Yusada model for the viscosity with ENULA=23., ENULB=0.5, ENULC=3., ENULD=4., ENULE=5., ENULF=0.010372, ENULG=240., IENULA=4 and IENULB=1.
Relevant "help-file" entries are :- ENUL Reserved Names.
The FORTRAN coding for the non-Newtonian laws may be found in REAL FUNCTION KINVISL which resides in the file gxknvsl.for.
M R Malin
2022.04.25