盾構機唇形密封流固熱耦合仿真研究*

紀佳馨 蘇世東 項 沖 郭 飛

（1. 中國石油大學（華東） 機電工程學院 山東青島 266580；2. 清華大學高端裝備界面科學與技術全國重點實驗室 北京 100084）

隨著國家城市化進程的發展， 盾構掘進技術在城 市地鐵、 鐵路交通等領域得到廣泛應用。 鑒于盾構機工作過程中特殊的隧道環境， 主驅動唇形密封系統[1-3］一旦失效， 會導致外界環境中的土砂倒灌進而破壞主軸承， 嚴重影響工程進度。 唇形密封件的性能直接決定了整臺盾構機的施工效率， 且在施工過程中應避免唇形密封件的更換， 為免影響施工效率， 其故在盾構機主驅動密封設計階段， 預測唇形密封性能是否滿足要求尤為重要。

國產盾構機主驅動密封件多依賴國外進口， 且價格高昂， 為提高國產大型盾構機的國際競爭力， 亟待開展盾構機主驅動密封研究。 譚鋒等人[4］利用正交試驗法針對盾構機主驅動唇形密封開展結構優化研究，根據密封工作要求確定最優密封結構， 優化后密封件接觸寬度降低23%。 張中華等[5］在主驅動密封失效分析的基礎上， 總結提出密封安裝結構和工藝優化方案， 改善盾構機主驅動密封性能。

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

1 數值仿真模型

唇封系統主要由密封唇、 旋轉軸和密封介質三部分組成， 其相互作用共同決定唇封系統的密封性能。如圖1 所示為唇封與軸接觸區的微觀示意圖， 可以看出建立盾構機主驅動唇形密封數值仿真模型， 需要綜合分析密封過程涉及的物理要素， 包括密封介質的流體力學分析， 微觀粗糙峰與軸之間的接觸力學分析和密封唇的固體力學分析。 圖中所示x方向為周向，y方向為軸向。

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

為方便模型求解， 忽略某些對計算結果影響不大的因素， 提出了以下合理的簡化和假設：

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

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

（3） 唇封兩側充滿密封介質， 以方便計算泄漏量；

（4） 轉速恒定， 唇封與軸長時間磨合， 密封狀態穩定；

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

1.1 流體力學分析

雷諾方程是解決流體黏性流動的基礎理論方程，對于潤滑脂這種典型的非牛頓流體[6-7］， 開展流體壓場分布的研究仍依靠雷諾方程的建立與求解。 區別于潤滑油等牛頓流體， 潤滑脂不滿足牛頓內摩擦定律，若沿用傳統推導方式， 其復雜的本構方程會導致雷諾方程在二維層面難以計算。

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

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

目前廣泛應用于潤滑脂研究的本構模型為Herschel-Bulkley 模型[9］， 然而密封過程中膜厚一般是微米級別的， 給定轉速下對應的剪應變率遠大于屈服應力， 為計算方便， 可將Herschel-Bulkley 本構模型簡化成Ostwald 模型[10-11］， 如下所示：

式中：τf為流體剪應力；η為稠度系數；u為流體周向速度；n為流變指數。

潤滑脂的平衡方程為

結合流速分離法可得Ostwald 模型對應的反函數為

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

預設量綱一化變量， 防止pf、h等參數量級上的差異導致方程求解誤差過大進而不收斂：

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

對潤滑脂二維雷諾方程進行量綱一化處理：

上述方程中A～D為系數項，Pf為未知項，H為常數項。 若計算網格節點為m×m， 為提高數值求解速度， 可利用矩陣法將上式轉換成系數項矩陣（m2×m2） ×Pf列向量（m2×1） ＝常數項列向量（m2×1）的形式， 通過矩陣求逆運算準確快速地得到Pf分布。

在已知pf和h后， 即可通過潤滑脂在軸向上的質量流量公式計算泄漏率：

式中：Qy為軸向質量流量。

1.2 接觸力學分析

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

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

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

粗糙峰接觸產生的摩擦剪應力計算公式如下：

式中：τc為粗糙峰剪應力；f為摩擦因數。

1.3 變形力學分析

實際裝配條件下， 唇封還會受預緊力作用發生宏觀變形， 產生靜態接觸壓力psc。 對于唇封來說， 正是流體壓力、 粗糙峰接觸壓力和靜態接觸壓力三者的共同影響使得唇口發生變形。

考慮到唇口變形量非常小， 符合小變形理論， 可以認為唇封的變形與施加的載荷呈線性關系， 因此可以采用影響系數法[15］描述變形與載荷的對應關系。其中接觸表面任一計算節點i法向變形的計算公式為

式中：δn為法向變形量；In為法向影響系數矩陣；psc為靜態接觸壓力。

可以看出節點i的法向變形量由該點所在軸向列m個節點的法向受力狀況共同決定。

同理， 接觸表面任一計算節點i切向變形的計算公式為

式中：Is為切向影響系數矩陣。

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

1.4 固體力學分析

文中利用有限元軟件ABAQUS 對唇封開展固體力學分析， 建立實際裝配受壓條件下唇封系統幾何模型， 結合材料單軸實驗數據， 獲取靜態接觸壓力和影響系數矩陣。

仿真具體流程： 首先參考真實尺寸構建唇封二維軸對稱模型， 模型包括唇封、 旋轉軸、 定位夾具等；其次根據丁腈橡膠單軸試驗測試數據定義材料屬性，材料評估后選定的應變能函數為Ogden，N＝3， 其中模型系數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； 然后模擬裝配過程， 保證唇封的過盈量準確， 并通過直接加壓法施加流體靜壓， 網格劃分后仿真計算， 得到唇封理論的結構變形和應力應變狀態， 進而提取唇封與軸之間的靜態接觸壓力。 為獲得影響系數矩陣， 對任一接觸區節點施加單位力， 獲得所有節點的位移值， 最終得到的位移矩陣即為影響系數矩陣。

1.5 熱力學與傳熱學分析

實際工作狀態下， 唇封和軸之間會發生摩擦生熱現象， 造成唇口局部溫升， 同時密封介質和軸之間進行對流換熱， 帶走部分熱量， 如圖2 所示。

圖2 唇封系統熱行為Fig.2 Thermal behavior of lip sealing system

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

接觸區域的表面熱通量計算公式如下：

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

摩擦熱按照一定比例分配到兩個物體上， 即熱流分配系數， 計算公式為

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

旋轉軸與流體發生強制對流換熱行為， 認為周向上軸與周圍流體相對速度保持一致， 故可以使用等效外掠平板模型[16］對其換熱系數進行計算。

式中：Re＝vL／νr為雷諾數， 其中L為特征尺度，νr為流體運動黏度；Pr＝νrρrCr／λr為普朗特數。

1.6 流固熱耦合模型

文中通過上述流體力學、 接觸力學、 固體力學和熱力學與傳熱學的相互耦合， 建立了盾構機主驅動唇形密封流固熱耦合仿真模型。 圖3 展示了所建立的數值仿真模型的具體計算流程。

圖3 盾構機主驅動唇形密封流固熱耦合數值仿真模型計算流程Fig.3 Calculation flow of numerical simulation model of fluidsolid-heat coupling for main drive lip seal of shield machine

主要步驟如下：

（1） 輸入唇封基本參數， 利用有限元軟件ABAQUS 進行固體力學分析， 獲得靜態接觸壓力和影響系數矩陣；

（2） 根據接觸區粗糙峰形貌分布， 定義初始膜厚， 結合有限元仿真結果， 分別對流體壓力及粗糙峰接觸壓力進行計算， 根據耦合關系進行迭代求解， 調整油膜厚度；

（3） 計算表面粗糙峰的切向變形， 更新接觸區粗糙峰形貌分布， 直至粗糙峰的切向變形誤差達到收斂標準；

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

2 結果與討論

根據提出的流固熱耦合仿真模型， 對盾構機主驅動唇形密封進行分析， 涉及的基本參數如表1 所示。

表1 仿真模型輸入參數Table 1 Input parameters of simulation model

初始有限元仿真結果如圖4 和圖5 所示， 分別是靜態接觸壓力分布和影響系數矩陣。

圖4 靜態接觸壓力分布Fig.4 Static contact pressure distribution

圖5 影響系數矩陣Fig.5 Influence coefficient matrix

數值求解受網格密度、 收斂準則和表面復雜性影響， 根據網格無關化條件， 計算網格設置為97×97；為提高計算結果的準確度， 設置收斂精度為5×10-4，同時利用確定性模型定義初始膜厚， 即膜厚由常數項以及唇和軸表面微變形產生的波動項組成， 膜厚公式如下：

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

結合上述基本參數構建的初始膜厚如圖6 所示。

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

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

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

多次迭代后流固熱耦合仿真模型精度滿足收斂要求， 得到考慮溫度場后摩擦力矩、 泄漏率等密封性能參數， 其中摩擦力矩為3.14 N·m， 泄漏率為-2.47 g／mL （即泵送率2.47 g／mL）。 考慮溫度場后， 摩擦力矩下降， 泵送率降低， 主要因為唇口溫度提高使得材料軟化， 導致唇封徑向力減小。

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

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

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

圖9 考慮溫度場前后唇封應力應變狀態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 結論

基于脂潤滑條件下唇形密封流固耦合模型， 結合有限元軟件實現溫度場與唇封力場的完全耦合， 建立盾構機主驅動唇形密封流固熱耦合仿真模型。 主要結論如下：

（1） 考慮溫度場后， 唇封的摩擦力矩下降， 泵送率降低， 主要因為唇口溫度提高使得材料軟化， 導致唇封徑向力減小。

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

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

（4） 溫度對唇封應力應變狀態及密封性能產生較大影響， 考慮溫度場后的預測結果更貼合實際， 這對盾構機主驅動唇形密封設計具有一定指導作用。