摘要：本文回顾了常用的几种岩石粘弹性本构模型，依据蠕变柔量的概念，用位移和蠕变柔量两步反分析法从粘弹性本构模型的一般表达式中辨识出岩石的本构模型.先由相应的位移实测值用解析法反演岩石的蠕变柔量，再由岩石的蠕变柔量运用非线性优化技术辨识出岩石本构模型表达式，最后进一步由模型参数反算出岩体的粘弹性参数.文中最后给出了工程实例.
关键词：粘弹性本构模型 蠕变柔量 两步反分析法 反问题 岩石
岩体属于各向异性流变介质，由于其客观复杂性，在岩体理论分析和数值模拟方面，参数和模型的正确给定是岩石力学研究中的两大难题，它们的合理性将大大增强岩体工程分析设计的可靠性.新奥法及目前广泛应用的现场监控法或信息反馈法是把岩石的变形观测作为对原设计及方案进行修正的依据[1，2].这些方法紧密结合工程实际指导修改原设计，收到了良好的效果，已广泛应用于岩土工程实践领域.参数取值，传统的方法大多凭经验、工程类比来进行，常常不能作出科学的判断.参数反演作为参数辨识的一种方法，是基于实测位移反求系统某种参数的一种逆问题，对确定认识系统和进一步正演分析相当重要.参数的辨识，包括模型参数的辨识，是在模型的式给定的情况下进行的，模型的给定与实际相符合的程度显然相当重要.工程岩体是复杂的不确定系统，岩体的流变力学模型较多，模型识别的问题有着重要的理论意义和实用价值.在岩体本构模型的辨识，尤其是与时间有关的粘弹性本构关系辨识方面当前还有许多问题值得研究.
本文回顾了工程应用和科学研究中几种常用的岩石粘弹性本构模型，然后依据蠕变柔量概念，用位移和蠕变柔量两步反分析法从粘弹性本构模型的一般表达式中辨识岩石本构模型，第一步由相应的位移实测值用解析法反演岩石的蠕变柔量，第二步由岩石的蠕变柔量运用非线性优化技术辨识出岩石本构模型的表达式，最后进一步由模型参数与粘弹性参数的关系反算出岩石的粘弹性参数.
1 粘弹性本构模型回顾
工程岩体的长期稳定性是当今岩土工程领域的一个十分重大的前沿问题[3].工程岩体的变形破坏是有时间过程的，具有时间效应，岩体流变是岩体的重要基本力学特性和行为，对于某些实际岩体工程，研究变形的时间效应有着极其重要的意义[8,9].真实岩体是一种非连续、非均质、各向异性的流变介质，目前很难用数学手段作出适当精确描述.在一定情况下，可将岩体视为似连续、均匀化、类各向同性介质，采用经验与理论相结合，定量分析定性使用的原则，对工程设计与起重要的定性指导作用[10,11].在这一前提下，在岩土工程应用和科学研究中提出的粘弹性本构模型主要包括经验模型和组合模型[5](弹性和粘性元件的不同组合而成).表1列出了常用的几种微分型组合模型的公式及其特征(H代表弹性元件，N代表粘性元件).在一般状态下，这些组合模型的本构方程微分形式的一般式可表示为[5]
表1 常用岩石粘弹性微分型组合本构模型
式中：
式(1)即线粘弹性微分型本构关系的一般表达式.
2 粘弹性模型识别
2.1岩体介质蠕变柔量的确定在弹性岩体中开挖任意形状断面的隧道，应用平面复变影射、保角变换方法可得围岩内任一点在曲线坐标中的位移为[7]
由弹性理论确定.其中：ω(ζ)为垂直隧道轴线的物理平面Z上的非圆形洞室外域到数学平面ζ内单位圆外域的映射函数；φ(ζ)和ψ(ζ)为满足应力边界条件的复势函数.
应用弹性—粘弹性对应原理，对照式(2)，得出任意断面的隧洞内任一点因开挖引起的粘弹性位移在曲线坐标中的表达式为[8]
D(x)为荷载逐步释放系数，若计算断面远离开挖面，则可作为平面应变问题处理；若接近开挖面，则严格说来是空间问题.为使问题简化，考虑了开挖面空间效应后，仍作为平面问题来处理.开挖面空间效应通过洞周释放荷载逐步释放代替瞬间完全释放来体现.据研究资料[4]表明，其取如下形式：D(x)=1-0.7exp(-3.15x/2a).其中，a为隧洞半径；p0为垂直向初始地应力；G为剪切模量；r为所考虑点到隧道轴线的距离；λ为侧压力系数；k为体积模量，D(x)为应力释放系数，x为计算断面距开挖面距离.
在流变岩体中开挖任意形状的隧道，在开挖之前布置量测仪器，可测到由于开挖引起的全部位移，在ti时刻可测到n个测点由于开挖引起的相对径向全部位移，记为Uk(ti),k=1,2,，…，n,i=1,2,,…，l(假定量测了l个时段)，代入式(4)，则有蠕变柔量Jl(ti)和广义蠕变柔量J2(ti)的线性方程组
2.2粘弹性应力应变关系的确定线粘弹性微分型本构关系的一般表达式如式(1)所示，对其进行拉氏变换，考虑光滑化假定，得[7]
对于
(s)=0有重根或虚根的情况，J(t)中含有时间t的阶次项和正余弦项，且未知参量个数增加，给反演优化大大增加了难度，收敛性差，不利于程序的统一编制，这里作了一些简化，在工程实际中，复杂的高阶次模型也不必要.
式(10)中，在模型阶次确定时，参数Ai(i=1～2m 1)与模型式(1)中的模型参数p0,p1,p2,…,pm,q0,q1,q2,…,qn有一定的对应关系，由式(9)和式(10)可推导得出，模型参数均可用Ai简单运算表达.
表2 常用岩石粘弹性模型参数表
利用实测位移由式(6)可求出众多不同时段ti的蠕变柔量Ji(ti)的最佳估计值，又由式(10)可得相对应的J(ti)(i=1～L)，因而构造如下的非线性优化目标函数：
（11）
式(10)是Ai的非线性函数，该问题是一非线性最小二乘优化问题，带约束的隐式非线性优化问题，这里采用了可变容差优化方法，它是在单纯形法和复合形法的基础上变化而来，具体见文献[6].
上述模型参数估计是在模型已经确定的情况下进行的，即偏微分方程的阶次m,n为已知的.这里模型的最佳阶次m,n的确定不是通过数学推导求出，而是通过试验来确定的.令模型阶次m,n的取值从1开始，分别以步长1递增，比较m,n取不同值时最优估计式(11)的值，即拟合误差.通常随着阶次增大，拟合误差下降.但当阶次为最佳阶次时，拟合误差达到最小；然后随着m,n的增大，拟合误差趋向增大.
3 粘弹性参数的反算
粘弹性力学本构模型中的模型参数为粘弹性参数的函数，若辨识出岩石的粘弹性微分型模型，则可反算出岩石的粘弹性参数，这里粘弹性参数包括弹性模量和粘性系数.在研究和工程应用中，常用的几种粘弹性模型已在表1中列出，它们的模型参数与粘弹性参数的关系如表2，可以反算出相应的岩体力学参数.
4 工程实例
某地下工程开挖一条半径为2.0m的试验洞，围岩体可视为均匀各向同性粘弹性体.围岩垂直初始应力为4.5MPa，侧压力系数为1.5，泊松比为0.26.实测数据由通过圆中心的水平和竖直两条测线L1和L2得到.以某断面被开挖瞬间作为时间起点进行观测，实测数据如表3所示.
表3 各测线的位移实测值(单位：mm)
将上述实测结果及洞室几何尺寸等数据输入到用本文方法编制的程序中，进行迭代计算.第一步利用位移量测值由式(6)得出不同时段的蠕变柔量值，第二步由计算出的不同时段的蠕变柔量值通过式(11)优化计算，得到如表4所示结果.由表4可得出，模型的阶次为1阶，相应辨识出的模型为0.329×10-8σ 0.366×10-4
=0.843×10-3ε 
由实测数据初始时刻位移及辨识的模型形式初步确定，该岩体粘弹性模型可作为三参量Kelvin-Voigt模型，由表2反算粘弹性参数为E1=2.73E 4，E2=4.14E 5，η=4.91E 8，弹模型单位为MPa，粘性系数单位为MPa·d.
上述结果同已有的勘测和试验成果较一致，室内试验的结果为E1=2.51E 4，E2=4.62E 5，η=5.36E 8.
如果量测数据较多，可多算出几组结果，统计出其中拟合误差最小的阶次所占的百分比，通过百分比最高的阶次来得出模型的阶次及模型参数.
表4 模型辨识计算结果
5 结语
岩体本构模型辨识及参数反演是岩石力学理论和工程实践的重要问题，本文依据蠕变柔量明确的物理意义，初步探讨了通过两步反分析辨识出岩体的粘弹性模型，并对于常用的粘弹性模型反算出其粘弹性参数，粘弹性模型的表达式取自于线粘弹性模型的一般形式.算例表明有较好的工程实用价值.对于工程岩体，完全符合某一种力学模型是难以做到的，但可以根据一定的匹配原则，在允许误差内选取最恰当的力学模型.岩体本构模型辨识当前研究相对较少，如何选取适当的匹配准则来获得最佳模型仍是今后需要进一步深入研究的课题.
减小字体
增大字体