有机玻璃观察窗的蠕变特性及数值模拟

2019-06-05郭大猷黄小平

上海交通大学学报 2019年5期

郭大猷，黄小平，王芳

(1.上海交通大学高新船舶与深海开发装备协同创新中心，上海 200240；2.上海海洋大学上海深渊科学工程技术研究中心，上海 201306)

载人深潜器的研发对于我国开发利用深海资源有着重要的意义，其中观察窗作为重要的开口结构，对于保障深潜器和作业人员的安全至关重要.观察窗的材料一般采用有机玻璃(PMMA)，其力学特性不同于通常的线性或者非线性材料，表现为对环境和加载方式敏感.

若严格按照ASME规范的设计要求，全海深潜水器载人舱的观察窗厚度可能会超过设计限度，严重地限制了载人深潜器的发展.现阶段，观察窗的设计需要大量的数值模拟和实验验证.基于有限元法(FEM)的模拟由于能够处理非线性问题(材料非线性、大变形、接触等)而得到广泛应用.

刘道启等[1]指出，观察窗与窗座之间的相对位移主要由两部分构成：① 观察窗与窗座在海水压力作用下发生的挤压变形；② 观察窗玻璃随时间的蠕变变形.由常温下PMMA蠕变的回归公式推导出反映蠕变变形的本构模型，采用有限元软件得到的观察窗随时间蠕变的计算结果，与实验分析结果接近.田常录等[2]针对潜水器观察窗的蠕变变形理论展开进一步的探索，获得了观察窗受力变形以及蠕变变形的计算公式,计算结果与实验结果较为吻合.可见，针对保载阶段蠕变的有限元模拟研究具有一定的准确性.

然而，随着潜深的增加，下潜过程中观察窗将受到更长时间、更高压力的静水载荷，这一阶段的蠕变影响将变得十分重要.因此，为了更贴切地模拟观察窗的作业流程，亟需一种能分析加压-保载全过程的有限元方法.

首先必须获得PMMA的材料模型，目前国内外的研究大多采用黏弹性本构关系模型.谢中秋等[3]在恒应变率加载单轴压缩实验中使用了含应变率效应的ZWT(朱王唐)黏弹性本构模型，并基于实验数据获取了PMMA材料的本构模型参数.Hu等[4]研究了温度和应变率对PMMA力学性能的影响，提出了考虑温度、应变率耦合影响的本构模型，研究结果与实验结果吻合良好.Wang和Uzair等[5-6]对ZWT模型进行改进，在考虑温度影响的情况下对PMMA的碰撞过程进行了建模分析.

考虑到观察窗的工作环境相当于一个准静态加载过程，本文借鉴文献 [7] 的研究思路，采用广义Maxwell模型表征PMMA的黏弹性材料特性，得出Prony级数形式的松弛弹性模量，进而结合有限元软件Abaqus，采用质量缩放技术和等效时长方法加快了显式动力计算过程，计算并校验了观察窗的蠕变变形.此外，针对观察窗模型局部应力集中问题，通过设置倒角有效降低了局部应力，同时使得网格的收敛性更好.

1 观察窗蠕变的黏弹性特性

当物体受到恒定外力时，其应力与变形随时间变化的现象称为蠕变.在海水压力作用下，观察窗变形主要包括与窗座间的挤压变形和蠕变变形.考虑蠕变特性时，PMMA的应力-应变并不是简单的线性关系，故不能仅以线弹性特征来描述其本构关系.

学界对PMMA等聚合物的认识经历了一个逐步变化的过程，在模拟研究聚合物变形行为的时候曾将聚合物建立为弹塑性材料模型、超弹性材料模型、黏塑性材料模型，而近几年随着对聚合物力学性能研究的逐步深入，发现黏弹性材料模型更符合聚合物的变形特点[7].

1.1 PMMA的黏弹性力学模型

黏弹性的最大的特性就是与应变率相关，有蠕变、应力松弛等现象，并对时间和温度有很大的依赖性[7].目前对PMMA黏弹性特征的表述可归结为两类：① 时温等效性；② 在不同应变率下的松弛或蠕变行为.所谓时温等效性，即可以通过确定位移因子用较高温度的松弛实验曲线平移得到低温下的松弛实验曲线，而位移因子可以由经典的Williams-Landell-Ferry(WLF)方程计算[7].另一方面，应变率的影响则表现为：恒定的实验温度下，随着应变率增大，PMMA的流动应力显著增加，且可以采用ZWT黏弹性本构模型进行拟合[3].

为表征聚合物这种黏弹性材料的力学行为，可以将其看成是由表征黏性的元件(用黏壶表示)和表征弹性的元件(用弹簧表示)组合起来的复合模型.由黏性和弹性元件的不同组合排列方式即可得出模拟黏弹性材料特性的各种模型[7].

最常用的表征黏弹性的模型是Maxwell模型和Kelvin模型，但由于其参数较少，难以准确表达较为复杂的蠕变规律.

黏弹性力学模型中的弹性和黏性单元越多，越能精确地表征出材料的黏弹性特性.本文所采用的广义Maxwell模型由n个Maxwell体和1个弹性元件并联而成，如图1所示.图中：Ei和ηi(i=1,2,…,n)分别为各Maxwell体中弹性元件的弹性模量及黏壶的黏度；Ee为单独弹性元件的弹性模量.

图1 广义Maxwell模型Fig.1 Generalized Maxwell model

考虑松弛问题，即从某初始时刻起对模型施加恒定应变，并且施加到并联的每个Maxwell体上的应变大小相同.对于Maxwell体，将τi=ηi/Ei称为松弛时间，基于Maxwell模型的应力松弛公式，各 Maxwell 元件的应力随时间的变化应满足

σi(t)=ε0Eie-t/τi

(1)

故广义Maxwell模型的总应力为

(2)

从而瞬时的松弛模量

(3)

由以上的分析可知，广义的Maxwell模型可以很好地表征聚合物的蠕变(或松弛)变形特征，且上述松弛模量的数学表达形式与Prony级数形式相近.Abaqus中的线性黏弹性分析正是基于松弛模量的Prony级数形式的积分式本构方程[7]的，据此可对观察窗蠕变特性进行分析.

1.2 Abaqus中的线性黏弹性

有限元软件Abaqus提供了线弹性、塑性、黏弹性等多种本构模型，结合静态或显式动态的求解方法，可模拟PMMA观察窗在加压-保载条件下的变形过程.本节重点介绍Abaqus中的黏弹性模型.

Abaqus提供了3种黏弹性的定义方式：① 基于频域的黏弹性函数——适用于稳态、微幅振动分析；② 基于时域的黏弹性函数——适用于依赖时间的分析；③ 基于用户指定的蠕变法则——适用于非线性黏弹性分析.本文采用基于时域的黏弹性本构关系来研究观察窗的蠕变问题.

1.2.1适用条件在定义连续弹性体材料的性质时，Abaqus的时域黏弹性模型假定切变和体积的变化行为在多轴应力状态下是相互独立的(不包括橡胶式泡沫(elastomeric foam)材料)，本构模型要能反映切变变形和体积变形的规律.黏弹性本构模型需要与线弹性，或类橡胶类物质的超弹性行为协同使用，并适用于大应变问题.

黏弹性主要在准静态分析、隐式直接积分法的动态分析、显式动态分析、稳态流动(transport)分析等分析类型中使用.黏弹性参数的输入主要有两种方法：① 按照归一化的Prony级数格式输入材料参数；② 直接输入经标准化的松弛或蠕变实验数据，由软件自行计算生成本构模型.

1.2.2数值计算 Abaqus通过对无量纲的松弛模量作Prony级数展开来定义黏弹性，如切变模量可表示为

(4)

对体积模量作Prony级数展开，可获得类似的数学形式：

(5)

且Abaqus假设

(6)

瞬时模量与初始模量、无量纲瞬时模量的关系为

GR(t)=G0gR(t)

(7)

KR(t)=K0kR(t)

(8)

在弹性范围内，切变模量和体积模量与弹性模量和泊松比的关系为

(9)

为计算黏弹性响应，需要在Abaqus中输入线弹性与线性黏弹性的相关参数.

2 观察窗加压-保载有限元模型

Wang等[8]开展的观察窗加压-保载实验具有缓慢加载(不超过 4.5 MPa/min)、高压下长时间稳压保载(5 h)的特点，且实验过程中温度保持不变.因此，可对观察窗PMMA材料的本构模型作出一定的简化，即在恒温条件下忽略WLF方程所考虑的时温等效性，并将整个实验过程视为准静态加载，不考虑在动态加载下较高应变率对本构关系的影响.

本文参考文献[8-9]的实验，建立用于观察窗加压-保载分析的有限元模型，并验证网格的收敛性.以观察窗的下表面圆心为观察点，对比线弹性、线性黏弹性两种材料特性下的观察窗变形情况，并与文献中的实测数据进行对比.

2.1 有限元模型的建立

锥台形观察窗结构如图2所示.图中：观察窗锥角为90°；下表面直径为130 mm；厚度为153 mm.本文使用两种模型：旋转轴对称平面模型(2D模型)和立体模型(3D模型).前者在对称轴截面上取1/2模型，后者取1/4模型，分别如图3和4所示，

图2 深潜器的PMMA观察窗设计尺寸 (mm)Fig.2 PMMA frustum design for deep-sea manned submersible (mm)

图3 观察窗2D模型Fig.3 2D model of observation window

图4 观察窗3D模型Fig.4 3D model of observation window

2.2 定义接触和约束

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

(10)

式(10)为非线性方程，可以通过Newton-Raphson法求解.在t+Δt时刻，主从接触面之间的距离可以表示为

t+ΔtS=[tPc-tPs]·t+ΔtN≥0

(11)

式中:t+ΔtN为t+Δt时刻的法向接触体内的单位法向向量.式(11)可用线性表达式表示为

t+ΔtS=tS+

[Δtu(Pc)-Δtu(Ps)]·tN≥0

(12)

式中：Δu(P)为P点的位移矢量增量.

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

(13)

本文以刚度大的窗座面作为主面，观察窗作为从面，窗座与窗体相接触并设置摩擦，取摩擦因数为 0.1.在模拟过程中，接触方向总是主面的法线方向，从面上的节点不会穿越到主面，但主面上的节点可以穿越从面.对窗座底部施加固定约束，并在截面上建立对称约束.

2.3 网格划分及收敛性分析

基于线弹性本构关系开展网格收敛性分析，PMMA的弹性模量取 2.76 GPa，泊松比取 0.35(供货商数据).模型外侧施加115 MPa的静水压力作为载荷输入.

对于3D模型，选用含有缩减积分和沙漏控制的8节点线性六面体单元(C3D8R)划分网格，遵循接触非线性分析中接触主面网格尺寸大于从面的原则.静力条件下试算得到的应力云图如图5所示.

图5 静力试算得到的应力云图 (MPa)Fig.5 Stress contour plot of trial static analysis (MPa)

通过观察应力分布情况，对窗体与窗座接触的应力集中区域进行局部网格细化，以期得到更稳定的计算结果，网格划分如图6所示.

图6 局部加密2 mm(过渡到4 mm)的3D模型Fig.6 Local grid refined 3D model (2—4 mm)

表1 无圆角3D模型在不同网格密度下的应力及位移
Tab.1 Stress and displacement of none-filleted 3D model with different mesh densities

网格尺寸/mmσ/MPaτ′12/MPau/mm6111.339.09-4.437489.5445.13-4.3472116.966.55-4.335

为了加快收敛分析的计算速度，转而采用2D模型，使用旋转轴对称(Axisymmetric)单元划分网格，并通过类似3D模型的方法对网格局部加密，如图7所示.对于无圆角过渡的模型，在网格加密的过程中，观察点位移值趋于稳定，但局部最大应力难以收敛，见表2.

于是尝试给出了含圆角过渡(半径r=20,40 mm)的模型，如图8所示.计算结果见表3.结果发现：① 圆角存在时，随着网格的逐步加密，局部应力更易趋于稳定，故能较好地解决应力集中问题；② 圆角半径越大，局部应力水平越低，观察点位移越大.

通过2D与3D模型计算结果的比较，不难发现：① 无倒角模型存在局部应力集中，Mises应力和切应力对网格尺寸、形状敏感性较强，难以横向比较；② 倒角相同的3D和2D模型，其观察点处位移趋近，误差在 1.5% 左右，并且3D模型位移较大；③ 随着网格的细化，3D和2D模型的位移都呈现出减小的趋势.

图7 局部网格加密的2D模型Fig.7 Local grid refined 2D model

图8 含圆角过渡的2D模型网格Fig.8 Filleted 2D model mesh

表2 无圆角2D模型在不同网格密度下的应力及位移Tab.2 Stress and displacement of none-filleted 2D model with different mesh densities

表3 存在圆角时2D模型在不同网格密度下的应力及位移
Tab.3 Stress and displacement of 20 mm-filleted 2D model with different mesh densities

r/mm网格尺寸/mmσ/MPaτ′12/MPau/mm202.0143.475.03-4.5401.0129.171.61-4.4310.5131.974.17-4.4250.2133.475.98-4.4200.1133.876.17-4.419402.0117.264.03-4.6831.0121.568.82-4.6790.5122.770.00-4.6670.2123.570.37-4.6620.1123.770.51-4.661

因此，为了节省计算资源，使用由较粗网格的3D和2D模型计算得到的观察点位移是偏安全的，且误差在可以接受的范围内.为了与实验相匹配，以下有限元建模都采用无倒角模型.

3 有限元计算及实验结果比较

根据文献 [8] 的实验描述，有限元建模时，设加压时长为 4.0×103s，且为准静态加载，当压力爬升至115 MPa时开始保载状态，保载时长 1.8×105s.为了加快计算，本例中采用质量缩放技术和等效时长方法.质量缩放是指通过增加非物理质量到模型单元上从而获得大的显式时间步长的技术，对于不计应变率的准静态加载问题，采用质量缩放系数加快计算是有效的手段[11].等效时长方法，是基于准静态加载的特点，对黏弹性松弛时间参数和加载-保载时间同时除以103，使得在较短的时间步内求得位移变形历程，而不改变位移变形的大小.

3.1 线弹性模型

基于线弹性本构关系的有限元分析，不计加载-保载时间的影响，难以反映观察窗在高压加载-保载过程中的蠕变变形.据此校核观察窗位移及由此带来的密封问题，会对观察窗产生的位移和变形估计不足.实验数据和有限元计算结果的比较如图9所示，无论采用显式动力计算还是隐式静力计算，基于线弹性模型得到的位移结果都明显小于实测数据.

图9 线弹性2D模型的位移曲线Fig.9 Displacement curves of linear elastic 2D model

3.1.1线弹性模型的误差分析由图9不难看出，线弹性模型在反映PMMA变形方面存在较大误差，主要原因有：① 加载过程中，随着压力增大，蠕变效应产生的位移越来越重要，线弹性模型由于忽略加载过程中的蠕变位移，将产生显著的误差；② 在保载阶段，由于载荷维持不变，线弹性模型的位移保持不变，但实验中PMMA会继续发生蠕变变形.

3.1.2显式动力学计算的准确性和稳定性采用Abaqus的显式动力学模块进行计算，分别考察计算结果的准确性和稳定性.

比较图9中显式动力学、隐式静力学的位移曲线可知，两者误差很小.据此认为，使用质量缩放系数加速的显式动力学分析可得到准确的计算结果.

一般认为，若模型动能与内能之比小于5%则计算结果稳定[11].取质量缩放系数为1×103进行有限元计算，得到模型能量(W)如图10所示.计算过程中，模型动能的数量级从10-5下降到10-9，除了在加载开始极短的时间内，模型的动能与内能之比均远小于5%，可认为显式动力计算结果稳定.

图10 模型的内能和动能Fig.10 Internal energy and kinetic energy for model

3.2 黏弹性模型

3.2.1黏弹性参数 Luo等[12]对PMMA进行了一系列不同温度和应力水平下的较长时间的蠕变实验，并提出了时间-温度-压力等效性理论.本文通过分析蠕变实验数据，推导出基于Abaqus线性黏弹性的Prony级数参数，据此开展有限元动态分析.

图11 Prony级数曲线与蠕变实验数据比较Fig.11 Prony series curves compared with creep test data

大量蠕变实验结果表明：在较低应力水平和较低温度下，松弛效应将放缓，即松弛时间参数将变大.由于文献[12]中蠕变实验所采用的应力普遍低于观察窗中主要的分布应力，且其实验温度也低于观察窗实验的环境温度，故需要适当调整黏弹性参数，尤其是要缩短松弛时间参数，以得到与观察窗实验数据较吻合的结果.

依据式(9)，在恒定泊松比下，弹性模量与切变模量、体积模量呈现等比例关系，所以弹性模量的Prony级数参数可以作为Abaqus中切变模量和体积模量的输入参数.运用等效时长法处理后，得到如表4所示的Prony级数参数.据此，瞬时弹性模量E(t)满足

(14)

式中：E0=2.76 GPa为初始时刻弹性模量.

表4 用于有限元分析的Prony级数参数Tab.4 Prony series parameters for FEM analysis

3.2.2计算结果分析图12为不同质量缩放系数(ζ)下2D和3D模型的计算结果与实验数据的对比.可见，其他条件相同时，2D与3D模型的位移相当接近.总的来说：在加载阶段，各个模型与实验结果吻合良好；而在保载阶段，ζ= 1×104时的计算结果与实验值吻合得更好.

图12 黏弹性模型显式动力计算位移Fig.12 Displacement of viscoelastic model using explicit dynamic FEM

3D模型相比2D模型位移偏大的可能原因：① 2D 模型的窗座采用刚体单元，而3D模型的窗座则用实体单元来模拟，窗座变形会导致位移偏大；② 在显式动力学计算中，时间差分会导致计算误差，3D模型的单元数量远多于2D模型，这种误差会更加明显，表现为位移值偏大；③ 对于一个轴对称的实体模型，即使受到轴对称的外载荷和约束，其在不同轴截面上的有限元计算结果仍会略有差异，显式动力学计算尤其如此，在黏弹性效应更明显的保载期间，截面差异会加大这种计算误差.

从安全性角度出发，考虑到PMMA材料非均匀等变动因素，建议在模拟计算时预留位移冗余.因此，3D模型的计算结果虽与实验结果有一定偏差，但在观察窗设计中，仍然具有较高的参考价值.