葛仁余， 張金輪， 姜忠宇， 韓有民， 索小永， 牛忠榮

(1.安徽工程大學(xué) 力學(xué)重點(diǎn)實(shí)驗(yàn)室， 安徽 蕪湖 241000； 2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院， 合肥 230009)

軸向功能梯度變截面Timoshenko梁自由振動(dòng)的研究

葛仁余1， 張金輪1， 姜忠宇1， 韓有民1， 索小永1， 牛忠榮2

(1.安徽工程大學(xué) 力學(xué)重點(diǎn)實(shí)驗(yàn)室， 安徽 蕪湖 241000； 2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院， 合肥 230009)

功能梯度材料可以提高結(jié)構(gòu)的強(qiáng)度、改善質(zhì)量分布和保證工程結(jié)構(gòu)的完整性，因此軸向功能梯度變截面梁已廣泛應(yīng)用于土木、機(jī)械和航空工程。提出了用插值矩陣法計(jì)算軸向功能梯度Timoshenko梁自由振動(dòng)固有頻率；基于Timoshenko梁理論，將軸向功能梯度Timoshenko梁自由振動(dòng)固有頻率的計(jì)算轉(zhuǎn)化為一組非線性變系數(shù)常微分方程特征值問(wèn)題；運(yùn)用插值矩陣法可一次性地計(jì)算出軸向功能梯度變截面梁各階振動(dòng)固有頻率，并可同時(shí)獲取相應(yīng)的振型函數(shù)。該方法對(duì)于材料梯度函數(shù)和截面幾何輪廓的具體形式無(wú)任何限制條件，計(jì)算結(jié)果與現(xiàn)有結(jié)果對(duì)比，發(fā)現(xiàn)吻合良好，表明了該方法的有效性。

變截面梁；橫向振動(dòng)；固有頻率；插值矩陣法；功能梯度材料

功能梯度材料(Functionally Graded Material, FGM)是一種新型復(fù)合材料，將多種性能各異的材料按照設(shè)計(jì)意愿，形成材料的物理性能連續(xù)變化的組織和結(jié)構(gòu)，同時(shí)使不同材料結(jié)合部位的界面消失，避免了材料物理性能的不連續(xù)性和應(yīng)力集中，材料的性能在空間沿某個(gè)方向連續(xù)變化，以滿足各種特殊工程結(jié)構(gòu)的需要[1]。

目前，材料性能沿厚度方向梯度變化的功能梯度梁關(guān)于彎曲變形、自由振動(dòng)和穩(wěn)定性已有了大量的研究，文獻(xiàn)[2]假設(shè)材料的彈性模量沿厚度方向按指數(shù)型函數(shù)連續(xù)變化，研究了受橫向荷載作用的功能梯度簡(jiǎn)支梁的彈性解。文獻(xiàn)[3]開(kāi)發(fā)了一種新型梁?jiǎn)卧芯苛瞬牧系牧W(xué)、熱學(xué)性能沿厚度方向梯度變化的熱彈性問(wèn)題；文獻(xiàn)[4]提出遺傳算法來(lái)優(yōu)化沿厚度方向梯度變化的功能梯度梁固有頻率；文獻(xiàn)[5]分析了沿厚度方向按冪律和指數(shù)律連續(xù)變化功能梯度簡(jiǎn)支梁的解析解。

相對(duì)各向同性均勻性材料，軸向功能梯度梁的振動(dòng)問(wèn)題研究比較復(fù)雜，因?yàn)椴牧系膹椥阅Ａ俊⒓羟心Ａ亢兔芏妊剌S線方向連續(xù)變化，這些性質(zhì)導(dǎo)致梁結(jié)構(gòu)振動(dòng)問(wèn)題數(shù)學(xué)模型的實(shí)質(zhì)就是求解變系數(shù)常微分方程問(wèn)題，通常情況下很難獲得軸向功能梯度梁自由振動(dòng)的解析解，因此國(guó)內(nèi)有關(guān)軸向功能梯度梁振動(dòng)的研究文獻(xiàn)相對(duì)較少；文獻(xiàn)[6]將非均質(zhì)錐形變截面Timoshenko梁看作成許多小的均質(zhì)梁組合進(jìn)行研究；文獻(xiàn)[7]研究了阻尼對(duì)非均質(zhì)變截面Timoshenko梁的影響；文獻(xiàn)[8-9]分別運(yùn)用微分變換單元法和有限元方法研究了軸向功能梯度變截面Timoshenko梁的自由振動(dòng)和穩(wěn)定性問(wèn)題；文獻(xiàn)[10]提出了利用Fredholm積分方程來(lái)研究軸向功能梯度變截面梁的自由振動(dòng)問(wèn)題，并確定了特定邊界條件下梁的固有頻率；文獻(xiàn)[11-12] 利用半逆解法研究了梁的材料性質(zhì)按特殊形式的多項(xiàng)式函數(shù)梯度變化的自由振動(dòng)問(wèn)題；文獻(xiàn)[13-14] 針對(duì)特殊梯度變化情況研究了不同邊界條件下功能梯度變截面Eular-Bernoulli梁的自由振動(dòng)問(wèn)題。

本文提出采用插值矩陣法(Interporlating Materix Method, IMM)[15]研究軸向功能梯度Timoshenko梁自由振動(dòng)的一個(gè)新途徑?；赥imoshenko梁理論推導(dǎo)出功能梯度梁自由振動(dòng)方程，將軸向功能梯度Timoshenko梁自由振動(dòng)固有頻率的計(jì)算轉(zhuǎn)化為一組變系數(shù)常微分特征方程的特征值問(wèn)題求解，運(yùn)用插值矩陣法求解該變系數(shù)常微分方程組，可獲得軸向功能梯度梁自由振動(dòng)前若干階固有頻率及其相應(yīng)的振型函數(shù)，且變系數(shù)常微分方程組里出現(xiàn)的所有振型函數(shù)及其各階導(dǎo)函數(shù)的計(jì)算值具有同等精度，在利用振型函數(shù)及其導(dǎo)函數(shù)計(jì)算梁橫向振動(dòng)的彎矩和剪力時(shí)，這是一個(gè)顯著優(yōu)點(diǎn)， 而差分法和有限元法求解振型函數(shù)導(dǎo)函數(shù)的精度是逐次降階的。

1 基本理論和計(jì)算方法

考慮一長(zhǎng)度為l且材料性能和截面面積沿軸向x任意連續(xù)變化的Timoshenko梁，發(fā)生自由振動(dòng)時(shí)，其撓度設(shè)為w(x,t)、轉(zhuǎn)角為θ(x,t)。假設(shè)材料的彈性模量為E(x)，材料的密度為ρ(x)，截面面積為A(x)，截面轉(zhuǎn)動(dòng)慣量為I(x)，均為x的函數(shù)，即E(x)=E0f1(x)，ρ(x)=ρ0f2(x)，A(x)=A0h1(x)，I(x)=I0h2(x)，G(x)為材料的剪切彈性模量，κ為剪切系數(shù)，v為泊松比，其中，E0,ρ0,A0,I0對(duì)應(yīng)于軸向功能梯度梁在左端邊界x=0位置材料的彈性模量、密度、截面積和截面慣性矩，如圖1所示。

Timoshenko梁的自由振動(dòng)方程為

(1)

圖1 軸向功能梯度變截面Timoshenko梁Fig.1 Schematic of an axially functionally graded Timoshenko beams with varying section

(2)

本文主要考慮梁的自由振動(dòng)以及諧波振動(dòng)問(wèn)題，則

w(x,t)=W(x)sinωt，θ(x,t)=Θ(x)sinωt

(3)

可得

(4)

(5)

式中：r為影響梁橫截面轉(zhuǎn)動(dòng)慣量的無(wú)量綱回轉(zhuǎn)半徑；s為影響梁剪切變形的無(wú)量綱參量；Ω為梁的無(wú)量綱固有頻率。為描述方便，令

代入式(4)、式(5)可得

Θ″(ξ)+g111(ξ)Θ′(ξ)-g110(ξ)Θ(ξ)+
g121(ξ)U′(ξ)+λq110(ξ)Θ(ξ)=0

(6)

U″(ξ)+g221(ξ)U′(ξ)-Θ′(ξ)-g210(ξ)Θ(ξ)+
λq220(ξ)U(ξ)=0

(7)

為方便描述梁的邊界條件，采用標(biāo)記C，H和F分別表示固定、鉸支和自由3種邊界條件，如C-F表示梁的邊界條件為左端固定、右端自由，本文Timoshenko梁自由振動(dòng)的邊界條件為

簡(jiǎn)支-簡(jiǎn)支梁(H-H)

Θ′(ξ0)=0,U(ξ0)=0

Θ′(ξn)=0,U(ξn)=0

固支-固支梁(C-C)

Θ(ξ0)=0,U(ξ0)=0

Θ(ξn)=0,U(ξn)=0

固支-自由梁(C-F)

Θ(ξ0)=0,U(ξ0)=0

Θ′(ξn)=0,U′(ξn)-Θ(ξn)=0

圖2為功能梯度梁IMM計(jì)算模型，將區(qū)間[0,1]劃分為n段，0=ξ0,ξ1,ξ2,…,ξn-1,ξn=1，Δli=ξi-ξi-1=1/n利用差分法將變系數(shù)常微分方程組中兩個(gè)函數(shù)Θ(ξ),U(ξ)的導(dǎo)數(shù)值用區(qū)間劃分點(diǎn)上的函數(shù)值表示。

圖2 軸向功能梯度梁插值矩陣法計(jì)算模型Fig.2 IMM computation mode of an axially functionally graded beams

(8)

式(3)中的Θ″(ξ)和U″(ξ)用插值函數(shù)來(lái)逼近，設(shè)

(9)

式中，Li(ξ)為拉格朗日插值基函數(shù)，所以

(10)

引入向量和矩陣符號(hào)

式中，矩陣D稱為積分矩陣，僅依賴于插值基函數(shù)Li(ξ)，采用二次拋物線插值，其基函數(shù)為

將式(6)和式(7)中的Θ(ξ)函數(shù)、U(ξ)及其各階導(dǎo)函數(shù)寫(xiě)成向量形式為

Θ′=τΘ(ξ0)+σΘ′(ξ0)+DΘ″=

U′=τU(ξ0)+σU′(ξ0)+DU″=

H1=[τ,σ,D](n+1)×(n+3)

(11a)

將低階導(dǎo)函數(shù)順次地采用高階導(dǎo)函數(shù)替換，逐步遞推可得

Θ=σΘ(ξ0)+DσΘ′(ξ0)+D2Θ″=

U=σU(ξ0)+Dσ·U′(ξ0)+D2U″=

H0=[σ,Dσ,D2](n+1)×(n+3)

(11b)

Θ″=τΘ(ξ0)+τΘ′(ξ0)+IΘ″=

U″=τU(ξ0)+τU′(ξ0)+IU″=

H2=[τ,τ,I](n+1)×(n+3)

(11c)

其中，

Θ={Θ(ξ0),Θ(ξ1),Θ(ξ2),…,Θ(ξn)}T，

Θ′={Θ′(ξ0),Θ′(ξ1),Θ′(ξ2),…,Θ′(ξn)}T，

Θ″={Θ″(ξ0),Θ″(ξ1),Θ″(ξ2),…,Θ″(ξn)}T，

U={U(ξ0),U(ξ1),U(ξ2),…,U(ξn)}T，

U′={U′(ξ0),U′(ξ1),U′(ξ2),…,U′(ξn)}T，

U″={U″(ξ0),U″(ξ1),U″(ξ2),…,U″(ξn)}T

將式(6)和式(7)振動(dòng)微分方程組中變系數(shù)寫(xiě)成對(duì)角陣形式

G111=diag(g111(ξ0),g111(ξ1),…,g111(ξn)),

G110=diag(g110(ξ0),g110(ξ1),…,g110(ξn)),

G121=diag(g121(ξ0),g121(ξ1),…,g121(ξn)),

Q110=diag(q110(ξ0),q110(ξ1),…,q110(ξn)),

G221=diag(g221(ξ0),g221(ξ1),…,g221(ξn)),

G210=diag(g210(ξ0),g210(ξ1),…,g210(ξn)),

Q220=diag(q220(ξ0),q220(ξ1),…,q220(ξn)),

(12)

將式(11)和式(12)代入式(6)和式(7)中，則功能梯度Timoshenko梁自由振動(dòng)方程寫(xiě)成向量形式為

(13a)

(13b)

不失一般性，以兩端固支梁(C-C)情況為例進(jìn)行討論，則相應(yīng)的邊界條件用向量可表示為

(14)

聯(lián)立式(13)和式(14)，將Timoshenko梁自由振動(dòng)控制微分方程和邊界條件合并寫(xiě)成矩陣形式為

(15)

2 數(shù)值算例與討論

考慮一長(zhǎng)度為l的Timoshenko梁，設(shè)無(wú)量綱回轉(zhuǎn)半徑、剪切系數(shù)和材料泊松比分別為r=0.01、κ=5/6、v=0.3，設(shè)梁的橫截面面積為以下兩種情形

(16)

式中：c為變截面錐度系數(shù)；x為從梁的左端起點(diǎn)沿軸線方向的坐標(biāo)；當(dāng)c=0時(shí)，表示為等截面梁，當(dāng)c=1時(shí)，表示為錐形梁，當(dāng)clt;0時(shí)，表示為梁的橫截面從左端到右端逐漸增大。設(shè)軸向功能梯度Timoshenko梁的材料由鋁和氧化鋯組成，它們的彈性模量和密度分別為Ea=70 GPa；ρa(bǔ)=2 702 kg/m3；Ez=200 GPa；ρz=5 700 kg/m3，材料的物理性能隨坐標(biāo)x變化關(guān)系為

(17)

式中，m為影響材料性能的非均勻性參數(shù)。

2.1 變截面錐度系數(shù)c對(duì)自由振動(dòng)固有頻率的影響

運(yùn)用IMM計(jì)算軸向功能梯度Timoshenko梁自由振動(dòng)固有頻率，在不同的邊界條件下，考慮非均勻性參數(shù)m=2，錐度系數(shù)c取不同的值時(shí)，Timoshenko梁自由振動(dòng)前4階固有頻率計(jì)算值如表1～表3所列，表1～表3中給出了網(wǎng)格密度n=10,20,40的IMM計(jì)算結(jié)果，從表中可知，隨著區(qū)間劃分點(diǎn)數(shù)n的加倍，IMM計(jì)算結(jié)果幾乎按數(shù)量級(jí)加速收斂，當(dāng)n=40時(shí)，計(jì)算結(jié)果與文獻(xiàn)[8]結(jié)果完全吻合，說(shuō)明了IMM計(jì)算軸向功能梯度Timoshenko梁自由振動(dòng)固有頻率給出了非常好的結(jié)果。由表1～表3計(jì)算結(jié)果可知，邊界條件為C-C、H-H時(shí)，隨著錐度系數(shù)c值的增大，梁的第1階固有頻率減小；而邊界條件為C-F時(shí)，隨著錐度系數(shù)值的c增大，梁的第1階固有頻率增大。

2.2材料非均勻性參數(shù)m對(duì)自由振動(dòng)固有頻率的影響

運(yùn)用本文的數(shù)值計(jì)算方法，在不同的邊界條件下，考慮錐度系數(shù)c=0.2，非均勻性參數(shù)m取不同的值時(shí)，給出了網(wǎng)格密度n=10,20,40,80時(shí)變截面Timoshenko梁自由振動(dòng)前4階固有頻率計(jì)算值如表4～表6所列，當(dāng)區(qū)間劃分點(diǎn)數(shù)n=40時(shí)計(jì)算結(jié)果開(kāi)始收斂，與區(qū)間劃分點(diǎn)數(shù)n=80時(shí)計(jì)算結(jié)果趨同，其中第1階固有頻率有5位有效數(shù)字相同，第2階、第3階固有頻率有4位有效數(shù)字相同。因此，本文插值矩陣法取區(qū)間劃分點(diǎn)數(shù)n=80時(shí)的計(jì)算結(jié)果可以作為變截面Timoshenko梁自由振動(dòng)固有頻率實(shí)際結(jié)果。由表4～表6計(jì)算結(jié)果可見(jiàn)，邊界條件為C-C、C-F時(shí)，隨著非均勻性參數(shù)m值的增大，梁的第1階固有頻率減??；而邊界條件為H-H時(shí)，隨著非均勻性參數(shù)m值的增大，梁的第1階固有頻率增大。

表1 錐度系數(shù)取不同值時(shí)，兩端固支軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

表2 錐度系數(shù)取不同值時(shí)，兩端鉸支軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

表3 錐度系數(shù)取不同值時(shí)，一端固支一端自由軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

表4 非均勻性參數(shù)取不同值時(shí)，兩端固支軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

表5 非均勻性參數(shù)取不同值時(shí)，兩端鉸支軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

表6 非均勻性參數(shù)取不同值時(shí)，一端固支一端自由軸向功能梯度變截面Timoshenko梁固有頻率計(jì)算值

2.3材料非均勻性參數(shù)m對(duì)等截面梁振型函數(shù)曲線的影響

本文的數(shù)值計(jì)算方法不但可計(jì)算出軸向功能梯度Timoshenko梁自由振動(dòng)前若干階固有頻率，同時(shí)相應(yīng)的振型函數(shù)Θ(x)和U(x)一并解出，圖3為軸向功能梯度等截面Timoshenko梁在不同的邊界條件下材料非均勻性參數(shù)m=1和m=5時(shí)，第1階固有頻率對(duì)應(yīng)的振型函數(shù)Θ(x)和U(x)的分布曲線圖，從圖3可知，材料的非均勻性參數(shù)m對(duì)梁自由振動(dòng)的振型函數(shù)的影響不容忽略。

圖3 不同的邊界條件下材料非均勻性參數(shù)m對(duì)梁的振型函數(shù)分布曲線的影響Fig.3 Variation of the mode shape with the material non-homogeneity parameter (m) for different BCs

4 結(jié) 論

本文基于Timoshenko梁理論，將軸向功能梯度變截面Timoshenko梁自由振動(dòng)固有頻率的計(jì)算轉(zhuǎn)化為一組變系數(shù)常微分方程組特征值問(wèn)題，然后，運(yùn)用IMM計(jì)算已建立的變系數(shù)常微分方程組，從而獲得軸向功能梯度變截面Timoshenko梁自由振動(dòng)前若干階固有頻率，本文具有以下優(yōu)點(diǎn)：

(1) 軸向功能梯度變截面梁振動(dòng)問(wèn)題的復(fù)雜性存在于振動(dòng)方程是一組變系數(shù)常微分方程，因此只有一些特殊情況才能獲得解析解。本文方法可避免迭代方法計(jì)算超越方程，一次性地計(jì)算出軸向功能梯度變截面Timoshenko梁的自由振動(dòng)固有頻率。本文方法計(jì)算量小、計(jì)算精度高和適應(yīng)性強(qiáng)，具有一定的工程應(yīng)用價(jià)值。

(2) 本文方法既可計(jì)算出軸向功能梯度Timoshenko梁自由振動(dòng)前若干階固有頻率，同時(shí)相應(yīng)的振型函數(shù)也一并解出，基于計(jì)算結(jié)果可知：材料的非均勻性參數(shù)對(duì)振型函數(shù)的分布曲線影響較大。

(3) 解析法分析軸向功能梯度變截面Timoshenko梁的自由振動(dòng)時(shí)，只能局限于特定形式的材料梯度函數(shù)和截面幾何性質(zhì)函數(shù)，而本文方法對(duì)于材料的梯度函數(shù)和截面幾何輪廓的具體形式無(wú)需任何限制條件。

[ 1 ] BIRMAN V, BYRD L W. Modeling and analysis of functionally graded materials and structures[J]. Applied Mechanics Reviews,2007, 60(5): 195-216.

[ 2 ] SANKAR B V. An elasticity solution for functionally graded beams [J].Composites Science amp; Technology, 2001, 61 (5): 689-696.

[ 3 ] CHAKRABORTY A, GOPALAKRISHNAN S, REDDY J N. A new beam finite element for the analysis of functionally graded materials [J]. International Journal of Mechanical Sciences，2003,45(3): 519-539.

[ 4 ] GOUPEE A J, SENTHIL S V. Optimization of natural frequencies of bidirectional functionally graded beams [J]. Structural amp; Multidisciplinary Optimization，2006,32(6): 473-484.

[ 5 ] AYDOGDU M, TASKIN V, Free vibration analysis of functionally graded beams with simply supported edges [J]. Materials amp; Design，2007, 28 (5): 1651-1656.

[ 6 ] TONG X, TABARROK B, YEH K Y. Vibration analysis of Timoshenko beams with non-homogeneity and varying cross-section [J]. Journal of Sound amp; Vibration，1995,186(5): 821-835.

[ 7 ] SORRENTINO S, FASANA A, MARCHESIELLO S. Analysis of non-homogeneous Timoshenko beams with generalized damping distributions [J]. Journal of Sound amp; Vibration，2007, 304(3): 779-792.

[ 8 ] RAJASEKARAN S, TOCHAEI E N. Free vibration analysis of axially functionally graded tapered Timoshenko beams using differential transformation element method and differential quadrature element method of lowest-order[J]. Meccanica, 2014, 49(4): 995-1009.

[ 9 ] SHAHBA A, ATTARNEJAD R, MARVI M T, et al. Free vibration and stability analysis of axially functionally graded tapered Timoshenko beams with classical and non-classical boundary conditions [J]. Composites Part B Engineering，2011,42(4): 801-808.

[10] HUANG Y, LI X F. A new approach for free vibration analysis of axially functionally gradedbeams with non-uniform cross-section [J].Journal of Sound and Vibration, 2010, 329(11): 2291-2303.

[11] ELISHAKOFF I, CANDAN S. Apparently first closed-form solutions for vibrating in homogeneous beams[J].International Journal of Solids and Structures,2001,38(19): 3411-3441.

[12] CALIO I, ELISHAKOFF I. Closed-form trigonometric solutions for inhomogeneous beam columns on elastic foundation [J]. International Journal of structural Stability and Dynamics, 2004, 4(1): 139-146.

[13] HEIN H, FEKLISTOVA L. Free vibrations of non-uniform and axially functionally graded beams using Haar wavelets [J].Engineering Structures, 2011, 33(12): 3696-3701.

[14] MABIE H H, ROGERS C B. Transverse vibrations of tapered cantilever beams with end loads [J]. Journal of the Acoustical Society of America, 1964, 36(3): 463-469.

[15] NIU Zhongrong, GE Dali, CHENG Changzheng, et al. Determining stress singularity exponents of plane V-notches in bonded bimaterial [J].Journal of University of Science and Technology of China, 2008, 38(3): 314- 319.

FreevibrationanalysisofaxiallyfunctionallyTimoshenkobeamswithanon-uniformcross-section

GERenyu1,ZHANGJinlun1,JIANGZhongyu1,HANYoumin1,SUOXiaoyong1,NIUZhongrong2

(1.KeyLaboratoryforMechanics,AnhuiPolytechnicUniversity,Wuhu241000,China; 2.SchoolofCivilEngineering,HefeiUniversityofTechnology,Hefei230009,China)

Non-uniform beams with varying axially material properties are widely used in civil, mechanical and aeronautical engineering, due to the fact that they can improve distribution of strength and weight, and guarantee structural integrity. In this paper, an interpolating matrix method (IMM) for determining the natural frequencies of free transverse vibration of axially functionally graded Timoshenko beams was proposed. Firstly, based on the Timoshenko beam theory, the governing equations of free vibration analysis of an axially functionally graded Timoshenko beam were transformed into a set of nonlinear characteristic ordinary differential equations with variable coefficients. Then, the interpolating matrix method (IMM) was adopted to solve the established equations. All the natural frequencies of free transverse vibration companying with the corresponding vibration mode functions of the axially functionally graded beam were calculated at a time. Furthermore, the present methods do not pose any restrictions on both the type of material gradation and the variation of the cross section profile. By comparing with the existing results of numerical examples, the validity of the present method was confirmed.

variable cross-section beam; transverse vibration; natural frequency; the interpolating matrix method； functionally graded material (FGM)

國(guó)家自然科學(xué)基金資助項(xiàng)目(11272111)；安徽省高校自然科學(xué)研究重點(diǎn)項(xiàng)目(KJ2016A055)

2016-11-08 修改稿收到日期： 2017-01-20

葛仁余 男，博士，副教授，1969年生

張金輪 男，博士生，講師，1984年生

O326

A

10.13465/j.cnki.jvs.2017.22.025