基于Matlab的花崗巖RHT模型參數確定研究

殷 達，趙 寧，劉雄峰，彭 俊，榮 冠

（1.中交廣州航道局有限公司，廣州 510290；2.武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

RHT 模型是由Riedel 等通過改進HJC 模型提出的混凝土模型［1］，RHT 模型在混凝土、陶瓷等脆性材料的沖擊損傷方面被廣泛應用［2，3］，后來又被引入描述巖石在動載作用下的力學特性［4］，如被應用與巖石的數值模擬研究（如霍普金森壓桿試驗的模擬或巖石受彈丸沖擊的模擬）中。RHT 模型考慮了偏應力張量第三不變量對破壞面的影響，引入了與壓力相關的彈性極限面、失效面和殘余強度面，分別描述了材料的初始屈服強度、失效強度及殘余強度的變化規律，更為全面地反映了巖石和混凝土等材料在不同應力狀態下的動態力學行為，在參考文獻中已經有關于RHT 模型的詳細描述［5］，該模型中需要確定的參數有30 余個。根據文獻研究成果，在RHT 模型參數當中，部分參數可通過室內靜力學和波速測試等試驗獲取得到，但是13個參數例如剪壓強度比f*

s、孔隙度指數N、拉壓強度比f*t、孔隙壓實時壓力pcomp等極難獲取，工程人員常常引用文獻中混凝土等類似材料的相關參數取值［6，7］，但該方式往往導致特定材料如花崗巖的RHT模型的估計偏差較大，給數值仿真計算帶來了不利影響，可能致使仿真結果不合理［7］。

在采用沖擊錘破巖或爆破破巖的工程中，為優化設計和施工方案，往往需要準確掌握巖石在動載作用下的力學特性。由于地質成因及環境不同，不同場地巖石材料的RHT模型參數不盡相同。為了準確模擬某一特定項目花崗巖在動載作用下的力學特性，對于特定工程花崗巖的RHT 模型參數，需要結合室內SHPB 試驗實施正交模擬從而確定其準確值。該過程需要多次進行數值模擬，從每次模擬得到的結果文件中提取特征值繪制曲線，進而優化確定RHT模型參數值。如果通過手動改寫輸入參數和提取結果文件［4，6］，則該過程較為繁瑣，效率極低。此外，花崗巖在動載作用下具有應變率效應［10］，有時在一個項目為了確定一類花崗巖的RHT模型參數，可能需要結合不同應變率情況下的應力-應變曲線進行模擬，進一步加大了工作量，給工程人員確定RHT 模型參數帶來了極大不便。Matlab 作為數據分析和建模等方面強大的商業數學軟件，在繪制函數和數據、實現算法、連接其他計算程序如LS-DYNA 等方面具有較大優勢［11］。本文運用Matlab 軟件對LS-DYNA 軟件的輸入文件、執行文件和結果文件進行修改和處理，以求提高該優化分析過程的效率。首先介紹了深中（深圳-中山）通道沉管隧道工程現場花崗巖RHT 模型參數的確定方法，而后介紹了基于Matlab 軟件和LS-DYNA 軟件的聯合數值模擬，利用正交優化方法在前文介紹內容的基礎上確定RHT模型參數。

1 RHT模型參數確定方法

強度模型、損傷模型和狀態方程組成的RHT 模型，在模擬巖石類材料沖擊分析中得到廣泛應用，運用過程中總共涉及34個參數。

1.1 結合靜力學和聲波測試試驗確定部分參數

首先，基于靜力學和聲波測試試驗確定深中通道沉管隧道工程現場花崗巖的部分RHT 模型參數。而后對于巖樣進行三軸、單軸、波速測定等試驗，得到了花崗巖的基本物理參數如表1所示。

表1 某項目微風化花崗巖基本物理參數Tab.1 Basic physical parameters of slightly weathered granite in a project

以花崗巖基本物理參數為基礎，根據文獻中結合靜力學和聲波測試試驗確定部分參數的常用公式［4，6，7］，可以計算出參數A1、A2、A3、B0、B1、T1、T2、Pel、G、βc、βt、fc、ρ0和α0。一般來說和D2是跟隨模型的定值，且B=0.010 5 和=0.7。故而得到了RHT模型中的21個參數取值如表2中所示。

表2 由靜力學試驗和理論分析確定的某項目花崗巖RHT模型部分參數取值Tab.2 Some parameters of granite RHT model determined by static test and theoretical analysis

表3 參考文獻研究確定的花崗巖RHT模型部分參數近似值Tab.3 Approximate values of some parameters of granite RHT model determined by references

1.2 結合SHPB試驗進行正交試驗確定剩余參數

以SHPB 沖擊試驗為基礎，以LS-DYNA 數值模擬SHPB 沖擊試驗為手段，利用正交試驗優化確定這13 個參數，下面簡要介紹該參數優化的過程。首先，需要基于室內SHPB 試驗獲取巖石應力-應變曲線，如圖1中藍色曲線所示，以供后續對比分析。

而后，在LS-DYNA 有限元軟件中對于SHPB 試驗進行模擬。其中，巖石試樣采用RHT 本構模型（模型參數為表2、3 中初步擬定的值）來描述其在沖擊作用下的力學行為。試樣和壓桿之間的接觸使用的是侵蝕面面接觸算法。入射桿和透射桿的材料模擬都是使用的線彈性模型。由此模擬實驗可分析得到巖石的應力-應變模擬曲線，如圖1中紅色曲線所示。若模擬結果和實驗曲線重合度較高，則確定的參數合理。

圖1 SHPB試驗應力-應變曲線和LS-DYNA軟件模擬的應力-應變曲線示意圖Fig.1 The stress-strain curve of SHPB test and the stress-strain curve simulated by LS-DYNA software

顯然，由于上述待定參數有13個，試驗量較大，可采用正交試驗進行參數確定。正交試驗的試驗因素即為13 個未確定的參數。試驗因素選定后，要對每個因素的水平進行確定，參考文獻研究，進行13 因素三水平試驗，即選擇標準正交表L27（313）進行模擬正交試驗。在完成數值模擬后，即可對正交試驗模擬結果與SHPB 沖擊巖石試驗曲線進行對比分析，計算各個試驗方案對應的試驗指標（即用“距離”表示的模擬曲線與SHPB 沖擊試驗曲線的接近程度），由此便能夠得到試驗因素的最優水平和參數之間的最優組合，從而得到最合理的參數取值。

該方法步驟雖然清晰，但需要重復進行多次數值模擬，從每次模擬得到的結果文件中提取特征值繪制曲線，進而優化確定RHT 模型參數值。如通過手動改寫輸入參數和提取結果文件，則該過程較為繁瑣，效率極低。此外，花崗巖在動載作用下具有應變率效應，導致需要結合不同應變率情況下的多條應力-應變曲線進行模擬，進一步加大了工作量，給工程人員確定RHT 模型參數帶來了極大不便。下一節介紹本文提出的基于Matlab 軟件和LS-DYNA 軟件的聯合數值模擬方法，以求提高該優化分析過程的效率。

2 基于Matlab 軟件和LS-DYNA 軟件的聯合數值模擬方法

本文提出了一種Matlab 與LS-DYNA 的聯合數值計算的方法，用以實施SHPB試驗數值模擬的正交試驗，從而確定RHT模型的參數，Matlab 作為數據分析和建模等方面強大的商業數學軟件，在繪制函數和數據、實現算法、連接其他計算程序如LSDYNA 等方面具有較大優勢。LS-DYNA 是適合求解非線性動力學問題，可用于模擬巖石在動態荷載作用下的損傷和破裂特征［12，13］。RHT 本構模型可用于描述巖石的動力學特性，它已被嵌入LS-DYNA 軟件中［14］。需要注意的是每一個實際工程項目的RHT 模型的參數需要進行標定。該標定過程往往需要采用正交試驗，也就是反復多次進行不同RHT 模型參數條件下的SHPB 試驗數值模擬工作。但是這個想法的實現十分的繁瑣，人工實現難度極高。本文所提方法可大大提高標定參數的效率，流程圖如圖2。

圖2 基于Matlab和LS-DYNA的聯合數值模擬流程圖Fig.2 Flow chart of joint numerical simulation based on MATLAB and LS-DYNA

首先，運用LS-DYNA 建立花崗巖SHPB 試驗的數值模型，包含幾何模型建立、網格劃分等，形成模型文件（k 文件）。k 文件中RHT模型參數的值可任意指定。而后，根據正交試驗的一般方法，制定正交試驗方案，確定每次試驗的RHT 模型參數。Matlab 可打開并修改k 文件，因此可將每一套參數代入Matlab中用于改寫k文件。下為基于Matlab編寫的核心代碼：

在完成所有k 文件的修改后，在Matlab 中通過執行LS-DY‐NA 的執行文件從而實現批量提交計算。本部分Matlab 核心代碼如下：

提交計算后，LS-DYNA 后臺自動進行計算，計算完成后，將產生每次試驗的結果文件數據庫。然后，采用Matlab 讀取結果文件，用于優化分析出RHT 模型參數。下為基于Matlab 編寫的核心代碼：

通過采用上述Matlab 與LS-DYNA 軟件結合的數值計算方法，用于下一節正交試驗中的數據處理、仿真計算以及結果的讀取和比對，使得繁雜的基于LS-DYNA 的模擬工作在Matlab中全部自動實現，提高RHT模型參數優化效率。

依據上述確定參數的框架，下面將介紹深中通道隧道項目花崗巖的RHT模型參數的確定過程。

3 基于正交試驗確定RHT模型參數

3.1 SHPB室內試驗結果

試驗中，撞擊桿以某一速度沖擊入射桿桿端，在入射桿桿端產生壓應力波并向入射桿另一端傳播。試驗選取了兩種撞擊桿的發射速度，對應的兩種入射波的應力持續時間均為約450 μs，應力峰值分別為75 MPa 和100 MPa，分別如圖3中藍色曲線所示。圖3中紅色和綠色曲線代表經過室內試驗獲得的反射波和透射波曲線。

圖3 項目室內試驗測得的“三波”曲線Fig.3 “Three wave”curve measured in indoor test of the project

運用“三波算法”公式中，根據圖3中的波形數據，求得試樣的平均應力、應變和應變率，試樣應力-應變曲線如圖4所示。

圖4 項目室內試驗測得的應力-應變曲線Fig.4 The stress-strain curve measured by the indoor test of the project

3.2 RHT模型參數的確定

表4 花崗巖RHT模型參數正交試驗的因素水平Tab.4 Factor level of orthogonal test for parameters of granite RHT model

最后是正交試驗的模擬。對每種應變率情況，分別各進行一組（共27 次模擬）模擬，得到共三組應力-應變曲線。在正交試驗完成之后，得到共27 次模擬的應力波信號，限于篇幅未予展示。依據每一次所得的“三波”推算得到應力-應變曲線見圖5。

圖5 進行正交試驗（共27次模擬）所得的應力-應變關系Fig.5 The stress-strain relationship was obtained by orthogonal test（27 simulations）

提取27 組模擬曲線上的控制點和實驗室SHPB 沖擊試驗曲線上的相關點，利用距離公式處理分析得到試驗指標值見表5。然后計算K1、K2和K3，確定試驗因素的優水平和優組合，然后對照表4可以得到正交試驗優化確定的花崗巖RHT 模型部分參數，具體取值見表6。

表5 進行正交試驗（共27次模擬）所得的試驗指標值Tab.5 The test index values were obtained by orthogonal test（27 times of simulation）

表6 正交試驗優化確定的花崗巖RHT模型部分參數值Tab.6 Some parameters of granite RHT model optimized by orthogonal test

本節確定RHT 模型參數的全部過程可通過采用Matlab 與LS-DYNA 軟件結合的數值計算方法實現，尤其是使得數十次的基于LS-DYNA 的模擬工作在Matlab 中全部自動實現，為實際項目中確定花崗巖的RHT模型參數提供了一種有效的手段。且所提方法可以被集成為一整套程序，在操作上方便快捷，極大降低了RHT模型參數確定的門檻。

4 RHT模型參數的有效性驗證

4.1 結合項目室內試驗進行驗證

利用前文中得到的兩組參數（具體可見表3、6）分別進行SHPB 沖擊花崗巖數值模擬試驗。將模擬曲線與實驗室SHPB沖擊試驗曲線進行對比，驗證所得參數的合理性。

圖6所示為未經過優化處理的參數得到的數值模擬結果與室內試驗所得的對比。其中紅色曲線對應的花崗巖RHT 模型近似參數中和pcomp是直接引用文獻中混凝土的相關參數取值，將試驗參數與優化前的模擬參數所得曲線進行比較，可以看出在關鍵點處數據相差較大，因而不加優化的引用他人確定的參數取值是不合理的。

圖6 優化前RHT模型近似參數模擬曲線與實驗室SHPB沖擊曲線的對比結果Fig.6 Comparison between the simulation curve before optimization and SHPB impact curve

然后，采用表6中的參數取值進行SHPB 沖擊花崗巖數值模擬試驗得到模擬曲線，而后將其與實驗室所得曲線對比，二者對比關系如圖7所示。由圖可知，試驗與模擬所得入射波曲線大致是一樣的。對于應力強化段，數值模擬得到的曲線與試驗得到的曲線重合度高，但是后半段中二者重合度較低。這是因為，花崗巖組成成分復雜且具有裂紋等微缺陷，表現出非均質性，這一點在數值模擬中無法被考慮進去，因而數值模擬結果不佳。由對比結果可知，花崗巖RHT模型參數經過正交試驗優化后得到的模擬曲線可以與室內試驗獲得的曲線可以很好地貼合，說明此時所采用的RHT模型參數適用于本項目微風化花崗巖。

圖7 優化后的花崗巖RHT模型參數模擬曲線與實驗室SHPB沖擊曲線對比結果Fig.7 Comparison between the optimized simulation curve and SHPB impact curve

4.2 結合不同應力峰值和荷載時長的影響進行驗證

采用常規動載SHPB 試驗數值模型，本節探索不同應力峰值和荷載時長影響的數值模型只需要在入射桿端加載不同峰值或不同時長的應力波。在試樣受壓后分析其應力應變以及損傷破壞形態，得出相關結果。

（1）不同應力峰值的影響。為了研究不同應力峰值情況下巖石的破壞規律，參考試驗中的加載波，數值模擬中的應力持續時間為250 μs，應力峰值分別為200、225、250、275、300、325、350、400、500、600 MPa。

圖8是通過模擬所得入射波、反射波和透射波數據計算得到的不同應力峰值下的應力-應變曲線。從中可以明顯看出隨著應力峰值的增加，試樣極限強度提高，說明巖石類材料存在明顯的應力峰值依賴性，隨著應力峰值的增大，巖石的動載抗壓強度增大趨勢更加明顯。

圖8 不同荷載峰值情況下的應力-應變曲線Fig.8 Stress-strain curves under different peak loads

從圖9應變率曲線中可以明顯看出隨著應力峰值的增加，試樣應變率提高。因此，圖9中反映的應力-應變曲線隨不同應力峰值變化的規律也可以認為是在不同應變率情況下的變化規律。

圖9 不同荷載峰值情況下的應變率曲線Fig.9 Strain rate curves under different peak loads

（2）不同加載時長的影響。為了研究不同應力峰值情況下巖石的破壞規律，參考試驗中的加載波，數值模擬中的應力持續時間分別為150、175、200、225、250、275、300、325和350 μs。

分別進行不同加載時長下的動載SHPB 試驗模擬，可得到入射波、反射波和透射波曲線。得到的應力-應變曲線圖見圖10。從圖10看出，加載時長對動載抗壓強度峰值基本無影響。從圖11應變率曲線可以看出，應力波加載時長對試樣應變率有明顯的影響，在同種應力幅值的條件下，試樣破壞時最大應變率隨載荷的持續時間加長而增大。當應力峰值為350 MPa的情況下，減小載荷的持續時間到200 μs，試樣不會出現破壞現象，可見，單次沖擊荷載作用下，若低于某一應力峰值和持續時間，則試樣不會被破壞。這是因為應力峰值較小時，試樣將仍處于彈性狀態，經歷沖擊荷載后試樣可以恢復至原有形態。

圖10 不同荷載時長情況下的應力-應變曲線Fig.10 Stress-strain curves under different load duration

圖11 不同荷載時長情況下的應變率曲線Fig.11 Strain rate curve under different load duration

5 結 論

（1）通過采用Matlab 與LS-DYNA 軟件結合的數值計算方法，可高效地實施優化RHT模型參數的正交試驗。研究人員或工程人員只需輸入室內SHPB 試驗結果等信息，即可由程序自動完成數據處理、仿真計算和優化分析，最終得出優化后的RHT 模型參數，提高RHT 模型參數確定的效率，具有一定的工程實際意義。

（2）在基于深中通道沉管隧道項目花崗巖的室內SHPB 試驗的基礎上，以Matlab 和LS-DYNA 數值模擬SHPB 沖擊試驗為手段，實現了不同應變率下RHT 模型參數的確定，得到了適用于該項目花崗巖的RHT 模型參數。所得模型參數的計算結果與實際試驗較為吻合，且該套參數能合理地反映應變率效應。