平面问题有限元应力分析程序

第六章 有限元程序设计中的若干问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 基本步骤： ⅰ.结构离散化,输入或生成 结点信息-结点坐标 单元信息-单元结点编号 ⅱ.计算单元刚度矩阵,形成总体刚度矩阵,包括计算 ⅲ.形成结点载荷向量 ⅳ.引入约束条件 ⅴ.解线性方程组 ⅵ.求出结点位移 ⅶ.计算单元的应力并输出 §6-1  约束条件的处理 1． 对称性与反对称性 (1) 对称结构承受对称载荷作用时 (2) 对称结构承受反对称载荷作用 2. 约束位移的引入 主元置1法 主元赋大值 §6-2  总刚度矩阵的存贮法 1． 半带宽存贮法 x x x             x     O       x x x x             x x             x             x x x x x             x               x x                 x                   x x                         x x x x x x x x x x x x x x x x x X　 X　 X　 X　   X　   O X　           2． 一维压缩存贮法 考虑到总体刚度矩阵中各行的带宽并不相等，有时由于结构的几何形状的原因，使总体刚度矩阵某些行的带宽特别大。这种情况下如采用半带宽存贮法，就可能把许多零元素也包含了进去，这对节省计算机的存贮量是很不利的。 一维压缩存贮法是将总体刚度矩阵的下三角形中每一行从第一个非零元素开始按行将元素排成一序列，存放于一维数组 中。但是为了确定SK中的元素在[K]中的行列号，还需要将[K]中各行对角线的元素在伊维数组中的序号存放于另一辅助数组KD（N2）中（N2是总刚度矩阵的阶数）。现举例说明这一存贮法： 设有一系数阵 在一维数组SK（13）中依次存放的是 而辅助数组KD（6）中存放的是 KD（6）其实就是[K]中对角元素在一维数组SK（13）中的地址。 将一结构离散化后，对结点进行编号，就能依据单元号确定出总刚度矩阵[K]各行的带宽，由它依次累加就可得出其对角线元素一维存贮中的序号。 显然，形成了数组Kd，就确定了[K]中被存贮的元素分布情况以及SK和[K]中元素的对应关系，例如可求出[K]中第I行带宽为 也可确定出[K]中第I行左边第一个非零元素在[K]中的列号 此外，也能立即确定出单元刚度矩阵 中的某子矩阵 组集到一维数组存贮总刚度矩阵SK中的地址 ， ， 结构总刚度矩阵的以上两种存贮 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 ，一般应用于直接解法中。 附录A 平面问题有限元应力分析程序 本程序是用FORTRAN语言编写的教学使用程序，采用常应变三角形单元，用于解决弹性力学平面应力问题.它极易由学员扩充为求解平面应力和平面应变问题的通用应力分析程序。 程序中总刚度矩阵按一维压缩存贮，线性代数方程组用三角分解直接法求解。 为使教学程序易读易懂，计算时输入计算机的载荷是结点载荷，任何其他载荷都要事先换算成等效的结点载荷.程序所处理的约束仅是X或Y轴方向上的零位移约束。 计算结果除输出结点位移外，还输出单元形心处的σx，σy， τxy，主应力σ1,σ2及主方向角θ。 程序的结构便于修改和扩充，易于连接图形库扩充为包含前、后处理程序（网格自动生成及计算结果的图形显示）的更完整的程序系统. 1.  输入数据说明（以输入的先后为序，自由 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 ） （1）    NN，NE，KU，KV，KRX，KRY，EO，PO————6个整型数，2个实型数，其中 NN      结点总数（≤500）； NE      单元总数（≤700）； KU      x方向位移受约束的结点数（≤50）； KV      y方向位移受约束的结点数（≤50）； KRX    在x方向受结点载荷作用的结点数（≤60）； KRY    在y方向受结点载荷作用的结点数(≤60)； EO      材料的弹性模量； PO      材料的泊桑比； （2） X（NN）  ————    NN个实型数，为结点的x坐标. （3） Y（NN）  ————    NN个实型数，为结点的y坐标. （4） IJM（NE，3）————  3*NE个整型数，按行输入，为单元按逆时针向的结点编号 （5） JU (KU)    ————  KU个整型数，为x方向位移受约束的结点号. （6） JV（KV）  ————  KV个整型数，为y方向位移受约束的结点号. （7） NRX(KRX)  ————  KRX个整型数，为在x方向受结点载荷作用的结点号. （8） NRY(KRY)  ————  KRY个整型数，为在y方向受结点载荷作用的结点号. （9） RX(KRX)    ————  KRX个实型数，为x方向结点载荷的大小. （10）RY(KRY)    ————  KRY个实型数，为y方向结点载荷的大小. 2.    其他标识符说明 NF（=2*NN）          方程的总数； LK              总刚度矩阵下三角形一维存贮的总长； D1,D2,D3        弹性矩阵中的元素； B(3),C(3)        Bi,Ci(i=1～3)的工作单元； DEL              三角形单元面积的工作单元； U（NF）          开始存放结点载荷，求解后存放结点位移； SK(LK)          按一维存贮的总刚度矩阵； EK(21)          单元刚度矩阵下三角形元素按一维存贮的数组； BM（3,6）        用于存放2Δ〔B〕矩阵的元素； CM（3,6）        用于存放2Δ﹝D﹞〔B〕矩阵的元素； JD（NF）        总刚度矩阵下三角形对角线元素在一维存贮中的序号指示数组； JLL（6）        单元刚度矩阵元素和总刚度矩阵元素行列对应关系的工作数组； SGM（3）        存放单元形心处的应力σx，σy，τxy； SGM 1，SGM 2    存放单元形心处的主应力σ1,σ2； CET A            存放主应力的方向角. 3.    各子程序功能 SKDD            计算总刚度矩阵下三角形对角线元在一维存 贮中的序号 SHAPE（N）      用于计算第N个单元的有关常数. FEK              计算单元刚度矩阵，并存放于EK（21）中. SKKE            单元刚度矩阵向总刚度矩阵传送. FIXD              用于处理约束条件. SOLVE            解线性代数方程组. STRESS            用于计算单元 形心处的应力，主应力和主 向，并输出. 4.    出错信息的处理 本程序中若干数组是半动态数组,如求解问题的规模超过所设 数组下标的界限，则停机，等待处理 （1） 求解问题网格中的结点数大于500时，则停机.此时需将数组 X（500），Y（500），U(1000),JD(1000)下标的大小作相应的扩大. （2） 求解问题网格中的单元数大于700时，则停机.此时需将 JE(700,3)数组中的前一个下标的大小作相应的扩大. （3） 当求解问题的x方向或y方向的约束自由度数超过50时，则 停机.此时需将JU（50）或JV（50）下标的大小作相应的扩充. （4） 当求解问题的x方向或y方向的结点载荷数超过60时，则停 机.此时需将NRX（60）,RX(60)或NRY(60),RY（60）下标的大小作相应的扩充. （5） 当一维存贮的总刚度矩阵长度大于12000时，则停机.此时应 将数组SK(12000)下标的大小作相应的扩大. 5．  主程序框图 开 始 读入原始数据 调用SKDD,计算总刚度矩阵存贮信息 输出原始数据供校对，保存 形成结点载荷向量 调用SHAPE和FEK，计算单元刚度矩阵EKEK(21) 结 束 调用SHAPE和STRESS计算单元形心处应力、主应力和主方向，并输出 输出结点位移 调用SOLVE,解方程组 调用FIXD，处理约束 调用SKKE,将EK（21）送总刚度矩阵SK(LK) 6.源程序文本 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C  P.FEM---A STRESS ANALYSIS FINITE ELEMENT  C C          PROGRAM FOR PLANE PROBLEM          C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C            MAIN PROGRAM                  C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /MAT/ D1,D2,D3 COMMON /ELEM/ B(3),C(3),DEL COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) DIMENSION JU(50),JV(50),NRX(60),NRY(60),RX(60),RY(60) OPEN(7,FILE= ′DATA1′,STATUS= ′OLD′ ) OPEN(8,FILE= ′DATA2′,STATUS= ′NEW′ ) READ(7,* ) NN,NE,KU,KV,KRX,KRY,EO,PO IF(NN.LE.500) GOTO 10 WRITE( *,901) STOP 10  IF(NE.LE.700) GOTO 20 WRITE( *,902) STOP 20  IF(KU.LE.50) GOTO 30 WRITE( *,903) STOP 30  IF(KV.LE.50) GOTO 40 WRITE( *,904) STOP 40  IF(KRX.LE.60) GOTO 50 WRITE( *,905) STOP 50  IF(KRY.LE.60) GOTO 55 WRITE( *,906) STOP 55  READ(7，*) (X(I),I=1,NN) READ(7，*) (Y(I),I=1,NN) READ(7，*) ((JE(I,J),J=1,3),I=1，NE) READ(7，*) (JU(I),I=1,KU) READ(7，*) (JV(I),I=1,KV) READ(7，*) (NRX(I),I=1,KRX) READ(7，*) (NRY(I),I=1,KRY) READ(7，*) (RX(I),I=1,KRX) READ(7，*) (RY(I),I=1,KRY) NF=2 * NN CALL SKDD WRITE(8,1000) WRITE(8,1100) NN,NE,KU,KV,KRX,KRY,EO,PO WRITE(8,* )′X(NN)= ′ WRITE(8,1200)(X(I),I=1,NN) WRITE(8,* )′Y(NN)= ′ WRITE(8,1200)(Y(I),I=1,NN) WRITE(8,* )′JE(NE,3)= ′ WRITE(8,1300)((JE(I,J),J=1,3),I=1，NE) WRITE(8,* )′JU(KU)= ′ WRITE(8,1300)(JU(I),I=1,KU) WRITE(8,* )′JV(KV)= ′ WRITE(8,1300) (JV(I),I=1,KV) WRITE(8，*)ˊNRX(KRX)=ˊ WRITE(8,1300) (NRX(I),I=1,KRX) WRITE(8,*)ˊNRY(KRY)=ˊ WRITE(8,1300) (NRY(I),I=1,KRY) WRITE(8,*)ˊRX(KRX)=ˊ WRITE(8,1200) (RX(I),I=1，KRX) WRITE(8,*)ˊRY(KRY)=ˊ WRITE(8，1200) (RY(I),I=1，KRY) D1=EO/(1.-PO*PO) D2=D1*PO D3=D1*(1.-PO)/2 DO 60 I=1，LK 60  SK(I)=0. DO 100 I=1，NF 100  U(I)=0. DO 120 I=1，KRX K1=2*NRX(I)-1 120  U(K1)=RX(I) DO 140 I=1,KRY K1=2*NRY(I) 140  U(K1)=RY(I) DO 160 M=1,NE WRITE(*,*)ˊELE=ˊ,M CALL SHAPE(M) CALL FEK 160  CALL SKKE DO 180 I=1,KU 180  CALL FIXD(2*JU(I)-1) DO 200 I=1,KV 200 CALL FIXD(2,*JV(I)) CALL SOLVE WRITE(*,*)ˊSOLVE PASSEDˊ WRITE(8,3000) WRITE(8,4000) (I,U(2*I-1),U(2*I),I=1，NN) WRITE(8,5000) DO 240 M=1，NE CALL SHAPE(M) CALL STRESS(M) 240  CONTINUE STOP 901  FORMAT(5X,ˊNUMBER OF NODES EXCEEDSˊ *          ˊ ALLOWABLEˊ) 902  FORMAT(5X,ˊNUMBER OF ELEMENTˊ *          ˊEXCEEDA  ALLOWABLEˊ) 903  FORMAT(5X,ˊNUMBER OF FIXED FREEDOM AT X-ˊ, *      ˊDIRECTION EXCEEDS  ALLOWABLEˊ) 904  FORMAT(5X,ˊNUMBER OF FIXED FREEDOM AT Y-ˊ, *      ˊDIRECTION EXCEEDS  ALLOWABLEˊ) 905  FORMAT(5X,ˊNUMBER OF NODAL FORCES AT X-ˊ, *        ˊDIRECTION EXCEEDS  ALLOWABLEˊ) 906  FORMAT(5X,ˊNUMBER OF NODAL FORCES AT Y-ˊ, *        ˊDIRECTION EXCEEDS  ALLOWABLEˊ) 1000 FORMAT (//10X,ˊPRIMARY DATAˊ//,1X,ˊCONˊ, *        ˊSTANTS=ˊ) 1100 FORMAT (6I5,2F12.2) 1200 FORMAT (5F12.4) 1300 FORMAT (15I5) 3000 FORMAT (//10X,ˊNODAL DISPLACEMENTˊ/, *ˊNODEˊ,9X,ˊUˊ,13X,ˊVˊ,4X,ˊNODEˊ,9X, *          ˊUˊ,13X,ˊVˊ/ˊ) 4000 FORMAT (2(I5，2E14.4)) 5000 FORMAT (//2X,ˊELE.ˊ,2X,ˊELEMENT STESSˊ,/8X, *  ˊSGMXˊ,8X,ˊSGMYˊ,9X,ˊTXYˊ,8X,ˊSGM1ˊ, *    8X,ˊSGM2ˊ,10X,ˊCETAˊ/) END C C SUBROUTINE SKDD CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      CALCULATE ORDERS OF DIAGONAL          C C          ELEMENTS OF MATRIX K              C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK                COMMON /DATA1/ X(500),Y(500),U(10001),JE(700， 3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) JD(1)=1 JD(2)=3 DO 200 N=2，NN KK=0 DO 100 I=1，NE IO=JE(I,1) JO=JE(I,2) KO=JE(I,3) IF(IO.NE.N.AND.JO.NE.N.AND.KO.NE.N)GOTO 100 M1=N-IO M2=N-JO M3=N-KO IF(M1.GT.KK) KK=M1 IF(M2.GT.KK) KK=M2 IF(M3.GT.KK) KK=M3 100  CONTINUE KK=2 * KK JD(2*N-1)=JD(2*N-2)+KK+1 JD(2*N)=JD(2*N-1)+KK+2 200  CONTINUE LK=JD(NF) WRITE(*,*)ˊLK=ˊ,LK IF(LK.LE.12000) GOTO 210 WRITE(*,1000) STOP 1000 FORMAT(5X,ˊNUMBER OF ELEMENTS OF ONE-DIˊ, *ˊMENSIONAL GLOBAL MATRIX EXCEEDS ALLOWˊ, *            ˊABLEˊ) 210  RETURN END C C SUBROUTINE SHAPE(N) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      EVALUATE ELEMENT PARAMETERS          C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /DATA1/ X(500),Y(500),U(1000),JE(700，3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /ELEM/ B(3),C(3),DEL DIMENSION XO(6),YO(6) DO 100 I=1，3 K1=JE(N,I) K2=I+3 XO(I)=X(K1) XO(K2)=XO(I) YO(I)=Y(K1) YO(K2)=YO(I) K2=2*I-1 K3=K2+1 JLL(K2)=K1*2-1 100  JLL(K3)=K1*2 DO 200 I=1,3 K1=I+1 K2=I+2 B(I)=YO(K1)-YO(K2) 200  C(I)=XO(K2)-XO(K1) DEL=(B(1)*C(2)-B(2)*C(1))/2. RETURN END C C SUBROUTINE FEK CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C        CALCULATE STIFFNESS MATRIX          C C                OF ELEMENT                  C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  COMMON /STIF/SK(12000),EK(21),JD(1000),JLL(6) COMMON /ELEM/B(3),C(3),DEL COMMON /MAT/D1,D2,D3 DIMENSION BM(3,6),CM(3,6) DO 100 I=1,3 DO 100 J=1,6 100  BM(I,J)=0. DO 200 I=1,3 K2=I+I K1=K2-1 BM(1,K1)=B(I) BM(3,K2)=B(I) BM(2,K2)=C(I) BM(3,K1)=C(I) CM(1,K1)=D1*B(I) CM(1,K2)=D2*C(I) CM(2,K1)=D2*B(I) CM(2,K2)=D1*C(I) CM(3,K1)=D3*C(I) 200  CM(3,K2)=D3*B(I) DO 400 I=1,6 K1=I*(I-1)/2 DO 400 II=1,I Z=0. DO 300 JJ=1,3 300  Z=Z+BM(JJ,I)*CM(JJ.II) EK(K1+II)=Z/DEL/4 400  CONTINUE RETURN END C C SUBROUTINE SKKE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      ASSEMBLE GLOBAL STIFFNESS MATRIX      C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  COMMON /STIF/SK(12000),EK(21),JD(1000),JLL(6) DO 200 I=1,6 K1=JLL(I) KK=I*(I-1)/2 DO 200 J=1,I K2=JLL(J) IF(K2.LT.K1) GOTO 100 K=JD(K2)-K2+K1 GOTO 200 100  K=JD(K1)-K1+K2 200  SK(K)=SK(K)+EK(KK+J) RETURN END C C SUBROUTINE FIXD(K) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      READ BOUNDARY CONDITIONS              C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) L=JD(K) IF(K.EQ.1) GOTO 200 M=JD(K-1) IA=M+1 IB=L-1 DO 100 I=IA,IB 100  SK(I)=0. 200  IA=K+1 DO 300 I=IA,NF IF(JD(I)-JD(I-1).LT.I-K+1) GOTO 300 M=JD(I)-I+K SK(M)=0. 300  CONTINUE U(K)=0 RETURN END C C SUBROUTINE SOLVE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      SOLVE FOR SOLUTION OF EQUATIONS      C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) DO 300 I=1,NF IG=JD(I)-I IF(I.EQ.1) GOTO 10 MI=JD(I-1)-IG+1 GOTO 20 10  MI=1 20  DO 200 J=MI,I IP=IG+J JG=JD(J)-J IF(J.EQ.1) GOTO 30 MJ=JD(J-1)-JG+1 GOTO 40 30  MJ=1 40  IJ=MI IF(MJ.GT.MI) IJ=MJ J1=J-1 DO 100 K=IJ,J1 JK=JD(K) 100  SK(IP)=SK(IP)-SK(IG+K)*SK(JK)*SK(JG+K) IF(I.EQ.J) GOTO 210 JJ=JD(J) SK(IP)=SK(IP)/SK(JJ) U(I)=U(I)-SK(IP)*SK(JJ)*U(J) 200  CONTINUE 210  JI=JD(I) U(I)=U(I)/SK(JI) 300  CONTINUE DO 400 I=1，NF J=NF-I+1 IG=JD(J)-J IF(J.EQ.1) GOTO 310 MI=JD(J-1)-IG+1 GOTO 320 310  MI=1 320  JI=J-1 DO 400 K=MI,J1 U(K)=U(K)-SK(IG+K)*U(J) 400  CONTINUE RETURN END C C SUBROUTINE STRESS(N) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C      EVALUATE STRESSES OF ELEMENTS        C C            AND OUTPUT THEM                C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3)        COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /MAT/D1,D2，D3 COMMON /ELEM/B(3),C(3),DEL DIMENSION SGM(3),UE(6),CM(3，6) DO 100 I=1，6 K1=JLL(I) 100  UE(I)=U(K1) DO 150 I=1，3 K2=I+I K1=K2-1 CM(1,K1)=D1*B(I) CM(1,K2)=D2*C(I) CM(2,K1)=D2*B(I) CM(2,K2)=D1*C(I) CM(3,K1)=D3*C(I) CM(3,K2)=D3*B(I) 150  CONTINUE DO 200 I=1，3 SGM(I)=0. DO 200 J=1，6 200  SGM(I)=SGM(I)+CM(I,J)*UE(J)/DEL/2 S1=SGM(1) S2=SGM(2) S3=SGM(3) A1=SQRT(((S1-S2)/2.)**2+S3**2) A2=(S1+S2)/2 SM1=A2+A1 SM2=A2-A1 IF(ABS(S3).GT.1.E-4) GOTO 500 IF(S1.GT.S2) GOTO 400 CETA=90. GOTO 600 400  CETA=0. GOTO 600 500  CETA=ATAN((SM1-S1)/S3)*57.29578 600  WRITE(8,1000) N,SGM,SM1,SM2,CETA RETUR 1000 FORMAT(15,5F12.4,E14.6) END