二、平面四邊形4結(jié)點(diǎn)等參有限單元法程序

1、罿芆蒈袂肁蒁莄袁芃芄螃袀羃薀蠆衿肅莂薅袈膇薈蒁袈芀莁蝿羇罿膃蚅羆肂荿薁羅膄膂蕆羄襖莇蒃羃肆芀螂羂膈蒅蚈羂芁羋薄羈羀蒄蒀肀肅芇螈聿膅蒂蚄肈芇芅薀肇肇蒀薆蚄腿莃蒂蚃芁蕿螁螞羈莁蚇蟻肅薇薃蝕膆莀葿螀羋膃螈蝿羈莈蚄螈膀膁蝕螇節(jié)蒆薆螆羂艿蒂螅肄蒅螀螄膇芇蚆螄艿蒃薂袃罿芆蒈袂肁蒁莄袁芃芄螃袀羃薀蠆衿肅莂薅袈膇薈蒁袈芀莁蝿羇罿膃蚅羆肂荿薁羅膄膂蕆羄襖莇蒃羃肆芀螂羂膈蒅蚈羂芁羋薄羈羀蒄蒀肀肅芇螈聿膅蒂蚄肈芇芅薀肇肇蒀薆蚄腿莃蒂蚃芁蕿螁螞羈莁蚇蟻肅薇薃蝕膆莀葿螀羋膃螈蝿羈莈蚄螈膀膁蝕螇節(jié)蒆薆螆羂艿蒂螅肄蒅螀螄膇芇蚆螄艿蒃薂袃罿芆蒈袂肁蒁莄袁芃芄螃袀羃薀蠆衿肅莂薅袈膇薈蒁袈芀莁蝿羇罿膃蚅羆肂荿薁羅膄膂蕆羄襖莇蒃羃

2、肆芀螂羂膈蒅蚈羂芁羋薄羈羀蒄蒀肀肅芇螈聿膅蒂蚄肈芇芅薀肇肇蒀薆蚄腿莃蒂蚃芁蕿螁螞羈莁蚇蟻肅薇薃蝕膆莀葿螀羋膃螈蝿羈莈蚄螈膀膁蝕螇節(jié)蒆薆螆羂艿蒂螅肄蒅螀螄膇芇蚆螄艿蒃薂袃罿芆蒈袂肁蒁莄袁芃芄螃袀羃薀蠆衿肅莂薅袈膇薈蒁袈芀莁蝿羇罿膃蚅羆肂荿薁 有限元程序設(shè)計(jì)平面四邊形4結(jié)點(diǎn)等參有限單元法程序設(shè)計(jì)1、程序功能及特點(diǎn)a.該程序采用四邊形4節(jié)點(diǎn)等參單元，能解決彈性力學(xué)的平面應(yīng)力應(yīng)變問題。b.前處理采用網(wǎng)格自動(dòng)劃分技術(shù)，自動(dòng)生成單元及結(jié)點(diǎn)信息。b.能計(jì)算受集中力、自重體力、分布面力和靜水壓力的作用。c.計(jì)算結(jié)點(diǎn)的位移和單元中心點(diǎn)的應(yīng)力分量及其主應(yīng)力。d.后處理采取整體應(yīng)力磨平求得各個(gè)結(jié)點(diǎn)的應(yīng)力分量。e.算

3、例計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果比較，并給出誤差分析。f.程序采用Visual Fortran 5.0編制而成。2、程序流程及圖框圖2- 程序流程圖圖2-子程序框圖其中，各子程序的主要功能為：INPUT輸入原始數(shù)據(jù)HUAFEN自動(dòng)網(wǎng)格劃分，形成COOR(2,NP),X,Y的坐標(biāo)值與單元信息CBAND形成主元素序號指示矩陣MA（*）SKO形成整體剛度矩陣KCONCR計(jì)算集中力引起的等效結(jié)點(diǎn)荷載RBODYR計(jì)算自重體力引起的等效結(jié)點(diǎn)荷載RFACER計(jì)算分布面力引起的等效結(jié)點(diǎn)荷載RDECOP支配方程LU三角分解FOBALU分解直接解法中的回代過程OUTDISP輸出結(jié)點(diǎn)位移分量STRESS計(jì)算單元應(yīng)力分

4、量OUTSTRE輸出單元應(yīng)力分量STIF計(jì)算單元?jiǎng)偠染仃嘑DNX計(jì)算形函數(shù)對整體坐標(biāo)的導(dǎo)數(shù)，1,2,3,4。FUN8計(jì)算形函數(shù)及雅可比矩陣JSFUN 應(yīng)力磨平-單元下的K=NCNSCN應(yīng)力磨平-單元下的右端項(xiàng)系數(shù)CNSUMSKN應(yīng)力磨平-單元下的右端項(xiàng)集成到總體的PSUMSTRS應(yīng)力磨平-單元下的集成到總體的KGAUSTRSS高斯消元求磨平后的應(yīng)力3、輸入數(shù)據(jù)及變量說明當(dāng)程序開始運(yùn)行時(shí)，按屏幕提示，鍵入數(shù)據(jù)文件的名字。在運(yùn)行程序之前，根據(jù)程序中INPUT需要的數(shù)據(jù)輸入建立一個(gè)存放原始數(shù)據(jù)的文件，這個(gè)文件的名字為INDAT.DAT。數(shù)據(jù)文件包括如下內(nèi)容：（1）總控信息。共一條，4個(gè)數(shù)據(jù)。L，B，

5、NNL，NNB，NM，NRL矩形體長度B-矩形體寬度NNL-L方向上劃分的結(jié)點(diǎn)數(shù)NNB-B方向上劃分的結(jié)點(diǎn)數(shù)NM單元材料類型數(shù)NR約束結(jié)點(diǎn)總數(shù)（2） 結(jié)點(diǎn)約束信息。共NR條，每條依次輸入：IP，IX，IYIP結(jié)點(diǎn)號IX、IY分別為IP結(jié)點(diǎn)在x，y方向的約束情況，如果約束填0，如果自由填1。（3）材料信息。共NM條，每條依次輸入：JJ,（AE（I，JJ），I=1，4）JJ材料類型號，（AE（I，JJ），I=1，4）該材料的材料參數(shù)，共4個(gè)參數(shù)，排列順序?yàn)椋簭椥阅Ａ俊⒉此杀?、容重、單元厚度。?） 荷載信息a. 荷載控制信息。共一條，3個(gè)數(shù)據(jù)NCP，IZNCP受集中力作用的結(jié)點(diǎn)數(shù)IZ面力批數(shù)b.

6、若NCP0，輸入IP，PX，PYIP結(jié)點(diǎn)號PX、PY分別為IP結(jié)點(diǎn)x，y方向的集中力分量。c. 若IZ0，輸入面力荷載信息，共IZ批，按批輸入：JS，NSE，（WG（I）I=1，4）JS面力批號NSE第JS批面力受到面力作用的單元個(gè)數(shù)，（WG（I），I1，4）該面力的特征參數(shù)共4個(gè)數(shù)據(jù)，第1個(gè)數(shù)為面力類型，填1表示受靜水壓力作用，填2表示受均布法向壓力作用；第2個(gè)數(shù)為水壓密度，如果是均布壓力情況，就填均布壓力的集度；第3個(gè)數(shù)為最高水位的y坐標(biāo)，如果是均布壓力情況，可以填任意數(shù)；第4個(gè)數(shù)為面力作用的單元面的面號，單元面號的規(guī)定見圖2-3。（IEW（M），M1，NSE），IEW（*）為受面力作用的

7、單元的單元號，共NSE個(gè)。 圖3-1 單元結(jié)點(diǎn)編碼與面號4、理論基礎(chǔ)和求解過程4.1、構(gòu)造插值函數(shù)：本有限元計(jì)算采取的是四邊形八結(jié)點(diǎn)等參元進(jìn)行插值計(jì)算的。直接調(diào)用教材115頁3.3.21的結(jié)果，寫出所有插值函數(shù)：; 4.2位移插值函數(shù)及應(yīng)變應(yīng)力求解：在有限元方法中單元的位移模式一般采用多項(xiàng)式作為近似函數(shù)，多項(xiàng)式的選取應(yīng)由低次到高次的完備多項(xiàng)式。位移模式的選取一般為：u。，為位移模式，為廣義坐標(biāo)向量。根據(jù)方程求解得出廣義坐標(biāo)，可將位移函數(shù)表示成結(jié)點(diǎn)位移的函數(shù)，即 ,，寫成矩陣的形式為： N稱為插值函數(shù)矩陣或形函數(shù)矩陣，為單元結(jié)點(diǎn)位移列陣。確定了單元位移后，可以很方便地利用幾何方程和物理方程求得單

8、元的應(yīng)變和應(yīng)力。單元應(yīng)變?yōu)椋築稱為應(yīng)變矩陣，L是平面問題的微分算子，其中：，根據(jù)物理方程可以求得單元應(yīng)力，其中，S稱為應(yīng)力矩陣，B是應(yīng)變矩陣。由于是平面應(yīng)力問題，E0和v0取為E和v，所以彈性矩陣。這部分內(nèi)容參考了教材第2.2節(jié)。4.3、等參元變換：為了將局部（自然）坐標(biāo)中幾何形狀規(guī)則的單元轉(zhuǎn)換成總體坐標(biāo)中幾何形狀扭曲的單元，以滿足對一般形狀求解域進(jìn)行離散化的需要，必須建立坐標(biāo)轉(zhuǎn)換：,最方便的方法是將上式中表示成插值函數(shù)的形式，, 在有限元分析中，為建立求解方程，需要進(jìn)行各個(gè)單元面積內(nèi)的積分，其一般形式為：而g中包含場函數(shù)對于總體坐標(biāo)x，y的導(dǎo)數(shù)。由于場函數(shù)是用自然坐標(biāo)表述的，又因?yàn)樵谧匀蛔鴺?biāo)

9、內(nèi)的積分限是規(guī)格化的，因此希望能在自然坐標(biāo)內(nèi)按規(guī)格化的數(shù)值積分方法進(jìn)行上述積分。為此需要建立兩個(gè)坐標(biāo)系內(nèi)導(dǎo)數(shù)之間的變換關(guān)系。設(shè), ，xi，yi，zi是結(jié)點(diǎn)在總體坐標(biāo)內(nèi)的坐標(biāo)值，Ni為形狀函數(shù)，實(shí)際上它也是用局部（自然）坐標(biāo)表示的插值函數(shù)。對于等參變換，結(jié)點(diǎn)數(shù)、結(jié)點(diǎn)號和插值函數(shù)都不變。函數(shù)Ni對的偏導(dǎo)數(shù)可以表示成： 集合成矩陣形式就是： 式中的J稱為雅克比矩陣，利用, ，J可以表示為自然坐標(biāo)的函數(shù)，即：,Ni對于x,y的偏導(dǎo)數(shù)可以用自然坐標(biāo)顯示地表示為:其中是的逆矩陣,它可以按照下式計(jì)算得到：,是J的伴隨矩陣，它的元素是J的元素的代數(shù)余子式。4.4、總體應(yīng)力磨平根據(jù)有限單元法第5章中的（5.3.

10、15），即 其另一種矩陣形式為：令并利用以下表達(dá)式，和以及就可以像單元元?jiǎng)偠染仃嚨娇傮w剛度矩陣一樣，求出，只不過，這里有Ke為12*12的矩陣，而“總剛”K為3*NP階的矩陣，NP為結(jié)點(diǎn)數(shù)。5、算例應(yīng)用本程序計(jì)算一個(gè)矩形土體受到均布荷載時(shí)的位移和應(yīng)力，如下圖所示，三面約束的土體尺寸為40m*10m，取一半為20m*10m，E=1.0104 kN/m2，在右端受均布荷載 kN/m2，不考慮自重體力。給定網(wǎng)格大小為，有限元網(wǎng)格自動(dòng)劃分如圖3-2所示，單元總數(shù)2000，節(jié)點(diǎn)總數(shù)2121。由于對稱性，右端約束為滑移支座，只限制x方向位移。雖然土體一般不為彈性，但是方法是一樣的，只是剛度矩陣改動(dòng)即可使用

11、彈塑性模型。圖5-1 算例模型圖5-2 ANSYS計(jì)算模型圖實(shí)際計(jì)算結(jié)果與ANSYS分析結(jié)果的比較（詳細(xì)計(jì)算結(jié)果見后面）：（1） 位移比較比較項(xiàng)目X方向最大位移Y方向最大位移本程序計(jì)算結(jié)果0.0863mm5.890mmANSYS分析結(jié)果0.1200mm6.000mm誤差28.0%1.83%（2）應(yīng)力比較比較項(xiàng)目最大主應(yīng)力本程序計(jì)算結(jié)果9.996ANSYS分析結(jié)果9.940誤差0.56%誤差分析：本程序計(jì)算得到的Y方向最大位移和ansys計(jì)算結(jié)果很相近，由于x方向上的位移不占主要部分，因此，其誤差雖然有些大，但對總體位移影響不大。 由于網(wǎng)格較密，因此，計(jì)算得到的單元應(yīng)力值（未磨平）與ansys結(jié)

12、果相近，誤差小于1%，不必要應(yīng)力磨平也可以達(dá)到較好的精度。 本程序可以進(jìn)行總體應(yīng)力磨平，但是由于網(wǎng)格數(shù)較多，因此，磨平矩陣階數(shù)較大，可能計(jì)算時(shí)間也很長，甚至無法計(jì)算。（3）實(shí)際計(jì)算結(jié)果MATLAB圖與ANSYS圖比較：(a)(b)圖5-3 X方向應(yīng)力圖，（a)ANSYS （b) MATLAB 圖5-4 主應(yīng)力計(jì)算結(jié)果的MATLAB 圖圖5-5 主應(yīng)力計(jì)算結(jié)果的ANSYS 圖圖5-6 剪應(yīng)力計(jì)算結(jié)果的MATLAB 圖圖5-7 剪應(yīng)力計(jì)算結(jié)果的ANSYS 圖圖5-8 ANSYS總位移圖圖5-9 ANSYS X方向位移圖圖5-10 ANSYS Y方向位移圖圖5-11 ANSYS Y方向應(yīng)力圖附錄：（

13、1）輸入文件數(shù)據(jù)：平面四邊形四結(jié)點(diǎn)單元輸入數(shù)據(jù)L B NNL NNB NM NR20.0 10.0 101 21 1 141*約束信息數(shù)據(jù)結(jié)點(diǎn)號 X向約束信息 Y向約束信息1002003004005006007008009001000110012001300140015001600170018001900200021002200430064008500106001270014800169001900021100232002530027400295003160033700358003790040000421004420046300484005050052600547005680058900610

14、006310065200673006940071500736007570077800799008200084100862008830090400925009460096700988001009001030001051001072001093001114001135001156001177001198001219001240001261001282001303001324001345001366001387001408001429001450001471001492001513001534001555001576001597001618001639001660001681001702001723

15、00174400176500178600180700182800184900187000189100191200193300195400197500199600201700203800205900208000210100210201210301210401210501210601210701210801210901211001211101211201211301211401211501211601211701211801211901212001212101*材料信息數(shù)據(jù)類型號 彈性模量 泊松比 容重 單元厚度1 10000.0 0.3 0.0 1.0 *荷載信息數(shù)據(jù)受集中力作用的結(jié)點(diǎn)數(shù) 面力批

16、數(shù)0 1集中力數(shù)據(jù)受集中力作用的結(jié)點(diǎn)號 X向集中力分量 Y向集中力分量*面力數(shù)據(jù)面力批號 受面力單元個(gè)數(shù) 面力類型 面力集度 最大集度 面力作用面號 1 15 2 -10 -10 4*面力作用的單元的單元號200019801960194019201900188018601840182018001780176017401720（3）源程序： PROGRAM FEM DIMENSION SK(200000),COOR(2,10000),AE(4,100),MEL(5,10000),WG(4),&JR(2,10000),MA(10000),R(10000),IEW(10000),STRE(3,100

17、00),&TOALFUN(10000,10000),SUMS(10000),GOODSTR(10000) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) OPEN(1,FILE=INDAT.DAT,STATUS=OLD) OPEN(2,FILE=OUT.DAT,STATUS=NEW) CALL INPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR) !輸入數(shù)據(jù) CALL HUAFEN(XL,B,NNL,NNB,COOR,MEL) CALL CBAND (M

18、A,JR,MEL) !形成主元素序號指示矩陣MA DO I=1,NH SK(I)=0.0 ENDDO CALL SK0(SK,MEL,COOR,JR,MA,AE) !形成整體剛度矩陣K DO I=1,N R(I)=0.0 ENDDO IF(NCP.GT.0) CALL CONCR(NCP,R,JR) CALL BODYR(R,MEL,COOR,JR,AE) READ(1,70) TL READ(1,70) TL READ(1,70) TL70 FORMAT(A) IF(iz.GT.0)then DO jj=1,iz READ (1,*)JS,NSE,(WG(I),I=1,4) READ(1,7

19、0) TL READ(1,70) TL READ(1,*)(IEW(M),M=1,NSE) CALL FACER(IEW,NSE,R,MEL,COOR,JR,WG) ENDDO ENDIF CALL DECOP (SK,MA) CALL FOBA (SK,MA,R) CALL OUTDISP(NP,R,JR) CALL STRESS (COOR,MEL,JR,AE,R,STRE) CALL GAUSTRSS(MEL,COOR,AE,TOALFUN,SUMS,GOODSTR,STRE,JR,R)!應(yīng)力磨平并輸出磨平后的應(yīng)力值 WRITE(2,(A) PROGRAM SAFF HAS BEEN E

20、NDED WRITE(*,(A) PROGRAM SAFF HAS BEEN ENDED CLOSE(1) CLOSE(2) END!INPUT為原始數(shù)據(jù)輸入子程序 SUBROUTINE INPUT(XL,B,IX,IY,AE,NCP,IZ,NNL,NNB,JR) DIMENSION AE(4,*),JR(2,*) CHARACTER*400TL COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH READ(1,70)TL READ(1,70)TL READ(1,*)XL,B,NNL,NNB,NM,NR WRITE(2,10)XL,B,NNL,NNB

21、,NM,NR10 FORMAT(5X,平面四邊形四結(jié)點(diǎn)單元有限元程序/5X,L=,F5.2,4X,&B=,F5.2,4X,NNL=,I3,4X,NNB=,I3,4X,NM=,I2,4X,NR=,I3) NP=NNL*NNB DO 15 I=1,NP DO 15 J=1,2 15 JR(J,I)=1 READ(1,70)TL READ(1,70)TL READ(1,70)TL WRITE(2,110)110 FORMAT(/5X,約束信息/5X,結(jié)點(diǎn)號,5X,X向約束,5X,Y向約束) DO 20 I=1,NR READ(1,*)IP,IX,IY WRITE(2,100)IP,IX,IY100

22、FORMAT(8X,I5,5X,I2,9X,I2) JR(1,IP)=IX JR(2,IP)=IY 20 CONTINUE N=0 DO 30 I=1,NP DO 30 J=1,2 IF (JR(J,I) 30,30,25 25 N=N+1 JR(J,I)=N 30 CONTINUE READ(1,70) TL READ(1,70) TL READ(1,70) TL DO 55 J=1,NM READ (1,*) JJ,(AE(I,JJ),I=1,4) WRITE(2,910) JJ,(AE(I,JJ),I=1,4)55 CONTINUE910 FORMAT (/5X,材料信息/5X,類型號,

23、5X,彈性模量,5X,泊松比,5X,容重& ,8X,單元厚度/5X,I5,4(5x,E8.3) READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)NCP,IZ WRITE(2,920)NCP,IZ920 FORMAT (/5X,荷載信息/5X,受集中力作用的結(jié)點(diǎn)數(shù),8X,面力批數(shù)&/(2(12X,I5) READ(1,70) TL READ(1,70) TL70 FORMAT(A) RETURN END !自動(dòng)劃分網(wǎng)格子程序，計(jì)算單元數(shù)，結(jié)點(diǎn)數(shù)，結(jié)點(diǎn)坐標(biāo)，每個(gè)單元包含的結(jié)點(diǎn)!單元的4個(gè)結(jié)點(diǎn)號與材料類型號，從左下角逆時(shí)針編號順序； SUBROU

24、TINE HUAFEN(XL,B,NNL,NNB,COOR,MEL) DIMENSION COOR(2,*),MEL(5,*) COMMON /CMN1/ NP,NE,NM,NR NP=NNL*NNB NE=(NNL-1)*(NNB-1) WRITE(2,*)總結(jié)點(diǎn)數(shù) NP=,NP WRITE(2,*)總單元數(shù) NE=,NE DL=XL/(NNL-1) DB=B/(NNB-1) N=1 M=1 T=NNL*NNB DO 10 K=NNB,T,NNB DO 20 I=N,K IP=I COOR(1,I)=DL*(K-NNB)/NNB20 CONTINUE N=K+110 CONTINUE DO

25、30 J=NNB,T,NNB DO 40 I=M,J COOR(2,I)=DB*(I-M)40 CONTINUE M=J+130 CONTINUE WRITE(2,100) (J,(COOR(I,J),I=1,2),J=1,NP)100 FORMAT(/7X,結(jié)點(diǎn)號,10X,X坐標(biāo),8X,Y坐標(biāo)/(8X,I5,5X,2F12.6) N=1 DO 50 K=1,(NNL-1) DO 60 I=N,(NNB-1)*K NNE=I MEL(1,I)=I+K-1 MEL(2,I)=MEL(1,I)+NNB MEL(3,I)=MEL(1,I)+NNB+1 MEL(4,I)=MEL(1,I)+160 CO

26、NTINUE N=N+(NNB-1)50 CONTINUE WRITE(2,111)(J,(MEL(I,J),I=1,4),J=1,NE)111 FORMAT(/7X,單元號,5X,結(jié)點(diǎn)1,5X,結(jié)點(diǎn)2,5X,結(jié)點(diǎn)3,5X,結(jié)點(diǎn)&4/(7X,I5,5X,I5,5X,I5,5X,I5,5X,I5) DO 70 IE=1,NE MEL(5,IE)=170 CONTINUE END!*! LU三角分解子程序! * SUBROUTINE DECOP (SK,MA) !方程LU分解 DIMENSION SK(*),MA(*) COMMON /CMN2/ N,MX,NH DO 50 I=2,N L=I-M

27、A(I)+MA(I-1)+1 K=I-1 L1=L+1 IF (L1.GT.K) GOTO 30 DO 20 J=L1,K IJ=MA(I)-I+J M=J-MA(J)+MA(J-1)+1 IF (L.GT.M) M=L MP=J-1 IF (M.GT.MP) GO TO 20 DO 10 LP=M,MP IP=MA(I)-I+LP JP=MA(J)-J+LP SK(IJ)=SK(IJ)-SK(IP)*SK(JP) 10 CONTINUE 20 CONTINUE 30 IF (L.GT.K) GOTO 50 DO 40 LP=L,K IP=MA(I)-I+LP LPP=MA(LP) SK(IP

28、)=SK(IP)/SK(LPP) II=MA(I) SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP) 40 CONTINUE 50 CONTINUE RETURN END!*! LU三角分解后的回代求解! * SUBROUTINE FOBA (SK,MA,R) !回代求解 DIMENSION SK(*),MA(*),R(*) COMMON /CMN2/ N,MX,NH DO 10 I=2,N L=I-MA(I)+MA(I-1)+1 K=I-1 IF (L.GT.K) GOTO 10 DO 5 LP=L,K IP=MA(I)-I+LP R(I)=R(I)-SK(IP)*R(

29、LP) 5 CONTINUE 10 CONTINUE DO 20 I=1,N II=MA(I)45 R(I)=R(I)/SK(II)20 CONTINUE DO 30 J1=2,N I=2+N-J1 L=I-MA(I)+MA(I-1)+1 K=I-1 IF (L.GT.K) GO TO 30 DO 25 J=L,K IJ=MA(I)-I+J55 R(J)=R(J)-SK(IJ)*R(I)25 CONTINUE30 CONTINUE RETURN END!*！ 計(jì)算形函數(shù)對整體坐標(biāo)的導(dǎo)數(shù)! * SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,IVEN,NEE) DIM

30、ENSION XYZ(2,*),DNX(2,4),RJAC(2,2),IVEN(*) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) CALL FUN8 (XYZ,R,S,DET)!求形函數(shù) IF (DET.LT.1.0E-5)THEN !判斷J的大小 WRITE(2,600) NEE,R,S,det WRITE(2,*) (IVEN(M),M=1,4) STOP ENDIF REC=1.00/DET !求J的逆 RJAC(1,1)=REC*XJAC(2,2) RJAC(2,2)=REC*XJAC(1,1) RJAC(2,1)=-REC*XJAC(2,1) RJA

31、C(1,2)=-REC*XJAC(1,2) DO 30 K=1,4 !求形函數(shù)的導(dǎo)數(shù)，坐標(biāo)轉(zhuǎn)換后 DO 20 I=1,2 DNX(I,K)=0. DO 25 M=1,2 DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K) 25 CONTINUE 20 CONTINUE 30 CONTINUE600 FORMAT(1X,ERR0R* NEGTIVE OR ZERO JACOBIAN DETERMINANT FOR &ELEMENT/ELE.=,I5,R=,F10.5,6X,S=,F10.5,det=,f12.5) RETURN END!*！ !形函數(shù)子程序，計(jì)算雅可比行列式!

32、 * SUBROUTINE FUN8 (XYZ,R,S,DET) DIMENSION XYZ(2,*),XI(4),ETA(4) COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) DATA XI/-1.0,1.0,1.0,-1.0/!雙線性自然坐標(biāo)值 DATA ETA/-1.0,-1.0,1.0,1.0/ DO 10 I=1,4 G1=(1.0+XI(I)*R) G2=(1.0+ETA(I)*S) FUN(I)=0.25*G1*G2 PN(1,I)=0.25*XI(I)*G2 !自然坐標(biāo)下的形函數(shù)導(dǎo)數(shù) PN(2,I)=0.25*ETA(I)*G1 10 CONTIN

33、UE DO 80 I=1,2 DO 75 J=1,2 DET=0.00 DO 70 K=1,4 DET=DET+PN(I,K)*XYZ(J,K) !求J矩陣的各項(xiàng)值 70 CONTINUE XJAC(I,J)=DET 75 CONTINUE 80 CONTINUE DET=XJAC(1,1)*XJAC(2,2)-XJAC(2,1)*XJAC(1,2)!求J矩陣的行列式值 RETURN END!*！ 計(jì)算單元?jiǎng)偠染仃囎映绦? *SUBROUTINE STIF(XYZ,AE,IVEN) DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),IVEN(*),RJAC(2,2) CO

34、MMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /CMN5/ FUN(4),PN(2,4),XJAC(2,2) COMMON /GAUSS/ RSTG(3),H(3) DO 40 I=1,8 RF(I)=0.00 DO 30 J=1,8 SKE(I,J)=0.00!賦初值 30 CONTINUE 40 CONTINUE E=AE(1,NME) !單元參數(shù)，E，v,重度 U=AE(2,NME) GAMMA=AE(3,NME)

35、D1=E*(1.00-U)/(1.00+U)*(1.00-2.00*U) !剛度系數(shù)矩陣的一些系數(shù) D2=E*U/(1.00+U)*(1.00-2.00*U) D3=E*0.50/(1.00+U) DO 120 I=1,4 II=2*(I-1) I1=II+1 I2=II+2 DO 115 J=1,4 JJ=2*(J-1) J1=JJ+1 J2=JJ+2 DXX=0 DXY=0 DYX=0 DYY=0 DO 99 IS=1,3 S=RSTG(IS) SH=H(IS) DO 98 IR=1,3 R=RSTG(IR) RH=H(IR) CALL FDNX (XYZ,DNX,DET,R,S,RJAC

36、,IVEN,NEE) DNIX=DNX(1,I) DNIY=DNX(2,I) DNJX=DNX(1,J) DNJY=DNX(2,J) DXX=DXX+DNIX*DNJX*DET*RH*SH !按三點(diǎn)高斯積分求解單元?jiǎng)偠戎?DXY=DXY+DNIX*DNJY*DET*RH*SH DYX=DYX+DNIY*DNJX*DET*RH*SH DYY=DYY+DNIY*DNJY*DET*RH*SH 98 CONTINUE 99 CONTINUE SKE(I1,J1)=DXX*D1+DYY*D3 !單元?jiǎng)偠?SKE(I2,J2)=DYY*D1+DXX*D3 SKE(I1,J2)=DXY*D2+DYX*D3

37、SKE(I2,J1)=DYX*D2+DXY*D3115 CONTINUE120 CONTINUE RETURN END!*！ 形成主元素序號指示矩陣MA（*）! * SUBROUTINE CBAND (MA,JR,MEL) DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH DO 65 I=1,N65 MA(I)=0 DO 90 IE=1,NE !對桿件循環(huán) DO 75 K=1,4 IEK=MEL(K,IE) !單元結(jié)點(diǎn)號 DO 95 M=1,2 JJ=2*(K-1)+M

38、NN(JJ)=JR(M,IEK) !NN八個(gè)元素，記錄單元每個(gè)結(jié)點(diǎn)自由度信息95 CONTINUE75 CONTINUE L=N DO 80 I=1,8 NNI=NN(I) IF(NNI.EQ.0) GOTO 80 IF(NNI.LT.L) L=NNI !找單元的最小自由度號 80 CONTINUE DO 85 M=1,8 JP=NN(M) IF(JP.EQ.0) GOTO 85 !排除約束的 JPL=JP-L+1 !把單元自由度號從1開始 IF(JPL.GT.MA(JP) MA(JP)=JPL !公共結(jié)點(diǎn)找最大的自由度號 85 CONTINUE 90 CONTINUE MX=0 MA(1)=

39、1 DO 10 I=2,N IF(MA(I).GT.MX) MX=MA(I) !找總體最大自由度號，也是半帶寬 MA(I)=MA(I)+MA(I-1) 10 CONTINUE NH=MA(N) WRITE (*,500) N,MX,NH WRITE (2,500) N,MX,NH500 FORMAT (/5X,總自由度數(shù) N=,I5,3X,半帶寬 MX=,I5,3X,一維存儲(chǔ)位數(shù) &NH=,I7) RETURN END!*！ !形成總體剛度矩陣子程序! * SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE) DIMENSION SK(*),MEL(5,*),COOR(2,

40、*),JR(2,*),MA(*),AE(4,*),&XYZ(2,4),IVEN(4) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) COMMON /CMN4/ NEE,NME COMMON /GAUSS/ RSTG(3),H(3) H(1)=0.5555555555555560 H(2)=0.8888888888888890 H(3)=H(1) RSTG(1)=-0.7745966692414830 !三點(diǎn)高斯積分 RSTG(2)=0.00 RSTG(3)=-RSTG(1) DO 10 IE=1,NE NEE=IE !單元號 NME=MEL(5,IE) !材料型號 DO 75 K=1,4 IEK=MEL(K,IE) !單元結(jié)點(diǎn)號碼 IVEN(K)=IEK DO 95 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK) !NN八個(gè)元素，記錄單元每個(gè)結(jié)點(diǎn)自由度號95 XYZ(M,K)=COOR(M,IEK) !結(jié)點(diǎn)坐標(biāo)75 CONTINUE CALL STIF(XYZ,AE,IVEN) !單剛 DO 60 I=1,8 DO 60 J=1,8 II=NN(I) JJ=NN(J) IF (JJ.EQ.0).OR.(II.L