Bush 有限单元原理及其在航天器结构建模中的应用

2010-01-08邹元杰

航天器工程 2010年1期

邹元杰

(北京空间飞行器总体设计部,北京 100094)

1 引言

航天器结构或机构中存在大量的杆系结构和连接机构[1-2],考虑这些结构的特点,在力学分析中有时采用弹簧元模拟比较方便。航空航天领域常用的有限元分析软件M SC.N AS TRAN,既提供了标量弹簧单元(Scalar Spring),也提供了6 个自由度广义弹簧单元(Bush)。其中,标量弹簧单元刚度系数物理意义明确,应用起来比较容易。而广义弹簧单元Bush 的单元原理和参数含义尚不明确。

作者在使用该单元的过程中遇到诸多问题,难以合理解释。首先,按照M SC.NAS TRAN 的用户手册,Bush 元的刚度阵K0(对应单个节点的6 个自由度,同下文(5)式中的Κ11或K22)形式为[3]:

其中, ku、kv、kw分别代表三个平动自由度的刚度系数, kθx、kθy、kθz分别代表三个转动自由度的刚度系数。按照这个解释,Bush 元的刚度阵(对应单个节点)其实是对角阵,其6 个自由度显然是不耦合的。而实际计算发现,在Bush 元节点上施加垂直于单元轴线方向的力,会产生弯曲方向的转角;同样,施加弯矩,也会在垂向产生平动位移。这说明,6个自由度的刚度并非解耦,即Bush 单元刚度阵不可能如(1)式所示,必然存在耦合项。其次,工程上对于杆系和连接机构建模,也常采用梁单元(主要包括欧拉梁元和剪切梁元),然而梁单元与Bush 元的对应关系尚不清楚。再次,在计算中发现,标量弹簧单元的刚度仅与输入的刚度系数有关,其长度可以为0,即定义单元的两个节点可以重合,但Bush 元的长度取为0 则出错。

欲解开上述困惑,唯一的途径是对Bush 元的原理进行研究。为此,我们开展了以下研究:1)设法找到Bush 元刚度阵与输入的6 个刚度系数的关系;2)搞清楚6 个刚度系数的真实物理意义,进而分析Bush 元与标量弹簧单元、梁单元的关系;3)结合工程算例,对上述Bush 有限元原理进行验证。

2 Bush 单元刚度阵表达式的获取

由于大多数商用软件对用户来说都是“黑箱子”, 很难确切知道内核程序和算法, 因此, 得到NAS TRAN 中Bush 元刚度阵的具体表达式是非常困难的。我们只有通过大量的计算、比对,并结合相关的有限元理论知识,由数值结果反推单元刚度阵的表达式。利用DMAP 语言,我们可以输出Bush元的刚度阵,然后,通过改变输入参数,观察刚度阵的变化,进而获取刚度阵的解析表达式。

图1 Bush 单元示意图Fig.1 Sketch of Bush element

对于图1 所示的Bush 单元,其两端节点的位移向量q 和对应的力向量R 为

单元的静力平衡方程为

其中,单元刚度阵K 可以按两个节点分解为4个分块矩阵

设Bush 元的6 个输入参数为ki(i =1, ...,6), 单元长度为L 。通过数值计算最终反推得到的各分块矩阵解析表达式如下:

上述表达式已经通过一些算例验证。从公式(6)～(8)以及相关算例的计算结果可以看出:1)对应单个节点的刚度阵(如K11)并非对角阵,存在垂向位移和转角的耦合刚度项,而且6 个对角线元素也并非对应输入参数ki(i =1, ...,6)。2)对于纯拉压或纯扭转问题,刚度阵中只有k1 或k4 起作用,此时,采用标量弹簧元和Bush 元建模是等效的。而对于更复杂的力学问题,由于存在耦合项,Bush元的刚度不可能由6 个相互独立的标量弹簧元(分别对应节点6 个自由度的刚度)来代替,二者是无法等效的。3)Bush 元刚度阵元素不仅和输入的刚度系数ki(i =1, ...,6)有关,而且是Bush 单元长度L 的函数。由于Nastran 程序默认Bush 元需要考虑耦合项影响,因此,耦合项不能为0(将k3或k2取为0 的情况除外),即式(6)～(8)中的L 不能取为0。而标量弹簧元的刚度阵与单元长度无关,因此,单元长度L 可以取为0,即两端节点可以重合。

3 Bush 单元与梁单元的对应关系

为了找到Bush 单元与梁单元的关系,对比了两类单元的刚度阵表达式,试图用梁单元的输入参数来表示Bush 元的输入参数。

梁单元的相关研究多年来持续不断[4-9]。一般单元刚度阵的推导都应用虚位移原理,但是对于梁,应用结构力学中的解析方法(直接法)推导更为方便,只要列出梁单元两端的位移要素与外力之间的关系即可求出。剪切梁(Timoshenko 梁)单元的刚度阵如下(节点位移和外力的定义同Bush 元,见图1)[10]:

在式(10)～(12)中, E 为材料的弹性模量, G为材料的剪切模量, A 为截面面积, L 为梁单元的长度, Iy、Iz分别为绕y 、z 轴的弯曲惯性矩, J 为扭转惯性矩,φy、φz分别为y 向和z 向的剪切修正系数。

比较式(10)～(12)和式(6)～(8)可以看出,Bush 元刚度阵的非零元素和剪切梁单元是一一对应的,说明二者所表征的物理含义有可能相同。将式(9)和式(5)逐项对应,可以将Bush 元的刚度参数用剪切梁的参数表达。

Bush 元各刚度系数与剪切梁参数的对应关系为

若将Bush 元的刚度参数表达为(13)～(18)式,则式(9)和式(5)中的刚度阵完全一样。进一步分析可以看出,两个刚度阵各有7 个独立参数表征:在Bush 元的刚度阵中,7 个参数分别为ki(i =1, ...,6)和L ;在剪切梁元中,7 个参数为EA 、EI y 、EIz 、GJ 、φy 、φz 和L 。虽然用于表达的参数不同,但Bush 元和剪切梁元的刚度阵却可以完全等价。

Bush 元各参数的物理意义也不难理解。k1 、k4分别是拉压刚度(轴向)和扭转刚度,这两个参数的测量也比较容易。k5、k6分别是z 向和y 向的弯曲刚度,只要使梁处于纯弯曲状态,这两个参数容易测出。k2、k3的物理意义解释稍为复杂,这里以k2为例加以说明。如图1 所示,对于Bush 单元i～j ,若i 端固支, j 端约束绕z 轴转角(θzj=0),同时施加y 向外力Fyj,此时, j 端y 向位移为uyj。由式(4)可以看出:也就是说, k2 对应悬臂梁自由端转角被约束时自由端的垂向(y 向)刚度。如果要直接测得k2或k3,则需要约束悬臂梁自由端的转角,测量垂向载荷和垂向位移,进而换算得到k2 或k3 。

对于欧拉梁,由于忽略其剪切效应影响,剪切修正系数φy、φz应取为0。此时,其对应的Bush 元6个刚度参数中, k2 与k6 、k3 与k5 是相互关联的。所以,与欧拉梁对应的Bush 元只有四个刚度参数是独立的。如果采用试验方法测得6 个刚度系数,则k2与k6(或k3与k5)可以相互校验测量结果的准确性。

4 Bush 单元的工程应用

本节针对典型的航天器结构, 应用Bush 单元进行建模和分析。通过对比采用Bush 元的分析结果与采用梁单元的计算结果,可以进一步验证Bush元参数与剪切梁参数对应关系的准确性。

4.1 太阳翼结构的分析

对太阳翼展开状态的动力学分析而言,铰链的建模十分关键。通常工程上采用梁单元进行建模。现改用Bush 元建立铰链模型,替换原模型中的梁单元,而其它结构单元属性均不改变。利用式(13)～(18),由梁单元参数确定Bush 元参数,研究计算结果的变化。由于Bush 单元没有质量属性,采用Bush 元时, 需要在节点增加集中质量(Lumped mass)单元,集中质量大小取为梁单元质量的一半。动力分析时,两种模型的质量阵均采用集中质量处理。图2 为某太阳翼的有限元模型。

图2 太阳翼有限元模型Fig.2 Finite element model for solar array

首先计算静力响应。将太阳翼根铰与太阳翼驱动机构(SADA)连接处固定支撑,在整个结构上施加1 gn的惯性加速度载荷(垂直于面板的方向,即Z向),计算太阳翼最外端,节点898 处的静力响应。对铰链分别采用Bush 元和剪切梁元建模时节点898 处的静力响应对比情况见表1。从表中数据看,两者的计算结果是完全一致的。T1, T2, T3分别对应于沿X、Y 、Z 轴方向的线位移,R1, R2,R3分别对应于绕X 、Y 、Z 轴的转角。

表1 太阳翼的静力分析结果Table 1 Static analysis results for solar array

然后,进行模态分析。边界条件同静力分析。分析结果表明:两种方法计算的各阶模态频率非常接近(见表2),振型也完全一致。

4.2 桁架结构的分析

下面计算一个更为复杂的桁架结构。针对某大型桁架结构(见图3),对所有的连接杆结构分别采用梁单元和Bush 元建模,对比两种方法的静力和模态分析结果。该桁架结构需要采用2 816 个Bush 元或梁单元建模。

表2 太阳翼的模态分析结果Table 2 Modal analysis results for solar array

首先在桁架左端固支,在右端某节点加静载荷(Z 向)1N,采用两种单元处理方法计算加载点的静力响应。加载点的静力响应见表3。T1,T2,T3分别对应于沿X、Y 、Z 轴方向的线位移,R1,R2,R3分别对应于绕X 、Y 、Z 轴的转角。从数据看,Bush 元的计算结果与剪切梁单元非常接近。另外,进行了桁架结构的模态分析,分析结果表明:两种方法计算的各阶模态频率(见表4)和振型均吻合较好。