999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

高速列車客室內呼吸道飛沫傳輸規律的數值模擬研究

2023-07-06 01:01:44徐任澤伍釩劉歡錢博森
中南大學學報(自然科學版) 2023年5期

徐任澤,伍釩,劉歡,錢博森

(1. 中南大學 交通運輸工程學院,軌道交通安全教育部重點實驗室,湖南 長沙,410075;2. 中南大學 軌道交通安全關鍵技術國際合作聯合實驗室,湖南 長沙,410075;3. 中南大學 軌道交通列車安全保障技術國家地方聯合工程研究中心,湖南 長沙,410075;4. 中南大學 湘雅三醫院,湖南 長沙,410013)

過去幾十年中,諸如肺結核、流感、埃博拉等呼吸道傳染性疾病一直威脅公共健康安全。隨著時間推移,全球范圍內的感染人數仍急劇上升,公共健康安全問題并沒有得到有效解決。空氣傳播和接觸傳播是呼吸道傳染性疾病傳播的主要途徑[1],患者通過不同的呼吸活動(說話、咳嗽、噴嚏)釋放出大量攜帶病原體的飛沫(飛沫核),這些飛沫(飛沫核)易懸浮在空氣中或者沉積在墻面上,致使人體暴露在含有病原體的公共環境中[2]。伴隨著城市化進程的加快,世界各國對軌道交通的需求急速上升。高速列車因其具有舒適性好、運行速度快、覆蓋范圍廣等優點而受到青睞[3-4],尤其是對于發展中國家,高速列車的重要性更加突出。例如,在每年重要節日時,列車內的乘客密度都會達到頂峰。然而,運營列車是一個高度封閉和擁擠的室內空間,這可能會進一步加劇呼吸道傳染病交叉感染現象[5]。因此,需要迫切研究高速列車客室內飛沫的流動特征,為設計有效的通風策略提供依據,這對于軌道交通持續發展和乘員健康安全極為關鍵。

迄今為止,許多學者通過理論分析[6-7]、數值模擬[8-10]和試驗[11-12]等方法研究了飛沫在不同環境下的擴散特征。YAMAKAWA 等[13]研究了攜帶病原體飛沫在教室中的長期擴散過程,并考慮病原體的衰減性,發現教室中的社交距離應大于建議值(2.0 m)的2倍以上。YAN等[14]發現咳嗽氣流和人體熱量均會顯著增加飛沫在空氣中的停留時間和移動距離,這導致高占用率的室內空間中人員被感染的風險性急劇增加。YANG等[15]通過數值模擬方法研究了室外環境條件對飛沫擴散的影響,發現室外風速和相對濕度增加會導致飛沫沉降到地面的數量明顯上升。QIAN等[12]以N2O示蹤氣體代替飛沫核,研究了3 種室內通風方案(混合、向下和置換通風)下飛沫核的擴散特性,發現與其他2種通風方案相比,置換通風方案更容易導致呼出的飛沫核聚集在局部空間中。總體而言,人們對醫院、房間、飛機等常規環境中的飛沫或者示蹤氣體擴散規律進行了大量研究,而對高速列車內的飛沫流動特征以及通風策略研究較少。為了深入研究高速列車客室內多場(速度、溫度等)耦合作用下飛沫的空間運動特性,本文作者在高還原度的客室內部空間中建立飛沫運動、擴散的數值模擬方法,分析釋放位置、初始粒徑和初始速度等因素對客室中飛沫傳播過程的影響。構建適用于列車客室流場與飛沫耦合運動的數值模擬方法,揭示不同邊界條件下飛沫的擴散規律,以便為降低客室內呼吸道傳染性疾病傳播和交叉感染風險提供參考。

1 客室內飛沫運動特性的數值模擬方法

1.1 主要假設

為模擬飛沫在高速列車客室內部空間中的傳輸過程,進行如下假設:1) 高速列車客室的氣密性良好,不考慮室內的滲風、漏風等;2) 飛沫顆粒之間不發生凝并、破碎現象;3) 飛沫顆粒在整個運動和蒸發過程中均保持球形;4) 僅考慮人體口腔產生的飛沫,即乘客口部為飛沫的噴射出口;5) 數值模擬中飛沫受力只考慮對其運動起主導作用的重力、氣動阻力、熱泳力以及氣動升力,忽略其他作用力[14-15];6) 人體衣服和客室材料具有高度的吸附性,當飛沫接觸到人體或者客室表面時,出現完全沉積現象,此時,停止該飛沫的追蹤計算。

1.2 幾何模型、邊界條件及網格設置

以某型號高速列車中間車為研究對象,采用比例為1∶1 的計算模型。該模型由客室乘坐區域、兩端連廊和完整的風道系統等構成,未建立空調單元的實際模型[16-18]。計算模型的主要幾何參數和結構分布如圖1(a)和圖1(b)所示,其中,H為包含風道結構在內的列車高度(H=2.78 m),列車總長為8.92H,寬度為1.02H。列車的風道系統由送風道、回風道和廢排風道組成,在圖1(c)中用不同的顏色標記風道系統各部分的構成。經過空調單元處理的空氣通過送風管道從2 種不同類型的出風口(第一種類型的出風口位于行李架上方,第二種類型的出風口位于座椅下方)進入客室中,客室中的部分空氣通過回風管道被吸入到空調單元混合腔,然后與來自外界環境的新風混合,而其余部分的空氣則通過廢排風道直接排放到外界環境中。此外,根據GB 10000—1988 標準[19],在5(列數)×18(行數)個座位上建立乘客坐姿模型(圖1(b)中數字1~18 為座位序號),該模型的適用性已在文獻[17]中得到驗證。

圖1 通風系統和客室幾何結構Fig.1 Geometric models of ventilation system and passenger compartment

在客室流場邊界條件方面,將送風口、回風出口以及廢排出口設置為體積流量邊界,總送風量為4 800 m3/h,回風量為3 000 m3/h,其余風量從廢排出口流出。以夏季作為列車運行環境(外環境溫度為35 ℃),空調送風溫度為12 ℃,送風濕度為70%。客室側壁設置為對流換熱邊界,對應的隔熱系數(K)為1.15 W·m-2·K-1。此外,在模擬中還考慮乘客熱量的影響,將36.5 ℃的溫度熱源應用于乘客表面[17]。在飛沫污染物邊界條件方面,將乘客口部看作1個面射流源,人體從口部呼出的飛沫污染物為連續相(空氣)和離散相(飛沫)的混合物,借鑒目前已得到充分驗證的飛沫組分模擬方法,將人體呼吸道飛沫假設為由98.2%的水和1.8%的非揮發性組成的固體化合物[20]。由于飛沫顆粒與空氣流之間存在熱量和質量傳遞過程,因此,飛沫從人體釋放后立即進行蒸發,當其粒徑不再變化時,即表明飛沫蒸發過程完成,將以固態顆粒物的形式繼續在客室流場中運動。在每個計算工況中共釋放15 000 個粒子模擬飛沫污染物的運動,15 000個粒子能有效降低由粒子數量導致的飛沫空間分布的差異性問題[21]。另外,為了精確地捕捉到客室流場變化和飛沫動力學特征,時間步長設置為0.001 s,每個計算工況的總模擬時間為15 s。

本研究采用四面體的非結構化網格對模型的幾何結構進行離散化,車身的基礎網格邊長為40 mm,風道區域的參考網格邊長為3.0 mm,并在速度梯度較大的人體周圍區域以及風道結構設置了加密的附面層。而對于速度梯度變化較小的區域(例如行李架附近),由于相關物理信息變化平穩,網格數量可以適當減少。同時,建立了3種不同數量的網格來進行網格無關性測試,3套網格的總數量分別為5 200萬個(細網格)、4 100萬個(中網格)和2 700 萬(粗網格)個。當網格數量達到4 100萬個時,客室流場計算結果的偏差小于1%。因此,為了節省存儲空間和提高計算效率,后續計算工況中將采用中等網格。

1.3 數值計算模型

客室流場由空氣和飛沫組成,本研究采用基于Eulerian-Lagrangian 方法的離散相模型進行數值模擬,空氣被設置為連續相,而飛沫被定義為離散相。在空氣構成的連續相方面,高速列車客室的空氣流馬赫數小于0.3[18],因此,選擇incompressible Navier-Stokes equations 作為流場的控制方程。Large eddy simulation(LES) and Reynolds average Navier-Stokes(RANS) turbulence models 是2種常見的計算模型。LES 模型的計算周期較長[5],而RANS 模型計算效率高、魯棒性強以及可靠性強,并廣泛應用于各種工程研究[16,22-23]。所以,本文采用RNGk-ε模型模擬高速列車客室內部的氣流組織。SIMPLE算法用于速度-壓力耦合,而k和ε方程采用second-order upwind scheme 進行離散以提高計算精度。在離散相方面,采用拉格朗日方法計算飛沫的釋放過程,可以準確捕捉到不同流場條件下飛沫的運動軌跡。

2 數值計算結果驗證

為了驗證本文數值模擬方法的準確性,將計算結果與相關文獻中的飛沫蒸發和運動距離進行對比[6-7,14]。具體而言,在CFD 模擬中進行相同的設置(如幾何模型和環境條件等)。人體位于空氣流速為0 m/s 的空間中,以10 m/s 的速度從人體嘴部釋放飛沫,嘴部的有效直徑為2 cm(目的是為了保持飛沫粒子的雷諾數相同)。在90%濕度和25 ℃恒溫的環境中分別模擬2 種典型初始粒徑dp0(10 μm和100 μm)飛沫的動態過程,時間步長為0.001 s。圖2(a)所示為飛沫粒徑隨時間的變化曲線,并與WEI 等[6]以及YAN 等[14]的研究結果進行比較,計算結果與文獻[6, 14]報道結果基本一致,從而驗證了本文使用的蒸發模型的準確性。此外,本文還模擬了粒徑為100 μm 的飛沫在咳嗽射流作用下的擴散特性,并與LIU等[7]的研究結果進行對比。圖2(b)中黑色虛線為咳嗽射流邊界,黑色點為文獻中理論計算結果,紅色點為本文數值計算結果,發現本文的數值模擬方法可以很好地預測大部分飛沫的擴散距離,部分誤差是湍流影響下飛沫運動的隨機性造成的。

圖2 飛沫蒸發模型和擴散模型的驗證Fig. 2 Validation of droplet evaporation and dispersion

總體而言,本研究的模擬結果基本準確,湍流條件下的飛沫擴散和蒸發模型是有效的,并應用于后續的客室環境分析。在客室內流場模擬方面,本文的研究方法已在實車試驗中得到了驗證[16-18]。因此,本文所建立的數值計算方法可以有效地模擬客室氣流組織與飛沫的兩相耦合運動過程。

3 客室內飛沫傳輸特征分析

3.1 客室流場分析

為了得到列車客室內部流場特征,選取客室典型位置的橫截面,包括送風出口截面(行李架上方以及座椅下方,具體位置可見圖1(c))和回風入口截面,各橫截面上的速度流場分布如圖3 所示。當客室中間過道的氣流向下運動時,逐漸往兩側壁面擴散,而兩側壁面的氣流從客室頂部往中間過道運動,從而在乘客前面形成2個渦流。客室左右兩側的速度流場均高度對稱,各位置局部流場存在差異的主要原因在于:1) 圖3(b)中座椅下方的送風出口會導致客室底部區域形成2個更加明顯的渦流;2) 圖3(c)中回風入口使得客室左右兩側的渦流中心在垂向上移動到更高位置。

圖3 客室內部空間中各橫截面上的速度和流線分布Fig. 3 Velocity distribution and streamlines of different locations in passenger compartment

3.2 釋放位置對飛沫擴散的影響

為了分析釋放位置對客室內飛沫顆粒流動的影響,分別選取3個位置作為飛沫釋放源。在縱向方向上(x方向,即車長方向),以座位3C 和座位7C 作為飛沫釋放源,而在橫向方向上(y方向,即車寬方向),以過道兩側的座位3C 和座位3D 作為飛沫釋放源。每個計算工況僅1 名乘客釋放飛沫,假設飛沫初始速度為10 m/s(近似于咳嗽的平均氣流速度),持續時間為0.5 s[7],釋放的飛沫初始粒徑為50 μm。

由于每次從乘客口部釋放的飛沫粒子數量較多,各粒子的流動狀態(如蒸發過程、位移特征以及各個方向上的最大傳輸距離等)可能出現差異。根據LI 等[24]的研究成果,本文以粒子蒸發時間、運動距離的中位數來表征飛沫群的整體運動狀態。

不同位置釋放的飛沫蒸發過程的影響如圖4所示。從圖4 可見:隨著飛沫在空氣中的運動時間(0~8 s)增加,其粒徑減小的速率增大;8.0~8.5 s時,各位置釋放的飛沫蒸發完成時間接近,這種略微差異是客室流場的不均勻性所致。

圖4 釋放位置對飛沫蒸發過程的影響Fig. 4 Effect of release location on evaporation of droplet

客室內飛沫的傳播軌跡隨時間的變化如圖5所示。釋放位置為座位3C,飛沫顆粒表面顏色對應它們在空氣中的運動時間(既從人體釋放后的運動時間)。從圖5可見:在前1 s內,乘客釋放的飛沫主要沿縱向傳播,飛沫的沉降運動相對不明顯(圖5(a)),這是由于咳嗽氣流增加了飛沫的初始速度,因此,在運動初始階段,重力和來自周圍流場的氣動力對飛沫傳輸的影響很小。然而,當飛沫進一步擴散和蒸發時,來自周圍流場的氣動力將對飛沫運動起主導作用,中間過道的向下氣流導致座位3C釋放的飛沫加速沉降(圖5(b))。隨著客室內部流場對飛沫傳播的進一步作用,座位3C 釋放的飛沫群出現了2個不同的擴散趨勢,一部分飛沫繼續以下降方式運動,而另一部分飛沫則跟隨周圍氣流向上運動(圖5(c))。并且當運動時間達到15 s時,約5.8%的飛沫流入回風管道中(圖5(d)),這部分飛沫可以通過空調單元和送風管道再次進入客室各區域中,從而威脅乘員健康。

圖5 客室內飛沫傳播軌跡隨時間的變化過程Fig. 5 Transmission processes of droplets in the passenger compartment

客室內不同位置釋放的飛沫擴散過程如圖6所示,可見各位置釋放的飛沫幾乎都分布在周圍的局部空間中。造成該現象的原因是:一方面,風道系統采用“分散式”回風方式,客室內各回風入口均勻分布在座位上方的行李架附近,減小了客室氣流的縱向流動;另一方面,飛沫在擴散過程時,蒸發導致重力作用減小,飛沫周圍的客室氣流開始主導其運動狀態,最終致使大部分飛沫在周圍的局部空間中流動。

圖6 客室內不同位置(座位)釋放的飛沫擴散過程Fig.6 Dispersion of droplets released from different locations(seats) in passenger compartment

為了進一步研究飛沫群在釋放之后的整體擴散特性,對飛沫的軌跡長度(L)和位移距離(D)的比值進行比較,發現約88%的飛沫的L/D大于2.0,這表明大部分飛沫在乘客附近循環運動,不會以長距離輸運方式遠離乘客。不同位置釋放的飛沫在客室內部空間中的數量占比見表1。據表1 可將客室內部空間在橫向上劃分為3個區域,即粒子釋放側的座位區域、中間過道區域以及非粒子釋放側的座位區域。非釋放側區域和過道區域內粒子數量之和的最高占比僅為4.85%,因此,大多數飛沫在釋放側一端流動,極少部分粒子會流動到過道以及非釋放側區域。綜合上述粒子軌跡以及空間數量分布特征,可以認為呼吸道飛沫主要以局部循環流動方式在該類型高速列車客室中運動。

表1 不同位置釋放的飛沫在客室內部空間中的數量占比Table 1 The proportion of droplets distribution in the passenger compartment

3.3 初始粒徑對飛沫運動的影響

飛沫初始粒徑與飛沫所受重力直接相關,在不同重力作用下,飛沫運動軌跡截然不同,因此,初始粒徑是影響飛沫傳輸規律的重要因素。本節采用普遍關注的粒徑(10、50、100、150 和200 μm)進行研究。不同初始粒徑的飛沫在咳嗽射流作用下從座位3C 釋放,其粒徑隨時間的變化過程如圖7所示。從圖7可見:小粒徑(10 μm)飛沫經歷極為快速的蒸發過程(0.1 s 內),而中等粒徑(50 μm)飛沫則需8.5 s 左右完成蒸發過程,這表明小、中粒徑飛沫能被上升氣流(人體熱羽流和客室氣流)帶入乘客呼吸區域,從而導致乘員感染風險顯著增加;而在更大的重力荷載作用下,大粒徑(大于100 μm)飛沫轉變為飛沫核前就會沉積在乘客和客室表面,表現出完全不同的蒸發特征,并且其減小的粒徑不足初始粒徑的15%。

圖7 不同初始粒徑飛沫的蒸發時間曲線Fig. 7 Variation of droplet diameter with different initial sizes

為了反映不同初始粒徑飛沫的沉降特性,研究客室內部空間中的飛沫運動軌跡,如圖8 所示。從圖8可見:對于初始粒徑小于50 μm的飛沫在經歷蒸發過程后,大部分飛沫會跟隨氣流聚集在乘客腿部高度以上的區域;而對于粒徑大于100 μm的飛沫,客室通風氣流失去了主導作用。不同初始粒徑的飛沫在垂向方向上的移動過程如圖9 所示,以乘客口部中心為相對坐標系的原點(乘客口部與客室地板之間的垂直距離為1.1 m)。從圖9可見:小粒徑飛沫和中等粒徑飛沫在垂向上的位移特征相似,在同樣的時間內,位移方式均為先下降再上升,位移方式改變的時刻發生在第3 秒左右,這是因為飛沫在橫向上是往靠近釋放側一端的車壁運動,從而更接近乘客前方渦流的上升側(圖3);大部分小、中粒徑飛沫長時間懸浮在乘客口部上方區域,其上升距離能夠達到0.6 m 以上,這表明乘客口部上方區域是病原體飛沫空氣傳播的高危區域。然而,當飛沫初始粒徑大于100 μm時,重力作用主導飛沫的傳輸過程,客室通風氣流的影響下降十分明顯,大粒徑飛沫無上升過程且快速(6 s內)沉積在乘客和座椅表面,導致病原體聚集,接觸感染的風險增加。

圖8 客室內不同初始粒徑飛沫的運動過程Fig.8 Dispersion processes of droplets with different initial diameters

圖9 不同初始粒徑飛沫隨時間變化的垂向位移Fig.9 Vertical traveling displacement of droplet group with different initial diameters with time

沉降特性對飛沫在客室中的傳播時間產生較大影響。飛沫在客室中的最大傳輸距離如圖10 所示。從圖10 可見:小粒徑飛沫和中等粒徑飛沫在三維方向(縱向、橫向、垂向)上的最大移動距離較接近,這是因為小、中等粒徑飛沫的擴散特征與客室氣流流動特征相似。大粒徑飛沫在縱向上的移動距離主要是由咳嗽氣流以及自身慣性所致,其最大移動距離為0.8 m左右,與前后相鄰座位之間的距離相近。而在橫向上,當初始粒徑從10 μm增加到200 μm 時,飛沫在空氣中的運動時間急劇減小(從15 s 下降到1.8 s,圖7),導致飛沫在橫向上的最大運動距離從1.73 m 下降到0.18 m。此外,由于大粒徑飛沫的重力作用大于客室通風氣流作用,大粒徑飛沫在垂向上僅存在下降過程。因此,大粒徑飛沫垂向上的最大移動距離(1.1 m)等于乘客嘴部到地板之間的垂直距離。

圖10 客室內不同初始粒徑的最大傳輸距離Fig.10 The maximum traveling distance of droplets with different initial diameters

3.4 初始速度對飛沫分布的影響

根據已有研究成果[2,7],本文假設人體不同呼吸活動(說話、咳嗽、噴嚏)釋放的氣流分別為1.0、10 和和20 m/s 的穩定射流。不同呼吸方式對飛沫蒸發過程的影響如圖11 所示,飛沫初始直徑為50 μm,其釋放位置為座位3C。從圖11 可見:在前2.4 s 內,當人體呼出氣流速度越快時,飛沫蒸發速率越大。這是因為更大的相對速度加速了飛沫和周圍空氣的傳熱傳質過程。然而,在后續的運動時間內,乘員以說話方式釋放的飛沫蒸發速率均比其他2種呼吸方式的飛沫蒸發速率大,飛沫蒸發曲線出現了2個拐點。導致該現象的主要原因是乘員以說話方式釋放的飛沫在人體附近區域運動,人體熱效應導致其局部空間溫度高于客室其他區域溫度(送風溫度為12 ℃),從而減少了飛沫蒸發完成的時間。

圖11 不同呼吸活動下飛沫蒸發過程的時間變化曲線Fig.11 Variation curves of droplet evaporation processes with time under different respiratory activities

初始階段下人體呼出氣流對飛沫傳輸的影響如圖12所示。從圖12可見飛沫以較低速度(1.0 m/s)被釋放時,飛沫聚集在乘客口部前方的0.15 m內;隨著釋放速度進一步增大,飛沫向釋放位置前方的座椅靠近,其縱向上的運動距離可以達到0.740 m;而當乘客以噴嚏方式釋放出氣流和飛沫時,受到釋放位置前方座椅表面的限制和客室通風氣流的影響,乘客噴嚏釋放的射流效應減弱,導致一部分飛沫向人體高度以下的空間運動。15 s內不同初始速度的飛沫在客室內部空間中的擴散軌跡如圖13所示。從圖13 可見:整體上看,當乘員以1.0 m/s(說話)釋放出飛沫時,飛沫向前傳播的距離較短,易受到客室局部氣流的影響,飛沫主要在乘客口部上方空間中運動;而當飛沫初始速度達到20 m/s(咳嗽)時,飛沫在人體高度以下的縱向傳播距離明顯增加,這表明人體呼吸活動釋放的氣流速度越大,會加劇攜帶病原體飛沫在客室內部空間中的擴散范圍。

圖12 1 s內不同呼吸活動對飛沫傳輸過程的影響Fig. 12 Effect of respiratory activity on droplets’transmission during 1 s

圖13 15 s內不同呼吸活動所釋放的飛沫擴散軌跡Fig.13 Dispersion of droplets released by different respiratory activities within 15 s

呼吸活動對飛沫最大傳輸距離的影響如14 所示。從圖14 可見:在運動初始階段,人體呼吸活動釋放的氣流速度對呼吸道飛沫傳輸速度有明顯影響,特別是縱向方向上;而在6 s 后,咳嗽和噴嚏的射流效應明顯減弱,飛沫傳輸速度主要受到客室通風氣流的影響;而在15 s時,隨著人體呼吸活動釋放的氣流速度從1.0 m/s(說話)增加到20 m/s(咳嗽),飛沫縱向移動的最大距離從0.89 m提高到5.23 m,橫向移動的最大距離從0.76 m提高到2.73 m,垂向移動的最大距離從1.21 m 提高到2.53 m。因此,相比于說話方式,噴嚏能夠導致飛沫在三維(縱向、橫向和垂向)方向上的最大移動距離分別增加4.87、2.59和1.09倍。

圖14 呼吸活動對飛沫最大傳輸距離的影響Fig.14 Effect of respiratory activity on the maximum transmission distance of droplet

4 結論

1) 在送風系統、客室內部結構以及乘客熱量的共同作用下,客室流場會產生輕微的不均勻特征,導致不同位置釋放的飛沫蒸發速率存在差異。

2) 對于本研究所涉及的高速列車結構,各位置釋放的飛沫在客室內主要以局部循環運動為主。在客室兩側分布的座位中,由于客室流場具有良好的對稱性,大多數飛沫分布在釋放者一側的空間中(數量占比超過95%)。建議客室每排座位的左右側均只開放1個座位,這樣可減少呼吸道傳染性疾病交叉感染風險。

3) 由于重力效應,大粒徑(大于100 μm)飛沫未完全蒸發就沉積在物體表面,而小粒徑(10 μm)飛沫則由于蒸發速率快(蒸發完成時間小于0.1 s),長時間內跟隨客室氣流一起運動。因此,與大粒徑飛沫相比,小粒徑飛沫的懸浮時間和擴散范圍均更大,從而導致病原體飛沫空氣傳播的感染風險更高。

4) 人體呼吸方式產生的氣流速度對飛沫傳播特性有很大影響。相比于說話釋放的飛沫傳播范圍(縱向、橫向、垂向的最大移動距離分別為0.89、0.76和1.21 m),噴嚏方式由于其強烈的射流效應,導致飛沫在三維(縱向、橫向和垂向)方向上的最大移動距離分別增加了4.87、2.59以及1.09倍。

主站蜘蛛池模板: 欧美成人午夜视频| 91蜜芽尤物福利在线观看| 国产一区亚洲一区| 国产三区二区| 日韩天堂在线观看| 亚洲黄色片免费看| 亚洲伦理一区二区| 成人午夜网址| 无码精油按摩潮喷在线播放 | 中国精品自拍| 亚洲无码熟妇人妻AV在线| 国产精品视频观看裸模| 激情亚洲天堂| 免费播放毛片| 日本人又色又爽的视频| 亚洲综合经典在线一区二区| 国产精品短篇二区| 成人福利在线视频| 国产成人三级| 黄色一级视频欧美| 精品超清无码视频在线观看| 无码又爽又刺激的高潮视频| 亚洲无码一区在线观看| 精品少妇人妻av无码久久| 伊人成人在线| 台湾AV国片精品女同性| 日韩AV手机在线观看蜜芽| 在线日韩日本国产亚洲| 欧美亚洲第一页| 激情视频综合网| 亚洲综合极品香蕉久久网| 免费一级毛片完整版在线看| 婷五月综合| 国产亚洲精品91| 四虎永久在线精品国产免费| 亚洲精品男人天堂| AV片亚洲国产男人的天堂| 成人精品午夜福利在线播放| 伊人久久大香线蕉aⅴ色| 一级香蕉人体视频| 国产成人91精品| 嫩草在线视频| 欧美成人日韩| 福利视频99| 狠狠久久综合伊人不卡| 在线观看网站国产| 97无码免费人妻超级碰碰碰| 丝袜亚洲综合| 国产成人精品免费视频大全五级| 青青青视频91在线 | 日韩视频免费| 久久夜色精品| 一本色道久久88综合日韩精品| 成人午夜免费观看| 日本一区二区不卡视频| 无码 在线 在线| 91丝袜乱伦| 色视频国产| 伊人久久福利中文字幕| 欧美曰批视频免费播放免费| 国产精品尹人在线观看| 亚洲狼网站狼狼鲁亚洲下载| 国产在线一二三区| P尤物久久99国产综合精品| 日韩免费视频播播| 日韩精品无码免费一区二区三区 | 日本欧美午夜| 中国毛片网| 亚洲人成人伊人成综合网无码| 亚洲无码高清视频在线观看| 成人小视频在线观看免费| 中文成人在线| 亚洲AV人人澡人人双人| 91成人免费观看| 国产一在线观看| 国产在线日本| 亚洲第一黄色网| 日韩在线播放中文字幕| 在线观看国产精品第一区免费| 九九久久精品免费观看| 伊人色综合久久天天| 全部无卡免费的毛片在线看|