Search

The dynamics of a solution of rodlike liquid-crystalline molecules are simulated for the related problems of isotropic–nematic spinodal decomposition and the coarsening of misaligned nematic grains. The Doi diffusion equation for the rod distribution function is coupled with the full Onsager intermolecular potential, discretized by the finite-element method, and integrated forward in time by using a parallel, semi-implicit scheme. The Onsager potential models rod interaction on the scale of a single rod length in order to resolve accurately defects and interfaces in structure. Simulation results show the effects of rotational and translational diffusivity ratios on the mechanisms for alignment and phase separation in spinodal decomposition. As rotational motion is restricted, individual grains become more aligned prior to coalescence events. When rods are restricted to diffusive motion along their axis, the spinodal decomposition process is arrested, and the system will reach a pseudo-steady state featuring misaligned nematic grains. These results mark the first dynamic computation of the Doi diffusion equation for spinodal decomposition in nonhomogeneous rigid-rod systems. Nematic coarsening simulations show the effects of misalignment between neighboring ordered domains on the coarsening time and director field around structured interfaces. Results show that the coarsening time is dependent not only on misalignment between grain directors, but also on the tilt angle of the directors into the interface.

The rheological phase diagrams for solutions of rigid rod molecules are computed for planar Couette and Poiseuille flow as a function of Deborah number, initial orientation, and wall separation; this analysis extends the seminal work of Larson and Ottinger [R.G. Larson, H.C. Ottinger, Effect of molecular elasticity on out-of-plane orientations in shearing flows of liquid-crystalline polymers, Macromolecules 24 (1991) 6270–6282] to nonhomogeneous flows. The Doi diffusion equation is solved by using a finite-element discretization of the rod distribution function and Onsager interaction potential. Simulations of planar Couette flow show the familiar logrolling–tumbling–wagging–flow-aligning cascade as Deborah number increases, and the critical Deborah numbers associated with transitions between states vary with wall separation. Defects are caused solely by wall interactions rather than artificial anchoring conditions. Simulations of planar Poiseuille flow show that the system tends toward either logrolling states or composite flow-aligning/logrolling attractors depending on the Deborah number, wall separation, and initial orientation.

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.

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.

A model and simulation method are developed for two-dimensional non-isothermal melt spinning of a visco elastic melt. The visco elastic stress is evaluated from a non-isothermal Giesekus constitutive equation developed by application of the pseudo-time method to the isothermal form of the model [J. Non-Newt. Fluid Mech. (2001)]. The crystallization kinetics is described with the model proposed by Nakamura et al. [J. Appl. Polym. Sci. 17 (1973) 1031], whereas the crystallization rate, which is a function of both temperature and molecular orientation, is evaluated according to the equation proposed by Ziabicki [Fundamentals of Fiber Formation, Wiley, New York, 1976]. The set of non-linear governing equations is solved by using the DEVSS-G/SUPG finite element method. Melt spinning is simulated for two different polymers: amorphous polystyrene and fast-crystallizing Nylon-6,6. The analysis demonstrates that although the kinematics in the thread-line are approximately one-dimensional, the radially non-uniform thermal history, caused by the leading order variation of the temperature gradient ∂T/∂r, gives rise to radially non-uniform visco elastic stresses. This stress gradient results in radially non-uniform molecular orientation and a strong radial variation in crystallinity for Nylon-6,6. The radially non-uniform stress profiles obtained from the simulations are in good agreement with experimental results for melt spinning of polystyrene. Simulations of Nylon-6,6 show that the thermally-induced crystallization depends strongly on the choice of the Avrami index n, and a sharp increase in crystallinity due to stress-induced crystallization is predicted only when the molecules are highly oriented in the drawing direction at high drawing speeds. The significant influences of visco elasticity, air drag, and operating conditions on non-isothermal melt spinning dynamics also are predicted.

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.

We present numerical methods for analysis of the stability of two-dimensional steady viscoelastic flows to small-amplitude, two-dimensional and three-dimensional disturbances based on finite element calculations of the steady base flow and the perturbation. Direct time integration of the linearized equations of motion and iterative calculation of the most dangerous components of the eigenspectrum are tested. Finite element discretizations based on the DEVSS-G finite element discretization with Newton’s method used to compute steady-state solutions. Two different time integration schemes are tested for computing the time evolution of general, random disturbances: a θ -method operator-splitting scheme and a fourth-order Runge–Kutta method. For both time integrators, time stepping is decoupled into a solution of a modified Stokes problem and an evaluation of the time-dependent constitutive equation. The overall efficiency of both methods is extremely high, as is the potential for implementation on parallel computers. An algorithm also is presented for calculating eigenvalues with the largest real part that combines time integration of the linearized equations with a Krylov subspace method to accelerate the calculation of the eigenvalues. Although this method does not dramatically reduce the computational cost over use of time integration alone, it does provide a more complete analysis of the eigenspectrum. For both direct time integration and the hybrid time integration/Krylov calculation, the stability results for cylindrical Couette flow show quantitative agreement with the eigenvalues calculated by using other methods of analysis [M. Avgousti, A.N. Beris, J. Non-Newtonian Fluid Mech. 50 (1993) 225–251]. Contrary to the results in our previous paper [R. Sureshkumar, M.D. Smith, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 82 (1999) 57–104], we find that the flow of an Oldroyd-B fluid through a closely-spaced cylinder array is stable to two-dimensional perturbations. However, allowing the perturbations to be three-dimensional and considering an isolated cylinder does not alter the conclusions of our earlier study [R. Sureshkumar, M.D. Smith, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 82 (1999) 57–104] of the two-dimensional stability of widely-spaced arrays of cylinders; the flow around an isolated cylinder is computed to be stable for all values of the Weissenberg number obtainable with these calculations, We≤0.75.

Accurate and robust finite element methods for computing flows with differential constitutive equations require approximation methods that numerically preserve the ellipticity of the saddle point problem formed by the momentum and continuity equations and give numerically stable and accurate solutions to the hyperbolic constitutive equation. We present a new finite element formulation based on the synthesis of three ideas: the discrete adaptive splitting method for preserving the ellipticity of the momentum/continuity pair (the DAVSS formulation), independent interpolation of the components of the velocity gradient tensor (DAVSS-G), and application of the discontinuous Galerkin (DG) method for solving the constitutive equation. We call the method DAVSS-G/DG. The DAVSS-G/DG method is compared with several other methods for flow past a cylinder in a channel with the Oldroyd-B and Giesekus constitutive models. Results using the Streamline Upwind Petrov–Galerkin method (SUPG) show that introducing the adaptive splitting increases considerably the range of Deborah number (De) for convergence of the calculations over the well established EVSS-G formulation. When both formulations converge, the DAVSS-G and DEVSS-G methods give comparable results. Introducing the DG method for solution of the constitutive equation extends further the region of convergence without sacrificing accuracy. Calculations with the Oldroyd-B model are only limited by approximation of the almost singular gradients of the axial normal stress that develop near the rear stagnation point on the cylinder. These gradients are reduced in calculations with the Giesekus model. Calculations using the Giesekus model with the DAVSS-G/DG method can be continued to extremely large De and converge with mesh refinement.

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.

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

Experimental measurements of the shear stress and normal stresses in the viscometric flow of a well-characterized polyisobutylene/polybutene solution between a cone-and-plate and between parallel plates show the onset of an “anti-thixotropic” flow transition, which results in a time-dependent, apparent shear-thickening of the viscosity and first normal stress difference. Measurements of the critical conditions for the onset of the flow transition, made by systematically varying the plate separation in the parallel-plate geometry, demonstrate that this rotational instability is a function of the angular velocity in the flow and initially grows exponentially in time, in agreement with theoretical predictions for an Oldroyd-B model. However, experiments show that the disturbance which causes this flow transition is nonaxisymmetric, over-stable in time and subcritical in γ. Spectral analysis also shows that the nonlinear flow which ultimately develops is temporally aperiodic. These features of the flow transition do not agree with the existing theoretical analysis. The experimentally constructed flow-stability diagram shows that, for a given fluid, the instability occurs at a critical Deborah number which is insensitive to the specific geometrical parameters, temperature, and shear-rate of the flow. This interpretation of the elastic instability clarifies recent measurements of the flow stability for the Ml viscoelastic test fluid.

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.

It is generally assumed in this work, as elsewhere, that the mass diffusion coefficient has an Arrhenius temperature dependence. This allows a combustion wave propagated by thermal conduction to develop similar to more conventional combustion systems. It is found that for realistic systems the problem must involve a three length scale analysis. Similar to the classic theory on premixed, gas-phase flames, the largest scale is identified with thermal conduction and a “reaction zone” scale proportional to the inverse of the mass diffusion activation energy is associated with a small length over which significant mass diffusion is possible. Unlike classic flame theory a still smaller scale is identified with the length over which mass diffusion takes place (i.e., the size of a typical domain of alloyable constituent). Flame speeds are derived for three different geometric configurations of a binary system using a singular perturbation analysis. First a system where the constituents are arrayed in alt...

A dynamic model of liquid propellant combustion is introduced which generalizes the classical models of Landau [ Acta Physicochim . 19:77–85 (1944)] and Levich [ Dokl. Akad. Nauk SSSR 109:975–978 (1956)]. These earlier studies assumed that the liquid/gas interface propagated normal to itself with constant velocity, and showed that steady, planar burning was intrinsically unstable unless gravitational effects, coupled with either the effects due to surface tension or (liquid) viscosity, were sufficiently strong. In contrast, the present model assumes a more general pyrolysis burning rate law, and preliminary results reported in Part I [ Twenty-Second Symposium (International) on Combustion , 1988] indicate that the resulting coupling of the temperature field can render steady, planar burning stable with respect to hydrodynamic (Landau) instability even in the absence of gravity and surface tension. At the same time, a new type of reactive/diffusive pulsating instability appears for sufficiently large values of a temperature sensitivity parameter. In Part II, a complete linear stability analysis of both hydrodynamic and reactive/diffusive instabilities in the limit of zero viscosity is presented, allowing for a pressure, as well as a temperature, dependence in the burning rate law, and nonzero gravity and surface tension. The results indicate still more types of hydrodynamic and reactive/diffusive instability phenomena, as well as a stable parameter region for steady, planar deflagration.

Time-dependent processes of solid propellant deflagration are modeled and a numerical simulation of ignition and extinction is undertaken. Since the analytic theory of ignition is well developed we rely on it to verify the ability of the numerical algorithm to track stiff, time-dependent features of an ignition event. The verification of our numerical method provides a degree of confidence for the subsequent study of deradiation, which is less well understood. An activation energy asymptotic theory for deradiation is developed, which predicts that an influx of radiation on the order of the steady heat release of the exothermic reaction will, at near term, extinguish the propellant. Beyond this, conduction processes inherited from the thermal profile at the instant of deradiation may, or may not, cause resumption of pyrolysis, and, presumably, the deflagration as a whole. The analytic theory ignores the subsequent impact of a possible gas-phase flame that is detached many conductive-diffusive lengths away from the gas solid surface. The assumption is based on the observation that for most physical systems the gas density—and hence its thermal capacitance—is many times smaller than that of the solid and thus lacks the potential to alter the surface temperature after deradiation. The simulation is used to verify more complicated and more physical situations than those accessible via the theory, including deradiation with a detached flame and in the presence of an instability.

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.

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.

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.

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.

A constitutive model for foams developed previously is extended here to study the influence of polydispersity on small deformations and the coupled effects of viscous and interfacial forces present in the foam films on both small and large deformations. A formalism for describing cell motion at large strains is also presented. To investigate the effects of polydispersity, the cell deformation and stress—strain behavior for a foam with a bimodal cell size distribution are calculated. The yield stress, critical (yield) strain and the stress—strain relation are found to be independent of the size distribution about a constant mean cell size and are coincident with that of foam consisting of monodisperse, regular hexagonal cells. The introduction of liquid viscosity into the model produces strong deviations from previous results in both the stress field and in the cell deformation for modified capillary numbers larger than 0.01. Further, these viscous effects are dependent on initial cell orientation. Although the magnitude of the viscous stress is not large, it affects the total stress by changing the cell structure. For smaller Ca′, however, viscosity has a minimal effect on the foam rheology. For large deformations, initial cell orientation and viscosity strongly affect the time periodicity of the system. For certain orientations, the cell structure and stresses showed aperiodic behavior regardless of the magnitude of the ratio of the viscous to the interfacial forces (Ca′), whereas for other initial orientations they were periodic. For large Ca′ the cells become highly elongated indicating the possibility of cell rupture.

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.

Theoretical analyses are developed for the deflagration of solids such as nitramines that experience exothermic reactions in liquid layers at their surfaces. Relative motion of gas and liquid in a two-phase region at the surface is considered, with influences of pressure gradients and of surface-tension gradients taken into account for the drops and bubbles. It is shown that these influences tend to produce gas velocities in excess of liquid velocities. Burning-rate expressions are derived by activation-energy asymptotics, with special attention paid to the role of interphase heat transfer.

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 two-dimensional model for foams having gas volume fraction approaching unity has been developed. A general expression for the stress tensor is obtained which gives the total stress in terms of the shape of the cells, interfacial tension, the initial cell orientation, and the rate of deformation in the liquid. A formalism for describing cell deformation is also presented. By assuming hexagonal, monodisperse foam cells, we are able to obtain analytic expressions for stresses for small shearing and elongational deformations (below the yield strain) as well as for steady shear flow. For strains below the yield point, the stress-strain relation is independent of initial cell orientation. However, the critical strain varies with orientation, and therefore, the yield stress, τ 0 , is a function of orientation. In steady shear flow, the shear stress has a contribution from the liquid which is proportional to the shear rate. Thus, τ yx = τ 0 − C μ γ . , where C is a constant determined from viscous dissipation in the thin liquid films, μ is the liquid film viscosity, and γ . is the shear rate. The yield stress is directly proportional to the liquid surface tension and inversely proportional to cell size.

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 .

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.

We present a model for nonhomogeneous solid propellant combustion that treats the nonhomogeneity as a small random perturbation of a homogeneous deflagration. Our model is based on our previous work (Margolis and Armstrong, 1986) and, unlike most previous models which are quasi-steady and one dimensional, allows for a fully nonsteady, multidimensional response of the deflagration to both spatial and temporal nonhomogeneities. In one calculation, we assume a random spatial distribution for the solid pyrolysis rate coefficient and obtain fluctuation statistics for the burning propellant in terms of the statistics of the random propellant medium. In another calculation, we obtain the burning rate response to a deterministic external pressure fluctuation. Our results are particularly applicable to nitrocellulose/nitroglycerin propellants, which are typically “almost” homogeneous.

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.