茄子冷冻干燥冻结过程模型的建立与求解.pdf

第 26 卷第 5 期农 业 工 程 学 报Vol.26No.5 2010 年5 月Transactions of the CSAEMay 2010335 茄子冷冻干燥冻结过程模型的建立与求解 郝新生，郭玉明，崔清亮 （山西农业大学工学院，太谷 030801） 摘要：该文以茄子球形样本为对象，研究果蔬物料冷冻干燥中预冻模型的建立及求解方法。通过能量衡算建立了果蔬 冻结过程的非稳态数学模型，依据冻结过程试验参数，给出了冻结过程中物料液相区和固相区温度变化所满足的定解问 题。采用“欧拉向后差分格式”数值计算方法，对冻结模型进行了数值求解，得出了茄子在冻结过程中冻结界面的移动 规律和冻结时间。试验结果表明，所建模型能够准确地描述茄子的冻结过程，采用的数值解法是可行的。该文研究方法 及结果为果蔬物料冷冻干燥预冻结过程模型的建立与求解提供了参考，并为果蔬冷冻干燥预冻结过程的优化控制提供了 科学依据。 关键词：冷冻，低温干燥，模型，数值分析，stefan 问题，果蔬 doi：10.3969/j.issn.1002-6819.2010.05.057 中图分类号：TQ028文献标识码：A文章编号：1002-6819(2010)-05-0335-07 郝新生，郭玉明，崔清亮. 茄子冷冻干燥冻结过程模型的建立与求解J. 农业工程学报，2010，26(5)：335341. Hao Xinsheng, Guo Yuming, Cui Qinliang. Modeling and numerical solution of eggplant freezing during freeze- drying processJ. Transactions of the CSAE, 2010, 26(5): 335341. (in Chinese with English abstract) 0引言 冷冻干燥是将含水物料冻结后，在适当的温度和真 空度下，使物料内部的冰直接升华为水蒸汽并逸出，从 而获得干制品的一种方法，其基本过程分为：冻结、升 华干燥和解吸干燥 3 个阶段。与其它干燥方法相比，冷 冻干燥得到的冻干制品能最大限度地保持物料原有的形 态与营养成分1- 2。但冷冻干燥所需时间长、能耗大，导 致产品成本高（干燥成本高出热风干燥的 48 倍） ，制 约了该技术在农产品加工领域的大规模应用。因此，缩 短冷冻干燥时间和降低能耗成为该领域的一个重要研究 方向，Sadall、Lichtfield 和 Liapis 等众多学者主要针对升 华干燥与解吸干燥阶段建立了相应的数学模型，对其传 热传质机理和工艺参数的操作控制进行了研究和探讨， 取得了许多有益成果3- 8。但用数学模型来描述和模拟冻 结阶段的研究则很少，本文采用球形果蔬物料物理模型， 建立了非稳态数学模型来描述和模拟果蔬物料的冻结过 程，为果蔬物料冻结过程优化控制提供理论依据和过程 控制模型。 1球形果蔬物料冻结模型的建立 果蔬物料冻结的最终温度以其共晶点为依据，物料 收稿时间：2009- 09- 26修订时间：2010- 04- 26 项目基金：国家自然科学基金资助项目（30771242） ；山西省科技攻关项目 （20080311051） 作者简介：郝新生（1968） ，男，山西武乡人，副教授，博士生，中国农 业工程学会会员（E041200466S） ，主要从事农产品加工过程的模拟和控制 的研究。太谷山西农业大学文理学院，030801。Email: haoxs001 通信作者：郭玉明（1954） ，男，山西平定人，教授，博士生导师，中 国农业工程学会高级会员（E041200306S） ，主要从事农业生物力学与物料 机械特性的研究。太谷山西农业大学工学院，030801。 Email: guoyuming99 的共晶点（或共晶温度）是指物料中的水分全部冻结时 的温度。为保证物料完全冻结，预冻温度要低于物料的 共晶点 510。若预冻温度过低，则会延长冻结时间， 增加能耗和生产成本；若预冻温度高于共晶点，则不能 保证物料中的水分完全冻结，物料内部水分将不能完全 以冰的形式升华，导致物料在干燥过程中发生收缩和失 形等问题。另外，未冻结的水分中所含溶质，在干燥过 程中可能随内部水分向物料表面迁移，出现冻干制品表 面硬化现象2。 物料冻结又可分为 3 个阶段：第 1 阶段为单纯的降 温过程，无相变发生，当物料表面温度降至水的结晶温 度时，开始形成冰晶，认为第 1 阶段结束；第 2 阶段为 结晶过程，完成相变，使物料内部水分几乎完全冻结； 第 3 阶段为再降温过程，直至物料中心处温度降至物料 共晶温度以下，可视为无相变发生。下面将根据热量守 恒定律来建立物料冻结阶段的温度变化所满足的偏微分 方程。为此需对物料作如下的基本假设： 1）物料被视为半径为 R 的理想球体； 2）物料内部各向同性，各物性参数均为常数； 3）传热沿半径方向； 4）物料在冷冻过程中体积无明显变化。 1.1球状果蔬物料降温阶段模型的建立 将物料置于温度为 T*的低温环境中，物料共晶点温 度为T。 取物料的中心点（球心）为坐标原点 O，建立空间 直角坐标系，如图 1 所示。 令T1(x,y,z,t)为t时刻物料内部点(x,y,z)处的温度， 1、c1、k1分别表示物料未冻结时的密度、比热容和导热 系数。 336农业工程学报2010 年 图 1球形物料直角坐标示意图 Fig.1Schematic diagram of rectangular coordinates of spherical materials 在物料内任取一块微元体 V，并用V表示其边界曲 面，S为V上的任一面积微元，其单位外法向量为n ， 则根据 Fourier 实验定律9，在时间段 t 内流经 S 的热 量为 1 1 T QkS t n 则从 t1到 t2时间内通过边界V流出的热量为 2 1 1 11 (d )d t t V T QkSt n 利用 Guass 公式 2 1 111 (d d d )d t t V QkT x y zt 另外，t1到 t2时间内 V 中温度降低放出的热量为 2111112 ( , , , )( , , ,)d d d V QcT x y z tT x y z tx y z 2 1 1 11 d d d d t t V T ct x y z t 由热量守恒律得，Q1=Q2，即 22 11 1 1112 (d d d )dd d d d tt tt VV T kT x y ztct x y z t 由 Fubini 积分顺序交换定理9及 t1,t2,V 的任意性，可得 222 21111 1 222 TTTT a txyz （1） 其中， 21 1 11 k a c 。 由于研究对象看作球体，故采用球面坐标（r, , ） 来表示方程及相应的边界条件，直角坐标（x,y,z）与球 面坐标（r, , ）的关系为 sincos sinsin cos xr yr zr 方程（1）在球面坐标系中可表示为 22 111 1 22 11 sin sin TTT ar trrrr 2 1 222 1 sin T r 根据假设条件，物料内温度分布是球对称的，即 11( , ) TT r t 这样在球面坐标系下，物料降温过程温度 T1所满足的微 分方程为 2 2111 1 2 2TTT a trrr (0,0)rR t（2） 式中 R球形物料半径。 1.2球状果蔬冻结阶段模型的建立 冻结过程中，物料内部从外向内逐渐生成冰晶，在 物料内部形成固、液两相区域，如图 2。在固相和液相区 中， 其传热机理是相同的， 只是对应的物性参数有所不同。 注：R球形物料半径；s(t)物料中心到冻结界面的距离 图 2球形物料冻结过程物理模型 Fig.2Physical model of freezing process of spherical material 1.2.1液相区温度变化所满足的热传导方程 记 T1(x,y,z,t)为 t 时刻物料液相区点(x,y,z)处的温 度，则(x,y,z,t)满足 2 2111 1 2 2TTT a trrr (0( ),0)rs tt（3） 其中，s(t)t 时刻冻结界面到物料中心的距离， 0( )s tR。 1.2.2固相区温度变化所满足的热传导方程 记 T2(x,y,z,t)为在 t 时刻物料固相区中点(x,y,z)处 的温度，则 T2(x,y,z,t)满足 2 2222 2 2 2TTT a trrr ( ( ),0)s trR t（4） 其中， 22 2 22 k a c ， 2 、 2 c、 2 k分别表示物料固相区的 密度、比热容和导热系数。 1.2.3冻结界面温度变化所满足的微分方程 冻结界面是液相区和固相区的过度区域，是水分开 始凝固的场所。通常情况下，冻结界面很薄，认为冻结 界面上不能积累能量，在任何时刻通过液相区和固相区 传递到冻结界面处的热量都等于水分凝结所放出的热 量。对于球形物料其冻结面可近似看作为一个球面。 从时刻 t 至 t+t，对冻结界面作热量衡算，可得 12 12 d ( ) ( ),0) d TTs t LCkkrs tt trr （5） 式中：L水的凝固潜热；C物料的含水率。 第 5 期郝新生等：茄子冷冻干燥冻结过程模型的建立与求解337 1.3初始条件和边界条件 初始条件 10200tt TTT （6） 边界条件 2 22 0 12 1 , ,( ) 00 T kTTrR r TTTrs t T r r （7） 式中：T0物料的初始温度； 物料表面与低温环 境的热交换系数；T0物料中水分的结晶温度。 2果蔬物料无相变降温阶段模型的求解 物料降温阶段温度变化所满足的定解问题 2 2111 1 2 100 1 11 2 ,0 t r R TTT arR trrr TT T kTT r 记 1 ( , )()u r tr TT ，则上述定解问题转化为 2 2 1 2 * 00 0 ,0 () 0 0 t r r R uu arR tr ur TT u u u r 其中， 1 1 kR . 令( , )( )( )u r tN tr，代入方程可得 2 1 Na N，即 2 1 N a N . 式中：N(t)关于 t 的待定函数；(r)关于 r 的待 定函数；比例系数。 由此可得 2 1 0 0 NaN 由边界条件可得 (0)0,( )( )0RR 解特征值问题 0 (0)0 ( )( )0RR 可得其特征值为 2 n ， 特征函数为 sin( nr)，n=1,2,3, 其中 n为方程 xcot(Rx)=-的第 n 个正根。 当 = n时，记 N=Nn， 则由 22 1 0 nnn NaN ，可得 22 1 e n at nn ND 式中 Dn依赖于 n 的系数。 这样， 22 1 1 ( , )expsin() nnn n u r tDatr 。 由初始条件可得 0 1 (*)sin() nn n r TTDr 其中， 0 0 2 0 ()sin()d sin ()d R n n R n r TTrr D rr 22 0 22 2()(1) () n nn TTR R 记 0 t为物料表面达到水的结晶温度所用时间，则 0*22 10 1 1 exp sin() nnn n TTDatR R 式中所含级数的收敛速度很快，在工程上取其第一项即 可达到精度要求，即： 22 0*11101 exp sin()DatR TT R 从而 0* 0 22 1111 1() ln sin() TTR t DRa （9） 在 0 t时刻物料内的温度分布为 22 1010 1 1 ( ,)expsin() nnn n T r tTDatr r 22 11101 expsin()Datr T r 32 221 11110 ()exp 6 r TDat 记 2 10 ( )( , )rT r tABr，（10） 其中， *22 1 1110 expATDat， 3 221 1 110 exp 6 D Bat 。 3球形果蔬物料冻结过程模型的求解 3.1移动界面的固化 在冻结过程物料温度变化所满足的定解问题中， 其求 解区域的部分边界是移动的， 它和问题的解有关并且要同 时被确定，这类问题在数学上被称为自由边界问题或 stefan 问题10- 14。在上面建立的冻结模型中液态区和固态 区温度变化所满足的定解问题一起构成 1 个两相的 stefan 问题，其精确解根本无法求出，下面采用区域边界固定 法10对其进行数值求解，得到满足一定精度的近似解。 338农业工程学报2010 年 首先引入 2 个无量纲变量将移动界面固定： 令, 0( ) ( ) r rs t s t ， ,( ) ( ) rR s trR s tR 在此变换下，模型方程中出现的偏导数相应地变为 11 , , 1 ( ) t r t TT rs t ， 22 11 222 , 1 ( ) t rt TT rst ， 111 , , ( ) ( ) t rt t TTTs t tts t ， 22 , , 1 ( ) t r t TT rs tR ， 22 22 222 , 1 ( ( ) t rt TT rs tR ， 222 , , ( ) ( ) t rt t TTTs t tts tR 。 这样，液相区和固相区温度变化所满足的微分方程 可化为 22 11111 22 2( ) ( )( ) TaTTTs t ts tst (01,0)t（11） 22 2222 22 2 ( ( ) ( ) TaTT R ts tR s tR 2 ( ) (01,0) ( ) Ts t t s tR （12） 冻结界面方程和边界条件转化为 1122 d ( ) d( )( ) kTkTs t LC ts ts tR (1,1,0)t（13） 22 2 ( ) kT TT s tR (0)（14） 1 0 T (0)（15） 0 12 TTT(1,1)（16） 3.2微分方程的离散化 对上述方程采用具有无条件稳定的向后Euler格式进 行求解 15：在固、液两相区，取空间步长均为 1 h m （m 为正整数） ，时间步长为 进行网格刨分，结点为( ,) ik t （液相区） 或(i,tk) （固相区） ， 其中 i=ih； i=ih； tk=k ； 1,2,3,1im，1,2,3,k 。 先在结点( ,) ik t处考虑方程（11） 22 1111 22 2 ( ,)( ,)( ,) ( ) ikikik ik TaTT ttt ts t 1 ( ) ( ,) ( ) ik ik k s tT t s t （17） 记 1( , ) k iki Ttu， 式中： k i u液相区结点( ,) ik t处温度。 则由 Taylor 公式可得 11 2 ( ,)( ,)(01) 2 k iktiik TT tD ut tt （18） 3 11 3 1 ( ,)()(,) 26 kk ikiiikk TTh tD uD ut 11 () iiki （19） 242 211 24 ( ,)(,) 12 k ikiikk TTh tut 11 () iiki （20） ( )( )( )() 2 ktkkk ikk s tD s ts tttt （21） 式中： t Dt的一阶导数向后差分格式；D的 一阶导数向后差分格式； 2 的二阶导数中心差分 格式。 将（18）（21）代入（17）中，并舍去无穷小量 项，可得如下差分格式 2 21 2 ( )1 ()() 2 ( )( ) kkkkkkitk tiiiiii kk D s ta D vuD uD uD uD u s ts t （22） 在结点(,) ik t处考虑方程（12） ，并记 2( ,) k iki Ttv， 可得方程（12）差分格式 2 2 2 2 ( ) () 2( ( )( ( ) 1 () ( ) kkkitk tiii kk kkk iii i k D s ta D vD vD v s tRs tR vD vD v R s tR （23） 依次将式（13）（16）离散化后可得 00 11 121 ()( ) ( )( ) kk mm kk kk TvTu s ts tkk LChs tRs t （24） 22 22 0 22 2 22 122 10 22 2 2211 1 ( ( )( ( ) 2211 ( ( )( ( ) k kk kk kk aa v kRh s tRh s tR aa vvT kRh s tRhs tR （25） 22 111 010 2222 66 1 ( )( ) kkk kk aa uuu h sth st （26） 0kk mm uvT（27） 将关于 k i u和(0,1,2,) k i vim的差分方程表示成矩 阵形式，可得 第 5 期郝新生等：茄子冷冻干燥冻结过程模型的建立与求解339 0 11 1 22 2 2,2, 2 1, 1 3233 k kk k kkk k kkk k mkkmk m k mkk m bb u abc u abc u abc u ab u 1 0 1 1 1 2 1 2 10 11, k k k k m k mmk u u u u ucT （28） 0 11 1 22 2 2,2, 2 1, 1 1 3 k kkk k kkk k kkk k mkkmk m k mkk m bdb v abc v abc v abc v ab v 1* 0 1 1 1 2 1 2 10 11, k k k k k m k mmk vd T v v v vcT （29） 其中， 2 11 22 ( ()( )1 (1) 2 ( )( ) kk ik kk i s ts ta a is th st ， 2 1 22 2 1 ( ) k k a b h st ， 2 11 22 ( ()( )1 (1) 2 ( )( ) kk ik kk i s ts ta c is th st ， 22 22 222 1 ( ( ) ( ( )( ( ) ( ()( ) 2( ) ik kkk kk k aa a hs tRh ih s tRR s tR i s ts t Rs t ， 2 2 22 2 1 ( ( ) k k a b hs tR ， 22 22 222 1 ( ( ) ( ( )( ( ) ( ()( ) 2( ) ik kkk kk k aa c hs tRh ih s tRR s tR i s ts t Rs t ， 2 1 2 211 ( ) k k a d kRh Rs t 。 3.3对离散后的差分方程采用 Matlab7.1 程序 数值求解16，算法如下： 1）取 0 ( )s tR， 11 ( )s tR（ 1 为任意的小正数） ， 022 1 ( ) () i uAB ihR，(0,1,2,1)im， 02 1 () i vAB Rih，(0,1,2,1)im， 由式 （25） 和 （26） 可求得 1 i u和 1 i v(0,1,2,1)im， 将 1 1m u 、 1 1m v 代入式（21） ，可求得 2 ( )s t； 2）将 2 ( )s t、 1 i u和 1 i v代入式（25） 、 （26）可求出 2 i u、 2 i v(0,1,2,1)im，再将 2 1m u 和 2 1m v 代入式（21） ，可 求出 s(t3),； 3）如此反复迭代，直到 s(tk)2（ 2 为任意的小正 数） ，迭代停止。 4）迭代次数与时间步长的乘积即为冻结阶段所需时 间，记为 1 t。 4球形果蔬冻结后继续降温阶段所需时间的计算 这一阶段温度变化所对应的定解问题为 2 2222 2 2 200 2 22 2 ,0 t r R TTT arR trrr TT T kTT r （30） 式中 0 T 物料中水分全部冻结时的温度分布。 采用分离变量法可得上述定解问题的解为 *22 22 1 1 ( , )exp sin() nnn n T r tTBatr r 其中， *22 0 22 2()(1) () n n nn TTR B R ， 2 1 kR ， n 为方程cot()xRx 的第 n 个正根。 物料中心温度随时间变化的关系为 22 0 *22 2 0 1 *22 2 1 (0, )lim( , ) 1 lim(exp sin() exp r nnn r n nnn n TtT r t TBatr r TBat 记 2 t为物料中心达到T 所用时间，其中T 为共晶温度T 以下的某一温度值，则 *22 22 1 exp nnn n TTBat 式中所含级数的收敛速度很快，取其第一项即可达 到精度要求，即： *22 11212 expTTBat 由此可得 * 2 22 1121 1 ln TT t Ba （31） 这样整个预冻结过程所用时间为 012 tttt 340农业工程学报2010 年 5计算结果与验证实验 5.1仪器、材料及茄子的物性参数 5.1.1主要测量仪器 SE1501F电子天平（奥豪斯上海有限公司）， DHG- 9023A电热恒温鼓风干燥箱（无锡三鑫精工电气设 备有限公司），DW- 40L188立式低温冰箱（青岛海尔集 团），JDG- 02试验冻干机（兰州科近真空冻干技术有限 公司），Hotdisk热物性测试仪（瑞典Hotdisk AB公司生 产）。 5.1.2实验材料及物性参数 选择成熟适度、直径大约10 cm的新鲜球形茄子。将 茄子去皮，测得茄子的物性参数，并作为模型计算参数： C=94.58%， =700 kg/m3，R=0.05 m，c1=4 000 J/(kgK)， c2=2 040 J/(kgK)， k1=0.58 J/(mKs)， k2=1.66 J/(mKs)， =2 J/(mKs)， L=3.344105J/kg， T0=269 K，261 TK， T*=243 K，T 253 K，T0=299 K。 5.2数值计算结果 5.2.1 0 t和 2 t的计算 首先采用二分法求解代数方程)cot(Rxx和 )cot(Rxx， 通 过 一 简 单 的 C 程 序 可 求 得 ： )cot(Rxx的前 5 个正根为 16.89，90.98，155.15， 218.54，281.68；)cot(Rxx的前 5 个正根为 18.59， 95.48，162.41，229.55，302.87。再由式（9）和（31） 可得， 0 53 mint ， 2 42 mint 。 5.2.2 1 t的计算 取 1=2=10- 5， =0.4，h=0.05，并将有关参数代入 Matlab 程序计算得，迭代次数为 k=16853，则 1 16853t 0.4s111 min。 因此冻结过程总时间min206t。由计算结果可绘 出冻结界面和物料中心温度随时间变化曲线， 分别见图 3 和图 4。 注：s(t)t 时刻冻结界面到物料中心的距离 图 3冻结界面随时间变化曲线 Fig.3Variation curve of freeze interface in accordance with time 5.3实验验证与分析 将低温冰箱冷冻室内的温度设定为-30，把温度传 感器（铂电阻式）插入茄子中心，然后置入低温冰箱内。 开启冷冻干燥试验机的调控系统，采集茄子中心温度变 化数据，采样时间间隔为 1 min。待茄子中心温度降至 -20（即低于共晶温度 8） ，整个冻结过程所用时间为 216 min。对采集数据进行处理及曲线拟合，绘制茄子的 冷冻曲线如图 4 所示。由拟合曲线解析表达式及其曲率 计算公式得 B、C 2 点为拟合曲线上曲率最大的点，曲线 段 AB、BC 和 CD 分别对应冻结过程的 3 个阶段，AB 段 表示物料表面温度降至结晶温度，为降温过程；BC 段表 示物料中水分结晶过程；CD 段表示物料再降温过程。这 3 个阶段所用时间分别为 56、115 和 45 min，冻结过程总 时间为 216 min。而数值计算得到冻结过程 3 个阶段所需 时间分别为53、 111和42 min， 冻结过程总时间为206min。 可见，模型计算结果与实验结果接近，并且由图 4 的冻 结曲线可看出，冻结过程各阶段的冻结规律也基本吻合。 注：A起点；B终点；C、D拟合曲线上曲率最大的点 图 4茄子的冻结曲线 Fig.4Curves of eggplant freezing 6结论与讨论 对球形果蔬物料冻结过程的数值模拟和实验研究表 明： 1）所建立的冻结过程的非稳态数学模型能较好地描 述物料冻结过程温度与冻结界面的变化规律，物料冻结 时间的计算结果与试验结果较吻合。 2）该模型对球形果蔬物料冻结过程的控制与冻结时 间的预测均具有重要指导意义。 由于实验设备所限，冻结界面随时间的变化曲线未 获得实验验证。从计算所得冻结界面的移动曲线来看， 其变化趋势与实际情况是一致的。 参考文献 1赵鹤皋， 郑效东， 黄良瑾， 等 冷冻干燥技术与设备M 武 汉：华中科技大学出版社，2004 2华泽钊，刘宝林，左建国食品和药品的冷冻干燥M北 京：科学出版社，2006 第 5 期郝新生等：茄子冷冻干燥冻结过程模型的建立与求解341 3Lichtfield R J, Liapis A I. An absorption- sublimation model for a freeze dryerJ. Chem Eng Sci, 1979, 34(9): 1085 1090. 4Edinara A B, Filho R M. Freeze drying process: real time modelandoptimizationJ.ChemicalEngineeringand Processing, 2004, 43: 14751485. 5Nastaj J F, Ambrozek B. Modeling of vacuum desorption of mutilcomponent moisture in freeze dryingJ. Transp Porous Med, 2007, 66(10): 201218. 6刘永忠，陈三强，孙皓冻干物料孔隙特性表征的分形模 型与分形维数J农业工程学报，2004，20(6)：4145 Liu Yongzhong, Chen Sanqiang, Sun Hao. Characterizing pores in freeze- dried materials by fractal models and fractal dimensionsJ. Transactions of the CSAE, 2004, 20(6): 41 45. (in Chinese with English abstract) 7郭玉明，姚智华，崔清亮，等真空冷冻干燥过程参数对 升华干燥能耗影响的组合实验研究J农业工程学报， 2004，20 (4)：180184 Guo Yuming, Yao Zhihua, Cui Qingliang, et al. Combined experimentalstudyontheeffectsoftheoperational parameters on energy consumption of sublimation- drying during vacuum freeze- dryingJ. Transactions of the CSAE, 2004, 20(4): 180184. (in Chinese with English abstract) 8罗瑞明，周光宏，乔晓玲干切牛肉冷冻干燥中高速率升 华条件的动态研究J农业工程学报，2008，24(2)： 226231 Luo Ruiming, Zhou Guanghong, Qiao Xiaoling. Dynamic study on operation conditions for the higher sublimation rate during freezing drying of cooked beef sliceJ. Transactions of the CSAE, 2008, 24(2): 226231. (in Chinese with English abstract) 9王明新 数学物理方程M 北京： 清华大学出版社， 2006 10 Crank J. Free and Moving Boundary ProblemsM. Oxford: Clarendon Press, 1984. 11 Natale M F, Tarzia D A. Explicit solutions to the two phase stefan problem for storm type materialsJ. J Phys, 2000, 33: 395404. 12 闫德宝一类新型 Stefan 问题局部解的存在唯一性J山 东理工大学学报：自然科学版，2007，21(1)：6164 Yan Debao. The existence and uniqueness of global solution to a new type of stefan problemJ. Journal of Shandong University of Technology: Science and Technology, 2007, 21(1): 6164. (in Chinese with English abstract) 13 申建中，易法槐一个反应- 扩散方程的自由边界问题 J数学物理学报，2003，23(2)：183192 ShenJianzhong,YiFahuai.AReaction- diffusions free- boundary problemJ. Acta Mathematiea Scientia, 2003, 23(2): 183192. (in Chinese with English abstract) 14 Kar A, Mazumder J. Analytic solution of the stefan problem in finite mediumJ. Quarterly of Applied Mathematics, 1994,