Re-Entrant Corner for a White-Metzner Fluid

: Viscoelastic ﬂuids can be difﬁcult to model due to the wide range of different physical behaviors that polymer melts can exhibit. One such feature is the viscous elastic boundary layer. We address the particular problem of a viscoelastic shear-dependent ﬂuid ﬂowing past a corner and investigate how the properties of the boundary layer change for a White-Metzner ﬂuid. The boundary layer equations are derived and the upstream layer is matched with the far-ﬁeld ﬂow. It was found that if the ﬂuid is sufﬁciently shear thinning then the viscoelastic boundary layer formulation fails due to the inertial forces becoming dominant. The depth of the boundary layer is controlled by the shear-thinning parameters. These effects are not a feature of other shear-thinning models, such as the Phan-Thien-Tanner model. This study provides insight in the different effects of some commonly used viscoelastic models in corner ﬂows in the upstream boundary layer, the downstream boundary layer is not addressed.


Introduction
Boundary layers have been well-studied for Newtonian fluids since the paradigm shifting work of Prandtl [1]. However, boundary layer structures are also ubiquitously found in non-Newtonian fluids due to the myriad of complex competing behaviors, such as those arising from the competing effects of viscoelasticity and viscosity [2,3].
The classic geometry in which to study viscoelastic boundary layer effects is the re-entrant corner, although wedge flows [4] have also been studied. Re-entrant flows arise in many rheological applications, for example, in cross-slot and upstream contraction geometries [5,6]. Initial numerical simulations of re-entrant corner flows for viscoelastic fluids were compromised due to poor numerical convergence [7]. Corner smoothing can be applied to overcome this, the effects of which have been quantified [8,9]. The first analytical approach for investigating the properties of viscoelastic stresses in corner flows was presented by Lipscomb et al. [10] for the case of a second order fluid, although this model is invalid in regimes with large velocity gradients. Other approaches have been based on a decoupling of the velocity field from the constitutive relation [11,12] to give insight into the structure of the stress singularity.
It is only since the seminal work of Renardy [13] that these boundary layer structures have started to become well-understood. Renardy argued that the local breakdown in the purely affine deformation is due to the no-slip condition sufficiently reducing the local Weissenberg number so that the viscous terms are of parity with the elastic terms. The boundary layer equations for an upper-convective Maxwell (UCM) fluid were derived and solved numerically upstream of the corner. The difficulties associated with analyzing the boundary layer downstream of the corner were addressed later [14,15].
The UCM model studied in Renardy [16] is a relatively simple viscoelastic model. Evans [17] generalized Renardy's analysis to consider an Oldroyd B fluid, thus physically accounting for solvent viscous effects not included in the UCM model. It was found that an Oldroyd B and UCM fluids exhibit similar features unless the retardation time is close to the relaxation time.
The UCM and Oldroyd B models are regarded as being accurate only for predicting the behavior of Boger fluids [18,19] as they do not exhibit the shear-dependent viscous behavior that generally occurs in polymer melts. The first studies to incorporate shearthinning effects with viscoelasticity for the re-entrant problem were performed by Sibley and Evans [20][21][22] for the case of a Phan-Thien-Tanner (PTT) fluid. It was found that as a result of the shear-thinning behavior, the velocity field decays at a slower rate compared to the Newtonian case. The stress and velocity singularities could still be expressed as functions of the corner angle, α, only. The Giesekus model, which also exhibits shearthinning viscoelastic behavior, does not permit the same boundary layer structure as the PTT model [22].
In this article, we derive the boundary layer equations for a White-Metzner (WM) fluid [23] and analyze the structure of the boundary layer. Unlike the PTT and Giesekus models, the WM model is linear with respect to the stress tensor and the parameters can be adjusted to model both shear-thinning and shear-thickening behaviors, although shearthinning effects are the more common response in viscous melts. Metzner and Whitlock [24] first reported viscoelastic shear-thickening behavior for a suspension of titanium dioxide in a sucrose solution. Viscoelastic shear thickening behavior has since been found to occur for a variety of suspensions [25,26]. Direct comparisons of the PTT, Giesekus, and WM models found that the WM model generally gave rise to the best predictions, although it did not accurately predict the uni-axial extensional viscosity [27].
The applicability of the WM model in re-entrant corner geometries has been investigated numerically in an upstream contraction [28,29]. Sousa et al. [30] studied bifurcations of polyacrylamide and polyethylene oxide in a water glycerol solvent, which was fitted to a WM model. However, a complete analytical description of the WM model around a re-entrant corner geometry has yet to be explicitly studied. This gives the motivation for this work. Determination of the scales and constraints for such boundary layers will be of importance for both further numerical and experimental investigations.
The paper is structured as follows. The governing equations are presented in Section 2. Results are presented in Section 3. Conclusions are discussed in Section 4.

Governing Equations
The geometry of the re-entrant corner problem is shown in Figure 1. Here two semiinfinite rigid plates meet at an angle φ = π α . We consider angles constrained by 1 2 < α < 1 which arises from the 2π periodicity and requirement that there is a reflex angle which will induce a stress singularity in the viscoelastic stresses. The origin of the two-dimensional plane is at the apex of the corner. The steady state mass and momentum conservation equations for an incompressible 85 fluid can be written as: The steady state mass and momentum conservation equations for an incompressible fluid can be written as: where ρ is the density, u is the velocity field, p is the pressure field and T denotes the dyadic stress tensor, the behavior of which is governed by a constitutive relation. For the WM fluid, this can be written as where ∇ T is the upper-convected derivative operator of the tensor field T, which for steady state is given by where D Dt is the substantive derivative, ∇u = ∂u j ∂x i is the tensor of the velocity derivatives for the fluid and superscript T denotes the transpose. The rate of strain tensorγ is given bẏ , where the shear-rateγ = 1 2γpqγpq with summation convention implied. In order to permit similarity solutions we will assume that the functional form of the relaxation time can be written as where K 1 and λ 0 are fixed constants which will be removed by scaling. n and q are the flow behavior indices. In order to reduce the parameter space a common assumption is to take q = n [30,31]. We refer to this as the "Bird equality" as it proves to be a useful reference point for our analysis. However, we do not need to restrict our model to this case. Several studies have found that q = n (Table 1).  [32]. The parameters for Fluid B are those fitted from a pseudo-Carreau Cross WM model to a low-density polyethylene melt [33]. Fluids C and D are polyethylene and polystyrene melts [34].
The parameters n and q can be readily obtained experimentally from shear-rheometry, with n being derived from the viscous shear relation and q from the first normal stress difference shear relation [35]. It should be noted that λ(γ) and µ(γ) could be fitted using a more generalized power-law model, such as a Carreau-type law. A Carreau viscosity model does not permit global similarity solutions in corner geometries [2], nor can it be readily approached analytically. However, these potential short-comings should have no bearing on this study. Due to the singularity in the shear rate which occurs at the apex of the corner, the Carreau model will be well approximated by a power-law relation in our regions of interest.
The WM model reduces to the upper-convective Maxwell model for n = q = 1 and to a power-law fluid model if λ 0 = 0. Thus, when both conditions are applied it reduces to a Newtonian model. To non-dimensionalize Equations (2) and (6) we introduce the scalings: As the re-entrant corner problem has no natural length or velocity scales one is free to 2 , which is equivalent to setting both the Reynolds and Weissenberg numbers to unity. In order to keep the solution as general as possible, throughout the remainder of this paper we will use scaled variables. However, the unscaled variables can be recovered using (5). Of course, a specific case study is needed in order to express the results in dimensional values. The system of equations reduces to where the only parameters which cannot be removed are n and q.

Boundary Layer Analysis Outer Solution
The upper convective derivative is assumed to be dominant in the limit as r ⇒ 0 but outside of the boundary layer [13], whereby the constitutive relation (2) reduces to ∇ T= 0.
Equation (7) has an exact solution T = h (ψ) 2 uu, where h is an arbitrary function and ψ is the stream-function [11], Renardy [13]. As the analysis is performed in the locality of the corner apex, it is assumed that h = Bψ 1 m , where B is a constant and the power m is to be determined through matching. Neglecting inertia and applying the laws of mass and momentum conservation leads to the classic inviscid Euler equations for the vector quantity v = h (ψ)u [13]: As the Euler equations permit potential flow solutions, we write v = ∂ ψ ∂y , − ∂ ψ ∂x , where ψ is a pseudo-streamfunction for v and x, y are defined as being parallel and perpendicular to the upstream boundary ( Figure 1). Using the known solution of Laplace's equation around a corner and combining it with the chain rule and the definition of v we obtain where A and B are constants. After manipulating the constants we can determine that for the outer solution where C 0 , C 1 are constants. The parameter m is still free and is determined through matching this solution to an inner boundary layer located along the upstream boundary θ = 0. As Cartesian co-ordinates provide a natural way of describing the upstream boundary, we can use the small angle approximations r ∼ x, θ ∼ y x which give the outer-behavior to match as The parameter m is subject to several constraints. The inertialess constraint is required for Equation (8) to hold. In addition, the upper-convective derivative needs to be the dominant term in (2). Both of these conditions impose restrictions on m. Analysis of the balance in the momentum equation gives u · ∇u ∼ r 2αm−3 , p ∼ h 2 u 2 ∼ r 2α−2 , and T ∼ r 2α−2 . For inertial effects to be sub-dominant results in constraint I, which is given below. Similarly, for (7) to be a valid outer approximation of from the constitutive relation (7), T T, which leads to the inequalities II and III. Constraint IV is the previously stated re-entrant condition on the angle α. In summary the conditions are given by We will a priori assume that the above constraints hold and, thus, match to the affine outer-solution. Once m has been determined we will validate the constraints I through III and discuss any limiting parameter regimes.

Upstream Boundary Layer
To analyze the boundary layer we follow the approach of Evans [14,36]. Firstly, we require that r be close to zero at the apex of the corner for (7) to be valid whereby we will assume r ∼ O(ε). However, there is a breakdown in this solution near the wall θ = 0, hence we will assume that θ ∼ O(δ/ε), from which it follows that 0 < δ ε 1. If we consider these scales in Cartesian co-ordinates we have x ∼ O(ε) and y ∼ O(δ), thus δ represents the scale of the boundary layer in the y-direction. Motivated by the outer solution, we assume the scalings: where X, Y,T ij are all O(1). Application of these scalings gives rise to the system of equations: where we employ the notation By assuming that the scaling is such that the boundary-layer equations retain the highest number of terms, it follows that from which we obtain In the Bird equality, the boundary layer size is the same as that for the UCM case, however, this is not generally so. The size of the layer is a monotonically decreasing function of n for a fixed value of q, and a monotonically decreasing function of q for fixed n. The size of the layer, δ, is bounded between ε 3−2α < δ < ε. The lower limit has a infimum of ε 2 , and, therefore, the boundary layer grows at most linearly with distance from the corner and at least parabolically. As α < 1 it is important to note that m is a monotonically increasing function of n. We need to verify that our a priori assumptions I to III still hold. The inertial constraint I reduces to Equation (21) can be used to simplify the above constraint to obtain the dependence on n solely in terms of α: From consideration of constraint IV, we see that (23) is only active for shear-thinning fluids with n < 0.25. Thus, for n > 0.25 the inertial effects are sub-dominant and are independent of the corner angle in the far field under the Bird equality. It should be noted that this condition will be different for the each of the fluids in Table 1.
Under the Bird equality constraint, II gives rise to the angle inequality α < 1, however, for n = q an additional constraint arises: This constraint is naturally satisfied for q ≥ n. In fact this condition only becomes active relative to IV if n < q + 2. (Note that constraint IV leads us to the condition (1 − α) −1 > 2.) This requires the fluid be shear-thickening. All of the fluids listed in Table 1 satisfy (24). To the best of our knowledge, there are no fluids for which (24) does not hold true. The final constraint to be verified is III, which can be written as q > −α 2(1−α) . This is naturally satisfied by q > 0.

Similarity Solution
As m has been determined and is consistent with the constraints, the system scaling can be written as It, therefore, follows that the boundary layer system can be written as As the system of partial differential equations is invariant under the scaling we seek a similarity solution of the form The similarity variable χ can be interpreted as the local Deborah number. Consideration of the local ratio of advection time to the characteristic relaxation time shows thaṫ To obtain the similarity solution, one can easily integrate the vertical component of the momentum equation to derive that the pressure is a constant, i.e.,p(χ) = p 0 . This gives the familiar result that the pressure gradient is independent of the angle θ, which drives the flow in the upstream region yet acts against the flow in the downstream region. The resulting boundary layer system becomes with matching conditions as χ → ∞ derived from (12): The relationship between p 0 and C 1 can be found from (29) and the matching conditions (33) to give p 0 = C 1 2 .

Near Wall Analysis
In order to derive consistent boundary conditions we must perform a local analysis of the system. By seeking a regular power series expansion around χ = 0 under the assumption that p 0 is fixed, one can recover the expressions where c, b ij and c ij are constants, the expressions for which are given in Appendix A.
The constant value of pressure p 0 can be expressed as a function of a and b by Equations (34)-(37) prescribe the initial conditions, however, we do not know how well-posed the problem remains. We proceed by considering an eigen-mode analysis around the regular solution. Let f = f (0) + ∆f and t ij = t ij are given by the expressions (34)-(37). Taking order ∆ leads to a system of linear ordinary differential equations inf andt ij . Three solutions can be found by searching for a solution in the form of a regular power-series: = χ 2 , t 11 = 4(n + q)|a| q+n−2 a, t 12 = 2n|a| n−1 , The remaining two independent solutions have an irregular singular point at χ = 0 and can be found by the method of dominant balance to give Here The constants S i are given in Appendix B. Let us assume that n + 3q + 1 − (2q − 1)α > 0, which imposes no additional constraints on the system (see Appendix C for details). We see that the irregular solutions mirror the same behavior as those found for the UCM case [14] for all n, q > 0, thus we conclude that the first and second modes are inconsistent with local viscometric behavior. The third mode is always consistent with the local viscometric behavior and can be determined by imposing a, thus giving one degree of freedom. The final two modes are consistent with the asymptotic condition for a < 0 and inconsistent when a > 0, independent of n and q. Thus, there are three degrees of freedom when a < 0 and one degree of freedom for a > 0. The consequence of this is that the upstream problem is well-posed with use of (34)-(37) as initial conditions near χ = 0, however, this is not the case for the downstream layer.

Upstream Boundary Layer
We now solve the parametric equations for the upstream region. As the flow is towards the corner apex a < 0, and, therefore, there is only one degree of freedom which can be accounted for by prescribing a value for a. Thus prescribing particular values for a and p 0 uniquely defines the initial value problem. We can reduce the number of parameters by noting that the system is invariant under the scaling Consequently it follows from (42) that Using the value of Γ = |a| allows us to set a = −1 for the upstream region without loss of generality. Thus, there is only one free parameter p 0 , or equivalently b. Solutions for a = −1 can then be recovered by the scalings (42). The governing Equations (29)- (32) were integrated using the MATLAB stiff ode23s solver. The asymptotic conditions were achieved by imposing the leading order term on a truncated domain [χ 0 , χ ∞ ], with the initial conditions f (χ 0 ) = 1 2 aχ 2 0 , f (χ 0 ) = aχ 0 t 11 (χ 0 ) = 2a 2 |a| n−1 , t 12 (χ 0 ) = a|a| n−1 , t 22 = 4(1−α)a|a| n−1 n+1 χ 0 . The ode23s solver requires explicit expressions for f , t 11 , t 12 and t 22 in terms of χ, f , f , t 11 , t 12 and t 22 . The latter three expressions (t 11 , t 12 , t 22 ) can be easily inverted from Equations (29)- (32) if the inversion f is given. After a little algebra it can be shown that f is given by the implicit equation where L and K can be written explicitly in terms of f , f , t 11 , t 12 and t 22 as The inversion is performed with a root finding algorithm that uses a combination of bisection, secant and inverse quadratic interpolation methods to achieve a relative tolerance of 10 −10 . The boundary layer equations are solved for both the shear-thinning (n = 0.2) and shear-thickening (n = 2) cases in Figure 2 for a = −1, p = 1, α = 2 3 , χ 0 = 10 −4 , χ ∞ = 10 4 . To show the behavior matches correctly, the solution scaled by the far-field conditions is shown in Figure 3.

Discussion
The upstream boundary layer near a re-entrant corner has been studied for the case of a WM fluid. A similarity solution possessing a similar structure to the UCM case was found. Numerical solutions were computed corresponding to both shear-thinning and shear-thickening fluids. Of particular focus was the comparison of the results of the WM model with those predicted by another shear-thinning model, the PTT model [20]. We found that for the WM model the fluid parameters n and q explicitly enter the expression for the size of the boundary layer, whereas for the PTT fluid the shear-thinning parameter does not. Similarly, the WM model predicts the existence of a critical corner angle where the inertial terms start to dominate. Such a feature does not occur in the PTT model in the upstream layer. Instead, in the case of the PTT model, the boundary layer is maintained for all re-entrant corner angles.
The analogous downstream boundary layer is more complex and poses a challenge for a future study. As the flow has reversed direction relative to the upstream problem, ψ → −ψ and a > 0. This latter constraint requires two additional degrees of freedom to account for the exponentially small terms in (41). Therefore, the downstream boundary layer cannot be well-posed as an initial value problem. In order to solve the correctly-posed problem conditions in the far-field that physically account for effects occurring upstream would need to be imposed. Such a problem would be better studied by reformulating the equations in the natural stress basis.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: