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

1、幾何非線性大作業(yè)荷載增量法和弧長(zhǎng)法程序設(shè)計(jì)系（所）：建筑工程系學(xué) 號(hào)：1432055姓 名：焦 聯(lián) 洪培養(yǎng)層次：專(zhuān)業(yè)碩士指導(dǎo)老師：吳 明 兒 2015年6月19日 一、 幾何非線性大作業(yè)（ Newton-Raphson法）用荷載增量法（Newton-Raphson法）編寫(xiě)幾何非線性程序：（1）用平面梁?jiǎn)卧?，可分析平面桿系（2）算例：懸臂端作用彎矩。懸臂梁最終變形形成周長(zhǎng)為懸臂梁長(zhǎng)度的圓。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 懸臂梁?jiǎn)卧畔冶哿悍殖?0個(gè)單元，如圖1.2所示2.1 MATLAB輸入信息 材料信息 單元信息 約束信息（0為約束，1為放松） 荷載信息（FX,FY,M） 節(jié)點(diǎn)信息2.2 求解過(guò)程梁彎成圓形：理論彎矩M=EIY"=24981.944N.m ,直徑為0.642m運(yùn)用ABAQUS和MATLAB進(jìn)行求解對(duì)比：圖1.3 加載圖 圖1.4 ABAQUS變形圖圖1.5 MATLAB變形曲線ABAQUS和MATLAB變形對(duì)比，最終在理論荷載作用下都彎成了一個(gè)圓，其直徑為0.64716m，與理論值相對(duì)比值為：（0.64716-0.642）/0.642=0.

3、00804.非常接近。2.3 加載點(diǎn)荷載位移曲線圖1.5 加載點(diǎn)Y方向的荷載位移曲線加載點(diǎn)的最大豎向位移分別為1.4525m和1.45246m，相對(duì)比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，說(shuō)明MATLAB的計(jì)算結(jié)果很好。二、 幾何非線性大作業(yè)（ 弧長(zhǎng)法）用弧長(zhǎng)法編寫(xiě)幾何非線性程序，分析荷載位移全過(guò)程曲線：1) 用平面梁?jiǎn)卧?，可分析平面桿系結(jié)構(gòu)2) 算例(1)受集中荷載的拱：考察拱的矢跨比、荷載位置對(duì)荷載位移曲線的影響。(2)其他有復(fù)雜平衡路徑的結(jié)構(gòu)3) 將結(jié)果與相關(guān)文獻(xiàn)進(jìn)行對(duì)比1.1 弧長(zhǎng)法基本思想圖2.1 弧長(zhǎng)法基本思想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個(gè)單元，如圖2所示圖2.2 拱單元信息2.1 MATLAB輸入信息材料信息 單元信息約束信息（0為約束，1為放松） 荷載信息（FX,FY,M）節(jié)點(diǎn)信息2.2 運(yùn)用ANSYS和MATLAB進(jìn)行求解對(duì)比（兩端鉸接）ANSYS中模型:圖2.3 ANSYS模型圖2.4 MATLAB和ANSYS變形圖2.3 加載點(diǎn)荷載位移曲線圖2.5 加載點(diǎn)荷載位移曲線ANSYS求得的極限承載力3042.53,對(duì)應(yīng)位移3.00142MATLAB求得的

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

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

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

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

9、105691.7829 誤差=100*（MATLAB-ANSYS）/ANSYS從表中可以初步得出：隨著加載點(diǎn)位置越靠近支座處，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力也大幅度提高。當(dāng)加載點(diǎn)位置靠近支座時(shí)，拱的承載力增加幅度最大，拱的穩(wěn)定性很強(qiáng)，不容易發(fā)生失穩(wěn)。2）MATLAB與ANSYS計(jì)算結(jié)果對(duì)比圖2.11 ANSYS和MATLAB對(duì)比荷載位移曲線表4 各加載點(diǎn)拱結(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計(jì)算的荷載位移曲線非常吻合，驗(yàn)證了MATLAB程序的可行性。從表中可以得出：MATLAB與ANSYS計(jì)算中拱的極限荷載和極限荷載時(shí)所對(duì)應(yīng)的位移非常接近。具體為：26加載點(diǎn)，荷載誤差為0.04，位移誤差為2.45；16加載點(diǎn)，荷載誤差為0.01，位移誤差為-1.87；11加載點(diǎn)，荷載誤差為-0.02，位移誤差為-1.34；4加載點(diǎn)，荷載誤差為-0.05

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

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

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

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

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

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

17、有對(duì)結(jié)構(gòu)的承載能力有益的初始缺陷。ANSYS的計(jì)算結(jié)果可以看出，當(dāng)初始缺陷為1/2000和-1/2000時(shí)，其承載能力不變。由于為對(duì)稱(chēng)結(jié)構(gòu)，所以缺陷的正負(fù)影響不大。圖2.26 ANSYS和ABAQUS引入初始缺陷加載點(diǎn)荷載位移曲線表6 對(duì)比ANNSYS和ABAQUS的極限荷載值和其對(duì)應(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

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

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子程序二：整體坐標(biāo)剛度矩陣function K=beam2d_stiffness530(E,A,I,L,cs,Ele_F1);F=Ele_F1(1,4);M1=Ele_F1

20、(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*E

21、*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子程序三：局部坐標(biāo)剛度矩陣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)數(shù)

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

27、(N_Node*3,1); %節(jié)點(diǎn)內(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é)點(diǎn)荷載向量讀入一個(gè)矩陣中,形成列向量=總的自由度個(gè)數(shù)enddel=;for i=1:N_Boundary; if Boundary(i,2)=0; m=3*Boundary(i,1)-2; del=del,m; %受約束節(jié)點(diǎn)位移為0時(shí)所對(duì)應(yīng)

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

29、t,6); %單元節(jié)點(diǎn)位移向量Ele_F1=zeros(N_Element,6); %單元節(jié)點(diǎn)荷載剛度矩陣中向量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é)點(diǎn)編號(hào) N2=Element(i,2); %j節(jié)點(diǎn)編號(hào) N_Section=Element(i,3); %截面的形狀控制參數(shù) C(i,:)=Update_Node(

30、N2,:)-Update_Node(N1,:); %a0下坐標(biāo)向量增量 L(i)=norm(C(i,:); %變形后的長(zhǎng)度 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%整體剛度矩陣 De

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

32、n(:,n1); % 剛度矩陣求逆后乘以荷載向量, Equation(:,n1)荷載向量，得到節(jié)點(diǎn)的位移（da0） %求解方程組 zz=1; %識(shí)別約束 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é)點(diǎn)內(nèi)力向量 Update_Node1=Update_Node; %C1狀態(tài) Update_Node=Node+All_Disp(:,1:2); %C2狀態(tài)坐標(biāo)位置更新（改）（迭代后的位置 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、取出當(dāng)前單元節(jié)點(diǎn)位移向量 end N1=Element(i,1); %i節(jié)點(diǎn)編號(hào) N2=Element(i,2); %j節(jié)點(diǎn)編號(hào) C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐標(biāo)向量增量 L1(i)=norm(C1(i,:); %變形后的長(zhǎng)度 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

35、(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_Element

36、=beam2d_stiffness520(E,A,I,L1(i),cs1,Ele_F1(i,:); %L（i）為a0對(duì)應(yīng)下的 Ele_F(i,:)=K_Element*ELEDISP(i,:)' %局部坐標(biāo)系下荷載（Q0）作用下的節(jié)點(diǎn)力 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,:); %變形后的長(zhǎng)度 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é)點(diǎn)荷載 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j節(jié)點(diǎn)荷載 Internal_F1=Internal_F1

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

39、對(duì)此迭代 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下坐標(biāo)向量增量 L(i)=norm(C(i,:); %變形后的長(zhǎng)度 cs=C(i,:)/L(i); E=Section(N_Section,1); A=Section(N_Section,2); I=Section(N_Section,3); K_Elemen

40、t=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; %在殘余應(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)Equa

41、tion(:,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é)點(diǎn)內(nèi)力向量 Update_Node1=Update_Node; Update_Node=N

42、ode+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),:); %取出當(dāng)前單元節(jié)點(diǎn)位移向量 end N1=Element(i,1); %i節(jié)點(diǎn)編號(hào) N2=Element(i,2

43、); %j節(jié)點(diǎn)編號(hào) C1(i,:)=Update_Node1(N2,:)-Update_Node1(N1,:); %a0下坐標(biāo)向量增量 L1(i)=norm(C1(i,:); %變形后的長(zhǎng)度 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對(duì)應(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,:); %變形后的長(zhǎng)度 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é)點(diǎn)荷載 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j節(jié)點(diǎn)荷載 Internal_F1=Internal_F1+F_Global; %a1對(duì)應(yīng)下的P1 end Val=Internal_F1-All_F; %荷載增量Q1-P2 Val(del,:)=0; e

47、nd 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('懸臂梁加載點(diǎn)荷載位移曲線');xlabel('位移（m）');ylabel('荷載（N）');legend('點(diǎn)11荷載位移曲線');grid主程序二：弧長(zhǎng)

48、法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é)點(diǎn)個(gè)數(shù) N_Element=size(Element,1); %單元個(gè)數(shù)N_Force=size(Force,1); %節(jié)點(diǎn)力個(gè)數(shù)N_Boundary=size(Boundary,1); %約束節(jié)點(diǎn)個(gè)數(shù)Displacement=zeros(N_