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

matlab連續(xù)梁程序的編制與使用matlab連續(xù)梁程序的編制與使用matlab連續(xù)梁程序的編制與使用xxx公司matlab連續(xù)梁程序的編制與使用文件編號：文件日期：修訂次數：第1.0次更改批準審核制定方案設計，管理制度連續(xù)梁程序的編制與使用結構力學中的矩陣位移法是隨著電子計算機進入結構力學領域中而產生的一種方法，而Matlab語言正是進行矩陣運算的強大工具，因此，用Matlab語言編寫結構力學程序有更大的優(yōu)越性。本章將詳細介紹如何利用Matlab語言編制連續(xù)梁結構的計算程序。矩陣位移法的解題思路是將結構離散為單元（桿件），建立單元桿端力與桿端位移之間的關系－單元剛度方程；再將各單元集成為原結構，在滿足變形連續(xù)條件和平衡條件時，建立整體剛度方程；在邊界條件處理完畢后，由整體剛度方程解出節(jié)點位移，進而求出結構內力。用矩陣位移法計算連續(xù)梁的步驟如下：整理原始數據，如材料性質、荷載條件、約束條件等，并進行編碼：單元編碼、結點編碼、結點位移編碼、選取坐標系。建立局部坐標系下的單元剛度矩陣。建立整體坐標系下的單元剛度矩陣。集成總剛。建立整體結構的等效節(jié)點荷載和總荷載矩陣邊界條件處理。解方程，求出節(jié)點位移。求出各單元的桿端內力。圖3-1程序流程圖實際上，圖3-1程序流程圖程序說明%*******************************************************************%矩陣位移法解連續(xù)梁主程序%*******************************************************************功能：運用矩陣位移法解連續(xù)梁的基本原理編制的計算主程序。基本思想：結點（結點位移）編碼默認為從左至右，從1開始順序進行；桿端彎矩的方向默認為逆時針。荷載類型：可計算結點荷載，每單元作用的跨中集中力和均布荷載。說明：主程序的作用是通過賦值語句、讀取和寫入文件、函數調用等完成算法的全過程，即實現程序流程圖的程序表達。%-----------------------------------------------------------------------------------------------------1程序準備 formatshorte %設定輸出類型 clearall %清除所有已定義變量clc %清屏說明：formatshorte －設定計算過程中顯示在屏幕上的數字類型為短格式、科學計數法；clearall －清除所有已定義變量，目的是在本程序的運行過程中，不會發(fā)生變量名相同等可能使計算出錯的情況；clc －清屏，使屏幕在本程序運行開始時%-----------------------------------------------------------------------------------------------------2打開文件 FP1=fopen('','rt'); %打開輸入數據文件存放初始數據 FP2=fopen('','wt'); %打開輸出數據文件存放計算結果說明：FP1=fopen('','rt'); －打開已存在的輸入數據文件，且設置其為只讀格式，使程序在執(zhí)行過程中不能改變輸入文件中的數值，并用文件句柄FP1來FP2=fopen('','wt'); －打開輸出數據文件，該文件不存在時，通過此命令創(chuàng)建新文件，該文件存在時則將原有內容全部刪除。該文件設置為可寫格式，可在程序執(zhí)行過程中向輸出文件寫入數據。%-----------------------------------------------------------------------------------------------------3讀入程序控制信息 NELEM =fscanf(FP1,'%d',1);%單元個數（單元編碼總數） NPOIN =fscanf(FP1,'%d',1);%結點個數（結點編碼總數） NVFIX =fscanf(FP1,'%d',1); ％約束個數（零位移總數） NFPOIN=fscanf(FP1,'%d',1); ％結點荷載個數（作用在結點上集中力偶總數） NFPRES=fscanf(FP1,'%d',1); ％非結點荷載數（作用在單元上分布荷載總數） YOUNG=fscanf(FP1,'%f',1); ％彈性模量說明：從輸入文件FP1中讀入單元個數，結點個數，約束個數，結點荷載個數，非結點荷載個數，彈性模量；程序中彈性模量僅輸入了一個值，表明本程序僅能求解一種材料構成的結構，如：鋼筋混凝土結構、鋼結構，不能求解鋼筋混凝土－鋼組合結構。采用了命令fscanf，其中：’%d’表示讀入整數格式；1表示讀取1個數。%------------------------------------------------------------------------------------------ fprintf(FP2,'\n結構初始數據\n\n'); fprintf(FP2,'單元總數=%d結點總數=%d約束總數=%d\n',NELEM,NPOIN,NVFIX); fprintf(FP2,'結點荷載個數=%d非結點荷載個數=%d彈性模量=%12.5g\n',NFPOIN,NFPRES,YOUNG);說明：在輸出文件FP2中顯示輸入的控制信息，便于程序執(zhí)行完畢后，從輸出文件中查找輸入錯誤。采用了命令fprintf，其中：\n為換行標志，表示換一行；%12.5g表示輸出12位且有5位小數的實數。括在引號中的信息'單元總數=%d結點總數=%d約束總數=%d\n'為輸出到FP2文件中時的格式，其后的變量表NELEM,NPOIN,NVFIX依次將其代表的數值輸出到%d所對應的位置。%-----------------------------------------------------------------------------------------------------4調用讀取初始數據函數，生成結構信息[LNODS,COORD,FPOIN,FPRES,FIXED]=…ele_INITIALDATA(FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX);%-----------------------------------------------------------------------------------------------------5調用總剛計算函數，生成單剛并集成總剛[HK]=ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG);%-----------------------------------------------------------------------------------------------------6調用荷載計算函數，由結點荷載與非結點荷載生成總荷載向量列陣[FORCE]=…ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,FPRES,LNODS,COORD);%-----------------------------------------------------------------------------------------------------7調用邊界條件處理函數，總剛、總荷載進行邊界條件處理[HK,FORCE]=ele_BOUNDARY(NVFIX,FIXED,HK,FORCE);%-----------------------------------------------------------------------------------------------------8方程求解DISP=zeros(NPOIN,1);%建立并初始化結構的結點位移列矩陣DISP=HK\FORCE;%計算出結構所有的結點位移%-----------------------------------------------------------------------------------------------------9調用單元桿端力計算函數，求單元桿端力[FE]=ele_MOMENTS(FP2,NPOIN,DISP,NELEM,LNODS,…COORD,YOUNG,NFPRES,FPRES);%*******************************************************************10關閉輸出數據文件fclose(FP2);%*******************************************************************%讀取初始數據函數ele_INITIALDATA%*******************************************************************%入口參數：FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX%出口參數：LNODS,COORD,FPOIN,FPRES,FIXED%------------------------------------------------------------------------------------------function[LNODS,COORD,FPOIN,FPRES,FIXED]=…ele_INITIALDATA(FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX)說明：FP1－文件句柄，指示輸入數據文件；FP2－文件句柄，指示輸出數據文件；%------------------------------------------------------------------------------------------%讀取結構信息 LNODS=fscanf(FP1,'%f',[3,NELEM])';說明：建立LNODS矩陣，該矩陣指出了每一單元的連接信息和慣性矩。矩陣的每一行針對每一單元，共計NELEM；每一列相應為單元左結點號（編碼）、（編碼）、慣性矩。命令中，[3,NELEM]表示讀取NELEM行3列數據賦值給LNODS矩陣。顯然，LNODS(i,1:3)依次表示i單元的左結點號、右結點號、慣性矩。從這種定義可見，每一單元的慣性矩均為常數，均為等截面直桿。%------------------------------------------------------------------------------------------fprintf(FP2,'單元號左結點號右結點號慣性矩\n'); fori=1:NELEMfprintf(FP2,'%3d%3d%3d%\n',i,LNODS(i,:)); end說明：在輸出文件FP2中顯示輸入的單元連接等信息。從1單元到NELEM單元進行循環(huán)，并依次輸出每一單元的連接等信息。%---------------------------------------------------------------------------------------------------- COORD=fscanf(FP1,'%f',[NPOIN])';%坐標：x坐標（共計NPOIN組）說明：建立COORD矩陣，該矩陣用來存儲各結點x方向的坐標值。連續(xù)梁結構各結點均分布在x軸，以1結點為起始點順序編碼，因此表征各結點位置時僅需存儲其x方向的坐標即可。從FP1文件中讀取全部結點個數NPOIN的坐標值。COORD(i)表示第i個結點的x坐標。%------------------------------------------------------------------------------------------ fprintf(FP2,'結點定義數據：結點號x坐標\n'); fori=1:NPOINfprintf(FP2,'%3d%6.2f\n',i,COORD(i)); end說明：在輸出文件FP2中顯示輸入的結點坐標信息。從1結點到NPOIN結點進行循環(huán)，并依次輸出每一結點的坐標信息。%-----------------------------------------------------------------------------------------------------FPOIN=fscanf(FP1,'%f',[2,NFPOIN])';說明：建立FPOIN矩陣，該矩陣用來存儲直接作用在結點上的荷載信息。從FP1文件讀取NFPOIN行2列數據，賦值給FPOIN矩陣。每一行的第一列表示荷載作用的結點號；第二列表示荷載的數值大小。連續(xù)梁結構每一結點僅有1個轉角位移，相應地直接作用在結點上的荷載為力偶，以順時針轉動為正。FPOIN(i,1)表示第i個直接結點荷載作用的結點號，FPOIN(i,2)表示第i個直接結點荷載的數值大小。若控制數據NFPOIN等于零，則此行命令不執(zhí)行。%------------------------------------------------------------------------------------------ ifNFPOIN>0fprintf(FP2,'結點荷載數據：結點號M力偶\n'); fori=1:NFPOINfprintf(FP2,'%3d%3d%6.2f\n',i,FPOIN(i,:)); endend說明：首先判斷控制數據NFPOIN是否為零，若為零則不需輸出數據；若不為零，表明本算例有直接作用在結點上的荷載，需要在輸出文件FP2中顯示輸入的結點荷載信息。從第1到第NFPOIN個直接結點荷載進行循環(huán)，并依次輸出每個直接結點荷載的順序號、結點號和荷載大小。%-----------------------------------------------------------------------------------------------------FPRES=fscanf(FP1,'%f',[3,NFPRES])';說明：建立FPRES矩陣，該矩陣用來存儲非結點荷載信息。從FP1文件讀取NFPRES行3列數據，賦值給FPRES矩陣。每一行的第一列表示非結點荷載作用的單元號；第二列表示荷載類型；第三列表示荷載的數值大小。連續(xù)梁結構的局部坐標系x軸以向右為正、y軸以向上為正、轉動方向以順時針為正。FPRES(i,1)表示第i個非結點荷載作用的單元號，FPRES(i,2)表示第i個非結點荷載的類型（類型定義在函數ele_FORCE中），FPRES(i,3)表示第i個非結點荷載的數值大小。若控制數據NFPRES等于零，則此行命令不執(zhí)行。%------------------------------------------------------------------------------------------ ifNFPRES>0fprintf(FP2,'非結點荷載數據：荷載號單元號荷載類型荷載大小\n'); fori=1:NFPRESfprintf(FP2,'%3d%3d%6.2f%6.2f\n',i,FPRES(i,:)); endend說明：首先判斷控制數據NFPRES是否為零，若為零則不需輸出數據；若不為零，表明本算例有非結點荷載，需要在輸出文件FP2中顯示輸入的非結點荷載信息。從第1到第NFPRES個非結點荷載進行循環(huán)，并依次輸出每個非結點荷載作用的荷載序號、單元號、荷載類型和荷載大小。%-----------------------------------------------------------------------------------------------------FIXED=fscanf(FP1,'%f',NVFIX)';說明：建立FIXED矩陣，該矩陣用來存儲零位移對應的位移編碼信息（即有限制轉動的約束使結點的轉角位移為零，對應的結點位移編碼。對于連續(xù)梁結構即固定支座對應的結點編碼）。從FP1文件讀取NVFIX行數據，賦值給FIXED列矩陣。每一行表示零位移對應的結點位移編碼。FIXED(i)表示第i個零位移對應的結點位移編碼。若控制數據NFPRES等于零，則此行命令不執(zhí)行。%----------------------------------------------------------------------------------------------- ifNVFIX>0 fprintf(FP2,'零位移約束數據：位移編號位移編號\n'); fprintf(FP2,'%3d%3d\n',FIXED(:));end說明：首先判斷控制數據NVFIX是否為零，若為零則不需輸出數據；若不為零，表明本算例有零位移，需要在輸出文件FP2中顯示零位移約束數據。從第1到第NVFIX個零位移進行循環(huán)，依次輸出每個零位移對應的結點號。%-----------------------------------------------------------------------------------------------------return說明：return表明主程序終止。%*******************************************************************%計算總剛矩陣函數ele_HK%*******************************************************************%入口參數：結點數、單元數、單元信息數組、結點坐標、彈性模量%出口參數：整體剛度矩陣function[HK]=ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG)%----------------------------------------------------------------------------------------------------- HK=zeros(NPOIN,NPOIN);%生成總剛矩陣并清零 %生成單剛并組成總剛 fori=1:NELEM%對單元個數循環(huán)%調用生成局部單剛(局部坐標)函數EK=ele_EK(i,LNODS,COORD,YOUNG);HK(i:i+1,i:i+1)=HK(i:i+1,i:i+1)+EK; end說明：對每一單元順序循環(huán)，依次計算每一單元的局部單元剛度矩陣。EK為一個臨時變量，計算得到i單元的單剛EK后，通過下一條語句集成入總剛矩陣中。在下一次循環(huán)中，則計算i+1單元的單剛，并賦值給EK，再進行集成。顯然，隨著循環(huán)的不斷進行，EK對應著不同單元的單剛，值不斷發(fā)生改變。單剛到總剛的集成是利用了連續(xù)梁結構的特點，以分塊的形式放入總剛矩陣中。例如：對于第i個單元，其四個單剛元素分別放入總剛矩陣的第（i,i）、（i,i+1）、（i+1,i）、（i+1,i+1）四個位置，因此通過HK的下標(i:i+1,i:i+1)來表示相應的元素位置。由于相鄰單元的單剛集成時，有元素疊加的情況，因此相應的總剛矩陣元素的疊加采用了HK(i:i+1,i:i+1)=HK(i:i+1,i:i+1)+EK的表達形式。%-----------------------------------------------------------------------------------------------------return%*******************************************************************%計算單元剛度矩陣函數ele_EK%*******************************************************************%入口參數：單元號、單元信息、結點坐標、彈性模量%出口參數：局部單元剛度矩陣說明：在ele_HK函數中調用本函數時，彈性模量的入口參數符號為YOUNG，而在本函數的輸入參數表中，符號為E。函數調用時可以采用此種用法。%-----------------------------------------------------------------------------------------------------functionEK=ele_EK(i,LNODS,COORD,E)%-----------------------------------------------------------------------------------------------------NL=LNODS(i,1);%左結點號NR=LNODS(i,2);%右結點號L=COORD(NR)-COORD(NL);%單元長度I=LNODS(i,3);%慣性矩i=E*I/L;說明：利用已知信息，計算每一單元的線剛度。LNODS(i,1)為i單元左結點的結點號賦值給NL；LNODS(i,2)為i單元右結點的結點號賦值給NR。COORD(NR)為第NR個結點的x坐標；COORD(NL)為第NL個結點的x坐標；二者的差為單元的長度。LNODS(i,1)為i單元的慣性矩。%----------------------------------------------------------------------------------------------%生成單剛(局部坐標)右手坐標系EK=[4*i2*i;2*i4*i];%----------------------------------------------------------------------------------------------return%*******************************************************************%生成總荷載向量函數ele_FORCE%*******************************************************************%入口參數：結點數，結點荷載個數，結點荷載信息，非結點荷載個數，非結點荷載信息，單元信息，結點坐標%出口參數：單元固端力—左右兩端的桿端彎矩function[FORCE]=ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,…FPRES,LNODS,COORD)%----------------------------------------------------------------------------------------------------FORCE=zeros(NPOIN,1);%生成總荷載向量并清零fori=1:NFPOIN%對結點荷載個數進行循環(huán)FORCE(FPOIN(i,1))=FORCE(FPOIN(i,1))+FPOIN(i,2);%取出結點荷載end說明：建立并初始化總荷載向量FORCE，其行數與總剛矩陣、總位移列矩陣相同，即與總結點數NPOIN相同（對于連續(xù)梁結構）。對結點荷載個數進行循環(huán)，將每個結點荷載的數值放入所在結點號對應的總荷載向量FORCE對應的行處FPOIN(i,1)為第i個結點荷載所作用的結點號；FPOIN(i,2)為第i個結點荷載的大小；FORCE(FPOIN(i,1))為總荷向量對應結點號處的荷載。%---------------------------------------------------------------------------------------------------- fori=1:NFPRES%對非結點荷載個數進行循環(huán)F0=ele_FPRES(i,FPRES,LNODS,COORD);%計算單元固端力%對單元局部桿端力要進行坐標轉換ele=FPRES(i,1);%取荷載所在的單元號NL=LNODS(ele,1);NR=LNODS(ele,2);%單元的左右結點號FORCE(NL)=FORCE(NL)-F0(1);%將固端力變成等效結點荷載FORCE(NR)=FORCE(NR)-F0(2);%固端力與等效荷載符號相反 end說明：計算由非結點荷載引起的等效結點荷載，并集成入已有的總荷向量中。對非結點荷載個數進行循環(huán)，即有幾個非結點荷載循環(huán)幾次，若沒有非結點荷載，則本循環(huán)語句不執(zhí)行。調用計算非結點荷載引起的單元固端力函數，計算出非結點荷載所在單元的兩端固端力。FPRES(i,1)為第i個非結點荷載所在的單元號賦值給ele；LNODS(ele,1)為ele單元左結點的結點號賦值給NL；LNODS(ele,2)為ele單元右結點的結點號賦值給NR。FORCE(NL)、FORCE(NR)分別為單元的左、右桿端的桿端力；-F0(1)、-F0(2)為非結點荷載引起的左、右端固端力的負值，即等效結點荷載。FORCE(NL)=FORCE(NL)-F0(1)為直接作用在單元左端結點上的荷載與左端等效結點荷載的疊加；FORCE(NR)=FORCE(NR)-F0(2)為直接作用在單元右端結點上的荷載與右端等效結點荷載的疊加。%----------------------------------------------------------------------------------------------------return%*******************************************************************%計算非結點荷載引起的單元固端力函數ele_FPRES%符號規(guī)定：正方向為：X向右Y向上M順時針%*******************************************************************%入口參數：荷載序號，非結點荷載信息，單元信息，結點坐標%出口參數：單元固端力（左右兩端固端彎矩）%----------------------------------------------------------------------------------------------------functionF0=ele_FPRES(i,FPRES,LNODS,COORD)%----------------------------------------------------------------------------------------------------ele=FPRES(i,1);%取荷載所在的單元號G=FPRES(i,3);%單元荷載大小NL=LNODS(ele,1);NR=LNODS(ele,2);%單元的左右結點號L=COORD(NR)-COORD(NL);%單元長度說明：將已知非結點荷載信息分別賦值給相應的變量。FPRES(i,1)為第i個非結點荷載作用的單元號；G=FPRES(i,3)為第i個非結點荷載的大??；LNODS(ele,1)為ele單元左結點的結點號賦值給NL；LNODS(ele,2)為ele單元右結點的結點號賦值給NR。COORD(NR)為第NR個結點的x坐標；COORD(NL)為第NL個結點的x坐標；二者的差為單元的長度。%----------------------------------------------------------------------------------------------------F0=[0;0];%單元固端彎矩清零switchFPRES(i,2)case1%橫向均布荷載F0(1)=-G*L^2/;F0(2)=G*L^2/;case2%橫向集中力F0(1)=-G*L/8;F0(2)=G*L/8;end說明：根據荷載類型，求解單元桿端的桿端彎矩。FPRES(i,2)為第i個非結點荷載的類型。該類型是在下面的switch判斷模塊中定義的。當FPRES(i,2)＝1時，為橫向均布荷載（橫向指垂直于桿件軸線）；當FPRES(i,2)＝2時，為橫向集中力。即荷載類型分別為1、2，這兩個數值需要在輸入非結點荷載信息時針對荷載的類型進行輸入。每種非結點荷載作用下的單元固端力的計算公式是按照結構力學教材得到的。其中，F0(1)表示單元左端的固端彎矩；F0(2)表示單元右端的固端彎矩。%----------------------------------------------------------------------------------------------------return%*******************************************************************%總剛、總荷載進行邊界條件處理函數ele_BOUNDARY%*******************************************************************%入口參數：零位移個數，約束信息，總剛矩陣，總荷載矩陣%出口參數：總剛矩陣，總荷載矩陣%----------------------------------------------------------------------------------------------------function[HK,FORCE]=ele_BOUNDARY(NVFIX,FIXED,HK,FORCE)%---------------------------------------------------------------------------------------------------- forj=1:NVFIX%對約束個數進行循環(huán)N1=FIXED(j);%找出約束所在的結點位移編碼HK(N1,N1)=HK(N1,N1)*1E10;%乘大數法 end說明：采用乘大數法處理邊界條件，僅需對總剛元素進行處理即可。FIXED(j)為第j個約束所在的結點位移編碼（連續(xù)梁即結點號）賦值給N1。對約束個數進行循環(huán)，對每j個約束處理時，需對總剛矩陣中的第j行第j列的主對角線元素HK(N1,N1)乘以一個大數，如。%----------------------------------------------------------------------------------------------------若采用置0置1法處理邊界條件，其相應程序為 forj=1:NVFIXN1=FIXED(j);%找出約束所在的結點位移編碼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行元素賦值為0 end%----------------------------------------------------------------------------------------------------return%*******************************************************************%計算單元桿端力函數ele_MOMENTS%*******************************************************************%入口參數：輸出文件句柄，結點數，結構位移列陣，單元數，單元信息，結點坐標，彈性模量，非結點荷載個數，非結點荷載信息%出口參數：單元桿端力—左右兩端的桿端彎矩function[FE]=ele_MOMENTS(FP2,NPOIN,DISP,NELEM,LNODS,…COORD,YOUNG,NFPRES,FPRES)%---------------------------------------------------------------------------------------------------- fprintf(FP2,'\n計算結果\n\n'); fprintf(FP2,'輸出結構結點位移結點號轉角\n'); fori=1:NPOINfprintf(FP2,'%3d%\n',i,DISP(i)); end說明：在輸出文件FP2中顯示提示信息：計算結果。通過循環(huán)按照結點位移編碼（連續(xù)梁結構即結點編碼）的順序，在文件FP2中輸出所有結點位移的數值。%----------------------------------------------------------------------------------------------------EDISP=zeros(2,1);%單元位移列向量清零fprintf(FP2,'輸出單元桿端內力，以順時針為正\n');fori=1:NELEM%對單元個數進行循環(huán)，依次計算每一單元的桿端力forj=1:2N1=LNODS(i,j);%單元桿端結點號EDISP(j)=DISP(N1);%取N1端的單元位移列向量end說明：通過循環(huán)，依次從結構位移列陣中對號，賦值給第i個單元的單元位移向量EDISP；j=1時N1為左端，j=2時N1為右端。%--------------------------------------------------------------------------------------EK=ele_EK(i,LNODS,COORD,YOUNG);%生成局部單剛(局部坐標)FE=EK*EDISP;%計算局部坐標系下的由結點位移產生的桿端力說明：計算第i個單元的單剛EK；由單剛乘以單元位移得到單元僅由結點位移產生的桿端力。%--------------------------------------------------------------------------------------forj=1:NFPRESifFPRES(j,1)==i%成立時，當前單元上有非結點荷載F0=ele_FPRES(j,FPRES,LNODS,COORD);%計算單元固端力FE=FE+F0;%考慮由非結點荷載引起的桿端力，直接疊加endend說明：對非結點個數進行循環(huán)。FPRES(j,1)為第j個非結點荷載作用的單元號，若等于當前循環(huán)的單元號i，則表明當前單元上有非結點荷載。調用計算單元固端力函數，得到非結點荷載引起的單元固端力。將桿端位移引起的桿端力和非結點荷載引起的單元固端力疊加，得到當前單元的全部桿端力。%------------------------------------------------------------------------------------------fprintf(FP2,'單元(%d)桿端力彎矩\n',i);fprintf(FP2,'左端%3d結點%7.2f\n',LNODS(i,1),FE(1));fprintf(FP2,'右端%3d結點%7.2f\n',LNODS(i,2),FE(2));說明：輸出單元的結點號，對應的桿端彎矩。end%對單元循環(huán)結束。%----------------------------------------------------------------------------------------------------return程序應用舉例【例1】利用連續(xù)梁程序計算圖2-9所示的連續(xù)梁結構。35353563圖2-9%-----------------------------------------------------------------------------------------------------輸入數據文件為：34112122334212314%-----------------------------------------------------------------------------------------------------說明：第一行：讀入程序控制信息NELEM =fscanf(FP1,'%d',1);%單元個數（單元編碼總數）NPOIN =fscanf(FP1,'%d',1);%結點個數（結點編碼總數）NVFIX =fscanf(FP1,'%d',1); ％約束個數（零位移總數）NFPOIN=fscanf(FP1,'%d',1); ％結點荷載個數（作用在結點上集中力偶總數）NFPRES=fscanf(FP1,'%d',1); ％非結點荷載數（作用在單元上分布荷載總數）YOUNG=fscanf(FP1,'%f',1); ％彈性模量第二、三四行：讀入單元連接信息：LNODS=fscanf(FP1,'%f',[3,NELEM])';%單元號，左結點號，右結點號，慣性矩第五、六、七、八行：讀入結點坐標COORD=fscanf(FP1,'%f',[NPOIN])';%結點坐標值第九行：讀入結點荷載信息FPOIN=fscanf(FP1,'%f',[2,NFPOIN])';%結點號，結點荷載數值（以順時針為正）第十、十一行：讀入非結點荷載信息FPRES=fscanf(FP1,'%f',[3,NFPRES])';%荷載號，單元號，荷載類型，荷載大小%荷載類型：1為均布荷載，2為集中力第十二行：讀入零位移信息FIXED=fscanf(FP1,'%f',NVFIX)';%零位移約束所在的結點位移編碼%-----------------------------------------------------------------------------------------------------輸出數據文件結構初始數據單元總數=3結點總數=4約束總數=1結點荷載個數=1非結點荷載個數=2彈性模量=1單元定義數據：單元號左結點號右結點號慣性矩112+000223+000334+000結點定義數據：結點號x坐標1234結點荷載數據：結點號M力偶12非結點荷載數據：荷載號單元號荷載類型荷載大小1123零位移約束數據：位移編號位移編號4計算結果輸出結構結點位移結點號轉角1+0012+0013+0014輸出單元桿端彎矩（順時針為正）單元(1)桿端力彎矩左端1結點右端2結點