结构及有限元教学课件：平面问题有限元-满阵.docVIP

平面问题有限元程序（部分）--满阵存储总体刚度矩阵 （形成总体刚度矩阵、解方程部分） //输入原始数据 （自编） //总刚矩阵的清零 for(i5 = 0; i5 NJF; i5++) for(i6 = 0;i6 NJF; i6++) { zk[i5][i6]=0; } //总刚生成程序 for(je=0;jeNE;je++) { //变量je表示单元编号 i1=ie[je][0]-1; i2=ie[je][1]-1; i3=ie[je][2]-1;// //z矩阵表示结点坐标矩阵 b[je][0]=z[i2][1]-z[i3][1];//bi=yj-ym b[je][1]=z[i3][1]-z[i1][1];//bj=ym-yi b[je][2]=z[i1][1]-z[i2][1];//bm=yi-yj c[je][0]=z[i3][0]-z[i2][0];//ci=xm-xj c[je][1]=z[i1][0]-z[i3][0];//cj=xi-xm c[je][2]=z[i2][0]-z[i1][0];//cm=xj-xi //aa向量表示面积,由b,c 矩阵的元素计算而得 aa[je]=0.5*(b[je][0]*c[je][1]-b[je][1]*c[je][0]); //a向量为中间量, E为弹性模量，T为板厚度，VM为泊松比 a[je]=E*T/(4.0*(1.0-VM*VM)*aa[je]); for(i=0;i3;i++) { for(j=0;j3;j++) { //计算单刚子块的四个元素，准备向总刚对应位置填写 g[0][0]=b[je][i]*b[je][j]+0.5*(1.0-VM)*c[je][i]*c[je][j]; g[1][1]=c[je][i]*c[je][j]+0.5*(1.0-VM)*b[je][i]*b[je][j]; g[1][0]=VM*b[je][j]*c[je][i]+0.5*(1.0-VM)*b[je][i]*c[je][j]; g[0][1]=VM*c[je][j]*b[je][i]+0.5*(1.0-VM)*c[je][i]*b[je][j]; for(i7 = 0;i7 2;i7++) { for(j7 = 0;j7 2;j7++) { i4=2*(ie[je][i]-1)+i7; j4=2*(ie[je][j]-1)+j7; zk[i4][j4]+=g[i7][j7]*a[je];//总刚元素的计算和填写 } } } } } //将载荷和约束条件引入 （自编） //高斯消元法解方程，求出节点位移 for(k = 0;k NJF-1;k++) { for(i = k+1;i NJF;i++) { cc=zk[i][k]/zk[k][k]; for(j = 0;j NJF;j++) zk[i][j]-=cc*zk[k][j]; P[i]=P[i]-cc*P[k];//计算前，p矩阵是边界条件向量，待输入；然后将线性方程组的系数矩阵转变为上三角矩阵 } } P[NJF-1]/=zk[NJF-1][NJF-1];//解最后一个方程 for(i = NJF-2;i = 0;i--) { for(j = i+1;j NJF;j++) P[i]-=zk[i][j]*P[j]; P[i]/=zk[i][i];//回带，解其它方程，最后解得的未知数放在P矩阵，得到节点位移向量 } printf(输出位移向量P:\n);