• <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ì)
    午夜福利视频1000在线观看| 啦啦啦观看免费观看视频高清| 久久精品国产亚洲av涩爱 | 亚洲在线自拍视频| 免费av不卡在线播放| 国产高潮美女av| 亚洲av.av天堂| 一个人免费在线观看电影| 日韩欧美精品免费久久| 国产大屁股一区二区在线视频| 婷婷亚洲欧美| 2021天堂中文幕一二区在线观| 国产精品日韩av在线免费观看| 丰满乱子伦码专区| 欧美xxxx黑人xx丫x性爽| 欧美xxxx黑人xx丫x性爽| 国产视频内射| 久久久久久久久久成人| 国产一区二区三区av在线 | 国产精品一及| 欧美激情在线99| 老司机午夜福利在线观看视频| 少妇的逼水好多| 国产熟女欧美一区二区| 高清日韩中文字幕在线| 久久久久久久精品吃奶| 中出人妻视频一区二区| 看黄色毛片网站| 1000部很黄的大片| 91久久精品国产一区二区三区| 高清毛片免费观看视频网站| 乱码一卡2卡4卡精品| .国产精品久久| 少妇丰满av| 嫩草影院新地址| 少妇人妻精品综合一区二区 | 亚洲黑人精品在线| 一本一本综合久久| 国产色婷婷99| 身体一侧抽搐| 亚洲国产欧洲综合997久久,| 日韩 亚洲 欧美在线| 成熟少妇高潮喷水视频| 国产黄片美女视频| 我的女老师完整版在线观看| 桃红色精品国产亚洲av| 免费人成在线观看视频色| 国产精品亚洲一级av第二区| 色噜噜av男人的天堂激情| 偷拍熟女少妇极品色| 亚洲成人久久爱视频| 国产又黄又爽又无遮挡在线| 亚洲专区国产一区二区| 久久99热这里只有精品18| 毛片一级片免费看久久久久 | 神马国产精品三级电影在线观看| 在线观看美女被高潮喷水网站| 日本撒尿小便嘘嘘汇集6| 少妇的逼水好多| 久久久午夜欧美精品| 九九热线精品视视频播放| 国产淫片久久久久久久久| av在线老鸭窝| 12—13女人毛片做爰片一| 99在线人妻在线中文字幕| 一本一本综合久久| 欧美潮喷喷水| 自拍偷自拍亚洲精品老妇| 欧美性感艳星| 国产精品永久免费网站| 日本黄大片高清| 国产高清视频在线观看网站| 国内精品久久久久久久电影| 韩国av一区二区三区四区| 国产三级在线视频| 最近中文字幕高清免费大全6 | 日韩欧美国产一区二区入口| 日本a在线网址| 亚洲精品456在线播放app | 免费观看的影片在线观看| 亚洲人成网站高清观看| 国内精品美女久久久久久| 免费高清视频大片| 日本撒尿小便嘘嘘汇集6| 伦精品一区二区三区| 久久人人精品亚洲av| 美女xxoo啪啪120秒动态图| 欧美丝袜亚洲另类 | 午夜老司机福利剧场| 乱码一卡2卡4卡精品| 日本-黄色视频高清免费观看| 亚洲av免费在线观看| 男人和女人高潮做爰伦理| 啦啦啦观看免费观看视频高清| 一区二区三区高清视频在线| 欧美性感艳星| 女的被弄到高潮叫床怎么办 | 真人一进一出gif抽搐免费| 人妻丰满熟妇av一区二区三区| 12—13女人毛片做爰片一| 亚洲av日韩精品久久久久久密| 51国产日韩欧美| 欧美高清成人免费视频www| 2021天堂中文幕一二区在线观| 国产黄色小视频在线观看| 成年女人毛片免费观看观看9| 黄色一级大片看看| 国产精品人妻久久久影院| 国产欧美日韩一区二区精品| 窝窝影院91人妻| 男人舔奶头视频| 黄色视频,在线免费观看| 99久国产av精品| 国产精品98久久久久久宅男小说| АⅤ资源中文在线天堂| 男人狂女人下面高潮的视频| 久久精品91蜜桃| 免费高清视频大片| 美女 人体艺术 gogo| 中文资源天堂在线| 一进一出好大好爽视频| 亚洲四区av| 亚洲男人的天堂狠狠| 老师上课跳d突然被开到最大视频| 午夜精品一区二区三区免费看| 少妇熟女aⅴ在线视频| 午夜a级毛片| 欧美成人一区二区免费高清观看| 国产精品野战在线观看| 无遮挡黄片免费观看| 免费人成视频x8x8入口观看| 日韩大尺度精品在线看网址| 美女大奶头视频| 男女视频在线观看网站免费| 国产精品电影一区二区三区| 黄片wwwwww| 18禁裸乳无遮挡免费网站照片| 又爽又黄a免费视频| 不卡一级毛片| 欧美一区二区亚洲| 国产 一区 欧美 日韩| 两个人的视频大全免费| 国产高清有码在线观看视频| 一本精品99久久精品77| 成人av一区二区三区在线看| 舔av片在线| 99热这里只有精品一区| 亚洲成人免费电影在线观看| 乱人视频在线观看| 日韩欧美在线二视频| 亚洲一区二区三区色噜噜| 亚洲精品国产成人久久av| 国产单亲对白刺激| 亚洲成人久久性| 免费大片18禁| 亚洲av成人精品一区久久| 老熟妇仑乱视频hdxx| 日韩欧美免费精品| 亚洲av成人精品一区久久| 亚洲精品在线观看二区| 最后的刺客免费高清国语| 国产精品电影一区二区三区| 国产高清激情床上av| 精品一区二区三区视频在线| 少妇高潮的动态图| 久久精品影院6| 欧美丝袜亚洲另类 | 麻豆成人午夜福利视频| 欧美精品啪啪一区二区三区| 美女高潮喷水抽搐中文字幕| 成年女人毛片免费观看观看9| 九九在线视频观看精品| 在线观看美女被高潮喷水网站| 精品久久久久久久久久久久久| 99久久久亚洲精品蜜臀av| 日韩欧美免费精品| 日日夜夜操网爽| 22中文网久久字幕| 一个人观看的视频www高清免费观看| 99热这里只有是精品50| 亚洲欧美日韩无卡精品| 在线免费十八禁| 国产精品亚洲美女久久久| 精品一区二区免费观看| 欧美中文日本在线观看视频| 麻豆成人av在线观看| 国产一区二区激情短视频| 91麻豆av在线| 97碰自拍视频| 天天躁日日操中文字幕| 国产亚洲精品综合一区在线观看| 美女高潮喷水抽搐中文字幕| 成人二区视频| 不卡一级毛片| 嫩草影院入口| 搡老熟女国产l中国老女人| 草草在线视频免费看| 男人和女人高潮做爰伦理| 亚洲无线观看免费| 日本精品一区二区三区蜜桃| 亚洲av不卡在线观看| 亚洲中文字幕一区二区三区有码在线看| 色综合站精品国产| 日韩欧美在线二视频| 日韩欧美 国产精品| 男女之事视频高清在线观看| 国产av一区在线观看免费| 欧美+亚洲+日韩+国产| 免费搜索国产男女视频| 在线观看舔阴道视频| 久久久久久伊人网av| 非洲黑人性xxxx精品又粗又长| 欧美不卡视频在线免费观看| 亚洲精品粉嫩美女一区| 尤物成人国产欧美一区二区三区| 无遮挡黄片免费观看| 亚洲四区av| 日本熟妇午夜| 亚洲狠狠婷婷综合久久图片| 国内精品久久久久精免费| 日本色播在线视频| 无人区码免费观看不卡| 精品一区二区三区人妻视频| 嫩草影院入口| 国产精品乱码一区二三区的特点| 少妇的逼好多水| 如何舔出高潮| 亚洲精品一区av在线观看| 最后的刺客免费高清国语| 尾随美女入室| 精品久久国产蜜桃| 亚洲在线自拍视频| 亚洲av.av天堂| 欧美3d第一页| 此物有八面人人有两片| 亚洲av免费高清在线观看| 少妇熟女aⅴ在线视频| 精品久久国产蜜桃| 淫妇啪啪啪对白视频| 美女免费视频网站| 桃色一区二区三区在线观看| 成人毛片a级毛片在线播放| 日本 欧美在线| 国产亚洲精品av在线| 国产精品一及| 日本在线视频免费播放| 国产 一区精品| 亚洲欧美日韩高清专用| 又紧又爽又黄一区二区| 欧美性猛交╳xxx乱大交人| 男女那种视频在线观看| 草草在线视频免费看| 两人在一起打扑克的视频| 欧美三级亚洲精品| 在线观看午夜福利视频| 欧美人与善性xxx| 亚洲av成人av| 99热这里只有精品一区| 美女免费视频网站| 国产69精品久久久久777片| 欧美日韩综合久久久久久 | a级一级毛片免费在线观看| 国产蜜桃级精品一区二区三区| 成年免费大片在线观看| 欧美+亚洲+日韩+国产| 午夜免费激情av| 老司机深夜福利视频在线观看| 有码 亚洲区| 尤物成人国产欧美一区二区三区| 精品国产三级普通话版| 99久国产av精品| 波多野结衣高清无吗| 免费观看的影片在线观看| 国产精品国产高清国产av| 露出奶头的视频| 又黄又爽又刺激的免费视频.| 97人妻精品一区二区三区麻豆| 人妻少妇偷人精品九色| 国产极品精品免费视频能看的| 亚洲欧美日韩卡通动漫| 一个人看视频在线观看www免费| 国产黄a三级三级三级人| 级片在线观看| 午夜免费成人在线视频| а√天堂www在线а√下载| 免费看日本二区| 欧美日韩国产亚洲二区| 丰满的人妻完整版| 午夜福利成人在线免费观看| 久久国产乱子免费精品| 麻豆av噜噜一区二区三区| 久久久久久九九精品二区国产| 性色avwww在线观看| 久久久久国产精品人妻aⅴ院| 国产一区二区三区在线臀色熟女| 赤兔流量卡办理| 亚洲精品色激情综合| 亚洲图色成人| 久久精品国产亚洲av香蕉五月| 欧美3d第一页| 国产精品久久久久久亚洲av鲁大| 在现免费观看毛片| 亚洲专区国产一区二区| 成人午夜高清在线视频| 久久久久久久久中文| 久久久成人免费电影| 亚洲人与动物交配视频| 午夜福利在线观看免费完整高清在 | 国产91精品成人一区二区三区| 美女高潮的动态| 国产免费一级a男人的天堂| 亚洲aⅴ乱码一区二区在线播放| 中文字幕高清在线视频| 亚洲国产欧洲综合997久久,| 尤物成人国产欧美一区二区三区| 少妇丰满av| 国产成人av教育| 欧美区成人在线视频| 我要搜黄色片| 91在线观看av| 有码 亚洲区| 午夜久久久久精精品| 91av网一区二区| 麻豆成人午夜福利视频| 我的老师免费观看完整版| 夜夜爽天天搞| 免费无遮挡裸体视频| 亚洲狠狠婷婷综合久久图片| 99久久精品一区二区三区| 久久久久久伊人网av| 如何舔出高潮| av国产免费在线观看| 亚洲人成伊人成综合网2020| 久久久久久九九精品二区国产| 午夜老司机福利剧场| 国产精品乱码一区二三区的特点| 18禁在线播放成人免费| 亚洲av中文字字幕乱码综合| 88av欧美| 国产色爽女视频免费观看| 3wmmmm亚洲av在线观看| 又爽又黄无遮挡网站| 又粗又爽又猛毛片免费看| 男插女下体视频免费在线播放| 色综合色国产| 亚洲黑人精品在线| 少妇人妻精品综合一区二区 | 真人做人爱边吃奶动态| 变态另类丝袜制服| 可以在线观看毛片的网站| 波多野结衣高清作品| 九色国产91popny在线| 久久久久久伊人网av| 精品久久久久久久久久久久久| 亚洲美女视频黄频| 日韩欧美国产在线观看| 在线免费观看的www视频| 国产大屁股一区二区在线视频| 亚洲综合色惰| 精品人妻一区二区三区麻豆 | 久久精品综合一区二区三区| 免费人成在线观看视频色| 久久人人爽人人爽人人片va| 美女被艹到高潮喷水动态| 哪里可以看免费的av片| 日日摸夜夜添夜夜添小说| 国产av不卡久久| 99精品在免费线老司机午夜| 女生性感内裤真人,穿戴方法视频| 大又大粗又爽又黄少妇毛片口| 99热只有精品国产| 午夜精品一区二区三区免费看| 亚洲国产欧洲综合997久久,| 国产激情偷乱视频一区二区| 久久久久久大精品| 日日摸夜夜添夜夜添av毛片 | 伊人久久精品亚洲午夜| 丝袜美腿在线中文| 国产精品亚洲一级av第二区| 亚洲av免费高清在线观看| 日本黄大片高清| 最近在线观看免费完整版| 成人av一区二区三区在线看| 中文字幕精品亚洲无线码一区| 久久久午夜欧美精品| 在线看三级毛片| 婷婷精品国产亚洲av在线| 99热6这里只有精品| 精品人妻1区二区| 免费在线观看成人毛片| 91av网一区二区| 老司机福利观看| 成人国产麻豆网| 久久精品国产亚洲av天美| 精品国产三级普通话版| 真实男女啪啪啪动态图| 天堂动漫精品| 人人妻人人看人人澡| 亚洲性夜色夜夜综合| 国产亚洲精品av在线| 日韩欧美免费精品| 免费无遮挡裸体视频| 久久午夜福利片| 女同久久另类99精品国产91| 男人和女人高潮做爰伦理| 一进一出抽搐gif免费好疼| 国产成人av教育| 婷婷精品国产亚洲av在线| 色视频www国产| 国产午夜精品久久久久久一区二区三区 | 啪啪无遮挡十八禁网站| 欧美区成人在线视频| 久久九九热精品免费| 亚洲一区二区三区色噜噜| 精品欧美国产一区二区三| 人妻夜夜爽99麻豆av| 乱系列少妇在线播放| 男女啪啪激烈高潮av片| 99热这里只有是精品在线观看| 日本与韩国留学比较| 老司机午夜福利在线观看视频| 国产一区二区在线观看日韩| 国产综合懂色| 床上黄色一级片| 波多野结衣巨乳人妻| 亚洲最大成人av| 久久久久免费精品人妻一区二区| 制服丝袜大香蕉在线| 真人做人爱边吃奶动态| 露出奶头的视频| 老司机深夜福利视频在线观看| 99国产精品一区二区蜜桃av| 久久香蕉精品热| 久久久国产成人免费| 成人国产一区最新在线观看| 男人狂女人下面高潮的视频| 色综合亚洲欧美另类图片| 国产色爽女视频免费观看| 九色成人免费人妻av| 小蜜桃在线观看免费完整版高清| 最近中文字幕高清免费大全6 | 观看美女的网站| 欧美日韩瑟瑟在线播放| 男插女下体视频免费在线播放| 12—13女人毛片做爰片一| 国内精品久久久久久久电影| 欧美bdsm另类| 最近最新中文字幕大全电影3| 亚洲熟妇熟女久久| 韩国av一区二区三区四区| 亚洲av熟女| 久久香蕉精品热| 人妻少妇偷人精品九色| 久久国产乱子免费精品| 国产大屁股一区二区在线视频| 亚洲熟妇熟女久久| 国产免费男女视频| 亚洲精品久久国产高清桃花| 国产91精品成人一区二区三区| 三级毛片av免费| 九九在线视频观看精品| 国产高清有码在线观看视频| 18+在线观看网站| 床上黄色一级片| 精品久久久久久久久久免费视频| 免费av毛片视频| 九九爱精品视频在线观看| 精华霜和精华液先用哪个| 日本黄色片子视频| av在线观看视频网站免费| 狂野欧美白嫩少妇大欣赏| ponron亚洲| av女优亚洲男人天堂| 午夜福利在线观看免费完整高清在 | 国产av麻豆久久久久久久| 国产欧美日韩精品亚洲av| 国产在线精品亚洲第一网站| 99久久中文字幕三级久久日本| 高清在线国产一区| 久久亚洲精品不卡| 欧洲精品卡2卡3卡4卡5卡区| 国产免费av片在线观看野外av| 亚洲三级黄色毛片| 校园人妻丝袜中文字幕| 在线天堂最新版资源| 狂野欧美白嫩少妇大欣赏| 婷婷色综合大香蕉| 十八禁国产超污无遮挡网站| 久久久精品欧美日韩精品| 色综合色国产| 午夜精品久久久久久毛片777| 无遮挡黄片免费观看| ponron亚洲| 免费电影在线观看免费观看| 一本久久中文字幕| 日韩精品青青久久久久久| 一卡2卡三卡四卡精品乱码亚洲| 高清日韩中文字幕在线| 国产成人a区在线观看| 亚洲性夜色夜夜综合| 免费无遮挡裸体视频| 1024手机看黄色片| 亚洲真实伦在线观看| 午夜福利成人在线免费观看| 亚洲av美国av| 在线播放无遮挡| 精品久久久久久久末码| 97碰自拍视频| 国内久久婷婷六月综合欲色啪| 哪里可以看免费的av片| 搞女人的毛片| 成年人黄色毛片网站| 国产午夜精品久久久久久一区二区三区 | 深爱激情五月婷婷| 久久香蕉精品热| 很黄的视频免费| 久久99热这里只有精品18| 看片在线看免费视频| 欧美在线一区亚洲| 成人av一区二区三区在线看| 色av中文字幕| 精品久久国产蜜桃| 精品福利观看| 欧美色视频一区免费| 亚洲aⅴ乱码一区二区在线播放| 亚洲第一区二区三区不卡| 国产69精品久久久久777片| 日日撸夜夜添| 人人妻人人看人人澡| 丝袜美腿在线中文| 日韩人妻高清精品专区| 国产久久久一区二区三区| 国国产精品蜜臀av免费| 亚洲av二区三区四区| 久99久视频精品免费| 一区福利在线观看| 国产激情偷乱视频一区二区| 亚洲精品影视一区二区三区av| 91午夜精品亚洲一区二区三区 | 变态另类丝袜制服| 白带黄色成豆腐渣| 欧美xxxx性猛交bbbb| 国产伦精品一区二区三区视频9| 欧美潮喷喷水| 性插视频无遮挡在线免费观看| 又粗又爽又猛毛片免费看| 日本在线视频免费播放| 国产一区二区亚洲精品在线观看| 成人国产综合亚洲| 亚洲国产色片| 亚洲 国产 在线| 99热这里只有精品一区| 九九热线精品视视频播放| 国产高清不卡午夜福利| 久久99热6这里只有精品| 在线免费观看的www视频| 熟妇人妻久久中文字幕3abv| 亚洲欧美日韩高清专用| 一a级毛片在线观看| 啦啦啦韩国在线观看视频| 免费看美女性在线毛片视频| 亚洲电影在线观看av| 91精品国产九色| 国产成人福利小说| 欧美最新免费一区二区三区| 亚洲成a人片在线一区二区| 在线免费观看的www视频| 日韩国内少妇激情av| 在线免费十八禁| 男女做爰动态图高潮gif福利片| 国产精品人妻久久久久久| 少妇猛男粗大的猛烈进出视频 | 久久精品91蜜桃| 91精品国产九色| 国产精品三级大全| 小说图片视频综合网站| 女人被狂操c到高潮| 国产精品一区二区三区四区久久| 国产精品一区二区三区四区免费观看 | 在线看三级毛片| 88av欧美| 少妇高潮的动态图| x7x7x7水蜜桃| 亚洲av美国av| 色综合色国产| 日本撒尿小便嘘嘘汇集6| 日韩欧美一区二区三区在线观看| av在线亚洲专区| 麻豆国产av国片精品| 国产 一区精品| 看片在线看免费视频| 99riav亚洲国产免费| 亚洲精品日韩av片在线观看| 日日摸夜夜添夜夜添小说| 亚洲精品久久国产高清桃花| 美女大奶头视频| 看片在线看免费视频| a在线观看视频网站| 亚洲欧美日韩无卡精品| 美女被艹到高潮喷水动态| 日韩av在线大香蕉| 色哟哟哟哟哟哟| 内地一区二区视频在线| 极品教师在线免费播放| 两个人的视频大全免费| 久久久久精品国产欧美久久久| 国产伦人伦偷精品视频| 亚洲自拍偷在线| 可以在线观看毛片的网站| 18禁裸乳无遮挡免费网站照片| av国产免费在线观看| 日韩欧美在线乱码| 久久午夜亚洲精品久久| 久久婷婷人人爽人人干人人爱| 成人永久免费在线观看视频|