有限单元法作业非线性分析 程序解读.docx
- 文档编号:14783240
- 上传时间:2023-06-27
- 格式:DOCX
- 页数:55
- 大小:886.56KB
有限单元法作业非线性分析 程序解读.docx
《有限单元法作业非线性分析 程序解读.docx》由会员分享,可在线阅读,更多相关《有限单元法作业非线性分析 程序解读.docx(55页珍藏版)》请在冰点文库上搜索。
有限单元法作业非线性分析程序解读
几何非线性大作业荷载增量法
和弧长法程序设计
系(所):
建筑工程系
学号:
*******
******
培养层次:
专业硕士
*******
2015年6月19日
一、几何非线性大作业(Newton-Raphson法)
用荷载增量法(Newton-Raphson法)编写几何非线性程序:
(1)用平面梁单元,可分析平面杆系
(2)算例:
悬臂端作用弯矩。
悬臂梁最终变形形成周长为悬臂梁长度的圆。
1.1Newton-Raphson算法基本思想
图1.1Newton-Raphson算法基本思想
1.2悬臂梁参数
基本参数:
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)
节点信息
2.2求解过程
梁弯成圆形:
理论弯矩M=EIY"=24981.944N.m,直径为0.642m
运用ABAQUS和MATLAB进行求解对比:
图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的计算结果很好。
二、几何非线性大作业(弧长法)
用弧长法编写几何非线性程序,分析荷载位移全过程曲线:
1)用平面梁单元,可分析平面杆系结构
2)算例
(1)受集中荷载的拱:
考察拱的矢跨比、荷载位置对荷载位移曲线的影响。
(2)其他有复杂平衡路径的结构
3)将结果与相关文献进行对比
1.1弧长法基本思想
图2.1弧长法基本思想
1.2拱基本参数
拱参数:
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.1MATLAB输入信息
材料信息单元信息
约束信息(0为约束,1为放松)荷载信息(FX,FY,M)
节点信息
2.2运用ANSYS和MATLAB进行求解对比(两端铰接)
ANSYS中模型:
图2.3ANSYS模型
图2.4MATLAB和ANSYS变形图
2.3加载点荷载位移曲线
图2.5加载点荷载位移曲线
ANSYS求得的极限承载力3042.53,对应位移3.00142
MATLAB求得的极限承载力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)
5mm
3043.8
3.0768
1765.2
7.0816
7.5mm
7623.3
4.0335
-595.82
11.21
10mm
14974
5.4026
-6408.1
14.886
20mm
39791
9.4831
-63049
30.513
从表中可以初步得出:
在一定随着矢跨比的增加,拱仍然呈现跳跃失稳的形式,拱结构的极限承载能力有大幅度的提高;在最低处的承载力呈现出反向,相当于有一个拉力在阻止拱结构发生跳跃失稳,矢跨比越大,拱越不容易发生跳跃失稳。
当拱的矢跨比超过一定范围后,拱将发生复杂的不同于跳跃失稳的失稳形式。
2)MATLAB与ANSYS计算结果对比
图2.8ANSYS和MATLAB对比荷载位移曲线
表2各矢跨比下拱结构的极限荷载对比
参数
矢高
F(N)MAT
位移(m)
F(N)ANA
位移(m)
误差(%)
误差(%)
5mm
3043.8
3.0768
3042.53
3.00142
0.04
2.45
7.5mm
7623.3
4.0335
7624.91
3.96303
-0.02
1.75
10mm
14974
5.4026
14974.3
5.3157
0.00
1.61
20mm
39791
9.4831
39695.7
9.59955
0.24
-1.23
从图中可以看出:
矢跨比在一定范围内,MATLAB与ANSYS计算的荷载位移曲线非常吻合,验证了MATLAB程序的可行性。
当矢跨比为0.15时,ANSYS中将跟踪不到失稳后复杂的平衡路径。
从表中可以得出:
MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近,加载点均为顶点26。
具体为:
矢高5mm,荷载误差为0.04,位移误差为2.45;矢高7.5mm,荷载误差为-0.02,位移误差为1.75;矢高10mm,荷载误差为0,位移误差为-1.61;矢高20mm,荷载误差为0.24,位移误差为-1.23。
实际误差相差很小,在工程允许的范围内是可以接受的。
2.5荷载位置对拱荷载位移曲线的影响
图2.9ANSYS模型简图
1)MATLAB中计算荷载位置对拱荷载位移曲线的影响
图2.10各加载点荷载位移曲线
表3改变加载点拱结构的极限荷载
参数
加载点
极值点
F(N)
位移(m)
最低点
F(N)
位移(m)
26
3043.8
3.0768
1765.2
7.0816
16
3570
3.1891
2258.8
6.116
11
4728
2.88
3220.5
4.7959
4
14317
1.2826
10569
1.7829
误差=100*(MATLAB-ANSYS)/ANSYS
从表中可以初步得出:
随着加载点位置越靠近支座处,拱结构的极限承载能力有大幅度的提高;在最低处的承载力也大幅度提高。
当加载点位置靠近支座时,拱的承载力增加幅度最大,拱的稳定性很强,不容易发生失稳。
2)MATLAB与ANSYS计算结果对比
图2.11ANSYS和MATLAB对比荷载位移曲线
表4各加载点拱结构的极限荷载
参数
矢高
F(N)MAT
位移(m)
F(N)ANA
位移(m)
误差(%)
误差(%)
26
3043.8
3.0768
3042.53
3.00142
0.04
2.45
16
3570
3.1891
3569.73
3.24865
0.01
-1.87
11
4728
2.88
4728.71
2.91862
-0.02
-1.34
4
14317
1.2826
14324.8
1.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.05,位移误差为-1.17。
实际误差相差很小,在工程允许的范围内是可以接受的。
2.6两端铰接和固接对拱荷载位移曲线的影响
矢高为5mm时,拱两端为固接和铰接时的荷载位移曲线如下:
图2.12ANSYS和MATLAB固接和铰接的荷载位移曲线
从图中可以看出:
拱的边界条件对其的失稳形式有很大影响。
两端固接拱的稳定性明显优于两端铰接拱,承载能力也大幅度提高。
固接拱不会发生跳跃失稳的形式,刚度在初始阶段会减小,待到达一定程度后刚度又会增加。
而两端铰接拱会发生跳跃失稳的形式。
2.7参数m对拱失稳形式的影响
文献中给出:
m是一个由拱截面在竖向平面内的回转半径r和拱的初始矢高h无确定的无量纲参数。
当m>=1时,扁拱不会出现跳跃屈曲,当0 2.13不同m值加载点的荷载位移曲线 2.14不同m值变形后拱曲线 从MATLAB的计算结果中可以验证: 不同的m系数对应拱不同的失稳形式。 当m>=1时,扁拱不会出现跳跃屈曲,当0 但拱的最终变形图非常接近,只是此时拱的失稳形式变了。 2.8具有复杂失稳形式的拱 铰支圆拱 该结构及其几何参数、物理性质均示于图4a中。 中心受集中荷载。 这个结构是研究分歧问题的经典题目。 将半跨结构划分为8个单元,得到图4b的基本路径和分歧路径,并与JChreseielewski和Rsehmiot的结果进行了比较。 本文对此结构也进行了缺陷分析。 拱的基本参数: L=254cm,R=381cm,I=41.62cm4,A=6.45cm2,E=6898kN/cm2。 文献中的计算结果。 采用MATLAB和ANSYS对其进行求解,得到其荷载位移曲线: 图2.15ABAQUS模型 图2.16ABAQUS变形图 图2.17ANSYS、MATLAB、ABAQUS加载点荷载位移曲线 从MATLAB、ANSYS、ABAQUS计算的荷载位移曲线可以看出: 两者的荷载位移曲线基本吻合。 MATLAB的计算结果可以看出在后期其承载能力会有较大提高,与文献中的荷载位移曲线趋势相同,所以验证出程序的可靠性。 ABAQUS不能模拟出后续段曲线也许是单元划分过少。 图2.18MATLAB加载点荷载位移曲线 第二次极值点会超过第一次极值点所对应的荷载,与文献一致,荷载点也接近。 加入初始缺陷: L/1000,L/2000初始缺陷后观察加载点的荷载位移曲线变化趋势。 图2.19ANSYS加入初始缺陷后的加载点荷载位移曲线 2.20初始缺陷为0.0001时的荷载位移曲线 加入初始缺陷后,拱的极限承载能力明显降低。 且失稳形式也发生了变化,初始缺陷的大小对其荷载位移曲线有明显影响。 所以在工程设计中应考虑结构或构件的初始缺陷,使结构的设计更加合理,安全。 为了研究初始缺陷对拱失稳路径的影响,应用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,对应的屈曲299915N和 288078N。 从屈曲模态中可以看出,两种软件的前二阶模态趋势吻合,屈曲系数和极限荷载也是吻合的较好。 第三阶模态出现不一样的变形趋势,屈曲系数和极限荷载也是也相差比较大,但计算时只引入一阶屈曲模态。 图2.24ANSYS、ABAQUS、MATLAB加载点荷载位移曲线 从图中可以看出: ANSYS对缺陷的模拟效果比较好,与文献中的比较接近ABAQUS模拟的极限荷载稍低于ANSYS,而MATLAB模拟的极限荷载远低于ANSYS,但是曲线最终都在位移为300多mm时交于一点。 还是有一定规律性。 图2.25ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线 始缺陷并一定都会降低承载力,也会有对结构的承载能力有益的初始缺陷。 ANSYS的计算结果可以看出,当初始缺陷为1/2000和-1/2000时,其承载能力不变。 由于为对称结构,所以缺陷的正负影响不大。 图2.26ANSYS和ABAQUS引入初始缺陷加载点荷载位移曲线 表6对比ANNSYS和ABAQUS的极限荷载值和其对应位移值 参数 矢高 F(N)ANS 位移(m) F(N)ABA 位移(m) 误差(%) 误差(%) 1/1000 58444.7 68.9799 55795.628 79.9318 4.53261 -15.8769 1/2000 60579.9 70.1384 57924.958 65.4542 4.38255 6.67851 误差=100*(ABAQUS-ANSYS)/ABAQUS 表中可以得出: ABAQUS求解出的机线荷载小于ANSYS,单对应的位移大于ANSYS对应的值。 这可能与单元划分的个数,求解精度有关,但在工程允许的范围内,还是可以接受的。 ABAQUS中迭代步的跳跃很快,位移增加速度很快,其路径不是很准确,可能是由于其单元划分比较少。 体会: 1)注意坐标系的转化和力、位移更新时所对应的状态(C1--C2) 2)拱是否发生跳跃失稳与矢跨比、矢高与截面旋转半径有关。 矢跨比太大不会发生跳跃失稳;m大于1时不能发生跳跃失稳。 3)有些拱在ANSYS中不能得出下降段,可见ANSYS中对拱的跨度矢高、截面参数的比值有一定要求。 内部计算和程序中有一些差别。 附录 子程序一: 刚度组装矩阵 functionK_G=Assemble(K_Element,Element,N_Node) K_G=zeros(N_Node*3,N_Node*3); fori=1: 2 n1=Element(1,i); forj=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); end end end 子程序二: 整体坐标刚度矩阵 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 子程序三: 局部坐标刚度矩阵 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法 clc clear Node=xlsread('Data.xls','Node'); Element=xlsread('Data.xls','Element'); Boundary=xlsread('Data.xls','Boundary'); Section=xlsread('Data.xls','Section'); Force=xlsread('Data.xls','Force'); %读取相关数据 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) %初始位移及转角为0 All_Disp=zeros(N_Node,3);%初始位移和转角为零(迭代后的节点位移) All_F=zeros(N_Node*3,1);%初始荷载向量为零(存放节点荷载向量) Internal_F=zeros(N_Node*3,1);%节点内力向量 ForceMatrix=zeros(N_Node*3,1);%总荷载向量 C=zeros(N_Element,2); L=zeros(N_Element,1); fori=1: N_Force ForceMatrix(Force(i,1)*3-2: Force(i,1)*3,1)=Force(i,2: 4)';%把节点荷载向量读入一个矩阵中,形成列向量=总的自由度个数 end del=[]; fori=1: N_Boundary; ifBoundary(i,2)==0; m=3*Boundary(i,1)-2; del=[del,[m]];%受约束节点位移为0时所对应的指标数值1 end ifBoundary(i,3)==0; m=3*Boundary(i,1)-1; del=[del,[m]];%受约束节点位移为0时所对应的指标数值2 end ifBoundary(i,4)==0; m=3*Boundary(i,1); del=[del,[m]];%受约束节点位移为0时所对应的指标数值3 end end %求出约束节点的标号,便于刚度、荷载矩阵归0 Update_Node=Node+Displacement(: 1: 2);%更新后的节点坐标向量(x,y坐标) Ele_F=zeros(N_Element,6);%单元节点荷载向量 ELEDISP=zeros(N_Element,6);%单元节点位移向量 Ele_F1=zeros(N_Element,6);%单元节点荷载刚度矩阵中向量 Ele1_F=zeros(1,6); ELEDISP1=zeros(1,6); qq(1,1)=0; pp(1,1)=0; forn=0: Allstep-1 n=n+1 K_Global=zeros(N_Node*3,N_Node*3);%总刚矩阵 fori=1: N_Element N1=Element(i,1);%i节点编号 N2=Element(i,2);%j节点编号 N_Section=Element(i,3);%截面的形状控制参数 C(i,: )=Update_Node(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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限单元法作业非线性分析 程序解读 有限 单元 作业 非线性 分析 程序 解读