-
QQ截图20200413155402.png (174.36 KB, 下载次数: 256)
QQ截图20200413155402.png (174.36 KB, 下载次数: 256)
Subroutine SDB(NEL)
!子程序SDB的功能是形成几何矩阵[B]和应力矩阵【S】
Real :: A,B,C,D,E,F,G,H
Real::p,q,r,s,t,u,v,w,x,y,z
Integer::I,J
Dimension NODES(3)
Common /BB/LNODS(350,3),MATNO(350),&
COORDS(200,2),ID(400)
Common /CC/NFIX(100,3),AFIX(100,2),SL(400),AMAT(5,4)
Common /WORK/D(4,4),B(4,6),DB(4,6),ELSTIF(6,6),ICN(6)
Common /WORK2/DISP(6),STRESS(6),NODE(80),NELNO(80)
Common /COUNT1/MAXNEL,MAXNOD,MAXLOD,MAXFIX,NCNTRAL,&
NLOAD,MATER,NN
! AMAT(i,1)——第i种材料的杨氏弹性模量；
! N 一维存放的总体刚度元素个数
! NCNTRAL 求解问题的类别控制信息；NCNTRAL=1,2,3分别表示平面
! 应力，平面应变，轴对称问题。
COMMON /COUNT2/GAMA,AREA,T
DO I=1,NN
DO J=1,6
B(I,J)=0.0
End Do
End Do
!295~305 取出单元的三个节点号，存放在工作数组NODES中，对轴对称
!！问题求出单元形心的R0和Z0
DO I=1,3
NODES(I)=LNODS(NEL,I)
!！ 取出单元的三个节点号，存放在工作数组NODES中，对轴对称问题求出
! 单元形心的R0和Z0&
End do
IF(NN==4) THEN
X0=0.0
Y0=0.0
DO I=1,3
J=NODES(I)
X0=X0+COORDS(J,1)/3.0
Y0=Y0+COORDS(J,2)/3.0
End Do
End if
!306~319 求得坐标参数Ai,Bi,Ci（i=i,j,m）,并形成集合矩阵【B】
!如为轴对称问题，形成【B】矩阵中第二行的非零元素
DO I=1,3
IA=I+1-I/3*3
IB=3-I/2*(4-I)
IA=NODES(IA)
IB=NODES(IB)
J=2*I
B(1,J-1)=COORDS(IA,2)-COORDS(IB,2)
B(NN,J-1)=COORDS(IB,1)-COORDS(IA,1)
B(NN,J)=B(1,J-1)
B(NN-1,J)=B(NN,J-1)
IF (NN==4) THEN
AI=COORDS(IA,1)*COORDS(IB,2)-COORDS(IB,1)*COORDS(IA,2)
B(2,J-1)=(AI+B(4,J-1)*Y0)/X0+B(1,J-1)
END if
End do
!321求三角形单元的面积
AREA=(B(1,3)*B(NN,5)-B(1,5)*B(NN,3))/2.0
!322~324最终形成集合矩阵【B】
DO I=1,NN
DO J=1,6
B(I,J)=B(I,J)/(2.0*AREA)
End do
End do
!325~328 取出单元的材料号，并得到相应的E
I=MATNO(NEL)
YM=AMAT(I,1)
PR=AMAT(I,2)
T=AMAT(I,4)
!329~333 对于平面应变问题，对材料常数做相应的处理
IF(NCNTRAL==2)THEN
YM=YM/(1.0-PR*PR)
PR=PR/(1.0-PR)
T=1.0
END if
DO I=1,NN
DO J=1,NN
D(I,J)=0.0
End do
End do
!337~357 形成弹性矩阵【D】，对平面问题是3*3阶，对轴对称问题是4*4阶
If(NCNTRAL==3)THEN
T=6.283185*X0
EC2=YM/(1.0+PR)/(1.0-2.0*PR)
EC=EC2*(1.0-PR)/PR
D(1,3)=EC2
D(2,3)=EC2
D(3,1)=EC2
D(3,2)=EC2
D(1,2)=EC2
D(2,1)=EC2
D(1,1)=EC
D(2,2)=EC
D(3,3)=EC
D(4,4)=0.5*EC*(2.*PR)/(1-PR)
ELSE
EC=YM/(1.0-PR*PR)
D(1,1)=EC
D(2,1)=PR*EC
D(2,2)=EC
D(1,2)=PR*EC
D(3,3)=0.5*EC*(1.0-PR)
END If
!359~364 应用公式[S]=[D][B]形成应力矩阵存放于数组DB
DO I=1,NN
DO J=1,6
DB(I,J)=0.0
DO K=1,NN
DB(I,J)=DB(I,J)+D(I,K)*B(K,J)
End do
End do
End do
Return
END Subroutine
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |
