文章快速检索  
  高级检索
基于二级声辐射模型的地面声场高效预测方法
王菲 , 史勇杰 , 徐国华     
南京航空航天大学 直升机旋翼动力学国家级重点实验室, 南京 210016
摘要: 建立了可用于直升机运动状态下噪声传播特性分析的高效地面声场预测方法。该方法包括计算旋翼气动、噪声的自由尾迹方法和时域FW-H方程;计入大气、地面等环境因素的声传播模型及地面声场计算模型;为提高计算效率,在声源计算与噪声传播之间引入基于“紧致球声源”的二级声辐射模型;在此基础上,还提出了通过建立“特征参数声辐射球库”以实现直升机运动状态下噪声实时预测的方法。以AH-1旋翼为算例,通过与地面声场直接计算法对比,说明了方法的有效性及高效性;此外,文中还分析了大气、地面等环境因素对噪声传播以及地面声场特性的影响。
关键词: 噪声     旋翼     二级声辐射球     大气吸声     地面吸声、反射    
Efficient ground noise prediction method based on secondary noise radiation model
WANG Fei , SHI Yongjie , XU Guohua     
National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
Received: 2015-11-17; Accepted: 2016-02-29; Published online: 2016-04-13
Foundation item: National Natural Science Foundation of China (11302103)
Corresponding author. Tel.:025-84892117,E-mail:shiyongjie@nuaa.edu.cn
Abstract: An efficient ground noise prediction method for the analysis of helicopter noise in flight based on the secondary noise radiation is developed. The method to predict the aerodynamic force and aeroacoustic characteristics based on the free-wake/FW-H, the noise propagation model that can take the effects of atmosphere absorption, ground absorption and reflection on noise radiation characteristics into account, as well as ground noise computation model are included in this method. In order to improve computational efficiency, the secondary sound radiation model built on the hypothesis of compact source is introduced between the noise source computation and noise radiation. To achieve the real time simulation of ground noise, the noise radiation model database is built according to the idea of parametric modeling. The efficiency and effectiveness are illustrated through the comparison of computed results and compute time with the direct computing method for AH-1 rotor. Furthermore, the effects of atmosphere and ground on noise radiation as well as ground noise characteristics are analyzed with this method.
Key words: noise     rotor     secondary noise radiation     atmosphere absorption     ground absorption and reflection    

随着直升机的广泛应用,直升机噪声大的缺点越来越引起人们的重视。对于军用直升机,噪声大使其过早地暴露目标,与其隐身要求相悖;对于民用直升机,由于常在人口稠密的地区起飞着陆以及低空慢速飞行,噪声污染已成为公害之一。旋翼气动噪声是直升机的主要噪声源,国内外在此方面已开展了很多研究工作。然而这些研究主要集中在旋翼气动噪声分析方法建模[1-2]、发声机理[3-5]、声学特性以及设计降噪[6-7]等方面,对于直升机飞行过程噪声的预测及其地面声场特性的研究则开展较少。

一般而言,地面噪声预测分为噪声源计算、噪声传播及地面噪声合成等部分[8]。若采用直接计算法,需要计算旋翼旋转一周内桨叶上所有辐射点对地面观测点的噪声。当观测点数目较多时,计算将非常耗时。这对于实时分析直升机运动过程中的动态声场特性显然是不现实的。20世纪80年代,Boxwell等[6]开展了飞行参数对旋翼气动噪声的影响研究,并对旋翼噪声进行了参数化表达。基于此,Fleming[9]和Conner[10]等开展了基于试验测量的直升机噪声实时预测方法研究。首先由外场测量不同飞行状态下的地面声场数据;通过噪声逆向传播技术,将地面声场映射形成声辐射球;建立关于飞行参数的声辐射球库,当计算直升机某一飞行过程中地面声场时,根据状态参数插值获取相应的声辐射球,再将其传播至地面。根据上述研究,目前NASA已经发展形成了直升机噪声实时分析系统——Rotor Noise Model(RNM)[11]

虽然实飞测量可以获得更丰富直升机噪声数据,但是也存在技术难度大、代价高等缺点。为此,本文借鉴RNM思想,建立基于数值计算的直升机噪声地面声场实时预测方法:其中声辐射球采用自由尾迹/FW-H方程进行求解;引入大气、地面等环境因素对噪声传播特性的影响噪声传播模型。以AH-1/OLS试验旋翼为算例,分析验证了方法的有效性。

1 计算方法

图 1是本文建立的方法与直接地面噪声预测方法的对比。地面声场的计算包括旋翼气动/噪声预测、声传播(包含环境因素对噪声传播特性影响)、地面声场合成等。本文在①、②过程间加入了声辐射球,由二级声辐射球来传递旋翼的噪声信息,由此建立了高效的地面声场预测方法。

图 1 2种地面噪声预测方法对比 Fig. 1 Comparison between two ground noiseprediction methods
1.1 旋翼气动/噪声预测模型

旋翼噪声预测包括气动力和噪声计算两部分。本文中旋翼气动力采用课题组发展的自由尾迹方法进行求解[12]。旋翼尾迹模型由近尾迹和远尾迹两部分构成。近尾迹由桨叶环量沿径向变化引起的尾随涡和沿方位角方向变化引起的脱体涡组成。远尾迹仅包括单根桨尖涡线,由沿桨叶后缘脱出的尾迹在转过30°~60°方位角后卷起形成,如图 2所示。桨尖涡线的控制方程为

(1)
图 2 旋翼自由尾迹模型示意图 Fig. 2 Schematic of rotor free-wake model

式中:r(ψ,ζ)为控制点的位置矢量,ψ为桨叶方位角,ζ为尾迹寿命角;V(r(ψ,ζ))为控制点空间运动速度,包括自由流速度、附着涡及其他涡线引起的诱导速度。诱导速度采用毕奥-萨伐尔定理进行求解,涡核模型采用了Scully代数涡核模型[13],初始涡核半径为0.2c(c为桨叶弦长)。

为更好地模拟桨叶三维效应的影响,采用了Weissinger-L升力面理论进行了新生涡的求解,该方法的具体实现可参见文献[14]

在获得桨叶气动力基础上,可以进行旋翼气动噪声的计算。噪声计算采用基于FW-H方程的Farassat 1A公式[15]。Farassat 1A公式将噪声分为两部分,即

(2)

式中:p′p′tpl分别为总噪声、厚度噪声和载荷噪声;X为观测位置;t为观测时间。

1.2 二级声辐射模型

当噪声传播距离较远时,可将发声体看作点声源[16],从而简化计算。然而先前研究表明,旋翼噪声辐射具有明显的方向性,即不同方位观测点接收到的旋翼噪声大小差别较大[17]。因此,将旋翼噪声看作单纯的点声源将失去传播的方向性特点。借鉴点声源思想,本文建立了基于球声源的二级声辐射模型,如图 3所示。在该模型中,以

图 3 二级声辐射模型示意图 Fig. 3 Schematic of secondary noise radiation model

位于桨毂处的声辐射球球心代替实际桨叶声源作为噪声辐射中心,声辐射球上等角度间隔的观测点噪声代表了不同方位旋翼噪声辐射的实际声源。旋翼噪声沿从辐射中心到观测点的直线传播,声辐射点为该直线与声辐射球面的交点,噪声信息可通过二维插值方法在相邻4个观测点上插值得到。由此将随桨叶旋转的多个运动声源转换为声辐射球面上单一静态辐射点,提高了计算效率。

1.3 噪声传播模型

噪声传播过程中不仅有随着传播距离增加造成的扩散损失,还有大气吸声效应对其造成的衰减和地面对其的反射与阻抗作用[18]。对于距地面一定高度的观测点,其接受到的噪声有直达声波和地面反射声波两部分(见图 4)。地面观测点接收到的噪声声压可表示为

(3)
(4)

式中:p(x,t)为观测点接收到的噪声;pff(x,t)为直达声波;pgr(x,t)为地面反射声波;r1为直达声波传播距离;r2为反射声波传播距离;1r1(2)为扩散损失;x为观测点位置;α和Π(Φ1,ω)分别为大气吸声系数和地面吸声衰减系数,其表达式如下:

(5)
(6)
图 4 观测点噪声接收示意图 Fig. 4 Schematic of noise receiving at observation point

式中:σe为地面等效流阻。其他参数意义可参见文献[17]

需要指出的是,现有的声传播理论一般是建立在频域内,而2.1节中旋翼气动噪声的计算则是在时域内开展的。为此,本文采用了傅里叶变换和反变换进行噪声的时、频变换。

1.4 声场计算流程

实际飞行过程中,直升机的飞行姿态(飞行速度、加速度、航迹角等)是不断变化的,每个飞行姿态对应着一个声辐射球。研究表明,定常或小机动飞行状态下旋翼噪声可归结为参数前进比和桨盘倾角的函数[6]。因此,本文基于“参数化建模”思想,建立了以前进比和桨盘倾角为特征参数的声辐射球库。这样在地面噪声预测中,可根据所需状态的特征参数,直接选取相应的声辐射球用于后续噪声传播计算,而不用实时计算相应的声辐射球,从而进一步提高计算效率。

图 5给出了采用本文所建方法进行地面噪声预测的流程图。计算流程主要分为3个部分。首先,建立了以前进比(μ)和桨盘倾角(αTPP)为特征参数的声辐射球库;其次,根据所需计算状态的特征参数,选取相应的声辐射球,确定相应的声辐射点并插值得到相应的噪声信息;最后,进行噪声传播计算并进行相应的噪声后处理。

图 5 快速地面声场计算流程图 Fig. 5 Flowchart of quick ground noise computing
2 结果讨论与分析 2.1 计算细节

本节以AH-1/OLS试验旋翼[6]作为算例,对所建地面声场预测方法的有效性进行了分析。OLS试验旋翼是AH-1直升机旋翼的1/7缩比模型,共有2片桨叶,采用BHT-540对称翼型,直径为1.916 m,展弦比为9.22,具有-10°线性扭转,桨叶根切比为18.2%。算例状态参数为拉力系数0.005 4,前进比0.164,总距6.2°,纵横向周期变距分别为-2.4°、1.0°。图 6给出了算例分析中旋翼空间位置及地面观测点分布示意图,其中O点为坐标系原点,桨榖中心位于地面中心正上方50 m处。参考直升机适航标准中测量点分布,文中设置了B、C、D 3个地面观测点,坐标分别为:(100,40,0),(100,0,0)和(100,-40,0)。此外,图中还给出了半径分别为25R、50R(R为旋翼半径)的声辐射球;B1C1D1B2C2D2分别为BCD观测点噪声传播路径与声辐射球的交点,即噪声传播至25R、50R时,3个地面观测点所在位置。

图 6 旋翼空间位置及地面观测点分布示意图 Fig. 6 Schematic of rotor spatial position anddistribution of ground observation point
2.2 方法验证

图 7给出了该状态下旋翼的自由尾迹几何形状。从图 7中看出,本文方法较好地模拟了旋翼尾迹畸变及旋翼尾迹的运动。

图 7 旋翼桨尖涡尾迹几何形状 Fig. 7 Wake geometry of rotor blade tip vortex

图 8给出了采用2种计算方法得到的半径为50R的噪声声压级SPL(单位为dB)等高图,其中快速计算法采用了半径为25R的二级声辐射球。该噪声等高图是将球面空间投影到二维平面图得出的。投影方法采用了具有良好保角性(即投影图形与原图形相比没有角度变形)的兰勃特投影[19],这样有利于准确反映出噪声辐射的方向性。图中周向角度表示旋翼旋转的方位角,径向角度表示观测点仰角(仰角为负表示观测点位于旋翼下方)。

图 8 2种计算方法噪声声压级计算结果对比 Fig. 8 Comparison of computed noise SPLresults of two different methods

对比图 8(a)~图 8(c)可以看出,相比于25R处的噪声声压级50R处的噪声声压级降低了约6 dB,这与传播距离增加一倍声压级降低6 dB的噪声衰减结果[20]一致。这说明半径为25R的二级声辐射球较好地模拟了旋翼远场噪声及噪声传播模型较好地模拟了噪声的衰减。图 8中还可以看出,该状态下旋翼最大噪声位于前行侧桨盘平面下方30°附近,体现出较强的载荷噪声辐射特性,这与文献[12]分析基本一致。此外,对比同一传播距离下由2种计算方法计算得出的声压级等高图可以看出,2种计算方法在噪声幅值及方向性两方面基本一致。

旋翼噪声包含近场噪声及远场噪声两部分。近场噪声随传播距离增加衰减较快而不会传播较远的距离,计算中应尽可能选择较大半径的声辐射球,使近场噪声充分衰减掉。然而过大的声辐射球半径将引起较大的大气吸声计算误差,因此合理选择二级声辐射球是非常重要的。图 8(d)给出了由半径为20R的声辐射球计算得到的50R处的噪声声压级等高图。与图 8(b)对比可知,虽然该图较好地体现出旋翼噪声传播的方向性,但是最大声压级高出约2 dB,反映出该半径的声辐射球并没有充分衰减掉近场噪声。综合上述分析可以看出本文所采用二级声辐射球是合理的。

图 9给出了由2种计算方法计算得到B2C2D2 3个球面观测点的声压时间历程。可以看出,二者在幅值及相位两方面基本一致。综合图 8图 9的计算分析说明了本文快速地面声场预测方法的有效性。

图 9 地面观测点声压时间历程 Fig. 9 Sound pressure time history of groundobservation points

表 1给出了2种方法的计算时间对比。算例中桨叶声源辐射点共有18×28×201(桨叶分段数×弦向分段数×观测时间)。在直接计算法中,需要计算所有声源点对地面观测点的噪声辐射;而快速地面声场预测模型中建立了旋翼各方位角噪声辐射特性的二级声辐射球,地面声场预测时,仅需将相应的二级声辐射点噪声映射至地面观测点即可。当地面观测点数目均为10万时,前者需要耗时3.2 h,而后者仅为305 s。需要指出的是,快速地面声场计算模型虽然可以大幅度提高地面声场计算效率,但是当传播距离较近时,噪声中含有近场项及远场项,噪声衰减并不随传播距离的增加线性衰减,此时应采用直接计算法进行地面声场预测。

表 1 2种方法的声辐射点数及计算时间对比 Table 1 Comparisons of radiation points and computation time between two methods
计算方法 声辐射点数目 观测点数目 计算用时
二级声源法 1 10 万 305 s
直接计算法 101 304 10 万 3.2 h

2.3 大气吸声影响分析

本节进行了大气吸声效应对噪声传播特性的影响分析。图 10对比了有、无大气吸声效应情况下地面观测面的噪声声压级(单位为dB)等高图。算例中直升机前飞方向(v)沿XX轴正向,旋翼投影位于地面中心,如图 10(a)所示。可以看出,该状态下最大噪声位于地面右前方、最小噪声出现在地面中心附近。此外,由于大气吸声作用地面旋翼噪声降低了约6.5 dB,但大气吸声效应并没有改变噪声辐射的方向性。图 11为地面观测点B、C、D的声压时间历程。可以看出大气吸声对不同地面观测点作用处处相似,仅改变了噪声的幅值大小而没有改变噪声的相位特征。

图 10 大气吸声对噪声传播特性影响 Fig. 10 Effects of atmosphere absorption on noiseradiation characteristics
图 11 大气吸声对噪声声压时间历程的影响 Fig. 11 Effects of atmosphere absorption on noise soundpressure time history
2.4 地面效应影响分析

图 12进一步分析了地面吸声效应对旋翼噪声传播特性的影响,给出了3种常见地面(地面类型及地面声学特性描述参数如表 2[21]所示)的噪声声压级(单位为dB)等高图。对比图 12图 10可以看出,地面吸声效应除进一步降低噪声强度外,对噪声的传播方向几乎没有影响。此外对比不同地面声压级等高图可以看出,流阻对地面吸声效应具有较大的影响,流阻较大的地面噪声声压级较低。因而,对直升机而言,应尽量选择旧土路等流阻较大的地面飞行以降低直升机的可探测性。图 13C观测点为例,给出了不同地面类型对旋翼噪声的阻抗作用。可以看出,不同地面下观测点接受到噪声的声压历程变化曲线相位一致,而幅值不同;这说明地面的阻抗作用对噪声幅值影响较大,几乎不会影响噪声的相位特征。

图 12 不同地面类型对噪声传播特性影响 Fig. 12 Effects of different grounds on noiseradiation characteristics
表 2 不同地面声学特性参数[21] Table 2 Acoustic characteristic parameters ofdifferent grounds[21]
编号 地面类型 空隙率 流阻/
(kg·s-1·m-3)
地面1 间隙填以碎石的旧土路 0.400 7 500
地面2 裸露并被雨水填实的土地 0.350 17 100
地面3 被尘封并很少使用的沥青路面 0.250 120 000

实际的地面噪声接收点(如人耳等)距地面有一定的高度,因此其接受到的噪声不仅含有旋翼直接辐射的直达声波,还有经地面反射后接收到的反射声波(如图 4所示)。图 14给出了噪声观测点位于地面观测点C点上方1.0R、2.1R、3.0R时,观测点接收到的声波时域信号。算例中由于观测点距噪声源较远,其高度变化对反射声波入射角影响较小,因此实际计算中采用了直达声波作为反射声波进行计算分析。由图 14可以看出,不同高度观测点接收到的反射声波与总声波均发生了较明显的变化。反射声波,由于传播距离较大,扩散衰减和大气吸声效应进一步增强,因而声压幅值较小;此外传播距离的改变引起观测点接收到的反射声波相位发生改变。总声波则由于直达声波与反射声波的叠加作用,发生了较明显改变,且当观测点高度为2.1R时,直达声波与反射声波相位变化一致,发生“声共振”现象,使得总声波声压幅值明显增大。图 15给出了观测点接收到的总噪声频域图,可以看出不同高度观测点接收到的总噪声频域信号发生了较为明显的变化。当发生“声共振”时,总声波低频段噪声明显增强,而低频噪声具有较强的穿透性,非常不利于直升机的声隐身。

图 13 观测点噪声声压时间历程 Fig. 13 Noise sound pressure time history of observation points
图 14 不同高度观测点接收到的声压时间历程 Fig. 14 Accepted sound pressure time history ofobservation point at different heights
图 15 不同高度观测点总声波频域特性 Fig. 15 Frequency-domain characteristics of total acousticwave of observation point at different heights
3 结 论

本文开展了直升机飞行过程中高效地面声场预测方法的研究并建立了相应计算模型,应用所建方法计算分析了大气、地面等环境因素对噪声传播特性的影响。

1) 在旋翼噪声预测及噪声传播之间,结合点声源假设和旋翼噪声辐射的方向性特点建立了二级声辐射球,从而将随桨叶旋转的多个运动声源转换为声辐射球面上单一静态辐射点,显著提高了计算效率。当地面观测点数量为10万时,传统计算方法耗时约为3.2 h,本文所建方法仅为305 s。

2) 为实时预测直升机运动过程中旋翼噪声对地面声场影响,文中提出了“特征参数声辐射球库”。在地面声场预测中,可根据旋翼所处状态的特征参数,直接选取相应的二级声辐射球进行地面噪声辐射计算。

3) 大气、地面吸声效应对旋翼噪声具有较强的衰减作用,但几乎不会改变旋翼噪声传播的方向性;地面流阻对地面吸声效应具有较大的影响。

4) 对于距地面具有一定的高度的观测点,其接收到的噪声有直达声波与反射声波;地面反射对观测点时、频域噪声均具有较大的影响,且会发生“声共振”现象,此时低频噪声将会增大,对地面人员产生较大的影响。

5) 应当指出,本文所建快速地面声场预测方法对二级声辐射球半径具有一定的要求,当噪声传播距离较近时,应采用直接计算法进行地面声场预测。

参考文献
[1] BRENTNER K S, LYRINTZIS A S, KOUTSAVDIS E K. Comparison of computational aeroacoustic prediction methods for transonic rotor noise prediction[J]. Journal of Aircraft, 1997, 34 (4) : 531 –538. DOI:10.2514/2.2205
[2] BOYD D D JR.HART-Ⅱ acoustic predictions using a coupled CFD/CSD method[C]//Proceedings of 65th Annual Forum of the AHS,2009. http://cn.bing.com/academic/profile?id=2119794446&encoded=0&v=paper_preview&mkt=zh-cn
[3] DESOPPER A, LAFON P, PHILIPPE J J, et al. Effect of an anhedral sweptback tip on the performance of a helicopter rotor[J]. Vertica, 1988, 12 (4) : 345 –355.
[4] YU Y H. Rotor blade-vortex interaction noise[J]. Progress in Aerospace Sciences, 2000, 36 (2) : 97 –115. DOI:10.1016/S0376-0421(99)00012-3
[5] EDWARDS B,COX C,BOOTH E R JR.Revolutionary concepts for helicopter noise reduction-S.I.L.E.N.T program:NASA-CR-2002-211650[R].Washington,D.C.BOOTH E R J R.: NASA,2002.
[6] BOXWELL D A, SCHMITZ F H, SPLETTSTOESSER W R, et al. Helicopter model rotor-blade vortex interaction impulsive noise:Scability and parameteric variations[J]. Journal of the American Helicopter Society, 1987, 32 (1) : 3 –12.
[7] SAKOSWKY P C,CHARLES B D.Noise measurement test results of AH-1G operational loads survey:Bell Helicopter Co.Report 299-099-831[R].Fort Worth,TX:Bell Helicopter Co.1976.
[8] GOPALAN G, SCHMITZ F H, SIM B W. Flight path management and control methodology to reduce helicopter blade-vortex(BVI) noise[J]. Journal of Aircraft, 2002, 39 (2) : 193 –205. DOI:10.2514/2.2923
[9] FLEMING G G,RICKLEY E J.HNM heliport noise model,version 2.2,user's guide:Dot/FAA/EE/94-01,DOT-VNTSC-FAA-94-3[R].Washington,D.C.:Office of Environment and Energy,1994.
[10] CONNER D A,PAGE J A.A tool for low noise procedures design and community noise impact assessment:The rotorcraft noise model(RNM)[C]//Heli Japan 2002,2002:T213-6.
[11] CONNER D A,BURLEY C L.Flight acoustic testing and data acquisition for the rotor noise model(RNM)[C]//Proceedings of 62nd Annual Forum of AHS,2006.
[12] 吕维梁, 招启军, 徐国华. 计入畸变修正的旋翼尾迹前飞状态稳定性分析[J]. 航空学报, 2012, 33 (11) : 1958 –1966. LU W L, ZHAO Q J, XU G H. Analysis of rotor wake s stability in forward flight associated with distortion revise[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33 (11) : 1958 –1966. (in Chinese)
[13] LEISHMAN J G. Principles of helicopter aerodynamics[M]. 2nd ed Cambridge: Cambridge University Press, 2006 .
[14] WEISSINGER J.The lift distribution of swept-back wings:NASA-TM-1120[R].Washington,D.C.:NASA,1947.
[15] FARASSAT F, SUCCI G P. The prediction of helicopter rotor discrete frequency noise[J]. Vertica, 1983, 7 (4) : 309 –320.
[16] DOWLING A P, FFOWCS WILLIAMS J E, WRIGHT W M. Sound and sources of sound[J]. American Journal of Physics, 1985, 53 (6) : 320 .
[17] SIM B W,SCHMITZ F H.Acoustic phasing and amplification effects of single-rotor helicopter blade-vortex interactions[C]//Proceedings of 55th Annua Forum of AHS,1998.
[18] WEIR D,JUMPER S,BURLEY C L.Aircraft noise prediction program theoretical manual:Rotorcraft system noise prediction system:NASA-TM-83199[R].Washington,D.C.:NASA,1995.
[19] SNYDER J P.Map projections:A working manual:USGS Paper 1395[R].Washington,D.C.:U.S.Government Printing Office,1987.
[20] 马大猷, 沈儫. 声学手册[M]. 北京: 科学出版社, 2001 : 31 -38. MA D Y, SHEN H. Manual of acoustics[M]. Beijing: Science Press, 2001 : 31 -38. (in Chinese)
[21] ATTENBOROUGH K. Acoustical impendance models for outdoor ground surfaces[J]. Journal of Sound and Vibration, 1985, 99 (4) : 521 –544. DOI:10.1016/0022-460X(85)90538-3
http://dx.doi.org/10.13700/j.bh.1001-5965.2015.0756
北京航空航天大学主办。
0

文章信息

王菲, 史勇杰, 徐国华
WANG Fei, SHI Yongjie, XU Guohua
基于二级声辐射模型的地面声场高效预测方法
Efficient ground noise prediction method based on secondary noise radiation model
北京航空航天大学学报, 2016, 42(11): 2445-2453
Journal of Beijing University of Aeronautics and Astronsutics, 2016, 42(11): 2445-2453
http://dx.doi.org/10.13700/j.bh.1001-5965.2015.0756

文章历史

收稿日期: 2015-11-17
录用日期: 2016-02-29
网络出版时间: 2016-04-13

相关文章

工作空间