闕仁波

(廈門大學嘉庚學院土木系 福建漳州 363105)

0 引言

在結構力學發(fā)展史上，為了解超靜定結構，19世紀建立了力法。但用于解之后大量出現(xiàn)的高次超靜定剛架，仍顯繁瑣。在20世紀初，以力法為基礎，發(fā)展求解高次超靜定結構比力法更具優(yōu)勢的位移法[1]，但隨著基本未知量的增多，位移法的手算法依然讓人不堪重負。隨著計算機的發(fā)展，位移法被以矩陣的方式來表達，即矩陣位移法，以適合編程計算，從而在傳統(tǒng)手算法的基礎上，發(fā)展了程序結構力學，極大地增強了求解高次超靜定結構的能力，并將設計者很大程度地從繁瑣的計算中解放出來，將更多的精力用于概念設計和定性分析[1-2]。

在適合桿件結構的矩陣位移法的啟示下，又發(fā)展了適合連續(xù)體的矩陣方法——有限元法，故把矩陣位移法稱為桿件有限元法[3-6]，但相對于連續(xù)體有限元法，門檻低，故若由此入門，將其核心思想、計算流程和程序編制掌握，然后，再拾階而上，在求同中類比而廣義化，在變異中對比而延拓，從一維桿件、線性(線性代數(shù)之線性、本構關系之線性)到多維連續(xù)體、非線性，不斷朝橫向和縱深發(fā)展，不僅難度梯度比直接由彈性力學有限元入門小得多，亦容易建構起合理而有序的知識體系。

此外，桿件有限元在土木工程的專業(yè)課中應用較多，比如從總體受力來看，橋梁的特點是長而不寬，它的受力特性與桿系結構相符，一般采用平面桿系有限元法來分析[7]。但對土木專業(yè)的本科生而言，有限元大多只是選修課，而在專業(yè)軟件如橋梁博士和PKPM中，卻要用到有限元的概念，故在必修的矩陣位移法教學中體現(xiàn)一點有限元的思想對后續(xù)課程學習有較大幫助。

而現(xiàn)在的矩陣位移法教學卻存在一些問題。一方面，可能受學時所限或?qū)W生沒修過編程課，教學只局限于原理的講解，而線性代數(shù)與結構力學的無“機”結合，使得缺乏程序?qū)崿F(xiàn)而采用手算的矩陣位移法，不僅沒顯示出其相對于傳統(tǒng)位移法和漸進法的優(yōu)越性，反而讓學生更加望而卻步，該部分內(nèi)容的內(nèi)在規(guī)律呼喚理論教學和程序?qū)崿F(xiàn)的有“機”結合。另一方面，某些觀念認為，都那么多現(xiàn)成的軟件，沒必要再自主開發(fā)。可在輸入和輸出之間，作為內(nèi)核的有限元對很多程序使用者而言卻是“黑盒子”[4]，盡管程序設計的智能化造就了軟件應用的機械化，但若不加理解，面對復雜問題時出錯的可能性會增大。而對于一個從事過開發(fā)的人而言，在應用別人開發(fā)的軟件時，會有一種同理心的優(yōu)勢，站在開發(fā)者的角度來審視軟件，透視其內(nèi)核而不覺其內(nèi)黑，從而能更快更好地掌握軟件的應用，并在欣賞中繼承，在批判中創(chuàng)新。在“計算機+”和人工智能的時代，這種同理心有助于更好地去理解和適應“機”智的內(nèi)在原理。此外，對一些不在“通用”之列的個性化問題，亦要靠自主開發(fā)程序來解決。故編程訓練仍有其存在的必要性。且矩陣位移法所蘊含的從一般到特殊的演繹(如從剛架到連續(xù)梁、桁架和組合結構)、化整為零的單元分析、集零為整的整體綜合思想，更具有方法論上的意義。通過編程，可提高一個人對復雜系統(tǒng)內(nèi)在邏輯認識的清晰度，將復雜系統(tǒng)分解為簡單子系統(tǒng)逐個加以處理，然后再加以集成，可將“知識單元”富有邏輯地集成為“知識整體”，從“構件”到“結構”，建構起牢固而“幾何不變的知識體系”。

基于矩陣位移法的程序開發(fā)，可采用不同的語言，如Fortran[2,8]、C++[5]和VB[9]等。而采用以擅長矩陣計算著稱的MATLAB(即矩陣實驗室)來處理“矩陣位移法”之“矩陣”，將會是一種非常具有優(yōu)勢的選擇，且該語言簡單易學，應用更廣泛，故目前采用MATLAB編寫的程序開始增多[4，6，10-13]，但這些資料中，要么從彈性力學和一般有限元開始講解[4,6,10]，對于本科生來說，難度相對較大；要么將連續(xù)梁、桁架和剛架分門別類地講[11，12]，篇幅較大。如何與結構力學教材直接配套，在一邊理論教學的同時，一邊快速切入程序?qū)崿F(xiàn)，而無需過多參考資料，以較小的跨度在傳統(tǒng)手算和程序計算之間架起一座溝通的橋梁，對大多數(shù)受學時所限的本科教學來說，是非常必要的。此外，以較短的時間讓學生一窺從理論到程序?qū)崿F(xiàn)的全過程，可更好地激起他們的興趣，然后再去拓展，教學效果會更好。而通過閱讀程序、修改程序到自主編程逐步進階，不失為一種入門的方法。

鑒于上述目的，本文將基于目前高校廣泛采用的結構力學教材—文獻[14](盡管它配套有求解器，為文獻[2]的程序?qū)崿F(xiàn)，但文獻[2]是以Fortran語言編寫)，采用MATLAB語言編寫一個程序，以供教學參考。

1 程序概況

(1)流程圖和源代碼請見附錄。

(2)該程序充分發(fā)揮了MATLAB擅長矩陣計算的優(yōu)勢，代碼十分簡潔，除注釋行外，才115行。

(3)該程序基于文獻[14]編寫，故對采用該教材的師生，參照非常方便。符號系統(tǒng)同文獻[14]，故限于篇幅，本文不再贅述。

(4)可自行控制輸出，比如輸出局部坐標系中的單元剛度矩陣、坐標轉(zhuǎn)換矩陣、整體坐標系中的單元剛度矩陣等階段性結果，以配合階段性演示或檢驗。

(5)為便于沒學過有限元英文專業(yè)詞匯的學生記憶，程序中的變量名大多采用漢化的表述模式，如JJHZ-結間荷載，JDHZ-結點荷載。

2 編寫說明

(1)該程序按分析剛架的一般模式編寫。

(2)對于邊界條件，均采用先處理法，即令結點位移為零的總碼為零，如圖3所示。

(3)將連續(xù)梁和桁架看作特殊的剛架，具體處理方式如下：

①對于連續(xù)梁(忽略軸向變形)，每個單元亦按6個桿端位移輸入，令4個線位移總碼為零即可；

②對于桁架，每個單元亦按6個桿端位移輸入，但同一鉸結點處各桿的結點線位移相同，故采用重碼；但結點角位移不同，故采用異碼。如圖7中單元①、②和⑤交點處的編碼。

③上述處理模式對于單純的連續(xù)梁和單純的桁架結構，可能比按單純計算連續(xù)梁和桁架的程序顯得麻煩。但本文中，一方面更注重程序的通用性，希望以一般而非特殊的思維來統(tǒng)一看待不同類型的結構；另一方面，本文的方法在處理諸如例四所示的組合結構時，就顯示出其優(yōu)越性。

(4)組合結點處，編碼時要注意線位移的重碼和轉(zhuǎn)角的異碼問題，如例四中A點和B點處。

3 數(shù)據(jù)準備和輸入

3.1 準備工作

(1)劃分單元

②對結點位移編碼，以確定單元定位向量。

(3)確定單元的結間荷載。

(4)確定單元的結點荷載。

3.2 輸入工作

(1)BM——編碼矩陣：每一行第1列元素代表單元號，第2～7列元素代表該單元定位向量(以行向量形式輸入)的6個元素。

(2)CLJ——材料參數(shù)和夾角矩陣：每一行第1列元素代表單元號，第2～6列元素分別代表該單元的彈性模量E、截面面積A、截面慣性矩I、長度l和夾角α。

(3)JJHZ——結間荷載矩陣：每一行第1列元素代表單元號，第2～5列元素分別代表該單元的荷載類型號、荷載大小(q、Fp、M)、a和b，如表1所示。

表1 荷載類型

該程序只考慮了表1所示的4種情況，更多可參閱文獻[14]表11-2，讀者可在求單元固端約束力的子程序function[GDL]=g(LX,F,a,b)中增加case即可擴展荷載類型、考慮支座移動或溫變的影響。

輸入JJHZ時要注意：

①有荷載的單元才需輸入；

②當同一單元有多種荷載，每一種均需輸入，如例一的單元1；

③若所有單元均無結間荷載(如桁架的情形)，則JJHZ為空矩陣，如例三。

(4)JDHZ—結點荷載矩陣：每一行第1列元素代表結點荷載所對應的結點位移編碼，第2列元素代表結點荷載大小。如例三中JDHZ=[1 10;2 -10]。若所有結點均沒荷載作用，則JDHZ為空矩陣，如例一、例四和例五。

4 結果輸出和控制

(1)一般控制輸出的結果為：

①結點位移向量Δ即程序中的Delta；

圖1 局部坐標系中的單元桿端力

(2)用戶可自行控制參數(shù)的輸出，只需在function[ ]=f(BM,CLJ,JJHZ,JDHZ) 左側(cè)方括號中填入或修改需要輸出的參數(shù)即可。如例二，[DeltaT,GDNLT]=f(BM,CLJ,JJHZ,JDHZ)。習慣于行向量輸出還是列向量輸出，只需通過在程序中轉(zhuǎn)置即可實現(xiàn)。亦可輸出階段性結果，以配合階段性演示或檢驗的需要。

5 算例

注：荷載作用下，超靜定結構的內(nèi)力分布只與各桿的相對剛度有關，與絕對剛度無關，而位移與絕對剛度有關。故在以下算例中，當各桿的E、A、I值以相對大小輸入時，算出的內(nèi)力屬于真實值，而位移真實值則需輸入實際的絕對值進行計算。

為體現(xiàn)程序的通用性，算例的選擇上著重體現(xiàn)結構類型、支座類型、截面變化情況和荷載類型的多樣化。

5.1 例一 用程序計算圖2所示剛架的桿端內(nèi)力

圖2 剛架算例

圖3 坐標建立、單元劃分和結點位移編碼

(1)輸入

>>BM=[1 1 2 3 4 5 6;2 1 2 3 0 0 0;3 4 5 7 0 0 0];

>>CLJ=[1 1 0.6 1/6 5 0;2 0.9 0.5 1/12 5pi/2;3 0.9 0.5 1/12 5pi/2];

>>JJHZ=[1 1 4.8 5 0;1 2 6 4 1;3 2 -8 2 3];

>>JDHZ=[];

(2)執(zhí)行

>>[GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

(3)輸出

Elapsed time is 0.133784 seconds.

GDNLT=

5.2 例二 用程序計算圖4所示連續(xù)梁的桿端內(nèi)力和桿端位移

圖4 連續(xù)梁算例

圖5 坐標建立、單元劃分和結點位移編碼

(1)輸入

>>BM=[1 0 0 0 0 0 1;2 0 0 1 0 2 3;3 0 2 3 0 0 4;4 0 0 4 0 5 6];

>>CLJ=[1 1 2 3 3 0;2 1 2 3 2 0;3 1 1.5 2 3 0;4 1 1.5 2 2 0];

>>JJHZ=[4 1 10 2 0];

>>JDHZ=[1 50];

(2)執(zhí)行

[DeltaT,GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

(3)輸出

Elapsed time is 0.024951 seconds.

DeltaT=

6.76512.0539-2.80277.874425.748814.5411

GDNLT=

5.3 例三 用程序計算圖6所示桁架的桿端內(nèi)力

圖6 桁架算例

圖7 坐標建立、單元劃分和結點位移編碼

(1)輸入

>>BM=[1 1 2 5 0 0 6;2 1 2 7 3 4 8;3 3 4 9 0 0 10;4 0 0 11 0 0 12;5 1 2 13 0 0 14 ;6 3 4 15 0 0 16 ];

>>CLJ=[1 1 1 1 1pi/2;2 1 1 1 1 0;3 1 1 1 1pi/2;4 1 1 1 1 0;5 1 1 1 2^0.5pi/4;6 1 1 1 2^0.5 3*pi/4];

>>JJHZ=[];

>>JDHZ=[1 10;2 -10];

(2)執(zhí)行

>>[GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

(3)輸出

Elapsed time is 0.026574 seconds.

GDNLT=

5.4 例四 用程序計算圖8所示組合結構的桿端內(nèi)力

圖8 組合結構算例

圖9 坐標建立、單元劃分和結點位移編碼

(1)輸入

>>BM=[1 0 0 0 1 2 3;2 1 2 3 4 5 6;3 4 5 6 0 0 0;4 0 0 7 1 2 8;5 4 5 9 0 0 10];

>>CLJ=[1 1 2 1 20 0;2 1 2 1 20 0;3 1 2 1 20 0;4 1 1/20 1 25 atan(3/4);5 1 1/20 1 25 -atan(3/4)];

>>JJHZ=[2 1 10 20 0];

>>JDHZ=[];

(2)執(zhí)行

>>[GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

(3)輸出

Elapsed time is 0.025520 seconds.

GDNLT=

5.5 例五 用程序計算圖10所示剛架的桿端內(nèi)力

圖10 帶滑動支座的剛架算例

圖11 坐標建立、單元劃分和結點位移編碼

(1)輸入

>>BM=[1 1 0 2 3 4 5;2 3 4 5 0 6 0;3 3 4 5 0 0 0];

>>CLJ=[1 1 1 1 4 0;2 1 1 2 2 0;3 1 1 1 4 pi/2];

>>JJHZ=[1 2 10 2 2];

>>JDHZ=[ ];

(2)執(zhí)行

>>[GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

(3)輸出

Elapsed time is 0.003276 seconds.

GDNLT =

若忽略軸向變形，則可在CLJ中令A為一個大數(shù)，比如A=106，則CLJ=[1 1 10^6 1 4 0;2 1 10^6 2 2 0;3 1 10^6 1 4 pi/2]，重新執(zhí)行可得：

GDNLT=

6 結語

(1)基于MATLAB語言和矩陣位移法思想，進行了平面桿系結構有限元程序設計。該程序注重通用性，對剛架、連續(xù)梁、桁架及組合結構，都采用統(tǒng)一的分析模式，這對分析組合結構尤其具有優(yōu)勢。

(2)該程序充分發(fā)揮了MATLAB——矩陣實驗室擅長矩陣計算的優(yōu)勢，代碼簡潔，數(shù)據(jù)準備簡單、計算流程清晰、結果輸出可選控，可作為矩陣位移法原理學習之外的一種實踐補充，與教材緊密配套，切入方便快捷，以較低的門檻，引領學生初窺有限元。

(3)該程序旨在拋磚引玉，讀者可在此基礎上進一步修改拓展，以處理更多實際情況，比如增加荷載類型、考慮支座位移和溫變的影響、考慮斜支座和彈性支座、考慮剪切變形和帶剛域桿[15]、基于MATLAB GUI開發(fā)可視化圖形用戶界面和內(nèi)力圖自動繪制功能等，進一步增加其通用性和智能化程度。

7 附錄：流程圖和源代碼

圖12 計算流程圖

function [GDNLT]=f(BM,CLJ,JJHZ,JDHZ)

tic %計時開始

DYSH=max(BM(:,1)); %計算單元數(shù)

ZDM=max(max(BM(:,2:7))); %計算最大總碼

K=zeros(ZDM,ZDM); %將整體剛度矩陣賦初值零

FeP=zeros(6,DYSH); %將由局部坐標系中的固端約束力向量按列連接而成的矩陣賦初值零

P=zeros(ZDM,1); %將整體結構的等效結點荷載向量賦初值零

GDWY=zeros(6,DYSH);%將由局部坐標系中的單元桿端位移向量按列連接而成的矩陣賦初值零

GDNL=zeros(6,DYSH);%將由局部坐標系中的單元桿端內(nèi)力向量按列連接而成的矩陣賦初值零

for N=1:DYSH

%計算局部坐標系中的單元剛度矩陣

kej=zeros(6,6);

EAL=CLJ(N,2)*CLJ(N,3)/CLJ(N,5);

kej([1 4 19 22])=[EAL -EAL -EAL EAL];

EIL3=12*CLJ(N,2)*CLJ(N,4)/CLJ(N,5)^3;

kej([8 11 26 29])=[EIL3 -EIL3 -EIL3 EIL3];

EIL2=6*CLJ(N,2)*CLJ(N,4)/CLJ(N,5)^2;

kej([9 12 14 32])=EIL2;

kej([17 27 30 35])=-EIL2;

kej([15 36])=4*CLJ(N,2)*CLJ(N,4)/CLJ(N,5);

kej([18 33])=2*CLJ(N,2)*CLJ(N,4)/CLJ(N,5);%將所有局部坐標系中的單元剛度矩陣按行連接生成一個大矩陣，以供后面求局部坐標系中的單元桿端內(nèi)力向量時調(diào)用

kejall((N*6-5):N*6,:)=kej;

%單元坐標轉(zhuǎn)換矩陣

T=zeros(6,6);

T([1 8 22 29])=cos(CLJ(N,6));

T([7 28])=sin(CLJ(N,6));

T([2 23])=-sin(CLJ(N,6));

T([15 36])=1;

Tall((N*6-5):N*6,:)=T; %將所有局部坐標系中的單元坐標轉(zhuǎn)換矩陣按行連接生成一個大矩陣，以供后面求局部坐標系中的單元桿端位移向量時調(diào)用

kezh=T'*kej*T; %整體坐標系中的單元剛度矩陣

%定位

FL=0;

for M=2:7

if BM(N,M)～=0

FL=FL+1; %統(tǒng)計單元中非零總碼的個數(shù)

BMFL(FL,:)=[M-1,BM(N,M)]; %找出局部碼(M-1)與總碼(B(N,M))對應關系

end

end

%累加

for I=1:FL

for J=1:FL

K(BMFL(I,2),BMFL(J,2))=K(BMFL(I,2),BMFL(J,2))...

+kezh(BMFL(I,1),BMFL(J,1));

end

end

%求局部坐標系中的單元固端約束力

SZJJ=size(JJHZ);

if SZJJ(1,1)～=0 %若結間作用荷載

for S=1:SZJJ(1,1)

if JJHZ(S,1)==N

[GDL]=g(JJHZ(S,2),JJHZ(S,3),JJHZ(S,4),JJHZ(S,5));

%調(diào)子函數(shù)

FeP(:,N)=FeP(:,N)+GDL; %單元固端約束力(局部坐標系中)

end

end

end

Pe(:,N)=-T'*FeP(:,N); %整體坐標系中的單元等效結點荷載向量

%將Pe在P中定位累加

for S=2:7

if BM(N,S)～=0

P(BM(N,S),1)= P(BM(N,S),1)+Pe(S-1,N); %整體結構的等效結點荷載向量

end

end

end

%將結點荷載在P中定位累加

SZJD=size(JDHZ);

if SZJD(1,1)～=0; %若作用結點荷載

for Z=1:SZJD(1,1)

P(JDHZ(Z,1),1)= P(JDHZ(Z,1),1)+JDHZ(Z,2);

end

end

Delta=KP; %(編碼非零)結點位移向量，按編碼從小到大的順序排列

DeltaT=Delta'; %轉(zhuǎn)置成行向量表示

for II=1:DYSH

for JJ=2:7

if BM(II,JJ)～=0

GDWY(JJ-1,II)=GDWY(JJ-1,II)+Delta(BM(II,JJ),1);

%由結點位移向量求整體坐標系中的單元桿端位移

end

end

DYWY=Tall((II*6-5):II*6,:)*GDWY(:,II); %局部坐標系中的單元桿端位移

GDNL(:,II)=GDNL(:,II)+...

kejall((II*6-5):II*6,:)*DYWY+ FeP(:,II); %局部坐標系中的單元桿端內(nèi)力

end

GDNLT=GDNL'; %轉(zhuǎn)置，以便輸出時以每單元一行而非每單元一列，節(jié)省文中顯示空間

%求解局部坐標系中單元桿端約束力的子函數(shù)

function [GDL]=g(LX,F,a,b)

L=a+b;

switch LX %判斷荷載類型

case 1

FxP1=0;

FxP2=0;

FyP1=-F*a*(1-a^2/L^2+a^3/(2*L^3));

FyP2=-F*a^3/L^2*(1-a/(2*L));

MP1=-F*a^2/12*(6-8*a/L+3*a^2/L^2);

MP2=F*a^3/(12*L)*(4-3*a/L);

GDL=[FxP1 FyP1 MP1 FxP2 FyP2 MP2]';

case 2

FxP1=0;

FxP2=0;

FyP1=-F*b^2/L^2*(1+2*a/L);

FyP2=-F*a^2/L^2*(1+2*b/L);

MP1=-F*a*b^2/L^2;

MP2=F*a^2*b/L^2;

GDL=[FxP1 FyP1 MP1 FxP2 FyP2 MP2]';

case 3

FxP1=0;

FxP2=0;

FyP1=6*F*a*b/L^3;

FyP2=-6*F*a*b/L^3;

MP1=F*b/L*(2-3*b/L);

MP2=F*a/L*(2-3*a/L);

GDL=[FxP1 FyP1 MP1 FxP2 FyP2 MP2]';

Otherwise

FxP1=0;

FxP2=0;

FyP1=-F*a/4*(2-3*a^2/L^2+1.6*a^3/L^3);

FyP2=-F/4*a^3/L^2*(3-1.6*a/L);

MP1=-F*a^2/6*(2-3*a/L+1.2*a^2/L^2);

MP2=F*a^3/(4*L)*(1-0.8*a/L);

GDL=[FxP1 FyP1 MP1 FxP2 FyP2 MP2]';

end

end

toc %計時結束

end