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

可靠性數值模擬的不確定度量化

2014-06-09 12:33:44馬智博李海杰殷建偉黃文斌
計算物理 2014年4期
關鍵詞:物理模型

馬智博, 李海杰, 殷建偉, 黃文斌

(1.北京應用物理與計算數學研究所,北京 100088;2.長沙機電產品研究開發中心,湖南長沙 410100;3.流體物理研究所,四川綿陽 621900)

可靠性數值模擬的不確定度量化

馬智博1, 李海杰2, 殷建偉1, 黃文斌3

(1.北京應用物理與計算數學研究所,北京 100088;2.長沙機電產品研究開發中心,湖南長沙 410100;3.流體物理研究所,四川綿陽 621900)

根據可靠性認證對不確定度量化的技術要求,結合復雜工程系統的層級結構和數值模擬從校準、驗證與確認到最終具有預測能力的歷程,對可靠性數值模擬不確定度量化的要求和方法進行探索,并結合爆轟算例對這些方法進行演示和驗證.

數值模擬;不確定度量化;驗證與確認;可靠性認證

0 引言

可靠性問題源自產品本身和任務環境中A、B兩類不確定度的存在,其中前者對應于客觀隨機性引起的偶然不確定度(aleatory uncertainty),后者對應于主觀認知缺陷帶來的認知不確定度(epistemic uncertainty).根據對不確定問題的存在是否預知,認知不確定度又分為已知的未知(known unknowns)和未知的未知(unknown unknowns)[1].為了保證產品工作可靠,必須在綜合考慮兩類不確定度的前提下,為定型的批量產品留下足夠的設計裕度.

實物試驗和統計分析可用于量化這兩類不確定度,也可用于消減認知不確定度并發現其中未知的未知,但偶然不確定度只能通過改進設計方案和制造工藝來實現其消減.

數值模擬包含建模和模擬(modeling and simulation,M&S)兩個技術環節,建模和模擬均受人類主觀意志的主導,只能在人的認識基礎上借助計算機完成對客觀現象的分析,所以數值模擬主要用于研究偶然不確定度以及認知不確定度中已知的未知,一般不直接用于研究認知不確定度中未知的未知.

當受到成本因素或其它條件制約時,很多大型復雜產品往往難以實施系統級可靠性試驗,此時,需要用M&S技術對可靠性特征量進行工程預測,而基于裕度和不確定度量化的QMU認證方法(quantification of margins and uncertainties,QMU)為M&S在可靠性認證中發揮巨大作用提供了科學的技術平臺,M&S預測結果以及預測結果的不確定度將是QMU認證的重要信息來源[2-4].

可靠性數值模擬的任務主要表現在兩個方面:其一是研究產品狀態和任務環境與產品性能之間復雜而又確定的物理關系(對隨機過程而言,統計意義下的物理關系是確定的),發現失效模式并量化可靠性特征量,其二是研究產品狀態和任務環境的不確定度所引起的產品性能的不確定度,促成不確定度的傳遞分析.因為后者的不確定度起源于工程因素而非數值模擬本身,所以數值模擬的不確定度主要體現在反映前者對應的確定性物理關系時所帶來的模擬結果的不確定度.

對M&S實施不確定度量化(uncertainty quantification,UQ),本身也是一項復雜的系統工程,相關技術正處在發展中[5-7].本文根據可靠性認證的要求,對M&S不確定度的基本概念、組成結構、量化原則、量化方法和方法檢驗進行探索,并以爆轟計算為例展示這些方法的思想基礎和應用效果.

1 數值模擬不確定度量化的基本原則

不確定度量化的源頭信息往往是誤差數據,綜合不同行業對誤差和不確定度的定義[8-13],可以從廣義的角度上總結出兩者的普適性內涵.本文認為,誤差是行為結果與期望目標之間可以辨識的差別,它來自有期望目標的行為,沒有期望目標的行為不產生誤差,不能辨識的差別形不成誤差,而不確定度是行為結果的分散性或背離期望目標的程度,對無期望目標的行為,不確定度主要表現為行為結果的分散性,對有期望目標的行為,不確定度主要表現為行為結果背離期望目標的程度.

可靠性數值模擬的期望目標是準確預測客觀物理過程中的確定性規律,所以既產生誤差,也產生不確定度,其中誤差為模擬結果減去參照解所得的差值,基于可靠性工程的需要,不確定度應被理解為模擬結果背離真值的程度.

無論可靠性認證(針對抽象的型號)還是可靠性評估(針對具體的產品),都需要對可靠性特征量的真值作出估計,根據QMU方法對應的數據格式,該估計結果包含兩個部分,其一是最佳估計結果,其二是針對該最佳估計的不確定度估計結果,兩者共同形成關于可靠性特征量估計的不確定度區間.

再考慮到其它方面的因素,數值模擬的不確定度量化應該遵循以下原則:

1)UQ的前提是M&S的技術狀態得到相對固化且計算過程不依賴人為干預;

2)UQ的目標是量化對應于工程預測且無參照解時數值模擬的不確定度;

3)可靠性數值模擬的不確定度應該理解為模擬結果背離真值的程度;

4)模擬結果和不確定度估計值所形成的不確定區間要立足于覆蓋真值;

5)基于給定輸入信息,方法上要盡量使不確定度的量化結果達到最小;

6)量化方法要有利于診斷和消除M&S存在的問題以促進不確定度消減;

7)量化方法要有利于保證量化信息的充分性和量化結果的可信性.

2 數值模擬不確定度的組成結構

對大部分工程問題,M&S的本質是求解物理系統初始狀態和邊界條件與其最終性能的映射關系,該映射關系可表示為一個泛函[14]

其中x為物質點的空間坐標;t為時間變量,下標“0”和“e”分別表示初始和最終時刻;ξ為物質點的物理狀態參數組,反映材料的物理狀態;ω為隨空間和時間而變化的邊界條件函數組,反映產品的工作環境;η為由物質行為屬性決定的描述物質行為規律的函數組,是物理狀態ξ的函數;ζ為隨時間變化的系統性能參數組;P為由普遍物理規律決定的聯結系統性能和眾多參數或函數的泛函關系式.

建模與模擬的不確定度,指的是用模擬計算結果來預測對應實體模型的真實行為所帶來的不確定度,根據式(1),該不確定度可分解為三個組成部分:U1為實體建模不確定度U01所引起的模擬結果不確定度,實體建模即構造ξ(x,t0)和ω(x,t);U2為物理建模不確定度U02所引起的模擬結果不確定度,物理建模即構造η(ξ)和P[·];U3為數值計算所引起的模擬結果不確定度,數值計算即用數值方法求解泛函P[·].

因此,M&S的總不確定度可表示為

其中符號“⊕”表示總不確定度由三個部分組成,在某些情況下,總不確定度可近似表示為右端三項簡單相加所得的總和.

對建模引起的不確定度U1和U2,需要首先量化建模的不確定度U01和U02,然后通過傳遞分析得到U1和U2,U01和U02均由兩部分構成,其一是模型形式的不確定度Uform,其二是模型輸入參數的不確定度Upara.數值計算引起的不確定度U3,主要源自數值離散、迭代計算和四舍五入等.

3 不確定度量化的技術環節

數值模擬不確定度量化所涉及到的技術環節及其對應方法主要有以下幾個方面.

1)物理系統的層級劃分

可根據物理過程從簡單到復雜的邏輯關系或從開始到結束的先后順序,將物理系統劃分為不同的層級,層級劃分有助于隔離和診斷M&S存在的問題,也有利于獲取更多的不確定度信息.

2)建模的不確定度量化

主要目的是量化U01和U02.對不同特性的模型,可基于相應的測量技術、科學理論和實驗統計方法來獲取U01和U02的量化信息.

3)不確定度的傳遞分析

常用于對建模的不確定度U01和U02進行傳遞以得到U1和U2,該過程可借助數值模擬程序用變動分析的方法來完成,兩者分別用變分簡記為

4)模擬結果與參照解的對比分析

該方法主要基于數值模擬的驗證與確認(verification&validation,V&V)技術,通過與高可信度參照解的對比來量化數值模擬的不確定度[1],基準模型的設計及其參照解的構造是V&V的重要技術支撐.

5)不確定度的預測推斷

數值模擬的實際價值體現在對新模型的工程預測,但此時的不確定度不能通過對比分析來獲取,需要基于已有的對比信息以及對不確定度隨模型參數變化規律的認識來綜合推斷預測結果的不確定度.

6)不確定度信息的整合

對一個新系統進行數值模擬,其不確定度信息至少有兩個來源,其一是建模不確定度傳遞到系統級模擬結果上的不確定度信息,其二是基于各層級對比分析和預測推斷得到的系統級模擬結果不確定度信息,另外還可能有不同程序之間計算結果的對比信息,這些信息可通過廣義信息理論進行整合.

4 應用域內不確定度的預測推斷

根據文獻[1],可將實體模型設計參數(包括環境參數)所對應的多維參數空間劃分為確認域和應用域,確認域即已經實施的物理試驗所對應的參數空間,應用域即數值模擬預期應用所對應的參數空間,由于應用域內一般沒有試驗數據用于對比分析,所以其不確定度信息需要根據確認域內的不確定度以及兩個區域之間不確定度的某種聯系來進行推斷.

為了揭示確認域和應用域之間不確定度的聯系方式,文獻[14]提出了不確定度與實體模型的設計參數存在函數關系的假設,并結合有理論解的黎曼問題對基于該假設的應用域不確定度預測量化方法進行了驗證分析,本文擬對該方法的思想基礎做進一步完善,并結合試驗模型對該方法的有效性進行檢驗.

本文就應用域內不確定度預測推斷的科學依據提出如下觀點:

1)用規定的工具和方法對物理量(如尺寸,質量,密度,溫度,速度等)進行測量,測量結果的不確定度是被測量的函數;

2)用規定的運算法則由觀測量推算其它導出量,導出量的不確定度是觀測量的函數;

3)計算程序固化后,計算環節引起的不確定度是實體模型設計參數的函數;

4)基于以上觀點,數值模擬不確定度可表示為實體模型設計參數的函數:

5)該函數既適用于確認域也適用于應用域.

根據以上五點,如果在確認域獲取了不確定度隨實體模型設計參數的變化規律,就可將該函數內推或外推以獲得應用域內M&S的不確定度.

5 數值模擬不確定度量化的步驟

根據數值模擬不確定度量化的原則和V&V的技術思想,建議按如下步驟實施數值模擬的不確定度量化:

1)根據可靠性認證的工程需求,確認數值模擬的不確定度要求;

2)根據物理過程從簡單到復雜的邏輯關系或從開始到結束的先后順序,將物理系統劃分為不同的層級,低層級的模型有利于診斷和解決M&S中的問題以及不確定度信息的獲取;

3)根據不同層級中物理現象的特點和耦合程度,分析對應的算法中需要重點考核的技術要素,基于這些考核要素為代碼驗證逐層建立基準模型并獲取基準模型的高可信度參照解(比如理論解、人為解或由其它程序得到的被公認為比較可信的數值解);

4)結合基準模型開展代碼驗證(code verification),發現并隔離算法設計和代碼編制中存在的錯誤,驗證離散方程的相容性以及數值求解的收斂性和穩定性,保證計算程序能夠忠實地表現數學模型.代碼驗證是解算驗證的前提與基礎,其中計算誤差來自基準模型數值解與參照解的對比,可通過對該誤差進行統計分析和區間估計以實現對計算不確定度的評估;

5)根據M&S的應用需求在不同層級上開展解算驗證(solution verification),以驗證數值計算是否滿足給定的準確度要求.解算驗證也要估計輸入數據和后處理結果的不確定度,它是確認活動的前提和技術基礎,所對應的模型較復雜,一般沒有用于對比分析的參照解,可通過數值解隨離散尺度的收斂趨勢來估計其收斂解——如Richardson外推方法等,并借助帶安全系數的網格收斂指標(grid convergence index,GCI)來估計數值計算的不確定度,并使該結果滿足不確定度量化的真值覆蓋原則;

6)針對所需的實體建模和物理建模,設計并實施各個層級的校準試驗,對建模和模擬的模型參數、計算參數和可調參數進行校準.可調參數指的是無物理意義或雖然有物理意義但存在較大認知缺陷的對M&S結果有明顯影響的建模或計算參數;

7)開展各個層級的確認試驗,用校準后的模型參數進行數值模擬,通過數值模擬結果與試驗結果的對比信息,量化各層級數值模擬的不確定度,并根據確認過程發現的問題,完成對建模和模擬的技術改進,以縮減數值模擬不確定度;

8)通過預測推斷由系統級確認域內的不確定度獲得應用域內的不確定度;

9)通過傳遞分析獲得由建模不確定度所引起的系統級數值模擬不確定度;

10)將預測推斷和傳遞分析得到的不確定度信息進行整合,實現應用域內數值模擬總體不確定度的量化.

6 實例演示

6.1 實體模型

該實體模型主要包含一個起爆器和一塊炸藥,在軸線上炸藥的厚度為δ,該模型及其試驗裝置主要研究炸藥起爆初期的爆轟波發展過程,用以建立該階段的起爆和傳爆物理模型.測量的參數為爆轟波傳到炸藥右端面時在軸線上貼近LiF窗口的炸藥峰值速度,該速度的相對測量誤差為±3.0%,如圖1所示.

圖1 實體模型及試驗設備Fig.1 The entitymodel and test facility

6.2 層級劃分

將該實體模型視為一個系統,并劃分為兩個層級:

1)起爆器層級,僅有起爆器,沒有炸藥件;

2)全系統層級,包含起爆器和炸藥件.

起爆器層級要先于系統層級開展校準和確認活動,用于建立和評價起爆器動作過程的數值模擬可信度,在此基礎上,再開展全系統層級的校準和確認活動.

6.3 模型校準

進行模型確認之前,要先進行模型校準,即通過試驗對物理模型的參數進行校準,在校準之前,要首先選定物理模型的形式,模型形式不在此細述.該模型的校準活動分兩步:

1)起爆器層級校準

通過起爆器外半球形金屬薄膜中心位置的峰值速度測量值來校準起爆器內數值模擬的物理模型和計算參數;

2)全系統層級校準

選擇不同的炸藥厚度,通過炸藥右端面峰值速度測量值來校準炸藥對應起爆和傳爆過程的物理模型和計算參數.

6.4 模型確認

該實體模型只有一個設計變量,即炸藥的厚度δ,炸藥密度不作為設計變量.

該例主要關心炸藥峰值速度V的數值模擬不確定度隨炸藥厚度的變化規律,炸藥厚度δ和炸藥密度均具有加工隨機性,數值模擬時要輸入它們的實測值.

確認試驗分4組,δ的設計值分別為5 mm、10 mm、15 mm和20 mm,每組共有2次重復性試驗,由于加工隨機性,即使重復性試驗,炸藥樣品的厚度和密度也不會嚴格等于其設計值.

在確認環節,固定校準后的模型參數(比如JWL狀態方程中的5個模型參數)不變,分別獲得8個實體模型的峰值速度模擬結果,如表1所示,對應的試驗結果如表2所示.

表1 確認域內峰值速度V的模擬結果(δ/mm,V/(m·s-1))Table 1 Simulation results of V in the validated domain(δ/mm,V/(m·s-1))

表2 確認域內峰值速度V的試驗結果(δ/mm,V/(m·s-1))Table 2 Test results of V in validated domain(δ/mm,V/(m·s-1))

本文將模擬誤差定義為模擬結果減去試驗結果:

由于每組的重復試驗只有兩次,我們根據這兩次試驗中模擬結果對試驗結果的最大背離來近似提取確認域內建模與模擬的不確定度信息,如表3所示.

6.5 應用域內的不確定度量化

確認域內炸藥厚度δ的設計范圍是5 mm~20 mm,本文將炸藥厚度δ分別為2 mm、3 mm、4 mm和25 mm的實體模型劃入應用域內,根據δ的取值范圍,該例的確認域和應用域不重合.

在現實情況下,應用域內沒有試驗值可以參照,不可能通過與試驗結果的對比來量化M&S不確定度,必須建立不確定度的預測模型.

基于式(5)所表達的基本思想,擬用二次多項式構造不確定度預測模型,獲取不確定度隨炸藥厚度δ的變化規律

表3 確認域內峰值速度V的模擬誤差和不確定度(δ/mm,EV/(m·s-1),UV/(m·s-1))Table 3 Simulation errors&uncertainties of V in validated domain(δ/mm,EV/(m·s-1),UV/(m·s-1))

方程(8)一般為超定方程,根據最小二范數原則進行求解,得到

至此,式(7)所對應的函數形式被完全確定下來,并由此得到應用域內的不確定度預測結果,如表 4所示.

表4 應用域內峰值速度V的M&S不確定度預測值(δ/mm,UV/(m·s-1))Table 4 Predicted M&S uncertainties of V in the applied domain(δ/mm,UV/(m·s-1))

7 對不確定度量化方法的檢驗

本文僅對基于預測推斷的應用域內不確定度量化方法進行檢驗.

基于上述算例模型,本文安排以下步驟進行檢驗:

1)根據預測模型,推斷應用域內對應不同炸藥厚度δ的M&S不確定度,如表4所示;

2)對應用域內每個δ設計值,加工2個試驗模型,共8個試驗模型,制定試驗計劃;

3)對8個試驗模型進行初態參數的測量,在試驗前完成數值模擬并凍結模擬結果;

4)完成8個試驗模型的試驗,記錄試驗的測試結果;

5)將模擬結果減去測試結果,得到8次試驗的模擬誤差;

6)對比模擬誤差和推斷結果以檢驗應用域內不確定度量化方法的準確性.

對比結果列于表5.將表5中每個δ設計值下兩次重復試驗的最大誤差絕對值與表4中由預測模型得到的不確定度進行對比,發現預測模型能基本上反映出不確定度隨δ的變化規律,具有一定的預測推斷效果.由于試驗數太少,統計結果的隨機性較大,應用域內的量化結果還不能完全達到真值覆蓋的要求,并且隨著設計點遠離確認域,預測的準確性也會逐漸降低.

表5 應用域內峰值速度V的模擬結果和模擬誤差(δ/mm,V/(m·s-1))Table 5 Simulation and test results and errors of V in applied domain(δ/mm,V/(m·s-1))

8 結論

在缺乏試驗的前提下,全系統級的可靠性認證主要依靠數值模擬結果以及該結果對應的不確定度,對不確定度過大或過小的估計將顯著增加產品認證決策中“棄真”和“納偽”的錯誤,為供需雙方帶來損失.本文根據復雜工程系統的可靠性認證將越來越多地依賴數值模擬這一發展趨勢,對可靠性數值模擬不確定度量化的概念、原則、思想方法與實施步驟進行了分析和歸納,重點對數值模擬用于預測任務時應用域內的不確定度量化方法進行了探索,并結合算例演示了不確定度量化所經歷的層級劃分、模型校準、驗證與確認、應用域內不確定度量化及其方法檢驗等過程.檢驗結果表明,以實體模型參數為自變量而建立的不確定度函數,能夠基本上反映數值模擬不確定度在該參數空間中的變化規律,并能夠用于應用域內數值模擬不確定度的推斷,但隨著實體模型遠離確認域,推斷的準確性會逐漸降低,另外,由于該方法的思想基礎源自對科學規律的主觀認識,所以還需要有更多的實踐檢驗.

[1]OberkampfWL,Roy C J.Verification and validation in science computing[M].Cambridge University Press,2010:1-767.

[2]Ma Zhibo,Ying Yangjun,Zhu Jianshi.QMU certifyingmethod and its implementation[J].Chinese Journal of Nuclear Science and Engineering,2009,29(1):1-9.

[3]Pilch M,Trucano T G,Helton J C.Ideas underlying the quantification of margins and uncertainties[J].Reliability Engineering and System Safety,2011,96:965-975.

[4]Helton JC.Quantification ofmargins and uncertainties:Conceptual and computational basis[J].Reliability Engineering and System Safety,2011,96:976-1013.

[5]Doostan A,Iaccarino G.A least-squares approximation of partial differential equations with high-dimensional random inputs[J].Journal of Computational Physics,2009,228:4332-4342.

[6]Helton JC,Johnson JD,Sallaberry C J,Storlie C B.Survey of sampling-based method for uncertainty and sensitivity analysis[J].Reliability Engineering and System Safety,2006,91(10-11):1175-1209.

[7]Roy C J,OberkampfWL.A comprehensive framework for verification,validation,and uncertainty quantification in scientific computing[J].Computer Methods in Applied Mechanics and Engineering,2011,200:2131-2144.

[8]Reston,VA,Guide for verification and validation of computational fluid dynamics[S].American Institute of Aeronautics and Astronautics,AIAA G-077-1998.

[9]Ye Depei,Song Zhenguo,Wang Xianzhi,Expression and evaluation of uncertainty inmeasurement[S].Military Standard of the People's Republic of China,1999,GJB3756-99.

[10]The American Society ofMechanical Engineers.Testuncertainty[S].American National Standard,2005,ASME PTC 19.1-2005.

[11]The American Society of Mechanical Engineers.Guide for verification and validation in computational solid mechanics[S]. American National Standard,2006,ASME V&V 10-2006.

[12]The American Society of Mechanical Engineers.Standard for verification and validation in computational fluid dynamics and heat transfer[S].American National Standard,2009,ASME V&V 20-2009.

[13]The American Society of Mechanical Engineers.An illustration of the concepts of verification and validation in computational solid mechanics[S].American National Standard,2012,ASME V&V 10.1-2012.

[14]Ma Zhibo,Zheng Miao,Yin Jianwei,et al.Quantification of uncertainties in detonation simulations[J].Chinese Journal of Computational Physics,2011,28(1):66-74.

Uncertainty Quantification of Numerical Simulation for Reliability Analysis

MA Zhibo1,LIHaijie2,YIN Jinwei1,HUANGWenbin3
(1.Institute of Applied Physics and Computational Mathematics,Beijing 100088,China;2.Electromechanical Product Research and Development Center,Changsha 410100,China;3.Institute of Fluid Physics,Mianyang 621900,China)

According to demands of uncertainty quantification for reliability certification,in hierarchy of complex engineering system and development from calibration,verification&validation to prediction ability ofmodeling&simulation(M&S),the principle and method of uncertainty quantification of M&S are investigated.An example on detonation process is shown in which the ideas are demonstrated.

numerical simulation;uncertainty quantification;verification&validation;reliability certification

date:2013-07-12;Revised date:2013-12-07

0242.1

A

1001-246X(2014)04-0424-07

2013-07-12;

2013-12-07

國家自然科學基金(11371066)、中國工程物理研究院科學技術發展基金(2013A0101004)及裝備預先研究(51319010207)資助項目

馬智博(1963-),男,河南,博士,研究員,主要從事計算力學和可靠性工程研究,E-mail:mazhibo@iapcm.ac.cn

猜你喜歡
物理模型
一半模型
只因是物理
井岡教育(2022年2期)2022-10-14 03:11:44
重要模型『一線三等角』
如何打造高效物理復習課——以“壓強”復習課為例
重尾非線性自回歸模型自加權M-估計的漸近分布
處處留心皆物理
我心中的物理
三腳插頭上的物理知識
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 99热最新网址| 国产在线97| 亚洲男人天堂久久| 久久国产精品娇妻素人| 免费人成网站在线高清| 九九久久99精品| 无遮挡一级毛片呦女视频| 九九热精品视频在线| 久久国产精品电影| 99er精品视频| a天堂视频在线| 在线观看免费黄色网址| 久久综合丝袜长腿丝袜| 伊人天堂网| 毛片卡一卡二| 午夜一区二区三区| 噜噜噜久久| 91精品伊人久久大香线蕉| 毛片在线播放a| 日本免费一级视频| 欧美不卡视频一区发布| 日本成人在线不卡视频| 青青久在线视频免费观看| 高清色本在线www| 国产91麻豆免费观看| 狠狠ⅴ日韩v欧美v天堂| 九月婷婷亚洲综合在线| 中文字幕久久波多野结衣| 亚洲国产精品国自产拍A| 中文一级毛片| 国产成人91精品| 国产在线精彩视频二区| 亚洲成人在线网| 欧美激情综合| 成人精品亚洲| 亚洲性视频网站| 国国产a国产片免费麻豆| 日本人又色又爽的视频| 欧美a在线看| 一本色道久久88| 国产高清不卡| 亚洲视屏在线观看| 亚洲中文精品久久久久久不卡| 99re热精品视频国产免费| 日韩在线中文| 国模视频一区二区| 波多野结衣在线一区二区| 免费又爽又刺激高潮网址| 激情无码字幕综合| 精品国产电影久久九九| 99在线视频免费观看| 在线精品亚洲国产| 国产自无码视频在线观看| 欧美一级夜夜爽www| 在线国产你懂的| 亚洲中文字幕无码爆乳| 日本爱爱精品一区二区| 天天做天天爱夜夜爽毛片毛片| 日本亚洲国产一区二区三区| 国产在线精彩视频二区| 亚洲丝袜中文字幕| 中文字幕在线播放不卡| 欧美日本激情| 欧美激情综合一区二区| 国产女人18毛片水真多1| 亚洲香蕉久久| 国产欧美日韩另类精彩视频| 婷婷亚洲最大| 亚洲国产精品成人久久综合影院 | 亚洲成人网在线播放| 欧美国产综合色视频| 亚洲熟女中文字幕男人总站| 国产综合精品日本亚洲777| 欧美自慰一级看片免费| 欧美无遮挡国产欧美另类| 狠狠v日韩v欧美v| 久久香蕉欧美精品| 日韩欧美中文在线| 国产成人8x视频一区二区| 久久久久久久蜜桃| 日韩第九页| 久久成人国产精品免费软件 |