基于微极理论的冻结砂土应变局部化数值模拟

doi:10.7522/j.issn.1000-0240.2023.0085

基于微极理论的冻结砂土应变局部化数值模拟

王冬勇^,¹^,², 邵博¹, 马玲³, 齐吉琳^,¹^,², 彭丽云¹^,²

1.北京建筑大学土木与交通工程学院，北京 100044

2.北京建筑大学城市交通基础设施建设北京市国际科技合作基地，北京 100044

3.北方民族大学土木工程学院，宁夏银川 750030

Numerical simulation of strain localization for frozen sand based on micropolar theory

WANG Dongyong^,¹^,², SHAO Bo¹, MA Ling³, QI Jilin^,¹^,², PENG Liyun¹^,²

1.School of Civil and Transportation Engineering，Beijing University of Civil Engineering and Architecture，Beijing 100044，China

2.Beijing International Science and Technology Cooperation Base for Urban Transportation Infrastructure Construction，Beijing University of Civil Engineering and Architecture，Beijing 100044，China

3.School of Civil Engineering，North Minzu University，Yinchuan 750030，China

通讯作者: 齐吉琳，教授，主要从事寒区岩土工程研究. E-mail： jilinqi@bucea.edu.cn

收稿日期: 2022-11-28 修回日期: 2023-01-26

基金资助:

国家自然科学基金项目. 42002277. 41972279. 42172299
北京建筑大学金字塔人才培养工程项目. JDYC20220807

Received: 2022-11-28 Revised: 2023-01-26

作者简介 About authors

王冬勇，讲师，主要从事计算岩土力学研究.E-mail：wangdongyong@bucea.edu.cn , E-mail：wangdongyong@bucea.edu.cn

摘要

应变局部化问题是冻土材料破坏的前兆，目前相关研究还较为薄弱。在微极连续体理论框架下，通过考虑温度对相关模型参数的影响，建立了可以模拟冻土应变局部化问题的微极弹塑性数值框架。通过对冻结砂土平面应变试验进行数值分析，发现所建立的数值框架可实现对冻土剪切带产生、发展和形成的全过程模拟，且模拟获得的应力-应变曲线及剪切带破坏形态与试验结果吻合较好；进而探究了温度、端部约束及初始缺陷对剪切带发展规律的影响，发现温度对剪切带倾角和宽度影响较为显著，倾角随温度的降低而逐渐增大，宽度随温度的降低而逐渐加宽，具体剪切带的破坏形态和破坏位置需要综合考虑端部约束及初始缺陷的共同作用；通过对不同有限元网格条件下的冻土应变局部化进行分析，发现粗网格和细网格获得的模拟结果基本一致，说明所建立的数值框架可以基本克服网格依赖性问题。

关键词： 冻结砂土 ; 微极理论 ; 应变局部化 ; 平面应变 ; 数值模拟

Abstract

Strain localization is the precursor of the failure of frozen soil. So far， the research on strain localization of frozen soil is relatively weak. A micropolar elastoplastic numerical framework is established to simulate the strain localization of frozen soil by considering the effect of temperature. Through numerical simulation of plane strain test of frozen sand， it is found that the process of the generation， development and formation of shear bands of frozen soil can be simulated by the established numerical framework， and the stress-strain curve and shear band are in good agreement with the test results. Besides， the effects of temperature， end restraints and initial defects on the shear band are investigated. The angle and the width of the shear band increase with the decrease of temperature， and the specific failure mode and location of the shear band are closely related to the initial defects and end restraints. Through the simulation analysis for different finite element mesh， it can be observed that the simulation results obtained based on the coarse mesh are consistent with the fine mesh， which indicates that the mesh-dependent problem can be basically avoided.

Keywords： frozen sand ; micropolar theory ; strain localization ; plane strain ; numerical simulation

PDF (9851KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

王冬勇, 邵博, 马玲, 齐吉琳, 彭丽云. 基于微极理论的冻结砂土应变局部化数值模拟[J]. 冰川冻土, 2023, 45(3): 1105-1115 doi:10.7522/j.issn.1000-0240.2023.0085

WANG Dongyong, SHAO Bo, MA Ling, QI Jilin, PENG Liyun. Numerical simulation of strain localization for frozen sand based on micropolar theory[J]. Journal of Glaciology and Geocryology, 2023, 45(3): 1105-1115 doi:10.7522/j.issn.1000-0240.2023.0085

0 引言

岩土材料在剪切破坏时，塑性应变通常会集中于某局部狭小带状区域内，即表现出应变局部化现象^［1-2］。在实际工程中，如边坡滑塌、路基失稳破坏、隧道围岩失稳等灾害通常是渐进破坏过程，普遍存在着应变局部化问题。岩土材料的应变局部化问题一直是岩土力学领域的研究热点，对局部应变发展和剪切带形成演化规律的研究是分析实际工程土体稳定性及渐进破坏过程的基础和关键。目前，国内外学者已经在室内试验、理论建模和数值模拟等方面对岩土材料应变局部化问题进行了相关研究^［3］。随着计算机技术和数值计算方法的发展，数值模拟已经逐渐成为研究土体应变局部化问题的一种有效且方便的手段，常采用的数值方法主要包括离散元法和有限元法^［3-6］。

目前，关于融土应变局部化的数值模拟已经相对较为成熟，许多科研人员通过数值计算对融土剪切带的发展规律进行了系统研究^［3-6］。由于冻土具有显著的温度相关性，目前关于冻土应变局部化的研究还相对较为薄弱。考虑到离散元法可以较好地从细观角度揭示冻土的力学行为，相关学者基于颗粒流理论对冻土应变局部化问题进行了离散元模拟。例如，周凤玺等^［7］利用二维颗粒流程序PFC^2D对冻结砂土的平面应变试验进行了模拟，研究表明离散元法可以较好地模拟冻结砂土的应力应变关系以及剪切带的发展过程，但颗粒流细观参数对温度具有显著的依赖性。袁伟等^［8］利用PFC^2D对冻结砂土的三轴试验进行了离散元模拟，分析了不同应变速率条件下剪切带的发展演化规律。为了更真实地反映三维情况，尹楠等^［9］基于三维颗粒流程序PFC^3D对冻结黏土的室内三轴试验进行了数值模拟，研究表明围压增大会使得接触黏结逐渐失效，剪切带中胶结冰的破坏区域将增大，细观参数对温度的依赖性也很明显。可见，利用离散元法研究冻土应变局部化问题具有一定的理论和应用价值，但在计算过程中细观参数的取值往往导致计算结果存在一定的偏差，不断改变参数反复调试则计算量较大。更为关键的问题是，数值计算中所涉及的细观参数目前很难通过试验测得，经常靠经验和假定。

与离散元法相比，有限元法中的宏观参数可以通过试验测定且无需反复调试，同样成为模拟土体应变局部化问题的一种有效方法。然而，直接采用基于经典连续体理论的有限元法模拟应变局部化时，会出现网格依赖性问题，即随着网格加密剪切带宽度逐渐变窄的非物理现象^［10］。为了解决这个问题，科研人员通过引入非局部理论、梯度塑性理论和微极理论等方法来克服网格依赖性^［3］。其中，微极理论是解决网格依赖性问题的一种较为有效且便于实施的方法之一。微极理论在经典连续体理论的基础上增加了转动自由度，同时在本构方程中通过引入内部特征长度参数作为正则化机制，可以有效规避经典连续体理论中的网格依赖性问题。目前，微极理论已经在融土应变局部化分析中应用较为广泛^［11-13］，例如，Tejchman等^［14］建立了基于微极理论的弹塑性本构模型，通过双轴压缩试验对剪切带模式进行了数值模拟，研究发现应变硬化和应变软化对剪切带模式影响显著。李锡夔等^［15］发展了压力相关弹塑性微极连续体模型，并应用于平面应变压缩试验模拟中，研究表明所发展的微极模型可以有效克服网格依赖性问题。常江芳等^［16］则在微极连续体框架下，通过平面应变试验研究了各向异性对剪切带发展规律的影响。考虑到在传统弹塑性计算框架下需要复杂的平衡迭代和应力积分算法，Wang等^［6］建立了基于二阶锥规划理论的微极弹塑性数值框架。以上研究工作主要针对融土，微极连续体理论还未应用于冻土的应变局部化研究中。因此，为了更好地揭示冻土材料的失稳破坏过程及其演化规律，有必要将微极理论应用于冻土的应变局部化数值模拟中，进而为冻土地区工程建设提供一定的理论依据。

本文在微极连续体理论框架下，通过考虑温度对相关模型参数的影响，建立可以模拟冻土应变局部化问题的微极弹塑性数值框架。进而对冻结砂土平面应变试验进行数值模拟，以期验证所建立数值框架分析冻土应变局部化问题的可行性和有效性，同时探究温度、端部约束和初始缺陷对冻结砂土应变局部化发展规律的影响。

1 考虑温度影响的微极连续体模型

对于微极连续体理论，求解域内每一点除位移自由度u_x 和u_y 外，还引入了独立的转动自由度w_z^［12］，即

u = {[\begin{matrix} u_{x}, & u_{y}, & w_{z} \end{matrix}]}^{T}

（1）

除Cauchy应力和应变外，还存在偶应力m_xz 和m_yz，以及对应的微曲率χ_xz 和χ_yz。其应力和应变矢量分别为

σ = {[\begin{matrix} σ_{x}, & σ_{y}, & σ_{z}, & τ_{x y}, & τ_{y x}, & m_{x z}, & m_{y z} \end{matrix}]}^{T}

（2）

ε = {[\begin{matrix} ε_{x}, & ε_{y}, & ε_{z}, & ε_{x y}, & ε_{y x}, & χ_{x z}, & χ_{y z} \end{matrix}]}^{T}

（3）

其中，

\begin{matrix} ε_{x} = \frac{\partial u_{x}}{\partial x}, ε_{y} = \frac{\partial u_{y}}{\partial y}, ε_{x y} = \frac{\partial u_{y}}{\partial x} - w_{z}, ε_{y x} = \frac{\partial u_{x}}{\partial y} + w_{z}, \\ χ_{x z} = \frac{\partial w_{z}}{\partial x}, χ_{y z} = \frac{\partial w_{z}}{\partial y} \end{matrix}\}

（4）

几何方程可表示为

ε = L u

（5）

其中，

L = {[\begin{matrix} \frac{\partial}{\partial x} & 0 & 0 & 0 & \frac{\partial}{\partial y} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 & \frac{\partial}{\partial x} & 0 & 0 & 0 \\ 0 & 0 & 0 & - 1 & 1 & \frac{\partial}{\partial x} & \frac{\partial}{\partial y} \end{matrix}]}^{T}

（6）

总应变增量为弹性应变增量与塑性应变增量之和，即

Δ ε = Δ ε_{e} + Δ ε_{p}

（7）

线弹性应力-应变关系可表示为

Δ σ = D_{e} Δ ε_{e}

（8）

考虑温度T对弹性参数的影响，则

D_{e} = [\begin{matrix} λ (T) + 2 G (T) & λ (T) & λ (T) & 0 & 0 & 0 & 0 \\ λ (T) & λ (T) + 2 G (T) & λ (T) & 0 & 0 & 0 & 0 \\ λ (T) & λ (T) & λ (T) + 2 G (T) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & G (T) + G_{c} (T) & G (T) - G_{c} (T) & 0 & 0 \\ 0 & 0 & 0 & G (T) - G_{c} (T) & G (T) + G_{c} (T) & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 4 G (T) l_{c}^{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 4 G (T) l_{c}^{2} \end{matrix}]

（9）

式中：G_c为微极连续体剪切模量，通常取G_c=0.5G^［10-13］；l_c为内部特征长度；λ和G为拉梅常数，可由弹性模量E和泊松比v计算获得，即

λ (T) = \frac{ν (T) E (T)}{[1 + ν (T)] [1 - 2 ν (T)]}, G (T) = \frac{E (T)}{2 [1 + ν (T)]}

（10）

采用Drucker-Prager强度准则，考虑温度对强度参数M和K的影响，则

F = q + p M (T) - K (T) \leq 0

（11）

D-P准则中的M和K参数可由黏聚力c和内摩擦角φ换算获得。其中，对于平面应变条件下的关联流动法则，采用内切圆形式^［17］，即

M (T) = \frac{3 t a n φ (T)}{\sqrt[]{3 + 4 t a n^{2} φ (T)}}, K = \frac{3 c (T)}{\sqrt[]{3 + 4 t a n^{2} φ (T)}}

（12）

式（11）中的p和q可分别表示为

p = \frac{1}{3} I_{1} = \frac{1}{3} (σ_{x} + σ_{y} + σ_{z}), q = \sqrt[]{3 J_{2}} = (\frac{1}{2} σ^{T} P_{σ} {σ)}^{\frac{1}{2}}

（13）

其中，

P_{σ} = [\begin{matrix} 2 & - 1 & - 1 & 0 & 0 & 0 & 0 \\ - 1 & 2 & - 1 & 0 & 0 & 0 & 0 \\ - 1 & - 1 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 / 2 & 3 / 2 & 0 & 0 \\ 0 & 0 & 0 & 3 / 2 & 3 / 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3 / l_{c}^{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 / l_{c}^{2} \end{matrix}]

（14）

当采用相关联流动法则时，塑性应变增量可表示为

Δ ε_{p} = Δ λ \frac{\partial F}{\partial σ}

（15）

由上可见，通过考虑温度对弹性参数和强度参数的影响，可建立考虑温度影响的微极连续体理想弹塑性框架，即式（7）~（15）。为了考虑应变软化，假设强度参数随等效塑性应变的发展逐渐弱化直到残余强度^［6，18］，如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 强度参数随等效塑性应变的变化关系

Fig. 1 Relationship between strength parameter and variation of equivalent plastic strain

图1中，c₀和φ₀分别为初始黏聚力和初始内摩擦角；c_r和φ_r分别为残余黏聚力和残余内摩擦角；ε_pr为强度软化阶段进入残余阶段的软化参数；ε_p为等效塑性应变，可由式（16）确定^［6］。

ε_{p}^{} = \sqrt[]{\frac{2}{3} ε_{p}^{T} P_{ε} ε_{p}}

（16）

其中，

P_{ε} = [\begin{matrix} 2 / 3 & - 1 / 3 & - 1 / 3 & 0 & 0 & 0 & 0 \\ - 1 / 3 & 2 / 3 & - 1 / 3 & 0 & 0 & 0 & 0 \\ - 1 / 3 & - 1 / 3 & 2 / 3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 / 2 & 1 / 2 & 0 & 0 \\ 0 & 0 & 0 & 1 / 2 & 1 / 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & l_{c}^{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & l_{c}^{2} \end{matrix}]

（17）

2 冻结砂土平面应变试验模拟

关于冻土应变局部化的试验研究还较为薄弱，目前只有Yao等^［19-20］对平面应变状态下的冻结砂土进行了应变局部化测试工作，利用PIV技术研究了温度为-2 ℃和-5 ℃条件下冻结砂土局部应变和剪切带的演化规律。为验证所建立数值框架模拟冻土应变局部化问题的可行性和有效性，本研究主要对Yao等^［20］开展的应变速率为1.0×10^-4 min^-1条件下冻结砂土平面应变试验进行数值模拟，同时探究温度、端部约束和初始缺陷等对冻土剪切带发展规律的影响。该平面应变试验采用的是冻结标准砂，颗粒平均粒径d₅₀为1.0 mm，土样宽度和高度分别为100 mm和200 mm。其几何及有限元网格如图2所示，求解域内共划分800个单元。左右边界为自由边界，不考虑围压；上下边界通过约束水平方向位移考虑端部约束作用，同时上边界采用增量加载的方式施加位移荷载。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 平面应变试样几何及有限元网格

Fig. 2 Geometry and finite element mesh of the samples for plane strain test

2.1　参数确定

冻土的弹性模量常用割线模量表示，可由峰值强度σ_p的一半除以所对应的应变ε_1/2获得，即E=σ_p/（2ε_1/2）^［21］。根据Yao等^［20］的试验结果，当温度分别为-2 ℃和-5 ℃时，弹性模量分别为153 MPa和315 MPa；由于缺少该平面应变条件下冻结砂土的强度参数，通过参考低围压条件下的三轴试验结果^［22-23］，确定两种温度条件下的黏聚力和内摩擦角，进而根据残余强度与峰值强度的比例关系确定残余黏聚力和残余内摩擦角。选取微极连续体剪切模量G_c=0.5G，内部特征长度l_c=d₅₀=1.0 mm。同时，假设两种温度条件下的泊松比分别是0.35和0.30，设置ε_pr均为1.5，具体的模型参数如表1所示。

表1 模型参数

Table 1 Parameters used in the model

T/℃	E/MPa	ν	c/MPa	φ/（°）	c_r/MPa	φ_r/（°）	ε_pr	G_c/MPa	l_c/mm
-2	153	0.35	1.6	15.6	1.06	11.5	1.5	0.5G	1.0
-5	315	0.30	1.8	28.3	1.24	20.4	1.5	0.5G	1.0

新窗口打开| 下载CSV

2.2　不同温度条件下的剪切带发展规律

为了较好地模拟文献［20］中冻结砂土平面应变试验的破坏形态，需考虑试样端部约束和初始缺陷的影响。初始缺陷是指在均匀试样内部设置弱单元，弱单元的强度和刚度都小于其他正常单元。相关科研人员主要通过弱化缺陷单元的相关模型参数来设置缺陷单元，但缺陷单元的参数具体如何取值目前还没有相对统一的定论^［24-25］。本研究主要通过设置缺陷单元的弹性模量和黏聚力是正常单元的一半来实现对缺陷单元的弱化，模型中共设置左上角和右上角两个单元（单元编号分别是781和800）为缺陷单元，如图2所示。模拟获得的不同温度条件下的应力-应变曲线如图3所示，可以看出两种温度下的冻结砂土轴向应力-轴向应变曲线均呈现应变软化，且模拟结果和试验结果吻合较好。同时发现两种温度下的体应变随着轴向应变的增加均呈现先剪缩后剪胀的趋势［图3（b）］，且随着轴向应变的增加，温度对体积变形的影响越来越明显。由于该平面应变试验没有考虑两侧围压，模拟获得的体应变较大，体积膨胀量随着温度的降低而显著增加；目前，尚缺少平面应变条件下体应变的实测数据，图3（b）只提供了体应变-轴向应变的模拟结果，与三轴条件下的规律基本相似^［26-27］。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 应力-应变曲线

Fig. 3 Stress-strain curves： axial stress-axial strain （a） and volumetric strain-axial strain （b）

为进一步研究剪切过程中试样局部应变发展和剪切带形成规律，图4~7分别给出了温度为-2 ℃和-5 ℃条件下等效塑性应变及体应变随轴向应变的发展规律。由于体应变随轴向应变的增加基本呈现先剪缩后剪胀的趋势，为了更直观地观察局部体应变的发展规律，每个轴向应变条件下的体应变均设置单独的应变区间（见图5和图7）。以T=-5 ℃为例进行说明，当轴向应变ε₁=2%时，试样中没有塑性应变产生，只在端部区域产生了体积收缩（体应变为正值），如图7（a）所示；随着加载的进行，当ε₁=5%时，试样端部和中间区域逐渐出现了塑性变形，此时靠近端部区域依然表现为体积收缩，而中间区域逐渐出现剪胀（体应变为负值）；当ε₁=10%时，试样中产生了明显的塑性变形和体积变形，局部应变集中的区域基本表现为剪胀；最后当ε₁=15%时，塑性区和剪胀区完全贯通，分别如图6（d）和图7（d）所示，形成了较为明显的剪切带。T=-2 ℃时的应变局部化发展规律与T= -5 ℃时基本相似，剪切过程中塑性应变集中的区域，同时伴随着局部化剪胀。可见，所建立的数值框架能够实现对冻土剪切带产生、发展和形成的全过程模拟。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 T=-2 ℃时的冻结砂土等效塑性应变发展规律

Fig. 4 Development law of equivalent plastic strain of frozen sand when T=-2 ℃

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 T=-2 ℃时的冻结砂土体应变发展规律

Fig. 5 Development law of volumetric strain of frozen sand when T=-2 ℃

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 T=-5 ℃时的冻结砂土等效塑性应变发展规律

Fig. 6 Development law of equivalent plastic strain of frozen sand when T=-5 ℃

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 T=-5 ℃时的冻结砂土体应变发展规律

Fig. 7 Development law of volumetric strain of frozen sand when T=-5 ℃

当轴向应变ε₁=15%时，基于模拟和试验获得的试样破坏形态如图8所示。可以看出，T=-2 ℃时模拟获得的剪切带破坏形态与试验结果较为吻合，而T=-5 ℃时模拟获得的破坏区域的位置与试验结果略有偏差，这与设置的端部约束和初始缺陷条件有关，具体关于端部约束和缺陷单元对剪切带破坏形态及破坏位置的影响将在下一节进行讨论分析。基于试验和模拟获得的剪切带倾角和宽度结果如表2所示，当温度分别为-2 ℃和-5 ℃时，本文模拟获得的剪切带倾角分别是51.8°和55.8°，与Yao等^［20］的试验结果50.5°和55.3°大致相同，同时也基本符合摩尔-库伦解45°+φ/2。对于剪切带宽度，需要说明的是，具体宽度值与塑性应变的可视化结果密切相关，科研人员也尝试给出了一些确定剪切带宽度的方法^［28-29］，但目前还没有相对统一的定论。本文主要以图8中观测到的剪切带测量其宽度，两种温度下剪切带宽度与颗粒平均粒径d₅₀的比值分别是18.0和20.4，与融土中剪切带宽度大致是d₅₀的10~20倍的规律类似；此外，随着温度的降低，剪切带宽度逐渐增大，这是由于冻结砂土强度越高，其承载能力越强，土样达到破坏时产生的塑性区也就越大。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 基于试验和模拟获得的冻结砂剪切带形态

Fig. 8 Shear band mode of frozen sand based on test and simulation： test result^［20］（T=-2 ℃）（a）， simulation result （T=-2 ℃）（b）， test result^［20］（T=-5 ℃）（c） and simulation result （T=-5 ℃）（d）

表2 剪切带倾角和宽度

Table 2 Inclination angle and width of the shear band

T/℃	剪切带倾角/（°）			剪切带宽度与d₅₀的比值
T/℃	摩尔-库伦解	试验结果^［20］	模拟结果	试验结果^［20］	模拟结果
-2	52.8	50.5	51.8	17.8	18.0
-5	59.2	55.3	55.8	22.7	20.4

新窗口打开| 下载CSV

2.3　端部约束和初始缺陷的影响

端部约束和初始缺陷使得冻土试样在加载过程中受力不均匀，进而产生不均匀的变形，是产生应变局部化问题的主要因素。以T=-2 ℃为例，说明端部约束和初始缺陷对剪切带发展规律的影响。分别考虑如下5种工况：工况1，不考虑端部约束，不考虑初始缺陷；工况2，不考虑端部约束，考虑非对称初始缺陷；工况3，不考虑端部约束，考虑对称初始缺陷；工况4，考虑端部约束，不考虑初始缺陷；工况5，考虑端部约束，考虑对称初始缺陷。具体不同工况下的端部约束和初始缺陷条件如表3所示。

表3 不同工况下的端部约束和初始缺陷条件

Table 3 End constraints and initial imperfection for different working conditions

工况编号	端部约束	初始缺陷
工况编号	端部约束	单元781	单元800
工况1	×	×	×
工况2	×	×	√
工况3	×	√	√
工况4	√	×	×
工况5	√	√	√

新窗口打开| 下载CSV

图9分别给出了5种工况下冻结砂土的破坏形态。可以看出，剪切带破坏形态主要包括单一型和交叉型，但具体破坏形态及破坏位置需要综合考虑端部约束和初始缺陷的共同作用。其中，当不考虑端部约束和初始缺陷时，试样材料均匀，在加载过程中各单元受力均匀，无局部化的剪切带产生，如图9（a）所示；当不考虑端部约束而只考虑非对称初始缺陷时，由于引入了缺陷单元使得试样材料非均匀，加载过程中受力不均匀，剪切带一般会沿着初始缺陷区域进行开展，具体破坏的位置跟初始缺陷的位置密切相关。工况2只在试样右上角设置了一个缺陷单元，此时局部应变沿着该缺陷单元逐渐开展直到形成宏观的单一型剪切带，如图9（b）所示；当不考虑端部约束而只考虑对称初始缺陷时，此时试样材料非均匀但具有对称性，一般会产生交叉型剪切带，且破坏区域靠近试样缺陷处，如图9（c）所示；当只考虑端部约束而不考虑初始缺陷时，虽然试验材料本身是均匀的，但由于端部约束会对上下边界的水平位移起到限制作用，使得端部区域和中间区域在加载过程中变形非均匀，此时通常会产生交叉型剪切带，并且剪切带位于试样中部，如图9（d）所示；当同时考虑端部约束和对称初始缺陷时，试验材料虽然非均匀但满足对称性，同时端部约束的作用易使得试样产生对称型的剪切带［图9（d）］。因此对称初始缺陷叠加端部约束作用一般也产生对称型剪切带，如图9（e）所示，但剪切带具体的破坏形态及破坏位置需要综合考虑两者的共同影响，例如当设置初始缺陷效果不明显时，其作用会被端部约束作用所覆盖。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 不同工况下的剪切带形态

Fig. 9 Shear band mode for Working Condition 1 （a）， Working Condition 2 （b）， Working Condition 3 （c）， Working Condition 4 （d） and Working Condition 5 （e）

2.4　网格依赖性分析

为了探究有限元网格密度对数值结果的影响，将整个求解域分别划分粗细两种有限元网格，粗网格和细网格的单元数分别为nels=800和nels=3 200。为了保证粗细网格有限元建模的一致性，细网格与粗网格在左上角和右上角相同区域设置了缺陷单元。图10给出了不同有限元网格下的应力-应变曲线，可以看出基于微极理论获得的粗细网格条件下的模拟结果具有较好的收敛性，未出现明显的网格依赖性问题。其中，当T=-2 ℃时，基于微极理论计算获得的不同网格密度条件下的剪切带如图11所示，可以看出粗网格与细网格的剪切带宽度基本一致，说明基于微极理论模拟冻土应变局部化问题时，可以基本克服经典连续体理论的网格依赖性问题。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 不同有限元网格下的应力-应变曲线

Fig. 10 Stress-strain curves for different finite element meshes

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 不同有限元网格下的剪切带形态

Fig. 11 Shear band mode for different finite element meshes

3 结论

在微极连续体理论框架下，通过考虑温度对相关模型参数的影响，建立了模拟冻土应变局部化问题的微极弹塑性数值框架，并应用于冻结砂土平面应变试验中，得到如下主要结论：

（1）基于微极理论所建立的弹塑性框架能够实现对冻土材料剪切带产生、发展和形成的全过程模拟，且模拟获得的应力-应变曲线及剪切带破坏形态与试验结果吻合较好，基本验证了所建立数值框架的合理性与有效性。

（2）通过对不同温度条件下的冻结砂土剪切带发展规律进行研究，发现冻土试样在剪切过程中，塑性应变集中的区域均伴随着局部化剪胀。温度对剪切带倾角和宽度影响较为显著，倾角随温度的降低而逐渐增大，宽度随温度的降低而逐渐加宽，这是由于冻结砂土强度越高，土样达到破坏时产生的塑性区越大。

（3）通过对不同端部约束和初始缺陷条件下的剪切带形态进行分析，发现端部约束和初始缺陷是产生应变局部化问题的主要因素，具体剪切带的破坏形态及破坏位置需要综合考虑端部约束和初始缺陷的共同作用。

（4）通过对不同有限元网格条件下的数值模拟结果进行分析，发现粗网格和细网格获得的应力-应变曲线及剪切带破坏形态大致相同，说明基于微极理论模拟冻土应变局部化问题时，可以基本克服经典连续体理论的网格依赖性问题。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Alshibli

K A

， Batiste

S N

， Sture

Strain localization in sand： plane strain versus triaxial compression

［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2003， 129（6）： 483-494.