高马赫数来流超燃冲压发动机燃烧流场分析
2017-03-27张时空黄志伟
张时空,李 江,黄志伟,秦 飞,薛 瑞
(西北工业大学燃烧、热结构和内流场重点试验室,西安710072)
高马赫数来流超燃冲压发动机燃烧流场分析
张时空,李 江,黄志伟,秦 飞,薛 瑞
(西北工业大学燃烧、热结构和内流场重点试验室,西安710072)
以模拟自由来流马赫数12的地面试验氢燃料超燃冲压发动机为研究对象,应用商用计算流体力学软件CFD++;针对高马赫数来流下的超燃冲压发动机的典型流场结构、空间释热分布、预混/非预混燃烧模式和火焰稳定机理开展了分析研究。计算中采用7组分、9反应步的氢气/氧气动力学模型,使用壁面函数结合两方程剪应力输运模型,基于雷诺时均化方法开展计算,数值结果与试验数据相符较好。1)验证了CFD++软件在高马赫数来流下的适用性和计算精度;2)分析了高超声速来流下的燃烧室流场特征;3)获得了高马赫来流条件下的发动机燃烧效率、释热区间、预混/非预混燃烧模式的空间分布规律;4)为进一步开展高马赫数下的发动机精细化流场计算和多尺度燃烧过程研究提供了重要依据。
超燃冲压发动机;高马赫数;燃烧模式;数值模拟
0 引 言
超声速燃烧是一种燃料在超声速气流中掺混与燃烧的复杂物理化学过程。超声速燃烧冲压发动机主要由进气道、燃烧室、尾喷管三部分组成:进气道将高超声速来流空气压缩到合适的温度与压力后送入燃烧室;燃料与空气在燃烧室内进行掺混、燃烧并释放热能;燃料化学能转变成燃气的动能后通过尾喷管膨胀做功产生推力。超燃冲压发动机工作范围较广,适用飞行速度低至Ma 5,高至Ma 10以上,是实现空天飞行器、临近空间飞行器和高超声速巡航导弹的最有效动力装置之一[1-3]。
当自由来流速度达到Ma 5以上时,发动机内为亚声速与超声速共存的流动状态,来流通过隔离段内的预燃激波串实现减速。当自由来流速度达到高马赫数(Ma 8以上)时,燃烧室内绝大部分区域为超声速流动状态[4]。与较低马赫数相比:发动机由双模态模式转变为纯超燃模式进行工作,流道内波系复杂,呈现明显的三维流场特征;飞行器的内外阻较大,Ma 8速度下,产生1份净推力需要用7份推力克服6份阻力,同时燃料掺混困难,增加掺混与减小阻力呈现突出的矛盾[2];发动机热环境较为严酷,必须采用具备较高热值与热沉的氢燃料代替碳氢燃料,对机体进行主动冷却的同时产生更高的推力[5];当来流总温在2500 K以上时,高温真实气体效应开始凸显[4]。伴随着高马赫数而出现的物理现象对该状态下超燃冲压发动机的设计提出了新的要求,原有计算分析手段在该条件下的适用性和准确性需要重新评判。
近年来,对Ma 8以上超燃冲压推进系统的研究日渐兴盛[2],试验器的最高飞行速度不断刷新[6]。由于超高速风洞的造价昂贵,目前仅有美国[6]、澳大利亚[7]、日本[8]等国开展了Ma 8以上的相关试验。国内对于Ma 8的超燃冲压发动机的研究刚刚兴起[9-10],更高马赫数的研究尚缺乏相关报道。相比其他方案,M12系列发动机模拟自由来流达到Ma 12,直连试验进气道入口速度达到Ma 6.72,这一数值已经处于双模态冲压发动机的自由来流速度区间,针对如此高速下的发动机燃烧流场的分析,国外鲜见,国内未见公开报道。
由于超声速燃烧过程在时空上具有强瞬变、强湍流、强压缩、各向异性和一些物理化学现象相互耦合的非线性特点,任何侵入式的物理传感器在对超声速流场带来激波干扰的同时也面临着生存问题,计算流体动力学(Computational fluid dynamics,CFD)的作用较为重要[11],CFD逐渐成为除试验外唯一可用的工具[2].
本文以验证数值模型对高马赫数超燃冲压发动机流场仿真的适用性为切入点,应用CFD++软件,针对日本宇宙航空研究开发机构(Japan Aerospace Exploration Agency,JAXA)的M12-02发动机开展研究。验证高马赫数下选取的湍流模型、燃烧模型、动力学模型以及差分格式的适用性和精确度,分析流场的基本结构、空间释热区间、预混/非预混模式燃烧的分布,对于研究高马赫数下的燃烧过程起到抛砖引玉的作用。
1 构型介绍
M12-02发动机隶属JAXA的M12系列发动机[8]。M12的研制目的为:获得设计点为Ma 12、工作范围为Ma 10~Ma 15的超燃冲压发动机。试验在自由活塞式高焓激波风洞(High Enthalpy Shock Tunnel,HIEST)进行。如图1所示,M12-02发动机为二元构型;使用二维侧压式进气道,等直燃烧室和直扩喷管;应用流向涡支板喷注燃料,以增强燃气掺混,改善发动机的火焰稳定性和点火能力。试验工况如表1所示。
2 数值方法与计算模型
2.1 计算方法
使用商用软件CFD++开展数值模拟工作。CFD++是Metacomp Technologies公司开发的一款流体计算软件平台[12]。本文使用二阶总变量衰减(Total Variation Diminishing,TVD)格式差分求解雷诺时均的N-S方程(RANS),应用有限速率 燃 烧 模 型,使 用 Harten-Lax-Van-LeerC (HLLC)类型的黎曼积分求解方式[13]。CFD++中可以使用一系列的湍流模型[12],本文使用SST模型封闭湍流项[13-15],所有的计算中,湍流生成和离散项都经过了压缩性效应的修正。目前,CFD++已经被应用于超燃流场数值模拟分析[7,16],但对Ma 10以上自由来流的燃烧内流场的性能计算尚不多见。
2.2 网格划分
因发动机为二元构型,选取发动机宽度方向一半的区域进行网格划分,如图2所示。计算区域包括两排共计六个支板。考虑可用计算资源,加密近壁面网格,保证全流道Y+小于5,整个构型网格总数约为1000万。
2.3 边界条件设置
进气道使用超声速入口条件,具体参数由设备喷管计算程序获得,见表1。计算过程中考虑高温真实气体效应(组分离解),设备喷管出口的NO被计入 N2中,因此,O2的质量分数仅为0.2031[8]。
尾喷管使用超声速出口条件,采用常温300 K、黏性无滑移壁面。使用经验壁面函数法求解近壁面湍流[12]。燃料支板使用质量流量与温度的入口条件,单个支板喷注燃料0.00752 kg,温度250 K。
表1 进气道入口状态[8]Table 1 Inlet entrance conditions[8]
2.4 化学动力学模型
使用7组分、9步的简化动力学模型[17],如表2所示。由于进气道入口气流的速度较高,气流通过与发动机等长度的距离所用时间在毫秒量级,已经接近化学动力学的时间尺度,故研究中使用有限速率燃烧模型,以考虑燃料的点火延迟。
为验证简化模型的适用性,比较简化模型与9组分27步详细模型[18]点火延迟时间,如图3所示,可见两者相符较好,简化模型可以用于流场计算。
表2 H2-O2化学反应模型[12]Table 2 Hydrogen-oxygen reaction mechanism[12]
2.5 数值结果/试验数据比对
计算结果中,状态方程的各项残差均小于10-5且不再变动,进出口质量流率之差小于10-6g,认为结果已经收敛。将数值结果与文献[8,19]中的试验数据进行对比(试验中所有测点处的压力数据均取试验时间内的平均值),图4为对比结果。
由图4可见,数值模型能较好预测燃烧室的壁面静压,但对燃烧室前半部分的压力预测稍差。
由壁面静压可见,与低来流马赫数亚燃冲压、双模态冲压燃烧室[16,20]相比,高马赫数来流下,燃烧流场中明显的波系贯穿整个流道。燃烧室的压力分布间接反映了空间释热量的分布。由此初步判定,高马赫数下发动机的燃烧模式与常规的亚燃、双模态冲压发动机相比具有显著不同的特点。
3 数值结果分析
以下从燃烧效率、流场结构、燃烧模式及燃烧释热区间入手,初步揭示高马赫数下超燃冲压发动机的燃烧、流动特征。
3.1 发动机燃烧效率
本文中将燃烧效率定义为,燃烧室任意横向截面上燃烧最终产物中H2O组分的实际质量分数与理论质量分数的比值[21]。即
式中:YH2O和YH2分别为H2O和H2的质量分数,ρ为气体密度,ux为气流轴向速度,Ayz为垂直于轴线的任意横向截面积,νH2O和νH2分别为H2单步反应时H2O和H2的化学计量数,WH2O和WH2分别为H2O和H2的分子质量。
图5为燃烧效率与释热量沿流道的分布曲线,图中可见,支板出口处 (x=0.200 m)至 x= 0.259 m处燃烧效率和释热量基本为0,化学反应很微弱,属于反应感应区(A区)。从x=0.259 m至x=0.579 m处燃烧效率快速增长,属于快速反应区(B区)。从x=0.579m至发动机出口处 (C区)燃烧效率沿程缓慢增加,气流处于短距离内的局部平衡状态。由释热率分布可见,流道内的强激波系对于燃气的影响明显,在x=0.299~0.429 m区间,燃气释放大量热量;在x=0.434~0.484 m区间,受波系影响,燃气吸热;在x=0.484~1.900 m区间,燃气持续受到上述波系加热和燃烧释热的双重影响,二者此起彼伏,造成流道内释热量的规律波动。在x=1.90 m处,燃烧效率达到0.810,此后超声速气流在喷管内加速减压,反应向放热方向进行,最终在喷管出口处(x=2.353 m)燃烧效率达到0.837,在高气流速度导致的燃料短驻留时间[22]与高气流总温导致的化学反应进度受到抑制等因素的共同作用下,发动机的燃烧效率较低。
3.2 流场结构分析
图6为冷流/燃烧状态下中心截面的压力分布云图,可见在燃烧条件下,燃烧室内仍存在着明显的波系结构。
高速来流先后受到进气道与支板型面的压缩作用,压力升高,在进气道内与燃烧室之前分别形成两道斜激波。进入燃烧室后,由于流道面积突扩,形成两道膨胀波;在支板喷嘴出口处,氢气的喷入对来流产生横向的剪切作用,形成两道弓形激波。此后,初始形成的激波和膨胀波在燃烧室壁面多次反射、相交并叠加,与燃烧释热形成强烈的耦合作用,最终建立起贯穿整个燃烧室的复杂波系。
在来流高超声速状态下,化学反应强度的相对减弱加上高超声速剪切层具有的较强稳定性,燃烧反应不足以充分破坏波系结构,流场呈现出与低超声速来流时[23]显著不同的特点。在燃烧室后段,激波系的强度减弱、厚度增加,燃烧释热在一定程度上“抹平”了波系,当气流进入到大膨胀比的尾喷管时,受膨胀减压的作用,波系遭到显著破坏,波系强度较弱。
图7为流道中心截面上冷流和燃烧工况下的局部马赫数云图,图中将亚音速区域隐藏。由图可见,在冷/热工况下,支板底部和近壁面边界层存在明显的亚声速区域,流道中其余部分为超声速流动,中心流道的气流速度变化较小。燃烧工况下,支板处的燃料与高速来流作用形成了弓形激波,故支板底部亚声速区域较大;近壁面气流速度较低,这反映出室压对于高速来流的减速作用。燃烧增大了流道中的激波角度。
图8为支板后方 H2质量分数的等值面分布,等值面以温度着色。从图8(a)~图8(d)可见,氢气射流受到高速气流的剪切与燃料支板的扰动作用,氢气在支板轴线上沿流向呈明显的带状分布。图中清楚地再现了复杂波系在壁面上的反射现象。H2的消耗主要发生在近壁面区域,沿壁面向高度方向扩散,然而由于超声速气流的刚性,H2直至燃烧室出口也无法实现高度方向的贯穿(图8(a))。从图8(c)~图8(d)可见,下游流道的燃气温度较高,燃料在燃烧室前半段放热较少,近支板处温度较低。这一方面是因支板附近为富燃环境,H2/O2掺混不均;另一方面,流道内气流速度较高,燃气在点火时间内需要流经一定的距离以完成自点火,活性基团产生后,其在燃气中的扩散同样需要一定的时间,故高温区间主要形成在燃烧室下游。
图9(a)~图9(d)分别为x=0.2~0.5 m截面处的湍动能云图,其中白色箭头表示速度矢量,箭头长短对应速度大小。由图9(a)可见,燃料支板扰动附近气流,形成了流向涡。在流道壁面附近(图右侧),高速来流受到壁面的摩擦减速作用,加之激波在流道壁面的反射,流道中心部分在激波相交点附近产生流动分离,形成了高湍动能区域(图9(a),9(c)),故壁面附近湍动能与流道中心相比较大。从图9各支板附近的速度矢量可见,在支板出口处,上下排流向涡互不干扰,当流动向下游发展时,流向涡的卷吸运动逐渐衰减,同时流向涡的扰动不断向外围扩展,在x=0.5 m贯穿了流道的高度方向。
3.3 燃烧模式分析
式中:YF、YO为燃料和氧化剂质量分数。当火焰指数为1时,为预混燃烧模式;当火焰指数为-1时,为非预混燃烧模式。
图10为流道中心截面火焰指数云图,可见流道内以非预混燃烧占主导,预混燃烧主要集中在x= 0.91 m之前的近壁面部分。在流道中心,燃料与空气受到激波的隔离作用,在燃烧室入口处到x= 0.41 m之前形成了不反应区域。
图10中,燃料沿燃烧室壁面向下游发展、与空气掺混,很快产生自点火,形成局部的预混燃烧区域;该区域在下游沿轴向和高度两方向进行扩展,起到了较好的加热、预燃、点火作用。而在燃烧室的其他区域,由于高超声速剪切层的强稳定性,加上高超声速来流条件下燃烧室具有较强的自点火特性,燃料在未实现充分掺混的条件下即迅速发生化学反应,因而以非预混燃烧为主导。
图11为流道中心截面释热率云图,图中隐藏了释热值为0以上(燃气吸热)的区间。由图可见,流道内释热区间呈菱形分布。在x=0.53 m之前,释热主要发生在近壁面区域;在x=0.53 m之后,随着燃气组分的扩散,流道中部也出现了释热。在x= 0.20 m(支板出口)到x=0.22 m之间,无热量释放;在x=0.22~0.50 m内,热量剧烈释放(颜色较深部分);在x=0.50 m之后,由于激波在流道内反射相交,在空白区域中,激波在交汇时释放的热量超过了燃气自身的化学反应释热;在尾喷管入口,高速气流形成了上下两道膨胀波,燃气加速降温,促进了释热反应进行,故在整个尾喷管内,燃气持续释放热量。
图12(a)~(d)分别对应局部流道中心截面上的H2组分,OH组分、火焰指数和静压(P)云图。由图12(a)可见,氢气在支板后方形成高浓度区间,而在该区间至x=0.23 m前的主流中,OH生成量很少(图12(b)),即燃料与空气反应较弱;由火焰指数与OH组分可见,流道近壁面区域(x=0.23~0.38 m)为预混燃烧模式,有明显的OH存在,OH生成后被输运至流道中部参与扩散燃烧;与静压云图类似,OH的分布呈现周期性。由OH、火焰指数、静压云图可见,x=0.28~0.38 m为预混燃烧区域,也是激波与燃烧室边界层相干涉的地方,这一区间压力较高,有大量OH组分生成;由图12(d)可见,x =0.72~0.76 m近壁面区域,同样为激波与边界层相干涉的位置,此区间亦是预混燃烧区域。反之,在x=0.38~0.48 m近壁面区域,虽然有一定浓度的H2存在,但这一区间OH生成较少。故认为激波一定程度上促进了预混燃烧的发生,同时预混燃烧生成的OH活性组分随着流动输运至流道中后部,进行扩散燃烧。流道内预混燃烧区间虽然较小,却是燃烧室整个化学反应的开始,强烈的预混燃烧改变了燃烧室入口来流状态,预混燃烧对整个发动机内的燃烧和流动过程十分重要。
图13为流道高度方向(y向)50%截面上的释热分布统计,可直观显示亚、超声速与释热的关系。x坐标为马赫数,y坐标为释热率,散点用OH的质量分数染色。由图13可见,燃烧在亚声速和超声速环境下同时进行;在亚声速区域有强烈的热量释放,而大部分的热量在超声速下释放,这也对应着OH分数较高的区域。
4 结束语
1)CFD++软件应用于高马赫数来流下的超声速燃烧流场分析,计算结果与实验数据相符较好。
2)高马赫数发动机内流场中,燃料的大部分化学能以超声速燃烧的形式释放;流道内以非预混燃烧为主,预混燃烧起到“自点火源”的作用;预混燃烧主要集中在燃烧室前半段的近壁面区域;激波诱导边界层分离对于预混燃烧模式具有重要影响。
3)高马赫数来流条件下,燃烧室内的波系结构贯穿整个流道,燃烧释热对于流场结构的改变作用下降,波系对整个燃烧过程的影响显著。
[1] 秦飞,何国强,刘佩进,等.圆形燃烧室支板火箭超燃冲压发动机数值模拟[J].固体火箭技术,2011,34(2):150-155.[Qin Fei,He Guo-qiang,Liu Pei-jin,et al.Numerical simulation of strut-rocket scramjet with circular combustor[J].Journal of Solid Rocket Technology,2011,34(2):150-155.]
[2] 俞刚,范学军.超声速燃烧与高超声速推进[J].力学进展,2013,43(5):449-471.[Yu Gang,Fan Xue-jun.Supersonic combustion and hypersonic propulsion[J].Adv.Mech.,2013,43(5):449-471.]
[3] Erik A,Ajay K,Alan W.Study of forebody injection and mixing with application to hypervelocity air breathing propulsion[C].48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference&Exhibit,Georgia,USA,July 30-August 1,2012.
[4] Heiser W,Pratt D,Daley D,et al.Hypersonic airbreathing propulsion[M].Washington DC:AIAA,1994.
[5] 黄伟,罗世彬,王振国.临近空间高超声速飞行器关键技术及展望[J].宇航学报,2010,31(5):1259-1265.[Huang Wei,Luo Shi-bin,Wang Zhen-guo.Key techniques and prospect of near-space hypersonic vehicle[J].Journal of Astronatics,2010,31(5):1259-1265.]
[6] Rogers R C,Shih A T,Hass N E.S cramjet development tests supporting the Mach 10 flight of the X-43[C].AIAA/CIRA 13th International Space Planes and Hypersonic Systems and Technologies,Capua,Italy,May 16-20,2005.
[7] Mathew G B,Laurie M B,Russell R B.Supersonic combustion processes in a premixed three-dimensional nonuniform compression scramjet engine[J].AIAA Journal,2014,52(8): 1670-1684.
[8] Masatoshi K,Vigor Y,Masahiro T,et al.Ignition transient phenomena in a scramjet engine at Mach 12 flight condition[C].43rd AIAA/ASME/SAE/SAEE Joint Propulsion Conference&Exhibit,OH,USA,July 8-11,2007.
[9] 周建兴,汪颖.高马赫数超燃冲压发动机性能数值研究[J].推进技术,2014,35(4):433-441.[Zhou Jian-xing,Wang Ying.Numerical investigation on performance of a high Mach number scramjet[J].Journal of propulsion technology,2014,35(4):433-441.
[10] 王培勇,陈明,邢菲,等.Hyshot超燃冲压发动机的数值模拟[J].航空动力学报,2014,29(5):1021-1023.[Wang Pei-yong,ChenMing,XingFei,etal.CFD numerical simulation of hyshot scramjet[J].Journal of Aerospace Power,2014,29(5):1021-1023.]
[11] 金亮,梁剑寒,罗世彬,等.超燃冲压发动机燃烧室流动问题的数值验证[J].宇航学报,2008,29(6):1922-1926.[Jin Liang,Liang Jian-han,Luo Shi-bin,et al.CFD validation for typical combustor flow fields in scramjet[J].Journal of Atronautics,2008,29(6):1922-1926.]
[12] METACOMP,CFD++User Manual,Version 11.1[M].Agoura Hills,California:Metacomp Technologies,2011.
[13] Batten P,Leschziner M A,Goldberg U C.Average-state jacobians and implicit methods for compressible and turbulent flows[J].Journal of Computational Physics.1997,137(1): 38-78.
[14] Menter F R,Kuntz M,Langtry R.Ten years of industrial experience with the SST turbulence model[J].Heat and Mass Transfer,2003,4(1):625-632.
[15] Batten P,Craft M,Leschziner M,Loyau H.Reynolds-stresstransport modelling for compressible aerodynamics applications[J].AIAA Journal,1999,37(7):785-797.
[16] Robert J Y,Datta V G.Numerical investigation of dual-mode operation in a rectangular scramjet flow path[J].Journal of Propulsion and Power,2014,30(2):474-489.
[17] Wang T S,Farmer R C,Edelman R B.Turbulent combustion kinetics forcomplexhydrocarbon fuels[C].26th AIAA Aerospace Sciences Meeting,Nevada,USA,January 11-14,1988.
[18] Marinov N M,Westbrook C K,Pitz W J.Detailed and global chemical kinetics model for hydrogen.[C].8th International Symposium on Transport Properties,San Francisco,USA,October,1995.
[19] Scott R,Masahiro T,Tetsuji S,et al.Viscous drag reduction in a mach 12 scramjet engine[C].AIAA/CIRA 13th International Space Planes and Hypersonics Systems and Technologies,Capua,Italy,May 16-20,2005.
[20] Sarah B,Axel V R.Numerical study of the hydrogen/air combustion with cedre on laerte dual mode ramjet combustion experiment[C].20th AIAA Space Planes and Hypersonic Systems and Technologies,Glasgow,Scotland,July 6-9,2015.
[21] Liu G,Xu X,Xie Y F.Numerical investigation on the supersonic combustion of liquid kerosene in a dual-staged strut based scramjet combustor[C].50th AIAA/ASME/ASE/SAEE Joint Propulsion Conference,OH,USA,July 28-30,2014.
[22] 李佩波,王振国,孙明波,等.超声速气流中液体横向射流的气液相互作用过程数值研究[J].宇航学报,2016,37 (2):209-215.[Li Pei-bo,Wang Zhen-guo,Sun Ming-bo,et al.Numerical simulation of the gas-liquid interaction of cross liquid jet in supersonic flow[J].Journal of Astronautics,2016,37(2):209-215.]
[23] Huang Z W,He G Q,Qin F,et al.Large eddy simulation of flame structure and combustion mode in a hydrogen fueled supersonic combustor[J].International Journal of Hydrogen Energy,2015,40(31):9815-9824.
[24] Yamashita H,Shimada M,Takeno T.A numerical study on flame stability at the transition point of jet diffusion flame[J].26th Symposium(International)on Combustion,1996,26(1): 27-34.
通信地址:西安友谊西路127号西北工业大学l64信箱
电话:(029)88460327
E-mail:qadr@mail.nwpu.edu.cn
(编辑:张宇平)
Combustion Flow Field Analysis of a Scramjet Engine at High Mach Number
ZHANG Shi-kong,LI Jiang,HUANG Zhi-wei,QIN Fei,XUE Rui
(Science and Technology on Combustion,Internal Flow and Thermo-structure Laboratory,Northwestern Polytechnical University,Xi’an 710072,China)
A hydrogen-fueled ground testing scramjet engine operating at flight Ma 12 is numerically studied.Computation of the combustion field is carried out based on the commercial software CFD++.The basic flow structure,the spatial distribution of heat release,the combustion mode of premixed/non-premixed,and the mechanism of flame stabilization of the scramjet engine which works at Ma 12 are investigated in detail.A hydrogen-air chemistry kinetic model consisting of 7 species,9 reaction steps is adopted in the Reynolds-Averaged Navier-Stokes simulation.The two-equation shear-stress-transport(SST) turbulence model which takes account of the wall functions is used to handle the turbulence-chemistry interactions.The results are validated by the experimentally measured wall pressure distribution with good agreement.1)The applicability and accuracy of CFD++ in the study of compressible reactive flows under high flight Mach numbers are verified;2)the combustion characteristics under hypersonic inflow conditions are studied;3)the combustion efficiency,the heat release,and the combustion mode are obtained for the very high Mach number;4)these observations offer insight into the potential for exploring the mechanisms of the combustion simulation and multi-scale combustion dynamics under hypersonic flow conditions.
Scramjet;High Mach number;Combustion mode;Numerical simulation
V235.21
A
1000-1328(2017)01-0080-09
10.3873/j.issn.1000-1328.2017.01.011
张时空(1986-),男,博士生,主要从事火箭冲压组合发动机总体设计。
2016-04-15;
2016-11-07
国家自然科学基金(91541110)