Multiscale simulations of viscoelastic fluids in complex geometries using a finitely extensible nonlinear elastic transient network model

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


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) scheme 3 or using the logarithm of the conformation tensor 4 are often used to extend the range of the Weissenberg numbers for which converged approximations are obtained.
Phan-Thien and Dou 5 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-Thien 6 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 Phillips 7 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-Thien 6 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 Miranda 9 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 Nithiarasu 11 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 Phillips 13 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 Ouyang 17 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 model 19 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.2][23] Rinc on 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. 25In 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 on 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 shearthinning.These effects agree with the predictions of Vaccaro and Marrucci. 23In 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 solutions 28 and a coupling of the model to a rheological constitutive equation of state within an irreversible thermodynamics framework. 29errer et al. 19 coupled the FENE model with the transient network of Rinc on 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 nondimensionalization 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
The mathematical statement of conservation of mass for an incompressible fluid is where v is the velocity field.The mathematical statement of conservation of momentum, also known as the Cauchy equation of motion, is where q is the density of the fluid, p is the pressure, and r is the extrastress tensor.The extra-stress tensor r is decomposed into solvent and polymeric contributions, r s and r p , respectively, where r ¼ r s þ r p .The solvent contribution is given by where g s is the solvent viscosity, and d is the rate-of-strain tensor defined by Using the decomposition of the stress just described, the momentum equation can be written in the form The polymeric contribution to the stress tensor is traditionally calculated by means of a constitutive equation.However, there are kinetic theory models that are able to provide an accurate representation of polymeric fluids but which do not possess an equivalent closed form constitutive equation.It is for this reason that in this paper the polymeric contribution to the stress is evaluated by means of a microscopic description of polymer dynamics.

B. Microscopic equations
Consider an elastic dumbbell comprising two identical beads and a connecting elastic spring immersed in a Newtonian solvent.The FENE dumbbell model is a particular case of the general elastic dumbbell model (see Fig. 1) in which the spring has finite extensibility.Suppose that the configuration of the dumbbell is represented by the end-to-end vector Q as shown in Fig. 1.The vector Q provides information about the stretch and orientation of the dumbbell.The equation of motion for the beads in the dumbbell can be derived by considering the following forces: the forces due to Brownian motion, the elastic spring force, and the viscous drag force.This is the starting point for the derivation of the Fokker-Planck equation, for the details of which the reader is referred to € Ottinger, 30 for example.The Fokker-Planck equation for the FENE dumbbell model is where w is the configuration probability density function (PDF) and wðQ; x; tÞdQ represents the probability of finding a dumbbell with configuration in the range Q to Q þ dQ at ðx; tÞ, where jðtÞ is the velocity gradient, 1 is the friction coefficient, k b is Boltzmann's constant, T is the absolute temperature, and f represents the spring force which for FENE dumbbells is given by where H is the spring constant and Q 0 is the maximum extension of each spring.The high-dimensional nature of the Fokker-Planck equation means that it is expensive to solve computationally.However, one can use the equivalence between the Fokker-Planck equation 5 and the stochastic differential equation, as the basis for deriving computationally tractable numerical schemes.In this equation, UðtÞ is a multi-dimensional Wiener process [see € Ottinger, 30 € Ottinger et al., 31 Lozinski and Chauviere, 32 for details].It is possible to establish a relation between the material constants at the microscopic and macroscopic levels: k ¼ 1=4H and g p ¼ nk b Tk, where k, g p , and n are the relaxation time, polymeric viscosity, and the number of dumbbells per unit of volume, respectively.Every dumbbell is characterized by k and Q 0 .Once the configuration vector field is defined, Kramers expression is used to determine r p : 25,31

C. Transient network model
In the FENE dumbbell model, a single dumbbell represents the behavior of a large number of chain molecules.Many of these chain molecules tend to get entangled or disentangled as a result of flow or due to energy conditions, which cause a local change in properties such as the maximum chain length, drag coefficient, elastic constant, and relaxation time, all of which vary with time.In the dilute regime, it is possible to consider a low number of interactions and low or negligible changes for every dumbbell.As the concentration increases, those changes are not negligible.In order to model these changes in behavior, it is possible to define some characteristic states termed "microstates." 24hese are shown in Fig. 1.These microstates can be defined by the number of chain segments and the number of nodes needed to produce it. 24n this framework, it is possible to consider different ways to determine what constitutes a microstate.For this reason, we only consider states that are more likely to appear and ones that provide information about the system.Every microstate, x i ; i ¼ 0; 1; 2; 3; 4, that defines a chain arrangement is identified as shown in Table I, where l i represents the fraction of maximum extension length for every microstate. 19,24,33As mentioned before, the formation and destruction of microstates depends on the available energy of the system to create or destroy a node.From a kinetic scheme analogous to that of a chemical reaction, the creation or destruction of each microstate can be determined as a function of free energy and viscous dissipation as follows: where k A i ; k B i represent the formation and destruction constants, respectively, for every microstate, and E 0 is the energy needed to create a node.The numerical values assigned to the kinetic constants are based on the assumption that the ratio between them is equal to the energy ratio of each kinetic equation, which means that ; which leads to the following simplified system of ordinary differential equations: x 0 0 1 1 1 x 2 2 7 3 3 7 where C i is the concentration of microstructure x i .Note that it is important that the total length of all microstates, L p , remains constant in time i.e., In the present paper, every microstate x i in a region of space defines a microstructure, which behaves as a dumbbell with its own properties: partial viscosity fraction i , relaxation time s i , and maximum extension length b i , as shown in Table II.Under these assumptions, the stress produced by each microstructure is obtained by the Kramers expression which can be expressed as where g i ¼ C i k b Tk i is the partial viscosity, and k i is the relaxation time associated with the microstate x i .The polymeric contribution to the stress is then obtained as

III. DIMENSIONLESS EQUATIONS
Let L and U denote the characteristic length and velocity scales.Define the following dimensionless variables, where g 0 denotes the zero shear viscosity defined by g 0 ¼ g s þ g p in the initial state.The dimensionless macroscopic equations are then where Re ¼ qUL=g 0 is the Reynolds number and b ¼ g s =g 0 is the viscosity ratio.The end-to-end vector Q is made dimensionless using ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k b T=H p so that the dimensionless stochastic differential equation 7 is where k H ¼ 1=ð4HÞ is the relaxation time and the dimensionless force is The polymeric stress is evaluated using the Kramers expression where Wi ¼ k H U=L is the Weissenberg number, and s i ¼ k i =k H and i ¼ g i =g p are the relaxation time and viscosity of the ith microstructure, respectively.The parameter a i in ( 26) is given by where d is the spatial dimension of the problem, and b i ¼ bl 2 i is the dimensionless maximum extension length of the ith microstructure x i , where b ¼ HQ 2 0 =k b T. The dimensionless form of the stochastic differential equation for the transient network model is where W represents the dimensionless Wiener process.The quantity s i Wi ¼ k i U=L may be viewed as the Weissenberg number associated with x i .
For the concentration, the dimensionless equations are in term of the kinetic constants for construction and destruction A and B, which are related to k A and k B in Eqs. ( 13)-( 17) by as follows: In order to simplify the notation, the symbol Ã is omitted from here onward.

IV. NUMERICAL SCHEME
The system is solved using a decoupled approach in which the velocity and pressure (macroscopic step) are solved separately from the stress (microscopic step) at each time step.Within each time step, Eq. ( 28) is solved using a predictor and corrector method to determine the stress while the conservation equations are solved using an implicit method.The spectral element method is used for the spatial discretization.The physical domain X is partitioned into The space of all polynomials of degree N on X k is denoted as P N ðX k Þ, and we define Every spectral element is mapped onto a parent element D ¼ ½À1; 1 Â ½À1; 1, using a transfinite mapping technique that associates each point ðn; gÞ 2 D with a point ðxðn; gÞ; yðn; gÞÞ 2 X k . 34agrangian interpolants of degree N in both spatial directions are used to approximate the dependent variables on D, based on the Gauss-Lobatto-Legendre (GLL) points creating a GLL grid inside the spectral elements. 13The general algorithm is summarized in the flow chart shown in Fig. 2.

A. Macroscopic equations
In the macroscopic stage of the time evolution of the equations, the polymeric contribution to the extra-stress is known.The velocity and pressure at the new time level are then obtained by solving the following semi-discrete problem: This is an implicit treatment of the conservation equations where ṽnþ1 corresponds to the solution of the following pure convection problem at time t ¼ t nþ1 : This is a first-order operator integration factor splitting (OIFS) technique. 35The pure convection problem is solved using a fourth-order Runge-Kutta method.The solution spaces at each time level are P ¼ L 2 0 ðXÞ for pressure, T ¼ ½L 2 ðXÞ 4 s for stress where s indicates the space of symmetric tensors and a subspace V of ½H 1 ðXÞ 2 for velocity. 13The subspace V contains vectors that satisfy the given Dirichlet boundary conditions for velocity.
The weak formulation of Eqs. ( 35) and ( 36) is discretized using appropriate approximation spaces for velocity, pressure, and stress (see Vargas et al. 12 for details).To guarantee that the corresponding discrete problem is well-posed, the LBB condition must be satisfied by the discrete approximation spaces.Maday et al. 35 have shown that for spectral elements the LBB condition is satisfied when the velocity approximation space is the polynomial space P N ðXÞ and the pressure approximation space is P NÀ2 ðXÞ.A Gauss-Lobatto quadrature rule is used to integrate the velocities, while a Gauss quadrature rule is used to integrate the pressure.The stress is approximated by polynomials in the space P N ðXÞ with the difference being that the stress components are allowed to be discontinuous across the element boundaries. 13The discrete version of the weak formulation of the problem is where D N and C N denote the discrete divergence operators acting on the velocity and stress, A N , D T N , and B N are the discrete Laplace, gradient, and velocity mass matrix operators, respectively, and g N is a vector that contains boundary data.Using the discrete continuity equation, the velocity is eliminated resulting in the pressure equation: where H N is a discrete Helmholtz-like operator, defined as The operator U N ¼ D N H À1 N D T N is the Uzawa operator.The pressure equation may be written as where b N stands for the right-hand side of Eq. ( 40).To solve Eq. ( 42), a nested solver is required for the pressure solution, because it is necessary to invert both the Uzawa operator and the Helmholtz operator.Schur's complement method has been used to reduce the size of the Helmholtz operator.An efficient preconditioner for the Uzawa operator, based on the finite element mass and stiffness matrices on the local finite elements, M FE k and E FE k , respectively, was proposed by Gerritsma and Phillips 13 and is given by where R k is the restriction operator, mapping a global vector to a vector of length equal to the number of GLL nodes of the spectral element, for full details of the method, see Ref. 12.

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,37The 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 Q m ðx; tÞ, with m ¼ 1; …; N f .The Eulerian version of Eq. ( 7) is obtained by adding a convective term to the right-hand side (see Lozinski et al. 38 ), Since the Wiener process second moment is OðtÞ and using the properties of UðtÞ described previously, Eq. (44) becomes At each point ðX; tÞ in space and time, a set of conformation vectors is generated that experience the same deformation history, with different stochastic processes.The polymeric stress r p can be determined from the conformation vectors as follows: As the fluid is initially stationary, then r p ¼ 0 at t ¼ 0. Every initial conformation field Q m ðX; t 0 Þ must satisfy lim for an initial state.For sufficiently large N f , the central limit theorem states that the statistical error is proportional to 1= ffiffiffiffiffi ffi N f p , which means the numerical convergence is order O 1= ffiffiffiffiffi ffi N f p .Equation ( 45) can be generalized to be applicable for each microstate: in which case the polymeric contribution to the extra-stress tensor at the new time level t nþ1 is then computed using where the subscript i denotes the microstructure and j denotes the trajectory so that i ¼ 0; …; 4; and j ¼ 1; …; N f .The temporal discretization of the evolution equation for the configuration fields (48) is performed using a second-order predictor-corrector scheme that is based on a forward Euler predictor and trapezoidal corrector: The right-hand side of Eq. ( 51) gives the orientation of the conformation vector Q nþ1 ij , with an arbitrary length.In order to obtain the length of Q nþ1 ij , Eq. ( 51) is expressed as or where L stands for the right-hand side of Eq. ( 51) and x ¼ jQ nþ1 ij j.The solutions of Eq. ( 53) give possible values for the length of the conformation vector Q nþ1 ij .€ Ottinger 30 has shown that there exists a unique solution of this cubic equation in the interval ½0; ffiffiffiffi b i p so there is precisely one physically admissible solution.This scheme circumvents difficulties that may arise using other temporal schemes in which unphysical extensions can be generated during the iterative process.The semidiscrete stochastic differential equations ( 50) and (51) are discretized in space using the discontinuous Galerkin (DG) method. 12The convection term is integrated by parts twice.In the forward step, the condition Q ¼ Q in is applied weakly along inflow portions, c in , of X k .The inflow portion of the boundary is characterized by the condition v Á n < 0, and Q in is the value of Q in the neighboring upwind spectral element or the prescribed condition at the inflow to the domain, and n is the unit outward normal to the boundary of X k .In the backward step, the new boundary terms are left unchanged.Therefore, in the weak formulation of Eqs. ( 50) or (51) over each spectral element X k , terms of the form ðv n Á rQ i ; SÞ X k are replaced by 12 The spatial discretization of the configuration fields is performed using discontinuous approximations in P N ðXÞ.

V. PROBLEM DESCRIPTION
The complex flow past a cylinder placed symmetrically in a channel is considered.The aspect ratio is defined by K ¼ R=H, where R is the cylinder radius and H is the half-width of the channel.In this paper, K ¼ 0:5.The computational domain is extended a distance 25 units upstream and 25 units downstream of the cylinder, thus the assumption of fully developed flow conditions at the inlet and outlet is valid.Figure 3 shows a close-up of the spectral element mesh used in this work.The number of spectral elements is K ¼ 20, the order of the spectral approximation is N ¼ 6, and the time step is Dt ¼ 2 Â 10 À3 .The dimensionless drag coefficient is used to verify the numerical accuracy of the solution of the flow through a cylinder by comparing with results published in the literature.The expression for the dimensionless drag on the cylinder is

A. Initial conditions
Initially, the fluid is assumed to be at rest, so that vð0Þ ¼ 0. Additionally, we require an equilibrium distribution for the configuration fields at time t ¼ 0. At equilibrium, the polymeric contribution to the extra-stress tensor should be zero.The equilibrium distribution, w eq , must satisfy ð where the integral is over configuration space.For the FENE model, the initial configuration fields are generated using the equilibrium distribution function which is a Gaussian distribution with zero mean and a covariance matrix which is b i =ðb i þ dÞ times the unit matrix.

B. Boundary conditions
A fully developed velocity profile with a mean velocity of unity is imposed at entry and exit: No-slip and no-penetration conditions are imposed on the surface of the cylinder.Furthermore, symmetry conditions are imposed along the axis of symmetry.

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 on et al. 24 and improved by Ferrer et al. 19 and G omez-L opez 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, Wi 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 N f reduces the fluctuations in the evolution of the drag and that there is not a significant improvement beyond N f ¼ 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 x 0 ¼ 1 and x i ¼ 0 for i ¼ 1; 2; 3; 4, (ii) entanglement, characterized by x 4 ¼ 1 and x i ¼ 0 for i ¼ 0; 1; 2; 3, and (iii) aleatory, where all the microstates are present (x i 6 ¼ 0 for i ¼ 0; 1; 2; 3; 4).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 and N f ¼ 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.

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 on 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 x 4 (see Table I).In Fig. 6, the concentrations of microstates x 2 and x 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 t % 3, and this is not shown in the figure.

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: 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 Wi 8.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 Wi % 1 for b ¼ 0:59. 7,10,40In 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. 8Afonso 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.

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 b ¼ 0:1 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,36Figure 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
For the flow around a sphere, the steady-state Stokes problem is solved.The drag is calculated by integrating the total stress over the surface of the sphere, F, given by the dimensionless drag factor F Ã , which is used for comparison of numerical solutions, defined as the ratio between the drag experienced by a sphere and the drag that the same sphere would experience in a Newtonian fluid without boundaries, given by Using the same mesh (20 spectral elements and N ¼ 6) as for the cylinder case, a dimensionless drag of F Ã ¼ 5:947 64 was determined, which compared to the value of F Ã ¼ 5:9474 reported in the literature 12,41 is a very good approximation.Figure 16 shows the evolution of drag force around a sphere by modifying the kinetic rate constants, the maximum extension length and elasticity, respectively, for  Re ¼ 0.01.For all cases, the same trends analyzed for the flow around the cylinder are present.

VII. CONCLUSIONS
This paper presents the first complex flow simulations of viscoelastic fluids using a FENE transient network model 24 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.

8 V
FIG.3.Sample mesh used to simulate the flow around a cylinder or a sphere with K ¼ 20 and N ¼ 6.

FIG. 4 .
FIG. 4. Evolution of the drag force on the cylinder (a) as a function of polynomial order N with N f ¼ 2000 configuration fields and (b) as a function of N f for N ¼ 6, with Re ¼ 0.01, Wi ¼ 1, b ¼ 0:1, b ¼ 50.

FIG. 6 .
FIG. 6. Influence of the kinetic rate constant A and B on the evolution of microstate concentrations for different initial networks around a cylinder with: Re ¼ 0.01, Wi ¼ 1, b ¼ 0:1, b ¼ 50, N ¼ 6, and N f ¼ 2000.

FIG. 7 .FIG. 9 .
FIG. 7. Effect of the kinetic rate constant A and B on the shear and normal stresses along the central line for different initial network configurations around a cylinder with: Wi ¼ 1, Re ¼ 0.01, b ¼ 0:1, b ¼ 50, N ¼ 6, and N f ¼ 2000.