計算結構力學PPT

1、第八章 結構分析程序設計81 概述:程序設計的基本概念與要點至此,我們已完成了結構矩陣分析桿系有限元法基本原理的教學內容,本章主要介紹計算機實現(xiàn)過程。程序設計:當人們需要計算機完成科學計算,數(shù)據(jù)處理等計算工作時,必須事先恰當?shù)匕才藕糜嬎悴襟E,每一步的計算內容以及在什么條件下進行這一部分計算,這就是所謂編制計算機程序。結構矩陣分析原理結構分析基礎數(shù)值計算方法程序設計算法語言計算結果分析上機調試程序設計應注意以下幾點：保證程序的正確性，可通過考題校驗使程序具有高效率,并盡可能降低費用,求解方程組最費機時(80%左右),應設計再解功能使程序便于調試、修改、擴充和完善,既要有通用性,又要留下可擴充修改

2、的窗口從結構矩陣分析原理到最終計算機實現(xiàn)解決具體問題,主要有以下過程模塊,用框圖可表示為：8.2 結構分析程序設計的框圖設計根據(jù)結構方程組的建立與求解來進行程序設計結構方程組的建立與求解是結構分析的核心問題,如建立結構方程組的有限元方法,解線性代數(shù)方程組的消元分解法,這兩個方法一經(jīng)確定,程序設計的大致思路與過程也就基本確定了。1. 編程要點由單元定位向量組織整體流程圖的運行實施 從形成結構方程組K=P過程中K與P的形成,到計算結構內力和反力的過程,都離不開MW的組織。應具有良好的通用性 本程序的編制特點是利用特殊結點的約束信息，可模擬六種類型的桿系結構，故具有廣泛的適用性。2. 設計流程圖主程

3、序,數(shù)據(jù)輸入形成總剛K對K消元分解形成荷載列陣P對P正消，回代求出結點位移求內力與反力 形成MW 形成單剛K1234567對各組數(shù)據(jù)循環(huán)對各組荷載循環(huán)3. 框圖設計主、數(shù)據(jù)輸入KJX:形成KNXFJ:消元分解YDX:形成 PNXFJ:回代求解YWY:打印位移QDL:求F與RQJW21222324252627 結束 開始 QJZHQCHQKDQXSDKXPGPDXL1112131415161718ZERO1ZERO2JZZQMW010203044. 本程序設計的模塊功能介紹模塊：具有單一的獨立的功能塊,由子程序或自定義函數(shù)所組成。模塊可根據(jù)其功能進一步劃分,依次分解成較低級的模塊,模塊之間通過調

4、用而組成一個協(xié)同的程序；這種可通過自上而下進行分解,并可通過自上而下的調用,一級一級地組成程序是程序設計的重要方法。各個模塊的功能在很大程度上是獨立的,因而不同的模塊可以由不同的人來完成。例如,不太了解矩陣力學的人也可以設計消元分解及正消回代的子程序SUB. NXFJ。模塊的相互獨立性不僅方便程序設計,也方便程序調試。調試時自下而上一塊一塊地進行。這時主要著重調試模塊接口和上一級模塊,而不必調試已通過的模塊。亦可以根據(jù)模塊功能組成其它程序。本程序的模塊設有三個級別:0、1、20級表示各個程序均可調用1級僅供2級調用。5. 本程序的靜力計算功能(1)結構形式:可對六種類型的桿系進行計算分析連續(xù)梁

5、桁架排架剛架框排架梁、桁組合結構?？捎蒍TX(4,NJT)進行分類。(2)材料:各向同性,按EI(EA)進行分組,分組數(shù)為NAI(3)荷載類型結點荷載信息NPJ:受載結點數(shù) 需輸入結點荷載信息數(shù)組PJZ(NPJ,2),NPJ行,2列:第一列輸JD.x,第二列輸大小(與坐標一致為正，反之為負)。JD是結點號，x是方向(1為X向, 2為Y向, 3為z向)。單元荷載信息單元荷載信息數(shù)組PMZ(NPM,3):NPM:受載單元數(shù)NPM行,3列第一列輸M.L:第M單元，第L類荷載；第二列輸荷載位置(距始端距離x)第三列輸荷載大小(與坐標一致為正，反之為負)L類：共六類,見講義附表(P63)或參考教材表5-

6、1(P131),程序PSTDY的子程序SUB.DXL中留有用戶入口,可接入這六種以外的荷載。83 單元定位向量的主線作用從程序設計框圖可以看出:程序設計的每一個環(huán)節(jié)都離不開單元定位向量,故它在程序設計中起到組織者的作用。我們稱之為主線作用?；仡?單元定位向量是按單元結點編號順序由結點的未知量編號所組成的一個列向量?？捎蒍W(3,NJ)直接生成,其作用主要有：1.裝備結構剛度矩陣（1）按方陣存貯SUB.KJX1DO 10 I=1,6L=MW(I)IF (L.LE.0) GOTO 10DO 20 J=1,6K=MW(J)IF (K.LE.0) GOTO 20ZK(L,K)=ZK(L,K)+DK(I

7、,J)20CONTINUE10CONTINUE思考題：若形成上三角陣如何改動（2）一維變帶寬上三角按行存貯SUB.KJXDO 10 I=1,6L=MW(I)IF (L.LE.0) GOTO 10II=KD(L)DO 20 J=1,6K=MW(J)IF (K.LE.L) GOTO 20IJ=II-L+KZK(IJ)=ZK(IJ)+DK(I,J)20CONTINUE10CONTINUE2.形成P見PSTDY中SUB.YDX3.單元位移的形成見PSTDY中SUB.YWY由(存在P中，存在D(NE,6)中4.FjRj亦由MW的第j個分量是否為零來判斷,見PSTDY中SUB.QDL5.在KD數(shù)組中求帶寬

8、亦用到MW見PSTDY中SUB.QKD帶寬公式：NDK=單元兩端未知量編號最大差值+16.更為重要的是單元定位向量還體現(xiàn)了單元間的相互聯(lián)接，以及對結構邊界條件的處理，如主從關系、無效未知量的處理等，可通過對特殊結點的約束信息數(shù)組JTX(4,NJT)來模擬實際結構中的復雜關系。換句話說，我們在程序的功能中所提到的六種類型的桿系結構，程序設計中最后的區(qū)別形式就是MW，這點尤為重要。對各種類型的結構，我們設計了統(tǒng)一的單元剛度矩陣形式，由MW來直接裝配總剛！六種類型桿系結構的JTX數(shù)組1.連續(xù)梁2.桁架3.排架4.框架5.框排架6.框排架根據(jù)各程序模塊的功能，本程序設計了下列控制變量與循環(huán)變量： NE

9、單元總數(shù) NJ結點總數(shù) NJT特殊結點數(shù) NJZ支座結點數(shù) NAIEA或EI的分組數(shù) NPJ結點荷載數(shù) NPM單元荷載數(shù)84 變量與數(shù)組設計 N未知量總數(shù) M單元序號 LOAD荷載分組序號 I、J單元剛度矩陣的行列號 L，K結構剛度矩陣的行列號 II，LL對角線元素地址 NI-KD數(shù)組的個數(shù)， NZY結構剛度矩陣的元素數(shù)整型數(shù)組有 JH(2，NE)單元兩端的結點號 MW(6)單元定位向量 JW(3，NJ)結構結點未知量編號 JTX(4，NJT)特殊結點信息， JZH(NJZ)支座維點號 NLX(2，NL)每組荷載信息，即NPJ數(shù)和NPM數(shù) NLX(2，NL)截面特性分組號 KD(NI)結構剛陣

10、一維存貯時主元素地址；雙精度實型數(shù)組 X(NJ)結點坐標值 Y(NJ)結點坐標值 SL(NE)單元長度 CX(NE)單元的cos SY(NE)單元的sin EA(NAI)單元的EA El(NAI)單元的EI XS(7)整體坐標系下各單元的7個常數(shù)， XSA(NE，7)整體坐標系下各單元的7個常數(shù) DK(6，6)單元剛度矩陣， T(6,6)坐標變換矩陣， ZK(NZY)結構剛度矩陣，一維存貯 P(N)荷載列陣 D(NE，6)單元兩端的結點位移，后存結點力 F(6)單元結點力 FE(6)等效結點力 FG(6)整體坐標系下的單元結點力 DG(6)整體坐標系下單元結點位移 FLZ(NJZ，3)支座反力

11、 PMZ(NPM，3)單元荷載信息 PJZ(NPJ，2)結點荷載信息 在采用FORTRAN語言編制程序時，一定要摘清楚數(shù)據(jù)如何傳遞。鑒于FORTRAN語言的模塊化性質，各程序的數(shù)據(jù)一般可通過下列三種方式進行傳遞：啞實結合 COMMON塊 數(shù)據(jù)文件方法還可運用在機器設備與外部設備(終端)的數(shù)據(jù)傳遞，如本程序設計就采用2個OPEN語句，建立了輸入、輸出數(shù)據(jù)文件。我們在這里介紹的程序設計各子程序間的數(shù)據(jù)傳遞均采用方法 。85 數(shù)據(jù)傳遞與動態(tài)數(shù)組設計子程序的一般形式為SUBROUTINE QJW(NJ，NJT，JTX，JW，N)RETURNEND其中括號內的NJ，N即為形式參數(shù)或稱為虛擬變量，或稱啞元

12、。它可以是變量名字、數(shù)組名字或數(shù)組元素。1、啞實結合的數(shù)據(jù)傳遞方式子程序中的形式參數(shù)沒有確切的數(shù)值，這就是虛擬變量的由來。只有在調用該子程序時，才對形式參數(shù)賦值，或賦予其實在的存貯空間。如：CALL QJW(NJ，NJT，JTX，JW，N)這時NJ，N稱為實在參數(shù)(實元)。注意：形式參數(shù)和實在參數(shù)的類型應一致，個數(shù)應相等，但參數(shù)名可以不相同。形式參數(shù)的作用可分為兩種：一種是從主程序或其它子程序來賦值的，也就是通過啞實結合接收從外面輸入的數(shù)據(jù)，作為本子程序計算的依據(jù)，好象是加工廠的原料一樣另一種是本程序模塊計算的結果，通過啞實結合傳遞到調用處，這是向外傳遞的數(shù)據(jù)，好象是加工廠的產(chǎn)品一樣。如：子程

13、序：SUB. QJW(NJ，NJT，JTX，JW，N)主程序：CALL QJW(NJ，NJT，JTX，JW，N)動態(tài)數(shù)組又稱為可調數(shù)組，如：JH(2,NE)、JW(3,NJ)JTX(4,NJT)、ZK(NZY)等，由于這里NE、NJ、NJT，NZY等均為數(shù)組變量，沒有確切的數(shù)字，即數(shù)組的大小未能得到確切的定義，因而可調數(shù)組在主程序中是不允許出現(xiàn)的。在主程序中只能出現(xiàn)確切定義的數(shù)組，如我們在前二個大作業(yè)里所介紹的數(shù)組如JH(2,20)、ZK(50,50)等，這樣才能在在DIMENSION語句中予以確切定義。 2、主程序中動態(tài)數(shù)組的設計根據(jù)子程序的形式參數(shù)(啞元)的定義，在子程序中可以出現(xiàn)動態(tài)數(shù)組

14、，但必須在調用時進行啞實結合，才能進行運算。由于在實際計算中大多數(shù)數(shù)組的大小是隨具體問題的不同而變化的，但在主程序中又不允許出現(xiàn)動態(tài)數(shù)組，這就給我們在DIMENSION語句中如何定義數(shù)組帶來不少困難：既受到計算機內存的限制又應使數(shù)組有足夠的存貯空間，這對于在微機上解決大型工程問題顯得更重要。我們知道，二維或高維數(shù)組在計算機內部都是按列存放的，即在計算機內部都是按一維數(shù)組 的方式來存貯的，這表示數(shù)組變量之間存在一定的關系，如對于：JH2*20數(shù)組：JH(2，3)JH(6)，JH(1，4)JH(7)，這說明數(shù)組一經(jīng)定義，這種關系便確定，計算機立即“了解”。此外，一個數(shù)組變量(下標變量)的下標值加1

15、，就是緊跟在它后面的下標值，這叫做數(shù)組變量(下標變量)的后繼函數(shù)，計算機的處理功能使程序會自動按后繼函數(shù)找到下一個元素，這叫做下標的自動后繼性質。如在主程序DIMENSION語句中定義說明了JH(2，20)，則在調用時實元用JH或JH(1，1)或JH(1)調用的效果是相同的，如在主程序中用一維定義說明了JH(40)，則用JH或JH(1)調用效果亦是相同。至于調用的JH具體形式，則由子程序對可調數(shù)組JH(2，NE)進行說明確定，換句話說，利用下標的自動后繼性質，采用首元素調用，即可在主程序中出現(xiàn)動態(tài)數(shù)組。根據(jù)算例的需要，我們在主程序中僅開設了整、實二個大數(shù)組： 整：IA(1000)，暫定其大小N

16、IA=1000 實：A(10000)，暫定其大小NA=10000然后按照程序中所要出現(xiàn)數(shù)組名字順序，定出每個數(shù)組的第個元素地址，由所輸入的變量或巳確定的變量確定這一點并不難，參考下圖：整：實：JH(2，NE)的第一個元素地址為IA(1)NLX(2，NL)的第一個元素地址為IA(K1)，這里(K1)=1+2*NEJW(3，NJ)的第一個元素地址為IA(K2)，這里(K2)=K1+2*NL對于各雙精度實型數(shù)組在A中第一個元素地址亦可由圖方便推出。如將某數(shù)組的第一個元素作為實元調用，通過啞實結合過程中的下標自動后繼，當子程序啞元表中相應的啞元為動態(tài)數(shù)組時，便可完全按照子程序中該數(shù)組的DIMENSIO

17、N語句的定義在主程序相應數(shù)組中得到反映，即在主程序中實現(xiàn)了動態(tài)數(shù)組。由此可知，只要找出控制各數(shù)組大小的一些變量，即可確定各數(shù)組的第一個元素。仔細研究這些變量，發(fā)現(xiàn)其中某些變量(如N，N1，NZY)等可通過另一些變量由程序計算確定，我們將后者稱為主控變量，須按其出現(xiàn)的次序在程序中首先輸入。本程序的主控變量為：NE、NJ、NJT、NJZ、NL、NAI。為促進同學們今后的工程應用能力，我們這里所提供的程序PSTDY并有動力分析內容，故主控變是還增加一個MJ。MJ：擬求振型數(shù)。 在PSTDY中，我們將MJ采用屏幕輸入，并兼作靜、動力分析的程序運行控制，在屏幕的提示下，若鍵人 MJ=0僅作靜力計算 MJ

18、=3做完靜力計算后，再求前三階頻率和振型，具體請參閱算例1，不另贅述。 另外，在主程序中還備有這兩大類型數(shù)組實際使用的元素個數(shù)顯示，如發(fā)生IA或A的溢出現(xiàn)象，只需在內存允許的范圍內修改這四條語句，見主程序的001、002、004、005句。1、根據(jù)設計框圖，各子程序及其功能匯總如下：22 KJX：形成結構剛陣ZK(NZY)23、25 NXFJ：消元分解法解線代方程組23：剛陣ZK的消元25：對P的正消回代24 YDX：形成P26 YWY：打印桿端位移27 QDL：計算F及R86 源程序設計28(補充) MJX：形成結構質量的矩陣29(補充) DYNA：逆迭代法計算結構前MJ階振型和頻率11 Q

19、JW：結點未知量編號數(shù)組12 QJZH：形成支座結點號13 DCH：形成單元常數(shù)14 QKD：形成主元地址數(shù)組15 QXS：形成單剛系數(shù)數(shù)組16 DKX：形成單剛矩陣17 PGP：疊加荷載，形成P的具體計算，18 DXL：形成FE01 ZER01：對向量充零02 ZER02：對二維數(shù)組充零03 JZZ：矩陣轉置04 QMW：形成單元定位向量。這里的順序由二位數(shù)表示，十位數(shù)上有0，1，2級，其中高級別的可調用低級別的模塊，個位數(shù)則表示被調用的順序。由 OPEN(1，F(xiàn)ILE=QAZ.TXT)建立數(shù)據(jù)文件 READ(1，*) NE，NJ，NJT,NJZ，NL，NAI(見主程序語句標號003) RE

20、AD(1，*) JH，NLX，JMH，JTX(見主程序語句標號006)READ(1，*) EA，EI，X，Y(見程序語句標號007)READ(1，*) PJZ(NPJ，2)(見SUBYDX語句標號008) READ(1，*) PMZ(NPM，3)(見SUBYDX浯句標號009)2、數(shù)據(jù)文件的形成在PSTDY程序中，我們用OPEN語句建立了數(shù)據(jù)輸入文件QAETXT及結果輸出文件 FCADTXT。全部的原始數(shù)據(jù)，除擬求振型數(shù)MJ兼作靜、動力計算的運行控制從屏幕輸入外，均可在QAETXT文件上從終端讀人。全部的數(shù)據(jù)均采用自由格式，由READ(1，*)語句讀人，這里“1”表示通道號，“*”表示自由格式

21、。自由格式要求兩個數(shù)之間用“，分隔，且每個數(shù)所占的位數(shù)不限，但要與對應變量的隱含數(shù)型一致。3. 算例1數(shù)據(jù)文件可在每行開始直接輸入數(shù)據(jù)，不用說明符號“C”，也沒有語句標號區(qū)及續(xù)行區(qū)，但每行不應超過72列2每個輸入語句的結尾不要加任何標點符號，但實型數(shù)據(jù)的小數(shù)點除外3下一個語句不要與上一個輸入語句接著輸入，而應別起一行建立數(shù)據(jù)文件時還要注意以下幾點算例1作圖示剛架的靜力計算，并求出其基本頻率與振型。任取l=4m,m=3.0*102kg/mP=2KN,E=2.1*106KN/m2I=4*105cm4=4*10-3m4EI=8.4*103KNm2則各結點標號與單元劃分如圖示。1第一個輸入語句在主程序

22、003句，要求輸入控制變量6個： NE， NJ， NJT， NJZ， NL， NAI 3， 4， 4， 2， 1， 2注意，這六個整形數(shù)成為一行，句尾不加任何符號。2第二個輸入語句在主程序006句，要求輸入整型數(shù)組IA，實際上輸入JH，NLX，JMH，JTX，其中JH有2*NE=6個數(shù)，NLX有2*NL=2個數(shù)，JMHE有NE=3個數(shù)，JTX有4*NJT=16個數(shù)，總共29個整型數(shù)，可分為兩行，但這兩行之間要用逗號“，”隔開。3,1,4,2,1,2,1,0,1,1,2,1,0,10001,0,2,1001,10001,0,0,3,1,1,1,4,1,1,1 注意：這里特殊節(jié)點約束信息JTX(4

23、，NJT)的填寫是一項非常細致的工作，需根據(jù)具體力學模型，參照17中的有關說明，認真填寫。3第3個輸入語句在主程序007句，要求輸入雙精度型數(shù)組A，實際上輸入EA，EI，X，y。亦采用自由格式，EA有2個數(shù)，EI也有2個數(shù)，X有4個數(shù)，Y也有4個數(shù)，總共有12個數(shù)，可一行輸入：0,0,8400,33600,0,8,0,8,4,4,0,04在主程序對荷載的循環(huán)中，要輸入每一組荷載的數(shù)據(jù)，PSTDY程序仍采用通常的荷載信息集約方法，在集中荷載作用時，要給出節(jié)點荷載的個數(shù)，荷載作用在哪號節(jié)點的哪個位移方向以及荷載的大小，這些信息可以定義一個數(shù)組PJZ(NPJ，2)來存放。PJZ(I，1)可填實型數(shù)J

24、D.x其中：JD一荷載作用的節(jié)點號； x一荷載的作用方向，可在總體坐標系下考慮。x=1，荷載作用沿X方向；x=2，荷載作用沿Y方向,x=3，繞Z軸正向作用的力矩。 PJZ(1，2)填荷載的大小，與坐標方向一致者為正，I表示節(jié)點荷載序號。 在單元荷載作用下，PSTDY的靜力計算可解決六種類型的單元荷載問題，我們只要給出如下信息： (1)在哪個單元M上有什么類型L的荷載作用；(2)在右手系下荷載離單元左節(jié)點的距離；(3)荷載大小Q，其符號規(guī)定與坐標方向一致為正。于是也可定義一個數(shù)組PMZ(NPM，3)來存放上述信息。 PMZ(I，1)填M.L； PMZ(I，2)填x； PMZ(I，3)填Q；其中，

25、L=1，表示左端有局部均布荷載Q作用；L=2，表示離左端x處有集中荷載Q作用；L=3，表示離左端x處有集中力矩Q作用；L=4，表示左端有局部三角形荷載作用，其x處為QL=5，表示左端有均布軸力荷載Q作用；L=6，表示離左端x處有集中軸力Q作用；在子程序YDX中，有兩個輸入語句，應分別輸入PJZ與PMZ。本例只有一組結點荷載NL=1，沒有單元荷載NPM=0，結點荷載數(shù)NPJ=1，所以第四個輸入語句在子程序SUBYDX中標號為008句，輸入結點荷載信息PJZ，共2*NPJ=2個數(shù)，也占行： 1.1，2000注意，由于NPM=0，單元荷載輸入信息PME(NPM，3)的輸入語句009，程序自動跳轉，該

26、語句可以不輸入。5至此，靜力計算的所有信息輸入完畢，如不做動力計算，可直接從屏幕上輸“0”，程序運行結束，并輸出結果。如需要進行動力計算時，應輸入所求振型數(shù)MJ，從屏幕上輸入。本例MJ=1，可直接輸入“1”。當MJ0時，還應在QAZTXT文件上繼續(xù)輸入動力計算的有關數(shù)據(jù)，這樣，第五個輸入語句在子程序SUBMJX的010句，輸入ERM(NAI)桿件質量線密度數(shù)據(jù)，本題有NAI=2個數(shù)，也占一行： 300，450解 根據(jù)圖55，按程序中輸入語句的順序依次填寫數(shù)據(jù)。由算例，可知： 1第一個輸入語句標號003，要求輸人控制變量6個：NE，NJ，NJT，NJZ，NL，NAI 9， 8， 2， 2， 1， 2這6個數(shù)恰好一行，在文件中叫做一個記錄。 算例2 做54節(jié)例3三層剛架的計算。2第二個輸入語句標號006，要求輸入整型數(shù)組IA，實際輸入JH，NLX，JMH，J