冻土–复合土工布界面复合幂–指数本构模型研究

Composite Power-Exponential Constitutive Modelling on Frozen Soil-Composite Geotextile Interface

摘要: 土工合成材料广泛应用于各种工程中，由于土体与土工合成材料的物理力学特性差异明显，且土体易受水热条件影响，不同条件下界面的剪切力学特性存在明显的差异，界面的应力应变曲线也存在较大的形态差异。为准确描述冻土–复合土工布界面剪切特性，基于在不同水热循环工况下进行直剪试验得到的试验数据，提出了一种改进的复合幂–指数界面本构模型，并对得到的模型参数进行了分析和讨论。结果表明该模型能准确地反映不同类型的软化型和硬化型曲线，模型中不同参数有其对应的明确物理意义。针对不同复杂工况条件下的应力应变关系进行有效的拟合，并结合具体试验工况给出典型的有代表性的拟合函数关系式，对冻土区基础工程建设具有指导和参考意义。

Abstract: Geosynthetics are widely used in a variety of projects, due to the obvious differences in the physical and mechanical properties of the soil and geosynthetics, and the soil is susceptible to hydrothermal conditions, there are obvious differences in the shear mechanical properties of the interface under different conditions, and there are also large differences in the shape of the stress-strain curve at the interface. To accurately describe the shear characteristics of the frozen soil-composite geotextile interface, an improved composite power-exponential interfacial constitutive model is proposed based on the test data obtained from direct shear tests carried out under different hydrothermal cycling conditions, and the obtained model parameters are analyzed and discussed. The results show that the model can accurately reflect different types of softening and hardening curves, and the different parameters in the model have their corresponding clear physical meanings. The stress-strain relationship under different complex working conditions is fitted effectively, and a typical and representative fitted functional relationship is given in combination with specific test conditions, which is of guiding and reference significance for the construction of infrastructure projects in frozen soil areas.

1. 引言

界面通常是由两种不同材料接触组成的不连续面，如岩石节理面、新旧混凝土界面及结构和土体界面。界面作为两种材料之间的薄弱环节，在外力作用下常首先出现破坏，从而引起一系列工程事故和灾害[1]。在冻土工程中，冻土与结构物接触界面的剪切力学特性相比融化环境下更加复杂，因为不同的水、热、力条件下界面水与胶结冰之间的动态变化使界面土体力学性质和土体与结构接触状态不断变化[2] [3]。在界面剪切过程中的不同阶段，冻土与结构物的接触力类型随着剪切阶段的进行不断发生变化。典型冻土–结构界面剪应力–剪切位移关系如图1所示。在剪切刚开始时，界面处的冰晶体、土体颗粒及接触状态是比较完整的，此时试样的变形近似为弹性变形阶段；随着剪切位移的增大，曲线斜率逐渐降低，此过程中冰晶体逐渐破坏，界面处的土颗粒、冰颗粒和部分冰包裹体开始与土工布产生摩擦作用；剪切强度达到峰值强度之后，界面处的胶结冰出现了大量的破碎，界面胶结力迅速下降，界面剪应力值也随之下降，此时破碎的冰土混合物与土工布之间的摩擦力是界面剪应力的主要组成部分，同时界面中也存在着一些新生成的冰晶颗粒；随着界面冰胶结力完全消失，界面剪切强度完全由冰土混合物和土工布之间的摩擦作用提供，剪应力值趋于稳定，此时的剪切强度称为界面的残余强度[4]。

界面的剪切行为受法向压力、材料特性、固结历史、表面粗糙度及各种耦合因素的影响[5]。许多学者通过剪切试验、数值模拟和理论分析对结构与土体界面的剪切行为开展了研究[6] [7]。由于界面剪切本构模型可以定量模拟界面的力学行为特征，因此其是接触面研究中的重要内容。在过去的研究中，界面本构模型可分为三类：基于试验现象的经验模型，统计损伤模型和弹塑性模型[8]-[13]。其中经验模型在实验室或原位试验结果的基础上建立，具有易用性和灵活性的优势。由于冻土特殊的水热力属性，其界面剪切力学特性也表现出复杂的强硬化、弱硬化、弱软化及软化特性，通过弹塑性理论或统计损伤理论将使得模型中参数过多，给其工程应用带来诸多不便。而经验本构模型通常参数较少，数据易于获得，能够很好地描述不同变化规律的界面剪切力学行为[14] [15] [16]。如对于强软化型曲线，使用驼峰曲线能够实现较好的拟和[17]，对于强硬化型曲线，用幂指数模型的拟合效果更好[18] [19]。王丽琴[20]等通过分析幂函数与指数函数，提出了一种新的非线性复合幂–指数模型，并给出了模型参数的确定方法，经测试该模型可描述各种型式的应力–应变曲线及不同形态的体变曲线，具有广泛的适用性和准确性。

对于冻结状态下的不同结构与土体界面，国内外学者基于不同理论建立了相应的界面本构模型。陈晓[21]等从界面细观力学特征出发，将剪切过程中的界面分为损伤和非损伤两部分，基于均匀化方法建立了冻土与构筑物界面的一维本构模型，拟合结果表明该模型对不同类型曲线都有较好的拟合效果。陈良致[22]等采用应变直剪仪开展了青藏冻结粉土–玻璃钢基础界面直剪试验，基于试验得到的接触面强度变化规律和受力变形特性建立了耦合温度效应的应力位移关系预测模型，拟合结果表明该模型能够较好的描述接触面的力学特性。He [23]，Zhang [24]，Shi [25]等学者开展了一系列界面直剪试验研究了冻土与不同材料的界面剪切特性，基于试验得到的界面强度变化规律和受力变形特性分别建立了不同形式的本构模型，参考原理包括扰动状态理论、二元介质理论和弹塑性理论，模型验证结果表明不同模型均可有效反映界面的力学行为。

综上所述，冻土与结构接触面的本构模型是分析冻土区工程结构与冻土相互作用、评价工程安全稳定性的基础和关键。土体与结构物界面的应力应变曲线反映了其受到荷载后的变形规律，准确描述界面的应力应变行为是建立本构模型需要考虑到核心问题，建立冻土与结构界面有效的本构模型对寒区工程的设计和建设有重要的指导意义。

Figure 1. Shear stress-shear displacement relationship at typical frozen soil-structure interface

图1. 典型冻土–结构界面剪应力–剪切位移关系

2. 复合幂–指数本构模型建立

2.1. 应力应变曲线形态

在讨论本构模型之前已经进行了一系列冻土–界面直剪试验，得到了不同水热循环条件下的试验数据。由于试验工况的不同，得到的试样应力应变曲线形态是多种多样的，包括稳定型、硬化型和软化型，其中硬化型又包括弱硬化型和强硬化型，软化型也同样包括弱软化型和强软化型，采用刘祖典[19]对曲线形态的划分标准，粉质黏土–复合土工布界面剪切试验中不同的应力位移曲线类型大致如图2所示。通过对不同的经验型模型进行分析，发现采用传统的单一形式的数学模型很难对冻土–土工布界面形态各异的应力位移曲线进行良好的拟合[26] [27] [28] [29]。因此，本论文采用一种改进的幂–指数函数模型对试验中的数据进行拟合，以下简称CPE模型(composite power exponential model)，拟合结果证实了复合函数模型的有效性。

Figure 2. Stress-displacement relationship at the interface of silty clay-composite geotextile

图2. 粉质黏土–复合土工布界面应力–位移关系

2.2. 模型函数式的建立

观察图2中不同形式的应力位移曲线，可以发现它们有两个共同特点：1) 过0点，2) 峰值点或无穷远处的函数值有界。对于不同条件下的试验结果，峰值前的曲线形态通常是相似的，包含初始的线性增加阶段及随后的非线性屈服阶段。结合图2，对于脆性破坏曲线的峰值前阶段，曲线形态与幂函数形式比较接近，取幂函数形式为：

$g\left(x\right)=a{x}^m-k$ (1)

式中，x为剪切位移，a，k为常数，假设x ≥ 0，m > 0，则a > 0时，g(x)为单调增函数，a < 0时，g(x)为单调减函数，a = 0时，g(x) = −k。

当界面剪切应力超过峰值强度后，不同的试验条件下峰值后阶段剪切应力发展趋势显著不同，包括剪切应力继续增大、剪切应力稳定和剪切应力衰减等。结合试验现象(图2)，对于塑性破坏曲线和脆性破坏曲线的峰值后阶段，曲线的形态与指数函数的形式比较接近，取指数函数的形式为：

$h\left(x\right)=e^{-b{x}^n}$ (2)

式中，b，n为常数，假设b，n > 0，则h(x)是单调递减函数且有界，当x = 0时，h(x) = 1，当x→+∞时，h(x)→0。

将式(1)与式(2)相乘得：

$p\left(x\right)=\left(a{x}^m-k\right)e^{-b{x}^n}$ (3)

容易观察到函数式(3)没有过坐标原点，因此在式中补充k值，使函数图形过原点。得到新的复合幂–指数模型表达式如式(4)所示[20]：

$f\left(x\right)=\left(a{x}^m-k\right)e^{-b{x}^n}+k$ (4)

当x→+∞时，f(x)→0，因此该模型符合应力位移曲线的两个特点。将式(4)乘以${\tau }_a$ 进行量纲统一，得到应力位移曲线的表达式如式(5)所示：

$\tau ={\tau }_a\left[\left(a{x}^m-k\right)e^{-b{x}^n}+k\right]$ (5)

式中，$\tau$ 为剪应力，${\tau }_a$ 为无量纲量，其值与平均剪应力值大小相等，考虑到本文试验数据值大多在60 kPa左右，因此将${\tau }_a$ 取为60，x为剪切位移，a，b，k，m，n均为模型参数，并且满足b，k ≥ 0，m，n > 0。

当x→+∞时，$\tau$ →$k{\tau }_a$ ，当应力位移曲线为软化型时，$k={\tau }_r/{\tau }_a$ ，${\tau }_r$ 的值为残余强度，这也说明k的值与界面的残余强度有关。当界面的应力位移曲线为硬化型时，$k={\tau }_u/{\tau }_a=\left({\sigma }_1-{\sigma }_3\right)/{\tau }_a$ ，其中，${\tau }_u$ 为界面的极限强度，${\sigma }_1$ 和${\sigma }_3$ 分别为该应力状态中的最大和最小主应力。残余强度${\tau }_r$ 和极限强度${\tau }_u$ 均可从应力位移曲线图中得出。值得一提的是，利用非线性拟合的方法可以快速的得到不同的模型参数值。

3. 结果和讨论

3.1. 冻融循环条件下曲线拟合分析

Figure 3. Interfacial shear characteristic curves under different water contents and freeze-thaw cycling conditions (σ = 50 kPa, T = −10℃)

图3. 不同含水率和冻融循环条件下界面剪切特性曲线(σ = 50 kPa, T = −10℃)

图3为不同含水率和冻融循环次数条件下应力位移曲线试验结果和拟合结果对比[4]。从图中可以看到，图3(b)中有弱硬化形态和稳定塑性形态的应力位移曲线，图3(d)中有不同软化程度的应力位移曲线，从拟合结果可以看出来，该模型能够较好地对不同形态的应力位移曲线进行拟合。表1为图3中拟合曲线的模型参数，其中R²为相关系数的平方值，从表1中可以看到，不同工况条件下的R²值均大于0.981，表明该模型值能够准确的描述不同类型的应力位移曲线。通过观察形态明显的软化和硬化曲线可以发现，脆性破坏型拟合曲线的a值都为正数，塑性破坏型拟合曲线的a值都为负数，除此之外，通过观察到为塑性稳定型的曲线可能有微弱的软化趋势或硬化趋势，从而导致部分塑性稳定型曲线a值出现正负不一的情况。b值与CPE模型函数表达式中的指数部分有关，能够更好的描述塑性破坏型曲线，因此硬化型曲线的b值大于软化型曲线的b值。k值与界面的残余强度有关，残余强度越大则k值越大，m和n值的大小与曲线的形态和应力位移关系有关。对于软化型曲线，当参数m值大于1，且应力位移曲线的下凹程度随着m值的增大而增大，对应土体裂隙发育情况，当参数m值小于1，且应力位移曲线的上凸程度随着m值的增大而减小，说明此时土体裂隙未发育[30]。

Table 1. Parameter values of CPE model fit for curves with different water contents and freeze-thaw cycle conditions

表1. 不同含水率和冻融循环条件下曲线的CPE模型拟合参数值

含水率/%	冻融循环次数/次	a	b	k	m	n	R^2
10	0	0.19	0.00	0.43	1.11	6.36	0.981
	5	−0.49	1.45	0.47	2.11	1.23	0.998
	10	0.02	1.15	0.44	1.11	1.10	0.996
	20	−0.29	1.29	0.48	3.02	1.44	0.998
13	0	−0.38	1.19	0.51	2.11	0.92	0.996
	5	−0.07	0.88	0.50	1.75	1.12	0.998
	10	−0.65	1.20	0.63	3.12	1.58	0.996
	20	−0.29	1.35	0.57	1.48	1.17	0.993
16	0	0.19	0.24	0.49	0.28	1.27	0.992
	5	0.25	0.11	0.47	0.63	4.88	0.997
	10	0.28	0.13	0.69	1.07	2.26	0.996
	20	−1.29	2.32	0.46	1.92	1.34	0.996
19	0	−0.42	1.29	0.49	3.25	1.46	0.982
	5	0.17	0.39	0.42	1.66	1.29	0.994
	10	0.03	0.86	0.66	5.45	1.38	0.994
	20	0.27	0.12	0.59	0.86	5.28	0.995

3.2. 湿干循环条件下曲线拟合分析

图4为不同法向应力和湿干循环条件下界面应力位移曲线试验结果和拟合结果对比[31]。从图中可以看到，低法向应力状态时剪应力值较小，曲线的拟合度也更高，高法向应力状态时界面的剪应力值更大，拟合曲线的拟合效果不如低法向应力状态，表现为虽然拟合曲线形态接近，但存在部分点不能完美对应，这与量纲统一时的选值有关。表2为图4所对应的湿干循环条件下应力位移曲线的CPE模型拟合参数，从表中可以看到，R²的值均大于0.979，说明模型拟合效果良好。通过观察a值发现，没有经历湿干循环的试样，拟合曲线的a值都为负，为硬化型，由前面的分析知这是由于摩擦作用不断累积造成的，经历若干次湿干循环后，应力位移曲线呈软化型，拟合曲线的a值都为正，这与冰胶结力的消失有关。由于k值与界面的残余强度有关，因此随着法向应力的增大，k值也在不断增大，从25 kPa到100 kPa，k值均值的增长幅度为116.5%，与实际残余强度值相对应。

Figure 4. Interfacial shear characteristic curves under different normal stress and wet-dry cycle conditions (T = −6℃)

图4. 不同法向应力和湿干循环条件下界面剪切特性曲线(T = −6℃)

3.3. 湿干冻融循环条件下曲线拟合分析

图5为不同法向应力和湿干冻融循环次数条件下应力位移曲线试验结果和拟合结果对比。从图中可以看到，k值的变化与前面湿干循环时的分析大致相同，即随着法向应力的增大而增大，法向应力从25 kPa增加到100 kPa，k值均值的增长幅度为197.4%，与实际残余强度值相对应。从图中也可以看出低法向应力时的拟合结果优于高法向应力条件。表3所示为湿干冻融条件下CPE模型的参数拟合值，从表中可以看到，相关系数R²的拟合值均在0.954之上，其中在75 kPa条件下，不同循环次数条件下的R²值均大于0.997，表明拟合结果非常准确。

Table 2. Parameter values of CPE model fit for curves with different normal stress and wet-dry cycle conditions

表2. 不同法向应力和湿干循环条件下曲线的CPE模型拟合参数值

法向应力/kPa	湿干循环次数/次	a	b	k	m	n	R^2
25	0	−0.12	1.06	0.35	3.01	1.25	0.999
	1	0.25	1.03	0.29	1.80	0.97	0.993
	3	0.24	0.46	0.25	0.96	1.08	0.998
	5	0.32	1.05	0.38	2.06	0.88	0.995
	10	0.20	1.32	0.27	3.15	0.91	0.992
50	0	−0.58	1.73	0.50	1.99	0.92	0.998
	1	0.06	0.78	0.45	3.35	1.30	0.997
	3	0.27	0.08	0.54	0.92	2.32	0.997
	5	0.48	0.63	0.52	1.74	1.20	0.996
	10	0.47	0.77	0.57	1.98	1.40	0.979
75	0	−0.17	0.95	0.65	2.20	1.26	0.999
	1	0.05	0.97	0.62	3.65	1.12	0.989
	3	0.32	0.16	0.75	1.06	1.70	0.991
	5	0.31	0.49	0.73	1.92	1.23	0.991
	10	0.37	0.63	0.72	2.09	1.21	0.997
100	0	−0.18	0.69	0.87	2.23	1.09	0.999
	1	−0.34	1.04	0.84	3.31	1.67	0.988
	3	0.23	0.84	0.87	3.14	1.27	0.991
	5	0.29	0.49	0.93	2.16	1.24	0.996
	10	0.06	0.54	0.92	3.72	1.41	0.979

Figure 5. Interfacial shear characteristic curves under different normal stress and wet-dry-freeze-thaw cycle conditions (T = −2℃)

图5. 不同法向应力和湿干冻融循环条件下界面剪切特性(T = −2℃)

Table 3. Parameter values of CPE model fit for curves with different normal stress and wet-dry-freeze-thaw cycle conditions

表3. 不同法向应力和湿干冻融循环条件下曲线的CPE模型拟合参数值

法向应力/kPa	湿干冻融循环次数/次	a	b	k	m	n	R^2
25	0	−0.13	1.20	0.36	3.18	1.23	0.997
	1	0.20	1.08	0.29	3.10	1.03	0.997
	3	0.36	0.99	0.31	1.64	1.12	0.954
	5	0.19	0.96	0.28	1.92	1.05	0.994
	10	0.15	0.96	0.31	1.83	0.91	0.987
50	0	−0.16	0.89	0.53	2.01	0.99	0.997
	1	0.00	0.36	1.05	4.61	1.08	0.983
	3	−0.42	1.28	0.50	2.93	1.62	0.999
	5	0.32	0.41	0.49	1.00	1.64	0.996
	10	0.28	0.83	0.49	1.72	1.22	0.993
75	0	−0.25	0.86	0.75	2.27	1.10	0.999
	1	0.35	0.62	0.80	2.01	1.19	0.999
	3	−0.30	0.90	0.69	3.17	1.69	0.998
	5	0.06	0.86	0.64	3.48	1.18	0.999
	10	0.22	0.53	0.70	1.64	1.46	0.997
100	0	−2.62	1.51	1.52	0.69	0.66	0.999
	1	0.08	0.75	1.64	3.45	1.12	0.979
	3	0.25	0.51	1.46	1.71	1.31	0.999
	5	−0.69	1.06	1.53	3.92	1.53	0.993
	10	0.12	0.63	1.42	2.74	1.16	0.998

3.4. 部分典型曲线形态及函数表达式

选用合适的本构模型可以准确的描述应力应变曲线的变化规律，能够详尽的拟合冻土复杂的曲线形态，对岩土工程相关的数值计算也有指导意义。从前面的拟合结果来看，CPE模型能够对软化型和硬化型中形态各异的曲线做到较好的拟合。基于不同水热循环工况试验数据，从中选取了部分有典型代表特征的曲线，结合前面拟合得到的适用于复杂水热循环作用下满足冻土–土工布界面剪切工况的模型参数，给出了拟合应力位移曲线对应的函数表达式，图2中不同应力位移曲线对应的函数表达式分别如式(6)~式(10)所示。

强软化脆性破坏型：

$f\left(x\right)=\left(0.03x^{5.45}-0.66\right)e^{-0.86x^{1.38}}+0.66$ (6)

弱软化脆性破坏型：

$f\left(x\right)=\left(0.27x^{0.92}-0.54\right)e^{-0.08x^{2.32}}+0.54$ (7)

塑性稳定型：

$f\left(x\right)=\left(-0.17x^{2.20}-0.65\right)e^{-0.95x^{1.26}}+0.65$ (8)

强硬化塑性破坏型：

$f\left(x\right)=\left(-0.18x^{2.23}-0.87\right)e^{-0.69x^{1.09}}+0.87$ (9)

弱硬化塑性破坏型：

$f\left(x\right)=\left(-0.16x^{2.01}-0.53\right)e^{-0.89x^{0.99}}+0.53$ (10)

4. 结论

冻土与结构接触面的本构模型是分析冻土区工程结构与冻土相互作用、评价工程安全稳定性的基础和关键。本文利用温控直剪仪在不同水热循环条件下进行了一系列冻土–复合土工布界面直剪试验，通过改进的复合幂–指数本构模型对部分试验数据进行了拟合验证，得到如下结论：

1) 针对不同水热循环条件下的软化型和硬化型曲线，改进的本构模型函数式均能实现较好的拟合，不同试验数据的拟合相关性R²均大于0.954，表明模型拟合效果比较准确。

2) 通过非线性拟合可以快速求得模型的参数值，模型参数值的变化有规律可循，且不同参数有其代表的物理意义，如a值可以反映曲线的类型，a > 0时曲线偏向塑性破坏型，a < 0时曲线偏向脆性破坏型。m值反映曲线上凸或下凹的程度。通过对湿干和湿干冻融条件下反映残余强度的k值进行分析，发现k值的变化幅度与残余强度实际值的变化幅度相吻合，且湿干冻融条件下k值的变化幅度大于湿干条件。

NOTES

^*通讯作者。