王 凱

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

0 引 言

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

1 特征值和特征向量的計(jì)算

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

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

的特征值是使方程組

(1)

即

(2)

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

1.2 計(jì)算特征值和特征向量的第1種方法

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

(3)

展開(kāi)方程(3)，即可進(jìn)一步得到方程：

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

(4)

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

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

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

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

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

(5)

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

1.3 計(jì)算特征值和特征向量的第2種方法

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

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

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

(6)

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

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

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

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

(7)

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

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

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

2 彈性力學(xué)中主應(yīng)力及其方向余弦的計(jì)算[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這樣的解答，因?yàn)樵摻獯鹋c式(9)相矛盾。要使方程組(8)有別的解答，則必須取方程組的系數(shù)行列式等于零，即：

展開(kāi)行列式得：

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

(10)

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

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

(11)

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

由參考文獻(xiàn)[6]，通過(guò)進(jìn)一步力學(xué)分析，得到：

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

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

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

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

1)計(jì)算公式法

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

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

2)雅可比方法

如果將式(8)改寫(xiě)成：

(12)

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

(13)

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

3 計(jì)算算例

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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