人工髖關節球形凹坑微織構表面摩擦學性能研究*

楊小洋 胡亞輝 毛璐璐 張春秋 鄭清春

(1.天津理工大學天津市先進機電系統設計與智能控制重點實驗室 天津 300384；2.機電工程國家級實驗教學示范中心(天津理工大學) 天津 300384)

隨著我國關節骨病發病率的增加，人工關節置換手術的臨床運用越來越廣泛。但相比于天然關節，由于一些主要結構的缺失，使得人工關節在減摩潤滑性能上與天然關節存在一定差距[1-3]。因此研究改善人工關節的摩擦性能具有十分重要的意義。

大量的研究發現，適當的表面織構可以有效地改善摩擦表面間的摩擦性能[4-7]。WANG等[8]通過在SiC試樣表面加工一系列直徑和深度的織構并進行摩擦學實驗發現，與無織構試樣相比，不同參數的微織構均起到了一定程度的減摩作用。汪家道等[9]利用激光刻蝕在試樣表面加工不同深度的微織構，通過對比其表面油膜壓力發現，選擇合適的微織構參數可以提高試樣表面的油膜壓力。 目前研究已經證明了可以通過加工微織構來改善人工髖關節的摩擦磨損性能。毛璐璐等[10]為改善鈦合金人工髖關節假體的耐磨性，通過實驗研究發現，加工合適的微織構參數可以有效降低鈦合金的摩擦因數。奚鐸聞等[11]利用數值方法研究球面紋理對人工髖關節潤濕性能的影響，結果表明相對于光滑表面，帶有微凸起的表面可以提升人工髖關節的承載能力。

MURAKAMI等[12]通過對橢圓形關節軟骨標本進行實驗發現，人在行走和跑步時關節處于流體潤滑狀態。由于人工關節在使用中大部分磨損發生在行走和跑步狀態下。因此本文作者提出一種在人工髖關節表面加工球形凹坑微織構的方法來改善關節表面摩擦性能，并利用CFD數值仿真方法，對流體潤滑狀態下表面微織構參數對摩擦學性能的影響進行研究。

1 人工髖關節表面微織構CFD模型建立

1.1 微織構流域單元幾何模型建立及參數設定

人工髖關節摩擦副表面間潤滑油膜厚度為微米級，髖關節球頭表面曲率半徑遠大于膜厚，故在此忽略摩擦副表面曲率的影響，將人工髖關節摩擦副表面簡化成平面進行研究。簡化后的微織構表面，如圖1所示。將微織構單元所在區域上下摩擦副之間的流體區域定義為流域，取其中一個微織構流域單元作為研究對象，如圖2所示。利用三維建模軟件建立微織構流域單元的幾何模型，微織構流域單元的長為a，寬為b，微織構沿速度方向的間距為c，織構直徑為D，摩擦副上下表面間初始油膜厚度為h0，球形凹坑微織構深度為hp，織構最大深度為l。其中，流域單元長度為

圖1 簡化后微織構表面示意Fig 1 Schematic of micro texture surface after simplification

圖2 微織構流域單元幾何模型Fig 2 Geometric model of micro textured watershed unit(a) parameter setting of micro texture element；(b) schematic of film thickness

a=D+c

(1)

此外定義織構面積密度S為控制單元內微織構圓形開口面積占微織構單元面積的比值，微織構深徑比hk為織構最大深度與織構直徑的比值，其公式為

S=πD2/(4ab)

(2)

hk=l/D

(3)

1.2 基于Reynolds方程的潤滑控制方程

建立人工關節摩擦副間潤滑控制方程前做以下假設：

(1)人工關節摩擦副間被潤滑液均勻分隔開，即處于全膜潤滑狀態；

(2)潤滑液為不可壓縮牛頓流體，采用層流方式；

(3)忽略流體的慣性力，忽略摩擦副間隙油膜兩端壓力差；

(4)摩擦副表面絕對光滑，流體與壁面接觸無相對滑動。

得到簡化后的Reynoleds方程[13]如下：

(4)

式中:h為潤滑油膜厚度；η為潤滑油黏度；U為滑動速度；p為油膜壓力。

1.3 織構流域單元膜厚方程

取織構上表面為坐標平面，如圖2所示，則摩擦副之間潤滑液膜厚方程為

(5)

式中:h0為人工髖關節摩擦副上下表面間初始油膜厚度；Φ為微織構區域;hp為球形凹坑微織構深度，

(6)

1.4 微織構流域單元邊界條件及網格劃分

利用ICEM CFD進行前處理，首先將各邊界單獨建立為Part，在ANSYS Fluent中對相應的塊設置邊界條件，如圖3所示。將垂直于速度U的入口、出口邊界設置為周期性邊界條件；將平行于速度U的前、后邊界設置為對稱邊界條件；控制單元的上、下邊界為無滑移壁面邊界條件，其中下邊界為固定壁面，上邊界以速度U運動。

圖3 微織構流域單元邊界條件設置Fig 3 Boundary condition setting of micro texture watershed unit

之后利用ICEM CFD進行網格劃分，采用有限體積法將計算流體域控制單元離散成為六面體結構性網格，網格將計算區域分割成35 000個足夠小的計算區域，然后在每一個計算區域內應用流體控制方程求解，從而獲得計算區域內的物理量分布。為保證計算的準確性，在計算前需檢查網格質量。采用不可壓縮牛頓流體作為連續相介質。壓力-速度耦合格式設為SMIPLE；壓力項選用PRESTO!離散，動量項采用二階迎風格式離散，以保證求解過程的精確和穩定。

1.5 微織構流域單元CFD仿真結果

當摩擦副上表面以一定速度相對運動，由于潤滑液黏性作用導致微織構表面間液體被帶動，在摩擦副間形成楔形間隙，產生流體動壓效應，同時由于潤滑液黏性導致潤滑液中剪切效應的產生。圖4(a)、(b)所示分別為面積密度為15%、微織構深度為10 μm、速度為0.25 m/s條件下，球形凹坑微織構流域單元的壓力分布和剪應力分布。

圖4(a)所示為球形凹坑微織構單元壓力分布云圖。沿速度方向從左到右，在凹坑織構區域內，左端流入處由于上下表面間距增大形成發散楔導致壓力下降產生負壓，右端流出處由于上下表面間距減小形成收斂楔導致壓力上升產生正壓，正壓最高值大于負壓最低值。由于織構凹坑兩側正負壓力差的存在使得流體流經織構摩擦副表面時形成動壓效應，這對增加油膜平均承載力與減小摩擦因數從而改善摩擦副之間的摩擦性能起到積極作用[12]。

圖4 球形凹坑微織構流域單元仿真結果Fig 4 Simulation results of watershed unit with micro texture ofspherical pit (a) pressure distribution(S=15%,l=10 μm,U=0.25 m/s)；(b) shear stress distribution(S=15%,l=10 μm,U=0.25 m/s)

1.6 CFD仿真實驗驗證

為驗證仿真模擬結果的可靠性，將膜厚h0為6 μm的無織構流域模型的剪應力理論值與模擬值進行比較。無織構壁面剪應力的理論公式[14]如下：

(7)

式中：τ為理論公式計算的剪應力；du/dz為速度梯度。

式中速度梯度等于速度U與膜厚h0之比。由表1可知，在不同速度下模擬值與理論值相對誤差約為1%，驗證模擬方法的正確性。

表1 理論剪應力與模擬剪應力及誤差對比Table 1 Comparison of theoretical shear stress and simulated shear stress and error

2 CFD仿真實驗設計及數據處理

2.1 微織構流域單元CFD仿真實驗設計

對單個微織構控制單元進行數值模擬計算，在此設置凹坑織構直徑D與單元寬度b為固定值，分別為100、200 μm，油膜初始厚度h0取6 μm。各組微織構流域單元數據如表2所示。

表2 微織構流域單元參數設置Table 2 Parameter setting of micro texture watershed unit

潤滑液黏度η設為0.05 Pa·s。在數值分析中將流體域單元上邊界的滑動速度U分別設置為0.25、0.3、0.35、0.4 m/s，來研究不同速度下微織構參數對摩擦學性能的影響。

2.2 仿真結果數據處理

利用ANSYS Fluent分別對各微織構單元在不同速度條件下進行求解，利用Fluent自帶后處理功能在織構控制單元內將壓力p對控制單元下壁面進行積分得到油膜沿z方向的承載力F，以及微織構表面剪應力τ，之后通過計算得到單位面積內平均承載力Fa和微織構表面沿x方向摩擦力Ff，公式為

(8)

(9)

Fa=F/(ab)

(10)

將摩擦因數f定義為摩擦力與承載力之比，公式為

f=Ff/Fa

(11)

以平均承載力Fa和摩擦因數f作為衡量潤滑性能的指標來分析微織構參數對摩擦學性能的影響。

3 結果與討論

3.1 深徑比hk對微織構表面摩擦學性能的影響

為討論深徑比hk對平均承載力Fa以及摩擦因數f的影響，選取微織構面積密度為10%的情況下，不同深徑比的微織構流域單元作為研究對象。

不同速度下微織構深徑比與平均承載力的關系，如圖5所示。在深徑比由0.06升高至0.12的過程中，各速度條件下的平均承載力均呈現出先下降后上升之后再下降的趨勢。在所設速度范圍內深徑比為0.06的微織構表面平均承載力最高，深徑比為0.12時的平均承載力最低。

圖5 不同速度下深徑比與平均承載力關系Fig 5 Relationship between depth diameter ratio andaverage bearing capacity at different speeds

圖6所示為不同速度下摩擦因數與深徑比的關系。可知，在同一深徑比條件下織構表面摩擦因數隨速度增加而降低。然而在相同速度下，在所設計的深徑比范圍內，摩擦因數隨深徑比增加呈現出先上升后下降再上升的波動趨勢，其中深徑比為0.06時微織構表面摩擦因數最低，當深徑比達到0.12時微織構表面摩擦因數最高。

圖6 不同速度下摩擦因數隨深徑比的變化Fig 6 Variation of friction coefficient with depthdiameter ratio at different speeds

一方面，隨著微織構深徑比的增加，由于流體剪切效應在微織構凹坑中產生的微渦流逐漸增強削弱了流體動壓，從而導致平均承載力下降；另一方面，隨著微織構深徑比的增加，微織構中儲存的潤滑液也隨之增加，可以為摩擦副提供更充足的潤滑，這可能是導致深徑比為0.10的表面織構比深徑比為0.08的表面織構的摩擦因數低且產生的平均承載力高的原因。當深徑比進一步增加時微織構中微渦流削弱作用的影響更加顯著，這可能是導致深徑比從0.10繼續提高到0.12時摩擦因數上升且平均承載力下降的原因。

3.2 面積密度S對微織構表面摩擦學性能的影響

為討論面積密度S對平均承載力Fa的影響，選取不同面積密度下深徑比為0.10的微織構流域單元作為研究對象。

圖7示出了在一定相對滑動速度下，在人工關節摩擦副表面設置球形凹坑微織構對摩擦副間油膜承載力的影響。可見，在所設速度范圍內，平均承載力隨著相對滑動速度的增加而增加；在相同速度條件下，隨著微織構面積密度由5%升至30%，平均承載力呈現出先升高后降低的趨勢，并在面積密度為10%時取得最高值。

圖7 不同速度下面積密度與平均承載力關系Fig 7 Relationship between area density and averagebearing capacity at different speeds

圖8反映出不同相對滑動速度下，摩擦因數與織構面積密度的關系。織構面積密度由5%增加到30%的過程中，摩擦因數呈現出先逐漸降低后緩慢上升的趨勢，并在織構面積密度為25%時摩擦因數取得最低值，這與吳霄[13]所得模型及實驗規律基本一致。

圖8 不同速度下摩擦因數隨面積密度的變化Fig 8 Variation of friction coefficient with areadensity at different speeds

隨著織構面積密度的增加，流體動壓潤滑效應逐漸增強，導致平均承載力提高以及摩擦因數逐漸降低。當面積密度到達一定程度且繼續增加時，織構沿相對滑動速度方向分布過于密集，導致難以建立穩定的潤滑油膜[15]，這導致油膜承載力的下降以及摩擦因數的上升，且摩擦因數與面積密度的這一趨勢隨相對滑動速度的降低逐漸明顯。對于平均承載力的最高值與摩擦因數的最低值對應的面積密度不一致這一問題，作者認為由于影響摩擦因數的因素有多種，平均承載力只是其中之一，由于面積密度的改變導致其他因素的變化，這可能是導致這一問題出現的原因。

4 結論

(1)在人工髖關節表面設置球形凹坑微織構可以提高油膜平均承載力降低摩擦因數，從而起到減小關節的摩擦磨損提高人工關節使用壽命的作用。

(2)摩擦副相對滑動速度對油膜平均承載力以及摩擦因數有明顯影響，在所設速度范圍內油膜平均承載力隨速度提高而增加，摩擦因數隨速度提高而減小。

(3)將深徑比hk作為研究織構深度的參考值，研究發現在所設計深徑比范圍內，平均承載力先下降再上升最后再下降，摩擦因數呈現與平均承載力相反的波動趨勢。當深徑比為0.06時平均承載力最高，摩擦因數最低。

(4)在所設計面積密度范圍內，平均承載力先上升再下降，摩擦因數先下降再上升。當面積密度為10%時平均承載力最高，面積密度為25%時摩擦因數最低。

(5)利用CFD數值仿真方法研究微織構參數對人工骨關節摩擦學性能影響的方法是可行的，為今后進一步開展相關研究奠定了一定基礎。