塑料食品包裝材料化學(xué)物遷移的數(shù)值方法

陳星玎

(北京工商大學(xué)理學(xué)院,北京 100048)

塑料薄膜這種高分子材料被廣泛用于食品包裝領(lǐng)域,在與食品的接觸過程中,塑料薄膜內(nèi)部存留的添加劑、加工助劑、聚合物單體等化學(xué)物會向食品發(fā)生遷移.這不僅會降低包裝材料對食品的保護功能,還會導(dǎo)致食品受污染而最終危害消費者健康.

20世紀(jì)70年代以來,歐美等國開展了大量包裝材料化學(xué)物遷移的實驗研究.由于遷移實驗過程比較復(fù)雜,需要耗費大量的時間和費用,測試分析更是需要昂貴的儀器設(shè)備,以致在一般的實驗室難以展開.因此,近年來國內(nèi)外的學(xué)者開始研究利用數(shù)學(xué)模型來對遷移進行預(yù)測.

1 化學(xué)物遷移模型

大量的科學(xué)實驗表明,當(dāng)化學(xué)物在塑料內(nèi)的擴散系數(shù)D 以及它在塑料薄膜 食品間的分配系數(shù)K已知時,可基于Fick第二定律對遷移過程進行預(yù)測.通常認為遷移僅發(fā)生在包裝材料的厚度方向上,可用一個一維二階偏微分方程來描述

當(dāng)擴散系數(shù)D與濃度無關(guān)時,式(1)可簡化為

對于模型問題(2),通用的數(shù)值方法包括有限差分法和有限元法.法國Vergnaud等人基于Crank-Nicolson有限差分法數(shù)值求解污染物向原生層遷移的數(shù)學(xué)模型,見文獻[1-3];Roduit等人采用有限元方法進行數(shù)值模擬,見文獻[4];本文將采用譜方法求解模型問題(2),并給出相應(yīng)的理論分析.

為了簡化分析同時考慮實際應(yīng)用情形,通常使用如下假設(shè)條件:1)初始時刻,化學(xué)物均勻分布在包裝薄膜中;2)化學(xué)物從薄膜一側(cè)進入食品,交界面處沒有傳質(zhì)阻力(傳質(zhì)系數(shù)很大);3)任一時刻食品中的化學(xué)物均勻分布;4)在整個遷移過程中,擴散系數(shù)D和分配系數(shù)(KP/F=CP,e/CF,e)為常數(shù);5)遷移過程任何時刻,在包裝薄膜和食品的界面上都是平衡的;6)忽略包裝材料邊界效應(yīng)及其與食品的相互作用.

2 譜方法

譜方法是20世紀(jì)80年代興起的一種數(shù)值求解方法.Quarteroni,Canuto等人對譜方法在理論上做了相關(guān)的研究,見文獻[5-6].譜方法不同于差分法和有限元法,在譜方法中,試探函數(shù)被取為無窮可微的整體函數(shù)(它們一般是奇異和非奇異的Sturm-Liouville問題的特征函數(shù));如果原問題的解充分光滑,譜方法具有“無窮階收斂性”,見文獻[7].

考慮薄膜分單層和多層(雙層、3層、5層等),相應(yīng)的初始條件有所差異.根據(jù)生產(chǎn)實際,假設(shè)初始時刻化學(xué)物均勻布于再生層,未使用的原生層(即原料層)內(nèi)不含污染物.初始條件(單層、多層通式,t=0時)

其中Cin(i)表示初始時刻第i層薄膜內(nèi)化學(xué)物濃度,mg·cm-3;LPi表示第i層薄膜的厚度.邊界條件(t＞0時)

下面給出模型問題(2)的Chebyshev配點法,當(dāng)采用坐標(biāo)變換

可將區(qū)間[LPi-1,Lpi]化為標(biāo)準(zhǔn)區(qū)間[-1,1].所以,在標(biāo)準(zhǔn)區(qū)間[-1,1]上考慮模型問題(2).

設(shè)Tk(x)=cos(karccosx),(k=0,1,…)是定義在[-1,1]上 k次 Chebyshev多項式,取 Tk(x),(k=0,1,2,…,N)為試探函數(shù),那么模型問題(2)的近似解CN(x,t)可表示為

其中ak(t)為待定系數(shù).

令 Tk(x)=cos(kθ),θ=arccosx,有

利用

可得

所以,從

可得

其中

類似的,令

可得

取檢驗函數(shù)φj(x)=δ(x-xj),j=1,2,…,N-1,其中δ(x)為符號函數(shù),xj是(-1,1)上的配置點,一般取.由

可得

因此,模型問題(2)的Chebyshev配點方程為

在式(5)中解出

3 理論分析

本節(jié)給出Chebyshev配點格式(5)的穩(wěn)定性和收斂性結(jié)果.

定義算子L為

其中(·,·)ω為帶權(quán)ω(x)的內(nèi)積.相應(yīng)的離散形式LN為

定理1 Chebyshev配點格式(5)是穩(wěn)定的.

證明:將式(5)寫成變分形式:求CN∈C1(0,+∞;XN),使得

由于

在(6)中取v=CN(t),利用算子L的強制性和連續(xù)性條件,則對 t＞0,有

其中α為強制性條件系數(shù).

于是,在式(7)的兩邊關(guān)于t積分,有

為了討論收斂性,令e(t)=CN(t)-RNC(t),其中RN:→PN上的投影算子,PN為N次正交多項式空間.由原方程和變分形式(6)的第一式,可得

利用算子L的強制性,有

這樣,由式(9)可得

于是,對 t＞0,有

其中c是與N無關(guān)的常數(shù),利用三角不等式有

根據(jù)RN的定義,在式(10)中有

從而可得到如下收斂性結(jié)果.

定理2 當(dāng)N→∞時,近似解CN(x,t)收斂于真解C(x,t).

[1] Perou A L,Vergnaud J M.Contaminant transfer during the coextrusion of food packages made of recycled polymer and virgin polymer layers[J].Computational and Theoretical Polymer Science,1997,7(1):1-6.

[2] Perou A L,Laoubi S,Vergnaud J M.Effect of the thickness of food packages made of recycled and virgin polymer layers coextruded in sandwich form on the time of food protection[J].Computational and Theoretical Polymer Science,1997,8:331-338.

[3] Perou A L,Vergnaud J M.Effect on the time of protection of the dood of the relative thicknesses of the layers[J].Journal of Applied Polymer Science,1999,73(10):1939-1948.

[4] Limm W,Hollifield H C.Modelling of additive diffusion in Polyolefins[J].Food Additives and Contaminants,1996,13(8):949-967.

[5] Canuto C,Hassaini M Y,Quarteroni A,et al.Spectral method in fluid dynamics[M].Berlin:Springer-Verlag,1998:1-568.

[6] Zhang Fa-yong.Spectral and pseudospectral approximations in time for parabolic equations[J].Journal of Computational Mathematics,1998,16(2):107 120.

[7] 向新民.譜方法的數(shù)值分析[M].北京:科學(xué)出版社,2000:1-327.