A New Progress in Random Hydrodynamic Lubrication for Movable Non-Rotational Curvilinear Biosurfaces with Phospholipid Bilayers

Krzysztof Wierzcholski ^*

WSG Bydgoszcz University Garbary 2, Bydgoszcz, Poland

* Correspondence: Krzysztof Wierzcholski

Academic Editor: Hossein Hosseinkhani

Received: October 08,2020 | Accepted: April 13,2021 | Published: June 02,2021

Recent Progress in Materials 2021, Volume 3, Issue 2, doi:10.21926/rpm.2102023

Recommended citation: Wierzcholski K. A New Progress in Random Hydrodynamic Lubrication for Movable Non-Rotational Curvilinear Biosurfaces with Phospholipid Bilayers. Recent Progress in Materials 2021; 3(2): 023; doi:10.21926/rpm.2102023.

© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

The subject of this paper reflects new progressive probabilistic estimation methods of lubrication solutions for movable, non-rotational, arbitrary, living biological surfaces. The superficial layer of living bio-surfaces and the overlaying phospholipid bilayer membrane, which is several nanometers thick, shape the height of the gap between two co-operating, curvilinear, non-rotational bio-surfaces [1,2,3]. Based on numerous recent studies conducted with the use of the Atomic Force Microscope (AFM), it can be established that the phospholipid layer and the superficial layer are not static [4]. The dimensions of the surface coated with phospholipids (PLs) as well as the total gap height ε_T between the co-operating surfaces is subject to several stochastic changes in relation to the nominal mean value. Abrupt changes can be caused by micro-vibrations, discrete load cumulated on the friction nod, or changes in the roughness geometry of the co-operating surfaces. Another cause of random variations in gap height is the genetic and volumetric growth of live cells on the living surface with a phospholipid layer. Such small random changes are significant [5,6,7]. Numerous previous experimental studies on the scope of influence of the phospholipid membrane on the hydrodynamic process of surface lubrication have chiefly focused on the theories and experiments only for the human rotational joints in the field of chemistry without considering any random changes [8,9,10,11,12,13,14,15,16,17,18,19]. There has been no previous stochastic research about the frictional forces, and the coefficients of friction based on experimental methods and analytical-numerical hydrodynamics for lamellar and laminar lubricant flows between two non-rotational biological surfaces [12,19]. Random changes in the gap height between two co-operating surfaces directly or indirectly influence the lubrication parameters, such as pressure, load-carrying capacity, and friction forces. The direct influence is demonstrated by an integral formula determining the lubrication parameters in the area of the lubricating layer. The indirect influence of random changes in the gap height on the value of lubrication parameters occurs through random changes in the dynamic viscosity of the lubricant biofluid, which randomly affects the changes in the friction forces and pressure [20]. In order to explain this phenomenon, it is necessary to consider the fact that with a random increase (decrease) in the gap height, the flow velocity of the fluid in the gap decreases (increases). This leads to a decrease (increase) in the shear rate. Owing to the non-Newtonian properties of the lubricant biofluid, dynamic viscosity would eventually increase (decrease). Moreover, the indirect influence of random changes in the region of an articular bio-surface coated with PL on the value of lubrication parameters occurs through random changes in the dynamic viscosity of the biofluid, which abruptly affects the changes in friction forces and pressure. In order to explain this phenomenon, it must be noted that the apparent viscosity of the biofluid η_T in the gap is directly proportional to the surface coating of PL [20,21,22,23,24]. The concentration of PL in the superficial layer of the surface flowed around, and the degree of wettability of these surfaces have a significant influence on the change in the viscosity of the physiological fluid accumulated in the boundary areas, values of the frictional forces, load-carrying capacity of the joint, and coefficient of friction [13]. Viscosity changes not only in the direction of the length and width of the flow but also in the direction of the nanometer thickness of the layer [20].

The objectives of the research presented in this work were as follows:

1.To study the direct and indirect influences of the randomly changing gap height between two co-operating, movable, non-rotational curvilinear, orthogonal bio-surfaces on the values of viscosity and flow velocity of the bio-lubricant, i.e., the bio-fluid, as well as on the values of hydrodynamic pressure and friction forces.

2.Using experimental measurements, define a new model of random variations of lubricant dynamic viscosity based on the features of the superficial layer of co-operating non-rotational bio-surfaces.

3.To present an estimation of the expected values of the main bio-tribology parameters based on the performed measurements, obtained probability density functions of the gap height, and semi-analytical and numerical solutions of stochastic hydrodynamic equations.

2. Materials and Tools for Measurements and Semi-Analytical Solutions: Random Gap Height, Density Function, and Standard Deviation

In the jump joint, the collar bone, or the blade bone, the lubrication regions are found between two movable, non-rotational but unparallel or seldom parallel co-operating surfaces (Figures 1a, 1b, and 2c). A new phenomenon of hydrodynamic sweat lubrication in a random form is observed in movable, non-rotational, curvilinear, non-parallel, or exceptionally parallel surfaces, e.g., in the external human skin surface and the tightly fitting dress surface during gymnastic training (Figures 1c, 1d, 1e, 2d, and 2e). Half-rotational, curvilinear, co-operating living surfaces are found in the human wrist joint and knee joint (Figures 1a, 1b, and 2b. Rotational co-operating biological bone surfaces are found in the hip joint and the elbow joint (Figures 1a, 1b, and 2a).

Figure 1 Biological human friction nods lubrication: a) Upper human limb: 1-collar bone, 2-blade bone, 3-radial bone, 4-wrist joint; b) Lower human limb: 5-hip bone, 6-sacral bone, 7-hip joint, 8-femoral bone, 9-knee cap, 10-knee joint, 11-calf bone, 12-tibia, 13-syndesmosis, 14-jump joint, 15-phalangeal joint; c), d), and e) human skin-tight dress (Based on author’s studies).

Figure 2 Random gap height ε_T between two co-operating biological surfaces: a) rotational; b) half-rotational; c) movable non-rotational; d), e) movable non-rotational (skin-sweat-sport dress); 1-phospholipid bilayer, 2-hydrated phospholipid, 3-biological liquid, 4-subchondral bone, 5-dress, 6-sweat, 7-opening of sweat duct, 8-human skin, 9-muscle, 10-gland, 11-hair follicle, 12-sweat duct, 13-sweat pressure, 14-sweat out, 15-heat flux; W-load, R-repulsion force, P_p-force of hydrodynamic pressure, ω-angular velocity (Based on author’s studies).

The variety in the curvature shapes of the analyzed non-rotational, co-operating bio-surfaces in the particular case of human joint surfaces dictates the description of the surface in a curvilinear orthogonal coordinate system: α₁, α₂, α₃. Hence, the α₁-width direction of the non-rotational surfaces or the circumferential direction is accepted in the case of a rotational surface, α₃-longitudinal direction, and α₂-gap height direction.

The joint gap is filled with physiological bio-fluid lubricant. A characteristic constant dimensional gap height value ε₀[m] and the dimensionless total gap height function ε_T1 based on variables α₁ and α₃ were accepted in this experiment. This function takes on the following equation [3]:

\[ \varepsilon_{T 1}=\frac{\varepsilon_{T}}{\varepsilon_{0}}=\varepsilon_{T 1 s}\left(\alpha_{1}, \alpha_{3}\right)\left[1+\delta_{1}\left(\alpha_{1}, \alpha_{3}\right)\right] \tag{1} \]

The dimensional symbol ε_T [m] denotes the total gap height. Symbol ε_T1s represents the dimensionless gap height limited by the nominally smooth biological movable and non-rotational surfaces. Random changes in the gap height are results of random hyper-elastic deformations of the articular superficial layer as well as random protrusions of the superficial layer or phospholipid bilayer surface roughness. The dimensionless random variable of corrections for the gap height is denoted by δ₁. The expected value for the random variable of corrections δ₁ is defined by equation (2a), and the expected function for the entire gap height is explained using equation (2b) [25,26]:

\[ E X\left(\delta_{1}\right) \equiv \int_{-\infty}^{+\infty}\left(\delta_{1}\right) \times f\left(\delta_{1}\right) d \delta_{1}, E X\left(\delta_{1}\right)^{2} \equiv \int_{-\infty}^{+\infty}\left(\delta_{1}\right)^{2} \times f\left(\delta_{1}\right) d \delta_{1} \tag{2a} \]

\[ E X\left(\varepsilon_{T 1}\right)=E X\left[\varepsilon_{T 1 s}\left(1+\delta_{1}\right)\right]=\varepsilon_{T 1 s}\left[1+E X\left(\delta_{1}\right)\right] \tag{2b} \]

The symbol EX denotes the operator of the expected function. Probability density function f assigns probability values as a function of the random variable of correction δ₁.

The ordinates of the density function f denote the probabilities established for the random variable of corrections δ₁ of the joint gap height. These values were determined experimentally, considering the articular superficial layer roughness. Standard deviation σ, from a random variable of corrections, is determined using the following equation [25,26]:

\[ \sigma \equiv \sqrt{E X\left(\delta_{1}\right)^{2}-E X^{2}\left(\delta_{1}\right)} \tag{2c} \]

Gap height ε_T, phospholipid-coated surfaces, apparent dynamic viscosity η_T, hydrodynamic pressure, temperature, and several other values are subject to random change corrections. The surface structure of the tested samples was irregular owing to the occurrence of random roughness (from 10 nm to 50 nm) or disease [20]. Based on the comparisons made between the random changes in the rough superficial layer surface structure measured in Cwanek’s research [1] and Dowson’s measurements [2], the probability density function was found to be unsymmetrical. This suggests that there is a difference in the probabilities of random increases and decreases in the gap height.

The measurements of the random height changes in the surface roughness on the co-operating, non-rotational bio-surfaces were conducted using a sample (10 mm × 10 mm) of diseased superficial layer taken from a human blade bone, as shown in Figure 3a and Figure 3b., and a sample (1 mm × 1 mm) of artificial friction nod cf. Figure 3c and Figure 3d. The experiments were conducted using a laser micro-sensor installed in a Rank Taylor Hobson-Talyscan 150 Apparatus. The results were compiled using Talymap Expert and Microsoft Excel computer software. The measured gap height limited by the rough surfaces of the articular superficial layer of the presented sample exhibited a variation in the range of 0.005 mm and 0.35 mm [1].

Figure 3 Friction nod surface structure and visualization of random changes δ₁ in total gap height ε_T caused by hyper-elastic deformations, unsteady load, and diseases, for sample: a), b), 10 mm × 10 mm, in the articular friction nod; and for sample: c), d), 1.0 mm × 1.0 mm, in the artificial friction nod; 1-surface limiting the gap height, 2-smooth surface limiting the gap without random changes, 3-random variable surface (Based on author’s studies).

Figure 3 illustrates the measured positive and negative values of random variable gap height variations δ₁, which represent the differences between the probabilistic, variable surface coordinate 3 and its primary localization on the theoretical smooth surface 2, limiting the gap height without any random changes.

3. Semi-Analytical Methods

3.1 Basic Lubrication Equations of Two Arbitrary Surfaces in Curvilinear Coordinates

The lubrication problem of the human movable, non-rotational friction nods is represented using the conservation of momentum, continuity, Maxwell equations, conservation of energy, and Young-Kelvin-Laplace’s equation for biological non-Newtonian liquid, particularly for a synovial or sweat liquid lubrication flow. A non-rotational, unsteady, non-isothermal, incompressible flow of viscoelastic biological liquid in an electromagnetic field has been assumed. The above-mentioned equations are in the following forms [21]:

The equation of motion:

\[ \operatorname{Div} \boldsymbol{S}+\mu_{m}(\boldsymbol{N} \nabla) \boldsymbol{H}+\boldsymbol{J} \times \boldsymbol{B}+\rho_{e} \boldsymbol{E}=\rho^{d \boldsymbol{v}} / d t \tag{3a} \]

The continuity equation:

\[ \operatorname{div}(\boldsymbol{v})=0 \tag{3b} \]

The reduced Maxwell equation:

\[ \nabla^{2} \boldsymbol{H} \equiv \mu_{m} \mu_{e}^{\partial^{2}} \boldsymbol{H} / \partial_{\partial t^{2}}, \boldsymbol{J}=\sigma_{e} \boldsymbol{E} \tag{3c} \]

The conservation of energy equation:

\[ \operatorname{div}(\kappa \operatorname{grad} T)+\Phi_{F}-\mu_{m} T \Xi(\boldsymbol{v} \nabla) \boldsymbol{H}=0, \Phi_{F}=\operatorname{div}(\boldsymbol{v} \boldsymbol{S})-\boldsymbol{v} \operatorname{Div} \boldsymbol{S} \tag{3d} \]

The Young-Kelvin-Laplace equation:

\[ \gamma=\gamma_{\max }+2 s R_{g} T \ln \left(\sqrt{\frac{K_{a}}{K_{b}}}+1\right)-s_{c} R_{g} T \ln \left[\left(\frac{K_{a}}{a_{H}^{+}}+1\right)\left(\frac{a_{H}^{+}}{K_{b}}+1\right)\right] \tag{3e} \]

where σ [Pa] is the stress tensor in the biological fluid, E [V/m] denotes the electric intensity vector, J [A/m²] represents the electric current density, B [T = kg/s²A] indicates the magnetic induction vector in bio-fluid, v [m/s] is the biological fluid velocity, H [A/m] represents the magnetic intensity vector with components (H₁, H₂, H₃), N [A/m] denotes the magnetization vector with components (N₁, N₂, N₃), σ_e [S/m] is the electrical conductivity of phospholipids bilayer, μ_m [mkgs^–2A^–2] signifies the magnetic permeability coefficient of biological fluid, μ_e [s⁴A² m^–3 kg^–1] is the electric permeability coefficient of the biological fluid, κ[W/mK] denotes the thermal conductivity coefficient for the biological fluid, Φ_F [W/m³] indicates the dissipation of energy, Ξ [A/mK] denotes the first derivative of the magnetization vector with respect to temperature, T [K] signifies the temperature, ρ_e [C/m³ = As/m³] is the electric space charge in the biological fluid, R_g indicates the gas constant (8.3144598 J/K·mol), s_c = (N_A∙A)^–1 [mol/m²] is the concentration of phospholipid particles, γ[J/m² = N/m] denotes the interfacial energy, γ_max is the maximum interfacial energy of the lipid membrane such that 0.1 mN/m < γ_max < 4 mN/m, K_a [J] is the acid equilibrium constant (denotes the amount of energy needed to stretch the bilayer), K_b [J] represents the base equilibrium constant (denotes the amount of energy needed to bend or flex the bilayer), α_H [J] indicates the energy activity of protons, A [m²] denotes the boundary regions between the areas of different concentrations of phospholipid molecules, N_A = 6.024∙10²³ [1/mol] represents the Avogadro number, and ρ[kg/m³] signifies the biological fluid density. Due to the presence of phospholipid bi-layers on the cartilage or superficial layer surface and the presence of liposomes, micelles, macromolecules, and lamellar aggregates in the biological fluid, this liquid has non-Newtonian as well as pseudo-plastic properties.

For the synovial fluid, the relation between stress tensor σ and displacement velocity tensor 2T_d = Θ, i.e., the constitutive equations, are accepted in the following form [21]:

\[ \boldsymbol{S}=-p \boldsymbol{\delta}+\eta_{T} \boldsymbol{\Theta} \tag{3f} \]

where unit tensor δ and strain tensor Θ [s^–1] possess the following components: δ_ij, and Θ_ij. The δ_ij denotes Kronecker Delta and p [Pa] is the pressure. For a non-Newtonian synovial fluid with Rivlin-Reiner type, the constitutive dependencies for apparent viscosity η_T [Pas] take the following form [21]:

\[ \eta_{T}=2^{n-1} M_{\text {con }} \cdot\left|\frac{1}{2} I_{1}^{2}(\Theta)-I_{2}(\Theta)\right|^{\frac{n-1}{2}}, I_{1}(\theta)=\Theta_{k k}, I_{2}(\Theta)=\frac{1}{2} e_{i j k} e_{i m n} \Theta_{j m} \Theta_{k n} \tag{3g} \]

where I₁(Θ)[s^–1], I₂(Θ)[s^–2] are the invariants of the shear rate components Θ_ij [s^–1], n is the dimensionless flow index dependent on the phospholipid concentration in the biological fluid (0.8 ≤ n ≤ 1.2), M_con (n, pH, T, We) [Pasⁿ] denotes the fluid consistency coefficient, e_ijk indicates tensor Levi-Civita, 50⁰ ≤ We ≤ 80⁰ represents the superficial layer wettability, and 2 ≤ pH ≤ 12 denotes the power hydrogen ion concentration. Geometrical non-linear relations between shear rate Θ_ij [s^–1] and biological fluid velocity component v_i [m/s] are as follows [21]:

\[ \Theta_{i j}=\frac{1}{2}\left(v_{i \mid j}+v_{j \mid i}\right), v_{i \mid j} \equiv \frac{1}{h_{i}}\left(\frac{\partial v_{i}}{\partial \alpha_{j}}-\frac{v_{j}}{h_{i}} \frac{\partial h_{j}}{\partial \alpha_{i}}+\delta_{i j} \sum_{k=1}^{3} \frac{v_{k}}{h_{k}} \frac{\partial h_{i}}{\partial \alpha_{k}}\right) \tag{3h} \]

where h_i denotes the Lamé coefficient. The lubrication equation in curvilinear coordinates for the conservation of momentum, continuity equation, energy conservation equation, and Young-Kelvin-Laplace equation (3a, b, c, d, e) are expanded and written in the following curvilinear directions: α₁, α₂, α₃. The expected functions (2a, b) of hydrodynamic pressure EX [p(α₁,α₃)], temperature EX [T(α₁,α₂,α₃)], the velocity of synovial fluid EX [v_i(α₁,α₂,α₃)], the viscosity of lubricating fluid EX [η_T(α₁,α₂,α₃)], and gap height EX [ε_T(α₁,α₃)] are taken into consideration. The incompressibility of biological fluid enables mitigating the influence of changes in the density of bio-fluid on the equation for biological fluid flow continuity. The known dependencies between interfacial energy γ[N/m] and power of hydrogen ion concentration, pH, and wettability We [14,15,16,18] are used. In the basic equations (3a, b, c, d), the influence of the electrostatic field on the viscosity of bio-fluids studied by the authors [21] is considered, while the magnetic field is neglected. Further, in equations (3a, b, c, d, e, f), the classic boundary simplification of the hydrodynamic equations is applied in the boundary layer by omitting the terms of the order of relative radial clearance ψ with a value of 10^–4. ψ is defined as the ratio of the thickness of the thin bio-fluid layer ε_T [m] to the length of the lubricated region of the movable, non-rotational surface, or curvature radius of the rotational surface flowed around. For two arbitrary co-operating biological surfaces, the curvilinear, orthogonal system of coordinates α₁, α₂, α₃ is applied to the respective Lamé coefficients h₁, h₂, h₃. The aforesaid boundary simplifications of the arbitrary, curvilinear, non-rotational, non-parallel thin-layer surfaces follows h₂ = 1 and h₁ = h₁(α₁,α₃), h₃ = h₃(α₁,α₃). This case is valid for the jump joint, collar bone, and blade bone coated with a thin liquid layer. Moreover, this case is applicable for the thin sweat layer between the non-rotational, movable, external skin surface and leggings or a tight training dress surface. For the thin liquid layer restricted by the two rotational surfaces in α₁ direction and the non-monotone-generating line in α₃ direction, the Lamé coefficients follow h₂ = 1, h₁ = h₁(α₃), h₃ = h₃(α₃), or h₃ = 1 for the monotone generating line. This case is valid for the human elbow joint and hip joint lubrication. After performing the above transformations with boundary simplifications and applying the probabilistic changes (2a, b, c), taking PL into account [20,21,22], the system of equations (3a, b, c, d) for the hydrodynamic lubrication of two arbitrary, non-rotational, biological co-operating surfaces tend to take the following form [26,27,28]:

\[ 0=-\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}+\frac{\partial}{\partial \alpha_{2}}\left[E X \eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \frac{\partial E X\left(v_{1}\right)}{\partial \alpha_{2}}\right]+\rho_{e} E_{1} \tag{4} \]

\[ 0=\frac{\partial E X(p)}{\partial \alpha_{2}} \tag{5} \]

\[ -\frac{\rho E X\left(v_{1}\right)^{2}}{h_{1} h_{3}} \frac{\partial h_{1}}{\partial \alpha_{3}}=-\frac{1}{h_{3}} \frac{\partial E X(p)}{\partial \alpha_{3}}+\frac{\partial}{\partial \alpha_{2}}\left[E X \eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \frac{\partial E X\left(v_{3}\right)}{\partial \alpha_{2}}\right]+\rho_{e} E_{3} \tag{6} \]

\[ \frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{1}}\left[h_{3} E X\left(v_{1}\right)\right]+\frac{\partial E X\left(v_{2}\right)}{\partial \alpha_{2}}+\frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left[h_{1} E X\left(v_{3}\right)\right]=0 \tag{7} \]

\[ \frac{\partial}{\partial \alpha_{2}}\left\{\kappa \frac{\partial E X\left[T\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}{\partial \alpha_{2}}\right\}+E X \eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\left\{\left[\frac{\partial E X\left(v_{1}\right)}{\partial \alpha_{2}}\right]^{2}+\left[\frac{\partial E X\left(v_{3}\right)}{\partial \alpha_{2}}\right]^{2}\right\}=J^{2} / \sigma_{e} \tag{8} \]

Combining equation (3e) with equation (3f-h) using the transformations among apparent dynamic viscosity, molecular velocity, and interfacial energy, it becomes evident that the expected function of viscosity η_T [Pas] [21,26] is as follows:

\[ E X\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=E X\left[\eta_{T}\left(n, p_{H}, W_{e}, T, \gamma, E\right)\right] \equiv \]

\[ =\frac{\gamma_{\max }\left(p_{H}, W_{e}\right)+k E X\left(A^{-1}\right) \cdot E X(T) \ln L}{\delta_{v} \cdot E X\left(v_{0}\right)} \]

\[ \left[1+\delta_{E}\left(p_{H}, E\right) \cdot E^{2}\right]\left(\sqrt{\left[\frac{\partial E X\left(v_{11}\right)}{\partial \alpha_{21}}\right]^{2}+\left[\frac{\partial E X\left(v_{31}\right)}{\partial \alpha_{21}}\right]}\right)^{n-1} \tag{9a} \]

$$0<L\equiv \frac{{\left(\sqrt{L_k}+1\right)}^2}{\left(L_a+1\right)\left(L_b+1\right)}<1,L_a\equiv \frac{K_a}{a_H^{+}},L_b\equiv \frac{a_H^{+}}{K_b},L_k\equiv L_a{L}_b$$

\[ \left(L_{a}+1\right)\left(L_{b}+1\right)>\left(\sqrt{L_{k}}+1\right)^{2} \tag{9b} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, and 0 ≤ α₂ ≤ EX(ε_T). Symbol 2a [m] and 2b [m] denote the width and length, respectively, of the friction region between the two co-operating surfaces. k = 1.38054∙10^–23 [J/K] represents the Boltzmann constant. It follows the formula (9a), according to which, apparent viscosity η_T of the biological liquid in the gap is directly proportional to the boundary surfaces A [m²] between the areas of different phospholipid (PL) concentration having different superficial layer susceptibility of Symbol δ_v occurring in formula (9a) represents the dimensionless random coefficient (0.2 < δ_v < 0.6). Suppose index δ_v increases from 0.2 to 0.6, the dynamic viscosity of the bio-fluid undergoes a significant reduction. Coefficient δ_v describes the concentration c_c of collagen fibers in the bio-fluid. For δ_v = 0.2, the value of c_c is 100,000 mol/mm³, while for δ_v = 0.6 concentration, the value of c_c is 100 mol/mm³. The dimension of parameters L_a and L_b was determined by the formula (9b). Since 0 < L < 1, lnL is a negative number, while it is always the case that γ_max >-kA^–1T lnL, where –100 < lnL < –10.

The symbol 0.03 m/s <v₀ < 0.04 m/s indicates the characteristic dimensional value of linear velocity for the bio-surface being flowed around. The derivatives in the dimensionless gap height direction α₂₁ of the dimensionless functions v₁₁ and v₃₁ of velocity vector components v₁[m/s] and v₃[m/s] in directions α₁ and α₃ in formula (9a) are the result of dimensionless transformations, and calculations of constitutive dependences of apparent viscosity (3g), where v₁ = v₀⋅ v₁₁, v₃ = v₀ ⋅v₃₁. Newtonian fluid is used if the dimensionless flow index n (0.8 < n < 1.2) has the value of ‘1’. The following points have been assumed: δ_E [m²/V²] is the experimental coefficient of the influence of electrical intensity E and concentration of hydrogen ions pH in the bio-fluid. The value of coefficient δ_E [m²/V²] for the bio-fluid has not yet been accurately measured. Though, it has been observed that for E = 10 V/m and pH = 8, the value of δ_E = 0.0003 m²/V² [29]. It thus follows that the dimensionless increase in the dynamic viscosity of the bio-fluid resulting from the influence of the electrostatic field in the phospholipid layer is 1+δ_EE² = 1.03, which is only a 3% increase. Based on the obtained estimations, it can be inferred that the direct influence of the electrostatic field in the tissue boundary layer is negligible.

Formula (9) was derived from known dependences of interfacial energy γ(pH, We), describing the relationship between the power of hydrogen ion concentration pH and wettability We. Interfacial energy was analytically transformed into apparent dynamic viscosity of the lubricating fluid η_T. Such dependencies are explained graphically in [21] for two types of phospholipids, PC and PS, with dimensionless values for acid pKa and base pKb equilibrium constants. An increase in pKa enhances the viscosity increase in the interval of 2 < pH < 4 and enhances the viscosity decrease in the interval of 4 < pH < 10. Bio-fluid dynamic viscosity for PC and PS lipids increases with the increasing power of hydrogen ion concentration pH to a certain iso-electric point IP (γ = 3.5 mJ/m²), with the established values of We, δ_v, T, and v₀. A further increase in pH leads to a decrease in dynamic viscosity. The dynamic viscosity of the bio-fluid decreases with decreasing wettability We at established values δ_v, T, and v₀. A drop in wettability from 70^° to 50^° indicates a transition from hydrophobic to hydrophilic properties in the bio-surfaces flowed around by the bio-fluid.

Similarly, for formula (1) and formula (2a), by specifying the random gap height changes, the expected functions of pressure, temperature, and lubricating fluid velocity components in the system of equations (4)-(8) take the following forms [26]:

\[ E X(p)=p \cdot\left[1+E X\left(\delta_{p}\right)\right], E X(T)=T \cdot\left[1+E X\left(\delta_{T}\right)\right] \]

\[EXv_{i}=v_{i}\cdot\left[\begin{array}{ll}\left.1+E X\left(\delta_{vi}\right)\right], \quad\text {for}\quad i=1,2,3 \end{array}\right.\tag{10} \]

where the symbols δ_p, δ_T, and δ_vi represent unknown random variables of pressure, temperature, and biological fluid velocity component corrections. If the random variable of gap height corrections δ₁ = 0, then δ_p = δ_T = δ_vi = 0. In this case EX(δ_p) = EX(δ_T) = EX(δ_vi) = 0; therefore, under the assumption (10), EX(p) = p, EX(T) = T, and EX(v_i) = v_i for i = 1, 2, 3. Thus, in this specific case, the system of equations (4)-(8) loses the stochastic properties.

The system of partial differential equations (4)-(8) determines the following expected functions of random variable unknowns: three bio-fluid velocity vector components EX [v_i(α₁,α₂,α₃)] [m/s] for i = 1, 2, 3; hydrodynamic pressure EX [p(α₁,α₃)] [Pa]; and temperature EX [T(α₁,α₂,α₃)] [K]. The term on the left side of equation (6) indicates the centrifugal forces occurring during the cooperation of two bio-surfaces. These forces occur only when Lamé coefficient h₁ is a function of α₃. This is applicable either to the non-rotational curvilinear bio-surface or to rotational surfaces with a non-monotone generating line, i.e., spherical, conical, parabolic, or elliptical shapes, but not a cylindrical shape where coefficient h₁ is a constant value.

3.2 Integrals for Random Hydrodynamic Lubrication of Non-Rotational Surfaces

Integration of the system of equations (4)-(8) describing the lubrication of two non-rotational bio-surfaces with the participation of phospholipids (PL) was carried out in curvilinear coordinates (α₁,α₂,α₃). These phospholipids are separated by a thin layer of biological fluid, with viscosity variations in the gap height direction. This lubricated bio-surface performs a non-rotational motion with linear velocity U₁(α₁,α₃) [m/s] in α₁ direction and U₃(α₁,α₃) [m/s] in α₃ direction, while the second bio-surface is stationary and limits a gap filled with a layer of bio-fluid with randomly varying height ε_T. The following boundary conditions are applied to the functional components of the expected randomly varying viscous bio-fluid velocities EX(v₁), EX(v₂), and EX(v₃) in the directions α₁, α₂, and α₃, respectively [25,28]:

\[ \operatorname{EX}\left(v_{1}\right)=U_{1} \text { for } \alpha_{2}=0, E X\left(v_{1}\right)=0 \text { for } \alpha_{2}=E X\left(\varepsilon_{T}\right) \tag{11} \]

\[ E X\left(v_{2}\right)=0 \text { for } \alpha_{2}=0, \text { and } E X\left(v_{2}\right)=0 \text { for } \alpha_{2}=E X\left(\varepsilon_{T}\right) \tag{12} \]

\[ E X\left(v_{3}\right)=U_{3} \text { for } \alpha_{2}=0, E X\left(v_{3}\right)=0 \text { for } \alpha_{2}=E X\left(\varepsilon_{T}\right) \tag{13} \]

By applying condition (11) to the general solution of equation (4), the following form of function of the expected randomly variable component of the bio-fluid velocity vector is obtained in direction α₁ of the non-rotational movable bio-surface [23,26]:

\[ \operatorname{EX}\left[v_{1}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right) A_{\eta}+\left(1-A_{\mathrm{s}}\right) U_{1} \tag{14} \]

where the subordinate functions A_s [1] and A_η [m⁴/Ns] are as follows:

\[ A_{s}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv \frac{\int_{0}^{\alpha_{2}} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}}{\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}}, A_{\eta}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv \]

\[ \int_{0}^{\alpha_{2}} \frac{\alpha_{2}}{E X\left(\eta_{T}\right)} d \alpha_{2}-\frac{\left(\int_{0}^{\alpha_{2}} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}\right) \cdot\left(\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{\alpha_{2}}{E X\left(\eta_{T}\right)} d \alpha_{2}\right)}{\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}} \tag{15} \]

Here, 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, 0 ≤ α₂ ≤ EX(ε_T), EX(ε_T) = EX [ε_T(α ₁,α₃)], h₁α₁[m], h₃α₃[m], α₂[m], a [m], and b [m], U₁ = U₁(α₁,α₃) [m/s]. Besides, A_s [α₁,α₂ =EX(ε_T),α₃] = 1, and A_η[α₁,α₂ =EX(ε_T),α₃] = 0. Substituting solution (14) into equation (6) and then applying condition (13), the function of the expected randomly variable component of the bio-fluid velocity vector is obtained in the longitudinal direction α₃ of the non-rotational bio-surface:

\[ E X\left[v_{3}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=\left(\frac{1}{h_{3}} \frac{\partial E X(p)}{\partial \alpha_{3}}-M_{3}\right) A_{\eta}+\left(1-A_{\mathrm{s}}\right) U_{3}-\frac{\rho}{h_{1} h_{3}} \frac{\partial h_{1}}{\partial \alpha_{3}} A_{p} \tag{16} \]

where U₃ = U₃(α₁, α₃) [m/s]. Moreover, M_i = ρ_eE_i [N/m³] for i = 1, 3 are the electrical terms.

The last term on the right-hand side of equation (16) represents the influence of the suction or inertial force of the non-rotational bio-surface on the distribution of the bio-fluid velocity in the gap. The characteristic function of the centrifugal effect A_p [m⁶/Ns³] was assumed to be as [22]:

\[ A_{p}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right)^{2} A_{\rho 1}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)+ \]

\[ -2 U_{1}\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right) A_{\rho 2}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)+U_{1}^{2} A_{\rho 3}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \tag{17a} \]

The auxiliary functions for non-rotational surfaces and the viscosity variations in gap height directions are represented as follows [21,22]:

\[ A_{\rho 1}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv \int_{0}^{\alpha_{2}}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}} A_{\eta}^{2} d \alpha_{2}\right) d \alpha_{2}-A_{s} \cdot \int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}} A_{\eta}^{2} d \alpha_{2}\right) d \alpha_{2} \]

\[ A_{\rho 2}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv \]

\[ \int_{0}^{\alpha_{2}}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}}\left(1-A_{s}\right) A_{\eta} d \alpha_{2}\right) d \alpha_{2}-A_{s} \cdot \int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}}\left(1-A_{s}\right) A_{\eta} d \alpha_{2}\right) d \alpha_{2} \]

\[ A_{\rho 3}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right) \equiv \]

\[ \int_{0}^{\alpha_{2}}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}}\left(1-A_{s}\right)^{2} d \alpha_{2}\right) d \alpha_{2}-A_{s} \cdot \int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}}\left(1-A_{s}\right)^{2} d \alpha_{2}\right) d \alpha_{2} \tag{17b} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, 0 ≤ α₂ ≤ EX(ε_T), EX(ε_T) = EX [ε_T(α ₁, α₃)], h₁α₁[m], h₃α₃[m], α₂[m], a [m], b [m], U₁ = U₁(α₁,α₃)[m/s], A_ρ1[m¹²/N³s³], A_ρ2[m⁸/N²s²], and A_ρ3[m⁴/Ns].

The influence of the centrifugal forces disappears when the Lamé coefficient h₁ is independent of α₃. Such a situation occurs if there are two co-operating plates or rotational bio-surfaces, like cylindrical co-operating surfaces where the generating line is a straight line parallel to the axis of rotation. Integration of the continuity equation (7) with the boundary condition (12), where v₂ = 0 for α₂ = 0, leads to the following form of function of the expected random variable component of bio-fluid velocity in the direction α₂ of gap height [20,26]:

\[ E X\left[v_{2}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=-\int_{0}^{\alpha_{2}} \frac{1}{h_{1} h_{3}} \frac{\partial\left[h_{3} E X\left(v_{1}\right)\right]}{\partial \alpha_{1}} d \alpha_{2}-\int_{0}^{\alpha_{2}} \frac{1}{h_{1} h_{3}} \frac{\partial\left[h_{1} E X\left(v_{3}\right)\right]}{\partial \alpha_{3}} d \alpha_{2} \tag{18} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], and b [m].

Further, the functions (14)-(16) are substituted into equation (18), and then boundary condition (12) is applied in the form of EX(v₂) = 0 for α₂ = EX(ε_T), i.e., to the component of the bio-fluid velocity vector in the direction of gap height α₂. The new forms (19a, 19b, 19c) of the hydrodynamic pressure equations are derived after intense calculations.

For arbitrary, cooperating, biological, non-rotational surfaces in α₁ and α₃, directions, where h₂ = 1, h₁ = h₁(α₁,α₃), and h₃ = h₃(α₁,α₃), for a movable surface with linear velocity U₁ in α₁ direction and U₃ in α₃ direction, the stochastically modified Reynolds type equation is derived. The bio-surfaces are lubricated with a thin liquid layer with dynamic viscosity η_T(a₁,α₂,α₃) varying in three directions. This equation determining the expected function EX [p(α₁,α₃)] of randomly variable hydrodynamic pressure has the following form:

\[ \frac{1}{h_{1}} \frac{\partial}{\partial \alpha_{1}}\left[\left(\frac{\partial E X(p)}{h_{1} \partial \alpha_{1}}-M_{1}\right) \cdot\left(\int_{0}^{E X\left(\varepsilon_{r}\right)} A_{\eta} d \alpha_{2}\right)\right]+\frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left[h_{1}\left(\frac{\partial E X(p)}{h_{3} \partial \alpha_{3}}-M_{3}\right) \cdot\left(\int_{0}^{E x\left(\varepsilon_{\tau}\right)} A_{n} d \alpha_{2}\right)\right] \]

\[ -\rho \frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left(\frac{\partial h_{1}}{h_{3} \partial \alpha_{3}} \int_{0}^{E x\left(\varepsilon_{r}\right)} A_{p} d \alpha_{2}\right)=\frac{U_{1}}{h_{1}} \frac{\partial}{\partial \alpha_{1}}\left[\left(\int_{0}^{E X\left(\varepsilon_{r}\right)} A_{s} d \alpha_{2}-E X\left(\varepsilon_{T}\right)\right)\right],\tag{19a} \]

where, 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], b [m], U₁ = U₁(α₁, α₃), and U₃ = U₃(α₁, α₃) [m/s]. The terms multiplied by the linear velocity components U₁ and U₃ represent the influence of surface motion on the hydrodynamic pressure. The integrals mentioned in equation (19a) are enclosed in Appendix 1. The mathematical derivations presented in the two appendices are crucial as they illustrate the influence of bio-fluid viscosity variations across the very thin gap height on the lubrication parameters. The aforementioned influences have never been referred to in any of the previous works on classical lubrication. Hence, such derivations have not been considered until now.

For rotational, biological, movable surfaces with angular velocity ω in α₁ direction and non-monotone generating line in α₃ direction, where h₂ = 1, h₁ = h₁(α₃), and h₃ = h₃(α₃), for the motionless surface in α₃ direction, where, U₁(α₃) = ω⋅h₁, U₃ = 0, and for the surfaces lubricated with a thin liquid layer with dynamic viscosity varying in three directions, η_T(α₁,α₂,α₃), the following stochastically modified Reynolds equation, which determine the expected function EX [p(α₁,α₃)] of randomly variable hydrodynamic pressure, is valid:

\[\frac{1}{h_{1}} \frac{\partial}{\partial \alpha_{1}}\left[\left(\frac{\partial E X(p)}{h_{1} \partial \alpha_{1}}-M_{1}\right) \cdot\left(\int_{0}^{E x\left(\varepsilon_{\tau}\right)} A_{n} d \alpha_{2}\right)\right]+\frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left[h_{1}\left(\frac{\partial \operatorname{EX}(p)}{h_{3} \partial \alpha_{3}}-M_{3}\right) \cdot\left(\int_{0}^{E X\left(\varepsilon_{T}\right)} A_{n} d \alpha_{2}\right)\right] \]

\[ -\rho \frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left(\frac{\partial h_{1}}{h_{3} \partial \alpha_{3}} \int_{0}^{E x\left(\varepsilon_{T}\right)} A_{p} d \alpha_{2}\right)=\frac{U_{1}}{h_{1}} \frac{\partial}{\partial \alpha_{1}}\left[\left(\int_{0}^{E X\left(\varepsilon_{r}\right)} A_{s} d \alpha_{2}-E X\left(\varepsilon_{T}\right)\right)\right],\tag{19b} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], b [m], and U₁ = U₁(α₃).

For rotational, cooperating, biological surfaces with angular velocity ω in α₁ direction and monotone motionless generating line in α₃ direction, where h₂ = 1 and h₁ = h₁(α₃), U₁ = ω⋅h₁, U₃ = 0, and for the surfaces lubricated with a thin liquid layer with dynamic viscosity varying in three directions, η_T(α₁,α ₂,α₃), the stochastically modified Reynolds equation (19b) is valid with h₃ = 1.

For rotational, co-operating, biological surfaces with angular velocity ω in α₁ direction and motionless in α₃ non-monotone direction, where h₂ = 1, h₁ = h₁(α₃), and h₃ = h₃(α₃), U₁ = U₁(α₃) = ω⋅h₁, U₃ = 0, and for the surfaces lubricated with a thin liquid layer with dynamic viscosity varying only in two directions and being constant in the gap height direction, η_T(α₁,α₃), the following stochastically modified Reynolds equation, which determine the expected function EX [p(α₁,α₃)] of random variable hydrodynamic pressure, is obtained:

\[ \frac{1}{h_{1}} \frac{\partial}{\partial \alpha_{1}}\left[\frac{E X\left(\varepsilon_{T}^{3}\right)}{E X\left(\eta_{T}\right)}\left(\frac{\partial E X(p)}{h_{1} \partial \alpha_{1}}-M_{1}\right)\right]+\frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left[\frac{E X\left(\varepsilon_{T}^{3}\right)}{E X\left(\eta_{T}\right)} h_{1}\left(\frac{\partial E X(p)}{h_{3} \partial \alpha_{3}}-M_{3}\right)\right] \]

\[ +\rho \frac{1}{h_{1} h_{3}} \frac{\partial}{\partial \alpha_{3}}\left(\frac{\partial h_{1}}{h_{3} \partial \alpha_{3}} \cdot 12 \cdot \int_{0}^{E X\left(\varepsilon_{T}\right)} A_{p} d \alpha_{2}\right)=6 \omega h_{1} \frac{\partial E X\left(\varepsilon_{T}\right)}{h_{1} \partial \alpha_{1}} \tag{19c} \]

where

\[ \text { 12. } \int_{0}^{E X\left(\varepsilon_{T}\right)} A_{p} d \alpha_{2}=-\frac{3}{10} \frac{E X\left(\varepsilon_{T}^{3}\right)}{E X\left(\eta_{T}\right)} \]

\[ \left[\frac{1}{28} \frac{E X\left(\varepsilon_{T}^{4}\right)}{E X\left(\eta_{T}^{2}\right)} \cdot\left(\frac{\partial E X(p)}{h_{1} \partial \alpha_{1}}-M_{1}\right)^{2}+\frac{\omega h_{1}}{3} \frac{E X\left(\varepsilon_{T}^{2}\right)}{E X\left(\eta_{T}\right)} \cdot\left(\frac{\partial E X(p)}{h_{1} \partial \alpha_{1}}-M_{1}\right)+\left(\omega h_{1}\right)^{2}\right] \tag{19d} \]

Here, 0 ≤ h₁α₁ < a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], and b [m].

The terms in equations (19a, 19b, 19c), when multiplied by the biological fluid density ρ, illustrate the effect of centrifugal forces on the hydrodynamic pressure.

The expected function of randomly varying pressure EX(p) is determined from equations (19a, 19b, 19c), assuming that the value of atmospheric pressure p_A at the edges of the lubrication area to be Ω(α₁,α₃):

\[ E X\left[p\left(\alpha_{1}, \alpha_{3}\right)\right]=p_{A} \text { for }\left(\alpha_{1}, \alpha_{3}\right) \in F r \Omega, \Omega \in\left(0 \leq h_{1} \alpha_{1} \leq {\rm a}\right) \times\left(-b \leq h_{3} \alpha_{3} \leq {\rm b} \right) \tag{20} \]

where Fr denotes the topological boundary set operator.

Now, the expected functions of the bio-fluid velocity vector components (14) and (16) are substituted into the energy equation (8). A constant value of bio-fluid thermal conductivity κ has been considered. For non-rotational, arbitrary, co-operating, biological surfaces in α₁ and α₃ directions, where h₂ = 1 and h₁ = h₁(α₁,α₃), and h₃ = h₃(α₁,α₃), movable in α₁ direction with linear velocity U₁ and in α₃ direction with linear velocity U₃, and for the surfaces lubricated with the liquid with dynamic viscosity varying in three directions:η_T(α₁,α₂,α₃), the following differential equations are obtained after the transformations, which enable the determination of the expected function of temperature EX [T(α ₁, α ₂,α₃)] in a random variable form [22,25,28]:

\[ \frac{\partial^{2} E X(T)}{\partial \alpha_{2}^{2}}+ \]

\[ +\frac{E X\left(\eta_{T}\right)}{\kappa}\left\{\left[\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right) \frac{\partial A_{\eta}}{\partial \alpha_{2}}-U_{1} \frac{\partial A_{S}}{\partial \alpha_{2}}\right]^{2}+\left[\left(\frac{1}{h_{3}} \frac{\partial E X(p)}{\partial \alpha_{3}}-M_{3}\right) \frac{\partial A_{\eta}}{\partial \alpha_{2}}-U_{3} \frac{\partial A_{S}}{\partial \alpha_{2}}\right]^{2}\right\}+ \]

\[ +\frac{E X\left(\eta_{T}\right)}{\kappa} \frac{\rho}{h_{1} h_{3}} \frac{\partial h_{1}}{\partial \alpha_{3}} \frac{\partial A_{p}}{\partial \alpha_{2}}\left[\frac{\rho}{h_{1} h_{3}} \frac{\partial h_{1}}{\partial \alpha_{3}} \frac{\partial A_{p}}{\partial \alpha_{2}}-2\left(\frac{1}{h_{3}} \frac{\partial E X(p)}{\partial \alpha_{3}}-M_{3}\right) \frac{\partial A_{\eta}}{\partial \alpha_{2}}\right]=\frac{M_{T}}{\kappa}, \tag{21} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, 0 ≤ α₂ ≤ EX(ε_T), h₁α₁[m], h₃α₃[m], α₂[m], a [m], b [m], U₃ = U₃(α₁, α₃), U₁ = U₁(α₁,α₃), and M_T/κ=J²/σ_ek[K/m²]. The derivatives of the terms A_s, A_η, and A_p are presented in Appendix 2. The terms in (21), when multiplied by the linear velocity components U₁ and U₃, represent the influence of surface motion on temperature.

In order to determine the expected function of randomly variable temperature EX [T(α₁,α₂,α₃)] from the second-order differential equation (21), two boundary conditions are required. The variations in the expected function of temperature below and above the characteristic temperature T₀ ultimately lead to a constant temperature value f_c on the first bio-surface (movable) and an unknown variable value of temperature changes f_p(α₁,α₃) on the second bio-surface (immovable). Hence, the two derived boundary conditions are as follows:

\[ E X\left[T\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=T_{0}+f_{c} \text { for } \alpha_{2}=0 \tag{22a} \]

\[ \operatorname{EX}\left[T\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]=T_{0}+f_{p}\left(\alpha_{1}, \alpha_{3}\right) \text { for } \alpha_{2}=E X\left(\varepsilon_{T}\right) \tag{22b} \]

In order to determine the unknown temperature function f_p(α₁,α₃) on the surface of the acetabulum, the condition of transportation of heat flux density q_c from the bonehead surface through the bio-liquid layer to the surface of acetabulum is used. This condition has the following form [21]:

\[ \kappa \frac{\partial E X(T)}{\partial \alpha_{2}}=-q_{c} \text { for } \alpha_{2}=0 \tag{22c} \]

3.3 Random Frictional Forces, Coefficient of Friction, and Load Carrying Capacity for Non-Rotational Surfaces in Curvilinear Form

The components of expected random functions of friction forces F_R1, F_R3 [N] are determined in curvilinear directions α₁ and α₃ occurring during the lubrication of two non-rotational bio-surfaces, with the participation of phospholipids (PL), separated by a thin layer of biological fluid with viscosity variations in the gap height directions, η_T(α₁,α₂,α₃₎.Mentioned friction forces are formulated in the curvilinear orthogonal coordinates (α₁,α₂,α₃) in the following forms [3,23,25,26]:

\[ E X\left(F_{R 1}\right) \iint_{\Omega}\left[E X\left(\eta_{T} \cdot \frac{\partial E X\left(v_{1}\right)}{\partial \alpha_{2}}\right)\right]_{\alpha_{2}=E X\left(\varepsilon_{T}\right)} h_{1} h_{3} d \alpha_{1} d \alpha_{3} \]

\[ E X\left(F_{R 3}\right) \iint_{\Omega}\left[E X\left(\eta_{T} \cdot \frac{\partial E X\left(v_{3}\right)}{\partial \alpha_{2}}\right)\right]_{\alpha_{2}=E X\left(\varepsilon_{T}\right)} h_{1} h_{3} d \alpha_{1} d \alpha_{3} \tag{23} \]

where 0 ≤ h₁α₁ ≤ 2, –b ≤ h₃α₃ ≤ b, 0 ≤ α₂ ≤ EX(ε_T), h₁α₁[m], h₃α₃[m], a₂[m], a [m], b [m], and Ω (α₁,α₃) lubrication surface, h₂ = 1, h₁ = h₁(α₁,α₃), and h₃ = h₃(α₁,α₃). When the derivatives of the expected values of random biological fluid velocity components (14) and (16) are substituted into expression (23), considering the derivatives of functions (14) and (17ab), the expected functions of the random friction force components are obtained in the following form of double integrals [21]:

\[ E X\left(F_{R 1}\right)=\iint_{\Omega}\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right)\left[E X\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]\right. \]

\[ \left.-\frac{\int_{0}^{E X\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]} \frac{\alpha_{2} d \alpha_{2}}{E X\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}}{\int_{0}^{E X\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]} \frac{d \alpha_{2}}{\operatorname{EX}\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}}\right] h_{1} h_{3} d \alpha_{1} d \alpha_{3}+ \]

\[ -\iint_{\Omega}\left[\frac{U_{1}\left(\alpha_{1}, \alpha_{3}\right)}{\int_{0}^{E X\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]} \frac{d \alpha_{2}}{\operatorname{EX}\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}}\right] h_{1} h_{3} d \alpha_{1} d \alpha_{3}\tag{24} \]

\[ E X\left(F_{R 3}\right)=\iint_{\Omega}\left(\frac{1}{h_{3}} \frac{\partial E X(p)}{\partial \alpha_{3}}-M_{3}\right)\left[\operatorname{EX}\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]-\frac{\int_{0}^{\operatorname{EX}\left[\varepsilon_{\tau}\left(\alpha_{1}, \alpha_{3}\right)\right]} \frac{\alpha_{2} d \alpha_{2}}{\operatorname{EX}\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}}{\int_{0}^{\operatorname{EX}\left[\varepsilon_{\tau}\left(\alpha_{1}, \alpha_{3}\right)\right]} \frac{d \alpha_{2}}{\operatorname{EX}\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}}\right] h_{1} h_{3} d \alpha_{1} d \alpha_{3}+\]

\[ -\iint_{\Omega}\left[\frac{U_{3}\left(\alpha_{1}, \alpha_{3}\right)}{\left[\int_{0}^{{E x\left[\varepsilon_{T}\left(\alpha_{1}, \alpha_{3}\right)\right]} } \frac{d\alpha_2}{E X\left[\eta_{T}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]}\right.}\right] h_{1} d \alpha_{1} h_{3} d \alpha_{3}-\rho \iint_{\Omega} \frac{1}{h_{1} h_{3}} \frac{\partial h_{1}}{\partial \alpha_{3}}\left\{\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right)^{2} \operatorname{EX}\left(F_{\rho 1}\right)+\right. \]

\[ \left.-2 U_{1}\left(\frac{1}{h_{1}} \frac{\partial E X(p)}{\partial \alpha_{1}}-M_{1}\right) \operatorname{EX}\left(F_{\rho 2}\right)+U_{1}^{2} E X\left(F_{\rho 3}\right)\right\} h_{1} h_{3} d \alpha_{1} d \alpha_{3}\tag{25}\]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], b [m], U₁ = U₁(α₁,α₃) [m/s], and U₃ = U₃(α₁,α₃) [m/s].

The terms, when multiplied by the linear velocity components U₁ and U₃, describe the influence of surface motion on the values of friction forces.

The last two terms on the right-hand side of equation (25), when multiplied by the biological fluid density ρ, indicate the effects of centrifugal forces on the friction forces. Considering equations (23), (25), and (17a, 17b), the EX(F_ρ1) [m⁹/N²s²], EX(F_ρ2)[m⁵/Ns], EX(F_ρ3) [m] form of auxiliary expectancy function occurring in equation (25) and presenting the centrifugal forces [23,26] can be derived as follows:

\[ E X\left(F_{\rho 1}\right) \equiv\left\{E X\left(\eta_{T}\right) \frac{\partial}{\partial \alpha_{2}}\left[A_{\rho 1}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]\right\}_{\alpha_{2}=E X\left(\varepsilon_{T}\right)}= \]

\[ \int_{0}^{E X\left(\varepsilon_{T}\right)} A_{\eta}^{2} d \alpha_{2}-\frac{\int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}} A_{\eta}^{2} d \alpha_{2}\right) d \alpha_{2}}{\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}} \tag{26a} \]

\[ E X\left(F_{\rho 2}\right) \equiv\left\{E X\left(\eta_{T}\right) \frac{\partial}{\partial \alpha_{2}}\left[A_{\rho 2}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]\right\}_{\alpha_{2}=E X\left(\varepsilon_{T}\right)}= \]

\[ \int_{0}^{E X\left(\varepsilon_{T}\right)} A_{\eta}\left(1-A_{s}\right) d \alpha_{2}-\frac{\int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}} A_{\eta}\left(1-A_{s}\right) d \alpha_{2}\right) d \alpha_{2}}{\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}} \tag{26b} \]

\[ E X\left(F_{\rho 3}\right) \equiv\left\{E X\left(\eta_{T}\right) \frac{\partial}{\partial \alpha_{2}}\left[A_{\rho 3}\left(\alpha_{1}, \alpha_{2}, \alpha_{3}\right)\right]\right\}_{\alpha_{2}=E X\left(\varepsilon_{T}\right)}= \]

\[ \int_{0}^{E X\left(\varepsilon_{T}\right)}\left(1-A_{s}\right)^{2} d \alpha_{2}-\frac{\int_{0}^{E X\left(\varepsilon_{T}\right)}\left(\frac{1}{E X\left(\eta_{T}\right)} \int_{0}^{\alpha_{2}}\left(1-A_{s}\right)^{2} d \alpha_{2}\right) d \alpha_{2}}{\int_{0}^{E X\left(\varepsilon_{T}\right)} \frac{1}{E X\left(\eta_{T}\right)} d \alpha_{2}} \tag{26c} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, 0 ≤ α₂ ≤ EX(ε_T), h₁α₁[m], h₃α₃[m], α₂[m], a [m], and b [m]. EX [η_T (α₁,α₂,α₃)] denotes the expected function of the random dynamic viscosity of synovial fluid, which also varies in the direction of gap height. EX [ε_T(α₁,α₃)] indicates the expected function of the variable total gap height.

The expected value of the load carrying capacity C [N] for non-rotational, cooperating, biological surfaces with h₂ = 1, h₁ = h₁(α₁,α₃), h₃ = h₃(α₁,α₃) can be formulated in the curvilinear coordinates in the following form of surface integrals [23,26]:

\[ E X(C)=\iint_{\Omega} E X\left[p\left(\alpha_{1}, \alpha_{3}\right)\right] d \Omega, \mathrm{d} \Omega=h_{1} h_{3} d \alpha_{1} d \alpha_{3} \tag{27} \]

where 0 ≤ h₁α₁ ≤ 2a, –b ≤ h₃α₃ ≤ b, h₁α₁[m], h₃α₃[m], α₂[m], a [m], and b [m]. The load carrying capacity is along the direction vertical to the biological friction nod surface Ω [m²] and in the direction opposite to load W [N].

For rotational, biological surfaces in α₁ direction and non-monotone generating line in α₃ direction, where h₂ = 1, h₁ = h₁(α₁,α₃), and h₃ = h₃(α₁,α₃), the expected value of load carrying capacity (27) in the curvilinear coordinates leads to the following double-integral form [25,28]:

\[ E X(\mathcal{C})=\sqrt{\left[\int_{-b}^{b}\left(\int_{0}^{\alpha_{k}} \operatorname{Ex}\left[p\left(\alpha_{1}, \alpha_{3}\right)\left(\sin \alpha_{1}\right) h_{1} d \alpha_{1}\right]\right) h_{3} d \alpha_{3}\right]^{2}+\left[\int_{-b}^{b}\left(\int_{0}^{\alpha_{k}} \operatorname{Ex}\left[p\left(\alpha_{1}, \alpha_{3}\right)\left(\cos \alpha_{1}\right) h_{1} d \alpha_{1}\right]\right) h_{3} d \alpha_{3}\right]^{2}} \tag{28} \]

for 0 ≤ α₁ ≤ α_k < 2π, –b ≤ h₃α₃ ≤ b, EX(ε_T) = EX [ε_T(α₁,α₃)], where α_k denotes the bio-fluid end coordinate in the circumferential direction α₁ of the femoral head, and 2b is the length in the longitudinal direction of the friction area.

The symbol EX(p) indicates the expected function of the randomly variable hydrodynamic pressure function. Based on Coulomb’s law of friction, and considering the expected function of adhesion force EX(A_D) [N], in a curvilinear orthogonal coordinate system, the dimensionless randomly variable coefficient of friction has the following form [23,26]:

\[ E X(\mu)=\frac{\left|\boldsymbol{e}_{\mathbf{1}} E X\left(F_{R 1}\right)+\boldsymbol{e}_{3} E X\left(F_{R 3}\right)\right|-E X\left(A_{D}\right)}{E X(C)} \tag{29} \]

where e₁ and e₃ are the unit vectors tangential to the circumferential α₁ and longitudinal α₃ directions.

4. Semi-Experimental and Theoretical Methods

Here, a general estimation of random lubricant parameters for two non-rotational, curvilinear, co-operating, biological surfaces is elaborated. First, an estimation of the value of the expected functions of random lubricant liquid velocity, random hydrodynamic pressure, load-carrying capacity, and friction forces for the lubricated, curvilinear, arbitrary, biological, non-rotational surfaces is presented based on the stochastic gap height analysis (Section 2). As mentioned in the analytical solutions presented in the sub-Section 3.2, a non-rotational surface moves with the linear velocity U₁(α₁,α₃)[m/s] in α₁ direction and U₃(α₁,α₃)[m/s] in α₃ direction, while the second bio-surface is stationary and limits a gap filled with a layer of bio-liquid with randomly varying height EX(ε_T).

Using the probability density function f for correction δ₁ of the gap height obtained experimentally, the expected correction value, i.e., EX(δ₁) = m, is determined from formula (2a). The standard deviation σ is calculated from equation (3). The inequality that enables the estimation of the expected values of gap height is based on equations (1) and (2b) and takes the following form [26,27,30]:

\[ (1+m-\sigma) \cdot \varepsilon_{T}\left(\delta_{1}=0\right) \leq E X\left(\varepsilon_{T}\right) \leq(1+m+\sigma) \cdot \varepsilon_{T}\left(\delta_{1}=0\right) \tag{30} \]

where ε_T(δ₁ = 0) = ε₀⋅ε_T1s(α₁,α₃) = ε_T(α₁,α₃,δ₁ = 0) represents the gap height without any random corrections.

For the stationary flow of biological fluid in the friction nod between surfaces or in the joint gap, there is an equal volumetric flow rate in places of the narrowing and expansion of gap height ε_T, in both classic hydrodynamic lubrication and lubrication by squeezing out. The flow rate, as the product of average flow velocity and flow surface, implies an increase (decrease) in velocity v, which is inversely proportional to the narrowing (widening) of the gap height [30]. Using the expected value of gap height correction m and its standard deviation σ, and considering the numerical estimations of the analytical velocity solutions (14) and (16), the following inequalities for the value of the expected velocity function are obtained [23]:

\[ \frac{1}{1+m+\sigma} \leq \frac{E X(v)}{v\left(\delta_{1}=0\right)} \leq \frac{1}{1+m-\sigma}, v\left(\delta_{1}=0\right)=O(U), U \equiv \sqrt{U_{1}^{2}+U_{3}^{2}} \tag{31} \]

where v(δ₁ =0) indicates the velocity without any random corrections of gap height.

The apparent viscosity variation is confirmed by formula (9a) that defines the apparent viscosity η_T. In this formula, velocity v is inversely proportional to the dynamic viscosity of the biological fluid. Thus, using the expected value of gap height correction m and its standard deviation σ, along with the numerical estimations for expression (9a), the following inequalities for the value of the expected viscosity function are observed:

\[ (1+m-\sigma) \leq \frac{E X\left(\eta_{T}\right)}{\eta_{T}\left(\delta_{1}=0\right)} \leq(1+m+\sigma), \eta_{T}\left(\delta_{1}=0\right)=O\left(\frac{\gamma_{\max }+k \cdot A^{-1} T \cdot \ln L}{\delta_{v} \cdot v_{0}}\right) \tag{32a} \]

where η_T(δ_{1 =} 0) indicates the biological liquid viscosity without random corrections of gap height and the susceptibility random changes on the superficial layer of surfaces.

Based on equation (32a), the apparent dynamic viscosity for synovial fluid is calculated, taking into account the following parameter values: A = 10^–15 m², lnL = –50, γ = 2.5 mJ/m², where 0.03 m/s < v₀ < 0.04 m/s; 0.2 < δ_v < 0.6; T = 310 K; and k = 1.380649×10^–23 J/K. Hence, kA^–1 T×lnL = –0.214 mN/m.

Eventually, the apparent viscosity has the following value [20]:

\[ \eta_{T}\left(\delta_{1}=0\right)=\frac{\gamma_{\max }+k \cdot A^{-1} T \cdot \ln L}{\delta_{v} \cdot v_{0}}=\frac{2.5 \frac{\mathrm{mJ}}{\mathrm{m}^{2}}-0.2140 \ldots \frac{\mathrm{m} N}{\mathrm{~m}}}{0.20 \cdot 0.0368 \frac{\mathrm{m}}{\mathrm{s}}}= \]

\[ =\frac{2.5 \frac{m N}{m}-0.2140 \ldots \frac{m N}{m}}{0.20 \cdot 0.0368 \frac{m}{s}}=0.3105 \frac{N \cdot s}{m^{2}}=0.3105 P a \cdot s \tag{32b} \]

Symbols A, L, γ, and k-Boltzmann constant are described and defined in Section 3. It is visible that the boundary surface A between the areas of different phospholipid concentrations has a noteworthy influence on the bio-fluid viscosity variations.

The derived formula (32a) for the apparent viscosity of biological liquid proves that the increments in phospholipid concentration on the joint cartilage surface corresponding to the decrements in the value 0.6 to value 0.2 of the dimensionless coefficient δ_v imply the increments in the apparent dynamic viscosity. The experimental studies [1,8] validate this feature of phospholipids.

The derived formula (32a) for the apparent viscosity of biological liquid proves that the decrements in the flow velocity v₀ of the biological liquid in the joint gap during the lubrication increase the apparent dynamic viscosity of the biological liquid.

According to the fundamental laws of fluid mechanics, the velocity distribution of a viscous liquid flow in the conduit gap has a parabolic shape, while the lowest values of velocity are located on the boundary surfaces of the gap and are restricted to the liquid layer. Experimental studies [6,12] confirm that the highest values of the dynamic viscosity of the lubricated biological liquid are obtained on the laminar boundary layer of the cartilage, which is flowed around by the biological liquid.

Analytical and numerical solutions of the system of equations (4)-(8) and (19a, 19b, 19c) prove that a smaller gap height implies a greater velocity and a lower viscosity of the non-Newtonian liquid. This fact applies to both hydrodynamic classic lubrication and co-operating surface lubrication by squeezing out. The expected function of the hydrodynamic pressure EX (p) estimated numerically from the differential equations (19a, 19b, 19c), with the simultaneous use of results (31) and (32a), leads to a dimensionless pressure quotient in the following interval [30]:

\[ \frac{(1+m-\sigma) \cdot \alpha_{p}}{(1+m+\sigma)^{2}} \leq\left(\frac{E X(p)}{p\left(\delta_{1}=0\right)} \equiv \xi_{p}\right) \leq \frac{(1+m+\sigma) \cdot \beta_{\mathrm{p}}}{(1+m-\sigma)^{2}} \]

\[ p\left(\delta_{1}=0\right)=O\left(\frac{U \cdot \eta_{T}\left(\delta_{1}=0\right)}{\psi \varepsilon_{T}\left(\delta_{1}=0\right)}\right), \psi \equiv \frac{\varepsilon_{T}}{L_{R}} \tag{33} \]

where p(δ_{1 =} 0) indicates the hydrodynamic pressure without random corrections of gap height and L_R [m] is the length of the lubrication region of the non-rotational surface. Dimensionless correction coefficients are dependent on the radial clearance between the two co-operating surfaces and the range of stochastic changes in the dimensionless fraction values in the pressure quotient ξ_p.

The expected values of the load-carrying capacity EX(C) = Ω∙EX(p) are dependent on the pressure p and the numerical estimations obtained from formula (27) for the lubricated surface dΩ = h₁h₃ dα₁dα_3, along with the correction coefficients , which are limited in the following interval [30]:

\[ \frac{(1+m-\sigma) \cdot \alpha_{C}}{(1+m+\sigma)^{2}} \leq\left(\frac{E X(C)}{C\left(\delta_{1}=0\right)} \equiv \xi_{C}\right) \leq \frac{(1+m+\sigma) \cdot \beta_{\mathrm{C}}}{(1+m-\sigma)^{2}} \]

\[ C\left(\delta_{1}=0\right)=\Omega \cdot O\left(\frac{U \cdot \eta_{T}\left(\delta_{1}=0\right)}{\psi \varepsilon_{T}\left(\delta_{1}=0\right)}\right), \Omega \equiv \iint_{\Omega} h_{1} h_{3} d \alpha_{1} d \alpha_{3} \tag{34} \]

where C(δ₁ = 0) represents the load-carrying capacity without random corrections of gap height.

The expected function of friction forces EX(F_R) with components estimated numerically from formulas (24) and (26), with the simultaneous use of inequalities (31–34), enables the representation of the upper and lower values of the friction quotient in the following form:

\[\frac{(1+m-\sigma) \cdot \alpha_{F}}{1+m+\sigma} \leq\left(\frac{E X\left(F_{R}\right)}{F_{R}\left(\delta_{1}=0\right)} \equiv \xi_{F}\right) \leq \frac{(1+m+\sigma) \cdot \beta_{\mathrm{F}}}{1+m-\sigma} \]

\[ F_{R}\left(\delta_{1}=0\right)=O\left(\frac{\Omega \cdot U \cdot \eta_{T}\left(\delta_{1}=0\right)}{\varepsilon_{T}\left(\delta_{1}=0\right)}\right) \tag{35} \]

where F_R(δ₁ = 0) indicates the friction force without random corrections of gap height, and α_F, β_F are the correction coefficients that are dependent on the stochastic changes in the fraction value ξ_F.

The expected function of the coefficient of friction EX(μ) = EX(F_R)/EX(C) expressed by formula (29), without the adhesion forces and with the simultaneous use of inequalities (34) and (35) and the correction coefficients α_μ, β_μ, leads to the estimation in the following interval:

\[ (1+m-\sigma) \cdot \alpha_{\mu} \leq\left(\frac{E X\left(\mu={ }^{F_{R} /} C\right)}{\mu\left(\delta_{1}=0\right)} \equiv \xi_{\mu}\right) \leq(1+m+\sigma) \cdot \beta_{\mu} \]

\[ \mu\left(\delta_{1}=0\right)=O\left(\frac{\varepsilon_{T}\left(\delta_{1}=0\right)}{L_{R}}\right) \tag{36} \]

where μ (δ₁ = 0) indicates the friction coefficient without any random corrections of gap height.

The data for estimations (33)-(36) will be considered in some selected surfaces.

4.1 Example 1

Here, the human hip joint with two co-operating half-spherical rotational surfaces is considered, as depicted in Figure (4a, b, c, d, e, f). The lower surface (bone head) moves in α₁ direction with angular velocity ω and the upper surface (acetabulum) is motionless. Both the spherical surfaces are motionless in the meridian α₃ direction. Hence, the values are U₁ = ω⋅h₁, U₃ = 0 and α₁ = φ, α₂ = r, α₃ = ϑ. Lubrication is performed only on the spherical surface region with the monotone generating line, i.e., inside the intervals: 0 < φ < π, πR/6 < ϑ < πR/2, and ϑ = ϑ₁R meridian coordinates. Thus, the Lamé coefficients are as follows: h₁ = Rsinθ₁, h₃ = 1, and R [m] = L_R which represents the radius of the sphere. Hence, U, Ω, and gap height have the following forms [30]:

\[ U \equiv \sqrt{U_{1}^{2}+U_{3}^{2}}=\sqrt{\left(\omega R \sin \vartheta_{1}\right)^{2}+0^{2}}=\omega R \sin \vartheta_{1}, \Omega \equiv \iint_{\Omega} h_{1} h_{3} d \alpha_{1} d \alpha_{3}=\frac{\sqrt{3} \pi R^{2}}{2} \tag{37} \]

\[ \varepsilon_{T}\left(\varphi, \vartheta_{1}, \delta_{1}\right)=\varepsilon_{0} \varepsilon_{T 1}\left(\varphi, \vartheta_{1}, \delta_{1}\right) \equiv \varepsilon_{0} \varepsilon_{T 1 s}\left(\varphi, \vartheta_{1}\right)\left[1+\delta_{1}\left(\varphi, \vartheta_{1}\right)\right] \tag{38a} \]

\[ \varepsilon_{T}\left(\delta_{1}=0\right)=\varepsilon_{0} \varepsilon_{T 1 s}\left(\varphi, \vartheta_{1}\right)=\varepsilon_{x 1} \cos \varphi \sin \vartheta_{1}+\Delta \varepsilon_{x 2} \sin \varphi \sin \vartheta_{1}-\Delta \varepsilon_{x 3} \cos \vartheta_{1}-R+ \]

\[ +\left[\left(\Delta \varepsilon_{x 1} \cos \varphi \sin \vartheta_{1}+\Delta \varepsilon_{x 2} \sin \varphi \sin \vartheta_{1}-\Delta \varepsilon_{x 3} \cos \vartheta_{1}\right)^{2}+\left(R+\varepsilon_{\min }\right)\left(R+2 D+\varepsilon_{\min }\right)\right]^{0.5} \tag{38b} \]

The symbol ε₀ denotes the characteristic dimensional value of the joint gap height. The center of the spherical bone head is assumed to be at O (0,0,0) and the center of the spherical surface of the acetabulum is at O₁ (x–Δε_x, y–Δε_y, z+Δε_z). Eccentricity has the value of D, as shown in Figures 4a and 4c.

Figure 4 Gap height between randomly deformable spherical cartilage surface: a) the geometrical simulation of hip bone head and sleeve, b) random variations correction δ of the gap height ε_T, c) eccentricity D and center of the acetabulum O₁, d) view of the hip joint gap, e) vertical section of the hip joint in the plane of ante distortion angle, f) phospholipid (PL) bilayer on the internal sleeve (acetabulum) and bonehead surface; R-radius of the bonehead, ε_min-the least value of the gap height, C-capacity, Ω-load (Figures are based on the author`s studies).

4.2 Example 2

Here, the human elbow joint with two co-operating, cylindrical, rotational surfaces is considered. The lower surface (bone head) moves in α₁ direction with angular velocity ω and the upper surface (acetabulum) is motionless. Both the cylindrical surfaces are motionless in the longitudinal α₃ direction.

Hence, the values are U₁ = ω⋅h₁, U₃ = 0 and α₁ = φ, α_{2 =} r, and α₃ = z. Lubrication is performed on the cylindrical surface region with the monotone generating line α₃ = z inside the intervals: 0 < φ < π/2, -b < z < +b, and z longitudinal coordinate. Thus, the Lamé coefficients are as follows: h₁ = R, h₃ = 1, where R [m] = L_R represents the radius of the cylindrical bone. Hence, U = ωR and Ω = πRb. Geometrical simulation of the elbow cylindrical joint is illustrated in Figure 5a.

4.3 Example 3

In this particular case, it is assumed that the thin biological liquid layer is restricted by the two cooperating non-rotational yet parallel biological planes in the rectangular Cartesian coordinates α₁ = x₁, α₂ = x₂, and α₃ = x₃. The Lamé coefficients have the form: h₂ = 1 and h₁ = 1, h₃ = 1, and 0 < x₁ < a, 0 < x₃ < b ≡ L_R. For example, this case is valid for the jump joint and in some regions of the collar bone, blade bone, and skin-tight dress regions. The lower bone (skin) surface moves in α₃ = x₃ direction with a linear velocity U₃ and the upper bone (tight fitting dress) surface is motionless. Both the parallel surfaces are motionless in the meridian α₁ = x₁ direction and, hence, U₁ = 0 and U = U₃, Ω = ab. The two non-rotational parallel plane surfaces are highlighted in Figure 5b.

4.4 Example 4

This case is valid for two movable, non-rotational and non-parallel surfaces represented in Figure 5c, where h₁ = h₁(α₁,α₃), and h₃ = h₃(α₁,α₃).

Figure 5 Geometry of two co-operating surfaces: a) the rotational region of the cylindrical elbow joint, b) movable, non-rotational parallel planes, c) movable, non-rotational and non-parallel surfaces; 1-collagen region, 2-PL bilayer, 3-hydrated phospholipid, 4-synovial liquid, 5-hydrated sodium ions, R-repulsion force, W-load, P_p-the force of the hydrodynamic pressure, ω- angular velocity (Figures are elaborated based on the author’s studies).

Further, it is important to experimentally obtain the density function f with expected value m and standard deviation σ of the gap height between two variable cooperation surfaces. The calculation also needs to be done for the numerical solutions of the stochastic partial differential equations (4)-(9) and the correction coefficients α_p, α_C, α_F, α_μ, β_p, β_C, β_F, β_μ. The realistic real estimation values of the dimensionless quotients of pressure, capacity, friction forces, friction coefficients, ξ_p, ξ_C, ξ_F, ξ_μ, are presented in estimations (33)-(36).

5. Results of Measurements

The real probability density function f of a random gap height change for the spherical hip (h) joint was measured (Refer to example 1 and Figure 4a, 4b, 4c, 4d) for the cylindrical elbow (e) gap height (Example 2 and Figure 5a), and for the two co-operating, non-rotational, biological planes occurring in the jump (j) joint (Refer to examples 3 and 4 and Figure 5b and 5c). The expected value m of the gap height and the standard deviations σ is considered to obtain an estimation of the analytical bio-tribology parameters (33–36) for each surface in concrete standard deviation intervals. The experimental measurement [1,5,12] of the aforementioned joint gap height with radial clearance ε₀ in the range of 2 µm to 10 µm proved that the dimensionless random variable of corrections for the gap height, marked with the symbol δ₁ for hip (h), elbow (e), and jump (j) surfaces, is most often manifested by two characteristic types of probability density functions, which typically are symmetrical and anti-symmetrical functions. Symmetrical density functions f_s of the correction parameters for (h), (e), and (j) surfaces occur less frequently than the anti-symmetrical functions (about 12 times in 100 measurements with a probability of P_s = 0.12). These are characterized by the symmetric distribution of random probability changes in terms of gap height variations. In this case, the probability of a random increase in the joint height gap is equal to the probability of a random decrease.

Among the frequently occurring unsymmetrical correction parameters, there are two types of density functions that describe the random variations in gap height. The first type concerns function f_N, where the probability of a random increase in the gap height dominates over the probability of random decrease in gap height. The second type is the function f_n, where the probability of random decrease in gap height dominates over the probability of the random increase in gap height. Both f_N and f_n are considered for surfaces (h), (e), and (j). In both cases of anti-symmetrical probability density distribution function f_N, each one of them occurs approximately 22 times out of 100 measurements with a probability of P_N = 0.22. In the next two cases of anti-symmetrical probability density distribution function f_n, each one of them occurs 22 times in 100 measurements with a probability of P_n = 0.22. The measurement multiplicity can be computed as 12 + 2 × 22 + 2 × 22 = 100, creating a probabilistically complete system of events [1,26].

In order to determine the final average expected value and the standard deviation of the hip gap height in the form m^h, σ^h, Appendix 3 (Figure A3.1 and figure A3.2a, b, c, d) highlights the probability density functions f_s, f_N, and f_n of random changes in the gap height between hip joint spherical surfaces (h) obtained in a set of experiments. The upper index h indicates that the considered value refers to the hip spherical surface.

The final measured values m^h, σ^h mentioned in Appendix 3 are considered for the current research work (Equations (A3.3). Based on the inequalities (30), (31), and (32a), the expected function for the spherical joint gap height at the speed of the synovial fluid and its viscosity together with probability 0.6708 ≤ P ≤ 1.000 possess values in the following intervals:

\[ 0.6120 \varepsilon_{T}\left(\delta_{1}=0\right)=\left(1+m^{h}-\sigma^{h}\right) \varepsilon_{T}\left(\delta_{1}=0\right) \leq E X\left(\varepsilon_{T}\right) \leq\left(1+m^{h}+\sigma^{h}\right) \varepsilon_{T}\left(\delta_{1}=0\right) \]

\[ =1.3879 \varepsilon_{T}\left(\delta_{1}=0\right) \tag{39} \]

\[ 0.7205 v\left(\delta_{1}=0\right)=\left(1+m^{h}+\sigma^{h}\right)^{-1} v\left(\delta_{1}=0\right) \leq E X(v) \leq\left(1+m^{h}-\sigma^{h}\right)^{-1} v\left(\delta_{1}=0\right) \]

\[ =1.6339 v\left(\delta_{1}=0\right) \tag{40} \]

\[ 0.6120 \eta_{T}\left(\delta_{1}=0\right)=\left(1+m^{h}-\sigma^{h}\right) \eta_{T}\left(\delta_{1}=0\right) \leq E X\left(\eta_{T}\right) \leq\left(1+m^{h}+\sigma^{h}\right) \eta_{T}\left(\delta_{1}=0\right) \]

\[ =1.3879 \eta_{T}\left(\delta_{1}=0\right) \tag{41} \]

Substituting the measured values of the expression 1+m^{h} \pm \sigma^{h} into the inequalities (33)-(36), the following upper and lower estimations of quotients \xi_{p}^{h}, \xi_{C}^{h}, \xi_{F}^{h}, \xi_{\mu}^{h}; are obtained, with upper index (h) for the hip spherical co-operating joint surfaces:

\[ 0.3177 \cdot \alpha_{p}^{h} \leq\left(\xi_{p}^{h} \equiv \frac{E X(p)}{p\left(\delta_{1}=0\right)}\right) \leq 3.7055 \cdot \beta_{p}^{h}, \quad 0.3177 \cdot \alpha_{C}^{h} \leq\left(\xi_{C}^{h} \equiv \frac{E X(C)}{C\left(\delta_{1}=0\right)}\right) \leq 3.7055 \cdot \beta_{C}^{h} \]

\[ 0.4409 \cdot \alpha_{F}^{h} \leq\left(\xi_{F}^{h} \equiv \frac{E X\left(F_{R}\right)}{F_{R}\left(\delta_{1}=0\right)}\right) \leq 2.2678 \cdot \beta_{F}^{h}, 0.6120 \cdot \alpha_{\mu}^{h} \leq\left(\xi_{\mu}^{h} \equiv \frac{E X(\mu)}{\mu\left(\delta_{1}=0\right)}\right) \leq 1.3879 \cdot \beta_{\mu}^{h} \tag{42}\]

Based on the probability density functions f_s, f_N, and f_n obtained experimentally through gap height random changes between elbow joint cylindrical surfaces (e) and plane jump surfaces (j), the average expected value of the gap height and standard deviation of the gap height are analogously obtained in the forms m^{e}, \sigma^{e}, m^{j}, \sigma^{j}, respectively. Substituting the measured values of the expressions 1+m^{e} \pm \sigma^{e} and 1+m^{j} \pm \sigma^{j} into the inequalities (33)–(36), the following upper and lower estimations of quotients \cdot \xi_{p}^{e}, \xi_{C}^{e}, \xi_{F}^{e}, \xi_{\mu}^{e} ; \xi_{p}^{j}, \xi_{C}^{j}, \xi_{F}^{j}, \xi_{\mu}^{j} for elbow and jump co-operating joint surfaces are obtained, with the upper index (e) and (j), respectively:

\[ 0.3277 \cdot \alpha_{p}^{e} \leq\left(\xi_{p}^{e} \equiv \frac{E X(p)}{p\left(\delta_{1}=0\right)}\right) \leq 3.6055 \cdot \beta_{p}^{e}, \quad 0.3277 \cdot \alpha_{C}^{e} \leq\left(\xi_{C}^{e} \equiv \frac{E X(C)}{C\left(\delta_{1}=0\right)}\right) \leq 3.6055 \cdot \beta_{C}^{e} \]

\[ 0.4509 \cdot \alpha_{F}^{e} \leq\left(\xi_{F}^{e} \equiv \frac{E X\left(F_{R}\right)}{F_{R}\left(\delta_{1}=0\right)}\right) \leq 2.1678 \cdot \beta_{F}^{e}, 0.6220 \cdot \alpha_{\mu}^{e} \leq\left(\xi_{\mu}^{e} \equiv \frac{E X(\mu)}{\mu\left(\delta_{1}=0\right)}\right) \leq 1.2879 \cdot \beta_{\mu}^{e} \tag{43a} \]

\[ \text { 0.3377. } \alpha_{p}^{j} \leq\left(\xi_{p}^{j} \equiv \frac{E X(p)}{p\left(\delta_{1}=0\right)}\right) \leq 3.505 \cdot \beta_{p}^{j}, \quad 0.3377 \cdot \alpha_{c}^{j} \leq\left(\xi_{c}^{j} \equiv \frac{E X(C)}{C\left(\delta_{1}=0\right)}\right) \leq 3.5055 \cdot \beta_{C}^{j} \]

\[ 0.4609 \cdot \alpha_{F}^{j} \leq\left(\xi_{F}^{j} \equiv \frac{E X\left(F_{R}\right)}{F_{R}\left(\delta_{1}=0\right)}\right) \leq 2.0678 \cdot \beta_{F}^{j}, 0.6320 \cdot \alpha_{\mu}^{j} \leq\left(\xi_{\mu}^{j} \equiv \frac{E X(\mu)}{\mu\left(\delta_{1}=0\right)}\right) \leq 1.1879 \cdot \beta_{\mu}^{j} \tag{43b} \]

Equations (42) and (43a, b) represent the upper and the lower estimations of the dimensionless quotients of the stochastic functions EX(p), EX(C), EX(F_R), and EX(μ) in the numerator, to the non-stochastic functions p(δ₁ = 0), C(δ₁ = 0), F_R(δ₁ = 0), and μ(δ₁ = 0) in the denominator, for pressure, load-carrying capacity, frictional forces, and the coefficients of friction. The quotient’s numerators are obtained from the numerical solutions of the stochastic system (4)-(9), with operator EX, assuming the stochastic changes of gap height, velocity, and dynamic viscosity. The quotient’s denominators are determined from the system of equations (4)-(9), excluding the random operator EX. The upper and the lower values of the quotients:

\[ \xi_{p}^{h}, \xi_{C}^{h}, \xi_{F}^{h}, \xi_{\mu}^{h} ; \xi_{p}^{e}, \xi_{C}^{e}, \xi_{F}^{e}, \xi_{\mu}^{e} ; \xi_{p}^{j}, \xi_{C}^{j}, \xi_{F}^{j}, \xi_{\mu}^{j} \tag{43c} \]

of pressure, capacity, friction forces, and friction coefficients described in (42) and (43a, b) need the value of corrections coefficients:

\[ \alpha_{p}, \alpha_{C}, \alpha_{F}, \alpha_{\mu}, \beta_{\mathrm{p}}, \beta_{\mathrm{C}}, \beta_{\mathrm{F}}, \beta_{\mu} \tag{43d} \]

for the surfaces h, e, and j. Therefore, the numerical solutions of the stochastic partial differential equations (4)-(9) for various radial clearances between two co-operating surfaces and the different ranges of stochastic changes are determined in the next section.

6. Sketch of Numerical Calculations

The diagram of the half-analytical system (14)-(19), as well as the numerical solutions of the system of equations (4)-(9), primarily concerns the stochastic and non-stochastic (excluding operator EX) non-linear partial differential integral hydrodynamic equations for spherical, cylindrical surfaces and planes. For this particular case, the numerical calculations are elaborated for equation (19), which determines two-dimensional hydrodynamic pressure p, and equation (21), which describes the three-dimensional temperature distribution T. These equations, when converted to spherical (φ, r, ϑ), cylindrical (φ, r, z), or rectangular (x, y, z) coordinates, are always mutually coupled through the three-dimensional dynamic viscosity function η_T represented by formula (9a and 9b), variable in the direction of the gap height and dependent on temperature. The coupling of the equations is based on the simultaneous dependence of pressure on temperature and temperature on pressure. The numerical integration of the coupled equations was carried out using the method of a convergent sequence of successive approximate solutions. In the first step of approximation, a constant dimensionless viscosity is accepted with the value of one. The hydrodynamic pressure for constant viscosity is determined, followed by the temperature distribution. In the second step of the calculation, the obtained pressure and temperature values are used to determine the variable 3D viscosity value. Hydrodynamic pressure and temperature distributions were calculated using the obtained variable viscosity. In the subsequent steps of the calculation, the described procedure is iterated until the obtained sequence of pressure and temperature functions is convergent with the boundary function of pressure and temperature. The convergence of such a sequence to the boundary function of pressure and temperature, regardless of the value of the accepted constant viscosity at the beginning of the calculation process, suggests that the solution of the coupled equation system has been obtained correctly. The calculations were carried out using the finite difference method, using some derived numerical procedures and professional Mathcad 15 software. The pressure and temperature from partial differential equations (19) and (21) for the selected coordinates are obtained by performing numerical calculations. The frictional force is calculated from formulas (24) and (25), followed by the load-carrying capacity C and the coefficient of friction µ, with the omission of adhesion forces (A_D = 0), based on formulas (27) and (29).

7. Final Results after Analytical, Numerical, and Experimental Research

In this section, the real and final estimations for the correction of upper and lower values of the dimensionless quotients (fractions) of pressure, capacity, friction forces, and coefficients of friction for the selected co-operating surfaces occurring in human joints are described.

The real dimensionless estimation values of fractions with upper indexes h, e, and j for human hip, elbow, and jump joint lubrication parameters, respectively, are presented in (42) and (43a, 43b, 43c), with the quotients of the aforementioned solutions obtained from the numerical solutions of (4)-(9). Based on the numerical solutions in Section 6, the correction coefficients for surfaces h, e, and j are determined. Hence, the lower and upper estimations of the aforementioned quotients are obtained in the following form:

\[ \xi_{p \min }^{h} \leq \xi_{p}^{h} \leq \xi_{p \max }^{h} ,\xi_{C \min }^{h} \leq \xi_{c}^{h} \leq \xi_{c \max }^{h}, \xi_{F m i n}^{h} \leq \xi_{F}^{h} \leq \xi_{F m a x}^{h}, \xi_{\mu m i n}^{h} \leq \xi_{\mu}^{h} \leq \xi_{\mu m a x}^{h} \tag{44a} \]

\[ \xi_{p \min }^{e} \leq \xi_{p}^{e} \leq \xi_{p \max }^{e}, \xi_{\operatorname{Cmin}}^{e} \leq \xi_{C}^{e} \leq \xi_{\operatorname{Cmax}}^{e} \xi_{F m i n}^{e} \leq \xi_{F}^{e} \leq \xi_{F \max }^{e}, \xi_{\mu m i n}^{e} \leq \xi_{\mu}^{e} \leq \xi_{\mu \max }^{e} \tag{44b} \]

\[ \xi_{p \min }^{j} \leq \xi_{p}^{j} \leq \xi_{p \max }^{j}, \xi_{C \min }^{j} \leq \xi_{C}^{j} \leq \xi_{C \max }^{j}, \xi_{F m i n}^{j} \leq \xi_{F}^{j} \leq \xi_{F \max }^{j} , \xi_{\mu \min }^{j} \leq \xi_{\mu}^{j} \leq \xi_{\mu \max }^{j} \tag{44c} \]

The fraction ξ_min (with lower index p, C, F, μ and upper index h, e, j) is calculated by assuming the stochastic variations of the lubrication parameters for which the fraction value is the lowest. The fraction ξ_max is calculated by assuming the stochastic changes for which the fraction value is the highest.

Based on the numerical calculations of the semi-analytical solutions of differential equations (19), (20), and (24)-(28) restricted for the rotational spherical hip surfaces, the values of correction coefficients (43d) α_p, α_C, α_F, α_μ, β_p, β_C, β_F, β_μ as well as the maximum and minimum extreme values \xi_{\max }^{h}, \xi_{\min }^{h} were calculated through an example. These values are listed in Table 1 as the quotients (42) of stochastic values to non-stochastic values for pressure, load carrying capacity, frictional forces, and the coefficients of friction.

The above-mentioned correction coefficients depend on the radial clearances (the characteristic values of gap height ε₀) and the dimensionless range of stochastic changes in the solutions of equations (4)-(9) that represent the random expected values of gap height for the specific numerical calculations of spherical hip joints.

Moreover, the numerical results data are restricted to the random hip gap height between two co-operating, rotational, spherical surfaces. Figure 6 (a, b, c) represents the variations in the dimensional values of load-carrying capacity and frictional forces as well as the changes in the dimensionless coefficient of friction for stochastic changes in gap height versus the radial clearance between two co-operating, spherical, bio-surfaces occurring in the hip joint.

Table 1 Extreme values: $\xi_{p \max ^{\prime}}^{h} \xi_{p \min }^{h}$, $\xi_{C \max ^{\prime}}^{h} \xi_{C \min }^{h}$, $\xi_{F \max ^{\prime}}^{h} \xi_{F \min }^{h}$, $\xi_{\mu \max ^{\prime}}^{h} \xi_{\mu \min }^{h}$ with upper index h and lower index (p, C, F, μ), describing the ratio of solution (pressure p, capacity C, friction forces F, and friction coefficient μ). The calculation is based on the system of equations (4–9) for stochastic changes in the spherical hip joint gap height 0.6120 < EX(ε_T)/ε_T (δ₁ = 0) < 1.3879 as mentioned in (A3.3), with a radial clearance 2 µm ≤ ε₀ ≤ 10 µm, to the same solutions, by assuming a gap height of hip joint (h) without random change.

Figure 6 Random lubrication parameters versus radial clearance 2 µm ≤ ε₀ ≤ 10 µm between two spherical hip joint (h) bio-surfaces, for stochastic changes in gap height 0.612·ε_T and 1.388·ε_T, in the absence of stochastic changes 1.000·ε_T; a) load-carrying capacity, b) friction forces and c) coefficients of friction. (Figures are based on the author’s numerical calculations).

Taking into account the density function f_N, Figure 7a suggests that the load-carrying capacity decreases about 44 percent (attained value 0.558) for PS and PC with a probability of P = 0.604, and increases about 239 percent (attained value 2.393) with a probability of P = 0.729 for PS and a probability of P = 0.788 for PC, compared to 100% for the dimensionless quotient ξ_C = 1 obtained without random changes. Inside the standard deviation interval on the axis δ₁, the expected value m₁ = 0.25 has been calculated for the corrected gap height (m₁+1)⋅ε_T = 1.25ε_T, where the expected capacity attained the increase of about 55% (1.55) with a probability of P = 0.916 for PC, compared to the capacity without random effects. The expected value m₂ = 0.125 for the corrected gap height (m₂+1)⋅ε_T = 1. 125ε_T is found, where the expected capacity attained increases of about 34% (1.34) with a probability of P = 0.937 for PS, compared to the capacity without random effects for ξ_C = 1.

Figure 7b illustrates that for the density function f_N in the radial clearance interval from 2 mm to 10 mm, the dimensionless expected value of capacity decreases from 1.55 to 1.50 with a constant probability of P = 0.916 for PC, and from 1.34 to 1.31 with a constant probability of P = 0.937 for PS, compared to the dimensionless capacity without random effects for .

Figure 7 Dimensionless random capacity quotient for spherical surface after the random density function f_N: a) versus the random variable gap height correction δ₁ and probability P, within the standard deviation interval of gap height (1+m₂-σ<1+m₂+σ) ε_T, i.e., (0.74 < 1.51)×ε_T for PS and (0.86 < 1.63)×ε_T for PC with a radial clearance of 2 μm; b) versus radial clearance from 2 μm to 10 μm, with a constant probability of P = 0.916 and gap height (m₁+1) ε_T = 1.25ε_T for PC and with a constant probability of P = 0.937 and gap height (m₂+1) ε_T = 1.125ε_T for PS.

For the density function f_n, for PS, the expected value m₂ = –0.125, the corrected value of gap height = 0,875ε_T, and the standard deviation interval is (0.49 < 1.26)⋅ε_T, while for PC, the expected value m₁ = –0.25, the corrected value of gap height = 0.75ε_T, and the standard deviation interval is (0.37<1.13)⋅ε_T. Variations in the expected capacity values obtained for density functions f_N and f_n are the same.

Numerical calculations are based on the influence of phospholipid concentration on the spherical hip surface of articular cartilage (A = 10^–15 m²); the electromagnetic field (M_i = 0, M_T = 0) was omitted. For pressure and temperature, boundary conditions (20) and (22a, b, c), and density functions of variable random joint gap height ε_T described by the formulas (A3.1) and (A.3.2a, b) were assumed. Besides, the following were assumed: ln(L) = –50, the characteristic value of synovial fluid flow velocity in the joint gap 0.03 m/s < v₀ < 0.04 m/s, the interfacial energy 0.1 mN/m < γ < 4 mN/m, the coefficient of collagen fibre concentration 0.2 < δ_v < 0.6, the characteristic value of synovial fluid dynamic viscosity η₀ = 0.300 Pas, radius of the sphere of femoral head R = 0.0265 m, the characteristic value of ambient temperature T₀ = 305 K, thermal conductivity of synovial fluid k = 0.13 W/mK, synovial fluid density ρ = 850 kg/m³, angular velocity of bonehead ω = 5s^–1, lubrication area Ω = 0.0004 m², flow index n = 1.1, and eccentricity of bone head: Δε_x1 = Δε_x2 = Δε_x3 = 5μm.

8. Conclusions for Rotational and Movable Non-Rotational Co-operating Biological Surfaces

This paper justified the following key results, which have not been previously obtained in the contemporary domain so far. A few papers on random lubrication theory are referred for the biological joints with co-operating movable, rotational, and non-rotational surfaces coated with phospholipids bilayers.

8.1 Capacity Increments after Bio-liquid Dynamic Viscosity Variations Across the Thin Layer

The variation, which has not been considered so far, of bio-liquid dynamic viscosity across the very thin lubricating bio-liquid layer between two co-operating multimode, rotational and non-rotational bio-surfaces in human joints accounts for about 12 percent of the variations in the joint load-carrying capacity.

8.2 Bio-fluid Dynamic Viscosity Increments after Phospholipids Concentration Increments

Results of the numerical calculation performed on the random apparent dynamic viscosity of the lubricated bio-liquid [(9a) and (32a, b)] lead to a conclusion that a 10 percent increment in the phospholipid concentration at surface A, with increments of different multimode boundary areas of the phospholipid concentration on the non-rotational surfaces, causes approximately 7 percent increments in the apparent dynamic viscosity of the biological liquid.

8.3 Random Estimation of Lubrication Parameters for Cooperating Non-rotational Bio-surfaces

The upper and lower estimations of the expected values of the joint gap height (30), the synovial fluid flow velocity (31), and its apparent viscosity (32), capacity, friction forces were determined based on the stochastic equations of the hydrodynamic lubrication theory of movable, non-rotational bio-surfaces and their semi-analytical-numerical solutions. These estimations are functions of the co-operating surface shapes, standard deviation σ, the expected value m of gap height, probability P, and values of the density function f of gap height.

8.4 Random Estimation of Bio-liquid Dynamic Velocity for Arbitrary Bio-surfaces Lubrication

The particular solutions (14) and (16) for biological fluid velocity obtained for arbitrary rotational and movable, non-rotational, co-operating biological surfaces imply the stochastic estimation (31), demonstrating about 28%–30% decrease (with P = 0.60) and 60%–63% increase (with 0.72 < P < 0.78) in the biological fluid flow velocity v, compared to the velocity obtained without any random changes in the gap height for δ₁ = 0. These changes are valid in the standard deviation interval of the gap height. However, the most realistic chances of the expected velocity attained a value of about 15% (for PS) with a probability of P = 0.93 and about 25% (for PC) with a probability of merely 0.91.

8.5 Random Estimation of Bio-liquid Dynamic Viscosity for Arbitrary Bio-surfaces Lubrication

The particular random solutions 9 (a, b) for biological fluid viscosity obtained for arbitrary rotational and movable, non-rotational, co-operating biological surfaces imply the stochastic estimation (32), demonstrating about 38% decrease (with P = 0.60) and 38% increase (with 0.72 < P < 0.78) in the fluid viscosity η_T, compared to the viscosity obtained without any random changes in the gap height for δ₁ = 0. These changes are valid within the deviation interval of the gap height. However, the most realistic chances of the expected viscosity attained the value of about 9% (for PS) with a probability of P = 0.93 and about 15% (for PC) with a probability of merely 0.91.

8.6 The Capacity and Friction Effects vs Random Gap Height Changes for Arbitrary Bio Surfaces

In the considered spherical, elbow, and jump joints with non-rotational parallel and non-parallel co-operating surfaces, the load-carrying capacity, friction forces, and friction coefficient decrease if the joint gap height increases after the probability effects. These decreases are about 5% smaller for the jump joints with non-rotational plane surfaces, compared to the decreases for the spherical surfaces. For spherical hip joint surfaces, the results are illustrated in Figure 6.

8.6.1 The Capacity vs random Gap Height Variations Between Cooperating Bio-surfaces

The joint capacity obtained without any random changes in the gap height is about 66 percent lower than the capacity for the minimum random gap height, and about 30 percent higher than the capacity for the maximum random gap height.

8.6.2 The Friction Force vs Random Gap Height Variations Between Cooperating Bio-surfaces

The friction forces determined in the joint without any random changes in the gap height are about 12% lower than the friction forces for the maximum random gap height and about 10% higher than the friction forces for the minimum random gap height.

8.6.3 The Friction Coefficient vs Random Gap Height Variations Between Cooperating Bio-surfaces

The friction coefficient obtained in the human joint without any random changes in the gap height is about 70% lower than the friction coefficient for the maximum random gap height, and approximately 50% higher than the friction coefficient for the minimum random gap height.

The differences between the numerical values of capacity, friction forces, and friction coefficients for spherical and cylindrical surfaces, along with the mentioned respective measured values [1,6,8] determined in this paper, attained the values in an interval from about 5 to 9 percent.

9. Discussion

This paper discusses a stochastic flow of synovial fluid with non-Newtonian properties in a biological gap between two isotonic, curvilinear, movable, non-rotational living bio-surfaces. Deformations of susceptible surfaces are caused by random changes in the superficial layer, stochastic loading of the surface, as well as the genetic and volumetric growth of the living articular tissue. The concept of iso-osmosis or isotonic refers to the layers of biological bodies remaining in an osmotic balance with respect to each other. Isotonic and impermeable biological membranes include the lipid bilayer, which eliminates the flow across the layer. Nanometer channels can be treated as an exception for transporting certain types of liquid [30,31].

For the case analyzed in this work, the influence of the physical properties of the articular tissue coated with a phospholipid bilayer and probabilistic gap height variations on the viscosity of the synovial fluid in the area of the boundary layer of the surface flowed around are elaborated. This paper concludes that the variations in the bio-fluid dynamic viscosity depend on the susceptibility of the co-operating movable, non-rotational surfaces, as well as on the probability effects. The research on random changes in the gap height and biological (synovial) fluid viscosity enables the determination of the upper and lower estimations of random changes in hydrodynamic pressure, frictional forces, and the coefficients of friction occurring within the standard deviation intervals. Within the deviation intervals, the localization of the real value of the expected hydrodynamic pressure, load carrying capacity, friction force, friction coefficient, and its probability are demonstrated.

The tools and methods represented in this paper refer to the hydrodynamic pressure obtained in the thin sweat layer between the skin and the tight-fitting gymnastic sports tracksuit. According to the author, this theoretically justifies some experimental results mentioned in the paper [31]. For example, in paper [31,32], it was proved that some training effects are better in the case when a tight dress was used, compared to the results obtained for a loose dress where the height gap between the skin and the dress was large. The experimental conclusion obtained in [32] is as follows: the training of running or walking implemented with a tight-fitting dress leads to increments in the dynamic viscosity of sweat and synovial fluid and alters the internal energy accumulated in the human muscles and cartilage. These results accelerate the slimming (weight losing) process along with reducing body weight.

Symbols and Notations

α_H [J]-protons energy activity,

a [m]-length in the width direction of the friction area,

b [m]-half length in the longitudinal direction of the friction area,

c_c [mol/mm³]-collagen fiber concentration in bio-fluid,

e_ijk [1]-tensor Levi-Civita for i, j, k = 1, 2, 3.,

e [1]-this symbol as upper index assigns the quantity of elbow cylindrical surface,

f [1]-probability density function,

f_s [1]-symmetrical probability density function,

f_n [1]-density function with probability decreases dominate over increases,

f_N [1]-density function for probability increases dominate over decreases,

h_i [1]-Lame Coefficients for i = 1, 2, 3.,

h_i⋅α_i [m]-length coordinate in α_i directions for i = 1, 2, 3.,

h [1]-this symbol as upper index assigns the quantity of the hip spherical surface,

j [1]-this symbol as upper index assigns the quantity of the jump joint plane surface,

k [J/K]-Boltzmann constant,

m [1]-expected value of dimensionless probabilistic corrections,

m_s [1]-expected value for f_s,

m_n [1]-expected value for f_n,

m_N [1]-expected value for f_N,

n [1]-dimensionless flow index,

p [Pa]-pressure,

p_A [Pa]-atmospheric pressure,

pH [1]-dimensionless power hydrogen ion concentration,

s_c [mol/m²]-concentration of phospholipid particles,

t [s]-time,

v [m/s]-biological fluid velocity vector,

v_i [m/s]-biological fluid velocity components in α_I directions for i = 1, 2, 3.,

v₀[m/s]-characteristic value of bio-fluid velocity,

x_i [m]-rectangular coordinates for i = 1, 2, 3.,

A [m²]-the region of areas of different phospholipids molecules concentration,

A_D [N]-adhesion force,

A_η[m⁴/Ns]-auxiliary function,

A_p [m⁶/Ns³]-auxiliary function,

A_ρ1[m¹²/N³s³]-auxiliary function,

A_p2[m⁸/N²s²]-auxiliary function,

A_p3[m⁶/Ns]-auxiliary function,

A_s [1]-dimensionless auxiliary function,

B [T =kg/s²A]-magnetic induction vector in bio-fluid,

C [N]-load carrying capacity of joint,

EX-operator of the expected function,

E [V/m]-electric intensity vector with components (E₁, E₂, E₃),

Fr-topological operator limiting boundary set,

F_Ri [N]-friction force components in α_i direction for i = 1, 2, 3.,

H [A/m]-the magnetic intensity vector with components (H₁, H₂, H₃),

I₁(Θ)[s^–1]-the first invariants of shear rate components,

I₂(Θ)[s^–2]-the second invariants of shear rate components

J [A/m²]-electric current density,

K_a [J]-acid equilibrium constant whereas pKa [1],

K_b [J]-base equilibrium constant whereas pKb [1],

L [1]-dimensionless auxiliary function,

L_k [1]-dimensionless auxiliary function,

M_con [Pasⁿ]-fluid consistency coefficient,

M_i [N/m³]-electrical terms for i = 1, 3.,

N [A/m]-the magnetization vector,

N_A [1/mol]-Avogadro number (= 6.024∙10²³),

P [1]-the probability,

P_n [1]-probability value for density function f_n ,

P_N [1]-probability value for density function f_N,

P_p [N]-the force of hydrodynamic pressure

PC-Phosphati-dylcholine,

PS-Phosphati-dylserine,

PL-Phospholipid bilayer,

R_g [J/Kmol]-gas constant (= 8.3144598),

R [m]-radius of the spherical bone head,

R[N]-repulsion force,

S [Pa]-the stress tensor in biological fluid,

T [K]-temperature,

U_i [m/s]-linear velocity components in α_i direction for i = 1, 3.,

We [°]-wettability,

W [N]-load,

Θ[1/s]-strain tensor,

Θ_ij [1/s]-strain tensor components,

Ξ[A/mK]-the first derivative of the magnetization vector with respect to temperature,

Φ_F [W/m³]-dissipation of energy,

Ω[m²]-lubrication region,

α₁-width (circumferential) direction of the non-rotational (rotational) surfaces,

α₃-longitudinal direction,

α₂-gap height direction,

α_p [1]-lower correction of pressure quotient estimation,

α_C [1]-lower correction of capacity quotient estimation,

α_F [1]-lower correction of friction forces quotient estimation,

α_μ[1]-lower correction of friction coefficient quotient estimation,

β_p[1]-upper correction of pressures quotient estimation,

β_C[1]-upper correction of capacity quotient estimation,

β_F[1]-upper correction of friction forces quotient estimation,

β_μ[1]-upper correction of friction coefficient quotient estimation,

γ[J/m^{2 =}N/m]-interfacial energy,

δ₁[1]-dimensionless random variable corrections,

δ_ν[1]-dimensionless coefficient determining the collagen fibers concentration,

δ_E [m²/V²]-coefficient describing influence of electrical intensity on dynamic viscosity,

δ_ij [1]-dimensionless Kronecker Delta symbol,

δ[1]-unit tensor,

ε₀[m]-dimensional, characteristic gap height value,

ε_T [m]-dimensional total gap height value,

ε_T1[1]-dimensionless total gap height value,

ε_T1s [1]-dimensionless gap height limited by the nominally smooth surfaces,

η_T [Pas]-apparent dynamic viscosity,

κ[W/mK]-thermal conductivity coefficient for biological fluid,

μ_m [mkgs^–2A^–2]-the magnetic permeability coefficient of biological fluid,

μ_e [s⁴A² m^–3 kg^–1]-the electric permeability coefficient of biological fluid,

μ [1]-the dimensionless friction coefficient,

ξ_p[1]-the dimensionless pressures quotient,

ξ_C[1]-the dimensionless capacities quotient,

ξ_F[1]-the dimensionless friction forces quotient,

ξ_μ[1]-the dimensionless friction coefficients quotient,

ρ[kg/m³]-bio-fluid density,

ρ_e [C/m³ =As/m³]-electric space charge in the biological fluid,

σ_e [S/m]-electrical conductivity of phospholipids bilayer,

σ[1]-standard deviation symbol,

σ_s [1]-standard deviation for f_s,

σ_n [1]-standard deviation for f_n,

σ_N [1]-standard deviation for f_N,

ϑ₁[radian]-meridian coordinate in spherical surface,

φ[radian]-circumferential coordinate in spherical and cylindrical surface,

ψ[1]-dimensionless relative radial clearance,

ω[1/s]-angular velocity

Acknowledgments

Acknowledge the WSG Bydgoszcz University Garbary 2, Bydgoszcz, Poland,and Prof.dr hab.A.Gadomski from University of Science and Technology Bydgoszcz, dr hab.med.J. Cwanek from University Rzeszów, and dr.hab.eng.A.Miszczak from Maritime University Gdynia in Poland that have technically supported this work, excluding fund provider. Author thanks for many discussions during the work preparation.

Additional Materials 

The following additional materials are uploaded at the page of this paper.

1. Appendix 1

2. Appendix 2

3. Appendix 3

Author Contributions

Dr.hab.eng.Andrzej Miszczak contributes numerical calculations on the ground of analytical solutions, dr.hab.med.Janusz Cwanek contributes the medical measurements of bio surfaces, Prof.dr.hab.eng.Adam Gadomski enables scientific discussions and disseminates primary scientific results during the paper preparation.

Competing Interests

The author has declared that no competing interests exist