地下结构方法，基于ANSYS的工程实例

1.工程背景

本次优化选取的工程项目为一公路隧道。隧道长度为520 m，属于曲墙式连拱隧道。隧道处于平曲线中，曲线半径为1 900 m，隧道纵坡变坡点桩号为K148+550，其前后纵坡分别为1.966%和-1.800%，左右线的坡率一致。隧道超高为2%，进、出洞门均为1∶1削竹式洞门。该隧道所处地段的围岩主要为强风化砂岩，初步鉴定为Ⅳ级和Ⅴ级围岩，隧道实际埋深约为35 m左右，隧道走向与山体走向斜交，属于典型的偏压连拱隧道。

该连拱隧道典型断面如图5-2所示，图5-3为隧道施工现场照片。

图5-2　中龙隧道入口处典型断面

根据实际工程问题，抽象出需要进行优化计算的模型。按照实现性强和具有实际工程意义的原则，本次优化设计选取连拱隧道的衬砌结构进行优化计算。在实际数值计算中选取的断面如图5-4所示。

图5-3　中龙隧道施工现场照片

图5-4　偏压连拱隧道整体模型示意（单位：m）

该隧道的断面几何参数见图5-5所示。

图5-5　隧道断面几何尺寸（单位：m）

2.材料计算参数的选取

对于地下工程的数值计算而言，材料计算参数的选取是一个常被简化处理但又十分重要的问题，因为计算参数选择的合理与否直接关系到数值计算的成败。同济大学地下建筑与工程系老教授侯学渊先生就曾经说过，数值计算材料参数选取得不对，数值计算就只是“Garbage in，garbage out”。然而，对于实际工程项目，现场地质勘察报告和岩土勘察报告给出的岩土材料的物理力学参数一般都只是一个大致的范围，而数值计算要求的材料计算参数必须是一个确定值，如何合理地解决这个矛盾呢？

基于普氏理论与公路隧道设计规范推荐公式在我国隧道工程界的广泛应用，可以通过该公式所得到的松动圈范围来反演数值模型中的计算参数。下面以《公路隧道设计规范》为依据，借助有限元软件ANSYS，以围岩的极限拉应变为判据给出模型材料的力学性能参数。技术路线如图5-6所示。

按照公路隧道设计规范

对于本节依托的双连拱隧道，单洞跨度为12 m左右，由此分别计算出Ⅳ，Ⅴ级围岩的荷载等效高度为：

图5-6　研究技术路线

Ⅳ级围岩：

Ⅴ级围岩：

根据隧道规范，对于Ⅳ，Ⅴ级围岩，深浅埋的分界高度为hq=2.5h，由此计算出按照隧道规范，Ⅳ，Ⅴ级围岩深埋隧道的最小埋深分别为：15.3 m与30.6 m。实际工程中，由于地表风化层的存在以及岩体结构节理裂隙的影响，深浅埋分界标准不能完全绝对化，一般而言，覆土厚度在30～50 m为宜。

结合石吉中龙隧道实际工程，对于Ⅴ级围岩，计算选取断面为K140+200，隧道围岩为强风化变余砂岩，结构破碎，隧道埋深35 m左右；对于Ⅳ级围岩，计算选取断面为K140+235，隧道围岩为弱风化变余砂岩，结构较破碎，隧道埋深35 m左右。

由此建立起单洞隧道模型如图5-7所示。

一般而言，岩石在力学作用下发生张拉、剪切或拉剪破坏。对于岩土工程，可以用最大拉应变准则解决岩体破坏问题。谢兴华（2004），T.R.Stacey（1986），杨林德（2001、2003），曲海峰（2007），严宗雪（2009）等人在这方面均进行了一定的研究。借鉴上述研究成果，这里采用极限拉应变值确定隧道开挖后围岩松动圈范围。

岩石的极限拉应变[ε]可以通过岩石的单轴抗拉强度R拉与弹性模量E的比值确定：

图5-7　单洞隧道有限元数值模型

在数值模拟中，通过ANSYS路径功能，在隧道顶部布置深度为20 m的参考线，等分为100段，计算相邻两测点的径向平均拉应变值ε：

式中，δ为相邻两测点的相对位移，ΔL在本数值模拟中为0.2 m。计算判断松动区深度的表达式为：

在本节中，[ε]值的选取根据石吉线地质勘查报告（此处参照石吉线五峰山详细地质勘查报告，围岩情况与中龙隧道类似），并参考公路隧道设计规范围岩参数的选取范围[ε]=0.2‰。

图5-8　隧道顶部设置的参考线

通过反复调整围岩参数，使得计算得到的单洞隧道松动区范围与公路隧道规范中的松动区范围相符合，0°参考线方向（图5-8）极限拉应变与松动区深度的关系如图5-9、图5-10所示。

图5-9　Ⅳ级围岩中0°参考线方向拉应变与围岩松动区深度关系曲线

图5-10　Ⅴ级围岩中0°参考线方向拉应变与围岩松动区深度关系曲线

通过数值模拟得到的Ⅳ，Ⅴ级围岩松动区深度如表5-1所示。

表5-1　数值模拟得到的Ⅳ，Ⅴ级围岩松动区深度

根据以上围岩松动区深度绘制出的松动区范围如图5-11、图5-12所示（其中红线为《公路隧道设计规范》给出的松动区范围，黑线为数值模拟的结果，可见二者是十分接近的）。

图5-11　Ⅳ级围岩松动区范围

图5-12　Ⅴ级围岩松动区范围

最终确定的围岩参数如表5-2所示。

表5-2　石吉中龙隧道围岩参数表

应该指出，围岩力学参数不是围岩的真实或者绝对的力学参数，由于岩土材料的复杂性，即使对于同一级别围岩，该值亦发生变化。这里将围岩力学参数设为定值，只是为了合理地减少设计变量，便于衬砌结构的优化。优化计算时选取Ⅴ级围岩计算。

3.基于ANSYS的偏压连拱隧道衬砌优化

1）有限元数值模型的建立

按照平面应变处理的有限元模型主要包括两种材料类型，模拟隧道围岩的平面单元和模拟隧道衬砌的梁单元。各材料的相关参数由中反演分析得到。值得指出的是，为了施工开挖时模拟的方便，将相同的材料赋予不同的材料号，此时只是材料的某一项指标稍有不同。数值模型的计算范围根据已有的经验确定。具体是，左右边界距离隧道外边缘均为连拱隧道三倍的单洞洞径；下边界距离隧道底边也是三倍单洞洞径，上部边界取自由边界。已有的研究无论是数值模拟或是现场实测均证明隧道开挖对周边的影响不会超过三倍的洞径。

计算模型选择为平面应变模式，模型的边界条件处理为位移边界条件，具体表述为：左右两侧竖直边界约束其水平方向的位移（即u=0）；底边水平边界约束其竖直方向的位移（即v=0）；上边界视作自由边界，不受约束。整体模型及隧道局部如图5-13—图5-16所示。

图5-13　有限元计算整体数值模型

图5-14　隧道局部详图

图5-15　隧道开挖后局部详图

图5-16　连拱隧道衬砌结构（平面应变实体模式显示）

2）施工工法选择

目前隧道工程界对于连拱隧道有多种施工工法，较为常用的是中导洞上下台阶法施工，本次数值模拟即采用中导洞上下台阶法施工，模拟先开挖埋深较浅一侧（左侧）的方案（图5-17）。

3）选取设计变量及优化目标

一般而言，工程结构的极限状态包括两方面：正常使用极限状态和承载能力极限状态。正常使用极限状态对应结构或者构件达到正常使用或者耐久性能的某项规定限制；承载能力极限状态对应结构或构件达到最大承载能力或不适于继续承载的状态。

对于公路隧道而言，为满足正常使用所需的建筑限界应该对隧道拱顶沉降和拱底隆起进行控制，使变形不至于太大从而影响正常营运。为了保证结构的安全可靠，应该对衬砌结构的最大拉、压应力进行控制从而不至于出现强度破坏。

据此，选取设计变量为连拱隧道的左洞衬砌厚度t1，右洞衬砌厚度t3，以及中导洞的衬砌厚度t2。考虑到基本的构造要求，t1、t3和t2均需满足一定的范围限制，具体见式（5.1.5）。其中，单位为m。

图5-17　先开挖浅埋一侧的施工方案

对该偏压连拱隧道进行优化设计的目标是在保证结构安全可靠并正常发挥功能的前提下，尽最大可能减少衬砌的钢筋混凝土用量，达到降低工程费用，节约经济的目的。优化目标的函数表达式为：

式中，C1为连拱隧道左洞平均周长；C2为中导洞平均周长；C3为连拱隧道右洞平均周长，均为常数。(www.xing528.com)

4）选取状态变量

状态变量（State Variables）：状态变量表征设计变量在变化过程中引起的工程控制参数的变化。在本工程实例中，选取不同的状态变量对应不同的极限状态。正常使用极限状态采用变形控制，对连拱隧道左右洞的洞周收敛进行控制。监测断面示意如图5-18所示。

图5-18　隧道断面监测示意图

正常使用极限状态采用变形控制，表达式为

式中，δlb，δlt，δrb，δrt分别表示左右洞的拱顶沉降和拱底隆起；δll，δlr，δrl，δrr表示左右隧道的水平位移。

按照公路隧道设计规范，取洞周收敛为隧道单洞跨径的0.6%，即：

承载能力极限状态采用应力控制，中墙和衬砌结构均采用C25钢筋混凝土。按照偏压连拱隧道的工程经验，中隔墙为受力的薄弱环节，因此先对中隔墙的应力进行控制。控制方程如下：

此外，还应分别对左右隧道衬砌结构和中导洞衬砌的应力进行控制，方程为

应力监测断面布置示意图如图5-19、图5-20所示。

图5-19　中墙应力监测

图5-20　衬砌结构应力监测

5）利用ANSYS进行衬砌结构的优化

（1）参数化建立模型

用设计变量作为参数建立模型的工作是在PREP7中完成的。本例中，设计变量为左边、中隔墙和右边隧道的衬砌厚度T1、T2和T3。

（2）参数提取并建立优化中的参数

建立模型并进行求解后，进行参数的建立和提取。这些参数一般为设计变量，状态变量和目标函数。提取数据的操作用*GET命令：（Utility Menu＞Parameters＞Get Scalar Data）实现。通常用POST1来完成本步操作，特别是涉及到数据的存储，加减或其他操作。

提取梁截面应力的命令流：

提取mises应力的命令流：

提取节点位移的命令流：

计算收敛位移：

这些参数的提取，都是为后面优化设计的状态变量做准备。可以通过LGWRITE命令（Utility Menu＞File＞Write DB Log File）生成命令流文件。LGWRITE将数据库内部的命令流写到文件Jobname.LGW中。内部命令流包含了生成当前模型所用的所有命令。

（3）进入OPT，指定分析文件（OPT）

进入优化处理器：

指定分析文件：

（4）声明优化变量

该步中，指定哪些参数是设计变量，哪些参数是状态变量，哪个参数是目标函数。ANSYS中，允许有不超过60个设计变量和不超过100个状态变量，但只能有一个目标函数。主要操作如下：

对于设计变量和状态变量可以定义最大和最小值，而目标函数不需要给定范围。每一个变量都有一个容差值，这个容差值可以由用户输入，也可以选择由程序计算得出。

（5）选择优化工具或优化方法

ANSYS程序提供了一些优化工具和方法。缺省方法是单次循环。指定后续优化的工具和方法用下列命令：

优化方法是使单个函数（目标函数）在控制条件下达到最小值的传统化的方法。ANSYS中提供了一阶方法（First Order）、零阶方法（Sub Problem）、随机方法（Random Design）、阶乘方法（Factorial Tool）、梯度法（Gradient）、等步长搜索（Sweep Tool）以及用户自定义方法（User）。选择合理的方法，对于得到合理的优化计算结果至关重要。

（6）进行优化分析及查看结果

选定优化循环控制后，即可进行分析，进行优化计算，并查看优化计算结果。

6）ANSYS优化方法对比

（1）零阶方法

零阶方法之所以称为零阶方法是由于它只用到因变量而不用到它的偏导数。在零阶方法中有两个重要的概念：目标函数和状态变量的逼近方法，由有约束的优化问题转换为非约束的优化问题。

程序用曲线拟合来建立目标函数和设计变量之间的关系。优化处理器开始通过随机搜索建立状态变量和目标函数的逼近。由于是随机搜索，收敛的速度可能很慢。需要合理的设计初值以加速收敛。每次优化循环生成一个新的数据点，目标函数就完成一次更新。实际上是逼近被求解最小值而并非目标函数。

状态变量也是同样处理的。每个状态变量都生成一个逼近并在每次循环后更新。

用户可以控制优化近似的逼近曲线。可以指定线性拟合，平方拟合或平方差拟合。缺省情况下，用平方差拟合目标函数，用平方拟合状态变量。

状态变量和设计变量的数值范围约束了设计，优化问题就成为约束的优化问题。ANSYS程序将其转化为非约束问题，因为后者的最小化方法比前者更有效率。转换是通过对目标函数逼近加罚函数的方法计入所加约束的。

搜索非约束目标函数的逼近是在每次迭代中用Sequential Unconstrained Minimization Technique（SUMT） 实现的，该方法在进行搜索时补偿与设计变量和状态变量无关，可能会导致搜索的越界。

一般来说，零阶方法可以解决大多数的工程问题，在基本确定设计空间的时候，一般可以先选择零阶方法进行优化。虽然零阶方法的精度相对较低，但是基本都可以得到全部最优设计序列（图5-21）。

图5-21　零阶方法得到的衬砌厚度和衬砌体积收敛曲线

（2）一阶方法

同零阶方法一样，一阶方法通过对目标函数添加罚函数将问题转换为非约束的。但是，与零阶方法不同的是，一阶方法将真实的有限元结果最小化，而不是对逼近数值进行操作。

一阶方法使用因变量对设计变量的偏导数。在每次迭代中，梯度计算（用最大斜度法或共轭方向法）确定搜索方向，并用线搜索法对非约束问题进行最小化。因此，每次迭代都有一系列的子迭代（其中包括搜索方向和梯度计算）组成。这就使得一次优化迭代有多次分析循环（图5-22）。

图5-22　一阶方法的衬砌厚度和体积优化设计序列

一阶方法精度高，但是计算代价也很大。同时，一阶方法还有一些需要特别注意的：

①一阶方法可能在不合理的设计序列上收敛。这时可能是找到了一个局部最小值，或是不存在合理设计空间。如果出现这种情况，可以使用零阶方法，因其可以更好的研究整个设计空间。也可以先运行随机搜索确定合理设计空间（如果存在的话），然后以合理设计序列为起点重新运行一阶方法。

②一阶方法更容易获得局部最小值。这是因为一阶方法从设计空间的一个序列开始计算求解，如果起点很接近局部最小值的话，就会选择该最小值而找不到全局最小值。一般可以用零阶方法或随机搜索验证得到的是否为局部最小值。

（3）随机搜索法

程序完成指定次数的分析循环，并在每次循环中使用随机搜索变量值。可以用OPRAND命令指定最大迭代次数和最大合理设计数。随机搜索法往往作为零阶方法的先期处理，它也可以用来完成一些小的设计任务。

（4）等步长搜索法

等步长搜索法用于在设计空间内完成扫描分析。将生成n*NSPS个设计序列，n是设计变量的个数，NSPS是每个扫描中评估点的数目。对于每个设计变量，变量范围将划分为NSPS-1个相等的步长，进行NSPS次循环。问题的设计变量在每次循环中以步长递增，其他的设计变量保持其参考值不变。

（5）乘子计算法

本工具[OPTYPE，FACT]用二阶技术生成设计空间上极值点上的设计序列数值。（这个二阶技术在每个设计变量的两个极值点上取值。）可以用OPFACT命令（Main Menu＞Design Opt＞Method/Tool）指定是完成整体的还是部分子的评估。对于整体评估，程序进行2n次循环，n是设计变量的个数。1/2部分的评估进行2n/2次循环，依此类推。

（6）最优梯度法

最优梯度法计算设计空间中某一点的梯度。梯度结果用于研究目标函数或状态变量的敏感性。梯度法在开始迭代时收敛较快，但越接近最优点，步长越小，逼近函数极小的过程是“之”形的，可能在可行域内无法得到最优的值。

在本例中，无论选取初始值为多少，发现都难以得到较好的优化设计值，优化搜索只是在初始值附近循环，这可能与模型本身有关，也说明用该方法可能有很大的模型依赖性（图5-23）。

图5-23　随机搜索法衬砌厚度和体积优化设计序列

表5-3　三种优化设计方法对比分析

从表5-3三种优化设计方法对比分析可以看出，随机搜索方法一般也能给出合理的优化设计值，但是该值的合理性与计算次数有关。零阶方法能得出较好的优化序列，运用一阶方法进行计算时，二者值比较接近，说明在本例的，零阶方法得到比较合理的值，一阶方法得到的是全局最优的较高精度的设计序列。