有限单元法作业非线性分析+程序

1、几何非线性大作业荷载增量法和弧长法程序设计系（所）：建筑工程系学 号：1432055姓 名：焦 联 洪培养层次：专业硕士指导老师：吴 明 儿 2015年6月19日 一、 几何非线性大作业（ Newton-Raphson法）用荷载增量法（Newton-Raphson法）编写几何非线性程序：（1）用平面梁单元，可分析平面杆系（2）算例：悬臂端作用弯矩。悬臂梁最终变形形成周长为悬臂梁长度的圆。1.1 Newton-Raphson算法基本思想图1.1 Newton-Raphson算法基本思想1.2 悬臂梁参数基本参数：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） 节点信息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的计算结果很好。二、 几何非线性大作业（ 弧长法）用弧长法编写几何非线性程序，分析荷载位移全过程曲线：1) 用平面梁单元，可分析平面杆系结构2) 算例(1)受集中荷载的拱：考察拱的矢跨比、荷载位置对荷载位移曲线的影响。(2)其他有复杂平衡路径的结构3) 将结果与相关文献进行对比1.1 弧长法基本思想图2.1 弧长法基本思想1.2 拱基本参数拱参数

4、：L=100m, A=0.32m2 ,I=1m4 ,E=1.0e7N/m2，F=-5000N，拱曲线 y=5×sin(3.1415926*x/L)将拱结构分成25个单元，如图2所示图2.2 拱单元信息2.1 MATLAB输入信息材料信息 单元信息约束信息（0为约束，1为放松） 荷载信息（FX,FY,M）节点信息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 各矢跨比下拱结构的极限荷载 参数矢高极值点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从表中可以初步得出：在一定随着矢跨比的增加，拱仍然呈现跳跃失稳的形式，拱结构的极限承载能力有大幅度的提高；在最低处的承载力呈现出反向，相当于有一个拉力在阻止拱结构发生跳跃失稳，矢跨比越大，拱越不容易发生跳跃失稳。当拱的矢跨比超过一定范围后，拱将发生复杂的不同于跳跃失稳的失稳形式。2）MATLAB与ANSYS计算结果对比图2.8 ANSYS和MATLAB对比荷载位移曲线表2 各矢跨比下拱结构的极限荷载对比 参数矢高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从图中可以看出：矢跨比在一定范围内，MATLAB与ANSYS计算的荷载位移曲线非常吻合，验证了MATLAB程序的可行性。当矢跨比为0.15时，ANSYS中将跟踪不到失稳后复杂的平衡路径。从表中可以得出：MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近，加载点均为顶点26。具体为：矢高5mm，荷载误差为0.04，位移误差为2.45；矢高7.5mm，荷载

8、误差为-0.02，位移误差为1.75；矢高10mm，荷载误差为0，位移误差为-1.61；矢高20mm，荷载误差为0.24，位移误差为-1.23。实际误差相差很小，在工程允许的范围内是可以接受的。2.5 荷载位置对拱荷载位移曲线的影响图2.9 ANSYS模型简图1）MATLAB中计算荷载位置对拱荷载位移曲线的影响图2.10 各加载点荷载位移曲线表3 改变加载点拱结构的极限荷载 参数加载点极值点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从表中可以初步得出：随着加载点位置越靠近支座处，拱结构的极限承载能力有大幅度的提高；在最低处的承载力也大幅度提高。当加载点位置靠近支座时，拱的承载力增加幅度最大，拱的稳定性很强，不容易发生失稳。2）MATLAB与ANSYS计算结果对比图2.11 ANSYS和MATLAB对比荷载位移曲线表4 各加载点拱结构的极限荷载 参数矢高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。实际误差相差很小，在工程允许的范围内是可以接受的。2.6 两端铰接和固接对拱荷载位移曲线的影响矢高为5mm时，拱两端为固接和铰接时的荷载位移曲线如下：图2.12 ANSYS和MATLAB固接和铰接的荷载位移曲线从图中可以看出：拱的边界条件对其的失稳形式有很大影响。两端固接拱的稳定性明显优于两端铰接拱，承载能力也大幅度提高。固接拱不会发生跳跃失稳的形式，刚度在初始阶段会减小，待到达一定程度后刚度又会增加。而两端铰接拱会发生跳跃失稳的形式。2.7 参数m对拱失稳形式的影响文献中给出:m是一个由拱截面在竖向平面内的回转半径r 和拱的初始矢高h无确定的无量纲参数。当m>

12、=1 时，扁拱不会出现跳跃屈曲, 当0<m<1时，扁拱有可能发生跳跃屈曲，而影响扁拱是否发生跳跃屈曲的主要因素是m值和荷载参数。 2.13 不同m值加载点的荷载位移曲线2.14 不同m值变形后拱曲线从MATLAB的计算结果中可以验证：不同的m系数对应拱不同的失稳形式。当m>=1 时，扁拱不会出现跳跃屈曲，当0<m<1时，扁拱有可能发生跳跃屈曲。但拱的最终变形图非常接近，只是此时拱的失稳形式变了。 2.8 具有复杂失稳形式的拱铰支圆拱该结构及其几何参数、物理性质均示于图4a 中。中心受集中荷载。这个结构是研究分歧问题的经典题目。将半跨结构划分为8个单元, 得到图4b

13、的基本路径和分歧路径, 并与JChreseielewski 和Rsehmiot的结果进行了比较。本文对此结构也进行了缺陷分析。拱的基本参数：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。文献中的计算结果。采用MATLAB和ANSYS对其进行求解，得到其荷载位移曲线：图2.15 ABAQUS模型图2.16 ABAQUS变形图图2.17 ANSYS、MATLAB、ABAQUS加载点荷载位移曲线从MATLAB、ANSYS、ABAQUS计算的荷载位移曲线可以看出：两者的荷载位移曲线基本吻合。MATLAB的计算结果可以看出在后期其承载能力会有较大提高

14、，与文献中的荷载位移曲线趋势相同，所以验证出程序的可靠性。ABAQUS不能模拟出后续段曲线也许是单元划分过少。图2.18 MATLAB加载点荷载位移曲线 第二次极值点会超过第一次极值点所对应的荷载，与文献一致，荷载点也接近。加入初始缺陷：L/1000, L/2000初始缺陷后观察加载点的荷载位移曲线变化趋势。图2.19 ANSYS加入初始缺陷后的加载点荷载位移曲线2.20 初始缺陷为0.0001时的荷载位移曲线加入初始缺陷后，拱的极限承载能力明显降低。且失稳形式也发生了变化，初始缺陷的大小对其荷载位移曲线有明显影响。所以在工程设计中应考虑结构或构件的初始缺陷，使结构的设计更加合理，安全。为了研

15、究初始缺陷对拱失稳路径的影响，应用ABAQUS和ANSYS以及MATLAB中加水平力模拟拱结构初始缺陷下的荷载位移曲线。为了探究ABAQUS和ANSYS计算结果，取其前三阶模态进行对比分析： 2.21 一阶屈曲模态 ABAQUS和MATLAB中的一阶屈曲系数为0.55884和0.564512，对应的屈曲荷载为74325.72N 和75080.096N。2.22 二阶屈曲模态 ABAQUS和MATLAB中的二阶屈曲系数为1.2259和1.253，对应的屈曲163044.7N 和166649N。2.23 三阶屈曲模态ABAQUS和MATLAB中的三阶屈曲系数为2.166和2.255，对应的屈曲29

16、9915N 和288078N。从屈曲模态中可以看出，两种软件的前二阶模态趋势吻合，屈曲系数和极限荷载也是吻合的较好。第三阶模态出现不一样的变形趋势，屈曲系数和极限荷载也是也相差比较大，但计算时只引入一阶屈曲模态。图2.24 ANSYS、ABAQUS、MATLAB加载点荷载位移曲线从图中可以看出：ANSYS对缺陷的模拟效果比较好，与文献中的比较接近ABAQUS模拟的极限荷载稍低于ANSYS，而MATLAB模拟的极限荷载远低于ANSYS,但是曲线最终都在位移为300多mm时交于一点。还是有一定规律性。图2.25 ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线 始缺陷并一定都会降低承载力，也

17、会有对结构的承载能力有益的初始缺陷。ANSYS的计算结果可以看出，当初始缺陷为1/2000和-1/2000时，其承载能力不变。由于为对称结构，所以缺陷的正负影响不大。图2.26 ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线表6 对比ANNSYS和ABAQUS的极限荷载值和其对应位移值 参数矢高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对应的值。这可能与单元划分的个数，求解精度有关，但在工程允许的范围内，还是可以接受的。ABAQUS中迭代步的跳跃很快，位移增加速度很快，其路径不是很准确，可能是由于其单元划分比较少。体会：1）注意坐标系的转化和力、位移更新时所对应的状态（C1-C2）2) 拱是否发生跳跃失稳与矢跨比、矢高与截面旋转半径有关。矢跨比太大不会发生跳跃失稳；m大于1时不能发生跳跃失稳。3）有些拱在ANSYS中不能得出下降段，可见ANSYS中对拱的跨度矢高、截面参数的比值有一定要求。内部计算和程序中有一些差

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');%读取相关

26、数据Allstep=1000; %迭代步数N_Node=size(Node,1); %节点个数 N_Element=size(Element,1); %单元个数N_Force=size(Force,1); %节点力个数N_Boundary=size(Boundary,1); %约束节点个数Displacement=zeros(N_Node,3); %位移向量（a0）%初始位移及转角为0All_Disp=zeros(N_Node,3); %初始位移和转角为零（迭代后的节点位移）All_F=zeros(N_Node*3,1); %初始荷载向量为零 （存放节点荷载向量）Internal_F=zero

27、s(N_Node*3,1); %节点内力向量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)' %把节点荷载向量读入一个矩阵中,形成列向量=总的自由度个数enddel=;for i=1:N_Boundary; if Boundary(i,2)=0; m=3*Boundary(i,1)-2; del=del,m; %受约束节点位移为0时所对

28、应的指标数值1 end if Boundary(i,3)=0; m=3*Boundary(i,1)-1; del=del,m; %受约束节点位移为0时所对应的指标数值2 end if Boundary(i,4)=0; m=3*Boundary(i,1); del=del,m; %受约束节点位移为0时所对应的指标数值3 endend%求出约束节点的标号，便于刚度、荷载矩阵归0Update_Node=Node+Displacement(:,1:2); %更新后的节点坐标向量（x，y坐标）Ele_F=zeros(N_Element,6); %单元节点荷载向量ELEDISP=zeros(N_Eleme

29、nt,6); %单元节点位移向量Ele_F1=zeros(N_Element,6); %单元节点荷载刚度矩阵中向量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节点编号 N2=Element(i,2); %j节点编号 N_Section=Element(i,3); %截面的形状控制参数 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); %建立调节位移矩阵的位移向量 Disp_Transefer(del,:)=0; %将约束节点的位移值定为0,其他的定位1 Equation(del,:)=;%把方程中约束节点所对应的行和列去掉 Equation(:,del)=; %引入约束条件修改方程组 n1=size(Equation,2); % 方程组中列数 dis1=(Equation(:,1:n1-1)Equati

32、on(:,n1); % 刚度矩阵求逆后乘以荷载向量, Equation(:,n1)荷载向量，得到节点的位移（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); %节点内力向量 Update_Node1=Update_Node; %C1状态 Update_Node=Node+All_Disp(:,1:2); %C2状态坐标位置更新（改）（迭代后的位置 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、 取出当前单元节点位移向量 end N1=Element(i,1); %i节点编号 N2=Element(i,2); %j节点编号 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）作用下的节点力 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节点荷载 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j节点荷载 Internal_F1=Internal_F

38、1+F_Global; %a1对应下的P1 end Val=Internal_F1-All_F %求出上次迭代完成时的残余应力Q1-P1 Correc_dis=zeros(N_Node,3); %构造新的节点位移向量，每次更新 Correc_dis1=zeros(N_Node,3); %构造新的节点位移向量，叠加位移增量 N_dis=size(dis1,1); %未受约束的节点位移数目,不为零的节点位移数目 dis3=zeros(N_dis,1); %构造一个向量 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); %节点内力向量 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),:); %取出当前单元节点位移向量 end N1=Element(i,1); %i节点编号 N2=Element(i,

43、2); %j节点编号 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节点荷载 F_Global(3*N2-2:3*N2,1)=Ele1_F(4:6); %j节点荷载 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&