KDP晶体生长的溶液流动和物质输运计算

KDP晶体生长的溶液流动和物质输运计算 第 39卷 第 1期 人工晶体学报 Vo .l 39 No. 1 2010年 2 月 JOURNAL O F SYN TH ET IC CR YSTAL S Feb rua ry, 2010 KD P晶体生长的溶液流动和物质输运计算 王晓丁 , 李明伟 , 曹亚超 , 刘玉姗 () 重庆大学动力工程学院 ,重庆 400030 摘要 :本文采用有限容积法 ,对 KD P晶体生长过程中溶液的流动和物质输运进行了数值模拟 。结果表明 :随着入口 溶液流动速度的增大 ,籽晶的上表面因自然对流而引起的抽吸作用减小 ,表面过饱和度的最小值沿 x 正向发生右 移 ,其上表面的剪切力先减小后增大 。随着入口溶液过饱和度的增大 ,籽晶上表面剪切力增大 。不同尺寸的籽晶 表面过饱和度的分布差异较大 。籽晶的生长边界层厚度与溶液流动密切相关 ,入口溶液流动速度越大 ,厚度越小 , 但其受入口溶液过饱和度的影响较小 。 关键词 : KD P晶体 ;数值模拟 ;表面过饱和度 ;剪切力 ;生长边界层 ( ) 中图分类号 : O614 文献标识码 : A 文章编号 : 10002985X 2010 0120088 207 Num er ica l S im u la t ion of the F low an d M a ss Tran sfer in the Grow th of KD P C ry sta ls from So lu t ion W AN G X iao2d ing, L I M ing 2w ei, CAO Ya 2chao, L I U Yu 2shan ( )Co llege of Power Enginee ring, Chongq ing U n ive rsity, Chongq ing 400030 , Ch ina ( )R eceived 16 M a rch 2009 A b stra c t: N um e rica l sim u la tion s of the flow and m a ss tran sfe r invo lved in the grow th of KD P c rysta ls had been p e rfo rm ed by the fin ite vo lum e m e thod. The re su lts showed tha t the p ump ing ac tion a t the upp e r su rface of the seed c rysta l cau se by the na tu ra l convec tion dec rea se s, the m in im um of the su rface sup e rsa tu ra tion ha s a righ t sh ift a long the po sitive x 2axis and the shea r stre ss a t the su rface firstly inc rea se s and then dec rea se s a s the in le t ve loc ity inc rea sing. The shea r stre ss inc rea se s w ith the inc rea se of the in le t sup e rsa tu ra tion. The d istribu tion of the su rface sup e rsa tu ra tion ha s an obviou s d iffe rence when changing the geom e trie s of the seed c rysta l. The th ickne ss of the grow th bounda ry laye r is a ssoc ia ted c lo se ly w ith the flu id flow , bu t it is influenced sligh tly w ith the in le t sup e rsa tu ra tion. Key word s: KDP crysta l; nume rical simula tion; surface supersaturation; shear stre ss; grow th boundary layer 1 引言 )( KD P KHPO晶体是 20世纪 30、40年代发展起来的一种综合性能非常优良的非线性光学材料 。由于 2 4 它具有较大的非线性光学系数和高激光损伤阈值的特性而被广泛应用于激光变频 、电光调制和光快速开关 等高科技领域 。随着惯性约束核聚变的发展 ,工业上对 KD P晶体的尺寸和质量提出了更高的要求 。为了更 收稿日期 : 2009 203 216 () 基金项目 :国家自然科学基金 No. 50676113 , 50976127 ( ) 作者简介 :王晓丁 1983 2,男 ,四川省人 ,博士 。 E2m a il: wxd_03 @ 163. com 通讯作者 :李明伟 ,教授 。 E2m a il: aoweixia@ 126. com 好地理解 KD P晶体在溶液中的结晶过程 ,除了传统的理论分析和实验研究外 ,数值模拟也成为研究晶体生 长必不可少的方法之一 。 [ 1 ][ 2 ]Zhou和 Yecke l等 完成了对加速 /减速旋转的 KD P晶体生长的三维时相关模拟 ,发现晶体生长受流 动的影响很大 ,但是由于并未考虑物质输运过程 ,因此很难将晶体的实际生长与计算联系起来 。B ra ilovskaya [ 3 ]等 考虑了自然对流和强制对流的影响 ,对不同入口溶液流动速度下的过饱和溶液流过 KD P晶体表面的情 况进行了二维非稳态数值模拟 ,得到了晶体的表面浓度和流函数的分布 ,对认识浓度边界层的均化作用和提 [ 4 ][ 5 ] 高晶体生长的质量有很好的指导作用 。 Robey和 M ayne s采用 Bou ssine sq假设模拟了三维条件下快速生 长大型 KD P单晶时相关的湍流流动情况 ,得到了不同尺寸 、不同生长速率和不同旋转条件下晶体表面剪切 力和过饱和度随空间和时间的变化规律 ,阐述了包裹物的形成与晶体表面剪切力和过饱和度之间的关系 。 [ 6 ]L iiri等 引入生长速率方程和表面反应方程 ,通过计算 KD P 晶体周围的速度和溶质边界层厚度 ,最终建立 了一个实际生长的自洽模型 。该模型对于研究晶体局部生长速率和平均生长速率受溶液浓度 、速度和晶体 方位的影响是有很大帮助的 。 [ 4 , 5 ] 晶体表面过饱和度和剪切力对晶体表面形貌稳定性以及包裹物的形成是有影响的 ,但究竟是如何 影响的 ,到目前为止 ,还没有较好的阐述 。该问题的难度主要是对表面过饱和度和剪切力的测量是难以实现 的 ,依靠数值计算可以得到上述值 ,但又该如何与表面形貌的实验结果相联系 ? 实验中提出了一种新的思 路 ,利用 A FM 实时原位扫描 ,获得表面形貌 ,同时进行与实验条件尽可能一致的数值模拟 ,然后对比实验获 得的表面形貌和数值计算得到的表面过饱和度和剪切力值 ,从而了解表面过饱和度和剪切力对晶体表面形 貌的影响情况 。 本文对 KD P晶体生长过程中溶液的流动和物质输运进行了稳态三维数值模拟 ,重点考察了不同入口溶 液流动速度 、不同入口溶液过饱和度和不同籽晶尺寸下籽晶上表面的过饱和度 、剪切力以及生长边界层厚度 的情况 。计算结果与液下 A FM 实时扫描实验相结合 ,将有利于更好地了解晶体生长的微观机制 。 2 物理模型和数学模型 2. 1 物理模型 ( )( ) 模拟采用与液下 A FM 实时扫描实验完全一致的物理模型 ,如图 1 a 、b 所示 ,籽晶放臵于内径为 12 mm ,溶液高度为 3. 5 mm 的圆柱型液池中间 ,坐标原点为液池底面的中心 ,即籽晶柱体下表面迎来流一侧中 ( ) ( ) 点 ,进口处中心点坐标为 - 13. 5 , - 20 , 1 ,出口处中心点坐标为 11. 5 , - 20 , 1 ,进出口管道的管径为 < = 1 mm ,籽晶的柱体尺寸分别为 size= 0. 5 mm ×0. 5 mm ×0. 5 mm、size= 1 mm ×1 mm ×1 mm 和 size= 1. 5 1 2 3 mm ×1. 5 mm ×1. 5 mm ,籽晶锥体高度为柱体边长的 1 /2。 ( )( )( )图 1 a生长液池正视图 ; b生长液池俯视图 ; c籽晶示意图 ( ) ( ) ( ) F ig. 1 aF ron t ou tline of grow th c rysta llize r; bV e rtica l view of grow th c rysta llize r; cGeom e try of the seed c rysta l 为方便作如下规定 :籽晶的上下表面分别称为 up 面和 down面 ,正对出口的面称为 righ t面 , y轴正 、负值 侧的面分别称为 beh ind面和 fron t面 。与柱体的 up 面 、down 面 、fron t面和 beh ind 面相连的锥面分别称为 39卷第 90 人 工 晶 体 学 报 cone2up 面 、cone2down面 、cone2fron t面和 cone2beh ind面 。在 up 面上取与 x 轴平行的中线 ,记为 up line,如图 ( ) 1 c所示 。 2. 2 数学模型 2. 2. 1 基本假设 ( )( ) 为简化可作如下假定 : 1 系统处于稳态的层流 ; 2 KD P溶液为可压缩流体 ,除密度外 ,所有物性参数 ( )( )均为常数 ; 3 溶液上表面是与大气相接的自由表面 ,其余固壁表面满足无滑移条件 ; 4 籽晶表面有生长源 项 ,即质量源项 S和动量源项 S。 1 2 2. 2. 2 控制方程及边界条件 在上述假定下控制方程可写为 : ( )(ρ) 1 A 〃 V = 0 ρρτ ( ) 2 V 〃A V = - A p - A 〃 - gez ( )(ρ) (ρ) 3 A 〃 V c= A 〃 D A c 边界条件为 : 溶液上表面 : ν 9u 9 ( )ω 4 = 0 , = 0 , = 09z 9z 容器各壁面 : νω ( )u = 0 , = 0 , = 0 5 籽晶表面 : νω ( )u = 0 , = 0 , = 0 6 9c (ρρ) ρ( )7 m = - cR =D | s 0 e 0 0 9n ρτz轴的单其中 , V 为速度矢量 , 为溶液密度 , p为压力 ,为粘性应力张量 , g 为重力加速度 , e为平行于 z - 9 2 () νω位矢量 , c是用 KD P质量分数表示的溶液浓度 , D 为物质扩散系数 D = 1 ×10 m/ s, u、、分别为 x 向 、y 3 ρ(ρρ) 向 、z向速度分量 , m 为籽晶表面质量流率 ,为 KD P晶体的密度 = 2338 kg /m , 为 50 ?时 KD P饱和 s s 0 3 (ρ( ) ( ) ( ) ) 溶液的密度 = 1200 kg /m , c为 50 ?时 KD P溶液的平衡浓度 c= 0. 2984 g KD P/ g so lu tion。 0 e e 3 网格划分及数值方法 籽晶锥体附近区域和进出口管道采用混合网格 ,其它区域采用六面体的结构网格 ,其中籽晶柱体近表面 区域进行 局部 加 密 , 最密 处网 格间 距 为 0. 05 mm。选用 size、size和 size时 所 构 造 的 结 点 数 目 分 别 为 1 2 3 355379 , 420942 , 483432。 [ 7 ] 模拟采用有限容积法对基本方程进行二阶迎风离散 ,使用 Simp le 算法 进行求解 。当质量残差达到- 5 10 时 ,认为离散方程的解收敛 。 4 结果与讨论 4. 1 入口溶液流动速度的影响 首先对籽晶尺寸为 size,入口溶液过的饱和度为 4 %时 ,计算了入口溶液的流动速度分别为 0. 008 m / s、 2 0. 011 m / s、0. 014 m / s、0. 017 m / s、0. 020 m / s时溶液的流动和物质输运情况 。 ( )( )图 2 a是入口溶液流动速度 V 为 0. 008 m / s时 , y = 0 截面在籽晶附近的速度矢量图 。从图 2 a中可 in 以看出 ,籽晶 up 面上的自然对流现象非常明显 ,这是由于籽晶表面附近区域存在生长源项 ,在生长溶液中形 ( ) ( )成了浓度差造成的 。在紧靠 up 面的水平截面 z = 1. 05 mm ,取其局部区域的流场 ,如图 2 b 所示 。该流 ( ) 场面积大小为 2 mm ×2 mm ,中心 1 mm ×1 mm 区域正对 up 面 。结合图 3 aup 面的等过饱和度线可以发现 σup 面上自然对流的抽吸作用使流体往过饱和度最低值 附近汇聚 ,在汇聚过程中由于籽晶不断消耗溶 m in 质 ,因而溶液过饱和度逐渐降低 。 ( ) ( )( ) 图 3 a、b、c分别为入口溶液流动速度 V 为 0. 008 m / s、0. 014 m / s、0. 020 m / s时 up 面的等过饱 in σσ和度线 。可以看出 :随着 V 的增大 , up 面的过饱和度 也随之增大 ,且最小值 也沿着 x 正向发生右移 。 inm in )σ(4 当 V = 0. 020 m / s时 ,几乎位于最右端 ;从 y = 0和 z = 1. 05 mm 两截面在籽晶附近的速度矢量图 图 in m in 还可以看出 ,此时位于 up 面上自然对流已很微弱 ,抽吸作用已很不明显 。 ( ) ( ) 图 2 ay = 0截面在籽晶附近的速度矢量图 ; bz = 1. 05 mm 截面在 up 面正上方的速度矢量图 σ( ) ( ) F ig. 2 aV e loc ity vec to r a t the sec tion of y = 0 nea r the seed c rysta l size, = 4% , V = 0. 008 m / s; 2 in in ( ) σ( )bV e loc ity vec to r a t the sec tion of z = 1. 05 mm above the up su rface size, = 4% , V = 0. 008 m / s 2 in in 图 3 up 面等过饱和度线 σσ σσ( ) ( ) F ig. 3 Con tou r line of sup e rsa tu ra tion on up su rface size, = 4% , = 0. 1 % a= 3. 1% , = 1. 6 % ,2 in m ax m in σσσσ( ) ( ) V = 0. 008 m / s; b= 3. 1 % , = 1. 8 % , V = 0. 014 m / s; c= 3. 8% , = 2. 1 % , V = 0. 020 m / sin m ax m in in m ax m in in ( ) ( ) 图 4 ay = 0截面在籽晶附近的速度矢量图 ; bz = 1. 05 mm 截面在 up 面正上方的速度矢量图 σ( ) ( ) F ig. 4 aV e loc ity vec to r a t the sec tion of y = 0 nea r the seed c rysta l size, = 4% , V = 0. 020 m / s; 2 in in ( ) σ( )bV e loc ity vec to r a t the sec tion of z = 1. 05 mm above the up su rface size, = 4 % , V = 0. 020 m / s 2 in in 图 5是不同入口溶液流动速度下 up line上的过饱和度值 。当 V = 0. 008 m / s时 ,在 x = 0. 725 mm 处有 in 极小值 ,说明在该点附近区域的自然对流最剧烈 。当 V = 0. 011 m / s和 0. 014 m / s时 ,极小值在 x = 0. 825in 39卷第 92 人 工 晶 体 学 报 mm 和 0. 925 mm 处 ,也对应着自然对流最剧烈的区域 。而 V 增大为 0. 017 m / s和 0. 020 m / s时 ,并未发现 in 过饱和度的极小值 ,此时 up 面上的抽吸作用和自然对流现象已不明显 。显然 ,低入口流速时 , up 面上因自 然对流的抽吸作用而引起的回流 ,使 up line上过饱和度存在一极小值 ,而随着 V 的增加 ,极小值将消失 。 in 不同 V 时 up line上的剪切力值如图 6 所示 。可以看出 ,在自然对流剧烈的区域 ,剪切力小 。由壁面剪in 9u 9u x yτττμ( )μ切力公式 : =+=+ 可知 ,在 一定时 ,决定 up 面剪切力大小的是 x 和 y 向速度分量在 z方zx zy 9z 9z 向的变化率 ,这一变化率在自然对流强的区域小 ,因此剪切力必然在自然对流最强烈的区域产生最小值 。此 外 ,图 6还表明 ,在各剪切力曲线的起始段 ,剪切力随入口溶液流动速度的增大先减小后增大 ,不具有单调 性 ,这应该是自然对流和强制对流综合作用的结果 。当流速增大到一定程度以后 ,自然对流的作用不明显 了 ,强制对流成为了决定流场的主要因素 ,这时流速越大 ,剪切力也就越大 。 图 5 不同 V 时 up line上的过饱和度值图 6 不同 V 时 up line上的剪切力值 inin F ig. 5 Sup e rsa tu ra tion on up line F ig. 6 Shea r stre ss on up line ( σ( σ))a t d iffe ren t V size, = 4 % a t d iffe ren t V size, = 4 % in 2 in in 2 in 图 7 是 y = 0 截面在籽晶附近区域的等过饱和度 σ线 。该区域有最大值 = 4 % ,与入口溶液过饱和度 m ax 相同 。位于 锥 面 底 部 区 域 的 溶 液 过 饱 和 度 最 小 , 为 σ= 1 % , 这是因为该区域溶液流动最弱 , 溶质难以 m in 补给 。在 up 面上方的等过饱和度线与前面提到的自 然对流的抽吸是对应一致的 。此外 ,从图 7 中不难看 出 ,在紧靠晶 /液界面存在溶质浓度发生剧烈变化的流 δ体层 ,也就是籽 晶的 生 长边 界层 , 其 厚 度 通常 定 义 图 7y = 0截面在籽晶附近的等过饱和度线 为 :生长边界层外缘处流体的浓度和晶面浓度差为最 F ig. 7 Con tou r line of sup e rsa tu ra tion a t the sec tion of y = 0 大浓度差的 99 %时的距离 。( σσnea r the seed c rysta l size, = 4 % , V = 0. 014 m / s, 2 in in m 入口溶液流动速度 V = 0. 014 m / s, 入口溶液过 in σσ )= 4 % , = 1. 0% , = 0. 1 % ax m in σ 饱和度 = 4 %时 , size籽晶 up line 上的生长边界层 in 2 δδ厚度 如图 8所示 ,沿着 x 的正方向 ,生长边界层厚度 是逐渐增大的 。在 x ?0. 7 mm ,增大的趋势并不明 显 ,当 x ?0. 7 mm 时 ,曲线陡增 。这主要是因为在 x = 0. 7 mm 附近有强烈的自然对流 ,从而导致在该点之后 的区域生长边界层厚度急剧增加 。 δ 图 9给出了 up 面上中心点生长边界层厚度 随入口溶液流动速度的变化 。入口溶液流动速度从 0. 008 δm / s增加至 0. 011 m / s,减小得较快 ,而当入口溶液流动速度逐渐从 0. 011 m / s增加到 0. 020 m / s的过程 中 ,生长边界层厚度虽然一直减小 ,但减小幅度较小 。需要指出的是 ,在有些情况下 ,由于 up 面上发生抽吸 作用的位臵比较靠前 ,若按定义计算的边界层厚度应为无穷大 ,在分析的时候 ,要注意区分这种情况 。 δδ9 up 面中心点的生长边界层厚度 随 V 的变化图 图 8 up line上生长边界层厚度 随坐标 x的变化 in F ig. 8 Th ickne ss of grow th bounda ry laye r on up line F ig. 9 Th ickne ss of grow th bounda ry laye r a t the cen te r po in t σ( )x coo rd ina te s size, = 4% a t d iffe ren t ( σ)of up su rface a t d iffe ren t V size, = 4% 2 in in 2 in 4. 2 入口溶液过饱和度的影响 σ图 10 和图 11 分别是入口溶液过饱和度 为 2 % 、3 % 、4 % 、5 % 、6 %时 , up line上的过饱和度和剪切力 。 in σσ从图中发现 ,尽管不同的 造成生长溶液的整体过饱和度不同 ,越大 ,生长溶液整体过饱和度越大 ,但是 inin 过饱和度和剪切力沿 up line的变化趋势是一致的 ,只是曲线极小值的位臵略有差别 。这是因为当 V 一定in时 ,对流物质输运一定 ,计算中假设了法向生长速度与入口溶液过饱和度存在线性关系 ,入口溶液过饱和度 越小 ,那么籽晶表面的法向生长速度就越小 ,溶质的消耗量就越小 ,整个生长液池中的浓度差也越小 ,自然对 σ流现象就越弱 ,因此出现极小值的位臵相对靠后 。当 = 2 %时 ,自然对流很不明显 ,因此在 up line上的过 in σ 饱和度和剪切力曲线均没有出现极小值 。此外 ,随着入口溶液过饱和度 的增大 ,自然对流的强度就越大 ,in () 因此在强制对流不变 入口溶液流动速度不变 的情况下 , up line上的剪切力也就更大一些 。 σσ图 10 不同 下 up line上的过饱和度值图 11 不同 下 up line上的剪切力值 inin F ig. 11 Shea r stre ss on up line F ig. 10 Sup e rsa tu ra tion on up line σ( )σ( ) size, V = 0. 014 m / sa t d iffe ren ta t d iffe ren tsize, V = 0. 014 m / s in 2 in in 2 in σup 面 取 up 面的中心点 ,研究不同入口溶液过饱和度 下 ,该点对应的生长边界层厚度 。图 12给出的 in σσ上中心点生长边界层的厚度随 而变化 。入口溶液流动速度 V 一定时 ,入口溶液过饱和度 对生长边界 ininin 层厚度的影响不大 ,这就说明生长边界层厚度与溶液流动密切相关 ,即物质输运受流动速度的影响很大 。 4. 3 籽晶尺寸的影响 σ( )( ) 图 13 为 = 4 % , V = 0. 014 m / s, size和 size的 up 面上的等过饱和度线 。比较图 13 a 和图 13 b in in 1 3 不难发现 ,尽管入口溶液的条件一样 ,籽晶尺寸不同 ,其 up 面上的过饱和度分布的差异也是较大的 。随籽晶 尺寸增加 ,表面过饱和度略有减小 ,溶质消耗增加 。 39卷第 94 人 工 晶 体 学 报 图 13 up 面等过饱和度线 σδF ig. 13 Con tou r line of sup e rsa tu ra tion on up su rface 图 12 不同 下 up 面中心点的生长边界层厚度 in (σσ ) = 4% , V = 0. 014 m / s, = 0. 1 % F ig. 12 Th ickne ss of grow th bounda ry laye r in in ( ) σσa t the cen te r po in t of up su rface a= 3. 0% , = 1. 4 % , size; m ax m in 1 σ( ) a t d iffe ren tsize, V = 0. 014 m / s( ) σσb= 3. 2 % , = 2. 2% , size in 2 in m ax m in 35 结论 本文对 KD P生长过程中溶液的流动和物质输运情况进行了数值模拟 ,得到以下结论 : ( )σ1 入口溶液过饱和度 不变 ,随着入口溶液流动速度 V 的增大 , up 面的抽吸作用和自然对流现象减 inin 弱 , up line上过饱和度的极小值也沿着 x 正向发生右移 。当 V 增大到 0. 020 m / s,在 up line上无极小值 ; in ( )σ2 在自然对流剧烈的区域 up line上的剪切力小 。入口溶液过饱和度 不变 ,剪切力随入口溶液流动 in速度 V 的增加先减小后增大 ,不具有单调性 ,这是因为自然对流和强制对流综合作用的结果 。V 不变 ,剪切 inin σ力随 的增大而增大 ; in ( )3 不同的籽晶尺寸 up 面上的过饱和度分布的差异较大 。随着籽晶尺寸的增加 , up 面过饱和度略有减 小 ,溶质消耗增加 ; ( )δσδδ 4 生长边界层厚度 与溶液流动密切相关 , V 越大 ,越小 。对 的影响很小 。inin 参 考 文 献 [ 1 ] Zhou Y M , Jeffrey J , D e rby. Th ree2d im en sional Comp u ta tion s of So lu tion H yd rodynam ic s D u ring the Grow th of Po ta ssium D ihyd rogen Pho sp ha te I. Sp in Up and Steady Ro tation [ J ]. J ou rna l of C rysta l G row th, 1997 , 180: 497 2509. [ 2 ] Yeckel A , Zhou Y M , D enn is M. Th ree2d im en siona l Comp u ta tion s of So lu tion H yd rodynam ics du ring the Grow th of Po ta ssium D ihyd rogen Pho sp ha te II. Sp in Down [ J ]. J ou rna l of C rysta l G row th, 1998 , 191: 206 2224. B ra ilovskaya V A , Zilbe rbe rg V V , Feok tistova L V. N um e rica l Inve stigation of N a tu ra l and Fo rced So lu ta l Convection above the Su rface of a [ 3 ] Grow ing C rysta l[ J ]. J ou rna l of C rysta l G row th, 2000 , 210: 767 2771. [ 4 ] Robey H F, M ayne s D. N um e rica l Sim u la tion of the H yd rodynam ic s and M a ss Tran sfe r in the L a rge Scale, R ap id Grow th of KD P C rysta ls. Pa rt1: Comp u ta tion of the Tran sien t, Th ree2d im en siona l F low F ie ld [ J ]. J ou rna l of C rysta l G row th , 2001 , 222: 263 2278. Robey H F. N um erica l Sim u la tion of the H yd rodynam ic s and M a ss Tran sfer in the L a rge Sca le, R ap id Grow th of KD P C rysta ls22: Comp u ta tion of [ 5 ] the M a ss Tran sfer[ J ]. J ou rna l of C rysta l G row th, 2003 , 259: 388 2403. ( ) L iiri M , Enqvist Y, Ka lla s J , e t a l. CFD Mode lling of Single C rysta l Grow th of Po ta ssium D ihyd rogen Pho sp hate KD P from B inary W ate r [ 6 ] So lu tion a t 30 ?[ J ]. J ou rna l of C rysta l G row th , 2006 , 286: 413 2423. () 陶文铨. 数值传热学 第二版 [M ]. 西安 : 西安交通大学出版社 , 2001: 207 2231.[ 7 ] ( ) Tao W Q. N um e rica l H eat Tran sfe r Second ed ition[M ]. X i’an: X i’an J iao tong U n ive rsity P re ss, 2001: 207 2231.