CN110991112B - 一种基于响应面法的随机缺陷模态叠加方法 - Google Patents

Abstract

本发明提供的一种基于响应面法的随机缺陷模态叠加方法，步骤为：利用多项式响应面法，将临界荷载拟合表达式表达为多项式的形式，根据临界荷载拟合表达式计算确定表达式待定系数，并由计算得到的待定系数确定临界荷载拟合表达式；确定临界荷载；计算结构临界荷载的失效概率。该方法通过建立结构缺陷模态组合系数的概率模型，利用多项式拟合出随机变量,对具有随机缺陷的经典网壳和网状壳进行稳定承载力分析，弥补了随机缺陷法及改进随机缺陷法样本计算量大的缺点，提高了随机缺陷模态叠加的效率，给出了设计临界荷载的可靠性问题。

Description

一种基于响应面法的随机缺陷模态叠加方法

技术领域

本发明属于工程结构计算分析领域，特别涉及一种基于响应面法的随机缺陷模态叠加方法。

背景技术

深单层网壳结构的典型破坏形态为失稳破坏，其失稳模式和稳定承载力分析方法的研究十分重要。目前，基于非线性有限元理论的荷载-位移全过程分析法是研究网壳结构非线性平衡路径的主要分析方法。该方法的关键是平衡路径的跟踪求解，目前的主要方法包括人工弹簧法、位移控制法、弧长控制法及自动增量求解技术。初始缺陷对单层网壳结构的失稳模式和稳定承载力有显著影响。

对有初始缺陷网壳结构的稳定性分析，除基于连续化方法的拟壳法外，目前应用较多的还有基于离散方法的解析法和数值计算法。在解析法研究方面，Koiter提出了缺陷敏感性分析的渐进理论，Tompson、Budiansky以及Hutchinson等发展了这一理论。渐进法以分支点的微小邻域作为研究对象，难以直接应用于复杂结构，因而数值计算方法开始发展起来，包括确定性方法和随机有限元方法。其中，忽略了缺陷随机性的确定性方法又包括优化方法、临界缺陷模态法、一致缺陷模态法及其改进方法。优化方法不具有通用性，即需针对具体问题进行求解。临界缺陷模态法认为非完善结构的平衡路径可看作完善结构在缺陷作用下发生的微小扰动，因而难以分析位移较大的复杂结构。

一致缺陷模态法认为最低阶屈曲模态是结构屈曲时的位移倾向，与结构屈曲模态相同的初始缺陷对结构产生不利影响，但是，目前尚没有相关研究成果或理论可以证明此说法，最低阶屈曲模态也可能不是最不利缺陷模态。关于初始缺陷随机性的研究相对较晚。赵惠麟等运用Monte Carlo法对随机稳定承载力进行了研究，假设结构的初始缺陷服从正态分布，得到承载力的特征值，进而得到统计意义上的结构稳定承载力。黄斌等采用随机缺陷模态法考虑节点缺陷的随机性，计算结果精度较高，但对多自由度复杂结构，随机变量为3倍节点个数，样本计算量大。

魏德敏基于概率可靠度理论对随机缺陷模态法进行了改进，认为样本的抽样数量应不少于90组。蔡健等提出n阶特征屈曲模态法，且认为使用N阶特征缺陷模态法需计算前20阶特征屈曲模态。有研究指出缺陷模态可看作是若干屈曲模态的耦合，对基于某一阶屈曲模态的缺陷网壳结构进行一致缺陷模态分析，定性得出缺陷的前4阶屈曲模态耦合系数，并未对缺陷模态耦合的随机性进行理论论证及数值分析。

进一步地，Luca Bruno等基于特征缺陷模态法对子午线网壳前20阶屈曲模态进行两两组合获得结构的初始缺陷，得出最不利的稳定荷载是由屈曲模态组合产生的，同时，最不利的缺陷模式受到组合系数不确定的影响。对此，有研究者将该方法推进到随机领域，提出了基于蒙特卡洛抽样法的随机模态叠加法，算例和试验数据均证明该方法在当时具有先进性。徐军改变了随机缺陷模态叠加法中的随机变量分布形式，对网壳随机初始缺陷下的稳定性做了分析，并试探性地做了基于变形限值的结构可靠度分析。

由上可知，现有渐进法只适用于简单结构，难以应用于在复杂工程中。而确定性方法由于忽略了缺陷的随机性，难以得出准确的结果，随机缺陷模态法及改进方法通常是基于Monte Carlo方法，样本计算量较大，工程应用受限。基于蒙塔卡罗方法的随机模态叠加法因受制于蒙特卡洛抽样的固有特性，仍需保证一定量的样本，以获得稳定和收敛的高精度解。而采用精细化程序编制的改进的随机模态叠加法在应用于大型复杂结构计算中时操作性欠佳。因此，如何保持该方法精度和操作性的同时，进一步提高计算速度，将是一个重要的研究方向。

发明内容

技术问题：为了解决现有技术中如何保持该方法精度和操作性的同时，进一步提高计算速度，本发明提供了一种基于响应面法的随机缺陷模态叠加方法。

技术方案：本发明提供的一种基于响应面法的随机缺陷模态叠加方法，步骤为：

(1)建立结构缺陷模态组合系数的概率模型，测定验算点相关数据；利用多项式响应面法，将临界荷载拟合表达式

表达为多项式的形式，根据临界荷载拟合表达式

计算确定表达式待定系数，并由计算得到的待定系数确定临界荷载拟合表达式

(2)确定临界荷载：临界荷载采用以下方法确定：

其中，

为临界荷载p_cr(r)的均值，

为临界载荷p_cr(r)的方差，p^*为最不利临界载荷；

(3)计算结构临界荷载的失效概率：随机缺陷下结构失效概率采用以下方法确定：

P_f＝P(p_cr(r)-5p₀≤0)      (2-14)

其中，P_f为随机缺陷下结构失效概率，p₀为网壳结构无缺陷时承载力设计值，p_cr为响应值-临界荷载。

其中，步骤(1)中，对大型结构来说，验算点的选择方式为：已知各基本变量的分布形式及分布参数，以均值点为中心展开,验算点的选择范围是

其中，

是随机变量r_i的均值，μ是均值，

是随机变量的方差，f是常数，根据工程上的3σ原则选择。

其中，步骤(1)中，临界荷载拟合表达式

为二次完全多项式，则

或者临界荷载拟合表达式

为一次多项式，则

式中，a₀，a_i，a_ii，a_ij是待定系数。

其中，步骤(1)中，待定系数的计算方法为：用最小二乘估算法求a_i,a_ij，假定得到k组数据为

(r)_i，其中i＝1,2,…,k，必须大于等于相应的未知参数的个数且为正整数，由如下(2-7)、(2-8)两式可以求出待定系数：

([r]^T[r]){A}＝[r]^T{P_cr}          (2-7)

{A}＝([r]^T[r])^-1[r]^T{P_cr}            (2-8)

式中，

10.{P_cr}＝(p_cr,1 p_cr,2 L p_cr,k)^T       (2-10)

{A}＝(a₀ a₁ L a_n a₁₁ L a_nn a₁₂ L a_n-1n)^T      (2-11)

上述式子中，r_ij为第i组数据的第j随机变量数值，p_cr,i为第i组数据的临界荷载值。

其中，步骤(1)中，具体而言，结构缺陷模态组合系数的概率模型的建立方法，包括以下步骤：

(一)获取线性屈曲的特征值和相应的屈曲模态：

根据结构刚度矩阵[K_T]，求解线性屈曲的特征值λ_i(i＝1,2,…,m)和相应的屈曲模态{U_i}，其中m是参与组合的模态阶数；

(二)获取随机有限元刚度方程：

(2.1)对结构节点进行编号，依次为1，2，L，n，L、n均为正整数；

(2.2)获得假定结构的任一缺陷模式向量{ΔX}'：其中

设定m为模态参与阶数，r₁，r₂，r_i……r_m为参与系数，是独立随机变量，{U_i}为结构第i阶线性屈曲模态；

(2.3)获得幅值为R的缺陷模式向量{ΔX}：对{ΔX}'进行幅值调整，以获得幅值为R的缺陷模式向量{ΔX}；

(2.4)获得节点坐标{X}的切线刚度矩阵[K_T]：先将(2.3)中获得{ΔX}计算到{X}＝{ΔX}+{X₀}中，设定：{X}引入缺陷后结构的节点坐标，{X₀}为3n维完善结构节点坐标向量(n为结构节点数)，{ΔX}为3n维结构节点位置缺陷向量；再采用引入缺陷后结构的节点坐标{X}形成；

(2.5)获得随机有限元刚度方程：采用(2.4)[K_T]获得[K_T]{Δa}＝{Q}-{F}，式中[K_T]为形成当前有缺陷结构节点坐标{X}的切线刚度矩阵，{Δa}、{F}均为随机变量r_i的函数，{Q}是荷载列向量；

(三)获得随机变量输入后的响应值-临界荷载p_cr：

利用弧长法，求解步骤(二)中随机有限元刚度方程，获得结构在随机变量输入后的响应值：临界荷载p_cr，及临界荷载的拟合值

(3.1)若取二次完全多项式，则

(3.2)若不考虑二次交叉项系数，则

(3.3)简单取

一次多项式，则

设定(3.1)-(3.3)中，a₀，a_i，a_ij，a_ii是待定系数；r_i，r_ir_j为参与系数；

为样本结构下临界荷载p_cr的近似值，i，j均为正整数；r为样本结构的参与系数。

其中，步骤(2)中，临界载荷的获得方法为：通过{r_i}进行Box-Behnken矩阵法、CCD法或2ⁿ析因法中任一种抽样方法，抽样次数为k次，结合有限元方程组，获得

为临界荷载p_cr(r)的均值，

为临界载荷p_cr(r)的方差，p^*为最不利临界载荷。

本发明还提供了一种基于响应面法的网壳结构随机缺陷模态叠加方法，方法步骤为：(1)对网壳结构的结构杆件离散为单元并且对其节点进行编号，定义结构杆的材料性能参数、荷载和边界条件，完成数值模型的建立；(2)针对数值模型，参与组合的模态为特征值分析的前n阶，即随机变量数为n，n为正整数，采用CCD配点方法，进行随机缺陷模态叠加法分析，获取所求临界荷载方差、均值以及所消耗的CPU时间。

其中，步骤(2)中所述随机缺陷模态叠加法采用权利要求1至6中任一种随机缺陷模态叠加方法。

有益效果：本发明提供的一种基于响应面法的随机缺陷模态叠加方法，在基于蒙塔卡罗的随机缺陷叠加法的基础上，通过建立结构缺陷模态组合系数的概率模型，采用试验设计方法，利用多项式拟合出随机变量，提出基于响应面法的改进的网壳结构稳定性随机缺陷模态迭加法。该方法对具有随机缺陷的经典网壳和网状壳进行稳定承载力分析，弥补了随机缺陷法及改进随机缺陷法样本计算量大的缺点，且提高了随机缺陷模态叠加的效率，给出了设计临界荷载的可靠性问题，并进行了算例验证。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。

附图说明

图1为本发明结构数据模型及节点表号，其中(a)为结构平面模型(m)，(b)为结构里面模型(m)，其中m为国际单位米。

图2为本发明实施例1中不同方法获得的最不利荷载柱形图。

图3为本发明实施例1中不同配点拟合时最不利载荷柱形图。

图4为本发明实施例2中不同方法获得的最不利载荷柱形图。

图5为本发明实施例2中不同配点拟合时最不利载荷柱形图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

对结构进行节点编号，例如结构选为网壳结构，以下进行具体的说明：

1.1有缺陷网壳结构稳定性分析

假定结构节点位置缺陷向量为{ΔX}，则引入缺陷后结构的节点坐标{X}为{X}＝{ΔX}+{X₀}                     (1-1)

(1-1)式中，{X₀}为3n维完善结构节点坐标向量(n为结构节点数)，{ΔX}为3n维结构节点位置缺陷向量。

则，当前有缺陷结构节点坐标{X}的切线刚度矩阵[K_T]，此时增量平衡方程为

[K_T]{Δa}＝{Q}-{F}                  (1-2)

(1-2)式中，{Δa}位移向量，{Q}外荷载向量，{F}为不平衡力。

采用Newton-Raphson方法结合柱面弧长法，可求解方程(1-2)，并可跟踪每个荷载增量下节点位移的增量，获得结构在整个荷载加载过程中的屈曲路径。

通过求解该方程(1-2)，获取结构在整个加载历史过程中每个增量步的结构切线刚度矩阵[K_T]，计算

Det([K_T])＝0                      (1-3)

式(1-3)第一次成立时，即结构第一次达到临界状态，获取临界荷载

和结构屈曲模态构型{U}。此时

λ_cl(1)＝{q}                       (1-4)

结构的临界变形{U_cr}

{U_cr}₁＝{U_s1}                     (1-5)

1.2稳定分析的随机缺陷模态迭加法

首先，根据结构刚度矩阵[K_T]，求解线性屈曲的特征值λ_i(i＝1,2,…,m)和相应的屈曲模态{U_i}，其中m是参与组合的模态阶数(模态参与阶数)。

由于节点位置缺陷具有随机性，考虑缺陷模式{ΔX}的随机性是，假定任一缺陷模式为

式中，m为模态参与阶数，r₁，r₂，…，r_m为对应阶数模态的参与系数，是独立随机变量，{U_i}为结构第i阶线性屈曲模态。

对{ΔX}'进行幅值调整，以获得幅值为R的缺陷模式向量{ΔX}。

则根据式(1-1)可得有缺陷结构的节点坐标向量

{X}＝{ΔX}+{X₀}                     (2-2)

由于{ΔX}为随机变量r_i的函数，则随机有限元刚度方程为[K_T]{Δa}＝{Q}-{F}  (2-3)

式中，[K_T]、{Δa}、{F}均为随机变量r_i的函数。

利用弧长法，求解非线性有限元，结构在随机变量输入后的响应值-临界荷载p_cr。

现基于多项式响应面法，将

近似表达为多项式的形式。若取二次完全多项式，则

式中，a₀，a_i，a_ii，a_ij是待定系数。若不考虑二次交叉项系数，则表达式为

若简单取

一次多项式，则有

为了得到待定系数，需要选择足够的展开点计算值，从而解方程组式(2-6)求出a_i,a_ij,得到临界荷载的拟合表达式

即临界载荷的近似值，为了得到更精确的解，可以引入一些数值计算的冗余度，用最小二乘估算法求a_i,a_ij，假定得到k组数据为

(r)_i，其中i＝1,2,…,k，必须大于等于相应的未知参数的个数且为正整数，由如下(2-7)、(2-8)两式可以求出待定系数：

([r]^T[r]){A}＝[r]^T{P_cr}                    (2-7)

{A}＝([r]^T[r])^-1[r]^T{P_cr}                   (2-8)

式中，

{P_cr}＝(p_cr,1 p_cr,2 L _pcr,k)^T          (2-10)

{A}＝(a₀ a₁ L a_n a₁₁ L a_nn a₁₂ L a_n-1n)^T           (2-11)

上述式子中，r_ij为第i组数据的第j随机变量数值，p_cr,i为第i组数据的临界荷载值。

通过对{r_i}进行Box-Behnken矩阵法、CCD法或2ⁿ析因法等抽样方法，抽样次数为k次，通过上述有限元方程组(2-8)，获得(2-12)。

最不利临界荷载为

p*＝μ_pcr-3σ_pcr                   (2-12)

对大型结构，当随机变量很多的时候，这一过程的计算量很大，实际上并不需要在整个空间上拟合出响应面和精确的失效界面相吻合。只需在验算点附近一致。因这一区域对总的失效概率贡献最大，因此，展开点应选在验算点的位置。但在计算时并不知道验算点的位置，如果在计算中展开点选取范围很宽，验算点较容易落在该范围内，但是所得到的多项式对实际的失效函数的拟合度就交叉，反之，若取值范围较窄，验算点可能不落在该范围内，从而使所得是多项式不能与实际上失效函数在该点处相拟合。

实际计算中，若已知个基本变量的分布形式及分布参数，首先，以均值点位中心展开点的选择范围是

是随机变量r_i的均值，μ是均值，

是随机变量的方差，f是常数，可以根据工程上的3σ原则进行选择。

获得一次拟合的响应面法后，采用MCS计算结构临界荷载的概率分布或失效概率。在计算中为了更好的拟合曲线，可以使用对数形式的变量ln(p_cr),lnr,lna_i,lna_ij。

若已知网壳结构无缺陷时承载力设计值为p₀，该结构临界状态是由可靠度功能函数表达的，其形式为

Z＝p_cr(r)-5p₀                     (2-13)

随机缺陷下结构失效概率为

P_f＝P(p_cr(r)-5p₀≤0)                (2-14)

采用上述同样的方法即可求出。

以上方法提出了修正的基于响应面法的随机缺陷模态叠加法，这新的修正拟合算法，可快速准确地用响应面法来代替结构系统的性能，这一响应面由基本随机变量的统计特征计算得到,从而提高了算法的有效性和精确性。获得了响应面后，结合MCS算法。能得到结构可靠性分析结果。因该指出，响应面法对大型结构系统的可靠性问题很有效，相对计算量也不大。

另外，本发明中还提供了一种基于响应面法的网壳结构随机缺陷模态叠加方法，方法步骤为：(1)对网壳结构的结构杆件离散为单元并且对其节点进行编号，定义结构杆的材料性能参数、荷载和边界条件，完成数值模型的建立；(2)针对数值模型，参与组合的模态为特征值分析的前三阶，随机变量数选为3或4，采用CCD配点方法，进行随机缺陷模态叠加法分析，获取所求临界荷载方差、均值以及所消耗的CPU时间。其中步骤(2)中所述随机缺陷模态叠加法采用上述任一种随机缺陷模态叠加方法。

现针对图1中经典算例1：6角扁网壳结构，进行算例分析详细的说明。

杆件离散为Timeshenko梁单元，其截面积为317mm²。节点编号如图1所示，周边为3向固定铰支座。材料的弹性模量E＝3030MPa，剪切模量为1.096×10³MPa。荷载P向下作用顶点处。

实施例1

针对随机变量数为3的情况下，进行不同多项式拟合结果、不同配点拟合结果分析。

(1)不同多项式拟合结果对比

采用CCD配点方法，进行随机模态叠加法分析，获取所求临界荷载方差、均值等因素以及所消耗的CPU时间。同时，采用蒙塔卡罗方法对该数值模型进行计算，抽样次数是10000。两种分析方法的计算结果和相对误差如表1所示。

表1随机变量数为3时算例1的不同拟合多项式计算结果(CCD，n＝3)

表1中给出了响应面法(CCD)中，分别采用不同的拟合多项式，根据

和最低原则计算的临界不稳定荷载p^*，并出了蒙特卡罗方法下的p^*，以及二者相对误差。由表1可知，在同样的抽样次数下，CCD法中线性一阶多项式拟合出的临界荷载最低，但响应面的拟合验算误差相对较大，比MCS所求的值大1.6％。分析为软件在对CCD法中含有交叉项的响应面函数进行回归分析优化过程中，过滤了二次多项式中交叉项的影响，即真实采用的是完全二次多项式。这一结论在图2的柱状图中更为直观的获得。从所消耗的时间看，CCD法的CPU时间仅为2分19秒，远小于MCS法所需的20小时32分9秒(同一台计算平台)。

由上可知，CCD方法中使用完全的二次多项式响应面函数是可以满足拟合验算点，且于MCS的误差在0.5％之内，是可以采用的。而线性多项式下的计算结果虽比MCS法小3.1％，但响应面函数的拟合精度相对偏低，因而不建议采用。

总体说，MCS方法比响应面法精度高，但响应面法效率在较小的精度折损下，大大提高了计算效率。

(2)不同配点拟合结果对比

参与组合的模态为特征值分析的前三阶，采用二阶含交叉项的多项式拟合，针对不同配点，对该模型进行随机模态叠加法分同时，采用蒙塔卡罗方法对该数值模型进行计算，抽样次数是10000。两种分析方法的计算结果和相对误差如表2所示。

表2随机变量数为3时算例1的不同配点的二阶交叉多项式

表1中给出了响应面法(无交叉项多项式)中，分别采用不同的配点，根据

和最低原则计算的临界稳定荷载p^*，并出了蒙特卡罗方法下的p^*，以及二者相对误差；图3用于直观表达各种方法临界荷载的相对大小。由表2可知，在同样的抽样次数下，BBM所获得的响应面拟合误差较小，所获得临界荷载比CCD大，BBM所消耗的时间最少。不论何种配点的响应面法，在较小的精度折损下，大大提高了计算效率。综合来看，CCD与BBM各有优缺点，但BBM法将适合于随机模态叠加法。

实施例2

(1)采用随机变量为4，进行不同多项式拟合结果对比，获得数据见表3。

表3随机变量数为4时算例1的不同拟合多项式计算结果(CCD，n＝4)

采用ccd二阶多项式拟合所计算出的临界荷载值p^*比ccd线性多项式拟合结果更接近于MCS法下的临界荷载，它们与MCS的误差分别为：0.39％、0.39％、0.66％；

与随机变量为3时一样，随机变量为4时，在对ccd法中含有交叉项的响应面函数进行回归分析优化过程中，过滤了二次多项式中交叉项的影响；ccd二阶多项式拟合且考虑交叉项的程序所需cpu时间最少，仅为3m7s远小于mcs法的19h59m16s。

因此，在随机变量为4时，ccd法下采用二阶多项式拟合且考虑交叉项的程序效率较高且误差较小，宜优先采用。

将采用不同方法得到的最不利荷载画成柱状图，获得图4，由图4中通过比较发现ccd法的线性多项式拟合和二次多项式拟合出的最不利荷载p^*均大于蒙特卡洛法抽样结果。ccd法线性多项式拟合出的最不利荷载与mcs法较为接近，仅比后者大了0.86％；虽然ccd法中二次多项式拟合出的误差最小，但拟合出的最不利荷载值比mcs法大1.6％。

(2)采用随机变量为4，进行二阶交叉多项式拟合结果对比，获得数据见表4。

表4随机变量数为4时算例1的不同配点的二阶交叉多项式

表4中将ccd、bbm与MCS对比发现：bbm抽样法拟合的响应面所计算出相比于MCS抽样法误差为0.30％，小于ccd法误差值0.39％；同时bbm与ccd的有限元分析次数、CPU时间均比MCS的少。在采用回归分析方法拟合响应面函数时，发现bbm与ccd都过滤了交叉项的影响。

因此，说明在随机变量数为4时，bbm法拟合出响应面效果最好，但是误差与ccd法相差仅为0.09％，响应面函数宜使用完全的二次多项式。

见图5所示，将考虑二阶交叉项的ccd法和bbm法与MCS法抽样得到的最不利荷载画成柱状图，结果显示：ccd与bbm所获得的最不利荷载值均大于mcs法，其比mcs法分别高5.08％、3.3％。综合而言，考虑二阶交叉项的bbm法抽样获得最不利荷载值。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (7)

1.一种基于响应面法的随机缺陷模态叠加方法，其特征在于：步骤为：

(1)建立结构缺陷模态组合系数的概率模型，测定验算点相关数据；利用多项式响应面法，将临界荷载拟合表达式

表达为多项式的形式，根据临界荷载拟合表达式

计算确定表达式待定系数，并由计算得到的待定系数确定临界荷载拟合表达式

(2)确定临界荷载：临界荷载采用以下方法确定：

其中，

为临界荷载p_cr(r)的均值，

为临界载荷p_cr(r)的方差，p^*为最不利临界载荷；

(3)计算结构临界荷载的失效概率：随机缺陷下结构失效概率采用以下方法确定：

P_f＝P(p_cr(r)-5p₀≤0)                (2-14)

其中，P_f为随机缺陷下结构失效概率，p₀为网壳结构无缺陷时承载力设计值，p_cr为响应值-临界荷载；

其中，步骤(1)中，结构缺陷模态组合系数的概率模型的建立方法，包括以下步骤：

(一)获取线性屈曲的特征值和相应的屈曲模态：

根据结构刚度矩阵[K_T]，求解线性屈曲的特征值λ_i(i＝1,2,…,m)和相应的屈曲模态{U_i}，其中m是参与组合的模态阶数；

(二)获取随机有限元刚度方程：

(2.1)对结构节点进行编号，依次为1，2，L，n，L、n均为正整数；

(2.2)获得假定结构的任一缺陷模式向量{△X}'：其中

设定m为模态参与阶数，r₁，r₂，r_i……r_m为参与系数，是独立随机变量，{U_i}为结构第i阶线性屈曲模态；

(2.3)获得幅值为R的缺陷模式向量{△X}：对{△X}'进行幅值调整，以获得幅值为R的缺陷模式向量{△X}；

(2.4)获得节点坐标{X}的切线刚度矩阵[K_T]：先将(2.3)中获得{△X}计算到{X}＝{△X}+{X₀}中，设定：{X}引入缺陷后结构的节点坐标，{X₀}为3n维完善结构节点坐标向量，n为结构节点数，{△X}为3n维结构节点位置缺陷向量；再采用引入缺陷后结构的节点坐标{X}形成；

(2.5)获得随机有限元刚度方程：采用(2.4)[K_T]获得[K_T]{△a}＝{Q}-{F}，式中[K_T]为形成当前有缺陷结构节点坐标{X}的切线刚度矩阵，{△a}、{F}均为随机变量r_i的函数，{Q}是荷载列向量；

(三)获得随机变量输入后的响应值-临界荷载p_cr：

利用弧长法，求解步骤(二)中随机有限元刚度方程，获得结构在随机变量输入后的响应值：临界荷载p_cr，及临界荷载的拟合值

(3.1)若取二次完全多项式，则

(3.2)若不考虑二次交叉项系数，则

(3.3)简单取

一次多项式，则

设定(3.1)-(3.3)中，a₀，a_i，a_ij，a_ii是待定系数；r_i，r_ir_j为参与系数；

为样本结构下临界荷载p_cr的近似值，i，j均为正整数；r为样本结构的参与系数。

2.根据权利要求1所述的一种基于响应面法的随机缺陷模态叠加方法，其特征在于：步骤(1)中，对大型结构来说，验算点的选择方式为：已知各基本变量的分布形式及分布参数，以均值点为中心展开,验算点的选择范围是

其中，

是随机变量ri的均值，μ是均值，

是随机变量的方差，f是常数，根据工程上的3σ原则选择。

3.根据权利要求1所述的一种基于响应面法的随机缺陷模态叠加方法，其特征在于：步骤(1)中，临界荷载拟合表达式

为二次完全多项式，则

或者临界荷载拟合表达式

为一次多项式，则

式中，a₀，a_i，a_ii，a_ij是待定系数。

4.根据权利要求1所述的一种基于响应面法的随机缺陷模态叠加方法，其特征在于：步骤(1)中，待定系数的计算方法为：用最小二乘估算法求a_i,a_ij，假定得到k组数据为

(r)_i，其中i＝1,2,…,k，必须大于等于相应的未知参数的个数且为正整数，由如下(2-7)、(2-8)两式可以求出待定系数：

([r]^T[r]){A}＝[r]^T{P_cr}                    (2-7)

{A}＝([r]^T[r])^-1[r]^T{P_cr}                   (2-8)

式中，

5. {P_cr}＝(p_cr,1 p_cr,2 L p_cr,k)^T       (2-10)

{A}＝(a₀ a₁ L a_n a₁₁ L a_nn a₁₂ L a_n-1n)^T    (2-11)

上述式子中，r_ij为第i组数据的第j随机变量数值，p_cr,i为第i组数据的临界荷载值。

5.根据权利要求1所述的一种基于响应面法的随机缺陷模态叠加方法，其特征在于：步骤(2)中，临界载荷的获得方法为：通过{r_i}进行Box-Behnken矩阵法、CCD法或2ⁿ析因法中任一种抽样方法，抽样次数为k次，结合有限元方程组，获得

为临界荷载p_cr(r)的均值，

为临界载荷p_cr(r)的方差，p^*为最不利临界载荷。

6.一种基于响应面法的网壳结构随机缺陷模态叠加方法，其特征在于：方法步骤为：(1)对网壳结构的结构杆件离散为单元并且对其节点进行编号，定义结构杆的材料性能参数、荷载和边界条件，完成数值模型的建立；(2)针对数值模型，参与组合的模态为特征值分析的前n阶，即随机变量数为n，n为正整数，采用CCD配点方法，进行随机缺陷模态叠加法分析，获取所求临界荷载方差、均值以及所消耗的CPU时间。

7.根据权利要求6所述的一种基于响应面法的网壳结构随机缺陷模态叠加方法，其特征在于：步骤(2)中所述随机缺陷模态叠加法采用权利要求1至5中任一种随机缺陷模态叠加方法。

Images