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

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):38; 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.

**Abstract**

This paper aims to highlight the result of a new progression of mathematical estimation methods of stochastic bio-hydrodynamic lubrication parameters for arbitrary, curvilinear, non-rotational, co-operating, living biological surfaces coated with phospholipid bi-layers. Movable, non-rotational, co-operating surfaces occur in various biological friction nods like the collar bone, the blade bone, the jump joint, and the wrist joint. Specifically, the author presents a synthetic and comprehensive estimation of stochastic bio-hydrodynamic lubrication parameters for co-operating, rotational cartilage bio-surfaces with phospholipid bi-layers occurring in human spherical hip joints and cylindrical elbow joints. The method of research discussed in this paper focuses on a review of stochastic analytical considerations performed by the author. This research is based on the measurements of the gap height between two movable, non-rotational bio-surfaces. The gap is restricted between two co-operating biological surfaces. After several experiments, it could be inferred that there are symmetric as well as asymmetric random increments and decrements in the gap height. Such changes are applicable to the hydrodynamic pressure, load-carrying capacity, friction forces, and wear of the co-operating biological surfaces in human friction nods and contacts. The prime purpose of this paper is to demonstrate the influence of variations in the expected values and standard deviation of the gap height on the hydrodynamic lubrication parameters that occur during the friction process. It can thus be concluded that the apparent dynamic viscosity of biological lubricant varies in the ultra-thin gap height direction, depending on the susceptibility of the superficial layer of the lubricated bio-surface. The results presented in this paper are obtained considering the 3D variations in the dynamic viscosity of the biological fluid, particularly the random variations crosswise the film thickness in non-Newtonian biological fluid properties.

**Keywords **

Non-Newtonian; bio-liquid; nano-layer; viscosity 3D-variations; stochastic solutions; theoretical; experimental results

**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

*η*in the gap is directly proportional to the surface coating of PL [

_{T}**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: *α _{1}, α_{2}, α_{3}*. Hence, the

*α*

_{1}-width direction of the non-rotational surfaces or the circumferential direction is accepted in the case of a rotational surface,

*α*

_{3}-longitudinal direction, and

*α*

_{2}-gap height direction.

The joint gap is filled with physiological bio-fluid lubricant. A characteristic constant dimensional gap height value *ε _{0}*[m] and the dimensionless total gap height function

*ε*based on variables

_{T1}*α*

_{1 }and

*α*

_{3}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

*ε*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

_{T1s }*δ*

_{1}. The expected value for the random variable of corrections

*δ*

_{1}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 *δ*_{1}.

The ordinates of the density function *f* denote the probabilities established for the random variable of corrections *δ*_{1 }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

*η*, 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 [

_{T}**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 *δ _{1}* in total gap height

*ε*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).

_{T}Figure 3 illustrates the measured positive and negative values of random variable gap height variations *δ _{1}*, 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*

*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,

**[V/m] denotes the electric intensity vector,**

*E***[A/m**

*J*^{2}] represents the electric current density,

**[T = kg/s**

*B*^{2}A] indicates the magnetic induction vector in bio-fluid,

**[m/s] is the biological fluid velocity,**

*v***[A/m] represents the magnetic intensity vector with components**

*H**(H*),

_{1}, H_{2}, H_{3}**[A/m] denotes the magnetization vector with components**

*N**(N*),

_{1}, N_{2}, N_{3}*σ*[S/m] is the electrical conductivity of phospholipids bilayer

_{e }*, μ*[mkgs

_{m }^{–2}A

^{–2}] signifies the magnetic permeability coefficient of biological fluid,

*μ*[s

_{e }^{4}A

^{2 }m

^{–3 }kg

^{–1}] is the electric permeability coefficient of the biological fluid,

*κ*[W/mK] denotes the thermal conductivity coefficient for the biological fluid,

*Φ*[W/m

_{F }^{3}] 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

^{3}= As/m

^{3}] is the electric space charge in the biological fluid,

*R*indicates the gas constant (8.3144598 J/K·mol),

_{g}*s*= (N

_{c}_{A}∙A)

^{–1}[mol/m

^{2}] is the concentration of phospholipid particles,

*γ*[J/m

^{2}= N/m] denotes the interfacial energy,

*γ*is the maximum interfacial energy of the lipid membrane such that 0.1 mN/m <

_{max}*γ*< 4 mN/m,

_{max }*K*[J] is the acid equilibrium constant (denotes the amount of energy needed to stretch the bilayer),

_{a }*K*[J] represents the base equilibrium constant (denotes the amount of energy needed to bend or flex the bilayer),

_{b }*α*[J] indicates the energy activity of protons,

_{H }*A*[m

^{2}] denotes the boundary regions between the areas of different concentrations of phospholipid molecules,

*N*= 6.024∙10

_{A}^{23 }[1/mol] represents the Avogadro number, and

*ρ*[kg/m

^{3}] 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

*2*

*T*

_{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:

*δ*. The

_{ij}, and Θ_{ij}*δ*denotes Kronecker Delta and

_{ij}*p*[Pa] is the pressure. For a non-Newtonian synovial fluid with Rivlin-Reiner type, the constitutive dependencies for apparent viscosity

*η*[Pas] take the following form [

_{T }**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*_{1}*(Θ)*[s^{–}^{1}], *I*_{2}*(Θ)*[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*[Pas

_{con }(n, pH, T, We)^{n}] denotes the fluid consistency coefficient

*, e*indicates tensor Levi-Civita, 50

_{ijk}^{0 }≤ We ≤ 80

^{0}represents the superficial layer wettability, and 2 ≤ pH ≤ 12 denotes the power hydrogen ion concentration. Geometrical non-linear relations between shear rate

*Θ*

_{i}_{j }[s

^{–1}] and biological fluid velocity component

*v*[m/s] are as follows [

_{i }**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

_{1}, α_{2}, α_{3}*EX*[

*p(α*

_{1},

*α*

_{3}

*)*], temperature

*EX*[

*T(α*], the velocity of synovial fluid

_{1},α_{2},α_{3})*EX*[

*v*], the viscosity of lubricating fluid

_{i}(α_{1},α_{2},α_{3})*EX*[

*η*], and gap height

_{T}(α_{1},α_{2},α_{3})*EX*[

*ε*] 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

_{T}(α_{1},α_{3})*γ*[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

*ε*[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

_{T }*α*is applied to the respective Lamé coefficients

_{1}, α_{2}, α_{3}*h*The aforesaid boundary simplifications of the arbitrary, curvilinear, non-rotational, non-parallel thin-layer surfaces follows

_{1}, h_{2}, h_{3}.*h*and

_{2}= 1*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

_{1}= h_{1}(α_{1},α_{3}), h_{3 }= h_{3}(α_{1},α_{3}).*α*direction and the non-monotone-generating line in

_{1}*α*direction, the Lamé coefficients follow

_{3}*h*or

_{2}= 1, h_{1}= h_{1}(α_{3}), h_{3}= h_{3}(α_{3}),*h*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 (

_{3}= 1*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 _{1}α*

_{1 }≤ 2a, –

*b*≤

*h*

_{3}α_{3 }≤

*b*, and 0 ≤

*α*

_{2 }≤

*EX(ε*. Symbol 2a [m] and 2

_{T})*b*[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

*η*of the biological liquid in the gap is directly proportional to the boundary surfaces

_{T }*A*[m

^{2}] between the areas of different phospholipid (PL) concentration having different superficial layer susceptibility of Symbol

*δ*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

_{v}*c*of collagen fibers in the bio-fluid. For

_{c}*δ*= 0.2, the value of

_{v}*c*is 100,000 mol/mm

_{c}^{3}, while for

*δ*= 0.6 concentration, the value of

_{v}*c*is 100 mol/mm

_{c}^{3}. The dimension of parameters

*L*and

_{a}*L*was determined by the formula (9b). Since 0 <

_{b }*L*< 1,

*lnL*is a negative number, while it is always the case that

*γ*>-

_{max }*kA*

^{–1}

*T*

*lnL,*where –

*100 < lnL < –10*.

The symbol 0.03 m/s <*v _{0 }*< 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 α

_{21}of the dimensionless functions

*v*and

_{11}*v*of velocity vector components

_{31}*v*[m/s] and

_{1}*v*[m/s] in directions

_{3}*α*in formula (9

_{1}and α_{3}*a*) are the result of dimensionless transformations, and calculations of constitutive dependences of apparent viscosity (

*3g*), where

*v*

_{1 }= v_{0}⋅*v*Newtonian fluid is used if the dimensionless flow index

_{11}, v_{3 }= v_{0 }⋅v_{31}.*n (0.8 < n < 1.2)*has the value of ‘1’. The following points have been assumed:

*δ*[m

_{E }^{2}/V

^{2}] is the experimental coefficient of the influence of electrical intensity

**and concentration of hydrogen ions**

*E**pH*in the bio-fluid. The value of coefficient

*δ*[m

_{E }^{2}/V

^{2}] 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

*δ*= 0.0003 m

_{E }^{2}/V

^{2}[

**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+δ*

_{E}E^{2 }=*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

^{2}), with the established values of

*We, δ*. A further increase in

_{v}, T, and v_{0}*pH*leads to a decrease in dynamic viscosity. The dynamic viscosity of the bio-fluid decreases with decreasing wettability

*We*at established values

*δ*. A drop in wettability from 70

_{v}, T, and v_{0}^{°}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}*,

*δ*, and

_{T}*δ*represent unknown random variables of pressure, temperature, and biological fluid velocity component corrections. If the random variable of gap height corrections

_{vi}*δ*

_{1 }= 0, then

*δ*=

_{p }*δ*=

_{T }*δ*= 0. In this case

_{vi }*EX(δ*0; therefore, under the assumption (10),

_{p}) = EX(δ_{T}) = EX(δ_{vi}) =*EX(p) = p*,

*EX(T) = T*, and

*EX(v*for i = 1, 2, 3. Thus, in this specific case, the system of equations (4)-(8) loses the stochastic properties.

_{i}) = v_{i}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}(*α*_{1},*α*_{2},*α*_{3})] [m/s] for i = 1, 2, 3; hydrodynamic pressure *EX *[*p*(*α*_{1},*α*_{3})] [Pa]; and temperature *EX *[*T*(*α*_{1},*α*_{2},*α*_{3})] [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*_{1 }is a function of *α*_{3}. 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*_{1 }is a constant value.

*3.2 Integrals for Random Hydrodynamic Lubrication of Non-Rotational Surfaces *

*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 (*α*_{1},*α*_{2},*α*_{3}). 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 _{1}(α*

_{1},

*α*

_{3}

*)*[

*m/s*] in

*α*

_{1}direction and

*U*

_{3}(α_{1},

*α*

_{3}

*)*[

*m/s*]

*in α*

_{3}direction, while the second bio-surface is stationary and limits a gap filled with a layer of bio-fluid with randomly varying height

*ε*. The following boundary conditions are applied to the functional components of the expected randomly varying viscous bio-fluid velocities EX(

_{T}*v*

_{1}), EX(

*v*

_{2}), and EX(

*v*

_{3}) in the directions

*α*

_{1},

*α*

_{2}, and

*α*

_{3}, 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 α_{1} 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^{4}/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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, 0 ≤ α_{2 }≤ EX(ε_{T}), EX(ε_{T}) = EX [ε_{T}(α_{ 1},α_{3})],_{ }h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], and b [m], U_{1 }= U_{1}(α_{1},α_{3}) [m/s]. Besides, A_{s }[α_{1},α_{2 }=EX(ε_{T}),α_{3}] = 1, and A_{η}[α_{1},α_{2 }=EX(ε_{T}),α_{3}] = 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 *α*_{3} 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_{3} = U_{3}(α_{1}, α_{3}) [m/s]. Moreover, M_{i} = ρ_{e}E_{i} [N/m^{3}] 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

^{6}/Ns

^{3}] 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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, 0 ≤ α_{2 }≤ EX(ε_{T}), EX(ε_{T}) = EX [ε_{T}(α_{ 1}, α_{3})],_{ }h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], b [m], U_{1} = U_{1}(α_{1},α_{3})[m/s], A_{ρ}_{1}[m^{12}/N^{3}s^{3}], A_{ρ}_{2}[m^{8}/N^{2}s^{2}], and A_{ρ}_{3}[m^{4}/Ns].

The influence of the centrifugal forces disappears when the Lamé coefficient *h _{1}* 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

_{3}*v*

_{2 }= 0 for

*α*

_{2 }= 0, leads to the following form of function of the expected random variable component of bio-fluid velocity in the direction α

_{2 }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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[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*_{2}) = 0 for *α*_{2 }= *EX*(*ε _{T}*), i.e., to the component of the bio-fluid velocity vector in the direction of gap height

*α*

_{2}. The new forms (19a, 19b, 19c) of the hydrodynamic pressure equations are derived after intense calculations.

For arbitrary, cooperating, biological, non-rotational surfaces in *α _{1} and α_{3}*, directions, where

*h*and

_{2 }= 1, h_{1 }= h_{1}(α_{1},α_{3}),*h*for a movable surface with linear velocity

_{3 }= h_{3}(α_{1},α_{3}),*U*in

_{1}*α*direction and

_{1}*U*in

_{3}*α*direction, the stochastically modified Reynolds type equation is derived. The bio-surfaces are lubricated with a thin liquid layer with dynamic viscosity

_{3}*η*(

_{T}*a*

_{1},

*α*

_{2},

*α*

_{3}) varying in three directions. This equation determining the expected function

*EX*[

*p*(

*α*

_{1},

*α*

_{3})] 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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], b [m], U_{1} = U_{1}(α_{1}, α_{3}), and U_{3} = U_{3}(α_{1}, α_{3}) [m/s]. The terms multiplied by the linear velocity components *U _{1} and U_{3}* represent the influence of surface motion on the hydrodynamic pressure. The integrals mentioned in equation (19

*a*) 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 *α _{1}* direction and non-monotone generating line in

*α*direction, where

_{3}*h*

_{2 }= 1, h_{1}*= h*and

_{1}(α_{3}),*h*for the motionless surface in

_{3 }= h_{3}(α_{3}),*α*direction, where,

_{3}*U*

_{1}(α_{3}) = ω⋅h_{1}, U_{3}*= 0*, and for the surfaces lubricated with a thin liquid layer with dynamic viscosity varying in three directions,

*η*the following stochastically modified Reynolds equation, which determine the expected function

_{T}(α_{1},α_{2},α_{3}),*EX*[

*p*(

*α*

_{1},

*α*

_{3})] 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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], b [m], and U_{1} = U_{1}(α_{3}).

For rotational, cooperating, biological surfaces with angular velocity ω in *α _{1}* direction and monotone motionless generating line in

*α*direction, where

_{3}*h*and

_{2 }= 1*h*,

_{1 }= h_{1}(α_{3}), U_{1 }= ω⋅h_{1}*U*and for the surfaces lubricated with a thin liquid layer with dynamic viscosity varying in three directions,

_{3 }= 0,*η*the

_{T}(α_{1},α_{ 2},α_{3}),*stochastically*modified Reynolds equation (

*19b*) is valid with

*h*

_{3 }= 1.For rotational, co-operating, biological surfaces with angular velocity *ω* in *α _{1}* direction and motionless in

*α*non-monotone direction, where

_{3 }*h*and

_{2 }= 1, h_{1 }= h_{1}(α_{3}),*h*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,

_{3 }= h_{3}(α_{3}), U_{1 }= U_{1}(α_{3}) = ω⋅h_{1}, U_{3 }= 0,*η*the following stochastically modified Reynolds equation, which determine the expected function

_{T}(α_{1},α_{3}),*EX*[

*p*(

*α*

_{1},

*α*

_{3})] 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_{1}α_{1 }< a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[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 (19*a*, *19b, 19c*), assuming that the value of atmospheric pressure *p _{A}* at the edges of the lubrication area to be

*Ω(α*

_{1},α_{3}):\[ 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 *α _{1} and α_{3}* directions, where

*h*and

_{2 }= 1*h*and

_{1 }= h_{1}(α_{1},α_{3}),*h*movable in

_{3 }= h_{3}(α_{1},α_{3}),*α*direction with linear velocity

_{1}*U*and in

_{1}*α*direction with linear velocity

_{3}*U*,

_{3}_{ }and for the surfaces lubricated with the liquid with dynamic viscosity varying in three directions:

*η*(

_{T}*α*

_{1},

*α*

_{2},

*α*

_{3}), the following differential equations are obtained after the transformations, which enable the determination of the expected function of temperature EX [

*T*(

*α*

_{ 1},

*α*

_{ 2},

*α*

_{3})] 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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, 0 ≤ α_{2 }≤ EX(ε_{T}), h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], b [m], U_{3} =* U _{3}(α_{1}, α_{3}),*

*U*and

_{1 }= U_{1}(α_{1},α_{3}),*M*[K/m

_{T}/κ=J^{2}/σ_{e}k^{2}]. The derivatives of the terms

*A*are presented in Appendix 2. The terms in (21), when multiplied by the linear velocity components

_{s}, A_{η}, and A_{p}*U*represent the influence of surface motion on temperature.

_{1}and U_{3},In order to determine the expected function of randomly variable temperature *EX *[*T(α _{1},α_{2},α_{3})*] 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*

_{0}ultimately lead to a constant temperature value

*f*on the first bio-surface (movable) and an unknown variable value of temperature changes

_{c}*f*(

_{p}*α*

_{1},

*α*

_{3}) 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}*(

*α*

_{1},

*α*

_{3}) on the surface of the acetabulum, the condition of transportation of heat flux density

*q*from the bonehead surface through the bio-liquid layer to the surface of acetabulum is used. This condition has the following form [

_{c}**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*

*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*[N] are determined in curvilinear directions

_{R3 }*α*

_{1}and

*α*

_{3}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,

*η*.Mentioned friction forces are formulated in the curvilinear orthogonal coordinates (

_{T}(α_{1},α_{2},α_{3)}*α*) in the following forms [

_{1},α_{2},α_{3}**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 _{1}α_{1 }≤*

*2*

*,*–

*b ≤*

*h*

_{3}α_{3 }≤*b, 0 ≤*

*α*

_{2 }≤*EX(ε*[m],

_{T}), h_{1}α_{1}*h*[m], a

_{3}α_{3}_{2}[m], a [m], b [m], and

*Ω*

*(α*lubrication surface,

_{1},α_{3})*h*and

_{2 }= 1, h_{1 }= h_{1}(α_{1},α_{3}),*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 [

_{3 }= h_{3}(α_{1},α_{3})**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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], b [m], U_{1 }= U_{1}(α_{1},α_{3}) [m/s], and U_{3 }= U_{3}(α_{1},α_{3}) [m/s].

The terms, when multiplied by the linear velocity components *U _{1}* and

*U*describe the influence of surface motion on the values of friction forces.

_{3},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 (17*a*, 17b), the *EX(F _{ρ}_{1}) *[m

^{9}/N

^{2}s

^{2}],

*EX(F*[m

_{ρ}_{2})^{5}/Ns]

*, EX(F*) [m] form of auxiliary expectancy function occurring in equation (25) and presenting the centrifugal forces [

_{ρ}_{3}**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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, 0 ≤ α_{2 }≤ EX(ε_{T}), h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], and b [m]. *EX *[*η _{T }(α_{1},α_{2},α_{3})*] denotes the expected function of the random dynamic viscosity of synovial fluid, which also varies in the direction of gap height.

*EX*[

*ε*(

_{T}*α*

_{1},

*α*

_{3})] 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 _{2 }= 1, h_{1 }= h_{1}(α_{1},α_{3}), h_{3 }= h_{3}(α_{1},α_{3}) *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_{1}α_{1 }≤ 2a, –b ≤ h_{3}α_{3 }≤ b, h_{1}α_{1}[m], h_{3}α_{3}[m], α_{2}[m], a [m], and b [m]. The load carrying capacity is along the direction vertical to the biological friction nod surface *Ω* [m^{2}] and in the direction opposite to load *W *[N].

For rotational, biological surfaces in *α _{1}* direction and non-monotone generating line in

*α*direction, where

_{3}*h*and

_{2 }= 1, h_{1 }= h_{1}(α_{1},α_{3}),*h*the expected value of load carrying capacity (27) in the curvilinear coordinates leads to the following double-integral form [

_{3 }= h_{3}(α_{1},α_{3}),**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 ≤ α_{1 }≤ α_{k }< 2π, –b ≤ h_{3}α_{3 }≤ b, EX(ε_{T}) = EX [ε_{T}(α_{1},α_{3})], where *α _{k}*

_{ }denotes the bio-fluid end coordinate in the circumferential direction α

_{1}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_{1}** and

*e*_{3}_{ }are the unit vectors tangential to the circumferential

*α*

_{1}and longitudinal

*α*directions.

_{3}**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 _{1}(α*

_{1},

*α*

_{3}

*)*[m/s] in

*α*

_{1}direction and

*U*

_{3}(α_{1},

*α*

_{3}

*)*[m/s] in

*α*

_{3}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 *δ _{1}*

_{ }of the gap height obtained experimentally, the expected correction value, i.e.,

*EX(δ*is determined from formula (2

_{1}) = m,*a*). 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) =

_{1}*ε*(

_{0}⋅ε_{T1s}*α*

_{1},

*α*

_{3}) =

*ε*(

_{T}*α*

_{1},

*α*

_{3},

*δ*0) represents the gap height without any random corrections.

_{1 }=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*(*δ*_{1 }=0) indicates the velocity without any random corrections of gap height.

The apparent viscosity variation is confirmed by formula (9*a*) 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 (9

*a*), 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^{2}, lnL = –50, γ = 2.5 mJ/m^{2},* where

*0.03 m/s < v*

_{0}*< 0.04 m/s; 0.2 < δ*and

_{v }< 0.6; T = 310 K;*k = 1.380649×10*. Hence,

^{–}^{23 }J/K*kA*.

^{–1 }T×lnL = –0.214 mN/mEventually, 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 (32*a*) 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 (32*a*) for the apparent viscosity of biological liquid proves that the decrements in the flow velocity v_{0} 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 _{1}h_{3} dα_{1}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(δ _{1 }= 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}*(

*δ*

_{1}= 0) indicates the friction force without random corrections of gap height, and

*α*,

_{F}*β*are the correction coefficients that are dependent on the stochastic changes in the fraction value

_{F}*ξ*.

_{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 *μ* (*δ*_{1 }= 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*

*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 *α _{1}* 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

_{3}*U*and

_{1 }= ω⋅h_{1}, U_{3 }= 0*α*Lubrication is performed only on the spherical surface region with the monotone generating line, i.e., inside the intervals:

_{1 }= φ, α_{2 }= r, α_{3}= ϑ.*0 < φ*

*< π, πR/6 < ϑ*

*< πR/2,*and

*ϑ = ϑ*meridian coordinates. Thus, the Lamé coefficients are as follows:

_{1}R*h*=

_{1 }*Rsinθ*= 1, and

_{1}, h_{3 }*R*[m] =

*L*which represents the radius of the sphere. Hence,

_{R}*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 *ε _{0}* 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) a*nd the center of the spherical surface of the acetabulum is at

*O*Eccentricity has the value of

_{1 }(x–Δε_{x}, y–Δε_{y}, z+Δε_{z}).*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

_{1}, 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,

*ε*-the least value of the gap height,

_{min}*C*-capacity,

*Ω*-load (Figures are based on the author`s studies).

*4.2 Example 2*

*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 *α _{1}* direction with angular velocity

*ω*and the upper surface (acetabulum) is motionless. Both the cylindrical surfaces are motionless in the longitudinal

*α*direction.

_{3}Hence, the values are *U _{1 }= ω⋅h_{1}, U_{3}* = 0 and

*α*and

_{1 }= φ, α_{2 = }r,*α*Lubrication is performed on the cylindrical surface region with the monotone generating line α

_{3}= z._{3}= z inside the intervals:

*0 < φ*

*< π/2, -b < z < +b,*and

*z*longitudinal coordinate. Thus, the Lamé coefficients are as follows:

*h*, where

_{1 }= R, h_{3}= 1*R*[m] =

*L*represents the radius of the cylindrical bone. Hence,

_{R}*U = ωR*and

*Ω = πRb*. Geometrical simulation of the elbow cylindrical joint is illustrated in Figure 5a.

*4.3 Example 3*

*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 *α _{1} = x_{1}, α_{2} = x_{2}, *and

*α*. The Lamé coefficients have the form:

_{3 }= x_{3}*h*and

_{2 }= 1*h*,

_{1 }= 1*h*

_{3 }=*1,*and

*0 < x*a

_{1 }<*, 0 < x*. 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

_{3 }< b ≡ L_{R}*α*direction with a linear velocity

_{3 }= x_{3}*U*and the upper bone (tight fitting dress) surface is motionless. Both the parallel surfaces are motionless in the meridian

_{3}*α*direction and, hence,

_{1 }= x_{1}*U*= 0 and

_{1 }*U = U*,

_{3}*Ω*= a

*b*. The two non-rotational parallel plane surfaces are highlighted in Figure 5b.

*4.4 Example 4*

*4.4 Example 4*

This case is valid for two movable, non-rotational and non-parallel surfaces represented in Figure 5c, where *h _{1}* =

*h*and

_{1}(α_{1},α_{3}),*h*

_{3}= h_{3}(α_{1},α_{3}).**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 *ε _{0}* 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

*δ*

_{1}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*of the correction parameters for

_{s}*(h), (e),*and

*(j)*surfaces occur less frequently than the anti-symmetrical functions (about 12 times in 100 measurements with a probability of

*P*= 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.

_{s }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*, where the probability of random decrease in gap height dominates over the probability of the random increase in gap height. Both

_{n}*f*and

_{N}*f*

_{n}_{ }are considered for surfaces

*(h), (e),*and

*(j).*In both cases of anti-symmetrical probability density distribution function

*f*each one of them occurs approximately 22 times out of 100 measurements with a probability of

_{N,}*P*= 0.22. In the next two cases of anti-symmetrical probability density distribution function

_{N }*f*, each one of them occurs 22 times in 100 measurements with a probability of

_{n}*P*= 0.22. The measurement multiplicity can be computed as 12 + 2 × 22 + 2 × 22 = 100, creating a probabilistically complete system of events [

_{n }**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*and

_{s}, f_{N},*f*of random changes in the gap height between hip joint spherical surfaces

_{n}*(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*obtained experimentally through gap height random changes between elbow joint cylindrical surfaces

_{n}*(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(δ*and

_{1}= 0), C(δ_{1}= 0), F_{R}(δ_{1}= 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

_{1}= 0)*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*= 0), based on formulas (27) and (29).

_{D}**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

*ξ*is calculated by assuming the stochastic changes for which the fraction value is the highest.

_{max}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 ε_{0}) 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 }*δ*

_{1}= 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

_{0 }*(h)*without random change.

**Figure 6 **Random lubrication parameters versus radial clearance 2 µm ≤ *ε _{0 }*≤ 10 µm between two spherical hip joint

*(h)*bio-surfaces, for stochastic changes in gap height 0.612·

*ε*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).

_{T}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

*ξ*= 1 obtained without random changes. Inside the standard deviation interval on the axis

_{C}*δ*, the expected value

_{1}*m*has been calculated for the corrected gap height

_{1}= 0.25*(m*, where the expected capacity attained the increase of about 55% (1.55) with a probability of

_{1}+1)⋅ε_{T }= 1.25ε_{T}*P = 0.916*for PC, compared to the capacity without random effects. The expected value

*m*for the corrected gap height

_{2}= 0.125*(m*is found, where the expected capacity attained increases of about 34% (1.34) with a probability of

_{2}+1)⋅ε_{T }= 1. 125ε_{T }*P = 0.937*for PS, compared to the capacity without random effects for

*ξ*= 1.

_{C}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 }*(1+m*, i.e., (0.74 < 1.51)×

_{2}-σ<1+m_{2}+σ) ε_{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

_{T}*P = 0.916*and gap height

*(m*for PC and with a constant probability of P = 0.937 and gap height

_{1}+1) ε_{T }= 1.25ε_{T}*(m*= 1.125

_{2}+1) ε_{T }*ε*for PS.

_{T}For the density function *f _{n}*, for PS, the expected value

*m*the corrected value of gap height =

_{2}= –0.125,*0,875ε*, and the standard deviation interval is (

_{T}*0.49 < 1.26)⋅ε*, while for PC, the expected value

_{T}*m*the corrected value of gap height =

_{1}= –0.25,*0.75ε*, and the standard deviation interval is (

_{T}*0.37<1.13)⋅ε*. Variations in the expected capacity values obtained for density functions

_{T}*f*and

_{N}*f*are the same.

_{n}Numerical calculations are based on the influence of phospholipid concentration on the spherical hip surface of articular cartilage (*A *= 10^{–15 }m^{2}); the electromagnetic field (*M _{i }*= 0,

*M*= 0) was omitted. For pressure and temperature, boundary conditions (20) and (22

_{T }*a*, b, c), and density functions of variable random joint gap height

*ε*described by the formulas (A3.1) and (A.3.2

_{T}*a*, 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 <

_{0 }*γ <*4 mN/m, the coefficient of collagen fibre concentration 0.2 <

*δ*< 0.6, the characteristic value of synovial fluid dynamic viscosity

_{v }*η*= 0.300 Pas, radius of the sphere of femoral head

_{0 }*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

_{0 }^{3}, angular velocity of bonehead ω = 5s

^{–}

^{1}, lubrication area Ω = 0.0004 m

^{2}, 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 *

*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 *

*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 [(9*a*) and (32*a*, 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*

*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*

*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

*δ*

_{1}= 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*

*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

*δ*

_{1}= 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*

*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

^{3}]-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

*α*directions for i = 1, 2, 3.,

_{i}*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

^{2}]-concentration of phospholipid particles,

*t *[s]*-*time,

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

*v _{i }*[m/s]-biological fluid velocity components in

*α*directions for i = 1, 2, 3.,

_{I}*v _{0}*[m/s]-characteristic value of bio-fluid velocity,

*x _{i }*[m]-rectangular coordinates for i = 1, 2, 3.,

*A *[m^{2}]-the region of areas of different phospholipids molecules concentration,

*A _{D }*[N]-adhesion force,

*A _{η}*[m

^{4}/Ns]-auxiliary function,

*A _{p }*[m

^{6}/Ns

^{3}]-auxiliary function,

*A _{ρ}_{1}*[m

^{12}/N

^{3}s

^{3}]-auxiliary function,

*A _{p2}*[m

^{8}/N

^{2}s

^{2}]-auxiliary function,

*A _{p3}*[m

^{6}/Ns]-auxiliary function,

*A _{s }*[

**1**]-dimensionless auxiliary function,

** B** [T =kg/s

^{2}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*),

_{1}, E_{2}, E_{3}*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*),

_{1}, H_{2}, H_{3}*I _{1}*(Θ)[s

^{–}

^{1}]-the first invariants of shear rate components,

*I _{2}*(Θ)[s

^{–}

^{2}]-the second invariants of shear rate components

*J *[A/m^{2}]-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

^{n}]

*-*fluid consistency coefficient,

*M _{i }*[N/m

^{3}]-electrical terms for i = 1, 3.,

** N** [A/m]-the magnetization vector,

*N _{A }*[1/mol]-Avogadro number (= 6.024∙10

^{23}),

*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

^{3}]-dissipation of energy,

*Ω*[m^{2}]-lubrication region,

*α*_{1}-width (circumferential) direction of the non-rotational (rotational) surfaces,

*α*_{3}-longitudinal direction,

*α*_{2}-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}[**1**]-dimensionless random variable corrections,

*δ _{ν}*[

**1**]-dimensionless coefficient determining the collagen fibers concentration,

*δ _{E }*[m

^{2}/V

^{2}]-coefficient describing influence of electrical intensity on dynamic viscosity,

*δ _{ij }*[

**1**]-dimensionless Kronecker Delta symbol,

*δ*[**1**]-unit tensor,

*ε _{0}*[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

^{–2}A

^{–2}]-the magnetic permeability coefficient of biological fluid,

*μ _{e }*[s

^{4}A

^{2 }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^{3}]-bio-fluid density,

*ρ*_{e }[C/m^{3 }=As/m^{3}]-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}ϑ_{1}[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 (if any)**

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

**References**

- Cwanek J. The usability of the surface geometry parameters for the evaluation of the artificial hip joint wear. Rzeszów: University Press; 2009. pp.10-70.
- Mow VC, Ratcliffe A, Woo S. Biomechanics of diarthrodial joints. Berlin-Heidelberg-New York: Springer Verlag; 1990. pp.30-80. [CrossRef]
- Wierzcholski K. Time depended human hip joint lubrication for periodic motion with stochastic asymmetric density function. Acta Bioeng Biomech. 2014; 16: 83-97.
- Andersen OS, Koeppe RE. Bilayer thickness and membrane protein function: An energetic perspective. Annu Rev Biophys Biomol Struct. 2007; 36: 107-130. [CrossRef]
- Bhushan B. Handbook of micro/nano tribology. 2nd ed. London and New York, Washington DC: CRC Press, Boca Raton; 1999. pp.150-200.
- Bhushan B. Nanotribology and nano-mechanics of MEMS/NEMS and BioMEMS/BioNEMS materials and devices. Microelectron Eng. 2007; 84: 387-412. [CrossRef]
- Chagnon G, Rebouah M, Favier D. Hyperelastic energy densities for soft biological tissues: A review.J Elast. 2015; 120: 129-160. [CrossRef]
- Gadomski A, Bełdowski P, Rubì JM, Urbaniak W, Augé II WK, Santamarìa-Holek I, et al. Some conceptual thoughts toward nanoscale oriented friction in a model of articular cartilage. Math Biosci. 2013; 244: 188-200. [CrossRef]
- Hills BA. Boundary lubrication in vivo. Proc Inst Mech Eng Part H J Eng Med. 2000; 214: 83-94. [CrossRef]
- Hills BA. Oligolamellar lubrication of joint by surface active phospholipid. J Rheumatol. 1989; 16: 82-91.
- Marra J, Israelachvili J. Direct measurements of forces between phosphatidylcholine and phosphatidylethanolamine bilayers in aqueous electrolyte solutions. Biochemistry. 1985; 24: 4608-4618. [CrossRef]
- Pawlak Z, Gadomski A, Sojka M, Urbaniak W, Beſdowski P. The amphoteric effect on friction between the bovine cartilage/cartilage surfaces under slightly sheared hydration lubrication mode. Colloids Surf B. 2016; 146: 452-458. [CrossRef]
- Pawlak Z, Urbaniak W, Hagner-Derengowska M, Hagner W. The probable explanation for the low friction of natural joints. Cell Biochem Biophys. 2015; 71: 1615-1621. [CrossRef]
- Pawlak Z, Figaszewski ZA, Gadomski A, Urbaniak W, Oloyede A. The ultra-low friction of the articular surface is pH-dependent and is built on a hydrophobic underlay including a hypothesis on joint lubrication mechanism. Tribol Int. 2010; 43: 1719-1725. [CrossRef]
- Pawlak Z, Urbaniak WI, Gadomski A, Yusuf KQ, Afara IO, Oloyede A. The role of lamellate phospholipid bilayers in lubrication of joints. Acta Bioeng Biomechs. 2012; 14: 101-106.
- Pawlak Z, Urbaniak W, Oloyede A. The relationship between friction and wettability in aqueous environment. Wear. 2011; 271: 1745-1749. [CrossRef]
- Pawlak Z, Petelska AD, Urbaniak W, Yusuf KQ, Oloyede A. Relationship between wettability and lubrication characteristics of the surfaces of contacting phospholipid-based membranes. Cell Biochem Biophys. 2013; 65: 335-345. [CrossRef]
- Ad P, Figaszewski ZA. Effect of pH on interfacial tension of bilayer lipid membrane. Biophys J. 2000; 78: 812-817. [CrossRef]
- Schwarz IM, Hills BA. Synovial surfactant: Lamellar bodies in type B synoviocytes and proteolipid in synovial fluid and the articular lining. Rheumatology. 1996; 35: 821-827. [CrossRef]
- Wierzcholski K, Miszczak A. Mathematical principles and methods of biological surface lubrication with phospholipids bilayers. Biosystems. 2019; 178: 32-40. [CrossRef]
- Wierzcholski K, Miszczak A. Electro-magneto-hydrodynamic lubrication. Open Phys. 2018; 16: 285-291. [CrossRef]
- Wierzcholski K. The joint cartilage lubrication with phospholipids bilayer. Tribologia. 2016; 2: 145-157. [CrossRef]
- Wierzcholski K. The topology of calculating pressure and friction coefficients for time-dependent human hip joint lubrication. Acta Bioeng Biomech. 2011; 13: 41-56.
- Yuan CQ, Peng Z, Yan XP, Zhou XC. Surface roughness evaluation in sliding wear process. Wear. 265: 341-348. [CrossRef]
- Fisz M. Probability calculation and mathematical statistics. In: Polish: Rachunek prawdopodobieństwa i statystyka matematyczna. Warsaw: PWN; 1967. pp.10-500.
- Helwig Z. Elements of probability calculations and mathematical statistics. In Polish: Elementy rachunku prawdopodobieństwa i statystyki matematycznej. Warsaw: PWN; 1977. pp.20-76.
- Truckenbrodt E. Fluidmechanik. In: German: Stroemungsmechanik. Berlin, New York: Springer Verlag; 1986.
- Schlichtling H. Boundary-Theorie. In: German: Grenzschicht-Theorie. Karlsruhe: Verlag G.Broun; 1989.
- Syrek P. Analysis of parameters of small dimensions space-applicators utilized in magneto-therapy-: Doctor thesis. In Polish-Analiza parametrów przestrzennych aplikatorów małogabarytowych, wykorzystywanych w magnetoterapii. Kraków: AGH University of Sciences and Technology; 2010.
- Fung VC. Biomechanice: Mechanical properties of living tissues. Berlin: Springer Verlag; 1993.
- Wierzcholski K. Nanotribology impact of run-walk, electromagnetic-hydrodynamic human joint and skin lubrication on the slimming and metabolic process. Ann Nanosci Nanotechnol. 2017; 1: 1-8.
- Wierzcholski K. The metabolic probabilistic effects of E-M Hydrodynamic lubrication for synovial fluid and sweat in the human joint and on the skin. In: Multiscale Locomotion. Bydgoszcz: Publishing House University of Science and Technology; Ed.A.Gadomski 2019. pp.81-119.