载人潜器观察窗通常采用高分子材料有机玻璃，恒定压力下将随时间发生蠕变，由接触挤压变形与蠕变变形引起的严重变形不协调，很可能导致观察窗结构突然失效。针对接触边界的观察窗结构设计与优化问题，已有学者进行了深入研究与分析，而有关观察窗蠕变特性的研究则不多见。文献[1-5]从理论和试验两方面分别研究了不同设计参数的观察窗在不同工况下的临界应力和轴向位移。刘道启等^[6]进行了观察窗蠕变试验和数值研究，结果表明，观察窗蠕变与负载、结构形式、接触条件等有关。田常录等^[7]运用弹塑性力学全量理论结合试验数据推导出了观察窗蠕变变形公式。杨青松等^[8]通过试验分析了半球形观察窗蠕变特性，并与数值模拟结果进行了对比，但未给出数值模拟的具体方法，也未深入研究观察窗设计参数对蠕变特性的影响。目前，对观察窗蠕变特性的理论与试验研究结果表明，观察窗几何参数和接触面摩擦系数等因素与观察窗蠕变有一定关系，但未对其影响规律和敏感性做深入研究，未形成一套有效的数值分析方法。为了准确把握观察窗蠕变特性，避免出现因观察窗设计不合理导致的结构突然失效、变形不协调等安全性问题，有必要对观察窗蠕变特性进行深入研究。因此，本文基于ABAQUS软件，给出一种观察窗蠕变有限元分析的方法，该方法能够有效模拟载人潜器观察窗在一定工作压力下，一定工作时间内的蠕变行为，并研究锥台形观察窗设计参数对蠕变变形和应力的影响规律及灵敏度。

1 观察窗蠕变有限元分析方法

蠕变是结构在恒定负载下应力与应变随时间变化的一种现象。深海中，观察窗变形主要由与窗座的受压变形及蠕变变形两部分组成。考虑到有机玻璃的粘弹性特性，其应力-应变并非一一对应关系，因此不能简单地按弹性材料来模拟计算。本文运用ABAQUS进行观察窗蠕变分析，分析流程如图 1所示。

1.1 有机玻璃本构关系

张志林^[9]通过系列试验拟合出常温下有机玻璃回归公式：

$ {\rm{lg}}\left( {\frac{\varepsilon }{{{\varepsilon _c}}}} \right) = 0.036\;8{\left[ {{\rm{lg}}\left( {\frac{t}{{{t_c}}}} \right)} \right]^2} + 0.27811{\rm{lg}}\left( {\frac{t}{{{t_c}}}} \right) $

(1)

式中:ε为应变；ε_c为临界断裂应变；t为时间；t_c为临界断裂时间；ε_c和t_c与应力水平相关，文献[7, 9]等在常温下回归出ε_c和t_c表达式：

$ {\varepsilon _c} = 1.4{\rm{exp}}\left( { - \frac{\sigma }{{15}}} \right) $

(2)

$ {t_c} = \left( {0.5E + 14} \right){\rm{exp}}\left( { - \frac{\sigma }{{1.68}}} \right) $

(3)

式中：E为弹性模量；σ为应力。将式(2)、(3)代入式(1)得出有机玻璃蠕变试验应力-应变关系式：

$ \begin{array}{l} {\rm{lg}}\varepsilon = 0.036{\rm{ }}8({\rm{lg}}t + 0.021{\rm{ }}94\sigma )2 + \\ \;\;\;\;\;\;\;\;0.278{\rm{ }}1({\rm{lg}}t + 0.021{\rm{ }}94\sigma ) - 0.028{\rm{ }}95 \end{array} $

(4)

式(1)~(4)可作为观察窗有机玻璃本构模型。

1.2 蠕变模型选取

ABAQUS软件中有3种蠕变模型可用来模拟粘弹性材料的蠕变行为，恒温与定负载下的蠕变行为通常选择时效硬化模型及应变硬化模型；变温条件下则使用双曲正弦法则来模拟温度相关的稳态蠕变行为。

根据与时间有关的Norton幂次法则，蠕变前2期可表示为：

$ \varepsilon \left( {\sigma , t} \right) = \frac{A}{{m + 1}}{\sigma ^n}{t^{m + 1}} $

(5)

式中A、m、n为系数，根据试验数据确定。

将式(5)对时间微分可获得应变率：

$ \frac{{{\rm{d}}\varepsilon }}{{{\rm{d}}t}} = A{\sigma ^n}{t^m} $

(6)

式(6)为时效硬化模型应变变化率的数学表达式。研究载人潜器作业时间内定负载下的观察窗蠕变特性则选择时效硬化模型。

1.3 接触定义

接触对之间的Kimematic接触条件^[10]为：对于主接触面Г_s上的任意一点P_s，从接触面Г_c上在变形方向上的最近接触点P_c可以通过它们之间的相对距离确定，其距离表达式为：

$ {P_c} - {P_s} = \mathop {{\rm{min}}}\limits_{{P_s} \in {\mathit{\Gamma }_s}} {(^t}{P_s}{ - ^t}{P_s}) $

(7)

该式为非线性方程，可以通过牛顿-拉菲森法求解。在t+Δt时刻，主从接触面之间的距离可以表达为

$ ^{t + \Delta t}G\left( P \right) = {[^{t + \Delta t}}{P_c}{ - ^{t + \Delta t}}{P_s}]{ \cdot ^{t + \Delta t}}\mathit{\boldsymbol{N}} \ge 0 $

(8)

式(8)可用线性表达式来表示：

$ ^{t + \Delta t}G\left( P \right){ = ^t}G\left( P \right) + [{\Delta ^t}\mathit{\boldsymbol{u}}({P_c}) - {\Delta ^t}\mathit{\boldsymbol{u}}({P_s})]{ \cdot ^t}\mathit{\boldsymbol{N}} \ge 0 $

(9)

式中：Δ^tu(P)为位移矢量增量；^tN为t单位法向向量。

根据Koulomb准则，总摩擦力与粘连接触及滑动接触相关，即：

$ f\left\{ \begin{array}{l} \alpha {{f'}_s}, \;\;\;\;粘连接触\\ {{f'}_s}, \;\;\;\;\;\;\;滑动接触 \end{array} \right. $

(10)

式中：α为系数，0＜α≤1；f_s′=μN，μ为滑动摩擦系数，N为接触正压力。最大粘连摩擦力f_s′与观察窗及窗座的材料、接触面润滑等有关。载人潜器观察窗与窗座为面-面接触，采用接触对算法，软质材料有机玻璃观察窗作为柔性体，硬质材料窗座作为刚性体，接触面摩擦系数取决于接触面加工打磨、涂润滑油脂的情况^[11]。

1.4 有限元设置

为保证单一变量原则，排除窗座的干扰，边界条件须约束窗座y方向自由度及转动自由度。计算采用多步分析，设置一个静力、通用分析步和多个粘性，依次进行接触非线性分析和蠕变分析，分析步时间为标准作业时间，根据所研究潜器任务要求进行设置。

1.5 网格划分

遵循接触非线性分析中接触主面网格尺寸大于从面的原则，主从面网格要有一定的疏密区别，防止接触穿透，尽量采用高精度网格单元进行模型划分，如图 2所示。

1.6 方法可信度验证

根据文献[8]试验模型建立有限元模型, 如图 3所示，以验证本方法可信度。施加6 MPa荷载，保压6 h所得蠕变位移结果与试验结果对比如图 4所示。

文献[8]试验测得观察窗下表面圆心处蠕变量为0.11 mm；计算得到观察窗顶部下表面中心点最大蠕变位移1.233 mm，蠕变量0.113 mm，与试验结果相比误差为2.7%，同时由图 4可以看出，蠕变位移变化曲线与试验结果变化曲线趋势基本一致，变化相差不大，位于误差范围内，充分证明了有限元计算方法的可信度。

2 锥台形观察窗蠕变计算与分析

锥台形是载人潜器观察窗广泛采用的结构形式，开孔较小、密封性能良好以及视野较广是其主要优点，因此，将运用有限元软件对锥台形观察窗进行蠕变特性分析。

2.1 有限元模型与材料参数

参照我国“蛟龙”号载人潜器锥台形观察窗结构形式进行有限元计算分析，观察窗结构和数值模型如图 5、6所示，观察窗倾角α=45°，下表面直径D₁=200 mm，上表面直径D₂=640 mm，厚度t=220 mm。根据结构对称性而采用轴对称模型。观察窗及窗座材料性能见表 1。

2.2 网格收敛性检验

计算结果随网格加密变化不大时，可认为已接近收敛。由表 2中数据可知：窗座网格大小为20 mm×20 mm时已经收敛；观察窗网格大小为3 mm×3 mm时接近收敛。为降低计算成本，选择窗座网格大小为20 mm×20 mm，观察窗网格大小为3 mm×3 mm。采用高精度四边形结构网格单元进行模型划分如图 7所示。

2.3 参数影响性分析

通常载人潜器观察窗与窗座接触面摩擦系数通过打磨、涂润滑脂安装可以控制为[0.001, 0.3]。参考中国船级社《潜器规范》^[12]，观察窗几何参数及接触面摩擦系数变化范围如下：

$ \begin{array}{l} 30^\circ < \alpha \le 60^\circ , 0.4 < t/{D_1} \le 1.2, \\ \;\;\;\;\;\;\;\;\;0.001 < \mu \le 0.3 \end{array} $

式中：ɑ为观察窗倾角；t/D₁为厚度与底部直径比；μ为接触面摩擦系数。固定观察窗底下表面半径D₁=100 mm，通过改变厚度-直径比t/D₁来实现厚度变化。设计参数α、t/D₁、μ分别选取20、7、7组，设计响应为观察窗下表面轴线一点的蠕变量、蠕变前后轴向位移及观察窗最大Mises应力值。计算结果如图 8~13所示。

由图 8可知，观察窗蠕变前后的轴向位移和最大Mises应力均随厚度-直径比的增大而减小，说明厚度-直径比的增加有利于减小观察窗的变形。对比图 9发现，厚度-直径比为0.4时，观察窗下表面接触尖角处有较大的应力集中，最大Mises应力为148.1MPa，最大剪应力为93.71 MPa；厚度-直径比为0.9时，尖角处应力集中较小，最大Mises应力99.7 MPa，最大剪应力59.20 MPa。比较发现，厚度-直径比由0.4增加至0.9，其下表面接触面尖角应力集中明显改善，应力分布更均匀，最大Mises应力减小48.4 MPa，最大剪应力减小34.51 MPa，因此，增大厚度-直径比使观察窗应力分布更均匀是改善应力集中现象的主要原因。

从图 10可知，观察窗蠕变前后的轴向位移随倾角的增大而增加，而观察窗应力随观察窗倾角的增大总体呈增大趋势，变化幅度不大，局部存在极小值。观察窗倾角从30°~32°时，应力减小；32°~44°时逐渐增大；44°~54°时逐渐减小，大于54°时逐渐增大。

由于观察窗应力整体变化范围不大，现仅对倾角30°和32°时的观察窗进行比较分析，如图 11所示。倾角为30°时，应力集中位于观察窗下表面与接触面相连的尖角处；倾角为32°时，应力集中的位置沿接触面向上偏移小段距离，该处相对于尖角处应力有所减小，同时观察窗剪应力有所减小，最大剪应力由54.58 MPa减小至50.32 MPa。综合来看，认为应力集中位置偏离尖角及剪应力减小是其应力变化的主要原因。

从图 12可知，随着接触面摩擦系数的增加，观察窗蠕变前后的轴向位移和最大Mises应力均呈减小的趋势，蠕变后的观察窗应力变化幅度不大。比较图 13中不同摩擦系数下的观察窗剪应力云图发现，接触面光滑粗糙与否对观察窗内部剪应力分布影响不大，摩擦系数增大，最大剪应力下降幅度较小，应力集中情况差别不大。

2.4 参数灵敏度分析

利用ISIGHT优化软件对多组参数值进行计算, 通过试验设计能够获得设计变量灵敏度，确认变量间的交互效应。统计学里，被研究的各项参数称为因子；因子被研究的范围称为试验范围；在DOE研究中定义的因子数值称为水平；试验结果称为响应或输出^[13]。试验设计方法包括参数试验、全因子试验、拉丁方、最优拉丁超立方等，根据不同的问题进行选择^[14]。由于本文设计变量较少，选用最优拉丁超立方法，该方法基于随机拉丁方法进行了改进，优化了试验设计矩阵每列中各个水平出现的次序，使得设计矩阵中各样本点在样本空间均匀取值^[15-16]。空间采样时最小样本数应满足:

$ n = \left( {N + 1} \right) \times \left( {N + 2} \right)/2 $

(11)

式中N为因子数。

考虑到精度要求及计算量，样本点数目取20。结果Pareto图如图 14所示，Pareto图反映设计参数对应的贡献程度百分比^[17]，包含各参数线性项、二阶项的主效应及各个参数间的交互效应。图中显示了各参数对观察窗蠕变位移、蠕变量和最大Mises应力的贡献率。可以看出，观察窗厚度-直径比对观察窗蠕变位移、蠕变量及最大Mises应力影响最大，贡献率均在20%以上。观察窗倾角及接触面摩擦系数存在较大的交互效应，二者同时变化对目标响应影响较大，贡献率均在10%左右。

3 结论

1) 观察窗蠕变变形和应力随厚度-直径比增大而减小，接触面应力分布不均和尖角处应力集中现象明显改善；观察窗蠕变变形随倾角增大而增大，应力总体呈增大趋势，变化幅度不大；观察窗蠕变变形和应力随接触面摩擦系数增大而减小，增大到0.150后应力变化趋于平缓。

2) 观察窗厚度-直径比对轴向位移、蠕变量和最大应力值影响最大，观察窗倾角及接触面摩擦系数存在较强的交互效应，二者同时变化时影响也较大，在优化时应主要考虑观察窗的厚度-直径比。