遺傳算法在動態黏彈性模型參數擬合中的應用研究

程軍飛，陳 興，鄧益民，李德群

(1．寧波大學機械工程與力學學院，浙江 寧波315211；2.華中科技大學材料科學與工程學院，武漢 430074)

0 前言

對黏彈性材料本構模型的研究已有很長一段時間[1-2]，但是通過廣義Maxwell模型研究動態黏彈性確定離散松弛時間譜卻鮮有報道[3-4]。在對黏彈性材料進行成型仿真過程中，準確描述黏彈性材料特性的關鍵是廣義Maxwell模型材料參數的準確性以及可靠性。在傳統方法中，通過對儲能模量與耗損模量的頻域曲線擬合模型參數的方法有最小二乘線性回歸法、正則法、非線性回歸法等等，這些傳統方法中都是基于回歸思想求解，在使用上諸多的限制導致產生相對較大的誤差[5-6]。遺傳算法作為一種全局隨機優化算法，能夠快速可靠地求解非常困難的問題，具有固有的并行性和并行計算的能力[7]。

本文通過對TOPAS?COC材料動態黏彈性模型的研究，首先利用遺傳算法對100、110、120、130、140 ℃溫度下的復數模量擬合出黏彈性材料模型參數，與最小二乘法擬合得到的模型參數進行對比，其吻合度更高，表明了遺傳算法在求解黏彈性模型參數上具有很好的優化性。然后利用平移原理,運用WLF方程將不同溫度的復數模量平移到參考溫度80 ℃下，得到TOPAS?COC材料模型的主曲線。最后利用Abaqus有限元分析軟件對滾動熱壓成型進行了仿真模擬，并分析了微結構表面節點在脫模階段中應力和位移變化情況。

1 動態黏彈性模型

1.1 彈性與黏性特性

高分子材料的彈性理論是應力與應變成正比而與應變率無關的理想彈性體的力學性能，對降低生熱、改善屈撓龜裂性以及儲存能量與釋放能量使變形恢復原狀都有很大關系(圖1)。該理論與彈性體Hooke定律相符[8]49-78[9]1-55[10]67-126：

σ=E·ε

(1)

式中σ——應力，MPa

E——彈性模量，MPa

ε——應變

圖1 彈性行為Fig.1 Elasticity behavior

流體力學中黏性研究的是應力與應變率成正比而與應變本身無關的理想黏性流體，使形變物體不能恢復原狀，使輸入的能量轉化為熱能或者摩擦力而消散，如圖2所示。該理論與黏性體牛頓定律[8]49-78[9]1-55[10]67-126相符：

(2)

式中η——黏度，Pa·s

圖2 黏性行為Fig.2 Viscous behavior

1.2 黏彈性模型

在實際的變形過程中，黏性與彈性總是共存的，物體受力時應力同時依賴于形變和形變速率，并且在玻璃化轉變溫度(Tg)附近以上時是介于理想彈性體和理想黏性體之間的黏彈性材料[8]49-78[9]1-55[10]67-126[11]66-69。

本文利用頻域下的黏彈性廣義Maxwell本構模型，結合復數模量的頻域曲線求解材料參數。聚合物黏彈性廣義Maxwell模型是由多支不同的Maxwell模型單元和1支Hooke彈簧并聯組成的，每支Maxwell模型單元具有不同的松弛時間。廣義Maxwell模型在頻域下的儲能模量Es(ω)和損耗能量El(ω)的Prony方程如式(3～4)[12]1-11所示：

(3)

(4)

式中ω——角頻域

Τi——第i支Maxwell模型單元對應的松弛時間

Ei——第i支Maxwell模型單元中Hooke彈簧的模量

(5)

(6)

1.3 時溫等效平移原理

通過大量實驗表明，作用力頻域ω(或時間t)和溫度T對聚合物力學性能的影響存在著某種等同的效果。在較高溫度下聚合物表現出猶如在較長作用時間或緩慢作用力頻域下相同的力學性能；反之，在較低溫度下聚合物表現出猶如在較短作用時間或較快作用力頻域下相同的力學性能。時溫等效的最大意義在于可以通過改變溫度的方法來擴大頻域ω的范圍，使我們用一種試驗方法得到反映聚合物力學性能全貌的整個時間譜或頻域譜，即可把不同溫度下測得的數據轉換成某一參考溫度的組合曲線，如圖3所示，圖中lgG和lgω分別表示模量G和頻域ω的對數函數。

圖3 復數模量曲線的平移Fig.3 Translation of complex modulus curves

描述時溫等效的模型主要使用的是WLF模型，而轉換中最重要的是確定移位因子aT，根據大量實驗數據總結得到aT的經驗公式：

(7)

式中T0——參考溫度，一般選取材料玻璃化轉變溫度50 ℃附近

C1、C2——經驗常數，不同的聚合物材料C1、C2也不同

2 遺傳算法

遺傳算法包括染色體編碼、種群初始化、適應度函數構造、遺傳操作等基本步驟[7]18-38[13]36-59。

2.1 初始化群體和編碼

對于模型參數的研究重點是在一定松弛譜時間(τ1τ2…τL)內對彈性模量E的優化求解，也即是尋求出一組最優的模量(E1,E2…EL),其中L為染色體長度，使{(E1,τ1)(E2,τ2)…(EL,τL)}代入E*后與實測結果吻合度最好。

編碼是把可行解從其解空間轉換到遺傳算法所能解決的搜索空間內的方法，本文采用的是實數編碼方法。

2.2 適應度函數的確定

本文的求解問題為最小值問題，那么作為量化評價的適應度函數通過目標函數F(x)得到:

(8)

2.3 遺傳操作

在選擇操作中，選擇算子采用截斷選擇[7]18-38[13]36-59，截斷閾值根據每次循環的適應度值決定。

在交叉操作中，本文采用算術交叉運算[13]36-59[14]10-11，通常情況下，經過交叉，M個父體產生M個子代，但由于截斷選擇操作過程中適應度在閥值以下的個體被淘汰掉，因此在交叉操作中需要每次隨機抽取2個個體進行交叉，并循環M/2次，直至恢復初始種群數量。

變異操作可以增加種群中染色體的多樣性，改善遺傳算法的局部搜索能力，本文采用非均勻性變異[15-16]，以變異概率對種群個體進行變異操作。

3 遺傳算法對黏彈性模型參數的優化擬合

3.1 遺傳算法的優化性

本文使用的材料為日本Topas Advanced Polymers Gmbh公司開發的一種環烯類共聚物聚高分子材料TOPAS?COC，這種材料具有高于聚碳酸酯(PC)的耐熱性，還由于低吸水性而具有比聚甲基丙烯酸甲酯(PMMA)和PC更加優良的尺寸穩定性，在市場上有很高的評價。圖4為TOPAS?COC材料在不同溫度(100、110、120、130、140 ℃)下的動態熱力學分析實驗得到的復數模量曲線[12]1-11。

溫度/℃：■—100 ●—110 ▲—120 ▼—130 ◆—140圖4 不同溫度下的復數模量曲線Fig.4 Complex modulus curve at different temperature

根據T=100、110、120、130、140 ℃溫度下復數模量 - 頻域的變化曲線，利用遺傳算法分別對5個溫度下的實測模量曲線進行尋優求解，表1為遺傳算法在優化過程中的運行參數。

表1 遺傳算法運行參數Tab.1 Operating parameters of genetic algorithm

圖5～9是在不同溫度下利用遺傳算法和最小二乘法擬合得到的結果與實測結果對比的曲線圖。從圖中可以看出，與最小二乘法相比，遺傳算法得到的結果與實測結果吻合度較好。如表2中所示，相對誤差總和在2.0左右，誤差平均值在15 %左右。而利用遺傳算法得到的結果，相對誤差總和均小于0.5，誤差平均值均在5 %范圍內。

表2 最小二乘法與遺傳算法的誤差結果對比Tab.2 Comparison of error between least square and genetic algorithm fitting

■—實驗 ●—最小二乘法 ▲—遺傳算法圖5 100 ℃溫度下實測與擬合的復數模量曲線Fig.5 Complex modulus curve of experimental and fitted results at 100 ℃

■—實驗 ●—最小二乘法 ▲—遺傳算法圖6 110 ℃溫度下實測與擬合的復數模量曲線Fig.6 Complex modulus curve of experimental and fitted results at 110 ℃

■—實驗 ●—最小二乘法 ▲—遺傳算法圖7 120 ℃溫度下實測與擬合的復數模量曲線Fig.7 Complex modulus curve of experimental and fitted results at 120 ℃

■—實驗 ●—最小二乘法 ▲—遺傳算法圖8 130 ℃溫度下實測與擬合的復數模量曲線Fig.8 Complex modulus curve of experimental and fitted results at 130 ℃

■—實驗 ●—最小二乘法 ▲—遺傳算法圖9 140 ℃溫度下實測與擬合的復數模量曲線Fig.9 Complex modulus curve of experimental and fitted results at 140 ℃

3.2 模量主曲線

利用前面提到的平移原理，對不同溫度下的復數模量數據通過移位因子進行時溫轉換，式(5)中的C1、C2是與材料屬性有關的常數，對于TOPAS?COC材料，可取C1=17.44，C2=51.6。選取TOPAS?COC的Tg=80 ℃為參考溫度，分別將100、110、120、130、140 ℃的曲線向80 ℃平移，如圖10所示。

溫度/℃：■—100 ●—110 ▲—120 ▼—130 ◆—140?—參考溫度(80)圖10 參考溫度80 ℃下的復數模量主曲線Fig.10 Complex modulus curve at 80 ℃

表3是利用遺傳算法求得的黏彈性模型參數(Ei,τi)，圖10為與之對應的復數模量主曲線示意圖，從圖中可以看出擬合得到的復數模量與參考溫度下的復數模量吻合度很高，如表4所示相對誤差總和為1.26，誤差平均值低于5 %。

表3 利用遺傳算法求得參考溫度80 ℃下的動態黏彈性模型參數Tab.3 Dynamic viscoelastic model parameters obtained by Genetic Algorithm at 80 ℃

注：E0=1.8E+1 MPa。

表4 參考溫度80 ℃下擬合結果的誤差Tab.4 Error of fitted results at 80 ℃

表5 參考溫度80 ℃下無量綱的動態黏彈性模型參數Tab.5 Dimensionless dynamic viscoelastic model parameters at 80 ℃

4 TOPAS?COC薄膜滾動式熱壓成型的數值模擬與分析

當聚合物在Tg附近以上時處于橡膠態階段。熱固性塑料在高溫下模壓成型后脫模冷卻至室溫，由于制品脫模后壓力下降，有彈性回復和塑性變形產生使制品的形狀發生變化[17-21]，通過模擬分析方法為熱壓成型提供參考，也為微流道成型質量提供預測幫助，本文針對這一問題研究脫模后彈性回復對聚合物成型圖像的影響規律。

4.1 數值模擬

本文在Abaqus有限元軟件中建立二維模型,如圖11所示，圓柱滾子半徑(R)為100 mm，滾子上的微結構為半圓形，其半徑(r)為0.5 mm，設置2個微弧結構的圓心距離為微弧半徑的6倍。驅動滾子的材料強度通常遠遠大于聚合物材料的強度，因此仿真中模具常常被定義為解析剛體，聚合物基板為可變形體[22-23]。

圖11 滾動熱壓成型制備工藝Fig.11 Rolling hot embossing process

熱滾壓成型過程中，當滾子模具與薄板表面分離時，聚合物處于脫模階段，此時聚合物表面受到的壓力必然減小，受到壓縮后的致密塑料也必然會發生彈性回復致使微結構發生變形，因此成型工藝控制參數對轉印圖形質量具有很大的影響。本文在成型流道上取不同位置的3個點A、B、C，如圖12所示，利用控制變量法研究溫度、預壓緊量及滾動速度對變形情況的影響。

圖12 微流道截面網格圖Fig.12 Grid chart of section of micro-channel

4.2 結果分析4.2.1 溫度對微流道彈性回復的影響

聚合物的物理狀態是某一溫度下的客觀表現，主要隨溫度而變化，因此溫度的高低是影響聚合物流動情況的直接因素。為了研究溫度大小對成型流道彈性回復變形的影響，在預壓緊量為0.45 mm、滾動速度為1.0 mm/s的條件下，選取100、110、120 ℃ 3個溫度分析其應力和位移的變化情況。

如圖13、14所示，相同預壓緊量的同一點位置，彈性回復起始值大小都相同，而衰減的初始值大小卻是溫度越高，初始值越小。在圖13中，A、B、C 3點位置在微結構回彈過程中，溫度越高，應力衰減得越快，并且應力衰減到的最低值也越小，相對應地在圖14中，應力衰減越快，彈性回復位移增長越快，然而位移上升的極限值卻越小。說明在一定范圍內增加溫度，雖然初期彈性回復變化較快，但最終的變化量卻越小，因此增加溫度可以增加聚合物材料在熱壓成型中的黏性變形，降低彈性回復對微結構變形的影響。

溫度/℃： ■—100 ●—110 ▲—120(a)A點 (b)B點 (c)C點圖13 不同溫度條件下A、B、C 3點的應力變化情況Fig.13 Stress at different temperature in A,B and C

溫度/℃： ■—100 ●—110 ▲—120(a)A點 (b)B點 (c)C點圖14 不同溫度條件下A、B、C 3點的位移變化情況Fig.14 Displacement at different temperature in A,B and C

4.2.2 預壓緊量對微流道彈性回復的影響

聚合物在外力作用下發生宏觀變形，材料內部產生附加內力以抵抗外力，當外力卸載時，根據虎克定理，彈性回復距離會隨著內力的變化而變化，由此可見預壓緊量的大小也是壓力的另一種表現形式。為了研究預壓緊量大小對成型流道彈性回復變形的影響，在滾動模具溫度為110 ℃、滾動速度為1.0 mm/s條件下，選取0.4、0.45、0.5 mm 3個預壓緊量分析其應力和位移變化情況。

從圖15中可以看出，在同一位置的預壓緊量越大，應力衰減的初始值越小，并且應力衰減也隨著預壓緊量的增加而加快。發生這一現象主要是由于預壓緊量不同，產生的應變不同，受到擠壓后的材料內部儲存的彈性能也不同，致使在彈性回復階段應力變化不同。從圖16中可以看出，A、B、C 3點位置的位移變化為：預壓緊量越大，彈性回復位移增長得越快，并且彈性回復位移的增長量也越大。

預壓緊量/mm： ■—0.5 ●—0.45 ▲—0.4(a)A點 (b)B點 (c)C點圖15 不同預壓緊量條件下A、B、C 3點的應力變化情況Fig.15 Stress at different preloading height in A,B and C

預壓緊量/mm： ■—0.5 ●—0.45 ▲—0.4(a)A點 (b)B點 (c)C點圖16 不同預壓緊量條件下A、B、C 3點的位移變化情況Fig.16 Displacement at different preloading height in A,B and C

4.2.3 滾動速度對微流道彈性回復的影響

平板熱壓成型過程中有保壓階段，保壓的存在可以使致密制品更好地定型。而在滾動熱壓成型中滾子速度的大小在一定程度上反映了保壓時間的長短，滾動速度越大，滾動模具與薄板的接觸時間越短，保壓時間越短；滾動速度越小，滾動模具與薄板的接觸時間越久，一定程度上起到了保壓的作用。為了研究滾動速度對成型流道彈性回復變形的影響，在滾動模具溫度為110 ℃、預壓緊量為0.5 mm的條件下，選取0.5、1.0、1.5 mm/s 3個滾動速度分析其應力和位移的變化情況。

如圖17所示，微結構脫模時A、B、C 3個節點的應力衰減是隨著滾動速度的增加而加快，出現這種情況主要是因為滾動速度越快，滾動模具與薄板表面的接觸時間越短，起不到保壓的作用，致使材料內部的應力迫切地釋放出去。另外，隨著應力的衰減，3點位置的最低值也逐漸趨向于等同；圖18為A、B、C 3點相對應的彈性回復位移變化情況示意圖，從圖中可以看出，滾動速度越大，位移上升的速度越快。然而在A、B 2點位置處位移上升的最大值隨著滾動速度的增加而增加，在C點位置位移上升的最大值卻是隨著滾動速度的增加而減小，由此可以看出滾動速度對彈性回復位移的影響是隨著應變的減小而減小。因此，當微結構尺寸增大時適當地減小滾動速度可以提高成型制品的質量。

滾動速度/mm·s-1：■—0.5 ●—1.0 ▲—1.5 (a)A點 (b)B點 (c)C點圖17 不同滾動速度條件下A、B、C 3點的應力變化情況Fig.17 Stress at different rolling speed in A,B and C

滾動速度/mm·s-1：■—0.5 ●—1.0 ▲—1.5 (a)A點 (b)B點 (c)C點圖18 不同滾動速度條件下A、B、C 3點的位移變化情況Fig.18 Displacement at different rolling speed in A,B and C

5 結論

(1)利用遺傳算法對黏彈性模型參數進行了尋優求解，以實測與擬合的復數模量相對誤差總和為目標函數，該算法的系統結構采用了實數編碼方法、截斷選擇、算術交叉運算和非均勻性變異，將求得的結果與最小二乘法進行了比較，其擬合結果誤差遠遠小于最小二乘法擬合得到的結果，表明了遺傳算法在求解該模型參數中具有很好的優化性；

(2)利用時溫等效平移原理將100、110、120、130、140 ℃溫度下的復數模量平移到參考溫度80 ℃下，并用遺傳算法對參考溫度下的模量曲線進行模型參數尋優求解，為Abaqus進行熱滾壓成型仿真分析提供了更為精確的模型參數；

(3)基于黏彈性模型在Abaqus有限元仿真軟件建立滾動熱壓成型的模型，并從不同節點處應力和位移變化方面對微流道脫模過程中彈性回復進行了分析，為更好地選擇工藝、控制參數提供了參考。