Leakage Threshold of a Saddle Point

The threshold condition for leakage inception is of great interest to many engineering applications, and it is essential for seal design. In the current study, the leakage threshold is studied by means of a numerical method for a mechanical contact problem between an elastic bi-sinusoidal surface and a rigid flat surface. The coalesce process of the contact patches is first investigated, and a generalized form of solution for the relation between the contact area ratio and the average applied pressure is acquired. The current study shows that the critical value of the average applied pressure and the corresponding contact area required to close the percolation path can be represented as a power law of a shape parameter, if the effect of the hydrostatic load from the pressurized fluid is ignored. With contact patches merged under a constant applied load, the contact breakup process is investigated with elevated sealed fluid pressure condition, and it is shown that the leakage threshold is a function of the excess pressure, which is defined as a ratio between the average applied pressure and the critical pressure under dry contact conditions.


S b
Critical solid-solid contact area when p f = p b S 0 Area of the solution domain p f Sealed fluid pressure Δ Geometry parameter of the initial gap g 0 (x, y) Side length of the square solution domain Material Poisson ratio Shape parameter of the initial gap g 0 (x, y) p Surface traction (surface stress orthogonal to the contact interface) p 0 Average applied pressure at vertical direction p saddle Surface traction at the saddle point p * 0 Critical pressure of dry contact p b Critical pressure of leakage p es Normalized excess pressure defined in Eq. (9) P t Total load at the vertical direction p s Pressure scaling factor E Material young's modulus g 0 Initial gap between the contacting surfaces g x The x-component of g 0 g y The y-component of g 0 l cx Half-width of the rectangular contact strip for g x l cy Half-width of the rectangular contact strip for g y p * x Pressure scaling factor for p x p *

Introduction
As a machine element that blocks the fluid flow, the seal achieves its function by compressing two surfaces together and closing the percolation paths between the high-and low-pressure sides. The contact patches are formed when two surfaces are pressed against each other. The initially isolated contact patches for rough surfaces can merge under increased contact load. The fluid leakage is stopped when the connected contact patches isolate the percolation paths between the high-and the low-pressure sides. Since the geometry connectivity of the contact patches decides the transition point of the leakage, it is of great interest to study the coalescence and breakup process of the contact patches.
The early contact mechanics models for the rough surfaces, such as the Greenwood-Williamson model [1] and the Greenwood-Tripp model [2,3], assumed that the contact between the micro asperities could be decomposed into a group of Hertzian-like contacts and the deformation of each asperity is independent of each other. Since then, both experiments and simulations [4,5] have shown that even for a simple contact problem between an elastic bi-sinusoidal and a rigid flat surface, the deformation of the contact patches is not independent of each other, suggesting that the contact patch interaction needs to be considered for the correct contact area, hence, also for the percolation threshold. Elastic contact mechanics accounting for the interaction between contact patches has been developed [6,7]. It has also been utilized to study the percolation threshold for both the isotropic fractal rough surface and the anisotropic rough surface [8][9][10]. The contact area ratio at the leakage threshold predicted by both contact mechanics theory and numerical simulations [11][12][13][14][15][16] is smaller than the value predicted by the bond-percolation theory [17] and the bearing area model [2,3], because the long-range elastic coupling between contact patches is ignored in the latter two.
Previous research on the leakage threshold of random rough surfaces usually builds on the assumption that the fluid pressure is much smaller than the contact pressure, and that the solid deformation generated by the fluid pressure is negligible. This assumption allows one to decouple the solid contact problem from the fluid leakage problem, and the leakage threshold is determined by the geometry connectivity of the deformed elastic surfaces without the fluid load. There are two situations where the fluid pressure is not negligible. First, the fluid can be trapped between the surface roughness due to the coalescing of the contact patches caused by a increase load. If the trapped fluid is further squeezed by the solid surfaces, its internal pressure will increase due to the volume decrease, thus providing additional load carrying capacity, and preventing further contact between the solid surfaces. Kuznetsov [18,19] studied the contact problem with fluid entraps in the valleys of a sinusoidal profile. His solution shows that the increased total load can be balanced by the compressed fluid pressure with a constant solid-solid contact length. The experimental findings in [20] have also confirmed that the fluid entrapped between a bi-sinusoidal and a flat surface prevents the further close-up of the gap between them after the adjacent contact patches touch each other. The further increase of the average applied pressure mainly contributes to the deformation near the saddle point, with little variation on the solid-solid contact areas, and the trapped fluid does not leak with its increased hydrostatic pressure [21].
Another situation where the fluid pressure must be considered is when the sealed fluid pressure is comparable with the contact load. It is not a rare case, especially near the boundaries of the contact zones adjacent to the highpressure side [22]. If there is no leakage, the pressurized fluid from the high-pressure side fills the gap between the contacting surfaces with uniform hydrostatic pressure. The hydrostatic fluid pressure will deform the contacting surfaces and change contact pressure distribution [23]. When the fluid is pressurized from the high-pressure side of an initially sealed system, leakage is observed when the sealed fluid pressure p f equals the maximum contact pressure [24]. The leakage threshold under this condition is governed by the initial geometries of the contact surfaces, the loading condition, and the sealed fluid pressure.
The contact mechanics problem between an elastic bisinusoidal and a rigid flat surface is probably the most simple one that may be used to show the coalescence of the contact patches and the disappearance of the percolation path. Hence, it serves as an important model problem involving the mechanics governing the leakage threshold, which has also been studied in the previous works [4,5,12], without considering the fluid load. In this work, we will revisit this problem and extend it to include the effect of the sealed fluid pressure at the high-pressure side. The main research objective here is to investigate the variability of leakage inception with the sealed fluid pressure in initially sealed contacts.

Problem Setup
The present study is devoted to the linear elastic contact problem between a non-planar smooth surface in contact with a rigid flat surface. The elastic surface has a initial bisinusoidal shape, i.e., with (x, y) ∈ [− ∕2, ∕2] 2 , and is the side length of the square solution domain. The geometry parameter Δ is chosen as Δ = 0.01 in the current study to satisfy the small slope requirement of the half-space approximation, while it is noteworthy that the solution in the current study is independent from Δ∕ when normalized with Δ∕ [25]. The shape parameter controls the initial gap size at the saddle point (x, y) = (0, 0) , and the gap equals Δ under zero loading condition. The initial geometry g 0 (x, y) reduces to a 1D cosine profile when = 0 . When = 2 , the geometry depicted in Eq. (1) is the bi-sinusoidal geometry studied in [4,5,12] under a rotation and a translation operation. In the current study, the shape parameter is chosen as = 0.1, 0.2, 0.3, 0.4 and 0.5 , as shown in Fig. 1.
The non-planar smooth surface g 0 (x, y) is brought into contact with a rigid flat surface under the total load P t , and the average applied pressure is p 0 = P t ∕S 0 , with S 0 = 2 the area of the solution domain. The solid-solid contact problem is first studied with the DCD-FFT method [26] (which is a modification of the DC-FFT method with duplicated padding), in which the periodic boundary condition is applied at y = ±0.5 . With the percolation path being sealed by the coalescence of the two contact patches under constant total load condition, the sealed fluid is filled in from the high-pressure side x = −0.5 with its pressure equaling p f . Under the elevated p f condition, the sealed fluid opens up the solid-solid contacts between the surfaces, allowing for the fluid front to move forward to the point that the percolation path is revealed once more, and the leakage is onset. The value of the sealed fluid pressure when the high-and low-pressure sides are first connected is defined as the critical pressure of leakage p b , and the critical solid-solid contact area when p f = p b is denoted as S b . It is worth to mention that the S b defined in the current study is not the same as the critical area defined in the studies of [5,12], in which the deformation of contacting surfaces by the sealed fluid pressure is ignored. The critical contact area under the dry condition is defined as S * in the current study and the corresponding average applied pressure when the solid-solid contact area S equals S * is defined as the critical pressure of dry contact, with the notation as p * 0 . The values of p * 0 are acquired by conducting the contact mechanics simulation under various p 0 conditions, and checking the geometry connectivity of the high-and lowpressure sides with a connected component labeling algorithm [23,27]. A binary connectivity identifier is defined as 1 when the high-and low-pressure sides are connected, and 0 otherwise. A bisection method is utilized to narrow down the range of the average applied pressure p 0 based on the connectivity identifier, and the simulation is stopped when the range of p 0 is smaller than a specified tolerance = 10 −3 p 0 , and the critical pressure of dry contact is taken as the lower endpoint of the range of p 0 .
The critical pressure of leakage p b is acquired with a boundary element method developed in [23]. The method is a further development of FFT-accelerated method [28] including the hydrostatic load at the contact interface. The mesh size is chosen as 1024 by 1024. The contact plane geometry tolerance is set to 10 −5 Δ , and the total load tolerance is set to 10 −9 p 0 for both the dry contact and the cases with the hydrostatic load.
It is worth to emphasize that the periodic boundary condition is only applied at the y direction in the DCD-FFT method. It differentiates the current study from the previous studies on the bi-sinusoidal surface [4,5,12], in which the periodic boundary condition is applied to both the x and y directions. The y-only periodic boundary condition leads to contact pressure concentration when the contact patches touch the corner of the solution domain due to the flat punch effect, and the contact pressure at the corner is unbounded with large value of shape parameter . So is set to be ∈ (0.0, 0.5] in the current study to avoid the contact pressure singularity under high loading conditions.

Results and Discussions
The merging of the two contact patches occurs continuously for non-adhesive solids, without considering the surface deformation by the sealed fluid pressure [12]. Both the critical pressure p * 0 and the critical area S * under dry contact conditions depend on the shape of the merging contact patches. In this section, the critical quantities p * 0 and S * are first studied for the initial geometries depicted in Eq. (1) and Fig. 1, followed by the analysis of the critical quantities p b and S b with hydrostatic load in consideration. The adhesion force between the contacting surfaces is ignored in the simulation, and the surface tension of the fluid is neglected. The surface stress orthogonal to the contact interface p is referred to as surface traction in the discussion below, following the naming convention in [29]. The surface traction p is the contact pressure under dry contact conditions. For the cases with hydrostatic load, the surface traction p equals the fluid pressure in the fluid-filled regions, and it is the contact pressure at the solid-solid contact regions. Both the force and pressure in the current study are normalized with a pressure scale factor p s , where E is the Young's modulus of the elastic surface, and is its Poisson ratio.  and Xu et al. [25] have shown that there exists a 2/3 power law between the contact area ratio S∕S 0 and the average applied pressure p 0 if the contact patches are semiellipses. The contact area ratio deviates from the 2/3 power law with higher p 0 , and it increases as a concave function of the average applied pressure. There is an inflection interval when p 0 → p * 0 , in which the contact area ratio becomes a convex function of the average applied pressure, as shown in Fig. 3. This inflection interval is firstly discovered by Yastrebov et al. [5] for the isotropic bi-sinusoidal surface, and being confirmed by Xu et al. [25] for the anisotropic bisinusoidal surfaces. The contact area ratio transits back to the concave function of the average applied pressure after p 0 passes p * 0 , and the two contact patches merge together into a dumbbell shape. The size of the dumbbell contact patch, as shown in Fig. 2c, is further increased with elevated p 0 .

Critical Quantities Under the Dry Contact Condition
The relationship between the contact area ratio S∕S 0 and average applied pressure p 0 ∕p s , as shown in Fig. 3   The power laws Eq. (5) give zero critical pressure and critical contact area ratio under dry contact condition when = 0 , this is coinciding with the 1D cosine contact problem, for which there is no percolation paths connecting the high-and low-pressure sides when the loading equals zero.
Substituting the power laws Eq. (5) into Eq. (4) gives The coefficient C can be fitted as with The fitting coefficients c 0 to c 3 are independent from the shape parameter , as shown in Fig. 6 for p 0 ∕p * 0 ∈ (0, 2.0] . The fitting Eq. (6) agrees well with the simulation results when p 0 is closed to p * 0 , and the average and the maximum relative error of Eq. (6) are shown in Table 1 for different ranges of p 0 ∕p * 0 . While Eq. (6) gives a poor estimation for the solid-solid contact area S when p 0 → 0 , the relative error decreases with increased applied load, and the maximum error is less than 2.5% when p 0 ∕p * 0 > 0.5. Equation (6), with the universal coefficients c 0 to c 3 in Eq. (8), suggests that the relationship between the contact area and the average applied pressure, when scaling with their critical values under the dry conditions, is weakly dependent on the initial gap shapes. Utilizing Eq. (6) can significantly simplify the computation of the solid-solid contact area at different loading conditions, and it is of a nature to scale the surface traction p with p * 0 . However, further analysis has shown that the surface traction p∕p * 0 at the saddle point (x, y) = (0, 0) , after the percolation path is closed, depends on the shape parameter (Fig. 7). Since the saddle point has the minimum non-zero surface traction after the percolation path closes, its value is closed relative to the hydrostatic pressure required to open up contact. Therefore, the critical pressure of leakage p b still depends on the initial gap shape. Figure 8 shows the surface traction distribution p∕p s for the initial geometry g 0 (x, y) depicted in Eq. (1), with the shape parameter = 0.4 . The average applied load is p 0 = 0.04p s in Fig. 8a-c, and the percolation path between the high-and low-pressure sides is disconnected under the dry contact condition. The sealed fluid enters from the low-pressure side x = −0.5 with its pressure equaling p f , and the value of p f is p f ∕p 0 = 2.5 × 10 −5 in Fig. 8a, p f ∕p 0 = 0.285 in Fig. 8b, and p f ∕p 0 = 0.6125 in Fig. 8c. The fluid from the high-pressure side filled the gaps between the two contacting surfaces  before being blocked by the contact patch. The surface traction distribution in the contact patch and the sealed fluid pressure together determined the fluid front before the contact is open, and the fluid front of the sealed fluid proceeds to the low-pressure side with increased p f . Since the total load P t is constant in Fig. 8a-c, the maximum solid-solid contact pressure decreases with increased p f , and the contact area S is decreased at the same time due to the proceeding of the sealed fluid. The loss of the solid-solid contact happens at both the high-and low-pressure sides, because of the global response of the elastic displacement with the surface traction distribution [23]. The contact patch has a dumbbell shape which is symmetric with x = 0 when p f = 0 , but loses the symmetry when the sealed fluid pressure increases, because the hydrostatic pressure is only applied at the high-pressure side. Since the neck position of the connected contact patches ( y = 0 ) has the smallest surface traction, the contact length at y = 0 shrinks the fastest with increased p f , and the contact is open along with it when the sealed fluid pressure reaches the critical pressure of leakage p b .

Critical Quantities with Hydrostatic Load
In the current study, the critical pressure of leakage p b is defined as the value of the sealed fluid pressure when the high-and low-pressure sides are first connected, so it can only be uniquely determined when p 0 > p * 0 . A normalized excess pressure p es is defined as and the critical pressure of leakage p b is represented as the function of p es for p es ≥ 0 , as shown in Fig. 9a. When p es = 0 , the two contact patches touch each other only at the saddle point (x, y) = (0, 0) , any infinitesimal fluid pressure can open up the contact and p b ≈ 0 . The ratio between p b and the average applied pressure p 0 firstly increases with p es , and slightly decreases with further growth of p es . The The critical pressure (a) and the critical contact area of leakage (b), as the function of the excess pressure defined in Eq. (9) decrease of p b ∕p 0 with the average applied pressure p 0 is consistent with the line contact problem studied in [30]. However, the critical pressure of leakage p b is always larger than the average applied pressure p 0 in the line contact problem, because the sealed fluid pressure must be increased to the level that all the contact area disappears to reveal the percolation path between the high-and low-pressure sides.
On the other hand, the sealed fluid only needs to open up the contact at the saddle point to create leakage for the geometries Eq. (1), so the critical pressure of leakage p b can be smaller than p 0 in the current study. Another consequence of the existence of the saddle point is that the solid-solid contact area does not need to decrease to zero to initialize the leakage. The critical contact area S b , which is the solid-solid contact area when the hydrostatic load opens the contact, is non-zero as depicted in Fig. 9b. The critical contact area S b equals S * when p es = 0 for all geometries, since any infinitesimal fluid pressure increases will lead to contact open at the saddle point when p 0 = p * 0 . When p 0 > p * 0 , the solid-solid contact area S is larger than S * , and it decreases with the increasing sealed fluid pressure p f . The remaining solid-solid contact area when the contact is opened is generally higher than S * in the current study, and it increases with the shape parameter . The curvatures of g 0 (x, y) at the two crest positions ( (x = 0, y = ±0.5 ) ) are and y increases with . The higher curvature of the initial gap results in lager maximum contact pressure and steeper contact pressure distribution under the same total load condition, and the fluid front is harder to proceed with elevated hydrostatic pressure [30].
The critical pressure of leakage p b illustrated in Fig. 9a is fitted into an empirical relationship: The relative error generated by Eq. (11) is 0.86% in average and maximize at 3.45%.
To further the understanding of the critical pressure ratio p b ∕p 0 , a dimensionless parameter Λ is introduced as the critical pressure factor as follows [30]: with S f the flooded area at the high-pressure side when the sealed fluid pressure p f = 0 and the average applied pressure p 0 > p * 0 . Due to the symmetry of the current problem under the dry contact condition, the S f can be written as The critical pressure factor Λ is plotted against the excess pressure p es in Fig. 10, and it is fitted into an empirical relationship, viz.
The fitting coefficients in Eq. (14) are independent from the shape parameter . Substituting Eq. (14) into Eq. (12) allows one to calculate the critical pressure ratio p b ∕p 0 with S and p es , which are quantities available from the dry contact simulation. The relative error of p b ∕p 0 generated from Eq. (14) is 1.48% in average, when comparing with the simulation results, and the maximum relative error is 4.61%.

Conclusions
A contact problem between a geometrically anisotropic bisinusoidal surface and a rigid flat surface is revisited in the current paper. The critical pressure and the critical solid-solid contact area are of main research interest since they are closely related to the leakage inception. The critical quantities are first studied without considering the hydrostatic load of the pressurized fluid. The relationship between the solid-solid contact area and the average applied load acquired from the current study qualitatively agrees with the previous studies from other researchers, and the detailed transition regime of the contact area is also observed. Instead of studying the asymptotic solution of this contact problem, a generalized form of solution Eq. (3) is constructed based on Westergaard's solution of the 1D cosine contact. By introducing a partition function (p 0 ; ) , 5.14 . Fig. 10 The critical pressure factor, as the function of the excess pressure Eq. (3) can fit the contact area with high accuracy for a wide range of shapes and loading conditions. Moreover, it is found that the relationship between the contact area and the average applied pressure, when normalized with their critical quantities, is independent of the shape parameter.
The critical pressure and critical contact area for leakage are further studied by first applying a load to close the saddle point, and then filling the contact interface with fluid from the high-pressure side. The contact is lost due to increased sealed fluid pressure, and it was discovered that the sealed fluid pressure required to open up the contact and the remaining solid-solid contact area before the revealing of the percolation path can be represented as a function of excess pressure, which is defined as the difference between the applied pressure and the critical pressure under dry contact conditions. On top of that, a dimensionless parameter, called the critical pressure factor, is defined based on the previous studies from the same authors. It is shown that the critical pressure factor, when represented as a function of the excess pressure, is independent from the shape of the initial geometries. Utilizing the critical factor allows one to acquire the critical pressure for leakage using the quantities from the dry contact simulations, which are well studied and easily accessible from the previous research.

Appendix A. The relation between the contact area ratio and the average applied pressure under dry contact condition
The contact area S is a continuous function of the average applied pressure p 0 for non-adhesive elastic contact [12], and it depends on the geometries of the contact surfaces. Hence, the ratio of the contact area to the nominal area can be formulated as where (the continuous function) f (p 0 ; ) should satisfy that: (A.2) f (0; ) = 0, Fig. 11 Top row: the initial elastic surface g 0 (x, y) with the shape parameter = 0.5 (a), the surface g x (x) depicted in eq. (A.4) (b), and the surface g y (y) depicted in eq. (A.5), with = 0.5 (c). Bottom row: the schematic view of the contact region for the initial geom-etry depicted in (a)-(c) contacts with a flat rigid surface. The contact region is colored in red, and the initial elastic surface is g 0 (x, y) in (d), g x (x, y) in (e) and g y (x, y) in (f) for any value of (the shape parameter that controls the initial gap size at the saddle point (x, y) = (0, 0) ) because the contact area is zero when there is no applied load.
In the current problem, the initial gap g 0 (x, y) in eq. (1) can be written as: with and When → 0 , g(x, y) → g x (x) , the contact area S is a rectangular strip of shape 2l cx × , as shown in Fig. 11 (e). The values of l cx , according to the 1D Westergaard's solution [31], is: with Therefore,