王 凱

(東南大學(xué) 交通學(xué)院，江蘇 南京 210096)

0 引 言

在經(jīng)典的彈性力學(xué)專著或教材中，均未給出主應(yīng)力及其方向余弦的計算公式。在參考文獻[1]和[2]中，筆者在前人工作的基礎(chǔ)上，對3個主應(yīng)力及其方向余弦的計算公式進行了詳細的推導(dǎo)和闡述，而國外的一些著名的力學(xué)計算程序例如BISAR、ANSYS等，則利用計算方法中求解實對稱方陣的全部特征值和特征向量的雅可比(Jacobi)方法，由應(yīng)力矩陣計算出3個主應(yīng)力及其方向余弦。這兩種計算主應(yīng)力及其方向余弦的方法在理論上有何聯(lián)系，在各種力學(xué)狀況下其計算結(jié)果是否一致，這是廣大讀者一直關(guān)心的問題。筆者擬對上述問題做進一步的研究和探討。筆者曾閱讀多本國內(nèi)出版的彈性力學(xué)教材或?qū)Ｖ?，在這些教材或?qū)Ｖ校磳煞N方法的數(shù)學(xué)原理和計算過程作深入細致地介紹和論述，并通過計算對兩種方法的計算結(jié)果進行對比和討論。因此僅憑這些教材與專著中的內(nèi)容，無法很好地回答上述讀者所關(guān)心的問題。

1 特征值和特征向量的計算

1.1 特征值和特征向量的定義

由線性代數(shù)可知[3]，n階方陣

的特征值是使方程組

(1)

即

(2)

有非零解的λ值。這個非零解(x1,x2,…,xn)即為該方陣的屬于特征值λ的特征向量。

1.2 計算特征值和特征向量的第1種方法

根據(jù)齊次線性代數(shù)方程組有非零解的充要條件是系數(shù)行列式的值為零，特征值λ滿足n次多項式方程：

(3)

展開方程(3)，即可進一步得到方程：

b0λn+b1λn-1+…+bn=0

(4)

方程(4)稱為方陣A的特征方程，該方程的n個根(λ1,λ2,…,λn)即為方陣A的n個特征值。

由上述可得到確定方陣A的特征值和特征向量的第1種方法，該方法分成如下兩步：

1)求解方陣A的特征方程(4)的全部根(λ1,λ2,…,λn)，得到全部特征值。

2)把所求得的特征值逐個代入方程組(2)，對于每一個特征值求解出方程組(2)的解(x1,x2,…,xn)，即可得到相應(yīng)的特征向量。

由式(2)可知，如果(x1,x2,…,xn)是方程組(2)的一個非零解，那么(kx1,kx2,…,kxn)(k≠0)也是滿足方程組(2)的非零解。這說明特征向量不是被特征值所唯一決定的。在實際應(yīng)用問題中，往往對特征向量的取值加以“規(guī)格化”限制，即令所求算的特征向量值必須滿足式(5)：

(5)

任意去掉方程組(2)中的一個方程式，并將剩余的(n-1)個方程式與式(5)聯(lián)立求解，即可得到“規(guī)格化”的特征向量。該特征向量有兩組解答，如其中一組為(x1,x2,…,xn)，則另一組為(-x1,-x2,…,-xn)。

1.3 計算特征值和特征向量的第2種方法

由線性代數(shù)和計算方法可知[3-4]，如果P為n階非奇異方陣，P-1為其逆矩陣，則利用矩陣相乘相似變換(B=P-1AP)得到的n階方陣B與n階方陣A有相同的特征值。

實際中方陣A為實對稱的情況較為常見，關(guān)于其特征值和特征向量的計算比一般情況要簡單得多。線性代數(shù)理論表明，對于實對稱方陣A,存在一個n階正交方陣Q，使得：

Q-1AQ=D=diag(λ1,λ2,…,λn)

(6)

式中：λ1,λ2,…,λn為n階方陣A的特征值；D為由n個特征值構(gòu)成對角線元素的對角方陣。

n階正交方陣Q的k列qk正好為n階方陣A之相應(yīng)于第k個特征值λk(k=1,2,…,n)的特征向量。

在計算數(shù)學(xué)中，雅可比方法是計算實對稱方陣的全部特征值和特征向量的最常用方法。雅可比方法的基本思想是經(jīng)過一系列正交變換(雅可比旋轉(zhuǎn)變換)將n階實對稱方陣A一步一步化成由A的n個特征值構(gòu)成對角線元素的對角方陣D，即：

(R1R2…Rm)-1A(R1R2…Rm)=D=diag(λ1,λ2,…,λn)

(7)

式中：R1,R2,…,Rm為旋轉(zhuǎn)方陣。

與此同時，n階方陣Q=R1R2…Rm為所有旋轉(zhuǎn)方陣的乘積，其k列qk正好為n階方陣A之相應(yīng)于第k個特征值λk(k=1,2,…,n)的特征向量。

應(yīng)用雅可比方法計算實對稱方陣全部特征值和特征向量，已編成專用的計算機程序[5]。在該計算機程序中，考慮到實際應(yīng)用的需要，一般還對所求的特征向量值規(guī)格化。因此當雅可比計算機程序計算結(jié)束時，即可得到實對稱方陣全部的特征值和各特征向量平方和等于1的規(guī)格化特征向量。不過該特征向量解答不是全部解答，而僅僅是上述兩組規(guī)格化特征向量解答中的一組解答。

2 彈性力學(xué)中主應(yīng)力及其方向余弦的計算[6]

彈性力學(xué)中，求解主應(yīng)力σ及其方向余弦(l,m,n)的方程組為：

(8)

又有：

l2+m2+n2=1

(9)

將式(8)和式(9)聯(lián)立求解，即可得到σ,l,m,n的解答。齊次方程組(8)不能有l(wèi)=m=n=0這樣的解答，因為該解答與式(9)相矛盾。要使方程組(8)有別的解答，則必須取方程組的系數(shù)行列式等于零，即：

展開行列式得：

σ3-I1σ2+I2σ-I3=0

(10)

式中：I1、I2、I3為應(yīng)力狀態(tài)不變量，其值不隨坐標系的改變而改變，分別為：

三次方程式(10)稱為該應(yīng)力狀態(tài)的特征方程，它的3個實根σ1,σ2,σ3即為所求的主應(yīng)力。運用一元三次方程的理論求解三次方程式(10)，即可得到3個主應(yīng)力的計算式[1]：

(11)

與先前國內(nèi)各種文獻中出現(xiàn)的主應(yīng)力計算公式相比，筆者所列的主應(yīng)力計算公式無論是公式本身，還是在參考文獻[1]中公式的推導(dǎo)過程都更完善，因而也更為學(xué)術(shù)界所接受。

由參考文獻[6]，通過進一步力學(xué)分析，得到：

1)在彈性體內(nèi)任意一點，一定存在3個相互垂直的主應(yīng)力及其微分作用面(主平面)，主平面的外法線方向稱為應(yīng)力的主方向，主平面外法線的方向余弦與主應(yīng)力的方向余弦相同。每一個主應(yīng)力σi(i=1,2,3)及其主平面都是成對出現(xiàn)的。

2)在彈性體內(nèi)任意一點，上述6個微分主平面構(gòu)成一個微分六面體，同一個主應(yīng)力的兩個主平面相互平行，在其上分別作用著與其相對應(yīng)的兩個大小相等、方向相反的主應(yīng)力。

由上述可知，彈性體內(nèi)任意一點的每一個主應(yīng)力σi(i=1,2,3)的方向余弦都有兩組解答，如果其中一組是(li,mi,ni)，則另一組是(-li,-mi,-ni)。

在實際計算時求解主應(yīng)力及其方向余弦，也有兩種方法：

1)計算公式法

首先由前述的主應(yīng)力計算式(11)即可計算出3個主應(yīng)力，然后將所求出的3個主應(yīng)力σi(i=1,2,3)依次代入式(8)中任意兩式并與式(9)聯(lián)立求解，即可得到相應(yīng)的方向余弦li,mi,ni(i=1,2,3)。由式(8)和式(9)可知，與σi相對應(yīng)的方向余弦有兩組解答，如果其中一組為(li,mi,ni)，則另一組為(-li,-mi,-ni)，這與上面的力學(xué)分析一致。

由于參考文獻[2]已將主應(yīng)力方向余弦的計算公式推導(dǎo)出，因此在實際計算時，可直接利用該計算公式計算主應(yīng)力方向余弦。由參考文獻[2]可知，該文中主應(yīng)力方向余弦的求解只列出正數(shù)解的表達式，因此采用該文中公式僅計算了兩組主應(yīng)力方向余弦解答中的一組解答。不過一旦該解答li,mi,ni(i=1,2,3)得到后，另一組解答(-li,-mi,-ni)也很容易得到。

2)雅可比方法

如果將式(8)改寫成：

(12)

由式(12)即可看出主應(yīng)力σ是應(yīng)力矩陣S的特征值，應(yīng)力矩陣S如式(13):

(13)

主應(yīng)力方向余弦為應(yīng)力矩陣S的屬于特征值σ的特征向量。因此，可以利用計算方法中求解實對稱方陣全部特征值及其相應(yīng)規(guī)格化特征向量的雅可比方法及相應(yīng)的計算機程序，計算出全部主應(yīng)力和主應(yīng)力方向余弦。如前所述，采用該方法只能計算出兩組主應(yīng)力方向余弦解答中的一組解答。與計算公式法一樣，一旦該解答li,mi,ni(i=1,2,3)得到后，另一組解答(-li,-mi,-ni)也很容易得到。

3 計算算例

在筆者編寫的NESCP程序[7]中，計算主應(yīng)力及其方向余弦采用計算公式法，而在BISAR程序[8]中則采用雅可比方法。作者利用上述兩個程序計算了各種應(yīng)力狀況下的幾十個算例。由于篇幅有限，僅列出其中4個算例，如表1、表2、表3。每個算例中NESCP程序的計算結(jié)果精確到4位有效數(shù)字， BISAR程序的計算結(jié)果精確到3位有效數(shù)字。

表1 算例1～算例4應(yīng)力分量計算結(jié)果匯總Table 1 Summary of calculation results of stress components in example 1～example 4 MPa

表2 算例1～算例4主應(yīng)力及其方向余弦計算結(jié)果匯總Table 2 Summary of calculation results of principal stresses and their direction cosines in example 1～example 4

表1、表2的計算結(jié)果表明，利用兩個程序得到的主應(yīng)力的計算結(jié)果十分接近；主應(yīng)力方向余弦的計算結(jié)果則分兩種情況：一種情況是計算結(jié)果的絕對值十分接近，符號也相同，另一種情況是計算結(jié)果的絕對值十分接近，而符號相反。由前面所述可知，兩種計算方法實際上均只計算了兩組主應(yīng)力方向余弦解答中的一組解答[假設(shè)該解答為li,mi,ni(i=1,2,3)]，如果將另一組解答(-li,-mi,-ni)補齊，則由表3[表3中的計算數(shù)據(jù)在表2中算例4主應(yīng)力及其方向余弦計算數(shù)據(jù)的基礎(chǔ)上將另一組主應(yīng)力方向余弦解答(-li,-mi,-ni)補齊后得到] 可知，在此情況下，不僅主應(yīng)力的計算結(jié)果十分接近，每一個主應(yīng)力方向余弦的計算結(jié)果也十分接近，僅計算結(jié)果中主應(yīng)力方向余弦的排列次序可能不一樣。

表3 算例4主應(yīng)力及其方向余弦完整計算結(jié)果Table 3 Complete calculation results of principal stresses and their direction cosines in example 4

4 結(jié) 論

綜上所述，可得到如下結(jié)論：

1)彈性力學(xué)中計算主應(yīng)力及其方向余弦的兩種方法即計算公式法和雅可比方法與線性代數(shù)中計算特征值及其特征向量的兩種方法一脈相承，或者說彈性力學(xué)中計算主應(yīng)力及其方向余弦的兩種方法具有相同的數(shù)學(xué)淵源。

2)彈性體內(nèi)任意一點的每一個主應(yīng)力σi(i=1,2,3)及其主平面均成對出現(xiàn)。同一個主應(yīng)力的兩個主平面相互平行，在其上分別作用著與其相對應(yīng)的兩個大小相等、方向相反的主應(yīng)力。因此每一個主應(yīng)力σi(i=1,2,3)的方向余弦都有兩組解答，如果其中一組解答為(li,mi,ni)，則另一組解答為(-li,-mi,-ni)。目前采用計算公式法或雅可比方法計算主應(yīng)力方向余弦，實際上均只計算了兩組主應(yīng)力方向余弦解答中的一組解答，只有將另一組解答補齊，才能得到完整的計算結(jié)果。

3)計算表明，在各種力學(xué)狀況下上述兩種方法的計算結(jié)果十分接近，僅計算結(jié)果中主應(yīng)力方向余弦的排列次序可能不一樣。無論是采用計算公式法還是采用雅可比方法均能得到足夠精確的主應(yīng)力及其方向余弦的計算結(jié)果。

由上述結(jié)論可知，筆者的研究成果不僅很好地回答了讀者的問題，也是對經(jīng)典的彈性力學(xué)基礎(chǔ)理論做了一點很好的補充。

致謝：筆者對東南大學(xué)數(shù)學(xué)學(xué)院趙業(yè)鑫副教授的幫助謹表示衷心感謝。