平面三角形3节点有限元程序Word格式.docx

1、TREAT由于有非零位移，对K和F进展处理；DEP整体劲度矩阵的分解运算；FOBA前代、回代求出未知结点位移；ERFAC计算约束结点的支座反力；KRS计算单元劲度矩阵中的子块Krs。4、输入数据与变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。数据文件包括如下内容：总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP结点总数；NE单元总数；NM材料类型总数；NR约束结点总数；NI问题类型标识，0为平面应力问题，1为平面应变问题；NL受荷载作用的结点的

2、数目；NG考虑自重作用为1，不计自重为0；ND非零位移结点的数目；NC要计算支座约束反力的结点数目。材料信息，共NM条，每条依次输入EO，VO，W，tEO弹性模量kN/m2；VO泊松比；W材料容重kN/m3；t单元厚度m。信息都存放在数组AE4，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条，每条依次输入一个

3、数IP Ix，Iy分别为该结点的约束情况，如果方向受约束时填0，如果自由如此填1。荷载信息，共NL条，每条依次输入yIP，Fx，FyFx，Fy分别为该结点的x，y方向的荷载分量kN。结点号存放在数组NFNL中，结点荷载分量存放在数组FV2，NL中。假如ND0，输入非零位移信息，共ND条，每条依次输入IP，ux，uyux，uy分别为该结点x，y方向位移分量m，假如其中某方向为自由，如此其相应分量为0。结点号存放在数NDIND中，位移分量存放在数组DV（2，ND）中。支座反力信息，共NC条，每条依次输入IP，M1，M2，M3，M4IP 支座结点号；M1，M2，M3，M4 为与该支座结点相关的单元号

4、，假如不足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）单

5、元劲度矩阵计算公式中的bi，bj，bm和ci，cj，cm；S三角形单元的面积；H11，H12，H21，H22单元劲度矩阵中子块Krs的4个元素。5、算例一个正方形弹性体，厚度为1m，四边受单位均布法向力作用，由于对称性，取其1/4进展计算，其有限元网格如图2-3所示，设，不考虑自重。该问题的准确解应力为1，1，0。图-3 有限元网格1输入文件数据6 41 5 0 3 0 0 52000.0 0.0 0.0 1.01 0.0 2.02 0.0 1.03 1.0 1.04 0.0 0.05 1.0 0.06 2.0 0.01 13 1 2 2 12 4 5 3 13 2 5 4 15 6 3101

6、 201 4005106101 -0.5 -0.53 -1.0 -1.06 -0.5 -0.51 1 0 0 02 1 2 3 04 2 0 0 05 2 3 4 06 4 0 0 02输出文件结果 NODAL DISPLACEMENTSNODE X-P. Y-P. 1 0.00000E+00 -0.10000E-02 2 0.00000E+00 -0.50000E-03 3 -0.50000E-03 -0.50000E-03 4 0.00000E+00 0.00000E+00 5 -0.50000E-03 0.00000E+00 6 -0.10000E-02 0.00000E+00 ELEM

7、ENT STRESSESELEMENT X-STRESS Y-STRESS XY-STRESS MAX-STRESS MIN-STRESS ANGLE 1 -1.000 -1.000 0.000 -1.000 -1.000 90.000 2 -1.000 -1.000 0.000 -1.000 -1.000 90.000 3 -1.000 -1.000 0.000 -1.000 -1.000 90.000 4 -1.000 -1.000 0.000 -1.000 -1.000 90.000 NODE STRESSESNODE X-STRESS Y-STRESS XY-STRESS MAX-ST

8、RESS MIN-STRESS ANGLE 5 -1.000 -1.000 0.000 -1.000 -1.000 90.000 6 -1.000 -1.000 0.000 -1.000 -1.000 90.000 NODAL REACTIONS NODE X-P Y-P 1 0.000 0.000 2 1.000 0.000 4 0.500 0.500 5 0.000 1.000 6 0.000 0.0006、源程序C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONALC TRIANGLE ELEMENTC DIMENSIONK（800000）,COOR（2

9、,3000）,AE（4,11）,* MEL（5,2000）,MA（6000） CHARACTER*32 dat MON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC WRITE（*,300） 300 FORMAT（/ * :*/+PLEASE INPUT FILE NAME OF DATA） READ（*,*）data OPEN （4,FILE=data,STATUS=OLD OPEN （7,FILE=OUT,STATUS=UNKNOWN READ （4,*） NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE （*,400） NP,NE,NM,NR,NI,

10、NL,NG,ND,NC C WRITE （7,400） NP,NE,NM,NR,NI,NL,NG,ND,NC CALL INPUT （JR,COOR,MEL,AE） CALL CBAND （MA,JR,MEL） CALL SK0（SK,R,COOR,MEL,MA,JR,AE）CALL LOAD （COOR,MEL,R,JR,AE） CALL DEP （SK,MA） CALL FOBA （SK,MA,R） WRITE（*,650） WRITE（7,650） CALL OUTPUT（JR,R） WRITE（*,700） WRITE（7,700） CALL CES （COOR,MEL,JR,R,AE

11、） 400 FORMAT （/2X,NP=,I3,2X,NE=NM= *,I3,2X,NR=NI=I3,2X,NL=,I3,2X, *NG=ND=NC=,I3）500 FORMAT（1X,TOTAL DEGREES OF FREEDOM N=, *I4,1X,TOTAL-STORAGE ,NH=,I5,1X,MAX-SEMI-BANDWIDTH MX=550 FORMAT（/20X,TOTAL STORAGE IS *GREATER THAN 50000600 FORMAT（30X,NODAL FORCES/8X,NODE *11X,X-P.,14X,Y-P.650 FORMAT（/30X,N

12、ODAL DISPLACEMENTS/8X,13X,12X,700 FORMAT（/30X,ELEMENT STRESSES/5X,ELEMENT,5X,X-STRESS,3X,Y-STRESS *2X,XY-STRESS,1X,MAX-STRESS,1X,MIN-STRESS,6X,ANGLE/） STOP END C * SUBROUTINE KRS （BR,BS,CR,CS） MON /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.0 V=（1.0-VO）/2.0 H11=ET*（BR*BS+V*CR*CS） H12=ET