平面三角形3结点有限元程序.docx

平面三角形3结点有限元程序.docx平面三角形3结点有限元程序1、 程序名:,、 程序功能本程序能计算弹性力学的平而应力问题和平而应变问题;能考虑自重和结点集中力两种荷载的作用,在计算口重吋y轴取垂直向上为止;能处理非零己知位移,如支座沉降的作用。主要输出的内容包括:结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。,输入数据全部采用口由格式。3、 程序流程及框图图1・1程序流程图图1-2程序框图其中,各子程序的功能如下:INPUT——输入结点坐标、单元信息和材料参数;MR——形成结点白由度序号矩阵;FORMMA——形成指标矩阵MA(N)并调用其他功能子程序,和当于主控程序;DIV——取出单元的3个结点号码和该单元的材料号并计算单-元的b,,5等;MGK——形成整体劲度矩阵并按一维存放在SK(NH)中;LOAD——形成整体结点荷载列阵F;OUTPUT——输出结点位移或结点荷载;TREAT——由于有非零己知位移,对K和F进行处理;P——整体劲度矩阵的分解运算;FOBA——前代、回代求出未知结点位移5;ERFAC——计算约束结点的支座反力;KRS——计算单元劲度矩阵中的子块Krso4、输入数据及变量说明当程序开始运行时,按屏幕提示,键入数据文件的名字。在运行程序Z前,必须根据程序中输入要求建立一个存放原始数据的文件,这个文件的名字由少于8个字符或数字组成。数据文件包括如下内容:(1)总控信息,共一条,9个数据NP,NE,NM,NR,NI,NL,NG,ND,NCNP——结点总数;NE——单元总数;NM——材料类型总数;NR——约束结点总数;NI——问题类型标识,0为平]何应力问题,1为平面应变问题;NL——受荷载作用的结点的数目;NG——考虑自重作用为1,不计自重为0;ND——非零已知位移结点的数H;NC——要计算支座约束反力的结点数Uo⑵材料信息,共NM条,每条依次输入EO,VO,W,tEO——弹性模量(kN/n?);VO——泊松比;W——材料容重(kN/m‘);t——单元厚度(m)。这些信息都存放在数组AE(4,NM)中。⑶坐标信息,共NP条,每条依次输入IP,X,YIP——结点号;X,Y 分别为结点的x坐标和y坐标。处标信息存放在数组X(2,NP)中。(4)单元信息,共NE条,每条依次输入JE, L,Io,Jo,MoJE——单元号;L——为该单元的材料类型号。Io,Jo,Mo——分别为该单元几,加的整休编码。单元信息存放在数组MEO(4,NE沖。⑸约束信息共NR条,每条依次输入一个数XXXXX' V '\\\IPI」IP——结点号;I”,I、,——分别为该结点的约束情况,如果方向受约束时填0,如果自由则填1。⑹荷载信息,共NL条,每条依次输入yIP, F”,FyIP——结点号;F、——分别为该结点的兀,y方向的荷载分量(kN)。结点号存放在数组NF(NL)中,结点荷载分量存放在数组FV(2,NL)中。⑺若ND>0,输入非零已知位移信息,共ND条,每条依次输入IP'叮jIP——结点号;uv,uv——分别为该结点x,y方向己知位移分量(m),若其中某方向为口由,则其相应分量为Oo结点号存放在数NDI(ND)中,己知位移分量存放在数组DV(2,ND)中。⑻支座反力信息,共NC条,每条依次输入IP,Ml,M2,M3,M4IP——支座结点号;Ml,M2,M3,M4——为与该支座结点和关的单冗号,若不足4个,则用0补充。支座结点号存放在数组NCI(NC)中,相关单元号存放在数组NCE(4,NC)中。以上数据须按如上顺序存放在数据文件中。除此Z外,程序中还用到其他一些主要变量和数组,说明如下:N——结构自由度总数;NH——按一维存贮的整体劲度矩阵的总容量;NX——最大半带宽;SK(10000)——维存贮的劲度矩阵;R(1000)——开始存放等效结点荷载,求解方程以后,用来存放结点位移;B(6) 存放单元应力ax, ,CT],a2,a:MA(1000)——主元索序号指标矩阵;JR(2,500)——结点自由度序号矩阵;M玖3)——存放单元结点i,j,m的整体编码;NN(6)——单元结点自由度序号;BI(3),S3)——单元劲度矩阵计算公式中的b「,b,,b,”和c,・,c厂c,”;S——三角形单元的面积;Hn,比,H21,H22——单元劲度矩阵中子块Krs的4个元素。5、算例一个止方形弹性休,厚度为1IB,四边受单位均布法向力作用,由于对称性,収其1/4进行计算,其有限元网格如图2・3所示,设E=2-4xl07^/m2,“=0」67,不考虑自重。该问题的精确解应力为crv=1kN/m2,ay=1kN/m2,rvv=0o(a)(c) (d)图1・3有限元网格(1)输入文件数据6