CN112966421B - 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 - Google Patents

Abstract

本发明涉及一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法，属于结构稳定性校核技术领域。本发明将p型有限元法应用于结构稳定性分析领域，求解在相应几何和约束条件下的应力场，得到相应的几何刚度矩阵并应用Lanczos迭代方法计算屈曲特征值，进而得到屈曲载荷因子和屈曲形状。该方法能减少计算成本，提高计算的收敛率和计算精度。

Description

使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形 状的计算方法

技术领域

本发明涉及一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法，属于结构稳定性校核技术领域。

背景技术

屈曲是结构常见的破坏形式之一，尤其是对于一些长、薄结构，具备这样性质的板壳结构广泛应用于航空、航天、航海、军工等领域中。飞机机翼、尾翼上的翼面壁板、梁腹板和机身上的蒙皮、隔框等常采用薄壁结构，这些结构都具有一个共同的特征，即一个方向的尺寸大大小于另外两个方向的尺寸，可以简化为板壳结构，当这种结构受到压缩、弯曲和剪切等外荷载单独或联合作用时，屈曲失稳是最常见的失效模式之一。薄板的屈曲分析在工程设计中至关重要，对于大多数薄板结构，屈曲失稳先于强度破坏。

对于屈曲分析，主要应用伽辽金法、摄动法和有限元法。伽辽金法在求解过程中不涉及变分问题，简单易行，应用较广，但是在微分方程算子不是线性算子时很难求解；摄动法在求解过程中需要解渐近展开式，由于渐近级数一般是发散的，在求解时可能在级数展开项的前几项表现为逐渐收敛形式，而在级数项增加到一定程度后又表现为发散形式，要解决这一问题就需要较高的数学专业技能，增加了求解难度；有限元法采用将复杂结构离散化为一定数量的单元求数值解，为了提高数值解的计算精度，必须划分更小的单元。传统有限元法划分单元较多，计算效率偏低。

本发明专利运用p型有限元法进行薄板结构的特征值屈曲分析，利用p型有限元法单元数量少、前处理少及收敛速度快的优点，提出了一种基于p型有限元法计算薄板结构屈曲载荷系数和相应屈曲形状的方法。

发明内容

本发明将p型有限元法应用于结构稳定性分析领域，求解在相应几何和约束条件下的应力场，得到相应的几何刚度矩阵并应用Lanczos迭代方法计算屈曲特征值，进而得到屈曲载荷因子和屈曲形状。该方法能减少计算成本，提高计算的收敛率和计算精度。本发明通过以下技术方案实现。

使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法，其步骤如下：

步骤1、确定薄板结构的几何尺寸和材料参数，建立三维有限元模型，该步骤包括：

步骤1.1、根据结构的实际受力情况，确定结构的受力特征；

步骤1.2、根据结构所使用的材料确定材料属性，如有必要应通过实验测定材料的弹性模量、泊松比。

步骤1.3、建立三维有限元模型。

步骤1.4、施加载荷和边界条件：根据确定的载荷和边界条件，在相应的边界上施加相对应的载荷和边界条件；

步骤2使用p型有限元法计算在对应荷载和约束条件下的线性解。具体步骤如下：

步骤2.1、结合所指定的插值多项式阶数，依据方程

Ka＝F (1)

求解结构位移列阵a，其中，K为整体刚度矩阵，K＝∑G^TK_eG,K_e为单元刚度矩阵，a为结构位移列阵，F为结构节点载荷列阵，K_e，G，表示结构节点自由度与单元节点自由度的转换矩阵。

步骤2.2、依据单元刚度矩阵公式

单元等效节点载荷列阵

其中：单元内部节点力

外部节点力

上式(2)至上式(5)中，Ω表示为单元内部；B＝LN_I，B表示应变矩阵，L表示微分算子，N_I表示高阶形函数矩阵，I表示单元的标号，Γ_t表示单元的外部边界；D表示应力矩阵，b表示体力，表示荷载边界条件，这些量由步骤1中确定的材料属性、载荷及边界条件计算得到；把上式(2)至上式(5)代入到上式(1)中，再对线性方程组(1)进行求解，则得到了结构位移列阵a；

步骤2.3、依据公式a_I＝Ga进一步得到了单元位移列阵a_I；

步骤2.4、依据公式u＝N_Ia_I，ε＝Lu，σ＝Dε＝DBa_I，求解得到应力场σ、应变场ε、位移场u。

本发明采用p型有限元法求解位移场和应力应变场。p型有限元法可在固定网格下，通过采用阶谱类型的高阶形状函数来提高计算精度，其主要特征在于在三维标准六面体单元上采用了一种阶谱类型的高阶形状函数。本发明采用六面体网格单元对薄板屈曲载荷因子和屈曲形状进行数值模拟分析。

阶谱类型的三维标准六面体单元高阶形状函数的构造形式如下：

三维标准六面体单元如图1所示，其形函数为：

(1)节点模式形状函数：

节点模式形状函数共有8个，其定义如下：

(2)边模式形状函数：

边模式形状函数共有12(p-1)个，其定义如下：

与连接节点1和节点2的边相关联的边模式形状函数为：

其中，

类似地，与连接节点2和节点3的边相关联的边模式形状函数为：

其余边模式形状函数以此类推可得。

(3)面模式形状函数：

面模式形状函数共有3(p-2)(p-3),p≥4个，其定义如下：

与连接节点1,2,5,6的面相关联的面模式形状函数为：

其中，下标m＝m(i,j)由编号约定决定。

其余面模式形状函数以此类推可得。

(4)内部模式形状函数：

内部模式形状函数共有p≥6个，其定义如下：

其中，下标m＝m(i,j,k),i,j,k＝2,3,…,p-4,i+j+k＝6,7,…,p。

插值多项式的阶次可依次从p＝1逐渐提高，提升插值多项式低阶的刚度矩阵可继续沿用，只需计算高阶部分，避免了刚度矩阵低阶部分的重复计算，具有良好的承袭性，节约了前处理的成本。

步骤3：通过所求得的线性解中的应力场，计算所对应的几何刚度矩阵；再使用几何刚度矩阵进行屈曲载荷因子特征值求解：具体为

步骤3.1、通过计算得到的应力场求得几何刚度矩阵；

步骤3.2、将所求得的几何刚度矩阵用于特征值计算，将广义的特征值问题Kφ＝λMφ转化为标准特征值问题Kφ＝λφ，即将M化为单位阵；

其中K为刚度矩阵，M质量矩阵，φ为特征矢量，λ即为所求得的屈曲载荷因子特征值；

采用Lanczos迭代方法进行特征值计算；标准特征值问题：Kφ＝λφ

任何初始向量U,||U||_k＝1，U₀＝0，其中||||_k为k范数

采用三项递推公式：{U_k+1}＝(K{U_k}-α_k{U_k}-β_k{U_k-1})/β_k+1

其中，β₁＝0

α_k＝{U_k}^TK{U_k}

β_k+1＝||K-{U_k}-α_k{U_k}-β_k{U_k-1}||₂

这里k＝1,2,3,...,m-1≤n，||||₂代表2范数，m为最高阶特征值个数，n为方阵阶数，完成迭代计算后，最终得到两组解：用来建立预屈曲应力状态的线性解和特征值屈曲载荷因子解；

步骤4：根据所求得的屈曲载荷因子是否满足精度要求，并判断能量范数误差是否满足一定的要求，如屈曲载荷因子未满足精度要求或者能量范数误差过大，提高插值多项式的阶数并增加网格剖分数量，返回步骤2；屈曲载荷因子满足精度要求，根据屈曲载荷因子获得相应屈曲形状。

所述步骤2中p型有限元法在单元数量确定的情况下，通过提升插值多项式的阶数来提高计算结果的精度。当单元厚度较小时，使用低阶插值函数会使剪切刚度占主导地位，发生剪切锁定，因此插值多项式的阶数应大于等于3。

上述p型有限元法的插值多项式是阶谱类型形状函数，提升插值多项式阶数时，低阶时的刚度矩阵可继续沿用，只需计算新增高阶项，装配形成新的总体刚度矩阵，具有较好的承袭性

上述相应屈曲形状是通过以下方式得到：屈曲载荷因子乘以步骤1中所施加荷载得到屈曲临界荷载，再以屈曲临界荷载作为外荷载以步骤2.1—2.4所述方式计算得到位移场，据此绘出薄板位移图即为相应的屈曲形状。

本发明的有益效果是：该方法能减少计算成本，提高计算的收敛率和计算精度。

附图说明

图1是三维有限元网格单元示意图；

图2是薄钢板结构示意图；

图3是薄钢板屈曲分析网格示意图；

图4是实施例1和实施例2屈曲形状示意图。

具体实施方式

下面结合附图和具体实施方式，对本发明作进一步说明。

实施例1

根据图2所示薄钢板构件，其中薄板长为a＝3000mm，宽b＝a/m，厚h＝10mm，薄板两端受压应力σ＝1.0MPa作用，边界条件为四边简支约束。薄钢板的材料参数：杨氏模量E＝205800MPa，泊松比ν＝0.3，m为变化参数，本实施例实例选择从1到2。

模型的网格划分如图3所示。当薄板的几何参数变化时，网格采用均匀剖分，并且单元数量均为9。

由采用平衡法求解的理论公式可得此时对应的屈曲载荷因子的理论解。

上式中：k为屈曲系数，与薄板的长宽比m＝a/b有关，当m为整数时，k＝k_min＝4，E为弹性模量，ν为泊松比，t为薄板厚度，b为薄板的宽度，下表为m＝1,…,5时临界屈曲载荷的理论值。

表1

基于p型有限元求解该薄板构件的线性解，并通过线性解得到几何刚度矩阵，再得到临界屈曲载荷因子，将屈曲载荷因子与所施加的外力相乘最终得到屈曲荷载。下表给出了p＝3,4,…,8且m＝1,2时有限元解的能量范数误差和对应的约束条件下得到的临界屈曲荷载，及与理论解相比的相对误差。

表2

基于p型有限元法计算临界屈曲载荷，随着插值多项式阶数p的增加，有限元解及有限元解能量范数误差逐渐收敛。基于上表可看出，本发明计算临界屈曲载荷时，具有前处理少，收敛速度快，采用较低的计算成本，便可获得较高精度的数值解。

实施例2

实例2模型与实例1类似，在相同网格和材料属性条件下，改变薄板的长宽比，继续计算m＝4,5时的临界屈曲载荷，具体数据如下表：

表3

由上表可看出改变薄板的长宽比，网格条件相同条件下，对比理论解，此时临界屈曲载荷的精度较高，且具有良好的收敛性。

实施例1与实施例2屈曲形状如图4所示。

以上对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (2)

1.一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法，其特征在于步骤如下：

步骤1：构建薄板构件三维有限元模型，该步首先创建薄板构件的几何模型，接着对所建立几何模型进行网格剖分；其次对要进行屈曲分析的薄板构件进行材料参数的设定；然后施加对应的荷载与边界条件，其步骤如下：

步骤1、确定薄板结构的几何尺寸和材料参数，建立三维有限元模型，该步骤包括：

步骤1.1、根据结构的实际受力情况，确定结构的受力特征；

步骤1.2、根据结构所使用的材料确定材料属性，通过实验测定材料的弹性模量、泊松比；

步骤1.3、建立三维有限元模型；

步骤2：使用p型有限元法计算在对应荷载和约束条件下的线性解中的应力场；

步骤3：通过所求得的线性解中的应力场，计算所对应的几何刚度矩阵；再使用几何刚度矩阵进行屈曲载荷因子特征值求解：具体为

步骤3.1、通过计算得到的应力场求得几何刚度矩阵；

步骤3.2、将所求得的几何刚度矩阵用于特征值计算，将广义的特征值问题Kφ＝λMφ转化为标准特征值问题Kφ＝λφ，即将M化为单位阵；

其中K为刚度矩阵，M质量矩阵，φ为特征矢量，λ即为所求得的屈曲载荷因子特征值；

采用Lanczos迭代方法进行特征值计算；标准特征值问题：Kφ＝λφ

任何初始向量U,||U||_k＝1，U₀＝0，其中|| ||_k为k范数

采用三项递推公式：{U_k+1}＝(K{U_k}-α_k{U_k}-β_k{U_k-1})/β_k+1

其中，β₁＝0

α_k＝{U_k}^TK{U_k}

β_k+1＝||K-{U_k}-α_k{U_k}-β_k{U_k-1}||₂

这里k＝1,2,3,...,m-1≤n，|| ||₂代表2范数，m为最高阶特征值个数，n为方阵阶数，完成迭代计算后，最终得到两组解：用来建立预屈曲应力状态的线性解和特征值屈曲载荷因子解；

步骤4：根据所求得的屈曲载荷因子是否满足精度要求，并判断能量范数误差是否满足一定的要求，如屈曲载荷因子未满足精度要求或者能量范数误差过大，提高插值多项式的阶数并增加网格剖分数量，返回步骤2；屈曲载荷因子满足精度要求，根据屈曲载荷因子获得相应屈曲形状。

2.根据权利要求1所述的使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法，其特征在于：所述步骤2中p型有限元法在单元数量确定的情况下，通过提升插值多项式的阶数来提高计算结果的精度，插值多项式的阶数应大于等于3。