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

鄭子君，王 洪

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

1 引 言

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

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

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

圖1 對稱式三輥滾彎工藝

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

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

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

2 虛擬載荷法

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

建立如圖1所示坐標(biāo)系。記t時刻某空間坐標(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)過微小時間步長Δt，位于x處的材料截面由原本位于x+s處的材料截面受輥輪帶動平移而來，即有Ct + Δt(x)=Ct(x+s)；又由步長Δt很小，可以假設(shè)任意材料點的應(yīng)力路徑是簡單的，則計算彎矩時，只需要知道該時間步內(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無關(guān)的常數(shù)，物理意義為該時間步內(nèi)的滾彎機的進料量。式(2)可進一步整理為

Δ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)第二項的意義是材料水平流動后，為維持t時刻構(gòu)型不變所需要的額外載荷，記為-ΔMa r t(x)。而式(3)第一項的物理意義是由變形引起的彎矩變化。若定義割線彎曲剛度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)的歐拉梁方程相比，差別僅是右端項除了實際彎矩增量ΔMe x(x)以外，還額外多了一項ΔMa rt(x)。它是空間構(gòu)型不變時，由于材料流動而導(dǎo)致的內(nèi)力不平衡量。由于該項在t+Δt時刻是已知的，故可以作為載荷項來處理。引入額外載荷項來考慮材料流動影響的方法可稱為虛擬載荷法。從定義容易看出，如果梁的截面參數(shù)C和塑性內(nèi)變量P在各處是均勻的，則虛擬載荷項自動消失。

3 有限元實施方案

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

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

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

(6)

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

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

(7)

(8)

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

(9)

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

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

(10)

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

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

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

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

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

4 算 例

取對稱式三輥滾彎機的輥距L=400 mm，輥徑Dt=Db=50 mm。中輥下壓量固定為d=3.8 mm。工件的橫截面為h=10 mm的矩形，初始時工件平直。材料采用理想彈塑性模型，楊氏模量E=200 GPa，屈服強度σs=200 MPa。總滾彎長度為600 mm。

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

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

求解得到的中輥的約束反力和出口處的曲率隨滾彎長度S的變化如圖3所示。可以看出各模型的結(jié)果曲線非常相似。大約在滾過1.5倍跨長后，各解基本達到穩(wěn)態(tài)。對比曲率積分法給出的穩(wěn)態(tài)解可知，虛擬載荷法與理論解最接近。板殼模型與其他模型相比，變形量明顯偏小，下壓力偏大，與文獻[11,12]結(jié)論相符，這是因為板殼模型中存在彎曲面外的正應(yīng)力，而其他模型均未計面外應(yīng)力。平面應(yīng)力模型的塑性變形偏大，因為該模型能考慮中輥下方的應(yīng)力集中，材料更容易發(fā)生塑性屈服。

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

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

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

達到穩(wěn)定狀態(tài)后變形區(qū)出口處的橫截面上的軸向殘余應(yīng)力分布如圖5所示。可以看到虛擬載荷法與理論解基本吻合；平面應(yīng)力模型在工件頂部的殘余應(yīng)力明顯大于理論解，在工件底部則略偏小；而板殼模型給出的殘余應(yīng)力則整體偏小。

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

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

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

5 結(jié) 論

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

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

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

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

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