平面三角形3节点有限元程序

平面三角形３结点有限元程序1、程序名：FEMT3.FOR，FEMT3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。程序采用VisualFortran5.0编制而成，输入数据全部采用自由格式。3、程序流程及框图图１-1程序流程图输入数据信息形成K形成R计算，输出应力解Kδ＝R，输出δ启动停机图１-2程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的bi，ci等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移δ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块Krs。4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR——约束结点总数；NI——问题类型标识，0为平面应力问题，1为平面应变问题；NL——受荷载作用的结点的数目；NG——考虑自重作用为1，不计自重为0；ND——非零已知位移结点的数目；NC——要计算支座约束反力的结点数目。⑵材料信息，共NM条，每条依次输入EO，VO，W，tEO——弹性模量（kN/m2）；VO——泊松比；W——材料容重（kN/m3）；t——单元厚度（m）。这些信息都存放在数组AE（4，NM）中。⑶坐标信息，共NP条，每条依次输入IP，X，YIP——结点号；X，Y——分别为结点的x坐标和y坐标。坐标信息存放在数组X(2，NP)中。⑷单元信息，共NE条，每条依次输入JE，L，Io，Jo，MoJE——单元号；L——为该单元的材料类型号。Io，Jo，Mo——分别为该单元i，j，m的整体编码。单元信息存放在数组MEO(4,NE)中。⑸约束信息共NR条，每条依次输入一个数××IPxyIP——结点号；Ix，Iy——分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。⑹荷载信息，共NL条，每条依次输入yIP，Fx，FyIP——结点号；Fx，Fy——分别为该结点的x，y方向的荷载分量（kN）。结点号存放在数组NF（NL）中，结点荷载分量存放在数组FV（2，NL）中。⑺若ND＞0，输入非零已知位移信息，共ND条，每条依次输入IP，ux，uyIP——结点号；ux，uy——分别为该结点x，y方向已知位移分量（m），若其中某方向为自由，则其相应分量为0。结点号存放在数NDI（ND）中，已知位移分量存放在数组DV(2，ND)中。⑻支座反力信息，共NC条，每条依次输入IP，M1，M2，M3，M4IP——支座结点号；M1，M2，M3，M4——为与该支座结点相关的单元号，若不足4个，则用0补充。支座结点号存放在数组NCI(NC)中，相关单元号存放在数组NCE(4，NC)中。以上数据须按如上顺序存放在数据文件中。除此之外，程序中还用到其他一些主要变量和数组，说明如下：N——结构自由度总数；NH——按一维存贮的整体劲度矩阵的总容量；NX——最大半带宽；SK(10000)——维存贮的劲度矩阵；R(1000)——开始存放等效结点荷载，求解方程以后，用来存放结点位移；B(6)——存放单元应力σx，σy，τxy，σ1，σ2，α；MA(1000)——主元素序号指标矩阵；JR(2，500)——结点自由度序号矩阵；ME(3)——存放单元结点i，j，m的整体编码；NN(6)——单元结点自由度序号；BI(3)，CI(3)——单元劲度矩阵计算公式中的bi，bj，bm和ci，cj，cm；S——三角形单元的面积；H11，H12，H21，H22——单元劲度矩阵中子块Krs的4个元素。5、算例一个正方形弹性体，厚度为1m，四边受单位均布法向力作用，由于对称性，取其1/4进行计算，其有限元网格如图2-3所示，设27/1042mknE，167.0，不考虑自重。该问题的精确解应力为x＝12kN/m，y＝12kN/m，xy＝0。图１-3有限元网格（1）输入文件数据6415030052000.00.00.01.010.02.020.01.031.01.040.00.051.00.062.00.0113122124531325415631012014005106101-0.5-0.53-1.0-1.06-0.5-0.51100021230420005234064000（2）输出文件结果NODALDISPLACEMENTSNODEX-COMP.Y-COMP.10.00000E+00-0.10000E-0220.00000E+00-0.50000E-033-0.50000E-03-0.50000E-0340.00000E+000.00000E+005-0.50000E-030.00000E+006-0.10000E-020.00000E+00ELEMENTSTRESSESELEMENTX-STRESSY-STRESSXY-STRESSMAX-STRESSMIN-STRESSANGLE1-1.000-1.0000.000-1.000-1.00090.0002-1.000-1.0000.000-1.000-1.00090.0003-1.000-1.0000.000-1.000-1.00090.0004-1.000-1.0000.000-1.000-1.00090.000NODESTRESSESNODEX-STRESSY-STRESSXY-STRESSMAX-STRESSMIN-STRESSANGLE1-1.000-1.0000.000-1.000-1.00090.0002-1.000-1.0000.000-1.000-1.00090.0003-1.000-1.0000.000-1.000-1.00090.0004-1.000-1.0000.000-1.000-1.00090.0005-1.000-1.0000.000-1.000-1.00090.0006-1.000-1.0000.000-1.000-1.00090.000NODALREACTIONSNODEX-COMPY-COMP10.0000.00021.0000.00040.5000.50050.0001.00060.0000.0006、源程序CFINITEELEMENTPROGRAMFORTWODIMENSIONALCTRIANGLEELEMENTCDIMENSIONK(800000),COOR(2,3000),AE(4,11),*MEL(5,2000),MA(6000)CHARACTER*32datCOMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NCWRITE(*,300)300FORMAT(///''*':****'/'+PLEASEINPUTFILENAMEOFDATA')READ(*,*)dataOPEN(4,FILE=data,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NR,NI,NL,NG,ND,NCCWRITE(*,400)NP,NE,NM,NR,NI,NL,NG,ND,NCCWRITE(7,400)NP,NE,NM,NR,NI,NL,NG,ND,NCCALLINPUT(JR,COOR,MEL,AE)CALLCBAND(MA,JR,MEL)CALLSK0(SK,R,COOR,MEL,MA,JR,AE)CALLLOAD(COOR,MEL,R,JR,AE)IF(ND.GT.0)CALLTREAT(SK,MA,JR,R)CALLDECOMP(SK,MA)CALLFOBA(SK,MA,R)WRITE(*,650)WRITE(7,650)CALLOUTPUT(JR,R)WRITE(*,700)WRITE(7,700)CALLCES(COOR,MEL,JR,R,AE)IF(NC.GT.0)callERFAC(COOR,MEL,JR,R,AE)400FORMAT(/2X,'NP=',I3,2X,'NE=',I3,2X,'NM='*,I3,2X,'NR=',I3,2X,'NI='I3,2X,'NL=',I3,2X,*'NG=',I3,2X,'ND=',I3,2X,'NC=',I3)500FORMAT(1X,'TOTALDEGREESOFFREEDOMN=',*I4,1X,'TOTAL-STORAGE','NH=',I5,1X,*'MAX-SEMI-BANDWIDTHMX=',I3)550FORMAT(/20X,'TOTALSTORAGEIS*GREATERTHAN50000')600FORMAT(30X,'NODALFORCES'/8X,'NODE',*11X,'X-COMP.',14X,'Y-COMP.')650FORMAT(/30X,'NODALDISPLACEMENTS'/8X,*'NODE',13X,'X-COMP.',12X,'Y-COMP.')700FORMAT(/30X,'ELEMENTSTRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',6X,'ANGLE'/)STOPENDC*********************************************SUBROUTINEKRS(BR,BS,CR,CS)COMMON/CB/EO,VO,W,T,A,H11,H12,H21,H22*,ME(3),BI(3),CI(3)ET=EO*T/(1.0-VO*VO)/A/4.0V=(1.0-VO)/2.0H11=ET*(BR*BS+V*CR*CS)H12=ET*(VO*BR*CS+V*BS*CR)H21=ET*(VO*CR*BS+V*BR*CS)H22=ET*(CR*CS+V*BR*BS)RETURNENDC*********************************************SUBROUTINEINPUT(JR,COOR,MEL,AE)DIMENSIONJR(2,*),COOR(2,*),AE(4,*),MEL(3,*)COMMON/CA/NP,NE,NM,NRCOMMON/CC/N,MX,NHDO70I=1,NPREAD(4,*)IP,X,YCOOR(1,IP)=XCOOR(2,IP)=Y70CONTINUEDO11J=1,NEREAD(4,*)NEE,NME,(MEL(I,NEE),I=1,3)MEL(3,NEE)=NME11CONTINUEDO10I=1,NPDO10J=1,210JR(J,I)=1DO20I=1,NRREAD(4,*)IP,IX,IYJR(1,IP)=IXJR(2,IP)=IY20CONTINUEN=0DO30I=1,NPDO30J=1,2IF(JR(J,I))30,30,2525N=N+1JR(J,I)=N30CONTINUEDO55J=1,NMREAD(4,*)jj,(AE(I,jj),I=1,4)CWRITE(*,910)jj