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

    氣液兩相分離的免方程多尺度模擬方法

    2015-08-21 07:01:48楊晨何航行
    化工學(xué)報 2015年6期
    關(guān)鍵詞:狀態(tài)方程氣液投影

    楊晨,何航行

    (重慶大學(xué)動力工程學(xué)院低品位能源利用技術(shù)及系統(tǒng)教育部重點實驗室,重慶400044)

    引 言

    在現(xiàn)實問題中,流體或多或少總會含有其他物質(zhì),絕大多數(shù)流動屬于多相流。氣液兩相流動現(xiàn)象廣泛存在于能源、動力、石油、航天等多個技術(shù)領(lǐng)域,如核反應(yīng)堆內(nèi)的流動沸騰現(xiàn)象等。過去幾十年中,各國學(xué)者對于一些典型的兩相流動傳熱過程進(jìn)行了大量的數(shù)值模擬研究。傳統(tǒng)的數(shù)值方法很難模擬相界面動力學(xué)行為,無法刻畫出兩相流動的細(xì) 節(jié)[1-3]。同時,隨著研究的深入,僅著眼于宏觀尺度的模擬研究已經(jīng)無法滿足研究和發(fā)展的需求[4]。目前,基于微細(xì)尺度(微觀尺度、介觀尺度等)的模擬方法如動力學(xué)Monte Carlo 方法(KMC)、分子動力學(xué)方法(MD)和格子Boltzmann 方法(LBM)等是深入研究這些復(fù)雜過程的有效手段。直接運(yùn)用這些方法在微細(xì)尺度上能夠建立反映過程內(nèi)在機(jī)理的微觀模型和介觀模型,但其模擬過程需要完全依賴計算量巨大的微細(xì)仿真器(fine scale simulator,F(xiàn)SS),經(jīng)常面臨計算機(jī)內(nèi)存開銷巨大、超出當(dāng)前 計算能力的問題,以致無法實現(xiàn)對某些復(fù)雜實際 問題的模擬,導(dǎo)致微細(xì)尺度模擬研究面臨極大的 困難[5]。

    近年來,有學(xué)者提出了一種計算機(jī)輔助分析復(fù)雜系統(tǒng)的新型多尺度模擬框架——免方程方法(equation-free method,EFM)[6-7]。其關(guān)鍵思路是繞過從微細(xì)尺度模型中獲得宏觀模型方程這一傳統(tǒng)過程,為復(fù)雜系統(tǒng)多尺度模擬提供了一種新的思路。目前免方程方法的研究在國外發(fā)展迅速,但國內(nèi)只有少數(shù)文獻(xiàn)提及。早在20世紀(jì)末,Theodoropoulos等[6]已涉及到粗?;瘯r間步進(jìn)(coarse time-stepper,CTS)這種跨尺度耦合思想;2001年Gear 等[7-9]開始從事粗?;队胺e分(coarse projective integration,CPI)的研究;之后Kevrekidis 和Gear 等[10-13]運(yùn)用免方程方法進(jìn)行了一系列研究。由此,該多尺度模擬框架逐漸受到復(fù)雜系統(tǒng)模擬研究者的廣泛關(guān)注,其應(yīng)用研究也逐漸增多。

    LBM 方法由于具有微觀粒子特性,可以方便地描述不同相之間的相互作用,清晰地捕捉氣液相界面[14]。因此,本研究充分運(yùn)用這一優(yōu)勢,提出了一種用于模擬氣液兩相分離過程的免方程多尺度方法。利用該方法模擬水在van der Waals(vdW)、Redlich-Kwong(R-K)和Peng-Robinson(P-R)狀態(tài)方程控制下的相變過程,發(fā)現(xiàn)多尺度模擬能夠從微細(xì)底層很好地模擬出氣液兩相分離過程中密度的分布,所得結(jié)果與LBM 單尺度模擬十分吻合;多尺度模擬計算產(chǎn)生的偽速度與LBM 單尺度模擬結(jié)果基本相同,不會導(dǎo)致偽速度增加。這說明該方法在提高了計算效率的同時保證了計算的穩(wěn)定性,驗證了該多尺度模擬方法的準(zhǔn)確性和高效性。通過免方程方法在氣液兩相分離過程中的應(yīng)用研究,可為進(jìn)一步從細(xì)微尺度出發(fā)對流動沸騰中氣泡行為等復(fù)雜多相流過程的多尺度模擬奠定基礎(chǔ)。

    1 免方程多尺度方法

    1.1 基本框架

    免方程方法(EFM)建立了一種基于閉包(closure)的系統(tǒng)辨識,同時構(gòu)建了微觀/隨機(jī)仿真和傳統(tǒng)連續(xù)科學(xué)計算、數(shù)值分析的橋梁。目前免方程方法中常見的應(yīng)用是時間尺度上的多尺度模擬,因此一般以時間尺度上的多尺度模擬為其基本框架。該基本框架由粗粒化時間步進(jìn)(CTS)和相應(yīng)的數(shù)值技術(shù)構(gòu)成,其中CTS 相當(dāng)于整個框架中的“內(nèi)核”。在EFM 中,粗?;瘯r間步進(jìn)包括“提升”(lifting,借助于適當(dāng)?shù)姆植己瘮?shù)從宏觀集中獲得分子狀態(tài))、演化和“約束”(restriction,通過平均、篩選、光滑處理從微觀或分子模型獲得宏觀狀態(tài)集)3 個過程[6]。其中,提升過程是該組合中首先需實現(xiàn)的過程,它是一種跨尺度過程,能夠把宏觀系統(tǒng)量離散轉(zhuǎn)化成與之相匹配的細(xì)小量離散分布,如在氣體動態(tài)模擬中提升過程把某氣體的密度離散成大量氣體粒子的隨機(jī)離散分布。演化過程把提升過程得到的細(xì)小量離散分布投入到細(xì)小尺度模型中進(jìn)行運(yùn)算,經(jīng)過一些短的演化時間后從中得到新的細(xì)小量離散分布。在CTS 中演化不需要尺度跨越,整個過程都在微細(xì)尺度上完成。最后執(zhí)行約束過程,其目的恰好與提升過程相反,即把新得到的細(xì)小量離散分布反饋到宏觀水平。接下來把從連續(xù)數(shù)次CTS中得到的一系列宏觀系統(tǒng)值抽離出來,運(yùn)用合理的數(shù)值技術(shù)如Euler 法、Runge-Kutta 法、Newton- Raphson 法等對它們進(jìn)行相應(yīng)的數(shù)值處理。重復(fù)進(jìn)行上述3 個過程,達(dá)到對系統(tǒng)分析研究的目的。其基本框架如圖1所示。需要強(qiáng)調(diào)的是,免方程方法并不是一種具體的模擬方法,而是一種多尺度模擬框架。

    圖1 基于EFM 的多尺度框架Fig.1 EFM-based multi-scale framework

    1.2 粗?;队胺e分

    粗粒化投影積分(CPI)是免方程方法中一種具有代表性的方法。其基本思想是在復(fù)雜系統(tǒng)的微細(xì)仿真器上運(yùn)用傳統(tǒng)的數(shù)值外推方法,達(dá)到對復(fù)雜系統(tǒng)進(jìn)行分析研究的目的,即將粗粒化時間步進(jìn)獲取的宏觀量實施外推投影以預(yù)測和捕獲系統(tǒng)的特 性[15-17]。

    這里簡要說明CPI 的基本原理(圖2):假設(shè)模擬過程中的某一小段總共有n個連續(xù)演化步,當(dāng)?shù)? 步的初始宏觀物理值提升至細(xì)微仿真器后,必須通過解析相對困難的細(xì)微仿真器依次計算這n個演化步,才能得到第n步的宏觀物理值。這樣雖然能計算出第n步的宏觀值,但其運(yùn)算效率低。而采用CPI 進(jìn)行多尺度計算,則只需通過該細(xì)微仿真器計算出前k個演化步的宏觀物理值,然后對這k個已獲得的宏觀物理值進(jìn)行有效的外推處理,就能快速準(zhǔn)確地預(yù)測出第n步的宏觀值。由于中間省去了h(即n k- )個演化步的計算量以及運(yùn)用了較為準(zhǔn)確的外推處理,CPI 能在保證系統(tǒng)模擬精度的同時有效地降低內(nèi)存消耗,提高CPU 的計算效率。CPI 模擬全過程如圖3所示。

    圖2 某一小段CPI 模擬過程Fig.2 CPI simulation diagram of a short part

    圖3 CPI 模擬全過程Fig.3 Entire CPI simulation process diagram

    由于CPI 無須獲取系統(tǒng)的宏觀模型方程,模擬研究是在現(xiàn)象及過程的細(xì)微仿真器上直接運(yùn)用傳統(tǒng)外推數(shù)值方法快速準(zhǔn)確地預(yù)測系統(tǒng)的宏觀性能,著 眼于過程細(xì)微尺度模型的多尺度模擬方法能更真實地反映實際過程的內(nèi)在特性。本文應(yīng)用伸縮式投影方法(TPM),選用具有二階精度的Adams-Bashforth預(yù)測校正數(shù)值方法進(jìn)行外推處理,其格式可表示為

    式中,F(xiàn)為變量的預(yù)報值,F(xiàn)為變量的校正值,k為演化數(shù),h為投影數(shù),α、β、χ、φ為Adams- Bashforth 方法中的二級系數(shù),以下仿真計算中α取0.4,β取0.6,χ取0.3,φ取0.7。

    2 基本模型

    2.1 格子Boltzmann 模型

    單組分多相LB 方程為[18-19]

    式中,x、t分別表示空間及時間坐標(biāo);i=0,1,2,…,b表示每個格點的不同離散速度方向;ci表示各個離散方向上的速度矢量;Δt表示格子離散的時間步長;τ為量綱1 松弛時間;為平衡態(tài)粒子分布函數(shù)。

    由于D2Q9 模型比D2Q7 模型的數(shù)值穩(wěn)定性好,本研究選取D2Q9 模型,因此Maxwell 平衡態(tài)分布函數(shù)為

    式中,權(quán)重系數(shù)ω0=4/9,ωi=1/9(i=1,2,3,4),1/36(i=5,6,7,8);c=Δx/Δt為粒子遷移速度;equ為平衡態(tài)速度。密度和動量的宏觀參數(shù)分別為

    根據(jù)文獻(xiàn),要使氣液兩相分離,必須考慮相間相互作用力。Shan 和Chen 提出了一種能夠直接刻畫粒子間相互作用的LB 模型,即利用一個偽勢函數(shù)來反映不同相之間的作用力

    式中,?k(x)、分別是x處第k相和x′處第相的有效密度。若只考慮臨近格點間的相互作用力,格林函數(shù)為

    參數(shù)的絕對值決定了相間相互作用的強(qiáng)度,其符號決定二者之間是相互吸引還是相互排斥,Δx為格子長度。在由Shan 和Chen 提出的單組分多相LB模型中,位置x和位置x′的粒子之間的相互作用力為[13]

    式中,c0為常數(shù),對于D2Q9 格子模型c0=6.0。

    Shan-Chen 偽勢模型通過平衡態(tài)速度來體現(xiàn)相間作用力的影響,新的平衡速度為

    2.2 提升算子和約束算子

    結(jié)合LB 模型氣液兩相流動過程的演變規(guī)律,經(jīng)過多次試驗分析得到提升算子的具體形式

    式中,1/48<λ<1/20。

    約束算子則選取為LB 模型密度的宏觀參數(shù)形式

    2.3 狀態(tài)方程

    粒子間相互作用力由式(8)給定,系統(tǒng)的狀態(tài)方程可表示為[20]

    式中,cs為LB 模型的格子聲速,對于D2Q9模型cs=1/。

    本研究選取了以下3 種狀態(tài)方程(equation of state,EOS):van der Waals(vdW)狀態(tài)方程,Redlich- Kwong(R-K)狀態(tài)方程以及Peng-Robinson(P-R)狀態(tài)方程。

    vdW 狀態(tài)方程

    式中,R為通用氣體常數(shù),格子Boltzmann 方法中取值為1.0。壓力對密度求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)并令其等于零,則可得到臨界點的數(shù)值:ρc=1/(3b),Tc=8a/(27bR)。a和b為常量,它們分別取值a=9/49、b=2/21。

    R-K 狀態(tài)方程

    式中,a和b的值由臨界參數(shù)表達(dá):a=

    P-R 狀態(tài)方程

    式中,a和b的值由臨界參數(shù)表達(dá):a=為偏心因子,與工質(zhì)本身有關(guān),水的偏心因子ω=0.344。

    3 仿真結(jié)果與分析

    本研究數(shù)值計算區(qū)域的格子劃分為256×256,平衡態(tài)分布函數(shù)采用D2Q9 模型,松弛因子 1τ=,上下邊界條件和左右邊界條件均為周期性邊界條件,偽勢密度函數(shù)為任意常數(shù)。初始密度值ρ=ρc,為了使方程演化,在初始密度附近取0.01 的擾動。模擬結(jié)果均在以內(nèi)存4GB、CPU 主頻3.1GHz 的四核Xeon 處理器為計算平臺的同一臺計算機(jī)上運(yùn)行得到。

    3.1 不同狀態(tài)方程的結(jié)果對比

    圖4給出了不同狀態(tài)方程控制下LB 模型模擬所得的密度曲線與實驗值的對比,圖中的實驗值根據(jù)水蒸氣性質(zhì)圖表繪制得出,離散點為數(shù)值模擬計算得出的結(jié)果。通過對比可以發(fā)現(xiàn),vdW 狀態(tài)方程的溫度取值范圍最小,T/Tc取值范圍約為[0.86,1.0];P-R 狀態(tài)方程的溫度取值范圍最大,約為[0.76,1.0];R-K 狀態(tài)方程中T/Tc取值范圍約為[0.8,1.0]。在可取的溫度范圍內(nèi),vdW 狀態(tài)方程控制下的模擬結(jié)果偏離實驗值較多,密度相對誤差較大,氣相密度比實驗值大,隨著溫度升高模擬結(jié)果與實驗值的相對誤差逐漸減小,而液相密度模擬結(jié)果比實驗值小,與實驗值的相對誤差基本維持在13%左右。R-K 和P-R 狀態(tài)方程模擬所得曲線與實驗值吻合較好,對于氣相密度R-K狀態(tài)方程得到的結(jié)果與實驗值最接近[圖4(b)],對于液相密度P-R 狀態(tài)方程得到的結(jié)果較好。

    圖5和圖6分別給出了vdW 和R-K、R-K 和P-R 狀態(tài)方程控制下最大偽速度隨密度比的變化曲線。通過對比可以看出,隨著密度比的增大,各模型的最大偽速度均會增大。對于某一給定的密度比,R-K 模型的最大偽速度明顯小于vdW 模型,而P-R 模型的最大偽速度小于R-K 模型;同時,P-R 模型密度比的取值范圍最大,可以處理較高密度比的問題。綜上所述,P-R 狀態(tài)方程能較好地模擬水的相變過程。

    圖4 不同狀態(tài)方程控制下的飽和密度曲線與實驗值的對比Fig.4 Comparison of coexistence curves obtained from simulations with experimental values for different EOS

    圖5 vdW 和R-K 狀態(tài)方程最大偽速度的變化曲線Fig.5 Maximum magnitude of spurious current changes with density ratio for vdW and R-K EOS

    圖6 R-K 和P-R 狀態(tài)方程最大偽速度的變化曲線Fig.6 Maximum magnitude of spurious current changes with density ratio for R-K and P-R EOS

    3.2 LBM 與TPM 模擬的對比

    圖7~圖12給出了LBM 單尺度模擬與TPM 多尺度模擬在不同狀態(tài)方程控制下所得結(jié)果的對比及其與實驗值的對比。總計算步數(shù)為20000,此部分TPM 多尺度模擬選取的演化數(shù)和投影數(shù)均為400。

    圖7 vdW 模型的飽和密度曲線對比Fig.7 Comparison of coexistence curves for vdW EOS

    圖8 vdW 模型最大偽速度對比Fig.8 Comparison of spurious current for vdW EOS

    圖9 R-K 模型的飽和密度曲線對比Fig.9 Comparison of coexistence curves for R-K EOS

    圖10 R-K 模型最大偽速度對比Fig.10 Comparison of spurious current for R-K EOS

    圖11 P-R 模型的飽和密度曲線對比Fig.11 Comparison of coexistence curves for P-R EOS

    圖12 P-R 模型最大偽速度對比Fig.12 Comparison of spurious current for P-R EOS

    由于vdW 狀態(tài)方程本身形式簡單,從圖7可以看出vdW 模型的飽和密度曲線偏離實驗值較大,但 其變化趨勢與實驗值保持一致。TPM 多尺度模擬結(jié)果與LBM 單尺度模擬結(jié)果十分接近,而且TPM 多尺度模擬結(jié)果比LBM 模擬結(jié)果更接近實驗值。

    圖8給出了vdW模型最大偽速度隨密度比的變化情況??梢钥闯?,TPM 模擬產(chǎn)生的偽速度與LBM模擬產(chǎn)生的偽速度基本相同,TPM 模擬產(chǎn)生的偽速度除了最后一個點之外都比LBM 模擬小,說明采用此種情況下的TPM 多尺度模擬不會導(dǎo)致模型偽速度增大,保證了模型計算的穩(wěn)定性。

    圖9給出了R-K 模型經(jīng)過20000 步迭代后的飽和密度曲線??梢钥闯?,TPM 多尺度模擬的結(jié)果與LBM 單尺度模擬的結(jié)果在氣相部分吻合良好,而且與實驗值十分接近;液相部分TPM 多尺度模擬結(jié)果比LBM 單尺度模擬結(jié)果更接近實驗值。模擬結(jié)果比vdW 模型更加準(zhǔn)確,模擬計算出的密度值準(zhǔn)確性較好。在低溫處,液相密度的相對誤差小于2%;氣相密度絕對誤差很小,但相對誤差較大。這是由狀態(tài)方程的性質(zhì)造成的,因為當(dāng)溫度遠(yuǎn)低于臨界溫度時密度比會變得很大。隨著溫度逐漸接近飽和溫度,氣相密度的相對誤差從27%減小到7%。

    圖10給出了R-K 模型最大偽速度隨密度比的變化趨勢??梢钥闯?,TPM 多尺度模擬結(jié)果與LBM單尺度模擬結(jié)果基本保持一致,TPM 產(chǎn)生的偽速度平均值小于LBM 單尺度模擬產(chǎn)生的偽速度平均值,說明采用此種情況下的TPM 多尺度模擬降低了模型的偽速度,提高了模型的穩(wěn)定性。

    圖11給出了P-R 模型經(jīng)過20000 步迭代后的飽和密度曲線。可以看出,TPM 多尺度模擬結(jié)果與LBM 單尺度模擬結(jié)果十分吻合,氣液兩相的密度變化趨勢與實驗值較為接近,但P-R 模型產(chǎn)生的誤差比R-K 模型產(chǎn)生的誤差大。在低溫處,液相密度的相對誤差小于10%;但氣相密度相對誤差還是較大,盡管絕對誤差很小。

    圖12給出了P-R 模型最大偽速度隨密度比的變化趨勢??梢钥闯?,P-R 模型的密度比變化范圍比R-K 模型大,能處理高密度比情況下的問題。TPM 多尺度模擬的偽速度小于LBM 單尺度模擬的偽速度,說明此種情況下采用TPM 多尺度模擬不僅沒有導(dǎo)致模型的偽速度增大,反而降低了模型的偽速度,提高了模型的穩(wěn)定性,充分體現(xiàn)了該多尺度模擬方法的有效性。

    3.3 不同投影步數(shù)的模擬結(jié)果對比

    改變投影和演化的步數(shù),分別選取演化400 步、投影400 步(TPM1),演化200 步、投影600 步(TPM2),演化100 步、投影700 步(TPM3),演化100 步、投影1900 步(TPM4),其他計算條件不變,比較TPM 多尺度模擬在不同投影情況下的模擬結(jié)果。

    圖13~圖18給出了選取不同的投影步長時3種狀態(tài)方程控制下取得的模擬結(jié)果。

    圖13和圖14為vdW 狀態(tài)方程控制下的模擬結(jié)果??梢钥闯?,當(dāng)投影步數(shù)分別取值400、600、700和1900 時,TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線基本重合。投影數(shù)為400、600 和700 時,模型的最大偽速度保持一致;投影數(shù)為1900時,模型的最大偽速度有所增大。

    圖13 不同投影步數(shù)下vdW 模型的飽和密度曲線對比Fig.13 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for vdW EOS

    圖14 vdW 模型最大偽速度的變化曲線對比Fig.14 Comparison of spurious current obtained from simulations using different extrapolation steps for vdW EOS

    圖15 不同投影步數(shù)下R-K 模型的飽和密度曲線對比Fig.15 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for R-K EOS

    圖15和圖16為R-K 狀態(tài)方程控制下的模擬結(jié)果??梢钥闯觯瑤追N情況下的TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線在氣相區(qū)域基本重合。在液相區(qū)域,投影步數(shù)為400 和600 時所得結(jié)果比LBM 單尺度模擬更接近實驗值,TPM 與LBM 模擬的最大偽速度基本相同;投影數(shù)為700 時,模擬結(jié)果在低溫部分與實驗值的偏差增大,模型的最大偽速度稍有增大;投影數(shù)為1900 時,模擬結(jié)果在低溫部分與實驗值最為接近,模型的偽速度與其他幾種情況基本相同。

    圖16 R-K 模型最大偽速度的變化曲線對比Fig.16 Comparison of spurious current obtained from simulations using different extrapolation steps for R-K EOS

    圖17 不同投影步數(shù)下P-R 模型的飽和密度曲線對比Fig.17 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for P-R EOS

    圖18 P-R 模型最大偽速度的變化曲線對比Fig.18 Comparison of spurious current obtained from simulations using different extrapolation steps for P-R EOS

    圖17和圖18為P-R 狀態(tài)方程控制下的模擬結(jié)果??梢钥闯?,TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線在氣相區(qū)域基本重合。在液相區(qū)域,TPM 多尺度模擬比LBM 單尺度模擬所得曲線在低溫部分與實驗值的偏差增大。當(dāng)選取700 和1900 的投影步長時,TPM 模擬的溫度取值范圍減小,而最大偽速度與LBM 模擬結(jié)果基本保持一致,投影步長為1900 的偽速度稍大。

    在模擬精度和運(yùn)算效率方面,以R-K 模型為例(見表1和表2)。由表1可以看出,對于氣相密度,TPM 多尺度模擬的平均誤差和均方根誤差都大于LBM 模擬的誤差。當(dāng)投影步數(shù)為400、600 和700時,氣相密度的平均誤差隨投影步長增大而增加;繼續(xù)增大投影步數(shù)至1900 時,相對誤差的平均值并沒有增大,而是保持為15%。液相密度方面,除投影步數(shù)為400 的情況之外,TPM 模擬的相對誤差均比LBM 單尺度模擬的相對誤差?。黄骄`差最小為4.0%,分別是投影步數(shù)為600 和1900 時取得的。而從液相的均方根誤差可以看出,TPM 多尺度模擬的誤差(10%左右)均小于LBM 單尺度模擬的誤差(12.2%),說明TPM 多尺度模擬在這一部分有效地提高了模型的準(zhǔn)確性。由表2可以看出,TPM多尺度模擬的CPU 耗時明顯小于LBM 單尺度模擬的CPU 耗時,可見TPM 多尺度模擬能夠大大提高模型的運(yùn)算效率;同時,由總體誤差可以看出TPM多尺度模擬的誤差均小于LBM 單尺度模擬的誤差,充分說明了該方法能夠在提高運(yùn)算效率的同時提高模型的精確度,體現(xiàn)了其在復(fù)雜計算中的優(yōu)勢。當(dāng)投影步數(shù)從700(TPM3)增大至1900(TPM4)時,誤差保持在8%左右,體現(xiàn)了該方法在相當(dāng)范圍內(nèi)的穩(wěn)定性。

    綜上所述,在合理的投影范圍內(nèi)運(yùn)用TPM 多尺度模擬方法對氣液兩相分離過程進(jìn)行模擬研究能夠大幅降低對細(xì)微仿真器的依賴性,說明所提多尺度模擬方法對現(xiàn)階段細(xì)微單尺度模擬中出現(xiàn)的巨大內(nèi)存開銷和CPU 耗時問題實現(xiàn)了實質(zhì)性的解決,體現(xiàn)了其有效性和高效性。當(dāng)在細(xì)微尺度上模擬更復(fù)雜的氣液兩相流動過程或更復(fù)雜的多尺度系統(tǒng)時,采用該多尺度方法能有效地解決細(xì)微單尺度模擬難以實現(xiàn)的難題,從而為復(fù)雜兩相流動傳熱問題的模擬研究提供新思路。

    表1 LBM 和TPM 模擬與實驗值的誤差Table 1 Error of simulation results obtained from LBM and TPM model verified by experimental values

    表2 模擬運(yùn)算效率Table 2 Operational efficiency of simulations

    4 結(jié) 論

    本研究提出了一種用于模擬氣液兩相分離的免方程多尺度模擬方法,該方法以氣液兩相的格子Boltzmann(LB)模型作為介觀仿真器,LB 模型進(jìn)行少量演化步后,運(yùn)用二階伸縮式投影方式對介觀仿真演化得到的氣液相密度進(jìn)行有效的外推處理,能夠快速準(zhǔn)確地獲取后續(xù)演化步的信息,從而實現(xiàn)對氣液兩相分離過程的多尺度模擬研究。

    對氣液兩相分離過程中密度分布和模型產(chǎn)生的偽速度的仿真分析表明,本研究所提出的多尺度模擬方法從介觀尺度出發(fā)真實地反映了氣液兩相分離的過程,在確保多尺度模擬的準(zhǔn)確性和真實性的基礎(chǔ)上有效地解決了細(xì)微尺度模擬中出現(xiàn)的CPU耗時過長、內(nèi)存開銷過大的問題。與以往的LBM單尺度模擬相比,本研究所提出的多尺度模擬方法具有明顯的計算優(yōu)勢以及多尺度特點,說明所提方法對于研究更為復(fù)雜的流動傳熱問題非常具有潛力。通過EFM 方法的改進(jìn)以及LBM 模型本身的改進(jìn),可為下一步工作實現(xiàn)流動沸騰中氣泡動力學(xué)行為等更復(fù)雜和更深層次的兩相流動傳熱過程的數(shù)值模擬研究打下堅實的基礎(chǔ)。另外,將該方法的理念運(yùn)用于其他能源動力系統(tǒng)的模擬也是重要的發(fā)展方向,以進(jìn)一步推動多尺度模擬的研究和發(fā)展。

    符 號 說 明

    a,b——狀態(tài)方程中的系數(shù)

    c——粒子遷移速度

    ci——離散方向上的速度矢量

    cs——LB 模型的格子聲速

    F——變量的校正值

    F——變量的預(yù)報值

    F?——粒子間相互作用力

    fi——粒子分布函數(shù)

    ——平衡態(tài)粒子分布函數(shù)

    ——格林函數(shù)

    h——伸縮式投影方法的微觀投影數(shù)

    i——離散速度方向

    k——伸縮式投影方法的微觀演化數(shù)

    p——壓力

    R——通用氣體常數(shù)

    T——溫度

    t——時間坐標(biāo)

    Δt——格子離散的時間步長

    δt——伸縮式投影方法的微觀演化步長

    u——宏觀速度

    ——平衡態(tài)速度

    ——偽勢函數(shù)

    x——空間坐標(biāo)

    Δx——格子離散的空間步長

    α,β,χ,?——Adams-Bashforth 方法中的系數(shù)

    λ——提升算子中的系數(shù)

    ρ——宏觀密度

    τ——量綱1 松弛時間

    φ——有效密度

    ω——偏心因子

    ωi——權(quán)重系數(shù)

    [1]Guo Zhaoli (郭照立),Zheng Chuguang (鄭楚光),Li Qing (李青),Wang Nengchao ( 王能超).Lattice Boltzmann Method for Hydrodynamics (流體動力學(xué)的格子 Boltzmann 方法)[M].Wuhan:Hubei Science & Technology Press,2002:158.

    [2]He Yaling (何雅玲),Wang Yong (王勇),Li Qing (李慶).Lattice Boltzmann Method:Theory and Applications (格子Boltzmann 方法的理論及應(yīng)用) [M].Beijing:Science Press,2008:156.

    [3]Zeng Jianbang (曾建邦),Li Longjian (李隆鍵),Liao Quan (廖全),Chen Qinghua (陳清華),Cui Wenzhi (崔文智),Pan Liangming (潘良明).Application of lattice Boltzmann method to phase transition process [J].Acta Physica Sinica(物理學(xué)報),2010,59 (1):178-185.

    [4]Deen N G,van Sint Annaland M,Kuipers J A M.Multi-scale modeling of dispersed gas-liquid two-phase flow [J].Chemical Engineering Science,2004,59 (8):1853-1861.

    [5]Armaou A,Kevrekidis I G,Theodoropoulos C.Equation-free gaptooth-based controller design for distributed complex/multiscale processes [J].Computers & Chemical Engineering,2005,29 (4):731-740.

    [6]Kevrekidis I G,Gear C W,Hyman J M,Kevrekidis P G,Runborg O,Theodoropoulos C.Equation-free coarse-grained multiscale computation:enabling microscopic simulators to perform system-level analysis [J].Communications in Mathematical Sciences,2003,1 (4):715-762.

    [7]Kevrekidis I G,Gear C W,Hummer G.Equation-free:the computer-aided analysis of complex multiscale systems [J].AⅠChE Journal,2004,50 (7):1346-1355.

    [8]Gear C W.Projective integration methods for distributions [J].NEC Trans.,2001,130:1-9.

    [9]Gear C W,Kevrekidis I G.Projective methods for stiff differential equations:problems with gaps in their eigenvalue spectrum [J].SⅠAM Journal on Scientific Computing,2003,24 (4):1091-1106.

    [10]Cisternas J,Gear C W,Levin S,et al.Equation-free modelling of evolving diseases:coarse-grained computations with individual-based models [J].Proceedings of the Royal Society of London.Series A:Mathematical,Physical and Engineering Sciences,2004,460 (2050):2761-2779.

    [11]Makeev A G,Kevrekidis I G.Equation-free multiscale computations for a lattice-gas model:coarse-grained bifurcation analysis of the NO + CO reaction on Pt (100) [J].Chemical Engineering Science,2004,59 (8):1733-1743.

    [12]Erban R,Kevrekidis I G,Adalsteinsson D,et al.Gene regulatory networks:a coarse-grained,equation-free approach to multiscale computation [J].Journal of Chemical Physics,2006,124 (8):084106.

    [13]Kevrekidis I G,Samaey G.Equation-free multiscale computation:algorithms and applications [J].Annual Review of Physical Chemistry,2009,60:321-344.

    [14]Zeng Jianbang (曾建邦),Li Longjian (李隆鍵),Liao Quan (廖全),Cui Wenzhi (崔文智),Chen Qinghua (陳清華),Pan Liangming (潘良明).Simulation of boiling process with lattice Boltzmann method [J].Journal of Xi’an Jiaotong University(西安交通大學(xué)學(xué)報),2009,43 (7):25-29.

    [15]Maluckov A M,Ishiguro S,?koric M M.A macro projective integration method in 2D microscopic system applied to nonlinear ion acoustic waves in a plasma [J].Communications in Computational Physics,2008,4 (3):556-574.

    [16]Gear C W,Li J,Kevrekidis I G.The gap-tooth method in particle simulations [J].Physics Letters A,2003,316 (3):190-195.

    [17]Yang Chen (楊晨),Peng Wei (彭偉).Telescoping projective methods for the multiscale simulation of reaction-diffusion processes [J].Journal of Engineering Thermophysics(工程熱物理學(xué)報),2011,32 (3):377-381.

    [18]Shan X,Chen H.Lattice Boltzmann model for simulating flows with multiple phases and components [J].Physical Review E,1993,47 (3):1815.

    [19]Shan X,Chen H.Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation [J].Physical Review E,1994,49 (4):2941.

    [20]Yuan P,Schaefer L.Equations of state in a lattice Boltzmann model [J].Physics of Fluids(1994-present),2006,18 (4):42101.

    猜你喜歡
    狀態(tài)方程氣液投影
    LKP狀態(tài)方程在天然氣熱物性參數(shù)計算的應(yīng)用
    煤氣與熱力(2021年6期)2021-07-28 07:21:30
    微重力下兩相控溫型儲液器內(nèi)氣液界面仿真分析
    解變分不等式的一種二次投影算法
    基于最大相關(guān)熵的簇稀疏仿射投影算法
    氣液分離罐液位計接管泄漏分析
    找投影
    找投影
    基于隨機(jī)與區(qū)間分析的狀態(tài)方程不確定性比較
    CO2 驅(qū)低液量高氣液比井下氣錨模擬與優(yōu)化
    用狀態(tài)方程模擬氨基酸水溶液的熱力學(xué)性質(zhì)
    国产精品久久久久久精品电影小说| 菩萨蛮人人尽说江南好唐韦庄| 久久精品国产自在天天线| av片东京热男人的天堂| 一区二区三区乱码不卡18| 晚上一个人看的免费电影| 美女午夜性视频免费| 黑丝袜美女国产一区| 婷婷色综合大香蕉| 国产一区二区 视频在线| 亚洲欧美一区二区三区国产| 国产午夜精品一二区理论片| 99精国产麻豆久久婷婷| 久久久亚洲精品成人影院| 少妇人妻久久综合中文| 免费人妻精品一区二区三区视频| 亚洲精品,欧美精品| 亚洲av.av天堂| 一区二区三区精品91| 视频在线观看一区二区三区| 26uuu在线亚洲综合色| 69精品国产乱码久久久| 日本色播在线视频| 最近最新中文字幕免费大全7| 成人影院久久| 91精品国产国语对白视频| 中文字幕精品免费在线观看视频| 少妇熟女欧美另类| 女人高潮潮喷娇喘18禁视频| 午夜福利视频在线观看免费| 免费观看无遮挡的男女| 日韩av不卡免费在线播放| 美女xxoo啪啪120秒动态图| 中国三级夫妇交换| 一级毛片电影观看| 又粗又硬又长又爽又黄的视频| 色婷婷av一区二区三区视频| 人人妻人人澡人人看| 国产黄色免费在线视频| 香蕉国产在线看| av片东京热男人的天堂| 亚洲少妇的诱惑av| 国产av一区二区精品久久| 国产免费福利视频在线观看| 亚洲精品国产一区二区精华液| 亚洲精品日本国产第一区| 亚洲第一av免费看| 一级,二级,三级黄色视频| 国产精品欧美亚洲77777| 国产 精品1| 秋霞伦理黄片| 电影成人av| 97在线视频观看| 亚洲欧美中文字幕日韩二区| 国产精品国产av在线观看| 日韩av不卡免费在线播放| 亚洲欧洲日产国产| 亚洲成人一二三区av| 亚洲,欧美,日韩| 欧美亚洲 丝袜 人妻 在线| 男女国产视频网站| 久久精品久久久久久久性| 欧美黄色片欧美黄色片| 中文字幕人妻熟女乱码| 午夜老司机福利剧场| 女人高潮潮喷娇喘18禁视频| 亚洲精品aⅴ在线观看| 亚洲av福利一区| 男人爽女人下面视频在线观看| 国产精品久久久av美女十八| 午夜福利在线观看免费完整高清在| 免费不卡的大黄色大毛片视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| 777米奇影视久久| 午夜福利视频精品| 另类精品久久| 久久久欧美国产精品| 精品国产乱码久久久久久男人| 欧美日韩视频精品一区| 国产一区亚洲一区在线观看| 黄片播放在线免费| 最近的中文字幕免费完整| 免费黄频网站在线观看国产| 亚洲精品美女久久av网站| 纵有疾风起免费观看全集完整版| 99香蕉大伊视频| 在线观看国产h片| 日本av手机在线免费观看| 超碰成人久久| 欧美 日韩 精品 国产| 另类亚洲欧美激情| 人体艺术视频欧美日本| 久久99蜜桃精品久久| 国产精品免费大片| 天堂俺去俺来也www色官网| 国产精品二区激情视频| 国产一区二区在线观看av| 国语对白做爰xxxⅹ性视频网站| 国产野战对白在线观看| 老鸭窝网址在线观看| 欧美另类一区| 中文字幕亚洲精品专区| 欧美日韩综合久久久久久| 爱豆传媒免费全集在线观看| 最近的中文字幕免费完整| 欧美老熟妇乱子伦牲交| 国产精品秋霞免费鲁丝片| 亚洲综合精品二区| 丝袜喷水一区| 国产精品 国内视频| 国产精品欧美亚洲77777| 麻豆精品久久久久久蜜桃| 日韩电影二区| 午夜av观看不卡| 国精品久久久久久国模美| 狠狠精品人妻久久久久久综合| 亚洲经典国产精华液单| 精品久久蜜臀av无| 国产精品国产三级国产专区5o| 99国产精品免费福利视频| 少妇被粗大猛烈的视频| 国产精品久久久久成人av| 美女午夜性视频免费| 久久99一区二区三区| 亚洲综合色惰| 少妇人妻精品综合一区二区| 国产精品久久久av美女十八| 制服诱惑二区| 中文字幕精品免费在线观看视频| 尾随美女入室| 亚洲国产av新网站| 欧美精品亚洲一区二区| 国产白丝娇喘喷水9色精品| 亚洲国产av新网站| 国产片特级美女逼逼视频| 少妇的逼水好多| 91aial.com中文字幕在线观看| 日日摸夜夜添夜夜爱| 国产乱人偷精品视频| 日韩熟女老妇一区二区性免费视频| 国产成人91sexporn| 亚洲国产色片| 欧美激情 高清一区二区三区| av福利片在线| 国产 精品1| 伦精品一区二区三区| 国产熟女午夜一区二区三区| 欧美精品av麻豆av| av国产精品久久久久影院| 国产精品熟女久久久久浪| 久久 成人 亚洲| 我要看黄色一级片免费的| 午夜激情av网站| 777久久人妻少妇嫩草av网站| 婷婷成人精品国产| 国产精品麻豆人妻色哟哟久久| 90打野战视频偷拍视频| 99香蕉大伊视频| 我的亚洲天堂| 黄网站色视频无遮挡免费观看| 男人添女人高潮全过程视频| 精品酒店卫生间| 哪个播放器可以免费观看大片| 国产视频首页在线观看| 日本vs欧美在线观看视频| 亚洲精品国产一区二区精华液| 99久国产av精品国产电影| 日日啪夜夜爽| 国产高清不卡午夜福利| 桃花免费在线播放| 少妇猛男粗大的猛烈进出视频| 波多野结衣av一区二区av| 午夜福利网站1000一区二区三区| 校园人妻丝袜中文字幕| 欧美av亚洲av综合av国产av | 日韩av在线免费看完整版不卡| 我的亚洲天堂| 国产淫语在线视频| 日本av免费视频播放| 人妻人人澡人人爽人人| 国产午夜精品一二区理论片| 成人国产av品久久久| 久久久久久人妻| 亚洲欧美精品自产自拍| 黄色毛片三级朝国网站| 1024香蕉在线观看| 大码成人一级视频| 亚洲国产成人一精品久久久| 视频在线观看一区二区三区| 精品国产国语对白av| 日韩欧美一区视频在线观看| 精品久久久久久电影网| 超碰97精品在线观看| 欧美日韩一区二区视频在线观看视频在线| 亚洲 欧美一区二区三区| 久久久a久久爽久久v久久| 汤姆久久久久久久影院中文字幕| 日本欧美国产在线视频| 国产成人精品在线电影| 在线观看美女被高潮喷水网站| 波多野结衣av一区二区av| freevideosex欧美| 性色avwww在线观看| 久久毛片免费看一区二区三区| 最新中文字幕久久久久| 亚洲精品日本国产第一区| 久久久久精品人妻al黑| 男女高潮啪啪啪动态图| 国产在线一区二区三区精| 男女边摸边吃奶| 国产免费视频播放在线视频| 黄片小视频在线播放| 国产成人91sexporn| 亚洲欧美精品综合一区二区三区 | 日本猛色少妇xxxxx猛交久久| 亚洲精品中文字幕在线视频| 亚洲精品自拍成人| 99久国产av精品国产电影| 亚洲在久久综合| 777米奇影视久久| 女的被弄到高潮叫床怎么办| 婷婷色麻豆天堂久久| 久久韩国三级中文字幕| 高清视频免费观看一区二区| 亚洲精品国产一区二区精华液| 妹子高潮喷水视频| a级毛片在线看网站| 亚洲国产欧美网| 亚洲欧美成人精品一区二区| 精品久久蜜臀av无| 欧美 日韩 精品 国产| 天天操日日干夜夜撸| 免费观看无遮挡的男女| 久久av网站| 久久这里只有精品19| 免费在线观看视频国产中文字幕亚洲 | 日产精品乱码卡一卡2卡三| 亚洲国产av影院在线观看| 一级毛片我不卡| 美女中出高潮动态图| 精品国产乱码久久久久久男人| 精品一区二区免费观看| 综合色丁香网| 成人二区视频| 日韩,欧美,国产一区二区三区| 婷婷色综合大香蕉| 久久99热这里只频精品6学生| 18禁观看日本| 免费高清在线观看日韩| 久久久久精品久久久久真实原创| 嫩草影院入口| 丁香六月天网| 久久人妻熟女aⅴ| 精品久久蜜臀av无| 午夜激情av网站| 亚洲精品一二三| 两性夫妻黄色片| av在线播放精品| 午夜福利一区二区在线看| 女人精品久久久久毛片| 精品99又大又爽又粗少妇毛片| 国产黄色免费在线视频| 777久久人妻少妇嫩草av网站| 大陆偷拍与自拍| 色视频在线一区二区三区| 涩涩av久久男人的天堂| 亚洲美女视频黄频| 天天躁夜夜躁狠狠躁躁| 波多野结衣一区麻豆| 国产黄色视频一区二区在线观看| 成年美女黄网站色视频大全免费| 一区二区三区乱码不卡18| 亚洲国产日韩一区二区| 色视频在线一区二区三区| 青春草视频在线免费观看| 日韩制服骚丝袜av| 纯流量卡能插随身wifi吗| 这个男人来自地球电影免费观看 | 天天躁夜夜躁狠狠躁躁| av免费观看日本| 亚洲五月色婷婷综合| 久久精品国产亚洲av高清一级| 黄色毛片三级朝国网站| 中文字幕人妻丝袜一区二区 | 亚洲精品第二区| 免费女性裸体啪啪无遮挡网站| 日本爱情动作片www.在线观看| 两个人免费观看高清视频| 国产精品亚洲av一区麻豆 | 久久 成人 亚洲| 免费高清在线观看视频在线观看| 一区在线观看完整版| 国产精品欧美亚洲77777| 日韩制服丝袜自拍偷拍| 国产精品 国内视频| 亚洲 欧美一区二区三区| 新久久久久国产一级毛片| 人体艺术视频欧美日本| 看非洲黑人一级黄片| 色吧在线观看| 久久久久久久精品精品| 久久精品国产亚洲av天美| av在线观看视频网站免费| 欧美日韩国产mv在线观看视频| 久久亚洲国产成人精品v| av免费在线看不卡| 国产精品久久久久久精品电影小说| 国产精品麻豆人妻色哟哟久久| 永久免费av网站大全| 国产男人的电影天堂91| 少妇熟女欧美另类| 欧美人与善性xxx| 99香蕉大伊视频| 免费日韩欧美在线观看| 欧美在线黄色| 伦理电影大哥的女人| 爱豆传媒免费全集在线观看| 卡戴珊不雅视频在线播放| 国产精品久久久av美女十八| 久久久亚洲精品成人影院| 啦啦啦中文免费视频观看日本| 免费少妇av软件| 亚洲精品,欧美精品| 久久午夜福利片| 日韩精品免费视频一区二区三区| 国产精品一国产av| 美国免费a级毛片| 蜜桃在线观看..| 搡女人真爽免费视频火全软件| 色婷婷久久久亚洲欧美| 亚洲成av片中文字幕在线观看 | 亚洲内射少妇av| 蜜桃国产av成人99| 中文字幕人妻丝袜制服| 国产日韩一区二区三区精品不卡| 乱人伦中国视频| 永久免费av网站大全| 看免费av毛片| 久久久久久久亚洲中文字幕| 少妇被粗大猛烈的视频| 女的被弄到高潮叫床怎么办| 国产极品粉嫩免费观看在线| 韩国精品一区二区三区| 色婷婷av一区二区三区视频| 日本91视频免费播放| 精品久久久久久电影网| 丁香六月天网| 精品久久蜜臀av无| 日韩精品免费视频一区二区三区| 欧美国产精品一级二级三级| 蜜桃在线观看..| www.熟女人妻精品国产| 另类亚洲欧美激情| 精品第一国产精品| 捣出白浆h1v1| 一二三四中文在线观看免费高清| 久久女婷五月综合色啪小说| 午夜老司机福利剧场| 大陆偷拍与自拍| 丰满少妇做爰视频| 午夜福利在线观看免费完整高清在| 春色校园在线视频观看| 中文字幕另类日韩欧美亚洲嫩草| 天天躁日日躁夜夜躁夜夜| 久久久国产精品麻豆| av卡一久久| 亚洲av电影在线进入| 91精品国产国语对白视频| 伦理电影大哥的女人| 少妇猛男粗大的猛烈进出视频| 女人久久www免费人成看片| 少妇猛男粗大的猛烈进出视频| 五月开心婷婷网| 午夜福利视频精品| 叶爱在线成人免费视频播放| 99久国产av精品国产电影| 久久精品国产亚洲av高清一级| 老司机亚洲免费影院| 高清黄色对白视频在线免费看| 久久 成人 亚洲| 18禁国产床啪视频网站| 国产精品亚洲av一区麻豆 | 久久97久久精品| 国产午夜精品一二区理论片| 色婷婷av一区二区三区视频| 侵犯人妻中文字幕一二三四区| 大陆偷拍与自拍| 又粗又硬又长又爽又黄的视频| 亚洲成人一二三区av| 天天躁夜夜躁狠狠久久av| 午夜福利影视在线免费观看| 中文字幕人妻丝袜一区二区 | 久久精品久久久久久噜噜老黄| 中文字幕精品免费在线观看视频| 人妻 亚洲 视频| 精品久久久久久电影网| 99久国产av精品国产电影| 99国产精品免费福利视频| 久久久国产一区二区| 伊人久久大香线蕉亚洲五| av在线观看视频网站免费| 丝袜脚勾引网站| 成人免费观看视频高清| 国产亚洲欧美精品永久| 国产爽快片一区二区三区| 老鸭窝网址在线观看| 精品久久久精品久久久| 亚洲男人天堂网一区| 亚洲av男天堂| 国产日韩欧美视频二区| 亚洲成人av在线免费| 精品一区在线观看国产| tube8黄色片| 亚洲一区二区三区欧美精品| 国产日韩欧美视频二区| 中文字幕人妻熟女乱码| 精品亚洲成a人片在线观看| 一级毛片 在线播放| 国产日韩欧美在线精品| √禁漫天堂资源中文www| 男女无遮挡免费网站观看| 国产精品亚洲av一区麻豆 | 久久久a久久爽久久v久久| 色网站视频免费| 久久99热这里只频精品6学生| 考比视频在线观看| 中文精品一卡2卡3卡4更新| 国产深夜福利视频在线观看| 最新的欧美精品一区二区| 五月天丁香电影| 人人妻人人澡人人爽人人夜夜| 一级毛片电影观看| 乱人伦中国视频| 国产精品国产三级专区第一集| 一区二区三区四区激情视频| 成年动漫av网址| 日韩精品有码人妻一区| 中文乱码字字幕精品一区二区三区| 中文字幕人妻丝袜一区二区 | 久久久久久伊人网av| 精品国产超薄肉色丝袜足j| 成人亚洲精品一区在线观看| 大码成人一级视频| 久久久久久久国产电影| 日产精品乱码卡一卡2卡三| 九色亚洲精品在线播放| 一区二区三区精品91| 国产欧美日韩一区二区三区在线| 最近最新中文字幕免费大全7| 成人黄色视频免费在线看| 婷婷色av中文字幕| 日韩伦理黄色片| 久久精品国产亚洲av天美| av又黄又爽大尺度在线免费看| h视频一区二区三区| 久久人人爽av亚洲精品天堂| 国产成人精品福利久久| 看免费av毛片| 免费女性裸体啪啪无遮挡网站| 亚洲欧美中文字幕日韩二区| 如日韩欧美国产精品一区二区三区| 五月天丁香电影| 中文字幕另类日韩欧美亚洲嫩草| 国产老妇伦熟女老妇高清| 日韩一本色道免费dvd| 欧美精品亚洲一区二区| 亚洲美女黄色视频免费看| 久久精品国产a三级三级三级| 自拍欧美九色日韩亚洲蝌蚪91| 久久毛片免费看一区二区三区| 久久人人爽av亚洲精品天堂| 岛国毛片在线播放| 咕卡用的链子| 亚洲一码二码三码区别大吗| 免费看av在线观看网站| 一级毛片 在线播放| 精品人妻在线不人妻| 久久久久久久精品精品| 亚洲av电影在线进入| 精品国产国语对白av| 丝袜美腿诱惑在线| 老司机亚洲免费影院| 精品人妻偷拍中文字幕| 成人二区视频| 性少妇av在线| 亚洲人成网站在线观看播放| 中文字幕av电影在线播放| 黄频高清免费视频| 黄网站色视频无遮挡免费观看| 亚洲欧美成人综合另类久久久| 国产一区有黄有色的免费视频| 男女国产视频网站| 国产免费现黄频在线看| 人妻少妇偷人精品九色| 波多野结衣一区麻豆| 久久久a久久爽久久v久久| 欧美精品av麻豆av| 国产精品一区二区在线观看99| 色婷婷久久久亚洲欧美| 国产 一区精品| 成人亚洲欧美一区二区av| 精品久久久精品久久久| 18在线观看网站| 熟妇人妻不卡中文字幕| 国产亚洲午夜精品一区二区久久| 岛国毛片在线播放| 久久久久精品性色| 少妇被粗大猛烈的视频| 精品视频人人做人人爽| 亚洲,一卡二卡三卡| 中文字幕另类日韩欧美亚洲嫩草| 青春草国产在线视频| 看非洲黑人一级黄片| 欧美在线黄色| 91在线精品国自产拍蜜月| 午夜福利在线观看免费完整高清在| 黄网站色视频无遮挡免费观看| 极品人妻少妇av视频| 成人午夜精彩视频在线观看| 卡戴珊不雅视频在线播放| 中文天堂在线官网| 天天躁夜夜躁狠狠躁躁| 国产亚洲午夜精品一区二区久久| 女性生殖器流出的白浆| 日韩 亚洲 欧美在线| 久久久久精品久久久久真实原创| 午夜福利乱码中文字幕| 国产av国产精品国产| 国产成人精品无人区| 另类精品久久| av又黄又爽大尺度在线免费看| 国产精品二区激情视频| 热99国产精品久久久久久7| 街头女战士在线观看网站| 日本91视频免费播放| 国产一区二区 视频在线| 在线天堂最新版资源| 国产乱来视频区| 亚洲综合色惰| 亚洲,欧美精品.| av电影中文网址| 国产黄频视频在线观看| av在线app专区| 亚洲欧洲精品一区二区精品久久久 | 天天影视国产精品| 最近中文字幕高清免费大全6| 国产福利在线免费观看视频| 老汉色∧v一级毛片| 一本—道久久a久久精品蜜桃钙片| 1024香蕉在线观看| 亚洲欧美清纯卡通| 男人操女人黄网站| 日韩精品免费视频一区二区三区| 午夜福利视频精品| 亚洲精品日韩在线中文字幕| 国产野战对白在线观看| 久久久久久久久免费视频了| 大话2 男鬼变身卡| 天天躁夜夜躁狠狠躁躁| 91久久精品国产一区二区三区| 久久久久久久久久人人人人人人| 一边摸一边做爽爽视频免费| av免费在线看不卡| 久久精品国产亚洲av涩爱| 久久久久网色| 欧美bdsm另类| 国产午夜精品一二区理论片| 亚洲综合精品二区| 一区在线观看完整版| 一区福利在线观看| 91精品伊人久久大香线蕉| 婷婷色麻豆天堂久久| 日本91视频免费播放| 日本av手机在线免费观看| 永久免费av网站大全| 丝瓜视频免费看黄片| 成人影院久久| 香蕉丝袜av| 免费看av在线观看网站| 高清av免费在线| 国产精品.久久久| 精品少妇一区二区三区视频日本电影 | 午夜日本视频在线| 久久精品国产亚洲av天美| 日本wwww免费看| 国产亚洲av片在线观看秒播厂| 久久av网站| 免费在线观看黄色视频的| 伊人久久大香线蕉亚洲五| 少妇的丰满在线观看| 日本欧美视频一区| 亚洲综合色网址| 日韩电影二区| 一边亲一边摸免费视频| 精品人妻在线不人妻| 久久午夜综合久久蜜桃| 久久鲁丝午夜福利片| 亚洲精品av麻豆狂野| 久久女婷五月综合色啪小说| av电影中文网址| 日本-黄色视频高清免费观看| 99久国产av精品国产电影| 久久99一区二区三区| 777久久人妻少妇嫩草av网站| av在线观看视频网站免费| 欧美成人精品欧美一级黄| 久久精品久久久久久噜噜老黄| 18禁动态无遮挡网站| 国产成人91sexporn| 精品一区二区三卡| 国产一区二区 视频在线| 日韩熟女老妇一区二区性免费视频| 又黄又粗又硬又大视频| 国产在线一区二区三区精| 成人国语在线视频|