Lubrication dynamics of a settling plate

Abstract If a flat, horizontal, plate settles onto a flat surface, it is known that the gap $h$ decreases with time $t$ as a power law: $h\sim t^{-1/2}$. We consider what happens if the plate is not initially horizontal, and/or the centre of mass is not symmetrically positioned: does one edge contact the surface in finite time, or does the plate approach the horizontal without making contact? The dynamics of this system is analysed and shown to be remarkably complex. We find that, depending upon the initial position of the plate and the position of the centre of force, the plate might either make contact in finite time or settle progressively without ever making contact. Our results show an excellent agreement between analytical exact solutions, asymptotic solutions and numerical studies of the lubrication equations.


Introduction
We discuss a foundational problem in lubrication theory.We consider the motion of a body with a flat lower surface settling onto a flat horizontal surface, impeded by a thin layer of viscous fluid.For simplicity, we consider the two-dimensional case, where the gap h at time t depends only upon one Cartesian coordinate of the plane (x, say), and is independent of y.It is also assumed that the plate is sufficiently wide that fluid motion in the y-direction can be neglected.If the plate is initially horizontal, and the weight acts through the centre, then it can be shown that the plate remains horizontal as it settles, with h ∼ t −1/2 , in accordance with the classic lubrication result that smooth solid objects take an infinite amount of time to make contact ([Brenner (1961)]).
If the force is applied off-centre, and/or the plate is not initially horizontal, then the system is much harder to analyse.There are three coupled degrees of freedom: the gap Z between the left-hand edge of the plate and the surface, the angle of the plate θ, and the horizontal displacement X of the left hand edge (both θ and Z must be assumed to be small, in order for the lubrication theory approximations to be valid).The geometry of the system is illustrated in figure 1.
The equations of motion are the Reynolds equations of lubrication theory ([Reynolds (1886)]), reviewed in [Michell (1950), Szeri (1998)].The plate may, or may not, contact the surface in a finite time.We find that there are two different types of solution, those which the edge of the plate makes contact in finite time, and those where the gap decreases as a power-law.
The motion of the plate is determined by three generalised forces: the vertical force F z , the horizontal force F x , and the torque (clockwise, about the left-hand edge), G. Section 2 derives the equations of motion, by determining the resistance matrix which relates the force vector (F x , F z , G) to the velocity vector, ( Ẋ, Ż, θ).A similar expression for the resistance matrix was obtained in [Cawthorn and Balmforth (2010)] for the problem of a falling wedge, but our Figure 1: The plate has width L and its centre of mass is displaced from the left-hand edge by sL.A coordinate x measures distance from the left-hand edge, and all variables are assumed to be independent of the other Cartesian coordinate in the plane, y.At a time t, the angle of the plate is θ(t), the gap at the left-hand edge is Z(t) and the horizontal displacement of the left-hand edge is X(t).
equations do contain an additional term, which makes a significant difference to their solution.Section 3 considers the reduction of the equations of motion to dimensionless form.
Sections 4, 5 and 6 discuss the solutions to these equations of motion.In section 4 we show that the Reynolds lubrication equations are exactly solvable, in terms of integrals of elementary functions, but that the general form of the solution is too complicated to be informative.This leads us to concentrate on understanding the qualitative form of the solutions, and their asymptotic approximations.We find that two qualitatively different types of solution exist.Section 5 considers solutions for which θ/Z approaches a constant, using the dimensionless form of the equations of motion.Section 6 discusses solutions where the plate contacts the surface on one edge in finite time.In this section we find it more convenient to use the original dimensional form of the equations.Section 7 details some numerical experiments, comparing our solutions with numerical integration of the Reynolds equations.Section 8 summarises our conclusions.

Equations of motion
A flat plate of width L, immersed in a viscous fluid of viscosity µ and constant density ρ, settles onto a flat horizontal surface.We consider the two-dimensional case, where the gap h between both surfaces depends only upon one Cartesian coordinate of the plane (x, say), and is independent of y.It is also assumed that the plate is sufficiently wide that fluid motion in the y-direction can be neglected.
The two surfaces move with velocities v 1 (t) (lower surface) and v 2 (t) (upper) and we assume that their separation h(x, t) is sufficiently small that lubrication theory applies (meaning that there is Poiseuille flow in the gap).The volume flux per unit depth is J(x, t), the flux and pressure p are related by and the continuity equation is Poiseuille flow also implies that the tangential stress on the upper surface is We can use these equations in any convenient Cartesian frame.Rather than using the laboratory frame, let us use a frame which is attached to the plate.The displacement of the left-hand edge of the plate relative to some fixed position in the laboratory frame is X(t).The velocity of the upper bounding surface in the plate frame is then v 2 = 0, and the velocity of the lower surface is v 1 = − Ẋ.The configuration of the plate is then described by specifying height at displacement x from the left edge: where Z(t) is the gap between the left hand edge of the plate and the surface, and θ(t) is the angle between the plate and the horizontal, see Fig. 1.We wish to determine the force components in the vertical and horizontal direction, F z and F x , respectively, and torque G on the plate.These are: where p 0 = p(0).The above expressions depend on the pressure gradient ∂p/∂x, which we can determine from Eq. (1) as follows.Inserting Eq. ( 4) into the continuity equation ( 2) and integrating once we obtain where J 0 is a constant.It will be useful to define We can now find the pressure gradient where J 0 is found by integrating the above expression and imposing that p(L) = p(0): Provided the background pressure is high enough to prevent cavitation, its value must be irrelevant, so we set p 0 = 0. Therefore, the force components are: We can write the relation between forces and velocities in terms of a matrix: where the matrix A is Equations of motion in which a generalised velocity vector is obtained from a generalised force vector by matrix multiplication occur quite generally in treatments of viscosity dominated flow, and the matrix A is termed the resistance matrix ( [Happel and Brenner (1983)]).The equations of motion ( 11) and ( 12) were obtained in a slightly different form in [Cawthorn and Balmforth (2010)].
In that work the contribution to the horizontal force which arises from tilting the plate was not included, and the first element of the force vector is F x , rather than F x + θF z .

Dimensionless equations
It will be efficient to define non-dimensional dynamical variables: and a non-dimensional version of the integrals I n m : with With these definitions the generalised force components are expressed in terms of the dimensionless dynamical variables as follows: It will also be helpful to introduce another dimensionless variable λ which can be used in place of ξ, defined by: λ = ln ξ .
Then the force equations are where the coefficients B ij (η) are obtained by comparison with ( 16), e.g.
The equations of motion can then be expressed in the form We are interested in settling of the body, due to a gravity force W which acts through a point x = sL from the left-hand edge, so the opposing forces on the body due to the fluid are (F x , F z , G) = (0, W, sW L).Introduce a transformed time variable τ , which satisfies For our settling plate we then have the following equations of motion: The matrix elements of B(η) and C(η) can be determined analytically: the expressions are given in Appendix A.
It is possible to find solutions in which Z(t) is increasing, but the height of the centre of gravity must always decrease monotonically.It is not obvious that Equations (23a)-( 23a) have this property.Appendix B demonstrates that this property holds.

Solution of dimensionless equations
In subsection 4.2 it will be shown that the equations of motion, ( 23), can be solved, expressing t, ξ and ζ as functions of η, defined as integrals of elementary functions.These exact solutions are, however, quite difficult to interpret.It is the qualitative behaviour of the solutions which is usually of more interest than precise expressions.In particular, it is desirable to understand whether the edge of the plate makes contact in finite time, or whether the angle of the plate decreases so that it settles without ever making contact.This is addressed in the next subsection 4.1, before we consider the exact solution.

Fixed points
The qualitative form of the solution is addressed by looking at the dynamics of the dimensionless variable η.If the left-hand edge makes contact, this corresponds to η → ∞ (and a contact of the right-hand edge is η → −1).
Consider equation (23c).This is a differential equation for η(τ ) which is independent of the other variables.Any fixed points of the η(τ ) dynamics are determined by the condition This determines a fixed point η * (s) as a function of the position of the centre of gravity, s.The fixed point is stable if (and, conversely, unstable if κ < 0).The matrix elements quoted in equations (72) of Appendix ?? enable us to express F η (η, s) in terms of elementary functions: we find that the values of s and η at a fixed point are related by i.e., η = −1 (right-hand edge approaches contact) is a fixed point for all s.By symmetry, the limit η → ∞ (i.e., left-hand edge approaches contact) can also be regarded as a fixed point.However, the function F η (η, s) has a singular form as η → −1, and the dynamics in the vicinity of these edge-contact fixed points is highly unusual.In section 6 we consider in some detail the case where the edge is close to contact.There we argue that the plate makes contact in finite time.
Panels 2(b, c, d) show the behaviour of F η for three different ranges of η as we increase its value above η = −1.We can see how fixed points are found at increasingly large values of η.It would be instructive to plot a phase-diagram which shows the locations of the stable and unstable fixed points as lines in the (η, s) plane, but the fact that η may approach infinity is inconvenient.There is a left-right symmetry of the system, such that a fixed point at (η, s) should be reflected as a fixed point with the same stability at (−η/(1 + η), 1 − s).If we were to Figure 3: Plot showing fixed points in phase space of model, using φ, defined by ( 28) and ( 13), to describe the configuration of the plate, and s to parametrise the position of its centre of gravity.The blue line indicates stable fixed points, the red line unstable.use an alternative dimensionless parameter, to describe the configuration of the plate, then the configuration space extends from φ = −1 (right-hand edge in contact) to φ = +1 (left-hand edge in contact), and a fixed-point at (φ, s) is mirrored by one at (−φ, 1 − s).In terms of the symmetrized configuration variable φ, the equation for the fixed points is The function S(φ) satisfies the symmetry relation Figure 3 shows the phase-diagram of the system using the (φ, s) coordinates.If we include the stable fixed points at φ = ±1, there are either three or five fixed points, depending upon the value of s (the center of mass position): 1.For 0 < s < 0.3764208755 . . .there is an unstable fixed point at a value of φ * which is very close to the point at which the right-hand side touches (note the dashed red line and solid horizontal blue lines are indistinguishable in Fig. 3).In this range of s, φ * lies in (−1, −0.999206 . ..].If the right-hand edge is exquisitely close to making contact, it does so.Otherwise, we expect that the left hand edge will make contact.In section 6 we argue that, in both cases, contact happens in finite time. 2. For 1 > s > 0.623579125 . .., the situation is the mirror image of the case 0 < s < 0.3764208755 . ... There is an unstable fixed point at a value of φ * which lies in the interval [0.999206 . . ., 1).If φ is initially greater than φ * , the left-hand edge contacts in finite time.If φ is less than φ * (which is the only case which is of practical relevance), the right-hand edge contacts in finite time.
3. There is a region in the interval 0.3764208755 . . .≤ s ≤ 0.6235791253 . . .which has two stable fixed points at φ = ±1, two unstable fixed points, and one other stable fixed point between the unstable fixed points.There is a basin of attraction for which the value of φ approaches the non-trivial stable fixed point φ * (s).This stable fixed point starts at φ * (0.3764208755 . ..) = 0.83774254 . .., passes through φ * (1/2) = 0, and disappears at φ * (0.6235791253 . ..) = −0.83774254 . ... The basin of attraction in the φ variable is at its largest when s = 1/2.At this point, the unstable fixed points are at φ = ±0.9956216,and all initial conditions in between these values converge to the plate settling in the flat, φ = 0, configuration.
In section 5 we show that the approach to the fixed point has a power-law dependence upon time t.

Exact solution
We define a dimensionless time variable where ξ 0 is the initial value of Z/L.Equations ( 21) and ( 23) then lead to the following differential equations for t(η), λ(η) and ζ(η): These three differential equations can be integrated to express t, λ = ln ξ and ζ as functions of η, for example (32a) gives where λ(η) is found by integrating (32b): These expressions give an exact expression for t in terms of integrals of elementary functions of η, but they are too complicated to be instructive.Furthermore, in order to express Z, θ and X as functions of t, it is necessary to invert the function t(η) defined by ( 33) and (34).

Asymptotic solution near stable fixed-point
If the initial conditions are in the basin of attraction of a stable fixed point η * (s), then the long-time behaviour of the solution is determined by the solution in the vicinity of this fixed point.
If η is close to η * then so, with η 0 = η(0), the solution of ( 23c) is approximated by where κ(s) was defined in Equation ( 25).When η(τ ) converges to this stable fixed point at η * , the functions F ζ and F λ in (23a), (23b) approach constant values.Equation (23b) then implies that ξ(τ and hence equation ( 21) can be integrated to give t as a function of τ and ξ 0 : so that, (noting that (37) implies that we must have F λ < 0), the leading term in the dependence of τ upon time is proportional to ln t.Using (37) in ( 38) we obtain We can now determine dependences upon the true time, valid for t → ∞.In the limit as t → ∞, Equation (39) then implies and Equation ( 36) can be written as The long-time evolution of ζ(t) is (assuming X(0) = 0): (where K and ζ are independent of t), so that the plate slides by an unbounded distance (unless In the special case where the centre of mass is symmetrically placed, so that s = 0, it is expected that the fixed point of (23c) will be η * (1/2) = 0.So let us expand the matrix C about η = 0.The leading order of K n m (η) when |η| ≪ 1 is: so that Using the Maple computer algebra system, we obtain up to second order: so that Hence κ(1/2) = 5 and the exponent of η for s = 1/2 is The motion of X in the symmetric case is different from Equation ( 42) because F ζ (0, 1/2) = 0.In this case, using equation ( 36) with τ ≫ 1, and η 0 ≪ 1, equation (23a) becomes Integrating this expression, we find that when η ≪ 1 and s = 1/2, ζ(τ, η 0 ) ≈ η 0 10 [exp(−5τ ) − 1].In the limit as τ → ∞ and η 0 → 0, the horizontal displacement therefore remains finite: 6 Asymptotic solution near and following edge contact In this section we first consider (section 6.1) the form of the solution of the lubrication equations when η ≫ 1, so that the left-hand edge is very close to contact, showing that contact occurs after a finite time.We then consider (section 6.2) what happens after the left-hand edge has contacted the plate.Because this latter calculation has to start from first principles, all of the calculations in this section will use the full, dimensional, equations of motion.The results of section 6.1 can be reproduced using the η → ∞ limit of the coefficients in equations ( 23).

Solution in the vicinity of left-edge contact
The equation of motion for η developed in section 3 always has a stable fixed point at η = −1 (right-hand edge contacts), and by symmetry there must be a corresponding motion in which η → ∞.Usually, a stable fixed point implies an exponential approach, but without the fixed point being reached for any finite time.However, the form of the equations on the vicinity of η = −1 and η → ∞ is so unusual that it is possible that η ≡ Lθ/Z → ∞ in a finite time.In the following, we consider that possibility.When Z is very small, the integrals I n m have simple asymptotic approximations: that is Then, equations ( 16) become We now set F x = 0, F z = W , G = W sL, and retain leading order terms in the variable Lθ/Z ≡ η: From Equations ( 53a) and ( 53c) we find Using the above equation with ( 53b) and (53a), and keeping only the leading order terms in η = Lθ/Z, we obtain Next substitute ( 54) and ( 55) into (53b) to obtain θ which has solution where θ 0 = θ(0).Equations ( 54), ( 55) and ( 56) are equations of motion for X, Z and θ, derived assuming that η = Lθ/Z is large.We must consider whether ( 55) and ( 56) imply that η → ∞ in finite time.Define a variable K as in (57), and a dimensionless pseudo-time τ (distinct from that defined by ( 21)) as follows: Noting that θ = ηZ/L, the equation of motion for θ is then Using (55), and (57), we obtain a differential equation for η(τ ): Figure 5: Schematic illustration of asperities on the surfaces.This simplified model assures that the pressure field is essentially unchanged across most of the surface, while the gap h cannot be less than the prominence of the asperities, ǫ.
We now write X = 1 + 4sθ 2 0 τ /3, and consider the limit as η → ∞, where this equation of motion is approximated by which has solution (with η = η 0 at τ = 0) When η 0 is very large, η → ∞ in a finite time t, which is approximated by

Motion after contact
We have seen that, if the initial value of η ≡ Lθ/Z is sufficiently large, the left hand edge contacts the surface in finite time.In order to complete our analysis of the problem, we need to consider what happens after contact is made.The motion after making contact depends upon the precise nature of the surfaces.The point of contact is able to tilt, and may also be able to slide.We assume that there are asperities on the contacting surfaces, such that the gap Z never falls below a microscopically small value, ǫ.We assume that the typical separation between the asperities, ∆L, satisfies L ≫ ∆L ≫ ǫ.We also assume that θ satisfies θ∆L ≪ ǫ, so that the size of the gap is nowhere reduced significantly due to tilt of the contacting surfaces.If these asperities are widely-spaced, as illustrated in figure 5, the equations of motion for the generalised forces which are mediated by the pressure in the gap, namely F z and G, may be assumed to be obtained by a simple modification of equations ( 53a) and (53b): after contact occurs, Z = ǫ, Ż = 0, and there is an upward reaction force per unit depth R at the left-hand edge.
Because this additional force acts through the axis which defines the torque G, it makes no contribution to the torque equation.The modified forms of equations (53a) and (53b) are, therefore: The equation for the horizontal force, F x , will be changed more fundamentally because this component is sensitive to the precise form of the asperities.However, we shall not need the modified version of equation (53c).
Equations ( 64) and ( 65) imply that so that the term proportional to Ẋ in ( 65) is negligible as ǫ → 0. After contact, the equation of motion for θ is very well approximated by θ and ( 66) implies that the sliding velocity, Ẋ, while being indeterminate because we have not determined the reaction force R, is very small.Combining ( 66) and ( 67), we find a relation between Ẋ and R: Note that the equation of motion for θ, (67) is the same as that which is applicable before contact (equation ( 56)), whereas (66) implies that the value of Ẋ changes upon contact.The values of Ẋ and R remain ambiguous because we have not specified how the equation for the horizontal force is modified when surfaces are in contact.We now assume that, whenever there is a non-zero reaction force between the plate and the surface, the plate is prevented from moving horizontally, so that Ẋ = 0 when R > 0. Then equation ( 69) implies that the reaction force is The condition that the reaction force R must be non-negative cannot, therefore, be satisfied, if s > 1/4 and Ẋ = 0, because the torque exerted by the downward force causes the left-hand edge of the plate to lift.If s > 1/4, the reaction force R is reduced to zero, and the plate is then able to slide.The sliding motion produces additional contributions to the force F z and the torque G, such that the left-hand edge remains in contact as it slides with reaction force R = 0.By symmetry, there must be sliding motion if 1/4 < s < 3/4.In this range of s, the slow sliding motion creates a downward force on the plate which counters the effect of the torque.After contact the plate can slide over the surface so that the reaction force is R = 0.The slip velocity when the left-hand edge is in contact is 7 Numerical studies We have compared our analytical predictions with numerical solutions of the equations of motion.In each case we compared the numerical solution with the exact solution, and with the appropriate asymptotic approximation.We used dimensionless variables in all of the simulations, equivalent to setting L = µ = ρ = W = 1.The exact solution given in subsection 4.2 was evaluated as follows.First we convert the initial values of the coordinates (X, Z, θ) to the dimensionless coordinates (ζ, λ, η).We then integrated equations (32a), (32b) and (32c) to obtain t, λ and ζ as functions of η.In the case where we expect convergence to a stable fixed point of η, these integrations can only be continued to a point close to the fixed-point value.The values of t were stored in a table along with the values of ζ, λ and η.The plots of the exact solution are generated by interpolating this table to plot η, ζ or ξ = exp(λ) as a function of t.
Figure 7 considers an asymmetric case, where s = 0.4, which has a stable fixed point at η * = 3.13266 (or φ * = 0.61033), and there are unstable fixed points at η = −0.99933and η = 33.20746 .We consider two different sets of initial conditions which are: (X, Z, θ) = (0, 0.001, 0.01) and (X, Z, θ) = (0, 0.005, 0.005), corresponding to initial values η = 10, (φ = 0.83333) and η = 1, (φ = 0.33333) respectively.Both initial conditions are within the basin of attraction of the stable fixed point.The double-logarithmic scale plot of ξ( t) shows the asymptotic power-law behaviour with exponent −1/2 also holds in the asymmetric case, but the exponent γ(s) for the decrease of η(t) takes a different value, plotted in figure 4(a).In figure 7, the horizontal displacement is plotted as a function of ln t, in accord with the prediction that X(t) increases logarithmically, with coefficient plotted in figure 4(b).
In both cases we plot η( t), ξ( t) and ζ( t) up to the point of contact.The contact time, t, is shown in each case.For the first set of initial conditions, η 0 is close to a fixed point, so that t is much larger than the prediction of Equation ( 63).Our second case has a relatively large initial η value of 100 and here the contact time is t = 2847, approaching the theoretical contact time of 2402 predicted by equation ( 63) for very large initial η values.

Concluding remarks
Our motivation for examining the lubrication theory for a settling plate came from the observation that there are two plausible types of solution: the plate might 'side-slip' onto the surface, making contact in finite time, or the motion might 'round out', so that the plate settles in a progressively flatter attitude, without ever making contact.We find that both of these types of solution may be realised, depending upon the initial position of the plate and the of the centre of force.We find that there is only one dimensionless coordinate which controls the nature of the trajectory: throughout most of the calculations we used η = Lθ/Z, but if we wish to represent the phase diagram of the system (Figure 3) in a manner which reflects its symmetry, the alternative choice φ = η/(2 + η) can be used.
The elements of both the resistance matrix (detailed in the Appendix) and its inverse can be expressed exactly as elementary functions of η.This observation was used to write an exact solution of the equations of motion, in parametric form, in terms of integrals of elementary functions.The exact solution (presented in section 4.2) is, however, difficult to interpret.This led us to investigating asymptotic approximations to the solution.We investigated two types of asymptotic solution: 1.There are solutions where η approaches a constant η * , with both Z and θ decreasing as t −1/2 .Numerical experiments show that these solutions only exist for s ∈ [s * , 1 − s * ], with s * ≈ 0.376 . ... The value of X(t) shows a logarithmic evolution, except for the case s = 1/2, and the relaxation of η(t) towards η * is a power-law, η − η * ∼ −γ(s) .The exponent γ(s) has a maximum γ(1/2) = 5/2, and approaches zero as s → s * or s → 1 − s * .

Figure 2 :
Figure 2: Plot of F η (η, s) as a function of η for different choices of the centre of mass position, s, and for four different ranges of η.Attractive fixed points of η occur where F η = 0 and F ′ η < 0.

Figure 2
Figure 2 shows plots of F η (η, s) as a function of η for different values of s.Because fixed points were found with widely differing values of η, the plots of Fig. 2 correspond to four different ranges of η.Panel 2(a) illustrates the limit: lim η→−1

Figure 4 :
Figure 4: (a) Plot of the exponent γ, which characterises the relaxation of η(t) towards its fixed point (equation (41)).(b) Coefficient of the logarithm, ζ, describing the slippage of the plate (equation (42)), as function of the centre of mass position parameter, s.

Figure 6 :
Figure6: Comparison of numerical solutions of the equations of motion with analytical predictions for the case where the centre of gravity is symmetrically placed, s = 1/2, and when the parameter η(t) = Lθ(t)/Z(t) approaches a stable fixed point.Panels (a,c,e) and (b,d,f) correspond to two different sets of initial conditions.Both η( t) (a,b) and the dimensionless plate height ξ( t) = Z( t)/L (c,d) show a power-law relaxation at large t, and the horizontal displacement ζ( t) = X( t)/L (e,f) approaches a constant.

Figure 7 :
Figure7: Comparison of numerical solutions of the equations of motion with analytical predictions for the case where the centre of gravity is displaced from the center, s = 0.4, and when the parameter η(t) = Lθ(t)/Z(t) approaches a stable fixed point.Panels (a,c,e) and (b,d,f) correspond to two different sets of initial conditions.The exponent for the relaxation of η, given by Equation (41), is different than in the symmetric case [compare panels (a,b) with Figure6(a,b)], but ξ( t) = Z( t)/L ∼ t−1/2 continues to hold in this asymmetric case (c,d).The horizontal displacement ζ( t) = X( t)/L (e,f) grows logarithmically, as given by Equation (42).