主要参考3DEC5.2 manual → OPTIONAL FEATURES → Section 2 : Dynamic Analysis 及《块体离散元数值模拟技术及工程应用》一书
1. 动力分析的主要问题
动力分析主要针对岩土工程中以下类型的问题:
(a) 地震荷载、
(b) 爆炸或冲击荷载
(c) 采矿引起的能量释放(例如岩爆)
(d) 重力作用下的颗粒(圆形或不规则形状)流动
这些问题与时间密切相关,求解持续时间可能只有1分钟,甚至只有1秒钟(爆破冲击等),在这些过程中需要考虑阻尼、网格、时间步的特殊设置,并关注块体的运动过程、应力波的波动效应等。
2. 3DEC动力分析阻尼的设置
由于块体在运动时动能转变为热能耗散掉,运动时并不发生弹跳,因此岩体的运动是不可逆过程。
为了避免弹性系统由于具有动能后在平衡位置振荡,就要采用加阻尼的办法来耗散系统在振动过程中的动能。
在任何动力分析程序中,必需说明物理系统中的能量损耗,而这在数值算法中并没有考虑。通常,在高度弹性系统中采用较小的阻尼,而土层等岩土介质中则取较大阻尼。
(1)Rayleigh阻尼(Rayleigh Damping)
用 UDEC/3DEC 解决动力问题时,阻尼常采用 Rayleigh阻尼。Rayleigh阻尼最初应用于结构和弹性体的动力计算中,以减弱系统的自然振动模式的振幅。在计算时,假设动力方程中的阻尼矩阵C 与刚度矩阵M 和质量矩阵K 有关。
对于线性系统而言,Rayleigh阻尼可写为:
式中,[C]为黏性阻尼矩阵,[M]为质量矩阵为刚度矩阵,α 和 β 分别为质量阻尼比例系数与刚度阻尼比例系数,由下式决定:
式中:ξmin为临界阻尼比,ωmin为圆频率,其意义如下图所示:
基频定义为:
对于多自由度系统,可以从系统的任何角频率ωi找到临界阻尼比ξi:
或者:
可以看出,对于低频率振动的系统,质量阻尼比比较有效,而刚度阻尼则适用于快速高频振动的系统。
由于非连续系统随着块体之间的滑移或分离,其振型是不确定的,因此Rayleigh阻尼对于非连续介质系统不完全适用。
设置瑞利阻尼应该尽量选择中间频率ωmin。对于地质体,阻尼一般是独立于频率的,中间频率ωmin选择出现在数值模拟中频率范围的中间值(自振频率或者输入频率的主频),这可以根据谱分析获得。
如果在 3DEC 中使用Rayleigh阻尼来研究动力问题,调用格式如下: Damp freq fcrit
其中 freq
对应基频 fmin , fcrit
对应ξmin, 这里不是直接输入α和β的值。 但是Rayleigh阻尼参数的给出是一个相当麻烦的问题,通常动力分析中的阻尼假定为频率不相关,而 Rayleigh阻尼属于频率相关模型,在频谱曲线中有一个水平区(如下图), 跨度约为频域的1/3,该图可以基于速度时程的傅立叶变换分析获得。
如果上限优势频率值大于下限优势频率值三倍,那么即可认为存在一个 3:1 (频域范围/优势频段)跨度,该范围可涵盖能量谱中的绝大多数动能。此时在进行动力分析时可调整Rayleigh阻尼中ωmin, 使其如上图中的优势频率平直段一致,ξmin与准确的物理阻尼比一致。这里优势频率既不是输入频率,也不是系统的自振模态,而是二者的综合结果。其目的是在问题求解时对重要频率设置正确的阻尼值,而忽略其他频率影响。
(2)局部阻尼(Local Damping)
Rayleigh阻尼属于频率相关模型,其参数的确定十分麻烦。通常来说动力分析中的阻尼假定为频率不相关,我原来求解边坡的地震响应问题时使用的是局部阻尼(原来参考的一篇文献是这么做的,文献在下面)。使用局部阻尼的命令是 damp local value
,默认值为0.8,该值可以根据下面的进行计算。局部阻尼最初用于静力学分析,但是它的一些特性可以使其在动力学分析中更具优势。它是通过在系统的振动循环中增加或减少节点的质量来实现,但是又在整个过程中保持质量守恒。当速度改变符号时使质量增加,当速度超过峰值点(最大值或最小值)后质量减少。局部阻尼可通过下式求解:
ξ 为临界阻尼比,αL为局部阻尼系数。对于岩土体材料而言,其临界阻尼比一般在2%~5%之间
【参考文献】肖克强. 地震荷载作用下顺层岩体边坡变形特征及稳定性研究[D]. 中国科学院研究生院(武汉岩土力学研究所), 2006.
3. 3DEC动力分析网格尺寸要求
在动态分析中,波在传播过程中可能会失真,而输入波形的频率成分和岩体介质的波速特性都会影响波传播的数值精度,由于频率与波长是相关的,所以波的传播精度也与输入波形的波长有关。
在对模型划分网格时,网格单元的尺寸划分的越大,计算结果会越粗糙,也会影响波的传播,造成波传播失真;网格单元的尺寸划分的越小,计算结果会越准确,然而计算的时间和耗费的资源就会越多。因此,合理的选取划分网格时网格单元尺寸的大小对于准确的进行数值分析至关重要。
Kuhlemeyer 和 Lysmer(1973)建议,为了准确描述通过模型的波传播,网格单元的尺寸大小必须小于输入波的最高频率分量对应波长的约十分之一到八分之一,即:
式中,Δl 为网格单元的最大长度,λ 为通过介质的峰值速度最高频率分量的波长,c为波的传播速度,f为波的频率。
由于任何离散单元都有其允许通过的峰值频率的限制,因而当输入动力荷载波形具有尖脉冲(即高峰值且短时间起跳)特征时,为满足Kuhlemeyer和 Lysmer对网格的要求, 就必须采用非常精细的网格划分以及非常小的时步,这就会耗费大量的计算时间以及占用大量的内存。在这种情况下,可以采用FFT (快速傅里叶变换)技术,过滤掉高频率波,而不会显著影响模拟结果,因为能量绝大部分集中在低频率波上,滤掉的高频率波携带的能量是非常少的。
4. 3DEC动力分析边界条件的设置
(1) 无反射边界
岩土力学问题的建模边界,从分析尺度上而言应该是无界的。因此深埋地下开挖问题,可以假定洞室被无限围岩介质空间包围,而地表或近地表结构则认为是处于半无限空间。数值方法依靠有限空间离散网格来分析问题,必须在人工截断边界上施加合适的边界条件。在静态分析时,可在关注区域一定范围外设置约束边界或者采用弹性边界单元技术实现。
但在动力问题分析时,无论是固定边界还是自由边界,都会产生应力波的反射,这些反射波在模型内相互叠加,无法耗散。虽然采用较大的模型可以减少边界条件的影响,因为材料阻尼会吸收远处边界反射的能量,但这会导致模型很大,计算速度很慢。此时最实际的是采用无反射边界条件,该方法通过设置独立的阻尼器,在边界上波入射角超过30°时基本能完全吸收体波能量,而对面波只能近似,能量吸收并不完全。但该方法在时域分析时不失为一种有效方法,已经广泛用于各类数值模拟技术。
这是根据Lysmer和 Kuhlemeyer (1969)提出的黏滞边界理论来实现的,在模型边界法向与切向分别设置独立的阻尼器,来提供黏性法向与切向牵引力,以此来抵消反射波的应力。具体的理论公式可以参考帮助文档。
(2)自由场边界
在进行地表结构(坝体、边坡等)地震响应分析数值模拟时,需要对相毗邻基础范围进行离散化。地震波可以用平面波自下部材料向上传播来模拟。
此时模型两侧的边界必须引起结构内存在但缺失的自由场运动。自由场边界提供的条件与无限域模型完全一样,平面波向上传播不会受边界条件带来的失真影响。
主网格模型的侧向边界与自由场网格通过黏性阻尼器耦合在一起,来模拟静止边界。自由场网格上的不平衡力将施加到主网格边界上(具体原理和理论公式可以查看帮助文档)。一个施加了自由场边界的边坡的示意图如下图所示:
在自由场施加前,模型需要首先计算至静态平衡。当调用 FF apply
命令时,在模型侧向边界施加自由场条件,所有的临近自由场网格的单元信息,包括模型类型和当前状态参量同时复制至自由场区域。自由场的应力则指定为相邻节点单元的平均应力。而模型底部的动边界必须在施加自由场边界前施加。自由场网格是连续的,如果主网格内含有接触面延伸到边界,自由场内并不出现。
使用FF apply
命令将自由场网格施加后,自由场的网格会自动显示于块体绘图中,并可利用 LIST ff
查看。下图即为施加自由场边界的一个边坡模型:
注意:为了在3DEC中施加自由场边界,模型的底部必需水平,法向为 y 轴,侧向边界必须垂直,且法向平行于坐标轴 x 或 z。如果入射地震波不是竖向,必须旋转坐标轴使得应力波传播方向与y轴方向一致。如果模型竖直方向为Z轴,不能采用3DEC默认的自由场边界条件。
一般在模型底部采用粘滞边界(无反射边界),在模型四周采用自由场边界。这样做的好处就在于粘滞边界可以用来施加地震波和吸收反射波,自由场边界可以提供与无限场相同的效果,从而在进行动力计算时,地震波的能量可以通过模型四周的自由场边界耗散掉。
未完待续:下一节将介绍动力荷载的输入方式及手册中一个混凝土大坝的地震响应分析的示例。
下一节:3DEC动力分析示例(2):动力荷载输入方式及地震加载实例