有限元求解器的案例
考慮熱源的瞬態熱傳導有限元求解器
關鍵詞：熱源，瞬態，熱傳導，有限元求解器，三角形單元，自研
在《瞬態熱傳導有限元求解器開發》一文中，我們介紹了自研的二維瞬態熱傳導求解器。
當時那個控制方程沒有考慮熱源，邊界條件中只涉及溫度、熱流、對流。然而在很多問題中，熱源才是最關鍵的邊界條件，比如電發熱、化學反應生熱。
熱源的處理
熱源是體熱，相對應的熱流是面熱。兩者處理方式類似，都是根據單位熱功率值和幾何尺寸計算熱功率，然后加到控制方程矩陣的右側，承擔類似于結構力學中的“載荷”的功能。
區別在于，熱源是作用在體上的，單位是W/m3，熱流是作用在面上，單位是W/m2。具體到編程上，熱源要分配到單元的三個節點上，熱流要分配到單元某個邊的兩個節點上。
從求解器編程的角度來說，這些邊界條件的處理方式都是固定和通用的。考驗一般出現在實際工程項目中使用自研求解器的時候。
在CAE軟件的開發中，交互端和求解器端永遠要解決的問題是，如何讓所有單元始終知道：
（1）它是誰？（材料參數，幾何參數）；
（2）它在哪？（和其他單元的相對位置）；
（3）它怎么了？（邊界條件）。
以熱源為例，在交互界面上，我們通過視口選擇單元，指定其體熱功率。那么前端數據在生成求解器輸入的時候，就要告知求解器所有單元的編號和其對應的體熱功率。
當求解器拿到單元編號以后，就需要索引或者計算其面積，并根據單元三個節點編號，將功率加到載荷列陣對應的位置。
驗證
設計案例如下，區域外部為20℃空氣，對流換熱系數取5W/(m2K)，時間總長18000s，每步時間間隔60s。
自研求解器得到模型中心最終溫度是84.6℃，與商用軟件結果完全一致。
展開 流體有限元求解器開發-二維斯托克斯方程
我在開發結構力學有限元求解器的時候，都是先去查資料，直接就把單元剛度矩陣拿過來用。
但是到流體這就不能完全這么干，原因是：
（1） 未必能找到直接可用的單元矩陣；
（2） 流體的邊界條件中，有很多第二類邊界條件（壓力、熱流、對流），這些邊界條件有的是放到右側載荷項，需要推導。有的載荷甚至會影響左側單元矩陣，部分放到左邊去修正。結構力學里面大部分情況下，載荷就是在右側，位移就是放到左邊修正剛度矩陣，清晰明確。
所以，開發流體求解器的時候，還是要從有限元的基本方法入手。這里采用加權余量法進行處理。有限元的教材里面講的很多了，這里簡單說一下流程：
（1） 根據單元類型，確定插值函數。此時速度、壓力等變量，都可以用權函數表達。
（2） 采用伽遼金方法，權函數=插值函數，控制方程與權函數相乘，積分取0。
（3） 在每個單元域內，方程轉換為權函數的積分形式，最終形成單元矩陣。
單元方程
最終得到單元的方程形式如下：
類比到結構有限元，左側第一大項，就是剛度矩陣。u、v、p相當于3個自由度，右側就是載荷列陣。
我的本意是開發一個三角形單元的斯托克斯求解器。使用三角形線性單元對應的插值函數：
無論是CFD還是結構有限元，只要單元類型一致，插值函數都是一樣的，區別只在單元方程。
但是這樣直接求解結果是震蕩的，原因是結構有限元中，節點的幾個自由度本質是同一種物理量，它們是“平級的”。但是速度u、v與壓力p不是同一種物理量，它們不平級。壓力和速度的降階是平級的。
一般的處理思路是，對速度用高階單元，對壓力用低階單元。還有一種思路是使用罰函數方法，將壓力用如下形式表達，只要λ足夠大即可。這樣就在原控制方程中消除了壓力。
展開 瞬態熱傳導有限元求解器開發
有限元思路
這部分在結構有限元教材中介紹的比較多，流程：
（1） 根據單元類型，確定插值函數。此時單元溫度用權函數表達。
（2） 采用伽遼金方法，權函數=插值函數，控制方程與權函數相乘，積分取0。
（3） 在每個單元域內，方程轉換為權函數的積分形式，最終形成單元矩陣。
單元方程
使用三角形線性單元對應的插值函數：
有些教材中，會把面積項提取出來，寫成以下這種形式，所以有的教材上剛度矩陣結果用a、b、c表達的時候，會存在差異，但是本質都是一樣的。
最終單元方程如下，其中M是熱容矩陣，K是傳導矩陣，F是熱載荷。
熱容矩陣乘的是溫度的導數。在瞬態問題的求解中，導數項可以寫成前后時間變量差值與時間間隔的比值：
代入后得到如下形式：
求解思路
在求解過程中，把Tn+1當作未知量，Tn作為已知量。這樣在每個時間點，求解方法和結構有限元方法一致。
初始時候，可以指定一個溫度作為全域已知初始溫度，然后在迭代過程中，Tn和Tn+1會逐漸接近，達到收斂狀態。
案例效果
設計案例如下，同時包含對流換熱邊界條件和熱流，時間總長10000s，每步時間間隔50s。
自研求解器和商用軟件結果對比如下，從結果可以看出，自研求解器結果與商用軟件結果一致。
自研求解器結果：最終溫度分布
商用軟件結果：最終溫度分布
自研求解器結果：平均溫度時間曲線
商用軟件結果：平均溫度時間曲線
展開 流體有限元求解器開發-SUPG迎風格式與SA湍流模型
關鍵詞：CFD，有限元，對流項，繞流，迎風格式，湍流模型
在《流體有限元求解器開發-不可壓定常流動模型》一文中，我們介紹了考慮對流項的不可壓流動求解器的實現。
然而正如所預料的那樣，一旦流速高一些，或者粘性小一些，仿真結果就容易發散，收斂性成為一大難題。
為了解決這個問題，CFD大神們想出了各種手段，有的嚴格按照理論去處理盡力彌合。有的則主打靈感修正，問就是人工粘性、人工擴散、人工穩定，實用至上。
SUPG（Streamline Upwind/Petrov-Galerkin，流線迎風/Petrov-Galerkin）迎風格式是計算流體力學和有限元方法中一種經典的穩定化技術，專門用于解決對流主導問題中的數值振蕩問題。
該方法是79年到82年Brooks 和 Hughes提出并確立的，目前廣泛用于流體有限元求解中。
Hughes T J R
SUPG的核心思想
我們前面文章介紹的伽遼金法，在推導過程中，令權函數=插值函數。在對流主導情況下，這種對稱處理無法捕捉流動的方向性特征，因此迭代過程中，速度場逐漸發散。
SUPG的核心思想，是修改權函數，引入迎風效應。增加的項一個只在流線方向上起作用的項。我的理解是人工給一個收斂的方向。
當這個項用的系數大，抹平振蕩的能力就越強，當然結果也可能偏離實際更多。用的系數小，就可能會發散。
在這個基礎上，我們進一步嵌入了SA湍流模型，這是因為高雷諾數流動求解中發現，上述方法收斂性還是差。SA湍流的引入，可以將N-S方程的擴散項系數增大，對流主導問題的病態程度降低，迭代求解更容易收斂。和無腦增加迎風項系數強制收斂比，這種方法得到的結果精度要好一些。
效果
圓柱繞流
設定圓柱半徑為0.05m，流體介質為空氣，來流速度100m/s。
展開 
Ansys Mechanical | 軟件介紹：業界一流的有限元求解器
wx_fmt=gif&from=appmsg"></p><p class="ql-align-center"><br></p><h1><strong>Ansys Mechanical 有限元分析軟件 </strong></h1><p><br></p><p>Ansys Mechanical 是業界一流的有限元求解器，具有結構、熱學、聲學、瞬態和非線性功能，可幫助改進建模。</p><p><strong>NO.1軟件概覽</strong></p><p><img src="https://mmbiz.qpic.cn/mmbiz_jpg/XUJorxDYRKRVcgS91JRTGvWkdD1mpibZW6lB58iaZcnic2aHAs3oJ83hXlX5LTucWQRvYVxAVlwHRZUD93pkpVP6A/640?wx_fmt=jpeg&wxfrom=5&wx_lazy=1&wx_co=1"></p><p><br></p><p>Ansys Mechanical 創建了一個使用有限元仿真分析軟件（FEA）進行結構分析的集成平臺。Mechanical 是具有完整分析工具的動態環境，從準備幾何結構進行分析到連接其它物理場以實現更高的保真度。可定制的直觀用戶界面使各個級別的工程師能夠快速而滿懷信心地獲得答案。</p><p>Ansys Workbench 可支持與商業 CAD 工具的可靠連接，提供點擊按鈕設計點更新，流體和電氣求解器可提供無縫集成的多物理場功能。
展開 工業軟件研發中處理超大模型(6)--有限元求解器
有限元基本原理和特點
這里簡要介紹一下有限元方法基本原理和特點，有限元本質上是求解偏微分方程的一種方法。針對偏微分方程整體求解解析解困難的情況，有限元方法是將分析目標對象本身進行細化，劃分成小的有限個網格單元，比如三角形，四面體，六面體等，有時也會根據求解特點將高緯度幾何降為低緯度幾何單元，將構造的基函數和形函數應用在每一個網格單元，利用加權余量方法求解未知量。
有限元方法最早應用于結構力學分析，后來逐步推廣到熱，聲，電磁，流體等領域，是解決工程數值分析的通用方法。最早基于三角面片，后來推廣到各種形狀網格單元，主流網格單元為三角形，四邊形，四面體和六面體。目前大部分介紹有限元的教程集中在結構領域，對多物理場以及多物理場耦合介紹的非常少，這是國內教材可以加強的一個方向。
有限元方法自誕生以來，在理論基礎，物理，數學，工程應用等都得到了驗證，是目前求解偏微分方程以及多物理場問題最完備，最有效，最廣泛的數值計算方法。
需要提及的是我國的胡海昌于1954年提出了廣義變分原理，錢偉長最先研究了拉格朗日乘子法和廣義變分原理之間關系。而馮康更是在20世紀60年代于國內獨立提出了有限元的方法和概念，并研究了有限元分析的精度，邊界條件以及收斂性問題，指導了國內多項重大工程實踐。（1965年馮康發表了論文“基于變分原理的差分格式”，這篇論文是國際學術界承認我國獨立發展有限元方法的主要依據）。知Zienkiewicz而不知馮康，是我們科普和教材需要改進的地方。
關于有限元更詳細的介紹，點擊鏈接查看
一篇文章入門多物理場有限元（全篇）
2.
展開 有限元理論基礎及Abaqus內部實現方式研究系列24： 顯式求解Step By Step
（原創，轉載請注明出處）
==概述==
本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過
（1） 基礎理論
（2） 商軟操作
（3） 自編程序
三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。
有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。
一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。
展開 有限元理論基礎及Abaqus內部實現方式研究系列22： 幾何非線性的剛度矩陣求解
（原創，轉載請注明出處）
==概述==
本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過
（1） 基礎理論
（2） 商軟操作
（3） 自編程序
三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。
有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。
一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。
展開 有限元理論基礎及Abaqus內部實現方式研究系列4：非線性問題的求解 ￥1
我們不研究有限元的新方法、新理論，只是研究商用有限元軟件的實現方式。有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用軟件在這些傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。
一方面我們查閱Abaqus軟件手冊得到修正方法的說明，另一方面我們自己編程實現簡單的結構有限元求解器，通過自研求解器和Abaqus的結果比較結合理論手冊如同管中窺豹一般來研究Abaqus的修正方法，從而猜測商用有限元軟件的內部計算方法。在研究的同時，準備將自己的研究成果記錄下來寫成一個系列文章，希望對那些不僅僅滿足使用軟件，而想了解軟件內部實現方法甚至是做自己的軟件的朋友有些幫助。由于水平有限，里面可能有許多錯誤，歡迎交流討論。
2 第四篇：非線性問題的求解
隨著分析對象越來越復雜，有限元軟件從線性分析向非線性分析（如材料為非線性、幾何大變形導致的非線性、接觸行為引起的邊界條件非線性等）發展，Abaqus的強大主要就在于能解決各種工程上的非線性問題。由于非線性理論和編程實現的困難性，研究起來也相對困難，只能一步步來，我們在本文中首先研究各種非線性都會遇到的求解問題和Abaqus的內部實現方式。很多人做非線性編程結果正確性驗證由于不知道商軟的內部實現方式，采用小模型或者有理論解的模型對比理論結果和商軟最終的迭代結果，可能最終的結果差不多，但這種方法很容易就忽略了非線性中細節問題，如果進一步對比迭代過程，就會發現迭代效率和精度和商軟相比差距較大，這種效率和精度當遇到復雜工程問題時就尤為重要了。
展開 有限元理論基礎及Abaqus內部實現方式研究系列35： 接觸求解算法
（原創，轉載請注明出處）
1 概述
本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過
（1） 基礎理論
（2） 商軟操作
（3） 自編程序
三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。
有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。
一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。
展開 有限元理論基礎及Abaqus內部實現方式研究系列46：約束關系（2）-Lagrange因子法求解
（原創，轉載請注明出處）
1 概述
本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過
（1） 基礎理論
（2） 商軟操作
（3） 自編程序
三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。
有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。
一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元軟件iSolver，通過自研CAE軟件和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。iSolver包括完整的前后處理和有限元求解器，功能如下，有興趣可直接在下面網址下載：
百度網盤鏈接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取碼: 6fdf
2 約束關系
我們在系列文章第35章介紹了接觸分析及常用的基本算法。
展開 
Samcef Mecano 非線性柔性多體有限元求解器
求解器集合了有限元和多體動力系統的優勢，使用戶可以不用在分析一系統時轉化分析平臺。
附件介紹了Samcef求解器Mecano的功能特點和優勢。
Mecano.pdf
P Flyer_SAMCEF_MecanoTemporarily.pdf
面向對象有限元編程|自定義求解器之共軛梯度法
共軛梯度法是方程組求解的一種迭代方法。這種方法特別適合有限元求解，因為該方法要求系數矩陣為對稱正定矩陣，而有限元平衡方程的系數矩陣正好是對稱正定矩陣(考慮邊界條件)。同時，共軛梯度法也適合并行計算。
FlexPDE 5.0和FlexPDE 7.0（偏微分有限元求解器）英文版安裝及破解 ￥35
FlexPDE是基于有限元方法的偏微分方程（Partial Differential Equations,PDE）求解軟件。FlexPDE屬于商業軟件，其免費的學生版本（Student Version）對網格數量以及求解方程數量都做出了限制（雖然做出了限制，但仍遠優于同樣基于有限元方法的的matlab PDE工具箱）。
FlexPDE可以進行化學、地下水、污染物、地熱開發、油藏開發等相關的數學模型求解。
下圖為FlexPDE 5.0的主界面：
下圖為FlexPDE 7.0的主界面：
然而，我們通常所下載的免費版有網格數目限制。
因此，此處就是要給出破解方法和相關文檔。
下圖是用于進行油藏開發兩相滲流的計算運行圖。未破解前無法運行相關模擬，但破解后就沒有網格數目限制了。
網格模型：
產油動態曲線：
儲層壓力分布云圖：
展開 一款FEMAP求解器TP2000(由德國奔馳公司有限元法奠基人Prof.Peter Groth領導開發)，免費下載
改求解器也適用于MEDINA.
