On the loading and unloading of metal-to-metal seals: A two-scale stochastic approach

Abstract

During operation, the mating surfaces of a metal-to-metal seal typically undergo significant plastic deformation, which in turn can have beneficial effect on its performance. In previous studies, it has, for instance, been shown that plastic deformation can provide for better sealing during unloading. Those studies did, however, only consider flow through unrealistically small domains. Therefore, it is possible that this might be a size effect, which would not be apparent in a real situation with a much larger domain. In this paper, we develop a model which can handle real-sized seal domains at the same time as fine details of the surface topography. More precisely, we construct a two-scale model, in which the global scale represents the seal domain and where the influence of the fine details at the local scale are represented by a stochastic element. By means of this stochastic two-scale model, we show that the beneficial effect associated with the plastic deformation persists also when real-sized seal domains are considered.

Keywords

Contact mechanics load cycle Reynolds equation seals stochastic two-scale modelling

Introduction

The metal-to-metal seal is an example of an application in which the performance is highly influenced by the evolution of the surface topography and, in particular, by the plastic deformation occurring during loading cycles. Kadin et al.¹ presented numerical simulation results showing that plastic deformation leads to a more intimate contact between the surfaces and this may improve the sealing performance. Experimental evidence supporting this can be found in the results by Murtagian et al.² Despite this, little is known about the mechanism through which plastic deformation affects the leakage. Therefore, a study of the effect of plastic deformation on sealing performance was undertaken by the present authors in Pérez-Ràfols et al.³ That study lead to an interesting observation: after loading the seal up to a certain applied load, one could reduce substantially the applied load without observing any significant increase in leakage. This is relevant to the robustness of the metal-to-metal seal, as it gives an indication of the amount of load a seal can lose before it fails. The main limitation of the work in Pérez-Ràfols et al.³ is that it considers unrealistically small domains.

In this work, we attempt to answer the question: will this behaviour still be prevalent when a domain large enough to cover the size of a realistic seal is considered? To this end, we extend the model presented by Pérez-Ràfols et al.⁴ to allow for studying the contact mechanics and leakage of seals during unloading as well. As in Pérez-Ràfols et al.,⁴ the present model couples a classical two-scale type of formulation with an explicit consideration to the stochastic nature of the roughness at the local scale. A novelty is a stochastic construction for the permeability, which considers the complete loading–unloading cycle.

The rest of the paper is organised as follows. First, a fully deterministic model (in opposition to a two-scale model) is introduced to serve as a reference. Then, the two-scale model is presented in the following section and its accuracy is evaluated in the Appendix. In ‘The two-scale stochastic model’ section, the two-scale stochastic model is presented and, in ‘Results and discussion’ section, the model is used to answer the question posed.

The deterministic model

A deterministic model is used to assess the accuracy and validity of the outcome of the stochastic two-scale model, which is the main development in this paper. This deterministic model is based on the one presented in Sahlin et al.⁵ and later in Pérez-Ràfols et al.³ It is, therefore, only briefly summarised here.

In this work, we assume that the deformation caused by fluid-structure interaction is small as compared to the deformation caused by direct mechanical contact between the two surfaces. It follows from this assumption that it is possible to separate the model into two distinct physics, i.e.

A contact mechanics model to determine the elasto-plastic deformation due to a total applied load W, and the elastic spring back when releasing the load.

A fluid flow model to compute the amount of fluid that flows through the rough aperture between the two deformed surfaces obtained from the contact mechanics model for a given a pressure drop driving the flow.

In what follows in this section, we present the contact mechanics model, then we explain how unloading is modelled and finally we present the fluid flow model. The same rectangular solution domain, Ω, defined as

Ω : = {x | 0 < x_{1} < L_{1}, 0 < x_{2} < L_{2}}

(1)

will be used when comparing the different models. For the discretisation of the deterministic model, a uniform rectangular mesh with

N_{1} \times N_{2}

elements of size

Δ x_{1} \times Δ x_{2}

is adopted. The coordinate for the node

({x_{1}}_{o}, {x_{2}}_{p})

is thus given by

{x_{1}}_{o} = o Δ x_{1}

and

{x_{2}}_{p} = p Δ x_{2}

, with

Δ x_{1} = L_{1} / (N_{1} - 1)

and

Δ x_{2} = L_{2} / (N_{2} - 1)

The contact mechanics model

The contact mechanics model utilised here is the same as the one presented in Sahlin et al.⁵ Therefore, only a summary of the model will be given. The model assumes that the behaviour of the surface can be described by a half-space in which friction is neglected. We can thus utilise the Boussinesq solution to model the relation between the contact pressure and the elastic deformation. Moreover, a linear elastic perfectly plastic material model is adopted. More precisely, when the contact pressure at a certain point exceeds the hardness of the softer material, H, the point floats at the surface and the pressure remains equal to the hardness. The hardness of the material is taken as 2.8 times the yield limit, see Johnson.⁶ The desired output for the contact mechanics model (in the context of this work) is the gap between the deformed surfaces, h, defined as

h (x) = h_{1} (x) + g_{00} + u_{e} (x) + u_{p} (x)

(2)

where g₀₀ is the rigid body separation, h₁ is the undeformed gap between the contacting surfaces and u_e and u_p are the elastic and plastic deformations, respectively. According to the Boussinesq solution, the elastic part of the deformation (u_e) is computed from the contact pressure, p_c, as

u_{e} (x) = \frac{1}{E^{*}} \int_{Ω} \frac{p_{c} (x')}{\sqrt{{(x_{1} - x'_{1})}^{2} + {(x_{2} - x'_{2})}^{2}}} d x'

(3)

where

E^{*}

is the reduced elastic modulus, defined as

E^{*} = (\frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - ν_{2}^{2}}{E_{2}})^{- 1}

(4)

Here E₁, E₂, ν₁ and ν₂ are the elastic modulus and Poisson ratios of the contacting surfaces. Since contact is defined by h = 0 and $0 < p_{c} \leq H$ and since h > 0 and p_c = 0 where there is no contact, we adopt the following Kuhn–Tucker complementary condition to model the mechanical contact between the surfaces.

\begin{matrix} h (x) > 0, p_{c} (x) = 0, x \in Ω_{c}, \\ h (x) = 0, 0 < p_{c} (x) \leq H, x \notin Ω_{c}, \end{matrix}

(5)

where Ω_c is the part of Ω where both surfaces are in contact. Also, since we want to obtain the deformed gap, h, for different nominal contact pressures, W, we impose

W = \frac{1}{A_{n}} \int_{Ω} p_{c} (x) d x

(6)

where A_n is the nominal contact area.

As in Sahlin et al.,⁵ the complementarity problem, (5), is numerically solved employing the variational principle, formulated for the total complementary potential energy and by using DC-FFT method to accelerate the computation of the convolution in (3).

Loading and unloading

From the way plasticity is implemented, it follows that the predicted unloading of the seal is equal to any subsequent reloading, all the way up to the load where the unloading previously started. This is obvious from the fact that $p \leq H$ and that unloading will only reduce the pressure. Note that, although being a result of our simple way of treating plastic deformation, this behaviour is actually not far from reality, as observed in Kadin et al.⁷ This observation allows us to use (5) for both loading and unloading. During loading, the original gap (h₁) is used when defining the deformed gap (h) in (2). For the unloading, u_p is set to $u_{p}^{max}$ , which is the plastic deformation occurred at the maximum load reached before unloading started. Notice that despite plasticity is included in the model, no extra plastic deformation besides $u_{p}^{max}$ will occur during unloading. The same will then apply to any subsequent reloading up to the maximum load reached during the previous loading cycle.

The fluid flow model

In this section, we present the model for the flow through the rough aperture between the deformed surfaces (h). This is a Reynolds equation based model that can be applied to study the flow situations for which the classical thin film assumptions are valid. If the fluid is Newtonian, incompressible and iso-viscous, it can be expressed as

\nabla \cdot (h^{3} \nabla p_{f}) = 0

(7)

where p_f is the fluid pressure. For convenience in the later derivation, p_f is considered dimensionless (the pressure is scaled by the pressure drop,

Δ P

, between the inlet and the outlet of the seal). This particular form of the Reynolds equation is solved together with periodic boundary conditions in x₂-direction and unitary pressure drop in x₁-direction, i.e.

p_{f} (x_{1}, 0) = p_{f} (x_{1}, L_{2})

(8)

p_{f} (0, x_{2}) = 1

(9)

p_{f} (L_{1}, x_{2}) = 0

(10)

Note that the unitary pressure drop comes from the non-dimensionalisation of p_f. In order to obtain the solution for a particular case, one simply needs to scale the non-dimensional pressure, p_f with $Δ P$ .

Once (7) has been solved for the fluid pressure, p_f, the leakage, Q_d, can be computed. In general, the leakage equals the mass flow integrated over the outlet boundary. However, due to the periodicity imposed in x₂-direction, the leakage can also be evaluated as the mass flow integrated over the whole domain, i.e.

Q_{d} = \frac{1}{L_{1}} \frac{Δ P}{12 η} \int_{Ω} h^{3} \frac{\partial p_{f}}{\partial x_{1}} d x

(11)

This turns out to be a better option for the numerical evaluation since the pressure gradient at the boundary is close to zero and the numerical error would thus be greater if the integration would be carried out over the outlet boundary only. In (11), η is the (constant) dynamic viscosity of the fluid.

The two-scale model

The two-scale model presented in this work is based on the same ideas as the model presented in Pérez-Ràfols et al.⁴ Therefore, detailed explanations will only be given when the two models differ from each other.

When presenting the two-scale model, we will follow the same structure we used to present the deterministic model. This means that we will first present the two-scale formulation of the contact mechanics model, then we will discuss the particularities of the modelling of loading and unloading when two scales are considered and finally we will present the two-scale formulation of the flow model. A schematic illustration of the outcome of a numerical simulation, obtained by means of the suggested two-scale model, the two-scale model is shown in Figure 1.

Figure 1.

A schematic presentation of the flow situation (left) and the contact pressure variation (right) at both the global scale and the local scale, depicted in the inserts. In the left figure, the red arrows represent the direction and magnitude of the flow and the yellow to blue colormap represents the fluid pressure (p_f) at each node of the global domain. Similarly, the pattern coloured in red in the insert depicts the local scale flow, whilst the corresponding fluid pressure distribution is represented by the colormap in yellow and blue. In the right figure, the nominal contact pressure ( ${\bar{p}}_{c}$ ) over each element of the global-scale domain is depicted by the yellow to blue colormap. The insert shows the contact pressure distribution (p_c) at the local scale.

The two-scale contact mechanics model

The two-scale contact mechanics model is built upon two main assumptions, i.e. (i) using a low-pass filtered surface as an input for the deterministic contact mechanics model permits obtaining a (sufficiently accurate) coarse representation of the global-scale distributions of contact pressure and average separation, and (ii) the local and global scales can be coupled by approximating the nominal contact pressure and average separation at the local-scale by the global scale pressure ${\bar{p}}_{c}$ and separation $\bar{h}$ . The latter assumption requires that the deformation caused by the pressure in the surrounding area only have a marginal influence on the shape of the local-scale topography. In Pérez-Ràfols et al.,⁴ these two assumptions were justified by applying a Gaussian filter to the equations governing the contact mechanics problem and by invoking St. Vennant's principle. In this work, a numerical justification of the two-scale contact mechanics model is given in the Appendix.

Based on the two previous assumptions, the two-scale contact mechanics model can be built following the steps, summarised in Figure 2. We start with a large measurement, which has been processed to avoid any artefact coming from the measurement. The height data are represented on a fine grid with $N_{1} \times N_{2}$ elements of size $Δ x_{1} \times Δ x_{2}$ (the grey grid in Figure 3). These height data are then filtered using a low-pass filter with a sharp cutoff at a frequency ω_c. Once filtered, these height data are interpolated onto a coarser grid, with grid nodes $({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ , where ${\tilde{X_{1}}}_{k} = k Δ \tilde{X_{1}}$ and ${\tilde{X_{2}}}_{l} = l Δ \tilde{X_{2}}$ (see black grid in Figure 3). We here define a coarsening factor $c = Δ \tilde{X_{1}} / Δ x_{1} = Δ \tilde{X_{2}} / Δ x_{2}$ . The coarser grid has thus $N_{1} / c \times N_{2} / c$ elements. We can now choose $ω_{c} = 2 π / (8 Δ \tilde{X}) = (2 π / (8 c Δ x)$ ). According to previous studies,⁸ this choice permits obtaining meaningful output from the contact mechanics model.

Figure 2.

Flow chart of the algorithm. On the left, the part corresponding to the computations on the local scale. On the right, the global scale contact mechanics part, described in ‘The two-scale contact mechanics model’ section. The yellow boxes correspond to intermediate outputs.

Figure 3.

Representation of the different domains and grids used in the contact mechanics model. The grey lines with nodes $({x_{1}}_{m}, {x_{2}}_{n})$ represent the grid in which the height data are measured. The global-scale contact mechanics problem is solved on the black grid, with nodes $({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ . Finally, the red grid, with nodes $({X_{1}}_{i}, {X_{2}}_{j})$ , is where the global-scale flow problem is solved. The domains $ω_{i + 1 / 2, j}^{1}$ and $ω_{i, j + 1 / 2}^{2}$ , presented in (12), are also depicted. The values of ${\bar{p}}_{c}$ and $\bar{h}$ are also averaged over these domains.

The deterministic contact mechanics model is used to obtain coarse grid representation of the contact pressure ${\tilde{p}}_{c} = {\tilde{p}}_{c} ({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ and the separation $\tilde{h} = \tilde{h} ({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ . Then, these coarse representation of contact pressure and separation are use to obtain representations on an even coarser grid $({X_{1}}_{i}, {X_{2}}_{j})$ , where ${X_{1}}_{i} = i Δ X_{1}$ and ${X_{2}}_{j} = j Δ X_{2}$ (see red grid in Figure 2). This is done by assigning to each point $({X_{1}}_{i}, {X_{2}}_{j})$ a nominal pressure, ${\bar{p}}_{c}$ , which is the average of the pressures, ${\tilde{p}}_{c}$ at the surrounding nodes $({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ . Similarly, an average separation $\bar{h}$ is defined. More precisely, the average is done over the local-scale domains $ω_{i + 1 / 2, j}^{1}$ and $ω_{i, j + 1 / 2}^{2}$ , around each pair of points $[({X_{1}}_{i}, {X_{2}}_{j}), ({X_{1}}_{i + 1}, {X_{2}}_{j})]$ and $[({X_{1}}_{i}, {X_{2}}_{j}), ({X_{1}}_{i}, {X_{2}}_{j + 1})]$ , respectively,

\begin{matrix} \begin{matrix} ω_{i + 1 / 2, j}^{1} = {x | {X_{1}}_{i} < x_{1} < {X_{1}}_{i + 1}, {X_{2}}_{j - 1 / 2} < x_{2} < {X_{2}}_{j + 1 / 2}} \\ ω_{i, j + 1 / 2}^{2} = {x | {X_{1}}_{i - 1 / 2} < x_{1} < {X_{1}}_{i + 1 / 2}, {X_{2}}_{j} < x_{2} < {X_{2}}_{j + 1}} \end{matrix} \end{matrix}

(12)

For clarity, these two domains are depicted in Figure 3. The reason for the particular choice of domains will be apparent in ‘The two-scale flow model’ section when discussing the two-scale flow model.

We can now solve the problem at the local scale using these domains following the procedure summarised in Figure 2. The contact mechanics model can be used to obtain the local-scale distributions of h and p_c on the local-scale domains $ω_{i + 1 / 2, j}^{1}$ and $ω_{i, j + 1 / 2}^{2}$ . The inputs to the local-scale problems are the material parameters E_i, ν_i and H, the undeformed gap between the contacting surfaces h₁. In order to decouple the solutions on the global and local scale, we solve for h and p_c for a wide range of nominal pressures ${\bar{p}}_{c}$ , corresponding to a wide range of separations $\bar{h}$ . Then, the results for the particular value of ${\bar{p}}_{c}$ obtained from the global scale solution, the result can be obtained via interpolation.

Loading and unloading in the two-scales model

The plastic deformation at both the local and global scales is accounted for in the same way as in the deterministic model (i.e. as described in ‘Loading and unloading’ section). However, since the plastic deformation is considered only point-wise, i.e. by assuming that for each point for which $p_{c} (x) = H$ ‘float’ to the contact surface, one can argue that the model only considers plastic deformations at frequencies of the order of $Δ x$ only. This implies that the plastic deformation predicted to occur at the global scale, in the two-scale model, will have another (lower) frequency than the plastic deformation predicted by the deterministic model. Together with the high frequency contribution to plastic deformation taking place at the local scale, one can argue that the two-scale model represents a more realistic plastic behaviour than the deterministic one. Nevertheless, we show in ‘Justification of the two-scale contact mechanics model’ section in the Appendix that the plastic deformation is small at the global-scale, which means that the results coming from the deterministic- and the two-scale model are comparable. Note that this is an implicit requirement for using the deterministic model as a reference representation of the contact mechanics in metal-to-metal seals.

The two-scale flow model

The two-scale flow model is constructed following the Heterogeneous Multi-scale Method (HMM) approach. In Engquist et al.,⁹ a comprehensive review of the method is given and in de Boer et al.,^10,11 it is applied to a related fluid flow problem including fluid structure interaction. The construction is based on the fundamental fact that the Reynolds equation is a mass flow continuity equation. Accordingly, in a discrete form, the global scale flow is described by imposing mass balance at each global-scale point $({X_{1}}_{i}, {X_{2}}_{j})$ (see the red grid in Figure 3). Indeed, let $j_{i + 1 / 2, j}^{1}$ represent the flow in X₁-direction (as indicated by the superscript) between the global-scale points $({X_{1}}_{i}, {X_{2}}_{j})$ and $({X_{1}}_{i + 1}, {X_{2}}_{j})$ and define other flows analogously. Mass balance at $({X_{1}}_{i}, {X_{2}}_{j})$ can then be expressed as

- j_{i - 1 / 2, j}^{1} + j_{i + 1 / 2, j}^{1} - j_{i, j - 1 / 2}^{2} + j_{i, j + 1 / 2}^{2} = 0

(13)

According to Darcy's law, the flow $j_{i + 1 / 2, j}$ is given by as

j_{i + 1 / 2, j} = \frac{1}{η} k_{i + 1 / 2, j}^{1} δ p_{f}^{1}_{i + 1 / 2, j}

(14)

where

k_{i + 1 / 2, j}^{1}

is the permeability of the gap h on

ω_{i + 1 / 2, j}^{1}

in the X₁-direction and

δ p_{f}^{1}_{i + 1 / 2, j} = p_{f}_{i + 1, j} - p_{f}_{i, j}

(15)

is the corresponding pressure difference. The other flows are computed in a similar manner.

Now, in order to obtain the permeability, we use the fluid flow model at the corresponding local-scale domain. Therefore, we solve Reynolds equation (7) imposing a unitary pressure drop in the direction in which the permeability is computed and periodic boundaries in the direction perpendicular to it. This requires mirroring the local-scale representation of h. The choice of boundary conditions is not trivial, as it must be consistent with the definition of the global-scale model. For a more detailed discussion of the actual choice, the reader is referred to Pérez-Ràfols et al.⁴ Once the Reynolds equation is solved, the permeability, $k_{i + 1 / 2, j}^{1}$ , can be computed as

k_{i + 1 / 2, j}^{1} = \frac{1}{12 Δ X_{1}} \int_{ω_{i + 1 / 2, j}^{1}} h^{3} \frac{\partial p_{f}}{\partial x_{1}} d x

(16)

Those readers more familiar with porous media literature will probably notice that the permeability presented here has units of m³ instead of m². This is because of the use of Reynolds equation, which simplifies the problem by integrating over the thin dimension. Once permeability is computed for all global-scale points, (13) can be solved for ${p_{f}}_{i, j}$ and then the leakage $Q_{2 s}$ , which equals the total flow over the global-scale's domain's outlet boundary $\partial Ω_{o}$ , can be computed as

Q_{2 s} = \frac{1}{η} \sum_{\partial Ω_{o}} k_{i + 1 / 2, j}^{1} δ p_{f}^{1}_{i + 1 / 2, j}

(17)

Note here that we have neglected the leakage in x₁-direction caused by a pressure drop in x₂-direction. For the specific surfaces studied in this work, this addition to the leakage is known to be small and this is an acceptable approximation.

The two-scale stochastic model

Once the two-scale model has been defined, we can use it to build a two-scale stochastic model that can be used to cover larger domains. Conceptually, the model is the same as the one presented in Pérez-Ràfols et al.⁴ In this work, we extend it to be able to asses seal performance during unloading.

In the two-scale model presented previously, the global domain is discretised into cells and the local-scale model is used to compute the permeability for each of those cells. However, as any other roughness related parameter, the permeability exhibits the properties of a random variable, following a given distribution. In principle, this implies that any value of permeability is possible. However, in reality only a certain range of values have a reasonably high probability of occurring and we can make use of this fact when constructing a random permeability map based on data from a limited number of local-scale cells.

We will now describe the procedure followed here to obtain the total flow $Q_{2 s}$ , by means of the two-scale stochastic model. First, we generate random permeability maps kⁱ in the X₁- and X₂-directions. Then we solve (13) for the pressure ${p_{f}}_{i, j}$ and finally we can compute $Q_{2 s}$ by means of (17). In order to generate these permeability maps, we start by computing the pressure and separation distributions ( ${\tilde{p}}_{c}$ and $\tilde{h}$ ) corresponding to the coarsened surface height data, i.e. using the grid points $({\tilde{X_{1}}}_{k}, {\tilde{X_{2}}}_{l})$ . The coarse grid pressure and separation distributions are then further coarsened into the points $({X_{1}}_{i}, {X_{2}}_{j})$ , and we obtain maps for nominal pressure ( ${\bar{p}}_{c}$ ) and average surface separation ( $\bar{h}$ ). These maps can then be extended to cover the whole global-scale domain, representing the full seal. This is done by extending the domain Ω by repetition. At this stage, long wavelength pressure variations, originating from e.g. out of flatness or out of roundness, can be added. Notice here that even when the pressure distribution is more or less constant at the global scale, the leakage through an extended domain can be different from the collective leakage through the individual domains that the extended one has been generated from. In Pérez-Ràfols et al.,⁴ it was shown that the flow paths through an extended domain may cover larger distances in the circumferential direction that what would be permitted by a smaller sized domain. Once the maps for ${\bar{p}}_{c}$ and $\bar{h}$ have been generated, they will be used to construct the random permeability maps.

In Pérez-Ràfols et al.,⁴ the permeability maps were constructed by describing the permeability, k, as a log-normal distributed random variable, with the parameters of the distribution being dependent on ${\bar{p}}_{c}$ (or $\bar{h}$ ). There are two reasons for this not being appropriate in this case. The first one is that the larger number of local-scale cells, for this study, revealed that a log-normal distribution is not adequate to describe the variability of k. Note that the log-normal distribution was also ruled out in a recent related work¹² by Waseem et al. The second one relates to the aim of this work. With the pressure study, we investigate leakage and contact mechanics during unloading, and it is therefore necessary that the permeability at a given location and load is correlated with the permeability at the same location but at all loads within the cycle. That is, we want to simulate random realisations of load cycles and not of separate load cases, which was the case in previous study.⁴ Doing otherwise might hide relevant information of the unloading process.

In the rest of this section, we will describe the behaviour of local-scale cells and how we apply this description to the construction of the permeability maps.

Local-scale behaviour and permeability map

Figure 4 depicts loading–unloading cycles computed for four different local-scale cells, selected at random. The results presented in the figure show the variability of the local-scale cells' loading–unloading behaviour. For instance, if we compare the blue (solid) and the purple (dotted) curves, we see a similar response at the initial stage of unloading but that the difference become larger as the unloading proceeds. The yellow (dashed) and orange (dot-dashed) depict a behaviour conceptually different than the blue and the purple, as they both exhibit zero permeability before the loading reaches 900 MPa.

Figure 4.

Loading (rightmost part of each curve) and unloading cycles for four local-scale cells selected at random. All cells have been loaded to a maximum nominal pressure of 900 MPa and subsequently unloaded. The absence of line indicates zero permeability.

Given the large variability of the local-scale cells' loading–unloading behaviour, it is clear that the permeability for every load in the loading–unloading cycle must be correlated with the permeability associated with all the other loads in the cycle. Otherwise, the resulting realisation would not show a physically representative behaviour.

In order to achieve the desired correlation, one would ideally describe the loading–unloading cycle at a local-scale level by means of a mathematical expression, and include the random variability of the parameters included in that expression. Unfortunately, there is no such expression available for the surfaces considered in this work, and developing one with physical meaning would be a matter of another study. Instead, the authors suggest a methodology related to the derivative of kⁱ with respect to ${\bar{p}}_{c}$ (or $\bar{h}$ ). An overview of this methodology is depicted in Figure 5.

Figure 5.

Flow chart of the algorithm for the stochastic realisations described in ‘The two-scale stochastic model’ section. The yellow boxes correspond to intermediate outputs.

In the following, the process for constructing the permeability maps for the loading part of the cycle will be detailed. The permeability map for the highest load in the loading cycle is randomised from an empirical cumulative density function, obtained from the local-scale results for different values of ${\bar{p}}_{c}$ (or $\bar{h}$ ). With this as a starting point, the permeability increases when the load reduces (or as the separation increases). This increase in permeability follows the derivative of kⁱ for a randomly chosen local-scale cell. In ‘Justification of the two-scale contact mechanics model’ section in the Appendix, we justify why this approach gives a more robust solution than starting at a low load. We found that, for higher loads, the relation between kⁱ and ${\bar{p}}_{c}$ captures this increase better than the relation between kⁱ and $\bar{h}$ . As seen in Figure 6, this is because, at higher loads, a very small change in h results in a large change in kⁱ and thus any error in h would change the results significantly. However, when the load approaches zero, the relation between kⁱ and $\bar{h}$ is the preferred one. Here we choose ${\bar{p}}_{c} = 0.01 H$ as a transition point for switching between the two relations. Numerically, it is more efficient to use the relation between $log (k^{i})$ and ${\bar{p}}_{c}$ (or $\bar{h}$ ) when calculating the derivative. The trapezoidal integration rule is then used to obtain the permeability increase for a given reduction of ${\bar{p}}_{c}$ or increase of $\bar{h}$ . In the case when the randomised permeability value at the highest load is zero, we cannot directly follow this procedure. Instead, we select the value of ${\bar{p}}_{c}$ at which the permeability of the chosen local-scale cell first becomes positive as the starting point.

Figure 6.

Permeability in x₁-direction as a function of the average separation (left) and nominal contact pressure (right). The orange color corresponds to loading and the blue one to unloading. The solid black line in the left represents the asymptote ${\bar{h}}^{3} / 12$ .

The permeability maps for the unloading part of the cycle are computed following the same procedure. When computing the unloading starting at a load W^u, we start from the permeability map computed during the loading part at the load $W = W^{u}$ . From there, the unloading curve is computed in the same manner as the loading part.

A final comment is made for when the separation is large. In such a case, the permeability approaches the asymptote ${\bar{h}}^{3} / 12$ . This is shown in Figure 6. In order to avoid unnecessary computations, it is convenient to set a threshold ${\bar{h}}_{0}$ above which the permeability is not computed randomly but defined as $k^{i} = {\bar{h}}^{3} / 12$ instead. According to the results presented in Figure 6, one can safely choose ${\bar{h}}_{0} = 5$ μm, at least for the surfaces used in this study.

Results and discussion

The main goal of this work is to assess whether the beneficial effect associated with the plastic deformation, observed in Pérez-Ràfols et al.,³ persists also when real-sized seal domains are considered. That is, the zone of more or less constant leakage, from the point when the unloading started, can be maintained during the release of a substantial part of the applied load. To test this hypothesis, we consider two different combinations of surfaces resulting in the two undeformed gaps (h₁), depicted in Figure 7. The topography depicted on the left is the gap between a measured turned surface contacting a flat one, and the topography to the right is a computer generated self-affine fractal surface, also contacting a flat one. After measuring the turned surface, a sample consisting of 2000 × 2000 points is selected. It is then low-passed filtered by truncation in the frequency domain. The reason to impose this low-pass filter is to ensure that even the shortest wavelength components are well resolved numerically. More precisely, as in Yastrebov et al.,⁸ we truncate frequencies higher than $1 / (8 Δ x)$ , resulting in eight elements resolving the shortest wavelength. In order to avoid problems caused by the periodic boundary conditions imposed in x₂-direction, we mirror the filtered topography.

Figure 7.

The unreformed gap, h₁, used for the evaluation of the two-scale model. On the left, corresponding to a turned surface (measured) against a flat one and, on the right, a self-affine fractal surface (generated) against a flat one.

The (periodic) self-affine fractal surface is generated using the algorithm described in Putignano et al.¹³ It has a size 1 mm × 1 mm, represented in a grid of size $2^{12} \times 2^{12}$ . It has an rms height of 1 μm, a lower cutoff frequency of $10 / L$ and a higher cutoff frequency of $1 / (8 Δ x)$ . Note that $Δ x$ is different for the two surfaces.

The material properties used are $E_{1} = E_{2} = 206$ GPa, $ν_{1} = ν_{2} = 0.3$ and H = 2.75 GPa. For the contact mechanics algorithm, we use a tolerance for the load balance of $10^{- 3}$ % and a tolerance for the contact plane of $10^{- 10}$ m, see Sahlin et al.⁵ for details. For the turned surface, we choose a local cell size of $132 Δ x \times 132 Δ x$ and a coarsening factor c = 4. The local-scale cell size in the x₁-direction is chosen to match (approximately) the turning grooves seen in Figure 7, and then this size is also used in the x₂-direction. For the fractal surface, the chosen local-scale cell size is $128 Δ x \times 128 Δ x$ . The reader can find a justification of the validity of the model and of the particular choices for the cell size and the coarsening factor in the Appendix.

In order to assess the influence of the seal size, we construct several domains of size $L_{1} \times n L_{2}$ (where L₂ is the size of the sampled surface in x₂-direction). Notice that no extension in x₁-direction is done, as the measured size is already large enough to be considered realistic. Figure 8 shows the leakage, $Q_{2 s}$ , for different domain sizes during loading and unloading. In all cases, the leakage is computed assuming unitary pressure drop and viscosity. It can be noticed how the mean value of the leakage does not change significantly when the domain size in x₂-direction increases. On the other hand, the uncertainty reduces and this is consistent with the observation made in Pérez-Ràfols et al.⁴ for the loading case.

Figure 8.

Leakage per unit length in x₂-direction. On the left, leakage for different domain of sizes $L_{1} \times n L_{2}$ . Solid lines represent the mean value of 500 realisations while the dashed lines correspond to the 95% confidence interval. On the right, the mean value corresponding to n = 4 for loading and various unloading curves. The results in the upper part correspond to the turned surface and the ones in the lower part to the self-affine fractal surface.

In Figure 8, the loading and several unloading curves are also shown. These results confirm the hypothesis that a zone of more or less constant leakage exists also when real-sized seal domains are considered. In Pérez-Ràfols et al.,³ some of the cases studied presented a much more clear transition at which the rapid leakage increase starts. Despite that the results are consistent with the general trends observed. An explanation for this behaviour can be made using the concept of the critical constriction. That is, at sufficiently small total leakage, the pressure drop occurs over a small number of critical constrictions. Because of this, the leakage at high loads is controlled by local features regardless of the size of the domain.

Concluding remarks

The model presented by Pérez-Ràfols et al.⁴ was extended to allow for studying the contact mechanics and leakage of seals during loading and unloading. The present model couples a classical two-scale type of formulation with an explicit consideration to the stochastic nature of the roughness at the local scale. A novelty is a stochastic construction for the permeability, which considers the complete loading–unloading cycle. We have assessed the validity of the model and we conclude that it can provide for trustworthy predictions of leakage during loading and unloading of metal-to-metal seals. Due to its two-scale stochastic nature, it is applicable for studying real-sized seal domains, returning the solution in reasonable times.

The presented model was used to test the hypothesis whether a zone of more or less constant leakage, which starts at the beginning of the unloading and persists during a substantial release of the load, can be observed. The results presented in this work support this hypothesis, which arose while studying the loading–unloading behaviour for small (local) seal domains in Pérez-Ràfols et al.³

Footnotes

Acknowledgements

The authors are grateful to Shell Global Solutions BV for permission to publish.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Shell Global Solutions BV and Swedish Research Council (VR).

Orcid ID

Francesc Pérez-Ràfols .

Appendix

References

Kadin

Kligerman

Etsion

. Unloading an elastic-plastic contact of rough surfaces. J Mech Phys Solids 2006; 54: 2652–2674.

Murtagian

Fanelli

Villasante

et al.

Sealability of stationary metal-to-metal seals. J Tribol 2004; 126: 591–596.

Pérez-Ràfols

Larsson

van Riet

et al.

On the flow through plastically deformed surfaces under unloading: A spectral approach. Proc IMechE, Part C: J Mech Eng Sci 2017. doi.org/10.1177/0954406217690180.

Pérez-Ràfols

Larsson

Lundström

et al.

A stochastic two-scale model for pressure-driven flow between rough surfaces. Proc Roy Soc A: Math Phys Eng Sci 2016; 472: 20160069.

Sahlin

Larsson

Almqvist

et al.

A mixed lubrication model incorporating measured surface topography. part 1: Theory of flow factors. Proc IMechE, Part J: J Eng Tribol 2010; 224: 335–351.

Johnson

. Contact mechanics, The Edinburgh building, Cambridge CB2 2RU, UK: University of Cambridge, 1985.

Kadin

Kligerman

Etsion

. Multiple loading-unloading of an elastic-plastic spherical contact. Int J Solids Struct 2006; 43: 7199–7127.

Yastrebov

Anciaux

Molinari

. From infinitesimal to full contact between rough surfaces: Evolution of the contact area. Int J Solids Struct 2015; 52: 83–102.

Engquist

WEB

Ren

et al.

Heterogeneous multiscale methods: A review. Commun Comput Phys 2007; 2: 367–450.

10.

de Boer

Hewson

Thompson

et al.

Two-scale EHL: Three-dimensional topography in tilted-pad bearings. Tribol Int 2014; 79: 111–125.

11.

de Boer

Gao

Hewson

et al.

A multiscale method for optimising surface topography in elastohydrodynamic lubrication (EHL) using metamodels. Struct Multidiscip Optim 2016; 54: 483–497.

12.

Waseem

Guilleminot

Temizer

. Stochastic multiscale analysis in hydrodynamic lubrication. Int J Numer Method Eng 2017; 112(8): 1070–1093.

13.

Putignano

Afferrante

Carbone

et al.

A new efficient numerical method for contact mechanics of rough surfaces. Int J Solids Struct 2012; 49: 338–343.

14.

Almqvist

Dasht

. The homogenization process of the Reynolds equation describing compressible liquid flow. Tribol Int 2006; 39: 994–1002.

15.

Almqvist

Essel

Fabricius

et al.

Variational bounds applied to unstationary hydrodynamic lubrication. Int J Eng Sci 2008; 46: 891–906.

16.

Almqvist

. Homogenization of the Reynolds equation governing hydrodynamic flow in a rotating device. J Tribol 2011; 133: 021705.