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

基于衛(wèi)星在軌溫度預(yù)示熱控涂層性能退化的方法

2015-06-15 19:05:42宋馨張有為劉自軍斯東波
中國空間科學(xué)技術(shù) 2015年6期
關(guān)鍵詞:設(shè)備研究

宋馨 張有為 劉自軍 斯東波

(北京空間飛行器總體設(shè)計部,北京 100094)

基于衛(wèi)星在軌溫度預(yù)示熱控涂層性能退化的方法

宋馨 張有為 劉自軍 斯東波

(北京空間飛行器總體設(shè)計部,北京 100094)

介紹了一種利用衛(wèi)星在軌溫度數(shù)據(jù)預(yù)示熱控涂層在軌性能退化情況的方法,建立并推導(dǎo)了適用于研究熱控涂層隨整數(shù)年變化規(guī)律的數(shù)學(xué)模型,該模型利用衛(wèi)星設(shè)備的溫度曲線通過反演計算得到設(shè)備表面熱控涂層的退化數(shù)據(jù)。為檢驗本文數(shù)學(xué)模型的計算效果,構(gòu)造了一組數(shù)值試驗,結(jié)果表明本文數(shù)學(xué)模型的計算數(shù)據(jù)與真實值非常接近,說明本文建立的數(shù)學(xué)模型能夠很好地根據(jù)衛(wèi)星在軌溫度計算出熱控涂層在軌退化數(shù)據(jù)。

熱控涂層;性能退化;數(shù)值模擬;在軌數(shù)據(jù)反演;衛(wèi)星

1 引言

熱控涂層在衛(wèi)星熱控制領(lǐng)域得到了非常廣泛的應(yīng)用,熱控涂層的兩個重要熱輻射性能參數(shù)太陽吸收比和發(fā)射率決定了暴露于空間環(huán)境中衛(wèi)星表面的溫度水平。由于衛(wèi)星在軌運(yùn)行期間會受到原子氧、紫外線、帶電粒子等空間輻射的作用使熱控涂層的性能發(fā)生退化,主要表現(xiàn)為太陽吸收比增大[1],導(dǎo)致衛(wèi)星溫度升高,影響衛(wèi)星的正常運(yùn)行和在軌安全。因此,研究衛(wèi)星熱控涂層性能的退化規(guī)律,對于衛(wèi)星在軌安全運(yùn)行和衛(wèi)星熱控設(shè)計都非常重要。NASA通過在軌飛行試驗和地面模擬試驗研究了熱控涂層性能退化規(guī)律[2-3];我國在“東方紅二號”衛(wèi)星上對熱控涂層(OSR)的在軌性能退化情況進(jìn)行了測試[4];文獻(xiàn)[5-7]中研制了一套空間輻射環(huán)境模擬試驗系統(tǒng),對常用熱控涂層的性能退化進(jìn)行了大量試驗研究;文獻(xiàn)[8]對熱控涂層性能退化的機(jī)理進(jìn)行了理論研究。

由于在衛(wèi)星上搭載熱控涂層監(jiān)視裝置所消耗的資源代價非常大,目前的熱控涂層性能退化研究工作仍然以地面試驗和理論分析為主。熱控涂層性能的退化會直接反映到衛(wèi)星溫度水平的變化上,而我國發(fā)射的各類衛(wèi)星在軌飛行過程中又積累了大量溫度數(shù)據(jù),因此可選取不同時間衛(wèi)星在軌溫度數(shù)據(jù)并通過計算比較得到熱控涂層性能退化規(guī)律。本文從一維傳熱問題出發(fā),探討了基于衛(wèi)星在軌溫度數(shù)據(jù)通過反演計算研究熱控涂層在軌性能退化情況的方法,并通過模擬算例對計算效果進(jìn)行了驗證。

2 數(shù)學(xué)模型

2.1 一維傳熱方程

本文研究的為一維傳熱問題,研究對象為安裝在某人造地球衛(wèi)星上的一臺設(shè)備,設(shè)備的特征尺寸為L,在x=0處為絕熱邊界,在x=L處與宇宙空間輻射換熱,一維瞬態(tài)傳熱方程為

(1)

(2)

(3)

(4)

式中αk為研究對象與宇宙空間換熱表面的太陽吸收比;Es、Ee、Er分別為太陽輻射強(qiáng)度、地球紅外輻射強(qiáng)度以及地球反照的太陽輻射強(qiáng)度。

2.2 熱控涂層性能退化計算公式推導(dǎo)

為研究熱控涂層在軌性能退化,總體思路是從某時刻τ1開始,截取τ1~τ1+Δτ時段的研究對象在軌溫度變化曲線,反演得到這段時間內(nèi)研究對象的吸收外熱流qk,然后由式(4)計算得到αk,下面對式(4)中各項進(jìn)行逐個分析:

1)qk為研究對象與宇宙空間換熱表面的吸收外熱流值,在已知研究對象在軌溫度變化曲線的情況下,采用導(dǎo)熱反問題方法由研究對象的溫度曲線反演得到,為已知量,導(dǎo)熱反問題反演過程將在第2.3節(jié)中給出。

2)Es為太陽輻射強(qiáng)度,對于給定的軌道及軌道位置,太陽輻射強(qiáng)度的量值有成熟的公式可以計算得到,太陽輻射與研究對象受照表面的角度關(guān)系也能由衛(wèi)星軌道與太陽輻射的位置關(guān)系計算得到,為已知量。

3)Ee為地球紅外輻射強(qiáng)度,與大氣溫度、地表溫度、地形地貌等很多因素有關(guān),地球紅外輻射的測量及計算非常困難,為未知量。

4)Er為地球反照的太陽輻射強(qiáng)度,與地表熱輻射特性、海洋湖泊分布、大氣傳輸和散射特性等很多因素有關(guān),地球反照的太陽輻射強(qiáng)度的測量及計算非常困難,為未知量。

5)εk為研究對象與宇宙空間換熱表面的發(fā)射率,已有熱控涂層性能退化研究工作表明,熱控涂層在軌性能退化主要表現(xiàn)為太陽吸收比的變化,εk退化很小可視為不變[1],因此εk可以在衛(wèi)星發(fā)射前通過地面試驗測量得到,為已知量。

通過分析看到,式(4)中除了αk這個未知量外,還有地球紅外輻射強(qiáng)度Ee和地球反照的太陽輻射強(qiáng)度Er為未知量,下面通過引入一些新條件逐個消去這些未知量。為消去未知量Ee,引入一臺參照設(shè)備,設(shè)備的安裝位置和角度應(yīng)滿足兩個條件:

條件1:參照設(shè)備與宇宙空間輻射換熱表面的到達(dá)外熱流條件需要與研究對象的相同,根據(jù)衛(wèi)星上設(shè)備的安裝位置和布局篩選出符合條件的設(shè)備。

條件2:參照設(shè)備與宇宙空間輻射換熱表面的熱控涂層性能參數(shù)為已知量,這需要找到熱控涂層性能在軌退化很少、基本不變的涂層(如鋁合金黑色陽極氧化[1]),則在衛(wèi)星發(fā)射前通過地面試驗測量得到太陽吸收比和發(fā)射率的值;或者在軌退化規(guī)律已經(jīng)研究清楚、有確定數(shù)據(jù)的熱控涂層,這樣在衛(wèi)星在軌飛行期間參照設(shè)備的太陽吸收比和發(fā)射率為已知量。

(5)

(6)

(7)

消去地球紅外輻射強(qiáng)度Ee后,式(7)中還需要消去地球反照的太陽輻射強(qiáng)度Er。為此,選取另外一個時間段τ2~τ2+Δτ在軌溫度變化曲線。選取τ2~τ2+Δτ時段應(yīng)滿足兩個條件:

條件3:τ2~τ2+Δτ時段與τ1~τ1+Δτ時段在衛(wèi)星飛行軌道上的位置相同。

條件4:τ2~τ2+Δτ時段與τ1~τ1+Δτ時段的太陽對地球的輻照條件相同,這種情況往往出現(xiàn)在τ2~τ2+Δτ時段與τ1~τ1+Δτ時段的時間間隔相差為整數(shù)年的情況,因此這個條件也使本文提出的數(shù)學(xué)模型適合于研究熱控涂層性能隨整數(shù)年的退化規(guī)律。

(8)

(9)

(10)

根據(jù)選取τ2~τ2+Δτ時段的限制條件4,τ2~τ2+Δτ時段與τ1~τ1+Δτ時段的太陽對地球的輻照條件相同,可得:

(11)

另外由于影響地球反照太陽輻射強(qiáng)度的地表熱輻射特性、海洋湖泊分布、大氣傳輸和散射特性等因素在幾年或幾十年周期內(nèi)可視為狀態(tài)變化不大、比較穩(wěn)定,因此可近似認(rèn)為

(12)

(13)

2.3 導(dǎo)熱反問題

導(dǎo)熱反問題是利用研究對象的測量溫度,通過一定反演算法計算得到熱物性參數(shù)、邊界條件等未知量的一類問題。對于本文研究的一維瞬態(tài)傳熱問題,傳熱方程為

(14)

(15)

(16)

初始條件為

(17)

從式(16)中可以看到,由于邊界條件中引入了表征輻射熱流的4次方非線性項,大大增加了導(dǎo)熱反問題的不適定性,因此采用抗不適定性能力較強(qiáng)的共軛梯度法進(jìn)行求解。

共軛梯度法的目標(biāo)函數(shù)為

(18)

式中Tcal,n為設(shè)備溫度計算值;Tmea,n為設(shè)備溫度在軌遙測值;n為不同時間點(diǎn)。吸收外熱流值q的迭代式為

(19)

式中b為迭代次數(shù),d為搜索方向,其計算式為

(20)

γ由式(21)計算

(21)

迭代步長β為

(22)

(23)

邊界條件式(15)、式(16)對qn求微分得到:

(24)

(25)

式(25)中

(26)

靈敏度方程的初始條件由式(17)對qn求微分得到:

(27)

(28)

共軛梯度法的收斂目標(biāo)是使式(29)成立:

(29)

式中δ為很小的正數(shù)。共軛梯度法的求解步驟如下:

1)求解傳熱方程式(14)得到溫度計算值Tcal,n。

3)根據(jù)溫度計算值和在軌遙測值,檢查收斂目標(biāo)式(29)是否達(dá)到;如果達(dá)到收斂目標(biāo),則停止迭代,否則,轉(zhuǎn)入第4)步。

2.4 計算熱控涂層在軌性能退化的計算過程

計算熱控涂層在軌性能退化的步驟為:

3)求解導(dǎo)熱反問題,反演出研究對象和參照設(shè)備的溫度曲線在τ1~τ1+Δτ和τ2~τ2+Δτ時段的吸收外熱流曲線qk,τ1~τ1+Δτ、qc,τ1~τ1+Δτ、qk,τ2~τ2+Δτ、qc,τ2~τ2+Δτ。

3 數(shù)值驗證

下面構(gòu)造一組數(shù)值算例對數(shù)學(xué)模型的效果進(jìn)行檢驗。首先假設(shè)一組熱控涂層在軌退化情況的虛擬算例,計算得到熱控涂層退化前以及退化后衛(wèi)星設(shè)備的在軌溫度曲線,虛擬算例中的計算結(jié)果視為真實值;然后根據(jù)虛擬算例中的部分?jǐn)?shù)據(jù),通過前面推導(dǎo)的數(shù)學(xué)模型進(jìn)行反演計算,計算結(jié)果與虛擬算例中的真實值比較,檢驗數(shù)學(xué)模型的效果。

3.1 虛擬算例

選取衛(wèi)星上某設(shè)備做為研究對象,其傳熱過程可視為一維傳熱,特征尺寸方向L=5 mm,材料為鋁合金,設(shè)備內(nèi)表面為絕熱邊界條件;外表面與宇宙空間環(huán)境通過輻射交換熱量,表面狀態(tài)為噴涂SR107熱控白漆。

選取某設(shè)備做為參照設(shè)備,其傳熱過程為一維傳熱,特征尺寸方向L=8 mm,材料為鋁合金,設(shè)備內(nèi)表面為絕熱邊界條件,外表面與宇宙空間環(huán)境通過輻射交換熱量,表面狀態(tài)為鋁合金黑色陽極氧化。研究對象和參照設(shè)備的安裝位置和布局使研究對象外表面與參照設(shè)備外表面到達(dá)外熱流情況相同,符合條件1和條件2。

熱物理參數(shù)取值為鋁合金密度ρ=2 700 kg/m3、比熱容cp=900 J/(kg·K)、導(dǎo)熱系數(shù)k=120 W/(m·K)。熱輻射性能參數(shù)取值為:SR107熱控白漆的地面測量值αk,τ0=0.18、εk=0.85,剛?cè)胲壓笠欢螘r間SR107熱控白漆熱輻射性能無退化即αk,τ1~τ1+Δτ=0.18,5年后退化為αk,τ2~τ2+Δτ=0.3;鋁合金黑色陽極氧化的地面測量值αc,τ0=0.87、εc=0.87,在軌飛行過程中鋁合金黑色陽極氧化熱輻射性能無退化,即αc,τ1~τ1+Δτ=αc,τ2~τ2+Δτ=0.87。

圖1 到達(dá)外熱流曲線

選取發(fā)射入軌后不久的12 000s時間段做為τ1~τ1+Δτ段,并根據(jù)條件3和條件4選擇發(fā)射5年后的12 000s時間段做為τ2~τ2+Δτ時段,研究對象和參照設(shè)備外表面到達(dá)的太陽輻射、地球紅外輻射、地球反射太陽輻射如圖1所示。

根據(jù)研究對象和參照設(shè)備外表面無退化和5年退化后的表面熱控涂層的熱輻射參數(shù)可由圖1中各項到達(dá)外熱流值計算出τ1~τ1+Δτ時段以及τ2~τ2+Δτ時段外表面的吸收外熱值,圖2和圖3分別給出了研究對象無退化和5年退化后的吸收外熱流曲線。

圖2 研究對象吸收外熱流

圖3 參照設(shè)備吸收外熱流

圖4 單元劃分示意

研究對象和參照設(shè)備沿特征尺寸方向劃分為5個單元,測溫點(diǎn)粘貼在第3個單元上,如圖4所示。

求解式(1)~式(3)組成的一維瞬態(tài)傳熱方程即可計算得到研究對象和參照設(shè)備分別在τ1~τ1+Δτ時段以及τ2~τ2+Δτ時段的在軌溫度曲線,如圖5和圖6所示。

圖5 研究對象在軌溫度曲線

圖6 參照設(shè)備在軌溫度曲線

3.2 數(shù)學(xué)模型的檢驗

通過本文建立的數(shù)學(xué)模型,根據(jù)衛(wèi)星設(shè)備在軌溫度計算熱控涂層在軌性能退化情況。進(jìn)行計算時能夠獲得的全部數(shù)據(jù)為

1)已知條件1:研究對象,傳熱過程可視為一維傳熱,特征尺寸方向L=5mm,材料為鋁合金,內(nèi)表面為絕熱邊界條件,外表面噴涂SR107熱控白漆。

2)已知條件2:參照設(shè)備,傳熱過程可視為一維傳熱,特征尺寸方向L=8mm,材料為鋁合金,內(nèi)表面為絕熱邊界條件,外表面為鋁合金黑色陽極氧化。

3)已知條件3:地面通過試驗測量得到鋁合金密度ρ=2 700kg/m3、比熱容cp=900J/(kg·K)、導(dǎo)熱系數(shù)k=120W/(m·K)。

4)已知條件4:SR107熱控白漆的地面測量值αk,τ0=0.18、εk=0.85;鋁合金黑色陽極氧化的地面測量值αc,τ0=0.87、εc=0.87,在軌飛行過程中鋁合金黑色陽極氧化熱輻射性能無退化。

5)已知條件5:研究對象和參照設(shè)備在軌溫度曲線,由衛(wèi)星飛行過程中通過地面遙測得到。

按照2.4節(jié)計算熱控涂層在軌性能退化的步驟進(jìn)行計算分析:

1)根據(jù)已知條件1~已知條件4,可以得到步驟1)所需的各項數(shù)據(jù)。

2)根據(jù)已知條件5,從衛(wèi)星飛行遙測數(shù)據(jù)中找出在τ1~τ1+Δτ時段以及τ2~τ2+Δτ時段的在軌溫度曲線,即圖5和圖6。

3)根據(jù)圖5和圖6,進(jìn)行4次導(dǎo)熱反問題反演求解,設(shè)備初始溫度取圖5和圖6中0時刻的值,計算出所要研究的設(shè)備和參照設(shè)備的溫度曲線在τ1~τ1+Δτ時段以及τ2~τ2+Δτ時段的吸收外熱流曲線;從圖2和圖3可以看到,研究對象以及參照設(shè)備在無退化和5年退化后的吸收外熱流曲線波動變化都非常劇烈,并且輻射邊界條件中引入的4次非線性項更加增大了導(dǎo)熱反問題反演的難度;導(dǎo)熱反問題反演得到的qk,τ1~τ1+Δτ、qc,τ1~τ1+Δτ、qk,τ2~τ3+Δτ、qc,τ2~τ2+Δτ曲線與真實值(圖2和圖3)比較如圖7~圖10所示,從圖中可以看到,導(dǎo)熱反問題較好地反演出了吸收外熱流真實值,最大相對偏差為3.9%。

圖7 qk,τ1~τ1+Δτ反演值與真實值比較

圖8 qc,τ1~τ1+Δτ反演值與真實值比較

圖9 qk,τ2~τ2+Δτ導(dǎo)熱反問題反演值與真實值比較

圖10 qc,τ2~τ2+Δτ導(dǎo)熱反問題反演值與真實值比較

4 結(jié)束語

本文介紹了一種利用衛(wèi)星在軌溫度數(shù)據(jù)計算熱控涂層在軌性能退化情況的方法,建立并推導(dǎo)了適用于研究熱控涂層隨整數(shù)年變化規(guī)律的數(shù)學(xué)模型,該模型利用衛(wèi)星設(shè)備的溫度曲線通過反演計算得到設(shè)備表面熱控涂層的退化數(shù)據(jù)。為檢驗數(shù)學(xué)模型的計算效果,構(gòu)造了一組數(shù)值試驗,先假設(shè)出熱控涂層退化前和退化后的真實值,分別計算出不同熱控涂層參數(shù)下衛(wèi)星設(shè)備的溫度曲線,然后采用本文模型對設(shè)備溫度曲線進(jìn)行反演計算,得到熱控涂層退化數(shù)據(jù),再與真實值比較。數(shù)值試驗結(jié)果表明數(shù)學(xué)模型計算數(shù)據(jù)與真實值非常接近,說明本文建立的數(shù)學(xué)模型能夠很好地根據(jù)衛(wèi)星在軌溫度計算出熱控涂層在軌退化數(shù)據(jù)。

[1]閔桂榮,郭舜.航天器熱控制[M].北京:科學(xué)出版社,1998:60-100.

[2]DONALD A J, GEORGE C T, DAVID T W, et al. Optical properties of thermal control coatings after weathering, simulated ascent heating and simulated space radiation exposure[R]. NASA 2008-215259.

[3]KRISTAN P G, RACHID D, GORDEN V. Reduced polarized thermal emission due to surface contamination by dew and atmospheric aerosols[J]. Journal of Quantitative Spectroscopy Radiative Transfer, 2004, 88(1):64-70.

[4]馮偉泉,丁義剛,閆德葵,等.地球同步軌道長壽命衛(wèi)星熱控涂層太陽吸收率性能退化研究[J].中國空間科學(xué)技術(shù),2005,25(2):34-40.

FENG WEIQUAN, DING YIGANG, YAN DEKUI, et a1.Study on long-term degradation of solar absorptance properties of GEO stationary satellite thermal control coating[J].Chinese Space Science and Technology, 2005,25(2):34-40.

[5]劉宇明,馮偉泉,丁義剛.SR107-ZK在空間輻照下的物性變化[J].宇航材料工藝,2007,37(4):61-64.

LIU YUMING, FENG WEIQUAN, DING YIGANG. Space radiation effects on SR107-ZK[J]. Aerospace Materials & Technology, 2007,37(4):61-64.

[6]馮偉泉,丁義剛,閆德葵.空間電子、質(zhì)子和紫外綜合輻照模擬試驗研究[J]. 航天器環(huán)境工程,2005,22(2):69-72.

FENG WEIQUAN,DING YIGANG, YAN DEKUI. Study on space electron, proton and ultraviolet combined irradiation simulation test[J]. Spacecraft Environment Engineering, 2005,22(2):69-72.

[7]馮偉泉,王榮,丁義剛.低能質(zhì)子對衛(wèi)星熱控涂層太陽吸收率的影響[J].中國空間科學(xué)技術(shù),2007,27(6):42-47.

FENG WEIQUAN, WANG RONG,DING YIGANG.Low energy proton effects on solar absorptance of thermal control coatings[J]. Chinese Space Science and Technology,2007,27(6):42-47.

[8]TAN H P, LUO J F, XIA X L. Transient coupled radiation and conduction in a three layer composite with semitransparent specular interfaces and surfaces[J].Journal of Heat Transfer, 2002, 124(3):470-481.

(編輯:楊嬋)

《中國空間科學(xué)技術(shù)》編輯部新動向

為提高論文發(fā)表速度、縮短科研成果公布的時間,《中國空間科學(xué)技術(shù)》與 “中國知網(wǎng)”簽署了優(yōu)先數(shù)字出版協(xié)議,文章一經(jīng)錄用即可實行優(yōu)先數(shù)字出版。

另外,從2016年第1期開始,本刊將采用雙欄排版。同時,新的版式每篇文章中添加了卷、通訊作者、引用格式等信息,并參照國標(biāo) GB/T 7714-2015對參考文獻(xiàn)的格式進(jìn)行了修改,方便閱讀、引用和檢索。

Prediction on Performance Degradation of Thermal Control Coating Based on Satellite Temperature on Orbit

SONG Xin ZHANG Youwei LIU Zijun SI Dongbo

(Beijing Institute of Spacecraft System Engineering,Beijing 100094)

A method was proposed to predict the performance degradation of thermal control coating based on temperature data on orbit.A mathematical model which was constructed to study the degradation of thermal control coating performance spanning integer years was verified by numerical simulations. The calculation data meet the true value well, indicating that the mathematical model is valid.

Thermal control coating; Performance degradation; Numerical simulation;Inversion of telemetry data;Satellite

2015-02-01。收修改稿日期:2015-06-01

10.3780/j.issn.1000-758X.2015.06.005

宋馨 1986年生,2010年獲北京航空航天大學(xué)人機(jī)與環(huán)境工程專業(yè)碩士學(xué)位。研究方向為航天器熱控技術(shù)。

猜你喜歡
設(shè)備研究
諧響應(yīng)分析在設(shè)備減振中的應(yīng)用
FMS與YBT相關(guān)性的實證研究
2020年國內(nèi)翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設(shè)計中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統(tǒng)研究
基于VB6.0+Access2010開發(fā)的設(shè)備管理信息系統(tǒng)
基于MPU6050簡單控制設(shè)備
電子制作(2018年11期)2018-08-04 03:26:08
新版C-NCAP側(cè)面碰撞假人損傷研究
500kV輸變電設(shè)備運(yùn)行維護(hù)探討
主站蜘蛛池模板: 久久国产亚洲偷自| 国产精品毛片一区| 国产网友愉拍精品| 99热精品久久| 日韩欧美一区在线观看| 69国产精品视频免费| 国产v精品成人免费视频71pao| 国产SUV精品一区二区6| 第一页亚洲| 国产日韩欧美在线视频免费观看| 亚洲国产中文在线二区三区免| 一级看片免费视频| 精品欧美视频| 国产成人综合在线视频| 久草热视频在线| 免费在线不卡视频| 国产精品蜜臀| 国产成人精品视频一区视频二区| 久久久久88色偷偷| 欧美人与性动交a欧美精品| 99精品国产电影| 在线视频亚洲色图| 国产精品亚洲一区二区三区z| 欧美福利在线播放| 国产丝袜无码精品| 日本欧美视频在线观看| 日韩一区二区在线电影| 午夜视频免费一区二区在线看| 99热这里都是国产精品| 又大又硬又爽免费视频| 亚洲色大成网站www国产| 国产精品偷伦视频免费观看国产| 无码专区第一页| 免费人成视网站在线不卡| 久久先锋资源| 免费高清毛片| 无码啪啪精品天堂浪潮av| 一本大道视频精品人妻| 免费国产不卡午夜福在线观看| 欧美国产日韩一区二区三区精品影视| 国产原创演绎剧情有字幕的| 无码乱人伦一区二区亚洲一| 亚洲精品另类| 蜜芽国产尤物av尤物在线看| 精品久久综合1区2区3区激情| 国产精品开放后亚洲| 久久中文无码精品| 日韩123欧美字幕| www亚洲精品| 成人午夜天| 亚洲成人精品| 亚洲午夜国产精品无卡| 直接黄91麻豆网站| 国产一级在线观看www色| 谁有在线观看日韩亚洲最新视频| 性做久久久久久久免费看| 国产不卡网| 网久久综合| 丰满人妻一区二区三区视频| 国产成人精品一区二区| 一级全免费视频播放| 国产色婷婷视频在线观看| 怡春院欧美一区二区三区免费| 亚洲一欧洲中文字幕在线| 中文字幕精品一区二区三区视频 | 国产杨幂丝袜av在线播放| 久久久久青草线综合超碰| 亚洲精品视频免费观看| 精品1区2区3区| 福利在线不卡| 国产人妖视频一区在线观看| 91无码人妻精品一区二区蜜桃| 国产三级韩国三级理| 亚洲国产成人综合精品2020| 国产91小视频在线观看| 日韩123欧美字幕| 九色综合视频网| 天天躁夜夜躁狠狠躁躁88| 中文字幕久久亚洲一区| 高h视频在线| 91精品视频网站| 日韩123欧美字幕|