梯形水坝的有限元分析.docx

梯形水坝的有限元分析.docx

《梯形水坝的有限元分析.docx》由会员分享，可在线阅读，更多相关《梯形水坝的有限元分析.docx（20页珍藏版）》请在冰豆网上搜索。

梯形水坝的有限元分析

一题目

下图为一梯形水坝，载荷沿厚度均匀分布，其尺寸如图所示，厚度t=1m,弹性模量E=79G,泊松μ=0.3，容重ρ=3200.

6m

9m

二划分单元及载荷处理

6m

9m

1.5m

1.5m

1.5m

1.5m

1.5m

1.5m

三输入数据

1.基本参数

⑴单元数NE=40;

⑵节点数NJ=30;

⑶支撑数NZ=14;

⑷节点载荷数NPJ=7;

⑸半宽带DD=14；

⑹节点位移数NJ2=NJ*2=60;

2.其他参数

⑴问题类型码LXM,LXM=

⑵弹性常数E,弹性模量EO=79，泊松比μ,MU=0.3;

⑶容重ρ,LOU=3200;

⑷板厚t,TE=1;

⑸节点坐标数组AJZ.

AJZ（NJ,2）=

AJZ（NJ,2）中的元素按节点整体码顺序输入，数组行号为节点整体码，每个节点的坐标值存一行，第一列存x的值.如

=AJZ（2,1）

=AJZ（4,2）.

⑹节点码数组JM

JM（NE,3）=

JM（NE,3）中的元素按单元输入。

单元号为行号，每个单元的整体码存一行，局部码为列号，局部码对应的整体码存放在相应的列。

⑺支撑数组NZC.

NZC（NZ）=

NZC（NZ）中的元素按支撑对应的位移数从小到大排列。

⑻节点载荷数组PJ

PJ（NPJ+1,2+1）=

PJ（NPJ,2）中的元素排列原则为：

一个节点载荷存一行，同一行中第一列是节点载荷值，第二列是载荷对应的位移数。

四程序及结果

#include

#include

#defineNE40

#defineNJ30

#defineNZ10

#defineNPJ7

#defineNJ260

#defineDD14

intLXM=1;

doubleEO=79000000000;

doubleMU=0.3;

doubleLOU=3200;

doubleTE=1.0;

//节点坐标数组

doubleAJZ[NJ+1][3]={{0,0,0},{0,0,0},{0,1.5,0},{0,3,0},{0,4.5,0},{0,6,0},\

{0,0,1.5},{0,1.5,1.5},{0,3,1.5},{0,4.5,1.5},{0,5.25,1.5},{0,0,3},{0,1.5,3},\

{0,3,3},{0,4.5,3},{0,4.85,3},{0,0,4.5},{0,1.5,4.5},{0,3,4.5},{0,4.5,4.5},\

{0,0,6},{0,1.5,6},{0,3,6},{0,4.275,6},{0,0,7.5},{0,1.5,7.5},{0,3,7.5},\

{0,3.75,7.5},{0,0,9},{0,1.5,9},{0,3,9}};

//单元节点码数组

intJM[NE+1][4]={{0,0,0,0},{0,1,7,6},{0,1,2,7},{0,2,8,7},{0,2,3,8},{0,3,9,8},\

{0,3,4,9},{0,4,10,9},{0,4,5,10},{0,6,12,11},{0,6,7,12},{0,7,13,12},{0,7,8,13},\

{0,8,14,13},{0,8,9,14},{0,9,15,14},{0,9,10,15},{0,11,17,16},{0,11,12,17},\

{0,12,18,17},{0,12,13,18},{0,13,19,18},{0,13,14,19},{0,14,15,19},{0,16,21,20},\

{0,16,17,21},{0,17,22,21},{0,17,18,22},{0,18,23,22},{0,18,19,23},{0,20,25,24},\

{0,20,21,25},{0,21,26,25},{0,21,22,26},{0,22,27,26},{0,22,23,27},{0,24,29,28},\

{0,24,25,29},{0,25,30,29},{0,25,26,30},{0,26,27,30}};

intNZC[NZ+1]={0,1,2,3,4,5,6,7,8,9,10};

doublePJ[NPJ+1][2+1]={{0,0,0},{0,207637.5,1},{0,339937.5,11},{0,229687.5,21},

{0,141487.5,31},{0,75337.5,41},{0,31237.5,47},{0,4900,55}};

doubleAE,KZ[NJ2+1][DD+1],P[NJ2+1],S[3+1][6+1],KE[6+1][6+1];

intIE,JE,ME;

voidDUGD（int,int）;

voidmain（）

{

intNJ1,k,IN,IM,jn,m,i,j,z,J0,ii,jj,h,dh,E,l,zl,dl;

doublePE,c,SIG1,SIG2,SIG3,PYL,RYL,MAYL,MIYL,CETA;

doubleWY[6+1],YL[3+1];

//?

÷?

?

?

?

?

?

DUGD?

¨?

?

?

?

ú?

?

KZ?

?

?

ó:

for（i=0;i<=NJ2;i++）

{

for（j=0;j<=DD;j++）

KZ[i][j]=0.0;

}

for（E=1;E<=NE;E++）

{

DUGD（E,3）;

for（i=1;i<=3;i++）

{

for（ii=1;ii<=2;ii++）

{

h=2*（i-1）+ii;

dh=2*（JM[E][i]-1）+ii;

for（j=1;j<=3;j++）

{

for（jj=1;jj<=2;jj++）

{

l=2*（j-1）+jj;

zl=2*（JM[E][j]-1）+jj;

dl=zl-dh+1;

if（dl>0）

KZ[dh][dl]=KZ[dh][dl]+KE[h][l];

}

}

}

}

}

for（i=1;i<=NJ2;i++）

P[i]=0.0;

if（NPJ>0）

{

for（i=1;i<=NPJ;i++）

{

j=（int）PJ[i][2];

P[j]=PJ[i][1];

}

}

if（LOU>0）

{

for（E=1;E<=NE;E++）

{

DUGD（E,1）;

PE=-LOU*（AE）*TE/3;

P[2*IE]=P[2*IE]+PE;

P[2*JE]=P[2*JE]+PE;

P[2*ME]=P[2*ME]+PE;

}

}

for（i=1;i<=NZ;i++）

{

z=NZC[i];

KZ[z][1]=1.0;

for（j=2;j<=DD;j++）

KZ[z][j]=0.0;

if（z!

=1）

{

if（z>DD）J0=DD;

elseJ0=z;

for（j=2;j<=J0;j++）

KZ[z-j+1][j]=0.0;

}

P[z]=0.0;

}

NJ1=NJ2-1;

for（k=1;k<=NJ1;k++）

{

if（NJ2>k+DD-1）

IM=k+DD-1;

elseIM=NJ2;

IN=k+1;

for（i=IN;i<=IM;i++）

{

l=i-k+1;

c=KZ[k][l]/KZ[k][1];

jn=DD-l+1;

for（j=1;j<=jn;j++）

{

m=j+i-k;

KZ[i][j]=KZ[i][j]-c*KZ[k][m];

}

P[i]=P[i]-c*P[k];

}

}

P[NJ2]=P[NJ2]/KZ[NJ2][1];

for（i=NJ1;i>=1;i--）

{

if（DD>NJ2-i+1）J0=NJ2-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（"JDUV\n"）;

for（i=1;i<=NJ;i++）

{

printf（"%2d%-9.6f%9.6f\n",i,P[2*i-1],P[2*i]）;

}

for（E=1;E<=NE;E++）

{

DUGD（E,2）;

for（i=1;i<=3;i++）

{

for（j=1;j<=2;j++）

{

h=2*（i-1）+j;

dh=2*（JM[E][i]-1）+j;

WY[h]=P[dh];

}

}

for（i=1;i<=3;i++）

{

YL[i]=0;

for（j=1;j<=6;j++）

YL[i]=YL[i]+S[i][j]*WY[j];

}

SIG1=YL[1];

SIG2=YL[2];

SIG3=YL[3];

PYL=（SIG1+SIG2）/2;

RYL=sqrt（pow（（SIG1-SIG2）/2.0,2）+pow（SIG3,2））;

MAYL=PYL+RYL;

MIYL=PYL-RYL;

if（SIG2==MIYL）

CETA=0;

else

CETA=90-57.29578*atan2（SIG3,（SIG2-MIYL））;

printf（"\n"）;

printf（"E=%2d\n",E）;

printf（"sx=%-18.6fsy=%-18.6ftou=%-18.6f\n",SIG1,SIG2,SIG3）;

printf（"s1=%-18.6fs2=%-18.6ftheta=%-18.6f\n",MAYL,MIYL,CETA）;

}

}

//DUGD?

?

?

?

:

voidDUGD（intE,intASK）

{

if（LXM!

=0）

{

EO=EO/（1.0-MU*MU）;

MU=MU/（1.0-MU）;

}

doubleB[3+1][6+1],D[3+1][3+1];

inti,j,k;

doubleCM,BM,CJ,BJ;

IE=JM[E][1];

JE=JM[E][2];

ME=JM[E][3];

CM=AJZ[JE][1]-AJZ[IE][1];

BM=AJZ[IE][2]-AJZ[JE][2];

CJ=AJZ[IE][1]-AJZ[ME][1];

BJ=AJZ[ME][2]-AJZ[IE][2];

AE=（BJ*CM-BM*CJ）/2.0;

if（ASK>1）

{

for（i=1;i<=3;i++）

for（j=1;j<=6;j++）

B[i][j]=0.0;

B[1][1]=（-BJ-BM）/（2.0*AE）;

B[1][3]=BJ/（2*AE）;

B[1][5]=BM/（2*AE）;

B[2][2]=（-CJ-CM）/（2*AE）;

B[2][4]=CJ/（2*AE）;

B[2][6]=CM/（2*AE）;

B[3][1]=B[2][2];

B[3][2]=B[1][1];

B[3][3]=B[2][4];

B[3][4]=B[1][3];

B[3][5]=B[2][6];

B[3][6]=B[1][5];

D[1][1]=EO/（1-MU*MU）;

D[1][2]=EO*MU/（1-MU*MU）;

D[1][3]=0;

D[2][1]=D[1][2];

D[2][2]=D[1][1];

D[2][3]=0;

D[3][1]=0;

D[3][2]=0;

D[3][3]=EO/（2*（1+MU））;

for（i=1;i<=3;i++）

{

for（j=1;j<=6;j++）

{

S[i][j]=0.0;

for（k=1;k<=3;k++）

S[i][j]=S[i][j]+D[i][k]*B[k][j];

}

}

if（ASK>2）

{

for（i=1;i<=6;i++）

{

for（j=1;j<=6;j++）

{

KE[i][j]=0.0;

for（k=1;k<=3;k++）

KE[i][j]=KE[i][j]+S[k][i]*B[k][j]*AE*TE;

}

}

}

}

}

五.结果分析：

1.输出变形图：

由运算结果中的u,v的即是各个节点在u,v方向的变形位移图，由CAD输出图形时，由于变形的位移很小，为了更清楚的显示出变形，应把u，v方向的位移相应的扩大一定的倍数，如图示：

6m

9m

1.5m

1.5m

1.5m

1.5m

1.5m

1.5m

2.输出等应力线图：

又运行结果中各个单元的s1,s2,tou计算出各个节点相应的s1,s2,tou，再次计算出给节点的主应力s1。

由c语言编程以计算各个节点的主应力的程序：

#include

#include

#defineMAX10

voidmain（）

{

intt;

inti;//循环变量

floata[MAX][3];

floatzx,zy,zz;

floatpjx,pjy,pjz;

floatPYL,RYL,MAYL;

zx=0;

zy=0;

zz=0;

printf（"输入该节点周围的单元数：

\n"）;

printf（"该节点周围的单元数t="）;

scanf（"%d",&t）;

printf（"输入该节点周围的单元的应力：

\n"）;

for（i=1;i<=t;i++）

scanf（"%f%f%f",&a[i][0],&a[i][1],&a[i][2]）;

for（i=1;i<=t;i++）

{

zx=zx+a[i][0];

zy=zy+a[i][1];

zz=zz+a[i][2];

}

pjx=zx/t;

pjy=zy/t;

pjz=zz/t;

PYL=（pjx+pjy）/2;

RYL=（float）sqrt（pow（（pjx-pjy）/2,2）+pow（pjz,2））;

MAYL=PYL+RYL;

printf（"s1=%-9.6f",MAYL）;

}//输出结果即是各个节点的主应力

计算各节点的主应力的结果

1.166733.1210942.133281.6718753.35427.1484384.8616.34375

5.34336.9531256.167154.3125007.118257.2031258.45992.125000

9.-1349.51562510.-8356.39062511.120053.78906312.55704.273438

13.33174.42968814.7007.82812515.-13444.14062516.11800.183594

17.19107.16796918.19210.58203119.996.2409420.-604.39113

21.4432.1758122.239223.-1597.296875

24.847.91992225.-1179.40820326.-4654.45459027.-4796.332031

28.-2451.62890629.-2070.97802730.-2764.493652

找出相应的应力相等的点，由CAD软件绘制等应力线的图，并且输出，如下图所示：