流固耦合作用下人體上呼吸道內(nèi)氣流運(yùn)動(dòng)特性數(shù)值仿真研究

孫 棟 李福生 徐新喜 趙秀國 譚樹林

（軍事醫(yī)學(xué)科學(xué)院衛(wèi)生裝備研究所 國家生物防護(hù)裝備工程技術(shù)研究中心，天津 300161）

引言

上呼吸道是呼吸系統(tǒng)的重要組成部分，是人體與外界環(huán)境進(jìn)行氣體交換的主要通道。隨著科學(xué)技術(shù)的迅猛發(fā)展以及大氣環(huán)境的不斷惡化，這條與外界物質(zhì)交換最為頻繁的通道越來越多地引起人們的關(guān)注。近年來，國內(nèi)外許多學(xué)者從生物流體力學(xué)的角度對(duì)其進(jìn)行了大量研究。徐新喜等對(duì)循環(huán)呼吸模式下人體口喉模型內(nèi)氣流組織特性進(jìn)行了數(shù)值仿真和實(shí)驗(yàn)研究［1］；Wang等對(duì)上呼吸道內(nèi)的呼吸氣流組織形式進(jìn)行了數(shù)值分析［2］，Nithiarasu等利用無結(jié)構(gòu)網(wǎng)格對(duì)上呼吸道中的層流和湍流計(jì)算進(jìn)行了數(shù)值模擬［3］，趙秀國等對(duì)包括前3級(jí)支氣管在內(nèi)的人體上呼吸道模型內(nèi)的氣流運(yùn)動(dòng)特性及氣溶膠沉積進(jìn)行了PIV實(shí)驗(yàn)研究［4－5］。目前的研究大多假定呼吸道壁面為剛性壁面，而實(shí)際上呼吸道內(nèi)的真實(shí)呼吸流是典型的流體-結(jié)構(gòu)交互作用下的流動(dòng)，這一流固耦合過程對(duì)呼吸道內(nèi)的氣流組織形式、氣溶膠擴(kuò)散沉積以及壁面剪切應(yīng)力的分布有著重要的影響。

流固耦合力學(xué)模型的發(fā)展和成熟，為研究人體上呼吸道內(nèi)的氣流運(yùn)動(dòng)提供了一種嶄新的思路。利用流固耦合力學(xué)模型，能較好地分析氣流和呼吸道壁面的相互耦合作用，避免以往單純假設(shè)呼吸道壁面為剛性壁面，從而使理論數(shù)值模擬分析與實(shí)際呼吸道內(nèi)的流動(dòng)情況更為接近，拓展了人們對(duì)呼吸道內(nèi)氣流運(yùn)動(dòng)的研究思路和研究方法。本研究通過構(gòu)建流固耦合力學(xué)模型，對(duì)流固耦合作用下包括前三級(jí)支氣管在內(nèi)的上呼吸道內(nèi)穩(wěn)態(tài)氣流運(yùn)動(dòng)特性進(jìn)行了數(shù)值模擬，系統(tǒng)分析了流固耦合作用下上呼吸道壁面的形變特點(diǎn)、壁面剪切應(yīng)力分布以及呼吸道內(nèi)的氣流運(yùn)動(dòng)特點(diǎn)。

1 人體上呼吸道模型和數(shù)值仿真方法

1.1 人體上呼吸道模型

研究所采用的上呼吸道模型主要由兩部分組成：口喉模型，包括口腔、咽、喉和氣管；前3級(jí)支氣管模型，包括G0～G3級(jí)支氣管模型。口喉模型的尺寸參考了ARLA（Aerosol Research Laboratory of Alberta）模型［6］和 Stapleton K W 模型［7］的尺寸，前 3級(jí)支氣管模型采用了Weibel的支氣管模型［8］，如圖1所示。

1.2 控制方程

圖1 人體上呼吸道模型Fig.1 Human upper respiratory tract model

1.2.1 流體控制方程上呼吸道中流體湍流模型采用耦合了Langtry-Menter轉(zhuǎn)捩模型的 Menter's SST兩方程模型。Langtry-Menter轉(zhuǎn)捩模型是Menter與Langtry提出的嚴(yán)格建立在局部變量基礎(chǔ)上的基于經(jīng)驗(yàn)關(guān)聯(lián)式的轉(zhuǎn)捩模型，因而能夠與CFD代碼耦合。該模型基于兩個(gè)輸運(yùn)方程，即間歇因子的輸運(yùn)方程和轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)輸運(yùn)方程。前者與SST湍流模型耦合，用以開啟轉(zhuǎn)捩點(diǎn)下游的湍動(dòng)能產(chǎn)生項(xiàng)，從而觸發(fā)當(dāng)?shù)氐霓D(zhuǎn)捩；后者用來捕捉非當(dāng)?shù)赝牧鲝?qiáng)度的影響，且能夠?qū)⒔?jīng)驗(yàn)關(guān)聯(lián)式與間歇因子輸運(yùn)方程中的轉(zhuǎn)捩起始準(zhǔn)則聯(lián)系起來，得到張量表達(dá)式［9－10］。

間歇因子γ的輸運(yùn)方程為

式中，生成項(xiàng)Pγ和耗散項(xiàng)Eγ的表達(dá)式分別為

當(dāng)?shù)剞D(zhuǎn)捩動(dòng)量厚度雷諾數(shù)～Reθt輸運(yùn)方程為

式中，源項(xiàng) Pθt定義為

式中，U為當(dāng)?shù)厮俣龋旌虾瘮?shù)Fθt的表達(dá)式為式中，θ?BL為邊界層動(dòng)量損失厚度，δ為邊界層厚度，混合函數(shù)Fθt用來關(guān)掉邊界層內(nèi)的生成項(xiàng)Pθt從而使得從邊界層以外擴(kuò)散得到，F(xiàn)wake的作用是保證混合函數(shù)Fθt在尾跡區(qū)不被激活。

1.2.2 壁面（固體）結(jié)構(gòu)控制方程

將上呼吸道壁材料假設(shè)為線彈性材料，而非傳統(tǒng)的剛性壁面假設(shè)，這樣考慮可以使其與實(shí)際氣流情況更為接近。根據(jù)文獻(xiàn)［11］，設(shè)定氣管為均勻各項(xiàng)同性的材料，即服從胡克定律的薄壁彈性模型，其楊氏彈性模量為9 MPa，泊松比為0.4。

固體模型控制方程可表達(dá)為式中，σs為是氣管壁應(yīng)力張量，ρs為氣管壁密度，as為氣管壁加速度。

在流固耦合交界面上，應(yīng)當(dāng)滿足

賞識(shí)教育作為幼兒教育中的首要教學(xué)方法，是在幼兒教育過程中通過鼓勵(lì)性的話語和肢體語言啟發(fā)幼兒，讓幼兒在受教育中得到更多的愛與關(guān)心。不僅能夠幫助幼兒建立自信心，而且有助于幼兒心智全面健康發(fā)展。以下分三個(gè)方面對(duì)賞識(shí)教育在幼兒教育中的應(yīng)用實(shí)踐進(jìn)行探索。

式中，d為位移，n為邊界法向，下標(biāo)s和 f分別表示固體和流體。

1.3 數(shù)值仿真方法

流固耦合問題一般分為兩類，一類為流-固單向耦合，一類為流-固雙向耦合。本研究對(duì)流固耦合采取雙向耦合的方式，其基本思路是：分別計(jì)算流場(chǎng)和固體結(jié)構(gòu)，然后通過中間平臺(tái)來交換耦合量。在每次大迭代中，分別進(jìn)行1次流體和固體計(jì)算，并交換兩次數(shù)據(jù)（每個(gè)方向各1次），直到最終收斂。其中，在固體結(jié)構(gòu)計(jì)算上，由于人體上呼吸道的生理結(jié)構(gòu)極其復(fù)雜，咽喉-氣管部位又受到周圍組織的非線性、非固定的力學(xué)約束，仿真非常困難。為簡(jiǎn)化起見，模型通過合理設(shè)置壁面的彈性系數(shù)以及在模型的口部添加固定位移約束，以達(dá)到仿真的目的。而在流場(chǎng)的計(jì)算上，設(shè)置其邊界條件為：在低強(qiáng)度呼吸條件（呼吸流量30 L/imn）下，口喉模型的入口邊界采用速度入口條件，入口邊界假定速度均勻分布，支氣管出口采用壓力出口邊界條件，出口相對(duì)壓力為0。

1.4 數(shù)值仿真模型驗(yàn)證

應(yīng)用粒子圖像速度儀（particleimage velocimetry，PIV），對(duì)人體在呼吸流量為30 L/min的低強(qiáng)度呼吸條件下上呼吸道內(nèi)的穩(wěn)態(tài)氣流運(yùn)動(dòng)特性進(jìn)行了試驗(yàn)研究；將試驗(yàn)結(jié)果與數(shù)值仿真結(jié)果進(jìn)行對(duì)比，驗(yàn)證數(shù)值仿真方法的準(zhǔn)確性。

圖2 人體上呼吸道喉部氣流組織形式試驗(yàn)結(jié)果與數(shù)值仿真結(jié)果對(duì)比。（a）試驗(yàn)結(jié)果；（b）仿真結(jié)果Fig.2 Comparison of airflow patterns between experimental results and large eddy simulation results of the larynx.（a）Experimental results；（b）Large eddy simulation results

圖2為上呼吸道喉部的試驗(yàn)結(jié)果與數(shù)值仿真結(jié)果的對(duì)比。圖2（a）為PIV試驗(yàn)各點(diǎn)氣流速度的測(cè)量結(jié)果，圖2（b）是仿真氣流的速度矢量?？梢钥闯觯瑑烧呔窃跉饬鬟M(jìn)入咽部時(shí)，截面積變小而氣流速度增大，整體氣流組織形式較為一致。同時(shí)試驗(yàn)速度最大為8.51 m/s，仿真速度最大為7.85 m/s，且測(cè)量數(shù)值與仿真數(shù)值誤差最大不超過10％，吻合較好，從而說明數(shù)值仿真方法的準(zhǔn)確性。

2 結(jié)果與討論

2.1 上呼吸道形變分析

圖3為上呼吸道壁面的位移變化云圖。圖3（a）為上呼吸道壁面的總體形變對(duì)比，顯示了人體上呼吸道壁面的總體位移特點(diǎn)。其中，總體位移是指上呼吸道壁面上各點(diǎn)受到氣流作用后相對(duì)初始時(shí)刻位置的變化，由3個(gè)方向的位移組成：前后位移dx、左右位移dy、縱向位移dz，即d=dx+dy+dz。圖3（b）顯示了上呼吸道壁面縱向形變特點(diǎn)，圖3（c）和圖3（d）顯示了上呼吸道壁面的前后形變特點(diǎn)。通過圖3（a）可以看出，在呼吸流量為30 L/min的穩(wěn)態(tài)氣流作用下，人體上呼吸道整體向后運(yùn)動(dòng)，咽喉部位向后的位移量較小，氣管支氣管向后的位移量較大，且氣管支氣管的位移量由上往下呈逐級(jí)遞增的狀態(tài)，氣管前壁面受到拉伸，后壁面受到壓縮。由圖3（c）和圖3（d）可以看出，咽喉部位受到氣流運(yùn)動(dòng)的作用，造成前壁面的中間前后位移小、兩邊前后位移大，后壁面的中間前后位移大、兩邊前后位移小，造成前壁凹、后壁凸的現(xiàn)象，在咽部的后壁面中心處形成一個(gè)小凸包。通過圖3（b）可以看出，在穩(wěn)態(tài)氣流的作用下，上呼吸道前壁面向下運(yùn)動(dòng)、后壁面向上運(yùn)動(dòng)，且咽喉部位后壁面向上運(yùn)動(dòng)的位移量最大，而氣管前壁面向下的位移量最大。

圖3 呼吸流量為30 L/min時(shí)上呼吸道壁面形變。（a）形變前后對(duì)比；（b）縱向位移；（c）和（d）前后位移Fig.3 The deformation of human upper respiratory tract model with breathing intensity Q=30 L/min.（a）Comparison of pre-and post-deformation；（b）Longitudinal displacement；（c）and（d）Lateral displacement

2.2 上呼吸道壁面剪應(yīng)力分布分析

圖4 呼吸流量為30 L/min時(shí)上呼吸道壁面的Von Mises stress分布。（a）口喉模型前壁正視圖；（b）口喉模型后壁正視圖；（c）氣管支氣管前壁正視圖；（d）氣管支氣管后壁正視圖Fig.4 The Von Mises stress distribution of human upper respiratory tract model with breathing intensity Q=30 L/min.（a）Front view of anterior wall in human mouth-throat model；（b）Front view of posterior wall in human mouth-throat model；（c）Front view of anterior wall in trachea to triple bifurcation region；（d）Front view of posterior wall in trachea to triple bifurcation region

圖4為人體上呼吸道壁面剪應(yīng)力分布云圖。其中，（a）和（b）為口喉模型上剪應(yīng)力的分布狀況，（c）和（d）為氣管支氣管上剪應(yīng)力的分布狀態(tài)?？梢钥闯觯诤砟Ｐ椭惺艿降谋诿婕羟袘?yīng)力較大，氣管支氣管受到的壁面剪切應(yīng)力較小，剪切應(yīng)力最大處可達(dá)30.34 Pa，且氣管支氣管所受壁面剪應(yīng)力由上往下呈現(xiàn)出逐級(jí)遞減的狀態(tài)。此外，口喉模型前壁面所受剪切應(yīng)力小于后壁面所剪應(yīng)力，氣管前壁面所受剪切應(yīng)力大于氣管后壁面所受剪切應(yīng)力。通過圖4中的（a）和（b）可以看出，氣流進(jìn)入咽部時(shí)，在咽部的后壁面中心處以及咽部兩側(cè)形成剪應(yīng)力集中區(qū)，剪應(yīng)力較大，而咽部的前壁面所受剪應(yīng)力較小。由于喉部結(jié)構(gòu)形狀復(fù)雜，后壁面中心處及兩側(cè)產(chǎn)生應(yīng)力集中區(qū)；同時(shí)，由于聲門處的截面減小和氣流轉(zhuǎn)向，在聲門上方的后壁面上產(chǎn)生兩個(gè)接近對(duì)稱的剪應(yīng)力集中區(qū)。通過圖4中的（c）和（d）可以看出，氣流到達(dá)聲門后，在聲門兩側(cè)產(chǎn)生較大的剪應(yīng)力，在中間部位產(chǎn)生較小的剪應(yīng)力，由于氣流在聲門處向氣管的前壁面產(chǎn)生噴射，造成氣管前壁面上的剪應(yīng)力明顯高于后壁面上的剪切應(yīng)力。在氣管、支氣管的分叉處出現(xiàn)剪切應(yīng)力集中區(qū)，該處的剪切應(yīng)力明顯大于其他部位。

2.3 上呼吸道內(nèi)流場(chǎng)分析

圖5為呼吸流量30 L/min時(shí)人體上呼吸道模型內(nèi)不同截面的軸向速度分布。其中，圖5（a）為口喉模型內(nèi)對(duì)稱垂直截面內(nèi)的氣流流場(chǎng)分布，圖5（b）為氣管支氣管對(duì)稱截面內(nèi)的氣流速度分布?？梢钥闯觯诳诤砟Ｐ蛢?nèi)，氣流分別在咽部和喉部形成兩個(gè)速度增長點(diǎn)，并且在咽部產(chǎn)生分離現(xiàn)象，在咽部外形成分離區(qū)，見圖5（c）中C-C＇截面。氣流在聲門位置受到幾何結(jié)構(gòu)的限制，產(chǎn)生湍流噴射的現(xiàn)象。聲門下方的氣管內(nèi)發(fā)生流動(dòng)分離，靠近氣管內(nèi)壁氣流速度較高，見圖5（c）中 D-D＇和E-E＇截面。隨著與聲門距離的增加，氣管內(nèi)、外壁的氣流速度差逐漸減少，但是氣管內(nèi)壁的氣流速度還是要高于外壁的氣流速度，如圖5（c）中F-F＇截面所示。氣流在支氣管分叉處發(fā)生分離，靠近支氣管內(nèi)壁速度較高，靠近外壁速度較低，見圖5（c）中 H-H＇截面。在 G2和 G3級(jí)支氣管內(nèi)，產(chǎn)生了非對(duì)稱流量的氣流流動(dòng)現(xiàn)象。

圖5 呼吸流量為30 L/min時(shí)上呼吸道不同截面內(nèi)的速度分布。（a）口喉垂直截面；（b）氣管支氣管垂直截面；（c）橫截面Fig.5 Axial velocity contours at different cross section in human upper respiratory tract model with breathing intensity Q=30 L/min.（a）The vertical section of human mouth-throat model；（b）The vertical section of trachea to triple bifurcation region；（c）The cross sections

3 結(jié)論

1）在呼吸流量為30 L/min的穩(wěn)態(tài)氣流作用下，人體的上呼吸道整體向后運(yùn)動(dòng)，咽喉部位向后的位移量較小，氣管、支氣管向后的位移量較大，且氣管、支氣管的位移量由上往下呈逐級(jí)遞增的狀態(tài)，氣管前壁面受到拉伸，后壁面受到壓縮。

2）口喉模型中受到的壁切面剪應(yīng)力較大，氣管、支氣管受到的壁面剪應(yīng)力較小，且氣管、支氣管所受壁面剪切應(yīng)力由上往下呈現(xiàn)出逐級(jí)遞減的狀態(tài)；口喉模型前壁面所受剪切應(yīng)力小于后壁面所受剪切應(yīng)力，氣管前壁面所受剪切應(yīng)力大于氣管后壁面所受剪切應(yīng)力。

3）在口喉模型內(nèi)，氣流分別在咽部和喉部形成兩個(gè)速度增長點(diǎn)，并且在咽部發(fā)生分離現(xiàn)象；氣流在聲門位置產(chǎn)生湍流噴射的現(xiàn)象，聲門下方的氣管內(nèi)發(fā)生流動(dòng)分離，靠近氣管內(nèi)壁氣流速度較高；氣流在支氣管分叉處發(fā)生分離，靠近支氣管內(nèi)壁速度較高。

［1］徐新喜，趙秀國，譚樹林，等.循環(huán)呼吸模式口喉模型內(nèi)氣流運(yùn)動(dòng)特性數(shù)值模擬［J］.力學(xué)學(xué)報(bào)，2010，42（2）：182 －190.

［2］Wang Ying， LiuYingxi， SunXiuzhen， etal. Numerical analysis of respiratory flow patterns within human upper airway［J］.Acta Mecheanica Sinica，2009，25（6）：737 －746.

［3］Nithiarasu P，Liu CB，Massarotti N.Laminar and turbulent flow calculations through a modelhuman upperairways using unstructured meshes［J］.Communications in Numerical Methods Engineering，2007，23（12）：1057－1069.

［4］趙秀國，徐新喜，譚樹林，等.人體上呼吸道內(nèi)穩(wěn)態(tài)氣流運(yùn)動(dòng)特性的PIV初步試驗(yàn)研究［J］.實(shí)驗(yàn)流體力學(xué)，2009，23（4）：60－64.

［5］徐新喜，趙秀國，譚樹林，等.人體上呼吸道內(nèi)氣流運(yùn)動(dòng)特性的數(shù)值模擬分析［J］.計(jì)算力學(xué)學(xué)報(bào)，2010，27（5）：881－886.

［6］Grgic B，F(xiàn)inlay WH，Heenan AF.Regional aerosol deposition and flow measurements in an idealized mouth and throat［J］.Aerosol Science，2004，35：21 －32.

［7］Stapleton KW，Guentsch E，Hoskinson MK，et al.On the suitability of k-εturbulence modeling for aerosol deposition in the mouth and throat：a comparison with experiment［J］.Aerosol Science，2000，31（6）：739 －749.

［8］Weibel ER.Morphometry of the human lung［M］.New York：Academic Press，1963：1 － 3.

［9］Menter FR，Langtry RB，Likki SR，et al.A correlation based transition model using local variables Part 1-Model Formulation［R］.ASME 2004-GT-53452.

［10］Langtry RB，Menter FR，Likki SR，et al.A correlation based transition model using local variables Part 2-Test Cases and Industrial Applications［R］.ASME 2004-GT-53434.

［11］Wlofgang AW，Timon R.Fluid-structure interaction in lower airways of CT-based lung geometries［J］.International Journal for Numerical Methods in Fluids，2008，57（5）：653 －675.