宋 浩， 曹聰慧2， 陳 偉， 郭 進(jìn)

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

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

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

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

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

圖1 相對(duì)于默認(rèn)方向的梁桿單元坐標(biāo)角

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

圖2 2個(gè)不同的空間直角坐標(biāo)系

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

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

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

(1)

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

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

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

(2)

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

(3)

將公共點(diǎn)的重心化坐標(biāo)代入式(1), 可得

(4)

(5)

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

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

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

(6)

對(duì)于n個(gè)公共點(diǎn), 可得λ的最小均方估計(jì)

(7)

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

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

(8)

式中，I為單位矩陣，S為反對(duì)稱矩陣。

設(shè)

(9)

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

(10)

展開(kāi)得

(11)

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

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

(12)

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

(13)

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

圖3 最終局部坐標(biāo)系和默認(rèn)局部坐標(biāo)系的相對(duì)空間位置

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

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

(14)

即

(15)

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

(16)

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

2 結(jié)果驗(yàn)證

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

表1 驗(yàn)證實(shí)驗(yàn)數(shù)據(jù) m

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

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

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

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

3 結(jié)論

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