基于MITC单元的加筋板动力稳定性分析*

2021-03-05裴志勇吴卫国

武汉理工大学学报（交通科学与工程版） 2021年1期

李念裴志勇吴卫国

(武汉理工大学交通学院1) 武汉 430063)(武汉理工大学绿色智能江海直达船舶与邮轮游艇研究中心2) 武汉 430063)

0 引言

为了确保结构具有足够的稳定性，临界载荷的计算非常重要.当所考虑的结构承受的载荷是动载荷时，则属于动力稳定性问题，在该问题上，由于载荷随着时间而变化，使得动力稳定性问题要比静力的情况复杂很多.

在弹性范围内，通常将动力稳定性问题转化为一个Mathiue方程进行求解，有限元方法也应用到动力稳定性问题中，Thana等[1]采用有限元法分析讨论了梁结构在轴向周期荷载作用下的动力稳定性问题.Loja等[2]研究了可变刚度的复合材料板的动稳性，并考虑了静态和动态载荷因数，以及边界条件的影响.Talimian等[3]采用有限差分法求解板的运动微分方程，分别从载荷系数、载荷类型，及其组合形式等方面进了讨论，总结了载荷参数对不稳定区的影响规律.相比于平板而言，加筋板由于其板和筋之间相互耦合，使得从解析角度难以求解.目前的研究方法，一类是将其简化为便于处理的计算模型，如正交异性板模型、交叉梁系模型[4-5]，但它们在使用上都具有一定的局限性；另一类方法则是对加筋板进行离散化处理，采用有限元法、有限差分法等数值方法进行求解，刘媛[6]基于Timoshenko梁和Mindlin板理论，建立了加筋板动力稳定性分析的有限元模型，具体分析了加强筋尺寸、种类及其数目对结构动力稳定性的影响.

MITC单元是一种基于假定应变法的张量分量混合插值板单元(mixed interpolation of tensorial components，MITC)，它具有完备的理论基础，不会出现剪切锁闭现象，且对单元的几何扭曲不敏感，能很好地适用于各类板壳结构，在收敛性和通用性方面具备优势[7-8].为了探究在工程上受到广泛应用的加筋板结构的动力稳定性，文中根据拉格朗日方程及哈密尔顿原理得到结构动稳性的矩阵方程，并利用Bolotin方法将动力稳定性问题转换为求解广义特征值问题，采用MITC4壳单元对加筋板进行离散化，推导了四节点24自由度壳单元特性矩阵，并在Matlab平台上开发了求解加筋板屈曲、模态以及动力稳定性的有限元计算程序，以周期性边界条件下加筋板为研究对象进行计算分析，先通过对比MITC4单元以及Abaqus软件的屈曲和模态计算结果，验证了MITC4壳单元的通用性和准确性，并能较好地预报加筋板的动力稳定性，最后考虑了静载系数、截面参数等因素对加筋板的动力稳定性的影响，为结构稳定性的改善提供设计参考.

1 结构动力稳定性控制方程

根据拉格朗日方程L=U+V-T及哈密尔顿原理，有如下关系.

(1)

式中：U为弹性体系的应变能；V为外力势能；T为结构的动能.

(2)

可以得到无阻尼情况下结构的平衡方程：

(3)

式中：M为质量矩阵；K为刚度矩阵；Kg为几何刚度矩阵.

现考虑动载荷为如下形式的周期性载荷，Pcr为结构的静态临界载荷.

P(t)=Ps+Pdcos (θt)=Pcr(α+βcos (θt))

(4)

根据Bolotin方法，形如式(3)带有周期系数的二阶Mathieu微分方程解的收敛性可以被周期为T和2T(T=2π/θ)的解为边界所确定，分别将两类解以傅里叶级数的形式展开.

(5)

将式(4)～(5)带入式(3)，利用正余弦项之间线性独立的关系，令其系数为零可以得到以下线性方程组.

(6b)

(6c)

(6d)

方程组(6)存在非零解的充要条件是其系数矩阵为奇异矩阵，通过求解前n阶行列式可以得到前n个动力不稳定区域，随着阶数的增大，频率区间迅速减小，其中第一阶不稳定区已包括参数平面的绝大部分，也称为主要动力不稳定区，它可以由式(6a)与(6b)的一阶解确定.

(7)

通过求解式(7)矩阵的广义特征值，从而可以得到由α，β，θ等参数决定的动力不稳定区域.

2 MITC4壳单元特性矩阵

2.1 平面应力四节点等参元

单元局部坐标系下，双线性四边形四节点等参元插值函数为

(8)

根据平面应力问题的应变-位移几何关系式，有：

(9)

(10)

则由膜应变对应的单元刚度矩阵Kme为

(11)

式中：ξm,ηk为等参元在自然坐标系下高斯积分点坐标；Wm,Wk为相应的高斯积分点权重系数；J为雅可比矩阵.

数值积分方法选用高斯积分，双线性四边形单元采用2×2阶高斯积分即可达到完全积分效果.

2.2 弯曲单元

目前常用的板弯曲壳单元是基于Mindlin板理论的板单元，这种考虑了板横向剪切应变的单元适合模拟厚板弯曲.但在有限元模拟时，如果划分单元的厚度与边长之比较小时，会导致单元的横向剪切应变被夸大，剪切刚度人为“刚化”，即发生剪切锁闭现象.尽管采用减缩积分或者选择积分等方式能解决剪切锁闭现象，但同时有可能会导致多余的零能模式，使得求解失真.因而，选用一种通用性更好，计算更为精确的单元是有必要的，MITC壳单元能较好地解决剪切自锁现象，且不必对数值积分方案做调整.

弯曲单元位移插值函数采用与3.1节平面应力单元相同的双线性函数：

(12)

由几何方程得到板面内弯曲应变关系：

(13)

(14)

由高斯积分得到单元弯曲刚度矩阵Kbe：

(15)

考虑板横向剪切应变的影响，基于Mindlin板理论，板单元的剪切应变为式(16)，结合插值函数式(12)可得到应变矩阵Bs.

(16)

而基于MITC板理论，为了消除剪切闭锁，使用与Mindlin板单元相同的插值函数，用插值点处的剪切应变值，对插值点处的横向位移、转角等协调应变进行混合插值，从而用假定的剪切应变来替换单元剪切应变部分，达到消除剪切自锁的目的.以单元四边中点的剪切应变进行插值，得到一组假定应变，见图1.

图1 MITC等参元和MITC4壳单元

(17)

(18)

2.3 第六自由度扭转单元

根据壳单元的公式，每个节点的刚度矩阵只耦合5个自由度，然而，采用5自由度壳单元离散加筋板结构存在困难.为此，文中基于Kanok[9]的方法，引入了面内扭转自由度分量.

定义罚能量函数，当壳体发生刚体转动，罚能量变为零，为

(19)

式中：G为材料剪切模量；κT为罚因子，文献[9]推荐取大于0.1的常数来抑制面内扭转刚度，随着κT增大，面内扭转模态的影响会迅速减小并趋于稳定，文中取κT=1.

扭转虚应变

(20)

采用减缩高斯积分计算单元扭转刚度矩阵Kte

(21)

2.4 单元质量矩阵与几何刚度矩阵

单元质量矩阵的推导类似于刚度矩阵，通过Hamilton原理即可得到单元质量矩阵的一般表达式.通常，质量矩阵分为一致质量阵和集中质量阵，考虑到动力稳定性问题会涉及到频率计算，文中采用在频率计算方面更为精确的集中质量矩阵.

(22)

广义几何应变：

(24)

由虚功原理得到几何刚度矩阵：

(25)

3 计算结果及分析

3.1 加筋板模型

现考虑三跨加筋板模型，图2.a=2 550 mm，b=850 mm，材料的弹性模量E=2.06×1011Pa，泊松比ν=0.3，密度ρ=7 850 kg/m3，受到纵向周期载荷的作用.约束条件为周期性边界条件，即在加筋板的四条边中，对边分别具有相同的位移模式，加载端和其中一条邻边用MPC约束，加筋板是构成大型结构的基本单元，采用周期性边界条件能较好地反映周边的影响，从而能够用部分的性质来表达整体的性质.加强筋尺寸见图3，h×bf×tw/tf=600 mm×150 mm×15/20 mm，板和加强筋均用壳单元进行离散，网格大小为50 mm.

图2 加筋板模型的尺寸和边界条件

图3 加强筋截面尺寸

3.2 屈曲及模态分析

表1 屈曲及模态(一阶)对比

3.3 动力稳定性分析

3.3.1静载系数的影响

取细长比为γ=2.21的加筋板分析计算，由3.2节的屈曲分析得到该加筋板的静态临界力Pcr=4 357 N/mm，动载荷P(t)=Ps+Pdcos (θt)=Pcr(α+βcosθt)，其中静载荷分量Ps=αPcr，为探讨静载分量对加筋板动力稳定性的影响，令α=0，0.2，0.4，0.6，计算不同静载荷分量下加筋板的动力不稳定区.计算结果见图4，在由β和动载荷频率θ构成的参数平面中，随着静载分量系数α的增大，加筋板的动力不稳定区失稳频率值越小，且不稳定频率区间范围变得越大，这说明静载分量越接近临界载荷，加筋板结构对动载荷的响应越敏感，即使在较低频率、较小动载荷幅值的动载荷作用下，也可能会发生失稳现象.

图4 不同静载分量下的动力不稳定区域

3.3.2加筋板截面参数的影响

最小剖面模数是衡量加筋板强度的重要指标，为研究加筋板的截面参数对动稳性的影响，设计了6组最小剖面模数基本相同，而细长比各不相同的加筋板，见表2.分别计算这6组加筋板的动力不稳定区(周期载荷静载分量均取α=0)，计算结果见图5，细长比是反应加筋板柔度的物理量，随着加筋板柔度的增大，失稳频率逐渐降低，同时不稳定区域的范围也逐渐缩小.因此，在最小剖面模数相同的前提下，加筋板的设计需要综合考虑外载荷的幅值和频率范围，选取合适的尺寸，从而提高加筋板的动力稳定性.