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

核島結(jié)構(gòu)PCS水箱FSI效應(yīng)簡(jiǎn)化方法研究

2019-01-23 10:17:58李小軍宋辰寧周國良
振動(dòng)與沖擊 2019年2期
關(guān)鍵詞:結(jié)構(gòu)質(zhì)量模型

李小軍, 宋辰寧, 周國良, 魏 超

(1. 北京工業(yè)大學(xué) 建筑工程學(xué)院,北京 100124;2.中國地震局地球物理研究所,北京 100081;3. 山東建筑大學(xué) 土木工程學(xué)院, 濟(jì)南 250101; 4. 環(huán)境保護(hù)部核與輻射安全中心,北京 100082)

AP1000是美國西屋公司設(shè)計(jì)研發(fā)的第三代先進(jìn)壓水堆核電廠,如圖1所示,核島結(jié)構(gòu)主要由屏蔽廠房,輔助廠房和鋼安全殼組成。屏蔽廠房頂部的冷卻水箱是非能動(dòng)安全殼冷卻系統(tǒng)(Passive Containment Cooling System, PCS)的重要組成部分,液體和屏蔽廠房之間的流固耦合(Fluid-Structure Interaction, FSI)效應(yīng)會(huì)對(duì)結(jié)構(gòu)的動(dòng)力特性和地震反應(yīng)造成影響。由于FSI效應(yīng)的復(fù)雜性,常用有限元軟件(ABAQUS,ANSYS,ADINA等)在求解該問題時(shí),通常需要花費(fèi)較長(zhǎng)時(shí)間。因此,在對(duì)核島結(jié)構(gòu)進(jìn)行動(dòng)力分析(地震、爆炸、撞擊等)時(shí),有必要對(duì)液體晃動(dòng)進(jìn)行一定的簡(jiǎn)化,在滿足計(jì)算要求的前提下,提高計(jì)算效率。

以地震反應(yīng)分析為例,針對(duì)儲(chǔ)液結(jié)構(gòu)的經(jīng)典理論,最早是Housner[1-2]針對(duì)規(guī)則儲(chǔ)液罐,基于液體無旋、無黏、不可壓縮及剛性罐壁、線彈性等假定提出的近似計(jì)算方法。核心理論是將液體動(dòng)力效應(yīng)分為兩部分,一部分是脈沖壓力,另一部分是對(duì)流壓力,使得流固耦合求解過程得到簡(jiǎn)化。在此基礎(chǔ)上,Haroun等[3-4]針對(duì)柔性結(jié)構(gòu),考慮儲(chǔ)液罐變形的影響,提出Haroun-Housner模型,進(jìn)一步將液體脈沖壓力分為柔性脈沖分量和剛性脈沖分量,該簡(jiǎn)化模型更適用于工程結(jié)構(gòu),被相關(guān)抗震設(shè)計(jì)規(guī)范和標(biāo)準(zhǔn)廣泛采用。

Housner模型和Haroun-Housner模型是針對(duì)規(guī)則儲(chǔ)液結(jié)構(gòu)地震分析的兩種常用的簡(jiǎn)化方法,對(duì)于不規(guī)則形狀儲(chǔ)液結(jié)構(gòu)的簡(jiǎn)化研究,大多集中在環(huán)形水箱[5]和純錐形水箱[6],核島結(jié)構(gòu)中涉及的PCS水箱為底部錐形的環(huán)形圓柱結(jié)構(gòu),對(duì)于該特殊的儲(chǔ)液結(jié)構(gòu)并沒有給出相關(guān)的簡(jiǎn)化計(jì)算方法。本文基于現(xiàn)有簡(jiǎn)化方法,針對(duì)PCS水箱,提出一種適用于核島結(jié)構(gòu)三向地震反應(yīng)分析的簡(jiǎn)化模型,并與FSI模型計(jì)算結(jié)果進(jìn)行對(duì)比,驗(yàn)證簡(jiǎn)化模型的合理性。

圖1 AP1000核電廠示意圖Fig.1 Diagram of AP1000 nuclear power plant

1 理論分析

考慮核島主體結(jié)構(gòu)和PCS水箱的幾何尺寸和結(jié)構(gòu)特性,可將其視為剛性儲(chǔ)液結(jié)構(gòu),滿足Housner模型的基本適用條件。因此,本文基于Housner模型的基本假定,開展簡(jiǎn)化計(jì)算模型研究,以提出工程分析計(jì)算精度的模擬PCS水箱液固耦合作用的簡(jiǎn)化模型。

1.1 Housner模型

1957年,Housner基于如下假定提出簡(jiǎn)化模型:①儲(chǔ)液結(jié)構(gòu)截面規(guī)則;②底部為平底;③適用于水平地震作用;④剛性結(jié)構(gòu)。其核心思想是將剛性儲(chǔ)液結(jié)構(gòu)中液體晃動(dòng)等效為“質(zhì)量-彈簧”模型,如圖2所示。圖2中:m0為脈沖質(zhì)量;m1為一階對(duì)流質(zhì)量;mn為n階對(duì)流質(zhì)量;k1為一階等效彈簧剛度;kn為n階等效彈簧剛度;h0為脈沖質(zhì)量高度;h1為一階對(duì)流質(zhì)量高度;hn為n階對(duì)流質(zhì)量高度。脈沖質(zhì)量附著在結(jié)構(gòu)內(nèi)壁,隨結(jié)構(gòu)作同步運(yùn)動(dòng),對(duì)流質(zhì)量依靠彈簧與箱壁相連,等效模擬液體晃動(dòng)效應(yīng)。

在具體計(jì)算和應(yīng)用時(shí),首先要確定脈沖質(zhì)量和對(duì)流質(zhì)量的比例關(guān)系以及其質(zhì)心相對(duì)于底面的高度,然后通過計(jì)算得到液體的晃動(dòng)特性以及液體晃動(dòng)產(chǎn)生的基底剪力和儲(chǔ)液結(jié)構(gòu)箱壁的傾覆彎矩。該模型簡(jiǎn)潔實(shí)用,至今仍為國內(nèi)外學(xué)者采用。

圖2 Housner模型Fig.2 Housner model

1.2 簡(jiǎn)化模型

一般的儲(chǔ)液結(jié)構(gòu)(如LNG儲(chǔ)罐、石油儲(chǔ)罐等)底部為平底,在采用Housner模型進(jìn)行簡(jiǎn)化計(jì)算時(shí)基于基本假定,通常只考慮水平地震作用。核島PCS水箱為底部錐形的環(huán)形圓柱結(jié)構(gòu),且相關(guān)規(guī)范條例要求必須進(jìn)行三向地震反應(yīng)分析,而Housner模型中并沒考慮液體豎向振動(dòng)。

考慮PCS水箱的特殊結(jié)構(gòu)形式和需要進(jìn)行三向地震反應(yīng)分析要求,本文參考Housner模型的計(jì)算方法,分兩步推導(dǎo)建立PCS水箱簡(jiǎn)化模型:①計(jì)算脈沖質(zhì)量和對(duì)流質(zhì)量的比例關(guān)系、對(duì)流質(zhì)量質(zhì)心相對(duì)于底面的高度;②計(jì)算液體的晃動(dòng)頻率,根據(jù)晃動(dòng)頻率和阻尼比確定連接對(duì)流質(zhì)量和箱壁的彈簧阻尼器參數(shù)。

基于Housner模型可以得出:液體高階對(duì)流壓力遠(yuǎn)小于一階對(duì)流壓力。因此,為了簡(jiǎn)化計(jì)算,本文模型中只考慮液體的一階晃動(dòng),即將PCS水箱中液體的總質(zhì)量分為脈沖質(zhì)量和一階對(duì)流質(zhì)量。設(shè)PCS水箱液體體積為V,液體總質(zhì)量為m,液面設(shè)計(jì)高度為h,內(nèi)半徑為Ri,外半徑為Ro,參考已有研究成果[7-8],將非規(guī)則環(huán)形水箱等效成圓柱形水箱,等效半徑R、一階對(duì)流質(zhì)量m1、一階對(duì)流質(zhì)量高度h1計(jì)算公式為

(1)

(2)

(3)

針對(duì)環(huán)形圓柱水箱,Meserole等給出了晃動(dòng)力學(xué)模型的參數(shù)表達(dá)式,考慮到PCS水箱為底部錐形的非規(guī)則環(huán)形結(jié)構(gòu),本文首先將其等效為平底環(huán)形結(jié)構(gòu)。如圖3所示,基本思路是保證V,Ri和Ro不變,將PCS水箱等效成環(huán)形圓柱水箱,計(jì)算等效液面高度ha,求解晃動(dòng)頻率,計(jì)算公式為

(4)

λ=Ri/Ro

(5)

ki=(ha/Ro)ξi

(6)

(a)PCS水箱

(b)環(huán)形圓柱水箱圖3 等效模型示意圖Fig.3 Diagram of equivalent model

在PCS水箱簡(jiǎn)化模型中,將脈沖質(zhì)量m0離散成質(zhì)量點(diǎn)附加于水箱底部和水箱壁,參考ASCE4-98規(guī)范[9],認(rèn)為液體豎向振動(dòng)近似于剛性,將一階對(duì)流質(zhì)量m1分成水平分量mh和垂直分量mv,將mv離散成質(zhì)量點(diǎn)附加于水箱底部。為了避免對(duì)流質(zhì)量在動(dòng)力計(jì)算過程中引起水箱壁局部應(yīng)力集中,對(duì)mh進(jìn)行分層離散處理,通過彈簧阻尼器連接于水箱壁,如圖4所示?;诮Y(jié)構(gòu)動(dòng)力學(xué)基本公式,彈簧阻尼器的剛度和阻尼系數(shù)計(jì)算公式為

(7)

圖4 液體晃動(dòng)效應(yīng)水箱簡(jiǎn)化模型Fig.4 Simplified sloshing model of water tank

式中:k和c為單個(gè)彈簧阻尼器的剛度和阻尼;K和C為彈簧阻尼體系的整體剛度和阻尼;n為彈簧阻尼器數(shù)目;ω1為一階晃動(dòng)頻率;γ為晃動(dòng)阻尼比。

2 結(jié)構(gòu)反應(yīng)數(shù)值分析

2.1 計(jì)算模型

采用ADINA軟件建立核島結(jié)構(gòu)有限元模型,如圖5(a)所示,X軸正方向?yàn)楸?,Y軸正方向?yàn)槲鳎捎?D-Solid單元模擬核島結(jié)構(gòu)混凝土底板,3D-Shell單元模擬核島結(jié)構(gòu)墻體和樓板,3D-Beam單元模擬空間鋼桁架。在流固耦合效應(yīng)模擬方面,如圖5(b)和圖5(c)所示,模型1中采用流體單元模擬PCS水箱中的液體,模型2中采用簡(jiǎn)化模型模擬PCS水箱的FSI效應(yīng),模型材料參數(shù)見表1。

表1 數(shù)值模型材料參數(shù)Tab.1 Material parameters of numerical models

圖5 核島結(jié)構(gòu)計(jì)算模型Fig.5 Numerical models of nuclear island building

2.2 結(jié)構(gòu)模態(tài)分析

振動(dòng)模態(tài)和頻率是結(jié)構(gòu)的固有特性,采用Lanczos法提取結(jié)構(gòu)的振動(dòng)模態(tài)和頻率,結(jié)構(gòu)前兩階自振頻率及液體一階晃動(dòng)頻率見表2。前兩階振動(dòng)模態(tài)為核島結(jié)構(gòu)整體沿兩個(gè)水平方向振動(dòng),兩個(gè)模型前兩階頻率相對(duì)誤差均小于2%,液體一階晃動(dòng)頻率相對(duì)誤差小于1%??梢钥闯觯疚腜CS水箱簡(jiǎn)化模型可以較好地反映結(jié)構(gòu)的動(dòng)力特性。

表2 結(jié)構(gòu)前兩階自振頻率及液體一階晃動(dòng)頻率Tab.2 The first two natural frequencies and 1st sloshing frequency

2.3 結(jié)構(gòu)反應(yīng)時(shí)程分析

選用典型的強(qiáng)震動(dòng)記錄El Centro地震動(dòng)時(shí)程和基于核電廠標(biāo)準(zhǔn)設(shè)計(jì)反應(yīng)譜擬合的人工地震動(dòng)時(shí)程分別進(jìn)行結(jié)構(gòu)反應(yīng)的時(shí)程分析,研究核島結(jié)構(gòu)地震反應(yīng)。地震動(dòng)采用三向輸入,按照SSE(Safe Shutdown Earthquake)要求,峰值加速度均取0.3g[10-12],El Centro地震動(dòng)和人工地震動(dòng)三向加速度時(shí)程曲線和反應(yīng)譜,如圖6和圖7所示。

圖6 El-Centro地震動(dòng)加速度時(shí)程和反應(yīng)譜Fig.6 Acceleration time histories and response spectra of El Centro ground motion

圖7 人工地震動(dòng)加速度時(shí)程和反應(yīng)譜Fig.7 Acceleration time histories and response spectra of artificial ground motion

考慮篇幅所限,本文選取屏蔽廠房西側(cè)水箱頂部、水箱根部、支撐斜坡與筒體交界處3個(gè)觀測(cè)點(diǎn)和輔助廠房頂部西南角1個(gè)觀測(cè)點(diǎn)進(jìn)行反應(yīng)的具體分析,分別驗(yàn)證兩個(gè)模型的地震反應(yīng)在時(shí)域、頻域以及液壓效應(yīng)方面的吻合情況。時(shí)域方面主要以觀測(cè)點(diǎn)的峰值加速度作為對(duì)比參數(shù),頻域方面主要采用觀測(cè)點(diǎn)的樓層反應(yīng)譜進(jìn)行對(duì)比分析,液壓效應(yīng)方面主要通過對(duì)比觀測(cè)點(diǎn)處的有效應(yīng)力時(shí)程曲線體現(xiàn)。4個(gè)觀測(cè)點(diǎn)的具體位置見圖5(a)。

2.3.1 結(jié)構(gòu)反應(yīng)的峰值加速度

計(jì)算結(jié)構(gòu)反應(yīng)并提取兩個(gè)模型4個(gè)觀測(cè)點(diǎn)在兩組地震動(dòng)下的峰值加速度。分析兩個(gè)模型計(jì)算結(jié)果的相對(duì)誤差,見表3,相對(duì)誤差的計(jì)算公式為

(8)

表3中數(shù)據(jù)表明,兩個(gè)模型各觀測(cè)點(diǎn)反應(yīng)的峰值加速度相對(duì)誤差基本控制在5%以內(nèi),少數(shù)超過5%。該誤差是由于簡(jiǎn)化模型本身采用了一些假定和簡(jiǎn)化處理,不能完全如實(shí)反應(yīng)真實(shí)液體晃動(dòng)情況造成的,但誤差在可接受的范圍內(nèi),可以認(rèn)為計(jì)算結(jié)果吻合良好。

表3 結(jié)構(gòu)反應(yīng)的峰值加速度相對(duì)誤差Tab.3 Relative errors of peak accelerations of structural responses

2.3.2 結(jié)構(gòu)反應(yīng)的樓層反應(yīng)譜

根據(jù)加速度反應(yīng)時(shí)程曲線計(jì)算得到樓層反應(yīng)譜,圖8列出了兩個(gè)模型在El Centro地震動(dòng)下4個(gè)觀測(cè)點(diǎn)的樓層反應(yīng)譜對(duì)比情況。在水平方向(X向、Y向),模型2的樓層反應(yīng)譜曲線與模型1基本重合,在各個(gè)頻段差異很小。在垂直方向(Z向),P1和P2兩個(gè)觀測(cè)點(diǎn)在結(jié)構(gòu)基頻對(duì)應(yīng)的周期附近(0.10~0.30 s),兩條曲線有一定的差異,在其他頻段差異很小。兩個(gè)模型的計(jì)算結(jié)果總體吻合良好。

圖9列出了兩個(gè)模型4個(gè)觀測(cè)點(diǎn)在人工地震動(dòng)下的樓層反應(yīng)譜對(duì)比情況。與El Centro地震動(dòng)計(jì)算結(jié)果相似,P1和P2兩個(gè)觀測(cè)點(diǎn)的垂直方向樓層反應(yīng)譜在某些頻段存在一定差異,兩個(gè)模型計(jì)算結(jié)果總體吻合良好。

通過兩個(gè)模型計(jì)算結(jié)果的對(duì)比分析可以進(jìn)一步看出,對(duì)于兩組特性(頻譜和持時(shí))差異很大的地震動(dòng),觀測(cè)點(diǎn)樓層反應(yīng)譜曲線的分析結(jié)果是一致的,說明提出的水箱簡(jiǎn)化模型在模擬FSI效應(yīng)方面具有較好的適用性。針對(duì)某些頻段(0.10~0.30 s)計(jì)算結(jié)果,為了更清楚地反映其差異大小,定量地分析兩者的誤差,提取加速度反應(yīng)譜曲線在0.10 s,0.20 s,0.30 s時(shí)刻的數(shù)值和反應(yīng)譜最大值,計(jì)算其相對(duì)誤差,統(tǒng)計(jì)結(jié)果見表4和表5。相對(duì)誤差采用式(8)計(jì)算。

表4 El Centro地震動(dòng)下加速度反應(yīng)譜值相對(duì)誤差Tab.4 Relative errors of spectral accelerations under El Centro ground motion

圖8 El Centro地震動(dòng)下觀測(cè)點(diǎn)反應(yīng)譜對(duì)比Fig.8 Comparison of response spectra under El Centro ground motion

圖9 人工地震動(dòng)下觀測(cè)點(diǎn)反應(yīng)譜對(duì)比Fig.9 Comparison of response spectra under artificial ground motion

表5 人工地震動(dòng)下加速度反應(yīng)譜值相對(duì)誤差Tab.5 Relative errors of spectral accelerations under artificial ground motion

從表4和表5中可以看出,兩個(gè)模型的樓層反應(yīng)譜相對(duì)誤差基本控制在10%以內(nèi),少數(shù)達(dá)到了12%左右。相對(duì)誤差較大位置出現(xiàn)在水箱結(jié)構(gòu)上的觀測(cè)點(diǎn)(P1和P2),特別是在垂直方向,屏蔽廠房和輔助廠房上的觀測(cè)點(diǎn)(P3和P4)相對(duì)誤差較小。

PCS水箱簡(jiǎn)化模型中,在考慮水平地震作用時(shí),基于Housner模型,將一階對(duì)流質(zhì)量m1的水平分量mh通過彈簧阻尼器連接于水箱壁來模擬水平晃動(dòng);同時(shí)在考慮豎向地震作用時(shí),參考相關(guān)規(guī)范采用了液體豎向振動(dòng)近似于剛性假定,將一階對(duì)流質(zhì)量m1的垂直分量mv附加在水箱底部。這些假定與真實(shí)豎向振動(dòng)并不完全相符,加上PCS水箱底部?jī)A斜,水箱外壁特別是根部(P2)受力情況復(fù)雜,簡(jiǎn)化模型并不能完全模擬真實(shí)情況。因此簡(jiǎn)化模型的計(jì)算結(jié)果在水平方向上更接近于流固耦合模型,垂直方向的誤差相對(duì)較大。

通過圖8、圖9、表4和表5中的數(shù)據(jù)可以看出,水箱簡(jiǎn)化模型可以用于模擬考慮FSI效應(yīng)的核島結(jié)構(gòu)動(dòng)力分析,特別是水箱以下的主體結(jié)構(gòu),誤差很小;對(duì)于水箱結(jié)構(gòu)本身,采用簡(jiǎn)化模型存在一些差異。針對(duì)結(jié)構(gòu)在三向地震動(dòng)作用下的動(dòng)力反應(yīng),水箱簡(jiǎn)化模型可以在保持計(jì)算結(jié)果合理的前提下,節(jié)省計(jì)算時(shí)間,提高計(jì)算效率。

2.3.3 結(jié)構(gòu)反應(yīng)的有效應(yīng)力

由于簡(jiǎn)化模型無法輸出液體動(dòng)水壓力,液體晃動(dòng)效應(yīng)作用在結(jié)構(gòu)上主要體現(xiàn)在結(jié)構(gòu)應(yīng)力、應(yīng)變的變化,因此,為了驗(yàn)證簡(jiǎn)化模型在模擬液體晃動(dòng)和液壓效應(yīng)方面的合理性,提取4個(gè)觀測(cè)點(diǎn)處的有效應(yīng)力時(shí)程進(jìn)行對(duì)比,如圖10、圖11所示??梢钥闯?,在兩組地震動(dòng)作用下,P1觀測(cè)點(diǎn)處兩個(gè)模型的有效應(yīng)力時(shí)程曲線形狀基本相同,具體數(shù)值上有一定的差異,P2,P3,P4觀測(cè)點(diǎn)處時(shí)程曲線基本重合。為了更清楚地反映差異大小,提取有效應(yīng)力最大值進(jìn)行比較,采用式(8)計(jì)算相對(duì)誤差,如表6所示。P1觀測(cè)點(diǎn)處有效應(yīng)力最大值相對(duì)誤差分別為8.30%和5.71%,P2,P3,P4觀測(cè)點(diǎn)處相對(duì)誤差均小于5.00%,4個(gè)觀測(cè)點(diǎn)處相對(duì)誤差控制在10%以內(nèi),認(rèn)為兩個(gè)模型計(jì)算結(jié)果總體吻合良好。

圖10 El Centro地震動(dòng)下觀測(cè)點(diǎn)有效應(yīng)力對(duì)比Fig.10 Comparison of effective stress under El Centro ground motion

圖11 人工地震動(dòng)下觀測(cè)點(diǎn)有效應(yīng)力對(duì)比Fig.11 Comparison of effective stress under artificial ground motion

表6 觀測(cè)點(diǎn)有效應(yīng)力最大值相對(duì)誤差Tab.6 Relative errors of maximum effective stress

由于P1觀測(cè)點(diǎn)位于水箱頂部,采用流固耦合模型計(jì)算時(shí),能考慮到水箱頂部受液體晃動(dòng)沖擊的影響,而采用簡(jiǎn)化模型計(jì)算時(shí),液體晃動(dòng)效應(yīng)通過“質(zhì)量-彈簧-阻尼”系統(tǒng)模擬,由于對(duì)流質(zhì)量高度有限,不能很好地考慮液體對(duì)水箱頂部的晃動(dòng)沖擊,存在一定的差異,造成有效應(yīng)力的計(jì)算結(jié)果小于流固耦合模型。對(duì)于水箱下部和結(jié)構(gòu)主體(P2,P3,P4觀測(cè)點(diǎn)),簡(jiǎn)化模型能夠很好地模擬結(jié)構(gòu)反應(yīng)的有效應(yīng)力。

3 結(jié) 論

本文考慮核島結(jié)構(gòu)PCS水箱FSI效應(yīng),基于Housner模型提出一種適用于核島結(jié)構(gòu)三向地震反應(yīng)分析的PCS水箱簡(jiǎn)化模型。在簡(jiǎn)化模型中,將液體分為脈沖質(zhì)量和一階對(duì)流質(zhì)量,一階對(duì)流質(zhì)量的水平分量通過彈簧阻尼器連接于水箱壁,脈沖質(zhì)量和一階對(duì)流質(zhì)量的垂直分量附加于水箱壁和底部。在此基礎(chǔ)上,采用ADINA軟件進(jìn)行了高置冷卻水箱核島結(jié)構(gòu)三向地震反應(yīng)分析,對(duì)比了FSI模型和簡(jiǎn)化模型的計(jì)算結(jié)果,提取觀測(cè)點(diǎn)結(jié)構(gòu)反應(yīng)的峰值加速度、樓層反應(yīng)譜和有效應(yīng)力,分析其相對(duì)誤差。計(jì)算結(jié)果對(duì)比顯示,PCS水箱簡(jiǎn)化模型可用于高置冷卻水箱核島結(jié)構(gòu)的三向地震反應(yīng)分析,能夠很好地模擬FSI效應(yīng),在保證計(jì)算精度的情況下提高計(jì)算效率。

猜你喜歡
結(jié)構(gòu)質(zhì)量模型
一半模型
“質(zhì)量”知識(shí)鞏固
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
質(zhì)量守恒定律考什么
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
做夢(mèng)導(dǎo)致睡眠質(zhì)量差嗎
論《日出》的結(jié)構(gòu)
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲性一区| 亚洲美女操| 国产女人18水真多毛片18精品| 成人a免费α片在线视频网站| 青青草一区| 婷婷伊人五月| 国产一区亚洲一区| 人妻无码中文字幕第一区| 亚洲国产精品一区二区第一页免 | 色综合五月婷婷| 日韩少妇激情一区二区| 国产欧美日韩视频一区二区三区| 久久狠狠色噜噜狠狠狠狠97视色 | 国产精品白浆无码流出在线看| 欧美成人在线免费| 亚洲日产2021三区在线| 国产精品所毛片视频| 亚洲精品视频网| 呦视频在线一区二区三区| 手机在线国产精品| 波多野结衣一区二区三视频| 国产成人综合亚洲欧洲色就色| 国产精品成人第一区| 国产成人精品2021欧美日韩| 久久国产精品麻豆系列| 午夜老司机永久免费看片| 有专无码视频| 久久国产黑丝袜视频| 亚洲欧美人成人让影院| 久久96热在精品国产高清| 国产香蕉在线| 99在线视频免费| 亚洲午夜国产精品无卡| 伊人久久大香线蕉aⅴ色| av在线人妻熟妇| 97人人做人人爽香蕉精品| 91精品国产麻豆国产自产在线| 狠狠做深爱婷婷综合一区| 妇女自拍偷自拍亚洲精品| 日本国产一区在线观看| 久久久久人妻一区精品| 日韩无码黄色网站| 中文字幕伦视频| 色偷偷一区| 精品国产aⅴ一区二区三区 | 一级毛片在线免费看| 无套av在线| 亚洲国产无码有码| 国内精品视频| 欧美日韩资源| 国产精品无码AV片在线观看播放| 免费在线国产一区二区三区精品| 国产一区二区三区免费观看| 亚洲视频欧美不卡| 在线欧美日韩国产| 国产精品午夜福利麻豆| 老司机午夜精品网站在线观看| 国产91精品久久| 久久九九热视频| 精品色综合| 2022精品国偷自产免费观看| 国产大片喷水在线在线视频 | 欧美a级在线| 国产91视频免费观看| 三区在线视频| 精品久久高清| 日韩欧美网址| 国产白浆在线观看| a色毛片免费视频| 国产美女久久久久不卡| 国产黑丝视频在线观看| 国产乱人伦精品一区二区| 中国美女**毛片录像在线| 欧美精品三级在线| 亚洲天堂网2014| 久青草国产高清在线视频| 国产精品亚洲一区二区三区z| 精品成人免费自拍视频| 中文字幕天无码久久精品视频免费 | 免费国产黄线在线观看| 无码内射在线| 国产精品网址你懂的|