弹性力学三结点三角形单元静力分析程序.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;jGlobalK[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;jfout<fout<}
GetLoadVector（LoadNodeNum,SurLoadNum,Density,ElementNum,Thick,EleNode,NodeDOF,Xcoord,Ycoord,NodeLoad,SurLoad,LoadorDisp）;
fout<<"结点等效荷载列阵"<for（i=0;i{
fout<}
fout<for（i=FreeDOF;igauss（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;ifor（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;jEffectA[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;jb[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
升级会员 