基于MATLAB和矩陣位移法的平面桿系結(jié)構(gòu)有限元程序設(shè)計(jì)

闕仁波

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

0 引言

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

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

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

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

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

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

1 程序概況

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

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

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

(4)可自行控制輸出，比如輸出局部坐標(biāo)系中的單元?jiǎng)偠染仃?、坐?biāo)轉(zhuǎn)換矩陣、整體坐標(biāo)系中的單元?jiǎng)偠染仃嚨入A段性結(jié)果，以配合階段性演示或檢驗(yàn)。

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

2 編寫說(shuō)明

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

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

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

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

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

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

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

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

3.1 準(zhǔn)備工作

(1)劃分單元

②對(duì)結(jié)點(diǎn)位移編碼，以確定單元定位向量。

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

(4)確定單元的結(jié)點(diǎn)荷載。

3.2 輸入工作

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

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

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

表1 荷載類型

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

輸入JJHZ時(shí)要注意：

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

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

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

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

4 結(jié)果輸出和控制

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

①結(jié)點(diǎn)位移向量Δ即程序中的Delta；

圖1 局部坐標(biāo)系中的單元桿端力

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

5 算例

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

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

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

圖2 剛架算例

圖3 坐標(biāo)建立、單元?jiǎng)澐趾徒Y(jié)點(diǎn)位移編碼

(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 例二 用程序計(jì)算圖4所示連續(xù)梁的桿端內(nèi)力和桿端位移

圖4 連續(xù)梁算例

圖5 坐標(biāo)建立、單元?jiǎng)澐趾徒Y(jié)點(diǎn)位移編碼

(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 例三 用程序計(jì)算圖6所示桁架的桿端內(nèi)力

圖6 桁架算例

圖7 坐標(biāo)建立、單元?jiǎng)澐趾徒Y(jié)點(diǎn)位移編碼

(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 例四 用程序計(jì)算圖8所示組合結(jié)構(gòu)的桿端內(nèi)力

圖8 組合結(jié)構(gòu)算例

圖9 坐標(biāo)建立、單元?jiǎng)澐趾徒Y(jié)點(diǎn)位移編碼

(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 例五 用程序計(jì)算圖10所示剛架的桿端內(nèi)力

圖10 帶滑動(dòng)支座的剛架算例

圖11 坐標(biāo)建立、單元?jiǎng)澐趾徒Y(jié)點(diǎn)位移編碼

(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(yù)為一個(gè)大數(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 結(jié)語(yǔ)

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

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

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

7 附錄：流程圖和源代碼

圖12 計(jì)算流程圖

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

tic %計(jì)時(shí)開始

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

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

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

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

P=zeros(ZDM,1); %將整體結(jié)構(gòu)的等效結(jié)點(diǎn)荷載向量賦初值零

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

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

for N=1:DYSH

%計(jì)算局部坐標(biāo)系中的單元?jiǎng)偠染仃?/p>

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);%將所有局部坐標(biāo)系中的單元?jiǎng)偠染仃嚢葱羞B接生成一個(gè)大矩陣，以供后面求局部坐標(biāo)系中的單元桿端內(nèi)力向量時(shí)調(diào)用

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

%單元坐標(biāo)轉(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; %將所有局部坐標(biāo)系中的單元坐標(biāo)轉(zhuǎn)換矩陣按行連接生成一個(gè)大矩陣，以供后面求局部坐標(biāo)系中的單元桿端位移向量時(shí)調(diào)用

kezh=T'*kej*T; %整體坐標(biāo)系中的單元?jiǎng)偠染仃?/p>

%定位

FL=0;

for M=2:7

if BM(N,M)～=0

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

BMFL(FL,:)=[M-1,BM(N,M)]; %找出局部碼(M-1)與總碼(B(N,M))對(duì)應(yīng)關(guān)系

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

%求局部坐標(biāo)系中的單元固端約束力

SZJJ=size(JJHZ);

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

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; %單元固端約束力(局部坐標(biāo)系中)

end

end

end

Pe(:,N)=-T'*FeP(:,N); %整體坐標(biāo)系中的單元等效結(jié)點(diǎ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); %整體結(jié)構(gòu)的等效結(jié)點(diǎn)荷載向量

end

end

end

%將結(jié)點(diǎn)荷載在P中定位累加

SZJD=size(JDHZ);

if SZJD(1,1)～=0; %若作用結(jié)點(diǎn)荷載

for Z=1:SZJD(1,1)

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

end

end

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

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);

%由結(jié)點(diǎn)位移向量求整體坐標(biāo)系中的單元桿端位移

end

end

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

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

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

end

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

%求解局部坐標(biāo)系中單元桿端約束力的子函數(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 %計(jì)時(shí)結(jié)束

end