计算力学基础有限元程序.docx
- 文档编号:4416266
- 上传时间:2022-12-01
- 格式:DOCX
- 页数:19
- 大小:69.09KB
计算力学基础有限元程序.docx
《计算力学基础有限元程序.docx》由会员分享,可在线阅读,更多相关《计算力学基础有限元程序.docx(19页珍藏版)》请在冰豆网上搜索。
计算力学基础有限元程序
计算力学有限元分析作业报告
机械设计及理论张冬明2111102138
如图1所示正方形薄板,沿其对角承受一对压力,载荷沿其厚度为均匀分布,设板的杨氏模量E=200GPa,泊松比v=0.3,厚度t=5mm。
由于薄板形状及所受载荷的对称性,故取其1/4作有限元分析。
有限元模型如图2所示,由6个节点,4个单元组成,节点坐标及单元编号如表1、2所示。
单元编号及边界条件所图3所示。
节点1和节点2在x方向上位移为0,节点5,和节点6在y方向上位移为0,节点4在x和y方向上位移均为0。
节点1受到沿y轴负方向100KN作用力。
图1正方形薄板
表1节点编号及坐标
NODE
X
Y
Z
1
0
2
0
2
0
1
0
3
1
1
0
4
0
0
0
5
1
0
0
6
2
0
0
表2单元编号及节点
ELEM
NODE
NODE
NODE
1
1
2
3
2
2
4
5
3
2
5
3
4
3
5
6
图2有限元模型
图3载荷及边界条件
程序说明:
本程序计算弹性力学平面应力问题(带厚度),包含主程序STRESS和四个子程序ELSTMX(KK)、MODIFY、DCMPBD和SLVBD。
主程序主要进行数据输入、结果输出、最大半带宽的计算、组装总体刚度矩阵以及计算单元内的应力、应变值。
子程序ELSTMX(KK)计算刚度系数举证,生成KK单元的刚度举证K,子程序MODIFY输入载荷节点处的载荷值、位移边界点处的位移值,对总体刚度矩阵、位移数组和节点力数组作修改,子程序DCMPBD用高斯消元法将对称等带宽总刚度矩阵的上半部分化为上三角矩阵,子程序SLVBD对DCMPBD得到的上三角矩阵进行回代计算,得到节点位移向量DD。
程序代码见附录A。
输入文件内容
输入信息包括节点数、单元数、弹性模量、泊松比、平板厚度、节点坐标、单元号及其节点编号、载荷及边界约束条件。
具体参见附录B。
计算输出结果
输出信息列于OUTPUT.DAT文件内,具体包括模型信息、节点位移及单元应力应变。
具体参见附录C。
ANSYS软件分析节点及单元形变如图4所示。
分析结果节点位移值如表3所示,与程序计算结果向吻合。
图4模型形变
表3ANSYS分析结果节点位移
NODE
UX
UY
1
0
-3.28E-04
2
0
-1.38E-04
3
2.65E-05
-3.61E-05
4
0
0
5
5.62E-05
0
6
6.70E-05
0
附录A:
STRESS程序
PROGRAMSTRESS
COMMON/ELMATX/ESM(6,6),X(3),Y(3),D(3,3)
COMMON/GRAD/B(3,6),AR2
COMMON/MTL/EM,PR,TH
COMMON/AV/A(8500),JGF,JGSM,NP,NBW,JEND
DIMENSIONNS(6),U(6),STRA(3),STRE(3)
DATAIN/60/,IO/61/
REAL,ALLOCATABLE:
:
XC(:
),YC(:
)
INTEGER,ALLOCATABLE:
:
NEL(:
:
)
CHARACTER(LEN=40)TITLE
OPEN(60,FILE='INPUT.DAT',STATUS='UNKNOWN')
OPEN(61,FILE='OUTPUT.DAT',STATUS='UNKNOWN')
READ(IN,1)TITLE
1FORMAT(900A)
WRITE(IO,1)TITLE
WRITE(IO,*)
!
输入计算模型的节点数NN和单元总数NE
READ(IN,*)NN
READ(IN,*)NE
WRITE(IO,2)NN,NE
2FORMAT(/,'输入数据为:
',/,'节点数=',I3,5x,'单元数=',I4,/)
NP=2*NN
ALLOCATE(NEL(1:
NE,1:
3),XC(NN),YC(NN))
!
输入材料的杨氏模量EM,泊松比PR,平板厚度TH,节点坐标XC(I)
!
YC(I)和组成单元的节点NEL(N,I)
READ(IN,*)EM
READ(IN,*)PR
READ(IN,*)TH
READ(IN,*)(XC(I),I=1,NN)
READ(IN,*)(YC(I),I=1,NN)
WRITE(IO,3)
3FORMAT('材料常数为:
')
WRITE(IO,4)EM,PR,TH
4FORMAT('弹性模量为:
',E12.5,5X,'泊松比为:
',E12.5,5X,'厚度为:
',E12.5)
WRITE(IO,*)
DO6I=1,NN
WRITE(IO,5)I,XC(I),YC(I)
5FORMAT('节点',I3,'坐标为:
X=',F8.3,3X,'Y=',F8.3)
6CONTINUE
WRITE(IO,*)
DO7KK=1,NE
READ(IN,*)N,(NEL(N,I),I=1,3)
WRITE(IO,8)N,(NEL(N,I),I=1,3)
8FORMAT('单元号码为:
',I3,'组成单元的节点号码为:
',I3,I3,I3)
7CONTINUE
WRITE(IO,*)
!
计算开始
!
计算最大半带宽,B=MAXe(De+1)*F
INBW=0
NBW=0
DO20KK=1,NE
DO25I=1,3
25NS(I)=NEL(KK,I)
DO21I=1,2
IJ=I+1
DO21J=IJ,3
NB=IABS(NS(I)-NS(J))
IF(NB.LE.NBW)GOTO21
INBW=KK
NBW=NB
21CONTINUE
20CONTINUE
NBW=(NBW+1)*2
!
初始化数组A
JGF=NP
JGSM=JGF+Np
JEND=JGSM+NP*NBW
write(*,*)jend
JL=JEND-JGF
DO24I=1,JEND
24A(I)=0.0
GOTO30
!
生成材料性质矩阵D
30R=EM/(1.0-PR**2)
D(1,1)=R
D(2,2)=D(1,1)
D(3,3)=R*(1.0-PR)/2
D(1,2)=PR*R
D(2,1)=D(1,2)
D(1,3)=0.0
D(3,1)=0.0
D(2,3)=0.0
D(3,2)=0.0
!
单元矩阵循环的开始
KK=1
!
节点自由度的生成,节点坐标的局部化
32DO31I=1,3
J=NEL(KK,I)
NS(2*I-1)=J*2-1
NS(2*I)=J*2
X(I)=XC(J)
31Y(I)=YC(J)
!
调用子程序ELSTMX计算单元刚度矩阵ESM(6,6)并输出
CALLELSTMX(KK)
!
单元刚度矩阵组装成总体刚度矩阵
DO33I=1,6
II=NS(I)
DO34J=1,6
JJ=NS(J)+1-II
IF(JJ.LE.0)GOTO34
J1=JGSM+(JJ-1)*NP+II-(JJ-1)*(JJ-2)/2
A(J1)=A(J1)+ESM(I,J)
34CONTINUE
33CONTINUE
KK=KK+1
IF(KK.LE.NE)GOTO32
!
调用子程序MODIFY输入载荷节点处的载荷值、位移边界节点处的位移值
!
对总体刚度矩阵、位移数组和节点力数组进行相应的修改
CALLMODIFY
WRITE(IO,35)
35FORMAT(/,'计算结果为:
',/)
WRITE(IO,36)NBW
36FORMAT(/,'最大半带宽为:
',I3)
WRITE(IO,37)JEND
37FORMAT(/,'总数组大小为:
',I3)
!
调用子程序DCMPBD,用高斯消元法将对称等带宽总刚度矩阵化为上三角阵
CALLDCMPBD
!
调用子程序SLVBD对子程序DCMPBD得到的上三角矩阵进行回代计算,得到节点位移向量DD
CALLSLVBD
!
输出各节点的位移向量
WRITE(IO,*)
DO45I=1,NP/2
WRITE(IO,43)I,A(2*I-1),A(2*I)
43FORMAT('节点号',I3,5X,'X方向位移UX=',E12.5,5X,'Y方向位移UY=',E12.5)
45CONTINUE
!
附加计算
DO96KK=1,NE
!
节点自由度的生成,节点坐标的局部化
DO51I=1,3
J=NEL(KK,I)
NS(2*I-1)=J*2-1
NS(2*I)=J*2
X(I)=XC(J)
51Y(I)=YC(J)
!
单元节点位移局部化
65DO73I=1,6,2
NS1=NS(I)
NS2=NS(I+1)
U(I)=A(NS1)
73U(I+1)=A(NS2)
!
计算矩阵B
DO200I=1,3
DO200J=1,6
200B(I,J)=0.0
B(1,1)=Y
(2)-Y(3)
B(1,3)=Y(3)-Y
(1)
B(1,5)=Y
(1)-Y
(2)
B(2,2)=X(3)-X
(2)
B(2,4)=X
(1)-X(3)
B(2,6)=X
(2)-X
(1)
B(3,1)=B(2,2)
B(3,2)=B(1,1)
B(3,3)=B(2,4)
B(3,4)=B(1,3)
B(3,5)=B(2,6)
B(3,6)=B(1,5)
AR2=X
(2)*Y(3)+X(3)*Y
(1)+X
(1)*Y
(2)-X
(2)*Y
(1)-X(3)*Y
(2)-X
(1)*Y(3)
!
计算单元应变
DO52I=1,3
STRA(I)=0.0
DO52K=1,6
52STRA(I)=STRA(I)+B(I,K)*U(K)/AR2
!
计算单元应力值
DO58I=1,3
STRE(I)=0
DO58K=1,3
58STRE(I)=STRE(I)+D(I,K)*(STRA(K))
!
计算主应力S1,S2,TM
AA=(STRE
(1)+STRE
(2))/2.0
AB=SQRT((ABS(STRE
(1)-STRE
(2))/2.0)**2+STRE(3)**2)
S1=AA+AB
S2=AA-AB
TM=AB
!
计算主应力方向与X轴的夹角
IF(ABS(STRE
(1)-STRE
(2)).LT.0.0001)GOTO93
AC=ATAN2(2*STRE(3),(STRE
(1)-STRE
(2)))
THM=(180/3.1415926*AC)/2
GOTO94
93THM=90
!
输出单元应变和应力的计算结果
94WRITE(IO,57)KK
57FORMAT(/,'单元',I4)
WRITE(IO,95)STRA
(1),STRA
(2),STRA(3)
95FORMAT('EPTOX=',E12.5,2X,'EPTOY=',E12.5,2X,'EPTOXY=',E12.5)
WRITE(IO,97)STRE
(1),STRE
(2),STRE(3)
97FORMAT('SX=',E12.5,5X,'SY=',E12.5,5X,'SXY=',E12.5)
WRITE(IO,98)S1,S2,TM
98FORMAT('S1=',E12.5,5X,'S2=',E12.5,5X,'TMAX=',E12.5)
WRITE(IO,99)THM
99FORMAT(44X,'ANGEL=',F8.2,'°')
96CONTINUE
CLOSE(60)
CLOSE(61)
STOP
END
!
计算刚度系数矩阵,计算第KK单元的刚度矩阵K
SUBROUTINEELSTMX(KK)
COMMON/MTL/EM,PR,TH
COMMON/ELMATX/ESM(6,6),X(3),Y(3),D(3,3)
COMMON/GRAD/B(3,6),AR2
DIMENSIONC(6,3)
DATAIN/60/,IO/61/
!
计算矩阵B
DO20I=1,3
DO20J=1,6
20B(I,J)=0.0
B(1,1)=Y
(2)-Y(3)
B(1,3)=Y(3)-Y
(1)
B(1,5)=Y
(1)-Y
(2)
B(2,2)=X(3)-X
(2)
B(2,4)=X
(1)-X(3)
B(2,6)=X
(2)-X
(1)
B(3,1)=B(2,2)
B(3,2)=B(1,1)
B(3,3)=B(2,4)
B(3,4)=B(1,3)
B(3,5)=B(2,6)
B(3,6)=B(1,5)
AR2=X
(2)*Y(3)+X(3)*Y
(1)+X
(1)*Y
(2)-X
(2)*Y
(1)-X(3)*Y
(2)-X
(1)*Y(3)
!
计算矩阵C
DO22I=1,6
DO22J=1,3
C(I,J)=0.0
DO22K=1,3
22C(I,J)=C(I,J)+B(K,I)*D(K,J)
!
计算矩阵ESM=CB
DO27I=1,6
DO27J=1,6
SUM=0.0
DO28K=1,3
28SUM=SUM+C(I,K)*B(K,J)
ESM(I,J)=SUM*TH/(2.0*AR2)
27CONTINUE
RETURN
END
!
输入载荷节点处的载荷值、位移边界节点处的位移值,对总体刚度矩阵、位移数组和节点力数组进行修改
SUBROUTINEMODIFY
COMMON/AV/A(8500),JGF,JGSM,NP,NBW,JEND
DATAIN/60/,IO/61/
!
输入节点的集中载荷,并放到A数组中节点力向量的相应位置
!
IB为施加外载荷节点的自由度,BV为该自由度上的载荷值
202READ(IN,*)IB
IF(IB.LE.0)GOTO208
READ(IN,*)BV
IF(MOD(IB,2).EQ.1)GOTO204
WRITE(IO,203)IB/2,BV
203FORMAT('节点',I3,'载荷为:
PY=',F8.3)
GOTO206
204WRITE(IO,205)(IB+1)/2,BV
205FORMAT('节点',I3,'载荷为:
PX=',F8.3)
206A(JGF+IB)=A(JGF+IB)+BV
GOTO202
!
输入节点的位移值,并放到A数组中节点力向量的相应位置
208READ(IN,*)IB
IF(IB.LE.0)RETURN
READ(IN,*)BV
IF(MOD(IB,2).EQ.1)GOTO214
WRITE(IO,213)IB/2,BV
213FORMAT('节点',I3,'位移约束为:
V=',F8.3)
GOTO209
214WRITE(IO,215)(IB+1)/2,BV
215FORMAT('节点',I3,'位移约束为:
PX=',F8.3)
209K=IB-1
DO211J=2,NBW
M=IB+J-1
IF(M.GT.NP)GOTO210
IJ=JGSM+(J-1)*NP+IB-(J-1)*(J-2)/2
A(JGF+M)=A(JGF+M)-A(IJ)*BV
A(IJ)=0.0
210IF(K.LE.0)GOTO211
KJ=JGSM+(J-1)*NP+K-(J-1)*(J-2)/2
A(JGF+K)=A(JGF+K)-A(KJ)*BV
A(KJ)=0.0
K=K-1
211CONTINUE
A(JGF+IB)=A(JGSM+IB)*BV
221CONTINUE
GOTO208
END
!
高斯消元子程序
SUBROUTINEDCMPBD
COMMON/AV/A(8500),JGF,JGSM,NP,NBW,JEND
NP1=NP-1
DO226I=1,NP1
MJ=I+NBW-1
IF(MJ.GT.NP)MJ=NP
NJ=I+1
MK=NBW
IF((NP-I+1).LT.NBW)MK=NP-I+1
ND=0
DO225J=NJ,MJ
MK=MK-1
ND=ND+1
NL=ND+1
DO225K=1,MK
NK=ND+K
JK=JGSM+(K-1)*NP+J-(K-1)*(K-2)/2
INL=JGSM+(NL-1)*NP+I-(NL-1)*(NL-2)/2
INK=JGSM+(NK-1)*NP+I-(NK-1)*(NK-2)/2
II=JGSM+I
225A(JK)=A(JK)-A(INL)*A(INK)/A(II)
226CONTINUE
RETURN
END
!
回代求解子程序
SUBROUTINESLVBD
COMMON/AV/A(8500),JGF,JGSM,NP,NBW,JEND
DATAIN/60/,IO/61/
NP1=NP-1
!
对子程序DCMPBD得到的上三角矩阵进行回代计算,得到节点位移向量DD
DO250I=1,NP1
MJ=I+NBW-1
IF(MJ.GT.NP)MJ=NP
NJ=I+1
L=1
DO250J=NJ,MJ
L=L+1
IL=JGSM+(L-1)*NP+I-(L-1)*(L-2)/2
250A(JGF+J)=A(JGF+J)-A(IL)*A(JGF+I)/A(JGSM+I)
!
求解节点自由度的回代计算,从下向上迭代
A(NP)=A(JGF+NP)/A(JGSM+NP)
DO252K=1,NP1
I=NP-K
MJ=NBW
IF((I+NBW-1).GT.NP)MP=NP-I+1
SUM=0.0
DO251J=2,MJ
N=I+J-1
IJ=JGSM+(J-1)*NP+I-(J-1)*(J-2)/2
251SUM=SUM+A(IJ)*A(N)
252A(I)=(A(JGF+I)-SUM)/A(JGSM+I)
RETURN
END
附录B:
输入文件INPUT内容
计算对角受压正方形薄板的应力分布
6!
节点数
4!
单元数
200000000000!
弹性模量
0.3!
泊松比
0.005!
厚度
001012!
节点坐标
211000
1123!
单元及其节点号
2245
3253
4563
2!
载荷
-100000
0
1!
约束条件
0
3
0
7
0
8
0
10
0
12
0
0
附录C:
输出文件OUTPUT内容
计算对角受压正方形薄板的应力分布
输入数据为:
节点数=6单元数=4
材料常数为:
弹性模量为:
0.20000E+12泊松比为:
0.30000E+00厚度为:
0.50000E-02
节点1坐标为:
X=0.000Y=2.000
节点2坐标为:
X=0.000Y=1.000
节点3坐标为:
X=1.000Y=1.000
节点4坐标为:
X=0.000Y=0.000
节点5坐标为:
X=1.000Y=0.000
节点6坐标为:
X=2.000Y=0.000
单元号码为:
1组成单元的节点号码为:
123
单元号码为:
2组成单元的节点号码为:
245
单元号码为:
3组成单元的节点号码为:
253
单元号码为:
4组成单元的节点号码为:
563
节点1载荷为:
PY=********
节点1位移约束为:
PX=0.000
节点2位移约束为:
PX=0.000
节点4位移约束为:
PX=0.000
节点4位移约束为:
V=0.000
节点5位移约束为:
V=0.000
节点6位移约束为:
V=0.000
计算结果为:
最大半带宽为:
8
总数组大小为:
120
节点号1X方向位移UX=0.00000E+00Y方向位移UY=-0.32789E-03
节点号2X方向位移UX=0.00000E+00Y方向位移UY=-0.13795E-03
节点号3X方向位移UX=0.26479E-04Y方向位移UY=-0.36053E-04
节点号4X方向位移UX=0.00000E+00Y方向位移UY=0.00000E+00
节点号5X方向位移UX=0.56226E-04Y方向位移UY=0.00000E+00
节点号6X方向位移UX=0.67042E-04Y方向位移UY=0.00000E+00
单元1
EPTOX=0.26479E-04EPTOY=-0.18994E-03EPTOXY=0.10190E-03
SX=-0.67042E+07SY=-0.40000E+08SXY=0.78383E+07
S1=-0.49513E+07S2=-0.41753E+08TMAX=0.18401E+08
ANGEL=12.61°
单元2
EPTOX=0.56226E-04EPTOY=-0.13795E-03EPTOXY=0.00000E+00
SX=0.32618E+07SY=-0.26612E+08SXY=0.00000E+00
S1=0.32618E+07S2=-0.26612E+08TMAX=0.14937E+08
ANGEL=0.00°
单元3
EPTOX=0.26479E-04EPTOY=-0.36053E-04EPTOXY=0.72151E-04
SX=0.34424E+07SY=-0.61778E+07SXY=0.55500E+07
S1=0.59767E+07S2=-0.87121E+07TMAX=0.73444E+07
ANGEL=24.54°
单元4
EPTOX=0.10816E-04EPTOY=-0.36053E-04EPTOXY=-0.29747E-04
SX=-0.42330E+00SY=-0.72105E+07SXY=-0.22883E+07
S1=0.66487E+06S2=-0.78754E+07TMAX=0.42701E+07
ANGEL=-16.20°
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 力学 基础 有限元 程序