CN104281730A - 一种大转动变形的板壳结构动响应的有限元分析方法 - Google Patents

Abstract

本发明公开了一种大转动变形的板壳结构动响应的有限元分析方法，属于计算非线性结构动力学领域。该方法包括以下步骤：1、在计算机辅助设计软件建立板壳结构的几何模型，然后导入到计算机辅助工程软件中进行网格划分；2、建立计入了大转动变形效应的板壳单元分析模型；3、引入时间积分算法和迭代技术，建立板壳结构动力学响应分析的计算分析模型；4、编制有限元分析程序；5、配置计算环境，进行板壳结构动力学响应的计算。本发明为涉及大转动变形板壳结构的动响应计算提供了一种有效的计算分析工具。本方法对于不受外力作用时的结构系统，不仅能够守恒结构能量，还能守恒结构动量，弥补了现有分析技术的不足。

Description

一种大转动变形的板壳结构动响应的有限元分析方法

技术领域

本发明涉及一种大转动变形的板壳结构动响应的有限元分析方法，属于计算非线性结构动力学领域。

背景技术

板壳结构以其优异的力学性能被广泛地应用于各结构工程中，由于它们在承受外部载荷时常伴随着大位移、大转动的非线性动力学现象，发展高效实用的分析技术以计算它们的非线性动响应对其设计及振动控制具有重要的研究意义和广泛的应用价值。

有限元分析是一种计算机辅助的求解方法，采用数值模拟技术来寻找偏微分方程的近似解，所述偏微分方程表述了结构受载下的力学规律。随着计算机技术的不断进步，有限元分析已经变成了辅助工程师进行结构设计时作决策的重要工具。 

对于涉及大转动变形的板壳结构的非线性动响应问题，结构中存在转角幅度超过1弧度（57.3°）的大转动变形，由于采用了小转动参数描述结点转动，导致基于传统的线性有限元分析方法和发展成熟的几何非线性有限元分析方法（完全的和更新的拉格朗日法）的商用有限元软件均无法准确地对此类问题进行求解。 

近年来，科学家提出了共旋有限元分析方法，由于该方法采用转动矩阵表述结点转动，使其具有大转动变形分析能力，但该方法目前主要应用于静力学分析中，在动力学响应分析中应用较少。在期刊《中国科学》2013年第43卷第1期的论文“薄壳结构非线性动力学响应的共旋有限元能量守恒与衰减算法”中，该文作者提出了一种涉及大转动变形的薄壳结构动响应分析的共旋有限元分析方法，建立了用于求解动响应的能量守恒与衰减算法。然而，对于不受外力作用的结构系统，除了结构能量满足守恒定律外，结构的线动量与角动量守恒也满足守恒定律，而在该论文中，这一重要守恒特性则没有得到体现。 

发明内容

本发明针对现有技术的不足，提出了一种大转动变形的板壳结构动响应的有限元分析方法。本方法对于不受外力作用时的结构系统，不仅能够守恒结构能量，还能守恒结构动量。

本发明为解决其技术问题采用如下技术方案： 

一种大转动变形的板壳结构动响应的有限元分析方法，它包括以下步骤：

        步骤1、壳结构的几何建模及有限元网格划分；

包括在计算机辅助设计软件中根据测量得到的板壳结构几何尺寸，建立结构几何模型，导入到网格划分软件中进行三角形网格划分，对每个单元的材料密度、弹性模量和泊松比及厚度进行赋值，生成单元结点坐标列表和单元—结点连结关系与材料属性、厚度数据列表；

        步骤2、建立计入了大转动变形效应的板壳单元分析模型；

包括确定单元惯性力矢量和内力矢量计算表达式，分别对单元惯性力矢量和内力矢量进行线性化，得到单元切向惯性矩阵和单元切向刚度矩阵，建立单元动力学平衡方程；

步骤3、建立板壳结构动响应分析的计算模型；

包括引入时间积分算法和迭代技术改写步骤2得到的单元动力学平衡方程，采用标准的有限元组集过程得到结构的迭代位移方程；

步骤4、编制有限元分析程序；

包括编写以下程序模块：矢量与矩阵运算模块、单元分析模块、结构组集模块、线性代数方程组求解模块以及主程序模块；

步骤 5、配置计算环境，进行板壳结构动响应的计算；

包括在计算程序中读入步骤1中生成的单元结点坐标列表、和单元—结点连结关系与材料属性、厚度数据列表；施加外载荷，引入边界条件和初始条件；设置时间步长和时间积分常数，进行响应计算，处理计算结果。

步骤2所述的大转动变形是指板壳结构中存在转角幅度超过1弧度的变形。 

步骤2所述的建立计入了大转动变形效应的板壳单元分析模型，是指采用了大转动运动的如下描述手段推导得到的：单元中面内任意点的总转动采用正交转动矩阵表述、单元中面内任意点处的角速度与角加速度是材料形式的。 

步骤2所述的单元内力矢量是指采用了三角形壳单元共旋列式理论建立的。 

步骤3所述的引入时间积分算法和迭代技术，指广义能量动量算法和牛顿-拉夫逊迭代技术。 

与现有技术相比，本发明具有以下有益效果： 

1）为涉及大转动变形板壳结构的动响应计算提供了一种有效的计算分析工具。

2）本发明用于求解板壳结构在不受外力作用时的非线性动响应时，得到的结构动响应除了能够守恒结构能量，还能守恒结构线动量与角动量，很好地符合重要的物理守恒定律。 

3）本发明不仅可用于涉及大转动变形的板壳结构的非线性动响应分析；也可用于仅存在小转动变形的板壳结构的线性动响应分析。 

附图说明

图1 为本发明实施流程图。

图2(a)为顶部受到阶跃集中载荷作用下的球盖壳的网格图；图2(b) 为顶部受到阶跃集中载荷作用下的球盖壳的剖面图。 

图3(a)、图3(b)、图3(c)、图3(d)依次为球盖壳在时刻 、、、的变形构形图。 

图4为球盖壳顶点A处的挠度与球盖壳高度之比的时间变化曲线图。 

图5为矩形薄壳的自由运动算例示意图。 

图6为矩形薄壳在不同时刻的变形构形图。 

图7为矩形薄壳结构的线动量随时间变化曲线图。 

图8为矩形薄壳结构的角动量随时间变化曲线图。 

图9为矩形薄壳结构的能量随时间变化曲线图。 

图10(a)为顶部开孔球盖壳的动态屈曲算例的网格图；图10(b)为顶部开孔球盖壳的剖面图。 

图11(a)、图11(b)、图11(c)、图11(d)、图11(e)、图11(f)、图11(g)、图11(h)、图11(i)、图11(j)依次为顶部开孔球盖壳在时刻、、、、、、、、的变形构形图。 

图12为顶部开孔球盖壳点A处的位移随时间变化曲线图。 

具体实施方式

下面结合附图对本发明创造做进一步详细说明。

本发明的实施流程如图1所示 

步骤1，板壳结构的几何建模及有限元网格划分；

在计算机辅助设计（CAD）软件中根据测量得到的板壳结构几何尺寸，建立结构几何模型，导入到网格划分软件中进行三角形网格划分，对每个单元的材料属性（包括密度、弹性模量和泊松比）及厚度进行赋值，生成单元结点坐标列表和单元—结点连结关系与材料属性、厚度数据列表。

步骤2，建立计入了大转动变形效应的板壳单元分析模型； 

对于任意的三角形平壳单元，采用如下的运动描述方法：矢量单元中面内任意点的总转动采用正交转动矩阵表述、单元中面内任意点处的角速度与角加速度是材料形式的。

利用动能变分建立单元惯性力矢量表达式： 

                              $f_{\mathrm{iner}}=\left\{\begin{array}{c}{\∫}_A{N}_1{M}_{\mathrm{\ρ}}\stackrel{\·\·}{u}\mathrm{dA}\\ {\∫}_A{N}_{1}R(I_{\mathrm{\ρ}}^0\stackrel{\·}{\mathrm{\Ω}}+\widetilde{\mathrm{\Ω}}{I}_{\mathrm{\ρ}}^0\mathrm{\Ω})\mathrm{dA}\\ {\∫}_A{N}_2{M}_{\mathrm{\ρ}}\stackrel{\·\·}{u}\mathrm{dA}\\ {\∫}_A{N}_{2}R(I_{\mathrm{\ρ}}^0\stackrel{\·}{\mathrm{\Ω}}+\widetilde{\mathrm{\Ω}}{I}_{\mathrm{\ρ}}^0\mathrm{\Ω})\mathrm{dA}\\ {\∫}_A{N}_3{M}_{\mathrm{\ρ}}\stackrel{\·\·}{u}\mathrm{dA}\\ {\∫}_A{N}_3(I_{\mathrm{\ρ}}^0\stackrel{\·}{\mathrm{\Ω}}+\widetilde{\mathrm{\Ω}}{I}_{\mathrm{\ρ}}^0\mathrm{\Ω})\mathrm{dA}\end{array}\right\}---\left(1\right)$                       

式中，为三角形壳单元的面积，为线性插值函数，，为壳单元材料密度，为单元厚度，、、、及分别为壳单元中面内任意点的平动加速度、转动矩阵、转动惯性张量、材料形式的角速度及角加速度。

基于共旋有限元方法确定单元内力矢量表达式： 

                                                                             

式中，为单元共旋转换矩阵，为单元局部内力矢量，为局部坐标系中的三角形线性平壳单元的材料刚度矩阵，为经共旋列式思想提取出的单元局部纯弹性变形位移矢量，上标T为矩阵转置符号。

 对单元惯性力矢量进行线性化处理，得到单元切向惯性矩阵： 

                                                                                  

式中， 为质量矩阵，为陀螺矩阵，为离心力刚度矩阵，与为纽马克时间积分参数，为计算时间步长。

对单元内力矢量进行线性化处理，得到单元切向刚度矩阵： 

                                                                                                          

式中，为材料刚度矩阵矩阵，为几何刚度矩阵。

根据上面确定单元惯性力及内力矢量表达式，建立单元动力学平衡方程： 

                                                                                   

式中为单元外力矢量。

步骤3，建立板壳结构动力学响应分析的计算模型； 

引入广义能量—动量算法积分算法和牛顿-拉夫逊迭代技术改写步骤2得到的单元动力学平衡方程，为：

                                                               

式中，与为时间积分算法中的积分常数， 、及分别为单元广义惯性力、内力及外力矢量，它们的计算表达式为：

                 

式中， 为广义平动加速度矢量，与分别为材料形式的广义角速度与角加速度矢量，为广义转动矩阵，为单元广义共旋矩阵，为单元广义局部内力矢量，为单元广义局部变形位移矢量，与分别为时刻和的单元结点外力矢量，它们的计算表达式如下：

                                                             

式中，与分别为时刻和的平动加速度，与分别为时刻和的转动加速度，与分别为时刻和的转动角速度，与分别为时刻和的转动矩阵，与分别为和的共旋转换矩阵，与分别为时刻和的局部纯弹性变形位移矢量。

仅对单元广义局部变形位移矢量的计算引入额外的数值阻尼，将的计算表达式改写为： 

                                                   

式中，为额外引入的数值阻尼系数，。此时，单元广义动力学平衡方程为：

                                                      

式中，为额外引入数值阻尼系数后的单元广义内力矢量。

将单元广义动力学平衡方程在第次迭代处作一阶Taylor级数展开进行线性化，可得： 

                       

式中

                                                             

式中，为单元迭代位移矢量，单元广义切向惯性矩阵与广义切向刚度矩阵分别为：

                        

式中，、及分别为单元质量矩阵、陀螺矩阵及离心刚度矩阵，与分别为单元广义材料刚度矩阵和几何刚度矩阵，与为纽马克公式的积分常数，为时间步长。将、、及采用谱半径进行表示：

                                               

采用上式表示后，当取时，得到的算法为能量动量守恒算法，当取时，得到的算法为能量衰减动量守恒算法，这两种算法合在一起称作广义能量动量算法。

将方程(20)(21)代入到方程(19)中，可得单元位移迭代方程为： 

                                                    

式中，为单元广义切向矩阵，为单元残余力矢量。

将单元位移迭代方程在整个结构上进行遍历与组集，得到结构位移迭代方程为： 

                                                                 

式中，为结构广义切向矩阵，为结构迭代位移矢量，为结构残余力矢量，符号“”表示标准的有限元分析组集过程。

步骤4， 编制有限元分析程序； 

该程序主要包括以下模块：

1）  矢量与矩阵计算模块，实现程序中涉及的矢量与矩阵运算；

2）  单元分析模块，实现三角形壳单元的惯性力、内力矢量以及切向惯性矩阵、切向刚度矩阵的计算；

3）   结构组集模块，实现单元矢量与矩阵向结构矢量与矩阵的组集过程；

4）   线性代数求解器模块，实现结构位移迭代方程的求解；

5）   主程序模块，实现配置计算环境，控制计算流程，调用其它模块完成计算。

步骤5，配置计算环境，进行板壳结构动力学响应的计算分析； 

将步骤1中生成的单元结点坐标列表和单元—结点连结关系与材料属性、厚度数据列表读入到步骤4编制的程序中，施加外载荷，引入边界条件和初始条件，设置时间步长及积分参数等；运行程序进行响应计算分析；计算完成后进行计算结果处理。

 实施实例 

下面结合具体的实例，对本发明的效果进行阐述。应该指出，此处所描述的具体实例仅仅用于解释本发明，并不用于限定本发明。此外，需要说明的是本发明的保护范围不限于下述的计算实例。

实例1：固支球盖壳的动力学响应 

球盖壳的几何形状及剖面如图2所示，几何尺寸为：半径，角度，厚度，球盖壳高度。材料弹性模量、泊松比及密度分别为、和。球盖壳周边固支，在其顶点A处受到阶跃集中力的作用，其中。利用对称性，选取球盖壳的四分之一结构进行建模，单元数量为200。时间步长为，计算时间为。选用的时间积分算法分别为能量动量守恒算法、能量衰减动量守恒算法1()能量衰减动量守恒算法2()。球盖壳在阶跃集中力的作用下，在时刻、、、的变形构形依次如图3(a)、图3(b)、图3(c)、图3(d)所示。点A处z方向上的挠度与球盖壳高度之比随时间变化曲线如图4所示。

实例2：矩形薄壳的自由运动 

矩形薄壳如图5所示，几何尺寸为：长度，宽度，厚度。材料弹性模量、泊松比及密度分别为、和。矩形薄壳不受约束，在图中所示三个位置的各结点上受到外力的作用，其中，的表达式为

                                             

计算单元数量为，时间步长为，计算时间为。所选用的时间积分算法为选用的时间积分算法分别为能量动量守恒算法、能量衰减动量守恒算法1()能量衰减动量守恒算法2()。

在外力的作用下，矩形薄壳产生了大的三维弯曲及扭转变形，同时伴随大的刚体运动， 在时间段内不同时刻的变形构形如图6所示，其中时间间隔为。结构的线动量与角动量分别见图7和8，可以看出，在外力为零后，能量动量守恒算法与两种能量耗散动量守恒算法均精确地守恒了结构的线动量与角动量。结构能量变化曲线如图9所示，由于在能量动量守恒算法中未引入任何数值阻尼，故而不会对结构能量产生耗散作用，算法精确地守恒了结构能量，而两种能量耗散算法均对结构能量产生了耗散作用。 

实例3：顶部开孔球盖壳的动态屈曲 

顶部开孔的球盖壳几何形状及剖面如图10所示，几何尺寸为：半径，角度，，厚度。弹性模量、泊松比、密度分别为、和。该球盖壳的底部边在z方向上的位移设置为零。顶部边的各结点上受到沿z轴向下的结点力作用，所采用的计算单元数量为，顶部边各结点上的作用力的大小随时间变化的关系式为

                                                                               

计算采用的时间步长为，计算时间为，所采用的时间积分算法为能量动量守恒算法、能量衰减动量守恒算法1()能量衰减动量守恒算法2()。  在外力由零线性增至恒定值的过程中，球盖壳经历了前屈曲—屈曲—后屈曲的动态变形过程，在时刻、、、、、、、、的变形构形依次如图11(a)、图11(b)、图11(c)、图11(d)、图11(e)、图11(f)、图11(g)、图11(h)、图11(i)、图11(j)所示。初始位置为在的A点处的位移变化曲线如图12所示。

Claims (5)

1.一种大转动变形的板壳结构动响应的有限元分析方法，其特征在于，它包括以下步骤：

     步骤1、壳结构的几何建模及有限元网格划分；

包括在计算机辅助设计软件中根据测量得到的板壳结构几何尺寸，建立结构几何模型，导入到网格划分软件中进行三角形网格划分，对每个单元的材料密度、弹性模量和泊松比及厚度进行赋值，生成单元结点坐标列表和单元—结点连结关系与材料属性、厚度数据列表；

     步骤2、建立计入了大转动变形效应的板壳单元分析模型；

包括确定单元惯性力矢量和内力矢量计算表达式，分别对单元惯性力矢量和内力矢量进行线性化，得到单元切向惯性矩阵和单元切向刚度矩阵，建立单元动力学平衡方程；

步骤3、建立板壳结构动响应分析的计算模型；

包括引入时间积分算法和迭代技术改写步骤2得到的单元动力学平衡方程，采用标准的有限元组集过程得到结构的迭代位移方程；

步骤4、编制有限元分析程序；

包括编写以下程序模块：矢量与矩阵运算模块、单元分析模块、结构组集模块、线性代数方程组求解模块以及主程序模块；

步骤 5、配置计算环境，进行板壳结构动响应的计算；

包括在计算程序中读入步骤1中生成的单元结点坐标列表、和单元—结点连结关系与材料属性、厚度数据列表；施加外载荷，引入边界条件和初始条件；设置时间步长和时间积分常数，进行响应计算，处理计算结果。

2.根据权利要求1所述的一种大转动变形的板壳结构动响应的有限元分析方法，其特征在于，步骤2所述的大转动变形是指板壳结构中存在转角幅度超过1弧度的变形。

3.根据权利要求1所述的一种大转动变形的板壳结构动响应的有限元分析方法，其特征在于，步骤2所述的建立计入了大转动变形效应的板壳单元分析模型，是指采用了大转动运动的如下描述手段推导得到的：单元中面内任意点的总转动采用正交转动矩阵表述、单元中面内任意点处的角速度与角加速度是材料形式的。

4.根据权利要求1所述的一种大转动变形的板壳结构动响应的有限元分析方法，其特征在于，步骤2所述的单元内力矢量是指采用了三角形壳单元共旋列式理论建立的。

5.根据权利要求1所述的一种大转动变形的板壳结构动响应的有限元分析方法，其特征在于，步骤3所述的引入时间积分算法和迭代技术，指广义能量动量算法和牛顿-拉夫逊迭代技术。