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

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

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

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

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

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

1 土壤本構(gòu)關(guān)系

土壤的應(yīng)力-應(yīng)變關(guān)系可以用固體力學(xué)中的彈塑性理論來描述。土壤的力學(xué)行為可以分為彈性行為和塑性行為。彈性行為是土壤在彈性變形范圍內(nèi)的應(yīng)力-應(yīng)變關(guān)系，可以用彈性本構(gòu)模型來描述。彈性模型主要包含兩個參數(shù)，彈性模量E和泊松比μ.常見的土壤彈性本構(gòu)模型包括線彈性模型、多孔介質(zhì)彈性模型和線粘彈性模型。筆者采用適用最廣泛的各向同性彈性模型描述土壤的彈性行為。

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

F=Rmcq-ptanφ-c=0，

(1)

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

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

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

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

表1 起伏路面土壤力學(xué)參數(shù)

2 土壤壓力沉陷關(guān)系非線性擬合

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

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

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

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

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

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

以表2中4種行駛速度所對應(yīng)的平板加載速率為例，在2.1節(jié)建模的基礎(chǔ)上，采用表1中2#土壤參數(shù)，分別設(shè)置不同的加載速率，得到v1～v44種加載速率下的土壤壓力沉陷關(guān)系曲線，并與準(zhǔn)靜態(tài)加載時的土壤壓力沉陷曲線對比，如圖2所示。

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

2.3 平板尺寸對土壤壓力沉陷關(guān)系的影響

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

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

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

3 土壤壓力沉陷關(guān)系擬合

3.1 Prony級數(shù)逼近方法

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

(2)

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

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

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

(3)

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

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

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

(4)

(5)

聯(lián)立式(2)～(5)可得

(6)

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

(7)

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

3.2 數(shù)值擬合

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

3.3 Kriging代理模型建立

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

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

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

表3 輸入?yún)?shù)范圍

為使得試驗樣本盡可能覆蓋整個樣本空間，采用最優(yōu)拉丁超立方的方法在樣本空間內(nèi)抽取80組試驗樣本?；贗sight-Abaqus聯(lián)合仿真，建立輸入?yún)?shù)與輸出參數(shù)之間的Kriging代理模型。并利用交叉驗證的方法對代理模型進行誤差分析，得到每個輸出參數(shù)響應(yīng)擬合，如圖5所示。

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

3.4 模型測試

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

表4 典型土壤壓力沉陷試驗參數(shù)

以表4中3種試驗條件為輸入?yún)?shù)，得到每種試驗條件下的Kriging代理模型輸出結(jié)果，并把仿真結(jié)果和Prony函數(shù)擬合結(jié)果對比，如圖6所示。

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

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

4 結(jié)束語

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