第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 动力时间步

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

??V???tcrit?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 了一个临界时间步的公式?t?，其中考虑了刚度比例阻尼的影响。

（11-3）

如果采用了刚度比例的阻尼，那么为了保持数值稳定性，时间步必须减小。Belytschko（1983）提出

?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变量，表示计算机目前的时间

setup ;执行变量赋值 def wave ;定义动荷载函数 wave = sin(omega * dytime) ;定义动荷载变量 end

apply xvel = 1 hist wave range z=-.1 .1 ;施加动荷载 apply zvel = 0 range z=-.1 .1 hist gp xvel 5,2,0 hist gp xvel 5,2,10 hist gp zvel 5,2,10 hist dytime

def tim ;估算程序运行的时间 tim = 0.01 * (clock - old_time） end

set dyn multi on ;设置动态多步 solve age 1.0

print tim ;输出计算时间 print dyn ;输出动力计算相关信息 save mult1.sav

注意：动力计算中必须设置材料的密度，若模型中存在结构单元，也必须设置结构单元的密度，否则会出错。采用FISH函数定义动力荷载时，FISH函数和变量应具有相同的名称。因为设置动力边界条件命令中的hist关键词后面必须要跟随一个FISH函数名，FISH变量需要调用FLAC3D中的内置标量dytime，该变量是动力计算的真实时间，通过调用可以给函数提供预订变化的数值。所以，一般FISH定义动荷载的方法如下（以定义函数xxx为例）：

def xxx xxx = …dytime end

app xvel = 1.0 hist xxx range …

3. 计算过程与输出结果

计算过程中命令窗口会提示动力计算的步数、动力时间和时间步，计算结束后可以将模型底部和墙体顶部节点的水平速度时程输出，使用以下的命令：

plot hist 1 2 v 4

输出结果见图11-2所示。读者可以采用设置动态多步和不设置动态多步两种情况分别进行计算，观察图11-2中的速度时程曲线以及FLAC3D命令窗口中输出结果（图11-3和图11-4）。可以发现，时程曲线、动力计算的时间步、迭代步数均一致，不同的是设置动态多步的情况下花费的计算时间较少。

注意：程序运行的时间将由于计算机配置的不同而存在差别。

从图11-3中的输出信息可以看出，动态多步将模型中的375个单元分成了两类，并提供了不同的时步乘子，其中只有小部分单元（墙体单元65个）的时步乘子为1，其他大部分单元（土体单元310个）拥有较高的时步乘子，这样可以大大加快计算的进度。

FLAC3D 3.00Step 821616:46:10 Wed Apr 02 2008 1.2History 1 X-Velocity Gp 55 Linestyle -1.000e+000 <-> 1.000e+000 2 X-Velocity Gp 688 Linestyle -1.192e+000 <-> 1.340e+000 1.0 0.8 0.6 Vs. 4 Dynamic Time 1.217e-003 <-> 9.993e-001 0.4 0.2 0.0-0.2-0.4-0.6-0.8-1.0Itasca Consulting Group, Inc.Minneapolis, MN USA 2.0 4.0x10^-1 6.0 8.0 图11-2 动力计算结束时模型底部和顶部的水平速度时程曲线

图11-3 设置动态多步情况下的输出信息

图11-4 未设置动态多步情况下的输出信息

11.4 动力荷载和边界条件

利用FLAC3D进行动力计算时，有以下3个问题需要读者认真考虑： ? 动力荷载和边界条件； ? 力学阻尼；

? 模型中波的传播。

本节及后续的两节将分别针对以上三个问题展开讨论。

11.4.1 动力荷载的类型与施加方法

FLAC3D可以在模型边界或内部节点施加动荷载来模拟材料受到外部或内部动力作用下的反应，程序允许的动力荷载输入可以是：（1）加速度时程，（2）速度时程，（3）应力（压力）时程，（4）集中力时程。

动力荷载的施加采用APPLY命令，另外，采用APPLY Interior命令可以将加速度、速度和力的时程施加到模型内部的节点上。动力荷载的形式主要有2种：

? FISH函数。FISH函数表达的动力荷载往往比较规则，也常用于试算阶段的动力输入，因为试

算时可以不用过多考虑荷载的频率、基线校正等问题。本章例11.1中已经给出了一个FISH函数作为动力荷载的例子，这里不再赘述。

? TABLE命令定义的表。常用于离散的动力荷载输入，包括地震波、实测振动数据、不规则动力输入等。下面简要介绍利用TABLE命令形成的表作为动力荷载的方法。

表是FLAC3D中的一种数据格式，表中的数据成对出现，相当于两列的表格。表建立的基本命令是：

TABLE n x1 y1 x2 y2 x3 y3

其中n表示表的ID号，（x1,y1）、（x2,y2）、（x3,y3）分别为表格中的三对数据，例如在命令行中输入如下命令：

table 1 1 1 2 4 3 6

表示建立了一个ID号为1的表，表中有3对数据。可以通过PLOT命令绘出该表的图形：

plot table 1 line

也可以通过PRINT命令打印该表的内容：

print table 1

在FLAC3D动力计算中，动力荷载往往数据点很多，用命令输入的方法显然不便，因此常用编辑文本文件的方法进行表的读入与调用，编辑文本文件的表有2种格式：

? x列均匀间隔的表，常用于等间隔的动力荷载形式。

第1行：表的名称

第2行：数据对的个数 空格 时间间隔（x列的数据间隔） 第3行：y列的第1个数据 第4行：y列的第2个数据 …… 空行

? 分别给出x，y数据对的表

第1行：表的名称 第2行：x1 空格 y1 第2行：x2 空格 y2 …… 空行

注意：在表的文本文件最后，需要有一个回车换行符，否则会出现“Error reading file xxx.dat”的错误；表的名称可以用英文、中文，也可以包含空格；表的文本文件可以保存成.dat、.txt等格式。 完成的文本文件需要进行读入操作才可以供FLAC3D调用。采用TABLE read命令进行读入引用：

table 1（ID号） read 文件名

读入后，读者可以使用PLOT table和PRINT table命令来查看生成的表文件是否读入正确。在进行FLAC3D动力边界条件设置时，使用APPLY命令调用已读入的表，下面的命令以施加水平向速度荷载为例。

app xvel 1.0 hist table 1 range …

命令中的1.0表示表格y向数据的乘子，可以方便地控制荷载幅值的大小。另外，动荷载的输入可以沿着x，y，z的任意方向施加或者模型边界的法向和切向施加。

注意：特定的边界条件不能在相同的边界上进行混合施加，否则程序会提示出错。

11.4.2 边界条件的设置

在动力问题中，模型周围边界条件的选取是一个主要内容，因为边界上会存在波的反射，对动力分析的结果产生影响。把分析模型的范围设置得越大，分析结果就越好，但较大的模型会导致巨大的计算负担。FLAC3D中提供了静止（粘性）边界和自由场边界两种边界条件来减少模型边界上的波的反射。

1. 静态边界

FLAC3D中允许采用静态边界（也称粘性边界，吸收边界）条件来吸收边界上的入射波。FLAC3D中的静态边界是Lysmer和Kuhlemeyer（1969）提出的，具体做法是在模型的法向和切向分别设置自由的阻尼器从而实现吸收入射波的目的，阻尼器提供的法向和切向粘性力分别为式（11-5）和（11-6）。

tn???Cpvn

（11-5）

ts???Csvs

（11-6）

其中，vn，vs分别为模型边界上法向和切向的速度分量，?为介质密度，Cp，Cs分别为p波和s波的波速。

这种静态边界对于入射角大于30度的入射波基本能够完全吸收，对于入射角较小的波，比如面波，虽然仍有一定的吸收能力，但吸收不完全。静态边界可以加在整体坐标系上，也可以加载在倾斜边界的法向和切向上。如果在倾斜边界的法向和切向施加静态边界，则需要同时使用nquiet，dquiet，squiet条件。

整体坐标系的静态边界条件设置使用命令为：

APPLY xquiet (yquiet, zquiet) range …

使用倾斜边界上的静态边界条件命令为：

APPLY nquiet dquiet squiet range …

注意：

（1）施加动力边界条件后，这些边界上原先的静力边界条件将被自动去掉（free），在动力荷载施加期间，程序始终自动计算边界上的作用力，用户不能将静态边界条件去掉；也不能在静态边界上施加加速度、速度边界条件，因为静态边界上的作用力是根据边界上的速度分量计算得到的，如果再施加速度荷载就会使静态边界失效。若需要在静态边界上输入动荷载，则只能输入应力时程。可以将加速度、速度时程通过转换公式（11-7）和（11-8）形成应力时程施加到静态边界上。

?n=-2(?Cp)vn

（11-7） （11-8）

?s?-2(?Cs)vs

式中的?n，?s分别为施加在静态边界上的法向应力和切向应力，公式中的系数2表示施加的能量中只有一半是向上传播作为动力输入的，另一半向边界下部传播。公式中的负号是为了使应力施加后节点的速度能与实际一致。

（2）对于动力荷载来源于模型内部（如隧道中的列车振动问题）的情况，可以将动力荷载直接施加在节点上，这种情况下使用静态边界可以有效减小人工边界上的反射，并且不需要再施加下面提到的自由场边界。

（3）动力计算过程中应避免静力荷载的变化。比如，在一个已经施加底部静态边界的模型上进行开挖，会造成整个模型向上移动。因为施加静态边界时程序自动计算了施加在模型底部边界上的反力，这些边界反力不能与开挖后的模型相平衡，会引起模型的整体上移。

下面给出一个简单的静态边界的例子。

问题描述：如图11-5所示，一根竖直的弹性杆，高50m，宽1m，体积模量和剪切模量分别为20MPa

3

和10MPa，材料密度为1000kg/m。在杆的底部边界的两个水平方向设置静态边界条件，而杆的顶部为自由表面，在杆的底部施加水平方向的应力冲击荷载。根据材料参数，可得计算得到剪切波速Cs为100 m/s，?Cs的乘积为105。应力冲击荷载的幅值设置为2×105，根据公式（11-7）可以得到等效速度荷载幅值为1m/s。计算中不考虑重力的作用，因此不需要进行初始应力的计算。

分析杆的底部、中心点、顶部三个典型节点的速度响应，计算结果如图11-6所示。可以发现模型底部产生的速度荷载幅值就是1m/s，图中最后两个脉冲是从模型顶部自由面反射回来的波。在模型顶部节点的速度幅值是底部输入幅值的两倍，而且顶部反射回来的波传到模型底部以后，模型基本结束了振动。这些现象说明采用静态边界的设置是合理有效的。

计算文件如下：

例11.2：静态边界的例子

new config dyn

gen zone brick size 1,1,50 model elas

prop shear 1e7 bulk 2e7 ini dens 1000 def setup

omega = 2.0 * pi * freq pulse = 1.0 / freq end set freq=4.0 setup def wave

if dytime > pulse wave = 0.0 else

wave = 0.5 * (1.0 - cos(omega * dytime)) endif end

range name bottom z=-.1 .1

fix z range z=.5 55 ;将上部网格都施加数值向约束 apply dquiet squiet range bottom

apply sxz -2e5 hist wave syz 0.0 szz 0.0 range bottom ;-2e5的系数来源于?Cs的值 apply nvel 0 plane norm 0,0,1 range bottom

hist gp xvel 0,0,0 hist gp xvel 0,0,25 hist gp xvel 0,0,50 hist dytime hist wave plot create hhh plot add hist 1 2 3 vs 4 plot show solve age 2

1 m

50 m

图11-5 静态边界条件示例

FLAC3D 3.00Step 147915:13:54 Thu Apr 03 2008 1.8History 1 X-Velocity Gp 1 Linestyle -3.836e-002 <-> 1.001e+000 2 X-Velocity Gp 101 Linestyle -3.348e-002 <-> 1.000e+000 3 X-Velocity Gp 201 Linestyle -5.125e-002 <-> 1.995e+000 1.6顶部 1.4 1.2底部 中心 Vs. 4 Dynamic Time 1.217e-002 <-> 1.789e+000 1.0 0.8 0.6 0.4 0.2 0.0Itasca Consulting Group, Inc.Minneapolis, MN USA 0.5 1.0 1.5 图11-6 在静态边界上输入应力荷载得到的响应

2. 自由场边界

对诸如大坝之类的地面结构进行动力反应分析时，在模型各侧面的边界条件须考虑为没有地面结构时的自由场运动。FLAC3D通过在模型四周生成二维和一维网格的方法来实现这种自由场边界条件（见图11-7所示），主体网格的侧边界通过阻尼器与自由场网格进行耦合，自由场网格的不平衡力施加到主体网格的边界上。由于自由场边界提供了与无限场地相同的效果，因此向上的面波在边界上不会产生扭曲。

图11-7 自由场边界示意图

对于FLAC3D而言，自由场边界模型包括4个平面网格和4个柱体网格，平面网格在模型边界上与主体网格是一一对应的，柱体网格相当于平面自由场网格的自由场边界。其中，平面自由场网格是二维计算，假设在面的法向无限延伸；柱体自由场网格是一维计算，假设在柱体两端无限延伸。

下面给出一个简单的自由场边界条件的例子，见例11.3。

例11.3：自由场边界条件的实例

new

;第一步：静力计算阶段 config dyn set dyn off

gen zone brick size 6 3 2 gen zone brick size 2 3 2 p0 0 0 2 gen zone brick size 2 3 2 p0 4 0 2 gen zone wedge size 1 3 2 p0 2 0 2

gen zone wedge size 1 3 2 p0 4 3 2 p1 3 3 2 p2 4 0 2 p3 4 3 4 & p4 3 0 2 p5 4 0 4 model elastic

prop bulk 66667 shear 40000 ini dens 0.0025 set grav 0 0 -10 fix x range x -0.01 0.01 fix x range x 5.99 6.01 fix y range y -0.01 0.01 fix y range y 2.99 3.01 fix z range z -0.1 0.1 hist unbal solve

save 11-3_1.sav ;第二步：动力计算阶段 set dyn on def iniwave per = 0.01 end iniwave def wave

wave = 0.5 * (1.0 - cos (2*pi*dytime/per)) end

free x y z ran z -0.1 0.1 ;去掉模型底部原有的静力条件 apply nquiet squiet dquiet ran z -0.1 0.1 ;静态边界条件 apply dstress 1.0 hist wave ran z -0.1 0.1 ;加动力荷载

apply ff ;施加自由场边界条件 group ff_corner

group ff_side ran x 0 6 group ff_side ran y 0 3 group main_grid ran x 0 6 y 0 3

set dyn time = 0 ;设置动力计算从0s开始 hist reset ;清空已有的历史信息

hist unbal hist dytime ; 主体网格 hist gp xvel 2 1 0 hist gp xvel 2 1 5.0 ; 柱体网格 hist gp xvel -1 -1 0 hist gp xvel -1 -1 5.0

; 平行于y方向的二维自由场网格 hist gp xvel -1 0 0 hist gp xvel -1 0 5.0

; 平行于x方向的二维自由场网格 hist gp xvel 2 -1 0 hist gp xvel 2 -1 5.0 solve age 0.015 save 11-3_2.sav

计算分两步，第一步进行重力作用下的初始应力计算，这里采用的是直接重力加载的方法生成初始应力，第二步进行动力计算。在设置动力边界条件时采用了4个步骤：

步骤1 去掉模型底部已有的静力约束条件； 施加静态边界条件；

步骤2 步骤3 施加动力荷载；

步骤4 施加自由场边界条件。

施加自由场边界条件以后，程序自动在主体网格周围生成一圈自由场网格。当计算模型的单元数较多时，这个生成过程需要花费较多的时间。可以对生成的自由场网格定义group，如例11.3中就将二维自由场网格和一维自由场网格分别赋值了group，这样便于模型的观察和主体网格结果的后处理。

图11-8 自由场边界施加前后的网格

计算中对主体网格、柱体自由场网格和平面自由场网格对应位置的水平向速度进行了监测，计算结果见图11-9所示。可以发现，主体网格与周围的自由场网格同步运动，达到了自由场网格的目的。

FLAC3D 3.00Step 99217:04:15 Thu Apr 03 2008x10^-2 9.0FLAC3D 3.00Step 99217:04:45 Thu Apr 03 2008x10^-1 1.0History 3 X-Velocity Gp 10 Linestyle 6.711e-006 <-> 9.060e-002 5 X-Velocity Gp 368 Linestyle 6.711e-006 <-> 9.049e-002 7 X-Velocity Gp 145 Linestyle 6.711e-006 <-> 9.049e-002 9 X-Velocity Gp 232 Linestyle 6.711e-006 <-> 9.063e-002 8.0History 4 X-Velocity Gp 102 Linestyle -6.480e-004 <-> 1.005e-001 6 X-Velocity Gp 380 Linestyle -3.849e-005 <-> 1.001e-001 8 X-Velocity Gp 176 Linestyle -3.849e-005 <-> 1.001e-001 10 X-Velocity Gp 274 Linestyle -6.919e-004 <-> 1.007e-001 0.8 7.0 6.0 0.6 5.0 Vs. 2 Dynamic Time 5.735e-005 <-> 1.497e-002 4.0 Vs. 2 Dynamic Time 5.735e-005 <-> 1.497e-002 0.4 3.0 2.0 0.2 1.0 0.0Itasca Consulting Group, Inc.Minneapolis, MN USA 0.2 0.4 0.6 0.8x10^-2 1.0 1.2 1.4 Itasca Consulting Group, Inc.Minneapolis, MN USA 0.2 0.4 0.6 0.8x10^-2 1.0 1.2 1.4 模型底部 模型顶部

图11-9 自由场边界条件的示例结果

讨论：读者可以修改例11.3中边界条件设置的命令，比如去掉某个命令或者调换命令的前后位置，可以得到以下结论：

? free命令在本例中可以删去，因为接下来设置静态边界本身就会将模型底部边界上已有的静力

条件设置为free。但是如果静力计算阶段模型底部的边界条件设置为fix x y z，而且施加的动力

荷载形式为加速度或者速度，那么就需要首先将这些fix速度的静力条件设置为free，否则会出现动力加载的错误。

? 本例中app ff的位置对结果没有影响，但是在动力计算中仍然要把app ff放在所有边界条件的

后面，换句话说，静态边界条件、动力荷载施加等必须在自由场边界条件设置之前完成。 另外，使用自由场边界条件还需要注意以下几点：

? 如果动力源只存在于模型的内部，那么可以不必设置自由场边界。

? 自由场边界设置对模型的形状有一定的要求：模型底部水平，重力方向为z向，四个侧面垂直，

法向分别为x、y向。如果地震波的传播方向不是竖直方向的，则应当进行坐标轴旋转，使得z

轴与地震波传播方向相同。这种情况下，重力方向将与z轴存在一定的夹角，而且模型边界也与水平面产生倾斜。

? 其他边界条件应该在APPLY ff之前进行设置。

? APPLY ff将边界上单元的属性、条件和变量全部转移至ff单元上，并且设置以后主体网格上的改动将不会被FF边界所响应，这一点与静态边界类似。

? 使用自由场边界对主体网格的本构模型没有要求，并可以进行竖向流体计算相耦合。 ? 自由场边界进行的是小变形计算，主体网格可进行大变形计算，因此自由场边界上的变形要相

对较小，如果自由场网格上的变形较大，那么需要扩大人工边界条件的选取。 ? 存在attach的边界将不能设置FF边界，而且主体网格边界上的Interface将不能连续到自由场

网格。 ? 可以对生成的自由场网格进行group赋值，这样在后处理的云图显示时去掉自由场网格，生成

只有主体网格的图形。

11.4.3 地震荷载的输入

地震反应分析是FLAC3D动力计算的主要应用领域，所以有必要对地震荷载的输入做单独介绍。地震荷载常用表的文本文件格式读入到FLAC3D，再施加到模型底部的边界上。对于地震荷载的动力边界

条件选择可以参考以下标准：

1. 刚性地基

如果模型底部为岩石等模量较大的材料，可以在底部直接施加加速度或速度荷载，并采用自由场边界条件，模型底部无需施加静态边界条件。

2. 柔性地基

如果模型底部的单元为土体，尤其是软土，则不能直接施加加速度和速度，而需要将加速度、速度转换成应力时程，再施加到模型底部。模型周围采用自由场边界条件，模型底部可采用静态边界条件（nquiet，dquiet，squiet）。

11.5 力学阻尼

阻尼的产生主要来源于材料的内部摩擦以及可能存在的接触表面的滑动。FLAC3D采用求解动力方程的方法解决两类力学问题：准静力问题和动力问题。在这两类问题中都要使用阻尼，但准静力问题需要更多的阻尼使得动力方程能够更快的收敛平衡。对于动力问题中的阻尼，需要在数值模拟中重现自然系统在动荷载作用下的阻尼大小。目前，FLAC3D动力计算提供了三种阻尼形式供用户选择，分别是瑞利阻尼、局部阻尼和滞后阻尼。

11.5.1 瑞利阻尼（Rayleigh damping）

瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关：

C??M??K （11-9）

其中，?为与质量成比例的阻尼常数，?为与刚度成比例的阻尼常数。

瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。图11-10为归一化的临界阻尼比与角频率之间的关系曲线，其中包含了仅设置刚度分量、仅设置质量分量以及二者叠加的三种结果。可以看出，采用叠加的方法得到的阻尼比在较大的频率范围内保持定值，因此瑞利阻尼可以近似地反映岩土体具有的频率无关性。

叠加 仅有刚度分量

?=0

仅有质量分量

?=0

图11-10 归一化的临界阻尼比与角频率之间的关系

1. 瑞利阻尼的两个参数 设置瑞利阻尼的命令为：

SET dyn damp rayleigh 参数1 参数2

根据图11-10叠加结果的阻尼比曲线最小值可以确定瑞利阻尼的两个参数，分别是最小临界阻尼比

，可以按照式（11-10）进行计算： ?min（参数1）和最小中心频率?min（参数2）

?min?(???)1/2?min?(?/?)1/2 （11-10）

注意：参数2中的单位是Hz，而不是角频率的单位。 2. 瑞利阻尼参数选择的原则

很多读者最十分关心瑞利阻尼参数的选择。对于岩土材料而言，临界阻尼比的范围一般是2～5%，而结构系统的临界阻尼比为2～10%。在使用弹塑性模型进行动力计算时，相当多的能量消散于材料发

生的塑性流动阶段，因此在进行大应变的动力分析时，只需要设置一个很小的阻尼比（比如0.5%）就能满足要求，而且达到塑形以后，随着应力-应变滞回圈的扩大，能量的消散也逐渐明显。

（1）中心频率的确定

瑞利阻尼是与频率相关的，但是在一定的频率范围内，瑞利阻尼基本与频率无关。这个频率范围的最大频率常是最小频率的3倍。对于任何动力问题，可以对速度时程进行谱分析得到速度谱与频率之间的关系，见图11-11所示。可以逐渐调整fmin，使得频率范围在fmin～3×fmin之间包含了动力能量的主要部分，这时的fmin就是瑞利阻尼的中心频率。在实际分析中，比如多种材料的土石坝动力分析，需要首先将分析模型假设为弹性材料进行动力计算，得到各种材料关键位置的功率谱，根据功率谱的分布来确定该区域的瑞利阻尼中心频率的大小。通过这种方法确定的中心频率fmin，既不是输入频率，也不是系统的自振频率，而是二者的叠加。

主频区域

速度谱

?min

3×?min

频率

图11-11 速度谱与频率之间的关系

对于简单的模型，也可以采用自振频率作为瑞利阻尼的中心频率，这需要对模型进行无阻尼的自振计算，下面给出一个实例。

自振频率的计算要求设置正确的边界条件，不设置阻尼，在重力作用下求解一定的步数，使模型产生振荡，分析模型中关键节点的响应（包括速度、位移随动力时间的变化曲线）。具体的求解步数可以根据输出结果，至少完成一个或几个周期的振荡，从而确定自振频率的大小。

例11.4：自振频率的计算

conf dy

gen zone brick size 3,3,3 model elas

prop bulk 1e8 shear 0.3e8 ini dens 1000 fix z range z -.1,.1

set dyn=on, grav 0 0 -10, hist_rep=12 hist gp zdisp 3.0,1.5,3.0 hist dytime plot create hh plot add his 1 vs 2

save 11-4_1.sav ;保存文件，在后续计算中需调用该文件 cyc 150

例11.4中记录了模型边界上一个节点的位移振荡曲线，见图11-12所示。可以发现完成一个振荡周期需要的时间约为0.0438s，据此可以计算出系统的自振频率为22.8Hz。读者在确定振荡周期时，可以

变量hist_rep的作用是将设置历史记录的间隔，默认为10，即每10步记录一次数据。由于本例中计算时间步较少，为了得到更平滑的记录曲线，将hist_rep取为1。

2

将计算结果输出到文本文件，再通过Excel操作得到更精确的周期值。

FLAC3D 3.00Step 15011:05:17 Fri Apr 04 2008Job Title: vertical displacement versus time (undamped)x10^-3History 1 Z-Displacement Gp 56 Linestyle -1.091e-003 <-> -3.527e-006-0.1-0.2 Vs. 2 Dynamic Time 5.939e-004 <-> 8.909e-002-0.3-0.4-0.5-0.6-0.7-0.8-0.9-1.0Itasca Consulting Group, Inc.Minneapolis, MN USA 2.0 4.0x10^-2 6.0 8.0 图11-12 自振频率计算得到的位移时程曲线

（2）阻尼比的确定

经验方法是，直接选取岩土体的阻尼比参数，比如0.5%。

另外一种，通过弹性阶段的动力计算，了解各关键位置的动应变幅值，并根据实验室得到的阻尼比-应变幅值曲线来选择阻尼比的大小。

3. 瑞利阻尼的实例

以例11.4中保存的11-4_1.sav文件为基础，采用如下3种参数形式的瑞利阻尼进行动力计算： ? 质量分量和刚度分量共同作用； ? 仅有质量分量； ? 仅有刚度分量。

为了更好地进行比较，上述三种参数均调试为临界阻尼状态，即模型刚好不做周期振动而回到平衡位置。达到临界阻尼状态必须将临界阻尼比设为1，使用自振频率作为瑞利阻尼的中心频率，而且同时使用质量分量和刚度分量，即例11.5的第（1）种情况。

如果仅采用质量阻尼或刚度阻尼，则需要将阻尼比设置为共同作用下的2倍，具体命令文件如下。

例11.5：瑞利阻尼的例子

;（1）质量分量和刚度分量共同作用 rest 11-4_1.sav

set dyn damp rayleigh 1 22.8 solve age=0.2 title

vertical displacement versus time (mass & stiffness damping) plot show pause3

;（2）只有质量分量

命令pause的作用是在命令文件执行中设置暂停，以便于用户观察暂停状态下的计算结果，在命令行中输入continue可以继续执行，若输入命令中存在错误，就停止执行。pause后面可以跟keyword关键词，表示暂停状态下按任意键即可继续。

3

rest 11-4_1.sav

set dyn damp rayleigh 2 22.8 mass solve age=0.08 title

vertical displacement versus time (mass damping only) plot show pause

;（3）只有刚度分量 rest 11-4_1.sav

set dyn damp rayleigh 2 22.8 stiffness solve age=0.08 title

vertical displacement versus time (stiffness damping only) plot show

三种计算结果得到的竖向位移时程曲线均如图11-13所示，表明三种情况都达到了临界阻尼状态。但是三种计算工况下的时间步存在较大差别，具体见表11-1所示。可以看出，瑞利阻尼中刚度分量的设置大大降低了FLAC3D中的动力时间步。

表11-1 瑞利阻尼三种参数格式下的计算时间步比较

计算工况 （1）质量分量和刚度分量共同作用 7.8E-5 （2）只有质量分量 （3）只有刚度分量 时间步 5.9E-4 3.9E-5 FLAC3D 3.00Step 202515:29:37 Fri Apr 04 2008Job Title: vertical displacement versus time (stiffness damping only)x10^-4History 1 Z-Displacement Gp 56 Linestyle -5.371e-004 <-> -1.561e-008-0.5-1.0 Vs. 2 Dynamic Time 3.951e-005 <-> 8.001e-002-1.5-2.0-2.5-3.0-3.5-4.0-4.5-5.0Itasca Consulting Group, Inc.Minneapolis, MN USA 1.0 2.0 3.0 4.0x10^-2 5.0 6.0 7.0 8.0 图11-13 竖向位移时程曲线的临界阻尼状态

11.5.2 局部阻尼（Local damping）

局部阻尼是FLAC3D静力计算中采用的阻尼形式，但是它的一些特性可以用来进行动力计算。它在振动循环中通过在节点或结构单元节点上增加或减小质量的方法达到收敛，由于增加的单元质量和减小的相等，因此总体来说，系统保持质量守恒。当节点速度的符号改变时增加节点质量，当速度达到最大值（或最小值）时减小节点质量。因此损失的能量?W是最大瞬时应变能W的一定比例，这个比值?W/W，是与频率无关和率无关的。因为?W/W是临界阻尼比D的函数：

局部阻尼的设置命令为：

SET dyn damp local 局部阻尼系数

?L??D

其中，?L为局部阻尼系数，临界阻尼比的取值可以参考瑞利阻尼中的?min。

（11-11）

局部阻尼系数不用求解系统的自振频率，而且相对于瑞利阻尼而言不会减小时间步，从这个意义上来说具有较大的优势。但局部阻尼只适合于简单问题的求解，实践证明设置局部阻尼不能有效地衰减复杂波形的高频部分，计算结果会产生一些高频“噪声”。因此，在使用时需慎重，最好将局部阻尼与瑞利阻尼的结果进行对比分析。

继续沿用例11.4中的保存文件11-4_1.sav，比较瑞利阻尼和局部阻尼的计算结果。瑞利阻尼中临界阻尼比设为5%，根据公式（11-11）可知局部阻尼的阻尼系数为0.1571。

例11.6：局部阻尼的例子

rest 11-4_1.sav

set dyn damp rayleigh 0.05 22.8 set hist_rep=5 solve age=0.5 title

vertical displacement versus time (5% Rayleigh damping) plot show pause rest 11-4_1.sav

set dyn damp local 0.1571 ; = pi * 0.05 set hist_rep=5 solve age=0.5 title

vertical displacement versus time (5% Local damping) plot show

计算结果见图11-14所示，可以看出两种阻尼形式具有相同的计算结果。另外，从时间步上看，瑞利阻尼的时间步为4.9E-4，而局部阻尼的时间步为5.9E-4，局部阻尼的时间步略大于瑞利阻尼，但二者相差并不明显。这是因为瑞利阻尼中质量分量的比例（?值）较小，读者可以通过PRINT dynamic输出动力计算信息了解瑞利阻尼中两种分量的比例。

FLAC3D 3.00Step 142016:02:52 Fri Apr 04 2008Job Title: vertical displacement versus time (5% Rayleigh damping)x10^-4FLAC3D 3.00Step 117916:02:35 Fri Apr 04 2008Job Title: vertical displacement versus time (5% Local damping)x10^-4-1.0-1.0History 1 Z-Displacement Gp 56 Linestyle -9.881e-004 <-> -3.623e-005-2.0History 1 Z-Displacement Gp 56 Linestyle -9.888e-004 <-> -4.737e-005-2.0 Vs. 2 Dynamic Time 2.466e-003 <-> 7.002e-001-3.0 Vs. 2 Dynamic Time 2.970e-003 <-> 6.978e-001-3.0-4.0-4.0-5.0-5.0-6.0-6.0-7.0-7.0-8.0-8.0-9.0-9.0Itasca Consulting Group, Inc.Minneapolis, MN USA 1.0 2.0 3.0x10^-1 4.0 5.0 6.0 7.0

Itasca Consulting Group, Inc.Minneapolis, MN USA 1.0 2.0 3.0 4.0x10^-1 5.0 6.0 瑞利阻尼 局部阻尼

图11-14 瑞利阻尼与局部阻尼的比较

11.5.3 滞后阻尼（Hysteretic damping）

FLAC3D将土动力学中岩土体的滞后特性用阻尼的形式加入到程序中。使用模量衰减系数Ms来描述土体的非线性特性。假设土体为理想粘弹性体，可以从模量衰减曲线上得到归一化的剪应力?：

示为：

Mt???Ms? dMsd? ?Ms??d?d?（11-12）（11-13）

式中，?为剪应变，Ms为归一化的割线模量，Mt为归一化的切线模量。则增量剪切模量G可以表

其中，G0为小应变下的剪切模量。

G?G0Mt

（11-14）

图11-15 模量衰减曲线（Seed & Idriss）

滞后阻尼是与材料无关的阻尼格式，在动力计算中，滞后阻尼可以满足Masing二倍法，从而构造土

体在动力作用下的滞回圈。另外，滞后阻尼的优点是：

? 可以直接采用动力试验中的模量衰减曲线；

? 相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步； ? 可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。

滞后阻尼与瑞利阻尼及局部阻尼的设置不同，采用的是初始条件INITIAL命令。

INITIAL damp hyst

FLAC3D中滞后阻尼提供了多种形式的割线模量衰减曲线模型，包括默认模型（default）、S型模型（包括三参数模型Sig3和四参数模型Sig4）以及哈丁模型（Hardin）。下面对这几种模型做简要介绍。

1. 默认模型

默认模型中Ms曲线可以用式（11-15）的三次方程来拟合。

其中

s?L2-L L2-L1Ms=s2(3?2s)

（11-15）

（11-16） （11-17）

L?log10(?)

L1和L2为默认模型的两个参数，表示Ms曲线的循环应变范围。比如L1=-3，L2=1表示Ms曲线中循

环应变的最小值是0.001%（10-3），最大值是10%（101），因此可见默认模型的参数确定较简单，设置方法是：

INITIAL damp hyst default L1 L2

2. S型模型

S型模型包括三参数模型和四参数模型，采用如下的公式来拟合Ms曲线。 Sig3模型包含a，b，x0三个参数：

设置命令为：

INITIAL damp hyst sig3 a b x0

Ms?a

1?exp(-(L-x0)/b)（11-18）

Sig4模型包含a，b，x0和y0四个参数：

设置命令为：

INITIAL damp hyst sig4 a b x0 y0

Ms?y0+a

1?exp(-(L-x0)/b)（11-19）

3. 哈丁模型

滞后阻尼中有一种Hardin / Drnevich模型，采用式（11-20）的双曲线公式来拟合Ms曲线：

Ms?1

1??/?ref（11-20）

其中，?ref为参考应变，一般取G/Gmax=0.5时的对应的应变值。 哈丁模型只有一个参数，设置命令为： INITIAL damp hyst hardin ?ref

表11-2提供了图11-15中模量衰减曲线的拟合结果，表11-3提供了黏土的模量衰减曲线拟合结果，

读者在使用时可参考。

表11-2 Seed & Idriss模量衰减曲线的拟合结果

数据来源 默认模型 三参数模型 四参数模型 a = 0.9762 a = 1.014 砂土 图11-15 L1 = -3.325 b = -0.4792 L2 = 0.823 x0 = -1.249 y0 = 0.03154 x0 = -1.285 b = -0.4393 哈丁模型 ?ref = 0.06 表11-3 Seed & Idriss模量衰减曲线的拟合结果（2）

数据来源 默认模型 三参数模型 四参数模型 a = 0.922 a = 1.017 L1 = -3.156 粘土 L2 = 1.904 x0 = -0.633 y0 = 0.0823 b = -0.587 x0 = -0.705 b = -0.481 哈丁模型 ?ref = 0.234 注意：滞后阻尼是FLAC3D新开发的一个技术，在使用过程中读者要反复调试，同时也要注意以下几点：

? 低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个较小分量的Rayleigh阻尼（比如0.2%刚度比例），这样也不会降低时步。

? 若初始剪应力不为0，剪应力-剪应变曲线可能不匹配，因此在生成初始应力时就要调用Hyst阻尼，这一点至关重要。

? 滞后阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，可能会导致动力响应的幅值偏大。

? 在设置滞后阻尼之前要做一次弹性无阻尼求解，获得各关键部位发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么使用滞后阻尼可能会存在问题。

? 即使应变较小，使用塑性模型也会增大应变，因此若模型存在广泛屈服的现象时，不能使用滞后阻尼。

11.5.4 关于阻尼设置的一些讨论

力学阻尼的设置是FLAC3D动力计算中讨论最多的话题，在这些阻尼形式中，瑞利阻尼由于其理论

与常规动力分析方法类似，而且实践证明，瑞利阻尼计算得到的加速度响应规律比较符合实际，因此最易为大家所接受，唯一也是最大的不足就是瑞利阻尼的计算时间步太小，导致动力计算时间过长，因此很多用户不得不使用局部阻尼来代替。

滞后阻尼是新版本FLAC3D的一个亮点，但在实际使用过程中读者可以发现使用时存在一定的困难，主要原因是滞后阻尼有过多的使用限制，且目前相关的参考资料极少。作者曾经尝试过一些滞后阻尼的算例，当模型较复杂时，很难得到满意的分析结果，尤其是在处理地震液化分析时更是如此。因此读者在使用滞后阻尼时需要慎重，从简单做起，逐步了解其功能后再应用于实践。

注意：不同的阻尼形式之间可以混合使用；对于不同的材料也可以按照初始条件的方式设置不同的阻尼形式和参数；动力分析中如果存在结构单元，需要采用SEL set damp命令，指定结构单元的阻尼，否则计算会提示出错。

11.6 网格尺寸的要求

输入波形的频率成分和土体的波速特性会影响波传播的数值精度。Kuhlemeyer和Lysmer（1973）的研究表明，要想精确描述模型中波的传播，那么网格的尺寸?l必须要小于输入波形最高频率对应的波长的1/8到1/10，也就是：

?11??l??~??

?810?

式中，?是最高频率对应的波长。

（11-21）

可见在动力计算中，土体的模量越小，即土体越软，最大网格尺寸越小，划分的网格数量越多。注意到，任何离散化的介质都存在能量传播的上限频率，只有当输入荷载的频率小于这个上限频率时，计算结果才有意义。因此上述公式不仅适用于FLAC3D，同样适用于其他基于时域的动力分析程序。

由于输入荷载的最大频率直接影响到单元的最大尺寸，所以对于脉冲荷载、爆炸荷载等这些频率范围很广的荷载形式，需要进行滤波处理，这部分内容在本章的11.7.1节中将会做介绍。

11.7 输入荷载的校正

这里的输入荷载一般指的是地震荷载，因为在地震反应分析中，常使用离散的荷载列表，因此在施加之前有必要进行滤波和基线校正。

11.7.1 滤波

滤波的目的是过滤掉原有波形中的高频分量，因为由式（11-21）可知，地震波的最大频率对网格尺寸的影响较大，最大频率越高，满足精度条件下的网格尺寸越小。采用滤波的方式，可以减小地震波的最大频率，从而增大计算所需的最小网格尺寸，减小单元数量，达到节约计算时间的目的。滤波可以通过OriginPro，SeismoSignal等软件进行，也可以使用FLAC3D提供的FFT.fis函数进行。

11.7.2 基线校正

在FLAC3D地震动力分析中，输入波通常为加速度时程，若将输入的加速度进行积分得到的最终速度和最终位移不为0，则在动力计算结束时模型底部会出现继续的速度和残余的位移，此时需要对加速度时程进行基线校正。即通过在原始加速度时程上增加一个低频率的波形（多项式或周期函数），使最终的速度和位移均为0。基线校正可以通过SeismoSignal软件进行，该软件提供了多种基线校正的方法，

建议读者使用。

11.8 动孔压模型与土体的液化

FLAC3D可以进行动力与渗流的耦合分析，能够模拟砂土在动力作用下的孔压积累直至土体的液化，FLAC3D采用了Finn模型来描述这种孔压积累的效应。Finn模型的实质是在Mohr-Coulomb模型的基础上增加了动孔压的上升模式，并假定动孔压的上升与塑形体积应变增量相关。

?时砂土的一维回弹模量为Er，则对于不排水条件下孔隙水压力的增量?u与塑性设在有效应力为?0体积应变增量??vd的关系为：

?u?Er??vd

（11-22）

FLAC3D提供了2种不同的塑形体积应变增量公式，包括Finn模式和Byrne模式。 1. Finn模式

Martin et al.（1975）[143]的试验表明，塑性体积应变与循环剪应变幅值之间的关系与固结压力无关。

为了实用目的，塑性体积应变增量??vd仅是总的累积体积应变?vd和剪应变?的函数：

??vd?C1(??C2?vd)?2C3?vd??C4?vd

（11-23）

其中，C1，C2，C3和C4为模型常数。对于相对密度为45%的结晶二氧化硅砂：C1 = 0.80，C2 = 0.79，C3 = 0.45，C4 = 0.73。

2. Byrne模式

Byrne（1991）提出了一种更简便的计算塑性体积应变增量的方法：

??vd?????C1exp??C2vd?

???0.4 C1（11-24）

其中，C1和C2为两个参数，大多数情况下两者存在如下的关系：

C2?（11-25）

参数C1与砂土的相对密度存在如下关系：

C1?7600（Dr)-2.5

（11-26）

另外，相对密度与标准贯入击数存在一定的经验关系：

Dr?15(N1)1/260

（11-27）

因此，参数C1也可以通过标准贯入击数来得到：

C1?8.7(N1)-1.2560

（11-28）

Byrne模型中还有一个参数C3表示剪应变阙值，即发生塑性体积应变的最小剪应变值。

由于Finn模型的基础是Mohr-Coulomb模型，因此Finn模型参数包含了Mohr-Coulomb的所有参数，

包括bulk，shear，cohesion，friction，tension，还包括用于计算孔压增量的参数。

? ff_switch表示孔压上升模式，0为Finn模式，1为Byrne模式。

? ff_c1，ff_c2，ff_c3，ff_c4分别为孔压上升模式公式中的C1，C2，C3和C4。 ? ff_latency表示两次应变反转之间的最小步数。

11.9 完全非线性动力耦合分析步骤

FLAC3D采用完全非线性的动力分析方法，可以考虑动力与渗流的耦合分析，模拟土体的液化。动力分析是在静力分析的基础上进行的，因此在动力计算之前要进行静力的力学计算和渗流计算，得到正确的应力场和渗流场。在动力计算前要考虑网格尺寸、边界条件、材料参数、阻尼类型、地震波调整等问题，一般动力分析过程见图11-16所示。

开始 静力分析

生成网格 否 网格尺寸满足要求？ 是 静力边界条件和初始条件 力学平衡 流体平衡 震前初始状态 选择材料模型及参数 动力分析

模型中刚度差异较大？ 是 设置动态多步计算 输入动力荷载 设置动力边界条件 选择阻尼类型 进行动力求解 地震波的调整 基线校正 滤波 结束 图11-16 完全非线性动力耦合问题求解流程图

11.10 应用实例——振动台液化试验模拟

下面用一个振动台液化试验的例子来说明动力分析的主要过程。在本节的分析中，比较了不同阻尼形式、流固耦合方式对计算结果的影响。读者在这个例子的基础上通过进一步细化改进，可以模拟其它条件下的振动台试验。

11.10.1 计算模型及参数

水平场地的液化模拟采用的计算模型见图11-17所示。为了减少计算时间，本例中没有采用真实的振动台尺寸，而采用了较大的模型（长30m，高11m）。其中，可液化砂土厚度为10m，其上部有一层厚度为1m的非液化层，垂直平面的方向取1m。计算中考虑了动力计算中的最大网格尺寸的要求，如表11-4所示，模型的网格划分见图11-18所示，共划分330个单元。为了便于不同计算工况的比较，在计算中对可液化层中的一个单元进行跟踪，分别记录该单元的孔压和平均有效应力随动荷载时间的变化曲线。

顶部非液化层可液化层监测单元3010单位：m1

图11-17 水平场地液化计算模型

图11-18 水平场地液化计算的网格划分 表11-4 水平场地液化计算中最大单元尺寸的计算 G (MPa) 30 cs?G/? f (Hz) 5 ? (m) 14.14 ?l (m/s) 70.71 (m) 1.4144 计算时，首先利用弹性模型使土体在重力作用下达到平衡，得到动力计算前的初始应力场，然后在模型底部及两侧的水平方向 (x方向) 施加正弦的速度边界，见图11-19所示。速度时程考虑了振动台试验激振的过程，从0 ~ 2 s速度逐渐增大达到最大值并稳定一定时间 (2 ~ 20s)，随后速度逐渐减小到0 (20

~ 30s)。计算中采用的速度时程曲线见图11-20所示。

计算工况为：

? 本构模型：分别采用FLAC3D自带的Finn模型；

? 阻尼形式：采用局部阻尼、瑞利阻尼和滞后阻尼，三种阻尼形式的参数见表111-5所示； ? 耦合形式：分别进行考虑动力‐渗流的耦合及不考虑耦合的计算。

表111-5 阻尼形式及参数选择

阻尼形式 局部阻尼 瑞利阻尼 滞后阻尼 阻尼参数 临界阻尼比D = 10% 阻尼比D = 10%，fmin = 3.67 Hz 参考应变 r = 6%

图11-19 水平场地液化计算采用的边界条件

0.10速度 / ms-10.05 0.00-0.05-0.100510152025303540时间 / s

图11-20 水平场地液化计算采用的速度荷载时程

11.10.2 计算过程

计算过程主要包括静力分析和动力分析，其中静力分析是动力分析的基础。 1. 静力分析

在进行动力分析之前要进行静力计算，获得振动施加前的初始应力状态。本例中初始应力状态的获得分为两步：

步骤1 设置弹性材料参数、流体模型参数、初始孔压场条件和边界条件，关闭流体模式和动力模式，并设置流体模量为0，在重力作用下达到静力平衡。

定义了一个设置水平应力的FISH函数，_k0是可以自行设定的侧压力系数，函数运行结束

后再次达到平衡，完成初始应力的计算。

步骤2

静力分析阶段的命令如例11.7。 例11.7：振动台模拟的静力分析

;振动台试验的例子 new

config dynamic fluid def model_dim h_R = 0 h_R1 = h_R + 1.0 end model_dim

gen zon bri p0 0 0 -10 p1 30 0 -10 p2 0 1 -10 p3 0 0 0 p4 30 1 -10 p5 0 1 0 p6 30 0 h_R p7 30 1 h_R size 30 1 10 group sand gen zon bri p0 0 0 0 p1 30 0 h_R p2 0 1 0 p3 0 0 1 p4 30 1 h_R p5 0 1 1 p6 30 0 h_R1 p7 30 1 h_R1 size 30 1 1 group top ;gen zon bri p0 0 0 -.5 p1 3 0 -.5 p2 0 1 -.5 p3 0 0 0 p4 3 1 -.5 p5 0 1 0 p6 3 0 0 p7 3 1 0 size 30 1 10 model elastic

prop bulk=3e7 shear=1e7 fric=35 ini dens 2000 model fl_iso

prop poro 0.5 perm 1e-8 ini fmod 2e8 fdens 1000

ini pp 0 grad 0 0 -10e3 ran z 0 -10.0 fix z ran z -9.9 -10.1 fix x ran x -.1 .1 fix x ran x 29.9 30.1 fix y set grav 10 set fluid off dyn off ini fmod 0 set mech rat 1e-6 solve def ini_conf _k0 = 1.0 pnt = zone_head loop while pnt # null

val = _k0 * z_szz(pnt) + (_k0-1.) * z_pp(pnt) z_sxx(pnt) = val z_syy(pnt) = val pnt=z_next(pnt) end_loop end ini_conf

solve save 11-8.sav

2. 动力分析

动力计算的命令文件见例11.8所示。在本例中对计算结果采取了分阶段保存的方法，定义了一个FISH函数，其功能是在动力计算完成特定的时刻即保存一个文件，这样便于计算结束后对中间结果的调用。

例11.8：振动台模拟的动力计算

rest 11-7.sav set dyn on fluid on ini fmod 2e8 set fluid pcut on model finn ran gro sand

prop bulk=3e7 shear=1e7 fric=35 ran gro sand ini dens 2000 ran gro sand

prop ff_latency=50 ff_switch = 0 ff_c1=0.8 ff_c2=0.79 ff_c3=0.45 ff_c4=0.73 ran gro sand ;扭剪试验结果 def setup freq=5.0 ampl=2

omega = 2.0 * pi * freq end setup

def sine_wave

vv = 9.36e-2 * sin(omega*dytime) if dytime < 2.0

sine_wave = dytime / 2.0 * vv else

if dytime < 20.0 sine_wave = vv else

if dytime <= 30.

sine_wave = (30.0 - dytime) / 10.0 * vv endif endif endif

if dytime > 30.0 sine_wave = 0.0 endif end free x

apply xvel = 1.0 hist sine_wave ran z -9.9 -10.1 apply xvel = 1.0 hist sine_wave ran x -.1 .1