薄壁圆筒壳初始几何缺陷不确定性量化的极大熵方法

李建宇，佘昌忠，张丽丽，杨 坤

(1.天津科技大学 机械工程学院，天津市轻工与食品工程机械装备集成设计与在线监控重点实验室，天津 300222；2.天津职业技术师范大学 理学院，天津 300222)

1 引 言

薄壁圆筒壳由于其优越的承载性能和简易的制造工艺在航空航天、军事导弹和化工液罐存储等现代工程中广泛应用。轴压条件下其结构形式非常简单，但表现出的力学特性却十分复杂，存在非弹性效应、大变形和屈曲等因素，从而导致该结构的轴压承载力总是由其稳定性控制。在关于轴压薄壁圆筒壳稳定性研究中，一个长期关注的课题就是屈曲载荷的实验结果与理论解之间存在较大的差异，关于其机理的研究促进了后屈曲理论[1]和缺陷敏感性理论[2]的发展。已有研究表明，初始缺陷的形状和大小均对轴压薄壁圆筒壳的屈曲载荷有明显影响，而该结构产品在其全生命周期内形成的真实几何缺陷又具有不确定性，如何量化这种不确定性并在设计环节给出既安全又经济的屈曲载荷目前仍是颇具挑战的研究课题。

关于薄壁圆筒壳初始几何缺陷量化表征的研究大致可分为确定性方法和概率方法[3]。确定性方法的思路一般是通过人为假设或某种搜索手段确定最不利缺陷形式，并据此获得屈曲载荷的下限值，如实测缺陷法[3]、特征模态缺陷法[4]、单点窝陷扰动法[5]以及各类缺陷组合和最优搜索类方法[6]。确定性方法由于其简洁性和便利性而得到广泛运用。但其不足是关于屈曲载荷不确定性量化的信息量太少而使得对承载力下限值的预测偏保守。概率方法针对具有不确定性特征的初始缺陷在统计意义下呈现的规律性，利用随机变量或随机场对初始缺陷进行建模，并通过随机分析方法给出屈曲载荷的统计预测结果[1]。在20世纪60年代，Budiansky等[7]提出应该对初始几何缺陷进行随机描述，以实现对缺陷高度敏感结构屈曲载荷的更可靠预测。利用概率方法量化表征缺陷的难点在于随机变量的选定及其概率密度函数的确定。Elishakoff等[8]最早将实测缺陷的傅立叶系数直接随机化为高斯型随机变量，并通过蒙特卡洛方法获得屈曲载荷的随机分布结果。近年来，随机场理论逐渐用于缺陷的建模表征，如文献[9]基于实测数据的功率谱分析，构造初始几何缺陷随机场的谱表征方法，文献[10]基于实测数据并假定随机场自相关结构，构造缺陷随机场的KL(Karhunen-Loève)分解方法。利用随机场理论对缺陷建模在数学上显得更自然，但目前的不足是，由于实测数据缺乏，已有研究通常对缺陷的随机场模型施加各类假设，如缺陷随机场的平稳性(Homogeneity)、高斯性和各态历经性等[11]。事实上，通过对已有实验数据的统计检验，发现真实缺陷大多不具备这些性质，这导致预测的屈曲载荷往往有偏[12]。

本文提出一种基于实测缺陷数据和极大熵建模的初始几何缺陷随机建模方法。基于实测缺陷数据的统计结果并引入极大熵原理对统计变量进行随机建模，给出已知信息条件下对初始几何缺陷随机分布的一个最无偏估计。该方法不仅无需事先假定随机场自相关结构，而且也放松了对随机场的平稳性、高斯性和各态历经性等假设。

2 初始几何缺陷的极大熵建模

薄壁圆筒壳的几何缺陷表现为实际与理想之间的形状偏差，包括圆柱度误差和表面局部凹陷等多种形式[13]。鉴于缺陷形式及分布具有随机性，学者们提出利用随机场理论对其进行建模表征。如图1所示，左侧是一高为L，截面半径为R，筒厚为t的薄壁圆筒壳结构。右侧矩形面为该结构沿其任意母线展平所得，记s=(x,y)为矩形面各点坐标，其中x为环向坐标(0≤x≤ 2πR)，y为轴向坐标(0≤y≤L)。

图1 薄壁圆筒壳及其展开示意图

将初始几何缺陷表征为薄壁圆筒壳与理想圆筒壳相比在半径方向的偏差，记为a(s,ω),且a(s,ω)∈Θ2(Ω,F,P)为一随机场，其中(Ω,F,P)为概率三元组，ω∈Ω为样本，F为事件域，P为概率测度，Θ2(Ω,F,P)为二阶范数有限的可测泛函空间[14]。

2.1 缺陷随机场的谱分解表征

假设已经获得M个薄壁圆筒壳径向初始几何缺陷的测量结果，令a1,a2,…,aM表示对缺陷a(ω)=[a(s1,ω),…,a(sN,ω)]的M次独立观测，其中N为每个样本筒的观测点数。基于观测数据，可获得缺陷随机场的均值和协方差矩阵等统计特征量，即

(1)

(2)

(3)

(4)

式中 〈·，·〉为在RN上的内积。由KL展开理论，缺陷随机场a(s,ω)还可降维表示为

(5)

2.2 KL随机向量的极大熵建模

已有的实测数据表明，薄壁圆筒壳结构的初始几何缺陷在统计上并不一定具有高斯性，这导致其KL随机向量ξ也不具有高斯性,其概率分布形式也往往是事先未知的。为此，本文利用极大熵方法给出其概率密度函数的一个估计。极大熵原理的基本思想是，在利用已知信息对随机向量的概率分布进行推断时，使得概率熵函数取得最大的概率分布是最合理的[15]。极大熵原理给出的概率分布认为是已知信息下概率分布的最无偏估计。为简化计算，假设KL随机向量ξ各分量之间是相互独立的。对于每一个分量ξj(j=1,…,r),其概率密度函数p(ξj)的极大熵估计为凸优化问题(6)的解。

(6)

式中μj 0=1,μj q为分量ξj的q阶统计中心矩，q=1，…，K。解凸优化问题(6)得

(7)

式中β={β0,…,βK}为优化问题(6)的拉格朗日乘子，是非线性方程组(8)的解，

(8)

在获得随机向量ξ的概率密度函数后，还需对其进行随机抽样，以便进行后续的随机屈曲载荷计算，本文利用马尔科夫链蒙特卡洛(MCMC)抽样[16]。

3 随机屈曲载荷分析

在获得缺陷随机场模型后，利用不确定性传播算法进行随机屈曲分析，最后统计各样本筒屈曲分析结果，并获取统计规律。

在工程运用中，折减因子KDF(Knockdown Factor)常用来度量薄壁圆筒壳轴压屈曲载荷的缺陷敏感性。令Pc r,perf为理想筒的临界载荷，Pc r,i mp为含缺陷筒的屈曲载荷，定义归一化因子αc r为

αc r=Pc r,perf/Pc r,i mp

(9)

在随机屈曲分析的基础上，定义随机变量Ac r，其概率分布函数FAc r(αc r)为

(10)

概率密度函数fAc r(αc r)可利用2.2节的极大熵方法获得。定义可靠度R(α)为

(11)

进一步给出基于可靠度的KDF估算方法，即可靠度为R(α)的α的估计值为KDF。

4 数值算例

表1 模型参数

图2 轴向加载圆筒壳的有限元模型

图3 圆筒壳的屈曲载荷收敛性分析

进一步利用NASA SP-8007法与特征模态缺陷法等经典含缺陷屈曲分析方法给出KDF的估计值分别为0.241和0.239。

4.1 数值实验1

为了测试极大熵方法的正确性，假定真实缺陷可表征为一个平稳高斯随机场，随机场的均值和协方差函数参见文献[14]。分别抽取5，10和50个样本作为初始几何缺陷数据，利用极大熵方法建模，并进行随机屈曲分析，进一步获得αc r的概率密度函数，考察KDF的估计值与理论值的关系。

表2对比了置信度为0.95时不同样本数对应KDF的估计值与理论值，可知随着样本筒数目的增加，KDF的估计值逐渐收敛于理论值。同时由图4可知αc r概率密度函数的估计和相关统计置信度都受观察的薄壁圆筒壳样本数的影响，并且随着更多的观察，两者的准确性都会提高，此时极大熵方法能给出更无偏的αc r概率密度函数估计。

表2 屈曲载荷预测结果对比

图4 不同样本的估计概率密度及附近95%的置信区间与理论值的对比

4.2 数值实验2

为验证本文方法对薄壁圆筒壳实际屈曲问题的有效性，利用文献[17]中基于实际测量的初始几何缺陷的傅里叶级数表示形式分别获取7个薄壁圆筒壳的缺陷数据。将这7个样本筒的缺陷数据视为对同一批次生产的薄壁圆筒壳缺陷的7次独立观测。

按照式(2)求出协方差矩阵C，并图形化显示为图5，发现协方差矩阵的数值不仅取决于两个节点的相对距离，而且与节点位置有关，说明缺陷随机场是非平稳的。

图5 缺陷随机场的协方差矩阵

对缺陷随机场进行KL分解，图6给出其协方差矩阵特征值λj的变化趋势，可知特征值自第6阶开始接近0。在确保包含实测缺陷数据主要特征信息的基础上进行降噪去冗余，截取前6阶特征值即可使累计方差贡献率达到99%。

图6 协方差矩阵的特征值

将截取的前6阶特征值及其对应的特征向量分别代入式(4)获取6维随机向量ξ，并利用极大熵方法分别计算该6维随机向量ξ的概率密度函数PDF，如图7所示。可以看出，随机向量ξ显然是非高斯型的。进一步利用 Shapro -Wilk 方法进行高斯性检验，由表3可知ξ2和ξ4两项显著拒绝高斯性假设(其显著性概率P≤0.05)。

利用MCMC法对随机向量ξ的概率密度函数抽取2000个样本，代入式(5)获取相应几何缺陷数据进行随机屈曲分析。

秦铁崖并不着急，一边砰砰砰击打对方，一边用跳步和滑步避开对方攻击，边打边思索，寻找对策。他的脑海里，老是想起黄河漂流所用的羊皮筏子。那些羊皮囊，经过河水冲激，经过岁月打磨，已变成深黄色乃至深棕色，坚韧异常。而现在，张万邦的躯壳，就是那皮囊，怎么打也不碍事，他的五脏六腑，被胀鼓鼓的皮囊保护得好好的，不仅隔着一层皮壳，而且隔着一层气。因为有这层气，五脏六腑根本不会受伤，想震荡它们一下，都很困难。

图7 随机向量ξ的频率分布直方图及极大熵原理近似的概率分布

表3 随机向量ξ的高斯性检验

表4列出了实验数据、相关结构法以及本文结果。可知文献[13，18]中事先假定的随机场自相关结构进行数值模拟所得的ac r均值都明显偏大，这往往会导致对薄壁圆筒壳承载能力偏危险的预测。而本文计算所得ac r的变异系数相较实验结果偏大，但其均值却与实验结果基本一致。

进一步分析本文预测结果分散性偏大的原因。一方面，极大熵方法给出的是已知信息条件下概率分布的最无偏估计，而已知信息依赖于已知的数据量及其统计矩信息，越少的数据和越低的统计矩信息会导致相应的极大熵分布越接近均匀分布。本文只用到了7个实测样本，4阶统计矩，因而所给出的缺陷随机分布模型分散性偏大；另一方面，实验结果给出的屈曲载荷不仅受初始几何缺陷的影响，还受材料不均匀性、非完美边界条件、荷载对中性偏差及结构尺寸偏差等非传统缺陷的影响，而本文在进行屈曲分析时并未考虑这些因素，这也导致屈曲载荷计算结果与实验结果的差异。

图8展示了实验结果与再次利用极大熵方法获取的含缺陷薄壁圆筒壳αc r的概率密度函数。

图8 各种方法估计薄壁圆筒壳αc r的概率密度对比

进一步给出基于可靠度的折减因子估算结果，从图9可以看出，当可靠度R(α)为0.95时对应的KDF为0.417，相较NASA SP-8007仅考虑径厚比R/t对屈曲承载能力下限的预测方法计算的 0.241 更为经济。

图9 αc r的可靠性函数

5 结 论

本文从随机场的角度较为系统地研究了轴压薄壁圆筒壳初始缺陷的建模方法。首先通过对实验数据进行统计分析，获得随机场模型的均值和协方差矩阵等统计特征；其次，基于KL展开法将缺陷随机场建模转化为KL随机向量的建模问题；然后，基于KL随机向量非高斯性的事实，提出利用极大熵方法构建其概率密度函数；最后，结合MCMC抽样和确定性屈曲分析方法，给出轴压薄壁圆筒壳屈曲分析的不确定性传播计算，获得屈曲载荷的随机分析结果，并给出基于可靠度指标的折减因子估算结果。

本文方法有两个特点，一是全面利用了实测缺陷数据的一阶和二阶统计信息，这使得构建的模型较为全面地涵盖了实际缺陷信息；二是从极大熵原理的角度进行随机建模，避免了文献中关于缺陷随机场高斯性、平稳性和各态历经性等假设，使得所构建的概率模型是已知信息条件下的一个最无偏概率分布估计。这些特征的具备，使得本文预测的屈曲载荷结果较好地覆盖了实验结果，并提高了折减因子估算的可靠度。