网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

超声速气流中复合材料壁板的动力学分析及其控制

  • 周凯
  • 倪臻
  • 黄修长
  • 华宏星
上海交通大学 高新船舶与深海开发装备协同创新中心,机械与系统国家重点实验室,振动冲击噪声研究所,上海,200240

最近更新:2021-01-05

DOI:10.6052/1672-6553-2020-080

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
目录contents

摘要

超声速飞行器外表具有大量的蒙皮壁板结构,其在巡航过程中会受到气动、热、声及机械等载荷的联合作用.本文针对超声速气流中的复合材料壁板结构,基于一阶剪切变形理论(First-Order Shear Deformation Theory,FSDT)和超声速活塞理论,推导得到了复合材料壁板的气热弹动力学的能量泛函,并运用Hamilton原理变分求得系统的控制方程,利用牛顿迭代法结合Newmark法,求解获得了壁板的临界颤振动压和时域动力学响应.通过变参数计算,分析了不同参数对于壁板动力学响应的影响.最后,应用非线性能量阱(Nonlinear Energy Sink,NES)对复合材料壁板的动力学响应进行颤振控制,结果表明,NES可以有效降低壁板的颤振极限环振幅,从而极大提高超声速飞行器的可靠性和寿命.

前言

超声速飞行器因其巡航速度快、突防性能好的特点,成为各大国的重点发展对象.蒙皮壁板作为超声速飞行器上重要的维形结构,其在巡航过程中会受到气动力、热载荷、声载荷和机械载荷的作

1,而这些载荷会显著影响壁板的动力学行为.因此,大量的学者们对超声速气流中的壁板动力学问题进行了深入的研究.

学者们在对超声速气流中的壁板进行动力学分析的过程中发展了大量的计算方法.对于一些简单的工况,壁板颤振问题的解析解可以得到.但是由于工程实际问题的复杂性,数值计算方法被应用得更多,如里茨法、伽辽金法、微分求积法、有限元法

2-5,它们被广泛地运用在各向同性壁板、正交各向异性壁板、功能梯度材料壁板、粘弹性壁板等结构的颤振分析.与此同时,研究者们进行了大量的参数分析,揭示了不同参数对于壁板的颤振特性影响规律.除此之外,壁板的颤振控制问题也得到了大量学者们的研究,并取得了一定的成6-10.

值得一提的是,尽管关于壁板颤振的文献较多,但是它们大多是基于经典边界下展开研究的.事实上,由于摩擦和蠕变的存

11,壁板的边界约束往往采用弹性边界更加符合实际.并且由于安装空间和经济条件的约束,传统的控制策略往往难以在工程中得到应用,因此,亟需发展新的壁板颤振控制策略.本文将建立一般边界下超声速气流中复合材料壁板的动力学模型,并采用简易便携的非线性能量阱(Nonlinear energy sink, NES)对壁板进行颤振控制.

1 模型描述

本文研究对象为复合材料层合壁板,如图1所示,其中复合材料壁板的尺寸为L1×L2×h.笛卡尔坐标系建立在复合材料层合壁板的中面上,假设超声速气流偏航角为θair.ZkZk+1分别为第k子层的下表面和上表面在z方向上的坐标.为了建立一般边界下复合材料壁板的动力学模型,根据改进傅里叶

12,五组人工虚拟弹簧将均布于复合材料层合壁板的边界来模拟壁板的约束条件.基于一阶剪切变形理论(First-Order Shear Deformation Theory,FSDT),三组位移弹簧(kukvkw)和两组扭转弹簧(kxky)将均布于壁板的四条边.

图1 超声速气流中的复合材料壁板

Fig.1 The composite laminated plate in supersonic airflow

图1所示的壁板为对象,基于FSDT,该壁板的位移可表示

13

U(x,y,z,t)=u(x,y,t)+zφx(x,y,t)V(x,y,z,t)=v(x,y,t)+zφy(x,y,t)W(x,y,z,t)=w(x,y,t) (1)

其中,uvw分别为壁板的中面在xyz方向的位移;φxφy为结构中面关于yx轴的转角;t为时间变量.

根据Von-Karman结构大变形理论,可以进一步得壁板应

13

εx=ux+zφxx+12(wx)2εy=vy+zφyy+12(wy)2γxy=(uy+vx)+z(φxy+φyx)+(wxwy)γxz=wx+φxγyz=wy+φy (2)

基于热弹性理论,可得壁板的应力如

14

σxkσykτyzkτxzkτxyk=Q11k¯Q12k¯00Q16k¯Q12k¯Q22k¯00Q26k¯00Q44k¯Q45k¯000Q45k¯Q55k¯0Q16k¯Q26k¯00Q66k¯εxk-α11k¯ΔTεyk-α22k¯ΔTγyzkγxzkγxyk-α12k¯ΔT (3)

其中,σxkσyk分别是第k子层xy方向的应力;τyzkτxzkτxyk为相应的剪应力;α11k¯α22k¯α12k¯为线膨胀系数;ΔT=T1-T0为温度变化量,其中T1为现在温度,T0为参考温度.

k子层的热应力可表示为:

σxTkσyTkτxyTk=-Q11k¯Q12k¯Q16k¯Q21k¯Q22k¯Q26k¯Q16k¯Q26k¯Q66k¯α11k¯ΔTα22k¯ΔTα12k¯ΔT (4)

k子层相应的热应变可表示

14

dxx=(ux)2+(vx)2+(wx)2+z2(φxx)2+z2(φyx)2dyy=(uy)2+(vy)2+(wy)2+z2(φxy)2+z2(φyy)2dxy=(ux)(uy)+(vx)(vy)+(wx)(wy)+z2(φxx)(φxy)+z2(φyx)(φyy) (5)

根据超声速活塞理论,气动力引起壁板的压差

14

Δp=-ρU2M2-1M2-2M2-11Uwt+wxcosθair+wysinθair (6)

其中,ρUM为来流的密度、速度和马赫数.

为了方便起见,定义无量纲气动力λ=ρU2L13/D11M2-1.其中,定义参考弯曲刚度D11=E11h3/[12(1-ν121ν211)]E11为壁板的第1子层的弹性模量.对于M1的工况,可以近似得到[(M2-2)/(M2-1)]μ/M2-1μ/M表达式.其中,μ=ρL1/ρh.

壁板的应变能可表示为:

Ue=12Ωk=1nzkzk+1σxk(εxk-α11k¯ΔT)+σyk(εyk-α22k¯ΔT)+τxyk(γxyk-α12k¯ΔT)+Ksτxzkγxzk+Ksτyzkγyzk+σxTkdxx+σyTkdyy+2τxyTkdxydzdxdy (7)

储存在虚拟弹簧里的势能为:

UBC=120L1kuy0(u)2+kvy0(v)2+kwy0(w)2+kxy0(φx)2+kyy0(φy)2y=0dx+120L1kuyL2(u)2+kvyL2(v)2+kwyL2(w)2+kxyL2(φx)2+kyyL2(φy)2y=L2dx+120L2kux0(u)2+kvx0(v)2+kwx0(w)2+kxx0(φx)2+kyx0(φy)2x=0dy+120L2kuxL1(u)2+kvxL1(v)2+kwxL1(w)2+kxxL1(φx)2+kyxL1(φy)2x=L1dy (8)

复合材料层合壁板的动能为:

Tk=120L10L2I0(ut)2+(vt)2+(wt)2+2I1(ut)(φxt)+2I1(vt)(φyt)+I2(φxt)2+(φyt)2dxdy (9)

其中,I0=k=1nρk(zk+1-zk)I1=12k=1nρk(zk+12-zk2)I2=13k=1nρk(zk+13-zk3).

采用有限带宽的高斯白噪声来模拟壁板受到的随机声载荷,其截断后的互谱密度函

7为:

Sp=S0=p0210SPL/10,0ffu,0,, (10)

其中,p0=20μPa为参考声压,SPL为声压分贝数,fu为截止频率.高斯白噪声的功率可通过其互谱密度函数和截止频率表

7为:

PW=S0fu (11)

因此,复合材料壁板所受的随机声载

7为:

p(t)=PWrandn([r,1]) (12)

气动力和声载荷对壁板做功之和为:

Wa=Ω(Δp+p)wdxdy (13)

将上述能量表达式带入到Hamilton原理中:

δt0t1(Tk-Ue-UBC+Wa)dt=0 (14)

变分可得系统的控制方程如下:

(K+KNL(t))Q(t)+CQ˙(t)+MQ¨(t)=F(t) (15)

NES有着质量轻、空间小、易于安装等特性,同时由于NES中非线性弹簧的存在,不仅使得可控频段大大加宽,也可以使得被控结构的能量实现靶向转移,从而控制性能远优于传统的线性吸振

15.如图2所示为背风面悬挂NES的复合材料层合壁板的示意图,这里假设NES悬挂于壁板背风面的(xnesynes)点.NES的设计参数中非线性刚度为knes,线性阻尼为cnes,小质量为mnes.

图2 超声速气流中含NES的复合材料壁板

Fig.2 The composite plate with NES in supersonic airflow

进一步建立含NES复合材料壁板的能量泛函.复合材料壁板的位移表达式、应力应变关系等与上述的一致,这里不再详细列出,重点给出NES的能量泛函.由于NES的存在,使得系统较上述复合材料壁板非线性动力学模型多出一个自由度,这里假设NES的位移为s,假设非线性弹簧的位移-力的关系为立方函数,则积分可得储存在NES非线性弹簧中的弹性势能

1617

Unes=14knes[w(xnes,ynes)-s]4 (16)

式中,w(xnes,ynes)为悬挂NES点处复合材料壁板的位移.

NES的阻尼项耗散能表达式为:

Wnes=-12cnes[w˙(xnes,ynes)-s˙]2 (17)

NES中质量块的动能表达式为:

Tnes=12mnes(s˙)2 (18)

将复合材料壁板和NES的能量表达式代入到Hamilton变分式中:

δt0t1(Tk+Tnes-Ue-Unes-UBC+Wa+Wnes)dt=0 (19)

进一步可写出如下的矩阵形式,即为含NES的复合材料壁板动力学控制方程:

(K+KNL(t))Q(t)+CQ˙(t)+MQ¨(t)=F(t) (20)

其中,上述控制方程(15)和(20)均可运用Newmark法联合牛顿迭代法求解.

2 验证与讨论

2.1 算例验证

为了验证本文模型的准确性,计算尺寸为0.305m×0.305m×0.00127m壁板在初始扰动下的动力学响应,该壁板的弹性模量为70GPa,泊松比为0.3,密度为2700kg/m3,边界约束为四边简支.如图3所示为本模型计算所得壁板(0.75L1,0.5L2)点处的极限环幅值与文献[

5]的对比(w˜=w/h),两条曲线吻合良好,可见本模型预测壁板气弹动力学响应的准确性.

图3 颤振极限环对比

Fig. 3 The comparison of the obtained LCO and that from literature

2.2 参数讨论

进一步进行参数讨论,以明确各参数对于壁板颤振特性的影响规律.这里选取的复合材料壁板尺寸为1m×1m,壁板共有五层铺层,每层的材料和厚度相同,为0.002m,壁板的材料特性如表1所示.

表1 复合材料层合板子层的材料特性
Table 1 Material property parameters of each single layer
E1(GPa)E2(GPa)G12(GPa)G13(GPa)G23(GPa)ρ(Kg/m3)v12α11(K-1)α22(K-1)α12(K-1)
175 32 12 12 5.7 1760 0.25 1.2e-6 2.3e-6 0

对于四边简支[0°/0°/0°/0°/0°]复合材料壁板,计算得其临界颤振动压为374.为了对比壁板在颤振前后的动力学响应情况,这里分别计算了气动压力为0,200和450下,壁板受到初始单位脉冲激励下的动力学响应,分别如图4-5和图6所示.图4(a)为壁板(0.75L1,0.5L2)点处的时域动力学响应,由于此处不考虑壁板自身的阻尼,故其在受到外界脉冲激励后,动力学响应不会衰减;图4(b)为该点时域动力学响应的相图,可以清晰看到壁板的响应没有收敛;图4(c)为该点响应的频谱图;图4(d)为壁板响应最大时刻的位移云图,可见在没有气动力时,壁板中点处的位移响应最大.进一步观察λ=200时,壁板(0.75L1,0.5L2)点处的动力学响应.由图5(a)可见,当气动压力小于临界颤振动压时,壁板在受到外界扰动后,由于系统气动阻尼的存在,最终动力学响应会收敛;绘制响应最大时刻壁板的位移云图如图5(d)所示,可知气动压力的存在会显著改变壁板最大位移点的位置.最后计算λ=450时,壁板(0.75L1,0.5L2)点处的动力学响应.如图6所示,当气动力大于临界颤振动压时,壁板在受到外界扰动后,壁板的响应不再会收敛,而是作等幅极限环运动;在图6(b)相图中可以清晰地看见封闭的环状图,说明了该运动为等幅极限环运动;图6(d)为壁板各点极限环振幅的云图,此时壁板最大变形点出现在(0.75L1,0.5L2).

图4 无气动力时壁板在脉冲激励下的响应(λ=200):(a)时域响应,(b)相图,(c)频谱图,(d)最大位移云图

Fig.4 The dynamic response of the composite plate with pulse excitation (λ=0): (a) time domain response, (b) phase plot, (c) frequency response, (d) the contour of the maximum displacement of the plate

图5 有气动力时壁板在脉冲激励下的响应(λ=200):(a)时域响应,(b)相图,(c)频谱图,(d)最大位移云图

Fig.5 The dynamic response of the composite plate with pulse excitation (λ=200): (a) time domain response, (b) phase plot, (c) frequency response, (d) the contour of the maximum displacement of the plate

图6 有气动力时壁板在脉冲激励下的响应(λ=450):(a)时域响应,(b)相图,(c)频谱图,(d)极限环响应云图

Fig.6 The dynamic response of the composite plate with pulse excitation (λ=450): (a) time domain response, (b) phase plot, (c) frequency response, (d) the contour of the LCO of the plate

如图7-9所示为复合材料壁板在λ=450时,不同温度下(0.75L1,0.5L2)点处的动力学响应.可见,随着温度的升高,壁板的动力学响应出现显著的变化.当温度上升一定范围时,壁板仍然作等幅极限环运动;当温度进一步升高时,壁板将做准周期颤振运动;当温度继续增高时,壁板将作混沌颤振运动,此时壁板的平衡位置不止一个且偏离原点,壁板将围绕其不同的平衡位置不断跳动.这说明了热效应会使得复合材料壁板的非线性动力学行为更加复杂.

图7 热环境下复合材料壁板的时域响应(ΔT=2.5Tcr):(a)时域响应,(b)相图

Fig.7 The dynamic response of the composite plate under thermal environment (ΔT=2.5Tcr): (a) time domain response, (b) phase plot

图8 热环境下复合材料壁板的时域响应(ΔT=5.0Tcr):(a)时域响应,(b)相图

Fig.8 The dynamic response of the composite plate under thermal environment (ΔT=5.0Tcr): (a) time domain response, (b) phase plot

图9 热环境下复合材料壁板的时域响应(ΔT=6.0Tcr):(a)时域响应,(b)相图

Fig.9 The dynamic response of the composite plate under thermal environment (ΔT=6.0Tcr): (a) time domain response, (b) phase plot

为了研究声载荷对复合材料壁板时域历程动力学响应的影响,计算不同声载荷作用下壁板的时域动力学响应结果.如图10所示,为不同声压下四边简支复合材料壁板中点处的时域动力学响应,这里分别计算了声压级为120dB和130dB下壁板的动力学响应.显而易见,复合材料壁板在声载荷激励下会在原点平衡位置附近振动,且声载荷的增大会导致壁板的响应显著增大.

图10 复合材料壁板在不同声压下的时域动力学响应

Fig.10 The dynamic response of composite plates with different acoustic loads

图11所示,为热声载荷联合作用下四边简支[0°/0°/0°/0°/0°]复合材料壁板中点处在不同温度下的时域动力学响应,这里所选取的声压级为120dB,两种温度分别为ΔT=0ΔT=0.5Tcr.显而易见,随着温度的增加,复合材料壁板的刚度由于热应力的增大而减小,从而使得壁板的响应显著增大,这意味着飞行器在巡航过程中由于气动加热引起的热载荷会导致结构的动力学响应增大.

图11 复合材料壁板在不同温度下的时域动力学响应

Fig.11 The dynamic response of composite plates with different temperatures

最后研究气动压力对声载荷激励下复合材料壁板动力学响应的影响规律.这里所选用的声载荷依旧为声压级120dB,如图12图13所示,为不同气动压力下复合材料壁板中点和1/4点(0.75L1,0.5L2)处的动力学响应对比.由结果可知,当气动压力上升一定范围时,壁板的动力学响应下降,而当气动压力进一步增大时,壁板的动力学响应显著增大.该现象主要是由于壁板的部分低阶固有频率随着气动压力的增加而增大,且气动阻尼的存在使得结构的动力学响应在气动力一定上升范围内出现下降,而当气动压力进一步增大时,结构出现了颤振失稳,壁板的动力学响应会随着气动力的增加而显著增大.故可以认为,当气动压力较小时,由于气动阻尼的存在和结构部分固有频率的上升,使得气动压力对于结构的气弹稳定性起“正面作用”;当气动压力较大时,由于结构发生颤振失稳,使得气动压力对于结构的气弹稳定性起“负面作用”.

图12 复合材料壁板在气声联合载荷下中点的时域动力学响应

Fig.12 The dynamic response at center point of composite plates with different aerodynamic pressures

图13 复合材料壁板在气声联合载荷下1/4点的时域动力学响应

Fig.13 The dynamic response at 1/4 point of composite plates with different aerodynamic pressures

2.3 颤振控制

由于NES中非线性弹簧初始处于未压缩状态,对壁板的动力学固有特性没有影响,故这里着重研究NES的质量和阻尼对于复合材料壁板临界颤振动压的影响.定义NES的无量纲参数m¯nes=mnes/mpc¯nes=cnes/(2mpω0)k¯nes=knesh2L1L2/D11,其中mp=I0L1L2ω0=D11/(ρhL12L22).出于轻量化的考虑,NES中质量不能太大,因此在后续的算例中限定NES中无量纲质量m¯nes选定在(0,0.15]范围内.如图14所示,为NES的无量纲质量和无量纲阻尼对四边简支[0°/0°/0°/0°/0°]复合材料壁板临界颤振动压的影响云图,此处NES安装在壁板背风面的中点上.显而易见,当NES的质量和阻尼同时取较大值时,壁板的临界颤振动压急剧下降,使得壁板的气弹稳定性下降.当NES质量和阻尼两个参数其中一个取较小值时,而剩余的一个设计参数在所选范围内无论其大小,壁板的临界颤振动压都将保持在一个相对较大的值.此外可以发现,当NES质量取较小值时,增大NES的阻尼值有利于提高壁板的临界颤振动压.因此,在NES参数选择时,应尽量选择较小的质量.故在接下来的壁板动力学响应控制中,选择的NES无量纲设计参数为m¯nes=0.06c¯nes=10k¯nes=1.

图14 简支复合材料壁板临界颤振动压随NES质量和阻尼的变化

Fig.14 The variation of the critical flutter aerodynamic pressure of the simply supported plate with different values of the mass and damping in NES

图15所示,为外界气动压力为400时壁板有无NES下的(0.75L1,0.5L2)点处动力学响应的对比.可见,当气动力超过临界颤振动压不是很多时,NES可以完全抑制住壁板的颤振极限环运动.

图15 有无NES复合材料壁板的响应对比 (λ=400)

Fig.15 The dynamic response of the composite plate with and without NES (λ=400)

图16所示,为外界气动动压为450时壁板有无NES下的(0.75L1,0.5L2)点处动力学响应的对比.虽然此时NES无法完全抑制壁板的颤振极限环运动,但是可以将颤振极限环振幅从5.787mm抑制至1.968mm,约降低了壁板的颤振极限环振幅的65.99%.这意味着可以大大降低壁板的动应力,从而提高壁板的疲劳寿命.

图16 有无NES复合材料壁板的响应对比 (λ=450)

Fig.16 The dynamic response of the composite plate with and without NES (λ=450)

3 结论

本文建立了一般边界下复合材料壁板的气热声弹非线性动力学模型,通过FSDT和超声速活塞理论建立了系统的能量泛函并利用Hamilton变分获得了系统的控制方程.运用改进傅里叶法获得了壁板动力学响应的解,并通过数值算例验证了模型的准确性.随后通过变参数计算,揭示了不同参数对于壁板颤振特性的影响规律.运用NES实现了壁板的颤振控制,数值计算结果表明,当外界气动力略大于壁板的临界颤振动压时,NES可以完全抑制壁板的颤振;当气动力进一步增大时,尽管NES无法完全抑制壁板的颤振,但是可以大幅度降低壁板的极限环振幅.因此,NES可在一定范围内显著抑制壁板的动力学响应.

参 考 文 献

1

吴振强程昊张伟.热环境对飞行器壁板结构动特性的影响.航空学报2013342):334~342 [百度学术

Wu Z QChen HZhang Wet al. Effect of thermal environment on dynamic properties of aerospace vehicle panel structures.Acta Aeronautica et Astronautica Sinica20133402):334~342(in Chinese) [百度学术

2

Dowell E H. Nonlinear oscillations of a fluttering plate.AIAA Journal196647):1267~1275 [百度学术

3

Navazi H MHaddadpour H. Nonlinear aero-thermoelastic analysis of homogeneous and functionally graded plates in supersonic airflow using coupled models. Composite Structures20119310):2554~2565 [百度学术

4

Xue D YMei C. Finite element nonlinear panel flutter with arbitrary temperatures in supersonic flow. AIAA Journal1993311): 154~162 [百度学术

5

Ibrahim H HYoo H HLee K S. Supersonic flutter of functionally grated panels subject to acoustic and thermal loads. Journal of Aircraft2009462): 593~600 [百度学术

6

杨智春杨飞张玲凌. 动力吸振器用于夹层壁板颤振抑制的研究. 振动与冲击2009282): 25~27 [百度学术

Yang Z CYang FZhang L L. Flutter suppression for sandwich panel using dynamic absorber. Journal of Vibration and Shock2009282): 25~27(in Chinese) [百度学术

7

赵海. 超声速流中壁板的颤振及其抑制[博士学位论文]. 哈尔滨工业大学2013(Zhao H. The panel flutter and its suppression in the supersonic flow [Ph.D Thesis]. Harbin: Harbin Institute of Technology2013(in Chinese)) [百度学术

8

Ibrahim H HTawfik MNegm H Met al. Aerothermoacoustic response of shape memory alloy hybrid composite panels. Journal of Aircraft2009465): 1544~1555 [百度学术

9

李丽丽赵永辉. 超音速下热壁板的颤振分析. 动力学与控制学报2012101): 67~70 [百度学术

Li L LZhao Y H. The flutter analysis of thermal panel under supersonic flow. Journal of Dynamics and Control2012101):67~70(in Chinese) [百度学术

10

李晋李鹏张德春. 低速轴向气流中曲壁板的稳定性及分岔分析.动力学与控制学报2018166):506~513 [百度学术

Li JLi PZhang D Cet al. Stability and bifurcation analysis of a curved plate in an axial low-speed flow. Journal of Dynamics and Control2018166):506~513(in Chinese) [百度学术

11

Beloiu D MIbrahim R APettit C L. Influence of boundary conditions relaxation on panel flutter with compressive in-plane loads. Journal of Fluids and Structures2005218): 743~767 [百度学术

12

Jin GYe TSu Z. Structural Vibration: A uniform accurate solution for laminated beams, plates and shells with general boundary conditions. BerlinSpringer2015 [百度学术

13

Liew K MXiang YKitipornchai Set al. Vibration of Mindlin plates: programming the p-version Ritz method. AmsterdamElsevier1998 [百度学术

14

Zhou KHuang XTian Jet al. Vibration and flutter analysis of supersonic porous functionally graded material plates with temperature gradient and resting on elastic foundation. Composite Structures201820463~79 [百度学术

15

Gourdon ELamarque C H. Energy pumping with various nonlinear structures: numerical evidences. Nonlinear Dynamics2005403): 281~307 [百度学术

16

Parseh MDardel MGhasemi M H. Investigating the robustness of nonlinear energy sink in steady state dynamics of linear beams with different boundary conditions. Communications in Nonlinear Science and Numerical Simulation2015291-3): 50~71 [百度学术

17

钟锐陈建恩葛为民. 串联非线性能量阱的高分支响应研究.动力学与控制学报2019173):251~257 [百度学术

Zhong RChen J EGe W Met al. Research on higher branch response of series nonlinear energy sink. Journal of Dynamics and Control2019173):251~257(in Chinese) [百度学术

微信公众号二维码

手机版网站二维码