有限元软件梁杆单元局部坐标系的简化确定方法

宋 浩， 曹聪慧2， 陈 伟， 郭 进

(1.石家庄铁道大学 土木工程学院，河北 石家庄 050043；2. 河北经贸大学 管理科学与工程学院，河北 石家庄 050061)

有限元软件在工程结构(特别是复杂三维结构)的反应分析中得到日益广泛的运用。梁杆单元具有概念简单、计算效率高和使用方便等特点，是使用最为广泛的单元之一。在梁杆单元的运用过程中，涉及到单元特性、单元荷载等属性的定义，还需对计算结果进行处理，这些都需要用到单元局部坐标系，所以理解和确定梁杆单元局部坐标系尤为重要。

以SAP2000为例，整体坐标系是三维直角坐标系，标记为OXYZ，相互垂直并满足右手准则。在整体坐标系下，每一梁杆单元都具有默认局部坐标系，标记为O′123。对于一般简单单元，默认坐标系即是梁杆单元的最终局部坐标系。但是对于某些具有特殊空间位置或特殊要求的单元，最终局部坐标系与默认局部坐标系并不重合。在SAP2000中，默认局部坐标系中的1轴沿单元长度方向，2、3轴位于和用户指定的单元方向相垂直的平面内。2、3轴的默认方向是根据局部1轴与整体Z轴的关系来确定的：

梁杆单元的最终局部坐标系是默认局部坐标系中2、3轴绕1轴旋转ang角度，旋转方向服从右手准则；默认情况下的ang等于0，即为默认局部坐标系，如图1所示。ang的大小根据杆件的空间位置和方向确定。

最终局部坐标系可由其原点O′和3个坐标轴上的点1′、2′、3′来表示。根据上述最终局部坐标系和默认坐标系的关系可知它们原点重合、轴O′1′与轴O′1也重合。因此点2′、3′在默认局部坐标系中非常容易确定，而其在整体坐标系中的表达需通过坐标转换实现。因此，上述最终局部坐标系的确定实际为坐标转换问题，即将点2′、3′在默认局部坐标系中的坐标转化到整体坐标系。

图1 相对于默认方向的梁杆单元坐标角

对于三维坐标转换问题，传统转换方法是基于泰勒级数展开的线性模型[1-3]和基于罗德格矩阵的三维坐标转换方法[4-7]等。基于泰勒级数展开的方法主要是利用泰勒级数展开的方法将模型线性化，然后解算坐标转换的旋转和平移参数，只有当2个坐标系间的旋转角为小角度时，才能对旋转角参数进行线性近似处理且计算繁琐，对于旋转角较大时采用线性模型会引起较大的模型误差；基于罗德格矩阵方法不需要进行三角函数的计算和迭代运算，利用反对称矩阵和罗德里格矩阵的性质，把传统旋转角参数用反对称矩阵的3个独立元素代替，计算过程相对于泰勒级数展开的方法更为简单，对坐标系间的旋转角大小没有限制，具有更好的适用性。两种方法都是需要求出坐标转换的旋转矩阵R、平移矩阵[ΔX,ΔY,ΔZ]T和尺度因子λ这7个参数。以上方法均涉及到较为复杂的数学变换，不利于理解和应用。而本文针对梁杆单元最终局部坐标系的特性，即默认情况下局部坐标系和最终局部坐标系共原点且一个轴重合，针对这一特殊性，简化传统三维坐标系的转换格式，指出了确定旋转矩阵R的简便方法，从而可较为简单地实现最终局部坐标系在整体坐标系的表示。

图2 2个不同的空间直角坐标系

1 2种三维坐标转换方法

设存在2个空间三维直角坐标系O-XYZ和O′-X′Y′Z′，2个坐标系空间位置如图2所示，直角坐标系O′-X′Y′Z′通过坐标原点的平移、分别以X′、Y′、Z′为旋转轴各旋转εX′、εY′、εZ′角度与坐标O-XYZ重合。

坐标转换关系如下

(1)

式中，λ为2个坐标系的尺度因子;R为旋转矩阵; [ΔX,ΔY,ΔZ]T为平移参数矩阵。

构造旋转矩阵R的方法[9]如下。设坐标系O′-X′Y′Z′绕Y′轴旋转的角度是εY′,绕X′轴旋转的角度是εX′，绕Z′轴旋转角度为εZ′。

整个旋转矩阵以各旋转轴旋转顺序对应的旋转矩阵右乘，就可以得到旋转矩阵R。

(2)

由式(1)、式(2)，求出ΔX、ΔY、ΔZ、εX′、εY′、εZ′、λ7个未知参数，就可以实现2个坐标系之间的相互转换。

(3)

将公共点的重心化坐标代入式(1), 可得

(4)

(5)

因而，转换参数可分两步来求解, 即先用式(4)求出旋转参数和比例因子, 再用式(5)求出平移参数。

1.1 传统方法基于罗德里格矩阵的转换

式(4)两边取2-范数, 由于λ> 0 及旋转矩阵为正交阵的特性，可得

(6)

对于n个公共点, 可得λ的最小均方估计

(7)

得到比例因子的最小均方估计后, 利用罗德里格矩阵的性质[10], 可将旋转矩阵R表示为

R=(I-S)-1(I+S)

(8)

式中，I为单位矩阵，S为反对称矩阵。

设

(9)

式中，a,b,c为罗德里格参数。将式(8)、式(9)代入式(4)得

(10)

展开得

(11)

因此，对于n个公共点, 可根据式(11)列出如下的总体误差方程[11]

V3n×1=A3n×3X3×1-L3n×1

(12)

利用式(12)根据最小二乘原理, 无需迭代即可直接求得罗德里格参数

(13)

求得罗德里格参数后，可按式(8)求得旋转矩阵R，然后再根据解算出的比例因子和旋转参数, 按式(5) 可求得平移参数。

图3 最终局部坐标系和默认局部坐标系的相对空间位置

1.2 基于本文简化方法的转换

此时的最终局部坐标系与默认情况下坐标系的坐标的转换关系为

(14)

即

(15)

显然旋转矩阵为

(16)

旋转矩阵R即为默认局部坐标系下1、2、3轴的坐标构成的矩阵。其中1轴为沿梁杆单元长度方向，后两轴位于和用户指定的单元方向相垂直的平面内。由此可见，本文简化方法避免了复杂的矩阵变换，直接可得出旋转矩阵。

2 结果验证

为了验证实验数据的正确性，现分3种情况验证，如表1所示。

表1 验证实验数据 m

设定梁杆单元不平与平行于3个轴，以1轴任意为例求得旋转矩阵、尺度因子、平移参数：

取ang角分别为 -30°，30°，45°，60°，90°，120°，易知3 ′在默认局部坐标系下的坐标3 ′=[0,-sin(ang),cos(ang)]。分别以传统方法基于罗德里格矩阵转化方法和本文简化方法求出3 ′在整体坐标系的坐标，如表2所示。

表2 传统转化方法和本文简化转化方法下的局部3 ′轴的转化结果 m

从验算结果可以看出，本简化方法和传统方法基于罗德里格矩阵转化方法计算结果完全一致，表明本文方法是准确和有效的。

3 结论

针对梁杆单元的三维坐标转换问题，列出了基于罗德里格矩阵的三维坐标转换方法，计算过程繁琐，不利于理解和应用。在SAP2000、Etabs、和Midas等有限元软件的梁杆单元的最终局部坐标系与整体坐标系的对应过程中，提出了利用默认情况下的梁杆单元1、2、3轴的坐标组成转换矩阵，从而简便地实现最终局部坐标系在整体坐标系中的表示。由表2可以看出，基于本文简化方法和基于罗德里格矩阵转化方法计算的最终局部坐标系中局部3 ′轴在整体坐标系的坐标，两种方法计算的结果完全一致，具有足够的精度，表明本文方法是准确有效的。本文简化方法理论严密，计算简单，无小角度限制，简化了旋转矩阵的求解过程，提高了运算效率。