龍 旭，毛明暉，盧昶衡，蘇天雄，賈馮睿

（1.西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院，西安 710072；2.遼寧石油化工大學(xué)土木工程學(xué)院，撫順 113001；3.浙江清華長(zhǎng)三角研究院，嘉興 314006）

作為一種質(zhì)優(yōu)價(jià)廉的建筑材料，混凝土材料已經(jīng)廣泛地運(yùn)用在日常的生產(chǎn)生活。然而在突發(fā)自然災(zāi)害以及軍事防護(hù)工程中，需要更加準(zhǔn)確且高效地確定混凝土材料的抗沖擊性能，因此需要進(jìn)一步了解高應(yīng)變率下混凝土材料力學(xué)性能及其破壞行為規(guī)律。傳統(tǒng)高應(yīng)變率測(cè)試實(shí)驗(yàn)，主要包括落錘試驗(yàn)、液壓伺服試驗(yàn)和霍普金森桿沖擊實(shí)驗(yàn)［1］。相比落錘試驗(yàn)和液壓伺服試驗(yàn)，霍普金森桿沖擊方法可以實(shí)現(xiàn)動(dòng)態(tài)加載條件下更加精確的應(yīng)變率控制以及更廣的應(yīng)變率測(cè)量范圍［2］，因此被廣泛地運(yùn)用于各種材料動(dòng)態(tài)性能測(cè)試［3?6］?；炷敛牧媳举|(zhì)上是一種多相的非均勻材料，但從宏觀角度而言，當(dāng)混凝土尺寸為最大骨料尺寸的4～5 倍時(shí)，混凝土可被認(rèn)為是均勻材料［7］。受限于集料的大小，混凝土霍普金森桿沖擊實(shí)驗(yàn)需要滿(mǎn)足一定的幾何尺寸要求，因此所用設(shè)備的入射桿及透射桿直徑需要相應(yīng)增大［8］。目前已有研究主要利用大直徑的分離式霍普金森壓桿實(shí)驗(yàn)裝置測(cè)量混凝土材料的動(dòng)態(tài)力學(xué)性能，同時(shí)運(yùn)用有限元方法對(duì)其進(jìn)行數(shù)值仿真，從而可以更加高效地揭示主導(dǎo)混凝土類(lèi)材料高應(yīng)變率下變形特性和失效規(guī)律。

當(dāng)前混凝土類(lèi)材料大直徑分離式霍普金森壓桿（Split Hopkinson pressure bar，SHPB）實(shí)驗(yàn)主要存在以下幾個(gè)問(wèn)題：端面摩擦效應(yīng)［9?11］、骨料尺寸［12?16］、慣性效應(yīng)［17?22］和溫度效應(yīng)［23］等。SHPB 實(shí)驗(yàn)研究的推進(jìn)，必須建立清晰準(zhǔn)確的入射波、反射波、透射波波形，因此如何準(zhǔn)確高效地確定波形，是混凝土類(lèi)材料動(dòng)態(tài)性能研究的關(guān)鍵所在。一般而言，SHPB 實(shí)驗(yàn)前需要進(jìn)行無(wú)試件空撞試驗(yàn)，確保試件、入射桿和透射桿位于同一直線。為保證桿件系統(tǒng)的一致性，還需要在標(biāo)準(zhǔn)試件上進(jìn)行試撞擊。因此，SHPB 實(shí)驗(yàn)準(zhǔn)備階段較為復(fù)雜［24］。此外，往往需要大量人工重復(fù)工作才能合理選取應(yīng)力波起始位置［25］，進(jìn)而進(jìn)行二波或三波對(duì)齊之后實(shí)驗(yàn)數(shù)據(jù)的處理，部分明顯失真數(shù)據(jù)仍需人工剔除，人為因素也一定程度上增大了實(shí)驗(yàn)誤差。因此，混凝土類(lèi)材料SHPB 實(shí)驗(yàn)亟待更為省時(shí)且更為客觀的波形確定方法。

混凝土作為一種多相的脆性復(fù)合材料，需要避免SHPB 實(shí)驗(yàn)中應(yīng)力波傳播時(shí)彌散效應(yīng)所致波形振蕩，且混凝土材料發(fā)生破壞時(shí)應(yīng)變較小，需要更加準(zhǔn)確地控制入射波的波形及作用時(shí)長(zhǎng)，上述客觀因素均會(huì)影響SHPB 實(shí)驗(yàn)測(cè)量數(shù)據(jù)的準(zhǔn)確性。宋博等［26］對(duì)比SHPB 實(shí)驗(yàn)中移波處理中簡(jiǎn)單時(shí)間法和不均勻時(shí)間應(yīng)力法，提供不同解耦方法，并且指出兩種方法適用材料范圍。Wang 等［4］利用SHPB系統(tǒng)測(cè)量不同含量鋼纖維混凝土動(dòng)態(tài)強(qiáng)度，研究鋼纖維對(duì)于動(dòng)態(tài)強(qiáng)度和破壞機(jī)理基本理論。胡亮亮等［27］搜集大量混凝土類(lèi)SHPB 高應(yīng)變率實(shí)驗(yàn)數(shù)據(jù)，認(rèn)為脆性材料直接用短平臺(tái)波表征其應(yīng)變率是不合理的，應(yīng)該選用其全段平均應(yīng)變率的1.38 倍。胡金生等［28］認(rèn)為大直徑SHPB 混凝土實(shí)驗(yàn)?zāi)芙档筒ㄐ握袷帉?duì)于數(shù)據(jù)擬合的影響，但是會(huì)限制混凝土類(lèi)脆性材料的高應(yīng)變率的測(cè)量。梁書(shū)鋒［29］通過(guò)大量巖石SHPB 實(shí)驗(yàn)，認(rèn)為脆性材料恒應(yīng)變率加載難度較大，指出三維應(yīng)力增量表達(dá)形式容易表示沖擊峰值前的動(dòng)態(tài)力學(xué)性能。王魯明等［9］針對(duì)脆性材料SHPB 實(shí)驗(yàn)討論萬(wàn)向頭、波形整形器、異形炮頭等方法降低波形振蕩，提高實(shí)驗(yàn)的準(zhǔn)確性。楊陽(yáng)等［30］討論不同霍普金森實(shí)驗(yàn)撞擊桿幾何形狀對(duì)波形的影響，針對(duì)不同材料特性合理選用撞擊桿，消除材料的Pochhammer?Chree 效應(yīng)。陳滔等［31］給出混凝土SHPB 實(shí)驗(yàn)壓縮性的改進(jìn)Gorham 慣性效應(yīng)公式，該公式可以更好地修正慣性效應(yīng)從而提高實(shí)驗(yàn)精度。

然而學(xué)者們只能根據(jù)實(shí)驗(yàn)波形進(jìn)行修正，很難做到新波形預(yù)測(cè)實(shí)現(xiàn)材料動(dòng)態(tài)力學(xué)性能預(yù)測(cè)。人工神經(jīng)網(wǎng)絡(luò)具有自學(xué)習(xí)的特點(diǎn)，適合處理非線性問(wèn)題，其中多層神經(jīng)元系統(tǒng)具有在理論上可以預(yù)測(cè)任何非線性的模型，適合預(yù)測(cè)SHPB 實(shí)驗(yàn)中反射波和透射波，從而實(shí)現(xiàn)混凝土類(lèi)材料的動(dòng)態(tài)力學(xué)性能預(yù)測(cè)。當(dāng)前的深度學(xué)習(xí)模型具有多種結(jié)構(gòu)，包括反向傳播神經(jīng)網(wǎng)絡(luò)（Back propagation，BP）、卷積神經(jīng)網(wǎng)絡(luò)（Convolutional neural network，CNN）、循環(huán)神經(jīng)網(wǎng)絡(luò)（Recurrent neural network，RNN）和長(zhǎng)短期記憶網(wǎng)絡(luò)（Long and short term memory，LSTM）模型等人工神經(jīng)網(wǎng)絡(luò)，上述神經(jīng)網(wǎng)絡(luò)均有具有自學(xué)習(xí)的特點(diǎn)。嚴(yán)曉明等［32］認(rèn)為BP 神經(jīng)網(wǎng)絡(luò)方法可根據(jù)誤差進(jìn)行所有節(jié)點(diǎn)權(quán)值的調(diào)整，相比較其他算法，該算法更加成熟穩(wěn)定，具有良好的適應(yīng)性。雖然學(xué)習(xí)速度較慢，但能夠降低誤差，能夠?qū)崿F(xiàn)輸入層和輸出層之間參數(shù)的非線性映射。Wang 等［33］將BP 神經(jīng)網(wǎng)絡(luò)技術(shù)運(yùn)用于波形處理，通過(guò)相關(guān)技術(shù)降低波的振蕩，實(shí)現(xiàn)實(shí)驗(yàn)波形和理想波形映射關(guān)系，從而提高了應(yīng)力?應(yīng)變響應(yīng)曲線的準(zhǔn)確性。徐明喬等［34］利用BP 神經(jīng)網(wǎng)絡(luò)技術(shù)研究了高聚合物SHPB 實(shí)驗(yàn)非線性率相關(guān)本構(gòu)行為。朱勵(lì)等［35］針對(duì)波形的隱式關(guān)系，利用反向誤差人工智能網(wǎng)絡(luò)技術(shù)，對(duì)波的彌散效應(yīng)進(jìn)行了修正，且修正結(jié)果具有良好的適用性。

綜上所述，已有研究中較多利用數(shù)值模擬和實(shí)驗(yàn)確定混凝土類(lèi)材料高應(yīng)變率下變形行為，尚未有運(yùn)用機(jī)器學(xué)習(xí)算法預(yù)測(cè)混凝土類(lèi)材料高應(yīng)變率下應(yīng)力?應(yīng)變力學(xué)性能的研究。因此，本文提出一種將SHPB 實(shí)驗(yàn)仿真模擬與BP 人工神經(jīng)網(wǎng)絡(luò)技術(shù)相結(jié)合的新方法，然后將其成功應(yīng)用于混凝土類(lèi)脆性材料動(dòng)態(tài)相應(yīng)研究，利用Zhang 等［36］和Li 等［37］的實(shí)驗(yàn)數(shù)據(jù)，建立有限元模型，對(duì)比SHPB 實(shí)驗(yàn)驗(yàn)證數(shù)值模擬準(zhǔn)確性。將數(shù)值模擬的入射波作為機(jī)器學(xué)習(xí)訓(xùn)練樣本的輸入層，反射波和透射波作為輸出層，建立兩組波形預(yù)測(cè)模型，從而預(yù)測(cè)材料動(dòng)態(tài)力學(xué)性能。該方法能夠代替人工量大且耗時(shí)有限元仿真、分析，處理等工作，建立的學(xué)習(xí)模型為人工智能預(yù)測(cè)混凝土類(lèi)材料動(dòng)態(tài)力學(xué)性能提供理論基礎(chǔ)和經(jīng)驗(yàn)借鑒。

1 SHPB 數(shù)值仿真

1.1 SHPB 實(shí)驗(yàn)介紹

SHPB 實(shí)驗(yàn)技術(shù)以?xún)蓚€(gè)基本假定［38］為基礎(chǔ)。一是彈性桿的一維應(yīng)力波假定，即壓桿在實(shí)驗(yàn)過(guò)程中始終處于彈性應(yīng)變范圍（壓桿剛度遠(yuǎn)大于試件剛度的情況下此假定可滿(mǎn)足）；二是試件中的應(yīng)力、應(yīng)變沿試件軸向方向均勻分布，即均勻性假定（試件的長(zhǎng)度遠(yuǎn)小于應(yīng)力波波長(zhǎng)的情況下此假定可滿(mǎn)足）。滿(mǎn)足這兩個(gè)基本假定，才能得到有效的實(shí)驗(yàn)結(jié)果。

常規(guī)的霍普金森桿如圖1 所示，主要包括撞擊桿、入射桿和透射桿等桿件。3 種不同的波形對(duì)應(yīng)于SHPB 實(shí)驗(yàn)中分別是入射波εi、反射波εr和透射波εt。入射波和反射波通過(guò)入射桿上的應(yīng)變片測(cè)得，透射波通過(guò)透射桿上的應(yīng)變片測(cè)得?；羝战鹕瓕?shí)驗(yàn)中需要通過(guò)移波處理，將3 種波形的起始位置移動(dòng)至同一個(gè)時(shí)間軸上，便于后期的數(shù)據(jù)處理。在數(shù)據(jù)處理中常用的方法為二波法和三波法，其中三波法更加適用于混凝土材料［39］，本文選用三波法處理應(yīng)變率，應(yīng)力?應(yīng)變之間關(guān)系。

圖1 SHPB 裝置示意圖Fig.1 Schematic diagram of SHPB setup

式中：C0代表?xiàng)U件波速，Eb代表?xiàng)U件彈性模量，ρ代表?xiàng)U件密度，Ls代表試件長(zhǎng)度，代表試樣應(yīng)變率，v1和v2分別表示入射桿和試樣以及透射桿與試樣接觸端面速度，εs代表試樣應(yīng)變，t代表整個(gè)沖擊過(guò)程時(shí)間，Ab代表?xiàng)U件橫截面面積，As代表試樣橫截面面積，σ1、σ2和σs分別代表試樣前后端面的應(yīng)力和試樣應(yīng)力。

1.2 SHPB 沖擊壓縮有限元計(jì)算

1.2.1 Drucker?Prager 線性擴(kuò)展理論

基于塑性和超彈性的本構(gòu)模型被廣泛應(yīng)用于混凝土類(lèi)材料SHPB 的沖擊數(shù)值模擬中。Drucker?Prager 模型可以描述因?yàn)閮?nèi)摩擦引起的靜水壓力對(duì)于單元體的綜合影響，ABAQUS 中線性Drucker?Prager 模型被廣泛運(yùn)用至混凝土類(lèi)材料的動(dòng)態(tài)力學(xué)行為模擬。本文選用的是線性屈服面子午線，對(duì)應(yīng)的Drucker?Prager 模型的屈服面函數(shù)可寫(xiě)為

式中：t為偏應(yīng)力的度量值；r為π平面上的半徑；J2為第二偏應(yīng)力不變量；J3為第三偏應(yīng)力不變量；p為等效圍壓效應(yīng)下的壓應(yīng)力，p=-l1/3；l1指的是第一應(yīng)力不變量；K為三軸拉伸屈服應(yīng)力與三軸壓縮屈服應(yīng)力之比，因此該值控制著屈服面對(duì)中間主應(yīng)力值的依賴(lài)性，當(dāng)K＝1時(shí)，t=q，屈服面在π平面上為Von Mises 圓上，為了保證屈服面外凸，要求0.778

當(dāng)材料的單軸受壓屈服的屈服應(yīng)力為σc時(shí)，黏聚力可表達(dá)為

Drucker?Prager 模型流動(dòng)性規(guī)律定義式為

式中：φ為p?t平面內(nèi)的膨脹角，p?t平面如圖2所示。

圖2 p?t 平面示意圖Fig.2 Schematic diagram of p?t plane

1.2.2 有限元模型建立和驗(yàn)證

本文有效性驗(yàn)證模擬Zhang 等［36］實(shí)驗(yàn)，該實(shí)驗(yàn)數(shù)值模擬中，用入射應(yīng)力波取代撞擊桿降低入射波振蕩。如表1 所示，SHPB 實(shí)驗(yàn)中入射桿長(zhǎng)度為3.0 m，入射桿應(yīng)變片在距離試樣與入射桿接觸面1.5 m 處，透射桿長(zhǎng)度為2.0 m，透射桿應(yīng)變片在距離試樣與透射桿接觸面0.5 m 處，入射桿和透射桿直徑均為100 mm，材料均為鋼。試樣的長(zhǎng)度21.5 mm，直徑為740 mm，基體材料為砂漿。如圖3 所示，本文采用ABAQUS 2017 商業(yè)軟件對(duì)上述模型進(jìn)行有限元仿真，取全結(jié)構(gòu)的1/4 進(jìn)行數(shù)值模擬，以降低計(jì)算時(shí)間。根據(jù)實(shí)際工況選擇合適的單元類(lèi)型，本文中針對(duì)的是混凝土類(lèi)材料下高應(yīng)變率力學(xué)性能實(shí)驗(yàn)屬于典型的動(dòng)力學(xué)分析，選擇顯示動(dòng)力學(xué)模塊中Explicit 且為3D Stress 單元，所有單元為C3D8 單元。入射桿和透射桿的布種間距為10 mm，試樣上布種間距為1 mm，入射桿共計(jì)4 848 單元，試樣共計(jì)28 710 單元，透射桿共計(jì)7 428 單元。自由端面入射應(yīng)力波與時(shí)間關(guān)系曲線見(jiàn)式（12），徑向邊界條件為對(duì)稱(chēng)固定以模擬全結(jié)構(gòu)中沖擊過(guò)程試樣徑向?qū)ΨQ(chēng)分布。仿真完成后，提取入射桿和透射桿上應(yīng)變片位置處的應(yīng)變與時(shí)間關(guān)系，選取合適的波形起始點(diǎn)和結(jié)束點(diǎn)，將應(yīng)變片上波形劃分成入射波、反射波和透射波。

表1 SHPB 實(shí)驗(yàn)裝置和試件幾何尺寸Table 1 Geometry of tested material and bars of SHPB

圖3 SHPB 有限元模型Fig.3 Finite element model of SHPB

在SHPB 實(shí)驗(yàn)當(dāng)中，試樣與端面之間摩擦已經(jīng)引起了人們的足夠重視，Li 等［37］開(kāi)始研究砂漿和混凝土材料中端面接觸摩擦的影響，因?yàn)槿鄙俣嗣婺Σ列?yīng)的具體數(shù)值，根據(jù)現(xiàn)有的石灰?guī)rSHPB 實(shí)驗(yàn)數(shù)值模擬，進(jìn)行端面摩擦估值0.3［40］。故對(duì)于模型中的接觸面接觸條件設(shè)置為面面接觸，法向接觸設(shè)置為硬接觸，切向接觸設(shè)置為罰函數(shù)，取摩擦因數(shù)估值為0.2。

在實(shí)際的數(shù)值模擬過(guò)程中，應(yīng)當(dāng)詳細(xì)地標(biāo)定相關(guān)的Drucker?Prager 的參數(shù)，根據(jù)前人［41?43］研究得知基體砂漿材料的β值介于40°至60°之間。通過(guò)Li 等［37］試算，在SHPB 實(shí)驗(yàn)中β=40°能夠很好地符合實(shí)驗(yàn)的數(shù)值結(jié)果，因此本文中的β=40°，同時(shí)Grote 等［44］發(fā)現(xiàn)φ在基體材料砂漿中影響有限，因此本文選用的φ=β=40°。Park 等［45］認(rèn)為砂漿和混凝土材料是一種K值不敏感性材料，故本文中的K選取為1。具體材料參數(shù)見(jiàn)表2。

表2 材料參數(shù)Table 2 Parameters of finite element model

入射桿應(yīng)力波輸入要保證能夠在較快時(shí)間內(nèi)達(dá)到指定的應(yīng)變率，從而保持試件前后端面的應(yīng)力平衡，將動(dòng)態(tài)測(cè)量轉(zhuǎn)換成準(zhǔn)靜態(tài)測(cè)量。常規(guī)的矩形波在測(cè)試巖石類(lèi)的混凝土試樣時(shí)，入射波的上升波形斜率很高，據(jù)應(yīng)力波的初等理論可知入射波可以分解為不同傳播速度和頻率的諧波分量組合而成，這些不同速度的波各自傳播使得波形出現(xiàn)明顯振蕩，入射波的振蕩必然導(dǎo)致反射波和透射波的振蕩從而影響整體試驗(yàn)效果。同時(shí)振蕩波相對(duì)于混凝土材料強(qiáng)度的相對(duì)值可能已經(jīng)超過(guò)混凝土強(qiáng)度本身，使得材料出現(xiàn)小幅度的加載和卸載，上述不利因素進(jìn)一步加劇波的振蕩，故需要對(duì)于入射波波頭緩和處理，以滿(mǎn)足有足夠的時(shí)間保持整體試樣達(dá)到應(yīng)力平衡。

本文選取直接在入射桿自由面添加入射應(yīng)力波法，入射桿自由端面的入射應(yīng)力和時(shí)間的關(guān)系見(jiàn)式（12）［36］，其中σ為入射波的應(yīng)力值。如圖4 所示，本文所采用入射波峰值應(yīng)力σi分別為60、80、100 和160 MPa，t對(duì)應(yīng)的是時(shí)間軸，單位為μs，應(yīng)力加載為500 μs，模擬總時(shí)長(zhǎng)1 400 μs。

圖4 不同入射波峰值所引起的入射應(yīng)力波形Fig.4 Incident stress waves induced by different inci?dent peak stresses

圖5（a）為160 MPa 峰值應(yīng)力SHPB 實(shí)驗(yàn)?zāi)M典型入射桿和透射桿圖，圖5（b）為60 MPa 峰值應(yīng)力前后端應(yīng)力對(duì)比圖，可以看出在整個(gè)過(guò)程中出現(xiàn)較好吻合，即保證材料在實(shí)驗(yàn)過(guò)程中的動(dòng)態(tài)平衡。圖5（c）為4 組不同的應(yīng)力波下的SHPB 實(shí)驗(yàn)應(yīng)力?應(yīng)變圖，圖5（d）為本文計(jì)算DIF 與文獻(xiàn)［36］實(shí)測(cè)DIF 的比較，其中DIF 為動(dòng)態(tài)峰值應(yīng)力與準(zhǔn)靜止態(tài)峰值應(yīng)力比值。可以看出，本文所建立SHPB 仿真預(yù)測(cè)結(jié)果可以較好地吻合實(shí)驗(yàn)結(jié)果，可以合理地復(fù)現(xiàn)此類(lèi)材料在高應(yīng)變率下變形特性和失效機(jī)理。

圖5 數(shù)值模擬所得動(dòng)態(tài)響應(yīng)Fig.5 Dynamic responses by numerical simulations

2 BP 神經(jīng)網(wǎng)絡(luò)及其預(yù)測(cè)結(jié)果

2.1 BP 神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)

反向傳播（Back propagation，BP）神經(jīng)元網(wǎng)絡(luò)系統(tǒng)是一種高效的人工智能預(yù)測(cè)方法，其主要包括輸入層、隱藏層和輸出層，其中單隱藏層BP 神經(jīng)元學(xué)習(xí)效率不高，不適用于較大數(shù)據(jù)的非線性耦合數(shù)據(jù)的預(yù)測(cè)。

BP 神經(jīng)網(wǎng)絡(luò)可以用數(shù)學(xué)表達(dá)式表示為［46］

式中：xi為輸入層參數(shù)；wij為神經(jīng)元權(quán)重；aj為閾值；f為激活函數(shù)。常見(jiàn)的激活函數(shù)有l(wèi)ogsig、tansig和purelin 等，本文選用的雙曲線正切函數(shù)tansig 能夠較好地實(shí)現(xiàn)輸入層到隱藏層以及隱藏層之間信號(hào)和誤差傳遞關(guān)系。

如圖6 所示，本文所開(kāi)展的BP 神經(jīng)網(wǎng)絡(luò)模型的計(jì)算流程為利用本文第1 節(jié)中參數(shù)進(jìn)行批量仿真提取相應(yīng)應(yīng)變，隨后進(jìn)行移波處理，形成入射波、反射波和透射波。將所有參數(shù)進(jìn)行數(shù)據(jù)歸一化以提高計(jì)算速度和避免信號(hào)之間不同數(shù)量級(jí)帶來(lái)的計(jì)算誤差。此外，需要對(duì)BP 神經(jīng)網(wǎng)絡(luò)中每一個(gè)權(quán)重和閾值進(jìn)行參數(shù)初始化，將入射波作為輸入層，反射波和透射波作為輸出層，建立兩組神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)模型。本文使用多層神經(jīng)網(wǎng)絡(luò)技術(shù)，可以實(shí)現(xiàn)輸入層到輸出層的非線性映射，從而實(shí)現(xiàn)反射波和透射波的波形非線性預(yù)測(cè)。其中神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程見(jiàn)圖7，訓(xùn)練過(guò)程主要分為如下步驟：

圖6 BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.6 BP neural network structure

圖7 BP 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程Fig.7 BP neural network training process

（1）初始化BP 神經(jīng)網(wǎng)絡(luò)，為BP 神經(jīng)網(wǎng)絡(luò)提供輸入x、目標(biāo)Y和目標(biāo)精度θ，學(xué)習(xí)效率η。

（2）計(jì)算根據(jù)給定輸入BP 神經(jīng)網(wǎng)絡(luò)實(shí)際輸出值，依次計(jì)算每層隱含層以及單元輸出層，式中f為雙曲線正切tansig 函數(shù)，i表示輸入層個(gè)數(shù)，l表示隱藏層節(jié)點(diǎn)個(gè)數(shù)。

（3）計(jì)算目標(biāo)值和神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值之間誤差ek，其中Ok為計(jì)算值，wjk表示隱藏層最后一層和輸出層之間的權(quán)值，bk表示該連接情況下的閾值。

（4）誤差反向傳遞誤差，更新輸入層和隱藏層權(quán)值wij和隱藏層和輸出層wjk，以及閾值aj和bk。

（5）判斷反饋后數(shù)據(jù)正向傳遞預(yù)測(cè)值是否滿(mǎn)足誤差條件，若滿(mǎn)足則停止模型訓(xùn)練；若不滿(mǎn)足，則重復(fù)步驟（2），直至滿(mǎn)足預(yù)定誤差值。

2.2 BP 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)

本研究所建立的模型為BP 多層神經(jīng)網(wǎng)絡(luò)模型，輸入層為入射桿上入射波，輸出層為反射波和透射波，建立兩組人工網(wǎng)絡(luò)模型。訓(xùn)練樣本為ABAQUS 仿真模擬中提取的入射桿和透射桿應(yīng)變，其輸入峰值應(yīng)力從5 MPa 開(kāi)始，峰值應(yīng)力每增加5 MPa 進(jìn)行一組SHPB 仿真實(shí)驗(yàn)，共計(jì)20 組。將2 μs 作為一個(gè)時(shí)間節(jié)點(diǎn)，共計(jì)252 個(gè)節(jié)點(diǎn)，記錄相應(yīng)的入射桿和透射桿上的應(yīng)變值。圖8 給出了20 組具體預(yù)測(cè)結(jié)果，數(shù)據(jù)具有良好的規(guī)律性，可在一定程度上保證機(jī)器學(xué)習(xí)模型訓(xùn)練的合理性。

圖8 20 組入射桿和透射桿應(yīng)變-時(shí)間響應(yīng)Fig.8 Strain?time responses of 20 groups of incident bar and transmission bar

BP 神經(jīng)網(wǎng)絡(luò)是一種信號(hào)向前傳播，誤差向后傳播的多層神經(jīng)網(wǎng)絡(luò)，因此選擇合適的層數(shù)、節(jié)點(diǎn)數(shù)以及相應(yīng)的訓(xùn)練函數(shù)就顯得尤為重要。多層神經(jīng)網(wǎng)絡(luò)可以避免單層神經(jīng)網(wǎng)絡(luò)線性逼近的缺點(diǎn)，可以解決任何線性和非線性的輸入輸出問(wèn)題。每個(gè)隱藏層的節(jié)點(diǎn)數(shù)對(duì)于非線性耦合相當(dāng)敏感。如果節(jié)點(diǎn)數(shù)目過(guò)低會(huì)出現(xiàn)適應(yīng)性不足的現(xiàn)象，而如果節(jié)點(diǎn)的數(shù)目過(guò)多則會(huì)出現(xiàn)過(guò)渡適應(yīng)的特點(diǎn)，即出現(xiàn)每一個(gè)點(diǎn)都能夠很好地?cái)M合，但是最終擬合的結(jié)果在真實(shí)數(shù)據(jù)點(diǎn)之間往復(fù)振蕩，因此需要選取合適的節(jié)點(diǎn)數(shù)進(jìn)行機(jī)器學(xué)習(xí)。當(dāng)前對(duì)于節(jié)點(diǎn)數(shù)目并沒(méi)有準(zhǔn)確的方法來(lái)確定，一般可通過(guò)經(jīng)驗(yàn)公式（21）和（22）作為參考［47］。

式中：n1、n、m代表隱藏層節(jié)點(diǎn)數(shù)、輸入節(jié)點(diǎn)、輸出節(jié)點(diǎn)數(shù)，c為1～10 之間的一個(gè)常數(shù)。本文輸入層和輸出層均為252，由計(jì)算公式可知25～35 之間，本文選取的節(jié)點(diǎn)數(shù)為35。由于本文訓(xùn)練樣本為20組，每組輸入層參數(shù)高達(dá)252 個(gè)，因此選用量化共軛梯度算法可以避免出現(xiàn)標(biāo)準(zhǔn)梯度算法訓(xùn)練速度過(guò)慢的不足。反射波和透射波是兩種類(lèi)型波形，所以構(gòu)建入射波和反射波一組訓(xùn)練，入射波和透射波一組訓(xùn)練。學(xué)習(xí)樣本和驗(yàn)證樣本應(yīng)當(dāng)具有一定的比例，過(guò)小的學(xué)習(xí)樣本會(huì)使得建立的非線性映射關(guān)系丟失部分關(guān)鍵信息；過(guò)多的學(xué)習(xí)樣本會(huì)降低神經(jīng)網(wǎng)絡(luò)的泛化能力。因此本文對(duì)于20 組樣本進(jìn)行分類(lèi)，其中19 組作為學(xué)習(xí)樣本（Train），1 組作為驗(yàn)證樣本（Validation），驗(yàn)證和測(cè)試過(guò)程中的均方誤差（MSE）隨網(wǎng)絡(luò)訓(xùn)練迭代次數(shù)（Epoch）的變化規(guī)律。學(xué)習(xí)與驗(yàn)證樣本占全體的比例分別是95%和5%。兩組學(xué)習(xí)模型循環(huán)次數(shù)均為3 000 次，學(xué)習(xí)速率為0.05，設(shè)置的目標(biāo)MSE 分別為5×10-5和5×10-3。驗(yàn)證誤差分別連續(xù)300 次和500 次上升時(shí)終止訓(xùn)練以降低模型的過(guò)擬合的可能性。上述任一條件達(dá)到，模型學(xué)習(xí)終止。

在輸入波形時(shí)選用式（12）所定義的曲線，然而由波的彌散效應(yīng)，在圖8（a）中入射波在波頭出現(xiàn)了一定程度的振蕩現(xiàn)象，如果直接選用傳統(tǒng)撞擊桿方式建模，波頭振蕩將更加明顯，故本文選用在入射桿自由端施加應(yīng)力波法代替?zhèn)鹘y(tǒng)撞擊桿法建模，所得計(jì)算模型中反射波和透射波應(yīng)變隨著輸入應(yīng)力波增加而增加，更加符合實(shí)驗(yàn)規(guī)律。

圖9（a）是反射波最佳性能驗(yàn)證，在第2 536 次時(shí)，驗(yàn)證集的MSE 最小為0.007 923。圖9（b）是透射波最佳性能驗(yàn)證，在第605 次時(shí)，驗(yàn)證集的MSE最小為0.003 063 7。上述組圖表明在MSE 最小值以后網(wǎng)絡(luò)訓(xùn)練繼續(xù)進(jìn)行，在指定的訓(xùn)練數(shù)據(jù)上對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行迭代訓(xùn)練，直到滿(mǎn)足驗(yàn)證集誤差連續(xù)上升終止條件滿(mǎn)足時(shí)停止訓(xùn)練，這從本質(zhì)上能夠避免網(wǎng)格過(guò)度收斂問(wèn)題以提高網(wǎng)絡(luò)的泛化能力。

圖9 均方誤差?訓(xùn)練迭代次數(shù)變化圖Fig.9 MSE?epochs relational graph

本文選用訓(xùn)練樣本入射波峰值應(yīng)力極大值以外一組入射波（峰值應(yīng)力強(qiáng)度為122 MPa 入射波）作為人工神經(jīng)網(wǎng)絡(luò)輸入層，使用建立的預(yù)測(cè)模型預(yù)測(cè)反射波、透射波、應(yīng)變率和應(yīng)力?應(yīng)變曲線，其中應(yīng)變?yōu)橹苯宇A(yù)測(cè)結(jié)果，應(yīng)力?應(yīng)變曲線和應(yīng)變率曲線為式（1～4）的計(jì)算結(jié)果。所有數(shù)值模擬和人工神經(jīng)網(wǎng)絡(luò)運(yùn)算均在小型計(jì)算工作站（配置：AMD（R）Ryzen（TM）3700X CPU@3.6 GHz 中央處理器、32 GB 內(nèi)存）上完成。對(duì)于入射波為122 MPa峰值應(yīng)力情況下透射波和反射波仿真結(jié)果如下，ABAQUS 2017 商業(yè)有限元軟件用時(shí)942.5 s，人工神經(jīng)網(wǎng)絡(luò)方法用時(shí)15 s，計(jì)算時(shí)間縮短60 倍以上。如果針對(duì)批量混凝土材料的動(dòng)態(tài)力學(xué)性能仿真，時(shí)間優(yōu)勢(shì)將更加明顯。圖10 給出了根據(jù)122 MPa 峰值應(yīng)變下機(jī)器學(xué)習(xí)與仿真模擬的主要結(jié)果。圖10（b）所有機(jī)器學(xué)習(xí)預(yù)測(cè)小于有限元計(jì)算的可能原因?yàn)閷W(xué)習(xí)樣本中的所有反射波峰值均低于122 MPa峰值應(yīng)力強(qiáng)度下反射波峰值，導(dǎo)致機(jī)器學(xué)習(xí)預(yù)測(cè)出現(xiàn)保守預(yù)測(cè)，但距45°回歸線偏離不遠(yuǎn)，表明反射波預(yù)測(cè)結(jié)果良好。圖10（d）中透射波機(jī)器學(xué)習(xí)預(yù)測(cè)值和有限元計(jì)算值對(duì)應(yīng)點(diǎn)均勻分布在45°回歸線附近，表明預(yù)測(cè)結(jié)果良好。如圖10（a～d）所示，基于機(jī)器學(xué)習(xí)方法所提出的預(yù)測(cè)模型可以很好地復(fù)現(xiàn)ABAQUS 有限元計(jì)算的反射波和透射波結(jié)果的主要趨勢(shì)，且誤差較小。圖10（e）對(duì)比了兩種方法所對(duì)應(yīng)的材料應(yīng)力?應(yīng)變曲線，在上升段和極值段的誤差控制在5%以?xún)?nèi)。類(lèi)似地，圖10（f）給出了應(yīng)變率的對(duì)比圖，在整個(gè)分析過(guò)程中，兩種方法預(yù)測(cè)得到的應(yīng)變率基本一致。因此可以認(rèn)為，機(jī)器學(xué)習(xí)預(yù)測(cè)模型可以相對(duì)可靠地預(yù)測(cè)材料動(dòng)態(tài)力學(xué)性能。

圖10 122 MPa 峰值應(yīng)變率下機(jī)器學(xué)習(xí)模型與數(shù)值仿真的結(jié)果對(duì)比Fig.10 Result comparison between machine learning model and numerical simulations for typical cases with the peak stress of 122 MPa

BP 算法使用反饋方法反復(fù)修正權(quán)重和閾值，故本文選取的35 個(gè)節(jié)點(diǎn)需要適應(yīng)全體的252 點(diǎn)。在應(yīng)變極小或應(yīng)力極小的點(diǎn)，較小的誤差都會(huì)形成極大的相對(duì)誤差率，因此本文選用反射波和透射波峰值應(yīng)變20%以上的點(diǎn)的平均相對(duì)誤差率E和Pearson 相關(guān)系數(shù)R來(lái)衡量?jī)山M數(shù)據(jù)的擬合程度，其計(jì)算公式分別為［46］

式中：Xi為ABAQUS 仿真計(jì)算值；Yi為人工神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值；N為滿(mǎn)足條件的預(yù)測(cè)點(diǎn)個(gè)數(shù)，對(duì)于相關(guān)系數(shù)的大小所表示的意義，目前在統(tǒng)計(jì)學(xué)研究中所得結(jié)論尚不一致，但通常認(rèn)為0.8～1.0 的取值代表兩種數(shù)據(jù)之間具有高度的相關(guān)性。從表3 可以看出：122 MPa 應(yīng)力入射波加載條件下，機(jī)器學(xué)習(xí)預(yù)測(cè)反射波與仿真模擬反射波曲線間相關(guān)系數(shù)高達(dá)0.987 0，平均相對(duì)誤差5.050%，預(yù)測(cè)透射波與仿真反射波相關(guān)系數(shù)0.967 4，平均相對(duì)誤差12.37%，預(yù)測(cè)應(yīng)變率與仿真預(yù)測(cè)應(yīng)變率相關(guān)系數(shù)0.979 5，平均相對(duì)誤差9.15%，上述3 類(lèi)預(yù)測(cè)均呈現(xiàn)強(qiáng)相關(guān)性。在動(dòng)態(tài)相應(yīng)峰值應(yīng)力之前，機(jī)器預(yù)測(cè)應(yīng)力?應(yīng)變和仿真應(yīng)力?應(yīng)變誤差在5%，表明本文所提出的基于人工神經(jīng)網(wǎng)絡(luò)SHPB 實(shí)驗(yàn)預(yù)測(cè)模型有效可靠且具有較高的精度，可以代替量大且耗時(shí)的有限元仿真建模、分析及后處理流程，實(shí)現(xiàn)更大應(yīng)變率范圍下材料應(yīng)力?應(yīng)變曲線預(yù)測(cè)。

表3 122 MPa峰值應(yīng)力下預(yù)測(cè)應(yīng)變、應(yīng)變率結(jié)果衡量系數(shù)Table 3 Correlation coefficient of strain and strain rate predictions for the input case with peak stress of 122 MPa

3 結(jié)論

隨著混凝土類(lèi)材料應(yīng)用場(chǎng)景的日益廣泛，研究各種應(yīng)變率下材料動(dòng)態(tài)本構(gòu)顯得尤為重要，本文通過(guò)機(jī)器學(xué)習(xí)和有限元仿真相結(jié)合的方法，得出如下研究成果：

（1）本文建立混凝土SHPB 實(shí)驗(yàn)?zāi)Ｐ陀邢拊抡婺Ｐ瞳@取訓(xùn)練樣本，所得的人工神經(jīng)網(wǎng)絡(luò)動(dòng)態(tài)混凝土力學(xué)性能預(yù)測(cè)模型有效地反映混凝土類(lèi)材料動(dòng)態(tài)入射波和反射波、入射波和透射波以及入射波和應(yīng)力應(yīng)變之間的非線性映射關(guān)系。

（2）所得的SHPB 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值與有限元仿真計(jì)算值相關(guān)性在0.95 以上，表明本文所建立的非線性映射關(guān)系具有很強(qiáng)的相關(guān)性，具有良好的精度。整個(gè)過(guò)程中各項(xiàng)參數(shù)設(shè)置合理，ABAQUS 計(jì)算值和文獻(xiàn)DIF 值具有良好的吻合度，證明將BP人工神經(jīng)網(wǎng)絡(luò)應(yīng)用于混凝土類(lèi)材料應(yīng)力應(yīng)變預(yù)測(cè)是可行的。

（3）該方法能預(yù)測(cè)給定訓(xùn)練樣本以外更大應(yīng)變率范圍下材料應(yīng)力?應(yīng)變曲線，實(shí)現(xiàn)高應(yīng)變率下混凝土類(lèi)材料應(yīng)力?應(yīng)變曲線形式本構(gòu)行為的高效準(zhǔn)確預(yù)測(cè)。與傳統(tǒng)的實(shí)驗(yàn)相比減少人工重復(fù)勞動(dòng)，節(jié)省時(shí)間物料成本，同時(shí)可代替量大且耗時(shí)的有限元仿真建模、分析及后處理流程，具有明顯優(yōu)勢(shì)。

本文所提的方法雖然能夠預(yù)測(cè)混凝土類(lèi)材料動(dòng)態(tài)力學(xué)性能，但是BP 人工神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)本身的權(quán)重和閾值調(diào)整是隨機(jī)無(wú)方向的，因此可以進(jìn)一步結(jié)合優(yōu)化算法，確定最優(yōu)權(quán)重和閾值，從而使所提出機(jī)器學(xué)習(xí)模型具有更好的普適性。