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

1、有限元编程的c+实现算例有限元编程的c+实现算例read.pudn./downloads76/doc/fileformat/290377/ganjian.cpp_.htm 1. #include 2. #include 3. 4. 5. #define ne 3 /单元数 6. #define nj 4 /节点数 7. #define nz 6 /支撑数 8. #define npj 0 /节点载荷数 9. #define npf 1 /非节点载荷数 10. #define nj3 12 /节点位移总数 11. #define dd 6 /半带宽 12. #define e0 2.1E8 /弹

2、性模量 13. #define a0 0.008 /截面积 14. #define i0 1.22E-4 /单元惯性距 15. #define pi 3.141592654 16. 17. 18. int jmne+13=0,0,0,0,1,2,0,2,3,0,4,3; /*gghjghg*/ 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

3、,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; 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坐标变换矩阵 35. 36. /*这是函数声明 37. void jdugd(int); 38.

4、void zb(int); 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; 47. double cl,wy7; 48. int im,in,jn; 49. 50. /* 51. / 52. /* 53. 54. if(npj0) 55. 56. for(i=1;i0) 63. 64. for(i=1;i=npf;i+) 65. /求固端反力F0 66. hz=i; 67. gdn

5、l(hz); 68. e=(int)pfhz3; 69. zb(e); /求单元 70. for(j=1;j=6;j+) /求坐标变换矩阵T 71. 72. pej=0.0; 73. for(k=1;k=6;k+) /求等效节点载荷 74. 75. pej=pej-tkj*f0k; 76. 77. 78. al=jme1; 79. bl=jme2; 80. p3*al-2=p3*al-2+pe1; /将等效节点载荷送到P中 81. p3*al-1=p3*al-1+pe2; 82. p3*al=p3*al+pe3; 83. p3*bl-2=p3*bl-2+pe4; 84. p3*bl-1=p3*

6、bl-1+pe5; 85. p3*bl=p3*bl+pe6; 86. 87. 88. 89. 90. /* 91. / 92. for(e=1;e=ne;e+) /按单元循环 93. 94. dugd(e); /求整体坐标系中的单元刚度矩阵ke 95. 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;jj0) 109. kzdhd

7、l=kzdhdl+kehl; /刚度集成 110. 111. 112. 113. 114. 115. 116. /*引入边界条件* 117. for(i=1;i=nz;i+) /按支撑循环 118. 119. z=zci; /支撑对应的位移数 120. kzzl=1.0; /第一列置1 121. for(j=2;jdd) 128. j0=dd; 129. else if(z=dd) 130. j0=z; /列（45度斜线）置0 131. for(j=2;j=j0;j+) 132. kzz-j+1j=0.0; 133. 134. pz=0.0; /P置0 135. 136. 137. 138.

8、139. 140. for(k=1;kk+dd-1) /求最大行码 143. im=k+dd-1; 144. else if(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; /修改KZ 151. jn=dd-l+1; 152. for(j=1;j=1;i-) 166. 167. if(ddnj3-i+1) 168. j0=nj3-i+1; 169. else j0=dd; /求最大列码j0 170. for(j=2;j=j0;j+) 171. 1

9、72. h=j+i-1; 173. pi=pi-kzij*ph; 174. 175. pi=pi/kzi1; /求其他位移分量 176. 177. printf(n); 178. printf(_n); 179. printf(NJ U V CETA n); /输出位移 180. for(i=1;i=nj;i+) 181. 182. printf( %-5d %14.11f %14.11f %14.11fn,i,p3*i-2,p3*i-1,p3*i); 183. 184. printf(_n); 185. /*根据E的值输出相应E单元的N,Q,M(A,B)的结果* 186. printf(E

10、N Q M n); 187. /*计算轴力N,剪力Q,弯矩M* 188. for(e=1;e=ne;e+) /按单元循环 189. 190. jdugd(e); /求局部单元刚度矩阵kd 191. zb(e); /求坐标变换矩阵T 192. for(i=1;i=2;i+) 193. 194. for(ii=1;ii=3;ii+) 195. 196. h=3*(i-1)+ii; 197. dh=3*(jmei-1)+ii; /给出整体坐标下单元节点位移 198. wyh=pdh; 199. 200. 201. for(i=1;i=6;i+) 202. 203. fi=0.0; 204. for(

11、j=1;j=6;j+) 205. 206. for(k=1;k0) 213. 214. for(i=1;i=npf;i+) /按非节点载荷数循环 215. if(pfi3=e) /找到荷载所在的单元 216. 217. hz=i; 218. gdnl(hz); /求固端反力 219. for(j=1;j=6;j+) /将固端反力累加 220. 221. fj=fj+f0j; 222. 223. 224. 225. printf(%-3d(A) %9.5f %9.5f %9.5fn,e,f1,f2,f3); /输出单元A（i）端力 226. printf( (B) %9.5f %9.5f %9.5fn,f4,f5,f6); /输出单元B（i）端力 227. 228. return; 229. 230. /*主程序结束* 231. 232. /* 233. / 234. /* 235. 236. void gdnl(int hz) 237. 238. int ind,e; 239. double g,c,l0,d; 240. 241. 242. g=pfhz1; /载荷值 243. c=pfhz2; /载荷位置 244. e=(int)pfhz3;