模擬滾彎非穩(wěn)態(tài)過程的虛擬載荷法

鄭子君，王 洪

(重慶理工大學(xué) 機(jī)械工程學(xué)院，重慶 400040)

1 引 言

滾彎是一種常見的型材和板料彎曲加工工藝。其涉及的機(jī)械結(jié)構(gòu)簡(jiǎn)單，便于連續(xù)生產(chǎn)，使用一套輥輪就可以加工多種曲率的產(chǎn)品，因此有著獨(dú)特的優(yōu)勢(shì)。以平面對(duì)稱式三輥滾彎為例，其基本工藝流程為[1,2],將工件一段置于滾彎?rùn)C(jī)中(圖1(a))，在靜態(tài)彎曲階段，調(diào)整輥輪至設(shè)計(jì)位置，夾緊工件(圖1(b))，然后在動(dòng)態(tài)彎曲階段，通過送料輥轉(zhuǎn)動(dòng)，帶動(dòng)工件在水平方向連續(xù)地通過滾彎?rùn)C(jī)并發(fā)生塑性彎曲(圖1(c))。經(jīng)過一段時(shí)間后，滾彎?rùn)C(jī)內(nèi)各空間點(diǎn)處的工件構(gòu)型將達(dá)到一個(gè)穩(wěn)定狀態(tài)，輥輪反力與輸出曲率均趨于定值。

對(duì)于滾彎的穩(wěn)定狀態(tài)的研究，已有較多的理論模型，如三點(diǎn)彎模型[3，4]、相切圓模型[5，6]、曲梁模型[7]及曲率積分法[8,9]等。此外，有限元模擬也是重要的研究手段，如Feng等[10]模擬了非對(duì)稱式三輥滾彎，并進(jìn)行了實(shí)驗(yàn)驗(yàn)證；Ktari等[11]通過顯式動(dòng)力學(xué)模擬得到了穩(wěn)態(tài)輸出曲率和輥輪位置的關(guān)系；Fu等[4]采用平面應(yīng)變單元來(lái)模擬穩(wěn)態(tài)滾彎，在小變形情形下結(jié)果與三點(diǎn)彎模型的理論解吻合良好；Shin等[12]比較了平面單元和梁?jiǎn)卧哪M結(jié)果，指出平面單元表現(xiàn)更好；Kim等[5]采用板殼單元模擬了薄板滾彎，提出采用這種單元時(shí)可按15%對(duì)輸出曲率進(jìn)行修正。

對(duì)滾彎非穩(wěn)態(tài)過程的研究，有助于控制直邊效應(yīng)和提高產(chǎn)品質(zhì)量，主要通過有限元模擬進(jìn)行，如Kagzi等[2]對(duì)圓錐滾彎過程中各輥輪受力的波動(dòng)情況進(jìn)行了仿真；Tran等[13]比較了有限元模擬和實(shí)際測(cè)量的板料表面應(yīng)變演化規(guī)律；Groth等[14,15]模擬了輥輪位置發(fā)生調(diào)整后，輸出曲率和工藝力出現(xiàn)波動(dòng)并最終達(dá)到穩(wěn)態(tài)的過程。

圖1 對(duì)稱式三輥滾彎工藝

上述研究在模擬滾彎時(shí)往往采用相似的有限元建模方案，即工件采用柔性體建模，輥輪采用剛體建模；工件和輥輪間定義接觸；通過指定輥輪的運(yùn)動(dòng)，帶動(dòng)工件產(chǎn)生位移和塑性變形。然而實(shí)踐中發(fā)現(xiàn)，這種方案的效率不高，通常明顯慢于同等尺度的有模彎曲的模擬。原因主要有兩點(diǎn)，(1) 滾彎過程中工件不斷平移通過變形區(qū)(圖 1(a)AC間)，必須在變形區(qū)外預(yù)留足夠長(zhǎng)的工件，導(dǎo)致網(wǎng)格數(shù)量過多；(2) 由于輥輪與工件間接觸面積小，且有切向相對(duì)滑動(dòng)，很容易發(fā)生穿透和接觸狀態(tài)突變，因此處理接觸時(shí)需采用小時(shí)間步長(zhǎng)、大接觸剛度和大檢索半徑，降低了接觸算法的效率[16,17]。

受計(jì)算流體力學(xué)的啟發(fā)，本文基于歐拉觀點(diǎn)描述三輥滾彎?rùn)C(jī)變形區(qū)內(nèi)的型材工件狀態(tài)變化。把滾彎?rùn)C(jī)內(nèi)的變形區(qū)視作一維流場(chǎng)，該空間內(nèi)的工件撓度和轉(zhuǎn)角視為流場(chǎng)變量。其優(yōu)點(diǎn)在于，工件未進(jìn)入滾彎?rùn)C(jī)的部分不需要建立網(wǎng)格，只需提供其在出入口處的狀態(tài)；網(wǎng)格的拓?fù)湫螒B(tài)和結(jié)點(diǎn)坐標(biāo)均是固定的，使得結(jié)點(diǎn)和輥輪間的接觸搜索得到簡(jiǎn)化。

然而，采用歐拉觀點(diǎn)時(shí)，單元方程形式往往和通常的拉格朗日列式有明顯的區(qū)別，增加了開發(fā)和使用的成本。因此本文針對(duì)小曲率滾彎工藝中工件移動(dòng)的特點(diǎn)，提出一個(gè)虛擬載荷的概念，使得基于歐拉觀點(diǎn)的滾彎控制方程可以直接套用較為成熟的拉格朗日梁?jiǎn)卧惺健?/p>

2 虛擬載荷法

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

建立如圖1所示坐標(biāo)系。記t時(shí)刻某空間坐標(biāo)x處的曲率為κt(x)，材料和截面參數(shù)集合為向量Ct(x)，塑性內(nèi)變量集合為向量Pt(x)。根據(jù)歐拉梁的本構(gòu)關(guān)系，彎矩由變形程度、截面參數(shù)以及塑性內(nèi)變量確定。由內(nèi)外力增量的平衡得

ΔMe x(x)=M[κt + Δ t(x),Ct + Δ t(x),Pt + Δ t(x)]-

M[κt(x),Ct(x),Pt(x)]

(1)

設(shè)經(jīng)過微小時(shí)間步長(zhǎng)Δt，位于x處的材料截面由原本位于x+s處的材料截面受輥輪帶動(dòng)平移而來(lái)，即有Ct + Δt(x)=Ct(x+s)；又由步長(zhǎng)Δt很小，可以假設(shè)任意材料點(diǎn)的應(yīng)力路徑是簡(jiǎn)單的，則計(jì)算彎矩時(shí)，只需要知道該時(shí)間步內(nèi)材料截面的初始塑性內(nèi)變量值Pt(x+s)。于是式(1)變?yōu)?/p>

ΔMe x(x)=M[κt + Δt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x),Pt(x)]

(2)

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

ΔMe x(x)={M[κt + Δt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x+s),Pt(x+s)]}+

{M[κt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x),Pt(x)]}

(3)

式(3)第二項(xiàng)的意義是材料水平流動(dòng)后，為維持t時(shí)刻構(gòu)型不變所需要的額外載荷，記為-ΔMa r t(x)。而式(3)第一項(xiàng)的物理意義是由變形引起的彎矩變化。若定義割線彎曲剛度ks為

ks(x)={M[κt + Δt(x),Ct + Δt(x),Pt + Δt(x)]-

M[κt(x),Ct + Δt(x),Pt + Δt(x)]}/

[κt + Δt(x)-κt(x)]

(4)

則平衡方程(3)可寫為

ks(x)Δκ(x)=ΔMe x(x)+ΔMa r t(x)

(5)

可以看出，式(5)雖然基于空間坐標(biāo)x得出，但與基于材料坐標(biāo)的歐拉梁方程相比，差別僅是右端項(xiàng)除了實(shí)際彎矩增量ΔMe x(x)以外，還額外多了一項(xiàng)ΔMa rt(x)。它是空間構(gòu)型不變時(shí)，由于材料流動(dòng)而導(dǎo)致的內(nèi)力不平衡量。由于該項(xiàng)在t+Δt時(shí)刻是已知的，故可以作為載荷項(xiàng)來(lái)處理。引入額外載荷項(xiàng)來(lái)考慮材料流動(dòng)影響的方法可稱為虛擬載荷法。從定義容易看出，如果梁的截面參數(shù)C和塑性內(nèi)變量P在各處是均勻的，則虛擬載荷項(xiàng)自動(dòng)消失。

3 有限元實(shí)施方案

考慮用有限元方法求解控制方程(5)。

(1) 網(wǎng)格方案。將兩個(gè)底輥間的變形區(qū)空間(圖 1中AC間)等分成若干單元。由于有限元程序通常只在高斯積分點(diǎn)處存儲(chǔ)截面數(shù)據(jù)，為了提高精度，取工件每步進(jìn)料長(zhǎng)度s恰等于單元的尺寸l。則t+ Δt時(shí)刻各單元積分點(diǎn)的參數(shù)即為其上游單元積分點(diǎn)在t時(shí)刻的參數(shù)；而入口處單元的參數(shù)為從進(jìn)口新流入工件的初始參數(shù)。

(2) 截面參數(shù)集合與本構(gòu)關(guān)系。簡(jiǎn)單起見，材料采用雙線性隨動(dòng)強(qiáng)化彈塑性模型。把橫截面在厚度上分為若干層，則高斯積分點(diǎn)處的截面參數(shù)集合C包括初始曲率κ0，楊氏模量E，屈服強(qiáng)度σs，硬化模量Et，各層形心到中性軸的距離ym(m為層號(hào))和慣性矩Im；塑性內(nèi)變量集合P則包含各層形心處的軸向塑性應(yīng)變?chǔ)舙 m。則彎矩與曲率和截面參數(shù)的關(guān)系為

(6)

式(6)在線彈性時(shí)退化為彈性歐拉梁方程M=EIκ。

(3) 單元平衡方程。原則上而言，虛擬載荷法的單元列式可基于現(xiàn)有的任意梁?jiǎn)卧匠掏茖?dǎo)[18-20]。簡(jiǎn)單起見，本文單元位移場(chǎng)采用兩結(jié)點(diǎn)三次Hermite插值，即取撓度w和轉(zhuǎn)角θ為結(jié)點(diǎn)變量，形函數(shù)N={1-3ξ2+2ξ3,(ξ-2ξ2+ξ3)l,3ξ2-2ξ3,(ξ3-ξ2)l}，則單元列式推導(dǎo)與經(jīng)典梁?jiǎn)卧猍20]完全相同。因此本文不加推導(dǎo)地給出單元平衡方程的增量格式：

(7)

(8)

式中ξg和Wg為積分域?yàn)?0,1)的第g個(gè)高斯積分點(diǎn)和對(duì)應(yīng)的權(quán)重。與結(jié)點(diǎn)外力載荷對(duì)應(yīng)的右端項(xiàng)為

(9)

式中 Δq，Δf和Δm分別為單元受分布力、集中力和集中力偶的增量。而與虛擬載荷對(duì)應(yīng)的右端項(xiàng)為

{M[κt(ξg),Ct(ξg),Pt(ξg)]-

(10)

(4) 方程的組裝和約束。組裝與約束過程和傳統(tǒng)有限元方法完全相同。由于歐拉網(wǎng)格中結(jié)點(diǎn)不會(huì)隨構(gòu)型改變而移動(dòng)，故接觸處理可以簡(jiǎn)化。更進(jìn)一步，如果忽略輥輪尺寸，即將模型簡(jiǎn)化為三點(diǎn)彎[3,4]，則模具對(duì)工件的作用簡(jiǎn)化為結(jié)點(diǎn)強(qiáng)制位移條件，可以通過矩陣代數(shù)變形直接施加[20]。

(5) 非線性方程的求解。當(dāng)工件發(fā)生彈塑性變形時(shí)，式(8)的割線彎曲剛度不是常數(shù)。本文采用割線迭代法，對(duì)割線剛度進(jìn)行迭代更新求解，直至位移增量收斂。迭代結(jié)束后，更新各高斯積分點(diǎn)處的塑性內(nèi)變量Pt +Δt(ξg)，并進(jìn)入下一時(shí)間步的求解。

在實(shí)際計(jì)算中，若遇到割線法迭代收斂困難，可以將每次平移產(chǎn)生的虛擬載荷分多個(gè)子步施加，并動(dòng)態(tài)調(diào)整子步步長(zhǎng)。

(6) 程序流程圖。用虛擬載荷法模擬滾彎過程的流程如圖2所示。

圖2 采用虛擬載荷法模擬滾彎流程

4 算 例

取對(duì)稱式三輥滾彎?rùn)C(jī)的輥距L=400 mm，輥徑Dt=Db=50 mm。中輥下壓量固定為d=3.8 mm。工件的橫截面為h=10 mm的矩形，初始時(shí)工件平直。材料采用理想彈塑性模型，楊氏模量E=200 GPa，屈服強(qiáng)度σs=200 MPa?？倽L彎長(zhǎng)度為600 mm。

采用Mathematica 11.0編寫虛擬載荷法程序；底輥AC間的空間劃分為 200個(gè)單元；單元內(nèi)用兩點(diǎn)高斯積分，每個(gè)積分點(diǎn)在厚度方向上等分為20層；動(dòng)態(tài)彎曲階段每步進(jìn)料長(zhǎng)度s=2 mm。作為對(duì)照，在ANSYS 18.0中分別用平面應(yīng)力單元和板殼單元建立兩個(gè)有限元模型，網(wǎng)格尺寸完全對(duì)照虛擬載荷法的模型；輥輪和工件間粗糙接觸，左底輥輪主動(dòng)轉(zhuǎn)動(dòng)帶動(dòng)工件進(jìn)入滾彎?rùn)C(jī)，其余輥輪被動(dòng)轉(zhuǎn)動(dòng)；設(shè)置動(dòng)態(tài)時(shí)間步長(zhǎng)上限使工件每次平移不超過2 mm。此外，還通過曲率積分法[10]得出靜態(tài)彎曲狀態(tài)和最終穩(wěn)定狀態(tài)下的理論解作為參考。

各模擬方案均在CPU為i7-6830HQ，內(nèi)存為16 GB的計(jì)算機(jī)上以串行模式進(jìn)行。ANSYS平面應(yīng)力模型耗時(shí)約575 min，板殼模型耗時(shí)約 61 min，而虛擬載荷法僅耗時(shí)15 min，明顯快于傳統(tǒng)方案。

求解得到的中輥的約束反力和出口處的曲率隨滾彎長(zhǎng)度S的變化如圖3所示?？梢钥闯龈髂Ｐ偷慕Y(jié)果曲線非常相似。大約在滾過1.5倍跨長(zhǎng)后，各解基本達(dá)到穩(wěn)態(tài)。對(duì)比曲率積分法給出的穩(wěn)態(tài)解可知，虛擬載荷法與理論解最接近。板殼模型與其他模型相比，變形量明顯偏小，下壓力偏大，與文獻(xiàn)[11,12]結(jié)論相符，這是因?yàn)榘鍤つＰ椭写嬖趶澢嫱獾恼龖?yīng)力，而其他模型均未計(jì)面外應(yīng)力。平面應(yīng)力模型的塑性變形偏大，因?yàn)樵撃Ｐ湍芸紤]中輥下方的應(yīng)力集中，材料更容易發(fā)生塑性屈服。

圖3 中輥工藝力與輸出曲率隨滾彎長(zhǎng)度的變化

圖4給出了達(dá)到穩(wěn)態(tài)后各模型在跨內(nèi)的塑性應(yīng)變分布?？梢钥闯?，平面應(yīng)力模型受應(yīng)力集中的影響，塑性應(yīng)變關(guān)于中性層不對(duì)稱，受壓側(cè)偏大。由于在計(jì)算中該模型的接觸誤差難以保持恒定，所得的結(jié)果始終有微小波動(dòng)，故應(yīng)變分布和輸出曲線不光滑(圖3和圖4(b))。此外該模型還偶爾會(huì)出現(xiàn)接觸穿透，如圖4的箭頭處，塑性應(yīng)變突然下降。板殼模型的塑性應(yīng)變分布雖然最為光滑，但數(shù)值明顯低于其他模型。而虛擬載荷法從原理上杜絕了偶然穿透現(xiàn)象，能給出光滑的結(jié)果曲線，并且數(shù)值上與曲率積分法給出的理論解最為接近。

圖4 穩(wěn)定狀態(tài)下變形區(qū)內(nèi)軸向塑性應(yīng)變的分布

達(dá)到穩(wěn)定狀態(tài)后變形區(qū)出口處的橫截面上的軸向殘余應(yīng)力分布如圖5所示?？梢钥吹教摂M載荷法與理論解基本吻合；平面應(yīng)力模型在工件頂部的殘余應(yīng)力明顯大于理論解，在工件底部則略偏?。欢鍤つＰ徒o出的殘余應(yīng)力則整體偏小。

圖5 穩(wěn)定狀態(tài)下殘余應(yīng)力在厚度方向上的分布

在上述算例的基礎(chǔ)上逐漸增加中輥下壓量，得到穩(wěn)態(tài)下的輥輪反力和曲率隨下壓量變化的關(guān)系如圖6所示。隨著變形區(qū)內(nèi)的工件構(gòu)型逐漸偏離小撓度情形，各有限元模型的誤差增加，預(yù)測(cè)的輸出曲率逐漸小于曲率積分得到的理論值。相比之下，平面應(yīng)力模型得到的輸出曲率誤差最小，而虛擬載荷法得到的下壓力與理論值吻合最好。

圖6 穩(wěn)定狀態(tài)下中輥工藝力和輸出曲率與下壓量的關(guān)系

5 結(jié) 論

針對(duì)三輥式滾彎的數(shù)值模擬，本文基于彈塑性歐拉梁模型提出了一種虛擬載荷法。相比通常的有限元模擬方案，其具有如下特點(diǎn)。

(1) 方法基于歐拉觀點(diǎn)，便于處理定義在空間坐標(biāo)上的力和位移載荷。

(2) 只對(duì)變形區(qū)內(nèi)的空間劃分網(wǎng)格，減少了網(wǎng)格數(shù)量，計(jì)算效率高。

(3) 不需要建立模具網(wǎng)格，能避免接觸面近似和偶然穿透導(dǎo)致的誤差。

本文方法的主要限制是只適于小曲率情形。對(duì)大曲率滾彎，適用的梁板理論、材料流動(dòng)追蹤和接觸處理都將變得更加復(fù)雜，仍需進(jìn)一步的研究。