This paper presents a novel implementation of a numerical scheme for predicting complex flows of viscoelastic fluids using a finitely extensible nonlinear elastic (FENE) transient network model. This model extends the FENE model by incorporating chain interactions and accounting for the way in which the maximum chain length, drag, and relaxation time are influenced by entanglement and disentanglement processes. Three different initial networks are considered (disentanglement, entanglement, and aleatory), and the influence of variables such as the kinetic rate constants, elasticity, and chain length on the microstate concentration, stresses, and drag force is investigated. It is shown that although the concentrations of the microstates are independent of the Weissenberg number and the maximum extension length, the stresses and hence the drag are influenced by them.
I. INTRODUCTION
The flow around a cylinder confined between two plates (see Fig. 1) and the flow around a sphere in a tube are two of the benchmark problems in computational rheology. The problems are used to compare various discretization methods and numerical solution techniques in terms of accuracy and stability and also to gain an improved understanding of some of the phenomena associated with viscoelastic fluids in complex flows. An early review of these two benchmark problems can be found in the monograph of Owens and Phillips.1 Most of the contributions to the development of computational models for viscoelastic flows have been based on a description of viscoelastic fluids using closed form constitutive equations such as the Oldroyd B, Giesekus, Phan-Thien Tanner (PTT), and finitely extensible nonlinear elastic (FENE)-P models. Stabilization techniques such as streamline upwind Petrov–Galerkin (SUPG)2 or a reformulation of the constitutive equation such as the discrete elastic viscous split stress (DEVSS) scheme3 or using the logarithm of the conformation tensor4 are often used to extend the range of the Weissenberg numbers for which converged approximations are obtained.
Phan-Thien and Dou5 developed a parallelized unstructured finite volume method (FVM) in conjunction with the elastic viscous split stress (EVSS) formulation for this problem. They obtained drag reduction for both the simplified PTT and upper convective Maxwell (UCM) constitutive models. In a subsequent paper, Dou and Phan-Thien6 studied an instability that occurs in the flow of an Oldroyd-B fluid past a circular cylinder in a channel using their parallelized FVM. They showed that the instability features stress oscillations on the top of the cylinder, which occurs at a critical value of the Deborah number in agreement with experiments. Claus and Phillips7 considered the problem of flow of a viscoelastic fluid (Oldroyd-B and Giesekus) around a confined cylinder using high-order spectral/hp element methods [a class of methods that combines both mesh refinement (h) with higher polynomial order (p)] and explored the causes for the breakdown in numerical convergence beyond a critical value of the Weissenberg number. They predicted the instability identified by Dou and Phan-Thien6 at the same critical value of the Weissenberg number.
Hulsen et al.8 implemented the log-conformation formulation of the constitutive equation within a finite element method using the DEVSS/Discontinuous Galerkin (DEVSS/DG) approach. They performed a 1D stability analysis to identify the failure of their original numerical scheme to balance exponential growth, which was a possible source of numerical instabilities at high Weissenberg numbers. The new version of the model generated results for an extensive range of Weissenberg numbers for the Oldroyd B and Giesekus models.
Oliveira and Miranda9 considered two-dimensional, inertia-free flow of a constant viscosity viscoelastic fluid and investigated the conditions under which the flow becomes unsteady using FVM. They showed that the failure to obtain converged steady solutions for moderate values of the Deborah number is due to the separation of the 2D flow that becomes unsteady and is characterized by a small pulsating recirculation region behind the cylinder. In a later publication, Afonso et al.10 investigated the uniform steady flow of viscoelastic fluids past a cylinder confined between two parallel moving plates using FVM for five constitutive models [UCM, Oldroyd-B, FENE-Chilcott Rallison (FENE-CR), PTT, and Giesekus]. Using very fine meshes, particularly in the wake of the cylinder where large normal stresses are generated at high Deborah numbers, they obtain converged solutions for Deborah numbers that were in excess of those that had been reported in the literature at that time.
Liu and Nithiarasu11 used a fully explicit, characteristic-based split (CBS) method for the flow of an Oldroyd B fluid past a cylinder. The combination of fractional solution stages used in the CBS scheme and the higher-order time step-based convection stabilization reduced the instability that was predicted above a Deborah number of around 0.6. They showed that by adding artificial diffusion to the discrete constitutive equation, the conformation tensor remains positive definite allowing higher Deborah numbers to be reached compared with other numerical schemes.
Vargas et al.12 performed numerical simulations of the viscoelastic flow past a cylinder in a channel and a sphere in a cylinder using a micro-macro technique and kinetic models. The basis of the numerical method employed is a micro-macro model in which the polymer dynamics is described by the evolution of an ensemble of Brownian configuration fields. The macroscopic equations are solved using an extension of the spectral element method developed by Gerritsma and Phillips13 for Stokes flow. The micro-macro approach can be used to simulate models that do not possess an equivalent closed form representation.
Ribeiro et al.14 performed a combined experimental and numerical study on the flow of a Newtonian fluid around a confined cylinder placed in a rectangular duct to assess 3D effects on the flow patterns. The flow visualizations were carried out using streak photography and the velocity measurements by particle image velocimetry for different aspect ratios (length to diameter of the cylinder) and Reynolds numbers until the flow becomes time-dependent. The numerical calculations were performed on 3D meshes using a finite volume scheme. For large aspect ratios, unexpected velocity peaks near the cylinder end walls downstream of the cylinder were observed, which were not diminished when the aspect ratio was increased. Reducing the aspect ratio eliminated flow separation as expected for Hele–Shaw type flows.
Garduño et al.15 present a new dissipative viscoelastic model in order to capture the levels of drag enhancement measured in experiments on flow past a sphere using Boger fluids. They showed that, for viscoelastic fluids with a large solvent component, an initial reduction in drag is followed by drag enhancement when hybrid dissipative models based on the exponential PTT (ePTT) and FENE-CR models are used. The dissipative ePTT model displays shear thinning and strain hardening/softening behavior while the dissipative FENE-CR model displays constant shear viscosity and strain hardening behavior that attains a plateau.
Peng et al.16 have studied two-dimensional viscoelastic flow past two side-by-side circular cylinders at a Reynolds number of 100 using the FENE-P model. The aim of this study is to provide new insight into multi-body dynamics. They studied the transitions in flow behavior as the Weissenberg number and the spacing between the cylinders is varied. In all cases, the drag on the cylinders increases with increasing Weissenberg number, while the repulsive force between the two cylinders gradually decreases for a lower spacing and increases for a higher spacing.
Ruan and Ouyang17 introduced a collocated FVM on unstructured meshes to simulate the viscoelastic flow of polymer melts with viscous dissipation past a confined cylinder, using as a constitutive equation a non-isothermal Peterlin approximation of the finitely extensible nonlinear elastic dumbbell (FENE-P) model.
Norouzi et al.18 performed a numerical investigation of the laminar viscoelastic shedding flow behind inclined square cylinders using the Giesekus model and a parallelized finite volume method. The Reynolds number is increased from Re = 60 to Re = 120 and different angles of incidence are considered. Their computational results indicate that the flow is destabilized either by increasing the Reynolds number or by increasing the Weissenberg number. Increasing polymer concentration accentuates the increase in the amplitude of the lift coefficient and shedding frequency.
This paper presents the complex flow simulations of viscoelastic fluids using a FENE transient network model19 to describe the polymer dynamics. The model can be viewed as an extension of the FENE model to more concentrated polymer solutions. Chain interactions are accounted for in the FENE transient network model, whereas they are ignored in the FENE model for dilute polymer solutions. As the polymer concentration increases so does the maximum chain length, drag, and relaxation time due to the processes of entanglement and disentanglement. These are accounted for in the transient network model. Network models are characterized by creation and destruction processes that determine the degree of entanglement of polymer chains and how these change in response to flow and changes in energy. Sim et al.20 studied the possible mechanisms of a network behavior composed of segments and junctions. Classical transient network models consider the creation of microstates as a thermal activation process, while flow destroys the network.21–23 Rincón et al.24 proposed the formation of a transient network through five microstates, using the moments of the distribution function to calculate the rheological functions.25 In the classical FENE model,26 the nonlinear dumbbell possesses a constant maximum extensibility. However, in the transient network model,24 the maximum extensibility is not a constant but a variable resulting from a kinetic process.
Rincón et al.24 show that the transient network model is able to predict effects observed in associative polymers and wormlike surfactants that cannot be reproduced by classical network models. For example, the model predicts a power-law region between Newtonian plateau regions at low and high shear rates. Furthermore, the model predicts a shear-thickening regime between the first Newtonian plateau and the shear-thinning region. In this regime, the rate at which structure is rebuilt is shorter than the characteristic time of the flow and so a flow-induced structure endures over a range of shear rates. At larger shear rates, the rate at which deformation of the network occurs is greater than the rate of reformation and this leads to shear-thinning. These effects agree with the predictions of Vaccaro and Marrucci.23 In predictions of stress relaxation after the cessation of steady shear flow, very close agreement is obtained concerning the stretched exponential region with experimental data for wormlike micelles (see Berret et al.,27 for example). Subsequent developments of the transient network model have included a study of the rheological behavior of micellar solutions28 and a coupling of the model to a rheological constitutive equation of state within an irreversible thermodynamics framework.29
Ferrer et al.19 coupled the FENE model with the transient network of Rincón et al.24 to describe the behavior of a complex fluid under simple shear flow. The novel contribution of this current paper is the development and implementation of a micro-macro numerical scheme for a FENE transient network model. This model is used to predict the complex flow around a confined cylinder (planar) and sphere (axisymmetric). Although many complex flow simulations have been performed using the FENE model and its variants,12,15,17 the numerical predictions presented in this paper are the first simulations based on the FENE transient network model. This is the first complex flow numerical simulation using this model for a benchmark problem in computational rheology. The effect of the maximum extension length, kinetic rate constants, and elasticity on the drag force and the shear and normal stresses is investigated and analyzed.
The rest of the paper is structured as follows. Section II introduces the mathematical formulation of the governing equations for the hybrid microscopic-macroscopic model and provides a description of the FENE transient network model. Section III describes the non-dimensionalization of the governing equations. Section IV provides a description of the spatial and temporal discretizations of these equations. Section V provides a mathematical description of the benchmark problems. Section VI presents results for planar flow past a cylinder in a channel and axisymmetric flow past a sphere in a cylinder. Concluding remarks are made in Sec. VII.
II. MATHEMATICAL MODEL
In this paper, a hybrid micro-macro model is used as the basis for numerical predictions since the approach provides greater flexibility in the modeling of polymeric fluids. Not all kinetic theory models have closed form macroscopic equivalents, and the micro-macro approach circumvents the need to implement closure approximations that are not universally accurate. We assume that the fluid is incompressible.
A. Macroscopic Equations
B. Microscopic equations
C. Transient network model
Microstate . | Number of nodes . | Number of segments . | Number of chains . | . |
---|---|---|---|---|
ω0 | 0 | 1 | 1 | 1 |
ω1 | 1 | 4 | 2 | |
ω2 | 2 | 7 | 3 | |
ω3 | 3 | 9 | 3 | |
ω4 | 4 | 12 | 4 |
Microstate . | Number of nodes . | Number of segments . | Number of chains . | . |
---|---|---|---|---|
ω0 | 0 | 1 | 1 | 1 |
ω1 | 1 | 4 | 2 | |
ω2 | 2 | 7 | 3 | |
ω3 | 3 | 9 | 3 | |
ω4 | 4 | 12 | 4 |
III. DIMENSIONLESS EQUATIONS
IV. NUMERICAL SCHEME
A. Macroscopic equations
B. Microscopic equations
The Brownian configuration fields method (BCFM) is a method for solving stochastic partial differential equation using a variance reduction technique.12,36,37 The main idea of BCFM is to assume that the same initial ensemble of dumbbells is located at every spatial discretization point instead of independent ensembles, which enables the same sequence of random numbers to be used to build the trajectories of many groups of dumbbells. The conformation fields are denoted by , with .
V. PROBLEM DESCRIPTION
A. Initial conditions
B. Boundary conditions
VI. NUMERICAL RESULTS
The numerical results obtained were generated using a code developed by the group of Phillips at the School of Mathematics at Cardiff University. The transient network model developed by Rincón et al.24 and improved by Ferrer et al.19 and Gómez-López et al.39 was incorporated into this numerical code. This represents the novel addition to this software from the current contribution.
A. Validation—Flow past a cylinder
To validate the numerical scheme, the benchmark problem of flow past a cylinder is considered, and a comparison is made with results in the literature. In this benchmark problem, fluid flows past a cylinder confined symmetrically in a channel. The dimensionless drag coefficient is calculated for a Newtonian fluid with Re = 0.01 and used to evaluate the accuracy of the numerical approximation. For example, Hulsen et al.8 calculated a drag of F = 132.358, which is in excellent agreement with the value of F = 132.351 determined by Vargas et al.12 and reproduced in this paper.
B. FENE model
The FENE model is used as a starting point to analyze the effect of the discretization parameters on the evolution of the drag on the cylinder for a viscoelastic fluid with Re = 0.01, , , b = 50 and the number of spectral elements remains constant. Since the transient network model has been developed for concentrated polymer solutions, a value of was chosen for the viscosity ratio throughout this numerical study since these fluids are characterized by small values of β. The dependence of the evolution of the drag on the cylinder on polynomial order N for Nf = 2000, and the number configuration fields Nf for N = 6 is shown in Figs. 4(a) and 4(b), respectively. Increasing the value of N with Nf = 2000 results in minor variability in the drag as can be seen in Fig. 4(a), and therefore, N = 6 is used in subsequent computations. In Fig. 4(b), we observe that increasing Nf reduces the fluctuations in the evolution of the drag and that there is not a significant improvement beyond Nf = 2000 and so this value is used in subsequent computations. For further results for the FENE, FENE-P, and Hookean dumbbell models, the reader is referred to Vargas et al.12 where mesh convergence of the evolution of the drag and the stress components along the axis and around the cylinder/sphere are presented. This validates the numerical solution procedure in the case of viscoelastic flows.
C. FENE transient network model
In this subsection, the numerical simulation of a complex fluid around a cylinder using the FENE-transient network model is analyzed for three networks with different initial conditions. The initial conditions are: (i) disentanglement, characterized by and for , (ii) entanglement, characterized by and for , and (iii) aleatory, where all the microstates are present ( for ). In addition to the mesh convergence studies showing the convergence of the drag force using the FENE model, we present a corresponding study for the FENE network model. In Fig. 5, we show the dependence of the evolution of the drag force on the cylinder on polynomial order N in the case of the disentanglement configuration with Re = 0.01, , b = 50, A = 1, B = 1, and Nf = 2000. This figure demonstrates that convergence with increasing polynomial order has been achieved and the close agreement between the results generated with N = 6 and N = 8 validates the choice of N = 6 used for the rest of the paper.
1. Influence of the kinetic rate constants A and B
The influence of the kinetic rate constants A and B on the three different initial entanglement networks, from highly structured (A = 100, B = 1) to weakly structured (A = 1, B = 100) on the principal variables, is analyzed. First, Fig. 6 shows the evolution of the concentration of the microstates for each of the starting conditions that were defined previously, i.e., disentanglement, entanglement, and aleatory. The network with disentanglement initial conditions is the most sensitive to changes in the kinetic rate constants. For this case, the available energy in the system is used to produce more entangled structures, as can be seen when B is kept constant and A is changed. Rincón et al.24 have shown that at low shear rates for both highly structured (A = 100, B = 1) and weakly structured (A = 1, B = 100) systems, l(t), the non-dimensional distance between nodes, is close to 1/3, which corresponds to a more entangled microstate ω4 (see Table I). In Fig. 6, the concentrations of microstates ω2 and ω4 increase in time. There is negligible effect on the evolution of the microstate concentrations in varying B keeping A = 1. For the network with entanglement initial conditions, the change in concentrations is very small for all the values of A and B considered, which indicates that a highly structured network requires more energy to produce simpler microstructures. Finally, for the network with aleatory initial conditions, the concentrations of the microstates exhibit minor modifications, and for all values of A and B, the behavior is basically the same. For this initial condition, the microstates are close to an equilibrium concentration, which requires the minimum amount of energy.
The effect of changes to the microstructural kinetics on the shear and normal stresses for networks characterized by the first two initial conditions, disentanglement and entanglement, is shown in Fig. 7. The network with disentanglement initial conditions is the most sensitive to changes in the kinetic rate constants giving rise to higher stresses. For the network with entanglement initial conditions, no significant changes are found. The drag force experienced by the cylinder is an overall macroscopic indicator of the different microscopic processes taking place during the flow of a complex fluid around the cylinder. Figure 8 shows the influence of the microstructural kinetics on the evolution of the drag force. For an initially disentangled network, the energy available from the system is used to produce more complex structures (see Fig. 6) increasing the shear and normal stresses and hence the drag force. For the entanglement and aleatory networks, there are no significant changes. Since there is very limited rebuilding of the structure in these networks, the corresponding stresses and drag force are smaller than for the initially disentangled configuration. For both cases, modifying the kinetic rate constants generates average drag values of approximately 50 and 45, respectively. For the aleatory network, there was a loss of convergence for the case A = 1 and B = 100 at , and this is not shown in the figure.
2. Effect of the elasticity (Wi)
In this subsection, the effect of elasticity (Wi) on the main variables is investigated for networks with the three different initial conditions with: A = 1, B = 1, Re = 0.01, , b = 50, N = 6, and Nf = 2000. Figure 9 presents the influence of Wi on the evolution of all microstates concentrations. No significant changes in the concentrations for any of the different initial network configurations are observed, i.e., the evolution of the microstructural kinetics appears to be independent of elasticity. Note that the expression for the polymeric stress in Eq. (26) involves the reciprocal of Wi, and therefore, one would expect a decrease in stress with increasing Wi as shown in Fig. 10. This is the case for the disentanglement and entanglement networks. The largest absolute values of stress are found for the disentanglement network. For the aleatory network, the stress profiles are not well-defined. Figure 11 shows the effect of elasticity on the evolution of the drag force. It is confirmed that increasing Wi decreases the drag force uniformly for the networks with disentanglement or entanglement initial conditions. For the network with disentanglement initial conditions, there is a sharp decrease in the drag force as Wi is increased from Wi = 1 to Wi = 2. This behavior has been reported by Vargas et al.12 In the aleatory network, the trend is not exactly the same, but in general terms, there is a reduction in the drag force. Further simulations were performed for higher Weissenberg numbers (up to Wi = 20) for the initial disentanglement configuration and no convergence problems were experienced. The dependence of the evolution of the drag force on Wi is shown on a log –log scale in Fig. 12. The drag decreases monotonically with respect to increasing Wi. The steady-state value of the drag is reached by t = 2 for . For Wi > 8, the steady-state drag value has not quite been reached by t = 5.
Since the flow past a cylinder is a benchmark problem in computational rheology, there is much available data showing the dependence of drag on Wi for different viscoelastic models. For the Oldroyd B model, extremely fine meshes are required to obtain convergence at the maximum attainable Weissenberg number for .7,10,40 In this, the drag decreases monotonically to around Wi = 0.7 before increasing ever so slightly. However, the lack of convergence of the stress in the wake of the cylinder for values of Wi much beyond this value means that the predictions are not useful or reliable.8 Afonso et al.10 showed that for viscoelastic models that exhibit shear-thinning behavior such as FENE-CR, PTT, and Giesekus, converged solutions can be obtained over a much larger range of values of the Deborah number, De, since the stress in the wake is weaker. For these models, the drag decreases monotonically with De. The steady-state predictions of drag shown in Fig. 11 are in agreement with this trend.
3. Effect of the maximum extension length b
The effect of the maximum extension length b for networks with different initial conditions with: A = 1, B = 1, Re = 0.01, Wi = 1, and on the main variables is presented. Figure 13 shows the influence of b on the evolution of all microstates concentrations. There are no appreciable changes in any of the entanglement scenarios, i.e., the microstate kinetics are also independent of extension length. In terms of magnitude, the disentangled network shows a higher evolution of all microstructure concentrations. This is in contrast to the entangled network, where the concentration of more entangled structure dominates. The effect of extension length is manifested on the shear and normal stresses for different initial network configurations as shown in Fig. 14. Increasing b results in an increase in the stresses presenting a greater magnitude in the disentangled network. For large values of the extensibility parameter, the behavior of the simplest microstructures described by the FENE model behaves like the Hookean model.12,25,30,36 Figure 15 shows the effect of b on the drag force. Drag fluctuations increase as b is smaller, due to the small molecular chains and the increased tension between the beads. The drag force increases with increasing b; a behavior that is present in all cases, with the largest drag attained for the network with the initial disentanglement configuration, confirming that an entanglement network requires more energy to destroy it.
D. Flow past a sphere
VII. CONCLUSIONS
This paper presents the first complex flow simulations of viscoelastic fluids using a FENE transient network model24 to describe the polymer dynamics. The model can be viewed as an extension of the FENE model to more concentrated polymer solutions for which chain interactions cannot be ignored as they are in the dilute regime. As the polymer concentration increases, so does the maximum chain length, drag, and relaxation time due to the processes of entanglement and disentanglement. These are accounted for in the transient network model. The numerical simulation of the flow of a complex fluid past a confined object is considered. In particular, two benchmark problems in computational rheology are considered, viz., flow past a cylinder and flow past a sphere. The blockage ratio is 2:1 in both cases. The spatial discretization is based on the spectral element method. This method is used to solve both the macroscopic governing equations and the microscopic equations based on the method of Brownian configuration fields that serves to determine the polymeric contribution to the extra stress tensor. The numerical approach is validated by making comparisons with the literature for the Oldroyd B and FENE models, and the convergence of the drag force is investigated with respect to polynomial order and the number of configuration fields. This information is used as the basis for choosing the discretization parameters for the remainder of the paper. For the FENE transient network model, three different initial network configurations were considered: disentanglement, entanglement, and aleatory. The influence of the main variables on these different initial networks is summarized as follows:
-
The influence of the kinetic rate constants A and B is strongly dependent on the initial network configurations, from highly structured (A = 100, B = 1) to weakly structured (A = 1, B = 100). The disentanglement network is the most sensitive network to changes in the kinetic rate constants leading to larger stresses, and a correspondingly larger drag force. The change in concentrations is very small for the entanglement network, which indicates that a highly structured network requires more energy to produce simpler microstructures. Finally, in the case of the aleatory network, the concentrations of the microstates exhibit minor modifications, since the microstates are close to an equilibrium concentration.
-
The evolution of the microstate concentrations is independent of elasticity. The influence of the elasticity is manifested on the shear and normal stresses; it is also confirmed that increasing Wi decreases the drag force.
-
The microstate kinetics are also independent of the maximum extension length. The effect of the extension length is manifested on the shear and normal stresses for different networks and drag fluctuations increase as b is reduced.
The same trends are observed in the investigation of flow past a sphere.
ACKNOWLEDGMENTS
A. Gómez-López gratefully acknowledges the projects CI2270 and PE1039323 and FES Cuautitlán UNAM for financial support. R. O. Vargas would like to acknowledge the financial support provided by the project SIP-20241726.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
A. Gómez-López: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing—original draft (equal); Writing—review & editing (equal). R. O. Vargas: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing—original draft (equal); Writing—review & editing (equal). A. Mil-Martinez: Formal analysis (equal); Investigation (equal); Writing—review & editing (equal). T. N. Phillips: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing—original draft (equal); Writing—review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.