有限元编程的c++实现算例.docx
- 文档编号:11398931
- 上传时间:2023-02-28
- 格式:DOCX
- 页数:20
- 大小:19.59KB
有限元编程的c++实现算例.docx
《有限元编程的c++实现算例.docx》由会员分享,可在线阅读,更多相关《有限元编程的c++实现算例.docx(20页珍藏版)》请在冰豆网上搜索。
有限元编程的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];
30.
31.//**kz[][]一整体刚度矩阵
32.//**ke[][]一整体坐标下的单元刚度矩阵
33.//**kd[][]一局部坐标下的单位刚度矩阵
34.//**皿]一坐标变换矩阵
35.
35.//**这是函数声明
36.voidjdugd(int);
37.voidzb(int);
38.voidgdnl(int);
39.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.
53.if(npj>0)
54.{
55.for(i=1;i<=npj;i++)
56.
〃把节点载荷送入P
{
57.j=pj[i][2];
58.p[j]=pj[i][1];
59.}
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];
//将等效节点载荷送到
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++)
//按单元循环
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.
165.
166.
167.
168.
169.
170.
171.
172.
173.
174.
175.
176.
177.
178.
179.
180.
181.
182.
183.
184.
185.
186.
187.
188.
189.
p[nj3]=p[nj3]/kz[nj3][1];
for(i=nj3-1;i>=1;i--)
(
if(dd>nj3-i+1)
j0=nj3-i+1;
elsej0=dd;
for(j=2;j<=j0;j++)
(
h=j+i-1;
p[i]=p[i]-kz[i][j]*p[h];
)
p[i]=p[i]/kz[i][1];
)
printf("\n");
printf("
//求最后一个位移分量
//求最大列码j0
//求其他位移分量
\n");
printf("NJUVCETA\n");//输出位移
for(i=1;i<=nj;i++)
(
printf("%-5d%14.11f%14.11f%14.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);
)
printf("\n");
//*根据E的值输出相应E单元的N,Q,M(A,B)的结果**
printf("ENQM\n");
//*计算轴力N,剪力Q,弯矩M*
for(e=1;e<=ne;e++)//按单元循环
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]);
⑴端力
226.
printf("(B)%9.5f
%9.5f
%9.5f\n",f[4],f[5],f[6]);
端力
〃输出单元A
//输出单元B⑴
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.
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.}
277.}
278.}
280.//*
279.
火味*味味*味味*味味*味味味*味味*味味*味味*味味味*味味*味味*味味*味味味*味味*味味*味味*味味味*味味**
281.//〈功能:
构成坐标变换矩阵>
282.//*
***********************************************************
283.voidzb(inte)
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;
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++)
298.
(
299.
t[i+3][j+3]=t[i][j];
300.
}
301.
}
302.}
284.(
//角度变弧度
//计算T右上角元素
//计算T的左下角元素
303.
304.
305.
306.//*
火味*味味*味味*味味*味味味*味味*味味*味味*味味味*味味*味味*味味*味味味*味味*味味*味味**
307.//<计算局部坐标下单元刚度矩阵kd[][]>
308.//*
****************************************************
309.voidjdugd(inte)
310.(
311.doubleA0,l0,j0;
312.inti;
313.intj;
314.
315.
316.A0=mj[e];
317.l0=gc[e];
318.j0=gx[e];
319.
320.
319.for(i=0;i<=6;i++)
320.for(j=0;j<=6;j++)
321.kd[i][j]=0.0;
324.
322.kd[1][1]=e0*A0/l0;
323.kd[2][2]=12*e0*j0/pow(l0,3);
324.kd[3][2]=6*e0*j0/pow(l0,2);
325.kd[3][3]=4*e0*j0/l0;
326.kd[4][1]=-kd[1][1];
327.kd[4][4]=kd[1][1];
328.kd[5][2]=-kd[2][2];
329.kd[5][3]=-kd[3][2];
330.kd[5][5]=kd[2][2];
331.kd[6][2]=kd[3][2];
332.kd[6][3]=2*e0*j0/l0;
333.kd[6][5]=-kd[3][2];
334.kd[6][6]=kd[3][3];
338.
335.for(i=1;i<=6;i++)
336.for(j=1;j<=i;j++)
337.kd[j][i]=kd[i][j];
338.}
343.
344.
//面积
//杆长
〃惯性镇
//kd清0
//计算kd左下角各元素
//将kd左下角元素按对称原则送到右下角
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
ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j];//计算刚度坐标单元刚度矩阵
363.}
364.}
365.}
366.}
367.}
368.
369.
368.//**程序结束**
371.
372.
373.
374.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 编程 c+ 实现