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

    基于完全耦合算法的繞水翼流固耦合特性研究

    2017-08-01 00:02:51王國玉曹樹良
    船舶力學(xué) 2017年7期
    關(guān)鍵詞:水翼慣性阻尼

    吳 欽,黃 彪,王國玉,曹樹良

    (1.清華大學(xué) 熱能工程系,北京 100084;2.北京理工大學(xué) 機(jī)械與車輛學(xué)院,北京 100081)

    基于完全耦合算法的繞水翼流固耦合特性研究

    吳 欽1,黃 彪2,王國玉2,曹樹良1

    (1.清華大學(xué) 熱能工程系,北京 100084;2.北京理工大學(xué) 機(jī)械與車輛學(xué)院,北京 100081)

    基于完全耦合算法對(duì)繞二維NACA0009水翼流固耦合特性進(jìn)行了數(shù)值模擬研究。采用Theodorsen模型和Munch模型對(duì)剛性和彈性水翼的水彈性響應(yīng)進(jìn)行了數(shù)值計(jì)算,分析了流體與結(jié)構(gòu)的相互作用關(guān)系,研究了影響結(jié)構(gòu)水彈性響應(yīng)和流固耦合特性的因素。研究結(jié)果表明:考慮了流體黏性的Munch模型與基于勢(shì)流理論的Theodorsen模型對(duì)氣動(dòng)彈性響應(yīng)的數(shù)值計(jì)算結(jié)果基本一致,而Theodorsen模型由于沒有考慮流體黏性在一定程度上低估了結(jié)構(gòu)的水彈性響應(yīng)。結(jié)構(gòu)的慣性、阻尼和剛度力矩與流體的相應(yīng)附加載荷均處于同一數(shù)量級(jí),故流體與結(jié)構(gòu)的相互作用不可忽略,尤其對(duì)于彈性水翼,流體的慣性、附加阻尼作用增大,流固耦合算法的數(shù)值穩(wěn)定性對(duì)流固耦合特性的計(jì)算結(jié)果影響將更大。外部激勵(lì)頻率為非共振頻率時(shí),結(jié)構(gòu)的剛度作用是影響水彈性響應(yīng)的主要因素,外部激勵(lì)頻率為共振頻率時(shí),流體的附加阻尼和附加剛度作用減弱,除結(jié)構(gòu)的剛度作用外,流體與結(jié)構(gòu)的慣性作用對(duì)水彈性響應(yīng)和流固耦合特性的影響也較大。

    完全耦合算法;水彈性;流固耦合

    0 引 言

    流固耦合振動(dòng)現(xiàn)象在工程和自然界中廣泛存在,對(duì)這種現(xiàn)象的研究在機(jī)械、航空、航天、海洋、建筑和生物等領(lǐng)域都具有十分重要的意義。人們對(duì)于流固耦合振動(dòng)現(xiàn)象的早期認(rèn)識(shí)源于機(jī)翼及葉片的氣動(dòng)彈性問題[1-4]。在與水動(dòng)力學(xué)相關(guān)的水力機(jī)械與船舶工程領(lǐng)域中,復(fù)雜水力環(huán)境下的水彈性響應(yīng)會(huì)導(dǎo)致系統(tǒng)動(dòng)力特性明顯改變,并與流速、壓力、湍動(dòng)能、渦湍粘性等流動(dòng)參數(shù)的變化緊密關(guān)聯(lián),隨流態(tài)變化而異,這顯著加劇了水動(dòng)彈性問題流固耦合研究的理論難度[5-6]。隨著對(duì)水力機(jī)械的安全穩(wěn)定性要求日益提高,國內(nèi)外學(xué)者對(duì)繞水翼流固耦合特性及其影響因素進(jìn)行了大量的研究工作。Gowing[7]針對(duì)兩種不同水彈性自適應(yīng)復(fù)合材料水翼展開實(shí)驗(yàn)研究,結(jié)果表明葉片變形可以有效減小葉片攻角,從而減小葉尖載荷,推延葉尖初生空化。Young[8]建立了一種低階基于勢(shì)流理論的三維邊界元方法用來求解空泡邊界和壓力波動(dòng),結(jié)合有限元軟件確定葉片的動(dòng)態(tài)響應(yīng),建立了邊界元與有限元的流固耦合計(jì)算方法。Ausoni[9]應(yīng)用數(shù)值計(jì)算與實(shí)驗(yàn)相結(jié)合的方法研究了繞水翼空化流動(dòng)和流固耦合特性對(duì)流場(chǎng)渦結(jié)構(gòu)產(chǎn)生機(jī)理的影響,結(jié)果表明,在翼型非共振條件下,初生空化數(shù)與雷諾數(shù)的平方根成正比,而在翼型共振條件下,渦脫落頻率與結(jié)構(gòu)的特征頻率一致,且渦結(jié)構(gòu)的空間相干性具有準(zhǔn)二維性質(zhì)。

    流固耦合問題的重要特征是流體和固體結(jié)構(gòu)之間的相互作用,包括固體結(jié)構(gòu)在流體載荷作用下產(chǎn)生的變形或運(yùn)動(dòng)、固體的變形和運(yùn)動(dòng)對(duì)流場(chǎng)的影響,因此,流固耦合問題的求解需要同時(shí)考慮流場(chǎng)和結(jié)構(gòu)場(chǎng)的求解及其耦合。耦合場(chǎng)的求解算法一般包括完全耦合法(或同步求解法)和分步求解法[10-14]。完全耦合法對(duì)流體和結(jié)構(gòu)建立統(tǒng)一的耦合方程,在一個(gè)時(shí)間步內(nèi)對(duì)流體域和固體域中所有的未知量進(jìn)行求解;分步求解法分別對(duì)流場(chǎng)和結(jié)構(gòu)場(chǎng)選擇合適的數(shù)值算法進(jìn)行獨(dú)立求解,流固耦合界面數(shù)據(jù)通過反復(fù)迭代求解并在界面間反復(fù)傳遞直至獲得收斂解。雖然完全耦合法需要自編程序?qū)崿F(xiàn),但由于其建立的流場(chǎng)與結(jié)構(gòu)場(chǎng)強(qiáng)烈的相關(guān)性而作為解決流固耦合問題的有效方法被廣泛應(yīng)用。Ryzhakov等[15]利用完全耦合法基于拉格朗日體系建立了流場(chǎng)和結(jié)構(gòu)場(chǎng)控制方程,對(duì)水流沖擊彈性平板和注水氣球彈性變形等算例進(jìn)行了分析研究,結(jié)果表明這種方法能精確求解出流體和結(jié)構(gòu)的響應(yīng),且數(shù)值收斂性較好。Michler[16]針對(duì)活塞與流體相互作用的一維模型問題,從算法穩(wěn)定性、計(jì)算精度和計(jì)算效率等方面對(duì)比了不同流固耦合算法在流固耦合問題求解中的應(yīng)用情況,結(jié)果表明完全耦合法無條件穩(wěn)定且計(jì)算精確度相對(duì)較高。

    目前,國內(nèi)外在水動(dòng)彈性力學(xué)的流固耦合數(shù)值計(jì)算研究中多基于理想有勢(shì)流動(dòng),考慮振蕩效應(yīng)、流體流動(dòng)等因素對(duì)物體振動(dòng)特性的影響較少,確定動(dòng)態(tài)參振質(zhì)量、動(dòng)態(tài)剛度和動(dòng)態(tài)阻尼的理論方法尚不完善。本文基于完全耦合算法,利用自編程軟件對(duì)二維NACA0009水翼模型在流固耦合作用下的水彈性響應(yīng)進(jìn)行了數(shù)值計(jì)算,通過對(duì)比不同流體介質(zhì)和不同結(jié)構(gòu)參數(shù)下的水彈性響應(yīng),分析了流體與結(jié)構(gòu)的相互作用關(guān)系,并基于計(jì)算所得彈性系統(tǒng)產(chǎn)生的附加質(zhì)量、附加阻尼和附加剛度,對(duì)繞水翼流固耦合特性進(jìn)行了研究。

    1 控制方程及其求解

    1.1 基本控制方程

    利用有限元對(duì)流場(chǎng)中的連續(xù)質(zhì)量結(jié)構(gòu)進(jìn)行瞬態(tài)分析,結(jié)構(gòu)動(dòng)力學(xué)方程定義為:

    其中:[Ms]、[Cs]和 [Ks]分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣,}分別為結(jié)構(gòu)的位移、速度和加速度,F(xiàn)EX是結(jié)構(gòu)的外部激振力,F(xiàn)HE是流體對(duì)結(jié)構(gòu)的附加作用力。

    計(jì)算采用單自由度的二維NACA0009水翼模型,考慮翼型僅在外部激勵(lì)作用下繞彈性軸發(fā)生振蕩運(yùn)動(dòng),忽略翼型弦長(zhǎng)方向和展向的變形。圖1給出了計(jì)算模型示意圖,圖中U為無窮遠(yuǎn)來流速度,b為半弦長(zhǎng),弦長(zhǎng)c=2b=0.1m,彈性軸位于翼型1/2弦長(zhǎng)處,與結(jié)構(gòu)的重心、壓心重合,θ代表翼型的扭轉(zhuǎn)變形角度。因此,(1)式可簡(jiǎn)化為:

    其中:Iθ為結(jié)構(gòu)關(guān)于彈性軸定義的轉(zhuǎn)動(dòng)慣量,Cθ為結(jié)構(gòu)扭轉(zhuǎn)變形的阻尼系數(shù),Kθ為結(jié)構(gòu)扭轉(zhuǎn)剛度系數(shù),θ、θ˙和θ¨分別為結(jié)構(gòu)的扭轉(zhuǎn)變形角度、角速度和角加速度,M0sinωe( )t是幅值為M0、頻率為ωe的外部激勵(lì)力矩,Mfluid是流體對(duì)結(jié)構(gòu)的附加作用力矩。本文分別采用Theodorsen模型[17]和Munch模型[18]對(duì)流體的附加作用力矩Mfluid進(jìn)行數(shù)值計(jì)算,分別記為

    圖1 單自由度二維NACA0009水翼模型Fig.1 NACA0009 hydrofoil with one degree offreedom

    Theodorsen模型[17]基于不可壓縮、無黏流動(dòng),流體對(duì)結(jié)構(gòu)的附加作用力矩MTfluid表述為:

    其中:MfT、CfT和KfT分別為根據(jù)Theodorsen模型確定的與結(jié)構(gòu)變形加速度、速度和位移有關(guān)的附加水動(dòng)力載荷,即結(jié)構(gòu)的附加質(zhì)量、附加阻尼和附加剛度,ρf為流體密度,k=ωeb/U為折合頻率,Theodorsen函數(shù)C()k是折合頻率的復(fù)函數(shù),表述為:

    C(k)是由于考慮了自由渦的作用而引起的修正項(xiàng),H21(k),H20(k)分別為第二類一階和零階Hankel函數(shù)[19-20]。

    為了考慮流體黏性對(duì)流體作用力的影響,Munch[18]基于繞NACA0009水翼不可壓縮、湍流黏性流動(dòng)的大量數(shù)值計(jì)算與實(shí)驗(yàn)結(jié)果進(jìn)行曲線擬合,提出Munch非定常動(dòng)力模型,將流體對(duì)結(jié)構(gòu)的附加作用力矩定義為:

    1.2 流固耦合算法

    本文采用完全耦合算法(Fully Coupled,以下簡(jiǎn)稱FC)對(duì)流固耦合場(chǎng)進(jìn)行數(shù)值求解。FC算法的優(yōu)點(diǎn)是計(jì)算無條件收斂,計(jì)算穩(wěn)定性和精確度較高,由于同時(shí)求解結(jié)構(gòu)場(chǎng)和流場(chǎng)方程組,計(jì)算效率較高,適用于流體和結(jié)構(gòu)之間具有強(qiáng)烈非線性特性的耦合問題求解。

    圖2給出了FC算法流程圖。FC算法對(duì)流體和結(jié)構(gòu)建立統(tǒng)一的耦合方程,同時(shí)對(duì)流場(chǎng)和結(jié)構(gòu)場(chǎng)進(jìn)行數(shù)值求解計(jì)算,根據(jù)解得的結(jié)構(gòu)變形進(jìn)行流場(chǎng)和結(jié)構(gòu)場(chǎng)的網(wǎng)格更新,對(duì)下一離散時(shí)間步重復(fù)上述求解過程。

    將(2)式在時(shí)域內(nèi)進(jìn)行離散可得:

    基于Theodorsen模型和Munch模型,將(3)式和(5)式分別代入(6)式,其中分別定義為因此(6)式可分別表述為:

    圖2 完全耦合算法流程圖Fig.2 Flowchart of the FC algorithm

    其中:下標(biāo)n表示時(shí)間增量,采用二階Crank-Nicholson方法[21]對(duì)方程進(jìn)行離散求解。

    2 結(jié)果與討論

    本文針對(duì)NACA0009水翼模型,編程實(shí)現(xiàn)流固耦合控制方程的求解,對(duì)繞翼型流固耦合特性進(jìn)行分析。表1給出了NACA0009剛/彈性水翼模型的物質(zhì)屬性,其中:ρs為結(jié)構(gòu)密度,ξs為結(jié)構(gòu)阻尼系數(shù),ωθ為結(jié)構(gòu)在真空中的特征頻率,ωn為結(jié)構(gòu)在水中的特征頻率。采用的計(jì)算工況為來流速度5 m/s,雷諾數(shù)5×105。

    表1 剛/彈性水翼模型物質(zhì)屬性Tab.1 Model parameters for the case of rigid and flexible hydrofoils subject to forced pitching motion

    圖3分別給出了在不同激勵(lì)條件和流體介質(zhì)中,數(shù)值計(jì)算的剛/彈性水翼俯仰角度隨時(shí)間的變化情況,當(dāng)ωe/ωn=0.01時(shí),表征外部激勵(lì)頻率遠(yuǎn)小于共振頻率,當(dāng)ωe/ωn=1時(shí)外部激勵(lì)等于共振頻率。

    對(duì)比圖3(a)、(b)與圖3(c)、(d)可知,當(dāng)流體介質(zhì)為空氣時(shí),采用Theodorsen模型和Munch模型數(shù)值計(jì)算得到的水翼俯仰角度隨時(shí)間的變化情況基本一致;當(dāng)流體介質(zhì)為水時(shí),采用Theodorsen模型和Munch模型對(duì)翼型水彈性響應(yīng)的數(shù)值計(jì)算結(jié)果出現(xiàn)差異,且共振激勵(lì)頻率下不同模型的數(shù)值計(jì)算結(jié)果差異比非共振激勵(lì)頻率下更為明顯。造成這種差異的原因在于,空氣黏性較小,因此,考慮了流體黏性的Munch模型與基于勢(shì)流理論的Theodorsen模型的數(shù)值計(jì)算結(jié)果基本一致,而水的黏性不可忽略,當(dāng)流體介質(zhì)為水時(shí),Theodorsen模型由于不考慮流體黏性而在一定程度上低估了結(jié)構(gòu)的彈性變形量,從而導(dǎo)致Theodorsen模型和Munch模型的數(shù)值計(jì)算結(jié)果差異較大。

    圖3 在不同激勵(lì)條件和流體介質(zhì)中,數(shù)值計(jì)算得到的剛性與彈性水翼俯仰角度隨時(shí)間的變化Fig.3 Comparison of the predicted pitching motion obtained using Theodorsen model and Munch model for the steel/flexible hydrofoil in air/water atωe/ωn=0.01 andωe/ωn=1

    對(duì)比圖3(a)、(c)與圖3(b)、(d)可以看出,當(dāng)結(jié)構(gòu)所受外部激勵(lì)的頻率為非共振頻率時(shí),剛性和彈性水翼的俯仰角度隨時(shí)間的變化情況基本一致,均隨外部激勵(lì)發(fā)生周期性變化;當(dāng)結(jié)構(gòu)所受外部激勵(lì)的頻率為共振頻率時(shí),振幅明顯大于非共振激勵(lì)頻率下的振幅,且剛性和彈性水翼發(fā)生共振時(shí)的振幅和頻率相差較大,水翼在水中的彈性變形量呈明顯發(fā)散趨勢(shì)。這反映了當(dāng)結(jié)構(gòu)所受外部激勵(lì)的頻率接近系統(tǒng)固有頻率時(shí),系統(tǒng)振幅顯著增大的共振現(xiàn)象,同時(shí)結(jié)構(gòu)的材料屬性對(duì)水彈性響應(yīng)的振幅和頻率影響較大。

    圖4給出了分別采用Theodorsen模型和Munch模型數(shù)值計(jì)算得到的剛/彈性水翼附加慣性力矩、阻尼力矩和剛度力矩隨時(shí)間的演變情況。在本小節(jié)的研究中,設(shè)定流體介質(zhì)為水,外部激勵(lì)頻率為非共振頻率,即ωe/ωn=0.01。表2給出了采用不同模型數(shù)值計(jì)算得到的流固慣性力矩、阻尼力矩和剛度矩系數(shù)之比。

    由圖4(a)-(b)可知,對(duì)于相同結(jié)構(gòu),基于Theodorsen模型和Munch模型計(jì)算的慣性力矩隨時(shí)間的演變規(guī)律相同,這是因?yàn)?,由?可知,基于Theodorsen模型和Munch模型計(jì)算所得附加慣性力矩系數(shù)(附加質(zhì)量)相等,即結(jié)合(3)式和(5)式可知,Theodorsen模型與Munch模型對(duì)附加慣性力矩系數(shù)項(xiàng)的定義相同,附加慣性力矩系數(shù)僅與流體密度和結(jié)構(gòu)幾何參數(shù)有關(guān),與折合頻率、雷諾數(shù)等無關(guān)。同時(shí),附加慣性力矩與結(jié)構(gòu)慣性力矩隨時(shí)間的演變呈同相分布,說明在水彈性響應(yīng)的影響下,由流體與結(jié)構(gòu)組成的彈性系統(tǒng)內(nèi)部慣性作用加強(qiáng)。此外,由圖4(a)-(b)不難發(fā)現(xiàn),初始時(shí)刻,剛性水翼的附加慣性力矩小于結(jié)構(gòu)的慣性力矩,而彈性水翼的附加慣性力矩大于結(jié)構(gòu)的慣性力矩。結(jié)合表2可知,剛性和彈性水翼的流固慣性力矩比值分別為0.245和1.718,剛性水翼的密度較大,流固密度比較小,因此流體對(duì)剛性結(jié)構(gòu)的慣性作用比流體對(duì)彈性結(jié)構(gòu)的慣性作用小。

    圖4 剛/彈性水翼的結(jié)構(gòu)慣性、阻尼和剛度力矩(黑色曲線)及流體附加慣性、附加阻尼和附加剛度力矩(紅色曲線為基于Munch模型的數(shù)值計(jì)算結(jié)果,藍(lán)色曲線為基于Theodorsen模型的數(shù)值計(jì)算結(jié)果)隨時(shí)間的演變情況(ωe/ωn=0.01)Fig.4 Comparison of the predicted fluid to solid inertial moments,damping moments and stiffness moments for the steel/flexible hydrofoil in water atωe/ωn=0.01

    由圖4(c)-(d)可知,采用Theodorsen模型數(shù)值計(jì)算得到的附加阻尼力矩小于采用Munch模型計(jì)算得到的附加阻尼力矩。結(jié)合表2可知,采用不同模型數(shù)值計(jì)算得到的剛性水翼流固阻尼力矩比值分別為0.186(Theodorsen模型)和8.914(Munch模型),而對(duì)于彈性水翼,數(shù)值計(jì)算得到的流固阻尼力矩比值則分別為0.905(Theodorsen模型)和18.429(Munch模型)。這是因?yàn)門heodorsen模型基于無黏勢(shì)流假設(shè),低估了流體的阻尼作用力。此外,采用Munch模型計(jì)算的附加阻尼力矩明顯大于結(jié)構(gòu)的阻尼力矩,兩者同相分布,說明水彈性響應(yīng)進(jìn)一步增強(qiáng)了彈性系統(tǒng)的阻尼作用,其中流體附加阻尼作用強(qiáng)于結(jié)構(gòu)阻尼作用。

    由圖4(e)-(f)可知,采用Theodorsen模型和Munch模型計(jì)算的剛度矩隨時(shí)間的演變情況基本一致,結(jié)合表2可知,采用不同模型計(jì)算得到的流固剛度矩比值也很接近。這是由于結(jié)構(gòu)的剛度系數(shù)較大,流固剛度矩比值較小,因此,流體的附加剛度作用對(duì)水彈性響應(yīng)影響較小。與慣性力矩和阻尼力矩不同的是,流體附加剛度矩與結(jié)構(gòu)剛度矩呈反相分布,這說明水彈性響應(yīng)在一定程度上削弱了彈性系統(tǒng)的剛度作用。

    綜合圖4(a)-(f)可以看出,當(dāng)外部激勵(lì)頻率為非共振頻率時(shí),結(jié)構(gòu)的剛度力矩遠(yuǎn)大于結(jié)構(gòu)的慣性力矩、阻尼力矩以及流體的相應(yīng)附加載荷,因此結(jié)構(gòu)的剛度作用是影響結(jié)構(gòu)水彈性響應(yīng)的主要因素。同時(shí),結(jié)構(gòu)的慣性、阻尼和剛度力矩與流體的相應(yīng)附加載荷均處于同一數(shù)量級(jí),因此流體與結(jié)構(gòu)的相互作用不可忽略,尤其對(duì)于彈性水翼,由于流體的慣性、附加阻尼作用增大,流固耦合算法的數(shù)值穩(wěn)定性將對(duì)彈性水翼流固耦合特性的數(shù)值計(jì)算結(jié)果產(chǎn)生更大的影響。

    表2 采用不同模型數(shù)值計(jì)算得到的流固慣性力矩、阻尼力矩和剛度力矩系數(shù)之比Tab.2 Comparison of the fluid-solid ratio of the inertial,damping and stiffness moments coefficients

    為了進(jìn)一步研究振蕩頻率對(duì)結(jié)構(gòu)水彈性響應(yīng)和流固耦合特性的影響,圖5對(duì)比了外部激勵(lì)頻率為共振頻率(ωe/ωn=1)時(shí),采用Theodorsen模型和Munch模型計(jì)算得到的剛/彈性水翼附加慣性力矩、阻尼力矩和剛度力矩隨時(shí)間的演變情況。

    圖5 剛/彈性水翼的慣性力矩、阻尼力矩和剛度矩(黑色曲線)及流體附加慣性、附加阻尼和附加剛度力矩(紅色曲線為基于Munch模型的數(shù)值計(jì)算結(jié)果,藍(lán)色曲線為基于Theodorsen模型的數(shù)值計(jì)算結(jié)果)隨時(shí)間的演變情況(ωe/ωn=1)Fig.5 Comparison of the predicted fluid to solid inertial moments,damping moments and stiffness moments for the steel/flexible hydrofoil in water atωe/ωn=1

    由于Theodorsen模型和Munch模型對(duì)附加慣性力矩系數(shù)定義相同,由圖5(a)—(b)可以看出,不同模型數(shù)值計(jì)算得到的附加慣性力矩隨時(shí)間的演變情況相同,且由于剛性水翼的密度遠(yuǎn)大于流體的密度,而彈性水翼的密度與流體的密度相當(dāng),因此剛性水翼的附加慣性力矩小于結(jié)構(gòu)的慣性力矩,而彈性水翼的附加慣性力矩大于結(jié)構(gòu)的慣性力矩。由于附加慣性力矩系數(shù)與折合頻率無關(guān),由表2可知,外部激勵(lì)頻率為共振頻率時(shí),數(shù)值計(jì)算得到的剛/彈性水翼流固慣性力矩系數(shù)之比分別為0.245和1.718,與外部激勵(lì)頻率為非共振頻率時(shí)的計(jì)算結(jié)果一致。值得注意的是,基于無黏勢(shì)流假設(shè)的Theodorsen模型是針對(duì)任意變形的翼型簡(jiǎn)諧運(yùn)動(dòng)推導(dǎo)結(jié)果,當(dāng)激勵(lì)頻率較高,結(jié)構(gòu)變形量較大時(shí),采用Theodorsen模型的數(shù)值計(jì)算結(jié)果可靠性大大降低。由圖5(c)—(f)不難發(fā)現(xiàn),基于Munch模型的數(shù)值計(jì)算結(jié)果表明,結(jié)構(gòu)阻尼力矩小于流體附加阻尼力矩,且二者同相分步;結(jié)構(gòu)剛度矩大于流體附加剛度矩,且二者呈反相分布,這與外部激勵(lì)頻率為非共振頻率時(shí)數(shù)值計(jì)算得到的水彈性響應(yīng)趨勢(shì)是一致的。由表2可知,與外部激勵(lì)頻率為非共振頻率時(shí)的數(shù)值計(jì)算結(jié)果相比,共振激勵(lì)頻率下,采用Munch模型計(jì)算所得的流固阻尼力矩系數(shù)比和剛度距系數(shù)比明顯減小,流體的附加阻尼和附加剛度作用減弱。此外,結(jié)合圖3(d)可知,當(dāng)結(jié)構(gòu)所受外部激勵(lì)頻率為共振頻率時(shí),結(jié)構(gòu)的變形量呈現(xiàn)發(fā)散趨勢(shì),流體的附加載荷也隨時(shí)間的發(fā)展逐漸增大。綜合圖5(a)—(f)可知,當(dāng)外部激勵(lì)頻率為共振頻率時(shí),彈性水翼的慣性、剛度力矩與流體附加慣性力矩大小相當(dāng),說明除結(jié)構(gòu)的剛度作用外,流體與結(jié)構(gòu)的慣性作用對(duì)結(jié)構(gòu)水彈性響應(yīng)和流固耦合特性的影響也較大。

    3 結(jié) 論

    本文基于完全耦合算法對(duì)二維NACA0009水翼模型進(jìn)行了數(shù)值模擬研究,分析了剛性和彈性水翼在流固耦合作用下的水彈性響應(yīng),并對(duì)影響繞水翼流固耦合特性的因素進(jìn)行了研究。主要結(jié)論如下:

    (1)由于空氣的黏度較小,采用Theodorsen模型和Munch模型對(duì)氣動(dòng)彈性響應(yīng)的數(shù)值計(jì)算結(jié)果基本一致;而水的黏性不可忽略,不考慮流體黏性的Theodorsen模型在一定程度上低估了結(jié)構(gòu)在水中的變形量。

    (2)結(jié)構(gòu)的慣性、阻尼和剛度力矩與流體的相應(yīng)附加載荷均處于同一數(shù)量級(jí),因此流體與結(jié)構(gòu)的相互作用不可忽略,尤其對(duì)于彈性水翼,流體的附加慣性、附加阻尼作用增大,流固耦合算法的數(shù)值穩(wěn)定性對(duì)彈性水翼流固耦合特性的數(shù)值計(jì)算結(jié)果影響也將增大。

    (3)當(dāng)外部激勵(lì)頻率為非共振頻率時(shí),結(jié)構(gòu)的剛度作用是影響結(jié)構(gòu)水彈性響應(yīng)的主要因素,而當(dāng)外部激勵(lì)頻率為共振頻率時(shí),流體的附加阻尼和附加剛度作用減弱,除結(jié)構(gòu)的剛度作用外,流體與結(jié)構(gòu)的慣性作用對(duì)結(jié)構(gòu)水彈性響應(yīng)和流固耦合特性的影響也較大。

    [1]Todd O N,Thomas W S.Aeroelastic response ofa rigid wing supported by nonlinear springs[J].Journal of Aircraft,1998, 35(4):616-622.

    [2]So R M C,Jadic I,Mignolet M P.Fluid-structure resonance produced by oncoming alternating vortices[J].Journal of Fluids and Structures,1999,13(4):519-548.

    [3]Kamakoti R,Shyy W.Fluid-structure interaction for aeroelastic applications[J].Progress in Aerospace Sciences,2005,40:535-558.

    [4]崔 鵬,韓景龍.新型運(yùn)輸機(jī)機(jī)翼的顫振特性分析[J].振動(dòng)工程學(xué)報(bào),2011,24(2):192-198.Cui Peng,Han Jinglong.Flutter analysis of new transport-type wings[J].Journal of Vibration Engineering,2011,24(2):192-198.

    [5]葉永林,吳有生,鄒明松,倪啟軍.基于水彈性力學(xué)的SWATH船結(jié)構(gòu)振動(dòng)與噪聲分析[J].船舶力學(xué),2013,17(4):430-438.Ye Yonglin,Wu Yousheng,et al.Analysis of the structural vibration and noise radiation of a SWATH ship based on hydroelastic method[J].Journal of Ship Mechanics,2013,17(4):430-438.

    [6]Chae E J.Dynamic response and stability of flexible hydrofoils in incompressible and viscous flow[D].University of Michigan,2015.

    [7]Gowing S,Coffin P,Dai C.Hydrofoil cavitation improvements with elastically coupled composite materials[C]//Proceeding of 25th American Towing Tank Conference,Iowa City,USA,1998.

    [8]Young Y L,Kinnas S A.Numerical modeling of supercavitating propeller flows[J].Journal of Ship Research,2003,47(1):48-62.

    [9]Ausoni P,Escaler X,Avellan F.Cavitation influenced on von karman vortex shedding and induced hydrofoil vibrations[J]. ASME Journal of Fluids Engineering,2007,129:966-973.

    [10]Hubner B,Walhorn E,Dinkler D.A monolithic approach to fluid-structure interaction using space-time finite elements [J].Computer Methods in Applied Mechanics and Engineering,2004,193:2087-2104.

    [11]Farhat C,Zee K vander,Geuzaine Ph.Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear aeroelasticity[J].Computer Methods in Appllied Mechanical Engineering,2006,195(17):1973-2001.

    [12]Wu Q,Huang B,Wang G.Experimental and numerical investigation of hydroelastic response of a flexible hydrofoil in cavitating flow[J].International Journalof Multiphase Flow,2015,74:19-33.

    [13]Matthies H G,Steindorf J.Partitioned strong coupling algorithms for fluid-structure interaction[J].Computers and Structures,2003,81:805-812.

    [14]Dowell E H,Hall K C.Modeling of fluid-structure interaction[J].Annual Review of Fluid Mechanics,2001,33:445-490.

    [15]Ryzhakov P B,Rossi R,Idelsohn S R,O?ate E.A monolithic Lagrangian approach for fluid-structure interaction problems[J].Computational Mechanics,2010,46(6):883-899.

    [16]Michler C,Hulshoff S J,Van Brummelen E H,De Borst R.A monolithic approach to fluid-structure interaction[J].Computers&Fluids,2004,33(5):839-848.

    [17]Theodorsen T H.Generaltheory ofaerodynamic instability and the mechanism offlutter[R].NACA Rep.496,1935.

    [18]Munch C,Ausoni P,Braun O,Farhat M,Avellan F.Fluid-structure coupling for an oscillating hydrofoil[J].Journal of Fluids and Structures,2010,26:1018-1033.

    [19]Abramowitz M,Stegun I A.Handbook of mathematical functions[D].National Bureau of Standards,Applied Math,Series #55,Dover Publications,1965.

    [20]Sears W R.Some aspects of non-stationary airfoil theory and its practical applications[J].Journal of the Aeronautical Sciences,1941,8:104-108.

    [21]Crank J,Nicolson P.A practical method for numerical evaluations of partial differential equations of the heat-conduction type[J].Mathematical Proceedings ofthe Cambridge Philosophical Society,Cambridge Press,1947,43(01):50-67.

    Fluid structure interaction analysis of a hydrofoil based on fully coupled algorithm

    WU Qin1,HUANG Biao2,WANG Guo-yu2,CAO Shu-liang1
    (1.Department of Thermal Engineering,Tsinghua University,Beijing 100084,China;2.School of Mechanical Engineering,Beijing Institute of Technology,Beijing 100081,China)

    A fully coupled algorithm for modeling of fluid structure interaction in viscous flow around the NACA0009hydrofoil is presented.The numerical simulations are performed by using the Theodorsen model and the Munch model to investigate the fluid structure interaction and characterize the factors that significantly affect the hydroelastic responses.It is revealed that Theodorsen’s approximation of the hydroelastic moment is based on inviscid,potential flow theory and hence it underestimated the hydroelastic responses to some extent.The fluid inertial,damping,and stiffness forces are not negligible compared to their structural counterparts.For the flexible hydrofoil,the fluid inertial and damping forces increase,suggesting that the numerical instability of fluid structure coupling algorithms may lead to a more significant impact on the numerical prediction of the fluid structure interactions.Meanwhile,the structural stiffness dominates the response of the hydrofoil when the external excitation frequency is set to be the non-resonance frequency,while as the external excitation frequency equals to the resonance frequency,the effect of the fluid damp-ing and stiffness is less significant and the structural response of this case is dominated by the fluid and solid inertial effects in addition to the structural stiffness.

    fully coupled algorithm;hydroelastic;fluid structure interaction

    TV131.2

    A

    10.3969/j.issn.1007-7294.2017.07.002

    1007-7294(2017)07-0804-10

    2017-01-03

    國家自然科學(xué)基金資助項(xiàng)目(51679005);北京市自然科學(xué)基金資助項(xiàng)目(3172029);高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(20131101120014)

    吳 欽(1989-),女,博士,助理研究員,E-mail:wuqin919@163.com;

    黃 彪(1985-),男,副教授,碩士生導(dǎo)師。

    猜你喜歡
    水翼慣性阻尼
    你真的了解慣性嗎
    沖破『慣性』 看慣性
    波浪滑翔機(jī)橢圓形后緣水翼動(dòng)力特性研究
    N維不可壓無阻尼Oldroyd-B模型的最優(yōu)衰減
    關(guān)于具有阻尼項(xiàng)的擴(kuò)散方程
    具有非線性阻尼的Navier-Stokes-Voigt方程的拉回吸引子
    袖珍水翼突防潛艇的設(shè)計(jì)構(gòu)想及運(yùn)用研究
    無處不在的慣性
    三維扭曲水翼空化現(xiàn)象CFD模擬
    普遍存在的慣性
    国产乱来视频区| 夜夜看夜夜爽夜夜摸| 国产男女内射视频| 少妇的逼水好多| 久久精品国产a三级三级三级| 亚洲欧美清纯卡通| 99九九线精品视频在线观看视频| av免费在线看不卡| 亚洲怡红院男人天堂| 亚洲精品日韩av片在线观看| 亚洲精品,欧美精品| 国产成人精品久久久久久| 亚洲美女搞黄在线观看| 一区二区三区免费毛片| 观看av在线不卡| 在线观看免费日韩欧美大片 | 国产视频内射| 精品人妻熟女av久视频| 日本欧美国产在线视频| 综合色丁香网| 9色porny在线观看| 亚洲欧洲日产国产| 激情五月婷婷亚洲| 欧美国产精品一级二级三级 | 国产欧美日韩一区二区三区在线 | 国产91av在线免费观看| 国产精品国产av在线观看| 亚洲成色77777| 伊人亚洲综合成人网| 亚洲美女视频黄频| 日韩成人av中文字幕在线观看| 夜夜骑夜夜射夜夜干| 国产精品一二三区在线看| 中文在线观看免费www的网站| 亚洲精华国产精华液的使用体验| 欧美精品一区二区大全| 亚洲一级一片aⅴ在线观看| 人人妻人人添人人爽欧美一区卜| 男人和女人高潮做爰伦理| 丝袜在线中文字幕| kizo精华| 交换朋友夫妻互换小说| 国产成人免费观看mmmm| 亚洲成人手机| 国产精品一二三区在线看| 国产一区有黄有色的免费视频| 日本黄色日本黄色录像| 免费人妻精品一区二区三区视频| 最近最新中文字幕免费大全7| 亚洲情色 制服丝袜| 我要看黄色一级片免费的| 国产在视频线精品| 香蕉精品网在线| 男女无遮挡免费网站观看| 人人妻人人澡人人爽人人夜夜| 亚洲真实伦在线观看| 18禁裸乳无遮挡动漫免费视频| 久久人人爽av亚洲精品天堂| 97精品久久久久久久久久精品| 国产色爽女视频免费观看| 一个人看视频在线观看www免费| 热re99久久精品国产66热6| 日韩成人av中文字幕在线观看| 人妻一区二区av| 赤兔流量卡办理| 欧美成人午夜免费资源| 十八禁网站网址无遮挡 | 一级毛片我不卡| 亚洲精品中文字幕在线视频 | 男女无遮挡免费网站观看| 女性生殖器流出的白浆| 日韩精品有码人妻一区| 亚洲欧洲精品一区二区精品久久久 | 亚洲国产精品成人久久小说| 中文字幕人妻熟人妻熟丝袜美| 久久这里有精品视频免费| 亚洲丝袜综合中文字幕| 久久这里有精品视频免费| 你懂的网址亚洲精品在线观看| 欧美精品一区二区免费开放| 亚洲av电影在线观看一区二区三区| av不卡在线播放| 国产精品久久久久久精品电影小说| 日本av手机在线免费观看| 少妇裸体淫交视频免费看高清| 久久青草综合色| 一级毛片久久久久久久久女| 日韩中文字幕视频在线看片| 九色成人免费人妻av| 亚洲精品国产成人久久av| 国产欧美另类精品又又久久亚洲欧美| 久热这里只有精品99| 久久久国产欧美日韩av| 伦理电影大哥的女人| av在线观看视频网站免费| 日韩中文字幕视频在线看片| 一本久久精品| 亚洲精品一区蜜桃| 伦精品一区二区三区| 一级,二级,三级黄色视频| 亚洲av中文av极速乱| 日韩中文字幕视频在线看片| 欧美日韩av久久| 日韩 亚洲 欧美在线| 一个人免费看片子| 国产av码专区亚洲av| 日韩不卡一区二区三区视频在线| 免费观看无遮挡的男女| 18禁动态无遮挡网站| 人人妻人人看人人澡| 日本-黄色视频高清免费观看| 99久久精品热视频| 一区二区三区免费毛片| 国产色爽女视频免费观看| 国产黄片美女视频| 日本vs欧美在线观看视频 | 美女cb高潮喷水在线观看| 狠狠精品人妻久久久久久综合| 91精品一卡2卡3卡4卡| 大话2 男鬼变身卡| 精品一区二区免费观看| 好男人视频免费观看在线| 最后的刺客免费高清国语| 美女视频免费永久观看网站| 成人亚洲精品一区在线观看| 欧美日韩在线观看h| 国产精品三级大全| 最近2019中文字幕mv第一页| 日本91视频免费播放| 国产黄片美女视频| 欧美精品一区二区免费开放| av免费观看日本| 日本av免费视频播放| 日产精品乱码卡一卡2卡三| 九九久久精品国产亚洲av麻豆| 久久久国产一区二区| 国产免费一区二区三区四区乱码| 中文乱码字字幕精品一区二区三区| 春色校园在线视频观看| 91精品国产九色| 免费少妇av软件| 香蕉精品网在线| 久久精品国产亚洲av天美| 99热这里只有是精品50| 国产淫语在线视频| a级片在线免费高清观看视频| 两个人免费观看高清视频 | 亚洲av二区三区四区| 欧美日韩综合久久久久久| 国产精品久久久久久久电影| 午夜福利在线观看免费完整高清在| 日本av免费视频播放| 成人国产麻豆网| 亚洲欧美成人综合另类久久久| 精品国产一区二区久久| 亚洲综合色惰| 亚州av有码| 国产亚洲最大av| 久久国产亚洲av麻豆专区| 欧美日韩国产mv在线观看视频| 2021少妇久久久久久久久久久| 国产精品久久久久久精品古装| 熟妇人妻不卡中文字幕| 97在线人人人人妻| 日韩欧美一区视频在线观看 | 久热久热在线精品观看| 国产日韩欧美在线精品| 国产精品一区二区三区四区免费观看| 亚洲精品日韩在线中文字幕| 日本vs欧美在线观看视频 | 五月开心婷婷网| 欧美精品一区二区免费开放| 噜噜噜噜噜久久久久久91| 色婷婷久久久亚洲欧美| 涩涩av久久男人的天堂| 女人久久www免费人成看片| 欧美日韩在线观看h| 亚洲国产精品一区三区| 欧美日韩av久久| 成人综合一区亚洲| 寂寞人妻少妇视频99o| 五月玫瑰六月丁香| 精品卡一卡二卡四卡免费| 香蕉精品网在线| 午夜福利在线观看免费完整高清在| 久久午夜综合久久蜜桃| av黄色大香蕉| 亚洲一区二区三区欧美精品| 亚洲av国产av综合av卡| 精品卡一卡二卡四卡免费| 如何舔出高潮| 国产精品不卡视频一区二区| 寂寞人妻少妇视频99o| 亚洲伊人久久精品综合| 在线观看美女被高潮喷水网站| 国产成人aa在线观看| 免费人成在线观看视频色| 人体艺术视频欧美日本| 亚洲欧洲国产日韩| 久久久久人妻精品一区果冻| 五月天丁香电影| av线在线观看网站| 欧美日本中文国产一区发布| 亚洲真实伦在线观看| 多毛熟女@视频| 2021少妇久久久久久久久久久| av在线app专区| 老司机影院毛片| 最黄视频免费看| 国产综合精华液| 国产中年淑女户外野战色| 国产色婷婷99| av天堂中文字幕网| 久久久久精品久久久久真实原创| 日本av手机在线免费观看| 久久精品国产亚洲网站| 国产淫语在线视频| 国产免费一级a男人的天堂| 一本—道久久a久久精品蜜桃钙片| 久久国产亚洲av麻豆专区| 免费看日本二区| 国产伦在线观看视频一区| 纵有疾风起免费观看全集完整版| 我的女老师完整版在线观看| 大片免费播放器 马上看| 亚洲精品456在线播放app| 中文资源天堂在线| a级毛片在线看网站| 777米奇影视久久| 男女边吃奶边做爰视频| 午夜91福利影院| 熟妇人妻不卡中文字幕| 亚洲中文av在线| 国产精品福利在线免费观看| 日韩三级伦理在线观看| 人人妻人人添人人爽欧美一区卜| 婷婷色综合www| 一级黄片播放器| 亚洲美女搞黄在线观看| 男人和女人高潮做爰伦理| 特大巨黑吊av在线直播| 嘟嘟电影网在线观看| 亚洲精品一区蜜桃| 国产又色又爽无遮挡免| 极品少妇高潮喷水抽搐| 日韩一本色道免费dvd| 免费看不卡的av| 肉色欧美久久久久久久蜜桃| 深夜a级毛片| 日韩av免费高清视频| 美女福利国产在线| 欧美精品亚洲一区二区| 妹子高潮喷水视频| 日韩精品免费视频一区二区三区 | 中文字幕精品免费在线观看视频 | 九草在线视频观看| 夫妻性生交免费视频一级片| 五月伊人婷婷丁香| 亚洲一级一片aⅴ在线观看| 黑人猛操日本美女一级片| 午夜福利,免费看| 国产成人午夜福利电影在线观看| 亚洲性久久影院| 日韩电影二区| 欧美激情极品国产一区二区三区 | 亚洲av成人精品一区久久| 熟女av电影| 天堂中文最新版在线下载| 少妇丰满av| 视频中文字幕在线观看| 最近的中文字幕免费完整| 欧美精品高潮呻吟av久久| 成人免费观看视频高清| av播播在线观看一区| av一本久久久久| 国产老妇伦熟女老妇高清| 极品少妇高潮喷水抽搐| 少妇高潮的动态图| freevideosex欧美| 99九九在线精品视频 | 亚洲国产av新网站| 一本大道久久a久久精品| 两个人免费观看高清视频 | 欧美丝袜亚洲另类| 精品人妻熟女毛片av久久网站| 一本—道久久a久久精品蜜桃钙片| 中文乱码字字幕精品一区二区三区| 尾随美女入室| 极品人妻少妇av视频| 2021少妇久久久久久久久久久| 成人午夜精彩视频在线观看| 日韩中字成人| 男人爽女人下面视频在线观看| 亚洲精华国产精华液的使用体验| 菩萨蛮人人尽说江南好唐韦庄| 成年人午夜在线观看视频| 国产精品99久久久久久久久| 亚洲av男天堂| 亚洲精品成人av观看孕妇| 日韩熟女老妇一区二区性免费视频| 日韩三级伦理在线观看| 亚洲精品中文字幕在线视频 | 亚洲美女黄色视频免费看| 亚洲成人手机| 国产欧美日韩精品一区二区| 日本午夜av视频| 热99国产精品久久久久久7| 久久av网站| 少妇人妻久久综合中文| 精品国产国语对白av| 日韩欧美一区视频在线观看 | 亚洲欧美清纯卡通| 国产一区二区三区综合在线观看 | a级片在线免费高清观看视频| 亚洲欧洲国产日韩| 成年av动漫网址| 亚洲欧美精品专区久久| 久久久久久人妻| 久热这里只有精品99| 欧美亚洲 丝袜 人妻 在线| 国产熟女午夜一区二区三区 | 欧美日韩视频精品一区| 久久久久人妻精品一区果冻| 最近的中文字幕免费完整| 国产色爽女视频免费观看| av网站免费在线观看视频| 视频中文字幕在线观看| av在线app专区| 日韩av在线免费看完整版不卡| av在线播放精品| 乱人伦中国视频| 天天躁夜夜躁狠狠久久av| 91在线精品国自产拍蜜月| 乱码一卡2卡4卡精品| 欧美亚洲 丝袜 人妻 在线| videossex国产| 精品一品国产午夜福利视频| 日本午夜av视频| 国产老妇伦熟女老妇高清| 中国美白少妇内射xxxbb| 日韩精品有码人妻一区| 亚洲欧美一区二区三区国产| 看十八女毛片水多多多| 一边亲一边摸免费视频| 乱码一卡2卡4卡精品| 亚洲怡红院男人天堂| 午夜精品国产一区二区电影| 国内揄拍国产精品人妻在线| av线在线观看网站| 99久久中文字幕三级久久日本| 亚洲不卡免费看| 婷婷色综合大香蕉| 中文字幕免费在线视频6| 国产一区二区三区av在线| 久久精品久久久久久久性| 亚洲真实伦在线观看| 一级黄片播放器| 亚洲精品视频女| 国产 精品1| 久久久a久久爽久久v久久| 自拍偷自拍亚洲精品老妇| 99热6这里只有精品| 日韩av不卡免费在线播放| 欧美日韩av久久| 久久久久久久久久久久大奶| 国产精品99久久99久久久不卡 | 日韩精品有码人妻一区| 亚洲图色成人| 七月丁香在线播放| 国产精品麻豆人妻色哟哟久久| 国产成人免费观看mmmm| 国产日韩欧美在线精品| 色哟哟·www| 精品亚洲成国产av| 亚洲综合色惰| 日本欧美视频一区| 97在线视频观看| 国产伦精品一区二区三区四那| 国国产精品蜜臀av免费| 国产伦理片在线播放av一区| 人人澡人人妻人| 国产成人一区二区在线| 日日爽夜夜爽网站| 欧美日韩综合久久久久久| 插阴视频在线观看视频| 两个人的视频大全免费| 99热6这里只有精品| 久久99精品国语久久久| 亚洲精品中文字幕在线视频 | 韩国av在线不卡| 欧美亚洲 丝袜 人妻 在线| 国产精品一二三区在线看| 黄片无遮挡物在线观看| 中文字幕免费在线视频6| 搡女人真爽免费视频火全软件| 亚洲av在线观看美女高潮| 亚洲怡红院男人天堂| 中文在线观看免费www的网站| 久久久久视频综合| 免费人成在线观看视频色| 国产精品欧美亚洲77777| 免费看不卡的av| 久久综合国产亚洲精品| 成人免费观看视频高清| 国产在视频线精品| 美女大奶头黄色视频| 亚洲欧美日韩另类电影网站| 日本与韩国留学比较| 亚洲国产欧美日韩在线播放 | 久久狼人影院| 精品一区在线观看国产| 久久久欧美国产精品| 中国三级夫妇交换| 99久久精品一区二区三区| 噜噜噜噜噜久久久久久91| 久久久久久久久久成人| 国产精品久久久久久av不卡| 亚洲国产日韩一区二区| 国语对白做爰xxxⅹ性视频网站| 永久网站在线| 国产美女午夜福利| 久久国产乱子免费精品| 超碰97精品在线观看| 日本午夜av视频| 丝瓜视频免费看黄片| 亚洲一级一片aⅴ在线观看| 大又大粗又爽又黄少妇毛片口| 国产男女内射视频| 中国美白少妇内射xxxbb| 精品少妇久久久久久888优播| 人妻制服诱惑在线中文字幕| 一区二区三区乱码不卡18| 免费观看在线日韩| 亚洲欧美日韩另类电影网站| 国产 一区精品| 久久久久久人妻| 国产精品国产三级国产av玫瑰| 日本欧美视频一区| 亚洲国产精品一区三区| h日本视频在线播放| 一级毛片我不卡| 国产女主播在线喷水免费视频网站| 亚洲欧美日韩另类电影网站| 日产精品乱码卡一卡2卡三| 欧美3d第一页| 日本vs欧美在线观看视频 | 色视频在线一区二区三区| 欧美激情国产日韩精品一区| 国产免费一区二区三区四区乱码| 亚洲精品aⅴ在线观看| 午夜激情福利司机影院| 日日啪夜夜撸| 中文字幕制服av| 国产成人91sexporn| 97在线人人人人妻| 多毛熟女@视频| 中文资源天堂在线| 久久久久久久精品精品| 国产高清不卡午夜福利| 99re6热这里在线精品视频| 久久ye,这里只有精品| 欧美+日韩+精品| 一级毛片久久久久久久久女| 黑丝袜美女国产一区| 国产精品久久久久久精品电影小说| 亚洲真实伦在线观看| 国产91av在线免费观看| 中文精品一卡2卡3卡4更新| 赤兔流量卡办理| 久久人人爽av亚洲精品天堂| 中国美白少妇内射xxxbb| 久久影院123| 99久久人妻综合| 国产成人精品一,二区| 亚洲精品久久久久久婷婷小说| 黄色视频在线播放观看不卡| 精品酒店卫生间| 美女国产视频在线观看| 国产成人aa在线观看| 男女免费视频国产| 女性被躁到高潮视频| 亚洲电影在线观看av| 亚洲,一卡二卡三卡| tube8黄色片| 99久久中文字幕三级久久日本| 国产乱人偷精品视频| 日韩制服骚丝袜av| 国产91av在线免费观看| 少妇丰满av| 亚洲不卡免费看| 精品久久久久久电影网| 中文字幕制服av| 最新中文字幕久久久久| 内地一区二区视频在线| 免费高清在线观看视频在线观看| 日本vs欧美在线观看视频 | 日韩欧美 国产精品| 熟女人妻精品中文字幕| 蜜臀久久99精品久久宅男| 免费大片18禁| 日韩欧美一区视频在线观看 | 80岁老熟妇乱子伦牲交| 我的女老师完整版在线观看| 新久久久久国产一级毛片| 我的女老师完整版在线观看| 精品亚洲成国产av| 中国三级夫妇交换| 久久久a久久爽久久v久久| 又粗又硬又长又爽又黄的视频| 国产高清三级在线| 日本午夜av视频| 亚洲精品中文字幕在线视频 | 毛片一级片免费看久久久久| 人妻制服诱惑在线中文字幕| 久久97久久精品| 女性生殖器流出的白浆| 嫩草影院入口| 精品人妻熟女毛片av久久网站| 美女xxoo啪啪120秒动态图| 日韩视频在线欧美| 大码成人一级视频| 老女人水多毛片| 日韩不卡一区二区三区视频在线| 高清毛片免费看| 久久精品久久久久久噜噜老黄| 国产深夜福利视频在线观看| 搡女人真爽免费视频火全软件| 免费看av在线观看网站| 日本午夜av视频| 嘟嘟电影网在线观看| 精品久久久久久久久亚洲| 内地一区二区视频在线| www.色视频.com| 91久久精品国产一区二区三区| 久久这里有精品视频免费| av免费在线看不卡| 精品视频人人做人人爽| av一本久久久久| 国产在线男女| 噜噜噜噜噜久久久久久91| 91精品一卡2卡3卡4卡| 综合色丁香网| 欧美日韩一区二区视频在线观看视频在线| 五月天丁香电影| 日韩中字成人| 99久久精品一区二区三区| 热re99久久国产66热| 久久免费观看电影| 伊人久久国产一区二区| 一级二级三级毛片免费看| 欧美日韩一区二区视频在线观看视频在线| 啦啦啦啦在线视频资源| 在线观看三级黄色| 国产在线免费精品| 国产成人aa在线观看| 九九久久精品国产亚洲av麻豆| 欧美高清成人免费视频www| 一级毛片我不卡| 永久免费av网站大全| 人妻系列 视频| 亚洲精品亚洲一区二区| 性色av一级| 搡女人真爽免费视频火全软件| 久久久a久久爽久久v久久| 一级a做视频免费观看| 尾随美女入室| 久久韩国三级中文字幕| 国产精品一区二区三区四区免费观看| 99九九在线精品视频 | 高清欧美精品videossex| 高清午夜精品一区二区三区| 国产片特级美女逼逼视频| 天堂中文最新版在线下载| 成人漫画全彩无遮挡| 韩国高清视频一区二区三区| 亚洲精品久久久久久婷婷小说| 嘟嘟电影网在线观看| 丰满人妻一区二区三区视频av| 赤兔流量卡办理| 亚洲精品乱久久久久久| 中文字幕久久专区| 蜜臀久久99精品久久宅男| 中文乱码字字幕精品一区二区三区| 亚洲欧美精品专区久久| 香蕉精品网在线| 亚洲伊人久久精品综合| 老司机亚洲免费影院| 国产午夜精品久久久久久一区二区三区| 免费黄网站久久成人精品| 亚洲内射少妇av| 国产欧美亚洲国产| 男女国产视频网站| 内射极品少妇av片p| 亚洲综合精品二区| 中文天堂在线官网| 建设人人有责人人尽责人人享有的| 欧美日韩视频精品一区| 黑人巨大精品欧美一区二区蜜桃 | 大又大粗又爽又黄少妇毛片口| 国产精品一区二区在线不卡| 国产成人freesex在线| 欧美日韩视频高清一区二区三区二| 人妻制服诱惑在线中文字幕| 最近最新中文字幕免费大全7| 日韩在线高清观看一区二区三区| 国产成人a∨麻豆精品| 人妻夜夜爽99麻豆av| 最近2019中文字幕mv第一页| av福利片在线| 久久狼人影院| 成人漫画全彩无遮挡| 黄色日韩在线| 日韩亚洲欧美综合| 精品熟女少妇av免费看| 麻豆成人午夜福利视频|