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

垂索-輔助索系統的建模與固有特性

2023-11-22 09:12:06孫測世焦德望趙碧航
工程力學 2023年11期
關鍵詞:模態振動模型

孫測世,焦德望,趙碧航,譚 超

(重慶交通大學土木工程學院,重慶 400074)

斜拉索是斜拉橋的主要承重結構,具有柔度大、阻尼低等特點,極易發生大幅振動,對橋梁的安全造成重大影響[1-2]。為了抑制斜拉索可能產生的大幅振動,工程上常采用氣動措施[3]、安裝阻尼器[4-7]或輔助索等措施。但隨著斜拉索的長度不斷增長,結構非線性和垂度效應的影響逐漸增大,阻尼器因受安裝位置的限制達不到理想的減振效果[8]。近年來,輔助索被認為是解決這一難題的有效手段[9],為此研究人員對拉索-輔助索系統進行了大量的理論研究和實驗研究。

早期研究中,為計算簡便,將拉索簡化為緊繃的弦。CARACOGLIA 等[10-11]建立了由輔助索連接水平弦的模型,基于解析和數值計算方法研究了系統的面內自由振動特性,并將該方法應用于一座既有斜拉橋的減振研究,證明了該方法的有效性。AHMAD 等[12-14]建立了水平弦-輔助索模型,理論分析了模型中各參數對系統面內剛度及阻尼特性的影響。GIACCU 等[15-17]采用基于能量的“等效線性化”方法研究了輔助索連接水平弦系統的非線性自由振動特性,并將模型應用到了既有橋梁的減振研究中。HE 等[18]研究了沿水平弦均勻布置彈性輔助索來抑制弦的振動,研究表明輔助索的剛度和間隙率顯著影響系統的模態頻率和振型。CHEN 等[19]建立了多個輔助索連接水平弦系統模型,研究了不同參數下彎曲剛度對系統自由振動特性的影響。值得一提的是,以孫利民和周??榻艹龃淼膰鴥葘W者在輔助索減振、輔助索加阻尼器復合減振的力學建模、理論分析和試驗研究中開展了一系列的系統研究[20-23],如:SUN 等[20]建立了由彈簧并聯阻尼器組成的輔助索連接兩根水平弦系統模型,并通過復模態分析得到系統特征方程。ZHOU 等[21]建立了采用負剛度阻尼器連接兩根水平弦的模型,在近似確定臨界粘性阻尼和負剛度的情況下,給出了附加模態阻尼比的漸近解,研究表明輔助索能顯著降低兩根水平弦的振動效應。除此以外,具有超彈性和高阻尼的形狀記憶合金(SMA)輔助索對弦-輔助索系統的振動特性影響也受到了廣泛關注[24-26]。

隨著拉索長度不斷增大,垂度效應不容忽視,若干最新研究均考慮了垂度影響。2019 年起,SUN 等[27]首先建立了考慮垂度的拉索-輔助索模型,研究表明垂度幾乎影響系統的所有振動模式,只有特定模式的頻率與垂度無關。其后,CHEN 等[28]提出了一種基于分量模態綜合法的索網自由振動和強迫振動分析的通用數值模型,并與SUN 的模型進行了對比分析。AHMAD[29]建立了雙懸索-輔助索系統模型,討論了考慮垂度和忽略垂度情況下不同參數對垂索-輔助索系統一階頻率的影響,研究表明忽略垂度影響會造成非常大的誤差。DI等[30-31]進一步考慮了含有預張力的輔助索連接兩根小垂度拉索的模型,研究了輔助索預張力對垂索-輔助索系統動力特性影響。

綜上所述,輔助索對斜拉索減振有著廣闊的應用前景,已受到國內外學者的廣泛關注。目前,垂索-輔助索模型多是在有量綱或部分無量綱化前提下進行的推導,或者推導過程并非直接從無量綱運動方程出發,不便更好地把握系統關鍵參數及其物理意義。另外,對于考慮垂度的拉索-輔助索系統,關鍵參數對系統固有特性影響的研究尚需進一步豐富。為此,本文推導了由N根垂索和M道輔助索組成的索網系統的無量綱運動方程,建立了垂索-輔助索動力學模型;然后,退化到雙索-輔助索系統,并求得其無量綱頻率方程和模態方程;最后,研究了系統關鍵參數對固有特性的影響,以期為后續輔助索減振設計提供理論參考。

1 動力學模型

1.1 垂索-輔助索動力學模型的建立

圖1 表示由輔助索(忽略質量)連接小垂度拉索[32]構成的一般垂索-輔助索系統,垂索總數為N,輔助索總道數為M,共計M(N+1)段。每根拉索的單位質量、長度、水平張力、軸向剛度分別為。每根輔助索剛度為,其中j=1, 2, ···,N+1,表示 第j根 拉 索,p=1, 2, ···,M+1,表示第p索段。整個垂索-輔助索系統中的拉索可分成N(M+1)個索段,各索段的長度為,其中j=1, 2,···,N,p=1, 2, ···,M+1。

圖1 垂索-輔助索系統Fig.1 The sagged-cable-crosstie system

如圖1 所示在各拉索j建立局部直角坐標系,xj從拉索左端開始,yj(xj)表示拉索j的初始位移。在自重作用下拉索的初始構型為:

式中,g*為重力加速度。為表示各索段的動態位移,以各索段左端為原點建立了如圖1 所示的索段局部直角坐標系。表示第j根索第p段的動態位移,其中j=1, 2, ···,N+1;p=1, 2, ···,M+1。每個索段考慮垂度的運動方程均可采用Irvine 的拉索經典理論方程表示[32],即:

為使運動方程式(4)更具一般性,引入無量綱量:

式中,L*為使垂索-輔助索系統在統一空間尺度下進行無量綱化而引入的任意長度。當L*=L*j時,即表示其他索的空間無量綱化均按第j根的索長進行。

利用式(5)對式(4)進行無量綱化處理,可得無量綱運動方程:

式中:

顯然,式(6)是考慮垂度的無量綱波動方程,其中:參數εj,p反映輔助索的安裝位置;為索j的Irvine 參數;1/αj是索j的無量綱波速。αj是對式(4)進行無量綱化處理的過程中產生的參數,本質上是將系統在統一時間尺度下進行無量綱化而引入的綜合參數,其中包含了整個垂索-輔助索系統的索力和與質量和之比(對應系統的波速),又包含了索j的索力與質量比(對應索j的波速)。由于式(6)為徹底無量綱化后動力學方程,具有普適性,物理世界中圖1 所示的所有索網結構均由該方程描述。由式(6)還可知,系統運動方程僅依賴于εj,p、和αj三個參數,這是后續開展參數分析的重要基礎,而若不進行無量綱化處理很難直接掌握這三個反映系統動力學本質的基礎參數。

1.2 垂索-輔助索動力學模型的通解

采用分離變量法將運動方程式(6)中的空間和時間進行分離,令:

將式(7)代入式(6)并求解可得任意索段的模態函數:

式 中:βj=αjω為 索j橫 向 振 動 的 無 量 綱 波 數;Bj,p和Dj,p為各索段模態函數的未知常數。

對式(8)兩邊從0~εj,p積分,并將各段索相加后再回代入式(8)可得:

式中,Qj為垂度產生的附加項,其表達式為:

是一個有量綱波數,而本文中的βj為無量綱波數,是一個反映了索j和整個結構之間關系的綜合參數。

1.3 垂索-輔助索動力學模型的定解

為了得到模態頻率和模態振型,依次引入拉索的邊界條件、位移連續條件和輔助索兩端節點處力的平衡條件。

將式(10)代入式(9)可得:

聯立整個索網N根拉索的邊界條件代入式(11)共可得到2N個方程。

由輔助索kj,p和拉索j的節點處力的平衡條件(假設輔助索拉力為正)可得:

將式(9)和式(11)代入式(12)可得:

聯立整個索網N×M個節點處的力平衡條件代入式(13)共可得到N×M個方程。

在輔助索kj,p和拉索j的節點左右需滿足位移連續條件,則有:

將式(9)代入式(14),得:

聯立整個索網N×M個節點處的位移連續條件,代入式(15)共可得到N×M個方程。

綜上,由式(11)、式(13)和式(15)聯立可以得到2N(M+1)個方程,整個索網中各索段模態函數的未知數Bj,p和Dj,p總數也為2N(M+1)個,因此,將方程聯立并寫成矩陣形式可得:

式中:S= {B1,1,B1,2, ···,B1,p, ···,B1,M+1,D1,1,D1,2, ···,D1,p, ···,D1,M+1, ···,Bj,1,Bj,2, ···,Bj,p, ···,Bj,M+1,Dj,1,Dj,2, ···,Dj,p, ···,Dj,M+1, ···,BN,1,BN,2, ···,BN,p, ···,BN,M+1,DN,1,DN,2, ···,DN,p, ···,DN,M+1}T為一個包含整個索網各索段模態函數的未知數Bj,p和Dj,p的向量;0 為零向量;R為2N(M+1)階的系數矩陣。顯然S不等于0,因此令系數矩陣R對應的行列式等于0,即可求得N根垂索和M道輔助索組成的索網系統的頻率方程,進而可以求得各階頻率和模態。

2 雙索-輔助索系統退化模型

實際應用時,一根輔助索可能僅連接少量拉索,如近年對蘇通大橋開展的輔助索減振實橋研究中便是連接3 根斜拉索[28]。因此,為便于討論,本文暫考慮最簡單的情況,將垂索-輔助索一般模型退化到一根輔助索連接兩根拉索的結構,即N=2、M=1。此時,方程中j=1, 2;p=1, 2。另外,令L*=L*1。此時式(16)中的S={B11,B12,D11,D12,B21,B22,D21,D22}T為一個包含8 個未知常數的向量;R為一個8 階的系數矩陣(見附錄)。

系統的頻率可以通過令系數矩陣R的行列式等于0 求得,展開行列式進行三角函數簡化,并合并同類項,可得系統的頻率方程:

其中:

式(17)中第一項為不考慮垂度的項,其余三項為考慮垂度影響的項,當γj=0 時便可以得到忽略垂度影響的頻率方程:

式(18)進行一定量綱處理可還原到AHMAD[12]的頻率方程;同時,若考慮兩所參數相同,則式(18)可進一步化簡為:

式(19)經過一定量綱處理亦與ZHOU[22]推導頻率方程忽略阻尼影響后的情形一致。

當式(17)中的k2,1取無窮大時方程可退化為:

式(20)亦可通過一定量綱處理還原到SUN 等[27]研究的剛性輔助索的情形。

3 模型驗證

3.1 雙索單輔助索退化模型驗證

為驗證方程推導的正確性,參照AHMAD[29]論文中拉索的參數進行計算,并與其計算結果對比。兩根拉索參數為:

輔助索放置在拉索的1/4 截面處,其剛度取無窮大,將上述參數代入式(17),分別求得系統考慮垂度(=2)和忽略垂度(=0)情況下前十階無量綱頻率,再將頻率轉化為有量綱頻率與AHMAD[12,29]論文進行對比;另外,利用上述拉索參數采用有限元軟件ANSYS 算得的前十階頻率結果如表1 所示。由表可知:本文求得的頻率有量綱化后與AHMAD 的解相同,同時與ANSYS結果也非常接近,說明方程推導的正確性。

表1 模態頻率對比Table 1 Comparison of frequencies

從表1 可以看到,各階模態振動均對應兩個頻率,這是由于各階模態振動均含有兩個模態振型,即同相振型和反相振型,較小的頻率對應同相振動模態振型,較大的頻率對應反相振動模態振型[24]。另外,還可以發現,當=0 時,第三階反相、第四階同相、第四階反相振動頻率相同,這與輔助索剛度和位置有關。

3.2 五索雙輔助索復雜模型驗證

為進一步驗證本文推導方程對更復雜索網的正確性,參照AHMAD[33]拉索參數(表2)建立了由2 道輔助索連接5 根拉索的模型。2 道輔助索布置在第一根拉索的L1/3 處,輔助索無量綱剛度取100。求得忽略垂度(λ=0)情況下前10 階無量綱頻率,再將頻率轉化為有量綱頻率與AHMAD[33]論文進行對比;同時,建立ANSYS 模型進行對比,結果如表3 所示。可見,本文所求頻率有量綱化后與AHMAD 論文[33]及ANSYS 結果均非常接近。

表2 拉索參數[33]Table 2 Cable parameters

表3 模態頻率對比Table 3 Comparison of frequencies

4 參數分析

從式(17)可以看出,垂索-輔助索系統的模態和頻率取決于4 個關鍵參數,即:式(6)所示的εj,p、和αj以及由平衡方程引入的kj,p。當然,拉索間無量綱波速1/αj的差異以波速比參數η 計入,即:

因此,本節以3.1 小節的退化模型為對象,針對Irvine 參數、輔助索剛度k2,1和位置ε2,1、波速比η 四個關鍵參數開展參數分析,分別討論它們對系統固有特性的影響。

4.1 Irvine 參數

本小節考慮兩根拉索參數相同,其索力、單位長度質量及長度均取3.1 小節中數值,輔助索剛度k2,1=10,輔助索位置ε2,1=0.3。直接改變參數大小,研究對垂索-輔助索系統頻率及模態的影響。

圖2 曲線Fig.2 Curvesof

綜上,圖2(a)由于是同相振動,輔助索不起作用,所以表現出單一懸索的固有特性,ω-曲線出現“穿越”(cross-over)現象,而圖2(b)的反相振動中,輔助索起作用,客觀上破壞了模態的對稱性。這和斜拉索中重力效應下的模態對稱性破壞類似,因此ω-曲線表現為“轉向”(veering)現象。但若輔助索恰好位于模態節點處,則仍表現為“穿越”現象。

為進一步考察同相振動和反相振動固有特性間的關系,將圖2(a)和圖2(b)中第一個“穿越”點附近的頻率曲線繪于圖3(a)。同時,考慮輔助索剛度對一階反相振動頻率的影響,增加k2,1=500(相對剛性)時的一階反相振動頻率曲線。圖3(b)為圖3(a)中的部分典型模態。由圖3(a)可知:<39.48 時一階反相振動頻率大于一階同相振動,其原因是,反相振動中輔助索發揮了作用,輔助索剛度越大垂索-輔助索系統反相振動頻率越大。隨著逐漸增加,一階同相振動頻率曲線增長速度最快,輔助索剛度對索網反相振動頻率的貢獻越來越低,同時,一階反相振動模態形狀逐漸向二階轉變(圖3(b))。當=39.48 時,3 條一階頻率曲線與二階同相振動的頻率曲線相交于“穿越”點。從模態圖可知:此時對于一階反相振動,輔助索兩端正好處于索的不動點處(模態的節點處),故失去了其作用。因此,同相和反相振動頻率相等,不同輔助索剛度下的固有頻率也相等,從而出現4 條曲線交匯于同一點的現象。當>39.48 時,3 條一階振動頻率曲線均越過二階同相振動頻率曲線。此后,一階同相振動頻率均大于反相振動頻率。輔助索對一階反相振動頻率和模態重新起作用,因此k2,1=500 的頻率曲線一直處于k2,1=10 之上。值得注意的是,一階同相振動模態呈現3 個半波,而一階反相振動為2 個半波。

圖3 第一個“穿越”點附近頻率曲線和典型模態圖Fig.3 Frequency curves and typical modes near the first cross-over point

為進一步分析系統模態的局部振動模態特性,利用ZHOU 等[23]提出的衡量模態局部化指標方法,結合GATTULLI 等[34]在研究斜拉橋索-梁-塔耦合振動時,提出的衡量拉索局部振動的指標,給出以下反映索網系統模態局部化的指標計算公式:

式中:

為第j根索第p段的模態勢能與整個索網系統模態勢能之比。顯然,Θ 取值范圍為[0,1],Θ 越大,模態局部化程度越高。從圖3(b)可以看出:隨著的增大,一階同相模態局部化程度增大;一階反相模態的局部化程度則隨增大而降低;二階同相模態不受影響,Θ 保持不變。

4.2 輔助索位置ε2,1 及剛度k2,1

由理論推導可以看出,輔助索的位置及剛度是影響其作用效果的關鍵因素,為此通過改變式(17)中的參數k2,1和ε2,1來研究輔助索剛度及位置對系統頻率及模態的影響。兩根拉索的索力、質量及長度均取3.1 節中數值,另外取=2。

圖4 是系統前四階頻率隨輔助索位置及剛度變化而改變的曲面圖。由圖可知,各階同相振動頻率不受k2,1及ε2,1變化的影響,主要原因在于,同步振動時兩索間無相對位移,輔助索無作用。反相振動時,兩索振動時存在相對位移,且不同位置處大小不同,故輔助索的影響不盡相同。當輔助索位置ε2,1=z/(n+1)時(n為階次,z為小于n的正整數),即輔助索安裝于兩索反相振幅最大處時對系統反相頻率的影響最大。當然,即使輔助索設置在這些位置,隨著剛度的增大,各階反相振動頻率并不會無限增大,而是接近于高一階次的同相振動頻率,即輔助索的安裝位置和剛度最多能使系統無量綱頻率增大1。

圖4 ω 的三維曲面圖Fig.4 The 3-D surfaces of ω

另外,對于反相模態,系統可視為一個以輔助索中點處為對稱軸的對稱結構,故取其中一半結構則正好可退化至周??〉萚35]所研究的拉索-接地輔助索結構。因此,圖4 中各階反相模態頻率隨k2,1及ε2,1變化的規律與文獻[35]完全一致。

圖5 是圖4 中頻率對應的典型模態,對比圖5(a)和圖5(b)中的第三階反相振動模態可以發現,若采用大剛度輔助索,其將抑制索振動位移,從而提高系統的剛度;另外,當ε2,1=0.25時,對于第四階反相振動模態,輔助索兩端點間未產生相對位移,輔助索無作用,因此其頻率與四階同相振動的相同,Θ 值亦相同,僅是方向不同而已。值得注意的是,由于剛度的增大各階反相振動頻率并不會無限增大,而是接近于高一階次的同相振動頻率,當=0 時,若k2,1取無窮大,三階反相振動頻率也無限接近于四階同相振動頻率,從而導致表1 中6 階、7 階、8 階三個頻率均為5.4247 Hz。

圖5 典型模態Fig.5 Typical modes

4.3 波速比η

波速的差異使拉索間產生相對位移,而輔助索的作用將改變垂索-輔助索系統的頻率及模態。為此,通過改變η 的大小來研究拉索間參數差異對系統頻率和模態的影響。計算時索1 的索力、質量及長度均取3.1 節中數值,索2 的索力及長度與索1 相同,但變化索2 質量,同時為僅研究η 單一參數的影響,保證兩根拉索=2 不變。輔助索剛度和位置參數分別取為:k2,1=10,ε2,1=0.3。

計算得到的系統各階頻率ω隨η 的變化曲線及典型模態見圖6。由圖可知,隨著η 減小,即兩索參數差異變大,各階頻率逐漸增大并出現了跳階,階次越高的頻率曲線受到η 變化的影響越大,二階及以上頻率隨著η 的減小甚至出現了連續跳階現象,且階次越高連續跳躍的次數也越多。結合圖6(b)模態圖也可以看出,隨著η 減小,兩索間模態位移差異逐步增大,安裝輔助索后系統整體剛度和頻率增大越顯著。

圖6 ω-η 曲線圖及典型模態圖Fig.6 Curves of ω-η and typical modes

從Θ 值來看,二階同相和三階反相模態的局部化程度隨η 減小而不斷增加。從方程上看,η 減小本質上可以認為是索2 質量的減小,在不變的情況下相當于拉索2 垂度降低,等效軸向剛度提高[36-38],因此系統的整體剛度和各階固有頻率變大,且對高階頻率的影響更明顯。其力學原理與禹見達等[36-38]的復合阻尼索設計一致,本質上均是利用副索的支撐作用減小主索垂度,使主索在幾何上近乎為一條直線,從而在第一時間承受可能來自端部或橫向的載荷,客觀上達到提高系統剛度和頻率的效果。斜拉橋中相鄰拉索的波速往往存在差異,因此,安裝輔助索后系統剛度和頻率將增大,且對高階頻率的影響大于低階。

5 結論

本文推導了包含N根垂索和M道輔助索的索網動力學模型,給出了求解頻率方程及模態函數的通用表達式,并進行了驗證?;陔p索-輔助索退化模型進行參數分析,研究了索網系統的固有特性,得到如下主要結論:

(1) 由于本文對垂索-輔助索系統運動方程進行了徹底的無量綱化,故方程更具普適性,且關鍵參數及其物理意義更明晰。系統運動方程僅依賴于輔助索位置εj,p、Irvine 參數和αj,以及由平衡方程引入的輔助索剛度kj,p四個關鍵參數。

(2) 當兩索參數相同時,系統的同相振動頻率曲線以及輔助索恰好位于模態節點處時的反相振動頻率曲線隨Irvine 參數的增加表現為“穿越”現象,而其它反相振動頻率則表現為“轉向”現象。特定的Irvine 參數下,任意輔助索剛度對應的系統一階反相振動頻率均等于一階同相頻率,因此出現各頻率曲線在第一個“穿越”點交匯的現象。

(3) 當兩索參數相同時,輔助索僅影響反相振動頻率,且其安裝位置和剛度最多能使系統各階反相頻率增大1,因此,僅依賴輔助索剛度和位置來提高系統頻率的方式效果有限。

(4) 系統同相和反相頻率均隨兩索波速差異增大而發生往高階的“跳階”現象(頻率增大),且頻率階次越高“跳階”次數越多。實際斜拉橋的拉索波速往往不同,因此安裝輔助索后的垂索-輔助索系統高階頻率可能有明顯增大,故合理調節兩索的參數可在一定程度上突破僅增大輔助索剛度這一手段的限制,進一步提高系統剛度和頻率。

猜你喜歡
模態振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
重要模型『一線三等角』
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
重尾非線性自回歸模型自加權M-估計的漸近分布
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 国产成人无码播放| 国产一级精品毛片基地| 成人福利在线观看| 国产一级二级在线观看| 亚洲无码视频一区二区三区| 亚洲黄色片免费看| 国产精品99久久久| 欧美日韩综合网| 国产欧美日韩91| 日韩欧美国产综合| 亚洲天堂视频网站| 91精品在线视频观看| 视频一本大道香蕉久在线播放| 毛片基地美国正在播放亚洲 | 1769国产精品视频免费观看| 亚洲天堂伊人| 91精品国产91久无码网站| 美女一级免费毛片| 思思热精品在线8| 91精品情国产情侣高潮对白蜜| 国产精品香蕉在线| 国产嫩草在线观看| 1级黄色毛片| 欧美成人免费一区在线播放| 亚洲综合久久一本伊一区| 国产精品主播| 重口调教一区二区视频| 欧美一区国产| 无码一区中文字幕| 日韩精品亚洲人旧成在线| 国产69精品久久久久孕妇大杂乱| 欧美、日韩、国产综合一区| 狠狠色噜噜狠狠狠狠色综合久| 成人一级免费视频| 97综合久久| 国产农村1级毛片| 欧美日韩免费| 日韩在线视频网站| 91在线国内在线播放老师| 久久77777| 亚洲精品免费网站| 亚洲精品国产精品乱码不卞| 国产在线精彩视频论坛| 99性视频| 国产丝袜91| 日韩黄色在线| аv天堂最新中文在线| 成人亚洲视频| 国产噜噜在线视频观看| 国产精品短篇二区| 丰满少妇αⅴ无码区| 欧美成人免费| 怡红院美国分院一区二区| 69国产精品视频免费| jizz在线观看| 国产在线观看第二页| 亚洲国产综合自在线另类| 亚洲激情99| 国产乱人免费视频| 亚洲乱码在线播放| 一区二区在线视频免费观看| 女人18毛片一级毛片在线| 国产精品视频第一专区| 91色爱欧美精品www| 精品1区2区3区| 波多野吉衣一区二区三区av| 国产成人在线无码免费视频| 久热这里只有精品6| 伊人网址在线| 亚洲人成在线免费观看| 欧美成人亚洲综合精品欧美激情 | 亚洲综合18p| 日本免费一区视频| 91亚瑟视频| 国产精品对白刺激| 日韩精品无码不卡无码| 国产乱人伦精品一区二区| 亚洲va欧美va国产综合下载| 18禁黄无遮挡网站| 91在线国内在线播放老师| 亚洲视频a| 国产黄网永久免费|