|
摘要：闡述了環(huán)件軋制三維有限元模擬中運用質(zhì)量縮放方法的意義及其理論依據(jù)，并對顯式動力學(xué)有限元模擬中各種質(zhì)量縮放方法進行了分類，給出了環(huán)件軋制模擬中質(zhì)量縮放方法的選取原則。最后基于ABAQUS/EXPLICIT操作平臺，建立了鉛環(huán)件的環(huán)件軋制三維有限元模型，通過具體模擬及分析比較計算結(jié)果對選取原則進行了驗證，從而總結(jié)出了適合于環(huán)件軋制模擬的有效的質(zhì)量縮放方法。
關(guān)鍵詞：環(huán)件軋制 有限元模擬 質(zhì)量縮放
中圖分類號：TG
1 前言
大變形動態(tài)非線性問題（nonlinear problem of dynamic large deformation）[1]的計算和分析可以說代表目前有限元分析的最高水平，也是目前還在處于迅速發(fā)展的領(lǐng)域,其包括的方面主要有：板料成型工藝的數(shù)值模擬（汽車覆蓋件），汽車安全性碰撞的全過程計算，材料的體積成形工藝（擠壓、鍛造、軋制）的全過程數(shù)值模擬等，環(huán)件軋制的三維有限元模擬就屬于其中的一種材料體積成形工藝的全過程數(shù)值模擬。對于此類問題中存在的幾何非線性和材料非線性問題，用于計算求解的方法主要有兩種：靜力隱式算法和動力顯式算法。靜力隱式算法的特征是迭代計算，在每一時間步內(nèi)需反復(fù)迭代，迭代收斂性會受許多因素的影響，在計算中須調(diào)整迭代以滿足收斂，而且其計算時間隨模型單元數(shù)量呈指數(shù)增長，所以其計算時間長。而動力顯式算法是遞推計算，無需迭代，就時間步長直接進行遞推計算，計算時間隨單元數(shù)量呈線性增長，計算時間短，而且可以通過質(zhì)量縮放來縮短計算時間，所以對于大規(guī)模計算和高度非線性問題[2]，動力顯式算法總體效率要比靜力隱式算法高，優(yōu)勢明顯。近年來，隨著在板料成型及碰撞沖擊領(lǐng)域得到日益廣泛的應(yīng)用，對于材料的體積成形工藝領(lǐng)域，人們也都開始采用動力顯式有限元法對其進行過程模擬研究[3]，而且取得了較好的成效。
動力顯式算法在有限元分析中通常被用來解決兩類典型的問題：結(jié)構(gòu)動力學(xué)響應(yīng)分析和準(zhǔn)靜態(tài)分析（其中包括復(fù)雜的非線性效應(yīng)和接觸條件），而環(huán)件軋制的模擬則是其中的一種準(zhǔn)靜態(tài)分析。由于動力顯式算法采用的中心差分法是通過對時間求導(dǎo)進行計算，引入了質(zhì)量矩陣[4]，所以在準(zhǔn)靜態(tài)分析中用于平衡方程式中的離散質(zhì)量矩陣對計算效率和精度都起著至關(guān)重要的作用。從而在這類問題的模擬中，如何正確適當(dāng)?shù)倪\用質(zhì)量縮放方法，使其能夠起到在保持計算精度和穩(wěn)定結(jié)果的同時提高計算效率的作用，也成為了一個至關(guān)重要的問題。目前，針對這個問題，人們也開始對其進行了相應(yīng)的分析研究，東北大學(xué)在利用顯式動力學(xué)有限元法模擬平板軋制過程中，通過計算實例研究了采用提高模擬軋制速度和質(zhì)量縮放技術(shù)來縮短計算時間等方法的有效性，在忽略慣性慣性力影響的前提下得出了與實驗結(jié)果相符合的計算結(jié)果[5]。燕山大學(xué)在用顯式動力學(xué)有限元法分析板帶軋制壓力分布時提出了為了加快計算速度而采用的虛擬速度和虛擬質(zhì)量帶來的虛擬慣性效應(yīng)會影響計算的結(jié)果[6]。而在環(huán)件軋制成形模擬方面，并沒有人通過具體的實例分析來驗證質(zhì)量縮放技術(shù)在環(huán)軋模擬中運用的有效性和合理性。本文運用Abaqus 軟件建立環(huán)件軋制的三維有限元模型，并采用動力顯式算法對鉛環(huán)件的軋制過程模擬進行了實例計算，驗證了質(zhì)量縮放技術(shù)在環(huán)軋模擬中運用的有效性，并通過對計算結(jié)果的分析比較，給出了適合環(huán)件軋制模擬的有效的質(zhì)量縮放方法。
2 質(zhì)量縮放在環(huán)軋模擬中的意義及理論依據(jù)
在Abaqus中利用Explict[7]來進行準(zhǔn)靜態(tài)模擬分析的求解，為了在較快的時間內(nèi)獲得可以接受的靜態(tài)解，以提高求解的效率，通常通過兩種方法來達到目的：一種是人為的提高成型模擬過程的速度；另一種則是通過質(zhì)量縮放來提高模擬速度。對于環(huán)件軋制，模型的材料性能隨應(yīng)變率的改變而變化[8]，人為的提高進給速度就會人為的改變整個變形的過程，使得與實際不相符，而且在實際的環(huán)軋生產(chǎn)工藝中，進給速度最大值受環(huán)件咬入條件的限制[9]，人為的提高進給速度，如果超過了最大值，即使模擬能夠得到穩(wěn)定解，在實際軋制中由于環(huán)件不能滿足咬入條件，軋制過程也是不可能正常進行的。所以人為的提高進給速度會使得模擬與實際脫離，失去了模擬用于驗證和指導(dǎo)實踐的作用，是不可取的。所以對于環(huán)軋的模擬，只能通過質(zhì)量縮放來提高模擬速度。在Abaqus中，質(zhì)量縮放是通過質(zhì)量縮放因子來完成的，質(zhì)量縮放因子取的太小達不到提高計算速度的效果，取的太大則會使得隨著質(zhì)量的增大而隨之產(chǎn)生的虛擬慣性力的增大影響到計算的結(jié)果，甚至使得結(jié)果不能收斂。所以質(zhì)量縮放因子的正確選取也就成為了環(huán)軋模擬中，質(zhì)量縮放方法運用的關(guān)鍵。
動力顯式算法的中心差分法是條件收斂的，其穩(wěn)定的收斂條件是：
其中△t為時間步長，△tcr為臨界時間步長，Tn是有限元系統(tǒng)的最小固有振動周期。在利用中心差分法求解具體問題時，時間步長必須小于由該問題求解方程性質(zhì)所決定的臨界時間步長△tcr，否則算法將是不收斂的。
臨界時間步長通常可以近似認為是應(yīng)力波通過任一單元網(wǎng)格所消耗時間的最小值，即：
式中Lmin為模型中單元的最小長度，與網(wǎng)格劃分的疏密相關(guān)，cd為應(yīng)力波在單元中傳播的速度，對于金屬材料可以表示為：
其中ρ為材料密度，λ和μ為拉梅常數(shù)，可以通過楊氏模量E和泊松比ν來表示：
代入上式則可以表示為成
在Abaqus中用Explict進行求解計算時，所消耗的時間長短與增量數(shù)n有直接聯(lián)系，n可以表示為n=T/△tcr ，其中T為整個模擬過程完成所需要的時間周期（對于環(huán)件軋制模擬，T即為軋制時間）。
 ， 則  ， 因為λ和μ與材料的物理性能相關(guān)，是常數(shù)，由此可以看出n值大小與時間周期T成正比，與單元長度L及ρ值成反比。n值越小，Explict在計算時處理這些增量所消耗的時間也就越短，計算效率也就越高。所以為了減小n值，可以通過減小時間周期和增大單元長度和密度來完成。在環(huán)軋模擬中，減小時間周期即為減小軋制時間，而軋制時間與進給速度成反比，減小軋制時間即意味著增大進給速度，這在前面已經(jīng)說過是不可取的，而單元長度的大小代表著網(wǎng)格的疏密，增大單元長度即為疏劃網(wǎng)格，網(wǎng)格劃的太疏則會降低計算精度，達不到進行精確分析的效果，一般來說也是不可取的。而對于增大密度，根據(jù)塑性變形的體積不變定理，增大密度也就是增大質(zhì)量，質(zhì)量增大倍，根據(jù)上面公式臨界時間步長值△tcr就會增大 倍，相應(yīng)的n值就會減小倍。所以在Abaqus中，通過質(zhì)量放大來達到等效于增大密度從而減小n值的作用，即能保證模擬的穩(wěn)定性和精確度，又能有效的提高計算效率，是最可取的一種方法。
3 質(zhì)量縮放方法的分類及選取原則
前面已經(jīng)提過，動力顯式算法主要用于解決結(jié)構(gòu)動力學(xué)響應(yīng)分析和準(zhǔn)靜態(tài)分析兩種類型的問題，雖然都是用同一種算法，但這兩類問題仍然有很大差別，Abaqus/Explict中用于動力顯式算法的質(zhì)量縮放方法有很多種，但這些方法對于兩種類型問題的分析并不是通用的，大多數(shù)適合于準(zhǔn)靜態(tài)分析的質(zhì)量縮放方法并不適合于結(jié)構(gòu)動力學(xué)響應(yīng)分析。由于環(huán)件軋制模擬屬于一種準(zhǔn)靜態(tài)分析。所以這里只介紹適合于準(zhǔn)靜態(tài)分析的質(zhì)量縮放方法。
Abaqus/Explict所提供的用于準(zhǔn)靜態(tài)分析的質(zhì)量縮放方法主要有兩種類型：
1.通過用戶給定的質(zhì)量縮放因子進行縮放。這是最簡單直接的一種質(zhì)量縮放方法，用戶輸入質(zhì)量想要放大的倍數(shù)即縮放因子值 ，則Abaqus/Explict會將臨界時間步長值放大倍，而且在計算過程的每一時間步中，都會將此時間步的DT值（Abaqus/Explict中以DT表示時間步長，以下均用DT表示）放大倍。
2.通過用戶給定的DT值進行縮放。這種方法就是用戶直接輸入所需求的DT值，而不需要通過給出質(zhì)量縮放因子進行縮放，但前提條件是所給出的DT值不能小于Abaqus/Explict根據(jù)用戶所給出的材料性能參數(shù)計算出的原始DT值，否則就起不到質(zhì)量放大的作用。這里需要說明的是Abaqus/Explict在進行求解計算時會對用戶建模所生成的求解文件進行檢查，并根據(jù)文件中給出的材料性能參數(shù)ρ，λ和μ從而估算給出一個原始DT值，也就是未進行質(zhì)量縮放前的DT值，這個值其實也就是通常所說的臨界步長值，但是Abaqus/Explict為了保證計算的穩(wěn)定和收斂，所給出的原始DT值一般都要比實際臨界步長值要小，對于三維模型，一般為實際值的 /3 ～ 1倍。然后Abaqus/Explict再根據(jù)用戶給出的質(zhì)量縮放方法對DT值進行放大，如果給出的是縮放因子a，則Abaqus/Explict將DT值增大 倍，如果用戶直接給出了DT值，Abaqus/Explict則直接將其替換成用戶給出的DT值，兩個DT值的比值也就是用戶運用這種縮放方法進行放大后，計算時間減少的倍數(shù)。所以用戶在運用這種方法時必須保證給定的DT值不能小于原始DT值。對于這種方法，根據(jù) Abaqus/Explict對DT值進行不同控制的操作，還可以將其分為三種方法。
a 直接根據(jù)給出的DT值進行縮放，對DT值不進行控制。這種方法其實和通過給出縮放因子進行縮放是同一原理。因為給出縮放因子，在表現(xiàn)方法上也是通過給出的因子值，將DT值進行相應(yīng)倍數(shù)的放大。如果用戶給出的DT值是原始DT值的4倍，效果就相當(dāng)于給出的縮放因子值為16，則Abaqus/Explic會將每一時間步的DT值放大4倍。
b 用戶給出DT值，在計算的每一時間步中，Abaqus/Explict自動對每一步的DT值進行控制，使得每一步的DT值在進行縮放后不超過用戶最初給定的DT值。因為在用Abaqus建立模型時，網(wǎng)格劃分采用的是自適應(yīng)網(wǎng)格劃分，那么在模擬過程中，隨著模型形狀變形，Abaqus會自動對模型進行自適應(yīng)網(wǎng)格重劃[10][11]。由于模型各個區(qū)域的變形程度不同，所以重劃后的網(wǎng)格與最初網(wǎng)格相比，其各個單元的單元長度也發(fā)生了變化，有的單元的長度會減小，有的則會增大，隨著模型變形的加劇，進行多次網(wǎng)格重劃后的模型，其各個單元的單元長度變化也會越來越大，有的單元的長度會逐漸減小，有的則會逐漸增大，Abaqus/Explict在確定每一時間步的DT值時，為了保證所有單元的計算收斂，是以這一時刻的應(yīng)力波穿過長度最大的一個單元所消耗的時間為這一時間步的DT值進行計算的。這樣，如果這個DT值不超過臨界時間步長則就能保證這個時間步計算的收斂性。所以在具體確定某一時間步的DT值時，隨著有些單元的長度的逐漸增大，應(yīng)力波穿過最大一個單元所消耗的時間也就越大，那么這一時間步的DT值就會增大，可能會超過原始DT值，甚至有可能超過了臨界時間步長值，那么這一時間步的計算就不收斂了。為了保證計算的精確度，Abaqus/Explict中的自適應(yīng)網(wǎng)格劃分在每一時間步都會對模型網(wǎng)格進行重劃，也就是說整個計算過程的每一時間步的DT值都是在變化的。所以這種縮放方法通過對每一時間步的DT值進行控制，使其不超過原始DT值，以避免出現(xiàn)DT值超出臨界時間步長值，使得計算結(jié)果不收斂的現(xiàn)象，對于大變形的準(zhǔn)靜態(tài)分析是非常合適的。
c 用戶給出DT值，在計算的每一時間步中，Abaqus/Explict自動對每一步的DT值進行控制，使得每一步的DT值在進行縮放后都等于用戶最初給定的DT值。這種方法通過縮放DT值使每一時間步的DT值都恒等于最初的DT值，這樣也就避免出現(xiàn)上面所說的隨著網(wǎng)格重劃而出現(xiàn)的逐漸增大的DT值超過臨界時間步長值使得計算不收斂的現(xiàn)象。所以這種方法對于保證準(zhǔn)靜態(tài)分析過程的收斂性也是很有效的。
以上幾種縮放方法都適合于準(zhǔn)靜態(tài)分析。但是，無論對于通過質(zhì)量縮放因子縮放，還是給出DT值縮放，如果縮放系數(shù)太大，則隨之產(chǎn)生的虛擬慣性力會影響到計算的精度和收斂性，縮放系數(shù)太小則達不到提高計算效率的目的。在準(zhǔn)靜態(tài)分析中[12]，為了保證模擬過程的穩(wěn)定性，一般要求變形的動能值較小，且變化平穩(wěn)，勢能值則較大，動能與勢能的比值較小，一般不能超過0.1，這樣使得能量大部分用于模型的變形。所以一般在準(zhǔn)靜態(tài)分析中，質(zhì)量縮放方法的基本選取原則就是：無論采用何種方法進行縮放，必須保證在整個過程中的動能值要較小，動能與勢能的比值不能超過0.1。在滿足這個前提條件的情況下進行縮放，一般都是可行的。
4 質(zhì)量縮放方法的實例比較
本文以Abaqus軟件為操作平臺，以工業(yè)純鉛為材料，建立鉛環(huán)件軋制的三維有限元模型。并用以動力顯式算法為基礎(chǔ)的求解器Abaqus/Explict進行求解計算。有限元模型如下：
圖1 環(huán)件軋制三維有限元模型圖
工業(yè)純鉛密度ρ=11340kg/m3 ,彈性摸量E=17GPa ,泊松比ν=0.3，最初網(wǎng)格劃分的單元長度為0.01m，根據(jù)上文所給出△tcr公式，將上述值代入求得臨界時間步長理論值為△tcr=4.955×10-6s，而Abaqus/Explict在檢查求解文件時所給出的原始DT值為2.511×10-6s
由此可以看出為了保證計算過程的穩(wěn)定性，Abaqus/Explict所給出原始DT值要比理論值小。
首先采用給出縮放因子進行縮放的方法，來比較不同數(shù)量級的縮放因子，對模擬結(jié)果的影響，本文分別取縮放因子
4、100、1000對環(huán)件軋制過程進行了模擬，環(huán)件最終的成形狀況如圖2所示:
圖2 不同數(shù)量級的縮放因子時的環(huán)件最終的成形狀況
圖2中從左至右依次分別為質(zhì)量縮放因子為4、100、1000時環(huán)件最終的成形網(wǎng)格圖。通過對比可以看出，最左邊的圖效果最好，為圓環(huán)形；中間的圖基本為圓環(huán)狀，但左中部分呈明顯的橢圓形，效果比要差；最右邊的圖效果最差，環(huán)件基本沒有變形左中部分就被壓扁了。由此，單從形狀上看，采用1000的縮放因子是不合理的。縮放因子分別為4和100時，環(huán)件端面寬展變形，環(huán)件變形的動能歷史，及動、勢能比值比較如下：
圖3 環(huán)件端面的寬展變形 圖4 環(huán)件變形的動能歷史 圖5 動、勢能比值
圖3橫軸是環(huán)件的徑向尺寸，從環(huán)件的內(nèi)環(huán)面開始為0，縱軸是環(huán)件端面的寬展距離。從圖中可以看出，縮放因子為4時，環(huán)件沿徑向的寬展變形較平緩，基本符合實際軋制中兩端寬展大，中間寬展小的情況，成形后的端面凹陷不會很明顯。而縮放因子為100時，環(huán)件沿徑向由內(nèi)向外呈明顯的上升趨勢，外端面的寬展要比內(nèi)端面大很多，那么成形后就會形成很明顯的端面凹陷[13]。再比較環(huán)件變形的動能歷史，從圖4可以看出，縮放因子為4時，環(huán)件在整個變形過程的動能值較小，且變化平緩，符合上文所述的準(zhǔn)靜態(tài)分析中，要求有較小且變化平穩(wěn)的動能值的原則，這樣環(huán)件變形過程就會較穩(wěn)定。縮放因子為100時，在同一時刻環(huán)件變形的動能值幾乎都要比前者大很多，而且整個變形過程中動能變化很不平穩(wěn)，升降幅度很大，不符和上述原則，這樣環(huán)件變形過程就會不穩(wěn)定，產(chǎn)生振動或搖擺的現(xiàn)象。最后比較環(huán)件的動、勢能比值，從圖四中可以看出，縮放因子為4時，在變形的初始階段，環(huán)件的動、勢能比值超過了0.1，最大值為0.3，且不穩(wěn)定，有波動現(xiàn)象，到大約0.4秒時比值下降為小于0.1，隨后趨于平穩(wěn)。而縮放因子為100時，動、勢能比值不僅明顯超過了0.1，而且比值較大，最大值達到了1.3，比同一時刻的縮放因子為4時的比值要大很多。在接近1秒時比值才小于0.1趨于穩(wěn)定，不穩(wěn)定時間持續(xù)較長，幾乎為整個過程的1/3，而在穩(wěn)定階段，比值也比縮放因子為4時高。由此，從以上3幅圖的比較可以看出：無論從工藝需求的寬展變形和端面凹陷方面，還是從Abaqus/Explict本身進行準(zhǔn)靜態(tài)分析時對動能輸出及動、勢能比值的要求方面來看，縮放因子取4比取100要更合理。
以上通過對比分析確定了縮放因子取4時符合本文所建鉛環(huán)件模型的模擬需求。下面比較直接給出縮放因子進行縮放即上文所述方法1和通過給出DT值進行縮放，并對DT值進行不同控制的縮放方法即上文所述的方法2中的方法b和方法c（前面已經(jīng)講過方法2中的方法a和方法一其實是一樣的，所以在下面的比較中以方法a代替方法1）這三種方法對模擬成形的影響。模擬過程中仍以縮放因子4為前提，轉(zhuǎn)化為DT值即為放大2倍。采用3種不同方法時，環(huán)件端面寬展變形，環(huán)件變形的動能歷史，及動、勢能比值比較如下：
圖6 環(huán)件端面寬展變形 圖7 環(huán)件變形的動能歷史
圖8 環(huán)件變形的動、勢能比值 圖9 環(huán)件變形的DT值比較
從圖6可以看出，三種方法中環(huán)件端面的寬展變形并無明顯差別，由此可以看出這三中縮放方法對寬展變形的影響不大。從圖7可以看出，雖然三種方法的動能變化都較平緩，但同一時刻方法b的動能值比方法a和c明顯要小，方法a和c的動能值變化則基本一致。圖8為三種方法的環(huán)件變形動、勢能比值比較，從圖中可以看出，方法a和c剛開始階段比值都超過了0.1，但c相對a要好，最大比值為0.17，而a最大值為0.31，c的比值波動的幅度要比a小，比值降為小于0.1時所經(jīng)歷的時間也要比a短。方法b的效果則最好，環(huán)件變形的整個過程中動、勢能比值都要小于0.1，最大值為0.095，而在隨后的穩(wěn)定階段，同一時刻方法b的動、勢能比值也要比方法a和c的小。圖8為計算過程中每個時間步的DT值輸出。由于Abaqus/Explict給出原始DT值為2.511×10-6s，而最初都是將DT值放大2倍即為5.022×10-6s，所以從圖中可以看出三種方法剛開始的DT值都為5.022×10-6s。方法a不對DT值進行任何控制，從圖中可以看出有些時間步的DT值明顯大于最初值，雖然在本模型的分析中，這些較大的DT值并沒有超過實際的臨界時間步長而導(dǎo)致計算結(jié)果不收斂，但在其它模型的分析中，如果Abaqus/Explict所給出的原始DT值與實際臨界時間步長值很接近時，那么運用方法a進行縮放時，較大的DT值就可能會超過實際臨界時間步長值而使得計算結(jié)果不收斂。方法b是控制DT值，使每步的DT值都不超過最初值，從圖中也可以看出，整個過程的DT值處在一個下降的趨勢中，且都不超過5.022×10-6s。方法c是控制DT值，使每步的DT值都等于最初值，從圖中明顯可以看出，整個DT值變化為一條直線都等于5.022×10-6s。另外，對于使用三種縮放方法進行模擬時，記錄CPU計算消耗時間分別為1小時39分，1小時40分和1小時42分。可見三種縮放方法計算效率幾乎一樣。
從以上四幅圖的分析比較中可以看出方法b既符合準(zhǔn)靜態(tài)分析對動能輸出及動、勢能比值的要求，又能保證有效的計算效率，所以是三種方法中最適合本模型分析的質(zhì)量縮放方法。
總結(jié)以上的分析比較可以得出：準(zhǔn)靜態(tài)分析中質(zhì)量縮放方法的選取原則對環(huán)件軋制的模擬是適用的。對于環(huán)軋模擬中縮放方法的運用，首先應(yīng)該從環(huán)件成形的形狀、寬展變形、動能變化及動、勢能輸出等方面進行比較，確定一個較合適的縮放因子，然后再運用方法b，通過對每一時間步的時間步長值進行控制，保證其不超過臨界時間步長值。這樣，即提高了計算效率又保證了計算過程的穩(wěn)定性和收斂性。
5 結(jié)論
本文采用不同的質(zhì)量縮放方法及放大系數(shù)對鉛環(huán)件三維有限元模型進行了實例分析，通過分析比較不同縮放系數(shù)時的模擬結(jié)果，驗證了在準(zhǔn)靜態(tài)分析中，為了保證模擬過程的穩(wěn)定性，質(zhì)量縮放方法的選取必須遵循：通過縮放后的整個模擬過程其變形的動能值應(yīng)較小，且變化平穩(wěn)，勢能值則較大，動能與勢能的比值應(yīng)較小，一般不能超過0.1，這一選取原則對環(huán)件軋制模擬是適用的；通過分析比較采用不同縮放方法時的模擬結(jié)果，總結(jié)了一套適用于環(huán)件軋制模擬的最有效的質(zhì)量縮放方法。
|