有限元軟件梁桿單元局部坐標系的簡化確定方法

宋 浩， 曹聰慧2， 陳 偉， 郭 進

(1.石家莊鐵道大學(xué) 土木工程學(xué)院，河北 石家莊 050043；2. 河北經(jīng)貿(mào)大學(xué) 管理科學(xué)與工程學(xué)院，河北 石家莊 050061)

有限元軟件在工程結(jié)構(gòu)(特別是復(fù)雜三維結(jié)構(gòu))的反應(yīng)分析中得到日益廣泛的運用。梁桿單元具有概念簡單、計算效率高和使用方便等特點，是使用最為廣泛的單元之一。在梁桿單元的運用過程中，涉及到單元特性、單元荷載等屬性的定義，還需對計算結(jié)果進行處理，這些都需要用到單元局部坐標系，所以理解和確定梁桿單元局部坐標系尤為重要。

以SAP2000為例，整體坐標系是三維直角坐標系，標記為OXYZ，相互垂直并滿足右手準則。在整體坐標系下，每一梁桿單元都具有默認局部坐標系，標記為O′123。對于一般簡單單元，默認坐標系即是梁桿單元的最終局部坐標系。但是對于某些具有特殊空間位置或特殊要求的單元，最終局部坐標系與默認局部坐標系并不重合。在SAP2000中，默認局部坐標系中的1軸沿單元長度方向，2、3軸位于和用戶指定的單元方向相垂直的平面內(nèi)。2、3軸的默認方向是根據(jù)局部1軸與整體Z軸的關(guān)系來確定的：

梁桿單元的最終局部坐標系是默認局部坐標系中2、3軸繞1軸旋轉(zhuǎn)ang角度，旋轉(zhuǎn)方向服從右手準則；默認情況下的ang等于0，即為默認局部坐標系，如圖1所示。ang的大小根據(jù)桿件的空間位置和方向確定。

最終局部坐標系可由其原點O′和3個坐標軸上的點1′、2′、3′來表示。根據(jù)上述最終局部坐標系和默認坐標系的關(guān)系可知它們原點重合、軸O′1′與軸O′1也重合。因此點2′、3′在默認局部坐標系中非常容易確定，而其在整體坐標系中的表達需通過坐標轉(zhuǎn)換實現(xiàn)。因此，上述最終局部坐標系的確定實際為坐標轉(zhuǎn)換問題，即將點2′、3′在默認局部坐標系中的坐標轉(zhuǎn)化到整體坐標系。

圖1 相對于默認方向的梁桿單元坐標角

對于三維坐標轉(zhuǎn)換問題，傳統(tǒng)轉(zhuǎn)換方法是基于泰勒級數(shù)展開的線性模型[1-3]和基于羅德格矩陣的三維坐標轉(zhuǎn)換方法[4-7]等?；谔├占墧?shù)展開的方法主要是利用泰勒級數(shù)展開的方法將模型線性化，然后解算坐標轉(zhuǎn)換的旋轉(zhuǎn)和平移參數(shù)，只有當2個坐標系間的旋轉(zhuǎn)角為小角度時，才能對旋轉(zhuǎn)角參數(shù)進行線性近似處理且計算繁瑣，對于旋轉(zhuǎn)角較大時采用線性模型會引起較大的模型誤差；基于羅德格矩陣方法不需要進行三角函數(shù)的計算和迭代運算，利用反對稱矩陣和羅德里格矩陣的性質(zhì)，把傳統(tǒng)旋轉(zhuǎn)角參數(shù)用反對稱矩陣的3個獨立元素代替，計算過程相對于泰勒級數(shù)展開的方法更為簡單，對坐標系間的旋轉(zhuǎn)角大小沒有限制，具有更好的適用性。兩種方法都是需要求出坐標轉(zhuǎn)換的旋轉(zhuǎn)矩陣R、平移矩陣[ΔX,ΔY,ΔZ]T和尺度因子λ這7個參數(shù)。以上方法均涉及到較為復(fù)雜的數(shù)學(xué)變換，不利于理解和應(yīng)用。而本文針對梁桿單元最終局部坐標系的特性，即默認情況下局部坐標系和最終局部坐標系共原點且一個軸重合，針對這一特殊性，簡化傳統(tǒng)三維坐標系的轉(zhuǎn)換格式，指出了確定旋轉(zhuǎn)矩陣R的簡便方法，從而可較為簡單地實現(xiàn)最終局部坐標系在整體坐標系的表示。

圖2 2個不同的空間直角坐標系

1 2種三維坐標轉(zhuǎn)換方法

設(shè)存在2個空間三維直角坐標系O-XYZ和O′-X′Y′Z′，2個坐標系空間位置如圖2所示，直角坐標系O′-X′Y′Z′通過坐標原點的平移、分別以X′、Y′、Z′為旋轉(zhuǎn)軸各旋轉(zhuǎn)εX′、εY′、εZ′角度與坐標O-XYZ重合。

坐標轉(zhuǎn)換關(guān)系如下

(1)

式中，λ為2個坐標系的尺度因子;R為旋轉(zhuǎn)矩陣; [ΔX,ΔY,ΔZ]T為平移參數(shù)矩陣。

構(gòu)造旋轉(zhuǎn)矩陣R的方法[9]如下。設(shè)坐標系O′-X′Y′Z′繞Y′軸旋轉(zhuǎn)的角度是εY′,繞X′軸旋轉(zhuǎn)的角度是εX′，繞Z′軸旋轉(zhuǎn)角度為εZ′。

整個旋轉(zhuǎn)矩陣以各旋轉(zhuǎn)軸旋轉(zhuǎn)順序?qū)?yīng)的旋轉(zhuǎn)矩陣右乘，就可以得到旋轉(zhuǎn)矩陣R。

(2)

由式(1)、式(2)，求出ΔX、ΔY、ΔZ、εX′、εY′、εZ′、λ7個未知參數(shù)，就可以實現(xiàn)2個坐標系之間的相互轉(zhuǎn)換。

(3)

將公共點的重心化坐標代入式(1), 可得

(4)

(5)

因而，轉(zhuǎn)換參數(shù)可分兩步來求解, 即先用式(4)求出旋轉(zhuǎn)參數(shù)和比例因子, 再用式(5)求出平移參數(shù)。

1.1 傳統(tǒng)方法基于羅德里格矩陣的轉(zhuǎn)換

式(4)兩邊取2-范數(shù), 由于λ> 0 及旋轉(zhuǎn)矩陣為正交陣的特性，可得

(6)

對于n個公共點, 可得λ的最小均方估計

(7)

得到比例因子的最小均方估計后, 利用羅德里格矩陣的性質(zhì)[10], 可將旋轉(zhuǎn)矩陣R表示為

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

(8)

式中，I為單位矩陣，S為反對稱矩陣。

設(shè)

(9)

式中，a,b,c為羅德里格參數(shù)。將式(8)、式(9)代入式(4)得

(10)

展開得

(11)

因此，對于n個公共點, 可根據(jù)式(11)列出如下的總體誤差方程[11]

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

(12)

利用式(12)根據(jù)最小二乘原理, 無需迭代即可直接求得羅德里格參數(shù)

(13)

求得羅德里格參數(shù)后，可按式(8)求得旋轉(zhuǎn)矩陣R，然后再根據(jù)解算出的比例因子和旋轉(zhuǎn)參數(shù), 按式(5) 可求得平移參數(shù)。

圖3 最終局部坐標系和默認局部坐標系的相對空間位置

1.2 基于本文簡化方法的轉(zhuǎn)換

此時的最終局部坐標系與默認情況下坐標系的坐標的轉(zhuǎn)換關(guān)系為

(14)

即

(15)

顯然旋轉(zhuǎn)矩陣為

(16)

旋轉(zhuǎn)矩陣R即為默認局部坐標系下1、2、3軸的坐標構(gòu)成的矩陣。其中1軸為沿梁桿單元長度方向，后兩軸位于和用戶指定的單元方向相垂直的平面內(nèi)。由此可見，本文簡化方法避免了復(fù)雜的矩陣變換，直接可得出旋轉(zhuǎn)矩陣。

2 結(jié)果驗證

為了驗證實驗數(shù)據(jù)的正確性，現(xiàn)分3種情況驗證，如表1所示。

表1 驗證實驗數(shù)據(jù) m

設(shè)定梁桿單元不平與平行于3個軸，以1軸任意為例求得旋轉(zhuǎn)矩陣、尺度因子、平移參數(shù)：

取ang角分別為 -30°，30°，45°，60°，90°，120°，易知3 ′在默認局部坐標系下的坐標3 ′=[0,-sin(ang),cos(ang)]。分別以傳統(tǒng)方法基于羅德里格矩陣轉(zhuǎn)化方法和本文簡化方法求出3 ′在整體坐標系的坐標，如表2所示。

表2 傳統(tǒng)轉(zhuǎn)化方法和本文簡化轉(zhuǎn)化方法下的局部3 ′軸的轉(zhuǎn)化結(jié)果 m

從驗算結(jié)果可以看出，本簡化方法和傳統(tǒng)方法基于羅德里格矩陣轉(zhuǎn)化方法計算結(jié)果完全一致，表明本文方法是準確和有效的。

3 結(jié)論

針對梁桿單元的三維坐標轉(zhuǎn)換問題，列出了基于羅德里格矩陣的三維坐標轉(zhuǎn)換方法，計算過程繁瑣，不利于理解和應(yīng)用。在SAP2000、Etabs、和Midas等有限元軟件的梁桿單元的最終局部坐標系與整體坐標系的對應(yīng)過程中，提出了利用默認情況下的梁桿單元1、2、3軸的坐標組成轉(zhuǎn)換矩陣，從而簡便地實現(xiàn)最終局部坐標系在整體坐標系中的表示。由表2可以看出，基于本文簡化方法和基于羅德里格矩陣轉(zhuǎn)化方法計算的最終局部坐標系中局部3 ′軸在整體坐標系的坐標，兩種方法計算的結(jié)果完全一致，具有足夠的精度，表明本文方法是準確有效的。本文簡化方法理論嚴密，計算簡單，無小角度限制，簡化了旋轉(zhuǎn)矩陣的求解過程，提高了運算效率。