弹性力学三结点三角形单元静力分析程序.docx
- 文档编号:25844961
- 上传时间:2023-06-16
- 格式:DOCX
- 页数:53
- 大小:27.27KB
弹性力学三结点三角形单元静力分析程序.docx
《弹性力学三结点三角形单元静力分析程序.docx》由会员分享,可在线阅读,更多相关《弹性力学三结点三角形单元静力分析程序.docx(53页珍藏版)》请在冰豆网上搜索。
弹性力学三结点三角形单元静力分析程序
弹性力学三结点三角形单元静力分析程序
1)先处理法程序PreFin
/**************************************************************************
***************
********平面问题三结点三角形单元有限元分析程序(先处理法)*******
***************
**************************************************************************/
/*************************************************************************
***************
********三结点三角形单元平面问题有限元分析程序(先处理法)*******
***************
**************************************************************************/
#include"Func.h"
#include
#include"math.h"
voidGetEleMainStress(intElementNum,double*Stress,double*MainStress);
voidGetElementStress(intEleNum,doubleMiu,boolIsPlainStrain,doubleElasticModule,int*EleNode,int*NodeDOF,double*Xcoord,double*Ycoord,double*Disp,double*Stress);
voidCalcuDOFIndex(intRestraintNodeNum,int*ResNodeDOF,intNodeNum,int&FreeDOF,int*NodeDOF);
voidGetStiffness(intEleNum,int*EleNode,int*ResNodeDOF,double*Xcoord,double*Ycoord,doubleElasticModule,doubleMiu,doubleThick,boolIsPlainStrain,doubleK[][6]);
voidgauss(double*a,double*b,intnTotal,intnFree);
voidStiffAssemble(intEleNum,intTotalDOF,int*EleNode,double*GlobalK,int*NodeDOF,doubleK[][6]);
voidInputControlPar(ifstream&fin,int&ElementNum,int&NodeNum,int&RestraintNodeNum,int&LoadNodeNum,int&SurLoadNum);
voidInputOtherPar(ifstream&fin,int&ElementNum,int&NodeNum,int&SurLoadNum,int&RestraintNodeNum,int&LoadNodeNum,double&ElasticModule,double&Miu,double&Density,double&Thick,double*Xcoord,double*Ycoord,int*ResNodeDOF,int*EleNode,double*NodeLoad,double*SurLoad);
voidOutPutPar(ofstream&fout,int&ElementNum,int&NodeNum,int&SurLoadNum,int&RestraintNodeNum,int&LoadNodeNum,double&ElasticModule,double&Miu,double&Density,double&Thick,double*Xcoord,double*Ycoord,int*ResNodeDOF,int*EleNode,double*NodeLoad,double*SurLoad,int*NodeDOF);
doubleGetArea(intEleNum,int*EleNode,double*Xcoord,double*Ycoord);
voidGetLoadVector(intLoadNodeNum,intSurLoadNum,doubleDensity,intElementNum,doubleThick,int*EleNode,int*NodeDOF,double*Xcoord,double*Ycoord,double*NodeLoad,double*SurLoad,double*LoadorDisp);
voidGetForce(intNodeNum,double*GlobalK,double*LoadorDisp,double*Force);
///////////主程序,输入参数,调用子函数进行计算
voidmain()
{
inti;
//申明变量
intElementNum,NodeNum,RestraintNodeNum,LoadNodeNum;//单元数、结点数、约束结点数、荷载数
doubleElasticModule,Miu;//材料弹性模量,泊松比
boolIsPlainStrain=false;//缺省是平面应力问题
doubleThick;//厚度
doubleDensity;//密度,不考虑可以自重,可以使其等于0;
doubleK[6][6];//单元刚度矩阵
intTotalDOF,FreeDOF;//总自由度、有效自由度数
ifstreamfin("例题5.1.txt");//定义输入文件;
ofstreamfout("result.txt");//定义输出文件;
////
intSurLoadNum;//分布荷载个数
///////////////////////////////////////////////////输入控制参数
InputControlPar(fin,ElementNum,NodeNum,RestraintNodeNum,LoadNodeNum,SurLoadNum);
TotalDOF=2*NodeNum;
/////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////计算过程
int*NodeDOF=newint[NodeNum*2];//单元自由度编码
int*ResNodeDOF=newint[RestraintNodeNum*3];//单元约束自由度
int*EleNode=newint[ElementNum*3];//单元结点编码
double*Xcoord=newdouble[NodeNum];//结点X坐标
double*Ycoord=newdouble[NodeNum];//结点Y坐标
double*GlobalK=newdouble[TotalDOF*TotalDOF];//单元结点编码
double*NodeLoad=newdouble[LoadNodeNum*3];//结点荷载
double*SurLoad=newdouble[SurLoadNum*5];//线荷载
double*LoadorDisp=newdouble[TotalDOF];//线荷载
double*Stress=newdouble[3];//应力
double*MainStress=newdouble[3];//主应力
double*Force=newdouble[NodeNum*2];//力
InputOtherPar(fin,ElementNum,NodeNum,SurLoadNum,RestraintNodeNum,LoadNodeNum,ElasticModule,Miu,Density,Thick,Xcoord,Ycoord,ResNodeDOF,EleNode,NodeLoad,SurLoad);
///////////////////////////重新进行结点编码
CalcuDOFIndex(RestraintNodeNum,ResNodeDOF,NodeNum,FreeDOF,NodeDOF);
OutPutPar(fout,ElementNum,NodeNum,SurLoadNum,RestraintNodeNum,LoadNodeNum,ElasticModule,Miu,Density,Thick,Xcoord,Ycoord,ResNodeDOF,EleNode,NodeLoad,SurLoad,NodeDOF);
//总刚度矩阵及荷载向量赋初值
for(i=0;i { for(intj=0;j GlobalK[i*TotalDOF+j]=0.0; LoadorDisp[i]=0.0; } for(i=0;i { GetStiffness(i,EleNode,ResNodeDOF,Xcoord,Ycoord,ElasticModule,Miu,Thick,IsPlainStrain,K); StiffAssemble(i,TotalDOF,EleNode,GlobalK,NodeDOF,K); //输出单元刚度矩阵 fout<<"输出单元刚度矩阵"< for(intii=0;ii<6;ii++) { for(intj=0;j<6;j++) fout< fout< } } //输出总刚度矩阵 fout<<"输出总刚度矩阵"< for(i=0;i { for(intj=0;j fout< fout< } GetLoadVector(LoadNodeNum,SurLoadNum,Density,ElementNum,Thick,EleNode,NodeDOF,Xcoord,Ycoord,NodeLoad,SurLoad,LoadorDisp); fout<<"结点等效荷载列阵"< for(i=0;i { fout< } fout< for(i=FreeDOF;i gauss(GlobalK,LoadorDisp,TotalDOF,FreeDOF); fout<<"结点位移为: "; for(i=0;i { fout< } fout< ///////////////////////////////////////////////////////////////////////////////////// for(i=0;i { GetElementStress(i,Miu,IsPlainStrain,ElasticModule,EleNode,NodeDOF,Xcoord,Ycoord,LoadorDisp,Stress); GetEleMainStress(i,Stress,MainStress); fout<<"单元"< "< "< "< fout<<"单元"< "< "< "< } /////////////////////////////输出力 GetForce(NodeNum,GlobalK,LoadorDisp,Force); fout<<"结点力"; for(intj=0;j<2*NodeNum;j++) fout< fout< /////////////////////////////////////////////////////////////////////////////////////输出参数 deleteNodeDOF;//单元自由度编码 deleteResNodeDOF;//单元约束自由度 deleteEleNode;//单元结点编码 deleteXcoord;//结点X坐标 deleteYcoord;//结点Y坐标 deleteGlobalK;//单元结点编码 deleteNodeLoad;//结点荷载 deleteSurLoad;//线荷载 deleteLoadorDisp;//总荷载或位移 deleteStress;//应力 deleteMainStress;//主应力 deleteForce;//力 fout.close(); fin.close(); } voidCalcuDOFIndex(intRestraintNodeNum,int*ResNodeDOF,intNodeNum,int&FreeDOF,int*NodeDOF) {//计算自由度重新编号 inti,j; intDOFindex=0;//自由度编码 /////////////////// for(i=0;i for(j=1;j<3;j++) { if(ResNodeDOF[i*3+j]>0)DOFindex++; } intn_TotalDOF=NodeNum*2; FreeDOF=n_TotalDOF-DOFindex; /////////////////求得无约束自由度数 DOFindex=0;//自由度编码 for(intiNode=0;iNode { boolisConstraint=false; for(intjResNode=0;jResNode { if(iNode==ResNodeDOF[jResNode*3+0]) { if(ResNodeDOF[jResNode*3+1]<=0) { NodeDOF[iNode*2+0]=DOFindex; DOFindex++; } if(ResNodeDOF[jResNode*3+2]<=0) { NodeDOF[iNode*2+1]=DOFindex; DOFindex++; } isConstraint=true; } } if(! isConstraint) { NodeDOF[iNode*2+0]=DOFindex;DOFindex++; NodeDOF[iNode*2+1]=DOFindex;DOFindex++; } } ////////////////////////////处理约束编号后 for(intjResNode=0;jResNode { if(ResNodeDOF[jResNode*3+1]>0) { NodeDOF[ResNodeDOF[jResNode*3+0]*2+0]=DOFindex; DOFindex++; } if(ResNodeDOF[jResNode*3+2]>0) { NodeDOF[ResNodeDOF[jResNode*3+0]*2+1]=DOFindex; DOFindex++; } } } ///////////////////////////////////////////// ////////////////////////////////集成总刚度矩阵********* voidStiffAssemble(intEleNum,intTotalDOF,int*EleNode,double*GlobalK,int*NodeDOF,doubleK[][6]) { intloop,loop1,iBuf,DOFIndex[6]; //求得单元三个结点的自由度 for(intj=0;j<3;j++) { intnNode=EleNode[EleNum*3+j];//单元三个结点的编号 for(inti=0;i<2;i++) { DOFIndex[i+2*j]=NodeDOF[nNode*2+i]; } } for(loop=0;loop<6;loop++) { iBuf=DOFIndex[loop]; for(loop1=0;loop1<6;loop1++) { intiBuf1=DOFIndex[loop1]; GlobalK[iBuf*TotalDOF+iBuf1]+=K[loop][loop1]; } } } /////////////// /////高斯消去法求解线性方程组 voidgauss(double*a,double*b,intnTotal,intnFree) { inti,j,k; doublea1; double*EffectA=newdouble[nFree*nFree];//单元结点编码 for(i=0;i { for(intj=0;j { EffectA[i*nFree+j]=a[i*nTotal+j]; } } for(k=0;k { for(i=k+1;i { a1=EffectA[k*nFree+i]/EffectA[k*nFree+k]; for(j=k+1;j EffectA[i*nFree+j]=EffectA[i*nFree+j]-a1*EffectA[k*nFree+j]; b[i]=b[i]-a1*b[k]; } } b[nFree-1]=b[nFree-1]/EffectA[(nFree-1)*nFree+(nFree-1)];//Base-0adjustment. for(i=nFree-2;i>-1;i--)//Base-0adjustment. { for(j=i+1;j b[i]=b[i]-EffectA[i*nFree+j]*b[j]; b[i]=b[i]/EffectA[i*nFree+i]; } deleteEffectA; //return(0); } //////////// doubleGetArea(intEleNum,int*EleNode,double*Xcoord,double*Ycoord) { doublexi,xj,xm,yi,yj,ym; doubleai,aj,am,bi,bj,bm,ci,cj,cm,area; intni=EleNode[EleNum*3+0]; intnj=EleNode[EleNum*3+1]; intnm=EleNode[EleNum*3+2]; /////////////////////对应单元的结点以及坐标 xi=Xcoord[ni];xj=Xcoord[nj];xm=Xcoord[nm]; yi=Ycoord[ni];yj=Ycoord[nj];ym=Ycoord[nm]; /////////////////////////// ai=xj*ym-xm*yj; aj=xm*yi-xi*ym; am=xi*yj-xj*yi; bi=yj-ym; bj=ym-yi; bm=yi-yj; ci=-(xj-xm); cj=-(xm-xi); cm=-(xi-xj); area=fabs(ai+aj+am)/2; returnarea; } voidGetStiffness(intEleNum,int*EleNode,int*ResNodeDOF,double*Xcoord,double*Ycoord,doubleElasticModule,doubleMiu,doubleThick,boolIsPlainStrain,doubleK[][6]) { ///////////////////////默认是平面应力问题,后面是平面应变问题 inti,j; doublexi,xj
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 弹性 力学 结点 三角形 单元 静力 分析 程序