有限元程序设计计算书附相关程序说明.docx

课程设计说明书有限元及程序设计课程设计专 业：工程力学学 生：某某某指导老师：许杨健河北工程大学土木工程学院2011年6月28日本有限元课程设计在有限元及程序设计程序例6的基础上，主要增加了一个单元，扩展 了计算集中力偶及指针输出方面的功能,并编辑了相应的程序语言.关键词：有限元法平面钢架单元刚度矩阵荷载AbstractThis design of Finite Element course is on the basis of example 6 which derive from Element Method and Programming ,A major increase a unit , expanded the calculation concentrated couples exerting on the annular gear and the function of pointer output and edit the corresponding programming language.Key words: Finite Element Method Planar steel element stiffness matrix load目录0引言31计算平面钢架的刚度矩阵31.1计算局部坐标系中的单元刚度矩阵41.2.1计算各单元的坐标变换矩阵51.2.2计算各单元的整体刚度矩阵51.2.3计算平面钢架的整体刚度矩阵62计算非节点荷载的等效节点荷载.72.1计算局部坐标系下等效节点荷载72.1.1单元有均布荷载72.1.2单元有集中荷载82.1.3单元有集中力偶92.2计算整体坐标系下等效节点荷载102.3计算钢架的整体等效节点荷载113计算各节点的位移113.1列平衡方程113.2引入边界条件124求单元内力134.1求节点力表达式134.2计算各单元的单元节点力144.2.1节点位移引起的单元节点力144.2.2非节点荷载引起的固端反力155画内力图165.1节点位移引起的内力图165.2非节点荷载引起的内力图175.3该有限元程序设计框架图19主要结论20心得体会20鸣谢20参考文献20河北工程大学土木工程学院课程设计说明书2011年附件1程序源代码.21有限元课程设计学生：某某某指导老师：许杨健河北工程大学土木工程学院工程力学专业0引言本课程设计是配合理论课教学的一个重要的实践教学环节，起到巩固课堂和书本上所学理论 知识、培养调试和编写有限元程序能力以及启发创新思维的效果。通过本课程设计能使我们系统 地，综合地和灵活地运用本课程学到的知识，并有效训练我们在应用有限元法去分析解决力学问 题、查阅文献和使用技术资料、计算、绘图、编写设计文件及独立工作等方面的能力，为将来的 毕业设计或毕业论文打下坚实的基础。1计算平面钢架的刚度矩阵划分单元，标单元号及节点号如图（1-1）功 9. 6kN/m 160kNI 必豆玉玉* *口一空 3K1OmK图（i-i）平面钢架及其单元节点取钢架的局部坐标系和整体坐标系如图（1-2）图（1-2 ）整体，局部坐标1.1计算局部坐标系中的单元刚度矩阵K'单元各杆件的尺寸相同，惯性矩I,弹性模量E和截面面积A相同，且A=0.5m2,I=0.4166m4,E=2.1X 107kN/m2。由单元刚度矩阵如式(1-1)(EA00lo012 EI6 EI1312006 EI4 EIKKe =121EA00l12 EI6 EI0131206 EI2 EI0 k121EA010012 EI6 EI131206 EI2 EI121EA001012 EI6 EI131206 EI4 EI121 )(1-1)(10.500-10.500、00.1050.5250-0.1050.52500.5253.50-0.5251.75K ®= K =K =-10.50010.500x 100-0.105-0.52500.105-0.525k 00.5251.750-0.5253.5 /将数据代入式（1-1）得(1-2)1.2计算整体坐标系下单元刚度矩阵1.2.1计算各单元的坐标变换矩阵由坐标变换矩阵T为:'cos 以sin以0000 )-sin 以cos以0000001000T=000cos以sin以0000-sin 以cos以0、000001 /(1-3)单元坐标变换矩阵：因为二者相同，且a =0，代入（1-3）式(1-4)单元的坐标变换矩阵：因为a =900代入（1-3）式得(010000，-100000001000000010000-100、000001 /T=T=I （单位矩阵）(1-5)1.2.2计算各单元的整体刚度矩阵由坐标变换式为：K e=TTK eT得(1-6)(10.5000.10500.525K=K=-10.500-0.1051 00.525将（1-2），（1-4）式代入（1-6）式，0-10.500、0.5250-0.1050.5253.50-0.5251.75X 105010.500-0.52500.105-0.5251.750-0.5253.5 /（1-7）将（1-2），（1-5）式代入（1-6）式，得'0.10500.5250.10500.525'010.50010.50一0.52503.50.52501.75K=0.10500.5250.10500.525X 105010.50010.50-0.52501.750.52503.5 J(1-8)1.2.3计算平面钢架的整体刚度矩阵由整体刚度矩阵的形式入下式:(KZKZKZKZ )11121314KZKZKZKZKZ =21222324KZKZKZKZ31323334KZKZKZKZ )41424344(1-9)三个单元的整体刚度矩阵可写成如下形式:单元 10.50010.500 100.1050.52500.1050.525'KZKZ、00.5253.500.5251.751112=KZKZ10.50010.500212200.1050.52500.1050.52500.5251.7500.5253.5 710.50010.500 100.1050.52500.1050.525'KZKZ)00.5253.500.5251.752223=KZKZ710.50010.500323300.1050.52500.1050.52500.5251.7500.5253.5 7 0.10500.5250.10500.525'010.50010.50'KZKZ、0.52503.50.52501.754442二KZKZ70.10500.5250.10500.5252422010.50010.500.52501.750.52503.5 7X 105单元X 105单元X 105(1-10)(1-11)(1-12)由KZ 5 = P得(KZKZKZKZ )(KZ 1KZ 100 I111213141112KZKZKZKZKZ1 KZ 1+ KZ 2 + KZ 3KZ 2KZ 3KZ =21222324 =212222222324KZKZKZKZ0KZ 2KZ 203132333432331 KZKZKZKZ )0KZ 30KZ 3 )V41424344 /4244 7由刚度集成法，将（1-10）式（1-11）式（1-12）式，代入（1-13）式得(1-13)(10.500-10.500.1050.525000.5253.50-10.50021.1050-0.105-0.5250KZ 00.5251.750.52500010.500000000000-0.10500000000.52500000000、-0.1050.525000000-0.5251.7500000000.525-10.500-0.10500.52510.7100-0.1050.5250-10.50010.50-0.5251.75-.52501.75x 1050010.500000-0.105-0.52500.105-0.5250000.5251.7500.5253.50000-0.5250000.1050-0.525-10.50000010.5001.75000-0.52503.5 /(1-14)2计算非节点荷载的等效节点荷载2.1计算局部坐标系下等效节点荷载2.1.1单元有均布荷载如图（2-1）所示9： 6kN/m图（2-1）可知1，2节点有等效节点荷载，取G=-9.6kN/m，均布荷载作用下固端反力计算公式公式（2-1），其中 c=l=10mU =U = 0V0iM0i一-c 2 c 3Gc(2气+履)=2= -Gc - V0i一 c c 2、Gc2 (6 8 F 3)l 12=12(2-1 )Gc 3cM 0 j =商(4 - 31)求得固端反力为:U 01 = U = 0V = 48kNV = 48kNM = 80kN CmM = 80kN Cm 02 则等效节点荷载为(2-2)式的负值.即 (七1'V 工MUIJ02(2-2)-48-80-4880(2-3)2.1.2单元有集中荷载160kN r图（2-2）可知2, 3节点有等效节点荷载，取G=-160kN，由集中荷载作用下固端反力计算公式公式 （2-4），其中 c=5m， d=5mU = U = 0G (/ + 2c)d 2V =0 i13G (l + 2d )c 2V =-0 j13Gcd 2M =0 i12Gc 2 dM = 0 j12(2-4)求的单元的固端反力为：U 02 = U 03 = 0、V02 =的kM = 200kNCm >02V03 =观NM = -200kN Cm03则单元的等效节点荷载为（2-5）式的负值，即(2-5)2.1.3单元有集中力偶如图（2-3）所U0、_02V-8002M-200F=020U003V-804IM 03 J"200 J(2-6)HORN, m个则2,4节点有等效节点荷载，G=80kN m, c=5m，由集中力偶作用下固端反力计算公式（2-7） 为：U = U = 00402匕=四 M = 20kNDm >V =-12kN02M = 20kN Dm02则单元的等效节点荷载为(2-7)式的负值，即（2-7）U（0、_04V-1204M-20F=04=0U00212V"M02)"-20 )（2-8）2.2计算整体坐标系下等效节点荷载单元局部坐标与整体坐标相同，所以单元的整体坐标下等效节点荷载分别为:(U _01V01F=M010U工V02"m )U _0?V02F=M020U_03V03"M03)(0 -48-800-48"80 )(0-80-2000-80"200 )（2-9）（2-10）单元a =90。，所以有F=TtF（2-11）由（1-5）式可求的TT为:0-10000100000001000TT =0000-10（2-12）000100000001J将（2-8）式（2-12）式代入（2-11）式得:I U_04V -04M 二 U 二V 二 M<12 )0一2°-120-20 J(2-13)022.3计算钢架的整体等效节点荷载将各个单元的整体坐标系下等效节点荷载，按节点对应相加。且无节点荷载作用所以由（2-9） 式（2-10）式和（2-13）式可得宕I-1IVIJ_MI-1lv2_ M 2一 U 3一lv3_ M13I M3计算各节点的位移3.1列平衡方程<0 -12 -128 -14( 0 -80 20012 0-20j(2-14)河北工程大学土木工程学院课程设计说明书2011年各节点位移列为一个矩阵如下式：8= C v 0 u v 0 u v 0 u v 0 )(3-1)111222333444乂刚度矩阵KZ，位移列阵。和节点荷载P仔任下列关系：KZ x8 = P(3-2)将（1-14）式，2-15）式和（3-1）代入（3-2）式得：r 10.500-10.50000000o,、/ 、00.105 0.5250-0.105 0.525000000U00.5253.50-0.525L75000000-V-10.50021.10500.525 -10.500-0.10500.5250-0.105 -0.525010.7100-0,105 0.5250-10.50。一 00.5251.750.525010.50-0.525L75-.52501.75-10xd00010.50010.500000M0000-0.105 -0.52500.105 -0.52500000000.525i.7500.5253.5000000-0.1050-0.5250000.1050-0,52500000-10.50000010.50-0000.5250i.75000-0.52503.5：'V6U°，V-48-80-12-128-14(°(3-3)-80200120<20/由图（1-1）知该平面钢架的边界条件为 =V =。= u = v =。= u = v =。= 0，并利111333444用对角线元素置“1”法对方程（2-17）进行处理后得方程式（3-4）：21.105厂 + 0.525耳= 12 x 10 -5、10.78 = -128 x 10-5（3-4）2 .0.525u +10.50 =-140 x 10-5由（3-4）求解后得：U = -0.237206 x 10-5= -11.951447 x10-5 '（3-5）二0 =-13.321462 x 10-52）程序输出各节点位移结果：节点码NJ=纵向位移U=横向位移V=角位移CETA=10.000000000000.000000000000.000000000002-0.00000237206-0.00011951447-0.0001332146230.000000000000.000000000000.0000000000040.000000000000.000000000000.000000000004求单元内力4.1求节点力表达式只要杆单元节点力，就能画出单元的内力图，并可以知单元各截面处的内力，由单元节点位 移引起的节点力的计算公式是Fe = KM(4-1)d(该公式是在局部坐标系中的表达式)又知非节点荷载引起的单元节点力为Fe =(UV M U V M)则单元总节点力的计算公式为Fe = Fe + Fe即（4-2）Fe = Ke5 e + Fe0因为由(3-3)输出的结果为整体坐标系内的位移，故需由式8。= T 6e其变换到局部坐标系中再13（4-3）河北工程大学土木工程学院课程设计说明书代入（4-2）式，即Fe = KeTe 8 e + Fe0式中Te是Q单元的坐标变换矩阵.4.2计算各单元的单元节点力 4.2.1节点位移引起的单元节点力 单元:f U1 1f 10.500 -10.50、00、1V100.1050.525-00.10500.5 2 51M100.5253 . 5-00.5250.75F=1=XdU1-1 0.50010.500-0.2372062V10-0.105-0.5250J 1 05-10.9512572顷1Jk 00.5251.75-00.52 J5k-13.321462JX 1032(2.490663F=K88=8d-5.7388657-17.038049-2.490663X 103(4-1)5.7388657k-40.350607； 单元：f U V 22M 22U 23V 23IM 2 )310.500-10.50000.1050.5250-0.1050.52500.5253.50-0.5251.75-10.50010.5000-0.105-0.52500.105-0.52500.5251.750-0.5253.58=8XF=K8df-0.237206 )-11.951447-13.321462X 103f-2.490663 )-8.2486785-52.8996272.490663X 103(4-2)8.2486785(-29.587068)单元：F=K55顼5d0005=0-10-0.237206-11.951447X 10-5 =-11.9514470.237206X 10-5'U 3 '(10.500-10.500 )'0、4V 300.1050.5250-0.1050.52504M 300.5253.50-0.5251.750F=4=XdU 3-10.50010.500-11.9514472V 30-0.105-0.52500.105-0.5250.2372062(M 3 )(00.5251.750-0.5253.5 )(-13.321462)X 1032-13.321462)(-13.321462)(125.49019-7.0186742-23.4372092-125.490197.0186742(-46.74965 )X 103(4-3)4.2.2非节点荷载引起的固端反力单元:"u 011，0、V4801M8001=U002V4802(M 02 )(-80 )F=0单元:f U )f 0 、02V8002M200F=02=0U003V8003顷03 )"-200) 单元:f U 04 .f 0、V0412M04=20U020V02-123 02)"20 JF=0程序输出各节点力结果:杆单元E=轴力N=剪力Q=弯矩M=1(A)2.4906642.2611362.96194(B)-2.4906653.73887-120.350642(A)-2.4906671.75132147.10033(B)2.4906688.24868-229.587093(A)125.490194.98132-3.43711(B)-125.49019-4.98132-26.749695画内力图5.1节点位移引起的内力图汽 I I II 错误!I I I 丁 I I 2 匚2490N125490N-295咿-m40350k m-23437N m5.2非节点荷载引起的内力图因为非节点作用时无纵向力（轴力），所以为048000N12000N将节点位移引起的内力图与非节点荷载引起的内力图对应叠加的该杆钢架结构的内力图2490N2490N125490N88248N5.3该有限元程序设计框架图程序输出结果：有限元及程序设计课程设计班级：工程力学（2）班姓名：某某某学号：00000节点码NJ=1234纵向位移U= 0.00000000000 -0.00000237206 0.00000000000 0.00000000000横向位移V=0.00000000000-0.000119514470.000000000000.00000000000角位移CETA= 0.00000000000 -0.00013321462 0.00000000000 0.00000000000杆单元E=轴力N=剪力Q=弯矩M=1(A)2.4906642.2611362.96194(B)-2.4906653.73887-120.350642(A)-2.4906671.75132147.10033(B)2.4906688.24868-229.587093(A)125.490194.98132-3.43711(B)-125.49019-4.98132-26.74969主要结论计算结果与程序输出结果基本一致，当然由于手算与机算之间有效数字取值不同产生部分微 小误差。但已经满足工程实际要求。心得体会经过一周多时间的不断努力终于完成了该课程设计。从开始的数据计算，中间的调试程序， 最后写计算书整个过程都是自己亲身亲为。过程很辛苦，结果很美好，也很有成就感，很好的锻 炼了个人独立完成工作的能力，也更加熟悉有限元法的应用，了解了 word更多的功能和进步掌 握了 cad等软件的实际操作。回头看看一路走来的每天感觉很充实，很有意义，更明白坚持就是 胜利。鸣谢特别要感谢许老师在课程设计中对我的帮助。许老师不仅教会我做学问要严谨的作风，更教 会我做人也要严格要求自己，做人比做学问更重要。参考文献1 李景V勇.有限元法.北京：北京邮电大学出版社，19992 吴银柱，吴丽萍.土建工程CAD （第二版）.北京：高等教育出版社，20063 谭浩强.C程序设.北京：清华大学出版社，19994 孟庆元.力学专业英语.哈尔冰：哈尔冰工业大学出版社，2002附件1该钢架有限元程序上机调试(含部分解释)# include<stdio.h># include<math.h># define NE 3# define NJ 4# define NZ 9# define NPJ 0# define NPF 3# define NJ3 12# define DD 9# define E0 2.1000E7# define A0 0.5# define I0 4.16667E-2# define PI 3.141592654int jmNE+13 = 0,0,0,0,1,2,0,2,3,0,4,2;double gcNE+1 = 0.0,10.0,10.0,10.0;double gjNE+1 = 0.0,0.0,0.0,90.0;double mjNE+1 = 0.0,A0,A0,A0;double gxNE+1 = 0.0,I0,I0,I0;int zcNZ+1 = 0,1,2,3,7,8,9,10,11,12;doublepfNPF+15 = 0.0,0.0,0.0,0.0,0.0,0.0,-9.6,10.0,1.0,1.0,0.0,-160,5.0,2.0,2.0, 0.0,80,5.0,3.0,3.0;double kzNJ3+1 NJ 3+1,pNJ 3+1;double pe7,f7,f07,t77;double ke77,kd77;void jdugd(int);void zb(int);void gdnl(int);void dugd(int);void main(void)int i,j,k,e,dh,h,ii,jj,hz,a1,b1,m,n,l,dl,zl,z,j0;double c1,wy7;int I虬 IN,jn;FILE*jg;/*声明指针文件*/jg=fopen（"某某某.txt，w+）；/*以追加的方式建立 某某某.txt，默认位置在你程序的目录下面*/for(i=1;i<=NPF;i+)hz=i;gdnl(hz);e=(int)pfhz3;zb(e);for(j=1;j<=6;j+)pej=0.0;for(k=1;k<=6;k+)pej=pej-tkj*f0k;a1=jme1;b1=jme