单位圆外域到任意开挖断面隧洞外域共形映射的计算方法

祝江鸿，杨建辉，施高萍，王 珺，蔡建平

（1. 浙江科技学院 建工学院，杭州 310023；2. 浙江水利水电专科学校，杭州 310018；3. 北京科技大学 土木与环境工程学院，北京 100083）

1 引 言

通过对地下隧洞围岩力学分析中的解析方法及其结果的分析，往往可以获得围岩应力和位移分布及稳定性规律的认识，因此它的理论价值是不容忽视的，众多学者为此展开了大量的研究。然而，受制于地下洞室开挖断面形状的制约，地下隧洞工程力学分析中的解析解成果主要集中在不同边界条件下的圆形和椭圆形洞室方面，如Fahimifar等[1]提出了静水应力场下伯格斯黏弹性体中圆形开挖断面隧洞的解析解，Gao等[2]开展了高寒区圆形隧道的弹塑性应力解析解研究，Li等[3]对深埋圆形隧道在地压和内衬压力作用下的弹性应力和位移场进行了分析，Verruijt[4]研究了无限半平面内圆形洞室的弹性变形机制，Zhang等[5]提出了横观各向同性岩石中深埋椭圆形隧道的位移解析解。相比较而言，对其他开挖断面形状的地下洞室围岩力学解析分析，则在数学和力学上要复杂得多，原则上可采用弹性力学的复变函数方法获得其平面弹性问题的解析解[6]。利用复变函数开展地下任意开挖断面隧洞围岩力学分析时，需获得符合精度要求的映射函数及用来表示围岩应力和位移的两个解析函数[7]。关于围岩应力复变函数分析法中的两个解析函数求解，笔者已在另文进行了研究[8]，本文对任意开挖断面隧洞的映射函数求解开展研究。

基于复变函数法的地下隧洞围岩力学分析，需将开挖边界线内域（外域）看成由单位圆内域（外域）共形映射而得[6－9]，将以复杂开挖断面为边界的隧洞力学分析问题转化成以单位圆周线为边界的力学问题。实际上，单位圆内（外）域到任意封闭曲线内（外）域之间的共形映射函数构建在数学理论上较完整，有 Schwarz-Christoffel积分公式[10－11]、有理分式[12－13]、泰勒级数[14－15]等模型。泰勒级数和洛朗级数等级数形式是隧洞映射函数求解的常用数学模型。皇甫鹏鹏等[16]以洛朗级数建立起单位圆外域-洞室外域的映射函数，提出了基于边界点搜索的映射函数求解方法，根据等距离比的原则解决了计算点的调整问题，然而对整个迭代过程实现及对应点组数和映射函数构成项数之间的关系交代不清晰。朱大勇等[17]同样以洛朗级数建立起单位圆外域-洞室外域的映射函数，利用三角函数系的正交性获取了函数中各项系数的表达式，阐明了计算点调整方法和系数求解的迭代过程，不足是对于复杂的开挖洞形，映射函数构成项数则多达200多项，计算量大且精度低。Zhu等[18]开展了基于梅林谢耶夫法的映射函数近似求解研究，然而给出的映射函数精度较低。王润富[19]以泰勒级数的有限项建立起单位圆内域-孔洞内域的映射函数，但并未明确对应点的调整方法及映射边界线与实际开挖边界线误差的界定。Deillo等[20]利用Faber级数法开展了单位圆内域-目标区域的共形映射的数值计算研究，然而映射函数中项数构成复杂、计算量偏大。Schwarz-Christoffel积分公式亦是地下隧洞映射函数的常用数学模型，但对应点的确定和被积函数中参数的计算非常复杂，有时候甚至难以完成。范广勤等[21]采用该积分公式求解了轴对称六角形洞室的映射函数，具体求解时用3个绝对收敛的幂级数在其收敛域内相乘来近似被积函数，这样处理的优点是将原积分转化成了对多项式积分，使积分求解成为可能，不足是 3个幂级数相乘展开后的系数求解依然繁琐，二是之所以可以用3个幂级数相乘代替被积函数，是因为其针对特定的轴对称六角形，适用范围有限。郑志强[22]也采用 Schwarz-Christoffel积分建立了单位圆与多边形的映射函数，研究了被积函数中未知参数的数值计算方法，解决了参数求解和参数精度问题，但获取参数只是解决问题的第一步，即明确了被积函数的表达式，被积函数的积分过程依然复杂。

综上所述，尽管地下任意开挖断面隧洞映射函数的求解原理在数学上较完整，但在计算方法的实现上，尤其在映射函数系数求解、对应点确定、计算点（映射点）调整及映射边界线误差控制等方面不甚成熟，所获得的映射函数存在项数构成复杂、计算量大及映射精度偏低等问题。本文根据黎曼存在定理，建立了以洛朗级数有限项表示的单位圆外域到任意开挖断面隧洞外域的映射函数，根据边界对应定理将域到域的映射函数求解问题转化为单位圆周线到任意开挖断面隧洞边界线的求解问题，并从三角插值理论出发，采用奇偶插值点反复相互迭代，对映射函数系数求解进行了研究。采用法线逼近法实现了计算点（映射点）的调整，并给出了映射函数系数迭代计算的收敛条件。利用本文的研究成果可方便地计算单位圆外域到任意开挖断面隧洞外域的共形映射函数。

2 映射函数建立

在复变函数理论中，规定了复平面中有一个惟一的“无穷远点”，且该点与复球面上的北极N相对应，因此，对于图1(a)和图1(b)中给出的区域D和D′均为单连通域。对两个给定的单连通域，由黎曼存在定理可知，总存在惟一映射函数将一个区域共形映射到另一个区域[23－24]。由此可知，将以单位圆周线Γ为边界的外部区域D共形映射到以复杂边界线Γ′为边界的隧洞外部区域D′的映射函数必存在（图1），且洛朗级数是满足该映射条件的数学模型[6－7,9]，可由式（1）给出。

图1 单位外域到任意开挖断面隧洞外域的共形映射Fig.1 Conformal mapping from the exterior of unit circle to the exterior of cavern with arbitrary excavation cross-section

式中： Cn=An+i Bn，An、Bn均为实常数。

根据边界对应定理，若式（1）能将图1(a)中的单位圆周线Γ共形映射到图 1(b)中的复杂边界线Γ′，即映射函数只需满足这两条边界线之间的共形映射关系，则该函数定能将D区域保角变换到D′区域[23－24]。据此，单位圆外域共形映射到任意开挖断面隧洞外域的映射函数求解转化成单位圆周线到任意开挖断面隧洞边界线的映射函数求解，即研究对象由域-域转换到单位圆周线-任意开挖断面隧洞边界线上。

单位圆周线Γ上点的坐标以三角函数表示，由式（2）给出：

复杂边界线Γ′上点的坐标可由式（3）给出：

由于式（1）是将单位圆周线共形映射到任意开挖断面隧洞边界线的映射函数，将式（2）、（3）代入式（1）后，由实部虚部相等原则可得

若式（1）已知，即An、Bn已知，则由式（4）、（5）可知，单位圆周线与任意开挖断面隧洞边界线之间的映射对应关系就完全确定，由此单位圆外域到任意开挖断面隧洞外域之间的映射对应关系也就完全确定。下面根据单位圆周线与任意开挖断面隧洞边界线之间的关系研究An和Bn的求解过程。

3 映射函数求解

映射函数是以洛朗级数的形式给出的，在实际求解时只能取有限项作近似。假设映射函数由洛朗级数的前m项构成，则式（1）可改写为

式中： j = 1, 2,… ,m。式（7）、（8）可根据式（4）、（5）及三角函数系的正交性推导而得，具体的推导过程可参见文献[6，14－15]，这里不再赘述。事实上，要知道单位圆周线Γ和复杂边界线Γ′上m对对应点是困难的，Aj和Bj的求解可采用单位圆上奇偶插值点相互反复迭代计算实现。

3.1 奇偶插值点相互迭代计算

将ζ平面中的单位圆周线Γ作 2m等分，2m个等分点分为奇数点和偶数点两组，相邻两奇数点的幅角差和相邻两偶数点幅角差Δφ均为2π/m，奇数点和偶数点的幅角φ2k和φ2k+1可由下式计算：

式中： k = 0, 1,… ,m - 1。在z平面中的Γ′上任意取m个点，假设这些点与ζ平面中Γ上的偶数点对应，则根据式（7）～（9）可得式（11）、（12），据此求得一组和（系数上标的“0”表示由偶数对应点求得的，“1”表示由奇数对应点求得）。

3.2 计算点（映射点）的调整

图2 计算点调整Fig.2 Calculating points adjustment

3.3 迭代计算的收敛条件

按前述步骤进行映射函数求解时，有迭代计算的收敛问题。这与工程的计算精度要求有关，也与迭代计算的有效性有关。一般来说，对由复杂边界线构成的区域，要获得精确的映射函数是有困难的，实际求解时只能获得符合计算精度要求的近似映射函数。因此，利用本文的解法获得的映射函数进行隧洞映射时，得到的映射断面与实际开挖断面是有误差的。只要误差控制在工程计算精度要求范围内，则利用该映射函数进行围岩应力和位移计算时，所得结果可满足围岩力学分析要求。因此，需对本文的映射函数求解设置精度要求，即收敛条件。本文以隧洞映射边界线与实际开挖边界线之间的最大距离小于某一数值ε作为迭代计算的收敛条件：

4 算 例

地下隧洞断面开挖边界线的设计线型通常有直线和弧线。为更一般性地验证本文研究结果的有效性，选取由直线段构成开挖边界线的矩形隧洞、以弧线段构成开挖边界线的公路两车道隧洞以及由直线段与弧线段共同组成开挖边界线的复杂形隧洞作为研究对象，开展将单位圆外域共形映射到上述 3种开挖断面隧洞外域的映射函数求解。利用Matlab语言将3个算例的映射函数求解过程编制成计算代码，整个计算过程由计算机自动完成，其中，迭代计算的收敛条件均采用式（16），按映射边界线与实际洞形边界线之间的绝对误差ε小于5 mm考虑。

4.1 矩形隧洞

矩形隧洞的几何尺寸及在坐标系中的位置如图3所示。映射函数的项数构成按9项考虑，因此在实际迭代计算时单位圆与实际边界线之间的对应点只需9对（考虑到奇偶插值点迭代计算的需要，单位圆周线作18等分）。图4给出了迭代1次、5次和16次后的映射断面图形，图形的绘制是按各轮迭代后求得的Aj和Bj及单位圆周线上的 100个点代入式（6）后得到的，由计算程序自动给出。其中图4(c)为符合精度要求的映射断面，图中各点均作了误差验算且满足精度要求，其映射函数系数见表1。

图3 矩形隧洞Fig.3 A rectangular tunnel

图4 矩形洞室的映射断面Fig.4 Mapping cross-sections of rectangular tunnel

表1 矩形隧洞映射函数系数Table 1 Coefficients of the mapping function for rectangular tunnel

4.2 两车道公路隧洞

以《公路隧道设计规范》(JTG D702004)中的两车道公路隧道标准内轮廓断面（图 5(a)）为例，进行映射函数求解。隧道实际开挖断面一般按径向方向比标准内轮廓线大50～500 mm左右考虑开挖尺寸，具体视围岩等级、埋深情况及支护设计等条件确定。本次计算的开挖尺寸按径向方向比设计的内轮廓线大300 mm设置计算尺寸（图5(b)）。映射函数的项数构成按16项考虑，为求解Aj和Bj，实际迭代计算时单位圆和实际边界线之间的对应点需16对，而单位圆线作了32等分。

图5 公路两车道隧洞断面Fig.5 Cross sections of two-lane road tunnel

图6给出了迭代1次、2次和6次后的映射断面，图形的绘制方法同算例 1。其中，图 6(c)为符合精度要求的映射断面，其映射函数系数见表2。

图6 公路两车道隧洞映射断面Fig.6 Mapping cross-sections of two-lane road tunnel

表2 公路两车道隧洞映射函数系数Table 2 Coefficients of the mapping function for two-lane road tunnel

4.3 复杂断面隧洞

以文献[17]中的复杂断面隧洞为计算对象开展映射函数求解，隧洞几何尺寸及在坐标系中的位置如图 7所示。映射函数的项数构成为24项。图8为迭代1次、4次和9次后的映射断面，图形的绘制方法同算例1和2。其中，图8(c)为迭代9次后满足精度要求的映射图形，其映射函数系数Aj和Bj见表3。

图7 复杂断面隧洞Fig.7 Cavern with complex cross-section

图8 复杂隧洞映射断面Fig.8 Mapping cross-sections of complex cavern

表3 复杂断面隧洞映射函数系数Table 3 Coefficients of the mapping function for cavern with complex cross-section

5 讨 论

上述给出了单位圆外域到任意开挖断面隧洞外域共形映射计算方法的研究过程和3个开挖断面隧洞映射函数的求解算例，下面就其中几个关键问题与文献中的方法作比较和讨论：

（1）初始对应点的确定。初始对应点的数量取决于映射函数的构成项数，与其他因素无关。映射函数中的构成项数为m时，则映射函数中有2m个待求的实系数，此时需要m对对应点才能求之（即式（7）、（8）、），考虑到采用奇偶点相互迭代计算，故需对单位圆周线作2m等分。初始对应点的对应关系，文献[16]假设有1对对应点是已知的，文献[17]将任意开挖断面隧洞边界线作m等分，得到m条线段，再将这些线段的中点作为与单位圆周线上确定的计算点相对应。本文在隧洞复杂开挖断面边界线上确定与单位圆周线上相对应的m个计算点时，是可任意的。从这个意义上来说，本文的方法更简单且更具一般性。

（2）计算点（映射点）的调整。文献[16]通过映射洞形和实际洞形上任意相邻两点的距离之比与相应的映射洞形周长和真实洞形周长之比相等的原则，重新确定实际洞形边界线上的计算点；文献[17]通过映射洞形上所有相邻2点间距离和映射洞形总周长之比，重新确定实际洞形边界线上的计算点，将这些点代替原映射点进行下一轮的映射函数系数计算。本文给出的法线逼近法，将映射点和坐标原点相连，可获得该连线与实际洞形边界线的交点，用该点代替原计算点（映射点）作为下一轮的计算点。相比较而言，该方法计算量少且更简单有效。

（3）迭代计算的收敛条件。文献[17－19]未对映射精度作讨论。文献[16]以映射洞室上所有映射点与实际洞形边界线上相应点距离的算术平均值与真实洞室边界总长之比作为计算控制条件，并给出比值取 5%时就能达到工程计算精度需要的结论；以标准两车道公路隧道为例，其开挖边界线周长在30 m以上，以该结论作为计算的收敛条件，映射洞形边界线与洞室实际开挖边线的总体误差将在1.5 m以上。本文给出隧洞映射边界线与开挖边界线之间的最大距离小于某一数值ε作为迭代计算的收敛条件，确保了映射函数的映射精度。另外，本文虽然给出了相对误差的迭代计算收敛条件，但在实际使用时应慎用。该条件只说明了迭代的有效性，即前一次迭代计算和后一次迭代计算结果之间的差值如果很小，说明继续迭代下去已经没有意义，却不能保证映射洞形接近于真实洞形。

（4）映射函数的构成项数。映射函数中的项数构成和隧洞洞形的复杂程度及对映射函数的精度要求有关。一般来说，洞形越复杂，则需要更多的项数去近似和逼近精确的映射函数；对映射函数的精度要求越高，则项数越多越容易获得符合精度要求的映射函数。另外，映射函数的项数构成与迭代的次数也有关系，对同一洞形，在相同精度条件下项数越多，则迭代计算收敛越快。

6 结 论

（1）单位圆外域到任意开挖断面隧洞外域共形映射函数的求解难点在于两条边界线上点的对应关系难以确定。从三角插值理论出发，采用奇偶插值点相互迭代计算可有效地克服此困难，方便地求得映射函数系数。

（2）映射点（计算点）的调整方法决定了映射函数系数迭代计算的快速性和有效性，提出的法线逼近法可较方便地解决该问题；以绝对误差即隧洞映射边界线与开挖边界线之间的最大距离小于某一数值，作为迭代计算的收敛条件，保证了计算映射函数的映射精度。

（3）给出的3个开挖断面隧洞映射函数求解算例及对研究结果的讨论表明，本文给出的映射函数求解算法、映射点（对应点）的调整方法及迭代计算的收敛条件更具可操作性（易于编程计算），两条边界线上对应点的数量和映射函数构成项数之间的关系更为明晰，整个求解过程更为简洁快速。

[1]FAHIMIFAR A, TEHRANI F M, HEDAYAT A, et al.Analytical solution for the excavation of circular tunnels in a visco-elastic Burger’s material under hydrostatic stress field[J]. Tunnelling and Underground Space Technology, 2010, 25(4): 297－304.

[2]GAO G Y, CHEN Q S, ZHANG Q S, et al. Analytical elasto-plastic solution for stress and plastic zone of surrounding rock in cold region tunnels[J]. Cold Regions Science and Technology, 2012, 72(3): 50－57.

[3]LI S C, WANG M B. Elastic analysis of stress–displacement field for a lined circular tunnel at great depth due to ground loads and internal pressure[J].Tunnelling and Underground Space Technology, 2008,23(6): 609－617.

[4]VERRUIJT A. Deformations of an elastic half plane with a circular cavity[J]. International Journal of Solids and Structures, 1998, 35(21): 2795－2804.

[5]ZHANG J U, HOU D Q, CHEN X P. Displacement analytical solution of a deep elliptical tunnel in transversely isotropic rock mass[C]//Advanced Materials Research. [S. l.]: [s. n.], 2012, 402: 593－597.

[6]陈子荫. 围岩力学分析中的解析方法[M]. 北京: 煤炭工业出版社, 1994.

[7]祝江鸿. 公路小净距隧道围岩稳定性与合理净距研究[博士学位论文D]. 北京: 北京科技大学, 2012.

[8]祝江鸿. 隧洞围岩应力复变函数分析法中的解析函数求解[J]. 应用数学和力学, 2013, 34(4): 345－354.ZHU Jiang-hong. Analytic functions in stress analysis of the surrounding rock for caverns with the complex variable theory[J]. Applied Mathematics and Mechanics, 2013, 34(4): 345－354.

[9]吕爱钟, 张路青. 地下隧道力学分析中的复变函数方法[M]. 北京: 科学出版社, 2007.

[10]胡龙飞, 刘全坤, 冯秋红. 基于共形映射的壁板型材挤压模具型腔建模[J]. 机械工程学报, 2008, 44(8): 180－184.HU Long-fei, LIU Quan-kun, FENG Qiu-hong. Modeling for shape extrusion of aluminum alloy flat-plate based on conformal mapping[J]. Chinese Journal of Mechanical Engineering, 2008, 44(8): 180－184.

[11]李成, 郑艳萍, 李大磊. 积分方程法对含复杂孔形复合材料板孔边应力分布的研究[J]. 机械强度, 2006, 28(6):931－936.LI Cheng, ZHENG Yan-ping, LI Da-lei. Research on hole-edge stress distribution of composite materials plate with complex holes by integral equations method[J].Journal of Mechanical Strength, 2006, 28(6): 931－936.

[12]MUSKHELISHVILI N I. Some basic problems of the mathematical theory of elasticity[M]. Groningen: P.Noordhoff Ltd., 1953.

[13]闻国椿. 共形映射和边值问题[M]. 北京: 高等教育出版社, 1985.

[14]徐趁肖, 齐红元, 朱衡君. 任意边界域映射建模理论及模具设计应用[J]. 机械工程学报, 2002, 38(9): 83－86.XU Chen-xiao, QI Hong-yuan, ZHU Heng-jun. Modeling theory of complicated contact region conformal mapping and application on die design[J]. Chinese Journal of Mechanical Engineering, 2002, 38(9): 83－86.

[15]徐趁肖, 朱衡君, 齐红元. 复杂边界单连通域共形映射解析建模研究[J]. 工程数学学报, 2002, 19(4): 135－138.XU Chen-xiao, ZHU Heng-jun, QI Hong-yuan.Analytical modeling of complicated boundary simply-connected region conformal mapping[J]. Journal of Engineering Mathematics, 2002, 19(4): 135－138.

[16]皇甫鹏鹏, 伍法权, 郭松峰, 等. 基于边界点搜索的洞室外域映射函数求解法[J]. 岩石力学, 2011, 32(5):1418－1424.HUNAGFU Peng-peng, WU Fa-quan, GUO Song-feng,et al. A new method for calculating mapping function of external area of cavern with arbitrary shape based on searching points on boundary[J]. Rock and Soil Mechanics, 2011, 32(5): 1418－1424.

[17]朱大勇, 钱七虎, 周早生, 等. 复杂形状洞室映射函数的新解法[J]. 岩石力学与工程学报, 1999, 18(3): 279－282.ZHU Da-yong, QIAN Qi-hu, ZHOU Zao-sheng, et al.New method for calculating mapping function of opening with complex shape[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 18(3): 279－282.

[18]ZHU J H, JI H G. The comparison of surrounding rock stress analytical solutions for three tunnel sections[C]//Proceedings of the 2nd ISRM International Young Scholars Symposium on Rock Mechanics. Leiden: CRC Press, 2011: 741－744.

[19]王润富. 一种保角映射法及其微机实现[J]. 河海大学学报, 1991, 19(1): 86－90.WANG Run-fu. A method of conformal mapping and its computer implementation[J]. Journal of Hohai University, 1991, 19(1): 86－89.

[20]DELILLO T K, ELCRAT A R, PFALZGRAFF J A.Numerical conformal mapping methods based on Faber[J].Journal of Computational and Applied Mathematics,1997, 83(2): 205－236.

[21]范广勤, 汤澄波. 应用三个绝对收敛级数相乘法解非圆形洞室的外域映射函数[J]. 岩石力学与工程学报,1993, 12(3): 255－263.FAN Guang-qin, TANG Deng-bo. Determination of the mapping function for the exterior domain of a non-circular opening by means of the multiplication of three absolutely convergent series[J]. Chinese Journal of Rock Mechanics and Engineering, 1993, 12(3): 255－263.

[22]郑志强. 单位圆到任意曲线保角变换的近似计算方法[J]. 应用数学和力学, 1992, 13(5): 449－457.ZHENG Zhi-qiang. An approximate method on the conformal mapping from a unit circle to an arbitrary curve[J]. Applied Mathematics and Mechanics, 1992,13(5): 449－457.

[23]SCHINZINGER R, LAURA P A A. Conformal mapping:methods and applications[M]. New York: Dover Publications, 2003.

[24]MCLACHLAN N W. Complex variable theory and transform calculus with technical application[M].Cambridge: Cambridge University Press, 1963.