Ignition process and propellant grains distribution of the two-module charge
-
摘要: 模块装药点传火过程中药粒堆积形态对膛内起始压力波特性有重要影响,而模块装药点传火过程中药盒破裂后药粒飞散过程决定了药粒最终堆积形态。为此设计了模块装药可视化点传火模拟试验装置,通过高速摄像系统,观测不同初始装填位置的两模块装药点传火、药盒破裂及药粒散布过程。试验结果表明,两模块初始装填位置远离底火端且两药盒装填间距增大时,药室内传火时间变长,两个模块药盒燃烧更充分,模块盒的破裂面增大。点传火试验结束后,药室内模拟药粒散布在以底火侧端面中心为起点的轴向195~500 mm区域。其中,药粒主要分布于药室右侧陡坡状堆积区域。基于试验建立了模块装药点传火过程中药盒破裂后药粒散布的三维非稳态气固两相流模型,并进行了模拟计算。计算得到的最终药粒散布与试验测得结果基本吻合,验证了模型的合理性。Abstract: The accumulation form of propellant grains has a great effect on the initial chamber pressure wave in the ignition and flame-spreading process of a modular charge. In this process, the grains distribution is determined by the dynamic characteristics of grains after the cartridge is broken. Therefore, a visualized ignition simulation experimental device was designed for the ignition test of the two-module charge with different initial loading positions. A high-speed camera system was used to observe the ignition and flame propagation, the rupture of combustible cartridge cases, and the moving process of the propellant grains. The experimental results show as follows. When the two-module charging position is far from the primer and the spacing between the two modules is increased, the time of flame propagation in the chamber is prolonged. And the cartridge cases are more completely burned and their rupture areas become larger. The propellant grains in the chamber are finally scattered in the axial 195–500 mm area starting from the end face of the primer side. The grains are mainly distributed in the steep-slope accumulation on the right side of the chamber. On the basis of the experiment, a three-dimensional unsteady gas-solid two-phase flow model for the modular charge was established. The dynamic process and distribution of the propellant grains was simulated. The calculation results are basically consistent with the test ones, which validate the rationality of the established model.
-
防护结构是指能够抵抗预定杀伤武器破坏作用的工程结构, 目前常见的野战阵地防护结构大多由钢筋、混凝土材料筑成, 并设置了各类尺寸的孔口供射击、观察、通风和出入使用。当常规武器命中带孔防护结构后, 接触爆炸产生的威力强大的冲击波经孔口进入防护结构内部, 可能对结构内部的人员及设备造成一定程度的毁伤, 甚至导致人员伤亡、设备瘫痪, 使野战阵地丧失作战能力。而在现代战争中, 带孔防护结构往往是精确制导武器打击的重点, 其生存能力直接关系着战争的胜负。因此, 关于常规武器对带孔防护结构毁伤效应的研究对防护工程的建设、综合作战效率的提高具有重大意义。
目前, 已对自由场空爆、坑道内爆炸等条件下冲击波的传播规律进行了深入的研究, 提出了诸如Henrych公式[1]、Sadovsk公式[2]、TM5-855-1公式[3]等具有广泛影响力的计算公式。C.R.Welch[4]在模型坑道中进行化爆实验, 得到了堵口爆炸情况下坑道内空气冲击波冲量的经验公式; LS2000设计规范[5]中提出了堵口爆炸时坑道内空气冲击波冲量的确定方法; 庞伟宾等[6]对爆炸冲击波在直坑道、T型坑道内的传播规律及其影响因素进行了原型坑道爆炸实验研究; 黄广炎等[7]对炸药装药在腔室中心、墙角等不同约束情况下爆炸冲击波的传播特性进行了数值计算。
虽然对冲击波传播规律问题的研究已取得了一定成果, 但既有研究均有其特定的适用范围[8], 尚不能准确反映接触爆炸荷载作用下带孔防护结构内冲击波的传播规律。随着精确制导技术的快速发展, 战斗部直接命中防护结构并发生接触爆炸的可能性越来越大, 这就要求对带孔防护结构内冲击波的传播规律有清晰的认识, 以便准确评估爆炸冲击波对防护结构的毁伤效应。本文中, 采用ATUODYN软件构建有限元模型, 对接触爆炸荷载作用下带孔防护结构内冲击波的传播过程进行计算分析, 以期获得冲击波各特征参数的变化规律。
1. 计算参数标定
为了确认数值计算所采用软件的有效性、所使用材料状态方程及本构模型的适用性, 对某直坑道化爆实验[9]进行相同工况下的数值分析, 通过计算结果与实测结果的比对, 判断误差是否在容许范围内。模型坑道总长21m, 由21节可拼装式钢结构单元组成。每节坑道单元长1m, 净截面尺寸为600mm×600mm。爆室壁厚80mm, 其余单元壁厚20mm。坑道实物如图 1所示, 为保持模型坑道在爆炸实验过程中的气密性, 坑道单元通过高强螺栓连接, 单元与单元之间设有密封圈, 坑道两端开口。
直通道测点布置如图 2所示。直坑道口部化爆实验测点1(测点距炸药262cm)和测点4(测点距炸药562cm)的冲击波超压Δp计算波形和实测波形[9]的比较如图 3所示。
从直坑道化爆实验的情况来看, 典型监测点计算波形与实验波形基本吻合, 峰值压力、到达时间、正压作用时间误差均不超过20%。因此, 本文计算所采用的方法和材料参数及求解设置均为合理的。
2. 计算概述
为了便于计算, 认为TNT炸药位于带孔防护结构中间孔口的正上方, 利用对称性, 建立1/2模型, 如图 4所示。带孔防护结构呈筒状, 组成材料为混凝土和钢筋, 混凝土标号为C35;结构高3 100mm, 壁厚d=500mm, 内部半径r=2 000mm; 结构正面开有3个孔口, 尺寸均为500mm×300mm, 背面开有一尺寸为500mm×1 500mm的出入口; 钢筋采用HRB335钢筋, 直径为12mm, 分3层布置, 如图 5所示。模型中, 混凝土采用Lagrange网格进行建模, 网格尺寸为50mm, 底部固接; 钢筋采用beam单元, 采用join方法将钢筋与钢筋、钢筋与混凝土之间的节点进行连接; 空气采用Euler算法, 网格尺寸为25mm, 边界面设为物质流出。由于带孔防护结构形状较复杂, 在AUTODYN中建立模型难度大, 对于该问题采用如下方法加以解决:在ANSYS主程序中建立带孔防护结构模型并对其进行网格划分, 将所建模型保存成.k文件的形式, 并导入AUTODYN。
炸药爆炸的反应区压力模型采用高能炸药燃烧函数因子模型, 状态方程采用JWL状态方程[10], 炸药爆轰过程中压力和比容的关系为:
p=A(1−ωR1V)e−R1V+B(1−ωR2V)e−R2V+ωe0V 式中:A、B、R1、R2、ω为材料常数, p为压力, V为相对体积, e0为初始比内能。
计算中使用的TNT材料参数为:密度ρ=1 630kg/m3, 爆速D=6 930m/s, A=381.2GPa, B=3.822GPa, R1=4.15, R2=0.9, ω=0.35, e0=6.12GJ/m3, pCJ=21.4GPa。
空气采用理想气体状态方程[11], 初始密度为1.225kg/m3, 比热比为1.4;混凝土采用修正后的RHT本构模型[12]; 钢筋采用Johnson-Cook热黏塑性材料模型[13]。
TNT炸药位于结构正面中间孔口的正上方, 到孔口上边缘的距离为x, 为了研究爆炸位置对空气冲击波毁伤效应的影响, 分别取x=0.1, 0.2, 0.4, 0.8m。TNT炸药质量为m, 为了研究炸药当量对毁伤效应的影响, 分别取m=1, 2, 4, 8kg。因此, 模拟的工况共有4×4=16个。
结构内布置2个观测点, 测点1位于结构内部中心处, 测点2位于孔口处, 如图 4所示。
3. 带孔防护结构内冲击波传播规律分析
3.1 冲击波传播过程分析
分别对不同药量、不同爆炸位置的16个工况进行计算, 得到了结构内部中心处(测点1)和孔口处(测点2)爆炸冲击波压力p随时间的变化过程。根据Δp=pt-p0即可得到对应的爆炸冲击波超压。
虽然在不同工况下炸药的质量和位置有所区别, 但测点处爆炸冲击波超压-时间曲线的大致线形都较相似。以x=0.4m, m=4kg的工况为例, 孔口处和结构内部中心处冲击波超压-时间曲线如图 6所示。爆炸冲击波超压Δp最初为零, 此时结构外TNT炸药虽然已经发生爆炸(炸药爆炸时为t=0时刻), 但由于冲击波尚未传播至测点处, 故孔口处和结构内部中心处的空气压力仍与大气压力p0相同。当t=0.616ms时, 爆炸冲击波波阵面到达孔口处, 测点2处压力急剧上升, 并在t=0.716ms时达到正压峰值, 此时孔口处峰值超压Δpmax=627.0kPa, 如图 6(a)所示。随后, t=2.138ms时, 冲击波到达结构内部中心处, 此处超压由最初的零迅速达到正超压峰值Δpmax=30.7kPa, 如图 6(b)所示。在爆炸冲击波波阵面经过结构内部中心处后, 测点1的超压值并没有迅速地稳定为零, 而是快速减小, 直至出现了负压。随后, 结构内部中心处出现了较长时间的超压震荡, 并且在超压震荡期间, 正压和负压的绝对值均较大, 如图 6(b)所示。
结合图 7所示的冲击波传播压力云图, 经分析可知, 造成上述超压震荡现象的原因是:爆炸冲击波波阵面通过孔口进入结构内部(如图 7(c)所示)后, 迅速到达测点1处(如图 7(d)所示), 此时结构内部中心处出现峰值超压, 随后, 此处压力由正压逐渐变为负压; 冲击波继续向前运动, 并与结构内壁发生碰撞, 产生的反射冲击波再次到达结构内部中心(如图 7(e)所示), 导致测点1处再次出现正超压; 之后, 冲击波与结构内壁发生了多次碰撞, 反射冲击波多次经过结构内部中心, 造成此处超压出现震荡, 并且正压作用时间也大大增加。
当冲击波在结构内壁产生多次反射后, 流场变得十分紊乱, 呈现复杂波特性, 部分位置压力将会增强, 如图 7(f)所示。
由图 6可知, 在相同的工况(x=0.4m, m=4kg)下, 孔口处(测点2)爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻显然要先于结构内部中心处(测点1), 并且孔口处的峰值超压远大于结构内部中心处(627.0kPa≫30.7kPa)。在爆炸冲击波波阵面第1次经过孔口处后, 测点2处的超压值虽然也会发生震荡, 但震荡幅值相对于超压峰值要小得多。
在研究空气冲击波毁伤判据时, 通常会考虑描述空气冲击波强弱的3个主要参数:峰值超压、正压区作用时间和正压冲量。冲击波峰值超压表示冲击波瞬间作用的量, 而正压冲量表示在正压区时间范围内超压的持续作用量, 在一定条件下两者均可作为目标毁伤程度的判定依据。在已获得冲击波超压-时间曲线的基础上, 通过对超压Δp对正压作用时间t+进行积分就可求得正压冲量I。由于图中冲击波超压在达到峰值后会出现一定时间的震荡, 同时考虑到第1股爆炸冲击波对目标的毁伤最严重, 所以在计算正压冲量时仅对第1段正压区进行积分。不同工况时带孔防护结构中心处冲击波的峰值超压和正压冲量如表 1所示。
表 1 不同工况下,结构内部中心处(测点1)冲击波的峰值超压和正压冲量Table 1. The maximum overpressures and the positive pressure impulse of the shock waves in the center of the structure(the gaging point 1)工况 x/m m/kg Δpmax/kPa I/(kPa·ms) 1 0.1 1 37.335 48.986 2 0.1 2 74.803 75.295 3 0.1 4 143.524 145.473 4 0.1 8 236.057 294.815 5 0.2 1 19.280 25.333 6 0.2 2 32.664 42.250 7 0.2 4 72.159 66.803 8 0.2 8 144.075 133.857 9 0.4 1 13.827 19.713 10 0.4 2 19.716 26.836 11 0.4 4 30.745 43.958 12 0.4 8 50.812 63.402 13 0.8 1 11.581 19.290 14 0.8 2 18.703 25.821 15 0.8 4 23.837 41.739 16 0.8 8 35.539 49.389 在爆炸位置固定的前提下, 比较不同TNT药量情况下的爆炸冲击波超压-时间曲线。如图 8(a)所示, 在工况1(x=0.1m, m=1kg)和工况2(x=0.1m, m=2kg)下, 炸药的爆炸位置相同, 均为x=0.1m, 工况2的药量是工况1的2倍, 两者的冲击波超压-时间曲线线形相似, 但工况2的峰值超压Δpmax, 2=74.803kPa高于工况1的峰值超压Δpmax, 1=37.335kPa, 并且工况2的爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻均要先于工况1。由此可知, 爆炸位置固定的情况下, 带孔防护结构内部中心处的峰值超压随着药量的增加而上升, 爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着药量的增加而提前。
在相同药量的前提下, 研究TNT炸药爆炸位置对冲击波传播规律的影响。如图 8(b)所示, 在工况3(x=0.1m, m=4kg)和工况7(x=0.2m, m=4kg)下, 炸药的药量相同, 均为m=4kg, 工况7的炸药到正面中间孔口上边缘的距离x是工况3的2倍, 两者的冲击波超压-时间曲线线形相似, 但工况7的峰值超压Δpmax, 7=72.159kPa低于工况3的峰值超压Δpmax, 3=143.524kPa, 并且工况7的爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻均要晚于工况2。由此可知, 在相同药量的前提下, 带孔防护结构内部中心处的峰值超压随着炸药至孔口距离的增加而下降, 爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着炸药至孔口距离的增大而推迟。
3.2 冲击波特征参数计算公式
为了建立接触爆炸荷载作用下带孔防护结构内部中心处冲击波特征参数的计算公式, 用量纲理论分析冲击波超压、正压冲量与时间的一般函数关系。
图 9为弹药在混凝土结构表面爆炸的示意图, 设装药爆炸能为Q, 混凝土结构的主要几何参数分别为正面孔口面积S、孔口下边缘至结构内地面的距离h、背面入口面积Sm、内部半径r、净高H和壁厚d。混凝土结构内部体积Vc=πr2 H, 可由r和H表示, 不作为主定参量。爆心到混凝土结构孔口上边缘的距离为x, 测点到混凝土结构墙壁内侧面的距离为l, 周围大气压力和密度分别为p0和ρ0, 测点处的峰值超压为Δpmax, 从爆炸时刻算起, 正压作用时间为t+, 正压冲量为I。
于是该问题的主定参量组为:Q, S, Sm, r, h, H, d, x, l, p0, ρ0; 待定参量组为:Δpmax, t+, I。首先确定测点处的峰值超压Δpmax, 如果用f表示爆炸冲击波参数关系, 那么有以下函数:
f(Δpmax,Q,S,Sm,r,h,H,D,x,l,p0,ρ0)=0 (1) 采用L-M-T度量单位系统作为基本量纲, 构建量纲矩阵:
(2) 并计算:
RankA=3 (3) 解齐次方程:
Ay=0 该齐次方程有12-3=9个基本解, 可取为:
{y1=(1000000000−10)Ty2=(0−1/3000000101/30)Ty3=(0010−20000000)Ty4=(0001−20000000)Ty5=(0000−10000100)Ty6=(00001−1000000)Ty7=(000010−100000)Ty8=(0000−10010000)Ty9=(00001000−1000)T 式(5)给出了9个相互独立的量纲一量
视p0、ρ0为常数, 为简化工程计算, 令其等于1, 则
将装药爆炸能Q仅视为与药量m有关时
有函数:
f(π1,π2,π3,π4,π5,π6,π7,π8,π9)=0 (6) 峰值超压Δpmax的显式表达式为:
Δpmax=f1(π2,π3,π4,π5,π6,π7,π8,π9) (7) 即:
Δpmax=f1(xm1/3,Sr2,Smr2,lr,rh,rH,Dr,rx) (8) 式中:f1表示一个未定函数。
由上面的计算结果可以得知峰值超压Δpmax的主要影响因素为比例距
, 带孔防护结构的几何尺寸以及结构内径与爆心到孔口上边缘的距离之比
运用同样的方法, 确定正压冲量:
Im1/3=f2(xm1/3,Sr2,Smr2,lr,rh,rH,Dr,rx) (9) 3.3 数据拟合
在用实验方法整理函数f(·)的形式时, 一般能在自变量的某一范围内采用指数关系式[14], 即π=k0π1k1 ·π2k2 ·...·πnkn。为了得到最终的拟合公式, 将式(8)、(9)采用指数关系式表示, 即:
Δpmax=α0(xm1/3)α1(Sr2)α2(Smr2)α3(lr)α4(rh)α5(rH)α6(Dr)α7(rx)α8 (10) Im1/3=β0(xm1/3)β1(Sr2)β2(Smr2)β3(lr)β4(rh)β5(rH)β6(Dr)β7(rx)β8 (11) 对数值计算数据进行多元线性回归, 得到式(10)、(11)中的待定参数αi, βi(i=0, 1, …, 8), 即可最终确定各待定函数f(·)的表达式, 从而可以分别求得峰值超压Δpmax、正压冲量I的计算公式。
式(10)、(11)中待定变量过多, 若在其基础上总结经验公式是比较困难的。对于本文所研究的防护结构, 其几何尺寸、开孔大小都是确定的, 此处可只保留结构内部半径r以反映结构尺寸对爆炸冲击波的影响。这样, 函数拟合的难度就大大降低。结合数值计算所得结果(表 1), 针对本文所研究的带孔防护结构, 结构内部中心处爆炸冲击波特征参数的拟合公式为:
Δpmax=19.2(xm1/3)−0.76(rm1/3)−1.42 (12) Im1/3=15.6(xm1/3)−0.59(rm1/3)−0.31 (13) 与数值计算结果相比, 由拟合公式计算所得的峰值超压和正压冲量的最大误差为25%(小于30%), 满足冲击波特征参数预估的精度要求。
4. 结论
(1) 当爆炸冲击波经孔口进入结构内部后, 结构内部中心处的超压会快速到达正超压的峰值, 并且该峰值大多出现在第1个波峰处。由于冲击波与结构内壁的反复碰撞, 在峰值超压出现后, 结构内部中心处的超压值会出现较长时间、较大幅值的震荡。
(2) 孔口处爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻先于结构内部中心处, 并且孔口处峰值超压约为结构内部中心处的20倍。
(3) 爆炸位置固定的情况下, 带孔防护结构内部中心处的峰值超压随着药量的增加而上升, 爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着药量的增加而提前; 在相同药量的前提下, 带孔防护结构内部中心处的峰值超压随着炸药至孔口距离的增加而下降, 爆炸冲击波波阵面的到达时刻和峰值超压所对应的时刻随着炸药至孔口距离的增大而后移。
(4) 利用数值计算的结果, 针对带孔防护结构, 拟合得到了结构内部中心处爆炸冲击波特征参数的预估公式。与数值计算结果相比, 由拟合公式计算所得的峰值超压和正压冲量的最大误差为25%。
-
表 1 模块药盒组合装填参数
Table 1. Filling parameters of the modular cartridge
工况 X1/mm L/mm m1/g m2/g m3/g 1 40 10 14.2 20 600 2 60 50 14.2 20 600 -
[1] 王泽山. 模块装药技术及其进展 [J]. 含能材料, 2004, 12(A01): 122. DOI: 10.3969/j.issn.1006-9941.2004.z1.031.WANG Z S. Modular charge technology and its progress [J]. Chinese Journal of Energetic Materials, 2004, 12(A01): 122. DOI: 10.3969/j.issn.1006-9941.2004.z1.031. [2] 余斌. 刚性组合装药可燃传火管性能试验研究 [J]. 火炸药学报, 2002, 25(4): 70–71. DOI: 10.3969/j.issn.1007-7812.2002.04.026.YU B. Experimental study on ignition and flame spreading of modular charge [J]. Chineses Journal of Explosives and Propellant, 2002, 25(4): 70–71. DOI: 10.3969/j.issn.1007-7812.2002.04.026. [3] 周瑶. 模块装药总体结构及装药弹道性能研究[D]. 南京: 南京理工大学, 2005. DOI: 10.7666/d.y1000320. [4] 韩博, 张晓志, 邢浴仁, 等. 大口径火炮发射装药点传火模拟试验装置的研究 [J]. 兵工学报, 2008, 29(3): 262–265. DOI: 10.3321/j.issn:1000-1093.2008.03.002.HAN B, ZHANG X Z, XING Y R, et al. Research on simulation testing device of the ignition system for large caliber propellant charging [J]. Acta Armamentarii, 2008, 29(3): 262–265. DOI: 10.3321/j.issn:1000-1093.2008.03.002. [5] 张洪林. 模块装药性能研究[D]. 南京: 南京理工大学, 2009. DOI: 10.7666/d.y1698703.ZHANG H L. The study of performance for modular propelling charge [D]. Nanjing: Nanjing University of Science and Technology, 2009. DOI: 10.7666/d.y1698703. [6] 陆中兵, 周彦煌. 模块装药膛内两相燃烧模型及压力波模拟 [J]. 爆炸与冲击, 1999, 19(3): 269–273.LU Z B, ZHOU Y H. Two-phase-combustion mode and numerical simulation of pressure wave in the gun with modular charges [J]. Explosion and Shock Waves, 1999, 19(3): 269–273. [7] 赵毅. 模块装药点传火过程的数值模拟 [J]. 火炸药学报, 2003, 26(2): 33–36. DOI: 10.3969/j.issn.1007-7812.2003.02.012.ZHAO Y. Numerical simulation for ignition and flame-spreading of modular charge [J]. Chineses Journal of Explosives and Propellant, 2003, 26(2): 33–36. DOI: 10.3969/j.issn.1007-7812.2003.02.012. [8] 王育维, 郭映华, 董彦诚, 等. 可燃容器对小号模块装药压力波影响的研究 [J]. 火炮发射与控制学报, 2016, 37(2): 31–35; 45. DOI: 10.19323/j.issn.1673-6524.2016.02.007.WANG Y W, GUO Y H, DONG Y C, et al. Study of combustible case effects on pressure waves for low-zone of bi-modular charge [J]. Journal of Gun Launch and Control, 2016, 37(2): 31–35; 45. DOI: 10.19323/j.issn.1673-6524.2016.02.007. [9] DONG Y C, WANG Y W, GUO Y H, et al. Effect of random factors of charge on interior ballistic performance [C] // 31st International Symposium on Ballistics. USA: DEStech Publication, Inc., 2019. DOI: 10.12783/ballistics2019/33088. [10] 曾思敏. 内弹道中的脉冲X射线摄影技术装备的研制及其应用的研究[D]. 南京: 南京理工大学, 1988. DOI: 10.7666/d.Y139890. [11] SOPER W G. Grain velocities during ignition of gun propellant [J]. Combustion and Flame, 1975, 24(2): 199–202. DOI: 10.1016/0010-2180(75)90148-0. [12] SOPER W G. Ignition waves in gun chambers [J]. Combustion and Flame, 1973, 20(2): 157–162. DOI: 10.1016/S0010-2180(73)80169-5. [13] EAST J L. Ignition and flame spreading phenomena in granular propellant gun charges [J]. Progress in Astronautics and Aeronautics, 1979, 66: 228–245. [14] 金志明. 火炮膛内压力波产生机理及其特征分析 [J]. 南京理工大学学报(自然科学版), 1992(1): 26–31. DOI: 10.14177/j.cnki.32-1397n.1992.01.006.JIN Z M. Generation mechanism and characteristic analysis of pressure wave in gun [J]. Journal of Nanjing University of Science and Technology, 1992(1): 26–31. DOI: 10.14177/j.cnki.32-1397n.1992.01.006. [15] 宋明, 金志明. 火炮膛内压力波的数值模拟 [J]. 兵工学报, 1992, 13(1): 12–19.SONG M, JIN Z M. Numerical simulation of pressure waves in guns [J]. Acta Armamentarii, 1992, 13(1): 12–19. [16] 王国强, 郝万军, 王继新. 离散单元法及其在EDEM上的实践[M]. 西安: 西北工业大学出版社, 2010: 16–18. 期刊类型引用(1)
1. 郑英,邹德波,赵铮. 爆炸切割冲击波防护仿真研究. 兵器装备工程学报. 2021(10): 243-248 . 百度学术
其他类型引用(6)
-