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

1、有限元编程的c+实现算例   1. #include<stdio.h>      2. #include<math.h>      3.       4.       5. #define ne 3               

2、;                                    /单元数      6. #define nj 4       &

3、#160;                                           /节点数      7. #define n

4、z 6                                                 

5、60; /支撑数      8. #define npj 0                                         

6、0;        /节点载荷数      9. #define npf 1                                  

7、60;               /非节点载荷数     10. #define nj3 12                            

8、;                     /节点位移总数     11. #define dd 6                      

9、60;                            /半带宽     12. #define e0 2.1E8               

10、                               /弹性模量     13. #define a0 0.008            &#

11、160;                                  /截面积     14. #define i0 1.22E-4        

12、60;                                    /单元惯性距      15. #define pi 3.141592654    

13、0;                                16.      17.      18. int jmne+13=0,0,0,0,1,2,0,2,3,0,4,3;   

14、          /*gghjghg*/                  19. double gcne+1=0.0,1.0,2.0,1.0;                 

15、                 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;    

16、; 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;     26. double kznj3+1dd+1,pnj3+1;     27. double pe7,f7,f07,t77;     28. double ke77,kd77;     29.      30. 

17、60;    31. /*kz整体刚度矩阵     32. /*ke整体坐标下的单元刚度矩阵     33. /*kd局部坐标下的单位刚度矩阵     34. /*t坐标变换矩阵     35.      36. /*这是函数声明     37. void jdugd(int);     38. void zb(int);    

18、; 39. void gdnl(int);     40. void dugd(int);     41.      42.      43. /*主程序开始     44. void main()     45.      46.  int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;   &#

19、160; 47.  double cl,wy7;     48.  int im,in,jn;     49.      50. /*     51. /<功能：形成矩阵P>     52. /*     53.      54.  if(npj>0)     55.   

20、;    56.      for(i=1;i<=npj;i+)     57.                                   &#

21、160;                         /把节点载荷送入P     58.          j=pji2;         

22、60;                    59.          pj=pji1;     60.           61.       62.  if(npf&

23、gt;0)     63.       64.      for(i=1;i<=npf;i+)     65.                            

24、0;                               /求固端反力F0     66.          hz=i;     67.

25、60;         gdnl(hz);     68.          e=(int)pfhz3;     69.          zb(e);          

26、0;                                  /求单元号码     70.          for(j=1;j<=6

27、;j+)                                  /求坐标变换矩阵T     71.           

28、60;  72.              pej=0.0;     73.              for(k=1;k<=6;k+)             &

29、#160;                /求等效节点载荷     74.                    75.           

30、;       pej=pej-tkj*f0k;     76.                   77.               78.       &

31、#160;  al=jme1;     79.          bl=jme2;     80.          p3*al-2=p3*al-2+pe1;                

32、         /将等效节点载荷送到P中     81.          p3*al-1=p3*al-1+pe2;     82.          p3*al=p3*al+pe3;     83.    

33、      p3*bl-2=p3*bl-2+pe4;     84.          p3*bl-1=p3*bl-1+pe5;     85.          p3*bl=p3*bl+pe6;     86.     

34、60;    87.       88.      89.      90.  /*     91.  /<功能：生成整体刚度矩阵kz>     92.  for(e=1;e<=ne;e+)             

35、60;                          /按单元循环     93.        94.      dugd(e);       

36、                                       /求整体坐标系中的单元刚度矩阵ke     95.     

37、for(i=1;i<=2;i+)                                     /对行码循环     96.      

38、0;    97.          for(ii=1;ii<=3;ii+)     98.                99.              h=3*(i-1)+ii; 

39、;                                /元素在ke中的行码    100.             

40、dh=3*(jmei-1)+ii;                         /该元素在KZ中的行码    101.              for(j=1;j<=2;j+) 

41、   102.                  103.                  for(jj=1;jj<=3;jj+)         &#

42、160;            /对列码循环    104.                      105.            &#

43、160;         l=3*(j-1)+jj;                         /元素在ke中的列码    106.         &#

44、160;            zl=3*(jmej-1)+jj;                 /该元素在KZ中的行码    107.             

45、60;        dl=zl-dh+1;                           /该元素在KZ*中的行码    108.        

46、60;             if(dl>0)    109.                          kzdhdl=kzdhdl+kehl;   /刚度集成&#

47、160;   110.                      111.                  112.        

48、60;     113.          114.      115.     116. /*引入边界条件*    117.   for(i=1;i<=nz;i+)               

49、0;                       /按支撑循环    118.      119.      z=zci;           

50、60;                                  /支撑对应的位移数    120.      kzzl=1.0;    &#

51、160;                                    /第一列置1    121.      for(j=2;j<=dd;j+) 

52、;   122.          123.          kzzj=0.0;                           &#

53、160;         /行置0    124.          125.      if(z!=1)    126.          127.        &#

54、160; if(z>dd)    128.              j0=dd;    129.          else if(z<=dd)    130.         

55、0;    j0=z;                                         /列（45度斜线）置0    131

56、.          for(j=2;j<=j0;j+)    132.              kzz-j+1j=0.0;    133.          134.      p

57、z=0.0;                                             /P置0    135.&#

58、160;     136.     137.     138.     139.     140.  for(k=1;k<=nj3-1;k+)    141.      142.      if(nj3>k+dd-1)    &

59、#160;                                   /求最大行码    143.          im=k+d

60、d-1;    144.      else if(nj3<=k+dd-1)    145.          im=nj3;    146.      in=k+1;    147.      for(i=in;i<=im;i

61、+)    148.          149.          l=i-k+1;    150.          cl=kzkl/kzk1;           &

62、#160;                 /修改KZ    151.          jn=dd-l+1;    152.          for(j=1;j<=jn;j+) 

63、  153.              154.              m=j+i-k;    155.              kzij=kzij-cl*kzkm

64、;    156.              157.          pi=pi-cl*pk;                      

65、;         /修改P    158.          159.      160.     161.     162.     163.     164.  pnj3=pnj3/kznj31

66、;                                /求最后一个位移分量    165.  for(i=nj3-1;i>=1;i-)    166.      

67、167.      if(dd>nj3-i+1)    168.          j0=nj3-i+1;    169.      else j0=dd;               

68、0;                          /求最大列码j0    170.      for(j=2;j<=j0;j+)    171.        

69、0; 172.          h=j+i-1;    173.          pi=pi-kzij*ph;    174.          175.      pi=pi/kzi1;  &#

70、160;                               /求其他位移分量    176.      177.      printf("n");&

71、#160;   178.      printf("_n");    179.      printf("NJ            U              

72、60;  V                CETA       n");     /输出位移    180.      for(i=1;i<=nj;i+)    181.   

73、60;      182.          printf(" %-5d %14.11f     %14.11f     %14.11fn",i,p3*i-2,p3*i-1,p3*i);    183.          184.  

74、60;   printf("_n");    185.      /*根据E的值输出相应E单元的N,Q,M(A,B)的结果*    186.      printf("E             N      

75、            Q                 M     n");    187.      /*计算轴力N,剪力Q,弯矩M*    188. &

76、#160;    for(e=1;e<=ne;e+)                                  /按单元循环    189.    

77、0;     190.          jdugd(e);                                   

78、    /求局部单元刚度矩阵kd    191.          zb(e);                               

79、;           /求坐标变换矩阵T    192.          for(i=1;i<=2;i+)    193.              194.     &

80、#160;        for(ii=1;ii<=3;ii+)    195.                  196.                  h=

81、3*(i-1)+ii;    197.                  dh=3*(jmei-1)+ii;                   /给出整体坐标下单元节点位移    198.

82、                  wyh=pdh;    199.                  200.            &#

83、160; 201.          for(i=1;i<=6;i+)    202.              203.              fi=0.0;    204.&#

84、160;             for(j=1;j<=6;j+)    205.                  206.             &

85、#160;    for(k=1;k<=6;k+)                      /求由节点位移引起的单元节点力    207.               &#

86、160;      208.                      fi=fi+kdij*tjk*wyk;    209.                &

87、#160;     210.                  211.              212.          if(npf>0)    213.&

88、#160;             214.              for(i=1;i<=npf;i+)                    

89、60;   /按非节点载荷数循环    215.                  if(pfi3=e)                       

90、 /找到荷载所在的单元    216.                      217.                      hz=i;  &

91、#160; 218.                      gdnl(hz);                         

92、; /求固端反力    219.                      for(j=1;j<=6;j+)                  /将固端反力累加   

93、; 220.                          221.                       &

94、#160;  fj=fj+f0j;    222.                          223.                  &#

95、160;   224.              225.          printf("%-3d(A)  %9.5f           %9.5f        

96、 %9.5fn",e,f1,f2,f3);         /输出单元A（i）端内力    226.          printf("   (B)  %9.5f           %9.5f    

97、0;    %9.5fn",f4,f5,f6);            /输出单元B（i）端内力    227.          228.       return;    229.     230. /*主程序结束*

98、    231.     232. /*    233. /<功能：将非节点载荷下的杆端力计算出来存入f0>    234. /*    235.     236. void gdnl(int hz)    237.     238.     int ind,e;   

99、; 239.     double g,c,l0,d;    240.     241.     242.     g=pfhz1;                      &#

100、160;                /载荷值    243.     c=pfhz2;                        &

101、#160;              /载荷位置    244.     e=(int)pfhz3;                         

102、         /作用单元    245.     ind=(int)pfhz4;                              

103、60; /载荷类型    246.     l0=gce;                                       &#

104、160;  /杆长    247.     d=l0-c;    248.     if(ind=1)    249.         250.         f01=0.0;    251.    

105、;     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.     &#

106、160;   f04=0.0;    254.         f05=-g*c-f02;    255.         f06=(g*c*c*c)*(4-3*c/l0)/(12*l0);    256.         257.  &#

107、160;  else    258.         259.         if(ind=2)                          

108、          /横向集中力的固端反力    260.             261.             f01=0.0;    262.      

109、;       f02=(-(g*d*d)*(l0+2*c)/(l0*l0*l0);    263.             f03=-(g*c*d*d)/(l0*l0);    264.             f04=0.0; 

110、;   265.             f05=(-g*c*c*(l0+2*d)/(l0*l0*l0);    266.             f06=(g*c*c*d)/(l0*l0);    267.      

111、       268.         else    269.             270.             f01=-(g*d/l0);     

112、;                    /纵向集中力的固端反力    271.             f02=0.0;    272.       

113、0;     f03=0.0;    273.             f04=-g*c/l0;    274.             f05=0.0;    275.     

114、        f06=0.0;    276.             277.         278.     279.     280. /*    281. /<功能：构成坐标变换矩阵> 

115、   282. /*    283. void zb(int e)    284.     285.     double ceta,co,si;    286.     int i,j;    287.     ceta=(gje*pi)/180;     

116、                        /角度变弧度    288.     co=cos(ceta);    289.     si=sin(ceta);    290. 

117、0;   t11=co;                                      /计算T右上角元素    291.    

118、; t12=si;    292.     t21=-si;    293.     t22=co;    294.     t33=1.0;    295.     for(i=1;i<=3;i+)    296.       &#

119、160; 297.         for(j=1;j<=3;j+)                            /计算T的左下角元素    298.   

120、0;         299.             ti+3j+3=tij;    300.             301.         302.     30

121、3.     304.     305.     306. /*    307. /<计算局部坐标下单元刚度矩阵kd>    308. /*    309. void jdugd(int e)    310.        311.     dou

122、ble A0,l0,j0;    312.     int i;    313.     int j;    314.         315.         316.     A0=mje;   

123、                                    /面积    317.     l0=gce;     &

124、#160;                                 /杆长    318.     j0=gxe;       &#

125、160;                               /惯性钜    319.     320.         321.     for(i=0;i<=6;i+)    322.         for(j=0;j<=6;j+)         &#