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

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

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

1 差異網(wǎng)格數(shù)據(jù)映射技術(shù)

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

如式（1）所示，其中Dp為單元內(nèi)任意點P的結(jié)果數(shù)據(jù)值（應(yīng)力、應(yīng)變或位移）。

Di為單元節(jié)點上的結(jié)果數(shù)據(jù)值，ηi為Di相應(yīng)的形函數(shù)參量。

計算8 節(jié)點的3D 六面體網(wǎng)格內(nèi)任意點P的數(shù)據(jù)值時，其形函數(shù)如式（2）[4]。

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

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

2 形函數(shù)法有限元數(shù)據(jù)傳遞算法流程

如圖1 所示網(wǎng)格模型A 的節(jié)點數(shù)據(jù)傳遞至模型B 時，其算法流程為：

圖1 數(shù)據(jù)映射的差異網(wǎng)格

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

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

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

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

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

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

7）步驟循環(huán)1）～6），直至獲取模型B 中所有節(jié)點的模擬數(shù)據(jù)值，退出循環(huán)。

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

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

3 六面體網(wǎng)格分解直接計算法

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

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

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

圖3 單元分解示意圖

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

表1 P點子域判斷規(guī)則

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

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

同理可得：

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

4 數(shù)據(jù)傳遞實例及分析

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

圖4 軟件用戶界面

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

圖5 孔板結(jié)構(gòu)Mise應(yīng)力云圖

圖6 T型焊接結(jié)構(gòu)Mise應(yīng)力云圖

表2 網(wǎng)格數(shù)據(jù)傳遞信息

圖7 T型結(jié)構(gòu)背側(cè)Mise應(yīng)力隨縱向焊接路徑

圖8 孔板中心至邊緣不同距離的Mise應(yīng)力曲線

5 結(jié)論

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