橢圓形封頭旋壓成形過(guò)程數(shù)值模擬與分析

田文爽，孫衛(wèi)國(guó)

(1. 中國(guó)海洋石油集團(tuán)有限公司節(jié)能減排監(jiān)測(cè)中心，天津 300452； 2. 艾西斯騰德國(guó)有限公司，德國(guó) 慕尼黑 80331-81929)

由于化工裝置規(guī)模不斷擴(kuò)大，工藝設(shè)備逐漸呈現(xiàn)大型化的趨勢(shì)。傳統(tǒng)的封頭沖壓成形技術(shù)在大型封頭生產(chǎn)中遇到了很多困難。旋壓成形技術(shù)以其獨(dú)特的優(yōu)點(diǎn)，已成為大型封頭的重要生產(chǎn)手段，國(guó)內(nèi)外對(duì)封頭旋壓成形工藝的研究也不斷深入【1】。隨著科學(xué)技術(shù)的飛速發(fā)展，封頭旋壓成形工藝的數(shù)值模擬研究得到了迅速發(fā)展并提高到了一個(gè)新的階段。

林波等人【2】提出兩模法和兩輪法的旋壓成形試驗(yàn)方案，并設(shè)計(jì)不同旋壓試驗(yàn)方案的芯模和旋輪工裝，提出了兩輪法旋壓超半球殼體的多道次拉深旋壓成形工藝。劉興家等人【3】建立了成形過(guò)程中各道次旋輪運(yùn)動(dòng)軌跡控制參數(shù)間的遞推關(guān)系， 求得各道次旋輪軌跡的控制參數(shù)值， 較系統(tǒng)地提出了一種合理分配各道次旋壓變形量的方法。Shimizu【4】和 Kawai等人【5】對(duì)旋壓工藝中一些旋輪的運(yùn)動(dòng)問(wèn)題進(jìn)行了討論。這些研究尚未涉及模型計(jì)算時(shí)間長(zhǎng)，旋輪軌跡復(fù)雜和回彈嚴(yán)重的問(wèn)題。

張艷秋等人【6】對(duì)薄壁鋁合金封頭冷旋壓成形過(guò)程進(jìn)行了彈塑性有限元數(shù)值模擬，揭示了成形時(shí)的塑性變形流動(dòng)規(guī)律，并分析了主要工藝參數(shù)對(duì)成形質(zhì)量的影響規(guī)律及失效問(wèn)題產(chǎn)生的原因。胡福泰等人【7】用有限元數(shù)值計(jì)算的方法，合理確定了封頭旋壓計(jì)算的力學(xué)模型。張晉輝等人【8】建立了錐形件剪切旋壓的三維有限元模型，進(jìn)而獲得了偏離率、旋輪圓角半徑、旋輪進(jìn)給量、芯模轉(zhuǎn)速及旋輪直徑對(duì)LY12M錐形件剪切旋壓力和壁厚差的影響規(guī)律。李文平等人【9】建立了碟形封頭冷旋壓的三維有限元數(shù)學(xué)模型，分析了封頭在冷旋壓成形過(guò)程中應(yīng)力應(yīng)變的分布規(guī)律，計(jì)算了旋壓力并分析了其變化的原因。

運(yùn)用數(shù)值模擬，可以對(duì)金屬壓力加工過(guò)程中金屬的流動(dòng)進(jìn)行分析，從而得到應(yīng)力應(yīng)變分布規(guī)律等，進(jìn)而預(yù)測(cè)金屬變形情況以及對(duì)工藝參數(shù)進(jìn)行優(yōu)化。有限元分析軟件方便的操作以及超強(qiáng)的仿真性促進(jìn)了數(shù)值模擬技術(shù)在鍛造和旋壓等行業(yè)中的廣泛應(yīng)用【10】，在保證工件質(zhì)量、減少材料消耗、提高生產(chǎn)效率、縮短試制周期方面顯示出無(wú)可比擬的優(yōu)越性。本文基于ANSYS/LS-DYNA有限元分析軟件，建立了封頭旋壓成形的三維有限元模型，利用有限元模型，對(duì)封頭成形過(guò)程中的應(yīng)力場(chǎng)、應(yīng)變場(chǎng)的變化進(jìn)行分析，為封頭旋壓成形實(shí)際生產(chǎn)提供了可靠參考。

1 大變形成形理論基礎(chǔ)

1.1 大變形動(dòng)力學(xué)數(shù)值計(jì)算方法

在實(shí)際工程中，經(jīng)常會(huì)遇到幾何大變形問(wèn)題，例如板殼結(jié)構(gòu)的大撓度、屈曲和過(guò)屈曲問(wèn)題等，這時(shí)需要使用大變形條件下的應(yīng)力與應(yīng)變度量。

對(duì)于與變形歷史不相關(guān)的彈性大變形問(wèn)題可以采用全量方法研究【11】，也就是直接求已知載荷與約束作用下的結(jié)構(gòu)變形與應(yīng)變、應(yīng)力。將變形前的構(gòu)形作為參考構(gòu)形，變形后的構(gòu)形及相應(yīng)的應(yīng)變、應(yīng)力作為待求量，采用牛頓法、擬牛頓法等方法進(jìn)行迭代【12】。

對(duì)與變形歷史有關(guān)的大變形問(wèn)題，如材料粘彈性模型以及慣性等時(shí)間效應(yīng)，必須采用增量方法，即將時(shí)間變量離散成某個(gè)時(shí)間序列：t0=0,t1,t2,…,tn,tn+1，然后求這些離散時(shí)間點(diǎn)上的數(shù)值解。求解方法依參考構(gòu)形選擇的不同可以分為更新拉格朗日格式(U.L.)以及完全拉格朗日格式(T.L.)【13】。更新拉格朗日格式在計(jì)算[t,t+Δt]區(qū)間的所有變量時(shí)，以t時(shí)刻的構(gòu)形作為參考，應(yīng)力、應(yīng)變描述主要采用Euler應(yīng)力和關(guān)于現(xiàn)時(shí)構(gòu)形的無(wú)限小應(yīng)變；完全拉格朗日以t0=0時(shí)刻的構(gòu)形作為參考構(gòu)形，應(yīng)力、應(yīng)變描述主要采用Kirchhoff應(yīng)力和關(guān)于初始構(gòu)形定義的Green應(yīng)變【14】。

1.2 大變形動(dòng)力學(xué)有限元基本解法與求解過(guò)程

求解彈性動(dòng)力學(xué)振動(dòng)響應(yīng)主要有以下3類方法：時(shí)域方法、頻域方法和響應(yīng)譜方法。其中時(shí)域方法根據(jù)解法的不同又可以分為直接積分法、模態(tài)疊加法與狀態(tài)空間法。直接積分法又可以分為中心差分法、Houbolt法、Wilsonθ法以及Newmark法等。中心差分法是LS-DYNA所采用的主要算法，在求解時(shí)，假定t0=0,t1,t2,…,tn時(shí)刻的節(jié)點(diǎn)位移、速度與加速度均為已知，求解tn+1(t+Δt)時(shí)刻的結(jié)構(gòu)響應(yīng)【15】。

2 標(biāo)準(zhǔn)橢圓形封頭有限元模型的簡(jiǎn)化及參數(shù)確定

2.1 封頭有限元模型的建立

由于封頭旋壓成形過(guò)程比較復(fù)雜，影響因素較多，在數(shù)值模擬中，為了能夠建立描述其成形過(guò)程的有限元模型，有必要對(duì)實(shí)際加工過(guò)程做出合理的簡(jiǎn)化。在模型中，旋壓輥與成形輥設(shè)為剛性體，并簡(jiǎn)化為2個(gè)直徑相同的球殼。板坯預(yù)壓成形部分為剛體， 需要旋壓部分為變形體， 變形體部分簡(jiǎn)化為圓錐面， 與預(yù)壓件剛體部分末端相切。預(yù)壓件剛體部分繞中心軸線轉(zhuǎn)動(dòng)， 旋壓輥和成形輥在所在平面內(nèi)做平面運(yùn)動(dòng)。為了使模擬條件與真實(shí)情況一致， 模型的工藝參數(shù)均與實(shí)際工作參數(shù)相同?？紤]是冷成形過(guò)程， 忽略溫度變化的影響。

標(biāo)準(zhǔn)橢圓形封頭尺寸通過(guò)查閱JB/T 4729—1994確定，選擇直徑為φ2 000 mm，厚度16 mm，直邊部分40 mm的封頭作為研究對(duì)象。板坯材料選用為Q345R，材料性質(zhì)見(jiàn)表1。板坯預(yù)壓件采用shell163單元，預(yù)壓件的尺寸應(yīng)用周長(zhǎng)法，通過(guò)MATLAB程序確定，LS-DYNA中建模需要的點(diǎn)的坐標(biāo)也通過(guò)MATLAB程序計(jì)算得到。板坯預(yù)壓件繞中心軸線旋轉(zhuǎn)的速度取75 r/min。旋壓輥和成形輥的直徑均設(shè)為φ0.2 m。為了提高計(jì)算效率，旋壓輥和成形輥均采用shell163單元，厚度設(shè)為5 mm，摩擦系數(shù)取0.25。

表1 封頭材料參數(shù)

在模型中，所有單元均采用四面體劃分網(wǎng)格，為了提高計(jì)算精度，同時(shí)控制計(jì)算時(shí)間，旋壓輥和成形輥的單元尺寸設(shè)為0.1，板坯的單元尺寸設(shè)為0.05。最終建立的模型及網(wǎng)格劃分如圖1所示。

圖1 封頭旋壓成形有限元模型及網(wǎng)格劃分示意

2.2 旋壓成形旋輪軌跡的確定

橢圓封頭由橢圓部分和直邊部分組成，旋壓成形工藝中兩部分要分別旋制而成，建立封頭旋壓成形工藝有限元模型的關(guān)鍵在于給出精確的旋輪運(yùn)動(dòng)軌跡。在模型中，毛坯旋轉(zhuǎn)的同時(shí)，成形輥和旋壓輥按照給定的運(yùn)動(dòng)軌跡運(yùn)動(dòng)，使板坯變形，逐漸成為所要加工的形狀。

對(duì)圖1模型，為了能正確地模擬結(jié)構(gòu)的響應(yīng)，必須定義與指定時(shí)間間隔相對(duì)應(yīng)的載荷，所以在ANSYS/LS-DYNA中用一對(duì)數(shù)組參數(shù)定義載荷：一個(gè)定義時(shí)間，另一個(gè)定義載荷。

與時(shí)間數(shù)組相對(duì)應(yīng)的是板坯旋轉(zhuǎn)的角速度數(shù)組，應(yīng)用MATLAB程序計(jì)算確定旋壓輥和成形輥在x方向和y方向運(yùn)動(dòng)的位移數(shù)組。由于旋壓過(guò)程中，毛坯繞主軸轉(zhuǎn)動(dòng)，故旋輪只需做平面運(yùn)動(dòng)。圖2為旋輪在旋制封頭橢圓部分時(shí)，在運(yùn)動(dòng)平面上軌跡的示意圖，因?yàn)樾龎耗Ｐ褪菍?duì)稱結(jié)構(gòu)，所以取其中的一半進(jìn)行計(jì)算并輸出旋壓輥和成形輥位移數(shù)組數(shù)據(jù)。

圖2 旋輪軌跡計(jì)算示意

3 橢圓形封頭旋壓成形數(shù)值模擬結(jié)果的分析

3.1 旋壓成形數(shù)值模擬結(jié)果

將數(shù)組計(jì)算結(jié)果導(dǎo)入ANSYS/LS-DYNA，并將載荷分別施加到相應(yīng)的part上。程序的計(jì)算時(shí)間根據(jù)旋輪的轉(zhuǎn)速和旋輪進(jìn)給比確定，并體現(xiàn)在時(shí)間數(shù)組里，載荷加載完成之后對(duì)模型進(jìn)行求解并得到結(jié)果。

圖3顯示的是在經(jīng)過(guò)600 s過(guò)程結(jié)束時(shí)的成形示意圖，從中可以看出，通過(guò)數(shù)值計(jì)算的手段，實(shí)現(xiàn)了封頭從旋壓初始毛坯結(jié)構(gòu)到實(shí)際產(chǎn)品的成形，數(shù)值計(jì)算的最終幾何結(jié)構(gòu)接近于實(shí)際標(biāo)準(zhǔn)橢圓形封頭結(jié)構(gòu)，表明本文所開(kāi)發(fā)的三維橢圓形封頭旋壓成形過(guò)程的有限元數(shù)值模型具有一定的可行性和可靠性。另外，從圖3中可以看到封頭成形的直邊段出現(xiàn)褶皺，這是封頭成形過(guò)程存在的一種主要缺陷形式，需要理論研究者或者實(shí)際操作者慎重考慮，并盡力避免褶皺缺陷生成或限制其在合理的范圍內(nèi)。

圖3 旋壓成形數(shù)值模擬結(jié)果示意

3.2 封頭旋壓成形應(yīng)力、應(yīng)變場(chǎng)分布

封頭旋壓成形過(guò)程中4個(gè)不同時(shí)期的Mises等效應(yīng)力云圖如圖4所示。由圖4可以看出，等效應(yīng)力最大值總是出現(xiàn)在旋壓輥?zhàn)饔命c(diǎn)附近的區(qū)域內(nèi)，由于本模型選用理想彈塑性材料，最大等效應(yīng)力為345 MPa，因此一旦材料等效應(yīng)力達(dá)到該值，金屬應(yīng)力將不再增加，而主要體現(xiàn)在塑性變形。同時(shí)，隨著變形量的變大，等效應(yīng)力達(dá)到最大值的區(qū)域范圍越來(lái)越大，這是因?yàn)榇藭r(shí)發(fā)生了更為劇烈的塑性變形。

圖4 不同階段mises等效應(yīng)力云圖

封頭旋壓成形過(guò)程中不同時(shí)期的等效塑性應(yīng)變?nèi)鐖D5所示，不同于等效應(yīng)力云圖，等效應(yīng)變?cè)茍D反映的是結(jié)構(gòu)的變形行為。由圖5可以看出，在旋壓成形的初始階段，也就是成形過(guò)程的小變形階段，等效塑性應(yīng)變很小，包含的區(qū)域也很小，隨后的幾個(gè)階段，即成形過(guò)程的大變形階段和直邊成形階段，毛坯的塑性變形較大，等效塑性應(yīng)變逐漸增大，云圖包含的區(qū)域也逐漸增大。

3.3 封頭旋壓成形厚度減薄率變化分布

封頭厚度減薄率是封頭成形質(zhì)量控制的一個(gè)非常重要的指標(biāo)，封頭旋壓成形過(guò)程中不同階段毛坯厚度的分布云圖如圖6所示。

圖5 不同階段塑性等效應(yīng)變分布云圖

圖6 不同時(shí)刻封頭厚度減薄率分布

由圖6可以看出，在旋壓初始階段，毛坯厚度沒(méi)有顯著變化。隨著旋壓成形的進(jìn)行，厚度減薄率峰值逐漸增加，成形中期出現(xiàn)減薄率最大部位。如在225 s時(shí)，最大厚度減薄率為3.05%，而在445 s時(shí)，最大厚度減薄率為11.82%，表明在該段時(shí)間內(nèi)，塑性變形較為劇烈，最大厚度減薄率增加較快。在600 s時(shí)，最大厚度減薄率峰值11.86%，且位置與445 s 時(shí)相比不再進(jìn)一步變化，說(shuō)明劇烈塑性變形階段之后變形過(guò)程趨緩，成形過(guò)程厚度減薄率不再繼續(xù)增加。

整個(gè)成形過(guò)程，最大厚度減薄率11.86%，滿足GB/T 25198—2010附錄J中對(duì)DN 2000、名義厚度16 mm的橢圓形封頭厚度減薄率不超過(guò)13%的要求【16】。

3.4 封頭褶皺形成分析

封頭直邊段某節(jié)點(diǎn)速度的時(shí)間歷程曲線如圖7 所示。由圖7可以看出，在旋壓成形400 s以前，速度趨于穩(wěn)定，但超過(guò)400 s以后，該節(jié)點(diǎn)速度出現(xiàn)一定的波動(dòng)。在宏觀上表現(xiàn)為此時(shí)結(jié)構(gòu)端部，即直邊段發(fā)生了褶皺，從而引起速度的變化。可見(jiàn)在封頭成形后半部分，封頭毛坯的金屬流動(dòng)過(guò)程復(fù)雜，易產(chǎn)生褶皺。同時(shí)也可以看出，隨著旋壓過(guò)程的結(jié)束，褶皺現(xiàn)象明顯減少，并逐漸趨于平穩(wěn)。

圖7 封頭直邊段某節(jié)點(diǎn)速度的時(shí)間歷程曲線

4 結(jié)語(yǔ)

1) 根據(jù)旋壓機(jī)的工作原理，結(jié)合生產(chǎn)實(shí)際，利用幾何方法求出封頭旋壓成形工藝中旋輪運(yùn)動(dòng)軌跡的計(jì)算公式。采用MATLAB程序?qū)崿F(xiàn)了旋輪軌跡的計(jì)算，并實(shí)現(xiàn)了與ANSYS/LS-DYNA軟件的對(duì)接， 為應(yīng)用數(shù)值模擬的方法研究封頭旋壓成形工藝提供了有效的改進(jìn)方法和可靠的依據(jù)。

2) 基于ANSYS/LS-DYNA有限元分析軟件，應(yīng)用MATLAB程序得到的旋輪軌跡，建立了封頭旋壓成形的三維有限元模型，對(duì)模型進(jìn)行了合理簡(jiǎn)化。利用有限元模型，對(duì)封頭成形過(guò)程中的應(yīng)力場(chǎng)、應(yīng)變場(chǎng)進(jìn)行了分析，得到了封頭旋壓過(guò)程中不同階段的應(yīng)力場(chǎng)、應(yīng)變場(chǎng)分布規(guī)律以及封頭壁厚的變化規(guī)律。

3) 由數(shù)值分析結(jié)果可知，等效應(yīng)力最大值出現(xiàn)在旋壓輥?zhàn)饔命c(diǎn)附近的區(qū)域內(nèi)，并且隨著變形量的增大，等效應(yīng)力最大值所在區(qū)域的范圍也越來(lái)越大；在旋壓成形過(guò)程中，隨著變形量的增大，等效塑性應(yīng)變逐漸增大，變形區(qū)域也逐漸增大；在旋壓初始階段，毛坯厚度減薄率較低，及至發(fā)生劇烈塑性變形的階段，毛坯的厚度減薄率逐漸增大并達(dá)到峰值，本文設(shè)置的參數(shù)下，減薄率能夠控制在12%以內(nèi)，符合相關(guān)標(biāo)準(zhǔn)要求。