在存在断层构造的富水地层，隧道施工常伴随着极大的突涌水风险，有必要开展隧道渗流场研究及涌水量预测．文献[1]通过保角坐标变换的方式将无限地层与隧道的关系转化为两个同心圆，优化边界条件，使渗流微分方程容易求解；文献[2]得到了半无限地层中渗流微分方程的通解；文献[3]总结了目前各国学者所提的隧道涌水量解析公式或经验公式，并通过对比计算分析了各公式的适用性．除了常规的半无限地层，针对在复杂构造地层条件下的隧道涌水量计算，学者们也做了许多尝试．文献[4]从地质条件出发，统计并总结了断层构造区域的地质特征，并将断层带分为三个区域：位于断层中心的断层破碎区Ⅰ；在中心区范围之外但存在节理裂隙的断层影响区Ⅱ；受断层影响较小的普通围岩区Ⅲ．文献[5]在点涌水模型的基础上，利用势函数叠加的方法得到隧道揭露断层破碎区Ⅰ时的掌子面涌水量公式；文献[6]针对巴塞罗那地铁9号线穿越断层破碎带面临的突涌水风险，通过三维地质特征调查，建立了包括各类地层特征和断层在内的精细化数值模型；文献[7]考虑到断层裂隙的影响，认为该区域附近渗流并非稳定的达西渗流，从而对渗透系数进行非达西修正，并通过有限元迭代的方法对深埋隧洞涌水进行预测；文献[8]同样利用有限元方法，对矿山隧道遭遇导水管道型岩溶时发生的突涌水演化特征进行了建模分析；文献[9-10]通过保角变换推导了双线平行隧道和上下分布的双洞隧道的地下渗流场，虽然没有特殊构造地质，但其叠加的思想仍然值得借鉴．对于普通地层下隧道渗流问题，许多成熟的模型将其视为半无限空间的地下渗流．在具有特殊地质条件下，若存在断层影响的渗流，介于复杂的边界条件，则学者们多选择有限元的方法来解答，而从理论解析的角度突破较少．本研究针对隧道在富水地层断层影响区Ⅱ可能面临的突涌水问题，利用断层等效和镜像叠加的方法来简化模型，并通过保角坐标变换进行求解，从而得到断层影响区Ⅱ内隧道地下水头重分布，提出了水头值计算和涌水量预测的半解析公式．1 等效断层的确定当隧道附近存在地层断裂面时，横断面内隧道与断层的关系如图1所示．图中：Ⅰ区为地层区域；Ⅱ区为隧道衬砌区域；h为隧道圆心到地下水位的距离；ω为横断面内断层的倾角；a为水位面到断层面的交点到隧道圆心的水平距离；d为横断面内隧道圆心到断层面的水平距离；r为隧道半径；R为衬砌圈的外半径．10.13245/j.hust.210115.F001图1横断面内隧道与断层关系图断层对隧道附近渗流场的影响体现在两方面．其一是断层影响范围内的节理裂隙，加速了渗流过程，该影响可以对岩土体的渗透系数KS进行非达西修正来表征[7]，KS=J-1(-μ/(2kρwβ)+(-μ/(2kρwβ))2+gJ/β), (1)式中：J为水力梯度；μ为液体的动力黏滞系数；k为岩土体的渗透率；ρw为液体的密度；β为非达西流影响系数[7]，取β=3×10-4k-0.613 4；g为重力加速度．其二是富水的断层面会起到补水边界的作用．此处将断层面作为水头值恒定的边界，水头值与地下水位线的水头值相等．在解析过程中，为了简化计算，将横断面内具有一定倾角的倾斜断层等效为倾角为90°的竖直断层．竖直断层与隧道的水平距离L决定了等效效果．将隧道涌水量Q作为衡量等效效果的指标，若等效工况的涌水量与实际工况相差不大，则等效效果较好．在图1所示的倾斜断层工况中，当h不变时，隧道与断层的相互位置关系由a和ω决定．当ω=90°时，断层竖直，此时L=a；当ω=0°时，断层水平，此时L趋向于正无穷，地层相当于半无限渗流空间．将图1中圆心到断层面的水平距离d作为等效距离，即L=a+h/tan ω．(2)式(2)保证了等效在物理上的真实性．为检验等效效果，利用Abaqus分别进行倾斜断层工况和等效断层工况的模拟计算．隧道离左边界的距离为200 m，而右边界主要由断层控制，隧道离模型底部边界185 m，埋深为65 m，各工况中r=6 m，R=6.5 m，h=55 m，a=50 m．ω=0°，10°，20°，30°，40°，50°，60°，70°，80°，90°时，L= ∞，362，201，145，116，96，82，70，60，50 m．以隧道轴线高程的水平面为基准面，在水位面和断层面上，设置总水头值为恒定值h；代表地层无限远处的模型水位面以下左右边界和底部边界水头值也恒定为h，而隧道开挖后其轮廓线孔隙水压为0 Pa．此外，为了重点研究地下水流情况，对所有土体施加了固定约束．假定地下模型为各向同性稳定渗流问题，同时参考文献[7]对土体渗透系数KS进行非达西修正，首先将ρwgk/μ作为初始渗透系数K0，带入模型运算，得到初始地下渗流场；然后根据式(1)中KS与J的关系，将隧道附近的初始水头梯度J带入式(1)进行修正．选取初始渗流场中隧道周围的平均水头梯度作为J的取值，式(1)中各参数大小及修正计算得到的渗透系数：k=5×10-9 m2，μ=1 Pa⋅s，β=37.07 m-1，KS=0.18 m/h，衬砌渗透系数KL取0.18 mm/h．对上述模型各工况进行数值求解，实际倾斜断层及等效竖向断层涌水量对比如图2所示．10.13245/j.hust.210115.F002图2实际倾斜断层及等效竖向断层工况涌水量对比可以看出：等效工况所计算的隧道涌水量与实际工况的涌水量十分接近，且变化趋势一致，差距最大的工况涌水量分别为0.788 9和0.788 4 m3/h，各点差距均在0.1%之内，说明将d作为L可行．2 渗流场解析推导及涌水量计算半无限地下空间的隧道渗流模型地下总水头ϕ满足渗流微分方程∂2ϕ/∂x2+∂2ϕ/∂y2=0．(3)文献[2]中利用保角变换得到半无限地层中式(3)的通解，ϕ=C0+C1ln ρ+C2(ρn-ρ-n)cos(nθ)，(4)式中：C0，C1和C2为常数；ρ=ε2+η2；θ=arctan(η/ε)．而ε，η与x，y关系为：ε=x2+y2-(h2-R2)/x2+(y-h2-R2)2;η=2h2-R2x/x2+(y-h2-R2)2. (5)为了后续表述方便，此处将ρ，θ与x，y的关系分别表示为ρ=f(x,y)，θ=g(x,y)．当地层中有断层存在时，为了考虑断层补给边界的影响，参考群井降水理论，采用镜像叠加的方法，将有限的地下渗流空间等效为半无限空间．断层地层中隧道渗流模型如图3所示．10.13245/j.hust.210115.F003图3断层地层中隧道渗流模型图3中，实际地层为Ⅰ区，衬砌区域为Ⅱ区．Ⅰ区采用正交坐标系统x-y，坐标原点位于隧道轴线上方的水平面上；而Ⅱ区采用极坐标系统ρ3-θ3，坐标原点位于隧道圆心．在Ⅰ区中，实际排水隧道以等效断层面作为镜像面，映射出一个虚拟的供水隧道．根据叠加原理，Ⅰ区的地下渗流场为各个排(供)水源影响下的子渗流场的线性叠加，所以依据式(4)，图3所示的模型中水头值分布为    ϕΙ=C0+C11ln ρ1+C12(ρ1n-ρ1-n)cos(nθ1)+C21ln ρ2+C22(ρ2n-ρ2-n)cos(nθ2), (6)式中：C11，C12，C21和C21均为常数；ρ1=ε2+η2=f(x,y)，(7)θ1=arctan(η/ε)=g(x,y)，(8)ρ2=(ε')2+(η')2=f(x',y')，(9)θ2=arctan(η'/ε')=g(x',y')，(10)其中，ε-η为Ⅰ区坐标系统x-y相应的变换坐标系统，x'-y'为虚拟地层的坐标系统，ε'-η'为相应的变换坐标系统，x-y与x'-y'坐标原点位置如图3所示．依据其位置关系可知：地层中某点在两个坐标系统中的坐标关系为x'=x-2L，y'=y．另一方面，由于Ⅱ区边界简单，在极坐标系统ρ3-θ3下，采用直接求解，可以得到满足渗流连续性方程的通解    ϕΙΙ=C3+C4lnρ3+(C5ρ3n+C6ρ3-n)∙cos(nθ3)+(C7ρ3n+C8ρ3-n)sin(nθ3), (11)式中C3~C8为常数．要求解式(6)和(11)中的常数，还须考虑模型中的边界条件．以水位线为基准面，首先考虑Ⅰ区．边界条件1 y=y'=0时，ϕ=0，此时由y=y'=0、式(7)和(9)可知，ρ1=ρ2=1，从而可以得出C0=0．(12)边界条件2 等效断层面水头值恒定为0，即当x=L，x'=-L时，ϕ=0，由x=-x'和式(7)~(10)可知，此时ρ1=ρ2，θ1=-θ2，从而可以得出C11=-C21,C12=-C22．(13)此时，式(6)变为    ϕΙ=C1lnρ1ρ2+C2(ρ1n-ρ1-n)cos(nθ1)-(ρ2n-ρ2-n)cos(nθ2). (14)边界条件3 考虑Ⅱ区，隧道内边界为自由出水界面，即ρ3=r时，其界面水头值等于其位置水头ϕΙΙρ3=r=-h+rsin θ3．(15)将式(11)带入式(15)，可以得到C3+C4lnr=-h;C5=C6=0;C7r+C8r-1=r, (16)此时，式(11)可以写成    ϕΙΙ=-h+C4ln(ρ3/r)+[C7ρ3+r2(1-C7)/ρ3]sinθ3. (17)要求式(14)和(17)的未知常数C1，C2，C4和C7，还须结合Ⅰ区和Ⅱ区交界面的情况，两侧区域在此边界上满足水头值相等，总流量相等，流速相等．边界条件4 当x2+(y+h)2=R2，ρ3=R时，ϕΙ=ϕΙΙ．(18)将式(14)带入式(18)．在Ⅰ区，衬砌外圈x2+(y+h)2=R2在ε-η坐标系上映射为圆，该圆的方程为ρ1=α=h/R-h2/R2-1，(19)而在ε'-η'坐标系上，则变为(x'+2L)2+(y'+h)2=R2的映射，由式(5)推导可知，此隧道面同样映射为圆，方程为    ε'+(A2+R2-4L2-h2)/(A2+2Ah+4L2+h2-R2)2+η'+4AL/(A2+2Ah+4L2+h2-R2)2=2AR/(A2+2Ah+4L2+h2-R2)2, (20)式中A=h(1-α2)/(1+α2)．因此在Ⅰ区衬砌外圈上，ρ1=α，ρ2则为与θ2有关的周期函数，此外衬砌外圈在ε'-η'坐标系中的映射圆半径很小，且R越小，ρ2的周期波动也越小，此处可将ρ2取值为一个波动周期内的平均值，即ε'-η'坐标系中原点到映射圆圆心的距离ρ2a，ρ2=ρ2a=(A2+R2-4L2-h2)2+16A2L2A2+2Ah+4L2+h2-R2 ．(21)同理，θ2也取一个波动周期内的平均值θ2a，根据式(20)，得θ2=θ2a=arctan4AL/(A2+R2-4L2-h2) ．(22)所以式(18)等号左边可以写成    ϕΙρ1=α=C1lnαρ2a+C2(αn-α-n)cos(nθ1)-C2(ρ2an-ρ2a-n)cos(nθ2a). (23)再将式(17)带入式(18)，同时，令C9=C7+r2(1-C7)/R2．(24)考虑ρ3-θ3与x-y之间的关系，式(18)等号右边可以写成    ϕΙΙρ3=R=-h+C4ln(R/r)+C9Rsin θ3=-h+C4ln(R/r)+C9(Rsinθ3-h)+C9h=(C9-1)h+C4ln(R/r)+C9y. (25)此外，根据保角坐标变换，在衬砌外圈上，y=A(α2-1)/(α2+1-2αcos θ1)，将其化成级数形式[2]则为y=-A-2A∑n=1∞αncos(nθ1)，(26)所以式(25)还可写成为    ϕΙΙρ3=R=(C9-1)h+C4ln(R/r)-C9A-2C9A∑n=1∞αncos(nθ1). (27)联立式(18)、(23)和(27)，可得到：    C2=-2C9Aαn/(αn-α-n);    C1=[C4ln(R/r)+(C9-1)h-C9A+∑n=1∞C2(ρ2an-ρ2a-n)cos(nθ2a)]/ln(α/ρ2a). (28)边界条件5 Ⅰ区和Ⅱ区在衬砌外圈边界上须满足总流量相等，即∫02πKSdϕΙdρ1ρ1=ααdθ1=∫02πKLdϕΙΙdρ3ρ3=RRdθ3．(29)将式(14)和(17)带入式(29)，式(14)中的ρ2和θ2分别由ρ2a和θ2a代替，则可以得到KSC1=KLC4．(30)边界条件6 同样，Ⅰ区和Ⅱ区在衬砌外圈边界上需满足流速相等，即KSdϕΙdρ1ρ1=α=KLdϕΙΙdρ3ρ3=R．(31)将式(14)和(17)带入式(31)，同时取一个特殊点辅助求解，即如图3中所示的点P，位于衬砌外圈边界上θ3=270°的位置，根据式(8)，此点在Ⅰ区ε-η坐标系统下，θ1=0°，此时式(31)可以写成    KSC1/α+KSC2(nαn-1+nα-n-1)=KLC4/R-KLC7+KL(1-C7)r2/R2. (32)联立式(24)、(28)、(30)和(32)，便可以求得各未知常数C1，C2，C4和C7，从而地层区域渗流场的水头解式(14)及衬砌区域渗流场的水头解式(17)便均可求得．为了得到涌水量公式，先考虑渗透系数的非达西修正．首先，将地层初始渗透系数KS取为ρgk/μ，然后利用式(14)、(24)、(28)、(30)和(32)得到初始地层渗流场分布，再计算衬砌外圈圆周的平均水力梯度J，J=∫02π∂ϕΙ∂ρ1ρ1=ααdθ1/(2πα), (33)将J带入式(1)中计算渗透系数KS，之后重复上述操作迭代求解，直至KS稳定．得到最终的水头解析式后，可得隧道涌水量Q=∫02πKL∂ϕΙΙ∂ρ3ρ3=rrdθ3=2πKLC4．(34)3 解析解退化验算及数值验算采用退化验算来检验所推导水头分布解析结果的正确性．在图1中，当a→∞或ω→0时，模型则变为半无限空间渗流模型，此时断层的等效距离L→∞．在地层水头分布式(14)中，当L→∞，根据式(9)和(21)，ρ2→1，ρ2a→1．同时，当不考虑衬砌作用，即R=r时，有：C9=1;C1=-A/ln α;C2=-2Aαn/(αn-α-n), (35)而地层的水头分布式(14)则变为    ϕΙ=-Aln ρ1/ln α-2Aα2nα2n-1(ρ1n-ρ1-n)cos(nθ1). (36)式(36)与文献[2]中关于隧道等水压边界条件下半无限空间水头分布式(13)一致，因此推导的隧道在断层影响区内的渗流解析式正确．利用Abaqus结合数值模拟，对隧道与断层的关系在几种不同的工况下进行了涌水量数值验算对比．分别对h=15，55 m；a=10，50 m；ω=10°，45°，80°，共12种工况进行了对比计算．各工况隧道涌水量的数值模拟计算结果和由式(34)所得的半解析计算结果如表2所示．10.13245/j.hust.210115.T001表2隧道涌水量的数值模拟计算与半解析计算结果h/ma/mω/(°)L/mQ/(m3·h-1)误差/%数值模拟半解析解151010950.217 6160.208 4304.4145250.217 9790.208 7004.4580130.2188520.209 1644.6350101350.217 5560.208 4194.3845650.217 6110.208 4704.3880530.217 6460.208 4894.395510103220.786 6000.751 2214.7145650.788 8150.753 3324.7180200.796 3740.760 5244.7150103620.786 4000.751 1974.69451050.786 8970.752 0924.6380600.788 4370.753 6484.62由表2可知：式(34)所计算的结果与数值模拟计算结果误差在5%以内．4 涌水量案例对比选取鸿图嶂隧道F0-6断层影响区域段落为分析对象，进行涌水量解析解与实际监测值对比．该区段网状裂隙发育，地下水主要为孔隙潜水和裂隙潜水两种．断层走向与隧道轴向夹角β=75°，倾角为80°，地下水至隧道的平均高差65 m，地层渗透率k=1.0×10-8 m2，液体的动力黏滞系数μ=1 Pa∙s，隧道半径r=5.9 m，衬砌厚度为0.5 m，衬砌渗透系数KL=0.18 mm/h．隧道与F0-6断层的平面位置示意图如图4所示．10.13245/j.hust.210115.F004图4F0-6断层与隧道平面示意图选择距离F0-6断层前50 m隧道范围内作为研究对象．在图4中，z为隧道某一断面距离断层面的距离，根据几何关系，断层等效距离则为L=ztanβ．将上述各参数带入式(14)、(24)、(28)、(30)、(32)及(34)中，得到鸿图嶂隧道在断层F0-6前50 m范围内的涌水量及现场实测涌水量分布如图5所示．10.13245/j.hust.210115.F005图5断层F0-6前50 m范围内涌水量及现场实测涌水量由图5可以看出：式(34)反映的涌水量变化规律与实测涌水量变化趋势一致，在靠近断层2倍洞径距离，涌水量急剧增加，在未施做防水措施前，极易发生突涌水害；随着断面距离断层越远，涌水量逐渐减小，在2倍洞径外变化基本稳定．本文方法计算得到的涌水量与实测值较为接近，最大相对误差为8.4%，说明了利用本文方法进行断层影响区内隧道涌水量预测具有较好的准确性．由式(28)、(30)、(32)和(34)可知：影响隧道涌水量的因素主要有隧道与断层的位置关系(由h和L表征)及地层或衬砌的水力学参数(由地层与衬砌的相对渗透系数KS/KL表征)．选取h=15，35，55，75 m；L=10，30，50，70 m及KS/KL=1 000，100，10，1，共64种工况进行涌水量计算．隧道涌水量在不同参数条件下的变化情况如图6所示．10.13245/j.hust.210115.F006图6隧道涌水量在不同参数条件下的变化情况1—KS/KL=1；2—KS/KL=10；3—KS/KL=100；4—KS/KL=1 000．由图6可以看出：相对渗透系数KS/KL较几何参数h和L对涌水量的影响更为显著，当KS/KL较小时，h和L的影响程度会被放大，说明衬砌可以很好地减少隧道涌水量；当KS/KL不变时，h越大，Q随L的变化斜率也会变大，说明h较大时，L的影响程度也会放大，所以在富水地层中进行隧道线路选择时，应尽可能考虑埋深设计，避免高水头压力．5 结论a. 建立了断层附近的隧道渗流简化模型，将隧道横断面内具有一定倾角的倾斜断层等效为倾角为90°的竖直断层，且其等效距离L为横断面内隧道圆心到断层面的水平距离d．b. 推导了富水断层影响区内，考虑衬砌作用的地下水头半解析公式，并得到了隧道涌水量计算公式，大大提高了计算效率．c. 经过多工况对比计算，发现在富水断层影响区内，地层与衬砌的相对渗透系数KS/KL，隧道与水位线的距离h，以及隧道与断层的距L均会影响涌水量大小，但相对来说，KS/KL较h和L对涌水量的影响更为显著．