有限元编程的c++实现算例
1.#include
2.#include
3.
4.
5.#definene3 //单元数
6.#definenj4 //节点数
7.#definenz6 //支撑数
8.#definenpj0 //节点载荷数
9.#definenpf1 //非节点载荷数
10.#definenj312 //节点位移总数
11.#definedd6 //半带宽
12.#definee02.1E8 //弹性模量
13.#definea00.008 //截面积
14.#definei01.22E-4 //单元惯性距
15.#definepi
16.
17.
18.intjm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/
19.doublegc[ne+1]={0.0,1.0,2.0,1.0};
20.doublegj[ne+1]={0.0,90.0,0.0,90.0};
21.doublemj[ne+1]={0.0,a0,a0,a0};
22.doublegx[ne+1]={0.0,i0,i0,i0};
23.intzc[nz+1]={0,1,2,3,10,11,12};
24.doublepj[npj+1][3]={{0.0,0.0,0.0}};
25.doublepf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};
26.doublekz[nj3+1][dd+1],p[nj3+1];
27.doublepe[7],f[7],f0[7],t[7][7];
28.doubleke[7][7],kd[7][7];
29.
30.
31.//**kz[][]—整体刚度矩阵
32.//**ke[][]—整体坐标下的单元刚度矩阵
33.//**kd[][]—局部坐标下的单位刚度矩阵
34.//**t[][]—坐标变换矩阵
35.
36.//**这是函数声明
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,al,bl,m,l,dl,zl,z,j0;
47. doublecl,wy[7];
48. intim,in,jn;
49.
50.//***********************************************
51.//<功能：形成矩阵P>
52.//***********************************************
53.
54. if(npj>0)
55. {
56. for(i=1;i<=npj;i++)
57. { //把节点载荷送入P
58. j=pj[i][2];
59. p[j]=pj[i][1];
60. }
61. }
62. if(npf>0)
63. {
64. for(i=1;i<=npf;i++)
65. { //求固端反力F0
66. hz=i;
68. e=(int)pf[hz][3];
69. zb(e); //求单元号码
70. for(j=1;j<=6;j++) //求坐标变换矩阵T
71. {
72. pe[j]=0.0;
73. for(k=1;k<=6;k++) //求等效节点载荷
74. {
75. pe[j]=pe[j]-t[k][j]*f0[k];
76. }
77. }
78. al=jm[e][1];
79. bl=jm[e][2];
80. p[3*al-2]=p[3*al-2]+pe[1]; //将等效节点载荷送到P中
81. p[3*al-1]=p[3*al-1]+pe[2];
82. p[3*al]=p[3*al]+pe[3];
83. p[3*bl-2]=p[3*bl-2]+pe[4];
84. p[3*bl-1]=p[3*bl-1]+pe[5];
85. p[3*bl]=p[3*bl]+pe[6];
86. }
87. }
88.
89.
90. //*************************************************
91. //<功能：生成整体刚度矩阵kz[][]>
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*(jm[e][i]-1)+ii; //该元素在KZ中的行码
101. for(j=1;j<=2;j++)
103. for(jj=1;jj<=3;jj++) //对列码循环
104. {
105. l=3*(j-1)+jj; //元素在ke中的列码
106. zl=3*(jm[e][j]-1)+jj; //该元素在KZ中的行码 107. dl=zl-dh+1; //该元素在KZ*中的行码 108. if(dl>0)
109. kz[dh][dl]=kz[dh][dl]+ke[h][l]; //刚度集成
110. }
111. }
112. }
113. }
114. }
115.
116.//**引入边界条件**
117. for(i=1;i<=nz;i++) //按支撑循环
118. {
119. z=zc[i]; //支撑对应的位移数
120. kz[z][l]=1.0; //第一列置1
121. for(j=2;j<=dd;j++)
122. {
123. kz[z][j]=0.0; //行置0
124. }
125. if((z!=1))
126. {
127. if(z>dd)
128. j0=dd;
129. elseif(z<=dd)
130. j0=z; //列（45度斜线）置0
131. for(j=2;j<=j0;j++)
132. kz[z-j+1][j]=0.0;
133. }
134. p[z]=0.0; //P置0
135. }
136.
138.
139.
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=kz[k][l]/kz[k][1]; //修改KZ
151. jn=dd-l+1;
152. for(j=1;j<=jn;j++)
153. {
154. m=j+i-k;
155. kz[i][j]=kz[i][j]-cl*kz[k][m];
156. }
157. p[i]=p[i]-cl*p[k]; //修改P
158. }
159. }
160.
161.
162.
163.
164. p[nj3]=p[nj3]/kz[nj3][1]; //求最后一个位移分量 165. for(i=nj3-1;i>=1;i--)
166. {
167. if(dd>nj3-i+1)
168. j0=nj3-i+1;
169. elsej0=dd; //求最大列码j0
170. for(j=2;j<=j0;j++)
171. {
173. p[i]=p[i]-kz[i][j]*p[h];
174. }
175. p[i]=p[i]/kz[i][1]; //求其他位移分量
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.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);
183. }
184. printf("_____________________________________________________________\n"); 185. //*根据E的值输出相应E单元的N,Q,M(A,B)的结果**
186. printf("E 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*(jm[e][i]-1)+ii; //给出整体坐标下单元节点位移
198. wy[h]=p[dh];
199. }
200. }
201. for(i=1;i<=6;i++)
202. {
203. f[i]=0.0;
204. for(j=1;j<=6;j++)
205. {
206. for(k=1;k<=6;k++) //求由节点位移引起的单元节点力
208. f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];
209. }
210. }
211. }
212. if(npf>0)
213. {
214. for(i=1;i<=npf;i++) //按非节点载荷数循环
215. if(pf[i][3]==e) //找到荷载所在的单元
216. {
217. hz=i;
218. gdnl(hz); //求固端反力
219. for(j=1;j<=6;j++) //将固端反力累加
220. {
221. f[j]=f[j]+f0[j];
222. }
223. }
224. }
225. printf("%-3d(A) %9.5f %9.5f %9.5f\n",e,f[1],f[2],f[3]); //输出单元A（i）端内力 226. printf(" (B) %9.5f %9.5f %9.5f\n",f[4],f[5],f[6]); //输出单元B（i）端内力 227. }
228. return;
229.}
230.//**主程序结束**
231.
232.//************************************************************
233.//<功能：将非节点载荷下的杆端力计算出来存入f0[]>
234.//************************************************************
235.
236.voidgdnl(inthz)
237.{
238. intind,e;
239. doubleg,c,l0,d;
240.
241.
243. c=pf[hz][2]; //载荷位置
244. e=(int)pf[hz][3]; //作用单元
245. ind=(int)pf[hz][4]; //载荷类型
246. l0=gc[e]; //杆长
247. d=l0-c;
248. if(ind==1)
249. {
250. f0[1]=0.0;
251. f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2; //均布载荷的固端反力 252. f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12;
253. f0[4]=0.0;
254. f0[5]=-g*c-f0[2];
255. f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);
256. }
257. else
258. {
259. if(ind==2) //横向集中力的固端反力
260. {
261. f0[1]=0.0;
262. f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);
263. f0[3]=-(g*c*d*d)/(l0*l0);
264. f0[4]=0.0;
265. f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);
266. f0[6]=(g*c*c*d)/(l0*l0);
267. }
268. else
269. {
270. f0[1]=-(g*d/l0); //纵向集中力的固端反力
271. f0[2]=0.0;
272. f0[3]=0.0;
273. f0[4]=-g*c/l0;
274. f0[5]=0.0;
275. f0[6]=0.0;
276. }
278.}
279.
280.//************************************************************ 281.//<功能：构成坐标变换矩阵>
282.//************************************************************ 283.voidzb(inte)
284.{
285. doubleceta,co,si;
286. inti,j;
287. ceta=(gj[e]*pi)/180; //角度变弧度
288. co=cos(ceta);
289. si=sin(ceta);
290. t[1][1]=co; //计算T右上角元素
291. t[1][2]=si;
292. t[2][1]=-si;
293. t[2][2]=co;
294. t[3][3]=1.0;
295. for(i=1;i<=3;i++)
296. {
297. for(j=1;j<=3;j++) //计算T的左下角元素
298. {
299. t[i+3][j+3]=t[i][j];
300. }
301. }
302.}
303.
304.
305.
306.//*****************************************************
307.//<计算局部坐标下单元刚度矩阵kd[][]>
308.//*****************************************************
309.voidjdugd(inte)
310.{
311. doubleA0,l0,j0;
313. intj;
314.
315.
316. A0=mj[e]; //面积
317. l0=gc[e]; //杆长
318. j0=gx[e]; //惯性钜
319.
320.
321. for(i=0;i<=6;i++)
322. for(j=0;j<=6;j++) //kd清0
323. kd[i][j]=0.0;
324.
325. kd[1][1]=e0*A0/l0;
326. kd[2][2]=12*e0*j0/pow(l0,3);
327. kd[3][2]=6*e0*j0/pow(l0,2);
328. kd[3][3]=4*e0*j0/l0;
329. kd[4][1]=-kd[1][1];
330. kd[4][4]=kd[1][1];
331. kd[5][2]=-kd[2][2]; //计算kd左下角各元素
332. kd[5][3]=-kd[3][2];
333. kd[5][5]=kd[2][2];
334. kd[6][2]=kd[3][2];
335. kd[6][3]=2*e0*j0/l0;
336. kd[6][5]=-kd[3][2];
337. kd[6][6]=kd[3][3];
338.
339. for(i=1;i<=6;i++)
340. for(j=1;j<=i;j++) //将kd左下角元素按对称原则送到右下角 341. kd[j][i]=kd[i][j];
342.}
343.
344.
345.//********************************************************** 346.//<计算整体坐标下单元刚度矩阵ke[][]>
348.voiddugd(inte)
349.{
350. inti,k,j,m;
351. jdugd(e); //计算局部单元刚度矩阵kd
352. zb(e); //计算坐标变换矩阵T
353. for(i=1;i<=6;i++)
354. {
355. for(j=1;j<=6;j++)
356. {
357. ke[i][j]=0.0;
358. for(k=1;k<=6;k++)
359. {
360. for(m=1;m<=6;m++)
361. {
362. ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; //计算刚度坐标内单元刚度矩阵ke 363. }
364. }
365. }
366. }
367.}
368.
369.
370.//**程序结束**
371.
372.
373.
374.