完全基于降阶模型的结构非线性屈曲平衡路径跟踪方法

摘要

一种完全基于降阶模型的结构非线性屈曲平衡路径跟踪方法，根据用户要求的结构承载响应曲线长度进行结构非线性屈曲分析，所得到的加载点沿载荷方向位移随载荷曲线变化。本发明充分考虑结构位移和降阶模型展开式中的残余项，并基于结构残余力将上述残余项迭代至零，从而实现对非线性预测解的精确修正，在解的修正阶段无需再返回有限元全模型进行计算，所提出的降阶技术应用在路径跟踪的预测环节，并在解的修正阶段得以实现，进而构造了一种完全“存粹”的路径跟踪降阶方法，具有结构非线性屈曲平衡路径分析效率高、精度好的特点，能够准确经济地评估结构的真实承载能力，并且计算时间仅为传统方法的1/5，大幅提高了结构非线性屈曲分析的计算效率。

说明书

技术领域

本发明属于结构力学分析技术领域，涉及一种完全基于降阶模型来跟踪结构非线性屈曲平衡路径的分析方法，特别适合于对存在屈曲失稳的结构的承载能力的精确高效评估。

背景技术

薄壁结构作为航空航天工程的典型结构构型具有重量轻、承载效率高及可设计性强等技术优势。此类结构的极限承载能力往往由其屈曲性能所决定，呈现出明显的几何非线性效应及工程应用的风险性，因此结构的非线性屈曲响应分析是轻质薄壁结构设计的重要技术环节。

为计及几何非线性效应，当前计算结构屈曲响应的常规方法是基于Newton增量迭代技术的非线性有限单元法(International Journal of Solids and Structures,2002,139(3):689-698)。该方法需要对有限元离散后的大规模非线性方程组进行迭代求解，存在求解效率低、计算规模大等问题。极大影响了飞行器全尺寸结构研发的周期、成本及其工程预测的可靠性，增大了研制过程的技术与经济风险。

Koiter-Newton降阶法是近些年所提出的、将Koiter摄动理论与Newton法的增量迭代思想相结合的、结构非线性路径跟踪计算方法(International Journal forNumerical Methods in Engineering,2013,96(12):763–786)。该方法基于模型降阶思想在计算效率上的优势以及在精算精度上的可靠性已得到充分验证。方法在每一个载荷步中采用非线性预测+线性修正的计算跟踪策略，其中非线性预测是通过对采用Koiter摄动理论所获得的降阶模型进行求解而得到的，而线性修正则还需要采用传统Newton法基于有限元全模型来开展。由此可以看出，此方法在路径跟踪的整个过程中并不是一种完全“纯粹”的降阶方法，其在计算效率上的优势还没有完全被挖掘。发明专利(申请号201510737131.2)所提出的“一种高效高精度地探测及跟踪结构屈曲分支路径的方法”就是完全基于Koiter-Newton降阶法来实现的，故也存在上述问题。

因此，本发明提出一种能够完全基于降阶模型来跟踪结构非线性屈曲平衡路径的新方法。

发明内容

为克服现有技术中承载的当前降阶方法在路径跟踪的修正阶段必须回到有限元全模型来计算的瓶颈，本发明提出了一种完全基于降阶模型的结构非线性屈曲平衡路径跟踪方法。

本发明的具体过程是：

步骤1，在已知平衡状态点处建立薄壁曲板结构的降阶模型；

所述已知平衡状态点处建立薄壁曲板结构在已知平衡状态点处的降阶模型为：

所建立的降阶模型是一个(m+1)阶的非线性方程组，其阶数一般小于10。上式中φ是载荷系数向量，它是第一个分量为μ其余分量均为0的单位向量，其中μ是状态增量所对应的外载荷向量的载荷系数。

所述建立薄壁曲板结构在已知平衡状态点处的降阶模型的具体过程是：

Ⅰ定义该薄壁曲板结构平衡路径上的已知平衡状态点为(u₀，λ₀)，其中u₀是结构在该已知平衡状态点处的结构位移，λ₀是结构在该已知平衡状态点处的外载荷向量的载荷系数。在第一个计算步中，选取从结构的未变形状态点处开始，此时u₀和λ₀分别对应的是结构在未变形状态处的位移和载荷系数；u₀和λ₀均为零。

在已知平衡点附近的定义一个未知平衡状态点为(q,λ)，其中q是结构在该未知平衡状态点处的结构位移，λ是相应的外载荷向量的载荷系数。所述的已知平衡点附近为载荷系数λ₀变化量小于10％的位置。

定义所述的已知平衡状态点和未知平衡状态点之间的状态增量为(u,μ)，其中u是结构状态增量所对应的结构位移，μ是相应的外载荷向量的载荷系数。所述的已知平衡状态点(u₀,λ₀)和未知平衡状态点(q,λ)之间分别满足以下的两个关系式：

q＝u₀+u(1)

λ＝λ₀+μ(2)

Ⅱ在已知平衡状态点(u₀,λ₀)处建立离散化的薄壁曲板结构非线性有限元平衡方程：

f(q)＝λf_ext(3)

其中：向量q是结构位移；f(q)是结构的非线性内力；向量f_ext是结构外载荷；λ是外载荷向量的载荷系数。

Ⅲ在已知平衡点(u₀,λ₀)处对薄壁曲板结构进行线性特征值屈曲分析：

K_tw＝zK_gw(4)

其中，矩阵K_t和K_g分别为薄壁曲板结构的切线刚度和几何刚度；w是薄壁曲板结构的屈曲模态；z是薄壁曲板结构的屈曲载荷。本实施例中，提取结构的m个密集屈曲模态构建降阶模型，并将K_g与屈曲模态w的乘积定义为摄动载荷向量

Ⅳ将所述离散化的结构非线性有限元平衡方程(3)在已知平衡状态点(u₀,λ₀)处关于位移u展开至三阶：

L(u)+Q(u,u)+C(u,u,u)＝Fφ(5)

其中：L、Q、C分别是展开式(5)中的一次算子、二次算子和三次算子；F是载荷向量矩阵，它的第1列为外载荷向量f_ext，其余各列由摄动载荷向量构成；φ是载荷系数向量，它的第一个分量为μ，其余各分量均为0的单位向量，其中μ是状态增量所对应的外载荷向量的载荷系数。

Ⅴ将位移u和载荷系数向量φ关于摄动参数向量ξ分别展开至三阶项：

u＝u_αξ_α+u_αβξ_αξ_β(7)

其中：变量下标α、β＝1,2,…,m+1，m是指该结构密集屈曲模态的个数；u_α是结构的一阶位移场，u_αβ是结构的二阶位移场；ξ_α、ξ_α是摄动参数向量ξ的各分量；分别为展开式(6)中待求解的一次算子、二次算子和三次算子。

Ⅵ将公式(5)中的位移u和载荷系数向量φ分别用式(6)和式(7)进行替换；令摄动参数向量ξ的各次幂的系数为零，即可获得式(6)和(7)中各算子的具体表达。

各算子分别为摄动参数向量ξ的一次、二次和三次函数表达式。步骤2，计算薄壁曲板结构非线性响应的预测解

采用Newton增量迭代技术求解该降阶模型，获得该方程的一个非线性解(ξ,μ)。将其分别代入到式(7)、(1)和(2)中，即得到此时薄壁曲板结构非线性响应的预测解(q,λ)。

步骤3，对薄壁曲板结构非线性响应的预测解进行修正

Ⅰ计算该非线性预测解所对应的结构残余力r^k：

r^k＝λf_ext-f(q)(9)

其中，上标k表示当前修正状态，f(q)是结构的非线性内力，向量f_ext是结构外载荷，λ是外载荷向量的载荷系数。

Ⅱ若在降阶模型(8)和位移展开式(7)中考虑残余项，则降阶模型和位移展开式分别重新写为：

其中，是降阶模型中的残余项；是位移展开式中的残余项。

Ⅲ计算降阶模型中残余项的增量和位移展开式中残余项的增量

其中，K_t是结构在已知平衡状态点处的切线刚度矩阵,F是摄动载荷向量，-F^T是指将摄动载荷向量取负值，r^k是结构残余力。求解线性方程组(12)即可分别得到降阶模型中残余项的增量和位移展开式中残余项的增量

Ⅲ采用计算得到的残余项增量更新降阶模型(10)中的残余项并采用计算得到的残余项增量更新位移展开式(11)中的残余项即：

上式中和是更新前的残余项，和是更新后的残余项，和是残余项增量。公式(12)和(13)对结构的残余项和进行了更新。

Ⅴ基于更新后的残余项通过公式(15)重新计算降阶模型(10)，

获得针对原先解(ξ,μ)的一个更新解(ξ^k+1,μ^k+1)，其中ξ^k+1是更新后的摄动参数解，μ^k+1是更新后的载荷系数解：

上式中是更新后的残余项(13)，φ是载荷系数向量，分别为展开式(6)中的一次算子、二次算子和三次算子。

Ⅵ基于更新后的摄动参数解ξ^k+1和更新后的残余项通过公式(16)重新计算位移u：

上式中u_α是结构的一阶位移场，u_αβ是结构的二阶位移场，是摄动参数更新解ξ^k+1的各分量，是更新后的残余项。

Ⅶ将更新后的载荷系数解μ^k+1和更新后的位移u代入式(1)和(2)即得到更新修正后的结构非线性预测解(q,λ)。到此，第一个修正步结束。

重复本步骤的修正过程，开始第二个修正步，直至式(12)中残余项增量和的二范数均小于使用者根据精度要求设定的某一给定阀值(0.0001)，则表示更新后的预测解(q,λ)满足精度要求，停止修正。返回步骤2，继续求解降阶模型获得下一个预测解，再按照步骤3对该解进行修正，得到第二个数据点。

重复所述修正步的过程，依次获得第1个计算步中的所有数据点。

至此，完成了第一计算步的循环过程。

步骤4，第二计算步的修正循环

Ⅰ 若步骤3中的的修正步次数大于10，则表明该降阶模型自身精度不足，则停止修正，开始第二个计算步。

Ⅱ 返回步骤1，以第一个计算步所获得的沿曲线的最后一个三角标记处重新建立降阶模型。重复步骤1～3即获得第二个计算步中的所有数据点。至此，完成了第二计算步的循环过程，进入下一个计算步；

步骤5，第三计算步的修正循环

Ⅰ 若步骤4中重复步骤3的修正步次数大于10，则表明该降阶模型自身精度不足，则停止修正，开始第三个计算步。

Ⅱ 返回步骤1，以第二个计算步所获得的沿曲线的最后一个空心圆标记处重新建立降阶模型。重复步骤1～3即获得第三个计算步中的所有数据点。至此，完成了第三计算步的循环过程，进入下一个计算步，直至获得满足用户要求的结构承载响应曲线长度。

本发明与现有技术相比的优点在于：

当前的Koiter-Newton降阶方法虽然在每一个载荷步的预测阶段只需要对降阶模型进行求解，从一定程度上显著减少了结构非线性平衡路径跟踪的计算效率，但该方法在对预测解进行修正时仍然需要返回有限元全模型来计算，从而并没有实现平衡路径跟踪中的全程降阶求解。为此，本发明充分考虑结构位移和降阶模型展开式中的残余项，并基于结构残余力将上述残余项迭代至零，从而实现对非线性预测解的精确修正。该技术在解的修正阶段不需要再返回有限元全模型进行计算，所提出的降阶技术不但应用在了路径跟踪的预测环节，也在解的修正阶段得以实现，进而构造了一种完全“存粹”的路径跟踪降阶方法。

本发明具有结构非线性屈曲平衡路径分析效率高、精度好的特点，可准确经济地评估结构的真实承载能力。

附图说明

图1为本发明的示意图

图2为受横向集中载荷作用的曲板示意图；

图3为结构承载响应曲线示意图；

图4为本发明的流程图。图中：

1.第一计算步；2.第二计算步；3.第三计算步；4.实线。

具体实施方式

本实施例是对图2所示薄壁曲板通过完全基于降阶模型以跟踪结构非线性屈曲平衡路径的分析方法。所述的薄壁曲板的直边长度为508mm，曲边所对应的曲率半径为2540mm，弧度为0.1rad，板厚为12.7mm，两直边简支约束，两曲边自由。板的弹性模量3102.75MPa，泊松比为0.3。该薄壁曲板上表面几何中心处受到一个横向集中点载荷。

本实施例根据用户要求的结构承载响应曲线长度，对该曲板进行结构非线性屈曲分析，所得到的加载点沿载荷方向位移随载荷变化曲线如图3所示。

本实施例具体步骤如下：

步骤1，在已知平衡状态点处建立薄壁曲板结构的降阶模型；

Ⅰ定义该薄壁曲板结构平衡路径上的已知平衡状态点为(u₀，λ₀)，其中u₀是结构在该已知平衡状态点处的结构位移，λ₀是结构在该已知平衡状态点处的外载荷向量的载荷系数。在第一个计算步中，选取从结构的未变形状态点处开始，如图3中的坐标原点，因而此时u₀和λ₀分别对应的是结构在未变形状态处的位移和载荷系数；u₀和λ₀均为零。

在已知平衡点附近的定义一个未知平衡状态点为(q,λ)，其中q是结构在该未知平衡状态点处的结构位移，λ是相应的外载荷向量的载荷系数。所述的已知平衡点附近为载荷系数λ₀变化量小于10％的位置。

定义所述的已知平衡状态点和未知平衡状态点之间的状态增量为(u,μ)，其中u是结构状态增量所对应的结构位移，μ是相应的外载荷向量的载荷系数。所述的已知平衡状态点(u₀,λ₀)和未知平衡状态点(q,λ)之间分别满足以下的两个关系式：

q＝u₀+u(1)

λ＝λ₀+μ(2)

Ⅱ在已知平衡状态点(u₀,λ₀)处建立离散化的薄壁曲板结构非线性有限元平衡方程：

f(q)＝λf_ext(3)

其中：向量q是结构位移；f(q)是结构的非线性内力；向量f_ext是结构外载荷；λ是外载荷向量的载荷系数。

Ⅲ在已知平衡点(u₀,λ₀)处对薄壁曲板结构进行线性特征值屈曲分析：

K_tw＝zK_gw(4)

其中，矩阵K_t和K_g分别为薄壁曲板结构的切线刚度和几何刚度；w是薄壁曲板结构的屈曲模态；z是薄壁曲板结构的屈曲载荷。本实施例中，提取结构的m个密集屈曲模态构建降阶模型，并将K_g与屈曲模态w的乘积定义为摄动载荷向量

Ⅳ将所述离散化的结构非线性有限元平衡方程(3)在已知平衡状态点(u₀,λ₀)处关于位移u展开至三阶：

L(u)+Q(u,u)+C(u,u,u)＝Fφ(5)

其中：L、Q、C分别是展开式(5)中的一次算子、二次算子和三次算子；F是载荷向量矩阵，它的第1列为外载荷向量f_ext，其余各列由摄动载荷向量构成；φ是载荷系数向量，它的第一个分量为μ，其余各分量均为0的单位向量，其中μ是状态增量所对应的外载荷向量的载荷系数。

Ⅴ将位移u和载荷系数向量φ关于摄动参数向量ξ分别展开至三阶项：

u＝u_αξ_α+u_αβξ_αξ_β(7)

其中：变量下标α、β＝1,2,…,m+1，m是指该结构密集屈曲模态的个数；u_α是结构的一阶位移场，u_αβ是结构的二阶位移场；ξ_α、ξ_α是摄动参数向量ξ的各分量；分别为展开式(6)中待求解的一次算子、二次算子和三次算子。

Ⅵ将公式(5)中的位移u和载荷系数向量φ分别用式(6)和式(7)进行替换；令摄动参数向量ξ的各次幂的系数为零，即可获得式(6)和(7)中各算子的具体表达。继而建立薄壁曲板结构在已知平衡状态点处的降阶模型：

所建立的降阶模型是一个(m+1)阶的非线性方程组，其阶数一般小于10。上式中φ是载荷系数向量，它是第一个分量为μ其余分量均为0的单位向量，其中μ是状态增量所对应的外载荷向量的载荷系数

步骤2，计算薄壁曲板结构非线性响应的预测解

薄壁曲板结构在已知平衡状态点处的降阶模型(8)是一个非线性方程组，因此采用经典的Newton增量迭代技术求解该降阶模型，获得该方程的一个非线性解(ξ,μ)。将其分别代入到式(7)、(1)和(2)中，即可得到此时薄壁曲板结构非线性响应的预测解(q,λ)。

步骤3，对薄壁曲板结构非线性响应的预测解进行修正

Ⅰ计算该非线性预测解所对应的结构残余力r^k：

r^k＝λf_ext-f(q)(9)

其中，上标k表示当前修正状态，f(q)是结构的非线性内力，向量f_ext是结构外载荷，λ是外载荷向量的载荷系数。

Ⅱ若在降阶模型(8)和位移展开式(7)中考虑残余项，则降阶模型和位移展开式分别重新写为：

其中，是降阶模型中的残余项；是位移展开式中的残余项。

Ⅲ计算降阶模型中残余项的增量和位移展开式中残余项的增量

其中，K_t是结构在已知平衡状态点处的切线刚度矩阵,F是摄动载荷向量，-F^T是指将摄动载荷向量取负值，r^k是结构残余力。求解线性方程组(12)即可分别得到降阶模型中残余项的增量和位移展开式中残余项的增量

Ⅲ采用计算得到的残余项增量更新降阶模型(10)中的残余项并采用计算得到的残余项增量更新位移展开式(11)中的残余项即：

上式中和是更新前的残余项，和是更新后的残余项，和是残余项增量。公式(12)和(13)对结构的残余项和进行了更新。

Ⅴ 基于更新后的残余项通过公式(15)重新计算降阶模型(10)，

获得针对原先解(ξ,μ)的一个更新解(ξ^k+1,μ^k+1)，其中ξ^k+1是更新后的摄动参数解，μ^k+1是更新后的载荷系数解：

上式中是更新后的残余项(13)，φ是载荷系数向量，分别为展开式(6)中的一次算子、二次算子和三次算子。

Ⅵ 基于更新后的摄动参数解ξ^k+1和更新后的残余项通过公式(16)重新计算位移u：

上式中u_α是结构的一阶位移场，u_αβ是结构的二阶位移场，是摄动参数更新解ξ^k+1的各分量，是更新后的残余项。

Ⅶ 将更新后的载荷系数解μ^k+1和更新后的位移u代入式(1)和(2)即得到更新修正后的结构非线性预测解(q,λ)。到此，一个修正步结束。

重复本步骤的修正过程，开始下一个修正步，直至式(12)中残余项增量和的二范数均小于使用者根据精度要求设定的某一给定阀值(0.0001)，则表示更新后的预测解(q,λ)满足精度要求，停止修正，该预测解对应图3中的第一个数据点，即从坐标原点起沿曲线的第一个三角标记，然后返回步骤2，继续求解降阶模型获得下一个预测解，再按照步骤3对该解进行修正，得到图3中的第二个数据点，即从坐标原点起沿曲线的第二个三角标记，这样就依次获得图三中第1个计算步中的所有数据点，即沿曲线的所有三角标记；

至此，完成了第一计算步1的循环过程。

步骤4，第二计算步的修正循环

Ⅰ 若所述第一计算步的修正步次数大于10，则表明该降阶模型自身精度不足，则停止修正，开始第二个计算步2。

Ⅱ 返回步骤1，以第一个计算步所获得的沿曲线的最后一个三角标记处重新建立降阶模型。重复步骤1～3即获得第二个计算步中的所有数据点。至此，完成了第二计算步2的循环过程。

所述的第二个计算步是图3中沿曲线的所有空心圆标记。

步骤5，第三计算步的修正循环

Ⅰ 若所述第二计算步的修正步次数大于10，则表明该降阶模型自身精度不足，则停止修正，开始第三个计算步3。

Ⅱ 返回步骤1，以第二个计算步所获得的沿曲线的最后一个空心圆标记处重新建立降阶模型。重复步骤1～3即获得第三个计算步中的所有数据点。至此，完成了第三计算步的循环过程。

经过第三计算步3的修正循环，获得满足用户要求的结构承载响应曲线长度。

所述的第三个计算步是图3中沿曲线的所有五角星标记。

通过上述重复循环过程，采用若干个计算步即可获得完整的结构非线性平衡路径。本实施例采用了三个计算步，即得到的平衡路径曲线，如图3所示。

本实施例中，采用传统算法，如基于全模型的Newton弧长法，获得实曲线4。曲线上的数据点均由本发明所提出的方法获得，三种不同形状的数据点表示本实施例采用了三个不同的计算步进行跟踪，每个计算步获得结构平衡路径长度中的一部分。

由图3可知本实施例对结构非线性平衡路径与采用当前传统算法得到的曲线完全吻合，其跟踪计算精度达到要求。为获得图3所示曲线，传统方法需要的CPU计算时间为92s，而本发明方法的计算时间仅为18s，计算时间仅为传统方法的1/5，大幅提高了结构非线性屈曲分析的计算效率。