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

基于阿里云的四維彈簧模型并行運算性能

2019-06-21 07:47:20趙高峰陳華
土木與環境工程學報 2019年3期
關鍵詞:裂紋計算機模型

趙高峰,陳華

(天津大學 建筑工程學院;水利工程仿真與安全國家重點實驗室,天津 300072)

由于具有可重復、經濟及參數可控等優點,數值模擬已經成為理論分析和物理實驗之外的第3種研究方法[1-3]。數值模擬不但被廣泛應用于隧道等地下工程在不同工況下的穩定性分析,而且被用來研究特定工程災變問題的內在力學機理。例如,唐春安等[1]采用RFPA軟件對隧道洞室周邊的分區破壞機理和演化規律進行了研究,吳順川等[2]采用三維離散元軟件研究了隧道巖爆的機理。上述兩個數值模擬案例分別采用了基于連續介質的有限元法和基于離散介質的離散元法。連續介質方法的基本思路是先建立求解對象的偏微分方程,然后,通過數學離散方法求解,是一種自頂向下的方法。由于其連續假設,連續介質方法在求解動態斷裂等問題時具有一定的局限性。將損傷模型引入連續介質方法可增強其求解漸進動態破壞的能力,例如RFPA軟件[1]及LSDYNA軟件[4]采用了類似單元生死法來實現對材料漸進破壞的描述,并已成功應用于巖石動態破壞及實際工程的計算。然而,這種關于破壞的處理方式沒有精確考慮破裂面的形態以及再接觸,因此,該方法對破壞后的描述可能會有偏差[5]。基于離散介質的數值方法則考慮了破裂面的分離和再接觸,更適合于巖石動態破壞問題,其中,最為著名的是Cundall等[6]提出的離散元和Shi[7]提出的非連續變形分析方法DDA。作為離散數值計算方法的一種,Lattice Spring Model (LSM)由亞歷山大博士于1941年最早提出[8],但由于泊松比限制問題,LSM發展一直很緩慢。為了解決該問題,Distinct Lattice Spring Model (DLSM)引入了多體剪切彈簧[8],該模型已被成功應用于巖石與煤的動態破壞研究[9-11]。最近,研究者提出了另一種解決方法,基于經典LSM引入額外維相互作用,稱之為四維彈簧模型(Four-lattice spring model,4D-LSM)[12]。4D-LSM的基本元件是由彈簧鍵連接的顆粒,通過彈簧鍵的變形和破壞來反映固體的宏觀力學響應。4D-LSM這種自底向上的建模方式與離散元類似,其單元(顆粒)數量必須達到一定規模才能得到足夠真實的模擬結果。因此,4D-LSM有龐大的計算需求,傳統的個人電腦已經不能滿足,高效的并行計算是唯一的解決辦法。

目前,并行計算器件主要分為多核CPU和GPU,而主要計算平臺是超級計算機。超級計算機擁有多個節點,每個節點一般是可以單獨實現CPU多核并行和(或)GPU并行的計算機,節點間通過網絡連接實現消息傳遞,從而將計算資源整合利用,并達到超高的計算性能。全球超級計算機Top500中,目前排名第一的超級計算機“Summit”由4 608個節點組成,每個節點搭載2個“Power9”CPU和6個“NVIDIA Tesla V100”GPU,CPU核心數量為202 752,GPU流處理器數量超過1.4億,其峰值性能為200PFLOPS[13]。近年來,GPU計算發展十分迅速[14],但是GPU計算并不能取代CPU計算,比如CPU更擅長處理邏輯控制密集任務,CPU多核并行仍然是一種便捷、可靠并且廣泛使用的高性能計算方式。實現CPU多核并行主要依靠應用程序接口,例如OpenMP(Open Multi-processing)等[15]。OpenMP是基于共享內存的應用程序接口,提供了對并行算法的高層抽象描述,非常適合多核CPU計算機的并行程序設計。OpenMP的顯著特點是精簡、易用,只需要在串行代碼中加入簡單的pragma指令即可實現并行,因此,OpenMP的使用非常普遍,例如4D-LSM和DLSM就采用OpenMP實現了多核并行[12,16]。

高性能計算通常以高性能計算機為依托,但超級計算機硬件的高昂費用和固定資產屬性常導致高性能計算的使用成本較高。近年來,計算領域中面向服務的“云計算”為解決這個問題提供了可能。云計算是指通過網絡按需提供虛擬計算資源和解決方案的有償服務,相對于傳統的計算模式,其主要優點是配置靈活、方便快捷、管理投入少等。例如,由“阿里云”提供的彈性云服務器類型有通用性、計算型、內存型等,CPU核數從2核到160核不等,內存從4 GB到1 920 GB不等[17],付費方式也有按量計費、按時間計費等不同選擇。筆者主要研究多核4D-LSM在云計算及常規多核工作站和個人電腦上的并行運算性能,通過大量數值模擬計算來研究線程數量、硬件配置、求解問題類型等對4D-LSM并行計算時間的影響,進而找到4D-LSM在“阿里云”計算環境下的極限規模和瓶頸,最后,通過4D-LSM求解脆性材料的三維破裂問題來展示離散數值計算方法和并行計算相結合帶來的優勢。

1 四維彈簧模型

1.1 基本原理

在經典物理學中,空間是三維的,時間作為第四維,它們共同構成了四維時空。有些研究者為了統一自然界的4種基本力,通過引入一個額外的空間維度,提出了五維時空。4D-LSM借鑒了五維時空理論。4D-LSM模型的構建過程如圖1所示。圖1(a)中,三維空間中的立方體晶格模型能夠再現各向同性彈性,其泊松比固定為0.25。該原始模型的彈簧鍵有兩種,即正彈簧(例如AB)和對角彈簧(例如AC),其剛度系數均為k。圖1(b)所示為原始模型在第四維的“平行體”,對于給定的質點A,其“平行體”即為A′,“平行體”模型的構造和彈簧剛度均與“本體”模型相同。然后,利用第四維相互作用(彈簧鍵)連接“本體”和“平行體”,如圖1(c),具體規則為:“本體”模型的一個彈簧(例如A-B)產生4個相應的四維彈簧(A-A′、B-B′、A-B′和A′-B)。

圖1 四維彈簧模型的構建過程[12]Fig.1 The model building process of

1.2 系統方程

在4D-LSM中,假定三維世界是一個四維超膜,離散的四維顆粒由彈簧鍵連接。4D-LSM的描述和證明詳見文獻[12],這里只關注對實現并行化必不可少的有關方程。四維顆粒的空間位置和運動參數表示為

xi=(xiyizi?i)T

(1)

(2)

(3)

(4)

式中:t為時間;Δt為時間增量。以相同的方法,可以得到顆粒的速度公式

(5)

顆粒i和顆粒j之間四維距離為

(6)

如果這些顆粒通過剛度為k的彈簧連接,那么顆粒j對顆粒i的作用力為

(7)

(8)

式中:mi為顆粒i的質量;gx、gy和gz為重力加速度。在4D-LSM中,假設牛頓第二定律也適用于第四維,則顆粒i的加速度為

(9)

式(1)~式(9)是實現4D-LSM并行化涉及的所有基礎性計算。

1.3 模型參數選取

對于立方體四維晶格,有3種類型的四維彈簧,剛度分別為kα、kβ、kγ。對于彈性各向同性體,它們的剛度值需滿足關系[12]

(10)

式中:λ4D為四維剛度系數;k為三維彈簧的剛度,k用式(11)計算。

(11)

式中:V為三維晶格模型的代表體積;E為彈性模量,li為三維晶格模型的初始彈簧長度;η為尺度參數。η可用式(12)計算[12]。

0.416 136 15λ4D+1.003 692 23

(12)

四維剛度系數λ4D也可以由泊松比得到[12]。

λ4D=-211.134 937 79v3+162.846 558 51v2-

55.424 497 19v+6.929 022 11

(13)

式中:v是泊松比。結合式(10)和式(13)可算出式(7)所需的力學參數(彈簧剛度)。這些參數都是預先計算的,與4D-LSM的計算循環無關,因此,參數計算部分不參與并行。更多細節和數學證明可以在4D-LSM的原始文獻[12]中找到。

1.4 OpenMP多核并行

圖2 4D-LSM的分叉并行策略Fig.2 Fork/Join parallel strategy of

1.5 計算模型

采用兩種4D-LSM計算模型,分別對應彈性問題和破壞問題。如圖3所示,模型外觀均為立方體,選用的顆粒直徑為1 mm。第1計算模型為立方體單軸壓縮試驗,底端在豎直方向被固定,頂端施加豎直向下的位移荷載,整個模擬過程不破壞,屬于線彈性問題。第2個問題是爆炸開裂模型,模型中心有球形空洞(綠色部分),沖擊荷載施加在球體的內表面,屬于動態破壞問題。每個計算模型采用不同規模,立方體邊長分別為20、50、100、150、200、250、300 mm,因此,最大的模型顆粒數達到2 700萬(300×300×300)。

圖3 兩種4D-LSM計算模型Fig.3 Two kinds of 4D-LSM computing

1.6 并行性能分析

云服務器是一種虛擬的計算機,根據客戶的需求可以有不同配置。如果用戶選用的云服務器操作系統與自己本地計算機的操作系統一樣(如Windows系統或者Linux系統等),那么,云服務器的操作體驗與本地計算機幾乎沒有區別,能夠在本地機運行的程序同樣可以在云服務器上運行,不需要做任何額外的更改,本文涉及的基于OpenMP的4D-LSM也是如此。如表1所示,選用的云服務器(CS)具有64核心、128 GB內存容量,CPU頻率是2.5 GHz。線程測試選用的是單軸壓縮模型,模型的大小有3種,邊長為50、100、150 mm,分別記為“Cube_50”、“Cube_100”和“Cube_150”,相應的顆粒數為12.5萬、100萬和337.5萬。將每一個模型在不同的計算機上采用不同的線程數進行重復計算,記錄每次計算的時長,并換算出加速比,加速比定義為單線程計算時間與多線程計算時間的比值。

表1 計算機主要參數Table 1 Main parameters of the computers

阿里云上的測試結果如圖4所示,最大加速比約為16.8×。單從最大加速比來看,云計算優于兩臺工作站WS-1和WS-2(見圖8,最大加速比約為9.0×),并且,云計算還有使用靈活、無需維護等優點。但是,相較于本地計算資源,云計算也有不足之處。首先是性價比的問題,測試用的云服務器按時間計費,費用約為350元/d,而工作站WS-2的一次性投入約為5萬元,該費用只能購買該云服務器5個月左右,但通常情況下,一臺工作站的性能至少可以在3年內保持相當的競爭力。數值計算方面的科研工作,經常需要修改模型的參數,這樣的重復計算是對云計算資源的浪費。因此,最合理的方式是利用本地計算資源調整數值模型,然后利用云服務一次性完成大規模計算。其次,大規模計算必定涉及到大量的數據存儲問題,由于云服務的存儲具有時效性,也不方便進行后處理工作,因此,如何快速將海量數據保存到本地存儲空間是云服務應用于數值分析計算面臨的另一個問題。

圖4 阿里云計算環境下4D-LSM并行效率測試結果Fig.4 Parallel efficiency test results of 4D-LSM in Alibaba

PC-1、PC-2、WS-1和WS-2的測試結果如圖5~圖8所示,最大加速比分別為2.6×、3.2×、10.8×和9.1×,隨著線程數的增加,計算速度總體上呈加快的趨勢。PC-1是4核8線程,由圖5可知,當線程數超過4以后,加速效果明顯下降,例如“Cube_150”的模型使用2線程、4線程和8線程時的加速比分別為1.84×、2.54×和2.6×,意味著加速比從2線程到4線程的增幅為38%,而從4線程到8線程的增幅僅為15%,這說明物理核心的加速效果遠遠超出超線程技術的加速效果。兩臺工作站WS-1和WS-2是20核40線程的雙路計算機,由圖7、圖8可知,當線程數超過20以后,加速效果的提升即開始放緩。更值得關注的是,對于WS-1和WS-2這兩臺雙路計算機(每個CPU有10個核心、兩個CPU共20核),10線程和20線程的加速效果幾乎相同,例如WS-1上“Cube_150”模型使用10線程和20線程時的加速比分別為7.1×和7.7×,當線程數介于10和20之間時,加速比呈現先降后升的“凹”型曲線,計算資源的增加卻適得其反。

圖5 4D-LSM在PC-1上并行效率測試結果Fig.5 Parallel efficiency test results of 4D-LSM on

圖6 4D-LSM在PC-2上并行效率測試結果Fig.6 Parallel efficiency test results of 4D-LSM on

圖7 4D-LSM在WS-1上并行效率測試結果Fig.7 Parallel efficiency test results of 4D-LSM on

圖8 4D-LSM在WS-2上并行效率測試結果Fig.8 Parallel efficiency test results of 4D-LSM on

1.7 并行性能影響因素及極限運算分析

1.7.1 求解類型的影響 4D-LSM模型中,破壞的顆粒在每一步的計算過程中都會進行動態接觸檢索,當有其他顆粒接觸到該破壞顆粒時,這兩個顆粒之間會產生一個新的特殊彈簧鍵,該彈簧鍵并不能受拉,其目的只是為了防止破壞顆粒由于運動而穿透其他顆粒,相對于非破壞模型,破壞模型的顆粒檢索將會消耗額外的時間。為了研究加速比與求解類型的關系,選用單軸壓縮模型和爆炸開裂模型進行對比,前者代表彈性(Elastic)問題,后者代表破壞(Fracture)問題,模型外觀均為立方體且邊長均為100 mm,計算機選用工作站WS-2,測試結果如圖9。對于線程數與計算效率的總體關系,破壞模型與前述彈性模型一致,但是,對比彈性模型與破壞模型的加速比發現,隨著線程數的增加,兩者的加速比差距越來越大,最終,使用40線程時彈性模型的加速比達到9.0×,而破壞模型相應的加速比為5.8×,僅為前者的64%。總之,并行化的4D-LSM求解破壞問題所獲得的加速效果要低于非破壞問題,使用的線程數越多,這種差距越明顯。

圖9 工作站WS-2對于彈性模型和破壞模型的并行效率測試結果Fig.9 Parallel efficiency test results for elastic model and failure model of 4D-LSM on

1.7.2 計算規模的影響 由圖7、圖8可知,對于大小不同的模型,在相同條件下,其加速比有一定的區別。例如工作站WS-1使用40線程時,“Cube_50”“Cube_100”和“Cube_150”的加速比分別為9.1×、8.0×和10.8×,而相同情況下,工作站WS-2上對應的加速比分別為9.0×、8.2×和7.7×。若僅就這3種大小的模型而言,則工作站WS-1上“Cube_150”加速效果最好(10.8×),而工作站WS-2上“Cube_50”才是加速效果最好的(9.0×)。因此,模型的大小對加速效果有一定的影響,但這種影響沒有普遍的規律,隨計算機硬件配置的不同而不同。

同時也發現,彈性模型的計算時間與模型的規模呈正比,而破壞模型則并非如此。對于破壞模型而言,顆粒檢索會消耗額外的時間,破壞的顆粒越多,每一步的計算時間就越長,但每一個破壞顆粒額外消耗多長時間還不明確,整個破壞過程目前也無法預測。由于這些“復雜性”,測試結果中破壞模型的計算時間與規模大小的關系曲線并不具有普適性,只能說明一般情況下是非線性的,從而區別于彈性模型的線性關系。

1.7.3 計算硬件的影響 CPU的主要參數包括頻率、核心數量和線程數量等,更高的CPU頻率、更多的核心或者線程都能夠獲得更快的計算速度。因此,由表1可知,一般情況下,擁有8線程3.6 GHz CPU的PC-1要比擁有4線程3.0 GHz CPU的PC-2更快,WS-1也會因為更高CPU的頻率而獲得比WS-2更好的性能。將圖4中關于計算時間的數據做進一步處理后得到表2,表中ΔPC、ΔWS分別為PC-1與PC-2、WS-1與WS-2計算同一模型所用時間之差。從表2來看,雖然有幾處Δ值為負數,但都是在模型較小、整個計算時間較短的情況下發生,不具有代表性,而絕大部分Δ值都為正數。因此,從統計的角度,對于同一個模型,可以認為PC-1比PC-2耗時更多,WS-1比WS-2耗時更多,也就是說,PC-2和WS-2計算速度更快,與之前的預測剛好相反,這說明4D-LSM的計算速度并非完全由CPU的性能決定。在表1中,對比4臺計算機的硬件,PC-2和WS-2唯一的優勢就是擁有更高的內存帶寬。由于計算時間不僅包括CPU處理數據的時間,也包括其他必要的時間消耗,如CPU和內存交換數據的時間,高內存帶寬意味著數據傳輸更快,最終的結果是PC-2和WS-2在計算時速度更快。因此,對于4D-LSM,若CPU性能差距不是很懸殊,則內存帶寬成為計算速度非常重要的影響因素。

表2 使用最大線程數計算不同大小模型的時間消耗表Table 2 Calculating time of different size models using maximum thread number

1.7.4 并行計算量瓶頸分析 將兩組模型(單軸壓縮、爆炸開裂)按從小到大的順序依次運行,記錄其計算時間以及消耗的物理內存,測試時,每臺計算機都使用最大線程,例如,PC-1使用8線程,而兩臺工作站WS-1、WS-2均使用40線程。從小到大的立方體模型的邊長為20、50、100、150、200、250、300、350、400 mm等,在此序列下,PC-1、PC-2、WS-1、WS-2能計算的最大模型邊長分別是150、150、250、300 mm,對應的顆粒數分別為337.5萬、337.5萬、1 562.5萬和2 700萬。由此可見,不論是彈性模型還是破壞模型,模型大小(顆粒數)與消耗的物理內存呈同一個線性關系。事實上,經過更進一步的數據分析發現,內存消耗與使用哪臺計算機也沒有聯系,即模型的顆粒數量與內存消耗存在一一對應的關系(如圖10所示),每100萬顆粒約需要1.8 GB內存,目前來看,4D-LSM的計算量由計算機的內存容量決定。例如,WS-1的內存容量是32 GB,當計算邊長為250 mm的模型時,4D-LSM消耗的內存約為28 GB,而下一個模型邊長是300 mm,顆粒數量2 700萬,按前述標準約需要48 GB的內存,因此,在WS-1上無法計算,最終,該模型在擁有64 GB內存的工作站WS-2上運行,而2 700萬顆粒也幾乎是WS-2的極限計算量。然而,對于阿里云來講,則可以最大運行10億單元的計算模型。從這點上來講,云計算為一些對顆粒規模要求十分龐大的問題提供了除傳統超級計算集群之外的可行解決途徑。相比傳統超級計算集群,云計算無需對代碼進行修改,也無需進行昂貴的硬件投資。然而,4D-LSM是采用自建的圖形交互界面建模,能夠建立多大的模型受制于顯存。例如,WS-2配備的“NVIDIA Quadro M5000”具有8 GB的顯存,其構建的最大模型是450×450×450(約9 000萬顆粒)。如果假設建模所需顯存與顆粒數成正比,構建10億顆粒的模型則至少需要大約88 GB的顯存,因此,目前4D-LSM大規模并行的瓶頸在于前處理。

圖10 模型大小與內存消耗的關系Fig.10 The relationship between model size and

2 應用案例

采用4D-LSM進行三維裂紋擴展分析。幣型裂紋試樣的尺寸及荷載條件如圖11(a)所示,裂紋形狀為圓形,直徑18 mm,厚度1 mm,中心位置與整個立方體試樣中心位置重合,裂紋面與試樣底面夾角θ=30°。建立兩個4D-LSM模型,一個解析度為110×110×110(約130萬顆粒),另一個解析度為220×220×220(顆粒數大約為1 060萬),除此之外,兩個模型并無其他任何差別。圖11給出了針對三維幣形裂紋的計算模型,采用并行4D-LSM進行了求解。圖12展示了三維幣型裂紋模型在不同解析度下的裂紋發展過程。通過對比,低解析度模型雖然能大致展現裂紋的擴展過程,但裂紋形態相對比較粗糙,裂紋擴展的對稱性遠不如高解析度模型。由此可見,更高解析度的4D-LSM模型對精準模擬三維裂紋擴展問題非常關鍵。基于云計算的并行計算技術可以求解更高解析度的計算模型,非常適合于求解三維裂紋擴展的計算。

圖11 三維幣型裂紋擴展模型Fig.11 The model of three-dimensional

圖12 裂紋擴展過程Fig.12 The crack propagation

3 結論

主要研究了4D-LSM在云計算環境下的并行性能,考慮了線程數量、硬件資源、模型大小、求解類型等因素。得到如下主要結論:

1)4D-LSM具有較好的并行性能,在20核的雙路計算機上的最大加速比接近11.0×,而在64核的云服務器上的加速比接近17×。

2)4D-LSM模型的規模對加速效果有一定的影響,并因使用的計算機不同而不同。

3)4D-LSM求解彈性問題的加速效果優于求解破壞問題,使用的線程數量越多,這種差別表現得越明顯。

4)非破壞模型的計算時間與顆粒數呈正比關系,而破壞模型由于其“復雜性”,通常情況下不是正比關系。

5)4D-LSM模型的顆粒數量與內存消耗呈正比,計算機的極限計算量由內存容量決定,每100萬顆粒大約需要1.8 GB的內存,若要求解10億顆粒的模型,理論上至少需要1.8 TB的內存。

6)對于雙路計算機應當注意,當線程數量介于單顆CPU的物理核心數和雙CPU的總物理核心數時,計算效率會下降,并且造成計算資源的浪費。

7)雖然云計算非常靈活且能提供強大的高性能計算能力,但其性價比也值得商榷,使用時應當綜合考慮,有的放矢。

另外,需要說明:

1)在測試極限計算量時只考慮了物理內存,實際上有些4D-LSM模型在內存需求超過計算機的物理內存時也可以計算,比如邊長為200 mm的單軸壓縮模型,顆粒數量是800萬,大約需要14.4 GB的內存,卻可以在內存容量8 GB的PC-2上運行,這是因為系統自動啟用了虛擬內存(此處虛擬內存是相對物理內存而言,并非編程模式下所指的虛擬地址空間),但此時計算速度非常緩慢,不在可接受的范圍,因此,未予以考慮。

2)對于雙路計算機,當使用的線程數量介于單顆CPU核心數和雙CPU總核心數時,不僅計算效率會下降,而且多次重復計算的結果表明:在此區間計算時間的離散程度也急劇增加,即計算效率不穩定,計算效率不穩定的情況與求解類型無關。相差40%的結論是因為統計了WS-2在10線程和20線程之間重復計算100次同一個模型的計算時間,求出了計算時間的變異系數,該變異系數最高為20%左右。因此,當線程數量處于該區間時,同樣模型的兩次計算時間有可能相差40%。

3)對于本地計算機,當使用的線程數量超出計算機的最大線程數時,其計算效率會下降20%左右,但對于云計算的虛擬服務器而言,超出最大線程后,計算效率并不會下降,而是保持在同一水平線上。

猜你喜歡
裂紋計算機模型
一半模型
裂紋長度對焊接接頭裂紋擴展驅動力的影響
計算機操作系統
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
基于計算機自然語言處理的機器翻譯技術應用與簡介
科技傳播(2019年22期)2020-01-14 03:06:34
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
信息系統審計中計算機審計的應用
消費導刊(2017年20期)2018-01-03 06:26:40
3D打印中的模型分割與打包
Fresnel衍射的計算機模擬演示
主站蜘蛛池模板: 日本精品视频| 视频国产精品丝袜第一页| 亚洲午夜国产精品无卡| 国产69精品久久| a毛片在线播放| 亚洲水蜜桃久久综合网站 | 五月婷婷综合在线视频| 男人天堂伊人网| 无码网站免费观看| 欧类av怡春院| 97久久人人超碰国产精品| 国产成人三级| 亚洲日韩高清在线亚洲专区| 丝袜久久剧情精品国产| 日本a级免费| 亚洲av无码人妻| 高清码无在线看| 韩国福利一区| 亚洲人成网站色7799在线播放| 国产区免费| 国产精品思思热在线| 亚洲日产2021三区在线| 婷婷综合亚洲| 91福利一区二区三区| 真实国产精品vr专区| 国产精品亚洲一区二区在线观看| 亚洲美女AV免费一区| 中文字幕在线观看日本| 久久91精品牛牛| 中文字幕亚洲乱码熟女1区2区| 乱人伦视频中文字幕在线| 九九精品在线观看| 高h视频在线| 亚洲中文字幕无码爆乳| 亚洲成aⅴ人在线观看| 欧美在线三级| 国内毛片视频| 国产在线一区视频| 992Tv视频国产精品| 伊人蕉久影院| 最新国产网站| 欧美天堂在线| 亚洲欧美成人在线视频| 精品国产成人国产在线| 影音先锋亚洲无码| 久久久亚洲色| 精品小视频在线观看| 亚洲福利片无码最新在线播放| 国产人成在线视频| 国产理论最新国产精品视频| 国产精品入口麻豆| 久久国产精品娇妻素人| 免费aa毛片| 日韩美毛片| 国产精品久久久久久久久kt| 黄色网站在线观看无码| 伊人久久福利中文字幕| 欧美精品v| 亚洲无码精彩视频在线观看| 亚洲福利视频一区二区| 欧美精品成人| 亚洲av无码专区久久蜜芽| 国产精品免费电影| 婷婷六月综合| 国产91成人| 无码丝袜人妻| 久久久久夜色精品波多野结衣| 亚洲视频免费在线| 国产精品福利社| 青青操国产| 亚洲黄网在线| 久久99国产综合精品女同| 日韩毛片在线播放| 午夜丁香婷婷| 青青草国产在线视频| 亚洲日韩图片专区第1页| 色爽网免费视频| 无码'专区第一页| 色综合日本| 精品久久国产综合精麻豆| 成人福利在线视频免费观看| h视频在线播放|