Tumbling-Snake Model for Polymeric Liquids Subjected to Biaxial Elongational Flows with a Focus on Planar Elongation

We have recently solved the tumbling-snake model for concentrated polymer solutions and entangled melts in the presence of both steady-state and transient shear and uniaxial elongational flows, supplemented by a variable link tension coefficient. Here, we provide the transient and stationary solutions of the tumbling-snake model under biaxial elongation both analytically, for small and large elongation rates, and via Brownian dynamics simulations, for the case of planar elongational flow over a wide range of rates, times, and the model parameters. We show that both the steady-state and transient first planar viscosity predictions are similar to their uniaxial counterparts, in accord with recent experimental data. The second planar viscosity seems to behave in all aspects similarly to the shear viscosity, if shear rate is replaced by elongation rate.


Introduction
Understanding the behavior of polymer liquids in shearfree (extensional) flows has attracted the interest of academic researchers and industrial companies alike, due to the capacity of such flows to align and stretch polymer chains at a preferred flow direction, such as in fiber spinning and film forming processes [1]. The reliable measurement of uniaxial extensional viscosity has been resolved more than two decades ago with the development of the filament stretching rheometer [2]. Today, this rheometer has reached a level of maturity that allows to demonstrate that systems with the same number of entanglements, and thus with identical linear rheology, have a drastically different nonlinear uniaxial extensional behavior [3][4][5].
On the other hand, the measurement of the planar or biaxial extensional viscosities is rather scarce and mainly unable to reach the steady-state (see e.g., [6,7]), while the such flow fields can be generated and controlled conveniently via optical birefringence in a cross-slot channel [8][9][10]. Rheooptics is then applied to interpret the data. The unavailability of reliable direct rheological data for planar elongation may be the reason for only a few works [11][12][13][14] devoted to testing the ability of rheological constitutive models to address this flow. Non-equilibrium molecular dynamics (NEMD) simulation of microscopic polymer chain models has helped in the past to clarify the applicability of constitutive relationships for simple flows, including uniaxial elongational and shear flows, while it is worthwhile recalling that steady-state planar elongation is easier to implement than uniaxial elongation (UE) in such a simulation setup [15][16][17].
Since the introduction of the tube/reptation concept by de Gennes and Doi & Edwards [18][19][20], this mean-field theory turned out to serve as the far-most capable starting point in an attempt to describe the dynamical nonequilibrium behavior of entangled (high molecular weight) polymer melts and concentrated polymer solutions. At equilibrium, the incorporation of additional mechanisms, such as contour length fluctuations and constraint release (CR) [19,21,22], lead to an accurate description of linear viscoelastic properties [21][22][23][24]; under flow, however, and despite numerous modifications such as the consideration of chain stretch [25], finite extensibility [26,27], and convective constraint release [27][28][29][30]), it still lacks consistency with available rheological data.
Another formalism that aims to address the rheological response of high molecular weight polymeric melts and concentrated solutions is the model developed by Curtiss & Bird [31,32] based on the phase-space formulation within the kinetic theory of undiluted polymers [33]. It invokes neither a mean field tube nor slip-links. It is also known as the tumbling-snake model [34], as it allows for both orientational and curvilinear diffusion of polymer segments. The model entails, as the original tube/reptation model, the solution of a Fokker-Planck (FP) for the single-link distribution function, f (σ, u, t), which describes the probability that at time t a chain segment at contour position σ ∈ [0, 1] along the chain is oriented in direction u, with u and σ independent dynamical variables, and u · u = 1. Segmental motion is not considered as a strict one-dimensional diffusion process ("reptation") along the polymer's backbone but the chain is also allowed to explore the surrounding space by moving perpendicular to its backbone (that may be identified as CR events) with the parameter ε controlling its significance. The strictly one-dimensional diffusion process of Doi & Edwards is recovered as a special case, when ε = 0. The extra stress tensor, see Equation (1) below, contains a term due to the anisotropy of the friction tensor ζ = ζ eq [δ − (1 − ε)uu] involving a link tension coefficient ε ∈ [0, 1]; if ε = 0 there is no friction against motion in the direction u, whereas for ε = 1 the friction tensor is isotropic as for an individual sphere. Despite the qualitatively different assumptions made by the two formalisms, the original tube/reptation model is obtained as a special case of the more general FP equation of the tumbling-snake model [31][32][33]35] when ε = ε = 0. Only the analytically tractable model with ε = 0 had been solved rigorously [31][32][33]36,37].
We have shown recently that the tumbling-snake model for ε > 0 can be analyzed conveniently via Brownian Dynamics simulations and applied this approach to both steady-state [34,35,38] and time-dependent shear flow [34,38], as well as to steady-state and time-dependent uniaxial elongation [39]. These works provided evidence that the tumbling-snake model is able to capture the damping behavior of the transient viscosity in start-up shear experiments at high rates [40][41][42], while preserving the absence of such undershoots in both normal stress coefficients, in line with experimental data [34,38]. The appearance of the undershoot has been associated with the shear-induced rotational motion of chains [38,42], further supported by non-equilibrium atomistic simulations [43][44][45]. As such, similar undershoots are not seen in elongational flows [39].
The qualitatively relevant and only modification to the original tumbling-snake model was the consideration of a variable link tension coefficient, that vanishes in the absence of flow, and is given by ε = ε 0 S 2 2 [34,38], where S 2 denotes the 2nd rank uniaxial nematic order parameter of polymer segments [46]. This adjustment of the Curtiss & Bird theory has eliminated certain disadvantages of their original model (exhibiting a constant link tension coefficient). Due to the refinement, the transient shear and elongational viscosities no longer approach constant values at small times, and spurious time oscillations of the transient second normal stress in startup of shear flow are absent. It has been demonstrated that the tumbling-snake model in its present form is able to qualitatively capture recent experimental evidence according to which the extensional viscosity of polymer solutions is seen to exhibit thinning below the inverse Rouse time and thickening above, whereas the extensional viscosity of polymer melts is monotonically decreasing for all strain rates [3][4][5], by having the strength ε 0 of the link tension coefficient increasing as the polymer concentration decreases [39].
In this work, we discuss the solution of the tumbling-snake for the more general case of biaxial elongational flows, with a focus on planar elongational flow. The structure of this manuscript is as follows: In Section 2 we revisit the stress tensor of the tumbling-snake model, parameterize the velocity gradient tensor and define the viscosities. Section 3 summarizes the Brownian dynamics method to solve the model. In Section 4, we provide the series expansion of the two biaxial elongational viscosities in the case of steady-state general shearfree flow for small dimensionless elongation rates for comparison with limiting results presented in Section 6. Similarly, in Section 5 we derive analytic expressions for the linear viscoelastic viscosities in the case of start-up flow, again for biaxial flows. In both sections, we discuss the cases of constant and variable link tension coefficient separately. This includes special cases such as the rigid dumbbell. In Section 6 we actually solve the model numerically using Brownian dynamics simulation for planar elongation, validate the analytic expressions, and further compare the predictions of the first planar viscosity with the uniaxial elongation viscosity. We conclude with Section 7 where we discuss the significance of this work.

Stress Tensor
In the case of a monodisperse polymer with "polymerization degree" N related to the number of entanglements per chain, Z, introduced by Doi and Edwards, via N = 3Z, and polymer number density n, the time-dependent (extra or polymeric) stress tensor τ of the tumbling-snake model subjected to a homogeneous flow field characterized by the transposed velocity gradient tensor κ is given by [33,35] with modulus G = nk B T(N − 1), temperature T, Boltzmann's constant k B , unit tensor I, coefficients ε and ε 0 interrelated via ε 0 ≡ ε (N − 1) 2 , and a link tension coefficient ε. The latter comes in two versions, the original one proposed by Curtiss & Bird, where ε is a constant, and the variable one proposed by us within the tumbling-snake model, ε = ε 0 S 2 2 with constant ε 0 and uniaxial order parameter S 2 determined by uu (1) (t). Note that we are adopting throughout the nomenclature of Bird et al. [33], while the τ in (1) is a pressure tensor, and thus the negative stress tensor, in the majority of scientific literature. This stress tensor (1) involves the following orientational averages calculated with the solution of the corresponding FP equation [34] for the single-link orientational distribution function f (σ, u, t) uu (1) where du denotes an integral over the unit sphere, λ a time constant proportional to ζ eq /k B T, squared bond length a, and N 3+β , where β is a chain constraint exponent, and κ : uu = (κ · u) · u stands for a two-fold contraction. The reptation or disengagement time is τ d ≡ λ/π 2 . In addition to the highlighted dependency on time t the stress tensor as well as the averages depend also on the flow field κ via f . For the case of a general shearfree elongational, homogeneous incompressible flow at rate˙ the transposed velocity gradient tensor κ is traceless and diagonal and thus of the form where it is sufficient to consider semipositive b for symmetry reasons. When b = 0, we retrieve homogeneous simple uniaxial elongation for˙ > 0 and biaxial stretching for˙ < 0, while b = 1 corresponds to planar elongation [33], b = 3 to equibiaxial extension, and b = 2 to so-called elliptical extension ( Figure 1). Arbitrary b can be realized experimentally via a multiaxial elongational rheometer [6]. While most results and methods to be presented below are valid for arbitrary b, we will focus on planar elongation (b = 1) in Section 6. The ratio of principle strain rates κ 22 /κ 33 , denoted by m by Demarmels and Meissner [6], is related to b by m = (b − 1)/2, and the projection of the motion of a material particle in the xzor yz-plane is given by x = Cz 1−m and y = Cz m with constants C, while y(t) = x(t)e bt in the xy-plane. Except for b = 0 (κ xx = κ yy ) and b = 3 (κ yy = κ zz ) there are two normal stress coefficients and corresponding viscosities that can be measured. The first, η 1 (˙ ) = −(τ zz − τ xx )/˙ , and the second, η 2 (˙ ) = −(τ yy − τ xx )/˙ , rate-dependent stationary elongational viscosity [33]. The corresponding transient viscosities we denote by η + 1 (˙ , t) and η + 2 (˙ , t). The transient viscosities in the linear viscoelastic regime do not depend on rate, are thus denoted by Because the elongation rate enters the stress tensor in the combination λ˙ , we introduce the dimensionless Weissenberg number

Brownian Dynamics Simulation
The Brownian dynamics algorithm that we employ in this work is identical, apart from the different choice of flow field, κ, with the one we had presented previously. The Fokker-Planck equation of the tumbling-snake model subject to isotropic boundary conditions at chain ends at all times, ∀ t f (0, u, t) = f (1, u, t) = 1/4π, can be cast in the form of two coupled Itô stochastic differential equations for variables U t (segment unit vector at time t) and σ t ∈ [0, 1] (relative contour position at time t) as follows where dW t is a three-dimensional Wiener process and dW t is another one-dimensional Wiener process, independent from the former ( Figure 2). The transverse projector operator I − U t U t guaranties that the stochastic dynamics preserves the U t property of being a unit vector. Note that the link tension coefficient ε affects the stress tensor, but not the dynamics of the orientational distribution function. From Equation (5) it is transparent that 1 − ε is related to a one-dimensional "reptation" diffusion coefficient, and ε 0 = ε (N − 1) 2 related to the orientational diffusion coefficient of polymer segments. The factor (N − 1) 2 appears because σ t is a relative rather than absolute contour position. For implementation details see [34,35].
The moments required to evaluate the stress tensor given by Equation (1), such as uu (1) (t), i.e., the left hand sides of Equation (2), are obtained during Brownian dynamics via plain averaging over an ensemble of stochastic trajectories at times t. To be specific, uu (1) involving the evaluation of a 2nd and 4th rank tensor, whose number of nonvanishing and independent components depends on κ (2 and 6 components, respectively, for the case of biaxial elongation).

Small Elongation Rate Expansion of the Stationary Viscosities for Biaxial Elongational Flows
Results for the stationary viscosities at small rates can be derived analytically. They are particularly useful because this limiting case can strictly not be accessed during Brownian dynamics, because the error bars increase with decreasing rate. They are furthermore useful to, e.g., test ratios between zero rate viscosities or to compare asymptotic behaviors for different types of flow. The approach to derive analytical results is based on a spherical harmonics expansion of the single-link distribution function around equilibrium.

Stationary Regime, Constant ε
To begin with, we assume a constant ε and we are interested in the rate-dependent steady-state values of the extensional viscosities. The methodology employed is described in detail in the Supplementary Section A and the final expression for the expansion, up to 2nd order in the dimensionless Weissenberg number Wi =˙ λ, is given in Supplementary Equation (S6). Upon inserting this expansion into the stress tensor Equation (1) we obtain both elongational viscosities up to second order in Wi in terms of the biaxiality parameter b or alternatively, if we normalize with the known zero-rate shear value, η 0 = 1 60 1+ 2 3 ε Gλ [33,35], the result can also be written as This way we see that the first and second zero-rate elongation viscosities lim Wi→0 η 1,2 (˙ ) are 3 + b and 2b times, respectively, the zero-rate shear viscosity, η 0 . The following abbreviations have been introduced for numerical prefactors appearing in (6) and (7) with the kernels k 1 (ν) = K 2 , k 2 (ν) = K 2 2 , and k 3 (ν) = K 2 K 4 that depend on both ε and ε 0 via Evaluating the ∆ j 's we can obtain more explicit predictions for two limiting cases (i) and (ii): (i) When ε = 0, implying K j (ν) = (πν) 2 , the ∆ j are readily evaluated, ∆ 1 = 1/40 and ∆ 2 = ∆ 3 = 17/6720, and Equation (6) reduces to up to order Wi 3 . These expressions (10) further include the Doi & Edwards results for ε = 0. (ii) In the second limit, ε = 1 with N = 2, thus ε 0 = ε , the chain reduces to a rigid dumbbell. For this case K j (ν) = j(j + 1) and all kernels k j are independent on ν, leading to ∆ 1 = 1/24, ∆ 2 = 1/144, and ∆ 3 = 1/480. We thus obtain from Equation (6) with Wi rd ≡˙ λ rd , G = 6G rd and λ = 6λ rd . To the best of our knowledge, this expansion for the rigid rod subjected to biaxial flows is provided here for the first time. Our result, Equation (11), also accounts for hydrodynamic interaction by identifying ε = λ

Stationary Regime, Variable ε
If, instead of a constant link tension coefficient, we consider the variable link tension coefficient of the tumbling snake model given by ε = ε 0 S 2 2 where S 2 2 = 3 2 tr( uu ani · uu ani ) is the uniaxial order parameter of the 2nd moment (here uu ani = uu (1) − 1 3 I) [34,38], then, up to second order in Wi we obtain where Γ 1 is a numerical coefficient that appeared in the Supplementary Equation (S1b) of Ref. [38]. The corresponding steady-state viscosities are given, up to O(Wi 3 ), by We refrain from writing down more explicit expressions for the special cases of (i) ε = 0, using the ∆ j 's mentioned above for this case, as well as Γ 1 = 1/8, and (ii) ε = ε 0 = 1, using Γ 1 = 1/4.

Transient Elongational Viscosities in the Linear Viscoelastic Regime
To study the transient viscosities after startup of flow we consider a time-dependent spherical harmonics expansion of the single-link distribution function around equilibrium, up to first order in Wi, to be able to obtain the linear viscoelastic analytical predictions; the procedure is described in the Supplementary Section A; the final expression for the expansion of the time-dependent single-link distribution function is given by Supplementary Equation (S4).

Transient Regime, Constant ε
Inserting this expansion into the stress tensor Expression (1), assuming a constant ε, we obtain analytical expressions for the time dependent viscosities, η + 1 (t), and η + 2 (t), which turn out to be 3 + b and 2b times the time-dependent shear viscosity, that was first presented in [34], where the following abbreviation has been introduced, with K 2 (ν) from Equation (9). An important property of ∆ 0 (t) is its initial value, ∆ 0 (0) = 1/4. It decreases monotonically with time, initially with a slope of −24[(2 − ε ) + ε 0 ]/λ, and ultimately vanishes as t → ∞. Taking the rigid dumbbell limit, ε = 1 and N = 2, the ∆ 0 (t) can be evaluated analytically and Equation (15) becomes To the best of our knowledge, this is the first time this result for a rigid dumbbell is presented.

Transient Regime, Variable ε
If, instead of a constant link tension coefficient, we consider a variable link tension coefficient given as ε = ε 0 S 2 2 , the least order expansion with respect to Wi gives where the dimensionless Γ 1 (t) = 12 ∑ ∞ ν=1,odd exp(−K 2 t/λ)/(πν) 2 K 2 with Γ 1 (0) = Γ 1 , cf. Equation (13), is taken from the Supplementary Equation (S3b) of Ref. [38]. Putting this together, ε at small rates adn times increases quadratically with rate and time. The precise expression is As the variable link tension coefficient thus vanishes in the linear viscoelastic regime, the full time dependent planar elongation viscosities are given by Equation (15) evaluated at ε = 0, i.e., For times t λ, this expression reduces, with ∆ 0 (t) given by Equation (16), to where use had been made of the already mentioned initial slope of ∆ 0 (t).

Brownian Dynamics Results for Planar Elongational Flow
Having derived analytical expressions for the various regimes, we now turn to the presentation of full rate-and time-dependent exact numerical results for the tumbling-snake model for the special case of planar elongational flow (b = 1).
where polymers tend to align in z-direction to the expense of a compression in x-direction, while y is the neutral direction. All figures presented in this manuscript are generated using the variable link tension coefficient; predictions for the case of a constant ε are available in the Supplementary Section B. All types of biaxial flows discussed above can be studied using an identical procedure. Mixed flows pose no additional problem and amount to consider a more general κ or even a time-dependent κ(t) in Equation (3) of the Brownian dynamics algorithm (5). The analytical results will be used to test the simulation results, and to extend their validity to "infinitely" small rates and times, where simulation results tend to become more difficult to sample.

Steady-State Planar Elongation
The steady-state value of the variable link tension coefficient ε = ε 0 S 2 2 as a function of dimensionless elongation rate Wi for polymerization degree N = 100 (Z ≈ 33 entanglements) and various values of ε 0 = ε (N − 1) 2 is shown in Figure 3. At small elongation rates, as expected from Equation (12) ε increases quadratically with the elongation rate, whereas at large Wi, and irrespective of the value of ε 0 , then ε → ε 0 . This is also expected since for a fully aligned sample the order parameter approaches unity, S 2 → 1. Thus, the model predictions for the case of a variable ε become identical to the ones for a constant ε at large rates, Wi 1.  (23) when ε 0 = 0 (dark blue) and 0.9 (dark yellow). The dependency of Γ 1 on ε 0 determines their offsets and is relatively weak for any N. For N > 2 the coefficient Γ 1 decreases monotonically with increasing ε , starting from Γ 1 = 0.125 at ε = 0. For any N ≥ 20 one has Γ 1 ≈ 0.081 at ε = 0.9 (and Γ 1 ≈ 0.078 at ε = 1).
The reduced steady-state first viscosity, η 1 (˙ ), as a function of the dimensionless elongation rate is presented in Figure 4. All solid lines in this figure are determined by Equation (14) evaluated at b = 1, i.e., with numerical prefactors ∆ j s given by Equation (8). The corresponding predictions of the simplified tumbling-snake model (ε = 0) read, where we have used the numerical values for ∆ 2 and ∆ 3 quoted after Equation (6), and Γ 1 = 1/8 for ε = 0.  Figure 4a shows the variation of η 1 (˙ ) upon changing ε 0 (or ε ) while keeping N = 100 and ε = 0 fixed. The prediction at small elongation rates is independent of the value of ε 0 and approaches the value 4η 0 . However, as Wi increases the first planar viscosity shear-thins after about Wi ≈ 1. This is in disaccord with the predictions for the uniaxial elongation viscosity of the tumbling-snake model, η E , which seems to be passing from a maximum when ε 0 > 0 (see inset of Figure 2a of Ref. [39]). However, both share the same power-law behavior at large elongation rates which is always equal to −1 irrespective of the value of ε 0 , as for the Doi & Edwards model [18].
In Figure 4b we show the same variation as in Figure 4a but now with ε 0 = 0.1. Again, η 1 (˙ ) follows Equation (24), or Equation (25) for the case of ε 0 = 0, at small elongation rates and, irrespective of the value of ε 0 , reaches monotonically the same value at large elongation rates. This value is simply equal to where η 0 = Gλ/60 is the zero rate shear viscosity for the case of variable , that vanishes at vanishing rate. The reduced first planar viscosity thus drops (or rises) from a value of 1/15 at small rates to ε 0 /6 at infinitely large ones. For the scenario ε 0 = 0.1 shown in Figure 4b, the reduced first planar viscosity reduces to η 1 (∞)/Gλ = 1/60. Equation (26) can be readily derived by noting that as Wi → ∞ then lim Wi→∞ ε = ε 0 as well as These expressions hold as long as κ zz is the largest diagonal component of κ, which is the case for any b < 3. Using similar arguments, η 2 (∞) = 0 for b < 3. Equation (27) are analogous to the case of uniaxial elongation [39], apply independently of the value of ε 0 , and originate from the dominance of the third term in the stress tensor expression, Equation (1), at large elongation rates leading to a leveling-off of the first viscosity at a value given by Equation (26). By further increasing the value of the parameter ε 0 , for given ε 0 and N, the small elongation rate predictions are unaffected (Figure 4c,d). At large Wi the curves always reach the value of the reduced η 1 (∞), Equation (26). When the value of ε 0 exceeds 2/5 then η 1 (∞) > 4η 0 (Figure 4d). The corresponding predictions when the link tension coefficient is considered constant, but non-vanishing, are given in Supplementary Figure S1.
From a visual comparison between planar η 1 (˙ ) and uniaxial η E (˙ ) (see Figure 2 of Ref. [39]) we note that their predictions are similar. To better quantify the similarities between η 1 and η E we compare the two in Figure 5 in a way so that they both have the same zero-rate prediction. We show the comparison upon changing ε 0 whilst keeping N = 100 and ε = 0 (Figure 5a) or ε 0 = 0.1 (Figure 5b) fixed. We find that η E (˙ )/3Gλ is slightly exceeding η 1 (˙ )/4Gλ in the intermediate flow regime, starting at about Wi ≈ 3. Further, η 1 (˙ )/4Gλ for ε 0 = 0.5 is basically coinciding with η E (˙ )/3Gλ for ε 0 = 0. For the case ε 0 > 0 in Figure 5b we note that at large elongation rates η E (∞)/3Gλ reaches a constant value larger than the corresponding one for planar elongation, η 1 (∞)/4Gλ, which stems from the way the two viscosities were made dimensionless; if both were made dimensionless with Gλ the corresponding value would be the same for both viscosities. The reduced steady-state second viscosity, η 2 (˙ ), as a function of the dimensionless elongation rate is presented in Figure 6. As for η 1 (˙ ), the prediction at small elongation rates is independent of the values of ε 0 and ε 0 in accord with Equation (24), or Equation (25) for the case of ε 0 = 0, and approaches the value 2η 0 . Also, like η 1 (˙ ), as Wi increases the second viscosity shear-thins after about Wi ≈ 1. It can be noted that the predictions when ε 0 is kept fixed (panels a and b) are almost the same irrespective of the value of ε 0 and the power-law behavior at large elongation rates when ε = 0 is equal to −3/2 as for the Doi & Edwards model. In panels Figure 6c,d we note that as the ε 0 increases under fixed ε 0 , η 2 (˙ ) increases after about Wi ≈ 20.

Transient Planar Elongation
Next we inspect the transient link tension coefficient, ε/ε 0 , as a function of dimensionless time t/λ for N = 100 and various values of the parameter ε 0 and dimensionless elongation rates Wi (Figure 7). At early times, this coefficient follows 12 25 (˙ t) 2 according to Equation (19) with b = 1, irrespective of Wi and ε 0 , whereas at larger times it monotonically approaches the steady-state values of the squared order parameter S 2 2 . A similar behavior was also noted for the transient behavior of the variable link tension coefficient in uniaxial elongation [39]. In Figure 8 we show the transient first viscosity η + 1 (˙ , t) as a function of the dimensionless time for various dimensionless elongation rates along with the linear viscoelastic prediction, that follows from Equation (20) The first term on the right hand side provides the full dependency on time, the second the initial slope. For all elongation rates we indeed notice from our Brownian dynamics results that as t → 0 the transient first viscosity η + 1 (˙ , t)/Gλ increases with increasing ε 0 , as suggested by the slope provided by Equation (28). As was the case for the shear viscosity [34,38] and the uniaxial elongation viscosity [39], by using the variable link tension coefficient ε = ε 0 S 2 2 we have amended the problematic predictions of the original tumbling-snake model, in which a constant link tension coefficient is considered, according to which both limiting η + 1 (t)/Gλ and η + 2 (t)/Gλ approach a constant value, ε 0 /15 and ε 0 /30, respectively, irrespective of the value of the parameter ε 0 and thus the mode, reptation versus orientational diffusion. Additional Figures S3 and S4 are provided in Supplementary Section B. By further increasing the elongation rate (Wi = 100) the transient first viscosity goes over the linear viscoelastic prediction only when ε 0 > 0, and reaches the steady-state value without reaching an overshoot, independently of the value of the parameters ε 0 > 0 and ε 0 . We observed a similar trend in the case of uniaxial η + E (˙ , t) [39]. Similarly for larger elongation rate (Wi = 1000) but now the curves depart much sooner and more intensely from the linear viscoelastic prediction, just like η + E (˙ , t) does depart from η + E (t). The non-appearance of an undershoot in shearfree flows is attributed to the absence of a rotational contribution to κ, and because the orientational diffusion of segments does therefore not lead to any deterministic rotation. On the other hand, under shear the tumbling behavior of polymer chains is not only due to the rotational component of κ, but triggered by the orientational diffusion term, present for ε 0 > 0. It produces undershoots at large shear rates [34,38]. Overall, the predictions for the transient first viscosity exhibit the same time-dependency with the uniaxial elongation viscosity except at large times as they approach their steady-state values ( Figure 9).
Finally, we investigate the transient second viscosity η + 2 (˙ , t) as a function of the dimensionless time for various dimensionless elongation rates along with the linear viscoelastic predictions given by Equation (28) in Figure 10. Like η + 1 (˙ , t), for all elongation rates and as t/λ → 0 the transient second viscosity η + 2 (˙ , t) increases with increasing ε 0 , following Equation (28). Also, at small times and irrespective of the elongation rate and the values of the parameters ε 0 and ε 0 , η + 2 (˙ , t) follows the linear viscoelastic prediction, Equation (28). It implies that the ratio between the two viscosities is two, initially. The most important distinctions between the transient behavior of η + 1 (˙ , t) and η + 2 (˙ , t) are the appearance of an overshoot, independently of the values of the parameters ε 0 and ε 0 , and the fact that for ε 0 > 0 and any ε 0 the curves are only slightly above the linear viscoelastic curve at small times. As for η + 1 (˙ , t), an undershoot is not predicted, as expected. Finally, it should be noted that, as was the case for shear flow and UE, model predictions for constant values of the parameters ε 0 and ε 0 but with a different number of Kuhn segments N are found to be identical, for both steady-state and transient quantities, for large values of N (N ≥ 10), since the two viscosities were scaled with the modulus G and the relaxation time λ, both of which do depend on N; thus, this comparison is not shown.

Conclusions
In this work, we discussed the features of the tumbling-snake model for concentrated solutions and entangled polymer melts subjected to both steady-state and transient biaxial elongation, focussing on planar elongation as an application of the more general framework. The model employs a variable link tension coefficient, given by ε = ε 0 S 2 2 [34,38], which, as for shear flow and uniaxial elongation, has amended several shortcomings of a constant link tension coefficient originally suggested by Bird et al. [33,37]. In particular, the two planar transient elongation viscosities η + 1,2 (˙ , t) do not approach a finite value as t → 0 upon using the variable link tension coefficient. We have shown that the predictions of the first planar viscosity η 1 and the uniaxial elongation viscosity are similar with respect to their transient and stationary behaviors, in accord with available experimental data [7]. On the other hand, the steady-state second planar viscosity η 2 always thins irrespective of the value of the ultimate strength of the link-tension coefficient, ε 0 . Overall the second viscosity seems to share many features with shear viscosity with respect to rate and time, except that η 2 , unlike the shear viscosity [34,38], does not (and should not) exhibit an undershoot in the course of time after startup of flow.
As a concluding remark, the tumbling-snake model with variable link-tension coefficient has been shown to provide a very adequate description of the available rheological measurements of entangled polymer melts and concentrated polymer solutions when subjected to shear [34,35,38], uniaxial elongation [39], and planar elongation. The model is tractable analytically at small and large rates, that are unaccessible or difficult to access by numerical inspection, and for all intermediate rates the model is solved conveniently by simple Brownian dynamics. Introducing further refinements, such as contour length fluctuations (see e.g., [23,47,48] and references therein), convective constraint release [27][28][29], flow-induced alignment of chain ends [46,49], and chain stretch [25,50], the latter being significant in strong elongation flows, could further improve the tumbling-snake's model capacity to quantitatively predict the rheological response of entangled polymer melts and concentrated polymer solutions. As for our works preceding the present study, no such refinements were attempted to present predictions of the tumbling-snake model for future reference. This includes our analytical results for biaxial elongation in the case of purely one-dimensional diffusion, ε 0 = 0, or the transient viscosities for a rigid dumbbell, as such results were apparently not available so far.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4360/10/3/329/s1, Section A: Methodology to obtain the real spherical harmonics expansion of the single-link distribution function, Section B: Results of the BD simulations when a constant link tension coefficient is employed (Figures S1-S4), Figure S1: Steady-state first planar viscosity, Figure S2: Steady-state second planar viscosity, Figure S3: Transient first planar viscosity, Figure S4: Transient second planar viscosity.