1、有限元講義 1 2.3 平面問(wèn)題有限元程序設(shè)計(jì) 2.3 平面問(wèn)題有限元程序設(shè)計(jì) 一、程序設(shè)計(jì)方法與結(jié)構(gòu)分析程序的特點(diǎn) 一、程序設(shè)計(jì)方法與結(jié)構(gòu)分析程序的特點(diǎn) 1程序設(shè)計(jì)方法論簡(jiǎn)述程序設(shè)計(jì)方法論簡(jiǎn)述 借助計(jì)算機(jī)來(lái)完成某項(xiàng)工作，通常都要先編寫相應(yīng)的計(jì)算機(jī)程序，或叫程序設(shè)計(jì)。完 成一個(gè)結(jié)構(gòu)分析或結(jié)構(gòu)CAD系統(tǒng)也必然要經(jīng)過(guò)程序設(shè)計(jì)才能實(shí)現(xiàn)。 程序設(shè)計(jì)要使用專門的程序語(yǔ)言。我國(guó)結(jié)構(gòu)程序設(shè)計(jì)中所采用的語(yǔ)言，在60年代和70 年代初以ALGOL語(yǔ)言ALGOL語(yǔ)言為主。此后逐步廣泛使用的主要是BASICBASIC語(yǔ)言和FORTRANFORTRAN語(yǔ)言，隨著CAD 和人工智能技術(shù)的發(fā)展，PASCAL、 C、LIS

2、P、 PROLOG等有著各自特長(zhǎng)的程序語(yǔ)言也逐步 進(jìn)入土木工程領(lǐng)域的計(jì)算機(jī)程序設(shè)計(jì)中。 過(guò)去人們通常認(rèn)為，程序設(shè)計(jì)的中心問(wèn)題就是學(xué)會(huì)使用一種程序語(yǔ)言，用以編寫程序。然 而學(xué)會(huì)用程序語(yǔ)言編程只是整個(gè)程序設(shè)計(jì)中的一部分。據(jù)有關(guān)資料介紹，編寫程序在整個(gè) 系統(tǒng)的研制過(guò)程中僅占1515的工作量。在一個(gè)大型程序系統(tǒng)的整個(gè)存在階段的工作量中， 在系統(tǒng)投入使用后的維護(hù)工作量為原來(lái)研制工作量總和的兩倍（這一點(diǎn)在作者所從事的軟 件開(kāi)發(fā)工作中也得到充分的證明）。 維護(hù)工作量是如此之高，這就使我們必須注意到，在程序研制階段便即應(yīng)當(dāng)考慮為以后的 維護(hù)工作提供方便，哪怕是為此要增加一些額外的工作量也是值得的。 要編制一個(gè)

3、好的程序系統(tǒng)并沒(méi)有一種絕對(duì)的規(guī)則，就象是工程設(shè)計(jì)沒(méi)有一種絕對(duì)規(guī)則 一樣。但對(duì)于程序設(shè)計(jì)的好壞現(xiàn)在已逐漸形成了一套評(píng)價(jià)的客觀標(biāo)準(zhǔn)。這些標(biāo)準(zhǔn)大致分為 以下幾個(gè)主要方面： (1) 程序的可讀性; (2) 正確性與可靠性; (3) 使用方便且效率高; (4) 軟件的可移置性; (5) 易于調(diào)試與維護(hù)。 直到1970年代中期人們才認(rèn)識(shí)到軟件的維護(hù)是軟件研究的一個(gè)關(guān)鍵領(lǐng)域。造成軟件維 護(hù)工作量大的原因之一是與程序研制過(guò)程中所采用的設(shè)計(jì)方法不夠科學(xué)化有關(guān)。為了解決 這一問(wèn)題，人們開(kāi)展了對(duì)于程序設(shè)計(jì)方法論的研究與實(shí)踐，其目標(biāo)是使軟件正確、可靠和 降低整個(gè)軟件研制活動(dòng)的費(fèi)用?？偟膩?lái)說(shuō)，程序設(shè)計(jì)已從強(qiáng)調(diào)靈活的技巧

4、和局部效率向著 強(qiáng)調(diào)程序結(jié)構(gòu)化和整體功能的方向發(fā)展。這實(shí)際上就是逐步發(fā)展起來(lái)的關(guān)于程序設(shè)計(jì)編寫 與調(diào)試的一套方法論，其要點(diǎn)可歸納為以下幾方面： (1) 編程結(jié)構(gòu)化。為了使程序設(shè)計(jì)者能按照一定的結(jié)構(gòu)形式，而不是隨心所欲地設(shè)計(jì) 和編寫程序，使編制的程序易讀、易修改，以提高程序設(shè)計(jì)和維護(hù)工作的效率，荷蘭學(xué)者 Dijkstra提出了“結(jié)構(gòu)化程序設(shè)計(jì)方法”。結(jié)構(gòu)化程序規(guī)定了三種基本的結(jié)構(gòu)形式，他們 是順序結(jié)構(gòu)、分支選擇結(jié)構(gòu)和循環(huán)結(jié)構(gòu)。編程結(jié)構(gòu)化又稱結(jié)構(gòu)化程序設(shè)計(jì)，他可使編寫的 程序?qū)哟畏置?，邏輯清楚，容易閱讀。 (2) 分層處理技術(shù)。為了解決現(xiàn)實(shí)世界中的許多復(fù)雜問(wèn)題，人們往往需要根據(jù)問(wèn)題的 內(nèi)在聯(lián)系將其

5、分割成有層次的一系列問(wèn)題來(lái)分別求解。對(duì)于一個(gè)大型程序系統(tǒng)設(shè)計(jì)來(lái)說(shuō)， 也需采用分層的辦法來(lái)處理， 在每一層里集中解決一個(gè)問(wèn)題， 并為下一層的執(zhí)行作好準(zhǔn)備。 分層處理技術(shù)的主要內(nèi)容是將程序劃分為多個(gè)層次的若干模塊，每個(gè)模塊完成一個(gè)（或幾 個(gè)）預(yù)定的功能。 為了保證模塊的獨(dú)立性，各模塊之間只能通過(guò)接口與其他模塊連接。另外，對(duì)于一個(gè) 較大的軟件系統(tǒng)要由多人合作才能完成，模塊化也為此提供較好的合作條件。模塊化的一 種典型結(jié)構(gòu)形式就是子程序結(jié)構(gòu)。 (3) 避免過(guò)多使用GOTO語(yǔ)句，特別是逆轉(zhuǎn)的GOTO語(yǔ)句。 這是結(jié)構(gòu)化程序設(shè)計(jì)的基本要求之一。 有限元講義 2 （4）可移植性 對(duì)于應(yīng)用軟件， 特別是大型的C

6、AD軟件， 可移植性高低同樣是衡量軟件的質(zhì)量的重要指標(biāo)。 因?yàn)檠兄埔粋€(gè)大型應(yīng)用軟件系統(tǒng)不但要耗費(fèi)大量的人力物力，而且還要花費(fèi)相當(dāng)長(zhǎng)的時(shí)間 才能研制成功。在其研制和應(yīng)用期間，不可避免地發(fā)生運(yùn)行環(huán)境的變化計(jì)算機(jī)硬件設(shè) 備的換代和系統(tǒng)軟件的更新?？梢浦残詮?qiáng)即等于延長(zhǎng)了軟件的生存期，從而可節(jié)省新的開(kāi) 發(fā)投資。可移植性主要體現(xiàn)在軟件對(duì)支撐環(huán)境的獨(dú)立性和軟件本身的封閉性。例如，用 FORTRAN、BASIC 、C等語(yǔ)言編寫，并盡量避免采用非標(biāo)準(zhǔn)語(yǔ)句和函數(shù)。此外，在軟件中采 用統(tǒng)一的I/O 模塊，也是提高可移植性的手段之一。 應(yīng)當(dāng)指出，程序設(shè)計(jì)方法論仍在發(fā)展探索之中，千萬(wàn)不能把上述有關(guān)內(nèi)容當(dāng)成一成不 變的教

7、條套用，而應(yīng)當(dāng)通過(guò)實(shí)踐來(lái)發(fā)展和豐富其內(nèi)容。然而程序設(shè)計(jì)發(fā)展到今天，已經(jīng)奠 定了很多必要的理論基礎(chǔ)。我們正在達(dá)到一個(gè)可以談?wù)摮绦蛟O(shè)計(jì)是一門科學(xué)而不僅僅是一 種技巧的階段。 2、結(jié)構(gòu)化程序設(shè)計(jì) 2、結(jié)構(gòu)化程序設(shè)計(jì) 結(jié)構(gòu)化程序設(shè)計(jì)，又稱結(jié)構(gòu)程序設(shè)計(jì)(Structured Programming)是 荷蘭學(xué)者E. W. dijkstra首先提出來(lái)的，簡(jiǎn)稱SP。 人們對(duì)SP有各種定義和解釋: 有人說(shuō)它是: 指導(dǎo)程序員編程的一般方法； 有人說(shuō)它是: 不使用goto語(yǔ)句的程序設(shè)計(jì); 有人說(shuō)它是: 自頂向下的程序設(shè)計(jì)。 有人把層次結(jié)構(gòu)、順序結(jié)構(gòu)、選擇結(jié)構(gòu)和重復(fù)結(jié)構(gòu)定義為SP的全部?jī)?nèi)容，而把一個(gè)程 序結(jié)構(gòu)模塊定

8、義為上述四種基礎(chǔ)結(jié)構(gòu)的某種集合，即： 結(jié)構(gòu)模塊(程序)(層次結(jié)構(gòu), 順序結(jié)構(gòu), 選擇結(jié)構(gòu), 重復(fù)結(jié)構(gòu)) 認(rèn)為SP就是定義為這些基本結(jié)構(gòu)在軟件開(kāi)發(fā)和維護(hù)中的嚴(yán)格應(yīng)用。 3、程序的正確性及其驗(yàn)證 3、程序的正確性及其驗(yàn)證 程序正確性(Program Correctness)被定義為一個(gè)程序和它打算要實(shí)現(xiàn)的功能之間的 一次符合(a correspondness) Gries說(shuō)，在程序設(shè)計(jì)的初期階段，人們很少看到程序正確性方面的問(wèn)題，那時(shí)人們 往往著重于調(diào)試(debugging),但是后來(lái)發(fā)現(xiàn)調(diào)試好的程序并不能表示錯(cuò)誤不存在。于是人 們就不得不考慮程序的正確性證明問(wèn)題。 什么是程序正確性證明呢？從詞義

9、上講，“證明”就是提供一種有力的證據(jù)使人們?cè)?思想上不得不接受一條真理或一件事實(shí)。但是在程序設(shè)計(jì)中要做到這一點(diǎn)并不容易，因?yàn)?現(xiàn)在還沒(méi)有找到一種象數(shù)學(xué)一樣非常完整的形式。目前主要還是通過(guò)選擇較好的算法，大 量考題，與已有結(jié)果比較等辦法來(lái)證明一個(gè)程序的正確性。但不管怎樣，程序正確性理論 已經(jīng)被提出來(lái)，并有了一些實(shí)際想法。 程序的驗(yàn)證(Program Verification)實(shí)際上就是檢驗(yàn)程序的正確性。 一個(gè)程序如果有錯(cuò)誤， 主要是兩方面的：一方面是語(yǔ)法錯(cuò)誤，這部分比較好解決，一般是在調(diào)試階段（編輯階段） 完成。但是一個(gè)語(yǔ)法完全沒(méi)問(wèn)題的程序并不一定是正確的。因?yàn)槌绦蛑械脑S多部分往往是 靠邏輯關(guān)系

10、來(lái)達(dá)到所要實(shí)現(xiàn)的目的。目前應(yīng)用程序的驗(yàn)證主要靠針對(duì)程序每一功能，每一 邏輯分支進(jìn)行各種類型的考題，包括考題的規(guī)模。 4、結(jié)構(gòu)分析程序的特點(diǎn) 4、結(jié)構(gòu)分析程序的特點(diǎn) 1. 規(guī)律性、 通用性好。 有限元講義 2. 計(jì)算工作量大、運(yùn)算時(shí)間長(zhǎng)（形成單、總剛、解方程、動(dòng)力、非線性）。 3. 矩陣階數(shù)高，與微機(jī)的存貯量發(fā)生矛盾。 二、平面問(wèn)題有限元程序設(shè)計(jì)（三角形單元） 二、平面問(wèn)題有限元程序設(shè)計(jì)（三角形單元） 三角形單元是平面問(wèn)題中最簡(jiǎn)單的單元形式，但其他單元形式的有限元程序與三角形 單元程序的結(jié)構(gòu)基本相同，并無(wú)本質(zhì)區(qū)別。故本節(jié)通過(guò)一個(gè)FORTRAN 源程序，詳細(xì)介紹三 角形單元的程序設(shè)計(jì)。 1.1.

11、程序總框圖 程序總框圖 輸入結(jié)構(gòu)控制參數(shù) 輸入其它數(shù)據(jù) 形成整體剛度陣 引入支承條件 解方程，輸出位移 求應(yīng)力，輸出應(yīng)力 形成節(jié)點(diǎn)荷載向量 開(kāi)始 結(jié)束 1 單元面積 求彈性矩陣 單元?jiǎng)偠染仃?位移-應(yīng)變矩陣 2 3 4 5 6 7 8 9 10 上述總框圖左邊主線對(duì)應(yīng)著源程序中的主程序流程（見(jiàn)源程序）。 現(xiàn)按照總框圖中的10個(gè)子框圖，結(jié)合詳細(xì)的子框圖設(shè)計(jì)，介紹附錄中給出的“平面問(wèn)題三 角形單元源程序”。 2子框圖1（SDATA2子框圖1（SDATA） 輸入其他數(shù)據(jù)并計(jì)算半帶寬 輸入5個(gè)控制參數(shù)后，程序運(yùn)行中所需的其他原始數(shù)據(jù)均放在該子程序中輸入。其框 圖如下： 3 有限元講義 4 有限元講義

12、半帶寬的計(jì)算： 如圖所示結(jié)構(gòu)，9個(gè)單元10個(gè)結(jié)點(diǎn)。將單剛按結(jié)點(diǎn)分塊，采用直接剛度法以子塊對(duì)號(hào)入座 的方式可形成圖示總剛。 圖中整體剛度矩陣K的非零元素分布在以主對(duì)角線為中心的斜帶形區(qū)域內(nèi)(圖中用 粗線標(biāo)明),這種矩陣稱作帶形矩陣。在半個(gè)斜帶形區(qū)域中(包括主對(duì)角線元素在內(nèi)),每行 具有的元素個(gè)數(shù)稱為半帶寬,用d表示。由圖中看出,在半帶中,每行有五個(gè)子塊,即十個(gè)元 素,因此半帶寬d=10。半帶寬d的一般公式是： 半帶寬d=(相鄰結(jié)點(diǎn)碼的最大差值+1)結(jié)點(diǎn)位移未知量數(shù) 半帶寬d=(相鄰結(jié)點(diǎn)碼的最大差值+1)結(jié)點(diǎn)位移未知量數(shù) 三角形單元： 三角形單元： 半帶寬d=(相鄰結(jié)點(diǎn)碼的最大差值+1)2 半帶寬

13、d=(相鄰結(jié)點(diǎn)碼的最大差值+1)2 圖中相鄰結(jié)點(diǎn)碼的最大差值是4,故d=(4+1)2=10 根據(jù)帶形矩陣的特點(diǎn),并利用矩陣的對(duì)稱性,則在計(jì)算機(jī)中可只存貯上半帶的元素。這種存 貯方式稱為半帶存貯。 如下面左圖總剛，半帶存貯時(shí),只從K中取出上半部斜帶中的元素,存貯在右圖的矩 陣K *的豎帶中,矩陣K *為n行d列,只有nd個(gè)元素。因此,實(shí)際存貯量與K中元素總數(shù) 之比為d/n。由此可見(jiàn),d值愈小,則存貯量愈省。 由K改成K *時(shí)，元素的行碼不變，新的列碼改為： 新列碼=原列碼-行碼+1 （存上三角） 新列碼= （存下三角）？ 5 有限元講義 3子框圖2（STE）3子框圖2（STE） 計(jì)算單元?jiǎng)偠染仃?/p>

14、KE 6 有限元講義 4子框圖3（AE） 計(jì)算單元面積AE 4子框圖3（AE） 計(jì)算單元面積AE 5子框圖4（DTE） 計(jì)算彈性參數(shù)矩陣D 5子框圖4（DTE） 計(jì)算彈性參數(shù)矩陣D 7 有限元講義 6子框圖5（BTE） 計(jì)算幾何矩陣B 6子框圖5（BTE） 計(jì)算幾何矩陣B 8 有限元講義 7子框圖6（STIFF） 形成總剛度矩陣KS 7子框圖6（STIFF） 形成總剛度矩陣KS (1) 結(jié)點(diǎn)局部碼與總碼的對(duì)應(yīng)關(guān)系（結(jié)點(diǎn)換碼） 圖示為單元3個(gè)結(jié)點(diǎn)的兩種編碼，某個(gè)結(jié)點(diǎn)的局部碼如果是i(1),則它在總剛中的對(duì)應(yīng)總碼 是LND(IE,1)。 9 有限元講義 （2）同一子塊在單剛KE和總剛KS的對(duì)應(yīng)位置

15、（子塊搬家） （3）同一元素在KE、總剛KS和帶狀總剛KS *中的對(duì)應(yīng)位置（元素搬家） 現(xiàn)在討論形成總剛中某個(gè)元素的搬家問(wèn)題 設(shè)該元素屬于單剛矩陣KE中的I行J列子塊,該元素在這個(gè)子塊中的位置是II行JJ列。則 此元素在單元?jiǎng)偠染仃嘖E的位置應(yīng)為： 單元行碼： IH=2*(I-1)+II 單元列碼: L =2*(J-1)+JJ 此元素在整體剛度矩陣KS中的位置為： 整體行碼： ID=2*(LND(IE,I)-1)+II 整體列碼: IL=2*(LND(IE,J)-1)+JJ 此元素在半帶存貯的整體剛度矩陣KS *中的 位置為 半帶行碼： IDH=ID 半帶列碼: IDL=IL-IDH+1 8子

16、框圖7（EQUPE） 形成結(jié)點(diǎn)荷載向量P 8子框圖7（EQUPE） 形成結(jié)點(diǎn)荷載向量P 第一步：將作用在結(jié)點(diǎn)上的荷載送入P中 10 有限元講義 第二步：將作用在單元上的荷載轉(zhuǎn)化為等效結(jié)點(diǎn)力后送入P 單元自重引起的等效結(jié)點(diǎn)力(豎向力)： PE=-V(材料容重)*AE(單元面積)*T(厚度)/3 11 有限元講義 9子框圖8（INSCD） 引入支承條件修改總剛 9子框圖8（INSCD） 引入支承條件修改總剛 子框圖8對(duì)應(yīng)于子程序INSCD(子程序8)。其功能是引入支承條件, 引入支承條件就是 對(duì)帶形矩陣KS *和向量P進(jìn)行修改。 設(shè)第I個(gè)支承結(jié)點(diǎn)的結(jié)點(diǎn)碼為IR,即IR =JR(I,1),又設(shè)支承結(jié)

17、點(diǎn)IR的約束位移所對(duì)應(yīng)的 位移分量碼為II?，F(xiàn)分別說(shuō)明對(duì)KS、KS *和P的修改內(nèi)容。 1.對(duì)矩陣KS來(lái)說(shuō),第II行的主對(duì)角線元素應(yīng)改為1,該行的其它元素改為0.此外,II列的 全部非主對(duì)角線元素也改為0,如圖所示。 2.對(duì)半帶矩陣KS *(程序中仍用KS表示)來(lái)說(shuō),第II行中第一個(gè)元素應(yīng)改為1,該行其它元 素都改為0。此外,以該行第1個(gè)元素為起點(diǎn),向右上方畫(huà)45度斜線,則此斜線上的元素也都 改為0,如圖所示。斜線上的元素如果列碼是JJ,則其行碼為II一JJ+1。 為了確定斜線元素的最大列碼JO,需分兩種情況考慮： 在圖a中,行碼II半帶寬NW,故JO =NW, 在圖b中,行碼II半帶寬NW,

18、 故JO =II。 總括起來(lái),最大列碼JO應(yīng)等于NW和II兩個(gè)數(shù)中的較小值 3.將P中對(duì)應(yīng)行的元素改為零。 其框圖如下： 12 有限元講義 10子框圖9 10子框圖9 消元法解方程(略) 11子框圖10（SIGME）求應(yīng)力 11子框圖10（SIGME）求應(yīng)力 子框圖10對(duì)應(yīng)于子程序SIGME（子程序10）。其功能是計(jì)算單元應(yīng)力。 上面解方程時(shí)，得出結(jié)構(gòu)的節(jié)點(diǎn)位移，并存于P中。首先從P中取出單元e的節(jié)點(diǎn)位移 a e（程序中用DE表示），然后根據(jù)單元的節(jié)點(diǎn)位移ae求單元應(yīng)力。求應(yīng)力的基本公式 是 eee aBDaS= 求矩陣D，B時(shí)，要分別調(diào)研子框圖4，5 關(guān)于主應(yīng)力和主平面角的計(jì)算公式，可以參考

19、下圖的應(yīng)力圓得出： 平均應(yīng)力： 2/ )( yx ASG+= 應(yīng)力圓半徑： 2 2 4 )( xy yx RSG + = 最大主應(yīng)力：SGMA=ASG+RSG 最小主應(yīng)力： SGMI=ASG-RSG 13 有限元講義 主平面角： = = SGMI arcctg SGMI arcctg y xy y xy 29578.5790 180 0 0 下面寫出詳細(xì)框圖。 14 有限元講義 算例： 1.輸入數(shù)據(jù) 1）主程序中輸入結(jié)構(gòu)控制參數(shù) 結(jié)點(diǎn)總數(shù) NJ=10； 單元總數(shù) NE=9 支承節(jié)點(diǎn)數(shù) NR=2; 結(jié)點(diǎn)荷載數(shù) NPJ=1 15 有限元講義 問(wèn)題類型碼 IPS=0 （2）子程序1中輸入數(shù)據(jù) 彈性模

20、量 E=1.0; 泊松比 PR=0.25 單元厚度 T=1.0； 容重 V=0.0； = 1098 896 976 865 562 673 743 632 521 39 LND 支承節(jié)點(diǎn)數(shù)組JR 節(jié)點(diǎn)荷載數(shù)組PJ 16 = 114 111 32 JR 0 . 100 0 . 10 31 = PJ 2. 計(jì)算半帶寬NW 相鄰節(jié)點(diǎn)碼最大差值為4，NW=2*(4+1)=10 3.程序標(biāo)識(shí)符說(shuō)明 (1)整型變量 (1)整型變量 NJ-結(jié)點(diǎn)總數(shù)； NE-單元總數(shù)； NS-支承結(jié)點(diǎn)數(shù) NPJ-有荷載作用的結(jié)點(diǎn)數(shù), IPS-平面問(wèn)題類型碼(平面應(yīng)力問(wèn)題IPS=0，平面應(yīng)變問(wèn)題IPS=1) NJ2-位移分量總

21、數(shù),NJ2=NJ*2, NW-半帶寬；IE-單元序號(hào)。 (2)實(shí)型變量 (2)實(shí)型變量 E-彈性模量； PR泊松比； T-單元厚度； V-材料容重 有限元講義 AE-單元面積 SGX,SGY,TXY-應(yīng)力分量 ASG,RSG-平均應(yīng)力,應(yīng)力圓半徑, SGMA,SGMI-最大與最小主應(yīng)力, CETA-主平面角。 （3）數(shù)組 （3）數(shù)組 X(NJ)，Y(NJ) -結(jié)點(diǎn)坐標(biāo)數(shù)組, LND(NE,3)單元的節(jié)點(diǎn)碼數(shù)組 PJ(NPJ,3) - 結(jié)點(diǎn)荷載數(shù)組, JR(NS,3)支承節(jié)點(diǎn)信息數(shù)組 D(3，3)- 彈性矩陣D B(3，6) - 位移-應(yīng)變轉(zhuǎn)換矩陣（幾何矩陣B） S(3.6) - 位移-應(yīng)力轉(zhuǎn)換

22、矩陣（應(yīng)力矩陣S） ST(3) - 單元應(yīng)力數(shù)組; DE(6) - 單元結(jié)點(diǎn)位移向量, KE(6,6) - 單元?jiǎng)偠染仃嘖 e KS(NJ2,NW) - 整體剛度矩陣半帶存寬 P(NJ2)- 荷載向量、解方程后存放結(jié)點(diǎn)位移。 2.4 矩形單元 2.4 矩形單元 如圖示深梁,設(shè)劃分成具有八個(gè)矩形單元的網(wǎng)格，從中取出一典形單元e，單元有個(gè) 結(jié)點(diǎn)(i,j,m,p)。 一、局部坐標(biāo) 一、局部坐標(biāo) 設(shè)單元的邊界平行于結(jié)構(gòu)整體坐標(biāo)x, y， 邊長(zhǎng)分別為2a, 2b, 結(jié)點(diǎn)編號(hào)從左下角開(kāi)始逆時(shí) 鐘編成i，j, m, p。 單元局部坐標(biāo)的原點(diǎn)取在單元的形心(x0, y0)上,并采用無(wú)因次坐標(biāo),。,軸分別平 行

23、于x, y軸(如圖),局部坐標(biāo),與結(jié)構(gòu)坐標(biāo)x, y的關(guān)系是: () () 0 0 yy xx 1 1 b a = = ()142 x0 , y0 為坐標(biāo)原點(diǎn)。由此得到個(gè)結(jié)點(diǎn)的坐標(biāo)分別是（-1, -1)、（1,-1）,(1,1),(-1,1)。 17 有限元講義 二、結(jié)點(diǎn)位移列陣和結(jié)點(diǎn)力列陣 二、結(jié)點(diǎn)位移列陣和結(jié)點(diǎn)力列陣 每個(gè)結(jié)點(diǎn)2個(gè)位移分量，共個(gè)位移分量，設(shè)結(jié)點(diǎn)位移和結(jié)點(diǎn)力列陣分別為： T p ()242 pmmjjii e VUVUVUVUd= T p Y()342 pmmjjii e XYXYXYXF= 三、單元位移函數(shù)和形函數(shù) 三、單元位移函數(shù)和形函數(shù) 單元共有個(gè)位移分量， 將結(jié)點(diǎn)位移分

24、量全部作為已知邊界條件， 則位移函數(shù)可取為： yxv yxu 765 321 xy xy 8 4 += += ()442 它是單值連續(xù)的。在平行于x軸的直線上,位移分量是x的線性函數(shù)；在平行于y軸的直 線上它是y的線性函數(shù)，故 2-4-4稱為坐標(biāo)x, y的雙線性函數(shù)。 將各結(jié)點(diǎn)坐標(biāo)式代入上式 （如u(-a， -b)=ui ） ， 通過(guò)解方程組便可求得待定常數(shù) 1 ， 2 8 ，將這些參數(shù)代回式（2-4-4），經(jīng)整理得： ()()()() p U pmmjjii NUNUNUNu,+= ()()()() p V pmmjjii NVNVNVNv,+= (2-4-5) 以矩陣形式表示： (2-4-

25、6) e pmji pmji d NNNN NNNN v u f 0000 0000 = = e dN= 形函數(shù)矩陣N中的元素： ()=1 4 1 i N()1 ()+=1 4 1 j N()1 ()+=1 4 1 m N()+1 ()=1 4 1 p N()+1 18 如果引進(jìn)參數(shù): i = 0 , i = 0i (i=i, j, m, p), ( , j )是矩形單元4個(gè)結(jié)點(diǎn)的 局部坐標(biāo)。 結(jié)點(diǎn)i( i , i )的坐標(biāo)值分別是 i(-1,-1), j(1,-1),m(1,1), p(-1,-1)。 )， 代 入上式，則可將上式簡(jiǎn)記成： ()()( iii N+=11 4 1 ,) ()7

26、32()pm, ji, 或() 0 11 4 1 +=() 0 + 矩形單元的形函數(shù)，具有與三角形單元類似的特點(diǎn)，即： 在結(jié)點(diǎn)上的形函數(shù)值： 有限元講義 19 i Ni( i )=1 Ni j ( j )=Ni( m m )=Ni p ( p )=0 即, Ni + Nj + Nm + Np =1 在單元形心上, Ni(0, 0)=1/4 (i,j,m,p) 四、幾何矩陣B,應(yīng)力矩陣S 四、幾何矩陣B,應(yīng)力矩陣S 在三角形單元中，我們已經(jīng)推知： 應(yīng)變： e dB= 應(yīng)力: e dS e dBD= 矩形單元自然具有相同關(guān)系，只是B, S的內(nèi)容有所區(qū)別。 仿照三角形單元： pm B 00 00 0

27、 0 ji pi pi BBB NN NN xy y x NHB= = L L (2-4-8) 式中: i i N 0 = i N ba a b ab B 0 0 0 1 () () ()() ji ba a b ab iiii ii ii + + + = 11 10 01 4 1 ()pm 單元應(yīng)力矩陣： p S mji SSSBDS= 式中： () ()() ()() ()() i ba ab ab ab E S iiii iiii iii i , 1 2 1 1 2 1 11 11 14 2 + + + + = ()pmj, dVB 五、單元?jiǎng)偠染仃?五、單元?jiǎng)偠染仃?將上述公式代入單剛

28、的一般表達(dá)式 DBK T v = 對(duì)于矩形單元, 有限元講義 20 dd BDBabtdxdyBDBtK T T a a b b = 1 1 1 1 計(jì)算可得按結(jié)點(diǎn)分塊后的子矩陣為: () sr 3 1 + + + + + + + + = srsrsr srsr srsrsrsrsrsr rs a b b a b a a b Et K 1 2 1 3 1 1 2 1 2 1 3 1 1 2 1 3 1 1 14 2 ()pmjisr, = 2.5 六結(jié)點(diǎn)三角形單元 2.5 六結(jié)點(diǎn)三角形單元 一、面積坐標(biāo) 一、面積坐標(biāo) 在開(kāi)始單元推導(dǎo)之前，先介紹有關(guān)面積坐標(biāo)的概念。 1面積坐標(biāo)的定義 在第二節(jié)中

29、，我們討論過(guò)三角形單元形函數(shù)的幾何意義，如果從三角形單元e內(nèi)任意 點(diǎn)p向三個(gè)頂點(diǎn)i,j,m 引連線,使其分割成三個(gè)小三角形便可得到這些小三角形面積與形函 數(shù)的關(guān)系。 A A N i i = A A N j j = A Am =Nm 反過(guò)來(lái)說(shuō)，只要給定了三角形面積 Ai , Aj , Am， p點(diǎn)即可確定。因?yàn)?Ai, Aj ,Am 中 只有兩個(gè)是獨(dú)立的(AiAjAm=A)，這與P點(diǎn)坐標(biāo)(x, y)對(duì)應(yīng)。因此，我們可以用三角形面 積來(lái)定義三角形內(nèi)的任意一點(diǎn)P的位置，并稱它為面積坐標(biāo)。P點(diǎn)的面積坐標(biāo)定義為以下三 個(gè)面積比: A A L i i = A Aj Lj= A A L m m = 由此定義

30、,必有：Li +Lj +Lm =1。由于面積坐標(biāo)Li , Lj , Lm 只限于用來(lái)確定三角形單 元上(內(nèi))點(diǎn)的位置（只有當(dāng)L(i,j,m)1才有定義）,并且只是在單元分析中使用,所以L(i,j,m) 屬于局部坐標(biāo)。 （直角坐標(biāo)與坐標(biāo)原點(diǎn)、平移轉(zhuǎn)動(dòng)等有關(guān)，可因人而異，故稱它為人為坐標(biāo)，而面積坐標(biāo) 具有不變性，故稱其為自然坐標(biāo)。） 2面積坐標(biāo)Li和整體坐標(biāo)x,y之間的關(guān)系 面積坐標(biāo)Li和整體坐標(biāo)x,y之間的關(guān)系為: () ()xxyyyxyx A yx yx yx AA A L mmjjmmj mm jj i i += 2 1 1 1 1 2 1 ()yxj ()mj, i 有限元講義 21 或

31、: ()ycxba A L iii += 2 1 i ()mji, ()152 將上式與(2-1-3) Ni =1/2A(ai+bix+ci y)比較得知,三角形 常應(yīng)變單元的形函數(shù)Ni , Nj , Nm 分別等于三角形面積坐 標(biāo)Ni , Nj , Nm ,因此在三角形的三個(gè)頂點(diǎn)上亦具有與形函 數(shù)相同的結(jié)點(diǎn)坐標(biāo): () = =yxL ssr 0 1 , sr sr 三角形內(nèi)部點(diǎn)和邊界點(diǎn)上的面積坐標(biāo)值可以用右圖 中的一些等值線表示（面積坐標(biāo)的等值線）。從圖中可看 出：在平行于jm邊的任何一條直線上各點(diǎn)都具有相同的Li 值(可從幾何作圖上看出)其值等于該直線到對(duì)應(yīng)jm邊界 的距離與結(jié)點(diǎn)i到j(luò)m邊

32、界距離的比值 其值等于該直線到對(duì)應(yīng)jm邊界 的距離與結(jié)點(diǎn)i到j(luò)m邊界距離的比值。其它亦然。 3面積坐標(biāo)與直角坐標(biāo)的關(guān)系 將式(2-5-1)寫成矩陣形式得到用x、y表示的面積坐標(biāo): y x 1 = cba cba cba A L L L mmm jjj iii m j i 2 1 ()252 其反變換式為: m j i L L L = mji mji yyy xxx y x 1111 ()352 上式的數(shù)學(xué)意義是:三角形ijm內(nèi)某點(diǎn)的坐標(biāo)x,y,可以用三角形頂點(diǎn)(3個(gè)結(jié)點(diǎn))的坐標(biāo) (xi , yi ), (xj , yj ), (xm , ym )進(jìn)行函數(shù)插值而獲得,其插值多項(xiàng)式是三角形的面積坐

33、 標(biāo)Li , Lj , Lm， 也就是三角形的形函數(shù)Ni , Nj , Nm 。這種坐標(biāo)（x, y）表達(dá)式的插值 函數(shù)和位移場(chǎng) （u, v） 表達(dá)式的插值函數(shù)采用同一函數(shù)的單元稱為等參單元 坐標(biāo)（x, y）表達(dá)式的插值 函數(shù)和位移場(chǎng) （u, v） 表達(dá)式的插值函數(shù)采用同一函數(shù)的單元稱為等參單元(簡(jiǎn)稱等參元)。 （又稱：坐標(biāo)變換函數(shù)和位移插值函數(shù)中的形函數(shù)相同）。 “等參元”的另一種定義是：（參考書(shū)：湖南科技社有限元法概論 ） 設(shè) 位移變量u的近似式: = = s i u 1 iiu ( ) i x u是用s-1次插值函數(shù)近似表示的。坐標(biāo)變換式： r i i x = = 1 通常,坐標(biāo)變換式的次

34、數(shù)與u 近似式的次數(shù)并不相等，由此可將單元分成三類。 有限元講義 、子參元 rs 又名 亞參元（Sub-parametric element） 、等參元 r =s (isoperimetric element) 、超參元 r (Super-parametric element) 亞參元概念： 有時(shí)為了減少結(jié)點(diǎn)坐標(biāo)數(shù)據(jù)的輸入量，用于位移描述的結(jié)點(diǎn)數(shù)和用于坐標(biāo)變換的結(jié)點(diǎn) 數(shù)可以不一樣。如圖示四邊形單元，可將位移取為8個(gè)結(jié)點(diǎn),即位移函數(shù)為 22 iiu N iiv N iix f iiy f T v3 T y3 2 6y 2 12y i u = = 8 1 i v = = 8 1 而坐標(biāo)變換式則取為

35、: ijmp x = ijmp y = 我們稱這種坐標(biāo)結(jié)點(diǎn)少,而位移結(jié)點(diǎn)多的單元為亞參元 二、六結(jié)點(diǎn)三角形單元 二、六結(jié)點(diǎn)三角形單元 在三角形單元三條邊界的中點(diǎn)上各增設(shè)一個(gè)結(jié)點(diǎn)，使三角 形單元共有12個(gè)位移自由度。 1、結(jié)點(diǎn)位移列陣： (2-5-4） mmii e vuvuvud 11 LL= 結(jié)點(diǎn)力陣: （2-5-5） mmii e yxyxyxF 11 LL= 、單元位移函數(shù)和形函數(shù) 設(shè)為次多項(xiàng)式： 2 54321 xxyyxu+= 2 1110987 xxyyxv+= 為了計(jì)算簡(jiǎn)單,采用面積坐標(biāo)Li, Lj ,Lm(Li=1/2A(ai+bix+ciy) (i,j,m) ,單元上六個(gè)結(jié)點(diǎn)

36、的 面積坐標(biāo)計(jì)算如圖示。 將六個(gè)結(jié)點(diǎn)的12個(gè)位移分量代入， 并采用面積坐標(biāo)， 按前面類似的推導(dǎo)過(guò)程，便可得到用形函數(shù)表達(dá)的單元 有限元講義 位移函數(shù)： 23 N 33v NvNvNv ijii33u uNuNu jjii +=L =+L e dN （2-5-6） 寫成矩陣形式: i ji ji v u NNN NNN v u f= = = 3 3 3 000 000 L L 式中形函數(shù): ()12= ii LN i Lmji, mjL L4()imj, () N1= , 3 , 2 , 1 、幾何矩陣B, 應(yīng)力矩陣S 應(yīng)變=Bd x u x = y v y = x v y u xy + = 面

37、積坐標(biāo)函數(shù)的運(yùn)算: 設(shè)函數(shù) Z= f(Li Lj Lm )， 由于面積坐標(biāo)Li, Lj Lm 都是x,y的函數(shù)。 因此將面積坐標(biāo)對(duì)直角坐標(biāo)求導(dǎo)時(shí)，應(yīng)按復(fù)合函數(shù)求導(dǎo)法則： x Lm m L z x L L z x L L z x z j j i i + + = 同理: y Lm m L z y L L z y L L z y z j j i i + + = 由(2-5-1) (Li=(ai+bix+ciy)/2A) (i j m) A b x Li 2 = ()mji i 所以，面積坐標(biāo)函數(shù)對(duì)直角坐標(biāo)的偏導(dǎo)數(shù)為: m m L z + + = j j i i b L z b L z b Ax z

38、2 1 ()752 + + = m m j j i i L z c L z c L z c Ay z 2 1 同理: () () i L A c y N L A b L z A b x N i ii i i i ii 14 2 14 22 = = = ()mj 有限元講義 () () ()imj A LcLc y N A LbLb x N jmmj jmmj ,321 2 4 2 4 1 1 + = + = 由此得： ()()() m u1 mmjjjiiix LbuLbuLb Ax u 41414 2 1 += =+ +()()() 3 uLi 21 444bLbuLbLbuLbLb jji

39、miimjmmj + 24 同理可得y和 xy 的類似表達(dá)式，于是得到六結(jié)點(diǎn)三角形單元的幾何方程: d()852BdBBBBBB mji = 321 = 式中： () () () () = 414 40 014 2 1 iiii ii ii i LbLc Lc Lb A B 1 1 ()mji () () () j j L L () + + + = mmjjmmj mmj jmmj bLbLcLc cLc LbLb A B 44 40 04 2 1 1 ()i m j 3 2 1 dSdBDD xy y x = = 單元應(yīng)力方程: ()952 應(yīng)力矩陣： 3 S 21 SSSSSS mji =

40、 式中: () () ()() i i b = i ii i i i c cb cb A LE S 11 22 22 14 14 2 ()mji () ()() ) () jm j j L L ()imj321 ()( ()()() + + + = mjjmmj mmjjmmj mmjjmmj cLbLcLc LcLcLbLb cLcLbLb A E S 1414 88 88 14 2 1 從（2-5-8、9）兩式可見(jiàn)，六結(jié)點(diǎn)三角形單元的應(yīng)力和應(yīng)變分量均為整體坐標(biāo) x,y或 面積坐標(biāo)L的一次函數(shù),因此它比常應(yīng)變單元精度好。而且，它們沿任何方向又都是線性變 化的，能滿足邊界的變形協(xié)調(diào)條件，即單元的兩條收斂準(zhǔn)則,完備性和協(xié)調(diào)性都能得到滿 足。 、單元?jiǎng)偠染仃?有限元講義 單剛為1212階，由一般公式K