盾構(gòu)機(jī)唇形密封流固熱耦合仿真研究*

紀(jì)佳馨 蘇世東 項(xiàng) 沖 郭 飛

（1. 中國(guó)石油大學(xué)（華東） 機(jī)電工程學(xué)院 山東青島 266580；2. 清華大學(xué)高端裝備界面科學(xué)與技術(shù)全國(guó)重點(diǎn)實(shí)驗(yàn)室 北京 100084）

隨著國(guó)家城市化進(jìn)程的發(fā)展， 盾構(gòu)掘進(jìn)技術(shù)在城 市地鐵、 鐵路交通等領(lǐng)域得到廣泛應(yīng)用。 鑒于盾構(gòu)機(jī)工作過程中特殊的隧道環(huán)境， 主驅(qū)動(dòng)唇形密封系統(tǒng)[1-3］一旦失效， 會(huì)導(dǎo)致外界環(huán)境中的土砂倒灌進(jìn)而破壞主軸承， 嚴(yán)重影響工程進(jìn)度。 唇形密封件的性能直接決定了整臺(tái)盾構(gòu)機(jī)的施工效率， 且在施工過程中應(yīng)避免唇形密封件的更換， 為免影響施工效率， 其故在盾構(gòu)機(jī)主驅(qū)動(dòng)密封設(shè)計(jì)階段， 預(yù)測(cè)唇形密封性能是否滿足要求尤為重要。

國(guó)產(chǎn)盾構(gòu)機(jī)主驅(qū)動(dòng)密封件多依賴國(guó)外進(jìn)口， 且價(jià)格高昂， 為提高國(guó)產(chǎn)大型盾構(gòu)機(jī)的國(guó)際競(jìng)爭(zhēng)力， 亟待開展盾構(gòu)機(jī)主驅(qū)動(dòng)密封研究。 譚鋒等人[4］利用正交試驗(yàn)法針對(duì)盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封開展結(jié)構(gòu)優(yōu)化研究，根據(jù)密封工作要求確定最優(yōu)密封結(jié)構(gòu)， 優(yōu)化后密封件接觸寬度降低23%。 張中華等[5］在主驅(qū)動(dòng)密封失效分析的基礎(chǔ)上， 總結(jié)提出密封安裝結(jié)構(gòu)和工藝優(yōu)化方案， 改善盾構(gòu)機(jī)主驅(qū)動(dòng)密封性能。

目前， 針對(duì)盾構(gòu)機(jī)主驅(qū)動(dòng)密封的研究集中于結(jié)構(gòu)和安裝優(yōu)化， 尚無預(yù)測(cè)密封性能的研究手段。 盾構(gòu)機(jī)主驅(qū)動(dòng)唇封密封介質(zhì)為潤(rùn)滑脂， 工作時(shí)唇口溫度可達(dá)50～60 ℃， 為更好地預(yù)測(cè)唇封的密封性能， 本文作者考慮潤(rùn)滑脂流變特性、 唇口溫度對(duì)流場(chǎng)分析、 密封材料的影響， 建立盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封流固熱耦合仿真模型； 然后利用流速分離法推導(dǎo)潤(rùn)滑脂二維雷諾方程， 采用赫茲接觸模型計(jì)算粗糙峰接觸壓力， 結(jié)合有限元軟件開展熱力耦合分析， 實(shí)現(xiàn)唇封溫度場(chǎng)及摩擦力矩、 泄漏率等關(guān)鍵性能參數(shù)的定量預(yù)測(cè)。 該模型對(duì)盾構(gòu)機(jī)主驅(qū)動(dòng)密封設(shè)計(jì)研究具有重要意義。

1 數(shù)值仿真模型

唇封系統(tǒng)主要由密封唇、 旋轉(zhuǎn)軸和密封介質(zhì)三部分組成， 其相互作用共同決定唇封系統(tǒng)的密封性能。如圖1 所示為唇封與軸接觸區(qū)的微觀示意圖， 可以看出建立盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封數(shù)值仿真模型， 需要綜合分析密封過程涉及的物理要素， 包括密封介質(zhì)的流體力學(xué)分析， 微觀粗糙峰與軸之間的接觸力學(xué)分析和密封唇的固體力學(xué)分析。 圖中所示x方向?yàn)橹芟?，y方向?yàn)檩S向。

圖1 接觸區(qū)微觀示意Fig.1 Microscopic schematic of the contact zone

為方便模型求解， 忽略某些對(duì)計(jì)算結(jié)果影響不大的因素， 提出了以下合理的簡(jiǎn)化和假設(shè)：

（1） 唇封宏觀性質(zhì)不受表面微觀變形影響；

（2） 唇封與軸接觸特征沿圓周方向周期性分布；

（3） 唇封兩側(cè)充滿密封介質(zhì)， 以方便計(jì)算泄漏量；

（4） 轉(zhuǎn)速恒定， 唇封與軸長(zhǎng)時(shí)間磨合， 密封狀態(tài)穩(wěn)定；

（5） 膜厚遠(yuǎn)小于軸徑， 故忽略油膜曲率的影響。

1.1 流體力學(xué)分析

雷諾方程是解決流體黏性流動(dòng)的基礎(chǔ)理論方程，對(duì)于潤(rùn)滑脂這種典型的非牛頓流體[6-7］， 開展流體壓場(chǎng)分布的研究仍依靠雷諾方程的建立與求解。 區(qū)別于潤(rùn)滑油等牛頓流體， 潤(rùn)滑脂不滿足牛頓內(nèi)摩擦定律，若沿用傳統(tǒng)推導(dǎo)方式， 其復(fù)雜的本構(gòu)方程會(huì)導(dǎo)致雷諾方程在二維層面難以計(jì)算。

文中利用流速分離法[8］推導(dǎo)得到的普適流體潤(rùn)滑方程可以有效解決一般非牛頓流體潤(rùn)滑問題， 普適流體潤(rùn)滑方程如下：

式中：h為流體膜厚；ρ為流體密度；t為時(shí)間；f-1為非牛頓流體本構(gòu)方程的反函數(shù)；pf為流體壓力；U0、Uh為周向上兩邊界速度；V0、Vh為軸向上兩邊界速度；W0、Wh為徑向上兩邊界速度。

目前廣泛應(yīng)用于潤(rùn)滑脂研究的本構(gòu)模型為Herschel-Bulkley 模型[9］， 然而密封過程中膜厚一般是微米級(jí)別的， 給定轉(zhuǎn)速下對(duì)應(yīng)的剪應(yīng)變率遠(yuǎn)大于屈服應(yīng)力， 為計(jì)算方便， 可將Herschel-Bulkley 本構(gòu)模型簡(jiǎn)化成Ostwald 模型[10-11］， 如下所示：

式中：τf為流體剪應(yīng)力；η為稠度系數(shù)；u為流體周向速度；n為流變指數(shù)。

潤(rùn)滑脂的平衡方程為

結(jié)合流速分離法可得Ostwald 模型對(duì)應(yīng)的反函數(shù)為

上述反函數(shù)代入到普適流體潤(rùn)滑方程， 即得到面接觸條件下簡(jiǎn)化的潤(rùn)滑脂二維雷諾方程：

預(yù)設(shè)量綱一化變量， 防止pf、h等參數(shù)量級(jí)上的差異導(dǎo)致方程求解誤差過大進(jìn)而不收斂：

式中：pr為參考流體壓力；Pf為量綱一流體壓力；hr為參考流體膜厚；H為量綱一流體膜厚；lx、ly為唇封和軸接觸寬度。

對(duì)潤(rùn)滑脂二維雷諾方程進(jìn)行量綱一化處理：

上述方程中A～D為系數(shù)項(xiàng)，Pf為未知項(xiàng)，H為常數(shù)項(xiàng)。 若計(jì)算網(wǎng)格節(jié)點(diǎn)為m×m， 為提高數(shù)值求解速度， 可利用矩陣法將上式轉(zhuǎn)換成系數(shù)項(xiàng)矩陣（m2×m2） ×Pf列向量（m2×1） ＝常數(shù)項(xiàng)列向量（m2×1）的形式， 通過矩陣求逆運(yùn)算準(zhǔn)確快速地得到Pf分布。

在已知pf和h后， 即可通過潤(rùn)滑脂在軸向上的質(zhì)量流量公式計(jì)算泄漏率：

式中：Qy為軸向質(zhì)量流量。

1.2 接觸力學(xué)分析

由于物體表面無法絕對(duì)光滑， 唇封與軸實(shí)際配合時(shí)， 會(huì)存在微觀層面的粗糙峰接觸擠壓行為， 其產(chǎn)生的接觸壓力和流體壓力共同起到支撐負(fù)載的作用。

文中選擇經(jīng)典的赫茲接觸模型[12-14］計(jì)算粗糙峰接觸壓力， 該模型為兩個(gè)彈性接觸的表面提供封閉形式的表達(dá)式。 其中接觸壓力計(jì)算公式如下：

式中：pc為粗糙峰接觸壓力；E′ ＝E／（1- ν2） 為等效彈性模量， 其中E為彈性模量，ν為泊松比；d*為粗糙峰的幾何重疊深度；R′為等效曲率半徑；r*為到粗糙峰接觸圓中心點(diǎn)的距離；a*為接觸區(qū)域半徑。

粗糙峰接觸產(chǎn)生的摩擦剪應(yīng)力計(jì)算公式如下：

式中：τc為粗糙峰剪應(yīng)力；f為摩擦因數(shù)。

1.3 變形力學(xué)分析

實(shí)際裝配條件下， 唇封還會(huì)受預(yù)緊力作用發(fā)生宏觀變形， 產(chǎn)生靜態(tài)接觸壓力psc。 對(duì)于唇封來說， 正是流體壓力、 粗糙峰接觸壓力和靜態(tài)接觸壓力三者的共同影響使得唇口發(fā)生變形。

考慮到唇口變形量非常小， 符合小變形理論， 可以認(rèn)為唇封的變形與施加的載荷呈線性關(guān)系， 因此可以采用影響系數(shù)法[15］描述變形與載荷的對(duì)應(yīng)關(guān)系。其中接觸表面任一計(jì)算節(jié)點(diǎn)i法向變形的計(jì)算公式為

式中：δn為法向變形量；In為法向影響系數(shù)矩陣；psc為靜態(tài)接觸壓力。

可以看出節(jié)點(diǎn)i的法向變形量由該點(diǎn)所在軸向列m個(gè)節(jié)點(diǎn)的法向受力狀況共同決定。

同理， 接觸表面任一計(jì)算節(jié)點(diǎn)i切向變形的計(jì)算公式為

式中：Is為切向影響系數(shù)矩陣。

通過變形力學(xué)分析， 計(jì)算不平衡應(yīng)力值下唇封的變形量， 用以修正膜厚， 反過來膜厚又會(huì)改變流體壓力和粗糙峰接觸壓力， 引起唇封變形， 故需經(jīng)過多次迭代， 最終使得應(yīng)力平衡、 膜厚穩(wěn)定。

1.4 固體力學(xué)分析

文中利用有限元軟件ABAQUS 對(duì)唇封開展固體力學(xué)分析， 建立實(shí)際裝配受壓條件下唇封系統(tǒng)幾何模型， 結(jié)合材料單軸實(shí)驗(yàn)數(shù)據(jù)， 獲取靜態(tài)接觸壓力和影響系數(shù)矩陣。

仿真具體流程： 首先參考真實(shí)尺寸構(gòu)建唇封二維軸對(duì)稱模型， 模型包括唇封、 旋轉(zhuǎn)軸、 定位夾具等；其次根據(jù)丁腈橡膠單軸試驗(yàn)測(cè)試數(shù)據(jù)定義材料屬性，材料評(píng)估后選定的應(yīng)變能函數(shù)為Ogden，N＝3， 其中模型系數(shù)Mu1 ＝-36.194 852 8， Mu2 ＝12.546 384 6，Mu3 ＝26.979 046 5， Alpha1 ＝2.965 641 8， Alpha2 ＝4.312 813 73， Alpha3＝1.569 973 01，D1＝0.012 090 53，D2＝D3＝0； 然后模擬裝配過程， 保證唇封的過盈量準(zhǔn)確， 并通過直接加壓法施加流體靜壓， 網(wǎng)格劃分后仿真計(jì)算， 得到唇封理論的結(jié)構(gòu)變形和應(yīng)力應(yīng)變狀態(tài)， 進(jìn)而提取唇封與軸之間的靜態(tài)接觸壓力。 為獲得影響系數(shù)矩陣， 對(duì)任一接觸區(qū)節(jié)點(diǎn)施加單位力， 獲得所有節(jié)點(diǎn)的位移值， 最終得到的位移矩陣即為影響系數(shù)矩陣。

1.5 熱力學(xué)與傳熱學(xué)分析

實(shí)際工作狀態(tài)下， 唇封和軸之間會(huì)發(fā)生摩擦生熱現(xiàn)象， 造成唇口局部溫升， 同時(shí)密封介質(zhì)和軸之間進(jìn)行對(duì)流換熱， 帶走部分熱量， 如圖2 所示。

圖2 唇封系統(tǒng)熱行為Fig.2 Thermal behavior of lip sealing system

為考慮溫度對(duì)唇封力場(chǎng)的影響， 文中在有限元軟件ABAQUS 中進(jìn)行了完全熱力耦合分析。 首先在初始流固耦合計(jì)算結(jié)果的基礎(chǔ)上考慮摩擦生熱作用得到熱載荷大小及作用位置； 然后在ABAQUS 中定義部件的材料屬性， 包括隨溫度變化的熱傳導(dǎo)率、 膨脹系數(shù)及應(yīng)變能函數(shù)； 接著設(shè)置對(duì)流換熱系數(shù)及環(huán)境溫度， 施加表面熱通量， 進(jìn)行溫度位移耦合分析， 最終得到唇封溫度場(chǎng)及應(yīng)力應(yīng)變狀態(tài)。 新得到的靜態(tài)接觸壓力重新代入流固耦合過程計(jì)算熱載荷， 重復(fù)上述熱力耦合過程， 直至唇封溫度及力學(xué)狀態(tài)不發(fā)生變化，完成流固熱耦合分析。

接觸區(qū)域的表面熱通量計(jì)算公式如下：

式中：Q為表面熱通量；Fs為流固耦合計(jì)算所得摩擦力；v為軸的線速度；S為密封接觸區(qū)面積。

摩擦熱按照一定比例分配到兩個(gè)物體上， 即熱流分配系數(shù)， 計(jì)算公式為

式中：Kq為熱流分配系數(shù)；λs、λr分別為固體、流體熱導(dǎo)率；ρs、ρr分別為固體、 流體密度；Cs、Cr分別為固體、 流體比熱容。

旋轉(zhuǎn)軸與流體發(fā)生強(qiáng)制對(duì)流換熱行為， 認(rèn)為周向上軸與周圍流體相對(duì)速度保持一致， 故可以使用等效外掠平板模型[16］對(duì)其換熱系數(shù)進(jìn)行計(jì)算。

式中：Re＝vL／νr為雷諾數(shù)， 其中L為特征尺度，νr為流體運(yùn)動(dòng)黏度；Pr＝νrρrCr／λr為普朗特?cái)?shù)。

1.6 流固熱耦合模型

文中通過上述流體力學(xué)、 接觸力學(xué)、 固體力學(xué)和熱力學(xué)與傳熱學(xué)的相互耦合， 建立了盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封流固熱耦合仿真模型。 圖3 展示了所建立的數(shù)值仿真模型的具體計(jì)算流程。

圖3 盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封流固熱耦合數(shù)值仿真模型計(jì)算流程Fig.3 Calculation flow of numerical simulation model of fluidsolid-heat coupling for main drive lip seal of shield machine

主要步驟如下：

（1） 輸入唇封基本參數(shù)， 利用有限元軟件ABAQUS 進(jìn)行固體力學(xué)分析， 獲得靜態(tài)接觸壓力和影響系數(shù)矩陣；

（2） 根據(jù)接觸區(qū)粗糙峰形貌分布， 定義初始膜厚， 結(jié)合有限元仿真結(jié)果， 分別對(duì)流體壓力及粗糙峰接觸壓力進(jìn)行計(jì)算， 根據(jù)耦合關(guān)系進(jìn)行迭代求解， 調(diào)整油膜厚度；

（3） 計(jì)算表面粗糙峰的切向變形， 更新接觸區(qū)粗糙峰形貌分布， 直至粗糙峰的切向變形誤差達(dá)到收斂標(biāo)準(zhǔn)；

（4） 根據(jù)摩擦力大小及作用范圍， 進(jìn)行密封系統(tǒng)熱力耦合分析， 得到摩擦生熱條件下靜態(tài)接觸壓力分布； 重新進(jìn)行上述過程， 直至唇口溫度誤差滿足收斂條件， 流固熱耦合計(jì)算完成， 輸出泄漏率等密封性能參數(shù)。

2 結(jié)果與討論

根據(jù)提出的流固熱耦合仿真模型， 對(duì)盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封進(jìn)行分析， 涉及的基本參數(shù)如表1 所示。

表1 仿真模型輸入?yún)?shù)Table 1 Input parameters of simulation model

初始有限元仿真結(jié)果如圖4 和圖5 所示， 分別是靜態(tài)接觸壓力分布和影響系數(shù)矩陣。

圖4 靜態(tài)接觸壓力分布Fig.4 Static contact pressure distribution

圖5 影響系數(shù)矩陣Fig.5 Influence coefficient matrix

數(shù)值求解受網(wǎng)格密度、 收斂準(zhǔn)則和表面復(fù)雜性影響， 根據(jù)網(wǎng)格無關(guān)化條件， 計(jì)算網(wǎng)格設(shè)置為97×97；為提高計(jì)算結(jié)果的準(zhǔn)確度， 設(shè)置收斂精度為5×10-4，同時(shí)利用確定性模型定義初始膜厚， 即膜厚由常數(shù)項(xiàng)以及唇和軸表面微變形產(chǎn)生的波動(dòng)項(xiàng)組成， 膜厚公式如下：

式中：A1為唇口表面粗糙峰半幅值；λ11為唇口表面形貌x方向波長(zhǎng)；λ12為唇口表面形貌y方向波長(zhǎng)；A2為轉(zhuǎn)軸表面粗糙峰半幅值；λ21為轉(zhuǎn)軸表面形貌x方向波長(zhǎng)；λ22為轉(zhuǎn)軸表面形貌y方向波長(zhǎng)；h0為膜厚中值。

結(jié)合上述基本參數(shù)構(gòu)建的初始膜厚如圖6 所示。

圖6 初始膜厚分布Fig.6 Initial film thickness distribution

將初始有限元仿真結(jié)果和基本參數(shù)代入到流固耦合模型中， 可求得無溫度場(chǎng)條件下唇封摩擦力等密封性能， 其中摩擦力矩為3.35 N·m， 泄漏率為-2.73 g／mL （代表無泄漏產(chǎn)生）。 利用摩擦力計(jì)算唇口熱通量， 結(jié)合有限元軟件ABAQUS 進(jìn)行熱力耦合分析，最終得到的溫度場(chǎng)分布， 如圖7 所示， 可以看出唇口溫度大約為55 ℃。

圖7 溫度場(chǎng)分布（℃）Fig.7 Temperature field distribution （℃）： （a） overall temperature field of lip seal and shaft； （b） local temperature field of lip

多次迭代后流固熱耦合仿真模型精度滿足收斂要求， 得到考慮溫度場(chǎng)后摩擦力矩、 泄漏率等密封性能參數(shù)， 其中摩擦力矩為3.14 N·m， 泄漏率為-2.47 g／mL （即泵送率2.47 g／mL）。 考慮溫度場(chǎng)后， 摩擦力矩下降， 泵送率降低， 主要因?yàn)榇娇跍囟忍岣呤沟貌牧宪浕?導(dǎo)致唇封徑向力減小。

圖8 所示為考慮溫度場(chǎng)前后唇封接觸壓力分布。 可以看出考慮溫度場(chǎng)后， 接觸壓力最大值由4.06 MPa 減小到3.52 MPa， 但接觸寬度由3.55 mm增加到3.61 mm， 考慮溫度的唇封接觸情況更符合實(shí)際。

圖8 考慮溫度場(chǎng)前后接觸壓力分布Fig.8 Distribution of static contact pressure before and after considering temperature field

圖9 展示了考慮溫度場(chǎng)前后唇封的應(yīng)力應(yīng)變狀態(tài)。 可見應(yīng)力應(yīng)變均有顯著變化， 其中最大Mises應(yīng)力值由2.3 MPa 減小到2.1 MPa， 最大主應(yīng)變由0.21 減小到0.19。 最大Mises 應(yīng)力和最大主應(yīng)變都集中于唇口部分， 受溫度場(chǎng)影響， 該區(qū)域材料變軟， 是導(dǎo)致最大Mises 應(yīng)力和最大主應(yīng)變減小的主要原因。

圖9 考慮溫度場(chǎng)前后唇封應(yīng)力應(yīng)變狀態(tài)Fig.9 Stress and strain state of lip seal before and after considering temperature field： （a） initial lip seal stress distribution （MPa）；（b） final lip seal stress distribution （MPa）； （c） initial lip seal strain distribution； （d） final lip seal strain distribution

3 結(jié)論

基于脂潤(rùn)滑條件下唇形密封流固耦合模型， 結(jié)合有限元軟件實(shí)現(xiàn)溫度場(chǎng)與唇封力場(chǎng)的完全耦合， 建立盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封流固熱耦合仿真模型。 主要結(jié)論如下：

（1） 考慮溫度場(chǎng)后， 唇封的摩擦力矩下降， 泵送率降低， 主要因?yàn)榇娇跍囟忍岣呤沟貌牧宪浕?導(dǎo)致唇封徑向力減小。

（2） 考慮溫度場(chǎng)后， 唇封接觸壓力最大值由4.06 MPa 減小到3.52 MPa， 但接觸寬度由3.55 mm增加到3.61mm， 考慮溫度的唇封接觸情況更符合實(shí)際。

（3） 考慮溫度場(chǎng)后， 最大Mises 應(yīng)力值由2.3 MPa 減小到2.1 MPa， 最大主應(yīng)變由0.21 減小到0.19。 最大Mises 應(yīng)力和最大主應(yīng)變都集中于唇口部分， 受溫度場(chǎng)影響， 該區(qū)域材料變軟， 是導(dǎo)致最大Mises 應(yīng)力和最大主應(yīng)變減小的主要原因。

（4） 溫度對(duì)唇封應(yīng)力應(yīng)變狀態(tài)及密封性能產(chǎn)生較大影響， 考慮溫度場(chǎng)后的預(yù)測(cè)結(jié)果更貼合實(shí)際， 這對(duì)盾構(gòu)機(jī)主驅(qū)動(dòng)唇形密封設(shè)計(jì)具有一定指導(dǎo)作用。