亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種再入飛行器結(jié)構(gòu)的熱流載荷反演方法

        2022-05-12 05:30:36劉子昂吳邵慶李彥斌費慶國
        宇航學報 2022年3期
        關(guān)鍵詞:共軛測量誤差熱流

        劉子昂,陳 強,吳邵慶,李彥斌,費慶國

        (1. 東南大學機械工程學院,南京 211189;2. 東南大學江蘇省空天機械裝備工程研究中心,南京 211189)

        0 引 言

        再入飛行器在進出大氣層時,其結(jié)構(gòu)與外側(cè)空氣相互摩擦,表面將產(chǎn)生大量熱能,形成強烈的氣動加熱效應(yīng)。結(jié)構(gòu)表面的熱流會沿飛行器結(jié)構(gòu)向內(nèi)傳遞,進而影響飛行器結(jié)構(gòu)靜、動力學特性,嚴重時可導致結(jié)構(gòu)發(fā)生強度破壞。因此,準確掌握再入飛行器結(jié)構(gòu)在實際服役過程中所面臨的熱流載荷,對飛行器結(jié)構(gòu)設(shè)計和評估至關(guān)重要。然而,直接測量服役狀態(tài)下再入飛行器結(jié)構(gòu)表面的熱流載荷存在諸多難點:一方面,再入飛行器服役于嚴酷的力/熱/振動/噪聲等復雜環(huán)境,將熱流傳感器布置于飛行器外側(cè)容易遭受多場載荷環(huán)境影響;另一方面,在再入飛行器表面的熱防護結(jié)構(gòu)中集成傳感器的難度較大。預埋測溫晶粒為一種對熱防護結(jié)構(gòu)影響小的測量方式,但此方法只能記錄服役過程中的最高溫度,難以準確反映溫度隨時間變化的特征。近年來所發(fā)展的載荷反演方法通過在再入飛行器內(nèi)部布置溫度傳感器并測量其內(nèi)壁面溫度,進而求解傳熱學反問題(Inverse heat conduction problem, IHCP)得到外壁面熱載荷,成為獲取再入飛行器熱流載荷的有效手段。

        傳熱學反問題作為一類典型的不適定性問題,對定解數(shù)據(jù)的擾動極為敏感,國內(nèi)外學者已通過各類數(shù)值方法開展了相關(guān)研究。錢煒祺等采用順序函數(shù)法對簡單平板結(jié)構(gòu)表面瞬變熱流進行反演,以當前時刻后若干時刻的結(jié)果對當前時刻熱流進行反演,反演結(jié)果依賴于用于反演的時刻數(shù)量。Lu等基于共軛梯度法,利用彎管外側(cè)布置的溫度傳感器獲得的數(shù)據(jù)對薄壁彎管結(jié)構(gòu)內(nèi)部流體溫度進行反演;在彎管內(nèi)外溫差不大的情況下共軛梯度法有著較好的反演結(jié)果。相比于順序函數(shù)法,共軛梯度法中沒有人為設(shè)定的參數(shù),并且反演精度和計算速度與順序函數(shù)法相當,因此共軛梯度法有著更好的適用性。黃少君等采用隨機優(yōu)化算法對熱源位置進行反演。隨機優(yōu)化算法相比于梯度方法在跳出局部最優(yōu)解上有著一定的優(yōu)勢,但由于其計算量大、計算時間長,難以用于復雜結(jié)構(gòu)上的熱流載荷反演。Lu等基于Tikhonov正則化方法以及共軛梯度法對三維瞬態(tài)熱傳導分析問題進行反演,通過正則化方法降低測量誤差對反演結(jié)果的影響。對于較厚的結(jié)構(gòu),Khajehpour等采用域分解法,將結(jié)構(gòu)沿厚度方向劃分為若干個區(qū)域,分區(qū)域?qū)崃鬟M行識別,識別結(jié)果相比于對整個結(jié)構(gòu)直接識別有更高的精度。現(xiàn)有研究大多針對簡單的厚板、圓管等結(jié)構(gòu),采用一維傳熱模型進行熱載荷反演,僅能考慮結(jié)構(gòu)沿厚度方向的熱量傳遞。而對于實際結(jié)構(gòu),熱載荷不僅會沿厚度方向傳遞,同樣會由于外表面溫度的非均勻性使得熱量沿面內(nèi)方向傳遞。再入飛行器尺寸較大,并且外形復雜。采用再入飛行器結(jié)構(gòu)簡化傳熱模型時,其在實際工況下的計算結(jié)果將存在較大誤差。

        此外,現(xiàn)有研究大多采用較為簡單的熱載荷進行理論研究和仿真分析,取得了較好的識別效果,但這些方法在復雜熱載荷工況下的穩(wěn)定性和收斂性有待研究。Nakamura等基于順序函數(shù)法以及奇異值分解法對飛行器再入階段表面熱流進行反演;通過在計算中引入奇異值分解法降低了測量噪聲對反演結(jié)果的影響,使得反演結(jié)果不會由于過度擬合從而偏離真實值。邵元培等針對飛行器中的燒蝕情況,考慮結(jié)構(gòu)表面幾何域的變化反演其表面熱流。

        對于再入飛行器而言,其服役過程中結(jié)構(gòu)處于上千攝氏度的高溫中。此時材料的熱導率、比熱容相比于常溫下有著很大的差異,材料熱物性參數(shù)隨溫度的變化會引起傳熱的非線性特性。薛齊文等基于T形板結(jié)構(gòu),在考慮材料溫變熱物性參數(shù)的基礎(chǔ)上,用共軛梯度法對熱物性參數(shù)以及邊界條件同時進行反演,基于數(shù)值仿真分析的手段探討了測點數(shù)目、測量誤差和變量初值對反演結(jié)果的影響。Bergagi等基于共軛梯度法求解邊界上的溫度,考慮了材料的溫變熱物性參數(shù)對反演結(jié)果的影響,證明了此方法可以較好地反演薄壁結(jié)構(gòu)的溫度邊界載荷。Xiong等結(jié)合共軛梯度法以及順序函數(shù)法的優(yōu)點,提出了順序共軛梯度法,當測量數(shù)據(jù)受較大噪聲影響時,此方法具有高于傳統(tǒng)順序函數(shù)法的計算精度以及高于傳統(tǒng)共軛梯度法的計算速度。

        本文針對再入飛行器的熱流載荷反演問題開展研究。首先,介紹基于共軛梯度法的熱流載荷反演理論基礎(chǔ)。進而,針對一維結(jié)構(gòu)開展熱流載荷反演的數(shù)值分析,討論材料溫變熱物性參數(shù)、溫度傳感器測量誤差等因素對反演結(jié)果的影響。然后,針對典型熱防護結(jié)構(gòu)和返回艙結(jié)構(gòu),結(jié)合有限元方法開展熱流載荷反演的數(shù)值分析,并討論測點位置對反演結(jié)果的影響。

        1 熱流載荷反演的理論基礎(chǔ)

        再入飛行器在服役過程中面臨復雜的熱流載荷,為了更好地擬合熱流載荷變化,將整個時程均勻離散為若干個時程節(jié)點,其中各個節(jié)點之間的熱流載荷由兩個相鄰節(jié)點的熱流載荷差值獲得。

        假定結(jié)構(gòu)所受熱流時程曲線函數(shù)為(),分布規(guī)律的空間函數(shù)為(),則全域上熱流載荷為:

        (,)=()()

        (1)

        式中:為結(jié)構(gòu)熱流載荷所處邊界上的位置矢量。

        將傳熱過程均勻劃分為個時間節(jié)點。在結(jié)構(gòu)上選取個位置作為模擬布置傳感器的位置,所受熱流載荷邊界上選取個位置。以此個位置的溫度時程曲線為已知值,以個位置的熱流時程曲線為反演值,建立優(yōu)化目標函數(shù):

        (2)

        式中:mea,,為第個傳感器測量點時刻的測量溫度,cal,,為第個傳感器測量點時刻的反演結(jié)果計算值。當小于極小值之后,認為其收斂,得到反演結(jié)果。

        ,為第個點的熱流時程曲線時刻的熱流值,則目標函數(shù)對未知參量,求偏導數(shù),可得到目標函數(shù)的梯度為:

        =1,2,…,;=1,2,…,

        (3)

        熱流時程曲線參量的迭代式為:

        (,)+1=(,)-(,),

        =1,2,…,;=1,2,…,

        (4)

        式中:表示迭代步數(shù),(,)表示第次迭代計算中獲得第個點的熱流時程曲線時刻的熱流值,表示迭代步長,,表示對第個熱流時程曲線時刻熱流值的迭代搜索方向。,由下式獲得:

        =1,2,…,;=1,2,…,

        (5)

        (6)

        式中:為共軛系數(shù);當=0時,=0。

        計算中的迭代步長為:

        =

        (7)

        2 熱流載荷反演驗證算例研究

        為驗證熱流反演方法的準確性,并探討影響反演結(jié)果的因素。針對如圖1所示的一維傳熱模型,采用有限體積法(Finite volume method, FVM)進行瞬態(tài)傳熱正問題分析計算。以結(jié)構(gòu)右側(cè)的仿真溫度結(jié)果,反演其左側(cè)的熱流載荷。

        2.1 有限體積法瞬態(tài)傳熱計算

        對于一維情況下采用有限體積法將總長度為的區(qū)域離散為長度為Δ的個控制體,如圖1所示,其中=0處為熱流入射面,=處為絕熱面,以此處的仿真結(jié)果模擬溫度傳感器的測量值。

        圖1 有限體積法中的控制體離散Fig.1 Discretization of the control volume in finite volume method

        依據(jù)能量守恒定律,對于每個控制體有:

        (-)Δ=Δ

        (8)

        式中:為單個控制體的輸入熱流,為單個控制體的輸出熱流,為控制體面積,為材料密度,為材料比熱容,Δ為時間間隔。

        對于控制體,,…,-1其輸入輸出熱流為:

        (9)

        式中:為材料熱導率,為第個控制體的溫度。

        聯(lián)立式(8)和式(9),可得:

        (10)

        定義式(11)所示的傅里葉數(shù),式(10)可簡化為:

        (11)

        ()(+1+-1-2)=1

        (12)

        依據(jù)實際工況設(shè)定第1個控制體與第個控制體為第二類邊界條件??刂企w1邊界為外界熱流輸入面,控制體為絕熱邊界,則邊界條件為:

        (13)

        式中:為環(huán)境輸入熱流。

        2.2 熱流載荷反演算例

        對于如圖1所示的一維結(jié)構(gòu)在=0處施加如式(14)所示熱流載荷:

        ()=-+300,0≤≤300

        (14)

        基于有限體積法計算熱傳導正問題獲取控制體的仿真值。以此仿真值模擬結(jié)構(gòu)內(nèi)部的實際測量溫度值。通過共軛梯度法對外側(cè)熱流載荷進行反演,計算流程圖如圖2所示。分別考慮材料熱物性參數(shù)為常數(shù)和溫變參數(shù)時的熱流反演結(jié)果,如圖3所示。

        反演熱流載荷的平均誤差為:

        (15)

        式中:為真實熱流值,為反演熱流值。

        圖2 共軛梯度法反演熱流載荷流程圖Fig.2 Flow chart of heat flux identification using conjugate gradient method

        由圖3可知,材料熱物性為常數(shù)時,反演熱流平均誤差為6.84 mW·mm;材料熱物性為溫變參數(shù)時,反演熱流平均誤差為3.88 mW·mm。當考慮材料熱物性參數(shù)為常數(shù)時,其反演誤差明顯高于考慮材料熱物性為溫變參數(shù)的結(jié)果。同時,這一誤差會隨著結(jié)構(gòu)整體溫度的升高而增大。這是因為材料在高溫時的熱導率和比熱容與常溫時相比有著很大差異。因此很有必要在熱流載荷反演中考慮材料熱物性隨溫度的變化特征。

        圖3 熱流載荷的反演結(jié)果Fig.3 Identification results of heat flux

        為模擬溫度傳感器的測量誤差,在仿真溫度值的基礎(chǔ)上施加方差為0.5、均值為0的隨機誤差。依據(jù)3σ準則,其誤差所在區(qū)間為[-2.12 K, 2.12 K]。由于反演結(jié)果對測量誤差的敏感性較高,且容易產(chǎn)生過擬合的情況。因此,需要對加噪聲之后的數(shù)據(jù)進行降噪處理。采用降噪處理后的數(shù)據(jù)與未降噪處理的模擬溫度測量數(shù)據(jù)對熱流載荷進行反演,反演得到的熱流載荷如圖4所示。

        圖4 測量數(shù)據(jù)含噪聲時熱流反演結(jié)果Fig.4 Identified results of heat flux with noisy measurement data

        表1中列出了在模擬測量結(jié)果上考慮不同測量誤差情況下熱流反演結(jié)果的誤差。由表1和圖4可知,當考慮溫度傳感器測量誤差時,其反演結(jié)果會受到較大的影響,特別是這種誤差會導致反演的結(jié)果在真值附近出現(xiàn)較大的震蕩。為了減小這種由于測量誤差所造成的反演誤差,一方面需要采用降噪方法,降低測量誤差對反演結(jié)果的影響;另一方面需要制定合適的收斂準則,使其不會因為迭代次數(shù)過多產(chǎn)生對測量數(shù)據(jù)的過度擬合。

        表1 測量數(shù)據(jù)含噪聲時熱流反演誤差Table 1 Error in identified heat flux with noisy measurement data

        以測量噪聲方差為0.5,均值為0的反演結(jié)果為例,分析反演過程中目標函數(shù)與熱流誤差的收斂過程,如圖5所示。其中熱流反演平均誤差以式(15)計算。

        圖5 反演結(jié)果隨迭代次數(shù)的變化Fig.5 Variation of the identified results with respect to the number of iterations

        當考慮溫度傳感器測量誤差時,隨著迭代次數(shù)的增加,目標函數(shù)值不斷減小,但熱流反演結(jié)果的誤差卻不會一直減小。這說明目標函數(shù)值越小并不等同于反演結(jié)果越好。由于測量誤差的存在,會導致共軛梯度法計算過程中各個時刻的搜尋方向出現(xiàn)誤差,從而影響最終的反演結(jié)果。因此,在制定共軛梯度法的收斂準則時,不僅要設(shè)置目標函數(shù)達到的最小值,也需要考慮下降速率。當目標函數(shù)下降速率較低時需要停止計算,防止過擬合。

        3 實際結(jié)構(gòu)算例研究

        為驗證熱流反演方法在再入飛行器上的運用,并探討影響反演結(jié)果的因素。本節(jié)首先針對夾芯式熱防護結(jié)構(gòu),反演二維結(jié)構(gòu)表面均勻分布的熱流載荷;進而,針對典型再入飛行器返回艙,反演三維結(jié)構(gòu)上所受熱流載荷。

        3.1 熱流載荷

        再入飛行器在服役過程中所受熱流載荷取決于飛行器的形狀、上升軌跡和再入條件等因素。參考文獻[21]中采用了如圖6中真實值曲線所示的一種兩級入軌飛行器在概念設(shè)計中采用的熱流載荷。為使熱流反演更貼近實際情況,在實際結(jié)構(gòu)熱流反演計算中采用此熱流載荷作為結(jié)構(gòu)的熱流載荷。

        圖6 含測量誤差時熱流反演結(jié)果Fig.6 Identified results of heat flux with measurement error

        3.2 夾芯式熱防護結(jié)構(gòu)熱流載荷反演

        本文所研究的夾芯式熱防護結(jié)構(gòu),其從外至內(nèi)分別包括復合材料面板結(jié)構(gòu)、氣凝膠結(jié)構(gòu)、復合材料面板結(jié)構(gòu)、膠層以及內(nèi)側(cè)金屬底板。各層的厚度依次為2 mm、28 mm、0.5 mm、2 mm、5 mm,結(jié)構(gòu)有限元模型如圖7所示。

        圖7 夾芯式熱防護結(jié)構(gòu)有限元模型Fig.7 Finite element model of the sandwich thermal-protection structure

        在實際工程應(yīng)用中,其結(jié)構(gòu)內(nèi)部難以布置傳感器,因此只能在其金屬底板內(nèi)側(cè)布置傳感器。其外表面受到空間上均勻分布的熱流。以底部溫度傳感器獲得的溫度值,對外部熱流載荷進行反演,獲取結(jié)構(gòu)外部施加的熱流載荷。

        當不考慮傳感器的測量誤差時,以底部金屬板獲得的數(shù)值仿真結(jié)果作為已知條件,反演其頂部輸入熱流,此時的熱流反演平均誤差為0.362 mW·mm。當使用不含噪聲的溫度測量數(shù)據(jù)反演時,熱流的反演結(jié)果不會出現(xiàn)過擬合的情況,因此可以在設(shè)置收斂準則時設(shè)置一個較小的值,增加迭代次數(shù)。最終的反演結(jié)果能夠很好的反映實際熱流載荷。

        在此基礎(chǔ)上,考慮測量中由于傳感器測量誤差所造成的數(shù)據(jù)波動。分別在原模擬測量溫度結(jié)果上加上方差分別為0.01、0.1,均值為0的正態(tài)分布隨機誤差。依據(jù)3σ準則,此兩工況下誤差分別在[-0.3 K, 0.3 K]與[-0.95 K, 0.95 K]之內(nèi)。將模擬測量值進行降噪處理前后的結(jié)果分別用于反演其外表面所受熱流,結(jié)果如圖6所示。在考慮測量誤差時,反演結(jié)果在整體的熱流趨勢上反演結(jié)果較好,但對于部分熱流峰值的反演存在誤差。

        由表2可見,當不考慮傳感器測量誤差時,共軛梯度法能夠有效地反演對復雜結(jié)構(gòu)的熱流載荷。而考慮傳感器測量誤差時,由于復合材料夾芯結(jié)構(gòu)本身良好的隔熱能力,溫度傳感器獲取到的溫度測量值較低,此時測量誤差對于反演的結(jié)果有著很大的影響。因此,為了保證熱流反演結(jié)果的精度,在布置傳感器位置時應(yīng)盡量選擇溫度較高的區(qū)域,以獲取較大的靈敏度。

        表2 含測量誤差時熱流反演平均誤差Table 2 Average error in heat flux identification with measurement errors

        3.3 再入飛行器結(jié)構(gòu)熱流反演

        返回艙結(jié)構(gòu)在再入大氣過程中,與空氣摩擦產(chǎn)生大量熱量。本文以日本宇宙航空研究開發(fā)機構(gòu)的PARTT返回艙結(jié)構(gòu)為例展開研究,其結(jié)構(gòu)如圖8所示,分為返回艙主體、熱防護結(jié)構(gòu)、支撐結(jié)構(gòu)和殼體四部分。并在結(jié)構(gòu)內(nèi)側(cè)布置模擬傳感器測點A與B,分別位于熱防護結(jié)構(gòu)的內(nèi)側(cè)與外側(cè)。

        圖8 返回艙結(jié)構(gòu)模型Fig.8 Structural model of the reentry capsule

        對于復雜三維結(jié)構(gòu),由于其內(nèi)側(cè)各部分溫差較大。在此采用不同位置的溫度測量值作為輸入,討論測點位置對熱流反演結(jié)果的影響。分別采用測點A與B的模擬測量結(jié)果作為輸入,反演其所受熱流載荷。反演結(jié)果如圖9所示。當采用測點A的溫度值作為反演輸入數(shù)據(jù)時,熱流反演平均誤差為2.22 mW·mm;當采用測點B的溫度值作為反演輸入數(shù)據(jù)時,反演平均誤差為1.30 mW·mm。

        圖9 測點不同時熱流反演結(jié)果Fig.9 Identified results of heat flux with different sensor arrangement

        由圖9可見,當測點位置越靠近熱流邊界位置時,反演效果也越好。這是由于越靠近熱載荷的位置,其結(jié)構(gòu)溫度越高,即靈敏度相對較高。此時誤差相對于有效測量值的比例較低,測量誤差對于反演結(jié)果的影響較小。因此在布置傳感器時,需要盡量使其位置靠近熱流邊界。

        在測點B處的結(jié)果上加上正態(tài)分布的隨機誤差模擬傳感器的測量誤差。分別加方差為1與2的誤差,依據(jù)3σ準則,其誤差所在區(qū)間分別為[-3 K, 3 K]與[-4.24 K, 4.24 K]。反演誤差見表3,反演結(jié)果如圖10所示。在考慮較大的測量誤差下,熱流反演結(jié)果與真實值相差較小,說明該反演方法具有良好的抗噪性。

        表3 測量誤差對熱流反演結(jié)果的影響Table 3 Influence of measured error on the heat flux identification results

        圖10 不同測量誤差時熱流反演結(jié)果Fig.10 Identified results of heat flux with different measurement errors

        4 結(jié) 論

        本文針對再入飛行器結(jié)構(gòu)開展熱流載荷反演問題研究,采用共軛梯度法求解傳熱學反問題。在已知熱流載荷作用區(qū)域和準確結(jié)構(gòu)傳熱模型的基礎(chǔ)上,通過結(jié)構(gòu)內(nèi)部區(qū)域布置的測溫傳感器反演了再入飛行器結(jié)構(gòu)外部熱流載荷變化過程。在此基礎(chǔ)上探討了材料熱物性參數(shù)的溫變特性、測點位置、傳感器誤差等因素對反演結(jié)果的影響,主要研究結(jié)論為:

        1)共軛梯度法能夠克服熱流反演問題中的不適定性,可以較好地反演典型熱防護結(jié)構(gòu)和返回艙結(jié)構(gòu)所受的熱流載荷;

        2)對于結(jié)構(gòu)內(nèi)部溫度差異較大的結(jié)構(gòu),其熱載荷反演過程中呈現(xiàn)出的非線性問題,通過修正靈敏度可以有效地提高識別精度。針對本文所研究的夾芯式熱防護結(jié)構(gòu),考慮材料熱物性為溫變參數(shù)時,熱流反演的平均誤差從6.84 mW·mm降至3.88 mW·mm;

        3)當反演復雜結(jié)構(gòu)的熱流載荷時,結(jié)構(gòu)溫度測點位置應(yīng)盡量靠近主要入射熱流位置,以提高計算中的靈敏度,從而提高反演精度。對于本文所研究的返回艙結(jié)構(gòu),當溫度測點布置于熱防護外側(cè)時相比于測點布置于熱防護結(jié)構(gòu)內(nèi)部,熱流反演誤差從2.22 mW·mm下降至1.30 mW·mm;

        4)為了降低測量誤差對反演結(jié)果的影響,需要采用合適的降噪手段并制定合適的收斂準則,從而防止反演結(jié)果出現(xiàn)過擬合。

        猜你喜歡
        共軛測量誤差熱流
        一個帶重啟步的改進PRP型譜共軛梯度法
        密度測量誤差分析
        一個改進的WYL型三項共軛梯度法
        巧用共軛妙解題
        一種自適應(yīng)Dai-Liao共軛梯度法
        縱向數(shù)據(jù)下變系數(shù)測量誤差模型的漸近估計
        內(nèi)傾斜護幫結(jié)構(gòu)控釋注水漏斗熱流道注塑模具
        空調(diào)溫控器上蓋熱流道注塑模具設(shè)計
        聚合物微型零件的熱流固耦合變形特性
        中國塑料(2017年2期)2017-05-17 06:13:24
        牽引變壓器功率測量誤差分析
        中文字幕一区二区人妻痴汉电车| 亚洲爆乳精品无码一区二区| 亚洲a∨无码一区二区| 国产精品麻豆A啊在线观看| 亚州中文字幕乱码中文字幕| 国产极品粉嫩福利姬萌白酱| 精品人妻人人做人人爽| 国产免费看网站v片不遮挡| 伊人久久综合狼伊人久久| 国产综合精品久久99之一| 日日摸天天摸人人看| 国产偷国产偷高清精品| 一本久道在线视频播放| 亚洲 欧美 偷自乱 图片| 思思久久96热在精品国产| 色综合久久加勒比高清88| 国产精品中文字幕日韩精品| 妺妺窝人体色www婷婷| 日韩好片一区二区在线看| 免费视频成人 国产精品网站| 少妇被粗大猛进进出出男女片| 免费人成年激情视频在线观看| 久久精品岛国av一区二区无码 | 日韩人妻无码精品二专区| 91九色国产老熟女视频| 精品视频无码一区二区三区| 亚洲va在线va天堂va手机| 国产在线拍91揄自揄视精品91| 国产一级内射视频在线观看| 9lporm自拍视频区| 久久综合给合久久狠狠狠9 | 美腿丝袜av在线播放| 中文字日产幕码三区国产| 国产精品自在线拍国产手机版| 欧美中文字幕在线看| 在线看亚洲一区二区三区| 一本色道久久88亚洲精品综合| jizz国产精品免费麻豆| 免费一区二区三区av| 亚洲av无码一区二区三区鸳鸯影院| www插插插无码免费视频网站|