房水是充满眼球前房和后房,夹在角膜和晶体之间的透明液体,由睫状体的无色素上皮细胞(睫状突)分泌,流经后房经瞳孔入前房,然后主要经前房角进入小梁网、Schlemm氏管、集合管以及进入巩膜内和巩膜上静脉.房水在压力差和温度差共同趋动下流动,涉及Poiseuille流动和自然对流现象.认识房水流动的详细流场对于眼睛疾病(如青光眼、白内障、Krukenberg条纹)的产生机理及治疗都有重要指导作用.
对房水流动的解析解研究,需要把眼内结构理想化为规则的几何形状.Friedland[1]通过把虹膜和晶体表面简化为球面,求解不可压缩流动流函数的双调和函数,并使用集总参数模型讨论眼压与小梁网、前房、后房以及睫状体等效阻力的关系.El-Shahed和Abd Elmaboud[2]简化了流体动量方程和能量方程,得到了温差驱动下的前房房水流场.Canning等[3]对前房房水在温差驱动下的流动进行了理论求解,并建立了颗粒在房水中的浓度微分方程模型,定量研究了前房积血、角膜内皮沉着物、前房积脓和Krukenberg条纹等现象.Fitt和Gonzalez[4]通过简化人眼空间形状,解析讨论了前房温差产生的浮力驱动流、房水分泌流、温差和重力共同作用下的流动,晶体震颤产生的房水运动以及睡眠时眼球快速运动时的房水流动.
随着计算流体动力学(CFD)的发展,数值求解眼内复杂流动已成为可能.Heys等[5]建立了人眼二维弹性模型,讨论了正常眼球的晶体曲率变化对前后房压差的影响,分析了虹膜和晶体的最小距离以及小梁网的渗透率等因素对眼压的影响.随后他们数值研究了虹膜变形后,前房房水堆积会导致前房压力比后房高1 mmHg[6](1 mmHg=133 Pa).Kumar等[7]数值模拟了兔眼前房内的房水流动,发现流场中的眼压差主要集中于小梁网处,瞳孔大小改变对于眼压影响不大.Villamarin等[8]建立了三维眼内流动数值模型,并研究了眼压与小梁网渗透率的关系.Crowder和Ervin[9]利用流体力学中的Navier-Stokes方程和渗流理论的Darcy方程对眼内流动进行建模,求得房水由后房晶体与虹膜通道到前房再通过小梁网的流动.郭竞敏等[10]基于房水生理学相关理论及流体力学基本原理,对不同宽窄的前房进行三维重建,并对不同瞳孔直径及体位的房水流动进行了数值模拟.陈伟等[11]根据房水生理学基本理论,在超声生物显微镜扫描结果的基础上构建了以人眼前房为主的二维模型,计算结果表明虹膜膨隆较为严重时,前后房的眼压差值急剧上升.梁书秀等[12]利用ANSYS FLUENT软件计算了虹膜至人工晶体不同间距下的白内障术后眼内前后房压力差,通过和术前对比,分析了白内障术后并发症的水动力因素,为后续白内障术后并发症的研究提供了新的途径.Wang等[13]利用有限元法数值研究了3D打印的人眼模型的房水流场以及虹膜变形,并用粒子图像测速(PIV)技术进行了实验验证,表明重力方向和温度对房水流动影响很大,从而影响虹膜的变形.
由上述文献可以看出,讨论前房流动的文献较多,这是因为前房容积大、角膜与虹膜的温差较大,导致自然对流主要发生在前房.对于前房后房流动的综合研究,由于虹膜与晶体之间距离很小(最小距离在10 μm左右[14]),此处的计算模型处理十分重要.许多学者的计算模型中虹膜-晶体间隙仅使用两层网格,这不足以充分表示虹膜-晶体间隙的流动.
本文对人眼房水流动进行CFD建模,研究不同体位考虑自然对流因素下的房水流动特性,计算两种不同虹膜-晶体间隙(30 μm和5 μm)情况,该间隙处使用11层计算网格,以捕捉该区域的流场物理量梯度,使得计算结果更为准确.近年来流动与传热理论在生物医学工程应用越来越多[15-16],本文的定量研究有助于深入认识眼内房水流动特性,为房水流动相关的眼科疾病如闭角型青光眼的诊断及手术提供参考.
区别于血液,房水为Newton流体,其流动受控于Navier-Stokes方程.同时由于房水流速很小(最大流速在mm/s量级),所以可视为不可压缩流动,其定常流动的连续方程及动量方程为
·v=0,
(1)
ρv·v=-p+μ2v+ρg,
(2)
式中v为速度,ρ,p,μ为房水密度、压力及动力黏度,g为重力加速度矢量.
由于眼内温度梯度的存在, 房水密度发生变化, 产生浮力, 形成自然对流.密度变化可由Boussinesq近似[17]得到
ρ=ρ0β(T-Tref),
(3)
式中ρ0为参考密度,β为房水体积膨胀系数,T和Tref分别为温度和参考温度.
对于不可压缩流动,由于房水流速低忽略黏性引起的耗散项后,描述稳态温度场的方程为
ρcpv·T=k2T,
(4)
其中cp为房水定压比热容,k为房水热传导系数.
本文使用的房水特性见表1.
表1 房水特性[14]
Table1 Properties of AH used in simulations
parametervaluedensity ρ/(kg/m3)990viscosity μ/(Pa·s)7.1×10-4heat conductivity k/(W/(m·K))0.58heat capacity cp/(J/(kg·K))4 200volume expansion coefcient β/K-13.2×10-4
判断房水流动的流态(层态或湍流),需使用特征Reynolds数和Rayleigh数:
(5)
其中D代表特征尺寸,如前房径向尺寸(10-2 m量级),v为特征速度,最大在10-3 m/s量级,本文研究流场的温度ΔT=3 ℃.可以估算出Re≈20≪2 300,Ra≈50≪109,所以可以判断房水流动为层流流动.
前房出口为小梁网,本文建模中使用了0.1 mm厚的各向同性多孔介质区域代表一部分小梁网,并把其内的压力梯度表示为黏性损失项(Darcy定理)和惯性损失项:
-
(6)
式中α为渗透率,C2为惯性阻力系数.它们可由下式估计[7]:
(7)
其中Dp是多孔介质填料床的平均粒径,它与孔径d之间的关系可以由(ε/(1-ε))Dp=d表示;ε为空隙率(定义为空隙体积除以填充床体积).
本小节主要介绍房水流道几何建模及网格划分.
本文的房水流道几何模型主要参考了文献[12,14]中的模型.取眼角膜的曲率半径为8.3 mm,晶体曲率半径为10 mm,在眼中心轴上晶体到眼角膜的距离为3 mm,瞳孔直径为3 mm.晶体与虹膜的缝隙分别设为30 μm和5 μm两个值,文献[14]认为10 μm基本上为该间隙的最小值,本文设置5 μm主要是为了研究极限情况.图1显示了5 μm的人眼房水模型.
图1 房水计算域及边界 图2 计算网格
Fig. 1 The intraocular flow domain and its boundaries Fig. 2 The computational mesh
在专业的网格生成软件ICEM-CFD中使用结构化六面体网格进行区域离散.六面体网格的优点是正交性好,且同样大小区域使用网格量少.为了更好捕捉房水在虹膜-晶体间隙处的详细流场,该处使用11层网格.眼前后房截面网格分布及虹膜-晶体间隙处的网格如图2所示,共计使用约200万个单元.文献[7]研究表明:30万个六面体网格已经可以达到房水流动仿真的网格无关性要求.
图1和图2中出口部分设置厚度为0.1 mm的多孔介质区域,代表小梁网的一小部分.完整小梁网分为三个特征性区域:与前房相接的为葡萄膜网,然后靠外侧是角巩膜网,占小梁网大部分,最后是紧连Schlemm管的邻管区薄层组织[18].前两个区域小梁网组织孔径较大,约为100 μm,而邻管区组织孔径为0.6 μm,小梁网空隙率ε≈0.5.本文多孔介质区域的参数取孔径d=1 μm,空隙率ε≈0.5,作为小梁网一小部分区域的近似.该区域网格层数为11层.设置该多孔计算区域,有助于解决在房角出口处使用压力出口边界引起的回流现象,提高计算精度.
房水流动进口设为速度进口条件,进口速度由房水分泌量2.4 μm/min除以进口面积得到,进口处房水温度设为37 ℃.出口设置为压力出口,设置压力为15 mmHg,即1 995 Pa.壁面设置为无滑移壁面,角膜壁面温度设为34 ℃,其余壁面温度设为37 ℃.
在通用CFD软件ANSYS FLUENT中,分析不同体位的流动,包括平视、仰视、俯视.采用有限体积法对房水流动的偏微分方程进行离散:动量方程和能量方程采用二阶迎风格式离散,压力速度耦合迭代采用SIMPLE算法[19].
迭代计算时,监测房水的体积平均压力及速度.迭代3 000步后,平均压力和速度均趋于常数,表明迭代趋于收敛,再继续迭代5 000步结束.
对虹膜-晶体间隙为5 μm和30 μm两种情况,求得仰视、平视、俯视三种体位的房水流动结果,对不考虑浮力和重力作用的情况也进行了相应计算.
图3显示了前后房内眼压的平均值.可以看出对于四种情况,压力基本相同,5 μm虹膜-晶体间隙眼压比30 μm高5 Pa~8 Pa,这是因为房水流过5 μm虹膜-晶体间隙需要克服更大的阻力.
图4显示了速度的平均值,可以看出整体上房水流动速度很小,在10-5 m/s~10-4 m/s量级,不同体位的房水平均速度可以相差一个数量级.平视时平均速度最大,由下文速度场分布图可知,此时前房存在从高温虹膜侧到低温角膜侧的整体自然对流.不考虑浮力时房水的平均速度值非常小,此时的房水流动只归因于压力趋动的Poiseuille流动,所以可以推断出自然对流将是前房房水的主要流动.
图3 不同体位的眼压平均值 图4 不同体位的房水平均速度
Fig. 3 Mean intraocular pressures of different eye orientations Fig. 4 Mean velocity values of different eye orientations
计算结果发现不同体位的压力场几乎相同,所以此处仅给出平视情况的压力分布,如图5所示,图中为眼球中截面.5 μm虹膜-晶体间隙处的压力变化可参见图5(b),可以看出后房压力比前房约高30 Pa.虹膜-晶体间隙为30 μm时,前后房压力几乎相同,压降集中在前房出口的小梁网处.
(a) 5 μm虹膜-晶体间隙情况 (b) 5 μm虹膜-晶体间隙处
(a) An eye with an iris-lens gap of 5 μm (b) At the gap between iris and lens
(c) 出口处 (d) 30 μm虹膜-晶体间隙情况
(c) At the exit of the anterior chamber (d) An eye with an iris-lens gap of 30 μm
图5 房水压力场
Fig. 5 Pressure distributions
注 为了解释图中的颜色,读者可以参考本文的电子网页版本,后同.
值得指出,虹膜-晶体间隙5 μm比文献[14]认为的最小值10 μm还小,用以研究虹膜间隙对前后房压差的影响.从计算结果可以看出,正常虹膜-晶体间隙情况下前后房的压差不会过大,除非虹膜发生病变使得其与晶体粘连,导致瞳孔阻滞,阻止后房房水进入前房.这也能从流体力学的Bernoulli方程来解释:正常虹膜间隙处房水流速低动压小,所以前后房压差不会过大.当虹膜贴近晶体前表面,后房压力增加,会趋使虹膜离开晶体:对于5 μm和30 μm虹膜-晶体间隙,计算得到虹膜表面x方向的力分别为3.72×10-3 N和5.473 16×10-5 N,可见5 μm时的排斥力比30 μm时高出两个数量级.因此,可以推断出后房压力的增高原因:当虹膜由于病理原因粘贴在晶体前表面时,导致后房房水不能进入前房,而房水仍然从睫状突分泌进入后房,导致后房压力骤升.
从压力分布图上还可以看出小梁网处压力降低明显,接近30 Pa.本文只设置了0.1 mm的厚度,如果设置真实的小梁网厚度,可以预见其内压降大为增加.因此,如果小梁网以及Schlemm管阻力系数提高(如堵塞),将会导致前房房水压力迅速提高.
为了研究房水压力场与速度场的联系及影响,图6显示了前房内微小压力场变化(压差变化在mPa量级).平视时前房对流速度场比其他体位强烈,因此图6(a)中的标尺设置与其余三图不同.从图中可以看出,对于平视前房从下房角到上房角压力逐渐增加,对于仰视从虹膜瞳孔侧向角膜侧依次增加,对于俯视瞳孔侧的压力稍高,对于不考虑自然对流的情况则前房各点压力几乎一样.对照下文的5 μm虹膜间隙不同体位的房水速度场,将发现压力高的地方为温度较高房水汇集之处.
(a) 平视 (b) 仰视
(a) The vertical orientation (b) The upward-facing position
(c) 俯视 (d) 不考虑自然对流
(c) The downward-facing position (d) No buoyancy
图6 5 μm虹膜间隙不同体位的房水详细压力场
Fig. 6 Detailed pressure distributions of aqueous humor (AH) flow with 5 μm gap
两个不同虹膜-晶体间隙情况的温度场类似,以5 μm虹膜-晶体间隙为例,图7显示了不同体位房水流动的温度场.可以看出对于平视,虹膜上半部分的高温区比下半部分的宽,这是典型的竖直平板自然对流热边界层特点,而角膜侧的低温区在下半部分较大;对于仰视,瞳孔附近的高温区伸向虹膜侧,产生羽流状的热特性,即房水受热沿眼球中心轴向虹膜侧流动(见下文5 μm虹膜间隙不同体位的房水速度场);对于俯视和不考虑自然对流,这两种情况的温度场非常相似,俯视情况角膜侧的低温区稍宽,对照下文的速度场可知,这是因为自然对流导致房水在角膜侧从房角向角膜中心流动,经较低温角膜壁的冷却,汇集在角膜中心部位.
(a) 平视 (b) 仰视
(a) The vertical orientation (b) The upward-facing position
(c) 俯视 (d) 不考虑自然对流
(c) The downward-facing position (d) No buoyancy
图7 5 μm虹膜间隙不同体位的温度场
Fig. 7 Temperature contours with 5 μm gap
后房各壁面的温度均为37 ℃,所以温差引起的自然对流现象几乎可以忽略,其流动主要为压差驱动的Poiseuille流动,因此不同体位的后房流动情况基本相同:房水由睫状突分泌进入后房,速度低,再流经虹膜-晶体间隙,此时通道横截面积减小流速提高.图8为虹膜-晶体间隙处的速度场,可以看出类似管道内流的黏性流动模式:壁面附近由于黏性作用流速低,中间流速最大,在10-3 m/s量级;当间隙增大到30 μm,最高流速下降明显,在10-4 m/s量级.
(a) 5 μm间隙 (b) 30 μm间隙
(a) A 5 μm gap (b) A 30 μm gap
图8 虹膜-晶体间隙处的房水速度场
Fig. 8 Streamlines and contours of the velocity magnitude at the iris-lens gap
两个不同虹膜-晶体间隙速度场情况类似,主要流动为前房内的自然对流,不同虹膜间隙时仅在间隙处差别较大.以5 μm虹膜-晶体间隙为例,说明不同体位的房水流动特性.图9显示了5 μm间隙情况的平视、 仰视、 俯视以及不考虑自然对流情况下的速度大小及流线分布图.虹膜-晶体间隙处的速度较大, 最大值在10-3 m/s量级.为了清楚起见,速度大小云图中标尺设为不同,并且速度最大值比标尺上限值大.
对于前房,不同体位速度场相差很大.对于平视(重力加速度(-y)方向,图9(a)),前房流动为整体性的旋涡流动:虹膜侧房水受热膨胀,密度减小,受浮力作用向上运动,流体质点在房角处遇到阻力大的小梁网后转到角膜侧,受到角膜低温壁的影响冷却,密度增加,向下沉落.
对于仰视(重力加速度(-x)方向,图9(b)),温度高的虹膜在下方,瞳孔处的房水向上运动,形成羽流状流动,遇到较冷的角膜壁面后下沉,但中间瞳孔处上来的房水浮力作用占主导地位,所以迫使下沉的房水向四周分散,有一部分房水从小梁网出口处流出,大部分沿虹膜回到瞳孔,形成图9(b)中的截面图中的一对反向旋涡.空间上讲,房水形成沿前房中心轴上升,然后沿角膜壁向四周运动,再沿虹膜回到瞳孔的空间轴对称对流模式.
对于俯视(重力加速度(x)方向,图9(c)),温度高的虹膜在上方,这时房水受热膨胀产生的浮力迫使房水沿虹膜向外侧流动,有一部分房水从小梁网出口处流出,另外一部分转到角膜侧,导致房水从角膜四周汇聚到中心,然后沿中心轴克服温度逆压梯度产生的反向浮力而重新回到瞳孔附近.显然,对于角膜温度比虹膜温度低时,就房水对流运动来说,仰视比俯视更为顺畅,从图中可以看出仰视时的整体速度比俯视时大.这里应该指出,如果角膜温度大于37 ℃(人体暴露在高温环境),那么仰视和俯视的流动情况将反过来.
(a) 平视 (b) 仰视
(a) The vertical orientation (b) The upward-facing position
(c) 俯视 (d) 不考虑自然对流
(c) The downward-facing position (d) No buoyancy
图9 5 μm虹膜间隙不同体位的房水速度场
Fig. 9 Streamlines and contours of the velocity magnitude with 5 μm gap
不考虑浮力时的速度场见图9(d),此时除了虹膜-晶体间隙处的流速大外,其余地方速度都很小.为了更清楚地显示流速场,图中的标尺数值设置在10-6 m/s~10-5 m/s量级,可以看出房水从瞳孔流出后弥散到前房中,其流线较不规则.
综上可知,对于前房流动主要为温差驱动的自然对流,温度较高壁面附近的流体受热后沿重力相反方向上升,受热房水汇聚地方压力稍有升高(变化在mPa量级).
本文利用计算流体动力学对5 μm,30 μm两种虹膜-晶体间隙的不同体位人眼房水流动进行数值研究,得到的结论如下:
1) 对于房水压力,整体上不同体位的平均压力基本相同.5 μm虹膜-晶体间隙的平均值比30 μm的平均值高出5 Pa~8 Pa,这归因于虹膜-晶体间隙处的压力损失.30 μm虹膜-晶体间隙时前后房压力几乎相同,而5 μm时后房比前房高出30 Pa(0.226 mmHg)左右,前后房压差驱动房水从后房流经虹膜间隙到前房.据此可以推断除非虹膜粘贴到晶体上(瞳孔阻滞)造成后房房水不能进入到前房,否则后房房水压力可以认为和前房压力基本一致.另外,小梁网处的压力下降明显,所以该处如果发生堵塞将导致房水压力增加.
2) 与压力不同,不同体位的房水平均速度值相差很大,这主要归因于虹膜与角膜温差所引起的自然对流.平视时平均流速最大,在10-4 m/s量级;其次是仰视和俯视,平均流速在10-5 m/s量级.如果不考虑自然对流,房水平均流速在10-5 m/s量级.虹膜-晶体间隙通道为类似于管道内流的黏性流动模式,5 μm间隙时中间最大流速在10-3 m/s量级,30 μm时为10-4 m/s量级.
3) 后房房水流动主要为压差驱动的Poiseuille流动,而前房房水流动主要为虹膜与角膜温差引起的自然对流.对于平视,前房流动为虹膜侧温度高的房水上升、在角膜侧遇冷下沉的整体循环流动.对于仰视,瞳孔处温度高的房水沿前房中心轴上浮,顺着角膜向四周下沉,然后从虹膜周边回到瞳孔,形成轴对称的对流模式.对于俯视,其对流循环模式与仰视正好相反,前房中心轴上房水从角膜运动到瞳孔处,由于要克服浮力,整体上流速比仰视小.前房自然对流引起的压差变化在mPa量级,温度较高房水汇集区域压力稍高.
我们后续将对房水-虹膜之间的流固耦合、虹膜膨隆时的流场特性进行研究.
[1] FRIEDLAND A B. A hydrodynamic model of aqueous flow in the posterior chamber of the eye[J]. Bulletin of Mathematical Biology, 1978, 40(2): 223-235.
[2] EL-SHAHED M, ABD ELMABOUD Y. On the fluid flow in the anterior chamber of a human eye with slip velocity[J]. International Communications in Heat and Mass Transfer, 2005, 32(8): 1104-1110.
[3] CANNING C R, GREANEY M J, DEWYNNE J N, et al. Fluid flow in the anterior chamber of a human eye[J]. Mathematical Medicine and Biology: a Journal of the IMA, 2002, 19(1): 31-60.
[4] FITT A D, GONZALEZ G. Fluid mechanics of the human eye: aqueous humour flow in the anterior chamber[J]. Bulletin of Mathematical Biology, 2006, 68(1): 53-71.
[5] HEYS J J, BAROCAS V H, TARAVELLA M J. Modeling passive mechanical interaction between aqueous humor and iris[J]. Journal of Biomechanical Engineering, 2001, 123(6): 540-547.
[6] HEYS J J, BAROCAS V H. Computational evaluation of the role of accommodation in pigmentary glaucoma[J]. Investigative Ophthalmology & Visual Science, 2002, 43(3): 700-708.
[7] KUMAR S, ACHARYA S, BEUERMAN R, et al. Numerical solution of ocular fluid dynamics in a rabbit eye: parametric effects[J]. Annals of Biomedical Engineering, 2006, 34(3): 530-544.
[8] VILLAMARIN A, ROY S, HASBALLA R, et al. 3D simulation of the aqueous flow in the human eye[J]. Medical Engineering & Physics, 2012, 34(10): 1462-1470.
[9] CROWDER T R, ERVIN V J. Numerical simulations of fluid pressure in the human eye[J]. Applied Mathematics and Computation, 2013, 219(24): 11119-11133.
[10] 郭竞敏, 张虹, 王军明. 人眼前房三维重建与房水流场数值模拟[J]. 眼科新进展, 2015, 35(4): 346-350.(GUO Jingmin, ZHANG Hong, WANG Junming. 3D model reconstruction and numerical simulation of fluid dynamics in anterior chamber[J]. Recent Advances in Ophthalmology, 2015, 35(4): 346-350.(in Chinese))
[11] 陈伟, 张向东, 余涵, 等. 虹膜膨隆对房水流动影响的数值模拟分析[J]. 眼科新进展, 2017, 37(1): 72-76.(CHEN Wei, ZHANG Xiangdong, YU Han, et al. Numerical simulation of aqueous humor hydrodynamics in human eye with iris bombe[J]. Recent Advances in Ophthalmology, 2017, 37(1): 72-76.(in Chinese))
[12] 梁书秀, 郝菲菲, 孙昭晨, 等. 基于Fluent的白内障术后房水流体力学数值分析[J]. 暨南大学学报(自然科学与医学版), 2014, 35(4): 350-356.(LIANG Shuxiu, HAO Feifei, SUN Zhaochen, et al. Numerical analysis of aqueous humor hydromechanics after cataract surgery based on Fluent[J]. Journal of Jinan University(Natural Science & Medicine Edition), 2014, 35(4): 350-356.(in Chinese))
[13] WANG W, QIAN X, SONG H, et al. Fluid and structure coupling analysis of the interaction between aqueous humor and iris[J]. Biomedical Engineering Online, 2016, 15(2): 570-586.
[14] HEYS J J, BAROCAS V H. A Boussinesq model of natural convection in the human eye and the formation of Krukenberg’s spindle[J]. Annals of Biomedical Engineering, 2002, 30(3): 392-401.
[15] 密斯让·J·C, 辛哈·A, 斯特·G·C, 等. 生物磁粘弹性流体的流动: 应用动脉电磁过热评估血液的流动,癌症治疗进程[J]. 应用数学和力学, 2010, 31(11): 1330-1343.(MISRA J C, SINHA A, SHIT G C, et al. Flow of a biomagnetic viscoelastic fluid: application to estimation of blood flow in arteries during electromagnetic hyperthermia,a therapeutic procedure for cancer treatment[J]. Applied Mathematics and Mechanics, 2010, 31(11): 1330-1343.(in Chinese))
[16] 李方方, 刘静, 乐恺. 生物组织在冻结过程中的三维相变传热问题精确解[J]. 应用数学和力学, 2009, 30(1): 64-74.(LI Fangfang, LIU Jing, YUE Kai. Exact analytical solution to three dimensional phase change heat transfer problems in biological tissues subject to freezing[J]. Applied Mathematics and Mechanics, 2009, 30(1): 64-74.(in Chinese))
[17] BERGMAN T L, LAVINE A S, INCROPERA F P, et al. Fundamentals of Heat and Mass Transfer[M]. Hoboken: John Wiley & Sons, Inc, 2017.
[18] AVTAR R, SRIVASTAVA R. Modelling aqueous humor outflow through trabecular meshwork[J]. Applied Mathematics and Computation, 2007, 189(1): 734-745.
[19] 陶文铨. 数值传热学[M]. 西安: 西安交通大学出版社, 2001.(TAO Wenquan. Numerical Heat Transfer[M]. Xi’an: Xi’an Jiaotong University Press, 2001.(in Chinese))
蔡建程(1980—),男,副教授,博士(E-mail: cai_jiancheng@foxmail.com);
曹月红(1979—),女,硕士(通讯作者. E-mail: 25258581@qq.com).