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

    艇槳耦合狀態(tài)螺旋槳水動(dòng)力與噪聲數(shù)值預(yù)報(bào)方法研究

    2021-11-26 03:44:22黃苗苗
    船舶力學(xué) 2021年11期
    關(guān)鍵詞:螺旋槳聲學(xué)潛艇

    張 楠,李 亞,黃苗苗,陳 默

    (中國船舶科學(xué)研究中心a.水動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室;b.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室,江蘇無錫 214082)

    0 引 言

    在水動(dòng)力學(xué)領(lǐng)域,潛艇模型快速性的數(shù)值模擬工作一般包括潛艇模型阻力數(shù)值模擬、推進(jìn)器敞水水動(dòng)力數(shù)值模擬以及帶推進(jìn)器潛艇自航水動(dòng)力數(shù)值模擬這三項(xiàng)主要內(nèi)容,流場數(shù)值模擬可伴隨這三項(xiàng)內(nèi)容開展,其中自航工況數(shù)值模擬的技術(shù)難點(diǎn)在于準(zhǔn)確預(yù)報(bào)艇-槳相互作用,艇-槳相互作用也可以簡稱為艇-槳干擾或艇-槳耦合。在完成上述阻力、敞水、自航三項(xiàng)數(shù)值模擬內(nèi)容之后,可按照基于模型試驗(yàn)的快速性預(yù)報(bào)方法,分析自航因子,外推預(yù)報(bào)實(shí)艇功率和航速。

    長期以來,潛艇艇-槳干擾數(shù)值模擬一直是國際水動(dòng)力學(xué)領(lǐng)域的重要內(nèi)容。國際上對(duì)于艇-槳干擾流動(dòng)與水動(dòng)力的數(shù)值模擬進(jìn)行了持續(xù)不斷的研究,經(jīng)歷了長久的發(fā)展與驗(yàn)證,模擬方法從早期的用分布體積力代替螺旋槳的方法已經(jīng)發(fā)展到目前常用的用滑移網(wǎng)格模擬艇后真實(shí)的螺旋槳運(yùn)轉(zhuǎn)的方法。近年來,國際上特別是美國,在此領(lǐng)域又取得了許多有價(jià)值的研究成果。

    Alin 等(2010)[1]采用大渦模擬方法(LES)對(duì)全附體艦船與潛艇流動(dòng)及流激噪聲進(jìn)行了計(jì)算嘗試,計(jì)算了SUBOFF 潛艇模型(DARPA AFF8)帶七葉大側(cè)斜槳(INSEAN E1619)的繞流,隨后又結(jié)合Light?hill聲學(xué)類比計(jì)算了流激噪聲。他對(duì)于螺旋槳的處理沒有采用常用的滑移網(wǎng)格方法,而是采用了一種動(dòng)網(wǎng)格類型的變形重生方法(D&R)。

    Chase 等(2012)[2]采用基于重疊網(wǎng)格的求解器CFDShip-Iowa(V4.5)對(duì)SUBOFF 潛艇拖曳狀態(tài)、自航狀態(tài)以及自航模(自由自航)超越機(jī)動(dòng)狀態(tài)的水動(dòng)力和流場進(jìn)行了數(shù)值模擬,作者將MIT 開發(fā)的螺旋槳?jiǎng)萘髑蠼獬绦騊UF-14也嵌入在該軟件中。在艇槳干擾模擬中,采用的計(jì)算模型是全附體潛艇帶E1619 螺旋槳。在自航與自航模操縱模擬中,對(duì)于螺旋槳采用滑移網(wǎng)格與PUF-14 計(jì)算這兩種方法,直線航行時(shí),兩種方法差別不大,但當(dāng)操縱運(yùn)動(dòng)時(shí),特別在尾部流動(dòng)分離與低進(jìn)速系數(shù)下,兩種方法的差別比較明顯。

    Liefvendahl等(2012)[3]采用大渦模擬方法計(jì)算研究了全附體SUBOFF潛艇模型帶E1619螺旋槳的水動(dòng)力和負(fù)荷脈動(dòng)問題,螺旋槳模擬采用動(dòng)網(wǎng)格方法。作者的數(shù)值模擬捕捉到了大尺度非定常相干結(jié)構(gòu),計(jì)算了圍殼與尾翼流動(dòng)結(jié)構(gòu)以及艇體邊界層在尾部的相互干擾流動(dòng)特征,分析了在尾部逆壓梯度作用下螺旋槳入流的形成過程,將這些流動(dòng)特征與螺旋槳推力、扭矩以及葉片上的負(fù)荷進(jìn)行了相關(guān),并與單獨(dú)艇體和單獨(dú)螺旋槳的計(jì)算結(jié)果進(jìn)行了對(duì)比分析。

    Kim 等(2014)[4]發(fā)展了一種耦合剛體運(yùn)動(dòng)(RBM)與非定常RANS(URANS)的數(shù)值模擬方法,對(duì)于SUBOFF潛艇帶E1619槳的拘束模狀態(tài)與自航模狀態(tài)繞流進(jìn)行了數(shù)值模擬。在拘束模偏航計(jì)算中,利用移動(dòng)參考系和滑移網(wǎng)格兩種方法來考慮螺旋槳運(yùn)轉(zhuǎn)的影響。在自航模回轉(zhuǎn)計(jì)算中,利用驅(qū)動(dòng)盤模型和滑移網(wǎng)格方法來考慮螺旋槳運(yùn)轉(zhuǎn)的影響,計(jì)算了三自由度的艇體運(yùn)動(dòng),分析了艇、槳、舵三者的相互干擾。

    Norrison 等(2016)[5]采用大渦模擬方法對(duì)實(shí)尺度全附體Joubert 潛艇的艇-槳干擾流動(dòng)進(jìn)行了數(shù)值模擬。雷諾數(shù)達(dá)到Re= 1.68 × 108與Re= 3.36 × 108。結(jié)果表明,繞實(shí)艇的流動(dòng)拓?fù)浣Y(jié)構(gòu)與模型尺度相似,符合隨雷諾數(shù)增高流動(dòng)結(jié)構(gòu)趨于緊致與收縮的經(jīng)典認(rèn)識(shí)。螺旋槳采用DSTG115-1 五葉槳與DSTG057-1七葉槳,與五葉槳相比,七葉槳的葉中負(fù)荷更高,梢渦破碎更早,尾流摻混更強(qiáng)。

    Carrica 等(2016)[6]對(duì)于Joubert BB2 潛艇自航模操縱工況進(jìn)行了試驗(yàn)和數(shù)值模擬研究,網(wǎng)格數(shù)達(dá)到3550 萬。自航模試驗(yàn)在MARIN 水池開展,包括近水面與深潛直航、回轉(zhuǎn)、垂直面與水平面Z 形機(jī)動(dòng)、應(yīng)急上浮等。他們將試驗(yàn)數(shù)據(jù)做成數(shù)據(jù)庫,用以驗(yàn)證操縱性預(yù)報(bào)方法的準(zhǔn)確性。作者采用兩種求解器進(jìn)行潛艇自航模操縱工況數(shù)值計(jì)算,一種是ReFRESCO(荷蘭MARIN牽頭,葡萄牙IST、巴西USPTPN、荷蘭DUT、荷蘭RuG、英國UoS、瑞典CUT、加拿大UNB 等大學(xué)和科研院所聯(lián)合開發(fā)的軟件),另一種是REX(CFDShip-Iowa V4.5+Magnus)。采用體積力法和滑移網(wǎng)格法模擬螺旋槳旋轉(zhuǎn),采用重疊嵌套網(wǎng)格模擬潛艇運(yùn)動(dòng)。對(duì)于深潛狀態(tài)三個(gè)航速下的潛艇直航運(yùn)動(dòng)、近水面垂向控制運(yùn)動(dòng)、操尾舵定深回轉(zhuǎn)運(yùn)動(dòng)、操圍殼舵與尾舵的Z形機(jī)動(dòng)、應(yīng)急上浮等都進(jìn)行了詳細(xì)的數(shù)值模擬、驗(yàn)證和分析,研究工作全面而且深入。

    Carrica 等(2018)[7]利用1700 萬到1 億3 千萬網(wǎng)格對(duì)于Joubert BB2 近水面波浪狀態(tài)下運(yùn)動(dòng)特性進(jìn)行了進(jìn)一步的數(shù)值模擬分析,波浪狀態(tài)主要考慮的是頂浪規(guī)則波。對(duì)于靜止與波浪狀態(tài)四個(gè)圍殼浸深的潛艇近水面自由自航特性進(jìn)行了數(shù)值模擬,利用滑移網(wǎng)格模擬螺旋槳旋轉(zhuǎn),采用重疊嵌套網(wǎng)格模擬潛艇運(yùn)動(dòng),潛艇尾舵為X 型舵,作者采用PID 控制器實(shí)現(xiàn)對(duì)于尾舵的操舵控制。Carrica 等(2016,2018)[6-7]的工作達(dá)到了很高的研究水平,可以稱為目前船舶領(lǐng)域CFD技術(shù)的標(biāo)桿。

    隨著國際上對(duì)艇-槳干擾流動(dòng)數(shù)值模擬研究的蓬勃開展,自航狀態(tài)船舶尾后螺旋槳輻射噪聲的數(shù)值預(yù)報(bào)問題也逐漸進(jìn)入研究者視野。近年來,也有一些研究問世,但相較艇-槳干擾流動(dòng)問題的模擬研究而言,艇-槳干擾流激噪聲的模擬研究還是比較少的,而且國外在計(jì)算螺旋槳噪聲時(shí)常用的是Lighthill所創(chuàng)立的聲學(xué)類比方法,近年來主要用的是FW-H聲學(xué)類比方法。

    Lighthill(1952)[8]創(chuàng)立了聲學(xué)類比方法,起初主要處理的是噴射噪聲問題。在聲學(xué)類比方法里,流動(dòng)特征(聲源)通過求解合適的非定常流動(dòng)方程得到,控制方程可以是可壓縮的流動(dòng)方程,也可以是不可壓縮的流動(dòng)方程,將求出的流動(dòng)聲源代入聲學(xué)類比方程的遠(yuǎn)場解中,通過格林公式來預(yù)報(bào)遠(yuǎn)場噪聲。聲學(xué)類比方法的研究思路就是使用CFD 技術(shù)求解非定常流場,作為等效聲源項(xiàng)輸入到聲場計(jì)算中,然后利用波動(dòng)方程的解,將聲源項(xiàng)產(chǎn)生的脈動(dòng)進(jìn)一步輻射到外場。從理論上講,只要知道流體運(yùn)動(dòng)的物理特性,無論是運(yùn)動(dòng)源還是力源,從Lighthill 方程出發(fā)都能求解其聲輻射,流體運(yùn)動(dòng)與聲輻射看似獨(dú)立的兩個(gè)物理現(xiàn)象從此聯(lián)系起來,流體動(dòng)力聲學(xué)作為一門學(xué)科從此發(fā)展起來。隨后,運(yùn)動(dòng)體引起的流動(dòng)噪聲問題也逐漸進(jìn)入學(xué)者的視野。1955年,Curle[9]采用Kirchhoff方法將Lighthill理論進(jìn)行推廣,導(dǎo)出了著名的Curle 方程,該方程可以處理靜止固體邊界的影響。從Curle 方程可知,一個(gè)四極子源與固體邊界的相互作用會(huì)產(chǎn)生新的低階的偶極子源,輻射效率增強(qiáng),這正反映了聲學(xué)的復(fù)雜性。隨后在1969 年,F(xiàn)fowcs Williams 和Hawkings[10]應(yīng)用極為有效的數(shù)學(xué)工具——廣義函數(shù)法將Curle 方程進(jìn)行拓展,使之可以處理固體邊界在流體中運(yùn)動(dòng)的發(fā)聲問題,得到了經(jīng)典的FW-H 方程。1974年,Gold?stein[11]用格林函數(shù)方法推導(dǎo)出了廣義的Lighthill方程。從這個(gè)方程我們可以清晰地知道Curle 方程和FW-H 方程均是該方程的特定表達(dá)。上述基于Lighthill 思想的各種方法統(tǒng)稱為聲學(xué)類比方法。在近代邏輯學(xué)中,類比法是根據(jù)兩個(gè)(或兩類)不同對(duì)象的部分屬性相似,而推出這兩個(gè)(或兩類)對(duì)象的其他屬性也可能是相似的一種推理方法。類比法是歸納法和演繹法的中間狀態(tài),其推理方式是從特殊到特殊,即從一個(gè)對(duì)象的特殊知識(shí)過渡到另一對(duì)象的特殊知識(shí)。

    在Lighthill方程提出之后,Powell(1964)[12]將渦量描述引入該方程,進(jìn)一步研究了渦運(yùn)動(dòng)和聲產(chǎn)生之間的聯(lián)系,推導(dǎo)出了Powell方程(渦聲方程),開創(chuàng)了渦聲理論,這一理論實(shí)質(zhì)上是Lighthill理論在低馬赫數(shù)下的一個(gè)演化。Powell 方程指出:渦是低馬赫數(shù)下等熵絕熱流動(dòng)發(fā)聲的根源。由于渦量分布往往集中在狹小的流動(dòng)區(qū)域,所以它是緊致的偶極子源。研究者普遍認(rèn)為,Lighthill 理論在預(yù)報(bào)流動(dòng)噪聲方面有實(shí)用價(jià)值,而在探索流動(dòng)發(fā)聲的內(nèi)部機(jī)制方面,Powell理論則以其簡潔深邃的內(nèi)涵顯示出了極大的優(yōu)勢。Howe(1975)[13]進(jìn)一步發(fā)展了渦聲理論,引入駐焓(stagnation enthalpy)的概念,考慮了熵變化和平均流對(duì)流動(dòng)發(fā)聲的影響,導(dǎo)出了描述聲音在氣流中傳播的非齊次方程——Howe 方程。Howe方程是描述由于渦量和熵的變化以及其相互作用而發(fā)聲的聲學(xué)普遍公式,它也描述了與聲音相關(guān)聯(lián)的非線性現(xiàn)象的相互作用。Howe方程指出,如不存在渦旋和熵梯度,則氣流不會(huì)發(fā)聲,聲源只集中在那些存在有渦量及熵梯度的區(qū)域。Howe 方程只有在無旋和等熵時(shí)才是封閉的,多數(shù)情況下,此方程不封閉,只有再引進(jìn)另外的旋度及熵的擾動(dòng)量方程才能求解。Howe(2003)[14]進(jìn)一步指出:氣體的非定常粘流運(yùn)動(dòng)是聲音、渦旋和熵運(yùn)動(dòng)成分的疊加和相互作用組成的,非線性效應(yīng)導(dǎo)致上述運(yùn)動(dòng)的相互轉(zhuǎn)變,這些問題是擺在連續(xù)介質(zhì)力學(xué)和氣動(dòng)聲學(xué)面前的艱巨任務(wù),當(dāng)然,這些內(nèi)容也遠(yuǎn)遠(yuǎn)超出一般流動(dòng)聲學(xué)的研究范疇。

    在FW-H 聲學(xué)類比方程建立之后,許多學(xué)者致力于求出遠(yuǎn)場解。Farassat(1981)[15]、Farassat 與Brentner(1988)等[16]、Francescantonio(1997)[17]、Brentner 與Holland(1997)等[18]、Casalino(2003)[19]在聲學(xué)遠(yuǎn)場解的推導(dǎo)研究方面做出了重要貢獻(xiàn)。在聲學(xué)計(jì)算結(jié)果處理技術(shù)方面,Wang、Lele 與Moin(1996)[20]提出了出口邊界修正方法。Mendez 等(2009)[21]研究了開口與閉口積分面的影響,認(rèn)為使用基于壓力的公式結(jié)合出流盤面平均技術(shù)給出的結(jié)果最好。Nitzkorski 與Mahesh(2014)[22]提出了“dy?namic end cap”處理方法。

    在國際水動(dòng)力學(xué)領(lǐng)域,Bensow等(2016)[23]采用OpenFOAM開源軟件中的大渦模擬與FW-H聲學(xué)類比方法對(duì)在船舶尾后運(yùn)轉(zhuǎn)螺旋槳的空泡噪聲進(jìn)行了計(jì)算研究,并與試驗(yàn)進(jìn)行了對(duì)比。利用VOF 結(jié)合顯式質(zhì)量傳遞模型計(jì)算空泡,模擬中忽略了水面的影響,螺旋槳運(yùn)轉(zhuǎn)采用滑移網(wǎng)格方法。計(jì)算模型包括單獨(dú)船體、單獨(dú)螺旋槳以及帶槳船體,計(jì)算工況包括兩個(gè)狀態(tài)的空泡噪聲,一個(gè)狀態(tài)的無空泡噪聲。對(duì)于某狀態(tài)的空泡噪聲,與試驗(yàn)結(jié)果相比,計(jì)算結(jié)果在400 Hz以下小了約20 dB,在400~1 100 Hz,計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合尚可。

    Ianniello等(2014)[24]采用不可壓縮RANS方法結(jié)合FW-H聲學(xué)類比方程對(duì)于某巡邏艇模型在定常直航下的螺旋槳噪聲與船體噪聲進(jìn)行了預(yù)報(bào)。作者通過數(shù)值計(jì)算明確指出,盡管水中船用螺旋槳轉(zhuǎn)速較低,但也不能忽視四極子噪聲。厚度噪聲、負(fù)荷噪聲分量對(duì)遠(yuǎn)場壓力的貢獻(xiàn)非常有限,主要噪聲源應(yīng)是流場中的速度脈動(dòng),特別是在無空泡產(chǎn)生的情況下,槳葉梢渦的發(fā)放與向下游傳播才是最主要的噪聲源。作者將水動(dòng)力與噪聲結(jié)果進(jìn)行了對(duì)比分析發(fā)現(xiàn),由于船體的存在,出現(xiàn)了明顯的散射。

    Ianniello等(2016)[25]采用FW-H 聲學(xué)類比方法計(jì)算了無空泡狀態(tài)下的螺旋槳噪聲,詳細(xì)分析了面積分對(duì)噪聲計(jì)算結(jié)果的影響。他經(jīng)過理論分析和數(shù)值計(jì)算發(fā)現(xiàn),與航空中旋翼的主要噪聲是厚度噪聲和負(fù)荷噪聲不同,水中螺旋槳的主要噪聲是螺旋槳的尾流/渦流所發(fā)放的噪聲,非定常非均勻的流場是主要噪聲源。

    在國內(nèi),朱錫清(2008)[26]深入研究了螺旋槳噪聲的機(jī)理、種類和預(yù)報(bào)方法,詳細(xì)闡釋了推進(jìn)器噪聲的成因和種類、非空泡螺旋槳離散譜(線譜)噪聲預(yù)報(bào)與寬帶噪聲的理論分析、螺旋槳空泡噪聲的預(yù)報(bào)以及推進(jìn)器(螺旋槳)噪聲的控制等方面的研究方法與成果。熊紫英等[27]采用面元法對(duì)于無空泡螺旋槳非定常推力脈動(dòng)及其誘導(dǎo)線譜噪聲進(jìn)行了計(jì)算分析。

    作者以往對(duì)于潛艇流場、水動(dòng)力與流激噪聲開展了一些數(shù)值模擬研究[28-38]。前期主要是利用RANS 方法對(duì)潛艇帶推進(jìn)器自航的流場與水動(dòng)力進(jìn)行了數(shù)值模擬,在流場與水動(dòng)力模擬的基礎(chǔ)上,近年來又利用LES結(jié)合FW-H聲學(xué)類比、Kirchhoff方法與Powell渦聲方程等聲學(xué)計(jì)算技術(shù)對(duì)潛艇流激噪聲進(jìn)行了數(shù)值預(yù)報(bào)計(jì)算探索,這些流聲耦合研究工作是本文研究的基礎(chǔ)。本文即采用LES結(jié)合Powell渦聲方程對(duì)SUBOFF潛艇帶AU5-65螺旋槳在敞水與自航工況下的流場和流激噪聲進(jìn)行了數(shù)值模擬,研究了LES 結(jié)合Powell 渦聲方程對(duì)螺旋槳水動(dòng)力與噪聲的預(yù)報(bào)能力。本文中的LES 方法利用商用軟件Fluent完成,Powell渦聲方程數(shù)值預(yù)報(bào)方法通過自編程序?qū)崿F(xiàn)。目前國內(nèi)外采用Powell渦聲方程計(jì)算流激噪聲的公開文獻(xiàn)還不多見,作者前期對(duì)于Powell渦聲方程預(yù)報(bào)方法的探索請(qǐng)見文獻(xiàn)[35]。

    1 數(shù)值計(jì)算方法

    1.1 大渦模擬方法

    將N-S方程進(jìn)行網(wǎng)格濾波,從而得到大渦模擬的控制方程。濾波之后,則大于網(wǎng)格體積的流動(dòng)結(jié)構(gòu)通過直接求解N-S方程得到,而那些小于網(wǎng)格體積的流動(dòng)結(jié)構(gòu)則通過亞格子渦模型來進(jìn)行模擬。濾波變量由下式定義:

    式中,φ(x)為流動(dòng)變量,D為流體域,G為濾波函數(shù)。取濾波函數(shù)為

    式中,V為計(jì)算網(wǎng)格的體積。

    濾波后的連續(xù)性方程和N-S方程可以表示為

    式中:σij為分子粘性引起的應(yīng)力張量;τij為亞格子應(yīng)力,需用亞格子模型進(jìn)行模擬。

    本文采用動(dòng)態(tài)Smagorinsky模型對(duì)亞格子應(yīng)力進(jìn)行模擬,這種模型是Germano在1991年提出的[39],通過對(duì)最小可解尺度的信息進(jìn)行采樣,然后利用這些信息來模擬亞格子尺度應(yīng)力,此模型在接近壁面邊界時(shí)給出了正確的漸近特性,而且并不需要阻尼函數(shù)或者間歇函數(shù)。此模型還能夠考慮逆散射的影響。Lilly(1992)[40]利用最小二乘法對(duì)此模型進(jìn)行了改進(jìn),這個(gè)辦法提供了一個(gè)自洽的動(dòng)態(tài)模型,可以在每個(gè)空間網(wǎng)格點(diǎn)、每一時(shí)間步上計(jì)算模型參數(shù)C。

    引入布西內(nèi)斯克(Boussinesq)假設(shè),

    1.2 Powell渦聲理論

    作者對(duì)于Powell 渦聲方程的研究請(qǐng)參見文獻(xiàn)[27],本文作簡要描述。在研究湍流誘發(fā)噪聲問題時(shí),關(guān)鍵的一步就是構(gòu)建流動(dòng)聲源數(shù)學(xué)公式,也就是如何在數(shù)學(xué)表達(dá)上將流體運(yùn)動(dòng)轉(zhuǎn)換為聲源。Lighthill建立的聲學(xué)類比方法是將雷諾應(yīng)力、壓力、剪應(yīng)力進(jìn)行組合作為聲源,通過面積分和體積分得到遠(yuǎn)場輻射噪聲,在工程實(shí)際中發(fā)揮了巨大作用,但聲學(xué)類比理論尚不足以深入了解流動(dòng)發(fā)聲的機(jī)理和細(xì)節(jié)。因?yàn)楸娝苤?,渦會(huì)產(chǎn)生噪聲,而在聲學(xué)類比理論中沒有辨識(shí)出渦動(dòng)力學(xué)特征,且Lighthiil應(yīng)力張量在空間分布比較疏散,不利于計(jì)算,但渦量的分布相對(duì)集中,屬于緊致聲源,便于計(jì)算。在低馬赫數(shù)流動(dòng)中,渦旋結(jié)構(gòu)一般而言都是緊致的,即渦旋結(jié)構(gòu)尺度相對(duì)聲波波長而言是小量。Powell深入研究了流體動(dòng)力與噪聲之間的關(guān)系,并將它們與渦運(yùn)動(dòng)聯(lián)系起來,得到了渦運(yùn)動(dòng)產(chǎn)生噪聲的機(jī)理與控制方程。

    本文求解的Powell渦聲方程表達(dá)為

    關(guān)于Powell渦聲方程的遠(yuǎn)場解,利用三維波動(dòng)方程自由空間格林函數(shù)

    可將Powell方程的遠(yuǎn)場解表達(dá)為

    注意到格林函數(shù)的下列關(guān)系:

    所以可將式(10)改寫為

    進(jìn)一步利用δ函數(shù)的性質(zhì)可得

    1.3 數(shù)值求解

    時(shí)間項(xiàng)采用二階隱式格式離散,動(dòng)量方程采用限界中心差分格式離散,壓力速度耦合采用SIMPLE算法。利用代數(shù)多重網(wǎng)格(AMG)方法加速收斂。計(jì)算中時(shí)間步長Δt=1×10-5s。壁面y+≈0.2~15。采用FFT結(jié)合Hanning窗處理非定常信號(hào)時(shí)間序列。

    1.4 邊界條件

    敞水工況采用邊界條件為

    速度入口:螺旋槳前方5倍槳徑處,設(shè)定來流速度的大小與方向。

    壓力出口:螺旋槳后方10倍槳徑處,設(shè)定相對(duì)于參考?jí)毫c(diǎn)的流體靜壓值。

    壁面:螺旋槳表面,設(shè)定無滑移粘附條件。

    外場:距離螺旋槳5倍槳徑,參數(shù)設(shè)置同速度入口。

    計(jì)算中采用全域模型,不含對(duì)稱面。見圖1(a)。

    自航工況采用邊界條件為

    速度入口:艇艏前方1倍艇長處,設(shè)定來流速度的大小與方向。

    壓力出口:艇艉后方2倍艇長處,設(shè)定相對(duì)于參考?jí)毫c(diǎn)的流體靜壓值。

    壁面:潛艇與螺旋槳表面,設(shè)定無滑移粘附條件。

    外場:距離艇身1倍艇長,參數(shù)設(shè)置同速度入口。

    計(jì)算中采用全域模型,不含對(duì)稱面。見圖1(b)。

    圖1 計(jì)算域Fig.1 Computational domain

    2 計(jì)算模型與網(wǎng)格

    SUBOFF 潛艇模型主體長為4.356 m,其前體長為1.016 m,平行中體長為2.229 m,后體長為1.111 m;艇身最大直徑為0.508 m;槳盤面位于距艇艏4.260 m處;指揮臺(tái)圍殼為一立柱體,其導(dǎo)邊位于距艇艏0.924 m 處,隨邊位于距艇艏1.292 m 處,其水平截面為兩橢圓相交而成;指揮臺(tái)上部為一外凸的頂蓋;四個(gè)尾翼形狀、大小都相同,其橫截面為NACA0020翼型,對(duì)稱布置于艇身上下左右。

    AU5-65 螺旋槳模型直徑為0.24 m,葉數(shù)為5,螺距比為0.782,盤面比為0.65,轂徑比為0.18,后傾角為10°。

    計(jì)算模型與網(wǎng)格如圖2 所示。敞水工況所用網(wǎng)格數(shù)為1 083 萬(靜止區(qū)域453 萬,旋轉(zhuǎn)區(qū)域630萬),自航工況所用網(wǎng)格數(shù)為1 937萬(靜止區(qū)域1 307萬,旋轉(zhuǎn)區(qū)域630萬),網(wǎng)格數(shù)量的確定參照了以往網(wǎng)格收斂性研究經(jīng)驗(yàn)[28-38],自航算例在本單位并行系統(tǒng)上利用10 個(gè)節(jié)點(diǎn)(200 核)計(jì)算,需要1 周時(shí)間達(dá)到穩(wěn)定收斂。包裹螺旋槳的圓柱體內(nèi)網(wǎng)格為四面體非結(jié)構(gòu)化網(wǎng)格,圓柱體之外的網(wǎng)格都為六面體結(jié)構(gòu)化網(wǎng)格。對(duì)于潛艇指揮臺(tái)圍殼、尾翼與螺旋槳周圍的網(wǎng)格進(jìn)行加密處理。

    圖2 潛艇與螺旋槳計(jì)算模型及網(wǎng)格Fig.2 Computational models and meshes of submarine and propeller

    3 計(jì)算結(jié)果分析

    3.1 水動(dòng)力與流場計(jì)算結(jié)果分析

    為了驗(yàn)證計(jì)算方法的有效性,先從計(jì)算水動(dòng)力入手,進(jìn)而計(jì)算渦旋流場,最后計(jì)算流激噪聲,采用層層遞進(jìn)的驗(yàn)證途徑。

    對(duì)于敞水工況,采用與試驗(yàn)一致的等轉(zhuǎn)速變航速的方法進(jìn)行螺旋槳水動(dòng)力數(shù)值模擬。利用滑移網(wǎng)格方法使螺旋槳以設(shè)定的轉(zhuǎn)速n旋轉(zhuǎn),改變多個(gè)來流水速V,從而計(jì)算多個(gè)進(jìn)速系數(shù)J0下的螺旋槳推力T0、轉(zhuǎn)矩Q0。推力系數(shù)、扭矩系數(shù)及敞水效率的定義如下:

    式中,D為螺旋槳直徑,ρ為流體介質(zhì)密度。

    計(jì)算中的轉(zhuǎn)速與敞水試驗(yàn)保持一致為25 r/s,使得螺旋槳在工作點(diǎn)處的雷諾數(shù)大于臨界值3× 105,滿足規(guī)程要求。槳葉雷諾數(shù)為

    式中,ν為流體運(yùn)動(dòng)粘性系數(shù)。

    敞水水動(dòng)力的數(shù)值模擬結(jié)果見圖3,圖中試驗(yàn)結(jié)果為MSU給出的拖曳水池測試結(jié)果。

    圖3 敞水工況水動(dòng)力與渦旋結(jié)構(gòu)計(jì)算結(jié)果Fig.3 Computed hydrodynamic forces and vortical structures in open water condition

    從圖3可以看到,不同進(jìn)速系數(shù)下,與試驗(yàn)結(jié)果相比,推力系數(shù)的計(jì)算誤差為2%~3%,扭矩系數(shù)的計(jì)算誤差為3%~4%,敞水效率的計(jì)算誤差都為1%~2%。從工程預(yù)報(bào)角度來看,大渦模擬方法對(duì)螺旋槳敞水水動(dòng)力(推力、扭矩)的預(yù)報(bào)達(dá)到了較高精度,計(jì)算效果是令人滿意的。由于在大渦模擬方法中,網(wǎng)格比較精細(xì),近壁面網(wǎng)格分辨率較優(yōu),可以直接解算到粘性底層,所以對(duì)于近壁面流動(dòng)的計(jì)算效果較好,推力與扭矩的預(yù)報(bào)精度可以達(dá)到工程實(shí)用要求。

    敞水工況下,J0= 0.7 時(shí),螺旋槳周圍的渦旋結(jié)構(gòu)計(jì)算結(jié)果見圖3,圖中渦旋結(jié)構(gòu)使用Q判據(jù)進(jìn)行捕捉,Q判據(jù)的定義為

    式中,ui、uj為三向速度,xi、xj為三向坐標(biāo)。

    從圖3可以看到,在敞水工況,螺旋槳周圍存在明顯的梢渦、葉根渦。由于在敞水工況計(jì)算中,槳軸延伸至出口,所以轂渦匯入葉根渦中沿槳軸發(fā)展,傳統(tǒng)意義上的轂渦計(jì)算結(jié)果見下面自航工況計(jì)算結(jié)果。梢渦以等螺距螺旋管形式向下游發(fā)展,由于螺旋槳對(duì)流動(dòng)的抽吸加速作用,螺旋管所在的柱面直徑略有減小。葉根渦的尺度較分散,但也以纏繞槳軸的螺旋管形式向下游發(fā)展。梢渦強(qiáng)度最高,形式最顯著。

    對(duì)于自航工況,采用與試驗(yàn)一致的等航速變轉(zhuǎn)速強(qiáng)迫自航方法進(jìn)行數(shù)值模擬。計(jì)算速度V選取為2.8~4.2 m/s,間隔為0.2 m/s。在每一自航速度下改變轉(zhuǎn)速n(4~5 點(diǎn)),每次記錄轉(zhuǎn)速n、強(qiáng)制力z、槳推力TB、艇阻力R和槳扭矩QB,結(jié)合敞水?dāng)?shù)據(jù)分析自航因子,并對(duì)各個(gè)速度下的自航因子進(jìn)行算術(shù)平均作為最終結(jié)果。通過自航計(jì)算可以得到如下的艇后螺旋槳水動(dòng)力特征值。

    用等推力法按照艇后推力系數(shù)從敞水曲線上查得J0、KQ0,進(jìn)而計(jì)算得到如下自航因子:

    潛艇帶槳自航狀態(tài)尾部渦旋結(jié)構(gòu)計(jì)算結(jié)果見圖4。自航因子計(jì)算結(jié)果與試驗(yàn)結(jié)果的對(duì)比見圖5,圖中試驗(yàn)結(jié)果為美國UM給出的拖曳水池測試結(jié)果。

    圖4 潛艇自航工況尾部渦旋結(jié)構(gòu)計(jì)算結(jié)果Fig.4 Computed vortical structures around the stern of submarine in self-propulsion condition

    圖5 自航因子計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.5 Comparison between computed results and measurements of self-propulsion factors

    從圖4 中可以看到潛艇自航工況下尾翼與螺旋槳周圍渦旋結(jié)構(gòu)的發(fā)展演化情況。圖4 中左圖為潛艇自航工況尾部渦旋結(jié)構(gòu)的起始階段,此時(shí)尾翼根部的馬蹄渦剛出現(xiàn)不久,渦腿尚未充分發(fā)展,螺旋槳梢渦、轂渦都已出現(xiàn),也在發(fā)展演化之中;右圖為潛艇自航工況尾部渦旋結(jié)構(gòu)的充分發(fā)展階段,經(jīng)過30 個(gè)旋轉(zhuǎn)周期渦旋流場的發(fā)展演化后,尾翼根部馬蹄渦、螺旋槳梢渦、轂渦、葉根渦都已達(dá)到成熟的形態(tài)。與敞水工況相比,螺旋槳周圍同樣存在明顯的梢渦、葉根渦、轂渦,而且此時(shí)螺旋槳在潛艇尾部與尾翼不均勻來流相互作用,使得流動(dòng)形式更為復(fù)雜。梢渦以等螺距螺旋管形式向下游發(fā)展,葉根渦非常明顯,以離散螺旋形式向下游泄落,轂渦如同兩股交替纏繞的水流相互作用,為雙螺旋結(jié)構(gòu),穩(wěn)定地向下游發(fā)展。梢渦強(qiáng)度最高,葉根渦與轂渦的形式也非常顯著。

    從圖5 可以看出,與試驗(yàn)結(jié)果相比,推力減額的計(jì)算誤差為2.8%~11.7%,伴流分?jǐn)?shù)的計(jì)算誤差為4.6%~9.5%,相對(duì)旋轉(zhuǎn)效率的計(jì)算誤差為1.2%~2.3%??梢姶鬁u模擬方法對(duì)自航因子與艇后螺旋槳水動(dòng)力(推力、扭矩)的預(yù)報(bào)達(dá)到了工程預(yù)報(bào)精度,計(jì)算效果是令人滿意的。

    3.2 螺旋槳噪聲計(jì)算結(jié)果分析

    朱錫清(2008)[26]對(duì)螺旋槳噪聲的成因與種類進(jìn)行了分析和闡釋。當(dāng)螺旋槳工作在非空泡狀態(tài),螺旋槳的噪聲主要由離散譜(線譜)噪聲、低頻寬帶噪聲和中高頻隨邊噪聲組成。一般認(rèn)為,離散譜噪聲主要是由于螺旋槳工作在船艉的非均勻流場中,當(dāng)螺旋槳葉片周期性旋轉(zhuǎn)時(shí),會(huì)和這非均勻流場相互作用產(chǎn)生非定常升力脈動(dòng),從而輻射出周期性的離散譜噪聲。螺旋槳的低頻寬帶噪聲是由于槳工作在船艉和附體形成的厚湍流邊界層中,當(dāng)這隨機(jī)的湍流和螺旋槳葉片相互作用時(shí)就會(huì)引起葉片攻角的脈動(dòng),即形成沿葉片葉展方向的升力脈動(dòng),從而產(chǎn)生了低頻寬帶噪聲。螺旋槳中高頻寬帶噪聲主要是由葉片上產(chǎn)生的湍流邊界層對(duì)流通過隨邊時(shí)發(fā)生散射引起的,因此也稱為隨邊噪聲。

    本文采用Powell 渦聲方程計(jì)算得到了AU5-65 螺旋槳敞水工況噪聲與自航工況噪聲。圖6 給出了兩個(gè)工況的典型計(jì)算結(jié)果,此時(shí)兩種工況下來流水速都為4 m/s,槳軸轉(zhuǎn)速為25 r/s。聲學(xué)積分區(qū)域?yàn)橐粓A柱體,上游在槳盤面之前1倍槳徑之處,下游在槳盤面之后5倍槳徑之處,圓柱體直徑為2倍槳徑。在敞水與自航工況下螺旋槳旋轉(zhuǎn)區(qū)域形式與網(wǎng)格數(shù)量完全相同,亦即由幾何形式與網(wǎng)格數(shù)量造成的計(jì)算誤差可以忽略。圖6中的螺旋槳噪聲結(jié)果為聲壓譜源級(jí),試驗(yàn)結(jié)果為MSU 給出的1/3 Oct.水槽測試結(jié)果,已經(jīng)扣除了尾翼噪聲,為了展示低頻線譜噪聲,計(jì)算結(jié)果給出的是連續(xù)譜,水槽中的試驗(yàn)一般難以辨識(shí)低頻線譜噪聲。從圖6中計(jì)算結(jié)果可以看出,無論是敞水工況,還是自航工況,在500 Hz以下的低頻段都存在窄帶線譜噪聲,自航工況的線譜噪聲更為明顯,這主要是由SUBOFF 潛艇尾部帶十字型尾翼造成的比較強(qiáng)烈的非均勻流場所致。由于軸頻為25 Hz,螺旋槳為5葉槳,所以理論上的1階葉頻為125 Hz,2階葉頻為250 Hz,3階葉頻為375 Hz。數(shù)值計(jì)算得到的1階葉頻為125.8 Hz,2階葉頻為251.7 Hz,3 階葉頻為375.6 Hz,與理論分析結(jié)果非常吻合。在1 階葉頻處,螺旋槳自航工況噪聲比敞水工況噪聲增大23.6 dB;在2 階葉頻處,自航工況噪聲比敞水工況噪聲增大19.9 dB;在3 階葉頻處,自航工況噪聲比敞水工況噪聲增大18.4 dB。可見艇尾的三維非均勻入流對(duì)低頻線譜噪聲有明顯影響。對(duì)于敞水工況,1階線譜噪聲比2階大5.2 dB,2階線譜噪聲比3階大2.4 dB;對(duì)于自航工況,1階線譜噪聲比2階大8.9 dB,2階線譜噪聲比3階大3.9 dB。計(jì)算結(jié)果從定性與定量兩方面來看也都與理論分析結(jié)果比較吻合。從圖6中的計(jì)算結(jié)果還可以看到,除了線譜噪聲之外,還存在隨頻率逐漸衰減的寬帶噪聲,試驗(yàn)結(jié)果表明0.5~10 kHz自航工況噪聲比敞水工況增大5.2~16.1 dB,數(shù)值計(jì)算所反映的差異與試驗(yàn)基本一致。在0.5~10 kHz 頻段,螺旋槳噪聲數(shù)值計(jì)算結(jié)果與試驗(yàn)差異為2.5~8.9 dB,從聲壓譜譜型和幅值來看,計(jì)算結(jié)果令人滿意,可以滿足工程設(shè)計(jì)中預(yù)報(bào)和優(yōu)化選型互比的需求。

    圖6 螺旋槳噪聲計(jì)算結(jié)果與試驗(yàn)對(duì)比Fig.6 Comparison of computed propeller noise with measurement

    4 結(jié) 論

    本文主要基于LES 結(jié)合Powell 渦聲理論對(duì)于螺旋槳水動(dòng)力與噪聲數(shù)值預(yù)報(bào)方法進(jìn)行了研究,闡釋了Powell 渦聲方程的物理內(nèi)涵,給出了遠(yuǎn)場解的數(shù)學(xué)表達(dá),對(duì)水動(dòng)力與噪聲計(jì)算結(jié)果進(jìn)行了分析,并與試驗(yàn)結(jié)果進(jìn)行了對(duì)比驗(yàn)證,證明了計(jì)算方法切實(shí)可行,計(jì)算結(jié)果可靠。將來可將此方法在螺旋槳工程初步設(shè)計(jì)與優(yōu)化設(shè)計(jì)中加以應(yīng)用,結(jié)合超級(jí)計(jì)算機(jī)并行計(jì)算技術(shù),以減少試驗(yàn)巨大的工作量,為實(shí)用設(shè)計(jì)服務(wù)。本文得到的主要結(jié)論如下:

    (1)基于LES 與Powell渦聲理論的艇槳耦合狀態(tài)螺旋槳水動(dòng)力與噪聲數(shù)值預(yù)報(bào)方法,計(jì)算過程穩(wěn)定,計(jì)算效果好,具有工程實(shí)用價(jià)值,可在螺旋槳初步設(shè)計(jì)階段以及優(yōu)化階段加以運(yùn)用。該方法能同時(shí)求解出水動(dòng)力和噪聲,可為艦艇航速預(yù)報(bào)與噪聲評(píng)估服務(wù)。

    (2)不同進(jìn)速系數(shù)下,與試驗(yàn)結(jié)果相比,敞水工況推力系數(shù)的計(jì)算誤差為2%~3%,扭矩系數(shù)的計(jì)算誤差為3%~4%,敞水效率的計(jì)算誤差為1%~2%;自航工況推力減額的計(jì)算誤差為2.8%~11.7%,伴流分?jǐn)?shù)的計(jì)算誤差為4.6%~9.5%,相對(duì)旋轉(zhuǎn)效率的計(jì)算誤差為1.2%~2.3%;在0.5~10 kHz頻段,螺旋槳噪聲數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果差異為2.5~8.9 dB。

    (3)結(jié)合前期研究發(fā)現(xiàn),使用上述計(jì)算方法時(shí),流動(dòng)渦旋結(jié)構(gòu)的準(zhǔn)確求解是正確計(jì)算流激噪聲的關(guān)鍵,否則會(huì)產(chǎn)生較明顯的虛假噪聲。合理的亞格子渦模型與網(wǎng)格數(shù)量要經(jīng)過系統(tǒng)的收斂性研究得到,對(duì)積分區(qū)域也要進(jìn)行收斂性研究。流動(dòng)計(jì)算時(shí)間一般要在30 個(gè)旋轉(zhuǎn)周期以上,然后再開始噪聲計(jì)算,噪聲計(jì)算時(shí)間一般不小于30個(gè)旋轉(zhuǎn)周期。

    致謝:作者在螺旋槳噪聲計(jì)算研究中得到了中國船舶科學(xué)研究中心朱錫清研究員的指教,在數(shù)值求解渦聲方程方面得到了張效慈研究員的指教,在此向兩位前輩表示衷心感謝!

    猜你喜歡
    螺旋槳聲學(xué)潛艇
    十分鐘讀懂潛艇史(下)
    潛艇哥別撞我
    十分鐘讀懂潛艇史(上)
    潛艇躍進(jìn)之黃金時(shí)代
    愛的就是這股Hi-Fi味 Davis Acoustics(戴維斯聲學(xué))Balthus 70
    基于CFD的螺旋槳拉力確定方法
    Acoustical Treatment Primer:Diffusion談?wù)劼晫W(xué)處理中的“擴(kuò)散”
    Acoustical Treatment Primer:Absorption談?wù)劼晫W(xué)處理中的“吸聲”(二)
    Acoustical Treatment Primer:Absorption 談?wù)劼晫W(xué)處理中的“吸聲”
    3800DWT加油船螺旋槳諧鳴分析及消除方法
    廣東造船(2015年6期)2015-02-27 10:52:46
    少妇人妻精品综合一区二区| 制服丝袜香蕉在线| 91aial.com中文字幕在线观看| 亚州av有码| 欧美zozozo另类| 麻豆精品久久久久久蜜桃| 黄色欧美视频在线观看| 日本-黄色视频高清免费观看| eeuss影院久久| 丰满乱子伦码专区| 国产黄片视频在线免费观看| 成人鲁丝片一二三区免费| 亚洲精品中文字幕在线视频 | 国产一级毛片在线| 又爽又黄a免费视频| 久久鲁丝午夜福利片| 亚洲精品乱码久久久久久按摩| 爱豆传媒免费全集在线观看| 日本欧美国产在线视频| 一级爰片在线观看| 永久免费av网站大全| 国产精品国产三级国产av玫瑰| av卡一久久| 成人亚洲精品av一区二区| 国产伦理片在线播放av一区| 少妇人妻一区二区三区视频| 免费观看无遮挡的男女| 国产精品久久久久久久久免| 中文资源天堂在线| 哪个播放器可以免费观看大片| 自拍欧美九色日韩亚洲蝌蚪91 | 嫩草影院精品99| 亚洲欧美中文字幕日韩二区| 看黄色毛片网站| 寂寞人妻少妇视频99o| 一级毛片 在线播放| 嫩草影院入口| 丝袜喷水一区| 夜夜看夜夜爽夜夜摸| 亚洲伊人久久精品综合| 禁无遮挡网站| 精品视频人人做人人爽| 婷婷色麻豆天堂久久| 亚洲国产欧美人成| 国产乱来视频区| 赤兔流量卡办理| 国产爱豆传媒在线观看| 国产精品三级大全| 青春草国产在线视频| 中文字幕久久专区| 国产爱豆传媒在线观看| 超碰av人人做人人爽久久| 建设人人有责人人尽责人人享有的 | 干丝袜人妻中文字幕| 欧美老熟妇乱子伦牲交| 亚洲成人中文字幕在线播放| 毛片一级片免费看久久久久| 久久久久久久亚洲中文字幕| 午夜福利在线观看免费完整高清在| 成人国产av品久久久| 人妻制服诱惑在线中文字幕| 97人妻精品一区二区三区麻豆| 久久人人爽av亚洲精品天堂 | 在线精品无人区一区二区三 | 亚洲精品久久午夜乱码| 亚洲久久久久久中文字幕| 欧美激情久久久久久爽电影| 搡女人真爽免费视频火全软件| 少妇猛男粗大的猛烈进出视频 | 国产亚洲91精品色在线| 国产中年淑女户外野战色| 久久久欧美国产精品| 亚洲av国产av综合av卡| 久久精品久久精品一区二区三区| 又爽又黄a免费视频| 高清午夜精品一区二区三区| 97在线视频观看| 麻豆乱淫一区二区| 国产成人精品福利久久| 亚洲国产精品专区欧美| 人妻少妇偷人精品九色| 91狼人影院| 欧美日韩在线观看h| 午夜亚洲福利在线播放| av线在线观看网站| 制服丝袜香蕉在线| 深爱激情五月婷婷| 亚洲图色成人| 永久免费av网站大全| 久久精品综合一区二区三区| 成人二区视频| 国产欧美亚洲国产| 亚洲精品色激情综合| 好男人在线观看高清免费视频| 成人一区二区视频在线观看| 91aial.com中文字幕在线观看| av国产精品久久久久影院| 欧美高清性xxxxhd video| 国产白丝娇喘喷水9色精品| 26uuu在线亚洲综合色| 国产精品久久久久久av不卡| 美女被艹到高潮喷水动态| 91午夜精品亚洲一区二区三区| 91久久精品国产一区二区三区| 99热这里只有是精品在线观看| 国产精品99久久久久久久久| 九九久久精品国产亚洲av麻豆| 一区二区三区免费毛片| 久久精品熟女亚洲av麻豆精品| 国内精品美女久久久久久| 大片电影免费在线观看免费| 水蜜桃什么品种好| 在线播放无遮挡| 欧美日韩视频精品一区| 最后的刺客免费高清国语| 久久久a久久爽久久v久久| 1000部很黄的大片| 一级二级三级毛片免费看| 干丝袜人妻中文字幕| 亚洲最大成人中文| 少妇人妻 视频| 精品久久久精品久久久| 看黄色毛片网站| 深夜a级毛片| 91久久精品国产一区二区成人| 亚洲国产精品国产精品| 天天躁夜夜躁狠狠久久av| 偷拍熟女少妇极品色| 中文字幕亚洲精品专区| 国产免费福利视频在线观看| 深爱激情五月婷婷| 丝袜喷水一区| 亚洲欧美精品专区久久| 卡戴珊不雅视频在线播放| 国产综合懂色| 一区二区三区精品91| av在线播放精品| 亚洲图色成人| 亚洲av国产av综合av卡| 97热精品久久久久久| 中文字幕久久专区| 3wmmmm亚洲av在线观看| 久久久色成人| 亚洲精品日韩av片在线观看| 热99国产精品久久久久久7| 午夜视频国产福利| 免费观看a级毛片全部| 18禁在线无遮挡免费观看视频| 免费黄频网站在线观看国产| 一区二区三区精品91| 国产免费视频播放在线视频| 国产淫语在线视频| 亚洲成人中文字幕在线播放| 男人爽女人下面视频在线观看| 国产伦理片在线播放av一区| 亚洲国产av新网站| 久久久久国产精品人妻一区二区| 伦精品一区二区三区| 一区二区av电影网| 97精品久久久久久久久久精品| 搞女人的毛片| 精品午夜福利在线看| 99热6这里只有精品| 国产欧美亚洲国产| 久久久久久久午夜电影| 特大巨黑吊av在线直播| 国产乱人视频| 美女xxoo啪啪120秒动态图| 成年av动漫网址| 亚洲欧美精品自产自拍| 免费观看无遮挡的男女| 另类亚洲欧美激情| 久久久久精品久久久久真实原创| 狂野欧美激情性bbbbbb| 搞女人的毛片| 男女边吃奶边做爰视频| 中文字幕免费在线视频6| 国产成人精品久久久久久| 嫩草影院精品99| 久久精品国产a三级三级三级| av在线蜜桃| 成人漫画全彩无遮挡| 久久精品国产鲁丝片午夜精品| 欧美xxⅹ黑人| 啦啦啦中文免费视频观看日本| 老女人水多毛片| 中文天堂在线官网| 六月丁香七月| 人妻一区二区av| 亚洲综合色惰| 香蕉精品网在线| 亚洲电影在线观看av| 舔av片在线| 国产欧美日韩一区二区三区在线 | 国产 精品1| 99re6热这里在线精品视频| 免费av毛片视频| 少妇人妻精品综合一区二区| 亚洲电影在线观看av| 噜噜噜噜噜久久久久久91| 特大巨黑吊av在线直播| 日韩av不卡免费在线播放| 一二三四中文在线观看免费高清| 可以在线观看毛片的网站| 亚洲成人中文字幕在线播放| 亚洲精品成人av观看孕妇| 国产白丝娇喘喷水9色精品| 九色成人免费人妻av| 欧美极品一区二区三区四区| 国产极品天堂在线| 久久久色成人| 久久亚洲国产成人精品v| 亚洲av男天堂| 欧美丝袜亚洲另类| 欧美三级亚洲精品| av在线播放精品| 成人一区二区视频在线观看| 久久久久精品久久久久真实原创| 国产亚洲一区二区精品| 在线天堂最新版资源| 日韩强制内射视频| 亚洲在线观看片| 男人舔奶头视频| 美女内射精品一级片tv| 亚洲人成网站在线观看播放| 亚洲av福利一区| 嘟嘟电影网在线观看| 国产午夜福利久久久久久| 一级毛片久久久久久久久女| 亚洲av中文字字幕乱码综合| 国产成人免费观看mmmm| 特大巨黑吊av在线直播| 亚洲av国产av综合av卡| 国产精品久久久久久精品电影小说 | 久久久久久久久久久免费av| 国产亚洲最大av| av在线app专区| 免费高清在线观看视频在线观看| 国产淫片久久久久久久久| 成人高潮视频无遮挡免费网站| 国产av不卡久久| 中文字幕制服av| 中文天堂在线官网| 99久国产av精品国产电影| 伦精品一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91 | 韩国高清视频一区二区三区| 边亲边吃奶的免费视频| 国产亚洲精品久久久com| 伊人久久精品亚洲午夜| 亚洲最大成人av| 国产综合懂色| 精品久久久久久电影网| 久热久热在线精品观看| 日产精品乱码卡一卡2卡三| 日韩免费高清中文字幕av| 国产精品99久久久久久久久| 舔av片在线| 欧美一级a爱片免费观看看| 少妇丰满av| 国产乱来视频区| 身体一侧抽搐| 男插女下体视频免费在线播放| 国产精品99久久久久久久久| 青春草亚洲视频在线观看| 亚洲av成人精品一区久久| 久久精品久久久久久噜噜老黄| 最近的中文字幕免费完整| 欧美区成人在线视频| 超碰av人人做人人爽久久| 国产一区二区三区综合在线观看 | 精品国产乱码久久久久久小说| 高清欧美精品videossex| 国产精品99久久久久久久久| 一级片'在线观看视频| 亚洲精品456在线播放app| 熟女电影av网| 啦啦啦中文免费视频观看日本| 色综合色国产| 午夜激情久久久久久久| 直男gayav资源| 国产爱豆传媒在线观看| 成人无遮挡网站| 国产精品人妻久久久久久| 一级毛片 在线播放| 久久精品国产亚洲网站| 亚洲国产高清在线一区二区三| 欧美一区二区亚洲| 精品久久久久久久久av| 国产一区有黄有色的免费视频| 亚洲欧美日韩无卡精品| 国产美女午夜福利| 你懂的网址亚洲精品在线观看| 欧美一区二区亚洲| av在线老鸭窝| 少妇人妻 视频| 99热全是精品| 日本-黄色视频高清免费观看| 久久99热这里只频精品6学生| tube8黄色片| 亚洲熟女精品中文字幕| 一本一本综合久久| 一个人看的www免费观看视频| 综合色av麻豆| 国产探花在线观看一区二区| 爱豆传媒免费全集在线观看| 男人狂女人下面高潮的视频| 婷婷色麻豆天堂久久| 美女cb高潮喷水在线观看| 国产毛片在线视频| 一级片'在线观看视频| 国产爽快片一区二区三区| 汤姆久久久久久久影院中文字幕| 国产人妻一区二区三区在| 久久久久久久久久久免费av| 日本爱情动作片www.在线观看| 精品视频人人做人人爽| 爱豆传媒免费全集在线观看| 涩涩av久久男人的天堂| 一区二区av电影网| 国产91av在线免费观看| 久久99精品国语久久久| 久久精品国产亚洲网站| 日韩人妻高清精品专区| 少妇 在线观看| 纵有疾风起免费观看全集完整版| 色婷婷久久久亚洲欧美| 1000部很黄的大片| 国产精品一区www在线观看| av.在线天堂| 美女内射精品一级片tv| 大陆偷拍与自拍| 欧美成人a在线观看| 岛国毛片在线播放| 色视频在线一区二区三区| 日日摸夜夜添夜夜添av毛片| 中国国产av一级| 亚洲av中文字字幕乱码综合| 国产亚洲最大av| 80岁老熟妇乱子伦牲交| 狂野欧美激情性xxxx在线观看| 成人欧美大片| 亚洲国产欧美人成| 免费看不卡的av| 夫妻午夜视频| 国产午夜精品一二区理论片| 亚洲精品国产色婷婷电影| 成人欧美大片| 一边亲一边摸免费视频| 久久这里有精品视频免费| 午夜激情久久久久久久| 一级毛片aaaaaa免费看小| 一级片'在线观看视频| 国产日韩欧美在线精品| 九色成人免费人妻av| 男人添女人高潮全过程视频| 国产男人的电影天堂91| 可以在线观看毛片的网站| 久久久国产一区二区| 欧美97在线视频| 国产男人的电影天堂91| 99久久精品一区二区三区| 麻豆国产97在线/欧美| 国产色婷婷99| 少妇的逼好多水| 联通29元200g的流量卡| 国产精品三级大全| 国产午夜精品一二区理论片| 一级二级三级毛片免费看| 国产伦理片在线播放av一区| 久久99热这里只频精品6学生| 国产探花极品一区二区| 在线观看免费高清a一片| 亚洲色图av天堂| 91精品伊人久久大香线蕉| 精品少妇黑人巨大在线播放| 免费观看无遮挡的男女| 国产男女超爽视频在线观看| 国语对白做爰xxxⅹ性视频网站| 人体艺术视频欧美日本| 在现免费观看毛片| 国产乱人视频| 美女高潮的动态| 2021天堂中文幕一二区在线观| 韩国高清视频一区二区三区| 欧美日韩精品成人综合77777| 亚洲av不卡在线观看| 少妇人妻一区二区三区视频| 亚洲人成网站在线播| 最新中文字幕久久久久| 国产日韩欧美亚洲二区| 国产一区亚洲一区在线观看| 99久久精品一区二区三区| 成人特级av手机在线观看| 午夜视频国产福利| 国产综合精华液| 日韩国内少妇激情av| 免费观看性生交大片5| 欧美激情久久久久久爽电影| 97热精品久久久久久| 全区人妻精品视频| 中文乱码字字幕精品一区二区三区| 交换朋友夫妻互换小说| 女人十人毛片免费观看3o分钟| 蜜桃久久精品国产亚洲av| 国产黄片美女视频| 麻豆国产97在线/欧美| 国产一级毛片在线| 最近中文字幕2019免费版| 综合色丁香网| 欧美一级a爱片免费观看看| av国产免费在线观看| 午夜爱爱视频在线播放| 99精国产麻豆久久婷婷| 国产黄片视频在线免费观看| 亚洲av中文字字幕乱码综合| 一本色道久久久久久精品综合| 婷婷色av中文字幕| 边亲边吃奶的免费视频| 午夜老司机福利剧场| 成人高潮视频无遮挡免费网站| a级毛片免费高清观看在线播放| 国产大屁股一区二区在线视频| 亚洲熟女精品中文字幕| 亚洲欧美一区二区三区国产| 国产有黄有色有爽视频| 最近中文字幕高清免费大全6| 99久久精品一区二区三区| 丰满人妻一区二区三区视频av| 亚洲人成网站高清观看| 久久久精品94久久精品| 亚洲,欧美,日韩| 最近中文字幕高清免费大全6| 国产黄片美女视频| 热99国产精品久久久久久7| 欧美精品人与动牲交sv欧美| 国产视频内射| 久久久久精品久久久久真实原创| 日本一本二区三区精品| 亚洲精品国产色婷婷电影| 在线精品无人区一区二区三 | 美女视频免费永久观看网站| 国产伦精品一区二区三区视频9| 男女那种视频在线观看| 国产视频首页在线观看| 国产成人freesex在线| 美女内射精品一级片tv| 成人毛片60女人毛片免费| 51国产日韩欧美| 国产精品成人在线| 亚洲精品一区蜜桃| 大片电影免费在线观看免费| 国产毛片在线视频| 欧美xxxx黑人xx丫x性爽| 亚洲第一区二区三区不卡| 简卡轻食公司| 亚洲精品日韩av片在线观看| 中文在线观看免费www的网站| 777米奇影视久久| kizo精华| 性插视频无遮挡在线免费观看| 肉色欧美久久久久久久蜜桃 | 亚洲成人中文字幕在线播放| 嘟嘟电影网在线观看| 亚洲欧美日韩无卡精品| 尤物成人国产欧美一区二区三区| 2018国产大陆天天弄谢| 伦精品一区二区三区| 嫩草影院入口| 亚洲av欧美aⅴ国产| 美女高潮的动态| 国产伦精品一区二区三区视频9| 午夜精品一区二区三区免费看| 九色成人免费人妻av| 麻豆成人av视频| 亚洲无线观看免费| 一个人观看的视频www高清免费观看| 熟女电影av网| 秋霞在线观看毛片| 亚洲av成人精品一二三区| 国产成人精品一,二区| 欧美少妇被猛烈插入视频| 亚洲最大成人中文| 亚洲最大成人手机在线| 自拍偷自拍亚洲精品老妇| 国产91av在线免费观看| 国产国拍精品亚洲av在线观看| 国产成人一区二区在线| 一级毛片我不卡| 日韩av不卡免费在线播放| 国产精品蜜桃在线观看| 国产精品人妻久久久影院| 午夜免费男女啪啪视频观看| av卡一久久| 纵有疾风起免费观看全集完整版| 国产一级毛片在线| 亚洲精品色激情综合| 国产精品不卡视频一区二区| 毛片一级片免费看久久久久| 欧美精品人与动牲交sv欧美| 中文字幕制服av| 亚洲国产欧美在线一区| av在线亚洲专区| 能在线免费看毛片的网站| 男女啪啪激烈高潮av片| 欧美高清性xxxxhd video| 高清av免费在线| 十八禁网站网址无遮挡 | 黄色怎么调成土黄色| 高清午夜精品一区二区三区| 国内精品宾馆在线| 亚洲欧美日韩无卡精品| 亚洲av成人精品一区久久| 深爱激情五月婷婷| 日本免费在线观看一区| 日韩一区二区三区影片| 国产成人91sexporn| 另类亚洲欧美激情| 好男人视频免费观看在线| 中文字幕av成人在线电影| 女的被弄到高潮叫床怎么办| 精品午夜福利在线看| 特级一级黄色大片| 婷婷色综合www| freevideosex欧美| 国产日韩欧美亚洲二区| 亚洲av电影在线观看一区二区三区 | 久久久久久久午夜电影| 亚洲经典国产精华液单| 日韩 亚洲 欧美在线| 免费播放大片免费观看视频在线观看| 青春草视频在线免费观看| 欧美日韩视频高清一区二区三区二| 亚洲人成网站在线播| 精品久久久久久久久av| 白带黄色成豆腐渣| 最近的中文字幕免费完整| 精品国产三级普通话版| 秋霞伦理黄片| 性插视频无遮挡在线免费观看| 精品人妻偷拍中文字幕| 亚洲图色成人| 舔av片在线| 在线天堂最新版资源| 热re99久久精品国产66热6| 爱豆传媒免费全集在线观看| 免费少妇av软件| 久久久久久国产a免费观看| www.av在线官网国产| 人妻少妇偷人精品九色| 国产精品av视频在线免费观看| 亚洲av中文字字幕乱码综合| 久久午夜福利片| 国产精品国产三级专区第一集| 尤物成人国产欧美一区二区三区| 亚洲av电影在线观看一区二区三区 | 青春草亚洲视频在线观看| 插逼视频在线观看| 欧美丝袜亚洲另类| 免费不卡的大黄色大毛片视频在线观看| 91精品一卡2卡3卡4卡| 国模一区二区三区四区视频| 国产黄频视频在线观看| 亚洲在线观看片| 亚洲欧美日韩卡通动漫| 日韩一本色道免费dvd| 尤物成人国产欧美一区二区三区| 国产精品女同一区二区软件| 欧美xxⅹ黑人| 成人亚洲精品av一区二区| 久久精品综合一区二区三区| 成人亚洲欧美一区二区av| 国产视频内射| 国产男人的电影天堂91| 成人免费观看视频高清| 少妇人妻精品综合一区二区| 我的女老师完整版在线观看| 高清视频免费观看一区二区| 毛片女人毛片| 日韩中字成人| 午夜免费观看性视频| 97热精品久久久久久| 国产午夜精品久久久久久一区二区三区| 免费少妇av软件| 国产免费视频播放在线视频| 亚洲精品久久久久久婷婷小说| 国产高清三级在线| 一区二区三区四区激情视频| 日韩av在线免费看完整版不卡| 日日啪夜夜爽| 少妇裸体淫交视频免费看高清| 亚洲欧美一区二区三区黑人 | 99热全是精品| 国产欧美另类精品又又久久亚洲欧美| 九草在线视频观看| 午夜激情福利司机影院| 黑人高潮一二区| 九草在线视频观看| 秋霞在线观看毛片| 国产精品人妻久久久影院| 色综合色国产| 国产精品国产三级国产专区5o| 777米奇影视久久| 三级经典国产精品| 久久精品久久久久久久性| 亚洲国产欧美在线一区| 人人妻人人澡人人爽人人夜夜| 欧美丝袜亚洲另类| 高清av免费在线| 99精国产麻豆久久婷婷| 国产v大片淫在线免费观看| 综合色av麻豆| 国内精品美女久久久久久| 国产在线一区二区三区精| 亚洲精品久久午夜乱码|