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

    擺動(dòng)河槽水動(dòng)力穩(wěn)定性特征分析1)

    2017-03-21 10:51:51白玉川冀自青徐海玨
    力學(xué)學(xué)報(bào) 2017年2期
    關(guān)鍵詞:波數(shù)雷諾數(shù)水流

    白玉川 冀自青 徐海玨,2)

    ?(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室,天津300072)

    ?(天津大學(xué)建筑工程學(xué)院,天津300072)

    擺動(dòng)河槽水動(dòng)力穩(wěn)定性特征分析1)

    白玉川?,?冀自青?徐海玨?,?,2)

    ?(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室,天津300072)

    ?(天津大學(xué)建筑工程學(xué)院,天津300072)

    河流形態(tài)與水動(dòng)力結(jié)構(gòu)息息相關(guān),形態(tài)約束水動(dòng)力結(jié)構(gòu),水動(dòng)力結(jié)構(gòu)則通過泥沙運(yùn)動(dòng)進(jìn)一步塑造形態(tài),在自然界河流中形成一對(duì)辯證互饋關(guān)系.天然河流形態(tài)形式多樣,大致可以分為順直、微彎、分叉和散亂游蕩幾種類型,其中微彎及多個(gè)彎曲構(gòu)成的河型為河流動(dòng)力演化中最重要的一環(huán).多個(gè)彎曲構(gòu)成的河型可用正弦派生曲線來描述,它也是天然河流主槽與水動(dòng)力結(jié)構(gòu)復(fù)雜相互作用的結(jié)果.作為探討這一過程的力學(xué)作用機(jī)理,構(gòu)建擺動(dòng)槽道并研究槽道擺動(dòng)與其內(nèi)部流動(dòng)結(jié)構(gòu)的互饋關(guān)系,既是流體力學(xué)研究的熱點(diǎn)內(nèi)容,也是目前河流動(dòng)力過程研究的基礎(chǔ)內(nèi)容.在此重點(diǎn)討論這一互饋關(guān)系前一部分,即水流對(duì)擺動(dòng)邊界的響應(yīng).文中建立了隨體坐標(biāo)系下擺動(dòng)河槽與內(nèi)部水流動(dòng)力響應(yīng)理論模型,通過給定擺動(dòng)彎曲槽道的不同特征參數(shù),研究討論了正弦派生型擺動(dòng)邊界下的槽道水流動(dòng)力穩(wěn)定性特征,明確了彎曲槽道擺動(dòng)對(duì)其內(nèi)部主流及擾動(dòng)水流結(jié)構(gòu)的影響,確定彎曲槽道擺動(dòng)波數(shù)、擺動(dòng)頻率對(duì)擾動(dòng)流發(fā)展影響的相應(yīng)參數(shù)定量關(guān)系,得到了槽道彎曲度和擺動(dòng)特征對(duì)其內(nèi)部水流不同尺度擾動(dòng)影響的閾值選擇性范圍.

    河槽擺動(dòng),槽道流動(dòng),正弦派生,水流特性,水動(dòng)力穩(wěn)定性特征

    引言

    復(fù)雜邊界與內(nèi)部水流的相互作用既是研究河流擺動(dòng)機(jī)理的基礎(chǔ),也是近來流體力學(xué)研究的熱門課題.

    平行槽道內(nèi)的流動(dòng)在小雷諾數(shù)時(shí)為層流,其流速分布為二次分布;然而,在復(fù)雜可動(dòng)邊界內(nèi)流速將產(chǎn)生變化,此時(shí)其水動(dòng)力穩(wěn)定性特征也會(huì)相應(yīng)變化.

    目前對(duì)于復(fù)雜可動(dòng)邊界內(nèi)流體流動(dòng)特征的研究主要沿著兩個(gè)方向進(jìn)行:一個(gè)是彈性固壁在水流作用下的自由振動(dòng).這方面主要代表學(xué)者有Davies和Carpenter[1-2],Hell和Waters[3],Guaus等[4],Pitman和Lucey[5],其核心思想是將邊壁構(gòu)造成無數(shù)彈性支承的薄壁梁,對(duì)彈性邊界和內(nèi)部水流分別建立控制方程,研究彈性邊壁對(duì)平面Poiseuille流、Taylor–Couette流水動(dòng)力穩(wěn)定性的影響.另一個(gè)是以水流減阻為主要目標(biāo),研究靜態(tài)或動(dòng)態(tài)復(fù)雜邊界內(nèi)層流運(yùn)動(dòng)的水動(dòng)力穩(wěn)定性和湍流擬序結(jié)構(gòu)的發(fā)展.例如,Hall[6]、Thomas等[7]探討了邊壁縱向振動(dòng)對(duì)平面Poiseuille流水動(dòng)力穩(wěn)定性特征的影響;Kuhlmann等[8]針對(duì)兩壁逆向旋轉(zhuǎn)的空腔,研究了其內(nèi)部水流的穩(wěn)定性;Hell和Waters[3]則針對(duì)彈性圓柱殼中的水流結(jié)構(gòu)進(jìn)行研究;Ren等[9]分析并得到了邊壁縱向勻速平移對(duì)內(nèi)部水流水動(dòng)力穩(wěn)定性的影響;Sen等[10]研究了固定邊壁和可動(dòng)邊壁交替出現(xiàn)時(shí)T-S波在邊界上的傳播.

    由于研究對(duì)象的不同,以上這些成果大多研究邊壁縱向運(yùn)動(dòng)對(duì)內(nèi)部水流結(jié)構(gòu)和水動(dòng)力穩(wěn)定性的影響,很少有模型考慮邊壁橫向擺動(dòng)對(duì)內(nèi)部水流的影響.然而對(duì)于彎曲河流的演化機(jī)理[11],槽道橫向擺動(dòng)邊界內(nèi)水流運(yùn)動(dòng)規(guī)律的研究更具有實(shí)際意義[12-13].

    因此,本文首先假設(shè)河槽的擺動(dòng)趨勢(shì),建立了隨體坐標(biāo)系下擺動(dòng)性槽道內(nèi)水流流動(dòng)的理論模型,討論了正弦派生型擺動(dòng)邊界下的槽道水流動(dòng)力穩(wěn)定性特征,其次確定了擺動(dòng)波數(shù)與擺動(dòng)頻率對(duì)主流流速和內(nèi)部擾動(dòng)發(fā)展影響的定量關(guān)系,最后分析并得到了槽道彎曲度和擺動(dòng)特征對(duì)其內(nèi)部水流不同尺度擾動(dòng)的選擇性影響.

    1 基本方程

    天然河道,尤其是下游河段一般是寬淺型的(河寬與水深比值不小于1).因此,河道在下游彎曲河段會(huì)呈現(xiàn)出典型的三維流動(dòng)特征.然而,河道的中上游河段,若經(jīng)過大峽谷,則一般會(huì)形成窄深型的河道(河寬與水深比值小于1).此時(shí),河道水深方向上變量的變化就可忽略不計(jì),表現(xiàn)出典型的平面二維水流流動(dòng)特征,這種河槽中的二維流動(dòng)在Parker等[14],Ikeda等[15]、Bai和Xu[16]以及Xu和Bai[17]的文中都有專門討論和分析.因此,本文也針對(duì)窄深型河槽建立平面二維坐標(biāo)系統(tǒng).

    采用如圖1中所示的坐標(biāo)系,河道中心線是沿著流向且處于河道中心位置的曲線,s?方向的其他曲線都是平行于河道中心線而n?方向是處處垂直于河道中心線并在水平面內(nèi)的方向,因此s?和n?線是處處正交的關(guān)系.這種坐標(biāo)系在Parker等[14]以及Ikeda等[15]的文章中首次使用,并在Bai等[18]的文章中多次討論.

    圖1 正弦派生曲線邊界平面示意圖Fig.1 Sketch of sine-generated boundary

    長(zhǎng)度尺度、速度尺度、時(shí)間尺度分別用半寬B?,零階基本流流速峰值無量綱化;河道曲率用河道最大曲率無量綱化,對(duì)于正弦派生曲線

    正弦派生曲線相關(guān)物理量的無量綱參數(shù):偏角幅值θm,擺動(dòng)波數(shù)αc,擺動(dòng)角頻率ωc,其與實(shí)際物理量直接的關(guān)系式為

    式中,uc為槽道擺動(dòng)的無量綱相速度,亦可寫為uc=ωc/αc.

    物理量無量綱化關(guān)系式為

    擺動(dòng)相位:φc=αcs-ωct;偏角寫為復(fù)數(shù)形式:θ=θmexp(iφc);無量綱曲率函數(shù)的復(fù)數(shù)形式:c=θ,s=iαcθ;無量綱曲率幅值

    無量綱形式的控制方程為

    無滑移邊界條件

    式中,S(s,t)為槽道中心線在拉格朗日坐標(biāo).槽道中心線伸縮比:當(dāng)槽道中心線無伸縮時(shí),?S/?s=1,亦即us0,s=0或us0=us0(t).

    值得注意的是,河道形態(tài)參數(shù)很多具有復(fù)數(shù)形式,其中包括:河道波數(shù)αc、河道頻率ωc、擺動(dòng)相位 φc、無量綱曲率函數(shù)c、Lame系數(shù)hs等.基本所有參數(shù)的實(shí)部具有物理含義,而虛部沒有具體物理含義,因此參數(shù)的取值一般為參數(shù)的實(shí)部.然而,在計(jì)算中需要用到參數(shù)的實(shí)部和虛部.如:偏角θ可以寫作θ=θmexp(iφc)=θmexp[i(φcr+iφci)]=θmexp(-φci)exp(iφcr).其中,φcr和φci分別是擺動(dòng)相位φc的實(shí)部和虛部.因此,擺動(dòng)相位實(shí)部φcr體現(xiàn)了河道周期性變化對(duì)偏角的影響,而其虛部φci則體現(xiàn)了河道非周期擺動(dòng)幅度改變對(duì)偏角的影響.

    2 基本流攝動(dòng)解

    攝動(dòng)法又稱小參數(shù)展開法,就是將系統(tǒng)視為理想模型的參數(shù)或結(jié)構(gòu)作了微小擾動(dòng)的結(jié)果來研究其運(yùn)動(dòng)過程的數(shù)學(xué)方法.這種方法最早應(yīng)用于天體力學(xué),用來計(jì)算小天體對(duì)大天體運(yùn)動(dòng)的影響,后來廣泛應(yīng)用于物理學(xué)和力學(xué)的理論研究.利用攝動(dòng)法求解方程,也稱攝動(dòng)分解,通常需要在無量綱方程中選擇一個(gè)能反映物理特征的無量綱小參數(shù)作為攝動(dòng)量,然后假設(shè)解可以按小參數(shù)展成冪級(jí)數(shù),將這一形式級(jí)數(shù)代入無量綱方程后,可得各級(jí)近似方程,依據(jù)這些方程可確定冪級(jí)數(shù)的系數(shù),對(duì)級(jí)數(shù)進(jìn)行截?cái)?,便得到原方程的漸近解.

    無量綱曲率幅值ψ=B?/r?,表示半河寬與最小曲率半徑之比,它代表了槽道的空間彎曲程度,是重要的彎曲特征參數(shù),這個(gè)幅值越大就證明彎道越尖銳.根據(jù)實(shí)際測(cè)量資料的結(jié)果顯示,天然形成河流槽道的ψ值大多在0.05~0.1的范圍內(nèi)[19-21],幾乎沒有超過0.3的.因此,在對(duì)方程的處理中,一般認(rèn)為ψ是可用于方程攝動(dòng)展開的小參數(shù)量.

    擺動(dòng)角頻率ωc反映擺動(dòng)槽道的時(shí)間特征,中心線橫向速度un0的量級(jí)為ucθm,槽道緩慢擺動(dòng)時(shí)可認(rèn)為其與ψ同量級(jí).

    2.1 基本流

    以無量綱曲率幅值ψ為小參數(shù),對(duì)Xψ進(jìn)行攝動(dòng)分解,只保留一階攝動(dòng)一次諧波量,可得

    式中,第 1項(xiàng)Xψ0為順直項(xiàng)部分 (ψ0),對(duì)應(yīng) ψ=0時(shí)的直槽流,為實(shí)向量函數(shù);第2項(xiàng)為一階彎曲修正項(xiàng),與槽道中心線擺動(dòng)同頻同波數(shù),表示為復(fù)數(shù)形式,其幅值:頂標(biāo)?表示波動(dòng)量的幅值,為復(fù)向量函數(shù)

    ψ0階對(duì)應(yīng)固定直槽道流動(dòng),其對(duì)應(yīng)基本流解為Poiseuille流流速分布

    ψ1階反映了槽道彎曲對(duì)水流流動(dòng)造成的影響,其解對(duì)應(yīng)的方程為

    對(duì)ψ1階方程(7)簡(jiǎn)化,可得一個(gè)O-S方程形式的常微分方程

    結(jié)合邊界條件,差分法求解方程,采用橫向網(wǎng)格數(shù)N=1000,由式(8)可解得橫向流速幅值將其代回方程(7)的連續(xù)方程和s向動(dòng)量方程,可解得一階縱向流速和壓力的幅值

    2.2 擺動(dòng)波數(shù)αc對(duì)一階基本流解的影響

    對(duì)于固定的彎曲槽道,擺動(dòng)角頻率ωc=0.雷諾數(shù)Re=2000時(shí),求解控制方程(7),可得一階流速分量和沿n方向分布圖,見圖2和圖3.

    圖2為一階流速分布在復(fù)數(shù)空間的分解,函數(shù)Re(),Im()分布表示復(fù)數(shù)的實(shí)部和虛部.實(shí)部為相位φ=0的正彎頂處的流速分布,虛部則為相位φ=-π/2的零曲率的直槽處流速分布,提前個(gè)相位.

    圖2顯示,一階s方向流速呈反對(duì)稱分布,而n方向流速則呈對(duì)稱分布;一階s方向流速在近兩壁處有最大值和最小值,極值大小和所處位置隨擺動(dòng)波數(shù)αc增變化而變化.

    隨擺動(dòng)波數(shù)αc增大,一階n方向流速隨擺動(dòng)波數(shù)αc增大明顯增大;一階s方向流速極值所處位置隨擺動(dòng)波數(shù)αc增大而逐漸靠近兩壁,在αc=(0,0.5)的范圍內(nèi)實(shí)部逐漸增大,而虛部則先增大后減小.

    圖2 一階流速沿n方向分布(復(fù)數(shù)空間)Fig.2 Distribution of first-orde velocities alongndirection(complex space)

    這是因?yàn)?,無量綱曲率幅值ψ是擺動(dòng)波數(shù)αc和偏角幅值θm的乘積,當(dāng)偏角幅值保持不變時(shí),ψ會(huì)隨著αc的增大而增大.因此,αc增大就意味著河槽的擺動(dòng)幅度增加,此時(shí),s方向的最大水流流速會(huì)由于慣性的作用偏向凹岸一側(cè),而且擺動(dòng)幅度越大,最大流速的偏離就越嚴(yán)重,甚至開始慢慢貼近凹岸邊界.

    圖3則為一階流速分布在幅值--相位空間的分解,Am()和Ph()分別表示復(fù)數(shù)的幅值和相位(弧度制).由圖3和圖2可知,縱向流速和橫向流速的幅值呈對(duì)稱分布,的相位呈對(duì)稱分布,但的相位在槽道兩側(cè)有180°的相位差,亦即

    圖3 一階流速沿n方向分布(幅值--相位空間)Fig.3 Distribution of first-orde velocities alongndirection(amplitude-phase space)

    2.3 擺動(dòng)角頻率ωc對(duì)一階基本流解的影響

    對(duì)于擺動(dòng)的彎曲槽道,擺動(dòng)角頻率ωc≠0.擺動(dòng)波數(shù)αc=0.3,雷諾數(shù)Re=2000時(shí),求解控制方程(1),所得一階流速分量的橫向分布繪圖,見圖4和圖5.

    圖4 一階流速沿n方向分布(復(fù)數(shù)空間)Fig.4 Distributions of first-orde velocities alongndirection(complex space)

    擺動(dòng)槽道的一階流速分布取決于兩方面的影響,一方面源自于槽道彎曲,即擺動(dòng)波數(shù)αc的影響,一方面源自于擺動(dòng)頻率ωc的影響.圖4中黑色實(shí)線為ωc=0的固定彎曲槽道中的一階流速分布,取ωc=-0.10,-0.05,0,0.05,0.10五種不同的擺動(dòng)情況,并以實(shí)線和虛線區(qū)分?jǐn)[動(dòng)傳播方向,實(shí)線對(duì)應(yīng)ωc>0,虛線對(duì)應(yīng)ωc<0.由圖4可見,擺動(dòng)向下游傳播時(shí),擺動(dòng)頻率ωc對(duì)一階流速影響比向上游傳播大,而且,整體來說,擺動(dòng)頻率ωc對(duì)一階流速的影響與擺動(dòng)波數(shù)對(duì)它的影響是相反的(圖2).這是因?yàn)楹硬壑械乃髁魉俜较蚴窍蛳掠蔚模虼?,水流的信息主要是向下游傳播的,擺動(dòng)向下游傳播間接推動(dòng)了水流信息的傳播效果.然而,又由于水流和河槽邊界同時(shí)向下游運(yùn)動(dòng),它們之間的相對(duì)速度就減小了,就相當(dāng)于間接減小了波數(shù)的影響.因此,擺動(dòng)頻率對(duì)一階流速的影響與擺動(dòng)波數(shù)對(duì)它的影響效果是相反的.

    隨著擺動(dòng)頻率ωc由零減小,一階縱向流速和橫向流速分布迅速反向,且ωc=0.05和ωc=0.10兩種擺動(dòng)頻率下相差不大.

    隨著擺動(dòng)頻率ωc由零增大,ωc=0.05時(shí),彎頂處兩壁附近出現(xiàn)一階縱向流速的反向分布 (與 ωc=0時(shí)的黑線分布對(duì)比),這時(shí)擺動(dòng)波數(shù)αc和擺動(dòng)頻率ωc的影響相當(dāng);擺動(dòng)頻率繼續(xù)增大ωc=0.10時(shí),反向流速分布擴(kuò)展到全斷面,這時(shí)擺動(dòng)頻率ωc的影響占主導(dǎo)地位.直槽處的一階橫向流速分布隨擺動(dòng)波數(shù)的變化與彎頂處分布相似.二者隨ωc的復(fù)雜變化規(guī)律是由于邊界流速非零,因此其近壁流速才呈現(xiàn)局部反向分布的現(xiàn)象.

    圖5為一階流速分布在幅值--相位空間的分解,一階縱向流速和橫向流速的幅值和相位均呈對(duì)稱分布.

    同時(shí),從圖5(b)中可以觀察到ωc=0.10對(duì)應(yīng)的藍(lán)色實(shí)線從1.0突然跳轉(zhuǎn)到-1.0,這是因?yàn)閳D5(b)中第2個(gè)圖代表的是unψ1的相位,而通過反正切計(jì)算相位的值域在(-π,π)之間.這種跳轉(zhuǎn)體現(xiàn)了相位的值域范圍的約束,其本身只有數(shù)學(xué)上的意義,在物理上Ph(unψ1)/π仍然是一個(gè)連續(xù)函數(shù).

    圖5 一階流速橫向分布(幅值--相位空間)Fig.5 Distributions of first-orde velocities(amplitude-phase space)

    3 水動(dòng)力穩(wěn)定性及擬序擾動(dòng)發(fā)展的邊界效應(yīng)

    按照描述湍流擬序結(jié)構(gòu)理論方法[23-25],方程的解X=[us,un,p]T分解為基本流解Xψ與擾動(dòng)解XT之和的形式,X=Xψ+XT.

    與平面 Poiseuille流相比,擺動(dòng)槽道基本流中的波動(dòng)項(xiàng)部分Xψ1的使內(nèi)部水流既有流動(dòng)的自然失穩(wěn)又有由于擺動(dòng)彎曲誘發(fā)的失穩(wěn),擾動(dòng)量XT=可表述為以下形式

    式中,擾動(dòng)相位φT=αTs-ωTt,由于基本量Xψ對(duì)ψ取一階近似,擾動(dòng)解中包含-1,0,1次諧波量,故式(9)在m=±1處截?cái)?,XT的幅值為n的復(fù)數(shù)函數(shù),

    將式(9)中的擬序擾動(dòng)XT、基本解Xψ代入原始控制方程(2)~方程(4),得到擬序擾動(dòng)的控制方程為

    邊界條件n=±1,usT=unT=0.

    式(10)~式(12)為擬序結(jié)構(gòu)擾動(dòng)項(xiàng)控制方程.方程中既有邊界波動(dòng)(拉梅系數(shù)hs)對(duì)擾動(dòng)解的直接影響,也有邊界波動(dòng)引起的一階基本流Xψ1對(duì)擾動(dòng)解的間接影響.擾動(dòng)方程的非線性使得各擾動(dòng)諧波量相互耦合,這種耦合有兩方面,一種是由高階基本流(usψi,unψi)的色散作用引起,一種是由拉梅系數(shù)hs的Taylor高階展開項(xiàng)的色散作用引起.由于本文只保留一階基本流和三種擾動(dòng)諧波量,前者只對(duì)相鄰擾動(dòng)諧波起作用,而拉梅系數(shù)hs則在各擾動(dòng)分量之間都能起到影響,故需保留hs的ψ的二階展開項(xiàng),即

    式中,c.c表示復(fù)共軛.將式(13)代入擬序擾動(dòng)方程(10)~方程(12),可得擾動(dòng)幅值方程

    式中,m=0,±1,式(14)~式(16)構(gòu)成9元耦合方程組.αT和ωT分別為無量綱形式的擾動(dòng)波數(shù)和擾動(dòng)角頻率;αc和ωc分別為無量綱邊界擺動(dòng)波數(shù)和擺動(dòng)角頻率.

    式 (14)~式 (16)代表彎曲河道水動(dòng)力穩(wěn)定性特征.方程組中未知量包含9個(gè)線性無關(guān)的特征變量:這些特征變量均為自變量n的特征函數(shù),且相互耦合,以矩陣形式可表述為

    系數(shù)矩陣A,B,C各元素的顯示表述可由方程(14)~方程(16)及附錄得到,限于篇幅不再贅述.系數(shù)矩陣中共有兩組參數(shù),一組為與直槽一樣的擬序擾動(dòng)參數(shù)組(αT,ωT,Re),還有一組為邊界彎曲參數(shù)組(θm,αc,ωc).

    特征值方程(17)采用時(shí)間模式,αT為實(shí)數(shù),ωT為復(fù)數(shù),用muller法求解擾動(dòng)幅值方程的特征值.擾動(dòng)角頻率ωT為(αT,Re)和(θm,αc,ωc)的函數(shù),亦即:ωT=ωT(αT,Re,θm,αc,ωc),其虛部即為擾動(dòng)增長(zhǎng)率:ωTi=ωTi(αT,Re,θm,αc,ωc).

    正弦派生曲線河槽內(nèi)的擾動(dòng)波,不僅與擾動(dòng)波數(shù)αT和雷諾數(shù)R的有關(guān),還受到邊界擺動(dòng)彎曲(θm, αc,ωc)的影響.直槽水流雷諾數(shù)Recr=5772.222 (Orszag[26]),而彎曲槽道流的臨界雷諾數(shù)Recr受到邊界擺動(dòng)彎曲參數(shù)(θm,αc,ωc)的影響,亦即:Recr=Recr(θm,αc,ωc),與此同時(shí),臨界雷諾數(shù)對(duì)應(yīng)的臨界擾動(dòng)參數(shù)(αTcr,ωTcr)均為(θm,αc,ωc)的函數(shù),亦即,αTcr=αTcr(θm,αc,ωc),ωTcr=ωTcr(θm,αc,ωc).

    由于 ψ=0時(shí)矩陣奇異,本文以θm=αc=1×10-6代替順直河道(θm=αc=0)來驗(yàn)證臨界雷諾數(shù)及中性曲線,采用目前國(guó)際最為公認(rèn)的理論和實(shí)驗(yàn)結(jié)果驗(yàn)證模型的準(zhǔn)確性,并與Reynold和Potter[27]的數(shù)值結(jié)果和Nishioka和Ichikawa[28]的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比,見圖6所示.可以看出理論計(jì)算結(jié)果與已發(fā)表的最為公認(rèn)的中性曲線結(jié)果吻合較好,說明程序計(jì)算順直河道計(jì)算結(jié)果的正確性.

    圖6 順直河道中性曲線計(jì)算結(jié)果驗(yàn)證Fig.6 Comparison of numerical results with data of straight channel

    3.1 邊界擺動(dòng)波數(shù)的影響

    圖7為擾動(dòng)頻率ωTr及擾動(dòng)增長(zhǎng)率ωTi隨擺動(dòng)波數(shù)αc的變化關(guān)系,對(duì)應(yīng)平面狀態(tài)參數(shù):θm=0.1, ωc=0,擬序擾動(dòng)參數(shù):αT=1,Re=6000;圖8臨界雷諾數(shù)隨擺動(dòng)波數(shù)αc的變化,對(duì)應(yīng)平面狀態(tài)參數(shù):θm=0.1,ωc=0.

    從圖7(a)中可以看出,隨著擺動(dòng)波數(shù)的增加,擾動(dòng)頻率有持續(xù)減少的趨勢(shì).造成這種變化趨勢(shì)的原因可能是擾動(dòng)頻率與擾動(dòng)波數(shù)、擺動(dòng)頻率和擺動(dòng)波數(shù)的一種適應(yīng)過程.即,擺動(dòng)波數(shù)在逐漸接近于擾動(dòng)波數(shù)過程中,擾動(dòng)頻率也逐漸接近于擺動(dòng)頻率(ωc=0).這種結(jié)果雖然與Bai等[18]以及Xu等[17]的結(jié)果由于采用不同坐標(biāo)系而產(chǎn)生一定的差別,但總體趨勢(shì)是一致的,體現(xiàn)出模型計(jì)算的合理性.

    圖7 擾動(dòng)頻率ωTr及擾動(dòng)增長(zhǎng)率ωTi隨擺動(dòng)波數(shù)αc的變化(平面狀態(tài)參數(shù):θm=0.1,ωc=0;擬序擾動(dòng)參數(shù):αT=1.02,Re=5772.222)Fig.7 Variation of disturbance frequency ωTrand growth rate ωTiwith swinging wave number αc(for:θm=0.1,ωc=0;αT=1.02,Re=5772.222)

    圖8 臨界雷諾數(shù)隨擺動(dòng)波數(shù)αc的變化(平面狀態(tài)參數(shù):θm=0.1,ωc=0)Fig.8 Variation of critical Reynolds number with swinging wave number αc(for:θm=0.1,ωc=0)

    從圖7(b)中可以看出,隨著擺動(dòng)波數(shù)的增加,擾動(dòng)增長(zhǎng)率出現(xiàn)先增大后減小的趨勢(shì).槽道彎曲對(duì)其內(nèi)部水流各種尺度的擾動(dòng)都有選擇性影響.擺動(dòng)波數(shù)αcm在(0,0.48)區(qū)間上,槽道擺動(dòng)加速擾動(dòng)增長(zhǎng),水流比直槽道水流更不穩(wěn)定.當(dāng)擺動(dòng)波數(shù)αcm=0.3,擾動(dòng)增長(zhǎng)最快,水流最不穩(wěn)定,此時(shí)臨界雷諾數(shù)不足直槽流的一半.

    直槽水流擾動(dòng)臨界參數(shù)Recr=5772.222,αTcr=1.020,ωTcr=0.27,擺動(dòng)槽道中心曲線和臨界擾動(dòng)顯然受擺動(dòng)參數(shù)(θm,αc,ωc)的影響.

    圖9 中性曲線隨擺動(dòng)波數(shù)αc的變化(平面狀態(tài)參數(shù):θm=0.1, ωc=0,0≤αc≤0.5)Fig.9 Variation of neutral curve with swinging wave number αc(for:θm=0.1,ωc=0,0≤αc≤0.5)

    圖9繪制了隨擺動(dòng)波數(shù)αc增大水流中心曲線的移動(dòng),圖中給出了αc=0,0.1,0.2,0.3,0.4,0.5時(shí)的中性曲線.由圖9可以看出,在擺動(dòng)波數(shù)αc較小時(shí),中性曲線變化緩慢,擺動(dòng)波數(shù)αc較大時(shí),中性曲線變化迅速,且趨于平坦化.隨擺動(dòng)波數(shù)αc增大,中性曲線向左移動(dòng),對(duì)應(yīng)臨界擾動(dòng)也隨之左移,臨界雷諾數(shù)減小,至αc=0.3時(shí)達(dá)到最小(亦可見圖8),然后中心曲線右移,臨界雷諾數(shù)增大.

    圖10只顯示了5條中心曲線的詳細(xì)信息,計(jì)算αc在(0,0.5)區(qū)間上共50條中性曲線,將各中性曲線臨界擾動(dòng)對(duì)應(yīng)的點(diǎn)列(Recr,αTcr)和(Recr,ωTcr)分別連接形成臨界擾動(dòng)曲線,繪制在圖8.由圖8可見,隨著擺動(dòng)波數(shù)αc增大,臨界擾動(dòng)對(duì)應(yīng)的臨界擾動(dòng)波數(shù)αTcr先增大后減小然后再次增大,并在臨近αc=0.3附近形成一個(gè)“繩套”;而臨界擾動(dòng)頻率ωTcr.則先增大后減小,形成類似于直槽中性曲線的“拇指”形狀.

    圖10 臨界擾動(dòng)點(diǎn)隨擺動(dòng)波數(shù)αc的變化(平面狀態(tài)參數(shù):θm=0.1,ωc=0)Fig.10 Variation of critical point with swinging wave number αc(for:θm=0.1,ωc=0)

    3.2 邊界擺動(dòng)頻率的影響

    圖11為擾動(dòng)頻率ωTr及擾動(dòng)增長(zhǎng)率ωTi隨擺動(dòng)頻率ωc的變化關(guān)系,對(duì)應(yīng)平面狀態(tài)參數(shù):θm=0.1, αc=0.3,擬序擾動(dòng)參數(shù):αT=1.016,Re=2307,該參數(shù)組合為圖8中最小臨界雷諾數(shù)對(duì)應(yīng)的擾動(dòng)情況.由圖可見,槽道擺動(dòng)頻率對(duì)其內(nèi)部水流各種尺度的擾動(dòng)有選擇性影響;擬序擾動(dòng)對(duì)擺動(dòng)頻率ωc非常敏感,這是因?yàn)閿_動(dòng)頻率的微小變化都能極大改變一階流速分布,如圖4和圖5所示,流速分布的改變進(jìn)一步影響擬序擾動(dòng)的發(fā)展.

    由圖11擾動(dòng)特征值隨擺動(dòng)頻率ωc的變化曲線可以看出,不同模態(tài)之間的交叉使得曲線偶有不可導(dǎo)的尖點(diǎn)和急劇變化的區(qū)段.這是因?yàn)?,任意一種流動(dòng)中都包含有各種尺度不同的擾動(dòng)波,這種擾動(dòng)波體現(xiàn)在特征方程的計(jì)算中則是可以計(jì)算出一系列的特征值.而湍流流動(dòng)通常是由其中最不穩(wěn)定/最少穩(wěn)定的特征值模態(tài)所控制,當(dāng)流動(dòng)條件發(fā)生改變時(shí)有時(shí)會(huì)出現(xiàn)第二甚至第三模態(tài)增長(zhǎng)率急劇增加,超過了第一模態(tài)的增長(zhǎng)率,從而發(fā)生模態(tài)交叉現(xiàn)象.這種現(xiàn)象在波紋狀底邊界[16]水流結(jié)構(gòu)的討論中首次被發(fā)現(xiàn),而在波狀側(cè)邊界[18]的流動(dòng)中被詳細(xì)討論,是目前彎曲型河槽水流結(jié)構(gòu)與順直型河槽水流結(jié)構(gòu)的重要區(qū)分之一.

    圖11 擾動(dòng)頻率ωTr和擾動(dòng)增長(zhǎng)率ωTi隨擺動(dòng)頻率ωc的變化(平面狀態(tài)參數(shù):θm=0.1,-0.1≤ωc≤0.1;擬序擾動(dòng)參數(shù):αT=1.016,Re=2307)Fig.11 Variation of disturbance frequency ωTrand growth rate ωTiwith swinging frequency ωc(for:θm=0.1,-0.1≤ωc≤0.1;αT=1.016,Re=2307)

    擾動(dòng)增長(zhǎng)率ωTi在ωc=0.085和ωc=-0.039處有極大值,在ωc=-0.006附近有極小值,該極小值處模態(tài)交叉;ωc=-0.009和ωc=0.064附近,擾動(dòng)頻率ωTr分別有極大值和極小值.

    層流轉(zhuǎn)捩為湍流,能量消耗會(huì)由于分子運(yùn)動(dòng)的加劇而急劇增加.然而根據(jù)最小耗能原理,自然的壁面邊界和水流結(jié)構(gòu)的相互適應(yīng)過程中會(huì)自然遵循使水流的能耗降低,即從自然現(xiàn)象上則表現(xiàn)為推遲層流向湍流的轉(zhuǎn)捩過程.Crosato[29]在論文中與VanBalen討論了蜿蜒型河道的擺動(dòng)規(guī)律,發(fā)現(xiàn)河道的發(fā)展不僅會(huì)向下游移動(dòng),在一些特定的條件下會(huì)向上游移動(dòng),其方向主要取決于沖刷池相對(duì)于河道頂點(diǎn)的位置[30].在相當(dāng)陡峻的(窄深型)河岸,最大近岸流速將出現(xiàn)在河道頂點(diǎn)的上游,上游的河岸在水流流速的作用下變?nèi)?,從而造成河道向上游方向遷移.在計(jì)算中就體現(xiàn)兩個(gè)重要的趨勢(shì):(1)模態(tài)出現(xiàn)交叉,控制湍流結(jié)構(gòu)的不穩(wěn)定模態(tài)發(fā)生交換,從而最大的近岸流速位置出現(xiàn)在河道頂點(diǎn)的上游而非通常的下游;(2)窄深型河道在水流結(jié)構(gòu)的影響下,開始向上游遷移.表現(xiàn)在方程特征值上則是向上游遷移的擺動(dòng)頻率在一定范圍內(nèi)減小了擾動(dòng)增長(zhǎng)率(如圖11(b)),即河道向上游遷移時(shí),能量耗散更小,水流結(jié)構(gòu)更加穩(wěn)定,如圖12中所顯示的趨勢(shì).

    圖12為臨界雷諾數(shù)Recr隨擺動(dòng)頻率ωc的變化關(guān)系,對(duì)應(yīng)平面狀態(tài)參數(shù):θm=0.1,αc=0.3.隨著擺動(dòng)頻率增大,臨界雷諾數(shù)減小,內(nèi)部水流更加容易失穩(wěn);與槽道擺動(dòng)向下游傳播(ωc>0)相比,槽道擺動(dòng)向上游傳播(ωc<0)時(shí)水流更加穩(wěn)定.

    圖12 臨界雷諾數(shù)Recr隨擺動(dòng)頻率ωc的變化(平面狀態(tài)參數(shù):θm=0.1,αc=0.3)Fig.12 Variation of critical Reynolds numberRecrwith swinging frequency ωc(for:θm=0.1,αc=0.3)

    圖13為中性曲線隨擺動(dòng)頻率ωc的變化關(guān)系,對(duì)應(yīng)平面狀態(tài)參數(shù):θm=0.1,αc=0.3,-0.02≤ωc≤0.02).圖中繪出了ωc=-0.015,0,0.015時(shí)的中性曲線.由圖13可以看出,隨擺動(dòng)頻率ωc增大,中性曲線左移,臨界點(diǎn)(Recr,αTcr)向左上方移動(dòng),即臨界雷諾數(shù)Recr減小,臨界擾動(dòng)波數(shù)αTcr增大,擾動(dòng)頻率ωTcr減小,該變化趨勢(shì)與圖11和圖12對(duì)應(yīng)一致.槽道擺動(dòng)向上游傳播時(shí),中性曲線隨擺動(dòng)頻率增大的速度大于向下游傳播時(shí)減小的速度.

    圖13 中性曲線隨擺動(dòng)頻率ωc的變化(平面狀態(tài)參數(shù):θm=0.1,αc=0.3,-0.02≤ωc≤0.02)Fig.13 Variation of neutral curve with swinging frequency ωc(for:θm=0.1,αc=0.3,-0.02≤ωc≤0.02)

    4 結(jié)論

    本文建立了隨體坐標(biāo)系下槽道內(nèi)水流流動(dòng)的理論模型,討論了正弦派生型擺動(dòng)邊界下的槽道水流動(dòng)力穩(wěn)定性特征,確定了擺動(dòng)波數(shù)與擺動(dòng)頻率對(duì)主流流速和內(nèi)部擾動(dòng)發(fā)展影響的定量關(guān)系,得到了槽道彎曲度和擺動(dòng)特征對(duì)其內(nèi)部水流不同尺度擾動(dòng)的選擇性影響.具體如下:

    (3)一階流速分布對(duì)擺動(dòng)頻率ωc的影響相當(dāng)敏感,由ωc=0向正負(fù)兩個(gè)方向增大時(shí)流速分布迅速變化,且開始出現(xiàn)與ωc=0時(shí)相反的分布規(guī)律;當(dāng)ωc<-0.05或ωc>0.1時(shí)反向分布占據(jù)主導(dǎo)地位.

    隨著擺動(dòng)頻率ωc由零減小,一階流速、相位分布均迅速反向,且在ωc<-0.05以后不再有明顯變化;隨著擺動(dòng)頻率ωc由零增大,一階流速首先在兩壁附近開始反向,ωc>0.1后反向分布占據(jù)主導(dǎo)地位,而流速相位則迅速反向,槽道中心處相位隨擺動(dòng)頻率ωc變化及其明顯.

    (4)槽道彎曲對(duì)其內(nèi)部水流各種尺度的擾動(dòng)有選擇性影響.擺動(dòng)波數(shù)αcm在(0,0.48)區(qū)間上,槽道擺動(dòng)加速擾動(dòng)增長(zhǎng),水流比直槽道水流更不穩(wěn)定.當(dāng)擺動(dòng)波數(shù)αcm=0.3,擾動(dòng)增長(zhǎng)最快,水流最不穩(wěn)定,此時(shí)臨界雷諾數(shù)不足直槽流的一半.大擺動(dòng)波數(shù)下中性曲線更加平坦.各中性曲線臨界擾動(dòng)對(duì)應(yīng)的點(diǎn)列(Recr,αTcr)和(Recr,ωTcr)分別呈現(xiàn)“繩套”和“拇指”形狀.

    (5)槽道擺動(dòng)對(duì)水流內(nèi)部各種尺度的擾動(dòng)有較大影響.擬序擾動(dòng)發(fā)展對(duì)擺動(dòng)頻率ωc變化十分敏感.隨著擺動(dòng)頻率ωc增大,擬序擾動(dòng)的增長(zhǎng)率增大,臨界雷諾數(shù)減小,水流更加容易失穩(wěn).槽道擺動(dòng)向上游傳播時(shí)水流更加穩(wěn)定.

    1 Davies C,Carpenter PW.Instabilities in a plane channel fl w be-tween compliant walls.Journal of Fluid Mechanics,1997,352: 205-243

    2 Davies C,Carpenter PW.Numerical simulation of the evolution of Tollmien–Schlichting waves over finit compliant panels.Journal of Fluid Mechanics,1997,335:361-392

    3 Hell M,Waters SL.Transverse fl ws in rapidly oscillating elastic cylindrical shells.Journal of Fluid Mechanics,2006,547:185-214

    4 Guaus A,Airiau C,Bottaro A,et al.E ff ects of wall compliance on the linear stability of Taylor–Couette fl w.Journal of Fluid Mechanics,2009,630:331-365

    5 Pitman MW,Lucey AD.Stability of plane-Poiseuille fl w interacting with a finit compliant panel.17th Australasian Fluid Mechanics Conference,Auckland,2010-12-5-9.New Zealand,2010

    6 Hall P.The stability of the Poiseuille fl w modulated at high frequencies//Proceedings of the Royal Society of London,London, 1975.England,1975.453-464

    7 Thomas C,Bassom AP,Blennerhassett PJ,et al.The linear stability of oscillatory Poiseuille fl w in channels and pipes//Proceedings of the Royal Society A,2011.Australia,2011.467:2643-2662

    8 Kuhlmann C,Wanschura M,Rath HJ.Flow in two-sided lid-driven cavities:non-uniqueness,instabilities,and cellular structures.Journal of Fluid Mechanics,1997,336:267-299

    9 Ren L,Chen JG,Zhu KQ.Dual role of wall slip on linear stability of plane Poiseuille fl w.Chinese Physics Letters,2008,25(2): 601-603

    10 Sen PK,Carpenter PW,Hegde S,et al.A wave driver theory for vortical waves propagating across junctions with application to those between rigid and compliant walls.Journal of Fluid Mechanics, 2009,625:1-46

    11 Friedkin JF.A laboratory study of the meandering of alluvial rivers. U.S.Army Engineer Waterways Experiment Station,Vicksburg, Mississippi,1945.USA

    12 Brice JC.Platform properties of meandering processes,In:River Meandering//Elliott CM.ed.Proceedings of the conference Rivers, Louisiana,1983-10-24-26.New Orleans:ASCE,1984.1-15

    13 Hooke JM,Redmond CE.Use of cartographic sources for analyzing river channel change with examples from Britain//Petts GE,ed. Historical of Large Alluvial Rivers:Western Europe,Wiley,1989

    14 Parker G,Sawai K,Ikeda S.Bend theory of river meanders.Part 2. Nonlinear deformation of finite-amplitud bends.Journal of Fluid Mechanics,1982,115:303-314

    15 Ikeda S,Parker G,Sawai K.Bend theory of river meanders.Part 1.Linear development.Journal of Fluid Mechanics,1981,112: 363-377

    16 Bai YC,Xu HJ.A study on the stability of laminar open-channel fl w over a sandy rippled bed.Science in China,Series E:Engineering&Material Science,2005,35:53-73

    17 Xu HJ,Bai YH.Theoretical analyses on hydrodynamic instability in narrow deep river with variable curvature.Applied Mathematics and Mechanics,2015,36(9):1147-1168

    18 Bai YC,Ji ZQ,Xu HJ.Instability and self-adaption character of turbulence coherent structure in narrow-deep river bend.Science China:Technology Science,2012,55:2990-2999

    19 HookeJM.Magnitude and distributionof rates of river bank erosion.Earth Surface Processes and Landforms,1980,5:143-157

    20 Hickin,E J.Lateral migration rates of rivers bends//Cheremisino ff PN,Chereminiso ffNP,Cheng SL eds.Handbook of Civil Engineering.Lancaster,Pennsylvinia:Technomic Publishing,1988

    21 Pyle CJ,Richards KS,Chandler JH.Digital photogrammetric monitoring of river bank erosion.Photogrammetric Record,1997, 15(89):753-764

    22 Luo JS,Wu XS.On the linear instability of a finit stokes layer: Instantaneous versus Floquet modes.Physics of Fluids,2010,22: 054106-1-054106-13

    23 Herbert T.Secondary of boundary layers.Annual Review of Fluid Mechanics,1988,20:487-526

    24 Zhang ZS,Lilley GM.A theoretical model of coherent structure in a plate turbulent boundary layer//Turbulent Shear Flow III.Berlin: Springer-Verlag,1981

    25 張兆順.湍流.北京:國(guó)防工業(yè)出版社,2002(Zhang Zhaoshun. Turbulence.Beijing:National DefenceIndustryPress,2002(inChinese))

    26 Orszag SA.Accurate solution of the Orr-Sommerfeld stability equation.Journal Fluid Mechanics,1971,50(4):689-703

    27 Reynolds WC,Potter MC.Finite-amplitude instability of parallel shear fl ws.Journal of Fluid Mechanics,1967,27:465-492

    28 Nishioka M,Ichikawa Y.An experimental investigation of the stability of plane Poiseuille fl w.Journal of Fluid Mechanics,1975, 72:731-751

    29 CrosatoA.AnalysisandModellingofRiverMeandering.[PhDThesis].Italia:University of Padua,2008

    30 Seminara G,Zolezzi G,Tubino M,et al.Downstream and upstream influenc in river meandering.Part 2.Planimetric development.Journal of Fluid Mechanics,2001,438:213-230

    附錄

    式(14)~式(16)中各參數(shù)為:

    HYDRODYNAMIC INSTABILITY CHARACTERISTICS OF LAMINAR FLOW IN A MEANDERING CHANNEL WITH MOVING BOUNDARY1)

    Bai Yuchuan?,?Ji Ziqing?Xu Haijue?,?,2)

    ?(State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin300072,China)

    ?(Department of Civil Engineering,Tianjin University,Tianjin300072,China)

    Configuratio of river is closely related with hydrodynamic structures of fl ws,for the shape o f a channel influence the fl w structures in it,and the fl w structures also a ff ect the developing trend of the channel through the movement of sediment,forming a pair of dialectical interactions in the river system.The natural rivers are di ff erent in configurations which can generally be divided into such types as straight,bending,branching and wandering.Among them,the bending river or the river composed of several curved channels,the result configuratio of interaction between the natural river and the complex hydrodynamic fl w structure in it,become one of the most important types in the study of river dynamics.As the basis of theoretical research,the establishment of model and the study on fl ws within a moving channel has become the focus not only from researchers of flui mechanics,but also from investigators of river dynamics. Therefore,this study firs established a theoretic model on the fl w in a meandering channel with a moving boundaryby using a streamwise-transverse coordinate system.It next discussed the hydro-dynamic instability characteristics of the laminar fl w within the sine-generated moving boundaries.Then it quantitatively analyzed the influence of various character parameters to the velocity distributions of main fl w.Finally,it obtained the selective influence from the curvature and meandering properties to the fl w structure.

    moving channel,fl w in channel,sine-generated curveform,fl w properties,hydrodynamic instability characteristics

    O352

    A

    10.6052/0459-1879-16-105

    2016–04–18收稿,2017–02–22錄用,2017–02–22網(wǎng)絡(luò)版發(fā)表.

    1)國(guó)家自然科學(xué)基金項(xiàng)目(41576093,51279124,51321065)和天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室基金(HESS-1606)資助.

    2)徐海玨,副教授,主要研究方向:流體力學(xué)、河流海岸動(dòng)力學(xué)、泥沙運(yùn)動(dòng)力學(xué).E-mail:xiaoxiaoxu_2004@163.com

    白玉川,冀自青,徐海玨.擺動(dòng)河槽水動(dòng)力穩(wěn)定性特征分析.力學(xué)學(xué)報(bào),2017,49(2):274-288

    Bai Yuchuan,Ji Ziqing,Xu Haijue.Hydrodynamic instability characteristics of laminar fl w in a meandering channel with moving boundary.Chinese Journal of Theoretical and Applied Mechanics,2017,49(2):274-288

    猜你喜歡
    波數(shù)雷諾數(shù)水流
    聲場(chǎng)波數(shù)積分截?cái)嗖〝?shù)自適應(yīng)選取方法
    一種基于SOM神經(jīng)網(wǎng)絡(luò)中藥材分類識(shí)別系統(tǒng)
    哪股水流噴得更遠(yuǎn)
    能俘獲光的水流
    我只知身在水中,不覺水流
    文苑(2020年6期)2020-06-22 08:41:56
    基于Transition SST模型的高雷諾數(shù)圓柱繞流數(shù)值研究
    失穩(wěn)初期的低雷諾數(shù)圓柱繞流POD-Galerkin 建模方法研究
    基于轉(zhuǎn)捩模型的低雷諾數(shù)翼型優(yōu)化設(shè)計(jì)研究
    民機(jī)高速風(fēng)洞試驗(yàn)的阻力雷諾數(shù)效應(yīng)修正
    重磁異常解釋的歸一化局部波數(shù)法
    99精品在免费线老司机午夜| 在线a可以看的网站| 中文字幕人妻熟人妻熟丝袜美 | 高清日韩中文字幕在线| 香蕉丝袜av| 熟女人妻精品中文字幕| 九色国产91popny在线| 久久久久国产精品人妻aⅴ院| 日韩精品青青久久久久久| 国产伦精品一区二区三区视频9 | 好看av亚洲va欧美ⅴa在| 悠悠久久av| 亚洲激情在线av| 成人午夜高清在线视频| 欧美一级毛片孕妇| 757午夜福利合集在线观看| 免费看a级黄色片| 亚洲国产欧美人成| 美女 人体艺术 gogo| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 老司机福利观看| 男女视频在线观看网站免费| www.999成人在线观看| 一进一出好大好爽视频| 日韩欧美三级三区| 久久九九热精品免费| 色尼玛亚洲综合影院| 一级毛片高清免费大全| ponron亚洲| 免费高清视频大片| 别揉我奶头~嗯~啊~动态视频| 亚洲国产精品合色在线| 色噜噜av男人的天堂激情| av在线天堂中文字幕| 无人区码免费观看不卡| 9191精品国产免费久久| 午夜免费成人在线视频| 国产黄色小视频在线观看| 国产中年淑女户外野战色| av专区在线播放| 三级国产精品欧美在线观看| 久久中文看片网| 99久久无色码亚洲精品果冻| 人妻久久中文字幕网| 欧洲精品卡2卡3卡4卡5卡区| 国产精品,欧美在线| 不卡一级毛片| 精品国产超薄肉色丝袜足j| 国产高清激情床上av| 日本免费a在线| 亚洲欧美激情综合另类| 欧美绝顶高潮抽搐喷水| 最好的美女福利视频网| 99热只有精品国产| 久久久久久久久大av| 中出人妻视频一区二区| 国产亚洲欧美98| 国产亚洲欧美98| 欧美一级a爱片免费观看看| 老熟妇乱子伦视频在线观看| 欧美日韩国产亚洲二区| 99久久九九国产精品国产免费| 一级黄色大片毛片| 不卡一级毛片| 亚洲欧美日韩东京热| 天堂网av新在线| 久久亚洲精品不卡| 亚洲18禁久久av| 一夜夜www| 一本久久中文字幕| 麻豆国产av国片精品| 国产成+人综合+亚洲专区| 可以在线观看毛片的网站| 看片在线看免费视频| 看黄色毛片网站| 成人国产一区最新在线观看| 舔av片在线| xxxwww97欧美| 国产精品亚洲美女久久久| 村上凉子中文字幕在线| 一个人观看的视频www高清免费观看| 国产97色在线日韩免费| 亚洲在线自拍视频| 久久精品国产综合久久久| 啦啦啦韩国在线观看视频| 9191精品国产免费久久| 99久久99久久久精品蜜桃| 无限看片的www在线观看| 老司机在亚洲福利影院| 色吧在线观看| 国产麻豆成人av免费视频| 麻豆一二三区av精品| 神马国产精品三级电影在线观看| 国产伦人伦偷精品视频| 精品久久久久久久久久久久久| 在线a可以看的网站| 精品乱码久久久久久99久播| 午夜精品在线福利| 中文资源天堂在线| 国内揄拍国产精品人妻在线| 午夜精品在线福利| 国产精品电影一区二区三区| 亚洲欧美日韩无卡精品| 国产精品 国内视频| www.熟女人妻精品国产| 每晚都被弄得嗷嗷叫到高潮| 亚洲国产欧美人成| 特大巨黑吊av在线直播| 全区人妻精品视频| 90打野战视频偷拍视频| 欧美日韩亚洲国产一区二区在线观看| 国产亚洲精品久久久久久毛片| 午夜福利高清视频| 欧美性猛交╳xxx乱大交人| 国产淫片久久久久久久久 | 麻豆国产av国片精品| 麻豆一二三区av精品| 免费看a级黄色片| 69人妻影院| 亚洲人成网站在线播| 在线观看66精品国产| 黄片小视频在线播放| 日韩欧美精品免费久久 | 身体一侧抽搐| 国产三级黄色录像| 夜夜爽天天搞| 亚洲欧美一区二区三区黑人| 色在线成人网| 中文字幕人成人乱码亚洲影| 草草在线视频免费看| 久久久久久大精品| 亚洲欧美一区二区三区黑人| 国产亚洲av嫩草精品影院| 欧美日韩福利视频一区二区| 成人精品一区二区免费| 男人舔女人下体高潮全视频| 久久久色成人| 无限看片的www在线观看| 在线观看一区二区三区| 久久久国产成人免费| 九色国产91popny在线| 网址你懂的国产日韩在线| 久久精品国产亚洲av涩爱 | 国产毛片a区久久久久| 亚洲精品在线美女| 国产欧美日韩一区二区三| 精品久久久久久久毛片微露脸| 久久久久九九精品影院| 成人性生交大片免费视频hd| 国产av麻豆久久久久久久| 久久精品人妻少妇| 欧美性猛交黑人性爽| 午夜精品在线福利| 国产综合懂色| www日本在线高清视频| 久久久久亚洲av毛片大全| 久久久国产精品麻豆| 非洲黑人性xxxx精品又粗又长| 欧美另类亚洲清纯唯美| 国产伦精品一区二区三区四那| 欧美日韩综合久久久久久 | 在线播放国产精品三级| 网址你懂的国产日韩在线| 久久久久久久久久黄片| 51国产日韩欧美| 伊人久久精品亚洲午夜| 免费高清视频大片| 亚洲avbb在线观看| 国内少妇人妻偷人精品xxx网站| 叶爱在线成人免费视频播放| 亚洲国产中文字幕在线视频| 国产乱人视频| 国产精品久久视频播放| 国产精品野战在线观看| 久久婷婷人人爽人人干人人爱| 99在线视频只有这里精品首页| 日韩欧美在线乱码| 高清日韩中文字幕在线| 中文字幕精品亚洲无线码一区| 免费看光身美女| 国产精品免费一区二区三区在线| 色视频www国产| 国语自产精品视频在线第100页| 五月玫瑰六月丁香| 一卡2卡三卡四卡精品乱码亚洲| 熟女电影av网| 欧美日韩乱码在线| 高潮久久久久久久久久久不卡| 女人高潮潮喷娇喘18禁视频| 最近最新免费中文字幕在线| 亚洲无线观看免费| 亚洲国产欧洲综合997久久,| 夜夜夜夜夜久久久久| 国产精品98久久久久久宅男小说| 最近视频中文字幕2019在线8| 日本三级黄在线观看| 黄色成人免费大全| 国产三级黄色录像| 久久久久精品国产欧美久久久| 国产午夜福利久久久久久| 久久久久久久久大av| 日韩欧美精品免费久久 | 观看美女的网站| 成年女人毛片免费观看观看9| 黄色日韩在线| 在线免费观看的www视频| 亚洲五月婷婷丁香| 中文字幕av在线有码专区| 男女午夜视频在线观看| 欧美性感艳星| svipshipincom国产片| 亚洲最大成人中文| www.色视频.com| 男插女下体视频免费在线播放| 九色国产91popny在线| 亚洲专区中文字幕在线| 男女视频在线观看网站免费| 亚洲中文日韩欧美视频| 村上凉子中文字幕在线| 美女高潮喷水抽搐中文字幕| 看免费av毛片| 法律面前人人平等表现在哪些方面| 级片在线观看| www.熟女人妻精品国产| 日日夜夜操网爽| 日韩免费av在线播放| 欧美一区二区亚洲| 最新美女视频免费是黄的| ponron亚洲| 久久久久久久精品吃奶| 美女被艹到高潮喷水动态| 亚洲国产精品999在线| 日本一二三区视频观看| 日本熟妇午夜| 久久精品国产清高在天天线| 狂野欧美白嫩少妇大欣赏| 在线十欧美十亚洲十日本专区| 18禁国产床啪视频网站| 亚洲精品456在线播放app | 午夜精品在线福利| www日本在线高清视频| 床上黄色一级片| 校园春色视频在线观看| 欧美一区二区精品小视频在线| 久久亚洲真实| 无限看片的www在线观看| 亚洲国产欧洲综合997久久,| 天堂影院成人在线观看| avwww免费| x7x7x7水蜜桃| 国产高清视频在线播放一区| 18禁裸乳无遮挡免费网站照片| 中文资源天堂在线| 久久精品综合一区二区三区| 欧美一级a爱片免费观看看| eeuss影院久久| 别揉我奶头~嗯~啊~动态视频| 国产乱人伦免费视频| 欧美日韩中文字幕国产精品一区二区三区| 国产色爽女视频免费观看| 搡老岳熟女国产| 亚洲人成电影免费在线| 亚洲人成网站在线播放欧美日韩| 99国产精品一区二区蜜桃av| 国产成+人综合+亚洲专区| 老汉色av国产亚洲站长工具| 男女那种视频在线观看| 搡老岳熟女国产| 一进一出抽搐动态| 国产99白浆流出| 欧美激情在线99| 日韩欧美在线乱码| 国产精品日韩av在线免费观看| 欧美国产日韩亚洲一区| av国产免费在线观看| 精品熟女少妇八av免费久了| 夜夜看夜夜爽夜夜摸| 亚洲欧美一区二区三区黑人| 午夜免费男女啪啪视频观看 | 亚洲av日韩精品久久久久久密| www.熟女人妻精品国产| 日本 av在线| 成年免费大片在线观看| 91久久精品电影网| 禁无遮挡网站| 亚洲av五月六月丁香网| 我要搜黄色片| 午夜免费男女啪啪视频观看 | 亚洲国产日韩欧美精品在线观看 | 欧美乱码精品一区二区三区| 女生性感内裤真人,穿戴方法视频| 午夜久久久久精精品| 国产av在哪里看| 久久6这里有精品| 两人在一起打扑克的视频| 他把我摸到了高潮在线观看| 女同久久另类99精品国产91| 欧美黄色片欧美黄色片| 18禁黄网站禁片午夜丰满| 国产伦在线观看视频一区| 99国产极品粉嫩在线观看| 男女之事视频高清在线观看| 亚洲精品456在线播放app | 男人的好看免费观看在线视频| 757午夜福利合集在线观看| 国产精品爽爽va在线观看网站| 夜夜躁狠狠躁天天躁| 国产精品久久久人人做人人爽| 免费看日本二区| 男人舔奶头视频| 日韩亚洲欧美综合| 精品一区二区三区人妻视频| 一区二区三区免费毛片| av天堂中文字幕网| 亚洲精品影视一区二区三区av| 特级一级黄色大片| 免费大片18禁| 好看av亚洲va欧美ⅴa在| 国产精品日韩av在线免费观看| 亚洲av成人精品一区久久| 婷婷六月久久综合丁香| 国产精品 欧美亚洲| 精华霜和精华液先用哪个| 夜夜夜夜夜久久久久| 午夜福利高清视频| 国产成人aa在线观看| 国产一级毛片七仙女欲春2| 最近最新免费中文字幕在线| 变态另类成人亚洲欧美熟女| 亚洲五月天丁香| 久久精品亚洲精品国产色婷小说| 国产蜜桃级精品一区二区三区| 在线观看日韩欧美| 在线十欧美十亚洲十日本专区| 天堂√8在线中文| 日本三级黄在线观看| 日韩欧美 国产精品| 啪啪无遮挡十八禁网站| 18+在线观看网站| 精品乱码久久久久久99久播| 美女高潮的动态| 亚洲成人免费电影在线观看| 深夜精品福利| 精品无人区乱码1区二区| 天堂网av新在线| 亚洲美女视频黄频| 国产真实乱freesex| 国产伦在线观看视频一区| 久久中文看片网| 亚洲成人久久爱视频| 亚洲最大成人中文| 国产精品精品国产色婷婷| 给我免费播放毛片高清在线观看| 757午夜福利合集在线观看| 超碰av人人做人人爽久久 | 久久中文看片网| 午夜久久久久精精品| www.色视频.com| 亚洲男人的天堂狠狠| 精华霜和精华液先用哪个| 欧美黄色淫秽网站| 美女免费视频网站| 欧美成人免费av一区二区三区| 在线视频色国产色| 色综合欧美亚洲国产小说| 亚洲精品成人久久久久久| 成年女人毛片免费观看观看9| 亚洲成人精品中文字幕电影| www国产在线视频色| 国产野战对白在线观看| www国产在线视频色| 亚洲欧美日韩无卡精品| 成人亚洲精品av一区二区| 国产真人三级小视频在线观看| 欧美一级毛片孕妇| 男女做爰动态图高潮gif福利片| 国产成人啪精品午夜网站| 欧美bdsm另类| 国产精品美女特级片免费视频播放器| 国语自产精品视频在线第100页| 日韩大尺度精品在线看网址| 精品国产亚洲在线| 91麻豆av在线| 精品一区二区三区视频在线观看免费| 淫秽高清视频在线观看| 蜜桃久久精品国产亚洲av| 亚洲电影在线观看av| 国产精品三级大全| 欧美中文综合在线视频| 99久久久亚洲精品蜜臀av| 亚洲黑人精品在线| 欧美精品啪啪一区二区三区| 69av精品久久久久久| 欧美精品啪啪一区二区三区| 最近最新免费中文字幕在线| 成年女人看的毛片在线观看| 一级作爱视频免费观看| 婷婷精品国产亚洲av在线| xxxwww97欧美| 国产一区在线观看成人免费| 欧洲精品卡2卡3卡4卡5卡区| 在线免费观看的www视频| 精品一区二区三区av网在线观看| 国产一区二区激情短视频| 国产伦精品一区二区三区四那| 欧美一区二区国产精品久久精品| 搡老熟女国产l中国老女人| 免费大片18禁| 日韩精品中文字幕看吧| 老鸭窝网址在线观看| 精品不卡国产一区二区三区| 最近在线观看免费完整版| 国产免费一级a男人的天堂| 一级作爱视频免费观看| 99国产极品粉嫩在线观看| 亚洲一区高清亚洲精品| 亚洲片人在线观看| 两个人看的免费小视频| 国产三级黄色录像| 无限看片的www在线观看| 最近最新中文字幕大全免费视频| 麻豆成人av在线观看| 久久性视频一级片| 精品电影一区二区在线| 成人鲁丝片一二三区免费| 美女大奶头视频| 久久久国产成人精品二区| 搡老岳熟女国产| 内地一区二区视频在线| 高清日韩中文字幕在线| 日日夜夜操网爽| 中文字幕熟女人妻在线| 国内精品一区二区在线观看| 99久久无色码亚洲精品果冻| 精品国产超薄肉色丝袜足j| 一级作爱视频免费观看| 国产高清videossex| 久久久国产成人精品二区| 久久精品国产99精品国产亚洲性色| av在线蜜桃| av视频在线观看入口| 国产精品亚洲美女久久久| 99国产综合亚洲精品| 亚洲第一电影网av| 国产一级毛片七仙女欲春2| 性色av乱码一区二区三区2| 少妇的逼水好多| 午夜免费观看网址| 欧美日韩亚洲国产一区二区在线观看| 国产三级在线视频| 好男人在线观看高清免费视频| 国产黄片美女视频| 成人午夜高清在线视频| www日本在线高清视频| 亚洲av熟女| 老熟妇乱子伦视频在线观看| 免费看美女性在线毛片视频| 人人妻人人看人人澡| 成年女人永久免费观看视频| 制服人妻中文乱码| 一级a爱片免费观看的视频| 国产美女午夜福利| 日本成人三级电影网站| 有码 亚洲区| 亚洲午夜理论影院| 女人十人毛片免费观看3o分钟| 三级男女做爰猛烈吃奶摸视频| 国产亚洲av嫩草精品影院| 91字幕亚洲| 欧美中文日本在线观看视频| 久久99热这里只有精品18| 欧美一级a爱片免费观看看| a级毛片a级免费在线| 欧美精品啪啪一区二区三区| 亚洲五月婷婷丁香| av天堂在线播放| 亚洲人成网站在线播| 99热只有精品国产| 久久精品国产综合久久久| 69av精品久久久久久| 欧美三级亚洲精品| 亚洲欧美激情综合另类| 婷婷六月久久综合丁香| 成人一区二区视频在线观看| 一夜夜www| 国产精品亚洲一级av第二区| bbb黄色大片| 国产一区在线观看成人免费| 国产成人av激情在线播放| 精品久久久久久久毛片微露脸| 色老头精品视频在线观看| 毛片女人毛片| 婷婷精品国产亚洲av在线| 亚洲avbb在线观看| 哪里可以看免费的av片| 99国产精品一区二区三区| 久久久久九九精品影院| 99国产极品粉嫩在线观看| 国产精品久久电影中文字幕| 国产欧美日韩精品亚洲av| 国产午夜精品论理片| 亚洲五月天丁香| 久久久久久人人人人人| 又爽又黄无遮挡网站| 国产91精品成人一区二区三区| 国产激情偷乱视频一区二区| 国产高清有码在线观看视频| 三级国产精品欧美在线观看| 久久精品国产99精品国产亚洲性色| 国产精品一及| 欧美日本视频| 免费看十八禁软件| 97碰自拍视频| 日韩欧美在线二视频| 亚洲内射少妇av| 亚洲欧美日韩东京热| 亚洲最大成人手机在线| 国产av在哪里看| 久久久久精品国产欧美久久久| 国产一区二区亚洲精品在线观看| 欧美不卡视频在线免费观看| 日韩有码中文字幕| 一区福利在线观看| 国产精品野战在线观看| 色播亚洲综合网| 一进一出抽搐gif免费好疼| 又爽又黄无遮挡网站| 精品乱码久久久久久99久播| 一本综合久久免费| 禁无遮挡网站| 欧美av亚洲av综合av国产av| 天堂网av新在线| 亚洲av第一区精品v没综合| 精品国产超薄肉色丝袜足j| 日本黄色片子视频| 久久久久久久久久黄片| 国产免费男女视频| 国产精品嫩草影院av在线观看 | 婷婷丁香在线五月| 波多野结衣高清作品| 99国产精品一区二区三区| 男人和女人高潮做爰伦理| 国产成人系列免费观看| 午夜精品一区二区三区免费看| 日本a在线网址| 999久久久精品免费观看国产| 国产单亲对白刺激| 亚洲熟妇中文字幕五十中出| 老司机深夜福利视频在线观看| 成年人黄色毛片网站| 美女大奶头视频| 人妻丰满熟妇av一区二区三区| 69人妻影院| 成人特级黄色片久久久久久久| 欧美3d第一页| 国产成年人精品一区二区| 熟妇人妻久久中文字幕3abv| 欧美成狂野欧美在线观看| 欧美av亚洲av综合av国产av| 看免费av毛片| 一个人观看的视频www高清免费观看| 中文字幕av成人在线电影| 高清在线国产一区| 俄罗斯特黄特色一大片| 精品国产超薄肉色丝袜足j| 成年女人毛片免费观看观看9| 亚洲熟妇中文字幕五十中出| 国产毛片a区久久久久| 国产美女午夜福利| 国产精品久久久久久精品电影| 在线国产一区二区在线| 俺也久久电影网| 日韩欧美国产一区二区入口| 精品人妻一区二区三区麻豆 | 丁香六月欧美| 天天躁日日操中文字幕| 国产av在哪里看| 国产成人a区在线观看| 叶爱在线成人免费视频播放| 亚洲中文字幕一区二区三区有码在线看| 欧美一级毛片孕妇| 久久久久久国产a免费观看| 久久国产精品影院| 免费在线观看亚洲国产| 亚洲精品成人久久久久久| 中文在线观看免费www的网站| 九九热线精品视视频播放| 亚洲成人免费电影在线观看| 蜜桃亚洲精品一区二区三区| 高清日韩中文字幕在线| 天堂av国产一区二区熟女人妻| 天堂√8在线中文| 日韩高清综合在线| 小说图片视频综合网站| 毛片女人毛片| 老汉色∧v一级毛片| 可以在线观看的亚洲视频| 99国产极品粉嫩在线观看| aaaaa片日本免费| 国产主播在线观看一区二区| 99精品在免费线老司机午夜| 国产一区二区激情短视频| 91九色精品人成在线观看| 99国产精品一区二区蜜桃av| 一级黄色大片毛片| 成人av在线播放网站| 免费看a级黄色片| 欧美极品一区二区三区四区| 午夜福利免费观看在线| 成年女人看的毛片在线观看| 欧美另类亚洲清纯唯美| 国产精品自产拍在线观看55亚洲| 老汉色∧v一级毛片| 国产av一区在线观看免费| 18禁美女被吸乳视频|