R134a喷雾冷却系统换热性能研究
2015-12-23钱春潮徐洪波邵双全田长青周光辉
钱春潮徐洪波邵双全田长青周光辉
(1中国科学院理化技术研究所热力过程节能技术北京市重点实验室 北京 100190;2中原工学院能源与环境学院 郑州 450007)
R134a喷雾冷却系统换热性能研究
钱春潮1,2徐洪波1邵双全1田长青1周光辉2
(1中国科学院理化技术研究所热力过程节能技术北京市重点实验室 北京 100190;2中原工学院能源与环境学院 郑州 450007)
本文建立了以R134a为冷却工质的封闭式喷雾冷却系统,研究了工质过冷度、质量流量和热流密度对喷雾冷却系统换热性能的影响。其中,工质过冷度由喷嘴入口前的过冷段控制,质量流量通过变频齿轮泵调节,热流密度通过改变加热电源电压和电流控制。实验结果表明,在热流密度和质量流量保持不变时,改变过冷度对热源表面温度和换热系数的影响并不明显;在热流密度和过冷度保持不变的条件下,系统存在一个临界质量流量值,在质量流量达到临界值之前,热源表面温度随质量流量的增大而降低,当质量流量高于临界值时,热源表面温度随质量流量的增大而升高;当质量流量和过冷度保持不变时,存在一个热流密度使液滴的蒸发量等于补充量,在此热流密度下热源表面系数能达到最大。
喷雾冷却;传热系数;过冷度;热流密度;质量流量;R134a
近年来,随着电子器件的小型化、集成化,局部的热流密度越来越高,散热问题日益突出,严重影响电子产品的可靠性以及使用寿命。以高功率固体激光器为例,其在使用过程中只有少部分能量转化为激光,绝大部分能量转化成废热,如果这些废热不能及时散发出去会严重降低激光光束质量,甚至损害激光介质[1-2]。所以解决高功率、高热流密度电子设备的散热问题尤其重要。很多研究表明,喷雾冷却具有较高的换热系数、较低的流量、冷却均匀等优点,因此喷雾冷却技术是解决诸如此类电子设备散热问题的一条有效途径[3-4]。
喷雾冷却换热过程中包括多种换热机制,如对流换热、液膜蒸发换热、核态沸腾和二次成核沸腾换热等,而且影响喷雾冷却性能的因素比较多,如喷嘴高度、喷射角度、雾滴速度、热源表面粗糙度、热流密度等[5-6]。Pautsch A G等[7]研究发现液膜的厚度对冷却性能的影响较大,液膜过厚,阻碍换热。王亚青等[8]发现当喷嘴处于最佳高度时,倾斜角度越大换热效果越佳,冷却效果越好。Rini D P等[9]研究发现一定热流密度下,增大液滴通量,可以降低热源表面温度。Lin L等[10]研究发现当热源表面过热度为定值,临界热流密度随体积通量的增大而增大。Si Chunqiang等[11]研究发现闭式喷雾冷却系统中,随着喷嘴进口压力的增加换热系数增加,热源表面温度降低。热源表面越粗糙换热特性越差,热源表面刻微槽道可以有效提高散热效果[12]。
R134a作为使用最广泛的中低温环保制冷剂,破坏臭氧潜能值为0,不过其全球变暖潜能值略高,但是由于良好的综合性能,它在相当长的时间内仍然将作为一种过渡工质广泛使用。因此,以R134a为冷却工质的喷雾冷却系统在国内外得到了较多的研究。
钱洋等[13]认为喷雾冷却系统以R134a作为冷却工质能够满足高功率电子元件散热需求且具有稳定的散热能力,其研究发现,当流量保持不变时,增大热流密度换热系数先快速升高然后有所下降。Tan Y B等[14]建立了以R134a为冷却工质的喷雾冷却系统,研究发现换热系数高达39000 W/(m2·℃),此时热流密度为145 W/cm2,热源表面温度为48℃。Hou Yu等[15]认为R134a在喷雾冷却系统中具有很好的换热效果,并研究了体积流量对喷雾冷却换热性能影响,同时指出增大工质体积流量可以增大临界热流密度,如体积流量为0.356 L/min时,最大临界热流密度为117.2 W/cm2且热源表面温度为46℃。Li Qiang等[16]研究了R134a喷雾冷却系统中工质体积流量对换热系数的影响,当流量小于1.6 L/min时,换热系数随着流量的增大而增大,当流量超过1.6 L/min后,继续增大流量对换热系数的影响并不明显。Hsieh S S等[17]研究了非沸腾区R134a的液膜厚度,当质量通量在1.33~1.4 kg/(m2·s)之间时,液膜厚度在0.95~1.35 mm之间。Eduardo Martínez⁃Galván等[18]研究了热源表面粗糙度对R134a喷雾冷却性能的影响,并指出液膜厚度与努塞尔特数存在一定的关系。
本文建立了以R134a为冷却工质的闭式喷雾冷却系统,主要研究过冷度、工质质量流量、热流密度影响喷雾冷却性能的规律,以便于喷雾冷却系统选择合适的运行工况。
1 实验系统与数据处理
1.1 实验系统
如图1所示,喷雾冷却系统主要由动力装置、喷雾换热装置、辅助制冷装置以及数据采集系统组成。其中动力装置主要由一齿轮泵提供循环动力,其转速可由变频器调节。喷雾换热装置为系统的关键部件,主要由喷嘴、喷雾室及模拟热源等构成,在喷嘴的进口及喷雾室分别设置热电偶及压力传感器测量相关温度及压力信息。辅助制冷装置包括冷水机组、冷凝段、保证泵进口制冷剂为全液态的过冷段1以及调节喷嘴进口前制冷剂过冷度的过冷段2、储液器、干燥过滤器等。数据采集系统由质量流量计、压力传感器、热电偶、数据采集器和电脑等组成。本实验所采用的喷嘴雾化锥角为60°,口径为0.83 mm,安装高度10.4 mm,雾化颗粒直径在70~90 μm范围内。其他主要器件的详细信息如表1所示。
图1 喷雾冷却系统示意图Fig.1 Schematic diagram of the experimental system
表1 主要实验器件参数Tab.1 The detail information of experimental devices
该系统的循环过程为:工质R134a从储液罐流经过滤器和过冷段1进入齿轮泵,然后工质通过质量流量计和过冷段2进入喷嘴,经喷嘴雾化后成为雾滴,进入喷雾室并喷射到热源表面进行换热,工质冷却热源后的成气液两相态或气态从喷雾室底部流出,经过冷凝段冷凝成液态后流入储液器,进行再次循环。实验过程中喷雾范围能完全覆盖热源表面。
模拟热源结构及其内部热电偶布置,如图2所示,制作热源材料为紫铜,圆柱体底座直径为42 mm,在底座内部均匀排布5根加热棒,底座通过45°锥台渐变成直径为12 mm的圆柱体,直径12 mm的圆柱体上端面为热源表面,该结构具有较好的一维导热性[14]。
图2 热源结构示意图Fig.2 Heat source structure
热源与喷雾室接触处添加聚四氟乙烯隔热垫,热源四周包裹20 mm厚绝热材料防止向环境漏热。为了测得温度参数,沿热源轴线方向,在距离换热表面3 mm、6 mm、9 mm处各布置一个热电偶,所测温度值分别为T1、T2、T3,并根据傅里叶导热定律计算出热流密度值。
1.2 数据处理及误差分析
热源表面温度及换热系数是衡量换热性能的重要指标,尤其是热源表面换热系数,所以本文主要研究工质过冷度,工质质量流量以及热流密度对热源表面温度和传热系数的影响。根据热源结构和隔热处理措施,忽略热源径向的温度变化,热源内部导热看作轴向的一维导热问题。通过热源内部热电偶测得温度数值结合傅里叶定律利用公式(1)计算出热源的热流密度,利用公式(2)计算出热源表面温度。通过压力传感器测出喷雾室内压力,然后得到R134a在此压力下对应的蒸发温度,利用公式(3)计算出换热系数。
式中:q为热流密度,W/cm2;λ为紫铜的导热系数,W/(m·℃);δ1为T1、T2测点之间的距离,mm。
式中:Tsur为热源表面温度,℃;δ0为T1测点和热源表面间的距离,mm。
式中:α为换热系数,W/(m2·℃);Tsat为喷雾室的蒸发温度,℃。
热电偶的测量误差为±0.1℃,相邻热电偶间距离的测量误差为±0.02 mm。通过方程(4)~方程(7)可计算得到热流密度q、热源表面温度Tsur及换热系数α的相对误差值[19]:
式中:φq为q的相对误差,%;ΔT1为T1、T2的测量误差,℃;ΔT2为T2的测量误差,℃;Δδ1为T1、T2两测点间距离的测量误差,mm。
式中:φT为Tsur的相对误差,%;Δδ0为测点T1与热源表面间距离的测量误差,mm。
式中:φα为α的相对误差,%;φTs为Tsat的相对误差,%。
φTs由方程(7)计算得出:
式中:n为测量次数;Tp为多次测量后得出的平均蒸发温度,℃;ΔTsat,i为每次测得蒸发温度所产生的相对偏差,℃。计算出φTs相对误差值为±3.7%。
通过计算,得到热流密度、热源表面温度、表面传热系数的误差分别为±2.9%、±5.5%、±6.4%。
2 实验结果与分析
在喷雾冷却过程中,喷雾室及热源本身向周围环境的散热量很小,可以忽略不计,因此热源的绝大部分热量被工质带走。本文用热源表面温度Tsur和传热系数α表征喷雾冷却系统换热性能的强弱,下面主要分析过冷度、质量流量、热流密度分别对热源表面温度和传热系数的影响。
2.1 过冷度对换热性能的影响
工质质量流量为3.0 kg/h,热流密度为48.9 W/cm2时,只改变过冷度Tsub,其他条件均保持不变。图3(a)和(b)分别表示热源表面温度和表面传热系数随过冷度的变化。过冷度由0℃增大到6.0℃,喷雾室压力为370 kPa,热源表面温度仅仅降低了1.4℃左右,换热系数略微升高约900 W/(m2·℃)。
实验结果表明单纯增大过冷度并不能明显提高喷雾冷却系统的换热性能。由于工质过冷度增大,提供的显冷量增加,且对流换热温差增大,增加了对流换热量。但是过冷度增大提供的显冷量远小于工质相变提供的潜冷量,即增大过冷度所增加的冷量在总冷量中所占的比例很小,所以增大工质过冷度不能显著提高换热性能。
图3 过冷度对热源表面温度和换热系数的影响Fig.3 The effect of subcooling degree on heat surface temperature and heat transfer coefficient
2.2 质量流量对换热性能的影响
热流密度和工质过冷度分别为38.1 W/cm2和3℃,只改变工质质量流量。实验过程中喷雾室内的压力是因变量,即喷雾室压力由热流密度,质量流量,过冷度等共同决定,其变化规律与热源表面温度变化规律一致,所以本实验不着重强调喷雾室压力。图4 (a)和(b)分别表示热源表面温度和热源表面换热系数变化规律。随着流量的增加,热源表面温度先降低随后略有回升,对应的换热系数先上升然后下降。
图4 质量流量对热源表面温度和换热系数的影响Fig.4 The effect of mass flow rate on heat surface temperature and heat transfer coefficient
当质量流量为1.0 kg/h时,工质流速较小且喷射到热源表面的工质量较少,因此工质经喷嘴雾化成液滴喷射到热源表面进行换热后全部呈气态,出现前面的液滴已经蒸干,后面的液滴还未补充上的现象,热源表面会出现干涸区域,导致临界热流密度出现。所以热源表面温度高达48.7℃,而换热系数仅有9.76 kW/(m2·℃)左右。随着工质质量流量从1.0 kg/h增大到4.0 kg/h,单位时间内喷射到热源表面的液滴量增加,逐渐满足热源的散热需求,热源表面的干涸区面积逐渐减小,同时流量增大使液滴喷射速度增大,增强了对热源表面液膜内工质流动的扰动,强化了对流换热,所以热源表面温度降低,换热系数升高。当质量流量达到4.0 kg/h时,液滴的补充量与蒸发量相当,得到最好的换热效果,热源表面温度降低至29.9℃,换热系数上升到最大值28.3 kW/(m2·℃)。当质量流量上升至6.0 kg/h,热源表面温度升高至33.5℃,换热系数降低至19.97 kW/(m2·℃)。这是因为随着质量流量的增大更多的液滴喷射到热源表面,液滴的补充量增大,而液滴的蒸发量不变,所以过量的液滴在热源表面堆积形成液膜,液膜厚度增加,蒸发换热相对减弱,换热方式以对流换热为主。根据湍流模型的层流底层理论推测出热源表面存在两层液膜:底层液膜和上层液膜。底层液膜是工质附着在热源表面形成一层很薄的液膜,而且这层液膜的流动速度很慢;上层液膜是处于底层液膜之上,厚度较厚且在液滴的冲刷下流动速度较快的一层液膜。上层液膜阻碍了底层液膜蒸发,液膜厚度越厚,阻碍作用越大,热源的热量以热传导的方式传递给底层液膜,然后底层液膜将热量传递给上层液膜,上层液膜以对流换热的方式将热量带走,所以此过程中对流换热代替了蒸发换热成为主要的散热方式。虽然质量流量增大强化了对流换热,但是换热方式的改变以及热阻的增大都阻碍了换热,且阻碍作用大于促进作用,所以换热性能降低。
2.3 热流密度对换热性能的影响
图5 热流密度对热源表面温度和换热系数的影响Fig.5 The effect of heat flux on heat surface temperature and heat transfer coefficient
在过冷度为3.0℃,质量流量为3.0 kg/h的条件下,热源表面温度和表面传热系数系数随热流密度的变化如图5(a)和(b)所示。热源表面温度随热流密度的增大而升高 (如图5(a)所示),热流密度从27.8 W/cm2增加到53.8 W/cm2,热源表面温度由25.9℃上升至43.0℃。因为冷却工质的质量流量是定值,在其他条件不变的情况下,增加热流密度导致热源内部的热量增多,在不改变其他条件的情况下,此部分增加的热量,必然导致热源表面温度上升。
随着热流密度的增大,换热系数先上升后下降(如图5(b)所示)。随着热流密度从27.8 W/cm2增大到43.1 W/cm2,换热系数从19010 W/(m2·℃)增大到23640 W/(m2·℃),热流密度继续增大到53.8 W/cm2,换热系数则降低到18700 W/(m2·℃)。因为在热流密度比较小的情况下,液滴的蒸发量小于液滴的补充量,热源表面形成一定厚度的液膜,阻碍了相变换热,而且增大了热阻。而随着热流密度的增大,液滴蒸发量增大导致液膜厚度变薄,阻碍作用减弱且热阻减小,换热系数增大。随着热流密度的增大,液滴的蒸发量逐渐接近补充量,此时热源表面上液膜足够薄且刚好润湿热源表面,换热系数最高。继续增大热流密度最终使液滴的蒸发量大于补充量,热源表面局部区域出现干涸现象,导致换热系数减小。
3 结论
本文建立了以R134a为工质的闭式喷雾冷却系统,实验研究了过冷度、质量流量以及热流密度对喷雾冷却系统换热性能的影响。主要结论如下:
1)在喷雾冷却系统中,当热流密度和质量流量为定值时,增大过冷度所提供的显冷量远小于液滴相变提供的冷量,因此单纯增大过冷度并不能明显提高喷雾的冷却性能。
2)当热流密度和过冷度为定值时,在质量流量增大的过程中,热源表面温度先降低然后又上升,表面传热系数先升高然后降低,所以存在一个临界质量流量。在质量流量未达到临界值之前,增大质量流量提供的液滴量渐渐满足换热需求,换热性能提高;当质量流量达到临界值时,刚好满足换热需求,换热性能最佳;但是质量流量超过临界值后,热源表面会形成较厚的液膜阻碍蒸发换热,换热性能降低。所以在喷雾冷却系统中并不是流量越大换热效果越好,而是要尽可能使质量流量达到临界值。
3)当质量流量和过冷度为定值时,热源表面温度随着热流密度增大而升高,然而表面传热系数先上升后下降。因为液滴的补充量不变,热流密度较小时,蒸发量小于补充量时热源表面液膜较厚,增大热流密度液膜厚度减小换热系数增大;持续增大热流密度最终会导致蒸发量大于补充量,热源表面出现干涸区域,临界热流密度出现,换热系数降低。
在此喷雾冷却系统中,工质的质量流量控制液滴的补充量,热流密度控制液滴的蒸发量。当液滴的补充量与蒸发量相当且热源表面恰好被完全润湿时,换热系数最大。所以根据热流密度适当地调节工质的质量流量能使喷雾冷却系统的换热性能达到最佳。
符号说明
n——测量次数
q——热流密度,W/cm2
Tp——某一蒸发压力下的平均蒸发温度,℃
Tsat——喷雾室的蒸发温度,℃
Tsur——热源表面温度,℃
ΔT1——T1的测量误差,℃
ΔT2——T2的测量误差,℃
ΔTsat——由蒸发压力得到蒸发温度产生的误差,℃
α——换热系数,W/(m2·℃)
δ0——T1测点和热源表面之间的距离,mm
δ1——T1、T2测点之间的距离,mm
Δδ0——测点T1与热源表面间距离的测量误差,mm
Δδ1——T1、T2两测点间距离的测量误差,mm
φq——q的相对误差,%
φT——Tsur的相对误差,%
φT s——Tsat的相对误差,%
φα——α的相对误差,%
λ——紫铜的导热系数,W/(m·K)
[1] 邵杰,李小莉,冯宇彤,等.激光二极管端面抽运Nd:YVO4板条激光器及其热效应[J].光学学报,2008,28(3):497⁃501.(Shao Jie,Li Xiaoli,Feng Yutong,et al.LD⁃end⁃pumped Nd:YVO4 slab laser and its thermal effects[J].Acta Optica Sinica,2008,28(3):497⁃501.)
[2] 司春强,邵双全,田长青,等.高功率固体激光器喷雾冷却技术[J].强激光与粒子束,2010,22(12):2789⁃2793.(Si Chunqiang,Shao Shuangquan,Tian Chan⁃gqing,et al.Spray cooling technology for high⁃power solid⁃state laser[J].High Power Laser and Particle Beams,2010,22(12):2789⁃2793.)
[3] Xu Hongbo,Si Chunqiang,Shao Shuangquan,et al.Ex⁃perimental investigation on heat transfer of spray cooling with isobutane(R600a)[J].International Journal of Ther⁃mal Sciences,2014,86(12):21⁃27.
[4] 杨波,高松信,刘军,等.高功率二极管激光器喷雾冷却实验研究[J].强激光与粒子束,2014,26(7):3⁃6. (Yang Bo,Gao Songxin,Liu Jun,et al.Spray cooling of high power diode laser[J].High Power Laser and Particle Beams,2014,26(7):3⁃6.)
[5] 王亚青,刘明侯,刘东,等.喷雾冷却换热机理和换热性能的因素[J].强激光与粒子束,2011,23(9):2277⁃2281.(Wang Yaqing,Liu Minghou,Liu Dong,et al.Heat transfer mechanism and influencing factors in spray cooling[J].High Power Laser and Particle Beams,2011,23(9):2277⁃2281.)
[6] Hou Yu,Liu Xiufang,Liu Jionghui,et al.Experimental study on phase change spray cooling[J].Experimental Thermal and Fluid Science,2013,46:84⁃88.
[7] Pautsch A G,Shedd T A.Adiabatic and diabatic measure⁃ments of the liquid film thickness during spray cooling with FC⁃72[J].International Journal of Heat and Mass Trans⁃fer,2006,49(15/16):2610⁃2618.
[8] 王亚青,刘明侯,刘东,等.倾斜喷射时喷雾冷却无沸腾区换热特性[J].化工学报,2009,60(8):1912⁃1919.(Wang Yaqing,Liu Minghou,Liu Dong,et al. Heat transfer performance of inclined spray cooling in non⁃boiling regime[J].CIESC Journal,2009,60(8):1912⁃1919.)
[9] Rini D P,Chen R H,Chow L C.Bubble behavior and nu⁃cleate boiling heat transfer in saturated FC⁃72 spray cooling [J].Journal of Heat Transfer,2002,124(1):63⁃72.
[10]Lin L,Ponnappan R.Heat transfer characteristics of spray cooling in a closed loop[J].International Journal of Heat and Mass Transfer,2003,46(20):3737⁃3746.
[11]Si Chunqiang,Shao Shuangquan,Tian Changqing,et al. Development and experimental investigation of a novel spray cooling system integrated in refrigeration circuit[J]. Applied Thermal Engineering,2012,33/34:246⁃252.
[12]王亚青,刘明侯,刘东,等.高功率激光器喷雾冷却的实验研究[J].强激光与粒子束,2009,21(12):1761⁃1766.(Wang Yaqing,Liu Minghou,Liu Dong,et al.Ex⁃perimental study on spray cooling for high⁃power laser[J]. High Power Laser and Particle Beams,2009,21(12):1761⁃1766.)
[13]钱洋,刘炅辉,李玫,等.采用R134a工质的相变喷雾冷却性能实验研究[J].西安交通大学学报,2015,49 (1):97⁃101.(Qian Yang,Liu Jionghui,Li Mei,et al. Study on heat transfer performance of phase change spray cooling with R134a[J].Journal of Xi’an Jiaotong Univer⁃sity,2015,49(1):97⁃101.)
[14]Tan Y B,Xie J L,Duan F,et al.Multi⁃nozzle spray cooling for high heat flux applications in a closed loop sys⁃tem[J].Applied Thermal Engineering,2013,54(2):372⁃379.
[15]Hou Yu,Liu Jionghui,Su Xuemei,et al.Experimental study on the characteristics of a closed loop R134⁃a spray cooling[J].Experimental Thermal and Fluid Science,2015,61:194⁃200.
[16] Li Qiang,Tie Peng,Xue Yimin.Investigation on heat transfer characteristics of R134a spray cooling[J].Experi⁃mental Thermal and Fluid Science,2015,60:182⁃187.
[17]Hsieh S S,Tien C H.R⁃134a spray dynamics and im⁃pingement cooling in the non⁃boiling regime[J].Interna⁃tional Journal of Heat and Mass Transfer,2007,50(3/4):502⁃512.
[18]Eduardo Martínez⁃Galván,Ramos J C,Raúl Antón,et al. Influence of surface roughness on a spray cooling system⁃with R134a.Part I:heat transfer measurements[J].Ex⁃perimental Thermal and Fluid Science,2013,46:183⁃190.
[19]Moffat R J.Describing the uncertainties in experimental re⁃ sults[J].Experimental thermal and fluid science,1988,1(1):3⁃17.
徐洪波,男,副研究员,中国科学院理化技术研究所,(010)82543433,E⁃mail:hbxu@mail.ipc.ac.cn。研究方向:制冷空调新技术、微通道换热、喷雾冷却及热泵技术。
About the corresponding author
Xu Hongbo,male,associate researcher,Technical Institute of Physics and Chemistry,CAS,+86 10⁃82543433,E⁃mail:hbxu @mail.ipc.ac.cn.Research fields:new technology of refrigera⁃tion and air conditioning,micro⁃channel heat transfer,spray cool⁃ing and heat pump technology.
Investigation on the Heat Transfer Performance of Spray Cooling System with R134a
Qian Chunchao1,2Xu Hongbo1Shao Shuangquan1Tian Changqing1Zhou Guanghui2
(1.Beijing Key Laboratory of Thermal Science and Technology and Key Laboratory of Cryogenics,Technical Institute of Physics and Chemistry,CAS,Beijing,100190,China;2.School of Energy and Environment Engineering,Zhongyuan University of Technology,Zhengzhou,450007,China)
Water In order to investigate the influences of subcooling degree,mass flow rate and heat flux on the performance of the spray cooling system,a closed spray cooling system with R134a was built.The subcooling degree was adjusted by the subcooled section located before the nozzle inlet,the mass flow rate was controlled by the variable gear pump,and the heat flux was dominated by changing the voltage and current of the heating power.The experimental results indicated that,the variation of subcooling degree did not have significantly influence on the heat surface temperature and heat transfer coefficient as the heat flux and the mass flow rate were fixed.When the heat flux and the subcooling degree were constant,there existed a critical value of the mass flow rate,as the mass flow rate was smaller than the critical val⁃ue the heat surface temperature decreased with the mass flow rate increasing,on the contrary,the heat surface temperature would increase with the mass flow rate increasing when the mass flow rate exceeded the critical value.When the mass flow rate and the subcooled temper⁃ature kept constant,there also existed a heat flux value which made the evaporation capacity equal to the supplement of droplets,therefore the heat transfer coefficient could reach the maximum value.
spray cooling;heat transfer coefficient;subcooling degree;heat flux;mass flow rate;R134a
TK124;TB61+2
A
0253-4339(2015)04-0001-07
10.3969/j.issn.0253-4339.2015.04.001
简介
国家自然科学基金(51106170)资助项目。(The project was supported by the National Natural Science Foundation of China(No. 51106170).)
2014年11月24日