Fluid flow characteristics within an oscillating lower spherical surface and a stationary concentric upper surface for application to the artificial hip joint.

The synovial fluid motion in an artificial hip joint is important in understanding the thermo-fluids effects that can affect the reliability of the joint, although it is difficult to be studied theoretically, as the modelling involves the viscous fluid interacting with a moving surface. A new analytical solution has been derived for the maximum induced fluid motion within a spherical gap with an oscillating lower surface and a stationary upper surface, assuming one-dimensional incompressible laminar Newtonian flow with constant properties, and using the Navier-Stokes equation. The resulting time-dependent motion is analysed in terms of two dimensionless parameters R and β, which are functions of geometry, fluid properties and the oscillation rate. The model is then applied to the conditions of the synovial fluid enclosed in the artificial hip joint and it is found that the motion may be described by a simpler velocity variation, whereby laying the foundation to thermal studies in the joint.

Introduction

In the extension of the classic viscous Newtonian Couette flow in fluid mechanics, the analytical solution for the case of an oscillating plate is often mentioned. However, they are usually referring to the case of a single plate oscillating upon an infinite initially stationary mass of viscous fluid, the so-called Stokes' second problem [1, 2, 3]. For instance, the Stokes' second problem was solved using Laplace transform method [4]. Later on, the same problem was extended to include the slip condition at the oscillating wall using Laplace transform method [5]. An arbitrary oscillatory Couette flow problem was solved for rarefied gas in infinite medium [6].

The case of an oscillation plate on a finite channel gap with a stationary parallel plate has been solved using the complex velocity field [7, 8, 9]. Studies have also been made in the so-called large-amplitude oscillatory shear flow, for Newtonian and non-Newtonian flows [10, 11]. Recently, viscous dissipation effect was studied numerically for Newtonian fluid induced by an oscillation plate using the basic conservation equations [12]. For molecular flow studies, an oscillating shear-driven flow with slip flow modelling was analysed [13]. The earlier rarefied gas flow in infinite medium was extended into channel flow, using an analytical-numerical approach [14].

For non-parallel channel flow, a study on flow in microtubes was done experimentally and numerically using the cylindrical coordinate system [15]. However, the case of viscous flow in a spherical gap solved in the spherical coordinate system is seldom encountered. Hence, the motivation is to investigate this particular case, the effect of an oscillating lower spherical surface upon a stationary concentric upper spherical surface, and this is covered in Sections 2 and 3.

The second motivation is related to the study of the artificial hip joint. The advancement in hip prosthesis implant has provided a key remedy to the aging population affected by chronic hip pain, and to victims' disability due to serious accidents. The most common forms of chronic hip pain include osteoarthritis, rheumatoid arthritis, post-traumatic arthritis, avascular necrosis and childhood hip disease. According to the National Joint Registry, there were 796,636 total hip replacements activities performed in the United Kingdom, within the period April 2003 to December 2015 [16]. Recent data from American Academy of Orthopaedic Surgeons has shown that more than 300,000 total hip replacements are performed each year in the United States alone [17], and the trend is expected to increase approximately up to 170 % by 2030 [18]. Currently, total hip replacements have an approximate 15 years of service life, hence is not satisfactory for active patients under 60 years old. In fact, 44 % of current active patients would require 20–25 years of life expectancy for their total hip replacements [19, 20]. Hence, each and every part of the artificial hip joint should be enhanced to prolong the service life.

In the artificial hip joint, the synovial fluid is enclosed between the moving femoral head and the stationary acetabular cup. As the gap is small and the motion of the femoral head is somewhat oscillatory, the concentric spherical surfaces may be applicable as a first model to derive an analytical solution to describe the flow characteristics in the fluid. Hence the general solution so obtained in the first part will be applied to the conditions of fluid motion within the joint. A brief review on the hip prosthesis and conditions for appropriate application will be presented in Section 4.

Theory

Problem statement

The physical model considered is the pair of concentric spherical surfaces distanced W apart, defined as . The gap contains a fluid initially at rest, and fluid motion is induced in the fluid by the sinusoidal oscillation of the lower surface, while the upper surface is fixed. It is assumed that the magnitude of oscillation is small compared to the whole spherical surface, so that the induced motion is confined to part of the gap. Furthermore, analysis is performed on a cross-section passing through the origin of the spheres, with the induced fluid velocity in θ-direction being a function of the radius only. This implies we are analysing the maximum induced flow in the gap.

As shown in Fig. 1, the r−θ spherical coordinate system is taken with the origin at the centre of rotation of the lower surface, radii and refer to the lower and upper surface, respectively. Thus the motion of the lower plate can be specified as , where is the velocity component in the θ-direction, is the angular frequency, is the time and is the magnitude of the velocity.

Fig. 1

Physical model.

The general equation of motion for Newtonian fluids, the Navier-Strokes equation with typical symbols, for incompressible, laminar, constant property, is applicable [21]:

For the case of no pressure gradient and body force, with one dimensional flow, the velocity components along the r and φ directions are neglected, and Eq. (1) reduces towhere is the kinematic viscosity. The general solution to Eq. (2) will yield , subjected to three conditions, assuming non-slip conditions at the surfaces.

The solution

First, for convenience and generality, the set of Eqs. (2), (3), (4), and (5) is re-cast into dimensionless form by defining the following 5 parameters in Eq. (6).

Then the dimensionless set to be solved is

One method to obtain a solution of the set Eqs. (7), (8), (9), and (10), involving the time-dependent BC1, is to apply the Duhamel's theorem [22, 23] to the solution of the same problem statement, but replacing BC1 with a time independent condition: at the dimensionless velocity is unity. Duhamel's principle provides a solution to a linear, inhomogeneous partial differential equation [24]. This principle required a step input and then superposing with Duhamel's integral. The simpler time dependent problem has the solutionwhere in Eq. (11) denotes this specific solution as a step input (see AppendixA).

Then, according to Duhamel's theorem, the solution to the original problem (Eqs. (7), (8), (9), and (10)) iswhere is the running time variable, and .where the integration involving the exponential term in Eq. (13) is shown in Eq. (14).

Inserting into Eq. (12), we have

Hence, it can be shown that solution to the set Eqs. (7), (8), (9), and (10) is

Eq. (15) will be examined by way of graphical plotting in the next section.

Results & discussion

From Eq. (15), it is clear that the unsteady solution depends on the two dimensionless parameters, R and β, which contain information on geometry, fluid properties and the oscillation rate. It is then of interest to understand how the velocity profiles vary with R and β. Since and are always positive, R and must always be positive. increases when or increases and when decreases; the converse is also true. Hence in applying the results to a given problem, the relative values of and will be relevant. However, without referring to specific values of and , we first discuss the solution with respect to the two single variables, R and .

Velocity profiles with varying β

Figs. 2, 3, 4, 5, and 6 show six different velocity profiles for fluid flow with and 100 in six different radius ratio and 0.01. Since the dimensionless velocity is oscillating sinusoidally with time and the reference origin is stationary, the range for is clearly within -1 to 1. Seven whole numbers of dimensionless time are selected as to illustrate the change in velocity profiles over periods of time and 20. At the top surface, the fluid velocity is always zero for all β and , satisfying the no-slip boundary condition. Meanwhile, at the bottom surface, the fluid velocity follows the surface. Also, the fluid velocity is always zero at time zero reflecting the assumed initial condition.

Fig. 2

Fluid velocity profile at varying dimensionless parameter when R = 0.99, for (a) (b) (c) (d) (e) (f) .

Fig. 3

Fluid velocity profile at varying dimensionless parameter when R = 0.90, for (a) (b) (c) (d) (e) (f) .

Fig. 4

Fluid velocity profile at varying dimensionless parameter when R = 0.50, for (a) (b) (c) (d) (e) (f) .

Fig. 5

Fluid velocity profile at varying dimensionless parameter when R = 0.10, for (a) (b) (c) (d) (e) (f) .

Fig. 6

Fluid velocity profile at varying dimensionless parameter when R = 0.01, for (a) (b) (c) (d) (e) (f) .

It is within the plates that the velocity profile depends on the values of R and β. The fluid motion is moving nearly with a linear velocity profile close to the bottom moving surface when β equal to 1 and R near to value of 1 as shown in Figs. 2(a) and 3(a). This is because the values of the terms under the summation sign in Eq. (15) are small enough not to influence the fluid velocity profile at all time. Also, the fluid velocities at different locations are varying due to the effect of the fluid inertia that influences the acceleration of the fluids and give rise to a velocity phase lag as time moves on. This effect is buried inside the general solution. For example, in Fig. 2(c), the fluid near to the bottom surface is having , while at , the .

The thickness of the stagnant layer from the top surface increases at two conditions, either the value of R is small or β is large. As observed from Figs. 2, 3, 4, 5, and 6, for β = 100, more than half of the fluid from the top surface remain stagnant for all R conditions and it is getting larger when R is going to a smaller value. Another observation from Figs. 2, 3, 4, 5, and 6 is that for and , the fluid still have a significant motion within the spherical gap. At the extreme case where , most of the fluid remains stagnant for all β conditions. Specifically, the motion only varies between .

Effects of the terms under the summation sign

From the general solutions, Eq. (15), it is seen that the fluid velocity compromises a sine term minus all the terms under the summation sign. The latter is named here as , as in Eq. (16). Then we have

Figs. 7, 8, 9, 10, and 11 show the effects of the terms under the summation sign for dimensionless parameters of and 100 while and 0.01, with selected varies of and 20. The first figures, labelled as (a) in Figs. 7, 8, 9, 10, and 11, show the variation of with directional sign, while the rest of the figures show only the magnitude, .

Fig. 7

The behaviour of at selected β at t* = 1 and with selected t* at R = 0.99, for (a) (b) (c) (d) (e) (f) .

Fig. 8

The behaviour of at selected β at t* = 1 and with selected t* at R = 0.90, for (a) (b) (c) (d) (e) (f) .

Fig. 9

The behaviour of at selected β at t* = 1 and with selected t* at R = 0.50, for (a) (b) (c) (d) (e) (f) .

Fig. 10

The behaviour of at selected β at t* = 1 and with selected t* at R = 0.10, for (a) (b) (c) (d) (e) (f) .

Fig. 11

The behaviour of at selected β at t* = 1 and with selected t* at R = 0.01, for (a) (b) (c) (d) (e) (f) .

To gauge the effect of R on , at constant β, it is seen that the peak of reduces significantly with the R value, at all dimensionless time conditions. For example, in Figs. 7(f) and 11(f), at , the peak of at and at . And when β is changing at constant R, the peak of the changes with both the β value and . Although the peaks for any particular value of β seem to vary with , however the maximum peak does seem to depend on the value of β. For example, the maximum peak found for , is always higher than the maximum peak for , and the maximum peak for is always higher than the maximum peak for , etc. Another observation is that at , the variation of the peak of can be considered to have negligible effect to the terms under the summation sign for all time.

There is only one peak of for all the selected at varying R and. However, there are two peaks occurring for some particular time at and 100, which is because there are actually two turnings for the fluid motion for at , as shown in Figs. 2(e) and 7(a). Moreover, the peak location of is varying within the spherical gap with time.

It can be concluded that the effect of would not vary much when the . On the other hand, when either R or β increases, the fluctuations of the become significant with time. Among the selected R and β, the magnitude of shows the largest difference when and 100.

Example

Physical characteristics of an artificial hip joint

The main components of an artificial hip joint include an acetabular component and the head of femur as shown in Fig. 12. There are several different combinations of materials for the components, including metal, plastic and ceramic materials [19]. In the 1960s, metal-on-plastic (MOP), also known as soft on hard couple is the most common. Later, it was replaced by hard on hard materials, called the metal-on-metal (MOM) combination introduced around 1997 for their lower wear rate [25, 26].

Fig. 12

Main components in the structure of an artificial hip joint [29].

The synovial fluid acts as a biological lubricant or as a shock absorber by filling the gap between the acetabular cup and the femur head [27]. This gap usually ranges from , particularly in MOM combination [28]. Although synovial fluid generally reduces the wear rate, they would at times contribute to the wear rate too. This section relates the derived maximum induced velocity profile applied to synovial fluids within the gap of an artificial hip joint.

Modelling the synovial fluid and the artificial hip joint motion

The synovial fluid in the artificial hip joint shows non-Newtonian behaviours at very low and high shear rate, described as rheopectic and shear-thinning effects, respectively [30, 31]. However, most of the synovial fluid cases can be treated as a simple incompressible fluid, Newtonian and isoviscous fluid model [32, 33], and hence the justification of applying the model as presented in Section 2.1.

The actual shape for artificial hip joint is hemi-spherical as shown in Fig. 12 and the parts of joint can be modelled in spherical geometry. Thus, in Fig. 1, the upper surface represents the acetabular cup and the lower surface refers to the head of femur, with the gap distance W.

There are some typical parameters that will be used to evaluate the fluid motion in the artificial hip joint as summarised in Table 1.

Table 1

Typical parameters used in deriving the velocity [34, 35].

Parameter	Value
Dynamic viscosity of synovial fluid, μ	10 mPa⋅s
Gap of the joint, W	30–150 μm
Angular frequency of bottom moving surface, ω	1 rad/s
Density of synovial fluid, ρ	1650 kg/m^3

Effects of the terms under the summation sign for artificial hip joint motion

As pointed out in Section 3, for the joint application, appropriate realistic values need to be considered, and inserted into Eqs. (15) and (16).

Fig. 13 shows the variation of at the assigned values according to the hip joint geometry and synovial fluid properties. The value selected for is , based on the average radius for the head of femur [28]. In Fig. 13(a) and (b) are shown the gap with and respectively. The angular frequency , which represents the average motion for the walking condition for an artificial hip joint. Hence, the value of R and β varies accordingly.

Fig. 13

The variation of with respect to and at different assigned values to and . (a) (b) .

Both the graphs in Fig. 13 show the occurrence of largest magnitude for synovial fluid motion located between to 0.5 and it is consistent throughout the selected t*. Thus, next section will reveal this range of r* for a longer period of time for both defined conditions where and .

The largest differences for Fig. 13(a) and (b) occurr at , and in term of magnitude are and , respectively. The collected data shows that the largest differences in the terms under the summation sign is insignificant and the contribution of these terms in Eq. (16) has negligible effect.

Fig. 14 illustrates the peak values of at three selected r*. The selected three r* values are 0.3, 0.4 and 0.5. An average value of for the three selected r* is plotted as well for both conditions where and , respectively.

Fig. 14

The peak values of for three selected and an average of r* at different assigned values to and . (a) (b) .

Recall that in Eq. (16), the consists of three types of time dependant terms, those that have the cosine term, followed by the sine term and lastly the exponential term. The exponential time dependant term vanishes as time increases. As observed from Fig. 14(a) and (b), the first peak cycle is lower compared to the later cycles. From the second cycle, the peaks of the remain the same with time, because the effect of the exponential term causes the first peak to be lower and the effect vanishes from the second cycle onwards, leaving only the effects of the cosine and sine terms in .

From Figs. 13 and 14, the approximate peaks of for R = 0.9894 and 0.9979 have magnitudes 0.00025 and 0.000010, respectively, over the range of fluid velocity , from 0 to 1. Thus, in terms of percentage of the maximum effect for both the minimum and maximum gap, are 0.001 and 0.025 %, respectively. It can be concluded that, for the maximum gap of artificial hip joint, the term, as given in Eq. (16), can be neglected when determining the velocity of the synovial fluid. This is because the maximum effect of to the overall is insignificant.

The velocity profile in the synovial fluid

Since the effect of terms under the summation sign is negligible, the circumstance reduces Eq. (15) to the simple Eq. (17),

Fig. 15 reflects the velocity profile for synovial fluid at different times. As observed from Fig. 15, the fluid velocity profile behaves proportionally with different gradients towards the position of the moving surface. The line gradient changes with time and the range of velocity varies to-and-fro between -1 to 1, reflecting the stationary cup and the sinusoidally movement of the femur head. Another observation from Fig. 15 is that, there are no significant differences in terms of fluid motion between the range of the different gap distance at the artificial hip joint.

Fig. 15

The change in velocity profile with respect to for six different times at different assigned values to and . (a) (b) .

The fluid motion for the range of gaps for an artificial hip joint have insignificant difference with the unsteady variation of fluid velocity being the same. Fig. 16 shows the time variation of fluid velocity in 0.2 steps, from to 1.0. The approximate maximum magnitudes for fluids velocity at 0, 0.2, 0.4, 0.6, and 1.0 are = 1.0, 0.8, 0.6, 0.4, 0.2 and 0, respectively. The maximum magnitude of the fluid velocity is decreasing as the position increasing and the fluid velocity varies sinusoidally henceforth.

Fig. 16

Unsteady variation of fluid velocity at 6 different positions in the gap.

Conclusion

This study reports a derived analytical expression for the maximum induced velocity in a spherical gap with the lower surface oscillating sinusoidally and upper surface remaining stationary. The behaviour of fluid motion is first analysed with respect to dimensionless property parameters at dimensionless time and dimensionless location . Results show that the fluid motion is dominated by a sinusoidal term, plus terms under the summation sign that have both sinusoidal and exponential decay characteristics. Single or double velocity peaks may occur in the channel and there can be a lag with respect to the activity moving surface.

When applied to realistic properties of the synovial fluid in a hip joint, it has been shown that the terms under the summation sign can be negligible. This is when or 0.01. It is of importance to further study into the thermal effect of the motion of synovial fluids in artificial hip joints.

Declarations

Author contribution statement

C. P. Tso: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

C. H. Hor: Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

G. M. Chen: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data.

C. K. Kok: Conceived and designed the analysis.

Funding statement

This work was supported by the Ministry of Higher Education, Malaysia, (FRGS/1/2016/TK03/MMU/01/1).

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.