有限元上机.docx
- 文档编号:3765703
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:30
- 大小:279.09KB
有限元上机.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;i scanf("%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;j tkz[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;i p[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;j tkz[mz-1][j]=0; if(mz j0=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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 上机