球形封头开孔接管结构屈曲特性研究

董乃逸，贺小华

(南京工业大学 机械与动力工程学院，南京 211816)

0 引言

球壳封头结构简单，受力均匀，广泛应用于高压容器、深海耐压潜水器等承压结构[1]。因工艺需要，耐压球壳需开设各类接管结构，开孔接管破坏了球壳结构的整体连续性，在外压作用下，导致球壳承压能力下降，易产生屈曲失稳破坏。对球形封头开孔接管结构屈曲特性及临界失稳压力进行研究意义重大。

对于外压球壳的屈曲行为研究由来已久。早期针对理想球壳，ZOELLY[2]基于小变形理论，推导了外压完整球壳的临界失稳压力公式。由于Zoelly公式计算得到的结果与试验值有较大差距，KARMAN[3]基于外压球壳的轴对称失稳问题，提出了大变形理论公式以计算外压球壳的临界失稳压力，结果与试验值相对较为接近。20世纪60年代，KRENZKE等[4]对200多个具有不同缺陷的耐压球壳进行试验，提出了临界失稳压力的计算公式。目前，国内外标准规范：GB/T 150—2011《压力容器》、ASME Ⅷ-1和ASME Ⅷ-2，以及EN 13445-3:2014分别基于上述小变形理论及修正公式，提出了临界失稳压力的计算公式。文献[5]分析描述了球壳在均匀外压下的前屈曲状态，使用有限元计算球壳的临界失稳压力和屈曲模式，显示了球壳的后屈曲行为；文献[6]研究了外压半球封头的弹塑性屈曲，提出了预测半球封头临界失稳压力的解析公式，并与数值结果进行对比。

由于球形封头开孔接管结构在其开孔处的应力状态较为复杂，难以通过解析方法求解球形封头开孔接管结构临界失稳压力。文献[7]基于外压球壳开孔区域的轴对称性和等面积补强法，对大深度潜水器载人球壳开孔加强形式进行几何变换，完善了球壳开孔加强理论计算方法，并对球壳进行外压试验，验证了计算方法的有效性和适用性。文献[8]运用有限元软件ANSYS，分析了侧壁加强、垫板加强、肘板加强这三种耐压球壳开孔加强方式的极限强度并进行对比，结果表明侧壁加强的极限强度最高；同时考虑到结构性能和工艺性方面，提出了梯形加强这一新型的开孔加强形式。文献[9]运用有限元软件ANSYS对潜水器开孔耐压壳半球封头承载能力进行非线性分析，讨论了开孔个数、孔径大小、围壁厚度3个参数对耐压壳半球封头承载力的影响，结果表明，开孔越多、孔径越大、围壁厚度越小都会使球形封头承载能力下降。

本文系统讨论开孔率di/Di、接管与封头厚度比δet/δe及径厚比Di/δe这3个无量纲结构参数对球形封头开孔接管结构临界失稳压力Pcr及屈曲特性的影响规律，为外压球形封头开孔接管结构稳定性研究及工程设计提供参考。

1 有限元模拟计算

1.1 材料性能参数

球形封头材料采用Q345R，接管材料采用16Mn，常温下的材料性能参数如表1所示。

表1 常温下材料性能参数

1.2 有限元计算方案

有限元分析结构为球形封头径向开孔接管结构，如图1所示。球形封头厚度取定值δe=10 mm，取筒体直边段长度L1=25 mm，接管外伸长度L2=200 mm，通过改变开孔率di/Di、接管与封头厚度比δet/δe和封头径厚比Di/δe的值进行有限元模拟计算。di/Di分别取值0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8；δet/δe分别取值0.2,0.4,0.6,0.8,1.0,1.25,1.5,1.75,2.0；Di/δe分别取值50,100,150,200。有限元分析模型共计288组。

1.3 材料本构关系

文献[10]指出本构关系对壳体结构的临界失稳压力影响甚少，因此,本文选取双线性材料本构模型进行有限元分析，塑性段斜率取E/100，其本构关系如图2所示。

图1 球形封头径向开孔接管结构

图2 双线性材料本构模型

1.4 单元类型及边界条件

分析结构选择Shell 181壳单元，该单元计算精度已在文献[10-11]中证实。在满足网格无关性要求下进行网格划分，分析结构的边界条件为封头直边段施加固支约束，开孔接管端面施加环向约束以保证接管端面的圆度，封头与接管均受外压作用。球形封头开孔接管有限元模型见图3。

图3 封头开孔接管有限元模型

2 有限元模型验证

2.1 完整球壳临界失稳压力值比较

对球壳屈曲行为进行的研究历程是从最初的不考虑缺陷，到后来发现初始缺陷对其屈曲耐压特性有巨大的影响。实际工程中半球封头是冲压成型的，由于模具和工艺不可避免地会造成一定的缺陷，从而使半球封头的实际形状和完美的半球封头存在一定的差距。本文选取文献[12]中的4个试验模型，模型参数分别为球壳直径D=117.54,117.98,117.54,117.5 mm，计算厚度分别为0.422,0.423,0.406,0.415 mm，弹性模量E=1.93×105MPa，泊松比μ=0.28，屈服强度σs=205 MPa。分别从理论公式、试验结果和有限元数值模拟三方面，对外压球壳的临界失稳压力进行计算对比。

小变形理论公式是基于完美几何形状且无缺陷的球壳，其临界失稳压力的计算公式[2]：

Pcr=1.21E(δe/R)2

(1)

式中，E为材料弹性模量,MPa；δe为球壳计算厚度,mm；R为球壳半径,mm。

大变形理论公式基于结构的非线性因素影响，其临界失稳压力的计算公式[3]：

Pcr=0.36516E(δe/R)2

(2)

本文有限元数值模拟时施加的缺陷方法有两种，一种是采用“一致缺陷模态法”[13]设置初始几何缺陷(一致缺陷模态法是基于最低阶屈曲模态，从而模拟结构的初始缺陷分布),通过ANSYS中的UPGEOM命令里的FACTOR值对模型统一施加10%的几何缺陷，然后进行非线性求解；另一种是施加局部厚度缺陷。GB/T 150.4—2011《压力容器 第4部分：制造、检验和验收》中6.4.2规定球形封头的最大形状偏差外凸不超过1.25%Di，内凹不超过0.625%Di。文献[14-15]表明内凹0.625%Di下临界失稳压力低于外凸1.25%Di下临界失稳压力，因此本文模型采用形状偏差内凹0.625%Di的球形封头模型，并在封头顶部最大壁厚减薄位置施加单个缺陷厚度减薄15%进行非线性分析。最终求解提取出位移最大点的载荷-位移图，并根据零曲率准则得到临界失稳压力值Pcr。

将文献[12]试验结果、公式(1)解、公式(2)解及本文有限元特征值、有限元非线性解1(一致缺陷模态法)以及有限元非线性解2(施加局部厚度缺陷)列于表2，并将其数据作成图4，以便于直观比较。

表2 球形封头临界失稳压力值Pcr比较

图4 球形封头临界失稳压力值Pcr比较

由表2和图4可以看出，有限元特征值与小变形公式(1)解结果较为接近，由于这两种方法均未考虑几何非线性以及缺陷的存在，导致其结果远大于其他解；大变形公式(2)解略大于试验值，计算结果偏于保守；有限元非线性解1及有限元非线性解2两者均考虑了结构的几何大变形和初始几何缺陷，因此其值均与试验结果较为接近。考虑到一致缺陷模态法结果的有效性及操作的方便性，后续有限元模型均采用此法设置初始几何缺陷。

2.2 开孔球壳临界失稳压力结果比较

选用文献[16]中按1∶10比例自制的不锈钢缩比试验模型，模型结构如图5(a)所示。该模型参数为：球壳外半径R0=100 mm，球壳厚度t=6.3 mm，接管内径d1=50 mm，接管外径d=77 mm，接管厚度t1=13.5 mm，内伸高度h1=32.2 mm，外伸高度h2=32.2 mm，弹性模量E=1.95×105MPa，泊松比μ=0.3，屈服强度σs=300 MPa。试验测得的缩比模型临界失稳压力为20.5 MPa，采用一致缺陷模态法，施加10%的初始缺陷进行有限元非线性求解，有限元模型如图5(b)所示，得到临界失稳压力为19.223 MPa，与试验值比较，误差为6.64%。两者误差基本在工程允许范围内，由此也说明本文非线性屈曲分析方法的有效性。

(a)开孔球壳围壁结构示意

(b)开孔球壳围壁有限元模型

3 球形封头开孔接管结构有限元非线性屈曲分析

3.1 临界失稳压力削弱系数K

为了分析讨论开孔接管结构对球形封头临界失稳压力的影响，通过引入削弱系数K来反映开孔接管的作用，K值为球形封头开孔接管结构的临界失稳压力Pcr与球形封头未开孔结构的临界失稳压力Pcrs(有限元非线性模拟值)之比，即K=Pcr/Pcrs。

3.1.1 接管与封头厚度比δet/δe与K值关系

图6示出了不同径厚比Di/δe及不同开孔率di/Di下，临界失稳压力削弱系数K与δet/δe的分布图。可以看出，不同径厚比及不同开孔率下，K值与δet/δe变化规律基本相同；同一δet/δe下，K值随di/Di的增大而减小；当δet/δe较小时，K值明显偏低，随着δet/δe增大，K值显著增大，当δet/δe≥1.0后，K值随δet/δe增大变化趋于平缓。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

3.1.2 开孔率di/Di与K值关系

图7示出了不同径厚比Di/δe及不同接管与封头厚度比δet/δe下，临界失稳压力削弱系数K与开孔率di/Di的分布图。可以看出，Di/δe与δet/δe一定时，K值随着di/Di增大而减小；当δet/δe较小时，K值随di/Di增大显著降低；而当δet/δe较大时，K值随di/Di增大、降低幅度较小；当δet/δe≥1.0时，同一di/Di下，δet/δe对临界失稳压力削弱系数K的影响较小。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

3.2 计算模型失稳模态分析

由图6可以看出，临界失稳压力削弱系数K与接管与封头厚度比δet/δe变化有个转折点，当δet/δe≥1.0时，K值随δet/δe增大变化不显著。同样由图7可以看出，当δet/δe≥1.0时，同一开孔率下，不同计算模型临界失稳压力削弱系数K值较为接近。为进一步讨论球形封头开孔接管结构屈曲特性，图8示出径厚比Di/δe=150时部分计算模型的失稳模态图。

图8 球形封头中心开孔接管结构失稳模态图(Di/δe=150)

由图8(a)可以看出，较小开孔率(di/Di=0.4)，当接管较薄(δet/δe=0.2)时，计算模型失稳位置位于接管；当δet/δe增大到0.6时，由图8(b)可知，计算模型失稳位置由接管转移到封头-接管连接处；当δet/δe增至2.0时，由图8(c)可知，计算模型失稳位置位于封头处。

由图8(d)～(f)可看出，当接管较薄(δet/δe=0.2)时，计算模型失稳位置位于接管，随着接管壁厚增大，计算模型失稳位置逐渐由接管转移到封头。由此，解释了图6中临界失稳压力削弱系数K与δet/δe分布图中转折点的由来。当接管与封头厚度比δet/δe低于图6中转折点时，接管较薄，承压能力较弱，计算模型失稳位置位于接管，此时分析结构临界失稳压力较小；当δet/δe增大到一定值时，接管承压能力超过封头，计算模型失稳位置转到封头，此时分析结构临界失稳压力较大，进一步由图6可以看出，随着开孔率di/Di增大，失稳位置由接管转到封头所需接管与封头厚度比越大。本文计算参数范围内，当δet/δe≥1.0时，计算模型失稳位置均位于封头上。

3.3 外压球形封头开孔接管结构临界接管与封头厚度比(δet/δe)c

由图6可以看出，临界失稳压力削弱系数K与δet/δe变化存在转折点。当δet/δe低于转折点对应值时，计算模型失稳位置位于接管；而当δet/δe高于转折点对应值时，计算模型失稳位置位于封头。本文将失稳位置由接管转向封头的最小接管与封头厚度比称为临界接管与封头厚度比 (δet/δe)c。表3列出本文计算参数范围下的临界接管与封头厚度比(δet/δe)c。

表3 临界接管与封头厚度比(δet/δe)c

由表3可以看出，临界接管与封头厚度比(δet/δe)c是3组定值(0.6,0.8,1.0)，其与径厚比及开孔率之间的关联式如下：

(3)

表3中列出的临界接管与封头厚度比(δet/δe)c基本对应于图6中曲线的转折点。由图6可以看出，当δet/δe<(δet/δe)c时，临界失稳压力削弱系数K明显较小，结构承载能力较弱；当δet/δe≥(δet/δe)c时，K值相对较大，且K值随δet/δe变化不显著。实际工程设计中，建议参考表3中给出的临界接管与封头厚度比，取δet/δe=(δet/δe)c，可以在满足结构轻量化设计基础上，获得较优的结构承载能力。

3.4 临界失稳压力削弱系数K值与结构参数的经验关系式

基于临界失稳压力削弱系数K随Di/δe，di/Di，δet/δe的变化曲线，进一步拟合出K值与3个结构参数的经验关系式。由于所研究参数δet/δe变化范围(0.2～2.0)较大，且K值随δet/δe的变化可分为两个阶段。为了保证拟合精度，采用分段拟合方法。当δet/δe位于[0.2,0.6]区间内，采用以下关系式进行拟合：

(4)

式中，a1～a7为待定系数。

基于有限元模拟所得数据点，最终拟合结果为：

(5)

当δet/δe位于(0.6,2.0]区间内，基于有限元模拟所得数据点，最终拟合结果为：

(6)

为了验证所得关系式的可靠性，进一步给出当开孔率di/Di=0.2,0.5,0.8时，利用以上经验关系式所得拟合K值与实际有限元K值数据点对比图，如图9所示。图9中补充计算了δet/δe=0.7时相应有限元模型的临界压力削弱系数K，并与第二区间段公式拟合值进行比较。由图9可以看出，不同区间内拟合所得K值与实际有限元K值两者较为接近，且拟合曲线能够反映出K值随不同结构参数的变化趋势。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

计算表明,除了δet/δe=0.2时个别点误差较大，其他点公式拟合值和有限元模拟值误差较小，所有分析计算点平均误差小于10%。

4 结论

分析讨论3个结构参数：径厚比Di/δe、开孔率di/Di、接管与封头厚度比δet/δe对外压球形封头开孔接管结构临界失稳压力Pcr的影响。研究结果如下。

(1)球形封头临界失稳压力削弱系数K值随着开孔率di/Di增大而减小。

(2)接管与封头厚度比较小时，临界失稳压力削弱系数K值越小，计算模型失稳位置位于接管，当接管与封头厚度比增至一定值时，K值显著增大，计算模型失稳位置位于封头。

(3)开孔率越大，失稳位置由接管转到封头所需接管与封头厚度比越大。本文计算参数范围内，当δet/δe≥1.0时，计算模型失稳位置均位于封头上。

(4)提出了临界接管与封头厚度比(δet/δe)c，实际工程设计中，当接管与封头厚度比取临界接管与封头厚度比时，可以在满足结构轻量化设计基础上，获得较优的结构承载能力。

(5)建立了临界失稳压力削弱系数K值与不同结构参数的经验关系式，对比发现所得关系式预测精度较高。