Nonisothermal two‐ and three‐dimensional flow simulations of inelastic and viscoelastic fluids by a finite‐volume method

This paper applies the finite‐volume method to computations of steady flows of viscous and viscoelastic incompressible fluids in complex two and three‐dimensional geometries. The materials adopted in the study obey different constitutive laws: Newtonian, purely viscous Carreau–Yasuda as also Upper‐Convected Maxwell and Phan‐Thien/Tanner differential models, with a Williams–Landel–Ferry (WLF) equation for temperature dependence. Specific analyses are made depending on the rheological model. A staggered grid is used for discretizing the equations and unknowns. Stockage possibilities allow us to solve problems involving a great number of degrees of freedom, up to 1 500 000 unknowns with a desk computer. In relation to the fluid properties, our numerical simulations provide flow characteristics for various 2D and 3D configurations and demonstrate the possibilities of the code to solve problems involving complex nonlinear constitutive equations with thermal effects. Copyright © 2009 John Wiley & Sons, Ltd.     

Viscoelastic flow analysis using the software OpenFOAM and differential constitutive equations

Viscoelastic fluids are of great importance in many industrial sectors, such as in food and synthetic polymers industries. The rheological response of viscoelastic fluids is quite complex, including combination of viscous and elastic effects and non-linear phenomena. This work presents a numerical methodology based on the split-stress tensor approach and the concept of equilibrium stress tensor to treat high Weissenberg number problems using any differential constitutive equations. The proposed methodology was implemented in a new computational fluid dynamics (CFD) tool and consists of a viscoelastic fluid module included in the OpenFOAM, a flexible open source CFD package. Oldroyd-B/UCM, Giesekus, Phan-Thien–Tanner (PTT), Finitely Extensible Nonlinear Elastic (FENE-P and FENE-CR), and Pom–Pom based constitutive equations were implemented, in single and multimode forms. The proposed methodology was evaluated by comparing its predictions with experimental and numerical data from the literature for the analysis of a planar 4:1 contraction flow, showing to be stable and efficient.     

Velocity overshoots in gradual contraction flows

This study reports the results of a systematic numerical investigation, using the upper-convected Maxwell (UCM) and Phan-Thien–Tanner (PTT) models, of viscoelastic fluid flow through three-dimensional gradual planar contractions of various contraction ratios with the aim of investigating experimental observations of extremely large near-wall velocity overshoots in similar geometries [R.J. Poole, M.P. Escudier, P.J. Oliveira, Laminar flow of a viscoelastic shear-thinning liquid through a plane sudden expansion preceded by a gradual contraction, Proc. Roy. Soc. Lond. Ser. A 461 (2005) 3827]. We are able to obtain good qualitative agreement with the experiments, even using the UCM model in creeping-flow conditions, showing that neither inertia, second normal-stress difference nor shear-thinning effects are required for the phenomenon to be observed. Guided by the numerical results we propose a simple explanation for the occurrence of the velocity overshoots and the conditions under which they arise.     

Effect of viscoelasticity on the rotation of a sphere in shear flow

When particles are dispersed in viscoelastic rather than Newtonian media, the hydrodynamics will be changed entailing differences in suspension rheology. The disturbance velocity profiles and stress distributions around the particle will depend on the viscoelastic material functions. Even in inertialess flows, changes in particle rotation and migration will occur. The problem of the rotation of a single spherical particle in simple shear flow in viscoelastic fluids was recently studied to understand the effects of changes in the rheological properties with both numerical simulations [D’Avino et al., J. Rheol. 52 (2008) 1331–1346] and experiments [Snijkers et al., J. Rheol. 53 (2009) 459–480]. In the simulations, different constitutive models were used to demonstrate the effects of different rheological behavior. In the experiments, fluids with different constitutive properties were chosen. In both studies a slowing down of the rotation speed of the particles was found, when compared to the Newtonian case, as elasticity increases. Surprisingly, the extent of the slowing down of the rotation rate did not depend strongly on the details of the fluid rheology, but primarily on the Weissenberg number defined as the ratio between the first normal stress difference and the shear stress.In the present work, a quantitative comparison between the experimental measurements and novel simulation results is made by considering more realistic constitutive equations as compared to the model fluids used in previous numerical simulations [D’Avino et al., J. Rheol. 52 (2008) 1331–1346]. A multimode Giesekus model with Newtonian solvent as constitutive equation is fitted to the experimentally obtained linear and nonlinear fluid properties and used to simulate the rotation of a torque-free sphere in a range of Weissenberg numbers similar to those in the experiments. A good agreement between the experimental and numerical results is obtained. The local torque and pressure distributions on the particle surface calculated by simulations are shown.     

Observability of viscoelastic fluids

We apply the observability rank condition to study the observability of various viscoelastic fluids under imposed shear or extensional flows. In this paper the observability means the ability of determining the viscoelastic stress from the time history of the observations of the first normal stress difference. We consider four viscoelastic models: the upper convected Maxwell (UCM) model, the Phan–Thien–Tanner (PTT) model, the Johnson–Segalman (JS) model and the Giesekus model. Our study reveals that all of the four models have observability for all stress components almost everywhere under shear flow whereas under extensional flow most of the models have no observability for the shear stress component. More specifically, for UCM and JS models under imposed shear flow, the observations of the first normal stress difference allow the reconstruction of all components of viscoelastic stress. For UCM and JS models under extensional flow, the two normal stress components can be determined from the measurements of the first normal stress difference; the shear stress component does not affect the evolution of the normal stress components and consequently it cannot be extracted from the observations. Under shear flow, the PTT and Giesekus models have observability almost everywhere. That is, all components of the viscoelastic stress can be determined from the observations when the vector formed by the components of viscoelastic stress does not lie on a certain surface. Under extensional flow, the PTT model has observability almost everywhere for normal stress components whereas the Giesekus model has observability almost everywhere for all stress components. We also run simulations using the unscented Kalman filter (UKF) to reconstruct the viscoelastic stress from observations without and with noises. The UKF yields accurate and robust estimates for the viscoelastic stress both in the absence and in the presence of observation noises.     

Viscous dissipation with fluid inertia in oscillatory shear flow

For liquids with high viscosity and low thermal conductivity, viscous dissipation can cause appreciable errors in rheological property measurements. Here, the influences of both viscous dissipation and fluid inertia on the property measurements in oscillatory sliding plate rheometry are investigated. For Newtonian fluids, Bird (1965) solved the combined problem analytically, but only for high frequencies. Here his solution is extended to any frequencies. Also, the equations of motion and energy are solved for linear viscoelastic fluids, and new analytical solutions for the velocity and temperature profiles are given. In both Newtonian and linear viscoelastic fluids, the temperature rise in the gap increases with frequency. The location of the maximum temperature shifts from the mid-plane at low frequency towards the moving wall at high frequency. The fluid inertia increases the viscous dissipation in both fluids. By solving the combined problem, this paper simplifies rheometer design by providing one unified criterion for avoiding measurement errors. Operating limits are presented graphically for minimizing the effects of both fluid inertia and viscous dissipation in oscillatory sliding plate rheometry.     

Finite element computations of viscoelastic two-phase flows using local projection stabilization

A three-field local projection stabilized (LPS) finite element method is developed for computations of a three-dimensional axisymmetric buoyancy driven liquid drop rising in a liquid column where one of the liquid is viscoelastic. The two-phase flow is described by the time-dependent incompressible Navier-Stokes equations, whereas the viscoelasticity is modeled by the Giesekus constitutive equation in a time-dependent domain. The arbitrary Lagrangian-Eulerian (ALE) formulation with finite elements is used to solve the governing equations in the time-dependent domain. Interface-resolved moving meshes in ALE allows to incorporate the interfacial tension force and jumps in the material parameters accurately. A one-level LPS based on an enriched approximation space and a discontinuous projection space is used to stabilize the numerical scheme. A comprehensive numerical investigation is performed for a Newtonian drop rising in a viscoelastic fluid column and a viscoelastic drop rising in a Newtonian fluid column. The influence of the viscosity ratio, Newtonian solvent ratio, Giesekus mobility factor, and the Eötvös number on the drop dynamics are analyzed. The numerical study shows that beyond a critical Capillary number, a Newtonian drop rising in a viscoelastic fluid column experiences an extended trailing edge with a cusp-like shape and also exhibits a negative wake phenomena. However, a viscoelastic drop rising in a Newtonian fluid column develops an indentation around the rear stagnation point with a dimpled shape.     

Numerical study on the dynamics of droplet passing through a cylinder obstruction in confined microchannel flow

The droplet dynamics passing through a cylinder obstruction was investigated with direct numerical simulations with FE-FTM (Finite Element-Front Tracking Method). The effect of droplet size and capillary number (Ca) was studied for both Newtonian and viscoelastic fluids. In the case of Newtonian droplet immersed in Newtonian medium, the droplet breakup induced by the geometric hindrance depends on the droplet size. As Ca increases, the short droplets (1.3 times longer than the channel width) break up while passing through the obstruction. However, the breakup does not occur for longer droplets (1.8 times longer than the channel width). When the viscoelastic fluid characterized by the Oldroyd-B model is considered, the Newtonian droplet immersed in viscoelastic medium breaks up into two smaller droplets while passing through the cylinder obstruction with increasing De_m (Deborah number of the medium). We also show that the normal stress difference plays a key role on the droplet breakup and the droplet extension. The normal stress difference is enhanced in the negative wake region due to the droplet flow, which also promotes droplet extension in that region. This numerical study provides information not only on underlying physics of the droplet flows passing through a cylinder obstruction but also on the useful guidelines for microfluidic applications.     

Perturbation theory for viscoelastic fluids between eccentric rotating cylinders

The circumferential and radial profiles of velocity, pressure and stress are derived for the flow of model viscoelastic liquids between two slightly eccentric cylinders with the inner one rotating. Singular perturbation methods are used to derive expansions valid for small gaps between the cylinders, but for all Deborah numbers. Results for Newtonian, second-order, Criminale-Ericksen-Filbey, upper-convected Maxwell, and White-Metzner constitutive equation separate the effects of elasticity, memory, and shear thinning on the development of the large stress gradients that hinder numerical solutions with these models in more complicated geometries. The effect of the constitutive equation on the critical Deborah number for flow separation is presented.     

Viscoelastic flow between eccentric rotating cylinders: unstructured control volume method

This paper reports a convergent numerical algorithm for the Upper-Convected Maxwell (UCM) fluid between two eccentric cylinders at various eccentricity ratios (?); the outer cylinder is stationary, and the inner one rotating. The problem is solved by an unstructured control volume method (UCV), which is designed for a general viscoelastic flow problem with an arbitrary computational domain. A self-consistent false diffusion technique and an iteration scheme are used in combination to solve the problem. The computations of the UCM fluid using the numerical algorithm are carried out to a higher value of the Deborah number (De) at each eccentricity tested than hitherto possible with previous numerical simulations. The solutions are compared with previous numerical results, confirming the effectiveness of the UCV method as a general technique for solving viscoelastic flow problems.     

A high-resolution finite-difference method for simulating two-fluid,viscoelastic gel dynamics

An important class of gels are those composed of a polymer network and fluid solvent. The mechanical and rheological properties of these two-fluid gels can change dramatically in response to temperature, stress, and chemical stimulus. Because of their adaptivity, these gels are important in many biological systems, e.g. gels make up the cytoplasm of cells and the mucus in the respiratory and digestive systems, and they are involved in the formation of blood clots. In this study we consider a mathematical model for gels that treats the network phase as a viscoelastic fluid with spatially and temporally varying material parameters and treats the solvent phase as a viscous Newtonian fluid. The dynamics are governed by a coupled system of time-dependent partial differential equations which consist of transport equations for the two phases, constitutive equations for the viscoelastic stresses, two coupled momentum equations for the velocity fields of the two fluids, and a volume-averaged incompressibility constraint. We present a numerical method based on a staggered grid, second order finite-difference discretization of the momentum equations and a high-resolution unsplit Godunov method for the transport equations. The momentum and incompressibility equations are solved in a coupled manner with the Generalized Minimum Residual (GMRES) method using a multigrid preconditioner based on box-relaxation. We present results on the accuracy and robustness of the method together with an illustration of the interesting behavior of this gel model for the four-roll mill problem.     

A new finite element scheme using the Lagrangian framework for simulation of viscoelastic fluid flows

A new numerical scheme for simulation of viscoelastic fluid flows was designed, making use of finite element algorithms generally regarded as advantageous for tackling the problem. This includes the Lagrangian approach for the solution of viscoelastic constitutive equation using the co-deformational frame of reference with a possibility of analytically solving the equation along the particles trajectories, which in turn allowed eluding the solution of any system of linear equations for the stress. Then, the full ellipticity of the momentum conservation equation was utilised thanks to a possibility of accurate determination of the stress tensor independently of the velocity field at the current stage of computation. The needed independent stress was calculated at each time step on the basis of the past deformation history, which in turn was determined on the basis of the past velocity fields, all incorporated into a modified Euler time stepping algorithm. Owing to explicit inclusion of the full viscous term from the viscoelastic model into the momentum conservation equation, no stress splitting was necessary. The trajectory feet tracking was done accurately using a semi-analytic solution of the displacement gradient evolution equation and a weak formulation of the kinematics equation, the latter at the expense of solving an extra symmetric system of linear equations.The error expressed in the form of the Sobolev norms was determined using a comparison with available analytical solution for UCM fluid in the transient regime or numerically obtained steady-state stress values for the PTT fluid in Couette flow. The implementation of the PTT fluid model was done by modifying the relative displacement gradient tensor so that a new convective frame was defined.The stability of the algorithm was assessed using the well-known benchmark problem of a sphere sedimenting in a tube with viscoelastic fluid. The stable numerical results were obtained at high Weissenberg numbers, with the limit of convergence Wi=6.6, exceeding any previously reported values. The robustness of the code was proven by simulation of the Weissenberg effect (the rod-climbing phenomenon) with the use of PTT fluid.     

Constitutive equation with fractional derivatives for the generalized UCM model

A fundamental problem on the constitutive equation with fractional derivatives for the generalized upper convected Maxwell model (UCM) is studied. The existing investigations on the constitutive equation are reviewed and their limitations or deficiencies are highlighted. By utilizing the convected coordinates approach, a mathematically rigorous constitutive equation with fractional derivatives for the generalized UCM model is proposed, which has an explicit expression for the stress tensor. This model can be reduced to the linear generalized Maxwell model with fractional derivatives, the UCM model and some other existing models. In addition, the rheological properties of this proposed model in the start-up of simple shear and elongation flows are investigated. It is shown that this generalized UCM model can describe the various stress evolution processes and the strain hardening effect of the viscoelastic fluids.     

Selective withdrawal of polymer solutions: Computations

This paper reports numerical simulations of selective withdrawal of Newtonian and polymeric liquids, and complements the experimental study reported in the accompanying paper (Zhou and Feng [2]). We use finite elements to solve the Navier–Stokes and constitutive equations in the liquid on an adaptively refined unstructured grid, with an arbitrary Lagrangian–Eulerian scheme to track its free surface. The rheology of the viscoelastic liquids are modeled by the Oldroyd-B and Giesekus equations, and the physical and geometric parameters are matched with those in the experiments. The computed interfacial deformation is in general agreement with the experimental observations. In particular, the critical condition for interfacial rupture is predicted to quantitative accuracy. Furthermore, we combine the numerical and experimental data to explore the potential of selective withdrawal as an extensional rheometer. For Newtonian fluids, the measured steady elongational viscosity is within 47% of the actual value, apparently with better accuracy than other methods applicable to low-viscosity liquids. For polymer solutions, an estimated maximum error of 300% compares favorably with prior measurements.     

A 2D boundary element method for simulating the deformation of axisymmetric compound non‐Newtonian drops

The boundary integral formulation of the solution to the Stokes equations is used to describe the deformation of small compound non‐Newtonian axisymmetric drops suspended in a Newtonian fluid that is subjected to an axisymmetric flow field. The non‐Newtonian stress is treated as a source term in the Stokes equations, which yields an extra integral over the domains containing non‐Newtonian material. By transforming the integral representation for the velocity to cylindrical co‐ordinates and performing the integration over the azimuthal direction analytically, the dimension of the problem can be reduced from three to two. A boundary element method for the remaining two‐dimensional problem aimed at the simulation of the deformation of such axisymmetric compound non‐Newtonian drops is developed. Apart from a numerical validation of the method, simulation results for a drop consisting of an Oldroyd‐B fluid and a viscoelastic material are presented. Moreover, the method is extended to compound drops that are composed of a viscous inner core encapsulated by a viscoelastic material. The simulation results for these drops are verified against theoretical results from literature. Moreover, it is shown that the method can be used to identify the dominant break‐up mechanism of compound drops in relation to the specific non‐Newtonian character of the membrane. Copyright © 1999 John Wiley & Sons, Ltd.     

Viscoelastic effects in double‐pipe single‐pass counterflow heat exchangers

The conducto‐convective heat loss from a viscoelastic liquid, in the core of a double‐pipe heat exchanger arrangement, to a cooler Newtonian fluid flowing in the outer annulus is investigated with direct numerical simulations. A numerical algorithm based on the finite difference method is implemented in time and space with the Giesekus constitutive model for the viscoelastic liquids. The flow of both the annulus and core‐fluids is considered to be Poiseuille flow, driven by respective pressure gradients. In general, the results show that a viscoelastic core‐fluid leads to slightly lower (albeit comparable) attainable temperatures in the core‐fluid stream as compared with a corresponding Newtonian fluid. Copyright © 2008 John Wiley & Sons, Ltd.     

Displacement of yield-stress fluids in a fracture

We consider a displacement of several yield-stress fluids in a Hele-Shaw cell. The topic is relevant to the development of a model for the flow of multiple phases inside a narrow fracture with application to hydraulically fracturing a hydrocarbon-bearing underground formation. Existing models for fracturing flows include only pure power-law models without yield stress, and the present work is aimed at filling this gap. The fluids are assumed to be immiscible and incompressible. We consider fluid advection in a plane channel in the presence of density gradients. Gravity is taken into account, so that there can be slumping and gravitational convection. We use the lubrication approximation so that governing equations are reduced to a 2D width-averaged system formed by the quasi-linear elliptic equation for pressure and transport equations for volume concentrations of fluids. The numerical solution is obtained using a finite-difference method. The pressure equation is solved using an iterative algorithm and the Multigrid method, while the transport equations are solved using a second-order TVD flux-limiting scheme with the superbee limiter. This numerical model is validated against three different sets of experiments: (i) gravitational slumping of fluids in a closed Hele-Shaw cell, (ii) viscous fingering of fluids with a high viscosity contrast due to the Saffman–Taylor (S–T) instability in a Hele-Shaw cell at microgravity conditions, (iii) displacement of Bingham fluids in a Hele-Shaw cell with the development of fingers due to the S–T instability. Good agreement is observed between simulations and laboratory data. The model is then used to investigate the joint effect of fingering and slumping. Numerical simulations show that the slumping rate of yield-stress fluid is significantly less pronounced than that of a Newtonian fluid with the same density and viscosity. If a low-viscosity Newtonian fluid is injected after a yield-stress one, the S–T instability at the interface leads to the development of fingers. As a result, fingers penetrating into a fluid with a finite yield stress locally decrease the pressure gradient and unyielded zones develop as a consequence.     

A FENE-P k–ε turbulence model for low and intermediate regimes of polymer-induced drag reduction

A new low-Reynolds-number k–ε turbulence model is developed for flows of viscoelastic fluids described by the finitely extensible nonlinear elastic rheological constitutive equation with Peterlin approximation (FENE-P model). The model is validated against direct numerical simulations in the low and intermediate drag reduction (DR) regimes (DR up to 50%). The results obtained represent an improvement over the low DR model of Pinho et al. (2008) [A low Reynolds number k–ε turbulence model for FENE-P viscoelastic fluids, Journal of Non-Newtonian Fluid Mechanics, 154, 89–108]. In extending the range of application to higher values of drag reduction, three main improvements were incorporated: a modified eddy viscosity closure, the inclusion of direct viscoelastic contributions into the transport equations for turbulent kinetic energy (k) and its dissipation rate, and a new closure for the cross-correlations between the fluctuating components of the polymer conformation and rate of strain tensors (NLT_ij). The NLT_ij appears in the Reynolds-averaged evolution equation for the conformation tensor (RACE), which is required to calculate the average polymer stress, and in the viscoelastic stress work in the transport equation of k. It is shown that the predictions of mean velocity, turbulent kinetic energy, its rate of dissipation by the Newtonian solvent, conformation tensor and polymer and Reynolds shear stresses are improved compared to those obtained from the earlier model.     

Numerical simulation of non-linear elastic flows with a general collocated finite-volume method

This paper reports the development and application of a finite-volume based methodology for the calculation of the flow of fluids which follow differential viscoelastic constitutive models. The novelty of the method lies on the use of the non-staggered grid arrangement, in which all dependent variables are located at the center of the control volumes, thus greatly simplifying the adoption of general curvilinear coordinates. The pressure–velocity–stress decoupling was removed by the development of a new interpolation technique inspired on that of Rhie and Chow, AIAA 82 (1982) 998. The differencing schemes are second order accurate and the resulting algebraic equations for each variable are solved in a segregated way (decoupled scheme). The numerical formulation especially designed for the interpolation of the stress field was found to work well and is shown to be indispensable for accurate results. Calculations have been carried out for two problems: the entry flow problem of Eggleton et al., J. Non-Newtonian Fluid Mech. 64 (1996) 269, with orthogonal and non-orthogonal meshes; and the bounded and unbounded flows around a circular cylinder. The results of the simulations compare favourably with those in the literature and iterative convergence has been attained for Deborah and Reynolds numbers similar to, or higher than, those reported for identical flow problems using other numerical methods. The application of the method with non-orthogonal coordinates is demonstrated. The entry flow problem is studied in more detail and for this case differences between Newtonian and viscoelastic fluids are identified and discussed. Viscoelasticity is shown to be responsible for the development of very intense normal stresses, which are tensile in the wall region. As a consequence, the viscoelastic fluid is more intensely decelerated in the wall region than the Newtonian fluid, thus reducing locally the shear rates and the role of viscosity in redeveloping the flow. A layer of high stress-gradients is formed at the wall leading edge and is convected below and away from the wall; its effect is to intensify the aforementioned deviation of elastic fluid from the wall.     

Influence of Heat and Mass Transfer on Newtonian Biomagnetic Fluid of Blood Flow Through a Tapered Porous Arteries with a Stenosis

Heat and mass transfer effects on Newtonian biomagnetic fluid of blood flow through a tapered porous artery with a stenosis is investigated. Governing equations have been modeled by treating blood as Newtonian biomagnetic fluid. The governing equations are simplified under the assumption of mild stenosis. Exact solutions have been evaluated for velocity, temperature, and concentration profiles. The effects of Newtonian nature of blood on velocity, temperature, concentration profile, wall shear stress, shearing stress at the stenosis throat and impedance of the artery are discussed graphically. Stream lines have been presented in last section of the article.