有限元程序設(shè)計課件

教學(xué)程序：FEP1李明瑞編FEP2李明瑞編ZFEP教學(xué)程序：1第一章緒論有限元程序的基本內(nèi)容有限元求解器第一章緒論有限元程序的基本內(nèi)容2§1.1有限元程序的基本內(nèi)容有限元法的解題步驟結(jié)構(gòu)離散化為有限元網(wǎng)格計算單元剛度矩陣引入約束條件求應(yīng)力等解方程組建立總剛度矩陣及外載§1.1有限元程序的基本內(nèi)容有限元法的解題步驟結(jié)構(gòu)離散化3有限元程序的基本內(nèi)容：數(shù)據(jù)輸入階段有限元矩陣的計算、組集和求解數(shù)據(jù)輸出階段前處理處理器后處理有限元程序的基本內(nèi)容：數(shù)據(jù)輸入階段有限元矩陣的計算、組集和求4一、數(shù)據(jù)輸入階段—前處理

讀入和生成數(shù)據(jù)，形成有限元網(wǎng)格，為有限元矩陣計算作準備。（數(shù)據(jù)或圖形輸入）標題和控制信息計算存貯分配節(jié)點坐標和約束信息單元信息材料信息載荷信息一、數(shù)據(jù)輸入階段—前處理讀入和生成數(shù)據(jù)，形成有限5二、有限元矩陣的計算、組集和求解—處理器計算插值函數(shù)矩陣[N

]、幾何矩陣[B

]、Jacob矩陣[J

]在高斯點進行數(shù)值積分，求得單剛、單元載荷組集成總剛、總載求解二、有限元矩陣的計算、組集和求解—處理器計算插值函數(shù)矩陣[6三、數(shù)據(jù)輸出階段—后處理輸出結(jié)果（場變量和場變量導(dǎo)數(shù)等），打印結(jié)果。場變量包括：位移、溫度、流場的速度勢等；場變量導(dǎo)數(shù)包括：應(yīng)力、應(yīng)變、熱流、流場速度勢等輸出形式：數(shù)據(jù)輸出和圖形輸出三、數(shù)據(jù)輸出階段—后處理輸出結(jié)果（場變量和場變量導(dǎo)數(shù)等），打7§1-2有限元求解器帶寬法：只存儲總剛矩陣的半帶寬內(nèi)的元素。等帶寬法：每行的帶寬相等，常采用二維數(shù)組存儲變帶寬法：每行的帶寬不等，常采用一維數(shù)組存儲一、帶寬法變帶寬分塊求解思想：一個變量的消元只影響剛度矩陣的一部分元素，較大的節(jié)點分量方程組可以分成較小的局部系數(shù)矩陣按節(jié)點逐步求解?！?-2有限元求解器帶寬法：只存儲總剛矩陣的半帶寬內(nèi)的元8

波前法是另一種分塊儲存的變帶寬法，其消元次序是按單元編號進行，邊組集邊消元。調(diào)入內(nèi)存的單元所保留的節(jié)點稱為波前節(jié)點，所消去的節(jié)點稱為消元節(jié)點，消元節(jié)點是與以后調(diào)入的單元沒有聯(lián)系的節(jié)點，即該節(jié)點的方程已組集完。波前法的回代按消元節(jié)點的反序進行。對內(nèi)存的需求取決與最大波前寬。二、波前法波前法是另一種分塊儲存的變帶寬法，其消9V-FortranFortran語言的基本格式變量基本語句子例子程序函數(shù)子程序其它功能、模塊V-FortranFortran語言的基本格式10有限元程序設(shè)計課件11有限元程序設(shè)計課件12有限元程序設(shè)計課件13第二章FEP2程序的總體設(shè)計與輸入數(shù)據(jù)

FEP2的設(shè)計任務(wù)

FEP2的結(jié)構(gòu)設(shè)計

FEP2的主程序

FEP2的主控程序

FEP2的數(shù)據(jù)輸入格式與程序?qū)崿F(xiàn)第二章FEP2程序的總體設(shè)計與輸入數(shù)據(jù)FEP214§2-1FEP2的設(shè)計任務(wù)1.FEP2的簡要題目說明FEP2是一個具有通用性的教學(xué)程序，可用于計算一般的線性靜力學(xué)問題。已設(shè)計了平面梁單元與平面3-9節(jié)點元兩種單元，但留下接口。2.支持軟件與硬件FORTRAN77以上編譯器、各種微機3.FEP2的功能§2-1FEP2的設(shè)計任務(wù)1.FEP2的簡要題目說153.FEP2的功能1）輸入文件名由用戶自由定義，但限制為4個字符，輸入文件擴展名為“DAT”，輸出文件擴展名為“OUT”。2）節(jié)點坐標、單元信息等具有線性自動生成功能。3）可以處理多種工況、多種類型單元組合結(jié)構(gòu)問題。4）可以處理固定約束和指定位移。5）采用變帶寬存儲、三角分解法求解剛度方程。6）為多種單元留下接口。3.FEP2的功能16§2-2FEP2的結(jié)構(gòu)設(shè)計FEP2：由形式主程序、主控程序、功能模塊組成。主程序的主要功能：定義輸入、輸出文件名，調(diào)用主控程序PCONTR主控程序PCONTR的主要功能模塊：1）內(nèi)存分配2）網(wǎng)格生成3）變帶寬存儲設(shè)置4）單剛的形成與組集5）剛度矩陣的分解6）形成節(jié)點載荷7）形成與組集單元載荷8）求解位移9）求單元應(yīng)力10）求結(jié)構(gòu)反力§2-2FEP2的結(jié)構(gòu)設(shè)計FEP2：由形式主程序、主控17FEP2程序框圖FEP2（主程序）PCONTR(主控程序)—SETMEM檢查內(nèi)存—PROFIL形成變帶寬剛度矩陣地址

—PFORM(3)形成單剛并組集—LDLT總剛的三角分解

—GENVEC形成節(jié)點載荷—PFOM(3)構(gòu)造并組集單元載荷

—FORBACK前消回代求出位移—PFORM(4)計算單元應(yīng)力

—PFORM(6)計算結(jié)構(gòu)反力FEP2程序框圖FEP2（主程序）PCONTR—S18§2-3FEP2的主程序一、文件名FEP2的文件名可由用戶自由定義（限制為4個字符），通過人機交互方式確定。設(shè)計中引入以下技巧：1)規(guī)定輸入文件名FI與輸出文件名FO為8個字符2)規(guī)定兩個字符數(shù)組NAMINP(8)，NAMOUT(8)3)用IQUIVALENCE語句使FI與NAMINP、FO與NAMOUT等價4)用DATA語句給FI和FO賦初值:“.DAT”、“.OUT”5)輸入NAMINP的前四個字符作為輸入輸出文件名§2-3FEP2的主程序一、文件名1)規(guī)定輸入文件19COMMON/PSIZE/MAXM,MAXACHARACTER*8FI,FOCHARACTERNAMINP(8),NAMOUT(8),HEAD*50COMMON/HEAD/HEAD1EQUIVALENCE(FI,NAMINP(1)),(FO,NAMOUT(1))DATAFI,FO/'.DAT','.OUT'/WRITE(*,'(A\)')'INPUTFILENAME(4LETTERSONLY):'READ(*,'(4A1)')(NAMINP(I),I=1,4)DO10I=1,410NAMOUT(I)=NAMINP(I)

COMMON/PSIZE/MAXM,MAXA20

OPEN(6,FILE=FO)OPEN(5,FILE=FI)MAXM=16000MAXA=16380CALLPCONTR(FO)END

OPEN(6,FILE=FO)21二、定義數(shù)組FEP2中設(shè)置了兩個數(shù)組M(1600)和A(16380)，數(shù)組M是作為存放坐標、單元信息、載荷、位移等，數(shù)組A則是存放剛度矩陣用的，其上界與FORTRAN（編譯器）版本和計算機內(nèi)存有關(guān)，確定了程序的計算規(guī)模。這兩個數(shù)組也可合并成一個。二、定義數(shù)組22§2-4FEP2的主控程序

PCONTR：實質(zhì)上的主程序。分以下幾大段：一、程序頭（前12語句）SUBROUTINEPCONTR(FO)CHARACTER*8FOLOGICALERR,ASEMCOMMON/M/M(6000)COMMON/A/A(16380)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/PSIZE/MAXM,MAXACOMMON/MDATA/NUMMATCHARACTERFD*25,FORC*6,PCOMMON/PRINT/PCOMMON/ASEM/ASEMDATAFORC/'FORC'/FD/'NODALFORCE/DISPL'/§2-4FEP2的主控程序PCONTR：23FO：輸出文件名ERR：邏輯變量，作出錯信息控制（.TRUE.為錯）ASEM：邏輯變量，作為判斷是否要組集總剛數(shù)組A：存放總剛的元素數(shù)組M：存放單剛、單元信息、節(jié)點坐標等NDF—最大自由度數(shù)NDM—空間維數(shù)NEN—單元的最多節(jié)點數(shù)NJ—節(jié)點總數(shù)NE—總單元數(shù)NEQ—總方程數(shù)NUMMAT—材料類型數(shù)FD,FORC：文字變量，作為輸出的標題用P：文字變量，判斷是否要輸出單剛FO：輸出文件名24二、輸入控制信息、確定主要數(shù)組的地址（14—34語句）NEN1=NEN+1每個單元信息的存儲量（節(jié)點號+材料號）NNEQ=NJ*NDF總節(jié)點節(jié)點自由度數(shù)NST=NDF*NEN單剛階數(shù)NXC=1節(jié)點坐標數(shù)組XC的首址NIX=NXC+NJ*NDM單元信息數(shù)組IX的首址NID=NIX+NEN1*NE結(jié)構(gòu)方程號編碼數(shù)組ID的首址ND=NID+NNEQ材料常數(shù)數(shù)組的首址NXL=ND+20*NUMMAT單元節(jié)點坐標數(shù)組XL的首址NLD=NXL+NEN*NDM單元節(jié)點自由度數(shù)組LD的首址NUL=NLD+NST單元節(jié)點位移數(shù)組UL的首址NS=NUL+NST*2單剛數(shù)組S的首址NP=NS+NST*NST總剛對角線地址向量JP的首址

……二、輸入控制信息、確定主要數(shù)組的地址（14—34語句）NE25READ(5,*)NJ,NE,NUMMAT,NDF,NDM,NENWRITE(*,2000)NJ,NE,NUMMAT,NDF,NDM,NENWRITE(6,2000)NJ,NE,NUMMAT,NDF,NDM,NENNEN1=NEN+1NNEQ=NJ*NDFNST=NDF*NENNXC=1NIX=NXC+NJ*NDM

……NS=NUL+NST*2NP=NS+NST*NSTNF=NP+NSTNU=NF+NNEQNJP=NU+NNEQNEND=NJP+NNEQREAD(5,*)NJ,NE,NUMM26§2-5FEP2的數(shù)據(jù)輸入格式與程序?qū)崿F(xiàn)FEP2的數(shù)據(jù)信息（自由格式）：1）控制信息：NJ，NE，NUMMAT，NDM，NEN2）坐標信息（一組）N，NG，（XL(I)，I=1，NDM）

節(jié)點號節(jié)點號增量坐標或載荷由GENVEC生成節(jié)點坐標或載荷向量，通過NDM,CD，F(xiàn)C進行識別生成的是節(jié)點坐標，還是載荷向量。CD：’NODALCOORDINATES’FC：’COOR’CD：’NODALFORCE/DISPL’FC：’FORC’§2-5FEP2的數(shù)據(jù)輸入格式與程序?qū)崿F(xiàn)FEP2的數(shù)據(jù)信27

SUBROUTINEGENVEC(NDM,X,CD,FC,NJ,ERR)CHARACTERCD*18,FC*6LOGICALERRREALX(NDM,*),XL(6)ERR=.FALSE.N=0NG=0102L=NLG=NGREAD(5,*,ERR=999)N,NG,(XL(I),I=1,NDM)101IF(N.LE.0.OR.N.GT.NJ)GOTO109DO103I=1,NDM103X(I,N)=XL(I)IF(LG)104,102,104104LG=ISIGN(LG,N-L)SUBROUTINEGENVEC(NDM,X,CD28LI=(IABS(N-L+LG)-1)/IABS(LG)DO105I=1,NDM105XL(I)=(X(I,N)-X(I,L))/LI106L=L+LGIF((N-L)*LG.LE.0)GOTO102IF(L.LE.0.OR.L.GT.NJ)GOTO108DO107I=1,NDMLMLG=L-LG107X(I,L)=X(I,LMLG)+XL(I)GOTO106WRITE(6,3000)L,CDWRITE(*,3000)L,CDERR=.TRUE.GOTO102CONTINUE……LI=(IABS(N-L+LG)-1)293）單元信息（一組）L，LX，LK，（IDL(K)，K=1，NEN）單元號單元號增量材料號節(jié)點號0，0，0，……結(jié)束由ELDA子程序生成單元信息，IX(NEN1,NE)二維數(shù)組。NEN=NEN+1NE—單元總數(shù)3）單元信息（一組）304）材料信息（一組或多組）對于每一類材料，要求輸入兩個記錄：a）材料類型號，單元類型號b）一批與單元類型相關(guān)的常數(shù)（不超過20個）對于梁單元，輸入兩個如下12個常數(shù)：NP,E,G,A或Ix,Iy或Iz，Rh0，As，N1，N2，Ni，W,a類型（NP=0：平面梁作用平面載荷；1：平面梁作用空間載荷）Rh0：單位長度的質(zhì)量密度N1：梁左端鉸為1，非鉸為0N2：右端鉸為1，非鉸為0Ni：單元載荷類型（0，1，2，3，4）（0為無載荷）單元載荷由EQLOAD處理成等效節(jié)點載荷4）材料信息（一組或多組）對于梁單元，輸入兩個如下12個31對于平面元，輸入6個常數(shù)：E,XNU,TH，L，K，ITH：單元厚度L：沿一個方向的高斯積分點數(shù)（2，3，4）K：沿一個方向應(yīng)力輸出的高斯點數(shù)I：=0平面應(yīng)變0平面應(yīng)力PMESH調(diào)用ELEMLIB，再調(diào)用BEAMorPLANE對于平面元，輸入6個常數(shù)：PMESH調(diào)用ELEM32SUBROUTINEPMESH(XC,IX,ID,D,NEN1)CHARACTERCOOR*6,XX*18COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/MDATA/NUMMATCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONXC(NDM,*),IX(NEN1,*),ID(NDF,*),D(20,*)LOGICALERRDATAXX/'NODALCOORDINATES'/COOR/'COOR'/1CALLGENVEC(NDM,XC,XX,COOR,NJ,ERR)

SUBROUTINEPMESH(XC,IX,ID331CALLGENVEC(NDM,XC,XX,COOR,NJ,ERR)IF(ERR)WRITE(*,3000)IF(ERR)STOP2CALLELDAT(IX,NEN1)3DO304N=1,NUMMATWRITE(6,2000)MA,IELWRITE(*,2000)MA,IELIF(MA.EQ.0)GOTO4CALLPZERO(D(1,MA),20)D(20,MA)=IELCALLLEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,NST,1)304CONTINUE2000FORMAT(/'MATERIALTYPE',I4,'ELEMENTTYPE',I4/)4CALLBOUN(ID)END1CALLGENVEC(NDM,XC,XX,CO345）約束信息（包括指定位移）（一組）PMESH調(diào)用BOUN輸入：N，NG，（IDL(I)，I=1，NDF）

節(jié)點號節(jié)點號增量坐標或載荷0，0，0，……（平面元4個0，梁元5個0）IDL(I)：約束信息，被約束的自由度給非零值（1或-1），其它為0。1：該自由度被約束，但不繼續(xù)生成-1：該自由度被約束，要求繼續(xù)生成5）約束信息（包括指定位移）（一組）356）載荷信息載荷輸入仍由GENVEC完成，維數(shù)：NDF

對于指定位移，生成方式與節(jié)點載荷相同，5）輸入約束信息，6）輸入其值。6）載荷信息對于指定位移，生成方式與節(jié)點載荷相36第三章總剛度矩陣的變帶寬存貯與求解總剛度矩陣的變帶寬存貯與對角線地址LDLT三角分解自由度指示矩陣ID與對角線地址矩陣第三章總剛度矩陣的變帶寬存貯與求解總剛度矩陣的變帶寬存37§3-1總剛度矩陣的變帶寬存貯與對角線地址總剛度矩陣的性質(zhì)：對稱、正定、稀疏只存貯半帶寬內(nèi)的元素（下三角、變帶寬）稀疏對稱矩陣A中的元素可用兩個一維數(shù)組B和JP完全確定。B的上界為A中的下三角、變帶寬內(nèi)元素總數(shù)，存貯A內(nèi)的元素。JP的上界為A的階數(shù)，存貯A的各行對角線元素的序號，稱為對角線元素地址數(shù)組。對應(yīng)關(guān)系：A(I,J)=B(JP(I)–I+J)每一行第一個非零元素的列號Mi:M1=1,Mi=I–JP(I)+JP(I-1)+1§3-1總剛度矩陣的變帶寬存貯與對角線地址總剛度矩陣的性質(zhì)38SUBROUTINEPROFIL(JP,ID,IX,NEN1,NK)

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

DIMENSIONJP(1),ID(NDF,1),IX(NEN1,*),IDL(18)

NEQ=0

NAD=0

DO50N=1,NJ

DO40I=1,NDF

J=ID(I,N)

IF(J)30,20,30

20NEQ=NEQ+1

ID(I,N)=NEQ

JP(NEQ)=0

GOTO40

30NAD=NAD-1

ID(I,N)=NAD

40CONTINUE

50CONTINUESUBROUTINEPROFIL(JP,ID39C.....COMPUTEROWBANDWIDTHDO500N=1,NEDO400I=1,NENII=IX(I,N)IF(II.EQ.0)GOTO400DO300K=1,NDFKK=ID(K,II)IF(KK.LE.0)GOTO300DO200J=1,NENJJ=IX(J,N)IF(JJ.EQ.0)GOTO200DO100L=1,NDFLL=ID(L,JJ)IF(LL.LE.0)GOTO100M=MAX0(KK,LL)JP(M)=MAX0(JP(M),IABS(KK-LL))100CONTINUE200CONTINUE300CONTINUECONTINUE500CONTINUEC.....COMPUTEROWBANDWIDTH40C.....COMPUTEDIAGONALPOINTERSFORPROFILNK=1JP(1)=1IF(NEQ.EQ.1)RETURNDO600N=2,NEQ600JP(N)=JP(N)+JP(N-1)+1NK=JP(NEQ)WRITE(*,2000)NK2000FORMAT('THEMAXIMUMSTORAGEOFGLOBALSTIFFNESS=',I6/)RETURNENDC.....COMPUTEDIAGONALPOINTER41

SUBROUTINEADDSTF(LD,JP,S,NST)CHARACTERYCOMMON/PRINT/YCOMMON/A/A(16380)COMMON/CDATA/NDF,NDM,NEN,NEJ,NE,NEQCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONLD(*),JP(*),S(NST,NST)IF(Y.EQ.'Y'.OR.Y.EQ.'y')WRITE(*,2000)N,SFORMAT('ELEMENTNUMBER',I4,'ELEMENTSTIFFNESS'/(6E12.5))SUBROUTINEADDSTF(LD,J42

DO200J=1,NSTK=LD(J)IF(K.LE.0)GOTO200L=JP(K)-KDO100I=1,NSTM=LD(I)IF(M.GT.K.OR.M.LE.0)GOTO100M=L+MA(M)=A(M)+S(J,I)100CONTINUE200CONTINUERETURNENDDO200J=1,NST43§3-2LDLT三角分解求解方程：Ax=b

A—對稱、正定矩陣=>Ux=yU—單位上三角矩陣對A進行三角分解=>A=LLTL—下三角矩陣，通過比較系數(shù)法確定，需要開方運算LDLT三角分解=>A=LDLTL—下三角矩陣，D—對角矩陣i=1,…,n;j=2,…,iLDLT(A,JP,N)DIMENSIONA(*),JP(*)§3-2LDLT三角分解求解方程：Ax=b44Ax=b=>LDLTx=b令LTx=y=>LDy=b

前消公式：i=2,…,n回代公式：i=n-1,…,1

SUBROUTINEFORBACK(A,B,JP,N)DIMENSIONA(*),B(*),JP(*)Ax=b=>45§3-3自由度指示矩陣ID與對角線地址矩陣ID(NDF,NJ)：非零值（1或-1）為被約束的自由度，0為活化自由度。在PROFIL中，改造數(shù)組ID，原來的零元素（活化自由度）用正數(shù)計，原來的非零元素（非活化自由度）用負數(shù)計，全部零元素（活化）的個數(shù)記為NEQ（總剛矩陣的階數(shù)）。§3-3自由度指示矩陣ID與對角線地址矩陣ID(NDF,46計算剛度矩陣的帶寬的方案比較：1）以各節(jié)點自由度為研究對象，利用單元信息IX，找出與該點聯(lián)系的節(jié)點，然后比較節(jié)點自由度之最大者。運算量：NJ*NDF*NE*NEN*NDF2）逐一研究各單元，分別取出各節(jié)點的各個自由度與該單元中其它節(jié)點的自由度比較，差值最大者再與所研究的自由度之帶寬值（初值為零）比較，取其中最大者代替原有的帶寬值。運算量：NE*NEN*NDF*NEN*NDF兩者運算量比：NJ:NE選方案2帶寬=自由度之最大差值+1暫存放于JP中計算剛度矩陣的帶寬的方案比較：1）以各節(jié)點自由度為研究對象，47若已知對角線地址數(shù)組JP，則第I行的帶寬：JP(1)=1,JP(I)–JP(I–1)

若JP存放的是帶寬,則對角線地址數(shù)組：JP(1)=JP(1)=1JP(2)=JP(2)+JP(1)……JP(I+1)=JP(I+1)+JP(I)總剛矩陣(按活化自由度存貯)全部存貯量:NK=JP(NEQ)若已知對角線地址數(shù)組JP，則第I行的帶寬：48第四章單剛與載荷的組集，指定約束的處理，單元分析的預(yù)處理子程序PFORM的功能子程序MODIFY完成將指定位移轉(zhuǎn)化成等價內(nèi)力將單剛組集成總剛子程序ADDSTF子程序BASBLY組集載荷與節(jié)點反力向量PRTREAC輸出節(jié)點反力PRTDIS輸出節(jié)點位移ELEMLIB單元庫管理程序第四章單剛與載荷的組集，指定約束的處理，單元分析的預(yù)處49§4-1子程序PFORM的功能子程序PFORM起到承上啟下的作用，為調(diào)用與單元運算有關(guān)的各種子程序做必要的前后處理。功能參數(shù)ISW、ASEM：ISW=3，ASEM=.TRUE.構(gòu)造單剛，并組集ISW=3，ASEM=.FALSE.構(gòu)造單元載荷，組集，將指定位移轉(zhuǎn)化為成等效載荷ISW=4計算單元內(nèi)力或應(yīng)力，并輸出ISW=6構(gòu)造單元內(nèi)力，并組集成結(jié)構(gòu)反力ISW=1讀入并構(gòu)造有關(guān)的材料常數(shù)§4-1子程序PFORM的功能子程序P50SUBROUTINEPFORM(F,LD,UL,XL,S,P,U,X,IX,D,ID,JP,NST,NEN1,ISW)

LOGICALASEMCOMMON/ASEM/ASEM

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL

DIMENSIONUL(NDF,*),U(*),F(NDF,*),S(NST,*),P(*),JP(*),+XL(NDM,*),X(NDM,*),D(20,*),IX(NEN1,*),ID(NDF,*),+LD(NDF,*)

IEL=0

MCT=0！設(shè)置打印輸出DO110N=1,NE

CALLPZERO(UL,NST*(NST+3))

MA=IX(NEN1,N)！材料號

DO108I=1,NEN

II=IX(I,N)

IF(II.NE.0)GOTO105SUBROUTINEPFORM(F,LD,UL,51DO102J=1,NDF

102LD(J,I)=0

GOTO108

105IID=II*NDF-NDF

NEL=I！確定單元的實際節(jié)點數(shù)DO106J1=1,NDM

106XL(J1,I)=X(J1,II)！確定單元的節(jié)點坐標DO107J=1,NDF

K=ID(J,II)！確定自由度標志數(shù)組IF(K.GT.0)UL(J,I)=U(K)！確定單元的活化自由度數(shù)INEN=I+NEN

IF(K.LT.0)UL(J,I)=U(NEQ-K)！確定單元的約束位移IF(K.LT.0)UL(J,INEN)=F(J,II)！確定單元的指定位移

IF(ISW.EQ.6)K=IID+J

107LD(J,I)=K！確定自由度標志數(shù)組108CONTINUEDO102J=1,NDF

102LD(J,I52IE20=D(20,MA)

IF(IE20.NE.IEL)MCT=0！設(shè)置打印輸出

IEL=IE20！材料號CALLELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)！確定聯(lián)系信息IX

IF(ASEM.AND.ISW.EQ.3)CALLADDSTF(LD,JP,S,NST)

130IF(ISW.EQ.6)CALLBASBLY(F,P,LD,NST)

120CONTINUE

IF(.NOT.ASEM.AND.ISW.EQ.3)CALLMODIFY(U,LD,S,UL(1,NEN+1),NST)

IF(.NOT.ASEM.AND.ISW.EQ.3)CALLBASBLY(U,P,LD,NST)

109CONTINUE

110CONTINUE

END

UL：單元位移S：單剛P：單元內(nèi)力IE20=D(20,MA)

IF(I53§4-2子程序MODIFY完成將指定位移轉(zhuǎn)化成等價內(nèi)力按活化自由度和被約束自由度，將總剛方程分塊：u1

：活化自由度的位移u2

：被約束自由度的位移（包括指定位移）F1

已知F2

：未知§4-2子程序MODIFY完成將指定位移轉(zhuǎn)化成等價內(nèi)力54

SUBROUTINEMODIFY(U,LD,S,DUL,NS)REALU,DUL,SDIMENSIONU(1),LD(1),S(NS,1),DUL(1)DO110I=1,NSK=LD(I)IF(K.LE.0)GOTO110DO100J=1,NS100U(K)=U(K)-S(I,J)*DUL(J)110CONTINUERETURNENDSUBROUTINEMODIFY(U,LD55§4-3將單剛組集成總剛子程序ADDSTFCALLELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST,ISW)PFOR35構(gòu)造了單剛，存放于S中，ADDSTF完成組集PFORM中，LD(NDF,NEN)ADDSTF中,LD(NST)一維數(shù)組與二維滿存數(shù)組的元素對應(yīng)關(guān)系：IJ=JP(I)–I+JIJ§4-3將單剛組集成總剛子程序ADDSTFCALLEL56SUBROUTINEADDSTF(LD,JP,S,NST)CHARACTERY……DIMENSIONLD(*),JP(*),S(NST,NST)IF(Y.EQ.'Y'.OR.Y.EQ.'y')WRITE(*,2000)N,S2000FORMAT('ELEMENTNUMBER',I4,'ELEMENTSTIFFNESS'/(6E12.5))DO200J=1,NSTK=LD(J)IF(K.LE.0)GOTO200L=JP(K)-KDO100I=1,NSTM=LD(I)IF(M.GT.K.OR.M.LE.0)GOTO100M=L+MA(M)=A(M)+S(J,I)100CONTINUE200CONTINUEENDSUBROUTINEADDSTF(LD,JP,57§4-4子程序BASBLY組集載荷與節(jié)點反力向量單元自由度指示數(shù)組LDISW=3組集載荷向量，按節(jié)點活化自由度次序排列ISW=6節(jié)點反力向量，按節(jié)點排列次序§4-4子程序BASBLY組集載荷與節(jié)點反力向量單元自由58130IF(ISW.EQ.6)CALLBASBLY(F,P,LD,NST)PFOR37IF(.NOT.ASEM.AND.ISW.EQ.3)CALLBASBLY(U,P,LD,NST)PFOR40SUBROUTINEBASBLY(B,P,LD,NS)BASB1DIMENSIONB(*),P(*),LD(*)BASB2DO100I=1,NSBASB4II=LD(I)BASB5IF(II.GT.0)B(II)=B(II)+P(I)BASB6100CONTINUEBASB7RETURNBASB8ENDBASB9130IF(ISW.EQ.6)CALLBASBLY59§4-5PRTREAC輸出節(jié)點反力節(jié)點反力兩種計算方式：1）由總剛方程求得2）分別計算每個單元內(nèi)力（單剛乘以單元節(jié)點位移），然后組集成總體內(nèi)力向量因總剛分解后，已不存在了，故選2）§4-5PRTREAC輸出節(jié)點反力節(jié)點反力兩種計算方式60SUBROUTINEPRTREAC(R)

COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ

REALR(NDF,*),RSUM(6)

NNEQ=NDF*NJ

DO50K=1,NDF

RSUM(K)=0.0

50CONTINUE

DO100N=1,NJ,50

J=MIN0(NJ,N+49)

WRITE(6,2000)(K,K=1,NDF)

DO100I=N,J

DO75K=1,NDF

R(K,I)=-R(K,I)

RSUM(K)=RSUM(K)+R(K,I)

75CONTINUESUBROUTINEPRTREAC(R)

61DO76K=1,NDF76IF(ABS(R(K,I)).GT.1.0E-3)GOTO77GOTO10077WRITE(6,2001)I,(R(K,I),K=1,NDF)100CONTINUEWRITE(6,2002)(RSUM(K),K=1,NDF)RETURNFORMAT(/5X,'NODALREACTIONS'/2X,'NODE',6(I8,'DOF'))2001FORMAT(I6,6F12.4)2002FORMAT(3X,'SUM',6F12.4)ENDDO76K=1,NDF62§4-6PRTDIS輸出節(jié)點位移節(jié)點位移按節(jié)點次序重新排列注意：1）解得得位移向量為活化自由度，編號小于NEQ2）約束(包括指定位移)得值存放在二維數(shù)組F（載荷）中§4-6PRTDIS輸出節(jié)點位移節(jié)點位移按節(jié)點次序重新63SUBROUTINEPRTDIS(ID,U,F)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQDIMENSIONID(NDF,*),U(*),F(NDF,*),UL(6)DO102II=1,NJ,50WRITE(6,2000)(K,K=1,NDF)JJ=MIN0(NJ,II+49)DO102N=II,JJDO100I=1,NDFK=ID(I,N)IF(K.LT.0)U(NEQ-K)=F(I,N)IF(K.LT.0)UL(I)=U(NEQ-K)100IF(K.GT.0)UL(I)=U(K)WRITE(6,2001)N,(UL(I),I=1,NDF)102CONTINUE2000FORMAT(//1X,'NODALDISPLACEMENTS'/2X,'NODE',6(I8,'DISP'))2001FORMAT(I6,6E12.4)ENDSUBROUTINEPRTDIS(ID,U,F64§4-7ELEMLIB單元庫管理程序FEP2中已有兩種單元：BEAM、PLANE若要增加單元，可修改第七句。例如，已編好一個板單元，名為：PLATE（形參），則程序修改如下：……GOTO(1,2,3)IEL……CALLPLATE（D,U,X,IX,S,P,I,J,K,ISW)RETURN……§4-7ELEMLIB單元庫管理程序FEP2中已有兩種65SUBROUTINEELEMLIB(D,U,X,IX,S,P,I,J,K,ISW)REALP(K),S(K,K),U(1)DIMENSIONIX(1),D(1),X(1)COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQCOMMON/ELDATA/LM,N,MA,MCT,IEL,NELIF(IEL.LE.0.OR.IEL.GT.2)GOTO400GOTO(1,2)IEL1CALLBEAM(D,U,X,IX,S,P,I,J,K,ISW)GOTO1002CALLPLANE(D,U,X,IX,S,P,I,J,K,ISW)GOTO100100RETURN400WRITE(*,4000)IELSTOPFORMAT(5X,'**FATALERROR**ELEMENTCLASSNUMBER',I5,'INPUT')ENDSUBROUTINEELEMLIB(D66第五章平面單元設(shè)計平面四節(jié)點單元的基本公式數(shù)值積分形函數(shù)及其導(dǎo)數(shù)四節(jié)點單剛的構(gòu)造3–9節(jié)點平面單元平面單元PLANE子程序第五章平面單元設(shè)計平面四節(jié)點單元的基本公式67§5-1平面四節(jié)點單元的基本公式在FEP2中，提供了3—9節(jié)點平面單元，其基本形式是四節(jié)點等參元（四邊形單元）。子程序：PLANE(D,UL,XL,IX,S,P,NDF,NDM,NST,ISW)

PLANE中包含的子程序PGAUSS(L,LINT,SG,TG,WG)SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.)SHAP2(SS,TT,SHP,IX,NEL)PSTRES(SIG,SIG(4),SIG(5),SIG(6))§5-1平面四節(jié)點單元的基本公式在FEP268§5-2數(shù)值積分

SUBROUTINEPGAUSS(L,LINT,R,Z,W)COMMON/ELDATA/LM,N,MA,MCT,IEL,NELDIMENSIONLR(9),LZ(9),LW(9),R(3),Z(3),W(3),G4(4),H4(4)DATALR/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/DATALW/4*25,4*40,64/LINT=L*LGOTO(1,2,3,4),L1R(1)=0.0E0Z(1)=0.0E0W(1)=4.0E0IF(NEL.EQ.3)THENZ(1)=-1.0E0/3.0E0W(1)=3.0E0ENDIFRETURN§5-2數(shù)值積分SUBROUTINEPG692G=1.0E0/SQRT(3.0E0)DO21I=1,4R(I)=G*LR(I)Z(I)=G*LZ(I)21W(I)=1.0E0RETURN3G=SQRT(0.6E0)H=1.0E0/81.0E0DO31I=1,9R(I)=G*LR(I)Z(I)=G*LZ(I)31W(I)=H*LW(I)RETURN2G=1.0E0/SQRT(3.0E0)704G=SQRT(4.8E0)H=SQRT(30.0E0)/36.0E0G4(1)=SQRT((3.0E0+G)/7.0E0)G4(4)=-G4(1)G4(2)=SQRT((3.0E0-G)/7.0E0)G4(3)=-G4(2)H4(1)=0.50E0–H；H4(2)=0.50E0+HH4(3)=0.50E0+H；H4(4)=0.50E0-HI=0DO41J=1,4DO41K=1,4I=I+1R(I)=G4(K)Z(I)=G4(J)W(I)=H4(J)*H4(K)41CONTINUEEND4G=SQRT(4.8E0)71§5-3形函數(shù)及其導(dǎo)數(shù)SUBROUTINESHAPE(SS,TT,X,SHP,XSJ,NDM,NEL,IX,FLG)DIMENSIONSHP(3,*),X(NDM,1),S(4),T(4),XS(2,2),SX(2,2),IX(*)DATAS/-0.5E0,0.5E0,0.5E0,-0.5E0/,T/-0.5E0,-0.5E0,0.5E0,0.5E0/DO100I=1,4SHP(3,I)=(0.5e0+S(I)*SS)*(0.5e0+T(I)*TT)SHP(1,I)=S(I)*(0.5e0+T(I)*TT)100SHP(2,I)=T(I)*(0.5e0+S(I)*SS)IF(NEL.GE.4)GOTO120C...FORMTRIANGLEBYADDINGTHIRDANDFOURTHTOGETHERDO110I=1,3110SHP(I,3)=SHP(I,3)+SHP(I,4)！3節(jié)點平面單元

C....ADDQUADRATICTERMSIFNECESSARY120IF(NEL.GT.4)CALLSHAP2(SS,TT,SHP,IX,NEL)§5-3形函數(shù)及其導(dǎo)數(shù)SUBROUTINESHAPE72C....CONSTRUCTJACOBIANANDITSINVERSEDO130I=1,NDMDO130J=1,2XS(I,J)=0.0DO130K=1,NEL130XS(I,J)=XS(I,J)+X(I,K)*SHP(J,K)XSJ=XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1)IF(FLG)RETURNSX(1,1)=XS(2,2)/XSJSX(2,2)=XS(1,1)/XSJSX(1,2)=-XS(1,2)/XSJSX(2,1)=-XS(2,1)/XSJC....FORMGLOBALDERIVATIVESDO140I=1,NELTP=SHP(1,I)*SX(1,1)+SHP(2,I)*SX(2,1)SHP(2,I)=SHP(1,I)*SX(1,2)+SHP(2,I)*SX(2,2)140SHP(1,I)=TPENDC....CONSTRUCTJACOBIANANDI73§5-4四節(jié)點單剛的構(gòu)造比較三種不同的算法，選運算量最少的一種。1）按矩陣直接運算2）預(yù)先算出DBj

，再與BiT相乘，展開3）按矩陣運算全部展開計算：BiT

DBj§5-4四節(jié)點單剛的構(gòu)造比較三種不同的算法，選運算量最少743L=D(5)IF(L*L.NE.LINT)CALLPGAUSS(L,LINT,SG,TG,WG)DO320L=1,LINTCALLSHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.)XSJ=XSJ*WG(L)*D(7)J1=1DO320J=1,NELW11=SHP(1,J)*XSJW12=SHP(2,J)*XSJK1=J1DO310K=J,NELS(J1,K1)=S(J1,K1)+W11*SHP(1,K)S(J1,K1+1)=S(J1,K1+1)+W11*SHP(2,K)S(J1+1,K1)=S(J1+1,K1)+W12*SHP(1,K)S(J1+1,K1+1)=S(J1+1,K1+1)+W12*SHP(2,K)310K1=K1+NDF320J1=J1+NDF3L=D(5)75§5-53–9節(jié)點平面單元

3節(jié)點平面單元SHAPE12~135–9節(jié)點平面單元§5-53–9節(jié)點平面單元3節(jié)點平面單元SHA76SUBROUTINESHAP2(S,T,SHP,IX,NEL)DIMENSIONIX(*),SHP(3,*)S2=(1.0E0-S*S)/2.0E0T2=(1.0E0-T*T)/2.0E0DO100I=5,9DO100J=1,3100SHP(J,I)=0.0IF(IX(5).EQ.0)GOTO101SHP(1,5)=-S*(1.0E0-T)SHP(2,5)=-S2SHP(3,5)=S2*(1.0E0-T)101IF(NEL.LT.6)GOTO107IF(IX(6).EQ.0)GOTO102SHP(1,6)=T2SHP(2,6)=-T*(1.0E0+S)SHP(3,6)=T2*(1.0E0+S)SUBROUTINESHAP2(S,T,SHP,IX,N77102IF(NEL.LT.7)GOTO107SHAP19IF(IX(7).EQ.0)GOTO103SHAP20SHP(1,7)=-S*(1.0E0+T)SHAP21SHP(2,7)=S2SHAP22SHP(3,7)=S2*(1.0E0+T)SHAP23103IF(NEL.LT.8)GOTO107SHAP24IF(IX(8).EQ.0)GOTO104SHAP25SHP(1,8)=-T2SHAP26SHP(2,8)=-T*(1.0E0-S)SHAP27SHP(3,8)=T2*(1.0E0-S)SHAP28C....INTERIORNODE(LAGRANGIAN)104IF(NEL.LT.9)GOTO107SHAP30IF(IX(9).EQ.0)GOTO107SHAP31SHP(1,9)=-4.0E0*S*T2SHAP32SHP(2,9)=-4.0E0*T*S2SHAP33SHP(3,9)=4.0E0*S2*T2102IF(NEL.LT.7)GOTO10778C....CORRECTEDGENODESFORINTERIORNODEDO106J=1,3SHAP36DO105I=1,4SHAP37105SHP(J,I)=SHP(J,I)-0.25E0*SHP(J,9)SHAP38DO106