有限單元法作業(yè)非線性分析+程序

1、幾何非線性大作業(yè)荷載增量法和弧長法程序設計系（所）：建筑工程系學 號：1432055姓 名：焦 聯(lián) 洪培養(yǎng)層次：專業(yè)碩士指導老師：吳 明 兒 2015年6月19日 一、 幾何非線性大作業(yè)（ Newton-Raphson法）用荷載增量法（Newton-Raphson法）編寫幾何非線性程序：（1）用平面梁單元，可分析平面桿系（2）算例：懸臂端作用彎矩。懸臂梁最終變形形成周長為懸臂梁長度的圓。1.1 Newton-Raphson算法基本思想圖1.1 Newton-Raphson算法基本思想1.2 懸臂梁參數(shù)基本參數(shù)：L=2m, D=0.03m, A=7.069E-4m2, I=3.976E-08m4

2、 ,E=2.0E11N/m2圖1.2 懸臂梁單元信息將懸臂梁分成10個單元，如圖1.2所示2.1 MATLAB輸入信息 材料信息 單元信息 約束信息（0為約束，1為放松） 荷載信息（FX,FY,M） 節(jié)點信息2.2 求解過程梁彎成圓形：理論彎矩M=EIY"=24981.944N.m ,直徑為0.642m運用ABAQUS和MATLAB進行求解對比：圖1.3 加載圖 圖1.4 ABAQUS變形圖圖1.5 MATLAB變形曲線ABAQUS和MATLAB變形對比，最終在理論荷載作用下都彎成了一個圓，其直徑為0.64716m，與理論值相對比值為：（0.64716-0.642）/0.642=0.

3、00804.非常接近。2.3 加載點荷載位移曲線圖1.5 加載點Y方向的荷載位移曲線加載點的最大豎向位移分別為1.4525m和1.45246m，相對比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，說明MATLAB的計算結(jié)果很好。二、 幾何非線性大作業(yè)（ 弧長法）用弧長法編寫幾何非線性程序，分析荷載位移全過程曲線：1) 用平面梁單元，可分析平面桿系結(jié)構(gòu)2) 算例(1)受集中荷載的拱：考察拱的矢跨比、荷載位置對荷載位移曲線的影響。(2)其他有復雜平衡路徑的結(jié)構(gòu)3) 將結(jié)果與相關(guān)文獻進行對比1.1 弧長法基本思想圖2.1 弧長法基本思想1.2 拱基本參數(shù)拱參數(shù)

4、：L=100m, A=0.32m2 ,I=1m4 ,E=1.0e7N/m2，F(xiàn)=-5000N，拱曲線 y=5×sin(3.1415926*x/L)將拱結(jié)構(gòu)分成25個單元，如圖2所示圖2.2 拱單元信息2.1 MATLAB輸入信息材料信息 單元信息約束信息（0為約束，1為放松） 荷載信息（FX,FY,M）節(jié)點信息2.2 運用ANSYS和MATLAB進行求解對比（兩端鉸接）ANSYS中模型:圖2.3 ANSYS模型 圖2.4 MATLAB和ANSYS變形圖2.3 加載點荷載位移曲線圖2.5 加載點荷載位移曲線ANSYS求得的極限承載力3042.53,對應位移3.00142MATLAB求得

5、的極限承載力3043.8, 對應位移3.0768相對誤差分別為0.0417%,2.45%，模擬效果比較好。2.4 拱的矢跨比a對拱荷載位移曲線的影響不同矢跨比（1/20，3/40，1/10，3/20）下加載點的荷載位移曲線1）MATLAB中計算拱的矢跨比a對拱荷載位移曲線的影響 圖2.6 荷載位移曲線圖2.7 荷載位移曲線表1 各矢跨比下拱結(jié)構(gòu)的極限荷載 參數(shù)矢高極值點F(N)位移（m）最低點F(N)位移（m）5mm3043.83.07681765.27.08167.5mm7623.34.0335-595.8211.2110mm149745.4026-6408.114.88620mm39791

6、9.4831-6304930.513從表中可以初步得出：在一定隨著矢跨比的增加，拱仍然呈現(xiàn)跳躍失穩(wěn)的形式，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力呈現(xiàn)出反向，相當于有一個拉力在阻止拱結(jié)構(gòu)發(fā)生跳躍失穩(wěn)，矢跨比越大，拱越不容易發(fā)生跳躍失穩(wěn)。當拱的矢跨比超過一定范圍后，拱將發(fā)生復雜的不同于跳躍失穩(wěn)的失穩(wěn)形式。2）MATLAB與ANSYS計算結(jié)果對比圖2.8 ANSYS和MATLAB對比荷載位移曲線表2 各矢跨比下拱結(jié)構(gòu)的極限荷載對比 參數(shù)矢高F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）5mm3043.83.07683042.533.001420.042.457.5

7、mm7623.34.03357624.913.96303-0.021.7510mm149745.402614974.35.31570.001.6120mm397919.483139695.79.599550.24-1.23從圖中可以看出：矢跨比在一定范圍內(nèi)，MATLAB與ANSYS計算的荷載位移曲線非常吻合，驗證了MATLAB程序的可行性。當矢跨比為0.15時，ANSYS中將跟蹤不到失穩(wěn)后復雜的平衡路徑。從表中可以得出：MATLAB與ANSYS計算中拱的極限荷載和極限荷載時所對應的位移非常接近，加載點均為頂點26。具體為：矢高5mm，荷載誤差為0.04，位移誤差為2.45；矢高7.5mm，荷載

8、誤差為-0.02，位移誤差為1.75；矢高10mm，荷載誤差為0，位移誤差為-1.61；矢高20mm，荷載誤差為0.24，位移誤差為-1.23。實際誤差相差很小，在工程允許的范圍內(nèi)是可以接受的。2.5 荷載位置對拱荷載位移曲線的影響圖2.9 ANSYS模型簡圖1）MATLAB中計算荷載位置對拱荷載位移曲線的影響圖2.10 各加載點荷載位移曲線表3 改變加載點拱結(jié)構(gòu)的極限荷載 參數(shù)加載點極值點F(N)位移（m）最低點F(N)位移（m）263043.83.07681765.27.08161635703.18912258.86.1161147282.883220.54.79594143171.282

9、6105691.7829 誤差=100*（MATLAB-ANSYS）/ANSYS從表中可以初步得出：隨著加載點位置越靠近支座處，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力也大幅度提高。當加載點位置靠近支座時，拱的承載力增加幅度最大，拱的穩(wěn)定性很強，不容易發(fā)生失穩(wěn)。2）MATLAB與ANSYS計算結(jié)果對比圖2.11 ANSYS和MATLAB對比荷載位移曲線表4 各加載點拱結(jié)構(gòu)的極限荷載 參數(shù)矢高F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）263043.83.07683042.533.001420.042.451635703.18913569.733.248650

10、.01-1.871147282.884728.712.91862-0.02-1.344143171.282614324.81.29764-0.05-1.17誤差=100*（MATLAB-ANSYS）/ANSYS從圖中可以看出：MATLAB與ANSYS計算的荷載位移曲線非常吻合，驗證了MATLAB程序的可行性。從表中可以得出：MATLAB與ANSYS計算中拱的極限荷載和極限荷載時所對應的位移非常接近。具體為：26加載點，荷載誤差為0.04，位移誤差為2.45；16加載點，荷載誤差為0.01，位移誤差為-1.87；11加載點，荷載誤差為-0.02，位移誤差為-1.34；4加載點，荷載誤差為-0.0

11、5，位移誤差為-1.17。實際誤差相差很小，在工程允許的范圍內(nèi)是可以接受的。2.6 兩端鉸接和固接對拱荷載位移曲線的影響矢高為5mm時，拱兩端為固接和鉸接時的荷載位移曲線如下：圖2.12 ANSYS和MATLAB固接和鉸接的荷載位移曲線從圖中可以看出：拱的邊界條件對其的失穩(wěn)形式有很大影響。兩端固接拱的穩(wěn)定性明顯優(yōu)于兩端鉸接拱，承載能力也大幅度提高。固接拱不會發(fā)生跳躍失穩(wěn)的形式，剛度在初始階段會減小，待到達一定程度后剛度又會增加。而兩端鉸接拱會發(fā)生跳躍失穩(wěn)的形式。2.7 參數(shù)m對拱失穩(wěn)形式的影響文獻中給出:m是一個由拱截面在豎向平面內(nèi)的回轉(zhuǎn)半徑r 和拱的初始矢高h無確定的無量綱參數(shù)。當m>

12、=1 時，扁拱不會出現(xiàn)跳躍屈曲, 當0<m<1時，扁拱有可能發(fā)生跳躍屈曲，而影響扁拱是否發(fā)生跳躍屈曲的主要因素是m值和荷載參數(shù)。 2.13 不同m值加載點的荷載位移曲線2.14 不同m值變形后拱曲線從MATLAB的計算結(jié)果中可以驗證：不同的m系數(shù)對應拱不同的失穩(wěn)形式。當m>=1 時，扁拱不會出現(xiàn)跳躍屈曲，當0<m<1時，扁拱有可能發(fā)生跳躍屈曲。但拱的最終變形圖非常接近，只是此時拱的失穩(wěn)形式變了。 2.8 具有復雜失穩(wěn)形式的拱鉸支圓拱該結(jié)構(gòu)及其幾何參數(shù)、物理性質(zhì)均示于圖4a 中。中心受集中荷載。這個結(jié)構(gòu)是研究分歧問題的經(jīng)典題目。將半跨結(jié)構(gòu)劃分為8個單元, 得到圖4b

13、的基本路徑和分歧路徑, 并與JChreseielewski 和Rsehmiot的結(jié)果進行了比較。本文對此結(jié)構(gòu)也進行了缺陷分析。拱的基本參數(shù)：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。文獻中的計算結(jié)果。采用MATLAB和ANSYS對其進行求解，得到其荷載位移曲線：圖2.15 ABAQUS模型圖2.16 ABAQUS變形圖圖2.17 ANSYS、MATLAB、ABAQUS加載點荷載位移曲線從MATLAB、ANSYS、ABAQUS計算的荷載位移曲線可以看出：兩者的荷載位移曲線基本吻合。MATLAB的計算結(jié)果可以看出在后期其承載能力會有較大提高

14、，與文獻中的荷載位移曲線趨勢相同，所以驗證出程序的可靠性。ABAQUS不能模擬出后續(xù)段曲線也許是單元劃分過少。圖2.18 MATLAB加載點荷載位移曲線 第二次極值點會超過第一次極值點所對應的荷載，與文獻一致，荷載點也接近。加入初始缺陷：L/1000, L/2000初始缺陷后觀察加載點的荷載位移曲線變化趨勢。圖2.19 ANSYS加入初始缺陷后的加載點荷載位移曲線2.20 初始缺陷為0.0001時的荷載位移曲線加入初始缺陷后，拱的極限承載能力明顯降低。且失穩(wěn)形式也發(fā)生了變化，初始缺陷的大小對其荷載位移曲線有明顯影響。所以在工程設計中應考慮結(jié)構(gòu)或構(gòu)件的初始缺陷，使結(jié)構(gòu)的設計更加合理，安全。為了研

15、究初始缺陷對拱失穩(wěn)路徑的影響，應用ABAQUS和ANSYS以及MATLAB中加水平力模擬拱結(jié)構(gòu)初始缺陷下的荷載位移曲線。為了探究ABAQUS和ANSYS計算結(jié)果，取其前三階模態(tài)進行對比分析： 2.21 一階屈曲模態(tài) ABAQUS和MATLAB中的一階屈曲系數(shù)為0.55884和0.564512，對應的屈曲荷載為74325.72N 和75080.096N。2.22 二階屈曲模態(tài) ABAQUS和MATLAB中的二階屈曲系數(shù)為1.2259和1.253，對應的屈曲163044.7N 和166649N。2.23 三階屈曲模態(tài)ABAQUS和MATLAB中的三階屈曲系數(shù)為2.166和2.255，對應的屈曲29

16、9915N 和288078N。從屈曲模態(tài)中可以看出，兩種軟件的前二階模態(tài)趨勢吻合，屈曲系數(shù)和極限荷載也是吻合的較好。第三階模態(tài)出現(xiàn)不一樣的變形趨勢，屈曲系數(shù)和極限荷載也是也相差比較大，但計算時只引入一階屈曲模態(tài)。圖2.24 ANSYS、ABAQUS、MATLAB加載點荷載位移曲線從圖中可以看出：ANSYS對缺陷的模擬效果比較好，與文獻中的比較接近ABAQUS模擬的極限荷載稍低于ANSYS，而MATLAB模擬的極限荷載遠低于ANSYS,但是曲線最終都在位移為300多mm時交于一點。還是有一定規(guī)律性。圖2.25 ANSYS和ABAQUS引入初始缺陷加載點荷載位移曲線 始缺陷并一定都會降低承載力，也

17、會有對結(jié)構(gòu)的承載能力有益的初始缺陷。ANSYS的計算結(jié)果可以看出，當初始缺陷為1/2000和-1/2000時，其承載能力不變。由于為對稱結(jié)構(gòu)，所以缺陷的正負影響不大。圖2.26 ANSYS和ABAQUS引入初始缺陷加載點荷載位移曲線表6 對比ANNSYS和ABAQUS的極限荷載值和其對應位移值 參數(shù)矢高F(N)ANS位移（m）F(N)ABA位移（m）誤差（%）誤差（%）1/100058444.768.979955795.62879.93184.53261-15.87691/200060579.970.138457924.95865.45424.382556.67851誤差=100*（ABAQU

18、S-ANSYS）/ABAQUS表中可以得出：ABAQUS求解出的機線荷載小于ANSYS，單對應的位移大于ANSYS對應的值。這可能與單元劃分的個數(shù)，求解精度有關(guān)，但在工程允許的范圍內(nèi)，還是可以接受的。ABAQUS中迭代步的跳躍很快，位移增加速度很快，其路徑不是很準確，可能是由于其單元劃分比較少。體會：1）注意坐標系的轉(zhuǎn)化和力、位移更新時所對應的狀態(tài)（C1-C2）2) 拱是否發(fā)生跳躍失穩(wěn)與矢跨比、矢高與截面旋轉(zhuǎn)半徑有關(guān)。矢跨比太大不會發(fā)生跳躍失穩(wěn)；m大于1時不能發(fā)生跳躍失穩(wěn)。3）有些拱在ANSYS中不能得出下降段，可見ANSYS中對拱的跨度矢高、截面參數(shù)的比值有一定要求。內(nèi)部計算和程序中有一些差

19、別。附錄子程序一：剛度組裝矩陣function K_G=Assemble(K_Element,Element,N_Node)K_G=zeros(N_Node*3,N_Node*3);for i=1:2 n1=Element(1,i); for j=1:2 n2=Element(1,j); K_G(3*n1-2:3*n1,3*n2-2:3*n2)=K_Element(3*i-2:3*i,3*j-2:3*j); endendend子程序二：整體坐標剛度矩陣function K=beam2d_stiffness530(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F

20、1(1,3);M2=Ele_F1(1,6);T=cs(1,1),cs(1,2),0,0,0,0; -cs(1,2),cs(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs(1,1),cs(1,2),0; 0,0,0,-cs(1,2),cs(1,1),0; 0,0,0,0,0,1;KE= E*A/L,0,0,-E*A/L,0,0; 0,12*E*I/L3,6*E*I/L2,0,-12*E*I/L3,6*E*I/L2;0,6*E*I/L2, 4*E*I/L, 0, -6*E*I/L2, 2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L3,-6*

21、E*I/L2,0,12*E*I/L3,-6*E*I/L2;0,6*E*I/L2,2*E*I/L,0,-6*E*I/L2,4*E*I/L;KG= F/L,0,-M1/L,-F/L,0,-M2/L; 0,12*F*I/A/L3+6*F/5/L, 6*F*I/A/L2+F/10, 0,-(12*F*I/A/L3+6*F/5/L), 6*F*I/A/L2+F/10;-M1/L, 6*F*I/A/L2+F/10, 4*F*I/A/L+2*F*L/15, M1/L, -(6*F*I/A/L2+F/10), 2*F*I/A/L-F*L/30;-F/L,0, M1/L,F/L,0,M2/L;0,-12*F*I

22、/A/L3-6*F/5/L,-6*F*I/A/L2-F/10,0,12*F*I/A/L3+6*F/5/L,-6*F*I/A/L2-F/10;-M2/L, 6*F*I/A/L2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L2+F/10), 4*F*I/A/L+2*F*L/15;K=T'*(KE+KG)*T;end子程序三：局部坐標剛度矩陣function K=beam2d_stiffness520(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F1(1,3);M2=Ele_F1(1,6);KE= E*A/L,0,0,

23、-E*A/L,0,0; 0,12*E*I/L3,6*E*I/L2,0,-12*E*I/L3,6*E*I/L2;0,6*E*I/L2, 4*E*I/L, 0, -6*E*I/L2, 2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L3,-6*E*I/L2,0,12*E*I/L3,-6*E*I/L2;0,6*E*I/L2,2*E*I/L,0,-6*E*I/L2,4*E*I/L;KG= F/L,0,-M1/L,-F/L,0,-M2/L; 0,12*F*I/A/L3+6*F/5/L, 6*F*I/A/L2+F/10, 0,-(12*F*I/A/L3+6*F/5/L), 6

24、*F*I/A/L2+F/10;-M1/L, 6*F*I/A/L2+F/10, 4*F*I/A/L+2*F*L/15, M1/L, -(6*F*I/A/L2+F/10), 2*F*I/A/L-F*L/30;-F/L,0, M1/L,F/L,0,M2/L;0,-12*F*I/A/L3-6*F/5/L,-6*F*I/A/L2-F/10,0,12*F*I/A/L3+6*F/5/L,-6*F*I/A/L2-F/10;-M2/L, 6*F*I/A/L2+F/10, 2*F*I/A/L-F*L/30, M2/L, -(6*F*I/A/L2+F/10), 4*F*I/A/L+2*F*L/15;K=(KE+KG

25、);End主程序一：Newton-Raphson法clcclearNode=xlsread('Data.xls','Node');Element=xlsread('Data.xls','Element');Boundary=xlsread('Data.xls','Boundary');Section=xlsread('Data.xls','Section');Force=xlsread('Data.xls','Force');%讀取相關(guān)

26、數(shù)據(jù)Allstep=1000; %迭代步數(shù)N_Node=size(Node,1); %節(jié)點個數(shù) N_Element=size(Element,1); %單元個數(shù)N_Force=size(Force,1); %節(jié)點力個數(shù)N_Boundary=size(Boundary,1); %約束節(jié)點個數(shù)Displacement=zeros(N_Node,3); %位移向量（a0）%初始位移及轉(zhuǎn)角為0All_Disp=zeros(N_Node,3); %初始位移和轉(zhuǎn)角為零（迭代后的節(jié)點位移）All_F=zeros(N_Node*3,1); %初始荷載向量為零 （存放節(jié)點荷載向量）Internal_F=zero

27、s(N_Node*3,1); %節(jié)點內(nèi)力向量ForceMatrix=zeros(N_Node*3,1); %總荷載向量C=zeros(N_Element,2);L=zeros(N_Element,1);for i=1:N_Force ForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4)' %把節(jié)點荷載向量讀入一個矩陣中,形成列向量=總的自由度個數(shù)enddel=;for i=1:N_Boundary; if Boundary(i,2)=0; m=3*Boundary(i,1)-2; del=del,m; %受約束節(jié)點位移為0時所對

28、應的指標數(shù)值1 end if Boundary(i,3)=0; m=3*Boundary(i,1)-1; del=del,m; %受約束節(jié)點位移為0時所對應的指標數(shù)值2 end if Boundary(i,4)=0; m=3*Boundary(i,1); del=del,m; %受約束節(jié)點位移為0時所對應的指標數(shù)值3 endend%求出約束節(jié)點的標號，便于剛度、荷載矩陣歸0Update_Node=Node+Displacement(:,1:2); %更新后的節(jié)點坐標向量（x，y坐標）Ele_F=zeros(N_Element,6); %單元節(jié)點荷載向量ELEDISP=zeros(N_Eleme

29、nt,6); %單元節(jié)點位移向量Ele_F1=zeros(N_Element,6); %單元節(jié)點荷載剛度矩陣中向量Ele1_F=zeros(1,6);ELEDISP1=zeros(1,6); qq(1,1)=0; pp(1,1)=0;for n=0:Allstep-1 n=n+1K_Global=zeros(N_Node*3,N_Node*3); %總剛矩陣 for i=1:N_Element N1=Element(i,1); %i節(jié)點編號 N2=Element(i,2); %j節(jié)點編號 N_Section=Element(i,3); %截面的形狀控制參數(shù) C(i,:)=Update_Node

30、(N2,:)-Update_Node(N1,:); %a0下坐標向量增量 L(i)=norm(C(i,:); %變形后的長度 cs=C(i,:)/L(i); %桿件的cos和sin值 E=Section(N_Section,1); A=Section(N_Section,2); I=Section(N_Section,3); K_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:); K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成總剛K0 end%整體剛度矩陣 D

31、elta_Force=ForceMatrix/Allstep;%初始荷載向量（Q0） Equation=K_Global,Delta_Force; %方程組 Disp_Transefer=ones(N_Node*3,1); %建立調(diào)節(jié)位移矩陣的位移向量 Disp_Transefer(del,:)=0; %將約束節(jié)點的位移值定為0,其他的定位1 Equation(del,:)=;%把方程中約束節(jié)點所對應的行和列去掉 Equation(:,del)=; %引入約束條件修改方程組 n1=size(Equation,2); % 方程組中列數(shù) dis1=(Equation(:,1:n1-1)Equati

32、on(:,n1); % 剛度矩陣求逆后乘以荷載向量, Equation(:,n1)荷載向量，得到節(jié)點的位移（da0） %求解方程組 zz=1; %識別約束 for i=1:N_Node*3; if Disp_Transefer(i,1)=1; Disp_Transefer(i,1)=dis1(zz,1); %總的位移向量 zz=zz+1; end end for i=1:N_Node Displacement(i,:)=Disp_Transefer(i*3-2:i*3,1); % 位移增量，形成n*3的位移向量（da0） end All_Disp=All_Disp+Displacement %

33、位移向量更新得到a1（總的位移增量） All_F=All_F+Delta_Force; %外荷載向量更新Q1 Internal_F1=zeros(N_Node*3,1); %節(jié)點內(nèi)力向量 Update_Node1=Update_Node; %C1狀態(tài) Update_Node=Node+All_Disp(:,1:2); %C2狀態(tài)坐標位置更新（改）（迭代后的位置 for i=1:N_Element F_Global=zeros(N_Node*3,1); %全局的荷載向量 for j=1:2 ELEDISP1(j*3-2:j*3)=Displacement(Element(i,j),:); %整體

34、 取出當前單元節(jié)點位移向量 end N1=Element(i,1); %i節(jié)點編號 N2=Element(i,2); %j節(jié)點編號 C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐標向量增量 L1(i)=norm(C1(i,:); %變形后的長度 cs1=C1(i,:)/L1(i); %桿件的cos和sin值 T1=cs1(1,1),cs1(1,2),0,0,0,0; -cs1(1,2),cs1(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs1(1,1),cs1(1,2),0; 0,0,0,-cs1(1,2),cs

35、1(1,1),0; 0,0,0,0,0,1; ELEDISP(i,:)=T1* ELEDISP1(:); X1=L1(i)+ ELEDISP(i,4)- ELEDISP(i,1); Z1= ELEDISP(i,5)-ELEDISP(i,2); LN=(X12+Z12)0.5; sin1=Z1/LN; cos1=X1/LN; Citaloca=atan2(sin1,cos1); Ub=LN-L1(i); THRA=ELEDISP(i,3)-Citaloca; THRB=ELEDISP(i,6)-Citaloca; ELEDISP(i,:)=0,0,THRA,Ub,0,THRB; K_Elemen

36、t=beam2d_stiffness520(E,A,I,L1(i),cs1,Ele_F1(i,:); %L（i）為a0對應下的 Ele_F(i,:)=K_Element*ELEDISP(i,:)' %局部坐標系下荷載（Q0）作用下的節(jié)點力 Ele_F1(i,:)= Ele_F1(i,:)+ Ele_F(i,:); C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐標向量增量 L2(i)=norm(C2(i,:); %變形后的長度 cs2=C2(i,:)/L2(i); %桿件的cos和sin值 T2=cs2(1,1),cs2(1,2),0

37、,0,0,0; -cs2(1,2),cs2(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs2(1,1),cs2(1,2),0; 0,0,0,-cs2(1,2),cs2(1,1),0; 0,0,0,0,0,1; Ele1_F(:)=T2'*Ele_F1(i,:)'%行向量 N1=Element(i,1); N2=Element(i,2); F_Global(3*N1-2:3*N1,1)=Ele1_F(1:3); %i節(jié)點荷載 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j節(jié)點荷載 Internal_F1=Internal_F

38、1+F_Global; %a1對應下的P1 end Val=Internal_F1-All_F %求出上次迭代完成時的殘余應力Q1-P1 Correc_dis=zeros(N_Node,3); %構(gòu)造新的節(jié)點位移向量，每次更新 Correc_dis1=zeros(N_Node,3); %構(gòu)造新的節(jié)點位移向量，疊加位移增量 N_dis=size(dis1,1); %未受約束的節(jié)點位移數(shù)目,不為零的節(jié)點位移數(shù)目 dis3=zeros(N_dis,1); %構(gòu)造一個向量 k=0;%修改位移矩陣形式 while norm(Val)>1e-7 & k<500 %在一個荷載增量步下進行

39、的對此迭代 k=k+1; K_Global=zeros(N_Node*3,N_Node*3); for i=1:N_Element N1=Element(i,1); N2=Element(i,2); N_Section=Element(i,3); C(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a1下坐標向量增量 L(i)=norm(C(i,:); %變形后的長度 cs=C(i,:)/L(i); E=Section(N_Section,1); A=Section(N_Section,2); I=Section(N_Section,3); K_Eleme

40、nt=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:); K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node); %形成總剛k1 end Equation=K_Global,Val; %在殘余應力下的位移求解 Disp_Transefer=ones(N_Node*3,1); Disp_Transefer(del,:)=0; Equation(del,:)=; Equation(:,del)=; n1=size(Equation,2); Dis2=-(Equation(:,1:n1-1)Equ

41、ation(:,n1) %荷位移增量da1 zz=1; for i=1:N_Node*3; if Disp_Transefer(i,1)=1; Disp_Transefer(i,1)=Dis2(zz,1); zz=zz+1; end end for i=1:N_Node Correc_dis(i,:)=Disp_Transefer(i*3-2:i*3,1); end Correc_dis1= Correc_dis1+ Correc_dis; Internal_F1=zeros(N_Node*3,1); %節(jié)點內(nèi)力向量 Update_Node1=Update_Node; Update_Node=

42、Node+All_Disp(:,1:2)+Correc_dis1(:,1:2);%為了求切線剛度矩陣（改）a2下 if abs(Update_Node(11,2)<=1e-4&&abs(Update_Node(11,1)<=1e-4 break end for i=1:N_Element F_Global=zeros(N_Node*3,1); for j=1:2 ELEDISP1(j*3-2:j*3)=Correc_dis(Element(i,j),:); %取出當前單元節(jié)點位移向量 end N1=Element(i,1); %i節(jié)點編號 N2=Element(i,

43、2); %j節(jié)點編號 C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐標向量增量 L1(i)=norm(C1(i,:); %變形后的長度 cs1=C1(i,:)/L1(i); %桿件的cos和sin值 T1=cs1(1,1),cs1(1,2),0,0,0,0; -cs1(1,2),cs1(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs1(1,1),cs1(1,2),0; 0,0,0,-cs1(1,2),cs1(1,1),0; 0,0,0,0,0,1; ELEDISP(i,:)=T1* ELEDISP1(:); X1

44、=L1(i)+ ELEDISP(i,4)- ELEDISP(i,1); Z1= ELEDISP(i,5)-ELEDISP(i,2); LN=(X12+Z12)0.5; sin1=Z1/LN; cos1=X1/LN; Citaloca=atan2(sin1,cos1); Ub=LN-L1(i); THRA=ELEDISP(i,3)- Citaloca; THRB=ELEDISP(i,6)- Citaloca; ELEDISP(i,:)=0,0,THRA,Ub,0,THRB; K_Element=beam2d_stiffness520(E,A,I,L1(i),cs1,Ele_F1(i,:);% L

45、（i）為a1對應下的 Ele_F(i,:)=K_Element*ELEDISP(i,:)' Ele_F1(i,:)= Ele_F1(i,:)+ Ele_F(i,:); C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:); %a0下坐標向量增量 L2(i)=norm(C2(i,:); %變形后的長度 cs2=C2(i,:)/L2(i); %桿件的cos和sin值 T2=cs2(1,1),cs2(1,2),0,0,0,0; -cs2(1,2),cs2(1,1),0,0,0,0; 0,0,1,0,0,0; 0,0,0,cs2(1,1),cs2(1,2),0

46、; 0,0,0,-cs2(1,2),cs2(1,1),0; 0,0,0,0,0,1; Ele1_F(:)=T2'*Ele_F1(i,:)'%行向量 N1=Element(i,1); N2=Element(i,2); F_Global(3*N1-2:3*N1,1)=Ele1_F(1:3); %i節(jié)點荷載 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j節(jié)點荷載 Internal_F1=Internal_F1+F_Global; %a1對應下的P1 end Val=Internal_F1-All_F; %荷載增量Q1-P2 Val(del,:)=0;

47、end All_Disp=All_Disp+Correc_dis1 qq(n+1,1)=All_F(N_Node*3,1); pp(n+1,1)=All_Disp(21,2);endplot(1.4525,qq,'g')text(1.3,1.5*104,'x=1.4525')hold onplot(pp,qq,'r')title('懸臂梁加載點荷載位移曲線');xlabel('位移（m）');ylabel('荷載（N）');legend('點11荷載位移曲線');grid主程序二：弧長法clcclearNode=xlsread('Data111.xls','Node');Element=xlsread('Data111.xls','Element');Boundary=xlsread('Data111.xls','Boundary&#