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

低頻超聲透皮給藥系統壓電-聲-熱計算模型*

2015-06-13 09:38:06彭瀚旻陳致鈞盧鵬輝冒林麗
振動、測試與診斷 2015年6期
關鍵詞:振動測量系統

彭瀚旻,陳致鈞,盧鵬輝,冒林麗

(1.南京航空航天大學超聲電機國家地方聯合工程實驗室 南京,210016)(2.南京航空航天大學機械結構力學及控制國家重點實驗室 南京,210016)

?

低頻超聲透皮給藥系統壓電-聲-熱計算模型*

彭瀚旻1,2,陳致鈞2,盧鵬輝2,冒林麗2

(1.南京航空航天大學超聲電機國家地方聯合工程實驗室 南京,210016)(2.南京航空航天大學機械結構力學及控制國家重點實驗室 南京,210016)

為分析低頻超聲透皮給藥過程中由超聲空化引起的發熱問題,基于壓電方程、熱平衡方程和聲壓平衡方程,利用COMSOL有限元軟件建立了低頻超聲透皮給藥過程中壓電-聲-熱耦合仿真計算模型。通過理論分析和FLIR熱成像儀的溫度測量實驗,獲得了系統在超聲輸入功率為5.5 W、頻率為21 kHz下的溫度場分布與外表面最大溫度變化曲線。超聲換能器與低頻超聲透皮給藥系統(空氣中)的溫度場分布及其最大溫度值的仿真分析結果都與實驗值一致,即在低頻超聲透皮給藥過程中,藥液中超聲空化會造成聲波大幅衰減,從而把部分聲能直接轉化為熱能,導致藥液溫度快速上升,15 min時外表面最高溫度可達40℃,內部輻射頭下方最大溫度的理論值可達41.3℃,達到了低頻超聲透皮給藥的溫度安全要求(≤42℃)。計算與實驗結果表明,所建立的壓電-聲-熱耦合計算模型可以預測系統溫度變化,可用于超聲作用時間、換能器結構尺寸和材料參數的優化設計。

低頻超聲; 透皮給藥; 超聲換能器; 多場耦合; 有限元

引 言

20世紀末,科學家們發現低劑量的超聲可以提高藥物滲透效率,其物理效應中的機械效應(波的傳播)和熱效應都可以改善生物組織的壓力、張力和血液流動速度,促進組織代謝,刺激神經系統等加速康復的作用。Mitrogotri等[1-2]發表了關于利用更低超聲頻率(20 kHz)在皮膚透皮給藥實驗中驗證了胰島素、γ-干擾素和促紅細胞生成素等蛋白質的低頻超聲透皮給藥方法,對傳統的治療性超聲頻段(1 MHz)和低頻超聲頻段(20 kHz)在4種不同滲透液中進行了對比測試,發現低頻超聲對滲透率的促進作用是治療性超聲的1 000倍[2]。以此為基礎,低頻超聲透皮給藥被應用于藥液滲透、麻醉等領域。Katz等[3]在臨床上對42名健康的志愿者進行了局部麻醉的低頻超聲透皮給藥,使麻醉時間從60 min縮短到5 min。Park等[4]使用超聲陣列結構的換能器(20 kHz)進行超聲透皮滲透,使葡萄糖水平在60 min后下降到72±5 mg/dl。Zhang等[5]使用低頻超聲透皮給藥技術來縮短局部鎮痛起效時間。呂川[6]發現低頻超聲透皮給藥方法可顯著提高豬超長寬比例筋膜皮瓣的毛細血管密度,延長給藥時間,顯著提高皮瓣的成活率。

為了進一步把低頻超聲透皮給藥技術應用于臨床,人們開始關注其對生物組織的安全問題。Boucaud等[7]發現當超聲在溶液中的強度低于2.5 W/cm2時皮膚沒有明顯變化,而強度到達5.2 W/cm2時則造成表皮脫落和真皮上部水腫。長時間且連續在4 W/cm2的超聲強度下,皮膚組織結構也發生改變,如表皮和真皮壞死脫落。Farzaneh等[8]在綜述中總結低頻超聲不能使生物組織的溫度超過42℃,過高的溫度將破壞藥液和生物組織的活性。但是,低頻超聲透皮給藥中的溫升與超聲換能器之間的關系卻鮮有研究,超聲換能器結構尺寸參數與材料的選取對溫度的影響規律未知。在超聲換能器設計領域,梁松等[9]利用COMSOL Multiphysics有限元軟件對換能器壓電耦合和聲場模型進行了分析計算,獲得了液體中的聲壓分布。筆者擬建立低頻超聲透皮給藥過程中壓電-聲-熱多物理場耦合計算數學模型,通過有限元計算獲得給藥系統的溫度分布及溫升規律,為設計滿足溫度安全性能要求的低頻超聲透皮給藥用換能器奠定理論基礎。

1 低頻超聲給藥系統結構

1.1 夾心式超聲換能器結構

低頻超聲透皮給藥系統中可以使用多種不同類型的超聲換能器來產生透皮給藥所需的超聲波,其中以夾心式超聲換能器和貼片式超聲換能器最為普遍。筆者以蘇州工業園區海納科技有限公司的夾心式超聲換能器為例,其結構如圖1所示。其工作原理為:當在壓電陶瓷上施加21 kHz的驅動電壓時,壓電陶瓷的逆壓電效應使整個結構產生共振(縱向振動),變幅桿的前端向外發送出21 kHz的低頻超聲波。

圖1 21 kHz超聲換能器(單位:mm)

1.2 Franz透皮擴散池

為測量藥液在皮膚表層的透皮滲透過程,通常采用Franz透皮擴散池,它是一種模擬皮膚表面藥物滲透過程的系統(即體外實驗)。本研究的低頻超聲透皮給藥系統采用天津正通公司TT-8透皮儀系統,用于模擬和測量皮膚上藥物透皮滲透的過程,系統如圖2所示。Franz透皮擴散池結構及主要尺寸如圖3所示,主要包括供給體、聚四氟乙烯壓蓋、皮膚(人造皮膚或者真實離體皮膚)和接受體等部件。此系統工作原理為:超聲換能器在藥液(圖2中擴散池的黃色液體)中激發出超聲波,超聲波破壞皮膚角質層,打開或加強滲透通道,使擴散池上方供給體中高濃度藥液加速透過皮膚進入下方的接受體內,從而增加接受體中藥液的濃度,達到低頻超聲透皮給藥的目的。整個擴散池系統進入在恒溫水域系統中,溫度一般設定在32~37℃之間(模擬人體體表或體內溫度)。

圖2 低頻超聲透皮給藥系統

圖3 低頻超聲透皮給藥用擴散池結構(單位:mm)

2 壓電-聲-熱耦合模型

2.1 理論模型

低頻超聲透皮給藥過程是一個多物理場耦合的復雜過程,主要包括電場、力場、聲場和熱場的綜合作用,主要的能量轉化關系如下:

1) 電能主要通過壓電陶瓷的逆壓電效應轉化為超聲換能器的振動能(此外還有很小一部分通過介電損耗和電路中的電損耗化為熱能);

2) 一小部分的振動能通過結構的阻尼損耗化為了熱能;

3) 大部分的振動能通過超聲輻射面轉化為藥液中的聲能;

4) 部分聲能又通過液體和皮膚中的介質耗散轉化為熱能。

本研究的理論模型主要考慮壓電、聲、熱幾個物理場的變化關系,通過有限元仿真計算,獲得最終藥液和皮膚內的溫度場分布。根據壓電平衡方程,即可獲得電場與力場之間的耦合關系(詳細可參見文獻[10])。

固體和流體傳熱時溫度場的瞬態方程為

(1)

其中:T為溫度;t為時間;ρ為密度;Cp為常壓熱容,k為導熱系數;為梯度;u為速度場(矢量矩陣);Q為總熱功率密度(總熱源)。

由于考慮了溶液和皮膚的聲衰減,其聲壓平衡方程可以利用復數i表示為

(2)

為建立整個系統壓電-聲-熱之間理論關系,首先建立超聲換能器的熱傳遞模型,認為換能器熱量[11]主要來自結構振動的耗散功率Qml與壓電陶瓷的介電損耗功率Qdl,可表示為

(3)

其中:ηm為結構阻尼損耗因子;ηe為介電損耗因子;1/2Re[Sconj(cS)]為應變能儲能密度函數;ω=2πf為角頻率(f為驅動頻率);εp為介電常數;E為電場強度。

將式(3)代入式(1)中的Q即可獲得振動與熱的耦合關系,并根據初始條件和邊界條件求得溫升。

超聲換能器在藥液中工作時將大部分能量轉化為聲能,聲能到達一定程度后引起聲空化,進一步導致藥液和皮膚溫度上升。其中,空化及皮膚引起的熱量主要來源于聲能產生的平均聲功率(也稱為聲強)I,平面波假設下熱功率密度QI、平均聲功率密度I與聲壓之間的關系[12]為

(4)

(5)

將式(5)與式(3)中的功率密度代入式(1),再考慮初始條件和邊界條件后即可獲得整個系統的溫升分布。其中,邊界條件中的外部空氣自然對流系數[13]為

(6)

(7)

其中:Gr為格拉曉夫準則;Pr為普朗特(Prandtl)準數;h為對流傳熱系數;β=1/T(T為絕對溫度)為容積膨脹系數;g為重力加速度;v為運動黏度;L為定型尺寸(垂直壁和水平比分別為高度和寬度);ΔT為壁與自然環境空氣的溫度差(由實驗值確定);k為空氣導熱系數;Nu=C(GrPr)n為努塞爾準則。

C與n的取值主要取決于層流、紊流、垂直壁和水平壁等,是實驗獲得的經驗值。

2.2 仿真計算

為求解以上壓電-聲-熱耦合方程,筆者以COMSOL有限元商務軟件為基礎,對整個低頻超聲透皮給藥系統進行分析和仿真計算。為便于計算上述復雜的物理場耦合,進行以下主要假設和簡化:

1) 忽略超聲換能器和擴散池對環境的熱輻射,因為實際測量的系統最高溫度與環境溫度差小于20℃,其熱交換過程中的熱傳遞與對流換熱起主要作用,熱輻射的影響可以忽略不計;

2) 擴散池的玻璃內表面為硬聲場邊界,即表面上聲能完全反射不存在吸收;

3) 皮膚的復雜結構簡化為具有一定聲衰減系數的薄層材料;

4) 由于主要考慮擴散池供給體及皮膚的溫度變化規律,為簡化3D模型,減少計算量,忽略擴散池接受體中支管對聲能和傳熱的影響。

仿真計算時的主要結構部件的部分尺寸參數見圖1、圖3和表1,零部件的材料及性能參數[14]見表2,部分參數來源于COMSOL軟件的材料數據庫。在COMSOL建模過程中,超聲換能器本體采用固體力學與靜電耦合(壓電模塊)的計算模塊,換能器與溶液邊界采用聲-結構邊界耦合(固體力學與壓力聲學耦合)的計算模塊,固體傳熱模塊則包括整個系統,內部空氣和溶液設為流體傳熱,其余為固體傳熱。

表1 零部件主要結構參數尺寸

表2 主要零部件材料及物理參數

Tab.2 Materials and physic parameters of main parts

零件名稱材料密度/(kg·m-3)導熱系數/(W·(m·K)-1)常壓熱容/(J·(kg·K)-1)壓電陶瓷PZT-87650.02.1420.0后擋蓋鈦合金4428.77.1545.4變幅桿鈦合金4428.77.1545.4套筒鋁合金2730.0155.0893.0螺栓鉻鉬合金鋼7850.044.5475.0溶液(供給)水997.60.624177.2皮膚皮膚1109.00.373391.0溶液(供給)水997.60.624177.2空氣空氣1.1470.02671006.0

初始條件如下:空氣中超聲換能器單獨溫升時起始溫度為26℃,環境溫度為25℃,壓電陶瓷施加電壓為58 Vpp(峰-峰值);超聲換能器作用在Franz擴散池上,起始溫度為26℃,環境溫度為25℃,電壓為83 Vpp。(以上值都與實驗時取值一致,有效功率都約為5.5 W)。

邊界條件如下:超聲換能器的鋁合金外壁設為固定邊界;溶液、皮膚與玻璃內壁接觸部分簡化為硬聲場邊界,與空氣接觸部分簡化為軟聲場邊界;整個系統外部設為自然對流邊界。

在計算壓電-聲-熱耦合場過程中,壓電陶瓷、螺栓及鋁合金的的結構阻尼損耗因子分別設為0.011,0.002和0.002,忽略鈦合金的結構阻尼損耗因子(鈦合金損耗因子遠小于鋼和鋁合金的值),根據式(3)即可獲得機械損耗和介電損耗。

溶液中由于超聲空化產生的氣泡會對聲波傳遞產生反射、散色、介質黏滯、熱傳導等因素[15-16],造成聲能的衰減,其溶液(供給體)平均聲衰減系數根據實驗值分別等效為25 NP/m[16]與20 NP/m,聲速取1 500 m/s。由于皮膚很薄且內部空隙一般小于低頻空化氣泡的直徑,故忽略皮膚內空化引起的聲衰減,認為皮膚的溫度主要由上下溶液的空化熱量傳遞而來。把上述等效衰減系數代入式(5)即可求得聲能在溶液內產生的熱能密度。

整個系統與外界的對流換熱可由式(6)和式(7)求得,ΔT對應15min時實驗測量的垂直壁與水平壁的溫度與環境溫度差分別為9.4 K與1.3 K,代入結構與材料物理參數后求得垂直壁與水平壁的h分別為4.6 W/(m2·K)與2.3 W/(m2·K)。

3 測試系統

為驗證壓電-聲-熱耦合模型的有效性和準確性,筆者采用了德國Polytec公司PSV-500型激光多普勒測振儀和美國的FLIR i7型熱成像儀(分辨率為0.1℃),分別測量超聲換能器變幅桿前端超聲輻射面的振動分布以及系統在15 min內的熱場分布。振動與溫度測量時,超聲換能器的輸入功率都維持在5.5 W左右,工作頻率約為21 kHz。其中,FLIR熱成像測量系統如圖4所示,熱成像儀與超聲換能器都被固定不可移動,每分鐘記錄一次熱成像數據,獲得整個低頻超聲透皮給藥系統外表面的溫度。實際上,Franz擴散池給藥時是在恒溫水浴系統中(圖2),但是由于水溫會阻斷擴散池下方的溫度測量,故筆者采用整個系統在空氣中的測量方法,僅用于本計算模型的測評。

圖4 FLIR熱成像測量系統

4 計算與實驗結果討論與分析

根據第2節中的理論和計算方法,在58 Vpp激勵下超聲換能器空氣中的振動幅值分布與速度分布由PSV-500型激光多普勒測振儀測量所得,結果如圖5(a)所示,其輻射面最大振幅與速度分別達到16.34 μm和2.16 m/s,與實驗測量值(見圖5(b))接近。因此,COMSOL中利用壓電方程獲得關于壓電耦合的振動場分布,特別是輻射面的振動分布與實際測量情況基本一致,從而驗證了壓電-聲-熱耦合計算中壓電耦合模塊的計算方法的有效性。

圖5 超聲換能器振動仿真計算與測量結果(空氣中)

結合式(3),利用壓電耦合模塊中計算所得的儲能密度函數就可以計算出機械損耗與介電損耗的能量密度,再代入式(1),并結合初始條件和邊界條件,即可計算出不同時間段內其內部及表面溫度分布場結果,如圖6(a)所示;FLIR i7測量的不同時間下超聲換能器表面溫度場分布如圖6(b)所示;表面溫度最大值隨時間的變化關系如圖6(c)所示。從圖中可以看出,不同時刻測量所獲得的表面溫度分布與仿真計算結果類似,溫度測量值則稍大于仿真計算結果,這是由于仿真時沒有考慮電路內部的電損耗引起的發熱等因素。因此,壓電-熱耦合計算方法適用于此類型超聲換能器。根據圖6(a),當超聲換能器空載時(或者近似認為空氣是負載),5.5 W的輸入功率通過結構的機械損耗和壓電材料的介電損耗轉化為自身的熱量,且機械損耗效果遠大于介電損耗,其中鈦合金(TC4)的機械損耗非常小,可忽略不計。此外,從圖6(c)可以看出,在15 min內,換能器的溫升并未達到穩態,還呈現遞增趨勢,但是超聲透皮給藥過程一般很短,通常幾分鐘至十幾分鐘即打開皮膚滲透通道,然后撤去超聲讓藥液自由滲透。因此,筆者只考慮超聲作用的前15 min對藥液的溫度影響。

圖6 超聲換能器溫度場仿真計算與測量結果(空氣中)

考慮超聲換能器在Franz擴散池透皮給藥工作中產生的熱量,先通過上述方法,獲得壓電耦合場的振動計算結果,再結合聲-結構邊界與式(2)獲得溶液及皮膚內的聲壓p分布,然后把p代入式(4)和式(5)獲得溶液和皮膚產生的熱量,最后通式(1)獲得整個系統的溫度分布,如圖7(a)所示。不同時間下超聲換能器表面溫度場分布如圖7(b)所示,表面溫度最大值隨時間的變化關系如圖7(c)所示。從圖中可以看出,外表面溫度場分布和溫升值的計算結果與實驗測量值基本一致,說明筆者提出的壓電-聲-熱多場耦合計算模型適用于低頻超聲透皮給藥系統,可以近似滿足預測要求。對比圖6(a)和圖7(a)可知,當溶液成為超聲換能器負載時,機械損耗與介電損耗所產生的熱量大幅下降。如圖7(a)中壓電陶瓷與螺栓位置處的溫度所示,而整個系統溫度最高處在超聲輻射頭下方和皮膚上方,這是由于聲能通過超聲空化轉化為熱量的結果,說明整個系統同樣為5.5 W的輸入功率時,能量主要是通過聲能中的聲衰減轉化為溶液和皮膚的熱量而耗散了,而不是如圖6(a)中所示結構振動的機械損耗占據主要地位。同時,15 min計算所得的最高溫度在皮膚上方,可達41.3℃,幾乎達到透皮給藥溫度的安全閥值[8](≤42℃)。從仿真結果可以得出,此換能器的輸入功率在不超過5.5W、工作時間不超過15 min時,其給藥操作是有效的。

圖7 超聲換能器及Franz擴散池溫度場仿真計算與測量結果(空氣中)

5 結束語

筆者基于COMSOL有限元軟件,對低頻超聲透皮給藥系統建立了壓電-聲-熱多場耦合模型,仿真計算出系統在空氣中工作15 min內的瞬態溫升,擴散池外表面溫升計算結果與實驗測量結果相差5%。超聲換能器在一定條件下滿足低頻超聲透皮給藥的溫度安全要求(≤42℃),表明所述模型可用于低頻超聲透皮給藥系統及其壓電換能器的優化設計。此外,通過模型分析發現,超聲給藥過程中,熱量主要集中在輻射頭下方,主要來源于超聲空化造成聲波衰減而產生的溫升,在結構設計時應予以足夠關注。

[1] Mitragotri S, Blankschrein D, Langer R. Ultrasound-mediated transdermal protein delivery[J]. Science. 1995, 269(5225): 850-853.

[2] Mitragotri S, Blankschrein D, Langer R. Transdermal drug delivery using low-frequency sonophoresis[J]. Pharmaceutical Research, 1996, 13(3): 411-420.

[3] Katz N, Shapiro D, Herrmann T, et al. Rapid onset of cutaneous anesthesia with EMLA cream after pretreatment with a new ultrasound-emitting device[J]. Anesthesia and Analgesia, 2004, 98(2): 371-376.

[4] Park E, Werner J, Smith N B. 3B-1 noninvasive insulin delivery in large pigs (100 lbs) using the lightweight cymbal array[C]∥Ultrasonics Symposium. New York, NY, United States: Institute of Electrical and Electronics Engineers Inc., 2007.

[5] Zhang Guoliang, Lu Rongjun, Chen Xin. Application of low frequency ultrasound to promote topical eutectic mixture of local anesthetics induced analgesia:a randomized controlled study[J]. Chinese Journal of Clinical Rehabilitation, 2004, 8(17): 3364-3365.

[6] 呂川. VEGF、FGF-2低頻超聲透皮給藥系統對豬超長筋膜瓣應用的初步研究[D]. 上海:第二軍醫大學, 2009.

[7] Boucaud A, Montharu J, Machet L, et al. Clinical, histologic, and electron microscopy study of skin exposed to low-frequency ultrasound[J]. Anatomical Record, 2001, 264(1): 114-119.

[8] Ahmadi F, Mcloughlin I V, Chauhan S, et al. Bio-effects and safety of low-intensity, low-frequency ultrasonic exposure[J]. Progress in Biophysics and Molecular Biology, 2012, 108(3): 119-138.

[9] 梁松,張義民. 超聲清洗換能器設計及性能分析[J]. 振動、測試與診斷, 2013, 33(S2): 87-90.

Liang Song, Zhang Yimin. Ultrasonic cleaning transducer design and vibration performance research[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(S2): 87-90. (in Chinese)

[10]王矜奉,蘇文斌,王春明,等. 壓電振動理論與應用[M]. 北京:科學出版社, 2011.

[11]Rajesh J T, Benjamin V V, Ramamoorthy V. Heat generation from dielectric loss and vibration using COMSOL multiphysics[C]∥Proceedings of the 2011 COMSOL Conference in Bangalore. Bangalore, India:[s.n.], 2011.

[12]杜功煥,朱哲民,龔秀芬. 聲學基礎[M]. 3版. 南京: 南京大學出版社, 2012: 126.

[13]張興中,黃文,劉慶國. 傳熱學[M]. 北京:國防工業出版社, 2011: 99.

[14]王榮津. 水聲材料手冊[M]. 科學出版社, 1983:142-149.

[15]Dahlem O, Reisse J, Halloin V. The radially vibrating horn: a scaling-up possibility for sonochemical reactions[J]. Chemical Engineering Science, 1999, 54(13-14): 2829-2838.

[16]Moran C M, Bush N L, Bamber J C. Ultrasonic propagation properties of excised human skin[J]. Ultrasound in Medicine & Biology, 1995, 21(9): 1177-1190.

歡迎訂閱《振動、測試與診斷》

《振動、測試與診斷》由工業和信息化部主管,南京航空航天大學和全國高校機械工程測試技術研究會聯合主辦,是反映振動、動態測試及故障診斷學科領域的科研成果及其應用情況的技術性刊物。主要刊登國內外以振動測試與故障診斷為中心的動態測試理論、方法和手段的研究及應用方面的技術文獻,包括實驗測試技術、測試儀器的研制、方法和系統組成、信號分析、數據處理、參數識別與故障診斷以及有關裝置的設計、使用、控制、標定和校準等,不拘泥于行業和測試項目。

本刊為EICompendex數據庫收錄期刊和中文核心期刊,雙月刊,每逢雙月末出版,每本定價20元,全年120元。歡迎訂閱和投稿,歡迎在本刊刊登各類廣告和科技信息。

編輯部地址:南京市御道街29號 郵政編碼:210016 電話:(025)84893332

傳真:(025)84893332 E-mail:qchen@nuaa.edu.cn 網址:http:∥zdcs.nuaa.edu.cn

10.16450/j.cnki.issn.1004-6801.2015.06.006

*國家自然科學基金青年基金資助項目(51405224);江蘇省科技計劃青年基金資助項目(BK20140818);中央高?;究蒲袠I務費專項資金資助項目(NJ20140027);江蘇省大學生創新創業訓練計劃資助項目(201510287010Y)

2015-06-26;

2015-07-20

TH113.1; TJ02; Q819

彭瀚旻,男,1984年4月生,講師。主要研究方向為壓電換能器理論分析、設計及應用。曾發表《Model study of IPMC beam response based on root deformation》(《Journal of Wuhan University of Technology-Mater》2013,Vol.28,No.1)等論文。 E-mail:penghm@nuaa.edu.cn

猜你喜歡
振動測量系統
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
中立型Emden-Fowler微分方程的振動性
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
主站蜘蛛池模板: 国内精品免费| 巨熟乳波霸若妻中文观看免费 | 久久综合九色综合97婷婷| 国产不卡在线看| 久久综合色88| 亚洲精品男人天堂| 亚洲欧美日韩久久精品| 欧美第二区| 伊人激情久久综合中文字幕| 国产成人永久免费视频| 国产黄色免费看| 中文字幕在线播放不卡| 尤物午夜福利视频| 综合久久久久久久综合网| 国产午夜人做人免费视频| 精品人妻一区二区三区蜜桃AⅤ| 99久久精品免费看国产免费软件| 久久精品嫩草研究院| 国产精品丝袜视频| 日韩精品一区二区深田咏美| 99精品视频九九精品| 波多野结衣的av一区二区三区| 欧美国产另类| 无码有码中文字幕| 亚洲欧美人成电影在线观看| 欧美精品成人| 免费在线a视频| 国产日韩欧美一区二区三区在线 | 久久久久无码精品| 1024国产在线| 免费99精品国产自在现线| 亚洲AV无码乱码在线观看裸奔 | 日韩色图在线观看| 久久久噜噜噜久久中文字幕色伊伊| 四虎成人精品在永久免费| 免费网站成人亚洲| 欧美精品亚洲二区| 久久精品国产免费观看频道| 91毛片网| 久久久久亚洲精品成人网| 最新国产精品第1页| 亚洲日韩Av中文字幕无码| 国产乱人伦AV在线A| 55夜色66夜色国产精品视频| 毛片免费在线视频| 久久久91人妻无码精品蜜桃HD| 国产成人综合网| 亚洲综合香蕉| 中文字幕一区二区视频| 国产草草影院18成年视频| 久草国产在线观看| 国内精品久久久久久久久久影视 | 欧美一级高清免费a| 国产成人一级| 国产91透明丝袜美腿在线| 国产区网址| 丰满人妻久久中文字幕| 国产在线拍偷自揄拍精品| 欧美精品伊人久久| 日韩精品免费在线视频| 一本久道久久综合多人| 亚洲综合专区| 国产精品成人久久| 无码国内精品人妻少妇蜜桃视频| 国产在线精品网址你懂的| 国产成人无码AV在线播放动漫| 四虎影视国产精品| 国产亚洲精久久久久久久91| 国产永久在线视频| 亚洲第一页在线观看| 国产在线观看成人91| 欧美午夜在线观看| 8090午夜无码专区| 日韩成人免费网站| 欧美自慰一级看片免费| 热久久国产| 国产麻豆va精品视频| 宅男噜噜噜66国产在线观看| 伊人久久大香线蕉综合影视| 99国产精品一区二区| 国产小视频在线高清播放| 国产精品乱偷免费视频|