一種基于Kriging代理模型的土壤壓力沉陷特性預測方法

馮文選，吳大林，馬吉勝，何 健

(1.陸軍工程大學石家莊校區，河北 石家莊 050003；2.西京學院，陜西 西安 710123)

在車輛地面力學中，土壤作為一種材料，其力學性能是影響履帶式自行武器越野行駛性能的重要因素。履帶板對土壤施加垂直壓力，使土壤在壓力作用下產生沉陷，從而影響履帶式自行武器的行駛阻力、牽引力和履帶板滑轉率等性能。在軍事行動中，還會影響履帶式裝備在特定區域的通過性能。因此，研究土壤在垂直壓力作用下的沉陷特性對提高武器裝備越野行駛性能和通過性意義重大。

由于土壤材料應力應變關系的復雜性，使得通過嚴格的力學推導來研究土壤承壓沉陷特性十分困難。因此，早期的土壤承壓沉陷特性研究方法主要以純經驗法和半經驗法為主[1]。20世紀70年代，Prumpral J V等首次用非線性彈性問題有限元解法分析了土壤在拖拉機車輪作用下的變形和應力分布。Wong J Y在專著中[2]闡述了土壤應力與破壞的試驗觀測與有限元法的研究。計算機技術的發展使車輛地面力學的研究方法更加豐富多樣，潘衛東等對軟地面車輛系統的建模仿真方法進行了研究[3]，呂唯唯等研究了快速加載條件下可變形土壤的特性[4]，趙家豐等研究了重塑土壤壓力沉陷特性[5]。

綜上所述，在研究方法上，雖然土壤承壓沉陷特性的研究手段不斷豐富，但以經驗法為主的經典土壤承壓沉陷特性理論仍然是各國學者開展研究的理論基礎。工程實際應用中，經典土壤承壓沉陷特性理論需要以大型原位土壤承壓試驗為依據，由于開展大型原位試驗受試驗環境、設備和資金等因素的影響，使其存在諸多困難。筆者利用有限元仿真的方法，基于數值擬合方法對土壤壓力沉陷特性進行分析，并基于Kriging插值法建立土壤力學參數和加載條件與擬合函數參數之間的關系。

1 土壤本構關系

土壤的應力-應變關系可以用固體力學中的彈塑性理論來描述。土壤的力學行為可以分為彈性行為和塑性行為。彈性行為是土壤在彈性變形范圍內的應力-應變關系，可以用彈性本構模型來描述。彈性模型主要包含兩個參數，彈性模量E和泊松比μ.常見的土壤彈性本構模型包括線彈性模型、多孔介質彈性模型和線粘彈性模型。筆者采用適用最廣泛的各向同性彈性模型描述土壤的彈性行為。

土壤的塑性行為主要包括屈服、硬化、剪脹和流動性，可以用塑性本構模型來描述。經典的土壤塑性本構模型主要有Mohr-Coulomb模型、Drucker-Prager模型和臨界狀態塑性模型[6]。筆者采用被廣泛應用的Mohr-Coulomb模型描述土壤的塑性行為。Mohr-Coulomb模型中的剪切屈服函數為

F=Rmcq-ptanφ-c=0，

(1)

式中：φ為q-p平面上的摩擦角；c為黏聚力；Rmc為屈服函數在π平面上的形狀參數。

在ABAQUS軟件中，為了避免π平面上屈服面存在尖角導致流動方向不唯一的現象，Mohr-Coulomb模型采用非相關聯的流動準則。并通過指定黏聚力c與等效塑性應變之間的關系，控制土壤的硬化規律。

土壤本構參數可以通過土工試驗來確定。土體抗剪強度試驗是測定土體剪切特性的一種試驗，分為室內試驗和現場試驗。室內試驗包括直剪試驗和三軸剪切試驗，直剪試驗儀器設備簡單、操作方便、試驗用土少，可以測量土壤的抗剪強度參數粘聚力c和內摩擦角φ.三軸壓縮試驗作為另一種常見的室內土工試驗，與直剪試驗相比，除了可以測定抗剪強度指標粘聚力c和內摩擦角φ，還可以測定土壤的彈性模量E、泊松比μ和壓縮系數等參數[7]。

采用上述兩種土工試驗對起伏路面松軟土壤進行參數測定[8]，得到土壤的力學參數如表1所示。

表1 起伏路面土壤力學參數

2 土壤壓力沉陷關系非線性擬合

2.1 土壤承壓沉陷試驗有限元建模

有限元仿真模型采取三維實體建模。載荷板為半徑R=0.1 m的圓板。為簡化計算，載荷板采用剛體約束，并設置參考點。綜合考慮土體模型邊界對仿真結果的影響和仿真計算效率，土體采用1 m×1 m×1 m的正方體模型。建立材料分別為鋼和土壤的截面屬性，鋼的材料參數為：密度7 800 kg/m3、彈性模量E=207 GPa、泊松比μ=0.3；土壤本構模型采用各向同性線彈性模型和Mohr-Coulomb模型，參數如表1所示。

約束土體模型4個側面水平方向位移自由度，底面設置為完全固定約束；約束載荷板參考點水平方向位移自由度和3個方向轉動自由度。整個模型采用八節點線性六面體單元，設置減縮積分、沙漏控制和網格重劃分區域。土體網格拓撲結構采用中間致密四周相對稀疏的方法。采用顯示動態分析步，設置載荷板參考點豎直方向載荷，加載載荷為50 kN，模型裝配體如圖1所示。

2.2 加載速率對土壤壓力沉陷特性的影響

車輛越野行駛過程中，隨著車速的不斷提高，車輛對地面施加載荷的速率可以達到每秒幾兆帕。以美軍M4坦克為例，其最大行駛速度42 km/h，總重32 t，平均最大接地壓力Pmm=282.2 kPa[9].按照M4坦克履帶接地長度為5 m，寬為0.25 m，計算不同行駛速度對應的對地加載速率，如表2所示。

表2 M4坦克行駛速度、加載速率對比

以表2中4種行駛速度所對應的平板加載速率為例，在2.1節建模的基礎上，采用表1中2#土壤參數，分別設置不同的加載速率，得到v1～v44種加載速率下的土壤壓力沉陷關系曲線，并與準靜態加載時的土壤壓力沉陷曲線對比，如圖2所示。

從圖2中可以看出，動態加載時的土壤壓力沉陷關系與準靜態加載時存在明顯區別。對照表2中M4坦克行駛速度與加載速率的關系，可以得出：M4坦克行駛速度大于20 km/h時，同一沉陷深度，動態加載條件下的土壤承載力明顯高于準靜態條件下的土壤承載力。此外，4種加載速率對應的土壤壓力沉陷關系之間也存在差異，土壤承載力隨著加載速率的增大呈現提高趨勢。

2.3 平板尺寸對土壤壓力沉陷關系的影響

土壤平板載荷實驗中，由于bekker公式所依據的土壤平板載荷試驗采用緩慢的加載方式，無法反映土壤在高速載荷作用下的承壓沉陷特性，因此有必要研究平板尺寸在高速加載條件下對土壤壓力沉陷關系的影響規律。

在2.1節建模的基礎上，采用表1中2#土壤參數和表2中v1加載速率，并分別設置平板半徑為0.05、0.10、0.15和0.20 m 4種尺寸，得到土壤壓力沉陷關系曲線，如圖3所示。

從圖3中可以看出，平板尺寸對土壤壓力沉陷特性有顯著的影響，同一沉陷深度時，土壤承載力隨著平板尺寸的增大而減小。且在尺寸越小時，承載力減小的趨勢越快。承載力與平板尺寸之間成負相關關系，這一趨勢與bekker公式的平板寬度對土壤承載力的影響趨勢是一致的。

3 土壤壓力沉陷關系擬合

3.1 Prony級數逼近方法

Prony方法是以構造形如式(2)的指數型函數為目的的一種非線性數值逼近算法[10]。

(2)

式中：A，s為待定參數；e為自然底數；t為自變量；f為函數。

Prony方法的主要原理為：令t=iT,使得式(3)成立。

f(iT)=fi，i=0,1,…,2n-1，

(3)

式中：T為步長；{fi}為給定型值。

定義變量zj、αi，且分別滿足：

zj=esjT，j=1,…，n,

(4)

(5)

聯立式(2)～(5)可得

(6)

zn+αn-1zn-1+…+α1z+α0=0，

(7)

從而求出αi(i=0,1,…,n-1)和Aj.

3.2 數值擬合

以上研究可以得出，土壤的壓力沉陷關系屬于典型的非線性關系，不僅與土壤的力學性質有關，而且與加載速率和平板尺寸有關。而Prony級數具有良好的非線性逼近能力，因此利用Prony級數，通過MATLAB軟件，對不同土壤力學參數和試驗條件下的壓力沉陷曲線進行擬合。以表1中1#土壤參數為例，在平板尺寸r=0.05 m，加載速率為1 254.22 kPa/s的條件下進行有限元數值仿真，并采用Prony指數函數擬合，如圖4所示。

3.3 Kriging代理模型建立

由3.2節可見，Prony級數能夠對不同力學性質和試驗條件下的土壤壓力沉陷曲線進行良好的擬合，且n=2時已具有很高的擬合精度，R2均大于0.9.但是其擬合過程較為復雜，工程實際操作性不強。無法通過已知的土壤力學參數和加載條件直接預測出土壤的壓力沉陷特性曲線。而Kriging插值法[11]以變異函數理論與結構分析為基礎，在無需確定輸入與輸出之間精確解析關系的條件下，可以較好的逼近輸入與輸出參數之間的非線性關系。

以彈性模量E、泊松比μ、內摩擦角φ、黏聚力c、平板半徑r和加載速率v為輸入參數，以n=2時Prony級數P=aebz+cedz的參數a、b、c、d為輸出，建立Kriging代理模型。

輸入參數取值范圍如表3所示，其中土壤力學參數以起伏路面為參考，平板半徑采用小尺寸平板，加載速率以M4坦克20 km/h以上的行駛速度為參考。

表3 輸入參數范圍

為使得試驗樣本盡可能覆蓋整個樣本空間，采用最優拉丁超立方的方法在樣本空間內抽取80組試驗樣本。基于Isight-Abaqus聯合仿真，建立輸入參數與輸出參數之間的Kriging代理模型。并利用交叉驗證的方法對代理模型進行誤差分析，得到每個輸出參數響應擬合，如圖5所示。

從圖5中可得，誤差分析樣本點對每個輸出參數的預測值與真值之間擬合精度良好。通過復相關系數計算，得出各輸出參數R2依次為0.96、0.98、0.94、0.99.

3.4 模型測試

利用上述代理模型對起伏路土壤壓力沉陷特性進行預測。選取典型試驗條件如表4所示。

表4 典型土壤壓力沉陷試驗參數

以表4中3種試驗條件為輸入參數，得到每種試驗條件下的Kriging代理模型輸出結果，并把仿真結果和Prony函數擬合結果對比，如圖6所示。

從圖6中可得，Kriging代理模型對Prony指數函數具有良好的逼近效果，并且能夠很好地對不同土壤力學參數和加載條件下的土壤壓力沉陷曲線進行擬合。其擬合復相關系數均在0.96以上。

綜上所述，Prony指數函數能夠對不同加載速率和平板尺寸的土壤壓力沉陷特性曲線進行擬合，且擬合精度較高，但是Prony函數在土壤壓力沉陷關系擬合中的應用較為繁瑣，擬合過程需要依據土壤壓力沉陷仿真結果或原位試驗結果。以土壤力學參數、加載速率和平板尺寸為輸入，以Prony指數函數參數為響應，建立輸入與輸出之間的Kriging代理模型，可以在不需要仿真和試驗的條件下，建立不同土壤力學參數和加載條件下的土壤壓力沉陷特性模型，具有一定的工程應用價值。

4 結束語

通過有限元仿真，研究了加載速率和平板尺寸對土壤壓力沉陷特性關系的影響。加載速率的提高會增大同一沉陷深度時的土壤承載力，而在仿真范圍內，平板半徑的增大會導致土壤單位面積承載力的減小。從數值逼近的角度對土壤壓力沉陷特性曲線進行了研究，發現Prony指數函數擬合方法可以對不同加載速率和平板尺寸條件下的土壤壓力沉陷特性曲線進行擬合，且擬合精度較高。從實際應用的角度出發，采用Kriging插值法建立了土壤力學參數、加載速率和平板尺寸與Prony指數函數參數之間的代理模型，通過輸入土壤參數和加載條件可以直接預測特定條件下的土壤壓力沉陷特性關系，避免了進行高精度數值仿真和原位試驗的過程。對研究土壤承壓沉陷特性具有一定的借鑒意義。