土与结构之间的相互作用一直以来都是岩土工程界关注的焦点之一，主要表现为基础、挡土墙、混凝土衬砌、土石坝防渗体等与岩土体之间的变形、传力及承载接触问题．接触面(带)两侧的介质存在刚度上的差异导致变形不同，在不同外荷载组合下，可能会产生法向张开或闭合等现象，使得土与结构之间的传力机制和结构受力响应变得更加复杂[1]．例如水工隧洞混凝土衬砌结构，混凝土衬砌自重较大，混凝土或灌浆浆液的冷凝、干缩等因素会导致围岩与结构之间可能产生局部初始缝隙．在内水压力作用下缝隙减小或闭合，在检修工况外水压力作用下缝隙又可能变大．只有将土与结构的接触条件与联合受力特性作为一个整体进行分析，才能更合理地揭示工程结构的传力机理．在计算过程中，有两个难题须解决：一是土与结构间接触边界约束或非连续模拟(几何模型)；二是接触或脱离力学行为模拟(本构模型)．目前主要有离散单元法、计算接触力学和接触面单元模拟法．离散单元法对介质运动学方面特性的描述比较直观、生动，但当须比较准确地确定接触面上的应力状态时，该法显得过于粗糙．现代计算接触力学方法在处理多体接触特性、非协调网格和不连续变形方面表现出很强的适应性[2-4]，但主要用于刚-刚接触分析，无法反映界面区域的剪缩和剪胀法向变形[5]．连续介质力学中的接触面单元以4节点Goodman单元和Desai薄层单元为代表，可以模拟接触面黏结、滑移、脱开等各种接触状态，且各大型商业数值模拟软件中都有所涉及．文献[6]采用接缝单元模拟衬砌与围岩之间的相互作用，研究了接缝大小、接触角变化对承载能力的影响；文献[7]开发了塑性-损伤-接触特征的开裂混凝土本构模型，研究衬砌结构在内水压与温度场联合作用下的开裂特征．文献[8-11]分别基于ANSYS，FLAC3D，ABAQUS和MIDAS中的接触单元研究隧道衬砌与围岩间的接触作用机理，分析衬砌结构的受力特性．可见：采用大型商业软件解决工程中的接触问题非常便捷，但共同的缺陷是处理两侧节点有初始空隙或法向相互嵌入时的接触精细算法考虑不周．本研究将在常规的连续介质力学模型基础上，引用FLAC3D中的Interface界面单元模拟水工隧洞混凝土衬砌与围岩的接触特性，对模拟初始缝隙功能及边界节点发生嵌入时的迭代计算过程进行改进．采用厚壁圆筒弹性和弹塑性解析解验证数值模拟结果的正确性及精度；分析水工隧洞在不同外荷载作用下衬砌与围岩有条件联合承载的变形和受力机制，为衬砌结构设计提供依据．1 接触面算法改进FLAC3D自带的Interface接触面模型是由m个3节点三角形平面构成的无厚度接触单元，可以设置在空间任意位置．每个时步循环时，通过计算得到目标面实体单元接触节点对之间的法向位移和剪切位移，再传递到Interface单元对应节点上，分别计算界面法向力和切向力[12]．Interface节点荷载取决于该节点所处的接触状态，在计算过程中，节点对可能存在黏结、滑移、张开和嵌入4种状态．1.1　初始缝隙的模拟及算法设界面的初始缝隙为δ0，当界面两侧节点对相对位移dun=|un1-un2|≤δ0(un1和un2分别为接触面两侧节点对的法向位移)时，界面处于张开状态，节点对间不传力，接触面两个方向上的接触应力为0 N，两侧节点在法向上可以自由伸缩，相当于不受约束的自由边界；当dun=|un1-un2|δ0时，界面两侧接触，节点对开始传递荷载，此时接触应力是由超过初始缝隙的那部分嵌入位移产生的．当dun≤δ0时，σn(t+Δt)=0， τs(t+Δt)=0；当dunδ0时，σn(t+Δt)=Kn(dun-δ0)，τs(t+Δt)=τsit+KsΔus(t+Δt/2)，式中：σn(t+Δt)和τs(t+Δt)分别为接触面某单元t+Δt时步法向应力和切向应力；τsit为t时步的初始切向应力；Δus(t+Δt/2)为t+Δt/2时步产生的相对剪切位移增量(Δt/2表示中心差分法计算)；Kn和Ks分别为接触面的法向刚度和切向刚度．1.2　接触边界节点嵌入修正方法界面可能存在的4种接触行为中，黏结、滑移和脱开状态在计算过程中可以预先得到并允许存在；而嵌入状态则表明接触体之间发生了相互侵入，这对于混凝土衬砌和岩体材料来说是不可能存在的．由于在当前计算步，不能预先判断物体受力变形后的状态，接触面两侧节点对的相对法向位移超过它们之间的初始距离产生而产生了嵌入．在计算过程中，应该在节点对上施加反向荷载，使相互侵入的位移退回去，最终将通过迭代计算使节点对处于黏结或滑移状态．设接触面两侧某节点对的法向相对嵌入量为dun，有dun=|un1-un2|，当dunδ0时，节点对发生嵌入，对节点对的位移进行修正．假定两个节点应分别调整法向位移量为dun1和dun2，且有：dun1+dun2=dun；dun1=k1dun；dun2=k2dun；k1+k2=1，式中k1和k2为两个节点位移应该调整的比例，取决于两接触体的相对刚度．修正后，节点对的法向位移相等，即处于法向黏结状态，un1N=un2N=k2un1+k1un2，式中：un1N和un2N为发生嵌入的节点对修正后法向位移；上标N表示修正后状态．节点对位移修正后，更新接触面两侧节点力，dFn1=Kndun1;  dFn2=Kndun2，式中dFn1和dFn2为节点对法向力增量．2 Interface模型验证根据文献[13]中的例子，考虑一个在均匀内压p作用下的厚壁圆形套筒，内筒外壁与外筒内壁之间的空隙间距为δ0=1.0 mm，内筒的内径a=1.0 m，外径c=1.3 m，外筒的外径b=2.0 m．厚壁圆筒的材料参数为：变形模量E=210 GPa，泊松比μ=0.3．Interface接触面单元法向和切向刚度均为560 GPa，内摩擦角和内聚力分别为40°和0.8 MPa．假设筒为无限长，可当平面应变问题处理，取1/4模型进行计算，如图1所示．10.13245/j.hust.220119.F001图1厚壁圆形套筒1/4模型2.1　弹性解验证2.1.1　解析解理论计算式厚壁圆筒仅受内水压力p作用时，其径向位移(u)及主应力(σr，σθ)计算公式[14]为u=1+μEpa2b2-a2(1-2μ)r+b2r;σr=a2pb2-a21-b2r2;σθ=a2pb2-a21+b2r2,式中：r为圆筒上任意一点至中心的距离；σr和σθ分别为最小和最大主应力．内筒外壁的径向位移为δ0时，所需径向压力为pδ=Eδ02(1-μ2)c2-a2ca2．在径向内压由0 Pa增至pδ的过程中(记为阶段1)，由于内外套筒间初始空隙的存在，只是内筒受力自由变形，不受外筒约束；而外筒处于无应力无变形状态，即内筒和外筒的应力分别为(上标i和e分别为内筒和外筒)：当ar≤c时，σr1i=a2pδc2-a21-c2r2,σθ1i=a2pδc2-a21+c2r2;当cr≤b时，σr1e=0,σθ1e=0.在径向内压由pδ 增至p的过程中(记为阶段2)，由于初始空隙在上一过程结束时已经闭合，内外筒达到联合承载条件，剩余荷载p-pδ由内外筒共同承担，相当于一个内径为a、外径为b的整体单筒，产生的内、外筒(a≤r≤b)应力可统一表示为：σr2=a2(p-pδ)b2-a21-b2r2;σθ2=a2(p-pδ)b2-a21+b2r2.2.1.2　数值解与解析解的比较初始缝隙δ0=1.0 mm弹性模型筒壁应力和位移数值解与解析解比较见图2．由图2可见：弹性数值解与解析解非常接近，初始缝隙的存在导致界面处的环向应力σθ不连续，接触法向应力分布均匀．初始缝隙的存在给内筒在初始受力阶段提供了临空自由面，不受外筒的约束，因此变形较大，而外筒处于无受力状态还来不及变形，只有在缝隙闭合后才联合受力共同变形，因此最终荷载加载完毕时界面处(r=1.3 m)两侧对应节点径向位移不连续．10.13245/j.hust.220119.F002图2弹性模型筒壁应力和位移数值解与解析解比较2.2　弹塑性解验证对于有初始缝隙的情况，塑性解比较复杂．此处不妨假定初始空隙δ0=0 mm，接触面两侧节点也有可能相互嵌入，不影响改进算法的验证.2.2.1　解析解假定材料服从Tresca屈服准则，屈服函数[14]为f(σr,σθ)=σθ-σr-σs，式中σs为屈服极限，取277 MPa．当f(σr,σθ)0时，单元处于塑性应力状态，塑性区半径ρ和内压p的关系为pσs=lnρa+121-ρ2b2．塑性区(a≤r≤ρ)内的应力分量、位移的表达式为：σr=-p+σslnra;σθ=-p+σs1+lnra;    u=(1-2μ)(1+μ)Eσsr2ρ2b2+rlnrρ-r2+(1-μ2)σsρ2Er.弹性区(ρr≤b)内的应力和位移计算式为：σr=σsρ221b2-1r2;σθ=σsρ221b2+1r2;u=1+μ2Eσsρ2b2(1-2μ)r+b2r.2.2.2　数值解与解析解的比较弹塑性模型筒壁应力和位移数值解与解析解对比如图3所示．应力和位移连续分布表明：接触面一直处于接触状态，接触法向应力约为106 MPa．数值解与解析解非常接近，塑性区内应力误差不超过3%，靠近外边界稍大但不大于10%，位移最大误差为1.5%.10.13245/j.hust.220119.F003图3弹塑性模型筒壁应力和位移数值解与解析解对比通过弹性模型及弹塑性模型数值解与解析解的对比，验证了Interface接触面单元可以模拟具有初始缝隙的接触问题，同时验证了接触后两侧节点对发生嵌入时位移及节点力修正方法的正确性，模拟精度较高．3 工程实例3.1　参数设定某水电站泄洪放空洞有压段通过的地层为孔王溪组灰岩∈2k1-3、泥质白云质灰岩～龙王庙组∈1l2-1钙质粉砂岩、板岩；∈2k1-2和∈1l2-2中岩溶较发育，以Ⅲ类围岩为主．该段隧洞为圆形，开挖断面直径为10.5 m，衬砌厚度为1.0 m，对溢0~010.500断面衬砌结构进行分析，内水水头为107.5 m，外水水头考虑折减后为37.5 m．岩体及主要结构面的物理力学参数见表1，Interface界面法向和切向刚度均为50 GPa，黏聚力为0.8 MPa，内摩擦系数为0.84．混凝土采用各向同性弹性损伤非线性模型计算[15]，假定衬砌与围岩在施工结束后是紧密接触的，分别考虑在自重、运行期内水压力、检修期外水压力工况下衬砌与围岩的接触运动状态、结构受力、破坏情况及安全系数[16]分布．安全系数β=Fu/F，式中：Fu为按材料极限强度计算的截面轴向或偏心受压(或受拉)构件的极限承载力；F为截面实际轴力．10.13245/j.hust.220119.T001表1岩体及主要结构面的物理力学参数介质E/GPaμ黏聚力 /MPa内摩擦系数抗拉强度/MPa粉砂岩、板岩(∈1l2-1)80.2501.01.120.70中～厚层灰岩(∈1l2-2)100.2401.21.120.90泥灰岩、白云质灰岩(∈2g/∈1l2-2)60.2600.80.840.70灰岩夹白云质灰岩(∈2k1-1/∈2k1-3)80.2400.80.930.30白云质灰岩(∈2k1-2)40.2600.40.650.10混凝土240.1671.753.2　施工期自重工况衬砌在自重作用下的位移见图4(a)，整个结构向下变形，总位移量为0.43~0.91 mm，上部大下部小，顶部与围岩产生张开脱离，最大脱离距离为1.29 mm(图4(c))，脱离范围约为顶拱圆心角120°，与布加耶娃模型假定围岩抗力分布的范围一致[17]．衬砌在自重作用下没有拉裂区，围岩有少量剪切塑性破坏，见图4(b)．最大和最小主应力较小(图4(d))，顶部和底部内侧受拉，两侧外部受拉，整个断面的安全系数较大(图4(e))，说明自重工况下结构很稳定.10.13245/j.hust.220119.F004图4自重工况计算结果3.3　运行期内水压力工况衬砌在内水压力作用下的位移见图5(a)，整个圆环结构有向外扩张的趋势，衬砌右下侧变形最大为0.94 mm，在较高内水压力作用下初始缝隙闭合，接触面处于黏结状态(图5(c))，衬砌与围岩联合受力．衬砌在内水压力作用下产生环向拉应力1.79~4.75 MPa(图5(d))，衬砌基本都被拉裂，围岩也分布着少量剪切塑性破坏区，见图5(b)．整个断面的安全系数都小于1，说明内水压力工况下结构不安全，应按照透水衬砌原理对结构进行配筋.10.13245/j.hust.220119.F005图5内水工况结果3.4　检修期外水压力工况衬砌在外水压力作用下的位移分布见图6(a)，整个结构向内收缩，衬砌最大变形发生在上部，最大值为0.99 mm，衬砌顶拱和围岩之间接触面脱离，最大张开缝隙为1.38 mm(图6(c))．此工况下几乎没有拉裂区，围岩有少量剪切塑性破坏，见图6(b)．最大和最小主应力均为压应力，值较小且分布较均匀(图6(d))，安全系数均大于1(图6(e))，结构处于安全状态.10.13245/j.hust.220119.F006图6外水工况结果综上，当隧洞开挖完成浇筑混凝土衬砌时，围岩应力释放及变形基本完成，衬砌和围岩只承受自重作用，它们之间仅存在简单的接触关系．由于施工方法、混凝土冷却收缩及衬砌自重等原因，围岩和衬砌不完全紧密接触，通常在隧洞顶部存在微小的缝隙．对于内水压力工况，初始缝隙闭合前内水压力完全由衬砌承担，衬砌单元有较大的变形空间，在与围岩接触之前其外表面相当于临空面，在较小的荷载作用下变形较大并很快能与围岩接触，导致初始缝隙闭合．闭合后围岩与衬砌接触，二者联合受力共同承担内水荷载，产生共同变形．外水压力工况类似于自重工况，对衬砌结构运行影响不大．通过以上分析可以得出：内水压力工况为控制工况；在内水压力作用下，衬砌开裂，结构处于欠安全状态，应按透水衬砌原理对结构进行配筋．4 结论a．引用FLAC3D中的Interface单元模拟衬砌与围岩之间的接触作用，推导了接触面两边存在初始缝隙和接触边界发生嵌入时的修正迭代算法，与厚壁圆筒弹性模型和弹塑性模型的解析解进行对比，精度较高.b．改进的Interface模型反映了衬砌结构在不同外荷载作用下与围岩有条件联合承载的各种受力机制．在自重和外水压力工况下，衬砌结构顶部与围岩存在脱离情况，结构安全性较高；在内水压力工况下，衬砌与围岩联合承载，混凝土开裂，应按透水衬砌原理对结构进行配筋.10.13245/j.hust.220119.F007