基于抛物线强度准则的冻土热‒弹性‒塑性柱孔扩张模型

Thermo‒Elasto‒Plastic Cavity Expansion Model of Frozen Soil Based on Strength Criterion with a Parabolic Curve

摘要

柱孔扩张作为典型的边值问题，可以用来分析冻土区的桩基工程在施工时土体受到径向压力时的应力和位移。在施工时随着孔内径向压力的增加，周围土体会逐渐由弹性变形状态转变为弹塑性变形状态，进而可将周围土体划分为弹性区和塑性区，其应力状态由不同的控制方程确定。现有的柱孔扩张模型适用于未冻土和岩石，其本构方程和强度准则仅与材料自身的力学特性有关，而冻土的柱孔扩张模型还与土体的温度分布有关。此外，冻土的强度准则还具有与未冻土明显不同的非线性特征，即随着平均应力的增大，冻土的强度先增大后减小，通常呈现抛物线形式。本文提出的冻土热‒弹性‒塑性柱孔扩张模型基于适用于冻土的抛物线强度准则，同时考虑了温度分布对冻土力学特性的影响，使用了包含温度变量的热弹性本构方程对土体的应力状态进行分析。在计算塑性区位移时，本文模型运用连续介质力学的运动方程推导了柱孔内侧位移的表达式。通过算例分析，可以发现弹塑性分界面上环向应力的不连续、塑性区的变形存在压硬性等，与冻土的实际力学特性相符合。当孔内温度变化，塑性区的范围和应力分布均会发生改变。本文模型对冻土地区的原位土工试验、桩基础设计和施工有指导意义。

Abstract

Objective The construction of pile foundation projects in cold regions requires quantitative analyses of displacement and stress in the surrounding soil, which can be achieved using cavity expansion methods. Classical cavity expansion methods fail to account for uneven temperature distribution, temperature variation, and the unique mechanical behavior of frozen soil. Therefore, this study proposes a thermo‒elasto‒plastic model based on a strength criterion with a parabolic curve, also referred to as the parabolic strength criterion, which is applicable to frozen soil and provides analytical expressions for displacement and stress. Methods Firstly, the geometry and boundary conditions were determined, and several fundamental assumptions of the proposed model were introduced. A Dirichlet boundary condition, defined as a fixed surface temperature, was considered when addressing the heat diffusion problem. Secondly, the governing equations of the thermo‒elastic cavity expansion model were introduced and systematically combined, through which the explicit expressions for displacement and stress in the elastic state were derived. With the increase in internal radial pressure, the state of the surrounding soil gradually transitioned from an elastic state to an elasto‒plastic state, whereby the surrounding soil was divided into elastic and plastic zones governed by different equations. Thirdly, to account for the nonlinear strength behavior of frozen soil, the parabolic strength criterion was introduced, indicating that the strength initially increased and decreased as the mean stress increased. The stress distribution within the plastic zone was derived by combining the equilibrium equations of continuum mechanics with the parabolic strength criterion of frozen soil. It was clearly impossible to derive an explicit expression for displacement using conventional plasticity theory because the parabolic strength criterion was based on triaxial experimental results. Therefore, the proposed model assumed that the compressive equations remained unchanged after the transition from the elastic state to the plastic state. The expression for the displacement of the internal wall surface was derived by substituting the compressive equations into the kinematic relations of continuum mechanics. Finally, an example involving pile foundation construction was analyzed using the proposed model, and the calculation results were compared to those obtained from a precedented model and a numerical method. Results and Discussions Comparing the modeling procedure with that of the classical cavity expansion model indicated that both stress and displacement were directly interconnected with the temperature distribution. A logarithmic temperature distribution was first obtained by solving the heat diffusion equation, and this distribution was used in the derivation of elastic stress and displacement. For stress and displacement in the elastic state, the influence of temperature was clearly reflected in the term containing the thermal expansion coefficient, which was generally small for solids such as ice, rock, and soil. A significant finding was that when compressive stress was defined as positive, the radial stress consistently represented the maximum principal stress, while the circumferential stress represented the minimum principal stress. This finding facilitated the determination of the plastic zone and enabled the substitution of principal stresses into the parabolic strength criterion. Within the plastic zone, the influence of temperature became considerably more complex because temperature variations not only altered the parameters of the parabolic strength criterion but also affected the subsequent calculations of stress and displacement. The application of the parabolic strength criterion made it challenging to obtain explicit expressions for stress and displacement in the plastic zone, whereas quantitative results were achievable through numerical analysis. The quantitative results of the example demonstrated several distinctive characteristics of the proposed model. For instance, with increasing internal radial pressure, the plastic boundary gradually moved outward, reaching 1 m when the pressure was 1 690 kPa and 2 m when the pressure was 6 306 kPa. The quadratic approximation of the plastic boundary coordinate exhibited good agreement with the results. The circumferential stress at the boundary surface of the plastic zone was found to be discontinuous, which was attributed to stress redistribution in the frozen soil. The displacement of the internal wall surface was calculated and was 83.5 mm when the pressure was 1 690 kPa and 170.0 mm when the pressure was 5 000 kPa. The characteristic of compressive hardening was clearly observed from the stress‒displacement curve. When the results were compared to those of the precedented model based on the Tresca strength criterion and the numerical method based on PFEM, the stress‒displacement relationship was considered acceptable from a geotechnical perspective. The sensitivity of model parameters was further analyzed and was primarily associated with the effects of temperature and the coefficient of plasticity-growing. Internal temperature groups of ‒5 ℃, ‒10 ℃, and ‒15 ℃ were respectively applied, resulting in smaller plastic zones, indicating higher frozen soil strength, and lower circumferential stress drops at the plastic boundary. In addition, the coefficient of plasticity-growing was shown to be effective in further adjusting the calculated results. Conclusions The proposed model is developed based on a strength criterion characterized by a parabolic curve applicable to frozen soil, with temperature variation explicitly considered. A thermoelastic constitutive equation is employed in which temperature is treated as a variable to derive the stress and displacement expressions of the soil. The model can incorporate more complex construction conditions and is verified to be effective and rational from a geotechnical perspective. The results provide valuable guidance for in situ soil testing, as well as for the design and construction of pile foundations in frozen soil regions.

Graphical abstract

关键词

柱孔扩张 / 冻土 / 温度场 / 强度 / 弹‒塑性分析

Key words

cavity expansion / frozen soil / temperature field / strength / elasto‒plastic analysis

冻土是一种温度低于0 ℃且含有冰的土岩，可根据冻融特征分为短时冻土、季节冻土、多年冻土（又称永久冻土）。其中，多年冻土的面积约占全球陆地面积的23%，占北半球裸露陆地面积的20%，主要分布在俄罗斯、加拿大、中国和美国等地^[1‒2]。由于冻土中的土颗粒为冰所胶结，其力学性质与一般工程中的三相土有较大的差异，尤其体现在温度、多晶冰的性质对冻土强度的影响上。

冻土地区的工程建设面临诸多挑战，包括在施工过程中冻土力学性能的改变^[3]、温度变化导致的冻融过程^[4]等。寒区工程建设实践表明，桩基础可以利用桩基嵌入基岩或多年冻土层，从而得到较高的承载力并且地温场变化较小，因此冻土区的建筑物基础多采用桩基^[5]。对于预制混凝土桩，由于在成桩过程中不断扰动、挤压周围土体，桩基与土体接触的界面温度会升高，同时周围土体还将受到逐渐增大的径向压力，如图1所示。对于多年冻土地区的桩基础工程，上述工况可以运用非均匀受热的柱孔扩张模型进行分析。因此，建立非均匀温度场下的冻土柱孔扩张模型对寒区桩基础工程建设具有重要意义。

对于边值问题，岩石或土体的柱孔扩张以厚壁圆筒的弹性分析为基础，在原位土工试验^[6]、桩基工程^[7‒8]、隧道工程^[9]等领域得到了广泛的应用。柱孔扩张理论的弹性解在许多文献中都有介绍，而随着内部荷载的进一步增加，土体进入塑性变形阶段，则需要引入材料的屈服准则、塑性本构方程进行求解。Hill^[10]推导了在弹塑性变形条件下Tresca材料应力和位移的解析解，也有学者采用了更适用于土体的Mohr‒Coulomb准则，推导出了土体在弹塑性变形下的解析解^[11‒12]。针对不同的土质，蒋明镜等^[13]考虑了材料应变软化特征，得到了土体柱孔扩张的理论解。也有学者考虑了土的固结，针对不同的排水条件提出了不同的模型，例如：Chen等^[14]提出的不排水条件下的柱孔扩张模型；Mo等^[15]提出的排水条件下的柱孔扩张模型。Rosenberg等^[16]根据von‒Mises屈服准则引入了归一化的压力与速度，推导了不同材料在弹塑性条件下的球形孔扩张理论解。张家奇等^[17]基于Mohr‒Coulomb准则，研究了考虑弹塑性卸荷的柱孔和球孔反向扩张解。上述推导过程中的弹性状态均将土体视为均匀受热的材料，然而，在实际工程中，施工机械散热或环境温度变化均会导致柱孔内部温度变化，使得土体处于非均匀受热状态，从而产生温度应力。此外，冻土的强度普遍较高，但受温度影响大^[3]，这是因为冻土中含有冰胶结物，其内聚力与温度有关^[5]。上述力学特性使得冻土的强度准则与适用于金属材料的Tresca准则、适用于一般岩土体的Mohr‒Coulomb准则有较大差异。因此，在建立冻土的柱孔扩张模型时，不仅需要考虑周围土体的温度场，还要引入适用于冻土的强度准则，以便对实际桩基工程中的应力和位移进行更精确的分析。

本文率先将适用于冻土的抛物线强度准则引入柱孔扩张问题的分析中，结合冻土受温度场影响的力学特性，构建了考虑非均匀温度分布的冻土热‒弹性‒塑性柱孔扩张模型，系统推导了模型的控制方程、应力与位移解析表达式，并通过工程算例验证模型的合理性与适用性，探究温度、内压等关键因素对冻土弹塑性分区、应力分布及位移变化的影响规律，以期为冻土地区桩基工程的设计、施工及原位土工试验提供理论支撑与技术参考。

1 问题分析

1.1 模型建立

桩基础施工中开挖形成的柱孔可视为从半径为b的圆柱体中挖去一个半径为a的圆柱体，如图2所示。因此选用柱坐标系(r,θ,z)描述周围土体的位移、应力较为方便。图2中，p_a 为内压，p_b 为外压，r为土体半径，s为塑性区外边界，T(a)为孔内温度，T(b)为孔外温度。

如果在柱坐标系下考察土体中任意一点的位置，则位移矢量 u 只是r的函数，且只有径向的分量，即 u = u(r) e_r。为了表述方便，下文都只用u来表示某一点的位移。由于桩身长度通常远大于桩的直径，故此类柱孔扩张问题可以视为平面应变问题，即竖向应变ε_z =0。

在柱坐标系描述的轴对称问题中，不为0的应力分量只有σ_r 、σ_θ 、σ_z，同时上述应力分量也是3个方向的主应力。

1.2 边界条件

在多年冻土区的桩基础施工过程中，以开挖形成的柱形孔为内边界，其周围的土体受到扰动从而发生应力状态的改变。扰动对边界条件的影响主要包括孔内温度T(a)的升高和孔内压力p_a 的增大。由于永久冻土区的地面温度总是低于冰的熔点，因此在自然条件下的永久冻土不发生相变。由于桩基的埋深一般不超过20 m，因此冻土随深度的温度变化并不明显^[1,18]。这里只考虑施工过程中不发生相变的情况，即T(a)总是小于0 ℃。同时可以看出这里的温度边界属于第一类边界条件。

在孔内压力p_a 逐渐增大的过程中，可以将周围土体由里到外分为3个区域：塑性区（a < r < s）、弹性区（s < r < b）、未扰动区（r > b）。弹性区的外边界具有自然条件下的温度T(b)以及自然条件下的地应力p_b，可以通过地质资料或现场试验测得。上述边界条件如图2所示。

1.3 计算假设

模型采用弹塑性力学、连续介质力学框架下的所有假定。除此之外的主要假定为：

1）应力与温度为单向耦合，即温度变化导致应力发生变化，而应力变化不影响土体温度分布，具体描述见第2节。对于冻土的温度分布，可以仅由固体的热传导方程确定^[1]，具体描述见第2节。

2）温度分布T(r)仅由柱孔的内外温度确定（第一类边界条件）。

3）土体处于平面应变状态，即ε_z = 0，后续所有公式中的应力和应变均以压为正。

4）塑性变形状态下，土体的本构方程发生改变，压缩速率与弹性变形状态时保持相同的形式。

2 柱孔扩张的热弹性基本方程

1）几何方程

平面应变条件下，轴对称问题的几何方程为：

{\epsilon }_r=\frac{\partial u}{\partial r},{\epsilon }_{\theta }=\frac{u}{r}（1）

式中，ε_r 为径向应变，ε_θ 为环向应变。

2）平衡方程

对于平面应变条件下的轴对称问题，各向同性物体的平衡方程用应力分量σ_r 、σ_θ 、σ_z 表示为：

\frac{\partial {\sigma }_r}{\partial r}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=\mathrm{0}（2）

3）热弹性本构方程

考虑非均匀受热对土体应力状态的影响，热‒弹性本构方程（平面应变条件下）用量σ_r 、σ_θ 、σ_z 表示为：

\left\{\begin{array}{l}{\sigma }_r=\frac{E}{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}((\mathrm{1}-\nu ){\epsilon }_r+\nu {\epsilon }_{\theta })-\frac{\alpha ET}{\mathrm{1}-\mathrm{2}\nu },\\ {\sigma }_{\theta }=\frac{E}{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}((\mathrm{1}-\nu ){\epsilon }_{\theta }+\nu {\epsilon }_r)-\frac{\alpha ET}{\mathrm{1}-\mathrm{2}\nu },\\ {\sigma }_z=\frac{E\nu }{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}({\epsilon }_r+{\epsilon }_{\theta })-\frac{\alpha ET}{\mathrm{1}-\mathrm{2}\nu }\end{array}\right.（3）

式中，α为热膨胀系数，E为弹性模量，ν为泊松比。更确切地说，由于不同埋深的土体固结程度不同，弹性模量和泊松比分别为E(h)、ν(h)，其中，h为深度。在实际工程应用中，上述参数可以通过现场试验确定。

4）热传导方程

如果土体中某一点的温度T(r)不随时间变化（稳定温度场），则土体中的温度分布符合固体的热传导方程：

{\nabla }^{\mathrm{2}}T(r)=\frac{{\mathrm{d}}^{\mathrm{2}}T}{\mathrm{d}{r}^{\mathrm{2}}}+\frac{\mathrm{1}}{r}·\frac{\mathrm{d}T}{\mathrm{d}r}=\mathrm{0}（4）

式中，▽为Nabla算符。对于第一类边界条件，如果内温为T(a)，外温为T(b)，则上述微分方程的解为：

T(r)=\frac{T(b)-T(a)}{\mathrm{l}\mathrm{n}b-\mathrm{l}\mathrm{n}a}\mathrm{l}\mathrm{n}r+\frac{T(a)\mathrm{l}\mathrm{n}b-T(b)\mathrm{l}\mathrm{n}a}{\mathrm{l}\mathrm{n}b-\mathrm{l}\mathrm{n}a}（5）

3 弹性变形状态下的位移与应力分析

在内压p_a 较小的情况下，周围土体整体处于弹性状态，故可用弹性理论的方程求解其位移和应力。弹性方程有多种求解方法，此处采用位移法求解可以同时解得位移和应力^[19‒20]。联立式（1）～（3），得到用位移表示的弹性平衡方程为：

\frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{\mathrm{1}}{r}·\frac{\mathrm{d}(ru)}{\mathrm{d}r}\right)=\frac{\alpha (\mathrm{1}+\nu )}{\mathrm{1}-\nu }·\frac{\mathrm{d}T}{\mathrm{d}r}（6）

在弹性理论中只考虑相对温度，这是因为均匀受热的物体在无外力的情况下不产生应变。此处可以定义起算温度T(b)=0，从而由式（5）和（6）解得弹性状态下位移u的表达式为：

u=\frac{\alpha (\mathrm{1}+\nu )}{\mathrm{1}-\nu }\left(C_{\mathrm{1}}r+\frac{C_{\mathrm{2}}}{r}+\frac{T(a)}{\mathrm{2}\mathrm{l}\mathrm{n}\frac{b}{a}}r\left(\mathrm{l}\mathrm{n}\frac{b}{r}+\frac{\mathrm{1}}{\mathrm{2}}\right)\right)（7）

式中，C₁与C₂为积分常数。通过几何方程（1）和本构方程（3）可以求出弹性状态下应力的表达式为：

{\sigma }_r=C_{\mathrm{3}}+\frac{C_{\mathrm{4}}}{r^{\mathrm{2}}}-\frac{\alpha ET(a)}{\mathrm{2}(\mathrm{1}-\nu )\mathrm{l}\mathrm{n}\frac{b}{a}}\left(\mathrm{l}\mathrm{n}\frac{b}{r}+\frac{\mathrm{1}}{\mathrm{2}}\right)（8）

{\sigma }_{\theta }=C_{\mathrm{3}}-\frac{C_{\mathrm{4}}}{r^{\mathrm{2}}}+\frac{\alpha ET(a)}{\mathrm{2}(\mathrm{1}-\nu )\mathrm{l}\mathrm{n}\frac{b}{a}}\left(-\mathrm{l}\mathrm{n}\frac{b}{r}+\frac{\mathrm{1}}{\mathrm{2}}\right)（9）

{\sigma }_z=\mathrm{2}{C}_{\mathrm{3}}-\frac{\alpha ET(a)}{(\mathrm{1}-\nu )\mathrm{l}\mathrm{n}\frac{b}{a}}\mathrm{l}\mathrm{n}\frac{b}{r}（10）

式（8）～（10）中，C₃、C₄为积分常数，由柱孔内外的应力边界条件σ_r (a)=p_a 和σ_r (b)=p_b 确定。

4 弹塑性分析

4.1 冻土的强度准则

在一般工程中，土的强度通常用线性的Mohr‒Coulomb准则描述，即屈服面为大主应力σ₁和小主应力σ₃的线性组合。如果采用广义的Mohr‒Coulomb准则，则屈服面为平均应力p与广义剪应力q的线性组合。而根据冻土的大量试验数据，随着围压的增大，冻土的强度会出现峰值，并且随着围压的进一步增大而降低，这一现象说明冻土的强度是非线性的^[21‒22]。马巍等^[23]根据兰州砂土的试验数据，通过引入平均应力的二次项p²来考虑冻土的非线性强度准则^[23]，即：

q=c-kp+\frac{k}{\mathrm{2}{p}_{\mathrm{m}}}{p}^{\mathrm{2}}（11）

式中：c为黏聚力；k=tan φ，φ为内摩擦角；p_m为当冻土的剪切强度q达到最大值时的平均应力值。由式（11）定义的强度准则在p‒q平面上的图形为抛物线，故称这种用于描述冻土的非线性强度准则为抛物线强度准则。其中，c、φ、p_m与温度T的函数关系如下：

\left\{\begin{array}{l}c=-\mathrm{0.147}T+\mathrm{2.255},\\ \phi =-\mathrm{0.607}T+\mathrm{27.27},\\ p_{\mathrm{m}}=\sqrt[]{-\mathrm{65.79}T+\mathrm{45.61}}\end{array}\right.（12）

式中，T的单位取℃。由式（12）可以看出：在温度较高的情况下（0 ℃左右），p_m对温度的变化最敏感，其次是内摩擦角，而黏聚力对温度变化最不敏感；随着温度进一步降低，p_m对温度的敏感性降低。

对于与砂土性质差别较大的土体（如粉质黏土），Wang等^[24]进行了一系列复杂应力路径的三轴试验，其结果显示冻结黏土的强度也具有非线性，但并不符合上述冻结砂土体现出的抛物线强度准则。

由于上述抛物线型屈服准则是通过三轴试验确定的，因此p=(σ₁+2σ₃)/3，q=σ₁-σ₃。而在α较小、p_a >p_b >0的条件下，由式（8）～（10）确定的弹性状态下土体各应力的大小关系始终为σ_r =σ₁>σ_z =σ₂>σ_θ =σ₃。因此在进行后续推导时，该屈服准则可以进一步写为：

({\sigma }_r-{\sigma }_{\theta })=c+\frac{k}{\mathrm{3}}({\sigma }_r+\mathrm{2}{\sigma }_{\theta })-\frac{k}{\mathrm{18}{p}_{\mathrm{m}}}({\sigma }_r+\mathrm{2}{\sigma }_{\theta }{)}^{\mathrm{2}}（13）

即认为试验得出的抛物线型屈服准则仅为σ_r 、σ_θ 的函数。

4.2 塑性区的确定

随着内压的进一步增大，当土体中某些区域的σ_r 与σ_θ 的组合达到式（13）所对应的临界值时，土体进入塑性变形状态形成塑性区，如图2所示。处于塑性变形状态的土体不再满足热弹性的本构方程，而是满足式（13）的抛物线型准则。

因此，通过式（13）可以确定在内压p_a 作用下的塑性区范围。当p_a >p_b 时，塑性区总是从柱孔处r=a开始发展。为简化记号，塑性区的范围可以表示为r\les，其中s为塑性区与弹性区的分界面所对应的半径。在柱孔的尺寸、材料参数、温度分布保持不变的情况下，s由内压p_a 唯一确定，即s=f(p_a )。

4.3 弹塑性状态下的应力

塑性区的应力除了满足抛物线型屈服准则（13），还应该满足平衡方程（2），这是因为应力的平衡方程适用于任意连续介质。将应力的平衡方程和抛物线型屈服准则联立，解得基于抛物线强度准则的塑性区应力分布为：

\frac{Y}{r}=\frac{{\left|k\left(\frac{{\sigma }_r+\mathrm{2}{\sigma }_{\theta }}{\mathrm{3}}-p_{\mathrm{m}}\right)+M\right|}^{\left(\frac{p_{\mathrm{m}}}{M}+\frac{\mathrm{2}}{\mathrm{3}}\right)}}{{\left|k\left(\frac{{\sigma }_r+\mathrm{2}{\sigma }_{\theta }}{\mathrm{3}}-p_{\mathrm{m}}\right)-M\right|}^{\left(\frac{p_{\mathrm{m}}}{M}-\frac{\mathrm{2}}{\mathrm{3}}\right)}}（14）

式中，Y为积分常数，M为已知材料常数。M计算如下：

M=\sqrt[]{k^{\mathrm{2}}{p}_{\mathrm{m}}^{\mathrm{2}}+\mathrm{2}k{p}_{\mathrm{m}}c}（15）

积分常数Y可由内壁r=a的应力边界条件确定：

{\sigma }_r\left|{}_{r=a}\right.=p_a（16）

{\sigma }_{\theta }\left|{}_{r=a}\right.=\psi (p_a)（17）

式中，ψ(p_a )为环向应力σ_θ 边界条件，将式（16）代入式（13）即可求得。

当给定了内压p_a 后，即可根据式（13）与（14）求出塑性区的σ_r 、σ_θ =ψ(σ_r )。由于塑性区的应力只由式（2）和（13）两个方程确定，因此无法直接求出竖向应力σ_z。由平面应变的条件，在弹性条件下，σ_z 来源于竖向应变ε_z 的约束，但冻土进入塑性变形后产生屈服，塑性区的进一步发展使得竖向约束变弱。因此可以认为塑性区σ_z 的值等于弹性区与塑性区的分界面处σ_z 的值。

对于弹性区的应力分布，根据分界面的应力边界条件有：

{\sigma }_r\left|{}_{r=s}\right.=p_s（18）

式中，p_s 为通过式（13）、（14）计算得出的r=s处的塑性区径向应力。进而可以将弹性区视为内径为s、外径为b的柱孔扩张问题，用弹性状态下的分析方式求出σ_r 、σ_θ 、σ_z。

4.4 弹塑性状态下的位移

弹性区的位移u(r,s)保持与式（7）相同的形式，仅积分常数发生变化，且依赖于塑性区发展的程度，即：

u=\frac{\alpha (\mathrm{1}+\nu )}{\mathrm{1}-\nu }\left(C_{\mathrm{1}}^{\text{'}}r+\frac{C_{\mathrm{2}}^{\text{'}}}{r}+\frac{T(s)}{\mathrm{2}\mathrm{l}\mathrm{n}\frac{b}{s}}r\left(\mathrm{l}\mathrm{n}\frac{b}{r}+\frac{\mathrm{1}}{\mathrm{2}}\right)\right)（19）

式中，C_{\mathrm{1}}^{\text{'}}、C_{\mathrm{2}}^{\text{'}}为积分常数。采用与第3节相同的分析方法可求出弹性区的位移。

由于抛物线强度准则是基于试验数据给出的，因此无法通过分析f(p,q,T)求出塑性区的位移。可以根据连续介质的运动方程，假定塑性区的压缩方程不变，从而给出塑性区位移的理论解^[8]。

在对塑性区的土体进行分析时，沿用连续介质力学中的一般规定，区分物质坐标 R 与空间坐标 r，即 R 代表某一物质颗粒， r 代表某一空间点，并且 R = r ( R, t₀)， r = r ( R,t)。而在弹性分析时采用的是小变形假定，故通常认为土体在受到外部荷载产生位移后，也不会距离初始位置太远，即 r ≈ R。考察塑性区土体中任意一个颗粒 R 的位移，可以认为其空间坐标仅与颗粒的物质坐标 R 和塑性区外边界s有关，如图3所示。

由于柱孔扩张模型具有对称性，向量 R、r 均只有一个分量，因此在后续分析中使用R、r代替 R 、 r。故对于柱孔扩张问题可以将塑性区的总位移写为：

u=u(R,s)=r(R,s)-R（20）

土体中某一颗粒R的速度为位置坐标r(R,t)对时间的1阶导数，同时也是位移的物质导数。如果形式上认为s与时间相关，则上述关系可写为：

\frac{Du}{Dt}={\left(\frac{\partial r}{\partial t}\right)}_R=\frac{\partial r}{\partial s}\dot{s}（21）

式中，D为物质导数求导运算，下标R表示对时间t求导时固定物质坐标R不变，\dot{s}为s对时间求导数。

位移的物质导数为：

\frac{Du}{Dt}=\frac{\partial u}{\partial t}+\frac{\partial r}{\partial t}·\frac{\partial u}{\partial r}=\left(\frac{\partial u}{\partial s}+\frac{\partial r}{\partial s}·\frac{\partial u}{\partial r}\right)\dot{s}（22）

联立式（21）、（22）可得：

V(R,s)=\frac{\partial r}{\partial s}=\frac{\frac{\partial u}{\partial s}}{\mathrm{1}-\frac{\partial u}{\partial r}}（23）

式中，V为物质颗粒的速度。

颗粒的应变速率{\dot{\epsilon }}_r、{\dot{\epsilon }}_{\theta }与应力速率{\dot{\sigma }}_r、{\dot{\sigma }}_{\theta }，实质上是各个分量的物质导数，即：

{\dot{\epsilon }}_r=\frac{D}{Dt}\left(\frac{\partial u}{\partial R}\right)=\frac{\partial }{\partial R}{\left(\frac{\partial u}{\partial t}\right)}_R=\frac{\partial }{\partial R}\left(\frac{Du}{Dt}\right)=\frac{\partial }{\partial R}{\left(\frac{\partial r}{\partial t}\right)}_R（24）

式（24）可进一步写为：

{\dot{\epsilon }}_r=\frac{\partial V}{\partial R}\dot{s}（25）

同理可得：

{\dot{\epsilon }}_{\theta }=\frac{D}{Dt}\left(\frac{u}{R}\right)=\frac{V}{R}\dot{s}（26）

对于应力速率{\dot{\sigma }}_r，根据物质导数的求导法则，有：

{\dot{\sigma }}_r=\left(\frac{\partial {\sigma }_r}{\partial s}+V\frac{\partial {\sigma }_r}{\partial r}\right)\dot{s}（27）

同理，对于{\dot{\sigma }}_{\theta }，有：

{\dot{\sigma }}_{\theta }=\left(\frac{\partial {\sigma }_{\theta }}{\partial s}+V\frac{\partial {\sigma }_{\theta }}{\partial r}\right)\dot{s}（28）

在塑性区，假定应变速率和应力速率与弹性阶段保持同样的形式，即塑性区的压缩方程为：

{\dot{\epsilon }}_{ii}=\frac{(\mathrm{1}-\mathrm{2}\nu )}{E}{\dot{\sigma }}_{ii}（29）

式中，{\dot{\epsilon }}_{ii}为体积应变速率，{\dot{\sigma }}_{ii}为体积应力速率。注意到压缩模量K=E/(3(1-2ν))，则对于平面应变问题，式（29）的分量形式可写为：

\frac{\partial V}{\partial R}+\frac{V}{R}=\left(\frac{\partial }{\partial s}+\mu V\frac{\partial }{\partial R}\right)\left(\frac{{\sigma }_{\mathrm{m}}}{K}+\mathrm{3}\alpha T\right)（30）

式中：σ_m为平均应力，σ_m=(σ_r +σ_θ +σ_z )/3；μ为表征塑性区发展程度的参数，μ=(\partialR/\partialr)。

由于塑性区内r≠R，因此在上式中引入了μ。后续推导中暂定μ=1，将在第6节的参数分析中进一步讨论。令P(R, s)=σ_m /K+3αT，其物理含义为塑性区的体积应变，由式（28）可解得V(R, s)为：

V(R,s)=\frac{{\mathrm{e}}^P}{R}\left({\int }_a^R\frac{\partial P}{\partial s}·\frac{R}{{\mathrm{e}}^P}\mathrm{d}R+g(s)\right)（31）

式中，g(s)为积分常数。再由V = (\partialr/\partials)进一步解得r为：

r(R,s)={\int }_a^s\mathrm{d}s{\int }_a^R\frac{\partial P}{\partial s}·\frac{R}{{\mathrm{e}}^P}\mathrm{d}R+{\int }_a^s\frac{g(s){\mathrm{e}}^P}{R}\mathrm{d}s+C（32）

式中，C为积分常数。对于柱孔扩张问题，内壁R=a处的位移对桩基的径向位移起决定性作用，内壁的位移可写为：

u(a,s)=r(a,s)-a={\int }_a^s\frac{g(s){\mathrm{e}}^P}{a}\mathrm{d}s+C-a（33）

作为对比，此处给出基于Tresca强度准则（σ₁-σ₃=2c）的塑性区位移控制方程（以拉为正）^[8]为：

\begin{array}{c}\frac{\mathrm{d}r(a,s)}{\mathrm{d}s}=\left((m-\frac{m}{\mathrm{1}-m})\frac{s^{\mathrm{2}}}{b^{\mathrm{2}}}+\frac{m}{\mathrm{1}-m}+\frac{\mathrm{2}c(\mathrm{1}+\nu )}{E}\right)·\\ {\left(\frac{r}{s}\right)}^{\mathrm{2}m-\mathrm{1}}+\frac{m}{\mathrm{1}-m}(\frac{s^{\mathrm{2}}}{b^{\mathrm{2}}}-\mathrm{1})\frac{r}{s}\end{array}（34）

式中，m为与土的性质有关的材料常数，m=2c(1+ν)(1-2ν)/E。

5 算例与分析

本模型的工程应用可通过下述算例体现：考虑多年冻土区（环境温度T_e=-15 ℃）的桩基础工程，桩径a=1.0 m；冻土容重γ=20 kN/m³，弹性模量E=25 MPa，泊松比ν=0.25；施工时孔内土体由于扰动升温至T(a)=-10 ℃。内压p_a 从0逐渐增加至R\le2.0 m范围内的土体全部屈服。下面考察1.0 m<R<5.0 m（5倍桩径）范围内土体的应力和位移。

5.1 外压的确定

外压通过静止侧压力系数K₀计算得出。当温度T>-10 ℃时，K₀随温度的变化较大，但在低温冻土（T<-10 ℃）的条件下，静止侧压力系数K₀趋于常数。这可能是由于在低温条件下，冻土各相的物理性质趋于稳定。根据姚晓亮等^[25]对冻结粉土和细砂土的研究结果，在竖向应力较小（<2 MPa）的情况下，可取K₀=0.8。

本算例中考虑深度z=5.0 m，故p_b =0.8×20×5=80 kPa。事实上，根据后续计算结果，土体达到屈服时的内压p_a 远大于外压p_b，且土体的变形主要集中在柱孔附近，因此在自重应力的条件下外压的影响不大。

5.2 热膨胀系数的确定

尚未有文献对冻土的热膨胀系数（此处的定义为线膨胀系数而非体积膨胀系数）进行研究，算例中根据冰的热膨胀系数α_ice=5×10^-5 K^-1，假设土颗粒的膨胀远小于冰，根据土颗粒与冰的夹杂关系加权平均，可取α=3×10^-5 K^-1。上述所有参数的取值汇总于表1。

5.3 温度场的确定

将温度边界条件代入式（5）计算冻土的温度场，计算结果如图4所示。

5.4 弹性变形状态下的应力与位移

通过柱孔内、外的应力边界条件求解积分常数C₃、C₄，以内压p_a =1 000 kPa为例，此时周围土体均处于弹性变形状态，其应力场σ和位移场u分布如图5所示。

5.5 弹塑性变形状态下的应力

首先，需要确定塑性区的边界。根据弹性变形状态的分析结果，在周围土体中搜索满足式（13）的区域，从而得到塑性区外边界s与内压p_a 的对应关系，如图6所示，仅用二次多项式就可以得到较好的拟合结果。

在搜索周围土体的塑性区时发现：柱孔附近的土体最先产生屈服，若要土体产生进一步屈服需要很大的内压，而此时冻土的整体结构早已发生破坏，不再满足弹塑性的假定。因此，在土体的整体结构不发生破坏的情况下，强度准则f中的温度近似认为是孔内温度T(a)。此时，各强度参数的取值由式（12）和（15）计算得出，结果见表2。

土体中r>s的区域为弹性变形状态。以p_a =3 000 kPa为例，距桩中心距离r>1.35 m的区域转化为内径为s=1.35 m、外径b=5.0 m的弹性问题。即弹性方程的通解不变，新的积分常数由边界条件式（18）确定。

r<s的区域为塑性变形状态，将边界条件式（16）和式（17）代入（14），解得积分常数Y，进而通过式（14）得到塑性区的应力，如图7所示。图8、9分别为内压为5 000 、6 306 kPa时的应力分布。

从图7～9中可以看出，当土体进入塑性变形状态后环向应力σ_θ 会逐渐由拉应力（-）调整为压应力（+），并且在塑性边界面上出现不连续的现象。从土体的破坏机理来看，可能是由于拉应力过大导致土体产生了拉应变，塑性区应力重分布导致环向应力不连续。

5.6 柱孔内侧的位移

柱孔R=a处的位移可以避开二重积分的计算，即只需要考察函数g(s)的值。以p_a =3 MPa为例计算柱孔的位移。

观察式（33），利用s=a时柱孔刚好达到塑性的条件，容易求出u(1,1)=C-a=-0.083 m。若要计算u(1, 1.35)，则至少需要求出g(1.35)，才能通过数值积分得到近似值。由式（31）可知，只要计算出V(1.35, 1.35)的值，即可通过数值积分得到g(1.35)的值。而V(1.35, 1.35)可以通过式（23）计算。因为：

\mathrm{d}u=\frac{\partial u}{\partial s}\mathrm{d}s+\frac{\partial u}{\partial R}\mathrm{d}R（35）

利用上式的差分形式，可以得到：

u(\mathrm{1.36,1.36})-u(\mathrm{1.35,1.35})=\left(\frac{\partial u}{\partial s}+{\epsilon }_r\right)\times \mathrm{0.01}（36）

式（36）中仅有\partialu/\partials是未知量，其他项均可以通过弹性状态下的位移和应变求出具体值。从而计算得出V (1.35, 1.35)=-0.189 2。

接下来计算式（30）中的积分项。由第5.5节中对塑性区应力的计算方法，逐渐增大内压p_a 并计算此时的塑性区边界s与塑性区体积应变P，进而计算\partialP/\partials的值。采用线性近似计算导数，再用梯形公式进行数值积分，代入式（30）计算得出g(1.35) = -0.275 8。采用相同的步骤可以分别计算得到g(1.00) = -0.004 83，g(1.57) = -0.360 8，g(1.77) = -0.437 2。最后通过式（33）计算得出随内压分别增加至3、4、5 MPa时内侧的位移变化如图10所示。

作为对比，根据算例参数和式（34）进行计算，得到基于Tresca准则计算的柱孔内侧的位移计算结果（图10）。可以看出，式（34）在高内压下的位移计算结果明显偏大，不适用于冻土。粒子有限元法（PFEM）是一种适用于模拟岩土材料大变形问题的数值模拟方法^[26‒28]，Hauser等^[29]根据该方法对黏土的柱孔扩张问题进行了数值计算，并在计算中综合考虑了大变形和CASM模型（clay and sand model）^[30]。对于超固结黏土（超固结比为16）排水条件下的柱孔扩张问题，根据归一化参数进行换算，若将先期固结压力视为算例中柱孔处刚好出现塑性变形时的压力，内压和柱孔内侧位移的情况如图10所示，用粒子有限元法进行数值计算的结果与本文的理论计算结果较为接近。

6 参数分析

6.1 孔内温度的影响

由于冻土的热膨胀系数α很小，柱孔内外的温度差对弹性状态下的应力场和位移场影响较小。受温度影响从而变化较大的力学指标主要为式（12）中的c、φ和p_m。在算例中已经指出，由于塑性区仅分布在距离柱孔很近的范围内，因此上述力学指标中的温度T= T(a)。对于不同的孔内温度T(a)，可以得到塑性区边界s曲线随温度变化如图11所示。

从图11可以看出，温度的升高导致冻土强度的降低。当内压相同时，塑性区会随着孔内温度的升高而向外扩展。根据图7～9的应力分布曲线可知，塑性区与弹性区的环向应力场σ_θ 有很大的不同，故塑性区范围的改变将会对σ_θ 的分布产生较大影响，如图12所示。在塑性区，孔内温度越高，对应的环向应力越小；而在弹性区，更高的温度则对应了更大的环向应力。此外，随着温度的升高，环向应力在塑性区边界面上的突变现象越来越明显。

6.2 塑性区发展系数的影响

在计算塑性区位移u时，引入参数μ=\partialR/\partialr来表示塑性变形与弹性小变形假定偏离的程度。如果简单地认为该参数是线性的，则有μ=R/r。将其代入塑性位移计算的表达式中，则柱孔内侧实际的塑性位移u_p (a,s)为：

u_p(a,s)=\frac{a}{a+u(a,s)}u(a,s)（37）

由式（37）可以看出，实际的塑性位移小于计算得出的位移，说明在塑性变形状态下，冻土仍有压硬性。

7 结 论

本文提出了基于抛物线强度准则的冻土热弹性‒塑性模型，并且针对冻土地区桩基础施工的工况进行了算例分析。主要结论如下：

1）对于非均匀温度场作用下的冻土，需要考虑温度对冻土位移、应力状态、强度参数的影响。因此本文提出的模型结合了热弹性本构方程、抛物线强度准则，将温度作为变量引入到了柱孔扩张模型中，进而得到了基于抛物线强度准则的热弹性‒塑性柱孔扩张模型。

2）随着内压的进一步增大，周围土体逐渐从弹性状态转变为弹塑性状态，并形成从柱孔内侧向外扩展的塑性区。根据抛物线强度准则，可以将整个计算区域分为塑性区和弹性区两部分，分别采用不同的方程进行求解。对于塑性区的应力，本文根据平衡方程和抛物线型准则给出了解析式；对于塑性区的位移，本文结合连续介质力学的运动方程进行了分析。

3）在冻土地区桩基础施工的过程中，桩身周围的土体处于非均匀温度场的作用下，其应力与位移的求解可视为温度非均匀分布情况下的柱孔扩张问题。本文据此进行了算例分析，得到了不同内压下周围土体的应力分布曲线，以及孔壁内侧的位移曲线。

4）随着孔内温度的升高，塑性区逐渐向外扩展。对于塑性区，孔内温度的升高导致环向应力减小；而对于弹性区，孔内温度的升高导致环向应力的增大，并且环向应力在塑性区边界面上的突变现象越来越明显。塑性变形与小变形假定偏离的程度可以用塑性区发展系数表示，其值随着塑性区范围的增大而减小，由此计算得出的实际塑性位移小于小变形假定下的塑性位移。

参考文献

[1] 徐学祖,王家澄,张立新.冻土物理学[M].北京:科学出版社,2001.

[2] Obu J.How much of the Earth's surface is underlain by permafrost?[J].Journal of Geophysical Research:Earth Surface,2021,126(5):e2021JF006123. doi:10.1029/2021jf006123

[3] Hjort J, Streletskiy D, Doré G,et al.Impacts of permafrost degradation on infrastructure[J].Nature Reviews Earth & Environment,2022,3(1):24‒38. doi:10.1038/s43017-021-00247-8

[4] De Guzman E M B, Stafford D, Alfaro M C,et al.Large-scale direct shear testing of compacted frozen soil under freezing and thawing conditions[J].Cold Regions Science and Technology,2018,151:138‒147. doi:10.1016/j.coldregions.2018.03.011

[5] 马巍,王大雁,等.冻土力学[M].北京:科学出版社,2014.

[6] Li Bo, Wang Tieliang, Luan Maotian.Incremental elastic-plastic solutions for cone tip resistance in sand[J].Journal of Wuhan University of Technology,2010,32(6):6‒10.

[7] 李波,王铁良,栾茂田.砂土静力触探锥头阻力增量弹塑性解[J].武汉理工大学学报,2010,32(6):6‒10.

[8] Han Tongchun, Dou Hongqiang.Application of space axisymmetric solution of cylindrical hole expansion theory in pile driving and soil squeezing effect[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(S1):3209‒3215. doi:10.3969/j.issn.1000-6915.2012.z1.083

[9] 韩同春,豆红强.柱孔扩张理论的空间轴对称解在沉桩挤土效应中的应用[J].岩石力学与工程学报,2012,31(S1):3209‒3215. doi:10.3969/j.issn.1000-6915.2012.z1.083

[10] Yu Haisui.Cavity Expansion Methods in Geomechanics[M].Dordrecht:Kluwer Academic,2000. doi:10.1007/978-94-015-9596-4

[11] Zeng Yingjun, Yang Min, Sun Qing.Elastic-plastic solution for cylinder cavity contraction and its application to tunnel engineering[J].Journal of Tongji University(Natural Science),2012,40(10):1480‒1485. doi:10.3969/j.issn.0253-374x.2012.10.008

[12] 曾英俊,杨敏,孙庆.圆柱孔收缩的弹塑性解及其在隧道工程中的应用[J].同济大学学报(自然科学版),2012,40(10):1480‒1485. doi:10.3969/j.issn.0253-374x.2012.10.008

[13] Hill R.The Mathematical Theory of Plasticity[M].Oxford:Clarendon Press,1998. doi:10.1093/oso/9780198503675.001.0001

[14] Vesić A S.Expansion of cavities in infinite soil mass[J].Journal of the Soil Mechanics and Foundations Division,1972,98(3):265‒290. doi:10.1061/jsfeaq.0001740

[15] Carter J P, Booker J R, Yeung S K.Cavity expansion in cohesive frictional soils[J].Géotechnique,1986,36(3):349‒358. doi:10.1680/geot.1986.36.3.349

[16] Jiang Mingjing, Shen Zhujiang.Expansion of cylindrical cavity of materials with strain-softening behaviour[J].Chinese Journal of Geotechnical Engineering,1995(4):10‒19.

[17] 蒋明镜,沈珠江.考虑材料应变软化的柱形孔扩张问题[J].岩土工程学报,1995(4):10‒19.

[18] Chen S L, Liu K.Undrained cylindrical cavity expansion in anisotropic critical state soils[J].Géotechnique,2019,69(3):189‒202. doi:10.1680/jgeot.16.p.335

[19] Mo Pinqiang, Yu Haisui.Drained cavity expansion analysis with a unified state parameter model for clay and sand[J].Canadian Geotechnical Journal,2018,55(7):1029‒1040. doi:10.1139/cgj-2016-0695

[20] Rosenberg Z, Dekel E.Analytical solution of the spherical cavity expansion process[J].International Journal of Impact Engineering,2009,36(2):193‒198. doi:10.1016/j.ijimpeng.2007.12.014

[21] Zhang Jiaqi, Zhao Chunfeng, Zhao Cheng,et al.Cylindrical and spherical cavity reverse expansion considering elastoplastic unloading[J].Rock and Soil Mechanics,2023,44(11):3224‒3234.

[22] 张家奇,赵春风,赵程,等.考虑弹塑性卸荷的柱孔和球孔反向扩张解[J].岩土力学,2023,44(11):3224‒3234.

[23] Wang Jiaeheng, Wang Shanling, Qiu Guoqing.Permafrost along the Qinghai—Xizang highway[J].Acta Geographica Sinica,1979,34(1):18‒32.

[24] 王家澄,王绍令,邱国庆.青藏公路沿线的多年冻土[J].地理学报,1979,34(1):18‒32.

[25] 徐秉业.应用弹塑性力学[M].北京:清华大学出版社,1995.

[26] 吴家龙.弹性力学[M].3版.北京:高等教育出版社,2016.

[27] 刘恩龙,刘星炎,张革.高寒地区冻结掺合土料的力学特性与本构模型[M].北京:中国水利水电出版社,2021.

[28] Liu Xingyan, Liu Enlong, Zhang De,et al.Study on strength criterion for frozen soil[J].Cold Regions Science and Technology,2019,161:1‒20. doi:10.1016/j.coldregions.2019.02.009

[29] Ma Wei, Wu Ziwang, Zhang Changqing.Strength and yield criterion of frozen soil[J].Journal of Glaciology and Geocryology,1993,15(1):129‒133.

[30] 马巍,吴紫汪,张长庆.冻土的强度与屈服准则[J].冰川冻土,1993,15(1):129‒133.

[31] Wang Baoxuan, Xu Xiangtian, Wang Xiaoshan,et al.Mechanical behavior and strength criterion of frozen silty clay under complex stress paths[J].Geoderma,2023,435:116506. doi:10.1016/j.geoderma.2023.116506

[32] Yao Xiaoliang, Qi Jilin, Yu Fan.The experimental study on coefficient of lateral earth pressure at rest for frozen soil[J].Chinese Journal of Underground Space and Engineering,2011,7(6):1108‒1113. doi:10.1088/0253-6102/55/2/31

[33] 姚晓亮,齐吉琳,余帆.冻土静止侧压力系数的试验研究[J].地下空间与工程学报,2011,7(6):1108‒1113. doi:10.1088/0253-6102/55/2/31

[34] Oñate E, Owen R.Particle-Based Methods:Fundamentals and Applications[M].Dordrecht:Springer Netherlands,2011. doi:10.1007/978-94-007-0735-1

[35] Monforte L, Carbonell J M, Arroyo M,et al.Performance of mixed formulations for the particle finite element method in soil mechanics problems[J].Computational Particle Mechanics,2017,4(3):269‒284. doi:10.1007/s40571-016-0145-0

[36] Monforte L, Arroyo M, Carbonell J M,et al.Coupled effective stress analysis of insertion problems in geotechnics with the Particle Finite Element Method[J].Computers and Geotechnics,2018,101:114‒129. doi:10.1016/j.compgeo.2018.04.002

[37] Hauser L, Schweiger H F.Simulation of cavity expansion with the clay and sand model using G-PFEM[C]//Challenges and Innovations in Geomechanics.Cham:Springer,2021:536‒543. doi:10.1007/978-3-030-64518-2_63

[38] Yu H S.CASM:A unified state parameter model for clay and sand[J].International Journal for Numerical and Analytical Methods in Geomechanics,1998,22(8):621‒653. doi:10.1002/(sici)1096-9853(199808)22:8<621::aid-nag937>3.3.co;2-#

基金资助

国家自然科学基金项目(U22A20596)