Continuous flow of a thixotropic elasto-viscoplastic-banded structured fluid: Exponential structure rheological constitutive equation

Ramírez-Torres, L. A.; Herrera-Valencia, E. E.; Soriano-Correa, C.; Sánchez-Villavicencio, M. L.; Campos-Fernández, L.; Ascanio, G.; Hernández-Abad, V. J.; Calderas, F.

doi:10.1063/5.0237760

In this work, the viscoelastic behavior of a complex structured liquid in a continuous squeeze flow is analyzed. This flow is simulated allowing a continuous flow of liquid into the narrow gap between two circular plates though the lower plate. The complex liquid is characterized by the exponential structure rheological (ESR) constitutive equation, which is a generalized exponential thixotropic-elasto-viscoplastic-banded model, chosen for this study due to its capacity to predict non-Newtonian and complex behavior such as (i) shear-thinning, (ii) shear-thickening, (iii) yield stress, (iv) thixotropy, (iv) rheopexy, and (v) shear banding flow behavior. The exponential rheological equation of state is defined by a class of exponential Phan-Thien–Tanner-type models, which includes specific cases documented in the literature. The viscoelastic, kinetic, and structural mechanisms in the ESR model are characterized by the association of non-dimensional numbers to each mechanism. To solve the set of non-linear partial differential equations, a perturbation scheme is suggested, based on a small parameter that represents the ratio between two characteristic lengths. At zeroth order (neglecting the inertial mechanisms of the momentum equation), it is found that the normal force on the upper disk is directly related to shear dependent viscosity (thixotropy, shear-thinning, shear-thickening, yield stress behavior, and concentration effects). At first order, the normal force is related to the effects of the elasticity, and it is parallel to the first normal stress difference associated with the elasticity of the ESR structured fluid.

I. INTRODUCTION

The analysis of fluid–structure interaction (FSI) is an interdisciplinary field that plays a key role in understanding the complex dynamics between fluid and structure.¹

A. Squeeze flow film

One of the most fundamental FSI phenomena in fluid mechanics and rheology is found in the squeeze flow.² The most common system consists of two parallel plates with a radius of r = a, separated by a distance z = h.³ The liquid is placed on the lower plate, and as the upper plate approaches the fluid, which induces radial flow due to the squeezing mechanism.⁴ This compression-driven flow system can have various geometric configurations, such as (i) disks with a radius of r = a and separated by a distance h,³ (ii) spheres separated by a distance h,^5,6 and (iii) a combination of spheres and disks.^3,4 From a rheological perspective, this mechanism is challenging to describe due to the presence of shear strain at the walls, extensional components at the center of the squeeze experiment, vorticity, and inertial mechanisms associated with the acceleration of the upper plate relative to the fluid.^7,8 The study of the squeeze flow film raises three important issues that deserve attention from rheologist worldwide: (A) The geometric configurations of the system,⁷ (B) The experimental and material conditions of the system,⁹ and (C) The mathematical and physical theoretical approaches, including computational methods, theoretical approximation, asymptotic analysis, and numerical methods.^10,11

B. Applications of the squeeze flow film

From an experimental perspective, one of the most important aspects of squeeze flow is the measurement of the material functions.⁷ In this context, the force required to maintain the flow of the liquid flow between the plates depends on the materials properties of the fluid and the velocity of the plate squeezing the flow.¹² Additionally, it has been shown that the inertial mechanisms are minimal, in contrast with the elastic effects associated with the load-bearing capacity of the squeeze flow film.¹³ The squeeze flow film in isothermal conditions for both Newtonian and non-Newtonian fluids has attracted significant interest due to its numerous applications, including (i) electroosmotic pumps,¹⁰ (ii) biorheology,^14,15 (iii) rheometry and rheology devices,¹¹ (v) enhanced oil recovery,¹⁶ (vi) electrorheological fluids,¹⁰ (vii) motors bearing,¹⁷ (viii) lubrication,¹⁸ (ix) fibers,¹⁹ (x) ceramics,²⁰ (xi) thermoplastic composites,²¹ (xii) dental composites resin,²² (xiii) food,²³ and (ix) coagulation system and medicine.²⁴ Remarkable examples of the squeeze flow theory can be found in the formation of foams, where bubble boundaries expand biaxially and shrink in thickness in a similar way as squeezing films (Ref. 7 and references therein). In medicine and biological sciences, the functioning of medical valves, diarthrodial joints, knee prostheses, and biological superfibers, such as spider silk are important.^25,26 Interesting transport phenomena and rheology occur during the formation of the collagenous film in tissue engineering.^27,28 The compression of food between the tongue and palate can be approximated (for some solid food) as a squeeze flow.²⁹ In biomimetics, the formation of the spider silk in its spinneret can be modeled as a double exponential geometry biological extruder.²⁶ In polymer science processing, squeeze flow phenomena are observed in fabrication operation, such as stamping, injection molding, sheet forming, and rheometric measurements using plastometers that incorporate parallel-disk squeeze flow geometry.^30,31

C. Mathematical and physical approaches

From a mathematical perspective, the aforementioned theoretical works involve several mathematical techniques, which include the following: (i) analytical approaches, (ii) regular perturbation techniques, (iii) numerical finite difference schemes, (iv) numerical finite element methods, (v) numerical finite volume methods, (vi) hybrid schemes combining perturbation techniques and numerical approaches, and (vii) computational software (e.g., COMSOL Multiphysics). Phan-Thien³² analyzes the small strain oscillatory squeeze film flow of simple fluids. The fluid is characterized by a functional of the strain history, involving the Cauchy–Green tensor, which reduces to a general integral linear viscoelastic rheological equation of state. The system is solved with a stream function, which reduces the problem to a biharmonic equation. Assuming, non-slip conditions, both analytical and asymptotic solutions were found. Phan-Thien and Tanner¹⁸ studied the lubrication mechanisms in a Boger fluid (nearly Newtonian), characterized by an Oldroyd-B constitutive equation. Assuming isothermal conditions and a lubrication approximation, the coupled partial differential equation system was solved using an exact flow kinematic solution and a variant of Galerkin method (One-Galerkin method). Gartling and Phan-Thien³³ studied a plastic fluid in a geometric in a parallel-plate plastometer configuration. The system was solved using a programing code based on a Galerking-finite element mesh for the numerical simulation. The in-silico experiments were calibrated with experimental data for complex fluid exhibiting yield stress mechanisms. Phan-Thien and Walsh³⁴ investigated similarity solutions in a squeeze flow film using an Oldroyd-B fluid. Assuming inertia-less flow and that the squeeze velocity varies exponentially with time, a similarity solution exists. Additionally, their results showed that there is a critical Weissenberg number above which at least one component of the stresses unbounded with time. Phan-Thien et al.³⁵ examined the squeeze flow film of an ideal elastic liquid using a constitutive equation similar to Oldroyd-B, and their numerical scheme was compared with rheometric data. According with their results, the non-linear model admits an exact solution. The non-linear case was solved using an explicit second-order finite difference scheme. The elastic mechanisms were well described with this theoretical-computation scheme. Phan-Thien et al.³⁶ studied the effect of a viscoelastic fluid with both solvent and polymer contribution, using a modified Phan-Thien-Tanner model with shear-dependent viscosity (Carreau model). The system was solved using the exact flow kinematic function and a boundary element method (BEM). The algebraic equations were solved using Gaussian elimination and predictor-corrector Rung–Kutta numerical methods. The squeeze flow film problem was addressed with a fourth order Runge–Kutta scheme. According to their simulations, load enhancements are related to the overshoots and the materials properties of the system. Phan-Thien and Low³⁷ investigated the squeeze flow of a viscoelastic fluid using a lubrication model. The system was modeled by splitting the solvent and polymer contribution, which allows for a Carreau-type viscosity and stress overshoot a startup of shear flow. This research compared the full numerical solution obtained with boundary element method, to the classical lubrication approximation, described by a fourth-order coupled partial differential equation in time and one spatial coordinate. The lubrication squeeze flow problem was solved using a finite difference scheme, and the results showed that the lubrication is a good alternative to full numerical analysis. Phan-Thien³⁸ studied the oscillatory squeeze flow film of a viscoelastic solid, using a neo-Hookean rubber-like finite deformation and the upper convected Maxwell models. The research, demonstrated that the load can exhibit a significant degree of asymmetry, which is largely due to the rubber-like elasticity in the response, which results in a larger force during the downward phase of the displacement. The system admits an analytical solution, and the general scheme was implemented using a central difference numerical scheme. The non-linear coupled differential equations were solved using explicit leap-frog and Euler methods. Phan-Thien et al.³⁹ analyzed large-amplitude oscillatory squeeze flow data for a complex biological material. The system was characterized by two rheological models: a bi-viscosity Newtonian model and a non-linear Maxwell model. The non-linear Maxwell model better matched the experimental data. The non-linear response of the material in large-amplitude oscillatory flows is primarily attributed to shear thinning.

D. Continuous squeeze flow film

In contrast to traditional squeeze flow systems, various approaches have been investigated to avoid the inertia of the upper plate approaching the fluid.^40,41 Oliver et al.⁴² proposed a new system called “continuous squeeze flow film.” In this system, the fluid is continuously injected at a constant rate into the gap between two parallel circular plates in close proximity, through a uniform array of the holes in the lower plate.⁴³ This system has the advantage of not considering the inertia of the plates, allowing the system to be treated as a steady-state problem. Oliver and Shahidullah⁴⁴ studied both normal and reverse squeeze flow in a continuous flow system with Newtonian and Non-Newtonian liquids.⁴⁵ The fluid is injected through the lower plate, with neither of the plates moving, and the liquid exudes from 1580 uniformly distributed holes in the plate surface. All tests were performed at a temperature of 24 °C.

Waters and Gooden⁴⁶ analyzed the flow in continuous-flow squeeze film from an analytical perspective. The inelastic liquid is characterized by the power-law model, and the system was solved using perturbation techniques at zero and first orders. They found that the force at zeroth and first orders depends on the volumetric flow and material

Properties, with inertial mechanisms contributing to the normal force on the upper plate. Waters and Gooden⁴⁷ studied the flow of a viscoelastic fluid in a continuous squeeze flow film. The liquid is characterized by a non-linear equation Oldroyd-B rheological equation of state.^48,49 The system was solved using a perturbation series at zeroth and first orders, and the normal force was obtained using a numerical scheme. At zeroth order, viscous mechanisms are relevant and depend on the material parameters of the constitutive equation. At first order, elastic (storage) mechanisms are significant, while inertial mechanisms do not contribute to the normal force.

E. Complex fluids

Structure fluids are complex systems that exhibit intriguing steady and unsteady behavior under flow.⁵⁰ These systems are found in various industries and are foundational in biomimetics.^51,52 The physics of these systems is a powerful tool for explaining biological phenomena, such as (i) polypeptides, (ii) cellulose, (iii) blood, (iv) DNA, (v) natural fibers, and (vi) synthetic materials (e.g., molten polymer nanocomposites, surfactant solutions, liquid crystals polymers, bio-composites). These materials often exhibit crystalline order, anisotropy, defects, viscoelasticity, and flow aligning structures.^53–55 Additionally, the mechanical response of the human ear has been studied through flexoelectric membranes embedded in non-Newtonian fluids,⁵⁶ and the effects of the solvent depletion on electrokinetic energy conversion in non-linear viscoelastic fluids.⁵⁷ Thixo-elasto-visco-plastic fluids can be described using constitutive equations that take into account build-up and break-down kinetics by flow.⁵⁶ These models incorporate sequences of self-assembly, structural changes, kinetic processes under flow, and mass transfer.⁵⁴ The characterization and modeling of these phenomena have been extensively studied by various authors using sophisticated rheological equations of state, such as (i) the Pom–Pom model,⁵⁸ (ii) the extended Pom-Pom model by Verbeeten et al.,⁵⁹ which incorporates a non-zero second normal stress difference in simple shear flow to avoid discontinuities at steady-state elongation and the unboundedness of the orientation equation at high strain rates. A sophisticated non-linear model introduced by Giesekus, accounts for mobility mechanisms associated with anisotropic Brownian motion and/or anisotropic hydrodynamic drag of the polymer molecules.⁶⁰ The Phan-Thien Tanner (PTT) and extended-generalized Phan-Thien–Tanner equation (g-PTT) have been studied in simple shear flow and classic inhomogeneous Poiseuille flow by Ferras and Afonso.⁶¹ These models serve as the basis for computer simulations based on relaxation time spectra or transient-network theory.^61,62 Boek et al.⁶³ proposed a model, which separates the contributions of solvent and polymer, consisting of a coupled system of non-linear equations, including the upper-convected Maxwell equation and a kinetic equation for structural breakdown. Building on the Phan-Thien–Tanner (PTT) formalism, Herrera-Valencia et al.⁶⁴ developed a new model based on structural energy changes the induced by flow. This model utilizes predator-prey dynamics and population balance principles biomathematics and bioengineering.⁶⁵ The kernel function is based on structural activation energy to alter microstructure, dependent on thermodynamic viscoelastic dissipation work that breaks down microscopic structures.⁶⁶

F. Wormlike micellar structures

Viscoelastic surfactants are characterized by an entangled network of large, worm-like micelle structures.⁶⁷ These structures break and re-form during flow, exhibiting complex rheological behavior.⁶⁸ Predicting the flow behavior of viscoelastic surfactants using constitutive equations has been a challenging issue.⁵³ These systems display Maxwell-type behavior in small–amplitude oscillatory shear flow and shear stress saturation in steady simple shear, which leads thixotropy and shear banding flow.⁶⁹ In the non-linear viscoelastic regime, elongated micellar solutions also exhibit notable features, such as a stress plateau in steady shear rate accompanied by slow transients to reach steady state⁷⁰ and normal stresses upon inception of flow.⁷¹ The plateau, where the stress becomes independent of shear rate has been interpreted as a purely mechanical instability leading to shear banding.⁷² Viscoelastic surfactants have been used as rheological modifiers in coating process⁷³ and in enhanced oil recovery operations,⁷⁴ particularly in the fracturing of subterranean formations.⁷⁵ Additionaly, oil extraction can be achieved by hydraulically inducing fractures in rock formations.⁷⁶ Water-based fracturing fluids have typically used high molecular weight water-soluble polymers.⁷⁷ Recently, polymer-free fracturing fluids, based on viscoelastic surfactants, have been developed for fracturing underground formations.⁷⁸ Fluids composed of viscoelastic surfactants can enhance fracture conductivity compared to polymer-based fluids.⁷⁹ Hydraulic fractures are characterized by having one dimension the width being very small compared to the lateral dimensions, height, and length.⁸⁰ This characteristic aligns with the lubrication, approximation, which assumes that flow is well approximated by a locally uniform flow between parallel plates separated by the local fracture width.⁸¹

Despite their significant technological applications, it is surprising that the rheological modeling of complex fluids, such as viscoelastic surfactants (e.g., worm-like micellar aqueous solutions), and their complex behavior (e.g., thixotropy, rheopexy, shear-banding flow) in lubrication squeeze films have been scarcely studied in the current literature. These rheological phenomena represent a test for new constitutive equations and enhancing rheometric techniques, which motivates the present investigation.^{56,61–64,70,82–84}

In this regard, the main objectives of this work are

To predict the flow behavior of a continuous squeeze film of a complex liquid modeled by the exponential structure rheological (ESR) constitutive equation of state.⁶⁴
To analyze the non-Newtonian mechanisms and their interactions with the material properties associated with the thixotropic-elasto-viscoplastic mechanisms of the fluid, through dimensionless groups that represent the physical characteristic in the system.
To study the effect of viscoelastic surfactant concentration by using rheometric data of an aqueous worm-like micellar solution (cetyl trimethyl ammonium tosylate) to predict the flow behavior for various micellar concentrations.⁸³

This study is organized as illustrated in Fig. 1.

FIG. 1.

View large Download slide

Flow chart of the paper organization.

Section I provides an introduction to the problem and a review of previous work. Section II covers continuity, momentum transport, and rheological equation of state (ESR model). Section III addresses the physical constraints and problem formulation. Section IV presents the scaling rules and dimensionless groups. Section V presents the regular perturbation technique applied to the ratio between axial and radial characteristic lengths, up to zero and first orders. In Sec. VI, analytical expressions and numerical solutions for radial and axial velocities, volumetric flow rate, pressure gradient, and normal force at the upper plate are provided, as a function of the shear-thinning/thickening, thixotropy/rheopexy, yield-stress, wormlike micellar data. Finally, Sec. VII includes concluding remarks and future work (Figs. 2–4).

FIG. 2.

View large Download slide

Illustration of the activation energy to change the initial configuration (structure), under the effect of an external tangential force.

FIG. 3.

View large Download slide

A schematic representation of the continuous squeeze flow film is shown. In the vicinity of the axes origin, where the flow can be approximated as pressure-driven Poiseuille flow in a narrow slit. The letter (a) shows the lower plates with random holes, where the complex structured liquid is flowing. Letter (b) shows a local approximation of the Poiseuille flow in a slit of gap H.

FIG. 4.

View large Download slide

Schematic representation of the material space for the dimensionless numbers of the system. The vertical axis represents the structural and viscoelastic mechanisms at high shear rates, the horizontal axis represents the structural and viscoelastic mechanisms at high shear rate, and the axis into the page represents viscoelasticity.

II. THEORETICAL EQUATIONS

A. Mass and momentum transfer equations

The conservation of momentum in a continuous medium is generally expressed by⁴⁸

\frac{\partial}{\partial t} ρ + \nabla \cdot (ρ v) = 0 .

(1)

In Eq. , v is the velocity vector field and ρ is the density of the liquid. This equation represents a mass balance without generation. The second equation is the momentum balance equation, which is given by the second Cauchy formalism⁴⁸

\frac{\partial}{\partial t} ρ v + \nabla \cdot [ρ v \otimes v] = \nabla \cdot T + f = \nabla \cdot T + ρ g .

(2)

In Eq. , T is the total stress tensor and f is the body forces, which can be defined as the product of the density and acceleration due to gravitational forces. The second-order stress tensor T can be decomposed into two contributions: the scalar pressure p and the deviatoric viscoelastic stress tensor σ

T = - p I + σ .

(3)

In Eq. , I is the unit tensor given by the following matrix:

I = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) .

(4)

B. Exponential structural rheological constitutive equation (ESR)

The exponential structure rheological constitutive equation (ESR) considers the build-up and the break-down dynamics of a complex fluid structure during flow and has proven to correctly predict the behavior of various complex systems in different flow situations. The ESR model describes the evolution of a structural parameter as a consequence of relaxation and dissipation mechanisms, respectively (Herrera-Valencia et al., 2024).⁵³

Φ σ + \frac{J_{0}}{φ_{\infty}} \overset{\nabla}{σ} {+ λ}_{D} \frac{J_{0}}{φ_{\infty}} f w (Tr (σ \cdot D)) \frac{\partial}{\partial t} \overset{\nabla}{σ} = \frac{2 D}{φ_{\infty}} + \frac{λ_{D}}{φ_{\infty}} f w (Tr (σ \cdot D)) \frac{\partial}{\partial t} 2 D .

(5)

In Eq. , Φ is the exponential kernel inspired by the Phan-Thien–Tanner formalism, σ is the viscoelastic stress tensor, J₀ is the compliance (the inverse of the elastic bulk modulus), φ_∞ is the fluidity at high shear rate, λ_D is the interfusion structure-momentum relaxation time,

\overset{\nabla}{σ}

is the co-deformational time derivative of σ, and D is the shear strain tensor, which is the symmetric part of the spatial velocity gradient tensor, i.e., D =∇V^S. The structural viscoelastic dissipation function Φ can be expressed in the following analytical form:

Φ = Exp [(φ_{r} - 1) f_{w} (Tr (σ \cdot D))] = {\sum_{n = 0}^{\infty} [(φ_{r} - 1) f_{w} (Tr (σ \cdot D))]}^{n} .

(6)

In Eq. , φ_r = φ₀/φ_∞ is a fluidity ratio related to the shear thinning and thickening mechanisms and f_w is an analytical function, which depends on the trace of viscoelastic dissipation work, i.e., Tr (σ ⋅D).

C. Exponential kernel viscoelastic dissipation function

The exponential kernel function, Φ, can be expanded using a Taylor series up to the second order as follows:

Φ ≅ 1 + (φ_{r} - 1) f_{w} (Tr (σ : D)) + \frac{1}{2!} {[(φ_{r} - 1) f_{w} (Tr (σ : D))]}^{2} + O ({[(φ_{r} - 1) f_{w} (Tr (σ : D))]}^{2}) .

(7)

Mathematically, the viscoelastic dissipation function fw can be defined such as

f w (Tr (σ \cdot D)) = Exp (- β Tr (σ \cdot D)) .

(8)

Equation can be interpreted as an Arrhenius equation associated with the minimum energy required to change the structure by flow. The parameter β(II_D) is associated with the inverse of the irreversible work needed, to change the structure due to the effect an external tangential force:

β ({I I}_{D}) {= β}_{0} Exp [λ_{B} ({I I}_{D})] = β_{0} \sum_{n = 0}^{\infty} \frac{{(λ_{B} ({I I}_{D}))}^{n}}{n!} .

(9)

In Eq. , the parameter β₀ is the irreversible work in the absence of flow, and λ_B is an intensive shear banding parameter. The upper convective derivative of the stress tensor σ is given in the following analytical form:

\overset{\nabla}{σ} = \frac{D}{Dt} σ - \nabla \otimes v^{T} \cdot σ - σ \cdot \nabla \otimes v .

(10)

In Eq. , (⋅)^T is the transpose of the spatial velocity gradient tensor, and the material time derivative of the shear stress Dσ/Dt can be defined in the following analytical form:

\frac{D}{Dt} σ = \frac{\partial σ}{\partial t} + v \cdot \nabla \otimes σ .

(11)

The shear rate tensor D [Eq. ] is defined as

D = \frac{\nabla \otimes v + \nabla \otimes v^{T}}{2} .

(12)

The second invariant of the shear strain, defined as II_D [Eq. ], is given by

{II}_{D} = \sqrt{2 Tr (D \cdot D)} .

(13)

Here, the reference state is defined as the high shear rate fluidity φ_∞ (inverse of the high shear rate viscosity). D is the shear rate tensor which is the symmetric part of the velocity gradient tensor, i.e., D = ∇v^S. II_D is the second invariant of the rate of deformation tensor, and J₀ is the compliance in the viscoelastic regime, which can be measured trough the inverse of the elastic modulus, i.e., J₀=1/G₀. In Eq. , {φ_r, φ_∞} are the fluidities at zero and very high shear rates, respectively, λ_D is an inter-diffusion momentum structural relaxation time; β₀ is the destruction rate constant in the absence of shear flow; and λ_B is the shear banding intensity parameter. All these phenomenological parameters can be calculated from steady and unsteady rheological experiments (Calderas et al., 2009).⁸⁵ Note that the structural parameter in this model is the fluidity function rather than the viscosity function. This is because the model includes yield stress mechanisms as a particular case (η → ∞, φ → 0). From a computational standpoint, using fluidity as the structural parameter is more efficient. Equations represent closed set of time evolution equations for all independent variables chosen to describe the behavior of complex fluids. Note that these equations are coupled to each other.

D. Particular cases

1. Small deformations of ESR model

In the regime of small deformations, the non-linear rheological equation of state reduces to the classical Burger's formalism, given by

σ = O_{η}^{B} (\frac{\partial}{\partial t}) 2 D,

(14)

where the viscosity operator is given by

O_{η}^{B} (\frac{\partial}{\partial t}) = \frac{1}{φ_{0}} \frac{1 {+ λ}_{D} \frac{\partial}{\partial t}}{1 + \frac{J_{0}}{φ_{0}} \frac{\partial}{\partial t} {+ λ}_{D} \frac{J_{0}}{φ_{0}} \frac{\partial^{2}}{\partial t^{2}}} .

(15)

2. Moderate and high shear strain deformations

Next, we consider the steady-state of Eqs. , which are given by the following expression:

Exp [(φ_{r} - 1) Exp (- β Tr (σ \cdot D))] σ = \frac{2}{φ_{\infty}} D + \frac{J_{0}}{φ_{\infty}} f w (Tr (σ \cdot D)) [- v \cdot \nabla \otimes σ + \nabla \otimes v^{T} \cdot σ + σ \cdot \nabla \otimes v] .

(16)

In steady-state and homogeneous flow, the system can be expressed as

Exp [(φ_{r} - 1) Exp (- β Tr (σ \cdot D))] σ = \frac{2}{φ_{\infty}} D + \frac{J_{0}}{φ_{\infty}} f w (Tr (σ \cdot D)) [\nabla \otimes v^{T} \cdot σ + σ \cdot \nabla \otimes v] .

(17)

In steady-state, homogeneous flow and neglecting the first normal stress contribution (inelastic regimen), the following non-linear equation is obtained:

Exp [(φ_{r} - 1) Exp (- β Tr (σ \cdot D))] σ = \frac{2}{φ_{\infty}} D .

(18)

Equation can be expanded in a Taylor series and the following non-lineal expression can be obtained as follows:

[φ_{r} + Exp (β Tr (σ \cdot D)) - 1] σ = Exp (β Tr (σ \cdot D)) \frac{2}{φ_{\infty}} D .

(19)

In the case where the fluidity ratio φ_r =1, the model reduces to the Newtonian equation, given by

σ = \frac{2}{φ_{0}} D .

(20)

The ESR model was selected for this study due to its ability to predict non-Newtonian phenomena, such as (i) shear thinning, (ii) shear thickening, (iii) yield stress, (iv) thixotropy, (iv) rheopexy, and (v) shear banding flow in the structured fluid (e.g., polymer melts, worm-like micellar solutions, triblock copolymers, biological fluids, dispersions of lamellar liquid crystals, and associative polymers). Additionally, analytical solutions for simple shear flow, as well as pulsating and oscillating flows, can be obtained using this model due to its simplicity compared to more complex models. Furthermore, all the material properties of the simplest ESR model (φ_r, φ_∞, J₀, λ_D, and β₀) are related to the fluid properties and can be estimated from independent rheological experiments in steady and unsteady flows. The fluidities at lower and upper shear rates (φ₀, φ_∞) can be estimated through experiments in steady shear flow. The structural time and elastic modulus (J₀ = 1/G₀, β₀) can be calculated using linear oscillatory flow.

III. PROBLEM FORMULATION

The physical system consists of two disks with a radius r = a and a gap z = h, as shown schematically in Fig. 1.^46,47 The lower disk has an array of holes uniformly distributed. The structured fluid (dilute worm-like micellar solution) flows from the lower to the upper plate, followed by a radial flow. The key objective is to calculate the normal force by the complex fluid on the upper plate.

To model this problem, the following assumptions are made:

It is assumed that the system is in a stationary state, meaning that the time derivatives of the kinematic and dynamic variables are set to zero, i.e., ∂v/∂t = ∂D/∂t = ∂σ/∂t = 0.
Gravitational forces are neglected compared to other forces, i.e., g ≅ 0.
The structured fluid is considered incompressible, i.e. the density does not vary with position or time. As a consequence, the material derivative of the density is zero, i.e., Dρ/Dt = 0⇒ ∇ ⋅ v = 0.
To describe the geometry of the flow system, cylindrical coordinates (r, θ, z) are chosen, defined with respect to an origin at the center of the lower disk.
Angular symmetry is assumed, meaning the angular derivatives of all kinematics and dynamic variables are set to zero, i.e., ∂v/∂θ = ∂D/∂θ = ∂σ/∂θ = 0.
Relaxation mechanisms are small in dilute and semi-dilute complex liquids, i.e., 1 + λ_D f(Tr (σ ⋅ D))∂/ ∂t ≅ 1
There is not mass-flux contribution. i.e., J = 0.
There are not instabilities associated to the shear banding flow, λ_B = 0
The liquid is characterized by the linearized ESR rheological equation of state.

The kinematic and dynamic tensors in the system are (i) spatial field vector, (ii) spatial velocity gradient tensor, and (iii) viscoelastic stress tensor. These mathematical objects can be represented in the following matrix form:

(21)

v (r, z) = (\begin{matrix} v_{r} (r, z) \\ 0 \\ v_{z} (r, z) \end{matrix}),

(21a)

\nabla \otimes v (r, z) = (\begin{matrix} \frac{\partial v_{r}}{\partial r} & 0 & \frac{\partial v_{r}}{\partial z} \\ 0 & v_{r} / r & 0 \\ \frac{\partial v_{z}}{\partial r} & 0 & \frac{\partial v_{z}}{\partial z} \end{matrix}),

(21b)

σ = (\begin{matrix} σ_{rr} & 0 & σ_{rz} \\ 0 & σ_{θ θ} & 0 \\ σ_{zr} & 0 & σ_{zz} \end{matrix}) .

(21c)

In Eq. , v_r and v_z represent the physical components of the velocity field r and z directions, respectively. It is assumed that the flow has axial symmetry with respect to r = 0. The motion is caused by the continuous flow of liquid, at a constant volumetric rate Q = 2πa ⟨V⟩ h through the lower disk into the gap h between the disks.

The boundary conditions can be written as

v_{r} (r, z) = 0 Q = \int_{0}^{2 π} \int_{0}^{a} vz (r, z) rdrd θ on z = 0,

(22a)

v_{r} (r, z) = v_{z} (r, z) on z = h,

(22b)

v_{r} (r, z) = 0 on r = 0,

(22c)

and the normal force F at the top plate is given by

F = - {\int_{S} (T \cdot e_{z}) |}_{z = h} \cdot d S = - {\int_{S} (T \cdot e_{z}) |}_{z = h} \cdot e_{z} dA = - \int_{0}^{2 π} {\int_{0}^{a} (T_{zz}) |}_{z = h} rdrd θ .

(23)

IV. SCALING RULES

A. Non-dimensional variables

To simplify the momentum and rheological equations, the following dimensionless variables for the axial and radial coordinates, radial and axial velocities, components of shear-stress, pressure gradient, total stress, fluidity, and second invariant of the shear-strain are introduced as follows:

(24)

z = \frac{z}{h},

(24a)

r = \frac{r}{a},

(24b)

v = \frac{v_{r}}{〈 V 〉},

(24c)

w = \frac{v_{z}}{(\frac{h}{a}) 〈 V 〉},

(24d)

σ_{(ik)} {= J}_{0} σ_{(ik)},

(24e)

p = \frac{p}{(\frac{a}{h}) \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h}},

(24f)

T_{(ik)} = \frac{T_{(ik)}}{(\frac{a}{h}) \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h}} φ = \frac{φ}{φ_{\infty}},

(24g)

I I_{D} = \frac{I I_{D}}{{(〈 V 〉 / h)}^{2}},

(24h)

where ⟨V⟩ is a characteristic velocity in the system, and it will be discussed in Sec. IV B.

B. Mass conservation

The flow inlet through the lower plate is given by the following expression:

Q_{r} = {\int_{0}^{2 π} \int_{0}^{h} v_{r} (r, z) \cdot r |}_{r = a} dzd θ = 2 π ah \frac{\int_{0}^{h} v_{r} dz}{\int_{0}^{h} dz} = 2 π ah \frac{\int_{0}^{1} v_{r} dz}{\int_{0}^{1} dz} = 2 π ah 〈 v_{r} (z, a) 〉 .

(25)

The volumetric flow entering the lower plate at z = 0 is defined as

Q_{z} = \int_{0}^{2 π} \int_{0}^{a} {v_{z} |}_{z = 0} rdrd θ = \int_{0}^{2 π} \int_{0}^{a} v_{z} (r, 0) {rdrd θ = π a}^{2} \frac{\int_{0}^{1} v_{z} (r, 0) (\frac{r}{a}) d (\frac{r}{a})}{\int_{0}^{1} (\frac{r}{a}) d (\frac{r}{a})} = {π a}^{2} 〈 v_{z} (r, 0) 〉 .

(26a)

By considering mass conservation, the following expression is obtained, Q_r = Q_z:

〈 v_{z} (r, 0) 〉 = 2 \frac{h}{a} 〈 v_{r} (a, z) 〉 ≪ 〈 v_{r} (a, z) 〉 .

(26b)

Equation indicates that the axial velocity is lower than the radial velocity. From Eq. , the following expression is derived:

Q_{r} = 2 π ah 〈 v_{r} (a, z) 〉 = {π a}^{2} 〈 v_{z} (r, 0) 〉 = {π a}^{2} \frac{\int_{0}^{1} v_{z} (r, 0) (\frac{r}{a}) d (\frac{r}{a})}{\int_{0}^{1} (\frac{r}{a}) d (\frac{r}{a})} .

(26c)

If the distribution of the volumetric flow rate in the lower plate follows a power low, i.e., v_z (r,0) = A(r/a)^m, with A positive real constant, A ∈ R⁺ and m ≥ 0, then, by substituting the axial velocity into the axial average integral, the constant A is given by

A = \frac{Q_{r}}{{2 π a}^{2}} (m + 2) .

(27)

So, the axial velocity in the lower plate (z = 0) can be approximated by the following analytical expression:

〈 v_{z} (r, 0) 〉 = \frac{Q_{r}}{2 {π a}^{2}} \underset{Non - homogeneity}{\underset{︸}{(m + 2) {(\frac{r}{a})}^{m}}} = 〈 V 〉 u (r; m) .

(28)

The axial velocity given in Eq. is expressed in terms of the radial volumetric flow Q_r. It is possible to obtain an expression for the average radial velocity, i.e.,

〈 v_{r} (a, z) 〉 = \frac{〈 v_{z} (r, 0) 〉}{h / a} = \frac{Q_{r}}{2 π ah} \underset{Non - homogeneity}{\underset{︸}{(m + 2) {(\frac{r}{a})}^{m}}} = 〈 V 〉 u (r; m) .

(29)

So, the characteristic average velocity ⟨V⟩ is given by the following expression, i.e., ⟨V⟩ = Q_r/2πah, and the non-homogeneity mechanisms are describe by the following analytical function u (r) = ⟨V⟩ (r/a)^m. It is important to note that the flow has the following constraint:

\int_{0}^{1} u (r) rdr = 1 .

(30)

Note that the component of the stress tensor was made dimensionless using the bulk-elastic mechanisms, which allows analyzing the case of a yield stress liquid in the ESR model.

C. Non-dimensional groups

Using the definitions from Eq. , the dimensionless components of the continuity, the momentum equation, the constitutive equations, and the normal force are obtained. Additionally, the following non-dimensional groups are defined:

α = \frac{h}{a} {≪ 1; Re}_{\infty} = ρ h 〈 V 〉 φ_{\infty}; {We}_{\infty} = \frac{J_{0}}{φ_{\infty}} \frac{〈 V 〉}{h}; {Jr}_{\infty} = \frac{\sqrt{β_{0} φ_{\infty}}}{J_{0}}; φ_{r} = \frac{φ_{0}}{φ_{\infty}}; {Ca = Re}_{\infty} {We}_{\infty} = ρ {〈 V 〉}^{2} J_{0} .

(31)

The first group α = h/a is the ratio of two characteristic lengths, the axial and radial lengths. It is important to note that α is lower than the unity, and for this reason, α is the perturbation parameter. The second group is the infinity Reynolds number Re_∞, which relates the inertial and high viscous mechanisms in the fluid. When Re_∞ ≫ 1, the system is dominated by the inertial forces over the dissipation, whereas Re_∞ ≪ 1, the opposite behavior is observed, with dissipation dominating over inertial contributions. The third group is the infinite Weissenberg number We_∞, which measures the viscoelasticity in the system at high shear rates. When We_∞ is less than one, i.e., We_∞ ≪ 1, the system is dominated by dissipation and structural mechanisms. However, if We_∞ ≫ 1, the storage processes control the system through elasticity. The fourth group J_r∞ is the ratio between two compliances associated with the yield stress mechanisms induced by the activation energy to change the structure and the fluidity at high shear rates, and the bulk-elastic modulus. It is important to note that the yield stress value is given by σ_y = (β₀φ_∞)^−1/2. This important number J_r∞ is related to the Bingham (Bi) and Papanastasiou numbers (Pa), trough the following analytical expression:

J_{r \infty} = \frac{\sqrt{β_{0} φ_{\infty}}}{J_{0}} = \frac{Jy}{J_{0}} = \frac{1 / J_{0}}{1 / \sqrt{β_{0} φ_{\infty}}} = \frac{σ_{E}}{σ_{Y}} = \frac{1}{Bi} = Pa .

(32a)

Note that the Bingam number can be interpreted as a plastic number associated with the yield stress mechanisms. In particular, we suppose that the J_r∞ number provides a better physical interpretation for structured fluids, as it leads to a better understanding of the resistance in the system associated with activation-structure, high-viscosity, and bulk elastic mechanisms. The fifth number φ_r measures shear thinning (φ_r < 1) and shear thickening (φ_r > 1) behavior. The last number is the Cauchy number Ca, which is the product between of the infinite Weissenberg and Reynolds numbers, i.e., Ca = Re_∞We_∞,

{Ca = Re}_{\infty} {We}_{\infty} = \frac{J_{0}}{1 / ρ {〈 V 〉}^{2}} = \frac{J_{0}}{J_{E}},

(32b)

when Ca is less than one, i.e., Ca ≪ 1, the elastic bulk compliance is small in compared to the kinetic energy compliance, whereas Ca ≫1 indicates that the system is dominated by elastic compliance and momentum compliance is the dominant mechanisms.

In terms of structure, shear thinning behavior can be conceived as a structure breakdown process, where the system transitions from higher to lower structure states, while shear thickening is the opposite mechanism. The three dimensionless groups [J_r∞, We_∞, φ_r] form a material parametric 3D space.

D. Continuity and momentum components

Substituting the dimensionless variables into the continuity equation, components of the momentum equation, second invariant of the rate of deformation tensor and components of the ESR rheology equations of state, the following equations are obtained:

\frac{\partial v_{z}}{\partial z} + \frac{1}{r} \frac{\partial}{\partial r} ({rv}_{r}) = 0,

(33)

α Ca (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z}) v_{r} = \frac{α}{r} (\frac{\partial}{\partial r} ({r σ}_{(rr)}) {- σ}_{(θ θ)}) + \frac{\partial σ_{(rz)}}{\partial z} {- We}_{\infty} \frac{\partial p}{\partial r},

(34)

α^{3} Ca (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z}) v_{z} = α (\frac{α}{r} \frac{\partial}{\partial r} ({r σ}_{(rz)}) + \frac{\partial σ_{(zz)}}{\partial z}) {- We}_{\infty} \frac{\partial p}{\partial z},

(35)

{I I}_{D} = {[{(\frac{\partial v_{r}}{\partial z})}^{2} {+ 2 α}^{2} {{(\frac{\partial v_{r}}{\partial r})}^{2} + (\frac{\partial v_{r}}{\partial z}) (\frac{\partial v_{z}}{\partial r}) + {(\frac{\partial v_{z}}{\partial z})}^{2}} + α^{4} ({(\frac{\partial v_{z}}{\partial r})}^{2} + \frac{1}{2} {(\frac{v_{r}}{r})}^{2})]}^{\frac{1}{2}} .

(36)

E. Total stress component

The components of the non-dimensional stress and are given by Eq. , and after using Eq. , they become

{We}_{\infty} T_{(zr)} = σ_{(zr)},

(37a)

{We}_{\infty} T_{(rr)} = - {We}_{\infty} p + α σ_{(rr)},

(37b)

{We}_{\infty} T_{(zz)} = - {We}_{\infty} p + α σ_{(zz)},

(37c)

{We}_{\infty} T_{(θ θ)} = - {We}_{\infty} p + α σ_{(θ θ)} .

(37d)

Similarly, the non-dimensional of the ESR constitutive equation can be obtained in the same way.

F. Rheological components

Substitution of dimensionless variables into the ESR model yields the following non-dimensional coupled non-linear partial differential equations:

\begin{matrix} Φ σ_{(rz)} {+ α We}_{\infty} (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z}) σ_{(rz)} {+ We}_{\infty} ({α σ}_{(rz)} \frac{\partial}{\partial r} {+ σ}_{(zz)} \frac{\partial}{\partial z}) v_{r} \\ {- α We}_{\infty} (σ_{(rr)} \frac{\partial}{\partial r} {+ σ}_{(rz)} \frac{\partial}{\partial z}) v_{z} {= We}_{\infty} (α^{2} \frac{\partial v_{z}}{\partial r} + \frac{\partial v_{r}}{\partial z}), \end{matrix}

(38)

Φ σ_{(rr)} + {α We}_{\infty} (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z}) σ_{(rr)} - {2 We}_{\infty} ({α σ}_{(rr)} \frac{\partial}{\partial r} {+ σ}_{(zr)} \frac{\partial}{\partial z}) v_{r} {= 2 We}_{\infty} α \frac{\partial v_{r}}{\partial r},

(39)

Φ σ_{(θ θ)} {+ α We}_{\infty} (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z} - 2 \frac{v_{r}}{r}) σ_{(θ θ)} {= 2 We}_{\infty} α \frac{v_{r}}{r},

(40)

Φ σ_{(zz)} {+ α We}_{\infty} (v_{r} \frac{\partial}{\partial r} {+ v}_{z} \frac{\partial}{\partial z} - 2 \frac{\partial v_{z}}{\partial z}) σ_{(zz)} {- 2 We}_{\infty} (α^{2} \frac{\partial v_{z}}{\partial r} σ_{(rz)}) {= 2 We}_{\infty} α \frac{\partial v_{z}}{\partial z},

(41)

(42)

Φ = Exp [(φ_{r} - 1) f w (Tr (σ \cdot D))],

(42a)

f w (Tr (σ \cdot D)) = Exp (- J_{r \infty}^{2} {We}_{\infty} Tr (σ \cdot D)),

(42b)

Tr (σ \cdot D) = σ_{(zr)} \frac{\partial v_{r}}{\partial r} + α (σ_{(rr)} \frac{\partial v_{r}}{\partial r} {+ σ}_{(θ θ)} \frac{v_{r}}{r} {+ σ}_{(zz)} \frac{\partial v_{z}}{\partial r}) {+ α}^{2} (σ_{(zr)} \frac{\partial v_{z}}{\partial r}) .

(43)

The boundary conditions can now be expressed as

v_{r} (r, z) = 0; v_{z} (z, r) = u (r) on z = 0,

(44a)

v_{r} (r, z) = v_{z} (r, z) = 0 on z = 1,

(44b)

v_{r} (r, z) = v_{z} (r, z) = 0 on r = 0,

(44c)

where u(r) is constant when the input velocity is uniform, but generally satisfies:

\int_{0}^{1} ru (r) dr = 1 .

(45)

Theoretical and experimental results indicate that the Cauchy number (Ca = Re_∞We_∞) given in Eqs. and is sufficiently small; therefore, these terms are neglected.^46,47

V. ASYMPTOTIC ANALYSIS

A. Regular perturbation scheme

It is evident that Eqs. and Eqs. and depend on the perturbation parameter α. Our primary hypothesis is that for a restricted range of αCa = αRe_∞We_∞, it is possible to expand the components of the velocity field, pressure, and components of the deviatoric stress tensor in terms the small parameter α, represented by the following power series:

v_{r} (r, z) = α^{0} v_{r 0} (r, z) + α^{1} v_{r 1} (r, z) + \cdot \cdot \cdot + α^{j} v_{r j} (r, z) + \cdot \cdot \cdot = \sum_{j = 1}^{\infty} α^{j} v_{r j} (r, z),

(46)

v_{z} (r, z) = α^{0} v_{z 0} (r, z) + α^{1} v_{z 1} (r, z) + \cdot \cdot \cdot + α^{j} v_{z j} (r, z) + \cdot \cdot \cdot = \sum_{j = 1}^{\infty} α^{j} v_{z j} (r, z),

(47)

p (r) {= α}^{0} p_{0} (r) + α^{1} p_{1} (r) + \cdot \cdot \cdot + α^{j} p_{j} (r) + \cdot \cdot \cdot = \sum_{j = 1}^{\infty} α^{j} p_{j} (r),

(48)

σ_{(ik) j} (r, z) = α^{0} σ_{(ik) 0} (r, z) + α^{1} σ_{(ik) 1} (r, z) + α^{2} σ_{(ik) 2} (r, z) + \cdot \cdot \cdot + = \sum_{j = 0}^{\infty} α^{j} σ_{(ik) j} (r, z) .

(49)

In Eq. , the indices i, k correspond to the components: [rz, zr, θθ, zz].

B. Taylor series and linear approximation

Using the Taylor theorem and the power series of the radial component, the components of the stress tensor can be expressed in the following form (Herrera et al., 2009, 2010),^86,87 where the shear strain derivatives are given by

σ_{rz} ({We}_{\infty} \frac{\partial v_{r}}{\partial z}) = σ_{rz} ({We}_{\infty} \frac{\partial v_{r 0}}{\partial z}) + ({We}_{\infty} \frac{\partial v_{r}}{\partial z} {- We}_{\infty} \frac{\partial v_{r 0}}{\partial z}) {\frac{\partial σ_{rz}}{\partial ({We}_{\infty} \frac{\partial v_{r}}{\partial z})} |}_{{We}_{\infty} \frac{\partial v_{r}}{\partial z} \to {We}_{\infty} \frac{\partial v_{r 0}}{\partial z}} + Rn .

(50a)

The linear approximation approach is

{We}_{\infty} \frac{\partial v_{r}}{\partial z} = {We}_{\infty} \frac{\partial v_{r 0}}{\partial z} + {α We}_{\infty} \frac{\partial v_{r 1}}{\partial z} .

(50b)

Substituting Eq. into Eq. , we obtain

σ_{rz} ({We}_{\infty} \frac{\partial v_{r}}{\partial z}) = σ_{rz} ({We}_{\infty} \frac{\partial v_{r 0}}{\partial z}) α^{0} + {We}_{\infty} \frac{\partial v_{r 1}}{\partial z} \frac{\partial σ_{rz}}{\partial ({We}_{\infty} \frac{\partial v_{r 0}}{\partial z})} α^{1} + Ο (α^{2}) .

(50c)

Equating the perturbation solution, with the Taylor approximation contribution

σ_{(rz) 0} = σ_{rz} ({We}_{\infty} \frac{\partial v_{r 0}}{\partial z}),

(50d)

σ_{(rz) 1} = {\frac{\partial σ_{rz}}{\partial ({We}_{\infty} \frac{\partial v_{r}}{\partial z})} |}_{{We}_{\infty} \frac{\partial v_{r}}{\partial z} \to {We}_{\infty} \frac{\partial v_{r 0}}{\partial z}} = \frac{\partial σ_{(rz) 0}}{\partial ({We}_{\infty} \frac{\partial v_{r 0}}{\partial z})} .

(50e)

Notice that the zeroth order in the parameter α⁰ corresponds to the steady-state and homogeneous case of the ESR constitutive equation.⁶⁴

It is important to highlight the following points:

The zeroth order corresponds to the steady-state and Poiseuille flow.
The first order in the parameter α¹ represents the local linear approximation, incorporating the elastic mechanisms through the first normal stress difference.

In Sec. V C and subsections V C 1–V C 3 the results from this analysis will be used to determine the zeroth and first-order terms for the normal force.

C. Zeroth-order theory O (α⁰)

Upon substituting Eq. into Eqs. and neglecting terms of order α and higher, the equations of continuity and the radial and axial components of the momentum equation are given by

\frac{\partial v_{z 0}}{\partial z} + \frac{1}{r} \frac{\partial}{\partial r} ({rv}_{r 0}) = 0,

(51)

\frac{\partial σ_{(rz) 0}}{\partial z} {= We}_{\infty} \frac{\partial p_{0}}{\partial r},

(52)

\frac{\partial p_{0}}{\partial z} = 0 .

(53)

The components of the ESR constitutive equation are given in the following form:

Φ_{0} σ_{(rz) 0} {+ σ}_{(zz) 0} X = X,

(54)

Φ_{0} σ_{(rr) 0} {= 2 σ}_{(zr) 0} X,

(55)

Φ_{0} σ_{(θ θ) 0} = 0,

(56)

Φ_{0} σ_{(zz) 0} = 0,

(57)

Φ_{0} = Exp [(φ_{r} - 1) f_{w} (Tr {(σ \cdot D)}_{0})],

(58)

f_{w} (Tr {(σ \cdot D)}_{0}) {= Exp [- J}_{r \infty}^{2} Tr {(σ \cdot D)}_{0}],

(59)

Tr {(σ \cdot D)}_{0} = σ_{(zr) 0} X .

(60)

In Eqs. , , and , X is the dimensionless shear stress, introduced as a mathematical assumption to solve the equations. It can be defined as

X (r, z) {= We}_{\infty} \frac{\partial v_{r 0}}{\partial z} .

(61)

Equations represent the usual lubrication approximation equations.^46,47 Equation implies that p₀ is independent of z. The zero-order boundary condition are given by

v_{r 0} (r, z) = 0; v_{z 0} (r, z) = u (r) on z = 0,

(62a)

v_{r 0} (r, z) = v_{z 0} (r, z) = 0 on z = 1,

(62b)

v_{r 0} (r, z) = v_{z 0} (r, z) = 0 on r = 0,

(62c)

and the constraint of the radial velocity and volumetric flow rate is given by the following normalized integral:

\int_{0}^{1} ru (r) dr = 1 .

(63)

1. Material functions

To find analytical expressions for the material properties of the system, Eqs. , , and are combined.

Exp [(φ_{r} - 1) {Exp [- J}_{r \infty}^{2} σ_{(zr) 0} X]] σ_{(rz) 0} = X .

(64)

Equation can be developed in a Taylor series in the following analytical form:

({Exp [J}_{r \infty}^{2} σ_{(zr) 0} {X] + φ}_{r} - 1) σ_{(rz) 0} {= XExp [J}_{r \infty}^{2} σ_{(zr) 0} X] .

(65)

Expanding can be developed in a Taylor series in the following analytical form:

J_{r \infty}^{2} X {(σ_{(zr) 0})}^{2} + (φ_{r} {- J}_{r \infty}^{2} X^{2}) σ_{(rz) 0} - X = 0 .

(66)

Notice that Eq. is quadratic in the shear stress and shear strain. The non-zero principal extra stresses, given by Eqs. , to zeroth order in α, are

σ_{(rz) 0} (X) = \frac{X}{Φ_{0} (X)},

(67)

σ_{(rr) 0} (X) = 2 {(σ_{(rz) 0} (X))}^{2} = 2 {(\frac{X}{Φ_{0} (X)})}^{2} .

(68)

The material functions viscosity and the first normal stress difference can be expressed as:

\frac{1}{Φ_{0} (X)} = \frac{\frac{J_{r \infty}^{2}}{φ_{r}} X^{2} - 1 + \sqrt{{(\frac{J_{r \infty}^{2}}{φ_{r}} X^{2} - 1)}^{2} + 4 \frac{J_{r \infty}^{2}}{φ_{r}^{2}} X^{2}}}{2 \frac{J_{r \infty}^{2}}{φ_{r}} X^{2}},

(69)

Ψ_{0} (X) = 2 {(\frac{1}{Φ_{0} (X)})}^{2} .

(70)

2. Radial component of the field velocity

Integration of Eq. with respect to z and using the condition that dp₀/dr is independent of z, we obtain

σ_{(rz) 0} (X) {= We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}),

(71)

where dp₀/dr is the pressure gradient in the physical system. Substituting Eq. into Eq. and solving for the shear rate to zero order in α, we get

\begin{matrix} X (r, z) = \frac{J_{r \infty}^{2} {({We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}))}^{2} - 1 + \sqrt{{(J_{r \infty}^{2} {({We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}))}^{2} - 1)}^{2} {+ 4 J}_{r \infty}^{2} φ_{r} {({We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}))}^{2}}}{{2 J}_{r \infty}^{2} {We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2})} . \end{matrix}

(72)

Equation contains to asymptotic cases at low and high Weissenberg numbers, respectively,

{We}_{\infty} \frac{\partial v_{r 0}}{\partial z} = {Lim}_{{We}_{\infty} \to 0} X (r, z) {= φ}_{r} {We}_{\infty} \frac{{dp}_{0}}{dr} z (z - 1) .

(72a)

The radial velocity at low Weissenberg number, using the boundary conditions given in Eq. , is

v_{r 0} (r, z) = φ_{r} \frac{1}{2} \frac{{dp}_{0}}{dr} z (z - 1) .

(72b)

At high Weissenberg number:

{We}_{\infty} \frac{\partial v_{r 0}}{\partial z} = {Lim}_{{We}_{\infty} \to \infty} X (r, z) {= We}_{\infty} \frac{{dp}_{0}}{dr} z (z - 1) .

(73a)

Integrating Eq. with respect to the radial coordinate and using the boundary condition, the following expression is obtained:

v_{r 0} (r, z) = \frac{1}{2} \frac{{dp}_{0}}{dr} z (z - 1) .

(73b)

The full analytical solution of Eq. can be calculated as

v_{r 0} (r, z) = \frac{1}{2} \frac{1}{{(J_{r \infty} {We}_{\infty})}^{2} \frac{{dp}_{0}}{dr}} I (z, \frac{{dp}_{0}}{dr}, {Jr}_{\infty} {, We}_{\infty}, φ_{r}),

(74a)

where the integral I, is given by the following analytical expression as a function of the dimensionless numbers and the pressure gradient:

\begin{matrix} I (z, \frac{{dp}_{0}}{dr}, φ_{r}, J_{r \infty} {, We}_{\infty}) = \frac{1}{2} (\begin{array}{l} u^{2} (z) {- u}_{0}^{2} + \sqrt{{(u (z) - 1)}^{2} {+ 4 φ}_{r} u (z)} - \sqrt{{(u_{0}^{2} - 1)}^{2} {+ 4 φ}_{r} u_{0}^{2}} \\ + 2 Ln | \frac{1 - u_{0}^{2} + \sqrt{{(u_{0}^{2} - 1)}^{2} {+ 4 φ}_{r} u_{0}^{2}}}{1 - u^{2} (z) + \sqrt{{(u^{2} (z) - 1)}^{2} {+ 4 φ}_{r} u^{2} (z)}} | \\ + 2 (φ_{r} - 1) Ln | \frac{1 - (u^{2} (z) {+ 2 φ}_{r}) \sqrt{{(u^{2} (z) - 1)}^{2} {+ 4 φ}_{r} u^{2} (z)}}{1 - (u_{0}^{2} {+ 2 φ}_{r}) \sqrt{{(u_{0}^{2} - 1)}^{2} {+ 4 φ}_{r} u_{0}^{2}}} | \end{array}) . \end{matrix}

(74b)

The functions u(z), u₀ are given by

u (z) = {(J_{r \infty} {We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}))}^{2}; u_{0} = u (z = 0) = {(J_{r \infty} {We}_{\infty} \frac{{dp}_{0}}{dr} (- \frac{1}{2}))}^{2} .

(74c)

Note that the radial velocity profile is a nonlinear function of the pressure gradient and depends on the product Jr_∞We_∞ and φ_r.

3. Volumetric flow rate

The radial volumetric flow rate exiting the continuous squeeze flow film can be expressed as

Q_{0} = \int_{0}^{1} v_{r 0} (r, z) {dz = v}_{r 0} (r, z) {z |}_{z = 0}^{z = 1} - \int_{0}^{1} \frac{\partial v_{r 0} (r, z)}{\partial z} z dz = - \frac{1}{{We}_{\infty}} \int_{0}^{1} X (r, z) zdz .

(75a)

Substituting the shear strain X, we obtain

Q_{0} = \frac{1}{{We}_{\infty}} \frac{1}{{We}_{\infty} (- \frac{{dp}_{0}}{dr})} \int_{- \frac{1}{2} {We}_{\infty} \frac{{dp}_{0}}{dr}}^{\frac{1}{2} {We}_{\infty} \frac{{dp}_{0}}{dr}} \frac{J_{r \infty}^{2} u^{2} - 1 + \sqrt{{(J_{r \infty}^{2} u^{2} - 1)}^{2} {+ 4 J}_{r \infty}^{2} φ_{r} u^{2}}}{{2 J}_{r \infty}^{2} u} (\frac{1}{2} + \frac{u}{{We}_{\infty} \frac{{dp}_{0}}{dr}}) du .

(75b)

In Eq. , the following variable changes have been used:

u (z, r) {= We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2}) .

(75c)

The radial volumetric flow rate can be expressed as

Q_{0} = \frac{1}{{We}_{\infty}} (\frac{1}{2} \frac{1}{{We}_{\infty} (- \frac{{dp}_{0}}{dr})} I_{10}),

(75d)

where the integral I₁₀ is given by

I_{10} = \frac{1}{2} (\begin{array}{l} J_{r \infty}^{2} u_{1}^{2} + \sqrt{{(J_{r \infty}^{2} u_{1}^{2} - 1)}^{2} + 4 φ_{r} J_{r \infty}^{2} u_{1}^{2}} - 1 \\ + 2 Ln | \frac{2}{1 - J_{r \infty}^{2} u_{1}^{2} + \sqrt{{(J_{r \infty}^{2} u_{1}^{2} - 1)}^{2} + 4 φ_{r} J_{r \infty}^{2} u_{1}^{2}}} | \\ + 2 (φ_{r} - 1) \\ Ln | \frac{1 - J_{r \infty}^{2} u_{1}^{2} - 2 φ_{r} + \sqrt{{(J_{r \infty}^{2} u_{1}^{2} - 1)}^{2} + 4 φ_{r} J_{r \infty}^{2} u_{1}^{2}}}{2 (1 - φ_{r})} | \end{array})

(75e)

and

u_{1} = u (r, z = 1) = \frac{1}{2} {We}_{\infty} \frac{{dp}_{0}}{dr} .

(75f)

4. Axial component of the field velocity

From the mass conservation [Eq. ], the axial velocity to zeroth-order in α is given by

v_{z 0} (z, r) = u (r) - \frac{1}{r} \frac{\partial}{\partial r} r (I_{w} (r, z)),

(76a)

where the Integral I_w(r,z) is given by

I_{w} (r, z) = \int_{0}^{z} v_{r 0} (r, β) d β = v_{r 0} (r, z) z - \frac{1}{{We}_{\infty}} \int_{0}^{z} X (r, β) dz = v_{r 0} (r, z) z + Q (z) .

(76b)

Q is the volumetric flow rate as a function of the orthogonal variable z, i.e., Q = Q(z). The mathematical function can be expressed as

Q (z) = \frac{1}{{We}_{\infty}^{2} (- \frac{{dp}_{0}}{dr})} \int_{- \frac{1}{2} {We}_{\infty} \frac{{dp}_{0}}{dr}}^{\frac{1}{2} {We}_{\infty} \frac{{dp}_{0}}{dr} (z - \frac{1}{2})} \frac{J_{r \infty}^{2} u^{2} - 1 + \sqrt{{(J_{r \infty}^{2} u^{2} - 1)}^{2} {+ 4 J}_{r \infty}^{2} φ_{r} u^{2}}}{{2 J}_{r \infty}^{2} u} (\frac{1}{2} + \frac{u}{{We}_{\infty} \frac{{dp}_{0}}{dr}}) du .

(76c)

Equation satisfies w₀ (z = 0, r) = u(r) and for w₀ (z =1, r) = 0 as follows:

\frac{1}{r} \frac{\partial}{\partial r} r (I_{w} (r, z)) = u (r) .

(76d)

Equation represents the axial velocity in our system to zeroth order in α. The asymptotic value of the axial velocity at low Weissenberg number is given by

v_{z 0} (z, r) = {Lim}_{{We}_{\infty} \to 0} v_{z 0} (z, r) = u (r) ({1 + 2 z}^{3} - {3 z}^{2}); \frac{1}{r} \frac{\partial}{\partial r} (r \frac{{dp}_{0}}{dr}) = - φ_{r}^{- 1} 12 u (r) .

(77a)

At high Weissenberg number:

v_{z 0} (z, r) = {Lim}_{{We}_{\infty} \to \infty} v_{z 0} (z, r) = u (r) ({1 + 2 z}^{3} - {3 z}^{2}); \frac{1}{r} \frac{\partial}{\partial r} (r \frac{{dp}_{0}}{dr}) = - 12 u (r) .

(77b)

5. Pressure gradient

To obtain the pressure gradient in the system as a function of the material properties and dimensionless numbers, Eq. is integrated with respect to the radial coordinate:

Q_{0} {= I}_{w} (r, z = 1) = \frac{1}{r} \int_{0}^{r} r' u (r') dr' .

(78)

The non-lineal Eq. is a function of the pressure gradient to zeroth order, i.e., Q₀ = Q₀ (dp₀/dr), dimensionless numbers (J_r∞, φ_r, and We_∞), and the non-homogeneous flow associated with the geometry (lower plate, m parameter). To solve it, the following steps are taken:

Choosing a numerical value for the dimensionless numbers {We_∞, J_r∞, φ_r} and m related to the non-homogeneity of the system.
Choosing the analytical form of the radial velocity distribution u(r), associated with the distribution of the holes in the lower plates. The continuous function u (r) must satisfy the normalized Eq. (21), so the mathematical u (r) function is given by

${u (r) = (m + 2) r}^{m} .$

The homogenous case appears when m = 0, i.e., u_H = u (r = 0) = 2. Physically, this means that the volume of the complex non-Newtonian fluid entering through the orifices of the lower plate is the same, and there is no distribution in the system.
Solve the non-analytical expression for the pressure gradient dp₀/dr. Given the steps A-B, and chosen and numerical value of r ε [0,1]⊆ R, the system can be solved with an iterative procedure such as the Newton–Raphson method.⁴⁷
The starting value, based on the Newtonian solution, is given by Eq. (41), convergence is found to be very rapid in this case.

$\frac{{dp}_{0}}{dr} = \frac{2}{{We}_{\infty}} X (r, z = 1) = - 12 \frac{1}{r} (\int_{0}^{r} α u (α) d α) = - 12 r^{m + 1} .$

The pressure gradient has the following analytical asymptotic values. At low Weissenberg number, the expression is

\frac{{dp}_{0}}{dr} = - φ_{r}^{- 1} \frac{1}{r} \int_{0}^{r} 12 u (r) rdr = - φ_{r}^{- 1} 12 r^{m + 1} .

(79a)

Equation can be separated to obtain the pressure in the system as follows:

p_{0} (r) = p_{0} (r = 1) {+ φ}_{r}^{- 1} \frac{12}{m + 2} (1 - r^{m + 2}) .

(79b)

For high Weissenberg numbers:

\frac{{dp}_{0}}{dr} = - \frac{1}{r} \int_{0}^{r} 12 u (r) rdr = - 12 r^{m + 1} .

(79c)

The linear differential Eq. can be separated to obtain the profile to pressure gradient as follows:

p_{0} (r) = p_{0} (r = 1) + \frac{12}{m + 2} (1 - r^{m + 2}) .

(79d)

It is important to note that, the asymptotic limits only depend on the dimensionless ratio φ_r. In Sec. V C 6, the normal force is obtained.

6. Normal force

The dimensional normal force on the upper disk by the complex fluid is given by the total zz component of the stress tensor evaluated at the upper plate as follows:

F = - \int_{0}^{2 π} \int_{0}^{a} ({T_{(zz)} |}_{z = h}) rdrd θ = - 2 π \int_{0}^{a} {T_{(zz)} |}_{z = h} rdr .

(80a)

To calculate the normal force induced by the liquid over the upper plate, we use Eq. as the basis. By substituting the dimensionless variables defined in Eq. into Eq. and integrating by parts, the following expression is obtained:

{F = - 2 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} {{{\frac{1}{2} T_{(zz)} |}_{z = 1} r^{2} |}_{0}^{1} - \frac{1}{2} \int_{0}^{1} \frac{\partial}{\partial r} {T_{(zz)} |}_{z = 1} r^{2} dr} = {π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} \int_{0}^{1} \frac{\partial}{\partial r} {T_{(zz)} |}_{z = 1} r^{2} dr .

(80b)

Applying the power expansion to Eq. and using the total stress tensor to zeroth order in α, we have We_∞T_(zz)0|_{z =1} = −We_∞p ₀. The non-dimensional force to zeroth order in α, denoted as (F₀), is

F_{0} = - {π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} {\int_{0}^{1} \frac{\partial p_{0}}{\partial r} |}_{z = 1} r^{2} dr = {- π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} {\int_{0}^{1} \frac{2}{{We}_{\infty}} σ_{(rz) 0} |}_{z = 1} r^{2} dr .

(80c)

Finally, Eq. can be expressed in terms of the shear strain evaluated at z =1, i.e., X₁(r) = X (r, z =1),

F_{0} (J_{r \infty} {, φ}_{r} {, We}_{\infty}) = - {π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} ({2 We}_{\infty}^{- 1} \int_{0}^{1} \frac{\frac{J_{r \infty}^{2}}{φ_{r}} X_{1} {(r)}^{2} - 1 + \sqrt{{(\frac{J_{r \infty}^{2}}{φ_{r}} X_{1} {(r)}^{2} - 1)}^{2} + 4 \frac{J_{r \infty}^{2}}{φ_{r}} \frac{1}{φ_{r}} X_{1} {(r)}^{2}}}{2 \frac{J_{r \infty}^{2}}{φ_{r}} X_{1} (r)} r^{2} dr) .

(80d)

Equation represents the dimensionless equation for the zeroth order contribution. The following mechanisms are studied in the continuous squeeze flow film: (i) shear thinning: φ_r = φ₀/φ_∞ < 1, (ii) shear thickening: φ_r = φ₀/φ_∞ > 1, (iii) Newtonian case: φ_r = 1 for any particular values of [J_r∞, We_∞], (iv) Pseudo-Thixotropy: φ_r = φ₀/φ_∞ < 1 and (We_∞)_min ≤ We_∞ ≤ (We_∞)_max, and (v) Pseudo-Rheopexy: φ_r = φ₀/φ_∞ > 1 and (We_∞)_min ≤ We_∞ ≤ (We_∞)_max, (vi) yield stress: φ_r → 0.

a. Lower shear strain

Equation has two asymptotic limits at low shear strain when z = 1, i.e., X₁(r):

{Lim}_{X_{1} (r) \to 0} F_{0} (J_{r \infty} {, φ}_{r} {, We}_{\infty}) {= - π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{{2 We}_{\infty}^{- 1}}{φ_{r}} \int_{0}^{1} X_{1} (r) r^{2} dr) .

(81a)

The shear strain given in Eq. is X₁ (r, z =1) = We_∞(−6r ^m+1). The non-homogeneous flow can be expressed as

{Lim}_{X_{1} (r) \to 0} F_{0} (J_{r \infty} {, φ}_{r} {, We}_{\infty}) {= π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} \frac{12}{φ_{r}} \int_{0}^{1} r^{m + 3} dr = \frac{3 {π a}^{3} 〈 V 〉}{φ_{\infty} φ_{r} h^{2}} (\frac{4}{m + 4}) .

(81b)

For homogeneous flow with m = 0, the following expression is:

F_{N 0} = F (m = 0) = \frac{3 {π a}^{3} 〈 V 〉}{φ_{\infty} φ_{r} h^{2}} .

(81c)

b. High shear strain

At high shear rates, the following expression is obtained as a function of the material properties:

{Lim}_{X_{1} (r) \to \infty} F_{0} (J_{r \infty} {, φ}_{r} {, We}_{\infty}) {= - π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} \frac{{2 We}_{\infty}^{- 1}}{φ_{r}} \int_{0}^{1} X_{1} (r) r^{2} dr .

(82a)

Similarly, the expression obtained is

{Lim}_{X_{1} (r) \to \infty} F_{0} (J_{r \infty} {, φ}_{r} {, We}_{\infty}) = \frac{3 {π a}^{3} 〈 V 〉}{φ_{\infty} h^{2}} (\frac{4}{m + 4}) .

(82b)

For homogeneous flow with m = 0, the expression is

F_{N \infty} {= F}_{0} (m = 0) = \frac{3 {π a}^{3} 〈 V 〉}{φ_{\infty} h^{2}} .

(82c)

7. Yield stress

The yield stress mechanism is relevant when φ_r→0. Thus, the integral of the force is given by

F_{0} (J_{r \infty, φ r,} {We}_{\infty}) ≅ {- π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} ({2 We}_{\infty}^{- 1} \int_{0}^{1} \frac{r^{2}}{J_{r \infty}} dr) = {π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (2 \frac{{We}_{\infty}^{- 1}}{J_{r \infty}} \int_{0}^{1} r^{2} dr) .

(83a)

The force corresponding to the yield stress is given by

F_{y - stress} {= - F}_{0} (J_{r \infty} {, φ}_{r} \to {0, We}_{\infty}) ≅ 3 {π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{2}{9} \frac{{We}_{\infty}^{- 1}}{J_{r \infty}}) .

(83b)

D. First-order theory O (α¹)

To first order in α, the equations of continuity and momentum are given by

\frac{\partial v_{z 1}}{\partial z} + \frac{1}{r} \frac{\partial}{\partial r} ({rv}_{r 1}) = 0,

(84)

0 = \frac{\partial σ_{(rz) 1}}{\partial z} {- We}_{\infty} \frac{\partial p_{1}}{\partial r} + \frac{1}{r} \frac{\partial}{\partial r} ({r σ}_{(rr) 0}) - \frac{σ_{(θ θ) 0}}{r},

(85)

{0 = - We}_{\infty} \frac{\partial p_{1}}{\partial z} + \frac{\partial σ_{(zz) 0}}{\partial z} .

(86)

The boundary conditions to first order in α are

v_{r 1} (r, z) = v_{z 1} (r, z) = 0 on z = 0,

(87a)

v_{r 1} (r, z) = v_{z 1} (r, z) = 0 on z = 1,

(87b)

v_{r 1} (r, z) = 0 on r = 0 .

(87c)

Equation can be simplified to

{We}_{\infty} \frac{\partial p_{1}}{\partial z} = \frac{\partial σ_{(zz) 0}}{\partial z} \Rightarrow \frac{\partial}{\partial z} [{We}_{\infty} p_{1} {- σ}_{(zz) 0}] = 0 \Rightarrow {We}_{\infty} p_{1} = Ξ {(r) + σ}_{(zz) 0} .

(87d)

From the zeroth order theory, the normal stress components σ_(zz)0 = σ₍_θθ₎₀ are zero, i.e., σ_(zz)0 = σ₍_θθ₎₀ = 0. Assuming cylindrical symmetry (∂p₁/∂θ = 0), and using Eq. , the pressure to first order, is a function exclusively of the radial coordinate r, i.e., We_∞p₁ = Ξ(r).

1. Radial component of the field velocity

Integrating the momentum equation to first order in α, the following expression is obtained:

σ_{(rz) 1} {= We}_{\infty} \frac{\partial p_{1}}{\partial r} z - \frac{{We}_{\infty}^{2}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) {(z - \frac{1}{2})}^{3} + C .

(88)

Assuming symmetry in the shear stress (σ_(rz)1 = 0 when z = 0), the following expression is obtained:

σ_{(rz) 1} {= We}_{\infty} \frac{\partial p_{1}}{\partial r} (z - \frac{1}{2}) - \frac{{We}_{\infty}^{2}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) {(z - \frac{1}{2})}^{3} .

(89)

Using a local approximation, for the shear stress at first order in the parameter α (Herrera et al., 2009, 2010; Herrera-Valencia et al., 2019)^83,86–88

\frac{\partial σ_{(rz) 0}}{\partial X} {We}_{\infty} \frac{\partial v_{r 1}}{\partial z} = We_{\infty} \frac{\partial p_{1}}{\partial r} (z - \frac{1}{2}) - \frac{{We}_{\infty}^{2}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) {(z - \frac{1}{2})}^{3} .

(90)

Integrating with respect to the axial coordinate z and using the boundary conditions and , the general solution is given by

v_{r 1} (r, z) = \frac{\partial p_{1}}{\partial r} \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} (z - \frac{1}{2}) dz - \frac{{We}_{\infty}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} {(z - \frac{1}{2})}^{3} dz .

(91a)

a. Low shear rate X

At small shear strain X, the following expression is obtained:

v_{r 1} (r, z) = φ_{r} \frac{1}{2} \frac{\partial p_{1}}{\partial r} ({(z - \frac{1}{2})}^{2} - {(\frac{1}{2})}^{2}) - φ_{r} \frac{{We}_{\infty}}{12} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) ({(z - \frac{1}{2})}^{4} - {(\frac{1}{2})}^{4}) .

(91b)

It is important to note that the effect of the viscosity is contained in the pressure gradient at zeroth order.

b. High shear rate

At high shear strain X, the following expression is obtained:

v_{r 1} (r, z) = \frac{1}{2} \frac{\partial p_{1}}{\partial r} ({(z - \frac{1}{2})}^{2} - {(\frac{1}{2})}^{2}) - \frac{{We}_{\infty}}{12} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) ({(z - \frac{1}{2})}^{4} - {(\frac{1}{2})}^{4}) .

(91c)

2. Axial component of the field velocity to first order

From the continuity equation, the general expression for the axial velocity is given as follows:

v_{z 1} (r, z) = v_{z 1} (r, z = 0) - \frac{1}{r} \frac{\partial}{\partial r} (r \int_{0}^{z} v_{r 1} (r, \bar{z}) d \bar{z}) = - \frac{1}{r} \frac{\partial}{\partial r} ({rQ}_{1} (r, z)),

(92a)

where the volumetric flow rate can be expressed in the following general form:

Q_{1} (r, z) = \int_{0}^{z} {vr}_{1} (r, \bar{z}) dz .

(92b)

From the second boundary condition of the axial velocity [v_z1(r, z =1) = 0], the following restriction is implied:

\frac{1}{r} \frac{\partial}{\partial r} ({rQ}_{1} (r, z)) = 0 .

(92c)

a. Low shear rate

At small shear strain, the volumetric flow rate, can be expressed in the following form:

\begin{matrix} Q_{1} (r, z) = φ_{r} \frac{1}{2} \frac{\partial p_{1}}{\partial r} {(\frac{1}{3} {(z - \frac{1}{2})}^{3} - z {(\frac{1}{2})}^{2}) - \frac{1}{3} {(- \frac{1}{2})}^{3}} \\ - φ_{r} \frac{{We}_{\infty}}{12} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) {(\frac{1}{5} {(z - \frac{1}{2})}^{5} - z {(\frac{1}{2})}^{4}) - \frac{1}{3} {(- \frac{1}{2})}^{5}} . \end{matrix}

(93)

b. High shear rate

At high shear strain, the volumetric flow rate can be expressed as

\begin{matrix} Q_{1} (r, z) = \frac{1}{2} \frac{\partial p_{1}}{\partial r} {(\frac{1}{3} {(z - \frac{1}{2})}^{3} - z {(\frac{1}{2})}^{2}) - \frac{1}{3} {(- \frac{1}{2})}^{3}} \\ - \frac{{We}_{\infty}}{12} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) {(\frac{1}{5} {(z - \frac{1}{2})}^{5} - z {(\frac{1}{2})}^{4}) - \frac{1}{3} {(- \frac{1}{2})}^{5}} . \end{matrix}

(94)

3. Pressure gradient

The pressure gradient can be obtained through the second boundary condition of the axial component of the velocity field:

Q_{1} (r, z) = \frac{C_{1}}{r} \Rightarrow \int_{0}^{z} v_{r 1} (r, \bar{z}) d \bar{z} = \frac{C_{1}}{r} .

(95)

Equation must be finite at any value of the radial coordinate r, in particular, if r = 0, the volumetric flow to first order goes to infinity, i.e., Q₁→∞. Therefore, the constant C₁, must be zero (C₁ = 0). The general expression for the pressure gradient to first order, can be expressed as

\frac{\partial p_{1}}{\partial r} = \frac{{We}_{\infty}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) \frac{\int_{0}^{1} \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} {(\bar{z} - \frac{1}{2})}^{3} d \bar{z} dz}{\int_{0}^{1} \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} (\bar{z} - \frac{1}{2}) d \bar{z} dz} .

(96a)

a. Low shear rate

At small shear rates X, the pressure gradient is obtained as follows:

\frac{\partial p_{1}}{\partial r} = \frac{{We}_{\infty}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) (\frac{3}{5}) {(\frac{1}{2})}^{2} = \frac{{We}_{\infty}}{3} φ_{r}^{- 2} 4^{2} \cdot 3^{2} (2 m + 3) r^{2 m + 1} (\frac{3}{5}) {(\frac{1}{2})}^{2} .

(96b)

b. High shear rate

At high shear rates X, the pressure gradient is

\frac{\partial p_{1}}{\partial r} = \frac{{We}_{\infty}}{3} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) (\frac{3}{5}) {(\frac{1}{2})}^{2} = \frac{{We}_{\infty}}{3} 4^{2} \cdot 3^{2} (2 m + 3) r^{2 m + 1} (\frac{3}{5}) {(\frac{1}{2})}^{2} .

(96c)

4. Normal force

To first order, the normal force can be calculated similarly to the zeroth order. Applying the power expansion to Eq. and using the total stress tensor to first order in α, We_∞T_(zz)1|_{z =1} = - We_∞p₁, the non-dimensionless force to first order in α (F₁) is

F_{1} {= 3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (- \frac{1}{3} {\int_{0}^{1} \frac{\partial p_{1}}{\partial r} |}_{z = 1} r^{2} dr) .

(97a)

Equation , is a function of shear strain evaluated in z =1, i.e., X₁(r) = X (r, z =1). Once the first-order pressure gradient [Eq. ] is used,

F_{1} = {3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (- \frac{{We}_{\infty}}{3} \int_{0}^{1} \frac{1}{r} \frac{\partial}{\partial r} (r {(\frac{\partial p_{0}}{\partial r})}^{2}) (\frac{\int_{0}^{1} \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} {(\bar{z} - \frac{1}{2})}^{3} d \bar{z} dz}{\int_{0}^{1} \int_{0}^{z} \frac{\partial X}{\partial σ_{(rz) 0}} (\bar{z} - \frac{1}{2}) d \bar{z} dz}) r^{2} dr) .

(97b)

Equation represents the contribution of elastic mechanisms in the continuous squeeze flow film.

a. Low shear rate

At small shear strain, the elastic contribution is given by

F_{1} ({m, We}_{\infty}) = {3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{1}{3} 4^{2} \cdot 3^{2}) (\frac{3}{5}) {(\frac{1}{2})}^{2} (\frac{2 m + 3}{2 m + 4}) ({- φ}_{r}^{- 2} {We}_{\infty}) .

(98a)

In the case of a homogeneous flow (m = 0), the following expression is obtained:

F_{1} = F_{1} ({0, We}_{\infty}) = {3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{1}{3} 4^{2} \cdot 3^{2}) (\frac{3}{5}) {(\frac{1}{2})}^{2} (\frac{3}{4}) φ_{r}^{- 2} (- {We}_{\infty}) .

(98b)

b. High shear rate

At high shear rates, the elastic contribution is

F_{1} (m) = {3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{1}{3} 4^{2} \cdot 3^{2}) (\frac{3}{5}) {(\frac{1}{2})}^{2} (\frac{2 m + 3}{2 m + 4}) ({- We}_{\infty}) .

(98c)

For homogeneous flow (m = 0), the expression is

F_{1} = {3 π a}^{3} \frac{1}{φ_{\infty}} \frac{〈 V 〉}{h^{2}} (\frac{1}{3} 4^{2} \cdot 3^{2} (\frac{3}{4})) (\frac{3}{5}) {(\frac{1}{2})}^{2} ({- We}_{\infty}) .

(98d)

Important points: (i) At small and high Weissenberg numbers, the force to first order depends linearly on the Weissenberg number We_∞ and (ii) the contribution of the elastic mechanisms is negative, meaning that the compression direction is opposite to that of the viscous mechanisms.

VI. SIMULATIONS AND RESULTS

In this section, the simulations of the ESR constitutive equation in a continuous squeeze flow film at zeroth and first orders are presented. The following mechanisms are explored: (i) shear thinning, (ii) shear thickening, (iii) thixotropy, (iv) rheopexy, (v) yields stress, and (vi) the weight of the worm-like micelles. The dimensionless numbers used in the simulations are given by: {We_∞, φ_r, J_r∞}. The analytical and numerical equations that are used in the simulations are (i) the axial velocity [Eq. (56)], (ii) the normal force to zeroth order in α (67), and (iii) the normal force to first order in α (67). The axial velocity is determined analytically, while the normal forces are computed numerically. The non-linear pressure gradient is calculated using the Newton–Raphson method, and the numerical integral is evaluated using a quadratic Gaussian technique. The results were programed in Wolfram Mathematica using institutional license codes. The numerical values of the material properties are taken from Herrera-Valencia et al. (2017).⁸³ The dimensionless numbers as a function of the weight concentration of the worm-like micelles solution are given in Table I. In what follows, the first analysis will be the zeroth-order in alfa, and the second analysis will be the first-order contribution. To facilitate the calculation, the normal force will be normalized by the Newtonian value, i.e., F₀/F_N∞ and F₁/F_N∞, where F_N∞ =3πa³⟨V⟩(φ_∞)⁻¹/h².

TABLE I.

Experimental material parameters of the ESR and its values.

Material properties vs mechanisms	Irreversible thermodynamical viscoelastic work constant β₀ (Pa⁻¹ s)	Interfusion structure-momentum relaxation time λ_D (s)	Elastic moduli J₀ =1/G₀ (a⁻¹)	Fluidity at low shear rate φ₀ =1/η₀ (Pa⁻¹ s⁻¹)	Fluidity at high shear rate φ_∞ =1/η_∞ (Pa⁻¹ s⁻¹)
Shear-thinning	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.20
Shear-thinning	5.46 × 10⁻⁶	0.14	0.0054	0.0053	1
Newtonian	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.0053
Shear-thickening	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.0002
Shear-thickening	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.002
Yield stress	5.46 × 10⁻⁶	0.14	0.0054	0	10.5
				0.0001
				0.001
				0.0053
				1
				10.5
Thixotropy	5.46 × 10⁻⁶	0.14	0.0054	0.0053	1
	3.9 × 10⁻⁶	0.10
	3.9 × 10⁻⁷	0.01
	3.9 × 10⁻⁸	0.001
Concentration CTAT (wt. %)
5%	30	0.12	0.0240	0.0275	19.8
10%	10	0.33	0.0057	0.0061	15.0
15%	4	0.38	0.0072	0.0050	12.6
20%	1.8	0.42	0.0016	0.0042	12.0

Material properties vs mechanisms	Irreversible thermodynamical viscoelastic work constant β₀ (Pa⁻¹ s)	Interfusion structure-momentum relaxation time λ_D (s)	Elastic moduli J₀ =1/G₀ (a⁻¹)	Fluidity at low shear rate φ₀ =1/η₀ (Pa⁻¹ s⁻¹)	Fluidity at high shear rate φ_∞ =1/η_∞ (Pa⁻¹ s⁻¹)
Shear-thinning	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.20
Shear-thinning	5.46 × 10⁻⁶	0.14	0.0054	0.0053	1
Newtonian	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.0053
Shear-thickening	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.0002
Shear-thickening	5.46 × 10⁻⁶	0.14	0.0054	0.0053	0.002
Yield stress	5.46 × 10⁻⁶	0.14	0.0054	0	10.5
				0.0001
				0.001
				0.0053
				1
				10.5
Thixotropy	5.46 × 10⁻⁶	0.14	0.0054	0.0053	1
	3.9 × 10⁻⁶	0.10
	3.9 × 10⁻⁷	0.01
	3.9 × 10⁻⁸	0.001
Concentration CTAT (wt. %)
5%	30	0.12	0.0240	0.0275	19.8
10%	10	0.33	0.0057	0.0061	15.0
15%	4	0.38	0.0072	0.0050	12.6
20%	1.8	0.42	0.0016	0.0042	12.0

A. Zeroth order mechanisms: O (α⁰)

1. Shear-thinning/thickening

Figure 5 shows the normalized force to zeroth order F_0/F_N∞ at the upper plate in the continuous squeeze flow to zero order in α vs the Weissenberg number We_∞, as a function of the dimensionless number φ_r for values: = (i) 0.001, (ii) 0.01, (1), (5), (10). The other dimensionless number used in the simulation is J_r∞ =1. It is interesting to note that the force is proportional to the viscosity of the system, so the shape of the curves observed in Fig. 5 is the same as the viscosity of the system (shear thinning, Newtonian, and shear thickening) as follows:

FIG. 5.

Normalized force F0/FN∞ to zeroth order in the parameter α vs Weissenberg number (We∞) as a function of shear-thinning/thickening mechanisms, through the dimensionless numbers Jr∞ and φr. The letters a-e shows the different values of the dimensionless numbers. Here, the material properties used in the simulation are provided in Table I.

View large Download slide

Normalized force F₀/F_N∞ to zeroth order in the parameter α vs Weissenberg number (We_∞) as a function of shear-thinning/thickening mechanisms, through the dimensionless numbers Jr_∞ and φ_r. The letters a-e shows the different values of the dimensionless numbers. Here, the material properties used in the simulation are provided in Table I.

At low Weissenberg number We_∞, all the curves show a constant plateau independent of the Weissenberg number (Newtonian zone), for φ_r equal to 1, the force remains constant independent of the Weissenberg number We_∞ (Newtonian case).
At intermediate Weissenberg number We_∞, a power law region is shown for for φ_r > 1 (shear thinning) and φ_r < 1 (shear thickening). In this power law region, the systems experience constant structural changes, with flow and rheology dominated by the dimensionless numbers {J_r∞,φ_r}, which are associated with the yield stress mechanism, bulk-elastic shear forces, and the fluidities at low and high infinity Weissenberg numbers. From a mathematical perspective, in the shear thinning mechanism, the system displays a monotonically decreasing behavior, while in the shear thickening forces, a monotonically increasing behavior is observed.
At large Weissenberg number We_∞ ≫ 1, the normal force (F₀/F_N∞) becomes constant and independent of the Weissenberg number (second Newtonian plateau). All the simulations (a-e) illustrate different material conditions, demonstrating rheological transitions from a constant structure to a more compact structure, influenced by the flow and rheology mechanisms through the dimensionless numbers. The partial conclusions from this plot are summarized as follows:
- The normal force (F₀/F_N∞) is proportional to the viscosity function in the case of viscometric flow, specifically in the case of homogeneous flow (m = 0).
- The shear-thinning/thickening mechanisms are characterized by the dimensionless groups {J_r∞, φ_r}.

2. Pseudo thixotropy/rheopexy

Figure 6 shows the normalized force F₀/F_N∞ to zeroth order vs Weissenberg number We_∞, as a function of the viscoelastic dissipation function at high shear rates and the bulk elastic mechanisms through the dimensionless number J_r∞. The other parameter fixed in the simulation is φ_r = 1.

FIG. 6.

Normalized force to zero order in α (F0/FN∞) vs Weissenberg number (We∞) of the number Jr∞ (pseudo thixotropy/rheopexy mechanisms). The parameters used in the simulation are given in Table I. In simulations (a) and (b), the dimensionless number that control the pseudo thixotropy and rheopexy mechanisms is the Jr∞.

View large Download slide

Normalized force to zero order in α (F₀/F_N∞) vs Weissenberg number (We_∞) of the number Jr_∞ (pseudo thixotropy/rheopexy mechanisms). The parameters used in the simulation are given in Table I. In simulations (a) and (b), the dimensionless number that control the pseudo thixotropy and rheopexy mechanisms is the J_r∞.

It is evident the competition of all these mechanisms causes the complex fluid not to recover after a period of deformation F₀/F_N∞. At zeroth order, different deformation histories are shown by the loops and they correspond to different relaxation times, as these curves are obtained in steady state, we have decided to call them pseudo thixotropy/rheopexy as the shape of the curves is similar to the unsteady case. Mathematically, Fig. 6 display two plateaus at small and high Weissenberg numbers and intermediate power low region with a slope close to unity. Unsteady state simulations may require more powerful numerical methods and they will render non analytical solutions which are outside the objectives of the present work and will be studied as future work.

The following important points are summarized as follows:

The normal force F₀/F_N as a function of the Weissenberg number We_∞ [i.e., F₀/F_N = F₀/F_N(We_∞)] is proportional to the apparent shear viscosity in a parallel plate system sheared by a pressure gradient dp₀/dr.
The normalized force and apparent shear stress loops are controlled by the two dimensionless numbers [J_r∞, φ_r], respectively.
At small and large Weissenberg numbers, both systems are independent of J_r∞, and its values is determined by the fluidity ratio φ_r.
The normalized plateau can be extended when the compliance associated to the viscoelastic irreversible work is small compared to bulk-elastic compliance. As J_r∞ → 0, the plateau extends, and the structured fluid require more energy to change the microscopic entanglements, bonds or links. In contrast, when J_r∞→ ∞, the plateau changes with less energy.

3. Yield stress

Figure 7 shows the normal force and the apparent shear stress vs Weissenberg number We_∞ as a function of the fluidity ratio φ_r = φ₀/φ_∞. The simulation (a) corresponds to the case when the system is Newtonian, and the fluidity ratio is equal to the unity, i.e., φ_r = 1. Here, the system is independent of the value of J_r∞ and the structure is completely disrupted, and it is equivalent to the Newtonian solvent viscosity. The simulations (b)–(d) show a constant behavior at small and moderate Weissenberg number We_∞ = J₀⟨V⟩/hφ_∞. At a critical value of the We_∞, the activation of the viscoelastic dissipation structure reaches its maximum, and the compact structure shows a monotonically decreasing behavior, which can be considered a small shear-thinning region followed by plateau behavior. Finally, the last corresponds to the highest value of the normal force at small and moderate We_∞ numbers (φ_r → 0). The following summary finding are

FIG. 7.

Normalized force to zeroth order in α (F0/FN∞) vs Weissenberg number (We∞) as a function of the fluidity ratio φr (yield stress mechanisms). Here, the value of the infinity compliance ratio is Jr∞ = 0.432 283. Panels (a)–(e) indicate different dimensionless numbers value as functions of the material physical properties are given in Table I.

View large Download slide

Normalized force to zeroth order in α (F₀/F_N∞) vs Weissenberg number (We_∞) as a function of the fluidity ratio φ_r (yield stress mechanisms). Here, the value of the infinity compliance ratio is Jr_∞ = 0.432 283. Panels (a)–(e) indicate different dimensionless numbers value as functions of the material physical properties are given in Table I.

The dimensionless value of φr controls the first plateau when the structure is constant and independent of the Weissenberg number (We_∞).
When φ_r → 0, the first plateau decreases drastically, and the system reaches the maximum compact structure, which its value determined by the force at zeroth Weissenberg value.
The value of the normal force in the maximum fluid structure is related with a simplified equation, given by the following expression:

$\frac{F_{0}}{F_{N \infty}} = \frac{2}{3^{2}} σ_{y} {({We}_{\infty})}^{- 1} = \frac{2}{3^{2}} {(J_{r \infty})}^{- 1 / 2} \cdot {({We}_{\infty})}^{- 1} .$

This equation can be linearized to theoretically obtain the yield stress value.

4. Wormlike micellar solutions data

Figure 8 shows the effect of CTAT content on the normal force and apparent shear stress as a function of the Weissenberg number. The material properties used in the simulations are provided in Table I. In both simulations, the effect of the concentration induces greater structural-interactions. Generally, Fig. 8 exhibits the same mathematical behavior as shown in the shear-thinning mechanisms (Fig. 5). At small and high Weissenberg numbers, the system displays two plateaus, while at moderate Weissenberg numbers, there is a monotonically decreasing and increasing behavior for the normal force and apparent shear stress. The main effects of the concentration on the load capacity of the continuous squeeze flow film are summarized as follow:

FIG. 8.

Normalized force at zeroth order in α (F0/FN∞) vs Weissenberg number (We∞) as a function of the CTAT content in the system. Panels (a)–(d) represent different micellar CTAT weight concentration. The material properties used in the simulation are provided in Table I.

View large Download slide

Normalized force at zeroth order in α (F₀/F_N∞) vs Weissenberg number (We_∞) as a function of the CTAT content in the system. Panels (a)–(d) represent different micellar CTAT weight concentration. The material properties used in the simulation are provided in Table I.

Increasing the CTAT content decreases the load-bearing capacity of the continuous squeeze flow film.
At high concentration (20 wt. %), the force is at a minimum compared to the other sample concentrations. Physically, as the concentration increases, the system exhibits more entanglement, links, or bonds, resulting in a more compact liquid structure.

B. Velocity profiles

In this section, the mechanisms of shear thinning/thickening, yield stress, and sample weight concentration are explored through the axial velocity profile. The material properties are obtained directly from Table I. In all cases, the system exhibits a parabolic profile, where the velocity is zero at the wall (non-slip boundary condition), and the maximum value of the velocity v_r0max = v_r0 (r, z =1/2) is reached and can be calculated using Eqs. (74a) and (74b). The following important points are highlighted in Fig. 9:

FIG. 9.

View large Download slide

(a) Radial velocity profile at zeroth order in α vs Weissenberg number (We_∞) as a function of the dimensionless number φ_r (shear thinning/thickening mechanisms). The parameters used in the simulation are provided in Table I. (b) Radial velocity profile at zero order in α vs Weissenberg number (We_∞) as a function of the dimensionless number φ_r (yield stress mechanisms). The parameters used in the simulation are provided in Table I. (c) Radial velocity profile at zeroth order in α vs Weissenberg number (We_∞) as a function of the weight sample concentration (worm-like micellar data). The parameters used in the simulation are provided in Table I.

Newtonian profile (red color) is observed when φ_r = 1, and is independent of Jr_∞, which represents the ratio between the yield stress compliance and bulk elastic mechanisms.
The shear thinning mechanisms display a larger parabolic profile compared to the Newtonian case (constant structure, simulations a and b). Here, the system is controlled by the dimensionless number φ_r ≪ 1.
The shear thickening mechanism shows a significant decrease in the parabolic profile (simulations d and e). Physically, the system exhibits a more compact structure due to internal mechanisms. Here, the dimensionless number φ_r is greater than one, i.e., φ_r ≫1.
Yield stress mechanisms are explained in Fig. 9(b). It is evident that the system shows a substantial decrease in the maximum velocity value. This behavior is observed when the φ_r → 0, i.e., φ₀ → 0 (infinite viscosity). In this case, the structured fluid has a highly compact structure and requires significant structural activation to flow approximating to plug flow (flat plateau, letters d and e).
Micellar concentration samples [Fig. 9(c)] show a parabolic velocity profile. As the concentration increases, the number or entanglements or links increases, and the velocity profile decreases considerably. The system is dominated by the two dimensionless numbers {J_r∞, φ_r}.

C. Non-homogeneous flow

Figure 10 shows the effect of the homogeneity through the parameter m on the normalized force at zeroth order. It is evident that the system exhibits the same rheological behavior as the viscosity shear-thinning mechanism (first Newtonian plateau, then power law zone, and finally second Newtonian plateau). The effect of the non-homogeneity, related to the parameter m, shifts the curve to lower values of the normalized force; however, the mathematical structure remains the same. Additionally, the non-homogeneity, represented by the parameter m, depends on the distribution of the holes in the lower plate and is completely determined by the experimental conditions.^73,74

FIG. 10.

Normalized force to zeroth order in α (F0/FN∞) vs Weissenberg number (We∞) as a function of the non-homogeneity parameter m (non-homogeneous flow). Pannel (a)–(e) shows the value of m employed in the simulation. The value of the dimensionless numbers Jr∞ = 0.193 323 and Jr∞ = 0.0265 (Letter b of Fig. 5). The parameters used in the simulation are provided in Table I.

View large Download slide

Normalized force to zeroth order in α (F₀/F_N∞) vs Weissenberg number (We_∞) as a function of the non-homogeneity parameter m (non-homogeneous flow). Pannel (a)–(e) shows the value of m employed in the simulation. The value of the dimensionless numbers Jr_∞ = 0.193 323 and Jr_∞ = 0.0265 (Letter b of Fig. 5). The parameters used in the simulation are provided in Table I.

D. First order mechanisms: O (α¹)

1. Elastic-shear-thinning/thickening

Figure 11(a) shows the force at first order in the parameter α vs the Weissenberg number (We_∞) as a function of the shear-thinning/thickening mechanisms through the dimensionless number φ_r. The first normalized force is related to the elastic contributions. For both mechanisms (shear thinning or thickening forces), the system exhibits three regions: (i) At low Weissenberg number, the system shows linear behavior (Power law). (ii) At intermediate Weissenberg numbers, the system shows a smooth, monotonically increasing behavior for the shear-thinning forces, while the shear-thickening mechanisms exhibit a pronounced transition. (iii) Finally, at high Weissenberg number, the system exhibits linear behavior again (Power law). Note that, in power law regions, the system undergoes constant changes in the elastic structure, and the flow and rheology are dominated by the numbers {J_r∞, φ_r}.

FIG. 11.

(a) Normalized force to first order in α (F1/FN∞) vs Weissenberg number (We∞) as a function of the non-dimensional number φr. The parameters used in the simulations are given in Table I. (b) Normalized force to first order in α (F1/FN∞) vs Weissenberg number (We∞) as a function of the number Jr∞. The parameters used in the simulations are given in Table I. (c) Normalized force to first order in α (F1/FN∞) vs Weissenberg number (We∞) as a function of φr. The parameters used in the simulations are given in Table I. (d) Normalized force to first order in α (F1/FN∞) vs Weissenberg number (We∞) as a function of the CTAT content in the system through the dimensionless numbers (Jr∞,φr). The parameters used are given in Table I.

View large Download slide

(a) Normalized force to first order in α (F₁/F_N∞) vs Weissenberg number (We_∞) as a function of the non-dimensional number φ_r. The parameters used in the simulations are given in Table I. (b) Normalized force to first order in α (F₁/F_N∞) vs Weissenberg number (We_∞) as a function of the number J_r∞. The parameters used in the simulations are given in Table I. (c) Normalized force to first order in α (F₁/F_N∞) vs Weissenberg number (We_∞) as a function of φ_r. The parameters used in the simulations are given in Table I. (d) Normalized force to first order in α (F₁/F_N∞) vs Weissenberg number (We_∞) as a function of the CTAT content in the system through the dimensionless numbers (Jr_∞,φ_r). The parameters used are given in Table I.

2. Elastic-pseudo thixotropy/rheopexy

Figure 11(b) shows the effect of the number Jr_∞ on the force at first order in the parameter α. It is clear that this number is related to the structural activation energy associated with the breakdown and buildup mechanisms through the shear-thinning/thickening mechanisms. At low Weissenberg numbers (We_∞ ≪ 1), the system shows linear behavior with a slope close the one. When 1 ≤ We_∞ < 10², it exhibits a smooth, monotonic transition, followed by linear behavior when 10² ≤ We_∞ < 10⁴. The elastic shear-thickening/thickening behavior shows a power law zone with an abrupt change in the internal microstructure, and for a critical Weissenberg number, it follows a linear behavior. These simulations show that elastic-shear-thinning/thickening mechanisms are coupled trough the pressure gradient to first order.

3. Elastic-yield/stress

Figure 11(c) shows the normal force to first order (F₁/F_N∞) vs Weissenberg number (We_∞) as a function of the material conditions associated with the yield stress forces. It is clear that, in the elastic-yield stress zone, the force at first order is independent of the Weissenberg number. At a critical value of the Weissenberg number, the system exhibits a linear behavior associated with the fluid zone (simulation a). The simulations (b-d) exhibit a mathematical behavior similar to shear-thinning mechanisms, with all these numerical results being controlled with the dimensionless number φ_r. The final case resembles Newtonian behavior.

4. Elastic-concentration wormlike micellar data

In Fig. 11(d) shows the effect of the force at first order vs the Weissenberg number as a function of the worm-like sample concentration in the system (see Table II). It is clear that the force at first order exhibits three basic zones. The first one, 0 < We_∞ < 10°, is associated with linear behavior. The second zone, 10° < We_∞ < 10², shows a monotonically increasing behavior, while the last zone shows a second linear behavior at high Weissenberg numbers. The simulations display the same mathematical behavior; however, as the worm-like micellar concentration increases, the elastic transition zone at moderate Weissenberg numbers decreases due to the increased structural interactions and elasticity.

TABLE II.

Dimensionless numbers vs sample concentration. Experimental results show that the gap between the plates h and the radius a are approximately of the order of 0.20 and 11 mm, respectively. The liquid density in this work is ρ =1000 kg/cm³ and the volumetric flow is of the order of Q = 5 × 10⁻⁶ m³/s. The average radial velocity ⟨V⟩ is given by ⟨Q⟩ = Q/2πah, so ⟨V⟩ ≅ 0.7234.

CTAT (wt. %)	Ca = J₀/J_m	We_∞ = J₀/Jv_∞	φr = φ₀/φ_∞	Jr = Jp/J₀	α = h/a
5%	0.076	720	0.000 105	3. 97 × 10⁻⁶	0.018
10%	0.015	2459	6.1 × 10⁻⁶	8.83 × 10⁻⁷
15%	0.0055	2520	2.18 × 10⁻⁶	7.23 × 10⁻⁷
20%	0.0024	2857	8.04 × 10⁻⁷	6.07 × 10⁻⁷

CTAT (wt. %)	Ca = J₀/J_m	We_∞ = J₀/Jv_∞	φr = φ₀/φ_∞	Jr = Jp/J₀	α = h/a
5%	0.076	720	0.000 105	3. 97 × 10⁻⁶	0.018
10%	0.015	2459	6.1 × 10⁻⁶	8.83 × 10⁻⁷
15%	0.0055	2520	2.18 × 10⁻⁶	7.23 × 10⁻⁷
20%	0.0024	2857	8.04 × 10⁻⁷	6.07 × 10⁻⁷

VII. SUMMARY AND CONCLUSIONS

The merits of this research focus on the study of structure–fluid interactions in a continuous squeeze-film flow of a complex liquid. The flow was simulated in a conventional squeeze film fixture by continuously injecting fluid into the narrow gap between two circular plates through the lower plate. The viscoelastic fluid is characterized by a new model based in the Phan-Thien–Tanner formalism.⁶⁴ The ESR non-linear equation of state describes the structure interactions through the activation energy associated with the minimal irreversible work to alter the structure under the influence of flow. This model incorporates five material properties related to shear-thinning/thickening, viscoelasticity, pseudo thixotropy/rheopexy, yield-stress, and banding mechanisms. The continuity, momentum, and rheology equations were scaled using characteristic variables, and the microscopic mechanisms are described by five dimensionless groups: (i) infinity Weissenberg number We_∞ = J₀/φ_∞h ⟨V⟩⁻¹ (viscoelastic mechanisms), (ii) fluidity ratio φ_r = φ₀/φ_∞ (shear thinning and thickening mechanisms), (iii) compliance ratio J_r∞ = Jp/J₀ (yield stress and bulk-elastic mechanisms), (iv) structure-interdiffusion Weissenberg We_λ_D = λ_D⟨V⟩/h (structure-interdiffusion and flow mechanisms), and finally (v) banding-flow Weissenberg We_λ_D = λ_B ⟨V⟩/h (flow instabilities). The scaled non-linear partial differential equations were solved with a regular perturbation-numerical technique to zeroth (viscosity forces) and first (elastic forces) orders, respectively. The study leads to the following conclusions.

A. Force to zeroth order: O (α⁰)

The force to zeroth order measures the effects of structure and is governed by two physical groups associated with the shear-thinning/thickening and yield/elastic-bulk/concentration mechanisms through the dimensionless groups (J_r∞, φ_r). The structure mechanisms can be evaluated through the numerical integral as follows:

$\frac{η}{η_{\infty}} \approx \frac{F_{0}}{F_{N \infty}} = \int_{0}^{1} (- \frac{{dp}_{0}}{dr}) r^{2} dr .$

Its magnitude is proportional to the variation in viscosity (inverse of the fluidity) with the rate of shear.
Shear thinning/thickening forces: The effects of the level of structure vs Weissenberg number can be measure through the dimensionless number φ_r = φ₀/φ_∞. The shear-thinning mechanisms are given by the following condition: φ_r = φ₀/φ_∞ < 1, and shear thickening mechanism: φ_r = φ₀/φ_∞ > 1. The Newtonian behavior is reached when φ_r = φ₀/φ_∞ = 1 (constant structure).
Pseudo thixotropy/rheopexy forces: The effect of the pseudo thixotropy and rheopexy (steady state, fully developed flow) can be evaluated using the dimensionless numbers (J_r∞, φ_r). This means that the competition between the yield-stress and bulk-elastic forces is related to the viscoelastic-irreversible work required to change the structure. When J_r∞→0, the plateau is extended, and the microstructure requires more structural activation energy to change. For very large J_r∞ values, i.e., J_r∞→∞, the first plateau decreases as a function of viscoelastic properties because the viscoelastic structural activation energy decreases, and the structure breaks down at lower Weissenberg numbers.
Yield-stress forces: Yield stress is maximized when the highest number of entanglements, links, or bonds. This effect is obtained when φ_r → 0, i.e, φ₀ → 0 (η₀ →∞). The yield stress force depends on the dimensionless numbers (Jr_∞, We_∞), and the scaling law is given by F_y-stress ≈ J_r∞⁻¹ We_∞⁻¹.
Weight concentration samples: The effect of different weight sample concentrations in the system is directly related to the numbers (J_r∞, We_∞). At low sample concentrations, the system promotes shear-thinning and pseudo thixotropy mechanisms. High micellar concentration induces more interactions through links, bonds or entanglements. The entire system is fully coupled, promoting the shear thickening/pseudo rheopexy/yield stress forces.

B. Force to first order: O (α¹)

Normal stress difference: The force to first order in the perturbation parameter α (F₁) measures the first normal stress difference, i.e., N₁ = σ_(rr)0-σ_(zz)0, and depends on the dimensionless numbers (J_r∞, φ_r, We_∞). It can be calculated through the following numerical integral F₁/F_N∞, where g (r) is a function of the radial coordinate:

$\frac{F_{1}}{F_{N \infty}} = - \int_{0}^{1} \frac{{dp}_{1}}{dr} r^{2} dr = - \int_{0}^{1} \frac{1}{r} \frac{d}{dr} (r {(\frac{{dp}_{0}}{dr})}^{2}) g (r) r^{2} dr \approx N_{1} = σ_{(rr) 0} - σ_{(zz) 0} .$
The valor of F₁/F_N∞ is negative, i.e., F₁/F_N∞ < 0. Physically, this is a consequence of the normal forces induced by the radial flow.
Small Weissenberg number We_∞: At small Weissenberg numbers We_∞, the system admits analytical solutions. The effects of elastic mechanisms decrease until a maximum value, determinate by a coupling effect between the dimensional groups. This demonstrates that although the overall effect is small and negative for small Weissenberg numbers, it is possible for normal stress effects to occur.
Load enhancement: The effects of the viscoelastic-thixotropic-elasto-plasto mechanisms demonstrate the marked load enhancement associate to normal-stress effects N₁ = σ_(rr)0-σ_(zz)0
Material properties: The material property (β₀, J₀, φ_r, λ_D, λ_D) associated with structure-activation energy are numerically fitted parameters of the models discussed. They implicitly depend on the underlying microscopic mechanism of micellar scissions and fusion.
Lubrication approximation: This system is equivalent to the traditional squeeze flow film, where the lubrication approximation is used, and dh(t)/dt is replaced by a dh(t)/dt = a <V>/2h, with a being the constant approach velocity of the disk in the conventional squeeze film.

It would be worthwhile to compare the theoretical predictions of the effect of pseudo thixotropy with experimental observations, for example, by using viscoelastic surfactants such as CTAT, EHAC, liquid crystalline suspensions or associative polymers. One of the most interesting effects of complex fluids is shear-banding flow, where mechanical and thermodynamic instabilities cause the system to separate into regions of different viscosities. The drastic shear-thinning behavior may lead to a significant enhancement in the context of lubrication theory.

ACKNOWLEDGMENTS

E.E.H.V. acknowledges the financial support of PAPIIT-DGAPA/UNAM project No. IN102823, PAPIME-DGAPA/UNAM project No. PE106224, and the theoretical discussion of James McGill Professor Alejandro Rey, chemical engineering department McGill University, Montreal, Quebec, Canada. L.A.R.T. gratefully appreciated the financial support from the Consejo Nacional de Ciencia y Tecnología (CONAHCyT) through CVU 860719. C.S.C. gratefully appreciated the financial support PAPIIT-DGAPA/UNAM project No. IN210023 and the Dirección General de Cómputo y de Tecnologías de Información y Comunicación (DGTIC) of the UNAM for allocation of computer time on the Miztli Supercomputer. V.J.H.A. acknowledges the financial support of PAPIIT-DGAPA/UNAM project No. IT-200323. This research is dedicated to the memory of my beloved father Emilio Herrera Caballero “el ave de las tempestades.”

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

L. A. Ramirez-Torres: Funding acquisition (equal); Methodology (equal); Software (lead); Visualization (equal). E. E. Herrera-Valencia: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). C. Soriano-Correa: Conceptualization (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal). M. L. Sánchez-Villavicencio: Data curation (equal); Project administration (equal); Software (equal); Validation (equal). L. Campos-Fernández: Formal analysis (equal); Software (equal); Visualization (equal). G. Ascanio: Formal analysis (equal); Software (equal); Supervision (equal); Validation (equal). V. J. Hernández-Abad: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). F. Calderas: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

REFERENCES

1.

A. A.

Dey

,

Y.

Modarres-Sadeghi

, and

J. P.

Rothstein

, “

Observation of lock-in for viscoelastic fluid-structure interactions

,”

J. Fluids Struct.

96

,

103025

(

2020

).

https://doi.org/10.1016/j.jfluidstructs.2020.103025

Google Scholar

Crossref

2.

K.

Zidi

,

B. D.

Texier

,

G.

Gauthier

, and

A.

Seguin

, “

Viscosimetric squeeze flow of suspensions

,”

Eur. Phys. J. E

47

,

17

(

2024

).

https://doi.org/10.1140/epje/s10189-024-00410-1

Google Scholar

Crossref

3.

Z.

Zheng

,

X.

Chen

, and

W.

Yang

, “

Squeeze flow of a Maxwell fluid between a sphere and a plate

,”

Phys. Fluids

36

,

013121

(

2024

).

https://doi.org/10.1063/5.0185335

Google Scholar

Crossref

4.

Z.

Zheng

,

H.

Xie

,

X.

Chen

,

X.

Liu

,

W.

Yang

,

Y.

Xu

, and

W.

Huang

, “

Squeeze flow of a Maxwell fluid between two parallel disks or two spheres

,”

Phys. Fluids

35

,

083105

(

2023

).

https://doi.org/10.1063/5.0161828

Google Scholar

Crossref

5.

R.

Dandekar

and

A. M.

Ardekani

, “

Nearly touching spheres in a viscoelastic fluid

,”

Phys. Fluids

33

,

083112

(

2021

).

https://doi.org/10.1063/5.0060222

Google Scholar

Crossref

6.

T. Y.

Zhou

,

D.

Lin

, and

Y. J.

Shen

, “

The squeeze flow of a bi-viscosity fluid between two rigid spheres with wall slip

,”

Particuology

79

,

153

–

160

(

2023

).

https://doi.org/10.1016/j.partic.2022.11.003

Google Scholar

Crossref

7.

J.

Engmann

,

C.

Servais

, and

A. S.

Burbidge

, “

Squeeze flow theory and applications to rheometry: A review

,”

J. Non-Newtonian Fluid Mech.

132

,

1

–

27

(

2005

).

https://doi.org/10.1016/j.jnnfm.2005.08.007

Google Scholar

Crossref

8.

R. J.

Poole

, “

Inelastic and flow-type parameter models for non-Newtonian fluids

,”

J. Non-Newtonian Fluid Mech.

320

,

105106

(

2023

).

https://doi.org/10.1016/j.jnnfm.2023.105106

Google Scholar

Crossref

9.

N.

Phan-Thien

,

S.

Nasseri

, and

L. E.

Bilston

, “

Oscillatory squeezing flow of a biological material

,”

Rheol. Acta

39

,

409

–

417

(

2000

).

https://doi.org/10.1007/s003970000093

Google Scholar

Crossref

10.

Y.

Aboelkassem

, “

Computational and theoretical model of electro-osmotic flow pumping in a microchannel with squeezing walls

,”

Phys. Fluids

35

,

052011

(

2023

).

https://doi.org/10.1063/5.0151554

Google Scholar

Crossref

11.

D.

K. Gartling

and

N.

Phan-Thien

, “

A numerical simulation of a plastic fluid in a parallel-plate plastometer

,”

J. Non-Newtonian Fluid Mech.

14

,

347

–

360

(

1984

).

https://doi.org/10.1016/0377-0257(84)80053-1

Google Scholar

Crossref

12.

J.

Lang

,

L.

Wang

, and

Q.

Wu

, “

Theoretical study of oscillating squeezing flow through a porous medium

,”

Trib. Int.

162

,

107110

(

2021

).

https://doi.org/10.1016/j.triboint.2021.107110

Google Scholar

Crossref

13.

C. Y.

Lee

and

C. Y.

Wen

, “

The oscillatory squeeze flow of electrorheological fluid considering the inertia effect

,”

Smart Mater. Struct.

11

,

553

–

560

(

2002

).

https://doi.org/10.1088/0964-1726/11/4/310

Google Scholar

Crossref

14.

N.

Phan-Thien

,

M.

Newberry

, and

R. I.

Tanner

, “

Non-linear oscillatory flow of a soft solid-like viscoelastic material

,”

J. Non-Newton Fluid Mech.

92

,

67

–

80

(

2000

).

https://doi.org/10.1016/S0377-0257(99)00110-X

Google Scholar

Crossref

15.

B.

Reynaud

and

T. M.

Quinn

, “

Anisotropic hydraulic permeability in compressed articular cartilage

,”

J. Biomech.

39

,

131

–

934

(

2006

).

https://doi.org/10.1016/j.jbiomech.2004.10.015

Google Scholar

Crossref

PubMed

16.

X.

Wang

,

M.

Gan

,

X.

Yang

,

P.

Zhang

,

X.

Peng

,

Y.

Ju

,

Y.

Kou

,

X.

Yu

,

L.

Zheng

, and

C.

Wang

, “

Mechanisms of enhanced oil recovery by fuzzy-ball fluid as a novel oil-displacement agent

,”

Energy Rep.

9

,

1447

–

1463

(

2023

).

https://doi.org/10.1016/j.egyr.2022.12.072

Google Scholar

Crossref

17.

T.

Karmakar

and

G. P.

Raja Sekhar

, “

Squeeze-film flow between a flat impermeable bearing and an anisotropic porous bed

,”

Phys. Fluids

30

,

043604

(

2018

).

https://doi.org/10.1063/1.5025882

Google Scholar

Crossref

18.

N.

Phan-Thien

and

R. I.

Tanner

, “

Lubrication squeeze film theory for the Oldroyd-B fluid

,”

J. Non-Newtonian Fluid Mech.

14

,

327

–

335

(

1984

).

https://doi.org/10.1016/0377-0257(84)80051-8

Google Scholar

Crossref

19.

K.

Rienesl

,

P. S.

Stelzer

, and

Z.

Major

, “

Squeeze flows rheometry and data analysis of carbon fiber sheet molding compounds

,”

Mater. Today: Proc.

62

,

2433

–

2435

(

2022

).

https://doi.org/10.1016/j.matpr.2022.02.583

Google Scholar

Crossref

20.

D. R.

Oliver

and

X.

Huang

, “

Squeeze film testing of ceramics pastes

,”

Br. Ceram. Trans.

99

,

101

–

108

(

2000

).

https://doi.org/10.1179/096797800680802

Google Scholar

Crossref

21.

R.

Arquier

,

G.

Miquelard-Garnier

,

I.

Iliopoulos

, and

G.

Régnier

, “

Assessing the shear viscous behavior of continuous carbon fiber reinforced PEKK composites with squeeze flow measurements

,”

Polym. Test.

123

,

108060

(

2023

).

https://doi.org/10.1016/j.polymertesting.2023.108060

Google Scholar

Crossref

22.

K.

Cho

,

G.

Rajan

,

P.

Farrar

,

L.

Prentice

, and

B. G.

Prusty

, “

Dental resin composites: A review on materials to product realizations

,”

Composites Part B

230

,

109495

(

2022

).

https://doi.org/10.1016/j.compositesb.2021.109495

Google Scholar

Crossref

23.

A.

Deblais

,

E.

den Hollander

,

C.

Boucon

,

A. E.

Blok

,

B.

Veltkamp

,

P.

Voudouris

,

P.

Verslus

,

H.-J.

Kim

,

M.

Mellema

,

M.

Stieger

,

D.

Bonn

, and

K. P.

Velikov

, “

Predicting thickness perception of liquid food products from their non-Newtonian rheology

,”

Nat. Commun.

12

,

6328

(

2021

).

https://doi.org/10.1038/s41467-021-26687-w

Google Scholar

Crossref

PubMed

24.

E.

Kucukal

,

Y.

Man

,

U. A.

Gurkan

, and

B. E.

Schmidt

, “

Blood flow velocimetry in a microchannel during coagulation using particle image velocimetry and wavelet-based optical flow velocimetry

,”

J. Biomech. Eng

143

,

091004

(

2021

).

https://doi.org/10.1115/1.4050647

Google Scholar

Crossref

PubMed

25.

A. D.

Rey

, “

Liquid crystals models of biological materials and processes

,”

Soft Matter

6

,

3402

–

3429

(

2010

).

https://doi.org/10.1039/b921576j

Google Scholar

Crossref

26.

A. D.

Rey

and

E. E.

Herrera-Valencia

, “

Liquid crystal models of biological material and silk spinning

,”

Biopolymers

97

,

374

–

396

(

2012

).

https://doi.org/10.1002/bip.21723

Google Scholar

Crossref

PubMed

27.

S. A.

Khadem

and

A. D.

Rey

, “

Thermodynamic modelling of acid collagenous solutions: From free energy contributions to phase diagrams

,”

Soft Matter

15

,

1833

(

2019

).

https://doi.org/10.1039/C8SM02140F

Google Scholar

Crossref

PubMed

28.

O. F.

Aguilar-Gutierrez

and

A. D.

Rey

, “

Theory and simulations of cholesteric film formation flows of dilute collagen solutions

,”

Langmuir

32

,

11799

(

2016

).

https://doi.org/10.1021/acs.langmuir.6b03443

Google Scholar

Crossref

PubMed

29.

M. A.

Nicolas

and

J. A.

Robbins

, “

The fluid mechanics of bolus ejection from the oral cavity

,”

J. Biochem.

34

,

1537

–

1544

(

2000

).

https://doi.org/10.1016/S0021-9290(01)00147-6

Google Scholar

Crossref

30.

S. J.

Lee

,

M. M.

Denn

,

M. J.

Crochet

, and

A. B.

Metzner

, “

Compressive flow between parallel disk: I. Newtonian fluid with a transverse viscosity gradient

,”

J. Non-Newtonian Fluid Mech.

10

,

3

–

30

(

1982

).

https://doi.org/10.1016/0377-0257(82)85002-7

Google Scholar

Crossref

31.

S. J.

Lee

,

M. M.

Denn

,

M. J.

Crochet

,

A. B.

Metz

, and

G. J.

Riggins

, “

Compressive flow between parallel disk: II. Oscillatory behavior of viscoelastic materials under a constant load

,”

J. Non-Newtonian Fluid Mech.

14

,

301

–

325

(

1984

).

https://doi.org/10.1016/0377-0257(84)80050-6

Google Scholar

Crossref

32.

N.

Phan-Thien

, “

Small strain oscillatory squeeze film flow of simple fluids

,”

J. Aust. Math. Soc. Ser. B

22

,

22

–

27

(

1980

).

https://doi.org/10.1017/S0334270000002514

Google Scholar

Crossref

33.

D. K.

Gartling

and

C. R.

Dohrmann

, “

Quadratic finite elements and incompressible viscous flows

,”

Comput. Methods Appl. Mech. Eng.,

195

,

1692

–

1708

(

2006

).

https://doi.org/10.1016/j.cma.2005.05.046

Google Scholar

Crossref

34.

N.

Phan-Thien

and

W.

Walsh

, “

Squeeze-film flow of an Oldroyd-B fluid: Similarity solution and limiting Weissenberg number

,”

Z. Angew. Math. Phys.

35

,

747

–

759

(

1984

).

https://doi.org/10.1007/BF00945440

Google Scholar

Crossref

35.

N.

Phan-Thien

,

J.

Dudek

, and

D. V.

Boger

, “

Squeeze film flow of ideal elastic liquids

,”

J. Non-Newtonian Fluid Mech.

18

,

227

–

254

(

1985

).

https://doi.org/10.1016/0377-0257(85)87001-4

Google Scholar

Crossref

36.

N.

Phan-Thien

,

F.

Sugeng

, and

R. I.

Tanner

, “

The squeeze film flow of a viscoelastic fluid

,”

J. Non-Newtonian Fluid Mech.

24

,

97

–

119

(

1987

).

https://doi.org/10.1016/0377-0257(87)85006-1

Google Scholar

Crossref

37.

N.

Phan-Thien

and

H. T.

Low

, “

Squeeze-film flow of a viscoelastic fluid a lubrication model

,”

J. Non-Newtonian Fluid Mech.

28

,

129

–

148

(

1988

).

https://doi.org/10.1016/0377-0257(88)85036-5

Google Scholar

Crossref

38.

N.

Phan-Thien

, “

Squeezing flow of a viscoelastic solid

,”

J. Non-Newtonian Fluid Mech.

95

,

343

–

362

(

2000

).

https://doi.org/10.1016/S0377-0257(00)00175-0

Google Scholar

Crossref

39.

N.

Phan-Thien

, “

N. Oscillatory torsional flow of a viscoelastic solid

,”

Comput. Mech.

29

,

143

–

150

(

2000

).

Google Scholar

Crossref

40.

F.

Avila

and

D. M.

Binding

, “

Normal and reverse squeezing flows

,”

J. Non-Newtonian Fluid Mech.

11

,

111

–

126

(

1982

).

https://doi.org/10.1016/0377-0257(82)85017-9

Google Scholar

Crossref

41.

D. R.

Oliver

and

R. C.

Ashton

, “

The flow of polymer-thickened oils in convergent jet thrust nozzles

,”

J. Non-Newtonian Fluid Mech.

2

,

367

–

384

(

1977

).

https://doi.org/10.1016/0377-0257(77)80022-0

Google Scholar

Crossref

42.

D. R.

Oliver

,

R. C.

Ashton

, and

G. D.

Wadelin

, “

The load bearing capacity of a continuous-flow squeeze film of liquid

,”

Appl. Sci. Res.

34

,

25

–

47

(

1978

).

https://doi.org/10.1007/BF00389274

Google Scholar

Crossref

43.

D. R.

Oliver

and

M.

Shahidullah

, “

Reverse squeeze film in a continuous flow system

,”

J. Non-Newtonian Fluid Mech.

15

,

3331

–

3339

(

1978

).

https://doi.org/10.1016/0377-0257(84)80017-8

Google Scholar

Crossref

44.

D. R.

Oliver

and

M.

Shahidullah

, “

Definitive load enhancement effects by polymer-thickened oils in a squeeze film experiment

,”

J. Non-Newtonian Fluid Mech.

9

,

257

–

267

(

1981

).

https://doi.org/10.1016/0377-0257(81)85004-5

Google Scholar

Crossref

45.

D. R.

Oliver

and

M.

Shahidullah

, “

Load enhancement effects by polymer thickened oils in strip squeeze film flow

,”

J. Non-Newtonian Fluid Mech.

13

,

93

–

102

(

1983

).

https://doi.org/10.1016/0377-0257(83)85024-1

Google Scholar

Crossref

46.

N. D.

Waters

and

D. K.

Gooden

, “

The flow of a power-law in a continuous-flow squeeze film

,”

Appl. Sci. Res.

40

,

169

–

179

(

1983

).

https://doi.org/10.1007/BF00386218

Google Scholar

Crossref

47.

N. D.

Waters

and

D. K.

Gooden

, “

The flow of an Oldroyd-B liquid in a continuous-flow squeeze film

,”

J. Non-Newtonian Fluid Mech.

14

,

361

–

376

(

1984

).

https://doi.org/10.1016/0377-0257(84)80054-3

Google Scholar

Crossref

48.

R. B.

Bird

,

R. C.

Armstrong

, and

O.

Hassager

,

Dynamics of Polymer Liquids

, Fluid Mechanics Vol. 1 (

Wiley

,

1987

).

Google Scholar

49.

G.

Larson

,

Constitutive Equations for Polymer Melts and Solutions

(

Butterworths

,

Boston, MA

,

1988

).

Google Scholar

50.

A. D.

Rey

, “

Capillary models for liquid crystals fibers, membranes, films, and drops

,”

Soft Matter

3

,

1349

–

1368

(

2007

).

https://doi.org/10.1039/b704248p

Google Scholar

Crossref

PubMed

51.

Z.

Wang

,

P.

Servio

, and

A.

Rey

, “

Geometry-structure models for liquid crystals interfaces, drops and membranes: Wrinkling, shape selection and dissipative shape evolution

,”

Soft Matter

19

,

9344

–

9364

(

2023

).

https://doi.org/10.1039/D3SM01164J

Google Scholar

Crossref

PubMed

52.

Z.

Wang

,

P.

Servio

, and

A.

Rey

, “

Pattern formation, structure and functionalities of wrinkled liquid crystals surfaces: A soft matter biomimicry platform

,”

Front. Soft Matter

3

,

1123324

(

2023

).

https://doi.org/10.3389/frsfm.2023.1123324

Google Scholar

Crossref

53.

E. E.

Herrera-Valencia

,

L. A.

Ramirez-Torres

,

C.

Soriano-Correa

,

M. L.

Sánchez-Villavicencio

,

O.

Bautista

,

V. J.

Hernández-Abad

, and

F.

Calderas

, “

Effects of multiple relaxation times in the annular flow of pulsatile electro-osmotic flow of a complex biological fluid: Blood with low and high cholesterol

,”

Front. Soft Matter

4

,

1385512

(

2024

).

https://doi.org/10.3389/frsfm.2024.1385512

Google Scholar

Crossref

54.

D. U.

Zamora-Cisneros

,

Z.

Wang

,

N. M.

Dorval Courchesne

,

M. J.

Harrington

, and

A. D.

Rey

, “

Geometric modeling of phase ordering for the isotropic smectic A phase transition

,”

Front. Soft Matter

4

,

1359128

(

2024

).

https://doi.org/10.3389/frsfm.2024.1359128

Google Scholar

Crossref

55.

A. D.

Rey

,

E. E.

Herrera-Valencia

, and

Y. K.

Murugesan

, “

Structure and dynamics of biological material liquid crystals

,”

Liq. Cryst.

41

,

430

–

451

(

2014

).

https://doi.org/10.1080/02678292.2013.845698

Google Scholar

Crossref

56.

L. A.

Ramírez-Torres

,

E. E.

Herrera-Valencia

,

M. L.

Sánchez-Villavicencio

,

C.

Soriano-Correa

,

V. J.

Hernández-Abad

, and

F.

Calderas

, “

Non-linear electrorheological model of a membrane immersed in Tanner-Power law fluids applied to outer hair cells: Shear-thinning mechanisms

,”

Phys. Fluids

36

,

033111

(

2024

).

https://doi.org/10.1063/5.0189731

Google Scholar

Crossref

57.

K.

Saha

,

P.-V.

Murthy

, and

S.

Chakraborty

, “

Effect of solvent depletion on electrokinetic energy conversion in viscoelastic fluids

,”

Phys. Fluids

36

,

062020

(

2024

).

https://doi.org/10.1063/5.0209974

Google Scholar

Crossref

58.

T. C. B.

McLeish

and

R. G.

Larson

, “

Molecular constitutive equations for a class of branched polymers: The Pom-Pom polymer

,”

J. Rheol.

42

,

81

–

110

(

1998

).

https://doi.org/10.1122/1.550933

Google Scholar

Crossref

59.

W. M. H.

Verbeeten

,

G. W. M.

Peters

, and

F. P. T.

Baaijens

, “

Differential constitutive equations for polymer melts: The extended Pom-Pom model

,”

J. Rheol.

45

,

823

–

843

(

2001

).

https://doi.org/10.1122/1.1380426

Google Scholar

Crossref

60.

L.

Abu-Farah

and

N.

Germann

, “

English translation of Giesekus's famous article on flows with constant velocity gradient and the motion of particles suspended therein. Part I. Spatial flows

,”

Phys. Fluids

36

,

031701

(

2024

).

https://doi.org/10.1063/5.0194331

Google Scholar

Crossref

61.

LL.

Ferras

and

A. M.

Afonso

, “

New insights into the extended and generalized PTT constitutive differential equations: Weak flows

,”

Fluid Dyn. Res.

55

,

035501

(

2023

).

https://doi.org/10.1088/1873-7005/acd0d4

Google Scholar

Crossref

62.

A. M.

Ribau

,

L. L.

Ferrás

,

M. L.

Morgado

,

M.

Rebelo

,

F. T.

Pinho

, and

A. M.

Afonso

, “

Analytical study of the annular flow of a generalized Phan-Thien-Tanner fluid

,”

Acta Mech.

235

,

1307

–

1317

(

2024

).

https://doi.org/10.1007/s00707-023-03784-z

Google Scholar

Crossref

63.

E. S.

Boek

,

J. T.

Padding

,

V. J.

Anderson

,

P. M. J.

Tardy

,

J. P.

Crawshaw

, and

J. R. A.

Pearson

, “

Constitutive equations for extensional flow of wormlike micelles: Stability analysis of the Bautista-Manero model

,”

J. Non-Newtonian Fluid Mech.

126

,

39

–

46

(

2005

).

https://doi.org/10.1016/j.jnnfm.2005.01.001

Google Scholar

Crossref

64.

E. E.

Herrera-Valencia

,

M. L.

Sánchez-Villavicencio

,

C.

Soriano-Correa

,

O.

Bautista

,

L. A.

Ramírez-Torres

,

V. J.

Hernández-Abad

, and

F.

Calderas

, “

Study of the electroosmotic flow of a structured fluid with a new generalized rheological model

,”

Rheol. Acta

63

,

3

–

32

(

2024

).

https://doi.org/10.1007/s00397-023-01418-8

Google Scholar

Crossref

65.

S.

Raczynski

, “

Extended prey–predator model

,” in

Catastrophes and Unexpected Behavior Patterns in Complex Artificial Populations

, Evolutionary Economics and Social Complexity Science Vol. 27 (

Springer

,

Singapore

,

2021

).

Google Scholar

Crossref

66.

Z. A.

Piskulich

,

O.

Mesele

, and

W. H.

Thompson

, “

Activation energies and beyond

,”

J. Phys. Chem. A

123

(

33

),

7185

–

7194

(

2019

).

https://doi.org/10.1021/acs.jpca.9b03967

Google Scholar

Crossref

PubMed

67.

Y.

Cai

,

J.

Qi

,

Y.

Lu

,

H.

He

, and

W.

Wu

, “

The in vivo fate of polymeric micelles

,”

Adv. Drug Deliv. Rev.

188

,

114463

(

2022

).

https://doi.org/10.1016/j.addr.2022.114463

Google Scholar

Crossref

PubMed

68.

R. A.

Abdel-Rahem

,

M.

Al-Remawi

,

A. Q.

Daraosheh

, and

H.

Hoffmann

, “

Rheological behavior of wormlike micelles (WLMs) in alcohol/water mixed solvent: Influence of alcohol chain length

,”

Colloid Polym. Sci.

299

,

1337

–

1351

(

2021

).

https://doi.org/10.1007/s00396-021-04852-3

Google Scholar

Crossref

69.

V. G.

Kulichikhin

and

A. Y.

Malkin

, “

The role of structure in polymer rheology: Review

,”

Polymers

14

,

1262

(

2022

).

https://doi.org/10.3390/polym14061262

Google Scholar

Crossref

PubMed

70.

P. A.

Vasquez

,

G. H.

McKinley

, and

L. P.

Cook

, “

A network scission model for wormlike micellar solutions: I. Model formulation and viscometry flow predictions

,”

J. Non-Newtonian Fluid Mech.

144

,

122

–

139

(

2007

).

https://doi.org/10.1016/j.jnnfm.2007.03.007

Google Scholar

Crossref

71.

C. J.

Pipe

,

N. J.

Kimp

,

P. A.

Vasquez

,

L. P.

Cook

, and

G. H.

McKinley

, “

Wormlike micellar solutions: II. Comparison between experimental data and scission model predictions

,”

J. Rheol.

54

,

881

–

913

(

2010

).

https://doi.org/10.1122/1.3439729

Google Scholar

Crossref

72.

N. A.

Spenley

,

X. F.

Yuan

, and

M. E.

Cates

, “

Non monotonic constitutive laws and the formation of shear banded flows

,”

J. Phys. II France

6

,

551

–

571

(

1996

).

https://doi.org/10.1051/jp2:1996197

Google Scholar

Crossref

73.

S.

Kheirandish

,

I.

Guybaidullin

,

W.

Wohlleben

, and

N.

Willenbacher

, “

Shear and elongational flow behavior of thickener solutions Part I: Effect of intermolecular aggregation

,”

Rheol. Acta

47

,

999

–

1013

(

2008

).

https://doi.org/10.1007/s00397-008-0292-1

Google Scholar

Crossref

74.

O.

Massarweh

and

A. S.

Abushaikha

, “

The use of surfactants in enhanced oil recovery: A review of recent advances

,”

Energy Rep.

6

,

3150

–

3178

(

2020

).

https://doi.org/10.1016/j.egyr.2020.11.009

Google Scholar

Crossref

75.

S.

Mihail

,

M.

Lyubov

,

M.

Denis

,

K.

Polina

,

B.

Sergei

, and

F.

Andrey

, “

Applicability assessment of viscoelastic surfactants and synthetic polymers as a base of hydraulic fracturing fluids

,”

Energies

15

,

2827

(

2022

).

https://doi.org/10.3390/en15082827

Google Scholar

Crossref

76.

Y.

Feng

,

G.

Zhaolong

,

Z.

Jinlong

, and

T.

Zhiyu

, “

Viscoelastic surfactant fracturing fluid for underground hydraulic fracturing in soft coal seams

,”

J. Pet. Sci. Eng.

169

,

646

–

653

(

2018

).

https://doi.org/10.1016/j.petrol.2018.06.015

Google Scholar

Crossref

77.

S.

Al-Hajri

,

B. M.

Negash

,

M.

Motiur Rahman

,

M.

Haroun

, and

T. M.

Al-Shami

, “

Perspective review of polymers as additives in water-based fracturing fluids

,”

ACS Omega

7

,

7431

–

7443

(

2022

).

https://doi.org/10.1021/acsomega.1c06739

Google Scholar

Crossref

PubMed

78.

X.

Li

,

X.

Zhang

,

L.

Wang

,

F.

Wen

,

Y.

Chen

,

Q.

Lv

,

H.

Ma

,

A.

Chen

,

R.

Wang

,

L.

Chen

,

Q.

Wang

,

D.

Dong

,

S.

Xu

, and

Q.

Niu

, “

Self-Assembled viscoelastic surfactant micelles with pH-responsive behaviour: A new fracturing-displacement integrated working fluid for unconventional reservoirs

,”

ACS Omega

9

(

21

),

22691

–

22702

(

2024

).

https://doi.org/10.1021/acsomega.4c00459

Google Scholar

Crossref

PubMed

79.

S.

Davoodi

,

M.

Al-Shargabi

,

D. A.

Wood

, and

V. S.

Rukavishnikov

, “

A comprehensive review of beneficial applications of viscoelastic surfactants in wellbore hydraulic fracturing fluids

,”

Fuel

338

,

127228

(

2023

).

https://doi.org/10.1016/j.fuel.2022.127228

Google Scholar

Crossref

80.

L. J.

Pyrak-Nolte

, “

Fracture specific stiffness: The critical link between the scaling behaviour of hydro-mechanical coupling in fractures and seismic monitoring

,” in

Science of Carbon Storage in Deep Saline Formations: Process Coupling Across Time and Spatial Scales

(

Springer

,

2019

), Chap. 14, pp.

311

–

335

.

Google Scholar

81.

M.

Marder

,

C. H.

Chen

, and

T.

Patzek

, “

Simple model of hydrofracture process

,”

Phys. Rev. E

92

,

062408

(

2015

).

https://doi.org/10.1103/PhysRevE.92.062408

Google Scholar

Crossref

82.

K.

Kamani

,

G. J.

Donley

, and

S. A.

Rogers

, “

Unification of the rheological physics of yield stress fluids

,”

Phys. Rev. Lett.

126

,

218002

(

2021

).

https://doi.org/10.1103/PhysRevLett.126.218002

Google Scholar

Crossref

PubMed

83.

E. E.

Herrera-Valencia

,

F.

Calderas

,

L.

Medina-Torres

,

M.

Pérez-Camacho

,

L.

Moreno

, and

O.

Manero

, “

On the pulsating flow behaviour of a biological fluid: Human blood

,”

Rheol. Acta

56

,

387

–

407

(

2017

).

https://doi.org/10.1007/s00397-017-0994-3

Google Scholar

Crossref

84.

D. C.

Vadillo

,

T. R.

Tuladhar

,

A. C.

Mulji

, and

M. R.

Mackley

, “

The rheological characterization of linear viscoelasticity for ink jet fluids using piezo axial vibrator and torsion resonator rheometers

,”

J. Rheol.

54

,

781

–

795

(

2010

).

https://doi.org/10.1122/1.3439696

Google Scholar

Crossref

85.

F.

Calderas

,

A.

Sánchez-Solís

,

A.

Maciel

, and

O.

Manero

, “

The transient flow of the PET-PEN-Montmorillonite clay nanocomposite

,”

Macromol. Symp.

283–284

,

354

–

360

(

2009

).

https://doi.org/10.1002/masy.200950942

Google Scholar

Crossref

86.

E. E.

Herrera

,

F.

Calderas

,

A. E.

Chávez

,

O.

Manero

, and

B.

Mena

, “Effect of random longitudinal vibrations on the Poiseuille flow of a complex fluid,”

Rheol. Acta

48

,

779

–

800

(

2009

).

https://doi.org/10.1007/s00397-009-0372-x

Google Scholar

Crossref

87.

E. E.

Herrera

,

F.

Calderas

,

A. E.

Chávez

, and

O.

Manero

, “

Study on the pulsating flow of a worm-like micellar solution

,”

J. Non-Newtonian Fluid Mech.

165

,

174

–

183

(

2010

).

https://doi.org/10.1016/j.jnnfm.2009.11.001

Google Scholar

Crossref

88.

E. E.

Herrera-Valencia

,

M. L.

Sanchez-Villavicencio

,

L.

Medina-Torres

,

D. M. Nuñez

Ramirez

,

Vicente Jesús

Hernández-Abad

,

F.

Calderas

, and

O.

Manero

, “New simple analytical method for flow enhancement predictions of pulsatile flow of a structured fluid,”

Phys. Fluids

31

,

063104

(

2019

).

https://doi.org/10.1063/1.5097867

Google Scholar

Crossref

Published open access through an agreement withCentro de Informacion Cientifica y Humanistica

Continuous flow of a thixotropic elasto-viscoplastic-banded structured fluid: Exponential structure rheological constitutive equation

I. INTRODUCTION

A. Squeeze flow film

B. Applications of the squeeze flow film

C. Mathematical and physical approaches

D. Continuous squeeze flow film

E. Complex fluids

F. Wormlike micellar structures

II. THEORETICAL EQUATIONS

A. Mass and momentum transfer equations

B. Exponential structural rheological constitutive equation (ESR)

C. Exponential kernel viscoelastic dissipation function

D. Particular cases

1. Small deformations of ESR model

2. Moderate and high shear strain deformations

III. PROBLEM FORMULATION

IV. SCALING RULES

A. Non-dimensional variables

B. Mass conservation

C. Non-dimensional groups

D. Continuity and momentum components

E. Total stress component

F. Rheological components

V. ASYMPTOTIC ANALYSIS

A. Regular perturbation scheme

B. Taylor series and linear approximation

C. Zeroth-order theory O (α0)

1. Material functions

2. Radial component of the field velocity

3. Volumetric flow rate

4. Axial component of the field velocity

5. Pressure gradient

6. Normal force

a. Lower shear strain

b. High shear strain

7. Yield stress

D. First-order theory O (α1)

1. Radial component of the field velocity

a. Low shear rate X

b. High shear rate

2. Axial component of the field velocity to first order

a. Low shear rate

b. High shear rate

3. Pressure gradient

a. Low shear rate

b. High shear rate

4. Normal force

a. Low shear rate

b. High shear rate

VI. SIMULATIONS AND RESULTS

A. Zeroth order mechanisms: O (α0)

1. Shear-thinning/thickening

2. Pseudo thixotropy/rheopexy

3. Yield stress

4. Wormlike micellar solutions data

B. Velocity profiles

C. Non-homogeneous flow

D. First order mechanisms: O (α1)

1. Elastic-shear-thinning/thickening

2. Elastic-pseudo thixotropy/rheopexy

3. Elastic-yield/stress

4. Elastic-concentration wormlike micellar data

VII. SUMMARY AND CONCLUSIONS

A. Force to zeroth order: O (α0)

B. Force to first order: O (α1)

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

C. Zeroth-order theory O (α⁰)

D. First-order theory O (α¹)

A. Zeroth order mechanisms: O (α⁰)

D. First order mechanisms: O (α¹)

A. Force to zeroth order: O (α⁰)

B. Force to first order: O (α¹)