幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計

幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計系（所）：建筑工程系學(xué)號：1432055姓名：焦聯(lián)洪培養(yǎng)層次：專業(yè)碩士指導(dǎo)老師：吳明兒2015年6月19日幾何非線性大作業(yè)（Newton-Raphson法）用荷載增量法（Newton-Raphson法）編寫幾何非線性程序：（1）用平面梁單元，可分析平面桿系（2）算例：懸臂端作用彎矩。懸臂梁最終變形形成周長為懸臂梁長度的圓。1.1Newton-Raphson算法基本思想圖1.1Newton-Raphson算法基本思想1.2懸臂梁參數(shù)基本參數(shù)：L=2m,D=0.03m,A=7.069E-4m2,I=3.976E-08m4,E=2.0E11N/m2圖1.2懸臂梁單元信息將懸臂梁分成10個單元，如圖1.2所示2.1MATLAB輸入信息材料信息單元信息約束信息（0為約束，1為放松）荷載信息（FX,FY,M）節(jié)點信息2.2求解過程梁彎成圓形：理論彎矩M=EIY"=24981.944N.m,直徑為0.642m運(yùn)用ABAQUS和MATLAB進(jìn)行求解對比：圖1.3加載圖圖1.4ABAQUS變形圖圖1.5MATLAB變形曲線ABAQUS和MATLAB變形對比，最終在理論荷載作用下都彎成了一個圓，其直徑為0.64716m，與理論值相對比值為：（0.64716-0.642）/0.642=0.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)其他有復(fù)雜平衡路徑的結(jié)構(gòu)3)將結(jié)果與相關(guān)文獻(xiàn)進(jìn)行對比1.1弧長法基本思想圖2.1弧長法基本思想1.2拱基本參數(shù)拱參數(shù)：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.1MATLAB輸入信息材料信息單元信息約束信息（0為約束，1為放松）荷載信息（FX,FY,M）節(jié)點信息2.2運(yùn)用ANSYS和MATLAB進(jìn)行求解對比（兩端鉸接）ANSYS中模型:圖2.3ANSYS模型圖2.4MATLAB和ANSYS變形圖2.3加載點荷載位移曲線圖2.5加載點荷載位移曲線ANSYS求得的極限承載力3042.53,對應(yīng)位移3.00142

MATLAB求得的極限承載力3043.8,對應(yīng)位移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.88620mm397919.4831-6304930.513從表中可以初步得出：在一定隨著矢跨比的增加，拱仍然呈現(xiàn)跳躍失穩(wěn)的形式，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力呈現(xiàn)出反向，相當(dāng)于有一個拉力在阻止拱結(jié)構(gòu)發(fā)生跳躍失穩(wěn)，矢跨比越大，拱越不容易發(fā)生跳躍失穩(wěn)。當(dāng)拱的矢跨比超過一定范圍后，拱將發(fā)生復(fù)雜的不同于跳躍失穩(wěn)的失穩(wěn)形式。2）MATLAB與ANSYS計算結(jié)果對比圖2.8ANSYS和MATLAB對比荷載位移曲線表2各矢跨比下拱結(jié)構(gòu)的極限荷載對比參數(shù)矢高 F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）5mm3043.83.07683042.533.001420.042.457.5mm7623.34.03357624.913.96303-0.021.7510mm149745.402614974.35.31570.001.6120mm397919.483139695.79.599550.24-1.23從圖中可以看出：矢跨比在一定范圍內(nèi)，MATLAB與ANSYS計算的荷載位移曲線非常吻合，驗證了MATLAB程序的可行性。當(dāng)矢跨比為0.15時，ANSYS中將跟蹤不到失穩(wěn)后復(fù)雜的平衡路徑。從表中可以得出：MATLAB與ANSYS計算中拱的極限荷載和極限荷載時所對應(yīng)的位移非常接近，加載點均為頂點26。具體為：矢高5mm，荷載誤差為0.04，位移誤差為2.45；矢高7.5mm，荷載誤差為-0.02，位移誤差為1.75；矢高10mm，荷載誤差為0，位移誤差為-1.61；矢高20mm，荷載誤差為0.24，位移誤差為-1.23。實際誤差相差很小，在工程允許的范圍內(nèi)是可以接受的。2.5荷載位置對拱荷載位移曲線的影響圖2.9ANSYS模型簡圖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.2826105691.7829誤差=100*（MATLAB-ANSYS）/ANSYS從表中可以初步得出：隨著加載點位置越靠近支座處，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力也大幅度提高。當(dāng)加載點位置靠近支座時，拱的承載力增加幅度最大，拱的穩(wěn)定性很強(qiáng)，不容易發(fā)生失穩(wěn)。2）MATLAB與ANSYS計算結(jié)果對比圖2.11ANSYS和MATLAB對比荷載位移曲線表4各加載點拱結(jié)構(gòu)的極限荷載參數(shù)矢高 F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）263043.83.07683042.533.001420.042.451635703.18913569.733.248650.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計算中拱的極限荷載和極限荷載時所對應(yīng)的位移非常接近。具體為：26加載點，荷載誤差為0.04，位移誤差為2.45；16加載點，荷載誤差為0.01，位移誤差為-1.87；11加載點，荷載誤差為-0.02，位移誤差為-1.34；4加載點，荷載誤差為-0.05，位移誤差為-1.17。實際誤差相差很小，在工程允許的范圍內(nèi)是可以接受的。2.6兩端鉸接和固接對拱荷載位移曲線的影響矢高為5mm時，拱兩端為固接和鉸接時的荷載位移曲線如下：圖2.12ANSYS和MATLAB固接和鉸接的荷載位移曲線從圖中可以看出：拱的邊界條件對其的失穩(wěn)形式有很大影響。兩端固接拱的穩(wěn)定性明顯優(yōu)于兩端鉸接拱，承載能力也大幅度提高。固接拱不會發(fā)生跳躍失穩(wěn)的形式，剛度在初始階段會減小，待到達(dá)一定程度后剛度又會增加。而兩端鉸接拱會發(fā)生跳躍失穩(wěn)的形式。2.7參數(shù)m對拱失穩(wěn)形式的影響文獻(xiàn)中給出:m是一個由拱截面在豎向平面內(nèi)的回轉(zhuǎn)半徑r和拱的初始矢高h(yuǎn)無確定的無量綱參數(shù)。當(dāng)m>=1時，扁拱不會出現(xiàn)跳躍屈曲,當(dāng)0<m<1時，扁拱有可能發(fā)生跳躍屈曲，而影響扁拱是否發(fā)生跳躍屈曲的主要因素是m值和荷載參數(shù)。2.13不同m值加載點的荷載位移曲線2.14不同m值變形后拱曲線從MATLAB的計算結(jié)果中可以驗證：不同的m系數(shù)對應(yīng)拱不同的失穩(wěn)形式。當(dāng)m>=1時，扁拱不會出現(xiàn)跳躍屈曲，當(dāng)0<m<1時，扁拱有可能發(fā)生跳躍屈曲。但拱的最終變形圖非常接近，只是此時拱的失穩(wěn)形式變了。2.8具有復(fù)雜失穩(wěn)形式的拱鉸支圓拱該結(jié)構(gòu)及其幾何參數(shù)、物理性質(zhì)均示于圖4a中。中心受集中荷載。這個結(jié)構(gòu)是研究分歧問題的經(jīng)典題目。將半跨結(jié)構(gòu)劃分為8個單元,得到圖4b的基本路徑和分歧路徑,并與JChreseielewski和Rsehmiot的結(jié)果進(jìn)行了比較。本文對此結(jié)構(gòu)也進(jìn)行了缺陷分析。拱的基本參數(shù)：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。文獻(xiàn)中的計算結(jié)果。采用MATLAB和ANSYS對其進(jìn)行求解，得到其荷載位移曲線：圖2.15ABAQUS模型圖2.16ABAQUS變形圖圖2.17ANSYS、MATLAB、ABAQUS加載點荷載位移曲線從MATLAB、ANSYS、ABAQUS計算的荷載位移曲線可以看出：兩者的荷載位移曲線基本吻合。MATLAB的計算結(jié)果可以看出在后期其承載能力會有較大提高，與文獻(xiàn)中的荷載位移曲線趨勢相同，所以驗證出程序的可靠性。ABAQUS不能模擬出后續(xù)段曲線也許是單元劃分過少。圖2.18MATLAB加載點荷載位移曲線第二次極值點會超過第一次極值點所對應(yīng)的荷載，與文獻(xiàn)一致，荷載點也接近。加入初始缺陷：L/1000,L/2000初始缺陷后觀察加載點的荷載位移曲線變化趨勢。圖2.19ANSYS加入初始缺陷后的加載點荷載位移曲線2.20初始缺陷為0.0001時的荷載位移曲線加入初始缺陷后，拱的極限承載能力明顯降低。且失穩(wěn)形式也發(fā)生了變化，初始缺陷的大小對其荷載位移曲線有明顯影響。所以在工程設(shè)計中應(yīng)考慮結(jié)構(gòu)或構(gòu)件的初始缺陷，使結(jié)構(gòu)的設(shè)計更加合理，安全。為了研究初始缺陷對拱失穩(wěn)路徑的影響，應(yīng)用ABAQUS和ANSYS以及MATLAB中加水平力模擬拱結(jié)構(gòu)初始缺陷下的荷載位移曲線。為了探究ABAQUS和ANSYS計算結(jié)果，取其前三階模態(tài)進(jìn)行對比分析：2.21一階屈曲模態(tài)ABAQUS和MATLAB中的一階屈曲系數(shù)為0.55884和0.564512，對應(yīng)的屈曲荷載為74325.72N和75080.096N。2.22二階屈曲模態(tài)ABAQUS和MATLAB中的二階屈曲系數(shù)為1.2259和1.253，對應(yīng)的屈曲163044.7N和166649N。2.23三階屈曲模態(tài)ABAQUS和MATLAB中的三階屈曲系數(shù)為2.166和2.255，對應(yīng)的屈曲299915N和288078N。從屈曲模態(tài)中可以看出，兩種軟件的前二階模態(tài)趨勢吻合，屈曲系數(shù)和極限荷載也是吻合的較好。第三階模態(tài)出現(xiàn)不一樣的變形趨勢，屈曲系數(shù)和極限荷載也是也相差比較大，但計算時只引入一階屈曲模態(tài)。圖2.24ANSYS、ABAQUS、MATLAB加載點荷載位移曲線從圖中可以看出：ANSYS對缺陷的模擬效果比較好，與文獻(xiàn)中的比較接近ABAQUS模擬的極限荷載稍低于ANSYS，而MATLAB模擬的極限荷載遠(yuǎn)低于ANSYS,但是曲線最終都在位移為300多mm時交于一點。還是有一定規(guī)律性。圖2.25ANSYS和ABAQUS引入初始缺陷加載點荷載位移曲線始缺陷并一定都會降低承載力，也會有對結(jié)構(gòu)的承載能力有益的初始缺陷。ANSYS的計算結(jié)果可以看出，當(dāng)初始缺陷為1/2000和-1/2000時，其承載能力不變。由于為對稱結(jié)構(gòu)，所以缺陷的正負(fù)影響不大。圖2.26ANSYS和ABAQUS引入初始缺陷加載點荷載位移曲線表6對比ANNSYS和ABAQUS的極限荷載值和其對應(yīng)位移值參數(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*（ABAQUS-ANSYS）/ABAQUS表中可以得出：ABAQUS求解出的機(jī)線荷載小于ANSYS，單對應(yīng)的位移大于ANSYS對應(yīng)的值。這可能與單元劃分的個數(shù)，求解精度有關(guān)，但在工程允許的范圍內(nèi)，還是可以接受的。ABAQUS中迭代步的跳躍很快，位移增加速度很快，其路徑不是很準(zhǔn)確，可能是由于其單元劃分比較少。體會：1）注意坐標(biāo)系的轉(zhuǎn)化和力、位移更新時所對應(yīng)的狀態(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)部計算和程序中有一些差別。附錄子程序一：剛度組裝矩陣functionK_G=Assemble(K_Element,Element,N_Node)K_G=zeros(N_Node*3,N_Node*3);fori=1:2n1=Element(1,i);forj=1:2n2=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子程序二：整體坐標(biāo)剛度矩陣functionK=beam2d_stiffness530(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F1(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/L^3,6*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2;0,6*E*I/L^2,4*E*I/L,0,-6*E*I/L^2,2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2;0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L];KG=[F/L,0,-M1/L,-F/L,0,-M2/L;0,12*F*I/A/L^3+6*F/5/L,6*F*I/A/L^2+F/10,0,-(12*F*I/A/L^3+6*F/5/L),6*F*I/A/L^2+F/10;-M1/L,6*F*I/A/L^2+F/10,4*F*I/A/L+2*F*L/15,M1/L,-(6*F*I/A/L^2+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/L^3-6*F/5/L,-6*F*I/A/L^2-F/10,0,12*F*I/A/L^3+6*F/5/L,-6*F*I/A/L^2-F/10;-M2/L,6*F*I/A/L^2+F/10,2*F*I/A/L-F*L/30,M2/L,-(6*F*I/A/L^2+F/10),4*F*I/A/L+2*F*L/15];K=T'*(KE+KG)*T;end子程序三：局部坐標(biāo)剛度矩陣functionK=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,-E*A/L,0,0;0,12*E*I/L^3,6*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2;0,6*E*I/L^2,4*E*I/L,0,-6*E*I/L^2,2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2;0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L];KG=[F/L,0,-M1/L,-F/L,0,-M2/L;0,12*F*I/A/L^3+6*F/5/L,6*F*I/A/L^2+F/10,0,-(12*F*I/A/L^3+6*F/5/L),6*F*I/A/L^2+F/10;-M1/L,6*F*I/A/L^2+F/10,4*F*I/A/L+2*F*L/15,M1/L,-(6*F*I/A/L^2+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/L^3-6*F/5/L,-6*F*I/A/L^2-F/10,0,12*F*I/A/L^3+6*F/5/L,-6*F*I/A/L^2-F/10;-M2/L,6*F*I/A/L^2+F/10,2*F*I/A/L-F*L/30,M2/L,-(6*F*I/A/L^2+F/10),4*F*I/A/L+2*F*L/15];K=(KE+KG);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)數(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=zeros(N_Node*3,1);%節(jié)點內(nèi)力向量ForceMatrix=zeros(N_Node*3,1);%總荷載向量C=zeros(N_Element,2);L=zeros(N_Element,1);fori=1:N_ForceForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4)';%把節(jié)點荷載向量讀入一個矩陣中,形成列向量=總的自由度個數(shù)enddel=[];fori=1:N_Boundary;ifBoundary(i,2)==0;m=3*Boundary(i,1)-2;del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值1endifBoundary(i,3)==0;m=3*Boundary(i,1)-1;del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值2endifBoundary(i,4)==0;m=3*Boundary(i,1);del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值3endend%求出約束節(jié)點的標(biāo)號，便于剛度、荷載矩陣歸0Update_Node=Node+Displacement(:,1:2);%更新后的節(jié)點坐標(biāo)向量（x，y坐標(biāo)）Ele_F=zeros(N_Element,6);%單元節(jié)點荷載向量ELEDISP=zeros(N_Element,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;forn=0:Allstep-1n=n+1K_Global=zeros(N_Node*3,N_Node*3);%總剛矩陣fori=1:N_ElementN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號N_Section=Element(i,3);%截面的形狀控制參數(shù)C(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a0下坐標(biāo)向量增量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);%形成總剛K0end%整體剛度矩陣Delta_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,其他的定位1Equation(del,:)=[];%把方程中約束節(jié)點所對應(yīng)的行和列去掉Equation(:,del)=[];%引入約束條件修改方程組n1=size(Equation,2);%方程組中列數(shù)dis1=(Equation(:,1:n1-1))\Equation(:,n1);%剛度矩陣求逆后乘以荷載向量,Equation(:,n1)荷載向量，得到節(jié)點的位移（da0）%求解方程組zz=1;%識別約束fori=1:N_Node*3;ifDisp_Transefer(i,1)==1;Disp_Transefer(i,1)=dis1(zz,1);%總的位移向量zz=zz+1;endendfori=1:N_NodeDisplacement(i,:)=Disp_Transefer(i*3-2:i*3,1);%位移增量，形成n*3的位移向量（da0）endAll_Disp=All_Disp+Displacement%位移向量更新得到a1（總的位移增量）All_F=All_F+Delta_Force;%外荷載向量更新Q1Internal_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)坐標(biāo)位置更新（改）（迭代后的位置fori=1:N_ElementF_Global=zeros(N_Node*3,1);%全局的荷載向量forj=1:2ELEDISP1(j*3-2:j*3)=Displacement(Element(i,j),:);%整體取出當(dāng)前單元節(jié)點位移向量endN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:);%a0下坐標(biāo)向量增量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=L1(i)+ELEDISP(i,4)-ELEDISP(i,1);Z1=ELEDISP(i,5)-ELEDISP(i,2);LN=(X1^2+Z1^2)^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（i）為a0對應(yīng)下的Ele_F(i,:)=K_Element*ELEDISP(i,:)';%局部坐標(biāo)系下荷載（Q0）作用下的節(jié)點力Ele_F1(i,:)=Ele_F1(i,:)+Ele_F(i,:);C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a0下坐標(biāo)向量增量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;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對應(yīng)下的P1endVal=Internal_F1-All_F%求出上次迭代完成時的殘余應(yīng)力Q1-P1Correc_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;%修改位移矩陣形式whilenorm(Val)>1e-7&k<500%在一個荷載增量步下進(jìn)行的對此迭代k=k+1;K_Global=zeros(N_Node*3,N_Node*3);fori=1:N_ElementN1=Element(i,1);N2=Element(i,2);N_Section=Element(i,3);C(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a1下坐標(biāo)向量增量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_Element=beam2d_stiffness530(E,A,I,L(i),cs,Ele_F1(i,:));K_Global=K_Global+Assemble(K_Element,Element(i,:),N_Node);%形成總剛k1endEquation=[K_Global,Val];%在殘余應(yīng)力下的位移求解Disp_Transefer=ones(N_Node*3,1);Disp_Transefer(del,:)=0;Equation(del,:)=[];Equation(:,del)=[];n1=size(Equation,2);Dis2=-(Equation(:,1:n1-1))\Equation(:,n1)%荷位移增量da1zz=1;fori=1:N_Node*3;ifDisp_Transefer(i,1)==1;Disp_Transefer(i,1)=Dis2(zz,1);zz=zz+1;endendfori=1:N_NodeCorrec_dis(i,:)=Disp_Transefer(i*3-2:i*3,1);endCorrec_dis1=Correc_dis1+Correc_dis;Internal_F1=zeros(N_Node*3,1);%節(jié)點內(nèi)力向量Update_Node1=Update_Node;Update_Node=Node+All_Disp(:,1:2)+Correc_dis1(:,1:2);%為了求切線剛度矩陣（改）a2下ifabs(Update_Node(11,2))<=1e-4&&abs(Update_Node(11,1))<=1e-4breakendfori=1:N_ElementF_Global=zeros(N_Node*3,1);forj=1:2ELEDISP1(j*3-2:j*3)=Correc_dis(Element(i,j),:);%取出當(dāng)前單元節(jié)點位移向量endN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:);%a0下坐標(biāo)向量增量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=L1(i)+ELEDISP(i,4)-ELEDISP(i,1);Z1=ELEDISP(i,5)-ELEDISP(i,2);LN=(X1^2+Z1^2)^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（i）為a1對應(yīng)下的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下坐標(biāo)向量增量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;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對應(yīng)下的P1endVal=Internal_F1-All_F;%荷載增量Q1-P2Val(del,:)=0;endAll_Disp=All_Disp+Correc_dis1qq(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*10^4,'x=1.4525')holdonplot(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');Section=xlsread('Data111.xls','Section');Force=xlsread('Data111.xls','Force');%讀取相關(guān)數(shù)據(jù)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é)點荷載向量）ForceMatrix=zeros(N_Node*3,1);%總荷載向量C=zeros(N_Element,2);L=zeros(N_Element,1);fori=1:N_ForceForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4)';%把節(jié)點荷載向量讀入一個矩enddel=[];fori=1:N_Boundary;ifBoundary(i,2)==0;m=3*Boundary(i,1)-2;del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值1endifBoundary(i,3)==0;m=3*Boundary(i,1)-1;del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值2endifBoundary(i,4)==0;m=3*Boundary(i,1);del=[del,[m]];%受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值3endend%求出約束節(jié)點的標(biāo)號，便于剛度、荷載矩陣歸0Update_Node=Node+Displacement(:,1:2);%更新后的節(jié)點坐標(biāo)向量（x，y坐標(biāo)）Ele_F=zeros(N_Element,6);%單元節(jié)點荷載向量ELEDISP=zeros(N_Element,6);%單元節(jié)點位移向量Ele_F1=zeros(N_Element,6);%單元節(jié)點荷載剛度矩陣中向量Ele1_F=zeros(1,6);%臨時存放向量ELEDISP1=zeros(1,6);%臨時存放向量flag=1;n=0;DetaL=0.1;%初始迭代荷載系數(shù)chushidian=zeros(1,N_Node*3+1);chushidian(del)=[];%存放m-1點的位移+荷載系數(shù)向量Deta_Disp=zeros(N_Node,3);%weiyiDetaVH=zeros(N_Node*3,1);DetaVH(del)=[];qq(1,1)=0;pp(1,1)=0;whileflag==1&&n<20000n=n+1DetaV=zeros(N_Node*3,1);DetaV(del)=[];K_Global=zeros(N_Node*3,N_Node*3);%總剛矩陣fori=1:N_ElementN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號N_Section=Element(i,3);%截面的形狀控制參數(shù)C(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a0下坐標(biāo)向量增量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);%形成總剛K0K_Global1=K_Global;end%整體剛度矩陣%%%%%%%%%%%%%%%%%弧長法第一步%%%%%%%%%%ifn==1Equation=[K_Global1,ForceMatrix];%方程組Equation(del,:)=[];%把方程中約束節(jié)點所對應(yīng)的行和列去掉Equation(:,del)=[];%引入約束條件修改方程組n1=size(Equation,2);%方程組中列數(shù)dis=(Equation(:,1:n1-1))\Equation(:,n1);%剛度矩陣求逆后乘以荷載向量,Equation(:,n1)荷載向量，得到節(jié)點的位移（da0）%求解方程組Lamuda=DetaL/(sqrt(dis'*dis)+1);dis1=Lamuda*dis;Delta_Force=Lamuda*ForceMatrix;%初始荷載向量（Q0）Disp_Transefer=ones(N_Node*3,1);%建立調(diào)節(jié)位移矩陣的位移向量Disp_Transefer(del,:)=0;%將約束節(jié)點的位移值定為0,其他的定位1zz=1;%識別約束fori=1:N_Node*3;ifDisp_Transefer(i,1)==1;Disp_Transefer(i,1)=dis1(zz,1);%總的位移向量zz=zz+1;endendfori=1:N_NodeDisplacement(i,:)=Disp_Transefer(i*3-2:i*3,1);%位移增量，形成n*3的位移向量（da0）endDetaV=DetaV+dis1;%用于求弧長的位移增量（第i荷載步的第一次迭代）DetaVH=DetaVH+dis1;Lamuda1=Lamuda;All_Disp=All_Disp+Displacement;%位移向量更新得到a1（總的位移增量）All_F=All_F+Delta_Force;%外荷載向量更新Q1end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ifn>=2Equation=[K_Global1,ForceMatrix];Equation(del,:)=[];Equation(:,del)=[];n2=size(Equation,2);%列向量的個數(shù)DetaVp=(Equation(:,1:n2-1))\Equation(:,n2);%更新位移Deta_Lamuda1=-DetaL/(sqrt(DetaVp'*DetaVp+1));Deta_Lamuda2=DetaL/(sqrt(DetaVp'*DetaVp+1));DetaVz1=Deta_Lamuda1*DetaVp;%代表位移增量DetaVz2=Deta_Lamuda2*DetaVp;%代表位移增量r1=[DetaVz1',Deta_Lamuda1]';%相對于上一個收斂點的相對向量r2=[DetaVz2',Deta_Lamuda2]';iffangxiangliang*r1>0Deta_Lamuda=Deta_Lamuda1;elseDeta_Lamuda=Deta_Lamuda2;endLamuda=Lamuda+Deta_Lamuda;%更新的荷載系數(shù)Lamuda1=Deta_Lamuda;Deta_V=Deta_Lamuda*DetaVp;%代表位移增量Disp_Transefer=ones(N_Node*3,1);%建立調(diào)節(jié)位移矩陣的位移向量Disp_Transefer(del,:)=0;%將約束節(jié)點的位移值定為0,其他的定位1zz=1;%識別約束fori=1:N_Node*3;ifDisp_Transefer(i,1)==1;Disp_Transefer(i,1)=Deta_V(zz,1);%總的位移向量zz=zz+1;endendfori=1:N_NodeDisplacement(i,:)=Disp_Transefer(i*3-2:i*3,1);%位移增量，形成n*3的位移向量（da0）endDetaV=DetaV+Deta_V;%用于求弧長的位移增量（第i荷載步的第一次迭代）DetaVH=DetaVH+Deta_V;All_Disp=All_Disp+Displacement;%位移向量更新得到a1（總的位移增量）Delta_Force=Deta_Lamuda*ForceMatrix;All_F=All_F+Delta_Force;%外荷載向量更新Q1end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Update_Node1=Update_Node;%C1狀態(tài)Update_Node=Node+All_Disp(:,1:2);%C2狀態(tài)坐標(biāo)位置更新（改）（迭代后的位置）Internal_F1=zeros(N_Node*3,1);%節(jié)點內(nèi)力向量fori=1:N_ElementF_Global=zeros(N_Node*3,1);%全局的荷載向量forj=1:2ELEDISP1(j*3-2:j*3)=Displacement(Element(i,j),:);%整體取出當(dāng)前單元節(jié)點位移向量endN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:);%a0下坐標(biāo)向量增量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=L1(i)+ELEDISP(i,4)-ELEDISP(i,1);Z1=ELEDISP(i,5)-ELEDISP(i,2);LN=(X1^2+Z1^2)^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（i）為a0對應(yīng)下的Ele_F(i,:)=K_Element*ELEDISP(i,:)';%局部坐標(biāo)系下荷載（Q0）作用下的節(jié)點力Ele_F1(i,:)=Ele_F1(i,:)+Ele_F(i,:);C2(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a0下坐標(biāo)向量增量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;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對應(yīng)下的P1endVal=All_F-Internal_F1;%求出上次迭代完成時的殘余應(yīng)力Q1-P1Val(del,:)=0;%受約束點的殘余應(yīng)力項歸零,便于下一次迭代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;%修改位移矩陣形式%whilenorm(Val)/norm(Delta_Force)>=1e-7&k<200%在一個荷載增量步下進(jìn)行的對此迭代whilenorm(Val)>=1e-4&k<200%在一個荷載增量步下進(jìn)行的對此迭代k=k+1;K_Global=zeros(N_Node*3,N_Node*3);%總剛矩陣fori=1:N_ElementN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號N_Section=Element(i,3);%截面的形狀控制參數(shù)C(i,:)=Update_Node(N2,:)-Update_Node(N1,:);%a0下坐標(biāo)向量增量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);%形成總剛K0endEquation=[K_Global,ForceMatrix];%求解方程組Equation(del,:)=[];Equation(:,del)=[];n1=size(Equation,2);%列向量的個數(shù)DetaV1=(Equation(:,1:n1-1))\Equation(:,n1);%zhengshuEquation=[K_Global,Val];%求解方程組Equation(del,:)=[];Equation(:,del)=[];n2=size(Equation,2);%列向量的個數(shù)DetaV2=(Equation(:,1:n2-1))\Equation(:,n2);%zhengshuifk==1Equation=[K_Global1,ForceMatrix];Equation(del,:)=[];%把方程中約束節(jié)點所對應(yīng)的行和列去掉Equation(:,del)=[];n1=size(Equation,2);%方程組中列數(shù)DetaV3=(Equation(:,1:n1-1))\Equation(:,n1);%剛度矩陣求逆后乘以荷載向量,Equation(:,n1)荷載向量，得到節(jié)點的位移（da0）Deta_Lamuda=-((DetaV3')*DetaV2)/(DetaV3'*DetaV1+1)elseDeta_Lamuda=-((DetaV')*DetaV2)/(DetaV'*DetaV1+Lamuda1)endDetaV3DetaVLamuda=Lamuda+Deta_Lamuda;%全過程荷載因子累積Lamuda1=Lamuda1+Deta_Lamuda;%每一荷載步下的荷載因子累積Deta_V=Deta_Lamuda*DetaV1+DetaV2;%代表位移增量DetaV=DetaV+Deta_V;%下一次迭代后的總位移增量DetaVH=DetaVH+Deta_V;All_F=Lamuda*ForceMatrix;%更新頂部荷載Disp_Transefer=ones(N_Node*3,1);Disp_Transefer(del,:)=0;zz=1;fori=1:N_Node*3;ifDisp_Transefer(i,1)==1;Disp_Transefer(i,1)=Deta_V(zz,1);zz=zz+1;endendfori=1:N_NodeCorrec_dis(i,:)=Disp_Transefer(i*3-2:i*3,1);endCorrec_dis1=Correc_dis1+Correc_dis;Internal_F1=zeros(N_Node*3,1);%節(jié)點內(nèi)力向量Update_Node1=Update_Node;Update_Node=Node+All_Disp(:,1:2)+Correc_dis1(:,1:2);%為了求切線剛度矩陣（改）a2下對應(yīng)的節(jié)點位置fori=1:N_ElementF_Global=zeros(N_Node*3,1);forj=1:2ELEDISP1(j*3-2:j*3)=Correc_dis(Element(i,j),:);%取出當(dāng)前單元節(jié)點位移向量endN1=Element(i,1);%i節(jié)點編號N2=Element(i,2);%j節(jié)點編號C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:);%a0下坐標(biāo)向量增量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=L1(i)+ELEDISP(i,4)-ELEDISP(i,1);Z1=ELEDISP(i,5)-ELEDISP(i,2);LN=(X1^2+Z1^2)^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（i）為a1對應(yīng)下的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下坐標(biāo)向量增量