基于ABAQUS的載人潛水器觀察窗結構蠕變行為分析

杜青海，江海港，胡曉康

1 上海海洋大學 深淵科學工程技術研究中心, 上海 201306

2 上海海洋大學 海洋科學學院, 上海 201306

3 上海海洋大學 工程學院, 上海 201306

0 引 言

目前，陸上資源的日益匱乏讓全世界的目光開始轉向海洋資源。隨著對海洋資源的深入研究，各國對海洋的關注重點開始從近淺海逐漸向深遠海發展[1]，而載人潛水器（HOV）在其中將發揮關鍵性作用。觀察窗是載人潛水器的關鍵結構之一，其安全性對于人員生命安全以及潛水器性能具有十分重要的作用。相對于金屬艙體結構而言，透明玻璃材質觀察窗結構的設計更具復雜性和不確定性[2]，故國際上的通用規范[3]采取提高安全系數(4～16倍)的方法來確保結構安全。

國內外學者對于觀察窗的研究始于上世紀中期，其中最為著名的是美國的Stachiw[4]，其畢生都在進行該方面的研究，成果已被寫入設計規范[3]并受到國際上相關行業的認可。值得注意的是，Stachiw的研究主要集中在7 000 m潛深范圍以內，對于更大的深度，其采取的是增大安全系數按比例進行保守設計的方法。近10年來，我國載人深潛裝備隨著7 000 m級“蛟龍”號的成功研制，也得到了蓬勃發展，4 500 m級“深海勇士”號以及萬米級的“奮斗者”號已相繼完成研制。針對載人潛水器的觀察窗結構，我國科研人員在結構應力[5]、蠕變特性[6-7]和結構協調性[8]等方面開展了研究。對于全海深大深度觀察窗結構，筆者所在課題組通過數值方法，針對觀察窗結構規范方法開展了探索性研究[9]，指出該設計方法過于保守；隨后，又針對觀察窗結構進行了優化計算，以尋找最優方案[10]。但對于觀察窗的長時間周期和蠕變壽命等問題，尚未開展相關研究。

因此，本文擬采用ABAQUS軟件構建錐臺型觀察窗的有限元計算模型，然后基于窗玻璃材料（丙烯酸塑料材料）的基本性能確定蠕變模型參數，進而分析、驗證觀察窗結構的蠕變特性，最后，基于有限元蠕變模型計算分析不同加載速度下觀察窗結構的蠕變特性及結構性能，為觀察窗結構的疲勞壽命研究奠定基礎。

1 錐臺型觀察窗

1.1 觀察窗幾何結構參數

錐臺型觀察窗的基本結構參數如下：承受的外載荷P、厚度δ、玻璃低壓面（即內側圓面）直徑Di、窗座內口徑Df、錐角α，以及厚?徑比δ/Di和窗座的懸伸比Di/Df。如圖1所示，取錐臺型觀察窗的內側圓面為AB面，外側高壓圓面為EF面，選擇EF面的圓心O點作為柱坐標系rOz的原點，觀察窗軸對稱旋轉方向，也即錐臺觀察窗圓面的圓周方向記為θ，AB面的圓心為C點，OF長記為r。

圖1 錐臺型觀察窗的基本結構參數Fig.1 Basic structural parameters of conical observation window

1.2 材料的蠕變特性

由于潛水器觀察窗大多數使用丙烯酸塑料（PMMA）[4]，這種材料具有比較明顯的蠕變特性，能使應力保持在較低水平，且不會導致嚴重的應力集中現象，有力地推動了觀察窗的飛躍式發展。本文錐臺型觀察窗主要結構材料的機械物理屬性[10]如表1所示。

表1 丙烯酸塑料的材料力學參數[10]Table 1 Mechanical parameters of PMMA[10]

圖2給出了丙烯酸塑料試件在不同外載荷作用下拉伸后的應變ε隨時間t的變化曲線[11]，揭示了此窗玻璃材料固有的蠕變特性。

圖2 PMMA拉伸試件試驗應變隨時間的變化曲線[11]Fig.2 Variation of strain level of PMMA with time in tensile specimen experimental test[11]

2 錐臺型觀察窗蠕變有限元理論

2.1 蠕變理論模型

由于采用傳統的ABAQUS線彈性觀察窗模型得出的結果忽略了丙烯酸塑料蠕變這一特性，導致計算結果與真實的試驗數據間存在較大差異，無法作為普遍的計算方法，因此需要尋找一種能體現蠕變特性的有限元仿真計算方法。為此，本文將在進一步研究觀察窗蠕變特性的基礎上，采用單步靜力?線彈性與多步?非線性相結合的分析步控制方法，同時引入加載和保壓的時間，重點得出觀察窗的軸向位移、應力和應變等隨時間的變化規律。

ABAQUS軟件中有3種蠕變模型：雙曲正弦模型、應變硬化模型和時效硬化模型，可用來仿真模擬材料的蠕變等特性。其中，恒溫與定負載條件下的蠕變特性選擇時效硬化模型或應變硬化模型，變溫條件下的蠕變特性則使用雙曲正弦法則模型。應變硬化的基本思想是：在蠕變發生的全過程中，起到強化作用的主要原因是蠕變變形與時間無任何關系，故其適用于短時間的試驗。然而，載人潛水器在工作深度下長時間作業的觀察窗結構蠕變特性[6-7]與時間有很大的關系，故觀察窗結構蠕變分析可選擇時效硬化模型。

時效硬化理論[12]的公式可以表述為：當外界溫度一定時，材料的應變ε、應力σ和時間t這3個參數滿足的狀態函數Φ為：

式中，ε˙為應變隨時間變化的應變率。

由文獻[12]可知，在蠕變發生階段，影響蠕變的是時效、擴散及恢復這幾個參數，且因蠕變曲線在第1、第2階段具有幾何相似的性質，故

其中，應力函數f1(σ)符合下面的冪函數形式，即

式中，n為應力冪函數的指數參數。

而時間響應函數f2(t)則用Ω(t)函數表示，即

式中，A，m分別為時間冪函數的系數和指數參數。

因此，將式(3)和式(4)代入式(2)，可得到時效硬化模型下應變的數學表達式，即

2.2 蠕變模型參數

根據對觀察窗結構進行的線彈性數值模擬仿真分析與優化[4,9-10]，可知規范規定的錐臺型觀察窗結構所對應的最大工作壓力作用下的整體應力水平基本集中在30～45 MPa之間。鑒此，為得到式(5)中的A，m，n這3個待定系數，根據圖2所示的試驗數據[11]，盡量選擇高應力條件下的試驗值，即選取30～40 MPa間的材料拉伸蠕變曲線。并且，選擇在拉伸試驗曲線的兩端取點，對拉伸試件試驗值與計算值進行對比，如圖3(a)和圖3(b)所示。由圖可見，在這2種取點情況下所得的計算值與試驗值相差較大，故應在試驗值的中間部分取點，并需滿足二者之間誤差很小，在此取點情況下的試驗值與計算值對比如圖3(c)所示。由圖可見，擬合出的計算值與試驗值基本吻合，證明該組參數有效。具體取點方式如表2所示。

表2 確定參數的取點方式及其解Table 2 Point selecting schemes and solution of determined parameters

圖3 PMMA拉伸試件應變試驗值與計算值對比Fig.3 Comparison of strain level between expreimental test and calculation of PMMA tensile specimen

2.3 有限元模型

根據錐臺型觀察窗結構和均勻靜水壓力載荷的軸對稱性特點，對觀察窗結構選取軸對稱形式的二維平面模型進行仿真模擬計算。根據模型尺寸，將選取大小合適的CAX4R單元對模型進行網格劃分，并設置模型部件之間的接觸形式、窗體與窗座之間的摩擦系數和罰函數因子來描述接觸[10]。觀察窗仿真計算模型的載荷邊界約束分布及網格單元劃分如圖4所示。

圖4 ABAQUS有限元模型的載荷邊界約束與網格劃分Fig.4 Boundary constraints and grid division of finite element model on ABAQUS

2.4 有限元計算方法驗證

現有的行業標準ASME PVHO是建立在Stachiw[4]的系列試驗模型研究基礎之上的，故選取其典型模型尺寸，即主尺度參數。具體如下：Di=25.4 mm，δ/Di= 0.75，1.0，1.25，α=90°。根據真實試驗情況將計算壓力加載到規范規定[3]的界限壓力68.9 MPa，然后保壓100 h。

在觀察窗結構數值模擬仿真研究[9-10]中，已針對上述尺度的模型驗證了線彈性計算環境下，摩擦系數μ=0.05時有限元計算結果隨網格粗細變化的準確性及收斂性，得出在該尺度下觀察窗和窗座分別采用1 和3 mm的網格比較合適，故此處不再贅述。

運用ABAQUS軟件對上述尺度的模型建立幾何模型并進行模擬仿真分析，獲得加壓和保壓階段的應變與變形等，然后將數值模擬仿真所得內側圓面中心軸向位移fc及其與試驗結果的相對誤差Δ進行對比分析，結果如圖5所示。由圖5可見，在前100 h，有限元計算結果與試驗結果[4]吻合，誤差絕對量在10%以內；在100～1 000 h間，有限元計算結果稍大于試驗結果，二者之間的誤差保持在11%左右，這說明在此階段采用數值模擬仿真方法具有保守性?？梢?，本文仿真計算結果與試驗結果較相符。而實際工程應用中保壓的時間一般都小于100 h，觀察窗體結構因其玻璃材料的黏彈特性，在保壓的起始階段(0～2 h)其結構性能主要由材料的黏彈屬性決定[4]，故本文將著重選擇10～100 h這一階段進行研究。

圖5 有限元計算結果與試驗結果對比Fig.5 Comparison of finite element calculation and experimental test results

3 全海深觀察窗結構蠕變分析

3.1 高壓作用下的蠕變特性

現有行業規范[3]中的設計方法是通過安全系數法則來確保結構疲勞壽命安全，但卻無法揭示長周期載荷下觀察窗所具有的蠕變現象，導致計算結果與實測結果存在較大誤差[7,9,10]，且不利于結構的優化。因此，本文將基于時效硬化模型和試驗確定參數的蠕變計算方法深入研究觀察窗結構在高壓作用下的蠕變等特性。“蛟龍”號觀察窗的結構尺寸參數和最大工作壓力恰好是規范[3]中方法——曲線設計參數分界[9-10]的典型結構尺寸參數和壓力范圍，其基本幾何結構參數如表3所示。

表3 “蛟龍號”觀察窗主尺度參數[5,7]Table 3 Main dimensions of conical observation window of Jiaolong ship[5,7]

利用第2節的數值模擬仿真方法對觀察窗結構進行有限元建模和模擬仿真，并在計算過程中分別選擇線彈性和時效硬化的蠕變模型。在數值計算過程中，同時采用單步靜力?線彈性與多步?非線性相結合的分析步控制的方法加入時間非線性響應。最終，將作用在觀察窗模型上的壓力加載至71.6 MPa并保壓10 h，觀察窗結構軸向位移的計算結果，如圖6所示。圖中，Uz為結構剖面的軸向位移。

圖6 保壓10 h后觀察窗的軸向位移云圖Fig.6 Axial displacement contours of conical observation window after holding pressure for 10 h

圖6中給出了觀察窗結構在保壓10 h后線彈性和蠕變這2種模型下的軸向位移云圖，而圖7則給出了這2種模型在加載和保壓過程中觀察窗內側圓面中心點軸向位移fc隨時間的變化曲線。由圖7可見，采用線彈性模型計算出的軸向位移在保壓階段無任何改變，而采用蠕變模型計算出的軸向位移fc在保壓階段則是隨著時間的增長逐步增大，這直接反映出了觀察窗結構蠕變這一特性，與實際結果比較契合。圖8進一步展示了最大應變εmax隨時間的變化趨勢。由圖8可見，在保壓階段，采用蠕變模型計算的最大應變隨著時間的推移逐漸增大，而采用線彈性模型計算的最大應變在保壓階段則無任何變化，可見觀察窗玻璃材料的固有蠕變特性直接導致了觀察窗結構在高壓下的蠕變特性。

圖7 觀察窗軸向位移隨時間的變化曲線Fig.7 Variation of the axial displacement of conical observation window with time

圖8 觀察窗最大應變隨時間的變化曲線Fig.8 Variation of the maximum strain level of conical observation window with time

根據對深海觀察窗結構強度優化的研究[10]，發現在主尺度參數下，錐角對錐臺型觀察窗的結構優化與安全評估具有重要意義。因此，本文選取與表3所示觀察窗主尺度參數一致的4組結構模型進行分析，即Di=220 mm，δ/Di=1.0，α=70°，80°，90°和100°，將作用壓力加載至71.6 MPa并保壓100 h，計算得到4組模型結構的基本力學特性，如圖9～圖12所示。圖中， σ為 最大等效應力，τmax為最大剪應力。

由圖9可見，在壓力載荷與摩擦系數都相同的情況下，相同厚?徑比觀察窗內側圓中心處的軸向位移fc與錐角α成反比；由圖10可見，觀察窗同一時間處的最大等效應力 σ隨著錐角α的增大呈現先增大后減小的趨勢；由圖11可見，觀察窗同一時間的最大剪應力τmax隨著錐角α的增大呈增大的趨勢；而圖12則揭示了觀察窗同一時刻的最大應變εmax是隨著錐角α的增大呈先增大后減小的規律。綜合上述分析，可得錐角α=70°時觀察窗盡管在軸向位移方面其數值是4個模型中最大的，但其最大等效應力、最大剪應力和最大應變卻是最小的，且應力集中現象不明顯。因此，α=70°時的觀察窗模型是相對優化的結果，這與線彈性數值模擬優化研究結論[10]一致。

圖9 軸向位移隨時間的變化趨勢Fig.9 The changing trend of axial displacement with time

圖10 最大等效應力隨時間的變化趨勢Fig.10 The changing trend of maximum equivalent stress with time

圖11 最大剪應力隨時間的變化趨勢Fig.11 The changing trend of maximum shear stress with time

圖12 最大應變隨時間的變化趨勢Fig.12 The changing trend of maximum strain with time

3.2 不同加載速度下的蠕變特性

經過大量的仿真計算，發現不同的加載（下潛）速度對觀察窗的位移、應變及最大等效應力等輸出值均有影響，會影響結構的蠕變特性。

為便于研究，針對3.1節中優化結果，設Di為固定值220 mm，α=70°，δ/Di=1.0，摩擦系數μ=0.05，采取3種不同的加載速度，即0.67，2和6 MPa/min，然后將觀察窗模型的作用壓力P加載至71.6 MPa并保壓100 h，最后研究采取不同加載速度時對觀察窗結構基本力學性能的影響。

圖13和圖14給出了加載和保壓過程中觀察窗結構的位移和最大應變特性。由圖可見，在將作用壓力加載至71.6 MPa的過程中，3種加載速度下的位移和應變值相差無幾，只是隨著加載時間的縮短，二者略有降低；而在保壓階段二者差距卻十分明顯，盡管保壓時間相同，但隨著加載時間的增加，位移和最大應變在保壓階段增長的幅度反而降低了。圖15顯示出在整個加載過程中不同加載速度下的最大等效應力值幾乎一致，其變化主要體現在保壓階段，即隨著加載時間的增大，最大等效應力的降幅在減小。由圖15可見，在相同的長時間保壓作用下，延長加載時間相對縮短加載時間而言，與正常加載速度（例如2 MPa/min）相比，觀察窗模型結構計算的等效應力變化幅度受到的影響較大。綜上所述，壓力加載速度不宜過快，否則不利于高壓保壓階段的結構蠕變應變及位移釋放。

圖13 不同加載速度下位移的變化曲線Fig.13 Variation of displacement under different loading speeds

圖14 不同加載速度下最大應變的變化曲線Fig.14 Variation of maximum strain under different loading speeds

圖15 不同加載速度下最大等效應力的變化規律Fig.15 Variation law of the maximum equivalent stress under different loading speeds

3.3 觀察窗蠕變特性下的結構強度

為了更好地研究觀察窗材料蠕變特性下其強度的力學結構屬性，需進一步分析觀察窗在長時間保壓作用下的結構強度特性，這對觀察窗的疲勞壽命分析來說意義重大。因此，選取前述優化的觀察窗結構主尺度參數，即Di=220 mm，δ/Di=1.0，α=70°，再按2 MPa/min加載速度將作用壓力加載至71.6 MPa，并保壓100 h。

圖16(a)所示為保壓后觀察窗上表面（EF面）在3個方向的應變曲線。由圖可見，在EF面中心附近存在高應變集中區域，數值上呈現從中心向四周遞減的趨勢，圖中虛線范圍所表示的區域隨著外部加載的不斷增大，出現了“凹坑”情況[4]。圖16(b)所示為內側圓面（AB面）在3個方向（即圖1所示柱坐標系下）的應變曲線。由圖可見，內側圓面在邊緣處易出現因幾何結構和應變突變所導致的裂紋等。圖中 ，εr，εθ和 εz分別為3個坐標軸方向上的主應變。

圖16 觀察窗外/內表面的應變曲線Fig.16 The strain distribution of outside/inside surface of conical observation window

圖17所示為觀察窗保壓后的3/4立體應變云圖。從圖中可以直觀地看出，EF面中心O點附近3個應變分量的絕對值明顯高于四周，這也進一步揭示了在高壓下觀察窗上表面會出現“凹坑”的破損失效現象。

圖17 觀察窗結構的應變云圖Fig.17 The strain contours of conical observation window structure

圖18給出了在保壓過程中，在0，10，100 h這3個時刻觀察窗中心軸線OC上沿z軸3個方向的應變變化情況，圖中縱坐標為軸線OC位置坐標z與厚度δ的比值，橫坐標為應變值大小。從圖中可以清晰地看出，虛線區域內的應變量會隨著保壓時間的增加而持續增加，相比其他中心軸線OC上區域，該區域受到的應變突變更嚴重，也更容易出現損傷破壞。

圖18 觀察窗中心軸線OC的應變曲線Fig.18 The strain distribution along OC at the center axis of conical observation window

4 結 論

本文基于丙烯酸塑料材料的蠕變特性，采用ABAQUS軟件對錐臺型觀察窗結構進行了數值模擬仿真計算，并通過時效硬化模型對材料和觀察窗結構的蠕變特性展開了分析，主要得到如下結論：

1) 基于ABAQUS軟件，構建了用于評估錐臺型觀察窗在長周期載荷下材料與結構蠕變特性的可靠數值方法。

2) 觀察窗結構有限元線彈性模型與蠕變模型的比較分析，以及筆者所在課題組之前對蠕變模型與典型參數模型試驗的對比分析表明，蠕變模型下的計算結果更符合試驗情況。通過多參數優化計算分析觀察窗結構的蠕變特性，并在特定主尺度參數下優化對比錐角對結構性能的影響，結果顯示錐角α=70°時觀察窗結構的蠕變特性更優。

3) 不同加載（下潛）速度下觀察窗的蠕變強度研究顯示，為避免蠕變特性下的結構損傷，不宜采取過快的加載速度。

本文針對觀察窗結構的蠕變特性分析可為今后觀察窗結構的優化設計提供參考。