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

1、幾何非線性大作業(yè)荷載增量法和弧長法程序設(shè)計系（所）：建筑工程系學(xué) 號：1432055姓 名：焦 聯(lián) 洪培養(yǎng)層次：專業(yè)碩士指導(dǎo)老師：吳 明 兒 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.00804.

3、非常接近。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)文獻進行對比1.1 弧長法基本思想圖2.1 弧長法基本思想1.2 拱基本參數(shù)拱參數(shù)：L=100

4、m, A=0.32m2 ,I=1m4 ,E=1.0e7N/m2，F(xiàn)=-5000N，拱曲線 y=5sin(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,對應(yīng)位移3.00142MATLAB求得的極限承載力3043.8

5、, 對應(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-63049

6、30.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.8 ANSYS和MATLAB對比荷載位移曲線表2 各矢跨比下拱結(jié)構(gòu)的極限荷載對比 參數(shù)矢高F(N)MAT位移（m）F(N)ANA位移（m）誤差（%）誤差（%）5mm3043.83.07683042.533.001420.042.457.5mm7623.34.03

7、357624.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，位移誤

8、差為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.2826105691.7829

9、 誤差=100*（MATLAB-ANSYS）/ANSYS從表中可以初步得出：隨著加載點位置越靠近支座處，拱結(jié)構(gòu)的極限承載能力有大幅度的提高；在最低處的承載力也大幅度提高。當(dāng)加載點位置靠近支座時，拱的承載力增加幅度最大，拱的穩(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.01-1.871147

10、282.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

11、。實際誤差相差很小，在工程允許的范圍內(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ù)。當(dāng)m=1 時，扁拱不會出現(xiàn)跳躍屈曲,

12、 當(dāng)0m=1 時，扁拱不會出現(xiàn)跳躍屈曲，當(dāng)0m1e-7 & k500 %在一個荷載增量步下進行的對此迭代 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,:); %變形后的長度 cs=C(i,:)/L(i); E=Section(N_Section,1); A=Sectio

13、n(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); %形成總剛k1 end Equation=K_Global,Val; %在殘余應(yīng)力下的位移求解 Disp_Transefer=ones(N_Node*3,1); Disp_Transefer(del,:)=0; Equation(del,:)=; Equation(:,del)=; n1

14、=size(Equation,2); Dis2=-(Equation(:,1:n1-1)Equation(:,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);

15、 %節(jié)點內(nèi)力向量 Update_Node1=Update_Node; Update_Node=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),:); %取出當(dāng)前單元節(jié)點位移向量 end N1=Elem

16、ent(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; EL

17、EDISP(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=beam2d_stiffness520(E

18、,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;

19、 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)下的P1 end Val=Internal_F1-All_F; %荷載增量Q1-

20、P2 Val(del,:)=0; 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=xlsr

21、ead(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); %位移向量

22、（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);for i=1:N_Force ForceMatrix(Force(i,1)*3-2:Force(i,1)*3,1)=Force(i,2:4); %把節(jié)點荷載向量讀入一個矩enddel=;for i=1:N_Boundary; if Boundary(i,2)=0; m=3*Boundary(i,1)-2; del=del,m; %受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值1 end if Boundary(i,3)=0; m=3*Boundary(i,1)-1; del=del,m; %受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值2 end if Boundary(i,4)=0; m=3*Boundary(i,1); del=del,m; %受約束節(jié)點位移為0時所對應(yīng)的指標(biāo)數(shù)值3 enden