平面四边形4结点等参有限单元法

-1-修改后的等参元有限元程序设计小组成员：-2-平面四边形4结点等参有限单元法程序设计1、程序功能及特点a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。b.能计算受集中力、自重体力、分布面力和静水压力的作用。c.计算结点的位移和单元中心点的应力分量及其主应力。d.后处理采取整体应力磨平求得各个结点的应力分量。e.程序采用VisualFortran5.0编制而成。2、程序流程及图框程序流程图-3-子程序框图其中，各子程序的主要功能为：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‘-4-SUMSTRS――应力磨平-单元下的集成到总体的‘K‘GAUSTRSS――高斯消元求磨平后的应力3、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，根据程序中INPUT需要的数据输入建立一个存放原始数据的文件，这个文件的名字为INDAT.DAT。数据文件包括如下内容：（1）总控信息。共一条，4个数据。L，B，NNL，NNB，NM，NRL——矩形体长度B--矩形体宽度NNL--L方向上划分的结点数NNB--B方向上划分的结点数NM—单元材料类型数NR——约束结点总数（2）结点约束信息。共NR条，每条依次输入：IP，IX，IYIP——结点号IX、IY——分别为IP结点在x，y方向的约束情况，如果约束填0，如果自由填1。（3）材料信息。共NM条，每条依次输入：JJ,（AE（I，JJ），I=1，4）JJ――材料类型号，（AE（I，JJ），I=1，4）――该材料的材料参数，共4个参数，排列顺序为：弹性模量、泊松比、容重、单元厚度。（4）荷载信息a.荷载控制信息。共一条，3个数据NCP，IZNCP——受集中力作用的结点数IZ——面力批数b.若NCP0，输入IP，PX，PYIP——结点号PX、PY——分别为IP结点x，y方向的集中力分量。c.若IZ0，输入面力荷载信息，共IZ批，按批输入：JS，NSE，（WG（I）I=1，4）JS——面力批号-5-NSE——第JS批面力受到面力作用的单元个数，（WG（I），I＝1，4）——该面力的特征参数共4个数据，第1个数为面力类型，填1表示受静水压力作用，填2表示受均布法向压力作用；第2个数为水压密度，如果是均布压力情况，就填均布压力的集度；第3个数为最高水位的y坐标，如果是均布压力情况，可以填任意数；第4个数为面力作用的单元面的面号，单元面号的规定见图2-3。（IEW（M），M＝1，NSE），IEW（*）为受面力作用的单元的单元号，共NSE个。单元结点编码与面号4、算例应用本程序计算一个矩形土体受到均布荷载时的位移和应力，如下图所示，三面约束的土体尺寸为40m*10m，取一半为20m*10m，E=1.0×104kN/m2，0.3，在右端受均布荷载10.0qkN/m2，不考虑自重体力。给定网格大小为0.20.5mm，有限元网格自动划分如图3-2所示，单元总数2000，节点总数2121。由于对称性，右端约束为滑移支座，只限制x方向位移。虽然土体一般不为弹性，但是方法是一样的，只是刚度矩阵改动即可使用弹塑性模型。-6-长度单位：5、输入文件数据：平面四边形四结点单元输入数据LBNNLNNBNMNR20.010.0101211141*********************************************************约束信息数据结点号X向约束信息Y向约束信息10020030040050060070080090010001100120013001400150016001700180019002000-7-21002200430064008500106001270014800169001900021100232002530027400295003160033700358003790040000421004420046300484005050052600547005680058900610006310065200673006940071500736007570077800799008200084100862008830090400-8-92500946009670098800100900103000105100107200109300111400113500115600117700119800121900124000126100128200130300132400134500136600138700140800142900145000147100149200151300153400155500157600159700161800163900166000168100170200172300174400176500178600180700182800-9-184900187000189100191200193300195400197500199600201700203800205900208000210100210201210301210401210501210601210701210801210901211001211101211201211301211401211501211601211701211801211901212001212101*********************************************************材料信息数据类型号弹性模量泊松比容重单元厚度110000.00.30.01.0*********************************************************荷载信息数据受集中力作用的结点数面力批数01集中力数据受集中力作用的结点号X向集中力分量Y向集中力分量*********************************************************-10-面力数据面力批号受面力单元个数面力类型面力集度最大集度面力作用面号1152-10-104*********************************************************面力作用的单元的单元号2000198019601940192019001880186018401820180017801760174017206、源程序：PROGRAMFEMDIMENSIONSK(200000),COOR(2,10000),AE(4,100),MEL(5,10000),WG(4),&JR(2,10000),MA(10000),R(10000),IEW(10000),STRE(3,10000),&TOALFUN(10000,10000),SUMS(10000),GOODSTR(10000)COMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHCOMMON/CMN3/RF(8),SKE(8,8),NN(8)OPEN(1,FILE='INDAT.DAT',STATUS='OLD')OPEN(2,FILE='OUT.DAT',STATUS='NEW')CALLINPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR)!输入数据CALLHUAFEN(XL,B,NNL,NNB,COOR,MEL)CALLCBAND(MA,JR,MEL)!形成主元素序号指示矩阵MADOI=1,NHSK(I)=0.0ENDDOCALLSK0(SK,MEL,COOR,JR,MA,AE)!形成整体刚度矩阵[K]DOI=1,NR(I)=0.0ENDDOIF(NCP.GT.0)CALLCONCR(NCP,R,JR)CALLBODYR(R,MEL,COOR,JR,AE)READ(1,70)TLREAD(1,70)TLREAD(1,70)TL70FORMAT(A)IF(iz.GT.0)thenDOjj=1,izREAD(1,*)JS,NSE,(WG(I),I=1,4)READ(1,70)TLREAD(1,70)TLREAD(1,*)(IEW(M),M=1,NSE)CALLFACER(IEW,NSE,R,MEL,COOR,JR,WG)ENDDO-11-ENDIFCALLDECOP(SK,MA)CALLFOBA(SK,MA,R)CALLOUTDISP(NP,R,JR)CALLSTRESS(COOR,MEL,JR,AE,R,STRE)CALLGAUSTRSS(MEL,COOR,AE,TOALFUN,SUMS,GOODSTR,STRE,JR,R)!应力磨平并输出磨平后的应力值WRITE(2,'(A)')'PROGRAMSAFFHASBEENENDED'WRITE(*,'(A)')'PROGRAMSAFFHASBEENENDED'CLOSE(1)CLOSE(2)END!INPUT为原始数据输入子程序SUBROUTINEINPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR)DIMENSIONAE(4,*),JR(2,*)CHARACTER*400TLCOMMON/CMN1/NP,NE,NM,NRCOMMON/CMN2/N,MX,NHREAD(1,70)TLREAD(1,70)TLREAD(1,*)XL,B,NNL,NNB,NM,NRWRITE(2,10)XL,B,NNL,NNB,NM,NR10FORMAT(5X,'平面四边形四结点单元有限元程序'//5X,'L=',F5.2,4X,'&B=',F5.2,4X,'NNL=',I3,4X,'NNB=',I3,4X,'NM=',I2,4X,'NR=',I3)NP=NNL*NNBDO15I=1,NPDO15J=1,215JR(J,I)=1READ(1,70)TLREAD(1,70)TLREAD(1,70)TLWRITE(2,110)110FORMAT(/5X,'约束信息'/5X,'结点号',5X,'X向约束',5X,'Y向约束')DO20I=1,NRREAD(1,*)IP,IX,IYWRITE(2,100)IP,IX,IY100FORMAT(8X,I5,5X,I2,9X,I2)JR(1,IP)=IXJR(2,IP)=IY20CONTINUEN=0DO30I=1,NPDO30J=1,2-12-IF(JR(J,I))30,30,2525N=N+1JR(J,I)=N30CONTINUEREAD(1,70)TLREAD(1,70)TLREAD(1,70)TLDO55J=1,NMREAD(1,*)JJ,(AE(I,JJ),I=1,4)WRITE(2,910)JJ,(AE(I,JJ),I=1,4)55CONTINUE910FORMAT(/5X,'材料信息'/5X,'类型号',5X,'弹性模量',5X,'泊松比',5X,'容重&',8X,'单元厚度'/5X,I5,4(5x,E8.3))READ(1,70)TLREAD(1,70)TLREAD(1,70)TLREAD(1,*)NCP,IZWRITE(2,920)NCP,IZ920FORMAT(/5X,