[Fortran] 纯文本查看 复制代码

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