平面问题有限元应力分析程序.docx
- 文档编号:5874464
- 上传时间:2023-01-01
- 格式:DOCX
- 页数:20
- 大小:69.31KB
平面问题有限元应力分析程序.docx
《平面问题有限元应力分析程序.docx》由会员分享,可在线阅读,更多相关《平面问题有限元应力分析程序.docx(20页珍藏版)》请在冰豆网上搜索。
平面问题有限元应力分析程序
第六章有限元程序设计中的若干问题
基本步骤:
ⅰ.结构离散化,输入或生成
结点信息-结点坐标
单元信息-单元结点编号
ⅱ.计算单元刚度矩阵,形成总体刚度矩阵,包括计算
ⅲ.形成结点载荷向量
ⅳ.引入约束条件
ⅴ.解线性方程组
ⅵ.求出结点位移
ⅶ.计算单元的应力并输出
§6-1约束条件的处理
1.对称性与反对称性
(1)对称结构承受对称载荷作用时
(2)对称结构承受反对称载荷作用
2.约束位移的引入
主元置1法
主元赋大值
§6-2总刚度矩阵的存贮法
1.半带宽存贮法
x
x
x
x
O
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
X
X
X
X
X
O
X
2.一维压缩存贮法
考虑到总体刚度矩阵中各行的带宽并不相等,有时由于结构的几何形状的原因,使总体刚度矩阵某些行的带宽特别大。
这种情况下如采用半带宽存贮法,就可能把许多零元素也包含了进去,这对节省计算机的存贮量是很不利的。
一维压缩存贮法是将总体刚度矩阵的下三角形中每一行从第一个非零元素开始按行将元素排成一序列,存放于一维数组
中。
但是为了确定SK中的元素在[K]中的行列号,还需要将[K]中各行对角线的元素在伊维数组中的序号存放于另一辅助数组KD(N2)中(N2是总刚度矩阵的阶数)。
现举例说明这一存贮法:
设有一系数阵
在一维数组SK(13)中依次存放的是
而辅助数组KD(6)中存放的是
KD(6)其实就是[K]中对角元素在一维数组SK(13)中的地址。
将一结构离散化后,对结点进行编号,就能依据单元号确定出总刚度矩阵[K]各行的带宽,由它依次累加就可得出其对角线元素一维存贮中的序号。
显然,形成了数组Kd,就确定了[K]中被存贮的元素分布情况以及SK和[K]中元素的对应关系,例如可求出[K]中第I行带宽为
也可确定出[K]中第I行左边第一个非零元素在[K]中的列号
此外,也能立即确定出单元刚度矩阵
中的某子矩阵
组集到一维数组存贮总刚度矩阵SK中的地址
,
,
结构总刚度矩阵的以上两种存贮方法,一般应用于直接解法中。
附录A
平面问题有限元应力分析程序
本程序是用FORTRAN语言编写的教学使用程序,采用常应变三角形单元,用于解决弹性力学平面应力问题.它极易由学员扩充为求解平面应力和平面应变问题的通用应力分析程序。
程序中总刚度矩阵按一维压缩存贮,线性代数方程组用三角分解直接法求解。
为使教学程序易读易懂,计算时输入计算机的载荷是结点载荷,任何其他载荷都要事先换算成等效的结点载荷.程序所处理的约束仅是X或Y轴方向上的零位移约束。
计算结果除输出结点位移外,还输出单元形心处的σx,σy,
τxy,主应力σ1,σ2及主方向角θ。
程序的结构便于修改和扩充,易于连接图形库扩充为包含前、后处理程序(网格自动生成及计算结果的图形显示)的更完整的程序系统.
1.输入数据说明(以输入的先后为序,自由格式)
(1)NN,NE,KU,KV,KRX,KRY,EO,PO————6个整型数,2个实型数,其中
NN结点总数(≤500);
NE单元总数(≤700);
KUx方向位移受约束的结点数(≤50);
KVy方向位移受约束的结点数(≤50);
KRX在x方向受结点载荷作用的结点数(≤60);
KRY在y方向受结点载荷作用的结点数(≤60);
EO材料的弹性模量;
PO材料的泊桑比;
(2)X(NN)————NN个实型数,为结点的x坐标.
(3)Y(NN)————NN个实型数,为结点的y坐标.
(4)IJM(NE,3)————3*NE个整型数,按行输入,为单元按逆时针向的结点编号
(5)JU(KU)————KU个整型数,为x方向位移受约束的结点号.
(6)JV(KV)————KV个整型数,为y方向位移受约束的结点号.
(7)NRX(KRX)————KRX个整型数,为在x方向受结点载荷作用的结点号.
(8)NRY(KRY)————KRY个整型数,为在y方向受结点载荷作用的结点号.
(9)RX(KRX)————KRX个实型数,为x方向结点载荷的大小.
(10)RY(KRY)————KRY个实型数,为y方向结点载荷的大小.
2.其他标识符说明
NF(=2*NN)方程的总数;
LK总刚度矩阵下三角形一维存贮的总长;
D1,D2,D3弹性矩阵中的元素;
B(3),C(3)Bi,Ci(i=1~3)的工作单元;
DEL三角形单元面积的工作单元;
U(NF)开始存放结点载荷,求解后存放结点位移;
SK(LK)按一维存贮的总刚度矩阵;
EK(21)单元刚度矩阵下三角形元素按一维存贮的数组;
BM(3,6)用于存放2Δ〔B〕矩阵的元素;
CM(3,6)用于存放2Δ﹝D﹞〔B〕矩阵的元素;
JD(NF)总刚度矩阵下三角形对角线元素在一维存贮中的序号指示数组;
JLL(6)单元刚度矩阵元素和总刚度矩阵元素行列对应关系的工作数组;
SGM(3)存放单元形心处的应力σx,σy,τxy;
SGM1,SGM2存放单元形心处的主应力σ1,σ2;
CETA存放主应力的方向角.
3.各子程序功能
SKDD计算总刚度矩阵下三角形对角线元在一维存
贮中的序号
SHAPE(N)用于计算第N个单元的有关常数.
FEK计算单元刚度矩阵,并存放于EK(21)中.
SKKE单元刚度矩阵向总刚度矩阵传送.
FIXD用于处理约束条件.
SOLVE解线性代数方程组.
STRESS用于计算单元形心处的应力,主应力和主
向,并输出.
4.出错信息的处理
本程序中若干数组是半动态数组,如求解问题的规模超过所设
数组下标的界限,则停机,等待处理
(1)求解问题网格中的结点数大于500时,则停机.此时需将数组
X(500),Y(500),U(1000),JD(1000)下标的大小作相应的扩大.
(2)求解问题网格中的单元数大于700时,则停机.此时需将
JE(700,3)数组中的前一个下标的大小作相应的扩大.
(3)当求解问题的x方向或y方向的约束自由度数超过50时,则
停机.此时需将JU(50)或JV(50)下标的大小作相应的扩充.
(4)当求解问题的x方向或y方向的结点载荷数超过60时,则停
机.此时需将NRX(60),RX(60)或NRY(60),RY(60)下标的大小作相应的扩充.
(5)当一维存贮的总刚度矩阵长度大于12000时,则停机.此时应
将数组SK(12000)下标的大小作相应的扩大.
5.主程序框图
开始
读入原始数据
调用SKDD,计算总刚度矩阵存贮信息
输出原始数据供校对,保存
形成结点载荷向量
调用SHAPE和FEK,计算单元刚度矩阵EKEK(21)
结束
调用SHAPE和STRESS计算单元形心处应力、主应力和主方向,并输出
输出结点位移
调用SOLVE,解方程组
调用FIXD,处理约束
调用SKKE,将EK(21)送总刚度矩阵SK(LK)
6.源程序文本
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CP.FEM---ASTRESSANALYSISFINITEELEMENTC
CPROGRAMFORPLANEPROBLEMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CMAINPROGRAMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/CONT/NN,NE,NF,LK
COMMON/MAT/D1,D2,D3
COMMON/ELEM/B(3),C(3),DEL
COMMON/DATA1/X(500),Y(500),U(1000),JE(700,3)
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
DIMENSIONJU(50),JV(50),NRX(60),NRY(60),RX(60),RY(60)
OPEN(7,FILE=′DATA1′,STATUS=′OLD′)
OPEN(8,FILE=′DATA2′,STATUS=′NEW′)
READ(7,*)NN,NE,KU,KV,KRX,KRY,EO,PO
IF(NN.LE.500)GOTO10
WRITE(*,901)
STOP
10IF(NE.LE.700)GOTO20
WRITE(*,902)
STOP
20IF(KU.LE.50)GOTO30
WRITE(*,903)
STOP
30IF(KV.LE.50)GOTO40
WRITE(*,904)
STOP
40IF(KRX.LE.60)GOTO50
WRITE(*,905)
STOP
50IF(KRY.LE.60)GOTO55
WRITE(*,906)
STOP
55READ(7,*)(X(I),I=1,NN)
READ(7,*)(Y(I),I=1,NN)
READ(7,*)((JE(I,J),J=1,3),I=1,NE)
READ(7,*)(JU(I),I=1,KU)
READ(7,*)(JV(I),I=1,KV)
READ(7,*)(NRX(I),I=1,KRX)
READ(7,*)(NRY(I),I=1,KRY)
READ(7,*)(RX(I),I=1,KRX)
READ(7,*)(RY(I),I=1,KRY)
NF=2*NN
CALLSKDD
WRITE(8,1000)
WRITE(8,1100)NN,NE,KU,KV,KRX,KRY,EO,PO
WRITE(8,*)′X(NN)=′
WRITE(8,1200)(X(I),I=1,NN)
WRITE(8,*)′Y(NN)=′
WRITE(8,1200)(Y(I),I=1,NN)
WRITE(8,*)′JE(NE,3)=′
WRITE(8,1300)((JE(I,J),J=1,3),I=1,NE)
WRITE(8,*)′JU(KU)=′
WRITE(8,1300)(JU(I),I=1,KU)
WRITE(8,*)′JV(KV)=′
WRITE(8,1300)(JV(I),I=1,KV)
WRITE(8,*)ˊNRX(KRX)=ˊ
WRITE(8,1300)(NRX(I),I=1,KRX)
WRITE(8,*)ˊNRY(KRY)=ˊ
WRITE(8,1300)(NRY(I),I=1,KRY)
WRITE(8,*)ˊRX(KRX)=ˊ
WRITE(8,1200)(RX(I),I=1,KRX)
WRITE(8,*)ˊRY(KRY)=ˊ
WRITE(8,1200)(RY(I),I=1,KRY)
D1=EO/(1.-PO*PO)
D2=D1*PO
D3=D1*(1.-PO)/2
DO60I=1,LK
60SK(I)=0.
DO100I=1,NF
100U(I)=0.
DO120I=1,KRX
K1=2*NRX(I)-1
120U(K1)=RX(I)
DO140I=1,KRY
K1=2*NRY(I)
140U(K1)=RY(I)
DO160M=1,NE
WRITE(*,*)ˊELE=ˊ,M
CALLSHAPE(M)
CALLFEK
160CALLSKKE
DO180I=1,KU
180CALLFIXD(2*JU(I)-1)
DO200I=1,KV
200CALLFIXD(2,*JV(I))
CALLSOLVE
WRITE(*,*)ˊSOLVEPASSEDˊ
WRITE(8,3000)
WRITE(8,4000)(I,U(2*I-1),U(2*I),I=1,NN)
WRITE(8,5000)
DO240M=1,NE
CALLSHAPE(M)
CALLSTRESS(M)
240CONTINUE
STOP
901FORMAT(5X,ˊNUMBEROFNODESEXCEEDSˊ
*ˊALLOWABLEˊ)
902FORMAT(5X,ˊNUMBEROFELEMENTˊ
*ˊEXCEEDAALLOWABLEˊ)
903FORMAT(5X,ˊNUMBEROFFIXEDFREEDOMATX-ˊ,
*ˊDIRECTIONEXCEEDSALLOWABLEˊ)
904FORMAT(5X,ˊNUMBEROFFIXEDFREEDOMATY-ˊ,
*ˊDIRECTIONEXCEEDSALLOWABLEˊ)
905FORMAT(5X,ˊNUMBEROFNODALFORCESATX-ˊ,
*ˊDIRECTIONEXCEEDSALLOWABLEˊ)
906FORMAT(5X,ˊNUMBEROFNODALFORCESATY-ˊ,
*ˊDIRECTIONEXCEEDSALLOWABLEˊ)
1000FORMAT(//10X,ˊPRIMARYDATAˊ//,1X,ˊCONˊ,
*ˊSTANTS=ˊ)
1100FORMAT(6I5,2F12.2)
1200FORMAT(5F12.4)
1300FORMAT(15I5)
3000FORMAT(//10X,ˊNODALDISPLACEMENTˊ/,
*ˊNODEˊ,9X,ˊUˊ,13X,ˊVˊ,4X,ˊNODEˊ,9X,
*ˊUˊ,13X,ˊVˊ/ˊ)
4000FORMAT(2(I5,2E14.4))
5000FORMAT(//2X,ˊELE.ˊ,2X,ˊELEMENTSTESSˊ,/8X,
*ˊSGMXˊ,8X,ˊSGMYˊ,9X,ˊTXYˊ,8X,ˊSGM1ˊ,
*8X,ˊSGM2ˊ,10X,ˊCETAˊ/)
END
C
C
SUBROUTINESKDD
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCALCULATEORDERSOFDIAGONALC
CELEMENTSOFMATRIXKC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/CONT/NN,NE,NF,LK
COMMON/DATA1/X(500),Y(500),U(10001),JE(700,
3)
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
JD
(1)=1
JD
(2)=3
DO200N=2,NN
KK=0
DO100I=1,NE
IO=JE(I,1)
JO=JE(I,2)
KO=JE(I,3)
IF(IO.NE.N.AND.JO.NE.N.AND.KO.NE.N)GOTO100
M1=N-IO
M2=N-JO
M3=N-KO
IF(M1.GT.KK)KK=M1
IF(M2.GT.KK)KK=M2
IF(M3.GT.KK)KK=M3
100CONTINUE
KK=2*KK
JD(2*N-1)=JD(2*N-2)+KK+1
JD(2*N)=JD(2*N-1)+KK+2
200CONTINUE
LK=JD(NF)
WRITE(*,*)ˊLK=ˊ,LK
IF(LK.LE.12000)GOTO210
WRITE(*,1000)
STOP
1000FORMAT(5X,ˊNUMBEROFELEMENTSOFONE-DIˊ,
*ˊMENSIONALGLOBALMATRIXEXCEEDSALLOWˊ,
*ˊABLEˊ)
210RETURN
END
C
C
SUBROUTINESHAPE(N)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CEVALUATEELEMENTPARAMETERSC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/DATA1/X(500),Y(500),U(1000),JE(700,3)
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
COMMON/ELEM/B(3),C(3),DEL
DIMENSIONXO(6),YO(6)
DO100I=1,3
K1=JE(N,I)
K2=I+3
XO(I)=X(K1)
XO(K2)=XO(I)
YO(I)=Y(K1)
YO(K2)=YO(I)
K2=2*I-1
K3=K2+1
JLL(K2)=K1*2-1
100JLL(K3)=K1*2
DO200I=1,3
K1=I+1
K2=I+2
B(I)=YO(K1)-YO(K2)
200C(I)=XO(K2)-XO(K1)
DEL=(B
(1)*C
(2)-B
(2)*C
(1))/2.
RETURN
END
C
C
SUBROUTINEFEK
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCALCULATESTIFFNESSMATRIXC
COFELEMENTC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
COMMON/ELEM/B(3),C(3),DEL
COMMON/MAT/D1,D2,D3
DIMENSIONBM(3,6),CM(3,6)
DO100I=1,3
DO100J=1,6
100BM(I,J)=0.
DO200I=1,3
K2=I+I
K1=K2-1
BM(1,K1)=B(I)
BM(3,K2)=B(I)
BM(2,K2)=C(I)
BM(3,K1)=C(I)
CM(1,K1)=D1*B(I)
CM(1,K2)=D2*C(I)
CM(2,K1)=D2*B(I)
CM(2,K2)=D1*C(I)
CM(3,K1)=D3*C(I)
200CM(3,K2)=D3*B(I)
DO400I=1,6
K1=I*(I-1)/2
DO400II=1,I
Z=0.
DO300JJ=1,3
300Z=Z+BM(JJ,I)*CM(JJ.II)
EK(K1+II)=Z/DEL/4
400CONTINUE
RETURN
END
C
C
SUBROUTINESKKE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CASSEMBLEGLOBALSTIFFNESSMATRIXC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
DO200I=1,6
K1=JLL(I)
KK=I*(I-1)/2
DO200J=1,I
K2=JLL(J)
IF(K2.LT.K1)GOTO100
K=JD(K2)-K2+K1
GOTO200
100K=JD(K1)-K1+K2
200SK(K)=SK(K)+EK(KK+J)
RETURN
END
C
C
SUBROUTINEFIXD(K)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CREADBOUNDARYCONDITIONSC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
COMMON/CONT/NN,NE,NF,LK
COMMON/STIF/SK(12000),EK(21),JD(1000),JLL(6)
COMMON/DATA1/X(500),Y(500),U(1000),JE(700,3)
L=JD(K)
IF(K.EQ.1)GOTO200
M=JD(K-1)
IA=M+1
IB=L-1
DO100I=IA,IB
100SK(I)=0.
200IA=K+1
DO300I=IA,NF
IF(JD(I)-JD(I-1).LT.I-K+1)GOTO300
M=JD(I)-I+K
SK(M)=0.
300CONTINUE
U(K)=0
RETURN
END
C
C
SUBROUTINESOLVE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CSOLVEFORSOLUTIONOFEQUATIONSC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 平面 问题 有限元 应力 分析 程序