田妍基,陳錦權(quán),陳梅英,3*
(1.寧德職業(yè)技術(shù)學院,福建 福安 355000;2. 福建農(nóng)林大學食品科學學院,福建 福州 350002; 3. 福建農(nóng)林大學管理學院旅游學院,福建 福州 350002)
?
界面厚度對果汁等溫結(jié)晶冰晶生長影響的相場法模擬
田妍基1,陳錦權(quán)2,陳梅英2,3*
(1.寧德職業(yè)技術(shù)學院,福建 福安 355000;2. 福建農(nóng)林大學食品科學學院,福建 福州 350002; 3. 福建農(nóng)林大學管理學院旅游學院,福建 福州 350002)
將果汁體系視為水和溶質(zhì)組成的二元系統(tǒng),采用相場法模擬二元系統(tǒng)的相變微觀結(jié)構(gòu),研究了等溫結(jié)晶過程中,界面厚度對冰晶生長的形貌及溶質(zhì)分布的影響,確定模型中界面厚度的合理取值。結(jié)果表明:界面厚度對冰晶的形貌有正負影響。在擾動干擾的范圍內(nèi),二次枝晶隨著界面厚度的增大而越來越發(fā)達,枝晶尖端也尖銳;在擾動干擾范圍外,晶粒不產(chǎn)生二次枝晶,同時主軸枝晶變細,部分主軸枝晶消失。界面厚度取值越大,晶粒生長速率越快,而枝晶尖端半徑減小。界面厚度取值在3 dx左右能較好地反應枝晶生長形貌及提高計算效率。
相場法;界面厚度;等溫結(jié)晶;冰晶生長
當前, 相場法已成為模擬凝固組織、界面形貌的重要方法[1]。相場法[2]的理論基礎(chǔ)是金茲堡-朗道相變理論,系統(tǒng)的模擬二元系統(tǒng)的相變微觀結(jié)構(gòu),研究了等溫結(jié)晶過程中界面厚度對冰晶生長的擴散、有序化勢以及熱力學驅(qū)動力的綜合作用,這是一種通過該理論的微分方程來反映的一種計算方法。近年來,國內(nèi)外研究者利用相場法模擬凝固微觀組織經(jīng)歷了從二維→三維[3-4]、 二元合金→多元合金[5]、單晶?!嗑Я6-8]、自由枝晶→定向凝固[9-10]、單相場→多相場[11]、無流場→包含流場[12-16]逐步深入的發(fā)展歷程,從而對凝固過程進行真實的模擬[17]。一般情況下,凝固體微觀組織的形成取決于溶質(zhì)擴散、界面曲率與熱擴散[1,18]。迄今關(guān)于利用相場法模擬在果汁冷凍濃縮過程中的冰晶生長情況的相關(guān)研究較少,目前仍處于初步摸索的階段[19-20],亟待進行深入研究。
如今,研究果汁冷凍濃縮尚未構(gòu)建冰晶生長特征的數(shù)學物理模型[21-22],還是主要以實驗為基礎(chǔ),冰晶夾帶的影響因素主要有溶液濃度、降溫速率、結(jié)晶時間等[23-25],另外,還包括處于非平衡態(tài)熱力學條件下的動力學規(guī)律[18,26]。根據(jù)目前已有研究結(jié)果表明,研究者何立群等[27]提出:異相成核后,研究呈現(xiàn)分形生長的冰晶形態(tài)結(jié)構(gòu)和冰晶的演變過程,并基于分形的相關(guān)理論,對水分子的結(jié)晶能力與降溫的速率、溶液的濃度之間的關(guān)系進行分析。陳梅英等[19-20,28-29]首次把相場模型引入食品領(lǐng)域,對冷凍濃縮過程冰晶生長機制進行初步探索,理論分析了冷凍濃縮過冷時間、過冷度等對冰晶生長的影響。
按理界面厚度應取其真實值,但在相場法模擬中需要一個能計算的界面厚度。本研究基于Kim S G等提出的KKS相場模型[26],結(jié)合Allen[30]提出的彌散型界面厚度的定義,探尋能反映真實情況且提高計算效率的界面厚度取值,將果汁的多元成分視為水和溶質(zhì)二元系統(tǒng),對二元系統(tǒng)冷凍濃縮等溫結(jié)晶過程進行相場法模擬,分析冰晶生長過程中界面厚度與冰晶形貌和溶質(zhì)分布之間的關(guān)系。
1.1 相場法的控制方程
冰晶生長相場模型以自由能減小原理為基礎(chǔ),基于kim S G等[26]提出的KKS相場模型,并結(jié)合近似稀溶液時的基礎(chǔ)方程,從而建立了果汁在冷凍濃縮過程中,冰晶生長的相場模擬模型。相場和溶質(zhì)場方程采用自由能密度形式的相場法分別推導出。自由能密度的定義是:利用固相與液相的自由能密度,分別乘以各自的分數(shù),之后再加上剩余自由能的和,可將其表示為:
f(c,φ)=h(φ)fS(cS)+[1-h(φ)]fL(cL)+Wg(φ)
(1)
式(1)中:f為自由能密度,J·m-2;fS表示固相,J·m-2;fL表示液相,J·m-2;且其固相分數(shù)為h(φ),h(φ)=φ3(6φ2-15φ+10),液相分數(shù)為1-h(φ);c為溶質(zhì)濃度;φ為相場變量,φ=-1或0時為液相,φ=1時為固相,在液-固的界面上,φ值在-1~1或0~1的范圍內(nèi)發(fā)生連續(xù)的變化;Wg(φ)表示剩余自由能,其中,W是相場參數(shù),g(φ)是剩余自由能函數(shù),g(φ)=φ2(1-φ)2。可將相場方程表示為:
(2)
式(2)中:M為相場遷移率參數(shù)[式(9)]; f(φ)是自由能密度對相場φ的一階導數(shù);ε為與界面能有關(guān)的相場參數(shù)。當利用稀溶液來近似時,則表達式為:
(3)
式(3)中: R表示氣體常數(shù), R的值為8.31J·kmol-1;T表示溫度,K;Vm表示摩爾體積,L·mol-1;ce表示平衡濃度;下標L和S分別表示液相和固相; h′(φ)為h(φ)的導函數(shù), h′(φ)=30φ2(1-φ)2; g′(φ)為 g(φ)的導函數(shù),g′(φ)=2φ(1-φ)(1-2φ)。
1.2 溶質(zhì)場擴散方程
采用自由能密度的形式,可以將溶質(zhì)場的擴散方程,改寫表示為:
(4)
式(4)中:D(φ) 表示溶質(zhì)擴散率;fc表示自由能對濃度的一階導數(shù);fcc表示自由能對濃度的二階導數(shù)。界面區(qū)域的溶質(zhì)濃度c,是指固、液相的摩爾分數(shù)之和,在固相和液相達到平衡時,在界面區(qū)域中,任意一點的固相和液相的化學勢相等,即:
c=h(φ)cS[1-h(φ)]cL
(5)
μS[cS(x,t)]=μL[cL(x,t)]
(6)
式(6)中μL和μS分別為液相和固相的化學勢。
1.3 加入擾動
計算時,在界面處施加微小的擾動也會影響真實的枝晶形貌,在相場方程中加入人為的隨機擾動:
(7)
式中χ在-1到1之間隨機取值, ?是與時間相關(guān)的相擾動強度因子,16g(φ)用以控制在固液界面中出現(xiàn)擾動,最大擾動可能出現(xiàn)在φ=0.5處,若遠離界面則擾動迅速減小。
2.1 參數(shù)的確定
相場參數(shù)ε和W與界面厚度及界面能都有關(guān)系,相場遷移率參數(shù)M則與界面動力學系數(shù)有關(guān),它們可以表示為:
(8)
(9)
式(8)和(9)中:σ表示固液兩相的界面能,J·m-2;R表示氣體常數(shù);T表示溫度;Vm表示摩爾體積;ke表示平衡常數(shù);μk表示動力系數(shù);me表示液相線斜率;ε表示與界面能有關(guān)的參數(shù), λ表示界面厚度,μm。
2.2 物性參數(shù)的取值
從已有的研究[31-34]得知,冰晶結(jié)構(gòu)在一般情況下,均表現(xiàn)為密排的六方結(jié)構(gòu)。模擬冰晶結(jié)構(gòu)各向異性模數(shù)k的取值為6。當控制其他相關(guān)的物性參數(shù)取值不發(fā)生改變時,只改變界面厚度λ,觀察及分析界面厚度對冰晶形貌的影響。選用質(zhì)量分數(shù)為9.6%的糖水溶液視為水和溶質(zhì)二元系統(tǒng),將其作為研究對象。在進行計算時,所使用的物性參數(shù),具體如表1所示。
表1 蔗糖水溶液的物性參數(shù)
注:表中數(shù)據(jù)參照文獻[35]計算得到。
2.3 數(shù)值計算
(1)邊界條件與初始條件。當計算界面區(qū)域的邊界時,相場與溶質(zhì)場的邊界條件都取絕熱邊界條件。r表示初始晶核半徑,則
(10)
式(10)中: x是模擬中的橫坐標;y是模擬中的縱坐標;T表示有量綱溫度,Tm表示模擬的初始溫度,ΔT表示過冷度。
(2)數(shù)值計算方法。相場方程利用顯示有限差分求解,時間步長Δt受到溶質(zhì)場計算的限制,即Δt<(Δx)2/4DL,DL表示液相中的溶質(zhì)擴散系數(shù),取值為:Δt=1.5×10-9s。晶體軸生長主軸與直角坐標系的x軸和y軸相對應。模擬計算網(wǎng)格數(shù)為1400×1400,空間步長為5×10-8m,初始晶核為半徑r=10個網(wǎng)格數(shù)的球。
在相變過程中, 固體與液體之間存在明顯分界面。將固液界面厚度擴大至可計算的尺度,在不同的界面厚度下利用相場法進行模擬,研究界面厚度這個因素對冰晶形貌與其冰晶溶質(zhì)分布的影響情況,這樣就能獲得穩(wěn)定的界面推進速度。
3.1 界面厚度對冰晶形貌的影響
圖1所示為不同界面厚度下的冰晶生長形貌。當λ=1.5dx時,晶粒形貌呈現(xiàn)等軸枝晶形貌,且二次分支不明顯;隨著界面厚度的增加,當λ=3dx時隨著界面厚度的增加,晶粒開始出現(xiàn)分支,二次分支更加發(fā)達;當λ=4dx時枝晶的主干變細,二次枝晶生長發(fā)達,枝晶整體形貌呈六邊形生長。當λ=4.5dx及λ=5dx時,界面厚度進一步增大,主軸枝晶變得更細,部分主軸枝晶消失,枝晶整體形貌呈球形。這是因為擾動對枝晶的影響存在一定的范圍。在擾動的影響范圍內(nèi),當界面厚度較小的時候,受到擾動的影響較小,二次枝晶不明顯;隨著界面厚度的增大,擾動的影響隨之增大,晶體的界面前沿不穩(wěn)定,二次枝晶開始產(chǎn)生,并隨著界面厚度的增加愈加發(fā)達。當界面厚度增大到一定程度時,能有效地限制擾動的對晶粒的影響減弱,二次枝晶減少最終晶粒沒有二次枝晶的產(chǎn)生,如圖1-d、e。這種現(xiàn)象分別與于艷梅[36]、王智平[37]等模擬的結(jié)果相吻合。
3.2 界面厚度對冰晶溶質(zhì)分布的影響
圖2 顯示了不同界面厚度條件下的溶質(zhì)分布情況,圖1和圖2對比可知相場和溶質(zhì)場相吻合。從述兩幅圖中可以發(fā)現(xiàn),枝晶的固相中夾帶著部分的溶質(zhì),凝固界面的溶質(zhì)濃度比周圍環(huán)境低,且枝晶前沿富集了大量的溶質(zhì)。這是溶質(zhì)擴散的過程。一方面水分子和溶質(zhì)分子都在不斷運動中,水分子之間容易形成氫鍵而結(jié)晶,同時溶質(zhì)分子與水分子相互作用,造成了夾帶;另外一方面,在結(jié)晶過程中,水結(jié)成冰是一個膨脹的過程,因為冰晶生長的速率大于溶質(zhì)擴散的速率,所以在溶質(zhì)來不及充分擴散之前,冰晶內(nèi)的濃縮溶質(zhì)會被迅速擠向尖端,從而在枝晶前端富集,形成一股流行尖端的溶質(zhì)流。由圖2可以發(fā)現(xiàn),隨著界面厚度的逐漸增大,枝晶的生長速率逐漸加快,而枝晶尖端半徑則逐漸減小。同時,從方程(9)可以看出,隨著界面厚度取值的增大,界面的遷移速率減小,溶質(zhì)的擴散較為充分,溶質(zhì)梯度較少,枝晶生長凝固時間越短,所以生長速度越快。所以,方程(9)的理論結(jié)果與相場模擬結(jié)果吻合。
由圖1和圖2可知,當界面厚度取值過大時,枝晶生長形貌變異。所以采用相場法模擬冰晶的生長時,界面厚度取值應在3dx左右,這樣不僅可以得到較可靠的冰晶生長形貌,又可以提高計算機的模擬效率。
4.1 將果汁視作水和溶質(zhì)的二元系統(tǒng),利用相場法模擬二元系統(tǒng)冷凍濃縮過程中的冰晶生長情況,探討界面厚度對冰晶形貌的影響。界面厚度對晶粒的形貌有正負影響,當界面厚度尺度達到一定取值時,其數(shù)值誤差隨之被放大,從而在界面處形成誤差噪音,誤差噪音在界面上形成明顯的擾動,使計算得到的冰晶形態(tài)發(fā)生變化。在擾動干擾的范圍內(nèi),隨著界面厚度的增大,二次枝晶越來越發(fā)達,枝晶尖端也尖銳;在擾動干擾范圍外,晶粒不產(chǎn)生二次枝晶,同時主軸枝晶變細,部分主軸枝晶消失。由于誤差噪音的性質(zhì)難以描述,故利用相場法模擬冰晶生長時,界面厚度大小應有適當?shù)娜≈?,從而達到有效控制誤差噪音。
4.2 相場和溶質(zhì)場相吻合,溶質(zhì)分布情況隨著界面厚度發(fā)生變化。隨著界面厚度取值越大,界面的遷移速率越小,導致溶質(zhì)的擴散速度與冰晶的生長速度相差較小,溶質(zhì)擴散有相對足夠的時間,避免大量溶質(zhì)聚集在固-液界面前沿,削弱了界面前沿出現(xiàn)溶質(zhì)截留現(xiàn)象。因此溶質(zhì)的擴散較為充分,溶度梯度較少,枝晶生長凝固時間越短,所以生長速度越快。
4.3 晶粒生長速率和枝晶尖端半徑受界面厚度取值的影響,界面厚度取值越大,晶粒生長速率越快,而枝晶尖端半徑減小。剛開始,晶粒生長速度呈指數(shù)上升,且隨著界面厚度取值的增大而迅速上升。這是由于開始結(jié)晶的時候伴隨著能量起伏波動,從而影響初始速度,界面厚度較小的時候,擾動對其的影響較少,在擾動干擾的范圍內(nèi),隨著界面厚度的增大,二次枝晶越來越發(fā)達,枝晶尖端半徑減小。因此界面厚度的取值應在3 dx左右,較能真實反映晶粒的生長形貌;界面厚度取值不恰當,將使冰晶形貌發(fā)生變異,也會受到誤差噪音的影響。
4.4 在果汁冷凍濃縮過程中,忽略冰晶結(jié)晶時釋放的潛熱,本研究采用相場法對冰晶生長進行等溫模擬,但模擬的結(jié)果很大程度上受到參數(shù)取值的影響。今后將繼續(xù)在探討熱擴散率、初始晶核半徑等其他他相關(guān)參數(shù)取值方面努力,以取得更優(yōu)的模擬結(jié)果。
[1]HUA HOU,YUHONG ZHAO,YUHUI ZHAO.Simulation of the precipitation process of ordered inter metallic compounds in binary and ternary Ni- Al based alloys by the phase-field model[J].Materials Science and Engineering,2009,499:204- 207.
[2]王潔玉,陳長樂,陳志,等.相場法模擬 Ni-Cu 合金枝晶生長中過冷度的影響[J].鑄造技術(shù),2007,28(4):538-540.
[3]KARMA A, RAPPEL W J. Quantitative phase-field modeling of dendritic growth in two and three dimensions[J].Phys Rev E,1998, 57(4): 4323-4349.
[4]朱昌盛, 馮力, 王智平, 等.三維枝晶生長的相場法數(shù)值模擬研究[J].物理學報,2009,58(11):8055-8061.
[5]ZHANG R J,JING T,JIE W Q,et al.Phase-field simulation of solidification in multi-component alloys coupled with thermodynamic and diffusion mobility databases[J].Acta Mater,2006,54(8):2235-2239.
[6]LI M E,XIAO Z Y,YANG G C,et al.Anisotropic growth of multigrain in equiaxial solidification simulated with the phase field method[J].Chin Phys Soc,2006,15(1):219-223.
[7]SUN Q,ZHANG Y T,CUI H X,et al.Phase field modeling of multiple dendrite growth of Al-Si binary alloy underisothermal solidification[J].China Foundry,2008, 5(4):265-267.
[8]FENG L,WANG Z P,ZHU C S,et al.Phase-field model of isothermal solidification with multiple grain growth[J].Chin Phys B,2009,18(5):1985-1990.
[9]TAKAKI T,F(xiàn)UKUOKA T,TOMITA Y.Phase-field simulation during directional solidification of a binary alloy using adaptive finite element method[J].J Cryst Growth,2005,283(1/2):263-278.
[10]WANG Z J,WANG J C,YANG G C.Phase field investigation on the initial planar instability with surface tension anisotropy during directional solidification of binary alloys[J].Chin Phys B,2010,19(1):492-496.
[11]TIADEN J,NESTLER B,DIEPERS H J,et al.The multiphase-field model with an integrated concept for modeling solute diffusion[J].Physica D,1998,115(1/2):73-86.
[12]TONG X,BECKERMANN C,KARMA A.Velocity and shape selection of dendritic crystals in a forced flow[J].Phys Rev E,2000,61(1):49-52.
[13]CHEN Z,CHEN C L,HAO L M.Numerical simulation of succinonitrite dendritic growth in a forced flow[J].Acta Metall Sin, 2008,21(6):444-450.
[14]龍文元,呂冬蘭,夏春,等.強迫對流影響二元合金非等溫凝固枝晶生長的相場法模擬[J].物理學報.2009.58(11):7802-7808.
[15]袁訓鋒,丁雨田,郭廷彪,等.對流作用下枝晶生長行為的相場法[J].中國有色金屬學報,2010,20(4):681-687.
[16]袁訓鋒,丁雨田,郭廷彪,等.強制對流作用下鎂合金枝晶生長的相場法數(shù)值模擬[J].中國有色金屬學報,2010,20(8):1474-1480.
[17]余德洋,劉寶林,王伯春.超聲波對蔗糖溶液中樹枝冰晶生長速度的影響[J].食品與發(fā)酵業(yè),2011,37(8):92-94.
[18]李剛,劉新田.相場模型及其在凝固組織模擬中的研究進展[J].鑄造技術(shù),2005,26(10):974-976.
[19]陳梅英,陳錦權(quán),陳永雪.荔枝汁冷凍濃縮冰晶生長的動力學分析和相場模擬[J].中國農(nóng)業(yè)大學學報:自然科學版,2011,16(2):153-157.
[20]陳梅英,陳永雪,王文成,等. 相場法模擬[J].福建農(nóng)林大學學報:自然科學版,2010,39(5):548-551.
[21]KIM S G,KIM W T,SUZUKI T.Phase-field model for bianry alloys[J].Phys Rev E,1999,60(6):71-86.
[22]呂???,劉寶林,李維杰.HA納米微粒對PEG-600低溫保護劑反玻璃化結(jié)晶的影響[J].低溫物理學報,2012,34 (4):315-320.
[23]LEVENT B.Mathematical analysis of freeze concentration of apple Juice[J].Journal of Food Engineering,1993(19):95-107.
[24]肖旭霖,李慧.蘋果汁冷凍濃縮工藝的研究[J].農(nóng)業(yè)工程學報,2006,22(1):192-194.
[25]方婷,陳錦權(quán),唐凌,等.橙汁冷凍濃縮動力學模型的研究[J].農(nóng)業(yè)工程學報,2008,24(12):243-248.
[26]GUILLERMO PETZOLD,JOSé M.AGUILERA.Ice Morphology:Fundamentals and technological applications in foods[J].Food,2009,4(4):378-396.
[27]何立群,張永鋒,羅大為,等.生命材料低溫保護劑溶液二維降溫結(jié)晶過程中的分形特征[J].自然科學進展,2002,12(11):1167-1171.
[28]陳梅英,陳永雪,陳威,等.過冷度對冷凍濃縮過程冰晶生長的宏微觀研究[J].西南大學學報:自然科學版,2010,32(5):140-145.
[29]陳梅英,王文成,陳錦權(quán).相場法模擬冷凍濃縮過程冰晶生長的可行性探討[J].江西農(nóng)業(yè)學報,2010,22(3):137-139.
[30]ALLEN S M,CAHN J W.A microscopic theory for anti-phase boundary motion and its application to anti-phase domain coarsening[J].Acta Metall,1979,27(6):1085-1092.
[31]陶樂仁,華澤釗.低溫保存時的氣泡形成現(xiàn)象試驗研究[J].食品科學,2000,21(6):12-15.
[32]陶樂仁,華澤釗.低溫保護劑溶液結(jié)晶過程的顯微實驗研究[J].工程熱物理學報,2001,22(4):481-484.
[33]劉占杰,華澤釗,陶樂仁.脂質(zhì)體懸浮液結(jié)晶對其凍干品質(zhì)影響的研究[J].青島海洋大學學報,2001,31(7):612-618.
[34]陶樂仁.水溶液凍結(jié)過程的顯微實驗研究[D].上海:上海理工大學,2000.
[35]李建國.凝固原理:第4版[M].北京:高等教育出版社,2013:171.
[36]于艷梅,楊根倉,趙達文,等.相場法模擬過冷熔體枝晶生長的界面厚度參數(shù)的取值[J].自然科學進展,2001,11(11):1192-1197.
[37]王智平,張殿喜,石可偉,等.多元合金等溫凝固相場法模擬[J].蘭州理工大學學報,2008,34(6):1-4.
(責任編輯:林海清)
Phase-field Simulation of Ice Crystal Growth by Interfacial Thickness Influence on Juice Isothermal Crystallization Process
TIAN Yan-ji1,CHEN Jin-quan2,CHEN Mei-ying2,3*
(1.NingdeVocationalandTechnicalCollege,Fu′an,Fujian355000,China; 2.CollegeofFoodScience,FujianAgricultureandForestryUniversity,Fuzhou,Fujian350002,China;3.ManagementCollege&Tourismcollege,FujianAgricultureandForestryUniversity,Fuzhou,Fujian350002,China)
The phase-field model which coupled the concentration field and phase field is applied to simulate microstructure evolution during crystal growth of juice binary in 2-D, exploring the process of fruit juice crystallization. In this paper, the influence of interfacial thickness on ice crystal growth morphology and solute distribution was investigated, and a reasonable value of interfacial thickness was estimated. The simulated results showed that the interfacial thickness had positive and negative effects on ice crystal morphology.In the range of disturbance, the secondary dendrite became more and more developed and the dendrite tip was also more sharper with the increase of the interfacial thickness. Outside the range of disturbance, the secondary dendrite didn't grow, and the principal axis of the dendrite became thin and even some disappeared. The growth speed of ice crystal would increase and the tip radius would decreased with the increase of interfacial thickness.When the interfacial thickness was around 3dx, the best simulation result of the dendrite morphology and high calculation efficiency could be obtained.
phase-field;interfacial thickness;isothermal crystallization;crystal growth
2016-05-12初稿;2016-06-05修改稿
田妍基(1981-),女,碩士,講師,主要從事食品加工與營養(yǎng)研究(E-mail:tyj6285268@126.com)
*通訊作者:陳梅英(1972-),女,副教授,博士,主要從事食品工程與旅游食品安全研究(E-mail:cmy2816@126.com)
國家自然科學基金項目(31101327);高等學校博士學科點專項科研基金博導類資助項目(20123515110016);食品生物技術(shù)校級示范專業(yè)和校級優(yōu)秀教學團隊建設(shè)項目(2015)
TS 275.5
A
1008-0384(2016)09-975-06
田妍基,陳錦權(quán),陳梅英.界面厚度對果汁等溫結(jié)晶冰晶生長影響的相場法模擬[J].福建農(nóng)業(yè)學報,2016,31(9):975-980.
TIAN Y-J,CHEN J-Q,CHEN M-Y.Phase-field Simulation of Ice Crystal Growth by Interfacial Thickness Influence on Juice Isothermal Crystallization Process[J].FujianJournalofAgriculturalSciences,2016,31(9):975-980.