-
摘要: 为了研究冲击载荷作用下考虑应力波效应弹性矩形薄板的动力屈曲,根据动力屈曲发生瞬间的能量转换和守恒准则,导出板的屈曲控制方程和波阵面上的补充约束条件,真实的屈曲位移应同时满足控制方程和波阵面上的附加约束条件。满足上述条件,建立了该问题的完整数值解法,对屈曲过程中冲击载荷、屈曲模态和临界屈曲长度之间的关系进行研究,定量计算了横向惯性效应对提高薄板动力屈曲临界应力的贡献。研究表明:板的厚宽比一定时,临界屈曲长度随冲击载荷的增大而减小;由于屈曲时的横向惯性效应,应力波作用下薄板一阶临界力参数是相应边界板的静力失稳临界力参数的1.5倍;随着边界约束逐渐减弱,板临界力参数逐渐减小,动力特征参数逐渐增大。Abstract: Dynamic buckling of thin rectangular plates under elastic compression waves caused by unaxial impact was investigated. Governing equations and a supplementary restraint-equation at compression wave front were obtained on the basis of energy transformation and conservation during buckling of plates; the real buckling displacement must fulfill the governing equations and the supplementary restraint-equation. Based on the requirements of above relations, a numerical method was established to solve the problem, and the relations of the amplitude of impulsive loads, buckling modes and critical buckling length were studied. Research results indicate that: critical buckling length decreases with the increase of impact load when the ratio of width to thickness of plate is fixed; the first order parameter of dynamic critical force is about 1.5 times of the same order parameter of static compression critical force; with weakening of boundary constraints, a smaller buckling load is required to cause the plate to buckle and transverse inertia increase gradually.
-
作为工程结构主要单元之一, 板在冲击载荷作用下的动力屈曲问题一直受到人们的普遍关注。已有的研究[1-4]中, 大多要求板具有某种形式的初缺陷, 且假定屈曲发生时结构中各截面处于均匀的轴向受力状态, 针对的是冲击载荷的后期效应, 此时考虑结构的整体屈曲, 忽略了应力波效应。近年来, 利用金属夹心构件的失稳吸能提高结构抗冲击性能的研究, 受到了重视[5-7], 在抗冲击机理研究中, 需考虑在应力波作用下的动力屈曲和后屈曲问题。D.G.Vaughn等[8-9]的研究表明, 高速冲击下板中应力波的传播与屈曲的发生是耦合的。考虑应力波效应, 施连会等[10]、韩大伟等[11]对阶跃载荷作用下薄板弹性动力屈曲进行了研究, 假设屈曲模态沿y方向为半个正弦波, 将问题转化为一维问题进行求解, 然而该方法只能求解垂直于应力波传播方向板的两个边界为简支约束的情况, 对于其他边界条件则无法求解。
本文中, 根据动力屈曲瞬间的能量率守恒准则, 导出应力波作用下矩形薄板动力屈曲控制方程和波前附加约束条件; 采用有限元离散建立求解冲击载荷作用下考虑应力波效应薄板弹性动力屈曲问题的完整解法, 定量计算横向惯性效应对提高薄板动力屈曲临界应力的贡献, 揭示屈曲模态与冲击载荷和临界屈曲长度之间的关系, 为薄板后屈曲问题的求解提供参考。
1. 动力屈曲平衡方程
在图 1所示矩形薄板, 建立直角坐标系。薄板长度为a, 宽度为b, 厚度为h, 材料密度为ρ, 弹性模量为E, 材料屈服极限为σs。假定在时刻t0=0, 平板左端作用一个面内冲击载荷Nx(Nx表示沿受载端y方向单位长度的线载荷, 并设N1 < σsh), 同时由载荷Nx引起的弹性压缩波开始向板的右端传播。弹性波波速
弹性波阵面到达固定端前, 板中压应力覆盖区域如图 1所示。假设压应力波在右端反射前结构发生动力屈曲, 在此阶段的任意时刻t, 弹性波阵面离开受载端的距离为l=ct, t≤a/c。
板在屈曲前保持平面状态, 与应力波传播方向垂直的截面上应力均匀分布。记临界屈曲的时间为t0, 用u(x, y, t0), w(x, y, t0)=0表示板在t0时刻的真实位移。记t1为屈曲发生的时间, t1=t0+Δt, Δt)为时间的一个微小增量, 屈曲时板的面内位移为u*(x, y, t1、v*(x, y, t1), 横向位移为w*(x, y, t1)。
由文献[12], 如果用广义坐标qst给定系统在真实运动中的位置, 若系统的真实运动与邻近运动的位形属于同一时刻, 用
表示广义坐标qst的邻近位形, 则有
其中δqs即是广义坐标qst的等时变分。若系统的真实运动与邻近运动的位形不属于同一时刻, 用
给定系统在邻近运动中无限接近, 为约束所允许的位置, 当只考虑一阶小量时, 有:
q∗s(t+Δt)=q∗s(t)+˙q∗s(t)Δt=qs(t)+δqs+˙q∗s(t)Δt (1) 以下, 将通过t1-t0时刻结构能量变化率守恒原理建立屈曲控制方程, 为方便计算, 假设板在同一时刻t0存在一个等时的邻近位形, 为u*(x, y, t0)、v*(x, y, t0)、w*(x, y, t0)。由等时变分原理可得δw=w*(t0)-w(t0), δu=u*(t0)-u(t0), δv=v*(t0)-v(t0)。以下记δw=w1, δu=u1, δv=v1。
应力波从冲击端传入后, 随着时间的增加, 板中承受轴压的部分增长, 在受压段达到临界长度时, 结构发生横向屈曲, 屈曲时板的动能、外力功和变形能分别为:
Kt1=ρh2∫b0∫ct10(˙w∗2+˙u∗2+˙v∗2)|t=t1 dx dy (2) Wt1=Nx∫b0u∗(0,y,t1)dy (3) Ut1=C2∫b0∫ct10(ε2x+2μεxεy+ε2y+1−μ2γ2xy)|t=t1 dx dy+D2∫b0∫ct10(χ2x+χ2y+2νχχy+2(1−ν)χ2xy)|t=t1 dx dy (4) 式中:C=Eh/(1-ν2), D=Eh3/12(1-ν2); 采用Karman平板理论中的几何关系, εx、εy、γxy表示板中面的应变分量, χx、χy、χxy表示板的曲率及扭率。
比较t1-t0时刻结构动能的变化情况, 可以得到:
ΔK=Kt1−Kt0=ρh2∫b0∫ct00(˙w∗2(x,y,t1)−˙w∗2(x,y,t0))dx dy+ρh2∫b0∫ct00(˙u∗2(x,y,t1)−˙u∗2(x,y,t0))dx dy+ρh2∫b0∫ct00(v∗2(x,y,t1)−˙u∗2(x,y,t0))dx dy+ρh2∫b0∫ct1ct0(˙w∗2(x,y,t1)+˙u∗2(x,y,t1)+v∗2(x,y,t1))dx dy (5) 式(5)右端第4项可以用中值定理简化, 认为位移满足一定的连续条件, 只考虑一阶小量, 有:
ρh2∫b0∫ct1ct0(˙w∗2+˙u∗2+v∗2)|t=t1 dx dy=ρh2cΔt∫b0(˙w∗2+˙u∗2+v∗2)|x=ct0t=t0 dy (6) 所以, 动能的变化率为:
˙K=limΔt→0ΔKΔt=ρh∫b0∫ct00(˙w∗¨w∗+˙u∗¨u∗+˙v∗¨v∗)|t=t0 dx dy+ρh2c∫b0(˙w∗2+˙u∗2+v∗2)|x=ct0t=t0 dy (7) 屈曲时面内惯性效应与横向惯性效应相比是小量, 忽略面内惯性效应, 式(7)可以简化为:
˙K=ρh∫b0∫ct00˙w∗¨w∗|t=t0 dx dy+ρh2c∫b0˙w∗2|x=ct0t=t0 dy (8) 忽略屈曲时面内新增微小载荷和位移, 同样可得, t1-t0时刻结构应变能的变化率为:
˙U=∫b0∫ct00(D(w∗,xx)∗,xx+μw∗,xxw∗,yy+μw∗,xxw∗,yy+w∗,yyw∗,yy+2(1−μ)w∗,xyw∗,xy)+Nx(u∗,x+w∗,xw∗,x))|t=t0 dx dy+12∫b0(D(w∗2,xx+2μw∗,xxw∗,yy+w∗2,yy+2(1−μ)w∗2,xy)+Nx(u∗,x+12w∗2,x))|∗=ct0t=t0 dx dy (9) t1-t0时刻外力功的变化情况为:
˙W=Nx∫b0u∗(0,y,t0)dy (10) 根据屈曲瞬间能量率守恒定律, 有:
˙W=˙U+˙K (11) 将式(8)~(10)代入式(11), 得
0=∫b0∫ct00(D(w∗⋅xx˙w∗⋅xx+μw∗,xxw∗⋅yy+μw∗⋅xx˙w∗⋅yy+w∗⋅yy˙w∗⋅yy+2(1−μ)w∗xy˙w∗⋅xy)+Nxw∗⋅x˙w∗⋅x+ρhw∗¨w∗)|t=t0dxdy+12c∫b0(D(w∗2⋅xx+2μw∗⋅xxw∗⋅yy+w∗2⋅yy+2(1−μ)w∗2⋅xy)+Nx(u∗,x+12w∗2⋅x)+cρhw∗2)|x=ct0t=t0dy (12) 将式(12)右端第1项积分, 得到关于位移w*(x, y, t0)的控制方程, 考虑到w(x, y, t0)=0, 可得:
D(∂4∂x4+∂4∂2∂y2+∂4∂y4)w1−Nx∂2w1∂x2+ρh∂2w1∂t2=0 (13) 式(12)右端第2项表示应力波波阵面处的边界条件, 根据波阵面处的连续条件, 有:
w1(ct0,y,t)=0w1,x(ct0,y,t)=0u∗,x(ct0,y,t)=0 (14) 将式(14)代入式(12)右端第2项, 得波阵面处的附加约束条件为:
w1,xx(ct0,y,t)=0 (15) 由以上推导可知, 要满足屈曲瞬间能量转换和守恒原理, 必须同时满足控制方程(式(13))和波前附加约束条件(式(15))。
2. 应力波作用下薄板弹性动力屈曲准则
由于临界屈曲时刻应力波传播距离是待定的, 目前广泛采用的商用有限元程序对该问题无法求解。而解析法仅对阶跃载荷作用下个别边界情况[10-11]可以实现, 当压应力波分布不均匀时, 式(13)中屈曲模态的解析形式不容易得到, 因而不易求解。本节中, 采用有限元离散, 寻找同时满足控制方程(13)和波前附加约束条件(15)的屈曲位移, 建立该问题的完整解法。
分离变量, 将板的屈曲模态函数写成:
w1(x,t)=ˉw(x,y)T(t) (16) 将式(16)代入式(12)右端第1项, 分离变量, 得:
¨T−λ2T=0 (17a) ∫b0∫ct00D(w2,xx+2μˉw−,xxw,yy+ˉw2,yy+2(1−μ)w2,xy)dx dy+∫b0∫ct00Nxw2,x dx dy+λ2∫b0∫ct00ρhˉw2 dx dy=0 (17b) 由式(17b)可以证明, 满足dλ/dl=0(l=ct0), 则同时满足了附加约束条件(15)。
将板用薄板矩形单元离散, 沿宽度方向划分为m+1个节点, 沿应力波传播长度方向划分为n+1个节点, 共有n×m个单元。单元几何刚度矩阵kσ=∭vN(x)CTCdv, 弯曲刚度矩阵k=∭vBTDBdv, 质量矩阵(一致质量矩阵)m=∭vρNTNdv。其中, N为形函数, C=N, x, D为薄板单元弹性矩阵, B为单元应变矩阵。
按照几何非线性有限元理论, 薄板动力屈曲控制方程可表述为:
\mathit{\boldsymbol{M\ddot \delta }} + \left( {\mathit{\boldsymbol{K}} + {\mathit{\boldsymbol{K}}_\sigma }} \right)\mathit{\boldsymbol{\delta }} = {\bf{0}} (18) 式中:M为结构整体的质量矩阵, K、Kσ分别为结构整体的弯曲刚度矩阵和几何刚度矩阵, δ为整体节点位移阵。
令整体节点位移阵:
\boldsymbol{\delta}=\boldsymbol{\bar \delta} \mathrm{e}^{\lambda\left(t-t_{\mathrm{cr}}\right)} (19) 代入式(18), 得:
\left(\lambda^{2} \boldsymbol{M}+\left(\boldsymbol{K}+\boldsymbol{K}_{\sigma}\right)\right) \boldsymbol{\bar \delta}=\bf{0} (20) 式(20)为压应力波作用下薄板发生动力分叉屈曲的特征方程, 若λ=0, 则转化为薄板静力失稳的特征方程。由式(19)可知, λ > 0时结构存在发散解。由以上分析可知, λ > 0, dλ/dl=0同时满足了屈曲控制方程(13)和附加约束条件(15), 因此对应长度l即为相应载荷作用下的临界屈曲长度。
可以看出, 若屈曲模态一定, 特征参数λ将直接决定屈曲模态放大的程度。λ是关于应力波传播距离l(l=ct0)和轴向分布载荷的函数, 临界屈曲时载荷在板中的分布一般是已知的, 所以λ将是关于长度l的函数, 由dλ/dl=0得到的是一定外载作用下, 结构屈曲发展最快模态对应的长度, 即由该准则得到的临界屈曲长度是最优模态对应的长度。
3. 薄板弹性动力屈曲求解
采用Matlab软件编制相关计算程序, 与文献结果进行对比, 验证本文理论和程序的正确性。计算时, 将应力波传播长度划分为10个单元, 板宽划分为5个单元, 改变l, 每一个长度可以得到若干特征值, 取其中最大的特征值, 绘制λ与l的关系图, 从图中找出满足条件的相关点。
3.1 算例
矩形薄板, 材料弹性模量E=71 GPa, 密度ρ=2 700 kg/m3。矩形截面:h=3 mm, b=250 mm。为方便表示, 从受载端开始, 沿逆时针方向, 对板的边界条件进行命名, 用C表示夹支边界、S表示简支边界、F表示自由边界, 如四边夹支板边界条件为CCCC。对应边界条件CCCC, 板中应力波覆盖部分应力为180、150和120 MPa的试件, 计算结果如图 2(a)所示。
在图 2(a)中标示了满足条件的相关点, 可以看出:对应应力180 MPa, 当l=120 mm时, λ=0, 该点与相应长度薄板静力失稳状态对应; 当120 mm < l < 184 mm时, 长度的微小增长将引起特征值λ迅速增大; 当l=184 mm时, 屈曲准则dλ/dl=0第一次得到满足, 即一阶临界屈曲长度为184 mm; 当l=256 mm时, 屈曲准则dλ/dl=0第二次得到满足, 即二阶临界屈曲长度为256 mm; 当l=338 mm时, 屈曲准则dλ/dl=0第三次得到满足, 即三阶临界屈曲长度为338 mm。同样可得:对应应力150 MPa, 一、二、三阶临界屈曲长度分别为204、288和376 mm; 对应应力120 MPa, 一、二阶临界屈曲长度分别为232、326 mm。由图 2(a)可以看出, 当板的厚宽比一定时, 临界屈曲长度随冲击载荷的增大而减小, 这与已有实验结论吻合。
限于篇幅, 本文中仅给出对应应力180 MPa的前三阶屈曲模态。由图 2(a)得到结构的临界屈曲长度和对应的特征值, 该特征值对应的特征参数即为结构屈曲模态。计算屈曲模态时, 将应力波传播长度和板宽划分为40×30的网格, 计算结果如图 3所示。
对应边界条件CSCS和CFCF, 几何物理参数同上, 板中应力波覆盖部分应力为180、150和120 MPa的试件, 计算所得λ与l关系如图 2(b)~2(c)所示。
3.2 与文献结果的对比
陈铁云等[13]对四边夹支矩形薄板静力失稳进行研究, 得到四边夹支矩形薄板的临界载荷近似为:
\sigma_{\mathrm{cr}}=\frac{4 \pi^{2} D}{3 h l^{2}}\left(3+3\left(\frac{l}{b}\right)^{4}+2\left(\frac{l}{b}\right)^{2}\right) (21) 由式(21)得到, 边界条件CCCC下相应薄板弹性静力失稳的临界载荷分别为193、162和129 MPa, 本文计算得到相应板的轴向压力为180、150和120 MPa, 最大误差为7.4%。
为方便对比, 引进量纲一量:
{\mathit{\Lambda}}_{1}^{(n)}=N_{x} l b /\left(\pi^{2} D\right), \quad {\mathit{\Lambda}}_{2}^{(n)}=\sqrt{\rho h / D}\left(\lambda l b / \pi^{2}\right) (22) 式中:
为n阶临界力参数, 表示屈曲载荷的大小;
为n阶动力特征参数, 表示结构横向惯性的影响。记静力失稳临界力参数为
将各种试件的数据代入式(22), 得到临界力参数和动力特征参数, 见表 1。
表 1 临界力参数和动力特征参数Table 1. Critical-load parameters and dynamic characteristic parameters边界条件 σ/MPa Λ1(0) Λ1(1) Λ2(1) Λ1(2) Λ2(2) Λ1(3) Λ2(3) 180 9.35 14.34 3.56 19.95 6.92 26.34 10.21 CCCC 150 8.77 13.25 3.01 18.70 6.14 24.41 9.03 120 8.36 12.05 2.26 16.93 5.02 - - 180 9.04 14.10 3.87 19.79 7.42 26.26 10.66 CSCS 150 8.31 13.12 3.48 18.57 6.65 24.22 9.58 120 7.64 11.90 2.99 16.88 5.78 - - 180 8.81 13.87 4.18 19.71 7.90 25.56 11.26 CFCF 150 7.99 12.73 3.83 17.99 7.20 23.44 10.32 120 7.17 11.12 3.34 15.95 6.38 - - 韩大伟等[11]对所述CSCS边界板进行研究, 得到了一阶临界力参数和动力特征参数的解析解, 换算到本文中为:
{\mathit{\Lambda}}_{1}^{(1)}=10+2(l / b)^{2}, \quad {\mathit{\Lambda}}_{2}^{(1)}=\sqrt{9-(l / b)^{4}}(b / l) (23) 由式(23)计算得到边界条件CSCS下一阶临界力参数分别为15.26、13.99和12.75, 本文计算得到的分别为14.10、13.12和11.90, 最大误差为7.6%;由式(23)计算得到一阶动力特征参数分别为4.08、3.62和3.14, 本文计算得到的分别为3.87、3.48和2.99, 最大误差为5.2%。这验证了本文理论和程序的正确性。
由表 1可以看出:(1)由于屈曲时的横向惯性效应, 应力波作用下薄板临界力参数远大于相应边界板的静力失稳临界力参数。一阶动力屈曲临界力参数是相应边界薄板弹性静力失稳临界力参数的1.5倍, 可见横向惯性效应大大提高了结构的屈曲临界力。(2)比较不同边界条件薄板同阶的临界力参数和动力特征参数, 可以发现对应边界条件CCCC、CSCS和CFCF, 临界力参数逐渐减小, 动力特征参数逐渐增大。这是因为, 板的边界约束逐渐减弱, 因而屈曲需要的作用载荷逐渐减小, 屈曲时惯性效应逐渐增大。
4. 结论
根据动力屈曲瞬间的能量率守恒准则, 导出了应力波作用下矩形薄板动力屈曲控制方程和波前附加约束条件; 采用有限元离散建立了求解考虑应力波效应时薄板弹性动力屈曲问题的屈曲准则和完整解法, 揭示了屈曲模态与冲击载荷和临界屈曲长度之间的关系, 定量计算了横向惯性效应对提高薄板动力屈曲临界应力的贡献。计算结果表明:板的厚宽比一定时, 临界屈曲长度随冲击载荷的增大而减小; 由于屈曲时的横向惯性效应, 应力波作用下薄板一阶临界力参数是相应边界板的静力失稳临界力参数的1.5倍; 随着边界约束逐渐减弱, 板临界力参数逐渐减小, 动力特征参数逐渐增大。
-
表 1 临界力参数和动力特征参数
Table 1. Critical-load parameters and dynamic characteristic parameters
边界条件 σ/MPa Λ1(0) Λ1(1) Λ2(1) Λ1(2) Λ2(2) Λ1(3) Λ2(3) 180 9.35 14.34 3.56 19.95 6.92 26.34 10.21 CCCC 150 8.77 13.25 3.01 18.70 6.14 24.41 9.03 120 8.36 12.05 2.26 16.93 5.02 - - 180 9.04 14.10 3.87 19.79 7.42 26.26 10.66 CSCS 150 8.31 13.12 3.48 18.57 6.65 24.22 9.58 120 7.64 11.90 2.99 16.88 5.78 - - 180 8.81 13.87 4.18 19.71 7.90 25.56 11.26 CFCF 150 7.99 12.73 3.83 17.99 7.20 23.44 10.32 120 7.17 11.12 3.34 15.95 6.38 - - -
[1] Ari-Gur J, Singer J, Weller T. Dynamic buckling of plates under longitudinal impact[J]. Israel Journal of Technology, 1981, 19(1): 57-64. http://www.researchgate.net/publication/307222378_Dynamic_buckling_of_plates_under_longitudinal_impact [2] Petry D, Fahlbusch G. Dynamic buckling of thin isotropic plates subjected to in-plane impact[J]. Thin-Walled Structures, 2000, 38(3): 267-283. doi: 10.1016/S0263-8231(00)00037-9 [3] Cui S, Hao H, Cheong H K. Numerical analysis of dynamic buckling of rectangular plates subjected to intermediate-velocity impact[J]. International Journal of Impact Engineering, 2001, 25(2): 147-167. doi: 10.1016/S0734-743X(00)00035-X [4] Cui S, Cheong H K, Hao H. Experimental study of dynamic buckling of plates under fluid-solid slamming[J]. International Journal of Impact Engineering, 1999, 22(7): 675-691. doi: 10.1016/S0734-743X(99)00022-6 [5] Xue Z, Hutchinson J W. Preliminary assessment of sandwich plates subject to blast loads[J]. International Journal of Mechanic Science, 2003, 45(4): 687-705. doi: 10.1016/S0020-7403(03)00108-5 [6] Xue Z, Hutchinson J W. A comparative study of impulse-resistant metallic sandwich plates[J]. International Journal of Impact Engineering, 2004, 30(10): 1283-1305. doi: 10.1016/j.ijimpeng.2003.08.007 [7] Fleck N A, Deshpande V S. The resistance of clamped sandwich beams to shock loading[J]. Journal of Applied Mechanics, 2004, 71(3): 386-401. doi: 10.1115/1.1629109 [8] Vaughn D G, Hutchinson J W. Bucklewaves[J]. European Journal of Mechanics A: Solids, 2006, 25(1): 1-12. doi: 10.1016/j.euromechsol.2005.09.003 [9] Vaughn D G, Canning J M, Hutchinson J W. Coupled plastic wave propagation and column buckling[J]. Journal of Applied Mechanics, 2005, 72(1): 139-146. doi: 10.1115/1.1825437 [10] 施连会, 王安稳.面内阶跃载荷下薄板弹性动力屈曲差分解[J].海军工程大学学报, 2008, 20(4): 25-29. http://d.wanfangdata.com.cn/Periodical/hjgcdxxb200804007Shi Lian-hui, Wang An-wen. Difference solution to elastic dynamic buckling of thin plates subject to in-plane step load[J]. Journal of Naval University of Engineering, 2008, 20(4): 25-29. http://d.wanfangdata.com.cn/Periodical/hjgcdxxb200804007 [11] 韩大伟, 王安稳, 毛柳伟, 等.弹性压应力波作用下矩形薄板动力屈曲[J].工程力学, 2012, 29(11): 12-15. doi: 10.6052/j.issn.1000-4750.2011.04.0204Han Da-wei, Wang An-wen, Mao Liu-wei, et al. Analytical resolution to dynamic buckling of rectangular thin plates under elastic compression wave[J]. Engineering Mechanics, 2012, 29(11): 12-15. doi: 10.6052/j.issn.1000-4750.2011.04.0204 [12] 梅凤翔, 刘端, 罗勇.高等分析力学[M].北京: 北京理工大学出版社, 1991. [13] 陈铁云, 沈惠中.结构的屈曲[M].上海: 上海科学技术文献出版社, 1993. 期刊类型引用(3)
1. 吴梦景,朱珏,黄栩浩. 开口截面拉胀超结构动力失稳分析. 宁波大学学报(理工版). 2025(02): 86-93 . 百度学术
2. 吴晓. 关于重物对简支梁冲击问题的讨论. 工程与试验. 2018(04): 1-3+49 . 百度学术
3. 韩大伟,王安稳. 流-固冲击载荷作用下直杆弹性动力屈曲研究. 海军工程大学学报. 2015(06): 1-5 . 百度学术
其他类型引用(9)
-