混凝土材料動態本構參數的分階段計算反求技術＊

陳 睿，劉 杰，韓 旭，畢仁貴

（1.湖南大學汽車車身先進設計制造國家重點實驗室，湖南 長沙410082；2.湖南大學機械與運載工程學院，湖南 長沙410082）

混凝土材料廣泛應用于防撞擊和防爆結構、機場跑道和核電站等民用和軍用設施，因此研究混凝土材料的動態力學性能具有十分重要的理論意義和現實指導作用。混凝土是一種復合材料，主要由水泥砂漿、骨料及黏結層組成。由于水泥砂漿和骨料的物理、力學性能差異明顯，骨料的形狀不規則和分布的不確定，以及含有大量不規則的初始裂紋和氣泡等因素，導致了混凝土材料的動態力學性能十分復雜。近幾十年來，學者們對混凝土材料的動態力學性能進行了廣泛的研究［1－3］，構建了許多關于混凝土材料的動態本構模型，常用的有 HJC模型［4］、RHT 模型［5］及 TCK 模型［6］等。

材料本構模型都含有很多參數，這些參數中大多數可以通過物理和力學實驗進行確定，但有些參數不能通過實驗直接獲取，常常需要大量不同的實驗數據，通過擬合來確定。然而對混凝土材料而言，由于混凝土試樣制備的困難和實驗設備的限制，要獲得大量有效的實驗數據是相當耗時費力的。因此，尋求其他快捷途徑來確定這些難以確定的參數是十分必要的。近些年來，基于材料的動態響應來反求材料特性取得了一定的進展［7－11］。這類參數反求方法中，材料的動力響應與材料參數之間的復雜關系一般由正問題的計算模型來表示。材料參數的反求是通過反復調整正問題計算模型中的材料參數，使得計算模擬響應與測量響應之間的誤差最小，最終獲得與測量響應吻合最好的材料參數。

通常，這類參數反求方法都是基于所有實驗測量響應數據反求出所有與測量響應具有一定敏感性的參數。這樣雖然能夠確定這些參數的取值，但由于未知參數對不同類型實驗測量響應的敏感性差異以及不同類型實驗測量響應的量級不同，會導致在一定程度上缺乏對實驗測量響應數據的綜合利用，并影響參數反求的精度。為了充分利用實驗測量響應信息，本文中提出了一種分階段參數計算反求技術以獲取混凝土HJC本構模型中的關鍵參數。結合不同應變率下的混凝土SHPB實驗，基于參數敏感性分析及分類，將透射波和反射波響應作為輸入信息，分階段反求出混凝土HJC本構參數，并將確定參數后的計算結果與實驗結果進行對比分析。

1 HJC動態本構模型及其參數

混凝土HJC本構模型［4］是一種綜合考慮大應變、高應變率及高壓效應的混凝土計算本構模型。該模型考慮了材料損傷、應變率效應以及靜水壓力對屈服應力的影響，主要包括3方面：等效強度模型、損傷模型和狀態方程。其中等效強度模型為：

式中：σ＊＝σ/fc，P＊＝P/fc和分別表示歸一化的等效強度、歸一化壓力和量綱一應變率。其中σ、P和分別為實際的等效強度，壓力和應變率，fc為準靜態單軸抗壓強度，為參考應變率。A、B、N、C是材料的強度模型參數，分別表示材料內聚強度、壓力硬化系數、壓力硬化指數、應變率硬化系數。D 為損傷因子（0≤D≤1.0），并引入一個最大量綱一強度Smax，σ＊≤Smax。

混凝土HJC動態本構模型［4］共含有20個參數，可以通過不同的途徑來確定。其中，可以直接通過簡單的物理或力學實驗獲取4個參數，即初始密度ρ0、彈性模量E、拉伸強度T和fc；可以間接通過實驗或經典的計算公式［4］確定6個參數，即2個損傷參數（D2，εf，min）和4個壓力參數（彈性體積模量K、彈性極限靜水壓力Pcrush及對應的體應變μcrush、壓實體應變μlock）；而其他10個參數（5個強度參數A、B、N、C、Smax，1個損傷參數D1和4個壓力參數（壓實靜水壓力Plock、K1、K2、K3））需要通過大量有效的實驗數據進行擬合來確定。然而要獲取大量有效的實驗數據成本較高，且費時耗力，想要快速獲取這些參數是較為困難的。因此本文中將采用分階段計算反求技術來快速獲取這些較難確定的參數。

2 混凝土本構參數反求

2.1 分階段計算反求技術

在利用計算反求技術進行參數確定過程中，通常都是基于所有實驗測量響應數據來同時反求所有的未知參數。這往往會由于未知參數對實驗測量響應的敏感程度不同，測量響應數據與未知參數類型不匹配，不同類型實驗測量響應數據量級水平不同等問題，使得在一定程度上缺乏對測量響應數據的綜合利用，并導致反求結果不夠理想。為此，本文中提出了一種分階段參數反求技術，該方法基于多種不同的實驗，將參數反求問題進行細化，對未知參數進行敏感性分析及分類，并在此基礎上分階段反求出各類參數，且每一次參數反求均利用前一次反求結果作為基礎，該參數反求方法的流程如圖1所示。

圖1 分階段參數反求流程圖Fig.1 Flow chart of the multi－stage inverse technique for parameters identification

分階段參數反求方法中最明顯的特點在于基于不同實驗建立的正問題，對所有未知參數進行敏感性分析，根據各參數的敏感性高低程度，將這些未知參數劃分為幾類（每類對應1組實驗）。在進行分階段參數反求時，首先對敏感性程度最高的一類參數，結合相應的實驗，利用反求方法確定這類參數。確定之后采用同樣的方法進行其他類參數反求，在每一次參數反求時都將上一次反求的結果作為當前參數反求的已知條件。這樣分階段地對各類參數進行反求，最后獲得所有參數的取值。在本文中每階段參數反求過程中，采用遺傳算法［12］作為反問題的求解器，該算法的最大優勢在于其具有較強的全局最優搜索能力和較快的全局收斂速度，能夠克服局部最優，并在全局范圍內尋找到最優解。

本文中將利用該分階段參數反求方法，結合不同應變率下的混凝土SHPB實驗，對混凝土HJC本構中關鍵參數進行反求。該SHPB實驗數據來源于文獻［3］對混凝土試樣進行3種不同應變率（case 1：17 s－1；case 2：25 s－1；case 3：200 s－1）加載的實驗。文獻［3］給出了較為詳細的SHPB實驗裝置的幾何尺寸和混凝土試件的基本材料參數，以及測試得到的響應曲線。

2.2 正問題的建立

基于文獻［3］中SHPB實驗裝置，建立相對應的有限元模型。由于撞擊桿、輸入桿、輸出桿和試件的幾何尺寸都是軸對稱的，同時為了提高計算效率，采用PLANE162平面軸對稱單元建立二維軸對稱有限元模型，如圖2所示。撞擊桿、壓桿和試件在徑向分別劃分15、15、30個單元，軸向分別劃分740、2 300、60個單元。在數值計算中，采用LS－DYNA軟件作為正問題的求解器，且兩兩接觸采用單面自動接觸算法。

圖2 SHPB的二維軸對稱有限元模型Fig.2 2D axisymmetric finite element model for the SHPB experiment

撞擊桿、輸入桿和輸出桿均為相同的彈性材料，因此都采用線彈性材料模型，材料參數為：密度7.83 g/cm3、彈性模量200 GPa、泊松比0.3。混凝土試件采用 HJC動態本構模型，其中ρ0、E、T 和fc等基本參數源自文獻［3］，材料的壓實密度ρgrain和Plock源自文獻［13］。鑒于第1節中所述，結合文獻［4］中的計算公式可獲取K、Pcrush、μcrush、μplock和μlock。由于在SHPB實驗的應變率范圍內混凝土的壓力僅能達到過渡段，使得狀態方程中壓實段的參數對數值計算結果影響較小［14］，因此該段方程的參數K1、K2和K3取原始文獻［4］的數值。另外，數值研究發現最大量綱一強度參數和損傷參數的變化對該問題的計算結果影響較小，因此本文中Smax、D1和D2也取原始文獻［4］中的數值。因此，除了強度參數A、B、C、N 外，混凝土HJC模型的其他參數如表1所示。

表1 混凝土HJC本構模型的基本參數Table 1 The material parameters for HJC concrete constitutive model

2.3 參數敏感性分析及分類

利用計算反求技術進行參數反求時，要求反求的參數必須對輸出的響應具有較強的敏感性。分別基于低應變率和高應變率的SHPB實驗，對強度參數A、B、C、N進行敏感性分析。基于低應變率工況（case 1）對這4個參數進行敏感性分析，結果如圖3所示。從圖3中可以發現，參數A、B、N 取不同值時，計算響應結果（反射波和透射波）變化很大，而參數C取值的變化對計算響應結果幾乎沒有影響。這說明在低應變率工況下，參數A、B、N 具有較高的敏感性，而參數C的敏感性較低。同樣基于高應變率工況（case 3）下對4個參數進行敏感性分析，發現計算響應結果隨著這4個參數取值的變化而變化，其中參數C的敏感性分析結果如圖4所示。這說明在高應變率工況下，4個參數都具有較高的敏感性，尤其參數C的敏感性很高。參數C表示應變率參數，且從公式（1）中可以看出應變率對強度的影響表現為應變率的對數，這就可以發現當處于低應變率時，（1＋C lnε′＊）的值幾乎在1.0左右；而當在高應變率ε·≥102s－1時，由于量級的影響，使得參數C的變化對（1＋C lnε′＊）的影響很大，這進一步說明在高應變率下，參數C具有較強的敏感性。

圖3 基于工況1對參數A、B、C、N 敏感性分析Fig.3 Sensitivity analysis for parameters A，B，C，N based on case 1

圖4 基于工況3對參數C敏感性分析Fig.4 Sensitivity analysis for parameter C based on case 3

基于上述敏感性分析及各參數的敏感性程度，將4個強度參數分為2類，即第1類包含3個參數A、B、N，這一類參數將基于低應變率工況的實驗進行參數反求。第2類為參數C，這類參數將基于高應變率工況實驗進行參數確定，在反求過程中，將第1類參數反求的結果作為已知條件。

2.4 結果及驗證

結合實驗測量的響應及參數敏感性分析，可以初步確定參數的取值范圍，即A∈［0.2，1.0］，B∈［0.4，1.5］，N∈［0.5，1.2］，C∈［0.05，0.2］，且在反求時這4個參數的初始值都設為取值范圍的平均值。在敏感變化的透射波和反射波響應790～1 100μs區域中采樣24個時刻的離散點以構造優化目標函數：

基于低應變率case 1實驗，利用上述的反求方法確定第1類參數A、B和N，計算反求結果如表2所示。

表2 分階段參數反求的結果Table 2 The inversed results by multi－stage inverse technique

為了驗證反求參數的適用性，分別將反求得到的參數代入正問題，在2種低應變率工況（case 1和case 2）下進行數值計算，并將計算結果與實驗數據進行對比。其中case 2實驗沒有參與參數反求，而作為反求參數適用性的驗證。圖5～6分別給出了2種低應變率下的計算結果，從圖中可以看出，2種應變率下的計算結果與實驗曲線較為吻合。

基于上述第1類參數反求的結果，結合較高應變率下的實驗數據，再利用反求技術確定第2類參數C的值，反求結果如表3所示。將得到的參數C值代入正問題計算，其結果如圖7所示。從圖中可以看出，數值計算的結果與實驗測量數據較為吻合。這說明基于本文分階段計算反求技術進行參數確定是可行的，所得到的結果較為準確。

圖5 工況1的計算結果與實驗結果比較Fig.5 The computational results for case 1 compared with experimental data

圖6 工況2的計算結果與實驗結果比較Fig.6 The computational results for case 2 compared with experimental data

圖7 工況3的計算結果與實驗結果比較Fig.7 The computational results for case 3 compared with experimental data

3 結 論

提出了一種分階段的計算反求技術以確定混凝土本構模型中的關鍵參數。在該方法中，基于2種不同應變率下的SHPB實驗，對未知參數進行敏感性分析及分類，在此基礎上分階段反求各類參數，且每次參數反求時均利用上一次反求結果作為基礎，從而分階段確定了混凝土本構參數。本文中每階段的參數反求都是根據透射波和反射波響應反求材料本構參數，這樣避免了應力波效應和應變率效應解耦問題，不必分離結構響應和材料響應。參數反求結果表明，該分階段計算反求方法充分利用不同的實驗響應，能較為快速、穩定地確定一些傳統方法較難以確定的參數。該計算反求方法充分結合實驗與數值計算，通過實驗響應確定參數，減少實驗次數，提高效率，在工程實際中具有重要意義。

［1］Ross C A，Thompson P Y，Tedesco J W.Split－Hopkinson pressure－bar tests on concrete and mortar in tension and compression［J］.ACI Materials Journal，1989，86（5）：475－481.

［2］Bischoff P H，Perry S H.Compressive behaviour of concrete at high strain rates［J］.Materials and Structures，1991，24（6）：425－450.

［3］Tedesco J W，Hughes M L，Ross C A.Numerical simulation of high strain rate concrete compression tests［J］.Computers ＆ Structures，1994，51（1）：65－77.

［4］Holmquist T J，Johnson G R，Cook W H.A computational constitutive model for concrete subjected to large strains，high strain rates，and high pressures［C］∥Jackson N.Proceedings of the 14th International Symposium on Ballistics.USA：American Defense Preparedness Association，1993：591－600.

［5］Riedel W，Thoma K，Hiermaier S，et al.Penetration of reinforced concrete by BETA－B－500 numerical analysis using a new macroscopic concrete model for hydrocodes［C］∥Proceedings of the 9th International Symposium on Interaction of the Effects of Munitions with Structures.Berlin，Germany，1999：315－322.

［6］Taylor L M，Chen E P，Kuszmaul J S.Microcrack－induced damage accumulation in brittle rock under dynamic loading［J］.Computer Methods in Applied Mechanics and Engineering，1986，55（3）：301－320.

［7］Han X，Liu G R.Computational inverse technique for material characterization of functionally graded materials［J］.American Institute of Aeronautics and Astronautics，2003，41（2）：288－295.

［8］Cooreman S，Lecompte D，Sol H，et al.Elasto－plastic material parameter identification by inverse methods：Calculation of the sensitivity matrix［J］.International Journal of Solids and Structures，2007，44（13）：4329－4341.

［9］Sedighi M，Khandaei M，Shokrollahi K H.An approach in parametric identification of high strain rate constitutive model using Hopkinson pressure bar test results［J］.Materials Science and Engineering：A，2010，527（15）：3521－3528.

［10］伍乾坤，韓旭，胡德安.一種陶瓷脆性材料動態本構參數的計算反求方法［J］.固體力學學報，2009，30（3）：280－285.Wu Qian－kun，Han Xu，Hu De－an.An inverse technique for identification of dynamic constitutive parameters of ceramic brittle material［J］.Acta Mechanica Solida Sinica，2009，30（3）：280－285.

［11］Chen R，Han X，Liu J，et al.A computational inverse technique to determine the dynamic constitutive model parameters of concrete［J］.Computers Materials ＆ Continua，2011，25（2）：135－157.

［12］Liu G R，Han X.Computational inverse techniques in nondestructive evaluation［M］.Boca Raton：CRC Press，2003.

［13］Hao Y F，Hao H，Li Z X.Numerical analysis of lateral inertial confinement effects on impact test of concrete compressive material properties［J］.International Journal of Protective Structures，2010，1（1）：145－167.

［14］巫緒濤，李耀，李和平.混凝土 HJC本構模型參數的研究［J］.應用力學學報，2010，27（2）：340－344.Wu Xu－tao，Li Yao，Li He－ping.Research on the material constants of the HJC dynamic constitutive model for concrete［J］.Chinese Journal of Applied Mechanics，2010，27（2）：340－344.