• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    粘彈性本構(gòu)對(duì)人耳動(dòng)力學(xué)特性影響的數(shù)值研究

    2016-01-11 00:35:49田佳彬,饒柱石,塔娜
    振動(dòng)與沖擊 2015年22期
    關(guān)鍵詞:粘彈性有限元分析耳蝸

    粘彈性本構(gòu)對(duì)人耳動(dòng)力學(xué)特性影響的數(shù)值研究

    田佳彬1,饒柱石1,塔娜1,許立富1,黃新生2

    (1.上海交通大學(xué)振動(dòng)、沖擊噪聲研究所 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海200240;2.復(fù)旦大學(xué)附屬中山醫(yī)院,上海200032)

    摘要:為分析中耳軟組織粘彈性材料特性對(duì)人耳系統(tǒng)動(dòng)力學(xué)特性影響,建立包括外耳道、中耳及耳蝸的整耳有限元模型。外耳道及中耳模型用微CT掃描與逆向成型技術(shù)建立,耳蝸采用雙腔導(dǎo)管形式簡(jiǎn)化模型?;谠撃P?,中耳部分軟組織材料屬性采用線性粘彈性,以表征動(dòng)態(tài)分析中能量損耗。在外耳道施加90 dB SPL聲壓模擬聲激勵(lì),并在計(jì)算中考慮外耳道氣體、中耳固體及耳蝸流體多場(chǎng)耦合作用。中耳結(jié)構(gòu)響應(yīng)包括鼓膜臍部與鐙骨底板位移及鐙骨底板速度傳遞函數(shù),耳蝸流體壓力響應(yīng)用于計(jì)算中耳壓力增益、耳蝸輸入聲阻抗及壓力逆向傳遞函數(shù)。結(jié)果表明,考慮粘彈性后,人耳系統(tǒng)動(dòng)態(tài)響應(yīng)參數(shù)較線彈性有一定程度改善,尤其在高頻段提升較明顯,與實(shí)驗(yàn)測(cè)量數(shù)據(jù)匹配效果更好。

    關(guān)鍵詞:粘彈性;動(dòng)態(tài)能耗;耳蝸;動(dòng)力學(xué)建模;多場(chǎng)耦合;有限元分析

    中圖分類號(hào):R318.01文獻(xiàn)標(biāo)志碼:A

    Numerical analysis on the effects of viscoelastic constitutive relation on dynamic characteristics of human ear

    TIANJia-bin1,RAOZhu-shi1,TANa1,XULi-fu1,HUANGXin-sheng2(1. Institute of Vibration, Shock and Noise, State Key Laboratory of Mechanical System and Vibration,Shanghai Jiao Tong University, Shanghai 200240, China; 2. Zhongshan Hospital A?liated to Fudan University, Shanghai 200032, China)

    Abstract:To analyze the effects of viscoelastic properties of middle ear soft tissues on the dynamic characteristics of human ear system, a finite element (FE) model of the human ear consisting of the external ear canal, middle ear and cochlea was developed. The geometric configuration of the external ear canal and middle ear was constructed via micro-CT scanning and reverse engineering technology, and the cochlea was simplified as an uncoiled, two-chambered and fluid-filled duct. The viscoelastic material effect was introduced into the behaviors of middle ear soft tissues to represent the energy dissipation in dynamic analysis. A multiphysics coupled analysis was conducted on the model in which the coupling effects among the air in the ear canal, the fluid in the cochlea and the middle ear structures were concerned. Then a sound pressure of 90 dB SPL was applied on the ear canal to simulate the sound stimulus on normal human ear. Middle ear structural responses such as movements of the tympanic membrane and stapes footplate in response to the sound stimulus were derived by this model. Meanwhile, based on calculating the pressure of the fluid in the cochlea, the sound pressure gain across the middle ear, the cochlear input impedance and the reverse pressure transfer function of cochlea were also obtained. The results show that taking into account the viscoelastic properties of middle ear soft tissues can improve the dynamic responses of the human ear system as compared with the results of a linear elastic model, especially at the high-frequency range. The better agreements between the model results and the experimental data in the literature illustrate the necessity of considering viscoelasticity for dynamic modeling of human ear.

    Key words:viscoelasticity; dynamic energy dissipation; cochlea; dynamic modeling; multiphysics coupling; finite element analysis

    人耳系統(tǒng)主要由耳廓、外耳道、中耳及耳蝸構(gòu)成。該系統(tǒng)具有復(fù)雜幾何形態(tài)、超微型結(jié)構(gòu)特征及非同質(zhì)材料特性,且在正常生理過(guò)程中伴隨聲固耦合、流固耦合等復(fù)雜機(jī)理?;诖耍捎糜邢迒卧ń⑷硕鷶?shù)值模型已成為人耳生物力學(xué)研究的重要方法[1]。

    中耳聲傳導(dǎo)模型主要用于模擬聲音經(jīng)外耳道、鼓膜、中耳聽(tīng)骨鏈至鐙骨底板的傳遞機(jī)制[1-2],其與耳蝸的耦合模型則可用于分析耳蝸內(nèi)部淋巴液的壓力響應(yīng)及基底膜沿長(zhǎng)度方向位移分布[3]等特性。人耳數(shù)值模型與實(shí)驗(yàn)測(cè)量數(shù)據(jù)的匹配程度及聲傳導(dǎo)機(jī)理研究已取得一定成果。力學(xué)實(shí)驗(yàn)[4-5]研究結(jié)果顯示,人耳部分組織與器官如鼓膜、中耳軟組織等均具有粘彈性材料特性,故動(dòng)態(tài)分析中的材料屬性普遍采用線彈性及瑞利阻尼的假設(shè)形式[6],與實(shí)驗(yàn)結(jié)果有出入,不能反映人耳各部分結(jié)構(gòu)真實(shí)力學(xué)特性。而準(zhǔn)確的人耳動(dòng)力學(xué)建模需考慮軟組織粘彈性特性的實(shí)際影響。

    人耳的主要功能在于將耳廓收集的動(dòng)態(tài)聲壓傳遞到耳蝸淋巴液轉(zhuǎn)化成電信號(hào)刺激耳蝸的聽(tīng)神經(jīng)。此動(dòng)態(tài)過(guò)程具有較大的頻率相關(guān)性,并取決于耳蝸基底膜的選頻特性[7],使人耳對(duì)高頻聲壓幅值較敏感。對(duì)具有粘彈性特性材料,動(dòng)態(tài)過(guò)程中能量耗散通常通過(guò)復(fù)模量虛部即損耗模量形式引入,而復(fù)模量則采用粘彈性本構(gòu)模型方式表達(dá)成頻率的函數(shù)。因此,人耳力學(xué)模型中考慮粘彈性材料屬性,必對(duì)人耳系統(tǒng)頻率相關(guān)性產(chǎn)生影響?,F(xiàn)有大部分人耳有限元模型的人耳結(jié)構(gòu)能量損耗通常用瑞利阻尼定義[8],阻尼系數(shù)設(shè)為α=0,β≠0,因而會(huì)導(dǎo)致系統(tǒng)阻尼隨頻率提高顯著增加,降低系統(tǒng)的高頻響應(yīng)。有限元模型中,中耳模型[9]及含耳蝸的整耳模型均出現(xiàn)位移響應(yīng)的高頻幅值相對(duì)較低現(xiàn)象,尤其鐙骨底板的高頻位移明顯低于實(shí)驗(yàn)結(jié)果,故在人耳動(dòng)力學(xué)建模中應(yīng)以更真實(shí)的本構(gòu)關(guān)系取代線彈性假設(shè)條件。

    完整的人耳聲傳導(dǎo)路徑應(yīng)含外耳道氣體、中耳固體、耳蝸流體及基底膜等結(jié)構(gòu)。耳蝸流體壓力響應(yīng)為耳蝸生物力學(xué)特性的重要參數(shù),Stieger等[10]研究中耳蝸底部隔膜壓力差被認(rèn)為耳蝸的輸入信號(hào),為聲傳導(dǎo)驅(qū)動(dòng)因素。基于正、逆向激勵(lì)的耳蝸隔膜壓力差相等方法,圓窗激勵(lì)式人工中耳在理論與實(shí)際應(yīng)用中也取得一定進(jìn)展[11]。因此,研究粘彈性特性引入對(duì)耳蝸流體壓力響應(yīng)等參數(shù)影響尚需建立外耳道、中耳及耳蝸的耦合模型。為此,本文基于耳蝸建模中普遍采用的雙腔導(dǎo)管模型對(duì)耳蝸實(shí)際螺旋形結(jié)構(gòu)進(jìn)行幾何簡(jiǎn)化。而外耳道、中耳的幾何模型基于逆向成型技術(shù)建立,并通過(guò)定義流固耦合面方式實(shí)現(xiàn)外耳、中耳、耳蝸整體集成。在此基礎(chǔ)上,考慮鼓膜、中耳軟組織的粘彈性材料屬性,計(jì)算中耳結(jié)構(gòu)位移及耳蝸流體壓力等動(dòng)力學(xué)參數(shù),比較線彈性與粘彈性響應(yīng)的異同點(diǎn),分析引入粘彈性本構(gòu)在人耳動(dòng)力學(xué)建模中的必要性??蔀榻⑼暾?、準(zhǔn)確的人耳動(dòng)力學(xué)模型提供思路。

    1方法

    1.1整耳網(wǎng)格模型

    用微CT掃描儀(GE Healthcare,eXplore Locus SP)掃描人體顳骨標(biāo)本(男,45歲,右耳),獲得882張醫(yī)療影像圖片,掃描層厚度43.5 μm。將所得圖片文件導(dǎo)入Simpleware軟件進(jìn)行逆向成型,建立外耳道及中耳的三維幾何模型。將模型尺寸與模型數(shù)據(jù)[12]進(jìn)行比對(duì),驗(yàn)證所建模型尺寸與結(jié)構(gòu)特征處于合理范圍內(nèi)。用有限元前處理軟件Hypermesh劃分網(wǎng)格,獲得人體外耳道及中耳模型。其中外耳道氣體采用四面體單元AC3D4,單元數(shù)22 366,中耳結(jié)構(gòu)包括聽(tīng)骨鏈、中耳韌帶與肌腱亦用四面體單元C3D4,總單元數(shù)為48 043。鼓膜、鼓膜環(huán)韌帶采用殼單元S3劃分,單元數(shù)804;鼓膜環(huán)韌帶及鼓膜張緊、松弛部厚度分別取0.2 mm、0.05 mm及0.1 mm。

    耳蝸采用非螺旋狀、充滿液體的雙腔導(dǎo)管模型,包括前庭階、鼓階、基底膜、卵圓窗、圓窗及蝸孔,模型整體尺寸見(jiàn)文獻(xiàn)[6],與人體耳蝸實(shí)際尺寸相近?;啄こ叽缬傻撞亢穸?.5 μm、寬度150 μm線性變化到頂部厚度2.5 μm、寬度500 μm[13],總長(zhǎng)度34 mm,劃分成482個(gè)S3、S4殼單元。前庭階及鼓階的液體容積分別為92.781 μL、93.270 μL,用六面體單元AC3D4及AC3D8劃分網(wǎng)格,單元數(shù)分別為17 577及13 802。圓窗直徑基于文獻(xiàn)[8]圓窗膜面積2.1 mm2確定為1.6 mm。所建人耳有限元模型見(jiàn)圖1、圖2。為顯示耳蝸內(nèi)部的基底膜結(jié)構(gòu),耳蝸流體設(shè)置成半透明形式。中耳幾何參數(shù)及見(jiàn)表1。

    圖1 外耳、中耳與耳蝸耦合模型 Fig.1 Coupled model of human ear consisting of external ear canal, middle ear and cochlea

    1.2材料屬性

    人耳有限元模型材料參數(shù)按外耳道氣體、中耳及耳蝸三部分設(shè)定。外耳道氣體材料屬性與空氣相同,體積模量0.142 MPa,密度1.21 kg/m3。中耳各部分結(jié)構(gòu)采用均勻及各向同性彈性材料,泊松比均為0.3。鼓膜、聽(tīng)骨鏈材料及韌帶、肌腱材料參數(shù)采用交叉調(diào)試[12]方法確定,獲得鼓膜、聽(tīng)骨鏈、中耳韌帶及肌腱材料屬性見(jiàn)表2。動(dòng)態(tài)分析中的能量損耗通過(guò)兩種形式實(shí)現(xiàn),鼓膜張緊部、鼓膜松弛部、砧錘關(guān)節(jié)、砧鐙關(guān)節(jié)、鐙骨環(huán)韌帶及圓窗用線性粘彈性材料模型,動(dòng)態(tài)能耗通過(guò)粘性部分表達(dá),材料參數(shù)見(jiàn)表3。中耳其它部分能量損耗采用瑞利阻尼形式,阻尼系數(shù)α=0 s-1,β=0.75×10-4s。

    圖2 中耳各部分結(jié)構(gòu)說(shuō)明 Fig.2 Illustration of components in the middle ear

    部件尺寸模型尺寸參考尺寸[14]鼓膜沿著錘骨柄方向的直徑/mm8.6678.0-10.0垂直錘骨柄方向的直徑/mm7.5747.5~9.0錐形高度/mm1.7742.0表面積/mm263.77755.8~85.0厚度/mm0.050.04~0.075錘骨錘骨柄底部到外側(cè)突底部的距離/mm4.0425.8總長(zhǎng)度/mm7.4537.6~9.1砧骨砧骨長(zhǎng)突的長(zhǎng)度/mm6.8437.0砧骨短突的長(zhǎng)度/mm4.7385.0鐙骨高度/mm3.4512.5~4.0鐙骨底板的長(zhǎng)度/nn2.7472.64~3.36鐙骨底板的寬度/mm1.4060.7~1.66外耳道長(zhǎng)度(鼓膜臍部到外耳道入口)/mm25.32925~31體積/mm3865.927830~1972

    耳蝸部分前庭階及鼓階流體材料特性與水相同,即體積模量2 250 MPa、密度1 000 kg/m3。為模擬耳蝸基底膜沿長(zhǎng)度方向剛度漸變的特性,用定義溫度場(chǎng)變量方式將基底膜楊氏模量設(shè)為由底部40 MPa線性減小中間15 MPa,再依次線性減小到頂部的3 MPa,瑞利阻尼系數(shù)取α=0 s-1,β=0.1×10-3s?;啄ぶ尾糠峙c圓窗的楊氏模量分別取1.41×104MPa及0.3 MPa。耳蝸流體粘度以間接方式引入,表達(dá)成流體的體積阻力,表達(dá)式為

    (1)

    式中:ω為圓頻率;ρ為耳蝸流體密度;K為耳蝸流體體積模量;η為體積粘度系數(shù),大小為0;μ為剪切粘度系數(shù),大小為1×10-3N·s/m2。

    表2 中耳結(jié)構(gòu)的材料參數(shù)

    表3 中耳軟組織的粘彈性材料參數(shù)

    1.3結(jié)構(gòu)與流體邊界設(shè)定

    人耳有限元模型中結(jié)構(gòu)邊界條件均采用節(jié)點(diǎn)位移固支形式,具體包括:鼓膜環(huán)韌帶環(huán)向節(jié)點(diǎn)、中耳韌帶及肌腱端部節(jié)點(diǎn)、基底膜支撐部分外側(cè)節(jié)點(diǎn)及圓窗外側(cè)周向節(jié)點(diǎn)。模型流體邊界按流體內(nèi)、外邊界分別設(shè)定。模型的流體外邊界設(shè)定中,外耳道氣體毗鄰耳道壁面流體單元表面及耳蝸流體毗鄰耳蝸骨壁的流體單元表面定義成Wall邊界,即流體壓力的法向梯度為零。外耳道的流體單元?jiǎng)t施加固定幅值的壓力載荷用于模擬正常人耳的聲激勵(lì)。對(duì)模型的流體內(nèi)邊界,通過(guò)定義流固耦合面實(shí)現(xiàn)流體壓力的法向梯度與結(jié)構(gòu)力的相互傳導(dǎo)。流固耦合面包括:外耳道氣體與鼓膜外側(cè)面的耦合、鐙骨底板表面與耳蝸流體的耦合、圓窗與耳蝸流體耦合、耳蝸隔膜(含基底膜)兩側(cè)與前庭階流體與鼓階流體的耦合。需指出的是,由于耳蝸流體與鐙骨底板表面及圓窗膜表面間網(wǎng)格節(jié)點(diǎn)并非一一對(duì)應(yīng),故實(shí)際建模時(shí)采用Abaqus中的Tie命令完成各結(jié)合面的連接處理。

    1.4粘彈性本構(gòu)關(guān)系

    為表達(dá)軟組織粘彈性材料特性,采用廣義的線性固體模型或Weichert模型[19]進(jìn)行模擬,見(jiàn)圖3。

    圖3 Weichert粘彈性模型 Fig.3 Weichert model for viscoelastic material

    基于圖3模型,軟組織的時(shí)域松弛模量E(t)可表達(dá)成Prony級(jí)數(shù)[20],即

    (2)

    式中:E∞為t→∞時(shí)長(zhǎng)期模量;Ei為第i個(gè)彈簧的松弛模量;τi為第i個(gè)阻尼器松弛時(shí)間。

    為計(jì)算松弛模量的傅里葉變換,定義無(wú)量綱松弛函數(shù)e(t)為

    (3)

    式中:ei=Ei/E∞。

    本文采用一階Prony級(jí)數(shù)(n=1)模擬人耳軟組織的粘彈性特性。對(duì)穩(wěn)態(tài)動(dòng)力學(xué)分析而言,時(shí)域的松弛模量轉(zhuǎn)換化成頻域復(fù)模量為

    E*(ω)=E′(ω)+iE″(ω)

    (4)

    式中:ω為圓頻率;E′(ω),E″(ω)分別為存儲(chǔ)模量及損耗模量,表達(dá)式為

    (5)

    復(fù)模量E*(ω)表達(dá)動(dòng)態(tài)分析中隨頻率變化的材料特性,能量損耗機(jī)制通過(guò)損耗模量E″(ω)反映到復(fù)模量虛部。

    對(duì)人耳系統(tǒng),由于動(dòng)態(tài)聲壓幅值較小,因此,采用穩(wěn)態(tài)小振幅振動(dòng)分析假設(shè),參數(shù)e1及τ1保持不變。而長(zhǎng)期模量E∞則對(duì)應(yīng)表2的楊氏模量。本文所建有限元模型中,中耳軟組織的5個(gè)部分均采用線性粘彈性本構(gòu)模型,即鼓膜張緊部、鼓膜松弛部、砧錘關(guān)節(jié)、砧鐙關(guān)節(jié)及鐙骨環(huán)韌帶。材料參數(shù)e1,τ1見(jiàn)文獻(xiàn)[21],圓窗材料參數(shù)見(jiàn)文獻(xiàn)[22](表3)。

    2結(jié)果分析

    用有限元分析軟件Abaqus 6.10于外耳道距離鼓膜約1.7 mm處,施加90 dB SPL(0.632 Pa,有效值)的聲壓激勵(lì),分析頻率范圍250~8 000 Hz。計(jì)算包括鼓膜臍部與鐙骨底板位移、鐙骨底板速度傳遞函數(shù)、中耳壓力增益及耳蝸輸入聲阻抗等動(dòng)態(tài)參數(shù)。在此基礎(chǔ)上,于圓窗膜中耳腔一側(cè)施加1 Pa的壓力載荷,計(jì)算耳蝸壓力逆向傳遞函數(shù)。在與相關(guān)實(shí)驗(yàn)測(cè)量數(shù)據(jù)對(duì)比同時(shí)將線彈性+瑞利阻尼及粘彈性計(jì)算結(jié)果進(jìn)行比對(duì)。其中瑞利阻尼系數(shù)取α=0 s-1,β=0.75×10-4s。

    2.1鼓膜臍部與鐙骨底板位移

    模型計(jì)算的鼓膜臍部與鐙骨底板位移響應(yīng)幅值及與10組人體顳骨實(shí)驗(yàn)數(shù)據(jù)平均值[23]比較見(jiàn)圖4。由圖4看出,模型計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)均值一致性較好,鼓膜臍部與鐙骨底板位移均在1 000 Hz處取得最大值,且頻率大于1 000 Hz時(shí)位移幅值隨頻率增大逐漸減小。通過(guò)粘彈性與線彈性比對(duì)知,用粘彈性本構(gòu)模型時(shí)鼓膜臍部與鐙骨底板位移響應(yīng)在250~1 000 Hz低頻段相差不大,隨頻率增大粘彈性計(jì)算結(jié)果逐漸趨近實(shí)驗(yàn)曲線,并高于線彈性計(jì)算結(jié)果。其中鼓膜臍部位移曲線從2 000 Hz處開(kāi)始高于線彈性計(jì)算曲線,而鐙骨底板則始于1 500 Hz頻率處。圖4(b)為模型計(jì)算結(jié)果,有限元模型基于外耳道、中耳及耳蝸的集成,且中耳軟組織采用線彈性本構(gòu)模型。粘彈性位移曲線高于1 000 Hz時(shí)明顯高于文獻(xiàn)[8]結(jié)果,與實(shí)驗(yàn)數(shù)據(jù)更接近。鐙骨底板作為連接中耳與耳蝸的關(guān)鍵部位,其位移響應(yīng)為衡量人耳模型有效性、評(píng)價(jià)植入式中耳助聽(tīng)裝置[24]的重要指標(biāo)。因此,較好匹配鐙骨底板位移實(shí)驗(yàn)曲線,會(huì)為模型進(jìn)一步應(yīng)用提供基礎(chǔ)條件。

    2.2鐙骨底板速度傳遞函數(shù)

    模型計(jì)算的鐙骨底板速度傳遞函數(shù)及文獻(xiàn)[25]的實(shí)驗(yàn)曲線均值及上、下限見(jiàn)圖5,用于表征中耳的聲傳遞特性,定義為鐙骨底板速度與外耳道內(nèi)鼓膜附近聲壓之比。實(shí)驗(yàn)數(shù)據(jù)取自12例人體顳骨樣本。用粘彈性本構(gòu)計(jì)算結(jié)果總體趨勢(shì)與實(shí)驗(yàn)曲線一致,尤其在250~4 500 Hz頻段內(nèi),計(jì)算、實(shí)驗(yàn)曲線一致性較好。較線彈性,粘彈性計(jì)算結(jié)果亦在高于1 500 Hz時(shí)出現(xiàn)較明顯的增幅現(xiàn)象,處于實(shí)驗(yàn)曲線范圍內(nèi)。

    圖4 中耳位移響應(yīng)Fig.4Displacementresponseofthemiddleear圖5 鐙骨底板速度傳遞函數(shù)Fig.5Stapesfootplatevelocitytransferfunction

    2.3中耳壓力增益

    中耳壓力增益定義為前庭階靠近卵圓窗處的壓力值與外耳道內(nèi)鼓膜附近的壓力值之比,并以dB形式表達(dá)。本文前庭階壓力提取點(diǎn)為耳蝸內(nèi)部距卵圓窗0.23 mm,距耳蝸隔膜1.5 mm,與文獻(xiàn)[26]實(shí)測(cè)點(diǎn)0.2 mm及2~3 mm較接近。計(jì)算結(jié)果、相關(guān)實(shí)測(cè)數(shù)據(jù)[27]及其它有限元模型的仿真曲線見(jiàn)圖6。由圖6看出,中耳軟組織采用粘彈性材料特性時(shí)模型計(jì)算的中耳壓力增益曲線與Puria等實(shí)驗(yàn)結(jié)果較接近,尤其在頻率高于2 000 Hz時(shí)。在250~850 Hz低頻段及3000~8 000 Hz高頻段內(nèi),模型計(jì)算結(jié)果均高于仿真曲線,與實(shí)驗(yàn)數(shù)據(jù)更吻合。

    對(duì)比圖6中線彈性與粘彈性結(jié)果曲線可知,二者差異主要集中在2 000~8 000 Hz頻率范圍內(nèi)。在8 000 Hz頻率點(diǎn)處,對(duì)應(yīng)線彈性壓力增益值為1.99 dB,而對(duì)應(yīng)粘彈性增益值為11.88 dB,相差量值達(dá)9.89 dB。結(jié)果表明,中耳軟組織采用粘彈性本構(gòu)關(guān)系亦會(huì)對(duì)耳蝸壓力響應(yīng)產(chǎn)生顯著影響。由模型計(jì)算知,僅考慮圓窗粘彈性特性時(shí)前庭階的壓力響應(yīng)與線彈性基本相同。研究表明,耳蝸前庭階及鼓階壓力差可作為耳蝸的輸入信號(hào),并用于圓窗激勵(lì)式人工中耳的植入性能評(píng)價(jià)[28],因此,圖6中前庭階壓力響應(yīng)的高頻改善會(huì)益于研究。

    2.4耳蝸輸入聲阻抗

    耳蝸輸入聲阻抗為理解聲能量由中耳傳遞到耳蝸的重要參數(shù),定義式[29]為

    (6)

    式中:PSV為前庭階近卵圓窗處壓力(提取點(diǎn)用與中耳壓力增益計(jì)算中相同位置);VS為鐙骨底板速度;AS=3.2 mm2為鐙骨底板面積。

    模型計(jì)算的耳蝸輸入聲阻抗及實(shí)驗(yàn)曲線[30]見(jiàn)圖7。由圖7看出,計(jì)算結(jié)果在總趨勢(shì)上與實(shí)驗(yàn)曲線相似,雖不同實(shí)驗(yàn)數(shù)據(jù)偏差較大,但結(jié)果曲線基本位于各組實(shí)驗(yàn)曲線范圍之內(nèi)。對(duì)應(yīng)線彈性與粘彈性不同材料屬性,聲阻抗最小值均位于1 000 Hz頻率點(diǎn)處,大小為17.89 GΩ。與圖4~圖6計(jì)算結(jié)果不同,在聲阻抗計(jì)算曲線高頻段未出現(xiàn)較明顯的粘彈性高于線彈性現(xiàn)象,且在250~750 Hz頻率范圍內(nèi)線彈性計(jì)算結(jié)果卻高于粘彈性。需進(jìn)一步指出的是,耳蝸輸入聲阻抗需綜合考慮中耳結(jié)構(gòu)與耳蝸流體的動(dòng)態(tài)響應(yīng),而在人耳數(shù)值仿真中不易準(zhǔn)確計(jì)算,人耳軟組織材料屬性變化對(duì)聲阻抗影響。研究表明[31],前庭階截面積也是影響聲阻抗的重要因素,故準(zhǔn)確仿真計(jì)算尚需實(shí)驗(yàn)、建模的深入研究。

    2.5耳蝸壓力逆向傳遞函數(shù)

    為進(jìn)一步分析引入粘彈性材料特性對(duì)人耳動(dòng)力學(xué)響應(yīng)影響,基于所建整耳模型計(jì)算耳蝸壓力逆向傳遞函數(shù)。此時(shí)激勵(lì)壓力加載于圓窗膜中耳腔一側(cè),大小為1 Pa,用于計(jì)算逆向激勵(lì)下耳蝸前庭階及鼓階的流體壓力響應(yīng)。其中,前庭階與鼓階壓力提取點(diǎn)分別為:距卵圓窗、圓窗中心點(diǎn)0.23 mm,距耳蝸隔膜1.5 mm。耳蝸壓力逆向傳遞函數(shù)為

    (7)

    模型計(jì)算耳蝸壓力逆向傳遞函數(shù)與三組實(shí)驗(yàn)結(jié)果[32](單位dB)見(jiàn)圖8。由圖8看出,模型計(jì)算曲線基本處于三組實(shí)驗(yàn)曲線范圍之內(nèi),且波動(dòng)主要在低于1 000 Hz頻率處,在1 000~8 000 Hz頻率內(nèi),結(jié)果曲線較平坦且接近0 dB。由線、粘彈性結(jié)果比對(duì)可知,引入粘彈性材料特性對(duì)耳蝸壓力逆向傳遞函數(shù)影響主要集中在250~1 000 Hz低頻段,最大差異量為250 Hz的3 dB;而頻率高于1 000 Hz時(shí)最大差異量?jī)H8 000 Hz的1.52 dB。此外,對(duì)計(jì)算結(jié)果,逆向激勵(lì)下鼓階壓力值普遍大于前庭階,因此,粘彈性曲線在低頻段下降現(xiàn)象表明,引入粘彈性會(huì)一定程度上提高逆向激勵(lì)的前庭階壓力響應(yīng)。需指出的是,現(xiàn)階段將人工中耳植入圓窗,采用逆向激勵(lì)形式進(jìn)行人耳聽(tīng)力補(bǔ)償,已取得較佳實(shí)驗(yàn)效果。在圓窗激勵(lì)評(píng)價(jià)中,耳蝸壓力逆向傳遞函數(shù)通??捎糜谟?jì)算等效驅(qū)動(dòng)力、驅(qū)動(dòng)位移等重要參數(shù)[32],故該計(jì)算結(jié)果有助于提高圓窗激勵(lì)的數(shù)值評(píng)價(jià)精度。

    圖6 中耳壓力增益Fig.6Middleearpressuregain圖7 耳蝸輸入聲阻抗Fig.7Cochlearinputimpedance圖8 耳蝸壓力逆向傳遞函數(shù)Fig.8Reversepressuretransferfunctionofcochlea

    3討論

    本文研究主要目的在于分析中耳軟組織引入粘彈性材料屬性對(duì)人耳動(dòng)力學(xué)特性影響。圖4、圖5分別從鼓膜臍部與鐙骨底板位移及鐙骨底板速度傳遞函數(shù)三方面,反映粘彈性本構(gòu)關(guān)系能改善中耳結(jié)構(gòu)的高頻響應(yīng),圖6的中耳壓力增益則表明,耳蝸流體壓力響應(yīng)同樣可通過(guò)考慮粘彈性達(dá)到與實(shí)驗(yàn)數(shù)據(jù)的較佳匹配。由此表明,用粘彈性表征人耳系統(tǒng)動(dòng)態(tài)過(guò)程中能量損耗機(jī)制更合適。而模型中采用粘彈性材料屬性的人耳結(jié)構(gòu)數(shù)量較有限,雖中耳韌帶、肌腱對(duì)人耳動(dòng)力學(xué)特性具有重要影響[33],但關(guān)于結(jié)構(gòu)的粘彈性材料實(shí)驗(yàn)仍無(wú)報(bào)道,因此模型中仍采用線彈性本構(gòu)模型。實(shí)驗(yàn)研究[34-35]表明,人體與豚鼠的鼓膜不僅具有頻率相關(guān)的粘彈性特性,在準(zhǔn)靜態(tài)域亦呈較明顯的超彈性特性。本文的人耳動(dòng)力學(xué)建模,未考慮軟組織超彈性材料屬性,在動(dòng)態(tài)分析中人耳系統(tǒng)的平衡位置始終保持不變。因此,綜合考慮人耳軟組織的超彈性、粘彈性特性,分析不同中耳靜壓加載下的人耳動(dòng)力學(xué)特性,將成為新的研究重點(diǎn)。

    圖9 基底膜峰值響應(yīng)位置與頻率關(guān)系 Fig.9 Place of maximum response on the basilar membrane vs. frequency

    本文的整耳有限元模型,耳蝸部分采用簡(jiǎn)化的雙腔導(dǎo)管形式,并未對(duì)蝸管及柯蒂氏器進(jìn)行實(shí)際建模,因此仍為被動(dòng)的耳蝸模型。基底膜沿長(zhǎng)度方向剛度漸變特性主要通過(guò)厚度尺寸、寬度尺寸及楊氏模量變化表達(dá)。目前關(guān)于基底膜厚度、寬度方向尺寸已有較多研究,并被相關(guān)耳蝸模型采用,但由于耳蝸內(nèi)部復(fù)雜的生理?xiàng)l件及基底膜螺旋形結(jié)構(gòu)特點(diǎn),關(guān)于基底膜材料屬性實(shí)驗(yàn)研究仍較少。為此,參考文獻(xiàn)[7]實(shí)驗(yàn)數(shù)據(jù),以基底膜楊氏模量為調(diào)節(jié)參數(shù),通過(guò)匹配實(shí)驗(yàn)曲線確定基底膜最終的模量參數(shù)?;啄し逯滴灰莆恢门c頻率的對(duì)應(yīng)關(guān)系見(jiàn)圖9,表示出基底膜的選頻特性。由圖9看出,沿基底膜長(zhǎng)度方向,共振頻率位置隨頻率增加逐漸由基底膜頂部32.2 mm移到底部7 mm處,并與實(shí)驗(yàn)曲線具有較好的一致性。需指出的是,由于基底膜厚度方向尺寸(2.5~7.5 μm)遠(yuǎn)小于另兩方向尺寸(0.15~34 mm),用實(shí)體單元?jiǎng)澐志W(wǎng)格會(huì)導(dǎo)致過(guò)多單元數(shù)量,導(dǎo)致計(jì)算困難。為此,基底膜采用殼單元進(jìn)行網(wǎng)格劃分,而沿長(zhǎng)度方向變化的厚度尺寸及楊氏模量則通過(guò)定義溫度場(chǎng)變量方式實(shí)現(xiàn)。圖9計(jì)算結(jié)果亦表明設(shè)定的合理性。

    4結(jié)論

    (1)通過(guò)所建含外耳道、中耳及簡(jiǎn)化耳蝸的集成人耳有限元模型,分析中耳軟組織粘彈性材料特性引入對(duì)人耳動(dòng)力學(xué)特性影響。并計(jì)算包括鼓膜臍部與鐙骨底板位移、鐙骨底板的速度傳遞函數(shù)、中耳壓力增益、耳蝸輸入聲阻抗及耳蝸壓力逆向傳遞函數(shù)等動(dòng)態(tài)響應(yīng)參數(shù)。

    (2)模型計(jì)算結(jié)果顯示,較線彈性,粘彈性材料屬性不僅能明顯提高中耳結(jié)構(gòu)的高頻響應(yīng),對(duì)耳蝸流體壓力響應(yīng)亦影響顯著。模型計(jì)算曲線與實(shí)驗(yàn)測(cè)量曲線比較表明,基于粘彈性本構(gòu)關(guān)系的動(dòng)態(tài)響應(yīng)能取得與實(shí)驗(yàn)結(jié)果更佳的匹配效果、用粘彈性表征人耳系統(tǒng)動(dòng)態(tài)能量損耗更合理。該研究有助于建立更貼近實(shí)際生理?xiàng)l件的人耳有限元模型。

    參考文獻(xiàn)

    [1]姚文娟,陳懿強(qiáng),葉志明,等. 耳聽(tīng)力系統(tǒng)生物力學(xué)研究進(jìn)展[J]. 力學(xué)與實(shí)踐, 2013, 35(6): 1-10.

    YAO Wen-juan, CHEN Yi-qiang, YE Zhi-ming, et al. Advance in biomechanics of human ear as hearing system [J]. Mechanics in Engineering, 2013, 35(6): 1-10.

    [2]Nie X, Liu H, Huang X, et al. Finite element model of human ear reconstruction through micro-computer tomography[J]. Acta Oto-Laryngologica, 2011, 131(3):269-276.

    [3]王學(xué)林,胡于進(jìn). 蝸窗激勵(lì)評(píng)價(jià)的有限元計(jì)算模型研究 [J]. 力學(xué)學(xué)報(bào), 2012, 44(3): 622-630.

    WANG Xue-lin, HU Yu-jin. Evaluation of round window stimulation by a FE model of human auditory periphery [J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(3): 622-630.

    [4]Cheng T, Gan R Z. Experimental measurement and modeling analysis on mechanical properties of tensor tympani tendon [J]. Medical Engineering & Physics, 2008, 30(3): 358-366.

    [5]Zhang X, Gan R Z. Dynamic properties of human tympanic membrane based on frequency-temperature superposition [J]. Annals of Biomedical Engineering, 2012, 41(1):205-214.

    [6]Kim N, Homma K, Puria S. Inertial bone conduction: symmetric and anti-symmetric components[J]. Journal of the Association for Research in Otolaryngology, 2011, 12(3): 261-279.

    [7]Greenwood D D. A cochlear frequency-position function for several species-29 years later[J]. Journal of the Acoustical Society of America, 1990, 87(6): 2592-2605.

    [8]Gan R Z, Reeves B P, Wang X L. Modeling of sound transmission from ear canal to cochlea[J]. Annals of Biomedical Engineering, 2007, 35(12): 2180-2195.

    [9]劉后廣,塔娜,饒柱石. 懸浮振子對(duì)中耳聲傳播特性影響的數(shù)值研究 [J]. 力學(xué)學(xué)報(bào), 2010, 42(1): 109-114.

    LIU Hou-guang, TA Na, RAO Zhu-shi. Numerical study on the effect of the floating mass transducer on middle ear sound transmission [J]. Chinese Journal of Theoretical andApplied Mechanics, 2010, 42(1):109-114.

    [10]Stieger C, Rosowski J J, Nakajima H H. Comparison of forward (ear-canal) and reverse (round-window) sound stimulation of the cochlea[J]. Hearing Research, 2013, 301: 105-114.

    [11]Nakajima H H, Dong W, Olson E S, et al. Evaluation of round window stimulation using the floating mass transducer by intracochlear sound pressure measurements in human temporal bones [J]. Otology & Neurotology, 2010, 31(3): 506-511.

    [12]Sun Q, Gan R Z, Chang K H, et al. Computer-integrated finite element modeling of human middle ear[J].Biomechanics and Modeling in Mechanobiology, 2002,2(1): 109-122.

    [13]Kwacz M, Marek P, Borkowski P, et al. A three-dimensional finite element model of round window membrane vibration before and after stapedotomy surgery [J]. Biomechanics and Modeling in Mechanobiology, 2013: 1-19.

    [14]Wang J, Zhao F, Li Y, et al. Effect of anterior tympanomeatal angle blunting on the middle ear transfer function using a finite element ear model[J]. Medical Engineering & Physics, 2011, 33(9): 1136-1146.

    [15]Koike T, Wada H, Kobayashi T. Modeling of the human middle ear using the finite-element method [J]. The Journal of the Acoustical Society of America, 2002, 111(3): 1306-1317.

    [16]Gan R Z, Sun Q, Feng B, et al. Acoustic-structural coupled finite element analysis for sound transmission in human ear-Pressure distributions[J]. Medical Engineering & Physics, 2006, 28(5): 395-404.

    [17]Zhang X, Gan R Z. Finite element modeling of energy absorbance in normal and disordered human ears[J]. Hearing Research, 2013,301:146-155.

    [18]Gan R Z, Yang F, Zhang X, et al. Mechanical properties of stapedial annular ligament [J]. Medical Engineering & Physics, 2011, 33(3): 330-339.

    [19]Machiraju C, Phan A V, Pearsall A, et al. Viscoelastic studies of human subscapularis tendon: Relaxation test and a Wiechert model [J]. Computer Methods and Programs in Biomedicine, 2006, 83(1): 29-33.

    [20]Park S, Schapery R. Methods of interconversion between linear viscoelastic material functions.part I-a numerical method based on prony series [J]. International Journal of Solids and Structures, 1999, 36(11): 1653-1675.

    [21]Zhang X, Gan R Z. A comprehensive model of human ear for analysis of implantable hearing devices [J]. Biomedical Engineering, IEEE Transactions on, 2011, 58(10): 3024-3027.

    [22]Zhang X, Gan R Z. Dynamic properties of human round window membrane in auditory frequencies running head: Dynamic properties of round window membrane[J]. Medical Engineering & Physics, 2013, 35(3): 310-318.

    [23]Gan R Z, Wood M W, Dormer K J. Human middle ear transfer function measured by double laser interferometry system [J]. Otology & Neurotology, 2004, 25(4): 423-435.

    [24]Wang X L, Hu Y J, Wang Z L, et al. Finite element analysis of the coupling between ossicular chain and mass loading for evaluation of implantable hearing device[J]. Hearing Research, 2011, 280(1): 48-57.

    [25]Aibara R, Welsh J T, Puria S, et al. Human middle-ear sound transfer function and cochlear input impedance [J]. Hearing Research, 2001, 152(1): 100-109.

    [26]Nakajima H H, Dong W, Olson E S, et al. Differential intracochlear sound pressure measurements in normal human temporal bones [J]. Journal of the Association for Research in Otolaryngology, 2009, 10(1): 23-36.

    [27]Puria S, Peake W T, Rosowski J J. Sound-pressure measurements in the cochlear vestibule of human-cadaver ears [J]. The Journal of the Acoustical Society of America, 1997, 101(5): 2754-2770.

    [28]王學(xué)林. 蝸窗激勵(lì)與外耳道激勵(lì)產(chǎn)生的耳蝸壓力差的比較分析[J].生物醫(yī)學(xué)工程學(xué)雜志,2012,29(6): 1109-1113.

    WANG Xue-lin. Comparison of differential intracochlear pressures between round window stimulation and ear canal stimulation[J]. Journal of Biomedical Engineering, 2012,29(6):1109-1113.

    [29]Zwislocki J. The role of the external and middle ear in sound transmission[J]. The Nervous System,1975, 3: 45-55.

    [30]Merchant S N, Ravicz M E, Rosowski J J. Acoustic input impedance of the stapes and cochlea in human temporal bones [J]. Hearing Research, 1996, 97(1): 30-45.

    [31]Wang X, Wang L, Zhou J, et al. Finite element modelling of human auditory periphery including a feed-forward amplification of the cochlea[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2012,17(10): 1096-1107.

    [32]Nakajima H H, Merchant S N, Rosowski J J. Performance considerations of prosthetic actuators for round-window stimulation [J]. Hearing Research,2010,263(1):114-119.

    [33]Gentil F, Parente M, Martins P, et al. The influence of muscles activation on the dynamical behaviour of the tympano-ossicular system of the middle ear[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2013, 16(4): 392-402.

    [34]Aernouts J, Dirckx J J. Static versus dynamic gerbil tympanic membrane elasticity: derivation of the complex modulus[J].Biomechanics and Modeling in Mechanobiolo-gy, 2012,11(6): 829-840.

    [35]Aernouts J, Aerts J R, Dirckx J J. Mechanical properties of human tympanic membrane in the quasi-static regime from in situ point indentation measurements [J]. Hearing Research, 2012, 290(1): 45-54.

    猜你喜歡
    粘彈性有限元分析耳蝸
    二維粘彈性棒和板問(wèn)題ADI有限差分法
    耳蝸微音器電位臨床操作要點(diǎn)
    時(shí)變時(shí)滯粘彈性板方程的整體吸引子
    不可壓粘彈性流體的Leray-α-Oldroyd模型整體解的存在性
    自錨式懸索橋鋼箱梁頂推施工階段結(jié)構(gòu)分析
    隨機(jī)振動(dòng)載荷下發(fā)射裝置尾罩疲勞壽命分析
    航空兵器(2016年4期)2016-11-28 21:54:01
    有限元分析帶溝槽平封頭的應(yīng)力集中
    飛機(jī)起落架支撐桿強(qiáng)度有限元分析
    科技視界(2016年18期)2016-11-03 22:31:14
    DR內(nèi)聽(tīng)道像及多層螺旋CT三維重建對(duì)人工耳蝸的效果評(píng)估
    豚鼠耳蝸Hensen細(xì)胞脂滴的性質(zhì)與分布
    国产福利在线免费观看视频| 午夜91福利影院| 亚洲视频免费观看视频| 97精品久久久久久久久久精品| 欧美久久黑人一区二区| 操美女的视频在线观看| 蜜桃在线观看..| 男人舔女人的私密视频| 另类精品久久| 卡戴珊不雅视频在线播放| 欧美xxⅹ黑人| 一个人免费看片子| 亚洲精品中文字幕在线视频| 亚洲精品国产色婷婷电影| 视频在线观看一区二区三区| 亚洲男人天堂网一区| 久久久国产欧美日韩av| 9色porny在线观看| av免费观看日本| 色网站视频免费| 国产精品香港三级国产av潘金莲 | 国产精品偷伦视频观看了| 欧美日本中文国产一区发布| 老熟女久久久| a 毛片基地| 精品少妇久久久久久888优播| 日韩一卡2卡3卡4卡2021年| 国产女主播在线喷水免费视频网站| 日韩,欧美,国产一区二区三区| 纯流量卡能插随身wifi吗| 水蜜桃什么品种好| 最近的中文字幕免费完整| 超碰成人久久| 日韩大片免费观看网站| 1024视频免费在线观看| 欧美日韩一区二区视频在线观看视频在线| 久久久久久久久久久免费av| 成人午夜精彩视频在线观看| 天天躁夜夜躁狠狠躁躁| 免费黄色在线免费观看| 男的添女的下面高潮视频| 女人被躁到高潮嗷嗷叫费观| 国产亚洲最大av| 少妇被粗大猛烈的视频| av在线观看视频网站免费| av.在线天堂| 国产日韩欧美在线精品| 亚洲熟女毛片儿| 最近的中文字幕免费完整| 一级黄片播放器| 国产日韩欧美视频二区| 久久人人爽av亚洲精品天堂| 精品午夜福利在线看| 超色免费av| 亚洲精品视频女| 一级片免费观看大全| 亚洲在久久综合| 国产伦理片在线播放av一区| 男女午夜视频在线观看| videos熟女内射| 女人高潮潮喷娇喘18禁视频| 日韩中文字幕欧美一区二区 | 亚洲欧洲日产国产| 欧美日韩视频精品一区| 91aial.com中文字幕在线观看| 汤姆久久久久久久影院中文字幕| 亚洲成色77777| tube8黄色片| 精品国产超薄肉色丝袜足j| 777久久人妻少妇嫩草av网站| 亚洲自偷自拍图片 自拍| 久久久久久久久免费视频了| 精品一区二区免费观看| 男女之事视频高清在线观看 | 狂野欧美激情性bbbbbb| 大片电影免费在线观看免费| 国产成人欧美在线观看 | 国产无遮挡羞羞视频在线观看| 黑人欧美特级aaaaaa片| 自拍欧美九色日韩亚洲蝌蚪91| 在线观看三级黄色| 亚洲欧洲日产国产| 日韩精品有码人妻一区| 麻豆av在线久日| 精品人妻熟女毛片av久久网站| 伊人久久大香线蕉亚洲五| 亚洲精品第二区| 9191精品国产免费久久| 欧美日韩视频高清一区二区三区二| 国产亚洲欧美精品永久| 免费观看a级毛片全部| 亚洲激情五月婷婷啪啪| 精品午夜福利在线看| 日韩不卡一区二区三区视频在线| 免费人妻精品一区二区三区视频| 精品免费久久久久久久清纯 | 人人妻人人添人人爽欧美一区卜| 国产成人欧美在线观看 | 国产精品久久久久久久久免| 麻豆精品久久久久久蜜桃| 亚洲av综合色区一区| 18禁观看日本| av不卡在线播放| 热re99久久国产66热| 9191精品国产免费久久| 亚洲欧美清纯卡通| 丝袜脚勾引网站| 麻豆精品久久久久久蜜桃| 国产黄色视频一区二区在线观看| 九色亚洲精品在线播放| 另类精品久久| 91成人精品电影| 亚洲欧美成人精品一区二区| 美女午夜性视频免费| 亚洲精品国产区一区二| 秋霞在线观看毛片| 欧美日韩精品网址| 午夜91福利影院| 亚洲精品国产色婷婷电影| 在线观看三级黄色| 美女扒开内裤让男人捅视频| 国产人伦9x9x在线观看| 亚洲国产欧美网| 亚洲在久久综合| 欧美黑人欧美精品刺激| 国产男女超爽视频在线观看| 国产成人一区二区在线| 欧美黄色片欧美黄色片| 国产av精品麻豆| 美女大奶头黄色视频| 男女之事视频高清在线观看 | 精品亚洲成国产av| 嫩草影视91久久| 国产亚洲精品第一综合不卡| 毛片一级片免费看久久久久| 黄色毛片三级朝国网站| 好男人视频免费观看在线| 日本爱情动作片www.在线观看| 天天躁日日躁夜夜躁夜夜| 成年人午夜在线观看视频| 自线自在国产av| 欧美av亚洲av综合av国产av | 国产福利在线免费观看视频| 人妻一区二区av| 男女床上黄色一级片免费看| 久久久国产精品麻豆| 在线观看免费视频网站a站| 日韩 欧美 亚洲 中文字幕| 欧美变态另类bdsm刘玥| 国产在线一区二区三区精| 亚洲av在线观看美女高潮| 免费看不卡的av| 精品人妻熟女毛片av久久网站| 国产成人一区二区在线| 人人妻人人爽人人添夜夜欢视频| 精品人妻熟女毛片av久久网站| 交换朋友夫妻互换小说| 日韩伦理黄色片| 欧美日韩视频精品一区| av有码第一页| 日韩,欧美,国产一区二区三区| 狠狠婷婷综合久久久久久88av| 1024视频免费在线观看| 中文欧美无线码| 中文天堂在线官网| 国产精品国产三级专区第一集| 久久性视频一级片| 99久久精品国产亚洲精品| 国产麻豆69| 久久久久精品久久久久真实原创| 少妇人妻 视频| 国产麻豆69| 国产成人欧美在线观看 | 如何舔出高潮| 欧美成人精品欧美一级黄| 一边摸一边做爽爽视频免费| 一边摸一边抽搐一进一出视频| 亚洲一区中文字幕在线| 国产精品久久久久久精品电影小说| 欧美日韩精品网址| 国产精品国产三级国产专区5o| 国产精品一区二区精品视频观看| 97人妻天天添夜夜摸| 国产精品国产三级国产专区5o| 亚洲av综合色区一区| 国产精品久久久久久久久免| 国产爽快片一区二区三区| 国产日韩欧美亚洲二区| 欧美av亚洲av综合av国产av | 欧美日韩一级在线毛片| 国产精品无大码| 观看av在线不卡| 欧美另类一区| 精品一区在线观看国产| 久久人人爽av亚洲精品天堂| 国产成人a∨麻豆精品| 高清欧美精品videossex| 丝袜在线中文字幕| 秋霞伦理黄片| 国产男女超爽视频在线观看| 亚洲av电影在线进入| 国产有黄有色有爽视频| 男人操女人黄网站| svipshipincom国产片| 日日摸夜夜添夜夜爱| 亚洲精品av麻豆狂野| 日本猛色少妇xxxxx猛交久久| av天堂久久9| 丝袜美足系列| av女优亚洲男人天堂| 日韩欧美精品免费久久| 亚洲精品国产av成人精品| 老司机在亚洲福利影院| av在线app专区| 免费观看av网站的网址| 最新在线观看一区二区三区 | 啦啦啦 在线观看视频| 精品国产一区二区三区久久久樱花| 国产熟女欧美一区二区| 欧美日韩亚洲综合一区二区三区_| 亚洲av电影在线进入| 天堂中文最新版在线下载| 国产 一区精品| 午夜av观看不卡| 国产黄色免费在线视频| 最黄视频免费看| 成人国产av品久久久| videosex国产| 中文欧美无线码| 视频区图区小说| 精品一区二区三区四区五区乱码 | 中文字幕人妻丝袜制服| 在线天堂中文资源库| 国产在线视频一区二区| 色精品久久人妻99蜜桃| 中文字幕色久视频| av国产久精品久网站免费入址| 精品久久久精品久久久| 国产乱来视频区| 久久久国产欧美日韩av| 久久亚洲国产成人精品v| 亚洲av欧美aⅴ国产| 99热网站在线观看| 欧美久久黑人一区二区| 精品卡一卡二卡四卡免费| 最近的中文字幕免费完整| 亚洲精品一区蜜桃| 亚洲精品美女久久久久99蜜臀 | 国产一级毛片在线| www.自偷自拍.com| 亚洲第一av免费看| 国产精品.久久久| 丝袜人妻中文字幕| 国产熟女午夜一区二区三区| 又粗又硬又长又爽又黄的视频| 久久 成人 亚洲| 黄网站色视频无遮挡免费观看| 国产一区有黄有色的免费视频| 国产不卡av网站在线观看| 女性生殖器流出的白浆| 欧美xxⅹ黑人| 国产 一区精品| 777米奇影视久久| av在线播放精品| av女优亚洲男人天堂| 国产欧美日韩一区二区三区在线| 国产免费又黄又爽又色| 久久亚洲国产成人精品v| 精品国产露脸久久av麻豆| 亚洲精品视频女| e午夜精品久久久久久久| 日韩av免费高清视频| 丁香六月欧美| 色婷婷久久久亚洲欧美| 狠狠精品人妻久久久久久综合| 精品国产一区二区三区四区第35| 精品人妻熟女毛片av久久网站| 国产成人精品无人区| 亚洲熟女精品中文字幕| 欧美日韩亚洲国产一区二区在线观看 | 纯流量卡能插随身wifi吗| 两个人免费观看高清视频| 欧美黑人欧美精品刺激| 成年美女黄网站色视频大全免费| 国产 一区精品| 9色porny在线观看| 国产精品国产三级国产专区5o| 热99久久久久精品小说推荐| 老鸭窝网址在线观看| 亚洲美女视频黄频| 丰满乱子伦码专区| 在线 av 中文字幕| 丁香六月天网| 免费观看人在逋| 精品国产一区二区久久| 成人国产av品久久久| 中文天堂在线官网| 大码成人一级视频| 少妇猛男粗大的猛烈进出视频| netflix在线观看网站| 我要看黄色一级片免费的| 亚洲国产精品国产精品| 99国产综合亚洲精品| 色吧在线观看| 国产在视频线精品| 乱人伦中国视频| 美女中出高潮动态图| 亚洲欧美一区二区三区黑人| 国产男人的电影天堂91| 卡戴珊不雅视频在线播放| 国产精品久久久久久精品古装| 亚洲精品美女久久久久99蜜臀 | 欧美激情极品国产一区二区三区| 男人操女人黄网站| 我的亚洲天堂| 国产精品熟女久久久久浪| 99久久人妻综合| 99精国产麻豆久久婷婷| 中文欧美无线码| 一边亲一边摸免费视频| 不卡av一区二区三区| 午夜影院在线不卡| 亚洲欧美精品自产自拍| 一区二区日韩欧美中文字幕| 高清视频免费观看一区二区| 另类亚洲欧美激情| 大香蕉久久网| 99热国产这里只有精品6| 高清不卡的av网站| 久热爱精品视频在线9| 欧美在线黄色| av.在线天堂| 操美女的视频在线观看| 亚洲婷婷狠狠爱综合网| 999久久久国产精品视频| 亚洲人成77777在线视频| 熟女少妇亚洲综合色aaa.| 久久精品国产综合久久久| 亚洲少妇的诱惑av| 免费女性裸体啪啪无遮挡网站| 亚洲免费av在线视频| 久久精品久久精品一区二区三区| 国产一级毛片在线| 欧美日本中文国产一区发布| 美女主播在线视频| 午夜福利视频精品| 日韩制服丝袜自拍偷拍| 国产成人精品在线电影| 久久久久久久久久久久大奶| 精品久久久久久电影网| 精品国产国语对白av| 日本欧美视频一区| 男女边摸边吃奶| 观看美女的网站| 如何舔出高潮| 免费在线观看完整版高清| 久久性视频一级片| 人人妻人人添人人爽欧美一区卜| 亚洲欧美激情在线| av有码第一页| 国产精品亚洲av一区麻豆 | 乱人伦中国视频| 五月开心婷婷网| 人人妻人人澡人人看| 日韩成人av中文字幕在线观看| 精品国产超薄肉色丝袜足j| 黄片小视频在线播放| 国产精品国产三级专区第一集| 精品一区在线观看国产| 免费黄色在线免费观看| 80岁老熟妇乱子伦牲交| 亚洲精品自拍成人| 久久ye,这里只有精品| 欧美日韩精品网址| 欧美变态另类bdsm刘玥| 1024视频免费在线观看| 在线观看国产h片| 日韩一区二区三区影片| 欧美黄色片欧美黄色片| 久久性视频一级片| 无限看片的www在线观看| 午夜免费观看性视频| 国产黄频视频在线观看| 中文字幕最新亚洲高清| 激情视频va一区二区三区| 久久亚洲国产成人精品v| 十八禁人妻一区二区| 久久99热这里只频精品6学生| 丝袜喷水一区| 99国产综合亚洲精品| 国产亚洲最大av| 国产日韩欧美亚洲二区| 中文天堂在线官网| 中文精品一卡2卡3卡4更新| 国产精品嫩草影院av在线观看| 男人添女人高潮全过程视频| 国产黄频视频在线观看| 热re99久久国产66热| 国产精品麻豆人妻色哟哟久久| 国产av码专区亚洲av| 亚洲一级一片aⅴ在线观看| 亚洲综合精品二区| 欧美在线黄色| 日韩av在线免费看完整版不卡| 国产1区2区3区精品| 国产日韩欧美在线精品| 丝袜美腿诱惑在线| 精品一区二区免费观看| 99九九在线精品视频| 女人被躁到高潮嗷嗷叫费观| 亚洲av日韩精品久久久久久密 | 久久鲁丝午夜福利片| 久久久久久久久久久久大奶| 中国国产av一级| 亚洲少妇的诱惑av| 这个男人来自地球电影免费观看 | 无遮挡黄片免费观看| 视频在线观看一区二区三区| 高清欧美精品videossex| a级片在线免费高清观看视频| 母亲3免费完整高清在线观看| 亚洲欧洲精品一区二区精品久久久 | 777久久人妻少妇嫩草av网站| 亚洲美女搞黄在线观看| 日韩熟女老妇一区二区性免费视频| 亚洲成国产人片在线观看| 亚洲欧美激情在线| 大香蕉久久网| 中文字幕色久视频| 国产乱来视频区| 亚洲成人av在线免费| 亚洲欧洲日产国产| 三上悠亚av全集在线观看| 久久人妻熟女aⅴ| 国产成人啪精品午夜网站| 99久久综合免费| 波多野结衣av一区二区av| 1024视频免费在线观看| 久久 成人 亚洲| 亚洲精品久久成人aⅴ小说| av福利片在线| 人人妻人人添人人爽欧美一区卜| 蜜桃国产av成人99| 老司机影院毛片| 国产成人啪精品午夜网站| 天堂俺去俺来也www色官网| 亚洲第一区二区三区不卡| 激情五月婷婷亚洲| 亚洲国产精品国产精品| 色精品久久人妻99蜜桃| 亚洲精品久久午夜乱码| 丁香六月天网| 伊人亚洲综合成人网| 日日爽夜夜爽网站| 国产av一区二区精品久久| 国产精品 国内视频| 成人18禁高潮啪啪吃奶动态图| 久热爱精品视频在线9| 欧美人与善性xxx| 日本av免费视频播放| 色播在线永久视频| 亚洲成人国产一区在线观看 | 制服丝袜香蕉在线| 精品少妇久久久久久888优播| 两个人看的免费小视频| 99久久精品国产亚洲精品| 国产精品嫩草影院av在线观看| 国产精品久久久人人做人人爽| 国产精品无大码| 亚洲国产欧美网| 一级爰片在线观看| 久久天躁狠狠躁夜夜2o2o | 亚洲av中文av极速乱| 成人影院久久| 如何舔出高潮| 在线天堂最新版资源| 在线观看国产h片| 久久精品国产a三级三级三级| 男女免费视频国产| 亚洲色图综合在线观看| 亚洲欧美色中文字幕在线| 婷婷色麻豆天堂久久| 日日爽夜夜爽网站| 亚洲美女搞黄在线观看| 精品国产一区二区三区久久久樱花| 国产一区亚洲一区在线观看| 欧美日韩精品网址| 亚洲中文av在线| svipshipincom国产片| 成年av动漫网址| 男的添女的下面高潮视频| 最近中文字幕高清免费大全6| 日韩视频在线欧美| 国产黄频视频在线观看| 国产成人a∨麻豆精品| 免费少妇av软件| 欧美日韩视频精品一区| 一级毛片我不卡| 久久av网站| 日日撸夜夜添| 激情视频va一区二区三区| 最近中文字幕2019免费版| 国产精品人妻久久久影院| 少妇 在线观看| 久久国产精品大桥未久av| 成人国语在线视频| 亚洲精品一区蜜桃| 操出白浆在线播放| 亚洲精品第二区| 欧美日韩一区二区视频在线观看视频在线| 综合色丁香网| 欧美 日韩 精品 国产| 老鸭窝网址在线观看| 老司机深夜福利视频在线观看 | 黑人欧美特级aaaaaa片| 久久精品熟女亚洲av麻豆精品| 精品福利永久在线观看| 久久精品久久精品一区二区三区| 蜜桃在线观看..| 国产精品欧美亚洲77777| 精品人妻一区二区三区麻豆| 制服诱惑二区| 2021少妇久久久久久久久久久| 老司机影院毛片| 少妇猛男粗大的猛烈进出视频| 欧美黄色片欧美黄色片| 一级毛片电影观看| 91成人精品电影| 婷婷色麻豆天堂久久| 在现免费观看毛片| 国产伦理片在线播放av一区| 国产精品免费大片| 无遮挡黄片免费观看| 99热网站在线观看| 在线亚洲精品国产二区图片欧美| 欧美精品人与动牲交sv欧美| 精品久久久精品久久久| 老汉色av国产亚洲站长工具| 亚洲情色 制服丝袜| 久久天堂一区二区三区四区| 日本猛色少妇xxxxx猛交久久| 美女午夜性视频免费| 成人国语在线视频| 乱人伦中国视频| 欧美人与性动交α欧美精品济南到| 在线免费观看不下载黄p国产| 国产av码专区亚洲av| 免费观看性生交大片5| kizo精华| 久久久久精品国产欧美久久久 | 国产有黄有色有爽视频| 91精品国产国语对白视频| 亚洲在久久综合| 国产99久久九九免费精品| kizo精华| 国产精品一区二区在线不卡| 午夜免费男女啪啪视频观看| 久久久精品94久久精品| 亚洲国产看品久久| 日本猛色少妇xxxxx猛交久久| 国产一区有黄有色的免费视频| 桃花免费在线播放| 男女午夜视频在线观看| 日韩一卡2卡3卡4卡2021年| 国产深夜福利视频在线观看| 中文字幕人妻熟女乱码| 999久久久国产精品视频| 在线观看免费午夜福利视频| 国产精品久久久久久久久免| 日韩大码丰满熟妇| 最近手机中文字幕大全| 国产av精品麻豆| 王馨瑶露胸无遮挡在线观看| 看免费成人av毛片| 亚洲色图综合在线观看| 啦啦啦 在线观看视频| 啦啦啦中文免费视频观看日本| 成年人免费黄色播放视频| 天堂8中文在线网| 亚洲,欧美精品.| 亚洲情色 制服丝袜| 欧美激情 高清一区二区三区| 色婷婷av一区二区三区视频| 嫩草影视91久久| 国产精品久久久久久精品古装| 哪个播放器可以免费观看大片| 超碰97精品在线观看| 香蕉国产在线看| 国产亚洲av高清不卡| 国产黄色免费在线视频| 亚洲欧美日韩另类电影网站| 亚洲av中文av极速乱| 女性被躁到高潮视频| 亚洲第一区二区三区不卡| 欧美人与性动交α欧美软件| 黄频高清免费视频| 国产精品熟女久久久久浪| 亚洲七黄色美女视频| 五月天丁香电影| a级毛片在线看网站| 亚洲精品成人av观看孕妇| 亚洲精品久久久久久婷婷小说| 韩国av在线不卡| 成人黄色视频免费在线看| 久久久久久人妻| 日本vs欧美在线观看视频| 少妇精品久久久久久久| 多毛熟女@视频| 精品人妻熟女毛片av久久网站| netflix在线观看网站| 国产成人一区二区在线| 国精品久久久久久国模美| 99香蕉大伊视频| 侵犯人妻中文字幕一二三四区| 18在线观看网站|