DOI: 10.3303/CET24111122 Paper Received: 22 July 2024; Revised: 01 August 2024; Accepted: 30 September 2024 Please cite this article as: Seethao T., Namkanisorn A., Samanmulya T., Wattananusorn S., 2024, A New Viscosity Model for Non–Newtonian Fluids : Part II - Anomalies and Discrepancies from Various Approaches, Chemical Engineering Transactions, 111, 727-732 DOI:10.3303/CET24111122 CHEMICAL ENGINEERING TRANSACTIONS VOL. 111, 2024 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Valerio Cozzani, Bruno Fabiano, Genserik Reniers Copyright © 2024, AIDIC Servizi S.r.l. ISBN 979-12-81206-11-3; ISSN 2283-9216 A New Viscosity Model for Non–Newtonian Fluids : Part II – Anomalies and Discrepancies from Various Approaches Tupthai Seethao*, Apinan Namkanisorn, Thachanan Samanmulya, Santi Wattananusorn Department of Chemical Engineering, School of Engineering, King Mongkut's Institute of Technology Ladkrabang, Bangkok 10520, Thailand tupthai.seethao@gmail.com This current study presents an explanation of shear–thinning fluid behavior outside the shear–thinning regions, where rheological responses typically exhibit like Newtonian fluids. Discussions on the alternative perspectives provided by a newly proposed numerical methodology are introduced. Additionally, the extension of this proposed model is included to offer another possibility for predicting fluid characteristics, especially in particular flowing zones. Previous numerical and experimental data are also unified to provide appropriate analysis and enhance the understanding of fluid flow phenomena in these specific regions. 1. Introduction This paper should be considered a direct continuation of a new viscosity model for non–Newtonian fluids: part I (Seethao T. et al., 2024), in which the characteristics were mathematically described. The development and performance of the proposed model were evaluated using Computational Fluid Dynamics (CFD) techniques, including grid independence analysis, experimental data validation, and comparison with commonly used viscosity models. Here, the same method is applied to demonstrate this mathematical approach, aiming to explain discrepancies between its predictions and other viscosity models, especially in the low and high shear– rate zones. Furthermore, an additional extension of the recently developed model is introduced for particular cases to generate more simulated data for predictions. In general, shear–thinning behavior is usually expressed within a specific range of applied shear rates. The fluid characteristics outside this range (the low and high shear rate regions) tend to be exceptions. The power law viscosity model (Ostwald, 1929) is the simplest and most often used for shear–thinning fluids. However, its performance is limited to capturing only the shear–thinning part of rheological data. Due to this limitation, several other complex constitutive models were later introduced, such as the Cross power law model (Cross M.M., 1965) and the Bird–Carreau model (Bird R.B. and Carreau P.J., 1968). These models can effectively reproduce the rheology of shear–thinning fluids and typically present fluid characteristics outside the shear–thinning region as Newtonian behavior. However, the depiction of “Newtonian” might be one of many conditions that fluid could exhibit. Considering this, the alternative perspectives from the newly proposed rheological model will provide an explanation for shear–thinning fluid behavior, particularly in these specific areas. 2. Numerical Method In this section, the mathematical formulation and the main aspects related to the proposed model are presented. The governing equations for laminar, isothermal, and incompressible flow of shear–thinning fluids are the mass (continuity) and the momentum conservations, respectively: 𝜕𝜕(𝜌𝜌 𝒖𝒖) 𝜕𝜕𝜕𝜕 + 𝛻𝛻 ⋅ (𝜌𝜌𝒖𝒖 𝒖𝒖) = −𝛻𝛻𝛻𝛻 + 𝛻𝛻 ⋅ 𝝉𝝉 (2) ∇ ⋅ (𝒖𝒖) = 0 (1) 727 together with a constitutive equation to describe the relation between the viscosity and shear rate �̇�𝛾 for the fluid of interest. In the above equations 𝒖𝒖 is the velocity vector, 𝜌𝜌 is fluids density, 𝛻𝛻 is the pressure, and 𝝉𝝉 is the stress tensor or shear stress which is given by 𝝉𝝉 = 𝜂𝜂[ ∇𝒖𝒖 + (∇𝒖𝒖)T ] = 2𝜂𝜂𝐃𝐃 (3) where 𝜂𝜂 is the dynamic viscosity of fluid and 𝐃𝐃 is the symmetric rate of strain tensor: 𝐃𝐃 = 1 2 [ ∇𝒖𝒖 + (∇𝒖𝒖)T ]. (4) Table 1 displays the constitutive equations or viscosity models used in this study. The upper part of the table includes well–established and widely employed viscosity models that are presented in numerous rheology publications. The lower part of the table features the proposed viscosity model (Eq(5)), which was presented in part (I) of this study, and the latest version of this model (Eq(6)). This extension has been developed in this part of the study to expand the explanation of fluid characteristics in the high shear–rate zone using different approaches. Table 1: Viscosity models used in this study. Model name Equation Power law (Ostwald–de Waele relation) 𝜂𝜂(�̇�𝛾) = 𝐾𝐾�̇�𝛾𝑛𝑛−1 Sisko (Sisko, A.W., 1958) 𝜂𝜂(�̇�𝛾) = 𝐾𝐾�̇�𝛾𝑛𝑛−1 + 𝜂𝜂∞ Cross power law 𝜂𝜂(�̇�𝛾) = 𝜂𝜂∞ + � (𝜂𝜂0 − 𝜂𝜂∞) 1 + (𝑚𝑚�̇�𝛾)𝑛𝑛 � Bird–Carreau 𝜂𝜂(�̇�𝛾) = 𝜂𝜂∞ + (𝜂𝜂0 − 𝜂𝜂∞)[1 + (𝑘𝑘�̇�𝛾)2] (𝑛𝑛−1) 2 The proposed viscosity model 𝜂𝜂(�̇�𝛾) = Ω 𝜇𝜇0 1 + � 1 𝐾𝐾∗ + (𝜆𝜆|�̇�𝛾|)𝑛𝑛−1� where Ω ≡ 𝐾𝐾∗ + 2 𝐾𝐾∗ + 1 (5) The proposed viscosity model (with its extension) 𝜂𝜂(�̇�𝛾) = Ω 𝜇𝜇0 1 + � 1 𝐾𝐾∗ + (𝜆𝜆|�̇�𝛾|)𝑛𝑛−1� + (1 − 𝜀𝜀) � 1 − �̇�𝛾𝑛𝑛−1 1 + 𝐾𝐾∞∗ � 𝜇𝜇∞ (6) Geometry, flow conditions, and mesh convergence analysis The procedure (and parameters associated with each viscosity model) used in the previous development (part I) is once again applied in this section. A planar channel with an aspect ratio of 100:1 was represented as a schematic of the computational domain to provide fully developed flow at the required flow rates, as shown in Figure 1. In the analysis, the following boundary conditions were applied for velocity: an imposed uniform velocity "U" at the entrance, a Neumann condition (setting the velocity derivative to zero in the flow direction) at the exit, no–slip at the wall, and symmetry at the centerline. For the boundary conditions for pressure in the momentum equation, a pressure value set to zero was also applied at the exit. Figure 1: Schematic of computational domain and boundary conditions. Wall Uniform Velocity Inlet U 100H H 2 x y Symmetry Plane 728 The system was solved by the finite volume method via the opensource OpenFOAM® computational package. The differencing schemes were discretized using a 2nd–order accurate scheme, specifically central differences for the diffusion term and the bounded linear upwind scheme (LUDS) for the convective terms. Coupling of pressure and velocity was achieved through the well–known semi–Implicit Method for Pressure–Linked Equations (SIMPLE), proposed by Patankar (1972). The scaled residuals for both pressure and velocity were observed to reach an asymptotic value of 10–6. For mesh convergence analysis, three consecutive hexahedral meshes, including 600 × 40 (Mesh 1), 1200 × 80 (Mesh 2), and 2400 × 160 (Mesh 3), were applied to determine a suitable mesh density and to investigate the accuracy of simulations. The analysis was carried out by using a Newtonian fluid at an intermediate generalized Reynolds number (Regen = 10.33) [Kozicki et al. (1966) and Nguyen H. & Nguyen N.–D. (2012)], following the same approach as the grid–independence analysis in Poole and Ridley (2007). The analysis also defined a relative error for each flow corresponding parameters (centreline velocity, pressure drop, and the development length) to verify the quantitative comparison among these meshes. Detailed analysis is provided in the first part of this study. From all analyzed data, it was indicated that the computational domain with the density of 1200 × 80 mesh cells (Mesh 2) ensure the desired accuracy with a respectable computational cost and time. 3. Results and discussion 3.1 Comparison with experimental data and various viscosity model The results obtained from a series of CFD simulations using various viscosity models with different Reynolds numbers were conducted to provide viscosity–shear rate relationship as illustrated in Figure 2. The simulated apparent viscosities using the proposed viscosity model (red line, Eq(5)) and this model with an extension (dark blue line, Eq(6)) were compared to previous experimental data for hydroxypropyl guar (HPG) of Guillot and Dunand (1985). These data were also previously simulated by other viscosity models and presented in the work of Moukhtari and Lecampion (2018). However, no experimental data are studied in high shear–rate zone. Figure 2: The viscosity – shear rate curve of HPG fluid from the experimental data and various viscosity models. As it can be seen, the prediction results from each viscosity model can be separately categorized into the shear– thinning part (moderate zone) and the external part of this relation. In the moderate shear–rate zone, prediction 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷ 10⁰ 10⁻¹ 10⁻² 10⁻³ 10⁻⁴ Vi sc os ity μ [P a ∙ s ] Shear Rate γ̇ [1/s] experimental (Guillot and Dunand, 1985) power law power law (truncated) Sisko Bird-Carreau Cross power law Eq(5) Eq(6) low shear–rate moderate shear–rate high shear–rate 729 from all viscosity models effectively demonstrates the shear–thinning behavior of the fluid and show good agreement with the experimental data (with similar simulated time). However, outside this specific region, some discrepancies among the considered models were noticed, including small differences in the low shear–rate zone and two distinct tendencies of the viscosity curve in a region beyond the informed (experimental) data. To clarify these differences, the explanation of various approaches is discussed in the following sections. 3.2 Discrepancies from various approaches in low and high shear–rate zones Low shear–rate zone To observe the variation of calculated viscosity in the low shear–rate zone, a detailed viscosity–shear rate relationship from all models in this zone (except for the power law and Sisko models due to their over–estimate predictions) is presented in Figure 3, along with the experimental data from Guillot and Dunand (1985). It is noted that the defined initial viscosity (𝜂𝜂0) from each viscosity model approaches their critical shear rate �̇�𝛾𝑐𝑐 (the shear rate at which fluid starts to exhibit non–Newtonian behavior) differently in terms of both value and mechanism. The predictions can be categorized into three pathways: Figure 3: Emphasis on the relationship between viscosity and shear rate for various viscosity models within the low shear–rate zone. 1. The initial viscosity gradually decreases from the beginning of the relationship (Cross power law model). 2. The value of initial viscosity remains the same until it reaches the critical point, then begins to decrease (truncated power law and Bird–Carreau model). 3. The value of the given initial viscosity slightly increases before decreasing, as characteristic of shear– thinning fluid, as demonstrated by the proposed model’s prediction. Differences between experimental and numerical results are observed. The unique pathway simulated by the proposed viscosity model may be explained by the responsive mechanism of the microstructure in shear– thinning fluid to the applied shear stress. The expansion of random coils in concentrated shear–thinning fluid can generate a small amount of viscosity before disentanglement. This increasing viscosity is similarly observed not only in Guillot and Dunand's experiment (1985) but also in the experimental investigation by Roumpea et al. (2017). Therefore, it can be stated that the non–linear function of proposed model reasonably illustrates a concealed behavior of shear–thinning fluids in this region. High shear–rate zone In the high shear–rate zone, as depicted in Figure 2, the predicted viscosity from Eq(5) continuously decreases, similar to the prediction made by the power law model. This decrease poses a limitation in numerous fluid flow calculations. However, concerning experimental data, the characteristics of shear–thinning fluid under high shear rate conditions are still evident. These data were reported in Selby and Miiller’s investigation (1995), where they measured the viscosity of multigrade engine oil using matched expansion metals as a wall of the viscometry within high shear rate conditions. According to the evidence from this previous experimental work 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10⁻¹ 10⁰ Vi sc os ity μ [P a ∙ s ] Shear Rate γ̇ [1/s] experimental (Guillot & Dunand, 1985) power law (truncated) Bird–Carreau Cross power law Eq(5) Eq(6) 730 (Selby T.W. and Miiller G.C., 1995), it can be noticed that the behavior of shear–thinning fluid in high shear rate flow is influenced in–between two main factors of the surrounding condition, which can be explained as follows: Wall–bounded effect – Typically, when measuring the viscosity of non–Newtonian fluids, the fluid is contained by a container. The viscous region generated by the presence of solid wall (boundary layer), can interfere with momentum flux transfer as shear rate increases. Consequently, the viscosity magnitude cannot decrease further but gradually approaches a certain value known as “viscosity at an infinite shear rate, 𝜇𝜇∞” Boundary–free effect – In cases where there is no (or flexible) restriction between shear rate and viscosity relation, as in boundary–free or adaptable boundary conditions, the characteristic of shear–thinning fluid can still be indicated in high shear rate flow. This is predicted by the power law model or Eq(5) of the proposed model. Under these conditions, it was determined that the simulated results from Eq(5) required an extension to broaden its capability in specific cases. This extension was carried out following a similar approach to the concept of the Sisko model, involving the addition of another term (corresponding to 𝜇𝜇∞ and its path function to a further infinite shear rate) to the proposed model. This addition was made to ensure that the predictions align with wall–bounded effect conditions, commonly encountered in many fluid flows. The function of 𝜇𝜇∞ in the additional term can be expressed as follows: � 1 − �̇�𝛾𝑛𝑛−1 1 + 𝐾𝐾∞∗ � 𝜇𝜇∞ where 𝐾𝐾∞∗ = � 𝐾𝐾 𝜇𝜇0 − 𝜆𝜆𝑛𝑛−1� �̇�𝛾𝑛𝑛−1. (7) Equations (5) and (7) can be mathematically combined into a more general form (Eq(6) in Table 1) as 𝜂𝜂(�̇�𝛾) = Ω 𝜇𝜇0 1 + � 1 𝐾𝐾∗ + (𝜆𝜆|�̇�𝛾|)𝑛𝑛−1� + (1 − 𝜀𝜀) � 1 − �̇�𝛾𝑛𝑛−1 1 + 𝐾𝐾∞∗ � 𝜇𝜇∞ where 𝜀𝜀 is a constant representing the wall–bounded effect ranging from 0 to 1, dependent on the given flow condition. Under the condition 𝜀𝜀 = 0, a fully wall–bounded effect is realized, whereas 𝜀𝜀 = 1 describes behaviors without the presence of any boundary, approaching an ideal situation. Therefore, choosing a value in the range of 0 < 𝜀𝜀 < 1 is suggested, providing a more realistic depiction for the incompressible flow of non–Newtonian fluids. Figure 4 illustrates the numerical predictions of viscosity for HPG fluid as a function of shear rate after the adaptation of the proposed viscosity model under a fully wall–bounded condition 𝜀𝜀 = 0. As can be seen, the additional term significantly affects the prediction in the high shear–rate zone. The results from modified model also predict the lowest possible value of viscosity compared to those other exiting models. Figure 5 emphasizes the difference between conditions with and without the presence of flow boundaries via Eq(6). Figure 4: Comparison between a fully wall–bounded (𝜀𝜀 = 0) and boundary–free (𝜀𝜀 = 1) flow conditions on the proposed viscosity model. 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷ 10⁰ 10⁻¹ 10⁻² 10⁻³ 10⁻⁴ Vi sc os ity μ [P a ∙ s ] Shear Rate γ̇ [1/s] experimental (Guillot & Dunand, 1985) Eq(6) with ε = 1 Eq(6) with ε = 0 C rit ic al sh ea r r at e, 𝛾𝛾 𝑐𝑐 B re ak th ro ug h sh ea r r at e, 𝛾𝛾 𝑏𝑏 731 The comparison shows that the simulated results in the low to moderate shear–rate zone for these two conditions are identical, even at the critical shear rate (�̇�𝛾𝑐𝑐) of each flow condition. However, when entering the high shear–rate zone, the additional term considerably affects the calculations. The solutions begin to diverge as shear rate continuously increases. For this turning point of fluid behavior, a “breakthrough shear rate” (�̇�𝛾𝑏𝑏), defined as a point at which a given fluid initially senses or “feels” the presence of surrounding walls, is also presented. 4. Conclusion This research investigated the discrepancies generated by various approaches, particularly the newly proposed viscosity model, to explain the characteristics of shear–thinning fluid in the exceptional zones. Using the proposed viscosity model in the low shear–rate zone allows the prediction to demonstrate the increment of fluid viscosity. This behavior can be observed in concentrated shear–thinning fluids, as reported in previous experimental data, and it does not behave as a Newtonian fluid. Moreover, considering the surrounding conditions in the high shear–rate zone, the prediction from the proposed viscosity model can provide a variety of apparent viscosities, depending on the flow boundary conditions (a fully wall–bounded and boundary–free). It could be clarified that perspectives of the recently proposed methodology provide a comprehensive explanation of fluid phenomenon across all shear rate ranges. Nomenclature 𝜌𝜌 – density of fluid, kg/m3 𝒖𝒖 – velocity vector, m/s 𝛻𝛻 – pressure, Pa 𝜂𝜂 – dynamic viscosity of fluid, kg/(m s) 𝝉𝝉 – stress tensor or shear stress, Pa 𝐃𝐃 – symmetric rate of strain tensor, 1/s 𝑛𝑛 – flow behavior index, - 𝐾𝐾 – flow consistency index, Pa sn �̇�𝛾 – shear rate of the fluid flow, 1/s 𝜇𝜇0 – initial viscosity, kg/(m s) 𝜇𝜇∞ – infinite viscosity, kg/(m s) 𝐾𝐾∗ – modified consistency index of 𝐾𝐾, - 𝜆𝜆 – time constant for proposed viscosity model, s 𝜃𝜃 – scale factor for time constant 𝜆𝜆, - Regen – generalized Reynolds number, - 𝑈𝑈 – imposed average velocity, m/s 𝑚𝑚 – time constant for Cross power law model, s 𝑘𝑘 – time constant for Bird–Carreau model, s �̇�𝛾𝑐𝑐 – critical shear rate, 1/s �̇�𝛾𝑏𝑏 – breakthrough shear rate, 1/s References Bird R.B. and Carreau P.J., 1968, A Nonlinear Viscoelastic Model for Polymer Solutions and Melts–I, Chemical Engineering Science, 23, 427 – 434. Cross M.M., 1965, Rheology of Non–Newtonian Fluids: A New Flow Equation for Pseudoplastic Systems, Journal of Colloid Science, 20(5), 417 – 437. Guillot D. and Dunand A., 1985, Rheological Characterization of Fracturing Fluids by Using Laser Anemometry, Society of Petroleum Engineers Journal, 25(1), 39 – 45. Kozichi W., Chou C.H. and Tiu C., 1966, Non–Newtonian Flow in Ducts of Arbitrary Cross–Sectional Shape, Chemical Engineering Science, 21, 665 – 679. Moukhtari F.–E. and Lecampion B., 2018, A Semi–Infinite Hydraulic Fracture Driven by a Shear–Thinning Fluid, Journal of Fluid Mechanics, 838, 573 – 605. Nguyen H. and Nguyen N.–D., 2012. Incompressible Non-Newtonian Fluid Flows, Continuum Mechanics– Progress in Fundamentals and Engineering Applications, InTech, Rijeka, Croatia. Ostwald W., 1929, Ueber die rechnerische Darstellung des Strukturgebietes der Viskosität, Kolloid–Zeitschrift, 47(2), 176 – 187. Patankar S.V. and Spalding D.B., 1972, A Calculation Procedure for Heat, Mass and Momentum Transfer in Three–Dimensional Parabolic Flows, International Journal of Heat and Mass Transfer, 115, 1787 – 1803. Poole R. J. and Ridley B. S., 2007, Development–Length Requirements for Fully Developed Laminar Pipe Flow of Inelastic Non-Newtonian Liquids, Journal of Fluids Engineering, 129, 1281 – 1287. Roumpea E., Chinaud M. and Angeli P., 2017, Experimental Investigations of Non–Newtonian / Newtonian Liquid–Liquid Flows in Microchannels, AIChE Journal, 63, 3599 – 3609. Seethao T., Namkanisorn A. and Wattananusorn S., 2024, A New Viscosity Model for Non–Newtonian Fluids : Part I – Physical Characteristics of Its Mathematical Description, Chemical Engineering Transactions, 109, 619 – 624. Selby, T.W. and Miiller, G.C., 1995, Very High Shear Rate Viscometry, Proceedings of the International Tribology Conference, Yokohama, Japan. Sisko, A.W. 1958, The Flow of Lubricating Greases, Industrial and engineering chemistry, 50(12), 17898 – 1792. 732 378Seethao-PAGATO.pdf A New Viscosity Model for Non–Newtonian Fluids : Part II – Anomalies and Discrepancies from Various Approaches