坯料隨機(jī)局部彎曲對(duì)滾彎成形結(jié)果影響的蒙特卡洛分析

　　摘 要：采用蒙特卡洛有限元模擬方法研究了坯料局部的隨機(jī)彎曲對(duì)滾彎成形結(jié)果的影響。為了提高滾彎模擬的效率，提出了基于歐拉網(wǎng)格和經(jīng)典梁?jiǎn)卧臐L彎模擬方案，并與傳統(tǒng)有限元模型和理論解對(duì)比驗(yàn)證了該方案的正確性。在此基礎(chǔ)上模擬了具有零均值正態(tài)分布的局部曲率的超長(zhǎng)坯料的滾彎過(guò)程，并統(tǒng)計(jì)產(chǎn)品曲率半徑的分布規(guī)律。結(jié)果表明：輸出曲率半徑分布近似滿足正態(tài)分布，且隨著坯料曲率標(biāo)準(zhǔn)差的增大，均值減小，方差增大，宏觀上表現(xiàn)為產(chǎn)品半徑減小。產(chǎn)品的目標(biāo)半徑越大，代表性單元長(zhǎng)度越長(zhǎng)，受初始彎曲的影響就越大;對(duì)于給定的目標(biāo)形狀，輥輪位置參數(shù)對(duì)實(shí)際輸出半徑的分布沒有影響。

　　鄭子君; 陶裕梅， 工程力學(xué) 發(fā)表時(shí)間：2021-08-02

　　關(guān)鍵詞：滾彎;歐拉網(wǎng)格;蒙特卡洛法;歐拉梁;隨機(jī)誤差

　　滾彎是對(duì)型材、板材進(jìn)行彎曲加工的一種常見工藝。以對(duì)稱式三輥滾彎?rùn)C(jī)為例，其工作方式為將坯料置于底輥和中輥之間，中輥以給定的下壓量壓緊坯料，通過(guò)輥輪或送料輥轉(zhuǎn)動(dòng)，帶動(dòng)坯料連續(xù)地通過(guò)滾彎?rùn)C(jī)并發(fā)生塑性彎曲。當(dāng)加工過(guò)程開始一段時(shí)間后，滾彎?rùn)C(jī)兩底輥之間的坯料將逐漸達(dá)到一個(gè)“定常”狀態(tài)，此時(shí)輥輪上的反力與流出出口的產(chǎn)品曲率不再隨著時(shí)間變化。

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

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

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

　　對(duì)坯料有隨機(jī)局部曲率的滾彎過(guò)程作蒙特卡洛模擬，可先建立整個(gè)坯料的模型并為其每個(gè)代表性單元長(zhǎng)度 (representative elementary length, REL，即曲率不發(fā)生明顯改變的長(zhǎng)度) 隨機(jī)生成初始曲率，然后進(jìn)行有限元模擬，最后分析流出滾彎?rùn)C(jī)的工件形狀。當(dāng)坯料足夠長(zhǎng)時(shí)，相當(dāng)于蒙特卡洛模擬的樣本數(shù)量足夠多。注意到本文考慮的是在局部 (單元級(jí)別) 的隨機(jī)初始曲率，模型參數(shù)的數(shù)量遠(yuǎn)遠(yuǎn)多于文獻(xiàn)中的情形，難以使用代理模型。而傳統(tǒng)的基于拉格朗日觀點(diǎn)的有限元模型進(jìn)行超長(zhǎng)坯料的滾彎模擬時(shí)，時(shí)間代價(jià)很高。這是因?yàn)椋阂环矫妫Ｐ托枰獮樵跐L彎?rùn)C(jī)外的坯料劃分網(wǎng)格，使得單元個(gè)數(shù)過(guò)多;另一方面，輥輪和坯料的接觸面積過(guò)小，使得接觸搜索算法效率不高;此外，初始構(gòu)型的計(jì)算也增大了建模的難度。

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

　　本文從計(jì)算流體力學(xué)受到啟發(fā)，提出基于歐拉觀點(diǎn)的滾彎模擬方案，并在經(jīng)典梁?jiǎn)卧夹g(shù)的基礎(chǔ)上，引入一個(gè)附加載荷項(xiàng)來(lái)處理材料在單元間流動(dòng)帶來(lái)的影響，從而提高有限元模擬的效率。在此基礎(chǔ)上采用蒙特卡洛法研究正態(tài)分布的局部初始曲率對(duì)對(duì)稱式三輥滾彎工藝的輸出形狀的影響。

　　1 基于歐拉網(wǎng)格的滾彎模擬方案

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

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

　　建立如圖 1 所示坐標(biāo)系。在歐拉觀點(diǎn)中，所有變量是一維空間坐標(biāo) 的函數(shù)。記 t 時(shí)刻 處的曲率為 ，將 處橫截面的幾何尺寸、材料屬性、初始曲率等截面參數(shù)合并記為向量 ，塑性內(nèi)變量集合記為向量 。根據(jù)歐拉梁的本構(gòu)關(guān)系，彎矩可以由變形程度、截面參數(shù)以及塑性內(nèi)變量確定，那么經(jīng)過(guò) 的時(shí)間，彎矩增量為：

　　設(shè)在 時(shí)刻，位于 處的材料截面由 t 時(shí)刻位于 處的材料截面被輥輪帶動(dòng)平移而來(lái) (圖 2(a)、 圖 2(b))， 于 是 有 。又當(dāng) 很短時(shí)，可認(rèn)為任意材料點(diǎn)的應(yīng)力路徑是簡(jiǎn)單的，于是計(jì)算 時(shí)刻 處的彎矩時(shí)只需知道 t 時(shí)刻該材料截面的塑性內(nèi)變量 ，于是式 (1) 變?yōu)椋?/p>

　　在小變形情形下， ?s 是與坐標(biāo) x 無(wú)關(guān)的常數(shù)，其物理意義為該時(shí)間步內(nèi)滾彎?rùn)C(jī)的進(jìn)料長(zhǎng)度。式 (2) 進(jìn)一步改寫為：

　　顯然，式 (3) 中第一項(xiàng)是由于空間 處的坯料構(gòu)型改變引起的，形式上可以寫為 ，其中為截面參數(shù)為 ，塑性內(nèi)變量為時(shí)的割線彎曲剛度。實(shí)際計(jì)算時(shí)， 需要用迭代法確定。

　　而式 (3) 第二項(xiàng)是 處構(gòu)型不變，只是材料發(fā)生流動(dòng)時(shí)引起的，類似于流體力學(xué)中的遷移項(xiàng)。該項(xiàng)涉及的所有變量在 時(shí)刻已經(jīng)求得，因此是已知的。記該項(xiàng)為 ，物理意義為 x 處的構(gòu)型不變而截面參數(shù)發(fā)生變化時(shí)的內(nèi)力不平衡量。則式 (3) 可寫為： kb (x)?κ (x)=?Mex (x)+ ?Ma (x) (4) 與材料坐標(biāo)下的經(jīng)典梁方程比較： kb (s)?κ (s)=?Mex (s) (5)

　　可以發(fā)現(xiàn)形式上僅僅相差了一項(xiàng) ，因此對(duì)式 (4) 進(jìn)行有限元離散時(shí)，只需要在拉格朗日列式的基礎(chǔ)上增加與內(nèi)力不平衡量 對(duì)應(yīng)的載荷項(xiàng)即可。如果采用經(jīng)典的 2 節(jié)點(diǎn) 4 變量 Hermite 插值，則式 (4) 的離散形式可寫為： Ke?u=?Fex + ?Fa (6) ?u = {?wI ,?θI ,?wJ ,?θJ } Ke ?Mex(x) ?Fex 其中：?jiǎn)卧灰葡蛄?由節(jié)點(diǎn)處的撓度和轉(zhuǎn)角的增量組成;單元割線剛度矩陣 和與 對(duì)應(yīng)的單元載荷向量 的公式與經(jīng)典梁?jiǎn)卧械耐茖?dǎo)完全相同，即[29]：

　　式中： 為單元局部坐標(biāo); 為形函數(shù)向量; 為單元長(zhǎng)度; 、 、 分別是分布力集度、作用在 處的集中力和作用在 處的集中力偶的增量。而與 對(duì)應(yīng)的單元附加載荷向量 為：

　　至此已得到式 (4) 的單元離散方程。在滾彎模擬中的每一時(shí)間步，總剛矩陣和載荷向量的組裝、約束的施加與非線性方程迭代求解過(guò)程均與常規(guī)的有限元方法基本相同，不同之處是該步開始前，需更新單元參數(shù)，即令 ，以模擬材料在滾彎?rùn)C(jī)內(nèi)的水平流動(dòng)。每一步模擬的原理性示意圖如圖 2 所示。

　　為了實(shí)現(xiàn)簡(jiǎn)單，實(shí)際計(jì)算時(shí)可取單元長(zhǎng)度 恰等于每步進(jìn)料長(zhǎng)度 ，并忽略輥輪的尺寸。此時(shí)具體的模擬步驟如下：

　　1) 初始化：將兩個(gè)底輥間的變形區(qū)空間等分成若干單元，初始化各單元高斯積分點(diǎn)處的塑性內(nèi)變量 和截面參數(shù) 。

　　2) 模擬調(diào)輥過(guò)程：約束底輥處的兩節(jié)點(diǎn)的垂向位移，為中輥下方節(jié)點(diǎn)施加指定的位移量，取附加載荷項(xiàng)為零，組裝式 (6) 并進(jìn)行非線性求解，更新塑性內(nèi)變量 。

　　3) 模擬材料流動(dòng)：從左到右各單元依次用右側(cè)相鄰單元的內(nèi)變量 和截面參數(shù) 覆蓋自身值 (圖 2(a)、圖 2(b));滾彎?rùn)C(jī)入口處的單元采用流入坯料的初始參數(shù);將出口處單元的數(shù)據(jù)寫入輸出文件 (圖 1)。

　　4) 平衡附加載荷：按式 (9) 計(jì)算單元內(nèi)附加載荷項(xiàng)，并組裝得到整體附加載荷向量;約束三個(gè)輥輪處的節(jié)點(diǎn)垂向位移，組裝式 (6) 并進(jìn)行非線性求解，更新變形區(qū)的構(gòu)型和各單元的塑性內(nèi)變量 (圖 2(c))。

　　5) 進(jìn)入下一時(shí)間步：返回第 3) 步繼續(xù)計(jì)算，直至總滾彎長(zhǎng)度達(dá)到預(yù)定目標(biāo)。

　　在進(jìn)行每一時(shí)間步的非線性方程的求解時(shí)，本文采用了割線迭代法更新 的方式[29]。求解方案采用 Mathematica 10.0 編程實(shí)現(xiàn)。

　　2 模擬方案的驗(yàn)證算例

　　為驗(yàn)證提出的基于歐拉網(wǎng)格的模擬方案，同時(shí)作為后續(xù)蒙特卡洛模擬的參照，考慮如下算例。

　　算例 1 中取對(duì)稱式三輥滾彎?rùn)C(jī)的底輥距，各輥直徑 ，中輥的下壓量為 ;坯料初始時(shí)完全平直，其橫截面為邊長(zhǎng) 的正方形，材料為理想彈塑性 ， 彈 性 模 量 為 ， 屈 服 強(qiáng) 度;總滾彎長(zhǎng)度為 。

　　使 用 本 文 方 案 時(shí) ， 底 輥 間 的 空 間 等 分 為 150 個(gè)單元;單元內(nèi)用兩點(diǎn)高斯積分，每個(gè)積分點(diǎn)在厚度方向上等分為 20 層;每步進(jìn)料長(zhǎng)度 mm (恰等于單元長(zhǎng)度)。進(jìn)行非線性求解時(shí)，收斂標(biāo)準(zhǔn)取為連續(xù)兩次迭代的位移差向量的二階范數(shù)，小于中輥下壓量的 10−7。作為對(duì)照，在 ANSYS 19.0 中建立平面應(yīng)力模型 (圖 3)，坯料長(zhǎng)度和厚度方向的網(wǎng)格尺寸分別為 和 ;設(shè)輥輪和坯料間粗糙接觸，左底輥主動(dòng)轉(zhuǎn)動(dòng)帶動(dòng)坯料進(jìn)入滾彎?rùn)C(jī)，其余輥輪被動(dòng)轉(zhuǎn)動(dòng);取動(dòng)態(tài)時(shí)間步長(zhǎng)上限使得坯料每步平移不超過(guò) 。此外，還通過(guò)曲率積分法得出最終定常狀態(tài)下的理論解作為參考。

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

　　中輥處的支座反力和滾彎?rùn)C(jī)出口處的坯料局部曲率隨著滾彎長(zhǎng)度變化關(guān)系如圖 4 所示。由于平面應(yīng)力模型的輸出曲率是用曲率圓擬合每 5 個(gè)相鄰節(jié)點(diǎn)的方式求得，這使曲率結(jié)果趨于平均，因此在滾彎長(zhǎng)度 處的峰值低于本文模型。除此之外本文方法和 ANSYS 平面模型得到的結(jié)果曲線幾乎重合。在滾彎開始時(shí)，受到直邊效應(yīng)的影響，曲線有較大的波動(dòng)，而當(dāng)滾彎長(zhǎng)度超過(guò)后，兩種模型的結(jié)果曲線均不再有明顯變化，指示均已基本達(dá)到定常狀態(tài)。這時(shí)兩種模型的結(jié)果均與曲率積分法得到的理論解吻合，誤差小于 2%。

　　提取定常狀態(tài)下產(chǎn)品橫截面上的殘余應(yīng)力分布如圖 5 所示。兩種有限元模型得到的殘余應(yīng)力均與理論解吻合良好，其中本文模型的吻合程度更高。這是因?yàn)楸疚哪Ｐ秃屠碚摻饩跉W拉 —伯努利梁理論，而 ANSYS 平面應(yīng)力模型并不遵循該理論的假設(shè)。

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

　　3 坯料有隨機(jī)局部彎曲時(shí)產(chǎn)品的曲率半徑分布

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

　　先考慮算例 2。當(dāng) 時(shí)，隨機(jī)生成的坯料初始曲率與由此重構(gòu)得的初始構(gòu)型如圖 6 (a) 所示。可以看到雖然坯料局部初始曲率分布范圍較大，但宏觀上看 200 m 長(zhǎng)度的坯料僅有 0.3 m 浪高，平均每米浪高僅 0.14 mm，在工程上屬于直線度很高的坯料，在通常的有限元模擬中可以當(dāng)作完全平直來(lái)處理。而用蒙特卡洛模擬計(jì)算輸出產(chǎn)品的各單元曲率和重構(gòu)得的構(gòu)型如圖 6 (b) 所示，各截面的輸出曲率雖然很不一致，但產(chǎn)品整體上仍可以用一個(gè)圓形很好地?cái)M合。由于該擬合圓的半徑就是工程實(shí)踐中觀測(cè)到的產(chǎn)品半徑，會(huì)直接影響到零件的安裝和配合，不妨稱其為宏觀曲率半徑;與之對(duì)應(yīng)，由各單元位移場(chǎng)求得的曲率半徑稱為局部曲率半徑;又稱當(dāng)坯料平直時(shí) ( ) 滾彎定常狀態(tài)得到的產(chǎn)品構(gòu)型為目標(biāo)形狀，其半徑為目標(biāo)曲率半徑。

　　圖 6 (b) 中顯示考慮了微小的隨機(jī)初始曲率后，宏觀曲率半徑略小于目標(biāo)曲率半徑，這并非偶然現(xiàn)象。圖 7 給出了不同偏差水平下的局部曲率半徑統(tǒng)計(jì)分布。可以看到局部曲率半徑基本符合正態(tài)分布，且隨著初始曲率偏差增大，擬合分布的均值變小，標(biāo)準(zhǔn)差近似線性增長(zhǎng);而宏觀曲率半徑則逐漸下降，偏離目標(biāo)曲率半徑。

　　采用相同的方法分析了其余算例，其偏差水平在 時(shí)的輸出曲率半徑分布如圖 8 所示，可見各算例的輸出半徑均近似滿足正態(tài)分布，都有宏觀半徑小于目標(biāo)半徑，但分布參數(shù)不同。

　　各算例的統(tǒng)計(jì)分布參數(shù)隨偏差水平的變化如圖 9 所示。從圖 9 (a) 看到相比曲率半徑的算術(shù)平均值，宏觀曲率半徑的下降趨勢(shì)更有規(guī)律，并且宏觀曲率半徑更符合工程使用的需要，更適合作為描述結(jié)果分布規(guī)律的位置參數(shù)。圖 9 (b) 顯示產(chǎn)品的曲率半徑標(biāo)準(zhǔn)差與坯料曲率的標(biāo)準(zhǔn)差基本成線性關(guān)系。圖 9 (c) 顯示當(dāng)坯料曲率標(biāo)準(zhǔn)差由小變大時(shí)，分布逐漸由負(fù)偏斜變?yōu)檎?峰度先減小后增大。在本文考慮的幾個(gè)算例中，偏度系數(shù)和峰度系數(shù)都很接近正態(tài)分布的對(duì)應(yīng)值 (0 和 3)，因此，輸出的曲率半徑分布可以用正態(tài)分布很好地近似。

　　對(duì)比算例 2 和算例 3 發(fā)現(xiàn)，兩者的目標(biāo)曲率半徑相近，此時(shí)在輥距相差 50% 的情況下，曲率半徑的分布規(guī)律和分布參數(shù)非常相似。這表明對(duì)于給定的目標(biāo)曲率半徑，局部輸出曲率的分布與輥輪位置參數(shù)幾乎無(wú)關(guān)。

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

　　算例 5 和算例 3 只有代表性單元長(zhǎng)度 REL 不同。這里 REL 表征隨機(jī)局部曲率的變化劇烈程度。當(dāng)坯料曲率的標(biāo)準(zhǔn)差不變時(shí)，該長(zhǎng)度越大，局部就越均勻，坯料的初始彎曲也就越明顯。對(duì)比結(jié)果發(fā)現(xiàn)，當(dāng) REL 增大時(shí)，產(chǎn)品曲率半徑的分布更加分散，而宏觀曲率半徑則沒有明顯變化。

　　4 結(jié)論

　　在用有限元模擬研究滾彎過(guò)程時(shí)，通常假定坯料是平直的或者具有一致的曲率。這種假設(shè)是否對(duì)模擬結(jié)果有偏向性的影響，則未見討論。為此本文在小曲率、細(xì)長(zhǎng)梁假設(shè)的基礎(chǔ)上，采用一種基于歐拉網(wǎng)格的滾彎模擬方案，使得超長(zhǎng)長(zhǎng)度的滾彎模擬成為可能，在此基礎(chǔ)上運(yùn)用蒙特卡洛方法，對(duì)坯料有零均值正態(tài)分布的隨機(jī)局部曲率的滾彎過(guò)程進(jìn)行了研究，結(jié)果發(fā)現(xiàn)：

　　(1) 產(chǎn)品宏觀曲率半徑隨著初始曲率偏差的增大而減小;這說(shuō)明要得到理想的產(chǎn)品形狀，坯料局部曲率波動(dòng)較大時(shí)，應(yīng)補(bǔ)償性地減小下壓量。

　　(2) 產(chǎn)品的局部曲率半徑大致服從正態(tài)分布，隨著初始曲率的標(biāo)準(zhǔn)差增大，輸出曲率半徑分布的標(biāo)準(zhǔn)差線性增加。

　　(3) 目標(biāo)曲率半徑相同時(shí)，輥輪位置參數(shù)對(duì)局部曲率半徑的分布沒有明顯影響。換言之，若不更換坯料，無(wú)法通過(guò)調(diào)整輥輪位置參數(shù)來(lái)使得輸出曲率更加集中到目標(biāo)曲率附近。

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