周　研,雙遠(yuǎn)華,趙春江,茍毓俊,劉邱祖

(太原科技大學(xué) a.機械工程學(xué)院,b.材料科學(xué)與工程學(xué)院,c.太原重型機械裝備協(xié)同創(chuàng)新中心,太原 030024)

?

無縫鋼管縱連軋無網(wǎng)格法熱力耦合分析

周研a,c,雙遠(yuǎn)華b,c,趙春江b,c,茍毓俊b,c,劉邱祖a

(太原科技大學(xué) a.機械工程學(xué)院,b.材料科學(xué)與工程學(xué)院,c.太原重型機械裝備協(xié)同創(chuàng)新中心,太原 030024)

摘要:采用混合交換法修正無網(wǎng)格移動最小二乘近似函數(shù),以便直接施加邊界條件;采用伽遼金法構(gòu)建速度場剛度矩陣;采用配點法構(gòu)建溫度場剛度矩陣;采用間接耦合法將速度場與溫度場耦合求解;最終推導(dǎo)出了剛塑性無網(wǎng)格法熱力耦合計算公式,實現(xiàn)了縱連軋過程的熱力耦合模擬。通過數(shù)值模擬與實驗研究對比表明,模擬結(jié)果與實驗結(jié)果吻合度較高,充分驗證了本文給出的無網(wǎng)格法的可靠性和正確性。

關(guān)鍵詞:無縫鋼管;縱連軋;無網(wǎng)格法;熱力耦合分析

無網(wǎng)格法(Meshless Method)對求解域進(jìn)行節(jié)點離散與節(jié)點近似,避免求解過程對網(wǎng)格的依賴,無需對網(wǎng)格進(jìn)行劃分或重構(gòu),保證計算精度的同時減小了計算難度[1-2],在工程數(shù)值模擬應(yīng)用中越來越受到人們的重視。國內(nèi)外許多學(xué)者均在嘗試使用無網(wǎng)格方法解決金屬塑性成形問題。最早，由美國學(xué)者CHEN et al將RKPM方法應(yīng)用于金屬環(huán)件延伸、冷墩粗和壓縮過程問題的研究[3-4]；趙國群等利用剛塑性材料假設(shè),并對本質(zhì)邊界條件采用變換法處理,對典型墩粗過程進(jìn)行了模擬[5];李光耀等利用RKPM 法對三維彈塑性金屬材料的成形問題進(jìn)行了模擬[6]；孫杰等基于無網(wǎng)格徑向點插值法(RPIM)與伽遼金法(EFG)對斜軋延伸過程進(jìn)行了數(shù)值仿真[7]。

筆者采用混合交換法[8]修正無網(wǎng)格移動最小二乘近似函數(shù),以便直接施加邊界條件;采用伽遼金法構(gòu)建速度場剛度矩陣,該算法穩(wěn)定性好、分析精度高;采用配點法構(gòu)建溫度場剛度矩陣,以簡化溫度場剛度矩陣構(gòu)建過程,并采用隱式Euler向后差分格式進(jìn)行溫度場時間域離散。在熱力耦合分析中,采用間接耦合法[9]，將速度場與溫度場兩個相對獨立的子系統(tǒng)耦合求解;借助本構(gòu)關(guān)系,將變形與傳熱相互影響同時考慮,最終推導(dǎo)出了剛塑性無網(wǎng)格法熱力耦合分析公式,實現(xiàn)縱連軋過程的熱力耦合分析。通過數(shù)值模擬與實驗研究對比發(fā)現(xiàn),模擬結(jié)果與實驗結(jié)果吻合度較高,驗證了本文給出的無網(wǎng)格法的可靠性和正確性。

1最小二乘近似與混合交換修正

求解域Ω內(nèi)，一點x的速度u(x)的移動最小二乘近似(MLS)為

(1)

式中：uI為節(jié)點速度向量；u為廣義速度向量；N為離散節(jié)點數(shù)目。MLS形函數(shù)矩陣為

由于MLS形函數(shù)在邊界處不滿足Kronecker delta條件,不能直接施加本質(zhì)邊界條件,需要進(jìn)行修正,修正的方法有完全交換法[9]與局部交換法[8]。當(dāng)節(jié)點數(shù)目較多時,使用完全交換法的計算工作量也隨之增加。因此，筆者采用混合交換法對u進(jìn)行局部修正,將N個離散節(jié)點分為NE個非約束節(jié)點和NB個約束節(jié)點,N=NE+NB。式(1)中,令x=xJ,LIJ=ΦI(xJ)，得到

(2)

式中，I為單位矩陣。

(3)

式中：u*為經(jīng)過混合交換后的速度向量；IE為NE×NE單位矩陣;Λ*為混合變換矩陣,

由式(3)可得

(4)

將式(4)代入式(2)得到局部修正的速度場函數(shù)向量

溫度場向量也需采用同樣的變換。

2連軋過程剛黏塑性無網(wǎng)格法

2.1速度場剛度矩陣的建立

基于剛黏塑性材料條件假設(shè),應(yīng)用不完全廣義變分原理,采用罰函數(shù)引入體積不變條件,在局部坐標(biāo)系下施加摩擦力邊界條件,摩擦模型為反正切模型,真實解使系統(tǒng)能量速率泛函的一階變分為零,即真實解滿足如下剛度方程

(5)

式中：ΠD為應(yīng)變能速率項;ΠP為引入體積不變條件的罰函數(shù)項;Πf為摩擦功率項。通過式(5)可以得到剛黏塑性無網(wǎng)格伽遼金法的速度場剛度方程

(6)

式中，α為懲罰因子。該方程為非線性方程,采用Newton-Raphson 迭代法對(6)式進(jìn)行求解,直到獲得收斂解。式(6)經(jīng)過Newton-Raphson迭代，進(jìn)行線性化處理后,并在局部坐標(biāo)系下總裝剛度方程,得到如下矩陣形式

(7)

2.2無網(wǎng)格熱力耦合分析

假設(shè)材料導(dǎo)熱各向同性,則節(jié)點xs瞬態(tài)溫度應(yīng)當(dāng)滿足如下平衡微分方程:

(8)

并且在給定初始條件后滿足如下邊界條件:

(9)

(10)

式中:T為某節(jié)點溫度;h為工件與環(huán)境的總體換熱系數(shù);Ta為工件表面溫度。工件接觸表面由于摩擦產(chǎn)生的熱通量為:

(11)

(12)

式中:hT為熱交換系數(shù);TT為工模具表面溫度。將局部修正后的無網(wǎng)格形函數(shù)代入式(8)、(9)得到如下方程,

(13)

(14)

(15)

(16)

其中,lxs,mxs,nxs為接觸點xs的外法矢方向余弦值,式(13)-(16)寫成矩陣形式得

(17)

式中:

KT=[HT(x1),HT(x2),…,HT(xN),MT(x1),

MT(x2),…,MT(xNf),NT(x1),NT(x2),…,

NT(xNC),OT(x1),OT(x2),…,OT(xNC)];

CT=[CT(x1),CT(x2),…，CT(xN),0,…,0] ;

QT=[q(x1),q(x2),…,q(xN),hTa,hTa,…,

C(xs)=[-ρCΨ1(xs),-ρCΨ2(xs),…,

-ρCΨN(xs)] ;

H(xs)=[H1(xs),H2(xs),…,HN(xs)];

M(xs)=[M1(xs),M2(xs),…,MNf(xs)];

N(xs)=[N1(xs),N2(xs),…,NNc(xs)] ;

O(xs)=[O1(xs),O2(xs),…,ONc(xs)] ;

HI(xs)=kΨI,xx(xs)+kΨI,yy(xs)+kΨI,zz(xs) ;

MI(xs)=klxsΨI,x(xs)+kmxsΨI,y(xs)+

knxsΨI,z(xs)-hΨI(xs) ;

NI(xs)=klxsΨI,x(xs)+kmxsΨI,y(xs)+

knxsΨI,z(xs) ;

OI(xs)=klxsΨI,x(xs)+

kmxsΨI,y(xs)+knxsΨI,z(xs)-hTΨI(xs) .

(18)

取β=1,得到溫度場剛度矩陣后,采用文獻(xiàn)[10]提出的間接耦合迭代法即可得到穩(wěn)定的溫度場與速度場。具體實施步驟:先在每一個迭代步中首先計算節(jié)點溫度場;再由所得溫度場計算節(jié)點速度場;然后重復(fù)上兩步,直到節(jié)點速度場與溫度場同時收斂；最后更新節(jié)點坐標(biāo)進(jìn)行下一迭代步。

3鋼管縱連軋過程熱耦合分析

3.1縱連軋工藝及參數(shù)

管坯材料選用20#無縫鋼管,不同溫度不同等效應(yīng)變率下的流變應(yīng)力、應(yīng)變曲線見圖1所示。

a-2 s-1;b-20 s-1;c-200 s-1圖1　不同應(yīng)變率等效應(yīng)力應(yīng)變曲線Fig.1　Equivalent stress and strain curves under different strain rates

本次模擬與實驗中設(shè)置縱連軋機架為3架,同一機架內(nèi)軋輥互成120°角布置,并且相鄰機架前后軋輥成60°角。其中，包含3組共9個軋輥,限動芯棒被插入中空的毛管中,并參與整個塑性變形過程如圖2所示。

圖2　縱連軋工藝幾何模型Fig.2　Geometrical model of the longitudinal rolling process

3.2無網(wǎng)格分析參數(shù)

管坯取1/6模型進(jìn)行分析,離散節(jié)點總數(shù)為31 156個,徑向分5層排列,8個相鄰節(jié)點構(gòu)成一個背景積分網(wǎng)格;采用完全高斯積分,積分階次2×2×2,基函數(shù)采用線性基函數(shù),權(quán)函數(shù)選取三次樣條函數(shù)。節(jié)點影響域為球體域,影響域系數(shù)選取2.5;時間步長為0.01s。管坯某部分無網(wǎng)格離散點模型如圖3所示。

圖3　無網(wǎng)格法節(jié)點圖Fig.3　Node points in meshless method

3.3結(jié)果分析

圖4　節(jié)點速度與軋輥線速度差異Fig.4　Velocity differences between node points and roller

圖4為管坯表面與軋輥接觸點相對切向速度分布,各區(qū)域邊界為大致分界。由圖可知，速度矢量方向與軋制方向相同,表示該節(jié)點速度快于所處位置軋輥線速度,為前滑區(qū),反之為后滑區(qū)。圖中標(biāo)示出了前滑區(qū)、接觸區(qū)域大致的邊界。矢量的長度說明了速度差異的大小,可以看出,在軋輥出口截面的輥頂處前滑趨勢最大,并沿輥頂至輥縫,接觸點軋輥線速度逐漸增加,前滑趨勢減弱;在Z軸位置495mm附近,金屬外表面與軋輥相對速度基本相等。圖中還可以看出，金屬在孔型入口處的后滑趨勢最大,在出口處前滑趨勢最大。

圖5為2#機架孔型內(nèi)部變形區(qū)外表面溫度場分布。從圖中可以看出,由于輥底兩側(cè)首先和鋼管外表面接觸,這兩區(qū)域率先發(fā)生溫度下降;并由于金屬內(nèi)部熱量傳導(dǎo),金屬通過孔型后外表面溫度逐漸趨于相同。

圖5　2#機架孔型內(nèi)變形區(qū)表面溫度場分布Fig.5　Surface temperature field distributions of the deformation area in 2# pass

4實驗驗證

圖6-a為第三機架孔型出口斷面圖,圖6-b為該仿真同工藝軋卡實驗鋸切后的截面。通過比對，兩圖形狀基本相似。經(jīng)過測量,管坯周向外徑絕對誤差最大值0.148mm,相對誤差0.32%;周向壁厚絕對誤差最大值為0.03mm,相對誤差為0.49%。無網(wǎng)格法仿真的結(jié)果與實驗實測值基本吻合,如圖7所示。

圖6　產(chǎn)品形狀對比Fig.6　Product shape contrasts

5結(jié)論

本文采用無網(wǎng)格配點法構(gòu)建了溫度場剛度矩陣,使方程構(gòu)建過程與兩物理場耦合過程得以簡化。由于實驗檢測條件限制,沒有測得管坯軋制過程中的溫度變化數(shù)據(jù),但縱連軋工藝數(shù)值模擬與實驗測量所得管坯尺寸數(shù)據(jù)對比表明,本文所建立的無網(wǎng)格模型是正確的、有效的。

圖7　外徑(a)、壁厚(b)仿真與實驗測量對比Fig.7　Contrasts between numerical simulations and experimental measurements of diameter (a) and thickness (b)

參考文獻(xiàn)：

[1]張雄,劉巖.無網(wǎng)格法[M].北京:清華大學(xué)出版社,2005.

[2]LIUGuirong.Meshfreemethods:movingbeyondthefiniteelementmethod[M].US:CRCPRESS,2003.

[3]CHENJS,ROQUECMOL,PANChunhui,etal.Analysisofmetalformingprocessbasedonmeshlessmethod[J].JournalMaterialsProcessingTechchnology,1998,80-81:642-646.

[4]YOONSP,WUCT,WANGHuiping,etal.Efficientmeshfreeformulationformetalformingsimulation[J].JounalofEngineeringMaterialsandTechnology,2001,123(4):462-467.

[5]趙國群,王衛(wèi)東,欒貽國.金屬塑性成形過程無網(wǎng)格伽遼金法數(shù)值模擬技術(shù)研究[J].機械工程學(xué)報,2004,40(11):13-16.

[6]LIGuangyao,SIDIBEK,LIUGuirong.Meshfreemethodfor3Dbulkforminganalysiswithlowerorderintegrationscheme[J].EngineeringAnalysiswithBoundaryElements,2004,28(10):1283-1292.

[7]孫杰,胡建華,雙遠(yuǎn)華,等.斜軋延伸過程的無網(wǎng)格RPIM方法數(shù)值模擬[J].四川大學(xué)學(xué)報(工程科學(xué)版),2013,45(1):67-73.

[8]CHENJS,WANGHuiping.Newboundaryconditiontreatmentsinmeshfreecomputationofcontactproblems[J].ComputerMethodsinAppliedMechanicsandEngineering,2000,187(3):441-468.

[9]REBELON,KOBAYASHIS.Acoupledanalysisofviscoelasticdeformationandheattransfer[J].InternationalJournalofMechanicalScience,1980,22(11):707-718.

(編輯：龐富祥)

Meshless Thermo-mechanical Analysis of Continual Mandrel Rolling Process for Seamless Steel Tube

ZHOU Yana,c,SHUANG Yuanhuab,c,ZHAO Chunjiangb,c,GOU Yujunb,c,LIU Qiuzua

(TaiyuanUniversityofScienceandTechnology,a.CollegeofMechanicalEng.,b.CollegeofMaterialSci.andEng.,c.CollaborativeInnovationCenterofTaiyuanHeavyMachineryEquipment,Taiyuan030024,China)

Abstract:In this paper, the meshless moving least square approximated functions were corrected by using the mixed transformation method in order to enforce the essential boundary conditions. The stiffness matrices of the velocity field were the temperature field were constructed using the meshless Galerkin method and the meshless collocation method,respectively. The solutions of the two coupled fields were calculated with the indirect coupling method. Finally a formula was derived for the rigid-plastic meshless thermo-mechanically coupled analysis for the simulation of the continual mandrel rolling process. The comparison between the simulation and the experiments shows good accordance, validating the reliability and precision of our meshless method.

Key words:seamless steel tube;continual mandrel rolling process;meshless method;thermo-mechanical analysis

文章編號:1007-9432(2016)02-0160-05

*收稿日期：2015-08-20

基金項目：山西省科技攻關(guān)項目：鎂合金管材可控張力熱連軋工藝與設(shè)備開發(fā)(20140321008-08);國家自然科學(xué)基金資助項目：薄壁管材高速旋壓工藝擬動力學(xué)特性及可控機理研究(51375325)

作者簡介：周研(1983-),男,湖南寧鄉(xiāng)人，博士生，主要從事無縫鋼管軋制工藝及設(shè)備、軋制過程的數(shù)值模擬研究,(E-mail)zy_harry@vip.163.com

中圖分類號:TG335.71

文獻(xiàn)標(biāo)識碼:A

DOI：10.16355/j.cnki.issn1007-9432tyut.2016.02.007