盾构隧道横断面反应位移法在ABAQUS软件中的开发及工程应用

2022-11-15梁建文于明东巴振宁杨贵生

自然灾害学报 2022年5期

梁建文，于明东，巴振宁，杨贵生

（1.天津大学建筑工程学院，天津 300350；2.中国铁路设计集团有限公司，天津 300142）

引言

反应位移法是一种简化的地下结构地震反应的计算方法，由日本学者于20世纪70年代提出［1-4］。近年来随着地下结构抗震设计需要，国内外学者对反应位移法进行了大量研究和探讨，完善了反应位移法的理论框架和应用细节［5-9］。目前，地下结构反应位移法已写入我国《城市轨道交通结构抗震设计规范》［10］和《地下结构抗震设计标准》［11］等多部国家规范。

关于结构横断面抗震分析的反应位移法，现行规范中针对不同场地类型提供了不同的应用模型。对于水平层状场地，需要求解一维地层地震响应，并将计算得到的地震响应信息赋予结构及地基弹簧远端，过程略显繁琐。在处理预制装配结构时，还需添加部分非线性弹簧单元以模拟接头的拉压异性、剪切及转动等行为，该过程难以在ABAQUS等通用有限元软件的可视化界面中操作。目前，尚未有专门针对地下结构抗震分析的可视化软件，工程技术人员需自行使用有限元软件建模求解，无疑增加了反应位移法的实际应用难度。

针对此问题，文中以盾构隧道横断面反应位移法为例，基于“梁-弹簧”模型模拟隧道结构及管片间接头在地震作用下的力学行为［12］，并借助有限元软件ABAQUS提供的二次开发功能，将一维地层地震响应［13］与盾构隧道有限元模型直接关联，实现了自动求解场地地震响应，并依据响应结果完善结构有限元模型、提交求解、输出结构内力图等环节的一体化及全自动化。同时，文中利用ABAQUS软件的GUI开发工具，设计出一套可视化软件，用户无需修改源代码，只需根据提示在可视化界面中输入场地、地震波、隧道结构等信息即可完成反应位移法计算，并直接输出隧道内力图。

文中软件适用于任意断面盾构隧道的横断面抗震分析，具有可视化界面输入、操作简便快捷、自由扩展模型前后处理的灵活性等优点，便于工程技术人员使用。此外，文中关于Python语言与Fortran语言在ABAQUS环境下联合开发、ABAQUS中非线性弹簧自动定义、ABAQUS直接输出隧道结构内力图等问题的讨论，对ABAQUS二次开发应用具有一定参考价值。

1 隧道横断面反应位移法

现行《城市轨道交通结构抗震设计规范》［10］规定，施加在结构上的地震作用包括地层相对位移、结构惯性力和结构与周围土层剪力，结构宜采用梁单元建模，同时通过在结构周边布置地基弹簧的方式模拟结构与周围土体的相互作用，并将土层位移施加在地基弹簧非结构连接端。隧道横断面反应位移法示意图如图1所示。

1.1 地基弹簧刚度系数

反应位移法通过地基弹簧传递土层位移的等效荷载，地基弹簧刚度系数计算公式如式（1）所示：

式中：k表示地基弹簧刚度系数；K表示基床系数；A表示单个地基弹簧反映的土体面积。

目前，基床系数的计算主要有3种方法：经验公式、静力有限元方法和规范提供的依据土体信息查表取经验值［10］。考虑到程序通用性以及兼顾求解精度，文中软件采用经验公式确定基床系数。

《城市轨道交通结构抗震设计规范》［10］中提供的基床系数计算公式如式（2）～式（3）所示：

图1 隧道横断面反应位移法示意图Fig.1 Schematic diagram of cross section response deformation method of tunnel

式中：Kn表示土层法向基床系数；Kt表示土层切向基床系数；G表示隧道处土体最大剪应变对应的剪切模量，R表示隧道半径。

1.2 土层相对位移

土层相对位移可由一维地层地震响应分析计算得到，文中通过编写Fortran程序完成场地响应计算。通过结构顶、底面处最大相对位移确定最不利时刻后，提取土层最大相对位移，并将其施加于相应地基弹簧非结构连接端。

1.3 结构与周围土层剪力

结构周围土体剪应力τA可由一维地层地震响应分析计算得到，圆形结构与周围土层剪力作用可转化为节点等效集中力，文中软件采用《城市轨道交通结构抗震设计规范》［10］提供的如式（4）～（5）所示：

式中：FAX表示作用于A点水平方向的节点力，FAY表示作用于A点竖直方向的节点力；τA表示A点处的剪应力；θ表示土与结构的界面A点处的法向与水平向的夹角。

文中采用以上公式计算等效节点力，并分解为节点径向力、切向力施加。

1.4 结构惯性力

反应位移法中结构惯性力可作为集中力作用于结构单元的形心。文中软件依据隧道结构形心处加速度计算结构惯性力，并以均布荷载形式施加于模型单元。

施加边界条件及等效地震荷载后的盾构隧道有限元模型如图2所示。

图2 盾构隧道有限元计算模型Fig.2 Finite element model of shield tunnel

2 ABAQUS GUI二次开发思路

2.1 基于Python语言的二次开发

ABAQUS二次开发分为子程序开发和图形用户界面程序（Graphical User Interface，简称GUI）2类：子程序开发基于Fortran语言，用户可根据需求编写材料本构关系、自定义单元等子程序；图形用户界面程序开发基于Python语言，用户可根据需求对原有ABAQUS/CAE界面程序进行扩展，开发专用的前后处理模块以及GUI工具［14］。文中程序为基于Python语言的用户图形界面程序。

ABAQUS软件在继承Python原有库函数的基础上，进一步扩展了专属的函数模块，通过调用这些函数，可直接调用ABAQUS/CAE内核程序，实现快速前后处理功能。Python脚本与ABAQUS/CAE内核程序进行交互主要有3种方式［15-16］：主窗口下部命令交互界面（CLI）、执行脚本程序（Script）以及图形用户界面（GUI）。具体交互命令通信模式如图3所示。

文中软件采用图形用户界面（GUI）方式，利用ABAQUS提供的RSG对话框构造器创建图形界面，设计后的整体界面下文介绍。

2.2 Python语言与Fortran语言在ABAQUS环境下的联合应用

图3 Python脚本与ABAQUS内核程序交互命令流Fig.3 The interactive command flow of Python script and ABAQUS kernel program

文中软件场地响应分析需由Fortran语言程序完成，而ABAQUS脚本基于Python语言，因此必须解决二者联合开发的问题。Python语言提供了大量的开源库，可实现与多语言混合编程，其中f2py.exe工具可将Fortran语言函数代码编译为Python语言的可扩展文件，从而实现二者的联合开发［17］。

文中软件Fortran语言的源代码文件为ARCS.f90，其中主函数为sub（…），选用GFortran编译器，于是编译指令为：

＞f2py-m ARCS-c ARCS.f90--fcompiler=gfortran--compiler=mingw32。

编译成功后获得可扩展文件ARCS.pyd，并将该文件存放至ABAQUS当前工作目录，于是在Python程序中调用该模块的具体方式为：

import ARCS#导入模块；

ARCS.sub（…）#调用主函数。

需要指出的是，ABAQUS二次开发接口是基于Python语言定制的，虽然保留了Python语言的大多数开源库，但并未提供f2py.exe工具以及其它Fortran语言混合编程所需的动态链接库文件（DLL）。因此使用者必须根据自己的软件版本、计算机系统配置以及所选编译器种类等补充相应的动态链接库文件或Python可扩展文件（PYD），才能在ABAQUS环境下使用Python语言正确调用Fortran语言程序。

3 计算系统的构建

图4 本文ARCS软件工作流程图Fig.4 The work flow chart of ARCS

基于上述反应位移法理论和软件开发思想，基于有限元软件ABAQUS二次开发，文中开发可以自动进行场地响应分析、建立盾构隧道横断面反应位移法模型、求解后输出隧道结构内力的高效自动化计算软件ARCS（ABAQUS-Response deformation method of Crosssection System，简称ARCS）。该计算软件可以实现场地响应分析与反应位移法模型自动关联，同时为用户提供可视化图形用户界面，用户只需在ABAQUS中建立初始状态的盾构隧道横断面梁单元模型、可视化界面中输入场地及模型信息，该计算软件自动识别模型信息、求解场地响应并完善模型、提交计算，最终输出内力结果及内力图。软件工作流程如图4所示。

首先是准备工作和信息输入：

（1）用户完成初步场地划分并在相应界面输入场地信息：在图5所示界面中输入地震波加速度时程并指定细分后的土层厚度；在图6所示界面中输入每一土层的厚度、土体类别编号、密度、剪切模型、泊松比等信息；在图7所示界面中按编号顺序输入土体剪切模量和阻尼非线性参数。

（2）用户完成隧道有限元模型初步工作，包括定义parts、定义material和section、定义step和划分网格单元。需要说明的是，本文软件允许使用者根据需求自定义网格划分情况，软件可自动识别划分后的单元信息。

（3）在图8所示界面输入已初步建立的有限元模型信息，以及隧道管片接头等效弹簧刚度、隧道埋深及重量等信息。其中，接头弹簧刚度系数的计算公式可参考相关文献资料［12，18］。

图5 本文ARCS软件图形用户界面-地震动信息输入界面Fig.5 The first page of ARCS GUI-earthquake information

图6 本文ARCS软件图形用户界面-地层参数信息输入界面Fig.6 The second page of ARCS GUI-site parameters information

图7 本文ARCS软件图形用户界面-土体非线性参数输入界面Fig.7 The third page of ARCS GUI-soil nonlinear parameters

图8 本文ARCS软件图形用户界面-有限元模型参数输入界面Fig.8 The fourth page of ARCS GUI-FEM parameters

至此，信息输入完毕，下面开始进行数据处理流程：

（4）调用ARCS.pyd模块完成场地响应分析，输出土体迭代模量和场地节点处的响应结果。

（5）根据隧道埋深，线性插值获得隧道节点处响应，以及结构顶、底处位移时程，并将土体剪切模量代入经验公式计算地基弹簧刚度系数。

（6）根据“结构顶、底面最大相对位移”准则确定隧道地震响应最不利时刻。已知结构顶、底处位移时程分别为u（t）、uB（t），结构处最大相对位移时刻tc可由式（6）确定：

（7）根据最不利时刻确定土层相对位移值、结构周围土层剪应力值以及隧道轴线处加速度值，并将结构周围剪力转化为模型节点集中力、将结构惯性力转化为模型单元均布荷载。

（8）为每一个节点新建一个局部坐标系、参考点，为模型节点和单元设置荷载边界条件、为参考点设置位移边界条件。文中软件在每一节点处定义一个局部坐标系，并在外法线方向创建一个参考点，为每一对“参考点-节点”和管片接头处“节点-节点”均定义一个Spring2类型单元，将上一步确定的土层相对位移以位移边界约束的形式施加于参考点水平方向，将剪应力转化后的节点力以集中荷载的形式施加于节点法向和切向，最后将惯性力均布荷载以线荷载的形式施加于所有单元的水平方向。

（9）根据已知的刚度系数为模型设置弹簧，包括管片间的接头弹簧和节点与参考点间的地基弹簧，因涉及在ABAQUS软件中设置非线性弹簧，详细流程下面介绍。在横断面模型中，隧道纵向取单位长度，文中软件将自动根据断面形状设置弹簧方向，具体表现为：关于管片间的接头弹簧，将依据接头处切线方向进行布置；关于地基弹簧，将依据节点处外法线方向进行布置。

（10）将完整的inp文件提交求解。因横断面反应位移法模型为拟静力问题，结构采用线性梁单元建模，因此采用常规ABAQUS/Standard求解器即可。

至此，模型求解结束，下面进入后处理阶段：

（11）等待计算完成后，读取输出数据库文件（ODB）并绘制隧道结构内力图，详细流程下面介绍。

（12）运行结束后，默认显示输出数据库文件（ODB），用户可自由进行模型后处理。

文中系统为实现程序运行连续，创新性地解决了2个问题：ABAQUS非线性弹簧自动定义问题和ABAQUS直接输出隧道结构内力图问题。下面就上述两点分别详细介绍。

3.1 ABAQUS软件中设置非线性弹簧

基于“梁-弹簧”模型的盾构隧道反应位移法计算模型中，管片接头等效弹簧包括与轴力相对应的拉压弹簧、与剪力相对应的剪切弹簧、与弯矩相对应的抗弯弹簧等3类，地基弹簧包括压缩弹簧、剪切弹簧2类。其中，管片接头拉压弹簧因衬砌混凝土抗压能力与螺栓抗拉能力差异，导致抗拉强度和抗压强度不一致，因此必须设置为非线性弹簧。由于土体只能受压不能受拉，只具有抗压刚度，同时，地基剪切弹簧具有明显的屈服位移，具体的屈服位移取值和土体性质有关，因此也需要设置为非线性弹簧。上述3类非线性弹簧的力学性质如图9所示。

图9 三类非线性弹簧力-位移曲线Fig.9 The force-displacement curve of three types of nonlinear springs

其中，管片接头弹簧抗拉刚度由单位宽度内的环向螺栓决定，抗压刚度由管片混凝土抗压强度决定，具体公式可见文献［12］；地基弹簧抗压、抗剪刚度可由文中式（1）～式（3）计算，剪切屈服位移的取值参考文献［18］。

ABAQUS/CAE提供了3种弹簧单元的定义方法，但都仅限于线性弹簧［19］。对于非线性弹簧，要求用户先完成模型其他部分并生成inp文件，然后在inp文件中添加设置非线性弹簧的关键字命令流，最后以inp文件创建作业提交计算。针对该问题，笔者提出一种可以自动设置非线性弹簧的方式，具体流程为：

（1）待完成模型除弹簧的其它内容后，生成xxx_0.inp文件，其中xxx为期望Job名称。

（2）逐行读取xxx_0.inp，找到“*End Assembly”所在行位置，将该行以上内容保存并记为A，将该行及该行以下内容保存并记为B。

（3）根据已知的弹簧信息，以字符串形式按序、逐行生成所有定义弹簧单元的关键字命令流，保存并记为C。对于每一个非线性弹簧，需先定义连接单元，再结合对应的局部坐标系定义弹簧力-位移曲线，并在命令流中添加“NONLINEAR”关键字。

（4）新建文档xxx.inp，依次将A、C、B写入。

（5）根据xxx.inp创建Job，并提交计算。

文中系统采用上述方法设置弹簧单元，可以实现系统的自动化和整体性，提高效率。

3.2 ABAQUS软件直接输出隧道结构内力图

ABAQUS软件求解功能强大，但其后处理功能，特别是数据结果展示功能不足，这也导致很多用户选择先利用ABAQUS求解，输出结果数据后再利用其它软件进行数据分析和展示。文中系统根据设计需求，提出一种可以直接在ABAQUS软件平台输出结构内力图的方法，具体流程为：

（1）Odb文件中读取所有单元场变量SF、SM，单元若有多个积分点则取积分点数值均值。

（2）准备一个以坐标原点为中心的结构形状曲线数据序列，保存为XY Data数据，记为数据a。

（3）将单元场变量数值按单元位置逆时针排序，根据已知节点外法线方向，并结合曲线序列a进行规范化，同样保存为XY Data数据，其中SF1记为数据b，SF3记为c，SM2记为d。

（4）新建3个窗口，每个窗口都绘制数据a，同时分别绘制数据b、c、d，分别为轴力图、剪力图、弯矩图，并在图中注释内力最值。

（5）将视图窗口保存为jpg图片。

文中软件采用上述方法由ABAQUS直接输出结构内力图，方便用户使用，特别是对ABAQUS后处理功能不熟悉的用户。

4 应用实例

以某地铁盾构隧道为例，选取典型横断面进行大震工况反应位移法分析，具体工程信息见文献［18］。

4.1 软件使用流程

ARCS计算软件使用流程如下：

（1）完成场地信息整理、ABAQUS有限元模型初步建立工作，并将插件程序文件复制到当前工作目录。其中，单元尺寸参考文献［18］进行处理，取圆形断面的1/60。

（2）打开模型数据库文件（CAE），找到主视图“Plug-ins”菜单下的“My-EERA-FORTRAN”选项，点击进入场地地震响应分析程序界面，输入场地信息及地震波信息，信息输入界面如图5～图7所示。

（3）点击“OK”或“Apply”按钮，进行场地地震响应分析计算，等待计算完成。

图10 大震作用下本文ARCS软件输出的隧道内力图（内力正负标准：正值在内，负值在外）Fig.10 Internal force chart of the tunnel output by the ARCS under the action of a large earthquake（The standard of internal force positive or negative：positive is inside，negative is outside）

（4）关闭场地地震响应分析程序界面，找到“Plug-ins”菜单下的“RDM-C”选项，点击进入盾构隧道反应位移法程序界面，输入隧道及有限元模型信息，输入后如图8所示。

（5）点击“OK”按钮进行反应位移法建模求解，等待求解完成后输出隧道结构内力图，如图10所示。

（6）关闭盾构隧道反应位移法程序界面，默认显示SM云图，用户可根据需要自由操作输出数据库文件（ODB）。

完成以上6步操作，即可实现盾构隧道横断面反应位移法模型从建立、到求解、最后输出内力响应结果等工作，简便高效。