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

1、有限元編程的 c+ 實(shí)現(xiàn)算例 n.cpp_.htm1. #include2. #include3.4.5. #define ne3/ 單元數(shù)6. #define nj4/ 節(jié)點(diǎn)數(shù)7. #define nz6/ 支撐數(shù)8. #define npj0/ 節(jié)點(diǎn)載荷數(shù)9. #define npf1/ 非節(jié)點(diǎn)載荷數(shù)10. #define nj312/ 節(jié)點(diǎn)位移總數(shù)11. #define dd6/ 半帶寬12. #define e02.1E8/ 彈性模量13. #define a00.008/ 截面積14. #define i01.22E-4/ 單元慣性距15. #define pi3.141592654

2、16.17.18. int/*gghjghg*/jmne+13=0,0,0,0,1,2,0,2,3,0,4,3;19. double gcne+1=0.0,1.0,2.0,1.0;20. double gjne+1=0.0,90.0,0.0,90.0;21. double mjne+1=0.0,a0,a0,a0;22. double gxne+1=0.0,i0,i0,i0;23. int zcnz+1=0,1,2,3,10,11,12;24. double pjnpj+13=0.0,0.0,0.0;25. double pfnpf+15=0,0,0,0,0,0,-20,1.0,2.0,2.0;

3、26. double kznj3+1dd+1,pnj3+1;27. double pe7,f7,f07,t77;28. double ke77,kd77;29.30.31. /*kz整體剛度矩陣32. /*ke整體坐標(biāo)下的單元?jiǎng)偠染仃?3. /*kd局部坐標(biāo)下的單位剛度矩陣34. /*t坐標(biāo)變換矩陣2.63.64./* 這是函數(shù)聲明void jdugd(int);void zb(int);void gdnl(int); void d

4、ugd(int);/* 主程序開始void main()int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;double cl,wy7;int im,in,jn;/*/*if(npj0)for(i=1;i0)for(i=1;i=npf;i+)65./ 求固端反力 F00.91.hz=i;gdnl(hz);e=(int)pfhz3;zb(e);/ 求單元號(hào)碼 for(j=1;j=6;j+)/ 求坐標(biāo)變換矩陣 T

5、pej=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;/*/0007.108.109./13.for(e=1;e=ne;e+)/ 按單元循環(huán)

6、dugd(e);/ 求整體坐標(biāo)系中的單元?jiǎng)偠染仃?ke for(i=1;i=2;i+)/ 對(duì)行碼循環(huán)for(ii=1;ii=3;ii+) h=3*(i-1)+ii;/ 元素在 ke 中的行碼 dh=3*(jmei-1)+ii;/ 該元素在 KZ中的行碼for(j=1;j=2;j+)for(jj=1;jj0)kzdhdl=kzdhdl+kehl剛度集成114. 115.116. /* 引入邊界條件 *38.139.for(i=1

7、;i=nz;i+)/ 按支撐循環(huán)z=zci;/ 支撐對(duì)應(yīng)的位移數(shù)kzzl=1.0;/ 第一列置 1for(j=2;jdd)j0=dd;else if(z=dd)j0=z;/ 列(45 度斜線)置 0 for(j=2;j=j0;j+)kzz-j+1j=0.0;pz=0.0;/P 置 067.for(k=1;kk+dd-1)/ 求最大行碼im=k+dd-1;else if(nj3=k

8、+dd-1)im=nj3;in=k+1;for(i=in;i=im;i+)l=i-k+1;cl=kzkl/kzk1;/ 修改 KZjn=dd-l+1;for(j=1;j=1;i-)if(ddnj3-i+1)168.j0=nj3-i+1;else169.j0=dd;/ 求最大列碼 j0170.for(j=2;j=j0;j+)172.171.173.h=j+i-1;pi=pi-kzij*ph;174.175.pi=pi/kzi1;/ 求其他位移分量177.176.printf(n);178.179.n);printf(NJ180.CETAn);/ 輸出位移for(i=1;i=nj;i+)181.1

9、82.printf( %-5d %14.11f%14.11f%14.11fn,i,p3*i-2,p3*i-1,p3*i);183.184.n);185.186./* 根據(jù) E的值輸出相應(yīng) E單元的 N,Q,M(A,B) 的結(jié)果*printf(En);187./* 計(jì)算軸力 N,剪力 Q,彎矩 M*188.for(e=1;e=ne;e+)/ 按單元循環(huán)190.189.jdugd(e);000015.216./ 求

10、局部單元?jiǎng)偠染仃?kdzb(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+)fi=0.0;for(j=1;j=6;j+)for(k=1;k0)for(i=1;i=npf;i+)/ 按非節(jié)點(diǎn)載荷數(shù)循環(huán)if(pfi3=e)/ 找到荷載所在的單元217.hz=i;218.gdnl(hz);25./ 求固端反力for(j=1;j233. / 功能：將非節(jié)點(diǎn)載荷下的桿端力計(jì)算

11、出來存入234. /* 235.236. void gdnl(int hz)238.237. int ind,e;239.double g,c,l0,d;240.241.242.g=pfhz1;/ 載荷值49.250.251./55.256.257.258.259.260.261.262.263.264.265.266.267.268.269./ 載荷位置e=(int)pfhz3;/ 作用單元ind=(int)pfhz4;/ 載荷類型l0=gce;/ 桿長d=l0-c;if(ind=1)f01=0.0;f02=-(g*

12、c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)/2;均布載荷的固端反力f03=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0)/12;f04=0.0;f05=-g*c-f02;f06=(g*c*c*c)*(4-3*c/l0)/(12*l0);elseif(ind=2)/ 橫向集中力的固端反力f01=0.0;f02=(-(g*d*d)*(l0+2*c)/(l0*l0*l0);f03=-(g*c*d*d)/(l0*l0);f04=0.0;f05=(-g*c*c*(l0+2*d)/(l0*l0*l0);f06=(g*c*c*d)/(l0*l0);else2

13、70.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. void zb(int e)284. 285.double ceta,co,si;286.int i,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.

14、t22=co;294.t33=1.0;295.for(i=1;i=3;i+)296.297.for(j=1;j=3;j+)/ 計(jì)算 T 的左下角元素298.299.300.301.302.303.304.305.306.307.308.309.310.311.312.313.314.315.316.317.318.319.320.321.322.323.324.325.ti+3j+3=tij;/*/*void jdugd(int e)double A0,l0,j0;int i;int j;A0=mje;/ 面積l0=gce;/ 桿長j0=gxe;/ 慣性鉅for(i=0;i=6;i+)for(

15、j=0;j=6;j+)/kd 清 0kdij=0.0;kd11=e0*A0/l0;326.327.328.329.330.331./332.333.334.335.336.337.338.339.340./341.342.343.344.345.346.347.348.349.350.351.352.353.kd22=12*e0*j0/pow(l0,3);kd32=6*e0*j0/pow(l0,2);kd33=4*e0*j0/l0;kd41=-kd11;kd44=kd11;kd52=-kd22;計(jì)算 kd 左下角各元素kd53=-kd32;kd55=kd22;kd62=kd32;kd63=2*e0*j0/l0;kd65=-kd32;kd66=kd33;for(i=1;i=6;i+)for(j=1;j=i;j+)將 kd 左下角元素按對(duì)稱原則送到右下角 kdji=kdij;/*/*void dugd(int e)int i,k,j,m;jdug