FEM之多物理場求解器(2)---C++開發(fā)熱-結構耦合求解器(詳細)

2017-01-15 by:CAE仿真在線來源:互聯網

在FEM之多物理場求解器(1)---C++開發(fā)熱-結構耦合求解器一文中介紹了分別利用 3D熱求解器和結構求解器求解三維穩(wěn)態(tài)熱結構耦合問題,本文介紹熱求解器開發(fā)的更多細節(jié)。雖然熱和結構有限元理論已經很完善,而且單元矩陣和算法公式也已經推導出來,但要開發(fā)出能和商業(yè)軟件相同的求解器還是有很多工作要做。工欲善其事必先利其器,列出需要做的基礎準備工作:

1. 與Ansys 做benchmark:

1.1. 需要解析Ansys有限元模型,開發(fā)Ansys2Fas.exe 以解析 *.cdb文件;

1.2. 求解的結果要與Ansys比較,需要導出Ansys計算結果;

1.3. 原始的結果數據通常都經過處理,所以需要設計統一的結果文件格;

1.4. 開發(fā)后處理器用于查看結果;

1.5. 精度性能比較程序;

2. 為了將求解器產品化,需要開發(fā)前處理器,僅這一項又是不少的工作:

2.1. 在前處理器中要能把有限元模型導出為*.cdb格式。

2.2.Ansys除了簡單的單元類型還有許多擴展的單元類型,前處理器中需要支持。

2.3. 為了驗證Ansys的準確性,能導出第三方軟件格式,比如Abaqus,Nastran

3. 在求解器設計之初需要考慮到單元類型,邊界條件,荷載等擴展需求,定義好接口。

熱求解器

1. 支持功能:

2. 理論:

穩(wěn)態(tài)熱的偏微分方程相對來說簡單,在求解開發(fā)中有些要注意的問題:

2.1. 對于熱分析,線性四面體單元通常就能產生較高精度的解,但對于很多結構,線性四面體的剛度過大,需要使用二階四面體單元。在兩者耦合計算時就存在單元兼容的問題,雖然可以溫度采用線性單元,結構采用二階單元,精度上會有影響。

2.2. 當存在多種材料時,材料接觸部分要保持連續(xù),即幾何需要使用非二維流體,網格劃分也要支持此種類型。對前處理器提出了更高的要求。Nastran要求不同材料接觸部分網格要保持連續(xù),觀察過Ansys的接觸部分網格可以不連續(xù),姑且不討論這樣對不對,如果Ansys支持這樣的網格,起碼在接觸部分求解器需要做額外處理。

2.3. 公式

四面體和六面體都使用常規(guī)單元即可。矩陣推導過程與彈性力學完全,而且因為一個節(jié)點只有一個自由度,比經典彈性力學的剛度矩陣要簡單。熱傳導的最后計算公式為:

K*t+C*T'=P;

該公式中K為熱傳導剛度矩陣,t為溫度分布,C為儲熱剛度矩陣,T'為對時間導數,P為溫度荷載

當P為0時,變?yōu)樽杂蓽囟确植?

當T'為0時,變?yōu)榉€(wěn)態(tài)方程;

當P,T'都為0時候,就是最簡單的穩(wěn)態(tài)溫度分布;

3. Benchmark:

以下四個case基本上覆蓋了穩(wěn)態(tài)熱求解器功能。

3.1. Case1: 線性四面體,第一類邊界條件,兩種材料;

3.2.Case2: 二階四面體,第一類邊界條件,兩種材料;

3.3.Case3: 一階四面體,第二類邊界條件,兩種材料;

3.4.Case4: 二階四面體,第一,二類邊界條件,三種材料;

3.5. Case5. 二階四面體,復雜幾何;

4. 實現:

第一個case詳細講一下分析實現過程,后面的case都類似。

4.1 Case1:上下兩個接觸的長方體,材料不同,上部長方體上表面溫度100度,熱導率20;下部長方體下表面溫度0度,熱導率4000,采用線性四面體單元,求溫度場分布。

4.1.1.幾何:

Ansys的 WorkBench中進行前處理:

將工程導入到傳統的Ansys,將模型導出為*.cdb文件。在HyperMesh中打開*.cdb文件可以觀察到有限元模型實際上分成了5部分:

如上圖,本來只有2個模型,但是多出4個模型,用來處理接觸部分的模型。這也證明了前面的猜想,Ansys支持接觸模型不連續(xù)的網格,但是在接觸部分做了特殊處理。

4.1.2.Ansys的計算結果:

4.1.3.一般來講,復雜幾何的六面體網格比較難劃分,市面上提供六面體網格引擎的不多。對三維實體熱單元,Ansys中提供了二階的四面體單元SOLID70。由于Ansys對網格的接觸單元做了特殊處理,為了簡化模型,在自己開發(fā)的前處理器中將兩個長方體做為非二維流體處理,劃分二階四面體,Mesh結果如下: