有限元编程的c++实现算例.docx
《有限元编程的c++实现算例.docx》由会员分享，可在线阅读，更多相关《有限元编程的c++实现算例.docx（11页珍藏版）》请在冰点文库上搜索。
有限元编程的c++实现算例
有限元编程的c++实现算例
read.pudn./downloads76/doc/fileformat/290377/ganjian.cpp__.htm
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.#definepi3.141592654
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;
67. gdnl（hz）;
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++）
102. {
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.
137.
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. {
172. h=j+i-1;
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++） //求由节点位移引起的单元节点力
207. {
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.
242. g=pf[hz][1]; //载荷值
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.
升级会员 