有限元分析课程设计.docx

有限元分析
课程设计

学 院：土木建筑工程学院
专 业： 工程力学
班 级： 力学131
学 号： 3121631023
姓 名： 崔晨露

一、 单元划分与节点编号

二、 C语言程序
#include<stdio.h>
#include<math.h>

#define NE 32 //单元数
#define NJ 45 //节点数
#define NZ 10 //支承数
#define NPJ 7 //节点荷载数
#define NJ2 90 //节点位移数
#define DD 14 //半带宽

int LXM=0; //类型判别码
double EO=170e9; //杨氏模量
double MU=0.3; //泊松比
double LOU=0; //容重
double TE=0.1; //厚度

double AJZ [NJ+1][3]= {{0,0,0},{0,0,0},{0,0,0.75},{0,0,1.5},{0,0,2.25},{0,0,3},{0,1.875,0},{0,1.875,0.75},{0,1.875,1.5},{0,1.875,2.25},{0,1.875,3},{0,3.75,0},{0,3.75,0.75},{0,3.75,1.5},{0,3.75,2.25},{0,3.75,3},{0,5.625,0},{0,5.625,0.75},{0,5.625,1.5},{0,5.625,2.25},{0,5.625,3},{0,7.5,0},{0,7.5,0.75},{0,7.5,1.5},{0,7.5,2.25},{0,7.5,3},{0,9.375,0},{0,9.375,0.75},{0,9.375,1.5},{0,9.375,2.25},{0,9.375,3},{0,11.25,0},{0,11.25,0.75},{0,11.25,1.5},{0,11.25,2.25},{0,11.25,3},{0,13.125,0},{0,13.125,0.75},{0,13.125,1.5},{0,13.125,2.25},{0,13.125,3},{0,15,0},{0,15,0.75},{0,15,1.5},{0,15,2.25},{0,15,3}};//共36个节点

int JM [NE+1][5]= {{0,0,0,0},{0,1,6,7,2},{0,2,7,8,3},{0,3,8,9,4},{0,4,9,10,5},{0,6,11,12,7},{0,7,12,13,8},{0,8,13,14,9},{0,9,14,15,10},{0,11,16,17,12},{0,12,17,18,13},{0,13,18,19,14},{0,14,19,20,15},{0,16,21,22,17},{0,17,22,23,18},{0,18,23,24,19},{0,19,24,25,20},{0,21,26,27,22},{0,22,27,28,23},{0,23,28,29,24},{0,24,29,30,25},{0,26,31,32,27},{0,27,32,33,28},{0,28,33,34,29},{0,29,34,35,30},{0,31,36,37,32},{0,32,37,38,33},{0,33,38,39,34},{0,34,39,40,35},{0,36,41,42,37},{0,37,42,43,38},{0,38,43,44,39},{0,39,44,45,40}};//共32个单元

int NZC [NZ+1]={0,1,2,3,4,5,6,7,8,9,10};//1-5号节点的x,y被约束

double PJ [NPJ+1][2+1]={{0,0,0},{0,-3.75e4,50},{0,-7.5e4,60},{0,-7.5e4,70},{0,-7.5e4,80},{0,8e4,83},{0,8e4,87},{0,-3.75e4,90}};

double AE,JS[NJ*4][4],js[NJ+1][4],KZ[NJ2+1][DD+1],P[NJ2+1],S[3+1][8+1],KE[8+1][8+1],SZ[3+1][32+1];

int IE,JE,ME,LE;

void DUGD(int,int);//生成S矩阵，KE矩阵

void main()
{
int NJ1,k,IN,IM,jn,m,i,j,z,JO,ii,jj,h,dh,E,l,zl,dl,n;
double PE,c,SIG1,SIG2,SIG3,PYL,RYL,MAYL,MIYL,CETA;
double WY[8+1],YL[3+1];

if(LXM!=0)//平面应力问题与平面应变问题的判别
{
EO=EO/(1.0-MU*MU);
MU=MU/(1.0-MU);
}

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<=4;i++)
{
for(ii=1;ii<=2;ii++)
{
h=2*(i-1)+ii;
dh=2*(JM[E][i]-1)+ii;
for(j=1;j<=4;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];
}
}
}
}
}
//********形成P矩阵********
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*9.8/4;
P[2*IE]=P[2*IE]+PE;
P[2*JE]=P[2*JE]+PE;
P[2*ME]=P[2*ME]+PE;
P[2*LE]=P[2*LE]+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)
JO=DD;
else
JO=z;
for(j=2;j<=JO;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;
else
IM=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)
JO=NJ2-i+1;
else
JO=DD;
for(j=2;j<=JO;j++)
{
h=j+i-1;
P[i]=P[i]-KZ[i][j]*P[h];
}
P[i]=P[i]/KZ[i][1];
}
printf("\n");
printf("JD U V\n");
for(i=1;i<=NJ;i++)
printf("%d %-9.9f %-9.9f\n",i,P[2*i-1],P[2*i]);

//*********求单元应力和主应力********
for(i=0;i<=NJ*4;i++)
for(j=0;j<=3;j++)
JS[i][j]=0.0;
k=0;
for(E=1;E<=NE;E++)
{
DUGD(E,2);
for(i=1;i<=4;i++)
{
for(j=1;j<=2;j++)
{
h=2*(i-1)+j;
dh=2*(JM[E][i]-1)+j;
WY[h]=P[dh];
}
}
for(n=1;n<=4;n++)
{
for(i=1;i<=3;i++)
{
YL[i]=0;
for(j=1;j<=8;j++)
YL[i]=YL[i]+SZ[i][8*(n-1)+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));
//将节点应力值记入节点应力矩阵
k++;
JS[k][0]=JM[E][n];
JS[k][1]=SIG1;
JS[k][2]=SIG2;
JS[k][3]=SIG3;
printf("\n");
printf("E=%d JD=%d\n",E,JM[E][n]);
printf("sx=%-9.2f sy=%-9.2f tou=%-9.2f\n",SIG1,SIG2,SIG3);
printf("s1=%-9.2f s3=%-9.2f theta=%-9.2f\n",MAYL,MIYL,CETA);
}
}
for(i=0;i<=NJ;i++)
for(j=0;j<=3;j++)
js[i][j]=0.0;
printf("\n");
printf("\n");
printf("平均节点应力:\n");
for(ii=1;ii<=NJ;ii++)
{
k=0;
for(jj=1;jj<=NJ*4;jj++)
{
if(JS[jj][0]==ii)
{
k++;
js[ii][1]=js[ii][1]+JS[jj][1];
js[ii][2]=js[ii][2]+JS[jj][2];
js[ii][3]=js[ii][3]+JS[jj][3];
}
}
js[ii][1]=js[ii][1]/k;
js[ii][2]=js[ii][2]/k;
js[ii][3]=js[ii][3]/k;
js[ii][0]=sqrt(((js[ii][1]-js[ii][2])*(js[ii][1]-js[ii][2])+js[ii][1]*js[ii][1]+js[ii][2]*js[ii][2]+js[ii][3]*js[ii][3]*6)/2);
printf("JD=%d\n",ii);
printf("sx=%-9.2f sy=%-9.2f tou=%-9.2f\nsi=%-9.2f\n\n",js[ii][1],js[ii][2],js[ii][3],js[ii][0]);
}
scanf("%d",&l);
}

/**************DUGD子程序******************/
void DUGD(int E,int ASK)
{
double B[3+1][8+1],D[3+1][3+1],x,y,a,b,C;
int i, j, k,n,w;
double xy[5][3]={{0,0,0},{0,-1,-1},{0,1,-1},{0,1,1},{0,-1,1}};

a=0.9375;//单元半长度
b=0.375; //单元半宽度

AE=4*a*b;

D[1][1]=EO/(1-MU*MU);
D[1][2]=MU*EO/(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));

if(ASK>1)
{
for(w=1;w<=4;w++)
{
for(i=0;i<=3;i++)
for(j=0;j<=8;j++)
B[i][j]=0.0;

x=xy[w][1]*a;
y=xy[w][2]*b;

B[1][1]=(y-b)/AE;
B[1][3]=(b-y)/AE;
B[1][5]=(b+y)/AE;
B[1][7]=-(b+y)/AE;
B[2][2]=(x-a)/AE;
B[2][4]=-(a+x)/AE;
B[2][6]=(a+x)/AE;
B[2][8]=(a-x)/AE;
B[3][1]=(x-a)/AE;
B[3][2]=(y-b)/AE;
B[3][3]=-(a+x)/AE;
B[3][4]=(b-y)/AE;
B[3][5]=(a+x)/AE;
B[3][6]=(b+y)/AE;
B[3][7]=(a-x)/AE;
B[3][8]=-(b+y)/AE;

for(i=1;i<=3;i++)
{
for(j=1;j<=8;j++)
{
S[i][j]=0;
for(k=1;k<=3;k++)
S[i][j]=S[i][j]+D[i][k]*B[k][j];
n=8*(w-1)+j;
SZ[i][n]=S[i][j];

}
}

}

}
if(ASK>2)
{
C=EO*TE/(1-MU*MU);

KE[1][1]=C*(b/(3*a)+(1-MU)*a/(6*b));
KE[1][2]=C*(1+MU)/8;
KE[1][3]=C*((-b)/(3*a)+(1-MU)*a/(12*b));
KE[1][4]=C*((3*MU-1)/8);
KE[1][5]=C*((-b)/(6*a)-(1-MU)*a/(12*b));
KE[1][6]=(-C)*(1+MU)/8;
KE[1][7]=C*(b/(6*a)-(1-MU)*a/(6*b));
KE[1][8]=C*((1-3*MU)/8);
KE[2][1]=KE[1][2];
KE[2][2]=C*(a/(3*b)+(1-MU)*b/(6*a));
KE[2][3]=KE[1][8];
KE[2][4]=C*(a/(6*b)-(1-MU)*b/(6*a));
KE[2][5]=KE[1][6];
KE[2][6]=C*((-a)/(6*b)-(1-MU)*b/(12*a));
KE[2][7]=KE[1][4];
KE[2][8]=C*((-a)/(3*b)+(1-MU)*b/(12*a));
KE[3][1]=KE[1][3];
KE[3][2]=KE[2][3];
KE[3][3]=KE[1][1];
KE[3][4]=KE[1][6];
KE[3][5]=KE[1][7];
KE[3][6]=KE[1][4];
KE[3][7]=KE[1][5];
KE[3][8]=KE[1][2];
KE[4][1]=KE[1][4];
KE[4][2]=KE[2][4];
KE[4][3]=KE[3][4];
KE[4][4]=KE[2][2];
KE[4][5]=KE[2][3];
KE[4][6]=KE[2][8];
KE[4][7]=KE[1][2];
KE[4][8]=KE[2][6];
KE[5][1]=KE[1][5];
KE[5][2]=KE[2][5];
KE[5][3]=KE[3][5];
KE[5][4]=KE[4][5];
KE[5][5]=KE[1][1];
KE[5][6]=KE[1][2];
KE[5][7]=KE[1][3];
KE[5][8]=KE[1][4];
KE[6][1]=KE[1][6];
KE[6][2]=KE[2][6];
KE[6][3]=KE[3][6];
KE[6][4]=KE[4][6];
KE[6][5]=KE[5][6];
KE[6][6]=KE[2][2];
KE[6][7]=KE[1][8];
KE[6][8]=KE[2][4];
KE[7][1]=KE[1][7];
KE[7][2]=KE[2][7];
KE[7][3]=KE[3][7];
KE[7][4]=KE[4][7];
KE[7][5]=KE[5][7];
KE[7][6]=KE[6][7];
KE[7][7]=KE[1][1];
KE[7][8]=KE[1][6];
KE[8][1]=KE[1][8];
KE[8][2]=KE[2][8];
KE[8][3]=KE[3][8];
KE[8][4]=KE[4][8];
KE[8][5]=KE[5][8];
KE[8][6]=KE[6][8];
KE[8][7]=KE[7][8];
KE[8][8]=KE[2][2];
}

}
三、C语言运行结果
JD U V
1 0.000000000 0.000000000
2 0.000000000 0.000000000
3 0.000000000 0.000000000
4 0.000000000 0.000000000
5 0.000000000 0.000000000
6 -0.000188086 -0.000164909
7 -0.000085891 -0.000145676
8 0.000005546 -0.000140236
9 0.000097076 -0.000147161
10 0.000199777 -0.000168014
11 -0.000342848 -0.000530025
12 -0.000162895 -0.000518568
13 0.000011559 -0.000514919
14 0.000186049 -0.000520018
15 0.000365862 -0.000532873
16 -0.000462248 -0.001078026
17 -0.000219051 -0.001068957
18 0.000017352 -0.001066485
19 0.000253834 -0.001070229
20 0.000497292 -0.001080578
21 -0.000545495 -0.001758929
22 -0.000258019 -0.001753258
23 0.000023292 -0.001752058
24 0.000304770 -0.001755558
25 0.000592898 -0.001763511
26 -0.000595214 -0.002527541
27 -0.000280247 -0.002524714
28 0.000029780 -0.002524792
29 0.000339881 -0.002527999
30 0.000655116 -0.002534183
31 -0.000617969 -0.003343643
32 -0.000289101 -0.003342750
33 0.000036320 -0.003343443
34 0.000361985 -0.003345918
35 0.000690729 -0.003349823
36 -0.000622391 -0.004177818
37 -0.000289316 -0.004178483
38 0.000043378 -0.004179828
39 0.000374731 -0.004182063
40 0.000708887 -0.004185090
41 -0.000620797 -0.005016007
42 -0.000282489 -0.005016335
43 0.000047901 -0.005017224
44 0.000385153 -0.005018842
45 0.000715651 -0.005020841

E=1 JD=1
sx=-18739706.59 sy=-5621911.98 tou=-5750681.88
s1=-3457881.11 s3=-20903737.46 theta=110.62

E=1 JD=6
sx=-17302482.56 sy=-831165.23 tou=3158628.43
s1=-246222.60 s3=-17887425.19 theta=79.51

E=1 JD=7
sx=-7120413.64 sy=2223455.45 tou=3829332.97
s1=3592279.42 s3=-8489237.61 theta=70.33

E=1 JD=2
sx=-8557637.66 sy=-2567291.30 tou=-5079977.34
s1=334757.93 s3=-11459686.89 theta=119.74

E=2 JD=2
sx=-8557637.66 sy=-2567291.30 tou=-5079977.34
s1=334757.93 s3=-11459686.89 theta=119.74

E=2 JD=7
sx=-8151125.39 sy=-1212250.39 tou=2891482.49
s1=-165310.18 s3=-9198065.59 theta=70.10

E=2 JD=8
sx=959114.42 sy=1520821.56 tou=3081188.22
s1=4333929.78 s3=-1853993.81 theta=47.60

E=2 JD=3
sx=552602.14 sy=165780.64 tou=-4890271.61
s1=5253286.21 s3=-4534903.43 theta=136.13

E=3 JD=3
sx=552602.14 sy=165780.64 tou=-4890271.61
s1=5253286.21 s3=-4534903.43 theta=136.13

E=3 JD=8
sx=35129.71 sy=-1559127.47 tou=3089199.22
s1=2428387.59 s3=-3952385.35 theta=37.77

E=3 JD=9
sx=9154524.94 sy=1176691.10 tou=2847712.09
s1=10066722.42 s3=264493.63 theta=17.76

E=3 JD=4
sx=9671997.38 sy=2901599.21 tou=-5131758.74
s1=12434523.13 s3=139073.46 theta=151.71

E=4 JD=4
sx=9671997.38 sy=2901599.21 tou=-5131758.74
s1=12434523.13 s3=139073.46 theta=151.71

E=4 JD=9
sx=8113710.99 sy=-2292688.73 tou=3821656.61
s1=9366388.14 s3=-3545365.88 theta=18.15

E=4 JD=10
sx=18346185.68 sy=777053.67 tou=3094456.29
s1=18875279.52 s3=247959.84 theta=9.70

E=4 JD=5
sx=19904472.06 sy=5971341.62 tou=-5858959.06
s1=22040677.42 s3=3835136.26 theta=159.97

E=5 JD=6
sx=-13982262.17 sy=164900.89 tou=-3822943.62
s1=1131869.94 s3=-14949231.21 theta=104.19

E=5 JD=11
sx=-14563359.43 sy=-1772089.99 tou=2955975.92
s1=-1122021.26 s3=-15213428.16 theta=77.60

E=5 JD=12
sx=-6816022.82 sy=552110.99 tou=2684797.20
s1=1426605.81 s3=-7690517.63 theta=71.96

E=5 JD=7
sx=-6234925.55 sy=2489101.88 tou=-4094122.34
s1=4109479.11 s3=-7855302.78 theta=111.59

E=6 JD=7
sx=-7265637.30 sy=-946603.96 tou=-5031972.82
s1=1835537.06 s3=-10047778.32 theta=118.94

E=6 JD=12
sx=-7399471.66 sy=-1392718.48 tou=2205311.70
s1=-670015.27 s3=-8122174.86 theta=71.86

E=6 JD=13
sx=871710.65 sy=1088636.22 tou=2142855.67
s1=3125772.32 s3=-1165425.45 theta=46.45

E=6 JD=8
sx=1005545.01 sy=1534750.73 tou=-5094428.85
s1=6371443.79 s3=-3831148.05 theta=133.51

E=7 JD=8
sx=81560.30 sy=-1545198.29 tou=-5086417.85
s1=4419222.89 s3=-5882860.88 theta=139.54

E=7 JD=13
sx=218004.27 sy=-1090385.06 tou=2146082.96
s1=1807387.71 s3=-2679768.51 theta=36.52

E=7 JD=14
sx=8483719.49 sy=1389329.50 tou=2209756.81
s1=9115713.11 s3=757335.88 theta=15.96

E=7 JD=9
sx=8347275.52 sy=934516.27 tou=-5022744.00
s1=10883107.62 s3=-1601315.83 theta=153.21

E=8 JD=9
sx=7306461.57 sy=-2534863.56 tou=-4048799.48
s1=8758058.99 s3=-3986460.98 theta=160.28

E=8 JD=14
sx=7904177.59 sy=-542476.81 tou=2673789.95
s1=8679414.72 s3=-1317713.94 theta=16.17

E=8 JD=15
sx=15587136.94 sy=1762410.99 tou=2952724.10
s1=16191378.35 s3=1158169.58 theta=11.57

E=8 JD=10
sx=14989420.92 sy=-229975.76 tou=-3769865.34
s1=15872036.37 s3=-1112591.20 theta=166.82

E=9 JD=11
sx=-11040187.35 sy=-715138.37 tou=-3421533.32
s1=315765.06 s3=-12071090.78 theta=106.77

E=9 JD=16
sx=-11218646.41 sy=-1310001.90 tou=2092094.30
s1=-886390.77 s3=-11642257.55 theta=78.55

E=9 JD=17
sx=-4917357.71 sy=580384.71 tou=2008813.40
s1=1236160.94 s3=-5573133.94 theta=71.92

E=9 JD=12
sx=-4738898.65 sy=1175248.24 tou=-3504814.21
s1=2803805.18 s3=-6367455.59 theta=114.92

E=10 JD=12
sx=-532234