自動微分在聚合物黏度擬合中的應用

石憲章，王松杰，趙振峰，田 中

（鄭州大學橡塑模具國家工程研究中心，河南 鄭州 450002）

研究開發

自動微分在聚合物黏度擬合中的應用

石憲章，王松杰，趙振峰，田 中

（鄭州大學橡塑模具國家工程研究中心，河南 鄭州 450002）

在聚合物七參數Cross-WLF黏度模型的非線性擬合中，需要解決的一個關鍵問題是如何精確計算梯度矢量Δf2(x)及漢森矩陣Δ f2(x)。因7個黏度參數的數量級差別巨大，傳統的數值微分算法因為容易產生歧義而不再適用，本文采用自動微分方法計算Δf2(x)及Δf2(x)，在忽略截斷誤差的情況下，計算結果即為解析解，該方法不僅有效地保障了Δf2(x)及Δf2(x)的計算精度，而且提高了擬合計算的穩定性，本文算例結果表明擬合誤差（RMS誤差）可達1×10－5量級。

自動微分；非線性擬合；黏度模型；漢森矩陣

聚合物黏度在注塑成型充填、流動模擬中起著非常重要的作用[1-3]，這是因為聚合物具有強剪切特性，其流變行為受黏度變化的影響較大。由于黏度的變化規律復雜，為了精確描述它，進而提高注塑成型充填、流動模擬的精度，黏度數學模型先后經歷了4次改進，由最初的冪律模型到Cross模型、修正的五參數Cross模型以及目前廣泛應用的七參數Cross-WLF黏度模型。七參數Cross-WLF黏度模型的定義如式（1）所示。

式中，η、η0、γ˙、T分別是黏度、零剪切黏度、剪切速率和溫度。n、τ*、D1、D2、D3、A1、A2、為表征黏度的7個參數。其中n為非牛頓指數，τ*為參考溫度（材料的玻璃化轉變溫度）；D1、A1、A2為模型常數；D2對應著低壓下的玻璃化轉變溫度；D3是壓力影響系數，表示黏度的壓力依賴性，對于高聚物材料其值為0。

當7個黏度參數已知時，可根據溫度、壓力和剪切速率來計算黏度值。對于充填、流動模擬來講，溫度、壓力及剪切速率的值是已知的，關鍵的問題是如何得到黏度模型中的7個黏度參數。可通過查閱專業技術手冊得到一些數據，但相當有限，實用性不強，主要是因為黏度隨材料的種類、組分的變化而變化，最有效的方法是通過做實驗，對實驗數據進行擬合得到7個參數的值。文獻[4-7]提出了幾種黏度測試的方法；文獻[8-11]介紹了黏度參數擬合的線性方法，擬合的誤差較大（約1%）；文獻[12-14]采用遺傳算法進行擬合，但遺傳算法的弱點是遺傳和變異的規則難以確定且不完善。 然而由式（1）可知這顯然是一個非線性問題，采用非線性方法進行擬合更為科學合理。

1 非線性擬合方法

為進行黏度參數的非線性擬合，需將式（1）中的溫度（T）和剪切速率（γ˙）當作參數而將其它參數當作未知量，這時黏度可表示為式（2）。

其中，x=[n，τ*，D1，D2，D3，A1，A2]；c=[T，γ˙]。

擬合的目的是獲取表征黏度的7個參數，從而使按式（1）計算得到的黏度與實驗黏度近似相等。為此，按照數值優化理論，可得殘值式（3）。

式中，ηi、φ(x，ci) 是第i次實驗黏度值及對應的擬合黏度值，m是實驗的總次數。

式（3）中ri(x)表示實驗黏度ηi和擬合黏度φ (x，ci)之間的誤差。黏度擬合的基本思想是給7參數選擇合適的值以使函數f(x)取最小值，見式（4）。

這是一個非線性最小二乘問題，對于任意向量p，根據泰勒原理有：

按照線性搜索原理[15]，當選擇了一個搜索方向pk后，在該方向上確定下一步迭代的Xk。這一步的關鍵問題是如何確定xk沿pk移動的距離，方法是通過近似求解式（6）的一維最小值問題計算出步長ak

然后根據式（7）進行迭代計算。

在計算αk時，為減少迭代次數和計算量，希望函數f（x）的值趨于減小并且降幅盡可能的大，為此在求解式（6）時要求αk滿足強Wolfe條件，見式（8-1）、（8-2）。

很明顯，有很多滿足式（6）的下降方向，其中最重要的一個是牛頓方向，通過求解方程得到式（9）。

梯度向量▽f(x)和漢森矩陣▽2f (x)中的元素按式（10）計算。

其中，x=[x1,…,x7]=∣n,τ*,D1,D2,D3,A1,A2∣，j，k=1，…，7。梯度向量▽f(x)和漢森矩陣▽2f (x)的計算在黏度擬合中非常重要，某種程度上它決定著擬合的成敗。這是因為7個黏度參數的數量級差別巨大，數值微分易產生歧義問題，為此本文采用自動微分方法進行計算，在忽略截斷誤差的情況下可以認為它即是解析解。

2 自動微分方法

黏度擬合中采用自動微分方法需要設9個臨時變量，這里用v1～v9進行標識。它們和黏度參數之間的聯系如圖1所示。

圖1 黏度參數和臨時變量的關系

臨時變量的計算公式見式（11）。

臨時變量 v9實質上是黏度η。 基于式（11）可以計算任何臨時變量的一階、二階導數，以v1、v2、v3為例，它們的計算式見式（12）～式（20）。

至此，梯度向量▽f (x)和漢森矩陣▽2f (x)已經得到，可采用數值迭代算法[6]計算求出7個黏度參數。

3 實驗算例

表1給出了一部分實驗黏度數據，試驗材料為PPF401，試驗設備為旋轉流變儀。

圖2所示的5條粉紅色曲線為5種不同溫度下的黏度/剪切速率關系曲線，每一條曲線上的數據都是在相同的溫度條件下測得的，因此5條曲線代表5種不同溫度下的黏度/速率關系。藍色曲線為擬合計算得到的黏度/剪切速率關系曲線，當前的擬合結果為7個參數n、τ*、D1、D2、D3、A1、A2、取初值 0.4丄40000、1×10、400、0、100、200時的結果，顯然和實驗值相差較大。圖3 顯示的是經過6步迭代的中間結果，可以看出擬合曲線與實驗曲線有一定的融合但誤差仍然很大。圖4為經過31步迭代后的最終結果，此時已很難在圖上看到擬合曲線及其數據點，說明擬合曲線已與實驗曲線高度重合，擬合結果與實驗結果高度一致，RMS擬合誤差僅為3.05421 ×10－6。最終得到的7個參數的值為0.3、52000、4.0260399×1012、382.14569、0、32.18244、62.14568。

表1 聚合物（PPF401）流變實驗數據

圖2 擬合前的黏度/剪切速率曲線

圖3 經過6次迭代擬合的黏度/剪切速率曲線

圖4 經過31次迭代后的最終結果

4 結 論

黏度是注塑CAE技術中的一個關鍵參數，由于注塑材料種類的多樣化和復合化，材料的黏度也復雜多變，要得到較精確的黏度參數是非常困難的，而通過實驗獲取黏度數據，再經擬合計算得到表征黏度的7個參數并應用到注塑CAE模型中不失為一種有效的方法。

本文建立了黏度擬合的二階非線性模型，迭代算法采用線性搜索法，較一階線性黏度模型更切合實際，求導計算中采用了自動微分技術，提高了求導的計算精度和擬合計算的穩定性。

算例表明本文所述方法是可行的，擬合結果與實驗結果高度一致，RMS誤差在10－5數量級。

[1] 許百平，瞿金平，任鴻烈. 動態加工過程中聚合物熔體流變行為的研究進展[J]. 塑料工業，2001，23（3）：34-36.

[2] 方傳，彭玉成. PET熔體流變特性的研究與應用[J]. 中國塑料，1998，12（2）：26-30.

[3] 谷諍巍，楊慎華，寇淑清，等. 注射成型中塑料熔體充填流動的黏度模型研究[J]. 模具工業，2002（11）：13- 15.

[4] 吳上生，李群華. 聚合物動態流變儀數據采集系統研制[J].工程塑料應用，2005，33（5）：60-62.

[5] 劉彩文，吳士軍，賈園，等. 基于PC機高聚物流變性能的測試[J].塑料，2006，35（2）：71-74.

[6] 茍剛，謝鵬程，丁玉梅，等. 基于注塑機的聚合物流變特性在線測試[J]. 塑料，2010，39（4）：1-3.

[7] Wagne M H，Bernnat A，Sehulze V. The rheology of the rheotens test[J]. Journal of Rheology，1998，42：917-928.

[8] 李寧，陳靜波，申長雨. 基于G-N法的高聚物粘度模型參數擬合[J]. 鄭州大學學報：工學版，2000，21（4）：55-57.

[9] 王松杰. 射模CAE中粘度模型的參數擬合[J]. 模具工業，2003（8）：13-16.

[10] 徐全軍，李德群.注塑材料流變參數擬合[J]. 塑料科技，1993（3）：19-24.

[11] 張杰，王克立，王鵬駒，等. 注射模CAD/CAE流動模型常數的數據擬合[J]. 模具工業，1998（11）：11-13.

[12] 陶俏，王敏杰，徐斌. 基于Pareto遺傳算法的聚合物黏度模型參數擬合[J]. 塑料，2010，39（2）：7-9.

[13] 周華民，燕立唐，黃棱. 塑料材料的流變實驗與流變參數擬合[J].中國塑料，2001，15（11）：49-52.

[14] 徐全軍，李德群.基因遺傳改進算法用于注塑材料流變參數擬合[J].機械與電子，1994（1）：23-24.

[15]Jorge Nocedal，Stephen J Wright. Numerical Optimization[M]. Berlin：Springer，2006.

[16] Gill P，Murray W，Wright M H. Practical Optimization[M]. New York：Academic Press，1982.

Automatic differentiation and its applications in polymer viscosity fitting technology

SHI Xianzhang，WANG Songjie，ZHAO Zhenfeng，TIAN Zhong
（National Engineering Research Center for Advanced Polymer Processing Technology，Zhengzhou University，Zhengzhou 450002，Henan，China）

The key problem in evaluating the parameters of the 7-parameter Cross-WLF viscosity model of polymer materials with nonlinear fitting method is how exactly to calculate the gradient vector f2(x) and the Hessian matrix. f2(x). Numerical method loses its efficiency here and singularity problem arises due to the significant difference of the parameters in magnitude. In order to overcome this problem，automatic differentiation method was adopted in this paper to calculate f2(x) and. f2(x). Its differentiation accuracy is up to that of analytical method when neglecting the truncation error. Further more，this method is of high stability and benefits much to the final result of the fitting program. Test case shows the RMS error between the fitting viscosity and the experimental viscosity is less than 1×10－5.

automatic differentiation；nonlinear fitting；viscosity model；Hessian matrix

TQ 32

A

1000-6613（2012）06-1298-04

2011-11-10；修改稿日期：2012-03-01。

及聯系人：石憲章（1967—），男，副教授，從事橡塑模具理論研究。E-mail sxely@zzu.edu.cn。