设为首页收藏本站--- 驰名中外的国际土木工程技术交流平台!

东南西北人

 找回密码
 注册

QQ登录

只需一步,快速开始

总共8822条微博

动态微博

中国土木工程师手册(上中下)东南西北人英文资料走马观花500多专业手册、工程手册如何获取积分和金币?
精彩施工和土木工程技术视频东南西北人英汉对照资料汇总各版块精彩讨论贴汇总! 
查看: 1699|回复: 4

《有限单元法基本原理和数值方法》一书的源程序

[复制链接]
鲜花(0) 鸡蛋(0)
wolf008 发表于 2008-10-9 23:50:02 | 显示全部楼层 |阅读模式
本帖最后由 三T上人 于 2015-11-8 18:06 编辑 <br /><br />《有限单元法基本原理和数值方法》一书的源程序  
--------------------------------------------------------------------------------
C*********************************************************************
C *
C PL ---- PROGRAM OF PLANE PROBLEM 96.1 *
C *
C*********************************************************************
C
C-------- 输入数据顺序--------
C 1.NG          1整型
C NG 结构结点总数
C NG=0 则停止运行
C 2.NE,MC,NX,NB,ND,EO,VO,T    5整型3实型
C NE 结构单元总数
C MC 计算控制类型参数
C MC=0 平面应力
C =1 平面应变
C NX 作用载荷组数
C NB 给定位移个数
C ND 结构刚度矩阵的半带宽
C EO 弹性模量
C VO 泊松比
C T 单元(结构)的厚度
C 3.NWA NWE,NWK,NWP,NWD      5整型
C 输出控制参数
C =1 输出
C =0 不输出
C NWA 单元参数的输出控制参数
C NWE 单元刚度矩阵的输出控制参数
C NWK 结构刚度矩阵的输出控制参数
C NWP 载荷向量的输出控制参数
C NWD 结点位移的输出控制参数
C 4.IJM(3,NE) 单元结点编码数组    3×NE整型
C IJM(1,I); IJM(2,I); IJM(3,I)
C 第I个三角形单元的结点编号,按结点编号顺序填写
C 5.XY(2,NG) 结构结点坐标数组    2×NG实型
C 6.MB(2,N,ZB(N      2×NB整型,NB实型
C MB(1,I)---第I个给定位移所在的结点号
C MB(2,I)=1--给定X方向位移
C =0--给定Y方向位移
C ZB(I)----给定位移值(以坐标正向为正)
C 7.NF,NP        2整型
C NF-----作用在结点上的集中载荷(坐标方向)的个数
C NP-----作用均布侧压的单元边数
C 若 NF>0 则填写
C 8.MF(2,NF),ZF(NF)      2×NF整型,NF实型
C MF(1,I)---第I个集中载荷所在的结点号
C MF(2,I)=1--给定X方向集中力
C =0--给定Y方向集中力
C ZF(I)-----作用的集中力值
C 若 NP>0 则填写
C 9.MP(2,NP),ZP  (NP)      2×NP整型,NP实型
C MP(1,I)----第I个载荷作用边的起始结点号
C MP(2,I)----第I个载荷作用边的起始结点号
C ZP(I)------第I个均布载荷值
C
C 若 NX>1 重复 7.-9. (NX-1) 次
C
C 最后 NG=0 表示数据结束
C
C-------- 输出数据顺序--------
C 1.IJM(3,NE) 单元结点编码数组    3×NE整型
C IJM(1,I); IJM(2,I); IJM(3,I)
C 第I个三角形单元的结点编号,按结点编号顺序填写
C 2.XY(2,NG) 结构结点坐标数组    2×NG实型
C
C 若NWA=1,则输出
C 3.I,B(7) 单元参数 NE行,1×NE整型,7×NE实型
C 每行结构为:'NE='+单元号+Bi+Bj+Bm+Ci+Cj+Cm+A
C
C 若NWE=1,则输出
C 4.IO,EK(6×6)单元刚度阵 NE行,1×NE整型,6×6×NE实型
C 每行结构为:'NE='+单元号+EK(单元刚度阵)
C
C 若NWK=1,则输出
C 5.SK(NT,ND) 结构刚度矩阵 NT×ND=2NG×ND实型
C
C 若NWD=1,则输出
C 6.I,B 结点位移数据    NG行,1×NG整型,2×NG实型
C 每行结构为:单元号+U+V
C
C 7.S1,S2,S3,X1,X2,CTA
C 单元应力数据    6×NE实型
C 分别代表σx,σy,τxy,σ1,σ2和主应力方向
C
C 若 NX>1 重复 6.-7. (NX-1) 次
C
C--------可调数组分配--------
C
C 实型数组 C(100000) 整型数组 IA(100000)
C C(1) XY(2,NG) IA(1) IJM(3,NE)
C C(N1) ZB(N IA(M1) MB(2,M
C C(N2) BCA(7,NE) IA(M2) MF(2,N
C C(N3) SK(NT,ND) IA(M3) MP(2,NP)
C C(N4) F(NT) IA(MEND) 下限
C C(N5) ZF(NF)
C C(N6) ZP(NP)
C C(NEND) 下限
C
C--------程序停止代码--------
C  0 正常停止
C 111 数组C越界
C 222 数组C/IA越界
C 333 单元面积非正
C 444 结构刚度矩阵主元非正
C----------------------------------------------------------------------
C
C 主程序
C
       DIMENSION C(500000),IA(50000),EK(36)
       CHARACTER*12 IN,OUT
C  IN和OUT为输入文件和输出文件的文件名
       WRITE(*,*)' '
       WRITE(*,*)' PLEASE INPUT THE INPUT-FILE NAME (A<12)'
       WRITE(*,*)' '
       READ(*,5) IN
C  输入输入文件的文件名
       WRITE(*,*)' '
       WRITE(*,*)' PLEASE INPUT THE OUTPUT-FILE NAME (A<12)'
       WRITE(*,*)' '
       READ(*,5) OUT
C  输入输出文件的文件名
  5 FORMAT(A12)
       OPEN(5,FILE=IN, STATUS='OLD')
       OPEN(6,FILE=OUT,STATUS='UNKNOWN')
C  打开对应的输入和输出文件
10 READ(5,*) NG
       IF(NG.EQ.0) STOP
C  输入结构结点数;如果结点数为0则停止运行
       READ(5,*)NE,MC,NX,NB,ND,EO,VO,T
C  按顺序输入结构单元数,问题类型参数,载荷组数,给定位移个数
C 结构刚度阵的半带宽,弹性模量,泊松比和结构厚度
       READ(5,*)NWA,NWE,NWK,NWP,NWD
C  按顺序输入各输出控制参数
       NT=2*NG
C  确定总刚度矩阵阶数NT
C
C 计算变界数组的下限
C
       M1=3*NE+1
       M2=M1+2*NB
       N1=2*NG+1
       N2=N1+NB
       N3=7*NE+N2
       N4=N3+NT*ND
       N5=N4+NT
C 得到各变界数组在一维大数组中的起始元素编号
C
C 检验实型数组C的下限
C
       NEND=N5
       IF(NEND.LE.500000) GOTO 35
       WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY C(IN THE MIDDLE)!! ***'
       WRITE(*,30) NEND
30 FORMAT(/,'******** NEND=',I6,1X,'>80000 ********')
       STOP 111
C  若C下限超出500000,则给出错误信息并停止运行
C
C 数据输入
C
35 CALL INPUT(NE,NG,NB,IA(1),C(1),IA(M1),C(N1))
C  调用INPUT子程输入数据
       WRITE(*,40)
40 FORMAT(/10X,'##### INPUT PASSED #####')
C  显示提示信息
       IF(MC.EQ.0) GOTO 45
C  检验是否是平面应力问题
C
C 平面应变问题
C
       E=EO/(1.0-VO*VO)
       V=VO/(1.0-VO)
C  平面应变问题时,先进行弹性常数替换
       GOTO 50
C
C 平面应力问题
C
45 E=EO
       V=VO
50 NX1=NX
       A1=E/(1.0-V*V)/4.0
       A2=0.5*(1.0-V)
C  初始化NX1,A1和A2 / NX1为剩余载荷的组数
C
C 计算单元参数
C
       CALL ABC(NE,NG,NWA,IA(1),C(1),C(N2))
       WRITE(*,55)
55 FORMAT(/10X,'##### ABC PASSED #####')
C  调用ABC子程计算单元参数并显示提示信息
C
C 集成结构刚度矩阵K
C
       DO 60 I=N3,N4
        C(I)=0.0
60 CONTINUE
C  初始化结构刚度矩阵SK
       DO 65 K=1,NE
C  遍历结构的所有单元
        IO=K
        CALL KE(IO,NE,NWE,T,A1,A2,V,EK(1),C(N2))
        CALL SUMK(IO,NE,ND,NT,IA(1),C(N3),EK(1))
C   调用KE子程计算出单元刚度阵并调用SUMK子程将其集成到结构刚度阵中
65 CONTINUE
       WRITE(*,70)
70 FORMAT(/10X,'##### SUMK PASSED #####')
C  显示提示信息
       CALL CHECK(NT,ND,NWK,C(N3))
C  调用CHECK子程检验结构刚度阵中的主元是否非正
       WRITE(*,75)
75 FORMAT(/10X,'##### CHECK PASSED #####')
C  显示提示信息
80 READ(5,*) NF,NP
C  输入集中载荷个数NF和均布载荷个数NP
C
C 再次计算变界数组的下限
C 并检验实型数组C和整型数组IA的下限
C
       M3=M2+2*NF
       N6=N5+NF
       NEND=N6+NP-1
       MEND=M3+2*NP-1
C  计算C和IA的下限
       NM=0
       IF(NEND.LE.500000) GOTO 85
       WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY C (AT THE END)!! ***'
       WRITE(*,30) NEND
       NM=1
85 IF(MEND.LE.50000) GOTO 95
       WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY IA (AT THE END)!! ***'
       WRITE(*,90) MEND
90 FORMAT(/,'******** MEND=',I6,1X,'>500 ********')
       STOP 222
95 IF(NM.EQ.1) STOP 222
C  检验两个数组的下限,若下限超出则给出错误信息并停止运行
C
C 集成结构结点载荷列阵P
C
       DO 100 I=N4,N5
        C(I)=0.0
100 CONTINUE
C  初始化结构结点载荷列阵F
       IF(NF.GT.0) CALL PF(NF,NP,NT,NWP,C(N4),IA(M2),C(N5))
       WRITE(*,105)
105 FORMAT(/10X,'##### PF PASSED #####')
C 若集中载荷个数>0,调用PF子程集成各结点集中力并显示提示信息
       IF(NP.GT.0) CALL PP(NP,NT,NG,NWP,C(1),C(N4),IA(M3),C(N6))
       WRITE(*,110)
110 FORMAT(/10X,'##### PP PASSED #####')
C 若均布载荷个数>0,调用PP子程集成各均布载荷并显示提示信息
C
C 引入给定位移
C
       CALL DBC(NT,ND,NB,NX,NX1,C(N3),C(N4),IA(M1),C(N1))
       WRITE(*,115)
115 FORMAT(/10X,'##### DBC PASSED #####')
C 调用DBC子程引入给定位移消除系数矩阵的奇异性,并显示提示信息
C
C 求解线性方程组KA=P
C
       CALL GAUSS(NT,ND,NWD,NX,NX1,C(N3),C(N4))
       WRITE(*,120)
120 FORMAT(/10X,'##### GAUSS PASSED #####')
C 调用GAUSS子程求解线性方程组并显示提示信息
C
C 计算单元应力
C
       CALL STRESS(NE,NT,A1,A2,V,IA(1),C(N2),C(N4))
       WRITE(*,125)
125 FORMAT(/10X,'##### STRESS PASSED #####')
C  调用STRESS子程输出各单元应力并显示提示信息
       NX1=NX1-1
       IF(NX1.GT.0) GOTO 80
C  剩余载荷组自减1;若还有载荷剩余则继续计算
       GOTO 10
C  重新输入
       END
C*********************************************************************
C 1 *
C 子过程名称: INPUT *
C 子过程功能: 按顺序输入并输出计算所需数据 *
C *
C*********************************************************************
       SUBROUTINE INPUT(NE,NG,NB,IJM,XY,MB,Z
C 形参说明
C 输入:
C NE 整型,结构单元总数
C NG 整型,结构结点总数
C NB 整型,给定位移的个数
C IJM(3,NE) 整型,单元结点编码数组
C XY(2,NG) 实型,结构结点坐标数组
C MB(2,N 整型,位移约束信息数组
C ZB(N 实型,位移约束数值数组
       DIMENSION IJM(3,NE),XY(2,NG),MB(2,N,ZB(N
C
       READ(5,*) ((IJM(I,L),I=1,3),L=1,NE)
C 输入单元结点编码数组
       READ(5,*) ((XY(I,J),I=1,2),J=1,NG)
C 输入结构结点坐标数组
       READ(5,*)((MB(I,L),I=1,2),L=1,N,(ZB(L),L=1,N
C 输入位移约束信息和数值数组
       WRITE(6,20)((IJM(M,I),M=1,3),I=1,NE)
20 FORMAT(1X,4(3I4,3X),3I4)
C 输出单元结点编码数组
       WRITE(6,40)((XY(M,I),M=1,2),I=1,NG)
40 FORMAT(1X,6E12.5)
C 输出结构结点坐标数组
       RETURN
       END
C*********************************************************************
C 2 *
C 子过程名称:ABC *
C 子过程功能:根据各单元的结点坐标, *
C 计算并输出所有单元的各参数 *
C *
C*********************************************************************
       SUBROUTINE ABC(NE,NG,NWA,IJM,XY,BCA)
C 形参说明
C 输入:
C NE 整型,结构单元总数
C NG 整型,结构结点总数
C NWA 整型,单元参数输出控制,0-不输出,1-输出
C IJM(3,NE) 整型,单元结点编码数组
C XY(2,NG) 实型,结构结点坐标数组
C 输出:
C BCA(7,NE) 实型,结构单元参数数组
C 变量说明
C X(2,5) 实型,当前计算单元的结点坐标数组
C  B(7) 实型,当前计算单元的单元参数数组
       DIMENSION IJM(3,NE),XY(2,NG),BCA(7,NE),X(2,5),B(7)
C
       IF(NWA.EQ.1) WRITE(6,5)
  5 FORMAT(/10X,'PARAMETERS OF ELEMENTS BCA(7,NE)'/)
C 若控制打开,则输出单元参数提示信息
       DO 80 I=1,NE
C 遍历所有单元
        DO 10 K=1,3
C   遍历单元内的3个结点
         K1=IJM(K,I)
C 取结点在结构中对应的结点号
         DO 10 J=1,2
         X(J,K)=XY(J,K1)
C 得到当前结点的坐标值
10 CONTINUE
        DO 20 J=1,2
         X(J,4)=X(J,1)
         X(J,5)=X(J,2)
20 CONTINUE
C 为编程方便,每单元多存两个结点坐标
        DO 30 K=1,3
         B(K)=X(2,K+1)-X(2,K+2)
C 计算Bm
         B(K+3)=X(1,K+2)-X(1,K+1)
C 计算Cm
30 CONTINUE
        B(7)=(B(1)*B(5)-B(4)*B(2))*0.5
C 计算单元面积A
        IF(NWA.GT.0) WRITE(6,40)I,B
40 FORMAT(1X,'NE=',I3,/3X,7E10.4)
C 若输出控制打开,则输出单元号和对应的参数
        IF(B(7).LE.0.0) GOTO 60
C 若当前单元面积为负,则出错
        DO 50 J=1,7
         BCA(J,I)=B(J)
50 CONTINUE
C 将当前单元参数数组中的值传送给输出数组
        GOTO 80
60 WRITE(6,70)I,(IJM(J,I),J=1,3)
70 FORMAT(/5X,'ELEMENT',I5,5X,'AREA IS NONPOSITIVE',5X,'IJM',3I5)
        STOP 333
C 显示出错信息并停止运行
80 CONTINUE
       RETURN
       END
C*********************************************************************
C 3 *
C 子过程名称:KE *
C 子过程功能:根据结构单元参数数组,计算出 *
C 指定单元的单元刚度矩阵 *
C *
C*********************************************************************
SUBROUTINE KE(IO,NE,NWE,T,A1,A2,V,EK,BCA)
C 形参说明
C 输入:
C IO 整型,计算单元编号
C NE 整型,结构单元总数
C NWE 整型,刚度矩阵输出控制,0-不输出,1-输出
C IJM(3,NE) 整型,单元结点编码数组
C XY(2,NG) 实型,结构结点坐标数组
C T 实型,单元厚度
C A1 实型,材料系数,A1=E/(4*(1-V**2))
C A2 实型,材料系数,A2=(1-V)/2
C V 实型,泊松比
C  BCA(7,NE) 实型,结构单元参数数组
C 输出:
C EK(6,6) 实型,单元刚度矩阵
       DIMENSION B(7),BCA(7,NE),EK(6,6)
       DO 10 I=1,7
        B(I)=BCA(I,IO)
10 CONTINUE
C 得到计算单元的参数
       A=A1/B(7)*T
       DO 20 I=1,3
        DO 20 J=I,3
C 将单元刚度矩阵分块成3×3个子矩阵
         I1=2*I
         J1=2*J
         EK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3))
         EK(I1-1,J1)=A*(V*B(I)*B(J+3)+A2*B(I+3)*B(J))
         EK(I1,J1-1)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3))
         EK(I1,J1)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J))
C   计算每个子矩阵各元素的值
20 CONTINUE
       DO 30 I=3,6
        DO 30 J=1,I
         EK(I,J)=EK(J,I)
30 CONTINUE
C 根据对称性得到左下角矩阵的值
       IF(NWE.EQ.0) GOTO 60
C 若输出控制关闭,则直接结束子过程
       WRITE(6,40) IO
40 FORMAT(/1X,'EK NE=',I5)
       WRITE(6,50)EK
50 FORMAT(1X,6E11.4)
C 输出单元号和对应的刚度矩阵
60 RETURN
       END
C*********************************************************************
C 4 *
C 子过程名称:SUMK *
C 子过程功能:将指定单元的单元刚度矩阵集成到结构刚度 *
C 矩阵中,结构刚度矩阵以二维等带宽方式存储 *
C *
C*********************************************************************
       SUBROUTINE SUMK(IO,NE,ND,NT,IJM,SK,EK)
C 形参说明
C 输入:
C IO 整型,计算单元编号
C NE 整型,结构单元总数
C ND 整型,结构刚度矩阵的半带宽
C NT 整型,结构刚度矩阵的阶数
C IJM(3,NE) 实型,单元结点编码数组
C EK(6,6) 实型,单元刚度矩阵
C 输出:
C SK(NT,ND) 实型,结构刚度矩阵
       DIMENSION IJ(3),SK(NT,ND),IJM(3,NE),EK(6,6)
C
DO 10 I=1,3
        IJ(I)=IJM(I,IO)
10 CONTINUE
C  取出计算单元的结点号
       DO 20 I=1,3
        DO 20 J=1,3
C 遍历单元刚度矩阵的3×3个子矩阵
         IF(IJ(I).GT.IJ(J)) GOTO 20
C   如果是下三角元素,则不储存
         M=2*IJ(I)-1
         N=2*(IJ(J)-IJ(I))+1
C   得到在结构刚度矩阵中等带宽储存的行列码
         MO=2*I-1
         NO=2*J-1
C   得到对应在单元刚度矩阵中的行列码
         SK(M,N)=SK(M,N)+EK(MO,NO)
         SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)
         SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1)
         IF(IJ(I).EQ.IJ(J)) GOTO 20
C   不存储主子块的下三角元素
         SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO)
C 将单元刚度矩阵的元素叠加到结构刚度矩阵中
20 CONTINUE
       RETURN
       END
C*********************************************************************
C 5 *
C 子过程名称:CHECK *
C 子过程功能:检验结构刚度矩阵的主元并输出结构刚度矩阵, *
C 如果主元非正,则输出错误信息并停止运行 *
C *
C*********************************************************************
       SUBROUTINE CHECK(NT,ND,NWK,SK)
C 形参说明
C 输入:
C NT 整型,结构刚度矩阵阶数
C ND 整型,结构刚度矩阵半带宽
C NWK 整型,结构刚度矩阵输出控制,0-不输出,1-输出
C SK(NT,ND) 实型,结构刚度矩阵
C 变量说明
C M 整型,错误主元个数
DIMENSION SK(NT,ND)
C
       IF(NWK.EQ.0) GOTO 30
       WRITE(6,20) ((SK(I,J),I=1,NT),J=1,ND)
20 FORMAT(1X,5E13.6)
C 若输出控制打开,则输出结构刚度矩阵
30 M=0
C  置错误个数为0
       DO 50 I=1,NT
C  遍历结构刚度矩阵各主元
        IF(SK(I,1).GT.1E-10) GOTO 50
        WRITE(6,40)I,SK(I,1)
40 FORMAT(/10X,'MAIN ELEMENT IS NONPOSITIVE NT=',I4,5X,E12.6)
        M=M+1
C   若主元非正,则输出错误信息,并使错误个数+1
50 CONTINUE
       IF(M.GT.0) GOTO 60
       GOTO 70
60 STOP 444
C  若错误个数>0则停止运行
70 RETURN
       END
C*********************************************************************
C 6 *
C 子过程名称F *
C 子过程功能:将结点集中力装入等效载荷列阵, *
C 同时输出等效载荷列阵 *
C *
C*********************************************************************
SUBROUTINE PF(NF,NP,NT,NWP,F,MF,ZF)
C 形参说明
C 输入:
C NF 整型,坐标方向上集中载荷的个数
C NP 整型,作用均布侧压的边数
C NT 整型,等效载荷列阵元素个数
C NWP 整型,载荷输出控制,0-不输出,1-输出
C MF(2,NF) 整型,作用于结点上集中载荷的信息数组
C MF(1,I):第I个载荷作用的结点号
C MF(2,I):第I个载荷作用的方向,0-Y向,1-X向
C ZF(NF) 实型,作用于结点上集中载荷的值
C 输出:
C F(NT) 实型,等效载荷列阵
       DIMENSION MF(2,NF),ZF(NF),F(NT)
C
       READ(5,*) ((MF(I,L),I=1,2),L=1,NF),(ZF(L),L=1,NF)
C 输入集中载荷信息和数值数组
       DO 40 I=1,NF
C  遍历所有集中载荷
        N=2*MF(1,I)-MF(2,I)
C 得到载荷在结构等效载荷列阵中的对应顺序
        F(N)=F(N)+ZF(I)
C   将集中载荷叠加到结构等效载荷列阵中
40 CONTINUE
RETURN
       END
C*********************************************************************
C 7 *
C 子过程名称P *
C 子过程功能:计算均布侧压的等效结点载荷并装入等效载荷列阵, *
C 同时输出等效载荷列阵 *
C *
C*********************************************************************
       SUBROUTINE PP(NP,NT,NG,NWP,XY,F,MP,ZP)
C 形参说明
C 输入:
C NP 整型,作用均布侧压的边数
C NT 整型,等效载荷列阵元素个数
C NG 整型,结构结点总数
C NWP 整型,载荷输出控制,0-不输出,1-输出
C MP(2,NF) 整型,作用于单元边上均布载荷的信息数组
C MF(1,I):第I个载荷作用边的起始结点号
C MF(2,I):第I个载荷作用边的终止结点号
C ZP(NF) 实型,作用于单元边上均布载荷的值
C 输出:
C F(NT) 实型,等效载荷列阵
       DIMENSION MP(2,NP),ZP(NP),XY(2,NG),F(NT)
C
       READ(5,*) ((MP(I,L),I=1,2),L=1,NP),(ZP(L),L=1,NP)
C 输入均布载荷信息和数值数组
       DO 40 I=1,NP
C  遍历所有均布载荷
        N1=MP(1,I)
C   得到载荷的起始结点号
        N2=MP(2,I)
C   得到载荷的终止结点号
        PX=XY(2,N1)-XY(2,N2)
        PY=XY(1,N2)-XY(1,N1)
        PX=.5*ZP(I)*PX
        PY=.5*ZP(I)*PY
C   得到等效结点载荷 PX=qt(Yi-Yj)/2,PY=qt(Xi-Xj)/2
        F(2*N1-1)=F(2*N1-1)+PX
        F(2*N1)=F(2*N1)+PY
        F(2*N2-1)=F(2*N2-1)+PX
        F(2*N2)=F(2*N2)+PY
C   将均布载荷的等效结点载荷叠加到结构等效载荷列阵中
40 CONTINUE
RETURN
       END
C*********************************************************************
C 8 *
C 子过程名称BC *
C 子过程功能:引入给定位移的边界条件,消除系数矩阵的奇异性 *
C *
C*********************************************************************
       SUBROUTINE DBC(NT,ND,NB,NX,NX1,A,B,MB,Z
C 形参说明
C 输入:
C NT 整型,系数矩阵阶数
C ND 整型,系数矩阵的半带宽
C  NB 整型,给定位移的个数
C  NX 整型,载荷的总组数
C  NX1 整型,载荷的剩余组数
C A(NT,ND) 实型,系数矩阵    (兼输出)
C B(NT) 实型,等效结点载荷列阵  (兼输出)
C MB(2,N 整型,给定位移的信息数组
C MB(1,I):第I个给定位移的结点号
C MB(2,I):第I个给定位移的方向,0-Y向,1-X向
C ZP(N 实型,给定位移的值
       DIMENSION MB(2,N,ZB(N,A(NT,ND),B(NT)
C
       DO 60 I=1,NB
C  遍历所有给定位移
        N=2*MB(1,I)-MB(2,I)
C   取要修改的方程的序数
        Z=ZB(I)
C   取对应位移的值
        IF(ABS(Z).LT.1E-10) GOTO 20
C   若位移为0,用对角元素改1法,否则用对角元素乘大数法
        IF(NX.NE.NX1) GOTO 10
C   若不是第1组载荷,则只在B中引入给定位移
        A(N,1)=A(N,1)*1E+15
10 B(N)=A(N,1)*Z
C   对角元素乘大数,并修改对应的等效结点载荷
        GOTO 60
C   以下为对角元素改1法
20 IF(NX.NE.NX1) GOTO 50
C   若不是第1组载荷,则只在B中引入给定位移
        A(N,1)=1.0
C   将对角元素置1
        DO 30 J=2,ND
         A(N,J)=0.0
30 CONTINUE
C   将对应行元素置0
        DO 40 K=2,ND
         IF(N.LT.K) GOTO 50
         M=N-K+1
         A(M,K)=0.0
40 CONTINUE
C   将对应列元素置0
50 B(N)=0.0
C   等效结点载荷置0
60 CONTINUE
       RETURN
       END
C*********************************************************************
C 9 *
C 子过程名称:GAUSS *
C 子过程功能:使用高斯消元法求解线性方程组 *
C *
C*********************************************************************
       SUBROUTINE GAUSS(NT,ND,NWD,NX,NX1,A,
C 形参说明
C 输入:
C NT 整型,结构刚度矩阵阶数
C ND 整型,结构刚度矩阵半带宽
C NWD 整型,载荷向量输出控制,0-不输出,1-输出
C  NX 整型,载荷的总组数
C  NX1 整型,载荷的剩余组数
C A(NT,ND) 实型,系数矩阵
C B(NT) 实型,等效结点载荷列阵  (兼输出)
C
       DIMENSION A(NT,ND),B(NT)
C
C
C 消去过程
C
       N=NT-1
       IF(NX.EQ.NX1) GO TO 10
       ASSIGN 50 TO M
       GO TO 20
10 ASSIGN 30 TO M
20 DO 60 K=1,N
        M1=ND-1
        IF((M1).GT.(NT-K)) M1=NT-K
        DO 60 L=1,M1
         C=A(K,L+1)/A(K,1)
         IF(ABS(C).LT.1E-18) GO TO 60
         GO TO M,(30,50)
30 M2=ND-L
         DO 40 J=1,M2
          A(K+L,J)=A(K+L,J)-C*A(K,J+L)
40 CONTINUE
50 B(K+L)=B(K+L)-C*B(K)
60 CONTINUE
C
C 回代过程
C
       B(NT)=B(NT)/A(NT,1)
       DO 80 K=1,N
        I=NT-K
        M1=ND
        C=B(I)
        IF((K+1).LT.(ND)) M1=K+1
        DO 70 J=2,M1
         L=I+J-1
         C=C-A(I,J)*B(L)
70 CONTINUE
       B(I)=C/A(I,1)
80 CONTINUE
C
C 输出求解结果
C
       IF(NWD.EQ.0)GOTO 150
       N=NT/2
       N11=N/2
       IF(FLOAT(N11)+.3-FLOAT(N)/2.0.GT.1E-7)N=N-1
       DO 140 I=1,N,2
        NT2=NT/4
        NT3=NT/2
        FL1=FLOAT(NT2)+.3-FLOAT(NT3)/2
        IF((FL1.LT.0).AND.(I.EQ.N)) GOTO 120
        J=2*I-1
        K=2*I
        I1=I+1
        J1=J+2
        K1=K+2
        WRITE(6,110) I,B(J),B(K),I1,B(J1),B(K1)
110 FORMAT(1X,2(I3,3X,E11.5,2X,E11.5,5X))
        GOTO 140
120 J=2*I-1
        K=2*I
        WRITE(6,130)I,B(J),B(K)
130 FORMAT(1X,I3,3X,E11.5,2X,E11.5)
140 CONTINUE
150 RETURN
       END
C*********************************************************************
C 10 *
C 子过程名称:STRESS *
C 子过程功能:根据结构各结点的位移,输出各单元的 *
C 单元应力,主应力和应力主方向 *
C *
C*********************************************************************
       SUBROUTINE STRESS(NE,NT,A1,A2,V,IJM,BCA,F)
C 形参说明
C 输入:
C NE 整型,结构单元总数
C NT 整型,结点位移列阵元素个数
C A1 实型,材料系数,A1=E/(4*(1-V**2))
C A2 实型,材料系数,A2=(1-V)/2
C V 实型,泊松比
C IJM(3,NE) 整型,单元结点编码数组
C  BCA(7,NE) 实型,结构单元参数数组
C F(NT) 实型,结构结点位移列阵
C 变量说明
C B(7) 实型,当前计算单元的单元参数数组
C  R(6) 实型,当前计算单元的结点位移数组
C  A 实型,A=E/(2*(1-V**2)*Ae)
C  S1,S2,S3 实型,当前计算单元的应力σx,σy,τxy
C  X1,X2 实型,当前计算单元的主应力σ1,σ2
C CTA   实型,当前计算单元的主应力方向
DIMENSION IJM(3,NE),BCA(7,NE),F(NT),B(7),R(6)
C
       WRITE(6,5)
  5 FORMAT(/,10X,' ',/)
       DO 60 I=1,NE
C  遍历所有单元
        S1=0.
        S2=0.
        S3=0.
C   当前单元的应力值初始化
        DO 20 J=1,7
         B(J)=BCA(J,I)
C   取当前单元的单元参数
20 CONTINUE
        A=2*A1/B(7)
C 得到E/(2*(1-V**2)*Ae)
        DO 30 J=1,3
C   遍历单元内的3个结点   
         N=IJM(J,I)*2
         R(2*J-1)=F(N-1)
         R(2*J)=F(N)
C 取对应结点的结点位移
30 CONTINUE
        DO 40 J=1,3
         K=2*J
         S1=S1+A*(B(J)*R(K-1)+V*B(J+3)*R(K))
         S2=S2+A*(V*B(J)*R(K-1)+B(J+3)*R(K))
         S3=S3+A*A2*(B(J+3)*R(K-1)+B(J)*R(K))
40 CONTINUE
C   计算当前单元的应力
C   σx=A*Σ(Bi*Ui+V*Ci*Vi),σx=A*Σ(V*Bi*Ui+Ci*Vi)
C τxy=A*Σ(A2*Ci*Ui+A2*Bi*Vi)
        P=.5*(S1+S2)
        Q=.5*(S1-S2)
        X1=P+SQRT(Q*Q+S3*S3)
        X2=2*P-X1
C   计算当前单元的主应力
        CTA=0.
        IF (S3.GT.0) CTA=ATAN((X1-S1)/S3)
C   计算当前单元的主应力方向
        WRITE(6,50) S1,S2,S3,X1,X2,CTA
50 FORMAT(1X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,2X,F8.4)
C   按顺序输出当前单元的σx,σy,τxy,σ1,σ2和主应力方向
60 CONTINUE
       RETURN
       END

本站强荐:185娱乐《城.足球《真_人.彩票齐全《手机可投《注任何游戏. 首次开户送10元.首存送58元.信誉绝对保证185.cc
鲜花(0) 鸡蛋(0)
chengshuren 发表于 2009-6-25 08:35:34 | 显示全部楼层
鲜花(0) 鸡蛋(0)
jasmine927 发表于 2010-3-10 15:43:51 | 显示全部楼层
鲜花(3) 鸡蛋(0)
nevada 发表于 2010-5-3 09:44:58 | 显示全部楼层
本帖最后由 三T上人 于 2015-11-8 18:06 编辑 <br /><br />能实现什么功能的不?

威尼斯人:wns185.com首存赠送58元м足球м真_人м各类彩票齐全м提现即时到账
您需要登录后才可以回帖 登录 | 注册

本版积分规则

QQ|关于我们|QQ即时充值|站点统计|手机版|小黑屋|百宝箱|留言|咨询|微信订阅|QQ189615688|东南西北人

GMT+8, 2024-3-28 20:15 , Processed in 0.086331 second(s), 34 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表