有限元编程的c++实现算例

1、有限元编程的C+实现算例n.cpp_.htm1. #include<stdio.h>2. #include<math.h>3.4.5. #definene3/单元数6. #definenj4/节点数7. #definenz6/支撑数8. #definenpj0/节点载荷数9. #definenpf1/非节点载荷数10. #definenj312/节点位移总数11. #definedd6/半带宽12. #defineeO2.1E8/弹性模量13. #definea00.008/截面积14. #definei01.22E-4/单元惯性距15. #definepi3.1415

2、9265416.17.18. int/*gghjghg*/jmne+13=0,0,0,0,1,2,023,0,4,3;19. doublegcne+1=0.0,1.0,2.0,1.0;20. doublegjne+1=0.0,90.0,0.0,90.0;21. doublemjne+1=0.0,a0,a0,a0;22. doublegxne+1=0.0,i0,i0,i0;23. intzcnz+1=0,1,2,3,10,11,12;24. doublepjnpj+13=0.0,0.0,0.0;25. doublepfnpf+15=0,0,0,0,0,0,-20,1.0,2.0,2.0;26.

3、doublekznj3+1dd+1,pnj3+1;27. doublepe7,f7,f07,t77;28. doubleke77,kd77;29.30.31. /*kz加一整体刚度矩阵32. /*kem一整体坐标下的单元刚度矩阵33. /*kd一局部坐标下的单位刚度矩阵34. /*t加一坐标变换矩阵35.36. /*这是函数声明37. voidjdugd(int);38. voidzb(int);39. voidgdnl(int);40. voiddugd(int);41.42.43. /*主程序开始44. voidmain()45. 46. inti,j,k,e,dh,h,ii,jj,hz,

4、al,bl,m,l,dl,zl,z,j0;47. doublecl,wy7;48. intim,in,jn;49.50. /*51. /<功能：形成矩阵P>52. /*53.54. if(npj>0)55. 56. for(i=1;i<=npj;i+)57. /把节点载荷送入P58. j=pj皿2;59. pj=pji1;60. 61. 62. if(npf>0)63. 64. for(i=1;i<=npf;i+)65./求固端反力F04.85.86

5、.0.91.hz=i;gdnl(hz)；e=(int)pfhz3;zb(e)；/求单元号码for(j=1;j<=6;j+)/求坐标变换矩阵Tpej=0.0;for(k=1;k<=6;k+)/求等效节点载荷pej=pej-tkj*f0k;al=jme1;bl=jme2;p3*al-2=p3*al-2+pe1;/将等效节点载荷送到P中p3*al-1=p3*al-1+pe2;p3*al=p3*al+pe3;p3*bl-2=p3*bl-2+pe4;p3*bl-1=p3*bl-1+pe5;p3*bl=p3*bl+pe6;/*/<功能：生成整体刚度矩阵kz>92

6、. for(e=1;e<=ne;e+)/按单元循环93. 94. dugd(e);/求整体坐标系中的单元刚度矩阵ke95. for(i=1;i<=2;i+)/对行码循环96. 97. for(ii=1;ii<=3;ii+)98. 99. h=3*(i-1)+ii;/元素在ke中的行码100. dh=3*(jmei-1)+ii;/该元素在KZ中的行码101. for(j=1;j<=2;j+)102. 103. for(jj=1;jj<=3;jj+)/对列码循环104. 105. l=3*(j-1)+jj;/元素在ke中的列码106. zl=3*(jmej-1)+jj

7、;/该元素在KZ中的行码107. dl=zl-dh+1;/该元素在KZ*中的行码108. if(dl>0)109. kzdhdl=kzdhdl+kehl/刚度集成110. 111. 112. 113.114. 115.116. /*引入边界条件*38.for(i=1;i<=nz;i+)/按支撑循环z=zci；/支撑对应的位移数kzzl=1.0;/第一列置1for(j=2;j<=dd;j+)kzzj=0.0；/

8、行置0if(z!=1)if(z>dd)j0=dd;elseif(z<=dd)j0=z;/歹I(45度斜线)置0for(j=2;j<=j0;j+)kzz-j+1j=0.0;pz=0.0;/P置0139.140. for(k=1;k<=nj3-1;k+)141. 142. if(nj3>k+dd-1)/求最大行码143.im=k+dd-1;144.elseif(nj3<=k+dd-1)145.im=nj3;146.in=k+1;147.for(i=in;i<=im;i+)148.149.l=i-k+1;150.cl=kzkl/kzk1;/修改KZ151.j

9、n=dd-l+1;152.for(j=1;j<=jn;j+)153.154.m=j+i-k;155.kzij=kzij-cl*kzkm;156.157.Pi=Pi-cl*Pk;/修改P64. pnj3=pnj3/kznj31;/求最后一个位移分量165. for(i=nj3-1;i>=1;i-)166. 167. if(dd>nj3-i+1)168.j0=nj3-i+1;169.elsej0=dd;/求最大列码j076. 177.178.n");179.180.181

10、.182.*i-2,p3*i-1,p3*i);183.184.n");for(j=2;j<=j0;j+)h=j+i-1;pi=pi-kzij*ph;pi=pi/kzi1;/求其他位移分量printf("n");printf("NJCETAfor(i=1;i<=nj;i+)printf("%-5d%14.11fUn");/输出位移%14.11f185.186.Q187.188.189./*根据E的值输出相应E单元的N,Q,M(A,B)的结果*printf("ENM'n");/*计算轴力N,剪力Q,

11、弯矩M*for(e=1;e<=ne;e+)按单元循环%14.11fn",i,p3190.jdugd(e)；000015./求局部单元刚度矩阵kdzb(e)；/求坐标变换矩阵Tfor(i=1;i<=2;i+)for(ii=1;ii<=3;ii+)h=3*(i-1)+ii;dh=3*(jmei-1)+ii;/给出整体坐标下单元节点位移wyh=pdh;for(i=1;i<=6;i+)

12、fi=0.0;for(j=1;j<=6;j+)for(k=1;k<=6;k+)/求由节点位移引起的单元节点力fi=fi+kdij*t皿k*wyk;if(npf>0)for(i=1;i<=npf;i+)/按非节点载荷数循环if(pfi3=e)/找到荷载所在的单元216.217.hz=i;218./求固端反力gdnl(hz);219./将固端反力累加for(j=1;j<=6;j+)220.221.fj=fj+f0j;25.printf("%-3d(A)%9.5f%9.5f%9.5fn",e,f1,f2,f3);/输出单元A

13、(i)端内力226.printf("(B)%9.5f%9.5f%9.5fn",f4,f5,f6);/输出单元B(i)端内力227. 228. return;229. 230. /*主程序结束*231.232. /*233./<功能：将非节点载荷下的杆端力计算出来存入f0>234. /*235.236. voidgdnl(inthz)237. 238. intind,e;239. doubleg,c,l0,d;42.g=pfhz1;/载荷值/载荷位置244. e=(int)pfhz3;/作用单元245. ind=(int)pfhz4;/载

14、荷类型246.l0=gce;/技247.d=l0-c;248.if(ind=1)249.250.f01=0.0;251.f02=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)/2/均布载荷的固端反力252.f03=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0)/12;253.f04=0.0;254.f05=-g*c-f02;255.f06=(g*c*c*c)*(4-3*c/l0)/(12*l0);256.257.else258.259.if(ind=2)/横向集中力的固端反力260.261.f01=0.0;262.f02=(-(g*d*d)

15、*(l0+2*c)/(l0*l0*l0);263.f03=-(g*c*d*d)/(l0*l0);264.f04=0.0;265.f05=(-g*c*c*(l0+2*d)/(l0*l0*l0);266.f06=(g*c*c*d)/(l0*l0);267.268.else269.270. f01=-(g*d/l0);/纵向集中力的固端反力271. f02=0.0;272. f03=0.0;273. f04=-g*c/l0;274. f05=0.0;275. f06=0.0;276. 277. 278. 279.280./*281./功能：构成坐标变换矩阵282./*283. voidzb(inte

16、)284. 285. doubleceta,co,si;286. inti,j;287. ceta=(gje*pi)/180;/角度变弧度288. co=cos(ceta);289. si=sin(ceta);290. t11=co;/计算T右上角元素291. t12=si;292. t21=-si;293. t22=co;294. t33=1.0;295. for(i=1;i<=3;i+)296. for(j=1;j<=3;j+)297./计算T的左下角元素298.299. ti+3j+3=tij;300. 301. 302. 303.304.305.306. /*307./计算

17、局部坐标下单元刚度矩阵kd叩308. /*309. voidjdugd(inte)310. 311. doubleA0,l0,j0;312. inti;313. intj;314.315.316. A0=mje;/面积317. l0=gce;/杆长318. j0=gxe;/惯性银319.320.321. for(i=0;i<=6;i+)322. for(j=0;j<=6;j+)/kd清0323. kdij=0.0;324.325. kd11=e0*A0/l0;326.327.328.kd33=4*e0*j0/l0;329.kd41=-kd11;330.kd44=kd11;331.k

18、d52=-kd22;/计算kd左下角各元素332.kd53=-kd32;333.kd55=kd22;334.kd62=kd32;335.kd63=2*e0*j0/l0;336.kd65=-kd32;337.kd66=kd33;338.339.for(i=1;i<=6;i+)340.for(j=1；j<=i；j+)/将kd左下角元素按对称原则送到右下角341.kd皿卜kdi皿kd22=12*e0*j0/pow(l0,3);342. 343.344.345. /*346. /<计算整体坐标下单元刚度矩阵ke>347. /*348. void dugd(int e)349. 350.int i,k,j,m;351