《有限元理論與應(yīng)用教案》遼寧工程技術(shù)大學(xué)YXY4三結(jié)點(diǎn)平面三角形單元的源程序設(shè)計

1、 用有限元法計算結(jié)構(gòu)的強(qiáng)度和剛度的 一個顯著特點(diǎn)，就是需要應(yīng)用電子計算機(jī) 和數(shù)學(xué)模型設(shè)計程序框圖，并編制計算機(jī) 程序。下面以簡單的平面元為例說明程序 框圖設(shè)計和源程序的編制方法。 進(jìn)行計算，這就需要根據(jù)一定的力學(xué)模型 該程序的適用范圍是：該程序的適用范圍是： A. 問題類型。（L*M=0是平面應(yīng)力問題， L*M=1是平面應(yīng) 變問題） B. 載荷類型。（包括結(jié)點(diǎn)載荷和自重，其它非結(jié)點(diǎn)載荷 需事先換算成等效節(jié)點(diǎn)力） C. 支承方式（約束方式），（平面體至少具有三個獨(dú)立 支承，以保證彈性體幾何位置） D. 材料性質(zhì)：彈性體由一種材料組成，各單元的厚度相 同。 E. 計算方法：采用位移法。整體剛度矩陣

2、采用剛度繼承 法，按半帶寬存儲，解方程時采用帶消去法。 一、總框圖一、總框圖 根據(jù)三角元計算平面問題的全過程，總框圖設(shè)計如下： 它共包括8個子框圖： 子框圖1和2是輸入原始數(shù)據(jù) 子框圖36是中間運(yùn)算 子框圖7和8是解剛度方程求單元應(yīng)力并輸出位移和應(yīng)力 總框圖總框圖 開 始 輸入基本參數(shù) 輸入基本參數(shù) 形成整體剛度 形成載荷向量 引入約束條件 解方程輸出位移 求應(yīng)力輸出應(yīng)力 結(jié) 束 形成單元剛度 1 2 4 5 6 7 8 3 （子程序） （主程序） 二、子框圖二、子框圖 按照總框圖中8個子框圖的次序,把各子框圖的 詳細(xì)框圖設(shè)計如下： 子框圖子框圖1，輸入基本數(shù)據(jù)。 此框圖由主程序完成，主要是輸

3、入5個基本數(shù)據(jù)， 而且這5個參數(shù)是來控制后面的數(shù)據(jù)個數(shù)，所以必須 先讀入。同時計算出位移分量個數(shù)NJ2，半帶寬IDD =（相鄰結(jié)點(diǎn)碼的最大差值+1） * 2，子框圖如下： 子框圖子框圖2 開 始 輸入節(jié)點(diǎn)個數(shù)NJ，單元個數(shù)NE支桿個數(shù) NZ，結(jié)點(diǎn)載荷個數(shù)NPJ，半帶寬IDD 位移分量個數(shù)NJ2=NJX2 輸入其他數(shù)據(jù) 輸入問題類型代碼L*M 輸入彈性模量EO，泊松比AMU， 容重ALOU，單元厚度TE 輸入結(jié)點(diǎn)坐標(biāo)數(shù)組AJZ（NJ,2） 單元結(jié)點(diǎn)碼數(shù)組JM（NE,3） 輸入約束數(shù)組IZC（NZ） 是 輸入結(jié)點(diǎn)載荷數(shù)組PJ（NPJ,2） L * U = I ？ 否 EO EO/（1-AMU *

4、AMU） AMU AMU/（1-AMU） 從原始數(shù)據(jù)中找出單元IE的三個結(jié)點(diǎn)碼 求出cm，bm，cj. ，bj和面積AE 形成幾何矩陣B 形成彈性矩陣D 形成矩陣S，S=DB I A S K 1 ？ I A S K 2 ？ 形成矩陣AkE ：AKE STBAE.TE 是 是 否 否 子框圖子框圖3 它是子程序DUGD的框圖， 這是一個獨(dú)立的程序段， 對于任一單元IE來說，它 有三個功能： 1. 當(dāng)計算信息IASK=1時， 求出單元面積AE； 2. 當(dāng)計算信息IASK=2時， 形成應(yīng)力矩陣S； 3. 當(dāng)計算信息IASK=3時， 形成單元剛度矩陣SKE； . . . . . . . . . . .

5、 d n n . d n . . . . . . . . . . 矩陣AK矩陣AK* 子框圖子框圖4：形成整體剛度矩陣 這部分主要是根據(jù)整體剛度矩陣的組集方法，把每個 單元剛度矩陣的子塊推到整體剛度矩陣（總剛）的對應(yīng)行 列中去，形成整體剛度矩陣。然后取其上三角陣AKz， 按半帶存貯方式存貯成AKz*。 由于半帶寬貯后，元素的行碼不變，新的列碼等于原來 的列碼減行碼加1，即公式： 新列碼=原列碼-行碼+1 子框圖子框圖4： 調(diào)用子程序DUGD(IE,3) 行碼I由1到NJ2循環(huán) 列碼J由1到NJ2循環(huán) 單元碼IE由1到NE循環(huán) AKE中子塊行碼I由1到3循環(huán) 該子塊中的元素行碼II由1到2循環(huán)

6、單元行碼IH 2*（I-1）+II 半帶行碼IDH 2*（JM（E,I）-1）+II AKE中子塊列碼J由1到3循環(huán) 該子塊中的元素列碼JJ由1到2循環(huán) 單元列碼IL 2*（J-1）+JJ 整體列碼IZL 2*（JM（E,J）-1）+JJ 半帶列碼IDL IZL-IDH+1 AKZ*（IDH,IDL） AKZ*（IDH,IDL） + AKZ*（IH,IL） I D L 0 ？ AKZ*(I,J) 0 是 否 I由1到NJ2循環(huán) I由1到NPJ循環(huán) IE由1到NE循環(huán) P（2 * IE） P（2 * IE）+PE P（2 * JE） P（2 *JE）+PE P（2 * ME） P（2 * ME）

7、+PE PE ALOU *AE*TE/3 N P J 0 ? P（I） 0 J PJ(I,2)；P（J） PJ（I,1） A L O U 0 ？ 調(diào)用DUGD（IE，1） p 置零 把結(jié) 點(diǎn)載 荷送 到 p 把自 重引 起的 等效 結(jié)點(diǎn) 載荷 累加 到 p 否 否 是 是 子框圖子框圖5：形成結(jié)構(gòu)載荷向量 對于平面問題，根據(jù)所 劃分的結(jié)點(diǎn)總數(shù)NJ和每個結(jié) 點(diǎn)的位移分量數(shù)2，可確定 總載荷列陣p的階數(shù)為 NJ2=2 *NJ。 p的形成可分兩步：首 先把直接作用在結(jié)點(diǎn)的載荷 送到p中去；其次把單元自 重引起的等效結(jié)點(diǎn)載荷累加 到p中去。設(shè)單元E的容重 為 ，面積為AE，厚度為tE ，則y向等效結(jié)點(diǎn)

8、載荷為PE =- AE Te/3，其框圖如下： 100 0 0 1 00 0 0 子框圖子框圖6：引入約束條件 按前述約束條件引入方法，將結(jié)構(gòu)的支桿引入對矩陣 AKz*和結(jié)點(diǎn)載荷向量P進(jìn)行修改。 設(shè)第I號支桿對應(yīng)的位移碼為IZ，現(xiàn)在分別說明對 AKz、 AKZ*和P的修改方法和內(nèi)容。IZ列 1. 對矩陣對矩陣AKz的修改：的修改： IZ行 IZ行 2. 對矩陣對矩陣AKZ*的修改：的修改： 3. 對向量對向量P修改：修改： 第IZ行的對角線元素應(yīng)改 名為1，該行其它元素改為0， 第IZ列的非對角線元素也改為0 第IZ行中第一個元素應(yīng)改 名為1，該行其它元素改為0， 此外，以該行第一個元素為起

9、點(diǎn)，向右上方畫45斜線，則此 斜線上的元素也都改為0。 第IZ行元素應(yīng)改為0。 子框圖子框圖6： 支桿碼I由1到NZ循環(huán) 列碼J由2到IDD循環(huán) 列碼J由2到J0循環(huán) AKZ*（IZ-J+1，J） 0 I Z I D D ？ 支桿相應(yīng)的位移分量IZ IZC（I） AKZ*（IZ，I） 1 AKZ*（IZ，J） 0 J0 IDDJ0 IZ P(IZ) 0 是 否 nnnnn n b b x x aa aaa 11 1 11211 nkkjib a a bb nka a a aa k kk ik i k i kj kk ik ij k ij 2, 1,; 12 , 1; )( )( 其中 子框圖子

10、框圖7：解結(jié)構(gòu)剛度方程式 PAKZ * P 前面已經(jīng)得到位移法基本方程 ，現(xiàn)在采用帶 消元法求解 。并將解出的 存在 中。 1. 高斯消元法公式 設(shè)方程組為： 根據(jù)“線性代數(shù)”公式，其消元公式為 ： （A） ii n ij jijii nn n n axabx nni a b x / )( 12, 1 1 其中 niijb a a bb nka a a aa k kk ki i k i kj kk ki ij k ij , 1,; 12 , 1:; )( )( 其中 nkki, 2, 1 回代公式為： （B） 作如下修改（以保證在消元中只存在上三角元素） 1） 為了使元素aij限制在上三角部分

11、，應(yīng)規(guī)定列碼j大于 和等于行碼i，即j i； 2） 屬于下三角部分的元素aik應(yīng)換成akj。 則（A）式變?yōu)椋?（C） )( 0 00 000 000 000 000 00 000 0000 00000 D a aaa A ij kjkikk )( 0 00 000 000 000 000 00 000 0000 00000 E a aaa A ij kmklkn 半帶寬di列 j列 k行 i行 半帶寬d 中的i列 i列 d列 D列 i行 Kd1行 n行 k行 2. 帶消去法公式 為省存貯單元，在帶形對稱矩陣A中，只需存貯主對 角線以上的元素。即把A的上半部分斜帶中的元素存貯在 豎帶A*中。（

12、注意：其行碼不變，新列碼原列碼行 碼） ; ; )( )( k kl kL i k i kM kl kL iJ k iJ b a a bb a a a aa kiJMLdJ ikkidki nkl mm ; 1, 2 , 1 , 2, 1; 1 1,2 , 1; 1 kiJkjMaa kiLaa laa ijJaa kMkj kLki klkk iJij 1 1 1; 1 il J HiJii nnn axabx abx / )( / 1 1; 1 3 , 2; 12, 1 iJHinJ JJnni m m 由（），（）兩式新舊列碼的對應(yīng)關(guān)系為： 所以（）式變?yōu)椋ㄏ剑?（）式變?yōu)椋ɑ卮?/p>

13、式） （） （） （） ; ; * * k kl kL ii kM kl kL iJiJ p Ak Ak pp Ak Ak Ak AkAk kiJMLdJ kiLlikki dkink m m ; 1, 2 , 1 1; 1;, 2, 1 1; 1,2 , 1 il J IJiJii nnn AkAkP AkP / )( / 1 1 1 3 , 2 12, 1 inJ JJ nni m m PAKZ * pbxAka ; 在程序中線性方程組是： 所以方便框圖設(shè)計將符號作如下變換： 略去輪碼（k）行消元公式為： 回代公式為： 在編框圖時，把Jm換成J0，d換成IDD， 換成P 子框圖子框圖7：

14、IDDNJ2-I+1 ? 行碼由NJ2到1,步長-1循環(huán) 消元輪碼K由1到NJ2-1循環(huán) 行碼I由K+1到IM循環(huán) 列碼J由1到IDD-L+1循環(huán) IM NJ2IM K+IDD-1 NJ2K+IDD-1 ? L I-K+1 C AKz*(K,L)/AKZ*(K,1) AKz*(I,J) AKZ*(I,J)-C * AKZ*(K,M) M J+i-K P(I) P(I) - C * P(K) P(NJ2) P(NJ2)/AKz*(NJ2,1) J0 IDDJ0 NJ2-I+1 列碼J由2到J0循環(huán) 打印位移 P(I) P(I) P(I) /AKZ*(I,1) IH J+I-1, P(I) P(I

15、)- AKZ*(I,J) * P(IH) 是 否 否 是 e s 0 xy 2 1 2 MAYL MIYL AMIYL arctg y xy 29578.5790 AMIYL arctg y xy /180 子框圖子框圖8：求單元應(yīng)力 上面解方程中，得出結(jié)構(gòu)的各結(jié)點(diǎn)位移分量，并存于P 中，現(xiàn)在根據(jù)各結(jié)點(diǎn)位移分量求單元應(yīng)力的基本公式是： 求應(yīng)力矩陣S時，可調(diào)用子框圖3中的子程序DUGD，并 令其中的計算信息IASK=2。關(guān)于主應(yīng)力和主平面角的計算公 式如下： 平均應(yīng)力： 應(yīng)力圓半徑： 最大主應(yīng)力： 最小主應(yīng)力： 主平面角： 或 2/ )( yx PYL 22 2/ xyyx RYL RYLPYL

16、AMAYL RYLPYLAMIYL 子框圖子框圖8： SIG2=AMIYL ? 調(diào)用子程序DUGD(IE,2) 單元碼IE由1到NE循環(huán) 結(jié)點(diǎn)局部碼I由1到3循環(huán) 應(yīng)力分量碼I由1到3循環(huán) 坐標(biāo)方向碼J由1到2循環(huán) 單元行碼IH 2*（I-1）+J 半帶行碼IDH 2*（JM（IE,I）-1）+J WY(IH) P(IDH) 位移分量碼J由1到6循環(huán) YL(I) YL(I) +S(I,J) *WY(J) YL(I) 0 SIG1 YL(1); SIG2 YL(2) SIG3 YL(3) PYL （SIG1+SIG2）/2 RYL AMAYL PYL+RYL AMIYL PYL-RYL 2 2

17、32/21SIGSIGSIG AMIYLarctg yxy / 29578.5790 0 結(jié) 束 打印IE，SIG1，SIG2，SIG3， AMAYL，AMIYL， 。 否 是 已知平面問題，共分10個結(jié)點(diǎn)，9個單元，4個支桿， 在結(jié)點(diǎn)10作用有y向荷載10，按平面應(yīng)力問題計算，設(shè)E=1， =0.25， =0，厚度TE=1 0 222 21 43 7 9 10 10 2 2 2 x y 8 56 6 9 8 7 4 321 5 根據(jù)源程序，首先準(zhǔn) 備子框圖1和2中的輸入數(shù) 據(jù)。 輸入結(jié)點(diǎn)數(shù)NJ=10； 單元數(shù) NE=9； 支桿數(shù) NZ=4； 結(jié)點(diǎn)載荷數(shù)NPJ=1； 半帶寬： IDD=2*（4+

18、1）=10 例：例： 解：解： 輸入問題類型碼L*M=0； 彈性模量 EO=1； 泊松比 AMU=0.25； 容量 ALOU=0； 厚度 TE=1； 結(jié)點(diǎn)坐標(biāo)數(shù)組 JZ； 單元結(jié)點(diǎn)碼數(shù)組 JM； 支承數(shù)組 IZC； 結(jié)點(diǎn)載荷數(shù)組 PJ； 39 210 1098 896 976 865 562 673 743 632 521 63 44 42 25 23 21 06 04 02 00 JZJZ 22 14 2010 00 8 7 2 1 JZJZ 位移分量號 下面按程序輸入語句填寫數(shù)據(jù)輸入清單。（*數(shù)組按列輸入 ） 平面元輸入文件（PLANE.DAT） 10,9,4,1,10 0,1.0,.25

19、,0.,1.0 0.0,2.0,4.0,6.0,1.0,3.0,5.0,2.0,4.0,3.0,0.0,0.0,0.0,0.0,2.0, 2.0,2.0,4.0,4.0,6.0 1,2,3,3,2,5,6,6,8,2,3,4,7,6,6,7,9,9,5,6,7,6,5,8,9,8,10 1,2,7,8 0.0,10.0,0.0,20.0 平面元輸出文件（PLANE.OUT計算結(jié)果文件） 10 9 4 1 10 0 1.000 .250 .000 1.000 .0 2.0 4.0 6.0 1.0 3.0 5.0 2.0 4.0 3.0 .0 .0 .0 .0 2.0 2.0 2.0 4.0 4.

20、0 6.0 1 2 3 3 2 5 6 6 8 2 3 4 7 6 6 7 9 9 5 6 7 6 5 8 9 8 10 1 2 7 8 .000 10.000 .000 20.000 平面元輸出文件（PLANE.OUT計算結(jié)果文件） OUTPUTING OF THE NODE DISPLACEMENTS JD= U= V= 1 .0000000000 .0000000000 2 1.0941080000 17.7564800000 3 -1.0941100000 17.7564800000 4 .0000000000 .0000000000 5 -1.6411640000 15.678530

21、0000 6 -.0000005743 20.9598900000 7 1.6411640000 15.6785300000 8 .8205854000 25.3126800000 9 -.8205772000 25.3126700000 10 .0000084227 44.4729700000 平面元輸出文件（PLANE.OUT計算結(jié)果文件） OUTPUTING OF THE STRESSES AND ANGLES E= 1 SX= 1.4902300000 SY= 3.7727030000 TAU= 3.1136530000 S1= 5.9476780000 S2= -.684744800

22、0 CET= 55.0646100000 E= 2 SX= -.7399293000 SY= 1.4167200000 TAU= -.0000000142 S1= 1.4167200000 S2= -.7399292000 CET= 90.0000000000 E= 3 SX= 1.4902300000 SY= 3.7727010000 TAU= -3.1136530000 S1= 5.9476770000 S2= -.6847451000 CET= 124.9354000000 E= 4 SX= .9503176000 SY= .5189403000 TAU= -.6733334000 S1

23、= 1.4416650000 S2= .0275932600 CET= 143.8809000000 平面元輸出文件（PLANE.OUT計算結(jié)果文件） E= 5 SX= .9503175000 SY= .5189423000 TAU= .6733328000 S1= 1.4416650000 S2= .0275950400 CET= 36.1191300000 E= 6 SX= 1.8077500000 SY= 3.9486720000 TAU= 1.3845050000 S1= 4.6282800000 S2= 1.1281420000 CET= 63.8551100000 E= 7 SX=

24、 1.8077500000 SY= 3.9486700000 TAU= -1.3845030000 S1= 4.6282770000 S2= 1.1281430000 CET= 116.1449000000 E= 8 SX= -.2949150000 SY= 2.1026650000 TAU= -.0000001708 S1= 2.1026650000 S2= -.2949150000 CET= 90.0000100000 E= 9 SX= 1.6794190000 SY= 10.0000000000 TAU= -.0000002427 S1= 10.0000000000 S2= 1.6794

25、190000 CET= 90.0000000000 平面元源程序（PLANE.FOR） dimension AKZ(700,100),IZC(100),P(700),PJ(100,2),JM(400,3), 1AJZ(350,2) CHARACTER*20 INFILE,OUTFILE write(*,1000) 1000 format(/1x, INPUT FILE NAME = ) READ(*,1001) INFILE 1001 FORMAT(A20) WRITE(*,1002) 1002 FORMAT(1X, OUTPUT FILE NAME = ) READ(*,1001) OUTF

26、ILE OPEN(6,FILE=OUTFILE,STATUS=NEW) OPEN(5,FILE=INFILE) READ(5,*) NJ,NE,NZ,NPJ,IDD write(*,28) nj,ne,nz,npj,idd write(6,28) nj,ne,nz,npj,idd 28 format(1x,5i5/) NJ2=NJ*2 NPJ=NPJ+1 CALL YFA(NJ,NE,NZ,NPJ,IDD,NJ2,AKZ,IZC,P,PJ,JM,AJZ) STOP END 總剛度矩陣組約束數(shù)組 單元節(jié)點(diǎn)碼數(shù)組 總載荷列陣 節(jié)點(diǎn)載荷數(shù)組 節(jié)點(diǎn)坐標(biāo)數(shù)組 (方程階數(shù))自由度數(shù)NZ-支桿個數(shù) NJ N

27、E 以SAP5P輸入和輸出文件命名 的形式來對此程序輸入輸出數(shù) 據(jù)文件命名,即使用TURBO編 輯器生成數(shù)據(jù)文件。 輸入節(jié)點(diǎn)數(shù)，單元數(shù)，支桿數(shù) ，節(jié)點(diǎn)載荷數(shù)，半帶寬數(shù)。（ “*”號表示按固定格式輸入， 按自由格式輸入）。 打印輸入數(shù)據(jù)文件 位移分量個數(shù) 節(jié)點(diǎn)載荷個數(shù) 調(diào)用YFA子程序 主程序段結(jié)束 打印輸出數(shù)據(jù)文件 單元矩陣應(yīng)力S SUBROUTINE YFA(NJ,NE,NZ,NPJ,IDD,NJ2,AKZ,IZC,P,PJ,JM,AJZ) DIMENSION IZC(NZ),PJ(NPJ,2),S(3,6),AKE(6,6),AKZ(NJ2,IDD) 1,P(NJ2),WY(6),YL(3

28、),JM(NE,3),AJZ(NJ,2) READ(5,*) LXM,EO,AMU,ALOU,TE WRITE(*,29) LXM,EO,AMU,ALOU,TE WRITE(6,29) LXM,EO,AMU,ALOU,TE 29 FORMAT(1X,I5,2X,4F10.3/) READ(5,*) AJZ WRITE(*,38) (AJZ(I,J),I=1,NJ),J=1,2) WRITE(6,38) (AJZ(I,J),I=1,NJ),J=1,2) 38 FORMAT(2X,9F8.1/) READ(5,*) JM WRITE(*,31) (JM(I,J),I=1,NE),J=1,3) WR

29、ITE(6,31) (JM(I,J),I=1,NE),J=1,3) 31 FORMAT(1X,20I4/) READ(5,*) IZC WRITE(*,37) (IZC(I),I=1,NZ) WRITE(6,37) (IZC(I),I=1,NZ) 37 FORMAT(1X,10I5/) READ(5,*) PJ WRITE(*,36) (PJ(I,J),I=1,NPJ),J=1,2) WRITE(6,36) (PJ(I,J),I=1,NPJ),J=1,2) 36 FORMAT(1X,4F10.3/) YFA子程序 單元剛度矩陣總剛度矩陣 自由格式輸入五個參數(shù)，問 題類型代碼L*M，彈性模量 E

30、O，泊松比AMU，容重 ALOU，單元厚度TE。 打印輸 入的數(shù) 據(jù)文件 輸入節(jié)點(diǎn)坐標(biāo)數(shù)組 輸入單元節(jié)點(diǎn)碼數(shù)組 輸入約束數(shù)組 輸入節(jié)點(diǎn)載荷數(shù)組 單元行碼 半帶行碼 形成 整體 剛度 矩陣 IF(LXM.NE.1) GO TO 200 EO=EO/(1.0-AMU*AMU) AMU=AMU/(1.0-AMU) 200 CONTINUE DO 22 I=1,NJ2 DO 22 J=1,IDD 22 akz(i,j)=0.0 DO 33 IE=1,NE CALL DUGD(IE,3,TE,AMU,EO,JM,AJZ,NE,NJ,ME,JE,AE,S,IIE,AKE) DO 33 I=1,3 DO 3

31、3 II=1,2 IH=2*(I-1)+II IDH=2*(JM(IE,I)-1)+II DO 33 J=1,3 DO 33 JJ=1,2 L=2*(J-1)+JJ IZL=2*(JM(IE,J)-1)+JJ IDL=IZL-IDH+1 IF (IDL.LE.0) GO TO 33 AKZ(IDH,IDL)=AKZ(IDH,IDL)+AKE(IH,L) 33 CONTINUE DO 44 I=1,NJ2 44 P(I)=0.0 LM=0平面應(yīng)力問題；LM=1平面應(yīng)變問題。 在LM=0平面應(yīng)力問題：E=EO，=AMU； 在LM=1平面應(yīng)變問題中：E=E/（1-2）； = /（1- ） 對位移分量

32、循環(huán)； 對半帶寬數(shù)循環(huán)； 總剛度矩陣允零 對所有單元循環(huán) 調(diào)用子程序DUGD （IE3）形成 單元剛度矩陣 AKE中子塊行碼I由1到3循環(huán) 該子塊中元素行碼II由1到2循環(huán) 找對應(yīng)的行碼 AKE中子塊列碼J由1到3循環(huán) 該子塊中元素行碼JJ由1到2循環(huán) 單元列碼 整體列碼 半帶列碼 找對應(yīng)的列碼 如果IDL 0時轉(zhuǎn)到33句即只存在上半帶寬 剛度集成 對位移分量循環(huán)； 總載荷列陣允零，P允零。 形成結(jié)構(gòu)載荷向量 引入 約束 條件 節(jié)點(diǎn) 載荷 IF (NPJ.LE.1) GO TO 66 DO 55 I=2,NPJ IAJ=PJ(I,2) 55 P(IAJ)=PJ(I,1) 66 CONTINUE

33、 IF(ALOU.LE.0.0) GO TO 88 DO 77 IE=1,NE CALL DUGD(IE,1,TE,AMU,EO,JM,AJZ,NE,NJ,ME,JE,AE,S,IIE,AKE) PE=-ALOU*AE*TE/3.0 P(2*IIE)=P(2*IIE)+PE P(2*JE)=P(2*JE)+PE 77 P(2*ME)=P(2*ME)+PE 88 CONTINUE DO 99 I=1,NZ IZ=IZC(I) AKZ(IZ,1)=1.0 DO 23 J=2 ,IDD 23 AKZ(IZ,J)=0.0 IF(IZ.GT.IDD)GOTO 34 JO=IZ GOTO 45 如果NPJ

34、1，則轉(zhuǎn)到66句 I由NPJ到1循環(huán) 把節(jié)點(diǎn)載荷送到P中 ALOU 0 ？成立轉(zhuǎn)到88句 單元由1到NE循環(huán) 調(diào)用子程序DUGD（ IE.1）形成單元面積AE 則y向等效節(jié)點(diǎn)力為 3/ EEE TAP 把自重引起的等效結(jié)點(diǎn)載荷累加到P中 對每個支桿循環(huán) 對支桿對應(yīng)的位移分量碼IZ 修改IZ行第一列元素 列碼J由2到IDD循環(huán) 若IZIDD 修改IZ行其他元素 確定最大列碼 求解結(jié)構(gòu) 剛度方程 34 JO=IDD 45 CONTINUE DO 56 J=2,JO 56 AKZ(IZ-J+1,J)=0.0 P(IZ)=0.0 99 CONTINUE IB=NJ2-1 DO 67 K=1,IB IF

35、(NJ2.GT.(K+IDD-1)GOTO 78 IM=NJ2 GOTO 89 78 IM=K+Idd-1 89 CONTINUE IC=K+1 DO 67 I=IC,IM L=I-K+1 C=AKZ(K,L)/AKZ(K,1) ID=IDD-L+1 DO 98 J=1,ID M=J+I-K 98 AKZ(I,J)=AKZ(I,J)-C*AKZ(K,M) 67 P(I)=P(I)-C*P(K) P(NJ2)=P(NJ2)/AKZ(NJ2,1) L=NJ2-1 列碼J由2到J0循環(huán) 修改P矩陣 修改斜線上元素 消元輪碼k由1到NJ2-1循環(huán) 若NJ2(k+IDD-1)則轉(zhuǎn)向78句 求最大行碼IM

36、 ) 1 ,(),(KAKLKAKc zz ),(),(),(MKAKcJIAKJIAK zzz 行碼I由k+1到IM循環(huán) L=I-k+1 列碼J由1到IDD-L+1 M=J+1-K 修改右邊項 求最后一個位移分量 修 改 系 數(shù) DO 87 I=L,1,-1 行碼由NJ2-1到1，步長為“-1”循環(huán) 消 元 過 程 IF(IDD.GT.(NJ2-I+1) GOTO 76 JO=IDD GOTO 75 76 JO=NJ2-I+1 75 CONTINUE DO 65 J=2,JO IH=J+I-1 65 P(I)=P(I)-AKZ(I,J)*P(IH) 87 P(I)=P(I)/AKZ(I,1)

37、 WRITE(*,50) WRITE(6,50) 50 FORMAT(5X,35HOUTPUTING OF THE NODE DISPLACEMENTS/) WRITE(*,54) WRITE(6,54) DO 43 I=1,NJ D1=P(2*I-1) D2=P(2*I) WRITE(*,32) I,D1,D2 43 WRITE(6,32) I,D1,D2 循環(huán)打印 打印位移 )(IP 如果IDDNJ2-I+1成立，轉(zhuǎn)到76句 J0=IDD J0=NJ2-I+1 列碼由J=2到J0循環(huán) IH=J+I-1 P（I）=P（I）-AKz*（I，J） P（IH） P（I）=P（I）/AKZ（I，1）

38、 求最 大列 碼J0 求 其 余 位 移 分 量 回 代 過 程 坐標(biāo) 方向分量J由1到2循環(huán)； 54 FORMAT(3X,3HJD=,5X,2HU=,16X,2HV=) 32 FORMAT(3X,I5,2(3X,F15.10) WRITE(*,53) WRITE(6,53) 53 FORMAT(/5X,36HOUTPUTING OF THE STRESSES AND ANGLES/) DO 21 IE=1,NE CALL DUGD(IE,2,TE,AMU,EO,JM,AJZ,NE,NJ,ME,JE,AE,S,IIE,AKE) DO 19 I=1,3 DO 19 J=1,2 IH=2*(I-1

39、)+J IDH=2*(JM(IE,I)-1)+J WY(IH)=P(IDH) 19 CONTINUE DO 18 I=1,3 YL(I)=0.0 DO 18 J=1,6 YL(I)=YL(I)+S(I,J)*WY(J) 18 CONTINUE SIG1=YL(1) SIG2=YL(2) SIG3=YL(3) y xy x 對每個單元循環(huán) 調(diào)用DUGD子程序（IE， 2）形成應(yīng)力矩陣S 節(jié)點(diǎn)局部碼I由1到3循環(huán)； 單元行碼 整體行碼； 形成單元E的結(jié)點(diǎn)位移向量 應(yīng)力分量碼I由1到3循環(huán)； 應(yīng)位分量碼由1到6循環(huán)； 求應(yīng)力 S 2 2 2 xy yx RYL PYL=(SIG1+SIG2)/2.0

40、 ryl=sqrT(SIG1-SIG2)/2.0)*(SIG1-SIG2)/2.0)+SIG3*SIG3) AMAYL=PYL+RYL AMIYL=PYL-RYL IF(SIG2.EQ.AMIYL)GOTO 17 SIG4=SIG2-AMIYL CETA=90.0-57.29578*ATAN2(SIG3,SIG4) GOTO 16 17 CETA=0.0 16 CONTINUE WRITE(*,15)IE WRITE(6,15)IE 15 FORMAT(5X,2HE=,I5) WRITE(*,14)SIG1,SIG2,SIG3 WRITE(6,14)SIG1,SIG2,SIG3 14 FORM

41、AT(1X,3HSX=,F15.10,5X,3HSY=,F15.10,5X,4HTAU=,F15.10) WRITE(*,101)AMAYL,AMIYL,CETA WRITE(6,101)AMAYL,AMIYL,CETA 101 FORMAT(1X,3HS1=,F15.10,5X,3HS2=,F15.10,5X,4HCET=,F15.10) 21 CONTINUE RETURN END 求平均應(yīng)力2/ yx PYL 如果 , 0?2AMIYLSIG )(29578.5790 AMIYL arctg y xy 打印 , , x y xy 最大主應(yīng)力 最小主應(yīng)力 求主應(yīng)力: 求主平面角: 打印該單元號 打印,最大主應(yīng) 力AMAYL，最 小主應(yīng)力 AMIYL，和主 平面角； 此子程序段到此結(jié)束。 求BD AKES陣 單元節(jié)點(diǎn)碼數(shù)組 DK當(dāng)IASK=1時求AE 當(dāng)IASK=2時求s 當(dāng)IASK=3時求AKE DUGD子程序（IE）對一任 意單元IE來說，它有三個功能： C subprgram for forming the element stiffess SUBROUTINE DUGD(IE,