基于近场动力学均匀化的混凝土防渗墙等效损伤模型

摘要

混凝土防渗墙在我国覆盖层高土石坝坝基防渗中应用非常广泛。受到巨大水推力、土压力及不均匀沉降作用的混凝土防渗墙可能会出现局部破损，该过程涉及到的尺度差异大、复杂破损精细化描述等问题通常是分析的难点。通过近场动力学多尺度均匀化分析，提出了一种高土石坝混凝土防渗墙破损等效损伤模型。首先给出了该等效损伤模型建立的基本思路，然后基于防渗墙破损的机理和应力状态规律的认识，采用近场动力学多尺度均匀化方法，建立了可用于工程三维分析的高土石坝防渗墙破损的等效损伤模型。该模型具有形式简单、参数少且物理意义明确、拉压统一描述、模型新引入的指数参数可以更好地模拟混凝土损伤发展速度等特点，便于工程应用。通过混凝土材料单轴受拉和受压试验及混凝土偏心三点弯曲梁算例验证了模型的有效性后，将该模型有效应用于某砾石土心墙堆石坝防渗墙的三维破损分析中。

关键词

混凝土防渗墙; 近场动力学; 均匀化多尺度方法; 等效损伤模型

由于能良好地适应基础条件、就地取材、充分利用建筑物开挖渣料、造价较低等优势，土石坝是坝工建设中最有发展前景的坝型之一。防渗结构是高土石坝安全的生命线^［1］，处于覆盖层内的混凝土防渗墙受力条件复杂，在巨大水推力、土压力及不均匀沉降作用下，可能会出现局部破损，对坝基渗透稳定性及坝体安全产生不利的影响。同时，高土石坝-混凝土防渗墙是一个多尺度相互作用的系统，存在着千米级的山体、百米级的大坝、米级的防渗体、厘米级的结构面乃至毫米级的裂缝及更小尺度的土体颗粒等多个尺度，这些不同尺度之间的相互耦联和相互作用共同决定了该复杂系统的性能演化和全寿命周期过程（图1）。因此需要发展多尺度数值分析方法，以更精确地刻画上述多个尺度之间的相互作用和更合理地揭示墙体复杂的破损机制。

图1  高土石坝-防渗结构多尺度相互作用系统

Fig. 1  Multi-scale system of the high embankment dam and its prevention structures

学者们基于现场监测、模型试验和数理分析等各种手段，对高土石坝防渗墙应力变形规律和破损机理等开展了研究，取得了很多进展。在现场监测方面，Rice和Duncan^［2］详细分析了美国境内多座坝体防渗结构的服役情况。在模型试验方面，主要以采用土工离心模型试验的方法为主^［3］。在数理分析方法方面，目前多采用有限元方法，在高土石防渗墙结构的应力变形^［4］、接触模拟^［5］、接头型式^［6］和水力耦合^［7］等多个方面取得了进展。本文主要研究防渗墙结构的破损开裂分析方法。何四清^［8］测定了塑性混凝土的断裂特性，并用弥散裂缝模型对册田水库副坝混凝土防渗墙进行了断裂破坏分析；Rice和Duncan^［9］通过安全度计算分析，指出防渗体最大可能开裂的位置为覆盖层与基岩接触的部位；Yu等^［10］采用塑性损伤模型模拟了墙体的损伤破坏；冯蕊等^［11］基于弥散裂缝模型分析了长河坝坝基廊道的开裂；Cen等^［12］使用细观随机损伤模型分析了面板挤压破损问题等。然而高土石坝中的混凝土防渗墙，因其受力条件复杂，所以破损的多尺度精细化分析方法及据此得到的破损机理规律等研究工作仍需加强。孙伟^［13］基于Peridynamics^［14-16］（PD，即近场动力学）理论，提出了高土石坝防渗墙破损的PD-FEM （Finite element method） 并发式多尺度精细化分析方法。该方法可以实现百米级大坝尺度中毫米级裂缝的精细模拟与扩展演化，较为真实地再现了混凝土防渗墙中拉剪耦合裂缝起裂与扩展形态^［17］。相比弥散裂缝模型或损伤模型等，能够得到更具体准确的分离式的裂缝位置和形态；相比DEM （Discrete element method）-FEM等方法^［18］，则可避免DEM细观参数可能难以确定以及用于实际工程时颗粒大小需要放大等问题。然而，在高土石坝混凝土防渗墙的破损分析中，局部破损计算采用非局部PD模型，虽然该局部区域可以选取得非常小，但在三维情形下，矩阵的稀疏性大大削弱，计算量陡增，面对实际复杂工程仍然是无能为力的。因此，此时并发式多尺度方法不再适用，有必要发展基于均匀化多尺度方法的等效简化分析模型，这种简化模型既能保留前期多尺度分析获得的某些破损信息，又能较为简单地为工程大尺度计算服务，这样就能更好地将多尺度方法应用于实际工程。

多尺度分析方法根据分析对象或问题的不同可大致划分为两类：均匀化递阶式方法和并发式方法^［19］。针对不同的问题，应根据两类方法的特点合理地加以选择。本文采用均匀化多尺度方法（Homogenization/Hierarchical multiscale method）沟通宏细观信息，建立可包含细观破损演化信息在内的宏观损伤模型。Fish等^［20-22］采用双尺度渐进展开法研究复合材料的破损，通过引入某些假设，缩减单胞RVE的自由度，发展了一种降阶多尺度方法（Reduced order multiscale method），其可根据分析问题的特点，引入某些假设，缩减单胞的自由度，这样耦合计算量即大为减少。Ghosh^［23-24］研究提出了另外一种简化多尺度方法的思路，在细观上采用VCFEM（Voronoi cell finite element method）模拟复合材料破损，均匀化后得到宏观上的各向同性或各向异性的损伤模型，直接使用该模型进行宏观破损分析。这种宏观损伤模型，被称为均匀化损伤模型（Homogenization-based continuum damage model）。在土木工程中，也有学者对多尺度等效分析模型进行了研究。Zhang^［25］根据复杂模型和简化模型动力模式等效来进行等效简化模型的参数识别，发展了桥梁隧道工程中大型正交各向异性结构的等效简化建模方法。梁诗雪等^［26］通过多尺度均匀化方法，对含有缺陷的混凝土细观RVE（Representative volume element）进行随机分析，得到了工程实用的混凝土随机损伤演化方程。李阳^［27］根据均匀化多尺度方法得到的刚度矩阵具有各向异性的特点，假设碎石桩复合地基为横观各向同性材料，据此推导了碎石桩复合地基等效简化模型，并给出了等效模型参数的取值方法。

本文借鉴Ghosh建立均匀化损伤模型的方法，采用近场动力学多尺度均匀化方法，建立一种包含破损演化信息在内的等效损伤模型，以供高土石坝工程三维分析使用。

1 近场动力学基本理论

近场动力学（Peridynamics， PD）理论^［14］是一种非局部计算力学理论。这种理论基于积分-微分方程，是对一般用偏微分方程描述的经典连续介质力学理论的一次重大改进，避免了位移导数项，获得了连续-非连续的统一表达，可直接采用键的断裂（bond breakage）来模拟裂缝扩展。它在描述复杂破损问题方面具有巨大的优势，从而可方便地分析防渗墙混凝土结构的破损问题。采用非常规态型（non-ordinary state-based）PD理论^［15］，该理论不再使用键的刚度来与连续介质力学模型等效，而是可以直接嵌入已有的连续介质力学中的本构模型，工程应用较为方便。

NOSBPD的运动方程^［15］为

$$\rho (\mathit{x})\ddot{\mathit{u}}(\mathit{x},t)={\int }_{{\mathrm{H}}_{\mathit{x}}}\left\{\underline{\mathit{T}}[\mathit{x},t]<\mathit{x}\mathrm{\text{'}}-\mathit{x}>-\underline{\mathit{T}}[\mathit{x}\mathrm{\text{'}},t]<\mathit{x}-\mathit{x}\mathrm{\text{'}}>\right\}\mathrm{d}{V}_{\mathit{x}\mathit{\text{'}}}+\mathit{b}(\mathit{x},t)$$

（1）

式中：$\rho$为质量密度；$\mathit{u}$为位移；$\mathit{b}$为给定的体积力；t 为时间；$H_{\mathit{x}}$为点$\mathit{x}$的邻域，它是指初始参考构型中截断半径$\delta$所包围的区域；$\mathit{f}$为对点力函数；$\underline{\mathit{T}}[\mathit{x},t]<\mathit{x}\mathrm{\text{'}}-\mathit{x}>-\underline{\mathit{T}}[\mathit{x}\mathrm{\text{'}},t]<\mathit{x}-\mathit{x}\mathrm{\text{'}}>$表示作用在点$\mathit{x}$上的由点${\mathit{x}}^{\mathrm{\text{'}}}$给的力密度向量；$\underline{\mathit{T}}[\mathit{x},t]$和$\underline{\mathit{T}}[\mathit{x}\mathit{\text{'}},t]$为力向量态。

初始位置参考构形向量态和某一时刻的变形向量态分别为

$$\underline{\mathit{X}}<\mathit{\xi }>=\mathit{\xi }={\mathit{x}}^{\mathrm{\text{'}}}-\mathit{x}$$

（2）

$$\underline{\mathit{Y}}<{\mathit{x}}^{\mathrm{\text{'}}}-\mathit{x}>=\mathit{\eta }+\mathit{\xi }=({\mathit{u}}^{\mathit{\text{'}}}+{\mathit{x}}^{\mathrm{\text{'}}})-(\mathit{u}+\mathit{x})$$

（3）

给定$\underline{\mathit{Y}}$，通过缩减技术，可以获得非局部的变形梯度，如式（4）：

$$\mathit{F}(\mathit{x})=\left[{\int }_{{\mathrm{H}}_{\mathit{x}}}\omega (‖\mathit{\xi }‖)(\underline{\mathit{Y}}(\mathit{\xi })\otimes \mathit{\xi })d{V}_{\mathit{\xi }}\right]\cdot \mathit{K}(\mathit{x})$$

（4）

其中，$\mathit{K}(\mathit{x})$为非局部形变张量。

$$\mathit{K}(\mathit{x})={\left[{\int }_{{\mathrm{H}}_{\mathit{x}}}\omega (‖\mathit{\xi }‖)(\mathit{\xi }\otimes \mathit{\xi })d{V}_{\mathit{\xi }}\right]}^{-\mathrm{1}}$$

（5）

给定某一经典连续介质本构关系，其应变、应力张量均可通过式（4）定义的非局部变形梯度求得。经过推导，力向量态可以用应力和变形梯度表示为

$$\begin{array}{l}\underline{\mathit{T}}[\mathit{x}]<\mathit{x}\mathrm{\text{'}}-\mathit{x}>=\\ \omega (‖\mathit{\xi }‖)\mathrm{d}\mathrm{e}\mathrm{t}(\mathit{F})\cdot {\mathit{F}}^{-\mathrm{1}}\cdot \mathit{\sigma }\cdot \mathit{K}(\mathit{x})\cdot \mathit{\xi }\end{array}$$

（6）

式（6）建立了NOSBPD中的力态与经典连续介质力学理论中的应力、变形梯度（或应变）等物理量之间的关系。即PD键发生变形后，通过式（4）可以计算出非局部变形梯度（也可以得到应变），通过材料本构可以得到应力。将变形梯度（或应变）与应力代入式（6），即可得到PD的力态，从而代入式（1），进行PD运动方程的计算。同时，材料损伤准则体现在影响函数$\omega (‖\mathit{\xi }‖)$中。

2 等效损伤模型建立的基本原理

2.1　总体思路

在采用精细化PD-FEM耦合方法分析防渗墙破损时，随着PD键的逐渐断裂，破损区材料性能不断退化。根据上述NOSBPD方法的基本原理可知，材料在受载过程中键的断裂对结构的作用是通过影响函数来体现的，而影响函数的值又会影响到变形梯度的计算（见式（4）），进而可认为影响到应变的计算，故可将应力表达为

$$\mathit{\sigma }=\mathit{E}:\stackrel{}{\overbrace{(\mathrm{1}-d)\mathit{\epsilon }}}$$

（7）

式中：为损伤变量（为表达方便，此处记其为标量，实际PD-FEM耦合分析时所对应的损伤应是各向异性的张量）。为方便说明式（7），假设某PD点邻域内所有的键都发生了同样的损伤，影响函数均从1.0变为了0.5（图2，图中δ = 2Δx），这样相当于

$$\mathit{\sigma }=\mathit{E}:\stackrel{}{\overbrace{\mathrm{0.5}\mathit{\epsilon }}}$$

（8）

图2  破损区某PD点邻域内键的状态示意图

Fig. 2  Damage status of the bonds in the horizon of a peridynamics material point located in the damage zone

即此时d = 0.5。

因此，精细化PD-FEM耦合分析时，防渗墙破损模拟是通过改造“形函数”来实现的，即某点邻域内的键，如果达到阈值，可以进行删除，这就相当于具有了丰富的“形函数”，可以表达多种变形模式（这与XFEM引入非连续变形的变形模式是类似的），同时键是有方向的，这就使得其依靠单个键的破损准则即能够模拟裂缝的走向，而无需再借助于其他判别标准。而对于基于有限元离散的经典连续介质力学理论，单元内部的变形模式是连续的，不具有表达非连续变形的能力，但可以用损伤变量来表达材料受荷过程中的性能衰退，达到模拟材料破损的目的，即可将式（7）改写成

$$\mathit{\sigma }=\stackrel{}{\overbrace{(\mathrm{1}-d)\mathit{E}}}:\mathit{\epsilon }$$

（9）

所以，首先可通过式（7）表达的PD-FEM耦合分析确定在给定应力路径下混凝土防渗墙破损过程，从而得到损伤变量的演化规律，再代入到式（9）中，即可得到与该精细化方法相等效的简化分析模型。

由精细化PD-FEM耦合分析可知，防渗墙结构的破损是一定应力状态下发生的复杂的材料性能退化过程。同时，防渗墙在特定应力路径下存在拉、压2种破损模式。因此，本文在对防渗墙破损应力路径进行分析和简化后，使用多尺度均匀化方法分析RVE拉压损伤演化过程，最终建立形式简单的损伤模型。

2.2　应力路径分析

（1）受拉破损模式。基于某堰塞坝防渗墙破损的PD-FEM耦合分析^［17］，研究防渗墙在受拉破损模式下的应力路径。该破损模式控制工况在蓄水期。蓄水过程中防渗墙破损区典型点应力路径如图3a所示。防渗墙施工完成后，水位较低时，存在着一定的卸载过程，而后随着墙前水位逐步淹没该点，防渗墙受到的围压和弯矩引起的拉应力同步增大，轴向从受压逐步过渡到受拉状态，最后在拉剪耦合应力作用下发生破损。根据破损区应力状态的分析^［13］，确定破损区围压${\sigma }_{\mathrm{1}}$为0~2.0MPa。

图3  防渗墙应力路径（仅显示了破损时刻附近）

Fig. 3  Stress path of cut-off wall (only near the damage instance is shown)

（2）受压破损模式。基于某砾石土心墙堆石坝覆盖层坝基防渗墙破损的多尺度分析^［13］，研究防渗墙在受压破损模式下的应力路径。该破损模式控制工况在施工期。有研究表明，土石坝填筑过程是一个坝体大小主应力增量比基本不变、各主应力逐渐增大的等应力比加载过程^［28］。防渗墙破损区典型点的应力路径如图3b所示。分析可知，在坝体填筑过程中，墙体的应力路径也基本符合等应力比加载的特点，随着大主应力的持续增大，墙体发生受压破损。根据破损区应力状态的分析^［13］，确定破损区围压${\sigma }_{\mathrm{3}}$为0~2.0MPa。

2.3　均匀化过程

首先取RVE，根据防渗墙破损应力路径的分析，施加相应的拉压和压压混合荷载（图4），同时为了引起不均匀破损，在RVE中引入了随机缺陷，采用PD分析该RVE破损过程。在PD模型中，采用的是基于能量的键的断裂准则^［17］。

图4  细观RVE及其边界条件

Fig. 4  RVE and its boundary conditions

这样，得到RVE的控制方程为

$$\left\{\begin{array}{l}\nabla \cdot {\mathit{\sigma }}_{\mathrm{L}}({\mathit{u}}_{\mathrm{L}})+\mathit{b}=\mathrm{0}\\ {\mathit{u}}_{\mathrm{L}}={\overline{\mathit{u}}}_{\mathrm{L}},\mathrm{在}\mathrm{边}\mathrm{界}{\partial }_u{\Omega }_{\mathrm{L}}\mathrm{上}\\ {\mathit{\sigma }}_{\mathrm{L}}\cdot \mathit{n}={\overline{\mathit{t}}}_{\mathrm{L}},\mathrm{在}\mathrm{边}\mathrm{界}{\partial }_t{\Omega }_{\mathrm{L}}\mathrm{上}\end{array}\right.$$

（10）

式中：下标L表示与RVE有关的量。破损信息由PD表达。

对由式（10）获得的应力场$\mathit{\sigma }$和应变场$\mathit{\epsilon }$（注意这里的非局部应变应是所有键都完好情形下计算出来的，即后处理时，对于断裂的键，应将其恢复为完整的键再进行应变计算）进行均匀化，如式（11）、（12）：

$$\overline{\mathit{\sigma }}=\frac{\mathrm{1}}{V}{\int }_V{\mathit{\sigma }}_{ij,L}^{}\mathrm{d}V$$

（11）

$$\overline{\mathit{\epsilon }}=\frac{\mathrm{1}}{V}{\int }_V{\mathit{\epsilon }}_{ij,L}^{}\mathrm{d}V$$

（12）

式中：$V$为RVE所占有的体积。在获得$\overline{\mathit{\sigma }}$和$\overline{\mathit{\epsilon }}$后，基于应变等价假设，按照式（13）计算损伤变量：

$$d=\frac{E-E\mathrm{\text{'}}}{E}$$

（13）

式中：$d$为损伤变量；$E$为无损时的弹性模量；$E\mathrm{\text{'}}$为损伤后的弹性模量。

PD方法可以自发地判断裂缝走向，因此得到的损伤应是各向异性的张量，但这里简单地通过一个标量d来定义损伤。即通过多尺度均匀化分析得到损伤因子$d$与总应变$\epsilon$的关系，从而建立包含破损信息在内的宏观各向同性损伤模型。

得到损伤变量的演化规律后，即可采用经典连续介质力学理论描述宏观尺度，并采用常规有限元方法进行离散。这样，得到总的宏观控制方程为

$$\left\{\begin{array}{l}\nabla \cdot {\mathit{\sigma }}_{{}^G}+\mathit{b}=\mathrm{0}\\ {\mathit{\sigma }}_{{}^G}={\mathit{\sigma }}_{{}^G}(\mathit{\epsilon }(\mathit{u}),\mathit{\xi },d)\\ {\mathit{\epsilon }}_{{}^G}={\mathit{\epsilon }}_{{}^G}(\mathit{u})=\frac{\mathrm{1}}{\mathrm{2}}(\nabla {\mathit{u}}_{{}^G}+\nabla {\mathit{u}}_{{{}^{G,}}^T})\\ {\mathit{u}}_{{}^G}={\overline{\mathit{u}}}_G,\mathrm{在}\mathrm{边}\mathrm{界}{\partial }_u{\mathrm{\Omega }}_G\mathrm{上}\\ {\mathit{\sigma }}_{{}^G}\cdot \mathit{n}={\overline{\mathit{t}}}_G,\mathrm{在}\mathrm{边}\mathrm{界}{\partial }_t{\mathrm{\Omega }}_G\mathrm{上}\end{array}\right.$$

（14）

式中：$\mathit{\xi }$为内变量；G表示与宏观有关的量。破损信息由损伤变量$d$表达。

3 等效损伤模型

3.1　损伤演化规律

对于普通混凝土防渗墙，按照泵送浇筑的要求，一般采用二级配混凝土，参考水工混凝土试验规程对于二级配混凝土试样尺寸的规定，选取RVE大小为$\mathrm{150}\mathrm{m}\mathrm{m}\times \mathrm{150}\mathrm{m}\mathrm{m}$。混凝土变形破损计算参数如表1所示。

表1  防渗墙RVE变形破损计算参数

Tab. 1  Computational parameters for the concrete cut-off wall

$\mathit{E}$/GPa	$\mathit{\upsilon }$	${\mathit{f}}_{\mathrm{t}}$/MPa	${\mathit{G}}_{\mathrm{I}}$/（N·m-1）	${\mathit{f}}_{\mathrm{c}}$/MPa	${\mathit{G}}_{\mathrm{I}\mathrm{I}}$/（N·m-1）
30.0	0.2	2.01	135	20.1	2 000

（1）受拉破损模式。根据前面应力路径的分析，这里研究${\sigma }_{\mathrm{1}}\mathrm{为}\mathrm{0}\mathrm{、}\mathrm{1}\mathrm{、}\mathrm{2}\mathrm{M}\mathrm{P}\mathrm{a}$不同围压情形下防渗墙RVE在拉压荷载作用下的破损规律。边界条件如图4a所示。最终得到不同围压下，承受拉压荷载作用的防渗墙RVE的损伤曲线如图5a所示。由分析可知，由于副轴方向压应力与抗压强度的比值比较小，破坏仍主要由主轴方向的拉应力控制，因此围压的增大对峰值强度和损伤曲线的影响较小。同时，损伤曲线均大致呈现外凸形状，即在损伤值0.8之前，损伤发展较快，在其之后则发展较慢，逐步趋近于完全损伤。

图5  不同围压荷载作用下RVE损伤曲线

Fig. 5  Damage curves at different confining pressure loads of the RVE

（2）受压破损模式。根据前面应力路径的分析，这里研究${\sigma }_{\mathrm{3}}\mathrm{为}\mathrm{0}\mathrm{、}\mathrm{1}\mathrm{、}\mathrm{2}\mathrm{M}\mathrm{P}\mathrm{a}$不同围压情形下防渗墙RVE在压压荷载作用下的破损规律。边界条件如图4b所示。最终得到不同围压下，承受压压荷载作用的防渗墙RVE损伤曲线如图5b所示。由分析可知，由于副轴/主轴应力比较小，所以围压对受压损伤发展影响比较小（特别是在损伤发展较快的初始阶段），并且仍然具有在损伤值小于0.8的曲线段发展较快而后发展趋缓的特点。

3.2　拉压统一损伤模型

由前面的分析可知，在防渗墙破损对应的应力路径下，围压对受拉和受压损伤发展的影响较小，同时受拉和受压损伤曲线均具有指数型的发展形式，因此这里以单轴受拉和单轴受压破损过程为研究对象，参考两类形式比较简单混凝土静动力损伤模型即双线性模型^［29］和指数型模型^［30］后，考虑多尺度均匀化得到的受拉和受压损伤曲线的特点，建立了如式（15）形式简单的拉压统一损伤模型：

$$d_{\mathrm{t}/\mathrm{c}}=\left\{\begin{array}{l}\mathrm{0}{,}_{}\epsilon \le {\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}\\ \mathrm{1}-{(\frac{{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}}{\epsilon })}^{n_{\mathrm{t}/\mathrm{c}}}{\mathrm{e}}^{(-\frac{\epsilon -{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}}{{\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}^{}-{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}})}{,}_{}\epsilon >{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}\end{array}\right.$$

（15）

式中：$\epsilon$为应变；${\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}$为损伤门槛应变；${\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}^{}$为接近完全损伤对应的应变；$n_{\mathrm{t}/\mathrm{c}}$为本文添加的参数，用来控制损伤发展的快慢。

使用式（15）拟合损伤曲线，受拉损伤取${\epsilon }_{\mathrm{0},\mathrm{t}}=\mathrm{6.67}\times {\mathrm{10}}^{-\mathrm{5}},{\epsilon }_{\mathrm{u},\mathrm{t}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}},n_{\mathrm{t}}=\mathrm{2.5}$，受压损伤取${\epsilon }_{\mathrm{0},\mathrm{c}}=\mathrm{6.67}\times {\mathrm{10}}^{-\mathrm{4}},{\epsilon }_{\mathrm{u},\mathrm{c}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}},n_{\mathrm{c}}=\mathrm{0.8}$。由图6可知：本文模型可以通过调节参数$n_{t/c}$来模拟相应的损伤发展速度，对受拉损伤模拟较好，而双线性模型是假设应力应变软化段为直线型而得到的简化损伤模型，其和指数型模型都有可能会低估混凝土损伤发展速度。为了使模型形式简单且参数少，本文没有对受压损伤进行更精细的描述，因此模型对受压损伤模拟有一定的偏差，但选择合适的参数后，模型总体上仍能近似模拟曲线的大致形状和演化过程。

图6  损伤因子演化规律拟合结果

Fig. 6  Fitting results of evolution laws of damage factors

3.3　模型参数的确定

模型拉压损伤$d_{t/c}$各具有3个参数：${\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}$、${\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}^{}$和$n_{\mathrm{t}/\mathrm{c}}$。

（1）${\epsilon }_{\mathrm{0},t/c}^{}$为损伤门槛应变。这里与Mazars模型^［31］一样，假设峰值前无损伤，近似取${\epsilon }_{\mathrm{0},t}=f_{\mathrm{t}}/E_{\mathrm{0}}$和${\epsilon }_{\mathrm{0},\mathrm{c}}=f_{\mathrm{c}}/E_{\mathrm{0}}$，其中$f_{\mathrm{t}}$为单轴抗拉强度，$f_{\mathrm{c}}$为单轴抗压强度，$E_{\mathrm{0}}$为初始弹性模量。${\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}$可通过混凝土抗拉和抗压强度试验确定。

（2）${\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}^{}$为接近完全损伤时的应变。它们可取为混凝土极限拉压应变，并可通过混凝土极限拉伸和混凝土柱承压试验确定。定义极限应变系数$\eta ={\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}^{}/{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}^{}$，若无试验数据，$\eta$可取10左右。

（3）$n_{\mathrm{t}/\mathrm{c}}$为控制损伤发展快慢的参数。$n_{\mathrm{t}/\mathrm{c}}$取值越大，则损伤越快。$n_{\mathrm{t}/\mathrm{c}}$可通过试验测定混凝土单轴受拉或受压应力应变全曲线，按照式（13）计算损伤值$d_{\mathrm{t}/\mathrm{c}}$，然后根据试验测点进行确定。假设有试验测点$({\epsilon }_{\mathrm{1},\mathrm{t}/\mathrm{c}},d_{\mathrm{1},\mathrm{t}/\mathrm{c}})$，其中$d_{\mathrm{1},\mathrm{t}/\mathrm{c}}\ne \mathrm{0},d_{\mathrm{1},\mathrm{t}/\mathrm{c}}\ne \mathrm{1}$，则

$$n_{\mathrm{t}/\mathrm{c}}=\left(\mathrm{l}\mathrm{n}(\mathrm{1}-d_{\mathrm{1},\mathrm{t}/\mathrm{c}})+\frac{{\epsilon }_{\mathrm{1},\mathrm{t}/\mathrm{c}}-{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}}{{\epsilon }_{\mathrm{u},\mathrm{t}/\mathrm{c}}-{\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}}\right)\frac{\mathrm{1}}{\mathrm{l}\mathrm{n}({\epsilon }_{\mathrm{0},\mathrm{t}/\mathrm{c}}/{\epsilon }_{\mathrm{1},\mathrm{t}/\mathrm{c}})}$$

（16）

在实际操作中，可选取多对试验点，最后取平均值加以确定$n_{\mathrm{t}/\mathrm{c}}$值，甚至还可以在曲线的不同阶段选取各自的测点，从而确定各自不同的$n_{\mathrm{t}/\mathrm{c}}$，这样可以更好地与试验曲线相匹配。若无试验数据，对于常见的刚性混凝土防渗墙，通过相应的多尺度均匀化计算分析，建议$n_{\mathrm{t}}$取2.5，而$n_{\mathrm{c}}$取0.8。

4 模型的验证

4.1　混凝土单轴试验

首先采用文献中的材料试验^［32］验证本文损伤模型模拟混凝土单轴受拉和受压应力-应变曲线及其损伤发展过程的有效性。混凝土材料参数为：$E=\mathrm{31}\mathrm{G}\mathrm{P}\mathrm{a}$，$f_{\mathrm{t}}=\mathrm{3.48}\mathrm{M}\mathrm{P}\mathrm{a}$，$f_{\mathrm{c}}=\mathrm{27.6}\mathrm{M}\mathrm{P}\mathrm{a}$。损伤模型参数根据试验数据确定为：${\epsilon }_{\mathrm{0},\mathrm{t}}=\mathrm{1.12}\times {\mathrm{10}}^{-\mathrm{4}},{\epsilon }_{\mathrm{u},\mathrm{t}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}},n_{\mathrm{t}}=\mathrm{2.5}$和${\epsilon }_{\mathrm{0},\mathrm{c}}=\mathrm{8.90}\times {\mathrm{10}}^{-\mathrm{4}},{\epsilon }_{\mathrm{u},\mathrm{c}}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{3}},n_{\mathrm{c}}=\mathrm{0.8}$。模型模拟效果如图7和图8所示，其中试验损伤值是根据试验数据采用式（13）计算而得到的。

图7  混凝土单轴受拉模拟

Fig. 7  Simulation results of uniaxial tension test of concrete

图8  混凝土单轴受压模拟

Fig. 8  Simulation results of uniaxial compression test of concrete

模拟结果表明，对于混凝土单轴受拉，本文模型模拟效果较好，而双线性模型（${\epsilon }_{\mathrm{u},\mathrm{t}}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$）和指数型模型（${\epsilon }_{\mathrm{u},\mathrm{t}}=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{4}},n_{\mathrm{t}}=\mathrm{1.0}$）则低估了损伤发展速度，模拟效果不佳；对于混凝土单轴受压，本文模型假设峰值前无损伤，这导致模拟结果与试验值存在一定的偏差，更精确地模拟可以考虑调整混凝土受压损伤门槛值，但这可能需要更多的模型参数，不便于工程应用。

总体上来说，本文模型新引入了指数参数$n_{\mathrm{t}/\mathrm{c}}$，相比于指数型或双线性模型，可较好地模拟混凝土材料单轴受拉损伤发展过程，同时也可近似地模拟混凝土材料单轴受压损伤发展过程。

4.2　偏心三点弯曲梁

采用混凝土偏心三点弯曲梁^［33］算例来验证等效模型在混凝土结构破损模拟方面的有效性。混凝土计算参数为：$E=\mathrm{38}\mathrm{G}\mathrm{P}\mathrm{a}$，$\upsilon =\mathrm{0.2}$，$f_{\mathrm{t}}=\mathrm{3.0}\mathrm{M}\mathrm{P}\mathrm{a}$，$G_{\mathrm{I}}=\mathrm{69}\mathrm{N}·{\mathrm{m}}^{-\mathrm{1}}$，损伤模型参数为：${\epsilon }_{\mathrm{0},\mathrm{t}}=\mathrm{7.89}\times {\mathrm{10}}^{-\mathrm{5}},{\epsilon }_{\mathrm{u},\mathrm{t}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}},n_{\mathrm{t}}=\mathrm{2.5}$（需要根据有限元单元尺寸大小，进行相应的正则化，以保持断裂能的恒定^［34］）。PD-FEM耦合模型中PD的离散参数为：$\mathrm{\Delta }x=\mathrm{1}\mathrm{m}\mathrm{m},\delta =\mathrm{3}\mathrm{\Delta }x$。

图9为试验布置和几何尺寸。图10为2种方法的计算模型。图11为最终裂缝形态。图12为2种方法得到的加载点承载力-位移曲线与试验测值的比较， 灰色带阴影部分为试验测得的包络线。由分析可知，在破损分布方面，损伤模型计算的损伤区域由于弥散效应，相比多尺度或者试验结果，稍微偏大；在承载力位移曲线方面，峰值及峰后行为基本都在试验测得的包络区间内，吻合较好。

图9  混凝土偏心三点弯曲梁试验（单位：mm）

Fig. 9  Concrete eccentric three-point bending beam test (unit: mm)

图10  计算模型

Fig. 10  Computational models for different methods

图11  变形放大100倍时的最终裂缝形态

Fig. 11  Final crack patterns (at a deformation magnification factor of 100)

图12  加载点荷载-位移曲线

Fig. 12  Curves of bearing capacities versus displacements

因此，本文建立的等效模型可以在避免多尺度复杂计算的前提下获得较好的混凝土结构破损模拟效果。

5 模型的工程应用

以某砾石土心墙堆石坝为工程背景，并进行一定的简化，使用本文等效模型进行该堆石坝防渗墙三维破损分析。该心墙堆石坝坝高123m，河床覆盖层（最大厚度72.4m）设置一道厚度为1.2m的全封闭式混凝土防渗墙防渗^［35］。计算模型如图13所示。计算参数^［36］如表2所列。另外，防渗墙损伤模型参数取自图6。

图13  某覆盖层上的砾石土心墙堆石坝工程计算模型

Fig. 13  Computational model for a high gravel earth core rockfill dam built on overburden

表2  计算参数

Tab. 2  Computational parameters

材料	密度/（g ·cm-3）	E/GPa	v	c/kPa	φ/ （°）	Δφ/ （°）	K	n	Rf	G	F	D	Kur
堆石料	2.25			0	54.0	11.0	1 200	0.38	0.70	0.3	0.03	3.0	2 200
接触黏土	1.60			40.0	14.0	0	150	0.20	0.86	0.45	0.10	2.0	300
心墙料	1.97			35.0	33.0	0	447	0.40	0.75	0.39	0.05	1.9	900
覆盖层	1.42			0	39.0	0	1 000	0.38	0.65	0.32	0.02	4.0	2 000
防渗墙	2.40	30.0	0.20
基岩	2.40	10.0	0.33

图14给出了基于等效损伤模型的墙体位移、应力和损伤分布。由分析可知，等效模型能够有效地模拟高土石坝坝基防渗墙的三维破损状况，其中受拉损伤主要分布在防渗墙顶部扩大端和防渗墙两岸岸坡约束区域，而河床中下部区域则出现较大范围的受压损伤。

图14  基于等效损伤模型的墙体位移、应力和损伤分布

Fig. 14  Distributions of displacements, stress, and damage of cut-off wall using the method proposed

6 结语

基于防渗墙破损的机理和应力状态规律的认识，采用近场动力学多尺度均匀化方法，建立了可用于工程三维分析的高土石坝防渗墙破损的等效模型，最后完成了模型的验证和工程应用。同时，该等效损伤模型建立的思路也可供其他类似工程参考。

近场动力学多尺度均匀化方法建立的拉压统一损伤模型，形式简单、参数少且具有明确的物理意义。模型新引入了指数参数$n_{\mathrm{t}/\mathrm{c}}$，相比于指数型或双线性模型，可以更好地模拟混凝土材料单轴受拉损伤发展过程，同时也可近似地模拟单轴受压损伤发展过程。采用偏心三点弯曲梁算例验证了等效模型可以在避免多尺度复杂计算的前提下获得较好的混凝土结构破损模拟效果，破损分布和结构承载力曲线计算结果与试验或多尺度方法较为一致。采用该模型进行了某砾石土心墙堆石坝工程的三维计算分析，结果表明建立的等效模型能够有效地模拟高土石坝坝基防渗墙的三维破损状况，其中受拉损伤主要分布在防渗墙顶部扩大端和防渗墙两岸岸坡约束区域，而河床中下部区域则出现较大范围的受压损伤。

作者贡献声明

孙 伟：完成论文主要计算和分析工作，撰写论文。

包世诚：整理论文初稿。

张 嘎：提出选题，设计论文框架。