有限元編程的c++實(shí)現(xiàn)算例

1、有限元編程的C+實(shí)現(xiàn)算例n.cpp_.htm1. #include<stdio.h>2. #include<math.h>3.4.5. #definene3/單元數(shù)6. #definenj4/節(jié)點(diǎn)數(shù)7. #definenz6/支撐數(shù)8. #definenpj0/節(jié)點(diǎn)載荷數(shù)9. #definenpf1/非節(jié)點(diǎn)載荷數(shù)10. #definenj312/節(jié)點(diǎn)位移總數(shù)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一整體坐標(biāo)下的單元?jiǎng)偠染仃?3. /*kd一局部坐標(biāo)下的單位剛度矩陣34. /*t加一坐標(biāo)變換矩陣35.36. /*這是函數(shù)聲明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. /把節(jié)點(diǎn)載荷送入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)；/求單元號(hào)碼for(j=1;j<=6;j+)/求坐標(biāo)變換矩陣Tpej=0.0;for(k=1;k<=6;k+)/求等效節(jié)點(diǎn)載荷pej=pej-tkj*f0k;al=jme1;bl=jme2;p3*al-2=p3*al-2+pe1;/將等效節(jié)點(diǎn)載荷送到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+)/按單元循環(huán)93. 94. dugd(e);/求整體坐標(biāo)系中的單元?jiǎng)偠染仃噆e95. for(i=1;i<=2;i+)/對(duì)行碼循環(huán)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+)/對(duì)列碼循環(huán)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+)/按支撐循環(huán)z=zci；/支撐對(duì)應(yīng)的位移數(shù)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;/求最后一個(gè)位移分量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./*根據(jù)E的值輸出相應(yīng)E單元的N,Q,M(A,B)的結(jié)果*printf("ENM'n");/*計(jì)算軸力N,剪力Q,

11、彎矩M*for(e=1;e<=ne;e+)按單元循環(huán)%14.11fn",i,p3190.jdugd(e)；000015./求局部單元?jiǎng)偠染仃噆dzb(e)；/求坐標(biāo)變換矩陣Tfor(i=1;i<=2;i+)for(ii=1;ii<=3;ii+)h=3*(i-1)+ii;dh=3*(jmei-1)+ii;/給出整體坐標(biāo)下單元節(jié)點(diǎn)位移wyh=pdh;for(i=1;i<=6;i+)

12、fi=0.0;for(j=1;j<=6;j+)for(k=1;k<=6;k+)/求由節(jié)點(diǎn)位移引起的單元節(jié)點(diǎn)力fi=fi+kdij*t皿k*wyk;if(npf>0)for(i=1;i<=npf;i+)/按非節(jié)點(diǎn)載荷數(shù)循環(huán)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)端內(nèi)力226.printf("(B)%9.5f%9.5f%9.5fn",f4,f5,f6);/輸出單元B(i)端內(nèi)力227. 228. return;229. 230. /*主程序結(jié)束*231.232. /*233./<功能：將非節(jié)點(diǎn)載荷下的桿端力計(jì)算出來存入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./功能：構(gòu)成坐標(biāo)變換矩陣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;/計(jì)算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./計(jì)算T的左下角元素298.299. ti+3j+3=tij;300. 301. 302. 303.304.305.306. /*307./計(jì)算

17、局部坐標(biāo)下單元?jiǎng)偠染仃噆d叩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;/計(jì)算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左下角元素按對(duì)稱原則送到右下角341.kd皿卜kdi皿kd22=12*e0*j0/pow(l0,3);342. 343.344.345. /*346. /<計(jì)算整體坐標(biāo)下單元?jiǎng)偠染仃噆e>347. /*348. void dugd(int e)349. 350.int i,k,j,m;351