三维化方法对土的平面应变强度和变形的影响分析

姚仰平; 武孝天; 崔文杰

doi:10.11779/CJGE20211389

三维化方法对土的平面应变强度和变形的影响分析 English Version

北京航空航天大学交通科学与工程学院, 北京 100083

基金项目:

国家重点研发计划项目 2018YFE0207100

国家自然科学基金项目 51979001

详细信息

作者简介:
姚仰平(1960—)，男，教授，博士生导师，“973”首席科学家，主要从事岩土本构理论研究工作。E-mail: ypyao@buaa.edu.cn

通讯作者:
崔文杰, E-mail: wcui21@buaa.edu.cn

中图分类号: TU473
计量
- 文章访问数: 431
- HTML全文浏览量: 133
- PDF下载量: 129
出版历程
- 收稿日期: 2021-12-24
- 网络出版日期: 2023-03-15
- 刊出日期: 2023-02-28

Influences of 3D model generalization approach on calculation of stress and strain of soils under plane strain

School of Transportation Science and Engineering, Beihang University, Beijing 100083, China

摘要

摘要: 临界状态本构模型需要通过三维化来反映材料在三维应力状态下的力学特性。采用不同三维化方法在描述应力罗德角和偏应力对应力应变的影响程度上具有差异性，在平面应变条件下，这种差异性对土的应力应变计算结果影响往往非常显著。以K₀固结土的统一硬化（UH）本构模型和空间滑动面（SMP）强度准则为例，分别采用变换应力（TS）方法和g(θ)方法这两种常用方法对模型三维化，推导了相应的三维弹塑性刚度矩阵[D_ep]，阐述了不同方法在三维应力应变计算过程上的本质区别。相比g(θ)方法，TS方法可以合理描述不同K₀状态下土的应力水平对偏平面上屈服曲线形状的影响规律，即由低应力比下的近圆形转变为剪切破坏时的SMP破坏准则形状，实现了土体从剪切屈服到剪切破坏的一致性。通过对一系列平面应变单元试验和边值问题的分析计算结果表明，采用TS方法三维化的UH模型预测结果与已有试验规律更为吻合。在不同K₀固结条件下，由于计算得到的破坏面呈现不规则形状，采用g(θ)方法往往过高或过低地估计平面应变条件下土体的应力水平。
- 三维化方法 /
- 本构模型 /
- 平面应变 /
- 土的应力应变
Abstract: To describe the shear yield and failure behavior of soils in the generalized 3D stress space, the critical state constitutive model needs to be combined with some strength criteria. The use of various 3D model generalization approaches is frequently highly important in characterizing the impacts of stress Lode angle and deviatoric stress on stress and strain, especially under plane strain. The unified hardening (UH) constitutive model for K₀ consolidated soils and the spatially mobilized plane (SMP) strength criterion are used as examples, and the transform stress (TS) and g(θ) methods are adopted for the 3D model generalization. The 3D elastoplastic stiffness matrix [D_ep] associated with each model generalization approach is deduced, and the key difference between various approaches for calculating 3D stress and strain is discussed. The TS method, in comparison to the g(θ) method, is able to transform the yield curve on the π plane from a nearly circular shape under a low stress ratio to the shape of the SMP strength criterion at failure, which demonstrates the consistency from yield to failure. Predictions of a series of plane strain element tests and boundary value problems are performed with the 3D UH model generalized by the TS and g(θ) methods. The results show that the predicted results from the TS method are more consistent with the existing experimental measurements. Due to the irregular shape of the failure surface when using the g(θ) method, the estimated stress level is often either too high or too low under plane strain for K₀-consolidated soils.
- 3D model generalization approach /
- constitutive model /
- plane strain /
- stress and strain of soil

HTML全文

0. 引言

在岩土工程数值分析中，采用弹塑性本构模型描述土的应力应变关系需要假定屈服面和塑性势面在偏平面上的形状，将基于三轴试验规律建立的本构关系扩展到三维应力空间中，实现本构模型的三维化^[1-4]。本构模型通过三维化可以描述罗德角的变化对土的应力应变关系和临界破坏强度的影响。大多数临界状态模型的建立都考虑了模型的三维化^[5-8]，目前已有的方法主要包括g(θ)方法^{[3, 9-11]}、M(θ)方法^[12]、t_ij方法^{[9, 13]}和变换应力（TS）方法^{[8, 14-17]}等。

采用不同的方法对本构模型三维化，计算得到的应力应变结果可能存在较大的差异性。孙德安等^[18]采用变换应力方法将空间滑动面（SMP）准则与K₀剑桥模型相结合，比较了模型三维化前后计算三轴压缩、三轴拉伸以及平面应变等问题上的差异。Sun等^[19]采用TS方法将莫尔库仑准则与修正剑桥（MCC）模型结合，比较了三维化前后计算三轴拉伸和真三轴等问题上的差异。Yao等^[14]基于TS方法将拉德准则与统一硬化模型结合，分析了土的三维强度和应力应变关系。Yao等^[20]基于g(θ)和TS方法将SMP准则与K₀-MCC结合，分析了土在三轴拉伸下的应力应变差异。然而，当前大多研究往往只针对某一特定的方法，将模型三维化后与未三维化模型在计算上的对比分析，但不同三维化方法之间的对比分析较少。此外，当前对采用不同的三维化方法影响的对比研究大多局限于对常规三轴或者真三轴等室内试验的预测分析，计算过程中应力罗德角往往保持不变，对罗德角发生变化的平面应变问题的影响分析较少。实际上，由于平面应变加载过程中土体的罗德角在不断变化，因此可以很好地反映出不同的三维化方法对应力应变计算时结果的影响程度。

本文以K₀固结土统一硬化（UH）本构模型和SMP强度准则为例，采用目前数值分析中最常用的两种三维化方法，即TS和g(θ)方法，实现模型的三维化。首先介绍了TS和g(θ)方法对K₀-UH模型实现三维化的过程，分析了不同方法下模型屈服面和破坏面演化规律的差异性及其对土的应力应变计算的影响；随后推导了三维化本构模型的弹塑性刚度矩阵[D_ep]，展示了二者在计算方法上的本质不同；最后分别以平面应变单元试验和边值问题分析为例，深入分析了不同应力路径下两种方法在计算分析中的差异，比较了两种方法的优劣。

1. 本构模型三维化方法

1.1 TS方法

TS方法是通过应力变换，使屈服面由真实应力空间在π平面上基于SMP准则的曲边三角形，转换为变换应力空间中的Mises圆，随后在变换应力空间将变换后的应力引入本构模型进行计算。该方法采用的变换公式如下^[16]：

${\tilde {\mathit{\boldsymbol{\sigma}}} _{ij}} = p{\mathit{\boldsymbol{\delta}} _{ij}} + \frac{{{q_{\text{c}}}}}{q}({\mathit{\boldsymbol{\sigma}} _{ij}} - p{\mathit{\boldsymbol{\delta}} _{ij}}) \text{，}$

(1)

式中： ${\mathit{\boldsymbol{\delta}} _{ij}}$ 为克朗内克符号； ${\tilde {\mathit{\boldsymbol{\sigma}}} _{ij}}$ 表示转换应力空间的应力；p为平均正应力；q为广义剪应力； ${q_{\text{c}}}$ 为土体在三轴压缩下的偏应力，计算式为^[14]

${q_{\text{c}}} = \frac{{2{I_1}}}{{3\sqrt {({I_1}{I_2} - {I_3})/({I_1}{I_2} - 9{I_3})} - 1}} 。$

(2)

式中：I₁，I₂和I₃分别为第一、二和三应力不变量。

在变换应力空间中，K₀-UH模型的屈服面表示为^[17]

$\tilde f = \tilde g = \ln \left( {\frac{{\tilde p}}{{{{\tilde p}_{x0}}}}} \right) + \ln \left[ {1 + {{\left( {\frac{{{{\tilde \eta }^ * }}}{{{{\tilde M}^ * }}}} \right)}^2}} \right] - \frac{{\tilde H}}{{{c_{\text{p}}}}} = 0 \text{，}$

(3)

其中， ${\tilde \eta ^ * }$ 和 ${\tilde M^ * }$ 分别为

${\tilde \eta ^ * } = \sqrt {3({{\tilde \eta }_{ij}} - {{\tilde \eta }_{ij0}})({{\tilde \eta }_{ij}} - {{\tilde \eta }_{ij0}})/2} \text{，}$

(4)

${\tilde M^ * } = \sqrt {{M^2} - \tilde \eta _0^2} \text{，}$

(5)

$\tilde H$ 为统一硬化参量，表达为

$\tilde H = \int {\frac{{\tilde M_f^4 - \eta _{\text{k}}^4}}{{{{\tilde M}^4} - \tilde \eta _{\text{k}}^4}}} {{{\rm{d}}}{\rm{\varepsilon }}}_{\text{v}}^{\text{p}} \text{，}$

(6)

${\tilde \eta _k} = \sqrt {3\left| {{{\tilde \eta }_{ij}}{{\tilde \eta }_{ij}} - {{\tilde \eta }_{ij0{\text{k}}}}{{\tilde \eta }_{ij0{\text{k}}}}} \right|/2} 。$

(7)

式中： ${\tilde \eta _{ij}} = ({\tilde \sigma _{ij}} - \tilde p{\delta _{ij}})/\tilde p$ ；c_p=( $\lambda$ -κ)/(1+e₀)，e₀为土体的初始孔隙比； $\lambda$ 和κ分别为e-lnp平面上等向固结加载曲线和卸载-再加载曲线的斜率；M为三轴压缩临界状态应力比。采用TS方法三维化的K₀-UH模型屈服面在一般应力空间中的形状如图 1所示。

图 1 采用TS方法三维化的K₀-UH屈服面

Figure 1. Three-dimensional yield surface of K₀-UH model using TS method

下载: 全尺寸图片幻灯片

1.2 g(θ)方法

g(θ)方法实现模型三维化时与TS方法保持一致，即真实应力空间中屈服面采用SMP准则，塑性势面保持Mises圆的形状。采用 $g(\theta )$ 方法三维化的K₀-UH模型的屈服面可表示为

$f = \ln \left( {\frac{p}{{{p_{x0}}}}} \right) + \ln \left[ {1 + {{\left( {\frac{{{\eta ^ * }}}{{{M^ * }g({\theta ^ * })}}} \right)}^2}} \right] - \frac{H}{{{c_{\text{p}}}}} = 0 \text{，}$

(8)

塑性势可表示为

$g = \ln \left( {\frac{p}{{{p_x}}}} \right) + \ln \left[ {1 + {{\left( {\frac{{{\eta ^ * }}}{{{M^ * }}}} \right)}^2}} \right] = 0 \text{，}$

(9)

式中：η^*和M^*的表达式分别为

${\eta ^ * } = \sqrt {3({\eta _{ij}} - {\eta _{ij0}})({\eta _{ij}} - {\eta _{ij0}})/2} \text{，}$

(10)

${M^ * } = \sqrt {{M^2} - \eta _0^2} 。$

(11)

基于SMP准则的 $g(\theta )$ 表达式^[4]可写为

$g({\theta ^ * }) = \frac{{\sqrt 3 (\sqrt {8 + {{\sin }^2}{\varphi _0}} - \sin {\varphi _0})}}{{4\sqrt {2 + {{\sin }^2}{\varphi _0}} \cos \psi }} \text{，}$

(12)

式中：SMP准则的材料参数ψ和θ^∗的表达式分别为

$\psi = \frac{1}{3}\arccos \left[ {{{\left( {\frac{3}{{2 + {{\sin }^2}{\varphi _0}}}} \right)}^{3/2}}\sin {\varphi _0}\sin 3{\theta ^ * }} \right] \text{，}$

(13)

${\theta ^ * } = - \frac{1}{3}\arcsin \left[ {\frac{{3\sqrt 3 }}{2}\frac{{J_3^ * }}{{{{(J_2^ * )}^{3/2}}}}} \right]{\text{ }}\left( { - \frac{{\text{π }}}{6} \leqslant {\theta ^ * } \leqslant \frac{{\text{π }}}{6}} \right) 。$

(14)

采用 $g(\theta )$ 方法三维化的K₀-UH模型屈服面在一般应力空间中的形状如图 2所示。

图 2 采用g(θ)方法三维化的K₀-UH屈服面

Figure 2. Three-dimensional yield surface of K₀-UH model using g(θ) method

下载: 全尺寸图片幻灯片

1.3 TS和 $g(\theta )$ 方法对三维屈服面及临界破坏强度的影响

图 3~6分别画出了K₀=1和K₀=0.625时采用TS和 $g(\theta )$ 方法时UH模型屈服面在π平面上从初始屈服到最终破坏的形状演化规律。从图 3和图 5可以看出，不论土体处于初始各向同性（K₀=1）还是不同K₀条件下的初始各向异性状态，TS方法均可以将屈服曲线在π平面上由低应力比下的近圆形变为剪切破坏时的SMP准则形状；从，可以看出，采用 $g(\theta )$ 三维化后的各向同性屈服面在加载过程中始终保持SMP破坏面形状不变，在初始应力各向异性条件下，屈服面最终变为内凹形状的破坏面。这种不规则的破坏面与已有试验规律相差较大，在计算过程中可能会错误地估计土的临界破坏强度。

图 3 采用TS方法的初始各向同性三维屈服面和破坏面（K₀=1）

Figure 3. Three-dimensional yield surface and failure surface using TS method with initial isotropic stress status (K₀=1)

下载: 全尺寸图片幻灯片

图 4 采用g(θ)方法的初始各向同性三维屈服面和破坏面（K₀=1）

Figure 4. Three-dimensional yield surface and failure surface using g(θ) method with initial isotropic stress status (K₀=1)

下载: 全尺寸图片幻灯片

图 5 采用TS方法的初始各向异性三维屈服面和破坏面（K₀=0.625）

Figure 5. Three-dimensional yield surface and failure surface using TS method with initial anisotropic stress status (K₀=0.625)

下载: 全尺寸图片幻灯片

图 6 采用g(θ)方法的初始各向异性三维屈服面和破坏面（K₀=0.625）

Figure 6. Three-dimensional yield surface and failure surface using g(θ) method with initial anisotropic stress status (K₀=0.625)

下载: 全尺寸图片幻灯片

2. 应力应变关系

2.1 计算方法

土体单元的弹塑性应力应变关系表示为

${\text{d}}{\sigma '_{ij}} = \mathit{\boldsymbol{D}}_{ijkl}^{{\text{ep}}}{\text{d}}{\varepsilon _{kl}} 。$

(15)

（1）采用TS方法三维化后，K₀-UH模型的弹塑性张量 $\mathit{\boldsymbol{D}}_{ijkl}^{{\text{ep}}}$ 计算式为

$\mathit{\boldsymbol{D}}_{ijkl}^{{\text{ep}}} = \mathit{\boldsymbol{D}}_{ijkl}^{\text{e}} - \frac{{\mathit{\boldsymbol{D}}_{ijmn}^{\text{e}}\frac{{\partial \tilde g}}{{\partial {{\tilde \sigma '}_{mn}}}}\frac{{\partial \tilde f}}{{\partial {{\sigma '}_{{\text{st}}}}}}\mathit{\boldsymbol{D}}_{stkl}^{\text{e}}}}{{ - \frac{{\partial \tilde f}}{{\partial \tilde H}}\frac{{\partial \tilde H}}{{\partial \varepsilon _{\text{v}}^{\text{p}}}}\frac{{\partial \tilde g}}{{\partial \tilde p'}} + \frac{{\partial \tilde f}}{{\partial {{\sigma '}_{uv}}}}\mathit{\boldsymbol{D}}_{uvwz}^{\text{e}}\frac{{\partial \tilde g}}{{\partial {{\tilde \sigma '}_{wz}}}}}} 。$

(16)

式（16）以矩阵形式表示为

$\left[{\mathit{\boldsymbol{D}}}_{\text{ep}}\right]=\left[\begin{array}{ccc}L+2G-\frac{{\alpha }_{r}{\beta }_{r}}{N}& L-\frac{{\alpha }_{r}{\beta }_{\theta }}{N}& L-\frac{{\alpha }_{r}{\beta }_{z}}{N}\\ L-\frac{{\alpha }_{\theta }{\beta }_{r}}{N}& L+2G-\frac{{\alpha }_{\theta }{\beta }_{\theta }}{N}& L-\frac{{\alpha }_{\theta }{\beta }_{z}}{N}\\ L-\frac{{\alpha }_{z}{\beta }_{r}}{N}& L-\frac{{\alpha }_{z}{\beta }_{\theta }}{N}& L+2G-\frac{{\alpha }_{z}{\beta }_{z}}{N}\end{array}\right]\text{，}$

(17)

式中， $L = K - 2G/3$ ，K和G分别为体积模量和剪切模量，表示为 $K=\frac{1+e}{\kappa }{p}^{\prime }，G=\frac{3(1-2\nu )}{2(1+\nu )}K$ 。其余各变量为

$\begin{array}{l}{\alpha }_{i}=L\left(\frac{\partial g}{\partial {{\tilde{\sigma }}^{\prime }}_{1}}+\frac{\partial g}{\partial {{\tilde{\sigma }}^{\prime }}_{2}}+\frac{\partial g}{\partial {{\tilde{\sigma }}^{\prime }}_{3}}\right)+2G\frac{\partial g}{\partial {{\tilde{\sigma }}^{\prime }}_{i}}\text{ }\text{，}\\ {\beta }_{i}=L\left(\frac{\partial f}{\partial {{\tilde{\sigma }}^{\prime }}_{r}}\frac{\partial {{\tilde{\sigma }}^{\prime }}_{r}}{\partial {{\sigma }^{\prime }}_{1}}+\frac{\partial f}{\partial {{\tilde{\sigma }}^{\prime }}_{s}}\frac{\partial {{\tilde{\sigma }}^{\prime }}_{s}}{\partial {{\sigma }^{\prime }}_{2}}+\frac{\partial f}{\partial {{\tilde{\sigma }}^{\prime }}_{t}}\frac{\partial {{\tilde{\sigma }}^{\prime }}_{t}}{\partial {{\sigma }^{\prime }}_{3}}\right)+2G\frac{\partial f}{\partial {\tilde{{\sigma }^{\prime }}}_{u}}\frac{\partial {{\tilde{\sigma }}^{\prime }}_{u}}{\partial {{\sigma }^{\prime }}_{i}}\text{，}\\ N=\frac{{\tilde{M}}_{f}^{4}-{\tilde{\eta }}_{\text{k}}^{4}}{{c}_{\text{p}}{\tilde{p}}^{\prime }({\tilde{M}}^{4}-{\tilde{\eta }}_{\text{k}}^{4})}\left(1-\frac{3\left({{\tilde{\eta }}^{\prime }}_{mn}{{\tilde{\eta }}^{\prime }}_{mn}-{{\tilde{\eta }}^{\prime }}_{mn}{{\tilde{\eta }}^{\prime }}_{mn0}\right)}{{\tilde{M}}^{\ast }{}^{2}+{\tilde{\eta }}^{\ast }{}^{2}}\right)+{\beta }_{i}\frac{\partial g}{\partial {\tilde{{\sigma }^{\prime }}}_{i}}。\end{array}\}$

(17a)

式中：

$\left. \begin{array}{l}\frac{\partial g}{\partial {{\tilde{\sigma }}^{\prime }}_{ij}}=\frac{\partial f}{\partial {{\tilde{\sigma }}^{\prime }}_{ij}}=\frac{1}{\tilde{{p}^{\prime }}}\left[\frac{{\delta }_{ij}}{3}+\frac{3({{\tilde{\eta }}^{\prime }}_{ij}-{{\tilde{\eta }}^{\prime }}_{ij0})-{{\tilde{\eta }}^{\prime }}_{mn}({{\tilde{\eta }}^{\prime }}_{mn}-{{\tilde{\eta }}^{\prime }}_{mn0}){\delta }_{ij}}{{\tilde{M}}^{\ast }{}^{2}+{\tilde{\eta }}^{\ast }{}^{2}}\right]\text{，}\\ \frac{\partial {{\tilde{\sigma }}^{\prime }}_{ij}}{\partial {{\sigma }^{\prime }}_{kl}}=\frac{1}{3}{\delta }_{kl}{\delta }_{ij}+\frac{{s}_{ij}}{q}\frac{\partial {q}_{\text{c}}}{\partial {{\sigma }^{\prime }}_{kl}}+\frac{{q}_{\text{c}}}{q}\left({\delta }_{ik}{\delta }_{jl}-\frac{1}{3}{\delta }_{kl}{\delta }_{ij}-\frac{{s}_{ij}{s}_{kl}}{{q}^{2}}\right)\text{，}\\ \frac{\partial {q}_{\text{c}}}{\partial {{\sigma }^{\prime }}_{kl}}=\left(\frac{2}{3F-1}+R{I}_{2}{I}_{3}\right){\delta }_{ij}+R{I}_{1}{I}_{3}({\sigma }_{mm}{\delta }_{ij}-{\sigma }_{ij})\\ \text{ }-R{I}_{1}{I}_{2}\left({\sigma }_{im}{\sigma }_{mj}-{\sigma }_{ij}{\sigma }_{mm}-\frac{1}{2}{\sigma }_{sr}{\sigma }_{sr}{\delta }_{ij}+\frac{1}{2}{\sigma }_{mm}^{2}{\delta }_{ij}\right)\text{，}\\ R=\frac{24{I}_{1}}{F{(3F-1)}^{2}{({I}_{1}{I}_{2}-9{I}_{3})}^{2}}\text{，}\\ F=\sqrt{({I}_{1}{I}_{2}-{I}_{3})/({I}_{1}{I}_{2}-9{I}_{3})}。\end{array} \right\}$

(17b)

（2）采用 $g(\theta )$ 方法三维化的K₀-UH模型的弹塑性张量 ${D}_{ijkl}^{{\text{ep}}}$ 计算式为

$\mathit{\boldsymbol{D}}_{ijkl}^{{\text{ep}}} = \mathit{\boldsymbol{D}}_{ijkl}^{\text{e}} - \frac{{\mathit{\boldsymbol{D}}_{ijmn}^e\frac{{\partial g}}{{\partial {{\sigma '}_{mn}}}}\frac{{\partial f}}{{\partial {{\sigma '}_{st}}}}\mathit{\boldsymbol{D}}_{stkl}^{\text{e}}}}{{ - \frac{{\partial f}}{{\partial H}}\frac{{\partial H}}{{\partial \varepsilon _{\text{v}}^{\text{p}}}}\frac{{\partial g}}{{\partial p'}} + \frac{{\partial f}}{{\partial {{\sigma '}_{uv}}}}\mathit{\boldsymbol{D}}_{uvwz}^{\text{e}}\frac{{\partial g}}{{\partial {{\sigma '}_{wz}}}}}} 。$

(18)

上式同样可以以矩阵形式（17）表示，式中各变量为

$\left. \begin{array}{l}{\alpha }_{i}=L\left(\frac{\partial g}{\partial {{\sigma }^{\prime }}_{1}}+\frac{\partial g}{\partial {{\sigma }^{\prime }}_{2}}+\frac{\partial g}{\partial {{\sigma }^{\prime }}_{3}}\right)+2G\frac{\partial g}{\partial {{\sigma }^{\prime }}_{i}}\text{ }\text{，}\\ {\beta }_{i}=L\left(\frac{\partial f}{\partial {{\sigma }^{\prime }}_{1}}+\frac{\partial f}{\partial {{\sigma }^{\prime }}_{2}}+\frac{\partial f}{\partial {{\sigma }^{\prime }}_{3}}\right)+2G\frac{\partial f}{\partial {{\sigma }^{\prime }}_{i}}\text{ }\text{，}\\ N=\frac{{M}_{f}^{4}-{\eta }_{\text{k}}^{4}}{{c}_{p}{p}^{\prime }({M}^{4}-{\eta }_{\text{k}}^{4})}\left[1-\frac{3({{\eta }^{\prime }}_{mn}{{\eta }^{\prime }}_{mn}-{{\eta }^{\prime }}_{mn}{{\eta }^{\prime }}_{mn0})}{{({M}^{\ast }g({\theta }^{\ast }))}^{2}+{\tilde{\eta }}^{\ast }{}^{2}}\right]+{\beta }_{i}\frac{\partial g}{\partial {{\sigma }^{\prime }}_{i}}。\end{array} \right\}$

(18a)

式中：

$\left. \begin{array}{l}\frac{\partial f}{\partial {{\sigma }^{\prime }}_{ij}}=\frac{1}{{p}^{\prime }}\left[\frac{{\delta }_{ij}}{3}+\frac{3({{\eta }^{\prime }}_{ij}-{{\eta }^{\prime }}_{ij0})-{{\eta }^{\prime }}_{mn}({{\eta }^{\prime }}_{mn}-{{\eta }^{\prime }}_{mn0}){\delta }_{ij}}{{[{M}^{\ast }g({\theta }^{\ast })]}^{2}+{\eta }^{\ast }{}^{2}}\right]+\\ \text{ }\frac{\partial f}{\partial g({\theta }^{\ast })}\frac{\partial g({\theta }^{\ast })}{\partial {\theta }^{\ast }}\left(\frac{\partial {\theta }^{\ast }}{\partial {J}_{2}^{\ast }}\frac{\partial {J}_{2}^{\ast }}{\partial {\sigma }_{ij}}+\frac{\partial {\theta }^{\ast }}{\partial {J}_{3}^{\ast }}\frac{\partial {J}_{3}^{\ast }}{\partial {\sigma }_{ij}}\right)\text{，}\\ \frac{\partial g}{\partial {{\sigma }^{\prime }}_{ij}}=\frac{1}{{p}^{\prime }}\left[\frac{{\delta }_{ij}}{3}+\frac{3({{\eta }^{\prime }}_{ij}-{{\eta }^{\prime }}_{ij0})-{{\eta }^{\prime }}_{mn}({{\eta }^{\prime }}_{mn}-{{\eta }^{\prime }}_{mn0}){\delta }_{ij}}{{[{M}^{\ast }g({\theta }^{\ast })]}^{2}+{\eta }^{\ast }{}^{2}}\right]\text{，}\\ \frac{\partial f}{\partial g({\theta }^{\ast })}=\frac{-2{\eta }^{\ast 2}}{{M}^{\ast 2}{g}^{3}({\theta }^{\ast })+{\eta }^{\ast }{}^{2}g({\theta }^{\ast })}\text{，}\\ \frac{\partial g({\theta }^{\ast })}{\partial {\theta }^{\ast }}=\frac{-\mathrm{tan}\psi }{\mathrm{sin}3\psi }{\left(\frac{3}{2+{\mathrm{sin}}^{2}{\phi }_{0}}\right)}^{3/2}g({\theta }^{\ast })\mathrm{sin}{\phi }_{0}\mathrm{cos}3{\theta }^{\ast }\text{，}\\ \psi =\frac{1}{3}\mathrm{arccos}\left[{\left(\frac{3}{2+{\mathrm{sin}}^{2}{\phi }_{0}}\right)}^{3/2}\mathrm{sin}{\phi }_{0}\mathrm{sin}3{\theta }^{\ast }\right]\text{，}\\ \frac{\partial {\theta }^{\ast }}{\partial {J}_{2}^{\ast }}=\frac{3\sqrt{3}{J}_{3}^{\ast }}{4{J}_{3}^{\ast 5/2}}\text{ }, \text{ }\frac{\partial {\theta }^{\ast }}{\partial {J}_{3}^{\ast }}=\frac{-\sqrt{3}}{2{J}_{2}^{\ast 3/2}}\text{ }。\end{array} \right\}$

(18b)

2.2 三维化方法对应力应变计算影响分析

由上文分析可知，采用TS和 $g(\theta )$ 方法三维化改变了本构模型的屈服面在π平面上的形状，从而改变了偏应力q在不同罗德角下对变形的影响。假设土体从各向同性p₀开始沿等p路径加载，则加载至任意应力比η时塑性体应变的增量为

$\Delta \varepsilon _{\text{v}}^{\text{p}} = {c_{\text{p}}}\ln \frac{{{p_x}}}{{{p_0}}} \text{，}$

(19)

式中：p_x为加载至应力比η时的当前硬化参数，参数p_x1（TS）和p_x2（ $g(\theta )$ ）的计算式分别为

$\left. \begin{array}{l}{p}_{x1}={p}_{0}\left(1+\frac{{\eta }^{\ast }{}^{2}}{{M}^{2}}\right)\text{ }\text{，}\\ {p}_{x2}={p}_{0}\left(1+\frac{{\eta }^{2}}{{M}^{2}{g}^{2}\left(\theta \right)}\right)\text{ }。\end{array} \right\}$

(20)

采用两种方法三维化后的UH模型加载至相同的应力比η，不同罗德角下屈服面的变化如图 7所示。

图 7 不同罗德角下三维本构模型屈服面在p-q面形状变化的对比分析

Figure 7. Comparison of variation in shape of yield surface for three-dimensional constitutive model under different Lode angles in p-q plane

下载: 全尺寸图片幻灯片

中，虚线为未三维化的UH模型屈服面，实线为采用TS或 $g(\theta )$ 方法三维化后的屈服面。从图 7中可以看出，当罗德角θ等于0°时，采用两种方法三维化后的屈服面在p-q平面上均与三维化前的屈服面相同，因而得到的体应变结果也相同。然而，当 $\theta$ 大于0°时，在相同的偏应力增量下，由 $g(\theta )$ 方法的计算得到的体应变值明显大于TS方法得到的计算值。这表明在上述加载条件下，采用TS方法计算的结果偏“硬”，采用 $g(\theta )$ 方法计算偏“软”。

3. 平面应变单元计算

3.1 加载路径

以文献[18]中已有的平面应变单元试验为例，对比分析两种三维化方法对土的应力应变计算结果的影响。图 8给出了平面应变单元试验采用4条应力路径AB、AB′、AC和AC′。其中，点A为土体的初始应力状态（K₀=0.5），路径AB的应力条件为 ${\sigma _x}$ =98 kPa，路径AB′为 ${\sigma _y}$ =196 kPa，路径AC和AC′为( ${\sigma _x}$ + ${\sigma _y}$ )/2 =147 kPa。

图 8 应力加载路径

Figure 8. Stress path under loading

下载: 全尺寸图片幻灯片

3.2 单元试验结果

文献[]中黏土采用的参数为： $\lambda$ /(1+e₀)=0.0508，κ/(1+e₀)=0.0112，ϕ=33.7°， $\nu$ =0.3，试验和预测结果如图 9所示。从中可以看出，采用TS方法三维化后的K₀-UH模型能够更好地预测平面应变条件下的应力应变特性，而采用 $g(\theta )$ 方法三维化后的模型计算的应力比偏低且在相同的应力比下变形较大。这一规律与本文2.2的分析结果相同。

图 9 平面应变单元试验的模型预测值和实测值

Figure 9. Predicted and measured results for plane strain element tests

下载: 全尺寸图片幻灯片

为了进一步分析两种三维化方法在预测平面应变单元试验时的差异，图 10展示了沿AB′应力路径加载时计算得到的π平面上的应力变化轨迹，从可以看出，相比采用TS方法得到的SMP破坏面，采用 $g(\theta )$ 方法计算的破坏面呈现不规则形状，且在该应力路径下低估了临界破坏时的应力比。

图 10 沿AB′应力路径加载时π平面应力路径

Figure 10. Predicted stress path on π plane when loading along AB'

下载: 全尺寸图片幻灯片

4. 平面应变边值计算——柱孔扩张

4.1 柱孔扩张计算模型

对土体中柱孔扩张问题的计算分析，已有文献[21, 22]常简化为一个圆柱形孔的平面应变扩张模型，其中竖向应变为零（dε_z=0），如图 11所示。土体初始孔径为a₀，扩孔后至a半径处时扩张压力为 ${\sigma _a}$ ；孔周土体单元的径向应力表示为 ${\sigma _r}$ ，切向应力表示为 ${\sigma _\theta }$ ；无穷远处土体的应力状态为原位应力状态，该处水平方向（径向）和切向应力大小均为 ${\sigma _{{\text{h0}}}}$ ，竖直方向应力大小为 ${\sigma _{{\text{v0}}}}$ 。在UH模型的假定中，加载时应力始终处于屈服面上，因此孔周只存在塑性区。

图 11 柱孔扩张模型

Figure 11. Sketch for modeling expansion of a cylindrical cavity

下载: 全尺寸图片幻灯片

4.2 不排水控制方程

式（17）的应力应变关系进一步表示为

$\left[ {\begin{array}{*{20}{c}} {{\text{d}}{{\sigma '}_r}} \\ {{\text{d}}{{\sigma '}_\theta }} \\ {{\text{d}}{{\sigma '}_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{D}_{rr}}}&{{{D}_{r\theta }}}&{{{D}_{rz}}} \\ {{{D}_{\theta r}}}&{{{D}_{\theta \theta }}}&{{{D}_{\theta z}}} \\ {{{D}_{zr}}}&{{{D}_{z\theta }}}&{{{D}_{zz}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\text{d}}{\varepsilon _r}} \\ {{\text{d}}{\varepsilon _\theta }} \\ {{\text{d}}{\varepsilon _z}} \end{array}} \right] 。$

(21)

式中：D_ij分别为式（17）中的各项。不排水条件下的孔周径向和切向应变为^[23]

${\varepsilon _r} = - {\varepsilon _\theta } = - \frac{1}{2}\ln \left( {1 - \frac{m}{{{\zeta ^2}}}} \right) 。$

(22)

式中： $\zeta$ =r/a，m=1-(a₀/a)²。将式（22）代入式（21），得到不排水土体应力分量随半径变化的微分方程组：

$\left. \begin{array}{l}\frac{\text{d}{{\sigma }^{\prime }}_{r}}{\text{d}\zeta }=({D}_{rr}-{D}_{r\theta })\frac{m}{{\zeta }^{3}-m\zeta }\text{ }\text{，}\\ \frac{\text{d}{{\sigma }^{\prime }}_{\theta }}{\text{d}\zeta }=({D}_{\theta r}-{D}_{\theta \theta })\frac{m}{{\zeta }^{3}-m\zeta }\text{ }\text{，}\\ \frac{\text{d}{{\sigma }^{\prime }}_{z}}{\text{d}\zeta }=({D}_{zr}-{D}_{z\theta })\frac{m}{{\zeta }^{3}-m\zeta }\text{ }。\end{array} \right\}$

(23)

式（23）即为不排水柱孔扩张控制方程。

4.3 排水控制方程

排水条件下孔周土体的体积在不断变化，因此新增加了一个比体积变量 $\upsilon$ ，而方程（21）只有3个，显然方程的数量不够。考虑在式（21）的基础上加入孔周土单元的平衡方程：

$\frac{{{\text{d}}{{\sigma '}_r}}}{{{\text{d}}r}} + \frac{{{{\sigma '}_r} - {{\sigma '}_\theta }}}{r} = 0 \text{，}$

(24)

引入辅助变量ϑ，令ϑ=r₀/r，孔周径向和切向应变计算式为^[25]

${\varepsilon _\theta } = \ln (\vartheta ) \text{，}$

(25)

${\varepsilon _r} = {\varepsilon _{\text{v}}} - {\varepsilon _\theta } = \ln \left( {\frac{{{\upsilon _0}}}{\upsilon }} \right) - \ln (\vartheta ) = \ln \left( {\frac{{{\upsilon _0}}}{{\upsilon \vartheta }}} \right) \text{，}$

(26)

式中： ${\upsilon _0}$ 为土体在初始应力下的比体积。将辅助变量ϑ微分，有如下关系：

${\text{d}}\vartheta = \left( {\frac{1}{r}\frac{{{\upsilon _0}}}{{\upsilon \vartheta }} - \frac{1}{r}\vartheta } \right){\text{d}}r 。$

(27)

结合式（21），（24）~（27），得到控制方程组：

$\left. \begin{array}{l}\frac{\text{d}{{\sigma }^{\prime }}_{r}}{\text{d}\vartheta }=-\frac{{{\sigma }^{\prime }}_{r}-{{\sigma }^{\prime }}_{\theta }}{\left(\frac{{\upsilon }_{0}}{\upsilon \vartheta }-\vartheta \right)}\text{ }\text{，}\\ \frac{\text{d}{{\sigma }^{\prime }}_{\theta }}{\text{d}\vartheta }=-\frac{{D}_{\theta r}}{{D}_{rr}}\left(\frac{{D}_{rr}-{D}_{r\theta }}{\vartheta }+\frac{{{\sigma }^{\prime }}_{r}-{{\sigma }^{\prime }}_{\theta }}{\frac{{\upsilon }_{0}}{\upsilon \vartheta }-\vartheta }\right)+\frac{{D}_{\theta \theta }-{D}_{\theta r}}{\vartheta }\text{，}\\ \frac{\text{d}{{\sigma }^{\prime }}_{z}}{\text{d}\vartheta }=-\frac{{D}_{zr}}{{D}_{rr}}\left(\frac{{D}_{r\theta }-{D}_{rr}}{\vartheta }+\frac{{{\sigma }^{\prime }}_{r}-{{\sigma }^{\prime }}_{\theta }}{\frac{{\upsilon }_{0}}{\upsilon \vartheta }-\vartheta }\right)+\frac{{D}_{z\theta }-{D}_{zr}}{\vartheta }\text{，}\\ \frac{\text{d}\upsilon }{\text{d}\vartheta }=\frac{\upsilon }{{D}_{rr}}\left(\frac{{D}_{r\theta }-{D}_{rr}}{\omega \vartheta }+\frac{{{\sigma }^{\prime }}_{r}-{{\sigma }^{\prime }}_{\theta }}{\frac{{\upsilon }_{0}}{\upsilon \vartheta }-\vartheta }\right)\text{ }。\end{array} \right\}$

(28)

式（28）前两式分别为排水柱孔扩张应力分量和比体积的控制方程。由于控制方程计算结果为应力分量与辅助变量ϑ的关系，需通过对式（27）积分，转化为应力分量与半径r的关系：

$\frac{r}{a} = \exp \left( {\int_{\vartheta (a)}^\vartheta {\frac{1}{{\frac{{{\upsilon _0}}}{{\upsilon \vartheta }} - \vartheta }}} {\text{d}}\vartheta } \right) \text{，}$

(29)

式中：ϑ(a)=a₀/a为孔壁处辅助变量ϑ的大小。

4.4 边界条件

不排水条件下，方程组（23）的应力和位移初始条件为^[23]

${\sigma '_r}({r_0}) = {\sigma '_{r0}}, {\text{ }}{\sigma '_\theta }({r_0}) = {\sigma '_{\theta 0}}, {\text{ }}{\sigma '_z}({r_0}) = {\sigma '_{z0}} \text{，}$

(30)

$\frac{{{r_0}}}{a} = \sqrt {{{\left( {\frac{r}{a}} \right)}^2} + {{\left( {\frac{{{a_0}}}{a}} \right)}^2} - 1} 。$

(31)

排水条件下，方程组（28）的应力和比体积边界条件为^[25]

${\sigma _r}(1) = {\sigma '_{r0}}, {\text{ }}{\sigma _\theta }(1) = {\sigma '_{\theta 0}}, {\text{ }}{\sigma _z}(1) = {\sigma '_{z0}}, {\text{ }}\upsilon (1) = 1 + {e_0} \text{，}$

(32)

位移边界条件为

${\vartheta _0} = 1 。$

(33)

4.5 不排水柱孔扩张计算结果

本节选取Boston Blue clay^[21]进行计算，土的具体参数见表 1。表中初始孔隙比计算式为

${\upsilon _0} = {\upsilon _{{\text{cs}}}} + (\lambda - \kappa )\ln 2 - \lambda \ln {p'_0} + \ln {\left( {1 + \frac{{q_0^2}}{{{M^2}p_0^2}}} \right)^{\kappa - \lambda }} 。$

(34)

表 1 Boston Blue clay参数

Table 1. Parameters for modeling Boston Blue clay

K₀	${\sigma '_{r{\text{0}}}}$ / kPa	${\sigma '_{\theta {\text{0}}}}$ / kPa	${\sigma '_{z{\text{0}}}}$ / kPa	${p'_{\text{0}}}$ / kPa	${q'_{\text{0}}}$ / kPa	${\upsilon _0}$
1	100	100	100	100	0	2.01
0.625	100	100	160	120	60	2.09
0.4	100	100	250	150	150	2.13
λ=0.15; κ=0.03; υ_cs=2.74

下载: 导出CSV

| 显示表格

不排水条件下柱孔扩张计算结果如图 12所示。从可以看出，当参数K₀等于1时，采用 $g(\theta )$ 方法或TS方法计算的应力分量在临界破坏时收敛于相同的值，且切向有效应力收敛于径向和竖向有效应力的平均值，即 ${\sigma '_\theta }$ =( ${\sigma '_r}$ + ${\sigma '_z}$ )/2。当参数K₀不等于1时，切向有效应力与径向和竖向有效应力无必然关系，该结论与文献[21, 22]等结论相同。需要特别指出的是，不排水条件下两种方法计算结果的差异不大。

图 12 不排水柱孔扩张的有效应力分量变化规律

Figure 12. Predicted variation in effective stresses for cylindrical cavity expansion under undrained conditions

下载: 全尺寸图片幻灯片

4.6 排水柱孔扩张计算结果

排水扩孔的应力分量和比体积的变化结果如图 13所示。当K₀=1时，不同方法计算的应力分量在临界破坏时收敛于相同的值。但从，可以看出，TS方法和 $g(\theta )$ 方法计算结果差异随初始各向异性的增大而增大，尤其在K₀=0.4的土体中，采用不同方法计算的孔壁竖向应力有很大差异。从比体积的变化可以看出， $g(\theta )$ 方法计算得到的比体积变化较大，TS方法计算较小，因此采用 $g(\theta )$ 方法计算偏“软”，采用TS方法计算的结果偏“硬”。

图 13 排水柱孔扩张有效应力分量和比体积变化规律

Figure 13. Predicted variation in effective stresses and specific volume for cylindrical cavity expansion under drained conditions

下载: 全尺寸图片幻灯片

利用柱孔扩张的自相似的特性^[24]，将K₀=0.4的孔壁土体在加载过程中π平面上应力路径表示在图 14中。

图 14 排水柱孔扩张孔壁土体在π平面上的应力路径

Figure 14. Predicted stress path on π plane for cylindrical cavity expansion under drained conditions

下载: 全尺寸图片幻灯片

从中可以看出，由于计算的破坏面呈现不规则形状，采用 $g(\theta )$ 方法过高地估计了孔壁土体临界破坏的应力比。因此，排水条件下在各向异性较大的土体中预测柱孔扩张的应力场时，建议采用TS方法三维化后的模型进行计算。

5. 结论

本文从机理上详细分析了采用TS方法和 $g(\theta )$ 方法实现K₀-UH模型三维化在描述强度和变形特性等方面的差异，并通过对平面应变单元试验和边值问题的分析，阐述了这一差异性对平面应变条件下土的应力应变计算的影响，得到以下3点结论。

（1）在强度方面，初始应力各向同性且塑性势面相同的情况下，TS方法和 $g(\theta )$ 方法在临界破坏时收敛于相同的值。当土体具有初始应力各向异性时，采用 $g(\theta )$ 方法三维化后各向异性屈服面由低应力比下的曲边三角形逐渐变为内凹形的不规则破坏面，预测平面应变试验的临界状态强度往往会偏高或偏低。TS方法是将屈服曲线在π平面上由低应力比下的近圆形转变为剪切破坏时的SMP破坏准则形状，临界破坏强度与实测值吻合。

（2）在变形方面，将土体加载至相同的应力状态时， $g(\theta )$ 方法计算出的应变结果偏大，导致计算结果偏“软”；TS方法计算的应变相对较小，计算结果偏“硬”。与文献中平面应变单元试验结果相比，TS方法更为合理。

（3）在柱孔扩张边值问题中，不排水条件下两种三维化方法计算的结果差异不大；而排水条件下，随着各向异性增大，由于计算的破坏面呈现不规则形状，采用 $g(\theta )$ 方法计算值过高地估计了孔壁土体临界破坏的应力比。因此，排水条件下在预测初始各向异性较大的土体柱孔扩张的应力场时，建议采用TS方法三维化后的模型进行计算。

图 1 采用TS方法三维化的K₀-UH屈服面

Figure 1. Three-dimensional yield surface of K₀-UH model using TS method

下载: 全尺寸图片幻灯片

图 2 采用g(θ)方法三维化的K₀-UH屈服面

Figure 2. Three-dimensional yield surface of K₀-UH model using g(θ) method

下载: 全尺寸图片幻灯片

图 3 采用TS方法的初始各向同性三维屈服面和破坏面（K₀=1）

Figure 3. Three-dimensional yield surface and failure surface using TS method with initial isotropic stress status (K₀=1)

下载: 全尺寸图片幻灯片

图 4 采用g(θ)方法的初始各向同性三维屈服面和破坏面（K₀=1）

Figure 4. Three-dimensional yield surface and failure surface using g(θ) method with initial isotropic stress status (K₀=1)

下载: 全尺寸图片幻灯片

图 5 采用TS方法的初始各向异性三维屈服面和破坏面（K₀=0.625）

Figure 5. Three-dimensional yield surface and failure surface using TS method with initial anisotropic stress status (K₀=0.625)

下载: 全尺寸图片幻灯片

图 6 采用g(θ)方法的初始各向异性三维屈服面和破坏面（K₀=0.625）

Figure 6. Three-dimensional yield surface and failure surface using g(θ) method with initial anisotropic stress status (K₀=0.625)

下载: 全尺寸图片幻灯片

图 7 不同罗德角下三维本构模型屈服面在p-q面形状变化的对比分析

Figure 7. Comparison of variation in shape of yield surface for three-dimensional constitutive model under different Lode angles in p-q plane

下载: 全尺寸图片幻灯片

图 8 应力加载路径

Figure 8. Stress path under loading

下载: 全尺寸图片幻灯片

图 9 平面应变单元试验的模型预测值和实测值

Figure 9. Predicted and measured results for plane strain element tests

下载: 全尺寸图片幻灯片

图 10 沿AB′应力路径加载时π平面应力路径

Figure 10. Predicted stress path on π plane when loading along AB'

下载: 全尺寸图片幻灯片

图 11 柱孔扩张模型

Figure 11. Sketch for modeling expansion of a cylindrical cavity

下载: 全尺寸图片幻灯片

图 12 不排水柱孔扩张的有效应力分量变化规律

Figure 12. Predicted variation in effective stresses for cylindrical cavity expansion under undrained conditions

下载: 全尺寸图片幻灯片

图 13 排水柱孔扩张有效应力分量和比体积变化规律

Figure 13. Predicted variation in effective stresses and specific volume for cylindrical cavity expansion under drained conditions

下载: 全尺寸图片幻灯片

图 14 排水柱孔扩张孔壁土体在π平面上的应力路径

Figure 14. Predicted stress path on π plane for cylindrical cavity expansion under drained conditions

下载: 全尺寸图片幻灯片

表 1 Boston Blue clay参数

Table 1 Parameters for modeling Boston Blue clay

K₀	${\sigma '_{r{\text{0}}}}$ / kPa	${\sigma '_{\theta {\text{0}}}}$ / kPa	${\sigma '_{z{\text{0}}}}$ / kPa	${p'_{\text{0}}}$ / kPa	${q'_{\text{0}}}$ / kPa	${\upsilon _0}$
1	100	100	100	100	0	2.01
0.625	100	100	160	120	60	2.09
0.4	100	100	250	150	150	2.13
λ=0.15; κ=0.03; υ_cs=2.74

下载: 导出CSV

参考文献(24)

[1]	WROTH C, HOULSBY G. Soil mechanics-property characterization and analysis procedures[C]// Proc of 11th Int Conf on SMFE, San Francisce, 1985.
[2]	ZIENKIEWICZ O. Some useful forms of isotropic yield surfaces for soil and rock mechanics[C]// Finite Elements in Geomechanics. London: John Wileye Sons, 1977.
[3]	POTTS D M, GENS A. The effect of the plastic potential in boundary value problems involving plane strain deformation[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(3): 259-286. doi: 10.1002/nag.1610080305
[4]	MATSUOKA H, YAO Y P, SUN D A. The Cam-clay models revised by the SMP criterion[J]. Soils and Foundations, 1999, 39(1): 81-95. doi: 10.3208/sandf.39.81
[5]	ROUAINIA M, MUIR WOOD D. A kinematic hardening constitutive model for natural clays with loss of structure[J]. Géotechnique, 2000, 50(2): 153-164. doi: 10.1680/geot.2000.50.2.153
[6]	LI X S. A sand model with state-dapendent dilatancy[J]. Géotechnique, 2002, 52(3): 173-186. doi: 10.1680/geot.2002.52.3.173
[7]	DAFALIAS Y F, MANZARI M T, PAPADIMITRIOU A G. SANICLAY: simple anisotropic clay plasticity model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(12): 1231-1257. doi: 10.1002/nag.524
[8]	YAO Y P, LIU L, LUO T, et al. Unified hardening (UH) model for clays and sands[J]. Computers and Geotechnics, 2019, 110: 326-343. doi: 10.1016/j.compgeo.2019.02.024
[9]	CHOWDHURY E Q, NAKAI T R, TAWADA M, et al. A model for clay using modified stress under various loading conditions with the application of subloading concept[J]. Soils and Foundations, 1999, 39(6): 103-116. doi: 10.3208/sandf.39.6_103
[10]	GRAMMATIKOPOULOU A, ZDRAVKOVIC L, POTTS D M. General formulation of two kinematic hardening constitutive models with a smooth elastoplastic transition[J]. International Journal of Geomechanics, 2006, 6(5): 291-302. doi: 10.1061/(ASCE)1532-3641(2006)6:5(291)
[11]	GRAMMATIKOPOULOU A, ZDRAVKOVIC L, POTTS D M. The effect of the yield and plastic potential deviatoric surfaces on the failure height of an embankment[J]. Géotechnique, 2007, 57(10): 795-806. doi: 10.1680/geot.2007.57.10.795
[12]	HASHIGUCHI K. A proposal of the simplest convex-conical surface for soils[J]. Soils and Foundations, 2002, 42(3): 107-113. doi: 10.3208/sandf.42.3_107
[13]	NAKAI T R, HINOKIO M. A simple elastoplastic model for normally and over consolidated soils with unified material parameters[J]. Soils and Foundations, 2004, 44(2): 53-70. doi: 10.3208/sandf.44.2_53
[14]	YAO Y P, SUN D A. Application of lade's criterion to cam-clay model[J]. Journal of Engineering Mechanics, 2000, 126(1): 112-119. doi: 10.1061/(ASCE)0733-9399(2000)126:1(112)
[15]	YAO Y P, LU D C, ZHOU A N, et al. Generalized non-linear strength theory and transformed stress space[J]. Science in China Series E: Technological Sciences, 2004, 47(6): 691-709. doi: 10.1360/04ye0199
[16]	YAO Y P, HOU W, ZHOU A N. UH model: three-dimensional unified hardening model for overconsolidated clays[J]. Géotechnique, 2009, 59(5): 451-469. doi: 10.1680/geot.2007.00029
[17]	YAO Y, TIAN Y, GAO Z. Anisotropic UH model for soils based on a simple transformed stress method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2017, 41(1): 54-78. doi: 10.1002/nag.2545
[18]	孙徳安, 松冈元, 姚仰平, 等. 初期异方性を考虑した黏土と砂の统一的な弾塑性构成式[C]// 土木学会论文集, 1999. SUN Dean, MATSUOKA H, YAO Yangping, et al. An anisotropic unified hardening elastoplastic model for clays and sands[C]// Preceedings of Japan Society of Civil Engineers, 1999. (in Japanese)
[19]	SUN D A, YAO Y P, MATSUOKA H. Modification of critical state models by Mohr-Coulomb criterion[J]. Mechanics Research Communications, 2006, 33(2): 217-232. doi: 10.1016/j.mechrescom.2005.05.006
[20]	YAO Y P, WANG N D. Transformed stress method for generalizing soil constitutive models[J]. Journal of Engineering Mechanics, 2014, 140(3): 614-629. doi: 10.1061/(ASCE)EM.1943-7889.0000685
[21]	CHEN S L, ABOUSLEIMAN Y N. Exact undrained elasto-plastic solution for cylindrical cavity expansion in modified Cam Clay soil[J]. Géotechnique, 2012, 62(5): 447-456. doi: 10.1680/geot.11.P.027
[22]	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
[23]	武孝天, 徐永福. 基于CSUH模型的砂/黏土不排水柱孔扩张统一解[J]. 岩土工程学报, 2021, 43(6): 1019-1028. doi: 10.11779/CJGE202106005 WU Xiaotian, XU Yongfu. Undrained unified solutions to cylindrical cavity expansion in soils and sands based on CSUH model[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(6): 1019-1028. (in Chinese) doi: 10.11779/CJGE202106005
[24]	武孝天, 徐永福. 超固结土中排水圆孔扩张弹塑性UH解[J]. 岩土工程学报, 2020, 42(10): 1903-1913. doi: 10.11779/CJGE202010016 WU Xiaotian, XU Yongfu. Elasto-plastic solution for drained cavity expansion in over-consolidated soil incorporating three-dimensional unified hardening model[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(10): 1903-1913. (in Chinese) doi: 10.11779/CJGE202010016

施引文献(4)

期刊类型引用(2)

1.	戴玉明，贺佐跃. 粉煤灰对水泥土三轴剪切峰值强度的影响性分析. 广东土木与建筑. 2024(03): 18-22 . 百度学术
2.	戴玉明，贺佐跃. 基于三轴试验的水泥土抗剪强度参数分析. 交通科学与工程. 2024(03): 65-73 . 百度学术