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

自由造型的寬速域二元進氣道優化設計

2022-12-25 08:21:48王健磊龔春林
火箭推進 2022年6期
關鍵詞:優化設計

王健磊,牟 桓,龔春林

(西北工業大學 航天學院 陜西省空天飛行器設計重點實驗室,陜西 西安 710072)

0 引言

自高超聲速飛行的概念提出以來,以超燃沖壓發動機為動力的吸氣式高超聲速飛行器的研究日益受到關注[1]。相比于火箭發動機,吸氣式發動機自身不用攜帶氧化劑,可以節約推進劑的質量,提供更多的動力[2]。進氣道為吸氣式發動機捕獲足夠的流量并提供相應的壓縮,其性能的優劣對于超燃沖壓發動機以及整個飛行器的正常工作都是至關重要的。二元進氣道是一種適用于吸氣式高超聲速飛行器的進氣道。

當前,二元進氣道多為多級壓縮的形式,以飛行器前體的下表面作為進氣道的壓縮面,根據等強度激波理論求出各級折轉角的大小使其滿足飛行要求。由于高超聲速飛行器的飛行空域廣,飛行速域大,為了保證二元進氣道可以滿足整個飛行包線的要求,需要進一步開展優化設計[3-4]。國內外學者對此進行了諸多研究。

Smart以最大總壓恢復系數為目標函數,對具有2~5個壓縮波的二元進氣道進行優化設計[5]。飛行器的飛行馬赫數為2~10,以燃燒室的流動平行于來流和總的靜壓比為0.7作為限制條件。但是其沒有考慮魯棒性和燃燒相關的限制,比如最小的燃燒溫度等相關因素的影響。Shukle等同樣以最大的總壓恢復系數為目標函數對二元進氣道進行優化,但是以出口流場具有較高的均勻性作為限制條件[6]。Markell以最大總壓恢復系數為目標函數,在二元進氣道的優化中考慮了非設計狀態下的情況,其結果修正了邊界層的厚度,但是沒有考慮激波-邊界層的干擾[7]。

國內研究者也對二元進氣道的設計開展了諸多研究。范曉檣等對二元進氣道的幾何參數進行制約關系分析和參數敏感性分析,建立了二元進氣道的參數化設計方法,為二元進氣道的優化設計提供了便捷的變參數路徑[8]。張曉嘉等研究了典型二元進氣道的設計和性能估算方法,給出了設計原則,在滿足進氣道設計要求的條件下,提出了對進氣道進口、外壓波系、內壓縮通道、唇罩及隔離段的快速設計方法[9]。王向轉等研究了二元進氣道的優化設計方法,用擬牛頓法作為優化方法、代數法生成結構化網格以及流體計算軟件,開發了二元高超聲速進氣道優化設計軟件[10]。孫菲等提出了一種高超聲速二元進氣道的參數化方法和性能計算方法,針對多目標優化問題,以設計點的總壓恢復系數和喉道馬赫數為優化目標,流量系數為約束條件,用遺傳算法進行優化[11]。

在這些研究中,所選擇的目標函數只能衡量進氣道在某一狀態下的性能,難以準確地反映整個飛行速域的要求。同時,多級壓縮的二元進氣道優化設計時,需要選擇合適的設計點設計出基準構型,再對基準構型開展優化設計。設計點和總壓縮角的選擇根據設計人員的主觀意愿選取,帶來較多的不確定性。

本文提出一種新的通用的目標函數,綜合考慮不同馬赫數下的性能參數,使其可以根據飛行要求覆蓋整個飛行速域;同時提出一種不考慮設計點的設計方法,直接采用優化設計,對進氣道壓縮面進行優化,壓縮面控制方程不同,優化變量的個數不同,壓縮面的自由度更大。

1 研究流程

本文中二元進氣道的設計均采用最佳波系理論,使壓縮面在設計狀態下產生的激波可以在唇口處封口,既保證了進氣道的壓縮效率又可以給燃燒室提供足夠的氣流。利用最佳波系理論設計二元進氣道時,需要首先選定設計點和總壓縮角以及壓縮面的數量,再根據斜激波關系式求出各級壓縮角的大小。

總壓縮角的大小和壓縮面的數量會直接影響到進氣道的性能,但是目前的研究中尚未提出一種合適的方法選擇總壓縮角和壓縮面的個數;設計點的選取會直接影響到進氣道能否滿足寬速域的工作要求。由于高速狀態下飛行高度較大,空氣稀薄,對進氣道的性能要求更高,因此設計時經常選擇高馬赫數作為設計點,但由此所設計的進氣道往往無法滿足進氣道在低速狀態下的工作要求。同時,最佳波系理論基于無黏流假設,在實際飛行過程中,由于黏性的影響會出現附面層,附面層對壓縮面產生的激波進行干擾,使其向唇口下方偏移,無法達到理想的設計狀態。在之前的研究中,需要進一步的優化設計和變幾何設計使進氣道滿足寬速域的要求[12]。

本文提出直接對二元進氣道的壓縮面進行優化的設計方法,初始壓縮面為一條任意曲線,通過形狀函數來控制曲線的形狀,控制參數作為優化變量。不同狀態下的性能指標作為約束條件,選擇合適的優化算法對模型進行優化,使優化后的構型可以滿足設計要求,該方法稱為自由造型的二元進氣道設計方法。

圖1是自由造型的二元進氣道設計流程圖。從圖中可以看出,首先要對二元進氣道進行參數化建模,主要是對壓縮面進行參數化。初始壓縮面為一條直線,選擇合適的曲線方程或控制點作為壓縮面的構型,將曲線的控制參數作為優化變量;根據飛行要求,選擇合適的狀態點,利用CFD計算不同狀態下的性能參數和目標函數,當目標函數收斂時結束優化,得到進氣道的構型。

圖1 自由造型的二元進氣道設計流程圖

2 目標函數

寬速域二元進氣道的優化設計中,需要考慮進氣道在不同飛行狀態下的諸多性能參數,包括流量系數φ、總壓恢復系數σ、壓升比π等。因此,其目標函數是多點多目標。對于多點多目標優化,如何權衡每一個馬赫數下的性能在總的飛行過程中所占據的比重,尚未有研究進行明確的說明。本文提出一種新的通用的目標函數,綜合考慮不同馬赫數下的性能參數,即

(1)

式中ωi是與馬赫數有關的權重系數。根據發動機所需流量與進氣道流量之間的匹配關系,建立了式(1)及權重系數ωi的求解方法。

發動機不僅需要足夠質量流量的空氣,以保證燃燒室內的燃料充分燃燒,還要求進入燃燒室的氣流具有一定的溫度及壓強。通常采用換算流量來表示這一特性,即將空氣的質量流量與其具有的總壓和總溫組合成一個總體量[13]。換算流量的定義式為

(2)

(3)

式中:K在空氣中為常數;p0為該截面位置的總壓;T0為該截面位置的總溫;A為該截面面積;q(Ma)為該截面處的流量函數,其表達式為

(4)

(5)

(6)

式(6)中,隨來流馬赫數改變的量有σ、φ、q(Ma∞)。

發動機推力F表達式如式(7)所示,可見發動機推力直接與進入發動機的空氣流量有關,即與進氣道的性能參數φ、σ有關。

(7)

(8)

(9)

式中c為當地聲速。

3 設計方法

3.1 參數化建模

自由造型的二元進氣道的幾何參數與多級壓縮的二元進氣道幾何參數基本一致,區別在于壓縮面的控制參數不同。自由造型的二元進氣道構型如圖2所示,基本幾何參數包括進氣道長度L、唇口距前緣長度Lf、隔離段長度Lin、進氣道迎風高度hc和喉道高度ht。AB段是進氣道的壓縮面,是進行優化設計的區域。

圖2 自由造型的二元進氣道構型

針對二元進氣道的基礎構型進行設計,需要先確定進氣道的基本幾何參數的大小,再對進氣道的壓縮面進行優化設計[14]。

根據二元進氣道的設計要求,進氣道的長度L可由前體的長度給定;進氣道的迎風高度hc決定了進氣道的捕獲面積,由發動機需求確定;喉道高度ht可以通過內收縮比的大小來計算[15];唇口距前緣長度Lf無法直接計算得出,可以先計算隔離段的長度Lin,再由Lf=L-Lin求出Lf。

隔離段的長度Lin影響了進氣道的抗反壓能力,增大隔離段的長度可以增強進氣道的抗反壓能力,但是過長的隔離段會增加飛行器的結構質量。一般用隔離段的長高比(Lin/hc)來衡量隔離段抗反壓能力的大小。相關研究發現,隔離段長高比與最大承受反壓之間存在著一定的關系,當隔離段的長度增加到一定值以后,繼續增加隔離段長度不會提高抗反壓能力[16]。根據所需抗反壓能力的大小,即可求出隔離段的長度Lin。

進氣道的壓縮面AB段是對氣流進行壓縮,是進行參數化的區域。在基本幾何參數確定后,可以在坐標系中確定A、B點的坐標。將壓縮面看做是經過兩點的曲線,選擇合適的曲線方程進行擬合,將曲線方程中的控制參數作為優化變量進行優化。本文采用類別形狀函數法對壓縮段進行參數化。

3.2 類別形函數法

當前,對自由造型的二元進氣道參數化方法研究較少,曲線的參數化方法集中在飛行器的氣動外形優化中,不同的參數化方式會影響到模型精度和優化效率。常用的參數化方法為了表達復雜的曲線往往需要增加控制參數,在優化過程中會出現不規則曲線,對優化過程的穩定性產生影響。Kulfan等提出了類別形狀函數法(class and shape transformation,CST)對飛行器的氣動外形進行參數化[17]。該方法將類別函數和形狀函數組合表示幾何外形,減少了設計變量,具有良好的可控性和表達精度[18]。自CST參數化方法提出以來,在很多方面都得到應用,包括乘波體的設計[19]、翼型的減阻設計[20]、增升裝置的設計[21]等。這些研究主要是對飛行器氣動外形的設計,還沒有應用到進氣道的優化設計中。本文選擇CST方法對二元進氣道的壓縮面進行參數化建模,簡化了進氣道的設計流程,為后續進氣道的設計提供了一種新的思路。

(10)

式中:ζ=y/L;ψ=x/L。其中:L為首末點的直線距離;x為橫軸坐標;y為縱軸坐標。

(11)

式中N1和N2的取值決定了曲線幾何外形的類別。

圖3是N1、N2取不同值時類別函數代表的幾何外形。

圖3 不同類別函數代表的幾何外形

由圖3可知,該類別函數曲線與二元進氣道的形狀大致相符。N1決定了曲線頭部的形狀,其值越小,頭部曲率越大;N2決定了曲線尾部的形狀,其值越小,尾部曲率越大。

形狀函數S(ψ)可以對類別函數形成的曲線基本形狀進行修正,最終生成設計過程所需要的曲線形狀。常用的曲線形狀有B樣條曲線(B-spline curve)、貝塞爾曲線(Bezier curve)等。其中,B樣條曲線在所選擇的系數振蕩劇烈時,會偏離曲線的基本形狀,而貝塞爾曲線仍然較為光滑[20]。因此,本文選擇貝塞爾曲線的方程作為形狀函數。

貝塞爾曲線是以空間逼近為基礎的曲線,一般通過Bernstein多項式得到。可以使用n階Bernstein多項式的加權和來表示形狀函數S(ψ),表達式為

(12)

在已知曲線的兩個端點坐標時,可以用ζ(ψ)來表示曲線方程,需要確定的參數為N1、N2、bi(i=0,1,2,…,n),可以將這些參數作為優化變量開展優化設計。

4 算例分析

4.1 設計要求

吸氣式飛行器將火箭發動機和沖壓發動機結合在一起,具有結構簡單、質量小、可靠性高的優點。發動機的工作模態可分為引射模態、亞燃模態、超燃模態和純火箭模態。其中引射模態為飛行器起飛時的狀態,此時只有火箭發動機工作,進氣道處于不啟動狀態。引射模態向亞燃模態過渡的馬赫數一般較低,以減少引射模態的工作時間,節約攜帶的燃料質量,減小飛行器質量。根據超燃沖壓發動機的總體設計要求,在Ma∞=2.5時,進氣道實現自啟動,發動機的工作狀態從引射模態過渡到亞燃模態,最大工作馬赫數為Ma∞=8.0。因此,進氣道的工作范圍為Ma∞=2.5~8.0在高超聲速飛行時,由于其飛行高度較高,空氣密度較小,需要更多的空氣流量,此時需要進氣道有較大的內收縮比;在低速時,要保證進氣道的啟動性能,要求進氣道有較小的內收縮比。根據飛行器總體方案的要求,確定了氣道設計的相關要求:

飛行高度:h=10~26 km;

飛行速度:Ma∞=2.5~8.0;

進氣道長度:L=12 m;

捕獲面積:S=2 m2;

捕獲高度:hc=2 m;

總壓恢復系數:

流量系數:

4.2 優化方法

本節中的優化變量是曲線方程的控制參數,對于不規則曲線而言,控制參數的數量較多,遺傳算法計算變量較多的優化問題時的效率低,計算代價大;梯度算法利用梯度信息計算局部最優解,對于變量較多的優化問題計算效率高;直接搜索法不需要計算函數梯度,可以直接搜索達到最優解。

本節選擇霍克-吉維斯直接搜索法(Hooke-Jeeves direct search method),該算法不需要連續的目標函數和線搜索,可以處理非連續參數空間,可以迅速收斂到局部最優解,優化過程中對目標函數的調用次數較少[22]。

4.3 優化模型

基于第3節章中的設計方法,建立優化模型。目標函數采用第2節中提出的目標函數,選擇Ma=2.5、5.0、6.0、8.0作為評估狀態點,建立優化模型如式(13)所示。

(13)

式中:CD為進氣道阻力系數;CDmax為進氣道阻力系數最大值。

4.3.1 優化變量

根據設計要求,進氣道的總長度L=12 m,迎風高度hc=2 m,喉道高度ht=500 mm。根據參考文獻[16],在Ma∞=2.5時,隔離段長高比取7就可以保證進氣道具有足夠的抗反壓能力,因此隔離段長度為3 500 mm。其基本構型如圖4所示。

圖4 二元進氣道基準構型

采用類別形狀函數法對進氣道的壓縮面進行參數化,壓縮面的長度為8 637.82 mm。形狀函數選擇3階多項式,得到類別形狀函數的表達式為

(14)

式中:N1、N2、b0、b1、b2、b3均為可控參數,可以作為優化變量;ζ=y/8 637.82;ψ=x/8 637.82。

各優化變量的取值范圍如表1所示。

表1 優化變量取值范圍

4.3.2 目標函數

根據設計要求,為了使進氣道能夠滿足整個飛行速域的要求,選擇Ma∞=2.5、 5.0、 6.0、8.0作為目標函數中的評估狀態點。

Ma∞=2.5是超燃沖壓發動機開始工作的最低速度,此時應保證進氣道可以正常啟動;Ma∞=5.0是超聲速和高超聲速的分界點,也是整個飛行速域中的中間狀態點和設計要求中需要考慮的設計狀態之一;Ma∞=6.0是吸氣式飛行器高速狀態下的中間狀態,可以作為評估狀態點;Ma∞=8.0是飛行器的最大飛行速度,此時的飛行高度較大,空氣較為稀薄,對進氣道吸入氣流的要求較高。因此選擇這4個飛行狀態作為評估狀態點,得到目標函數為

(15)

4.3.3 約束條件

對寬速域二元進氣道進行約束時,需要考慮進氣道在不同設計點下的狀態的性能滿足一定的要求,涉及的性能參數有流量系數、總壓恢復系數、壓升比。所選擇的飛行狀態點與目標函數中保持一致,便于優化過程中的計算。

Ma∞=2.5時產生的溢流較大,無法滿足發動機的需求,在優化設計時應使其流量系數和總壓恢復系數滿足設計要求,即φMa∞=2.5≥0.4,σMa∞=2.5≥0.8;同時,也要使進氣道出口的氣流品質較為均勻,具有足夠大的壓升比保證進氣道能夠正常啟動,約束πMa∞=2.5≥2。

Ma∞=5.0是飛行器飛行速域的中間點,對其評估主要是為了保證在飛行過程中進氣道均有良好的性能。總壓恢復系數是衡量進氣道氣流品質的重要參數。總壓恢復系數越大,進氣道出口氣流的能量越多,更能保證燃燒室內的燃燒充分。根據設計要求,約束總壓恢復系數σMa∞=5≥0.45。

Ma∞=6.0是吸氣式飛行器在高速狀態下的中間狀態,發動機需要進氣道提供均勻穩定的流量。根據設計要求,約束此時的流量系數φMa=6≥0.9。

Ma∞=8.0是飛行器的最大飛行速度,此時的飛行高度較大,空氣密度較小,超燃沖壓發動機既需要足夠流量的氣流保證其穩定工作,也需要氣流具有較高的能量來提高燃燒效率并保證足夠抗反壓能力。在Ma∞=6.0時已經對流量系數進行約束,根據斜激波理論,Ma∞越大,激波角越小,氣流溢流更少,因此其流量有所保證。同時,由于空氣密度較小,過大的阻力系數會嚴重影響飛行器凈推力的大小,因此,需要對飛行器的阻力系數進行約束。根據設計要求,約束σMa∞=8≥0.2,πMa∞=8≥10,CD≤0.02。

4.4 優化結果

利用Isight軟件搭建優化流程,圖5是Isight求解過程。優化變量由Data Exchanger模塊傳遞給Catia模塊和優化器,依次調用Pointwise、Fluent軟件進行數值分析計算。所選擇的湍流模型仍為k-ωSST模型,來流的邊界條件為壓力遠場。采用并行計算的方式,同時計算3個狀態點,計算結果輸出到Calculator模塊計算出性能參數。

圖5 Isight求解過程

采用Hooke-Jeeves算法對模型進行優化,相對步長(relative step size)是指算法最初尋優時,變量初值擾動的范圍,在本問題取值0.1。步長縮減因子(step size reduction factor)的值越大,對于高度非線性問題來說,收斂的可能性也就越大,但是用于函數評估的花費也會相應增大;其值越小,函數評估花費和程序運行的時間也會越小,但是不收斂的風險也會增加,在本問題取值0.5。考慮Fluent存在計算精度的數值誤差,優化時運行終止的步長(termination step size)取值0.001。優化100步之后達到收斂條件,目標函數收斂過程如圖6所示。

圖6 目標函數收斂過程

由圖6可知,目標函數I收斂于0.8附近,最小值為0.827,此時的各變量值如表2所示。

表2 優化變量的取值

對優化出來的構型進行全速域的性能分析,選取的狀態點如表3所示。使用Pointwise進行結構網格劃分,在靠近壁面處加密,如圖7所示。

圖7 氣動計算結構網格

表3 全速域內的計算狀態

采用FLUENT軟件對進氣道流場進行數值求解,湍流模型選用k-ωSST模型,使用二階迎風格式離散。來流設為理想氣體,分子黏性系數采用Sutherland公式計算,計算中選擇了壓力遠場、壓力出口邊界條件。圖8是其在典型狀態下的馬赫數云圖。

圖8 優化構型的馬赫數云圖

優化得到的構型壓縮面會產生3道激波將氣流壓縮,主要是因為在構型前半部分曲線的曲率較大,此后的曲線較為平緩。同時,曲線壓縮所產生的附面層較薄,減少了溢流損失。隨著馬赫數的增大,激波逐漸向唇口處移動,在Ma∞=6.0時,激波尚未封口,但流量系數已經達到要求,為90.1%。

在Ma∞=2.5時依然會產生較多的溢流,但是流量系數依然滿足要求,為41.2%。隔離段內的氣流可以正常通過,說明進氣道可以正常啟動。隔離段內的分離區隨著馬赫數的增大逐漸減小,出口處的氣流總體趨于均勻。Ma∞=5.0時總壓恢復系數為54.3,滿足設計要求。在Ma∞=8.0時激波可以封口,流量系數為1。出口氣流較為均勻,出口馬赫數為3.69,總壓恢復系數為21.02,均滿足設計要求。

5 結論

本文以寬速域吸氣式飛行器為研究對象,研究了自由造型的二元進氣道優化設計方法,得到以下結論。

1)基于進氣道和發動機的流量匹配關系,采用歸一化方法,提出了一種新的目標函數,將多點多目標問題轉化為單一問題。對不同馬赫數下的φ/σ進行加權求和,當其最小時,發動機產生的推力最大。權重系數與馬赫數有關,選擇不同的馬赫數作為評估狀態點,可以使目標函數覆蓋整個飛行速域。

2)提出了一種針對寬速域吸氣式飛行器二元進氣道設計的新型的設計方法,建立了研究流程。該方法直接對二元進氣道的壓縮面開展優化設計,使進氣道的壓縮面有更大的自由度,為寬速域吸氣式飛行器的二元進氣道設計提供了新的研究思路。

3)研究了自由造型的二元進氣道的參數化建模方法。首先根據設計要求確定二元進氣道的主要設計參數,將壓縮段看做一條直線,利用類別形狀函數法進行參數化建模,得到二元進氣道的優化模型。

4)對自由造型的二元進氣道開展優化設計并分析其性能。使用本文提出的目標函數,以類別形狀函數的控制參數為控制變量,采用Hooke-Jeeves算法對其開展優化設計,得到二元進氣道的構型。通過數值分析發現,二元進氣道在Ma∞=2.5時的流量系數為41.2;Ma∞=5.0時總壓恢復系數為54.3;Ma∞=6.0時流量系數為90.1;Ma∞=8.0時總壓恢復系數為21.02,均滿足設計要求。

猜你喜歡
優化設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
主站蜘蛛池模板: 免费A级毛片无码免费视频| 亚洲精品片911| 免费国产高清视频| 色AV色 综合网站| 高潮毛片免费观看| 亚洲欧美在线综合图区| 色噜噜狠狠狠综合曰曰曰| 国产精品无码制服丝袜| 中国一级特黄大片在线观看| 国产网站免费观看| 亚洲欧美在线综合一区二区三区| 69av免费视频| 在线观看精品自拍视频| 国产精品人人做人人爽人人添| 中文字幕日韩丝袜一区| 亚洲欧美在线综合一区二区三区 | 欧美怡红院视频一区二区三区| 一级毛片在线播放免费观看| 日韩区欧美区| 日韩毛片免费观看| 国产精品手机视频| 尤物在线观看乱码| 日韩精品资源| 国产精品女在线观看| 91精品人妻一区二区| 四虎影视8848永久精品| 1024你懂的国产精品| 国产成人亚洲欧美激情| www中文字幕在线观看| 国产美女在线观看| 国产麻豆精品久久一二三| 99热这里只有免费国产精品 | 97国产在线视频| 女人毛片a级大学毛片免费| 国产99精品久久| 中文字幕欧美日韩| 成人日韩欧美| 午夜日本永久乱码免费播放片| 国产成人精品日本亚洲| 青草娱乐极品免费视频| 麻豆国产原创视频在线播放| 色爽网免费视频| 国产男女免费视频| 亚洲欧美日韩视频一区| 国产激爽大片在线播放| 国产成人精品一区二区秒拍1o| 国产亚洲欧美在线中文bt天堂| 亚洲中文字幕手机在线第一页| 看看一级毛片| 亚洲色图欧美在线| 免费人成网站在线高清| aaa国产一级毛片| 国产在线观看第二页| 亚洲日本中文字幕天堂网| 亚洲精品无码日韩国产不卡| 国产精品永久在线| 免费国产不卡午夜福在线观看| 中文字幕在线欧美| 亚洲成AV人手机在线观看网站| 91九色最新地址| 亚洲精品手机在线| 丁香五月亚洲综合在线| 日本免费高清一区| 精品精品国产高清A毛片| 在线观看精品国产入口| a级毛片免费看| 亚洲欧美人成电影在线观看| 99热国产这里只有精品无卡顿"| 中文国产成人精品久久| 天天综合天天综合| 欧美激情二区三区| 人妻丰满熟妇AV无码区| 国产乱人乱偷精品视频a人人澡| 精品三级在线| 丁香综合在线| 欧美成人A视频| 在线国产你懂的| 国产a网站| 精品综合久久久久久97超人该| 国产成人你懂的在线观看| 99er这里只有精品| 日韩视频精品在线|