考慮剪脹性的面層基層混合銑刨料模型及其參數敏感性分析

中圖分類號：TU448 文獻標志碼：A

Abstract： In order to investigate the mechanical properties of face-based mixed milling planes and their ontological model that isnot unspecified yet.Basedonthe triaxial test，a nonlinear shearexpansion constitutive model was constructed based on the Duncan-Chang model and the Pastor-Zienkiewicz shear expansion equation， and thestifness matrix of this constitutive model was derived by applying the theoryof infra-elasticity.Moreover， the model was developed twice with the help of the UMAT subroutine of the ABAQUS finite element software platform，which verified the modelaccuracyandreliabilityof the model.Inaddition，theoverallsensitivityanalysis and the response surface analysis of the keyparameters of the model were also carried out.The results show thatat about 10% axial strain，the strength of the two materials gradually reaches the residual strength under dierent circumferential pressures，and the volumetric strain and axial strain show a tendency to decrease and then increase， and exhibit obvious shear expansion characteristics within acertain range.The triaxial tests were simulated using a second developed UMAT subroutine，and the results showed a high degre of agreement when compared to the test results. The study index deviatoric stresses （σ₁-σ₃） are particularly affected by the parameters φ₀ ， Δφ ，and K ， while they are relatively weakly affected by the parameters n ， φ_cr ， α ，and R_f ，and the axial strain becomes significantly affected by the parameters R_f ， K ，and φ₀ ，while the effects of n ， Δφ ， φ_cr ，and α appear to be more limited. The degree of sensitivity to {σ₁-σ₃} is φ₀ ， φ₀² ， Δφ ， φ₀Δφ ， KΔφ ， K ， K² ， Δφ² ， Kφ₀ in descending order， and the degree of sensitivity to ε₁ is R_f ， φ₀ ， R_f² ， K ， R_fφ₀ ， KR_f ， φ₀² ， Kφ₀ ， K² in descending order. The interactions of thekey parameters were analyzed through the response surface methodology for comparing the interaction effects between them，and it is found that they are not significant for （ σ₁-σ₃ ）is： φ₀gt;Δφgt;K ； the specific influence on （204號 is： R_fgt;φ₀gt;K .The research results can provide a theoretical basis for the engineering application of surface mixed milling planer and a reference for parameter selection.

Key Words：intrinsic model;shear expansion;hypoelasticity theory ;sensitivity analysis;response surface

隨著道路交通流量的激增與車輛荷載的加重，造成道路面層和基層損壞問題日益突出，給道路的安全、暢通帶來了嚴峻挑戰[1-4].在道路的后續翻新修補中，機械銑刨作業會產生面層瀝青與基層碎石混合的銑刨料，這種面基混合銑刨料因含有基層碎石料，在荷載作用下具有非線性、壓硬性、摩擦性及應力路徑相關性等特性[5-7]，此外，碎石料在低密度和高圍壓條件下具有剪縮性，而在高密度和較低圍壓條件下還具有剪脹性[8-9].因此，在研究面基混合銑刨料的力學響應及構建其本構模型時，必須既考慮其非線性，也反映剪縮與剪脹的復雜特性.

目前在巖土工程領域中應用廣泛的本構模型有：Duncan等[10-11]提出的 E-ν 型和用體積模量代替泊松比的E-B型增量彈性的Duncan-Chang模型;K-G模型通過將應力和應變分解為球張量和偏張量，建立了平均應力 dp 與體積應變 dε_v 廣義剪應力dq與廣義剪應變 之間的關聯；在K-G模型的基礎上，研究者們改進和擴展出了Domaschuk-Valliappan模型[13]、Naylor模型[14]、Izumi-Verruijt耦合模型（三參數模型）[15]、四川大學K-G模型[16]和清華非線性解耦K-G模型[17]等；Roscoe等[18]基于深入的試驗觀察、臨界狀態概念的引入以及加工硬化原理和能量方程等提出劍橋模型;Burland[19]對劍橋模型引入橢圓屈服面、相適應的流動準則，以塑性體應變為硬化參數以及基于試驗數據的分析和比較提出修正劍橋模型；沈珠江2在吸收了鄧肯模型和劍橋模型的優點后提出了沈珠江雙屈服面彈塑性模型（南水模型）；殷宗澤2在深入研究土體復雜應力應變關系的基礎上，對剪脹性與剪縮性、硬化與軟化特性以及工程實踐的需求等多方面因素的綜合考慮提出殷宗澤雙屈服面彈塑性模型;姚仰平等[22-23]依據土的臨界狀態理論、現代塑性力學理論和試驗數據分析提出了統一硬化（UH）本構模型;程展林等[24-25]提出了一種新的“三參量\"非線性剪脹模型，該模型可視為鄧肯-張模型的改進型模型等.現階段對于土體本構模型的研究已取得了豐碩的成果，但每種模型均存在一定的適用范圍和限制條件，而針對適用于面基混合銑刨料的本構模型還不明確，找到合適的本構模型對基混合銑刨料的工程再應用具有重要的意義.

本文在碎石料和面基混合銑刨料大型三軸試驗的基礎上，基于Duncan-Chang模型和Pastor-Zien-kiewicz剪脹方程構建出非線性剪脹本構模型，通過應用次彈性理論推導出該本構模型的剛度矩陣，借助ABAQUS有限元軟件平臺的UMAT子程序對該非線性剪脹模型進行了二次開發，并將模擬結果與試驗結果進行對比，以驗證該模型的準確性和可靠性.此外，還對非線性剪脹模型的參數進行了全體敏感性分析及關鍵參數響應面分析.研究成果可為面基混合銑刨料的工程應用提供理論依據及為參數選取提供參考．

1材料力學性能試驗

1.1試驗材料

試驗材料選用碎石料（作參照）和面基混合銑刨料（面層瀝青銑刨料和基層碎石料按照質量1：1比例混合），通過篩分比試驗測定，碎石料粒徑主要集中在 5～40mm ，面基混合銑刨料粒徑主要集中在0.075～40mm ，其顆粒相對密度分別為2.75、2.62，限制粒徑 d₆₀ 分別為 15.5mm 和 10.0mm ，不均勻系數分別為1.61、19.23，曲率系數分別為1.05、1.73.根據《土工試驗方法標準》（GB/T50123—2019）的規定，測得碎石料和面基混合銑刨料的最大干重度分別為17.39kN/m³ 和 18.31kN/m³ ，最小于重度分別為13.3kN/m³ 和 13.2kN/m³ ，含水率分別為 12% 和 10% 試驗材料樣品及其顆粒級配曲線如圖1所示.

圖1材料顆粒級配曲線及材料樣品Fig.1 Material grading curves and material samples

1.2試驗儀器及方案

試驗使用的儀器為大型三軸剪切試驗系統，主要有三軸壓力室、軸壓系統、側壓系統和孔隙水壓力測讀系統等.其主要技術指標為：圍壓 （σ_c） 0～3.0MPa ，軸向應力 （σ_a）0～21MPa ，軸向行程 （U_a）0～ 300mm 試驗采用固結排水剪切試驗（CD），試樣為圓柱形，其尺寸采用直徑 300mm× 高度 600mm

針對碎石料、面基混合銑刨料，分別開展圍壓為 100kPa.200kPa.300kPa.400kPa 的大型三軸剪切試驗.試樣制備、試驗前準備和加載試驗工作嚴格按照《土工試驗方法標準》（GB/T50123—2019）的規定執行.將級配均勻的試驗材料裝填入預先準備好的筒內，分5層均勻攤鋪并振動夯實，確保壓實度均達到 96% ；扎緊橡皮膜，安裝加壓帽并將其抽至真空；施加圍壓開始試驗，剪切過程中以 0.4mm/min 的速率進行加載，直到試樣剪切至軸向應變的 25% 時或試樣完全破壞，停止試驗.其部分試驗過程如圖2所示.

圖2試樣制備及試驗 Fig.2 Specimen preparation and testing

1.3試驗結果

圖3為碎石料和面基混合銑刨料的偏差應力-軸向應變曲線及體積應變-軸向應變曲線.由圖可知，不同圍壓下兩材料偏差應力-軸向應變關系曲線的變化趨勢較為相似，均具有明顯偏差應力峰值特征，隨著軸向加載不斷增加，表現出彈性變形階段、塑性變形階段和試樣破壞階段；而體積應變-軸向應變關系曲線在加載過程中呈現出體積應變隨軸向應變增加而先減小后增大的特征，即在一定范圍內先出現剪脹效應，隨后受到外部作用力影響逐漸恢復并增大體積應變，整體呈現出非線性的曲線特征.

從圖3（a）、（c）可知，從試驗開始加載到試樣軸向應變 ε₁ 達到0.02左右時，兩材料的偏差應力與軸向應變成正比，其關系曲線基本呈直線，在此階段可近似為線彈性材料；當兩材料的變形超過彈性軸向應變臨界點后，偏差應力的增長速率要小于軸向應變增加速率，表現為 $＼mathscr { E } _ { 1 } { ＼sim } ＼Big （ ＼sigma _ { 1 } - ＼sigma _ { 3 } ＼Big ）$ 關系曲線開始不斷向軸向應變軸偏移，兩材料強度接近極限強度時，偏差應力的增長速率迅速變小，相反軸向應變增加速率急劇增大，塑性特性更顯著；當兩材料抗剪強度達到極限時，試樣將不能承受更大的加載，此時試樣發生剪切破壞，達到的強度值即為峰值抗壓強度，與之相對應的應變值則為峰值應變.此外，當軸向應變達到0.01左右時，兩材料強度趨于一個定值，即殘余強度，且碎石料的殘余強度略大于面基混合銑刨料（約 5% ），低圍壓的殘余強度小于高圍壓的殘余強度（25%～50%） 兒

從圖3（b）3（d）可觀察到，隨著圍壓從 100kPa 逐漸增大至 400kPa ，碎石料及面基混合銑刨料的體積應變隨軸向應變的增大而呈現不同程度的增長，且兩材料均展現出了典型的先剪縮后剪脹行為.圍壓的增加導致剪切初期體積應變值顯著上升，同時剪縮段的長度也相應增加.值得注意的是，當圍壓增大時，碎石料及面基混合銑刨料在后期的剪脹段呈現減小趨勢，其體積應變與軸向應變的關系曲線逐漸向剪縮方向轉換，特別在圍壓達到 400kPa 時，這兩種材料在剪切后期表現出顯著的剪縮特性，這與傳統粗粒土剪切試驗所揭示的規律相吻合.同時，表明面基混合銑刨料與碎石料的力學特性基本一致.

圖3不同圍壓下兩材料三軸試驗結果

Fig.3Triaxial test results of two materials under different enclosurepressures

2非線性剪脹本構模型

2.1 Duncan-Chang模型

為了本文理論模型的完整性，仍對Duncan-Chang模型作簡單介紹[1o].Duncan-Chang模型是一種建立在增量廣義胡克定律基礎上的非線性彈性模型，由鄧肯等人根據雙曲線應力應變關系提出.它能夠反映土體變形的非線性特征，模型參數物理意義明確，且可通過靜三軸試驗確定.其切線變形模量 E_t 的表達式為：

式中： K 為彈性模量基數； n 為彈性模量指數； R_f 為破壞比； p_a 為大氣壓， p_a=101.4kPa;c 和 φ 分別為土體的黏聚力與內摩擦角，對于粗粒土，內摩擦角對限制壓力的依賴關系可用式（2）表示.

式中： φ₀ 為參考圍限壓力下的內摩擦角； Δφ 為內摩擦角隨圍壓的減小量.

2.2Pastor-Zienkiewicz剪脹方程

在巖土工程領域，盡管Duncan-Chang模型在描述土體強度變形特性方面具有諸多優點，但其無法充分反映土體的剪脹性質.通過對碎石料和面基混合銑刨料進行的力學性能試驗研究表明，這2種材料在低圍壓條件下均表現出明顯的剪脹性.為了更準確地描述和預測兩材料的剪脹行為，有必要引入專門的剪脹方程，從而更精確地預測其在受剪過程中的力學響應.

關于土體剪脹模型及理論研究成果較多.其中幾種在工程實踐中尤為常用，Nova等[26]基于土體應力狀態和顆粒結構，引入剪脹系數等參數，建立了剪切與體積應變間的定量關系.Rowe[27]研究了剪脹比與應力比、初始孔隙比的關系，提出剪脹方程.徐舜華等[28]在Rowe基礎上進行了改進和方程拓展.Man-zari等[29]和Gajo等[30]基于臨界狀態下土體特性提出了剪脹方程.這些方程在描述特定條件下的土體剪脹行為時具有優勢，但也存在局限性.Nova剪脹理論未充分考慮巖土材料達到臨界狀態時的特性，在模擬接近或達到臨界狀態時會出現偏差.Rowe剪脹方程及徐舜華修正方程在高圍壓下可能低估材料壓縮性，而在剪脹階段可能高估其剪脹性，導致計算結果精確度偏低.Manzari和Gajo的剪脹方程包含較多參數，確定這些參數需依賴大量試驗數據和復雜擬合過程，不利于實際工程應用.

為了方便模型參數計算能精確反映模型剪脹性，本文選用Pastor等學者[31-32]提出的Pastor-Zien-kiewicz剪脹方程，該方程基于彈塑性理論并引入相關狀態參數，使其具有廣泛的適用性和更高的預測精度，且該方程中僅有2個參數，均可通過常規三軸試驗得到.Ravanbakhsh等[33]、Cola等[34]和Song等[35]應用該方程精準預測了砂土的三軸試驗結果.Heidarzadeh等[36]和Mira等[37]基于該方程的理論計算結果與試驗數據的對比驗證，認為該剪脹方程能夠準確反映體積應變和廣義剪應變增量之間的動態關系.該剪脹方程的一般表達式為：

式中： dε_v 和 dε_s 分別是剪脹引起增量體積應變和廣義剪應變； α 為試驗確定的模型參數； M 為 p-q 平面上臨界狀態線的斜率，在三軸壓縮應力狀態下可通過土體的殘余內摩擦角估算：

式中： φ_cr 為土體的殘余內摩擦角.

在三軸應力狀態下，平均應力 p 和廣義剪應力 q 的表達式分別為：

2.3非線性剪脹本構模型

程展林等[24]學者為揭示土體在剪切過程中的剪脹行為特性，將土的壓硬性描述為 p 對偏應變 ε_s 的影響，同時，將剪脹性闡釋為 q 對體應變 ε_v 的影響.鑒于這種交叉影響，傳統的廣義胡克定律無法準確描述這種復雜的應力-應變關系，因此，提出了粗粒土非線性剪脹模型，其應力-應變關系增量形式為：

式中： dε_v 和 dε_s 分別為體積應變和廣義剪應變的增量； dp，dq 分別為平均應力增量和廣義剪應力增量; 分別稱為體變模量、剪脹模量、剪切模量，以下對參數， G$ 進行具體推導.

2.3.1 體變模量 K_p

基于胡克定律的基本原理，可以得出體變模量K_p 的具體表達式為：

式中 ?：μ^Γ 為材料彈性階段的泊松比，假定其為常數[38]；E_ur 為土體的切線彈性模量，可由式（8）確定，

式中： K_ur 為試驗確定的模型參數.

2.3.2剪脹模量 K_q

通過Pastor-Zienkiewicz剪脹方程并結合三軸應變關系，由剪脹作用引起的體積應變增量 dε_v^q 可表示為：

在三軸試驗應力條件下，有：

聯立式（10）和式（11），可得剪脹作用引起的軸向應變 dε₁^q 為：

dq=dσ₁

聯立式（9）式（12）和式（13），可得剪脹作用引起產生的剪脹模量 K_q

2.3.3剪切模量 G

將式（9）代入式（3），可以得到剪脹作用引起的廣義剪應變增量 dε_s^q 的表達式為：

在三軸應力條件下，可以推導出由彈性作用引起的廣義剪應變增量 dε_s^e 的表達式為：

廣義剪應變增量 dε_s 的表達式為：

dε_s=dε_s^e+dε_s^q

將式（15）～式（17）代入式（6b），可以推導出剪切模量 G 的具體表達式為：

3模型三軸應力空間推廣

3.1剛度矩陣推導

為了推廣該模型在三軸應力空間下的更廣泛應用，充分利用非線性剪脹模型具有應力增量張量和應變增量張量之間一一對應的彈性關系特性，引入了次彈性理論[12（也稱作最小彈性模型）.假設研究材料為各向同性，且變形是瞬時的，不隨時間累積.則次彈性理論的三維關系張量形式為：

式中： dσ_ij，dε_kl 分別為增量應力張量和增量應變張量； 為次彈性理論的剛度矩陣.次彈性理論轉換成矩陣形式可寫為：

假設所研究的材料具備各向同性特性，則 的一般形式為：

式中： β₁，β₂，…，β₁₂ 為12個剛度矩陣參數，它們僅受應力張量的不變量所影響，其他因素對其并無直接作用.δ為kroneckerdelta函數，當 i=j 時， δ_ij=1 ；當i≠j 時 δ_ij=0

由于常規試驗無法全面確定所有12階剛度矩陣參數，潘家軍等[38提出了假設：認為參數 β₅～β₁₂ 與高階應力張量不變量有關，其對于剛度矩陣的影響相對較小，因此可以將參數 β_s～β_?12 設為0.基于這一假設，式（21）得以簡化為：

β₃σ_ijδ_kl+β₄δ_ijσ_kl

將式（22）代入式（19），增量應力張量的表達式可寫為：

dσ_ij=β₁δ_ijdε_kk+2β₂dε_ij+β₃σ_ijdε_kk+β₄δ_ijσ_ijdε_kl

對比式（20）和式（23），可以推導出剛度矩陣的各個元素具體表達式為：

在上述矩陣中，未明確寫出的其余元素均默認為0.

3.2矩陣參數

求解式（6a）（6b）的逆表達式，可以表示為：

dq=Gdε_s

式中： dp 和 dq 分別為土體的平均主應力增量和廣義剪應力增量.將增量應力張量的表達式（23）展開，可寫為：

dσ₁₁=β₁dε_kk+2β₂dε₁₁+β₃σ₁₁dε_kk+

dσ₂₂=β₁dε_kk+2β₂dε₂₂+β₃σ₂₂dε_kk+

dσ₃₃=β₁dε_kk+2β₂dε₃₃+β₃σ₃₃dε_kk+

在三軸應力條件下，由于 dε₂₂=dε₃₃ 和σ₂₂=σ₃₃ ，故 dp 和dq有：

dq=dσ₁₁-dσ₃₃=β₃qdε_v+3β₂dε_s

針對式（27）參數 β₁ 中 p_kl ，當下標 k=l=1 時，p_kl=σ^₂₂+σ₃₃-2σ^₁₁ 當下標 k=l=2 時， p_ij=σ₁₁+ σ₃₃-2σ₂₂ ；當下標 k=l=3 時 ，p_ij=σ₁₁+σ₂₂-2σ₃₃ 此時便可以得出次彈性理論的剛度矩陣，從而使非線性剪脹模型在三軸應力空間下得到進一步推廣.

4理論與試驗結果的驗證

4.1理論模型參數

針對不同圍壓條件下兩材料的大型三軸剪切試驗結果，進行數據處理與分析得出Duncan-Chang模型非線性彈性參數：彈性模量基數 K， 卸載再加載時的彈性模量基數 K_ur. 彈性模量指數 n 、彈性泊松比 μ 及破壞比 R_f. 強度參數：黏聚力 ∣c∣ 、參考圍限壓力下的內摩擦角 φ₀ 、內摩擦角隨圍壓的減小量 Δφ 及殘余內摩擦角 φ_cr. 為了獲取Pastor-Zienkiewicz剪脹方程中與剪脹相關的模型參數 α ，采用式（9）進行反推計算獲得，見式（28）.具體的理論模型參數，如表1所示.

聯立求解式（24）式（26），可得各項剛度矩陣表達式的參數分別為：

其中， dε， 和 iε₁ 分別為大型三軸剪切試驗中測量的體積應變與軸向應變.

表1理論模型參數

Tab.1 Theoretical model parameters

4.2模型驗證

為了驗證所推導非線性剪脹模型的有效性和三軸應力空間下次彈性理論剛度矩陣的適用性，基于ABAQUS有限元軟件平臺的UMAT子程序，針對具有剪脹特性的材料（如碎石料和面基混合銑刨料），開發出了能在三軸應力狀態下運行的非線性剪脹模型計算子程序.在驗證該模型和UMAT子程序準確性的過程中，嚴格模擬了大型三軸剪切試驗中的加載條件、模型尺寸和加載路徑，確保模擬結果的精準度和真實性.具體細節包括：計算模型為直徑300mm 、高度 600mm 的圓柱體，分析流程遵循先固結后剪切的步驟；隨后，運用該UMAT子程序進行模擬計算，并將模擬結果與試驗結果進行細致比對，如圖4所示.

值得注意的是，該模型在針對具有剪脹特性的粗粒土材料時，展現出廣泛的適用性，能夠精確預測其強度特征及變形特性.然而，該模型也存在一定的局限性，即該模型主要用于預測粗粒土峰值強度及其之前階段的應力-應變關系，未能涵蓋從峰值強度衰減至殘余強度的材料軟化過程.具體而言，當粗粒土的強度條件符合Mohr-Coulomb準則時，即標志著其達到峰值強度，此時該模型的UMAT子程序將自動停止計算過程，這意味著該模型對于峰值強度之后的狀態無法進行預測.然而，這一限制對實際工程應用的影響相對有限，因為出于結構安全性的首要考慮，工程分析往往僅關注峰值強度及其之前階段的應力-應變關系.此階段內土體的力學特性是評估土體結構變形、穩定性以及承載能力的關鍵依據：

圖4的對比結果表明，非線性剪脹模型在模擬不同圍壓條件下，碎石料和面基混合銑刨料的偏應力-軸應變及體積應變-軸應變關系曲線表現出色.該模型精確地預測了兩材料在受力過程中的非線性剪脹特性.與試驗數據相比，對于面基混合銑刨料，模型在 100kPa?200kPa?300kPa 和 400kPa 圍壓下，峰值偏應力的相對誤差分別為 8.8%4.6%?3.7% 和4.4% ，而峰值體積應變的相對誤差則在 4.7%～16.5% 之間波動.對于碎石料，模型在不同圍壓下的峰值偏應力和峰值體積應變的相對誤差均值分別為 6.3% 和 19.1%. 這些數據充分證明了模擬計算結果與實驗結果的良好一致性，進一步驗證了非線性剪脹模型的準確性、剛度矩陣的合理性以及子程序編制的正確性.

此外，鑒于該模型預測試驗結果所存在誤差的原因，主要歸結于以下幾個方面：1）土體的各向異性和非均質性.土體本身具有各向異性和非均質性，而該模型假設土體是均質和各向同性的，因此會引發一定的誤差；2）試驗本身和模型參數確定誤差.三軸試驗過程中可能存在操作誤差，如試樣制備、飽和度控制、試驗操作等，且模型參數如 K?n?μ?K_ur?R_f?c 、φ₀，Δφ，α 和 φ_cr 等，均通過三軸試驗確定，但這些參數的測定可能受到試驗方法、試驗條件等因素的影響，這些因素都可能引入誤差，影響模型參數的準確性；3）土體應力路徑的多樣性.實際土體在受力過程中經歷復雜的應力路徑，而該模型主要基于軸向應力增加的情況，可能無法準確模擬側向應力變化下的土體行為.

5模型參數敏感性分析

選用面基混合銑刨料為例作為參數敏感性分析對象.在進行非線性剪脹模型參數敏感性分析時，考慮到涉及的模型參數較多（共10個參數），難以全面考量所有因素間的交互作用39，為此，首先將全體參數運用方差分析法對其進行初步篩選，以識別出對模型具有顯著敏感性的關鍵參數；隨后，針對這些篩選出的關鍵參數，進一步實施響應面分析，從而評估出模型的參數敏感度.

5.1全體參數敏感性分析

本節以面基混合銑刨料的強度為例，作參數敏感性分析.在該模型的全體參數的選擇上，Pastor-Zienkiewicz模型包含2個參數 α 和 φ_cr ，而Duncan-Chang模型涉及8個參數，分別是 K，n，μ，K_ur，R_f°，φ₀， Δφ .需要注意的是，由于面基混合銑刨料的 c 為0kPa ，而卸載再加載時的彈性模量基數 K_ur 和彈性泊松比 μ 均為常數，且經大量計算發現 K_ur 和 μ 對模型強度影響微小，因此不把參數 c，K_ur 和 μ 作為全體參數敏感性分析的討論對象，則在全局參數敏感性分析中，共需考慮7個參數.

針對 K，n，R_f，φ₀，Δφ，φ_cr，α 共7個模型參數進行研究因子設計，研究指標為圍壓 100kPa 時所對應的（σ₁-σ₃） 和 ε₁. 在研究因子設計水平中，7個模型參數被視為7個研究因子，以考察其在不同取值下對研究指標的影響.設置低（-1）和高（1）兩個水平因子進行研究，低水平設定為模型參數的 90% ，而高水平則設定為模型參數的1.1倍.各研究因子及水平設計如表2所示．

表2研究因子水平設計Tab.2 Research factordesign levell

為了有效地分析出7個參數對研究指標的敏感性程度，設計出了一個包含14種工況的部分因子數值計算方案.詳細方案及計算結果見表3.

在模型全體參數的敏感性分析中采用了方差分析法，這是一種能夠系統而全面評估各因子對研究指標影響程度的有效方法，通過該方法能夠精確地確定模型參數的顯著性，以及因子間交互作用的顯著性，方差分析結果見表4.在分析過程中，將 P 值作為顯著性水平的度量，提供了一個概率框架來評估研究因子的交互效應.具體來說，當 P 值小于0.01時，認為該研究因子具有高度的顯著性，并用**進行標記；若 P 值在0.01～0.05之間，則判定該研究因子顯著，用*進行標注；當 P 值介于0.05和0.1之間時，稱該研究因子為邊緣顯著；而 P 值大于或等于0.1時，則意味著該研究因子不顯著.

從表4中可以看出，全體模型參數對 （σ₁-σ₃） 的敏感性排序為 φ₀，Δφ，K，n，φ_cr，α，R_f， 其中 φ_?0 具有高度的顯著性， Δφ 具有顯著性.對 ε₁ 的敏感性排序為R_f?K?φ₀?π?Δφ?φ_cr?α ，其中對 ε₁ 影響高度顯著的是 R_f. 綜合來看，各參數對 （σ₁–σ₃） 和 ε₁ 的敏感性排序有較大差異，對于研究指標 （σ₁-σ₃） ，受到參數 φ₀，Δφ，K 的影響尤為突出，相比之下，參數 n，φ_cr，α，R_f 的影響則相對較弱；然而，當關注研究指標 ε₁ 時，參數 R_f，K， φ₀ 的影響變得顯著，而 n，Δφ，φ_cr，α 的影響則顯得較為有限.

5.2關鍵參數響應面分析

通過全體參數敏感性分析，可知各參數對 （σ₁–σ₃） （2號和 ε₁ 的敏感性排序存在顯著差異，鑒于這種差異性，分別針對 （σ₁–σ₃） 和 ε₁ 的關鍵參數開展響應面分析，以更精確地揭示它們之間的相互作用和影響機制.

表3研究因子設計方案

Tab.3 Research factor design program

表6研究因子設計方案（偏差應力）

5.2.1偏差應力

針對 φ₀，Δφ，K 這3個關鍵參數進行響應面分析，以探究它們之間兩兩交互作用對 （σ₁-σ₃） 的影響，在此不研究三階及更高階的交互效應（影響極小）.每個關鍵參數設計三個水平：高水平（1）為模型參數的1.1倍，而低水平（-1）則為模型參數的 90% 研究因子與水平取值見表5.

表5研究因子與水平取值（偏差應力）

Tab.5Researchfactorsandlevelstaking（deviationstress）

根據所設定的水平值制定了研究因子設計方案，并將相應的計算結果匯總在表6中。

根據各試驗方案的計算結果計算 K，φ₀，Δφ 這3個參數以及參數之間二階交互效應的方差大小.對表6數據結果進行方差分析，偏差應力的方差分析結果如表7所示.可以看出，對偏差應力的敏感性程度由高到低次序為 φ₀，φ₀²，Δφ，φ₀Δφ，KΔφ，K，K²， Δφ²，Kφ₀ ，其中 φ₀，φ₀²，Δφ 為高度顯著， φ₀Δφ 為顯著.

表4部分研究因子理論方差分析結果

Tab.6Research factor design program （deviation stress）

為了深入剖析參數 φ₀，Δφ 和 K 之間交互作用對（2號 （σ₁–σ₃） 的具體影響程度，通過繪制各參數間響應面圖及等值線圖（如圖5所示）來分析其內在關系.由圖5（a）（b）可知，當保持 Δφ=8.35^° 時，無論參數 K 處于高水平還是低水平時， （σ₁-σ₃） 隨著 φ₀ 的增大而增大，并且在這一過程中，等值線的分布也變得更加密集；而在 φ₀ 處于任意水平時， （σ₁–σ₃） 幾乎不會隨著 K 的增大而發生顯著變化，說明 φ₀ 的變化對 （σ₁-σ₃） 的變化影響更加顯著，顯著性隨著 φ_?0 的減小而減弱，相比之下， K 的變化對 （σ₁-σ₃） 的影響幾乎可以忽略不計.

表7偏差應力的方差分析

Tab.7Analysisofvariancefordeviatoricstress

由圖5（c）、（d）可知，在保持 K=365.59 時，當 φ₀ 處于任意水平時， （σ₁–σ₃） 幾乎不受 Δφ 變化的影響；而當 Δφ 處于任意水平時， φ₀ 的變化對 （σ₁-σ₃） 的變化影響更加顯著，顯著性隨著 φ₀ 的減小而減弱.由圖5（e）、（f）可知，當保持 φ₀=46.47^° 時， K，Δφ 的水平的變化對偏差應力 （σ₁-σ₃） 的變化影響均不顯著，隨著 K，Δφ 的增大或減小， （σ₁-σ₃） 基本無變化.

5.2.2軸向應變

采用響應面分析法探究參數 R_f，K 和 φ₀ 之間兩兩交互作用對 ε₁ 的具體影響.同樣，設計三個水平，高水平（1）取模型參數的1.1倍，低水平（-1）取模型參數的 90% .研究因子與水平取值見表8，研究因子設計方案見表9.

根據研究因子設計方案對各關鍵參數以及參數之間二階交互效應的方差進行計算，方差分析結果詳見表10.可以看出，對軸應變的敏感性程度由高到低次序為 R_f°φ₀°R_f²°K°φ₀°KR_f°φ₀ ²°Kφ₀°K² ，其中 R_f? φ₀，R_f²?K，R_fφ₀ 高度顯著， KR_f、φ₀² 顯著.

表8研究因子與水平取值（軸向應變）

Tab.8Researchfactorsandlevelvalues（axialstrain）

表9研究因子設計方案（軸向應變）

Tab.9Research factor design program （axial strain）

表10軸向應變的方差分析

Tab.1oAnalysisofvarianceforaxialstrain

鑒于參數 R_f，K 和 φ₀ 之間對 ε₁ 有顯著的交互作用，通過繪制響應面圖及等值線圖（如圖6所示）來深入剖析其內在關系.由圖6（a）（b）可知，當 φ₀ 保持在 46.47^° 時，若 K 處于低水平，觀察到 ε₁ 會隨著 R_f 的增大而增大，且此時等值線分布密集，說明此時 R_f 的變化對 ε₁ 的影響更為顯著，而這種顯著性會隨著 K 的減小而逐漸減弱；若 R_i 處于高水平時，隨著 K 的增大， ε₁ 呈現減小的趨勢，且此時等值線的分布變得相對稀疏，說明此時 K 的變化對 ε₁ 的影響不顯著，顯著性隨著 R_f 的減小而減弱.此外，值得注意的是，云圖呈現出弧度過渡的形態，而非簡單的線性過渡，證明參數 R_f 和 K 間存在顯著的交互作用，并共同影響著ε₁ 的變化.

圖6關鍵參數兩兩交互作用對軸向應變的影響

Fig.6 Influence of pairwise interaction of key parameters on axial strain

由圖6（c）（d）可知，當 K 保持在365.59時，當 φ₀ 處于高水平時， R_f 的變化對 ε_?1 的影響更為顯著，這種顯著性隨著 φ₀ 的減小而逐漸減弱.然而，在 R_f 處于低水平的情況下， φ_?0 的變化對 ε₁ 的影響則相對不顯著，并且隨著 R_f 的增大，其顯著性反而增強.

由圖6（e）（f）可知，當 R_f 保持在0.81時，在 K 處于低水平時， φ₀ 的變化對 ε₁ 的影響更加顯著，顯著性隨著 K 的增大而減弱；在 φ₀ 處于高水平時， K 的變化對 ε₁ 的影響更加顯著，顯著性隨著 φ₀ 的減小而減弱，K 與 φ_?0 之間存在一定的交互作用，且共同影響 ε₁ 的變化.

6結論

1）在軸向應變為 10% 左右時，不同圍壓下兩種材料的強度都逐漸達到殘余強度，其中碎石料的殘余強度略高出面基混合銑刨料約 5% 此外，低圍壓條件下的殘余強度相較于高圍壓條件要低，差距在25%～50% 之間.體積應變與軸向應變之間的關系呈現出一種先減小后增大的趨勢，且在一定范圍內表現為明顯的剪脹特性.

2）在三軸應力狀態下，利用ABAQUS有限元軟件平臺的UMAT子程序對文中非線性剪脹模型進行二次開發，該模型計算結果與試驗結果相比展現出高度的吻合度，驗證了模型的準確性和有效性.

3）通過全體參數敏感性分析發現，對于研究指標 （σ₁-σ₃） ，受到參數 φ₀，Δφ，K 的影響尤為突出，相比之下，參數 的影響則相對較弱;然而，對于研究指標 ε₁ ，參數 R_f?K，φ₀ 的影響變得顯著，而 n 、Δφ，φ_cr，α 的影響則顯得較為有限.

4）通過關鍵參數及其二階效應敏感性分析，發現對 （σ₁-σ₃） 的敏感性程度由高到低次序為 Δφ，φ₀Δφ，KΔφ，K，K²，Δφ²，Kφ₀ ，其中 φ₀，φ₀²，Δφ 為高度顯著， φ₀Δφ 為顯著.對軸應變的敏感性程度由高到低次序為 R_t?φ₀?R_f²?K?R_fφ₀?KR_f?φ₀²?Kφ₀?K²， 其中R_f、φ₀、R_f²、K、R_fφ₀ 高度顯著， KR_f、φ₀² 顯著.

5）無論 K，Δφ 的增大或減小， φ₀ 的變化對 （σ₁- σ₃ 的影響更加顯著.當 K 處于低水平時， R_i 的變化對ε₁ 的影響更為顯著；當 φ₀ 處于高水平時， R_i 的變化和K 的變化對 ε₁ 的影響更加顯著.

6文中開發出的面基混合銑刨料的非線性剪脹模型，不僅適用于具有剪脹性巖土材料的數值分析，還能為處理復雜的實際工程問題提供高效、實用的解決方案.另外，該模型有望通過等效剛度和等效強度的方法進一步拓展至土工合成材料加筋粗粒土結構，從而推動加筋土結構數值計算領域的發展，

參考文獻

[1] 趙曉華，楊海益，姚瑩，等．高速公路橋梁路段交通安全風險評 價及影響要素挖掘[J]．同濟大學學報（自然科學版），2022，50 （11） ：1637-1646. ZHAOXH，YANGHY，YAOY，etal.Traffic safety evaluation and feature analysis of freeway bridge sections[J].Journal of Tongji University （Natural Science），2022，50（11） ：1637-1646. （in Chinese）

[2] 楊曉光，胡仕，張夢雅．智能高速公路交通應用技術發展綜述 [J]．中國公路學報，2023，36（10）：142-164. YANG XG，HU S，ZHANG MY.Development of intelligent motorwaytrafficapplicationtechnologies：areview[J].China JournalofHighway and Transport，2023，36（10）：142-164. （in Chinese）

[3] 裴玉龍，何慶齡，侯琳，等.基于文獻計量分析的交通量預測 研究綜述[J].公路交通科技，2023，40（2）：171-181. PEI YL，HEQL，HOUL，et al.Review of studies on traffic volumepredictionbasedonbibliometricanalysis[J].Journal of Highway and Transportation Research and Development，2023， 40（2）：171-181.（in Chinese）

[4] 陳昌富，李偉，張嘉睿，等．山區公路邊坡工程智能分析與設 計研究進展[J]．湖南大學學報（自然科學版），2022，49（7）： 15-31. CHENCF，LI W，ZHANG JR，et al．State-of-the-art of intelligent analysis and design in slope engineering of highways in mountainous areas[J]. Journal of Hunan University（Natural Sciences），2022，49（7）：15-31.（in Chinese）

[5] YUL T，XIEJ，LIR，et al. Study on the performance of emulsified asphalt recycled subgrade based on the evaluation of semi-rigid milling material[J].Construction and Building Materials，2022， 324：126614.

[6] 童海濤，郭知濤，秦凱強．高速公路改擴建銑刨料再生利用技 術案例分析[J]．公路，2020，65（3）：243-247. TONG H T，GUO Z T，QIN K Q.Case study on recycling technologyofmillingmaterialsinexpresswayreconstructionand expansion[J].Highway，2020，65（3）：243-247.（in Chinese）

[7] 田源，牛冬瑜.銑刨料摻量及摻配方式對再生水穩碎石路用性 能的影響[J].長安大學學報（自然科學版），2020，40（4）： 39-49. TIANY，NIUD Y.Efect of milling material content and blending method on road performance of regeneration cement stabilized gravel[J].Journal of Chang’an University（Natural Science Edition），2020，40（4）：39-49.（in Chinese）

[8] VERMA G，KUMAR B． Prediction of compaction parameters for fine-grained and coarse-grained soils：a review[J]. International Journal ofGeotechnical Engineering，2020，14（8）：970-977.

[9] 蔣明杰，吉恩躍，王天成，等．粗粒土抗剪強度的縮尺效應規律 試驗研究[J]．巖土工程學報，2023，45（4）：855-861. JIANGMJ，JIEY，WANGTC，et al.Experimental studyon laws of scale effects of shear strength of coarse-grained soils[J]. Chinese Journal ofGeotechnical Engineering，2023，45（4）： 855-861.（in Chinese）

[10]DUNCANJM，CHANGCY.Nonlinear analysis of stress and strain insoils[J].Journal of the Soil Mechanics and Foundations Division，1970，96（5）：1629-1653.

[11]DUNCANJM，BYRNEP，WONG K S，et al.Strength，stressstrainand bulkmodulusparametersfor finiteelementanalysis of stress and movement in soil masses[R].In：Report No. UCB/GT/ 80-01，Berkerley：University of California，1980.

[12］李廣信．高等土力學[M]．北京：清華大學出版社，2004. LIGX.Advanced soil mechanics[M].Beijing：Tsinghua UniversityPress，2004.（in Chinese）

[13]DOMASCHUK L，VALLIAPPAN P.Nonlinear setlement analysis byfinite element[J].Journal of the Geotechnical Engineering Division，1975，101（7）：601-614.

[14]NAYLOR D J. Stress-strain laws for soils [M]. Essex：Applied Science Publishers Ltd.，1978.

[15］IZUMI N，VERRUIJT A.The dynamic interaction between water flow and soil deformation[J].Soilsand Foundations，1982，22（1）： 29-44.

[16］曾以寧，屈智炯，劉開明，等.土的非線性K-G模型的試驗研究 [J].成都科技大學學報，1985，17（4）：143-149. ZENGYN，QU ZJ，LIU KM，et al．Experimental study on nonlinear K-G model of soil[J]．Journal of Chengdu University of Science and Technology，1985，17（4）：143-149.（in Chinese）

[17］張丙印，于玉貞，張建民．高土石壩的若干關鍵技術問題[C]// 第八屆土力學及基礎工程學術研討會論文集.北京：中國建筑 工業出版社，2003：163-186. ZHANGBY，YUY Z，ZHANG JM.The key questions of high earth-rock dam[C]//Proceedings of the 8th Symposium on Soil Mechanics and Foundation.Beijing：China Architecture and Building Press，2003：163-186.（in Chinese）

[18]ROSCOE KH，SCHOFIELD AN，THURAIRAJAH A. Yielding of clays in states weter than critical[J].Géotechnique，1963， 13（3）：211-240.

[19］BURLAND JB.The yielding anddilation ofclay [J]. Geotechinique，1965，15（2）：211-214.

[20］沈珠江.理論土力學[M].北京：中國水利水電出版社，2000. SHENZJ.Theoretical soil mechanics[M].Beijing：China Water amp; Power Press，2000.（in Chinese）

[21］殷宗澤．一個土體的雙屈服面應力-應變模型[J]．巖土工程 學報，1988，10（4）：64-71. YIN Z Z. A double yield surface stress-strain model of soil[J]. Chinese Journal ofGeotechnical Engineering，1988，10（4）：64- 71.（in Chinese）

[22]姚仰平．UH模型系列研究[J].巖土工程學報，2015，37（2）： 193-217. YAO YP. Advanced UH models for soils[J].Chinese Journal of Geotechnical Engineering，2015，37（2）：193-217.（in Chinese）

[23]ZHU EY，YAO Y P.Structured UH model for clays[J]. Transportation Geotechnics，2015，3： 68-79.

[24］程展林，姜景山，丁紅順，等．粗粒土非線性剪脹模型研究[J]. 巖土工程學報，2010，32（3）：460-467. CHENG Z L，JIANG JS，DING H S，et al. Nonlinear dilatancy modelforcoarse-grainedsoils [J].Chinese Journalof Geotechnical Engineering，2010，32（3）：460-467.（in Chinese）

[25］程展林，陳鷗，左永振，等．再論粗粒土剪脹性模型[J]．長江 科學院院報，2011，28（6）：39-44. CHENG ZL，CHEN O，ZUO Y Z，et al.Further discussion on dilatancy model for coarse-grainedsoils[J]. Journal of Yangtze River Scientific Research Institute，2011，28（6）：39-44.（in Chinese）

[26]NOVAR，WOOD D M.A constitutive model for sand in triaxial uu AnalyticalMethods inGeomechanics，1979，3（3）：255-278.

[27]ROWE PW. The stress-dilatancy relation for static equilibrium of an assemblyof particles in contact[J].Proceedings of the Royal Society of London Series A Mathematical and Physical Sciences， 1962，269（1339）： 500-527.

[28］徐舜華，鄭剛，徐光黎．考慮剪切硬化的砂土臨界狀態本構模 型[J]．巖土工程學報，2009，31（6）：953-958. XUS H，ZHENG G，XU G L. Critical state constitutive model of sand with shear hardening[J].Chinese Journal of Geotechnical Engineering，2009，31（6）：953-958.（in Chinese）

[29]MANZARI M T，DAFALIAS Y F.A critical state two-surface plasticitymodel forsands[J].Geotechnique，1997，47（2）： 255-272.

[30]GAJO A，WOOD M. Severn-Trent sand：a kinematic-hardening constitutive model： the q-p formulation[J].Geotechnique，1999， 49（5）：595-614.

[31]PASTOR M，ZIENKIEWICZ O C，LEUNG K H. Simple model for transient soil loadingin earthquake analysis II.Non-associative models for sands[J]. International Journal for Numerical and AnalyticalMethodsinGeomechanics，1985，9（5）：477-498.

[32]PASTOR M，ZIENKIEWICZ O C，CHAN A H C. Generalized plasticity and the modelling of soil behaviour[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 1990，14（3）：151-190.

[33]RAVANBAKHSH E，HAMIDI A.Development of a generalized plasticity constitutive model for cemented sands using critical state concepts].IteatialJoalofotehalEig 2013，7（4）：364-373.

[34]COLA S，TONNI L. Adapting a generalized plasticity model to reproduce the stress-strain response of silty soils forming the Venice lagoon basin[M]/Solid Mechanics and Its Aplications. Dordrecht： Springer Netherlands，2007： 743-758.

[35]SONGF，JINYT，LIU HB，et al．Analyzing the deformation and failure of geosynthetic-encased granular soil in the triaxial stress condition[J].Geotextilesand Geomembranes，2O20，48（6）： 886-896.

[36]HEIDARZADEH H，OLIAEI M. Development of a generalized model using a new plastic modulus based on bounding surface plasticity[J]. Acta Geotechnica，2018，13（4）：925-941.

[37]MIRA P，TONNI L，PASTOR M，et al.A generalized midpoint algorithm for the integration of a generalizedplasticity model for sands[J].International Journal for Numerical Methodsin Engineering，2009，77（9）：1201-1223.

[38］潘家軍，程展林，饒錫保，等．一種粗粒土非線性剪脹模型的 擴展及其驗證[J]．巖石力學與工程學報，2014，33（增刊2）： 4321-4325. PAN JJ，CHENG Z L，RAO XB，et al.A nonlinear dilatancy model’s extention and validation for coarsed-grained soils[J]. Chinese Journal of Rock Mechanicsand Engineering，2014， 33（Sup.2）：4321-4325.（in Chinese）

[39］馮亞新，江兆強，孫一清，等．基于改進響應面的鄧肯-張E-B 模型參數全局敏感性分析[J]．水利水電科技進展，2023， 43（2）：44-50. FENGY X，JIANG ZQ，SUN YQ，et al.Global sensitivity analysisof parametersin Duncan-Chang E-B model basedon improved response surface [J]. Advances in Scienceand TechnologyofWaterResources，2023，43（2）：44-50.（in Chinese）