鲜花( 0) 鸡蛋( 0)
|
本帖最后由 三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 |
|