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

半懸臂式燃料元件在間隙限制約束下非線性振動(dòng)的等效方法研究

2023-06-05 06:46:16楊國威凡天娣陳建偉
核安全 2023年3期
關(guān)鍵詞:振動(dòng)模型

楊國威,張 勇,宋 勇,凡天娣,陳建偉,林 峰

(1.中國科學(xué)院合肥物質(zhì)科學(xué)研究院,合肥 230031;2.中國科學(xué)技術(shù)大學(xué),合肥 230026;3.中子科學(xué)國際研究院,青島 266041)

與大型核電站相比,小型模塊化反應(yīng)堆具有投資更低、應(yīng)用場景靈活、模塊化設(shè)計(jì)和固有安全性高等優(yōu)勢[1]。在各種第四代反應(yīng)堆中,鉛基堆具有中子性能良好、傳熱能力優(yōu)越、燃料增殖性能優(yōu)良和固有安全性高[2]等特點(diǎn),更具有小型化應(yīng)用前景[3-5]。

近年來,各國提出了多種鉛基堆設(shè)計(jì)方案,包括俄羅斯的SVBR-75/100[6]和BRESTOD-300項(xiàng)目[7],比利時(shí)的MYRRHA項(xiàng)目[8],歐盟的ALFRED項(xiàng)目[9],以及中國的CLEAR-I項(xiàng)目[10-12]。目前,鉛基堆中燃料元件均采用了下端固定上端簡支的結(jié)構(gòu),該結(jié)構(gòu)可以利用液態(tài)鉛合金浮力,同時(shí)還提供了燃料元件軸向熱膨脹裕量。為了提高燃料元件穩(wěn)定性,添加了中間支撐結(jié)構(gòu),其中SVBR、BREST-300與MYRRHA采用繞絲結(jié)構(gòu)提供中間支撐,ALFRED采用格架提供中間支撐,代表性鉛基堆燃料元件參數(shù)見表1。

表1 代表性鉛基堆燃料元件設(shè)計(jì)方案Table 1 Representative lead-based reactor fuel pin design

為了盡量減小反應(yīng)堆堆芯體積,小型模塊化反應(yīng)堆的燃料元件采用密集型排列。相較于傳統(tǒng)鉛基堆,小型模塊化鉛基堆的燃料元件長度大大縮短,可以省略中間支撐結(jié)構(gòu),并采用上端固定,下端間隙配合的半懸臂結(jié)構(gòu)。與傳統(tǒng)的固定方式相比,半懸臂式燃料元件的穩(wěn)定性略有降低[13],但減少了安裝空間,解決了安裝難題。

1 半懸臂式燃料元件非線性振動(dòng)問題

半懸臂式燃料元件底部間隙導(dǎo)致了燃料元件底端接觸變化。接觸是一種很普遍的非線性行為,面對接觸非線性問題,傳統(tǒng)的線性模態(tài)分析技術(shù)無法得到準(zhǔn)確的結(jié)果。由于系統(tǒng)非線性因素的控制難度較大,利用實(shí)驗(yàn)手段進(jìn)行的非線性模態(tài)研究較少,目前采用的方法大多是尋求非線性模態(tài)的近似解解析,因此離散系統(tǒng)的自由度不超過3個(gè),否則,計(jì)算量過于龐大。Rosenberg[14]等引入非線性模態(tài)理論,主要研究雙自由度離散、無阻尼非線性系統(tǒng)的自由振動(dòng)。堆芯內(nèi)存在上百根燃料元件,若在多元件抗震分析中完全采用真實(shí)的間隙單元模擬,工作量與計(jì)算量巨大,為了減少工作量與計(jì)算量,根據(jù)燃料元件整體動(dòng)態(tài)特性,將下部間隙采用等效彈簧模擬,將接觸非線性結(jié)構(gòu)等效為彈簧線性結(jié)構(gòu)。

國際上完成的分析與試驗(yàn)[15]證明,第二階及更高階頻率對整個(gè)元件在地震情況下的總的響應(yīng)和貢獻(xiàn)比例很小。故采用第一階特征頻率等效原則,通過單根燃料元件自由振動(dòng)分析,將間隙配合懸臂梁等效為彈支梁模型。底端間隙配合半懸臂梁與等效彈簧彈支梁的振動(dòng)分析模型,如圖1所示。常規(guī)的等效彈簧剛度計(jì)算過程如下:根據(jù)懸臂梁模態(tài)計(jì)算得到的特征向量施加初始平動(dòng)和轉(zhuǎn)動(dòng)位移,然后釋放,可以得到半懸臂式燃料元件第一階特征碰撞運(yùn)動(dòng)頻率,調(diào)整代替真實(shí)間隙的等效彈簧剛度,當(dāng)兩種不同模型的第一階特征碰撞運(yùn)動(dòng)頻率相等時(shí),此時(shí)所對應(yīng)的彈簧剛度即為等效彈簧剛度,此時(shí)半懸臂式燃料元件的間隙碰撞效應(yīng)與等效彈支梁彈簧效應(yīng)等效[16]。

圖1 模型線性化過程Fig.1 Model linearization process

上述方法需要進(jìn)行多次迭代,無法快速準(zhǔn)確地得到等效彈簧剛度。因此,本文推導(dǎo)得到了一種快速準(zhǔn)確地計(jì)算等效彈簧剛度的方法。該方法具體步驟如下:第一步,結(jié)合彈支梁運(yùn)動(dòng)方程,建立運(yùn)動(dòng)周期與等效剛度關(guān)系,確定等效彈簧剛度值;第二步,從單擺碰撞模型入手,逐步等效得到懸臂梁第一階碰撞運(yùn)動(dòng)周期計(jì)算公式,進(jìn)而求解彈簧等效剛度。

2 半懸臂式燃料元件等效剛度求解方法

2.1 彈支梁運(yùn)動(dòng)方程推導(dǎo)

第一步,推導(dǎo)得到彈支梁振動(dòng)方程,確定等效彈簧剛度K的計(jì)算公式。

考慮具有彈性控制邊界條件的均勻歐拉懸臂梁,如圖2所示。其中y(x,t) 是梁在時(shí)間t點(diǎn)x處的橫向位移,E是梁的楊氏模量,I表示面積的二階矩,ρ是線性質(zhì)量。

圖2 歐拉彈支梁示意圖Fig.2 Schematic diagram of Euler’s bullet support beam

忽略剪切變形和轉(zhuǎn)動(dòng)慣量效應(yīng)的均勻彈性,梁的自由彎曲振動(dòng)運(yùn)動(dòng)方程為:

彎曲梁的自由振動(dòng)方程是一個(gè)四階偏微分方程。求解微分方程,可采用分離變量法,式(1)可以表示為:

其中T(t)是一個(gè)簡諧函數(shù),可以表示為:

其中ω是系統(tǒng)的固有頻率。

式(2)可以表示為:

將式(4)帶入式(2),可以得到:

將微分方程式(5)的解設(shè)為:

微分方程中有4個(gè)參數(shù)需要確定,邊界條件表示如下。

(1)在固定端,梁的撓度和拐角均為0。固定端的邊界條件寫為式(7)、式(8)。

(2)在彈性支撐端,梁的彎矩為0,邊界條件為式(9),剪力與彈簧力平衡,邊界條件為式(10)。

通過上述4個(gè)邊界條件,解得彈支梁振動(dòng)方程:

當(dāng)彈支梁結(jié)構(gòu)確定后,E、I、l也隨之確定,根據(jù)式(11)可知,只需要確定λ值,即可以求出等效彈簧剛度。是一個(gè)與彈支梁運(yùn)動(dòng)周期相關(guān)的參數(shù),即確定半懸臂式燃料元件的一階固有頻率,便可以求解等效彈支梁彈簧的剛度。

2.2 半懸臂式燃料元件一階固有頻率計(jì)算

在半懸臂式燃料元件一階固有頻率求解過程中,采用如下簡化方法:第一步,通過假設(shè)一個(gè)簡單的單擺模型來初步估計(jì)碰撞運(yùn)動(dòng)周期;第二步,根據(jù)單擺與懸臂梁的差異,對上一步得到的計(jì)算公式進(jìn)行優(yōu)化,進(jìn)而得到帶間隙限制的懸臂梁系統(tǒng)一階固有頻率的計(jì)算公式。與傳統(tǒng)的懸臂梁一階碰撞運(yùn)動(dòng)法相比,該方法在保證計(jì)算精度的同時(shí)提高了速度。

2.2.1 單擺模型碰撞運(yùn)動(dòng)

假設(shè)梁質(zhì)量集中在底端,且梁是剛性的,則懸臂梁模型簡化為單擺模型。以單擺模型為研究對象,研究在碰撞過程中單擺小球的運(yùn)動(dòng)軌跡,單擺角度為θ0,當(dāng)左側(cè)轉(zhuǎn)角為α?xí)r發(fā)生碰撞,模型如圖3所示。

圖3 單擺模型Fig.3 Pendulum model

正常情況下,單擺運(yùn)動(dòng)質(zhì)量點(diǎn)的速度曲線與運(yùn)動(dòng)曲線如圖4所示,當(dāng)左側(cè)存在擋板時(shí),底端發(fā)生碰撞,能量不發(fā)生損失,質(zhì)量點(diǎn)速度在碰撞瞬間方向反向,速度大小不變。

圖4 單擺運(yùn)動(dòng)軌跡Fig.4 Pendulum motion trajectory

不發(fā)生碰撞時(shí),單擺質(zhì)量點(diǎn)運(yùn)動(dòng)周期可以通過式(12)求得,周期僅與單擺長度有關(guān)。

其中,l為單擺長度,g為重力加速度。

單擺運(yùn)動(dòng)角度與時(shí)間成余弦關(guān)系:

碰撞狀態(tài)下的運(yùn)動(dòng)周期與夾角α相關(guān)。若小球不在擋板左側(cè)運(yùn)動(dòng),則運(yùn)動(dòng)的周期減少量即為該部分運(yùn)動(dòng)所需要的時(shí)間,該部分用時(shí)為:

其中ω可以表示為:

此時(shí)的運(yùn)動(dòng)周期為:

當(dāng)α=0時(shí),即在中間處設(shè)置有碰撞擋板,周期為,即T/2;當(dāng)α=-θ0時(shí),相當(dāng)于沒有擋板,即為原有周期。在兩個(gè)極端角度下,根據(jù)式(15)可以準(zhǔn)確計(jì)算單擺碰撞運(yùn)動(dòng)周期,初步判斷式(15)可以計(jì)算不同碰撞角度的單擺碰撞運(yùn)動(dòng)周期。

為了進(jìn)一步驗(yàn)證式(15)的正確性,使用ANSYS中CONTA174&TARGE170單元建立點(diǎn)線接觸,構(gòu)建單擺碰撞模型。單擺模型碰撞點(diǎn)的角度α分別為2°、4°、6°、8°,使用AYSYS與理論公式計(jì)算得到單擺碰撞運(yùn)動(dòng)周期,見表2。

表2 不同碰撞角度單擺周期Table 2 Pendulum cycles for different collision angles

與通過式(15)理論計(jì)算得到的結(jié)果進(jìn)行對比,發(fā)現(xiàn)誤差小于2%,碰撞角度越大,誤差越大,因此得知:式(15)計(jì)算結(jié)果可靠。

2.2.2 鉸接剛性梁模型碰撞運(yùn)動(dòng)

為了更接近間隙配合懸臂梁系統(tǒng),引入鉸接剛性梁模型。與懸臂梁模型相比,鉸接剛性桿不會(huì)發(fā)生形變,相較于單擺模型,桿模型質(zhì)量分布均勻。可以簡化為長度為固定點(diǎn)到中心的單擺模型,如圖5所示。

圖5 鉸接剛性梁模型Fig.5 Articulated rigid beam model

該結(jié)構(gòu)在小角度時(shí),運(yùn)動(dòng)形式與單擺一致。

不發(fā)生碰撞時(shí),剛性桿末端點(diǎn)運(yùn)動(dòng)周期可以通過式(16)求得,周期僅與長度相關(guān)。

發(fā)生碰撞的運(yùn)動(dòng)周期為:

使用LINK1&MASS&TARGE170單元建立質(zhì)量均勻剛性桿碰撞模型,桿長為1 m,材料為結(jié)構(gòu)鋼,擺角為10°。

碰撞點(diǎn)的角度α分別為2°、4°、6°、8°、10°,使用AYSYS與理論公式計(jì)算得到的運(yùn)動(dòng)周期見表3。

表3 不同角度碰撞單擺周期Table 3 Collision pendulum cycles at different angles

與通過式(17)理論計(jì)算得到的結(jié)果進(jìn)行對比,誤差小于3%,驗(yàn)證了式(17)計(jì)算結(jié)果可靠。

2.2.3 懸臂梁一階自由運(yùn)動(dòng)模型

與上兩個(gè)模型相比,懸臂梁一階形變自由運(yùn)動(dòng)更為復(fù)雜,該系統(tǒng)不僅存在動(dòng)能,還存在懸臂梁的彎曲勢能。懸臂梁底端與擋板發(fā)生彈性碰撞,碰撞后不發(fā)生能量交換,碰撞發(fā)生時(shí)間極短,可以忽略該過程的懸臂梁形變量。碰撞后,懸臂梁底端速度方向相反,速度大小不發(fā)生變化,即碰撞后懸臂梁的動(dòng)能與勢能均不發(fā)生變化,碰撞前后懸臂梁端點(diǎn)運(yùn)動(dòng)規(guī)律可以與單擺等效。

在不考慮底端間隙配合時(shí),本文對懸臂梁進(jìn)行模態(tài)分析,得到第一階特征頻率所對應(yīng)的特征向量,并將特征向量成比例添加在各單元節(jié)點(diǎn)上,作為懸臂梁自由運(yùn)動(dòng)初始位移(本文取最大平動(dòng)位移特征向量為1 mm),初始位移1 mm自由振動(dòng)的情況下,元件的底端節(jié)點(diǎn)位移—時(shí)間曲線如圖6所示。可以看到其振動(dòng)基本體現(xiàn)出正弦曲線特性,由于并未設(shè)置阻尼,振動(dòng)無衰減趨勢。自由振動(dòng)的運(yùn)動(dòng)周期為0.147 s,頻率為6.80 Hz,該頻率即為假設(shè)懸臂固定時(shí)燃料元件的第一階固有頻率。

圖6 懸臂一階自由振動(dòng)底端運(yùn)動(dòng)Fig.6 Cantilever first-order free vibration bottom end movement

懸臂梁一階形變自由運(yùn)動(dòng)底端依舊為余弦運(yùn)動(dòng),其振動(dòng)頻率與懸臂梁一階模態(tài)頻率一致。

間隙配合懸臂梁結(jié)構(gòu)如圖7所示,桿長度為L,間隙為a。

圖7 單擺運(yùn)動(dòng)模型Fig.7 Pendulum motion model

初始角度θ近似為:

碰撞角度約為:

與單擺模型和鉸接剛性梁模型進(jìn)行類比,得到發(fā)生碰撞的懸臂梁運(yùn)動(dòng)周期:

其中,T為懸臂梁一階固有振動(dòng)周期。懸臂梁一階固有頻率可以由式(21)計(jì)算得到。

只需確定懸臂梁一階固有頻率,即可以求得間隙配合懸臂梁系統(tǒng)一階振動(dòng)頻率。

本文使用BEAM3&TARGE170單元建立點(diǎn)線接觸,構(gòu)建碰撞模型。間隙距離為0.1 mm、0.3 mm、0.5 mm,得到自由端運(yùn)動(dòng)軌跡,如圖8所示,碰撞點(diǎn)存在位移擾動(dòng),可以得到半懸臂式燃料元件運(yùn)動(dòng)周期。

圖8 間隙配合懸臂系統(tǒng)自由端運(yùn)動(dòng)軌跡Fig.8 Gap with cantilever system free-end motion trajectory

本文使用AYSYS與理論公式計(jì)算得到運(yùn)動(dòng)周期,見表4。

表4 不同角度碰撞單擺周期Table 4 Collision pendulum cycle at different angles

ANSYS計(jì)算過程中,發(fā)現(xiàn)碰撞時(shí)存在抖動(dòng)現(xiàn)象,且周期均比理論計(jì)算長,兩者之間誤差小于3%,可以通過上述公式快速得到間隙配合懸臂系統(tǒng)碰撞運(yùn)動(dòng)周期。

2.3 間隙配合懸臂結(jié)構(gòu)等效剛度確定

得到間隙配合懸臂梁系統(tǒng)振動(dòng)周期后,結(jié)合式(11)可以快速得到彈支梁等效剛度K。計(jì)算得到間隙距離為0.1 mm、0.3 mm、0.5 mm時(shí),間隙配合懸臂系統(tǒng)等效剛度的計(jì)算結(jié)果見表5。

本文使用BEAM3建立懸臂梁模型,并使用COMBIN40&TARGE170對底端進(jìn)行彈支固定,模型如圖9所示。

圖9 彈支梁模型Fig.9 Projectile beam model

本文通過ANSYS計(jì)算如圖9所示彈支梁模型的固有頻率,得到不同等效剛度下彈支梁一階固有頻率,見表6。

表6 不同等效剛度的彈支梁一階固有頻率Table 6 First-order natural frequencies of elastic beams with different equivalent stiffnesses

該計(jì)算結(jié)果與間隙配合懸臂梁一階自由振動(dòng)頻率對比見表7。

表7 彈支梁與半懸臂梁式燃料元件振動(dòng)頻率Table 7 Comparison of vibration frequencies of bullet support beams and semi-cantilever beam fuel pins

由表7可以看出,兩個(gè)振動(dòng)頻率基本一致,因此,間隙配合懸臂系統(tǒng)可以使用該等效剛度下的彈支梁進(jìn)行等效替換。

3 半懸臂式燃料元件抗震分析案例驗(yàn)證

為了進(jìn)一步確定本方法的正確性,本文對半懸臂式燃料元件與等效彈支梁進(jìn)行時(shí)程分析,分析模型如圖10所示。半懸臂式燃料元件底端間隙分別為0.1 mm、0.3 mm、0.5 mm,對應(yīng)等效彈支梁彈簧剛度分別為1.22×107N/m、1.41×108N/m、5.47×108N/m。

圖10 半懸臂式燃料元件與等效彈支梁分析模型Fig.10 Analysis model of semi-cantilever fuel pin and equivalent bullet support beam

本文分別給兩個(gè)模型在X方向添加RG1.60加速度譜,RG1.60加速度時(shí)程譜如圖11所示。

圖11 RG1.60加速度時(shí)程譜Fig.11 RG1.60 acceleration timer spectrum

本文分別對圖10中的兩個(gè)模型進(jìn)行25 s的時(shí)程分析,時(shí)間步長為0.01 s,其中等效彈支梁模型為線性模型,收斂速度快,僅需要17 min;半懸臂式燃料元件模型存在接觸非線性,需要用直接積分的方法進(jìn)行求解,由于所建立的動(dòng)力學(xué)方程的每個(gè)對象是最基本的單元,故求解速度和收斂時(shí)間較慢,結(jié)果文件龐大,計(jì)算耗時(shí)1258 min(20.9 h)。等效彈支梁模型計(jì)算耗時(shí)僅為半懸臂式燃料元件模型耗時(shí)的1/74。

本文分析得到兩個(gè)模型最底端位移軌跡,如圖12所示。對比兩個(gè)模型,半懸臂式燃料元件模型與等效彈支梁模型運(yùn)動(dòng)軌跡有所不同。由于底端存在碰撞,所以半懸臂式燃料元件底端位移軌跡相較于等效彈支梁模型更復(fù)雜,出現(xiàn)大量小幅振動(dòng)。若忽略這些小幅振動(dòng),則兩者位移軌跡整體趨勢具有一定相似性,且兩個(gè)模型的最大位移差距較小。

圖12 半懸臂式燃料元件模型與等效彈支梁模型底端位移軌跡Fig.12 Semi-cantilever fuel pin model and equivalent elastic support beam model bottom end displacement trajectory

半懸臂式燃料元件模型位移略大,最大位移見表8,不同間隙下兩者之間最大位移誤差均小于5%,在抗震分析中,等效彈支梁模型可以準(zhǔn)確地模擬出半懸臂式燃料元件底端最大位移情況。

表8 彈支梁與半懸臂式燃料元件與等效彈支梁底端最大位移Table 8 Max displacement at the bottom end of the ammunition support beam and semi-cantilever fuel pin and equivalent elastic support beam

本文分析得到兩個(gè)模型中間位置位移軌跡,如圖13所示。對比兩個(gè)模型,間隙為0.5 mm時(shí),半懸臂式燃料元件整體呈現(xiàn)懸臂梁陣型,兩者存在較大差距。誤差隨著間隙減小而減小,當(dāng)間隙小于0.3 mm,且兩個(gè)模型振型基本一致時(shí),振型更接近簡支梁,半懸臂式燃料元件振動(dòng)依舊更為復(fù)雜,最大位移基本相同,等效彈支梁模型可以準(zhǔn)確地模擬出半懸臂式燃料元件中間運(yùn)動(dòng)情況。

圖13 半懸臂式燃料元件模型與等效彈支梁模型中間位置位移軌跡Fig.13 Displacement trajectory in the middle of the semi-cantilever fuel element model and the equivalent elastic support beam model

結(jié)合底端位移軌跡分析,可以基本確定,彈支梁與半懸臂梁式燃料元件模型與等效彈支梁模型的運(yùn)動(dòng)軌跡基本一致,且在不同位置處最大位移基本一致。

在整個(gè)時(shí)程分析中,半懸臂式燃料元件模型與等效彈支梁模型最大應(yīng)變均集中在頂端固定處,不同參數(shù)下最大應(yīng)變見表9,對比得到兩模型之間誤差小于5%。半懸臂式燃料元件模型與等效彈支梁模型在地震波作用下,受力情況基本一致。

表9 彈支梁與半懸臂梁式燃料元件與等效彈支梁最大應(yīng)變Table 9 Maximum strain of the ammunition support beam and semi-cantilever beam fuel element and equivalent elastic support beam

綜上,等效彈支梁與半懸臂梁式燃料元件模型與等效彈支梁模型位移軌跡基本一致,且受力情況基本一致,兩者之間可以相互等效替代。

4 結(jié)論

(1) 質(zhì)量均勻分布的等截面間隙配合懸臂梁系統(tǒng)可以等效為彈支梁結(jié)構(gòu),進(jìn)而消除碰撞非線性,大大簡化計(jì)算過程。

(2) 結(jié)合彈支梁振動(dòng)方程與間隙配合懸臂梁一階振動(dòng)周期,可以快速得到彈支梁等效彈簧剛度K,等效彈支梁一階固有頻率與間隙配合懸臂梁系統(tǒng)一階振動(dòng)頻率基本一致,誤差小于3%,通過等效公式,可以快速進(jìn)行半懸臂式燃料元件在間隙限制約束下的模態(tài)分析。

(3) 通過RG1.60時(shí)程譜25 s的時(shí)程分析,確定等效彈支梁計(jì)算耗時(shí)僅為半懸臂梁式燃料元件模型的1/74,且關(guān)鍵部位的位移軌跡基本一致,受力情況基本一致,兩者之間可以相互等效替代。

猜你喜歡
振動(dòng)模型
一半模型
振動(dòng)的思考
噴水推進(jìn)高速艇尾部振動(dòng)響應(yīng)分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
振動(dòng)攪拌 震動(dòng)創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動(dòng)性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国内老司机精品视频在线播出| 久久久91人妻无码精品蜜桃HD| 亚洲一区无码在线| 国产精品尤物铁牛tv| 欧美日韩中文国产va另类| 一本视频精品中文字幕| 亚洲精品手机在线| 国产成人久视频免费| 日韩av无码DVD| 久久人搡人人玩人妻精品一| 精品视频一区在线观看| 一级毛片免费高清视频| 亚洲va在线观看| 欧美日韩一区二区三| 亚洲综合狠狠| 亚洲人网站| 伊人久久精品无码麻豆精品| 亚洲最新在线| 欧美中文字幕在线播放| 99在线观看精品视频| 99视频在线观看免费| 国产噜噜噜视频在线观看| 国产精品亚洲一区二区三区z| 国产在线精品人成导航| 久久国产精品娇妻素人| 国产高清无码麻豆精品| 欧美在线一二区| 国产区91| 精品国产一二三区| 国产v欧美v日韩v综合精品| 国产精品毛片一区| 国产精选小视频在线观看| 亚洲另类国产欧美一区二区| 亚洲精品欧美日本中文字幕| 999精品视频在线| 暴力调教一区二区三区| 女人av社区男人的天堂| 精品无码一区二区三区在线视频| 日韩欧美中文亚洲高清在线| 波多野结衣在线一区二区| 国产性精品| 欧美精品色视频| 在线精品亚洲国产| 日本一区中文字幕最新在线| 中文字幕乱码二三区免费| 国产精品亚洲一区二区三区在线观看| 中文国产成人精品久久| 久久超级碰| 欧美亚洲第一页| 久久香蕉欧美精品| 综合人妻久久一区二区精品| 欧美日韩精品一区二区视频| 在线国产91| 国产毛片基地| 99在线视频免费| 成人蜜桃网| 久久精品国产在热久久2019| 亚洲国产成人精品无码区性色| 欧美日本激情| 欧美激情网址| 亚洲视频在线青青| 亚洲人成人无码www| 九色最新网址| 日本午夜三级| 日韩成人在线一区二区| 欧美色图第一页| 亚洲无线视频| 综合网久久| 国产美女精品一区二区| 精品国产免费人成在线观看| 亚洲国产理论片在线播放| 亚洲精品在线观看91| 日韩在线第三页| 久久这里只精品国产99热8| 国产福利大秀91| 国产免费高清无需播放器| 色婷婷色丁香| 久久精品国产亚洲麻豆| 国产成人亚洲欧美激情| 色成人综合| 一级福利视频| 国产成人1024精品下载|