平面四节点等参单元分析程序

合集下载

四边形四结点等参元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）用虚功原理求得单元刚度矩阵，最后用高斯积分法计算完成。

平面四节点等参单元有限元程序

exit(1);
}
//控制信息读入
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)

四边形四节点等参元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.何谓等参单元?等参单元具有哪些优越性?①等参数单元(简称等参元)就是对坐标变换和单元内的参变量函数(通常是位移函数)采用相同数目的节点参数和相同的插值函数进行变换而设计出的一种单元。

①优点:可以很方便地用来离散具有复杂形体的结构。

由于等参变换的采用使等参单元特性矩阵的计算仍在单元的规则域内进行，因此不管各个积分形式的矩阵表示的被积函数如何复杂，仍然可以方便地采用标准化的数值积分方法计算。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。