有限元程序.docx
- 文档编号:29011607
- 上传时间:2023-07-20
- 格式:DOCX
- 页数:13
- 大小:48.44KB
有限元程序.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受载节点矩阵
ZXX方向已知载荷向量
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)