?

主動脈穿刺型軸流血泵折邊結構葉輪的數值模擬及溶血分析

2024-02-02 08:22荊騰潘愛娣顧發東王秀禮
排灌機械工程學報 2024年2期
關鍵詞:血泵葉頂葉輪

荊騰,潘愛娣,顧發東,王秀禮

(江蘇大學國家水泵及系統工程技術研究中心,江蘇 鎮江 212013)

心力衰竭是指人體心臟機能減弱,泵血能力不足,無法產生足夠的血壓和流量導致人體血液循環出現病理情況[1].目前治療心力衰竭最有效的方法是心臟移植,但由于供體缺乏導致每年有大量的心衰患者在等待移植的過程中死亡.人工心臟泵能夠對衰竭心臟的功能進行補償,為血液循環提供動力,從而可為等待心臟移植的患者贏得一定時間,也為不適合進行心臟移植的患者帶來了生存的希望[2].

軸流血泵因其結構簡單、體積小、效率高等優點,成為血泵研究發展的主要方向之一[3].但傳統軸流血泵的轉速較高,且存在葉頂間隙,會產生葉頂間隙泄漏流[4],從而在葉頂區產生較高的切應力,影響血泵血液相容性,阻礙其臨床應用和發展.研究發現,影響旋轉式血泵血液相容性問題的關鍵是葉輪結構[2],其設計不當會直接導致血泵中溶血和血栓的形成[5].柳光茂等[6]設計的紡錘形轉子葉輪帶有分流葉片和懸臂葉片的尾導結構,利用數值模擬及粒子成像測速(particle image velocity,PIV)方法分析了血泵的水力性能及溶血特性,結果表明其流道內血流無較大分流,溶血性能良好,壓力流量性能滿足臨床需要.KANNOJIYA等[7]對三彎葉片葉輪和螺旋葉片葉輪的軸流血泵進行了三維模擬研究,發現螺旋葉片葉輪能夠改善血液流場并大大降低溶血.許多學者研究了葉輪與泵殼間隙大小對溶血和血小板活化的影響,間隙大小被認為是影響血液損傷的重要參數之一[8-9].WIEGMANN等[8]研究發現潛在的破壞性切應力與間隙尺寸和葉輪葉片數有關.停滯區和再循環區的范圍隨著葉片數量的減少和半開放式葉輪的出現而減少,但隨著間隙的縮小而增加,且小間隙會對紅細胞產生高切應力并降低沖刷性能.GRAEFE等[9]研究表明,間隙長度通過影響紅細胞的停留時間,在數量上顯著影響血液創傷的風險.大間隙會導致顯著的流動擾動,從而增加渦量和暴露時間.同時,研究還發現,保持相同的周向速度和間隙,小直徑泵有較好的血液相容性.

文中以主動脈穿刺型軸流血泵為研究對象,設計一種折邊葉片結構葉輪[10],使得血液進入葉輪后,在封閉的葉輪腔體內流動,直到從出口流出,避免因葉頂間隙泄漏造成血液破壞從而引發溶血和血栓的產生,并通過CFD數值模擬技術,對折邊不等間距葉輪、折邊等間距葉輪及非折邊不等間距葉輪的內部流場及溶血情況進行比較分析,驗證折邊葉輪在提高血液相容性方面的優勢,為軸流血泵的后續研究奠定基礎.

1 模型與方法

1.1 軸流血泵模型

以主動脈軸流血泵為研究對象,采用NX 11.0三維軟件建模,其三維模型如圖1所示,主要包括前導葉、葉輪、后導葉、轉軸、泵殼以及驅動裝置,該泵由主動脈弓位置微創植入,入口位于左心室,出口位于升主動脈,直接將血液從左心室泵入主動脈,減少了進出口導管,降低了血栓的產生.為了提高軸流血泵的血液相容性,文中提出了折邊葉片結構,設計了折邊等間距葉輪和折邊不等間距葉輪,并與非折邊不等間距葉輪進行模擬分析比較,驗證折邊葉輪葉片在提高血液相容性方面的有效性.3種不同結構葉輪模型如圖2所示,其中,非折邊不等間距葉輪采用傳統葉輪設計方法設計,葉片間距從入口到出口逐漸增大,葉頂間隙為0.2 mm;折邊不等間距葉輪和非折邊不等間距葉輪具有相同的葉片包角和安放角,只增加折邊結構以抑制葉頂間隙泄漏;在折邊不等間距葉輪基礎上設計得到折邊等間距葉輪,它與折邊不等間距葉輪擁有相同的折邊結構,但其葉片間距相等.折邊不等間距葉輪血泵主要設計參數中,葉輪直徑D2=20 mm,葉輪葉片數Z=4,葉頂弦長c=51 mm,葉輪葉片輪轂半徑hl=15.5 mm,葉高h=1.75 mm,葉片厚度b=0.5 mm,前導葉葉片數Zqd=6,后導葉葉片數Zhd=5,轉速n=8 000 r/min,額定流量Qopt=5 L/min,通過數值模擬得到血泵對血液的揚程約為2.14 m.為使3種血泵模擬后的揚程工況一致,非折邊不等間距葉輪血泵轉速n1=10 300 r/min,折邊等間距葉輪血泵轉速n2=11 450 r/min.

圖1 主動脈穿刺型軸流血泵的結構示意圖Fig.1 Schematic diagram of structure of aortic puncture type axial flow blood pump

圖2 葉輪的三維模型Fig.2 Three-dimensional models of impellers

1.2 計算域與網格劃分

圖3為主動脈穿刺型軸流血泵計算域示意圖,分別由進口延長段、前導葉計算域、葉輪計算域、后導葉計算域及出口段組成.

圖3 主動脈穿刺型軸流血泵計算域示意圖Fig.3 Schematic diagram of computational domain of aortic puncture type axial flow blood pump

由于葉輪結構復雜,采用結構化網格劃分比較困難,所以文中采用非結構多面體網格對各計算域進行劃分,其整體網格與局部計算域網格如圖4所示,各計算域之間采用Interface面進行連接.為了便于捕捉葉片表面區域的血液流動情況以及更精確求解壁面切應力的分布情況,對各計算域壁面網格進行邊界層網格加密處理,如圖5所示.

圖4 主動脈穿刺型軸流血泵的網格示意圖Fig.4 Schematic diagram of meshes of aortic puncture type axial flow blood pump

圖5 血泵壁面邊界層網格加密示意圖Fig.5 Schematic diagram of boundary layer meshes encryption on static wall surface of blood pump

1.3 物性參數及邊界條件設定

采用CFD軟件對血泵進行數值仿真模擬.由于血液屬于非牛頓流體,數值仿真時計算過程相當復雜,因此假設血液為黏性不可壓縮牛頓流體[11].血液參數參考正常成年人指標,即血液密度設置為1 055 kg/m3,黏度設置為0.003 5 Pa·s.參照正常成年人心臟的生理狀況,軸流血泵入口設置為質量流量入口條件,質量流量為0.087 917 kg/s,對應體積流量為5 L/min.血泵出口設置為壓力出口條件,其值參照正常成年人左心室與動脈連接處的血壓,約為100 mmHg,即13 330 Pa.湍流模型選用k-ε標準雙方程模型,求解方法采用計算收斂速度更快的SIMPLEC算法,計算時所有殘差的收斂標準均為10-5.

1.4 溶血模型

人工血泵在運行過程中會對人體血液造成一定的破壞,導致溶血狀況的發生,嚴重時將會導致血液的氧供給缺乏,威脅人體健康,因此需要對所設計血泵開展溶血預測分析.TETSUYA等[12]采用拉格朗日粒子跟蹤法對溶血模型進行了驗證,結果表明,拉格朗日粒子跟蹤法能夠實現對溶血值的預測.文中使用GIERSIEPEN等[13]創建的冪函數模型來建立主動脈穿刺型軸流血泵溶血預測模型,并結合粒子追蹤法來獲取單個紅細胞的流動軌跡,提取其中的相關物理量數據,進行溶血性能的預測.

數學公式為

(1)

式中:溶血指數HI定義為血漿游離血紅蛋白濃度的增量ΔHb與全血血紅蛋白濃度Hb的比值;τ為作用于血液的有效切應力,Pa;t為紅細胞暴露時間,s.

切應力由2部分組成,即湍流切應力(雷諾應力)和分子切應力,其表達式為

(2)

式中:μ為血液黏度;μt為湍流黏度;ρ為液體密度;k為湍動能;σij為Kronecker數.

BLUDSZUWEIT[14]提出切應力的標量計算關系式,把三維湍流的切應力轉換為一個代表切應力大小的標量,公式為

(3)

式中:τii,τjj,τij為黏性切應力分量.

文中采用有效切應力來進行計算,即

(4)

式中:σ1,σ2,σ3表示正應力.

文中基于拉格朗日粒子跡線追蹤法對軸流血泵進行溶血值預測,其原理如下:假設單個紅細胞在進入血泵入口時均未受到任何破壞,即溶血值為0[15].當血液流經血泵時,紅細胞所受到的切應力隨時都在發生改變.把單個紅細胞視為一個質點,記錄每個質點從進入到離開血泵的運動軌跡及其在該軌跡上不同時刻受到的切應力,通過積分得到該質點所受的破壞程度大小即紅細胞流經血泵的溶血預測值.具體方法如下:假設第p個質點在t0時刻進入血泵,tm時刻離開血泵,記錄此質點在t0,t1,…,ti,…,tm時刻受到的切應力,將紅細胞在Δti內的時均切應力τi代入,則其在ti-1時刻至ti時刻受力時間內的破壞概率dp,i[12]可表示為

(5)

其中:Δti=ti-ti-1.

通過上式可以得到質點p在第i個單位時刻的累積血液破壞程度Dp,i,Dp,i-1是ti-1時刻前的累積血液破壞程度,則在ti-1時刻至ti時刻還有(1-Dp,i-1)可以被破壞,且破壞概率為dp,i,公式為

Dp,i=Dp,i-1+(1-Dp,i-1)dp,i.

(6)

則以紅細胞為主體的流體質點的最終累積破壞程度為Dp,m.當分析大量流體質點時,重復上述步驟即可得到待研究質點的最終破壞程度,再對所有跡線的結果求取平均值,得到整體的血液損傷程度,即血泵的溶血指數HI.

1.5 網格無關性驗證

為考察計算中使用網格單元數的合理性,選擇折邊不等間距葉輪血泵進行計算域網格無關性驗證[16-17].以軸流血泵的揚程參數隨網格數量的變化作為判斷標準,采用5套網格對血泵性能進行模擬,計算結果顯示,當網格數超過354萬時,隨著網格數量的增加,軸流血泵的外特性數值變化率逐漸減小,且變化率小于1%,綜合考慮計算資源和時間成本,在后續軸流血泵的相關研究中,折邊不等間距葉輪血泵、非折邊不等間距葉輪血泵與折邊等間距葉輪血泵的網格數分別為354.123 6萬、374.833 2萬與365.688 9萬.

2 結果與分析

2.1 水力性能

2.1.1不同結構血泵的流場

圖6為不同葉輪結構血泵整體和局部流線及血液切應力分布,其中流線顏色代表軸流血泵流道內血液所受切應力值的大小.由圖6可知,3種葉輪血泵的血液流動趨勢大致相同,只在葉輪入口和出口區域有較大差別.非折邊不等間距葉輪在與前、后導葉的動靜交界面流域中流線較為紊亂,存在回流和渦流現象,同時伴有局部極高切應力;折邊等間距葉輪與前導葉的動靜交界區域存在大量的渦流,在與前、后導葉的交界處的切應力均較高,且范圍較廣;折邊不等間距葉輪在與前、后導葉交界的流域中有極少數流線存在渦流現象,且在與后導葉交界的流域中存在少量高于周邊切應力的區域.

圖6 不同血泵整體及局部流線圖和血液切應力Fig.6 Overall and local flow diagrams and blood shear stress of different blood pumps

流場內回流和渦流的產生使得紅細胞在切應力中的暴露時間變長,增大了血細胞的破碎概率,會引發溶血現象[1].與此同時,流場中渦流的產生會使該位置流域內的速度梯度增大,造成較大的切應力以及較強的湍動效應[18],使血泵能量損耗增大,效率降低.

由于葉輪與導葉之間的動靜干涉作用會引起局部流域血液運動狀態發生急劇變化[19],使得該位置速度梯度增大,從而引起較大的切應力.對比3種血泵流線圖可知,折邊不等間距葉輪血泵在葉輪入口與出口處的流線紊亂現象較少,這大大降低了血液的停留時間,使得流態得以改善.同時折邊不等間距葉輪無葉頂間隙,可有效減小葉頂流域中的較高切應力,因此流場中的切應力普遍較低,內部流場較為平穩,內流損耗相對較小.從流場分布來看,折邊不等間距葉輪的流場優化效果較好.

2.1.2不同結構葉輪葉頂間隙流域內的壓力場、速度矢量場

葉頂間隙血液泄漏是由于葉輪高速旋轉對流體做功,將機械能轉化為動能與壓能,即血液在葉輪流道中是一個增壓過程,葉輪流道中下游壓力要大于上游壓力.在壓差的作用下,流體會從高壓區流向低壓區.在葉輪流道中表現為血液從葉片壓力側經過葉頂間隙流向吸力側.葉頂間隙泄漏流的產生容易在葉頂區靠近葉片吸力側處引起葉頂泄漏渦以及誘導渦,造成較大速度梯度,并產生較高的切應力,對紅細胞造成損傷.由圖6的速度矢量分布可知,非折邊不等間距葉輪葉頂間隙處存在血液泄漏現象,與此同時,在葉輪出口段壓力側存在明顯的渦流,這在一定程度上會影響血泵的外特性;折邊等間距葉輪和折邊不等間距葉輪可以有效解決葉頂間隙泄漏問題,血液只在葉輪入口段附近存在渦流和回流現象,其中折邊等間距葉輪入口段吸力側存在較大渦流;折邊不等間距葉輪除了有效解決葉頂間隙的泄漏問題,其入口段吸力側的血液回流也明顯減少.由圖7的壓力分布云圖可知,折邊等間距葉輪的壓力分布云圖與其他2個葉輪有較大差別.折邊等間距葉輪的流場壓力從入口到出口先升高后降低,同時,在葉輪入口段存在較高壓力梯度,易引起紅細胞膜因內外壓差過大而破碎;而另外2個葉輪的流場壓力從入口至出口處于持續增加的過程.

圖7 不同結構葉輪軸截面葉頂間隙流域的壓力分布云圖及速度矢量分布Fig.7 Cloud diagram of pressure distribution and velocity vector distribution of impeller top clearance basin on impeller shaft section with different structures

2.1.3不同葉輪結構血泵內的軸向速度、切應力分布

由圖6流線圖可知,3種葉輪血泵流場中均會出現相對較高的切應力,且流道內流動較為復雜,單從流線圖上只能做粗略的定性分析,為了更為精確地評估3種葉輪血泵的內流場性能,提取血液在3種葉輪血泵流道內的軸向速度、切應力進行定量分析.圖8為不同葉輪血泵前導葉、葉輪和后導葉各過流斷面的軸向速度分布曲線圖.其中,A為非折邊不等間距葉輪結構,B為折邊等間距葉輪結構,C為折邊不等間距葉輪結構;Ⅰ,Ⅱ,Ⅲ分別代表前導葉、葉輪和后導葉區域;λ為量綱一化血泵軸向長度的位置,計算過程如下:

圖8 泵內各過流斷面軸向速度分布Fig.8 Axial velocity distribution on each overflow cross-section in pumps

設ΔZ為軸向長度增量;LⅠ,LⅡ和LⅢ分別為前導葉、葉輪、后導葉的軸向長度,對于前導葉區域,計算公式為λⅠ=ΔZ/LⅠ,將前導葉軸向長度量綱一化;對于葉輪區域,計算公式為λⅡ=(ΔZ-LⅠ)/LⅡ+1;對于后導葉區域,計算公式為λⅢ=(ΔZ-LⅠ-LⅡ)/LⅢ+(LⅡ-LⅠ)/LⅡ+1.

由圖8整體對比來看,3種葉輪血泵內流場軸向速度vr分布趨勢大體一致,僅在葉輪區域有一定區別.其中血泵B在葉輪和前、后導葉的交界處產生了較大的軸向速度梯度,而在葉輪區域內基本保持不變.血泵A和C在葉輪前、后導葉交界處的軸向速度梯度基本一致且小于血泵B,血泵A的速度變化比血泵C略微平緩.血液在葉輪入口與出口處產生較大軸向速度梯度是由于葉輪區域為動計算域,葉片旋轉對流體做功,因此血液在進入和離開葉輪區域時,軸向速度會迅速增大或降低從而產生較大的速度梯度,由此引發較大的切應力,增加了血液破壞的可能性.血液在流入或流出前導葉、后導葉的流道時,受到葉片自身體積阻塞流道的影響,也會在導葉的前緣與后緣處出現軸向速度迅速增大或減小的情況.

圖9為不同葉輪血泵前導葉、葉輪和后導葉各過流斷面切應力分布.

由式(2)—(4)可知,切應力與速度梯度成線性關系,即與圖8所示的軸向速度斜率分布趨勢相對應.其中,血泵C在葉輪入口與出口流域的切應力峰值最小,血泵B在葉輪葉片流道內的切應力最小,血泵A在葉輪葉片流道內的切應力大于血泵C.對比圖8可知,血泵A的軸向速度梯度略小于血泵C,由此推斷其切應力也應該略小于血泵C,然而血泵A在葉頂間隙泄漏流的影響下,流道中整體切應力升高,這也體現了折邊不等間距葉輪結構的優越性.

2.2 溶血性能

2.2.1跡線無關性驗證

采用拉格朗日方法進行溶血估算需要提取足夠多的跡線才能使計算得到的溶血指數HI不會隨跡線數量的改變而發生變化,計算結果才能夠更真實地反映血泵流域的溶血情況[2].文中對折邊不等間距葉輪血泵開展跡線無關性驗證,結果表明,當跡線數N取到1 163時,溶血指數HI開始趨于平緩,溶血值的變化率在5%以內,因此文中在血泵溶血預測計算中跡線數目取值為1 163條,此時溶血指數HI為4.009 9×10-3%.

血液流經某一區域后所產生的溶血破壞程度受很多因素的影響,包括流場切應力的大小、曝光時間的長短、血液本身的性質(例如物種、年齡、紅細胞比容等)、剪切力加載歷史以及湍流等[20],其中,定量研究主要集中于建立溶血破壞量與剪切力大小和曝光時間之間的對應關系[2].對于其他影響因素的研究雖然有助于探究細胞破壞的過程和機理、完善建立準確溶血模型的理論基礎,但是目前尚未引入定量溶血模型中.所以文中對血液在3種血泵模型中的曝光時間和切應力分布分別進行分析,研究其對溶血值的影響.

2.2.2不同結構血泵的紅細胞曝光時間

由冪律公式可知,曝光時間t和切應力是溶血評定的重要參數.曝光時間越長,切應力越大,對血泵溶血性能越不利.從血泵入口選取1 163條跡線檢測紅細胞流經血泵的曝光時間,運用數理統計方法,通過眾數反映試驗所取數據的一般水平,由此得到每種血泵紅細胞曝光時間的大致范圍.圖10所示為跡線數所對應曝光時間眾數在其附近密集區域80%內的區間,非折邊不等間距葉輪血泵、折邊等間距葉輪血泵與折邊不等間距葉輪血泵的曝光時間眾數所在區間分別為(0.060 2,0.200 1),(0.060 6,0.180 7)與(0.110 9,0.410 0).即折邊不等間距葉輪血泵的紅細胞曝光時間最長,非折邊不等間距葉輪血泵的曝光時間次之,折邊等間距葉輪血泵的曝光時間最短,且曝光時間均在合理范圍內.曝光時間與血泵運行轉速有很大關系,轉速越高則血液流過血泵的速度越快、曝光時間就越短,與文中研究設置的轉速相對應.

圖10 流經不同血泵的紅細胞曝光時間Fig.10 Exposure time of red blood cells flowing through different blood pumps

2.2.3不同結構血泵的壁面切應力分布

由流體力學理論可知,固體與液體之間的切應力要大于同一條件下液體與液體之間的切應力,因此最大切應力通常分布在血液與泵體直接接觸的壁面上[21].圖11為3種血泵的壁面切應力分布云圖,3種血泵壁面切應力均在前導葉和后導葉處較低,在葉輪入口和出口處較高,同時受動靜干涉影響,前導葉葉片后緣與后導葉葉片前緣處也出現了局部切應力較高的現象.

圖11 不同血泵壁面切應力分布云圖Fig.11 Cloud diagrams of shear stress distribution on the wall of different blood pumps

通過圖11可知,非折邊不等間距葉輪血泵、折邊等間距葉輪血泵與折邊不等間距葉輪血泵的最大切應力分別為896,945和512 Pa,同時折邊不等間距葉輪血泵的整體切應力分布值明顯小于其他2種血泵,尤其在葉輪區域顯示出較大的優越性.

根據NIIMI等[22]的研究,當泵內紅細胞所受切應力在0~150 Pa時,紅細胞可以很好地保持其生理形態不受破壞;當所受切應力在150~1 000 Pa且受力時間大于1 s時,紅細胞很可能被破壞;當所受切應力超過1 000 Pa時,盡管暴露時間極短,紅細胞也會被破壞.圖12為葉輪血泵在不同區間段內壁面切應力的占比圖.由圖可知,3種血泵流場中的切應力主要集中在0~200 Pa.比較3種葉輪血泵中大于150 Pa的切應力占比δ可知,非折邊不等間距葉輪血泵、折邊等間距葉輪血泵與折邊不等間距葉輪血泵中的占比分別約為22.76%,21.76%和3.16%,即剪應力在0~150 Pa的占比分別為77.24%,78.24%和96.84%.由此可見,折邊不等間距葉輪血泵的切應力遠低于其他2種葉輪血泵,能夠有效降低對血液的損傷程度.

圖12 不同葉輪血泵在不同區間段內壁面切應力的占比Fig.12 Percentage of wall shear stress in different interval segments of different blood pumps with different impellers

文中采用拉格朗日方法對血泵進行溶血值估算,結果如圖13所示.

圖13 不同血泵溶血指數圖Fig.13 Graph of hemolysis index of different blood pumps

3種血泵的溶血值均能滿足血泵的溶血設計指標,其中折邊不等間距葉輪血泵的溶血指數HI最低,約為4.010 0×10-3%,較非折邊不等間距葉輪血泵下降了近3.50%,較折邊等間距葉輪血泵下降了近12.50%.由此可見,折邊不等間距葉輪血泵大大減少了溶血,提高了血泵的血液相容性.

3 結 論

1) 對比3種血泵內流場流線及切應力分布可知,折邊葉片結構血泵避免了葉頂間隙泄漏流的產生,同時折邊不等間距葉輪血泵減少了葉輪入口及出口處的回流和渦流,且流道內的切應力低于其他2種血泵,流場最穩定,內流損耗相對較小.

2) 軸流血泵中葉輪區域是血泵血液相容性問題的重點關注區域,葉輪入口和出口處均存在較高的速度梯度和切應力.折邊不等間距葉輪在葉輪入口及出口處的速度梯度和切應力較其他2種血泵小,同時折邊不等間距葉輪血泵的內流場壓力是逐漸增大的,可使紅細胞保持良好的生理性質.

3) 3種血泵中,非折邊不等間距葉輪血泵與折邊等間距葉輪血泵在葉輪入口與出口處均存在較大的壁面切應力,且切應力在0~150 Pa的占比分別約為77.24%和78.24%;而折邊不等間距葉輪血泵的壁面切應力較小,在0~150 Pa的占比約為96.84%.

4) 3種血泵中,折邊不等間距葉輪血泵的紅細胞曝光時間最長,非折邊不等間距葉輪血泵的曝光時間次之,折邊等間距葉輪血泵的曝光時間最短,這主要是由于同等工況下轉速不同引起的,但曝光時間均在合理范圍內.

5) 在同等工況下,折邊不等間距葉輪可有效降低血泵的溶血值.其中,折邊不等間距葉輪血泵的溶血指數較非折邊不等間距葉輪血泵下降了近3.50%,較折邊等間距葉輪血泵下降了近12.50%,溶血性能最優.

猜你喜歡
血泵葉頂葉輪
平面葉柵多凹槽葉頂傾斜圓柱孔氣膜冷卻與氣動特性研究
分離渦模擬在葉頂間隙流模擬中的可行性分析
透平動葉多凹槽葉頂氣膜冷卻特性的研究
透平動葉葉頂氣膜冷卻設計方案研究
1.4317 QT2鋼在高能泵葉輪上的應用
基于BP神經網絡的旋轉血泵生理控制
應用石膏型快速精密鑄造技術制造葉輪
離心泵葉輪切割方法
基于CFD/CSD耦合的葉輪機葉片失速顫振計算
Flow field CFD analysis of axial flow blood pump*
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合