人工边界转换方法

静-动力分析中人工边界转换方法的研究
摘要：通过将粘弹性动力人工边界应用于同时考虑静力效应和动力效应的工程算例，阐明了此类问题静-动力分析人工边界转换时保证模型为静力平衡体的必要性。

通过将粘弹性静-动力统一人工边界应用于半无限空间体有限元模型的静力分析中，验证了静力计算中的误差将使模型动力分析的稳态反应出现相近的误差。

在此基础上，系统阐述了适用于同时考虑静力效应和动力效应的工程问题的静-动力分析人工边界转换方法。

关键词：人工边界，静力分析，动力分析，边界转换
Abstract:Though the application of dynamic viscous-spring artificial boundary to an engineering case with a consideration of both static and dynamic effect, and the application of the unified viscous-spring boundary for static and dynamic analysis to static analysis of a finite modal of half space, the problems of the applications of viscous-spring artificial boundary to this kind of engineering calculation was pointed out, and its corresponding solving method was proposed. On the base, a systematic switching method of these artificial boundaries was specified.
Keywords: artificial boundary, static analysis, dynamic analysis, switching of boundaries
1 前言
人工边界从广义上可分为静力人工边界和动力人工边界。

静力人工边界由来已久，通常有固定边界、滚轴边界等。

动力人工边界经过几十年的研究发展，已形成具有全局人工边界和局部人工边界的两大类别，并应用于各自适应的工程计算中[1]。

动力人工边界发展到现在已有透射边界、粘性边界、粘弹性边界等几种类型。

1994年，Deeks 提出粘弹性人工边界[11]。

1998年，刘晶波等人发展了二维的黏弹性人工边界[3]，又于2005年将其发展为三维时域黏弹性人工边界[4]。

2006年，刘晶波等人再将二维黏弹性边界发展成一致粘弹性人工边界及其对应的粘弹性边界单元[5]，并于2007年推导了三维一致粘弹性人工边界及等效粘弹性边界单元[6]。

目前对静-动力分析的普遍做法是采用静力人工边界和动力人工边界分别对静力问题和动力问题进行计算，将计算结果进行叠加后得到完整的结果[1]。

但由于叠加原理仅在线弹性小变形范围内适用，原则上不能应用于涉及非线性或大变形问题的分析。

目前对涉及非线性或大变形问题的静-动力分析，常用的人工边界转换方法主要有以下几种：（1）静力分析和动力分析都采用滚轴边界或固定边界；（2）静力分析中采用滚轴边界或固定边界，动力分析采用粘弹性边界、透射边界、粘性边界等人工边界；（3）静力分析和动力分析都采用静-动力统一边界，如粘弹性静-动力统一人工边界。

对第（1）种方法，由于固定边界使波动全部反射，已有许多文献证明其具有放大振动效应的作用，目前已经使用得不多。

刘晶波等人基于黏弹性动力人工边界和半无限空间中静力问题的基本解，建立了对动力问题和静力问题均适用的三维黏弹性静-动力统一人工边界，从而上述第（3）种方法得以解决[1]。

然而，在使用人工边界对地下结构进行动力分析时，还存在一些问题。

如第（2）种方法，由于在静-动力分析的人工边界转换时的方法存在问题，致使产生错误的结果。

在第（3）种方法中，将粘弹性静-动力统一人工边界应用于地下结构的静力分析时，其解与准确值存在误差。

本文将就此两问题进行论证和分析，并阐述合理的地下结构静-动力分析人工边界转换方法。

2 静力和动力有限元分析原理
2.1 静力分析原理
地铁等地下工程初始应力场的确定须先计算未开挖状态下围岩的自重应力场，进而根据施工步
骤，采用释放荷载法，计算出衬砌结构和围岩的静应力场[13]。

许多地下结构的自重应力场模型可以假设为半无限空间体，根据经典围岩压力理论和弹性力学理论，半无限空间体中距地表面任一深度h 处的应力状态可定义为[7][8]：
h V γσ= (1) h V H λγλσσ== (2)
式中，σV 为竖向应力，γ为围岩重度，σH 为横向应力，λ为侧压力系数。

对浅层围岩，可假设其为各向同性介质，侧压力系数可用泊松比表示如下[7]：
μ
μλ-=
1 (3) 对有限元计算中所取的有限区域，可以据此确定有限域边界条件。

2.2静力分析边界条件
根据上述静力分析方法，在静-动力共同作用问题的计算中，在进行动力分析之前须先确定地下结构模型的静应力场。

地下结构开挖前可将大地假设为半无限空间体，其在重力作用下的静力计算，根据对称性，模型中任一处的水平位移0=H u 。

因此，计算该应力场时，有限区域模型两侧可用水平约束即法向约束，底部可用全约束或仅约束竖直方向，顶面即地面应为自由边界。

2.3 动力分析原理
在动力荷载作用下，有限元体系在t +∆t 时刻的运动平衡方程为：
t t t t t t t t F u K u C u
M ∆+∆+∆+∆+=++ (4) 式中M 为体系的总质量矩阵；C 为体系的总阻尼矩阵；K 为体系的总刚度矩阵； u
t t +∆为体系的节点加速度向量； u
t t +∆为体系的节点速度向量；u t t +∆为体系的节点位移向量；t t F ∆+为外荷载向量。

体系的总阻尼矩阵采用瑞利阻尼[9]
：
C M K =+αβ (5)
式中α、β为常数，可按两种不同的振动频率下测得的阻尼比ξ加以确定。

计算中常数α、β可由
α+βω2i =2ωi ξi 和α+βω2j =2ωj ξj 求得[9]。

则α和β可表示为：
j i j i j i j i i j ωωωωωωωξωξα))(()
(2-+-=
(6)
)
)(()
(2j i j i j j i i ωωωωωξωξβ-+-=
(7) ωi 、ξi 分别为振型向量φi 对应的自振圆频率和阻尼比。

根据振型分析结果可求得ωi 和ωj ，阻尼比在计算中取ξi =ξj =0.05。

2.4 粘弹性人工边界
粘弹性人工边界从用途上可分为动力人工边界和静-动力统一人工边界；从具体实现方法上可分为弹簧-阻尼器边界单元和一致粘弹性边界单元。

下面先介绍动力人工边界。

具体计算中模型边界材料参数由其相邻的围岩介质材料决定。

则当人工边界采用等效的弹簧和阻尼器物理元件来模拟时，其弹簧系数和阻尼系数的计算分别如下[3][4]：
法向边界：
)1(2μαα+==R E
R G K T
T
T ， s T c C ρ= (8) 切向边界：
)
1(2μαα+==R E
R G K N
N
N ， p N c C ρ= (9) 式中K T 、K N 分别为法向与切向的弹簧刚度；R 为波源至人工边界点的距离；c s 和c p 分别为S
波和P 波波速；E 和G 分别为介质弹性模量和剪切模量；μ为介质泊松比；ρ为介质质量密度；αT 与αN 分别为切向与法向粘弹性人工边界参数，具体取值情况见表1。

表1 粘弹性动力人工边界中参数α的取值
模型类型 方向 α 二维人工边界
平面内法向 2.0 平面内切向 1.5 出平面切向 0.5 三维人工边界
法向 4.0 切向
2.0
若人工边界采用一致粘弹性边界单元来模拟，其边界单元等效剪切模量、等效弹性模量和阻尼系数可分别用如下几式计算[5][6]：
)
1(2~
μαα+==R E
h
R G h
G T T (10) ~
~
~
~~
~
~
)
1()21)(1()1(2)1()21)(1(νννμανννα--++=--+=R E h R G h
E N N (11)
])1[()
1(2])
1[(~
N
p T s N
p
T
s
c c n nE R c c n nG
R
ααμρααρη+-+=
+
-=
(12)
其中等效阻尼系数取的是法向和切向的平均值。

式中h 为等效边界单元厚度，~
ν为等效泊松比，
n 为计算模型维数。

其余符号意义同上。

表2 粘弹性静-动力统一人工边界中参数α的取值
模型类型 人工边界位置 方向 α 二维人工边界
模型底面
法向 α*/3 切向 2/3 模型侧面 法向 α*/4 切向 1/8 三维人工边界
模型底面
法向 α* 切向 2 模型侧面
法向 α* 切向
1/2
考虑到实施的方便，实际计算中边界材料常采用各向同性材料，这时上述等效剪切模量和等效
弹性模量之间存在隐含关系式~
E =2（1+~ν）~G 。

考虑到普通有限元材料泊松比应限制在0～0.5范围内，则该等效泊松比可按如下取值[5]。

⎪⎩
⎪⎨⎧≥--=其它，，T N T N T N 02/)
1/(22
/~
ααααααν (13) 即计算等效泊松比时应先确定T N αα/大小范围，再确定其计算式。

静-动力统一人工边界是在上述动力人工边界的基础上对人工边界参数α进行调整，其具体取值见表2[1]。

其中参数α*具体计算公式：
]
)(1[])(1)[1(26
222*R
d
R d +++-=
μα (10)
式中μ为围岩泊松比，d 为位置坐标，R 为荷载作用点到人工边界点的距离。

对于底面人工边界，式中D 取荷载作用点至边界单元的水平距离；对于侧面人工边界，式中D 取荷载作用点至边界单元的垂直距离。

3 静-动力分析中的人工边界转换的几种方法
在引言中已经介绍了地下结构静-动力分析中常用的三种人工边界转换方法，下面采用算例对第二种和第三种方法存在的问题进行分析。

3.1算例一
某双线铁路隧道，其净空跨度10.2m ，为三心圆断面；衬砌初期支护15cm ，二次衬砌30cm ；隧道埋深45m ，Ⅳ类围岩。

采用ANSYS 大型有限元分析软件对该隧道进行建模和计算，模型区域竖向取100m ，横向取100m 。

隧道及围岩物理力学参数见表4。

静力分析采用的边界条件如图1所示，其中g 为重力加速度。

动力分析采用如图2如示动力模型，其最外一层即为粘弹性动力人工边界，厚度2米。

图中P (t )为某实测的列车振动荷载，其时程曲线见图3。

第（2）方法在静力分析中采用滚轴边界或固定边界（可称之为传统静力边界条件），在动力分析中采用粘弹性动力人工边界，在静力分析和动力分析中的荷载和边界条件见表3。

图1隧道静力分析边界条件 图2 动力分析边界条件
表3 计算步骤及条件表
表4 围岩物理力学参数表
图3 列车荷载时程曲线
由于隧道结构是对称的，整个模型有限域以及荷载也是对称的，监测点可以取自其对称轴及其一侧即可；取隧道衬砌的拱顶、拱腰、拱脚以及仰拱底等几个点监测其动力响应情况。

各监测点位移时程曲线见图4，第一主应力时程曲线见图5。

从图4和图5可以看出，各监测点的竖向位移和第一主应力均呈现首先大幅波动，然后逐渐趋于稳定的情形。

而从图2可以看出，列车振动荷载在整个时间段为振幅相近的周期荷载，且其大小不足以使结构应力和位移发生如此大的波动。

据此，可以推断模型在零时刻就受到一个或一组与自重相当的不平衡力，追溯前面的计算过程，可知这不平衡力就是静力计算结果中的约束反力。

这说明这种方法存在很大问题。

图4 算法1各监测点竖向位移时程曲线
（a）拱顶
（b）拱腰
（c）拱脚
（d）仰拱底
图5 算法1各监测点第一主应力时程曲线
3.2 方法三的算例
下面将粘弹性静-动力统一人工边界应用于方法二算例中所述的隧道模型动力分析，检验方法三的适用性。

隧道模型如图6所示，具体的计算条件见表5，其中动力荷载P(t)为与方法二的算例相同的列车振动荷载。

监测点的设置同该算例，计算结果见图7和图8。

表5 方法三计算步骤和条件表
图6 双线铁路隧道统一人工边界模型
图7 模型几何中心点横向应力和竖向应力图
图8 模型底部中点横向应力和竖向应力图
从图中可以看出，对双线铁路隧道施加粘弹性静-动力统一人工边界进行动力分析时，模型节点的第一主应力和竖向位移在振动初期的波动仍为整个过程的最大波动，但其波动幅值与振动荷载引起的稳态反应幅值相近或仅有几倍之差。

这说明该方法中从静力分析到动力分析的应力场和荷载等过渡得较好。

3.3 粘弹性静-动力统一人工边界在静力分析中的应用
方法三在静-动力分析中人工边界转换时静应力场和荷载等过渡得较好，但静力分析时沿须考虑粘弹性静-动力统一人工边界的精度，以下通过某静力分析算例检验其精度。

由于半无限空间体具有经典弹性力学解，即使在岩土力学中也有理论解，前面已有论述；因此，对半无限空间体的有限元模型进行静力分析，将结果与理论解进行对比，可以检验其计算方法的有效性。

该模型的理论解可以通过(1)式、(2)式和(3)式算得。

为了比较，进行了理论解、固定边界的数值解和粘弹性静-动力统一人工边界的数值解三种工况的计算，其计算模型分别如图9、图10和图11所示。

具体的计算步骤和计算条件见表6。

在模型对称轴上设置观测点，由地面向下每10m设置一个观测点，并将三种模型计算结果进行比较，比较结果见图12和图13。

图9 半无限空间体理论解模型（工况一）
图10 传统静力边界模型（工况二）
图11 粘弹性静-动力统一人工边界模型（工况三）
图12 竖向应力对比图
图13 横向应力对比图
从图中可以看出，该有限元模型在施加粘弹性静-动力统一人工边界时，其静力计算结果与理论解相比有较大误差，特别是横向应力结果误差较大；而其施加前面所述静力边界条件时，竖向应力和横向应力结果都与理论解符合得相对较好。

这说明粘弹性静-动力统一人工边界应用于受重力作用的半无限空间体的静力计算时，其误差偏大；此结果可推广到许多地下结构的静力分析。

4 方法二的改进
根据以上分析，上述三种静-动力分析中的人工边界转换方法，方法三虽然从静力分析到动力分析的应力场和荷载等过渡得较好，但其相应人工边界静力分析误差偏大，本文将方法二算例中的计算程序加以改进，将表3中动力分析的荷载条件再加上静力分析中的约束反力。

静力计算模型仍采用图1所示模型，动力计算采用图14所示模型。

监测点设置与前面相同，各监测点动力响应见图15和图16。

图14 方法二的改进动力分析模型
图15 改进后动力分析模型的监测点竖向位移响应
(a)拱顶第一主应力响应
(b)拱腰第一主应力响应
(c)拱脚第一主应力响应
(d)仰拱底第一主应力响应
图16 改进后动力分析模型的应力响应
从图中可以看到，各监测点的最大竖向位移与该方法改进前的结果相比大为减小，其后其的稳态振动与初期相对较大的振动差值比前面也小许多，第一主应力也有类似的现象。

初期出现相对较大的波动是因为这里的列车振动荷载是作为突加荷载施加于模型的，上述动力模型改进后其初期动力响应波动值与稳态反应值接近，是人工边界转换更为合理的表现。

5 几种方法的比较
从以上算例中可以看到，几种方法下各监测点在振动初期的振幅都是最大的，这是因为这些计算中列车振动荷载都是在突加荷载的形式施加于模型的。

但各种方法下振动幅大小有所不同，表7列出了各种方法下监测点的第一主应力振幅大小。

从表中可以看到，使用方法二时多数监测点应力的振幅与使用另外两种方法时的结果比较为最大，其主要原因是方法二在静-动力分析的人工边界转换中未能较好地将静力状态施加于模型。

因此，方法二存在缺陷。

表7 各种方法下的第一主应力振幅
表8 各种方法下第一主应力的振动稳定值与静力解
从方法三的算例以及粘弹性静-动力统一人工边界在静力分析中的应用可以看到，方法三虽然从静力分析到动力分析的应力场和荷载等过渡得较好，但其相应人工边界应用于静力分析时误差偏大。

表8为各种方法下各监测点下第一主应力的振动稳定值与静力解的对比，此处的振动稳定值为前面各监测点第一主应力响应图中呈水平线的时间段的应力值，取最后一段稳定值。

从中可看到方法三算出的应力稳定值与静力解差值为几种方法中最大的，这再次证实了上述结论。

从方法二改进后的算例中可以看到，该方法改进后从静力分析到动力分析的应力场和荷载等过渡更加合理。

而表8显示改进后的方法其应力振动稳定值与静力解相对较为接近。

因此，本文推荐改进后的方法为地下结构静-动力分析方法。

5 结论
本文通过将粘弹性动力人工边界应用于同时考虑静力效应和动力效应的工程问题，阐明了此类问题静-动力分析人工边界转换时保证模型为静力平衡体的必要性。

通过将粘弹性静-动力统一人工边界应用于半无限空间体有限元模型的静力分析中，验证了该边界条件及其相应计算方法在地下结构静-动力分析中的局限性。

通过对几种方法下对同一隧道结构模型的动力分析结果的比较，本文推荐了改进的方法为最适用于地下结构静-动力分析的方法。

在此基础上，本文提出适用于地下结构静-动力分析的人工边界转换方法，其具体步骤如下：
(1) 采用静力方法和静力约束条件对模型进行静力分析，计算结果应包括应力场和约束反力。

(2) 将静力模型中的约束条件去掉，施加粘弹性人工边界，在边界单元外边缘施加约束。

(3) 在原静力约束被去掉的地方施加约束反力，对模型施加静力场和静力荷载（如重力加速度或其他恒力）。

(4) 对模型施加动力荷载，进行动力计算。

该方法为此类工程问题在静力分析阶段和动力分析阶段都使用精度较高的边界条件提供了条件，从而使静力人工边界和动力人工边界可以相对独立地发展，加速人工边界条件的发展。

参考文献
[1]刘晶波, 李彬. 三维黏弹性静-动力统一人工边界[J]. 中国科学(E辑). 2005, 35(9): 966—980
[2]高峰. 静应力场对隧道列车振动响应的影响分析[J]. 岩土工程学报. 2009, 31(7): 1105—1109
[3]刘晶波, 吕彦东. 结构-地基动力相互作用问题分析的一种直接方法[J]. 土木工程学报. 1998,
31(3): 55—64
[4]刘晶波, 王振宇, 杜修力, 杜义欣. 波动问题中的三维时域粘弹性人工边界[J]. 工程力学. 2005:
22(6): 46—51
[5]刘晶波, 谷音, 杜义欣. 一致粘弹性人工边界及粘弹性边界单元[J]. 岩土工程学报. 2006, 28(9):
1070—1075
[6]谷音, 刘晶波, 杜义欣. 三维一致粘弹性人工边界及等效粘弹性边界单元[J]. 工程力学. 2007,
24(12): 31—37
[7]关宝树. 隧道力学概论[M]. 成都：西南交通大学出版社. 1993
[8]徐芝纶. 弹性力学(上册)[M]. 北京: 高等教育出版社. 2006
[9]潘昌实. 隧道力学数值方法[M]. 铁道出版社. 1995
[10]杜修力. 有缝拱坝-地基系统非线性地震波动反应分析方法[J]. 地震工程与振动. 2000, (1)
[11]Deeks A J. Randolph M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of
Engineering Mechanics. 1994, 120 (1) : 25 - 42.
[12]蒋树屏, 蒋华, 王晓雯, 陈杨. 高烈度地震区公路隧道洞口段地震响应分析[J]. 公路. 2009,
(2):194—198
[13]高峰、关宝树. 深圳地铁地震反应分析[J]. 西南交通大学学报.2001，36(4) 355～359。