坯料隨機局部彎曲對滾彎成形結果影響的蒙特卡洛分析

鄭子君，陶裕梅

(重慶理工大學機械工程學院，重慶400040)

滾彎是對型材、板材進行彎曲加工的一種常見工藝。以對稱式三輥滾彎機為例，其工作方式為將坯料置于底輥和中輥之間，中輥以給定的下壓量壓緊坯料，通過輥輪或送料輥轉動，帶動坯料連續地通過滾彎機并發生塑性彎曲。當加工過程開始一段時間后，滾彎機兩底輥之間的坯料將逐漸達到一個“定?！睜顟B，此時輥輪上的反力與流出出口的產品曲率不再隨著時間變化。

由于滾彎工藝屬于無模成型，回彈量大，實踐中要得到目標曲率往往需要對工藝參數進行大量的試錯[1]。為減少工藝開發成本，已有大量文獻對坯料參數、工藝參數和產品構型之間的關系進行了探索。在理論研究方面，Basset 等[2]、Fu 等[3]基于三點彎曲模型推導了輥輪位置；Kim 等[4]、黃世軍等[1]提出采用圓弧來近似滾彎機內的坯料構型；王安恒等[5]提出在大曲率滾彎中應考慮中性層的偏移；劉志芳[6]、張子騫等[7]采用曲率積分的方式求解變形區內的坯料構型。在數值模擬方面，Fu 等[3]采用平面單元來模擬穩態滾彎，結果與三點彎模型得到的理論解吻合良好；Kim 等[4]采用板殼單元模擬薄板的滾彎，提出可按15%對輸出曲率進行修正；Feng 等[8]模擬了非對稱式三輥滾彎，并進行了實驗驗證；Ktari等[9]通過顯式動力學分析得到設計形狀和輥輪位置的關系圖譜；Shin 等[10]比較了平面單元和梁單元的模擬結果；Kagzi 等[11]對圓錐滾彎過程中非定常階段各輥輪受力的波動情況進行了仿真；Tran 等[12]比較了有限元模擬和實驗測量的板料表面應變演化規律；Groth 等[13]模擬了輥輪位置發生調整后，輸出曲率和工藝力的變化過程，并對調整速率提出了建議。

前述研究都假設坯料初始是平直的。實際生產中，由于制造精度的限制和偶然因素的影響，坯料雖然宏觀上整體平直，但若取很短的長度來觀測，則常常有隨機的局部曲率，這種隨機性在通常的研究中被忽略了。這樣的近似處理是否會對分析結果產生影響則未見討論。

要分析隨機參數對力學過程的影響，將有限元模擬與蒙特卡羅法結合起來是一種直觀而有效的辦法，其基本思路是：根據模型輸入參數的統計分布規律，隨機生成大量的有限元模型，求解后再對結果進行統計分析，從而得到輸出的分布規律。如朱健等[14]在碳纖維布加固的廠房結構中考慮了尺寸、材料強度、載荷的隨機誤差對抗震能力的影響；陳力波等[15]在簡支梁橋中考慮結構參數的隨機誤差，分析了其在地震中的易損性；金路等[16]分析了鋼架中梁柱的初始側移和直線度對整體剛度的影響；楊智勇等[17]對土質邊坡的多種失效模式進行了概率分析，指出次風險滑面也應予以重視；Rafiee 等[18]、Abebe 等[19]和吳永強[20]采用蒙特卡洛方法實現了機械制造的誤差六西格瑪分析和魯棒性設計。然而，基于有限元法的蒙特卡洛模擬，始終受到計算效率低的困擾[14?16,21?23]。通過簡化有限元模型[15]，建立代理模型[15,18?21]，優化抽樣方案[16]，引入矩方法近似計算[22]，以及利用模型對參數的梯度信息[23]等方式，可以一定程度上提高分析的效率。

對坯料有隨機局部曲率的滾彎過程作蒙特卡洛模擬，可先建立整個坯料的模型并為其每個代表性單元長度(representative elementary length, REL，即曲率不發生明顯改變的長度)隨機生成初始曲率，然后進行有限元模擬，最后分析流出滾彎機的工件形狀。當坯料足夠長時，相當于蒙特卡洛模擬的樣本數量足夠多。注意到本文考慮的是在局部(單元級別)的隨機初始曲率，模型參數的數量遠遠多于文獻中的情形，難以使用代理模型。而傳統的基于拉格朗日觀點的有限元模型進行超長坯料的滾彎模擬時，時間代價很高。這是因為：一方面，模型需要為在滾彎機外的坯料劃分網格，使得單元個數過多；另一方面，輥輪和坯料的接觸面積過小，使得接觸搜索算法效率不高；此外，初始構型的計算也增大了建模的難度。

近年大量文獻嘗試了采用歐拉或任意拉格朗日—歐拉網格進行金屬塑性加工的模擬，發現歐拉網格在接觸判斷、網格數量和反畸變方面具有優勢[24?28]。但歐拉網格在塑性加工中的應用主要集中在鍛造[25]、擠壓[26]、厚板軋制[27?28]等體積成型領域，采用三維或平面實體單元。而對于以彎曲變形為主的滾彎過程，顯然基于梁/板理論建立歐拉觀點的結構單元計算效率更高，但目前未見相關研究。

本文從計算流體力學受到啟發，提出基于歐拉觀點的滾彎模擬方案，并在經典梁單元技術的基礎上，引入一個附加載荷項來處理材料在單元間流動帶來的影響，從而提高有限元模擬的效率。在此基礎上采用蒙特卡洛法研究正態分布的局部初始曲率對對稱式三輥滾彎工藝的輸出形狀的影響。

1 基于歐拉網格的滾彎模擬方案

模擬塑性加工過程時，通常采用拉格朗日網格，這種網格會隨著材料點的運動而移動，便于追蹤材料點的受力路徑。然而，對于待加工坯料很長的滾彎過程，實際新增塑性變形的卻只有滾彎機內的一小部分，采用拉格朗日網格是不經濟的。若將滾彎機內的空間視為一維流場，左、右底輥視為流場的入口和出口，坯料在滾彎機內的撓度視為流場變量，則可以采用流體力學中常用的歐拉網格，在任意時刻只分析滾彎機內的部分坯料。此外，歐拉網格中節點沒有水平位移，坯料節點與輥輪的可能接觸位置是確定的，接觸搜索變得容易；在接觸點處，輥輪對坯料的推動作用可簡化為歐拉網格節點上的指定位移，用代數方法直接施加。

為簡單起見，沿用文獻[2?3,6]中對平面滾彎模型作的假設：在變形區內的坯料始終處于小撓度狀態；底輥間距顯著大于坯料橫截面的尺寸；重力、剪力和慣性力對坯料的變形影響可以忽略。此時坯料適用歐拉—伯努利直梁模型。

圖1 對稱式三輥滾彎工藝的示意圖Fig.1 Schematic sketch of symmetrical three-rolled bending

設在t+?t時刻，位于x處的材料截面由t時刻位于x+?s處的材料截面被輥輪帶動平移而來(圖2(a)、圖2(b))，于 是有C t+?t(x)=C t(x+?s)。又當 ?t很短時，可認為任意材料點的應力路徑是簡單的，于是計算t+?t時刻x處的彎矩時只需知道t時刻該材料截面的塑性內變量P t(x+?s)，于是式(1)變為：

圖2 基于歐拉網格的滾彎模擬原理示意圖Fig.2 Schematic sketch of Eulerian-grid-based simulation of roll bending process

可以發現形式上僅僅相差了一項?Ma，因此對式(4)進行有限元離散時，只需要在拉格朗日列式的基礎上增加與內力不平衡量?Ma對應的載荷項即可。如果采用經典的2節點4變量Hermite 插值，則式(4)的離散形式可寫為：

為了實現簡單，實際計算時可取單元長度l恰等于每步進料長度 ?s，并忽略輥輪的尺寸。此時具體的模擬步驟如下：

1)初始化：將兩個底輥間的變形區空間等分成若干單元，初始化各單元高斯積分點處的塑性內變量P和截面參數C。

2)模擬調輥過程：約束底輥處的兩節點的垂向位移，為中輥下方節點施加指定的位移量，取附加載荷項為零，組裝式(6)并進行非線性求解，更新塑性內變量P。

3)模擬材料流動：從左到右各單元依次用右側相鄰單元的內變量P和截面參數C覆蓋自身值(圖2(a)、圖2(b))；滾彎機入口處的單元采用流入坯料的初始參數；將出口處單元的數據寫入輸出文件(圖1)。

4)平衡附加載荷：按式(9)計算單元內附加載荷項，并組裝得到整體附加載荷向量；約束三個輥輪處的節點垂向位移，組裝式(6)并進行非線性求解，更新變形區的構型和各單元的塑性內變量P(圖2(c))。

5)進入下一時間步：返回第3)步繼續計算，直至總滾彎長度達到預定目標。

在進行每一時間步的非線性方程的求解時，本文采用了割線迭代法更新kb的方式[29]。求解方案采用Mathematica 10.0編程實現。

2 模擬方案的驗證算例

為驗證提出的基于歐拉網格的模擬方案，同時作為后續蒙特卡洛模擬的參照，考慮如下算例。

算例1中取對稱式三輥滾彎機的底輥距L=300 mm，各輥直徑Dt=Db=50 mm，中輥的下壓量為d=2.7 mm；坯料初始時完全平直，其橫截面為邊長10 mm的正方形，材料為理想彈塑性，彈性模量為E=200 GPa ，屈服強度YS=200 MPa ；總滾彎長度為600 mm。

使用本文方案時，底輥間的空間等分為150個單元；單元內用兩點高斯積分，每個積分點在厚度方向上等分為20層；每步進料長度?s=2 mm(恰等于單元長度)。進行非線性求解時，收斂標準取為連續兩次迭代的位移差向量的二階范數，小于中輥下壓量的10?7。作為對照，在ANSYS 19.0中建立平面應力模型(圖3)，坯料長度和厚度方向的網格尺寸分別為2 mm 和1 mm；設輥輪和坯料間粗糙接觸，左底輥主動轉動帶動坯料進入滾彎機，其余輥輪被動轉動；取動態時間步長上限使得坯料每步平移不超過2 mm。此外，還通過曲率積分法得出最終定常狀態下的理論解作為參考。

圖3 ANSYS中的滾彎平面應力模型Fig.3 Plane stress model of roll bending process in ANSYS

算例在CPU i7-6820HQ，內存32 GB的工作站上以串行方式計算。ANSYS平面應力模型的模擬耗時近10 h，而本文方法僅耗時10 min，計算效率提升明顯。

中輥處的支座反力和滾彎機出口處的坯料局部曲率隨著滾彎長度變化關系如圖4所示。由于平面應力模型的輸出曲率是用曲率圓擬合每5個相鄰節點的方式求得，這使曲率結果趨于平均，因此在滾彎長度150 mm處的峰值低于本文模型。除此之外本文方法和ANSYS平面模型得到的結果曲線幾乎重合。在滾彎開始時，受到直邊效應的影響，曲線有較大的波動，而當滾彎長度超過400 mm后，兩種模型的結果曲線均不再有明顯變化，指示均已基本達到定常狀態。這時兩種模型的結果均與曲率積分法得到的理論解吻合，誤差小于2%。

圖4 中輥反力和輸出曲率隨滾彎長度的變化Fig.4 Top roller reaction forceand output curvature versus roll bending length

提取定常狀態下產品橫截面上的殘余應力分布如圖5所示。兩種有限元模型得到的殘余應力均與理論解吻合良好，其中本文模型的吻合程度更高。這是因為本文模型和理論解均基于歐拉—伯努利梁理論，而ANSYS平面應力模型并不遵循該理論的假設。

圖5 定常狀態下輸出坯料的殘余應力在厚度方向上的分布Fig.5 Distribution of residual stress through the thickness of the crosssection in steady state

此外，由圖4還可以看出ANSYS模型的輥輪反力和輸出曲率曲線始終有微小波動，這主要是由于接觸算法的數值誤差引起的。這種誤差具有一定的隨機性，因此可能在蒙特卡洛模擬中與坯料局部曲率導致的結果波動相混淆；但要縮小接觸誤差將導致接觸算法收斂困難。而本文模型能從原理上避免這種誤差，輸出曲線非常光滑，因此更適于進行蒙特卡洛模擬。

3 坯料有隨機局部彎曲時產品的曲率半徑分布

機械制造中，包括初始曲率半徑在內的坯料誤差常被認為服從正態分布。那么運用本文方法進行滾彎的蒙特卡洛模擬時，只需在前述模擬流程的基礎上，在每個時間步為入口處的REL生成正態分布的偽隨機初始曲率κ0即可。為研究模型參數對產品的曲率半徑分布的影響，在算例1的坯料模型參數和網格設置的基礎上，選取了4組不同的輥輪位置和REL 長度，并將四種情形分別記為算例2～算例5，如表1所示。對每個算例，考慮了4種偏差水平(σκ0=0.005 m?1～0.02 m?1)的零均值正態分布隨機初始曲率。對每種偏差水平，進行滾彎長度為2 00 m的蒙特卡洛模擬(即樣本數量為5萬～10萬個)，并在統計時去掉初始的1000個樣本以避免直邊效應的影響。

表1 各數值算例的參數設置Table 1 The parameter sets of the numerical examples

先考慮算例2。當σκ0=0.02 m?1時，隨機生成的坯料初始曲率與由此重構得的初始構型如圖6(a)所示?？梢钥吹诫m然坯料局部初始曲率分布范圍較大，但宏觀上看200 m 長度的坯料僅有0.3 m 浪高，平均每米浪高僅0.14 mm，在工程上屬于直線度很高的坯料，在通常的有限元模擬中可以當作完全平直來處理。而用蒙特卡洛模擬計算輸出產品的各單元曲率和重構得的構型如圖6(b)所示，各截面的輸出曲率雖然很不一致，但產品整體上仍可以用一個圓形很好地擬合。由于該擬合圓的半徑就是工程實踐中觀測到的產品半徑，會直接影響到零件的安裝和配合，不妨稱其為宏觀曲率半徑；與之對應，由各單元位移場求得的曲率半徑稱為局部曲率半徑；又稱當坯料平直時(σκ0=0)滾彎定常狀態得到的產品構型為目標形狀，其半徑為目標曲率半徑。

圖6(b)中顯示考慮了微小的隨機初始曲率后，宏觀曲率半徑略小于目標曲率半徑，這并非偶然現象。圖7給出了不同偏差水平下的局部曲率半徑統計分布?？梢钥吹骄植壳拾霃交痉险龖B分布，且隨著初始曲率偏差增大，擬合分布的均值變小，標準差近似線性增長；而宏觀曲率半徑則逐漸下降，偏離目標曲率半徑。

圖6 當 σκ0=0.02 m?1時算例2的工件形狀Fig.6 Theconfigurations ofthework-piece in Case2when σκ0=0.02m?1

圖7 算例2產品的局部曲率半徑分布Fig.7 The distributionsof theoutput local radiusin Case 2

采用相同的方法分析了其余算例，其偏差水平在σκ0=0.02 m?1時的輸出曲率半徑分布如圖8所示，可見各算例的輸出半徑均近似滿足正態分布，都有宏觀半徑小于目標半徑，但分布參數不同。

圖8 當 σκ0=0.02 m?1時算例3～算例5的局部曲率半徑分布Fig.8 The distributions of the output local radius in Case 3～Case 5 when σκ0=0.02 m?1

各算例的統計分布參數隨偏差水平的變化如圖9所示。從圖9(a)看到相比曲率半徑的算術平均值，宏觀曲率半徑的下降趨勢更有規律，并且宏觀曲率半徑更符合工程使用的需要，更適合作為描述結果分布規律的位置參數。圖9(b)顯示產品的曲率半徑標準差與坯料曲率的標準差基本成線性關系。圖9(c)顯示當坯料曲率標準差由小變大時，分布逐漸由負偏斜變為正偏斜；峰度先減小后增大。在本文考慮的幾個算例中，偏度系數和峰度系數都很接近正態分布的對應值(0和3)，因此，輸出的曲率半徑分布可以用正態分布很好地近似。

圖9 產品的局部曲率半徑分布參數隨坯料曲率標準差的變化規律Fig.9 Distribution parametersof output local radiusversus standard deviation of input curvature

對比算例2和算例3發現，兩者的目標曲率半徑相近，此時在輥距相差50%的情況下，曲率半徑的分布規律和分布參數非常相似。這表明對于給定的目標曲率半徑，局部輸出曲率的分布與輥輪位置參數幾乎無關。

而算例4具有和算例2相同的跨度，但目標曲率半徑增大了近30%，結果宏觀曲率半徑下降幅度增大了近1倍，局部曲率半徑的分布更加分散，標準差增大近90%。這表明隨著設計目標曲率半徑增大，產品受到坯料初始隨機曲率的影響逐漸變大。

算例5和算例3只有代表性單元長度REL不同。這里REL 表征隨機局部曲率的變化劇烈程度。當坯料曲率的標準差不變時，該長度越大，局部就越均勻，坯料的初始彎曲也就越明顯。對比結果發現，當REL 增大時，產品曲率半徑的分布更加分散，而宏觀曲率半徑則沒有明顯變化。

4 結論

在用有限元模擬研究滾彎過程時，通常假定坯料是平直的或者具有一致的曲率。這種假設是否對模擬結果有偏向性的影響，則未見討論。為此本文在小曲率、細長梁假設的基礎上，采用一種基于歐拉網格的滾彎模擬方案，使得超長長度的滾彎模擬成為可能，在此基礎上運用蒙特卡洛方法，對坯料有零均值正態分布的隨機局部曲率的滾彎過程進行了研究，結果發現：

(1)產品宏觀曲率半徑隨著初始曲率偏差的增大而減?。贿@說明要得到理想的產品形狀，坯料局部曲率波動較大時，應補償性地減小下壓量。

(2)產品的局部曲率半徑大致服從正態分布，隨著初始曲率的標準差增大，輸出曲率半徑分布的標準差線性增加。

(3)目標曲率半徑相同時，輥輪位置參數對局部曲率半徑的分布沒有明顯影響。換言之，若不更換坯料，無法通過調整輥輪位置參數來使得輸出曲率更加集中到目標曲率附近。

(4)想要得到的產品半徑越大，代表性單元長度越長，滾彎結果越容易受到坯料局部初始彎曲的影響。