平面問題有限元應力分析程序

1、第六章 有限元程序設計中的若干問題基本步驟：.結構離散化,輸入或生成結點信息-結點坐標單元信息-單元結點編號.計算單元剛度矩陣,形成總體剛度矩陣,包括計算.形成結點載荷向量.引入約束條件.解線性方程組.求出結點位移.計算單元的應力并輸出§6-1 約束條件的處理1 對稱性與反對稱性(1) 對稱結構承受對稱載荷作用時(2) 對稱結構承受反對稱載荷作用2. 約束位移的引入主元置1法主元賦大值§6-2 總剛度矩陣的存貯法1 半帶寬存貯法xxxxOxxxxxx x xxxxxxxxxxxxxxxxxxxxxxxxxxxxXXXXXOX2 一維壓縮存貯法考慮到總體剛度矩陣中各行的帶寬并

2、不相等，有時由于結構的幾何形狀的原因，使總體剛度矩陣某些行的帶寬特別大。這種情況下如采用半帶寬存貯法，就可能把許多零元素也包含了進去，這對節(jié)省計算機的存貯量是很不利的。一維壓縮存貯法是將總體剛度矩陣的下三角形中每一行從第一個非零元素開始按行將元素排成一序列，存放于一維數組中。但是為了確定SK中的元素在K中的行列號，還需要將K中各行對角線的元素在伊維數組中的序號存放于另一輔助數組KD（N2）中（N2是總剛度矩陣的階數）?，F舉例說明這一存貯法： 設有一系數陣在一維數組SK（13）中依次存放的是而輔助數組KD（6）中存放的是KD（6）其實就是K中對角元素在一維數組SK（13）中的地址。 將一結構離散

3、化后，對結點進行編號，就能依據單元號確定出總剛度矩陣K各行的帶寬，由它依次累加就可得出其對角線元素一維存貯中的序號。顯然，形成了數組Kd，就確定了K中被存貯的元素分布情況以及SK和K中元素的對應關系，例如可求出K中第I行帶寬為也可確定出K中第I行左邊第一個非零元素在K中的列號此外，也能立即確定出單元剛度矩陣中的某子矩陣組集到一維數組存貯總剛度矩陣SK中的地址，結構總剛度矩陣的以上兩種存貯方法，一般應用于直接解法中。附錄A平面問題有限元應力分析程序本程序是用FORTRAN語言編寫的教學使用程序，采用常應變三角形單元，用于解決彈性力學平面應力問題.它極易由學員擴充為求解平面應力和平面應變問題的通用

4、應力分析程序。 程序中總剛度矩陣按一維壓縮存貯，線性代數方程組用三角分解直接法求解。 為使教學程序易讀易懂，計算時輸入計算機的載荷是結點載荷，任何其他載荷都要事先換算成等效的結點載荷.程序所處理的約束僅是X或Y軸方向上的零位移約束。 計算結果除輸出結點位移外，還輸出單元形心處的x，y，xy，主應力1,2及主方向角。 程序的結構便于修改和擴充，易于連接圖形庫擴充為包含前、后處理程序（網格自動生成及計算結果的圖形顯示）的更完整的程序系統(tǒng).1. 輸入數據說明（以輸入的先后為序，自由格式） （1） NN，NE，KU，KV，KRX，KRY，EO，PO6個整型數，2個實型數，其中NN 結點總數（500）；

5、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個整型數，為

6、在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=13)的工作單元；DEL 三角形單元面積的工作單元；U（NF） 開始存放結點載荷，求解后存放結點位移；SK(LK) 按一維存貯的總剛度矩陣；EK(21) 單元剛度矩陣下三角形元素按一維存貯的數組； B

7、M（3,6） 用于存放2B矩陣的元素； CM（3,6） 用于存放2DB矩陣的元素； 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 用于處理

8、約束條件.SOLVE 解線性代數方程組.STRESS 用于計算單元 形心處的應力，主應力和主向，并輸出.4. 出錯信息的處理本程序中若干數組是半動態(tài)數組,如求解問題的規(guī)模超過所設數組下標的界限，則停機，等待處理（1） 求解問題網格中的結點數大于500時，則停機.此時需將數組X（500），Y（500），U(1000),JD(1000)下標的大小作相應的擴大.（2） 求解問題網格中的單元數大于700時，則停機.此時需將JE(700,3)數組中的前一個下標的大小作相應的擴大.（3） 當求解問題的x方向或y方向的約束自由度數超過50時，則停機.此時需將JU（50）或JV（50）下標的大小作相應的擴充.

9、（4） 當求解問題的x方向或y方向的結點載荷數超過60時，則停機.此時需將NRX（60）,RX(60)或NRY(60),RY（60）下標的大小作相應的擴充.（5） 當一維存貯的總剛度矩陣長度大于12000時，則停機.此時應將數組SK(12000)下標的大小作相應的擴大.5 主程序框圖開 始讀入原始數據調用SKDD,計算總剛度矩陣存貯信息輸出原始數據供校對，保存形成結點載荷向量對單元進行循環(huán)調用SHAPE和FEK，計算單元剛度矩陣EKEK(21)對單元進行循環(huán)結 束調用SHAPE和STRESS計算單元形心處應力、主應力和主方向，并輸出輸出結點位移調用SOLVE,解方程組調用FIXD，處理約束調用

10、SKKE,將EK（21）送總剛度矩陣SK(LK)6.源程序文本CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC P.FEM-A STRESS ANALYSIS FINITE ELEMENT CC PROGRAM FOR PLANE PROBLEM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAIN PROGRAM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /MAT/ D

11、1,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)

12、 GOTO 10 WRITE( *,901) STOP10 IF(NE.LE.700) GOTO 20 WRITE( *,902) STOP20 IF(KU.LE.50) GOTO 30 WRITE( *,903) STOP30 IF(KV.LE.50) GOTO 40 WRITE( *,904) STOP40 IF(KRX.LE.60) GOTO 50 WRITE( *,905) STOP50 IF(KRY.LE.60) GOTO 55 WRITE( *,906) STOP55 READ(7，*) (X(I),I=1,NN) READ(7，*) (Y(I),I=1,NN) READ(7，*)

13、 (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)

14、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)

15、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，LK60 SK(I)=0. DO 100 I=1，NF100 U(I)=0. DO 120 I=1，KRX K1=2*NRX(I)-1120 U(K1)=RX(I) DO 140 I=1,KRY K1=2*NRY(I)140 U(K1)=RY(I) DO 160 M=1,NE WRITE(*,*)E

16、LE=,M CALL SHAPE(M) CALL FEK160 CALL SKKE DO 180 I=1,KU180 CALL FIXD(2*JU(I)-1) DO 200 I=1,KV200 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 STOP901 FORMAT(5X,NUMBER O

17、F 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 ALLOWA

18、BLE)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 F

19、ORMAT (/2X,ELE.,2X,ELEMENT STESS,/8X, * SGMX,8X,SGMY,9X,TXY,8X,SGM1, * 8X,SGM2,10X,CETA/) ENDCC SUBROUTINE SKDDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALCULATE ORDERS OF DIAGONAL CC ELEMENTS OF MATRIX K CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /DATA

20、1/ 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=M3100 CONTINU

21、E KK=2 * KK JD(2*N-1)=JD(2*N-2)+KK+1 JD(2*N)=JD(2*N-1)+KK+2200 CONTINUE LK=JD(NF) WRITE(*,*)LK=,LK IF(LK.LE.12000) GOTO 210 WRITE(*,1000) STOP1000 FORMAT(5X,NUMBER OF ELEMENTS OF ONE-DI, *MENSIONAL GLOBAL MATRIX EXCEEDS ALLOW, * ABLE)210 RETURN ENDCC SUBROUTINE SHAPE(N)CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

22、CCCCCCCCCCCCCCCCCC EVALUATE ELEMENT PARAMETERS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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

23、)=Y(K1) YO(K2)=YO(I) K2=2*I-1 K3=K2+1 JLL(K2)=K1*2-1100 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 ENDCC SUBROUTINE FEK CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALCULATE STIFFNESS MATRIX CC OF ELEMENT CCCCCCCCCCCCCCCCC

24、CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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,6100 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

25、)=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,3300 Z=Z+BM(JJ,I)*CM(JJ.II) EK(K1+II)=Z/DEL/4400 CONTINUE RETURN ENDCC SUBROUTINE SKKECCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ASSEMBLE GLOBAL STIFFNESS M

26、ATRIX CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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 200100 K=JD(K1)-K1+K2200 SK(K)=SK(K)+EK(KK+J) RETURN ENDCC SUBROUTINE FIXD(K)CCCCCCCCCCCCCCCCCCCCC

27、CCCCCCCCCCCCCCCCCCCCCCCCCCC READ BOUNDARY CONDITIONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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,IB100 SK(I)=0.

28、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 ENDCC SUBROUTINE SOLVECCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SOLVE FOR SOLUTION OF EQUATIONS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /CONT/ NN,NE,NF,LK COMMON /STIF/

29、 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 2010 MI=120 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=140 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 CONTINUE210 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