文章快速检索  
  高级检索
两电极等离子体合成射流性能及出口构型影响仿真研究
周岩, 刘冰, 王林, 罗振兵, 夏智勋    
国防科学技术大学高超声速冲压发动机技术重点实验室, 湖南长沙 410073
摘要: 通过将火花放电的物理效应等效为气体焦耳加热的过程,在能量方程中引入能量源项,进行了单次能量沉积下两电极等离子体合成射流的唯象模拟。提出将等离子体合成射流对外界流场的动能和热能注入分别作为表征射流"冲击效应"和"热效应"的参数。研究表明在单次放电条件下射流建立的自维持振荡过程中,射流动能和热能主要集中于主射流阶段且射流的"冲击效应"相比"热效应"衰减更快,在一个大气压下两电极激励器总的能量转化效率约为2.3%。分析了出口构型对射流的影响,研究表明收缩孔结构可以有效提高射流速度,但将导致射流动量和饱和频率的降低。
关键词: 等离子体合成射流     唯象模拟     能量效率     出口构型     性能    
Numerical simulation of performance characteristics of two-electrode plasma synthetic jet and the influence of different actuator orifice shapes
Zhou Yan , Liu Bing, Wang Lin, Luo Zhenbing , Xia Zhixun     
Science and Technology on Scramjet Laboratory, National University of Defense Technology, Changsha 410073, China
Abstract: Plasma synthetic jet (some called "spark jet" or "pulsed-plasma jet") is a new type of plasma aerodynamic actuation. It is a synthetic jet that is generated by striking an electrical discharge in a small cavity and the gas in the cavity spurts out through a small orifice in a high speed after pressurization owing to the heating caused by electrical discharge. To study the performance characteristics of two-electrode plasma synthetic jet and the influence of different actuator orifice shapes, a phenomenological simulation of two-electrode plasma synthetic jet built by a single energy deposition was accomplished by equating the physical effects of the spark discharge with gas Joule heating and adding source term in energy equation. The kinetic energy and heat injection into outer flow field of plasma synthetic jet were used as characterizations of "impact effect" and "thermal effect". The results show that the kinetic energy and heat injection of plasma synthetic jet are mostly concentrating in primary jet during the self-sustained oscillation process established by a single discharge, and "thermal effect" of jet works longer than "impact effect". Overall energy efficiency of two-electrode plasma synthetic jet actuator is about 2.3% at 1atm. The influence of different actuator orifice shapes was studied and it shows that shrinking orifice can effectively improve the jet velocity, but will reduce jet momentum and saturation frequency simultaniously.
Keywords: plasma synthetic jet     phenomenological simulation     energy efficiency     orifice shape     performance characteristics    

 0 引 言

对流场进行有效操控具有重要的应用价值。新型流动控制技术的研究对于改善飞行器的气动性能、提高飞行器的安全性和可操作性具有重要意义。目前,主动流动控制技术受到越来越广泛的关注。

主动流动控制激励器的设计和研究是主动流动控制发展的核心问题之一,等离子体气动激励器是目前最引人关注的一种主动流动控制激励器[1, 2, 3]。根据放电特性和工作原理的不同,等离子体气动激励器主要分为介质阻挡放电(DBD)[4, 5, 6, 7]、直流/准直流电弧放电[8, 9]和等离子体合成射流[10, 11]等几种类型。其中,等离子体合成射流激励器又称为火花放电式合成射流激励器,它融合了合成射流与等离子体激励器两者的优势,克服了常规等离子体激励器诱导气流速度较低的不足,因而在高速流动控制领域表现出良好的应用前景[1, 12]

两电极等离子体合成射流激励器的结构和工作示意图如图 1所示,它由开有射流出口的绝缘腔体和一对电极组成。其工作过程分为三个阶段:1)能量沉积阶段,通过在两电极之间加载高压,在腔体内产生火花放电,电加热作用使得腔内气体的温度和压力快速升高;2)射流喷出阶段,腔内气体从出口高速喷出,产生等离子体射流;3)吸气复原阶段,由于高速射流的引射导致腔体内形成负压,外部气体重新充填腔体,准备进入下一工作周期。

图 1 两电极等离子体合成射流激励器示意图 Fig. 1 Schematic of two-electrode plasma synthetic jet actuator

美国约翰霍普金斯大学[10, 13, 14]于2003年首先开始了两电极等离子体合成射流激励器的研究,随后德克萨斯大学[12, 15, 16, 17]、法国宇航研究中心[18, 19]、新泽西州立大学[20, 21]及空军工程大学[22, 23]、国防科学技术大学[11, 24]、南京航空航天大学[25]等单位也相继开展了相关试验和仿真研究。目前,试验研究主要是采用纹影和PIV技术对射流瞬时流场结果进行观察,以及对激励器放电时的电参数进行测量,但是由于等离子体合成射流激励器工作空间狭小、电磁干扰大、射流流场变化剧烈等原因,射流的动量、质量流量、温度等参数不易精确测得。因此,数值仿真对于指导激励器的设计、研究不同工作参数对激励器工作性能的影响规律十分必要。

目前,针对等离子体合成射流的数值仿真方法主要分为两类:一是将激励器气体放电过程简化为对腔体的等容瞬时加热,并将按理想气体状态方程计算的激励器腔体达到的理论高温、高压值作为射流形成和发展的初始条件,进行理想气体条件下等离子合成射流的模拟;二是将激励器气体放电过程等效为功率密度按一定时间和空间分布的能量注入过程,通过在能量方程中添加源项的方法进行模拟,此方法相较于前一种方法更能反应激励器的真实工作过程,因而计算结果具有更高的可信度。文献[23]采用第一种方法进行了射流演化过程的模拟,并分析了放电持续时间的影响。文献[24, 25]采用第二种方法研究了放电频率、放电电能、激励器尺寸等参数对射流性能的影响。射流出口是激励器结构中一个关键部分,出口构型的不同对于激励器的工作性能具有重要影响,文献[24]等研究了不同出口直径的影响,结果表明小的激励器出口直径可以产生速度更高、饱和频率较小的射流,但是针对孔的收缩角不同对激励器性能的影响还缺乏相关研究。

本文采用文献[24]在能量方程中添加源项的方法,进行了等离子体合成射流的唯象模拟,针对三维模型网格量较大的问题,采用了简化的二维轴对称模型,缩短了计算周期。通过与文献[16]试验数据的对比,验证了计算方法的可行性。通过对射流参数进行详细分析,结合等离子体流动控制的机理,提出将等离子体合成射流对外界流场的动能和热能注入分别作为表征射流“冲击效应”和“热效应”的参数,并由此定义了激励器总能量效率的计算方法。同时针对之前研究中存在的不足,进一步研究了激励器射流出口收缩角对工作性能的影响。

1 物理模型和计算方法 1.1 控制方程

由于等离子体合成射流涉及流体力学、电磁学、等离子体物理学等多个学科,对其建立精确的物理仿真模型十分困难。因此本文采用了文献[24]中的物理模型和简化假设,控制方程为非定常可压缩粘性N-S方程组,通过在N-S方程组的能量方程中添加能量源项的方法来模拟放电过程中的热量注入,采用有限体积法对控制方程进行离散,空间离散采用二阶迎风格式,时间离散为二阶隐式格式,计算时间步长取为2×10-9s,每个时间步内迭代20步,使得所有变量迭代计算残差小于10-6以保证收敛。

1.2 计算域及网格划分

计算域包括激励器腔体、射流出口和外部流场三部分,其网格划分如图 2所示,考虑到计算域的对称性,为减小计算网格数,仅选取射流流场的1/2进行计算。其中激励器的结构与文献[16]试验中所使用的激励器结构相同,具体尺寸如图 2中所示。为了消除外部边界对计算结果的影响,使外边界侧向距激励器出口35 mm,纵向距激励器出口80 mm。

图 2 计算区域及网格划分 Fig. 2 Computational domain and mesh

为获得精细的射流流场结构,对激励器腔体和射流出口部分的网格进行局部加密,网格尺寸约为0.026 mm×0.026 mm。外部流场第一层网格高度取为0.026 mm,网格增长率1.005。计算网格数约为8万。

1.3 边界条件

计算域外边界设置为压力出口边界条件,根据文献[16]中的试验条件,总温为300 K,总压约为4666 Pa。激励器腔体和射流出口边界设置为固体壁面,壁面与外部环境的热交换可表示为

其中φ为有效传热热流密度(W/m2),传热系数h取为8 W/(m2·K)[22]Tw为固体壁面温度,T=300 K为环境大气温度。

根据文献[16],激励器工作过程中可以将整个腔体视为放电通道,并作为能量注入区域,在放电电流为3.5 A条件下单次放电注入电能约为40 mJ,能量注入过程主要集中在放电前5 μs内,电能到气体热能的转换效率约10%。根据基本假设认为气体加热在时间和空间上为均匀分布,可以得到能量注入区域的功率密度Q为:

其中气体加热的效率ηe=10%,注入电能E=40 mJ,腔体体积V=90.5 mm3,注入时间τ=5 μs,计算可得气体加热的功率密度Q=4.42×1010 W/m3

2 计算方法验证及结果分析 2.1 计算方法验证

图 3所示为放电开始后30 μs的流场密度及速度云图。由图 3可知激励器腔体内的气体由于受热快速膨胀而从射流出口高速喷出,射流的最高速度达到358 m/s,并且在射流外形成一道前驱激波,射流两侧形成对称旋涡。

图 3 射流流场密度及速度云图(放电开始后30 μs) Fig. 3 Contours of density and velocity of plasma synthetic jet (30 μs after the start of the discharge trigger)

激励器轴线上流场的压力、温度分布曲线如图 4所示,其中横坐标表示至激励器出口的距离。压力曲线中,气压在至激励器出口约11 mm处出现由P1点突降到P2点,这代表了前驱激波前后压力的变化。温度曲线中,激励器出口温度呈现先降低、后略有回升的变化趋势,如图 4(b)T1至T3点所示,其中T1点到T2点为高温射流主流所在位置,射流锋面大约位于至激励器出口约7 mm的T2点,T2点到T3点温度的微升可能来自于前驱激波造成的气动加热。

图 4 激励器轴线上压力与温度分布曲线 Fig. 4 Pressure and temperature on the axis of actuator

计算得到的不同时刻前驱激波和射流锋面位置与相同条件下试验结果[16]的对比如图 5所示,其中横坐标表示放电开始后的时间,纵坐标表示至激励器出口的距离。由图 5可知,前驱激波位置随时间基本呈线性变化,这表明激波运动的速度保持恒定,计算得到的激波推进速度约为412.7 m/s;而射流锋面的移动速度逐渐减慢。图 5的结果表明计算得到的前驱激波及射流特性与试验结果较为一致,所采用计算方法能满足本文计算要求。分析认为产生误差的主要原因有:1) 计算时假设能量在整个腔体内均匀注入,而实际的放电电弧区域并不会充满整个腔体,因而导致计算得到的前驱激波和射流锋面位置要略大于试验结果;2) 假设能量注入过程在时间上为均匀分布;3) 对计算流场的二维轴对称简化。

图 5 放电开始不同时刻前驱激波与射流锋面位置 Fig. 5 Trajectory of the plasma synthetic jet precursor shock and front jet
2.2 射流性能及出口构型的影响 2.2.1 计算算例

激励器结构参数是等离子体合成射流性能的重要影响因素。选取如图 6所示的三种激励器出口构型作为计算算例,研究不同出口构型对等离子体合成射流性能的影响。其中注入电能E=50 mJ及气体加热的效率ηe=10%保持不变,激励器腔体尺寸均为4 mm×4 mm,环境温度和压力分别为300 K和1 atm。激励器出口喉道高为0.7 mm,底径为Φ 1 mm,收缩角分别为0°、12°和23°,如图 6所示。

图 6 三种激励器出口构型 Fig. 6 Three kinds of actuator orifice shapes
2.2.2 性能评价参数

文献[26]指出,等离子体流动控制的作用主要包含动力效应、冲击效应、热效应三个方面。其中,动力效应指等离子体在电磁力作用下的定向运动,而等离子体合成射流的流动控制主要在无外加强磁场的情况下进行,因此其动力效应较小,主要依靠射流的冲击效应和热效应。

冲击效应取决于等离子体合成射流动量或动能大小,射流动能Ekinetic可以表示为[24]:

其中t1t2为射流喷射的时间,$\dot{m}$为射流的质量流率,ujet为射流的速度。

热效应来自于气体放电加热引起的流场物性参数的改变。在传统航空条件下,这方面的作用很小,但是在临近空间稀薄空气高电离率和高超声速飞行条件下,这方面的作用可能增大。热效应可以通过等离子体合成射流向外部流场的热能注入来衡量,等离子体合成射流向外部流场的热能注入Ethermal可表示为:

其中Tjet为射流的温度,T=300 K为流场的初始温度,射流气体的定压比热CpTjet的函数。

为了对激励器的工作性能进行评价,定义激励器总能量转化效率η为:

其中ηkineticηthermal分别表示电能向射流动能和热能的转化效率。

为保证激励器腔体内有足够的气体工质,优化连续脉冲工作的射流特性,必须合理选择激励器的放电频率。文献[24]定义了能够实现腔体充分回填的激励器最大工作频率为等离子体合成射流激励器的饱和频率fsat,当放电频率超过fsat时,容易造成激励器腔体无法充分回填而出现“哑火”[16]

此外,为了评价激励器的快速响应能力,定义激励器射流动量到达最大值的时间为射流的峰值时间tpeak,峰值时间越短表明激励器的响应越快。

2.2.3 射流性能分析

图 7所示为孔型1结构激励器对应的射流出口速度、射流动能、腔体压强、腔体内剩余气体质量、射流出口温度、射流注入热能等参数随时间的变化。

图 7 孔型1结构激励器参数变化曲线 Fig. 7 Jet and actuator parameters varying with time for orifice shape 1

图 7可知,激励器射流出口速度存在明显振荡,当ta=196 μs时射流出口速度开始为负值,此时主射流喷射结束,外部气体开始回填腔体;当tb=283 μs时回填过程结束,至此激励器完成一个工作周期,并开始进入一个振幅不断衰减的自维持周期性工作过程,直至下一次气体放电开始。激励器的饱和频率fsat=1/tb=3.54 kHz。射流的动能和注入热能绝大部分集中于主射流,在主射流结束后基本降为0。腔体内压强在气体放电阶段由101 325 Pa急剧升高到140 124 Pa,腔体内剩余气体质量由初始时刻的5.92×10-8 kg降低到ta时刻的4.63×10-8 kg,在主射流阶段腔体内最多有21.8%的气体喷出。对比射流出口速度与温度可以发现,射流速度和动能在达到峰值后开始逐渐衰减,而射流温度却有一个高温平台区,其持续时间约为整个主射流喷出阶段,直至接近ta时刻才开始快速衰减,表明射流可以在较长时间内维持高温。同时对比射流动能与射流热能曲线可以发现,射流热能曲线在达到峰值后的下降速度更慢,这表明射流的“热效应”相比“冲击效应”衰减相对缓慢。

2.2.4 出口构型的影响

三种不同出口构型的激励器各项性能参数比较如图 8所示,其中Vmax、$\dot{m}$maxpmax分别表示出口处射流的最大速度、最大质量流率和最大动量。由图 8可知,随着射流出口收缩角的增大,可以明显达到提高射流最大速度的目的,并且射流的峰值时间略有缩短,表明收缩孔结构仍能具有比较好的快响应特性。但是射流最大速度的提高存在一定的限度,由孔型1到孔型2出口顶部面积收缩51%,射流最大速度提高28.3 m/s;而由孔型2到孔型3出口顶部面积进一步收缩33%,射流最大速度仅仅提高了1.4 m/s。说明当超过一定限度时,出口收缩导致的阻塞作用将使射流增速大大放缓。并且相比射流速度的提高,射流最大质量流率的减小更为显著,这导致了射流最大动量的降低。

图 8 三种出口孔型激励器性能参数比较 Fig. 8 Performance comparison of three kinds of actuator orifice shapes

随着射流出口收缩角的增大,射流到达峰值的时间虽然变化不大,但是峰值过后射流的衰减速度却显著放缓,即射流由“短时间大流量”喷射转变为“小流量长时间”喷射。这造成两方面的影响:一是主射流的持续时间显著增加,这表明射流的有效作用时间延长;二是射流的饱和频率降低,限制了激励器工作频率的提高,如图 8(c)所示。通过能量效率的计算可以发现,收缩孔结构由于会造成更多的动能损失因此总能量转化效率略有降低,但总体而言各个孔型效率差别不大,大约为2.3%。

综上所述,为了实现激励器整体性能的最优,需要根据不同的使用要求合理设计射流出口构型。如果射流出口收缩过度,不仅对于射流增速作用不大,而且会造成射流动量和效率的降低,更重要的是导致激励器工作频率的提高受限,而激励器的高频工作特性对于特征频率较高的高速流动控制至关重要。

3 结 论

(1) 通过在能量方程中添加源项的方法进行了两电极等离子合成射流的唯象模拟,与相同条件下试验得到的前驱激波与射流锋面等数据对比表明,本文所采用的等离子合成射流唯象仿真模型可以与试验数据获得较好的吻合,仿真结果能够比较真实地反映等离子体合成射流的特性。

(2) 基于等离子体合成射流流动控制机理,提出将等离子体合成射流对外界流场的动能和热能注入分别作为表征射流“冲击效应”和“热效应”的参数,进而给出了射流的总能量转化效率和反映射流响应特性的峰值时间等评价参数。研究表明射流动能和热能主要集中于主射流阶段且射流的“冲击效应”相比“热效应”衰减更快,在一个大气压下两电极激励器总的能量转化效率约为2.3%。

(3) 在一定程度上采用收缩孔结构可以有效增大射流的最大速度,但是当出口面积收缩过大(收缩率由51%进一步收缩为84%)时,由于阻塞作用显著增强,增速的效果出现明显降低。同时,收缩孔结构会造成射流动量的降低,更重要的是导致激励器工作频率的提高受限,而激励器的高频工作特性对于特征频率较高的高速流动控制至关重要。因此,在激励器实际工作过程中要根据不同的使用要求对激励器的出口结构进行合理优化。

参考文献
[1] Wang L, Luo Z B, Xia Z X, et al. Review of actuators for high speed active flow control[J]. Sci. China Tech. Sci., 2012, 55:2225-2240.
[2] Wang J J, Choi K S, Feng L H, et al. Recent developments in DBD plasma flow control[J]. Progress in Aerospace Sciences, 2013, 62:52-78.
[3] Nie W S, Cheng Y Fa, Che X K. Review on dielectric barrier discharge plasma flow control[J]. Advances in Mechanics, 2012, 42(6):722-734. (in Chinese)聂万胜,程钰锋,车学科.介质阻挡放电等离子体流动控制研究进展[J].力学进展, 2012, 42(6):722-734.
[4] Meng X S, wang J L, cai J S, et al. Flow control over a slender conical forebody by different plasma actuations[J]. Acta Aerodynamica Sinica, 2013, 31(5):647-651. (in Chinese)孟宣市,王健磊,蔡晋生,等.不同形式等离子体激励对细长体分离涡的控制[J].空气动力学学报, 2013, 31(5):647-651.
[5] Zhang P F, Dai C F, Liu A B, Wang J J. The effect of actuator strength on the plasma synthetic jet[J]. Acta Aerodynamica Sinica, 2012, 30(2):228-232. (in Chinese)张攀峰,戴晨峰,刘爱兵,等.激励强度对等离子体合成射流的影响[J].空气动力学学报, 2012, 30(2):228-232.
[6] Xue B M, Yang Y. Numerical research on airfoil leading edge separation control using plasma actuator[J]. Acta Aerodynamica Sinica, 2009, 27(1):1-4. (in Chinese)薛帮猛,杨永.应用等离子体进行翼型前缘分离控制的数值模拟研究[J].空气动力学学报, 2009, 27(1):1-4.
[7] Li Y H, Wu Y, Zhang P, et al. Experimental investigation on airfoil stall separation suppression by plasma actuation[J]. Acta Aerodynamica Sinica, 2008, 26(3):373-377. (in Chinese)李应红,吴云,张朴,等.等离子体激励抑制翼型失速分离的实验研究[J].空气动力学学报, 2008, 26(3):373-377.
[8] Leonov S B, Bityurin V A, Yarantsev D A, et al. High-speed flow control due to interaction with electrical discharges[R]. AIAA 2005-3287.
[9] Cheng Y F, Nie W S, LI G Q. Numercial study of plasma aerodynamic actuation mechanism[J]. Acta Phys. Sin., 2012, 61(6), 060509. (in Chinese)程钰锋,聂万胜,李国强.等离子体气动激励机理数值研究[J].物理学报, 2012, 61(6), 060509.
[10] Grossman K R, Cybyk B Z, Vanwie D M. Sparkjet actuators for flow control[R]. AIAA 2003-0057.
[11] Wang L, Xia Z X, Luo Z B, et al. A three-electrode plasma synthetic jet actuator for high speed flow control[J]. AIAA Journal, 2014, 52(4):879-882.
[12] Narayanaswamy V, Clemens N T, Raja L L. Investigation of a pulsed-plasma jet for shock/boundary layer control[R]. AIAA 2010-1089.
[13] Haack S J, Taylor T, Emhoff J, et al. Development of an analytical sparkJet model[R]. AIAA 2010-4979.
[14] Grossman K R, Cybyk B Z, Rigling M C, et al. Characterization of sparkjet actuators for flow control[R]. AIAA 2004-0089.
[15] Narayanaswamy V, Shin J, Clemens N T, et al. Investigation of plasma-generated jets for supersonic flow control[R]. AIAA 2008-0285.
[16] Narayanaswamy V, Raja L L, Clemens N T. Characterization of a high-frequency pulsed-plasma jet actuator for supersonic flow control[J]. AIAA Journal, 2010, 48(2):297-305.
[17] Narayanaswamy V, Clemens N T, Raja L L. Method for acquiring pressure measurements in presence of plasma-induced interference for supersonic flow control applications[J]. Measurement Science and Technology, 2011, 22:125107.
[18] Caruana D, Barricau P, Hardy P, et al. The "plasma synthetic jet" actuator aero-thermodynamic characterization and first flow control applications[R]. AIAA 2009-1307.
[19] Hardy P, Barricau P, Belinger A, et al. Plasma synthetic jet for flow control[R]. AIAA 2010-5103.
[20] Anderson K. Characterization of spark jet for flight control[D]. New Brunswick:Rutgers, The State University of New Jersey, 2012.
[21] Golbabaei M, Knight D, Anderson K. Spark Jet efficiency[R]. AIAA 2013-0928.
[22] Jia M, Liang H, Song H M, et al. Characteristic of spark discharge plasma jet driven nanosecond pulses[J]. High Voltage Engineering, 2011, 37(6):1493-1498. (in Chinese)贾敏,梁华,宋慧敏,等.纳秒脉冲等离子体合成射流的气动激励特性[J].高电压技术, 2011, 37(6):1493-1498.
[23] Liu P C, Li J, Jia M, et al. Investigation on flow field of the plasma synthetic jet device[J]. Journal of Air Force Engineering University, 2011, 12(6):22-25. (in Chinese)刘朋冲,李军,贾敏,等.等离子体合成射流激励器的流场特性分析[J].空军工程大学学报, 2011, 12(6):22-25.
[24] Wang L, Luo Z B, Xia Z X, et al. Energy efficiency and performance characteristics of plasma synthetic jet[J]. Acta Phys. Sin., 2013, 62(12):125207. (in Cinese)王林,罗振兵,夏智勋,等.等离子体合成射流能量效率及工作特性研究[J].物理学报, 2013, 62(12):125207.
[25] Shan Y, Zhang J Z, Tan X M. Numerical study of the flow characteristics and excitation parameters for the sparkjet actuator[J]. Journal of Aerospace Power, 2011, 26(3):551-557. (in Chinese)单勇,张靖周,谭晓茗.火花型合成射流激励器流动特性及其激励参数数值研究[J].航空动力学报, 2011, 26(3):551-557.
[26] Wu Y. Study on the mechanism of plasma aerodynamic actuation and its application in compressor stability extension[D]. Xi'an:Air Force Engineering University, 2007. (in Chinese)吴云.等离子体气动激励及其扩大压气机稳定性的原理研究[D].西安:空军工程大学, 2007.
http://dx.doi.org/10.7638/kqdlxxb-2014.0037
中国空气动力学会主办。
0

文章信息

周岩, 刘冰, 王林, 罗振兵, 夏智勋
Zhou Yan, Liu Bing, Wang Lin, Luo Zhenbing, Xia Zhixun
两电极等离子体合成射流性能及出口构型影响仿真研究
Numerical simulation of performance characteristics of two-electrode plasma synthetic jet and the influence of different actuator orifice shapes
空气动力学学报, 2015, 33(06): 799-805
ACTA Aerodynamica Sinica, 2015, 33(06): 799-805.
http://dx.doi.org/10.7638/kqdlxxb-2014.0037

文章历史

收稿日期: 2014-05-19
修订日期: 2014-08-21

相关文章

工作空间