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

1、有限元编程的 c+ 实现算例 n.cpp_.htm1. #include2. #include3.4.5. #define ne3/ 单元数6. #define nj4/ 节点数7. #define nz6/ 支撑数8. #define npj0/ 节点载荷数9. #define npf1/ 非节点载荷数10. #define nj312/ 节点位移总数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整体坐标下的单元刚度矩阵33. /*kd局部坐标下的单位刚度矩阵34. /*t坐标变换矩阵2.63.64./* 这是函数声明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);/ 求单元号码 for(j=1;j=6;j+)/ 求坐标变换矩阵 T

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

6、dugd(e);/ 求整体坐标系中的单元刚度矩阵 ke for(i=1;i=2;i+)/ 对行码循环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+)/ 按支撑循环z=zci;/ 支撑对应的位移数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./* 根据 E的值输出相应 E单元的 N,Q,M(A,B) 的结果*printf(En);187./* 计算轴力 N,剪力 Q,弯矩 M*188.for(e=1;e=ne;e+)/ 按单元循环190.189.jdugd(e);000015.216./ 求

10、局部单元刚度矩阵 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+)fi=0.0;for(j=1;j=6;j+)for(k=1;k0)for(i=1;i=npf;i+)/ 按非节点载荷数循环if(pfi3=e)/ 找到荷载所在的单元217.hz=i;218.gdnl(hz);25./ 求固端反力for(j=1;j233. / 功能：将非节点载荷下的杆端力计算

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;/ 计算 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+)/ 计算 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;计算 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 左下角元素按对称原则送到右下角 kdji=kdij;/*/*void dugd(int e)int i,k,j,m;jdug