用伪代码方式写出有限元求解步骤

作业：一1用伪代码方式写出有限元求解步骤解：有限元分析流程图：形成K消除K奇异性求解Ka=PI形成K消除K奇异性求解Ka=PNGNG:结构的结点总数NE,MC,NX,NB.ND,EO,VO,TNE:结构单元总数MC:计算控制类型参数=0平面应力=1平面应变NX:作用荷载组数NB：给定位移的个数NB:给定位移的个数ND：结构刚度阵的半带宽EO:弹性模量VO：泊松比T：单元的厚度3.NWA.NWE,NWKNVP,NWDNWA：单元参数的输出控制参数NWE:单元刚度矩阵的输出控制参数NWK：结构刚度矩阵的输出控制参数NWP：荷载向量的输出控制参数NWD：结点位移的输出控制参数输出控制参数=1输出=0不输出4.IJM(3,NE):单元结点编码数组UM(1,I),UM(2,I)JJM(3,I),第I个三角形单元的节点号，按单元编号顺序填写。XY(2,NG)：结点坐标数组XY(1,I):第I个结点的坐标，XY(2,I):第I个结点的Y坐标按结点编号顺序填写。MB(2,NB).ZB(NB):给定位移约束的信息数组与值数组MB(1,I):第I个给定位移所在的结点号NB(2,1)=1：给定X方向位移=2：给定Y方向位移ZB(NB):给定位移值(以坐标正向为正)NF,NPNF：作用于结点上的集中荷载的个数NP：作用于均布侧压的单元边数若NF>0,填8MF(2,NF),ZF(NF):作用于结点上集中荷载的信息组与值数组MF(1,I):第I个集中荷载作用的结点号MF(2,1)=1:作用于x方向的集中力=2:作用y方向的集中力ZF(NF):作用的集中力值若NP>0,填9MP(2,NP),ZP(NP):作用于单元边上均布荷载的信息数组与值数组MF(LI):第I个均布荷载作用边的起始结点号MF(2,I):第I个均布荷载作用边的终止结点号，逆时针排列ZP(NP)：第I个均布荷载值输入数据格式，建立数据文件，文件名小于12个字符：NGNE,MC,NX,NB.ND,E,P,T3.NWA,NWE,NWKNWP,NWD4.IJM(3,NE)XY(2,NG)MB(2.NB),ZB(NB)NF,NP当NF>0测填8MF(2,NF),ZF(NF)当NP>0,则填9MP(2„NP),ZP(NP)若NX>1,即多组荷载情况，重复7〜9，若计算多个结构则重复1〜9。结束(2)平面问题主程序框图~——HNGNE.MC.NX.NB.NDEP.NVANWE.NWKNWRNWDL输入子程序INPUT伪代码DIMENSIONGUM(3,NE),XY(2,NG)・MB(2,NB),ZB(NB)READ(5,*)((IJM(I,L),I=1,3),L=I,NE)READ(5,*)((XY(I,J),I=1,2),J=LNG)READ(5,*)(MB(I,L),I=1,2),L=1,NB),(ZB(L),L=1,NB)WRITE(6,10)10FORMAT(/10X,4ELEMENTCODEBLOCKSUM(2,NE)'/)WRITE(6,20)((UM(MJ)M=l,3),1=1,NE)20FORMAT(IX,4(314,3X),3I4)WRITE(6,30)30FORMATE。」OX,4COORDINATESOFNODESXY(2,NG)/)WRITE(6,40)((XY(M,I),M=l,2),1=1,NG)40FORMATE(1X,6E12.5)WRITE(6,50)50FORMATE(/1OX.INRMATIONANDVALUESOFGIVENDISPLACEMENTS/20X,'MBL15X,^ZB\/)DO701=1.NBWRITE(6,60)(MB(K,I),K=l,2),ZB(I)60FORMATE(1OX,2I5,5X.E14.6)70CONGTINUERETURNEND计算单元参数子程序伪代码1、X(2,5)当前计算单元的结点坐标2、B(7)当前计算单元的单元参数3、BCA(7,NE)将B(7)按单元存储于该数组中DIMENSIONIJM(3,NE),XY(2,NG),BCA(7,NE)X(2,5),B(7)IF(NEWA・EQ』)WRITE(6,5)5FORMAT(/1OX,TARAMETERSOFELEMENTSBCA(7,NE)7)DO801=1.NEDO10K=l,3K1=UM(KJ)DO10J=l,2X(J,K)=XY(J・K1)10CONTINUEDO20J=l,2X(J,4)=X(J,1)X(J,5)=X(J,2)20CONTINUEDO30K=l,3B(K)=X(2,K+l).X(2,K+2)B(K+3)=X(LK+2).X(LK+1)30CONTINUEB(7)=(B(1)*B(5).B(4)*B(2))*0.5IF(NWA.GT.O)WRITE(6,40)I,B40FORMAT。X；NE=',13,/3X,7E10.4)IF(B(7).LE.0.0)GOTO60DO50J=l,7BCA(J,I)=B(J)50CONTINUEDOTO8060WIRTE(6,70)I,(UM(J,I),J=l,3)70FORMAT(/5X,ELEMENT\I5,5X,5AREAISNONPOSITIVE\5X.1JM\315)STOP80CONTINUERETURNEND计算单元刚度矩阵程序KE伪代码1、BCA(7,NE)单元参数组2、EK(6,6)单元刚度矩阵DIMENSIONB(7).BCA(7,NE),EK(6.6)DO101=1,7B(I)=BCA(I,IO)10CONTINUEA=A1/B(7)*TDO201=1,3DO20J=L311=2*1J1=2*JEK(Il.l,Jl-l)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3))EK(Il.l,Jl)=A*(V*B(I)*B(J+3)+A2*B(I+3*B(J))EK(Il,Jl-l)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3)EK(ILJl)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J))20CONTINEDO301=3,6DO30J=1,IEK(IJ)=EK(J,I)30CONTINUEIF(NEW.EQ.O)GOTO60WRITE(6,40)IO40FORMATflX,4EKNE=\I5)WRITE(6,50)EK50FORMAT(1X,6E1L4)60RETURNEND4、总刚集成子程序伪代码IO按单元循环的单元号EK(6,6)单元刚度矩阵IJM(3.NE)单元结点编码数组SK(NT,ND)按二维等宽存储的总刚度矩阵DOMENSIONIJ(3),SK(NT,ND),IJM(3.NE),EK(6,6)DO101=1,3IJ(I)=UM(I,IO)10CONTINUEDO201=1,3DO20J=l,3IF(IJ(I).IT(J))GOTO20M=2*IJ(I)-1N=2*IJ(J).U(I)+1MO=2*I-1NO=2*J-1SK(M,N)=SK(M,N)+EK(MO.NO)SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)SK(M+1,N+1)=SK(M+1,N+1TEK(MO+1,NO+1)IF(IJ(I).EQ・IJ(J))GOTO20SK(M+1,N.1)=SK(M+1,N.1)+EK(MO+1,NO)20CONTINUERETURNEND5、集中荷载子程序PF伪代码参数及数组说明：NF集中荷载的个数MF(2,NF)荷载信息组数ZF(NF)集中荷载数组F(NT)等效结点荷载列阵DEMENSIONMF(2,NF),ZF(NF),F(NT)READ(5,*)(MF(I,L),I=1,2),L=1,NF),(ZF(L),L.1,NF)WRITE(6,10)10FORMAT^1OX/CONCEENTRATEDLOADSMF(2,NF)ZF(NF)7)DO30L=LNFWRITE(6,20)(MF(I,L),I=1,2),ZF(L))20FORMAT(IX,2110,10X,E12.6)30CONTNUEDO401=1,NFN=2*MF(1,I).MF(2,I)F(N)=F(N)+ZF(I)40CONTINUEIF(NP.GT.O)GOTO70IF(NWP.EQ・0)GOTO70WRITE(6,50)50FORMAT^1OX.TOADVECTOR?)WRITE(6,60)F60FORMAT(1X,5E14.6)70RETURNEND6、均布侧压等效结点荷载子程序PP伪代码参数及数组说明NP作用侧压的单元边数MP(2,NP)作用侧压信息数组ZP(NP)均布侧压值数组XY(2,NG)结点坐标数组F(NT)等效结点荷载阵DIMENSIONMP(2,NP),ZP(NP),XY(2,NG),F(NT)READ(5,*)((MF(I,L),I=I,2),L=1,NF),(ZF(L),L=1,NF)WRITE(6,10)10FORMAT^1OX.^UNIFORMPRESSUREMP(2,NP)ZP(NP)7)DO30L=l,NPWRITE(6,20)(MP(I,L),I=1,2),ZP(L)20FORMAT(IX,2110,10X.E12.6)30CONTINUEDO401=1,NPNl=MP(bI)N2=MP(2,I)PX=XY(2,N1).XY(2,N2)PY=XY(1,N2).XY(LNl)PX=5*ZP(I)*PXPY=5*ZP(I)*PYF(2*N1.1)=F(2*N1.1)+PXF(2*N1)=F(2*N1)+PYF(2*N2-1)=F(2*N2.1)+PXF(2*N2)=F(2*N2)+PY40CONTINUEIF(NWP.EQ.O)GOTO70WRITE(6,50)50FORMAT(/,10X,IOADVECTOR?)WRITE(6,60)F60FORMAT(1X,5E14.6)70RETURNEND7、定位移子程序DBC伪代码参数说明：A(NT,ND)系数矩阵B(NT)等效结点荷载矩阵NB给定位移个数MB(2,NB)给定位移信息数组ZB(NB)给定位移数组DIMENSIONMB(2,NB),ZB(NB).A(NT,ND),B(NT)DO601=1,NBN=2*MB(1,I).MB(2,I)Z=ZB(I)IF(ABS(Z).LT.1E-1O)GOTO20IF(NX.NE.NXl)GOTO10A(N,1)=A(N,1)*1E+1510B(N)=A((N,1)*ZGOTO6020IF(NX.NE.NX1)GOTO50A(N,l)=1.0DO30J=2,NDA(NJ)=0.030CONTINUEDO40K=2,NDIF(N.LT.K)GOTO50M=N-K+]A(M,K)=0.040CONTINUE50B(N)=0.060CONTINUERETURNEND单元应力子程序stress伪代码参数说明：1.A1,A2,V2.IJNI(3,NE)单元结点编码数组BCA(7,NE)单元参数数组F(NT)结点位移DIMENSIONIJM(3,NE),BC(7,NE).F(NT),B(7),R(6)WRITE(6,5)5FORMAT^1OX,ELEMENTSTRESSES?)WRITE(6,10)10FORMAT(2X；NE\5X/S-X\8X；S-XY,,10X；S1\10X；S2\8X/AG7)DO601=1.NESl=0.S2=0.S3=0.DO20J=l,7B(J户BCA(J,I)20CONTINUEA=2*AUB(7)DO30J=l,3N=UM(J,I)*2R(2*J-1)=F(N.1)R(2*J)=F(N)30CONTINUEDO40J=l,3K=2*JS1=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))40CONTINUEP=.5*(S1+S2)Q=.5*(S1・S2)X1=P+SQRT(Q*Q+S3*S3)X2=2*P-X1CTA=0.IF(S3.GT.O)CTA=ATAN((X1-S1)/S3)WRITE(6,50)I,SI,S2,S3,XLX2,CTA50FORMAT(1X,I3,1X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,,2X,F8.4)60CONTNUERETURNEND2求集中力P=20作用点的位移，己知各竖直边边长为1,水平边长为2,E=2,u=0oft?：(i)先求单元刚度矩阵［KY<1000-10)(200、[B]=l/200020-2[D]=020k0120-2-1；ob<1000-10、[S]=[D][B]=00020-200.510-1-0.5r1000-10、00.510-1-0.5[K]e=[BY[D][B]tA=[BY[S]tA=^20120-2-100040-4-1-1-2031(