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.

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.

FIG. 1.

Complex flow around a confined object. Cylinder radius R placed symmetrically in a two-dimensional channel of half width Hi.

FIG. 1.

Complex flow around a confined object. Cylinder radius R placed symmetrically in a two-dimensional channel of half width Hi.

Close modal

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.

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.

The mathematical statement of conservation of mass for an incompressible fluid is
(1)
where v is the velocity field. The mathematical statement of conservation of momentum, also known as the Cauchy equation of motion, is
(2)
where ρ is the density of the fluid, p is the pressure, and σ is the extra-stress tensor. The extra-stress tensor σ is decomposed into solvent and polymeric contributions, σ s and σ p, respectively, where σ = σ s + σ p. The solvent contribution is given by
where η s is the solvent viscosity, and d is the rate-of-strain tensor defined by
(3)
Using the decomposition of the stress just described, the momentum equation can be written in the form
(4)
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.
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 Öttinger,30 for example. The Fokker-Planck equation for the FENE dumbbell model is
(5)
where ψ is the configuration probability density function (PDF) and ψ ( Q , x , t ) d Q represents the probability of finding a dumbbell with configuration in the range Q to Q + d Q at ( x , t ), where κ ( t ) is the velocity gradient, ς 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
(6)
where H is the spring constant and Q0 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 equation5 and the stochastic differential equation,
(7)
as the basis for deriving computationally tractable numerical schemes. In this equation, Φ ( t ) is a multi-dimensional Wiener process [see Öttinger,30 Öttinger 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: λ = ς / 4 H and η p = n k b T λ, where λ, η p, and n are the relaxation time, polymeric viscosity, and the number of dumbbells per unit of volume, respectively. Every dumbbell is characterized by λ and Q0. Once the configuration vector field is defined, Kramers expression is used to determine σ p:25,31
(8)
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.”24 These 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.24 In 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, ω 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,33 As 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:
(9)
(10)
(11)
(12)
where k i A , k i B 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
and
which leads to the following simplified system of ordinary differential equations:
(13)
(14)
(15)
(16)
(17)
where Ci is the concentration of microstructure ωi. Note that it is important that the total length of all microstates, L p, remains constant in time i.e.,
(18)
TABLE I.

Microstate properties.

Microstate Number of nodes Number of segments Number of chains l i
ω0 
ω1  1 2 
ω2  3 7 
ω3  1 3 
ω4  12  1 3 
Microstate Number of nodes Number of segments Number of chains l i
ω0 
ω1  1 2 
ω2  3 7 
ω3  1 3 
ω4  12  1 3 
In the present paper, every microstate ωi in a region of space defines a microstructure, which behaves as a dumbbell with its own properties: partial viscosity fraction νi, relaxation time τi, and maximum extension length bi, as shown in Table II. Under these assumptions, the stress produced by each microstructure is obtained by the Kramers expression
(19)
which can be expressed as
(20)
where η i = C i k b T λ i is the partial viscosity, and λi is the relaxation time associated with the microstate ωi. The polymeric contribution to the stress is then obtained as
(21)
TABLE II.

Microstructure coefficients.

Microstructure νi τi li
ω0  1 4 
ω1  1 2  1 2  1 2 
ω2  3 4  1 3  3 7 
ω3  3 4  1 3  1 3 
ω4  1 4  1 3 
Microstructure νi τi li
ω0  1 4 
ω1  1 2  1 2  1 2 
ω2  3 4  1 3  3 7 
ω3  3 4  1 3  1 3 
ω4  1 4  1 3 
Let L and U denote the characteristic length and velocity scales. Define the following dimensionless variables,
(22)
where η0 denotes the zero shear viscosity defined by η 0 = η s + η p in the initial state. The dimensionless macroscopic equations are then
(23)
(24)
where Re = ρ U L / η 0 is the Reynolds number and β = η s / η 0 is the viscosity ratio. The end-to-end vector Q is made dimensionless using k b T / H so that the dimensionless stochastic differential equation7 is
(25)
where λ H = ς / ( 4 H ) is the relaxation time and the dimensionless force is
The polymeric stress is evaluated using the Kramers expression
(26)
where W i = λ H U / L is the Weissenberg number, and τ i = λ i / λ H and ν i = η i / η p are the relaxation time and viscosity of the ith microstructure, respectively. The parameter αi in (26) is given by
(27)
where d is the spatial dimension of the problem, and b i = b l i 2 is the dimensionless maximum extension length of the ith microstructure ωi, where b = H Q 0 2 / k b T. The dimensionless form of the stochastic differential equation for the transient network model is
(28)
where W represents the dimensionless Wiener process. The quantity τ i Wi = λ i U / L may be viewed as the Weissenberg number associated with ω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 kA and kB in Eqs. (13)–(17) by
as follows:
(29)
(30)
(31)
(32)
(33)
In order to simplify the notation, the symbol * is omitted from here onward.
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 Ω is partitioned into K nonoverlapping spectral elements Ωk, 1 k K, such that k = 1 K Ω k = Ω. The space of all polynomials of degree N on Ωk is denoted as P N ( Ω k ), and we define
(34)
Every spectral element is mapped onto a parent element D = [ 1 , 1 ] × [ 1 , 1 ], using a transfinite mapping technique that associates each point ( ξ , η ) D with a point ( x ( ξ , η ) , y ( ξ , η ) ) Ω k.34 Lagrangian 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.13 The general algorithm is summarized in the flow chart shown in Fig. 2.
FIG. 2.

Flow chart of the micro-macro numerical algorithm.

FIG. 2.

Flow chart of the micro-macro numerical algorithm.

Close modal
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:
(35)
(36)
This is an implicit treatment of the conservation equations where v ̃ n + 1 corresponds to the solution of the following pure convection problem at time t = t n + 1:
(37)
This is a first-order operator integration factor splitting (OIFS) technique.35 The pure convection problem is solved using a fourth-order Runge–Kutta method. The solution spaces at each time level are P = L 0 2 ( Ω ) for pressure, T = [ L 2 ( Ω ) ] s 4 for stress where s indicates the space of symmetric tensors and a subspace V of [ H 1 ( Ω ) ] 2 for velocity.13 The 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 ( Ω ) and the pressure approximation space is P N 2 ( Ω ). 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 ( Ω ) with the difference being that the stress components are allowed to be discontinuous across the element boundaries.13 The discrete version of the weak formulation of the problem is
(38)
(39)
where DN and CN denote the discrete divergence operators acting on the velocity and stress, AN, D N T, and BN are the discrete Laplace, gradient, and velocity mass matrix operators, respectively, and gN is a vector that contains boundary data. Using the discrete continuity equation, the velocity is eliminated resulting in the pressure equation:
(40)
where HN is a discrete Helmholtz-like operator, defined as
(41)
The operator U N = D N H N 1 D N T is the Uzawa operator. The pressure equation may be written as
(42)
where bN 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 k F E and E k F E, respectively, was proposed by Gerritsma and Phillips13 and is given by
(43)
where Rk 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.

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 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),
(44)
Since the Wiener process second moment is O ( t ) and using the properties of Φ ( t ) described previously, Eq. (44) becomes
(45)
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 σ p can be determined from the conformation vectors as follows:
(46)
As the fluid is initially stationary, then σ p = 0 at t = 0. Every initial conformation field Q m ( X , t 0 ) must satisfy
(47)
for an initial state. For sufficiently large Nf, the central limit theorem states that the statistical error is proportional to 1 / N f, which means the numerical convergence is order O ( 1 / N f ). Equation (45) can be generalized to be applicable for each microstate:
(48)
in which case the polymeric contribution to the extra-stress tensor at the new time level t n + 1 is then computed using
(49)
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:
(50)
(51)
The right-hand side of Eq. (51) gives the orientation of the conformation vector Q i j n + 1, with an arbitrary length. In order to obtain the length of Q i j n + 1, Eq. (51) is expressed as
(52)
or
(53)
where L stands for the right-hand side of Eq. (51) and x = | Q i j n + 1 |. The solutions of Eq. (53) give possible values for the length of the conformation vector Q i j n + 1. Öttinger30 has shown that there exists a unique solution of this cubic equation in the interval [ 0 , b i ] 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.12 The convection term is integrated by parts twice. In the forward step, the condition Q = Q i n is applied weakly along inflow portions, γin, of Ωk. The inflow portion of the boundary is characterized by the condition v · n < 0, and Q i n 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 Ω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 Ωk, terms of the form ( v n · Q i , S ) Ω k are replaced by
(54)
for i = 1 , , N f,.12 The spatial discretization of the configuration fields is performed using discontinuous approximations in P N ( Ω ).
The complex flow past a cylinder placed symmetrically in a channel is considered. The aspect ratio is defined by Λ = R / H, where R is the cylinder radius and H is the half-width of the channel. In this paper, Λ = 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 Δ t = 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
(55)
FIG. 3.

Sample mesh used to simulate the flow around a cylinder or a sphere with K = 20 and N = 6.

FIG. 3.

Sample mesh used to simulate the flow around a cylinder or a sphere with K = 20 and N = 6.

Close modal
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, ψeq, must satisfy
(56)
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.
A fully developed velocity profile with a mean velocity of unity is imposed at entry and exit:
(57)
No-slip and no-penetration conditions are imposed on the surface of the cylinder. Furthermore, symmetry conditions are imposed along the axis of symmetry.

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.

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.

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, W i = 1, β = 0.1, 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 β = 0.1 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.

FIG. 4.

Evolution of the drag force on the cylinder (a) as a function of polynomial order N with Nf = 2000 configuration fields and (b) as a function of Nf for N = 6, with Re = 0.01, Wi = 1, β = 0.1, b = 50.

FIG. 4.

Evolution of the drag force on the cylinder (a) as a function of polynomial order N with Nf = 2000 configuration fields and (b) as a function of Nf for N = 6, with Re = 0.01, Wi = 1, β = 0.1, b = 50.

Close modal

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 ω 0 = 1 and ω i = 0 for i = 1 , 2 , 3 , 4, (ii) entanglement, characterized by ω 4 = 1 and ω i = 0 for i = 0 , 1 , 2 , 3, and (iii) aleatory, where all the microstates are present ( ω i 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 Re = 0.01, β = 0.1, 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.

FIG. 5.

Dependence of the evolution of the drag force on the cylinder on polynomial order N for the disentanglement configuration with Re = 0.01, β = 0.1, b = 50, A = 1, B = 1, and Nf = 2000.

FIG. 5.

Dependence of the evolution of the drag force on the cylinder on polynomial order N for the disentanglement configuration with Re = 0.01, β = 0.1, b = 50, A = 1, B = 1, and Nf = 2000.

Close modal

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.

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, β = 0.1, b = 50, N = 6, and N f = 2000.

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, β = 0.1, b = 50, N = 6, and N f = 2000.

Close modal

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.

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, β = 0.1, b = 50, N = 6, and N f = 2000.

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, β = 0.1, b = 50, N = 6, and N f = 2000.

Close modal
FIG. 8.

Effect of the kinetic rate constants A and B on the drag force evolution for different initial network configurations around a cylinder with: Wi = 1, Re = 0.01, β = 0.1, b = 50, N = 6, and Nf = 2000.

FIG. 8.

Effect of the kinetic rate constants A and B on the drag force evolution for different initial network configurations around a cylinder with: Wi = 1, Re = 0.01, β = 0.1, b = 50, N = 6, and Nf = 2000.

Close modal

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, β = 0.1, 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 W i 8. For Wi > 8, the steady-state drag value has not quite been reached by t = 5.

FIG. 9.

Effect of Wi on the evolution of microstate concentrations for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

FIG. 9.

Effect of Wi on the evolution of microstate concentrations for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

Close modal
FIG. 10.

Effect of Wi on the shear and normal stresses along the centerline for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

FIG. 10.

Effect of Wi on the shear and normal stresses along the centerline for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

Close modal
FIG. 11.

Effect of Wi on the evolution of the drag force for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

FIG. 11.

Effect of Wi on the evolution of the drag force for different initial network configurations around a cylinder with: Re = 0.01, β = 0.1, A = 1, B = 1, b = 50, N = 6, and N f = 2000.

Close modal
FIG. 12.

Effect of Wi on the evolution of the drag force for flow past a cylinder for the disentanglement configuration for Re = 0.01, β = 0.1, A = 1, B = 1, N = 6, and Nf = 2000 (log-log plot).

FIG. 12.

Effect of Wi on the evolution of the drag force for flow past a cylinder for the disentanglement configuration for Re = 0.01, β = 0.1, A = 1, B = 1, N = 6, and Nf = 2000 (log-log plot).

Close modal

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 W i 1 for β = 0.59.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 β = 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,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.

FIG. 13.

Effect of b on the evolution of microstate concentrations with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6, and N f = 2000.

FIG. 13.

Effect of b on the evolution of microstate concentrations with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6, and N f = 2000.

Close modal
FIG. 14.

Effect of b on the shear and normal stresses along the central line with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6, and N f = 2000.

FIG. 14.

Effect of b on the shear and normal stresses along the central line with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6, and N f = 2000.

Close modal
FIG. 15.

Effect of b on the drag force evolution with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6 and N f = 2000.

FIG. 15.

Effect of b on the drag force evolution with different entanglement scenarios around a cylinder for: Re = 0.01, Wi = 1, β = 0.1, A = 1, B = 1, N = 6 and N f = 2000.

Close modal
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
(58)
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
(59)
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 literature12,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.
FIG. 16.

Dependence of the evolution of the drag on a sphere on (a) kinetic rate constants, (b) b, (c) Wi, with Re = 0.01 , β = 0.1, b = 50, A = 1, B = 1, Nf = 2000, N = 6, and disentanglement initial conditions.

FIG. 16.

Dependence of the evolution of the drag on a sphere on (a) kinetic rate constants, (b) b, (c) Wi, with Re = 0.01 , β = 0.1, b = 50, A = 1, B = 1, Nf = 2000, N = 6, and disentanglement initial conditions.

Close modal

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
R. G.
Owens
and
T. N.
Phillips
,
Computational Rheology
(
Imperial College Press
,
London
,
2002
).
2.
A. N.
Brooks
and
T. J.
Hughes
, “
Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations
,”
Comput. Methods Appl. Mech. Eng.
32
,
199
259
(
1982
).
3.
R.
Guénette
and
M.
Fortin
, “
A new mixed finite element method for computing viscoelastic flows
,”
J. Non-Newtonian Fluid Mech.
60
,
27
52
(
1995
).
4.
R.
Fattal
and
R.
Kupferman
, “
Constitutive laws for the matrix-logarithm of the conformation tensor
,”
J. Non-Newtonian Fluid Mech.
123
,
281
285
(
2004
).
5.
N.
Phan-Thien
and
H.-S.
Dou
, “
Viscoelastic flow past a cylinder: Drag coefficient
,”
Comput. Methods Appl. Mech. Eng.
180
,
243
266
(
1999
).
6.
H.-S.
Dou
and
N.
Phan-Thien
, “
Viscoelastic flow past a confined cylinder: Instability and velocity inflection
,”
Chem. Eng. Sci.
62
,
3909
3929
(
2007
).
7.
S.
Claus
and
T. N.
Phillips
, “
Viscoelastic flow around a confined cylinder using spectral/hp element methods
,”
J. Non-Newtonian Fluid Mech.
200
,
131
146
(
2013
).
8.
M. A.
Hulsen
,
R.
Fattal
, and
R.
Kupferman
, “
Flow of viscoelastic fluids past a cylinder at high Weissenberg number: Stabilized simulations using matrix logarithms
,”
J. Non-Newtonian Fluid Mech.
127
,
27
39
(
2005
).
9.
P. J.
Oliveira
and
A. I.
Miranda
, “
A numerical study of steady and unsteady viscoelastic flow past bounded cylinders
,”
J. Non-Newtonian Fluid Mech.
127
,
51
66
(
2005
).
10.
A.
Afonso
,
M. A.
Alves
,
F. T.
Pinho
, and
P. J.
Oliveira
, “
Uniform flow of viscoelastic fluids past a confined falling cylinder
,”
Rheol. Acta
47
,
325
348
(
2008
).
11.
C.-B.
Liu
and
P.
Nithiarasu
, “
The characteristic-based split (CBS) scheme for viscoelastic flow past a circular cylinder
,”
Numer. Methods Fluids
57
,
157
176
(
2008
).
12.
R. O.
Vargas
,
O.
Manero
, and
T. N.
Phillips
, “
Viscoelastic flow past confined objects using a micro–macro approach
,”
Rheol. Acta
48
,
373
395
(
2009
).
13.
M. I.
Gerritsma
and
T. N.
Phillips
, “
Spectral element methods for axisymmetric Stokes problems
,”
J. Comput. Phys.
164
,
81
103
(
2000
).
14.
V.
Ribeiro
,
P.
Coelho
,
F. T.
Pinho
, and
M. A.
Alves
, “
Viscoelastic fluid flow past a confined cylinder: Three-dimensional effects and stability
,”
Chem. Eng. Sci.
111
,
364
380
(
2014
).
15.
I. E.
Garduño
,
H. R.
Tamaddon-Jahromi
, and
M. F.
Webster
, “
Flow past a sphere: Predicting enhanced drag with shear-thinning fluids, dissipative and constant shear-viscosity models
,”
J. Non-Newtonian Fluid Mech.
244
,
25
41
(
2017
).
16.
S.
Peng
,
Y.-L.
Xiong
,
X.-Y.
Xu
, and
P.
Yu
, “
Numerical study of unsteady viscoelastic flow past two side-by-side circular cylinders
,”
Phys. Fluids
32
,
083106
(
2020
).
17.
C.
Ruan
and
J.
Ouyang
, “
Numerical simulation of the non-isothermal viscoelastic flow past a confined cylinder
,”
Chin. J. Chem. Eng.
18
,
177
184
(
2010
).
18.
M.
Norouzi
,
S. R.
Varedi
, and
M.
Zamani
, “
Wake instability of viscoelastic flows past an unconfined inclined square cylinder
,”
Phys. Fluids
28
,
023101
(
2016
).
19.
V. H.
Ferrer
,
A.
Gómez
,
J. A.
Ortega
,
O.
Manero
,
E.
Rincón
,
F.
López-Serrano
, and
R. O.
Vargas
, “
Modeling of complex fluids using micro-macro approach with transient network dynamics
,”
Rheol. Acta
56
,
445
459
(
2017
).
20.
H.
Sim
,
K.
Ahn
, and
S.
Lee
, “
Large amplitude oscillatory shear behavior of complex fluids investigated by a network model: A guideline for classification
,”
J. Non-Newtonian Fluid Mech.
112
,
237
250
(
2003
).
21.
M. S.
Green
and
A. V.
Tobolsky
, “
A new approach to the theory of relaxing polymeric media
,”
J. Chem. Phys.
14
,
80
92
(
1946
).
22.
A. S.
Lodge
, “
A network theory of flow birefringence and stress in concentrated polymer solutions
,”
Trans. Faraday Soc.
52
,
120
130
(
1956
).
23.
A.
Vaccaro
and
G.
Marrucci
, “
A model for the nonlinear rheology of associating polymers
,”
J. Non-Newtonian Fluid Mech.
92
,
261
273
(
2000
).
24.
E.
Rincón
,
A. E.
Chávez
,
R.
Herrera
, and
O.
Manero
, “
Rheological modelling of complex fluids: A transient network model with microstates
,”
J. Non-Newtonian Fluid Mech.
131
,
64
77
(
2005
).
25.
R. B.
Bird
,
F. C.
Curtiss
,
R. C.
Amstrong
, and
O.
Hassager
,
Dynamics of Polymeric Liquids. Volume 2: Kinetic Theory
, 2nd ed. (
Wiley-Interscience
,
1987
).
26.
H. R.
Warner
, “
Kinetic theory and rheology of dilute suspensions of finitely extendible dumbbells
,”
Ind. Eng. Chem. Fundam.
11
,
379
387
(
1972
).
27.
J.-F.
Berret
,
D. C.
Roux
, and
G.
Porte
, “
Isotropic-to-nematic transition in wormlike micelles under shear
,”
J. Phys. II France
4
,
1261
1279
(
1994
).
28.
J. P.
García-Sandoval
,
A.
Martín del Campo
,
F.
Bautista
,
O.
Manero
, and
J. E.
Puig
, “
Nonhomogeneous flow of micellar solutions: A kinetic–network theory approach
,”
AIChE J.
64
,
2277
2292
(
2018
).
29.
J. P.
García-Sandoval
,
E.
Hernández
,
F.
Bautista
,
J. E.
Puig
, and
O.
Manero
, “
A simple kinetic model for complex rheological fluids based on irreversible thermodynamics
,”
AIChE J.
66
,
e16766
(
2020
).
30.
H. C.
Öttinger
,
Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing Simulation Algorithms
(
Springer
,
1996
).
31.
H. C.
Öttinger
,
B. H.
Van Den Brule
, and
M. A.
Hulsen
, “
Brownian configuration fields and variance reduced CONNFFESSIT
,”
J. Non-Newtonian Fluid Mech.
70
,
255
(
1997
).
32.
A.
Lozinski
and
C.
Chauvière
, “
A fast solver for Fokker–Planck equation applied to viscoelastic flows calculations: 2D FENE model
,”
J. Comput. Phys.
189
,
607
625
(
2003
).
33.
O.
Manero
,
J. E.
Puig
,
F.
Bautista
, and
J. P.
Garcia-Sandoval
, “
Nonlinear viscoelasticity of complex fluids: A kinetic network model
,”
Rheol. Acta
54
,
53
67
(
2015
).
34.
W. J.
Gordon
and
C. A.
Hall
, “
Construction of curvilinear co-ordinate systems and applications to mesh generation
,”
Numer. Methods Eng.
7
,
461
477
(
1973
).
35.
Y.
Maday
,
A. T.
Patera
, and
E. M.
Rønquist
, “
An operator-integration-factor splitting method for time-dependent problems: Application to incompressible fluid flow
,”
J. Sci. Comput.
5
,
263
292
(
1990
).
36.
M. A.
Hulsen
,
A. P. G.
Van Heel
, and
B. H. A. A.
Van Den Brule
, “
Simulation of viscoelastic flows using Brownian configuration fields
,”
J. Non-Newtonian Fluid Mech.
70
,
79
(
1997
).
37.
T. N.
Phillips
and
K. D.
Smith
, “
A spectral element approach to the simulation of viscoelastic flows using Brownian configuration fields
,”
J. Non-Newtonian Fluid Mech.
138
,
98
110
(
2006
).
38.
A.
Lozinski
,
R. G.
Owens
, and
T. N.
Phillips
, “
The Langevin and Fokker–Planck equations in polymer rheology
,” in
Numerical Methods for Non-Newtonian Fluids
, Handbook of Numerical Analysis Vol.
16
, edited by
R.
Glowinski
and
J.
Xu
(
Elsevier
,
2011
), pp.
211
303
.
39.
A.
Gómez-López
,
V. H.
Ferrer
,
E.
Rincón
,
J. P.
Aguayo
,
A. E.
Chávez
, and
R. O.
Vargas
, “
Large-amplitude oscillatory shear flow simulation for a FENE fluid
,”
Rheol. Acta
58
,
241
260
(
2019
).
40.
W.
Doherty
,
T. N.
Phillips
, and
Z.
Xie
, “
A stabilised finite element framework for viscoelastic multiphase flows using a conservative level-set method
,”
J. Comput. Phys.
477
,
111936
(
2023
).
41.
R. G.
Owens
and
T. N.
Phillips
, “
Steady viscoelastic flow past a sphere using spectral elements
,”
Int. J. Numer. Methods Eng.
39
,
1517
1534
(
1996
).
Published open access through an agreement withJISC Collections