[Matlab科學計算] 1.你想要了解的尤拉角
[問題由來]
在計算鐵磁材料多晶體的有效模量時,需要考慮晶粒在多晶體中的方向分佈,一般用三個尤拉角(, , )來表示晶粒在多晶體中的方向,用方向分佈函式來表示某個方向的分佈密度。基於此,迫使我要掌握尤拉角,但是在閱讀眾多教材和部落格文章中發現,大家對尤拉角的說法不是很統一,所以,基於我的理解整一下尤拉角的相關概念及使用注意事項。
一、幾個概念
1.1經典尤拉角(Proper Euler angles)和泰勒布萊恩角(Tait-Bryan angles)
這兩種尤拉角是按照旋轉軸的順序定義的,經典尤拉角是按照Z-X-Z,Y-X-Y,X-Y-X,Z-Y-Z,X-Z-X,Y-Z-Y這樣的軸序旋轉,即第一個旋轉軸和最後一個旋轉軸相同;而泰勒布萊恩角是按照X-Y-Z,Y-Z-X,Z-X-Y,X-Z-Y,Z-Y-X,Y-X-Z這樣的軸序旋轉,即每次旋轉軸都不相同。
1.2內在旋轉(intrinsic rotations)和外在旋轉(extrinsic rotations)
這兩個概念通過下圖進行說明:即內在旋轉每次旋轉圍繞的軸是上次旋轉之後座標系的某個軸,外在旋轉每次旋轉的軸是固定座標系中的軸。
內在旋轉與外在旋轉的轉換關係:互換第一次和第三次旋轉的位置則兩者結果相同。例如Z-Y-X旋轉和內部旋轉和X-Y-Z旋轉的外部旋轉的旋轉矩陣相同。
1.3旋轉矩陣(很重要)
主動旋轉和被動旋轉:主動旋轉是指將向量或座標系逆時針圍繞旋轉軸旋轉,被動選是對座標軸進行的逆時針旋轉,相當於主動旋轉的逆操作。
旋轉矩陣有二維旋轉和三維旋轉,研究需要這裡僅討論三維旋轉(被動旋轉)。
假設OP為單位向量且在XY平面內,則OP向量在O-XYZ座標系中座標為,在新座標系O-X’Y’Z’中的座標為,展開之後為,寫成矩陣的形式如下:
將上式中的3*3的變換矩陣稱為旋轉矩陣。從Fig. 3中可以看出,被動旋轉相當於將向量OP繞Z軸主動順時針旋轉θ角。
同理,可以得到繞X軸和繞Y軸被動旋轉
,
相應的向量繞座標軸主動旋轉的旋轉矩陣如下所示:
, ,
二、舉例計算
計算旋轉矩陣一定要明確旋轉順序和旋轉模式(內在旋轉還是外在旋轉,主動旋轉還是被動旋轉),由於內在旋轉易於圖解表達,通常使用內在旋轉說明。如圖4所示,旋轉順序為Z-X-Z的被動內在旋轉。假設向量OP在O-xyz座標系中的座標為OP0=(x, y, z),在O-XYZ座標系中的座標為OP=(X, Y, Z),分解計算步驟如下:
- 繞Z軸旋轉之後的座標為;
- 繞第一次旋轉之後的X座標旋轉之後的座標為;
- 繞第二次旋轉之後的Z座標旋轉之後的座標為;
可以得到OP向量在新座標系中的座標與在舊座標系中的座標之間的關係為 ,所以總的旋轉矩陣為,採用Matlab的符號計算可以得到該旋轉矩陣的展開形式。為了便於Matlab中書寫,令,,。
Matlab程式碼:
syms a b c
Rza = [cos(a) sin(a) 0;-sin(a) cos(a) 0;0 0 1];
Rxb = [1 0 0;0 cos(b) sin(b);0 -sin(b) cos(b)];
Rzc = [cos(c) sin(c) 0;-sin(c) cos(c) 0;0 0 1];
T = Rzc*Rxb*Rza;
disp(T);
計算結果:
[ cos(a)*cos(c) - cos(b)*sin(a)*sin(c), cos(c)*sin(a) + cos(a)*cos(b)*sin(c), sin(b)*sin(c)]
[-cos(a)*sin(c) - cos(b)*cos(c)*sin(a), cos(a)*cos(b)*cos(c) - sin(a)*sin(c), cos(c)*sin(b)]
[ sin(a)*sin(b), -cos(a)*sin(b), cos(b)]
三、在晶體織構學中的應用
根據尤拉角,可以得到旋轉矩陣為:
根據旋轉變換,則有如下公式:
進而,得到晶粒座標系C中的向量在試樣座標系S中的表達:
假設所有晶粒在試樣中的方向分佈函式(ODF)為,而且有 ,多晶體中常計算等效彈性剛度張量(四階張量),所以得出,一個在試樣座標系S中的四階張量 能夠從晶粒座標系C中的計算得到,
四階張量的體積平均變成取向平均,計算方法如下:
四、由旋轉前後的向量值計算旋轉矩陣(擴充套件)
- 兩個向量點乘計算夾角;
- 兩個向量叉乘計算旋轉軸;
- 利用羅德里格旋轉公式(Rodrigues' rotation formula)計算旋轉矩陣;
這不詳細內容和計算步驟以及程式見參考資料4.
致謝
感謝所有參考的部落格作者!