FLAC动力分析

第11章 非线性动力反应分析

FLAC / FLAC3D可以进行非线性动力反应分析，而且具有强大的动力分析功能。本章以FLAC3D为例，详细介绍了动力分析过程中的边界条件、阻尼形式、荷载要求等，并通过一些实例对个别问题做了详细解答。

本章要点：

? FLAC动力分析与等效线性方法的差别 ? 动力分析时间步的确定方式及影响因素 ? 动态多步的概念

? 动力荷载的形式及施加方法 ? 动力边界条件的类型及适用条件

? 地震荷载输入的要点

? 三种阻尼形式的概念、参数确定及适用条件 ? 网格尺寸的要求 ? 输入荷载的校正 ? 地震液化的模拟

? 完全非线性动力分析的步骤 Equation Section 11

11.1 概述

FLAC / FLAC3D可以进行二维或三维的完全动力分析，FLAC/FLAC3D中的动力分析功能是可选模块，需要在程序中添加动力分析模块才可以进行。

FLAC3D中在动力分析前需要采用以下的命令：

CONFIG dynamic

对于FLAC，在程序开始时的Model Options对话框中选择Dynamic复选框。

FLAC / FLAC3D中的动力分析并不是只能孤立进行的，还可以与其他FLAC/FLAC3D元素进行耦合。 （1）与结构单元相耦合，可以用来进行土与结构的动力相互作用。

（2）与流体计算相耦合，可以模拟动力作用下土体孔隙水压力的上升直至土体液化。 （3）与热力学计算相耦合，可以计算热力荷载和动力荷载的共同作用。

（3）采用大变形计算模式，可以分析岩土体在动力荷载作用下发生的大变形。

FLAC和FLAC3D可以模拟岩土体在外部（如地震）或内部（如风、爆炸、地铁振动）荷载作用下的完全非线性响应，因此可以适用于土动力学、岩石动力学等学科的计算。

本章将以FLAC3D为例讨论动力计算的相关内容，FLAC的动力分析可以参照执行。

注意：FLAC和FLAC3D的动力计算十分复杂，读者在阅读本章内容之前要对FLAC3D的静力计算、流体计算十分熟悉，具体可以参阅本书的第7章和第12章的内容。

对于初次接触FLAC3D动力计算的读者，大多数都会提以下2个问题： （1）FLAC3D动力分析与一般的等效线性方法有什么区别？

（2）FLAC3D动力分析怎么会采用静力本构模型，比如Mohr-Coulomb模型？ 下面就这两个问题展开初步的讨论。

11.1.1 与等效线性方法的关系

在岩土地震工程中，等效线性方法广泛应用于计算地基土体中波的传播及土与结构的动力相互作用。该方法已被工程师、科研人员广泛接受。而FLAC3D采用的完全非线性方法没有获得广泛使用，因此需要对这两种方法之间的差异做简要介绍。

1. 等效线性方法的特点

等效线性方法的基本原理是，假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。目前用于土动力分析的等效线性模型已有数种，根据骨干曲线的形状可以分为双直线模型、Ramberg-Osgood模型、Hardin-Drnevich模型等，其中又以Hardin模型使用最多。等效线性方法有如下的特点：

? 使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。

? 不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。

? 塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。

? 大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。

? 本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。

2. FLAC3D非线性方法的特点

FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。相对于等效线性方法而言，完全非线性分析方法主要有以下优点：

? 可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。

? 采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。

? 由于采用了弹塑性模型，因此程序可以自动计算永久变形。 ? 采用合理的塑性方程，使得塑性应变增量与应力相联系。

? 可以方便地进行不同本构模型的比较。

? 可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦

合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。

另外，FLAC3D3.0已将等效线性方法中的模量衰减曲线以阻尼的形式嵌入到程序当中（见本章11.6.3节），使得FLAC3D的动力分析结果更易于被岩土地震工程师们所接受。

11.1.2 FLAC3D动力计算采用的本构模型

FLAC3D的动力计算可以采用任意的本构模型，比如弹性模型、Mohr-Coulomb模型。这一点很多读者都不能接受，他们普遍认为Mohr-Coulomb是静力本构模型，不适合用于动力分析，而应当采用更合适的Hardin模型。

其实这是对FLAC3D动力计算的误解。FLAC3D的原理是求解动力方程，所以从其算法上来说，不管是进行静力分析还是动力分析，其实质都是求解运动方程。只是对于静力分析而言，采用了特定的阻尼方式以达到快速收敛的目的。所以，有的场合将FLAC3D的静力分析方法称为“拟动力方法”。相应的，FLAC3D在进行动力分析时，通过求解动力方程理所当然地可以得到合适的动力问题解答。对于本构模型的选择，主要是描述单元的应力-应变关系，如果是弹塑性的，则考虑的是单元的屈服准则、流动法则等。

等效线性方法考虑土体的滞后性常常是通过将骨干曲线进行变换，比如Masing二倍法，而在FLAC3D

的动力分析中，滞后性是通过阻尼来考虑，通过设置合适的阻尼形式和阻尼参数，同样可以描述土体在动力作用下的滞回曲线和滞回圈。

因此，FLAC3D动力分析中采用的本构模型可以选取任意模型，其参数也是对应静力本构模型的参数，关键是要设置合适的阻尼形式、阻尼参数、边界条件等，这些内容将在本章的后续内容中进行讲解。

11.2 动力时间步

动力计算中临界计算时间步的计算如下：

?tcrit??V???min?? fCA??pmax??（11-1）

其中，Cp为p波波速，与材料的体积模量K和剪切模量G有关，可以表示为：

K?4/3G

Cp?? （11-2）

f为与四面体子单元相关的最大表面积，min{}表示遍历V为四面体子单元（sub-zone）的体积，Amax所有的单元，包括结构单元和接触面单元。

由于式（11-1）只是临界时间步的一个估计值，因此在使用中采用了一个安全系数，乘以0.5。因此，当采用无刚度比例的阻尼时，动力分析的时间步为：

?td??tcrit2

（11-3）

如果采用了刚度比例的阻尼，那么为了保持数值稳定性，时间步必须减小。Belytschko（1983）提出了一个临界时间步的公式?t?，其中考虑了刚度比例阻尼的影响。

?2??t??????max??1??2???

（11-4）

其中，?max为系统的最高特征频率，?为该频率下的临界阻尼比。

注意：FLAC3D在动力计算中，程序会根据数值计算的稳定性自动设置动力计算时间步，一般不建议读者对这个默认的时间步进行放大。甚至，在大应变计算过程中，如果出现很大的网格变形并导致网格的几何错误时，还要对默认的时间步进行折减，降低动力时间步，以达到数值稳定的目的。

11.3 动态多步

由式（11-1）可知，FLAC3D动力计算中时间步需要遍历所有单元，取所有单元临界时间步中的最小值，因此时间步是由几何尺寸较小、模量较大的单元来确定的。因此，在计算中，尤其是在试算期间，要尽量避免较小的单元尺寸及刚度很大的材料，比如用实体单元来模拟较薄的混凝土墙，这种情况下必然会使动力时间步非常小，从而造成计算时间过长。可以通过采用结构单元或暂时不考虑混凝土墙的办法来进行试算，等到有关参数调试完成后再进行细化计算。

当计算模型中存在刚度差异较大、模型网格尺寸不均匀的情况时，FLAC3D可以采用“动态多步”（Dynamic Multi-stepping）的过程来有效减少计算所需要的时间。在此过程中，模型单元和节点按照相近最大时步进行分组和排序，然后每个组在特定的时步下运行，信息在适当的时候在单元之间进行交换。

动态多步的调用采用如下命令：

SET dyn multi on

下面用一个简单的例子来描述动态多步的应用效果。同时读者可以从例子中了解到利用FISH函数来编写简单的动力荷载的方法。

1. 问题描述

如图11-1所示，土体的深度为10 m，挡土墙的高度为5 m，两者的模量差异为20倍。动力荷载从模型底部输入，主要分析目的是了解动态多步对计算时间的影响。

墙体 K = 4000 MPa

5 m

G = 2000 MPa

10 m

土体 K = 200 MPa G = 100 MPa

图11-1 动态多步作用的实例

本例计算中不考虑重力的影响，因此不用进行初始应力设置和平衡，直接进行动力计算。动力荷载采用正弦函数，采用FISH函数的方法进行定义，可以方便修改荷载的频率（freq）、幅值等。

2. 命令流

例11.1：动态多步实例

new

conf dyn ;打开动力计算功能 gen zone brick size 10 5 10 mod elas

mod null range x=0,5 z=5,10 ;删除部分网格 fix z range x=-.1 .1 z=.1 10.1 ;设置静力边界条件 fix z range x=9.9,10.1 z=.1 10.1 fix y range y=-.1 .1 fix y range y=4.9 5.1

prop bulk 2e8 shear 1e8 ;设置土体参数

prop bulk 4e9 shear 2e9 range x=5,6 z=5,10 ;设置墙体参数（土体参数的20倍） ini dens 2000 ;设置密度

def setup ;动荷载中的变量赋值 freq = 1.0

omega = 2.0 * pi * freq old_time = clock end

1

1

clock是FISH变量，表示计算机目前的时间

本文来源：https://www.bwwdw.com/article/yakv.html