小通道氮氣-水兩相流三譜切片波動特性及其流型表征*

李洪偉 周云龍 王世勇 孫斌

1)(東北電力大學(xué)能源與動力工程學(xué)院,吉林 132012)

2)(撫順礦業(yè)中機熱電有限責(zé)任公司,撫順 113001)

(2013年1月22日收到;2013年4月7日收到修改稿)

1 引言

小通道兩相流廣泛地存在于石油、化工和能源等現(xiàn)代工業(yè)領(lǐng)域.與常規(guī)通道相比,微尺度通道的尺寸減小2—3數(shù)量級,使其兩相流動特性表現(xiàn)出與常規(guī)通道尺度下的流動特性有很大區(qū)別,這主要是因為在微小通道內(nèi)兩相間的表面張力和潤濕性等處于主導(dǎo)地位.由于氣液兩相間內(nèi)在的復(fù)雜性,有關(guān)微小通道內(nèi)氣液兩相流動特性的研究僅僅處于初級階段,但是因其在化學(xué)、化工和能源方面得到廣泛的應(yīng)用,所以一直以來都受到人們的廣泛關(guān)注[1,2].

兩相流是一個復(fù)雜的非線性動力學(xué)系統(tǒng),自20世紀(jì)90年代以來,基于混沌、分形、復(fù)雜網(wǎng)絡(luò)以及時頻域分析等方法進行流型識別的研究日趨增多[3],文獻[4—7]采用多個混沌參數(shù)指標(biāo)、有限穿越可視圖以及相空間多元圖重心軌跡等方法分析了油水和氣液兩相流型電導(dǎo)信號的動力學(xué)特性,在揭示兩相流型流動機理方面取得了一定進展.Sun等[8,9]應(yīng)用自適應(yīng)最優(yōu)核時頻分析方法對水平文丘里管內(nèi)氣液兩相流型與流動特性進行了識別與分析.Du等[10]同樣應(yīng)用自適應(yīng)最優(yōu)核時頻表征方法對氣液兩相流電導(dǎo)信號進行了分析,通過提取總能量和時頻熵兩個特征值準(zhǔn)確區(qū)分出不同流型.Manfredo[11]基于氣液兩相流圖像信號相密度時頻分布特征對其流動特性進行了表征.Ommen等[12]對兩相流時頻域分析方法進行了綜述,分析了在表征兩相流流型動力學(xué)特性方面的各自優(yōu)勢與不足.

高階譜估計技術(shù)是分析非高斯信號的有力工具,它從更高階概率結(jié)構(gòu)表征隨機信號,可以彌補二階統(tǒng)計量(功率譜)不包含相位信息的缺陷.高階譜有很強的消噪能力,能夠辨識非因果、非最小相位系統(tǒng),檢測和表征系統(tǒng)的非線性等[13].高階譜中Welch,AR和EV功率譜近年來被廣泛應(yīng)用,但幾種分析方法或多或少存在一些問題,最為普遍的缺點是抗噪能力不強且喪失了相位信息,這在非線性系統(tǒng)分析中是致命的弱點.而由文獻[14]演化而來的三譜估計可以有效地規(guī)避上述問題,其1.5維切片對系統(tǒng)內(nèi)在動力學(xué)特性揭示更加直觀及具體.

目前流型研究大都集中在常規(guī)管道和微米級通道,而介于兩種尺度之間的毫米級小通道(尤其是1.0—5.0 mm)的研究成果相對匱乏,相關(guān)文獻報道很少.基于此,本文應(yīng)用基于三譜切片的波動特性理論,針對內(nèi)徑2.0 mm×0.81 mm矩形垂直玻璃管道內(nèi)的氣液兩相流流型動力學(xué)特性進行了研究.

2 理論基礎(chǔ)

2.1 高階譜

高階譜分析是一種新的信號處理的有效技術(shù),它可以定量描寫信號間的非線性耦合,可以抑制噪聲、保留信號的相位信息,是處理非線性、非高斯信號的有力工具[15,16].

設(shè)系統(tǒng)的傳遞函數(shù)為H(X),而a(t)為零均值的非高斯白噪聲,由于采樣信號常受到加性高斯噪聲u(t)的干擾,用傳統(tǒng)的時間序列建模無法排除加性的高斯噪聲干擾.設(shè)y1(t)為y(t)中去掉確定性信號的零均值的有色隨機信號,利用y1(t)的高階累積量進行建模就可以抑制u(t)的干擾.對(2)式中的高階累積量 cy,k(τ1,τ2,···,τk-1) 進行 k-1 維的離散Fourier變換,得到高階譜的一般定義為

(1)式中,ω為頻率,H(ω)為系統(tǒng)的傳遞函數(shù),H?(ω) 為 H(ω) 的共軛函數(shù),Sy,k(ω1,ω2,···,ωk-1)為高階譜,也被稱為高階累積量譜.當(dāng)k=2,3,4時,它分別表示功率譜 P(ω)、雙譜B(ω1,ω2)和三譜T(ω1,ω2,ω3).

2.2 自回歸三譜

自回歸(AR)模型,

式中,ai(i=1,2,···,p)為自回歸系數(shù),p為自回歸模型的階數(shù).對于穩(wěn)定的線性物理過程h(t),考慮到系統(tǒng)為最小相位系統(tǒng),根據(jù)(1)式得到基于AR模型的雙譜和三譜表達式分別為

對y1(t)運用AR模型的參數(shù)法估計出模型系數(shù)a,結(jié)合(4)式將系數(shù)a代人(2)式和(3)式,得到AR雙譜和三譜幅值為

2.3 切片譜

二次相位耦合是常見的一種現(xiàn)象.對于這種非線性耦合現(xiàn)象,僅用功率譜很難從根本上解決問題,而高階譜就可以定量地描述二次相位耦合程度.當(dāng)凍結(jié)三譜 TAR(ω1,ω2,ω3)中的一個頻率 ω3,即當(dāng)ω3=C(const)時,三譜切片就變成了雙譜[17].定義三譜二維切片為

其規(guī)一化后的幅值為

用三譜切片譜函數(shù)計算二次相位耦合(信號兩種主要振蕩模式相互耦合),存在計算量大、三維圖形不夠直觀等缺點.為了改善圖形顯示不夠直觀的缺點,采用二維的對角切片譜函數(shù)表征系統(tǒng)的1二次相位耦合程度.當(dāng)ω1=ω2時,定義雙譜的12維(對角)切片譜函數(shù)為

2.4 波動特征參數(shù)

概率密度分布為非正態(tài)分布的隨機信號統(tǒng)稱非高斯信號,在工程中通常用偏斜因子S,平坦因子K和波動強度σ/〈p〉三個參數(shù)來描述.本文將這三個參數(shù)應(yīng)用到三譜切片的分析中,實現(xiàn)了三譜切片譜的定量分析.

式中x為離散時間序列,N為序列長度,E(·)為括號內(nèi)均值.波動強度反映了序列幅值變化梯度的強烈程度,偏斜因子是衡量隨機信號的分布偏離對稱分布的歪斜程度,偏斜因子不等于零的信號必定服從非對稱分布;而平坦因子表征統(tǒng)計頻率曲線接近分布中心時的大致狀態(tài).這三個參數(shù)都從不同方面反映出了序列的復(fù)雜程度.

3 方法評價

為了評價本文提出的方法在時間序列的處理分析中的有效性,本文主要將抗噪性與復(fù)雜度表征性作為評價指標(biāo),從具有代表性的分形、周期、混沌時間序列著手考察該方法在處理不同類型時間序列所表現(xiàn)出的抗噪與表征能力.

3.1 分形時間序列

本文對Brown分形序列進行分析,考察三譜切片波動方法在處理分形序列時抗噪與辨識方面所表現(xiàn)出的能力.

分形Brown運動(fractional Brownian motion,fbm)是典型的隨機性分形,其方程 f(x)是一個實值隨機方程{：}

其中,x是E維歐空間RE中的一點;F(t)為隨機變量t的分布函數(shù),該隨機變量服從標(biāo)準(zhǔn)正態(tài)分布 N(0,σ2);‖Δx‖ 為采樣間隔;H 為 Hurst指數(shù)[6].本文采用Matlab中的wfbm函數(shù)生成長度分別為10000,5000,1000,500的分形Brown序列各5組,分別取Hurst指數(shù)H=0.1,H=0.3,H=0.5,H=0.7,H=0.9,結(jié)果示于圖1.

本文對Brown分形序列的分析有兩個目的：一是驗證三譜切片的抗噪性能,二是討論三譜切片計算過程中數(shù)據(jù)的長度與快速傅里葉變換(FFT)長度對辨識結(jié)果的影響.

圖1 分形Brown運動時間序列

圖2 是Brown序列三譜圖,Brown分形序列雖然表面上看雜亂無章,但從圖中還是可以看到Brown分形序列的內(nèi)在規(guī)律性.二次耦合點主要集中在中心部位,接近于高斯分布.圖3是分形Brown序列在不同噪聲強度干擾下的三譜切片圖.從圖3中可以看出,五種強度噪聲干擾下的功率譜線趨勢基本相同,只存在微小的細節(jié)波動,再次驗證了三譜切片良好的抗噪性能.

在計算同一種時間序列的三譜切片及波動特征值過程中,我們發(fā)現(xiàn)有兩個參數(shù)對計算結(jié)果的影響比較大,一個是所取的序列長度N,另外一個是在進行FFT變換過程中所選擇的長度nFFT.本文對分形Brown運動時間序列在取不同H值情況下進行了N與nFFT的正交計算,觀察每相鄰兩組H值的三個波動特征參數(shù)類間距值,長度為10000時的具體數(shù)據(jù)如表1所示.

圖2 分形Brown序列三譜圖(H=0.1,N=10000)

圖3 Brown序列不同加噪下三譜圖(H=0.1,N=10000)

從表1中可以很清楚地看到,在沒有噪聲干擾下的分形Brown序列隨著nFFT的減小三個波動特征參數(shù)間的差值逐漸減小,當(dāng)nFFT為64時,每相鄰兩組序列的間距值只有一組不為0,且這一組中也只有波動強度值不為0,平坦因子和偏斜因子也為0.這表明在同一長度N下nFFT值越大,波動特征參數(shù)間的差值越大,辨識度越高.對于不同長度的分形Brown序列又有什么規(guī)律,在加入噪聲干擾情況下規(guī)律是否發(fā)生變化,圖4給出了詳細的說明.從圖4中可以看出,在沒有噪聲干擾的情況下不同長度序列在取得最大類間距時的nFFT長度都固定為1024,是一條平行于x軸的直線,而在有噪聲干擾的情況下,呈現(xiàn)出近似反比關(guān)系,當(dāng)N越大nFFT取值越小.

表1 分形Brown序列不同n FFT時三個參數(shù)間距值

圖4 Brown序列取得最大類間距值下的N與n FFT曲線

表明在時間序列嚴格符合高斯正態(tài)分布的情況下,N與nFFT兩者越大辨識效果越好,在有不規(guī)則噪聲干擾從而破壞正態(tài)分布結(jié)構(gòu)的情況下,N與nFFT成反比關(guān)系.這也為后面流型差壓序列分析時的參數(shù)選擇提供了理論依據(jù).

3.2 周期、混沌時間序列

本文為了考察三譜切片及其波動理論適用的廣泛性,對正弦、正弦+噪聲混合、Lorenz以及Lorenz+噪聲混合四種典型信號進行分析計算.幾種時間序列的產(chǎn)生條件如下：

1)正弦信號,y1=sin(x),采樣間隔為π/50;

2)正弦+噪聲,y=y1+p×y2,其中y1為正弦序列,y2為白噪聲序列,p為隨機成分的比例,這里取p=0.2;

3)Lorenz混沌信號,由Lorenz方程

初始條件x=2,y=2,z=20,采用四階龍格-庫塔方法迭代,取變量x為仿真序列[18];

4)Lorenz+噪聲,3)中所得變量x+40 dB噪聲.

圖5 Lorenz與sin信號的加噪前后三譜切片圖

圖5 為抗噪性能檢測,從圖中很容易看出,三譜切片在處理混沌序列時抗噪效果非常好,加噪后與加噪前波動趨勢幾乎如出一轍,而在處理周期序列時受噪聲擾動影響十分明顯,抗噪性能相對較差.而圖6則為四種信號三譜切片的波動特征值,從圖中可以看出,每種信號的三個波動參數(shù)值排序為：

Lorenz+噪聲＞Lorenz＞sin+噪聲＞sin,這也正表征出四種信號的內(nèi)在復(fù)雜程度.

圖6 Lorenz與sin信號的加噪前后波動特征值分布

通過對分形序列、周期序列和混沌序列的應(yīng)用分析,得出結(jié)論：三譜切片結(jié)合波動特性理論在處理分形和混沌序列時具有較好的抗噪能力,處理周期信號時抗噪能力相對較差.在復(fù)雜程度表征方面,無論應(yīng)用在哪種類型的序列中都具有很好的效果.此外,在計算三譜切片過程中,序列長度N與nFFT的選擇應(yīng)遵循如下原則：在時間序列嚴格符合高斯正態(tài)分布的情況下,N與nFFT兩者越大辨識效果越好,在有不規(guī)則噪聲干擾從而破壞正態(tài)分布結(jié)構(gòu)的情況下,N與nFFT近似成反比關(guān)系.

4 流型辨識及其動力學(xué)機理分析

實驗系統(tǒng)如圖7所示.該實驗系統(tǒng)主要包括兩部分,即氣液兩相流動與控制系統(tǒng)、數(shù)據(jù)采集系統(tǒng).實驗段采用透光性較好的有機玻璃制成.在通道相距200 mm處取兩個直徑為0.5 mm的取壓小孔,上游小孔距通道入口100 mm,下游小孔距通道出口80 mm,矩形通道截面寬為2.00 mm,縫隙為0.81 mm.

實驗中數(shù)據(jù)的采集主要有壓差信號、單相氣相流量和單相液相流量.其中：壓差信號采用羅斯蒙特3051C差壓變送器進行測量,測量精度為±0.075%,更新響應(yīng)時間為100 ms/次,差壓校驗量程為0.5inH2O-2000psi,實驗中采集壓差信號的采樣頻率為256 Hz;液體流量測量采用HW5智能金屬管浮子流量計,此流量計以MCU微處理器為核心,瞬時值和累計值的精度為±0.5%;氣體流量測量采用型號為D-600MD數(shù)字式氣體流量計,該流量計的量程范圍為0—5 SLM,精度等級為±1%.實驗中實驗環(huán)境壓力為0.1 MPa±0.0005 MPa,實驗環(huán)境溫度為25°C±0.5°C,表觀氣速為0.1—30 m/s,表觀液速為0.01—5 m/s.

在圖7中實驗段上采集到四種典型流型的差壓信號,并進行了三譜切片計算,結(jié)果如圖8所示.

圖7 氮氣-水兩相流實驗系統(tǒng)

圖8 氮氣-水兩相流典型流型差壓信號、三譜圖、切片圖 (a)毛細泡狀流;(b)分散泡狀流;(c)段塞流;(d)環(huán)狀流

從圖8中可以看出,不同流型的三譜圖表現(xiàn)差異較為明顯,這是由于不同流型氣液兩相界面分布不同引起壓降變化不同的結(jié)果.相較于三譜3D立體圖的復(fù)雜等頻對角切片更加直觀,更能體現(xiàn)不同流型的相分布特點.主要通過對二次耦合譜線的判斷挖掘氣相與液相兩種振蕩頻率與相位的耦合特性.小通道內(nèi)氮氣-水分布情況最為復(fù)雜的是環(huán)狀流,其對角切片呈現(xiàn)出頻率峰值分布混亂的特點,主要頻率峰值出現(xiàn)在-π/2,-7π/16,-π/4,π/4以及π/16處,其他很多處也存在微小的耦合現(xiàn)象,說明氣液兩相振蕩并無明顯耦合處,對壓降分布的影響也比較隨機,主要是由于環(huán)狀流是由段塞流演變而來,在段塞流中氣流量繼續(xù)增加使得氣塞被擊碎,形成不規(guī)則振蕩式運動模式,附在小通道表面的液膜在氣相推動下呈現(xiàn)出混沌、分形特征的抖動,而這種抖動的不規(guī)則性類似于信號當(dāng)中加入了白噪聲,其二次相位耦合也相對最為復(fù)雜.雖然耦合峰值不穩(wěn)定,但從圖8(d)可以看到環(huán)狀流在不規(guī)則的切片譜上還是存在幾處相對較為明顯的峰值.這種不穩(wěn)定峰值形成的原因有兩種：一方面是實驗系統(tǒng)中氮氣瓶輸出氣體時的波動以及氣體流經(jīng)三通時分配不均勻形成的波動.而這種波動反映到實驗段中存在一定的延遲,使得實驗段中環(huán)狀流的液膜被氣相沖擊產(chǎn)生振蕩存在一定的間歇性,這樣就形成了切片譜中的耦合雙峰值;另一方面是環(huán)狀流的波動近似混沌、分形并摻有白噪聲,這種不規(guī)則振蕩信號本身存在一定的相位周期性,這種周期性使得環(huán)狀流切片譜中出現(xiàn)類似圖3和圖5的幾處峰值.

從圖8(a)和圖8(b)中可以看到毛細泡狀流在-π/4到π/4之間存在四處明顯峰值,而分散泡狀流外側(cè)兩個峰值與毛細泡狀流較相近,中間兩個峰值相對較小,這與兩種不同類型的泡狀流氣泡尺寸大小與間隔頻率不同有關(guān),毛細泡狀流中氣泡尺寸分布均勻,間隔頻率具有很好的規(guī)律性,所以其兩相振蕩頻率耦合譜線明顯,而分散泡狀流氣泡尺寸分布不如毛細泡狀流均勻.同時兩種泡狀流的運動特征存在明顯的差異.氣液兩相流體以毛細泡狀流型流動時,生成的均勻氣泡與矩形管道的兩壁相貼,氣泡的截面積和小通道的橫截面積近似等同,液相不再連續(xù).隨著氣流量的波動,所生成的氣泡時間間隔大概在平均值的上下波動,這樣氣泡與液塞以幾種周期間歇性運動.由于三譜切片主要考察的是信號本身兩種主要振蕩模式的耦合現(xiàn)象,將高于平均值的時間間隔定位一種周期,那么平均值以下的是另外一種周期,所以從切片圖中可以看到明顯四峰值(兩種周期).而分散泡狀流由于不同尺寸的氣泡均單獨且分散地存在于連續(xù)液相之中,不具備毛細泡狀流那樣良好的周期間歇性,但由于氣體流量的波動性,使得微小氣泡在隨著液相的運動過程所出現(xiàn)的量不完全一致,時多時少,正是這種原因使得分散泡狀流的差壓波動信號主振蕩模式被氣泡數(shù)量大和小的兩種群體的間歇性出現(xiàn)所主宰,所以其切片譜中的雙峰值現(xiàn)象明顯,由于周期時間較長,所以峰值出現(xiàn)在外側(cè).而靠內(nèi)側(cè)的兩個小的峰值是散落在液相當(dāng)中小氣泡連續(xù)不間斷出現(xiàn)使得差壓信號所產(chǎn)生的微小波動所致,由于周期性不明顯,所以峰值不夠突出.而段塞流由于其特殊的流動結(jié)構(gòu)決定了它的氣液兩相振蕩模式耦合特別明顯(如圖8(c)所示),而且由此可以給出段塞流動的頻率,在π/3左右(一種周期).

圖9給出了四種典型流型差壓信號的三譜切片波動特征值分布情況,通過特征值的分布可以實現(xiàn)對流型的辨識.通過圖9(d)可以發(fā)現(xiàn)S與K相結(jié)合對段塞流的分類效果要明顯好于其他兩組組合,說明針對周期性較強的時間序列,由于其與高斯序列差距明顯,應(yīng)用S與K處理分析時會與其他接近分形、混沌的時間序列有較明顯的差距,而在毛細泡狀流、分散泡狀流與環(huán)狀流的識別上效果不是很理想,主要原因是由于這三種流型差壓信號波動周期性不強,具有分形與混沌特征,在這種情況下,σ/〈p〉就顯示出了它在處理接近高斯分布序列的優(yōu)勢,從圖9(b),(c)可以看到,環(huán)狀流的分類效果相對還是比較理想,但這三個波動特征參數(shù)任意兩兩組合的效果都不是十分完美,始終在毛細泡狀流與分散泡狀流的分類上存在一些瑕疵,總有部分流型樣本混雜在一起,難以準(zhǔn)確識別.而將三者結(jié)合到一起得到了對四種流型的滿意識別結(jié)果(圖9(a)).

5 結(jié)論

本文提出了一種基于三譜切片的波動特征分析理論,為了驗證該方法的適用性,將其應(yīng)用到分形、周期以及混沌時間序列,主要驗證了該方法的抗噪能力與對復(fù)雜度表征能力.結(jié)果表明該方法在處理分形與混沌時間序列時抗噪能力比較突出,而對于周期序列相對較差,基于三譜切片的波動特征值 (σ/〈p〉,S,K)能夠很好地表征出系統(tǒng)的復(fù)雜程度.同時發(fā)現(xiàn)在時間序列嚴格符合高斯正態(tài)分布的情況下,N與nFFT兩者越大辨識效果越好,而在有不規(guī)則噪聲干擾從而破壞正態(tài)分布結(jié)構(gòu)的情況下,N與nFFT近似成反比關(guān)系.

通過將該方法應(yīng)用到小通道氮氣-水兩相流差壓信號中,發(fā)現(xiàn)四種典型流型三譜切片差異明顯,主要振蕩模式的耦合現(xiàn)象可以揭示出不同流型的內(nèi)在流動機理,環(huán)狀流較泡狀流與段塞流更為復(fù)雜,其氣相與液膜相互作用顯示出強烈混沌特征,段塞流由于明顯的周期性使得對角切片分布最為簡單明了.

流型差壓信號的三譜切片結(jié)合三個波動特征值實現(xiàn)了對四種典型流型的準(zhǔn)確識別.其中S,K對非高斯信號識別效果比較突出,而σ/〈p〉對高斯信號的辨識度更高,將三者結(jié)合到一起可以使各自優(yōu)勢互補,實現(xiàn)小通道內(nèi)氮氣-水兩相流型的準(zhǔn)確識別.如將其應(yīng)用于其他多相流流型識別領(lǐng)域,以及在其他模式識別領(lǐng)域進行拓展,也將是有益的探索與補充.

[1]Ding GL,Huang D P,Zhang CL 2002 Journal of Refrigeration 2 7(in Chinese)[丁國良,黃東平,張春路2002制冷學(xué)報2 7]

[2]Li H W,Zhou Y L,Liu X,Sun B 2012 Acta Phys.Sin.61 030508(in Chinese)[李洪偉,周云龍,劉旭,孫斌2012物理學(xué)報61 030508]

[3]Gao Z K,Jin N D 2008 Acta Phys.Sin.57 6909(in Chinese)[高忠科,金寧德2008物理學(xué)報57 6909]

[4]Zong Y B,Jin N D,Wang Z Y 2009 Acta Phys.Sin..57 7544(in Chinese)[宗艷波,金寧德,王振亞2009物理學(xué)報57 7544]

[5]Dong F,Jin ND,Zong Y B 2008 Acta Phys.Sin.57 6145(in Chinese)[董芳,金寧德,宗艷波2008物理學(xué)報57 6145]

[6]Zhou T T,Jin N D,Gao Z K 2012 Acta Phys.Sin.61 030506(in Chinese)[周婷婷,金寧德,高忠科2012物理學(xué)報61 030506]

[7]Zhao JY,Jin N D 2012 Acta Phys.Sin.61 094701(in Chinese)[趙俊英,金寧德2012物理學(xué)報61 094701]

[8]Sun B,Wang EP,Ding Y 2011 Chinese Journal of Chemical Engineering 19 243

[9]Sun B,Wang E P,Zheng Y J 2011 Acta Phys.Sin.60 014701(in Chinese)[孫斌,王二朋,鄭永軍2011物理學(xué)報60 014701]

[10]Du M,Jin N D,Gao Z K 2012 Chem.Engin.Sci.82 144

[11]Manfredo G 2013 Int.J.Multiphase Flow 51 1

[12]van Ommen JR,Sasic S,van der Schaaf J,Gheorghiu S,Johnsson F,Coppens M 2011 Int.J.Multiphase Flow 37 403

[13]Barker RW,Klutke G,Hinich M J1993 ASMETrans.115 23

[14]Zhang Y,Wang S X,Li S H 1996 Acta Electron.Sin.24 109(in Chinese)[張嚴,王樹勛,李生紅1996電子學(xué)報24 109]

[15]Lutes L D,Chen D CK 1991 Int.J.Non-Linear Mech.6 893

[16]Hinich M J1994 Circuits Sys.Signal Proc.13 391

[17]Zhang W J,Huang Y J 2009 Machine Tool and Hydraulics 37 52(in Chinese)[張文俊,黃宜堅2009機床與液壓37 52]

[18]Zheng GB,Jin N D 2009 Acta Phys.Sin.58 4485(in Chinese)[鄭桂波,金寧德2009物理學(xué)報58 4485]