基于统一状态参数模型的砂中柱孔扩张分析.pdf

DOI：10.16030/ki.issn.1000-3665.202210044刘贯飞，雷胜友.基于统一状态参数模型的砂中柱孔扩张分析 J.水文地质工程地质，2023，50(6):99-111.LIU Guanfei,LEI Shengyou.An analysis of cylindrical cavity expansion in sand based on a unified state parameter modelJ.Hydrogeology&Engineering Geology,2023,50(6):99-111.基于统一状态参数模型的砂中柱孔扩张分析刘贯飞1,2，雷胜友1（1.长安大学公路学院，陕西 西安710064；2.中国十九冶集团有限公司，四川 成都610031）摘要：目前关于不同初始状态砂土砂中柱孔扩张的研究结果还缺乏更深层次的分析，并且由于没有考虑砂土屈服面形状因素的影响，造成许多研究成果难以在不同类型砂土中推广。采用统一状态参数模型（clay and sand model，CASM）和Rowe 剪胀方程来描述砂的弹塑性变形特点，结合大变形理论并引入辅助变量，推导了基于拉格朗日描述的弹塑性区内砂土体积和有效应力的一阶偏微分方程组，在此基础上结合弹塑性区的边界条件和柱孔扩张弹性解，建立了饱和砂中的排水柱孔扩张半解析解。结果表明，CASM 可以通过改变应力状态参数 n 和间距比 r*的值使砂的屈服面形状发生改变，进而使文中解答能够用于不同类型饱和砂中的排水柱孔扩张计算，其中 n、r*值越大，松砂初次屈服时的偏应力和后续砂中的扩孔压力越大，但中密、密实砂土中的情况与松砂完全相反。极限扩孔压力随砂土初始状态参数的减小而增大，相应的砂土体积也从一直剪缩变为先剪胀后剪缩，弹塑性区半径先减小后增大，硬化行为从一直硬化变为先软化后硬化。静止侧压力系数增大时，极限扩孔压力也增大，但对砂的体积变化规律影响不大。本研究可为相关岩土工程问题分析提供可靠理论支持。关键词：柱孔扩张；砂土；统一状态参数模型；屈服面形状；拉格朗日描述中图分类号：TU473 文献标志码：A 文章编号：1000-3665（2023）06-0099-13An analysis of cylindrical cavity expansion in sand based on aunified state parameter modelLIU Guanfei1,2，LEI Shengyou1（1.School of Highway,Changan University,Xian,Shaanxi710064,China；2.China 19th Metallurgical Corporation,Chengdu,Sichuan610031,China）Abstract：The current studies of the cylindrical cavity expansion in sand still lack a deeper analysis for the resultsobtained in sands with different initial states.In addition,most of previous studies have not considered theinfluence of the shape of yield surface of sand,which made their results difficult to be popularized in differenttypes of sand.In order to obtain a general solution to expansion of cylindrical cavity in sand under the drainedcondition,a unified state parameter model-clay and sand model(CASM)with Rowes stressdilatancy relation isused to describe the characteristics of elastic-plastic deformation of sand.By employing that large strain occurs insand and introducing an auxiliary variable,several partial differential equations to calculate the effective stress andspecific volume of sand in the elastoplastic zone are derived on the basis of the Lagrangian description.Under theelastic-plastic boundary conditions and the elastic solution of cylindrical cavity expansion,a semi-analytical 收稿日期：2022-10-19；修订日期：2023-02-13投稿网址：基金项目：国家自然科学基金项目（59479017）；铁道部科技研究开发项目（2010G003-F）第一作者：刘贯飞（1989-），男，博士研究生，主要从事土力学与深基础承载力的研究工作。E-mail：通讯作者：雷胜友（1965-），男，教授，博士，主要从事土力学、加筋土强度等方面的研究工作。E-mail： 第 50 卷 第 6 期水文地质工程地质Vol.50 No.62023 年 11 月HYDROGEOLOGY&ENGINEERING GEOLOGYNov.，2023solution for drained cylindrical cavity expansion in sand is obtained by solving the governing equationsnumerically.The results show that the solution of cylindrical cavity expansion established in this paper can beused in many types of sand by changing the values of stress-state parameter n and spacing ratio r*to select anappropriate shape of yield surface of sand,and the greater the values of n,r*are,the greater the initial yielddeviatoric stress of loose sand and the subsequent expansion pressure,but these situations will reverse in mediumdense and dense sand.The ultimate expansion pressure increases with the decreasing initial state parameter of thesand,and together with the volume variation rule of sand changes from always contraction into dilatation first andthen contraction,the radius of the elastic-plastic zone decreases first and then increases,and the hardeningresponse of sand changes from always hardening into softening first and then hardening.The ultimate expansionpressure increases with the increasing coefficient of earth pressure at rest of sand,but the volume change law ofsand has little change.This study provides a reliable theoretical support for the analysis of related geotechnicalengineering problem.Keywords：cylindrical cavity expansion；sand；unified state parameter model；the shape of yield surface；Lagrangian description 孔的扩张理论是岩土工程中的一个经典课题，因其形式简单、力学原理明确的优点，在桩基承载力1 3、静力触探4 5和旁压试验6分析中得到了广泛应用。为了合理描述材料弹塑性变形和剪切破坏等力学行为，临界状态模型近年来被广泛用于圆孔扩张的解答。例如 Collins 等7、Cao 等8分别使用原始剑桥（originalCam clay，OCC）模型和修正剑桥（modified Cam clay，MCC）模型，推导了饱和黏土中不排水圆孔扩张的近似解析解；Chen 等9使用大应变理论和拉格朗日分析法，基于 MCC 模型推导了初始应力各向异性条件下的饱和土中排水柱孔扩张半解析解；李镜培等10、Liu 等11、Castro 等12分别考虑 SMP 屈服准则、土的结构性和自然堆积各向异性等特点，采用类似 Chen 等9的方法，推导了不同情况下的饱和土中不排水或排水柱孔扩张解答。以上基于剑桥模型或改进型剑桥模型建立的饱和土中圆孔扩张解答，对孔扩张理论的发展具有较大的理论意义。然而，由于剑桥模型所用的相关联流动法则和形状固定的屈服面形式（OCC 模型为子弹头形，MCC 模型为椭圆形）13 14，导致其仅对正常固结和弱超固结饱和黏土具有较好的适用性，但不能准确模拟密砂的剪胀变形和松砂不排水剪切时加载初期的峰值偏应力15 16，因此，上述柱孔扩张解答均不适用于砂土。为了分析饱和砂中的柱（球）孔扩张力学过程，一些学者有针对地使用砂土类本构模型来描述砂的弹塑性变形特点。例如李林等17使用 Yao 等18提出的砂土临界状态模型推导了剪胀性砂土中的球孔排水扩张半解析解；Su 等19使用一个剪胀依赖状态的边界面模型推导了饱和砂中柱孔排水扩张半解析解；武孝天等20基于 CSUH 模型推导了同时适用于饱和黏土和砂土的柱孔不排水扩张解答。上述针对饱和砂的柱孔或球孔扩张解答分别考虑了密砂剪胀、松砂液化、砂土相对密度和静止侧压力系数等因素的影响，但忽略了饱和砂屈服面形状因素的影响，进而使以上解答尽管可以合理反映某一特定类别砂土在圆孔扩张时的剪切变形，但无法推广至其他类别砂土。研究表明，与黏土不同，砂土的屈服面形状通常不是固定的，而是要随土质情况的不同不断发生改变15,21 22。统一状态参数模型（clay and sand model，CASM）15可以通过改变土的间距比 r*和应力状态参数 n 的取值，使砂的屈服面形状发生改变，进而可用于计算不同类型饱和砂中的排水柱孔扩张力学过程。基于此，本文选择使用 CASM 和 Rowe 剪胀方程23来描述砂土的弹塑性变形，通过引入辅助变量，使用拉格朗日分析法建立饱和砂中的排水柱孔扩张半解析解。1 问题描述和本构模型介绍 1.1 柱孔扩张过程描述v0h0v0h0a如图 1 所示，无限均质饱和砂具有初始竖向应力（）、水平应力（=K0）和孔隙水压力（u0），以及体积（v0）。假设砂中有一柱孔在孔壁压力的作用下发生排水扩张，当孔壁应力从增大到时，柱孔孔径从 a0扩张至 a。在此过程中，孔壁处砂土最先屈服，之后随着砂土屈服区域越来越大，孔周围逐渐形成内侧为已屈服的弹塑性区和外侧未屈服的弹性区，100 水文地质工程地质第 6 期h0h0v0rzrz弹塑性区的半径假定为 rp。由于砂的渗透系数较大，一般认为饱和砂中的柱孔扩张属严格排水的情况。针对孔周某一砂颗粒的应力变化，在上述扩孔过程中，砂颗粒从其初始位置 rx0移动到了 rx，其中 x 为柱坐标系的 x 轴方向，应力状态从（,）变为（,），其中，、和分别表示砂的径向、环向和竖向有效应力，砂的体积从 v0变为 v。v0h0h0aa0arprp0rx0rxrp弹性区弹塑性区 图 1 砂中排水柱孔扩张示意Fig.1 Schematic diagram of drained cylindrical cavityexpansion in sand 将柱孔扩张视为平面应变过程，在柱坐标系中建立土的应力平衡方程。严格排水条件下扩孔时，饱和砂的孔隙水压力不发生改变，故应力平衡方程可以写为：drdr+rr=0（1）式中：d物理量的空间导数（欧拉描述），下同；r砂颗粒的径向位置/m。1.2 CASM如图 2 所示，根据 Been 等24的研究，砂的状态参数 的定义为 v-lnp平面内饱和砂当前状态与临界状态线（critical state line，CSL）之间的差值，表达式为：=v+ln(p)（2）式中：CSL 线的斜率；p=(r+z)/3p砂的平均有效应力/kPa，；p=1 kPa 时砂的临界体积。基于砂的状态参数概念和临界状态土力学理论，Yu15提出了一个能灵活描述黏性土和砂土受力变形特点的 CASM，近年来被许多学者用于分析饱和土的旁压试验、静力触探和桩基沉降等工程问题1,4,25 27。CASM 的屈服面函数为：F=qMp(1R)1n=0（3）或F=(qMp)n+ln(p/pm)ln(r)=0（4）式中：R参考状态参数；n应力状态参数；pm砂的硬化参数/kPa；M破坏应力比；r*间距比参数；q砂的偏应力/kPa，q=(r)2+(z)2+(rz)2。如图 2、图 3 所示，参数 n 控制着饱和土的屈服面形状，参数 r*控制着正常固结线（normal consolidationline，NCL）与 CSL 线之间的竖向距离。CASM 一般使用非相关联流动法则，本文选择基 NCLCSL（p,v）11Rv=1+eNlnp=0状态参数：=v+lnp间距比：r*=pm/pcs参考状态参数：R=（）lnr*当前应力状态Olnpcslnpmlnp图 2 土体状态参数的描述Fig.2 Description of state parameter of soilrpmpcs注：为土体回弹曲线的斜率；参数表示硬化参数与临界状态平均有效应力之间的比值；NCL 为正常固结线。On=4n=1n=2n=3r*=3q/kPapmpcs=pm/3p/kPa临界状态线屈服面图形图 3 不同 n 值情况下的 CASM 模型屈服面形状Fig.3 Yield surface shape of CASM with different n values2023 年刘贯飞，等：基于统一状态参数模型的砂中柱孔扩张分析 101 于 Rowe 剪胀方程23的塑性势函数来描述土的受力变形方向，表达式为：G=3Mln(p)+(3+2M)ln(3+2qp)(3M)ln(3qp)（5）式中：塑性势函数形状尺寸参数。根据式（5），土的塑性应变增量计算式为：Dpij=Gij=(Gppij+Gqqij)（6）式中：D物理量的物质导数（拉格朗日描述）；塑性因子；Dpij用张量表示的塑性应变增量。根据各向同性体积硬化规则，得到Dpm=vpmDpp=vpmGp（7）pp式中：砂的塑性体应变增量。图 2 中土体压缩回弹线的斜率。1.3 材料的本构关系根据广义胡克定律，砂的弹性本构关系为：Deij=1+EDijEDmmij（8）Deij式中：用张量表示的弹性应变增量；Dij用张量表示的有效应力增量；ij克朗内克符号；pv0E弹性模量/kPa，E=3(1)/；泊松比。根据砂土屈服面上的塑性一致性条件得到：DF=FijDij+FpmDpm=0（9）将式（7）代入式（9），整理后可得到塑性因子 和塑性模量 H 计算式为：=FijdijHH=FpmvpmGpFpm=1pmln(r)Gp=3(3+2M2q+3p3M3pq)在柱坐标系中建立砂的弹塑性应力应变关系式，将式（6）与式（8）相加，整理后可得到 CASM 的弹塑性应力应变关系为：Di=1+EDiEDp+1HFiGjDj（10）式中：Di柱坐标系 3 个主方向上砂的应变增量，i=r、z；DiDj（或）柱坐标系 3 个主方向上砂的应 力增量，i 或 j=r、z；式中其他物理量的计算式分别为：Fi=13(1plnrnqnMnpn+1)+3(i p)2q(nqn1Mnpn)Gj=3+2M2q+3p3M3pq+3(j p)2q2(3+2M)2q+3p+3M3pq 2 砂中柱孔排水扩张半解析解 2.1 柱孔扩张弹性解柱孔周围土体发生纯弹性变形时，由于变形量很小，可认为土的体积和平均有效应力不发生改变，根据式（8）和式（1）容易解得砂中柱孔扩张弹性解为9 12：r=h0+(rph0)(rpr)2（11a）=h0(rph0)(rpr)2（11b）z=v0（11c）ure=rph02G0r2pr（11d）rp式中：弹塑性区交界面处砂的径向有效应力/kPa；rp弹塑性区的半径/m；ure弹性区内砂的径向位移/m；G0剪切模量/kPa，G0=E/2（1+）。2.2 柱孔扩张弹塑性解 2.2.1 砂的有效应力和体积计算式弹塑性区内的砂土发生大变形，用对数应变表示砂的径向应变（r）和环向应变（），表达式分别为9：r=ln(dr/dr0)（12a）=ln(r/r0)（12b）柱孔扩张是平面应变问题，竖向应变 z=0，体应变 v=r+=ln(v/v0)，因此 r也可表示为：102 水文地质工程地质第 6 期r=v=ln(v/v0)ln(r/r0)（13）孔周砂土的应变表达式（式 13、12b）和应力平衡方程（式 1）均为欧拉坐标描述，而弹塑性应力应变关系（式 10）为拉格朗日坐标描述，为便于推导，需使用辅助变量将式（13）、式（12b）和式（1）也转化为用拉格朗日坐标描述的形式。参考 Chen 等9的方法，引入辅助变量，定义为：=urpr=rr0r（14）式中：urp弹塑性区内任一位置点处砂颗粒的径向 位移/m。将式（14）代入式（12b）、式（13），可以得到：=ln(rr0)=ln(1rr0r)=ln(1)（15a）r=v=ln(v/v0)ln(1)（15b）对式（15a）（15b）两边微分，得到：D=D1，Dr=DvvD1将 Dr、D和 Dz=0 代入式（10），可解得：DrD=1c11fv+c12c111（16a）DD=1c21fv+c22c211（16b）DzD=1c31fv+c32c311（16c）DvD=f（16d）式中：c11=12+Ehbza+Ehazb+Ehazbz+Ehab；c12=Ehar(b+bz)+(1+Ehbaz+Ehazbz)；c21=Ehar(a+az)+(1+Ehabz+Ehazbz)；c22=12+Eharbr+Eharbz+Ehbraz+Ehazbz；c31=Ehbr(a+az)+(1+EhabEhbaz)；c32=+2+EharbrEhbazEharbEhbraz；=1+E122+Eh(1)(arbr+ab+azbz)+Eh(abr+bar+azb+bza+azbr+bzar)塑性势函数形状尺寸参数；pp砂的塑性体应变增量。rz式（16）所示的微分方程组有 5 个未知量（、v、f）和 4 个控制方程，为了求解式（16），还需增加一个控制方程。令 对 r 求导，得到：ddr=1r(durpdr)（17）将式（17）代入式（1），可以得到：DrD(1r(durpdr)+rr=0（18）根据式（13a），式（18）中的 durp/dr 可表示为：durpdr=1er（19）将式（15b）代入式（19），而后将得到的 durp/dr 代入式（18），之后联立式（16a）可以解得：f=vc11(r1v0/v(1)c12c11(1)（20）将式（20）代入式（16），得到：DrD=rN1（21a）DD=c21c11(rN1c12c11N2)c22c21N2（21b）DzD=c31c11(rN1c12c11N2)c32c31N2（21c）DvD=vc11(rN1c12c11N2)（21d）N1=1v0/v(1)N2=(1)其中，。r(p)(p)z(p)v(p)r()()z()v()p fp f式（21）为弹塑性区内基于拉格朗日描述的砂土体积和有效应力求解方程组，适用于弹塑性区内任意径向位置处的饱和砂土。结合初始条件、和，可计算得到不同 值情况下砂的、和。的取值范围为，对应砂颗粒刚开始发生塑性变形时的，对应砂颗粒扩孔结束时的。2.2.2 边界条件r(p)=(p)=h0z(p)=v0v(p)=v0对于正常固结砂，砂的初始应力状态即为其发生弹塑性变形的起始条件，、。对于超固结砂，根据式（4）、式（11）和偏应力 q 的定义，可以解得：(p)=h013(q(p)2(h0v0)2)（22a）r(p)=h0+13(q(p)2(h0v0)2)（22b）z(p)=v0（22c）2023 年刘贯飞，等：基于统一状态参数模型的砂中柱孔扩张分析 103 v(p)=v0（22d）p=(urr)r=rp=r(p)h02G0（22e）q(p)=Mp0ln(Rpm0/p0)ln(r)1n=Mp010R1n（23）q(p)式中：砂的初始屈服偏应力/kPa。R=pC/pm0R砂土的超固结比，；pC砂的历史最大硬化参数/kPa。pC=rp0exp(0)（24）pm0式中：砂土在当前应力状态对应的虚拟硬化参 数/kPa，pm0=p0exp(lnr(q0Mp0)n)（25）2.2.3 计算结果的转换rz由式（21）解得的、均是关于辅助变量 的拉格朗日描述。为便于分析，还须将计算结果变换为关于孔周径向坐标 r 的欧拉描述，根据式（15b）（17）（19）可得：drr=d1v0/(v()(1)（26）式中：dr弹塑性区内相邻 2 个砂土颗粒之间的径 向距离/m。p从到 对式（26）积分，可以得到欧拉坐标 r 与拉格朗日坐标 之间的对应关系为：ra=expwad1v0/v()(1)（27）a式中：孔壁处的 值。pa从到对式（26）积分，可得到孔周的弹塑性区半径 rp为：rpa=expwpad1v0/v()(1)（28）在式（14）中，令 r=a、r0=a0，变化后可得到 与实时孔径 a 之间的关系为：a=a0/(1)（29）3 参数分析r(p)(p)z(p)v(p)p假设柱孔初始孔径 a0=0.01 m，下同；使用 R 语言编写程序计算饱和砂中的排水柱孔扩张力学过程。根据砂的初始应力状态，首先使用 2.2.2 节中所述办法计算砂初始屈服时的、和；之后使用四阶龙格库塔算法求解式（21），并且用式（26）（29）对所得结果进行转化；最后结合由式（11）所得的弹性区内的应力解答，给出完整的砂中柱孔扩张解答。cscsv0以 Ticino 砂为例，分析屈服面形状参数 n、r*，初始状态参数 0和静止侧压力系数 K0对砂中排水柱孔扩张结果的影响。根据 Yu15，Ticino 砂的临界状态参数为 =1.986，=0.024，=0.008，=0.3，M=1.29，=32o；实验室测得砂的最大孔隙比 emax和最小孔隙比emin分别为 0.89 和 0.6028；另外，下文无特殊说明时，饱和砂的静止侧压力系数取 K0=1 sin=0.47，=200 kPa，n=2，r*=108.6。3.1 初始状态参数 0v0针对多种不同初始状态参数 0的 Ticino 砂，取 n=2、r*=108.6、=200 kPa，用文中方法计算不同 0值饱和砂中柱孔扩张时砂的力学指标变化，所得结果如图 4、图 5 所示。a/z0a/z0a/z0由图 4（a）可知，5 种不同 0值情况下扩孔所得均随 a/a0的增大而增大，当 1a/a02 时，增大的较快，之后越来越慢并且在 a/a0=4 以后达到稳定；0越小，砂土越密实，同样扩孔半径a/a0时的越大。由图 4（b）可知，当 0=0.044 6 时，饱和砂在扩孔过程中一直剪缩并逐渐靠近临界状态线；当 0=0 时，砂的剪缩过程为先远离临界状态线而后再逐渐靠近后者；当 0=0.1，0.2，0.3 时，砂土在剪切初期会发生剪胀，这是密实砂土受剪时的砂颗粒错位引起的，并且 0越小，砂的剪胀量越大，而当体积路径越过临界状态线后，初始的中密砂和密砂此时已近似达到了临界状态，在后续扩孔过程中，体积路径为紧贴临界状态线发生剪缩。扩孔结束时，砂的体积 v 随 0的减小而减小，对应的p随0的减小而增大。另外，在图4（b）中，5 种砂的初始孔隙比分布为 0.570.91，这覆盖了大多数的 Ticino 砂天然孔隙比28。pmpmpm如图 4（c）所示，针对 6 种饱和砂，用于表示砂土屈服面大小的硬化参数在扩孔过程中均发生了较大变化。当 0=0.044 6 和 0.02 时，随 a/a0的增大而增大，说明这 2 种砂土在扩孔过程中一直在硬化；当 0=0，0.1，0.2，0.3 时，随 a/a0的增大先减小后增大，说明这 4 种砂土在扩孔过程中先软化而后再硬化。如图 4（d）所示，不同情况下的 rp/a 均随 a/a0的增大而增大，增速先快后慢，并且在 a/a0=6 以后几乎不 104 水文地质工程地质第 6 期再变化。当 0从 0.04 降低到 0 时，最终的 rp/a 随 0的减小而减小，而后当 0从 0 降低到0.3 时，rp/a 又随 0的减小而增大。这是因为当 0 0 时，0越大，砂土越易在低应力时屈服，进而弹塑性区的面积越大；而当 0 0 时，0越小，扩孔时砂的剪胀量越大，弹塑性区的面积越大。rzr图 5 给出了 3 种 Ticino 砂中柱孔排水扩张 a/a0=2 时，孔周砂土的径向、环向、竖向有效应力（、）和体积（v）的分布。可以看到，3 种砂的、和zrz在孔壁附近均随 r/a 的增大而减小，说明砂土未达到临界状态。当 0=0.044 6 时，扩孔以后的孔周只有弹塑性区，联立式（23）（24）可知，此时砂土的初始状态对应为正常固结土，其从扩孔一开始即发生弹塑性变形9；当 0=0.1 和0.3 时，扩孔以后的孔周可分为弹塑性区和弹性区 2 个区。和在弹塑性区内随r/a 的增大快速减小，而的变化相对较小，这是因为柱孔扩张时的孔周土体会因为环向受拉而进入主动状态，进而引起了的减小；同时，因为砂中超孔隙水 1234567 8 9 100.1110100a/a0a/a0a/a00=00=0.10=0.20=0.3为abaqus结果（0=0.044 6）为abaqus结果（0=0.1）为abaqus结果（0=0.3）为abaqus结果（0=0.044 6）为abaqus结果（0=0.1）为abaqus结果（0=0.3）456789101.451.501.551.601.651.701.751.801.851.901.95lnpv临界状态线12345678910024681012141618rp/a0=0.04 0=0.020=00=0.10=0.20=0.31234567 8 9 101091081071061051041031021011010pm/z0a/z0（a）孔壁压力a/z0随a/a0的变化0=0.044 60=00=0.10=0.20=0.30=0.044 60=00=0.10=0.20=0.30=0.044 6（b）孔壁处砂土的体积路径（c）硬化参数pm/z0随a/a0的变化（d）弹塑性区半径rp/a随a/a0的变化图 4 0值对 Ticino 砂中柱孔扩张结果的影响Fig.4 Effects of 0 value on the expansion results of cylindrical cavity in Ticino sand 11010000.51.01.52.02.53.03.54.04.51.801.821.841.861.881.901.921.94r/ar/ar/av土体应力/z0土体应力/z0土体应力/z0实线为本文解答结果“”为abaqus结果实线为本文解答结果“”为abaqus结果实线为本文解答结果“”为abaqus结果v弹塑性区11010002468101214v1.721.741.761.781.801.821.84v弹性区弹塑性区110100505101520253035v1.451.501.551.601.651.701.751.801.85v弹性区弹塑性区（a）0=0.044 6（b）0=0.1（c）0=0.3rrrzzz图 5 a/a0=2 时孔周的径向、环向、竖向有效应力分布Fig.5 Distributions of radial,tangential and vertical stresses and specific volume around the cavity when a/a0=22023 年刘贯飞，等：基于统一状态参数模型的砂中柱孔扩张分析 105 rz压力的完全消散，土中应力向砂颗粒转移，这也会引起、和发生相应的改变。弹性区内的砂土体积不发生变化，在弹塑性区内，松砂（0=0.044 6）均发生了剪缩，最终的砂土体积随 r/a 的增大而增大，中密砂（0=0.1）和密砂（0=0.3）在扩孔过程中先剪胀后剪缩，最终砂土的体积随 r/a 的增大先增大后减小。为了验证文中方法所得结果的合理性与正确性。使用 abaqus 软件对 0=0.044 6、0.1 和0.3 的 Ticino砂中柱孔扩张力学过程进行了计算。为了能够在 abaqus软件中使用 CASM 模型同时又便于程序实现，基于显示积分法27编写了相应的 umat 子程序。在轴对称坐标系中建立有限元模型，柱孔初始孔径 a0=0.01 m，模型的径向最大计算区域（右边界）取 400a0来消除边界效应，采用长方形网格，网格仅在靠近孔壁处加密以节约计算时间，采用 CAX8RP 单元。柱孔扩张分为初始、地应力平衡和扩孔 3 个分析步，其中在扩孔分析步中对左边界施加径向位移增量 9a0，同时将右边界设为排水条件，扩孔分析步的时间取 1，砂的渗透系数取 0.01 m/s。图 4（a）、4（b）和图 5 中给出了使用 abaqus软件所得结果与本文结果的对比，可以看到，本文方法针对上述 3 种饱和砂中柱孔扩张的计算结果与 abaqus软件所得是非常接近的。3.2 应力状态参数 n 和间距比 r*针对 0=0.02、0.1 和0.3 的 3 种 Ticino 砂，取 n=1、2、3、4，r*=108.6，用文中方法计算得到不同 n 值情况下砂中排水柱孔扩张结果，如图 6、图 7 所示；取 r*=10，108.6，1 000，10 000，n=2，用本文方法计算得到不同 r*值时的砂中柱孔扩张结果，如图 8、图 9 所示。1234567 8 9 1010100n=1n=2n=3n=4n=1n=2n=3n=40=0.020=0.020=0.10=0.10=0.30=0.3a/a0（a）孔壁压力a/z0随a/a0的变化456789101.551.601.651.701.751.801.851.901.95lnpv临界状态线（b）孔壁处砂土的体积路径a/z0 图 6 n 值对 Ticino 砂中柱孔扩张结果的影响（r*=108.6）Fig.6 Effects of n value on the expansion results of cylindrical cavity in Ticino sand(r*=108.6)110110临界状态线有效应力路径1100110100临界状态线有效应力路径11010100110100临界状态线有效应力路径n=1n=2n=3n=4n=1n=2n=3n=4n=1n=2n=3n=4（a）0=0.02（b）0=0.1（c）0=0.3q/p0q/p0q/p0p/p0p/p0p/p0 图 7 不同 n 值时孔壁处砂土的 p-q 应力路径（r*=108.6）Fig.7 Stress path of sand at cavity wall in p-q space for different n values(r*=108.6)a/z0a/z0a/z0图 6（a）和图 8（a）分别给出了不同 n 值和不同 r*值情况下柱孔扩张所得随 a/a0的变化。可以看到，对于上述 3 种 Ticino 砂，改变 n、r*的取值对扩孔时的均有较大影响，而且单独改变 n 和 r*的值对所得值的影响规律是相似的。对于同样的 a/a0，a/z0a/z0a/z0当 0=0.02 时，随 n、r*的增大而增大；当 0=0.1 和0.3 时，随 n、r*的增大而减小；0越小，砂土越密实，改变 n、r*的取值对所得的值影响越大。图 6（b）和图 8（b）分别给出了不同 n 值和不同 r*106 水文地质工程地质第 6 期值饱和砂中排水柱孔扩张时孔壁处砂土的体积路径。如图所示，对于每一种 Ticino 砂，n、r*的不同取值使饱和砂具有不同的体积路径。当 0=0.02 时，不同n、r*值时的饱和砂在柱孔扩张时一直处于剪缩状态，砂的状态参数 在此过程中不断减小，直至无限接近于 0；n、r*值越小，饱和砂在低 p时的剪缩量越大，进而使砂土可以在低应力条件下近似达到临界状态，扩孔至 a/a0=10 时，砂土体积 v 随 n、r*的增大而减小，对应的 p随的 n、r*增大而增大。当 0=0.1 和0.3 时，不同 n、r*值时的饱和砂在柱孔扩张时表现为先剪胀，而后在其体积路径越过临界状态线后紧贴后者不断剪缩，状态参数 在扩孔过程中先增大至稍大于 0，而后再不断减小；n、r*值越小，饱和砂在低 p时的剪胀量越大，但随着 p的增加，剪胀量的增加幅度不断减小，进而导致总的剪胀量与 n、r*值较大时的相等甚至被后者超越；在砂的剪缩阶段，n、r*值越小，剪缩量越大，相应 p增大的也越多，扩孔至 a/a0=10 时，砂的体积 v 随 n、r*的增大而增大，对应的 p随 n、r*的增大而减小。3 种 Ticino 砂中排水柱孔扩张时，不同 n 值和不同r*值情况下孔壁处砂土的 p-q 应力路径分别如图 7、图 9 所示。可以看到，对于每一种 Ticino 砂，n、r*的不同取值使砂土具有不同的应力路径，而且单独改变 n和 r*的值对应力路径的影响规律是相似的。3 种 Ticino砂的 p-q 应力路径初始阶段均为一条向上的竖线，p不变，q 一直增大，说明此时砂土发生了纯弹性变形。当 0=0.02 时，不同 n、r*值情况下的砂土初始屈服偏应力均位于临界状态线以下，之后的弹塑性变形应力路径从临界状态线的下方逐渐接近后者，并且 p、q 一直增大；当 0=0.1 和0.3 时，不同 n、r*值时的砂土初始屈服偏应力均位于临界状态线上方，后继的弹塑性变形应力路径从临界状态线的上方逐渐接近并越过后者，之后再紧贴后者不断上升，在此过程中，砂的 p、q 值总体上也一直增大，但 0=0.3 时的砂土在 n=1、r*=108.6 和 n=2、r*=10 两种情况下会因为剪切初期的偏应力 q 过大，造成应力路径在越过临界状态线之前的 q 值保持不变或者先减小而后再保持不变。另外，结合式（3）、式（23）可知，当 0=0.02 时，1234567 8 9 1010100a/a0r*=10 r*=108.6 r*=1 000r*=10