薄板件銑削顫振穩定性的非線性判據實驗研究

吳　石，邊立健，劉獻禮，宋盛罡，姜彥翠

(哈爾濱理工大學 機械動力工程學院，哈爾濱　150080)

?

薄板件銑削顫振穩定性的非線性判據實驗研究

吳石，邊立健，劉獻禮，宋盛罡，姜彥翠

(哈爾濱理工大學 機械動力工程學院，哈爾濱150080)

銑削過程中非線性動力學行為一直伴隨整個切削過程，為準確地判定和預測加工過程的顫振穩定性，基于實驗方法，研究了兩端固定薄板件銑削顫振穩定性的非線性判據。實驗中以薄板件振動信號為研究對象，基于相平面法、龐加萊法和頻譜分析了不同加工參數時的振動信號，繪制并討論了最大Lyapunov指數與主軸轉速和銑削深度的變化關系。最后以最大Lyapunov指數作為判據，通過等高線法確定銑削顫振穩定域，并和基于全離散法得出的銑削顫振穩定域進行比較分析，實驗得出了航空鋁合金7075-T6薄板件顫振穩定域的非線性判據。

薄板件；銑削顫振；非線性判據；李雅譜諾夫指數

現代航空航天對零件的重量、強度的要求逐漸提高，使薄壁結構件被廣泛應用于飛機的大梁、壁板等眾多主承力結構件上，有的成品零件壁厚僅僅1 mm。薄壁件雖然具有比強度高、相對重量較輕等眾多優點，但是薄壁件的加工卻存在諸多問題，其中切削加工過程中的振動失穩(如顫振)就嚴重制約著薄壁工件的加工質量，為實現薄板件銑削加工參數的優化，確保銑削過程的穩定性，必須采用方便有效的方法對加工過程的顫振穩定性進行準確的判定、預測。

在實際加工中主軸-刀具-工件，以及夾具系統是復雜的非線性系統[1]，傳統線性理論對于切削厚度較小時的非線性系統采用線性理論近似方法，可以滿足低速要求，但是隨著高速銑削的發展，采用傳統線性理論已經不能保證所建模型的精度，難以準確預測顫振臨界值和工件的表面位置誤差、粗糙度等。于是GRADISEK等[2-3]分析了不同切削深度對加工中非線性振動的影響，并基于分岔模型進行論證，提出顫振發生時振動信號中具有低維混沌振動現象。DAVID等[4]基于非線性動態切削力建立了二階時滯的再生型顫振模型，分析了系統振動由Hopf分岔進入混沌振動狀態，得到了Hopf分岔穩定性邊界條件以及表征值。SZALAI等[5]將刀具對工件切削簡化為碰撞行為，建立了高速銑削的非線性動力學模型，認為高速銑削時存在混沌振動，并對穩定性邊界進行了亞臨界分岔分析。STEFANSKI等[6-7]針對非線性銑削振動指出，動態系統的李雅普諾夫指數(Lyapunov Exponent，LE)是分析混沌運動的一種有效方法。這是由于動態切削力的非光滑特性等導致動力學模型的相空間是多維的，雖然計算兩個軌跡間距離及相鄰軌跡間發散度是比較困難的，但是可以采用最大Lyapunov指數來進行分析銑削系統的振動信號。

李忠群等[8]通過數值方法求解非線性銑削動力學方程，并根據穩定性判據得到了時域顫振穩定域。孔繁森等[9]基于Lyapunov指數和Kolmogorov熵分析了切削系統從無顫振狀態到顫振狀態這一轉變過程中的非線性行為的變化特征。李西彬等[10]研究中發現加工過程振動信號的非線性特征明顯，在相同工藝參數條件下卻不一定獲得完全相同的加工效果。同時，在銑削過程實驗中，發現薄壁工件振動信號的非線性特征比基于理論模型的振動信號要更為明顯和復雜[11]，目前薄壁工件的多場耦合振動的分岔與混沌運動的研究，也取得了一些成果，YEH等[12-14]對熱彈耦合矩形板振動的分岔與混沌進行了研究。

雖然非線性動力學模型在顫振預測方面占有巨大優勢，但是如何基于非線性理論判定薄壁件切削的穩定性，并在實驗的基礎上進行準確的顫振穩定域預測，建立一個針對具體工件的銑削顫振穩定性判據仍是個挑戰。本文將以航空鋁合金7075-T6兩端固定的薄板件為研究對象，探討以工件振動信號的最大Lyapunov指數作為閾值，通過實驗判定銑削加工過程顫振穩定域，確定兩端固定的薄板件加工的顫振穩定域的非線性判據，并為從理論上研究熱、力耦合作用下的薄板件切削振動的非線性穩定性奠定基礎。

1　顫振過程的分析方法

由于銑削系統的非線性(動態銑削力的非光滑性、系統的非線性阻尼和剛度等)，使切削振動具有明顯的非線性特性，而最大李雅譜諾夫(Maximum Lyapunov Exponent)指數、龐加萊(Poincare)映射及相平面圖法等分析方法，是有效地分析銑削顫振非線性動力學行為的有效手段。

1.1最大Lyapunov指數

最大Lyapunov指數是描述初始時刻兩個無限靠近的點隨時間演化而分離的特征，作為混沌運動的特征參量，表示相軌跡的最大發散程度，或對初值的最大敏感程度[15]。

設變切深銑削振動信號x1，x2，… xN(單變量時間序列)，式中N是時間序列的總數。根據PACKED等[16]提出的時間延遲思想，可以重構出所觀察到的動力學系統的相空間，基于這個思想，對時間序列進行相空間重構，得到重構的軌跡X，表示為

X=[X1,X2,…,XM]T

(1)

式中，M為相空間重構后軌跡點的個數，Xi為銑削振動系統在間斷時間點i的狀態，它可以表示為

Xi=[xi,xi+τ,…,xi+(m-1)·τ]

(2)

式中，τ為時間延遲，m為嵌入維數，其中M=N-(m-1)τ。對重構相空間X分為n段：[X1,X2,…,XT]，[XT+1,XT+2,…X2T]，…，[X(n-1)T+1,X(n-1)T+2,…,XnT]。其中每段的長度T=M/n，稱為演化時間。

(3)

式中:Δt是采樣時間間隔。

在不同的相空間維數下，對銑削時薄壁件的振動信號進行分析計算，進而研究不同加工參數下的銑削顫振穩定性。

1.2Poincare映射和相平面圖

Poincare映射的分析方法是橫坐標為系統響應的位移值，縱坐標為系統響應的速度值，表現為一個運動點，其計算方式為間隔每一個時間段進行數據提取。若銑削系統為周期振動，則對應于Poincare圖上的一個孤立點，若振動信號的周期為N時，在圖中出現N個獨立點，且周期數同孤立點數相同。如果振動信號變化劇烈發生顫振時，則在圖中顯示為離散點堆積圖[16]。

相平面圖是振動系統的解軌跡投影在相空間的運動曲線，對于確定性的運動系統，有下列幾種形式：振動信號周期為1或N時，對應相圖中出現1條或N條自封閉的曲線，振動信號具有強顫振特性時封閉的曲線會逐漸拓展為雜亂的曲線群。

2　薄板件銑削加工實驗

實驗在VDL-1000E型機床上進行，工件模型為兩端固定的薄板件，該件長120 mm，寬80 mm，壁厚為5 mm，實驗裝置與傳感器布置如圖1所示。測力儀型號為KISTLER9257B，振動加速度傳感器為PCB356A25，置于加工面的背側，其靈敏度為10.42 mv/g。刀具為山特維克(R216.64-08030-AO09G 1610)直徑8 mm、螺旋角30°的四刃硬質合金球頭銑刀。數采分析系統為東華DH5922的信號采集系統，采樣頻率為5 kHz。加工參數：順銑，徑向銑削寬度ae=0.2 mm，進給速度fz=0.1 mm/tooth，銑削實驗中工件與刀具的相對傾角設為30°，走刀路徑為沿水平方向來回進給。

表1　薄板件模態參數

假設刀具為剛性、工件為柔性，工件前二階模態參數如表1所示。首先基于刀具-工件模型雙自由度系統，根據全離散法[17]獲得的兩端固定薄板件的銑削顫振穩定域如圖3所示，根據圖中所標示的銑削加工參數點進行銑削加工實驗，即每一個參數點對應一定的切削深度和主軸轉速，這樣總共進行121個加工參數的切削振動測試。如圖2所示的兩端固定的薄壁樣件一共加工了二個，針對第一個樣件銑削實驗時將每一轉速作為一組實驗，一共進行13組實驗，轉速從3 000 r/min增加到6 000 r/min。針對第二個樣件銑削實驗時將每一轉速作為一組實驗，也一共進行13組實驗，轉速從6 000 r/min減少到3 000 r/min。

加工區域在工件上如圖2所示。每組包括10個或9個不同的切深。同組實驗過程中忽略上一加工參數實驗后由于工件材料去除對下一參數實驗所產生的影響。同時進行下一組實驗前，先對工件進行加工修模，使每組實驗工件薄壁部分的厚度及高度等尺寸相同，這樣可以盡量確保各組實驗工件的模態特性相對一致性。這兩個試驗樣件在相同參數切削時，薄壁件振動情況基本相同，總體誤差小于9%。然后針對每一種工況，基于相平面法、龐加萊法和頻譜分析了其實測的工件振動信號，判定在此加工參數下是否發生了顫振。并計算該振動信號的最大Lyapunov指數。

圖1　薄板件銑削加工實驗Fig.1　The milling experiments of thin wall

圖2　兩端固定薄板件的13組加工區域分布Fig.2　The13 groups at both ends of the fixed thin wall processing regional distribution

同時，在實驗中發現，在薄壁件加工中瞬時銑削力信號中的混沌(不確定性)特性，沒有對振動信號進行混沌(不確定性)分析明顯，故薄板件銑削加工實驗以薄板件的振動信號為研究對象。

圖4分別是兩端固定薄板件在不同加工參數時(A(3 500 r/min，0.3 mm)、B(3 500 r/min，0.5 mm)、C(4 500 r/min，0.3 mm)及D(5 000 r/min，0.7 mm))的工件加速度振動信號以及其相平面圖、Poincare截面圖及頻譜分析結果。

圖3　兩端固定薄板件的銑削顫振穩定域(○為銑削加工參數點)Fig.3　Milling chatter stability region of fixed at both ends of the thin wall(○ is the parameter point of milling)

由圖4可知，當以加工參數A(3 500 r/min，0.3 mm)進行銑削實驗時，其工件振動的加速度時域信號幅值為5 m/s2左右，研究中取振動信號的4 096個采樣點，通過濾波、離散傅里葉變換后發現能量在刀齒切削頻率(3 500×4/60=233 Hz)及諧波頻率(467 Hz)附近，且能量分布較均勻，如圖4(d)所示。從振動信號的相平面圖4(b)可知，工件的振動響應逐漸收斂為圍繞中心點的封閉曲線，而Poincare截面圖4(c)中的點也較集中且僅有數個點，此加工參數下振動信號的最大Lyapunov指數(嵌入維數等于2)雖然為0.277 3，但經測試其表面粗糙度、以及相平面、Poincare截面圖、頻譜分析后，綜合判定該加工參數時兩端固定的薄板件銑削處于無顫振的穩定切削狀態。

當以加工參數B(3 500 r/min，0.5 mm)進行銑削實驗時，工件振動加速度時域信號如圖4(e)所示，振動幅值達到11.2 m/s2左右。從振動信號的相平面圖4(f)可知，工件的振動響應逐漸發散為多條曲線且沒有規律，而Poincare截面圖4(g)中的點也變得分布不均，此加工參數下振動信號的最大Lyapunov指數(嵌入維數等于2)為0.613 1。經濾波、傅里葉變換后得到的頻譜如圖4(h)所示，發現峰值不僅僅處于刀齒切削頻率(3 500×4/60=233 Hz)及諧波頻率(467 Hz)附近，而且能量有向刀具結構頻率(1 278 Hz)處積聚的趨勢，經測試其表面粗糙度、以及相平面、Poincare截面圖、頻譜分析后，綜合判定該加工參數時兩端固定的薄板件銑削處于顫振切削狀態。

(e) 振動信號(3500r/min,0.5mm)(f) 相平面圖(3500r/min,0.5mm)(g) Poincare截面圖(3500r/min,0.5mm)(h) 頻譜分析(3500r/min,0.5mm)

(i) 振動信號(4500r/min,0.3mm)(j) 相平面圖(4500r/min,0.3mm)(k) Poincare截面圖(4500r/min,0.3mm)(l) 頻譜分析(4500r/min,0.3mm)

(m) 振動信號(5000r/min,0.7mm)(n) 相平面圖(5000r/min,0.7mm)(o) Poincare截面圖(5000r/min,0.7mm)(p) 頻譜分析(5000r/min,0.3mm)

圖4不同加工參數下工件振動信號及其相平面圖,Poincare截面圖,頻譜分析
Fig.4 Workpiece vibration signal and its phase plane，Poincare-sectional view， spectrum analysis under different processing parameters

當以加工參數C(4 500 r/min，0.3 mm)進行實驗時，振動加速度時域信號如圖4(i)所示，工件出現了明顯振動，最大振幅甚至達到13.3 m/s2左右。其振動信號的相平面圖4(j)雜亂無章，Poincare截面圖4(k)中的點也變得發散，且此參數下振動信號的最大Lyapunov指數(嵌入維數等于2)為0.611。將振動加速度信號通過濾波、傅里葉變換后發現能量向刀具結構頻率(1 278 Hz)處積聚，并且能量增大非常多，如圖4(l)所示，符合非穩定切削顫振情況。

同理，當以加工參數D(5 000 r/min，0.7 mm)進行實驗時，振動加速度時域振動信號如圖4(m)所示，工件發生輕微顫振，加速度的振動幅值在9.7 m/s2左右。從相平面圖4(n)可知，工件的振動響應逐漸發散為多條曲線，振動信號的整個相平面圖比較混亂，而圖4(o)中Poincare截面圖中的點也分布不均，該加工參數下對振動加速度信號進行最大Lyapunov指數(嵌入維數等于2)計算，其值為0.609。經傅里葉變換后得到的頻譜如圖4(p)所示，發現峰值不僅僅處于刀齒切削頻率(5 000×4/60=333 Hz)及諧波頻率(666 Hz)附近，而且能量向刀具結構頻率(1 278 Hz)處積聚，因此在該銑削參數加工時發生顫振，處于非穩定切削狀態。

通過振動信號的相平面圖、Poincare截面圖及頻譜分析，發現當振動信號的最大Lyapunov指數超過0.59以后銑削處于非穩定切削狀態。

3　非線性穩定域預測

基于相平面法、龐加萊法和頻譜分析判定在該加工參數下是否發生了顫振后，對所有實驗工況的工件振動加速度信號進行最大Lyapunov指數(嵌入維數等于2)計算，結果如圖5所示。由圖5可知，在每一轉速下隨著銑削深度的變化，其最大Lyapunov指數值也在隨之變化，而在最大Lyapunov指數(嵌入維數等于2)為0.3～0.6附近發生突變。但是，當轉速為5 250 r/min時折線在銑削深度0.9 mm處突然降低，其相平面圖、Poincare截面圖，頻譜分析也說明此加工參數處于穩定切削狀態，該點處于hof分叉葉瓣內，顫振穩定域邊界附近，切削穩定穩定性邊界存在一定的不確定性。

由于刀具轉速和銑削深度是影響銑削穩定性的兩個較主要的因素，那么通過插值方法可獲得最大Lyapunov指數(嵌入維數等于2)與主軸轉速和銑削深度的變化關系如圖6所示，這樣可以較直觀看到主軸轉速和銑削深度的變化對振動信號的非線性特征的影響。

圖5　不同銑削深度時的振動加速度信號的最大Lyapunov指數Fig.5　The maximum Lyapunov exponent of vibration acceleration signals of different milling depth

圖6　振動加速度信號的最大Lyapunov指數與主軸轉速、銑削深度的變化關系Fig.6　Changes in the relationship between the maximum Lyapunov exponent of vibration acceleration signals with the speed of spindle and milling deep

由圖6可知，在主軸轉速4 000 r/min、5 500 r/min左右分別存在一個Lyapunov指數較小的區域，選用此處加工參數易獲得較好的加工質量(該加工參數下的兩端固定薄壁件的表面粗糙度小于其它工況下的表面粗糙度)。根據圖6分別選定最大Lyapunov指數為0.6、0.59、0.61時作為閾值，畫出等高線。以此作為穩定域邊界，來分析銑削加工的穩態臨界銑削深度與主軸轉速的關系，圖7為以不同最大Lyapunov指數值閾值時銑削穩定預測圖與基于全離散法的銑削顫振穩定域的比較。

在圖7中，其中標注“×”的表示在該切削參數實驗中會導致顫振發生，而標注“○”則表示在該切削參數實驗中切削處于穩定切削狀態。其中，曲線1為采用全離散法得到的穩定性臨界曲線；曲線4為振動信號最大Lyapunov指數值為0.61時所獲得的穩定性臨界曲線，這兩條曲線的重合度較好，整體誤差度小于9%，這也同時證明了此方法的有效性。曲線2、3為不同最大Lyapunov指數值閾值時銑削穩定預測圖。本文所獲的穩定性臨界曲線在極限切深上低于全離散法曲線，這是由于加工系統非常復雜，存在著各種非線性因素所導致。

圖7　銑削加工顫振穩定域臨界曲線及實驗結果(×非穩定切削參數；○穩定切削參數)Fig.7　Milling chatter stability region critical curve and experimental results (×The non stable cutting parameters；○ The stable cutting parameters)

本文基于實驗所獲得的穩定性曲線與一些實驗結果有一定的誤差，出現誤差的原因可能由于以下的因素：① 由于動力學參數的變化銑削穩定域邊界是在一定范圍內變化的。② 利用插值技術所獲得的最大Lyapunov指數與主軸轉速及銑削深度的關系圖本身具有一定的誤差。③ 還需大量實驗數據來進一步精確此顫振穩定性臨界曲線。④ 實際加工中銑削穩定域的邊界存在一定的不確定性。

4　結　論

本文以航空鋁合金7075-T6兩端固定薄板件為研究對象，研究了以最大李雅普諾夫指數作為閾值確定銑削顫振穩定性，并以此作為穩定判據。具體結論如下：

(1)以鋁合金7075-T6兩端固定薄板件振動信號為研究對象，基于相平面法、龐加萊法和頻譜分析了不同加工參數時的振動信號，以此確定銑削過程是否發生顫振。

(2)Lyapunov指數是衡量系統動力學不確定性特征的重要判據，不規則振動中含有混沌，混沌是由銑削系統的非線性引起的，會導致不規則振動。穩定切削時與顫振切削時的不規則振動的程度不同。

(3)提取不同加工參數情況下兩端固定薄壁件振動信號的最大Lyapunov指數，繪制并討論了最大Lyapunov指數與主軸轉速和銑削深度的變化關系。最大Lyapunov指數隨銑削深度的增大而增大，實驗中發現該指數會在某一切深時發生數值突變，由此可以將振動信號的最大Lyapunov指數作為穩定性的非線性判據。

(4)針對兩端固定薄板件加工實驗，以振動信號的最大Lyapunov指數值為0.61時作為銑削顫振的非線性判據，通過等高線法得到銑削顫振穩定性臨界曲線。雖然本文所獲得的穩定性臨界曲線與基于全離散法所獲取的臨界曲線有一定的誤差，但這種誤差是可以接受的。另外，加工時可以以工件振動信號的最大Lyapunov指數作為加工參數優化的一種手段。

[1]GAO S H,LONG X H,MENG G.Noniinear response and nonsmooth bifurcations of an unbalanced machine-tool spindle-bearing system[J].Nonlinear Dynamics,2008,54:365-377.

[2]GRADISED J,GOVEKAR E,GRABED I.A chaotic cutting process and determining optimal cutting parameter values using neural networks[J].Institution Journal Machine Tools and Manufacture,1996,36(10):1161-1172.

[3]GRADISED J,GOVERKAR E,GRABED I.Time series analysis in metal cutting:Chatter versus chatter-free cutting[J].Mechannical Systems and Signal Processing,1998,12(6):839-854.

[4]DAVID E,GILSINN D E.Estimating critical hopf bifurcation parameters for a second-order delay differential equation with application to machine tool chatter[J].Nonlinear Dynamics,2002，30(2)：103-154.

[5]SZALAI R,STEPAN G,HOGAN S J.Global dynamics of low immersion high-speed milling[J].Chaos,2004,14(4):1069-1077.

[6]STEFANSKI A,KAPITANIAK T.Using chaos synchronization to estimate the largest Lyapunov exponent of non-smooth systems[J].Discrete Dynamics in Nature Sociences,2000,4(3):207-215.

[7]STEFANSKI A,KAPITANIAK T,DABROWSKI A.The largest lyapunov exponent of dynamical systems with time delay[C]//Proceedings of IUTAM Symposium on Chaotic Dynamics and Control of Systems and Processes in Mechanics,Rome,Italy,2003.6.

[8]LI Z Q,LIU Q.Solution and analysis of chatter tability for end milling in the time-domain[J].Chinese Journal of Aeronautics,2008,21(2):169-178.

[9]孔繁森，劉鵬，王曉明.切削振動加速度時間歷程演化過程的動力學特征[J].振動與沖擊,2011,30(7):10-15.

KONG Fansen,LIU Peng,WANG Xiaoming.Dynamic characteristics of evolution process of cutting vibration acceleration time history[J].Journal of Vibration and Shock,2013,30(7):10-15.

[10]劉志兵，王西彬. 微細銑削振動信號非線性特征的試驗研究[J]. 兵工學報,2010,31(1):84-87.

LIU Zhibing,WANG Xibin.An expermental study on nonlinear character istic of vibration signal in Micro-end-milling Process[J].Acta Armamentar,2010,31(1):84-87.

[11]吳石，劉獻禮，肖飛.銑削顫振過程中的振動非線性特征試驗[J].振動測試與診斷,2012,32(6):935-940.

WU Shi,LIU Xianli,XIAO Fei.Experimental study of the nonlinear characteristics of vibration in milling chatter[J].Journal of Vibrat ion,Measurement and Diagnosis,2012,32(6):935-940.

[12]王平，陳蜀梅，王知人.大撓度簡支矩形薄板受熱力磁耦合作用分岔與混沌[J].振動與沖擊,2013,32(7):129-134.

WANG Ping,CHEN Shumei,WANG Zhiren.Bifurcation and chaos of a thin rectangular plate simply upported with large deflection in a coupled environment of heating,force and magnetic field[J].Journal of Vibration and Shock,2013,32(7):129-134.

[13]YEH Y,CHEN C K,LAI H.Chaotic and bifurcation dynamics for a simply supported rectangular plate of themo-mechanical coupling in large deflection[J].Chaos,Solitons Fractals,2002,13(7):1493-1506.

[14]YEH Y.The effect of thermo-mechanical coupling for a simply supported orthotropic rectangular plate on nonlinear dynamics[J].Thin-Walled Structures,2005,43(8)：1277-1295.

[15]楊智春,張蕊麗.基于最大李雅普諾夫指數的壁板熱顫振特性分析[J].西北工業大學學報,2009,27(6):770-776.

YANG Zhichun,ZHANG Ruili.Analysis of panel therm al flutter using max imum Lyapuno v ex ponent[J].Journal of Northwestern Polytechnical University,2009,27(6):770-776.

[16]聞邦椿，李以農，徐培民,等.工程非線性振動[M].北京：科學出版社，2007.

[17]DING Y,ZHU L M,ZHANG X J,et al.A full-discretization method for prediction of milling stability[J].International Journal of Machine Tools and Manufacture,2010,50(5):502-509.

Tests for milling chatter stability nonlinear criterion of thin parts

WU Shi,BIAN lijian,LIU Xianli,SONG Shenggang,JIANG Yancui

(School of Mechanical Engineering,Harbin University of Science and Technology,Harbin 150080,China)

The nonlinear dynamic behavior of milling process accompanies the entire cutting process.In order to accurately determine and predict the chatter stability of machining process,a nonlinear criterion for milling shatter stability of a thin part fixed at its two ends was studied with tests.Vibration signals of the thin part were taken as the study object in test,based on the phase plane method,Poincare method and spectra the vibration signals under different processing parameters were analyzed and plotted.The relationships between the maximum Lyapunov exponent and the spindle rotating speed,the former and the milling depth were discussed.Finally,taking the maximum Lyapunov exponent as the criterion,the milling chatter stability domain of the thin part was determined with the contour method,compared with the milling chatter stability domain obtained with the full discrete method,the nonlinear criterion for the milling chatter stability domain of a thin part made of aviation aluminum alloy 7075-T6 was obtained with tests.

thin part; milling chatter; nonlinear criterion; lyapunov exponent

國家自然科學基金資助項目(51275139);哈爾濱市科技創新人才研究基金項目(2011RFLXG019)

2015-04-07修改稿收到日期：2015-07-17

吳石 男，博士,教授，1971年8月生

TB122;TG506

A DOI:10.13465/j.cnki.jvs.2016.17.032