差異六面體網格模型間數據映射傳遞研究

劉 安，胡廣旭，苗玉剛，鞏慶濤

（1.哈爾濱商業大學機械系，黑龍江哈爾濱 150028；2.哈爾濱工程大學船舶工程學院，黑龍江哈爾濱 150001；3.魯東大學蔚山船舶與海洋學院，山東煙臺 264025）

近些年來，由于制造行業對智能制造與數值模擬技術的廣泛認可，制造工藝鏈模擬技術也隨之誕生。制造工藝鏈模擬技術是工藝模擬技術向多工序制造流程的延伸。隨著工藝模擬研究的逐漸深入，采用多工藝制造的復雜零部件在不同工藝下的物理化學作用以及相互影響,漸漸受到工藝研究者的重視。歐洲先進制造領域研究者首先探索了制造工藝鏈模擬技術的相關研究。針對汽車結構制造優化問題，Papadakis L 等基于有限元方法開展了車體薄板結構沖壓、焊接以及強度與碰撞計算的連續模擬，系統地分析了各工藝環節對結構制造最終質量的影響，進而全面協調各工藝環節進行整體優化，提高了產品的綜合制造性[1]。Alexander Govik，Zaeh 等進一步開展了輕質合金的成形與焊接制造工藝鏈模擬研究[2-3]。在此基礎上，為了保證工藝模擬數據傳遞，Afazov 開發了用于不同有限元計算軟件之間的數據格式轉換系統FEDES[4]，并實現航空發動機圓盤組件的制造工藝鏈模擬，以保證應力應變數據在不同制造工藝模型間的繼承與傳遞，進而實現了熱處理、機械加工、噴丸處理等多種工序的制造工藝鏈模擬[5]。然而，盡管上述制造工藝鏈模擬技術可以“透視”整個制造過程，但實現制造工藝鏈模擬需有效地將有限元模擬數據在不同工藝模型間進行映射與傳遞，Afazov 相關程序代碼尚未公開和商業化，限制這項技術進一步普及應用。因此，國內連續模擬仍處于單工藝連續模擬階段，主要用于板料塑性成形工藝優化研究[6-7]。其中限制這項技術普及應用的關鍵問題是不同工藝模擬間的數據映射與傳遞技術，尤其是采用形函數法實現差異六面體網格間數據傳遞時，需求解非線性方程組，涉及迭代循環計算，計算數據量大，限制了該項技術的普及應用[8]。因此，文中將在前期研究基礎上[9-10]，通過直接計算方法代替迭代循環計算，實現了差異六面體網格的有效數據映射，并開發其用戶程序。

1 差異網格數據映射技術

開展制造工藝鏈模擬時，為了保證模擬物理數據在不同工藝階段數值模型的連續性,需將前道工序的模擬結果數據傳遞至后道工序的模型中。如果前后工序采用不同網格模型,還需將模擬結果數據映射至后道工序網格模型對應網格節點。常用差異網格間的映射插值技術有最近點法、點場距離法、單元距離法及形函數法[4]。其中,形函數法符合有限元網格原理,傳遞精度最高。該方法利用有限元單元形函數將單元節點的結果數據插值獲得單元內任意位置點的結果數據。

如式（1）所示，其中Dp為單元內任意點P的結果數據值（應力、應變或位移）。

Di為單元節點上的結果數據值，ηi為Di相應的形函數參量。

計算8 節點的3D 六面體網格內任意點P的數據值時，其形函數如式（2）[4]。

預求解六面體網格單元內任意點P對應映射局部坐標點(g,h,r)位置時,需求解式（3）的非線性方程組[4]。

此時計算式（3）時，需采用牛頓拉夫森進行循環迭代計算，且計算的精度取決于收斂精度。但采用上述方法實現三維六面體網格的差異網格數據傳遞時,求解局部坐標點(g,h,r)位置增加了網格傳遞循環計算量。

2 形函數法有限元數據傳遞算法流程

如圖1 所示網格模型A 的節點數據傳遞至模型B 時，其算法流程為：

圖1 數據映射的差異網格

1）在模型B 中循環獲取任意節點的全局坐標值，即為P(x,y,z)；

2）查找P(x,y,z)點在模型A 中所在單元EAi；

3）獲取單元EAi的節點坐標數據Nj(x,y,z)；

4）基于單元形函數，求解非線性方程組，如式（3），獲得P(x,y,z)點的局部坐標映射值P(g,h,r)；

5）在模型A 的模擬結果數據(應力、應變、溫度等)中查找單元EAi的節點或積分點結果數據；

6）采用形函數插值法，如式（2），通過單元EAi的節點模擬數據，計算獲取P(g,h,r)的插值映射值DP；

7）步驟循環1）～6），直至獲取模型B 中所有節點的模擬數據值，退出循環。

按照上述算法[11-13]，實現差異網格數據傳遞時，算法總循環為模型B 的總數CNB，傳遞其每一個節點數據時，都要循環獲取EAi，對于六面體網格為8節點，此時計算所需循環數為CEA×8；進行算法步驟4）時，求解非線性方程組進行牛頓拉夫森迭代時，仍需循環計算，循環數為nC，取決于迭代步長和收斂精度；獲取模型A 指定單元EAi模擬結果數據時，仍需進行循環A 單元數次εCEA，ε表示當獲取指定單元后剩余單元數無需循環的系數，0<ε<1，可近似為0.5。基于上述可知總算法循環數為式（4）：

假設模型A 六面體網格數量為10 萬，網格B 節點數量為7 萬，假設牛頓拉夫森循環數平均為50 次，此時總循環數按上式計算為2.803 5×1012次，且在每個循環下都進行大量空間計算，因此網格數據傳遞計算時間較長。數據傳遞的算法中，針對模型A 與B 的循環是不可避免的，因此為提高數據傳遞效率，本研究預實現步驟4）的簡化直接計算，進而使nC值為1 次，從總循環數將降為5.95×1010。

3 六面體網格分解直接計算法

數據傳遞算法步驟4）計算主要實現單元內P點對應的局部坐標系值P(g,h,r)。該局部坐標系是用于定義六面體網格形函數建立，如圖2 所示為局部坐標系定義[10，14-16]。該局部坐標系內，單元平面上局部坐標系值分別為-1≤g,h,r≤1。計算中,只要得到實際網格內任意點P(g,h,r)對應的局部坐標系點P(g,h,r)即可實現形函數方法映射。

圖2 六面體單元及局部坐標

為了避免求解式（3）非線性方程組，分解的比例系數法獲取P(g,h,r)值。首先進行如圖3 所示的單元分解，點O(xO,yO,zO)為當前單元中心點，即為：

圖3 單元分解示意圖

構建G(xg,yg,zg),H(xh,yh,zh),R(xr,yr,zr)，形成確定向量3 個方向向量以及面GOH，面HOR,面GOR，從而將單元劃分為8 個子域空間。根據向量與面GOH法向量nGOR的夾角θnGOR值可以判斷P點位于面GOR的正方向或負方向。

表1 P點子域判斷規則

根據P點所在子域位置，進一步計算P(g,h,r)。按照圖3 所示，此時，P點位于子域六面體Esub{O,H,M37,G,T,M23,N3,M43} 內。由于實際六面體網格為非理想等邊六面體,此時計算P點局部坐標值(g,h,r)時，g不僅僅是沿的投影，還需協調考慮相關向量。此時,為了近似求解(g,h,r)，該研究將的向量和平均值近似視為理想六面體的方向值，即為,如式（7）所示。

因此可以得到近似局部坐標系值如式（8）所示。

同理可得：

基于上述計算，即可獲得單元內任意點P對應的局部坐標系映射值(g,h,r)，實現了直接計算代替了迭代計算。

4 數據傳遞實例及分析

基于VC 對話框程序設計了用于網格數據傳遞的用戶界面，如圖4 所示。利用軟件用戶界面可以通過簡單操作實現網格數據映射傳遞。具體流程如下，首先讀取所需映射計算的A 與B 網格數據，然后讀取A 網格的模擬結果數據，進而在程序中建立映射所需的網格與模擬數據相關的數組與類信息，最后點擊“開始數據映射計算”按鈕實現數據映射計算并監控運算進度。

圖4 軟件用戶界面

如圖5 和圖6 所示分別為將孔板與T 型焊接結構的細網格A 模擬數據傳遞至粗網格B。由云圖結果可知，兩種結構網格A 與網格B 的應力分布基本一致，由表2 數據可知，孔板結構差異網格數比為ηk=12.7，而T 型焊接件的差異網格數比為ηt=62.6 。由圖7 和圖8 所示的Mise 應力曲線對比可知，對比范圍內數據傳遞獲得粗網格的應力值基本與細網格的值一致，僅在結構邊緣附近，由于網格尺寸差距，傳遞精度稍低，而在大部分位置傳遞數據基本與源數據一致，說明該研究采用直接計算法獲取局部坐標映射點(g,h,r)是有效的。

圖5 孔板結構Mise應力云圖

圖6 T型焊接結構Mise應力云圖

表2 網格數據傳遞信息

圖7 T型結構背側Mise應力隨縱向焊接路徑

圖8 孔板中心至邊緣不同距離的Mise應力曲線

5 結論

基于對單元形函數分析，設計并開發了實現差異六面體網格間模擬數據傳遞的程序，并成功實現了孔板與T 型焊構件由粗網格向細網格的模擬數據傳遞，對比數據傳遞前后結果可知，其應力值傳遞結果數據與源結果模擬數據基本一致，且其傳遞精度與差異網格的網格數比值無關。與此同時，在差異六面體網格數據映射計算過程中，為了避免迭代循環計算單元內任意點的映射局部坐標值，提出了單元分解直接計算代替迭代求解計算方法，從而顯著減少了網格數據傳遞的程序循環次數，提高了差異六面體網格間數據傳遞技術的可用性。