有限元上机.docx
《有限元上机.docx》由会员分享，可在线阅读，更多相关《有限元上机.docx（30页珍藏版）》请在冰豆网上搜索。
有限元上机
弹性力学及有限元基础
上机实践报告
指导老师：
单丽君
班级：
机械093班
学号：
0904010321
姓名：
尉骚狗
用常应变三角形单元解弹性力学平面问题
本程序只适用于三结点三角形单元，可计算平面应力问题也可计算平面应变问题，这两类问题用IND来区别，IND输入0的数据为解决平面应力问题，IND输入1的数据为解决平面应变问题。
本程序运行环境为VisualC++。
一、程序说明
1.程序中记号说明
nj—结点个数；ne—单元个数；nz—支杆个数；
ndd—半带宽；ind—问题类型码；nj2—位移分量个数；
eo—弹性模量；un—泊松比；gama—材料容重；
te—单元厚度；ae—单元面积；JM—单元结点码数组；
NZC—支撑数组；CJZ—结点坐标数组；PJ—结点载荷数组；
b[]—几何矩阵；d[]—弹性矩阵；s[]—应力矩阵；
tkz—半带存储的整体刚度矩阵；eke—单元刚度矩阵；p—载荷向量；
npj1—结点载荷个数加1；meo—任一单源码；iask—子程序的形参；
iask=1求单元面积AE；iask=2求应力矩阵S；iask=3求单元刚度K；
ie、je、me—单元meo的三个角点码；cm、bm、cj、bj—计算形函数中的常数；
i、j、k、ii、jj—循环参数；lh—单元行码；ldh—半带存储的整体刚度矩阵中的行码；
l—单元列码；ld—半带存储的整体刚度矩阵中的列码；
i1—工作单元；pe—自重的等效结点载荷；mz—支杆相应的位移分量码；
jo—最大列码；j1—工作单元；im—最大行码；
c—系数比值；ld1、m—工作单元；wy—单元结点位移量；
yl—应力向量；sigx、sigy、toxy—应力的三个分量；pyl—平均应力；
ryl—应力圆半径；sig—工作单元；sig1、sig2—最大最小主应力；
ceta1—工作单元；ceta—主平面角。
2.子程序名
DATA—输入数据子程序；
ELEST（meo，iask）—求单元meo的面积A；
TOTSTI—形成整体刚度矩阵子程序；
LOAD—形成载荷向量子程序；
SUPPOR—引入支撑子程序；
SOLVEQ—解方程输出位移子程序；
STRESS—求应力输出应力子程序。
二、程序流程图
1.程序总框图
图1程序总框图
整个程序由一个主调主程序（主函数main（））和七个子程序组成，其中数据子程序DATA（）用来接受输入的参数和变量。
单刚子程序ELEST（）为含有参数的函数，功能控制参数iask可取1、2、3，分别计算相应单元的面积、应力矩阵和单元刚度矩阵，主程序不直接调用它，而通过总刚子程序、载荷子程序和求应力子程序间接调用。
总刚子程序TOTSTI（）用来合成总刚矩阵，载荷子程序LOAD（）用来计算合成载荷，支承子程序SUPPOR（）用来引入约束，解方程子程序SOLVEQ（）用来求解并输出各个结点的位移，求应力子程序STRESS（）用来计算和输出应力、主应力及主平面角。
2.主函数及各子程序流程图
图2voidmain（）和voidDATA（）
图3voidELEST（intmeo,intiask）
图4voidTOTSTI（）图5voidLOAD（）
图6voidSUPPOR（）
图7voidSOLVEQ（）
图8voidSTRESS（）
三、程序源代码
#include"stdio.h"
#include"math.h"
intnj,ne,nz,ndd,ind,nj2;
intjm[100][3],nzc[200],npj1,npj;
floateo,un,gama,te,ae;
floatcjz[100][2],pj[100][2];
floatb[3][6],d[3][3],s[3][6],eke[6][6],tkz[200][20],p[200];
voidDATA（）/*数据输入函数*/
{
inti,j;
printf（"§请输入6个基本参数，结点个数、单元个数、支杆个数、半带宽、结点载荷个数、问题类型码:
\n"）;
scanf（"%d,%d,%d,%d,%d,%d",&nj,&ne,&nz,&ndd,&npj,&ind）;
nj2=nj*2;
npj1=npj+1;
printf（"\n"）;
printf（"§请输入4个基本常量，弹性模量、泊松比、材料容重、单位厚度:
\n"）;
scanf（"%f,%f,%f,%f",&eo,&un,&gama,&te）;
printf（"\n"）;
printf（"nj=%d,ne=%d,nz=%d,ndd=%d,npj=%d,ind=%d\n",nj,ne,nz,ndd,npj,ind）;
printf（"eo=%f,un=%f,gama=%f,te=%f\n",eo,un,gama,te）;
printf（"\n"）;
printf（"§请输入JM矩阵:
\n"）;
for（i=0;i{
for（j=0;j<3;j++）
scanf（"%d",&jm[i][j]）;
getchar（）;
}
printf（"\n"）;
printf（"§请输入CJZ矩阵:
\n"）;
for（i=0;i{
for（j=0;j<2;j++）
scanf（"%f",&cjz[i][j]）;
getchar（）;
}
printf（"\n"）;
printf（"§请输入NZC矩阵:
\n"）;
for（i=0;iscanf（"%d",&nzc[i]）;
getchar（）;
printf（"\n"）;
printf（"§请输入PJ矩阵:
\n"）;
for（i=0;i{
for（j=0;j<2;j++）
scanf（"%f",&pj[i][j]）;
getchar（）;
}
printf（"\n"）;
printf（"CJZ矩阵如下:
\n"）;
for（i=0;i{
for（j=0;j<2;j++）
printf（"%f\t",cjz[i][j]）;
printf（"\n"）;
}
}
voidELEST（intmeo,intiask）/*单元刚度矩阵*/
{
intie,je,me,i,j,k;
floatcm,bm,cj,bj;
ie=jm[meo-1][0];
je=jm[meo-1][1];
me=jm[meo-1][2];
cm=cjz[je-1][0]-cjz[ie-1][0];
bm=cjz[ie-1][1]-cjz[je-1][1];
cj=cjz[ie-1][0]-cjz[me-1][0];
bj=cjz[me-1][1]-cjz[ie-1][1];
ae=（bj*cm-bm*cj）/2;
if（iask>1）
{
for（i=0;i<3;i++）
{
for（j=0;j<6;j++）
b[i][j]=0;
}
b[0][0]=-bj-bm;
b[0][2]=bj;
b[0][4]=bm;
b[1][1]=-cj-cm;
b[1][3]=cj;
b[1][5]=cm;
b[2][0]=-cj-cm;
b[2][1]=-bj-bm;
b[2][2]=cj;
b[2][3]=bj;
b[2][4]=cm;
b[2][5]=bm;
for（i=0;i<3;i++）
{
for（j=0;j<6;j++）
b[i][j]=b[i][j]/（2*ae）;
}
d[0][0]=eo/（1-un*un）;
d[0][1]=eo*un/（1-un*un）;
d[0][2]=0;
d[1][0]=eo*un/（1-un*un）;
d[1][1]=eo/（1-un*un）;
d[1][2]=0;
d[2][0]=0;
d[2][1]=0;
d[2][2]=eo/（2*（1+un））;
for（i=0;i<3;i++）
{
for（j=0;j<6;j++）
{
s[i][j]=0;
for（k=0;k<3;k++）
s[i][j]=s[i][j]+d[i][k]*b[k][j];
}
}
if（iask>2）
{
for（i=0;i<6;i++）
{
for（j=0;j<6;j++）
{
eke[i][j]=0;
for（k=0;k<3;k++）
eke[i][j]=eke[i][j]+s[k][i]*b[k][j]*ae*te;
}
}
}
}
}
voidTOTSTI（）/*总体刚度矩阵*/
{
inti,j,ii,jj,meo,lii,l,lz,ldh,ld;
for（i=0;i{
for（j=0;jtkz[i][j]=0;
}
for（meo=1;meo<=ne;meo++）
{
ELEST（meo,3）;
for（i=1;i<4;i++）
{
for（ii=1;ii<3;ii++）
{
lii=2*（i-1）+ii;
ldh=2*（jm[meo-1][i-1]-1）+ii;
for（j=1;j<4;j++）
{
for（jj=1;jj<3;jj++）
{
l=2*（j-1）+jj;
lz=2*（jm[meo-1][j-1]-1）+jj;
ld=lz-ldh+1;
if（ld>0）
tkz[ldh-1][ld-1]=tkz[ldh-1][ld-1]+eke[lii-1][l-1];
}
}
}
}
}
}
voidLOAD（）/*添加载荷*/
{
inti,j,meo,ie,je,me;
floatp0;
for（i=0;ip[i]=0;
if（npj>0）
{
for（i=1;i{
j=（int）（pj[i][1]）;
p[j-1]=pj[i][0];
}
}
if（gama>0）
{
for（meo=1;meo<=ne;meo++）
{
ELEST（meo,1）;
p0=-gama*ae*te/3;
ie=jm[meo-1][0];
je=jm[meo-1][1];
me=jm[meo-1][2];
p[2*ie-1]=p[2*ie-1]+p0;
p[2*je-1]=p[2*je-1]+p0;
p[2*me-1]=p[2*me-1]+p0;
}
}
}
voidSUPPOR（）/*添加约束*/
{
inti,j,mz,j0;
for（i=1;i<=nz;i++）
{
mz=nzc[i-1];
tkz[mz-1][0]=1;
for（j=1;jtkz[mz-1][j]=0;
if（mzj0=mz;
else
j0=ndd;
for（j=2;j<=j0;j++）
tkz[mz-j][j-1]=0;
p[mz-1]=0;
}
}
voidSOLVEQ（）/*求解过程*/
{
intk,i,j,im,l,m,ii,jo,lh;
floatc;
for（k=1;k<=nj2-1;k++）
{
if（nj2>=k+ndd-1）
im=nj2;
else
im=k+ndd-1;
for（i=k+1;i<=im;i++）
{
l=i-k+1;
c=tkz[k-1][l-1]/tkz[k-1][0];
for（j=1;j<=ndd-l+1;j++）
{
m=j+i-k;
tkz[i-1][j-1]=tkz[i-1][j-1]-c*tkz[k-1][m-1];
}
p[i-1]=p[i-1]-c*p[k-1];
}
}
p[nj2-1]=p[nj2-1]/tkz[nj2-1][0];
for（ii=1;ii<=nj2-1;ii++）
{
i=nj2-ii;
if（ndd<=nj2-i+1）
jo=ndd;
else
jo=nj2-i+1;
for（j=2;j<=jo;j++）
{
lh=j+i-1;
p[i-1]=p[i-1]-tkz[i-1][j-1]*p[lh-1];
}
p[i-1]=p[i-1]/tkz[i-1][0];
}
printf（"运算结果：
\n"）;
printf（"uv\n"）;
for（i=0;i{printf（"%16f",p[i]）;
if（i%2==1）
printf（"\n"）;
}
printf（"\n"）;
}
voidSTRESS（）/*应力计算*/
{
doubleceta,cetal,sigx,sigy,toxy,pyl,ryl,sig1,sig2;
intmeo,i,j,lh,ldh;
floatwy[6],yl[3];
for（meo=1;meo<=ne;meo++）
{
ELEST（meo,2）;
for（i=1;i<=3;i++）
{
for（j=1;j<=2;j++）
{
lh=2*（i-1）+j;
ldh=2*（jm[meo-1][i-1]-1）+j;
wy[lh-1]=p[ldh-1];
}
}
for（i=1;i<=3;i++）
{
yl[i-1]=0;
for（j=1;j<=6;j++）
yl[i-1]=yl[i-1]+s[i-1][j-1]*wy[j-1];
}
sigx=yl[0];
sigy=yl[1];
toxy=yl[2];
pyl=（sigx+sigy）/2;
ryl=sqrt（（sigx-sigy）*（sigx-sigy）/4+toxy*toxy）;
sig1=pyl+ryl;
sig2=pyl-ryl;
if（sigy==sig2）
ceta=0;
else
{
cetal=toxy/（sigy-sig2）;
ceta=90-57.29578*atan（cetal）;
}
printf（"meo=%d\nsigx=%fsigy=%ftoxy=%f\nsig1=%fsig2=%fceta=%f\n",
meo,sigx,sigy,toxy,sig1,sig2,ceta）;
}
}
voidmain（）/*主程序*/
{
printf（"★★★★★★★●JIEGUOSHUCHU●★★★★★★★\n"）;
printf（"\n"）;
DATA（）;
if（ind!
=0）
{
eo=eo/（1-un*un）;
un=un/（1-un）;
}
TOTSTI（）;
LOAD（）;
SUPPOR（）;
SOLVEQ（）;
STRESS（）;
printf（"\n"）;
printf（"\n"）;
}
四、工程算例
如图所示钢梁，受
、
压力作用，以平面问题计算取
将结构分成18个单元，按图中所示的网格划分，用程序计算结点位移及各单元中的应力，其中弹性模量
。
输入数据：
基本参数
结点数nj
单元数ne
支杆数nz
半带宽ndd
结点载荷数npj
问题类型码ind
36
34
9
6
3
0
其它数据
弹性模量eo
泊松比un
材料密度gama
单元厚度te
1.0
0.167
0.0
1.0
单元结点码数组JM矩阵,见程序运行结果。
结点坐标数组CJZ矩阵，见程序运行结果。
支承数组NZC矩阵，见程序运行结果。
结点载荷数组PJ矩阵，见程序运行结果。
运行过程：
★★★★★★★●JIEGUOSHUCHU●★★★★★★★
§请输入6个基本参数，结点个数、单元个数、支杆个数、半带宽、结点载荷个数、问题类
型码:
36,34,9,12,3,0
§请输入4个基本常量，弹性模量、泊松比、材料容重、单位厚度:
1.0,0.167,0.0,1.0
nj=36,ne=34,nz=9,ndd=12,npj=3,ind=0
eo=1.000000,un=0.167000,gama=0.000000,te=1.000000
§请输入JM矩阵:
232724
192324
192420
151920
152016
111516
111612
71112
7128
478
485
245
253
123
356
596
6910
91310
101314
131714
141718
172118
182122
212622
212526
252926
252829
283129
283031
303331
303233
323533
323435
343635
§请输入CJZ矩阵:
0.00.0
4.00.0
0.03.0
8.00.0
4.03.0
0.06.0
12.00.0
8.03.0
4.06.0
0.09.0
16.00.0
12.03.0
4.09.0
0.012.0
20.00.0
16.03.0
4.012.0
0.015.0
24.00.0
20.03.0
4.015.0
0.018.0
28.00.0
24.03.0
8.015.0
4.018.0
28.03.0
12.015.0
8.018.0
16.015.0
12.018.0
20.015.0
16.018.0
24.015.0
20.018.0
28.018.0
§请输入NZC矩阵:
1211144346566768
§请输入PJ矩阵:
0.00.0
-250.025.0
-240.040.0
-245.058.0
CJZ矩阵如下:
0.0000000.000000
4.0000000.000000
0.0000003.000000
8.0000000.000000
4.0000003.000000
0.0000006.000000
12.0000000.000000
8.0000003.000000
4.0000006.000000
0.0000009.000000
16.0000000.000000
12.0000003.000000
4.0000009.000000
0.00000012.000000
20.0000000.000000
16.0000003.000000
4.00000012.000000
0.00000015.000000
24.0000000.000000
20.0000003.000000
4.00000015.000000
0.00000018.000000
28.0000000.000000
24.0000003.000000
8.00000015.000000
4.00000018.000000
28.0000003.000000
12.00000015.000000
8.00000018.000000
16.00000015.000000
12.00000018.000000
20.00000015.000000
16.00000018.000000
24.00000015.000000
20.00000018.000000
28.00000018.000000
运算结果：
uv
0.0000000.000000
464.764801135.639404
238.640747-89.400124
797.274963285.189545
391.376770152.356354
0.000000-155.526627
1015.2744750.000000
951.806152251.767960
-78.154236200.602188
-456.948700-70.193687
1439.958374-1622.574097
1672.331665-188.007919
-605.509399135.125992
-425.82330382.762360
2264.695313-2831.118652
2212.583740-1774.606567
-409.131073-9.852731
-184.436066121.669426
3292.391602-2175.435059
2351.984619-3015.543945
-127.921700-66.468155
0.00000098.286995
4165.24853547.616310
2266.479736-2228.899414
56.311920-75.690674
-43.825157-83.433342
2302.066162138.425858
381.356018-70.290672
-43.175243-86.627029
43.7331622414.275635
-206.516098287.141296
85.072670431.341309
-95.6151202164.874023
14.9687812704.031982
871.73400997.624718
687.739441-503.969238
meo=1
sigx=14.351885sigy=32.666618toxy=-12.523359
sig1=39.023500sig2=7.995003ceta=116.912438
meo=2
sigx=221.413055sigy=19.154507toxy=91.599068
sig1=256.729812s
升级会员 