Aerospace Engineering

Analysis of the flutter characteristics of sandwich functionally graded material panels

  • Wuchao QI ,
  • Deming CHEN ,
  • Sumei TIAN ,
  • Weitao ZHAO
Expand
  • Key Laboratory of Liaoning Province for Aircraft Composite Structural Analysis and Simulation,Shenyang Aerospace University,Shenyang 110136,China

Received date: 2025-03-23

  Revised date: 2025-04-25

  Accepted date: 2025-04-29

  Online published: 2026-06-15

Abstract

To address the delamination failure issue in composite laminates caused by interlaminar stress concentration during supersonic flight,a novel sandwich functionally graded material (FGM) panel structure was proposed. First,the nonlinear geometric relationship of the sandwich FGM panel was formulated based on the Kirchhoff thin panel theory and the von Kármán large deformation theory. The nonlinear aerodynamic forces acting on the sandwich FGM panel were simulated using the third-order piston theory. The Hamilton principle was employed to derive the differential equations of motion for the sandwich FGM panel during supersonic flight. Subsequently,the Galerkin method was introduced to discretize these differential equations of motion spatially in both the streamwise and spanwise directions,yielding the corresponding system of ordinary differential equations. Finally,the dynamic response of the sandwich FGM panel was obtained by solving the derived ordinary differential equations using the Runge-Kutta method. The results indicate that,for a constant thickness of the sandwich core,a smaller gradient index n of the core material leads to a higher flutter critical dynamic pressure; When the gradient index n of the sandwich core is less than 1.0,the flutter critical dynamic pressure exhibits a monotonic increase as the thickness ratio of the sandwich core rises; Furthermore,as dimensionless dynamic pressure increases,the motion of the sandwich FGM panel transitions progressively from static stability to limit cycle oscillation.

Cite this article

Wuchao QI , Deming CHEN , Sumei TIAN , Weitao ZHAO . Analysis of the flutter characteristics of sandwich functionally graded material panels[J]. Journal of Shenyang Aerospace University, 2026 , 43(2) : 1 -9 . DOI: 10.3969/j.issn.2095-1248.2026.02.001

复合材料层压板已被广泛应用于航空航天领域,但在实际工程中发现其受外载作用时往往会出现层间应力集中现象,进而影响板的整体性能,尤其是在将复合材料层压板运用到热、力共同作用的服役环境中时会表现得更为明显。为此,研究人员提出了功能梯度材料板1,该材料由多种组分混合而成,其材料分布在空间上呈现连续变化,从而避免了层间应力集中问题2
近年来,功能梯度材料板的力学性能及功能特性备受关注。Gholami等3研究了石墨烯片增强多孔芯金属夹层板的屈曲与后屈曲行为。Nguyen等4探讨了多向FGM多孔夹层板的非线性自由振动,分析了梯度指数与孔隙率系数的影响。Pham等5首次研究了斜向加劲肋增强泡沫金属芯FG夹层板的振动和动力响应。Singh等6重点建立了功能梯度结构中孔隙率的模型。Tran等7则基于有限元程序研究了多功能梯度夹层板的非线性振动行为。
FGM板应用于飞行器壁板时需承受气动载荷,且材料组分变化会改变其气动弹性性能,进而影响颤振特性。在这方面,Li等8研究了热-电-机械载荷下功能梯度压电板的非线性气动弹性特性与主动控制。Muc等9运用解析法和瑞利里兹法分析了多孔FGM板的颤振特性。Zanussi等10采用等几何方法分析了任意形状FGM板的非线性气动弹性行为。Zanussi等11则探讨了超音速流下碳纳米管增强 FGM板的气动弹性。Su等12研究了弹性约束加筋功能梯度板的振动和颤振行为。Tao等13研究了三向功能梯度矩形板和楔形板在静态平衡状态下的热后屈曲、气动弹性颤振及热致后屈曲颤振。Zhang等14提出了一种用于不规则板结构的离散耦合建模方法对功能梯度材料结构进行主动颤振控制。Zhong等15研究了在任意偏航角的超音速气流作用下,磁-电-热-弹性功能梯度板的颤振行为。Khorshidi等16提出了一种用于热环境中具有功能梯度面层的夹层板颤振分析的解析模型。Abdollahi等17研究了在超声速气流作用下FGM多孔斜板的气动热弹性颤振特性。Khalafi等18研究了超声速流场中含裂纹的FGM板的气动热弹性板颤振特性。
从上述研究来看,目前关于FGM板的颤振行为的研究还停留在传统FGM板。相比之下,本文所提出的夹层功能梯度板在上面板中使用了金属陶瓷材料,增强了壁板表面的硬度和耐磨性,而在下面板中使用纯金属材料,增强了壁板在气动载荷作用下的强度和刚性。同时,夹芯层的存在保证了材料性能在上、下面板之间的均匀过渡,从而减小或消除了层间应力。另外,在服役过程中,不同位置壁板受到的载荷工况不同。此时,可通过简单调整夹芯层厚度、材料成分及材料所占比例等参数,将其铺设至相应位置。
目前,在工程实践中,通常使用金属材料板或者复合材料层压板。与金属材料板相比,使用复合材料层压板可以在保证颤振速度的前提下,通过调整纤维铺层角、铺层顺序等参数获取最大刚度,从而实现大幅减重。但是,在实际服役过程中,复合材料层压板往往会因为层间应力梯度大而出现纤维与基体脱离及分层现象,同时还可能伴随基体开裂等现象。相比复合材料层压板,本文所提出的夹层功能梯度材料板可以在继承复合材料板优势的前提下,最大程度地减小层间应力所带来的不利影响。
基于以上分析,本文研究了所提出的夹层FGM板的颤振特性,考虑了不同层的厚度比、夹芯层材料的梯度指数及板的几何尺寸等因素对夹层FGM板颤振特性的影响规律,并分析了在不同无量纲动压下,夹层FGM板的非线性气动弹性响应。

1 夹层FGM板的几何模型

图1为夹层FGM板的几何模型。其中,沿顺气流方向的长度为 a;展向宽度为 b;厚度为 h。以夹层FGM板上的o为原点,以顺气流方向为 x轴,以板的展向为 y轴,以垂直于气流方向向上为z轴,建立 o x y z坐标系。夹层FGM板的表层材料为氧化锆(ZrO2)金属陶瓷,底层金属材料为钛合金 (Ti-6Al-4V),中间层为由金属陶瓷和金属材料混合而成的功能梯度材料。假设金属陶瓷层、功能梯度材料层和金属层的厚度分别为 h 1 h 2 h 3,且假设材料性能只沿板厚方向发生变化。夹层FGM板各组分材料的参数值如表1所示。
图1 夹层FGM板的几何模型
表1 夹层FGM板各组分材料的参数值
材料 E /(N·m-2 v ρ /(kg·m-3
ZrO2 168.0 × 10 9 0.28 3 800
Ti-6Al-4V 106.7 × 10 9 0.28 4 420
对于由功能梯度材料形成的夹芯层,在不同z向坐标处,其力学属性可表示为
P ( z ) = ( P m - P c ) V ( z ) + P c
式中:Pz)为功能梯度层在厚度坐标z处的等效材料属性;上角标c和m分别为陶瓷相和金属相对应的材料属性; V ( z )为夹层板中厚度 z处金属材料所占的体积比,可表示为
V ( z ) = 0                 , z 1 < z h 2 z - z 1 z 2 - z 1 n , z 2 < z z 1 1                  , - h 2 z z 2
式中: n为梯度指数,n≥0; z 1 z 2分别为 h 2 -h 1 h 2 -h 1-h 2

2 夹层FGM板运动方程的建立

基于广义Hooke定律,夹层FGM板在厚度 z处的应力—应变关系为
σ x σ y τ x y z = Q ( P ( z ) ) ε x ε y γ x y
式中:刚度矩阵 Q ( P ( z ) )表示为
Q ( P ( z ) ) = Q 11 Q 12 0 Q 21 Q 22 0 0 0 Q 66
Q i j为刚度系数,可表示为
Q 11 = Q 22 = E ( z ) 1 - v ( z ) 2 Q 12 = Q 21 = E ( z ) v ( z ) 1 - v ( z ) 2 Q 66 = E ( z ) 2 ( 1 + v ( z ) )
基于Kirchhoff薄板理论和von Kármán大变形理论,夹层FGM板的运动微分方程为
I 0 2 u 0 ( t * ) 2 = N x x x + N x y y I 0 2 v 0 ( t * ) 2 = N y y y + N x y x I 0 2 w 0 ( t * ) 2 - I 2 ( 4 w 0 x 2 ( t * ) 2 + 4 w 0 y 2 ( t * ) 2 ) = ( N x x x + N x y y ) w 0 x + ( N y y y + N x y x ) w 0 y + N x x 2 w 0 x 2 + N y y 2 w 0 y 2 + 2 N x y 2 w 0 x y + 2 M x x x 2 + 2 M y y y 2 + 2 2 M x y x y - Δ p
式中: u 0 v 0 w 0分别为夹层FGM板沿 x y z轴方向的位移; t *为描述夹层FGM板随时间变化的动力学响应过程; I 0 I 2分别为壁板截面的静矩和极惯性矩; N x x N y y N x y为板的面内合力; M x x M y y M x y为板的面内合力矩; Δ p为采用三阶活塞气动力理论模拟壁板飞行时的气动力。
对于四边简支边界条件,采用Galerkin法将夹层FGM板的3个方向的位移展开为
u 0 ( x , y , t * ) = i = 1 l j = 1 m a i j ( t * ) Φ i j ( x , y ) v 0 ( x , y , t * ) = i = 1 l j = 1 m b i j ( t * ) Φ i j ( x , y ) w 0 ( x , y , t * ) = i = 1 l j = 1 m c i j ( t * ) Φ i j ( x , y )
式(7)中模态函数为
Φ i j ( x , y ) = s i n ( i π x a ) s i n ( j π y b )
式(7)代入式(6),并沿流向和展向分别取4阶模态和1阶模态,可得四边简支夹层FGM板的运动微分方程组为
a ¨ i = F 0 0 a 0 b ( N x x x + N x y y ) Φ i j ( x , y ) d x d y b ¨ i = F 0 0 a 0 b ( N y y y + N x y x ) Φ i j ( x , y ) d x d y c ¨ i = F 1 0 a 0 b ( N x x x w 0 x + N x y y w 0 x + N y y y w 0 y + N x y x w 0 y + 2 M x x x 2 + 2 M y y y 2 + 2 2 M x y x y + N x x 2 w 0 x 2 + N y y 2 w 0 y 2 + 2 N x y 2 w 0 x y - Δ p ) Φ i j ( x , y ) d x d y
式中: F 0 = 4 a b I 0 F 1 = 4 a b a 2 b 2 I 0 + ( i 2 b 2 + a 2 ) π 2 I 2。令 y = a 11 , a ˙ 11 , , a 41 , a ˙ 41 , b 11 , b ˙ 11 , , b 41 , b ˙ 41 , c 11 ,…,c 41 c ˙ 41 T,则式(9)可以重写为
y ˙ = A y + g ( y )
式中: A为线性系数矩阵; g ( y )为非线性向量。
定义无量纲动压
λ = 2 q a 3 M a D
式中: M a为马赫数; D为夹层FGM板的弯曲刚度。对线性系数矩阵 A进行特征值计算,假设得到的特征值为 φ j = α j ± i β j。当 λ较小时,矩阵 A的所有特征值的实部均为负值。随着 λ的增大,矩阵 A的特征值实部的最大值由负变正。当 A =0时,与之对应的颤振临界动压 λ即为颤振临界动压,标志着系统从稳定状态向不稳定状态转变。本文取 μ = m a x ( α j )为矩阵 A的特征值实部的最大值。

3 算例分析

考察一个夹层FGM板的颤振特性和非线性动力响应特性,取大气密度 ρ a 0.3   k g / m 3。板沿顺气流方向的尺寸 a 1   m,沿展向方向的尺寸 b 1   m,即考察一个方形夹层FGM板的颤振特性。

3.1 模型验证

为验证本文所提方法的正确性,选取表2所示的Al/Al2O3的材料参数制成的传统结构的功能梯度板,并基于文献[19]的结果,对本文建立的动力学模型进行验证。
表2 Al/Al2O3 的材料参数
参数 给定值 参数 给定值
E A l /(N·m-2 7.0 e 10 a/m 1
E A l 2 O 3 /(N·m-2 3.8 e 11 b/m 0.5
ρ A l /(kg·m-3 2   707 h/m 0.002   5
ρ A l 2 O 3 /(kg·m-3 3   800 n 2
ν 0.3 ρ a /(kg·m-3 0.3
采用本文的理论对Al/Al2O3功能梯度板的颤振特性进行研究,得到Al/Al2O3功能梯度板的颤振临界动压 λ c r 1   078.5。与此同时,文献[19]中颤振临界动压 λ c r 1   051。将本文计算结果和文献[19]计算结果进行对比,误差为2.6%,这充分证明了本文中所提方法的正确性。

3.2 夹层FGM板的颤振特性分析

本节分别考虑不同层的厚度比、夹芯层梯度指数的变化及长宽比的变化对夹层FGM板的线性颤振特性的影响。首先,考虑夹芯层厚度的变化对夹层FGM板的颤振特性的影响,针对6 种不同厚度比的夹层FGM板的颤振临界动压展开研究。在梯度指数分别为 0.2和5.0这两种情况下,详细探究上面板、夹芯层及下面板的不同厚度比对夹层 FGM 板线性颤振特性所产生的影响 。
图2为在6种不同厚度占比下,特征值实部最大值随无量纲动压的变化曲线。从图2a可知,6种不同厚度比下的夹层FGM板的无量纲颤振临界动压分别为656.31、659.23、668.78、677.40、689.28、720.84。因此,当 n = 0.2时,增大夹芯层在厚度方向的占比可以提升夹层FGM板的颤振临界动压。
图2 在6种不同厚度比下,特征值实部最大值随无量纲动压的变化曲线
图2b可知,6种不同厚度比下的夹层FGM板的无量纲颤振临界动压分别为656.31、655.52、651.18、646.73、640.28、582.40。因此,当 n = 5.0时,增大夹芯层在厚度方向的占比会导致夹层FGM板的颤振临界动压下降。
图3为不同梯度指数下,沿厚度方向上金属材料的体积占比的变化曲线。从图3可以看出,当梯度指数 n = 1.0时,金属材料和金属陶瓷材料在夹芯层体积占比为1∶1。当 n < 1.0时,夹芯层中金属材料的体积占比高于金属陶瓷材料。在此情形下,随着夹芯层厚度增加,金属材料在夹层FGM板内的总占比上升,进而致使夹层 FGM 板的颤振临界动压增大。当 n > 1.0时,夹芯层内金属材料的体积占比低于金属陶瓷材料。此时,随着夹芯层变厚,金属材料于夹层 FGM 板内的总体占比降低,使得夹层 FGM 板的颤振临界动压随之减小 。
图3 不同梯度指数下,沿厚度方向上金属材料的体积占比的变化曲线
综上所述,当 n < 1.0时,适当增加夹芯层厚度,能够提升壁板的颤振临界动压;当 n > 1.0时,适当减小夹芯层厚度,可以有效提高壁板颤振临界动压。
进一步研究不同梯度指数对夹层FGM板的颤振临界动压的影响。首先需要明确材料沿厚度方向的分布情况。此次,选取了h 1h 2h 3=1∶2∶1、h 1h 2h 3=1∶3∶1、h 1h 2h 3=1∶4∶1 这3种材料分布模式。然后,分析5种不同梯度指数对夹层FGM颤振临界动压的影响。
图4为不同厚度比条件下,特征值实部最大值随无量纲动压的变化曲线。从图4a中可以得出,当梯度指数 n 0.2 0.5 1.0 2.0 5.0时,对应的颤振临界动压分别为668.78、663.10、658.14、654.38、651.18。从图4b中可以得出,当梯度指数 n 0.2 0.5 1 . 0 2.0 5.0时,对应的颤振临界动压分别为677.40、667.45、658.94、652.51、646.73。从图4c中可以得出,当梯度指数 n 0.2 0.5 1.0 2.0 5.0时,对应的颤振临界动压分别为689.28、673.33、659.89、649.74、640.28。
图4 不同层厚比条件下,特征值实部最大值随无量纲动压的变化曲线
综上所述,在夹层FGM板各层厚度已确定的前提下,随着夹芯层梯度指数的增大,该夹层 FGM 板的颤振临界动压呈现出减小的趋势。其原因是随着梯度指数的不断增大,夹芯层内金属材料的体积占比相应减少,这种材料组成上的变化进一步致使夹层 FGM 板的颤振临界动压随之降低。
在不改变壁板面积的情况下,考虑不同长宽比对夹层FGM板的颤振边界的影响。选定h 1h 2h 3=1∶2∶1,在考虑两种不同的梯度指数的情况下,分析长宽比对夹层FGM板颤振特性的影响。图5为不同长宽比下,特征值实部最大值随来流速度的变化曲线。从图5可以看出,当梯度指数 n = 0.2时,在 a / b = 0.5条件下,夹层FGM板的颤振速度达到了405 m/s,而在 a / b = 2.0条件下,夹层FGM板的颤振速度则为155 m/s;当梯度指数 n = 5.0时,在 a / b = 0.5条件下,夹层FGM板的颤振速度达到了395 m/s,而在 a / b = 2.0条件下,夹层FGM板的颤振速度则为150 m/s。由此可知,当夹层FGM板的长宽比 a / b = 0.5时,其颤振特性优于长宽比 a / b = 2.0。这是由于在来流方向上,随着夹层 FGM 板尺寸的增大,使得夹层 FGM 板在流体作用下更易发生颤振,即来流方向上尺寸的增加促使夹层 FGM 板的颤振更容易被激发,从而导致长宽比较大时,其颤振特性相对较弱。
图5 不同长宽比下,特征值实部最大值随来流速度的变化曲线

3.3 夹层FGM板的非线性响应分析

所选取夹层FGM板的结构参数具体为h 1h 2h 3=1∶1∶1, n = 0.2 a / b = 1.0。
图6 λ = 640.58时,夹层FGM板横向振动的时间历程图。从图6可以看出,当 λ = 640.58时,夹层FGM板在受初始扰动的作用后,随着时间的不断推进,其横向振动的振幅逐渐减小,最终夹层FGM板收敛于零点。这表明,当前动压小于夹层FGM的颤振临界动压,在当前动压下夹层FGM板处于静稳定状态。
图6 λ=640.58时,夹层FGM板横向振动的时间历程图
图7 λ = 671.67时,夹层FGM板的时间历程图和相平面图。通过对图7进行分析可知,当 λ = 671.67时,夹层FGM板在受初始扰动的作用后,随着时间的不断推进,夹层FGM板进入极限环运动的状态。这表明当 λ = 671.67时,壁板发生颤振,其振动状态为一种相对稳定的周期性振荡模式。
图7 λ=671.67时,夹层FGM板的时间历程图和相平面图
图8 λ = 746.30时,夹层FGM板的时间历程图和相平面图。从图8可以看出,夹层FGM板在初始扰动的作用下进入极限环运动状态。从图8a可以看出,夹层FGM板在 λ = 746.30时的振幅相比于夹层FGM板在 λ = 671.67时的振幅有所增大。
图8 λ=746.30时,夹层FGM板的时间历程图和相平面图
综上所述,夹层FGM板的运动状态由其所受动压与颤振临界动压的相对大小决定。当夹层FGM板所受动压小于颤振临界动压时,在初始扰动作用下,经过一定时间后,夹层FGM板最终处于静稳定运动状态,即振动逐渐衰减直至停止。当夹层FGM板所受动压超过颤振临界动压时,在初始扰动作用下,夹层FGM板发生颤振,且最终处于极限环运动状态,呈现出稳定的等幅振动模式。此外,在夹层FGM板发生颤振的情况下,其振幅与动压之间存在正相关关系,即夹层FGM板的振幅随着动压的增大而增大。

4 结论

通过对夹层FGM板的颤振特性进行研究,得到如下主要结论:
1)通过改变夹层FGM板夹芯层的厚度可以有效提升颤振临界动压。当梯度指数 n < 1.0时,增加夹芯层厚度 h 2的占比可以有效提高颤振临界动压。当梯度指数 n > 1.0时,减小夹芯层厚度 h 2的占比可以有效提高颤振临界动压。这表明,夹芯层厚度占比与颤振临界动压之间的关系并非简单的线性关系,而是会受到梯度指数这一因素的显著影响。
2)在不改变夹层FGM板夹芯层厚度的条件下,夹层FGM板的颤振临界动压随着梯度指数的增大而减小。因此,选取较小的梯度指数 n,可以提升夹层FGM板的颤振临界动压。
3)当夹层 FGM 板顺着来流方向的边长小于临边边长时,其在抵抗颤振方面表现出更好的性能,更不易发生颤振。
4)当动压小于颤振临界动压时,夹层FGM板在初始扰动作用后,夹层FGM板处于静稳定运动状态;当动压大于颤振临界动压时,夹层FGM板进入极限环运动状态。此外,当夹层FGM板已经处于发生颤振的状态下,继续增大动压,夹层FGM板的振幅也会变大。
[1]
Do D M Gao K Yang W,et al.Hybrid uncertainty analysis of functionally graded panels via multiple-imprecise-random-field modelling of uncertain material properties[J].Computer Methods in Applied Mechanics and Engineering2020368: 113116.

[2]
Ren S H Cheng C Z Yu B,et al.New refined higher-order shear deformation theories for functionally graded panels conforming to graded variations of material properties[J].European Journal of Mechanics-A/Solids202294: 104621.

[3]
Gholami R Ansari R Aghdasi P,et al.Buckling and postbuckling of embedded sandwich mode-rately thick panels with functionally graded graphene nanopanellet-reinforced porous core and metallic face sheets[J].Thin-Walled Structures2025210: 113063.

[4]
Nguyen V C Tran H Q Tran M T.Nonlinear free vibration analysis of multi-directional functionally graded porous sandwich panels[J].Thin-Walled Structures2024203: 112204.

[5]
Pham T T Tran H Q Nguyen V L,et al.Free and forced vibration of functionally graded porous sandwich panels reinforced by arbitrarily oblique stiffeners[J].Thin-Walled Structures2025210:113039.

[6]
Singh H Bhardwaj G Grover N.Modeling and static analysis of porous functionally graded and FG-sandwich panels[J].Structures202468: 107034.

[7]
Tran T T Nguyen V C Zenkour A M,et al.The nonlinear vibration analysis-based an enhanced finite element procedure of multi-functionally graded sandwich panels[J].Thin-Walled Structures2025211: 113042.

[8]
Li P Q Yang Z C Tian W.Nonlinear aeroelastic analysis and active flutter control of functionally graded piezoelectric material panel[J].Thin-Walled Structures2023183: 110323.

[9]
Muc A Flis J.Flutter characteristics and free vibrations of rectangular functionally graded porous pa-nels[J].Composite Structures2021261: 113301.

[10]
Zanussi P V Shahverdi H Khalafi V,et al.Nonli-near flutter analysis of arbitrary functionally graded panels using Isogeometric approach[J].Thin-Wall-ed Structures2023182: 110236.

[11]
Zanussi P V Shahverdi H Khalafi V,et al.Nonli-near flutter analysis of quadrilateral panels consis-ting of functionally graded carbon nanotubes reinforced composites using Isogeometric Analysis[J].Thin-Walled Structures2024198: 111701.

[12]
Su Z Wang L F Sun K P,et al.Vibration characteristic and flutter analysis of elastically restrained stiffened functionally graded panels in thermal environment[J].International Journal of Mechanical Sciences2019157-158: 872-884.

[13]
Tao C Dai T Chen Y.Thermal postbuckling and thermally induced postbuckled flutter of tri-directional functionally graded panels in yawed supersonic flow[J].Aerospace Science and Techno-logy2024154: 109491.

[14]
Zhang H Sun W Zhang Y,et al.Dynamic mode-ling and active aeroelastic flutter control of functionally graded irregular panels with piezoelectric layers based on the discrete-coupling method[J].Thin-Walled Structures2024205:112421.

[15]
Zhong R Qin B Wang Q S,et al.Investigation on flu-tter instability of magnetic-electric-thermo-elastic functionally graded panels in the supersonic airflow with any yawed angle[J].International Journal of Mechanical Sciences2021198: 106356.

[16]
Khorshidi K Karimi M.Flutter analysis of sandwich panels with functionally graded face sheets in thermal environment[J].Aerospace Science and Technology201995: 105461.

[17]
Abdollahi M Saidi A R Bahaadini R.An investigation of aero-thermo-elastic flutter and divergence of functionally graded porous skew panels[J].Composite Structures2022286: 115264.

[18]
Khalafi V Fazilati J.Panel flutter analysis of cracked functionally graded panels in yawed supersonic flow with thermal effects[J].Applied Mathematical Modelling2022101: 259-275.

[19]
谭小俊.复杂环境下梁板结构屈曲与颤振研究[D].哈尔滨: 哈尔滨工业大学,2017.

Outlines

/