塑料包裝材料化學物向食品遷移的模擬

李友愛

(北京工商大學 計算機與信息工程學院, 北京 100048)

1 問題的提出

由于具有抗腐蝕能力強、不與酸堿反應、制造成本低、防水、質輕、容易被塑制成不同形狀等眾多優點，塑料已經被廣泛地應用于食品和藥品包裝. 但是這些塑料包裝材料中所含有的化學物可能遷移到被包裝的食品和藥品中，因而會危害人體健康. 如何檢測和分析塑料包裝材料化學物的遷移并形成相關的理論，已經成為食品和藥品包裝領域一個重要的研究課題. 目前，主要檢測和分析手段是實驗的方法. 然而，化學物遷移實驗對實驗儀器要求非常高，且實驗復雜、費用昂貴，很難在普通的實驗室完成. 因此，很多國內外的研究者將研究的重點向理論研究轉移，其中的一個研究熱點就是將化學物的遷移轉化為一個微分方程，然后通過求解這一微分方程模擬化學物的遷移.

大量科學實驗表明(Figge[1]，Chang[2]， Till等[3])， 當化學物在塑料內的擴散系數D以及它在塑料薄膜—食品間的分配系數K已知時，可基于Fick第二定律對遷移過程進行預測. 通常認為遷移僅發生在包裝材料厚度方向上，可用一維二階偏微分方程式

(1)

描述. 若擴散系數D與濃度無關，則上式簡化為

(2)

對于上述遷移模型，只有在如下特殊的情況下才有可能求得解析解[4]:1)初始時刻，化學物均勻分布在包裝薄膜中；2)化學物從包裝薄膜一側進入食品，交界面處沒有傳質阻力(傳質系數 很大)；3)任一時刻食品中的化學物均勻分布；4)在整個遷移過程中，擴散系數D和分配系數KP/F=CP,e/CF,e為常數；5)遷移過程任何時刻，在包裝薄膜和食品的界面上都是平衡的； 6)忽略包裝材料邊界效應及其與食品的相互作用. 一般情況下，我們無法獲得解析方法而只能借助于數 值方法.

對于模型(2),常用的空間離散方法有有限差分[5]、有限元[6]和有限體積方法. 對空間離散后，在時間方向的離散常用的有顯式格式和隱式格式. 用顯式格式計算簡單，但是為了保證數值穩定性，時間步長必須較小；而用隱式格式計算時，可取較長的時間步長，但計算量要大得多. 并且這些方法有個共同的缺點，即時間步長相對較小，且精度不是很高. 鐘萬勰等在文獻[7]中提出了一種精細時程積分法. 這種方法的一個顯著特點就是具有高精度. 本文的目的是將精細時程積分法推廣應用到化學物遷移模型(2)并給出誤差分.

2 精細時程積分法

一般而言，化學物遷移模型(2)經過空間半離散后將得到如下關于時間的常微分方程組：

(3)

其中，v為濃度態變量,H為常系數矩陣,r為非齊次項，表示系統外部作用量. 方程(3)的解v因非齊次項r(t) 而分為兩部分：v1和v2，即

v(t)=v1(t)+v2(t),

(4)

v1(t)=exp(Ht)·v0,

(5)

(6)

其中v0為初始條件.則該系統的解在時刻tk為

(7)

根據方程(7), 解在時刻tk+1(τ=tk+1-tk, 其中τ為時間步長) 為

(8)

通過推導vk和vk+1的關系，可得逐步積分公式

(9)

式(9)為精細積分法的基本公式，是基于解析方法導出的，將初始條件代入上式，即可逐步求出系統在各時刻的解，精細積分法的關鍵是計算上式中的指數矩陣

(10)

為簡單起見，應用二階Taylor展開計算上式中括號內的部分，得

(11)

其中

在式(11)中，高階的Taylor展開值得推薦，會提高結果的精度. 注意到 Δt很小，故Tα相對于I也很小. 如果直接將他們累加會導致有效數位的很大損失. 有鑒于此，精確積分算法先將值較小的矩陣Tα累加，然后將之加到矩陣I. 注意到

(I+Tα)2=I+2Tα+Tα×Tα,

(12)

算法設計為

for (iter=0;iter

(13)

當循環結束后，

T=I+Tα.

(14)

方程(13)和(14)是PTI計算指數矩陣(10)的關鍵步驟.

3 精細時程積分法的誤差估計

若r(t)=0，則方程(3)簡化為

(15)

H為常系數矩陣,其通解可形式地寫成

v=exp(H·τ)v0,

(16)

其中v0為初始時刻的值. 令T=exp(H·τ)=[exp(H·τ/m)]m，又令Δt=τ/m，m=2N(在文獻[7]中，取N=20)，則

T=[exp(H·Δt)]m.

令A=exp(H·Δt)≈I+HΔt+(HΔt)2/2!. 下面考察用精細時程積分法來近似精確解的誤差. 不妨假定H對稱負定，λ1,…,λn為H的特征值，且λi≤0. 又設|λ1|≤|λ2|≤…≤|λn|， 則

H=Qdiag(λ1,…,λn)QT,

(17)

且QQT=I. 令vh=Amv0，于是誤差為

v-vh=(exp(H·τ)-Am)v0=
Qdiag(Δ1,Δ2,…,Δn)QTv0,

(18)

exp(λ1Δt)=1+λ1Δt+(λ1Δt)2/2!+
exp(λ1Δξ1)(λ1Δt)3/3!,
0 <Δξ1<Δt,

(19)

即

1+λ1Δt+(λ1Δt)2/2! =exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!.

(20)

由方程(20)及(18)可知

Δ1=exp(λ1τ)-(1+λ1Δt+(λ1Δt)2/2!)m=
exp(λ1τ) -(exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!)m=
exp(λ1τ) -exp(λ1Δtm)(1-
exp(λ1(Δξ1-Δt))(λ1Δt)3/3!)m=
exp(λ1τ)-exp(λ1τ)(1+x)m=
exp(λ1τ)[1-(1+x)m],

(21)

其中x=-exp(λ1(Δξ1-Δt))(λ1Δt)3/3!. 利用一階Taylor展開可得

(1+x)m=1+m(1+y1)m-1x, 0

(22)

從式(22)可知和(21)式, 可得

Δ1=-exp(λ1τ)m(1+y1)m-1x=
-exp(λ1τ)m(1+y1)m-1(-exp(λ1(Δξ1-
Δt))(λ1Δt)3/3!)=
exp(λ1τ)exp(λ1(Δξ1-Δt))m(1+
y1)m-1(λ1Δt)3/3!=
exp(λ1(Δξ1-Δt))mexp(λ1τ)(1+
y1)m-1(λ1Δt)3/3!.

(23)

下面依次考察式(23)中各項. 由于 0

exp(λ1τ)(1+y1)m-1=
exp(λ1Δt)(exp(λ1Δt)+exp(λ1Δt)y1)m-1≤
exp(λ1Δt)(exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!)m-1=
exp(λ1Δt)(1+λ1Δt+(λ1Δt)2/2!)m-1.

(24)

因為Δt=τ/m,λ1<0， 假設|λ1|≤|λ2|≤…≤|λn|

由上知,(24)式為

(25)

將(25)式代入(23)式得

Δ1≤exp(λ1Δξ1)|(λ1Δt)3/3!|.

(26)

定理1 設v為方程(15)精確解，vh為用精細時程積分法得到的近似解， 則v-vh的l2誤差為

‖v-vh‖l2≤max1≤i≤ncexp(λiΔξi)|λi/m|3τ3‖v0‖l2,

(27)

其中λi,i=1,2,…,n為矩陣A的特征值，m=2N，τ為時間步長，v0為初始時刻的值.

證明：在式(18)的兩邊左乘QT，令u=QT(v-vh)，并注意到Q為正交矩陣，則

u=QT(v-vh)=diag(Δ1,Δ2, …,Δn)QTv0.

(28)

令

y=QTv0=(y1,…,yn)T,

(29)

則

u=diag(Δ1,Δ2, …,Δn)y.

(30)

于是u的離散l2范數為

(31)

因Q為正交矩陣，故有

(32)

‖v-vh‖l2=‖u‖l2,

(33)

即

‖v-vh‖l2≤max1≤i≤n(Δi)‖v0‖l2.

(34)

于是有

‖v-vh‖l2≤
max1≤i≤ncexp(λiΔξi)|λi/m|3τ3‖v0‖l2,

(35)

其中i=1,…,n，c表示不同的常數.