胡 純，韓明哲，鄭德智,3，彭 鵬

(1.北京航空航天大學(xué) 電子信息工程學(xué)院，北京 100191；2.北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院，北京 100191；3.北京航空航天大學(xué) 前沿科學(xué)技術(shù)創(chuàng)新研究院，北京 100191；4.重慶德新機(jī)器人檢測(cè)中心有限公司，重慶 400700)

流量測(cè)量、壓力測(cè)量與溫度測(cè)量被稱為工業(yè)測(cè)量領(lǐng)域的三大參數(shù)測(cè)量，其中以流量測(cè)量最為復(fù)雜、精度要求最高，其測(cè)量的經(jīng)濟(jì)性、安全性、可靠性對(duì)生產(chǎn)過(guò)程有著重要影響，在工業(yè)中具有重要地位[1]。流量傳感器分為體積流量傳感器與質(zhì)量流量傳感器。體積流量會(huì)受到壓力、溫度、黏度與密度等因素的影響，并且標(biāo)定過(guò)程煩瑣；質(zhì)量流量則不受上述因素的影響，并且適用于氣液兩相流，因此日益得到人們的青睞[2]。

目前針對(duì)靈敏度的理論較為成熟，對(duì)測(cè)量管振動(dòng)模型、集中質(zhì)量的影響、最佳檢測(cè)位置等的研究較為充分，針對(duì)U型管和直管都進(jìn)行了較多研究。文軍浩等[3]關(guān)注了科氏質(zhì)量流量傳感器中的非線性因素，并研究了其對(duì)U型管質(zhì)量流量傳感器的影響。樊剛等[4]就溫度對(duì)科氏質(zhì)量流量傳感器進(jìn)行了分析，并建立了流量傳感器溫度補(bǔ)償模型，對(duì)溫度補(bǔ)償進(jìn)行了仿真。曹勝?gòu)?qiáng)等[5]將污垢按集中質(zhì)量處理，通過(guò)微分方程和ANSYS仿真的方法分析了污垢對(duì)單直管科氏質(zhì)量流量傳感器靈敏度的影響。部分科氏質(zhì)量傳感器采用了在測(cè)量管兩端安裝波紋管的結(jié)構(gòu)，波紋管剛度遠(yuǎn)小于鋼管，因此此類測(cè)量管狀態(tài)與一般測(cè)量管不同，徐偉國(guó)等[6]和紀(jì)彩虹等[7]分別就兩段波紋管和三段波紋管的直管科氏質(zhì)量流量傳感器進(jìn)行了靈敏度計(jì)算，但尚沒(méi)有用于科氏質(zhì)量流量傳感器靈敏度分析的軟件。為彌補(bǔ)此空缺，本文設(shè)計(jì)了一款具有靈敏度分析、頻率計(jì)算、最佳檢測(cè)位置分析等功能的軟件，為科氏質(zhì)量流量傳感器的設(shè)計(jì)與開(kāi)發(fā)提供便利。

首先從科氏質(zhì)量流量傳感器的數(shù)學(xué)模型入手，對(duì)科氏質(zhì)量流量傳感器進(jìn)行了分析，建立了測(cè)量管的歐拉梁靜力模型，并分別計(jì)算了測(cè)量管在激振力與科氏力作用下的撓度曲線；建立了測(cè)量管的微分方程模型，并通過(guò)求解其振動(dòng)微分方程，計(jì)算出了測(cè)量管在激振力下的諧振頻率和科氏質(zhì)量流量傳感器的靈敏度；接下來(lái)分析了測(cè)量管的壓力損失，給出了在層流和湍流下的壓損計(jì)算公式；研究了溫度對(duì)科氏質(zhì)量流量傳感器的影響，對(duì)當(dāng)前測(cè)量管的常用材料進(jìn)行了調(diào)研，采用線性模型估計(jì)了溫度對(duì)彈性模量的影響，并用以對(duì)分析結(jié)果進(jìn)行補(bǔ)償。

在理論分析的基礎(chǔ)上，通過(guò)Matlab軟件的APP Designer模塊制作了直管科氏質(zhì)量流量傳感器的分析軟件，實(shí)現(xiàn)了計(jì)算靈敏度、頻率、壓損等數(shù)據(jù)的功能；并使用ANSYS Workbench進(jìn)行有限元分析；對(duì)其進(jìn)行了靜態(tài)結(jié)構(gòu)分析，分別得到其在激振力和科氏力作用下的靜撓度曲線，由此計(jì)算出傳感器的最佳安裝位置；進(jìn)行了瞬態(tài)動(dòng)力學(xué)分析，計(jì)算測(cè)量管在激振力作用下的振動(dòng)；進(jìn)行了模態(tài)分析得到其各階模態(tài)及諧振頻率，并分析了溫度對(duì)諧振頻率的影響；最后通過(guò)諧響應(yīng)分析得到測(cè)量管位移振幅曲線。有限元仿真的結(jié)果被用于驗(yàn)證模型與所設(shè)計(jì)軟件的計(jì)算結(jié)果，以確保所設(shè)計(jì)軟件的正確。

1 直管型科氏質(zhì)量流量傳感器的數(shù)學(xué)模型

為設(shè)計(jì)靈敏度分析軟件，首先要建立科氏質(zhì)量流量傳感器的數(shù)學(xué)模型。為此分別采用了靜力學(xué)和微分方程的方法進(jìn)行分析，得到了靈敏度等傳感器指標(biāo)的計(jì)算公式，為GUI軟件的設(shè)計(jì)提供支持。

1.1 直管型科氏質(zhì)量流量傳感器的靜力學(xué)模型

在科氏質(zhì)量流量傳感器工作時(shí)，其激振力頻率應(yīng)當(dāng)接近測(cè)量管的固有頻率，此時(shí)為簡(jiǎn)便分析過(guò)程，可以將測(cè)量管在恒定激振力下的靜撓度曲線近似為測(cè)量管撓度曲線進(jìn)行分析[8]。直管型測(cè)量管的雙端固定且長(zhǎng)細(xì)比大，符合雙端固支梁模型，可視為雙端固支梁進(jìn)行計(jì)算。

雙端固支梁為超靜定結(jié)構(gòu)，為對(duì)其進(jìn)行分析，需要解除多余的約束條件并采用等效的作用力和力矩建模。將激振力F=Ffsin(ωt+φ)視為恒定的力F，設(shè)測(cè)量管兩端約束力與約束力矩分別為Fa,Fb,Ma,Mb，由于測(cè)量管結(jié)構(gòu)對(duì)稱，因此有

(1)

Ma=Mb

(2)

為計(jì)算兩固支端的約束力矩，考慮其兩端和中心的轉(zhuǎn)角。由于兩端是固支端，其3個(gè)自由度都受到約束，其轉(zhuǎn)角和撓度均為0；又由于雙端固支梁為對(duì)稱結(jié)構(gòu)，故在作用于中點(diǎn)的激振力作用下，梁的形變同樣具有對(duì)稱性。因此梁的中點(diǎn)轉(zhuǎn)角也為0。

通過(guò)材料力學(xué)的方法分析梁的彎曲變形，通常借助梁的撓曲線近似微分方程：

(3)

對(duì)于集中力P作用下的測(cè)量管，其彎矩為

(4)

使用積分法求解梁的中點(diǎn)轉(zhuǎn)角，將式(3)從梁的端點(diǎn)到中點(diǎn)積分，得

(5)

(6)

因此，結(jié)合式(3)～式(6)，得到在激振力F作用下梁的撓度微分方程：

(7)

對(duì)其進(jìn)行積分得到梁的轉(zhuǎn)角方程和撓度方程：

(8)

(9)

激振力F作用下梁的靜撓度曲線如圖1所示。

圖1 激振力F作用下梁的靜撓度曲線

圖1中對(duì)橫坐標(biāo)x和撓度y都做了去量綱處理。

在工作過(guò)程中，測(cè)量管所受科氏力載荷為

dFk=2Qω(x)dx

(10)

式中，w(x)為角速度，即角度對(duì)時(shí)間的微分：

(11)

為對(duì)測(cè)量管上的科氏力載荷進(jìn)行分析，考慮在測(cè)量管上橫坐標(biāo)a處施加集中力dFk的情況。可以使用與分析激振力時(shí)相同的方法，將其視為恒定的力，解除其一邊約束，并加以等效的約束力F′(a)和約束力矩M′(a)進(jìn)行求解。再將求解出的撓度曲線關(guān)于dFk進(jìn)行積分，得到科里奧利力載荷作用下的撓度曲線。經(jīng)計(jì)算，前半段梁的撓度為

5.64×10-4x3L3+1.95×10-4x2L4)

(12)

后半段梁的撓度可由前半段關(guān)于中點(diǎn)對(duì)稱得到。其撓度曲線如圖2所示。

圖2 科里奧利力載荷作用下梁的靜撓度曲線

分析式(12)，可以得到在x=0.287L處，梁的撓度取得最大值。因此，最佳檢測(cè)位置在0.287L處和0.713L處。

1.2 直管型科氏質(zhì)量流量傳感器的微分方程模型

對(duì)于等截面歐拉梁振動(dòng)力學(xué)中已有成熟的理論，根據(jù)達(dá)朗貝爾原理可以列出其動(dòng)力學(xué)方程：

(13)

式中，X,Y分別為測(cè)量管橫縱坐標(biāo)；T為時(shí)間；V為待測(cè)流體流速；E為測(cè)量管的彈性模量；I為測(cè)量管的截面慣性矩；m1,m2,m3分別為科氏質(zhì)量流量傳感器的檢測(cè)器和激振器質(zhì)量，其中m1,m3為檢測(cè)器，m2為激振器，兩檢測(cè)器質(zhì)量相等，安裝在關(guān)于中點(diǎn)對(duì)稱的位置上；mp,mf分別為測(cè)量管線密度和待測(cè)流體的線密度。

為計(jì)算諧振頻率，需要對(duì)該微分方程進(jìn)行求解。首先為方便分析，將其變形為無(wú)量綱形式：

(14)

式中，x,y,t,σk,v,β分別為無(wú)量綱的橫坐標(biāo)、縱坐標(biāo)、時(shí)間、集中質(zhì)量、流速、流體密度。

記測(cè)量管震動(dòng)角頻率為Ω，測(cè)量管無(wú)量綱振動(dòng)角頻率為ω：

對(duì)于此類微分方程，已經(jīng)有成熟的求解理論，其解為[9]

y(x,t)=RΦ(x)eiωt

(15)

式中，R為常系數(shù)；Φ(x)為振型函數(shù)；i為單位虛數(shù)。對(duì)測(cè)量管的振型函數(shù)，可以各階振型函數(shù)構(gòu)成的無(wú)窮級(jí)數(shù)來(lái)表述：

(16)

其振型函數(shù)為[10]

φr(x)=cosh(krx)-cos(krx)-λr(sinh(krx)-sin(krx))

(17)

本文關(guān)注的是梁在激振力作用下的振動(dòng)和在科里奧利力作用下的振動(dòng)，分別對(duì)應(yīng)一階和二階振型，振型曲線如圖3與圖4所示，其參數(shù)k和λ分別為：k1=4.730,k2=7.853,λ1=0.9825,λ2=1.001。

圖3 一階振型函數(shù)曲線

圖4 二階振型函數(shù)曲線

注意到微分方程解的一、二階振型函數(shù)，與前文所推導(dǎo)的測(cè)量管在激振力與科氏力作用下的撓度方程具有相同的物理意義，由于在分析過(guò)程中采用了不同的方法去量綱以簡(jiǎn)化分析，因此二者應(yīng)當(dāng)具有線性關(guān)系。在此采用了相關(guān)系數(shù)的方法進(jìn)行驗(yàn)證。

首先對(duì)計(jì)算出的函數(shù)進(jìn)行采樣，將其離散化。采樣間隔設(shè)為0.001，取1000個(gè)采樣點(diǎn)，然后分別計(jì)算一階振型與激振力撓度、二階振型與科氏力撓度的皮爾遜相關(guān)系數(shù)，分別為r1=0.999052,r2=0.999849，這說(shuō)明二者高度相關(guān)，驗(yàn)證了上述分析結(jié)果。二階振型函數(shù)在x=0.290處取得最大值，與2.1節(jié)中分析得到的最大值位置0.287相接近。歸一化后的對(duì)比圖如圖5、圖6所示。

圖5 一階振型與激振力撓度對(duì)比圖

圖6 二階振型與科氏力撓度對(duì)比圖

由圖5和圖6可以看出，兩種計(jì)算方法所得到的振型與撓度曲線高度相似，此結(jié)果對(duì)上述分析的正確性提供了支持。

使用伽遼金法求解微分方程(14)，并取其前2階振型作為該微分方程的近似解，解得：

(18)

(19)

在科氏質(zhì)量流量傳感器工作時(shí)，測(cè)量管同時(shí)存在激振力引起的振動(dòng)和科氏力引起的振動(dòng)，這使得測(cè)量點(diǎn)存在振動(dòng)的相位差。因?yàn)榭剖腺|(zhì)量流量傳感器通過(guò)測(cè)量此相位差以得到待測(cè)流體的質(zhì)量流量，所以將科氏質(zhì)量流量傳感器的靈敏度定義為單位流量的相位差，即：

(20)

由式(19),可以得到測(cè)量管的振動(dòng)函數(shù)：

(21)

其相位差Δφ為

(22)

進(jìn)一步計(jì)算出科氏質(zhì)量流量傳感器的靈敏度K：

(23)

其中，

1.3 壓損計(jì)算

壓損按照產(chǎn)生的機(jī)理可以分為沿程壓力損失和局部壓力損失。對(duì)于直管型測(cè)量管，不存在閥口、彎管等通流截面變化或流動(dòng)方向改變的結(jié)構(gòu)，因此局部壓力損失可以忽略不計(jì)。

圓管層流是流體運(yùn)動(dòng)中較為簡(jiǎn)單的一種情況，也是能夠得出流速分布及壓力損失解析解的為數(shù)不多的情況之一。經(jīng)計(jì)算，在層流中沿程壓力損失系數(shù)為[11]

(24)

湍流運(yùn)動(dòng)的情況較為復(fù)雜，沿程壓力損失系數(shù)的計(jì)算往往依賴經(jīng)驗(yàn)公式或半經(jīng)驗(yàn)公式。1913年德國(guó)水力學(xué)家布拉修斯(Blasius)提出了計(jì)算湍流沿程壓力損失系數(shù)的經(jīng)驗(yàn)公式:

(25)

布拉修斯公式適用于Re<105的紊流光滑區(qū)，計(jì)算簡(jiǎn)單方便，且具有較高精度。

沿程壓力損失的理論計(jì)算公式為

(26)

式中，l為測(cè)量管長(zhǎng)度；ρ為測(cè)量管密度；v為流體流速，且滿足

(27)

在層流和湍流下測(cè)量管的壓損計(jì)算公式為

(28)

式中，Q為流體質(zhì)量流量(kg/s)；η為流體動(dòng)力黏度(Pa·s)；ρ為流體密度(kg/m3)。

1.4 溫度對(duì)科氏質(zhì)量流量傳感器的影響

科氏質(zhì)量流量傳感器在工業(yè)領(lǐng)域中的應(yīng)用逐步廣泛，其所面臨的工作環(huán)境也多種多樣。因此，在各種環(huán)境下都能保持較高的精度已成為科氏質(zhì)量流量傳感器的迫切需要。為此，科氏質(zhì)量流量傳感器通常會(huì)設(shè)置溫度傳感器，測(cè)量被測(cè)介質(zhì)的溫度，通過(guò)溫度補(bǔ)償?shù)姆椒ǖ窒麥囟葞?lái)的影響，以保證測(cè)量數(shù)據(jù)的精度。

科氏質(zhì)量流量傳感器最常用的材料是316L不銹鋼，除此之外哈氏合金c-22也被用于科氏質(zhì)量流量傳感器，在直管科氏質(zhì)量流量傳感器上，有時(shí)以硅合金作為測(cè)量管材料。筆者選取了最常見(jiàn)的這3種金屬作為備選材料，對(duì)其屬性進(jìn)行了調(diào)研，結(jié)果如表1～表4所示。

表1 316L不銹鋼(Acerinox公司)屬性(自20 ℃起)

表2 哈氏合金c-22(Hastelloy公司)屬性(自20 ℃起)

表3 鈦合金Ti 6Al-4V(NeoNickel公司)屬性(自20 ℃起)

表4 常見(jiàn)材料密度

上述各材料的平均線膨脹系數(shù)都在10-5量級(jí)，因此在本研究中，可以忽略溫度對(duì)密度和測(cè)量管結(jié)構(gòu)參數(shù)的影響，溫度的影響主要體現(xiàn)在彈性模量上。對(duì)溫度與彈性模量進(jìn)行線性回歸分析，其結(jié)果如圖7～圖9所示。

圖8 哈氏合金c-22彈性模量與溫度的線性回歸模型

圖9 鈦合金Ti 6Al-4V彈性模量與溫度的線性回歸模型

上述回歸分析相關(guān)系數(shù)R2均大于0.997，下面使用F檢驗(yàn)法進(jìn)行回歸方程顯著性檢驗(yàn)。

相關(guān)系數(shù)R2與統(tǒng)計(jì)量F有如下關(guān)系：

(29)

由式(29)計(jì)算得，上述回歸分析統(tǒng)計(jì)量F分別為

F1=5795.89，F(xiàn)2=4340.53，F(xiàn)3=2896.18

對(duì)于316L不銹鋼、哈氏合金c-22和鈦合金Ti 6Al-4V的數(shù)據(jù)，其樣本容量分別為n=6,5,4。對(duì)F檢驗(yàn)進(jìn)行查表，得

F0.01(1,4)=21.20

F0.01(1,3)=34.12

F0.01(1,2)=98.49

因此上述回歸均在0.01的水平上高度顯著，故接受回歸分析的結(jié)果，并在靈敏度計(jì)算中采用上述模型進(jìn)行彈性模量的計(jì)算。在軟件設(shè)計(jì)中還加入了自定義材料的功能，支持直接輸入材料的彈性模量和密度，以適應(yīng)使用其他材料或用戶具有準(zhǔn)確的材料屬性的情況。

2 GUI軟件設(shè)計(jì)

本GUI軟件基于Matlab的APP designer進(jìn)行設(shè)計(jì)，由參數(shù)設(shè)置、材料屬性設(shè)置、集中質(zhì)量設(shè)置、最佳檢測(cè)位置分析、指定位置靈敏度分析、壓損分析和曲線繪制7個(gè)面板構(gòu)成，根據(jù)以上各節(jié)所推導(dǎo)計(jì)算的公式，完成了研究目標(biāo)中所提出的靈敏度分析、諧振頻率分析、壓損分析等目標(biāo)，實(shí)現(xiàn)了輸入測(cè)量管結(jié)構(gòu)參數(shù)、材料屬性參數(shù)、傳感器集中質(zhì)量等數(shù)據(jù)，進(jìn)行分析計(jì)算最佳檢測(cè)位置、靈敏度、頻率、壓損等數(shù)據(jù)的功能，并繪制相關(guān)曲線圖。軟件還具有根據(jù)溫度對(duì)材料屬性進(jìn)行補(bǔ)償?shù)墓δ?。設(shè)計(jì)的GUI軟件功能架構(gòu)圖如圖10所示，GUI軟件界面如圖11所示。

圖10 軟件功能架構(gòu)圖

圖11 GUI軟件界面

Matlab的APP Designer提供了豐富而強(qiáng)大的回調(diào)函數(shù)，極大地便利了軟件的設(shè)計(jì)。本軟件充分利用了APP Designer的回調(diào)函數(shù)，基本上全部功能都通過(guò)回調(diào)函數(shù)實(shí)現(xiàn)。在程序運(yùn)行的初期，先進(jìn)行GUI界面的初始化，生成GUI界面、初始化數(shù)據(jù)并繪制表格和坐標(biāo)軸；然后等待用戶對(duì)界面進(jìn)行操作，根據(jù)用戶的操作執(zhí)行相對(duì)應(yīng)的回調(diào)函數(shù)。程序流程如圖12所示。

3 科氏質(zhì)量流量傳感器的仿真分析

使用ANSYS Workbench對(duì)測(cè)量管進(jìn)行了有限元仿真。首先建立了直管型測(cè)量管的模型，然后分別對(duì)其進(jìn)行了靜態(tài)結(jié)構(gòu)分析、瞬態(tài)動(dòng)力學(xué)分析、模態(tài)分析和諧響應(yīng)分析，以驗(yàn)證理論計(jì)算結(jié)果。所使用的測(cè)量管結(jié)構(gòu)為單圓直管，并忽略了集中質(zhì)量，其結(jié)構(gòu)參數(shù)如表5所示。

表5 測(cè)量管結(jié)構(gòu)尺寸

3.1 靜態(tài)結(jié)構(gòu)分析

于測(cè)量管兩端施加固定約束，施加作用于測(cè)量管外表面上、作用位置為測(cè)量管中央的正上方、方向豎直向下的激振力F，以進(jìn)行激振力作用下的靜撓度仿真；為計(jì)算科氏力作用下的撓度，采用了分段施加載荷的方式，將測(cè)量管分為20段，在每段下施加相應(yīng)的科氏力載荷，以模擬實(shí)際的科氏力載荷分布。在網(wǎng)格劃分方面，采用了由ANSYS軟件智能劃分的自由網(wǎng)格。其仿真結(jié)果如圖13、圖14所示。

圖13 激振力作用下的靜態(tài)結(jié)構(gòu)仿真

圖14 科氏力作用下的靜態(tài)結(jié)構(gòu)仿真

通過(guò)將激振力作用下測(cè)量管的靜撓度曲線和以靜力學(xué)模型計(jì)算出撓度曲線、微分方程的振型函數(shù)進(jìn)行對(duì)比，驗(yàn)證了前文所進(jìn)行分析的正確性。

3.2 模態(tài)分析和諧響應(yīng)分析

同樣于兩端施加固定約束，無(wú)預(yù)應(yīng)力輸入，取20 ℃常溫下的材料屬性，同樣采用自由網(wǎng)格劃分。激振力作用下的振動(dòng)為一、二階模態(tài)，科氏力作用下的振動(dòng)為三、四階模態(tài)，因此選取前六階模態(tài)進(jìn)行模態(tài)分析。得到其前六階模態(tài)諧振頻率如表6所示。

表6 模態(tài)分析諧振頻率

直管型測(cè)量管結(jié)構(gòu)關(guān)于其中軸線對(duì)稱，因此其一二階、三四階、五六階的諧振頻率相同，其振型函數(shù)在空間上正交，且在繞測(cè)量管中軸線旋轉(zhuǎn)90°后重合，可以看作是一個(gè)二自由度振動(dòng)在空間上的分解。

上述仿真得到測(cè)量管的諧振頻率為131.53 Hz，運(yùn)行所設(shè)計(jì)的GUI軟件，輸入表1所示的材料屬性數(shù)據(jù)和表5所示的結(jié)構(gòu)尺寸參數(shù)，選擇測(cè)量管內(nèi)流體為空，將溫度設(shè)置為20 ℃，計(jì)算得到的諧振頻率為132.2 Hz，其相對(duì)誤差為+0.5094%。

施加在測(cè)量管外表面中央正上方，方向向下的簡(jiǎn)諧激振力作為激勵(lì)，進(jìn)行諧響應(yīng)分析結(jié)果見(jiàn)圖15。

圖15 幅值響應(yīng)

通過(guò)修改模態(tài)分析中的溫度條件，對(duì)軟件的溫度補(bǔ)償功能進(jìn)行驗(yàn)證。計(jì)算與仿真結(jié)果如表7所示。

表7 諧振頻率仿真驗(yàn)證數(shù)據(jù)

上述兩組數(shù)據(jù)關(guān)于溫度的線性回歸分析都在0.01的水平上顯著，此結(jié)果支持了溫度對(duì)諧振頻率具有線性影響的結(jié)論，驗(yàn)證了前文中進(jìn)行的溫度對(duì)科氏質(zhì)量流量傳感器的影響的正確性。

4 結(jié)束語(yǔ)

設(shè)計(jì)了直管型科氏質(zhì)量流量傳感器的靈敏度分析軟件，通過(guò)理論分析得到的公式進(jìn)行靈敏度、頻率、壓損的計(jì)算，并具有對(duì)溫度進(jìn)行補(bǔ)償?shù)墓δ?。并使用ANSYS進(jìn)行有限元仿真，對(duì)所設(shè)計(jì)軟件的計(jì)算結(jié)果進(jìn)行了驗(yàn)證，表明其具有較高的可靠性。