Viscous flows in corner regions: Singularities and hidden eigensolutions by James Sprittles | Papers by James

Viscous flows in corner regions: Singularities and hidden eigensolutions James E. Sprittles and Yulii D. Shikhmurzaev June 18, 2009 arXiv:0906.3462v1 [physics.flu-dyn] 18 Jun 2009 Abstract Numerical issues arising in computations of viscous flows in corners formed by a liquidfluid free surface and a solid boundary are considered. It is shown that on the solid a Dirichlet boundary condition, which removes multivaluedness of velocity in the ‘moving contact-line problem’ and gives rise to a logarithmic singularity of pressure, requires a certain modification of the standard finite-element method. This modification appears to be insufficient above a certain critical value of the corner angle where the numerical solution becomes mesh-dependent. As shown, this is due to an eigensolution, which exists for all angles and becomes dominant for the supercritical ones. A method of incorporating the eigensolution into the numerical method is described that makes numerical results mesh-independent again. Some implications of the unavoidable finiteness of the mesh size in practical applications of the finite-element method in the context of the present problem are discussed. 1 Introduction The ability of a numerical scheme to accurately approximate a physical problem in a domain containing corners is critical for the description of a number of phenomena, ranging from electromagnetic wave propagation in a waveguide [1] to die-swell effects in polymer extrusion [2]. Often, an analysis of such problems reveals singular behaviour of variables as the corner is approached, which requires special numerical treatment. Such problems are well known and have been thoroughly investigated in the setting of fracture mechanics, where one considers the propagation of a crack into a material [3], and have also been studied in some fluid dynamics problems [4]. Our interest here is the viscous flow in a corner formed between a liquid-fluid free surface and a solid boundary. Although a free surface is generally bent, to leading order as the corner is approached, it is often possible for the purpose of a local analysis to consider the flow domain as having a wedge shape. This is the case, for example, in dynamic wetting flows [5], where the free surface and the solid boundary form what is referred to as the ‘contact angle’ and the liquid-fluid-solid ‘contact line’ moves with respect to the solid surface. This differs from the situation considered in some previous investigations on flow in a corner where the contact line is stationary with respect to the solid and motion is generated by disturbances in the far field [6, 2]. It is well known that the classical fluid-mechanical approach when applied to dynamic wetting problems fails to provide an adequate description of the flow [7]. The conventional remedy to the problem is to relax the no-slip boundary condition on the solid surface and allow for ‘slip’ between the liquid and solid. A number of different forms for this slip behaviour have been examined in the literature (for a recent review see Ch. 3 of [5]). Broadly, these split into 1 conditions which (i) relate tangential stress to the slip velocity, such as the Navier condition [8], or (ii) explicitly prescribe the velocity along the solid surface. In this paper, we consider numerical problems arising in the second case where the velocity along the solid surface is a priori prescribed in a form which ensures that a solution exists, and that in the far field the usual no-slip condition is restored. This approach is appealing to some users due to its mathematical simplicity, and our goal here is to show numerical pitfalls one comes across in its numerical implementation and give a method of overcoming them which provides a framework for modelling this class of problems. A number of functions prescribing the fluid velocity on the solid surface have been proposed in the literature [9, 10, 11] and here we consider just one of these which captures all the main features of the problem. 2 Problem formulation The problem is most easily formulated in a polar coordinate system (r, θ) in a frame moving with the contact line (now referred to in our two-dimensional domain as the corner point). The wedge is formed by a solid surface at θ = 0 which moves at speed U parallel to itself, a flat free surface at θ = α and a ‘far field’ boundary which is placed at an arc of a sufficiently large radius r = R. The liquid is Newtonian and incompressible, with density ρ and viscosity µ. Near the corner the flow is characterized by a small length scale so that the Reynolds number Re based on this scale is small. Then as Re → 0, to leading order in Re we have the Stokes flow1 . The nondimensional Stokes equations for the bulk pressure p and the radial and azimuthal components of velocity (u, v) take the form: 1 ∂(ru) 1 ∂v + = 0, r ∂r r ∂θ ∂p u 2 ∂v = ∆u − 2 − 2 , ∂r r r ∂θ where ∆= (0 < r < R, 0 < θ < α), 1 ∂p v 2 ∂u = ∆v − 2 + 2 , r ∂θ r r ∂θ (1) (2) ∂2 1 ∂ 1 ∂2 + + 2 2. ∂r2 r ∂r r ∂θ On the solid surface, for a solution not to have multivalued velocity at the corner point [12, 5], we replace the no-slip condition (u = 1) with a prescribed velocity that has free-slip at the corner point and attains no-slip in the far field, that is: u = 0, at r = 0 and u → 1, as r → ∞. (3) Following [11], we use an exponential form for this function and, to complete the boundary conditions on the solid surface, it is combined with the usual impermeability condition for the component of velocity normal to the surface: u = 1 − exp(−r/s), v = 0, (0 < r < R, θ = 0). (4) The region in which the velocity deviates from no-slip is characterized by the value of s, which is a (non-dimensional) ‘slip length’. On the free surface, we have the standard boundary conditions of zero tangential stress and impermeability: ∂u = 0, v = 0, (0 < r < R, θ = α). (5) ∂θ The analysis remains valid for the full Navier-Stokes problem since it considers the r → 0 limit, and here we consider the Stokes equations in the whole region as a convenient way to illustrate the idea. 1 2 In the far field, we assume that the flow is fully developed and apply ‘soft’ conditions: ∂u ∂v = = 0, ∂r ∂r (r = R, 0 < θ < α), (6) which imply that the influence of slip has attenuated; these conditions are satisfied by the (multivalued at the corner point) solution obtained using the no-slip condition all along the solid surface [13]. Equations (1)–(6) fully specify the problem of interest. 3 Local asymptotics Consider the leading-order asymptotics for the solution of (1)–(6) as r → 0 [5] that we will later need to use in the numerical code and to provide a test of accuracy of the numerical results presented in the next section. After introducing the stream function ψ by u= 1 ∂ψ , r ∂θ v=− ∂ψ , ∂r (7) equations (1)–(2) are reduced to a biharmonic equation ∆2 ψ = 0 with boundary conditions (4)–(5) taking the form ∂ψ = r (1 − exp (−r/s)) , ∂θ ψ = 0, (θ = 0, 0 < r < R), (8) ∂ 2ψ = 0, ψ = 0, (θ = α, 0 < r < R), (9) ∂θ2 where, for definiteness, we assign the value zero to the streamline coinciding with the wedge’s boundary. Condition (8) is the only inhomogeneous boundary condition in the problem, i.e. the condition that drives the flow. To leading order as r → 0, it has the form ∂ψ ∂θ = ar2 + O r3 , θ=0 (10) where a = 1/s. An alternative prescribed velocity that satisfies (3) and (10) known in the literature [9] is given by u = (r/s)/(1 + r/s) and the asymptotic analyses throughout this paper are equally valid for this function as well. The form (10) suggests looking for the leadingorder term of the local asymptotics in the form ψ = r2 F (θ), which is a particular case from a known family of separable solutions of the biharmonic equation of the form ψ = rλ F (θ). After substituting ψ = r2 F (θ) into ∆2 ψ = 0, one arrives at ψ = r2 (B1 + B2 θ + B3 sin 2θ + B4 cos 2θ) , where the constants of integration Bi (i = 1, . . . , 4), found from (8)–(9), are given by2 : B1 = −B4 = aα sin 2α , 2α cos 2α − sin 2α B2 = −B1 /α, B3 = B1 cot 2α. (12) (11) The pressure field obtained from (2) using (7) and (11) has the form p = 4B2 ln r + p0 . 2 (13) Here, we correct a typographical error in B1 on p. 126 of [14] and on p. 153 of [5]. 3 where p0 is a constant which sets the pressure level. It is immediately obvious from (12) that the coefficients are singular when 2α cos 2α − sin 2α = 0, which occurs at a critical value αc determined by tan(2αc ) = 2αc . In the range of interest, i.e. for 0 < αc < 180◦ , we have αc ≈ 128.7◦ . It is noteworthy, that in the limit r → 0, the velocity scales linearly with r whilst the pressure is logarithmically singular at the corner and is independent of the angular coordinate θ. 4 Numerical results From a numerical viewpoint, the steady fixed-boundary problem considered in this paper is complicated only by the presence of a singularity in the pressure which, according to (13), is logarithmic as r → 0. The simplicity of the rest of the problem and the availability of asymptotic results, which not only give the behaviour of the velocity and pressure near the corner, but also provide the coefficients, make this a perfect testing ground for a numerical method’s ability to approximate flows in corner regions formed by boundaries on which different types of boundary conditions are applied. In the standard implementation of the finite-element, as well as finite-difference, algorithm, one assigns an a priori unknown finite value to the pressure at the corner point. If such a code attempts to approximate a solution where the pressure at the corner point is singular, like the one whose local asymptotics we considered earlier, the nodal value of the pressure, as well as the pressure at the neighbouring nodes, will vary as one refines the mesh. In other words, an attempt to approximate a singular analytic solution using regular numerical representations of the unknown function on each element will lead to a numerical ‘solution’ that is mesh-dependent and hence, strictly speaking, it is not a solution to the original problem formulated in terms of PDEs that ‘do not know’ about any mesh. In order to achieve a uniformly valid solution in the framework of a finite-element method, one approach is to use singular elements, i.e. to redefine the pressure interpolation in the elements that contain the corner point node in such a way that the pressure is allowed to be infinite at the corner point and behave as described by the local asymptotics. A simple implementation of this idea is described in [4]. In the present work, the problem formulated in Section 2 has been considered using part of a finite-element-based numerical platform which has been developed to simulate a range of microfluidic capillary flows and has already been used to obtain new results for the flow of liquids over surfaces of varying wettability [15, 16]. The idea here, is to modify the finiteelement’s basis function Φp associated with the pressure at the corner point. In the standard triangular Taylor-Hood element with six velocity nodes and three pressure nodes, Φp is linear and takes the value 1 at the corner point and 0 at all other pressure nodes. Now, instead of using Φp and determining the coefficient in front of it, we will be using and determining the coefficient in front of a singular basis function Φs = Φp ln r. We denote the computed coefficient of Φs as Bp . Instead of Φp one could use other functions to pre-multiply the logarithm, for example sin(πΦp /2); such functions have been tried and it was found that they do not offer noticeable advantages over Φp . It is pointed out in [4] that the usual Gaussian quadrature is not well suited to the integration of a singular function, such as ln r, and a special Simpson quadrature routine has been suggested to provide a very accurate approximation of the integrals. This routine has also been incorporated into our numerical platform; however, we found that using Gaussian quadrature with enough integration points provided an accurate enough estimation of the integral and was significantly quicker and easier to implement. Sixteen Gauss integration points were found to be more than sufficient. 4 1.0 θ=75° 0.1 u 1 0.5 0 2 θ=0° 0.0 0.5 1.0 −0.1 0 2 r × 103 4 Figure 1: Left: streamlines in increments of ψ = 0.05, with the θ = 0, α interfaces corresponding to ψ = 0. Right: comparison of the computed velocity with the analytic prediction (dashed line) along the liquid-solid (1) and liquid-fluid (2) interfaces. At the critical angle αc ≈ 127.8◦ , more complex asymptotic analysis should be considered to resolve the singular behaviour and the results should be incorporated into the code in a way similar to what is being described. However, since this paper is concerned with the general numerical approximation of corner flows which contain singularities and their numerical treatment, we shall consider angles away from this critical value, i.e. the range where, as one might expect at this stage, our asymptotics of Section 3 can be used. First, we shall consider subcritical angles α < αc , in particular α = 75◦ as a representative case. We take s = 0.1 and R = 10 in all the simulations that we present as an investigation into the variation of these parameters would be about the physical problem rather than its numerics and would detract from the main emphasis of this paper. 4.1 Numerical approximation at subcritical corner angles In Fig. 1, we show the streamlines generated by the exponential slip model. As expected, the prescribed velocity on the solid draws fluid near the solid out of the corner, thus reducing the pressure there which then sucks in fluid from the far field. In the same figure, the components of the radial velocity along the interfaces are compared to the asymptotic predictions; the agreement between the calculated velocity along the liquid-fluid interface and the asymptotic prediction is visibly excellent (agreement along the liquid-solid interface near the corner is guaranteed as the velocity is prescribed, so we give it here just to show the range in which the velocity varies). The plot of pressure along the two interfaces in Fig. 2 shows that the pressure is indeed θ-independent, and it is almost graphically indistinguishable from the asymptotic result. To confirm the mesh-independence of the result, we also consider how well the singular behaviour of pressure is captured as the mesh is resolved over ten orders of magnitude using ten different meshes, characterized by the width of the smallest element r1 . To do so, we show both the coefficient Bp in front of the basis function Φs (curve 1 in Fig. 2) and the appropriate gradient determined from the pressures p1 , p2 at the two pressure nodes on the solid surface, for simplicity, closest to the corner point at radial distances r1 , r2 , which is given by (p2 −p1 )/(ln r2 −ln r1 ) (curve 1g). This second method allows us to compare the code with the singular elements to one without (curve 2g in the plot) where Bp is not explicitly calculated. The clearest conclusions are drawn from the results of the second method (curves 1g and 5 −50 p 9 Bp 8.5 8 7.5 1 1g −75 −100 r 10 −6 7 6.5 −4 −5 2g −10 −5 r1 10 0 10 10 10 10 Figure 2: Left: comparison of computed pressure along the θ = 0, α interfaces, which are graphically indistinguishable, as the corner is approached with the analytic pressure (dashed line). Right: convergence of the pressure coefficient of ln r with (1) and without (2) singular elements as the mesh is resolved. Bp is calculated from the singular element’s coefficient (1) and from the local gradient (1g) and (2g). The analytic value is 7.23 and r1 is the size of the smallest element. 0.1 0 −0.1 −0.2 θ=0° 0.5 0.0 0.5 1.0 u θ=175 1.0 ° −0.3 0 2 r × 103 4 Figure 3: Left: streamlines in increments of ψ = 0.1 with the θ = 0, α interfaces corresponding to ψ = 0. Right: comparison of the velocity along the liquid-fluid interface with the analytic prediction (dashed line). 2g). Here, the plot shows that, by using the singular element, the local gradient converges to the asymptotic value of Bp = 7.23 (the dashed line in the figure): without these singular elements the code converges to the incorrect value as the mesh is resolved. This is strong support for the inclusion of singular elements in dynamic wetting codes. Without them, the behaviour of pressure is wrong not only in the element adjacent to the corner point, but, by continuity, also in a neighbourhood of this element. For α < αc the asymptotics and numerics are in excellent agreement, and the special treatment of the corner singularity was a success. Now we consider a supercritical angle α = 175◦ > αc . 4.2 Numerical approximation at supercritical corner angles The streamlines in Fig. 3 are as one may intuitively expect, but when we compare the numerical and asymptotic results for the velocity along the liquid-fluid interface there is no agreement. In fact, the asymptotic result predicts that the flow should be up the liquid-fluid interface, which is clearly not the case in the computed solution. 6 5000 4000 3000 2000 1000 0 −1000 p 0 −2 2 B ×10− 4 p −4 1 r 10 −6 −6 −4 r1 10 −10 10 −5 10 10 −5 10 0 Figure 4: Left: computed pressure in the vicinity of the corner along the liquid-solid (1) and liquid-fluid (2) interfaces compared to the analytic prediction (dashed line). Right: pressure coefficient compared to the analytic value 1.122 (dashed line) plotted against r1 the size of the smallest element. The computed pressure along the liquid-solid interface is given as curve 1 in Fig. 4. It is not only that the pressure strongly deviates from the analytic result (dashed line) as the corner point is approached; one can see that there also appear huge oscillations: the pressure decrease in the element adjacent to the special corner elements (the line to the left of the point 10−6 in the plot) is followed by a steep increase in the element comprising the corner point (not shown in the semilogarithmic plot). Such mesh-dependence of the numerical result indicates that the obtained solution cannot be regarded as a valid approximation of the solution to the original set of partial differential equations. This conclusion is re-enforced when we study the value of Bp , the coefficient to the logarithm in the singular element, as we refine the mesh: there is no convergence. A similar trend is observed if we study Bp using the local gradient method. The same conclusions may be drawn for all supercritical angles. Thus, the standard FEM coupled with the local-asymptotics-based approximation of the pressure does not allow one to obtain an acceptable numerical approximation for solutions of the Stokes equations in a corner with a combination of Dirichlet and Neumann boundary conditions on the interfaces for all angles greater than 127.8◦ . The robustness of the obtained numerical ‘solution’ suggests that there is a fundamental numerical problem. It should be emphasized that we arrived at this difficulty just by varying the wedge angle in a code that produces excellent results for smaller wedge angles, and then cannot provide any mesh-independent solution after a critical angle. An immediate (tentative) explanation for this situation is that, besides the analytical solution whose asymptotics has been considered in Section 3 and incorporated into the code, there exists a ‘local eigensolution’, i.e. a solution satisfying zero boundary conditions on the sides of the wedge, that becomes dominant for the supercritical angles. We will now examine this conjecture. 4.3 Asymptotics of an eigensolution ∂ψ = 0, ∂θ ∂ 2ψ = 0, ∂θ2 The near-field asymptotics of the eigensolution to the biharmonic equation satisfying conditions ψ= is given by ψe = rλ [A1 sin (λθ) + A2 cos (λθ) + A3 sin ((λ − 2) θ) + A4 cos ((λ − 2) θ)] . 7 (15) on θ = 0, and ψ= on θ = α, (14) Using (14) and noting that λ = 1, 2, we find that λ is determined by the equation: 2 sin (λα) cos ((λ − 2) α) = λ sin (2α) . Defining the degree of freedom by A ≡ A1 , the boundary conditions (14) give: A2 = −A4 = −A λ(λ − 2)−1 sin ((λ − 2) α) − sin (λα) , cos ((λ − 2) α) − cos (λα) A3 = −A λ . λ−2 (17) (16) and the pressure has the form: pe = 4A(λ − 1)rλ−2 [a3 cos((λ − 2)θ) − a4 sin((λ − 2)θ)] + p1 (18) where ai = Ai /A, i = 3, 4 and p1 is a constant setting the pressure level. Promisingly, equation (16) has roots λ ∈ (1, 2) for α ∈ (αc , 180◦ ), with λ → 2 as α → αc , λ → 3/2 as α → 180◦ and λ as a function of α varying monotonically between these limiting values. This eigensolution has been derived in a number of other works considering the flow in a corner formed between an impermeable no-slip boundary and an impermeable, sometimes free, zero-tangential stress boundary, e.g. in [17]; these are sometimes referred to as ‘stick-slip phenomena’ [18, 6]. The difference between our problem and the aforementioned flows with static corner points is that, unlike these situations where the flow is driven by the far field, here the fluid motion is generated by the movement of the solid and analytically this behaviour is captured in the asymptotics of Section 3. The eigensolution comes on top of this solution and, in the near field, it ‘does not know’ about the motion of the solid, although, ultimately, it is the solid’s motion that generates the flow in the far field that gives rise to this solution. The eigensolution exists in the range of subcritical angles as well, but there it is regular in all variables and therefore causes no problem for numerical computations; it is only for α > αc that the eigensolution becomes both singular and dominant. For α < αc the pressure at the corner point can be referred to as single-valued: the coefficient in front of the logarithm is independent of θ. In contrast, the solution for α > αc is manifestly multivalued as predicted by the eigensolution (18) and as seen numerically: if one takes a vicinity of the corner point, then, no matter how small this vicinity is, there will be points which are equidistant from the corner point with an arbitrarily large pressure difference. Here, being interested in the numerical side of the problem, we set aside physical arguments that might arise in connection with the obtained solution. The existence of an eigensolution and its dominance for α > αc suggests that our initial attempt at computing the flow at large angles were flawed because, given the asymptotics of Section 3, we assumed that the pressure scaled as ln r whereas in fact the most singular term (i) has order r−k where k = 2 − λ ∈ (0, 1) and (ii) is dependent on θ. This suggests a generalisation of our approach: we need to incorporate the new singular behaviour into the special elements adjacent to the corner point. Due to the presence of the eigensolution, we now have an unknown constant A in our asymptotics which will prevent us from comparing a priori determined analytic curves with our numerical results. However, once A is determined numerically, we may use it to extrapolate the asymptotic behaviour outside the singular element in which it is calculated, i.e. we may then compare, a now semi-analytic, asymptotic prediction to the computed solution globally. An alternative method that we used to verify our singular element solution and do not describe in detail here is to analytically remove the eigensolution, which is the cause of numerical difficulties, prior to computation and then superimpose it back on after. This approach that has been shown to be essential for the simulation of flows using the Navier slip condition, a Robin-type boundary condition, on the solid surface is described elsewhere [19]. In more complex problems, where the corner is just one element, this method of removing the eigensolution everywhere is overly complex and a local method should be used which removes the 8 5000 4000 3000 2000 1000 0 −1000 p 0.1 0.08 2 1 Ap 0.06 0.04 0.02 r 0 10 −10 r1 10 −5 10−6 10−5 10−4 10 0 Figure 5: Left: computed pressure distributions in the vicinity of the corner along the liquidsolid (1) and liquid-fluid (2) interfaces compared with the semi-analytic result (dashed lines). Right: convergence of pressure coefficient Ap to 0.092, plotted against r1 the size of the smallest element. eigensolution near to the corner; this method is also described in [19] using examples of full scale dynamic wetting simulations. 4.3.1 Modified singular elements In the limit as r → 0, the first two terms in the expansion of pressure are p = Ap rλ−2 g(θ) + Bp ln r. For α > αc both terms are singular. We will begin by using only the leading order term in r in our singular elements and will return to the two-term approximation later. Taking the leading term, our new singular elements have a basis function of the form: Φs = Φp ln r, α < αc , Φp rλ−2 g(θ), α > αc , (19) where the unknown coefficients are Bp and Ap , respectively. In Fig. 5 with α = 175◦ and hence, from (16), λ = 1.529, we see that the implementation of the new singular elements solves our previous problems by (i) removing the oscillations in pressure as the corner point is approached, and (ii) converging as the mesh is refined. The value of A = Ap , which is the coefficient of the singular basis functions and is determined by the finite element method, may be used after computation with the supplementary asymptotics of Section 4.3 to produce a fully determined asymptotic solution for the velocity and pressure. Then we may compare the analytic results of Section 4.3 with our numerical results globally. It should be pointed out that using the value of A to extrapolate the analytic behaviour of the eigensolution well outside the first elements provides a quick check to see if A is in the correct range, this value does not in any way actually determine the velocity field outside the first elements. The comparison with pressure in Fig. 5 shows good agreement between numerics and asymptotics; however, very close to the corner point along the free surface the numerical solution is not as smooth as the asymptotic result. This is no surprise given the huge gradients in pressure which are being approximated by linear basis functions both in the radial and angular directions. In Fig. 6, we compare the velocity along the interfaces of the computed numerical solution to the asymptotic result, showing in particular how the full asymptotic solution is a superposition of the eigensolution ue and the supplementary solution ua . We see that, although the supplementary asymptotics predicts that flow will be reversed near the contact line, this is blown away by the strength of the eigensolution, which restores what one would intuitively think is the correct direction for the flow. The agreement we see in this figure is sufficient when 9 0.1 ua 0 −0.1 −0.2 −0.3 0 2 u ue u +u a e r ×10 4 3 Figure 6: Comparison of the computed velocity along the liquid-fluid interface u with the semi-analytic result which is shown decomposed into a contribution from the supplementary asymptotics ua and the eigensolution contribution ue with A=0.092. we consider that it is determined by the coefficient of the singular pressure which only plays a role in elements adjacent to the corner point. 4.3.2 Numerics incorporating two-term asymptotics For an angle of α = 175◦ we have completely resolved the situation. However, we have observed for smaller angles, roughly αc < α < 150◦ , that, in terms of mesh refinement, convergence is very slow, i.e. a mesh-independent regime is only realised for exceptionally well resolved meshes which are well outside the scope of most numerical platforms where the corner is just one part of a larger problem. In this range of angles, although asymptotically in the limit r → 0 the logarithmic pressure behaviour (p ∼ ln r) is overshadowed by the eigensolution, where p ∼ rλ−2 g(θ), in reality, unavoidable finiteness of the resolution of the mesh means that one could be sufficiently far away from the corner point for the logarithm to be still dominant in the numerics. To see if this is indeed the case, we consider the relative size of the two pressure terms at the edge opposite the corner point in the first elements. Taking this element to have size 10−n , we see that the two singular functions are comparable, at a critical value λc , at the edge of the first element when ln 10−n = 10−n(λc −2) , (20) ln | − n ln 10| . (21) n ln 10 Taking n = 6, from (21) we have that λc = 1.81, which means, using (16), that the logarithmic behaviour dominates a numerical scheme, with the stated spatial resolution, in the range of angles αc < α < 143◦ . Thus, we have seen both numerically and analytically, that even for an extremely well resolved mesh, we should expect that there is a range of angles in which the logarithmic solution cannot be neglected; approximately for α < 150◦ . This is a very serious issue which must be resolved as almost all numerical schemes will not be able to afford the required resolution to accurately capture the singular behaviour. The solution to this issue is to include the second term in the asymptotic expansion of pressure, so that p = Ap Φp rλ−2 g(θ) + Bp Φp ln r for α > αc . The coefficient to this logarithm Bp is known exactly from the supplementary asymptotics (13) and we have already shown in Section 4.1 that this value is numerically reproduced. Therefore, the simplest solution is to prescribe the value of Bp and see if this improves the speed of convergence of Ap . In Fig. 7 we show how nicely this method works with, as one would expect, the most improvement occurring for smaller angles where the logarithmic pressure is strongest. For example, for λc = 2 − 10 that is when 3 Ap 140 ° 2 1 150 ° 160° 170° −10 −5 r1 10 10 0 0 10 Figure 7: A comparison of the convergence of the pressure coefficient Ap with mesh size r1 . With (dashed line) and without (full line) the supplementary asymptotics prescribed, for a range of different angles. α = 150◦ , where the converged value is Ap = 1.007, if we neglect the logarithm, then we require 5000 elements to get within 5% of this value, but, if we include the asymptotic logarithmic behaviour then we only need 2500 elements to attain the same accuracy. Computationally this is a terrific saving and, given the simplicity of its implementation, we conclude that the logarithmic asymptotic behaviour should be included for all angles α > αc . 5 Conclusion We have shown that accurate numerical calculation of viscous flows near corners of the flow domain in the framework of the finite-element method requires special treatment of the elements adjacent to the corners. In the case of zero-stress/prescribed velocity boundary conditions on the sides of the corner, the range of corner angles is split into two distinct regions. For the angles below the critical angle αc ≈ 128.7, it is sufficient to use a logarithmic basis function for the pressure, whereas for supercritical angles there appears a ‘hidden’ eigensolution which considerably complicates numerics. In order to obtain an acceptable (i.e. mesh-independent) solution, one has to incorporate this eigensolution into the code by altering the basis functions for the pressure in the elements adjacent to the corner point. Close to the critical angle it becomes necessary to use a two-term asymptotics in the numerical algorithm, including the leading terms of both the eigensolution and the solution of the inhomogeneous problem, as the finiteness of the mesh size could result in the algebraically and logarithmically singular terms having comparable values. An issue that is now opened up for numerical and analytic investigation is how to generalize the developed methods for a three-dimensional case, i.e. in a situation where both the contact angle and the direction of velocity of the solid vary along a contact line. The first of these aspects becomes particularly challenging when the angle varies from subcritical to supercritical. Acknowledgements The authors kindly acknowledge the financial support of Kodak European Research and the EPSRC via a Mathematics CASE award. 11 References [1] J. Juntunen and T. Tsiboukis. On the fem treatment of wedge singularities in waveguide problems. IEEE Transactions on Microwave Theory and Techniques, 48:1030–1037, 2000. [2] G.C. Georgiou, W.W. Schultz, and L.G. Olson. Singular finite elements for the suddenexpansion and die-swell problems. International Journal for Numerical Methods in Fluids, 10:357–372, 1990. [3] M. Duflot. A meshless method with enriched weight functions for three-dimensional crack propagation. International Journal for Numerical Methods in Engineering, 65:1970–2006, 2006. [4] M.C.T. Wilson, J.L. Summers, Y.D. Shikhmurzaev, A. Clarke, and T.D. Blake. Nonlocal hydrodynamic influence on the dynamic contact angle: Slip models versus experiment. Physical Review E, 83:041606, 2006. [5] Y.D. Shikhmurzaev. Capillary Flows with Forming Interfaces. Taylor & Francis, London, 2007. [6] G.C. Georgiou, L.G. Olson, W.W. Schultz, and S. Sagan. A singular finite element for Stokes flow: the stick-slip problem. International Journal for Numerical Methods in Fluids, 9:1353–1367, 1989. [7] C. Huh and L.E. Scriven. Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. Journal of Colloid and Interface Science, 35:85–101, 1971. [8] C.L.M.H. Navier. M´moire sur les lois mouvement des fluides. M´m. de l’Acad. de Sciences e e l’Inst. de France, 6:389–440, 1823. [9] E.B. Dussan V. The moving contact line: The slip boundary condition. Journal of Fluid of Mechanics, 77:665–684, 1976. [10] M. Zhou and P. Sheng. Dynamics of immiscible-fluid displacement in a capillary tube. Physical Review Letters, 64:882–885, 1990. [11] S. Somalinga and A. Bose. Numerical investigation of boundary conditions for moving contact line problems. Physics of Fluids, 12:499–510, 2000. [12] E.B. Dussan V and S.H. Davis. On the motion of a fluid-fluid interface along a solid surface. Journal of Fluid Mechanics, 65:71–95, 1974. [13] H.K. Moffatt. Viscous and resistive eddies near a sharp corner. Journal of Fluid Mechanics, 18:1–18, 1964. [14] Y.D. Shikhmurzaev. Singularities at the moving contact line. Mathematical, physical and computational aspects. Physica D, 217:121–133, 2006. [15] J.E. Sprittles and Y.D. Shikhmurzaev. Viscous flow over a chemically patterned surface. Physical Review E, 76:021602, 2007. [16] J.E. Sprittles and Y.D. Shikhmurzaev. A continuum model for the flow of thin liquid films over intermittently chemically patterned surfaces. The European Physical Journal Special Topics, 166:159–163, 2009. 12 [17] D.M. Anderson and S.H. Davis. Two-fluid viscous flow in a corner. Journal of Fluid Mechanics, 257:1–31, 1993. [18] S. Richardson. A stick-slip problem related to the motion of a free jet at low Reynolds numbers. Journal of Fluid Mechanics, 67:477–489, 1970. [19] J.E. Sprittles and Y.D. Shikhmurzaev. Viscous flow in domains with corners: Numerical artifacts, their origin and removal. Submitted to the Journal of Fluid Mechanics, 2009. 13