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

平面三角形３结点有限元程序 1、程序名：FEMT3.FOR，FEMT3.EXE 2、程序功能 本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种 荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。 主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结 点的支座反力。 程序采用Visual Fortran 5.0编制而成，输入数据全部采用自由格式。 3、程序流程及框图 图１-1 程序流程图 图１-2 程序图 INPUT——输入结点坐标、单元信息和材料参数； MR——形成结点自由度序号矩阵； FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序； DIV——取出单元的3个结点号码和该单元的材料号并计算单元的bii等； MGK——形成整体劲度矩阵并按一维存放在SK（NH）中； LOAD——形成整体结点荷载列阵F； OUTPUT——输出结点位移或结点荷载； TREAT——由于有非零已知位移，对K和F进行处理； DECOMP——整体劲度矩阵的分解运算； FOBA——前代、回代求出未知结点位移δ； ERFAC——计算约束结点的支座反力； KRS——计算单元劲度矩阵中的子块Krs。 4、输入数据及变量说明 9个数据 NP，NE，NM，NR，NI，NL，NG，ND，NC NP——结点总数； NE——单元总数； NM——材料类型总数； NR——约束结点总数； NI——问题类型标识，0为平面应力问题，1为平面应变问题； NL——受荷载作用的结点的数目； NG——考虑自重作用为1，不计自重为0； ND——非零已知位移结点的数目； NC——要计算支座约束反力的结点数目。 ⑵　材料信息，共NM条，每条依次输入 EO，VO，W，t EO——弹性模量（kN/m2）； VO——泊松比； W——材料容重（kN/m3）； t——单元厚度（m）。 这些信息都存放在数组AE（4，NM）中。 ⑶　坐标信息，共NP条，每条依次输入 IP，X，Y IP——结点号； X，Y——分别为结点的x坐标和y坐标。 坐标信息存放在数组X(2，NP)中。 ⑷　单元信息，共NE条，每条依次输入 JE，L ，Io，Jo，Mo JE——单元号； L——为该单元的材料类型号。 Io，Jo，Mo——分别为该单元i，j，m的整体编码。 单元信息存放在数组MEO(4,NE)中。 ⑸　约束信息共NR条，每条依次输入一个数 ×× IP IP——结点号； Ix，Iy——分别为该结点的约束情况，如果方向受约束时填0，如果自由则填1。 ⑹　荷载信息，共NL条，每条依次输入y IP，Fx，Fy IP——结点号； Fx，Fy——分别为该结点的x，y方向的荷载分量（kN）。 结点号存放在数组NF（NL）中，结点荷载分量存放在数组FV（2，NL）中。 ⑺　若ND＞0，输入非零已知位移信息，共ND条，每条依次输入 IP，ux，uy IP——结点号； ux，uy——分别为该结点x，y方向已知位移分量（m），若其中某方向为自由，则其相应分量为0。 结点号存放在数NDI（ND）中，已知位移分量存放在数组DV(2，ND)中。 ⑻　支座反力信息，共NC条，每条依次输入 IP，M1，M2，M3，M4 IP ——支座结点号； 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所示，设，，不考虑自重。该问题的精确解应力为＝1，＝1，＝0。 图 １-3 有限元网格