連續(xù)體的有限單元法課程理論與上機(jī)課程建設(shè)

張一鳴，王雪雅，馮 春，叢俊余

（1.河北工業(yè)大學(xué) 土木與交通學(xué)院，天津 300401；2.中國科學(xué)院力學(xué)研究所，北京 100190；

3.北京極道成然科技有限公司，北京 100085）

有限單元方法自上世紀(jì)50 年代誕生后，伴隨計算數(shù)學(xué)與計算機(jī)科學(xué)的蓬勃發(fā)展，已成為工業(yè)界的首選數(shù)值分析理論與方法。同時其數(shù)學(xué)完備，相對直觀易理解，計算框架擴(kuò)展性與可靠性強(qiáng)，是工科本科生與研究生入門數(shù)值計算的最佳選擇。

另一方面，針對土木工程學(xué)生，傳統(tǒng)的有限單元法往往更強(qiáng)調(diào)結(jié)構(gòu)單元的掌握與教學(xué)。例如梁、板、殼類型單元的構(gòu)造、特性與使用。然而，與建筑工程及橋梁工程不同，在巖土與水利工程中，計算分析中時常會采用連續(xù)體單元，即實體單元來模擬巖土體介質(zhì)，并結(jié)合精細(xì)化建模來研究工程問題。同時伴隨計算機(jī)算力的提升，應(yīng)用連續(xù)體單元、考慮非線性過程同時結(jié)合多尺度方法，來解決大規(guī)模復(fù)雜工程問題這一路徑也得到了越來越多科研人員的青睞?？梢韵胂笤诮酉聛砗荛L一段時期內(nèi)，連續(xù)體單元在工程與科研中會得到更為廣泛的應(yīng)用。

本文即針對如何開展、教授連續(xù)體的有限單元方法展開具體的規(guī)劃與討論。關(guān)于課程的結(jié)構(gòu)設(shè)計，我們參考了面向?qū)ο蟮木幊踢^程，從多個“關(guān)鍵知識點”出發(fā)，注重“關(guān)鍵點”之間的鏈接與邏輯關(guān)系，尤其是可能存在的繼承關(guān)系。同時與其他課程例如數(shù)值計算與有限元加以區(qū)別和類比。在理論授課時深入淺出，以線彈性問題為主，材料非線性問題為輔，幾何非線性問題為參考；透徹一維問題，主講二維問題，選講三維問題。另一方面，雖然本課程題目是關(guān)于“連續(xù)體”，但是也會討論有非連續(xù)，即裂縫的情況，也會討論一些特殊的計算方法與計算策略，如內(nèi)嵌不連續(xù)單元及破裂單元法構(gòu)造。上機(jī)實踐是本課程教學(xué)的一大重點，本文認(rèn)為應(yīng)盡可能采用較為基礎(chǔ)高效的編程語言與計算軟件，并以開源編譯器與開源軟件為主，雖然有一定的門檻，但是對于學(xué)生后期成長幫助極大，也更適用于團(tuán)隊化協(xié)作編程，在當(dāng)前國際形勢下也顯得尤為重要。本文后續(xù)章節(jié)將對課程的具體路徑、授課內(nèi)容與比例、上機(jī)操作與考試考察展開討論。

一、課程路徑

本文理論授課包含的關(guān)鍵點如圖1 所示，圖中實線部分表明關(guān)鍵點間存在較強(qiáng)的邏輯關(guān)系，授課過程中不宜忽略，而虛線指向的內(nèi)容應(yīng)了解概念，不作為主要授課內(nèi)容。部分授課關(guān)鍵點與其他課程可能有重疊，授課老師應(yīng)在第一節(jié)課時關(guān)注多課程的排課次序，避免授課內(nèi)容重復(fù)，而應(yīng)基于遞進(jìn)關(guān)系更新授課計劃，從而加深學(xué)生理解與認(rèn)識。

圖1 本文理論授課涉及的關(guān)鍵點

圖1 中部分較為基礎(chǔ)的內(nèi)容應(yīng)重點闡述，例如“實體單元”與“拉格朗日插值方法”以及“高斯積分方法”的內(nèi)在關(guān)聯(lián)。作者經(jīng)過多年教學(xué)認(rèn)為，有限單元法方法的基石可高度濃縮于兩大方法：弱形式條件?F=0→∫|?F-?G|dx=0 與形函數(shù)的引入F（x）=∑N（x）F。前者理論性較強(qiáng)，應(yīng)注重推導(dǎo)與表達(dá)，此處不再贅述。而后者實踐性較強(qiáng)，應(yīng)結(jié)合案例教學(xué)。在一維情況下的案例包括計算一個沿厚度方向的熱傳導(dǎo)問題，二維情況則應(yīng)分析平面彈性力學(xué)問題，側(cè)重于分析采用不同單元條件下計算結(jié)果的差異。實體單元雖然計算流程較復(fù)雜，但其整個計算流程又極為重要，通過紙面計算幾個案例可大幅度強(qiáng)化學(xué)生對于方法的理解。作者此處建議可采用較少的單元來計算較為基本的案例，例如年代較早的教科書提供了不少經(jīng)典案例與習(xí)題，作者將文獻(xiàn)[4]中的一道習(xí)題進(jìn)行擴(kuò)展，作為案例示于圖2 中。圖2 所示為一個一端固定，另一端受集中荷載或已知位移的深梁結(jié)構(gòu)，對于這個案例采用有限單元法進(jìn)行計算時，考慮3 種情況：（1）剖分為2 個三角形常應(yīng)變單元；（2）剖分為1 個四邊形雙線性單元；（3）剖分為一個四邊形二次單元。通過計算情況（1）讓學(xué)生初步了解連續(xù)體有限單元法計算流程，并且考慮集中荷載或已知位移來研究邊界條件的施加；通過計算情況（2）來引入高斯積分，分析高斯點的數(shù)量與對應(yīng)位置、權(quán)重，以及其對結(jié)果的影響；通過計算情況（3）來介紹高階單元，與程序計算結(jié)果進(jìn)行對比，闡明高階單元的適用性。另一方面，也可考慮采用具有對稱性的計算算例來進(jìn)一步簡化計算過程。

圖2 本文推薦的一個二維例題

該案例還可以進(jìn)一步擴(kuò)展到三維情況，采用單個三維六面體單元，考慮不同階次，如圖3 的單元示例。通過計算三維情況來分析平面應(yīng)變問題與平面應(yīng)力問題分別對應(yīng)的是哪種厚度方向的約束條件，加深理解。

圖3 算例中推薦采用的三維單元

當(dāng)涉及材料非線性問題時，可退回至一維算例或二維常應(yīng)變單元的算例?？紤]一個給定的應(yīng)力-應(yīng)變（σ-ε）關(guān)系如理想彈塑性，結(jié)合加載時間步，通過隱式（時間）積分，引入牛頓迭代方法，迭代前兩步來觀察收斂性質(zhì)。若學(xué)生基礎(chǔ)較好，也可介紹徑向返回算法（return mapping）等高級方法來研究分析材料非線性問題（建議一維）。通過顯式（時間）積分，進(jìn)一步介紹穩(wěn)定性條件以及動態(tài)松弛技術(shù)（由于動態(tài)松弛技術(shù)的適用性與特殊性，本文未將動態(tài)松弛方法歸類于數(shù)值計算工具包）。

最后，當(dāng)涉及非連續(xù)即裂縫的計算時，可闡述基于結(jié)點富集的擴(kuò)展有限元方法（XFEM）的處理方法或者基于單元富集的破裂單元法。其中擴(kuò)展有限元方法相對較為復(fù)雜，可能需要耗費較多學(xué)時，以應(yīng)用及上機(jī)操作為主，破裂單元法則相對較為簡單，由于有之前有限單元類型的鋪墊，授課時可將破裂單元法認(rèn)為是一種構(gòu)造類似9結(jié)點平面單元的特殊單元來展開，相關(guān)內(nèi)容可參考文獻(xiàn)[5]。

二、具體內(nèi)容比例

建議總課時數(shù)約為32 學(xué)時，其中理論課課時16 學(xué)時，上機(jī)課課時16 學(xué)時。若將課程作為本科生課程，建議作為三年級（下）課程，作為研究生課程時，建議作為碩士研究生課程。

理論課時分配：數(shù)值計算方法1 課時；有限單元方法理論1 課時；實體單元理論與算例計算4 課時；材料非線性理論與算例6 課時，其中隱式（時間）積分計算3課時，顯式（時間）積分計算3 課時；擴(kuò)展有限元2 課時；破裂單元法2 課時。

本文重點闡述上機(jī)課程分配：（1）數(shù)值計算工具箱4 學(xué)時，包括插值函數(shù)（一維與二維）2 學(xué)時，高斯積分計算（一維與二維）1 學(xué)時，牛頓迭代1 學(xué)時；（2）二維實體有限單元法（建議采用8 結(jié)點二維等參單元）上機(jī)指導(dǎo)編程4 學(xué)時，其中建模1 學(xué)時，單元設(shè)計1 學(xué)時，組裝與邊界條件設(shè)置1 學(xué)時，求解器調(diào)用與計算獲得結(jié)果1 學(xué)時；（3）非線性分析4 學(xué)時，包括隱式時間積分2 學(xué)時，顯式時間積分2 學(xué)時；（4）破裂單元法1 學(xué)時；（5）其他軟件如Abaqus，LS-Dyna，Ansys 等演示操作2 學(xué)時；（6）上機(jī)考試1 學(xué)時。

上機(jī)課程的具體內(nèi)容流程：插值函數(shù)僅教授拉格朗日插值，完成一維插值算例后，進(jìn)階二維三角形單元的面積坐標(biāo)插值，再完成四邊形雙線性插值函數(shù)構(gòu)造，同時引入等參單元，獲得N、?N 以及雅可比矩陣的表達(dá)。然后可以較為自然地進(jìn)入高斯積分環(huán)節(jié)，再轉(zhuǎn)而二維實體有限單元法，完成線彈性問題分析后進(jìn)入非線性分析，順序教授牛頓迭代與隱式積分，然后與顯式積分作對比，最后過渡到非連續(xù)問題與破裂單元法，上機(jī)部分的授課流程如圖4 所示。

圖4 上機(jī)操作內(nèi)容授課流程

三、上機(jī)編程要點

上一章節(jié)對于課程的上機(jī)操作考慮了學(xué)習(xí)的循序漸進(jìn)，數(shù)值計算工具箱部分可以讓學(xué)生了解各類子函數(shù)、全局變量、結(jié)構(gòu)體的書寫格式與調(diào)用方法。同時掌握“算法”的程序?qū)崿F(xiàn)方法，而具體的有限元程序書寫則能夠鍛煉學(xué)生對于“工程”的設(shè)計與規(guī)劃能力，即如何將具體功能的子函數(shù)組裝到模塊，實現(xiàn)特定功能，然后將模塊再拼裝，實現(xiàn)一個任務(wù)。

另一方面，如前所述，本文建議應(yīng)盡可能采用高效基本的編程工具，而非Matlab 一類的商業(yè)數(shù)值計算軟件。而事實上，一旦掌握了較為基礎(chǔ)的編程工具后，也可以大幅度降低學(xué)習(xí)Matlab 的難度。我們建議可選擇面向?qū)ο蟮腃++或者Python，面向過程的C 或者Fortran。事實上，隨著編程軟件的發(fā)展，現(xiàn)在兩種分類有一定模糊性，兩類編程語言也在互相借鑒。目前實體有限元程序有大量開源代碼供學(xué)習(xí)，如Calculix，Opensees，但是大部分起點較高，直接對源代碼進(jìn)行理解存在一定難度。無論采用面向?qū)ο筮€是面向過程編程語言，作者認(rèn)為實體有限元程序中存在諸多“要素”，我們將主要“要素”列于圖5 中。

圖5 實體有限元程序設(shè)計中的一些要點與功能

其中的幾何部分，與建模程序的幾何部分有一定相似性，但又有所不同。我們進(jìn)行了簡單的展開，每個要素可以根據(jù)不同情況進(jìn)行進(jìn)一步擴(kuò)展，其所包含的信息、函數(shù)、數(shù)據(jù)方法在不同條件下也存在差異，其他部分也可作類似展開。對于學(xué)生的考查可側(cè)重于對相關(guān)功能的要素歸類，或者子功能的具體實現(xiàn)。對于一些較為復(fù)雜的“黑箱式”功能，例如大規(guī)模線性方程組求解，可通過熟悉開源代碼的輸入和輸出格式，直接書寫API 接口取用，這樣可以鍛煉學(xué)生對于既有代碼的調(diào)用能力。

本課程的考試建議以上機(jī)操作為主，可設(shè)置多類必須通過編程來求解的選擇題與計算題，要求學(xué)生必須通過現(xiàn)場書寫代碼，調(diào)試，運行后獲得結(jié)果答題。例如提供結(jié)點位移求解某個單元某個高斯點的應(yīng)變值，提供材料信息計算某個單元的剛度矩陣，提供裂縫法方向求解破裂單元的整體剛度矩陣等。

四、結(jié)論與展望

CAE 即計算機(jī)輔助工程分析、設(shè)計、優(yōu)化，是一種新時代下的國之重器，各國對于CAE 核心技術(shù)都極為重視。而CAE 技術(shù)的自主研發(fā)離不開優(yōu)秀CAE 開發(fā)人員的發(fā)掘與培養(yǎng)，這是典型的工程-軟件設(shè)計交叉學(xué)科。本文對于連續(xù)體的有限單元法課程理論與上機(jī)課程進(jìn)行了設(shè)計與探討，對于授課內(nèi)容、邏輯結(jié)構(gòu)、課時分配、上機(jī)編程進(jìn)行了規(guī)劃。相信本課程的順利建設(shè)對于培養(yǎng)中國制造2025 背景下的優(yōu)秀工程人才有一定的推動作用。