1、C 有限单元法基本原理和数值方法一书的源程序 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,NW
2、P,NWD 5整型 C 输出控制参数 C =1 输出 C =0 不输出 C NWA 单元参数的输出控制参数 C NWE 单元刚度矩阵的输出控制参数 C NWK 结构刚度矩阵的输出控制参数 C NWP 载荷向量的输出控制参数 C NWD 结点位移的输出控制参数 C 4.IJM(3,NE) 单元结点编码数组 3NE整型 C IJM(1,I); IJM(2,I); IJM(3,I) C 第I 个三角形单元的结点编号,按结点编号顺序填写 C 5.XY(2,NG) 结构结点坐标数组 2NG实型 C 6.MB(2,N,ZB(N 2NB整型,NB实型 C MB(1,I)-第I个给定位移所在的结点号 C MB
3、(2,I)=1-给定X方向位移 C =0-给定Y方向位移 C ZB(I)-给定位移值(以坐标正向为正) C 7.NF,NP 2整型 C NF-作用在结点上的集中载荷( 坐标方向)的个数 C NP-作用均布侧压的单元边数 C 若 NF0 则填写 C 8.MF(2,NF),ZF(NF) 2NF整型,NF实型 C MF(1,I)-第I个集中载荷所在的结点号 C MF(2,I)=1-给定X方向集中力 C =0-给定Y方向集中力 C ZF(I)-作用的集中力值 C 若 NP0 则填写 C 9.MP(2,NP),ZP (NP) 2NP整型,NP实型 C MP(1,I)-第I个载荷作用边的起始结点号 C M
4、P(2,I)-第I个载荷作用边的起始结点号 C ZP(I)-第I个均布载荷值 C C 若 NX1 重复 7.-9. (NX-1) 次 C C 最后 NG=0 表示数据结束 C C- 输出数据顺序- C 1.IJM(3,NE) 单元结点编码数组 3NE整型 C IJM(1,I); IJM(2,I); IJM(3,I) C 第I 个三角形单元的结点编号,按结点编号顺序填写 C 2.XY(2,NG) 结构结点坐标数组 2NG实型 C C 若NWA=1,则输出 C 3.I,B(7) 单元参数 NE行,1NE整型,7NE实型 C 每行结构为:NE=+ 单元号+Bi+Bj+Bm+Ci+Cj+Cm+A C
5、C 若NWE=1,则输出 C 4.IO,EK(66)单元刚度阵 NE行,1NE 整型,66NE实型 C 每行结构为:NE=+ 单元号+EK(单元刚度阵) C C 若NWK=1,则输出 C 5.SK(NT,ND) 结构刚度矩阵 NTND=2NGND实型 C C 若NWD=1,则输出 C 6.I,B 结点位移数据 NG行,1NG整型,2NG 实型 C 每行结构为:单元号+U+V C C 7.S1,S2,S3,X1,X2,CTA C 单元应力数据 6NE实型 C 分别代表x,y,xy,1,2 和主应力方向 C C 若 NX1 重复 6.-7. (NX-1) 次 C C-可调数组分配- C C 实型数
6、组 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 结构刚度矩阵主元非正
7、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 (A80000 *) 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,#
8、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
9、(/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,#
10、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) GOT
11、O 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
12、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子程集
13、成各均布载荷并显示提示信息 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 计算单元应力
14、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
15、 输入: 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=
16、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,B
17、CA) 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(/1
18、0X,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
19、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)
20、,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
21、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 将
22、单元刚度矩阵分块成33个子矩阵 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
23、(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 N
24、E 整型 ,结构单元总数 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 遍历单元刚度矩阵的33个子矩阵 IF(IJ(I).GT.IJ(J) GOTO 20 C 如果是下三角元
25、素,则不储存 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
26、 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 WRI
27、TE(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 6
28、0 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(
29、2,I):第 I个载荷作用的方向,0Y向,1X 向 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
30、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)
31、实型 ,作用于单元边上均布载荷的值 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 得到等效
32、结点载荷 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
33、 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个给定位移的方向,0Y向,1X 向 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 取要
34、修改的方程的序数 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
35、)=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 整
36、型,载荷向量输出控制,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).
37、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
38、 输出求解结果 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
39、,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 实型
40、 ,材料系数 ,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 实型 ,当前计算单元的主应力方向 DIMENSIO
41、N 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
42、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
