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

    電磁力控制翼型三維繞流場(chǎng)特性的數(shù)值研究

    2016-05-04 18:39:28尹紀(jì)富林忠義尤云祥胡天群
    船舶力學(xué) 2016年11期
    關(guān)鍵詞:渦量電磁力攻角

    尹紀(jì)富,林忠義,李 巍,尤云祥,胡天群

    (1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240;2.中交疏浚技術(shù)裝備國(guó)家工程研究中心,上海 201208;3.嘉興南陽(yáng)職業(yè)技術(shù)學(xué)院,浙江 嘉興 314003)

    電磁力控制翼型三維繞流場(chǎng)特性的數(shù)值研究

    尹紀(jì)富1,2,林忠義3,李 巍1,尤云祥1,胡天群1

    (1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240;2.中交疏浚技術(shù)裝備國(guó)家工程研究中心,上海 201208;3.嘉興南陽(yáng)職業(yè)技術(shù)學(xué)院,浙江 嘉興 314003)

    在弦長(zhǎng)雷諾數(shù)ReL=2.97×106下,采用脫體渦模擬方法對(duì)弱電解質(zhì)中電磁力作用下翼型繞流場(chǎng)特性進(jìn)行了數(shù)值模擬,研究了電磁力作用控制翼型失速攻角時(shí)繞流場(chǎng)中三維流動(dòng)特性及失速問(wèn)題的規(guī)律和機(jī)理。結(jié)果表明:電磁力作用可有效改善翼型周?chē)牧鲌?chǎng)結(jié)構(gòu),顯著減小翼型繞流場(chǎng)三維特性,并使翼型表面流體動(dòng)能增加,當(dāng)電磁力作用足夠大時(shí),其表面渦量轉(zhuǎn)變?yōu)檎驕u量。同時(shí),電磁力作用可增加翼型升力,減小阻力,顯著減小升阻力脈動(dòng)特性,提高翼型升阻比,還可以顯著延緩翼型失速特性,增加失速攻角,提高工作性能。

    電磁力;翼型繞流;失速;數(shù)值模擬

    0 引 言

    翼型結(jié)構(gòu)在船舶工程中應(yīng)用廣泛,在其隨船舶以一定攻角航行的過(guò)程中,翼型表面的流動(dòng)分離會(huì)使其升力減小而阻力增加,甚至失速,直接影響船舶的推進(jìn)效率和操控性能。失速是翼型攻角增加到某個(gè)臨界值時(shí)升力突然減小的一種狀態(tài),這是一種在船舶航行過(guò)程中不希望出現(xiàn)的現(xiàn)象[1]。因此如何有效抑制和延緩翼型結(jié)構(gòu)的失速問(wèn)題,在船舶工程中具有十分重要的工程價(jià)值和學(xué)術(shù)研究意義。

    翼型表面流動(dòng)分離是產(chǎn)生失速現(xiàn)象的本質(zhì)原因之一,因此有效地抑制翼型表面的流動(dòng)分離是抑制和延緩翼型失速現(xiàn)象的重要途徑。目前對(duì)翼型結(jié)構(gòu)表面流動(dòng)分離的控制方法一般分為被動(dòng)控制和主動(dòng)控制兩種。被動(dòng)控制方法主要有采用柔壁型結(jié)構(gòu)、安裝端板、波浪形表面結(jié)構(gòu)、尾刷、導(dǎo)流板和襟翼結(jié)構(gòu)等方法[2-5]。由于被動(dòng)控制方法屬于預(yù)先設(shè)定的附加裝置,當(dāng)實(shí)際流動(dòng)狀態(tài)與設(shè)計(jì)狀態(tài)不符時(shí),難以達(dá)到預(yù)想的流動(dòng)控制效果。主動(dòng)控制方法主要有合成射流技術(shù)、吸吹氣方法、等離子體技術(shù)以及電磁力控制,主動(dòng)控制方法可根據(jù)流動(dòng)狀態(tài)調(diào)整控制參數(shù),達(dá)到最佳的流動(dòng)控制效果[6-10]。

    海水作為天然的弱電解質(zhì)給電磁力主動(dòng)控制在船舶工程中的應(yīng)用提供了先天的優(yōu)勢(shì),而且電磁力控制具有響應(yīng)速度快、無(wú)動(dòng)作部件以及易操作等特點(diǎn),因此電磁力控制方式在船舶工程中具有寬廣的應(yīng)用前景。Gailitis等[10]最早提出了電磁力主動(dòng)控制方法,將相間分布的電極與磁極的激活板放置于弱電解質(zhì)溶液,弱電解質(zhì)中產(chǎn)生的Lorentz力作用改變了平板邊界層的流動(dòng)。上世紀(jì)末期,對(duì)此主動(dòng)控制方法的研究逐漸開(kāi)始趨向于機(jī)理和應(yīng)用層面,研究者們開(kāi)始研究電磁力對(duì)平板、道槽及圓柱繞流等流動(dòng)的影響[11-14]。

    近年來(lái),利用電磁力控制翼型表面流動(dòng)分離,進(jìn)而延緩或抑制翼型失速現(xiàn)象的問(wèn)題,成為國(guó)際上研究的熱點(diǎn)。目前研究成果表明:電磁力作用增加了翼型近壁面流動(dòng)速度,能夠提高翼型最大失速攻角,在減阻方面也有較為顯著的作用,電磁力作用可有效抑制和延緩翼型失速現(xiàn)象的發(fā)生,提高翼型的舵效。當(dāng)電磁力作用與流動(dòng)方向相反時(shí),翼型上表面形成較顯著的渦街結(jié)構(gòu)[15-21]。在此基礎(chǔ)上,研究者還對(duì)流向的洛倫茲力控制翼型繞流場(chǎng)的機(jī)理進(jìn)行了深入的研究,且對(duì)洛倫茲力的控制效率的問(wèn)題進(jìn)行了討論[22]。

    由于受到實(shí)驗(yàn)條件的限制,目前的研究主要集中于翼型結(jié)構(gòu)表面層流邊界層分離時(shí)的電磁力控制問(wèn)題,大多采用二維電磁力的近似公式對(duì)亞臨界區(qū)低雷諾數(shù)時(shí)二維翼型結(jié)構(gòu)的繞流場(chǎng)進(jìn)行控制,對(duì)翼型結(jié)構(gòu)失速后,電磁力控制繞流場(chǎng)中顯著的三維特性和失速問(wèn)題的控制規(guī)律和機(jī)理尚不十分清楚。

    隨著計(jì)算流體力學(xué)CFD技術(shù)的發(fā)展,數(shù)值模擬已成為研究電磁力流動(dòng)主動(dòng)控制的重要手段之一,其中針對(duì)高雷諾數(shù)的湍流問(wèn)題主要采用RANS、LES和DES等方法[23]。RANS方法即雷諾平均方法,該方法模擬湍流脈動(dòng)的平均流動(dòng)作用,忽略了小尺度脈動(dòng)信息,因此對(duì)復(fù)雜的非定常湍流模擬存在困難,對(duì)翼型結(jié)構(gòu)附面層內(nèi)大范圍流動(dòng)分離等現(xiàn)象的模擬結(jié)果并不十分精確。DNS方法即直接數(shù)值模擬方法,直接求解流場(chǎng)中的所有湍流尺度,但由于計(jì)算機(jī)能力的限制,目前還只能研究相對(duì)較低雷諾數(shù)的湍流流動(dòng)問(wèn)題。LES方法即大渦模擬方法,是一種介于RANS和DNS之間的湍流數(shù)值模擬方法,該方法也可以獲得交警錫的真實(shí)瞬時(shí)流場(chǎng),但在近壁面附近需要足夠精細(xì)的網(wǎng)格,網(wǎng)格數(shù)量和質(zhì)量要求較高,目前在實(shí)際工程中很難廣泛應(yīng)用。DES方法是一種介于LES和RANS之間的方法,綜合了兩者之間對(duì)流場(chǎng)數(shù)值模擬的優(yōu)勢(shì),能夠很好地處理高雷諾數(shù)時(shí)附面層內(nèi)的大范圍非穩(wěn)態(tài)流動(dòng)分離,不僅可以提高計(jì)算效率,還能獲得比RANS方法更加精確的流動(dòng)信息[24-25]。

    有鑒于此,本文利用Maxwell方程直接數(shù)值分析表面包覆電極和磁極的翼型結(jié)構(gòu)表面的三維電磁力,將其加入到流動(dòng)控制方程中,在弦長(zhǎng)雷諾數(shù)ReL=2.97×106時(shí),采用DES數(shù)值模擬方法數(shù)值模擬電磁力作用下翼型失速攻角時(shí)繞流場(chǎng)特征,重點(diǎn)分析電磁力控制翼型失速后三維流場(chǎng)特性的規(guī)律及失速問(wèn)題的機(jī)理。

    1 數(shù)值模擬方法

    圖1 表面覆蓋電極與磁極翼型模型Fig.1 Configuration of the hydrofoil covered with electrodes and permanent magnets

    在翼型上表面布置相間分布的電極和磁極,將其放置于弱電解質(zhì)的海水中,對(duì)電極施加電壓后,電極產(chǎn)生的電場(chǎng)與永磁體產(chǎn)生的磁場(chǎng)相互作用,在海水中形成流向的電磁力作用,如圖1所示。其中:u∞為來(lái)流速度,F(xiàn)為電磁力。

    將電磁力作為體積力加入到動(dòng)量方程中,可得不可壓縮流動(dòng)Navier-Stokes方程為

    式中:t為時(shí)間;ρ為電解質(zhì)密度;ν為運(yùn)動(dòng)黏度;u=(u1,u2,u3)為流場(chǎng)瞬時(shí)速度;p為擾動(dòng)壓力;F=(f1,f2,f3)為洛倫茲力;i,j為坐標(biāo)軸方向。

    流動(dòng)的弱電解質(zhì)在電磁場(chǎng)中受Lorentz力的作用,其表達(dá)式為

    式中:J為電流密度,B為磁感應(yīng)強(qiáng)度。電流密度J由兩部分組成,一部分為電極產(chǎn)生的電流,另一部分為流體以速度u在磁場(chǎng)中運(yùn)動(dòng)產(chǎn)生的感應(yīng)電流,可以表示為

    式中:E為電場(chǎng)強(qiáng)度,σ為弱電解質(zhì)電導(dǎo)率。在弱電解海水中,電導(dǎo)率σ很小,因此即便磁場(chǎng)強(qiáng)度很強(qiáng),感應(yīng)電流也很小,將Lorentz力表達(dá)式簡(jiǎn)化為

    采用DES方法對(duì)流動(dòng)控制方程進(jìn)行數(shù)值模擬,DES方法主要思想是在近壁面處采用RANS方法模擬小尺度的脈動(dòng)運(yùn)動(dòng),而在遠(yuǎn)離壁面區(qū)域采用LES方法模擬大尺度脫體分離渦。在翼型壁面邊界層內(nèi),將控制方程時(shí)間平均,采用基于SST k-ω湍流模型的RANS方法對(duì)邊界層內(nèi)流場(chǎng)進(jìn)行數(shù)值模擬。在遠(yuǎn)離翼型壁面區(qū)域,將控制方程空間濾波,采用基于Smagorinsky亞格子模型的LES方法對(duì)附面層外的流場(chǎng)進(jìn)行模擬[26]。

    在SST k-ω湍流模型中k方程耗散項(xiàng)中的湍流尺度參數(shù)lk-ω的表達(dá)式為

    式中:βk為模型常數(shù),取值為βk=0.09。

    在DES方法中,RANS和LES的分辨尺度lDES由下式定義:

    其中,為了防止過(guò)早地開(kāi)啟亞格子應(yīng)力模型引起附面層內(nèi)雷諾應(yīng)力的急劇下降,導(dǎo)致非物理性質(zhì)的分離,產(chǎn)生網(wǎng)格誘導(dǎo)分離現(xiàn)象,系數(shù)CDES=0.61,Δ是網(wǎng)格尺度,對(duì)于非均勻網(wǎng)格有

    在緊貼物面的區(qū)域,由于lk-ω≤CDESΔ,因此lDES=lk-ω,這時(shí)在DES中的湍流模型即為RANS中的SST k-ω湍流模型。在其它區(qū)域,由于lk-ω>CDESΔ,因此SST k-ω湍流模型生成源項(xiàng)和消散源項(xiàng)中的lk-ω被替換為CDESΔ,當(dāng)生成與消散作用達(dá)到平衡時(shí),SST k-ω湍流模型中的湍動(dòng)黏度與Smagorinsky模型中的亞格子渦黏系數(shù)在形式相似,這時(shí)DES方法與LES方法相似。

    記yd為中截面上離開(kāi)翼型壁面的垂向距離,電極磁極寬度a=0.01 m,由于流動(dòng)分離主要發(fā)生在翼型的背壓面,因此電極和磁極布置在翼型背壓面上。利用Maxwell方程,采用有限元方法,對(duì)表面包覆電極與磁極的翼型表面的電磁力分布進(jìn)行直接數(shù)值求解,結(jié)果如圖2所示,其中(a)為法向電磁力密度為某展向斷面處等勢(shì)線(xiàn)分布。x為翼型弦向,y為翼型法向?yàn)殡姶帕?。結(jié)果表明,電磁力在法向方向上呈指數(shù)衰減的現(xiàn)象,在翼型法向約0.015倍的弦長(zhǎng)時(shí)已基本為零,沿展向呈周期性變化特性,最大值出現(xiàn)在電磁極交界處,其變化規(guī)律與文獻(xiàn)[26]中圓柱表面的電磁力分布特征相類(lèi)似。

    圖2 翼型表面電磁力法向(a)和展向(b)分布的數(shù)值結(jié)果Fig.2 Three dimensional distribution of the Lorentz force around the hydrofoil (a)Normal direction(b)Spanwise direction

    圖3 計(jì)算區(qū)域平面示意圖Fig.3 Schematic of a plane of the computational domain

    圖4 數(shù)值計(jì)算平面及三維網(wǎng)格Fig.4 Mesh for numerical computations

    本文針對(duì)NACA0018翼型進(jìn)行研究,弦長(zhǎng)c=0.3 m,攻角在5°~30°之間。數(shù)值模擬時(shí),展向長(zhǎng)度對(duì)翼型繞流場(chǎng)特性有明顯的影響,因此為模擬翼型繞流場(chǎng)中完整的三維特性,翼型展向長(zhǎng)度H選為5c。以翼型前緣起點(diǎn)為坐標(biāo)原點(diǎn),來(lái)流上游為12.5c,下游為21c,數(shù)值模擬計(jì)算平面區(qū)域如圖3所示。入口邊界為速度入口,出口邊界為自由出流,上下邊界為壁面邊界條件,三維方向上為對(duì)稱(chēng)邊界,翼型上下表面為無(wú)滑移壁面邊界條件,以進(jìn)口邊界為計(jì)算初始條件。

    數(shù)值模擬網(wǎng)格劃分采用C型和H型網(wǎng)格對(duì)流場(chǎng)區(qū)域進(jìn)行網(wǎng)格劃分,翼型上下表面等分240份網(wǎng)格,展向在每個(gè)弦長(zhǎng)高度等分30層網(wǎng)格,總網(wǎng)格為2.38×106。數(shù)值網(wǎng)格的平面及三維視圖如圖4所示。為了保證翼型壁面處第一層網(wǎng)格滿(mǎn)足y+≈1,采用下式估算第一層網(wǎng)格控制點(diǎn)離開(kāi)翼型壁面的距離Δyˉ:

    式中:弦長(zhǎng)雷諾數(shù)ReL=u∞c/ν。采用有限體積法離散動(dòng)量和連續(xù)性方程,對(duì)流項(xiàng)離散格式采用QUICK格式,擴(kuò)散項(xiàng)采用中心差分格式,壓力速度耦合迭代采用PISO算法,計(jì)算初始時(shí)間步長(zhǎng)為Δt=0.001 s。

    2 結(jié)果與分析

    對(duì)翼型繞流場(chǎng)進(jìn)行三維數(shù)值模擬時(shí),定義翼型總阻力系數(shù)Cd、總升力系數(shù)Cl及電磁力作用參數(shù)N分別為

    式中:Fd、Fl分別為翼型總阻力和升力,u∞為來(lái)流速度,S為浸濕面積;j0=σE0,E0為平均電場(chǎng)強(qiáng)度,B0為平均磁場(chǎng)強(qiáng)度。

    2.1 無(wú)電磁力作用時(shí)數(shù)值結(jié)果

    在數(shù)值模擬時(shí)弦長(zhǎng)雷諾數(shù)ReL=2.97×106,弱電解質(zhì)密度ρ=1 025 kg/m3,流體導(dǎo)電率為σ=5 s·m-1,運(yùn)動(dòng)黏度ν=10-6m2/s。首先對(duì)沒(méi)有電磁力作用下翼型不同攻角時(shí)的繞流場(chǎng)特性進(jìn)行數(shù)值模擬。當(dāng)數(shù)值模擬時(shí)間t=20 s時(shí),此時(shí)翼型繞流場(chǎng)已經(jīng)達(dá)到穩(wěn)定狀態(tài),然后再對(duì)不同攻角時(shí)的翼型繞流場(chǎng)施加電磁力作用,電磁力作用時(shí)間為5 s,計(jì)算結(jié)束時(shí)刻為t=25 s。

    圖5 當(dāng)t=20 s時(shí),翼型升阻力系數(shù)時(shí)均值數(shù)值與實(shí)驗(yàn)結(jié)果對(duì)比Fig.5 Comparisons of the time-average value of lift/drag coefficients results with the experimental values

    結(jié)果表明:當(dāng)翼型攻角小于失速攻角時(shí),升力系數(shù)隨著攻角的增加而增大,阻力系數(shù)緩慢增大,當(dāng)攻角大于失速攻角時(shí),升力系數(shù)急劇減小,阻力系數(shù)急劇增加。由此可知:翼型失速時(shí),其升力減小,阻力增加,在實(shí)際應(yīng)用中是不期望出現(xiàn)的一種現(xiàn)象,因?yàn)檫@會(huì)導(dǎo)致翼型效率降低,甚至不能夠產(chǎn)生足夠的有效升力,直接影響其穩(wěn)定性和操控性能。

    圖6 當(dāng)t=20 s時(shí),翼型失速前后繞流場(chǎng)渦量場(chǎng)分布 (a)α=12°;(b)α=20°Fig.6 Numerical results for the vorticity distribution at different attack angles(a)α=12°;(b)α=20°

    圖6中給出了翼型失速前和失速后兩個(gè)不同攻角時(shí)渦量場(chǎng)數(shù)值結(jié)果,圖中Ωz為展向渦量。圖中結(jié)果表明,當(dāng)攻角α=12°時(shí),流體流過(guò)翼型時(shí),在翼型背壓面后緣產(chǎn)生小尺度流動(dòng)分離現(xiàn)象,僅在翼型上下表面的尾緣處有狹長(zhǎng)的渦量場(chǎng)分布,此時(shí)翼型繞流場(chǎng)中的渦量場(chǎng)特性顯現(xiàn)為典型層流分離中長(zhǎng)泡分離現(xiàn)象[28],渦量場(chǎng)分布具有層流流動(dòng)特性。當(dāng)攻角α=20°時(shí),流體流經(jīng)翼型時(shí),流體在翼型上表面產(chǎn)生大尺度流動(dòng)分離,流場(chǎng)中有明顯的漩渦產(chǎn)生,翼型繞流場(chǎng)中具有顯著的三維特性。翼型失速后大尺度流動(dòng)分離增加了翼型前后的壓差阻力,升力急劇減小,阻力增加。

    2.2 電磁力作用對(duì)翼型流場(chǎng)特性的影響

    首先考慮電磁力作用對(duì)翼型繞流場(chǎng)流動(dòng)特性的影響規(guī)律。當(dāng)攻角一定時(shí),隨著電磁力作用參數(shù)的改變,翼型繞流場(chǎng)特性隨之改變。圖7和圖8分別給出了ReL=2.97×106下翼型失速后,當(dāng)α=20°時(shí)電磁力作用下流場(chǎng)中三維渦量場(chǎng)及不同展向斷面處渦量變化特性。

    圖中結(jié)果表明,電磁力作用可以抑制渦量場(chǎng)展向大尺度脈動(dòng),隨著電磁力作用的增加,電磁力可抑制翼型表面的流動(dòng)分離,當(dāng)電磁力作用足夠大時(shí),翼型繞流場(chǎng)中渦量場(chǎng)首先轉(zhuǎn)變?yōu)閷恿鞣蛛x泡模型中長(zhǎng)泡分離狀態(tài),隨電磁力作用再增加,在翼型后方電磁力作用推動(dòng)翼型表面的流體流動(dòng),產(chǎn)生類(lèi)似射流狀狹長(zhǎng)渦結(jié)構(gòu)。具體變化規(guī)律如下:

    圖7 當(dāng)α=20°和t=25 s時(shí),電磁力作用下翼型繞流場(chǎng)中渦量場(chǎng)數(shù)值結(jié)果(a)N=0;(b)N=0.1;(c)N=0.25;(d)N=0.5;(e)N=1;(f)N=2Fig.7 Numerical results for the vorticity distribution due to the electromagnetic force when t=25 s and α=20° (a)N=0;(b)N=0.1;(c)N=0.25;(d)N=0.5;(e)N=1;(f)N=2

    圖8 當(dāng)α=20°和t=25 s時(shí),翼型不同展向斷面處渦量場(chǎng)數(shù)值結(jié)果Fig.8 Numerical results for the vorticity distribution of different spanwise sections due to the electromagnetic force when t=25 s and α=20°

    當(dāng)N=0,即沒(méi)有電磁力作用時(shí),翼型后方產(chǎn)生具有顯著三維特性的大尺度渦街結(jié)構(gòu),展向不同斷面處渦量分布有明顯區(qū)別,不同斷面處流動(dòng)分離的位置亦具有顯著差異。當(dāng)N=0.1時(shí),此時(shí)電磁力作用較小,電磁力作用僅改變了翼型前緣處的流動(dòng)狀態(tài),不同展向斷面處渦量分布依舊具有顯著的差異,尾流場(chǎng)中渦街結(jié)構(gòu)的變化并不明顯。當(dāng)N=0.25時(shí),電磁力作用增加,在翼型展向方向上,大尺度渦街結(jié)構(gòu)有所減少,由于展向不同斷面處渦量分布的相位并不相同,在翼型展向中截面附近區(qū)域依舊存在大尺度渦結(jié)構(gòu)。當(dāng)N=0.5時(shí),電磁力作用將翼型表面流動(dòng)分離完全抑制,不同展向斷面處渦量場(chǎng)分布特征逐漸一致,流動(dòng)狀態(tài)類(lèi)似于翼型未失速時(shí)典型層流分離中長(zhǎng)泡分離特征,翼型上表面渦量逐漸增加,開(kāi)始出現(xiàn)正向的渦量。當(dāng)N=1時(shí),翼型上表面渦量由負(fù)渦量轉(zhuǎn)變?yōu)檎郎u量,此時(shí)翼型表面速度場(chǎng)特征發(fā)生明顯變化,在翼型尾部渦量場(chǎng)中出現(xiàn)正負(fù)渦量交叉分布特征,渦量場(chǎng)分布不再由流動(dòng)分離而導(dǎo)致,而是由電磁力推動(dòng)了翼型表面流動(dòng)加速形成狹長(zhǎng)狀渦結(jié)構(gòu)。隨電磁力作用再增加,當(dāng)N=2時(shí),翼型上表面渦量全部轉(zhuǎn)變?yōu)檎郎u量,由電磁力作用產(chǎn)生的渦結(jié)構(gòu)影響范圍逐漸增加,并向流動(dòng)方向的下游發(fā)展。

    圖9中給出了ReL=2.97×106,α=20°時(shí),翼型展向z/c=2.5中截面處上表面速度u隨yd的變化,其中x/c=0和x/c=1.0分別為翼型的前緣和后緣。結(jié)果可知:當(dāng)沒(méi)有電磁力作用時(shí),在翼型展向中截面處翼型前緣速度最大,隨著向翼型后緣移動(dòng),速度值逐漸減小,而且出現(xiàn)速度方向上的變化,由此可知在翼型上表面具有較明顯的速度梯度變化,當(dāng)N=0.1和N=0.25時(shí),此時(shí)電磁力作用較小,由流場(chǎng)結(jié)果可知翼型繞流場(chǎng)中存在顯著的三維特性,此時(shí)展向不同斷面處分離渦脫落的相位存在差異,在展向中截面處的速度場(chǎng)分布并沒(méi)有明顯的規(guī)律性,但當(dāng)yd>0.6c時(shí),不同斷面處的速度基本一致。當(dāng)N=0.5時(shí),翼型上表面各處速度顯著增加,隨著yd的增加,流體速度首先增加然后再逐漸與外流場(chǎng)速度相一致,從翼型前緣到后緣不同垂向斷面處上表面的速度逐漸減小,當(dāng)x/c=0.8斷面時(shí),邊界層內(nèi)的速度與外流場(chǎng)幾乎相當(dāng)。當(dāng)N=1時(shí),翼型上表面各處的速度值繼續(xù)增大,當(dāng)在斷面x/c=1.0時(shí),即在翼型后緣處,其邊界層內(nèi)的速度大小與外流場(chǎng)速度基本達(dá)到一致,當(dāng)yd大于電磁力作用范圍時(shí),其速度值又有所減小,因此電磁力作用使得翼型上表面接近尾緣處的速度梯度方向發(fā)生改變,導(dǎo)致電磁力作用范圍內(nèi)的渦量場(chǎng)由負(fù)渦量轉(zhuǎn)變?yōu)檎郎u量。N=2時(shí),壁面邊界層內(nèi)的速度明顯增大且超過(guò)外流場(chǎng)的速度值,因此當(dāng)yd大于電磁力作用的影響范圍時(shí),其速度場(chǎng)又明顯減小,此時(shí)電磁力作用使上表面正向渦量更加顯著,在流向x方向上的影響范圍亦增加。

    圖9 當(dāng)α=20°和t=25 s時(shí),翼型展向中截面上表面速度u隨yd變化的數(shù)值結(jié)果(a)N=0;(b)N=0.1;(c)N=0.25;(d)N=0.5;(e)N=1;(f)N=2Fig.9 Effect of the electromagnetic force on the tangential velocity near the up-surface of hydrofoil at the middle spanwise section when t=25 s and α=20°(a)N=0;(b)N=0.1;(c)N=0.25;(d)N=0.5;(e)N=1;(f)N=2

    2.3 電磁力作用對(duì)翼型失速特性的影響

    圖10中給出了當(dāng)攻角α=20°時(shí),不同電磁力作用下,升力系數(shù)時(shí)歷及時(shí)均值的數(shù)值結(jié)果。由圖10(a)升力系數(shù)時(shí)歷數(shù)值結(jié)果可知,在無(wú)電磁力作用時(shí),由于翼型繞流場(chǎng)中存在顯著的大尺度分離渦結(jié)構(gòu),其升力系數(shù)時(shí)歷沒(méi)有明顯周期性特征,存在有較大的隨機(jī)性。當(dāng)施加電磁力作用時(shí),隨著電磁力作用的增加,升力系數(shù)隨之增加,且其隨機(jī)脈動(dòng)幅值逐漸減小。當(dāng)N=0.5時(shí),已不存在明顯的脈動(dòng)量,這說(shuō)明此時(shí)翼型繞流場(chǎng)中已沒(méi)有三維特性存在,翼型表面不存在大尺度的流動(dòng)分離現(xiàn)象。由圖10(b)升力系數(shù)時(shí)均值的結(jié)果可知:當(dāng)N≤0.5時(shí),升力系數(shù)時(shí)均值隨著電磁力作用的增加近似線(xiàn)性增加,當(dāng)N>0.5時(shí),升力系數(shù)時(shí)均值隨電磁力作用增加亦近似線(xiàn)性增加,但增加速率有所減緩。

    圖10 當(dāng)α=20°時(shí),電磁力作用對(duì)總升力系數(shù)Cl時(shí)歷及其時(shí)均值的影響 (a)Cl;(b)Fig.10 Effect of the electromagnetic force on the time history and the time average value of the lift coefficients at α=20°(a)Cl;(b)

    圖11 當(dāng)α=20°時(shí),電磁力作用對(duì)總阻力系數(shù)Cd時(shí)歷及其時(shí)均值的影響 (a)Cd;(b)Fig.11 Effect of the electromagnetic force on the time history and the time average value of the drag coefficients at α=20°(a)Cd;(b)

    圖12 當(dāng)α=20°時(shí),電磁力作用對(duì)摩擦阻力系數(shù)時(shí)均值和壓差阻力系數(shù)時(shí)均值的影響 (a);(b)Fig.12 Effect of the electromagnetic force on the total average pressure and friction coefficients at α=20°(a);(b)

    圖11中給出了當(dāng)攻角α=20°時(shí),不同電磁力作用下阻力系數(shù)時(shí)歷及其時(shí)均值的數(shù)值結(jié)果。由圖11(a)阻力系數(shù)時(shí)歷結(jié)果可知:當(dāng)N=0.1時(shí),電磁力作用較小,阻力系數(shù)的時(shí)歷變化并不明顯,隨著電磁力作用的增加阻力系數(shù)減小,其脈動(dòng)幅值也顯著減小。當(dāng)電磁力作用N>0.5時(shí),電磁力作用穩(wěn)定后,其阻力系數(shù)變化并不明顯。由圖11(b)中的時(shí)均值的結(jié)果可知:當(dāng)N≤0.5時(shí),阻力系數(shù)時(shí)均值隨電磁力作用參數(shù)的增加總體上呈現(xiàn)近似線(xiàn)性顯著減小;當(dāng)N>0.5時(shí),阻力系數(shù)的變化趨于平緩,此時(shí)阻力系數(shù)大小基本與翼型未失速時(shí)阻力系數(shù)大小相當(dāng)。

    圖12中給出了當(dāng)攻角α=20°時(shí),電磁力作用對(duì)摩擦阻力系數(shù)和壓差阻力系數(shù)時(shí)均值的影響。圖12(a)中摩擦阻力系數(shù)時(shí)均值數(shù)值結(jié)果可知:摩擦阻力系數(shù)隨電磁力作用的增加在N≤0.5和N>0.5時(shí)均呈現(xiàn)線(xiàn)性增加的趨勢(shì),且當(dāng)N>0.5時(shí),摩擦阻力系數(shù)增加的速率要比N≤0.5時(shí)要小。圖12(b)為壓差阻力系數(shù)時(shí)均值的結(jié)果表明:當(dāng)N≤0.5時(shí),其變化特征與總阻力系數(shù)時(shí)均值變化特征相類(lèi)似;當(dāng)N>0.5時(shí),壓差阻力系數(shù)隨電磁力作用的增加近似線(xiàn)性減小,但減小的速率與摩擦阻力系數(shù)增加的速率相差不大,因此造成了圖11(b)中阻力系數(shù)變化并不明顯的現(xiàn)象。

    升阻比是衡量翼型工作性能的重要標(biāo)志,圖13中給出了當(dāng)攻角α=20°時(shí),電磁力作用對(duì)翼型升阻比的影響特性。結(jié)果可知:翼型升阻比隨著電磁力作用的增加而增大,當(dāng)N≤0.25,電磁力作用較小時(shí),升阻比略有增加但并不顯著;但N=0.5時(shí),翼型升阻比突然增大,隨著電磁力作用再增加,升阻比增加的速率逐漸平緩。

    圖13 當(dāng)α=20°時(shí),電磁力作用升阻比的影響特性Fig.13 Effect of the electromagnetic force on the lift-drag ratio at α=20°

    圖14 電磁力作用下對(duì)升力系數(shù)時(shí)均值及最大升力角的影響Fig.14 Effect of the electromagnetic force on the time-averaged value of the lift coefficients and the angles of maximum lift

    圖15 電磁力作用對(duì)阻力系數(shù)時(shí)均值的影響Fig.15 Effect of the electromagnetic force on the time-averaged of the drag coefficients

    為研究電磁力作用對(duì)翼型失速特性的影響,圖14分別給出了電磁力作用對(duì)翼型升力系數(shù)時(shí)均值及失速攻角的影響。圖14(a)升力系數(shù)時(shí)均值的結(jié)果表明:隨著電磁力作用的增加,翼型升力系數(shù)時(shí)均值增大,電磁力作用下,翼型依舊產(chǎn)生失速現(xiàn)象,翼型失速后,其升力依然急劇減小,但相對(duì)于未受電磁力控制時(shí)依舊增加。圖14(b)中最大失速攻角隨電磁力作用的變化結(jié)果,結(jié)果顯示在給出的電磁力作用范圍內(nèi),最大失速攻角隨著電磁力作用的增加近似線(xiàn)性增加。

    圖15給出了電磁力作用對(duì)阻力系數(shù)時(shí)均值的影響。結(jié)果表明:電磁力作用下,阻力系數(shù)時(shí)均值突然增加的攻角顯著提高,即翼型發(fā)生失速時(shí)的攻角增加,當(dāng)翼型失速后,隨著攻角的增加,阻力系數(shù)時(shí)均值急劇增加,翼型發(fā)生失速現(xiàn)象。

    3 結(jié) 論

    利用Maxwell方程直接數(shù)值計(jì)算表面包覆電極與磁極翼型產(chǎn)生的三維電磁力分布,將其加入到流體控制方程的動(dòng)量方程中,在弦長(zhǎng)雷諾數(shù)ReL=2.97×106時(shí),采用脫體渦模擬方法對(duì)弱電解質(zhì)中電磁力作用下翼型繞流場(chǎng)的流場(chǎng)特性進(jìn)行了數(shù)值模擬,研究了電磁力作用控制失速攻角α=20°時(shí),翼型繞流場(chǎng)中顯著的三維特性的規(guī)律和失速問(wèn)題機(jī)理。

    結(jié)果表明:當(dāng)電磁力作用參數(shù)N<0.25時(shí),電磁力作用對(duì)翼型繞流場(chǎng)的控制效果并不顯著;當(dāng)N= 0.5時(shí),電磁力作用可以抑制翼型表面的大尺度流動(dòng)分離,將翼型繞流場(chǎng)轉(zhuǎn)變?yōu)槲词贂r(shí)的層流分離狀態(tài);當(dāng)N>0.5時(shí),翼型表面流動(dòng)速度顯著增加,速度梯度方向發(fā)生轉(zhuǎn)變,表面渦量顯著增加,渦量由負(fù)渦量轉(zhuǎn)變?yōu)檎郎u量,翼型后方狹長(zhǎng)渦由電磁力作用產(chǎn)生,隨電磁力作用的增加,電磁力作用產(chǎn)生的狹長(zhǎng)渦向下游影響范圍增加。

    當(dāng)電磁力作用參數(shù)N<0.5時(shí),電磁力作用顯著減小翼型升阻力系數(shù)的脈動(dòng)幅值,使得翼型升力系數(shù)增加,阻力系數(shù)減小,升阻比增加,但升阻比增加速率較緩慢;當(dāng)N=0.5時(shí),翼型升力和阻力系數(shù)的脈動(dòng)幅值基本減小到零,翼型升阻比顯著增大;當(dāng)N>0.5時(shí),翼型升阻比增加的速率逐漸減緩。電磁力作用可有效延緩翼型失速問(wèn)題,提高翼型工作性能。

    [1]Ekaterinaris John A,Platzer Max F.Computational prediction of airfoil dynamic stall[J].Aerospace Sci,1997,33:759-846.

    [2]馬 俠,肖林奎,蔣金貴等.用NS方程求解柔壁翼型繞流[J].航空動(dòng)力學(xué)報(bào),1993,8(1):11-18. Ma Xia,Xiao Linkui,Jiang Jingui,et al.Computation of flexible-wall airfoil flow using N-S equations[J].Journal of Aerospace Power,1993,8(1):11-18.

    [3]Tang D,Dowell E H.Aerodyamic flow control of an airfoil with small trailing-edge strips[J].Journal of Aircraft,2006,43 (6):1854-1866.

    [4]郝禮書(shū),喬志德,宋文萍.基于渦流發(fā)生器和Gurney襟翼的翼型繞流流動(dòng)控制試驗(yàn)研究[J].航空學(xué)報(bào),2011,32(8): 1429-1434. Hao Lishu,Qiao Zhide,Song Wenping.Experimental research on control flow over airfoil based on vortex generator and gurney flap[J].Acta Aeronautica ET Astrnautica Sinica,2011,32(8):1429-1434.

    [5]Staelens Y D,Blackwelder R F,Page M A.Novel pitch control effectors for a blended wing body airplane in takeoff and landing configuration[J].AIAA Paper,2007:68-229.

    [6]Hao Feng L,Wang J J,Pan C.Effect of novel synthetic jet on wake vortex shedding modes of a circular cylinder[J].Journal of Fluids and Structures,2010,26(6):900-917.

    [7]Corke T C,Enloe C L,Wilkinson S P.Dielectric barrier discharge plasma actuators for flow control[J].Annual Review of Fluid Mechanics,2010,42:505-529.

    [8]Li Z,Navon I M,Hussaini M Y,et al.Optimal control of cylinder wakes via suction and blowing[J].Computers&Fluids, 2003,32(2):149-171.

    [9]Henoch C,Stace J.Experimental investigation of a salt water turbulent boundary layer modified by an applied streamwise magnetohydrodynamic body force[J].Phys Fluids,1995,7(6):1371-1383.

    [10]Gailitis A,Lielausis O.On a possibility to reduce the hydrodynamic resistance of a plate in a electrolyte[C].Applied Magneto hydrodynamics Report of the Physics Institute Riga 12,1961:143-146.

    [11]Mutschke G,Gerbeth G,Shatrov V,et al.Two-and three-dimensional instabilities of the cylinder wake in an aligned magnetic field[J].Physics of Fluids,1997,9:3114.

    [12]Weier T,Fey U,Gerbeth G,et al.Boundary layer control by means of wall parallel Lorentz forces[J].Magnetohydrodynamics,2001,37(1/2):177-186.

    [13]梅棟杰,范寶春,黃樂(lè)萍等.槽道湍流的展向振蕩電磁力壁面減阻[J].物理學(xué)報(bào),2010,59(10):6786-6792. Mei Dongjie,Fan Baochun,Huang Leping,et al.Drag reduction in turbulent channel flow by spanwise oscillating Lorentz force[J].Acta Phys.Sin,2010,59(10):6786-6792.

    [14]Zhang H,Fan B C,Chen Z H.In-depth study on cylinder wake controlled by Lorentz force[J].Chinese Physics Letters, 2011,28(12):124701.

    [15]Weier T,Gerbeth G.Control of separated flows by time periodic Lorentz forces[J].European Journal of Mechanics B/Fluids,2004,23:835-849.

    [16]Cierpka C,Weier T,Gerbeth G.Synchronized force and particle image velocimetry measurements on a NACA 0015 in poststall under control of time periodic electromagnetic forcing[J].Phys.Fluids,2010,22(7):075109.

    [17]陳耀慧,范寶春,陳志華等.翼型繞流電磁控制的實(shí)驗(yàn)和數(shù)值研究[J].物理學(xué)報(bào),2008,57(2):648-653. Chen Yaohui,Fan Baochun,Chen Zhihua,et al.Experimental and numerical investigations on the electro-magnetic control of hydrofoil wake[J].Acta Phys.Sin,2008,57(2):648-653.

    [18]Chen Y H,Fan B C,Chen Z H,et al.Influences of Lorentz force on the hydrofoil lift[J].Acta Mechanica Sinica,2009, 25(5):589-595.

    [19]劉宗凱,周本謀,劉會(huì)星等.電磁流體表面推進(jìn)機(jī)理與效果分析[J].物理學(xué)報(bào),2011,60(8):84701. Liu Zongkai,Zhou Benmou,Liu Huixing,et al.The analysis of effects and theories for electromagnetic hydrodynamics propulsion by surface[J].Acta Phys.Sin,2011,60(8):84701.

    [20]Albrecht T,Weier T,Gerbeth G,et al.Numerical and experimental investigation of electromagnetic separation control using different wave forms[C]//5th AIAA Flow Control Conference.Chicaco,IL,USA,.2010.

    [21]Mutschke G,Gerbeth G,Albrecht T,el al.Separation control at hydrofoil using Lorentz forces[J].European Journal of Mechanics-B/Fluids,2006,25(2):137-152.

    [22]陳耀慧,董祥瑞,陳志華等.翼型繞流的洛倫茲力控制機(jī)理[J].物理學(xué)報(bào),2014,63(3):034701. Chen Yaohui,Dong Xiangrui,Chen Zhihua,et al.Control of flow around hydrofoil using the Lorentz force[J].Acta Phys. Sin,2014,63(3):034701.

    [23]Breuer M,Jovicˇic'N,Mazaev K.Comparison of DES,RANS and LES for the separated flow around a flat plate at high incidence[J].International Journal for Numerical Methods in Fluids,2003,41(4):357-388.

    [24]STRELETS M.Detached eddy simulation of massively separated flows[R].AIAA,USA,2001:0879.

    [25]Dong L,Men’shov I,Nakamura Y.Numerical prediction for airfoil Stall[J].ICAS2004,Yokohama,Japan,2004.

    [26]尹紀(jì)富,尤云祥,李 巍等.電磁力控制湍流邊界層分離圓柱繞流場(chǎng)特性數(shù)值分析[J].物理學(xué)報(bào),2014,63(4): 044701. Yin Jifu,You Yunxiang,Li Wei,et al.Numerical analysis for the characteristics of flow control around a circular cylinder with a turbulent b oundary layer separation using the electromagnetic force[J].Acta Phys.Sin,2014,63(4):044701.

    [27]Jacobs E N,Sherman A.Airfoil section characteristics as affected by variations of the Reynolds number[M].National Advisory Committee for Aeronautics,1937.

    [28]Horton H P.Laminar separation bubbles in two-and three-dimensional incompressible flow[D].PhDThesis.London:University of London,1968.

    Numerical investigation on the three-dimensional flow characteristics around a hydrofoil influenced by the electromagnetic force

    YIN Ji-fu1,2,LIN Zhong-yi3,LI Wei1,YOU Yun-xiang1,HU Tian-qun1
    (1.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China;2.CCCC National Engineer Research Center of Dredging Technology and Equipment,Shanghai 201208,China;3.School of Jiaxing Nanyang Profession and Technology,Jiaxing 314003,China)

    This paper numerically researches the flow structures around a hydrofoil and the lift/drag characteristics influenced by the electromagnetic force in a weakly conductive fluid at a chord length Reynolds number2.9×106with the method of detached eddy simulation(DES).The effects of the electromagnetic force on the three-dimensional flow characteristics and stall problem of a hydrofoil are investigated at a stall attack angle α=20°.The results show that the electromagnetic force can improve the flow structures around a hydrofoil and significantly reduce the three-dimensional flow characteristics.The electromagnetic force increases the kinetic energy near the boundary layer of the hydrofoil surface,delay the flow separation and change the vorticity field characteristics.Meanwhile,the electromagnetic force can decrease the drag force and the fluctuation amplitudes of the lift,increase the lift,reduce the drag,and raise the lift-drag ratio. Moreover,the electromagnetic force can significantly delay the stall problem and increase the stall angle, improve the working performance.

    electromagnetic force;flow around a hydrofoil;stall;numerical simulation

    O361

    A

    10.3969/j.issn.1007-7294.2016.07.005

    1007-7294(2016)11-1388-12

    2016-09-22

    國(guó)家自然科學(xué)基金(11272211);上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室研究基金資助課題(GP010819)

    尹紀(jì)富(1984-),男,博士研究生,E-mail:yinjifuytu@163.com;尤云祥(1963-),男,教授,博士生導(dǎo)師,通訊作者,E-mail:youyx@sjtu.edu.cn。

    猜你喜歡
    渦量電磁力攻角
    對(duì)真空及電磁波傳播機(jī)理的重新認(rèn)識(shí)
    某型異步感應(yīng)電機(jī)電磁力計(jì)算及分析
    含沙空化對(duì)軸流泵內(nèi)渦量分布的影響
    風(fēng)標(biāo)式攻角傳感器在超聲速飛行運(yùn)載火箭中的應(yīng)用研究
    自由表面渦流動(dòng)現(xiàn)象的數(shù)值模擬
    大攻角狀態(tài)壓氣機(jī)分離流及葉片動(dòng)力響應(yīng)特性
    被動(dòng)電磁裝甲對(duì)金屬射流箍縮電磁力的計(jì)算及驗(yàn)證
    附加攻角效應(yīng)對(duì)顫振穩(wěn)定性能影響
    民用飛機(jī)攻角傳感器安裝定位研究
    航態(tài)對(duì)大型船舶甲板氣流場(chǎng)的影響
    久久99一区二区三区| 日日夜夜操网爽| 在线十欧美十亚洲十日本专区| 999精品在线视频| 日韩欧美一区二区三区在线观看| 免费搜索国产男女视频| 免费不卡黄色视频| 亚洲中文日韩欧美视频| 成人国产一区最新在线观看| 国产成年人精品一区二区 | 欧美最黄视频在线播放免费 | 男女之事视频高清在线观看| 国产成人一区二区三区免费视频网站| 一级毛片女人18水好多| 欧美日韩亚洲高清精品| 中文字幕人妻丝袜一区二区| 香蕉国产在线看| 国产av一区在线观看免费| 老司机深夜福利视频在线观看| 好男人电影高清在线观看| 国产成年人精品一区二区 | 好男人电影高清在线观看| 日本精品一区二区三区蜜桃| 成年人黄色毛片网站| 夫妻午夜视频| 亚洲av五月六月丁香网| 一区福利在线观看| 亚洲黑人精品在线| 欧美成人免费av一区二区三区| 精品国产乱子伦一区二区三区| 亚洲午夜理论影院| www.www免费av| 一区在线观看完整版| 丰满人妻熟妇乱又伦精品不卡| 香蕉久久夜色| 国产又爽黄色视频| 好看av亚洲va欧美ⅴa在| 久久九九热精品免费| 免费高清视频大片| 久久中文字幕一级| 嫩草影视91久久| 久久影院123| 啦啦啦在线免费观看视频4| 欧美黄色片欧美黄色片| 午夜激情av网站| 欧美大码av| 午夜老司机福利片| 久久久久久久久久久久大奶| 一进一出好大好爽视频| 男人舔女人的私密视频| 咕卡用的链子| 午夜久久久在线观看| 国产熟女午夜一区二区三区| 天天添夜夜摸| 久久精品国产亚洲av香蕉五月| 777久久人妻少妇嫩草av网站| 色播在线永久视频| 美女扒开内裤让男人捅视频| 免费不卡黄色视频| 久99久视频精品免费| 天堂中文最新版在线下载| 国产免费av片在线观看野外av| 又紧又爽又黄一区二区| 日韩免费高清中文字幕av| 国产高清国产精品国产三级| 久久这里只有精品19| 免费人成视频x8x8入口观看| 一级a爱片免费观看的视频| 免费在线观看日本一区| 亚洲欧美一区二区三区黑人| 琪琪午夜伦伦电影理论片6080| 少妇被粗大的猛进出69影院| 亚洲精品国产区一区二| 老熟妇乱子伦视频在线观看| 亚洲自偷自拍图片 自拍| 精品久久久久久久久久免费视频 | 极品人妻少妇av视频| 大型av网站在线播放| 80岁老熟妇乱子伦牲交| 亚洲一区二区三区色噜噜 | 9191精品国产免费久久| 热99re8久久精品国产| 亚洲一区二区三区欧美精品| 国产av精品麻豆| 欧美日韩中文字幕国产精品一区二区三区 | 妹子高潮喷水视频| 久久青草综合色| 午夜久久久在线观看| 性少妇av在线| 午夜福利在线观看吧| 亚洲性夜色夜夜综合| 女性生殖器流出的白浆| 日本黄色日本黄色录像| avwww免费| 国产精品久久久av美女十八| 村上凉子中文字幕在线| 亚洲成人久久性| 国产精品1区2区在线观看.| 欧美乱码精品一区二区三区| 久久国产亚洲av麻豆专区| 亚洲男人的天堂狠狠| aaaaa片日本免费| 久久久久九九精品影院| 香蕉丝袜av| 亚洲熟妇中文字幕五十中出 | 69精品国产乱码久久久| 午夜福利在线观看吧| 天堂中文最新版在线下载| 叶爱在线成人免费视频播放| 欧美午夜高清在线| 美女午夜性视频免费| 欧美日韩av久久| 国产高清videossex| 国产麻豆69| 一个人观看的视频www高清免费观看 | 老司机亚洲免费影院| 日本三级黄在线观看| 国产欧美日韩综合在线一区二区| av超薄肉色丝袜交足视频| 国产黄a三级三级三级人| 亚洲一码二码三码区别大吗| 波多野结衣av一区二区av| 激情视频va一区二区三区| 亚洲欧美一区二区三区久久| 中出人妻视频一区二区| 久久九九热精品免费| 91av网站免费观看| 国产一区二区激情短视频| 国产91精品成人一区二区三区| 在线观看一区二区三区激情| 国产精品亚洲一级av第二区| 精品久久久精品久久久| 日本三级黄在线观看| 成人手机av| 丝袜美腿诱惑在线| 免费高清视频大片| 久热这里只有精品99| 丁香欧美五月| 欧美成人午夜精品| 精品高清国产在线一区| 在线观看日韩欧美| 波多野结衣一区麻豆| 久久午夜亚洲精品久久| 在线观看午夜福利视频| 久久中文字幕一级| 一个人免费在线观看的高清视频| 一进一出好大好爽视频| 丰满的人妻完整版| a级片在线免费高清观看视频| 精品国产美女av久久久久小说| 精品国产美女av久久久久小说| 午夜福利欧美成人| 日本欧美视频一区| 午夜影院日韩av| 亚洲国产精品sss在线观看 | 免费高清视频大片| 亚洲人成77777在线视频| 亚洲av熟女| 久久精品亚洲av国产电影网| 久久午夜亚洲精品久久| 女人精品久久久久毛片| 搡老乐熟女国产| 中文字幕av电影在线播放| 一区二区三区国产精品乱码| 欧美成人午夜精品| 亚洲美女黄片视频| 黄色女人牲交| 国产av一区二区精品久久| 亚洲激情在线av| 国产精品一区二区三区四区久久 | 午夜亚洲福利在线播放| 午夜福利在线观看吧| 国产极品粉嫩免费观看在线| 女人被狂操c到高潮| 国产精品永久免费网站| 久久国产亚洲av麻豆专区| 一本大道久久a久久精品| 亚洲色图av天堂| 亚洲成人免费电影在线观看| 999久久久国产精品视频| 国产av一区二区精品久久| 午夜福利免费观看在线| 一个人观看的视频www高清免费观看 | 美女午夜性视频免费| 男女做爰动态图高潮gif福利片 | 夜夜看夜夜爽夜夜摸 | 欧美+亚洲+日韩+国产| 97碰自拍视频| 午夜福利一区二区在线看| 国产一区二区三区综合在线观看| 免费高清在线观看日韩| 中国美女看黄片| 免费在线观看亚洲国产| 精品久久久久久久久久免费视频 | 女性被躁到高潮视频| 三上悠亚av全集在线观看| 精品一区二区三卡| 99在线人妻在线中文字幕| cao死你这个sao货| 人人妻人人澡人人看| 久久精品国产99精品国产亚洲性色 | 天天躁夜夜躁狠狠躁躁| 午夜a级毛片| 免费高清在线观看日韩| 在线观看免费日韩欧美大片| 色婷婷av一区二区三区视频| 夫妻午夜视频| 精品国内亚洲2022精品成人| 欧美日韩av久久| 波多野结衣高清无吗| 欧美黄色片欧美黄色片| 一级a爱视频在线免费观看| 侵犯人妻中文字幕一二三四区| 国产男靠女视频免费网站| 久久中文看片网| 天堂中文最新版在线下载| 99在线视频只有这里精品首页| 国产单亲对白刺激| 久久 成人 亚洲| 少妇 在线观看| 午夜精品国产一区二区电影| 精品电影一区二区在线| 欧美在线黄色| 国产91精品成人一区二区三区| 亚洲一码二码三码区别大吗| 久久久精品国产亚洲av高清涩受| 50天的宝宝边吃奶边哭怎么回事| 免费在线观看视频国产中文字幕亚洲| 女同久久另类99精品国产91| 好男人电影高清在线观看| 精品国产美女av久久久久小说| 丝袜美腿诱惑在线| 欧美成狂野欧美在线观看| 国产精品久久电影中文字幕| 亚洲国产欧美网| 纯流量卡能插随身wifi吗| 久久精品国产清高在天天线| 亚洲欧美一区二区三区黑人| 很黄的视频免费| 久久热在线av| 水蜜桃什么品种好| 男女午夜视频在线观看| 亚洲av美国av| 日本a在线网址| 欧美黄色淫秽网站| 午夜日韩欧美国产| 国产精品综合久久久久久久免费 | 久久人人爽av亚洲精品天堂| 亚洲欧美日韩高清在线视频| 成人18禁高潮啪啪吃奶动态图| 91精品国产国语对白视频| 欧美乱码精品一区二区三区| 午夜两性在线视频| 搡老岳熟女国产| 国产有黄有色有爽视频| 黄色片一级片一级黄色片| 久久人妻av系列| 亚洲精品成人av观看孕妇| avwww免费| 久久人妻熟女aⅴ| 精品国产一区二区久久| 日韩大尺度精品在线看网址 | 亚洲一区二区三区不卡视频| 黄网站色视频无遮挡免费观看| av欧美777| 日日干狠狠操夜夜爽| 久久精品国产99精品国产亚洲性色 | 国产亚洲av高清不卡| 午夜免费成人在线视频| 亚洲精品一二三| 国内毛片毛片毛片毛片毛片| 久久九九热精品免费| 免费久久久久久久精品成人欧美视频| 国产精品乱码一区二三区的特点 | 成人特级黄色片久久久久久久| 免费看a级黄色片| 天堂√8在线中文| 99riav亚洲国产免费| 黑人猛操日本美女一级片| 无限看片的www在线观看| 日韩免费av在线播放| 亚洲,欧美精品.| av超薄肉色丝袜交足视频| www.熟女人妻精品国产| 国产精品一区二区精品视频观看| 欧美av亚洲av综合av国产av| 国产精品自产拍在线观看55亚洲| 99久久国产精品久久久| a级毛片黄视频| 日韩精品青青久久久久久| 天堂√8在线中文| 中文欧美无线码| 亚洲国产中文字幕在线视频| 99久久99久久久精品蜜桃| 亚洲少妇的诱惑av| 国内毛片毛片毛片毛片毛片| 亚洲成a人片在线一区二区| 精品久久久久久久久久免费视频 | 成年版毛片免费区| 欧美日韩视频精品一区| 国产视频一区二区在线看| 少妇粗大呻吟视频| 久久久精品欧美日韩精品| 国产一区二区三区在线臀色熟女 | av电影中文网址| 日本五十路高清| 露出奶头的视频| 18禁黄网站禁片午夜丰满| 在线国产一区二区在线| 乱人伦中国视频| 91成年电影在线观看| 精品一区二区三区四区五区乱码| 美女福利国产在线| 人人妻,人人澡人人爽秒播| 日韩视频一区二区在线观看| av福利片在线| 久久久久久久久中文| 97碰自拍视频| 黑丝袜美女国产一区| 精品国产一区二区久久| 成人永久免费在线观看视频| 精品一区二区三区av网在线观看| 9热在线视频观看99| 国产精品国产av在线观看| 99riav亚洲国产免费| 黑人猛操日本美女一级片| 女同久久另类99精品国产91| 欧美激情 高清一区二区三区| 黄片大片在线免费观看| 级片在线观看| 亚洲精品在线美女| 99国产精品免费福利视频| 午夜福利在线观看吧| 法律面前人人平等表现在哪些方面| 又紧又爽又黄一区二区| 久久天堂一区二区三区四区| 满18在线观看网站| 色尼玛亚洲综合影院| 亚洲熟女毛片儿| 欧美激情 高清一区二区三区| 大香蕉久久成人网| 91老司机精品| 国产成人欧美| 丰满的人妻完整版| 久9热在线精品视频| 亚洲色图综合在线观看| 国产欧美日韩一区二区三区在线| 淫秽高清视频在线观看| 中文字幕色久视频| 又紧又爽又黄一区二区| 中文字幕人妻丝袜一区二区| 别揉我奶头~嗯~啊~动态视频| 久久热在线av| 老汉色∧v一级毛片| 中文字幕人妻丝袜一区二区| 无人区码免费观看不卡| 一区二区三区国产精品乱码| 村上凉子中文字幕在线| 亚洲精华国产精华精| xxx96com| 天天躁夜夜躁狠狠躁躁| 叶爱在线成人免费视频播放| www国产在线视频色| av免费在线观看网站| 午夜福利一区二区在线看| 欧美乱妇无乱码| 一级毛片女人18水好多| 日韩免费av在线播放| 免费观看人在逋| 欧洲精品卡2卡3卡4卡5卡区| 亚洲欧美激情在线| 免费在线观看完整版高清| 青草久久国产| 黄片播放在线免费| www.自偷自拍.com| 免费高清在线观看日韩| 亚洲欧美激情综合另类| 欧美日韩瑟瑟在线播放| 亚洲成人国产一区在线观看| 欧美乱妇无乱码| 国产亚洲精品综合一区在线观看 | 9热在线视频观看99| 亚洲视频免费观看视频| 在线观看免费午夜福利视频| 91精品国产国语对白视频| 国产不卡一卡二| 18禁裸乳无遮挡免费网站照片 | 在线观看免费视频日本深夜| 99国产精品免费福利视频| 99久久国产精品久久久| 国产亚洲精品一区二区www| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲黑人精品在线| av天堂久久9| 自线自在国产av| 18禁美女被吸乳视频| 久久久久精品国产欧美久久久| 99久久久亚洲精品蜜臀av| 黄色视频,在线免费观看| av网站在线播放免费| 成人特级黄色片久久久久久久| 久久久国产精品麻豆| 国产成人啪精品午夜网站| 久久久久久免费高清国产稀缺| 黄色视频不卡| 成熟少妇高潮喷水视频| 色哟哟哟哟哟哟| 日本免费a在线| 国产高清videossex| 999久久久国产精品视频| 狠狠狠狠99中文字幕| 50天的宝宝边吃奶边哭怎么回事| 男女午夜视频在线观看| 在线看a的网站| 国产精品电影一区二区三区| 国产男靠女视频免费网站| 午夜免费激情av| 欧美人与性动交α欧美软件| 他把我摸到了高潮在线观看| 美女 人体艺术 gogo| www日本在线高清视频| 天堂√8在线中文| 欧美日韩亚洲国产一区二区在线观看| 在线观看免费视频日本深夜| 一个人免费在线观看的高清视频| 黑人猛操日本美女一级片| 18禁观看日本| 久久精品91无色码中文字幕| 看黄色毛片网站| 亚洲成人精品中文字幕电影 | 电影成人av| 久久国产精品人妻蜜桃| 97碰自拍视频| 亚洲精品中文字幕一二三四区| √禁漫天堂资源中文www| 久久国产精品影院| 91精品三级在线观看| 免费在线观看亚洲国产| 不卡av一区二区三区| 黄片小视频在线播放| av片东京热男人的天堂| 亚洲va日本ⅴa欧美va伊人久久| 1024香蕉在线观看| 亚洲人成电影免费在线| 一本综合久久免费| 亚洲片人在线观看| 欧美乱妇无乱码| 老司机深夜福利视频在线观看| 巨乳人妻的诱惑在线观看| 十八禁人妻一区二区| 亚洲欧美一区二区三区久久| 国产精品久久久av美女十八| 国产精华一区二区三区| 黄色毛片三级朝国网站| 亚洲精品久久成人aⅴ小说| 淫妇啪啪啪对白视频| 中文字幕色久视频| 亚洲国产欧美日韩在线播放| 精品国产乱码久久久久久男人| 黄色a级毛片大全视频| avwww免费| 日本vs欧美在线观看视频| 国产亚洲精品久久久久久毛片| 女性生殖器流出的白浆| 免费久久久久久久精品成人欧美视频| 少妇粗大呻吟视频| av超薄肉色丝袜交足视频| 看黄色毛片网站| 午夜福利,免费看| 欧美一级毛片孕妇| 自线自在国产av| 久久天躁狠狠躁夜夜2o2o| 亚洲一区二区三区不卡视频| 亚洲精品一二三| 欧美日韩亚洲国产一区二区在线观看| 99久久久亚洲精品蜜臀av| 国产成人精品无人区| 亚洲在线自拍视频| 欧美av亚洲av综合av国产av| 99热国产这里只有精品6| 欧美不卡视频在线免费观看 | 国产欧美日韩精品亚洲av| 后天国语完整版免费观看| 国产在线精品亚洲第一网站| 国产野战对白在线观看| 久久人妻熟女aⅴ| 成人影院久久| 亚洲伊人色综图| tocl精华| 国产精品一区二区精品视频观看| 亚洲狠狠婷婷综合久久图片| 91老司机精品| 国产精品久久久人人做人人爽| 97超级碰碰碰精品色视频在线观看| 黄色 视频免费看| 人成视频在线观看免费观看| 9191精品国产免费久久| 欧美激情 高清一区二区三区| 老司机福利观看| av中文乱码字幕在线| 日本一区二区免费在线视频| 妹子高潮喷水视频| 亚洲精华国产精华精| 免费在线观看视频国产中文字幕亚洲| 久久久久久久午夜电影 | 一a级毛片在线观看| av片东京热男人的天堂| 精品人妻1区二区| 亚洲欧美激情综合另类| 国产精品亚洲一级av第二区| 亚洲av熟女| 91在线观看av| 日本精品一区二区三区蜜桃| 多毛熟女@视频| 丁香欧美五月| 久久久久久亚洲精品国产蜜桃av| 91麻豆av在线| 亚洲精品在线观看二区| 久久九九热精品免费| 欧美不卡视频在线免费观看 | 美女高潮到喷水免费观看| 黑人欧美特级aaaaaa片| 国产精品久久久人人做人人爽| 大码成人一级视频| 日韩欧美国产一区二区入口| 欧美人与性动交α欧美精品济南到| 亚洲精品成人av观看孕妇| 久久99一区二区三区| 无限看片的www在线观看| 露出奶头的视频| 日日夜夜操网爽| 亚洲成av片中文字幕在线观看| 制服人妻中文乱码| 久久天堂一区二区三区四区| 免费av中文字幕在线| 美女午夜性视频免费| 欧美丝袜亚洲另类 | 国产亚洲精品第一综合不卡| 国产精品美女特级片免费视频播放器 | 男女下面进入的视频免费午夜 | 欧美不卡视频在线免费观看 | 91国产中文字幕| 麻豆成人av在线观看| 欧美在线一区亚洲| a级毛片黄视频| 亚洲一码二码三码区别大吗| 岛国在线观看网站| 大型av网站在线播放| 色婷婷久久久亚洲欧美| 精品免费久久久久久久清纯| 午夜91福利影院| 国产成人精品在线电影| 精品高清国产在线一区| 色综合站精品国产| 看黄色毛片网站| 黄色视频,在线免费观看| 精品久久久精品久久久| 国产精品av久久久久免费| 免费在线观看日本一区| 天堂√8在线中文| 一级毛片精品| 美女高潮到喷水免费观看| 精品高清国产在线一区| 淫秽高清视频在线观看| 久久精品国产综合久久久| 午夜福利在线免费观看网站| 在线观看免费高清a一片| 啦啦啦 在线观看视频| 久久青草综合色| 欧美日本中文国产一区发布| xxxhd国产人妻xxx| 国产精品日韩av在线免费观看 | 久久国产亚洲av麻豆专区| 国产成人精品在线电影| netflix在线观看网站| 久久中文字幕一级| 欧美乱色亚洲激情| а√天堂www在线а√下载| 久99久视频精品免费| 超碰97精品在线观看| 国产熟女午夜一区二区三区| 一区二区三区激情视频| 真人一进一出gif抽搐免费| 每晚都被弄得嗷嗷叫到高潮| 亚洲人成电影免费在线| 伊人久久大香线蕉亚洲五| 久久这里只有精品19| 无遮挡黄片免费观看| 首页视频小说图片口味搜索| 色在线成人网| 精品午夜福利视频在线观看一区| 日本欧美视频一区| svipshipincom国产片| 亚洲成人国产一区在线观看| 国产aⅴ精品一区二区三区波| 日韩欧美三级三区| 黄色a级毛片大全视频| 久久精品影院6| 91精品国产国语对白视频| 黑丝袜美女国产一区| 国产熟女xx| 人成视频在线观看免费观看| 一级黄色大片毛片| 如日韩欧美国产精品一区二区三区| 动漫黄色视频在线观看| 成人黄色视频免费在线看| 国产深夜福利视频在线观看| 国产单亲对白刺激| 一区二区日韩欧美中文字幕| 成人影院久久| 19禁男女啪啪无遮挡网站| 咕卡用的链子| 在线看a的网站| 国产精品亚洲一级av第二区| 国产熟女午夜一区二区三区|