有限元程序.docx
《有限元程序.docx》由会员分享，可在线阅读，更多相关《有限元程序.docx（13页珍藏版）》请在冰豆网上搜索。
有限元程序
有限元程序
设
计
报
告
1.程序功能说明
本程序包含四个子程序，其主要功能分别是：
1.MODPS（DMATX，YOUNG，POISS）
该子程序是形成弹性矩阵[D]的，[D]阵的平面应力问题表达式为：
（4-56）
如果要计算平面应变问题，可将上式中弹性模量、泊松比值分别用、值代替。
2、BMATPS（INELE，COORD，LNODE，BMATX，AREA，NPOIN，NELEM）
该子程序用来计算应变矩阵[B]，[B]的表达式为：
（4-57）
3、DBE（DMATX，BMATX，DMATX）
该子程序用来形成应力矩阵[S]，[S]阵的表达式为
（4-58）
4、GUASS（NZERO，LZERO，ASTIF，NHBW，ALOAD，NP2，NPOIT）
该子程序是用来进行支座处理和解方程的，最后输出节点处位移分量。
支座处理的基本原理是如果位移受到约束，则将其对应刚度对角元素变为1，该对角元素所对应的行和列的其余元素变为0，并且将右端项约束位移对应元素充为0。
解方程应用带消去法，解出位移后，将节点位移分量存放在原来存放右端项的标识符ALOAD中，最后将ALOAD（节点位移）输出。
2.程序框图设计
3.程序标识符及数组的说明
1.数组及变量说明
改正COORD节点坐标数组
LNODE单元节点数组
ALOAD荷载（右端项）数组，解出位移后，将节点处的位移充于其中
ESTIF单元刚度矩阵，为6×6阶方阵
ASTIF总刚度矩阵，带状存贮
DBMAT应力矩阵[S]阵，为3×6阶
DMATX弹性矩阵[D]阵，为3×3阶
BMATX应变矩阵[S]阵，为3×6阶
STRES单元应力矩阵，为3×1阶
JJS受载节点矩阵
ZXＸ方向已知载荷向量
ZYY方向已知载荷向量
DISP单元节点位移矩阵
LZERO约束位移矩阵
2.变量的说明
NPOIT最大节点数
NELEM最大单元数
NNODE单元节点数（三角元中为3）
NLOAD最大受载节点数
NZERO最大约束位移个数
YOUNG弹性模量
POISS泊松比
TE板厚（假设t=1）
NHBW半带宽
NP2位移总数
AREA单元面积
4.源程序
FORFINITEELEMENT
DIMENSIONCOORD（6,2）,LNODE（4,3）,ALOAD（12）,ESTIF（6,6）,
1ASTIF（12,8）,DBMAT（3,6）,STRES（3）,BMATX（3,6）,DMATX（3,3）,
1jjS
（1）,ZX
（1）,ZY
（1）,DISP（6）,LZERO（6）
c输入已知数据
OPEN（1,FILE='D:
\NMX\DATA.DAT'）
DATANPOIN,NELEM,NNODE,NLOAD,NZERO/6,4,3,1,6/
DATAYOUNG,POISS,TE/210000.0,0.3,1.0/
DATAZX/0.0/
DATAZY/-1.0/
DATAJJS/1/
DATALNODE/1,4,2,6,2,5,5,3,3,2,3,5/
DATACOORD/0.0,0.0,0.5,0.0,0.5,1.0,1.0,0.5,0.5,0.0,0.0,0.0/
DATALZERO/1,3,7,8,10,12/
C计算半带宽
NHBW=0
DO11INELE=1,NELEM
DO11I=1,NNODE
DO11J=1,NNODE
LN=IABS（LNODE（INELE,I）-LNODE（INELE,J））
IF（LN.GT.NHBW）NHBW=LN
11CONTINUE
NHBW=（NHBW+1）*2
WRITE（1,*）"半带宽"
WRITE（1,12）NHBW
12FORMAT（1x,'NHBW=',I3）
NP2=2*NPOIN
DO50I=1,NP2
ALOAD（I）=0.0
DO50J=1,NHBW
ASTIF（I,J）=0.0
50CONTINUE
C对单元循环
DO70INELE=1,NELEM
CALLMODPS（DMATX,YOUNG,POISS）
CALLBMATPS（INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM）
CALLDBE（DMATX,BMATX,DBMAT）
DO30I=1,6
DO30j=1,6
ESTIF（I,J）=0.0
DO30K=1,3
ESTIF（I,J）=ESTIF（I,J）+DBMAT（K,I）*BMATX（K,J）*AREA*TE
30CONTINUE
DO40ID=1,NNODE
DO40II=1,2
IH=2*（ID-1）+II
IDH=2*（LNODE（INELE,ID）-1）+II
DO35JD=1,NNODE
DO35JJ=1,2
IL=2*（JD-1）+JJ
IDL=2*（LNODE（INELE,JD）-1）+JJ-IDH+1
IF（IDL.LE.0）GOTO35
ASTIF（IDH,IDL）=ASTIF（IDH,IDL）+ESTIF（IH,IL）
35CONTINUE
40CONTINUE
70CONTINUE
C求右端项
DO90I=1,NLOAD
IL=JJS（I）*2
ALOAD（IL-1）=ALOAD（IL-1）+ZX（I）
ALOAD（IL）=ALOAD（IL）+ZY（I）
90CONTINUE
C支座处理、解方程
CALLGAUSS（NZERO,LZERO,ASTIF,NHBW,ALOAD,NP2,NPOIN）
WRITE（1,*）"单元应力"
WRITE（1,*）"单元号σxσyτxy"
DO400INELE=1,NELEM
CALLMODPS（DMATX,YOUNG,POISS）
CALLBMATPS（INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM）
CALLDBE（DMATX,BMATX,DBMAT）
DO410I=1,NNODE
DO410J=1,2
LH=2*（I-1）+J
MH=2*（LNODE（INELE,I）-1）+J
DISP（LH）=ALOAD（MH）
410CONTINUE
DO420I=1,NNODE
STRES（I）=0.0
DO420J=1,6
STRES（I）=STRES（I）+DBMAT（I,J）*DISP（J）
420CONTINUE
WRITE（1,430）INELE,（STRES（I1）,I1=1,NNODE）
430FORMAT（1X,I5,1X,3F13.5,1X）
400CONTINUE
STOP
END
C求弹性矩阵D
SUBROUTINEMODPS（DMATX,YOUNG,POISS）
DIMENSIONDMATX（3,3）
DMATX（1,1）=YOUNG/（1.0-POISS*POISS）
DMATX（1,2）=YOUNG*POISS/（1.0-POISS*POISS）
DMATX（2,1）=DMATX（1,2）
DMATX（2,2）=DMATX（1,1）
DMATX（3,3）=YOUNG/（2.0*（1.0+POISS））
RETURN
END
C求应变矩阵B
SUBROUTINEBMATPS（INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM）
DIMENSIONCOORD（NPOIN,2）,LNODE（NELEM,3）,BMATX（3,6）
IE=LNODE（INELE,1）
JE=LNODE（INELE,2）
ME=LNODE（INELE,3）
BI=COORD（JE,2）-COORD（ME,2）
BJ=COORD（ME,2）-COORD（IE,2）
BM=COORD（IE,2）-COORD（JE,2）
CI=COORD（ME,1）-COORD（JE,1）
CJ=COORD（IE,1）-COORD（ME,1）
CM=COORD（JE,1）-COORD（IE,1）
AREA=（BJ*CM-BM*CJ）/2.0
DO3I=1,3
DO3J=1,6
3BMATX（I,J）=0.0
CH=2.0*AREA
BMATX（1,1）=BI/CH
BMATX（1,3）=BJ/CH
BMATX（1,5）=BM/CH
BMATX（2,2）=CI/CH
BMATX（2,4）=CJ/CH
BMATX（2,6）=CM/CH
BMATX（3,1）=BMATX（2,2）
BMATX（3,2）=BMATX（1,1）
BMATX（3,3）=BMATX（2,4）
BMATX（3,4）=BMATX（1,3）
BMATX（3,5）=BMATX（2,6）
BMATX（3,6）=BMATX（1,5）
RETURN
END
C求应力矩阵DB
SUBROUTINEDBE（DMATX,BMATX,DBMAT）
DIMENSIONDBMAT（3,6）,DMATX（3,3）,BMATX（3,6）
DO3I=1,3
DO3J=1,6
DBMAT（I,J）=0.0
DO3K=1,3
DBMAT（I,J）=DBMAT（I,J）+DMATX（I,K）*BMATX（K,J）
3CONTINUE
END
C支座处理、解方程
SUBROUTINEGAUSS（NZERO,LZERO,ASTIF,NHBW,ALOAD,NP2,NPOIN）
DIMENSIONLZERO（NZERO）,ASTIF（NP2,NHBW）,ALOAD（NP2）
DO260I=1,NZERO
IZ=LZERO（I）
ASTIF（IZ,1）=1.0
DO210J=2,NHBW
ASTIF（IZ,J）=0.0
210CONTINUE
J0=NHBW
IF（IZ-NHBW.LE.0）J0=IZ
DO250J=2,J0
M=IZ-J+1
ASTIF（M,J）=0.0
250CONTINUE
ALOAD（IZ）=0.0
260CONTINUE
KK=NP2-1
DO290K=1,KK
IM=NP2
IF（NP2-K-NHBW+1.GT.0）IM=NHBW+K-1
IK=K+1
DO285I=IK,IM
L=I-K+1
C=ASTIF（K,L）/ASTIF（K,1）
IN=NHBW-L+1
DO280J=1,IN
M=J+I-K
ASTIF（I,J）=ASTIF（I,J）-C*ASTIF（K,M）
280CONTINUE
ALOAD（I）=ALOAD（I）-C*ALOAD（K）
285CONTINUE
290CONTINUE
ALOAD（NP2）=ALOAD（NP2）/ASTIF（NP2,1）
DO315IB=1,KK
I=NP2-IB
J0=NHBW
IF（NHBW-NP2+I-1.GT.0）J0=NP2-I+1
DO310J=2,J0
IH=J+I-1
ALOAD（I）=ALOAD（I）-ASTIF（I,J）*ALOAD（IH）
310CONTINUE
ALOAD（I）=ALOAD（I）/ASTIF（I,1）
315CONTINUE
WRITE（1,*）"节点位移"
WRITE（1,*）"节点号x方向位移y方向位移"
DO111I=1,NPOIN
II=2*I-1
IL=II+1
WRITE（1,580）I,ALOAD（II）,ALOAD（IL）
111CONTINUE
580FORMAT（I5,7X,F10.7,7X,F10.7）
RETURN
END
4.算例
已知一对角受压的正方形薄板，厚度
为1cm，荷载沿厚度均匀分布，为2Kg/cm2,
泊松比
，求板内点的应力与位场。
题中由于XZ面及YZ面均为该板的对称面，所以只选取1/4部分作为计算模型
如图4-12
计算模型的输入数据有以下各量:
NPOINNELEMNNODENLOADNZEROYOUNGPOISSTE
643162100000.31.0
JJS1
ZX0
ZY-1.0
I
COORD123456
COORD（I,1）0.00.00.50.00.01.0
COORD（I,2）1.00.50.50.00.00.0
编号
数组123456
LZERO13781012
LNODEIⅡⅢⅣ
LNODE（I，1）1426
LNODE（I，2）2553
LNODE（I，3）3235
结果的输出
NPOIN
NELEM
NNODE
NLOAD
NZERO
YOUNG
POISS
TE
6
4
3
1
6
210000
0.3
1.0
LNODEDISP-xDISP–y
10.00000E+000.000000E+00
2-3.0954687E+00-9.970900E+00
32.6363658E+00-1.075648963E+00
41.3829005E+00-1.2653849E+00
ELEMENTX-STRY-STRZ-STR
1.38838094E+05-.363800E+-10.108680E+06
2.18348094E+04-.5649400E+05.600930E+05
3-.49598094E+05-.723520E+05-.335670E+05
4.45560238E+05-.728480E+05-.856190E+04
升级会员 