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

細(xì)觀尺度下推進(jìn)劑/襯層/絕熱層粘接界面參數(shù)反演研究

2023-05-31 13:43:50肖云東王玉峰李何龍賴帥光韓永恒
兵器裝備工程學(xué)報(bào) 2023年5期
關(guān)鍵詞:界面模型

肖云東,王玉峰,李何龍,賴帥光,伍 鵬,韓永恒

(1.海軍航空大學(xué), 山東 煙臺(tái) 264001; 2. 91526部隊(duì), 廣東 湛江 524000;3. 91468部隊(duì), 海南 陵水 572400;4.海軍裝備部駐北京地區(qū)軍事代表局, 北京 100071)

0 引言

粘接界面的力學(xué)性能是影響固體火箭發(fā)動(dòng)機(jī)壽命的關(guān)鍵因素之一[1]。隨著對(duì)粘接界面損傷研究的逐漸深入,眾多學(xué)者從細(xì)觀角度出發(fā),運(yùn)用細(xì)觀均勻化理論和有限元法定量分析界面的力學(xué)性能,這有利于掌握界面破壞的本質(zhì)。

有限元法中通常采用內(nèi)聚力模型描述界面的力學(xué)性能,準(zhǔn)確的內(nèi)聚力模型是模擬界面失效過程的前提。關(guān)于如何確定內(nèi)聚力模型的問題中,大多預(yù)先定義模型種類,再通過經(jīng)驗(yàn)法、反演識(shí)別(有限元模型修正技術(shù))等方法確定模型參數(shù)。針對(duì)細(xì)觀尺度下的界面脫粘問題,由于推進(jìn)劑內(nèi)部顆粒夾雜,粘接界面處的組成復(fù)雜,采用經(jīng)驗(yàn)法則很難保證內(nèi)聚力模型的準(zhǔn)確性。

有學(xué)者采用反演識(shí)別方法,即基于仿真與實(shí)測(cè)結(jié)果構(gòu)造目標(biāo)函數(shù),通過最小化目標(biāo)函數(shù),迭代識(shí)別細(xì)觀尺度下的內(nèi)聚力模型參數(shù)。封濤[2]與職世君[3]建立了復(fù)合固體推進(jìn)劑二維細(xì)觀填充模型,均以仿真與實(shí)測(cè)的應(yīng)力-應(yīng)變數(shù)據(jù)信息構(gòu)造目標(biāo)函數(shù),分別通過Hooke-Jeeves優(yōu)化算法與分步迭代法對(duì)顆粒/基體界面之間的內(nèi)聚力模型參數(shù)進(jìn)行了反演識(shí)別。但其目標(biāo)函數(shù)信息中僅能反映出試件的力學(xué)行為,無法反映出試件的局部信息。

通過數(shù)字圖像相關(guān)技術(shù)(digital image correlation,DIC)可獲取試件全場(chǎng)二維/三維的變形信息,可實(shí)現(xiàn)位移和應(yīng)變的量化分析[4]。李高春[5]使用SEM對(duì)復(fù)合固體推進(jìn)劑拉伸變形破壞過程進(jìn)行了觀察,并通過DIC技術(shù)獲取了推進(jìn)劑細(xì)觀表面位移場(chǎng)。基于DIC技術(shù)的反演識(shí)別方法中包含了大量的材料響應(yīng)數(shù)據(jù),在統(tǒng)計(jì)學(xué)角度上,該反演識(shí)別方法的準(zhǔn)確性更高。Fedele[6]采用DIC技術(shù)與有限元法相結(jié)合的反演識(shí)別方法對(duì)眩光層壓板的脫膠實(shí)驗(yàn)結(jié)果擬合,確定了指數(shù)型內(nèi)聚力模型參數(shù)。楊思滿[7]通過DIC技術(shù)獲取試件表面的位移場(chǎng),并與仿真位移場(chǎng)構(gòu)造目標(biāo)函數(shù),采用遺傳算法和序列二次規(guī)劃法相結(jié)合的策略對(duì)材料參數(shù)進(jìn)行反求。

為提高細(xì)觀尺度下粘接界面內(nèi)聚力模型的準(zhǔn)確性,通過SEM原位拉伸實(shí)驗(yàn),記錄粘接試件在拉伸過程中的變形與破壞,根據(jù)界面附近顆粒分布情況建立了簡(jiǎn)化后的粘接界面細(xì)觀模型,采用DIC技術(shù)結(jié)合Hooke-Jeeves 優(yōu)化算法對(duì)界面所采用的雙線性內(nèi)聚力模型相關(guān)參數(shù)開展反演識(shí)別研究。

1 基于DIC技術(shù)的反演識(shí)別方法

1.1 反演思路

通過相關(guān)文獻(xiàn)[8],預(yù)先定義粘接界面所采用的模型為雙線性內(nèi)聚力模型,經(jīng)過大量的試算工作確定初始模型的相關(guān)參數(shù);采用有限元法和DIC技術(shù)分別獲取ROI的位移場(chǎng)信息;基于仿真與實(shí)測(cè)位移場(chǎng)信息構(gòu)建目標(biāo)函數(shù),并采用Hooke-Jeeves優(yōu)化算法沿著目標(biāo)函數(shù)減小的趨勢(shì)逐步修正待求參數(shù)。反演流程如圖1所示。

圖1 反演流程框圖

1.2 Hooke-Jeeves優(yōu)化算法

Hooke-Jeeves優(yōu)化算法又稱模式搜索算法,通過目標(biāo)函數(shù)來探索設(shè)計(jì)空間,可對(duì)待求的參數(shù)按照正確的“趨勢(shì)”賦值,具有適定性強(qiáng)、收斂速度快等優(yōu)點(diǎn)。算法基于所賦的初值點(diǎn)開始,由“探測(cè)移動(dòng)”沿著不同軸不同方向進(jìn)行探測(cè)性移動(dòng),從而確定目標(biāo)函數(shù)值下降的方向,再通過所設(shè)定的“模式移動(dòng)”,使目標(biāo)函數(shù)值快速沿著“山谷”方向,逐漸逼近目標(biāo)函數(shù)最小值。恰當(dāng)?shù)某跏蓟c(diǎn)可大幅度提高反演精度,降低計(jì)算成本。通過“探測(cè)移動(dòng)”與“模式移動(dòng)”多次交替運(yùn)行,直至最終符合收斂條件時(shí),計(jì)算終止,從而實(shí)現(xiàn)尋優(yōu)目的[9]。

2 細(xì)觀拉伸實(shí)驗(yàn)

制作了厚度為3 mm的微型非標(biāo)準(zhǔn)試件,如圖2所示。其兩端為有機(jī)玻璃,中間部分是通過切取矩形粘接試件得到的推進(jìn)劑/襯層/絕熱層粘接件,使用302膠將3部分黏結(jié),并自然固化2 d。

圖2 微形非標(biāo)準(zhǔn)試件及其尺寸

通過日本島津JSM-5410LV型試驗(yàn)系統(tǒng)(全數(shù)字電液伺服加載系統(tǒng)和高精度的掃描電子顯微鏡)進(jìn)行細(xì)觀拉伸實(shí)驗(yàn),并實(shí)時(shí)記錄試件在拉伸載荷作用下的力學(xué)行為與表面細(xì)觀結(jié)構(gòu)的變化情況。

拉伸過程由伺服控制器自動(dòng)加載,加載速率設(shè)定為1.2 mm·min-1。SEM放大倍數(shù)為50,調(diào)節(jié)對(duì)比度和亮度后,可清晰成像,且圖片特征明顯。拉伸位移與推進(jìn)劑/襯層/絕熱層粘接件邊長的比值定義為拉伸應(yīng)變?chǔ)?襯層附近區(qū)域的破壞過程如圖3所示。

圖3 粘接界面破壞過程(×50)

其中圖3(b)、圖3(c)是ε為8%、18%時(shí),通過開源數(shù)字圖像相關(guān)分析軟件Ncorr獲取圖3(a)虛線內(nèi)沿著x方向的實(shí)測(cè)位移(Uexpx)云圖。圖3(a)中,2條白色虛線分別位于絕熱層/襯層界面與襯層/推進(jìn)劑界面處,d0為初始襯層厚度。Δdi表示拉伸過程中的界面張開位移:

Δdi=di-d0

(1)

式中,i=1、2、3分別表示ε為8%、18%、30%的受載狀態(tài)。拉伸過程中的襯層厚度與界面張開位移情況,見表1所示。

表1 襯層厚度與界面張開位移Table 1 Interlayer thickness and opening displacement μm

3 有限元建模

3.1 細(xì)觀模型建立

推進(jìn)劑/襯層/絕熱層粘接試件的代表性體積單元(representative volume element,RVE)中包括HTPB推進(jìn)劑、HTPB/IPDI襯層、EPDM絕熱層,HTPB推進(jìn)劑中又包含了大量的固體燃料等顆粒。RVE在宏觀上尺寸需足夠小,可看成一個(gè)物質(zhì)點(diǎn),滿足連續(xù)介質(zhì)力學(xué)的基本假定;在細(xì)觀上尺寸需足夠大,包含足夠多的細(xì)觀元素和細(xì)觀結(jié)構(gòu)信息,可代表局部連續(xù)介質(zhì)的統(tǒng)計(jì)平均性質(zhì)。模型中的各部分尺寸,如表2所示。

表2 模型尺寸

其中,基體尺寸滿足推進(jìn)劑代表性體積單元最小尺寸680 μm×680 μm的要求[10]。Al顆粒粒徑一般為10~30 μm,AP顆粒粒徑總體呈雙峰分布規(guī)律,小顆粒平均粒徑約為20 μm,大顆粒平均粒徑在100~300 μm,各顆粒的粒徑相差大,若將其全部構(gòu)造,會(huì)極大程度上增加計(jì)算量,且不易收斂。因此,采用多步等效法將粒徑小于80 μm的顆粒等效為基體的一部分,其余顆粒簡(jiǎn)化為圓形或橢圓形,通過二維圖像重構(gòu)的方式[11]建立粘接界面細(xì)觀模型,如圖4所示。

圖4 粘接界面細(xì)觀模型

粘接界面的脫粘行為與顆粒/基體的脫濕行為分別通過黏結(jié)單元、黏結(jié)接觸開展仿真計(jì)算。黏結(jié)單元較為常見,不過多介紹。黏結(jié)接觸是通過從屬面與主控面之間的接觸從而定義損傷模型,由從屬面上的點(diǎn)對(duì)主控面進(jìn)行投影,從而確定接觸點(diǎn)對(duì),若對(duì)應(yīng)的接觸點(diǎn)不在網(wǎng)格節(jié)點(diǎn)處,可通過鄰近節(jié)點(diǎn)的計(jì)算結(jié)果進(jìn)行插值獲取,具有建模簡(jiǎn)單、易收斂等優(yōu)點(diǎn)[12]。

3.2 本構(gòu)關(guān)系

雙線性內(nèi)聚力模型的形式簡(jiǎn)單、且易于實(shí)現(xiàn),在界面脫粘問題中被廣泛采用[13]。粘接界面脫粘行為與顆粒/基體脫濕行為均采用雙線性內(nèi)聚力模型,如圖5所示。

圖5 雙線性內(nèi)聚力模型

n為垂直界面的法向方向,s為沿界面方向的切向方向,δ0為界面損傷起始位移,δf為界面失效位移。

當(dāng)δ≤δ0時(shí),界面未發(fā)生損傷,界面力T隨界面張開位移的增大而線性增大,其關(guān)系如式(2):

(2)

式中,K為界面剛度。

界面損傷起始準(zhǔn)則采用最大名義應(yīng)力準(zhǔn)則,認(rèn)為界面力達(dá)到界面的最大名義應(yīng)力時(shí),界面開始發(fā)生損傷:

(3)

式中:σmax、τmax為法向與切向的最大名義應(yīng)力;〈〉為Macaulay括號(hào),表示純壓應(yīng)力狀態(tài)不會(huì)引發(fā)損傷。

采用基于位移損傷演化準(zhǔn)則,當(dāng)δ0≤δ≤δf時(shí),界面發(fā)生損傷,界面力隨著界面張開位移增加而線性減小:

(4)

式中,D表征界面損傷狀態(tài),對(duì)應(yīng)仿真軟件的輸出量為SDEG。當(dāng)D=0時(shí)表明界面未損傷,D介于0~1時(shí)表明界面處于損傷階段,D=1時(shí)表明界面完全失效開裂:

(5)

式中,δm表示混合模式下的界面張開位移:

(6)

可以看出,雙線性內(nèi)聚力模型的形狀可由界面剛度、最大名義應(yīng)力、界面失效位移所決定,為了方便后續(xù)仿真計(jì)算,模型簡(jiǎn)化為各向同性,將各物理量法向與切向的數(shù)值設(shè)置相同。基于文獻(xiàn)[2]中的部分界面模型參數(shù),又經(jīng)過大量的有限元試算,最終各界面所采用的模型參數(shù),如表3所示。

表3 各界面的模型參數(shù)

將絕熱層、基體、顆粒簡(jiǎn)化為彈性材料,對(duì)絕熱層試件進(jìn)行單向拉伸實(shí)驗(yàn),取其初始模量作為彈性模量,其泊松比以及基體和顆粒的力學(xué)性能參數(shù)參考文獻(xiàn)[14],如表4所示。

表4 粘接結(jié)構(gòu)各部分力學(xué)性能參數(shù)

3.3 網(wǎng)格劃分與邊界條件

絕熱層與顆粒的楊氏模量較基體相比要大得多,載荷較小的情況下可近似認(rèn)為二者不發(fā)生變形,對(duì)其采用四結(jié)點(diǎn)雙線性平面應(yīng)變四邊形單元CPE4;基體作為近似不可壓材料,受載時(shí)體積不發(fā)生變化,采用四節(jié)點(diǎn)四邊形線性積分的雜交單元CPE4H來模擬其在大變形下的力學(xué)響應(yīng);襯層采用四節(jié)點(diǎn)二維粘接單元COH2D4。各界面附近的網(wǎng)格需適當(dāng)細(xì)化以提高計(jì)算的準(zhǔn)確性和收斂性。

粘接試件在拉伸載荷的作用下,Y方向變形小,對(duì)邊界條件進(jìn)行簡(jiǎn)化,仿真過程中不考慮試件沿Y方向的變形。將A、C邊設(shè)置為沿Y方向固定,D邊設(shè)置為完全固定,B邊沿X正方向受均布拉伸載荷,速率為1.2 mm/min。網(wǎng)格劃分與邊界條件的設(shè)置情況,如圖6所示。

圖6 網(wǎng)格劃分與邊界條件

4 參數(shù)反演

4.1 反演過程與結(jié)果

將部分襯層與推進(jìn)劑區(qū)域作為ROI,為圖3(a)藍(lán)色虛線內(nèi)的區(qū)域。文獻(xiàn)[15]中指出界面模型會(huì)對(duì)試件的力學(xué)行為產(chǎn)生影響,界面本構(gòu)關(guān)系也勢(shì)對(duì)影響推進(jìn)劑表面的變形特征,同時(shí),通過細(xì)觀拉伸實(shí)驗(yàn)發(fā)現(xiàn):界面的本構(gòu)關(guān)系主要影響ROI沿x方向的位移場(chǎng)分布,對(duì)沿y方向的位移場(chǎng)分布的影響可忽略。

以ε為8%、18%時(shí)的位移場(chǎng)(Ux)信息共同構(gòu)造目標(biāo)函數(shù),構(gòu)造方法參考文獻(xiàn)[16]并改進(jìn):

(7)

式中:p為所有界面參數(shù)組成的向量集合;k為拉伸狀態(tài);i為目標(biāo)點(diǎn)的編號(hào);m為目標(biāo)點(diǎn)的數(shù)量。

反演程序中的增量步長設(shè)為0.05、步長縮減因子設(shè)為0.5、步長加速因子設(shè)為1、增量步長閾值設(shè)為5×10-3,粘接界面的內(nèi)聚力模型參數(shù)初值,見表1所示。

迭代過程中,界面剛度、最大名義應(yīng)力、界面失效位移以及目標(biāo)函數(shù)值的變化情況,如圖7所示。

圖7 迭代過程曲線

迭代過程共計(jì)100次賦值計(jì)算,目標(biāo)函數(shù)曲線在第83次計(jì)算后,開始逐漸平穩(wěn),各曲線的收斂情況較好。在第16輪的“探測(cè)移動(dòng)”中,增量步長為3.12×10-3小于設(shè)定的增量步閾值,完成此輪“探測(cè)移動(dòng)”后,計(jì)算終止。第98次賦值試算中的目標(biāo)函數(shù)值最小,認(rèn)為此界面模型下的ROI的仿真位移與實(shí)測(cè)位移最為接近,即最終的反演結(jié)果,如表5所示。

表5 反演結(jié)果值

4.2 反演效果分析

反演獲取的模型參數(shù)的準(zhǔn)確性需通過相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證。基于反演結(jié)果開展仿真計(jì)算,將仿真獲取的ROI位移場(chǎng)與實(shí)測(cè)數(shù)據(jù)進(jìn)行比較。ε為8%、18%時(shí),ROI的仿真位移(Usimx)云圖,如圖8所示,相對(duì)應(yīng)的實(shí)測(cè)位移(Uexp x)云圖,見圖3(b)、圖3(c)。可見2種邊界條件下的位移(Ux)云圖分布特征均基本一致,相似度極高。

各節(jié)點(diǎn)處的位移誤差ηU:

(8)

ε為8%、18%時(shí),ROI的位移(Ux)誤差云圖,如圖9所示。

圖9 位移(Ux)誤差云圖

ROI沿x方向的位移誤差統(tǒng)計(jì)情況,見表6所示。

主要誤差因素有4個(gè)方面:① 界面附近存在大量微裂紋與微孔洞等缺陷,部分顆粒存在破碎與包裹不完全等缺陷,導(dǎo)致ROI的下區(qū)域平均位移(Uexpx)略高于上區(qū)域,而細(xì)觀模型中未對(duì)其考慮;② 各圖片亮度和對(duì)比度在獲取過程中很難保證完全一致,會(huì)影響圖片的灰度分布,進(jìn)而影響Ncorr計(jì)算結(jié)果;③ (1~7)顆粒附近區(qū)域的變形梯度大,Ncorr軟件中出現(xiàn)匹配相關(guān)性差、測(cè)量精度低的問題;④ 仿真過程中,將各材料簡(jiǎn)化為彈性,部分材料的力學(xué)性能參數(shù)直接引用相關(guān)文獻(xiàn),與實(shí)際力學(xué)性能存在偏差。

表6 ROI位移誤差情況(Ux)

5 界面損傷分析

采用界面單元的SDEG對(duì)界面的損傷程度進(jìn)行表征,使用Python腳本語言對(duì)仿真軟件進(jìn)行2次開發(fā),提取內(nèi)聚力界面單元高斯積分點(diǎn)的SDEG數(shù)據(jù),繪制SDEG變化曲線,如圖10所示。

圖10 損傷因子變化曲線

起初在外載荷的作用下,界面力隨著界面張開位移的增加呈線性增長,SDEG均為0,未進(jìn)入損傷階段;界面力達(dá)到最大值后,界面處損傷開始萌生并擴(kuò)展;隨著界面張開位移的逐漸增加,裂紋逐步形成擴(kuò)展,界面承受載荷的能力降低;當(dāng)界面力減小至0時(shí),界面損傷達(dá)到高峰,SDEG為1,此時(shí)界面完全損傷。

結(jié)合粘接界面破壞過程(見圖3),以及獲取的法向內(nèi)聚力模型可知:當(dāng)ε=8%時(shí),由于推進(jìn)劑的楊氏模量較小,其在載荷作用下的變形明顯,此時(shí)界面張開位移為75 μm,SDEG為0,界面力為0.48 MPa;當(dāng)界面張開位移為112.5 μm時(shí),界面力達(dá)到最大值0.72 MPa;當(dāng)ε=18%時(shí),(1~7)顆粒脫濕現(xiàn)象明顯,且6號(hào)、7號(hào)顆粒脫濕情況額外嚴(yán)重,脫濕處的承載能力降低,界面附近區(qū)域產(chǎn)生微裂紋,導(dǎo)致(Uexpx)云圖中產(chǎn)生像素點(diǎn)不匹配的無位移輸出區(qū)域,此時(shí)界面張開位移為150 μm,大于界面損傷起始位移,界面進(jìn)入損傷階段,表現(xiàn)出線性軟化行為,SDEG為0.26,界面力為0.70 MPa;當(dāng)ε=30%時(shí),界面與顆粒/基體附近的微裂紋匯合產(chǎn)生宏觀裂紋,且宏觀裂紋的擴(kuò)展方向幾乎與拉伸方向垂直,此時(shí)界面張開位移為375 μm,SDEG為0.75,界面力為0.61 MPa。

6 結(jié)論

1) 根據(jù)雙線性內(nèi)聚力模型中的參數(shù)特性,采用基于DIC技術(shù)與Hooke-Jeeves優(yōu)化算法相結(jié)合的反演識(shí)別方法,實(shí)現(xiàn)了對(duì)細(xì)觀尺度下粘接面損傷模型(雙線性內(nèi)聚力模型)中的界面剛度、最大名義應(yīng)力、界面失效位移的求解。當(dāng)加載速率為1.2 mm·min-1時(shí),模型參數(shù)分別為6.40 MPa·mm-1、0.72 MPa、1.80 mm。

2) 基于反演識(shí)別方法獲取的界面模型開展仿真計(jì)算。ε為8%、18%時(shí),仿真與實(shí)測(cè)ROI位移(Ux)的平均誤差分別為7.1%、4.9%、最大誤差分別為16.5%、12.8%、誤差小于10%的區(qū)域占比分別為78.6%、85.7%,均表明仿真結(jié)果與實(shí)驗(yàn)吻合較好,所建立的界面模型可以在一定程度上表征界面真實(shí)的本構(gòu)關(guān)系。

3) 獲取了界面單元的損傷因子變化曲線,與粘接試件拉伸過程的損傷情況有較好地吻合,驗(yàn)證了所建立的界面模型可正確地反映出粘接界面在拉伸過程中脫粘規(guī)律。

猜你喜歡
界面模型
一半模型
重要模型『一線三等角』
國企黨委前置研究的“四個(gè)界面”
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
基于FANUC PICTURE的虛擬軸坐標(biāo)顯示界面開發(fā)方法研究
空間界面
金秋(2017年4期)2017-06-07 08:22:16
電子顯微打開材料界面世界之門
人機(jī)交互界面發(fā)展趨勢(shì)研究
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 日韩欧美中文字幕在线韩免费| 国产凹凸视频在线观看| 亚洲精品大秀视频| 99久久99这里只有免费的精品| 久久久久国产精品免费免费不卡| 国产高清毛片| 国产尤物jk自慰制服喷水| 久久香蕉欧美精品| 日韩二区三区无| 日韩成人在线网站| 在线观看亚洲天堂| 亚洲国产天堂久久综合226114| 91福利片| 精品视频在线观看你懂的一区| 九色在线观看视频| 999福利激情视频| 午夜国产精品视频黄| 国产日韩AV高潮在线| 狠狠操夜夜爽| 在线播放精品一区二区啪视频| 亚洲人成电影在线播放| 国产一区二区在线视频观看| 国产精品香蕉在线观看不卡| 欧美翘臀一区二区三区| 婷婷亚洲最大| 亚洲天堂高清| 久久综合成人| 手机在线看片不卡中文字幕| 色呦呦手机在线精品| 一本久道热中字伊人| 99re视频在线| 欧美一区二区精品久久久| 久久性妇女精品免费| 欧美日韩亚洲国产主播第一区| 成人国产一区二区三区| 国产99精品视频| 在线一级毛片| 国产成人做受免费视频| 国产精品欧美在线观看| 亚洲天堂2014| 99精品国产电影| 最新国产网站| 成人蜜桃网| 亚洲成人精品| 久久久久免费精品国产| 蜜桃视频一区| 97影院午夜在线观看视频| 亚洲天堂区| 91久久精品国产| 亚洲色图欧美一区| 欧美一级在线| 丁香婷婷久久| 毛片免费网址| 综合天天色| 波多野结衣第一页| 欧美日韩精品在线播放| 国产一二视频| 国产精品久线在线观看| 午夜视频www| 久久久久青草线综合超碰| 久久国产V一级毛多内射| 国产亚洲视频播放9000| 青青草91视频| 亚洲天堂在线免费| 国产精品亚洲精品爽爽| 青青青视频免费一区二区| 91精品国产自产在线老师啪l| 免费99精品国产自在现线| 婷五月综合| 免费福利视频网站| 国产日韩欧美精品区性色| 婷婷六月综合网| 熟女成人国产精品视频| 国产激爽爽爽大片在线观看| 国产精品久久久久久久久| 国产精女同一区二区三区久| 一区二区三区四区精品视频| 久久久久无码国产精品不卡| 香蕉精品在线| 亚洲天堂自拍| 中文字幕免费在线视频| 国产免费久久精品99re丫丫一|