断层错动数值模拟中的几个关键影响因素

张炳焜

(四川省建筑科学研究院，四川成都610081)

针对断层错动类问题的有限元数值模拟，有几个关键问题值得注意，全面的考虑这些问题无疑能使计算结果更为准确。本项研究分别考虑了以下几个影响因素：准静态法、网格尺寸、应变软化和材料阻尼。为了能够说明这四个影响因素的重要性，本文建立简单均一上覆砂土层二维有限元模型，深度H=20 m，宽度L=80 m，土层内摩擦角φ=35°，剪胀角为ψ=15°，弹性模量E=60 MPa，泊松比ν=0.3，不考虑内聚力的影响。以各影响因素为变量展开计算，模拟逆断层60°，垂直位错量d=1 m时上覆土层应变与位移变形情况。

1 准静态分析

准静态法是一种用静力学方法近似解决动力学问题的简易方法，该方法能在有限程度上反映荷载的动力特性，它物理概念清晰，与全面考虑结构物动力相互作用的分析方法相比，计算方法较为简单，计算工作量很小、参数易于确定，并积累了丰富的使用经验，易于设计工程师所接受。在进行数值模拟时，有两个方面值得注意，一是在对地震导致的地面震动问题进行分析时，动力分析法是必不可少的。而另一方面，对于断层错动问题，若是综合考虑断层错动产生的波动和位移这两种因素的影响，则模拟结果会缺乏针对性。为了能够有针对性的说明断层错动时上覆土层的响应情况，选择准静态法进行模拟。ABAQUS/EXPLICIT模块常常在解决准静态类型问题是非常有效的。

本文对断层垂直错动速度分别为v=1 m/s，v=0.5 m/s，v=0.1 m/s和v=0.05 m/s的情况进行了计算和对比，从上覆土层塑性应变云图(图1)可以发现，当v=1 m/s时，塑性应变带未曾贯通上覆土层，而是较为集中的分布在了上覆土下部，由于高速的下部断层错动，导致上覆土层下部，即靠近断层处的土体受到高速的撞击而产生局部的应变，该区域的应变集中，消耗了本应传递至上覆土上部的应变能量。随着错动速度的降低，当v=0.1 m/s时，塑性应变带较为完整的贯通了上覆土层。根据以上分析得出结论：相对高速的断层错动速度会导致上覆土下部应变的集中，并且地面倾斜值低于相对低速的情况。而当断层错动速度减小到一定程度后，若是继续降低断层错动速度，覆土的塑性应变和地面位移将不再发生明显的变化，反而增加了计算所需的时间。因此，本文推荐采用的断层错动速度为0.1 m/s。

(a)v=1 m/s

(b)v=0.5 m/s

(c)v=0.1 m/s

(d)v=0.05 m/s图1 塑性应变云图

2 网格尺寸

对有限元模拟计算而言，除了材料参数的选取以及模型的建立外，模型的网格划分也是非常关键的，网格划分的好坏直接影响计算的精度。断层错动情况的网格尺寸是个关键因素，与其他模拟不同，断层错动伴随着较大的位移和应变的发生，若是模型网格尺寸较大，则可能导致塑性应变区的宽度发生变化，即破裂带可能变宽，而影响了计算精度，对后续的研究造成误导。本文有限元模型中部断层区域是研究的重点，在该区域内的网格尺寸更应该进行严格的控制，力求达到均匀和精细。

本文对模型中部网格尺寸分别为L=10 %H，L=7.5 %H，L=5 %H和L=2.5 %H的情况进行了计算和对比，H为覆土深度。从上覆土层塑性应变云图(图2)可以发现，随着网格尺寸的减小，塑性应变带的宽度逐渐减小，当L=2.5 %H时，塑性应变带宽度达到合理水平。因此，本文推荐网格尺寸L=2.5%H。

(a)L=10 %H

(b)L=7.5 %H

(c)L=5 %H

(d)L=2.5 %H图2 塑性应变云图

3 应变软化

应变软化是指材料试件经一次或多次加载和卸载后，进一步变形所需的应力比原来的要小，即应变后材料变软的现象。在断层发生错动后，涉及到较大的上覆土体变形和应变，土体强度的降低是应该考虑的，若是不考虑此因素的影响而进行数值模拟，得出的模拟结果与实际可能会产生一定的差别，之前学者们对此得出了较好的对比结果[1-3](Roth et al.1981,1982;Loukidis 1999)，针对路堤边坡稳定性而言，亦有类似结论[4-5](Potts et al.1990,1997)。针对摩尔库伦模型，内摩擦角、剪胀角和内聚力是该模型最主要的强度参数，将上覆土层的塑性剪应变e作为控制参数，考虑内摩擦角、剪胀角和内聚力随塑性剪应变的增加而数值逐渐降低，当塑性剪应变达到某一个临界值后，内摩擦角、剪胀角和内聚力将不再发生变化，保持其残余值。此处，内摩擦角和剪胀角的峰值分别记为φmax和ψmax，残余值分别记为φres和ψres，塑性剪应变临界值记为e0。则内摩擦角和剪胀角与塑性剪应变e之间的关系为：

(1)

(2)

本文分别对考虑应变软化和未考虑应变软化的情况进行数值模拟，当考虑应变软化时，内摩擦角和剪胀角的残余值分别为φres=25°，ψres=0°，塑性剪应变临界值e0=0.1。

从以上两种情况的塑性应变云图可以看出(图3)，二者的塑性应变带均贯通了上覆土层，并且分布较为均匀，但是，当未考虑应变软化时，破裂带和水平面的夹角明显较小。从某种意义上降低了断层的倾角。为了能够充分的估计断层错动对上覆土层的影响，考虑应变软化模型是合理的，有必要的。

(a)考虑应变软化

(b)未考虑应变软化图3 塑性应变云图

4 材料阻尼

对于准静态问题而言，虽然该方法大大的降低了动力波动对模型的影响，但是材料的阻尼仍然是不能忽视的问题，本文考虑Rayleigh阻尼的影响。

在有限元分析中，当采用振型分解法时，即把多自由度结构的复杂振动分解为按各个振型的独立振动的叠加，则常假定阻尼矩阵[C]满足振型的正交条件，即:

(3)

瑞利(Rayleigh)曾提出阻尼矩阵采用质量矩阵与刚度矩阵的线性组合，必满足正交条件，即：

[C]=α[M]+β[K]

(4)

阻尼矩阵用转置振型矩阵及振型矩阵前后相乘后，将得: 如令α=0，即C=β[K]，此时β=2ξn/ωn

ξn=(α/ωn+βωn)/2

(5)

如令β=0，即C=α[M]，此时α=2ξnωn

可见，当阻尼矩阵仅与质量矩阵有关时，阻尼比与振动频率成反比，高阶振型的阻尼就非常小；同样，当阻尼矩阵仅与刚度矩阵有关时，阻尼比正比于频率，高阶振型的阻尼非常大，这就使结构的高阶振型对结构动力反应的影响减弱甚至消除了。

系数α、β可由给定的第i、第j阶振型的阻尼比ξi、ξj反算得到：

(6)

通常结构振动分析中广泛采用的是简单的Rayleigh阻尼模型。

在工程中常采用简化求法，即：

α=ξ1·ω1β=ξ1/ω1

(7)

本文分别对考虑材料阻尼和未考虑材料阻尼的情况进行数值模拟，从塑性应变云图可以看出(图4)，当考虑材料阻尼时，塑性应变带贯通了上覆土层，并且分布较为均匀，而未考虑材料阻尼时，破裂带除了贯通上覆土层外，在破裂带左侧地面出现了一个局部的塑性应变区域，上覆土层破裂的不规律。因此，虽然采用准静态法进行数值模拟，但是材料的阻尼仍是不容忽视的关键因素之一。

(a)考虑材料阻尼

(b)未考虑材料阻尼图4 塑性应变云图

5 结论

通过以上计算分析，得出以下结论：对于地震引发的地面震动问题，动力分析方法是切实有效的。不过对于下覆断层错动类问题，为了能够较有针对性的对问题进行分析，建议采用准静态法进行研究，并综合考虑有限元网格尺寸、土体应变软化特性以及材料阻尼等影响因素。避免对上覆土层及结构的影响估计不足，使分析结果与实际贴切，为日后的工程设计和研究提供帮助。

[1] Roth WH, Scott RF, Austin I (1981). Centrifugemodelling of fault propagation through alluvial soils.Geophys Res Lett 8(6):561-564

[2] Roth WH, Sweet J, Goodman RE (1982) Numerical and physical modelling of flexural Slip phenomena and potential for fault movement. Rock Mech 12:27-46 (Suppl.)

[3] Loukidis D (1999). Active fault propagation through soil. Diploma Thesis, School of Civil Engineering, National Technical University, Athens

[4] Potts DM, Dounias GT, Vaughan PR (1990). Finite element analysis of progressive failure of Carsington Embankment. Géotechnique 40(1):79-101

[5] Potts DM, Kovacevic N, Vaughan PR (1997). Delayed collapse of cut slopes in stiff clay. Géotechnique 47(5):953-982