介老師第一課有限元計(jì)算探討

1、PAGE PAGE 6有限元編程和計(jì)算中的若干問題探討介玉新（清華大學(xué)水利水電工程系，北京，100084）【摘要】針對(duì)有限元編程和計(jì)算中的若干問題進(jìn)行探討，指出了一些應(yīng)當(dāng)注意的概念和細(xì)節(jié)，比如坐標(biāo)軸和符號(hào)約定、流量邊界和流速邊界名稱的辨析、固結(jié)分析中流量邊界的處理等。也指出了數(shù)值積分中的一些錯(cuò)誤之處。結(jié)合編程和計(jì)算經(jīng)驗(yàn)，給出了程序不能正常運(yùn)算時(shí)的一些解決方法，以及人們對(duì)有限元認(rèn)識(shí)的誤區(qū)與原因分析。【關(guān)鍵詞】有限元；固結(jié)分析；數(shù)值積分；程序1 前言筆者近年來主要從事土工數(shù)值計(jì)算等工作，進(jìn)行有限元編程和數(shù)值分析。在科研和教學(xué)中有一些簡(jiǎn)單的經(jīng)歷和體會(huì)。在此談一下，希望能夠拋磚引玉，與從事相同工作的同

2、仁共勉。文中談到的問題相信大多數(shù)人都曾經(jīng)發(fā)現(xiàn)過，只是沒有說出來而已。2 坐標(biāo)軸及符號(hào)約定xyz圖1 坐標(biāo)軸一般彈塑性力學(xué)分析中約定坐標(biāo)軸為右手螺旋方式，如圖1。三個(gè)坐標(biāo)軸方向滿足，。在編程中應(yīng)當(dāng)遵循這種隱含約定，否則可能會(huì)使面積計(jì)算出現(xiàn)負(fù)值以及其它一系列的問題。xyz圖1 坐標(biāo)軸在堤壩二維滲流分析中，對(duì)xy坐標(biāo)系常常隱含約定y坐標(biāo)為豎直方向。在一般地下水二維滲流分析的文獻(xiàn)中，xy坐標(biāo)則常常定義在水平面內(nèi)。注意到這些細(xì)微差別，對(duì)閱讀相關(guān)文獻(xiàn)會(huì)有好處。土力學(xué)中應(yīng)力應(yīng)變的符號(hào)規(guī)定與一般彈塑性力學(xué)符號(hào)規(guī)定相反。在編程時(shí)，建議采用一般彈塑性力學(xué)的符號(hào)規(guī)定。雖然這種規(guī)定使固結(jié)分析時(shí)土力學(xué)中習(xí)慣為正的孔隙水

3、壓力變?yōu)樨?fù)號(hào)，但能夠使之在整個(gè)有限元計(jì)算中與一般彈塑性力學(xué)中的定義保持一致。如果將孔隙水壓力按土力學(xué)符號(hào)定義，則容易造成公式推導(dǎo)和編程中的很大困擾。原因是在連續(xù)介質(zhì)力學(xué)中坐標(biāo)軸、邊界荷載、應(yīng)力分量以及位移、應(yīng)變等等是在統(tǒng)一的正負(fù)號(hào)約定下定義的，牽一發(fā)而動(dòng)全身。要盡量保持符號(hào)約定的一致性。3 流量邊界還是流速邊界在滲流計(jì)算中，第二類邊界條件為流量邊界，常寫為 ， 或 （1）式中q（或）為單位面積上的滲流量。、為邊界外法線的方向余弦。一些人覺得既然是“單位面積”上的流量，為何不稱作“流速”邊界。原因在于“流速”是矢量，“流量”是標(biāo)量。邊界條件中“單位面積上的滲流量”其實(shí)是流速在邊界外法線方向上的分

4、量，等于，即與外法線方向單位矢量的點(diǎn)積。所以將此邊界條件稱為“流量”邊界比叫做“流速”邊界更貼切些。 很多貌似合理的直觀感覺其實(shí)是經(jīng)不起推敲的。比如人們總以為制造衣服或地毯的纖維絲比土軟，其實(shí)其楊氏模量是土的幾百至上千倍。又如為了說明自然界存在很大的吸力，人們常拿上百米高的大樹舉例，但如果樹木等將水份輸運(yùn)到頂部只是依靠力學(xué)上的吸力，那么樹上的枯枝為什么沒有得到水分？還有一旦植物死去，為什么哪怕一棵小草也會(huì)干枯呢？4 數(shù)值積分中的錯(cuò)誤筆者在編程時(shí)，發(fā)現(xiàn)二維三角形的Hammer積分當(dāng)積分點(diǎn)數(shù)為4時(shí)，權(quán)系數(shù)應(yīng)分別為-27/96和25/96，而不應(yīng)當(dāng)是書中所寫的-27/48和25/48。兩者相差2倍。

5、筆者曾就此詢問文獻(xiàn)1的第一作者，作者的回答是在寫書時(shí)其實(shí)注意到了這個(gè)問題，但是看到國內(nèi)外很多書中都是如此，所以就照直采用。這個(gè)事情讓人費(fèi)解。筆者核對(duì)了許多書，發(fā)現(xiàn)幾乎所有的書本都在這里犯了錯(cuò)誤。尊重原作者固然是一種美德，但容忍其錯(cuò)誤卻并不值得提倡。其實(shí)這個(gè)錯(cuò)誤很容易就可以發(fā)現(xiàn)。不需推導(dǎo)公式，只需將被積函數(shù)取為1，就能輕易地發(fā)現(xiàn)問題。相信第一個(gè)得到這個(gè)公式的人應(yīng)當(dāng)是正確的，但在輾轉(zhuǎn)反抄的過程中何時(shí)引用錯(cuò)了就不得而知。相信這許多年來也有不少人發(fā)現(xiàn)了這個(gè)錯(cuò)誤，但不知為何大家都沒有在書中將它改過來。文獻(xiàn)2中作者用“*”做注，指出了正確的權(quán)系數(shù)計(jì)算值，但仍沒有直接在原表格中進(jìn)行更正。文獻(xiàn)中的錯(cuò)誤在所難免

6、。ASCE的Journal of Geotechnical and Geoenvironmental Engineering 中也有不少打印錯(cuò)誤。在程序編寫時(shí)，將用到的公式親自推導(dǎo)一遍，就能避免很多麻煩。另外，一些簡(jiǎn)單方法也有助于方便地發(fā)現(xiàn)錯(cuò)誤，比如，a 量綱分析：檢查量綱是否匹配；b 檢查對(duì)稱性：公式是否對(duì)稱，x、y、z方向三個(gè)分量輪換是否改變公式含義；c 用特殊數(shù)值進(jìn)行檢驗(yàn)：比如取變量為0、1、或等進(jìn)行分析，就可以檢驗(yàn)公式是否有誤。 另外，多找?guī)妆緯容^著學(xué)習(xí)，也有利于發(fā)現(xiàn)問題，避免錯(cuò)誤。5 固結(jié)計(jì)算中的問題比奧固結(jié)理論的平衡方程為： （2）其中p為孔隙水壓力。上述及以下的符號(hào)規(guī)定均按一般

7、彈塑性力學(xué)的符號(hào)規(guī)定。連續(xù)性方程， （3）其中vx、vy、vz為流速，v為體應(yīng)變。根據(jù)達(dá)西定律，有， ， 其中h為總水頭，（式中負(fù)號(hào)是由于孔隙水壓力p按一般彈塑性力學(xué)定義所引起）。由此代入式（3）可以得到以孔壓表示的連續(xù)性方程： （4）以式（2）、（4）為基礎(chǔ)可以推導(dǎo)出相應(yīng)的有限元格式： （5a） （5b）注意連續(xù)方程的流量邊界條件，用孔壓表示變?yōu)?（6）有些書中漏掉了上式中的項(xiàng)，將有限元格式（5b）的右邊項(xiàng)寫為 （7）其中，為形函數(shù)。一些程序中也犯了上述錯(cuò)誤。這個(gè)錯(cuò)誤雖然幸運(yùn)地尚不影響對(duì)超靜孔壓?jiǎn)栴}的計(jì)算，但將導(dǎo)致在考慮總孔壓的計(jì)算中在不透水邊界處有流量不為0的錯(cuò)誤。事實(shí)上，不但在不透水邊界

8、上，而且在不同的內(nèi)部單元邊界上也會(huì)出錯(cuò)。相信很多人都曾發(fā)現(xiàn)了這個(gè)問題。但在修改上一些人可能走了彎路。考慮項(xiàng)后，式（7）應(yīng)當(dāng)改寫為 （8）于是，除了邊界流量外，需要針對(duì)每個(gè)單元邊界計(jì)算。除了在整個(gè)計(jì)算域的邊界上要對(duì)該式進(jìn)行計(jì)算外，在不同的內(nèi)部單元邊界上也要進(jìn)行這項(xiàng)操作（在相同的內(nèi)部單元邊界上倒不必處理）。這種計(jì)算是相當(dāng)繁瑣的，且極易在編程中出現(xiàn)失誤。其實(shí)只要利用高斯公式將邊界積分改寫為體積分，即將式（8）進(jìn)一步寫為 （9）其中 ， 就能簡(jiǎn)便地解決這個(gè)問題。因?yàn)轶w積分計(jì)算要比邊界積分計(jì)算方便許多。6 程序不能正常運(yùn)算的解決辦法筆者也常遇到有人詢問程序算不下去的情況。此時(shí)往往是輸入文件有誤。在確認(rèn)輸

9、入文件正確的情況下，應(yīng)當(dāng)檢查單元和節(jié)點(diǎn)數(shù)目是否超過相關(guān)數(shù)組所定義的數(shù)值。還有一種情況就是被0除。遇到分母為0的情況，程序往往中斷運(yùn)行。編程中消除被0除的情況可以采用如下方法：對(duì)于這種算式，如果分母b 0時(shí)，同時(shí)有a 0和y 0，此時(shí)可以在上述算式前加上一句“時(shí)，取的語句”（其中為一小量，式中正負(fù)號(hào)根據(jù)b本身的正負(fù)號(hào)確定）。這樣，將b取為一小量，可以避免被0除的情況；如果b 0時(shí)，y趨于一常數(shù)或，程序中可以按和兩種情況用“ifelseendif”等選擇語句分別加以處理；如果時(shí)處理比較繁瑣，且b=0的情況非常罕見，也可以不予處理，只在算式前加上諸如“如果，輸出警告信息”等語句，以便于一旦程序算不下

10、去時(shí)排查出錯(cuò)原因。對(duì)于等算式，在計(jì)算前加上“時(shí)，取”及“時(shí)，取”的語句。將x限制在其定義域內(nèi)，也能避免程序算不下去的困擾。同時(shí)將x超出其定義域的情況輸出出來，便于計(jì)算完成后分析其對(duì)整個(gè)計(jì)算的影響。對(duì)于等算式，也有類似對(duì)定義域要求的情況需要處理。對(duì)于應(yīng)力計(jì)算中有時(shí)遇到的，表面上是一個(gè)簡(jiǎn)單算式，只要分母不為0即可正常計(jì)算，但考慮到該式其實(shí)等于，其中為L(zhǎng)ode角，所以也應(yīng)當(dāng)限制其值域在-11之間，以避免導(dǎo)致其它意外錯(cuò)誤。總之，編程時(shí)只要注意處理好分母為0、定義域、值域、數(shù)組越界等情況，相信能夠解決絕大多數(shù)程序算不下去的困擾。養(yǎng)成良好的編程和應(yīng)用程序的習(xí)慣，按部就班，一般就不會(huì)出現(xiàn)疏漏。7 對(duì)有限元認(rèn)

11、識(shí)的誤區(qū)及分析高等數(shù)學(xué)中微積分的本質(zhì)是由連續(xù)體中一個(gè)點(diǎn)或單元體的信息去推求二維（面）或三維（體）的特征，其紐帶是偏微分方程。求解偏微分方程的方法可以是解析法，也可以是數(shù)值解法。在計(jì)算機(jī)出現(xiàn)之前，解析方法占有主導(dǎo)地位，被數(shù)學(xué)家和力學(xué)家發(fā)揮得淋漓盡致。計(jì)算機(jī)的發(fā)展帶動(dòng)了數(shù)值方法的發(fā)展。應(yīng)當(dāng)說有限元的發(fā)明是數(shù)值方法的最大成就。對(duì)絕大多數(shù)力學(xué)問題，有限元法是當(dāng)前惟一有效的方法。微積分是解決由點(diǎn)到面（體）推廣的最有效的數(shù)學(xué)解析工具，有限元?jiǎng)t是解決這一問題的最成功的數(shù)值解算手段。有限元的本質(zhì)是求解偏微分方程問題，即用數(shù)值方法解偏微分方程。對(duì)土工數(shù)值計(jì)算來說，除了好的程序之外，正確的計(jì)算還需要把握好以下各個(gè)

12、環(huán)節(jié)：對(duì)問題本身的理解是否清楚？對(duì)程序的特點(diǎn)是否真正掌握？試驗(yàn)和參數(shù)是否恰當(dāng)，可信度有多少？對(duì)問題的簡(jiǎn)化是否抓住了主要矛盾且符合實(shí)際？有限元網(wǎng)格的劃分是否合理？對(duì)邊界條件的處理是否恰當(dāng)？數(shù)據(jù)文件是否嚴(yán)格按程序要求準(zhǔn)備，并正確輸入？輸出成果是否違背自己的預(yù)期？任何一個(gè)環(huán)節(jié)出錯(cuò)，都有可能導(dǎo)致求解錯(cuò)誤。目前對(duì)有限元有一些誤區(qū)，一種是事事都希望用有限元算一下，另一種則是對(duì)有限元加以排斥，說它是糊弄湊數(shù)。但總的來講后一種情況可能更多一些。分析其原因，可能在于：（1）因有限元公式符號(hào)復(fù)雜，涉及內(nèi)容太多，令人望而生畏。由畏生怨，因此予以拒絕排斥。這在有限元發(fā)展之初可能較多，隨著有限元在本科和研究生教學(xué)中的普

13、及，人們大都受過這方面的良好教育，在接受有限元方面已不成問題；（2）不恰當(dāng)?shù)氖褂谩Ｓ捎谑褂谜邔?duì)有限元以及待解決的問題認(rèn)識(shí)不清，或參數(shù)問題，或?qū)τ?jì)算結(jié)果不能正確解釋，因此導(dǎo)致錯(cuò)誤。這可能是有限元招致批評(píng)的主要原因；（3）過分期待。理論上有限元可以解決巖土力學(xué)中的所有問題，但事實(shí)上并不能事事都能解決或解決得很好。土不是鋼材那樣的彈性材料，單單勘察取樣和土的本構(gòu)關(guān)系就會(huì)制約數(shù)值計(jì)算的準(zhǔn)確性。過分要求，讓有限元做力所不能及的事情，就容易使有限元因出錯(cuò)而受到批評(píng)；（4）直覺錯(cuò)誤。計(jì)算結(jié)果與預(yù)期或直覺不符，因此產(chǎn)生對(duì)有限元的不信任感。這種事情更多發(fā)生在生產(chǎn)單位。他們實(shí)踐經(jīng)驗(yàn)豐富，所以很多情況下他們的直覺或

14、經(jīng)驗(yàn)是對(duì)的。有限元計(jì)算不準(zhǔn)往往與參數(shù)選擇有關(guān)，不能怪有限元本身。有的時(shí)候，人的直覺以及經(jīng)驗(yàn)其實(shí)也是錯(cuò)誤的。另外，工程師們已習(xí)慣于對(duì)強(qiáng)度問題的計(jì)算，對(duì)與變形相關(guān)的有限元計(jì)算，總習(xí)慣上要與根據(jù)強(qiáng)度計(jì)算得到變形結(jié)果進(jìn)行對(duì)照。土工有限元計(jì)算中施工順序?qū)ψ冃蔚挠绊懞艽螅瑥?qiáng)度問題計(jì)算卻不必考慮施工順序。概念和理解以及對(duì)直觀現(xiàn)象解釋上的錯(cuò)誤常常使工程師們非常困惑，這也導(dǎo)致有限元受到一些誤解。8 結(jié)束語隨著計(jì)算機(jī)以及以有限元為代表的各種計(jì)算手段的發(fā)展，數(shù)值計(jì)算已深入人心。人們對(duì)有限元計(jì)算也不再認(rèn)為是單純的數(shù)字游戲。事實(shí)上，有限元分析從邊界條件簡(jiǎn)化、模型參數(shù)選擇、輸出結(jié)果判讀等各個(gè)環(huán)節(jié)，也同工程設(shè)計(jì)施工一樣，需要深厚的理論基礎(chǔ)和敏銳的經(jīng)驗(yàn)判斷，決不是簡(jiǎn)單的數(shù)據(jù)輸入和結(jié)果輸出。把數(shù)值計(jì)算理想化或簡(jiǎn)單化都是不對(duì)的。近年來，隨著商業(yè)軟件在國內(nèi)的傳播，使用商業(yè)軟件漸漸成為一種趨勢(shì)。但目前的商業(yè)