1種通過聲帶空氣動力學模型求解發(fā)聲閾值壓力的方法

張莉麗，王 琰，陳莉媛，肖仲喆，張曉俊，陶 智

(蘇州大學 光電科學與工程學院，江蘇 蘇州 215000)

嗓音的產(chǎn)生來源于肺部氣流引起的聲帶振動，其中發(fā)聲閾值壓力是嗓音功能的1個重要指標.當肺部氣壓達到閾值壓力時，氣流向聲帶傳遞足夠的能量克服聲帶組織耗散所帶來的損失，使聲帶產(chǎn)生大幅度振蕩，再經(jīng)過口腔輻射最終產(chǎn)生嗓音輸出.從功能上講，發(fā)聲閾值壓力是能夠產(chǎn)生聲帶振動所需要的最低能量指數(shù)，在主觀感受上表現(xiàn)為發(fā)聲的力度和發(fā)聲的容易程度.因此發(fā)聲閾值壓力的研究對于更好地理解嗓音產(chǎn)生機理具有重要意義.

發(fā)聲閾值壓力最早由Titze[1]于1988年提出.在臨床上，其可作為聲帶功能健康的指標，也可作為量化嗓音情感和聲帶疲勞度的客觀指標，反映了聲帶黏性程度的變化[2].研究人員利用物理模型[3-4]、離體喉[5]和真人測量[6]等方式來探究與發(fā)聲閾值壓力相關(guān)的理論和臨床意義.Titze[1]根據(jù)二質(zhì)量模型在發(fā)聲時的聲帶表面壓力分布推導出發(fā)聲閾值壓力與預發(fā)聲聲門直徑、跨聲門壓系數(shù)成正比，與聲帶長度成反比.此外，為了測試閾值壓力與聲門半徑的關(guān)系，Titze等[7]采用樹脂玻璃搭建半喉模型，在這個模型中通過螺旋調(diào)節(jié)聲門半徑，并在聲帶表面覆蓋1層薄硅膠膜，硅膠膜中可以循環(huán)不同粘度的液體來模擬聲帶的不同種情況，實驗證實了閾值壓力與預發(fā)聲聲門直徑成正比，與聲帶長度成反比.但實驗結(jié)果在聲門直徑等于零時，測得的閾值壓力并不為零，與通過二質(zhì)量模型得到的結(jié)果不相符.為了探究這一原因，Lucero[8]考慮在聲門直徑較小時增加黏性效應(yīng)，對發(fā)聲閾值壓力表達式修正,通過泊肅葉效應(yīng)求解聲門直徑較小時聲帶內(nèi)的黏性效應(yīng)，這種修正與Titze等[7]得到的實驗結(jié)果基本吻合.后來，Chan等[9]考慮將粘膜波和聲道的慣性阻抗作為自激振動的主要能量傳遞機制，實驗得到的結(jié)果與閾值壓力表達式和聲門半徑之間的關(guān)系在圖像的上升和下降趨勢上保持一致，但斜率上存在差異，推測可能是由于未考慮聲門直徑比較小時聲帶內(nèi)的黏性效應(yīng)所導致的，但與Lucero[8]獲得的表達式在聲門直徑較小的情況仍然存在很大的差異.Fulcher等[10]根據(jù)入口損失系數(shù)和出口系數(shù)計算聲帶表面壓力分布，并且將計算得到的值與實驗結(jié)果進行對比，得出了這2種系數(shù)與聲門直徑和跨聲門壓之間的關(guān)系.近年來，隨著數(shù)值方法的發(fā)展和計算能力的提高，數(shù)值計算模型被廣泛地應(yīng)用于發(fā)聲過程中聲門內(nèi)流場變化和聲帶振動的分析[11-12].數(shù)值計算模型通過離散化求解流體力學中的Navier-Stokes控制方程，考慮肺部產(chǎn)生的動態(tài)激勵源對發(fā)聲的作用，將氣流運動和聲帶振動相結(jié)合，更符合嗓音的實際產(chǎn)生過程[13].

現(xiàn)有的發(fā)聲閾值壓力推導出的結(jié)果并不能完全預測發(fā)聲閾值壓力的表達形式，與實驗結(jié)果均有較大差異，因此有必要從發(fā)聲原理角度研究閾值壓力的具體表現(xiàn)形式.本文從聲帶振動特性的角度出發(fā)，建立聲帶的空氣動力學模型對聲門內(nèi)流場進行數(shù)值求解，得出聲帶表面壓力和聲門內(nèi)氣流的體積流量，進而推導出聲門直徑與入口損失系數(shù)的關(guān)系，對發(fā)聲閾值壓力的表達式進行修正，并將求解出的發(fā)聲閾值壓力表達式與前人的實驗結(jié)果擬合，驗證修正后發(fā)聲閾值壓力表達式的準確性.

1 聲帶空氣動力學建模的閾值壓力修正

1.1 發(fā)聲閾值壓力表達式的求解缺陷

發(fā)聲閾值壓力即為能使聲帶產(chǎn)生振蕩的最小肺部氣壓[1]，根據(jù)二質(zhì)量塊運動的穩(wěn)定性，求得發(fā)聲閾值壓力表達式為

(1)

Chan等[9]于2006年建立了1個半喉的硬塑模型，將閾值壓力與聲帶振動部分的線性黏彈性剪切特性和上聲道的慣性電抗聯(lián)系起來，并在物理模型上進行了額外的實驗，將粘彈性生物材料植入人工聲道內(nèi)，包括透明質(zhì)酸(Hyaluronic Acid, HA)、纖維連接蛋白(Hyaluronic Acid with FibroNectin, HA with FN)，對聲帶組織粘彈性在發(fā)聲閾值壓力上的效應(yīng)進行評估.實驗過程中，聲門下壓力從0開始逐漸增大，直至聲帶出現(xiàn)穩(wěn)定的自激振蕩.對比不同濃度的透明質(zhì)酸和纖維連接蛋白添加至聲帶表皮得到的發(fā)聲閾值壓力結(jié)果，發(fā)現(xiàn)當預發(fā)聲聲門直徑很小時，發(fā)聲閾值壓力并不為0，而式(1)的結(jié)果與實驗得到的結(jié)果在預發(fā)聲聲門直徑等于0時不一致.

上述實驗結(jié)果說明僅從發(fā)聲驅(qū)動力的角度并不能給發(fā)聲閾值壓力以很好的修正.而聲帶的空氣動力學模型可以依據(jù)聲帶和聲門的精確尺寸制成，并可以精確地控制聲帶的跨聲門壓和聲帶運動，通過聲帶氣動仿真模型研究聲帶內(nèi)的壓力變化可以對發(fā)聲閾值壓力的表達式進行修正.

1.2 聲帶空氣動力學模型和數(shù)值求解過程

根據(jù)聲帶的生理組織結(jié)構(gòu)，以聲門為基本氣流通道建立聲帶的空氣動力學模型，利用流體動力學計算軟件Fluent對聲門管內(nèi)流場進行數(shù)值計算，得出聲帶表面受力和聲門內(nèi)氣體的體積流量，以探究預發(fā)聲聲門直徑的改變對聲帶受力和聲門內(nèi)氣流分布的影響.

根據(jù)Scherer等[14]建立的M5聲帶模型，聲門入口的形狀采用以聲門角度為變量的半徑來定義，聲門出入口曲線用直線連接起來構(gòu)成聲帶表面.由于在剛性模型中壓力剖面的經(jīng)驗性結(jié)果表明聲門形狀不會在前后逐漸變尖，相對較尖銳的聲門入口與聲帶上皮的實際情況不太一致，因此將其改為圓形入口，以更符合聲帶的實際生理結(jié)構(gòu).

圖1 聲帶的空氣動力學模型圖Fig.1 Aerodynamic model of the vocal cords

本文建立的聲帶的空氣動力學模型如圖1所示，聲門入口半徑Rin=0.15cm，出口半徑Rout=0.108cm，聲帶軸向長度lg=0.308cm，聲帶寬度Lg=0.5cm，聲門入口偏轉(zhuǎn)角φin=50°，出口偏轉(zhuǎn)角φout=90°，入口部分的軸向長度Lin=0.2cm，出口部分的軸向長度Lout=1.5cm.氣流出口壓力設(shè)為標準大氣壓，入口壓力根據(jù)實驗需求設(shè)置不同的數(shù)值.

流體在運動過程中遵循機械運動普遍適用的守恒定律: 質(zhì)量守恒、動量守恒、能量守恒.由于本模型中氣流雷諾數(shù)小于4000，可認為流體為均質(zhì)、不可壓縮.在忽略重力的情況下，聲門內(nèi)流場可由Navier-Stokes方程表示:

(2)

同時對于定常、不可壓縮流體，密度ρ為常數(shù)，其質(zhì)量守恒方程(連續(xù)性方程)表示為

(3)

其中:vx，vy，vz分別為x，y，z方向上的速度分量；ρ為空氣密度；P為空氣壓力；μ為空氣的黏性系數(shù).

在求解氣流控制方程時，將計算區(qū)域劃分為一系列不重復的控制體積，并使每個網(wǎng)格點周圍有1個控制體積，將待解的微分方程對每個控制體積積分，得到1組離散方程組.在對模型進行網(wǎng)格劃分離散化后，利用SIMPLE(Semi-Implicit Method for Pressure Linked Equations)算法迭代求解壓力-速度關(guān)聯(lián)方程.首先給定單元面壓力場的估算值，并使用此估算值求解動量方程得到該單元面的速度場，然后通過求解壓力修正方程得到壓力的修正值，基于壓力修正值對壓力場和速度場進行修正并檢查結(jié)果是否收斂，若不收斂，則以修正后的壓力作為新的估算值繼續(xù)迭代計算，直至滿足收斂條件.

1.3 發(fā)聲閾值壓力表達式的修正

(4)

式中:ρ=1.3×10-3g/cm3，為空氣密度；v為平均速度，與體積流量U和橫截面積A之間的關(guān)系為U=vA.

但是由于聲門內(nèi)黏性損失力的存在，實際的聲帶壓力測量實驗中Pin≠Pout[15].設(shè)聲門入口處黏滯效應(yīng)的損失系數(shù)為kloss，則肺部氣壓與聲門入口之間的壓力差為

(5)

式中:U是聲門內(nèi)氣流的體積流量；lg是聲帶長度；d為聲門直徑；lg與d的乘積即為聲門的橫截面積A.

(6)

(7)

利用聲帶空氣動力學模型對聲門內(nèi)流場進行數(shù)值求解，可以得到聲帶兩側(cè)受力和聲門內(nèi)氣流的體積流量，同時設(shè)置不同的聲門直徑參數(shù)，并對其施加不同的聲門下壓，代入式(7)中，可得聲門入口損失系數(shù)和聲門直徑之間的關(guān)系.對聲帶空氣動力學模型中的d分別設(shè)置為0.005，0.0075，0.010，0.020，0.040，0.080和0.160cm，每種聲門直徑下均采用300，500，1000，1500Pa的聲門下壓作為聲門入口的邊界條件.圖2(a)展示了不同聲門下壓時聲門直徑與入口損失系數(shù)的曲線關(guān)系.從圖2(a)中可以看出，對于較低的聲門下壓，聲門入口損失系數(shù)隨著聲門直徑的減小有顯著地增加，并且聲門入口損失系數(shù)與預發(fā)聲聲門直徑呈現(xiàn)出1種反比例函數(shù)的趨勢，表明聲門入口損失系數(shù)為非恒定數(shù)值.將入口損失系數(shù)與聲門入口直徑相乘，可以得到圖2(b)的結(jié)果.對其中的曲線進行多項式擬合，得到聲門入口損失系數(shù)與預發(fā)聲聲門直徑之間的關(guān)系表達式，即

(8)

其中α，β，ε為由實驗決定的常數(shù).由于入口損失系數(shù)僅與聲帶形狀相關(guān)[10]，α值的選取取決于具體的聲帶模型尺寸.此外，對于式(7)，設(shè)置不同的聲門直徑和聲門下壓，得出不同聲門下壓時入口損失系數(shù)和聲門直徑之間的曲線關(guān)系，在圖2中，不同的聲門下壓得到的曲線關(guān)系基本是一致的，說明入口損失系數(shù)與聲門直徑之間的關(guān)系不取決于聲門下壓的改變.

圖2 不同聲門下壓時入口損失系數(shù)和聲門直徑之間的關(guān)系Fig.2 The relationship between entrance loss coefficient and glottal diameter under different glottal pressures

將式(8)代入式(1)中，得到發(fā)聲閾值表達式.在計算過程中，為了使參數(shù)達到最小，將其進行簡化，令ε=0，最終修正后的閾值壓力表達式為

(9)

通過這種形式可以直接確定發(fā)聲閾值壓力表達式的斜率和截距，將閾值壓力表達式的斜率和截距與聲帶模型本身配置的參數(shù)建立聯(lián)系，解決了預發(fā)聲聲門直徑等于零時而閾值壓力不為零的問題.

2 發(fā)聲閾值壓力表達式修正的實驗驗證

為了驗證修正后發(fā)聲閾值壓力表達式的有效性，將求解出的發(fā)聲閾值壓力的表達式與前人的實驗結(jié)果進行擬合.Chan等[9]搭建了在硅膠膜下植入生物力學材料的聲帶物理模型，通過在聲帶模型中加入不同濃度的透明質(zhì)酸(HA)來改變聲帶黏度，以測得不同聲門直徑下的發(fā)聲閾值壓力.將修正后的發(fā)聲閾值壓力的表達式和Chan等的實驗數(shù)據(jù)進行比較時，首先要確定入口損失參數(shù)α和β.Chan等的實驗是基于半喉模型而非對稱的M5模型，但聲門具有較大寬度時半喉聲帶與對稱聲帶的入口損失系數(shù)沒有太大的區(qū)別，因此α值與Chan等的實驗中的相同，考慮Chan等使用的實驗裝置尺寸，當α值為0.385cm時，可以較好地擬合出0.01%透明質(zhì)酸(P-=12880g/s2)填充聲帶時的情況(其中聲帶寬度Lg=2.22cm，聲帶垂直厚度T=1.11cm，下質(zhì)量塊厚度d1=5T/6)，同時為簡化計算，將β值設(shè)為1，結(jié)果如圖3所示，其中Eq.為式(9).

圖4是在0.01%濃度的透明質(zhì)酸中加入纖連蛋白(P-=17600g/s2)，然后填充聲帶所得的實驗結(jié)果.可以發(fā)現(xiàn)在相同的聲門寬度下，纖連蛋白的填充會使聲帶發(fā)聲的閾值壓力增大，如圖4中虛線所示.與式(9)保持一致，增大P-需要更大的截距和更陡的斜率.

圖3 2種濃度的透明質(zhì)酸中發(fā)聲閾值 壓力與聲門半徑之間的關(guān)系Fig.3 The relationship between the phonation threshold pressure and glottal radius in two concentrations of hyaluronic acid

圖4 加入纖連蛋白的0.01%濃度的透明質(zhì)酸中發(fā)聲 閾值壓力與聲門半徑之間的關(guān)系Fig.4 The relationship between the phonation threshold pressure and glottal radius in 0.01% concentration of hyaluronic acid added to the fibronectio

表1是濃度為0.01%，0.1%的透明質(zhì)酸和加入纖連蛋白的0.01%濃度的透明質(zhì)酸中閾值壓力的預測誤差.表1中，在所有閾值壓力的預測中，誤差最小可以達到0.36%，最大不超過30%.上述實驗結(jié)果可以驗證修正后的閾值壓力表達式能夠較準確地預測聲帶發(fā)聲時的閾值壓力，同時能夠解決當聲門半徑為零，發(fā)聲閾值壓力不為零的問題.

表1 聲帶模型中加入不同填充物時閾值壓力計算值與 測量值的誤差Tab.1 The error between the calculated value and the measured value of threshold pressure when dif- ferent fillings are added to the vocal cord model

3 結(jié) 語

為了對Titze提出的閾值壓力表達式與其實驗結(jié)果不相符的問題進行修正，本文從發(fā)聲原理的角度建立聲帶的空氣動力學模型，考慮聲門入口損失系數(shù)與發(fā)聲時的聲門直徑的關(guān)系，修正閾值壓力的表達式，解決在聲門直徑較小時發(fā)聲閾值壓力不為零的問題.將修正后的閾值壓力表達式與前人的實驗結(jié)果進行對比，說明本文算法可以較準確地預測發(fā)聲閾值壓力的變化趨勢，并且能夠滿足閾值壓力在聲門直徑較小時不為零的情況，更符合實際上的嗓音的發(fā)聲條件.本文主要考慮了聲門半徑與聲門入口損失系數(shù)之間的變化關(guān)系，進一步的研究將考慮聲帶厚度、下質(zhì)量厚度等相關(guān)因素對聲門入口損失系數(shù)的影響.