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

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

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

?

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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

式中，I為單位矩陣。

(3)

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

由式(3)可得

(4)

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

溫度場(chǎng)向量也需采用同樣的變換。

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

2.1速度場(chǎng)剛度矩陣的建立

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

(5)

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

(6)

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

(7)

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

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

(8)

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

(9)

(10)

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

(11)

(12)

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

(13)

(14)

(15)

(16)

其中,lxs,mxs,nxs為接觸點(diǎn)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,得到溫度場(chǎng)剛度矩陣后,采用文獻(xiàn)[10]提出的間接耦合迭代法即可得到穩(wěn)定的溫度場(chǎng)與速度場(chǎng)。具體實(shí)施步驟:先在每一個(gè)迭代步中首先計(jì)算節(jié)點(diǎn)溫度場(chǎng);再由所得溫度場(chǎng)計(jì)算節(jié)點(diǎn)速度場(chǎng);然后重復(fù)上兩步,直到節(jié)點(diǎn)速度場(chǎng)與溫度場(chǎng)同時(shí)收斂；最后更新節(jié)點(diǎn)坐標(biāo)進(jìn)行下一迭代步。

3鋼管縱連軋過(guò)程熱耦合分析

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

管坯材料選用20#無(wú)縫鋼管,不同溫度不同等效應(yīng)變率下的流變應(yīng)力、應(yīng)變曲線見(jiàn)圖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í)驗(yàn)中設(shè)置縱連軋機(jī)架為3架,同一機(jī)架內(nèi)軋輥互成120°角布置,并且相鄰機(jī)架前后軋輥成60°角。其中，包含3組共9個(gè)軋輥,限動(dòng)芯棒被插入中空的毛管中,并參與整個(gè)塑性變形過(guò)程如圖2所示。

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

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

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

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

3.3結(jié)果分析

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

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

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

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

4實(shí)驗(yàn)驗(yàn)證

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

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

5結(jié)論

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

圖7　外徑(a)、壁厚(b)仿真與實(shí)驗(yàn)測(cè)量對(duì)比Fig.7　Contrasts between numerical simulations and experimental measurements of diameter (a) and thickness (b)

參考文獻(xiàn)：

[1]張雄,劉巖.無(wú)網(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]趙國(guó)群,王衛(wèi)東,欒貽國(guó).金屬塑性成形過(guò)程無(wú)網(wǎng)格伽遼金法數(shù)值模擬技術(shù)研究[J].機(jī)械工程學(xué)報(bào),2004,40(11):13-16.

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

[7]孫杰,胡建華,雙遠(yuǎn)華,等.斜軋延伸過(guò)程的無(wú)網(wǎng)格RPIM方法數(shù)值模擬[J].四川大學(xué)學(xué)報(bào)(工程科學(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

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

*收稿日期：2015-08-20

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

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

中圖分類號(hào):TG335.71

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

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