matlab連續(xù)梁程序的編制與使用

1、第三章 連續(xù)梁程序的編制與使用結(jié)構(gòu)力學(xué)中的矩陣位移法是隨著電子計算機進入結(jié)構(gòu)力學(xué)領(lǐng)域中而產(chǎn)生的一種方法，而Matlab語言正是進行矩陣運算的強大工具，因此，用Matlab語言編寫結(jié)構(gòu)力學(xué)程序有更大的優(yōu)越性。本章將詳細(xì)介紹如何利用Matlab語言編制連續(xù)梁結(jié)構(gòu)的計算程序。矩陣位移法的解題思路是將結(jié)構(gòu)離散為單元（桿件），建立單元桿端力與桿端位移之間的關(guān)系單元剛度方程；再將各單元集成為原結(jié)構(gòu)，在滿足變形連續(xù)條件和平衡條件時，建立整體剛度方程；在邊界條件處理完畢后，由整體剛度方程解出節(jié)點位移，進而求出結(jié)構(gòu)內(nèi)力。用矩陣位移法計算連續(xù)梁的步驟如下：1） 整理原始數(shù)據(jù)，如材料性質(zhì)、荷載條件、約束條件等，并進

2、行編碼：單元編碼、結(jié)點編碼、結(jié)點位移編碼、選取坐標(biāo)系。2） 建立局部坐標(biāo)系下的單元剛度矩陣。3） 建立整體坐標(biāo)系下的單元剛度矩陣。4） 集成總剛。5） 建立整體結(jié)構(gòu)的等效節(jié)點荷載和總荷載矩陣6） 邊界條件處理。7） 解方程，求出節(jié)點位移。8） 求出各單元的桿端內(nèi)力。圖3-1 程序流程圖實際上，上述步驟也是編制Matlab程序的基本步驟，在求出計算結(jié)果后，還可以利用Matlab的繪圖功能繪制結(jié)構(gòu)圖、內(nèi)力圖、變形圖等等。3.1 程序說明%*% 矩陣位移法解連續(xù)梁主程序%* l 功能：運用矩陣位移法解連續(xù)梁的基本原理編制的計算主程序。l 基本思想：結(jié)點（結(jié)點位移）編碼默認(rèn)為從左至右，從1開始順序進行

3、；桿端彎矩的方向默認(rèn)為逆時針。l 荷載類型：可計算結(jié)點荷載，每單元作用的跨中集中力和均布荷載。l 說明：主程序的作用是通過賦值語句、讀取和寫入文件、函數(shù)調(diào)用等完成算法的全過程，即實現(xiàn)程序流程圖的程序表達(dá)。%-1 程序準(zhǔn)備format short e %設(shè)定輸出類型clear all%清除所有已定義變量clc%清屏l 說明：format short e 設(shè)定計算過程中顯示在屏幕上的數(shù)字類型為短格式、科學(xué)計數(shù)法； clear all 清除所有已定義變量，目的是在本程序的運行過程中，不會發(fā)生變量名相同等可能使計算出錯的情況； clc 清屏，使屏幕在本程序運行開始時%-2 打開文件FP1=fopen(

4、input.txt,rt); %打開輸入數(shù)據(jù)文件 存放初始數(shù)據(jù)FP2=fopen(output.txt,wt); %打開輸出數(shù)據(jù)文件 存放計算結(jié)果l 說明：FP1=fopen(input.txt,rt); 打開已存在的輸入數(shù)據(jù)文件input.txt，且設(shè)置其為只讀格式，使程序在執(zhí)行過程中不能改變輸入文件中的數(shù)值，并用文件句柄FP1來FP2=fopen(output.txt,wt);打開輸出數(shù)據(jù)文件，該文件不存在時，通過此命令創(chuàng)建新文件，該文件存在時則將原有內(nèi)容全部刪除。該文件設(shè)置為可寫格式，可在程序執(zhí)行過程中向輸出文件寫入數(shù)據(jù)。%-3 讀入程序控制信息NELEM=fscanf(FP1,%d,1

5、); %單元個數(shù)（單元編碼總數(shù)） NPOIN=fscanf(FP1,%d,1); %結(jié)點個數(shù)（結(jié)點編碼總數(shù)） NVFIX=fscanf(FP1,%d,1);約束個數(shù)（零位移總數(shù)） NFPOIN=fscanf(FP1,%d,1);結(jié)點荷載個數(shù)（作用在結(jié)點上集中力偶總數(shù)） NFPRES=fscanf(FP1,%d,1);非結(jié)點荷載數(shù)（作用在單元上分布荷載總數(shù)） YOUNG=fscanf(FP1,%f,1);彈性模量l 說明：從輸入文件FP1中讀入單元個數(shù)，結(jié)點個數(shù)，約束個數(shù)，結(jié)點荷載個數(shù)，非結(jié)點荷載個數(shù)，彈性模量；程序中彈性模量僅輸入了一個值，表明本程序僅能求解一種材料構(gòu)成的結(jié)構(gòu)，如：鋼筋混凝土結(jié)

6、構(gòu)、鋼結(jié)構(gòu)，不能求解鋼筋混凝土鋼組合結(jié)構(gòu)。采用了命令fscanf，其中：%d表示讀入整數(shù)格式；1表示讀取1個數(shù)。%-fprintf(FP2,n 結(jié)構(gòu)初始數(shù)據(jù)nn);fprintf(FP2, 單元總數(shù)=%d 結(jié)點總數(shù)=%d 約束總數(shù)=%d n,NELEM,NPOIN,NVFIX); fprintf(FP2, 結(jié)點荷載個數(shù)=%d 非結(jié)點荷載個數(shù)=%d 彈性模量=%12.5g n,NFPOIN,NFPRES,YOUNG); l 說明：在輸出文件FP2中顯示輸入的控制信息，便于程序執(zhí)行完畢后，從輸出文件中查找輸入錯誤。采用了命令fprintf，其中：n為換行標(biāo)志，表示換一行；%12.5g表示輸出12位

7、且有5位小數(shù)的實數(shù)。括在引號中的信息 單元總數(shù)=%d 結(jié)點總數(shù)=%d 約束總數(shù)=%d n為輸出到FP2文件中時的格式，其后的變量表NELEM,NPOIN,NVFIX依次將其代表的數(shù)值輸出到%d所對應(yīng)的位置。%-4 調(diào)用 讀取初始數(shù)據(jù) 函數(shù)，生成結(jié)構(gòu)信息 LNODS,COORD,FPOIN,FPRES,FIXED=ele_INITIALDATA (FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX);%-5 調(diào)用 總剛計算 函數(shù)，生成單剛并集成總剛 HK = ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG);%-6調(diào)用 荷載計算 函數(shù)，由結(jié)點

8、荷載與非結(jié)點荷載生成總荷載向量列陣 FORCE = ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,FPRES,LNODS,COORD);%-7調(diào)用 邊界條件處理 函數(shù)，總剛、總荷載進行邊界條件處理 HK,FORCE = ele_BOUNDARY(NVFIX,FIXED,HK,FORCE);%-8 方程求解 DISP=zeros(NPOIN,1); %建立并初始化結(jié)構(gòu)的結(jié)點位移列矩陣 DISP=HKFORCE; %計算出結(jié)構(gòu)所有的結(jié)點位移%-9調(diào)用 單元桿端力計算 函數(shù)，求單元桿端力 FE = ele_MOMENTS(FP2,NPOIN,DISP,NELEM,LNODS

9、,COORD,YOUNG,NFPRES,FPRES);%*10 關(guān)閉 輸出數(shù)據(jù)文件 fclose(FP2);%*% 讀取初始數(shù)據(jù) 函數(shù)ele_INITIALDATA%* % 入口參數(shù)：FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX% 出口參數(shù)：LNODS,COORD,FPOIN,FPRES,FIXED%-function LNODS,COORD,FPOIN,FPRES,FIXED=ele_INITIALDATA(FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX)l 說明：FP1 文件句柄，指示輸入數(shù)據(jù)文件；FP2 文件句柄，指示輸出

10、數(shù)據(jù)文件；%-% 讀取結(jié)構(gòu)信息 LNODS=fscanf(FP1,%f,3,NELEM); l 說明：建立LNODS矩陣，該矩陣指出了每一單元的連接信息和慣性矩。矩陣的每一行針對每一單元，共計 NELEM；每一列相應(yīng)為單元左結(jié)點號（編碼）、（編碼）、慣性矩。命令中，3,NELEM表示讀取NELEM行3列數(shù)據(jù)賦值給LNODS矩陣。顯然，LNODS(i,1:3)依次表示i單元的左結(jié)點號、右結(jié)點號、慣性矩。從這種定義可見，每一單元的慣性矩均為常數(shù)，均為等截面直桿。 %- fprintf(FP2,單元號 左結(jié)點號 右結(jié)點號 慣 性 矩 n); for i=1:NELEM fprintf(FP2, %3

11、d %3d %3d %10.2e n,i,LNODS(i,:); endl 說明：在輸出文件FP2中顯示輸入的單元連接等信息。從1單元到NELEM單元進行循環(huán)，并依次輸出每一單元的連接等信息。%-COORD=fscanf(FP1,%f,NPOIN); % 坐標(biāo)： x坐標(biāo)（共計 NPOIN 組）l 說明：建立COORD矩陣，該矩陣用來存儲各結(jié)點x方向的坐標(biāo)值。連續(xù)梁結(jié)構(gòu)各結(jié)點均分布在x軸，以1結(jié)點為起始點順序編碼，因此表征各結(jié)點位置時僅需存儲其x方向的坐標(biāo)即可。從FP1文件中讀取全部結(jié)點個數(shù)NPOIN的坐標(biāo)值。COORD(i)表示第i個結(jié)點的x坐標(biāo)。%-fprintf(FP2, 結(jié)點定義數(shù)據(jù)：結(jié)

12、點號 x 坐 標(biāo) n);for i=1:NPOIN fprintf(FP2, %3d %6.2fn,i,COORD(i);endl 說明：在輸出文件FP2中顯示輸入的結(jié)點坐標(biāo)信息。從1結(jié)點到NPOIN結(jié)點進行循環(huán)，并依次輸出每一結(jié)點的坐標(biāo)信息。%-FPOIN=fscanf(FP1,%f,2,NFPOIN); l 說明：建立FPOIN矩陣，該矩陣用來存儲直接作用在結(jié)點上的荷載信息。從FP1文件讀取NFPOIN行2列數(shù)據(jù)，賦值給FPOIN矩陣。每一行的第一列表示荷載作用的結(jié)點號；第二列表示荷載的數(shù)值大小。連續(xù)梁結(jié)構(gòu)每一結(jié)點僅有1個轉(zhuǎn)角位移，相應(yīng)地直接作用在結(jié)點上的荷載為力偶，以順時針轉(zhuǎn)動為正。FP

13、OIN(i,1)表示第i個直接結(jié)點荷載作用的結(jié)點號，F(xiàn)POIN(i,2)表示第i個直接結(jié)點荷載的數(shù)值大小。若控制數(shù)據(jù)NFPOIN等于零，則此行命令不執(zhí)行。%-if NFPOIN0 fprintf(FP2, 結(jié)點荷載數(shù)據(jù)：結(jié)點號 M力偶 n); for i=1:NFPOIN fprintf(FP2, %3d %3d %6.2fn,i, FPOIN(i,:); endendl 說明：首先判斷控制數(shù)據(jù)NFPOIN是否為零，若為零則不需輸出數(shù)據(jù)；若不為零，表明本算例有直接作用在結(jié)點上的荷載，需要在輸出文件FP2中顯示輸入的結(jié)點荷載信息。從第1到第NFPOIN個直接結(jié)點荷載進行循環(huán)，并依次輸出每個直接結(jié)

14、點荷載的順序號、結(jié)點號和荷載大小。%-FPRES=fscanf(FP1,%f,3,NFPRES);l 說明：建立FPRES矩陣，該矩陣用來存儲非結(jié)點荷載信息。從FP1文件讀取NFPRES行3列數(shù)據(jù)，賦值給FPRES矩陣。每一行的第一列表示非結(jié)點荷載作用的單元號；第二列表示荷載類型；第三列表示荷載的數(shù)值大小。連續(xù)梁結(jié)構(gòu)的局部坐標(biāo)系x軸以向右為正、y軸以向上為正、轉(zhuǎn)動方向以順時針為正。FPRES(i,1)表示第i個非結(jié)點荷載作用的單元號，F(xiàn)PRES(i,2)表示第i個非結(jié)點荷載的類型（類型定義在函數(shù)ele_FORCE中），F(xiàn)PRES(i,3)表示第i個非結(jié)點荷載的數(shù)值大小。若控制數(shù)據(jù)NFPRES等

15、于零，則此行命令不執(zhí)行。%-if NFPRES0 fprintf(FP2, 非結(jié)點荷載數(shù)據(jù)：荷載號 單元號 荷載類型 荷載大小n); for i=1:NFPRES fprintf(FP2, %3d %3d %6.2f %6.2fn,i,FPRES(i,:); endendl 說明：首先判斷控制數(shù)據(jù)NFPRES是否為零，若為零則不需輸出數(shù)據(jù)；若不為零，表明本算例有非結(jié)點荷載，需要在輸出文件FP2中顯示輸入的非結(jié)點荷載信息。從第1到第NFPRES個非結(jié)點荷載進行循環(huán)，并依次輸出每個非結(jié)點荷載作用的荷載序號、單元號、荷載類型和荷載大小。%-FIXED=fscanf(FP1,%f,NVFIX);l 說

16、明：建立FIXED矩陣，該矩陣用來存儲零位移對應(yīng)的位移編碼信息（即有限制轉(zhuǎn)動的約束使結(jié)點的轉(zhuǎn)角位移為零，對應(yīng)的結(jié)點位移編碼。對于連續(xù)梁結(jié)構(gòu)即固定支座對應(yīng)的結(jié)點編碼）。從FP1文件讀取NVFIX行數(shù)據(jù)，賦值給FIXED列矩陣。每一行表示零位移對應(yīng)的結(jié)點位移編碼。FIXED(i)表示第i個零位移對應(yīng)的結(jié)點位移編碼。若控制數(shù)據(jù)NFPRES等于零，則此行命令不執(zhí)行。%-if NVFIX0 fprintf(FP2, 零位移約束數(shù)據(jù)： 位移編號 位移編號n); fprintf(FP2, %3d %3dn,FIXED(:); endl 說明：首先判斷控制數(shù)據(jù)NVFIX是否為零，若為零則不需輸出數(shù)據(jù)；若不為零

17、，表明本算例有零位移，需要在輸出文件FP2中顯示零位移約束數(shù)據(jù)。從第1到第NVFIX個零位移進行循環(huán)，依次輸出每個零位移對應(yīng)的結(jié)點號。%- returnl 說明： return表明主程序終止。%*% 計算總剛矩陣 函數(shù)ele_HK%* % 入口參數(shù)：結(jié)點數(shù)、單元數(shù)、單元信息數(shù)組、結(jié)點坐標(biāo)、彈性模量% 出口參數(shù)：整體剛度矩陣function HK = ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG)%-HK=zeros(NPOIN,NPOIN); % 生成總剛矩陣并清零% 生成單剛并組成總剛for i=1:NELEM % 對單元個數(shù)循環(huán) % 調(diào)用生成局部單剛(局部坐標(biāo))

18、 函數(shù) EK=ele_EK(i,LNODS,COORD,YOUNG); HK(i:i+1,i:i+1) = HK(i:i+1,i:i+1)+EK;endl 說明：對每一單元順序循環(huán)，依次計算每一單元的局部單元剛度矩陣。EK為一個臨時變量，計算得到i單元的單剛EK后，通過下一條語句集成入總剛矩陣中。在下一次循環(huán)中，則計算i+1單元的單剛，并賦值給EK，再進行集成。顯然，隨著循環(huán)的不斷進行，EK對應(yīng)著不同單元的單剛，值不斷發(fā)生改變。單剛到總剛的集成是利用了連續(xù)梁結(jié)構(gòu)的特點，以分塊的形式放入總剛矩陣中。例如：對于第i個單元，其四個單剛元素分別放入總剛矩陣的第（i,i）、（i,i+1）、（i+1,i）

19、、（i+1,i+1）四個位置，因此通過HK的下標(biāo)(i:i+1,i:i+1)來表示相應(yīng)的元素位置。由于相鄰單元的單剛集成時，有元素疊加的情況，因此相應(yīng)的總剛矩陣元素的疊加采用了HK(i:i+1,i:i+1) = HK(i:i+1,i:i+1)+EK的表達(dá)形式。%-return%*% 計算單元剛度矩陣 函數(shù)ele_EK%* % 入口參數(shù)：單元號、單元信息、結(jié)點坐標(biāo)、彈性模量% 出口參數(shù)：局部單元剛度矩陣l 說明：在ele_HK函數(shù)中調(diào)用本函數(shù)時，彈性模量的入口參數(shù)符號為YOUNG，而在本函數(shù)的輸入?yún)?shù)表中，符號為E。函數(shù)調(diào)用時可以采用此種用法。%-function EK=ele_EK(i,LNOD

20、S,COORD,E)%- NL=LNODS(i,1); %左結(jié)點號 NR=LNODS(i,2); %右結(jié)點號 L=COORD(NR)-COORD(NL); %單元長度 I=LNODS(i,3); %慣性矩i=E*I/L;l 說明：利用已知信息，計算每一單元的線剛度。LNODS(i,1)為i單元左結(jié)點的結(jié)點號賦值給NL；LNODS(i,2)為i單元右結(jié)點的結(jié)點號賦值給NR。COORD(NR)為第NR個結(jié)點的x坐標(biāo)；COORD(NL)為第NL個結(jié)點的x坐標(biāo)；二者的差為單元的長度。LNODS(i,1)為i單元的慣性矩。 %- % 生成單剛(局部坐標(biāo)) 右手坐標(biāo)系EK = 4*i 2*i; 2*i 4

21、*i; %-return%*% 生成總荷載向量 函數(shù)ele_FORCE%* % 入口參數(shù)：結(jié)點數(shù)，結(jié)點荷載個數(shù)，結(jié)點荷載信息，非結(jié)點荷載個數(shù)，非結(jié)點荷載信息，單元信息，結(jié)點坐標(biāo)% 出口參數(shù)：單元固端力左右兩端的桿端彎矩function FORCE = ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,FPRES,LNODS,COORD)%- FORCE=zeros(NPOIN,1); % 生成總荷載向量并清零 for i=1:NFPOIN % 對結(jié)點荷載個數(shù)進行循環(huán) FORCE(FPOIN(i,1)=FORCE(FPOIN(i,1)+FPOIN(i,2); % 取出結(jié)點荷

22、載 endl 說明：建立并初始化總荷載向量FORCE，其行數(shù)與總剛矩陣、總位移列矩陣相同，即與總結(jié)點數(shù)NPOIN相同（對于連續(xù)梁結(jié)構(gòu)）。對結(jié)點荷載個數(shù)進行循環(huán)，將每個結(jié)點荷載的數(shù)值放入所在結(jié)點號對應(yīng)的總荷載向量FORCE對應(yīng)的行處FPOIN(i,1)為第i個結(jié)點荷載所作用的結(jié)點號；FPOIN(i,2) 為第i個結(jié)點荷載的大??；FORCE(FPOIN(i,1)為總荷向量對應(yīng)結(jié)點號處的荷載。%-for i=1:NFPRES % 對非結(jié)點荷載個數(shù)進行循環(huán) F0=ele_FPRES(i,FPRES,LNODS,COORD); % 計算單元固端力 % 對單元局部桿端力要進行坐標(biāo)轉(zhuǎn)換 ele=FPRES(

23、i,1); % 取荷載所在的單元號 NL=LNODS(ele,1); NR=LNODS(ele,2); %單元的左右結(jié)點號 FORCE(NL)=FORCE(NL)-F0(1); %將固端力變成等效結(jié)點荷載 FORCE(NR)=FORCE(NR)-F0(2); %固端力與等效荷載符號相反endl 說明：計算由非結(jié)點荷載引起的等效結(jié)點荷載，并集成入已有的總荷向量中。對非結(jié)點荷載個數(shù)進行循環(huán)，即有幾個非結(jié)點荷載循環(huán)幾次，若沒有非結(jié)點荷載，則本循環(huán)語句不執(zhí)行。調(diào)用計算非結(jié)點荷載引起的單元固端力函數(shù)，計算出非結(jié)點荷載所在單元的兩端固端力。FPRES(i,1)為第i個非結(jié)點荷載所在的單元號賦值給ele；L

24、NODS(ele,1)為ele單元左結(jié)點的結(jié)點號賦值給NL；LNODS(ele,2)為ele單元右結(jié)點的結(jié)點號賦值給NR。FORCE(NL)、FORCE(NR)分別為單元的左、右桿端的桿端力；-F0(1)、-F0(2)為非結(jié)點荷載引起的左、右端固端力的負(fù)值，即等效結(jié)點荷載。FORCE(NL)=FORCE(NL)-F0(1)為直接作用在單元左端結(jié)點上的荷載與左端等效結(jié)點荷載的疊加；FORCE(NR)=FORCE(NR)-F0(2) 為直接作用在單元右端結(jié)點上的荷載與右端等效結(jié)點荷載的疊加。%-return%*% 計算非結(jié)點荷載引起的單元固端力 函數(shù) ele_FPRES% 符號規(guī)定：正方向為：X向

25、右 Y向上 M順時針%* % 入口參數(shù)：荷載序號，非結(jié)點荷載信息，單元信息，結(jié)點坐標(biāo)% 出口參數(shù)：單元固端力（左右兩端固端彎矩）%-function F0=ele_FPRES(i,FPRES,LNODS,COORD) %- ele=FPRES(i,1); % 取荷載所在的單元號 G=FPRES(i,3); % 單元荷載大小 NL=LNODS(ele,1); NR=LNODS(ele,2); % 單元的左右結(jié)點號L=COORD(NR)-COORD(NL); % 單元長度l 說明：將已知非結(jié)點荷載信息分別賦值給相應(yīng)的變量。FPRES(i,1)為第i個非結(jié)點荷載作用的單元號；G=FPRES(i,3)

26、為第i個非結(jié)點荷載的大小；LNODS(ele,1)為ele單元左結(jié)點的結(jié)點號賦值給NL；LNODS(ele,2)為ele單元右結(jié)點的結(jié)點號賦值給NR。COORD(NR)為第NR個結(jié)點的x坐標(biāo)；COORD(NL)為第NL個結(jié)點的x坐標(biāo)；二者的差為單元的長度。%- F0=0;0; % 單元固端彎矩清零 switch FPRES(i,2) case 1 % 橫向均布荷載 F0(1)=-G*L2/12.0; F0(2)= G*L2/12.0; case 2 % 橫向集中力 F0(1)=-G*L/8; F0(2)= G*L/8;endl 說明：根據(jù)荷載類型，求解單元桿端的桿端彎矩。FPRES(i,2)為

27、第i個非結(jié)點荷載的類型。該類型是在下面的switch判斷模塊中定義的。當(dāng)FPRES(i,2)1時，為橫向均布荷載（橫向指垂直于桿件軸線）；當(dāng)FPRES(i,2)2時，為橫向集中力。即荷載類型分別為1、2，這兩個數(shù)值需要在輸入非結(jié)點荷載信息時針對荷載的類型進行輸入。每種非結(jié)點荷載作用下的單元固端力的計算公式是按照結(jié)構(gòu)力學(xué)教材得到的。其中，F(xiàn)0(1)表示單元左端的固端彎矩；F0(2)表示單元右端的固端彎矩。%-return%*% 總剛、總荷載進行邊界條件處理 函數(shù)ele_BOUNDARY%* % 入口參數(shù)：零位移個數(shù)，約束信息，總剛矩陣，總荷載矩陣% 出口參數(shù)：總剛矩陣，總荷載矩陣%-functi

28、on HK,FORCE = ele_BOUNDARY(NVFIX,FIXED,HK,FORCE)%-for j=1:NVFIX % 對約束個數(shù)進行循環(huán) N1=FIXED(j); %找出約束所在的結(jié)點位移編碼 HK(N1,N1)=HK(N1,N1)*1E10; %乘大數(shù)法endl 說明：采用乘大數(shù)法處理邊界條件，僅需對總剛元素進行處理即可。FIXED(j)為第j個約束所在的結(jié)點位移編碼（連續(xù)梁即結(jié)點號）賦值給N1。對約束個數(shù)進行循環(huán)，對每j個約束處理時，需對總剛矩陣中的第j行第j列的主對角線元素HK(N1,N1)乘以一個大數(shù)，如。%-若采用置0置1法處理邊界條件，其相應(yīng)程序為for j=1:NV

29、FIX N1=FIXED(j); %找出約束所在的結(jié)點位移編碼HK(:,N1)=0; %將總剛矩陣HK(:,N1)的第N1列所有元素賦值為0。 HK(N1,:)=0; %將總剛矩陣HK(N1,:)的第N1行所有元素賦值為0。 HK(N1,N1)=1; %將總剛矩陣HK(N1,N1)主對角線元素賦值為1。 FORCE(N1)=0; %將總荷載矩陣FORCE(N1)的第N1行元素賦值為0end%-return%*% 計算單元桿端力 函數(shù)ele_MOMENTS%* % 入口參數(shù)：輸出文件句柄，結(jié)點數(shù)，結(jié)構(gòu)位移列陣，單元數(shù)，單元信息，結(jié)點坐標(biāo)，彈性模量，非結(jié)點荷載個數(shù)，非結(jié)點荷載信息% 出口參數(shù)：單元

30、桿端力左右兩端的桿端彎矩function FE = ele_MOMENTS(FP2,NPOIN,DISP,NELEM, LNODS,COORD, YOUNG, NFPRES, FPRES)%-fprintf(FP2,n 計 算 結(jié) 果nn); fprintf(FP2, 輸出結(jié)構(gòu)結(jié)點位移 結(jié)點號 轉(zhuǎn)角n);for i=1:NPOIN fprintf(FP2, %3d %11.3e n,i,DISP(i);endl 說明：在輸出文件FP2中顯示提示信息：計算結(jié)果。通過循環(huán)按照結(jié)點位移編碼（連續(xù)梁結(jié)構(gòu)即結(jié)點編碼）的順序，在文件FP2中輸出所有結(jié)點位移的數(shù)值。%-EDISP=zeros(2,1); %

31、 單元位移列向量清零fprintf(FP2, 輸出單元桿端內(nèi)力，以順時針為正n);for i=1:NELEM % 對單元個數(shù)進行循環(huán)，依次計算每一單元的桿端力 for j=1:2 N1=LNODS(i,j); % 單元桿端結(jié)點號 EDISP(j)=DISP(N1); % 取N1端的單元位移列向量 endl 說明：通過循環(huán)，依次從結(jié)構(gòu)位移列陣中對號，賦值給第i個單元的單元位移向量EDISP；j=1時N1為左端，j=2時N1為右端。 %- EK=ele_EK(i,LNODS,COORD,YOUNG); % 生成局部單剛(局部坐標(biāo)) FE=EK*EDISP; %計算局部坐標(biāo)系下的由結(jié)點位移產(chǎn)生的桿端

32、力l 說明：計算第i個單元的單剛EK；由單剛乘以單元位移得到單元僅由結(jié)點位移產(chǎn)生的桿端力。 %- for j=1:NFPRES if FPRES(j,1) = i %成立時，當(dāng)前單元上有非結(jié)點荷載 F0=ele_FPRES(j,FPRES,LNODS,COORD); %計算單元固端力 FE=FE+F0; %考慮由非結(jié)點荷載引起的桿端力，直接疊加 end end l 說明：對非結(jié)點個數(shù)進行循環(huán)。 FPRES(j,1)為第j個非結(jié)點荷載作用的單元號，若等于當(dāng)前循環(huán)的單元號i，則表明當(dāng)前單元上有非結(jié)點荷載。 調(diào)用計算單元固端力函數(shù)，得到非結(jié)點荷載引起的單元固端力。 將桿端位移引起的桿端力和非結(jié)點荷載

33、引起的單元固端力疊加，得到當(dāng)前單元的全部桿端力。 %- fprintf(FP2, 單元(%d)桿端力 彎 矩n,i); fprintf(FP2, 左端%3d結(jié)點 %7.2fn,LNODS(i,1),FE(1); fprintf(FP2, 右端%3d結(jié)點 %7.2fn,LNODS(i,2),FE(2); l 說明：輸出單元的結(jié)點號，對應(yīng)的桿端彎矩。end %對單元循環(huán)結(jié)束。%-return3.2 程序應(yīng)用舉例【例1】利用連續(xù)梁程序計算圖2-9所示的連續(xù)梁結(jié)構(gòu)。353563圖2-9%-輸入數(shù)據(jù)文件input.txt為：3 4 1 1 2 1.01 2 1.02 3 1.03 4 1.00.06.012.017.02 6.0