許嘉誠(chéng)

(上海交通大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 上海 200240)

0 引言

血壓在臨床上被認(rèn)為是診斷和預(yù)測(cè)心血管疾病的重要參考指標(biāo)[1-4].在測(cè)量血壓的各種技術(shù)中，示波法是目前主流方法之一.示波法最早是由Marey[5]提出的，目前的示波算法主要分為基于振幅的方法和基于斜率的方法.基于振幅的方法通常采用最大振幅算法(MAA).目前人們普遍認(rèn)為，這種方法可以準(zhǔn)確估計(jì)平均血壓，相關(guān)文獻(xiàn)利用動(dòng)脈順應(yīng)性來(lái)解釋[6-9].研究人員[10-13]通常是利用特征比值來(lái)獲得收縮壓和舒張壓，但由于特征比值是利用經(jīng)驗(yàn)獲得并且往往固定不變，因此準(zhǔn)確性不高，誤差較大.基于斜率的方法[14-15]是通過(guò)示波波形振幅包絡(luò)線的斜率獲得收縮壓和舒張壓的方法.這種方法在干凈的數(shù)據(jù)中表現(xiàn)良好，但在實(shí)際的“臟”數(shù)據(jù)中，噪聲的影響極易產(chǎn)生較大的誤差[14].

事實(shí)上，要想提高示波算法的精度，理論研究是必不可少的.由于理論分析的重要性，越來(lái)越多的學(xué)者投入到示波算法理論研究中.例如，Babbs[16]將動(dòng)脈壁模型與袖帶模型相結(jié)合，以模擬袖帶壓力振蕩，并通過(guò)回歸算法來(lái)提取收縮壓和舒張壓.由于組織的力學(xué)性能對(duì)應(yīng)力傳遞效率有著重要的影響，一些學(xué)者在前人研究的基礎(chǔ)上增加組織的影響[17-19].Ursino等[17]提出袖帶-手臂-動(dòng)脈系統(tǒng)的數(shù)學(xué)模型，以研究生物力學(xué)因素對(duì)測(cè)量的影響.梁等[19]構(gòu)建了與骨結(jié)合的三維袖帶-手臂-動(dòng)脈系統(tǒng)，獲得其應(yīng)力分布.

基于前人的研究，本文將對(duì)袖帶-組織-動(dòng)脈系統(tǒng)做相應(yīng)的數(shù)值和理論分析.本文的結(jié)構(gòu)如下：第1節(jié)介紹數(shù)學(xué)模型；第2節(jié)介紹基于虛功原理和基于有限元法的數(shù)值算法；第3節(jié)是量綱分析和單元分析；第4節(jié)是基于平移不變性的理論驗(yàn)證；最后一節(jié)是未來(lái)展望.

1 模型建立

在建立袖帶-組織-動(dòng)脈系統(tǒng)的數(shù)學(xué)模型之前，我們首先提出適當(dāng)?shù)募僭O(shè)，以簡(jiǎn)化模型并了解系統(tǒng)的機(jī)理.在這個(gè)問(wèn)題中，系統(tǒng)是軸對(duì)稱的，即動(dòng)脈、組織和袖帶是軸對(duì)稱的，具有相同的軸(參見(jiàn)圖1(a)).脂肪、肌肉和其他軟組織雖然具有不同的力學(xué)性質(zhì)，如彈性模量、泊松比等，但它們可以看作是具有相同力學(xué)性質(zhì)的各向同性均勻彈性體.因此，彈性模量和泊松比保持不變.此外，忽略了骨骼和身體力量的影響.

(a) (b)圖1 (a)袖帶-組織-動(dòng)脈系統(tǒng)結(jié)構(gòu)圖；(b)計(jì)算區(qū)域Fig.1 (a)Geometry of cuff-tissue-artery system;(b)The computational region

1.1 袖帶模型

在袖帶充氣和放氣的過(guò)程中，袖帶的外壁會(huì)發(fā)生輕微的變形.因此，外壁視為剛性.在袖帶變形過(guò)程中，袖帶內(nèi)壁緊貼組織，因此認(rèn)為袖帶的內(nèi)壁順應(yīng)性足夠大.袖帶壓力始終與外部組織的應(yīng)力相同.袖帶內(nèi)的空氣被認(rèn)為是理想氣體.此外，袖帶充氣和放氣的過(guò)程都被視為絕熱過(guò)程：

(1)

式中：Pc是袖帶壓強(qiáng)；Patm是標(biāo)準(zhǔn)大氣壓；Vc是袖帶體積；參數(shù)k=1.4，代表絕熱過(guò)程；Q代表在特定時(shí)間袖帶內(nèi)的剩余空氣量.

1.2 動(dòng)脈壁模型

實(shí)驗(yàn)發(fā)現(xiàn)，動(dòng)脈壁的橫截面積隨著血壓的升高而增大，直至達(dá)到最大橫截面積.然而，隨著跨壁壓力的增加，橫截面積的膨脹先快后慢.也就是說(shuō)，當(dāng)動(dòng)脈擴(kuò)張時(shí)，隨著跨壁壓力的增大，動(dòng)脈壁順應(yīng)性迅速降低.此外，當(dāng)血管壁塌陷時(shí)，跨壁壓力變小，血管順應(yīng)性也顯著下降.根據(jù)這些實(shí)驗(yàn)結(jié)果，文獻(xiàn)[20-22]提出了由動(dòng)脈順應(yīng)性描述的動(dòng)脈壁模型，見(jiàn)式(2).

(2)

式中：A0是跨壁壓力為零時(shí)動(dòng)脈的橫截面積；Am是動(dòng)脈的最大橫截面積，參數(shù)Ca和Cb用來(lái)描述動(dòng)脈力學(xué)性質(zhì).在后面的模擬中，我們使用Ca=0.09 mmHg-1,Cb=0.03 mmHg-1，該參數(shù)數(shù)值來(lái)自于論文[20-21].

1.3 組織模型

組織通常被視為彈性體，組織的變形近似為線彈性.因此，組織的變形可以用一組張量形式的彈性方程來(lái)描述，即平衡方程式(3)、幾何方程式(4)和物理方程式(5).

平衡方程：

σij,j+Fi=0,

(3)

式中：σij=σji，是應(yīng)力張量；Fi是體積力.

幾何方程：

(4)

式中：εij=εji，是應(yīng)變張量.

物理方程：

σij=Cijklεkl,

(5)

式中，Cijkl是剛度張量.

在這個(gè)問(wèn)題中，組織的幾何形狀、外力分布和約束條件都是對(duì)稱的.因此，該問(wèn)題可視為空間軸對(duì)稱問(wèn)題.一般來(lái)說(shuō)，極坐標(biāo)r,θ,z比笛卡爾坐標(biāo)x,y,z能更方便地解決這個(gè)問(wèn)題.由于對(duì)稱性，任意一點(diǎn)的應(yīng)力、應(yīng)變和位移都是獨(dú)立于θ的.也就是說(shuō)，它們是r和z的函數(shù).更重要的是，為了保持對(duì)稱，我們可以得到τrθ=τzθ=uθ=0.因此，對(duì)于空間軸對(duì)稱問(wèn)題，可以用柱坐標(biāo)形式具體地寫(xiě)出一組彈性方程.

平衡方程：

(6)

幾何方程：

(7)

物理方程：

(8)

式中：θ=εr+εθ+εz；E是彈性模量；μ是泊松比.

2 迭代算法

由于邊界條件的非線性特性，直接用有限元法求解彈性方程是困難的.為了克服這一問(wèn)題，提出了利用虛功原理的迭代法.

由于先前的假設(shè)忽略了物體的體積力，因此彈性體只受表面力的作用.對(duì)于具有體積τ和表面積S的彈性體，S分為力面域Sσ和位移面域Su.表面力，其分量是Xi，作用于Sσ.位移，其分量為ui，作用在Su.假定物體中存在任意的虛位移，由于Su上位移是固定的，物體表面力在Su上所做的功為零.不失一般性，可以將表面力作用的區(qū)域設(shè)置為S.因此，表面力所做功的變分是

(9)

由于變形，物體的彈性勢(shì)能變分為

(10)

應(yīng)力張量具有對(duì)稱性，將式(7)代入式(10)，可將物體彈性勢(shì)能的變分改寫(xiě)為

(11)

總能量變分為

(12)

由于平衡方程σij,j=0，有

(13)

利用高斯散度定理，可以得到

(14)

其中，lj是表面法方向的方向余弦.

因此，式(12)可以改寫(xiě)為

(15)

當(dāng)彈性體處于平衡狀態(tài)時(shí)，總能量最小.為了使總能量達(dá)到最小，我們選擇δui=-c(Xi+σijlj),c>0來(lái)更新ui使得δE<0.

3 數(shù)值實(shí)驗(yàn)

在進(jìn)行數(shù)值實(shí)驗(yàn)之前，我們預(yù)先設(shè)定了一些參數(shù)的值：E=5×104Pa,μ=0.45,Ca=0.09 mmHg-1,Cb=0.03 mmHg-1,r0=0.02 cm,rm=0.252 31 cm, 收縮壓Ps=130 mmHg,舒張壓Pd=70 mmHg, 袖帶壓強(qiáng)Pc=100 mmHg ，袖帶外壁和內(nèi)壁的半徑Rm=2.7 cm，R0=2.2 cm,組織和袖帶長(zhǎng)度L=20 cm，Lc=4 cm.

由于幾何結(jié)構(gòu)的對(duì)稱性(參見(jiàn)圖1(a))，我們只分析了一半組織的縱切面.當(dāng)使用有限元法時(shí)，我們使用正方形單元.

3.1 量綱分析

(16)

(17)

通過(guò)忽略ΔVcΔPc，式(17)可以被進(jìn)一步簡(jiǎn)化為

(18)

在實(shí)際測(cè)量中，當(dāng)袖帶壓力為100 mmHg時(shí)，袖帶壓力振蕩幅度約為1.5 mmHg.在這種情況下，我們可以估計(jì)相對(duì)體積變化率為

(19)

ΔuR0≈6.68×10-4cm.

(20)

3.1.1 評(píng)價(jià)由式(20)可知，外部組織位移變化的數(shù)量級(jí)約為O(10-4)cm.因此，在數(shù)值計(jì)算中，網(wǎng)格細(xì)化所產(chǎn)生的數(shù)值誤差應(yīng)小于它.這一發(fā)現(xiàn)之所以如此重要，是因?yàn)樵跓o(wú)法達(dá)到精度的情況下，數(shù)值誤差會(huì)隨著迭代次數(shù)的增加而擴(kuò)大，這在以前的文獻(xiàn)中被忽略了.

3.2 單元分析

實(shí)際上，在使用有限元法時(shí)，需要選擇合適的單元和單元長(zhǎng)度，以保證計(jì)算精度，減少計(jì)算時(shí)間和內(nèi)存開(kāi)銷(xiāo).在這里，我們比較兩個(gè)元素，低階單元和高階單元.它們之間的主要區(qū)別在于高階單元具有二次位移行為，因?yàn)閮蓚€(gè)節(jié)點(diǎn)之間有一個(gè)中間節(jié)點(diǎn).

3.2.1 收斂階在數(shù)值計(jì)算中，我們一直關(guān)注收斂階.利用收斂階可以計(jì)算出合適的單元長(zhǎng)度.假設(shè)u*是彈性方程的真解，uh是當(dāng)單元長(zhǎng)度為h時(shí)的數(shù)值解.如果存在常數(shù)C使得|uh-u*|?Chα成立，則稱該數(shù)值方法的收斂階為α.于是有

(21)

對(duì)式(21)兩邊取對(duì)數(shù)，有

(22)

圖2 (a)、(b) 舒張壓和收縮壓下網(wǎng)格細(xì)化后的外部組織位移變化；(c) 收斂階Fig.2 (a) (b) Displacement change of external tissue by mesh refinement under diastolic and systolic pressure；(c) Convergence order

事實(shí)上，在我們的研究中，我們更關(guān)注血壓變化時(shí)的袖帶壓力振蕩或袖帶容積振蕩.在這里，我們研究外部組織從舒張壓到收縮壓的位移振蕩，即外部組織的位移振幅.

圖3 (a) 不同單元長(zhǎng)度下的外部組織位移振幅;(b) 網(wǎng)格細(xì)化后外部組織的位移幅度；(c) 收斂階Fig.3 (a) Displacement amplitude of external tissue under different element lengths; (b) Displacement amplitude change of external tissue by mesh refinement；(c) Convergence order

3.2.3 高階單元現(xiàn)在用高階單元來(lái)研究特定血壓下網(wǎng)格細(xì)化產(chǎn)生的位移變化.從圖4(a)、(b)可以看出，當(dāng)單元的長(zhǎng)度變小時(shí)，由網(wǎng)格細(xì)化引起的從袖帶中心到組織末端區(qū)域內(nèi)的外部組織的位移變化的數(shù)量級(jí)迅速減小到O(10-4) cm，甚至更小.

本文對(duì)參加“2018年江蘇省定向錦標(biāo)賽暨江蘇定向邀請(qǐng)賽”中南京普通高校的部分學(xué)生進(jìn)行隨機(jī)問(wèn)卷調(diào)查，共發(fā)放150份調(diào)查問(wèn)卷，回收150份調(diào)查問(wèn)卷，刪除無(wú)效問(wèn)卷25份后，有效問(wèn)卷為125份。

圖4 (a)、(b)舒張壓和收縮壓下網(wǎng)格細(xì)化后的外部組織位移變化；(c)收斂階Fig.4 (a)、(b) Displacement change of external tissue by mesh refinement under diastolic and systolic pressure；(c) Convergence order

圖5 (a) 不同單元長(zhǎng)度下的外部組織位移振幅；(b) 網(wǎng)格細(xì)化后外部組織的位移幅度；(c) 收斂階Fig.5 (a) Displacement amplitude of external tissue under different element lengths； (b) Displacement amplitude change of external tissue by mesh refinement；(c) Convergence order

3.2.4 評(píng)價(jià)無(wú)論是特定血壓下的位移問(wèn)題還是一段時(shí)間內(nèi)的位移變化問(wèn)題，高階單元所需的單元長(zhǎng)度是低階單元的8倍以上，高階單元的收斂階是低階單元的2倍左右.這些都是高階單元的優(yōu)點(diǎn)，我們將在后面的數(shù)值模擬中使用它們.

從表1可以看出，網(wǎng)格生成對(duì)于特定血壓下的位移問(wèn)題非常重要.需要大量的網(wǎng)格才能達(dá)到預(yù)期的精度.事實(shí)上，這一問(wèn)題在以往的文獻(xiàn)中并沒(méi)有得到重視.然而，對(duì)于位移變化問(wèn)題，有一個(gè)有趣的結(jié)果是我們只使用較大長(zhǎng)度的單元就可以滿足精度要求.實(shí)際上，我們研究的重點(diǎn)是袖帶壓力振蕩.通過(guò)分析可知，在數(shù)值模擬中采用較大長(zhǎng)度的單元能獲得較高的精度和較少的計(jì)算成本.

表1 單元種類和問(wèn)題種類的影響

4 理論驗(yàn)證

本節(jié)給出一些假設(shè)來(lái)簡(jiǎn)化模型，以相關(guān)的理論分析來(lái)驗(yàn)證迭代算法的準(zhǔn)確性.由于在袖帶下組織外側(cè)的邊界條件幾乎保持不變，在這里對(duì)模型進(jìn)行了一些簡(jiǎn)化：

1) 組織是軸對(duì)稱的，無(wú)限長(zhǎng).

2) 平移不變性，即徑向應(yīng)力和位移不隨縱向坐標(biāo)的變化而變化.

4.1 Love函數(shù)

在1.3節(jié)中，我們給出了空間軸對(duì)稱問(wèn)題的平衡方程、物理方程和幾何方程.因此，我們可以得到位移控制方程，如式(23)所示

(23)

(24)

?4φ=0.

(25)

式(25)說(shuō)明φ(r,z)是雙調(diào)和函數(shù).

此外，通過(guò)將位移分量代入物理方程和幾何方程，可以得到應(yīng)力分量

(26)

因此，一旦知道了雙調(diào)和函數(shù)是什么，就可以得到所有的位移和應(yīng)力分量.

由徑向位移的變量分離和平移不變性，可以假設(shè)φ=f0(r)+f1(r)z+f(z).由于?4φ=0，可以得到

(27)

由于r和z的任意性,可以得到

(28)

因此，φ(r,z)可以寫(xiě)成φ=f0(r)+f1(r)z+C1z2+C2z3.于是，可以得到縱向位移

(29)

由于w|z=0=0，有2(1-μ)Lf0(r)+2C1(1-2μ)=0.因此式(29)可以改寫(xiě)為

(30)

同時(shí)，可以得到

(31)

由于ur,σr,σz,w獨(dú)立于f0,C1，可以忽略它們的影響，令f0=0,C1=0，于是φ(r,z)可以重寫(xiě)為φ=f1(r)z+C2z3.

Lf1(r)=C3,

其中C3是常數(shù).如果Lf1(r)=C3成立，則L2f1(r)=0顯然也成立.

另一方面，知道應(yīng)該有兩個(gè)邊界條件和兩個(gè)未知參數(shù)C2,C3.然而，可發(fā)現(xiàn)，當(dāng)f1有一個(gè)常數(shù)之間的差異，位移和應(yīng)力分量不發(fā)生改變，因此，只需要三個(gè)條件.

由Lf1(r)=C3，可知

(32)

由于f1顯示表達(dá)，則有

(33)

4.2 位移和應(yīng)力表達(dá)式

在該問(wèn)題下給出三個(gè)條件，在內(nèi)部組織上給出位移邊界條件，在外部組織上設(shè)置應(yīng)力邊界條件.此外，橫截面上的縱向應(yīng)力積分等于0，即

因此，我們有

(34)

4.3 理論和數(shù)值結(jié)果的比較

如前所述，只考慮位移在一段時(shí)間內(nèi)的變化.因此，有必要對(duì)位移變化問(wèn)題的理論結(jié)果和數(shù)值結(jié)果進(jìn)行比較.

根據(jù)先前的估計(jì)，當(dāng)袖帶壓力振蕩的振幅為1.5 mmHg時(shí)，外部組織的位移變化為O(10-4) cm(見(jiàn)式(20)).此外，根據(jù)式(34)，可以得到外部組織的位移

進(jìn)而，有

在量綱分析中，可知道ΔuR0≈6.68×10-4cm，ΔPc≈1.5 mmHg，這樣就可以得到內(nèi)部組織的位移變化Δu≈0.043 cm，這與數(shù)值解(≈0.04 cm)幾乎是一致的.這在一定程度上說(shuō)明了我們理論的正確性，這一結(jié)果也說(shuō)明動(dòng)脈壁變形達(dá)到O(10-2) cm或超過(guò)一定程度時(shí)，會(huì)出現(xiàn)明顯的袖帶壓力振蕩幅度.

5 未來(lái)展望

今后，我們將進(jìn)一步對(duì)非軸對(duì)稱問(wèn)題進(jìn)行理論分析和數(shù)值模擬.一方面，動(dòng)脈和組織的幾何形狀不是軸對(duì)稱的.另一方面，動(dòng)脈在受到外力壓迫時(shí)會(huì)發(fā)生塌陷，很難保持圓形.經(jīng)驗(yàn)動(dòng)脈壁模型就不再適合.因此，在非軸對(duì)稱問(wèn)題中有許多有趣的問(wèn)題值得探討.