移動線載荷板局部非線性響應(yīng)的網(wǎng)格重劃分?jǐn)?shù)值方法

李元泰，趙耀，袁華，王波

1華中科技大學(xué)船舶與海洋工程學(xué)院，湖北武漢430074

2武昌船舶重工集團(tuán)有限公司，湖北武漢430064

移動線載荷板局部非線性響應(yīng)的網(wǎng)格重劃分?jǐn)?shù)值方法

李元泰1，趙耀1，袁華1，王波2

1華中科技大學(xué)船舶與海洋工程學(xué)院，湖北武漢430074

2武昌船舶重工集團(tuán)有限公司，湖北武漢430064

在船舶建造中，常涉及板的焊接、線加熱、局部滾壓成型等加工方式。該加工方式中的載荷作用區(qū)域相對于板的尺度較窄且具有高度的非線性特征，同時載荷穩(wěn)定移動。把握加工參數(shù)與響應(yīng)的關(guān)系在加工過程中很重要。由于載荷移動導(dǎo)致板局部非線性響應(yīng)的特性不斷隨之變化，故其相應(yīng)的計算非常復(fù)雜。常規(guī)的數(shù)值方法是主要采用的手段，但由于該方法在載荷全局作用區(qū)域上需預(yù)先將網(wǎng)格細(xì)密劃分，導(dǎo)致計算效率低下。如果考慮到在該加工方式中載荷作用區(qū)域之外的大部分區(qū)域處于弱非線性或線彈性狀態(tài)，在數(shù)值模擬中將載荷當(dāng)前作用區(qū)域的單元劃分細(xì)密，遠(yuǎn)離載荷作用區(qū)域劃分得相對較粗疏，同時載荷移動時細(xì)密網(wǎng)格隨著載荷一起移動，即采用網(wǎng)格重劃分?jǐn)?shù)值方法來分析上述問題就能節(jié)省大量計算時間。針對該方法，重點探究細(xì)密網(wǎng)格尺寸以及網(wǎng)格重劃分頻度對數(shù)值計算的影響及其與計算效率和精度的關(guān)系，并以板局部滾壓成型加工為例進(jìn)行數(shù)值計算。結(jié)果表明，該方法可作為移動載荷作用下板局部非線性響應(yīng)的一種高效的計算手段。

移動線載荷；網(wǎng)格重劃分頻度；細(xì)密網(wǎng)格尺寸

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/42.1755.TJ.20160921.1339.022.html期刊網(wǎng)址：www.ship-research.com

引用格式：李元泰，趙耀，袁華，等.移動線載荷板局部非線性響應(yīng)的網(wǎng)格重劃分?jǐn)?shù)值方法［J］.中國艦船研究，2016，11（5）：63-70.

LI Yuantai，ZHAO Yao，YUAN Hua，et al.Rezoning method in nonlinear simulation of plate under moving load［J］. Chinese Journal of Ship Research，2016，11（5）：63-70.

0　引言

在船舶金屬板結(jié)構(gòu)的加工成型過程中，常伴隨著載荷作用區(qū)域相對于被加工板面尺度較窄的移動線載荷作用（如熱載荷、力載荷等），如金屬板的焊接和線加熱成型；通過一對上下的凹凸?jié)L輪加壓滾動，對放置于凹凸?jié)L輪之間的板曲面成型的過程。這類加工的特點是，在載荷作用區(qū)域內(nèi)，不論是熱載荷還是力載荷，其載荷量均非常大，且作用過程穩(wěn)定移動。從而在載荷作用的狹窄區(qū)域附近，被加工板殼產(chǎn)生較大的塑性變形，局部響應(yīng)呈強非線性特征，而在遠(yuǎn)離這一區(qū)域的大部分被加工板殼則處于弱非線性或線彈性狀態(tài)，同時隨著載荷的移動，被加工板經(jīng)歷著加載、卸載的變化。

實際情況中，由于載荷大小及路線的變化，加工精度和效率與加工參數(shù)之間存在著復(fù)雜的對應(yīng)關(guān)系，因此高效率和高精度地計算預(yù)報不同加工參數(shù)與加工要求如變形量之間的關(guān)系是非常必要的。對于這樣一個涉及彈塑性同時包括載荷移動甚至接觸變化的復(fù)雜問題，使用完全的理論求解分析幾乎是不可能的，目前一般會借助于非線性有限元方法來解決。

常規(guī)的彈塑性數(shù)值方法對于處理這樣一類問題通常依據(jù)經(jīng)驗在非線性特征明顯的載荷區(qū)域周圍劃分十分細(xì)密的網(wǎng)格，以保證計算結(jié)果的準(zhǔn)確性。這樣導(dǎo)致在移動載荷尚未作用前和已經(jīng)作用后的整個載荷移動線路上均需要預(yù)先劃分為細(xì)密單元，使得網(wǎng)格數(shù)量過于龐大，導(dǎo)致非線性迭代計算十分耗時，往往模擬一次簡單的載荷計算過程需要很長時間，計算效率非常低下。因此，針對該問題的高效彈塑性計算數(shù)值方法研究就顯得尤為重要。如果能夠僅僅在載荷作用區(qū)域范圍內(nèi)加密計算網(wǎng)格尺度，而在其他地方適當(dāng)粗化計算網(wǎng)格尺度，減少節(jié)點數(shù)，縮小剛度矩陣規(guī)模，計算效率就會提高。網(wǎng)格重劃分?jǐn)?shù)值方法可以有效解決該類問題。其主要做法是根據(jù)分析過程中載荷當(dāng)前的位置和移動路線，在載荷作用的一段區(qū)域內(nèi)將網(wǎng)格進(jìn)行細(xì)密劃分，而對其他區(qū)域則適當(dāng)將網(wǎng)格粗疏劃分；而且在細(xì)密網(wǎng)格隨著載荷一起移動時，對載荷已經(jīng)離開的區(qū)域，細(xì)密網(wǎng)格重新劃分為粗疏網(wǎng)格，而對載荷即將進(jìn)入的區(qū)域，原粗疏網(wǎng)格劃分為細(xì)密網(wǎng)格。在網(wǎng)格重劃分前后，新舊網(wǎng)格模型的節(jié)點占據(jù)相同的幾何空間，其他場變量，如應(yīng)力應(yīng)變實現(xiàn)相應(yīng)的傳遞，以使分析求解得以連續(xù)進(jìn)行下去。國內(nèi)外學(xué)者對于該方法已發(fā)表了相應(yīng)的研究成果。Brown等［1］結(jié)合網(wǎng)格重劃分和動態(tài)子結(jié)構(gòu)技術(shù)，采用殼體單元分析了平板激光焊接后的變形和應(yīng)力，所使用的時間與常規(guī)的熱彈塑性計算相比，僅為后者的1/7，在保證精度的前提下大幅提高了效率，其分析局限于熱載荷作用。Huang等［2］以熱彈塑性有限元法為基礎(chǔ)開發(fā)了網(wǎng)格重劃分技術(shù)來預(yù)測線加熱過程中板的變形，并應(yīng)用于三維問題的分析，但載荷形式單一，也未涉及網(wǎng)格劃分的分析。Alsamhan等［3］利用實時網(wǎng)格重劃分技術(shù)，數(shù)值模擬了薄板局部滾壓成型過程，并通過實驗對比了薄膜應(yīng)力分布，二者吻合較好，成功提升了預(yù)測效率，雖然在接觸非線性下采用實時網(wǎng)格進(jìn)行了處理，但對細(xì)密網(wǎng)格劃分與變形結(jié)果尚未開展討論。另外，關(guān)于網(wǎng)格重劃分方法中的網(wǎng)格重劃分頻度，即在載荷作用路線上，一次計算的細(xì)密網(wǎng)格區(qū)域預(yù)先需要選擇劃分的次數(shù)及其在數(shù)值計算中的影響，在目前的研究中相對討論較少。

綜上所述，且考慮到關(guān)于熱載荷移動的研究結(jié)果相對較多，本文將重點考慮力載荷接觸變化作用條件下的網(wǎng)格重劃分?jǐn)?shù)值方法，并主要研究采用不同細(xì)密網(wǎng)格尺寸和不同頻度的網(wǎng)格劃分對計算結(jié)果的影響。這里，細(xì)密網(wǎng)格尺寸均指載荷作用區(qū)域的細(xì)密網(wǎng)格區(qū)域中網(wǎng)格大小；而網(wǎng)格重劃分頻度則指這一細(xì)密網(wǎng)格預(yù)先劃分的數(shù)量。在此討論的基礎(chǔ)上，本文選擇船舶金屬板材滾壓加工成型進(jìn)行計算分析，研究接觸非線性計算、計算效率及精度與細(xì)密網(wǎng)格尺寸和頻度的關(guān)系。通過滾壓加工成型實例，驗證細(xì)密網(wǎng)格尺寸和網(wǎng)格重劃分頻度對結(jié)果的影響。

1　網(wǎng)格重劃分?jǐn)?shù)值方法

網(wǎng)格重劃分?jǐn)?shù)值方法的基本思想如前所述。假想存在一塊區(qū)域，它與載荷速度一致，當(dāng)網(wǎng)格單元落入該封閉區(qū)域時將實現(xiàn)細(xì)密的網(wǎng)格劃分，而遠(yuǎn)離封閉區(qū)域則劃分得相對粗疏（圖1），從而能夠達(dá)到縮小整體剛度矩陣規(guī)模的效果。那么在兩次網(wǎng)格劃分之間，上一步的子模型形狀和相關(guān)的場變量等數(shù)據(jù)需要傳遞到下一步計算的子模型中?？紤]到分析的連續(xù)性與精確性，變化的網(wǎng)格形狀必須要準(zhǔn)確模擬，即保證前后網(wǎng)格節(jié)點占據(jù)相同的幾何空間，而相關(guān)的場變量要順利地實現(xiàn)映射，這樣才能保證計算結(jié)果的精度。具體示意圖如圖1所示。

網(wǎng)格疏密位置的變化、節(jié)點單元編號的改變導(dǎo)致有限元網(wǎng)格模型需不斷重新構(gòu)建。為了匹配下一步網(wǎng)格模型節(jié)點在上一步網(wǎng)格模型中的位置，判斷節(jié)點具體處于哪一個單元內(nèi)部成為問題的關(guān)鍵。不論是二維模型，還是三維模型均可采用Murti的判斷方法［4］。以八節(jié)點三維六面體單元為例，如果點在單元內(nèi)部則滿足

式中：i為單元面數(shù)；ai，bi，ci如圖2所示。如果點不在單元內(nèi)部，則Vi為負(fù)向量，Vi為0則表示點在單元邊界上。

圖1　網(wǎng)格重劃分和數(shù)據(jù)傳遞示意圖Fig.1Rezoning and data transfer

圖2　空間點的位置Fig.2A spatial point either inside or outside an element

根據(jù)等參元的性質(zhì)和六面體單元形函數(shù)Ni，單元中任意一點處的整體坐標(biāo)滿足

式中：x'y'z為整體坐標(biāo)系的坐標(biāo)值；ξ，η，ζ為局部坐標(biāo)系的坐標(biāo)值。已知上一步網(wǎng)格模型單元八節(jié)點坐標(biāo)值和位于該單元中一點的坐標(biāo)，采用Newton-Raphson方法求解得到局部坐標(biāo)系中的坐標(biāo)［5］。

式中，ai，bi，ci為各節(jié)點坐標(biāo)的函數(shù)。最后利用形函數(shù)對節(jié)點位移進(jìn)行插值，從而得到網(wǎng)格模型的節(jié)點在整體坐標(biāo)系中的坐標(biāo)值。相關(guān)場變量的傳遞則采用有限元軟件的變量映射功能，映射算法中采用插值技術(shù)將上一步網(wǎng)格中的殘余應(yīng)力和塑性應(yīng)變等變量外插到每一個單元節(jié)點上，然后對所有相鄰?fù)还?jié)點的單元變量取平均；新網(wǎng)格的積分點位置由舊網(wǎng)格的積分點來確定，先找到該點所屬舊網(wǎng)格的單元，這樣可以確定該點在單元中的位置，變量由舊網(wǎng)格的節(jié)點內(nèi)插到新網(wǎng)格的積分點上［6］。這樣就可以保證分析的連續(xù)性和準(zhǔn)確性。

對于力載荷接觸變化作用在板上時，由于非線性響應(yīng)特征明顯，一般需要將網(wǎng)格細(xì)密劃分才能保證計算的精確性；當(dāng)力載荷作用深度加大且接觸區(qū)域面積變大時，細(xì)密網(wǎng)格面積也相應(yīng)地要增加且密度也要加大，才能有利于接觸點對的建立，使計算分析成功。因為在接觸問題的有限元計算中，通常將接觸界面上的面稱為接觸塊，被動接觸塊上的節(jié)點和主動接觸塊與其接觸的點構(gòu)成一個接觸對。那么接觸點對的搜尋與合理建立是保證有限元分析結(jié)果可靠的關(guān)鍵。如果施載體，即主接觸體是剛體，那么作為變形體板殼的受載體，即從接觸體的網(wǎng)格必須細(xì)密劃分，以適應(yīng)剛體施載過程中接觸形狀的不斷變化。倘若網(wǎng)格比較粗疏，對于較大的載荷量，施載體進(jìn)入受載體的深度以及接觸面積就會相應(yīng)加大。主接觸面單元就會侵入從接觸面，導(dǎo)致計算失準(zhǔn)，甚至無法繼續(xù)分析下去。由于施載條件的變化，主、從接觸面上的變化如圖3所示，圖3（a）所示為由于從接觸面的網(wǎng)格比較粗糙，主接觸面上的單元侵入了從接觸面；圖3（b）所示為從接觸面的網(wǎng)格細(xì)化后，就防止了主接觸面的侵入［7］。

圖3　主、從接觸面上的變化Fig.3The change of master and slave contact surface

從另一方面講，即使選擇了保證不同載荷深度和大小都適當(dāng)?shù)木W(wǎng)格密度，也還存在一個預(yù)先細(xì)密網(wǎng)格區(qū)域劃定多大的問題。

由于密網(wǎng)格區(qū)域長度決定了需要多少次對網(wǎng)格模型重劃分才能完成一次數(shù)值模擬，所以頻度大小會影響數(shù)值模擬計算效率，頻度越小，實際上就越接近于完整模型。另外，頻度過大，則由于傳遞次數(shù)過多將導(dǎo)致精度損失，影響計算的精確性；反之則減少精度損失。

通過以上分析，對于板殼而言：

1）載荷增加，接觸面積增大，那么載荷區(qū)域的網(wǎng)格就需要細(xì)化，倘若細(xì)密網(wǎng)格尺寸過大，會造成計算失效，但細(xì)密網(wǎng)格尺寸太小又會導(dǎo)致計算時間冗長，效率低下；

2）載荷增加，接觸面積增大，每一步網(wǎng)格重劃分模型的加載區(qū)域長度就要相應(yīng)增加。如果劃分頻度過小，會造成計算效率降低，而劃分頻度過大又會造成精度損失。

2　計算結(jié)果與分析

本文選擇了船舶金屬板成型加工方式中常用的局部滾壓成型進(jìn)行計算。具體滾壓成型過程是將板材于一對上下凹凸的滾輪中通過，通過給上滾輪施壓并隨著下滾輪的回轉(zhuǎn)，依靠摩擦將板材向前送進(jìn)，從而形成變形。在該成型加工過程中涉及接觸非線性和滾輪持續(xù)穩(wěn)定的移動載荷作用，該工程對象從受力特征上是典型的在移動線載荷作用下非線性響應(yīng)的板。選擇此工程對象進(jìn)行網(wǎng)格重劃分?jǐn)?shù)值方法應(yīng)用計算分析具有代表性。矩形鋼板模型的尺寸為2 000 mm×1 000 mm×20 mm，材料性質(zhì)設(shè)為理想彈塑性模型，彈性模量為2.09× 105MPa，泊松比為0.3，屈服應(yīng)力σY=400MPa。

局部滾壓加載裝置主要包括上滾輪和下滾輪。其中上滾輪為下壓輪，尺寸最大半徑120 mm，輪轂厚度為200 mm；下滾輪尺寸最大半徑160 mm，輪轂厚度為240 mm；上滾輪和下滾輪的具體形狀如圖4所示。

圖4　上滾輪截面和下滾輪截面形狀Fig.4The section shape of upper roller and lower roller

滾壓過程的有限元模型如圖5所示，上滾輪和下滾輪簡化為剛體，滾壓線位于板寬中間位置。加載方式選擇施加位移，即給上滾輪一定的下壓量進(jìn)行加載。

圖5　常規(guī)有限元計算滾壓模型Fig.5Conventional FEM roll model

數(shù)值模擬步驟為：從板的一端接觸→下壓→連續(xù)勻速滾壓→板的另一端卸載。

即接觸過程為上滾輪下壓建立接觸，然后給定下壓量，滾壓過程是通過下滾輪順時針旋轉(zhuǎn)帶動板材來實現(xiàn)，同時為方便計算選擇讓板材勻速移動。整個滾壓過程邊界條件為：開始前平板采用遠(yuǎn)端位移約束，滾壓過程中解除位移約束，滾壓過程利用平板與滾輪之間的摩擦約束。

下面，分別對細(xì)密網(wǎng)格尺寸與網(wǎng)格重劃分頻度的影響進(jìn)行計算。其主要研究方案是對細(xì)密網(wǎng)格尺寸進(jìn)行計算，通過不同的下壓量和細(xì)密網(wǎng)格尺寸大小的組合進(jìn)行相應(yīng)計算，而對其中某一種下壓量的計算結(jié)果就計算效率和計算精度進(jìn)行綜合分析；對于網(wǎng)格重劃分頻度的計算，選擇某一細(xì)密網(wǎng)格尺寸，通過不同下壓量和頻度大小的組合進(jìn)行相應(yīng)計算，而對其中一種下壓量的完整計算結(jié)果就計算效率和計算精度進(jìn)行綜合分析。計算還對常規(guī)有限元方法進(jìn)行了建模計算，以此作為網(wǎng)格重劃分計算的比較分析依據(jù)。常規(guī)有限元計算模型如圖5所示。

2.1細(xì)密網(wǎng)格尺寸的影響

根據(jù)滾輪尺寸選重劃分頻度為8，而選擇4種不同的細(xì)密網(wǎng)格尺寸。據(jù)前文定義，頻度8表示隨著載荷的移動，8次重劃分網(wǎng)格模型并分別依次通過載荷作用進(jìn)行計算。根據(jù)前期計算調(diào)查，對于本次計算模型，細(xì)密網(wǎng)格區(qū)域的寬度均選為0.3 m較為合適，與完整模型的加載區(qū)域?qū)挾纫恢?，過渡區(qū)域和遠(yuǎn)離加載區(qū)域的網(wǎng)格也與完整模型保持一致。而沿載荷作用路線方向上的細(xì)密網(wǎng)格區(qū)域的長度在重劃分頻度為8的條件下取0.4 m。在上述區(qū)域內(nèi)，分別用4種細(xì)密網(wǎng)格進(jìn)行劃分，網(wǎng)格尺寸如圖6所示，網(wǎng)格尺寸依次逐漸變大。

Ⅰ：5 mm×10 mm×5 mm

Ⅱ：10 mm×10 mm×5 mm

Ⅲ：12 mm×10 mm×5 mm

Ⅳ：16 mm×10 mm×5 mm

圖6　4種細(xì)密網(wǎng)格尺寸示意圖Fig.6Four fine mesh size

對4種加載區(qū)域細(xì)密網(wǎng)格尺寸分別采取不同下壓量計算，計算結(jié)果如表1所示。

表1　不同下壓量下4種網(wǎng)格模型的計算結(jié)果Tab.1The results of four models under different amount of depression

從計算結(jié)果可知，當(dāng)細(xì)密網(wǎng)格尺寸為Ⅰ時，下壓量在1～15 mm都能完成計算；當(dāng)細(xì)密網(wǎng)格尺寸增大至Ⅱ時，下壓量在2～15 mm能完成計算；當(dāng)細(xì)密網(wǎng)格尺寸增大至Ⅲ時，下壓量在10～15 mm能完成計算；當(dāng)細(xì)密網(wǎng)格尺寸增大至Ⅳ時，下壓量在1～15 mm內(nèi)均不能完成計算。上述結(jié)果驗證了當(dāng)載荷區(qū)域細(xì)密網(wǎng)格尺寸過大時，會造成接觸點對的建立無法實現(xiàn)而造成計算失效。

2.1.1計算結(jié)果精度

計算結(jié)果精度是選取細(xì)密網(wǎng)格區(qū)域相同的3種不同細(xì)密網(wǎng)格尺寸Ⅰ，Ⅱ，Ⅲ與常規(guī)有限元模型在下壓量為10 mm的計算結(jié)果比較來進(jìn)行的。這里選擇滾壓0.2 m后板長方向0.1 m處橫向中點縱向截面的變形來進(jìn)行比較，結(jié)果如圖7所示。

圖7　不同細(xì)密網(wǎng)格尺寸模型的變形結(jié)果Fig.7Final deformation of different fine mesh size models

通過對比發(fā)現(xiàn)，3種細(xì)密網(wǎng)格尺寸的計算結(jié)果略有差別。在細(xì)密網(wǎng)格區(qū)域相同的情況下，細(xì)密網(wǎng)格尺寸越小，變形越接近于常規(guī)有限元模型計算變形結(jié)果。

2.1.2計算效率

當(dāng)加載區(qū)域細(xì)密網(wǎng)格尺寸為Ⅰ時，計算時間平均為60 min；當(dāng)加載區(qū)域細(xì)密網(wǎng)格尺寸為Ⅱ時，計算時間平均為26 min；當(dāng)加載區(qū)域細(xì)密網(wǎng)格尺寸板長方向為Ⅲ時，計算時間約為20 min。選取下壓量為5和10 mm，分別就3種網(wǎng)格尺寸Ⅰ，Ⅱ，Ⅲ的其中一步重劃分模型計算時間進(jìn)行比較，結(jié)果如圖8所示。

圖8　不同細(xì)密網(wǎng)格尺寸模型的計算時間Fig.8Computational time of different fine mesh size models

由圖可知，在下壓量為5 mm時，Ⅱ的計算時間相比Ⅰ縮短了50%，而Ⅲ的計算由于網(wǎng)格尺寸過大而失效；在下壓量為10 mm時，Ⅱ的計算時間相比Ⅰ縮短了50%，而Ⅲ的計算時間相比Ⅱ縮短了不到30%。這表明適當(dāng)增大細(xì)密網(wǎng)格尺寸可提高計算效率。如果考慮加工要求，同時又保證數(shù)值計算相對高效，細(xì)密網(wǎng)格尺寸選擇Ⅱ比較合適。

綜上所述，加載區(qū)域如何選擇合適的細(xì)密網(wǎng)格尺寸對計算精度和效率都有著重要的影響。細(xì)密網(wǎng)格尺寸太小會導(dǎo)致計算時間冗長，效率低下，但精度相對較高；隨著細(xì)密網(wǎng)格尺寸逐漸增大，會提高計算效率，但會損失精度；如果網(wǎng)格尺寸過大，會導(dǎo)致計算無法進(jìn)行。

2.2網(wǎng)格重劃分頻度的影響

本次模擬過程選擇細(xì)密網(wǎng)格尺寸為10 mm× 10 mm×5 mm的3種不同重劃分頻度進(jìn)行計算。如前所述，細(xì)密網(wǎng)格區(qū)域?qū)挾染赃x定為0.3 m，而沿載荷作用路線方向上的細(xì)密網(wǎng)格區(qū)域長度根據(jù)頻度的不同分別取為2，0.8，0.6，0.4 m。有限元網(wǎng)格模型如圖9所示，由圖9可知模型a，b，c，d分別為頻度0，4，6，8的4種重劃分模型，其周圍網(wǎng)格沿橫向和縱向2個方向由密逐步變疏，其中頻度為0的重劃分模型即為常規(guī)有限元計算模型。網(wǎng)格模型參數(shù)如表2所示，其中每種頻度的重劃分模型節(jié)點數(shù)是相同的。

圖9　不同頻度的有限元網(wǎng)格模型Fig.9FEM models of different frequencies

表2　4種不同頻度的網(wǎng)格模型參數(shù)Tab.2Four models'parameters of different frequencies

2.2.1計算結(jié)果的精度

計算結(jié)果的精度是通過與頻度0完整模型，即常規(guī)有限元模型在下壓量為5 mm的計算結(jié)果進(jìn)行比較。這里選擇板長邊方向0.1 m處和中點處的2個橫向截面變形以及板橫向中點縱向截面的變形來進(jìn)行比較，如圖10所示。

圖10　下壓量為5 mm時的板材最終變形圖Fig.10Final deformation of the plate under depression of 5 mm

由于在板滾壓過程中空間位置會變化，為便于比對，將重劃分模型的最終變形位置按照完整模型的坐標(biāo)系進(jìn)行相應(yīng)調(diào)整。先選取距離加載初始位置0.1和1 m處橫截面上的變形結(jié)果進(jìn)行對比，再選取中心線方向縱截面上的變形結(jié)果進(jìn)行對比。計算結(jié)果如圖11和圖12所示。

圖11　板長邊方向橫截面上的變形結(jié)果比較Fig.11Comparisons of final deformation in the cross section of the plate

圖12　板橫向中點縱截面上的變形結(jié)果比較Fig.12Comparisons of final deformation in the longitudinal section of the plate

由于滾壓線為矩形板中心線，所以滾輪沿著滾壓線加工矩形板，其變形的形狀是對稱的。由橫向結(jié)果對比可以看出：不論是在距離初始位置0.1 m處還是1 m處，完整模型與重劃分模型的變形趨勢一致，即角變形大小吻合很好，但模型b的變形結(jié)果相比于模型c，d來說更接近于完整模型a，不過計算時間要大于模型c，d。兩者之間的協(xié)調(diào)取決于個人需求。中心線方向縱向截面的變形趨勢可以說明重劃分模型與完整模型的板材最終變形形狀是一致的。從變形場的計算結(jié)果來看，頻度越小，計算效率越低，但精度較高；反之，頻度越高，計算效率相對越高，但精度會變差。

平板經(jīng)過滾壓載荷作用，由于存在卸載過程，載荷區(qū)域一定會存在較大的殘余應(yīng)力。以下壓量5 mm的計算結(jié)果為例，完整模型和頻度為8的網(wǎng)格重劃分模型d的最終應(yīng)力場計算結(jié)果的對比如圖13所示。

從圖中可以看出，網(wǎng)格重劃分?jǐn)?shù)值方法仍然再現(xiàn)了殘余應(yīng)力結(jié)果，殘余應(yīng)力分布與常規(guī)有限元的計算結(jié)果相近，不過相對而言網(wǎng)格重劃分?jǐn)?shù)值方法所得塑性應(yīng)力略微偏大，這是由于網(wǎng)格變化過程中力學(xué)信息的傳遞存在一定的計算誤差。

圖13　最終殘余應(yīng)力對比Fig.13Comparisons of final residual stress of the plate

2.2.2計算效率

完整模型與網(wǎng)格重劃分模型的計算時間如圖14所示。將3種不同下壓量下完成該次滾壓模擬的計算時間進(jìn)行對比。

圖14　不同頻度下模型的計算時間Fig.14Computational time of different frequencies'models

從圖中可以看出，使用網(wǎng)格重劃分?jǐn)?shù)值方法后，相對于完整模型，重劃分模型b的計算時間減少了約57%；網(wǎng)格重劃分模型c減少了約66%；重劃分模型d減少了約70%。其原因是網(wǎng)格重劃分模型中的節(jié)點數(shù)明顯少于完整模型的節(jié)點數(shù)，也即模型中總的自由度數(shù)前者少于后者，那么總體剛度矩陣規(guī)模也減少，內(nèi)存的需求也必將因網(wǎng)格重劃分?jǐn)?shù)值方法而得以減少，計算效率因此大幅提高。上述計算結(jié)果表明，網(wǎng)格重劃分?jǐn)?shù)值方法可以很好地解決效率問題，頻度越小，計算效率越低；反之頻度越高，計算效率相對越高。

3　結(jié)論

本文對網(wǎng)格重劃分?jǐn)?shù)值方法中細(xì)密網(wǎng)格尺寸與網(wǎng)格重劃分頻度對移動線載荷作用于板的數(shù)值模擬效率和精度的影響進(jìn)行了研究。以板滾壓加工成型數(shù)值模擬為例，討論了4種細(xì)密網(wǎng)格尺寸在不同下壓量下對效率和精度的影響，并采用不同頻度的網(wǎng)格重劃分?jǐn)?shù)值方法進(jìn)行了對比。主要結(jié)論如下：

1）由于接觸深度和接觸區(qū)域受載荷作用的影響，要保證獲得計算結(jié)果，載荷量增加，則網(wǎng)格劃分要細(xì)密且細(xì)密網(wǎng)格區(qū)域要增大；當(dāng)細(xì)密網(wǎng)格尺寸變小，計算精度變高但計算效率變低時，細(xì)密網(wǎng)格劃分的尺寸大小要根據(jù)實際加載情況綜合考慮。

2）在保證精度的前提下加大網(wǎng)格重劃分頻度會大幅提高效率，重劃分頻度越高，會在一定程度上提升效率，但會損失精度，所以頻度的大小也要根據(jù)綜合要求來選擇。

3）在適當(dāng)?shù)木W(wǎng)格劃分條件下，網(wǎng)格重劃分?jǐn)?shù)值方法總的來說可以實現(xiàn)，其在保證精度的前提下相比常規(guī)數(shù)值計算方法在計算效率上有大幅度提高，可以作為在移動線載荷作用下局部非線性響應(yīng)問題的有效計算手段。

［1］BROWN S B，SONG H.Rezoning and dynamicsubstructuringtechniques in FEM simulations of welding processes［J］.Journal of Engineering for Industry，1993，115（4）：415-423.

［2］HUANG H，SERIZAWA H，WANG J C，et al.Development of thermal elastic-plastic fem for line heating with remeshingtechnique［J］.Quarterly Journal of the Japan Welding Society，2013，31（4）：134s-137s.

［3］ALSAMHAN A，HARTELY P，PILLINGER I.The computer simulation of cold-roll-forming using FE methods and applied real time re-meshing techniques［J］.JournalofMaterialsProcessingTechnology，2003，142（1）：102-111.

［4］MURTI V，VALLIAPPAN S.Numerical inverse isoparametricmapping in remeshing and nodal quantity contouring［J］.Computers&Structures，1986，22（6）：1011-1021.

［5］黃輝，趙耀，袁華，等.三維網(wǎng)格重劃分技術(shù)及其在焊接數(shù)值模擬中的應(yīng)用［J］.艦船科學(xué)技術(shù)，2011，33（5）：120-125. HUANG Hui，ZHAO Yao，YUAN Hua，et al.Research on three-dimensional rezoning technique and its application in simulation of welding［J］.Ship Science and Technology，2011，33（5）：120-125.

［6］HIBBITT D，KARLSON B，SORENSON P.ABAQUS：analysis user's manual：V6.11［M］.Dassault systemes，2011.

［7］王勖成.有限單元法［M］.北京：清華大學(xué)出版社，2003：666-694.

Rezoning method in nonlinear simulation of plate under moving load

LI Yuantai1，ZHAO Yao1，YUAN Hua1，WANG Bo2
1 School of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology，Wuhan 430074，China
2 Wuchang Shipbuilding Industry Co.Ltd，Wuhan 430064，China

Ship building includes line heating，welding，local roll forming and other processing methods for plates.Compared to the plate dimensions，the region under steady moving linear load in the processing method is narrow with highly nonlinear characteristics.How to coordinate the relationship between the process parameters and the response is one of the main contents of the process.Since the linear load moves on the plate，the characteristics of local nonlinear response keep changing continuously，demanding complicated calculations.Because the mesh for the global region under linear load should be fine，conventional numerical methods result in low efficiency.When considering that other regions away from the loading area in the processing remain weak，nonlinear or elastic，the mesh for the region under current loads can be fine，while the mesh for another region away from the load may be relatively coarser in numerical simulations，with fine mesh moving along with the moving loads，i.e.the rezoning numerical method，which can save much computational time.In this paper，the impact of fine mesh size and rezoning frequency on the numerical computation，and the relationship with numerical accuracy and computational efficiency were researched using this method.The numerical computation example of the local rolling forming process was carried out，resulting in the impact of fine mesh size and rezoning frequency on the computational efficiency and numerical accuracy of the rezoning numerical method.The efficient numerical computation means of plates under moving linear loads with local nonlinear response verify the effectiveness of the method.

moving load；rezoning frequency；fine mesh size

U671.3

A

10.3969/j.issn.1673-3185.2016.05.010

2016-03-04網(wǎng)絡(luò)出版時間：2016-9-21 13:39

國際科技合作專項項目（2012DFR80390）

李元泰，男，1990年生，碩士生。研究方向：數(shù)值計算方法研究。

E-mail：xyt_lee@hust.edu.cn

趙耀（通信作者），男，1958年生，教授，博士生導(dǎo)師。研究方向：結(jié)構(gòu)靜動態(tài)響應(yīng)和企業(yè)信息化研究。E-mail：yzhaozzz@hust.edu.cn