Search

We present a methodology to analyze the stationary states and stability of complex fluid flows by using hybrid, discrete, and/or continuum multi-scale simulations. Building on existing theories, our scheme extracts dynamical and equilibrium characteristics from carefully chosen time integrations of these multi-scale evolution equations. Two canonical problems are presented to demonstrate the ability and accuracy of the formalism. The first is an investigation of flow-induced transitions seen in homogeneous, hard- rod liquid crystal suspensions subjected to a linear shear flow. In the second problem, we study the phenomenon of draw resonance, a dynamical instability in an isothermal fiber-spinning process, by using a multi-scale hybrid simulation that incorporates both stochastic and continuum models.

Kinetic theory models exhibit dynamics that depend on a few low-order moments of the underlying conformational distribution function. This dependence is exhibited in a compact spectrum of eigenvalues for the Jacobian matrix associated with the dynamical system. We take advantage of this spectrum of eigenvalues through Newton-GMRES iterations to enable dynamic viscoelastic simulators (time-steppers) to obtain stationary states and perform stability/bifurcation analysis. Results are presented for three example problems: (1) the equilibrium behavior of the Doi model with the Onsager excluded volume potential, (2) pressure-driven flow of non-interacting rigid dumbbells in a planar channel, and (3) pressure-driven flow of non-interacting rigid dumbbells through a planar channel with a linear array of equally spaced cylinders.

Sharp-interface and thin-interface asymptotic analyses are presented for a generalization of the Beckermann–Karma phase-field model for solidification of a dilute binary alloy when the interface curvature is macroscopic. The ratio of diffusivities, R m ≡ D s ′ / D m ′ , in the solid and melt is arbitrary with 0 ⩽ R m ⩽ 1 . Discrepancies between this diffuse-interface model and the classical, two sided solutal model (TSM) description are quantified up to second order in the small parameter e that controls the interface thickness. We uncover extra terms in the interface species flux balance and in the Gibbs–Thomson equilibrium condition introduced by the finite width of the interface. Asymptotic results in the limit of rapid-interfacial kinetics are presented for both finite phase-field mobility and a quasi-steady state approximation for the phase-field wherein the phase-field responds passively to the concentration field. The possibility of adding additional terms to the phase-field version of the species conservation equation is explored as a means of achieving O ( e 2 ) consistency with the classical model. Our results naturally lead to a generalization of the anti-trapping solutal flux suggested by Karma [Phase-field-formulation for quantitative modeling of alloy solidification, Phys. Rev. Lett. 87(11) (2001) 115701] for the limit R m = 0 . Achieving second order accuracy for arbitrary R m requires judicious choices for the interpolating functions; these are calculated a posteriori using the functional forms of the error terms as a guide.

This study compares the predictions of several “suspension temperature” models of particle migration to laser Doppler velocimetry measurements in a concentrated suspension of noncolloidal spheres. We compare the shear rate, concentration, and suspension temperature profiles in narrow-gap Couette flow. The models predict the observed macroscopic shear rate and concentration profiles well at moderate bulk particle concentration but diverge from one another and from the data at high concentrations. In addition, the predictions of the models compare poorly with suspension temperature measurements. Most of the models greatly underpredict the magnitude of the scalar temperature, capturing instead only the magnitude of the smaller two diagonal components of the temperature tensor. Also, the models do not predict the observed variation of the suspension temperature with particle concentration. Our investigation shows that both neglect of suspension temperature anisotropy and qualitative choices of model coefficie...

A numerical method based on wavelet approximations is developed to solve the diffusion equations that arise from kinetic theory models of polymer dynamics and for coupling with continuum calculations for simulating complex flows of polymer solutions. The dynamics of a liquid crystalline polymer solution in a simple, homogeneous shear flow is computed by using the Doi model with orientational-dependent rotational diffusivity as a function of Deborah number, De, and different initial conditions. Stochastic methods were previously used for computing this flow, and aperiodic behavior was found at high De. Calculations with the wavelet-Galerkin method demonstrate a stable limit cycle in the same parameter regime. The existence of the time-periodic solution is traced to a supercritical Hopf bifurcation at De=De Hopf ; the value of De Hopf scales with ( U − U c ) α , where U is the strength of the intermolecular potential. The exponent in this correlation, α , depends on the form of the mean-field intermolecular potential and can be used to characterize the anisotropic environment experienced by the polymer molecules.

The linear stability of the flow of an Oldroyd-B fluid through a linear array of cylinders confined in a channel is analyzed by computing both the steady-state flows and the linear stability of these states by finite element analysis. The linear stability of two-dimensional base flows to three-dimensional perturbations is computed both by time integration of the linearized evolution equations for infinitesimal perturbations and by an iterative eigenvalue analysis based on the Arnoldi method. These flows are unstable to three-dimensional perturbations at a critical value of the Weissenberg number that is in very good agreement with the experimental observations by Liu [Viscoelastic Flow of Polymer Solutions around Arrays of Cylinders: Comparison of Experiment and Theory, Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, MA, 1997] for flow of a polyisobutylene Boger fluid through a linear periodic array of cylinders. The wave number of the disturbance in the neutral or spanwise direction also is in good agreement with experimental measurements. For closely spaced cylinders, the mechanism for the calculated instability is similar to the mechanism proposed by Joo and Shaqfeh [J. Fluid Mech. 262 (1994) 27] for the non-axisymmetric instability observed in viscoelastic Couette flow, where perturbations to the shearing component of the velocity gradient interact with the polymer stresses in the base flow. However, when the cylinder spacing is increased, we discover a new mechanism for instability that involves the coupling between the elongational component of the velocity gradient and the streamwise normal stress in the base state. In addition, the most unstable disturbance in the Couette geometry is over-stable (leading to time-periodic states) whereas the instability calculated for closely spaced cylinders grows monotonically with time. This apparent inconsistency is resolved by the observation that in a complex flow, the evolution of the perturbations to the base flow can be transient in a Lagrangian frame of reference, while steady in an Eulerian frame.

Recent statistical constitutive models of suspensions of neutrally buoyant, noncolloidal, solid spheres in Newtonian fluids suggest that the particles migrate in response to gradients in “suspension temperature,” defined as the average kinetic energy contained in the particle velocity fluctuations. These models have not yet been compared systematically with experimental data. In addition the “temperature'’ models assume isotropic particle velocity fluctuations, since the suspension temperature is given as a scalar. However, highly anisotropic particle velocity fluctuations have been observed experimentally, which suggests that a suspension temperature tensor is more realistic. We use laser Doppler velocimetry to measure particle velocity fluctuations arising from interparticle collisions in a concentrated suspension under nearly homogeneous shear flow in a narrow-gap concentric cylinder Couette device. We compare the relative sizes of the fluctuating velocity components and determine the variation of each...

A highly parallel time integration method is presented for calculating viscoelastic flows with the DEVSS-G/DG finite element discretization. The method is a synthesis of an operator splitting time integration method that decouples the calculation of the polymeric stress by solution of a hyperbolic constitutive equation from the evolution of the velocity and pressure fields by solution of a generalized Stokes problem. Both steps are performed in parallel. The discontinuous finite element discretization of the hyperbolic constitutive equation leads to highly-parallel element-by-element calculation of the stress at each time step. The Stokes-like problem is solved by using the BiCGStab Krylov iterative method implemented with the block complement and additive levels method (BCALM) preconditioner. The solution method is demonstrated for the calculation of two-dimensional (2D) flow of an Oldroyd-B fluid around an isolated cylinder confined between two parallel plates. These calculations use extremely fine finite elements and expose new features of the solution structure.

A wide variety of bead-spring kinetic theory models have been proposed to explain the stress growth and hysteretic behavior of dilute polymer solutions in uniaxial extension. We analyze the Kramers chain, a fine-scale model for polymer dynamics, in order to assess the validity of the coarser-grained bead-spring models in these deformations. Whereas the spring force is a simple function of the dumbbell length for the FENE spring, we find that the relationship between the ensemble-averaged end-to-end force and the extension for a Kramers chain depends on the kinematic history to which it has been subjected. In a quiescent fluid, the Kramers chain force–extension relationship is identical to the FENE force law. However, during start up of elongational flow, the ensemble-averaged end-to-end force for a given (end-to-end) length of the molecule increases with strain until steady state is reached. If the extensional flow is suddenly stopped, the Kramers chain force–extension relationship relaxes back to the FEN...

Ultimately, numerical simulation of viscoelastic flows will prove most useful if the calculations can predict the details of steady-state processing conditions as well as the linear stability and non-linear dynamics of these states. We use finite element spatial discretization coupled with a semi-implicit θ -method for time integration to explore the linear and non-linear dynamics of two, two-dimensional viscoelastic flows: plane Couette flow and pressure-driven flow past a linear, periodic array of cylinders in a channel. For the upper convected Maxwell (UCM) fluid, the linear stability analysis for the plane Couette flow can be performed in closed form and the two most dangerous, although always stable, eigenvalues and eigenfunctions are known in closed form. The eigenfunctions are non-orthogonal in the usual inner product and hence, the linear dynamics are expected to exhibit non-normal (non-exponential) behavior at intermediate times. This is demonstrated by numerical integration and by the definition of a suitable growth function based on the eigenvalues and the eigenvectors. Transient growth of the disturbances at intermediate times is predicted by the analysis for the UCM fluid and is demonstrated in linear dynamical simulations for the Oldroyd-B model. Simulations for the fully non-linear equations show the amplification of this transient growth that is caused by non-linear coupling between the non-orthogonal eigenvectors. The finite element analysis of linear stability to two-dimensional disturbances is extended to the two-dimensional flow past a linear, periodic array of cylinders in a channel, where the steady-state motion itself is known only from numerical calculations. For a single cylinder or widely separated cylinders, the flow is stable for the range of Deborah number (De) accessible in the calculations. Moreover, the dependence of the most dangerous eigenvalue on De≡ λV / R resembles its behavior in simple shear flow, as does the spatial structure of the associated eigenfunction. However, for closely spaced cylinders, an instability is predicted with the critical Deborah number De c scaling linearly with the dimensionless separation distance L between the cylinders, that is, the critical Deborah number De L c ≡ λV / L is shown to be an O (1) constant. The unstable eigenfunction appears as a family of two-dimensional vortices close to the channel wall which travel downstream. This instability is possibly caused by the interaction between a shear mode which approaches neutral stability for De ≫ 1 and the periodic modulation caused by the presence of the cylinders. Nonlinear time-dependent simulations show that this secondary flow eventually evolves into a stable limit cycle, indicative of a supercritical Hopf bifurcation from the steady base state.

Local similarity solutions are presented for the stress field of a fluid described by the Oldroyd-B viscoelastic constitutive equation near the singularity caused by the intersection of a planar free surface and a solid surface along which Navier’s slip law holds, the partial-slip/slip problem. For the case where the velocity field is given by Newtonian kinematics, the elastic stress field is predicted to have a logarithmic singularity as the point of attachment of the free surface is approached. Asymptotic analysis for the fully-coupled flow, where the stress and flow fields are determined simultaneously, results in a local form for the flow and elastic stress fields that is similar in form to that for the decoupled case. For both the coupled and decoupled flow problems, the strength of the singularity depends on the dimensionless solvent viscosity and the slip coefficient, but not upon the Deborah number. The asymptotic results for the coupled flow differ from the predictions with Newtonian kinematics i...

The local solution behavior near corners formed by the intersection of a slip surface with either a no-slip or a shear-free boundary is analyzed by finite element calculations of the two-dimensional flow of an inertialess Newtonian fluid in several model flow geometries; these flows are the flow in a tapered contraction, a sudden expansion and the extrudate swell from a planar die. Local finite element mesh refinement based on irregular, embedded elements is used to obtain extremely fine resolution of the velocity and pressure fields near the region where there is a sudden change in boundary condition. The calculations accurately reproduce the expected asymptotic behavior for a shear-free surface intersecting a no-slip boundary, where the solution is given by a self-similar form for the velocity and pressure fields. Replacing the shear-free condition with a slip condition yields a similar form for the local velocity and pressure fields and indicates that the slip boundary behaves, to leading order, as a shear-free surface. Calculations for a slip boundary intersecting a shear-free surface yield similar results, with the local behavior being given by asymptotic analysis for two shear-free surfaces intersecting to form a wedge. These results suggest that replacing the no-slip boundary condition in planar Newtonian die swell with a slip boundary condition can give rise to local behavior of velocity gradients and pressure which is more singular than the flow created with no-slip boundary conditions. This prediction is confirmed by calculations of Newtonian die swell with slip. These calculations also demonstrate that the local solution in Newtonian die swell is sensitive to the details of the numerical method.

The two‐dimensional, free‐surface flow of a Newtonian fluid exiting from a planar die is computed by finite element analysis using quasiorthogonal mesh generation and local mesh refinement with irregular, embedded elements to obtain extreme resolution of the velocity and pressure fields near the die edge, where the fluid sheet attaches to the solid boundary. Calculations for the limit of large surface tension, the stick‐slip problem, reproduce the singular behavior near the die edge expected from asymptotic analysis using a self‐similar form for the velocity field. Results for finite capillary number (Ca) predict that the meniscus separates from the die at a finite contact angle and suggest that the capillary force enters the dominant normal stress balance at the die edge through an infinite curvature, as previously suggested by Schultz and Gervasio. The size of this region with large positive curvature increases with increasing Ca, and the strength of the singularity is in good agreement with theoretical predictions for a straight meniscus attached to the die at the appropriate contact angle predicted by the simulations. The contact angle appears to be determined from matching of the inner solution structure valid near the singularity with the bulk flow, in agreement with arguments made by Ramalingam; increasing the Reynolds number decreases the contact angle, corroborating this effect. Introducing fluid slip along the surface of the die changes the structure of the singularity in the pressure and stresses, but does not alleviate the singular behavior. In fact, the calculations with slip coefficients small enough not to change the bulk solution are more difficult than calculations with the no‐slip boundary condition.

A new mixed finite element method is presented that has improved numerical stability and has similar numerical accuracy compared to the EVSS/FEM developed by Rajagopalan et al. (J. Non-Newtonian Fluid Mech. 36 (1990) 159–192). The new method, denoted EVSS-G/FEM, is based on separate bilinear interpolation of all relevant components of the velocity gradient tensor. Lagrangian bilinear approximations of this variable are formed by least-squares interpolation of the derivatives of a Lagrangian biquadratic approximation to the velocity field that is generated by solution of the momentum and continuity equations using a standard mixed formulation for velocity and pressure. Elastic components of the stress are computed by solving the constitutive equation using either the streamline upwinding (SU) or streamline upwind Petrov-Galerkin (SUPG) methods with bilinear interpolation. The motivation for the interpolation of the velocity gradients and for the stress interpolation is to force the compatibility between these variables that is needed to satisfy differential viscoelastic constitutive equations in the limit where the velocity field vanishes, as it does near solid boundaries. The enhanced numerical stability of the EVSS-G/FEM over the EVSS/FEM formulation is demonstrated for the calculation of the linear stability of planar Couette flow of a UCM fluid, where closed form results exist for the linear stability of the steady-state flow; hence, any instability detected in the numerical solution is a result of either the finite element approximation or the time integration method. The second possibility is eliminated by using a fully implicit time integration method. Calculations with the EVSS-G/FEM are stable for De >100, whereas calculations with the EVSS/FEM become numerically unstable for De >5. The EVSS-G formulation is combined with finite element discretizations that incorporate local adaptive mesh refinement for increasing the accuracy of the calculations in regions where the stress and velocity gradients change rapidly. The accuracy of the EVSS-G/FEM is demonstrated by steady-state calculations for the flow between eccentric rotating cylinders, flow through a wavy-walled tube, and for flow through a square array of cylinders; each has been used before as a test problem for the accuracy of viscoelastic flow calculations. As in previous calculations, discretization of the constitutive equation by SUPG gives superior accuracy at high values of De compared to the application of SU.

Finite‐amplitude waves propagating at constant speed down an inclined fluid layer are computed by finite element analysis of the Navier–Stokes equations written in a reference frame translating at the wave speed. The velocity and pressure fields, free‐surface shape and wave speed are computed simultaneously as functions of the Reynolds number Re and the wave number μ. The finite element results are compared with predictions of long‐wave, asymptotic theories and boundary‐layer approximations for the form and nonlinear transitions of finite‐amplitude waves that evolve from the flat film state. Comparisons between the finite element calculations and the long‐wave predictions for fixed μ and increasing Re agree well for small‐amplitude waves. However, for larger‐amplitude waves the long‐wave results diverge qualitatively from the finite element predictions; the long‐wave theories predict limit points in the solution families that do not exist in the finite element solutions. Comparisons between the finite ele...

Laser Doppler velocimetry (LDV) and video flow visualization are used to investigate the creeping motion of a highly elastic, constant-viscosity fluid flowing past a cylinder mounted centrally in a rectangular channel. A sequence of viscoelastic flow transitions are documented as the volumetric flow rate past the cylinder is increased and elastic effects in the fluid become increasingly important. Velocity profiles clearly show that elasticity has almost no effect on the kinematics upstream of the cylinder, but that the streamlines in the wake of the cylinder are gradually shifted further downstream . Finite element calculations with a nonlinear constitutive model closely reproduce the evolution of the steady two-dimensional velocity field. However, at a well defined set of flow conditions the steady planar stagnation ow in the downstream wake is experimentally observed to become unstable to a steady, three-dimensional cellular structure. The Reynolds number at the onset of the flow instability is less than 0.05 and inertia plays little role in the flow transition, LDV measurements in the wake close to the cylinder reveal large spatially periodic fluctuations of the streamwise velocity that extend along the length of the cylinder and more than five cylinder radii downstream of the cylinder. Fourier analysis shows that the characteristic spatial wavelength of these flow perturbations scales closely with the cylinder radius R . Flow visualization combined with LDV measurements also indicates that the perturbations in the velocity field are confined to the narrow region of strongly extensional flow near the downstream stagnation point. A second flow transition is observed at higher flow rates that leads to steady translation of the cellular structure along the length of the cylinder and time-dependent velocity oscillations in the wake. Measurements of the flow instability are presented for a range of cylinder sizes, and a stability diagram is constructed which shows that the onset point of the wake instability depends on both the extensional deformation of the fluid in the stagnation flow and the shearing flow between the cylinder and the channel.

Finite element analysis is used to compute the inertialess flow of a viscoelastic fluid around a sphere falling in a cylindrical tube. Calculations are reported for three differential constitutive models with constant viscosities: the Upper-Convected Maxwell model (UCM), the Oldroyd-B model (OLDB) and the dumbbell model of Chilcott and Rallison (CR). Calculations are based on the Explicitly Elliptic Momentum Equation (EEME) and the Elastic Viscous Split Stress (EVSS) methods for calculations with models without and with a Newtonian solvent contribution, respectively. The calculations converged with mesh refinement for all three models for values of De below 1.6. Calculations with the UCM and OLDB models are limited below De = 2 by loss of convergence of the Newton iteration. Mesh refinement does not seem to alleviate these limits. Calculations with the CR model and moderate values of the maximum extensibility of the dumbbell L converge to higher values of De ; however, computations for large values of L require fine meshes in the stagnation region behind the sphere.

A decoupled finite-element method (INT/FEM) is presented for calculation of two-dimensional viscoelastic flows with integral constitutive models. The momentum and continuity equations are solved by Galerkin's method with the viscoelastic stress treated as a fixed body force. The viscoelastic stress is computed by using the stream function to track fluid particles upstream, integrating a system of ordinary differential equations that govern the displacement-gradient tensor, and evaluating the integral constitutive equation by numerical quadrature. The quasi-linear upper-convected Maxwell and Oldroyd-B models, as well as the nonlinear model of Papanastasiou, Scriven and Macosko (PSM), are used in the simulations. The efficiency of the integral method is compared to that of the recently developed finite-element method (EVSS/FEM) for differential constitutive models. Convergence and accuracy of the INT/FEM are shown by calculations for flow between eccentric cylinders. The upper limit in De attainable by using the INT/FEM is comparable to values for the EVSS/FEM only for constitutive models with a shear-thinning ratio of the first normal stress difference to the shear stress and a large solvent contribution to the solution viscosity. The INT/FEM becomes the more efficient technique for simulation with this type of constitutive equation when three or more relaxation modes are included in the memory function.

The EEME/finite-element method is used to compute steady flows of fluids with constitutive behavior described by the Modified Upper Convected Maxwell Model of Apelian and a modified form of the dumbbell model of Chilcott and Rallison through abrupt, axisymmetric contractions. Both constitutive models predict constant viscosity and shear thinning first normal stress coefficient, but differ qualitatively in the behavior of the elongational viscosity. Asymptotic analysis for both models predicts that the solution has Newtonian-like spatial structure near the reentrant corner and integrable stresses and velocity gradients there. With the Newtonian-like asymptotics, the stress field can be approximated by conventional Lagrangian finite elements and computed by the streamline upwind Petrov Galerkin (SUPG) method. The finite element calculations are stable and convergent: higher values of Deborah number are reached with increasing mesh refinement. Moreover, the predicted asymptotic structure of the stress and velocity fields is recovered near the corner in the calculations. Calculations with both constitutive equations show the stretching of the Newtonian corner vortex toward the reentrant corner and its growth upstream with increasing Deborah number for 4:1 and 8:1 contraction ratios. The characteristics of the computed vortex are in semi-quantitative agreement with experiments for Boger fluids for which the flow is axisymmetric and steady.

Time-dependent simulations using the EEME finite-element method are reported for calculation of the linear stability and nonlinear dynamics of the transition to time-periodic, viscoelastic flow in axisymmetric TaylorCouette flow between parallel cylinders. The linear stability analysis is based on time integration of the linearized finite-element equations and reproduces the oscillatory linear instability recently analyzed by Larson et al. past a critical Deborah number Dec. Linear analysis presented here shows that the viscoelastic instability corresponds to a secondary flow structure composed of multiple toroidal flow cells nested radially, and that these cells travel across the gap between the cylinders. Nonlinear simulations demonstrate the existence of supercritical, i.e. for De - Dec > 0, time-periodic flow states. For only small changes in De >Dec, the flatness of the neutral stability curve with variation in the flow cell height leads to nonlinear interactions between flows that are closely spaced in De. These interactions will make the observation of simple oscillations near onset extremely difficult.

A constitutive equation for computing particle concentration and velocity fields in concentrated monomodal suspensions is proposed that consists of two parts: a Newtonian constitutive equation in which the viscosity depends on the local particle volume fraction and a diffusion equation that accounts for shear‐induced particle migration. Particle flux expressions used to obtain the diffusion equation are derived by simple scaling arguments. Predictions are made for the particle volume fraction and velocity fields for steady Couette and Poiseuille flow, and for transient start‐up of steady shear flow in a Couette apparatus. Particle concentrations for a monomodal suspension of polymethyl methacrylate spheres in a Newtonian solvent are measured by nuclear magnetic resonance (NMR) imaging in the Couette geometry for two particle sizes and volume fractions. The predictions agree remarkably well with the measurements for both transient and steady‐state experiments as well as for different particle sizes.

The finite element method based on the explicitly elliptic momentum equation (EEME) formulation for solution of viscoelastic flow problem is extended to the calculation of flows using a multimode upper convected Maxwell (UCM) model. The generalized EEME for the multimode model makes explicit the elliptic character of momentum transport in creeping flows. Mixed finite-element approximations for the velocity, pressure and stress are used to discretize the EEME and continuity equations by Galerkin's method and the constitutive equation by the streamline-upwind Petrov-Galerkin (SUPG) technique. Sample calculations are described for a two-mode model for flow between eccentric rotating cylinders. The numerical method is stable and the solutions convergent with mesh refinement. The flow structure predicted for a single-mode UCM model is largely unchanged by the addition of the second mode, although the magnitude of the stresses is affected significantly.

The explicitly elliptic momentum equation (EEME) formulation for the upper-convected Maxwell model is extended to the analysis of time-dependent flows. The formulation makes explicit the elliptic operator, which controls the spatial dependence of the velocity field, and the hyperbolic-like character of the time-dependent equations, which gives rise to the wave-like character of transients when inertia is included. Numerical solutions are computed by using finite element approximations for the spatial dependence of the velocity, stress and pressure, coupled with finite difference discretizations of time derivatives. Numerically stable and accurate solutions are demonstrated for the flows between concentric and eccentric rotating cylinders when the condition for change of type of the steady momentum equation from elliptic to hyperbolic is not satisfied. A semi-implicit algorithm is presented in which the normal stresses are computed implicitly from the corresponding components of the constitutive equation and the velocities, pressure and shear stress are computed implicitly from the EEME-continuity pair and the shear component of the constitutive equation at each time step. The coupling between these two calculations is explicit; however, calculations are not constrained to small time steps. An algorithm that decouples all stress components from velocities and pressure is shown to be severely limited to small time steps. These results are consistent with the constraints caused by the mathematical type of the governing equations.

Adding a Newtonian solvent to most differential viscoelastic constitutive equations mathematically regularizes the coupled set formed by the momentum, continuity and constitutive equations. The momentum and continuity equations form an elliptic saddle point problem for velocity and pressure, and the constitutive equation is hyperbolic in stress. Three finite element algorithms are presented that exploit this elliptic behavior by expressing the momentum equation in different differential forms. The first is based on the viscous elliptic operator that arises naturally with the introduction of a Newtonian solvent viscosity; the second is based on the explicitly elliptic momentum equation formulation developed for the upper-convected Maxwell model; and the third is based on an elastic-viscous splitting of the momentum equation. Finite element discretizations are created by using Galerkin's method for the momentum and continuity equations and the streamline-upwind Petrov-Galerkin method for the components of the constitutive equation. In the latter two methods, additional interpolants are introduced so as to maintain continuous representations of velocity derivatives across element boundaries; this is a requirement if only those higher-order velocity terms which are explicitly elliptic are integrated by parts. Calculations for flow between eccentric cylinders and through a corrugated tube demonstrate the numerical stability and accuracy of each of the formulations. The robustness of each algorithm for calculation of flows at high Deborah numbers depends on the value of the ratio β ≡ ηS/gh0, where ηs is the solvent viscosity and η0 is the viscosity of the solution. The algorithm based on the elastic-viscous splitting is the most robust across the full range 0 ≤ β ≤ 1.

We present a simple linear stability analysis of the diffusion equation for nematic polymers that delivers the equilibrium bifurcation diagram for rigid rod, excluded volume potentials. Symmetry properties of the diffusion equation yield insight into invariant subspaces of solutions and allow for interpretation of recently published simulation results.

The dynamics of concentrated solutions of rodlike molecules in nonhomogeneous shear flow are explored using a consistent numerical simulation of the Doi diffusion equation and the nonhomogeneous Onsager model of excluded-volume rod interactions. Simulations of planar, wall-driven shear flow show that out-of-plane structure instabilities occur when nematic anchoring constraints at the boundaries are removed. A new composite state with misaligned logrolling and flow-aligning domains is observed for pressure-driven flow in a planar channel. These results mark the first use of the Doi diffusion equation to show how a nonhomogeneous flow field generates sharp inter-domain interfaces analogous to those observed in rheological experiments.

The initial stage of spinodal decomposition is investigated for a rigid-rod system. Spinodal decomposition proceeds through either of two mechanisms: (1) The randomly aligned rods rotate toward a common director with no inherent length scale. (2) The rods diffuse axially and segregate into regions of common alignment with a selected length scale [script-l]. Previous studies on spinodal decomposition yielded radically different conclusions about which mechanism is dominant. A computational method is employed to analyze the growth rate and properties of the dominant fluctuation mode through an eigenvalue analysis of the linearized Doi diffusion equation in Fourier space. The linearized operator is discretized in Fourier mode and orientation space (k,theta,phi) space, and the maximum eigenvalue and corresponding eigenvector of the operator are calculated. Our analysis generalizes the results of previous studies and shows that the conflicts seen in those studies are due to differences in the diffusivities for rotational motion, motion perpendicular to the rod axis, and motion along the rod axis. For a given system density, a plot of the dominant perturbation wave number k(max) as a function of the diffusivity ratios shows two separate regions corresponding to mechanisms (1) and (2). High rotational diffusivity corresponds to mechanism (1), whereas high axial diffusivity corresponds to mechanism (2). The transition between the two mechanisms depends on the ratio of diffusivities and on system density. Also, the dominant wave number increases with increasing density indicating that a deeper quench into the spinodal regime leads to a smaller characteristic length scale.

Isotropic-Nematic and Nematic-Nematic transitions from a homogeneous base state of a suspension of high aspect ratio, rod-like magnetic particles are studied for both Maier-Saupe and the Onsager excluded volume potentials. A combination of classical linear stability and asymptotic analyses provides insight into possible nematic states emanating from both the isotropic and nematic non-polarized equilibrium states. Local analytical results close to critical points in conjunction with global numerical results (Bhandar, 2002) yields a unified picture of the bifurcation diagram and provides a convenient base state to study effects of external orienting fields., Comment: 3 Figures

A crystallization model appropriate for application in continuum modeling of complex processes is presented. As an extension to the previously developed Schneider equations [ W. Schneider, A. Köppel, and J. Berger, "Non-isothermal crystallization of polymers," Int. Polym. Proc. 2, 151 (1988) ], the model presented here allows one to account for the growth of crystals of various shapes and to distinguish between one-, two-, and three-dimensional growth, e.g., between rod-like, plate-like, and sphere-like growth. It is explained how a priori knowledge of the shape and growth processes is to be built into the model in a compact form and how experimental data can be used in conjunction with the dynamic model to determine its growth parameters. The model is capable of treating transient processing conditions and permits their straightforward implementation. By using thermodynamic methods, the intimate relation between the crystal shape and the driving forces for phase change is highlighted. All these capabilities and the versatility of the method are made possible by the consistent use of four structural variables to describe the crystal shape and number density, irrespective of the growth dimensionality.

The steady-state and time-dependent flow transitions observed in a well-characterized viscoelastic fluid flowing through an abrupt axisymmetric contraction are characterized in terms of the Deborah number and contraction ratio by laser-Doppler velocimetry and flow visualization measurements. A sequence of flow transitions are identified that lead to time-periodic, quasi-periodic and aperiodic dynamics near the lip of the contraction and to the formation of an elastic vortex at the lip entrance. This lip vortex increases in intensity and expands outwards into the upstream tube as the Deborah number is increased, until a further flow instability leads to unsteady oscillations of the large elastic vortex. The values of the critical Deborah number for the onset of each of these transitions depends on the contraction ratio β, defined as the ratio of the radii of the large and small tubes. Time-dependent, three-dimensional flow near the contraction lip is observed only for contraction ratios 2 [les ] β [les ] 5, and the flow remains steady for higher contraction ratios. Rounding the corner of the 4:1 abrupt contraction leads to increased values of Deborah number for the onset of these flow transitions, but does not change the general structure of the transitions.

The steady-state, similarity solutions of the flow of an upper-convected Maxwell fluid through a tube with a porous wall are constructed by asymptotic and numerical analyses as functions of the direction of flow through the tube, the amount of elasticity in the fluid, as measured by the Deborah number De, and the degree of fluid slip along the tube wall. Fluid slip is assumed to be proportional to the local shear stress and is measured by a slip parameter β that ranges between no-slip (β = 1) and perfect slip (β = 0). The most interesting results are for fluid injection into the tube. For β = 1, the family of flows emanating from the Newtonian limit (De = 0) has a limit point where it turns back to lower values of De. These solutions become asymptotic to De = 0) and develop an O(De) boundary layer near the tube wall with singularly high stresses matched to homogeneous elongational flow in the core. This solution structure persists for all nonzero values of the slip parameter. For β ≠ 1, a family of exact solutions is found with extensional kinematics, but nonzero shear stress convected into the tube through the wall. These flows differ for low De from the Newtonian asymptote only by the absence of the boundary layer at the tube wall. Finite difference calculations evolve smoothly between the Newtonian-like and extensional solutions because of approximation error due to under-resolution of the boundary layer. The radial gradient of the axial normal stress of the extensional flow is infinite at the centerline of the tube for De > 1; this singularity causes failure of the finite difference approximations for these Deborah numbers unless the variables are rescaled to take the asymptotic behavior into account.

Approximate analytical solutions of the local temperature variations, due to viscous heating, in channel flows of non-Newtonian fluids with small Nahme-Griffith numbers were obtained by the WKB-J method. The results are general and easy to use so that extensive numerical computations are not required. Comparisons between the approximate analytical solution and complete numerical solution showed that the relative differences were very small.

A solid sphere falling through a Bingham plastic moves in a small envelope of fluid with shape that depends on the yield stress. A finite-element/Newton method is presented for solving the free-boundary problem composed of the velocity and pressure fields and the yield surfaces for creeping flow. Besides the outer surface, solid occurs as caps at the front and back of the sphere because of the stagnation points in the flow. The accuracy of solutions is ascertained by mesh refinement and by calculation of the integrals corresponding to the maximum and minimum variational principles for the problem. Large differences from the Newtonian values in the flow pattern around the sphere and in the drag coefficient are predicted, depending on the dimensionless value of the critical yield stressYgbelow which the material acts as a solid. The computed flow fields differ appreciably from Stokes’ solution. The sphere will fall only whenYgis below 0.143 For yield stresses near this value, a plastic boundary layer forms next to the sphere. Boundary-layer scalings give the correct forms of the dependence of the drag coefficient and mass-transfer coefficient on yield stress for values near the critical one. The Stokes limit of zero yield stress is singular in the sense that for any small value ofYgthere is a region of the flow away from the sphere where the plastic portion of the viscosity is at least as important as the Newtonian part. Calculations For the approach of the flow field to the Stokes result are in good agreement with the scalings derived from the matched asymptotic expansion valid in this limit.

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

For a dilute solution of macromolecules idealized as elastic dumbbells with Brownian motion, expressions are obtained giving the stress tensor in terms of rate‐of‐strain tensors for an arbitrary, steady, homogeneous flow. The rheological equations of state are fit into a retarded motion expansion. For dumbbells with arbitrary, nonlinear, elastic connectors, the result is given up through terms of second‐order. Stress tensor equations for two specific models, Hookean dumbbells and finitely extendible non‐linear elastic (FENE) dumbbells, are given up through terms of third order.

Results are presented for the flow of an upper-convected Maxwell (UCM) model through an axisymmetric corrugated tube with sinusoidally varying cross-section. The calculations are based on the explicitly elliptic momentum equation (EEME) which makes explicit the elliptic character of the momentum equation when inertia is absent. Mixed finite element approximations for velocity, pressure, and stress coupled with Galerkin's method for the EEME/continuity pair and streamline upwind Petrov-Galerkin (SUPG) method for the hyperbolic constitutive equations are shown to converge with mesh refinement. The maximum value of the Weissenberg number that can be computed for a particular discretization is limited only by inability to resolve the spatial structure that forms in both stress and velocity fields. Calculations are reported for the Weissenberg number up to 15 and are in excellent agreement with the spectral/finite-difference calculations reported recently by Pilitsis and Beris.

Under the influence of a uniform and parallel magnetic field, a ferromagnetic fiber suspended in a Newtonian fluid rotates to align with the field direction. This study examines the field-induced rotation process for an individual non-Brownian axisymmetric ellipsoid suspended in a stagnant Newtonian fluid. Theoretical predictions are derived by a perturbation analysis for the limiting case where the strength of the applied magnetic field far exceeds the saturation magnetization of the ellipsoid. Numerical calculations are performed for the more general problem of an ellipsoid with known isotropic, non-hysteretic magnetic properties, using nickel and a stainless steel as examples. The analysis encompasses materials with field-induced, nonlinear magnetic properties, distinguishing these results from the simpler cases where the particle magnetization is either independent of, or linearly dependent on, the strength of the applied external field. In this study, predictions indicate that when the ellipsoid is magnetically saturated, the particle rotation is governed by the magnetoviscous time constant,τ MV = ηs/e0 M s 2 . It is found that the rotation rate depends strongly on the aspect ratio,a/b, of the ellipsoid, but only weakly on the dimensionless magnetization,M s/H 0.

For an arbitrary, time‐dependent, homogeneous flow, an expression is obtained for the stress tensor describing the behavior of a dilute solution of macromolecules idealized as elastic dumbbells with Brownian motion. The result is given in terms of the relaxation modulus from linear viscoelasticity. The general relaxation modulus found for an arbitrary, nonlinear, elastic dumbbell is specialized to that for a finitely extendible nonlinear elastic (FENE) dumbbell. Good agreement is found when the result is compared with Warner's numerical calculations for the FENE dumbbell.

Laser Doppler velocimetry (LDV) measurements of the axial and radial velocity components are reported for the axisymmetric contraction flow of a viscoelastic solution of polyisobutylene dissolved in tetradecane. This fluid exhibits shear-thinning of both the viscosity and first normal stress difference; the viscosity is well described by the Carreau-Yasuda model. Shear-rate-dependent Reynolds and Deborah numbers are used to described the flow conditions. Measurements at low flow rates for contraction ratios of 2:1, 4:1, and 8:1 indicate that the flow near the contraction plane scales with the small tube radius. A viscoelastic corner vortex is detected at moderate flow rates, yet vortex growth is not observed with increased flow rate, because inertial forces press the vortex into the outer corner. At high flow rates, where both the Reynolds number and Deborah number are large, the flow near the contraction plane is divided into an accelerating core and an outer region where the flow retains its upstream profile. The magnitude of the local rate-of-strain is calculated from the axial and radial velocity profiles. In the core, the magnitude of >50s −1 and its value is dominated by the axial gradients ∂ v /∂ z and ∂ v r /∂ z . Large extension rates along the centerline suggest that extensional thinning occurs in this core region. Finite element calculations for a Carreau-Yasuda fluid with the same viscosity as the PIB/C14 solution and under the same flow conditions do not predict the two-regions of the flow.

Calculations of viscoelastic flows using the upper-convected Maxwell (UCM) model in geometries which include sharp corners or moving and free liquid/fluid contact lines are known to be non-convergent with mesh refinement. A modified upper-convected Maxwell (MUCM) model is proposed which partially alleviates this difficulty. The MUCM model is derivable from network theory and allows the fluid relaxation time to decrease at increasing values of the trace of the stress tensor. The MUCM model yields stress fields that reduce to the a symptotic expressions for a Newtonian fluid near singularities at non-deformable boundaries. Calculations using a Galerkin finite element method are presented for the planar stick-slip problem of the flow between two no-slip surfaces joined to two shear-free surfaces. Results for fine meshes show the correct asymptotic behavior near the singularity for the MUCM model and converge to much higher values of the Deborah number than for the UCM model. However, the results for the MUCM model are still constrained by numerical instabilities related to approximating the stress behavior near the singularity.

The two-dimensional flows between two slightly eccentric cylinders with the inner one rotating are calculated by finite-element methods for five viscoelastic constitutive equations each derived as a limit of the Giesekus model. Comparisons with exact results for Newtonian, second-order, and corotational Maxwell-like fluids set the accuracy of the calculations as a function of eccentricity and Deborah number (De). Computer-implemented perturbation methods are used to demonstrate bifurcation and turning points in De for an upper-convected Maxwell fluid. The locations of the limit points are moderately stable to extensive mesh refinement and, therefore, seem to be an intrinsic property of this constitutive equation. Similar solution pathology is demonstrated for the three-constant Oldroyd-B model, but no limiting value of De is found for calculations with the Leonov-like version of the Giesekus fluid.

Finite element calculations of two-dimensional flows of viscoelastic fluids in a journal bearing geometry reported in an earlier paper (J. Non-Newt. Fluid Mech. 16 (1984) 141-172) are extended to higher eccentricity (ρ = 0.4); at this higher eccentricity flow separation occurs in the wide part of the gap for a Newtonian fluid. Calculations for the second-order fluid (SOF), upper-convected Maxwell (UCM), and the Giesekus models are continued in increasing Deborah number for each model until either a limit point is reached or oscillations in the solution make the numerical accuracy too poor to warrant proceeding. No steady solutions to the UCM model were found beyond a limit point De c , as was the case for results at low eccentricities. The value of De c was moderately stabel to mesh refinement. A limit point also terminated the calculations with a SOF model, in contradiction to the theorems for uniqueness and existence for this model. The critical value of De increased drastically with increasing refinement of the mesh, as expected for solution pathology caused by approximation error. Calculations for the Giesekus fluid with the mobility parameter α ≠ O showed no limit points, but failed when irregular oscillations destroyed the quality of the solution. The behavior of the recirculation region of the flow and the load on the inner cylinder were very sensitive to the value of α used in the Giesekus model. The recirculation disappeared at low values of De except when the mobility parameter α was so small that the viscosity was almost constant over the range of shear rates in the calculations. The recirculation persisted over the entire range of accessible De for the UCM fluid, the limit of α = O of the Giesekus model. The behavior of the recirculation is coupled directly to the viscosity by calculations with an inelastic fluid with the same viscosity predicted by the Giesekus model.

Numerical calculations of complex, two-dimensional flows of viscoelastic fluids fail when the elastic contribution to the stresses, measured by a Deborah number, becomes comparable to the viscous contribution. Reasons for the limit on Deborah number are explored by a sequence of finite element calculations with contravariant convected Maxwell and second-order fluid models. The possibilities of bifurcation or loss of a steady, two-dimensional flow field are ruled out by employing continuation methods for flow through a planar contraction and by existence and uniqueness proofs for a second-order fluid in a driven cavity. Error in the finite element approximation to steep gradients of stress causes the breakdown of all calculations. This is most clearly seen in the calculations for a second-order fluid, because the exact flow field for any Deborah number is known.

Finite element calculations for viscoelastic flows are reported that use a restructured form of the equation of motion that makes explicit the elliptic character of this equation. We call this restructured equation the Explicity Elliptic Momentum Equation, and its use is illustrated for flow of an upper convected Maxwell (UCM) model between eccentric and concentric rotating cylinders and also for a modified upper convected Maxwell (MUCM) model in the stick-slip problem. Sets of mixed-order approximations for velocity, stress, and a modified pressure are used to test the algorithm in both problems. Both sets of calculations are shown to converge with mesh refinement and are limited at high values of Deborah number by the formation of elastic boundary layers that are identified in the momentum equation by the growth of low-order derivative terms that involve the local velocity gradient and divergence of stress. Similar convergence properties are observed for bilinear and biquadratic Lagrangian approximations to the stress components. However, calculations with the more accurate basis for stress converge to higher values of De and are sensitive to the weighted residual method used for the constitutive equation, particularly for the eccentric cylinder problem. Streamline-upwind Petroy-Galerkin (SUPG) and artificial diffusivity (AD) formulations of the constitutive equation are tested for solution of both problems by calculations of the stress fields with fixed kinematics and by solution of the coupled problem. The SUPG method improves the performance of the calculations with the biquadratic basis set for the eccentric cylinder problem. For the UCM model, adding artificial diffusion to the constitutive equation in the stick-slip problem changes the dominant balance for the stress field near the singularity, making it appear as an integrable stress approximation for fixed mesh. For the MUCM model the Newtonian-like behavior of the stress near this point is unaffected by the AD method and calculations converge to moderate De .

A new, automated, two-color laser Doppler velocimetry apparatus is described which has been constructed especially for mapping flow fields for viscoelastic fluids. This system is capable of measuring both slow, secondary flows and time-dependent motions that are characteristic of viscoelastic flow in complex geometries. Application of this apparatus to flow between eccentric cylinders and through an axisymmetric sudden contraction for a Newtonian fluid and for polyisobutylene (PIB) and polyacrylamide (PAC) solutions is described. In both geometries excellent agreement is found between calculations and experiments for Newtonian fluids, as expected. The Newtonian recirculation region which is present for flow between cylinders with large eccentricity disappears for the PAC solution and shrinks with increasing De for the PIB solution. These changes are attributed to the substantial shear thinning of the PAC solution and the slight shear thinning of the PIB mixture under the conditions studied. New flow transitions were discovered in the sudden contraction flow of the polyisobutylene solution. At a critical value of Deborah number De 1 ≈ 0.80 the flow changed from a steady, two-dimensional motion, nearly identical with the Newtonian result, to a time-periodic flow with a tangential velocity component that fluctuated about zero. Multiple, time-periodic motions and flow hysteresis were found for De 1 De De 2 ≈ 1.20. At De 2 the flow suddenly reverts to a two-dimensional, time-independent motion. A vortex emanating from the lip of the contraction was observed for De > De 2 .

A dilute solution of flexible macromolecules is modeled as a dilute suspension of dumbbells with slightly nonlinear elastic connectors (nearly‐Hookean dumbbells) and subject to Browian motion. The molecular variables are eliminated directly from the kinetic theory by means of a perturbation scheme to give an expression for the stress tensor which is good for flows that do not deform the dumbbells drastically. The stress tensor is expressed in terms of several structure tensors, each of which is given by an auxiliary differential equation. This Hand form of the rheological equation of state can also be expressed in an Oldroyd form or as a memory integral expansion. Calculations of rheological properties for the model show substantial improvements over Hookean dumbbell predictions. In an appendix the distribution function for the nearly‐Hookean dumbbell is derived, and it can be used to obtain the same rheological equation of state by traditional methods.

The steady-state, two-dimensional creeping flow of an Upper-Convected Maxwell fluid between two eccentric cylinders, with the inner one rotating, is computed using a spectral/finite-element method (SFEM). The SFEM is designed to alleviate the numerical oscillations caused by excessive dispersion error in previous finite-element calculations and to resolve the stress boundary-layers that exist for high elasticity, as measured by the Deborah number De . Calculations for cylinders with low eccentricity (ϵ = 0.1) converged to oscillation-free solutions for De ≈ 90, extending the domain of convergence over traditional finite-element methods by a factor of thirty. The results are confirmed by extensive refinement of the discretization. At high De , steep radial boundary layers form in the stress, which match closely with those predicted by asymptotic analysis. Calculations at higher eccentricity require extreme refinement of the discretization to resolve the variations in the stress field in both the radial and azimuthal directions associated with the existence of the recirculation region. Results for ϵ = 0.4 show that the recirculation region present for the Newtonian fluid ( De = 0) shrinks and then grows with increasing De . Calculations for ϵ = 0.4 are terminated by a limit point near DeL ≈ 7.24 for the finest discretization used. The Fourier series approximations are not convergent for this mesh, so the limit point must be considered to be an artifact of the discretization.

For a suspension of rigid macromolecules with Brownian motion, expressions are obtained for the components of the stress tensor for any time‐dependent shearing flow. The results are expressed up to terms of the fourth order as a series of memory integrals. The first two terms in the series enable one to calculate the explicit expressions for the kernel functions for second‐order viscoelasticity. It is established that the Oldroyd six‐constant model does not give kernel functions which are of the same form as those given by the molecular theory of rigid macromolecules. The general time‐dependent elongational flow results are also given through terms of the third order.

The use of the retarded motion expansion as a constitutive equation has long been known to be valid only for small Deborah number flow. We demonstrate that ‘solution’ which cannot be obtained as a regular perturbation about the Newtonian solution is physically meaningless. This is discussed in connection with a recent paper in this journal.