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

    多層速度格子Boltzmann 方法進(jìn)展及展望

    2022-07-13 01:55:12單肖文
    關(guān)鍵詞:平衡態(tài)等溫流動(dòng)

    楊 鯤,單肖文

    (1. 南方科技大學(xué) 前沿與交叉科學(xué)研究院,深圳 518055;2. 南方科技大學(xué) 工學(xué)院 力學(xué)與航空航天系,深圳 518055)

    0 引 言

    誕生于格子氣自動(dòng)機(jī)[1]的格子Boltzmann 方法(LBM),與傳統(tǒng)直接求解宏觀量的計(jì)算流體力學(xué)(CFD)方法相比,由于避開(kāi)了Navier-Stokes 方程非線性項(xiàng)的處理,且沒(méi)有求解壓力泊松方程這一繁重過(guò)程,使得LBM 具有程序編制簡(jiǎn)單、并行計(jì)算效率高的優(yōu)點(diǎn)[2],在湍流、多相流、多孔介質(zhì)流動(dòng)及滲流、顆粒流、微流動(dòng)、氣動(dòng)聲學(xué)等領(lǐng)域的模擬中得到了廣泛應(yīng)用[3-6]。早期的LBM 離散模型都是基于半經(jīng)驗(yàn)理論推導(dǎo),其約束條件往往局限于滿足質(zhì)量和動(dòng)量守恒方程,并以低馬赫數(shù)假設(shè)為基礎(chǔ),能量方程無(wú)法正確還原,使得LBM 方法只適用于等溫低速弱可壓流動(dòng)成為了較普遍的觀點(diǎn)[2,7]。

    隨著航空航天科技發(fā)展,高速流動(dòng)的數(shù)值模擬在飛行器設(shè)計(jì)驗(yàn)證中起到越來(lái)越重要的作用[8],CFD 領(lǐng)域涌現(xiàn)了多種可壓縮流動(dòng)的處理格式與方法,例如黎曼問(wèn)題求解器和TVD 格式[9]、用于高精度激波捕捉的WENO 格式[10],以及計(jì)算激波-湍流相互作用的Compact-WENO 混合格式[11];這些數(shù)值方法在理論研究和工程應(yīng)用領(lǐng)域都取得了可觀的成果[12-13]。不可否認(rèn)的是,盡管傳統(tǒng)CFD 方法已經(jīng)過(guò)了多年的發(fā)展和完善,在具體問(wèn)題中要同時(shí)準(zhǔn)確解析流動(dòng)小尺度結(jié)構(gòu)和激波間斷,實(shí)現(xiàn)模擬方法的精度、耗散、穩(wěn)定性和計(jì)算效率等因素的最佳平衡,仍然存在著諸多挑戰(zhàn)[14]。另一方面,工程傳熱流動(dòng)問(wèn)題中,高效換熱器普遍具有形狀復(fù)雜的管腔結(jié)構(gòu),處理復(fù)雜幾何邊界、并兼顧精度和計(jì)算效率也對(duì)流動(dòng)模擬方法提出了更高的要求。

    LBM 源自于氣體動(dòng)理學(xué)理論,對(duì)氣體運(yùn)動(dòng)描述的物理背景清晰,且在處理復(fù)雜幾何外形、特別是多孔結(jié)構(gòu)時(shí)具有獨(dú)特優(yōu)勢(shì),理論上LBM 在模擬可壓縮流動(dòng)和復(fù)雜結(jié)構(gòu)的傳熱流動(dòng)問(wèn)題方面擁有巨大潛力。研究者們對(duì)傳統(tǒng)LBM 模型進(jìn)行了諸多改進(jìn),以期突破LBM 在計(jì)算可壓縮流動(dòng)和熱力學(xué)流動(dòng)時(shí)的局限性,其基本思想是增加LBM 離散模型的自由度和約束條件以還原能量守恒方程。改進(jìn)方案大致分為兩類,一類是額外增加一個(gè)能量(內(nèi)能或總能)分布函數(shù)以滿足能量方程約束[15-20],稱之為雙分布函數(shù)模型。然而,根據(jù)氣體動(dòng)理學(xué)理論,氣體的密度、動(dòng)量和能量都是氣體分子熱運(yùn)動(dòng)的統(tǒng)計(jì)結(jié)果,其動(dòng)力學(xué)狀態(tài)也應(yīng)當(dāng)只需由單一的分布函數(shù)描述(對(duì)于單組分、單原子氣體而言),因此引入能量的分布函數(shù)必然需要與密度分布函數(shù)進(jìn)行耦合,該耦合關(guān)系帶有一定的經(jīng)驗(yàn)性,在增加復(fù)雜度的同時(shí)也對(duì)模型的廣泛適用性帶來(lái)了一定的限制。另一種方案是構(gòu)造更高階的單一平衡態(tài)分布函數(shù),在平衡態(tài)分布函數(shù)中引入更多的待定系數(shù),并使用更多的離散速度點(diǎn),通過(guò)增加離散點(diǎn)數(shù)量(未知數(shù)個(gè)數(shù))的方案來(lái)滿足額外的能量方程約束,典型做法是在每個(gè)離散速度方向上使用2 個(gè)以上的速度大小模態(tài),稱之為多層速度模型。與雙分布函數(shù)模型類似,多層速度模型在能量方程約束條件、分布函數(shù)的形式和離散速度點(diǎn)選取上也具有一定的經(jīng)驗(yàn)性,自Qian[21]和Alexander 等[22]提出以來(lái),出現(xiàn)了形式多樣的方案且各有優(yōu)缺點(diǎn)。

    多層速度模型按平衡態(tài)分布函數(shù)展開(kāi)式的構(gòu)造方法主要分為三類:第一類遵循傳統(tǒng)LBM 模型的思路,將Maxwell-Boltzmann 分布函數(shù)泰勒展開(kāi)至更高階數(shù)[23-26]。第二類由Shan 等[27]提出,他們發(fā)展了Grad[28-29]的早期工作思路,認(rèn)為平衡態(tài)分布函數(shù)可由Hermite 多項(xiàng)式級(jí)數(shù)展開(kāi),只需要將展開(kāi)級(jí)數(shù)保留足夠的截?cái)嚯A數(shù),即可分別還原到Navier-Stokes 方程的質(zhì)量、動(dòng)量和能量方程,而密度、動(dòng)量、能量等宏觀統(tǒng)計(jì)量均可由Hermite 展開(kāi)式的各階系數(shù)組合表達(dá),并在進(jìn)一步的工作中闡明傳統(tǒng)等溫弱可壓LBM是Hermite 多項(xiàng)式模型的低階形態(tài),且離散速度是五階精度Gauss-Hermite 積分公式的積分點(diǎn)的事實(shí)[30]。第三類由Xu 等[31]提出,其基本思想是不預(yù)設(shè)平衡態(tài)分布函數(shù)的展開(kāi)形式,而是將離散模型還原質(zhì)量、動(dòng)量和能量方程的約束條件視為平衡態(tài)分布函數(shù)離散點(diǎn)的線性代數(shù)方程組,直接求解方程組得出離散分布函數(shù)。

    由于Maxwell-Boltzmann 分布函數(shù)的各階Hermite多項(xiàng)式展開(kāi)是確定的,Gauss-Hermite 積分公式的積分點(diǎn)也具有確定性,使得上述第二類基于Hermite 多項(xiàng)式構(gòu)造的多層速度模型不依賴于經(jīng)驗(yàn)而具有普適性。然而,由于該模型的數(shù)學(xué)基礎(chǔ)理論較為晦澀,各部分細(xì)節(jié)的闡述散落于跨度接近20 年的多篇文獻(xiàn)中,使得科研工作者對(duì)該模型進(jìn)行深入理解存在一定的困難。為了彌補(bǔ)這個(gè)缺點(diǎn),使讀者連貫、系統(tǒng)地理解Hermite 多項(xiàng)式模型的構(gòu)造原理和方法,同時(shí)也在與其他LBM 模型的對(duì)比中全面認(rèn)識(shí)各類方案的優(yōu)劣,我們對(duì)Hermite 多項(xiàng)式模型重新梳理歸納的同時(shí),也一并列出了其他多層速度模型的構(gòu)造思想,作為L(zhǎng)BM 多層速度模型的一個(gè)簡(jiǎn)要總結(jié)性工作。

    本文的各部分內(nèi)容為:第1 節(jié)為Boltzmann 方程和LBM 的基礎(chǔ)理論。由于多層速度模型起源于經(jīng)典的等溫弱可壓LBM,且經(jīng)典LBM 模型的很多思想在多層速度模型中仍然適用,故第2 節(jié)對(duì)等溫弱可壓LBM 進(jìn)行簡(jiǎn)單介紹。第3 節(jié)介紹各類多層速度LBM模型(除Hermite 多項(xiàng)式模型以外)的基本思想。第4 節(jié)詳細(xì)描述Hermite 多項(xiàng)式模型的理論和平衡態(tài)分布函數(shù)的具體形式,以及基于Gauss-Hermite 積分公式的離散速度模型及其構(gòu)造方法。第5 節(jié)對(duì)LBM 和傳統(tǒng)CFD 方法的結(jié)合進(jìn)行了簡(jiǎn)要介紹,例如LBM 有限差分、LBM 有限體積和LBM 有限元方法。第6 節(jié)對(duì)現(xiàn)有的多層速度模型存在的問(wèn)題進(jìn)行總結(jié),并提出未來(lái)需要完善和發(fā)展的方向。

    1 Boltzmann-BGK 方程

    由于無(wú)量綱化的Boltzmann-BGK 方程和平衡態(tài)分布函數(shù)分別與式(10)和式(11)形式上完全一致。因此

    2 等溫弱可壓LBM 模型

    等溫弱可壓LBM 模型是最早提出的LBM 模型,并在不可壓流動(dòng)計(jì)算中得到廣泛應(yīng)用。嚴(yán)格意義上,包括氣體、液體在內(nèi)的所有流體都具有可壓縮性,例如液體的壓縮性可用Tait 方程描述[34]。傳統(tǒng)求解壓力泊松方程的不可壓流計(jì)算方法將導(dǎo)致無(wú)限大聲速,也就是壓力擾動(dòng)的無(wú)限速度傳播,這與物理事實(shí)是相悖的。采用弱可壓方法計(jì)算“不可壓流動(dòng)”理論上更加接近物理本質(zhì)。本節(jié)對(duì)等溫弱可壓LBM 模型進(jìn)行簡(jiǎn)要陳述,這也將成為多層速度LBM 模型修正和拓展的基礎(chǔ)。

    2.1 速度空間離散模型

    圖1 等溫弱可壓流動(dòng)LBM 的一維和二維離散速度模型Fig. 1 One- and two-dimensional discrete velocity models for isothermal weakly compressible flows

    2.2 時(shí)間和空間離散

    3 多層速度LBM 模型

    從第2 節(jié)可以看到,等溫弱可壓流動(dòng)LBM 模型,其離散速度在每個(gè)方向都只有一個(gè)值(分布函數(shù)信息只傳播到相鄰的網(wǎng)格點(diǎn)),因此稱之為單層速度模型。推導(dǎo)該模型時(shí),只約束了離散模型的質(zhì)量、動(dòng)量和動(dòng)量輸運(yùn)量,故無(wú)法準(zhǔn)確還原能量方程。受此啟發(fā),業(yè)內(nèi)學(xué)者開(kāi)始對(duì)離散模型施加能量相關(guān)(f(0)的ξ三階以上矩)的約束條件,伴隨而來(lái)的是更高階的f(0)展開(kāi)式和更多的離散速度點(diǎn),其中普遍做法是在各離散速度方向上采用2 個(gè)以上的速度值(模態(tài)),因此這些模型稱之為多層速度模型。本節(jié)對(duì)幾類較為典型的多層速度模型進(jìn)行簡(jiǎn)要陳述。

    3.1 早期多層速度模型

    圖2 Qian 的二維離散速度模型[21]Fig. 2 Two-dimensional discrete velocity model of Qian[21]

    圖3 Alexander 等[22]的二維離散速度模型Fig. 3 Two-dimensional discrete velocity model of Alexander et al[22]

    3.2 Watari-Tsutahara 模型

    3.3 比熱可變模型

    圖4 Chen 等的二維多層速度模型[23]Fig. 4 Two-dimensional discrete velocity model of Chen et al.[23]

    圖5 Watari-Tsutahara 二維多層速度模型[24-25]Fig. 5 Two-dimensional discrete velocity model of Watari-Tsutahara[24-25]

    圖6 Kataoka-Tsutahara[26]二維多層速度模型,圖中的坐標(biāo)代表離散速度點(diǎn)ξ a/cs在第一象限的坐標(biāo)Fig. 6 Two-dimensional discrete velocity model of Kataoka-Tsutahara[26]. Coordinates represent the discrete velocities in the first quadrant

    圖7 離散Boltzmann 方法的幾種二維離散速度模型[39]Fig. 7 Several two-dimensional discrete velocity models in DBM[39]

    4 Hermite 多項(xiàng)式模型

    4.1 速度空間展開(kāi)

    從第3 節(jié)的介紹可以看出,大部分多層速度模型還是沿襲了單層速度模型“預(yù)設(shè)平衡態(tài)分布函數(shù)離散形式—預(yù)設(shè)離散速度模型—待定系數(shù)法求解”的思路。這種思路可以針對(duì)每一類特定問(wèn)題都構(gòu)造相對(duì)優(yōu)化的多層速度模型,但經(jīng)驗(yàn)性較強(qiáng),對(duì)于開(kāi)始接觸多層速度模型的學(xué)者而言,面對(duì)種類繁多的離散模型、平衡態(tài)分布函數(shù)形式和約束條件,如何選擇最優(yōu)方案往往無(wú)所適從;即便獲得了對(duì)某一問(wèn)題的優(yōu)化組合,當(dāng)流動(dòng)參數(shù)改變后,計(jì)算結(jié)果不一定同樣令人滿意,無(wú)法保證模型的普適性。

    Shan 等[27,30]在Grad[28-29]工作的基礎(chǔ)上,發(fā)展了Hermite 多項(xiàng)式模型,該模型基于Hermite 多項(xiàng)式級(jí)數(shù)展開(kāi)的數(shù)學(xué)性質(zhì),使平衡態(tài)分布函數(shù)的展開(kāi)形式和離散速度模型都不再依賴于經(jīng)驗(yàn)而具有確定性,且可以證明單層速度模型事實(shí)上是多層速度模型的低階形態(tài)。這正是本節(jié)將詳述的內(nèi)容。

    參照Grad 的定義,D維空間的n階Hermite 多項(xiàng)式定義為n階張量:

    表1 不同流動(dòng)分布函數(shù)的Hermite 截?cái)嗍阶畹碗A數(shù)Table 1 Lowest order of truncated Hermite distribution function for different types of flows

    4.2 時(shí)間和空間離散

    4.3 離散速度模型

    圖8 2 種D2V17 離散速度模型Fig. 8 Two types of D2V17 discrete velocity models

    在之前的討論中,有一個(gè)重要的問(wèn)題沒(méi)有展開(kāi)——如何獲得具有K階精度的數(shù)值積分公式(式點(diǎn)序列要保留表2 中列出的至少2 個(gè)8 積分點(diǎn)組,其總積分點(diǎn)數(shù)最少為 1+4×5+8×2=37。若保留表2中的前8 組積分點(diǎn)組與權(quán)重,則得出與Pillippi[51]相同的D2V37 模型,如圖9 所示。

    表2 二維積分公式的積分點(diǎn)與權(quán)重Table 2 Quadrature abscissas and weights of the 2D integral formula

    圖9 D2V37 離散速度模型Fig. 9 D2V37 discrete velocity models

    對(duì)于三維情況下的積分公式,也可以遵循同樣的思路求解。相關(guān)的細(xì)節(jié)過(guò)程可參考文獻(xiàn)[47]。

    4.4 等溫弱可壓模型的推導(dǎo)

    5 與傳統(tǒng)CFD 方法的結(jié)合

    經(jīng)典的LBM 中時(shí)間和空間離散均為沿特征線積分(參考2.2 節(jié)和4.2 節(jié)),并設(shè)置歸一化的 Δt=1,這極大簡(jiǎn)化了LBM 的計(jì)算程序。然而這種時(shí)空離散方法也帶來(lái)了幾個(gè)弊端:

    1)由于 Δx為固定值,只能采用均勻正方形(體)網(wǎng)格。對(duì)于邊界層流動(dòng)計(jì)算而言,為了解析邊界層內(nèi)的流動(dòng)而加密全計(jì)算域網(wǎng)格會(huì)極大增加計(jì)算量。在計(jì)算量和解析度之間取平衡的辦法是使用局部加密網(wǎng)格,但是粗網(wǎng)格和細(xì)網(wǎng)格之間的插值會(huì)帶來(lái)新誤差[52]。

    2)無(wú)法使用貼體網(wǎng)格。在傾斜方向邊界,或者曲面邊界上,需要同時(shí)使用笛卡爾網(wǎng)格和浸沒(méi)邊界法[53],以鋸齒形邊界近似斜邊界或者曲面邊界。這種非貼體網(wǎng)格由于精度較低,且在對(duì)流擴(kuò)散問(wèn)題中易造成收斂性困難[54]和穩(wěn)定性問(wèn)題,往往需要與局部加密網(wǎng)格和多松弛時(shí)間模型(MRT)結(jié)合[55]。

    3)計(jì)算穩(wěn)定性。經(jīng)驗(yàn)表明,經(jīng)典LBM 方法在雷諾數(shù)較大時(shí)只能使用較小的物理時(shí)間步長(zhǎng)才能獲得穩(wěn)定的數(shù)值解;當(dāng)使用Zhou-He 非平衡態(tài)反彈邊界條件[56]時(shí)穩(wěn)定性問(wèn)題更加突出,這對(duì)實(shí)際工程湍流的計(jì)算會(huì)形成較大的障礙。

    由于傳統(tǒng)CFD 方法可以使用邊界層加密貼體網(wǎng)格,使用合適的時(shí)間與空間格式時(shí),數(shù)值穩(wěn)定性良好,將傳統(tǒng)CFD 方法與LBM 結(jié)合成為了解決以上問(wèn)題的一種思路。Reider 和Sterling[57]、Cao 等[58]將有限差分法與LBM 模型結(jié)合,提出了FDLBM 方法。FDLBM 的基本思路是將LBM 方程(式(104))視為關(guān)于各個(gè)離散分布函數(shù)fa的雙曲型方程,通過(guò)有限差分法的時(shí)間離散格式,例如Euler 格式或者Runger-Kutta 格式,將 ?fa/?t離散化;并使用迎風(fēng)型或中心型格式將對(duì)流項(xiàng) ξa·?fa進(jìn)行空間離散[59]。在可壓縮問(wèn)題中,需要將多層速度模型和激波捕捉差分格式結(jié)合。Kataoka 和Tsutahara[26]在多層速度模型的基礎(chǔ)上,使用二階迎風(fēng)格式離散對(duì)流項(xiàng);Wang 等[60]將二階TVD 和五階WENO 格式[10]與Kataoka-Tsutahara多層速度模型結(jié)合使用;許愛(ài)國(guó)等[39]則使用了NND 格式[61]來(lái)處理對(duì)流項(xiàng),并在一維和二維黎曼問(wèn)題的計(jì)算中取得較為滿意的結(jié)果。

    使用隱式時(shí)間離散是傳統(tǒng)CFD 方法中提高數(shù)值穩(wěn)定性的常用方法,伴隨而來(lái)的是需要迭代求解線型代數(shù)方程組,這往往使得計(jì)算量和實(shí)施難度顯著增加。在FDLBM 中,Guo 和Zhao[62]使用了與經(jīng)典LBM 類似的隱式時(shí)間處理辦法,即類似于式(30)定義一個(gè)中間時(shí)刻的fˉa代入FDLBM 方程中,可將隱式方法轉(zhuǎn)換為顯式方法。Wang 等[60]沿用類似的思路,將隱式-顯式Runger-Kutta 格式[63]引入FDLBM 并消除了隱式格式需要的迭代求解過(guò)程。

    為了適應(yīng)復(fù)雜幾何構(gòu)型和不規(guī)則形狀網(wǎng)格的計(jì)算需求,Nannell 和Succi[64]、Benzi 等[65]將有限體積方法引入LBM,提出了FVLBM 方法。FVLBM 的基本思想是,將LBM 方程(式(104))對(duì)流場(chǎng)中任意一個(gè)控制體單元 Ω (邊界為 ?Ω ,此處 Ω不代表碰撞模型)進(jìn)行體積分得到:

    則式(123)可寫為關(guān)于分布函數(shù)單元體平均量和單元面對(duì)流通量的方程。求解該方程的過(guò)程即與傳統(tǒng)有限體積方法類似。Chen[66]和Peng 等[67]將FVLBM的理論和其在非均勻網(wǎng)格上的應(yīng)用進(jìn)行了完善。為了提高穩(wěn)定性,Stiebler 等[68]使用迎風(fēng)型格式重構(gòu)單元面通量。Patil 和Lakshmisha[69]則將TVD 格式應(yīng)用于單元面通量的重構(gòu)。

    隨著計(jì)算流體力學(xué)有限元方法日趨成熟[70-71],Lee 和Lin[72]、Shi 等[73]將有限元方法引入LBM,發(fā)展了 有 限 元LBM 方 法。近 年 來(lái),MacMeccan 等[74]、Krüger 等[75]使用混合型有限元LBM 方法模擬了可變形顆粒流動(dòng)并應(yīng)用于生物流計(jì)算領(lǐng)域。由于主題和篇幅所限,此處不對(duì)有限元LBM 方法作進(jìn)一步展開(kāi)。有興趣的讀者可參閱有關(guān)文獻(xiàn)了解相關(guān)細(xì)節(jié)。

    6 總結(jié)與展望

    LBM 多層速度模型起源于單層速度模型,在熱流和可壓縮流動(dòng)應(yīng)用需求下不斷發(fā)展,經(jīng)歷了平衡態(tài)分布函數(shù)展開(kāi)從二階提高到四階、離散速度模型從單層網(wǎng)格提高到2 層或者更多層網(wǎng)格、約束條件也逐步增多的過(guò)程。Hermite 多項(xiàng)式模型則提供了一個(gè)既理論又實(shí)用的判定框架:要想準(zhǔn)確還原能量方程,平衡態(tài)分布函數(shù)必須展開(kāi)到速度四階量。Watari-Tsutahara預(yù)設(shè)的多層速度模型[24-25]與該結(jié)論殊途同歸。Hermite 多項(xiàng)式模型也同時(shí)指出,離散速度模型事實(shí)上為相應(yīng)階數(shù)積分公式的積分點(diǎn)和權(quán)重,只需查表(或用程序包計(jì)算)找出這些積分點(diǎn)和權(quán)重,選出一組應(yīng)用即可,而不是根據(jù)經(jīng)驗(yàn)預(yù)設(shè)和優(yōu)化。盡管Hermite 多項(xiàng)式模型給出了較為簡(jiǎn)潔和確定的形式,但是依照經(jīng)驗(yàn)理論進(jìn)行馬赫展開(kāi)或待定系數(shù)得到的其他多層速度模型,在特定范圍和領(lǐng)域內(nèi)仍具有一定的應(yīng)用價(jià)值。

    必須承認(rèn),相對(duì)于計(jì)算等溫弱可壓流動(dòng)的經(jīng)典單層速度模型而言,多層速度模型仍然是一種尚不成熟的模型,在很多方面仍然需要改進(jìn)和完善。這主要在以下幾個(gè)方面:

    1)模型的誤差、色散與耗散、穩(wěn)定性分析。目前對(duì)單層速度模型的誤差和穩(wěn)定性分析已有不少工作,例如:Sterling 和Chen[76]對(duì)D2Q7(六邊形模型)、D2Q9 和D3Q15 進(jìn)行線性穩(wěn)定性分析后得出 τ >0.5才能穩(wěn)定的結(jié)論。Marie 等[77]對(duì)D3Q19 模型進(jìn)行色散和耗散分析,并與低色散高階差分格式進(jìn)行了比較。Silva 等[78]和Bauer 等[79]對(duì)最常用的D3Q15、D3Q19和D3Q27 模型進(jìn)行了截?cái)嗾`差分析,指出這三種模型都能在低馬赫數(shù)下還原動(dòng)量方程,但是動(dòng)量輸運(yùn)項(xiàng)的高階截?cái)嗾`差不同,這將導(dǎo)致在各向異性流動(dòng)的計(jì)算中出現(xiàn)較大差別;在將平衡態(tài)分布函數(shù)展開(kāi)式根據(jù)連續(xù)平衡態(tài)分布函數(shù)的二階矩約束條件和D3Q19 速度模型進(jìn)行系數(shù)修正后,能獲得與D3Q27 模型相當(dāng)?shù)母飨蛲运健?/p>

    對(duì)多層速度模型而言,McNamara 等[80]對(duì)D3Q21模 型(D3Q15 加 上 (±2c,0,0)、 (0,±2c,0) 和 (0,0,±2c)這6 個(gè)離散點(diǎn))進(jìn)行了時(shí)間離散格式分析,指出使用Lax-Wendroff 時(shí)間格式將在高瑞利數(shù)下具有更好的穩(wěn)定性。Wissocq 等[81]則對(duì)D2Q9 和D2V17 進(jìn)行了模態(tài)和穩(wěn)定性分析,指出D2V17 能準(zhǔn)確解析線性等溫Navier-Stokes 方程的3 個(gè)特征值模態(tài),而D2Q9 則在這些模態(tài)的高波數(shù)上產(chǎn)生一定的誤差。這一點(diǎn)與Hermite 多項(xiàng)式模型的結(jié)論是一致的,即只有使用七階積分公式的積分點(diǎn)(即D2V17 模型)作為離散速度點(diǎn)才能準(zhǔn)確還原動(dòng)量方程;否則由于速度三階量誤差的存在而必須使用低馬赫數(shù)假設(shè)。

    值得注意的是,在Hermite 多項(xiàng)式模型框架下,經(jīng)典D3Q15、D3Q19 和D3Q27 模型從積分精度角度來(lái)講,并不存在差別,因?yàn)槎季哂形咫A積分精度;從計(jì)算經(jīng)濟(jì)性原則來(lái)講,應(yīng)該選擇積分點(diǎn)最少的D3Q15 模型。然而Silva 等[78]和Bauer 等[79]的工作提示我們,在選擇離散模型時(shí)還需要考慮各向同性、截?cái)嗾`差等更多方面的因素,即計(jì)算最經(jīng)濟(jì)的模型不一定是最優(yōu)的計(jì)算結(jié)果。Sieber 等[82]將Hermite 多項(xiàng)式模型與D2Q9 模型、Chen 等[23]的多層速度模型進(jìn)行了線性穩(wěn)定性分析對(duì)比,發(fā)現(xiàn)對(duì)等溫弱可壓流動(dòng)而言,將平衡態(tài)分布函數(shù)展開(kāi)到三階Hermite 多項(xiàng)式級(jí)數(shù),或使用D2V17、D2V37 模型都能顯著提高穩(wěn)定性;對(duì)于可壓縮熱流而言,使用四階Hermite 展開(kāi)的f(0)與D2V37 模型,穩(wěn)定性都優(yōu)于Chen 等[23]提出的多層速度模型。該結(jié)論在Hermite 多項(xiàng)式模型的理論框架下是易于被推斷的。然而,對(duì)于完全還原了能量方程、且離散速度點(diǎn)數(shù)量最少的D2V37、D3V103 模型[47]而言,與同階積分精度、但具有更多離散速度點(diǎn)的模型橫向?qū)Ρ?,例如D3V107 模型,是否也存在各項(xiàng)同性、穩(wěn)定性、色散和耗散差異?這一方面需要在后續(xù)工作中展開(kāi)詳盡地理論分析;另一方面也需要大量的數(shù)值試驗(yàn),例如方管Poiseuille 流動(dòng)、旋轉(zhuǎn)管道流、強(qiáng)可壓縮流動(dòng)、大溫差熱流等實(shí)際流動(dòng)計(jì)算,來(lái)驗(yàn)證這些多層速度模型之間的性質(zhì)差異。

    2)邊界條件。LBM 單層速度模型在經(jīng)過(guò)多年發(fā)展后,其邊界條件處理方法已經(jīng)日趨完善,有關(guān)工作可謂汗牛充棟[33,36,83],對(duì)于壁面反彈邊界條件[56,84]、周期邊界條件[85]、適用于流動(dòng)出入口的無(wú)反射邊界條件[86-87]以及其他類型的邊界條件等都已有充分的研究工作,可滿足絕大多數(shù)計(jì)算工況的需求。然而,當(dāng)使用多層速度模型時(shí),一方面由于離散速度點(diǎn)跨越多層網(wǎng)格,使得邊界附近的內(nèi)部結(jié)點(diǎn)與邊界結(jié)點(diǎn)區(qū)分變得模糊;另一方面平衡態(tài)分布函數(shù)的形式也變得更為復(fù)雜,如何恰當(dāng)處理邊界條件成為了一個(gè)難點(diǎn),相關(guān)的工作也屈指可數(shù)[83]。Malaspinas 等提出[88],將邊界節(jié)點(diǎn)的分布函數(shù)離散點(diǎn)按內(nèi)部和邊界分為已知量和未知量,將離散分布函數(shù)各階矩的約束條件構(gòu)成方程組,通過(guò)求解方程組構(gòu)造一種“通用”的邊界條件處理方法。他們?cè)贒2Q9 和D3Q19 模型中進(jìn)行了測(cè)試,并指出該方法可適用于多層速度模型的邊界條件(但沒(méi)有后續(xù)驗(yàn)證)。Meng 和Zhang[89]提出了適用于多層速度模型的擴(kuò)散-反彈邊界條件格式,并應(yīng)用到D2V16、D2V17、D2V37 和D3V121[45]這四種模型,對(duì)等溫流動(dòng)和熱流進(jìn)行了測(cè)試計(jì)算。Lee 等[90]將Hecht 和Harting 針對(duì)D3Q19 的非平衡態(tài)壁面反彈邊界條件[91],推廣到了用多層速度模型計(jì)算等溫弱可壓流動(dòng)的情況,并用D2V17 模型進(jìn)行了驗(yàn)算。Klass等[92]進(jìn)一步延續(xù)了Lee 等的工作,并將其應(yīng)用到了弱可壓熱流中,使用D2V17 和D2V37 模型進(jìn)行了驗(yàn)證,計(jì)算結(jié)果比Meng 和Zhang 的擴(kuò)散-反彈邊界條件表現(xiàn)更優(yōu)??傮w來(lái)看,針對(duì)多層速度模型的邊界條件研究工作還處于起步階段,且大部分局限于弱可壓熱流和Dirichlet 邊界類型;對(duì)于壓縮性較強(qiáng)的跨聲速、超聲速流動(dòng),以及流動(dòng)出入口涉及的無(wú)反射邊界條件,尚無(wú)相關(guān)工作涉足。將單層速度模型中的各類邊界條件處理方法擴(kuò)展到多層速度模型,并在各類流動(dòng)中廣泛驗(yàn)證測(cè)試,仍是今后需要較多投入的方向。

    3)高性能算法。LBM 最大優(yōu)勢(shì)是計(jì)算速度和高度并行特性,然而它也有一個(gè)眾所周知的缺點(diǎn):由于需要存儲(chǔ)數(shù)量眾多的離散分布函數(shù),對(duì)計(jì)算機(jī)內(nèi)存的使用量比傳統(tǒng)CFD 方法更多,因此有不少工作聚焦于如何改進(jìn)算法,減少內(nèi)存消耗并提高計(jì)算效率[93-95]。另外,正如上一節(jié)所述,經(jīng)典的LBM(相對(duì)于FDLBM 和FVLBM 及其他混合方法而言)一般使用均勻網(wǎng)格,當(dāng)涉及到邊界層計(jì)算時(shí),將LBM 與自適應(yīng)局部加密網(wǎng)格(AMR)結(jié)合[96-98]是更加經(jīng)濟(jì)有效的方案。這些算法與目前流行的GPU 高性能計(jì)算相結(jié)合,能發(fā)揮LBM 的最大計(jì)算潛力[99-101]。然而,目前這些高性能算法都是集中應(yīng)用在單層速度模型中。對(duì)于多層速度模型而言,特別是三維情況,準(zhǔn)確模擬可壓縮熱流的最經(jīng)濟(jì)模型是D3V103,但其分布函數(shù)離散點(diǎn)的數(shù)量已經(jīng)是D3Q27 的4 倍。這一方面使其內(nèi)存占用問(wèn)題更加突出;另一方面,在并行化和局部自適應(yīng)加密網(wǎng)格時(shí),不同塊之間的界面數(shù)據(jù)交互量也更大,需要相當(dāng)謹(jǐn)慎地處理。然而這些往往是工程實(shí)際應(yīng)用中決定算法效率的制約因素。發(fā)展適用于多層速度模型的高性能算法,包括內(nèi)存節(jié)約算法、高效率通信和AMR 方法、CPU+GPU 異構(gòu)高效算法,并在各類復(fù)雜幾何外形、復(fù)雜流動(dòng)中廣泛測(cè)試,應(yīng)引起業(yè)內(nèi)學(xué)者的廣泛關(guān)注。

    總體來(lái)說(shuō),LBM 模型,特別是多層速度模型,仍然需要諸多補(bǔ)充和完善才能成為一種可靠、高效的流體力學(xué)研究手段和應(yīng)用技術(shù)。正如1998 年Chen 和Doolen[2]在文章結(jié)尾中提到,LBM 多相流和復(fù)雜反應(yīng)系統(tǒng)模型主要集中于等溫弱可壓流問(wèn)題,適應(yīng)于可壓縮熱流的LBM 模型仍待開(kāi)發(fā)。20 余年后的今天,這些方面雖然取得了一些進(jìn)步,但是離可靠性要求還有不少差距。我們期望多層速度模型能夠百花齊放,特別是具有較完備數(shù)學(xué)理論框架的Hermite 多項(xiàng)式模型,能在這些方面繼續(xù)發(fā)展壯大并做出應(yīng)有的貢獻(xiàn)。

    附錄A

    猜你喜歡
    平衡態(tài)等溫流動(dòng)
    從平衡態(tài)到非平衡態(tài)
    物理與工程(2024年6期)2024-12-16 00:00:00
    初析固體物理學(xué)中平衡態(tài)的熱力學(xué)條件
    EPDM/PP基TPV非等溫結(jié)晶行為的研究
    流動(dòng)的光
    流動(dòng)的畫
    為什么海水會(huì)流動(dòng)
    快速檢測(cè)豬鏈球菌的環(huán)介導(dǎo)等溫?cái)U(kuò)增方法
    納米CaCO3對(duì)FEP非等溫結(jié)晶動(dòng)力學(xué)的影響
    “三態(tài)”模型:化學(xué)平衡移動(dòng)教學(xué)有效的教學(xué)思維模型
    流動(dòng)的光線
    亚洲精品第二区| 1000部很黄的大片| 色吧在线观看| 日本av手机在线免费观看| 丝袜脚勾引网站| 久久久久久久亚洲中文字幕| 欧美日韩综合久久久久久| 少妇被粗大猛烈的视频| 亚洲色图综合在线观看| 国产精品久久久久久精品古装| av视频免费观看在线观看| 99视频精品全部免费 在线| 观看av在线不卡| 亚洲av电影在线观看一区二区三区| 下体分泌物呈黄色| 国产av一区二区精品久久 | 男女边摸边吃奶| 美女国产视频在线观看| 夜夜骑夜夜射夜夜干| 亚洲美女搞黄在线观看| 成人亚洲欧美一区二区av| 日韩精品有码人妻一区| 精品久久久久久久末码| 国产日韩欧美亚洲二区| 欧美日韩一区二区视频在线观看视频在线| 精品久久久久久电影网| 欧美成人a在线观看| 六月丁香七月| 蜜桃亚洲精品一区二区三区| 伊人久久国产一区二区| 久久久久国产网址| 亚洲av成人精品一二三区| 色视频在线一区二区三区| 免费在线观看成人毛片| 91狼人影院| 在线观看av片永久免费下载| 精品久久国产蜜桃| 国产精品久久久久成人av| av黄色大香蕉| 网址你懂的国产日韩在线| 国产黄色免费在线视频| 丰满乱子伦码专区| 人妻系列 视频| 国产精品人妻久久久久久| 交换朋友夫妻互换小说| 噜噜噜噜噜久久久久久91| 女性被躁到高潮视频| 国产大屁股一区二区在线视频| 三级经典国产精品| 亚洲国产av新网站| 美女cb高潮喷水在线观看| 99热这里只有精品一区| 在线天堂最新版资源| av免费在线看不卡| 高清欧美精品videossex| 久久精品国产亚洲av天美| 五月开心婷婷网| 中文欧美无线码| 国产精品女同一区二区软件| 内地一区二区视频在线| 美女脱内裤让男人舔精品视频| 日韩制服骚丝袜av| 99久久精品国产国产毛片| 国产亚洲av片在线观看秒播厂| 精品国产一区二区三区久久久樱花 | 亚洲激情五月婷婷啪啪| 亚洲综合色惰| 中文字幕精品免费在线观看视频 | 国产成人91sexporn| 亚洲成人手机| 高清在线视频一区二区三区| av在线蜜桃| 又粗又硬又长又爽又黄的视频| 亚洲四区av| 国产成人精品一,二区| 久久久久久久久久久免费av| 99久久中文字幕三级久久日本| 午夜老司机福利剧场| 中文字幕人妻熟人妻熟丝袜美| 老司机影院毛片| 国产成人91sexporn| 美女cb高潮喷水在线观看| 国产日韩欧美在线精品| 三级国产精品欧美在线观看| 日韩强制内射视频| 国产精品伦人一区二区| 国产高清不卡午夜福利| 有码 亚洲区| 亚洲欧美一区二区三区国产| 亚洲无线观看免费| 99久久精品热视频| 久久久久网色| 80岁老熟妇乱子伦牲交| 91aial.com中文字幕在线观看| 综合色丁香网| 欧美极品一区二区三区四区| 亚洲欧美一区二区三区黑人 | 最近2019中文字幕mv第一页| 精品少妇黑人巨大在线播放| 久久久久久久久大av| 一级毛片我不卡| 国产成人a区在线观看| 在线天堂最新版资源| 久久久国产一区二区| 国产人妻一区二区三区在| 赤兔流量卡办理| 日日啪夜夜撸| 最近中文字幕高清免费大全6| 国产爽快片一区二区三区| 国产男女超爽视频在线观看| 有码 亚洲区| 国产乱来视频区| 精品少妇黑人巨大在线播放| 亚州av有码| 女性被躁到高潮视频| 人妻 亚洲 视频| 久久99精品国语久久久| 色视频在线一区二区三区| 99热这里只有精品一区| 久久久久久九九精品二区国产| 精品一区在线观看国产| 久久久久久人妻| 人体艺术视频欧美日本| av一本久久久久| 男人狂女人下面高潮的视频| 欧美xxⅹ黑人| 伦理电影免费视频| 九九爱精品视频在线观看| 欧美三级亚洲精品| 日韩强制内射视频| 精品99又大又爽又粗少妇毛片| 欧美日韩综合久久久久久| 大码成人一级视频| 色视频www国产| www.av在线官网国产| av免费在线看不卡| 美女高潮的动态| 黄色一级大片看看| 国产亚洲午夜精品一区二区久久| 欧美日韩亚洲高清精品| 人妻夜夜爽99麻豆av| 黑人猛操日本美女一级片| 国产69精品久久久久777片| 人妻制服诱惑在线中文字幕| 美女国产视频在线观看| 国产人妻一区二区三区在| 人人妻人人爽人人添夜夜欢视频 | 午夜激情久久久久久久| 一级av片app| 人人妻人人澡人人爽人人夜夜| 亚洲精品乱久久久久久| 免费在线观看成人毛片| 日本免费在线观看一区| 亚洲av电影在线观看一区二区三区| 欧美极品一区二区三区四区| 久久精品人妻少妇| 亚洲成人手机| 精品视频人人做人人爽| 少妇被粗大猛烈的视频| 国产精品一区www在线观看| 99国产精品免费福利视频| 国产精品秋霞免费鲁丝片| 精品一区二区三区视频在线| 不卡视频在线观看欧美| 高清av免费在线| 妹子高潮喷水视频| 日本黄大片高清| 中国国产av一级| 亚洲精品乱码久久久v下载方式| 精品熟女少妇av免费看| 久久久精品94久久精品| 中文字幕人妻熟人妻熟丝袜美| 精品人妻偷拍中文字幕| 大片免费播放器 马上看| 免费人成在线观看视频色| 久久99精品国语久久久| 毛片女人毛片| 亚洲性久久影院| 一区二区av电影网| 久久人人爽人人片av| 亚洲经典国产精华液单| 亚洲精品国产成人久久av| 成年av动漫网址| 国产女主播在线喷水免费视频网站| 丰满迷人的少妇在线观看| 老熟女久久久| 国产爱豆传媒在线观看| 亚洲成人手机| 久久久欧美国产精品| 草草在线视频免费看| a 毛片基地| 免费大片黄手机在线观看| av又黄又爽大尺度在线免费看| 99热这里只有是精品在线观看| 欧美最新免费一区二区三区| 成人漫画全彩无遮挡| 极品少妇高潮喷水抽搐| av播播在线观看一区| 插阴视频在线观看视频| 久久久欧美国产精品| www.色视频.com| 久久精品夜色国产| 国产有黄有色有爽视频| 亚洲欧美日韩卡通动漫| 小蜜桃在线观看免费完整版高清| 亚洲美女视频黄频| 99热这里只有是精品在线观看| 欧美日韩综合久久久久久| 午夜激情久久久久久久| 亚洲不卡免费看| 午夜福利高清视频| 一本久久精品| 亚洲性久久影院| 国产精品国产三级国产av玫瑰| 青春草国产在线视频| 色综合色国产| 国产国拍精品亚洲av在线观看| 嘟嘟电影网在线观看| 日本-黄色视频高清免费观看| 偷拍熟女少妇极品色| 免费黄频网站在线观看国产| 91aial.com中文字幕在线观看| 夜夜骑夜夜射夜夜干| 伦精品一区二区三区| 精品国产三级普通话版| 韩国高清视频一区二区三区| 纵有疾风起免费观看全集完整版| 好男人视频免费观看在线| 久久久精品免费免费高清| 日韩不卡一区二区三区视频在线| 久久久久性生活片| 纵有疾风起免费观看全集完整版| 亚洲婷婷狠狠爱综合网| 日本黄色日本黄色录像| 一级毛片电影观看| 精品久久久噜噜| 精华霜和精华液先用哪个| 国产淫片久久久久久久久| 九草在线视频观看| 偷拍熟女少妇极品色| 国产日韩欧美在线精品| 少妇丰满av| 国产 精品1| 亚洲成色77777| 中国国产av一级| 国产成人精品一,二区| 看非洲黑人一级黄片| 熟女人妻精品中文字幕| 在线天堂最新版资源| 校园人妻丝袜中文字幕| 国产成人午夜福利电影在线观看| 欧美精品国产亚洲| 国产日韩欧美亚洲二区| 亚洲美女视频黄频| 国产v大片淫在线免费观看| 纵有疾风起免费观看全集完整版| 春色校园在线视频观看| 女性生殖器流出的白浆| 伦理电影免费视频| 国产欧美日韩精品一区二区| 777米奇影视久久| 中文字幕亚洲精品专区| 久久99精品国语久久久| 少妇 在线观看| 精品国产一区二区三区久久久樱花 | 麻豆精品久久久久久蜜桃| 亚洲电影在线观看av| 久久99热这里只频精品6学生| 精品一品国产午夜福利视频| 亚州av有码| 欧美高清性xxxxhd video| av在线观看视频网站免费| 我要看黄色一级片免费的| 国产在线免费精品| 视频中文字幕在线观看| 欧美高清成人免费视频www| 爱豆传媒免费全集在线观看| 亚洲精品一区蜜桃| 黄片无遮挡物在线观看| 日韩一区二区三区影片| 久久精品人妻少妇| 国产 精品1| 插逼视频在线观看| 欧美成人一区二区免费高清观看| 精品少妇黑人巨大在线播放| 欧美一区二区亚洲| 日韩在线高清观看一区二区三区| 一级a做视频免费观看| 高清在线视频一区二区三区| 日本wwww免费看| 啦啦啦啦在线视频资源| 只有这里有精品99| 日韩一区二区三区影片| 边亲边吃奶的免费视频| 插阴视频在线观看视频| 亚洲中文av在线| 中国国产av一级| 成人一区二区视频在线观看| 一边亲一边摸免费视频| 中文字幕久久专区| 自拍欧美九色日韩亚洲蝌蚪91 | 国产av码专区亚洲av| 日韩成人av中文字幕在线观看| 小蜜桃在线观看免费完整版高清| 国产国拍精品亚洲av在线观看| videos熟女内射| 欧美日韩精品成人综合77777| 麻豆精品久久久久久蜜桃| 国产爱豆传媒在线观看| 熟女人妻精品中文字幕| 精品人妻偷拍中文字幕| 国产黄色视频一区二区在线观看| 交换朋友夫妻互换小说| 在线免费十八禁| 国产av码专区亚洲av| 久久韩国三级中文字幕| www.av在线官网国产| 在线观看美女被高潮喷水网站| 80岁老熟妇乱子伦牲交| a级毛片免费高清观看在线播放| 日本欧美视频一区| 自拍欧美九色日韩亚洲蝌蚪91 | 大话2 男鬼变身卡| 国产精品熟女久久久久浪| av天堂中文字幕网| 大片免费播放器 马上看| 久久鲁丝午夜福利片| 丰满迷人的少妇在线观看| 午夜视频国产福利| 久久ye,这里只有精品| 大又大粗又爽又黄少妇毛片口| 亚洲精品乱久久久久久| 国产成人精品一,二区| 高清日韩中文字幕在线| av黄色大香蕉| 美女国产视频在线观看| 国产 一区 欧美 日韩| 日本wwww免费看| 99热这里只有精品一区| 久久99热这里只有精品18| 成年av动漫网址| 欧美bdsm另类| 亚洲av不卡在线观看| 深夜a级毛片| 国产亚洲一区二区精品| 少妇人妻一区二区三区视频| 美女视频免费永久观看网站| 久久国产乱子免费精品| 中文字幕制服av| 国产真实伦视频高清在线观看| 国产乱来视频区| 国产黄色免费在线视频| 日本黄色片子视频| 日本色播在线视频| 亚洲国产成人一精品久久久| 亚洲av成人精品一二三区| 国产91av在线免费观看| 精品一区二区免费观看| 老师上课跳d突然被开到最大视频| 日韩不卡一区二区三区视频在线| 亚洲精品一区蜜桃| 中国三级夫妇交换| 免费少妇av软件| 国产亚洲91精品色在线| 精品一品国产午夜福利视频| 国产精品爽爽va在线观看网站| 99久久精品一区二区三区| 2018国产大陆天天弄谢| 成人午夜精彩视频在线观看| 成人影院久久| 尾随美女入室| 亚洲精品久久午夜乱码| 精品酒店卫生间| av在线播放精品| 国产成人freesex在线| 综合色丁香网| 欧美激情极品国产一区二区三区 | 亚洲欧美一区二区三区黑人 | h视频一区二区三区| 国产男人的电影天堂91| 国产伦精品一区二区三区四那| 欧美国产精品一级二级三级 | 国产精品国产三级专区第一集| 欧美日韩亚洲高清精品| 国产精品国产三级专区第一集| 国产精品99久久99久久久不卡 | 人体艺术视频欧美日本| 一本一本综合久久| 欧美三级亚洲精品| a 毛片基地| av网站免费在线观看视频| 一本一本综合久久| 亚洲精品亚洲一区二区| 99热全是精品| 日本色播在线视频| 亚洲av在线观看美女高潮| 亚洲美女视频黄频| 中文在线观看免费www的网站| 亚洲综合色惰| 又黄又爽又刺激的免费视频.| 在线观看三级黄色| 久久毛片免费看一区二区三区| 亚州av有码| 国产色爽女视频免费观看| 国产成人免费观看mmmm| 少妇精品久久久久久久| 夫妻午夜视频| 久热这里只有精品99| 精品一区在线观看国产| 内射极品少妇av片p| 最近中文字幕高清免费大全6| 亚洲国产日韩一区二区| 18禁在线无遮挡免费观看视频| 国产黄片美女视频| 成人无遮挡网站| 一级毛片黄色毛片免费观看视频| 丝袜脚勾引网站| 男女边摸边吃奶| 成年免费大片在线观看| 久久精品国产自在天天线| 久久久久久久久久久丰满| 少妇人妻久久综合中文| 精品一区二区三卡| 91久久精品国产一区二区三区| 另类亚洲欧美激情| 91在线精品国自产拍蜜月| 伦理电影大哥的女人| 六月丁香七月| 欧美xxⅹ黑人| 少妇熟女欧美另类| 深夜a级毛片| 蜜桃久久精品国产亚洲av| 亚洲国产精品专区欧美| 日本欧美视频一区| 激情 狠狠 欧美| 欧美人与善性xxx| 一二三四中文在线观看免费高清| 人体艺术视频欧美日本| av女优亚洲男人天堂| 亚洲一区二区三区欧美精品| 男人狂女人下面高潮的视频| 99久久精品热视频| 黄色欧美视频在线观看| 国产精品爽爽va在线观看网站| 久久精品国产亚洲av涩爱| 一二三四中文在线观看免费高清| 99视频精品全部免费 在线| 亚洲av免费高清在线观看| 高清黄色对白视频在线免费看 | 一级黄片播放器| 日日啪夜夜撸| 超碰97精品在线观看| 亚洲av中文av极速乱| 亚洲精品日本国产第一区| 在线观看av片永久免费下载| 国产精品久久久久久久久免| 99久久精品热视频| 精品一区二区三卡| 久久人人爽av亚洲精品天堂 | 五月天丁香电影| 亚洲精品亚洲一区二区| 欧美成人精品欧美一级黄| www.色视频.com| 夜夜爽夜夜爽视频| 久久综合国产亚洲精品| 五月开心婷婷网| 国产亚洲精品久久久com| 久久久久久九九精品二区国产| 建设人人有责人人尽责人人享有的 | 一级毛片久久久久久久久女| 少妇高潮的动态图| av不卡在线播放| 成人毛片a级毛片在线播放| 亚洲精品aⅴ在线观看| 51国产日韩欧美| 久久99精品国语久久久| 成年av动漫网址| 精品午夜福利在线看| 精品一品国产午夜福利视频| 日韩成人av中文字幕在线观看| 网址你懂的国产日韩在线| 亚洲va在线va天堂va国产| 国产女主播在线喷水免费视频网站| 在线免费观看不下载黄p国产| 一区在线观看完整版| 国产亚洲av片在线观看秒播厂| 黑人猛操日本美女一级片| 亚洲无线观看免费| 欧美日韩国产mv在线观看视频 | 美女视频免费永久观看网站| 日日撸夜夜添| 国产精品人妻久久久久久| av在线播放精品| 欧美三级亚洲精品| 日韩一区二区三区影片| 国产精品av视频在线免费观看| 欧美一级a爱片免费观看看| 国产精品一区二区在线不卡| 免费av不卡在线播放| 国产精品嫩草影院av在线观看| 日本与韩国留学比较| 亚洲图色成人| 多毛熟女@视频| www.av在线官网国产| 2021少妇久久久久久久久久久| 免费高清在线观看视频在线观看| 欧美极品一区二区三区四区| 777米奇影视久久| 国产成人一区二区在线| 美女内射精品一级片tv| 中文字幕av成人在线电影| 如何舔出高潮| 一级片'在线观看视频| 在线播放无遮挡| 天天躁日日操中文字幕| 久久久久久九九精品二区国产| 18禁动态无遮挡网站| 日韩三级伦理在线观看| 久久av网站| 亚洲综合精品二区| 蜜臀久久99精品久久宅男| 高清毛片免费看| 免费久久久久久久精品成人欧美视频 | 免费高清在线观看视频在线观看| 丝瓜视频免费看黄片| 国国产精品蜜臀av免费| 伦理电影大哥的女人| 简卡轻食公司| 少妇丰满av| 在线观看免费视频网站a站| 久久久久久久精品精品| 寂寞人妻少妇视频99o| 天堂俺去俺来也www色官网| 国产精品一及| 日韩在线高清观看一区二区三区| 国产免费一区二区三区四区乱码| 成人毛片60女人毛片免费| 国产无遮挡羞羞视频在线观看| 97超视频在线观看视频| 91精品一卡2卡3卡4卡| 国产免费一区二区三区四区乱码| 亚洲精品视频女| 蜜桃久久精品国产亚洲av| 欧美少妇被猛烈插入视频| 岛国毛片在线播放| 欧美精品亚洲一区二区| 又大又黄又爽视频免费| 秋霞伦理黄片| 美女福利国产在线 | 美女cb高潮喷水在线观看| 免费高清在线观看视频在线观看| 观看免费一级毛片| 美女xxoo啪啪120秒动态图| 免费观看av网站的网址| 麻豆精品久久久久久蜜桃| 日日撸夜夜添| 在线精品无人区一区二区三 | 国产黄频视频在线观看| videos熟女内射| 欧美日韩视频高清一区二区三区二| 国产有黄有色有爽视频| 一级av片app| 精品午夜福利在线看| 国产欧美亚洲国产| 亚洲国产精品成人久久小说| 少妇人妻 视频| 99久久精品一区二区三区| 国产免费一区二区三区四区乱码| 日日啪夜夜撸| 中文字幕人妻熟人妻熟丝袜美| 一级二级三级毛片免费看| 精品国产露脸久久av麻豆| 日本黄色片子视频| 97超视频在线观看视频| 日韩,欧美,国产一区二区三区| 女性被躁到高潮视频| 成人国产麻豆网| 国产在线男女| 午夜福利影视在线免费观看| 自拍偷自拍亚洲精品老妇| av在线老鸭窝| 色视频在线一区二区三区| 身体一侧抽搐| 亚洲精品国产av蜜桃| 一级a做视频免费观看| 久久久久久久久久人人人人人人| 秋霞在线观看毛片| 亚州av有码| 日本黄色日本黄色录像| 国产美女午夜福利| 菩萨蛮人人尽说江南好唐韦庄| 亚洲性久久影院| 免费大片18禁| 99精国产麻豆久久婷婷| 男人狂女人下面高潮的视频| 国产男女内射视频| 亚洲成人av在线免费| 国产精品一二三区在线看| 午夜激情久久久久久久| 国产精品.久久久| 久久精品国产亚洲网站| 国产男人的电影天堂91| 亚洲中文av在线| a级毛色黄片| 久久久久久久久大av| 欧美人与善性xxx| 黄色视频在线播放观看不卡| 亚洲欧洲国产日韩| 人妻一区二区av| 在线 av 中文字幕| 国产精品av视频在线免费观看| 亚洲欧美成人精品一区二区| 大香蕉97超碰在线| 日本午夜av视频| 午夜视频国产福利|