• <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ì)
    欧美精品国产亚洲| 乱人伦中国视频| 中文字幕久久专区| 91久久精品电影网| 国产精品国产三级国产av玫瑰| 亚洲色图综合在线观看| a级毛色黄片| 免费少妇av软件| 国产精品一区二区三区四区免费观看| 国产 一区精品| 欧美精品国产亚洲| 国产精品99久久99久久久不卡 | 一区二区三区精品91| 一级毛片 在线播放| 毛片一级片免费看久久久久| 亚洲美女搞黄在线观看| 老司机影院毛片| 伊人久久国产一区二区| 99国产精品免费福利视频| 日韩av在线免费看完整版不卡| 国产男人的电影天堂91| 下体分泌物呈黄色| 97超碰精品成人国产| 伊人久久精品亚洲午夜| 我的老师免费观看完整版| 最新中文字幕久久久久| 中文字幕av电影在线播放| 美女内射精品一级片tv| 国产69精品久久久久777片| 精品99又大又爽又粗少妇毛片| 国产日韩欧美在线精品| 你懂的网址亚洲精品在线观看| 国产一区二区在线观看av| 日日爽夜夜爽网站| 飞空精品影院首页| 精品酒店卫生间| 午夜av观看不卡| av一本久久久久| 欧美日韩av久久| 中文字幕久久专区| 国产黄色免费在线视频| 视频中文字幕在线观看| 高清毛片免费看| 伦理电影大哥的女人| 婷婷色麻豆天堂久久| 国产欧美另类精品又又久久亚洲欧美| 99久久精品国产国产毛片| 新久久久久国产一级毛片| 久久国产精品大桥未久av| 久久精品久久久久久噜噜老黄| 国精品久久久久久国模美| 蜜桃国产av成人99| 色94色欧美一区二区| 视频区图区小说| 亚洲,一卡二卡三卡| 人妻系列 视频| 午夜免费观看性视频| 满18在线观看网站| 男女国产视频网站| 少妇熟女欧美另类| 成年人免费黄色播放视频| 国产精品一国产av| 卡戴珊不雅视频在线播放| 国产精品蜜桃在线观看| www.av在线官网国产| 婷婷色综合大香蕉| 丰满少妇做爰视频| 午夜福利视频在线观看免费| 久久热精品热| av国产精品久久久久影院| 婷婷色av中文字幕| 老司机影院毛片| 桃花免费在线播放| av女优亚洲男人天堂| 男女无遮挡免费网站观看| 两个人的视频大全免费| 国产精品女同一区二区软件| 黑人巨大精品欧美一区二区蜜桃 | 99热国产这里只有精品6| 亚洲欧美日韩卡通动漫| 国产成人精品一,二区| 免费观看在线日韩| 午夜激情久久久久久久| 99热6这里只有精品| 国产午夜精品一二区理论片| 欧美最新免费一区二区三区| 日本-黄色视频高清免费观看| 国产女主播在线喷水免费视频网站| 欧美日韩成人在线一区二区| 日韩亚洲欧美综合| 亚洲综合精品二区| 麻豆乱淫一区二区| 亚洲成人一二三区av| 久久精品久久久久久久性| 蜜桃久久精品国产亚洲av| 亚洲欧美色中文字幕在线| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 看免费成人av毛片| 在线观看一区二区三区激情| 午夜激情久久久久久久| 一级毛片电影观看| 亚洲成人手机| 久久久久精品性色| 新久久久久国产一级毛片| 久久精品国产亚洲av涩爱| 国产在线视频一区二区| 亚洲第一区二区三区不卡| 一区在线观看完整版| 国产男人的电影天堂91| 亚洲四区av| 成人二区视频| 国产日韩欧美在线精品| 成年人免费黄色播放视频| a级毛片在线看网站| 亚洲精品自拍成人| 色婷婷av一区二区三区视频| 免费不卡的大黄色大毛片视频在线观看| 亚洲国产av影院在线观看| av一本久久久久| 日本欧美国产在线视频| 熟女人妻精品中文字幕| 人妻少妇偷人精品九色| videossex国产| 各种免费的搞黄视频| 大码成人一级视频| 日日摸夜夜添夜夜爱| 午夜激情福利司机影院| 老女人水多毛片| 九草在线视频观看| av在线观看视频网站免费| av又黄又爽大尺度在线免费看| 亚洲高清免费不卡视频| 99久久精品一区二区三区| 免费高清在线观看视频在线观看| 精品国产国语对白av| 亚洲欧美色中文字幕在线| 91午夜精品亚洲一区二区三区| 99九九线精品视频在线观看视频| 老女人水多毛片| 久久热精品热| 好男人视频免费观看在线| 欧美一级a爱片免费观看看| 久久免费观看电影| 国产精品久久久久久av不卡| 少妇的逼水好多| 如日韩欧美国产精品一区二区三区 | 一边亲一边摸免费视频| 飞空精品影院首页| 亚洲av欧美aⅴ国产| 欧美性感艳星| 少妇精品久久久久久久| 久久99精品国语久久久| 精品久久久精品久久久| 97超碰精品成人国产| 国产欧美日韩综合在线一区二区| 国产精品一区www在线观看| 麻豆成人av视频| 亚洲伊人久久精品综合| 韩国高清视频一区二区三区| 色网站视频免费| 精品久久久久久久久av| 九九爱精品视频在线观看| av.在线天堂| 九色成人免费人妻av| 亚洲精品日本国产第一区| 久久人人爽av亚洲精品天堂| 欧美精品国产亚洲| 18禁在线无遮挡免费观看视频| 亚洲五月色婷婷综合| av专区在线播放| 看免费成人av毛片| 观看av在线不卡| 2018国产大陆天天弄谢| 伦理电影免费视频| 伊人亚洲综合成人网| 国产又色又爽无遮挡免| 不卡视频在线观看欧美| 久久久久久久久久久久大奶| .国产精品久久| 亚洲精品视频女| 日韩av免费高清视频| 国精品久久久久久国模美| 大话2 男鬼变身卡| 午夜日本视频在线| 国产男人的电影天堂91| 精品少妇内射三级| 久久久久久久大尺度免费视频| 欧美激情 高清一区二区三区| 边亲边吃奶的免费视频| 十分钟在线观看高清视频www| 热re99久久国产66热| 欧美 日韩 精品 国产| 亚洲精品一二三| 22中文网久久字幕| 自拍欧美九色日韩亚洲蝌蚪91| 最后的刺客免费高清国语| av线在线观看网站| .国产精品久久| 亚洲国产av新网站| 亚洲精品久久午夜乱码| 婷婷成人精品国产| 亚洲,欧美,日韩| 蜜臀久久99精品久久宅男| 99久久综合免费| 欧美变态另类bdsm刘玥| 日韩大片免费观看网站| 久久久久久久大尺度免费视频| 欧美xxⅹ黑人| 男男h啪啪无遮挡| 日韩一区二区视频免费看| 国产熟女午夜一区二区三区 | 麻豆成人av视频| 七月丁香在线播放| 婷婷成人精品国产| 婷婷色麻豆天堂久久| 亚洲国产最新在线播放| 人妻制服诱惑在线中文字幕| 高清av免费在线| 国产精品久久久久久精品古装| 午夜福利在线观看免费完整高清在| 综合色丁香网| 国产av一区二区精品久久| 日韩伦理黄色片| 欧美精品亚洲一区二区| 制服丝袜香蕉在线| 91精品国产国语对白视频| 午夜视频国产福利| 国产精品人妻久久久影院| 999精品在线视频| 国产免费一级a男人的天堂| 美女大奶头黄色视频| 青青草视频在线视频观看| 一个人免费看片子| 黄片无遮挡物在线观看| 狂野欧美白嫩少妇大欣赏| av卡一久久| 午夜免费观看性视频| 观看av在线不卡| 日韩欧美精品免费久久| 成人二区视频| 午夜激情av网站| 日韩成人av中文字幕在线观看| 中文乱码字字幕精品一区二区三区| av视频免费观看在线观看| 亚洲精品日韩在线中文字幕| 91国产中文字幕| 国产男人的电影天堂91| 亚洲情色 制服丝袜| 午夜福利影视在线免费观看| 在线观看人妻少妇| 成年美女黄网站色视频大全免费 | 国产 精品1| 亚洲精品久久成人aⅴ小说 | 免费久久久久久久精品成人欧美视频 | 久久女婷五月综合色啪小说| 欧美日韩精品成人综合77777| 黄色毛片三级朝国网站| 伊人久久国产一区二区| av.在线天堂| 欧美国产精品一级二级三级| 国产日韩欧美在线精品| 国产精品秋霞免费鲁丝片| 最近的中文字幕免费完整| www.av在线官网国产| 91久久精品国产一区二区三区| 国产亚洲精品第一综合不卡 | 人人妻人人添人人爽欧美一区卜| 国产不卡av网站在线观看| 九九久久精品国产亚洲av麻豆| 麻豆精品久久久久久蜜桃| 午夜福利在线观看免费完整高清在| 中文字幕av电影在线播放| 热99久久久久精品小说推荐| 各种免费的搞黄视频| 久久久久久久国产电影| 大香蕉久久成人网| 亚洲欧美日韩另类电影网站| 国产免费福利视频在线观看| 18禁在线无遮挡免费观看视频| 亚洲无线观看免费| 国产高清国产精品国产三级| 久久婷婷青草| 在线观看三级黄色| 成年av动漫网址| 亚洲精品,欧美精品| 亚洲国产色片| 人妻夜夜爽99麻豆av| 精品卡一卡二卡四卡免费| 久久久久久久久久成人| 久久午夜综合久久蜜桃| 午夜免费观看性视频| 日韩大片免费观看网站| 涩涩av久久男人的天堂| 草草在线视频免费看| 韩国高清视频一区二区三区| 日产精品乱码卡一卡2卡三| 在现免费观看毛片| kizo精华| 国产一区二区在线观看日韩| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 永久网站在线| 亚洲国产av影院在线观看| 国产精品久久久久成人av| 人妻 亚洲 视频| 一区二区三区精品91| 美女福利国产在线| 青春草视频在线免费观看| 国产黄色免费在线视频| 久久97久久精品| 欧美丝袜亚洲另类| 亚洲精品成人av观看孕妇| 久久精品人人爽人人爽视色| 国产男女内射视频| 免费黄色在线免费观看| 精品酒店卫生间| 久久久久久久久大av| 日韩一区二区三区影片| 成人亚洲精品一区在线观看| 亚洲av男天堂| 一本大道久久a久久精品| 午夜视频国产福利| 精品久久久久久久久亚洲| 最近手机中文字幕大全| 在线观看免费视频网站a站| 成人国产麻豆网| 久久国内精品自在自线图片| 日韩一区二区视频免费看| 男人操女人黄网站| 最近手机中文字幕大全| 国产午夜精品一二区理论片| 久久精品久久久久久噜噜老黄| 99久国产av精品国产电影| 在线播放无遮挡| 搡女人真爽免费视频火全软件| 老司机影院成人| 免费看光身美女| 免费av中文字幕在线| 精品久久久噜噜| 最近最新中文字幕免费大全7| 亚洲精品乱码久久久v下载方式| 欧美日韩精品成人综合77777| 欧美3d第一页| av免费观看日本| 韩国高清视频一区二区三区| 人妻夜夜爽99麻豆av| 久久精品久久精品一区二区三区| 中文字幕人妻熟人妻熟丝袜美| 日本免费在线观看一区| 久久综合国产亚洲精品| av在线app专区| 精品人妻熟女毛片av久久网站| 一级毛片电影观看| 久久久国产欧美日韩av| 制服诱惑二区| 天堂8中文在线网| 亚洲av成人精品一区久久| 亚洲精品乱码久久久久久按摩| 亚洲国产欧美在线一区| 色婷婷久久久亚洲欧美| 日韩视频在线欧美| 欧美日韩一区二区视频在线观看视频在线| 中国三级夫妇交换| 黑人欧美特级aaaaaa片| 亚洲国产成人一精品久久久| 亚洲人成77777在线视频| 少妇人妻久久综合中文| 一级片'在线观看视频| 高清不卡的av网站| 在线观看www视频免费| 免费看光身美女| 国产精品秋霞免费鲁丝片| 久久韩国三级中文字幕| 熟妇人妻不卡中文字幕| 性高湖久久久久久久久免费观看| 国产亚洲午夜精品一区二区久久| 十八禁高潮呻吟视频| 永久免费av网站大全| 亚洲精品av麻豆狂野| 成人综合一区亚洲| 亚洲精品日韩av片在线观看| 午夜激情福利司机影院| 久久国产精品大桥未久av| 久久人人爽av亚洲精品天堂| 亚洲三级黄色毛片| 最新中文字幕久久久久| 国产精品国产三级国产专区5o| av.在线天堂| 最后的刺客免费高清国语| 性高湖久久久久久久久免费观看| 黄色怎么调成土黄色| 97超碰精品成人国产| 桃花免费在线播放| 自线自在国产av| 婷婷色综合大香蕉| 超碰97精品在线观看| 母亲3免费完整高清在线观看 | 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 乱码一卡2卡4卡精品| 欧美bdsm另类| 欧美成人精品欧美一级黄| 日韩av不卡免费在线播放| 欧美日韩av久久| 国产女主播在线喷水免费视频网站| 亚洲欧洲日产国产| 九九爱精品视频在线观看| 亚洲熟女精品中文字幕| 青春草国产在线视频| av网站免费在线观看视频| 我要看黄色一级片免费的| 国产免费现黄频在线看| 一二三四中文在线观看免费高清| 麻豆成人av视频| 成人漫画全彩无遮挡| 久久这里有精品视频免费| 少妇人妻久久综合中文| 中文字幕制服av| 中文天堂在线官网| 热99国产精品久久久久久7| 波野结衣二区三区在线| 久久99热这里只频精品6学生| 精品一区二区三卡| 自拍欧美九色日韩亚洲蝌蚪91| av在线老鸭窝| 日韩精品有码人妻一区| 建设人人有责人人尽责人人享有的| 国产女主播在线喷水免费视频网站| 一本久久精品| 欧美最新免费一区二区三区| 一二三四中文在线观看免费高清| 国产精品99久久久久久久久| 综合色丁香网| 91精品三级在线观看| 制服人妻中文乱码| 国产成人午夜福利电影在线观看| 又粗又硬又长又爽又黄的视频| 欧美一级a爱片免费观看看| 国产精品人妻久久久影院| 久久久国产欧美日韩av| 色婷婷av一区二区三区视频| 欧美成人午夜免费资源| 婷婷色综合大香蕉| 亚洲精品456在线播放app| 欧美日韩成人在线一区二区| 欧美精品亚洲一区二区| 欧美精品人与动牲交sv欧美| 熟女人妻精品中文字幕| 22中文网久久字幕| 黄色视频在线播放观看不卡| 国产女主播在线喷水免费视频网站| 青春草国产在线视频| 国产色爽女视频免费观看| 少妇人妻久久综合中文| 免费日韩欧美在线观看| 成人亚洲精品一区在线观看| 综合色丁香网| 亚洲精华国产精华液的使用体验| 大片电影免费在线观看免费| 亚洲av成人精品一区久久| 久久久久视频综合| 亚洲国产毛片av蜜桃av| 桃花免费在线播放| 自线自在国产av| 亚洲国产av影院在线观看| 伊人久久国产一区二区| 婷婷色麻豆天堂久久| 亚洲精品成人av观看孕妇| 婷婷色麻豆天堂久久| 97精品久久久久久久久久精品| 大话2 男鬼变身卡| 婷婷色综合www| 波野结衣二区三区在线| 国产免费一区二区三区四区乱码| 亚洲人成77777在线视频| 亚洲国产最新在线播放| 久久久久久久久久人人人人人人| 日本爱情动作片www.在线观看| 欧美激情国产日韩精品一区| 成人18禁高潮啪啪吃奶动态图 | 母亲3免费完整高清在线观看 | 国产免费一级a男人的天堂| 国产精品嫩草影院av在线观看| 99精国产麻豆久久婷婷| 成人无遮挡网站| 99精国产麻豆久久婷婷| 国产成人a∨麻豆精品| 欧美激情极品国产一区二区三区 | 一本一本综合久久| 全区人妻精品视频| 插阴视频在线观看视频| 99国产综合亚洲精品| 99久久综合免费| 夫妻性生交免费视频一级片| 一二三四中文在线观看免费高清| 亚洲在久久综合| 青青草视频在线视频观看| 日韩成人伦理影院| 亚洲av电影在线观看一区二区三区| 97超碰精品成人国产| 国产视频首页在线观看| 欧美变态另类bdsm刘玥| 日本猛色少妇xxxxx猛交久久| 久久狼人影院| 久久久国产一区二区| 夜夜爽夜夜爽视频| 久久久国产一区二区| 日韩电影二区| 色哟哟·www| 男女无遮挡免费网站观看| 国产精品免费大片| 亚洲图色成人| 欧美亚洲 丝袜 人妻 在线| 国产片特级美女逼逼视频| 两个人的视频大全免费| 欧美精品亚洲一区二区| 久久久亚洲精品成人影院| 2021少妇久久久久久久久久久| 最近2019中文字幕mv第一页| 国产国语露脸激情在线看| 精品视频人人做人人爽| 啦啦啦视频在线资源免费观看| 成人毛片60女人毛片免费| 亚洲国产欧美日韩在线播放| 欧美精品国产亚洲| 精品人妻在线不人妻| 国产视频首页在线观看| 在线播放无遮挡| 黑人欧美特级aaaaaa片| 麻豆成人av视频| 亚洲,一卡二卡三卡| 国产成人一区二区在线| 国产高清国产精品国产三级| 欧美变态另类bdsm刘玥| 97在线视频观看| 午夜免费男女啪啪视频观看| 久久影院123| 亚洲在久久综合| 精品久久国产蜜桃| 两个人免费观看高清视频| 男女高潮啪啪啪动态图| 少妇被粗大的猛进出69影院 | videosex国产| 毛片一级片免费看久久久久| 夫妻性生交免费视频一级片| 中文字幕精品免费在线观看视频 | 纵有疾风起免费观看全集完整版| 国产综合精华液| 岛国毛片在线播放| 亚洲成人手机| 永久网站在线| 亚洲人与动物交配视频| 在线观看免费日韩欧美大片 | 免费观看a级毛片全部| 欧美日韩亚洲高清精品| 日本vs欧美在线观看视频| 在线观看一区二区三区激情| 久久久国产一区二区| 精品人妻在线不人妻| 国产日韩一区二区三区精品不卡 | 国产欧美亚洲国产| 高清不卡的av网站| 夫妻性生交免费视频一级片| 一级毛片我不卡| 九九在线视频观看精品| 久久国产亚洲av麻豆专区| 91精品三级在线观看| 国产av国产精品国产| 色94色欧美一区二区| 老司机亚洲免费影院| 亚洲精品aⅴ在线观看| 最近2019中文字幕mv第一页| 日产精品乱码卡一卡2卡三| 欧美亚洲日本最大视频资源| 免费人妻精品一区二区三区视频| 成人二区视频| 国产精品秋霞免费鲁丝片| 精品亚洲成国产av| 精品午夜福利在线看| 精品人妻在线不人妻| 国产综合精华液| 久久久国产精品麻豆| 国产精品久久久久成人av| 伊人亚洲综合成人网| 免费看av在线观看网站| 人人妻人人澡人人看| 午夜影院在线不卡| 国产视频首页在线观看| 日本黄色日本黄色录像| 人妻人人澡人人爽人人| 永久免费av网站大全| 丰满饥渴人妻一区二区三| 看非洲黑人一级黄片| 18禁观看日本| 爱豆传媒免费全集在线观看| 考比视频在线观看| 日韩成人av中文字幕在线观看| 成人亚洲精品一区在线观看| 免费日韩欧美在线观看| 国产极品天堂在线| 纯流量卡能插随身wifi吗| 大片电影免费在线观看免费| .国产精品久久| 成人二区视频| 黑人猛操日本美女一级片| 夜夜看夜夜爽夜夜摸| 嘟嘟电影网在线观看| a级毛片免费高清观看在线播放| 全区人妻精品视频| 狠狠精品人妻久久久久久综合| 制服人妻中文乱码| 国产成人精品无人区| 97在线人人人人妻| 黄色怎么调成土黄色| 亚洲精品乱码久久久久久按摩| 国产亚洲av片在线观看秒播厂| 在现免费观看毛片|