考虑截面扁平化和翘曲的土-隧相互作用模型

摘要

基于广义梁理论与三维空间土抗力，将隧道视为在弹性地基梁上的具有弯曲、剪切、翘曲和扁平化4个变形模态的薄壁圆环结构，推导了横纵向相互耦合的土隧相互作用模型。结合案例验证了模型的有效性，并通过参数化分析得到了隧道厚径比、土抗力系数对隧道横纵向结构响应的影响规律。研究表明：高应力压载条件下大直径隧道翘曲和扁平化效应显著，截面呈现典型的“横鸭蛋”形式。翘曲将加剧隧道的纵向应变，增幅可达30%以上。隧道厚径比越小，其横纵向变形越大，增幅均可达30%左右。相较于法向土抗力系数，隧道纵向受力性能对切向土抗力系数更为敏感而横截面受力性能对双向抗力系数的敏感性一致。

关键词

土隧相互作用; 广义梁理论; 椭圆变形; 翘曲; 隧道横纵向力传递机制

地下隧道在受到地层长期沉降及外部荷载扰动作用下易发生纵向不均匀变形和截面的压扁^［1-2］。过大的纵向不均匀变形和横向收敛将引发隧道发生开裂、渗漏及螺栓脱落等一系列病害^［3-4］。这些病害产生的机制极为复杂且重要，为此，不少研究学者针对隧道横纵向变形的响应机制及效应的展开了深入的研究^［5-6］。

近年来，研究学者针对隧道纵向性能提出了大量的纵向结构模型。Koizumi等^［7］提出的梁-弹簧模型和Shiba等^［8］提出的纵向等效连续化模型则是研究盾构隧道纵向性能最常用的方法。纵向等效连续化模型将隧道视为刚度折减的纵向均质梁，因其计算简便且满足工程精度的特点得到了广泛的应用^［9-10］。与纵向等效连续模型相比，纵向梁-弹簧模型具有明确的物理意义，分别采用梁和弹簧来模拟管片和接头，能够充分考虑接头与管片之间的差异^［11-12］。除了盾构隧道的纵向受力性能，其横向受力性能同样是分析隧道结构安全的重要一环，且隧道横纵向性能往往是相互影响的^［13］。然而截至目前，关于隧道横纵向受力性能的研究，往往是将横向和纵向单独考虑，并未能充分考虑其相互耦合以及截面翘曲的效应。

大直径圆形隧道或异形隧道在纵向上发生弯曲与剪切变形的同时，在理论上必然伴随着横断面压弯和翘曲等复杂变形的出现^［11］，这对隧道结构安全是不利的。在高应力压载条件下，大直径隧道的横截面往往会产生压弯等畸变，进而呈现出典型的“鸭蛋”式扁平化变形。而由于变形协调关系，横截面的扁平化在纵向上会进一步引发显著的翘曲变形，即隧道中轴和两端分别沿纵向的不同方向变形，这势必将恶化隧道纵向变形状态。且随着隧道埋深和断面的不断加大、穿越地层的复杂化、隧道截面异型化以及建筑材料的高强薄壁化，由纵向不均匀变形导致的横截面效应将愈发显著。因此，该效应应得到更进一步的关注和研究。

扁平化和翘曲是薄壁结构稳定性和纵向弯曲设计的核心问题，在大型箱型桥梁等工程结构中得到了充分重视，而在隧道工程界尚认识不足。其原因可能是岩土介质中的中小直径隧道的扁平化和翘曲变形并不明显所致。但对于深埋的大直径隧道而言，翘曲和扁平化所产生的次生应力则可足以导致隧道开裂甚至局部失稳。不少地下隧道，特别是大直径隧道（D ≥ 11m），其管片厚度远小于其直径（t/D ≤ 1/20）往往可视为薄壁结构。针对薄壁圆环结构的特点，Schardt等^［14］在Vlasov梁的基础上提出了广义梁理论，该理论巧妙地用多个独立的截面变形模式来表征梁的总变形，进而快速求解荷载引发的结构响应。基于广义梁理论，不少研究学者针对直线型或曲线型的薄壁圆环结构在弯剪扭以及屈曲受力条件下的结构响应进行了研究，总结了圆形薄壁结构的多模态变形及受力状态，并论证了翘曲和扁平化的显著性^［15-16］。综上所述，为准确分析隧道的受力状态，考虑翘曲以及扁平化等变形对结构的影响极其重要。

为揭示复杂受力状态下隧道结构横纵向响应特点，本文提出了考虑截面扁平化和翘曲的横纵向土隧相互作用模型解析方法。基于广义梁理论，将隧道变形分解为弯曲、剪切、扁平和翘曲4个变形模态，同时为考虑土体和隧道之间的三维相对位移，引入三维地基模型。基于最小势能原理，推导了模型的有限元解。进一步结合两个工程案例验证了提出模型的有效性。最后通过参数化分析得到了关键设计参数对结构横纵向受力性能的影响。本文研究可为隧道精细化设计与调控提供理论基础。

1 横纵向土隧相互作用模型

1.1　基本假设

大量实测分析表明，隧道在周围土体长期变形、盾构穿越扰动、地裂缝以及隧道渗漏水的影响下，不仅仅会产生显著的纵向差异弯曲和剪切变形，同时在横断面会发生常见的压扁式“竖鸭蛋”和“横鸭蛋”的扁平化变形。此外，横断面的这种畸变将进一步引发纵向产生明显的翘曲效应。以上这类变形使得土隧相互作用存在复杂性，且过量的横纵向变形甚至会引起螺栓的脱漏和混凝土的剥落。为更精细的对复杂受力条件下的隧道结构响应进行分析，对这一系列三维空间的横纵向变形研究存在必要性。

本文基于广义梁理论^［14］，将隧道变形主要分解为4类，如图1所示：图1a弯曲：表现为隧道截面绕其中性轴旋转；图1b翘曲：表现为隧道中轴和两端沿纵向相背变形；图1c剪切：表现为隧道截面沿竖向产生位移；图1d扁平化：表现为隧道顶、底被挤扁，腰部向外变形。并作以下假设：

图 1  隧道变形的分解

Fig. 1  Decomposition of tunnel deformation

（1）隧道截面不可拉伸；

（2）隧道结构视为圆环薄壁结构（t<<D）；

（3）不考虑截面径向切应力和法向应力的影响；

（4）纵向纤维始终垂直于隧道轮廓中面；

（5）材料具有线性弹性和各向同性。

基于上述假设，根据广义梁理论，可将隧道变形分解为法向、切向和纵向分量来表示：

$$\left.\begin{array}{c}{u}_{\mathrm{s}}(\theta ,z)={\psi }_{\mathrm{s}}^{\beta }(\theta )\beta (z)+{\psi }_{\mathrm{s}}^U(\theta )U(z)+{\psi }_{\mathrm{s}}^V(\theta )V(z)+{\psi }_{\mathrm{s}}^F(\theta )F(z)\\ u_{\mathrm{n}}(\theta ,z)={\psi }_{\mathrm{n}}^{\beta }(\theta )\beta (z)+{\psi }_{\mathrm{n}}^U(\theta )U(z)+{\psi }_{\mathrm{n}}^V(\theta )V(z)+{\psi }_{\mathrm{n}}^F(\theta )F(z)\\ u_{\mathrm{z}}(\theta ,z)={\psi }_{\mathrm{z}}^{\beta }(\theta )\beta (z)+{\psi }_{\mathrm{z}}^U(\theta )U(z)+{\psi }_{\mathrm{z}}^V(\theta )V(z)+{\psi }_{\mathrm{z}}^F(\theta )F(z)\end{array}\right\}$$

（1）

式中：u_s、u_z、u_n分别为切向变形、纵向变形、法向变形；θ 是计算点与圆心的连线和竖直线的夹角；z为隧道纵向坐标；β、U、V、F分别为隧道弯曲、翘曲、剪切以及扁平化变形；${\psi }_j^i$为单位i变形引起的隧道j方向的变形，其表达式如下：

$$\left.\begin{array}{l}{\psi }_{\mathrm{s}}^{\beta }(\theta )=\mathrm{0}\\ {\psi }_{\mathrm{s}}^U(\theta )=\mathrm{0}\\ {\psi }_{\mathrm{s}}^V(\theta )=\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ {\psi }_{\mathrm{s}}^F(\theta )=\frac{\mathrm{1}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta \end{array}\right\}$$， $$\left.\begin{array}{l}{\psi }_{\mathrm{n}}^{\beta }(\theta )=\mathrm{0}\\ {\psi }_{\mathrm{n}}^U(\theta )=\mathrm{0}\\ {\psi }_{\mathrm{n}}^V(\theta )=-\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ {\psi }_{\mathrm{n}}^F(\theta )=-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta \end{array}\right\}$$， $$\left.\begin{array}{l}{\psi }_{\mathrm{z}}^{\beta }(\theta )=R\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ {\psi }_{\mathrm{z}}^U(\theta )=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta \\ {\psi }_{\mathrm{z}}^V(\theta )=\mathrm{0}\\ {\psi }_{\mathrm{z}}^F(\theta )=\mathrm{0}\end{array}\right\}$$

（2）

1.2　模型建立

基于前文假设，结合最小势能原理推导隧道横纵向土隧相互作用模型。将隧道视为弹性地基梁上的具有4变形模态的薄壁梁结构，结构受到荷载和三维土抗力的作用，如图2所示。利用C₀薄壁圆环单元，在静力分析中隧道势能表达式为

$$\mathrm{\Pi }=U_e+U_s+U_p$$

（3）

式中：U_e为隧道应变能；U_s为土隧相互作用做功；U_p为外荷载做功。

图 2  隧道纵向模型

Fig. 2  Longitudinal structural model of tunnel

则最小势能原理表达式为

$$\frac{\partial \prod }{\left\{\partial \mathrm{\Delta }\right\}}=\frac{\partial U_{\mathrm{e}}}{\left\{\partial \mathrm{\Delta }\right\}}+\frac{\partial U_{\mathrm{s}}}{\left\{\partial \mathrm{\Delta }\right\}}+\frac{\partial U_{\mathrm{p}}}{\left\{\partial \mathrm{\Delta }\right\}}=\mathrm{0}$$

（4）

式中：｛Δ｝为隧道结点变形矩阵。

进一步可以得到有限元的矩阵表达式为

$$\left\{\left[k\right]+\left[k_{\mathrm{s}}\right]\right\}\left\{\mathrm{\Delta }\right\}=F_{\mathrm{0}}+F_{\mathrm{p}}$$

（5）

式中：［k］为隧道结构刚度矩阵；［k_s］为土隧相互作用刚度矩阵；F₀为分布荷载矩阵；F_p为结点荷载矩阵。后文给出各参数的计算方法。

1.2.1　隧道结构刚度矩阵

根据所假设的隧道变形模式及几何方程，可进一步得到隧道的应变表达式为

$$\left.\begin{array}{l}{\epsilon }_{\mathrm{z}\mathrm{z}}=\frac{\partial u_{\mathrm{z}}}{\partial z}\\ {\epsilon }_{\mathrm{\theta }\mathrm{\theta }}=\mathrm{0}\\ {\gamma }_{\mathrm{\theta }\mathrm{z}}=\frac{\mathrm{1}}{R}\frac{\partial u_{\mathrm{n}}}{\partial \theta }+\frac{\partial u_{\mathrm{s}}}{\partial z}\\ {\kappa }_{\mathrm{\theta }\mathrm{\theta }}=\frac{\mathrm{1}}{R^{\mathrm{2}}}\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{n}}}{\partial {\theta }^{\mathrm{2}}}-\frac{\partial u_{\mathrm{s}}}{\partial \theta }\right)\end{array}\right\}$$

（6）

式中：ε_zz为纵向应力；ε_θθ为切向应力；γ_θz为剪应力；κ_θθ为截面内曲率；R为隧道半径。

则进一步可得到隧道应变能表达式为

$$U_{\mathrm{e}}=\frac{E_{\mathrm{1}}}{\mathrm{2}}\underset{l}{\int }\underset{A}{\int }{\epsilon }_{\mathrm{z}\mathrm{z}}^{\mathrm{2}}\mathrm{d}A\mathrm{d}z+\frac{G}{\mathrm{2}}\underset{l}{\int }\underset{A}{\int }{\gamma }_{\mathrm{\theta }\mathrm{z}}^{\mathrm{2}}\mathrm{d}A\mathrm{d}z+\frac{\mathrm{1}}{\mathrm{2}}{E}_{\mathrm{1}}I\underset{l}{\int }\underset{L}{\int }{\kappa }_{\mathrm{\theta }\mathrm{\theta }}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z$$

（7）

式中：E₁为折算弹性模量，$E_{\mathrm{1}}=\frac{E}{\mathrm{1}-{\mu }^{\mathrm{2}}}$，μ为泊松比；E为弹性模量；G为剪切模量；A为截面面积范围；L为截面轮廓范围，用θ表示其范围则为［0， 2π］；l为单元长度；I为惯性矩，I=t³/12。

令C₀单元的结点位移表达为

$$\left\{\mathrm{\Delta }\right\}={\left\{\begin{array}{cccccccc}{\beta }_{\mathrm{1}}& V_{\mathrm{1}}& U_{\mathrm{1}}& F_{\mathrm{1}}& {\beta }_{\mathrm{2}}& V_{\mathrm{2}}& U_{\mathrm{2}}& F_{\mathrm{2}}\end{array}\right\}}^{\mathrm{T}}$$

（8）

式中：T表示转置。

根据隧道变形表达式（1），隧道的应变表达式可进一步表达为

$$\left.\begin{array}{l}{\epsilon }_{\mathrm{z}\mathrm{z}}={\psi }_z^{\beta }{N}^{\text{'}}(z)\left\{\beta \right\}+{\psi }_z^{\beta }{N}^{\text{'}}(z)\left\{U\right\}\\ {\gamma }_{\mathrm{\theta }\mathrm{z}}=\frac{\mathrm{1}}{R}\left({\psi }^{\text{'}}{}_z^{\beta }N(z)\left\{\beta \right\}+{\psi }^{\text{'}}{}_z^{\beta }N(z)\left\{U\right\}\right)+{\psi }_s^V{N}^{\text{'}}(z)\left\{V\right\}+{\psi }_s^F{N}^{\text{'}}(z)\left\{F\right\}\\ {\kappa }_{\mathrm{\theta }\mathrm{\theta }}=\frac{\mathrm{1}}{R^{\mathrm{2}}}\left({\psi }^{″}{}_n^F-{\psi }^{\text{'}}{}_s^F\right)N(z)\left\{F\right\}\end{array}\right\}$$

（9）

式中：$\left\{\beta \right\}=\left[N(z)\right]\left\{{\beta }_i\right\}$，［N（z）］为形函数，$N(z)=\left[\begin{array}{cc}\mathrm{1}-\frac{z}{l}& \frac{z}{l}\end{array}\right]$；同理｛U｝、｛V｝、｛F｝；${\psi }^{\text{'}}{}_i^j$表示${\psi }_i^j$对θ求导。

将式（9）代入到式（7），进一步根据最小势能原理对位移求偏导可得到隧道结构的刚度矩阵为

$$\left[ks\right]=\left[\begin{array}{cc}{k}_{\mathrm{11}}& k_{\mathrm{12}}\\ k_{\mathrm{21}}& k_{\mathrm{22}}\end{array}\right]$$

（10）

式中：

$$\begin{array}{l}\left[k_{\mathrm{11}}\right]=\left[\begin{array}{cccc}{E}_{\mathrm{1}}{a}_{\mathrm{1}}/l+G{b}_{\mathrm{1}}l/\mathrm{3}& -G{b}_{\mathrm{1}}/\mathrm{2}& & \\ -G{b}_{\mathrm{1}}/\mathrm{2}& G{b}_{\mathrm{1}}/l& & \\ & & E_{\mathrm{1}}{a}_{\mathrm{2}}/l+G{b}_{\mathrm{2}}l/\mathrm{3}& -G{b}_{\mathrm{3}}/\mathrm{2}\\ & & -G{b}_{\mathrm{3}}/\mathrm{2}& E_{\mathrm{1}}{c}_{\mathrm{1}}l/\mathrm{3}+G{b}_{\mathrm{4}}/l\end{array}\right]\\ \left[k_{\mathrm{12}}\right]={\left[k_{\mathrm{21}}\right]}^T=\left[\begin{array}{cccc}G{b}_{\mathrm{1}}l/\mathrm{6}-E_{\mathrm{1}}{a}_{\mathrm{1}}/l& G{b}_{\mathrm{1}}/\mathrm{2}& & \\ -G{b}_{\mathrm{1}}/\mathrm{2}& -G{b}_{\mathrm{1}}/l& & \\ & & G{b}_{\mathrm{2}}l/\mathrm{6}-E_{\mathrm{1}}{a}_{\mathrm{2}}/l& G{b}_{\mathrm{3}}/\mathrm{2}\\ & & -G{b}_{\mathrm{3}}/\mathrm{2}& E_{\mathrm{1}}{c}_{\mathrm{1}}l/\mathrm{6}-G{b}_{\mathrm{4}}/l\end{array}\right]\\ \left[k_{\mathrm{22}}\right]=\left[\begin{array}{cccc}{E}_{\mathrm{1}}{a}_{\mathrm{1}}/l+G{b}_{\mathrm{1}}l/\mathrm{3}& G{b}_{\mathrm{1}}/\mathrm{2}& & \\ G{b}_{\mathrm{1}}/\mathrm{2}& G{b}_{\mathrm{1}}/l& & \\ & & E_{\mathrm{1}}{a}_{\mathrm{2}}/l+G{b}_{\mathrm{2}}l/\mathrm{3}& G{b}_{\mathrm{3}}/\mathrm{2}\\ & & G{b}_{\mathrm{3}}/\mathrm{2}& E_{\mathrm{1}}{c}_{\mathrm{1}}l/\mathrm{3}+G{b}_{\mathrm{4}}/l\end{array}\right]\end{array}$$

式中：a₁=πR³t；a₂=πRt；b₁=πRt；b₂=4πt/R；b₃=πt；b₄=Rπt/4；c₁=3πt³/R³。

1.2.2　土隧相互作用刚度矩阵

隧道的4变形模态将引起土隧之间产生纵向、切向以及法向等变形，必然使得土隧相互作用呈现显著的空间分布特性。本文在推导过程中，引入法向和切向土抗力来考虑三维土隧相互作用力。并假设土抗力与土隧相对位移服从温克尔地基梁理论。则土抗力做功可分为3部分，如图3所示：

图 3  土抗力分解图

Fig. 3  Decomposition of soil reaction

（1）由隧道沉降或隆起和拱腰水平变形引起的法向土抗力作用；主要影响区域为隧道与上/下卧层接触的区域，假设超载条件下计算区间为外轮廓线的2/3（见图 3a中（1）所示）、卸载条件为外轮廓线的1/3（见图 3b中（1）所示），令其土抗力引起的能量表达为U_s，1；

（2）由隧道沉降或隆起和扁平化引起的切向土抗力作用，如图 3a和3b中（2）所示，令其土抗力引起的能量表达为U_s，2；

（3）隧道纵向弯曲和翘曲引起的纵向土抗力作用，如图 3a和图 3b中（3）所示，令其土抗力引起的能量表达为U_s，3。

则可以将土抗力引起的总能量表达为

$$U_{\mathrm{s}}=U_{\mathrm{s},\mathrm{1}}+U_{\mathrm{s},\mathrm{2}}+U_{\mathrm{s},\mathrm{3}}$$

（11）

根据3类能量的定义，以超载形式为例，可进一步得到各分量表达式为

$$\left.\begin{array}{l}{U}_{\mathrm{s},\mathrm{1}}=\frac{k_{\mathrm{n}}}{\mathrm{2}}\underset{l}{\int }\underset{\mathrm{2}L/\mathrm{3}}{\int }{u}_{\mathrm{n}}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z=\frac{k_{\mathrm{n}}}{\mathrm{2}}\underset{l}{\int }\underset{\mathrm{2}L/\mathrm{3}}{\int }{\left({\psi }_{\mathrm{n}}^{V}V+{\psi }_{\mathrm{n}}^{F}F\right)}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z\\ U_{\mathrm{s},\mathrm{2}}=\frac{k_{\mathrm{t}}}{\mathrm{2}}\underset{l}{\int }\underset{L}{\int }{u}_{\mathrm{s}}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z=\frac{k_t}{\mathrm{2}}\underset{l}{\int }\underset{L}{\int }{\left({\psi }_{\mathrm{s}}^{V}V+{\psi }_{\mathrm{s}}^{F}F\right)}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z\\ U_{\mathrm{s},\mathrm{3}}=\frac{k_{\mathrm{t}}}{\mathrm{2}}\underset{l}{\int }\underset{L}{\int }{u}_{\mathrm{z}}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z=\frac{k_t}{\mathrm{2}}\underset{l}{\int }\underset{L}{\int }{\left({\psi }_{\mathrm{z}}^{\beta }\beta +{\psi }_{\mathrm{z}}^{U}U\right)}^{\mathrm{2}}\mathrm{d}s\mathrm{d}z\end{array}\right\}$$

（12）

式中：k_n为法向土弹簧刚度；k_t为切向土弹簧刚度。

将式（1）代入到（12），进一步根据最小势能原理对位移求偏导可得到土抗力刚度矩阵为

$$\left[ks\right]=\left[\begin{array}{cc}k{s}_{\mathrm{11}}& k{s}_{\mathrm{12}}\\ k{s}_{\mathrm{21}}& k{s}_{\mathrm{22}}\end{array}\right]$$

（13）

其中

$$\begin{array}{l}\left[k{s}_{\mathrm{11}}\right]=\left[k{s}_{\mathrm{22}}\right]=\frac{l}{\mathrm{3}}\left[\begin{array}{cccc}{k}_{\beta }& & & \\ & k_{v\mathrm{1}}+k_{v\mathrm{2}}& & k_{vF}\\ & & k_U& \\ & k_{vF}& & k_{F\mathrm{1}}+k_{F\mathrm{2}}\end{array}\right]\\ \left[k{s}_{\mathrm{12}}\right]={\left[k{s}_{\mathrm{21}}\right]}^T=\frac{\mathrm{1}}{\mathrm{2}}\left[k{s}_{\mathrm{11}}\right]\end{array}$$

式中：k_&#x1d6fd; = πR³k_t；k_U = πRk_t；k_v1=$\frac{\mathrm{8}\pi -\mathrm{3}\sqrt[]{\mathrm{3}}}{\mathrm{12}}R{k}_{\mathrm{n}}$；k_v2= πRk_t；k_vF=$\frac{-\sqrt[]{\mathrm{3}}}{\mathrm{2}}R{k}_{\mathrm{n}}$；k_F2=$\frac{\mathrm{16}\pi +\mathrm{3}\sqrt[]{\mathrm{3}}}{\mathrm{24}}R{k}_{\mathrm{n}}$；k_F2=Rπk_t/4。

1.2.3　荷载及边界条件

荷载做功所产生的系统能量可表达为

$$U_{\mathrm{p}}=-\left(\underset{l}{\int }\underset{L}{\int }q{u}_{\mathrm{n}}\mathrm{d}s\mathrm{d}z+\underset{l}{\int }\underset{L}{\int }q{u}_{\mathrm{s}}\mathrm{d}s\mathrm{d}z+\underset{l}{\int }\underset{L}{\int }q{u}_{\mathrm{z}}\mathrm{d}s\mathrm{d}z\right)$$

（14）

式中：q为截面内的分布荷载，式中荷载和位移为矢量点乘。

由此，可进一步得到荷载矩阵的表达式为

$${\left\{F\right\}}_i=-\frac{\partial U_{\mathrm{p}}}{\partial \left\{i\right\}}=\underset{l}{\int }\underset{L}{\int }p{\left[N(z)\right]}^{\mathrm{T}}{\psi }_i(s)\mathrm{d}s\mathrm{d}z=\left[\begin{array}{c}\underset{l}{\int }\underset{L}{\int }\stackrel{\rightharpoonup }{q}(\mathrm{1}-\frac{z}{l}){\phi }_i(s)\mathrm{d}s\mathrm{d}z\\ \underset{l}{\int }\underset{L}{\int }\stackrel{\rightharpoonup }{q}\frac{z}{l}{\phi }_i(s)\mathrm{d}s\mathrm{d}z\end{array}\right]$$

（15）

式中：i为任一自由度i。

将隧道远端车站视为固定约束，联立各刚度矩阵（10）、（13）和荷载矩阵（15），则可得到隧道横纵向变形和内力的分布情况，包括沉降以及收敛变形等。进一步可得到隧道纵向应变和由于隧道截面扁平化而导致的截面弯矩，其表达式如下

$$\left.\begin{array}{l}\left\{{\sigma }_{\mathrm{z}}\right\}=E_{\mathrm{1}}\frac{\partial u_{\mathrm{z}}}{\partial z}=E_{\mathrm{1}}(N^{\text{'}}(z)\left\{\beta \right\}{\phi }_{\mathrm{z}}^{\beta }(s)+\\ N^{\text{'}}(z)\left\{U\right\}{\phi }_{\mathrm{z}}^U(s))\\ \left\{{\tau }_{\mathrm{s}}\right\}=G(\frac{\partial u_{\mathrm{z}}}{R\partial \theta }+\frac{\partial u_{\mathrm{s}}}{\partial z})=G[\frac{\mathrm{1}}{R}({\psi }^{\text{'}}{}_{\mathrm{z}}^{\beta }(N(z)\left\{\beta \right\}+\\ {\psi }^{\text{'}}{}_{\mathrm{z}}^U(N(z)\left\{U\right\})+{\psi }_{\mathrm{s}}^V(N^{\text{'}}(z)\left\{V\right\}+\\ {\psi }_{\mathrm{s}}^F(N^{\text{'}}(z)\left\{F\right\}]\\ M(\theta ,z)=E_{\mathrm{1}}I{\kappa }_{\mathrm{\theta }\mathrm{\theta }}=\frac{E_{\mathrm{1}}{t}^{\mathrm{3}}}{\mathrm{4}{R}^{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\theta \right)N(z)\left\{F\right\}\end{array}\right\}$$

（16）

2 模型验证

为进一步验证本文提出的有限元解的优越性和合理性，本节结合两例工程案例，将提出的方法与现有的弹性地基梁上的欧拉伯努利梁解以及ABAQUS模型进行对比，分析在给定工程案例下隧道的受力和变形的性能。

2.1　盾构穿越既有隧道案例

（1）工程案例

新隧道施工易导致附近既有隧道产生纵向不均匀变形和截面的扁平化。Lin 等^［17］基于长沙地铁案例，通过数值模拟研究了不同角度盾构下穿对既有隧道的影响。案例中，既有隧道外径6.0 m，厚度0.3 m，新建隧道与即有隧道之间的最小垂直距离为2.9m，穿越地层为中粗砂层。以该案例中盾构30°斜角下穿既有隧道为例，验证本文方法的有效性。

（2）对比验证

图4为案例结果、弹性基础上的欧拉伯努利梁法和本文所提出方法的隧道变形对比图。本文方法与案例实际结果吻合度高，说明了本文方法的有效性。此外，从图中可以看出，随着计算截面距新隧道中心线距离的增加，既有隧道的断面由“竖鸭蛋”（隧道底部变形大于顶部变形）变为“横鸭蛋”（隧道顶部变形大于底部变形），本文提出的方法能很好的体现出这种截面的变化，且与案例结果吻合度高。相较于始终将隧道截面视为圆形的欧拉伯努利梁解，新方法可以充分考虑截面的畸变形式，使得计算结果更符合实际状态，说明了新方法的优越性和可靠性。

图 4  结果对比图

Fig. 4  Comparison of tunnel deformation

2.2　受地面超载影响的隧道

（1）工程案例

以上海某公路隧道为工程背景，隧道外径为11m，管片厚度0.48n。隧道埋深6m，地表存在局部堆载 （如图5所示）。隧道上方超载作用区域大小为10m×10 m，荷载压力为100 kPa。该处地层法向抗力系数为4MPa·m^-1，切向抗力系数为1MPa·m^-1。

图 5  盾构隧道受地面超载简图

Fig. 5  Shield tunnel subjected to ground surcharge

（2）ABAQUS模型对比验证

以该公路隧道为工程背景，分别利用温克尔地基上欧拉伯努利梁解析解、ABAQUS有限元模型和本文提出的方法对（1）中超载作用下的隧道内力和变形进行计算，两端均取为固定约束，将三者的结果进行对比分析，进一步验证所提出模型的有效性。其中，在利用本文所提出的有限元解计算过程中，在综合计算的效率和精度的基础上，将隧道被分割为单位长度取为0.5m的圆环单元，该单元长度相较于更小计算单元具有更快的计算效率和足够的精度。隧道总计算长度取为200m（如图6a所示）。ABAQUS有限元模型中隧道采用壳单元，土抗力用三维土弹簧来模拟，荷载则直接施加于隧道的外表面上（如图6b所示）。

图 6  有限元模型图

Fig. 6  FEM model

图7为3种方法的计算结果对比图，可以发现本文方法的计算结果和ABAQUS输出解吻合较好，验证了本文方法的有效性。对比欧拉伯努利梁解，可发现该传统方法能得到和本文方法精确一致的内力解，所计算的沉降变形与本文提出方法计算的平均位移一致。但传统方法不能全面地计算出隧道的截面变形情况，特别是无法解释隧道截面压扁和翘曲的效应。而本文提出的方法能很好的反映了隧道的截面变形情况，能得到与有限元模型输出结果一致的全截面的变形情况，且计算相较于ABAQUS模型更为简便，说明了该方法的优越性。同时本文提出的方法相较于有限元模拟软件存在以下优势：避免繁琐的有限元软件建模过程且物理含义清晰，仅需在MATLAB等计算软件在已搭建刚度和荷载矩阵的基础上，改变隧道参数和荷载参数，即可快速得出结果。

图 7  既有方法与本文方法的对比图

Fig. 7  Comparison of published solution and proposed results

根据本文的计算结果可知在地面超载的作用下，隧道截面呈现典型的“横鸭蛋”的形式，这将显著增大隧道顶、底部的变形和前、后截面的错台。弯矩、截面水平收敛和沉降呈峰状的对称分布，剪切力呈双峰反对称分布，拐点出现在加载区边缘附近。

2.3　翘曲和扁平化效应

现有研究往往将隧道视为弹性地基上的欧拉伯努利梁或者铁木辛柯梁^［9-10］，该类模型不能考虑翘曲和扁平化的效应。为了进一步探究翘曲和扁平化对结构的影响。提取了案例2的隧道纵向应变的纵向分布以及截面弯矩的分布情况，如图8所示。由图8a可以看出，隧道的纵向应变主要由纵向弯曲和纵向翘曲两部分组成，其中翘曲约占总应变的1/3，说明了翘曲的显著性。最大应变出现在载荷区的中心，约为极限拉应变的2/3。图8b展示了最大扁平化截面的截面弯矩分布情况，隧道的扁平化引起截面呈现沿腰部对称的分布，最大截面弯矩可达约80kN·m·m^-1。这说明扁平化不仅仅引起额外的沉降变形和错台，同时引起较大的截面收敛变形和截面弯矩，对隧道安全极其不利。因此在高受载条件下土-隧相互作用中需进一步考虑截面扁平化和翘曲的效应

图 8  翘曲和扁平化引起的隧道响应

Fig. 8  Structural responses of tunnel induced by warping and flattening

2.4　讨论

从结构刚度矩阵（10）中可以发现隧道横纵向的耦合关系表现为纵向翘曲和截面畸变的相互耦合以及纵向弯曲与截面剪切的耦合，因此可将外荷载的形式，将隧道的变形分为3种模式：

（1） 当隧道受到的为截面均布的竖向力作用时（如重力作用），截面会同时发生显著的弯曲和剪切变形，而不会发生翘曲和畸变；

（2） 当隧道受到截面上下反对称均布的竖向力作用时（如上、下地层以及侧向土体的对等挤压变形），截面在产生显著的扁平化变形的同时会在纵向引起翘曲变形，而不会发生弯曲和剪切；

（3） 当隧道受到截面除以上例外的任意荷载时（附近开挖扰动以及土地长期沉将等），将同时发生显著的截面剪切和扁平化，同时伴随者纵向弯曲和翘曲的发生。

除此之外，之前的研究中往往将隧道纵向和横向割裂开来而单独考虑，而从刚度矩阵（10）可以看出纵向翘曲和截面扁平化、纵向弯曲和截面剪切是相互影响并耦合的，因此对隧道横纵向性能的研究应充分考虑横纵向相互耦合的作用。

为进一步对隧道横纵向变形和内力的传递机制进行分析，提取了隧道纵横向变形和内力分布情况：从图 9a可以得到隧道纵向变形由翘曲和弯曲组成，其中弯曲导致的纵向变形约为翘曲导致的2倍，且呈现为沿加载中心的反对称分布，但各自的最大变形点错位；而隧道截面变形由剪切和扁平化组成，其中剪切导致的截面变形约为扁平化导致的1.5倍，且均表现为沿加载中心的对称分布。同时可以发现各纵向变形的峰值位置对应了截面变形的反弯点位置；截面变形的峰值位置对应了纵向变形的零值位置，即横纵向变形传递的规律类似于原函数和导数关系。从图 9b中可以看出，横截面最大弯矩与纵向弯矩分布规律一致，均表现为沿加载中心对称的中间大，两边小的对称分布，即横纵向内力传递表现为同步变化的规律。

图 9  横纵向内力和变形传递

Fig. 9  Load and deformation transfer in longitudinal and transverse direction

3 参数化分析

为探析隧道横纵向受力性能的变化规律，便于工程应用，以下进一步对隧道关键设计参数的影响进行分析。参数化分析均基于2.1中的案例2为背景。

3.1　隧道厚径比的影响

图10为不同隧道宽厚度比隧道的计算变形和内力。

图 10  不同隧道厚径比的影响

Fig. 10  Effect of different thick-to-diameter ratios

计算过程中仅改变隧道厚度，计算厚径比从0.03到0.05的隧道结构响应变化。其中隧道沉降、纵向弯矩、纵向剪切力表示了隧道的弯剪效应，收敛变形和截面弯矩表示了扁平化效应的变化，纵向应变代表了隧道的弯曲和翘曲的变化情况。

可以看出，随着隧道厚度的增加，横纵向内力均增大，其中纵向弯矩和剪切力增幅最大约为10%，但截面弯矩增幅可达2倍以上，说明截面的扁平化与隧道的厚径比密切相关，这也与厚度决定截面的惯性矩大小有关。除此之外，随着隧道厚径比的减小，隧道横纵向变形随之增大，增幅均可达30%左右，说明纵向翘曲在厚径比越小的隧道如大直径隧道或异型隧道越为明显，因而在隧道设计中需重视翘曲的效应。

3.2　土抗力系数的影响

图11为不同土抗力系数下隧道的变形和内力计算结果。计算中分别改变切向土抗力系数和法向土抗力系数，其中刚度变化范围为2MPa·m^-1至10MPa·m^-1。从图中可以看出，切向刚度和法向刚度对隧道横纵向的受力性能的影响程度是不同的。于纵向受力性能而言，隧道纵向沉降以及内力对切向土抗力系数的变化更为敏感；随着切向土抗力系数增大5倍，变形将缩减为原来的50%。而对横向受力性能而言，横截面变形以及弯矩对双向抗力系数的敏感性是一致的，均表现为随着抗力系数的增大，收敛变形和截面弯矩均减小，减幅可达30%。因此在隧道设计中，除了需掌握法向土抗力系数的大小，也需要考虑切向土抗力系数的影响。

图 11  不同土抗力系数的影响

Fig. 11  Effect of different soil resistance coefficient

4 结论

本文通过理论推导得到了考虑截面扁平化和翘曲的土-隧相互作用模型，经过案例验证及其参数化分析，得到如下结论：

（1）基于广义梁理论，将隧道视为在弹性地基梁上的具有弯曲、翘曲、剪切和扁平化4变形模态的薄壁圆环结构。进一步采用最小势能原理，推导了任意外荷作用下的隧道横纵向内力和变形的计算模型。该模型可充分考虑隧道的横纵向相互耦合作用，其中隧道横纵向的耦合关系表现为纵向翘曲和截面扁平化的相互耦合以及纵向弯曲与截面剪切的耦合。

（2）新模型考虑了三维土隧相互作用，通过与既有模型比较，验证了本文方法的有效性和优越性。在高应力受载条件下，大直径隧道翘曲和扁平化效应不可忽略。扁平化导致截面呈现典型的“横鸭蛋”和“竖鸭蛋”形式；而翘曲将加剧隧道的纵向应变，增幅可达30%以上。

（3）基于所提出的新模型，可根据外荷载的形式将隧道的受力状态分为3类：弯剪；翘曲扁平；弯剪翘曲扁平。并总结出了横纵向变形传递的规律类似于原函数和导数关系以及横纵向内力传递表现为同步变化的规律。

（4）隧道的横纵向性能与截面厚径比密切相关。随着隧道厚度的增加，横纵向内力增大而变形减小，其中纵向内力增幅最大可达约10%，截面弯矩增幅可达2倍以上，变形减幅均可达30%。厚径比越小的隧道，纵向畸变和扁平化效应越为明显，应在设计中尤其注意。

（5）土抗力系数对隧道结构响应影响显著，且法向和切向抗力系数影响程度不同。隧道纵向受力性能对切向土抗力系数的变化更为敏感；随着切向土抗力系数增大5倍，变形将缩减为原来的50%。横截面受力性能对双向抗力系数的敏感性一致，随着抗力系数的增大，收敛变形和截面弯矩均减小，减幅可达30%。

值得指出的是，本文仅对圆形隧道横断面翘曲变形的影响等进行了探究，而异型隧道（矩形、类矩形、马蹄形等）隧道横断面的扁平化和翘曲变形将更为显著，但其本质相同，将在后续论文中进行详细探讨。

作者贡献声明

刘颖彬：论文构思、写作和修改，理论推导；

廖少明：论文构思、写作和修改，基金获取；

肖明清：基金获取、论文修改；

门燕青：基金获取、论文修改。