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

    顆??s放理論標(biāo)定未篩選煙煤離散元仿真參數(shù)

    2024-09-24 00:00:00梅瀟吳葦榮劉祥偉
    中國(guó)粉體技術(shù) 2024年2期
    關(guān)鍵詞:煙煤

    摘要: 【目的】為了給未篩選煙煤的仿真研究提供參數(shù)依據(jù), 分析未篩選煙煤的離散元仿真參數(shù), 保證仿真與實(shí)際顆粒的幾何、 材料及運(yùn)動(dòng)學(xué)相似, 實(shí)現(xiàn)未篩選煙煤的可靠仿真研究。 【方法】采用實(shí)驗(yàn)和仿真相結(jié)合的方法, 測(cè)得未篩選煙煤的粒徑分布、 密度、 靜摩擦因數(shù)、 堆積密度、 休止角等基本參數(shù); 基于顆??s放理論, 建立不同粒徑范圍內(nèi)典型顆粒的放大模型,通過(guò)Plackett-Burman、最陡爬坡、 Box-Behnken試驗(yàn)對(duì)未篩選煙煤的泊松比、 切變模量、 滾動(dòng)摩擦因數(shù)、 恢復(fù)系數(shù)、 Johnson-Kendall-Roberts(JKR)表面能等仿真參數(shù)進(jìn)行標(biāo)定。 【結(jié)果】以實(shí)驗(yàn)和仿真休止角相對(duì)誤差最小為優(yōu)化目標(biāo), 得到最優(yōu)參數(shù)組合下的未篩選煙煤的仿真休止角為37.59°, 與休止角實(shí)驗(yàn)值37.81°的誤差為0.58%; 仿真堆積密度為717 kg/m3, 與堆積密度實(shí)驗(yàn)值722 kg/m3的誤差為0.69%。 【結(jié)論】煤-煤恢復(fù)系數(shù)與休止角呈負(fù)相關(guān), 煤-煤和煤-鋼滾動(dòng)摩擦因數(shù)與休止角呈正相關(guān),濺射現(xiàn)象會(huì)阻礙顆粒堆積。

    關(guān)鍵詞: 煙煤; 離散元; 休止角; 顆??s放理論; 參數(shù)標(biāo)定

    中圖分類(lèi)號(hào): TB44; TP391.9; O347.7文獻(xiàn)標(biāo)志碼:A

    引用格式:

    梅瀟, 吳葦榮, 劉祥偉. 顆??s放理論標(biāo)定未篩選煙煤離散元仿真參數(shù)[J]. 中國(guó)粉體技術(shù), 2024, 30(2): 67-81.

    MEI X, WU W R, LIU X W. Discrete elemental simulation parameters of unscreened bituminous coal calibrated by particle scaling theory[J]. China Powder Science and Technology, 2024, 30(2): 67-81.

    港口散料輸送多采用螺旋、 抓斗式卸船機(jī)等大型港口裝卸設(shè)備,利用離散元法(discrete element method,DEM)對(duì)不同類(lèi)型散料的運(yùn)輸狀態(tài)進(jìn)行仿真分析是裝卸設(shè)備結(jié)構(gòu)尺寸設(shè)計(jì)的重要步驟。在對(duì)粒徑較小的散料顆粒進(jìn)行仿真分析時(shí),為了滿足計(jì)算機(jī)的運(yùn)算能力要求,通常需要對(duì)顆粒尺寸放大處理,并進(jìn)行參數(shù)標(biāo)定,保證與真實(shí)顆粒有相同的物理特性。已有研究對(duì)包括植物種子、 養(yǎng)料、 土壤、 混凝土、 礦石等不同顆粒模型進(jìn)行仿真參數(shù)標(biāo)定,結(jié)果表明仿真參數(shù)標(biāo)定的應(yīng)用對(duì)象多、 領(lǐng)域廣[1-9]。

    煤炭是港口散料輸送的重點(diǎn)對(duì)象,中國(guó)經(jīng)濟(jì)高速發(fā)展依賴(lài)于煤炭能源的安全穩(wěn)定供應(yīng)[10-12]。仿真分析作為輸送研究的重要方法,仿真結(jié)果的準(zhǔn)確性直接影響研究者對(duì)輸送機(jī)制的判斷,因此,為了使輸送仿真分析結(jié)果可信,應(yīng)對(duì)煤炭的仿真參數(shù)進(jìn)行研究。夏蕊等[13]對(duì)3種不同粒徑的煤炭顆粒進(jìn)行了休止角實(shí)驗(yàn)并基于DEM進(jìn)行仿真分析,探究了落料高度對(duì)休止角的影響,研究發(fā)現(xiàn),落料高度與休止角度呈負(fù)相關(guān)關(guān)系;李鐵軍等[14-15]針對(duì)粒徑為6~8 mm的煤炭顆粒,探究了物料特性參數(shù)在單獨(dú)或互相作用下對(duì)休止角的影響,得到了干燥、含水煤炭的最優(yōu)仿真參數(shù)組合;Mei等[16]通過(guò)在較長(zhǎng)時(shí)間內(nèi)隨機(jī)調(diào)整且不斷試錯(cuò)的方法得到了煤炭的最優(yōu)仿真參數(shù)組合;Zhang等[17]采用正交試驗(yàn)仿真計(jì)算和多因素非線性聯(lián)合反演的方法,以三軸拉伸實(shí)驗(yàn)數(shù)據(jù)為響應(yīng)目標(biāo)對(duì)煤炭進(jìn)行仿真參數(shù)標(biāo)定,得到了煤炭仿真的最優(yōu)參數(shù)組合,為基于DEM的仿真參數(shù)標(biāo)定提

    供了新思路;Li等[18]選取頂煤塌落形狀的幾何參數(shù)作為響應(yīng)目標(biāo),通過(guò)與實(shí)驗(yàn)結(jié)果對(duì)比,對(duì)數(shù)值阻尼進(jìn)行仿真參數(shù)標(biāo)定,反映了頂煤塌落過(guò)程中煤炭顆粒的多尺度破碎和顆粒間的摩擦作用,得到了數(shù)值阻尼的合理仿真參數(shù)取值;Ma等[19]基于響應(yīng)面法研究了長(zhǎng)焰煤的黏結(jié)顆粒模型參數(shù)對(duì)力學(xué)性能的影響,以顆粒破碎特性為響應(yīng)目標(biāo),基于DEM進(jìn)行了仿真參數(shù)標(biāo)定;Xia等[20]以實(shí)驗(yàn)休止角為響應(yīng)目標(biāo),基于DEM對(duì)含

    水量較大的煤炭顆粒進(jìn)行仿真分析,得到了含水煤炭的最優(yōu)仿真參數(shù)組合。由于煤炭種類(lèi)過(guò)多且物料特性差異較大, 因此目前關(guān)于煤炭物料特性參數(shù)的仿真研究仍不充分,尤其缺少對(duì)未篩選煙煤的相關(guān)研究。

    未篩選煙煤常用于電廠發(fā)電或者作為工業(yè)鍋爐的燃料,使用量極大,但粒徑分布范圍較廣且顆粒不規(guī)則度較大,難以獲得準(zhǔn)確的仿真參數(shù)。本文中針對(duì)粒徑小于12 mm的未篩選煙煤,在考慮粒徑分布的情況下,通過(guò)實(shí)驗(yàn)和仿真獲得對(duì)未篩選煙煤顆粒的休止角影響較顯著的參數(shù),建立休止角與影響顯著參數(shù)的回歸方程,并以實(shí)驗(yàn)和仿真休止角相對(duì)誤差最小為優(yōu)化目標(biāo),確定未篩選煙煤顆粒的最優(yōu)參數(shù)組合,為該類(lèi)煤炭的仿真分析提供基礎(chǔ)參數(shù)。

    1 實(shí)驗(yàn)

    1.1材料

    實(shí)驗(yàn)材料為未篩選煙煤(上海松桿煤炭有限公司), 呈黑色不規(guī)則顆粒狀, 微具黏結(jié)性, 未經(jīng)過(guò)篩分。

    1.2方法

    實(shí)驗(yàn)參考國(guó)家標(biāo)準(zhǔn)GB/T 35017—2018[21],測(cè)得未篩選煙煤的粒徑分布、 密度、 靜摩擦因數(shù)、 堆積密度、 休止角等基本參數(shù)。

    1.2.1 粒徑分布

    使用孔徑分別為2、 4、 6、 8、 10 mm的篩子對(duì)未篩選煙煤樣本進(jìn)行篩選, 重復(fù)10次, 不同粒徑顆粒的質(zhì)量分布曲線與質(zhì)量分?jǐn)?shù)如圖1所示。 由圖可以看出, 粒徑為0~2、 gt;2~4、 gt;4~6、 gt;6~8、 gt;8~10、 gt;10 mm的未篩選煙煤的質(zhì)量分?jǐn)?shù)分別為28.4%、 20.2%、 13.9%、 11.1%、 8.9%、 17.5%。

    1.2.2 密度

    由于未篩選煙煤顆粒形狀為不規(guī)則多面體,因此采用排水法測(cè)定未篩選煙煤顆粒的質(zhì)量和體積。在容積為250 mL的量筒中注入體積為160 mL的自來(lái)水后放在電子秤上去皮,加入未篩選煙煤顆粒至水面凹處與量筒容積180 mL刻度線平齊后,讀取電子秤讀數(shù),重復(fù)7次,計(jì)算平均值。實(shí)驗(yàn)樣本的質(zhì)量、 體積與密度如表1所示。由表1中的數(shù)據(jù)計(jì)算得到未篩選煙煤顆粒密度為1 326 kg/m3。

    1.2.3 靜摩擦因數(shù)

    在測(cè)量煤-煤、 煤-鋼靜摩擦因數(shù)時(shí)應(yīng)使被測(cè)未篩選煙煤顆粒呈純滑動(dòng)狀態(tài),使用強(qiáng)力膠和尺寸相同的2塊亞克力板分別制作被測(cè)煤板和材料煤板。在測(cè)量煤-煤靜摩擦因數(shù)時(shí),將材料煤板固定在靜摩擦因數(shù)測(cè)量?jī)x的滑道上,并把被測(cè)煤板放置在材料煤板上,抬起滑道使傾斜角度逐漸增大,直至被測(cè)煤板滑動(dòng)的瞬間停止傾斜,讀取傾斜角度。煤-鋼靜摩擦因數(shù)的實(shí)驗(yàn)過(guò)程相同,區(qū)別在于被測(cè)煤板放置在材料鋼板上滑動(dòng)。重復(fù)20次,計(jì)算置信水平為0.90的置信區(qū)間下的平均值。摩擦因數(shù)實(shí)驗(yàn)裝置如圖2所示。

    靜摩擦因數(shù)計(jì)算公式[21]為

    f=tan θ ,(1)

    式中: f為靜摩擦因數(shù); θ為靜摩擦角,即角度刻度尺讀數(shù)。

    置信水平為0.90的置信區(qū)間下的靜摩擦角測(cè)量值如圖3所示。由圖可知,煤-煤、 煤-鋼靜摩擦角分別為35.65°、 22.08°,計(jì)算得到煤-煤、 煤-鋼靜摩擦因數(shù)分別為0.72、 0.41。

    1.2.4 堆積密度

    在測(cè)量未篩選煙煤的堆積密度時(shí)應(yīng)保證顆粒自由下落,向出料口直徑為25 mm的漏斗中加入足量未篩選煙煤顆粒后,撤除漏斗出料口處的擋板,未篩選煙煤顆粒落至長(zhǎng)度、 寬度、 高度分別為103、 103、 101 mm的長(zhǎng)方體容器中。當(dāng)未篩選煙煤顆粒堆積量超過(guò)容器容積后停止落料,用毛刷輕輕掃平容器表面的未篩選煙煤顆粒,再用電子秤測(cè)量未篩選煙煤顆粒的質(zhì)量,重復(fù)10次,計(jì)算平均值。實(shí)驗(yàn)樣本的質(zhì)量與堆積密度如表2所示。由表2中的數(shù)據(jù)計(jì)算得到未篩選煙煤顆粒的堆積密度為722 kg/m3。

    1.2.5 休止角

    在測(cè)定未篩選煙煤的休止角時(shí)應(yīng)注意漏斗出料口與托盤(pán)距離合適(過(guò)大造成未篩選煙煤顆粒飛濺, 過(guò)小限制未篩選煙煤顆粒堆積高度, 均不利于堆積), 經(jīng)過(guò)多次預(yù)實(shí)驗(yàn), 選取漏斗出料口與托盤(pán)距離為100 mm,休止角實(shí)驗(yàn)裝置如圖4所示。 在實(shí)驗(yàn)時(shí)向漏斗中加入足量未篩選煙煤顆粒, 然后撤除漏斗出料口處的擋板, 讓未篩選煙煤顆粒自由下落堆積, 堆積完成后在正向和側(cè)向拍照記錄堆積情況。 利用MATLAB軟件對(duì)所得圖片依次進(jìn)行灰度、 二值化、 輪廓提取處理, 再測(cè)量左、 右坡面的休止角, 重復(fù)20次, 計(jì)算置信水平為0.90的置信區(qū)間下休止角的平均值。

    置信水平為0.90的置信區(qū)間下的休止角測(cè)量值如圖5所示。由圖可知,未篩選煙煤的休止角為37.81°。

    2 顆??s放理論

    2.1相似理論及量綱分析

    為了滿足計(jì)算機(jī)的計(jì)算能力要求并且節(jié)省仿真時(shí)間,在仿真過(guò)程中通常需要對(duì)顆粒模型進(jìn)行放大處理。根據(jù)黃松元[22]和Feng等[23]提出的仿真相似理論,經(jīng)過(guò)放大處理后的仿真顆粒模型應(yīng)滿足顆粒應(yīng)變能和運(yùn)動(dòng)方程,即

    E(s, r)=∫s0F(s, r)ds ,(2)

    ma+F(s, r)=F(t) ,(3)

    式中: E(s, r)為顆粒的應(yīng)變能, s為顆粒重疊、 分離、 滑動(dòng)的位移; r為顆粒的半徑; F(s, r)為顆粒之間接觸力的合力; m為顆粒的質(zhì)量; a為顆粒的加速度; F(t)為顆粒所受的外力合力,t為外力作用的時(shí)間。

    根據(jù)相似理論,仿真顆粒密度與真實(shí)顆粒密度應(yīng)成比例;但由于在該條件下難以保持顆粒接觸力相似,因此在實(shí)際應(yīng)用時(shí)選取密度ρ的縮放比λρ為1,選擇長(zhǎng)度L的縮放比λL和時(shí)間t的縮放比λt為h。獲得以下數(shù)量關(guān)系[24]

    λm=λρλ3L=h3gt;0

    λv=λLλ-1t=1

    λF=λρλ4Lλ-2t=h2gt;0

    λE=λρλ2Lλ-2t=1

    λσ=λρλ2Lλ-2t=1

    λε=λρλ2Lλ-2tλ-1ρλ-2Lλ2t=1 ,(4)

    式中: λm為質(zhì)量m的數(shù)量關(guān)系; λv為速度v的數(shù)量關(guān)系; λF為力F的數(shù)量關(guān)系; λE為切變模量E的數(shù)量關(guān)系; λσ為應(yīng)力σ的數(shù)量關(guān)系; λε為應(yīng)變?chǔ)诺臄?shù)量關(guān)系。

    2.2接觸理論

    由于未篩選煙煤微具黏結(jié)性,因此選用Johnson-Kendall-Roberts(JKR)接觸模型,JKR接觸模型可以有效體現(xiàn)顆粒間的黏結(jié)情況。

    2個(gè)顆粒在外載荷FN和表面黏結(jié)共同作用下的等效載荷FD和接觸面半徑r1方程[25-26]為

    FD=FN+3πreEγ+(3πreEγ)2+6πreEγFN ,(5)

    r31=3re4EeFN+3πreEγ+(3πreEγ)2+6πreEγFN ,(6)

    式中:" re為有效顆粒半徑; Eγ為黏結(jié)能; Ee為有效切變模量。

    在顆粒接觸變形的重疊位移增量為Δs時(shí),顆粒的載荷法向力增量ΔFN方程[25-26]為

    ΔFN=2r1EeΔs3FN-3Fc

    3FN-Fc ,(7)

    滿足

    FNFc-r31r3c2=4r1rc3 ,(8)

    其中

    Fc=3πreEγ2 ,(9)

    r3c=3reFc4Ee=9πre2Eγ8Ee ,(10)

    式中: Fc為使2個(gè)顆粒分開(kāi)的臨界拉力; rc為相應(yīng)接觸面的半徑。

    3 仿真參數(shù)標(biāo)定

    3.1顆粒模型

    由于顆粒形狀尺寸直接影響仿真結(jié)果,因此在建立放大處理后的仿真顆粒模型時(shí),應(yīng)以各粒徑范圍內(nèi)的典型顆粒形狀為依據(jù),以實(shí)驗(yàn)堆積密度為評(píng)判標(biāo)準(zhǔn),對(duì)指定放大倍數(shù)時(shí)的仿真顆粒模型形狀進(jìn)行調(diào)整。在仿真堆積密度與實(shí)驗(yàn)堆積密度結(jié)果誤差在合理范圍內(nèi)時(shí),即可認(rèn)為所構(gòu)建的仿真顆粒模型具有代表性。本文中將各粒徑范圍內(nèi)的典型顆粒粒徑放大5倍,并對(duì)顆粒形狀進(jìn)行調(diào)整得到未篩選煙煤的仿真顆粒模型,未篩選煙煤的實(shí)際顆粒和仿真模型如圖6所示。

    3.2仿真參數(shù)

    實(shí)驗(yàn)僅獲得了未篩選煙煤的基本參數(shù),應(yīng)使用仿真對(duì)泊松比、 切變模量、 滾動(dòng)摩擦因數(shù)、 恢復(fù)系數(shù)、 JKR表面能等仿真參數(shù)進(jìn)行標(biāo)定。首先通過(guò)查找相關(guān)文獻(xiàn)[14,27]和離散元仿真軟件EDEM自帶的物料數(shù)據(jù)庫(kù)GEMM確定取值范圍,然后通過(guò)休止角仿真標(biāo)定確定具體取值,未篩選煙煤顆粒的仿真參數(shù)如表3所示。

    3.3仿真休止角的測(cè)量

    根據(jù)相似理論,在EDEM中建立放大漏斗模型,放大比例與未篩選煙煤顆粒的放大比例一致,在漏斗模型中按比例共生成質(zhì)量為150 kg的未篩選煙煤顆粒,待顆粒完全生成后靜置時(shí)間為1 s,隨后撤除漏斗出料口處擋板,使顆粒自由下落,在堆積的顆粒不再發(fā)生運(yùn)動(dòng)時(shí)即可認(rèn)為堆積完成,休止角實(shí)驗(yàn)的仿真過(guò)程如圖7所示。

    由于顆粒的形狀不規(guī)則,堆積坡面凹凸不平,因此為了更準(zhǔn)確地測(cè)量仿真休止角,采用基于Python開(kāi)發(fā)的EDEMpy庫(kù)進(jìn)行仿真結(jié)果分析[28]。

    3.4參數(shù)標(biāo)定試驗(yàn)設(shè)計(jì)

    3.4.1 Plackett-Burman(P-B)試驗(yàn)

    為了確定標(biāo)定參數(shù)對(duì)休止角的影響顯著性,使用Design-Expert軟件對(duì)表3中的范圍參數(shù)進(jìn)行P-B參數(shù)顯著性分析試驗(yàn)。將未篩選煙煤的泊松比(X0),切變模量(X1),煤-煤(X2)、 煤-鋼(X3)恢復(fù)系數(shù),煤-煤(X4)、 煤-鋼滾動(dòng)摩擦因數(shù)(X5),JKR表面能(X6)等7個(gè)試驗(yàn)變量和4個(gè)虛擬變量分別分為低、 高水平,低、 高水平分別為表3中仿真參數(shù)取值范圍的下、 上限值。P-B試驗(yàn)設(shè)計(jì)如表4所示。

    試驗(yàn)變量顯著性分析如表5所示。由表可知,決定系數(shù)R2為0.997 4,非常接近1,表示擬合情況很好;離散系數(shù)為3.04%,表示模型離散概率較低;測(cè)量信噪比Rsn為35.611 7,表示該模型非常適用。

    P-B試驗(yàn)分析圖如圖8所示。由圖8(a)可知,正態(tài)殘差圖的殘差點(diǎn)大致分布在一條直線附近,說(shuō)明殘差符合正態(tài)分布;由圖8(b)可知,預(yù)測(cè)及實(shí)際對(duì)比圖的試驗(yàn)點(diǎn)基本沿著與X軸角度為45°線分布,且離直線較近,表示預(yù)測(cè)值與實(shí)際觀測(cè)值之間存在較好的線性關(guān)系。以上現(xiàn)象均說(shuō)明模型的擬合程度較高。

    3.4.2 最陡爬坡試驗(yàn)

    為了使仿真結(jié)果快速向?qū)嶒?yàn)結(jié)果逼近,對(duì)P-B試驗(yàn)中影響顯著的3個(gè)試驗(yàn)變量煤-煤恢復(fù)系數(shù)、 煤-煤滾動(dòng)摩擦因數(shù)、 煤-鋼滾動(dòng)摩擦因數(shù)進(jìn)行最陡爬坡試驗(yàn)。選擇P-B試驗(yàn)第12組次(最接近真實(shí)休止角)的相關(guān)變量取值,即未篩選煙煤的泊松比為0.4, 未篩選煙煤的切變模量為0.67 GPa, 煤-鋼恢復(fù)系數(shù)為0.15, JKR表面能為14 J/m2。影響顯著變量取值見(jiàn)表6。由表可知,試驗(yàn)編號(hào)3的相對(duì)誤差最小,結(jié)果趨于真實(shí)值。

    3.4.3 Box-Behnken(B-B)試驗(yàn)

    為了確定休止角和影響顯著性試驗(yàn)變量的函數(shù)關(guān)系,以最陡爬坡試驗(yàn)中試驗(yàn)編號(hào)3的變量取值為中間水平,使用Design-Expert軟件對(duì)變量進(jìn)行B-B試驗(yàn),B-B試驗(yàn)水平和試驗(yàn)設(shè)計(jì)如表7所示。

    B-B試驗(yàn)得到二次回歸方程模型的方差分析見(jiàn)表8。由表可知,決定系數(shù)R2為0.987 6,預(yù)測(cè)系數(shù)R2p與校正系數(shù)R2a取值接近(差值小于0.2),說(shuō)明擬合情況很好;離散系數(shù)Cv為0.945%,模型離散概率較低;測(cè)量信噪比Rsn為24.553 7,模型的F值較大且p值較小,說(shuō)明該模型準(zhǔn)確度較高。

    為了進(jìn)一步提高模型的擬合度,舍去表8中的影響不顯著的試驗(yàn)變量X2X5、 X4X5、 X25以?xún)?yōu)化模型,再對(duì)優(yōu)化后的模型進(jìn)行方差分析,結(jié)果見(jiàn)表9。由表可知,優(yōu)化模型的決定系數(shù)R2幾乎不變,但是預(yù)測(cè)系數(shù)R2p、 校正系數(shù)R2a、 測(cè)量信噪比Rsn、 模型的F值均有顯著的增大,離散系數(shù)Cv和模型的p值略有減小,說(shuō)明優(yōu)化后的模型準(zhǔn)確度更高。

    表9對(duì)應(yīng)的優(yōu)化后的二次回歸方程為

    β=-62.21X2+482.41X4+26.88X5-233.75X2X4+84.94X22-1 351.44X24+21.49

    。(11)

    優(yōu)化后的B-B試驗(yàn)分析圖如圖9所示。 由圖9(a)、 (b)可知, 殘差點(diǎn)和試驗(yàn)點(diǎn)均沿直線分布。 由圖9(c)、 (d)可知, 試驗(yàn)點(diǎn)在區(qū)間內(nèi)均勻分布, 并沒(méi)有表現(xiàn)出明顯的規(guī)律或趨勢(shì), 說(shuō)明模型擬合程度較高。

    以變量的高、 低水平為邊界條件,實(shí)驗(yàn)休止角為優(yōu)化目標(biāo),對(duì)式(11)進(jìn)行求解,得到影響顯著的試驗(yàn)變量的最優(yōu)水平組合是:X2為0.500, X4為0.107, X5為0.095。

    4 結(jié)果分析

    4.1顯著試驗(yàn)變量分析

    式(11)中變量X2、 X4、 X5對(duì)休止角的影響曲線如圖10所示。由圖可知,變量X2與休止角呈負(fù)相關(guān),變量X4、 X5與休止角呈正相關(guān)。由圖10(a)可知,圖中曲線尾端有略微上升的趨勢(shì),這是由仿真誤差導(dǎo)致的,處于合理范圍之內(nèi),不影響對(duì)整體趨勢(shì)的判斷。

    式(11)中的交互項(xiàng)X2X4響應(yīng)曲面如圖11所示。由圖11(b)可知,X2對(duì)休止角的影響隨著X4取值增大而增大。原因是漏斗出料口與堆積平面的距離是固定的,X4越小,未篩選煙煤的休止角也越小,即堆積坡峰高度越低,此時(shí)下落顆粒與堆積顆粒碰撞瞬間的動(dòng)能與顆粒下落距離(出料口至堆積坡峰的距離)呈正相關(guān),無(wú)論X2取范圍內(nèi)的任何值,顆粒在碰撞后仍有較大的動(dòng)能,造成嚴(yán)重的濺射現(xiàn)象,不利于顆粒的堆積;X4對(duì)休止角的影響隨著X2取值增大而減小,原因是X2越大,濺射現(xiàn)象越嚴(yán)重,導(dǎo)致顆粒越難堆積,無(wú)論滾動(dòng)摩擦因數(shù)X4取值如何,都難以形成休止角較大的坡面。可見(jiàn),濺射現(xiàn)象是評(píng)判顆粒休止角大小的重要指標(biāo),二者呈負(fù)相關(guān)關(guān)系,而這種現(xiàn)象隨著X4增大或X2減小而緩解。

    4.2仿真參數(shù)組合的驗(yàn)證

    輸入最優(yōu)參數(shù)組合,得到仿真堆積密度為717 kg/m3,與堆積密度實(shí)驗(yàn)值的誤差為0.69%,顆粒模型設(shè)置合理;仿真休止角度為37.59°,與實(shí)驗(yàn)值的誤差為0.58%。仿真休止角與實(shí)驗(yàn)休止角對(duì)比如圖12所示。由圖12可知,由于仿真與實(shí)驗(yàn)休止角的輪廓高度重合,因此認(rèn)為仿真參數(shù)設(shè)置合理,仿真結(jié)果與實(shí)際情況相符。

    5 結(jié)論

    1)通過(guò)實(shí)驗(yàn)得到未篩選煙煤的基本參數(shù),根據(jù)顆??s放理論將仿真顆粒模型的粒徑放大5倍,并通過(guò)P-B試驗(yàn)確定出對(duì)休止角影響顯著的試驗(yàn)變量的排序依次為: 煤-煤滾動(dòng)摩擦因數(shù)、 煤-煤恢復(fù)系數(shù)、 煤-鋼滾動(dòng)摩擦因數(shù),其中煤-煤滾動(dòng)摩擦因數(shù)的影響尤為明顯。

    2)通過(guò)最陡爬坡試驗(yàn)使仿真結(jié)果快速向?qū)嶒?yàn)結(jié)果逼近,通過(guò)B-B試驗(yàn)建立并優(yōu)化了影響顯著的試驗(yàn)變量的二次回歸方程模型,發(fā)現(xiàn)3個(gè)影響顯著的試驗(yàn)變量的一次項(xiàng)、 煤-煤恢復(fù)系數(shù)和煤-煤滾動(dòng)摩擦因數(shù)的二次項(xiàng)以及煤-煤恢復(fù)系數(shù)與滾動(dòng)摩擦因數(shù)的交互項(xiàng)均對(duì)休止角有顯著影響,其中煤-煤恢復(fù)系數(shù)和煤-煤滾動(dòng)摩擦因數(shù)的一次項(xiàng)的影響尤為明顯。同時(shí)通過(guò)分析交互項(xiàng)響應(yīng)曲面的變化趨勢(shì)確定濺射現(xiàn)象會(huì)阻礙未篩選煙煤顆粒的堆積。

    3)確定密度為1 326 kg/m3的未篩選煙煤顆粒的最優(yōu)仿真參數(shù)組合為泊松比為0.4,切變模量為0.67 GPa,煤-煤、 煤-鋼靜摩擦因數(shù)分別為0.72、 0.41,煤-煤、 煤-鋼滾動(dòng)摩擦因數(shù)分別為0.107、 0.095,煤-煤、 煤-鋼恢復(fù)系數(shù)分別為0.5、 0.15,JKR表面能為14 J/m2。該參數(shù)組合下的仿真堆積密度與實(shí)驗(yàn)堆積密度、 仿真休止角度與實(shí)驗(yàn)休止角度的誤差均在合理范圍內(nèi)。

    4)以上仿真參數(shù)組合僅適用于顆粒粒徑放大5倍后的未篩選煙煤的仿真研究。 由于已知煤-煤恢復(fù)系數(shù)和煤-煤滾動(dòng)摩擦因數(shù)的一次項(xiàng)對(duì)休止角的影響尤為明顯, 因此在確定顆粒粒徑其余倍數(shù)下的未篩選煙煤的仿真參數(shù)組合時(shí), 僅需采用同樣的方法重新標(biāo)定煤-煤恢復(fù)系數(shù)和煤-煤滾動(dòng)摩擦因數(shù)即可。

    利益沖突聲明(Conflict of Interests)

    所有作者聲明不存在利益沖突。

    All authors disclose no relevant conflict of interests.

    作者貢獻(xiàn)(Author’s Contributions)

    梅瀟參與了實(shí)驗(yàn)設(shè)計(jì),吳葦榮參與了論文的寫(xiě)作和仿真,劉祥偉參與了論文的修改。所有作者均閱讀并同意了最終稿件的提交。

    The study was designed by MEI Xiao, the manuscript and discrete element simulation completed by WU Weirong, the manuscript modification completed by LIU Xiangwei. All authors have read the last version of paper and consented for submission.

    參考文獻(xiàn)(References)

    [1]黃文景. 砂質(zhì)土壤顆粒離散元模型參數(shù)標(biāo)定[J]. 福建建材, 2023(2): 1-5.

    HUANG W J. Parameter calibration of discrete element model of sandy soil particles[J]. Fujian Building Materials, 2023(2): 1-5.

    [2]洪波, 樊志鵬, 烏蘭圖雅, 等. 揉碎玉米秸稈螺旋輸送仿真離散元模型參數(shù)標(biāo)定[J]. 中國(guó)農(nóng)業(yè)科技導(dǎo)報(bào), 2023, 25(3): 96-106.

    WANG H B, FAN Z P, WU L, et al. Parameter calibration of discrete element model for simulation of crushed corn stalk screw conveying[J]. Journal of Agricultural Science and Technology, 2023, 25(3): 96-106.

    [3]MA G G, SUN Z J, MA H, et al. Calibration of contact parameters for moist bulk of shotcrete based on EDEM[J]. Advances in Materials Science and Engineering, 2022, 6072303.

    [4]FANG M, YU Z H, ZHANG W J, et al. Friction coefficient calibration of corn stalk particle mixtures using Plackett-Burman design and response surface methodology[J]. Powder Technology, 2022, 396: 731-742.

    [5]WANG Y K, WANG G Q, WU S W, et al. A calibration method for ore bonded particle model based on deep learning neural network[J]. Powder Technology, 2023, 420: 118417.

    [6]LI X Y, DU Y F, LIU L, et al. Parameter calibration of corncob based on DEM[J]. Advanced Powder Technology, 2022, 33(8): 103699.

    [7]SONG X F, DAI F, ZHANG F W, et al. Calibration of DEM models for fertilizer particles based on numerical simulations and granular experiments[J]. Computers and Electronics in Agriculture, 2023, 204: 107507.

    [8]ZHU J Z, ZOU M, LIU Y S, et al. Measurement and calibration of DEM parameters of lunar soil simulant[J]. Acta Astronautica, 2022, 191: 169-177.

    [9]DING X T, WANG B B, HE Z, et al. Fast and precise DEM parameter calibration for Cucurbita ficifolia seeds[J]. Biosystems Engineering, 2023, 236: 258-276.

    [10]孟麗英. 新形勢(shì)下港口運(yùn)輸經(jīng)濟(jì)問(wèn)題探析[J]. 商展經(jīng)濟(jì), 2022(5): 89-92.

    MENG L Y. Analysis of the economic problems of port transportation under the new situation[J]. Trade Fair Economy, 2022(5): 89-92.

    [11]賈大山, 徐迪, 蔡鵬. 2021年沿海港口發(fā)展回顧與2022年展望[J]. 中國(guó)港口, 2022(1): 3-14.

    JIA D S, XU D, CAI P. Review of coastal port development in 2021 and prospect in 2022[J]. China Ports, 2022(1): 3-14.

    [12]高宏杰. 煤炭行業(yè)發(fā)展現(xiàn)狀和供需形勢(shì)分析[J]. 中國(guó)煤炭工業(yè), 2022(3):75-77.

    GAO H J. Analysis of development status and supply and demand situation of coal industry[J]. China Coal Industry, 2022(3): 75-77.

    [13]夏蕊, 楊兆建, 李博, 等. 基于離散元法的煤散料堆積角試驗(yàn)研究[J]. 中國(guó)粉體技術(shù), 2018, 24(6): 36-42.

    XIA R, YANG Z J, LI B, et al. Experimental study for repose angle of coal bulk material based on discrete element method[J]. China Powder Science and Technology, 2018, 24(6): 36-42.

    [14]李鐵軍. 煤顆粒離散元模型宏細(xì)觀參數(shù)標(biāo)定及其關(guān)系[D]. 太原:太原理工大學(xué), 2019.

    LI T J. Calibration of DEM model parameters for coal particles and research on relationships between macro and micro parameters[D]. Taiyuan: Taiyuan University of Technology, 2019.

    [15]李鐵軍, 王學(xué)文, 李博, 等. 基于離散元法的煤顆粒模型參數(shù)優(yōu)化[J]. 中國(guó)粉體技術(shù), 2018, 24(5): 6-12.

    LI T J, WANG X W, LI B, et al. Optimization method for coal particle model parameters based on discrete element method[J]. China Powder Science and Technology, 2018, 24(5): 6-12.

    [16]MEI L, HU J Q, YANG J G, et al. Research on parameters of EDEM simulations based on the angle of repose experiment[C]//International Conference on Computer Supported Cooperative Work in Design. Proceedings of the 2016 IEEE 20th International Conference on Computer Supported Cooperative Work in Design. United States: Institute of Electrical and Electronics Engineers Inc., 2016: 570-574.

    [17]ZHANG H, ZHENG M X, RONG Y, et al. Calibration and sensitivity analysis of macro and meso parameters of discrete element model for coal measure soil[C]//International Conference on Applied Mechanics, Materials Physics, and Engineering Structures. Journal of Physics: Conference Series. Institute of Physics, 2023, 2519(1): 012036.

    [18]LI H B, LI D Y, ZHANG W Y, et al. Numerical damping calibration study of particle element method-based dynamic relaxation approach for modeling longwall top-coal caving[J]. Energies, 2021, 14(9): 2348.

    [19]MA H Z, WANG X W, LI B, et al. Calibration of discrete element microparameters of coal based on the response surface method[J]. Particulate Science and Technology, 2022, 40(5): 543-557.

    [20]XIA R, LI B, WANG X W, et al. Measurement and calibration of the discrete element parameters of wet bulk coal[J]. Measurement, 2019, 142: 84-95.

    [21]中國(guó)機(jī)械工業(yè)聯(lián)合會(huì). 連續(xù)搬運(yùn)設(shè)備散狀物料分類(lèi)、 符號(hào)、 性能及測(cè)試方法: GB/T 35017—2018[S]. 北京: 中國(guó)標(biāo)準(zhǔn)出版社, 2018.

    China Machinery Industry Federation. Continuous handling equipments—Classification, symbols, performance and test methods for bulk materials: GB/T 35017—2018[S]. Beijing: Standards Press of China, 2018.

    [22]黃松元. 散體力學(xué)[M]. 北京:機(jī)械工業(yè)出版社, 1993: 159-161.

    HUANG S Y. Mechanics of granular media[M]. Beijing: China Machine Press, 1993:159-161.

    [23]FENG Y T, HAN K, OWEN D R J, et al. On upscaling of discrete element models: similarity principles[J]. Engineering Computations, 2009, 26(6): 599-609.

    [24]任建莉, 周龍海, 韓龍, 等. 基于顆??s放理論的垂直螺旋輸送離散模擬[J]. 過(guò)程工程學(xué)報(bào), 2017, 17(5): 936-943.

    REN J L, ZHOU L H, HAN L, et al. Discrete simulation of vertical screw conveyor based on particle scaling theory[J]. The Chinese Journal of Process Engineering, 2017, 17(5): 936-943.

    [25]E SILVA B B, DA CUNHA E R, DE CARVALHO R M, et al. Modeling and simulation of green iron ore pellet classification in a single deck roller screen using the discrete element method[J]. Powder Technology, 2018, 332: 359-370.

    [26]孫其誠(chéng), 王光謙. 顆粒物質(zhì)力學(xué)導(dǎo)論[M]. 北京: 科學(xué)出版社, 2009: 20-22.

    SUN Q C, WANG G Q. Introduction to granular material mechanics[M]. Beijing: Science Press, 2009: 20-22.

    [27]謝仁海, 渠天祥, 錢(qián)光謨. 構(gòu)造地質(zhì)學(xué)[M]. 2版. 徐州: 中國(guó)礦業(yè)大學(xué)出版社, 2007: 47-49.

    XIE R H, QU T X, QIAN G M. Geotectonics[M]. 2nd ed. Xuzhou: China University of Mining and Technology Press, 2007: 47-49.

    [28]鄒洋, 湯佟, 高自成, 等. 基于顆??s放理論的生石灰粉離散元參數(shù)標(biāo)定[J]. 中國(guó)粉體技術(shù), 2023, 29(2): 81-91.

    ZOU Y, TANG T, GAO Z C, et al. Discrete element parameter calibration of quicklime powder based on particle scaling theory[J]. China Powder Science and Technology, 2023, 29(2): 81-91.

    Discrete elemental simulation parameters of unscreened bituminous coal

    calibrated by particle scaling theory

    MEI Xiao, WU Weirong, LIU Xiangwei

    (Logistics Engineering College, Shanghai Maritime University, Shanghai 201306, China)

    Abstract

    Objective In recent years, the proportion of coal in port traffic has been increased. Unscreened bituminous coal, as an important fuel for power plants or industrial boilers, is one of the widely utilized coal types. Therefore, conducting simulation studies on the transport process of unscreened bituminous coal holds extreme importance. However, according to different simulation conditions, it is usually necessary to scale the particle size of unscreened bituminous coal to meet various simulation requirements. The validity of setting particle discrete element simulation parameters directly affects the accuracy of the simulation results. To obtain discrete elemental simulation parameters for unscreened bituminous coal particles after scaling processing, the simulation parameters should be calibrated to ensure that it meets the geometric similarity, material similarity and kinematic similarity with the actual particles.

    Methods In this study, optimal parameter combinations for the discrete element simulation of unscreened bituminous coal were determined through a combination of physical experiments and simulation. Firstly, the basic parameters of bituminous coal, such as particle size distribution, density, static friction coefficients, bulk density, and angle of repose, were determined through physical experiments. Subsequently, typical particle models within different particle size ranges were established and the particle sizes were magnified by a factor of five using particle scaling theory principles. This approach aimed to shorten simulation time and reduce the computational demands on the computer performance during simulations. Following this, the Plackett-Burman (P-B) test was employed to analyze the significance of calibration parameters, including Poisson's ratio, shear modulus, rolling friction coefficient, restitution coefficient and Johnson-Kendall-Roberts (JKR) surface energy. The steepest ascent test was then utilized to quickly determine the range of optimal parameter combinations for the simulated angle of repose. Subsequently, a quadratic regression equation linking the significance parameters to the angle of repose was established by Box-Behnken (B-B) test.

    Results and Discussion With the optimization objective of minimizing the relative error between the experimental and simulated angle of repose, the optimal parameter combinations are as follows: density at 1 326 kg/m3, Poisson's ratio at 0.4, shear modulus at 0.67 GPa, coal-coal static friction coefficient at 0.72, coal-steel static friction coefficient at 0.41, coal-coal rolling friction coefficient at 0.107, coal-steel rolling friction coefficient at 0.095, coal-coal restitution coefficient at 0.5, coal-steel restitution coefficient at 0.15 and the JKR surface energy is 14 J/m2. The simulated angle of repose under the optimal parameter combination is determined to be 37.59°, with a mere 0.58% relative error from the experimental value of 37.81°. Additionally, the simulated bulk density is 717 kg/m3, with a mere 0.69% relative error from the experimental value of 722 kg/m3.

    Conclusion 1) Among the calibrated parameters, the parameters that have a significant influence on the angle of repose are: coal-coal restitution coefficient, coal-coal rolling friction coefficient and coal-steel rolling friction coefficient in sequence. The influence of the coal-coal rolling friction coefficient on the angle of repose is particularly significant. 2) The primary terms of the three significant parameters, the quadratic terms of coal-coal restitution coefficient and coal-coal rolling friction coefficient, and the interaction terms of coal-coal restitution coefficient and coal-coal rolling friction coefficient have significant influence on the angle of repose. The influence of the primary terms of the coal-coal restitution coefficient and coal-coal rolling friction coefficient on the angle of repose are particularly significant. 3) The coal-coal restitution coefficient is negatively correlated with the angle of repose, while coal-coal and coal-steel rolling friction coefficients show a positive correlation with the angle of repose. Additionally, sputtering phenomena are observed to hinder particle packing. 4) The above simulation parameter combinations are only applicable to the simulation study of unscreened bituminous coals at a 5-fold enlargement of particle size. Since the primary terms of the coal-coal restitution coefficient and coal-coal rolling friction coefficient are known to have a particularly significant influence on the angle of repose, only the coal-coal restitution coefficient and coal-coal rolling friction coefficient need to be recalibrated using the same method when determining the simulation parameter combinations for unscreened bituminous coal at the remaining multiples of the particle size.

    Keywords: bituminous coal; discrete element; angle of repose; particle scaling theory; calibration of parameter

    (責(zé)任編輯:武秀娟)

    收稿日期: 2023-11-08,修回日期:2023-11-30,上線日期:2024-01-14。

    基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目,編號(hào):52105466; 上??莆纸瞬庞?jì)劃項(xiàng)目,編號(hào):21PJ1404600。

    第一作者簡(jiǎn)介:梅瀟(1974—),女,副教授,博士,碩士生導(dǎo)師,研究方向?yàn)楦劭跈C(jī)械的結(jié)構(gòu)設(shè)計(jì)與安全評(píng)估、焊接結(jié)構(gòu)的現(xiàn)代設(shè)計(jì)方法和多相場(chǎng)仿真。E-mail: xiaomei@shmtu.edu.cn。

    猜你喜歡
    煙煤
    煙煤特性及其對(duì)噴吹安全性的影響研究
    冶金能源(2024年1期)2024-02-05 06:28:54
    高爐混合噴吹上灣高鈣煙煤的性能研究
    冶金能源(2023年4期)2023-08-15 06:34:44
    百分百高爐煙煤比例成為國(guó)內(nèi)“首家”
    2020年泰國(guó)動(dòng)力煤進(jìn)口量同比增長(zhǎng)8.48%
    2月份泰國(guó)動(dòng)力煤進(jìn)口量環(huán)比增長(zhǎng)5.43%
    煙煤氧解耦化學(xué)鏈氣化及氮氧化物生成機(jī)理
    氣氛及后置催化劑對(duì)平朔煙煤熱解特性的影響
    能源工程(2021年2期)2021-07-21 08:39:46
    煙煤黏結(jié)指數(shù)測(cè)定準(zhǔn)確性的影響因素及操作要點(diǎn)探究
    山西化工(2021年6期)2021-01-26 12:56:21
    煙煤煙氣吸附劑脫汞技術(shù)的現(xiàn)狀及展望
    9月份泰國(guó)動(dòng)力煤進(jìn)口量同比增長(zhǎng)25.58%
    成人亚洲精品av一区二区 | 高清毛片免费观看视频网站 | 天堂俺去俺来也www色官网| 老鸭窝网址在线观看| 好看av亚洲va欧美ⅴa在| 欧美成人性av电影在线观看| 黄网站色视频无遮挡免费观看| 亚洲一区高清亚洲精品| 亚洲精品久久成人aⅴ小说| 国产成人av激情在线播放| 天天添夜夜摸| 美国免费a级毛片| 99热国产这里只有精品6| 精品乱码久久久久久99久播| 亚洲av成人一区二区三| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲男人的天堂狠狠| 老汉色av国产亚洲站长工具| 欧美精品啪啪一区二区三区| 波多野结衣一区麻豆| 两个人免费观看高清视频| 国产精品国产高清国产av| 成人三级做爰电影| 丝袜美腿诱惑在线| 久久久久久人人人人人| 桃色一区二区三区在线观看| 天堂动漫精品| 在线国产一区二区在线| 一区二区三区精品91| 日本欧美视频一区| 国产亚洲精品一区二区www| 久久这里只有精品19| a级片在线免费高清观看视频| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美日韩黄片免| 欧美黄色淫秽网站| 成年人黄色毛片网站| 搡老岳熟女国产| 国产单亲对白刺激| 男女做爰动态图高潮gif福利片 | 夜夜躁狠狠躁天天躁| 午夜福利在线免费观看网站| 精品人妻1区二区| 国产成人精品久久二区二区免费| 亚洲人成电影免费在线| 亚洲久久久国产精品| 久久精品亚洲熟妇少妇任你| svipshipincom国产片| 视频在线观看一区二区三区| 国产一区在线观看成人免费| 亚洲五月色婷婷综合| 黄片播放在线免费| 18禁黄网站禁片午夜丰满| 精品国产乱码久久久久久男人| 啦啦啦免费观看视频1| 女性被躁到高潮视频| 久久影院123| 精品无人区乱码1区二区| 电影成人av| 99精国产麻豆久久婷婷| 美女午夜性视频免费| 欧美日韩黄片免| 最新美女视频免费是黄的| 电影成人av| 亚洲精品一二三| 成人精品一区二区免费| 日韩一卡2卡3卡4卡2021年| 国产伦人伦偷精品视频| 巨乳人妻的诱惑在线观看| 中文字幕高清在线视频| 两人在一起打扑克的视频| 欧美激情高清一区二区三区| 日韩 欧美 亚洲 中文字幕| 高清欧美精品videossex| 中文字幕最新亚洲高清| 色精品久久人妻99蜜桃| 日本 av在线| 变态另类成人亚洲欧美熟女 | 亚洲精华国产精华精| 男女之事视频高清在线观看| 80岁老熟妇乱子伦牲交| 亚洲专区国产一区二区| 可以免费在线观看a视频的电影网站| 宅男免费午夜| 久久天躁狠狠躁夜夜2o2o| 99久久人妻综合| 国产精品香港三级国产av潘金莲| 多毛熟女@视频| 亚洲狠狠婷婷综合久久图片| 九色亚洲精品在线播放| 成人黄色视频免费在线看| 亚洲男人天堂网一区| 一边摸一边抽搐一进一出视频| 男人操女人黄网站| 在线天堂中文资源库| 亚洲国产精品合色在线| www.精华液| 色综合站精品国产| 久久久久久免费高清国产稀缺| 好男人电影高清在线观看| 亚洲中文日韩欧美视频| 18美女黄网站色大片免费观看| 色在线成人网| 国产熟女xx| 99精品欧美一区二区三区四区| 黄片小视频在线播放| 久久青草综合色| 狂野欧美激情性xxxx| 一级a爱片免费观看的视频| xxx96com| 欧美黄色片欧美黄色片| 免费看a级黄色片| 91国产中文字幕| 50天的宝宝边吃奶边哭怎么回事| 777久久人妻少妇嫩草av网站| 国产亚洲欧美精品永久| 精品电影一区二区在线| 亚洲国产中文字幕在线视频| 免费高清视频大片| 天堂中文最新版在线下载| 欧美精品啪啪一区二区三区| 久久99一区二区三区| 欧美av亚洲av综合av国产av| 欧美 亚洲 国产 日韩一| 丰满迷人的少妇在线观看| 久久久久九九精品影院| 亚洲成av片中文字幕在线观看| 成人永久免费在线观看视频| 精品国产一区二区久久| 男女下面进入的视频免费午夜 | 中文字幕最新亚洲高清| 亚洲一码二码三码区别大吗| 免费在线观看影片大全网站| 成人三级做爰电影| 久久久水蜜桃国产精品网| 又黄又粗又硬又大视频| bbb黄色大片| 国产一区二区三区综合在线观看| 无限看片的www在线观看| 18禁黄网站禁片午夜丰满| 日韩精品免费视频一区二区三区| 极品人妻少妇av视频| 看免费av毛片| 亚洲国产欧美一区二区综合| 涩涩av久久男人的天堂| 亚洲av成人一区二区三| 18禁黄网站禁片午夜丰满| 亚洲狠狠婷婷综合久久图片| 亚洲一区二区三区色噜噜 | 亚洲精品在线观看二区| 亚洲五月婷婷丁香| 久热这里只有精品99| 亚洲午夜精品一区,二区,三区| 国产亚洲欧美98| 精品久久久久久久毛片微露脸| 国产欧美日韩一区二区三区在线| 法律面前人人平等表现在哪些方面| 成年人免费黄色播放视频| 少妇 在线观看| 久久精品国产99精品国产亚洲性色 | 国产成人av激情在线播放| 最近最新中文字幕大全免费视频| 女人高潮潮喷娇喘18禁视频| 丁香六月欧美| 亚洲一卡2卡3卡4卡5卡精品中文| 精品免费久久久久久久清纯| 亚洲少妇的诱惑av| 国产1区2区3区精品| 国产精品 国内视频| 久久欧美精品欧美久久欧美| 91精品国产国语对白视频| 久久伊人香网站| 欧美成人性av电影在线观看| 久热爱精品视频在线9| 精品久久久久久电影网| 啦啦啦免费观看视频1| 丰满迷人的少妇在线观看| videosex国产| 亚洲精品国产区一区二| 久久影院123| 一区在线观看完整版| 国产精品自产拍在线观看55亚洲| 国产亚洲av高清不卡| 久久午夜综合久久蜜桃| 亚洲一码二码三码区别大吗| 美女高潮到喷水免费观看| 国产精品1区2区在线观看.| 婷婷丁香在线五月| 黄频高清免费视频| 丰满的人妻完整版| 成人国产一区最新在线观看| 在线观看舔阴道视频| 国产亚洲欧美98| 亚洲国产欧美一区二区综合| 成人av一区二区三区在线看| 国产欧美日韩一区二区精品| av在线天堂中文字幕 | 精品国产国语对白av| 久久久国产成人精品二区 | 国产成年人精品一区二区 | 黄色女人牲交| 纯流量卡能插随身wifi吗| 99久久99久久久精品蜜桃| 午夜成年电影在线免费观看| 欧美精品一区二区免费开放| 亚洲av日韩精品久久久久久密| 自线自在国产av| 国产一区二区三区视频了| 久久国产亚洲av麻豆专区| 嫁个100分男人电影在线观看| 亚洲视频免费观看视频| 久久精品亚洲精品国产色婷小说| 午夜久久久在线观看| 久久久久久久久久久久大奶| 最好的美女福利视频网| 搡老熟女国产l中国老女人| 亚洲人成电影免费在线| 中文字幕另类日韩欧美亚洲嫩草| 国产成人系列免费观看| 777久久人妻少妇嫩草av网站| 久久伊人香网站| 国产成人一区二区三区免费视频网站| 欧美成人性av电影在线观看| 岛国视频午夜一区免费看| 人人妻,人人澡人人爽秒播| av在线播放免费不卡| 精品高清国产在线一区| 国产高清videossex| 欧美人与性动交α欧美精品济南到| 啦啦啦 在线观看视频| 久久国产乱子伦精品免费另类| 女性被躁到高潮视频| 老熟妇乱子伦视频在线观看| 久久久久国产一级毛片高清牌| 美女高潮到喷水免费观看| 日本五十路高清| 热99国产精品久久久久久7| 成年女人毛片免费观看观看9| 午夜亚洲福利在线播放| 搡老岳熟女国产| 日韩免费高清中文字幕av| 国产av一区在线观看免费| 1024视频免费在线观看| 亚洲熟妇熟女久久| 一a级毛片在线观看| 曰老女人黄片| 国产在线精品亚洲第一网站| 成人手机av| 露出奶头的视频| 天堂中文最新版在线下载| 亚洲av电影在线进入| 国产不卡一卡二| 亚洲 欧美一区二区三区| 久久久久久亚洲精品国产蜜桃av| 99riav亚洲国产免费| 露出奶头的视频| 1024香蕉在线观看| 亚洲精品一区av在线观看| 免费观看精品视频网站| 亚洲成国产人片在线观看| 日韩大码丰满熟妇| 久久中文字幕人妻熟女| 日韩三级视频一区二区三区| 久久天堂一区二区三区四区| 韩国精品一区二区三区| 精品少妇一区二区三区视频日本电影| 精品久久蜜臀av无| www.www免费av| 亚洲av熟女| 午夜福利免费观看在线| 精品乱码久久久久久99久播| 国产不卡一卡二| 女人被躁到高潮嗷嗷叫费观| 中文亚洲av片在线观看爽| 在线观看66精品国产| 一边摸一边抽搐一进一出视频| 亚洲欧美日韩高清在线视频| 精品一区二区三区四区五区乱码| 久久精品91无色码中文字幕| 免费高清在线观看日韩| 深夜精品福利| 精品国产一区二区三区四区第35| 国产精品影院久久| 老司机靠b影院| 又紧又爽又黄一区二区| 一二三四在线观看免费中文在| 国产乱人伦免费视频| 欧美日韩瑟瑟在线播放| 中文字幕高清在线视频| 久久久久久久久中文| 国产高清国产精品国产三级| 视频区图区小说| 操出白浆在线播放| 亚洲精品美女久久av网站| 美女大奶头视频| 国产成人欧美| 色婷婷久久久亚洲欧美| 人妻丰满熟妇av一区二区三区| 老司机福利观看| 久久性视频一级片| 高清黄色对白视频在线免费看| 精品乱码久久久久久99久播| 日韩欧美一区二区三区在线观看| 在线观看www视频免费| 国产精品一区二区在线不卡| 国产97色在线日韩免费| 国产精品秋霞免费鲁丝片| 中文欧美无线码| 久久亚洲真实| 一级,二级,三级黄色视频| 天天添夜夜摸| 日本黄色日本黄色录像| 精品国产一区二区三区四区第35| 丝袜美足系列| 国产人伦9x9x在线观看| 亚洲一区高清亚洲精品| 热99re8久久精品国产| 国产国语露脸激情在线看| 日本vs欧美在线观看视频| 很黄的视频免费| 婷婷丁香在线五月| 免费久久久久久久精品成人欧美视频| 伦理电影免费视频| av天堂久久9| 久久久久久久久免费视频了| 精品久久久精品久久久| 亚洲一区中文字幕在线| 妹子高潮喷水视频| 日本免费一区二区三区高清不卡 | 99在线视频只有这里精品首页| 久久性视频一级片| 在线观看日韩欧美| 久久久久国产一级毛片高清牌| 久久草成人影院| 99久久国产精品久久久| 99精品久久久久人妻精品| 操出白浆在线播放| 三级毛片av免费| 久久久久国产一级毛片高清牌| 日本三级黄在线观看| 黄片播放在线免费| 国产精品 欧美亚洲| 亚洲欧美激情综合另类| 国产单亲对白刺激| 最近最新免费中文字幕在线| 久久99一区二区三区| 伦理电影免费视频| 国产视频一区二区在线看| 高清毛片免费观看视频网站 | 久久国产精品人妻蜜桃| 麻豆国产av国片精品| 国产色视频综合| 黄色片一级片一级黄色片| 女人被狂操c到高潮| 99久久综合精品五月天人人| 久久国产精品人妻蜜桃| 人人澡人人妻人| 成人影院久久| 人人澡人人妻人| 亚洲情色 制服丝袜| 一级毛片精品| 激情在线观看视频在线高清| 一级黄色大片毛片| 色哟哟哟哟哟哟| 好男人电影高清在线观看| 亚洲美女黄片视频| 美女午夜性视频免费| 男女床上黄色一级片免费看| 国产成人系列免费观看| 亚洲一区二区三区不卡视频| 9热在线视频观看99| 超碰97精品在线观看| 日韩免费高清中文字幕av| 午夜福利影视在线免费观看| 91老司机精品| 真人做人爱边吃奶动态| 欧美久久黑人一区二区| 熟女少妇亚洲综合色aaa.| www国产在线视频色| 热re99久久国产66热| a级毛片黄视频| 欧美一区二区精品小视频在线| 男人舔女人下体高潮全视频| 黄片小视频在线播放| 99在线人妻在线中文字幕| 欧美精品一区二区免费开放| 极品人妻少妇av视频| 亚洲精品国产色婷婷电影| 一级毛片精品| 亚洲午夜理论影院| 亚洲午夜精品一区,二区,三区| 国产精品综合久久久久久久免费 | 国产精品乱码一区二三区的特点 | 乱人伦中国视频| 精品福利观看| 日韩免费av在线播放| 午夜a级毛片| 国产精品一区二区三区四区久久 | 丰满饥渴人妻一区二区三| 好看av亚洲va欧美ⅴa在| 一本大道久久a久久精品| 欧美激情久久久久久爽电影 | 黄色视频,在线免费观看| 99国产精品免费福利视频| 午夜福利在线免费观看网站| 日本a在线网址| 国产欧美日韩一区二区三| 国产精品日韩av在线免费观看 | 中文字幕av电影在线播放| av国产精品久久久久影院| 丝袜美足系列| 久久香蕉精品热| 亚洲久久久国产精品| 欧美在线一区亚洲| 中文字幕人妻熟女乱码| 亚洲精品国产一区二区精华液| 免费一级毛片在线播放高清视频 | 最新在线观看一区二区三区| 精品久久蜜臀av无| 丰满迷人的少妇在线观看| 亚洲自偷自拍图片 自拍| 亚洲精品在线美女| av在线播放免费不卡| 91成年电影在线观看| 国产成人影院久久av| e午夜精品久久久久久久| 成人手机av| 欧美激情 高清一区二区三区| 欧美日本中文国产一区发布| 人成视频在线观看免费观看| 91精品三级在线观看| 黄频高清免费视频| 亚洲精品中文字幕一二三四区| 国产av又大| 国产精品永久免费网站| 韩国av一区二区三区四区| 国产视频一区二区在线看| 国产精品久久久av美女十八| 国产精品日韩av在线免费观看 | 亚洲一区二区三区欧美精品| 日韩视频一区二区在线观看| 交换朋友夫妻互换小说| 欧美丝袜亚洲另类 | 啦啦啦免费观看视频1| 中文字幕人妻丝袜一区二区| 久久99一区二区三区| 中文欧美无线码| 日韩欧美在线二视频| 看片在线看免费视频| 国产精品永久免费网站| 热re99久久国产66热| 久久久久久免费高清国产稀缺| 日韩欧美一区二区三区在线观看| 色精品久久人妻99蜜桃| 自拍欧美九色日韩亚洲蝌蚪91| 黄色a级毛片大全视频| 两性午夜刺激爽爽歪歪视频在线观看 | 在线观看日韩欧美| 国产亚洲欧美98| 午夜免费鲁丝| 欧美日韩精品网址| 国产深夜福利视频在线观看| a级片在线免费高清观看视频| 久久午夜亚洲精品久久| 色综合欧美亚洲国产小说| 亚洲欧美一区二区三区久久| 日韩免费av在线播放| 一区二区日韩欧美中文字幕| 国产91精品成人一区二区三区| 纯流量卡能插随身wifi吗| 真人做人爱边吃奶动态| 一级黄色大片毛片| 亚洲av片天天在线观看| 国产欧美日韩一区二区三区在线| 丝袜在线中文字幕| 国产精品1区2区在线观看.| 亚洲精品在线观看二区| 国产精品自产拍在线观看55亚洲| 久久久久久人人人人人| 午夜a级毛片| 啪啪无遮挡十八禁网站| 亚洲欧美精品综合久久99| 性欧美人与动物交配| 中文字幕高清在线视频| 久久久久久人人人人人| 极品教师在线免费播放| 成人三级做爰电影| 欧美乱色亚洲激情| 一区二区日韩欧美中文字幕| 久久久久国内视频| 91老司机精品| 女同久久另类99精品国产91| 免费av毛片视频| 制服诱惑二区| 可以在线观看毛片的网站| 欧美日韩黄片免| 精品一品国产午夜福利视频| 一级作爱视频免费观看| 两人在一起打扑克的视频| 欧美日韩国产mv在线观看视频| 国产99白浆流出| 男女下面进入的视频免费午夜 | 欧美日本中文国产一区发布| 一级片免费观看大全| 制服人妻中文乱码| 一边摸一边做爽爽视频免费| 午夜福利一区二区在线看| 18禁裸乳无遮挡免费网站照片 | 妹子高潮喷水视频| 久久天堂一区二区三区四区| 一夜夜www| 啪啪无遮挡十八禁网站| 亚洲精品中文字幕一二三四区| 麻豆av在线久日| 最新在线观看一区二区三区| 香蕉久久夜色| 色播在线永久视频| 亚洲专区中文字幕在线| 国产一区在线观看成人免费| 免费在线观看视频国产中文字幕亚洲| 在线观看免费视频日本深夜| 亚洲欧美日韩无卡精品| 啪啪无遮挡十八禁网站| 母亲3免费完整高清在线观看| 自拍欧美九色日韩亚洲蝌蚪91| √禁漫天堂资源中文www| 午夜激情av网站| 国产蜜桃级精品一区二区三区| 色精品久久人妻99蜜桃| 欧洲精品卡2卡3卡4卡5卡区| 一级毛片女人18水好多| 国产精品亚洲一级av第二区| 人人妻,人人澡人人爽秒播| 欧美激情 高清一区二区三区| 这个男人来自地球电影免费观看| 91av网站免费观看| 老司机午夜福利在线观看视频| 99精品在免费线老司机午夜| 精品高清国产在线一区| 亚洲成人久久性| 午夜视频精品福利| 美女高潮到喷水免费观看| 热99re8久久精品国产| 亚洲国产精品合色在线| 日韩成人在线观看一区二区三区| 国产亚洲av高清不卡| 亚洲片人在线观看| 欧美日韩亚洲综合一区二区三区_| 高清av免费在线| 亚洲成人免费电影在线观看| 淫秽高清视频在线观看| 国产精品日韩av在线免费观看 | 久久 成人 亚洲| 精品久久久久久,| 最近最新免费中文字幕在线| 国产av一区在线观看免费| 变态另类成人亚洲欧美熟女 | 女人被躁到高潮嗷嗷叫费观| 亚洲 国产 在线| 久久久久久久久久久久大奶| 丰满人妻熟妇乱又伦精品不卡| 黄色成人免费大全| 一二三四社区在线视频社区8| 国产麻豆69| 久久青草综合色| 国产精品1区2区在线观看.| 欧美在线黄色| 他把我摸到了高潮在线观看| 日日干狠狠操夜夜爽| 90打野战视频偷拍视频| 天天添夜夜摸| 久久久久久大精品| 亚洲人成电影免费在线| 免费在线观看日本一区| 男女午夜视频在线观看| 久久精品国产亚洲av高清一级| 国产av一区在线观看免费| 精品国产亚洲在线| 亚洲 国产 在线| 久久人妻熟女aⅴ| 国产国语露脸激情在线看| 搡老乐熟女国产| av免费在线观看网站| 亚洲成av片中文字幕在线观看| 精品国产一区二区久久| 黑人猛操日本美女一级片| 不卡av一区二区三区| 91麻豆精品激情在线观看国产 | 在线视频色国产色| 在线观看一区二区三区| 黄色视频不卡| 宅男免费午夜| 午夜福利在线观看吧| 热99re8久久精品国产| 又紧又爽又黄一区二区| 天堂动漫精品| 精品无人区乱码1区二区| 亚洲精品一卡2卡三卡4卡5卡| 美女高潮到喷水免费观看| 91国产中文字幕| 99久久国产精品久久久| 亚洲精品一卡2卡三卡4卡5卡| 99国产综合亚洲精品| 久9热在线精品视频| 国产成+人综合+亚洲专区| 一区二区日韩欧美中文字幕| 高清毛片免费观看视频网站 | 国产精品亚洲一级av第二区| 亚洲成av片中文字幕在线观看| 中文字幕色久视频| 女人被躁到高潮嗷嗷叫费观| 欧美在线黄色| 日韩人妻精品一区2区三区| 久久久久久久久免费视频了| 女人被躁到高潮嗷嗷叫费观|