999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于CTA的椎動脈瘤血流動力學分析

2020-10-09 11:01:23劉芳趙思淇盧立彬
軟件 2020年8期
關鍵詞:分析

劉芳 趙思淇 盧立彬

摘 ?要: 分析椎動脈的動脈瘤血流動力學指標在動脈瘤發生、發展及治療后的作用,判斷引起動脈瘤發生與治療后復發的特定血流動力學因素,并為動脈瘤的預防、治療提供理論依據。選取一例顱內動脈瘤患者的CTA影像數據三維建模和仿真計算獲取血流動力學指標:time average wall shear stress、time average wall shear stress grade、oscillatory shear index 、aneurysm formation index、relative retention time等參數作為觀察指標分析。結果顯示:1. 動脈瘤TAWSS及TAWSSG的不穩定,栓塞手術能夠降低動脈瘤破裂的風險性,而在栓塞術后的血管交叉處,血管壁則較易受損;2. OSI值較高,改變瘤體內震蕩水平導致血流紊亂,OSI值減低,血流趨于穩定;3. 隨著動脈瘤AFI值逐漸升高,血液流動可逐漸平穩,可降低動脈瘤破裂危險性。

關鍵詞: 顱內動脈瘤;計算流體力學

中圖分類號: TP319 ? ?文獻標識碼: A ? ?DOI:10.3969/j.issn.1003-6970.2020.08.028

本文著錄格式:劉芳,趙思淇,盧立彬,等. 基于CTA的椎動脈瘤血流動力學分析[J]. 軟件,2020,41(08):97-102

【Abstract】: To analyze the role of vertebral artery aneurysm parameters in the occurrence, development and post treatment of aneurysms, ? ?to determine the specific Hemodynamics that cause the occurrence and recurrence of aneurysms after treatment, and to provide a theoretical basis for the prevention and treatment of aneurysms. Three-dimensional modeling and simulation of CTA images of a patient with cerebral aneurysm were used to obtain Hemodynamics parameters: Time average wall shear stress, time average wall shear stress, Oscillatory Shear Index, Eurasian SM Formation Index and relative retention time. 1. The instability of TAWSS and Tawssg, embolization can reduce the risk of aneurysm rupture, but the vessel wall is more vulnerable at the cross-section after embolization. 2. Osi Value is higher, changes in the level of turbulence in the tumor lead to blood flow disorder, Osi value is reduced, blood flow tends to be stable; 3. With the gradual increase of AFI, the blood flow could be stabilized and the risk of aneurysm rupture could be reduced.

【Key words】: Cerebral aneurysm; Computational fluid dynamics

0 ?引言

顱內動脈瘤(Intracranial aneurysm,IA)是多種因素導致的動脈壁的異常瘤樣擴張,常發生于顱內大動脈的分叉及彎曲處,破裂會導致蛛網膜下腔出血,具有極高的致死率及致殘率。影響動脈瘤生長和破裂的因素主要包括先天生理性、病理性及血流動力學因素。對于IA的治療主要包括開顱夾閉術及血管內介入栓塞兩種方法[1-3]。無論是哪種方法,由于動脈瘤自身的復雜性及不完全閉塞的發生,即使是有經驗的臨床醫生,術后的復發率依舊很高[4-7]。栓塞手術由于創傷小,操作相對簡單等優勢和栓塞材料及技術的不斷發展進步逐漸被廣泛應用于臨床。但同時由于彈簧圈具有可壓縮性使得動脈瘤復發的可能性大大增加,有研究表明栓塞程度是動動脈瘤復發的重要影響因素,Brzegowy等人回顧性分析破裂與未破裂前交動脈瘤的栓塞治療,同樣得出影響顱內動脈瘤復發的最大因素就是栓塞密度,初始栓塞的不完全極易引起動脈瘤的復發[8-9]。栓塞程度低,彈簧圈會隨著血流的沖擊逐漸壓縮,進而向遠側移位、復發。趙慶平等提出瘤腔內的血流速度與瘤腔大小呈負相關,即栓塞程度越低,腔內血流速度加快時,血液對壁面產生的力就可能導致動脈瘤的復發[10]。近年來隨著計算機的發展及有限元軟件的開發,尤其是計算流體動力學數值模擬方法的應用,使得血流建模能更好解釋血流動力學在IA發病機制中的作用[11-16]

1 ?材料與方法

原始影像數據采集:采集解放軍第78集團軍醫院一椎動脈瘤患者CTA影像數據,男性患者,62歲,患者主訴頭部持續頭痛,臨床表現為:行走不穩三個月。經臨床診斷為頭部椎最動脈瘤。經患者本人知情同意并簽署意見書與醫院倫理委員會批準。

圖像后處理工作站:DELL圖像向工作站:DELL 7810/CPU E5/16G內存/英偉達K2200顯卡;

圖像后處理軟件:醫學交互式影像控制系統(Materialis Interative Medical Image Control System,MIMICS,比利時Materialise公司)、醫學建模軟件3-matic medical(比利時Materialis公司);

計算機仿真軟件:ANSYS 19.2:流體仿真軟件CFX,網格劃分軟件FLUENT MESHING

計算結果分析軟件:ENSIGHT10.6。

椎動脈瘤三維重建:將頭部影像DICOM數據導入MIMICS軟件,使用MIMICS分割工具:閾值分割等算法,最后三維計算生成動脈瘤三維初步模型以stl格式導入3-matic medical軟件中,使用光順表面、去除細小分支、切好出、入口平面,最后形成動脈瘤三維模型,如圖1所示。

網格劃分:由于模型結構復雜,使用非結構化的四面體網格劃分,為保證計算精度,在動脈瘤管壁進行五層加密。

邊界條件:本計算不考慮能量的傳遞,不考慮重力。血液密度為1056 kg/m3,動力粘度為0.0035 。計算采用瞬態計算,兩個入口,兩個出口,壁面無滑移。入口采用速度入口,出口采用壓力出口。為保證盡快收斂,入口速度采用極小的速度差。出口采用壓力出口,壓力曲線如圖3所示。

血流作用在內皮細胞上的力的血流動力學參數GON,壁面切向和正交向上的向量,如果空間梯度G數值變化,代表對內皮細胞產生震蕩張力和壓縮力,在一個心動周期內,如果某個點發生較大的梯度變化,單位面積內發生強烈的震蕩張力或者壓縮力作用于內皮細胞上,GON是用來量化震蕩張力和壓縮力的程度。

CFX無法直接實現上述參數指標,使用CFX ccl語言編程,如下為子程序的部分內容:

IBRARY:

CEL:

EXPRESSIONS:

DOMAIN: FLUIDdom

Coord Frame = Coord 0

Domain Type = Fluid

Location = Assembly

BOUNDARY: INLET1

Boundary Type = INLET

Location = INLET1

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Normal Speed = invel1

Option = Normal Speed

END

END

END

BOUNDARY: INLET2

Boundary Type = INLET

Location = INLET2

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Normal Speed = invel2

Option = Normal Speed

END

END

END

BOUNDARY: OUTLET1

Boundary Type = OPENING

Location = OUTLET1

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW DIRECTION:

Option = Normal to Boundary Condition

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Option = Opening Pressure and Direction

Relative Pressure = OUTLET1f

END

END

END

BOUNDARY: OUTLET2

Boundary Type = OPENING

Location = OUTLET2

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW DIRECTION:

Option = Normal to Boundary Condition

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Option = Opening Pressure and Direction

Relative Pressure = OUTLET2f

END

END

END

BOUNDARY: OUTLET3

Boundary Type = OPENING

Location = OUTLET3

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW DIRECTION:

Option = Normal to Boundary Condition

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Option = Opening Pressure and Direction

Relative Pressure = OUTLET3f

END

END

END

BOUNDARY: OUTLET4

Boundary Type = OPENING

Location = OUTLET4

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW DIRECTION:

Option = Normal to Boundary Condition

END

FLOW REGIME:

Option = Subsonic

END

MASS AND MOMENTUM:

Option = Opening Pressure and Direction

Relative Pressure = OUTLET4f

END

END

END

BOUNDARY: WALL_VESSEL

Boundary Type = WALL

Location = WALL_PARENT_VESSEL

BOUNDARY CONDITIONS:

ADDITIONAL VARIABLE: WSSField

Additional Variable Value = WallShearMag

Option = Value

END

ADDITIONAL VARIABLE: WSSxF

Additional Variable Value = Wall Shear X

Option = Value

END

ADDITIONAL VARIABLE: WSSyF

Additional Variable Value = Wall Shear Y

Option = Value

END

ADDITIONAL VARIABLE: WSSzF

Additional Variable Value = Wall Shear Z

Option = Value

END

MASS AND MOMENTUM:

Option = No Slip Wall

END

END

END

DOMAIN MODELS:

BUOYANCY MODEL:

Option = Non Buoyant

END

2.2 ?TAWSS云圖分析

從平均壁面切應力(TAWSS)的云圖(圖4)上來看,云圖顏色越偏紅,代表平均壁面切應力越大,越接近藍色則平均壁面切應力越小。載瘤動脈由于其具有較高的流速,因此載瘤動脈的平均壁面切應力要大于動脈瘤。

2.3 ?動脈瘤的TAWSSG云圖分析

從TAWSSG云圖中(圖5)可以看出,瘤體的TAWSSG開始時低于載瘤動脈,而后逐漸接近。可能是由于開始時動脈瘤體積較大[18-20],血液流速較慢,壁面切應力數值變化不變明顯,所示動脈瘤偏藍色,TAWSSG值較低。

2.4 ?動脈瘤的OSI云圖分析

圖6為動脈瘤OSI云圖,在動脈瘤頂端存在小部分高OSI區域。動脈瘤頂端的高OSI區域較之前擴大,振蕩剪切系數代表的是整個心動周期內壁面切應力方向變化快慢的量,OSI不同是反應震蕩水平,即流動的強度和方向的改變,越大表示震蕩越強,流體在周期內流動的方向不穩定,導致動脈瘤內的血流運動趨于紊亂。

2.5 ?動脈瘤AFI云圖分析

圖7為動脈瘤AFI云圖,瘤體側壁上存在部分AFI低區域,即偏藍色區域。流增多形成渦流并不斷的沖擊著動脈瘤管壁,壁面切應力的方向變化明顯,血液流動不穩定。

2.6 ?動脈瘤GON云圖分析

動脈瘤GON云圖(圖8)表明在動脈瘤表面存在強烈的震蕩力和壓縮力,原因是血液在進入瘤腔后形成渦流,導致動脈瘤壁震蕩,這種沖擊對瘤壁造成膨脹或者擴張。

3 ?討論與結論

通過血流動力學計算分析,發現完全栓塞手術可以阻斷進入動脈瘤內的血液[21-24],提高TAWSS及降低OSI等,降低了破裂出血的風險。通過本實驗對最動脈瘤血流動力學的參數的變化分析可得出以下結論:

(1)動脈瘤TAWSS及TAWSSG的不穩定,栓塞手術能夠降低動脈瘤破裂的風險性[25-30],而在栓塞術后的血管交叉處,血管壁則較易受損;

(2)OSI值較高,改變瘤體內震蕩水平導致血流紊亂,OSI值減低,血流趨于穩定;

(3)隨著動脈瘤AFI值逐漸升高,血液流動可逐漸平穩,可降低動脈瘤破裂危險性。

參考文獻

[1] Wang Hua-Wei,Sun Zheng-Hui,Wu Chen et al. Surgical management of recurrent aneurysms after coiling treatment[J]. Br J Neurosurg, 2017, 31: 96-100.

[2] Phan Kevin,Huo Ya R,Jia Fangzhi et al. Meta-analysis of stent-assisted coiling versus coiling-only for the treatment of intracranial aneurysms[J].J Clin Neurosci, 2016, 31: 15-22.

[3] Gawlitza Matthias,Soize Sebastien,Januel Anne-Christine et al. Treatment of recurrent aneurysms using the Woven EndoBridge (WEB): anatomical and clinical results[J]. J Neurointerv Surg, 2018, 10: 629-633.

[4] Huang De-Zhang,Jiang Bin,He Wei et al. Risk factors for the recurrence of an intracranial saccular aneurysm following endovascular treatment[J]. Oncotarget, 2017, 8: 33676-33682.

[5] Yu Le-Bao, Fang Zhi-Jun, Yang Xin-Jian, et al. Management of Residual and Recurrent Aneurysms After Clipping or Coiling: Clinical Characteristics, Treatments, and Follow-Up Outcomes[J].World Neurosurg, 2019, 122: e838-e846.

[6] Zhang Donghuan, Wang Honglei, Liu Tianyi et al. Re-Recurrence of Intracranial Aneurysm with Proximal Vascular Stenosis After Primary Clipping and Secondary Endovascular Embolization: A Case Report and Literature Review[J]. World Neurosurg, 2019, 121: 28-32.

[7] Kim S-T, Baek J W, Jin S-C et al. Coil Embolization in Patients with Recurrent Cerebral Aneurysms Who Previously Underwent Surgical Clipping[J]. AJNR Am J Neuroradiol, 2019, 40: 116-121.

[8] Brzegowy Pawe?, Kucyba?a Iwona, Krupa Kamil, et al. Angiographic and clinical results of anterior communicating artery aneurysm endovascular treatment[J]. Wideochir Inne Tech Maloinwazyjne, 2019, 14: 451-460.

[9] Chueh Ju-Yu, Vedantham Srinivasan, Wakhloo Ajay K, et al. Aneurysm permeability following coil embolization: packing density and coil distribution[J]. J Neurointerv Surg, 2015, 7: 676-81.

[10] Zhao Qingping, Chen Guangzhong, Li Tielin, Zhao Wei, Yuan Yuan, Feng Yanqiu. Hemodynamic analysis of embolization density and recurrence of intracranial aneurysm [J]. Chinese Journal of neuropsychiatric diseases, 2013,39 (06): 339-343

[11] Liang Li,Steinman David A,Brina Olivier et al. Towards the Clinical utility of CFD for assessment of intracranial aneurysm rupture - a systematic review and novel parameter- ranking tool[J]. J Neurointerv Surg, 2019, 11: 153-158.

[12] Valen-Sendstad Kristian,Bergersen Aslak W,Shimogonya Yuji et al. Real-World Variability in the Prediction of Intracranial Aneurysm Wall Shear Stress: The 2015 International Aneurysm CFD Challenge[J]. Cardiovasc Eng Technol, 2018, 9: 544-564.

[13] Roloff Christoph, Stucht Daniel, Beuing Oliver et al. Comparison of intracranial aneurysm flow quantification techniques: standard PIV vs stereoscopic PIV vs to-mographic PIV vs phase-contrast MRI vs CFD[J]. J Neurointerv Surg, 2019, 11: 275-282.

[14] Frolov S V,Sindeev S V,Liepsch D et al. Experimental and CFD flow studies in an intracranial aneurysm model with Newtonian and non-Newtonian fluids[J]. Technol Health Care, 2016, 24: 317-33.

[15] Botti Lorenzo, Paliwal Nikhil, Conti Pierangelo et al. Modeling hemodynamics in intracranial aneurysms: Comparing accuracy of CFD solvers based on finite element and finite volume schemes[J]. Int J Numer Method Biomed Eng, 2018, 34: e3111.

[16] Xiang J, Tutino V M, Snyder K V et al. CFD: computational fluid dynamics or confounding factor dissemination? The role of hemodynamics in intracranial aneurysm rupture risk assessment[J]. AJNR Am J Neuroradiol, 2014, 35: 1849-57.

[17] Resnick Nitzan, Yahav Hava, Shay-Salit Ayelet et al. Fluid shear stress and the vascular endothelium: for better and for worse[J]. Prog. Biophys. Mol. Biol., 2003, 81: 177-99.

[18] Davies Peter F, Hemodynamic shear stress and the endothelium in cardiovascular pathophysiology[J]. Nat Clin Pract Cardiovasc Med, 2009, 6: 16-26.

[19] Shojima Masaaki, Oshima Marie, Takagi Kiyoshi et al. Magnitude and role of wall shear stress on cerebral aneurysm: computational fluid dynamic study of 20 middle cerebral artery aneurysms[J]. Stroke, 2004, 35: 2500-5.

[20] Boussel L, Rayz V, Mcculloch C, et al. Aneurysm Growth Occurs at Region of Low Wall Shear Stress: Patient-Specific Correlation of Hemodynamics and Growth in a Longitudinal Study[J]. Stroke, 2008, 39(11): 2997-3002.

[21] Reneman Robert S, Arts Theo, Hoeks Arnold P G, Wall shear stress-an important determinant of endothelial cell function and structure-in the arterial system in vivo. Discrepancies with theory[J]. J. Vasc. Res., 2006, 43: 251-69.

[22] Yeow Siang Lin,Leo Hwa Liang,Is Multiple Overlapping Uncovered Stents Technique Suitable for Aortic Aneurysm Repair?[J].Artif Organs, 2018, 42: 174-183.

[23] Suzuki Tomoaki, Stapleton Christopher J, Koch Matthew J et al. Decreased wall shear stress at high-pressure areas predicts the rupture point in ruptured intracranial aneurysms[J].J. Neurosurg., 2019, undefined: 1-7.

[24] Omodaka Shunsuke, Sugiyama Shin-Ichirou, Inoue Takashi et al. Local hemodynamics at the rupture point of cerebral aneurysms determined by computational fluid dynamics analysis[J]. Cerebrovasc. Dis., 2012, 34: 121-9.

[25] Dabagh Mahsa, Nair Priya, Gounley John et al. Hemodynamic and morphological characteristics of a growing cerebral aneurysm[J]. Neurosurg Focus, 2019, 47: E13.

[26] Zhang Yisen, Jing Linkai, Zhang Ying, et al. Low wall shear stress is associated with the rupture of intracranial aneurysm with known rupture point: case report and literature review[J]. BMC Neurol, 2016, 16: 231.

[27] Liu Xiuxian. Hemodynamic study of surgical treatment of intracranial aneurysm[D]. Southeast University, 2017.

[28] Chang Yu, Liu Chang, Zhang Qi, Shi Yue, Gao Bin. Hemodynamic study on the influence of Subarachnoid Aneurysm on the risk of cerebral aneurysm rupture[J]. Journal of Beijing University of technology, 2017, 43(07): 1079-1085. Computational fluid dynamics analysis[J]. Cerebrovasc. Dis., 2012, 34: 121-9.

[29] Murayama Yuichi, Fujimura Soichiro, Suzuki Tomoaki, et al. Computational fluid dynamics as a risk assessment tool for aneurysm rupture[J]. Neurosurg Focus, 2019, 47: E12.

[30] Sheng Bin, Wu Degang, Yuan Jinlong, et al. Hemodynamic Characteristics Associated With Paraclinoid Aneurysm Recurrence in Patients After Embolization[J]. Front Neurol, 2019, 10: 429.

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗證分析
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發展趨勢分析
經濟危機下的均衡與非均衡分析
對計劃生育必要性以及其貫徹實施的分析
現代農業(2016年5期)2016-02-28 18:42:46
GB/T 7714-2015 與GB/T 7714-2005對比分析
出版與印刷(2016年3期)2016-02-02 01:20:11
中西醫結合治療抑郁癥100例分析
偽造有價證券罪立法比較分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 久久这里只有精品免费| 91麻豆精品视频| 亚洲午夜福利精品无码| 少妇极品熟妇人妻专区视频| 色婷婷亚洲综合五月| 国产美女久久久久不卡| 亚洲欧美自拍视频| 亚洲Av激情网五月天| 99这里只有精品在线| 国产高清在线观看91精品| 国产一区成人| 精品国产成人三级在线观看| 在线无码九区| 亚洲美女一级毛片| 欧美一区二区丝袜高跟鞋| 亚洲人成电影在线播放| 在线播放精品一区二区啪视频 | 中文字幕乱码中文乱码51精品| 99热这里只有精品免费国产| a级毛片免费播放| 九九热精品视频在线| 欧美成人精品高清在线下载| 美女免费黄网站| 亚洲一区毛片| 五月六月伊人狠狠丁香网| 欧美日韩北条麻妃一区二区| 国产成人盗摄精品| 激情视频综合网| 在线观看国产精美视频| 国产欧美专区在线观看| 国产精品视频a| 欧美日韩一区二区在线免费观看| 免费在线色| 在线看免费无码av天堂的| 国产福利免费在线观看| 久久国产精品嫖妓| 国产鲁鲁视频在线观看| 久久永久精品免费视频| 日本午夜在线视频| 播五月综合| 国产三级成人| 亚洲美女一级毛片| 嫩草在线视频| 亚洲有无码中文网| 97在线公开视频| 亚洲国产一成久久精品国产成人综合| 欧美一区二区人人喊爽| 国产99视频精品免费观看9e| 黄色三级网站免费| 国产91在线|日本| 国产亚洲现在一区二区中文| 最新国产麻豆aⅴ精品无| 在线免费看片a| 国产簧片免费在线播放| 日本少妇又色又爽又高潮| 1769国产精品视频免费观看| 韩日午夜在线资源一区二区| 五月激情婷婷综合| 亚洲日本www| 精品少妇人妻av无码久久| a级免费视频| 58av国产精品| 午夜国产不卡在线观看视频| 国产又色又刺激高潮免费看| 亚欧美国产综合| 久久黄色小视频| 免费在线a视频| 狠狠v日韩v欧美v| 国产欧美专区在线观看| 欧美激情福利| 中文字幕在线不卡视频| 国产十八禁在线观看免费| 国产精品成人AⅤ在线一二三四 | 国产精品jizz在线观看软件| 日本亚洲最大的色成网站www| 午夜不卡福利| 午夜影院a级片| 日韩国产精品无码一区二区三区| 午夜三级在线| 国产亚洲欧美在线人成aaaa | 国产资源免费观看| 91色老久久精品偷偷蜜臀|