摘要：顆粒離散元的模擬精度主要取決于所標(biāo)定細(xì)觀參數(shù)的準(zhǔn)確性，為提升細(xì)觀參數(shù)的標(biāo)定精度和效率，本文首先采用抗轉(zhuǎn)動(dòng)線性接觸模型模擬砂土三軸壓縮試驗(yàn)；再利用正交試驗(yàn)中的方差分析法篩選出宏觀參數(shù)的控制性細(xì)觀參數(shù)；然后建立宏觀參數(shù)與控制性細(xì)觀參數(shù)的響應(yīng)面方程，依據(jù)室內(nèi)三軸試驗(yàn)確定的砂土宏觀參數(shù)，代入所建的響應(yīng)面方程求解得到控制性細(xì)觀參數(shù)；最后結(jié)合三軸試驗(yàn)結(jié)果驗(yàn)證所標(biāo)定細(xì)觀參數(shù)的準(zhǔn)確性。結(jié)果表明：宏觀參數(shù)包括彈性模量、泊松比、峰值摩擦角和剪脹角，其中彈性模量的控制性細(xì)觀參數(shù)為有效模量和摩擦系數(shù)，泊松比的控制性細(xì)觀參數(shù)為剛度比與摩擦系數(shù)，峰值摩擦角和剪脹角的控制性細(xì)觀參數(shù)都為摩擦系數(shù)；通過(guò)三軸試驗(yàn)得出宏觀參數(shù)取值，彈性模量為83.83 MPa、泊松比為0.45、峰值摩擦角為34.47°、剪脹角為8.93°；進(jìn)一步標(biāo)定細(xì)觀參數(shù)取值，有效模量為153.35 MPa，剛度比為2.16，摩擦系數(shù)為0.45；利用標(biāo)定的細(xì)觀參數(shù)進(jìn)行三軸試驗(yàn)?zāi)M，結(jié)果與室內(nèi)試驗(yàn)結(jié)果在低圍壓（100 kPa）下峰值強(qiáng)度相差7.2%，在高圍壓（300、500 kPa）下誤差都控制在了±15.0%以內(nèi)，證明了所提出的標(biāo)定方法具有準(zhǔn)確性和可靠性。

關(guān)鍵詞：砂土；離散元法；細(xì)觀參數(shù)；正交試驗(yàn)；響應(yīng)面法；宏觀參數(shù)；彈性模量；泊松比

doi：10.13278/j.cnki.jjuese.20230132

中圖分類號(hào)：TU4

文獻(xiàn)標(biāo)志碼：A

Supported by the Natural Science Foundation of Hebei Province （E2020201017），the Scientific Research Fund of Institute of Engineering Mechanics， China Earthquake Administration （2019EEEVL0202） and the Science and Technology Research Project of Higher Education Institutions in Hebei Province （ZD2020157）

Calibration of Sand Mesoscopic Parameters Based on Orthogonal-Response Surface Method

Liu Hongshuai"" Zhang Dongtao1，2

1. College of Civil Engineering and Architecture， Hebei University，Baoding 071000，Hebei，China

2. Institute of Geotechnical Engineering， Hebei University，Baoding 071000，Hebei，China

Abstract： The simulation accuracy of the particle discrete element method （DEM） mainly depends on the accuracy of calibrated mesoscopic parameters. To improve the calibration accuracy and efficiency of mesoscopic parameters， firstly，the rolling-resistance linear contact model was used to simulate sand triaxial compression tests. Then， the variance analysis method in orthogonal experiments is used to select the controlling mesoscopic parameters （CMP）， and response surface equations are established to relate the macro parameters to the CMP. By substituting the macro parameters determined from triaxial tests into the response surface equations， the values of the CMP are solved. Finally， the accuracy of the calibrated micro parameters is verified by comparing the simulation results with the triaxial test results. The results show that the macro parameters include elastic modulus， Poisson’s ratio， dilatancy angle， and peak friction angle. The CMP for the elastic modulus are effective modulus and friction coefficient， while the controlling micro parameters for Poisson’s ratio are stiffness ratio and friction coefficient. The CMP of the peak friction angle and dilatancy angle are the friction coefficient. The elastic modulus， Poisson’s ratio， peak friction angle and dilatancy angle from the triaxial tests， are 83.83 MPa， 0.45， 34.47° and 8.93°， respectively. The corresponding mesoscopic parameters are determined as follows： Effective modulus of 153.35 MPa， stiffness ratio of 2.16， and friction coefficient of 0.45. Comparing simulation results from the calibrated mesoscopic parameters with the triaxial tests， the difference in peak strength is 7.2% at a low level of confining pressure （100 kPa）， and the errors are within 15.0% at higher level of confining pressures （300 and 500 kPa）.

Key words： sand;" discrete element method; mesoscopic parameter; orthogonal test; response surface method； marco parameter； elastic modulus； Poisson’s ratio

0 引言

砂土作為淺表地層的天然散粒體或建筑材料，廣泛應(yīng)用于各種工程建設(shè)中。砂土顆粒間基本無(wú)黏結(jié)，在地震荷載、降雨和高壓差水頭等作用下，可能會(huì)發(fā)生液化、滲透變形、滑坡等非連續(xù)大變形破壞，造成大量工程破壞、人員傷亡等嚴(yán)重后果[14]。離散元法基于牛頓第二定律模擬砂土顆粒的運(yùn)動(dòng)，引入接觸模型處理顆粒單元之間的相互作用[5]，通過(guò)顆粒接觸行為，得出材料內(nèi)部受力變形的演化過(guò)程，為深入理解砂土液化、沖刷、管涌和應(yīng)變局部化的力學(xué)機(jī)理提供了有效方法。離散元的模擬精度主要取決于細(xì)觀參數(shù)標(biāo)定的準(zhǔn)確性（即細(xì)觀參數(shù)能否真實(shí)地反映砂土顆粒間的接觸力學(xué)行為）。因此，如何快速、準(zhǔn)確地標(biāo)定砂土細(xì)觀模型參數(shù)對(duì)于研究砂土力學(xué)行為及其災(zāi)變機(jī)理具有重要的理論價(jià)值和工程應(yīng)用意義。

目前，離散元細(xì)觀模型參數(shù)的標(biāo)定方法主要分為兩類：一為“試錯(cuò)法”。即通過(guò)大量試算不斷修改模型參數(shù)，待與室內(nèi)試驗(yàn)結(jié)果相匹配后結(jié)束[610]。該方法計(jì)算量大、隨機(jī)性高、耗時(shí)長(zhǎng)、標(biāo)定精度較低，且多數(shù)研究主要基于偏應(yīng)力軸向應(yīng)變關(guān)系曲線直接得到最終的細(xì)觀參數(shù)，但未明確參數(shù)標(biāo)定的詳細(xì)過(guò)程。二為“優(yōu)化法”?？梢苑譃閱我环ê蛷?fù)合法。其中單一法常用的有正交試驗(yàn)法

和PlackettBurman（PB）試驗(yàn)設(shè)計(jì)法。此類方法是

將細(xì)觀參數(shù)作為自變量，宏觀參數(shù)看作因變量，設(shè)計(jì)不同自變量在不同水平下的正交組合。先對(duì)每一組合進(jìn)行單元試驗(yàn)?zāi)M（如直接剪切試驗(yàn)、三軸壓縮試驗(yàn)等），分析各細(xì)觀參數(shù)對(duì)宏觀參數(shù)的影響規(guī)律；再通過(guò)比較各細(xì)觀參數(shù)對(duì)宏觀參數(shù)影響的相對(duì)大小，篩選出各宏觀參數(shù)的控制性細(xì)觀參數(shù)，從而確定調(diào)整哪些細(xì)觀參數(shù)能快速擬合室內(nèi)試驗(yàn)結(jié)果[1115]。雖然正交試驗(yàn)法和PB設(shè)計(jì)法明確了細(xì)觀參數(shù)對(duì)宏觀參數(shù)的影響規(guī)律，在很大程度上縮小了試算范圍，但是并不能直接確定控制性細(xì)觀參數(shù)的取值，仍需縮小范圍進(jìn)一步試算。為此，部分學(xué)者在單一法的基礎(chǔ)上，結(jié)合等值線法[16]、隨機(jī)森林法[17]等提出了復(fù)合優(yōu)化法。

相較于正交試驗(yàn)法和PB設(shè)計(jì)法，復(fù)合優(yōu)化法在明確宏細(xì)觀參數(shù)之間的敏感性后不需繼續(xù)試錯(cuò)，而是通過(guò)建立數(shù)學(xué)關(guān)系式，直接求出細(xì)觀參數(shù)的具體值，這不僅降低了計(jì)算成本，且模擬精度也可以得到保障。

本文基于離散元軟件PFC3D，選取抗轉(zhuǎn)動(dòng)線性接觸模型，采用正交試驗(yàn)確定控制性細(xì)觀參數(shù)，運(yùn)用響應(yīng)面法構(gòu)建各宏觀參數(shù)與控制性細(xì)觀參數(shù)的預(yù)測(cè)方程，提出一種新的砂土細(xì)觀參數(shù)標(biāo)定方法——正交響應(yīng)面復(fù)合優(yōu)化法，并利用三軸試驗(yàn)結(jié)果驗(yàn)證所提出方法的有效性。

1 砂土離散元建模及試驗(yàn)指標(biāo)

1.1 接觸模型

砂土離散元模型若考慮顆粒不規(guī)則形狀，則需用顆粒簇來(lái)模擬，會(huì)造成巨大的計(jì)算代價(jià)。已有研究結(jié)果[1819]表明，只要能合理考慮顆粒間的轉(zhuǎn)動(dòng)阻力，球形顆粒亦能很好地模擬不規(guī)則顆粒的力學(xué)行為；因此，本文選取抗轉(zhuǎn)動(dòng)線性接觸模型模擬砂土的力學(xué)行為。

抗轉(zhuǎn)動(dòng)線性接觸模型包括法向、切向和抗轉(zhuǎn)動(dòng)三種接觸類型（圖1），其共同點(diǎn)是都包含彈簧及與之并聯(lián)的阻尼器，彈簧和阻尼器分別反映顆粒接觸面的力學(xué)行為和能量耗散特性。

法向接觸同時(shí)還包含分離器，用于模擬法向只受壓不受拉的力學(xué)行為。通常假定法向壓力與顆粒間的法相位移重疊量成正比，若忽略阻尼力，法向力Fn計(jì)算公式為：

Fn=KnUn，Un≥0；Fn=0，Unlt;0。（1）

式中：Kn為法向剛度；Un為顆粒間法向重疊量。

切向接觸與法向接觸不同的是滑片代替了法向接觸中的分離器，用于模擬接觸點(diǎn)超過(guò)剪切強(qiáng)度后的塑性滑動(dòng)。切向力Fs的計(jì)算公式為：

Fs，max=μFn； Fs=Fs0+KsΔUs。（2）

式中：Fs，max為切向力閾值，超過(guò)該值后切向力不再增大；μ為摩擦系數(shù)；Fs0為計(jì)算開始前的切向力；Ks為切向剛度；ΔUs為顆粒間切向位移重疊量的變化量。

抗轉(zhuǎn)動(dòng)接觸模型由抗轉(zhuǎn)動(dòng)彈簧和阻尼器組成，用于等效彌補(bǔ)圓形顆粒不能模擬真實(shí)顆粒的抗轉(zhuǎn)動(dòng)能力?？罐D(zhuǎn)動(dòng)力矩Mr的計(jì)算公式為：

Mr，max=μrR-Fn；Mr=Mr0+KrΔθr；1R-=1RA+1RB。（3）

式中：Mr，max為抗轉(zhuǎn)動(dòng)力矩閾值，超過(guò)該值后抗轉(zhuǎn)動(dòng)力矩不再增大；μr為轉(zhuǎn)動(dòng)阻力系數(shù)；R-為等效半徑；Mr0為計(jì)算開始前的抗轉(zhuǎn)動(dòng)力矩；Kr為扭轉(zhuǎn)剛度，與接觸點(diǎn)的切向剛度有關(guān)；Δθr為相對(duì)轉(zhuǎn)角的變化量；RA和RB分別為相接觸兩個(gè)球的半徑。

1.2 離散元模擬

離散元模擬分為數(shù)值制樣、固結(jié)、加載三個(gè)階段。1）數(shù)值制樣：以砂土室內(nèi)三軸壓縮試驗(yàn)為研究對(duì)象，試樣高100 mm，直徑49.9 mm，所使用砂土的級(jí)配如圖2所示（粒徑范圍為0.1～ 0.5 mm）。考慮到若按實(shí)際粒徑大小生成試驗(yàn)顆粒，數(shù)目過(guò)于龐大，計(jì)算代價(jià)難以接受；為減小計(jì)算量，本文采用與土工試驗(yàn)粒徑相當(dāng)、縮小數(shù)值試樣尺寸的方法。根據(jù)相關(guān)研究成果，只要滿足試驗(yàn)的容器直徑（本文為49.9 mm）與最大粒徑的比值大于5[2021]，縮小試樣直徑或擴(kuò)大顆粒直徑所得到的結(jié)果即能滿足精度要求。本文試驗(yàn)顆粒的最大粒徑為0.5 mm，故將土樣尺寸縮小為高15 mm、直徑7.5 mm，可滿足室內(nèi)三軸試驗(yàn)高徑比為2的條件。采用剛性墻模擬三軸橡皮膜，上、下底部各生成一個(gè)無(wú)限大的正方形墻體模擬加載板。按照試驗(yàn)級(jí)配生成9 913個(gè)顆粒，顆粒密度設(shè)置為2 650 kg/m3，細(xì)觀接觸模型為抗轉(zhuǎn)動(dòng)

線性模型。在生成試樣過(guò)程中，有效模量和剛度比取自正交試驗(yàn)設(shè)計(jì)方案，先將顆粒間、顆粒與墻體的摩擦系數(shù)和抗轉(zhuǎn)動(dòng)系數(shù)均設(shè)為0，以便生成均勻的數(shù)值試樣；之后控制墻體應(yīng)力為10 kPa，改變顆粒間的摩擦系數(shù)和抗轉(zhuǎn)動(dòng)系數(shù)，用于生成不同密實(shí)度的試樣。2）固結(jié)：設(shè)置固結(jié)應(yīng)力分別為100、300、500 kPa進(jìn)行雙向固結(jié)，待達(dá)到指定圍壓后關(guān)閉上下墻伺服。3）加載：對(duì)固結(jié)完成后的試樣進(jìn)行加載，按5％的軸向應(yīng)變率加載，待軸向應(yīng)變達(dá)到20%后完成加載。在加載過(guò)程中記錄試驗(yàn)的偏應(yīng)力、軸向應(yīng)變和體積應(yīng)變。為模擬理想的三軸應(yīng)力狀態(tài)，整個(gè)數(shù)值試驗(yàn)忽略了墻體和顆粒間的摩擦力和顆粒重力。

1.3 宏觀力學(xué)參數(shù)確定方法

本文選擇工程應(yīng)用最為廣泛的彈性模量E50、泊松比υ、剪脹角和峰值摩擦角θ為試驗(yàn)宏觀參數(shù)。相關(guān)參數(shù)的選取規(guī)則如圖3所示。

各宏觀力學(xué)參數(shù)計(jì)算公式如下。

1）彈性模量：E50=σ50εa，50。（4）

2）泊松比：υ=-εr，50εa，50。（5）

3）剪脹角：=arcsinεp，v-2εp，a+εp，v。（6）

4）峰值摩擦角：θ=arcsin3Mp6+Mp。（7）

其中：

Mp=qfpf[22]；（8）

qf=σ1f—σ3f；（9）

pf=σ1f+2σ3f3。（10）

式中：σ50為峰值強(qiáng)度取一半時(shí)對(duì)應(yīng)的偏應(yīng)力；

εa，50為彈性模量所對(duì)應(yīng)的軸向應(yīng)變；

εr，50為彈性模量所對(duì)應(yīng)的徑向應(yīng)變；

εp，v為峰值強(qiáng)度對(duì)應(yīng)的體積應(yīng)變；

εp，a為峰值強(qiáng)度對(duì)應(yīng)的軸向應(yīng)變；

Mp為峰值強(qiáng)度比；σ1f和σ3f分別為極限狀態(tài)下第一和第三主應(yīng)力。

2 細(xì)觀參數(shù)正交試驗(yàn)及其方差分析

每個(gè)細(xì)觀參數(shù)都會(huì)對(duì)宏觀參數(shù)有不同程度的影響，同時(shí)，各細(xì)觀參數(shù)之間的交互效應(yīng)對(duì)宏觀力學(xué)參數(shù)的影響也不可忽略。為此，本節(jié)采用正交試驗(yàn)中的方差分析法對(duì)細(xì)觀參數(shù)進(jìn)行敏感性分析，給出各宏觀參數(shù)的控制性細(xì)觀參數(shù)，為建立宏觀參數(shù)和細(xì)觀參數(shù)的非線性關(guān)系提供依據(jù)。

2.1 細(xì)觀參數(shù)取值范圍初步確定

本文采用的細(xì)觀本構(gòu)模型為抗轉(zhuǎn)動(dòng)線性接觸模型，主要包括有效模量、摩擦系數(shù)、抗轉(zhuǎn)動(dòng)系數(shù)、剛度比四個(gè)細(xì)觀參數(shù)。為減小搜索范圍，本論文收集分析了上述細(xì)觀參數(shù)的經(jīng)驗(yàn)范圍值[2324]：有效模量一般在64～300 MPa之間，在滿足彈性球體剛度模型的基礎(chǔ)上剛度比應(yīng)該在1.0～1.5之間。參照已有的研究結(jié)果，本次正交試驗(yàn)設(shè)計(jì)取值范圍見表1。

2.2 正交試驗(yàn)方案及結(jié)果

本節(jié)在前述基礎(chǔ)上設(shè)計(jì)了四因素五水平的正交試驗(yàn)，共25組，具體試驗(yàn)方案如表1所示。四因素表示不同的細(xì)觀參數(shù)，五水平表示細(xì)觀參數(shù)的不同取值大小。根據(jù)正交試驗(yàn)設(shè)計(jì)進(jìn)行三軸壓縮試驗(yàn)?zāi)M，結(jié)合所給出的宏觀參數(shù)確定方法，計(jì)算出四因素五水平所對(duì)應(yīng)的彈性模量、泊松比、峰值摩擦角和剪脹角。正交試驗(yàn)匯總結(jié)果列于表2，可得出各宏觀參數(shù)范圍：彈性模量為24～281 MPa；泊松比為0.35～0.53；峰值摩擦角為25.3°～39.2°；剪脹角為5.3°～16.3°。

2.3 方差分析

方差分析能夠同時(shí)分析多個(gè)細(xì)觀參數(shù)對(duì)宏觀參數(shù)影響的顯著性，并能按照影響的重要程度進(jìn)行排序，從而確定控制宏觀力學(xué)參數(shù)的細(xì)觀參數(shù)（控制性細(xì)觀參數(shù)）。F統(tǒng)計(jì)量是方差分析法的重要組成部分，表示組間方差與組內(nèi)方差的比值。 F統(tǒng)計(jì)量通常用伴隨概率p來(lái)表示，本文取檢驗(yàn)顯著水平因子α=0.05，當(dāng)plt;α?xí)r認(rèn)為細(xì)觀參數(shù)對(duì)宏觀參數(shù)影響顯著（p值越小表示越顯著），否則認(rèn)為影響不顯著，可忽略不計(jì)。表2中各細(xì)觀參數(shù)對(duì)宏觀參數(shù)影響的正交試驗(yàn)方差統(tǒng)計(jì)結(jié)果見圖4。分析圖4可以得出：彈性模量的控制性細(xì)觀參數(shù)分別為有效模量和摩擦系數(shù)，其余兩個(gè)參數(shù)的影響可忽略不計(jì)；泊松比的控制性細(xì)觀參數(shù)分別為摩擦系數(shù)和剛度比，其余兩個(gè)細(xì)觀參數(shù)可忽略不計(jì)；峰值摩擦角和剪脹角的控制性細(xì)觀參數(shù)只有摩擦系數(shù)，其余三個(gè)參數(shù)可忽略不計(jì)。

3 響應(yīng)面法

響應(yīng)面法是一種多元統(tǒng)計(jì)分析方法，為用于建立因變量和自變量之間非線性關(guān)系的數(shù)學(xué)模型。因此本節(jié)采用響應(yīng)面法建立宏觀參數(shù)和控制性細(xì)觀參數(shù)間的預(yù)測(cè)方程。

運(yùn)用Design-Expert軟件，對(duì)控制性細(xì)觀參數(shù)和宏觀參數(shù)進(jìn)行中心復(fù)合設(shè)計(jì)（CCD）試驗(yàn)，采用二階模型來(lái)近似表示響應(yīng)。為便于表示，將細(xì)觀參數(shù)有效模量、摩擦系數(shù)、抗轉(zhuǎn)動(dòng)系數(shù)和剛度比分別表示為x1、x2、x3和x4，宏觀參數(shù)彈性模量、泊松比、峰值摩擦角和剪脹角分別表示為y1、y2、y3和y4。

3.1 彈性模量響應(yīng)曲面

由前述的方差分析結(jié)果可知，抗轉(zhuǎn)動(dòng)系數(shù)和剛度比對(duì)彈性模量的影響可忽略，因此將這兩個(gè)參數(shù)取值分別設(shè)為0.5和1.5。設(shè)計(jì)13組CCD試驗(yàn)，計(jì)算出彈性模量對(duì)控制性細(xì)觀參數(shù)（有效模量和摩擦系數(shù)）的響應(yīng)值，根據(jù)試驗(yàn)設(shè)計(jì)模擬三軸壓縮試驗(yàn)，結(jié)果列于表3。采用二次回歸法擬合得到宏觀參數(shù)彈性模量與細(xì)觀參數(shù)有效模量、摩擦系數(shù)間的響應(yīng)面方程：

y1=109.28742+0.61084x1-338.21678x2-0.22381x1x2-0.0001624x21+258.44835x22。

（11）

運(yùn)用表3的數(shù)據(jù)對(duì)式（11）進(jìn)行檢驗(yàn)，所得誤差均在±10%以內(nèi)，可見響應(yīng)面方程具有較高的精度，可滿足工程要求。

依據(jù)式（11）計(jì)算得到了宏觀參數(shù)彈性模量與其控制性細(xì)觀參數(shù)有效模量、摩擦系數(shù)間的響應(yīng)曲面。圖5a為彈性模量的等值線圖，可以看出當(dāng)有效模量一定時(shí)，隨著摩擦系數(shù)增大，彈性模量逐漸減??；當(dāng)摩擦系數(shù)一定時(shí)，彈性模量隨有效模量的增大而增大。

由圖5b可知：彈性模量和有效模量存在較強(qiáng)的正相關(guān)性，當(dāng)有效模量增大時(shí)，彈性模量也隨著增大；而彈性模量有隨摩擦系數(shù)增大而減小的趨勢(shì)，尤其是在摩擦系數(shù)小于0.50的情況下，彈性模量的下降速度要比摩擦系數(shù)大于0.50時(shí)更大，其對(duì)彈性模量的影響比大于0.50時(shí)更為顯著；此外，有效模量方向的變化梯度要明顯大于摩擦系數(shù)方向，這表明彈性模量的響應(yīng)值對(duì)有效模量更加敏感。該特征在方差分析中的p值大小上也有體現(xiàn)。

3.2 泊松比響應(yīng)曲面

泊松比的控制性細(xì)觀參數(shù)為摩擦系數(shù)和剛度比，非控制性細(xì)觀參數(shù)有效模量取150 MPa，抗轉(zhuǎn)動(dòng)系數(shù)取0.5，設(shè)計(jì)13組CCD試驗(yàn)計(jì)算泊松比對(duì)剛度比和摩擦系數(shù)的響應(yīng)值，模擬三軸壓縮試驗(yàn)，結(jié)果列于表4。

采用二次回歸法擬合得到宏觀泊松比與細(xì)觀剛度比、摩擦系數(shù)間的響應(yīng)面方程：

y2=0.2593+0.17943x2+0.078588x4+0.030833x2x4-0.15404x22-0.012173x24。

（12）

用表4的數(shù)據(jù)對(duì)式（12）進(jìn)行檢驗(yàn)，所得誤差均在±10%以內(nèi)，可見響應(yīng)面方程具有較高的精度，滿足工程要求。

依據(jù)式（12）的響應(yīng)面方程計(jì)算可得宏觀參數(shù)泊松比與細(xì)觀參數(shù)摩擦系數(shù)、剛度比的響應(yīng)曲面，見圖6。

圖6a為泊松比的等值線圖，可見：當(dāng)剛度比一定時(shí)，隨著摩擦系數(shù)增大，泊松比逐漸增大；當(dāng)摩擦系數(shù)一定時(shí)，泊松比隨剛度比的增大而增大。

如圖6b所示：泊松比與摩擦系數(shù)和剛度比都呈現(xiàn)較強(qiáng)的正相關(guān)性，當(dāng)摩擦系數(shù)小于0.50時(shí)，曲線更加陡峭，表明摩擦系數(shù)較小時(shí)對(duì)泊松比的影響較大；剛度比也有類似的規(guī)律，當(dāng)剛度比小于2.00時(shí)，泊松比的變化梯度比剛度比大于2.00時(shí)更大；此外剛度比方向上的變化梯度要明顯大于摩擦系數(shù)方向的，說(shuō)明泊松比對(duì)剛度比的變化更加敏感。

3.3 峰值摩擦角、剪脹角與摩擦系數(shù)擬合

由于影響峰值摩擦角和剪脹角的控制性細(xì)觀參數(shù)只有摩擦系數(shù)，因此只需對(duì)摩擦系數(shù)取不同值進(jìn)行重復(fù)試驗(yàn)，得出結(jié)果后進(jìn)行多項(xiàng)式擬合即可。將摩擦系數(shù)設(shè)為0.20、0.40、0.60、0.08和1.00。進(jìn)行三軸壓縮試驗(yàn)?zāi)M得到峰值摩擦角分別為29.82°、33.93°、36.17°，37.41°和38.15°，剪脹角分別為7.57°、8.94°、9.06°、8.07°和7.20°。通過(guò)二次曲線擬合得到響應(yīng)面方程：

y3=25.19576+26.8282x2-13.80355x22;

（13）

y4=5.856+11.17357x2-9.98214x22。

（14）

依據(jù)式（13）與式（14）響應(yīng)面曲線方程分別計(jì)算得到峰值摩擦角、剪脹角與摩擦系數(shù)的擬合曲線圖（圖7、8）。由圖7可見，峰值摩擦角與摩擦系數(shù)呈正相關(guān)性，且摩擦系數(shù)在0.60以下時(shí)斜率更大，可見在摩擦系數(shù)較小時(shí)，峰值摩擦角對(duì)摩擦系數(shù)更加敏感。由圖8可知，剪脹角隨摩擦系數(shù)的增大呈現(xiàn)

先增大后減小的趨勢(shì)，在摩擦系數(shù)為0.55時(shí)剪脹角

達(dá)到峰值。這可能是由于：剪脹角是描述試樣剪脹規(guī)律的參數(shù)，在較小的摩擦系數(shù)下，顆粒間的摩擦阻力較小，試樣受到剪切時(shí)更容易發(fā)生滑移，導(dǎo)致剪脹現(xiàn)象增強(qiáng)，剪脹角增大；隨著摩擦系數(shù)的增大，顆粒間摩擦阻力增加，導(dǎo)致剪脹現(xiàn)象減弱，剪脹角增長(zhǎng)速度減??；當(dāng)摩擦系數(shù)進(jìn)一步增大時(shí)，顆粒間摩擦阻力達(dá)到一定程度，導(dǎo)致顆粒緊密堆積且發(fā)生互鎖現(xiàn)象使土體體積變化受到限制，剪脹角開始減小。

4 參數(shù)標(biāo)定及試驗(yàn)驗(yàn)證

首先依據(jù)室內(nèi)試驗(yàn)獲取宏觀力學(xué)參數(shù)，將其代入響應(yīng)面預(yù)測(cè)方程，求解得到控制性細(xì)觀參數(shù)；再結(jié)合已有數(shù)值模擬經(jīng)驗(yàn)，綜合確定控制性細(xì)觀參數(shù)取值；然后根據(jù)標(biāo)定出的細(xì)觀參數(shù)進(jìn)行三軸壓縮試驗(yàn)?zāi)M，得出宏觀應(yīng)力應(yīng)變曲線，并與試驗(yàn)曲線進(jìn)行對(duì)比分析；最后完成細(xì)觀參數(shù)標(biāo)定的驗(yàn)證工作。

將圍壓設(shè)為100 kPa對(duì)福建標(biāo)準(zhǔn)砂進(jìn)行室內(nèi)三軸試驗(yàn)，所選福建砂的物理指標(biāo)如表5所示，圖9為室內(nèi)三軸試驗(yàn)相應(yīng)的偏應(yīng)力軸向應(yīng)變和體積應(yīng)變軸向應(yīng)變曲線。

依據(jù)圖9可分別確定彈性模量、泊松比、峰值摩擦角和剪脹角，結(jié)果如表6所示。細(xì)觀參數(shù)中，摩擦系數(shù)作為單因素變量分別與峰值摩擦角（式（13））、剪脹角（式（14））成對(duì)應(yīng)關(guān)系，代入相關(guān)公式計(jì)算得到的摩擦系數(shù)分別為0.45、0.63，二者差別顯著。考慮到砂土為摩擦材料，因此本次標(biāo)定所使用的摩擦系數(shù)取0.45。由式（11）和式（12）可知，摩擦系數(shù)確定后，彈性模量與有效模量、泊松比與剛度比成對(duì)應(yīng)關(guān)系。綜合表6的彈性模量、泊松比和式（11）（12），便可得到有效模量和剛度比。由前述可知，抗轉(zhuǎn)動(dòng)系數(shù)為所有宏觀參數(shù)的非控制性細(xì)觀參數(shù)，本文采用Gu等[24]建議值，抗轉(zhuǎn)動(dòng)系數(shù)取為0.5。細(xì)觀參數(shù)標(biāo)定的最終結(jié)果列于表7。

根據(jù)標(biāo)定的細(xì)觀參數(shù)建立100、300、500 kPa三種圍壓條件下的三軸壓縮試驗(yàn)?zāi)M，給出應(yīng)力應(yīng)變曲線，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比，以此來(lái)驗(yàn)證所標(biāo)定細(xì)觀參數(shù)的可靠性。對(duì)比模擬結(jié)果和試驗(yàn)結(jié)果（圖10）可看出，兩者的總體趨勢(shì)基本一致，能夠較合理地預(yù)測(cè)應(yīng)變軟化現(xiàn)象。與標(biāo)定圍壓相一致的情況下（100 kPa），擬合效果最好，預(yù)測(cè)峰值強(qiáng)度比試驗(yàn)值低7.2%；隨圍壓的增大，即距標(biāo)定圍壓的差別越大，預(yù)測(cè)峰值強(qiáng)度的誤差越大（圍壓為300 kPa時(shí)，峰值預(yù)測(cè)誤差12.1%；圍壓為500 kPa時(shí)，峰值預(yù)測(cè)誤差13.4%），但峰值預(yù)測(cè)誤差均小于15.0%。由此可知，當(dāng)使用所標(biāo)定出的細(xì)觀參數(shù)進(jìn)行不同于標(biāo)定圍壓的試驗(yàn)時(shí)，其預(yù)測(cè)的強(qiáng)度峰值誤差會(huì)變大。因此若模擬圍壓變化較大時(shí)，為提高模擬精度，建議采用模擬圍壓上下限的試驗(yàn)數(shù)據(jù)聯(lián)合標(biāo)定來(lái)保證結(jié)果的可靠性。

5 結(jié)論與建議

1）彈性模量的控制性細(xì)觀參數(shù)為有效模量和摩擦系數(shù)，泊松比的控制性細(xì)觀參數(shù)為摩擦系數(shù)和剛度比，峰值摩擦角和剪脹角的控制性細(xì)觀參數(shù)為摩擦系數(shù)。

2）所提出的正交響應(yīng)面細(xì)觀參數(shù)標(biāo)定法所標(biāo)定出的細(xì)觀參數(shù)能夠較合理地反映砂土的應(yīng)力應(yīng)變特性。與標(biāo)定圍壓相一致的應(yīng)力情況下 （100 kPa），預(yù)測(cè)峰值強(qiáng)度比試驗(yàn)值低7.2%。

3）當(dāng)使用所標(biāo)定出的細(xì)觀參數(shù)進(jìn)行不同于標(biāo)定圍壓的試驗(yàn)時(shí)，其預(yù)測(cè)的峰值強(qiáng)度誤差會(huì)變大，圍壓為300 kPa時(shí)，峰值預(yù)測(cè)誤差12.1%；圍壓為500 kPa時(shí)，峰值預(yù)測(cè)誤差13.4%。

4）若模擬試驗(yàn)圍壓變化較大時(shí)，為提高模擬精度，建議采用模擬圍壓上下限的試驗(yàn)數(shù)據(jù)聯(lián)合標(biāo)定來(lái)保證結(jié)果的可靠性。

參考文獻(xiàn)（References）：

[1] He S Y，Wang， X L， Fan H B， et al. The Study on Loess Liquefaction in China： A Systematic Review[J]. Natural Hazards， 2020， 103（2）： 16391669.

[2] 蘇德辰，楊戰(zhàn)兵，孫愛萍，等. 山西云崗石窟侏羅系地震液化砂巖柱的發(fā)現(xiàn)及其大地構(gòu)造意義[J]. 地質(zhì)學(xué)報(bào)，2019， 93（8）： 18141830.

Su Dechen， Yang Zhanbing， Sun Aiping， et al. Discovery of the Jurassic Seismic Liquefaction Sandstone Pipes in the Yungang Grottoes， Shanxi Province and Its Tectonic Significance[J]. Acta Geologica Sinica， 2019， 93（8）：" 18141830.

[3] 劉紅帥，宋東松，張東濤，等. 干砂小應(yīng)變動(dòng)力學(xué)參數(shù)的彎曲元試驗(yàn)[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2023，53（6）：18261834.

Liu Hongshuai， Song Dongsong， Zhang Dongtao， et al. Small Strain Dynamic Parameters of Dry Sand by Bender Element[J]. Journal of Jilin University （Earth Science Edition）， 2023，53（6）：18261834.

[4] 劉鑫，苗雪青，黃良，等. 5種排水條件下飽和重塑黃土三軸剪切特性[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版）， 2023， 53（5）：14991509.

Liu Xin，Miao Xueqing，Huang Liang，et al. Triaxial Shear Behavior of Saturated Remolded Loess Subjected to Five Drainage Conditions[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（5）：14991509.

[5] Cundall P A， Strack O D L. A Discrete Numerical Model for Granular Assemblies[J]. Geotechnic， 1979， 29（1）： 4765.

[6] 鄧志穎，范祥，何忠明，等. L形邊裂隙單圓孔砂巖單軸壓縮力學(xué)特性[J]. 吉林大學(xué)學(xué)報(bào)（地球科學(xué)版）， 2023， 53（2）： 541554.

Deng Zhiying， Fan Xiang， He Zhongming， et al. Uniaxial Compressive Mechanical Properties of Sandstone with Single Circular Hole and L-Shaped Side Fissures [J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（2）： 541554.

[7] 孫佳政，傅翼，仇雅詩(shī)，等. RB模式下砂土非極限主動(dòng)土壓力的離散元模擬與理論研究[J]. 巖土力學(xué)，2023， 44（2）： 603614.

Sun Jiazheng， Fu Yi， Qiu Yashi， et al. Discrete Element Simulation and Theoretical Study on Non-Limit Active Earth Pressure of Sand Under RB Mode[J]. Rock and Soil Mechanics， 2023， 44（2）：" 603614.

[8] 葉斌，宋思聰，倪雪倩. 制樣方法對(duì)砂土液化力學(xué)性質(zhì)影響的離散元模擬[J]. 同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版），2022， 50（7）： 9981008.

Ye Bin， Song Sicong， Ni Xueqian. Effects of Sample Preparation Method on Sand Liquefaction Mechanical Properties with Discrete Element Simulation[J]. Journal of Tongji University （Natural Science）， 2022， 50（7）： 9981008.

[9] 楊文劍，劉懷忠，謝紅強(qiáng)，等. 基于直剪試驗(yàn)的堆積碎石類材料細(xì)觀參數(shù)標(biāo)定方法及模擬驗(yàn)證[J]. 工程科學(xué)與技術(shù)，2022， 54（3）： 4654.

Yang Wenjian， Liu Huaizhong， Xie Hongqiang， et al. Mesoscopic Parameter Calibration Method of Accumulated Debris Materials Based on Direct Shear Test and Simulation Verification[J]. Advanced Engineering Sciences， 2022， 54（3）： 4654.

[10] 蔣明鏡，石安寧，劉俊，等. 結(jié)構(gòu)性砂土力學(xué)特性三維離散元分析[J]. 巖土工程學(xué)報(bào)，2019， 41（增刊2）： 14.

Jiang Mingjing， Shi Anning， Liu Jun， et al. Three-Dimensional Distinct Element Analysis of Mechanical Properties of Structured Sands[J]. Chinese Journal of Geotechnical Engineering， 2019， 41（Sup.2）： 14.

[11] 許江波，曹寶花，余洋林，等. 基于PFC3D的黃土三軸試驗(yàn)細(xì)觀參數(shù)敏感性分析[J]. 工程地質(zhì)學(xué)報(bào)，2021， 29（5）： 13421353.

Xu Jiangbo， Cao Baohua， Yu Yanglin， et al. Sensitivity Analysis of Meso-Parameters in Loess Triaxial Test Based on PFC3D[J]. Journal of Engineering Geology， 2021， 29（5）：" 13421353.

[12] 張寶玉，張昌鎖，王晨龍，等. PFC2D平節(jié)理模型細(xì)觀參數(shù)標(biāo)定方法[J]. 計(jì)算力學(xué)學(xué)報(bào)，2021， 38（5）： 665673.

Zhang Baoyu， Zhang Changsuo， Wang Chenlong， et al. Calibration Method of Meso-Parameters of PFC2D Flat-Joint Model[J]. Chinese Journal of Computational Mechanics， 2021， 38（5）： 665673.

[13] 陳坤，廖良健，朱瑞，等. 煤系土離散元模型宏細(xì)觀參數(shù)標(biāo)定及其敏感性分析[J]. 南昌工程學(xué)院學(xué)報(bào)，2023， 42（1）： 3944.

Chen Kun， Liao Liangjian， Zhu Rui， et al. Calibration and Sensitivity Analysis of Macro and Meso Parameters of Discrete Element Model for Coal Measure Soil[J]. Journal of Nanchang Institute of Technology， 2023， 42（1）： 3944.

[14] 崔熙燦，張凌凱，王建祥. 高堆石壩砂礫石料的細(xì)觀參數(shù)反演及三軸試驗(yàn)?zāi)M[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2022， 38（4）： 113122.

Cui Xican， Zhang Lingkai， Wang Jianxiang， et al. Inversion of Meso Parameters and Triaxial Test Simulation of the Gravel Materials for High Rockfill Dam[J]. Transactions of the Chinese Society of Agricultural Engineering， 2022， 38（4）： 113122.

[15]

汪洋，孟祥林，高進(jìn)，等. 基于PFC2d的平行粘結(jié)模型宏細(xì)觀參數(shù)標(biāo)定研究[J]. 四川建筑，2022， 42（6）： 137140.

Wang Yang， Meng Xianglin， Gao Jin， et al. Research on Calibration of Macro and Micro Parameters of Linear Parallel Bond Model Based on PFC2d[J]. Sichuan Architecture，2022， 42（6）： 137140.

[16] 王晉偉，遲世春，邵曉泉，等. 正交–等值線法在堆石料細(xì)觀參數(shù)標(biāo)定中的應(yīng)用[J]. 巖土工程學(xué)報(bào)，2020， 42（10）： 18671875.

Wang Jinwei， Chi Shichun， Shao Xiaoquan， et al. Application of Orthogonal-Contour Method in Calibration of Microscopic Parameters of Rockfill Materials[J]. Chinese Journal of Geotechnical Engineering， 2020， 42（10）： 18671875.

[17] 馬少坤，黃驍，韋榕寬，等. 考慮滾動(dòng)阻抗的線性接觸模型離散元宏細(xì)觀參數(shù)敏感性研究[J]. 中國(guó)安全生產(chǎn)科學(xué)技術(shù)，2021， 17（6）： 104110.

Ma Shaokun， Huang Xiao， Wei Rongkuan， et al. Study on Sensitivity of Macro and Micro Parameters of Discrete Element in Linear Contact Model Considering Rolling Resistance[J]. Journal of Safety Science and Technology， 2021， 17（6）： 104110.

[18] Wensrich C M， Katterfeld A. Rolling Friction as a Technique for Modelling Particle Shape in DEM[J]. Powder Technology， 2012，217： 409417.

[19] Jiang M J， Shen Z F， Wang J F. A Novel Three-Dimensional Contact Model for Granulates Incorporating Rolling and Twisting Resistances[J]. Computers and Geotechnics， 2015， 65： 147163.

[20] 劉海濤，程曉輝. 粗粒土尺寸效應(yīng)的離散元分析[J]. 巖土力學(xué)，2009， 30（增刊1）： 287292.

Liu Haitao， Cheng Xiaohui. Discrete Element Analysis for Size Effects of Coarse-Grained Soils[J]. Rock and Soil Mechanics， 2009， 30（Sup.1）： 287292.

[21] 尹小濤，鄭亞娜，馬雙科. 基于顆粒流數(shù)值試驗(yàn)的巖土材料內(nèi)尺度比研究[J]. 巖土力學(xué)，2011， 32（4）： 12111215.

Yin Xiaotao， Zheng Yana， Ma Shuangke. Study of Inner Scale Ratio of Rock and Soil Material Based on Numerical Test of Particle Flow Code[J]. Rock and Soil Mechanics， 2011， 32（4）： 12111215.

[22] 童森杰，黃茂松，時(shí)振昊，等. 基于初始狀態(tài)參數(shù)的砂土峰值摩擦角和極限強(qiáng)度[J]. 巖石力學(xué)與工程學(xué)報(bào)，2022， 41（2）： 389398.

Tong Senjie， Huang Maosong， Shi Zhenhao， et al. Peak Frictional Angle and Ultimate Drainage Strength of Sand Based on Initial State Parameter[J]. Chinese Journal of Rock Mechanics and Engineering， 2022， 41（2）： 389398.

[23] Goldenberg C， Goldhirsch I. Friction Enhances Elasticity in Granular Solids[J]. Nature， 2005， 435： 188191.

[24] Gu X， Zhang J， Huang X. DEM Analysis of Monotonic and Cyclic Behaviors of Sand Based on Critical State Soil Mechanics Framework[J]. Computers and Geotechnics， 2020， 128： 103787.