The tube model requires just two parameters for each choice of local polymer chemistry, reflecting the physics at the coarse-grained level of the tube diameter: The entanglement modulus G_e, (or alternatively the plateau modulus G) and the Rouse relaxation time of an entanglement segment _e [McLeish (2002)]. The tube diameter a is not an additional parameter, but is related to G_e via the definition of M_e [Larson et al. (2003)]. At 170 °C, we find that these are 2.0×10⁵  Pa and 7.1×10^–4  s, respectively, for PS, by fitting the molecular theory of Likhtman and McLeish (2002) to a wide data set. Treating all of the physics in the greatest detail requires solving a partial differential equation for the tensor correlation function f(s,s;t) defined in terms of the arc coordinate of the polymers R(s,t), averaged over chains as

Both the viscoelastic stress and the single chain scattering function may be calculated from a knowledge of f (the latter within a Gaussian approximation). The full partial differential equation (PDE) for f(s,s;t) contains terms that arise from advection, reptation, CLF, CR, and retraction [Milner et al. (2001); Graham et al. (2003)]:

In Eq. 2, the first term describes advection by the flow, and the second contains both reptation and CLF (from the effective local diffusion constant D^*(s,s). Here, Z is the equilibrium number of entanglement segments comprising the chain, and Z^*(t) is the time dependent instantaneous value that may differ from Z because of chain stretch. The third term arises from CR, and models the tube as a free polymerlike object with a local hopping rate . In the language of polymer physics, this is equivalent to a "Rouse theory" of the tube. The dynamics of f(s,s;t) arising from CR alone then take the form of a diffusion equation in the two variables s and s. Here, a is the "tube diameter," related to the plateau modulus directly, (s) is the local mean stretch of the chains, given by (s) = . The CR rate is calculated in turn self-consistently from averages of both reptation (diffusive CR) and retraction (CCR) over the ensemble. The only additional parameter not precisely known at present is the O(I) dimensionless constant that counts the number of local hops of a tube segment generated by one CR event. But comparison of the full theory to date, as well as theoretical considerations of the tube arising as a many-body effect (many chains intersect the volume of a single tube segment), suggests a value for this number, termed c, as c = 0.1, which we use in all calculations. For full details, see Graham et al. (2003). The final term accounts for free-Rouse chain retraction along its deforming tube. The constant R_s accounts for the decoupling approximation used in the retraction term. Graham et al. (2003) demonstrated that a universal value of R_s = 2.0 produces uniform agreement over a wide range of experimental data and we employ this value for all calculations.

Although this detailed theory is necessary to make contact with the scattering experiment, it is at present not feasible to use it locally in a finite-element calculation of the flow. Fortunately, calculation of the flow and stress fields requires only a simpler constitutive equation that is able to calculate the stress generated by a given flow history on an element of fluid. Recent work by Likhtman and Graham (2003) has shown that the full treatment at the level of f(s,s;t) may be projected onto a simpler equation for the stress only, that preserves the quality of rheological prediction. Based on the "ROuse-CCR tube model for LInear Entangled POLYmers," the "RoliePoly" constitutive equation takes the form

Here, the four terms correspond to: (1) Advection, (2) reptation, (3) chain retraction, and (4) CCR. Likhtman and Graham (2003) found that the parameter values = 1 and = –0.5 gave the closest fit to the predictions of the full model with the preferred CR parameter of c = 0.1. The negative value of is instructive: It implies that strong stretching flow suppresses CCR. One physical way in which this might arise is that the longer path length of the stretched molecules simply pick up more entanglements, so that the CCR-generated drag is greater. This is, however, not obvious: An alternative picture might view the entanglement structure simply convecting with chain stretch. The corresponding predictions of the latter are, however, not consistent with transient shear rheology. For a full discussion, see Graham et al. (2003). To recover the details of the linear spectrum of our materials, up to six modes carrying the nonlinear structure of Eq. 3 were used to make computable models of the series of monodisperse melts that approximated closely to the rheology of both the material and the full model in the flow rates of the experiment.

The time dependent flow of polymer in the MPR was simulated using the Lagrangian–Eulerian finite-element solver developed at Leeds University [Bishko et al. (1999)], using the RoliePoly constitutive model described above. In this method, the fluid velocity at each time step is calculated using a standard Eulerian finite-element technique, while the evolution of the stress is calculated in a Lagrangian frame by allowing each element to deform with the local velocity gradient. Details of the numerical method are given in Bishko et al. (1999).

The flow domain is divided into triangular elements, with velocities and pressures held at the vertices in a continuous interpolation, and material constitutive parameters ( in the RoliePoly equations, which incorporates chain orientation and stretch) piecewise constant on each element. Modes with relaxation rates much faster than the fluid velocity gradients behave as a Newtonian fluid with viscosity G_ibi, and so to reduce computational time, modes for which _di0.01 were treated as Newtonian solvent. Here = Q/d² is the typical shear rate in the constriction where Q is the area flow rate and d is the width of the channel.

At each time step, the velocity (u) and pressure (p) are found from a finite-element solution of the equations of mass and momentum conservation

Here, µ is the effective viscosity formed from the short relaxation time modes. The calculated velocities are then used to advect the grid with the flow and the corresponding triangle deformations are used to update the internal constitutive parameters for the triangles. The full RoliePoly equation was used for only the slowest mode, other modes used the nonstretching version of the equations [Likhtman and Graham (2003)]. The vertices are then reconnected as necessary to maintain a Delaunay triangulation. The reconnection process introduces a small degree of stress diffusion, but this is controlled by convergence under grid refinement. The resolution of the grid is maintained by an automated adaptive routine which divides any element whose side length is greater than a prescribed maximum length, _max (which is a function of position, so that finer grids may be specified in regions of high gradients).

At the start of the simulation, the entire domain is filled with undeformed fluid, thus the transient effects of start-up flow are simulated. Since the flow domain is symmetric, only one-half of the domain needs to be calculated. In order to reduce the size of the calculation further, a full two-dimensional (2D) model is only employed in the region from two channel widths upstream to two channel widths downstream of the contraction. The region far upstream of the contraction is modeled as start-up flow in a uniform channel using a one-dimensional (1D) simulation employing the same RoliePoly constitutive equations. These velocities are then transferred to vertices in the flow entry region of the 2D simulation. In order to mimic the small, but non-negligible compression of the upstream reservoir, a single-exponential time dependence of growth toward the steady-state upstream velocity was introduced at the upstream boundary (see below). Observation of the results indicates that deviations from uniform channel flow do not start until well into the modeled region. Thus, any upstream effects of the piston in the experimental geometry are not modeled. At the bottom of the simulated region, vertices are removed from the simulation and velocities from the uniform 1D channel flow model are imposed on the flow exit boundary points.

Following previous work [Lee et al. (2001)], the convergence of the numerical scheme was tested by using three different levels of mesh refinement in the region near the contraction. A satisfactory degree of convergence was obtained using a mesh that contained approximately ten elements spanning the narrow constriction in the center of the simulation. The stress singularity at the sharp right-angled corners at the entrance to the contraction was avoided by rounding these corners with a radius of 1% of the channel width [Bishko et al. (1999)]. Each simulation took approximately 50 h to run on a 633 MHz Pentium III processor.