平面四节点等参单元分析程序
四边形四结点等参元matlab程序广州大学
%---------------四边形四节点等参元matlab计算程序----------------------------% 2013年% 13级建筑与土木工程PHILDER% B15-404% 本程序只输出了节点位移、单元中心点的应力%******************************************************************* clear all; %清楚内存变量clc; %清屏format short e %设定数据输出类型FP1=fopen('input.txt','rt');%打开初始数据文件NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); % 总结点数NFPOIN=fscanf(FP1,'%d',1); %受力结点数NVFIX=fscanf(FP1,'%d',1); %约束结点个数E=fscanf(FP1,'%f',1); %弹性模量v=fscanf(FP1,'%f',1); % 泊松比t=fscanf(FP1,'%f',1); %厚度LNODS=fscanf(FP1,'%f',[4,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号x,y坐标(整体坐标下)FPOIN=fscanf(FP1,'%f',[3,NFPOIN])'; % 节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号(n,2),(n,3)约束对应的位移编码HK=zeros(2*NPOIN ,2*NPOIN); %生成总刚度矩阵并清零FORCE=zeros(2*NPOIN,1); %生成总荷载列阵并清零for i=1:NELEM %对单元个数进行循环EK=ele_EK(i, E, v, LNODS, COORD ,t); %生成单刚矩阵for j=1:4 %对行循环N1=LNODS(i,j)*2; %j结点对应的第二个位移编码for k=1:4 %对列循环N2=LNODS(i,k)*2; %k结点对应的第二个位移编码HK(N1-1:N1,N2-1:N2)= HK(N1-1:N1,N2-1:N2)+EK(j*2-1:j*2,k*2-1:k*2);endendendfor i=1:NFPOIN %对作用荷载的结点数进行循环N1=FPOIN(i,1)*2-2; %作用荷载的节点对应的第一个位移编码-1for j=1:2 %X, Y方向FORCE(N1+j)=FPOIN(i, j+1);endendfor i=1:NVFIX %对约束数进行循环N1=FIXED(i); %约束对应的位移编码HK(:,N1)=0;HK(N1,:)=0;HK(N1,N1)=1;FORCE(N1,1)=0; %将零位移约束对应的行、列变成零，主元变成1endDISP=HK\FORCE; %方程求解，得结点位移a=E/(1-v*v);D=[a,a*v,0;a*v,a,0;0,0,a*(1-v)/2];%弹性矩阵B=ele_B(i, -1, -1, COORD,LNODS); %计算B矩阵，并赋值（-1，-1）用以验算EDISP=zeros(8,1); %生成单元位移列向量，并清零S=D*B; %计算应力矩阵Sfor i=1:NELEM %对单元数进行循环for j=1:4 %对结点数进行循环N1=LNODS(i, j)*2; %结点所对应的最后一个位移编码EDISP(2*j-1:2*j)=DISP(N1-1:N1);endFE=S*EDISP; %计算结点力FE(:,i)=FE; %将同一单元的节点力集合endFEfunction B=ele_B(i, s, t, COORD,LNODS)%(3*4)J=ele_J(i,s,t,LNODS,COORD);dN=ele_dN(s,t);A=[J(2,2) -J(1,2);-J(1,2) J(1,1)]*[J(1,1)*J(2,2)-J(1,2)*J(2,1)]*dN; % 形函数在整体坐标系下的导数N_x=A(1,1:4);N_y=A(2,1:4);B=zeros(3,4);for j=1:4B(1:3,2*j-1:2*j)=[N_x(1,j),0;0,N_y(1,j);N_y(1,j), N_x(1,j)];endreturnfunction dN=ele_dN(s,t)dN=zeros(2,4);dN1s=- (1-t) *0.25;dN1t= - (1-s) *0.25;dN2s = (1-t) *0.25;dN2t = (1+s) *0.25;dN3s = (1+t) *0.25;dN3t = (1+s) *0.25;dN4s= - (1+t) *0.25;dN4t = - (1-s) *0.25;returnfunction EK=ele_EK(i, E, v, LNODS, COORD ,t) EK=zeros(8,8);a=E/(1-v*v);D=[a,a*v,0;a*v,a,0;0,0,a*(1-v)/2];%弹性矩阵%高斯积分采用3 x 3 个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W=[W1 W2 W3];r=15^(0.5)/5;x=[-r 0 r];%积分点for j=1:3for k=1:3B=ele_B(i,x(j),x(k),COORD,LNODS);J=ele_J(i,x(j),x(k),LNODS,COORD);EK=EK+W(j)*W(k)*B'*D*B*det(J)*t;endendreturnfunction J=ele_J(i,s,t,LNODS,COORD)x=zeros(4,1);y=zeros(4,1);for j=1:4N=LNODS(i,j);%取结点号x(j,1)=COORD(N,1);y(j,1)=COORD(N,2);enddN=ele_dN(s,t);x_s=dN(1,1:4)*x;x_t=dN(2,1:4)*x;y_s=dN(1,1:4)*y;y_t=dN(2,1:4)*y;J=[ x_s, y_s; x_t, y_t]; Returnfunction N=shape(s,t) %ξ,ηN(1) = (1-s)*(1-t)*0.25;N(2) = (1+s)*(1-t)*0.25;N(3) = (1+s)*(1+t)*0.25;N(4) = (1-s)*(1+t)*0.25; returninput.3 8 1 2 3E11 0.3 11 2 3 42 5 8 34 6 7 80 01 01 10 12 02 13 03 14 0 1001 1 14 1 1。
有限元考试精彩试题及问题详解——第一组
有限元考试试题及答案一、简答题（5道，共计25分）。
1.有限单元位移法求解弹性力学问题的基本步骤有哪些？（5分）答：（1）选择适当的单元类型将弹性体离散化；（2）建立单元体的位移插值函数；（3）推导单元刚度矩阵；（4）将单元刚度矩阵组装成整体刚度矩阵；（5）代入边界条件和求解。
2. 在划分网格数相同的情况下，为什么八节点四边形等参数单元精度大于四边形矩形单元？（5分）答：在对于曲线边界的边界单元，其边界为曲边，八节点四边形等参数单元边上三个节点所确定的抛物线来代替原来的曲线，显然拟合效果比四边形矩形单元的直边好。
3.轴对称单元与平面单元有哪些区别？（5分）答：轴对称单元是三角形或四边形截面的空间的环形单元，平面单元是三角形或四边形平面单元；轴对称单元内任意一点有四个应变分量，平面单元内任意一点非零独立应变分量有三个。
4.有限元空间问题有哪些特征？（5分）答：（1）单元为块体形状。
常用单元：四面体单元、长方体单元、直边六面体单元、曲边六面体单元、轴对称单元。
（2）结点位移3个分量。
（3）基本方程比平面问题多。
3个平衡方程，6个几何方程，6个物理方程。
5.简述四节点四边形等参数单元的平面问题分析过程。
（5）分）答：（1）通过整体坐标系和局部坐标系的映射关系得到四节点四边形等参单元的母单元，并选取单元的唯一模式；（2）通过坐标变换和等参元确定平面四节点四边形等参数单元的几何形状和位移模式；（3）将四节点四边形等参数单元的位移模式代入平面问题的几何方程，得到单元应变分量的计算式，再将单元应变代入平面问题的物理方程，得到平面四节点等参数单元的应力矩阵；（4）用虚功原理求得单元刚度矩阵，最后用高斯积分法计算完成。
二、论述题（3道,共计30分）。
1. 简述四节点四边形等参数单元的平面问题分析过程。
（10分）答：（1）通过整体坐标系和局部坐标系的映射关系得到四节点四边形等参单元的母单元，并选取单元的唯一模式；（2） 通过坐标变换和等参元确定平面四节点四边形等参数单元的几何形状和位移模式；（3）将四节点四边形等参数单元的位移模式代入平面问题的几何方程，得到单元应变 分量的计算式，再将单元应变代入平面问题的物理方程，得到平面四节点等参数单元的应力矩阵；（4）用虚功原理求得单元刚度矩阵，最后用高斯积分法计算完成。
平面四节点等参单元有限元程序
}
//控制信息读入
fscanf(fp1,"%d",&np);
fscanf(fp1,"%d",&ne);
fscanf(fp1,"%d",&nr1);
fscanf(fp1,"%d",&nr2);
fscanf(fp1,"%d",&nd);
exit(1);
}
fprintf(fp2,"位移为：\n"); //输出位移
for(i=0;i<n;i+=2)
fprintf(fp2,"u%d=%f v%d=%f\n",i/2+1,p[i],i/2+1,p[i+1]);
fprintf(fp2,"\n各单元应力：Sx,\tSy,\tSxy,\tS1,\tS2,\tSmises\n"); //输出单元应力
float *x,*y,*p,**pp,u1,u2; //x,y,me,nrr,p为输入信息
float coords[2][4],det,fun[4],pn[2][4],xjac[2][2],rjac[2][2],dnx[2][4]; //求单刚定义的变量
int *LD,**me,*nrr,*nu,*IS,nn; //nn为变带宽一维组中元素个数
void output() //结果输出函数
{ int i,j;
FILE *fp2;
fp2=fopen("output.dat","w"); //运算结果放在output.dat中
平面四边形四节点等参单元Fortran源程序
C ************************************************C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C ************************************************PROGRAM ELASTICITYcharacter*32 dat,cchDIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)WRITE(*,*)'PLEASE ENTER INPUT 'READ(*,'(A)')DATOPEN(4,'OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',npWRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',neWRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',NrCALL INPUT (JR,COOR,AE,MEL)CALL CBAND (MA,JR,MEL)DO I=1,NHSK(I)=0、0enddoCALL SK0(SK,MEL,COOR,JR,MA,AE)do I=1,NR(I)=0、0enddopause 'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,izIF(NCP、GT、0)CALL CONCR(NCP,R,JR)IF(NBE、GT、0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)IF(iz、GT、0)thendo jj=1,izREAD (4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendifCALL DECOP (SK,MA)CALL FOBA (SK,MA,R)CALL OUTDISP(NP,R,JR)CALL STRESS (COOR,MEL,JR,AE,R,STRE)WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'STOPc RETURNENDC *********************************************SUBROUTINE INPUT (JR,COOR,AE,MEL)DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I)) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8、3))) RETURNENDC **********************************************SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,4IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,2*4NNI=NN(I)IF(NNI、EQ、0) GO TO 80IF(NNI、LT、L) L=NNI80 CONTINUEDO 85 M=1,2*4JP=NN(M)IF(JP、EQ、0) GO TO 85JPL=JP-L+1IF(JPL、GT、MA(JP)) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I)、GT、MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',NWRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH 500 FORMAT (/5X,'FREEDOM N='*,I5,3X,'SEMI-BANDWI、MX=',I5,3X,* 'STORAGE NH=',I7)RETURNENDC **********************************************SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN4/ NEE,NMEMON /GAUSS/ RSTG(3),H(3)H(1)=0、5555555555555560H(2)=0、8888888888888890H(3)=H(1)RSTG(1)=-0、7745966692414830RSTG(2)=0、00RSTG(3)=-RSTG(1)DO 10 IE=1,NENEE=IENME=MEL(5,IE)DO 75 K=1,4IEK=MEL(K,IE)iven(k)=IEKDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUECALL STIF(XYZ,AE,iven)DO 60 I=1,8DO 60 J=1,8II=NN(I)JJ=NN(J)IF ((JJ、EQ、0)、OR、(II、LT、JJ)) GO TO 60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60 CONTINUE70 CONTINUEwrite(7,1111) ((ske(i,j),j=1,8),i=1,8)1111 format(2x,8f12、2)10 CONTINUERETURNENDC *********************************************SUBROUTINE STIF(XYZ,AE,iven)DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),* RJAC(2,2)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN4/ NEE,NMEMON /GAUSS/ RSTG(3),H(3)DO 40 I=1,8RF(I)=0、00DO 30 J=1,8SKE(I,J)=0、0030 CONTINUE40 CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(3,NME)D1=E*(1、00-U)/((1、00+U)*(1、00-2、00*U))D2=E*U/((1、00+U)*(1、00-2、00*U))D3=E*0、50/(1、00+U)DO 120 I=1,4II=2*(I-1)I1=II+1I2=II+2DO 115 J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO 99 IS=1,3S=RSTG(IS)SH=H(IS)DO 98 IR=1,3R=RSTG(IR)RH=H(IR)CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE)DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SHDYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUERETURNENDC *********************************************SUBROUTINE CONCR(NCP,R,JR)DIMENSION R(*),JR(2,*),XYZ(2)DO 100 I=1,NCPREAD (4,*) IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO 95 J=1,2L=JR(J,IP)IF(L、EQ、0) GO TO 95R(L)=R(L)+XYZ(J)95 CONTINUE100 CONTINUERETURNENDC **********************************************SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE)DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),& AE(4,*),XYZ(2,4),iven(4)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)MON /GAUSS/ RSTG(3),H(3)H(1)=1、0H(2)=1、0RSTG(1)=-0、57735RSTG(2)=-RSTG(1)DO 10 IE=1,NBEDO I=1,8RF(I)=0、00ENDDOc READ(4,*)NEENEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO 75 K=1,4IEK=MEL(K,NEE)iven(k)=iekDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUEDO 99 IS=1,2S=RSTG(IS)SH=H(IS)DO 98 IR=1,2RR=RSTG(IR)RH=H(IR)CALL FUN8 (XYZ,RR,S,DET)DO 30 I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30 CONTINUE98 CONTINUE99 CONTINUECALL ASLOAD (R)10 CONTINUERETURNENDC *********************************************SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG)DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*)* ,XYZ(2,4),iew(*),PR(2)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN4/ NEE,NMEMON /GAUSS/ RSTG(3),H(3)H(1)=1、0H(2)=1、0RSTG(1)=-0、57735RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0、1if(ir、eq、1)nwf=1if(ir、eq、2)nnf=1DO 510 IE=1,NSEDO I=1,8RF(I)=0、00ENDDOnee=iew(ie)DO 575 K=1,4IEK=MEL(K,NEE)DO 595 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUEIF(NWF、EQ、1) thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0、1CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF、EQ、1) thenq=WG(2)NSU=WG(4)+0、1do j=1,2PR(J)=qenddoCALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)endifCALL ASLOAD (R)510 CONTINUERETURNENDC *********************************************SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FV AL(4),NODES(2),FACT(4)MON /CMN1/ NP,NE,NM,NRMON /CMN2/ N,MX,NHMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN4/ NEE,NMEMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)MON /GAUSS/ RSTG(3),H(3)DATA KCRD/1,1,2,2/DATA KFACE/1, 4,* 2, 3,* 1, 2,* 4, 3/DATA FVAL/-1、00,1、00,-1、00,1、00/FACT(1)=1、0FACT(2)=-1、0FACT(3)=-1、0FACT(4)=1、0FACTNUS=FACT(NSU)DO I=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF (NSI、EQ、1) THENDO I=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0、00IF (Z、GT、0、00) PR(I)=Z*GAMAENDDOENDIFML=KCRD(NSU)IF(ML、EQ、1)MM=2IF(ML、EQ、2)MM=1RST(ML)=FVAL(NSU)DO 70 LX=1,2RST(MM)=RSTG(LX)CALL FUN8 (XYZ,RST(1),RST(2),DET)PXYZ=0、00DO 25 I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)25 CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)30 DO 60 I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A260 CONTINUE70 CONTINUERETURNENDC *********************************************SUBROUTINE ASLOAD (R)DIMENSION R(*)MON /CMN1/ NP,NE,NM,NRMON /CMN3/ RF(8),SKE(8,8),NN(8)DO 20 I=1,8L=NN(I)IF (L、EQ、0) GO TO 20R(L)=R(L)+RF(I)20 CONTINUERETURNENDC ***********************************************SUBROUTINE DECOP (SK,MA)DIMENSION SK(*),MA(*)MON /CMN2/ N,MX,NHDO 50 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF (L1、GT、K) GO TO 30DO 20 J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF (L、GT、M) M=LMP=J-1IF (M、GT、MP) GO TO 20DO 10 LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10 CONTINUE20 CONTINUE30 IF (L、GT、K) GO TO 50DO 40 LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP)=SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40 CONTINUE50 CONTINUERETURNENDC *************************************************************SUBROUTINE FOBA (SK,MA,R)DIMENSION SK(*),MA(*),R(*)MON /CMN2/ N,MX,NHDO 10 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF (L、GT、K) GO TO 10DO 5 LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5 CONTINUE10 CONTINUEDO 20 I=1,NII=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUEDO 30 J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF (L、GT、K) GO TO 30DO 25 J=L,KIJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUERETURNENDC *****************************************************************SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),& B(3,8),R(*),iven(4)MON /CMN1/ NP,NE,NM,NRMON /CMN3/ RF(8),SKE(8,8),NN(8)MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DO 106 IE=1,NENME=MEL(5,IE)DO 300 K=1,4IEK=MEL(K,IE)DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK)DO 320 M=1,2JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1、00-U)/((1、00+U)*(1、00-2、00*U))D2=E*U/((1、00+U)*(1、00-2、00*U))D3=0、50*E/(1、00+U)SS=0、0RR=0、0CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE)DO 30 I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0、B(3,J1)=CIB(1,J2)=0、B(2,J2)=CIB(3,J2)=BI30 CONTINUEDO 55 II=1,3SIG(II)=0、0055 CONTINUEDO 70 K=1,8NA=NN(K)IF (NA、EQ、0) GO TO 70DO 60 L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60 CONTINUE70 CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106 CONTINUECALL OUTSTRE(NE,STRE)RETURNENDC *********************************************SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE)DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)CALL FUN8 (XYZ,R,S,DET)IF (DET、LT、1、0E-5)THENWRITE(7,600) NEE,R,S,detWRITE(7,*) (iven(m),m=1,4)STOPENDIFREC=1、00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO 30 K=1,4DO 20 I=1,2DNX(I,K)=0、DO 25 M=1,2DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)25 CONTINUE20 CONTINUE30 CONTINUE600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '* 'JACOBIAN DETERMINANT FOR '* 'ELEMENT'/'ELE、=',I5,' R=',F10、5,6X,'S=',F10、5,* 'det=',f12、5)RETURNENDC *********************************************SUBROUTINE FUN8 (XYZ,R,S,DET)DIMENSION XYZ(2,*),XI(4),ETA(4)MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DATA XI/-1、0,1、0,1、0,-1、0/DATA ETA/-1、0,-1、0,1、0,1、0/DO 10 I=1,4G1=(1、0+XI(I)*R)G2=(1、0+ETA(I)*S)FUN(I)=0、25*G1*G2PN(1,I)=0、25*XI(I)*G2PN(2,I)=0、25*ETA(I)*G110 CONTINUEDO 80 I=1,2DO 75 J=1,2DET=0、00DO 70 K=1,4DET=DET+PN(I,K)*XYZ(J,K)70 CONTINUEXJAC(I,J)=DET75 CONTINUE80 CONTINUEDET=XJAC(1,1)*XJAC(2,2)* -XJAC(2,1)*XJAC(1,2)RETURNENDC **********************************************SUBROUTINE OUTDISP(NP,R,JR)DIMENSION R(*),JR(2,*),U(2)WRITE(*,650)WRITE(7,650)DO I=1,NPDO M=1,2L=JR(M,I)IF(L、EQ、0)U(M)=0、0IF(L、GT、0)U(M)=R(L)ENDDOWRITE(*,'(5X,I5,10X,2E14、3)') I,UWRITE(7,'(5X,I5,10X,2E14、3)') I,UENDDO650 FORMAT(/25X,'NODAL DISPLACEMENTS'/8X, * 'NODE',13X,'X-P、',8X,'Y-P、')RETURNENDC **********************************************SUBROUTINE OUTSTRE(NE,STRE)DIMENSION STRE(3,*),ST(6)WRITE(*,700)WRITE(7,700)DO IE=1,NESX=STRE(1,IE)SY=STRE(2,IE)SXY=STRE(3,IE)ST(1)=SXST(2)=SYST(3)=SXYH1=SX+SYH2=SQRT((SX-SY)*(SX-SY)+4、0*SXY*SXY)ST(4)=(H1+H2)/2、0ST(5)=(H1-H2)/2、0IF(ABS(SXY)、LT、1、0E-4)THENIF (SX、GT、SY) ST(6)=0、0IF (SX、LE、SY) ST(6)=90、0ELSEST(6)=ATAN((ST(4)-SX)/SXY)*57、29578ENDIFWRITE(*,'(6X,I4,3X,6F11、3)') IE,STWRITE(7,'(6X,I4,3X,6F11、3)') IE,STENDDO700 FORMAT(/30X,'ELEMENT STRESSES'/5X, *'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',4X,'ANGLE'/)RETURNENDC *********************************************。
第九章 平面广义四节点等参元 GQ4
18.24 23.00 23.24 23.06 23.43 23.78 23.02 23.04 23.69 23.67 23.82 23.9
0.2113 0.2209 0.2287 0.2221 0.23337 0.2226 0.222 0.2661 0.2261 0.2393 0.2377 0.2360
式中 N e 为单元 e 的插值函数矩阵 值函数对应的广义形函数形式为
De 为单元 e 的广义自由度向量 一阶广义节点插 y , (i = 1,2,3,4) 0 x2 0 y2 0
1 0 x 0 y N ei = ϕ ei 0 1 0 x 0 1 0 x 0 y N =ϕ 0 1 0 x 0
9.3.3
剪切自锁考查
MacNeal 细长梁问题
弹性模量为 E = 106 泊松比为
计算简图见图 9-3
材料参数选为
ν = 0.3 纯
弯和端部受剪两种工况 计算网格及工况如图 9-3 中(a) 格计算结果列于表 9-3 中
(b)和(c)
不同工况下各种网
0.5 50 0.2 1 工况 1 0.5 50 工况 2
( ) (E )(B )
i T e g
(9-10)
j e g
式中 n g 为单元内高斯点个数 取值
t 表示材料的厚度
下标 g 表示该表达式在高斯点处的
W g 为高斯点积分权系数 具体的数值实现步骤如下
(1) 按照所选用的高斯积分阶次 (2) 按单元节点循环 i. 形成该单元的所有高斯点局部坐标 (ξ i ,η i )
vC
网格 a
ε A (max) ε B (min)
vC
网格 b
ε A (max) ε B (min)
最新平面四边形4结点等参有限单元法
有限元程序设计平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术，自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较，并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框图2-１程序流程图图2-２子程序框图其中，各子程序的主要功能为：INPUT――输入原始数据HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）SKO――形成整体刚度矩阵[K]CONCR――计算集中力引起的等效结点荷载{R}eBODYR――计算自重体力引起的等效结点荷载{R}eFACER――计算分布面力引起的等效结点荷载{R}eDECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数TiiyNxN⎥⎦⎤⎢⎣⎡∂∂∂∂，=i1,2,3,4。
FUN8――计算形函数及雅可比矩阵[J]SFUN ――应力磨平-单元下的‘K’=NCN‘SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。
四边形四节点等参元matlab程序
悬臂钢梁，尺寸如图一所示；v=0.3。
h=1，E=2.1e11.图一悬臂钢梁图二单元划分与结点编号附录：%---------------四边形四节点等参元matlab计算程序----------------------------% 2013年% 13级建筑与土木工程Bruce% B15-405% 本程序只输出了节点位移、单元中心点的应力%******************************************************************* % 变量说明% E v h% 弹性模量泊松比厚度% NPOIN NELEM NVFIX NNODE NFPOIN% 总结点数, 单元数, 约束结点个数, 单元节点数,受力结点数% COORD LNODS% 结构节点整体坐标数组, 单元定义数组,% FPOIN FORCE FIXED% 结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%******************************clear allformat short eFP1=fopen(' sjd.txt','rt'); %打开数据文件%%读入控制数据E=fscanf(FP1,'%f',1); %弹性模量v=fscanf(FP1,'%f',1); % 泊松比h=fscanf(FP1,'%f',1); %厚度NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); % 总结点数NNODE=fscanf(FP1,'%d',1); %单元节点数NFPOIN=fscanf(FP1,'%d',1); %受力结点数NVFIX=fscanf(FP1,'%d',1); %约束结点个数LNODS=fscanf(FP1,'%f',[NNODE,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号x,y坐标(整体坐标下) FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';% 节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0% 刚度矩阵的生成%计算刚度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零%调用子程序生成单元刚度矩阵for m=1:NELEM %m为单元号Ke=K(E,v,h,...COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2)); %调用单元刚度矩阵a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理for j=1:4for k=1:4HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+...Ke(j*2-1:j*2,k*2-1:k*2);endendend %—————————————————————————————————%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOINb1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)为作用点FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力end %—————————————————————————————————%将约束信息加入总刚,总荷载for i=1:NVFIXif FIXED(i,2)==1c1=2*FIXED(i,1)-1;HK(c1,:)=0; %将一约束序号处的总刚列向量清0HK(:,c1)=0; %将一约束序号处的总刚行向量清0HK(c1,c1)=1; %将行列交叉处的元素置为1FORCE(c1)=0;endif FIXED(i,3)==1c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0;endendDISP=HK\FORCE %计算节点位移向量%———————————求解单元应力————————————————stress=zeros(3,NELEM);for m=1:NELEMu(1:8)=0;d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号for i=1:NNODEu(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);%从总位移向量中取出当前单元的节点位移endD=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%形成应变矩阵BMBM=zeros(3,8);for i=1:NNODEJ=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),...COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...COORD(LNODS(m,4),1),COORD(LNODS(m,4),2),0,0);[N_s,N_t]=DHS(0,0);B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);endstressm=D*BM*u';stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x2,y2,x3,y3,x4,y4)%=========单元刚度矩阵=============== % E 弹性模量% v 泊松比% h 厚度% x1,y1,x2,y2,x3,y3,x4,y4 为4个角结点的坐标%矩阵尺寸为8 x 8Ke=zeros(8,8);D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵%高斯积分采用3 x 3 个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W=[W1 W2 W3];r=15^(1/2)/5;x=[-r 0 r];%积分点for i=1:3for j=1:3B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j));Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;endendfunction B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,s,t)%调用导函数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t);%求应变矩阵BB=zeros(3,8);for i=1:4B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);B(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i];endB=B/det(J);function J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t) %-------Jacobi-----------%单元坐标x=[x1 x2 x3 x4];y=[y1 y2 y3 y4];%%调用形函数对局部坐标的导数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;for i=1:4x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);endJ=[x_s y_s;x_t y_t];function N=shape(s,t)%ξ,ηN(1)=(1-s)*(1-t)/4;N(2)=(1+s)*(1-t)/4;N(3)=(1+s)*(1+t)/4;N(4)=(1-s)*(1+t)/4;function [N_s,N_t]=DHS(s,t)%形函数求导%ξ,ηN_s(1)=-1/4*(1-t);N_s(2)=1/4*(1-t);N_s(3)=1/4*(1+t);N_s(4)=-1/4*(1+t);N_t(1)=-1/4*(1-s);N_t(2)=-1/4*(1+s);N_t(3)=1/4*(1+s);N_t(4)=1/4*(1-s);sjd.txt 文件数据2.1E11 0.3 1 5 12 4 1 21 2 8 72 3 9 83 4 10 94 5 11 105 6 12 110.0 0.01.0 0.02.0 0.03.0 0.04.0 0.05.0 0.00.0 1.01.0 1.02.0 1.03.0 1.04.0 1.05.0 1.06 0 -100001 1 17 1 1。
有限元考试复习资料(华东交通大学)
有限元考试复习资料（含习题答案）1试说明用有限元法解题的主要步骤。
（1）离散化：将一个受外力作用的连续弹性体离散成一定数量的有限小的单元集合体，单元之间只在结点上互相联系，即只有结点才能传递力。
（2）单元分析：根据弹性力学的基本方程和变分原理建立单元结点力和结点位移之间的关系。
（3）整体分析：根据结点力的平衡条件建立有限元方程，引入边界条件，解线性方程组以及计算单元应力。
（4）求解方程，得出结点位移（5）结果分析，计算单元的应变和应力。
2.单元分析中,假设的位移模式应满足哪些条件,为什么?要使有限元解收敛于真解，关键在于位移模式的选择，选择位移模式需满足准则：（1）完备性准则：（2）连续性要求。
P210面简单地说，当选取的单元既完备又协调时，有限元解是收敛的，即当单元尺寸趋于0时，有限元解趋于真正解，称此单元为协调单元；当单元选取的位移模式满足完备性准则但不完全满足单元之间的位移及其导数连续条件时，称为非协调单元。
3.什么样的问题可以用轴对称单元求解？在工程问题中经常会遇到一些实际结构，它们的几何形状、约束条件和外载荷均对称某一固定轴，我们把该固定轴称为对称轴。
则在载荷作用下产生的应力、应变和位移也都对称此轴。
这种问题就称为轴对称问题。
可以用轴对称单元求解。
4.什么是比例阻尼?它有什么特点?其本质反映了阻尼与什么有关?答:比例阻尼:由于多自由度体系主振型关于质量矩阵与刚度矩阵具有正交性关系，若主振型关于阻尼矩阵亦具有正交性,这样可对多自由度地震响应方程进行解耦分析。
比例阻尼的特点为具有正交性。
其本质上反应了阻尼与结构物理特性的关系。
5.何谓等参单元?等参单元具有哪些优越性?①等参数单元(简称等参元)就是对坐标变换和单元内的参变量函数(通常是位移函数)采用相同数目的节点参数和相同的插值函数进行变换而设计出的一种单元。
①优点:可以很方便地用来离散具有复杂形体的结构。
由于等参变换的采用使等参单元特性矩阵的计算仍在单元的规则域内进行，因此不管各个积分形式的矩阵表示的被积函数如何复杂，仍然可以方便地采用标准化的数值积分方法计算。
matlab四节点四边形等参元的刚度矩阵计算程序
-1.2150 0.1440 2.1968 -0.8933 0.5212 -0.0151 -1.5030 0.7645
0.0616 -0.0419 -0.8933 1.8399 0.0673 -0.5250 0.7645 -1.2729
x2=2;y2=0;
x3=2.25;y3=1.5;
x4=1.25;y4=1;
材料参数：
E=30e12;
NU=0.3;
h=1;
ID=1;
Ai=1
Aj=1;
3、程序（由附件给出）
4、计算结果及讨论（需有图表来说明）
计算结果如下：
k =
1.0e+012 *
1.4619 0.2211 -1.2150 0.0616 -0.3716 -0.4236 0.1248 0.1409
-0.3716 -0.4236 0.5212 0.0673 1.1645 0.2763 -1.3141 0.0800
-0.4236 -0.2582 -0.0151 -0.5250 0.2763 0.9061 0.1624 -0.1229
0.1248 0.0585 -1.5030 0.7645 -1.3141 0.1624 2.6923 -0.9854
%输入平面问题性质参数ID（1为平面应力，2为平面应变）
%输出单元刚度矩阵
%-------------------------------------------------------------------
symsst;
a=[-(1-t)*x1+(1-t)*x2+(1+t)*x3-(1+t)*x4]/4;
平面四边形四节点等参单元Fortran源程序复习过程
平面四边形四节点等参单元F o r t r a n源程序C ************************************************C * FINITE ELEMENT PROGRAM *C * FOR Two DIMENSIONAL ELASticity PROBLEM *C * WITH 4 NODE *C ************************************************PROGRAM ELASTICITYcharacter*32 dat,cchDIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)WRITE(*,*)'PLEASE ENTER INPUT FILE NAME'READ(*,'(A)')DATOPEN(4,FILE=dat,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NRWRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',npWRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL)CALL CBAND (MA,JR,MEL)DO I=1,NHSK(I)=0.0enddoCALL SK0(SK,MEL,COOR,JR,MA,AE)do I=1,NR(I)=0.0enddopause 'aaa'stopREAD(4,*)NCP,NBE,izWRITE(*,'(5i8)')NCP,NBE,izWRITE(7,'(5i8)')NCP,NBE,izIF(NCP.GT.0)CALL CONCR(NCP,R,JR)IF(NBE.GT.0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)IF(iz.GT.0)thendo jj=1,izREAD (4,*)Js,nse,(WG(I),I=1,4)read(4,*)(iew(m),m=1,nse)CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)enddoendifCALL DECOP (SK,MA)CALL FOBA (SK,MA,R)CALL OUTDISP(NP,R,JR)CALL STRESS (COOR,MEL,JR,AE,R,STRE)WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'STOPc RETURNENDC *********************************************SUBROUTINE INPUT (JR,COOR,AE,MEL)DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 70 I=1,NPREAD(4,*) IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70 CONTINUEDO 11 J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,4)MEL(5,NEE)=NME11 CONTINUEDO 10 I=1,NPDO 10 J=1,210 JR(J,I)=1DO 20 I=1,NRREAD(4,*) IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20 CONTINUEN=0DO 30 I=1,NPDO 30 J=1,2IF (JR(J,I)) 30,30,2525 N=N+1JR(J,I)=N30 CONTINUEDO 55 J=1,NMREAD (4,*)JJ,(AE(I,JJ),I=1,4)WRITE(*,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8.3))) RETURNENDC **********************************************SUBROUTINE CBAND (MA,JR,MEL)DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHDO 65 I=1,N65 MA(I)=0DO 90 IE=1,NEDO 75 K=1,4IEK=MEL(K,IE)DO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUEL=NDO 80 I=1,2*4NNI=NN(I)IF(NNI.EQ.0) GO TO 80IF(NNI.LT.L) L=NNI80 CONTINUEDO 85 M=1,2*4JP=NN(M)IF(JP.EQ.0) GO TO 85JPL=JP-L+1IF(JPL.GT.MA(JP)) MA(JP)=JPL85 CONTINUE90 CONTINUEMX=0MA(1)=1DO 10 I=2,NIF(MA(I).GT.MX) MX=MA(I)MA(I)=MA(I)+MA(I-1)10 CONTINUENH=MA(N)WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH500 FORMAT (/5X,'FREEDOM N='*,I5,3X,'SEMI-BANDWI. MX=',I5,3X,* 'STORAGE NH=',I7)RETURNENDC **********************************************SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*), * AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=0.5555555555555560H(2)=0.8888888888888890H(3)=H(1)RSTG(1)=-0.7745966692414830RSTG(2)=0.00RSTG(3)=-RSTG(1)DO 10 IE=1,NENEE=IENME=MEL(5,IE)DO 75 K=1,4IEK=MEL(K,IE)iven(k)=IEKDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUECALL STIF(XYZ,AE,iven)DO 60 I=1,8DO 60 J=1,8II=NN(I)JJ=NN(J)IF ((JJ.EQ.0).OR.(II.LT.JJ)) GO TO 60JN=MA(II)-(II-JJ)SK(JN)=SK(JN)+SKE(I,J)60 CONTINUE70 CONTINUEwrite(7,1111) ((ske(i,j),j=1,8),i=1,8)1111 format(2x,8f12.2)10 CONTINUERETURNENDC *********************************************SUBROUTINE STIF(XYZ,AE,iven)DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),* RJAC(2,2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)DO 40 I=1,8RF(I)=0.00DO 30 J=1,8SKE(I,J)=0.0030 CONTINUE40 CONTINUEE=AE(1,NME)U=AE(2,NME)GAMA=AE(3,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=E*0.50/(1.00+U)DO 120 I=1,4II=2*(I-1)I1=II+1I2=II+2DO 115 J=1,4JJ=2*(J-1)J1=JJ+1J2=JJ+2DXX=0DXY=0DYX=0DYY=0DO 99 IS=1,3S=RSTG(IS)SH=H(IS)DO 98 IR=1,3R=RSTG(IR)RH=H(IR)CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DNIX=DNX(1,I)DNIY=DNX(2,I)DNJX=DNX(1,J)DNJY=DNX(2,J)DXX=DXX+DNIX*DNJX*DET*RH*SHDXY=DXY+DNIX*DNJY*DET*RH*SHDYX=DYX+DNIY*DNJX*DET*RH*SHDYY=DYY+DNIY*DNJY*DET*RH*SH98 CONTINUE99 CONTINUESKE(I1,J1)=DXX*D1+DYY*D3SKE(I2,J2)=DYY*D1+DXX*D3SKE(I1,J2)=DXY*D2+DYX*D3SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUERETURNENDC ********************************************* SUBROUTINE CONCR(NCP,R,JR)DIMENSION R(*),JR(2,*),XYZ(2)DO 100 I=1,NCPREAD (4,*) IP,PX,PYXYZ(1)=PXXYZ(2)=PYDO 95 J=1,2L=JR(J,IP)IF(L.EQ.0) GO TO 95R(L)=R(L)+XYZ(J)95 CONTINUE100 CONTINUERETURNENDC ********************************************** SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),& AE(4,*),XYZ(2,4),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)DO 10 IE=1,NBEDO I=1,8RF(I)=0.00ENDDOc READ(4,*)NEENEE=ieNME=MEL(5,NEE)GAMA=AE(3,NME)DO 75 K=1,4IEK=MEL(K,NEE)iven(k)=iekDO 95 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)95 XYZ(M,K)=COOR(M,IEK)75 CONTINUEDO 99 IS=1,2S=RSTG(IS)SH=H(IS)DO 98 IR=1,2RR=RSTG(IR)RH=H(IR)CALL FUN8 (XYZ,RR,S,DET)DO 30 I=1,4J=2*IRF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA30 CONTINUE98 CONTINUE99 CONTINUECALL ASLOAD (R)10 CONTINUERETURNENDC *********************************************SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG) DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*) * ,XYZ(2,4),iew(*),PR(2)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /GAUSS/ RSTG(3),H(3)H(1)=1.0H(2)=1.0RSTG(1)=-0.5773502691896260RSTG(2)=-RSTG(1)nwf=0nnf=0ir=wg(1)+0.1if(ir.eq.1)nwf=1if(ir.eq.2)nnf=1DO 510 IE=1,NSEDO I=1,8RF(I)=0.00ENDDOnee=iew(ie)DO 575 K=1,4IEK=MEL(K,NEE)DO 595 M=1,2JJ=2*(K-1)+MNN(JJ)=JR(M,IEK)595 XYZ(M,K)=COOR(M,IEK)575 CONTINUEIF(NWF.EQ.1) thenGAMA=WG(2)Z0=WG(3)NSU=WG(4)+0.1CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)endifIF(NNF.EQ.1) thenq=WG(2)NSU=WG(4)+0.1do j=1,2PR(J)=qenddoCALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)endifCALL ASLOAD (R)510 CONTINUERETURNENDC *********************************************SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4), & FVAL(4),NODES(2),FACT(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN2/ N,MX,NHCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN4/ NEE,NMECOMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)COMMON /GAUSS/ RSTG(3),H(3)DATA KCRD/1,1,2,2/DATA KFACE/1, 4,* 2, 3,* 1, 2,* 4, 3/DATA FVAL/-1.00,1.00,-1.00,1.00/FACT(1)=1.0FACT(2)=-1.0FACT(3)=-1.0FACT(4)=1.0FACTNUS=FACT(NSU)DO I=1,2J=KFACE(I,NSU)NODES(I)=JENDDOIF (NSI.EQ.1) THENDO I=1,2J=NODES(I)Z=Z0-XYZ(2,J)PR(I)=0.00IF (Z.GT.0.00) PR(I)=Z*GAMAENDDOENDIFML=KCRD(NSU)IF(ML.EQ.1)MM=2IF(ML.EQ.2)MM=1RST(ML)=FVAL(NSU)DO 70 LX=1,2RST(MM)=RSTG(LX)CALL FUN8 (XYZ,RST(1),RST(2),DET) PXYZ=0.00DO 25 I=1,2J=NODES(I)PXYZ=PXYZ+FUN(J)*PR(I)25 CONTINUEA1=XJAC(MM,2)A2=-XJAC(MM,1)30 DO 60 I=1,2J=NODES(I)K2=2*JK1=K2-1Q=PXYZ*FUN(J)*H(LX)*FACTNUSRF(K1)=RF(K1)+Q*A1RF(K2)=RF(K2)+Q*A260 CONTINUE70 CONTINUEENDC ********************************************* SUBROUTINE ASLOAD (R)DIMENSION R(*)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)DO 20 I=1,8L=NN(I)IF (L.EQ.0) GO TO 20R(L)=R(L)+RF(I)20 CONTINUERETURNENDC *********************************************** SUBROUTINE DECOP (SK,MA)DIMENSION SK(*),MA(*)COMMON /CMN2/ N,MX,NHDO 50 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1L1=L+1IF (L1.GT.K) GO TO 30DO 20 J=L1,KIJ=MA(I)-I+JM=J-MA(J)+MA(J-1)+1IF (L.GT.M) M=LMP=J-1IF (M.GT.MP) GO TO 20DO 10 LP=M,MPIP=MA(I)-I+LPJP=MA(J)-J+LPSK(IJ)=SK(IJ)-SK(IP)*SK(JP)10 CONTINUE20 CONTINUE30 IF (L.GT.K) GO TO 50DO 40 LP=L,KIP=MA(I)-I+LPLPP=MA(LP)SK(IP)=SK(IP)/SK(LPP)II=MA(I)SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)40 CONTINUE50 CONTINUEENDC *************************************************************SUBROUTINE FOBA (SK,MA,R)DIMENSION SK(*),MA(*),R(*)COMMON /CMN2/ N,MX,NHDO 10 I=2,NL=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 10DO 5 LP=L,KIP=MA(I)-I+LPR(I)=R(I)-SK(IP)*R(LP)5 CONTINUE10 CONTINUEDO 20 I=1,NII=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUEDO 30 J1=2,NI=2+N-J1L=I-MA(I)+MA(I-1)+1K=I-1IF (L.GT.K) GO TO 30DO 25 J=L,KIJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUERETURNENDC ***************************************************************** SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),& B(3,8),R(*),iven(4)COMMON /CMN1/ NP,NE,NM,NRCOMMON /CMN3/ RF(8),SKE(8,8),NN(8)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DO 106 IE=1,NENME=MEL(5,IE)DO 300 K=1,4IEK=MEL(K,IE)DO 310 M=1,2310 XYZ(M,K)=COOR(M,IEK)DO 320 M=1,2JRR=2*(K-1)+M320 NN(JRR)=JR(M,IEK)300 CONTINUEE=AE(1,NME)U=AE(2,NME)D1=E*(1.00-U)/((1.00+U)*(1.00-2.00*U))D2=E*U/((1.00+U)*(1.00-2.00*U))D3=0.50*E/(1.00+U)SS=0.0RR=0.0CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE) DO 30 I=1,4II=2*(I-1)J1=II+1J2=II+2BI=DNX(1,I)CI=DNX(2,I)B(1,J1)=BIB(2,J1)=0.B(3,J1)=CIB(1,J2)=0.B(2,J2)=CIB(3,J2)=BI30 CONTINUEDO 55 II=1,3SIG(II)=0.0055 CONTINUEDO 70 K=1,8NA=NN(K)IF (NA.EQ.0) GO TO 70DO 60 L=1,3SIG(L)=SIG(L)+B(L,K)*R(NA)60 CONTINUE70 CONTINUESX=D1*SIG(1)+D2*SIG(2)SY=D2*SIG(1)+D1*SIG(2)SXY=D3*SIG(3)STRE(1,IE)=SXSTRE(2,IE)=SYSTRE(3,IE)=SXY106 CONTINUECALL OUTSTRE(NE,STRE)RETURNENDC *********************************************SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE) DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)CALL FUN8 (XYZ,R,S,DET)IF (DET.LT.1.0E-5)THENWRITE(7,600) NEE,R,S,detWRITE(7,*) (iven(m),m=1,4)STOPENDIFREC=1.00/DETRJAC(1,1)=REC*XJAC(2,2)RJAC(2,2)=REC*XJAC(1,1)RJAC(2,1)=-REC*XJAC(2,1)RJAC(1,2)=-REC*XJAC(1,2)DO 30 K=1,4DO 20 I=1,2DNX(I,K)=0.DO 25 M=1,2DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)25 CONTINUE20 CONTINUE30 CONTINUE600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '* 'JACOBIAN DETERMINANT FOR '* 'ELEMENT'/'ELE.=',I5,' R=',F10.5,6X,'S=',F10.5,* 'det=',f12.5)RETURNENDC *********************************************SUBROUTINE FUN8 (XYZ,R,S,DET)DIMENSION XYZ(2,*),XI(4),ETA(4)COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)DATA XI/-1.0,1.0,1.0,-1.0/DATA ETA/-1.0,-1.0,1.0,1.0/DO 10 I=1,4G1=(1.0+XI(I)*R)G2=(1.0+ETA(I)*S)FUN(I)=0.25*G1*G2PN(1,I)=0.25*XI(I)*G2PN(2,I)=0.25*ETA(I)*G110 CONTINUEDO 80 I=1,2DO 75 J=1,2DET=0.00DO 70 K=1,4DET=DET+PN(I,K)*XYZ(J,K)70 CONTINUEXJAC(I,J)=DET75 CONTINUE80 CONTINUEDET=XJAC(1,1)*XJAC(2,2)* -XJAC(2,1)*XJAC(1,2)RETURNENDC ********************************************** SUBROUTINE OUTDISP(NP,R,JR)DIMENSION R(*),JR(2,*),U(2)WRITE(*,650)WRITE(7,650)DO I=1,NPDO M=1,2L=JR(M,I)IF(L.EQ.0)U(M)=0.0IF(L.GT.0)U(M)=R(L)ENDDOWRITE(*,'(5X,I5,10X,2E14.3)') I,UWRITE(7,'(5X,I5,10X,2E14.3)') I,UENDDO650 FORMAT(/25X,'NODAL DISPLACEMENTS'/8X, * 'NODE',13X,'X-COMP.',8X,'Y-COMP.')RETURNENDC ********************************************** SUBROUTINE OUTSTRE(NE,STRE)DIMENSION STRE(3,*),ST(6)WRITE(*,700)WRITE(7,700)DO IE=1,NESX=STRE(1,IE)SY=STRE(2,IE)SXY=STRE(3,IE)ST(1)=SXST(2)=SYST(3)=SXYH1=SX+SYH2=SQRT((SX-SY)*(SX-SY)+4.0*SXY*SXY)ST(4)=(H1+H2)/2.0ST(5)=(H1-H2)/2.0IF(ABS(SXY).LT.1.0E-4)THENIF (SX.GT.SY) ST(6)=0.0IF (SX.LE.SY) ST(6)=90.0ELSEST(6)=ATAN((ST(4)-SX)/SXY)*57.29578ENDIFWRITE(*,'(6X,I4,3X,6F11.3)') IE,STWRITE(7,'(6X,I4,3X,6F11.3)') IE,STENDDO700 FORMAT(/30X,'ELEMENT STRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',4X,'ANGLE'/)RETURNENDC *********************************************。
平面有限元分析-等参单元
等参元变换的条件为 J ≠ 0 ，因此在有限元网格划分时，要 特别注意这一点。
等参单元等效节点力（4节点）
（1）集中力引起的单元节点载荷
单元内某点受到集中载荷P=[Px Py]T，移置到单元节 点上的等效节点力为：
j
同理 得
dη = ∂x dηi + ∂y dη j
∂η
∂η
∂x dξ
dA =
dξ × dη
=
∂ξ
∂x
dη
∂µ
∂y dξ ∂x
∂ξ
∂y
dη
=
∂ξ
∂x
∂η ∂µ
∂y
∂ξ
∂y
dξdη
∂η
等参单元刚度（4节点）
因为
∂x ∂y
J
=
∂ξ
∂x
∂ξ
∂y
∂η ∂η
雅可比行列式 Jacobi
曲边面积元dA:
dA = J dξdη
8 平面问题有限元分析 等参单元
8.1等参曹单国元华刚度（4节点） 8.2等参单元等效节点力（4节点） 8.3矩形单元（8节点） 8.4等参单元（8节点） 8.5高斯积分法
等参单元刚度（4节点）
4
4
=u ∑= Ni (ξ ,η )ui ,v ∑ Ni (ξ ,η )vi
=i 1 =i 1
4
4
=x ∑= Ni (ξ ,η ) xi , y ∑ Ni (ξ ,η ) yi
i =1
4
y = ∑ Ni (ξ ,η ) yi
i =1
平面四节点等参单元分析程序
\变分原理与有限元大作业平面四节点等参单元分析程序.姓名：潘清学号：SQ完成时间：2011-4-26：一、概述通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。
有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。
因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。
因此，必须寻找新的编程语言。
随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。
现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。
C语言最初是为操作系统、编译器以及文字处理等编程而发明的。
随着不断完善，它已应用到其它领域，包括工程应用软件的编程。
近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。
用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。
C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。
有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析
2.6 四结点四边形单元(The four-node quadrilateral element)前面介绍了四结点的矩形单元其位移函数:xy y x U 4321αααα+++=xy y x V8765αααα+++=为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。
但它对边界要求严格。
本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。
对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。
可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。
正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。
正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。
如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。
当然, 局部坐标上的A 点与整体坐标的A 点对应。
一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V引入边界条件，即可得位移函数：∑=ijmpi i U N Ui ijmpi V N V ∑==写成矩阵形式：{}{}[]{}ee p i p i ed N d N N N N V U f =⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=000 式中形函数: ()()()ηηξξηξi i i N ++=1141, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmpi x N x N x N x N x N x +++==∑p p m m j j i i i ijmpi y N y N y N y N y N y +++==∑ ()162-- 式中形函数N 与位移函数中的完全一致。
第5章 平面问题有限元分析-等参单元
2020/6/30
平面问题有限元分析-等参单元
2
5.1四节点矩形单元位移函数
如图所示的矩形单元，不失一般性，令矩形单元的长、
宽分别为2a、2b。矩形单元有4个节点，共8个自由度，即 共有8个节点位移，采用类似三角形单元的分析方法，同样 可以完成对矩形单元的力学特性分析。
y
4 3
2b
o
1
2a
2
o
x
2020/6/30
14
5.2四节点矩形单元应变与应力矩阵
由前面的讨论可以发现，四边形单元的位移模式比常
应变三角形单元所采用的线性位移模式增添了项（即相
当于xy项），把这种位移模式称为双线性模式。在这种模 式 下 ， 单 元 内 的 应 变 分 量 将 不 再 是 常 量 ， 这 一 点 可 以 从Be 的表达式中看出。另外四边形单元的位移模式中的1 ~ 7 与 三角形单元相同，它反映了刚体位移和常应变，而且在单
求出α1, α2, α3, α4；α 5, α 6 , α7 , α8
u1 1 1 1
u2 1 1 1
u3 1 1 1
1
u4 1
1 1 1 1
1 1
1 1 1 1
11 1 1
1 1 1 1
1 1 1 1 1 u1
23 4
1 4
1
1
1
1 1 1
1 1 1
1
1 1
uu23 u4
5 1 1 1 1 v1
K
e rs
4ab
Et 1
2
K1
K3
K2
K4
式中：
K1
b2rs
1
rs
3
1
平面四边形4结点等参有限单元法
有限元程序设计平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术，自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较，并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框图2-１程序流程图图2-２子程序框图其中，各子程序的主要功能为：INPUT――输入原始数据HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息CBAND――形成主元素序号指示矩阵MA（*）SKO――形成整体刚度矩阵[K]CONCR――计算集中力引起的等效结点荷载{R}eBODYR――计算自重体力引起的等效结点荷载{R}eFACER――计算分布面力引起的等效结点荷载{R}eDECOP――支配方程LU三角分解FOBA――LU分解直接解法中的回代过程OUTDISP――输出结点位移分量STRESS――计算单元应力分量OUTSTRE――输出单元应力分量STIF――计算单元刚度矩阵FDNX――计算形函数对整体坐标的导数TiiyNxN⎥⎦⎤⎢⎣⎡∂∂∂∂，=i1,2,3,4。
FUN8――计算形函数及雅可比矩阵[J]SFUN ――应力磨平-单元下的‘K’=NCN‘SCN――应力磨平-单元下的右端项系数‘CN‘SUMSKN――应力磨平-单元下的右端项集成到总体的‘P‘SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。
04-等参数单元
第四章 等参数单元为了方便应用和提高计算精度，目前许多实用程序采用了等参数单元，取得了较好的效果。
本章从平面问题的任意四边形单元入手，介绍等参数单元的一些基本概念，并根据工程实际应用需要，重点介绍空间六面体等参数单元分析。
§4-1 等参数单元的概念一、平面等参数单元1. 四节点四边形等参数单元在平面问题的有限单元法中，最简单和最常用的是三节点三角形单元，其次是四节点矩形单元。
由于三节点三角形单元采用的是线性位移模式，是实际位移分布的最简单逼近形式，求解精度受到限制。
而四节点矩形单元，由于它的位移函数是坐标的二次函数，单元内的应力不是常量而是线性变化的，所以能比简单三角形单元较好地反映实际应力变化情况。
但是矩形单元不能适应曲线边界和非正交的直线边界，也不便随意改变大小，应用范围受到很大限制。
如果采用任意四边形单元(如图4-1所示)，而仍采用矩形单元的位移模式，则基本上能够克服矩形单元的上述不足之处。
但是此时在不平行于坐标轴的边界上，由于y=kx(k ≠0)，位移函数为C Bx Axxy y x u ++=+++=24321αααα是坐标的二次函数，这样就不能由边界上二个节点的位移来唯一地确定该边界上的位移，故位移的连续性将得不到保证，其变形协调条件就得不到满足。
采用坐标变换可解决这一矛盾，现说明如下。
在图4-1所示的任意 四边形单元上，用等分四边的两族直线分割该四边形。
以两族直线的中心为原点(ξ=η=0)，并令四边上的ξ值、η值分别为±1,这样就得到一个新的坐标系，单元上的任一点都取一个新的坐标(ξ，η)。
这里的ξ，η是一种局部坐标，只适用于一个单元的范围内。
与此相反，原坐标x ，y 则是一种整体坐标，和以前一样地通用于所有单元的整体。
为确保局部坐标ξ，η和原坐标x,y 有一一对应关系，即存在确定的坐标变换关系，应使任意四边形不能大歪斜，它的任意 一条边的延伸线不能再分割单元(如图4-2)。
土木工程中的数值方法-6-有限单元法-等参单元
x y
J
x
y
Ni x
N
i
x
y Ni Ni
y
x N i
J
x N i
y y
Ni
Ni
x N i
J
1
N
i
y
4.单元刚度矩阵
代入坐标变换关系，得到
x1
y1
x y
N1 0
0 N1
N2 0
0 N2
N3 0
0 N3
N3 0
0 N3
N4 0
0 N4
e
B
e
y x
[B]为应变转换矩阵。
单元内应力：
x
y
DB
e
z
4.单元刚度矩阵
有虚功原理得到单元刚度矩阵
ke e BT DBt dΩ
面积元转化 x
dΩ x
y
dd J dd
y
单元刚度矩阵公式：
k e
1 1 B T 1 1
DBt
J
dd
4.单元刚度矩阵
0 N3
N4 0
0 N4
x2 y2 x3
y3
x4
y4
利用划线法构造局部坐标系里的形函数：
N1
1 4
(1
)(1 )
N2
1 4
(1 )(1)
N3
1 4
(1 )(1)
N4
1 4
(1 )(1)
2.坐标变换与平面四结点等参单元
用同样的形状函数来插值单元内任意一点(x, y)的位移。
u v
N1 0
0 N1
N2 0
N1
1 4
(1
)(1 )