Revisiting boundary layer flows of viscoelastic fluids

In this article we reconsider high Reynolds number boundary layer flows of fluids with viscoelastic properties. We show that a number of previous studies that have attempted to address this problem are, in fact, incomplete. We correctly reformulate the problem and solve the governing equations using a Chebyshev collocation scheme. By analysing the decay of the solutions to the far-field we determine the correct stress boundary conditions required to solve problems of this form. Our results show that both the fluid velocity within the boundary layer and the stress at the solid boundary increase due to the effect of viscoelasticity. As a consequence of this, we predict a thinning of the boundary layer as the value of the dimensionless viscoelastic flow parameter is increased. These results contradict a number of prominent studies in the literature but are supported by results owing from an asymptotic analysis based on the assumption of the smallness of the non-dimensional viscoelastic flow parameter.


Introduction
Interest in the theory of viscoelastic boundary layer flows first developed in late 1950s and early 1960s, some 50 years, or so, after the fundamental contributions of [1] and [2]. The early works of [3], [4], and [5] sought to develop the theory of, and approximate analytical solutions to, the stagnation point boundary layer flow of viscoelastic fluids described using Rivlin-Ericksen tensors [6].
Ultimately, each of these early studies suffered either from a misunderstanding of the governing physics and/or mathematical shortcomings. It wasn't until the slightly later study of [7] that the first complete, correct, boundary layer analysis concerning viscoelastic fluids was completed. In this study the authors begin by considering an Oldroyd-B constitutive viscosity law and then choose to restrict their attention to fluids with short relaxation times. This restriction reduces their analysis to consider what we would now refer to as a 'second-order' fluid model. In the absence of advanced numerical techniques, the authors develop approximate analytical solutions using an asymptotic expansion based on their assumption of the smallness of the non-dimensional viscoelastic flow parameter. They conclude that as viscoelasticity is increased, so does the stress at the solid boundary and, subsequently, the boundary layer is thinned when compared to its Newtonian counterpart.
Following the work of [7] there have been a plethora of studies concerning boundary layer flows of viscoelastic fluids characterised by the second-order fluid model (sometimes referred to as 'fluids of second grade'). Rajagopal and co-authors have produced numerous studies of this nature, see for example [8], [9], [10], [11], [12]. In these works the authors cover the theory of second-order fluid boundary layer flows, use the perturbation method introduced by [7] to analyse flows over wedges, investigate the effects of suction through a permeable surface, use the idea of an augmented boundary condition to arrive at new solutions to the stagnation point flow problem, and, lastly, show that their stagnation point analysis can be extended, in a local sense, to flows over wedges.
The use of Garg and Rajagopal's augmented boundary condition for problems of this nature has been discussed at some length in the literature. Indeed, we will subsequently be required to revisit these discussions as part of our analysis. At this stage, we refer the interested reader to the work of [13], wherein the vast majority of suitable references on this discussion point can be found.
Seeking to improve upon the aforementioned second-order fluid analyses, in the mid 2000s, Sadeghy and co-authors focussed their efforts on the study of viscoelastic boundary layer flows with an upperconvected Maxwell (UCM) constitutive viscosity law ( [14], [15]). Given the relative limitations of the second-order fluid model, the motivation for these studies came from a want to capture the correct physics of these problems for larger values of the viscoelastic flow parameter. However, in both studies, the authors relied upon the problem formulation presented in Harris' book 'Rheology and Non-Newtonian Flow ' ( [16]). Careful examination of this reference shows that Harris made a rather a significant oversight when deriving the equations relevant to UCM boundary layers. As such, the analysis presented by Sadeghy et al. fails to correctly capture the effects of viscoelasticity for high Reynolds number flows. This realisation is important, not only because of the need to produce mathematically correct solutions, but also because of the number of subsequent studies that have used this formulation type as a basis for their investigations.
To date, both [14] and [15] have received well over 100 citations each, with studies involving flows over porous surfaces ( [17]), flows involving the transfer of mass ( [18]), magnetohydrodynamic flows ( [19]) and rotating flows ( [20]), to name but a few, relying, via Sadeghy et al., on the mathematical shortcomings first presented by [16]. Interestingly, we note that there is a far smaller readership for the work of [21], published well in advance of both [14] and [15], who use the correct equation set, contrary to those described by [16]. This study conducted in the mid 1980s considers the stagnation point flow UCM problem, and is overlooked by Sadeghy and co-workers. The corresponding problem of the flow of an Oldroyd-B fluid is reported in [22], and we note that a similar case study concerning the Giesekus fluid model is conducted by [23]. Comparisons of our work with all these investigations is presented in §4.
Studies of boundary layer flows of other common viscoelastic models have been undertaken more recently. For example, [24] and [25] consider a FENE-P constitutive viscosity law, which seeks to improve on the Oldroyd-B model by introducing an elasticity which is finitely-extensible. These investigations focus on a forced convection boundary layer flow, and the calculation of local self-similar solutions for the flat plate boundary layer problem, respectively. We also note the work of [26], who seek to extend the analysis presented by [25]. They combine semi-analytical results, in the sense of local self-similarity, with full numerical results owing from OpenFOAM simulations, showing a good agreement between the two sets of solutions. We go on to show in §2 that the flat plate geometry does not yield self-similar solutions for the Oldroyd-B model, and this is indeed the case for the FENE-P model, as was originally shown by [25].
The motivation for this study stems from our work concerning the injection of non-Newtonian fluids into otherwise Newtonian boundary layer flows. In an effort to validate the results of our injection analyses, for flows involving viscoelastic fluids, we explored the large injection limit expecting our solutions to tend towards the known published results for entirely viscoelastic boundary layers. When this was not the case, we analysed the literature regarding flows of this nature and noted the mathematical shortcomings mentioned above. Before presenting the results of our much broader injection study we felt it important to revisit, and correct, the current state of the literature.
In this article, we have improved on the aforementioned studies in a number of ways. Firstly, where previously the necessity of a stagnation point flow has been assumed for flows of a viscoelastic fluid, we have worked independently of this assumption, and show analytically that one is restricted to this case study in order to obtain fully self-similar solutions. Secondly, we investigate the choice of stress boundary conditions in a full and rigorous manner, where in other articles the conditions are often simply stated or assumed without proper reference. In addition to this, we provide a detailed explanation of our numerical scheme, which we will show is highly accurate 1 . Lastly, our solution method benefits from the possibility of evaluating moderate levels of viscoelastic effects. Previous studies [15], [21], [22], [23] only show results for a relatively low ratio of elastic to viscous effects, in some cases due to the imposition of a boundary condition which restricts the validity of the range of their dimensionless parameters.
The outline of this article is as follows, in §2 we formulate the problem, firstly in a general sense, and then derive the equations relevant to the flow of viscoelastic fluids described using the Oldroyd-B model, which similarly describes the UCM model by fixing a single parameter. In §3 we outline the numerical scheme used to solve the governing equations and present a range of results for fluids with varying levels of viscoelasticity. In §4 we compare our results to those in the literature and discuss our findings in this context. Lastly, in §5 we provide a brief summary and outline some future directions for work on problems such as this.

Formulation
Consider the steady flow of viscoelastic fluid over an impermeable, semi-infinite, flat plate inclined at an angle of mπ/(m + 1), from the horizontal where m is a constant that will be defined in due course. The streamwise coordinate is x * and the wall normal coordinate is y * (asterisks denote dimensional quantities throughout). The flow is governed by the continuity and Cauchy momentum equations where ρ * is the fluid density, p * is the pressure, u * = (u * , v * ) is the two-dimensional velocity field with u * and v * being the streamwise and wall normal velocity components respectively. Lastly, we note that the definition of the stress tensor τ * , changes according to the form of the viscoelastic model in question. In this study we will focus our attention on three prominent viscoelastic models, these being the Oldroyd-B model, the upper-convected Maxwell model and the second-order fluid model. For the types 1 Data created during this research can be found online at a location provided in C of flows discussed in this article all, three models can be described using the following constitutive viscosity law where µ * 0 is the total viscosity, E * = ∇ * u * + (∇ * u * ) T /2, is the rate of strain tensor, λ * 1 is the relaxation time, and λ * 2 is the retardation time. We note that the velocity gradient is defined with the following convention in index notation: (∇ * u * ) ij = ∇ * i u * j . The upper-convected derivative for a (steady) tensor q * is defined as follows which describes how a quantity is translated and rotated under the influence of the flow field. In the case when the total viscosity, relaxation time and retardation time are all non-zero, (2) returns the Oldroyd-B model. If the retardation time is set equal to zero then this model reduces to the upper-convected Maxwell model. In the limit of short relaxation times the equation of state (2) reduces to [7]). Now, the second-order fluid model is defined like so (see [27]) While we consider the use of the upper-convected derivative, it is noted that other sources in the literature [28], [29], [30] use either the lower-convected or corotational derivative in order to describe the second-order fluid model. All definitions are equivalent so long as the material constants for the model are defined in an appropriate manner. Upon comparison of the definitions of the stress tensors above, we see that the limiting case of an Oldroyd-B fluid in short relaxation times does not capture the full general second-order fluid model [27], rather a specific case study. The simplest way of quantifying this is to set the constants α * 2 = 0, with −κ * 0 = α * 1 . However, we show in A that the full second-order fluid model is captured in our boundary layer analysis, since the quadratic rate of strain term, E * · E * , results in no extra contribution to the terms associated with either the velocities or viscoelastic stresses. Rather, it provides a contribution to the isotropic pressure term. Therefore, the final set of boundary layer equations are independent of α * 2 under the relevant limit of time scales. Here we solve the coupled system of equations ((1) -(2)) by explicitly considering the independent stress components of τ * in two-dimensions. However, [16] uses a separate and incomplete technique (in the UCM case study where λ * 2 = 0), which involves the decoupling of velocity and stress. Firstly, Harris re-arranges equation (1b) to define ∇ * · τ * , the stress divergence, precisely in terms of the velocity and pressure. Then, the divergence operator is applied to (2), with the aim of capturing all stress dependent terms as a function of the same vector function, the stress divergence. One can see that, if this were possible, the definition of ∇ * · τ * could be inserted into this new equation such that we are left with one equation in velocity and pressure only. However, it is not possible to describe the viscosity law as such, since the upper-convected derivative and divergence are not commutative operators: where Q * = ∇ * · q * , for any general quantity q * . We note that the vector equivalent of the (steady) upper-convected derivative is defined as To better visualise the inequivalence, we highlight a simple case study, where the quantity q * represents an isotropic tensor I * . The upper-convected derivative of I * is some multiple of the rate of strain tensor; since the identity tensor is homogeneous it has no spacial or temporal dependence, and therefore the convective term depends entirely on the velocity gradient. It must also be the case that the upperconvected derivative is frame invariant, which reduces the quantity to the symmetric component of the velocity gradient i.e., the rate of strain tensor. After the explicit calculation, we confirm that this is the case, and further show that the divergence of ∇ * q * reduces to a velocity Laplacian The other side of our non-commutative equation is represented by the upper-convected derivative of Q * . From our definitions above clearly Q * = ∇ * · I * = 0, and thus the upper-convected derivative of this quantity must also be zero. This shows conclusively that the two operators are not commutative.
We start the correct analysis with governing equations (1) and (2), which are made dimensionless via the introduction of the following variables where δ = (ρ * U * L * /µ * 0 ) −1/2 = Re −1/2 , is the standard boundary layer length scale, and L * and U * are characteristic length and velocity scales, respectively. This choice of scaling leads to the definition of two viscoelastic dimensionless parameters where Wi, the ratio of relaxation time and the flow time scale, is the Weissenberg number and 1−β is the ratio of retardation and relaxation times respectively. We note that, under the definitions µ * with µ * s and µ * p the solvent and polymer viscosities respectively, β also represents a ratio of viscosities, which must remain bounded: 0 ≤ β ≤ 1.
In the first instance we will consider an Oldroyd-B (OB) analysis, this being the case when Wi = 0, and β = 1. In this instance, given the definitions of the dimensionless variables, it is relatively straightforward to show that the boundary layer flow is governed by the following system of equations where the notation ∂ k denotes the partial derivative with respect to k. At this stage of our investigation, we restrict the OB model to, at largest, moderate values of the Weissenberg number i.e., Wi ∼ O(δ 0 ), and we note that different scalings are introduced if Wi is assumed to be large. All that remains now is to determine the relevant boundary layer scalings for the three stress components (T xx , T xy , T yy ). In the case when Wi → 0, we must return the standard Newtonian boundary layer equations. As such, it follows immediately that the correct scaling for the shear stress component is, T xy = δ −1 τ xy . Given this result, when considering (3f), we see that in order to ensure the correct leading order balance it must follow that T yy = δ 0 τ yy . Lastly, from (3d), we observe that a non-zero shear stress will only be predicted in the case when Having now determined all the correct scales we replace T ij with τ ij throughout (3b)-(3f), and take the limit as Re → ∞. We then arrive at the OB boundary layer equations noting that to leading order the pressure is a function of x only; this follows immediately from (3c). In order to remove the pressure from the problem we consider the flow in the free-stream where the Figure 1: Schematic diagram of the stagnation point boundary layer flow over a wedge inclined at angle of π/2 from the horizontal. Sketches of the streamlines are indicated in blue. A typical boundary layer velocity profile with free-stream velocity, U * ∞ , is highlighted in black. All stated quantities are dimensional.
streamwise velocity is a function of x only, and the wall normal velocity is identically zero. Outside of the boundary layer it must also be true that the stress varies only in the streamwise direction. Therefore, in the limit as y → ∞, we have from (4b) that where the superscript F denotes a free-stream quantity. As a first attempt to solve system (4) we seek a self-similar solution of the form where, at this stage, we do not specify any restrictions on the three stress components. From (4b) we have that where the primes indicate differentiation with respect to the similarity variable η. In order to achieve a self-similar solution, we must guarantee that all the terms in above equation have the same x dependence. Given the form of the left-hand side of the (5), this concisely forces the definition of the free-stream velocity, such that, u F = x m , where m dictates the angle of inclination of the flat plate as noted at the outset of this section. In turn, this implies that τ xy = x (3m−1)/2 s xy (η), and that τ xx = x 2m s xx (η). Therefore, upon inserting these three functions of x into (5), we arrive at It must then follow that τ F xx = Cx 2m , where the constant C is determined from the solution of (4c) in the limit as y → ∞. One finds that, in this limit, a general solution of (4c) can only be found in the instance when C is identically zero and we must then have that s xx → 0, as η → ∞.
In order to determine the range of values of m for which self-similar solutions can exist, we must investigate the form of the governing stress equations. Analysis of (4c)-(4e) reveals that τ yy = x m−1 s yy (η), and that x independence can only be ensured in the specific case when m = 1, i.e., when one considers stagnation point flow over a wedge inclined at an angle of π/2 from the horizontal. We provide a schematic diagram of this flow configuration in Figure 1. In this specific case we have that u = xf (y), v = −f (y), τ xx = x 2 s xx (y), τ xy = xs xy (y), and that τ yy = s yy (y). We note that primes now indicate differentiation with respect to y.
In order to simplify the forthcoming analysis it proves useful to introduce the following transformations T xx = s xx , T xy = s xy − f , and T yy = s yy + 2f . Having done so, we are able to express the problem as a system of six coupled first order ordinary differential equations This system must be closed subject to six boundary conditions. Given that the surface of the plate is impermeable we must have that v(y = 0) = 0, thus f (y = 0) = 0. The no-slip condition at the surface implies that u(y = 0) = 0, and hence g(y = 0) = 0. Given that the boundary layer flow must match with the free-stream above, we have the condition u(y → ∞) → u F , which implies that g(y → ∞) → 1.
Our final three conditions must relate to the three stress components.
We have already shown that T xx (y → ∞) → 0, it remains, therefore, for us to arrive at physically relevant conditions for the functions T xy and T yy . We note that as we approach the limit of infinite distance from the wall, both τ xy and τ yy must tend to functions of x, the variable of flow direction, and be independent of y. This is not to say that the functions must be zero, in fact we will go on to show that this is not the case. However, we must have no change in the stress contributions with respect to y, otherwise there will be a discontinuity in our solutions at the free-stream. It follows that our final two boundary conditions are then s xy (y → ∞) → 0, and s yy (y → ∞) → 0. In summary, we must solve system (6) with respect to the following six boundary conditions where the infinity subscript denotes evaluation in the limit as y → ∞. In order to determine the values for the constants h ∞ , and h ∞ , we look to the decay of the solutions in the far-field.
Given that the streamwise velocity function g tends to unity as y → ∞, then f → (y − δ 1 ), where is the displacement thickness. It proves useful, in what follows, to introduce the shifted coordinate Y = y − δ 1 . Given the far-field conditions on T xx , and T xy , we have from (6c) that the function h must satisfy in the limit as y → ∞. This ODE has the solution where c is a constant of integration and D is Dawson's Integral function. Then, from (6f) it immediately follows that as y → ∞. In order to ensure matching with the free-stream it must then be true that h ∞ = 0: if this is not the case then the solution for T yy grows linearly with respect to y at the outer edge of the boundary layer. Given this form for the solution for T yy , in the limit as y → ∞, from (6e), one can show that where a n = nWi (1 + nWi) .
Again, in order to ensure matching with the free-stream it must then be true that h ∞ = 0, if this is not the case then the absolute value of the solution for T xy grows linearly with respect to y at the outer edge of the boundary layer (note that this conclusion could also have been drawn from analysis of the expression for T yy ). Thus, the function h decays to zero exponentially, h → c e −(Y 2 −δ 2 1 )/2 , irrespective of the value of the Weissenberg number. Therefore, as y → ∞, the three stress functions decay to their respective constant values as follows The result for the function T xx follows directly from integrating (6d) with the the known large-y expression for T xy . This ODE has a decaying solution that can be expressed in terms of the function h, the upper incomplete gamma function and the complementary error function. To leading order, the terms involving upper incomplete gamma function and the complementary error function can be approximated in terms of powers and h and Y , and the result above follows. In order to ensure some level of brevity, the details of this calculation are included for the interested reader in the §B.
Given the above analysis, we are now in a position to solve the system of six coupled first order ODEs that govern OB stagnation point flow, (6), subject to the six physically correct boundary conditions

Numerical Method and Results
We solve (6) subject to (8) using a spectral method, in particular a variation on the Chebyshev collocation scheme. We split the one dimensional domain, y = [0, y ∞ ], into N + 1 points, including the surface of the flat plate and the far-field location, which we define to be at y 0 = 0 and y N = y ∞ , respectively. At this stage, we choose not to associate a specific value with y ∞ , instead we keep the formulation general and determine a value for y ∞ based on the outcome of our numerical testing procedure. The solution points are linearly related to the Chebyshev collocation points χ n = cos (πn/N ), with n = 0, 1, ..., N , such that y n = y ∞ (1 − χ n ) /2, covers the entire domain. The majority of the interesting dynamics associated with this problem occur in a close proximity to the flat plate, which can be better approximated by these collocation points as opposed to a uniform grid, given the higher concentration of points at both ends of the y-domain. We approximate all the flow quantities, q, as a finite sum where denotes that the first and last terms of the sum are halved, a j represents the set of constants to be fitted from our discrete solutions, and T * j is related to the jth Chebyshev polynomial of the first kind T j , applied at location χ n , via This form for the yet unknown functions is applied to all the quantities present in both the governing equations (6), and boundary conditions (8), including the stress components and their derivatives. The one exemption being h , which is obtained post-hoc. The polynomials themselves then form a set of orthonormal basis functions which are widely used in interpolation and optimisation problems, not least due to their property of having a maximum magnitude of 1 in the range χ ∈ [−1, 1]. One can define these functions using their orthogonality condition with respect to the inner product in integral form, but a simpler definition follows directly from the recurrence relation Using the discrete orthogonality relation we can then calculate the unknown constants where q(y n ) is the value of a flow quantity at a point in the discrete domain. These definitions allow us to return a continuous function for any quantity q, which is known exactly at the N + 1 distinct points. Furthermore, we may find the derivative of a generic quantity, q , which is described by way of a differentiation matrix which satisfies q (y i ) = N j=0 A ij q(y j ), where T * n represents the derivative of Chebyshev polynomial T * n with respect to y. This statement holds true for general choice of flow quantity, which is particularly useful when q itself is a derivative. The governing equations (6) are converted into a set of 9 closed equations, taking care to account for the number of boundary conditions noted in (8). The first 4 equations represent the conservation of momentum and constitutive stress relations. These are implemented at every point in the y-domain, which therefore provides 4N + 4 equations to solve. Where a flow quantity, for example, g, is known at either boundary due to our imposed conditions, this value is implemented directly in the numerical scheme and is not solved for. We also note that the only function that is not directly determined as part of our numerical scheme is h . Therefore, it is necessary to make explicit use of the differentiation matrix A ij , in the determination of this function only.
The remaining governing equations describe the first order derivatives of the 5 primary flow functions f, g, T xx , T xy , and T yy . We do so by making use of the differentiation matrix as described overleaf, solving, respectively, for where δ iN represents the Kronecker delta. The y i locations at which we solve the first order system above have been chosen to remove degrees of freedom equal to the number of boundary conditions, such that we have a total of 9N + 3 equations to solve numerically. We choose to remove either one or two equations from the system based on whether there is a fixed boundary condition on a particular quantity or its integral. The first study involves a relation where the quantity summed over has at least one boundary condition attached to it. Take, for example, the function h(y i ), as described above. We evaluate this function at all solution points apart from at the flat plate (y = 0) and the far-field location (y = y ∞ ), since we constrain its integral g(y i ) at these points. The other case is considered in, for example, the calculation of T xy (y i ), where the boundary condition is applied to the quantity on the left-hand side. In this case, since no boundary conditions are placed on T xy , we are required to solve for it at all points y n . However, since we know the value of T xy (y = y ∞ ), the location i = N is excluded from the calculations. This is equivalent to the statement that we know the value of q at a certain location, so we should not solve for it. We note however that this is not true of the first equation for g(y i ) above, which does constrain f with the known condition g(y ∞ ) → 1. We justify this by recognising that f does have an attributed boundary condition, and so falls into the first category of equation removals described above. We use MATLAB to solve the system as described above, with the Levenberg-Marquardt algorithm up to a norm mean square error of 10 −18 . The numerical solutions presented here are calculated with a value of N + 1 = 200 solution points for a variety of both Wi and β values.
In order to determine a suitable value for the location of the far-field boundary y ∞ , we solved the Newtonian flow problem (Wi = 0) for a range of combinations of both N and y ∞ . In each case we maintain the norm mean square error tolerance of 10 −18 . These Newtonian results were also compared to solutions obtained from a fourth-order Runge-Kutta integration scheme twinned with a secant shooting method. The difference between the results obtained via these schemes was numerically indistinguishable. Moving forward, we set y ∞ = 5 to be the free-stream location, which captures the dynamics in the farfield in precisely the same way as any larger value for each of our viscoelastic models.
In the first instance, we will consider solutions to (6) subject to (8) for an arbitrary fixed value of the dimensionless retardation parameter, β = 0.8. We note that the special case when β = 1, corresponding to a system with an upper-convected Maxwell constitutive viscosity law, will be covered, in detail, in §4.
friction coefficient is defined as follows: Given that f (0) = g(0) = s yy (0) = 0, neither Wi, nor β appear explicitly in the expression for the Oldroyd-B skin friction coefficient. Therefore, C f = 2x Re −1/2 h(0), which is identical to the standard Newtonian expression. The other important quantity for flows such as these is the first normal stress difference at the wall: We observe that the first normal stress difference at the wall is directly proportional to the square of the skin friction coefficient,

Small Wi ∼ O(1)
In the first instance we consider low Weissenberg number flows, which is in keeping with previous analyses available in literature. A comparison between our solutions and those of preceding studies can be found in §4. It is evident, from the results presented in Figure 2 (a), that as the value of the Weissenberg number is increased, so does the streamwise velocity within the boundary layer. Similarly, as evidenced in Table  1 and Figure 3, the value the skin friction coefficient, and therefore also the first normal stress difference at the wall, increases with increasing Wi. This is an indication that the effects of viscoelasticity act to thin the boundary layer profile. Indeed, the fluid particles accelerate near the wall and obtain the value of the free-stream velocity closer to the surface of the plate, primarily because of the increased difference in the first normal stress difference at the wall, with increasing Weissenberg number. This correlates directly to a decreasing of the displacement thickness, δ 1 , as the value of Wi increases.
From Figure 2 (b), (c), and (d) we are able to correlate our approximate predictions for the decay of stress functions, in the limit as y → ∞, with our numerical solutions. We observe, as predicted by (7), the rapid decay of the streamwise stress function τ xx /x 2 , when compared to both the shear stress function τ xy /x, and the wall normal stress function τ yy . For unconfined boundary layer flows, such as the flows analysed here, the shear rate is largest at the wall. As one moves away from the wall the shear rate dissipates and the free-stream conditions are attained. The introduction, therefore, of fluid elasticity generates, at the wall, a non-zero stress contribution parallel to the direction of the flow. This explains the form of the τ xx /x 2 profiles for increasing values of the Weissenberg number. In a similar fashion, as we transition from analysing a purely viscous fluid (Wi = 0), to a fluid with elastic properties (Wi > 0), the shear stress at the wall increases as a result of the introduction of elastic forces, hence the form of the profiles observed in Figure 2 (c). We note that the boundary conditions imply that the wall normal stress function, τ yy , will always be zero at the wall. From (4e) we see that τ yy is proportional to ∂v/∂y = −f = −u/x, it must, therefore, be the case that this function attains a constant negative value at the outer edge of the boundary layer. From Figure 2  We note that, although not plotted here, as the value of β decreases from 1, for a fixed value of the Weissenberg number, the thickness of the boundary layer increases, and the values of both the skin friction coefficient and first normal stress difference at the wall decrease. Physically, these results match with ones' intuition. Decreasing the value of β effectively reduces the ability of the fluid to overcome retardation effects and thus, for a fixed free-stream velocity, one would expect to observe a thickening of the boundary layer.

Moderate Wi ∼ O(1)
The spectral method we use to calculate the flow profiles is not limited to relatively low values of Weissenberg number, as is frequently considered in the literature. We find that, given the number of collocation points as prescribed, the system can be well described up to a value of Wi ≈ 3. This value falls well within the O(1) restriction which was imposed in the non-dimensionalisation procedure outlined in §2. Figure 4 shows the dependence of our horizontal velocity and stress components, given a variety of moderate O(1) Weissenberg numbers and a constant value of the retardation parameter, β = 0.5. We can see that, similarly to the smaller Wi cases, an increase in the Weissenberg number yields a return to free-stream profiles closer to the wall. As noted previously, an increase in viscoelastic effects within the fluid phase therefore leads to a decrease in the boundary layer thickness. Figure 4 (a) shows a seeming convergence to fixed flow profile as Wi increases. This likely indicates that there is a threshold for Wi at which the viscoelastic contributions to the governing equations become dominant. We note that the dependence is clearly non-linear, and we confirm that intermediate Weissenberg number velocities converge in a monotonic fashion. The investigation into Wi ∼ O(δ −1 ) simulations, given the same scalings provided in §2, is left as an open problem.
Stress functions τ xx , τ xy and τ yy are plotted in Figure 4 (b), (c), and (d), respectively. Qualitatively similar results are obtained to those presented in § 3.1. The main difference between the two sets of results is, however, the magnitude of effects. For example, we observe that the stress function τ xx /x 2 attains a much higher peak at the wall when compared to Figure 2 (b). Furthermore, we note that the   free-stream value in all the stress components is attained closer to the wall, which mimics the behaviour of the streamwise velocity.
For fixed Wi = 2, we evaluate the variation of the streamwise velocity profile u/x as the retardation parameter β decreases. We observe in Figure 5 (a) that as β reduces in value, which represents a decrease in polymer viscosity (compared to solvent viscosity), the streamwise velocity decreases accordingly. This results in a thickening of the boundary layer at fixed moderate values of the Weissenberg number.
In Figure 5 (b) we plot the variation of the displacement thickness, δ 1 , for a range of moderate Weissenberg numbers and our familiar choices of the retardation parameter. Naturally, the same behaviour as described above in Figure 5 (a) is confirmed here: the boundary layer thickens as the value of β decreases. Furthermore, we observe that δ 1 appears to approach a fixed constant value as Wi increases towards the upper end of the spectrum of moderate values. We note that this result is not immediately evident from an analysis of the governing equations. Independent of the value of β, all curves show a monotonic decreasing nature, which has strongest gradient at the Newtonian limit. They also show clearly that any introduction of viscoelastic effects, or indeed a polymer viscosity will produce a thinning of the boundary layer, when compared to the Newtonian constant value δ 1 = 0.6479. Indeed, when β = 0.7, and Wi = 3, we observe an approximate halving the thickness of the boundary layer when compared to the Newtonian counterpart.

Comparisons With Previous Studies
As noted in §1, numerous previous studies have considered the stagnation point boundary layer flow problem for viscoelastic fluids, with particular emphasis being focused on the upper-convected Maxwell model (β = 1). To the best of our knowledge, all previous studies are confined to the cases when Wi ≤ 1/2. For this reason, we restrict comparisons to our small Wi case study. We note that, in contrast to this high Reynolds number study, in the limit of creeping flow (Re 1), many previous investigation have not been limited to this restrictive range of values of the Weissenberg number and, instead, choose to scale Wi with some inverse power of the boundary layer thickness [31], [32], [33].

Upper-convected Maxwell fluid model
In the case when β = 1, the governing ODEs for this problem, (6), reduce to those representing an upper-convected Maxwell constitutive viscosity law. In this instance we are able to compare our results with those of [15] who, following the analysis of [16], arrived at the following ODE that is claimed to model the stagnation point flow of a fluid exhibiting an upper-convected Maxwell viscosity where k is referred to as the elasticity number. Given our notation, we will consider this to be equivalent to our definition for the Weissenberg number. The authors solve this ODE subject to the standard boundary layer conditions noted previously f (y = 0) = f (y = 0) = 0, f (y → ∞) → 1.
Given the relatively simple form of this third order ODE one could use a variety of different numerical methods to arrive at accurate solutions. Indeed, [15] (from this point onwards referred to as SHT), opted to use a Chebyshev spectral scheme, somewhat similar in nature to our own scheme described in §3. SHT produce results for Wi in the range Wi = [0, 0.3] and analyse both the streamwise velocity profile and the displacement thickness. We reproduce their results using a simplified version of our own spectral scheme. Qualitatively, our reproduction of SHT's results appears to be exact. However, we are unable to make a quantitative statement due to the fact that SHT did not report any numerical values from which we could benchmark.
In Figure 6 we compare our results with those of SHT for both of these quantities. We observe that as the value of the dimensionless viscoelastic parameter, Wi, increases from 0 (Newtonian flow) our results differ in every respect when compared to those of SHT. As the Weissenberg number increases from zero our results indicate that the free steam velocity will be attained closer to the surface of the flat plate. This result is a direct consequence of the predicted increase of shear stress at the wall (see Table 1). Naturally, this is reflected in a thinning of the boundary layer and a reduction in the value of the displacement thickness δ 1 . On the other hand, SHT's analysis predicts that the shear stress at the wall will decrease and the boundary layer itself will become broader when compared to its Newtonian counterpart. These results are contrary not only to those owing from our analysis of the full governing   [7] is represented by the dash-dotted black curve, and that of [15] is represented by the dashed-line curve.
equations but also to those of [7] who predict, in the limit when Wi 1, that both the velocity in the boundary layer, and stress at the solid boundary, will increase because of the effect of viscoelasticity.
In Figure 7 we compare our solution for the streamwise velocity component with those of SHT and also [7] (from this point onwards referred to as BW). In order to reproduce the results of BW we solve their governing equations ((23) -(24)) subject to the relevant boundary conditions ((25) -(26)) using a fourth-order Runge-Kutta (RK4) integrator alongside a Newton-Raphson shooting routine. In the case when Wi = 0.05, BW predict that τ xy (0) = 1.2895. When we attempt to reproduce their analysis we find that τ xy (0) = 1.2896. The relative closeness of these results suggest that our RK4 scheme correctly reproduces the results of BW (one would not necessarily expect the two values to match identically given that we are able to produce results to a much higher numerical tolerance). Given that BW's analysis is only strictly valid for small values of the Weissenberg number, and in this case we have set Wi = 0.05, it is reaffirming that their result match almost identically with our solution obtained from the governing equations (6) and boundary conditions (8). In contrast, SHT's result does not closely match that of BW's, something one would expect given the relative smallness of the Weissenberg number.
SHT argue that the difference between their solutions and those of BW are because BW are essentially analysing a second-order fluid model problem. Whereas, they claim to be solving the equations relevant to a UCM boundary-layer flows. However, as we have shown, for flows of this nature, the UCM model collapses down to the second-order fluid model in instances when the Weissenberg number is small. Thus, one would not expect there to be a vast disparity between the results of the two studies for the range of values of the Weissenberg number that were analysed Wi = [0, 0.3]. Indeed, SHT note that their 'work serves to demonstrate quite clearly that the constitutive equation of a fluid may be of crucial importance in boundary layer studies of viscoelastic fluids'. However, having followed the incorrect formulation presented by [16], SHT have, in fact, arrived at erroneous conclusions that fail to capture the correct physics of the problem.
Aside from the work of SHT, we can compare our results to those reported by [21]. The governing equations used in [21] cannot be reduced to a single governing ODE, as they have been in the case of SHT's analysis, and it is a relatively straightforward task to show that Phan-Thien's system of ODEs is directly equivalent to those presented in §3. We note, however, that there is a discrepancy between the boundary conditions, for the three stress functions, that [21] employs at the far-field, and our own. Due to a slightly differing formulation to the problem, Phan-Thien's analysis is also restricted to flows whereby the Weissenberg number must satisfy the condition that Wi < 1/2. The only reported values in that study, to which we can make any comparison, are those for the thickness of the boundary layer. This quantity is defined in [21] to be the point at which the streamwise velocity function attains 99% of the free-stream velocity (from this point onwards we will denote this constant value ∆ 99 ). We find that, to the same order of accuracy as given in [21] (2 d.p.), the results owing from our spectral approach match precisely with those obtained from Phan-Thien's central difference scheme.
In an attempt to extend the study of SHT, [23] return to the stagnation point flow problem, and consider a Giesekus governing viscosity law. This model is only slightly more complex than the upperconvected Maxwell model, and it should be noted that under the correct limiting choice of constants, the Giesekus model reduces in complexity to UCM. However, upon comparison of both our own UCM results (and those of [21]) with those presented by [23], we note a discrepancy between the reported values for the boundary layer thickness, ∆ 99 . Again, this is the only tangable reported value that we can make any comparison to. Having analysed the formulation presented in [23] with our own (and Phan-Thien's) it would appear that the reason for this discrepancy is due to the relative accuracy of their numerical scheme. [23] utilise a Keller box method twinned with a Newton linearisation procedure. Very few details about this linearisation process are forthcoming in the Mirzadeh and Sadeghy's article and it would seem likely that it is this numerical procedure that causes the disparity between the results, as opposed to the aforementioned work of SHT, where clearly the erroneous governing equations are the direct cause of the incorrect solutions.

Oldroyd-B fluid model
In the more general case when β = 1, the Oldroyd-B fluid equations do not reduce to an explicit upperconvected Maxwell model. We can therefore compare our results to those of [22] (referred to from now on as P-T), who consider a similar, though not identical, stagnation point flow over a flat plate. As previously mentioned in §1, due to P-T's formulation of the problem he is restricted to considering only a small finite range of the Weissenberg number, Wi = [0, 0.5]. Indeed, P-T chooses not to present any results, either tabulated or graphical, for the case when Wi > 0.4. This is in contrast to our study where no such restriction applies. [22] uses an equivalent set of governing equations as (6), the only difference being a reformulation of the governing viscosity equation in terms of τ * N N = τ * − 2E * , the non-Newtonian contribution to the stress. This difference in formulation is expected to be of no consequence. However, upon comparing results for the boundary layer thickness ∆ 99 , there is a clear difference between our results and those of P-T. We provide our calculated values in Table 2 to an accuracy of 5 significant figures, along with the literature values (3 s.f.) as well as the relative difference between the results for the range of Wi and β values quoted by P-T.
We observe that there is a marked increase in the difference between the results as β decreases from its UCM value, β = 1. At this specific value, we match the literature ( [21] & [22]) precisely to the given  Table 2: A comparison between our results for the boundary layer thickness, ∆ 99 , and those presented by [22]. The relative difference between the solutions is noted in each of the third sub-columns.
accuracy of that paper, which cannot be said in any of the cases when β = 1. At its largest, we observed a relative difference between our results and those of P-T of 1.75%. We note that we do not see a uniform change in the percentage error for increasing Wi or decreasing β, as one might expect. This is attributed to the relative accuracy of P-T's stated results. Having said that, in general, we observe that for larger values of both the Weissenberg number and the retardation parameter the error between the two solution sets grows. The disparity between these two sets of results cannot be due to P-T's differing formulation of governing equations: we have used our in-house spectral code to model exactly P-T's system of ODEs. We find that the values of ∆ 99 from these calculations match precisely with our own results tabulated in Table 2, i.e., when we mimic P-T's analysis we obtain our results. The finite difference scheme employed by P-T is based on a central difference approach, and provides accuracy up to 3 s.f.. However, this cannot account for the source for difference between our results and those of P-T, since even the EG values to the same order of accuracy do match those of P-T. We must conclude, therefore, that there is some small error in the code used by P-T.

Conclusions
In this paper we have investigated viscoelastic boundary layer flows of fluids described using one of three different models. We have derived the full set of coupled equations in two dimensions which represent all of these fluid models, and further rigorously determined the boundary conditions that must be imposed at the far-field. In the process of doing so, we have highlighted the inadequacies of a number of previous studies in the literature, which fail to capture the correct physics of the problem. The coupled system of governing equations was solved using a Chebyshev collocation method, to a high degree of accuracy.
Results for the streamwise velocity profile and the individual stress components are presented for specific choices of viscoelastic and viscosity parameters. We find that, contrary to the findings of previous investigations, in particular [15], an increase in viscoelasticity results in a thinning of the boundary layer. Our results tend to agree qualitatively with the work of [21], [22] for lower Weissenberg numbers. However, with reference to the determination of the boundary layer thickness, there is a clear disparity between the specific values quoted by these studies and our own. We provide context for the variation of the three stress components across the boundary layer, and show that all of them behave well in the Newtonian limit, where the solution is well known. Furthermore, we evaluate both the skin friction coefficient and first normal stress difference at the wall for fluids with viscosity captured by the Oldroyd-B model. Both of these quantities show a monotonic increase with increasing viscoelastic effects, with the same statement being true as the ratio of the solvent to total viscosity increases.
For moderate values of the Weissenberg number, we present results that retain the high degree of accuracy associated with our spectral method, which would typically be very challenging for finite difference schemes to reproduce. The streamwise flow velocity, along with the stress functions can be seen to approach invariant profiles with increasing values of the Weissenberg number. This behaviour is likely indicative of a transition change from the regime where Wi ∼ O(1) to higher orders. We note that at a value of Wi = 3, with β = 0.7, the displacement thickness approximately halves when compared to the Newtonian reference.
Given the relatively general nature of this work, there is a reasonably large range of natural extensions that one may wish to consider. Firstly, the framework presented here could be extended to encapsulate other viscoelastic models. One could consider, for example, extensions to the Phan-Thien Tanner, Giesekus or FENE-P models.
Conversely, one may choose to analyse different types of boundary layer flows whilst considering only the viscoelastic models discussed here. One obvious candidate, given its relative prominence in the literature, would be the boundary layer flow induced by the stretching of a solid surface. In effect, a study of this nature would serve as a correction to the analysis presented by [14].
At the outset of this study we restricted our attention to moderate Weissenberg number flows. However, there is no reason why the methodology presented here could not be reformatted to solve boundary layer flows associated with large values of the Weissenberg number. Indeed, the governing equations for that problem, as noted by [33], are simply a modified subset of the governing boundary layer equations presented here.
With regards to our original motivation to tackle this problem, we intend to use the results of this study to form the basis for a much wider body of work focusing on the injection of a non-Newtonian fluid in to an otherwise Newtonian boundary layer flow. We expect to be able to report on the results of this study in the relatively near future.
with the easiest comparison between the two coming from choice of constants α * 2 = 0, −κ * 0 = α * 1 . The apparent missing dynamics, which originate from the quadratic strain rate tensor term, can be shown to play no role in the solution of velocity or extra stress contributions for flat plate boundary layer flows. We investigate this idea first by assuming a simplified version of the stress tensor: τ * = 2µ * 0 E * + 4α * 2 (E * · E * ), so chosen to capture the minimum amount of required physics. The first term is the Newtonian contribution, which we know must be dominant in the limit of zero viscoelasticity [1], and the second is our non-Newtonian term of interest. This stress profile can be inserted directly into the momentum equations, reducing the number of governing equations from 6 to 3. Under the same non-dimensionalisation and scaling process as outlined in §2, with the definition of a new viscoelastic parameter k = α * 2 U * /µ * 0 L * , we derive the following set of coupled differential equations: ∂ x u + ∂ y v = 0, u∂ x u + v∂ y u = −∂ x p + δ 2 ∂ xx u + ∂ yy u − 2k ∂ y u∂ xy u + δ 2 ∂ x v∂ xy u + 4δ 2 ∂ x u∂ xx u + δ 2 ∂ y u∂ xx v + δ 4 ∂ x v∂ xx v , u∂ x v + v∂ y v = −δ −2 ∂ y p + δ 2 ∂ xx v + ∂ yy v − 2k δ −2 ∂ y u∂ yy u − 4∂ x u∂ yy v + ∂ x v∂ yy u − ∂ y u∂ xx u − δ 2 ∂ x v∂ xx u .
Here we note that the momentum equations above have been simplified via use of the continuity equation.
In a similar fashion to our previous analysis, we solve the y momentum equation for pressure, leading to the succinct formula p(x, y) = P (x) − k[∂ y u(x, y)] 2 , where P is a yet unknown function of x. One can see that upon inserting this definition into the x momentum equation, all terms dependent on the viscoelastic parameter k will naturally cancel out: u∂ x u + v∂ y u = −∂ x P + ∂ yy u.
Since we know k to be a factor of every term in our boundary layer equations which comes from the E * · E * contribution, it is clear that both velocities u * , and stresses τ * , act independently of it. Instead, we find that the quadratic rate of strain tensor contributes to the pressure function only, in the form of the Newtonian shear derivative squared. We also note that this result is true independent of the angle of inclination of our boundary layer flow.

B Decay of the function T xx in the far-field
Using the results for the far-field decay of the functions T xy and T yy , presented in §2, one can then determine the form of the decay of the function T xx , to zero, as y → ∞. From analysis of (6d). The relevant ODE to solve is then where b = 1 + a 1 (2a 2 − 3). This differential equation has the decaying solution , where Γ is the upper incomplete gamma function, erfc is the complementary error function, and the constant of integration has been chosen such that the free-stream boundary condition is satisfied. Now, to leading order