• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于深度學(xué)習(xí)和細(xì)觀力學(xué)的顆粒材料本構(gòu)關(guān)系研究1)

    2021-11-10 09:48:44瞿同明馮云田王孟琦趙婷婷狄少丞
    力學(xué)學(xué)報(bào) 2021年9期
    關(guān)鍵詞:深度模型

    瞿同明 馮云田 ,2) 王孟琦 趙婷婷 狄少丞

    * (斯旺西大學(xué)辛克維奇工程計(jì)算中心,英國(guó)斯旺西 SA1 8EP)

    ? (太原理工大學(xué)機(jī)械與運(yùn)載工程學(xué)院,太原 030024)

    ** (哈爾濱工程大學(xué)船舶工程學(xué)院,哈爾濱 150001)

    引言

    顆粒材料廣泛存在于各類(lèi)自然環(huán)境和工程活動(dòng)中,人類(lèi)幾乎所有的基礎(chǔ)設(shè)施都建造于巖土材料這類(lèi)典型的顆粒材料上,理解和預(yù)測(cè)此類(lèi)材料在外載作用下的力學(xué)響應(yīng)具有重要的意義.巖土材料通常由粒徑不一,形狀各異的礦物顆粒以及顆粒之間的空隙組成.顆粒間的互相滑動(dòng)會(huì)引起不可恢復(fù)的塑性變形,因而,顆粒材料在剪切作用下極易發(fā)生體積收縮或膨脹(剪脹性或者剪縮性).在外載荷作用下,顆粒材料主要通過(guò)顆粒之間的接觸力來(lái)傳遞和平衡外力,內(nèi)部常常形成高度復(fù)雜的力鏈網(wǎng)絡(luò)結(jié)構(gòu)(非均勻性)[1].這種力鏈網(wǎng)絡(luò)的穩(wěn)定性來(lái)源于顆粒間的摩擦強(qiáng)度和形狀引起的互鎖效應(yīng),宏觀上體現(xiàn)出顆粒材料的抗剪強(qiáng)度和剛度隨著壓應(yīng)力的增加而增大(壓硬性和摩擦性).除此以外,顆粒材料的力學(xué)響應(yīng)還和應(yīng)力路徑,材料初始各向異性等因素相關(guān)[2-3].這些復(fù)雜的特性使得預(yù)測(cè)顆粒材料的應(yīng)力?應(yīng)變關(guān)系成為一個(gè)極具挑戰(zhàn)的科學(xué)問(wèn)題[4-7].

    在過(guò)去幾十年間,基于唯象假設(shè)建立應(yīng)力?應(yīng)變?cè)隽恐g的解析關(guān)系一直是巖土等顆粒材料本構(gòu)問(wèn)題中最主流的研究方式,然而發(fā)展一個(gè)統(tǒng)一的應(yīng)力?應(yīng)變理論模型仍然存在許多困難.并且,許多現(xiàn)有的本構(gòu)模型數(shù)學(xué)形式復(fù)雜,需要借助精細(xì)的室內(nèi)實(shí)驗(yàn)來(lái)標(biāo)定數(shù)十個(gè)自由參數(shù),對(duì)一般工程師而言學(xué)習(xí)和使用成本過(guò)高.由于理論發(fā)展和實(shí)踐應(yīng)用的不相適應(yīng),盡管各類(lèi)高級(jí)模型層出不窮,大量的巖土從業(yè)人員仍然使用基于線(xiàn)彈性理想塑性的摩爾?庫(kù)倫等簡(jiǎn)單本構(gòu)模型開(kāi)展工程問(wèn)題分析.

    除唯象模型以外,一些學(xué)者也在尋找新的研究范式來(lái)研究本構(gòu)問(wèn)題.在層級(jí)多尺度模擬方法中,有限元等連續(xù)計(jì)算模型不再使用唯象的本構(gòu)模型,而是通過(guò)低尺度的離散元模擬來(lái)得到材料的應(yīng)力?應(yīng)變關(guān)系[8].在有限元的計(jì)算過(guò)程中,每個(gè)高斯點(diǎn)的位移梯度或者應(yīng)變信息實(shí)時(shí)傳遞給離散元作為邊界條件;離散元計(jì)算后,又將對(duì)應(yīng)的應(yīng)力響應(yīng)或者切向剛度矩陣實(shí)時(shí)傳遞給有限元.然而,此種多尺度計(jì)算方法中每個(gè)高斯點(diǎn)均需要有對(duì)應(yīng)的離散元模型,計(jì)算成本過(guò)高,很難真正應(yīng)用到工程尺度問(wèn)題的研究中.

    另外一種代表性的方案是采用神經(jīng)網(wǎng)絡(luò)從大量實(shí)驗(yàn)數(shù)據(jù)中訓(xùn)練代理本構(gòu)模型.這種模型無(wú)需標(biāo)定自由參數(shù),在有限元等宏觀計(jì)算方法中使用方便,并且隨著新數(shù)據(jù)的不斷增加,其預(yù)測(cè)精度可以不斷提高,相比傳統(tǒng)本構(gòu)模型有獨(dú)特的優(yōu)勢(shì).使用神經(jīng)網(wǎng)絡(luò)來(lái)訓(xùn)練巖土材料本構(gòu)模型的思想可以追溯到20 世紀(jì)90 年代[9-12].近年來(lái),隨著深度學(xué)習(xí)技術(shù)的突破性進(jìn)展,基于神經(jīng)網(wǎng)絡(luò)的本構(gòu)模型再次得到廣泛的重視.Wang 等[13-15]結(jié)合博弈論和強(qiáng)化學(xué)習(xí)模型研究了顆粒材料的界面張力和位移的本構(gòu)關(guān)系以及應(yīng)力?應(yīng)變關(guān)系;Zhang 等[16]使用長(zhǎng)短期記憶模型(LSTM)預(yù)測(cè)了顆粒材料的循環(huán)剪切行為;Karapiperis 等[17]在有限元模型中實(shí)施了基于離散元數(shù)據(jù)訓(xùn)練得到的神經(jīng)網(wǎng)絡(luò)本構(gòu)模型.Qu 等[18]基于離散元數(shù)據(jù)探索了3 類(lèi)不同的顆粒材料本構(gòu)訓(xùn)練策略.盡管深度學(xué)習(xí)模型已經(jīng)被證明能有效捕捉特定工況下的顆粒材料力學(xué)響應(yīng),但是要將此類(lèi)純粹的數(shù)據(jù)驅(qū)動(dòng)模型發(fā)展成為可靠可用的,適用于多種加載工況的本構(gòu)模型仍有較遠(yuǎn)的距離.

    解析模型和數(shù)據(jù)模型作為兩種不同的科學(xué)研究范式具有不同的優(yōu)缺點(diǎn).解析模型基于一些假設(shè)條件,利用理論方法構(gòu)造整個(gè)預(yù)測(cè)模型.此種范式的優(yōu)點(diǎn)是模型在自由參數(shù)得到標(biāo)定后具有一定的預(yù)測(cè)能力,但是受限于假設(shè)條件的約束,面對(duì)復(fù)雜荷載工況時(shí)往往預(yù)測(cè)與實(shí)際存在不少出入.數(shù)據(jù)模型不需要任何假設(shè),完全利用經(jīng)驗(yàn)數(shù)據(jù)來(lái)構(gòu)造預(yù)測(cè)模型.然而一個(gè)純粹的數(shù)據(jù)驅(qū)動(dòng)模型需要大量且覆蓋全面的訓(xùn)練樣本,對(duì)于本構(gòu)問(wèn)題來(lái)說(shuō),獲取大量的材料應(yīng)力?應(yīng)變數(shù)據(jù)并不容易.另外,當(dāng)數(shù)據(jù)模型遇到超出經(jīng)驗(yàn)樣本以外的情況時(shí),往往預(yù)測(cè)能力有限.本文擬嘗試結(jié)合理論模型和數(shù)據(jù)驅(qū)動(dòng)模型兩種不同的研究范式,通過(guò)將客觀的力學(xué)規(guī)律與應(yīng)力?應(yīng)變數(shù)據(jù)經(jīng)驗(yàn)結(jié)合起來(lái),發(fā)展顆粒材料的數(shù)據(jù)驅(qū)動(dòng)本構(gòu)模型訓(xùn)練方法.

    本文首先基于Vogit 假設(shè)推導(dǎo)小應(yīng)變條件下顆粒材料的應(yīng)力?應(yīng)變解析關(guān)系,利用此關(guān)系識(shí)別出一組與顆粒材料本構(gòu)響應(yīng)相關(guān)的重要變量.通過(guò)有向圖的方式,將所有與本構(gòu)響應(yīng)相關(guān)的重要變量包含到從應(yīng)變到應(yīng)力的預(yù)測(cè)中,變量間的映射關(guān)系通過(guò)時(shí)序的深度學(xué)習(xí)模型來(lái)描述.采用離散元常規(guī)三軸以及真三軸數(shù)值試驗(yàn)生成顆粒材料應(yīng)力?應(yīng)變數(shù)據(jù)用于訓(xùn)練,驗(yàn)證和測(cè)試本文提出的深度學(xué)習(xí)模型.

    1 基于深度學(xué)習(xí)的數(shù)據(jù)驅(qū)動(dòng)本構(gòu)模型

    1.1 基本原理

    深度學(xué)習(xí)是機(jī)器學(xué)習(xí)算法中的一類(lèi),主要指基于多層人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的代理模型方法.數(shù)學(xué)形式上,神經(jīng)網(wǎng)絡(luò)是由線(xiàn)性的矩陣運(yùn)算和非線(xiàn)性的激活函數(shù)構(gòu)成的復(fù)合函數(shù),這類(lèi)代理函數(shù)具有強(qiáng)大的表示能力,已經(jīng)被證明能夠以任意的精度來(lái)估計(jì)任何一個(gè)定義在實(shí)數(shù)空間中的有界閉集函數(shù)[19-20].

    一個(gè)問(wèn)題是否適合深度學(xué)習(xí)方法,主要有3 個(gè)基本前提:(1)問(wèn)題本身包含某種確定的模式或者規(guī)律;(2)無(wú)法在數(shù)學(xué)上找到精確的解析方法來(lái)描述這種模式;(3)有包含這個(gè)問(wèn)題所有自變量和因變量的完備數(shù)據(jù).對(duì)于顆粒材料而言,本構(gòu)關(guān)系由材料自身物理性質(zhì)所決定,這種關(guān)系必然是受物理材料規(guī)律支配的;并且,已有的大量研究表明:建立一個(gè)精確、統(tǒng)一的解析方程來(lái)描述這種關(guān)系極為困難;再者,通過(guò)各類(lèi)顆粒材料加載試驗(yàn),我們可以得到反映本構(gòu)關(guān)系的大量應(yīng)力?應(yīng)變數(shù)據(jù).這些因素構(gòu)成了發(fā)展基于深度學(xué)習(xí)框架的數(shù)據(jù)驅(qū)動(dòng)材料本構(gòu)關(guān)系的前提.

    一般而言,在準(zhǔn)靜態(tài)變形情況下,材料內(nèi)部一點(diǎn)的應(yīng)力張量σij和應(yīng)變張量εmn可以用6 個(gè)獨(dú)立的分量來(lái)表示,即

    如圖1 所示,基于深度學(xué)習(xí)的本構(gòu)模型即是利用神經(jīng)網(wǎng)絡(luò)來(lái)建立材料應(yīng)力和應(yīng)變之間的映射關(guān)系.通常輸入?yún)?shù)為應(yīng)變張量,輸出參數(shù)為應(yīng)力張量,中間隱藏層的激活函數(shù)可按如下方式計(jì)算

    圖1 基于深度學(xué)習(xí)的本構(gòu)模型示意圖Fig.1 Diagram of deep learning-based constitutive models

    假設(shè)第j層為最后一個(gè)隱藏層,第(j+1)層為輸出層,那么應(yīng)力的預(yù)測(cè)值可按照如下公式計(jì)算

    其中,e為最后一個(gè)隱藏層的單元數(shù)量,W和b分別為神經(jīng)網(wǎng)絡(luò)的權(quán)重和偏置系數(shù).

    在一個(gè)確定的神經(jīng)網(wǎng)絡(luò)架構(gòu)下尋找描述對(duì)應(yīng)映射關(guān)系的網(wǎng)絡(luò)權(quán)重和偏置的過(guò)程稱(chēng)為訓(xùn)練或者學(xué)習(xí).如圖2 所示,一個(gè)神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程主要包括前向傳播、反向傳播和迭代更新網(wǎng)絡(luò)參數(shù)(權(quán)重和偏置).在訓(xùn)練前,網(wǎng)絡(luò)參數(shù)通常是隨機(jī)初始化的,此時(shí)通過(guò)前向傳播計(jì)算得到的預(yù)測(cè)值必然是偏離實(shí)際值的.這個(gè)預(yù)測(cè)值和實(shí)際值的差別通常用損失函數(shù)L來(lái)衡量.訓(xùn)練一個(gè)可靠的預(yù)測(cè)模型的過(guò)程可以轉(zhuǎn)化為一個(gè)以最小化損失函數(shù)L為目標(biāo)的優(yōu)化問(wèn)題,即

    圖2 人工神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程Fig.2 Basic procedures of training artificial neural networks

    其中,argmin 是指使得目標(biāo)函數(shù)取最小值時(shí)的自變量值.

    常規(guī)神經(jīng)網(wǎng)絡(luò)一般通過(guò)基于一階梯度的優(yōu)化算法,如自適應(yīng)矩估計(jì)算法(Adam)等來(lái)不斷迭代,從而找到使得L達(dá)到極小值時(shí)的網(wǎng)絡(luò)參數(shù).現(xiàn)代神經(jīng)網(wǎng)絡(luò)中,這些梯度信息都使用反向傳播(BP)算法通過(guò)自動(dòng)微分和鏈?zhǔn)椒▌t計(jì)算得到.對(duì)于一個(gè)已經(jīng)訓(xùn)練好的模型,也即網(wǎng)絡(luò)中的權(quán)重和偏置參數(shù)能很好地描述映射關(guān)系時(shí),我們可以通過(guò)向前傳播,計(jì)算得到較為可靠的應(yīng)力預(yù)測(cè)值.

    1.2 材料本構(gòu)模型的時(shí)序本質(zhì)與深度學(xué)習(xí)模型的選用

    顆粒材料應(yīng)力響應(yīng)具有明顯的應(yīng)變歷史依賴(lài)和應(yīng)變路徑依賴(lài)特征,一個(gè)能反映復(fù)雜應(yīng)力?應(yīng)變路徑條件下的材料本構(gòu)響應(yīng)模型,必須具備“記憶”應(yīng)力或者應(yīng)變歷史的能力[21],因此本構(gòu)模型本質(zhì)上是一種時(shí)序模型或者過(guò)程模型.在當(dāng)前深度學(xué)習(xí)模型中,LSTM 和GRU 等循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)模型主要適用于包含時(shí)間序列的預(yù)測(cè)問(wèn)題.由于幾種RNN 模型預(yù)測(cè)能力近似,同等網(wǎng)絡(luò)規(guī)模條件下,GRU 相比LSTM 需要更少的訓(xùn)練參數(shù),因此本文主要采用GRU 作為基本的深度學(xué)習(xí)模型來(lái)訓(xùn)練材料的本構(gòu)行為.

    2 基于Vogit 均質(zhì)化假設(shè)的顆粒材料小應(yīng)變應(yīng)力?應(yīng)變關(guān)系

    本文嘗試引入顆粒物質(zhì)力學(xué)理論來(lái)指導(dǎo)深度學(xué)習(xí),以期降低深度學(xué)習(xí)在訓(xùn)練階段對(duì)數(shù)據(jù)量的需求,提高模型在預(yù)測(cè)階段的泛化能力.對(duì)于顆粒材料而言,目前主要有兩種唯象本構(gòu)模型:連續(xù)模型和細(xì)觀力學(xué)模型.基于連續(xù)假設(shè)的彈塑性模型采用屈服面來(lái)描述材料的塑性,借助流動(dòng)法則來(lái)描述塑性的發(fā)展.細(xì)觀力學(xué)模型雖然不再需要假想的屈服面,但仍需要通過(guò)預(yù)測(cè)細(xì)觀組構(gòu)的演化來(lái)描述應(yīng)力?應(yīng)變關(guān)系.不管是基于連續(xù)假設(shè)還是細(xì)觀力學(xué)的唯象模型,在預(yù)測(cè)材料隨外部應(yīng)變引起的應(yīng)力響應(yīng)時(shí),常常需要使用一些從試驗(yàn)數(shù)據(jù)中擬合得來(lái)的參數(shù)(如硬化參數(shù)等)來(lái)描述主應(yīng)力?主應(yīng)變關(guān)系.為避免使用此類(lèi)擬合參數(shù),在準(zhǔn)靜態(tài)加載條件下,本研究將時(shí)序的應(yīng)力?應(yīng)變過(guò)程簡(jiǎn)化成單個(gè)時(shí)刻的應(yīng)力?應(yīng)變響應(yīng)來(lái)考慮,即小應(yīng)變條件下的顆粒材料應(yīng)力?應(yīng)變關(guān)系.

    散粒體在小變形情況下,顆粒體系內(nèi)的彈性勢(shì)能U可按照如下方式計(jì)算

    當(dāng)顆粒接觸為線(xiàn)性彈性接觸時(shí),儲(chǔ)存在接觸k內(nèi)的彈性勢(shì)能Uk為

    式中,kn和ks為顆粒之間的法向和切向剛度.

    如圖3 所示,通過(guò)將顆粒位移等效為相應(yīng)連續(xù)體內(nèi)對(duì)應(yīng)點(diǎn)處的位移,則顆粒相對(duì)位移與對(duì)應(yīng)連續(xù)體內(nèi)的局部應(yīng)變(顆粒A 與B 之間的等效應(yīng)變)可表示為如下關(guān)系

    圖3 顆粒接觸與接觸位移Fig.3 A contact between particles and contact displacements

    根據(jù)彈性力學(xué)理論,連續(xù)體的應(yīng)力張量可通過(guò)應(yīng)變能密度u對(duì)應(yīng)變張量求取偏導(dǎo)得到,即

    假設(shè)顆粒試樣內(nèi)的變形是統(tǒng)計(jì)均勻的(Vogit’s 假設(shè)),即在對(duì)應(yīng)連續(xù)體內(nèi)任意位置處的整體應(yīng)變 εij等于材料中沿著顆粒接觸方向的局部應(yīng)變

    基于式(13),假設(shè)散粒體和對(duì)應(yīng)連續(xù)體在小應(yīng)變情況下的彈性勢(shì)能相等,則式(12)可表示為

    其中,V為顆粒試樣的體積.

    將應(yīng)力張量對(duì)應(yīng)變張量求偏導(dǎo),可得到對(duì)應(yīng)的彈性剛度矩陣為

    式中,δin為克羅內(nèi)克函數(shù),當(dāng)i=n時(shí),δin=1;當(dāng)i≠n時(shí),δin=0.若令

    式(17) 表明Cijmn由顆粒間的法向接觸剛度kn、切向剛度ks以及兩個(gè)單位體積內(nèi)的組構(gòu)張量aijmn和bijmn構(gòu)成.在顆粒試樣承受外界荷載作用時(shí),假設(shè)顆粒未發(fā)生破碎或其他材料劣化,則kn和ks保持不變;隨著試樣在宏觀上發(fā)生不可恢復(fù)的塑性變形(微觀上表現(xiàn)為顆粒間發(fā)生相互滑移)時(shí),組構(gòu)張量aijmn和bijmn隨之不斷演化.

    小應(yīng)變情況下的顆粒材料的宏觀應(yīng)力?應(yīng)變關(guān)系可表示為當(dāng)切應(yīng)力和切應(yīng)變?yōu)? 時(shí),3 個(gè)主應(yīng)力主應(yīng)變的增量關(guān)系展開(kāi)后為

    其中,C1122=C2211,C1133=C3311,C2233=C3322.因而彈性剛度矩陣中共有6 個(gè)獨(dú)立分量描述主應(yīng)力與主應(yīng)變的關(guān)系.對(duì)于準(zhǔn)靜態(tài)加載而言,式(19)的關(guān)系在試樣加載過(guò)程中總是成立,但隨著顆粒材料組構(gòu)的不斷演化,彈性剛度矩陣Cijmn中各個(gè)分量的值也會(huì)隨之變化.基于細(xì)觀力學(xué)的唯象模型通常需要借助流動(dòng)法則和硬化函數(shù)等來(lái)擬合得到試樣內(nèi)部組構(gòu)的演化過(guò)程,本文中,這些微觀結(jié)構(gòu)的演化過(guò)程則直接由深度學(xué)習(xí)來(lái)預(yù)測(cè)得到.

    3 力學(xué)知識(shí)指導(dǎo)下的深度學(xué)習(xí)策略

    由于顆粒材料內(nèi)部描述細(xì)觀結(jié)構(gòu)的組構(gòu)張量屬于內(nèi)變量,在有限元等連續(xù)介質(zhì)方法中計(jì)算應(yīng)力?應(yīng)變關(guān)系時(shí),此類(lèi)內(nèi)變量不能直接作為輸入.本文嘗試引入有向圖理論,將微觀結(jié)構(gòu)的演化過(guò)程包含在應(yīng)力?應(yīng)變關(guān)系的描述中.

    有向圖是圖論中用于描述物與物之間聯(lián)系的一種方法,主要由結(jié)點(diǎn)和含有方向的邊(連接兩個(gè)結(jié)點(diǎn))構(gòu)成.結(jié)點(diǎn)表示某一個(gè)物理量或者對(duì)象,邊線(xiàn)表示兩個(gè)物理量之間存在某種聯(lián)系并且這種聯(lián)系是從源結(jié)點(diǎn)(尾部)指向目標(biāo)結(jié)點(diǎn)(頭部)的.根據(jù)有向圖理論和小應(yīng)變條件下得到的顆粒材料應(yīng)力?應(yīng)變關(guān)系,本文設(shè)計(jì)的本構(gòu)關(guān)系訓(xùn)練策略如圖4 所示.

    圖4 基于有向圖包含組構(gòu)演化的本構(gòu)訓(xùn)練方式Fig.4 A directed graph-based constitutive training approach incorporating fabric evolution

    圖4 中,初始結(jié)點(diǎn)為主應(yīng)變,作為輸入?yún)?shù);目標(biāo)結(jié)點(diǎn)為主應(yīng)力,作為輸出參數(shù);中間節(jié)點(diǎn)為6 個(gè)相互獨(dú)立的彈性剛度張量分量.每條邊之間的映射關(guān)系通過(guò)深度神經(jīng)網(wǎng)絡(luò)來(lái)描述.將該有向圖展開(kāi)后,可用兩條子神經(jīng)網(wǎng)絡(luò)來(lái)描述整個(gè)有向圖結(jié)構(gòu):

    在網(wǎng)絡(luò)學(xué)習(xí)階段,兩條子網(wǎng)絡(luò)單獨(dú)訓(xùn)練;在使用階段,兩條子網(wǎng)絡(luò)合并成一條完整的信息流:子網(wǎng)絡(luò)1 的輸出連同3 個(gè)主應(yīng)變序列一起作為子網(wǎng)絡(luò)2 的輸入,共同預(yù)測(cè)材料的應(yīng)力響應(yīng).因此,基于有向圖的神經(jīng)網(wǎng)絡(luò)模型雖然充分利用了顆粒材料內(nèi)部的微觀結(jié)構(gòu)演化信息,但在使用過(guò)程中并不違背本構(gòu)模型僅僅將主應(yīng)變和主應(yīng)力作為最終輸入和輸出信息的基本原則.

    4 數(shù)據(jù)準(zhǔn)備與深度學(xué)習(xí)模型的實(shí)施

    4.1 基于離散元的數(shù)值三軸試驗(yàn)與模擬工況

    細(xì)觀尺度的離散元仿真已經(jīng)被證明能夠有效模擬顆粒材料的本構(gòu)行為[22],并且在層級(jí)多尺度計(jì)算模型中,離散元模擬直接代替了傳統(tǒng)的唯象本構(gòu)關(guān)系進(jìn)行有限元計(jì)算[23].本文采用離散元模擬三軸實(shí)驗(yàn)來(lái)產(chǎn)生顆粒材料的應(yīng)力?應(yīng)變數(shù)據(jù),用于深度學(xué)習(xí)訓(xùn)練.如圖5 所示是本文采用的數(shù)值三軸試樣,每個(gè)試樣共包含4037 個(gè)顆粒,半徑為2~ 4 mm 均勻分布;顆粒之間采用線(xiàn)性接觸模型,法向和切向接觸剛度分別為105N/m 和5 × 104N/m;顆粒密度為2600 kg/m3,局部阻尼系數(shù)為0.5.為生成一個(gè)較密實(shí)的均勻試樣,在指定空間內(nèi)生成顆粒后,采用一個(gè)較小的制樣顆粒摩擦系數(shù)(如 0.05),使用伺服方法通過(guò)移動(dòng)試樣的6 個(gè)邊界來(lái)施加三向等壓壓力,直到試樣穩(wěn)定地達(dá)到200 kPa 的目標(biāo)圍壓值,再將制樣階段的接觸摩擦系數(shù)改成實(shí)際樣本摩擦系數(shù)0.5.

    圖5 離散元三軸試驗(yàn)?zāi)P虵ig.5 Triaxial compression models via discrete element modelling

    本文主要考慮了3 類(lèi)典型的三軸試驗(yàn)工況,具體包括:

    (1)包含應(yīng)力?應(yīng)變加卸載循環(huán)的常規(guī)三軸壓縮試驗(yàn)291 組;

    (2)平均有效應(yīng)力p恒定的真三軸加卸載試驗(yàn)122組(p=σ11+σ22+σ33,σ11=σ22);

    (3)主應(yīng)力系數(shù)b恒定的真三軸壓縮試驗(yàn)20 組

    所有工況的最大軸向應(yīng)變?yōu)?2%,在每組加卸載工況中,卸載應(yīng)變和再加載應(yīng)變的取值各不相同,在等b真三軸加載中,b以0.05 的間隔在(0,1]之間取值.在440 組所有三軸試驗(yàn)中,288 組用于訓(xùn)練,72 組用于驗(yàn)證(validation),剩下的80 組用于測(cè)試(test).注意驗(yàn)證不同于測(cè)試,驗(yàn)證組用于神經(jīng)網(wǎng)絡(luò)模型在訓(xùn)練過(guò)程中調(diào)整超參數(shù)組合,而測(cè)試集則僅用于測(cè)試目的,測(cè)試數(shù)據(jù)在模型訓(xùn)練過(guò)程中從未接觸到神經(jīng)網(wǎng)絡(luò),這樣的設(shè)計(jì)可避免在測(cè)試神經(jīng)網(wǎng)絡(luò)模型時(shí)發(fā)生數(shù)據(jù)泄露的可能性,從而保證神經(jīng)網(wǎng)絡(luò)測(cè)試結(jié)果的可靠性.

    4.2 深度學(xué)習(xí)模型實(shí)施過(guò)程

    在構(gòu)建深度學(xué)習(xí)模型之前,離散元模型中計(jì)算得出的數(shù)據(jù)需要預(yù)處理.一方面,現(xiàn)實(shí)問(wèn)題中的輸入變量往往具有不同的數(shù)據(jù)大小、單位以及分布特征,而差異較大的原始輸入數(shù)據(jù)會(huì)降低訓(xùn)練效率,阻礙神經(jīng)網(wǎng)絡(luò)在訓(xùn)練過(guò)程中的收斂.因此,原始數(shù)據(jù)需要通過(guò)標(biāo)準(zhǔn)化將其縮放到平均值為0、標(biāo)準(zhǔn)差為1 的分布空間,這將有利于降低訓(xùn)練過(guò)程中網(wǎng)絡(luò)陷入局部最優(yōu)的可能性.另一方面,循環(huán)神經(jīng)網(wǎng)絡(luò)模型的輸入數(shù)據(jù)需要滿(mǎn)足某些特定特征:所有的輸入數(shù)據(jù)必須是三維矩陣的形式,其中第一維表示樣本,一段應(yīng)力?應(yīng)變響應(yīng)的序列數(shù)據(jù)即是一個(gè)樣本;第二維是時(shí)間步(time steps),也稱(chēng)移動(dòng)窗口的長(zhǎng)度,這是指循環(huán)神經(jīng)網(wǎng)絡(luò)在給出下一個(gè)預(yù)測(cè)輸出時(shí)需要依據(jù)在此之前的輸入變量的數(shù)量,若時(shí)間步取20,則在預(yù)測(cè)當(dāng)前應(yīng)力時(shí),需要輸入當(dāng)前應(yīng)變值以及在此之前的19 個(gè)應(yīng)變值;第三維則是數(shù)據(jù)特征,比如3 個(gè)主應(yīng)變分量可當(dāng)作3 個(gè)數(shù)據(jù)特征.

    本文的深度學(xué)習(xí)模型通過(guò)Tensorflow 和keras來(lái)構(gòu)建和訓(xùn)練.神經(jīng)網(wǎng)絡(luò)的大小(層數(shù)和神經(jīng)元數(shù)量)決定了神經(jīng)網(wǎng)絡(luò)所能達(dá)到的上限表示能力、網(wǎng)絡(luò)訓(xùn)練過(guò)程中學(xué)習(xí)速率等超參數(shù)則決定了模型能在當(dāng)前數(shù)據(jù)集下將其預(yù)測(cè)能力發(fā)揮到何種水平.本文采用的網(wǎng)絡(luò)結(jié)構(gòu)和超參數(shù)組合通過(guò)試錯(cuò)法得到:兩個(gè)子網(wǎng)絡(luò)除了輸入層和輸出層以外,還包括兩個(gè)含有120 個(gè)神經(jīng)元的GRU 隱藏層和20 個(gè)神經(jīng)元的全連接隱藏層.全連接層采用Sigmoid 激活函數(shù),輸出層采用linear 激活函數(shù).移動(dòng)窗口的長(zhǎng)度為40,采用Adam 優(yōu)化器訓(xùn)練神經(jīng)網(wǎng)絡(luò);學(xué)習(xí)速率為0.01,批量數(shù)據(jù)的尺寸(batch size)為128.損失函數(shù)為預(yù)測(cè)值與實(shí)際值的平均絕對(duì)誤差函數(shù)(MAE).模型在訓(xùn)練了500 個(gè)epochs 后基本達(dá)到收斂狀態(tài),兩個(gè)神經(jīng)網(wǎng)絡(luò)在訓(xùn)練集和驗(yàn)證集上的預(yù)測(cè)表現(xiàn)分別如圖6 中的學(xué)習(xí)曲線(xiàn)所示.

    圖6 學(xué)習(xí)曲線(xiàn)Fig.6 Learning curves

    需要說(shuō)明的是,選擇40 的移動(dòng)窗口長(zhǎng)度并不意味著不關(guān)注每條應(yīng)力-應(yīng)變曲線(xiàn)上前39 步應(yīng)變引起的應(yīng)力值,整個(gè)序列中的應(yīng)力響應(yīng)都應(yīng)考慮進(jìn)去.實(shí)踐中為了確保總能得到長(zhǎng)度為40 的應(yīng)變序列,對(duì)于前40 組應(yīng)力值,我們通過(guò)在第一個(gè)應(yīng)變值前面不斷充填零值,使得應(yīng)變序列的長(zhǎng)度總為選定的移動(dòng)窗口長(zhǎng)度.這樣的措施確保了整個(gè)應(yīng)力序列都能夠得到預(yù)測(cè),實(shí)踐表明具有較好的預(yù)測(cè)效果.此外,充零的方式還具有較合理的物理意義,即在荷載作用之前,試樣并不產(chǎn)生應(yīng)變響應(yīng).

    4.3 深度學(xué)習(xí)模型的評(píng)價(jià)指標(biāo)

    預(yù)測(cè)模型的好壞通過(guò)量化預(yù)測(cè)值與實(shí)際值的平均差異得到.將平均絕對(duì)誤差函數(shù)作為損失函數(shù)是評(píng)價(jià)模型預(yù)測(cè)能力的一種有效方式,但此種誤差評(píng)價(jià)指標(biāo)往往與訓(xùn)練對(duì)象的值有關(guān),且缺乏對(duì)預(yù)測(cè)能力的直觀體現(xiàn).本文除了采用平均絕對(duì)誤差函數(shù)外,還采用另外一種評(píng)分的方式來(lái)對(duì)模型進(jìn)行評(píng)價(jià).該方法首先需要計(jì)算得到第j條應(yīng)力?應(yīng)變曲線(xiàn)上第i個(gè)點(diǎn)的平方誤差損失SSEij

    式中,Nj為第j條應(yīng)力?應(yīng)變曲線(xiàn)上的樣本點(diǎn)數(shù)目,將曲線(xiàn)上所有點(diǎn)的SSEij按照升序排列,可以為每條曲線(xiàn)計(jì)算一個(gè)對(duì)應(yīng)的得分

    式中,εP%是累計(jì)分布函數(shù)上對(duì)應(yīng)P%時(shí)的SSE值,這個(gè)值被用來(lái)作為一個(gè)代表值來(lái)評(píng)價(jià)整條應(yīng)力?應(yīng)變曲線(xiàn)上的預(yù)測(cè)得分;當(dāng)預(yù)測(cè)應(yīng)力?應(yīng)變序列中εP%對(duì)應(yīng)的SSE小于εcrit時(shí),可認(rèn)為預(yù)測(cè)足夠準(zhǔn)確,得分為1.在本文中,P%取90%,εcrit為0.001.

    以上介紹了單條應(yīng)力?應(yīng)變預(yù)測(cè)曲線(xiàn)的打分標(biāo)準(zhǔn),在評(píng)價(jià)每個(gè)神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)能力時(shí),一個(gè)更合理的指標(biāo)是判斷該模型在測(cè)試集上所有樣本的平均預(yù)測(cè)表現(xiàn),即使用已經(jīng)訓(xùn)練好的模型,得到對(duì)所有測(cè)試樣本上的平均MAE 值或者平均得分.

    5 預(yù)測(cè)結(jié)果

    5.1 總體預(yù)測(cè)表現(xiàn)

    在80 組測(cè)試樣本中,該深度學(xué)習(xí)模型在測(cè)試集上的平均MAE 值為0.013;平均得分為0.986;87.5%的測(cè)試樣本(70 組)獲得了滿(mǎn)分,其中最差的預(yù)測(cè)得分為0.639,對(duì)應(yīng)MAE 值為0.051;次差的一組預(yù)測(cè)得分為0.823,對(duì)應(yīng)MAE 值為0.031;在70 組滿(mǎn)分預(yù)測(cè)中,最低的MAE 值為0.004 2 和0.004 5.為直觀體現(xiàn)模型的預(yù)測(cè)能力,兩組最好和最差的預(yù)測(cè)情況如圖7 所示.

    圖7 兩組最佳與最差預(yù)測(cè)Fig.7 Examples of the two best and worst predictions

    在得分最低的兩組預(yù)測(cè)中,盡管在加卸載處的應(yīng)力?應(yīng)變預(yù)測(cè)值與實(shí)際值存在差別,但總體的應(yīng)力響應(yīng)趨勢(shì)已經(jīng)被神經(jīng)網(wǎng)絡(luò)模型捕捉到.在損失函數(shù)MAE 最低的兩組預(yù)測(cè)中,深度學(xué)習(xí)模型則幾乎重現(xiàn)了實(shí)際的應(yīng)力?應(yīng)變響應(yīng).

    5.2 內(nèi)插與外推預(yù)測(cè)表現(xiàn)

    根據(jù)已有數(shù)據(jù)以及模型預(yù)測(cè)未知區(qū)域的目標(biāo)對(duì)象時(shí),若預(yù)測(cè)點(diǎn)在已有數(shù)據(jù)范圍內(nèi)通常稱(chēng)為內(nèi)插預(yù)測(cè);若數(shù)據(jù)點(diǎn)完全在已有數(shù)據(jù)范圍之外,則為外插預(yù)測(cè).本文所用的訓(xùn)練數(shù)據(jù)集僅在常規(guī)三軸和等p真三軸中包含了一次應(yīng)力?應(yīng)變加卸載循環(huán),在等b真三軸壓縮工況中僅包含單調(diào)加載工況.圖8 分別給出了在測(cè)試集內(nèi)的常規(guī)三軸和等p真三軸一次加卸載,以及等b真三軸加載工況中的幾個(gè)代表性?xún)?nèi)插預(yù)測(cè).

    圖8 幾組代表性?xún)?nèi)插預(yù)測(cè)Fig.8 Some representative interpolation predictions

    注意到在常規(guī)三軸加載工況下,由于水平圍壓在加載過(guò)程中保持穩(wěn)定,σ22和σ33均保持不變;在平均有效應(yīng)力p恒定的真三軸加載工況中,σ22和σ33保持了相同的變化趨勢(shì);在中主應(yīng)力系數(shù)b恒定的真三軸加載工況中,σ33維持大小不變.經(jīng)過(guò)訓(xùn)練的神經(jīng)網(wǎng)絡(luò)模型均準(zhǔn)確預(yù)測(cè)到了這些趨勢(shì).

    本研究在測(cè)試數(shù)據(jù)集中準(zhǔn)備了常規(guī)三軸和等p真三軸單調(diào)加載工況,也準(zhǔn)備了6 組包含2 次及以上加卸載循環(huán)的常規(guī)三軸加載工況.這些工況完全在模型訓(xùn)練過(guò)程中所用的訓(xùn)練集和測(cè)試集范圍之外,因此被用來(lái)測(cè)試神經(jīng)網(wǎng)絡(luò)的外推預(yù)測(cè)能力.8 組測(cè)試樣本中除了一組預(yù)測(cè)得分為0.98 以外,其余預(yù)測(cè)都得到了滿(mǎn)分.幾組代表性預(yù)測(cè)工況如圖9 所示,結(jié)果表明深度學(xué)習(xí)模型能準(zhǔn)確地適應(yīng)單調(diào)加載和多次加卸載的外推工況.

    圖9 幾種代表性外推預(yù)測(cè)Fig.9 Some representative extrapolation predictions

    6 討論

    顆粒材料的多次加卸載等復(fù)雜工況本質(zhì)上都體現(xiàn)為細(xì)觀結(jié)構(gòu)的往復(fù)演化.在材料參數(shù)不發(fā)生較大變化的情況下,預(yù)測(cè)試樣應(yīng)力?應(yīng)變響應(yīng)的關(guān)鍵在于預(yù)測(cè)試樣組構(gòu)的演化.由于采用了基于時(shí)序的循環(huán)神經(jīng)網(wǎng)絡(luò)模型,在預(yù)測(cè)當(dāng)前應(yīng)力值時(shí),包含了足夠長(zhǎng)的應(yīng)變歷史信息來(lái)推斷當(dāng)前的組構(gòu)演化和應(yīng)力響應(yīng),因此深度學(xué)習(xí)模型能有效區(qū)分加載和卸載,做出準(zhǔn)確的應(yīng)力預(yù)測(cè).另外,盡管包含多次加卸載的工況未曾直接出現(xiàn)在模型訓(xùn)練過(guò)程中,但是類(lèi)似的細(xì)觀組構(gòu)演化過(guò)程可能已經(jīng)被訓(xùn)練數(shù)據(jù)所覆蓋,這也是深度學(xué)習(xí)模型對(duì)多次加卸載工況具備較好外推能力的原因.外推預(yù)測(cè)結(jié)果表明當(dāng)訓(xùn)練集中包含豐富的一次加卸載樣本時(shí),模型將具有一定的兩次或三次等加卸載工況的外推預(yù)測(cè)能力,這在實(shí)踐中能有效降低準(zhǔn)備訓(xùn)練樣本的數(shù)據(jù)規(guī)模.

    相比于顆粒材料室內(nèi)三軸試驗(yàn),基于離散元的虛擬三軸試驗(yàn)所需成本更低,效率更高,且方便模擬各類(lèi)加載工況.采用基于離散元的虛擬三軸試驗(yàn)產(chǎn)生顆粒材料應(yīng)力?應(yīng)變數(shù)據(jù)主要包含兩個(gè)方面的意義:(1) 離散元方法已經(jīng)表明具備定性再現(xiàn)顆粒材料基本響應(yīng)特征的能力[24],通過(guò)低成本的離散元模擬,可為探索數(shù)據(jù)驅(qū)動(dòng)本構(gòu)模型提供經(jīng)驗(yàn).(2) 高質(zhì)量的離散元模擬本身有能力對(duì)實(shí)際顆粒材料的宏觀力學(xué)響應(yīng)做出定量預(yù)測(cè).在離散元的參數(shù)得到合理標(biāo)定[25-26],使用更加真實(shí)合理的接觸模型,并考慮復(fù)雜形狀因素的離散元模擬中[27-28],實(shí)際顆粒材料的本構(gòu)行為能很好的被離散元所捕捉到,此時(shí)基于離散元數(shù)據(jù)基本能代表真實(shí)材料數(shù)據(jù)[22,29-30].

    本文的理論推導(dǎo)主要基于線(xiàn)彈性接觸模型,由于從線(xiàn)彈性接觸到非線(xiàn)性彈性接觸的映射關(guān)系并不復(fù)雜,深度神經(jīng)網(wǎng)絡(luò)完全有能力捕捉到這個(gè)映射關(guān)系的轉(zhuǎn)變.因此本文提出的深度學(xué)習(xí)框架也適用于預(yù)測(cè)非線(xiàn)性彈性接觸模型的應(yīng)力?應(yīng)變響應(yīng).此外,若需嚴(yán)格地發(fā)展基于非線(xiàn)性彈性接觸模型的深度學(xué)習(xí)策略,可在本文提出的訓(xùn)練框架下,參考文獻(xiàn)[25]中推導(dǎo)的基于非線(xiàn)性彈性接觸模型的顆粒材料應(yīng)力?應(yīng)變關(guān)系進(jìn)行深度學(xué)習(xí)訓(xùn)練.

    數(shù)據(jù)驅(qū)動(dòng)本構(gòu)模型相比于傳統(tǒng)基于唯象假設(shè)的解析模型有幾個(gè)突出優(yōu)點(diǎn):(1)可以根據(jù)新數(shù)據(jù)的補(bǔ)充而不斷改進(jìn)和優(yōu)化,而唯象模型在參數(shù)標(biāo)定好以后,成為一個(gè)封閉的體系,無(wú)法充分利用工程活動(dòng)中產(chǎn)生的海量數(shù)據(jù)資源;(2)復(fù)雜的唯象模型往往需要標(biāo)定多達(dá)十多個(gè),甚至數(shù)十個(gè)自由參數(shù),并且有些模型在有限元中實(shí)施較為復(fù)雜;然而數(shù)據(jù)驅(qū)動(dòng)的本構(gòu)模型不論考慮何種復(fù)雜的荷載條件,在有限元等連續(xù)方法中實(shí)施時(shí)均能夠做到標(biāo)準(zhǔn)化和模塊化.

    數(shù)據(jù)驅(qū)動(dòng)模型和理論驅(qū)動(dòng)模型具有不同的優(yōu)點(diǎn)和缺點(diǎn),數(shù)據(jù)驅(qū)動(dòng)本構(gòu)模型并不能替代理論本構(gòu)模型,兩種研究范式是相互補(bǔ)充和促進(jìn)的.一個(gè)典型的例子是,盡管離散元模擬可以得到代表性試樣中的剪應(yīng)力分量,在室內(nèi)試驗(yàn)條件下我們一般僅能測(cè)到試樣中的3 個(gè)主應(yīng)力和主應(yīng)變分量;深度學(xué)習(xí)模型在主應(yīng)力和主應(yīng)變數(shù)據(jù)集下也只能訓(xùn)練得到主應(yīng)力與主應(yīng)變空間的映射關(guān)系.一個(gè)完整的本構(gòu)模型需要進(jìn)一步將主應(yīng)力與主應(yīng)變空間的應(yīng)力關(guān)系變換到一般應(yīng)力?應(yīng)變空間當(dāng)中.由于顆粒材料在加載過(guò)程中主應(yīng)力和主應(yīng)變?cè)隽糠较虼嬖诜枪草S特征[31],因此純粹的坐標(biāo)變換方法將不可避免存在誤差,考慮顆粒材料主應(yīng)力和主應(yīng)變關(guān)系的非共軸特性是唯象模型和數(shù)據(jù)模型共同的挑戰(zhàn)[32].

    7 結(jié)論

    本文結(jié)合小應(yīng)變應(yīng)力?應(yīng)變理論關(guān)系和深度神經(jīng)網(wǎng)絡(luò)模型,通過(guò)有向圖發(fā)展了一種基于深度循環(huán)神經(jīng)網(wǎng)絡(luò)的顆粒材料應(yīng)力?應(yīng)變數(shù)據(jù)驅(qū)動(dòng)模型訓(xùn)練方法.通過(guò)常規(guī)三軸加卸載、等p真三軸加卸載和等b真三軸加載等多軸加載工況對(duì)模型進(jìn)行訓(xùn)練和測(cè)試.結(jié)果表明使用的深度神經(jīng)網(wǎng)絡(luò)模型能有效預(yù)測(cè)訓(xùn)練和驗(yàn)證過(guò)程從未見(jiàn)過(guò)的測(cè)試數(shù)據(jù),在僅含有常規(guī)三軸和等p真三軸一次加卸載循環(huán)的訓(xùn)練集中,深度學(xué)習(xí)模型有足夠的泛化能力預(yù)測(cè)到包含單調(diào)加載和多個(gè)加卸載循環(huán)的應(yīng)力響應(yīng),具有良好的內(nèi)插和外推預(yù)測(cè)能力.結(jié)果表明通過(guò)結(jié)合先驗(yàn)的應(yīng)力?應(yīng)變解析范式和經(jīng)驗(yàn)的數(shù)據(jù)驅(qū)動(dòng)范式,可以發(fā)展出適應(yīng)多種加載工況的顆粒材料主應(yīng)力?主應(yīng)變預(yù)測(cè)模型.

    猜你喜歡
    深度模型
    一半模型
    深度理解一元一次方程
    重要模型『一線(xiàn)三等角』
    重尾非線(xiàn)性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    深度觀察
    深度觀察
    深度觀察
    深度觀察
    3D打印中的模型分割與打包
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    成人黄色视频免费在线看| 精品一区二区三区四区五区乱码 | 黄色配什么色好看| 老熟女久久久| 热99国产精品久久久久久7| 久久久国产一区二区| 制服诱惑二区| 天天操日日干夜夜撸| 欧美变态另类bdsm刘玥| 日韩一区二区三区影片| 午夜日韩欧美国产| 国产精品欧美亚洲77777| 最黄视频免费看| 熟女电影av网| 久久毛片免费看一区二区三区| 999久久久国产精品视频| 女性被躁到高潮视频| 一二三四中文在线观看免费高清| 精品一区二区免费观看| 人成视频在线观看免费观看| 欧美激情 高清一区二区三区| 熟妇人妻不卡中文字幕| 不卡av一区二区三区| 制服人妻中文乱码| 性少妇av在线| 满18在线观看网站| 在线观看免费视频网站a站| 黑人巨大精品欧美一区二区蜜桃| 纯流量卡能插随身wifi吗| 97人妻天天添夜夜摸| 欧美中文综合在线视频| 欧美精品一区二区免费开放| 国产成人精品在线电影| 亚洲视频免费观看视频| 三上悠亚av全集在线观看| 狂野欧美激情性bbbbbb| 最近中文字幕高清免费大全6| 最近最新中文字幕免费大全7| 日韩伦理黄色片| 中文字幕人妻丝袜制服| 亚洲经典国产精华液单| 亚洲视频免费观看视频| 伊人久久国产一区二区| 国产成人aa在线观看| 中文字幕av电影在线播放| 亚洲内射少妇av| 美女主播在线视频| 男女免费视频国产| 亚洲国产精品999| av线在线观看网站| 久久久久国产精品人妻一区二区| xxxhd国产人妻xxx| 色哟哟·www| av福利片在线| 这个男人来自地球电影免费观看 | 蜜桃在线观看..| 久久人人97超碰香蕉20202| 国产精品久久久久久精品古装| 亚洲四区av| 久久韩国三级中文字幕| 97在线视频观看| 精品久久蜜臀av无| 日本色播在线视频| 亚洲欧美一区二区三区黑人 | 搡老乐熟女国产| 成人国产av品久久久| 婷婷色综合www| 精品少妇久久久久久888优播| 制服人妻中文乱码| 日本猛色少妇xxxxx猛交久久| 国产一区亚洲一区在线观看| 国产欧美日韩综合在线一区二区| 国产精品av久久久久免费| 欧美+日韩+精品| 日韩精品免费视频一区二区三区| 免费看av在线观看网站| 少妇猛男粗大的猛烈进出视频| 亚洲av电影在线观看一区二区三区| 一本—道久久a久久精品蜜桃钙片| 亚洲av成人精品一二三区| av在线观看视频网站免费| 天堂中文最新版在线下载| 亚洲经典国产精华液单| 老汉色av国产亚洲站长工具| 午夜福利一区二区在线看| 精品亚洲乱码少妇综合久久| 亚洲美女搞黄在线观看| 午夜精品国产一区二区电影| 在线观看www视频免费| 国产 精品1| 人妻人人澡人人爽人人| av在线播放精品| 久久ye,这里只有精品| 成人国语在线视频| 日韩伦理黄色片| 国产毛片在线视频| 亚洲图色成人| 色播在线永久视频| 国产精品一区二区在线不卡| 2018国产大陆天天弄谢| 一本色道久久久久久精品综合| 国产欧美日韩综合在线一区二区| 免费人妻精品一区二区三区视频| 99国产综合亚洲精品| 捣出白浆h1v1| 亚洲国产精品999| 久久av网站| 国产成人午夜福利电影在线观看| 亚洲国产看品久久| 久久久久久久亚洲中文字幕| 久久精品久久久久久久性| 97人妻天天添夜夜摸| 曰老女人黄片| 不卡av一区二区三区| av国产精品久久久久影院| 91精品伊人久久大香线蕉| 国产1区2区3区精品| 久久综合国产亚洲精品| 夫妻午夜视频| 国产欧美日韩综合在线一区二区| 五月伊人婷婷丁香| 亚洲精品久久久久久婷婷小说| 超色免费av| 国产白丝娇喘喷水9色精品| 高清视频免费观看一区二区| 亚洲人成网站在线观看播放| 赤兔流量卡办理| h视频一区二区三区| 新久久久久国产一级毛片| 国产又爽黄色视频| av有码第一页| 亚洲欧美一区二区三区黑人 | 老司机影院成人| 亚洲精品成人av观看孕妇| 免费看av在线观看网站| av在线老鸭窝| 蜜桃在线观看..| 午夜福利在线免费观看网站| 精品国产露脸久久av麻豆| 亚洲一区二区三区欧美精品| 看非洲黑人一级黄片| 精品久久久久久电影网| 国产男人的电影天堂91| 美女xxoo啪啪120秒动态图| 丰满迷人的少妇在线观看| 啦啦啦在线免费观看视频4| 99久久综合免费| 在现免费观看毛片| 欧美成人午夜精品| 国产欧美日韩一区二区三区在线| 免费黄网站久久成人精品| 丝袜喷水一区| 考比视频在线观看| 99国产精品免费福利视频| 日韩熟女老妇一区二区性免费视频| 日本免费在线观看一区| 老鸭窝网址在线观看| 亚洲美女黄色视频免费看| 极品人妻少妇av视频| 一级片'在线观看视频| 亚洲av电影在线观看一区二区三区| 亚洲欧美一区二区三区久久| 亚洲国产精品成人久久小说| 亚洲欧美日韩另类电影网站| 国产精品蜜桃在线观看| 亚洲第一av免费看| 国产一区二区三区综合在线观看| 免费高清在线观看日韩| 汤姆久久久久久久影院中文字幕| 精品久久久久久电影网| 成人18禁高潮啪啪吃奶动态图| av国产精品久久久久影院| 亚洲综合精品二区| 亚洲五月色婷婷综合| 韩国高清视频一区二区三区| 如何舔出高潮| 少妇被粗大的猛进出69影院| 97在线人人人人妻| 亚洲国产毛片av蜜桃av| 国产探花极品一区二区| 成人亚洲欧美一区二区av| 母亲3免费完整高清在线观看 | 国产精品三级大全| 丝袜美足系列| 国产精品99久久99久久久不卡 | 狠狠精品人妻久久久久久综合| 秋霞伦理黄片| 国产欧美亚洲国产| 免费观看在线日韩| 中文欧美无线码| 亚洲综合色惰| 亚洲欧洲精品一区二区精品久久久 | 久久精品夜色国产| 亚洲少妇的诱惑av| 久久久欧美国产精品| 亚洲国产毛片av蜜桃av| www.av在线官网国产| 久久久久久久国产电影| 亚洲精品一区蜜桃| 免费黄频网站在线观看国产| 日韩电影二区| 蜜桃在线观看..| 精品人妻熟女毛片av久久网站| av有码第一页| 美女脱内裤让男人舔精品视频| 免费高清在线观看视频在线观看| av网站免费在线观看视频| 美女高潮到喷水免费观看| 久久久久久久久久久免费av| 国产免费福利视频在线观看| 一级片免费观看大全| av一本久久久久| 国产乱来视频区| 十八禁网站网址无遮挡| 免费高清在线观看视频在线观看| 国产精品欧美亚洲77777| 春色校园在线视频观看| 晚上一个人看的免费电影| 精品国产一区二区三区四区第35| 国产免费福利视频在线观看| av免费在线看不卡| 国语对白做爰xxxⅹ性视频网站| 丰满饥渴人妻一区二区三| av视频免费观看在线观看| 亚洲欧美一区二区三区黑人 | 不卡av一区二区三区| 久久精品久久久久久久性| 国产成人aa在线观看| 少妇熟女欧美另类| 久久久久久免费高清国产稀缺| 欧美成人午夜精品| 蜜桃国产av成人99| 国产亚洲精品第一综合不卡| 亚洲欧美色中文字幕在线| 久久久久久久亚洲中文字幕| 91成人精品电影| 国产日韩欧美亚洲二区| 国产av一区二区精品久久| 亚洲人成电影观看| 精品一区在线观看国产| 国产精品久久久久成人av| 韩国高清视频一区二区三区| 80岁老熟妇乱子伦牲交| 精品国产乱码久久久久久男人| 2022亚洲国产成人精品| 亚洲av福利一区| 人妻人人澡人人爽人人| 丝瓜视频免费看黄片| 久久久国产精品麻豆| 日韩精品有码人妻一区| 看免费av毛片| 国产精品嫩草影院av在线观看| 亚洲美女搞黄在线观看| 免费在线观看完整版高清| 亚洲av日韩在线播放| 日韩不卡一区二区三区视频在线| av一本久久久久| 国产精品嫩草影院av在线观看| 亚洲伊人色综图| 赤兔流量卡办理| 免费女性裸体啪啪无遮挡网站| 制服诱惑二区| 黑人猛操日本美女一级片| 最近手机中文字幕大全| 精品99又大又爽又粗少妇毛片| 亚洲男人天堂网一区| 欧美精品一区二区大全| 各种免费的搞黄视频| 9191精品国产免费久久| 亚洲国产欧美日韩在线播放| 少妇被粗大的猛进出69影院| 男女啪啪激烈高潮av片| 久久久久精品久久久久真实原创| 69精品国产乱码久久久| av又黄又爽大尺度在线免费看| 欧美日韩亚洲国产一区二区在线观看 | 色播在线永久视频| 少妇熟女欧美另类| 极品少妇高潮喷水抽搐| 看免费成人av毛片| 街头女战士在线观看网站| 在线观看国产h片| av线在线观看网站| a级片在线免费高清观看视频| 电影成人av| 国产精品欧美亚洲77777| 国产片特级美女逼逼视频| 国产又爽黄色视频| 在线看a的网站| 一级毛片 在线播放| 9191精品国产免费久久| 国产精品蜜桃在线观看| 亚洲精品久久成人aⅴ小说| 中文字幕另类日韩欧美亚洲嫩草| 国产黄色免费在线视频| 菩萨蛮人人尽说江南好唐韦庄| 综合色丁香网| 亚洲av免费高清在线观看| 涩涩av久久男人的天堂| 国产成人精品在线电影| 18在线观看网站| 又粗又硬又长又爽又黄的视频| 热re99久久国产66热| 妹子高潮喷水视频| 成人午夜精彩视频在线观看| 欧美日韩一区二区视频在线观看视频在线| 又大又黄又爽视频免费| www.熟女人妻精品国产| 亚洲国产色片| 午夜免费观看性视频| 亚洲一码二码三码区别大吗| 爱豆传媒免费全集在线观看| 日日摸夜夜添夜夜爱| 国产精品久久久av美女十八| 欧美激情高清一区二区三区 | 香蕉丝袜av| 人人妻人人澡人人看| 中文字幕av电影在线播放| 亚洲精品成人av观看孕妇| 美女午夜性视频免费| 99国产综合亚洲精品| 五月伊人婷婷丁香| 青春草亚洲视频在线观看| 97人妻天天添夜夜摸| 午夜福利视频在线观看免费| videosex国产| 亚洲天堂av无毛| 国产人伦9x9x在线观看 | 成人亚洲欧美一区二区av| 99热国产这里只有精品6| av国产精品久久久久影院| 国产一级毛片在线| 女性生殖器流出的白浆| av福利片在线| 日日摸夜夜添夜夜爱| 另类精品久久| kizo精华| 哪个播放器可以免费观看大片| 观看美女的网站| 黄色 视频免费看| 熟女电影av网| 成人国语在线视频| 老司机亚洲免费影院| 久久久a久久爽久久v久久| 国产精品av久久久久免费| 亚洲伊人久久精品综合| 99久久精品国产国产毛片| 久久国产精品男人的天堂亚洲| 男女啪啪激烈高潮av片| 国产成人av激情在线播放| 99久久人妻综合| 欧美xxⅹ黑人| 亚洲精品国产一区二区精华液| 亚洲国产最新在线播放| 黄色 视频免费看| 久久精品国产亚洲av涩爱| 一区二区三区精品91| 久久人人爽av亚洲精品天堂| 久久久久久人人人人人| 欧美激情高清一区二区三区 | 亚洲一区二区三区欧美精品| 国产亚洲午夜精品一区二区久久| 人妻系列 视频| 看免费成人av毛片| 成人国产av品久久久| 岛国毛片在线播放| 少妇人妻 视频| 欧美 亚洲 国产 日韩一| 午夜免费男女啪啪视频观看| 亚洲精品,欧美精品| 精品少妇久久久久久888优播| av女优亚洲男人天堂| 一个人免费看片子| 久久久久国产网址| freevideosex欧美| 国产精品欧美亚洲77777| av国产久精品久网站免费入址| 精品99又大又爽又粗少妇毛片| 久久亚洲国产成人精品v| 丝袜美腿诱惑在线| 99热网站在线观看| av电影中文网址| 丁香六月天网| 天堂中文最新版在线下载| 又粗又硬又长又爽又黄的视频| 国产乱来视频区| 国产在线免费精品| 各种免费的搞黄视频| a级毛片黄视频| 高清视频免费观看一区二区| 日韩制服骚丝袜av| 丝袜脚勾引网站| 99re6热这里在线精品视频| 一区在线观看完整版| 99九九在线精品视频| 国产欧美日韩一区二区三区在线| 亚洲在久久综合| 如日韩欧美国产精品一区二区三区| 午夜老司机福利剧场| 天天躁狠狠躁夜夜躁狠狠躁| 熟女电影av网| 十分钟在线观看高清视频www| 免费日韩欧美在线观看| 亚洲精品久久成人aⅴ小说| 国产精品国产三级专区第一集| 黑人欧美特级aaaaaa片| 男女免费视频国产| 亚洲欧洲国产日韩| 夜夜骑夜夜射夜夜干| 午夜免费观看性视频| 欧美精品一区二区免费开放| 男男h啪啪无遮挡| 精品视频人人做人人爽| 欧美日韩成人在线一区二区| av在线播放精品| 国产xxxxx性猛交| 亚洲激情五月婷婷啪啪| 黑人欧美特级aaaaaa片| 午夜精品国产一区二区电影| 人体艺术视频欧美日本| 亚洲内射少妇av| 欧美成人午夜精品| 免费大片黄手机在线观看| 久久久久久伊人网av| 午夜激情av网站| 中文字幕人妻丝袜一区二区 | 午夜福利在线免费观看网站| 国产成人精品婷婷| 少妇熟女欧美另类| 黄色一级大片看看| av片东京热男人的天堂| 老司机影院成人| 国产毛片在线视频| 久久这里只有精品19| 在线 av 中文字幕| 中文字幕人妻丝袜一区二区 | 欧美日韩视频精品一区| 亚洲国产最新在线播放| 久久久精品免费免费高清| av福利片在线| 我要看黄色一级片免费的| 久久这里有精品视频免费| av免费观看日本| av电影中文网址| 亚洲精品美女久久久久99蜜臀 | 亚洲国产成人一精品久久久| 午夜福利网站1000一区二区三区| 午夜免费鲁丝| 久久国产精品男人的天堂亚洲| 久久精品久久久久久噜噜老黄| 99久国产av精品国产电影| 国产亚洲午夜精品一区二区久久| 国产欧美日韩综合在线一区二区| 亚洲精品日本国产第一区| 中文字幕人妻丝袜一区二区 | 少妇人妻 视频| 亚洲精品成人av观看孕妇| 最近最新中文字幕大全免费视频 | 亚洲国产精品国产精品| 大陆偷拍与自拍| 爱豆传媒免费全集在线观看| 久久久久精品性色| 精品酒店卫生间| 七月丁香在线播放| 国产黄色视频一区二区在线观看| 国产不卡av网站在线观看| 成人国产麻豆网| a级毛片黄视频| 在线观看免费视频网站a站| 婷婷色麻豆天堂久久| 大片免费播放器 马上看| 免费看不卡的av| 久久久久视频综合| 中文字幕亚洲精品专区| 一级片免费观看大全| 一级毛片我不卡| 不卡视频在线观看欧美| 国产精品 欧美亚洲| 精品亚洲成a人片在线观看| 国产 一区精品| 亚洲欧美精品综合一区二区三区 | 亚洲一码二码三码区别大吗| av网站在线播放免费| 免费黄色在线免费观看| 日韩,欧美,国产一区二区三区| 精品国产乱码久久久久久小说| 久久韩国三级中文字幕| 国产亚洲av片在线观看秒播厂| 国产在线视频一区二区| 91aial.com中文字幕在线观看| 日本欧美国产在线视频| 国产探花极品一区二区| 久久久久视频综合| 中文字幕亚洲精品专区| a级毛片在线看网站| av卡一久久| 国产成人午夜福利电影在线观看| 日韩制服骚丝袜av| 日韩一区二区三区影片| 寂寞人妻少妇视频99o| 久久久国产一区二区| 国产一级毛片在线| 色网站视频免费| 亚洲国产精品一区三区| 新久久久久国产一级毛片| 一区在线观看完整版| 欧美精品亚洲一区二区| 欧美日韩综合久久久久久| 美女xxoo啪啪120秒动态图| 日日啪夜夜爽| 国产黄色免费在线视频| 最近手机中文字幕大全| 尾随美女入室| 自拍欧美九色日韩亚洲蝌蚪91| 久久久久久人人人人人| 亚洲一区二区三区欧美精品| av不卡在线播放| 国产一区二区激情短视频 | 亚洲av男天堂| 成人二区视频| 国产精品 国内视频| 婷婷色综合www| 五月开心婷婷网| 久久综合国产亚洲精品| 精品一品国产午夜福利视频| av片东京热男人的天堂| 欧美日韩亚洲高清精品| 亚洲国产欧美日韩在线播放| 一级片免费观看大全| 亚洲av免费高清在线观看| 黄色配什么色好看| 久久人人爽av亚洲精品天堂| 国产免费又黄又爽又色| 91国产中文字幕| 中文乱码字字幕精品一区二区三区| 精品人妻熟女毛片av久久网站| 国产男人的电影天堂91| 高清在线视频一区二区三区| 人人妻人人澡人人看| 熟妇人妻不卡中文字幕| 国产熟女午夜一区二区三区| 下体分泌物呈黄色| 满18在线观看网站| 久久国内精品自在自线图片| 免费高清在线观看日韩| 毛片一级片免费看久久久久| 久久久国产精品麻豆| 美女福利国产在线| 大香蕉久久网| 精品一区在线观看国产| 亚洲精品一区蜜桃| 国产熟女午夜一区二区三区| 麻豆乱淫一区二区| 男女边摸边吃奶| 免费黄色在线免费观看| a级毛片在线看网站| 国产精品免费视频内射| 男女啪啪激烈高潮av片| 水蜜桃什么品种好| 久久久久久久久久久久大奶| 一本大道久久a久久精品| 国产有黄有色有爽视频| 美国免费a级毛片| 在线天堂最新版资源| 亚洲成人手机| 精品人妻在线不人妻| 十分钟在线观看高清视频www| 亚洲av在线观看美女高潮| 精品人妻偷拍中文字幕| 丝袜美足系列| 免费高清在线观看日韩| 午夜老司机福利剧场| 你懂的网址亚洲精品在线观看| 99九九在线精品视频| 亚洲人成77777在线视频| 大片免费播放器 马上看| 老司机影院毛片| 9191精品国产免费久久| 黑丝袜美女国产一区| 99久久精品国产国产毛片| 叶爱在线成人免费视频播放| 卡戴珊不雅视频在线播放| 亚洲精品美女久久久久99蜜臀 | 久久久国产欧美日韩av| 久久久精品国产亚洲av高清涩受| 香蕉精品网在线| 少妇人妻 视频| 纯流量卡能插随身wifi吗| 成人手机av| 欧美日韩亚洲高清精品| 国产 一区精品| 免费看av在线观看网站| 免费观看a级毛片全部| 国产精品蜜桃在线观看| 国产熟女午夜一区二区三区| 伊人久久大香线蕉亚洲五| 亚洲美女视频黄频| 伊人久久国产一区二区| 少妇人妻久久综合中文| 欧美av亚洲av综合av国产av | 国产97色在线日韩免费| 男人添女人高潮全过程视频| 免费av中文字幕在线| 免费看不卡的av| 国产av国产精品国产| 欧美最新免费一区二区三区| 精品久久久精品久久久| 亚洲欧美一区二区三区黑人 | 熟女电影av网| 五月天丁香电影| 最新的欧美精品一区二区| 精品国产超薄肉色丝袜足j| 欧美日韩一级在线毛片| av在线app专区| 国产男女内射视频|