1、第 11 章 非线性动力反应分析FLAC / FLAC3D 可以进行非线性动力反应分析，而且具有强大的动力分析功能。本章以 FLAC3D 为例，详细介绍了动力分析过程中的边界条件、阻尼形式、荷载要求等，并通过一些实例对个别问题做了详细解答。本章要点： FLAC 动力分析与等效线性方法的差别 动力分析时间步的确定方式及影响因素 动态多步的概念 动力荷载的形式及施加方法 动力边界条件的类型及适用条件 地震荷载输入的要点 三种阻尼形式的概念、参数确定及适用条件 网格尺寸的要求 输入荷载的校正 地震液化的模拟 完全非线性动力分析的步骤Equation Section 1111.1 概述FLAC / F
2、LAC3D 可以进行二维或三维的完全动力分析， FLAC/FLAC3D 中的动力分析功能是可选模块，需要在程序中添加动力分析模块才可以进行。FLAC3D 中在动力分析前需要采用以下的命令：CONFIG dynamic对于 FLAC，在程序开始时的 Model Options 对话框中选择 Dynamic 复选框。FLAC / FLAC3D 中的动力分析并不是只能孤立进行的，还可以与其他 FLAC/FLAC3D 元素进行耦合。（1）与结构单元相耦合，可以用来进行土与结构的动力相互作用。（2）与流体计算相耦合，可以模拟动力作用下土体孔隙水压力的上升直至土体液化。（3）与热力学计算相耦合，可以计算热
3、力荷载和动力荷载的共同作用。（3）采用大变形计算模式，可以分析岩土体在动力荷载作用下发生的大变形。FLAC 和 FLAC3D 可以模拟岩土体在外部（如地震）或内部（如风、爆炸、地铁振动）荷载作用下的完全非线性响应，因此可以适用于土动力学、岩石动力学等学科的计算。本章将以 FLAC3D 为例讨论动力计算的相关内容， FLAC 的动力分析可以参照执行。注意：FLAC 和 FLAC3D 的动力计算十分复杂，读者在阅读本章内容之前要对 FLAC3D 的静力计算、流体计算十分熟悉，具体可以参阅本书的第 7 章和第 12 章的内容。对于初次接触 FLAC3D 动力计算的读者，大多数都会提以下 2 个问题：
4、（1）FLAC 3D 动力分析与一般的等效线性方法有什么区别？（2）FLAC 3D 动力分析怎么会采用静力本构模型，比如 Mohr-Coulomb 模型？下面就这两个问题展开初步的讨论。11.1.1 与等效线性方法的关系在岩土地震工程中，等效线性方法广泛应用于计算地基土体中波的传播及土与结构的动力相互作用。该方法已被工程师、科研人员广泛接受。而 FLAC3D 采用的完全非线性方法没有获得广泛使用，因此需要对这两种方法之间的差异做简要介绍。1. 等效线性方法的特点等效线性方法的基本原理是，假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新
5、赋值。目前用于土动力分析的等效线性模型已有数种，根据骨干曲线的形状可以分为双直线模型、Ramberg-Osgood 模型、Hardin-Drnevich 模型等，其中又以 Hardin 模型使用最多。等效线性方法有如下的特点： 使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。 不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。 塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量
6、是应力张量的函数，称之为“流动法则” 。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。 大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。 本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。2. FLAC3D 非线性方法
7、的特点FLAC3D 采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。相对于等效线性方法而言，完全非线性分析方法主要有以下优点： 可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用 Rayleigh 阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。 采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。 由于采用了弹塑性模型，因此程序可以自动计算永久变形。 采用合理的塑性方程，使得塑性应变增量与应力相联系。 可以方便
8、地进行不同本构模型的比较。 可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。另外，FLAC 3D3.0 已将等效线性方法中的模量衰减曲线以阻尼的形式嵌入到程序当中（见本章 11.6.3节） ，使得 FLAC3D 的动力分析结果更易于被岩土地震工程师们所接受。11.1.2 FLAC3D 动力计算采用的本构模型FLAC3D 的动力计算可以采用任意的本构模型，比如弹性模型、Mohr-Coulomb 模型。这一点很多读者都不能接受，他们普遍认为 Mohr-Coulomb 是静力本构
9、模型，不适合用于动力分析，而应当采用更合适的 Hardin 模型。其实这是对 FLAC3D 动力计算的误解。 FLAC3D 的原理是求解动力方程，所以从其算法上来说，不管是进行静力分析还是动力分析，其实质都是求解运动方程。只是对于静力分析而言，采用了特定的阻尼方式以达到快速收敛的目的。所以，有的场合将 FLAC3D 的静力分析方法称为“拟动力方法” 。相应的，FLAC3D 在进行动力分析时，通过求解动力方程理所当然地可以得到合适的动力问题解答。对于本构模型的选择，主要是描述单元的应力-应变关系，如果是弹塑性的，则考虑的是单元的屈服准则、流动法则等。等效线性方法考虑土体的滞后性常常是通过将骨干曲
10、线进行变换，比如 Masing 二倍法，而在FLAC3D 的动力分析中，滞后性是通过阻尼来考虑，通过设置合适的阻尼形式和阻尼参数，同样可以描述土体在动力作用下的滞回曲线和滞回圈。因此，FLAC 3D 动力分析中采用的本构模型可以选取任意模型，其参数也是对应静力本构模型的参数，关键是要设置合适的阻尼形式、阻尼参数、边界条件等，这些内容将在本章的后续内容中进行讲解。11.2 动力时间步动力计算中临界计算时间步的计算如下：（11-1）maxincrit fpVCA其中， 为 p 波波速，与材料的体积模量 K 和剪切模量 G 有关，可以表示为：C（11-2）4/3pC为四面体子单元（sub-zone
11、）的体积， 为与四面体子单元相关的最大表面积， 表示遍历VmaxfA min所有的单元，包括结构单元和接触面单元。由于式（11-1）只是临界时间步的一个估计值，因此在使用中采用了一个安全系数，乘以 0.5。因此，当采用无刚度比例的阻尼时，动力分析的时间步为：（11-3）2dcritt如果采用了刚度比例的阻尼，那么为了保持数值稳定性，时间步必须减小。Belytschko（1983）提出了一个临界时间步的公式 ，其中考虑了刚度比例阻尼的影响。t（11-4）2max1t其中， 为系统的最高特征频率， 为该频率下的临界阻尼比。max注意：FLAC 3D 在动力计算中，程序会根据数值计算的稳定性自动设置
12、动力计算时间步，一般不建议读者对这个默认的时间步进行放大。甚至，在大应变计算过程中，如果出现很大的网格变形并导致网格的几何错误时，还要对默认的时间步进行折减，降低动力时间步，以达到数值稳定的目的。11.3 动态多步由式（11-1）可知，FLAC 3D 动力计算中时间步需要遍历所有单元，取所有单元临界时间步中的最小值，因此时间步是由几何尺寸较小、模量较大的单元来确定的。因此，在计算中，尤其是在试算期间，要尽量避免较小的单元尺寸及刚度很大的材料，比如用实体单元来模拟较薄的混凝土墙，这种情况下必然会使动力时间步非常小，从而造成计算时间过长。可以通过采用结构单元或暂时不考虑混凝土墙的办法来进行试算，等
13、到有关参数调试完成后再进行细化计算。当计算模型中存在刚度差异较大、模型网格尺寸不均匀的情况时，FLAC 3D 可以采用“动态多步”（Dynamic Multi-stepping）的过程来有效减少计算所需要的时间。在此过程中，模型单元和节点按照相近最大时步进行分组和排序，然后每个组在特定的时步下运行，信息在适当的时候在单元之间进行交换。动态多步的调用采用如下命令：SET dyn multi on下面用一个简单的例子来描述动态多步的应用效果。同时读者可以从例子中了解到利用 FISH 函数来编写简单的动力荷载的方法。1. 问题描述如 图 11-1 所示，土体的深度为 10 m，挡土墙的高度为 5 m
14、，两者的模量差异为 20 倍。动力荷载从模型底部输入，主要分析目的是了解动态多步对计算时间的影响。图 11-1 动态多步作用的实例本例计算中不考虑重力的影响，因此不用进行初始应力设置和平衡，直接进行动力计算。动力荷载采用正弦函数，采用 FISH 函数的方法进行定义，可以方便修改荷载的频率（freq） 、幅值等。2. 命令流例 11.1：动态多步实例newconf dyn ;打开动力计算功能gen zone brick size 10 5 10mod elasmod null range x=0,5 z=5,10 ;删除部分网格fix z range x=-.1 .1 z=.1 10.1 ;设置
15、静力边界条件fix z range x=9.9,10.1 z=.1 10.1fix y range y=-.1 .1fix y range y=4.9 5.1prop bulk 2e8 shear 1e8 ;设置土体参数prop bulk 4e9 shear 2e9 range x=5,6 z=5,10 ;设置墙体参数（土体参数的 20 倍）ini dens 2000 ;设置密度def setup ;动荷载中的变量赋值freq = 1.0omega = 2.0 * pi * freqold_time = clock1end1 clock 是 FISH 变量，表示计算机目前的时间土体K = 20
16、0 MPaG = 100 MPa墙体K = 4000 MPaG = 2000 MPa5 m10 msetup ;执行变量赋值def wave ;定义动荷载函数wave = sin(omega * dytime) ;定义动荷载变量endapply xvel = 1 hist wave range z=-.1 .1 ;施加动荷载apply zvel = 0 range z=-.1 .1hist gp xvel 5,2,0hist gp xvel 5,2,10hist gp zvel 5,2,10hist dytimedef tim ;估算程序运行的时间tim = 0.01 * (clock - o
17、ld_time）endset dyn multi on ;设置动态多步solve age 1.0print tim ;输出计算时间print dyn ;输出动力计算相关信息save mult1.sav注意：动力计算中必须设置材料的密度，若模型中存在结构单元，也必须设置结构单元的密度，否则会出错。采用 FISH 函数定义动力荷载时， FISH 函数和变量应具有相同的名称。因为设置动力边界条件命令中的 hist 关键词后面必须要跟随一个 FISH 函数名，FISH 变量需要调用 FLAC3D 中的内置标量 dytime，该变量是动力计算的真实时间，通过调用可以给函数提供预订变化的数值。所以，一般
18、FISH 定义动荷载的方法如下（以定义函数 xxx 为例）：def xxxxxx = dytimeendapp xvel = 1.0 hist xxx range 3. 计算过程与输出结果计算过程中命令窗口会提示动力计算的步数、动力时间和时间步，计算结束后可以将模型底部和墙体顶部节点的水平速度时程输出，使用以下的命令：plot hist 1 2 v 4输出结果见 图 11-2 所示。读者可以采用设置动态多步和不设置动态多步两种情况分别进行计算，观察 图 11-2 中的速度时程曲线以及 FLAC3D 命令窗口中输出结果（ 图 11-3 和 图 11-4） 。可以发现，时程曲线、动力计算的时间步、
19、迭代步数均一致，不同的是设置动态多步的情况下花费的计算时间较少。注意：程序运行的时间将由于计算机配置的不同而存在差别。从 图 11-3 中的输出信息可以看出，动态多步将模型中的 375 个单元分成了两类，并提供了不同的时步乘子，其中只有小部分单元（墙体单元 65 个）的时步乘子为 1，其他大部分单元（土体单元 310 个）拥有较高的时步乘子，这样可以大大加快计算的进度。FLAC3D 3.0Itasca Consulting Group, Inc.Mineapolis, MN USAStep 821616:46:10 Wed Apr 02 08History2.0 4.0 6.0 8.0x10-
20、1-1.0-0.8-0.6-0.4-0.20.00.20.40.60.81.01.21 X-Velocity Gp 5 Linestyle -1.00e+00 1.00e+00 2 X-Velocity Gp 68Linestyle -1.92e+00 1.340e+00 Vs.4 Dynamic Time 1.217e-03 9.93e-01图 11-2 动力计算结束时模型底部和顶部的水平速度时程曲线图 11-3 设置动态多步情况下的输出信息图 11-4 未设置动态多步情况下的输出信息11.4 动力荷载和边界条件利用 FLAC3D 进行动力计算时，有以下 3 个问题需要读者认真考虑： 动力荷
21、载和边界条件； 力学阻尼； 模型中波的传播。本节及后续的两节将分别针对以上三个问题展开讨论。11.4.1 动力荷载的类型与施加方法FLAC3D 可以在模型边界或内部节点施加动荷载来模拟材料受到外部或内部动力作用下的反应，程序允许的动力荷载输入可以是：（1）加速度时程， （2）速度时程， （3）应力（压力）时程， （4）集中力时程。动力荷载的施加采用 APPLY 命令，另外，采用 APPLY Interior 命令可以将加速度、速度和力的时程施加到模型内部的节点上。动力荷载的形式主要有 2 种： FISH 函数。FISH 函数表达的动力荷载往往比较规则，也常用于试算阶段的动力输入，因为试算时可以
22、不用过多考虑荷载的频率、基线校正等问题。本章例 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
23、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 数据对的
24、表第 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 动力边界条件设置时，使
25、用 APPLY 命令调用已读入的表，下面的命令以施加水平向速度荷载为例。app xvel 1.0 hist table 1 range 命令中的 1.0 表示表格 y 向数据的乘子，可以方便地控制荷载幅值的大小。另外，动荷载的输入可以沿着 x，y，z 的任意方向施加或者模型边界的法向和切向施加。注意：特定的边界条件不能在相同的边界上进行混合施加，否则程序会提示出错。11.4.2 边界条件的设置在动力问题中，模型周围边界条件的选取是一个主要内容，因为边界上会存在波的反射，对动力分析的结果产生影响。把分析模型的范围设置得越大，分析结果就越好，但较大的模型会导致巨大的计算负担。FLAC 3D 中提供
26、了静止（粘性）边界和自由场边界两种边界条件来减少模型边界上的波的反射。1. 静态边界FLAC3D 中允许采用静态边界（也称粘性边界，吸收边界）条件来吸收边界上的入射波。FLAC 3D 中的静态边界是 Lysmer 和 Kuhlemeyer（1969）提出的，具体做法是在模型的法向和切向分别设置自由的阻尼器从而实现吸收入射波的目的，阻尼器提供的法向和切向粘性力分别为式（11-5）和（11-6 ） 。（11-5）npntCv（11-6）ss其中， ， 分别为模型边界上法向和切向的速度分量， 为介质密度， ， 分别为 p 波和 s 波nvs pCs的波速。这种静态边界对于入射角大于 30 度的入射波
27、基本能够完全吸收，对于入射角较小的波，比如面波，虽然仍有一定的吸收能力，但吸收不完全。静态边界可以加在整体坐标系上，也可以加载在倾斜边界的法向和切向上。如果在倾斜边界的法向和切向施加静态边界，则需要同时使用 nquiet，dquiet，squiet 条件。整体坐标系的静态边界条件设置使用命令为：APPLY xquiet (yquiet, zquiet) range 使用倾斜边界上的静态边界条件命令为：APPLY nquiet dquiet squiet range 注意：（1）施加动力边界条件后，这些边界上原先的静力边界条件将被自动去掉（free） ，在动力荷载施加期间，程序始终自动计算边界上
28、的作用力，用户不能将静态边界条件去掉；也不能在静态边界上施加加速度、速度边界条件，因为静态边界上的作用力是根据边界上的速度分量计算得到的，如果再施加速度荷载就会使静态边界失效。若需要在静态边界上输入动荷载，则只能输入应力时程。可以将加速度、速度时程通过转换公式（11-7）和（11-8）形成应力时程施加到静态边界上。（11-7）=-2()npnCv（11-8）ss式中的 ， 分别为施加在静态边界上的法向应力和切向应力，公式中的系数 2 表示施加的能ns量中只有一半是向上传播作为动力输入的，另一半向边界下部传播。公式中的负号是为了使应力施加后节点的速度能与实际一致。（2）对于动力荷载来源于模型内部
29、（如隧道中的列车振动问题）的情况，可以将动力荷载直接施加在节点上，这种情况下使用静态边界可以有效减小人工边界上的反射，并且不需要再施加下面提到的自由场边界。（3）动力计算过程中应避免静力荷载的变化。比如，在一个已经施加底部静态边界的模型上进行开挖，会造成整个模型向上移动。因为施加静态边界时程序自动计算了施加在模型底部边界上的反力，这些边界反力不能与开挖后的模型相平衡，会引起模型的整体上移。下面给出一个简单的静态边界的例子。问题描述：如 图 11-5 所示，一根竖直的弹性杆，高 50m，宽 1m，体积模量和剪切模量分别为20MPa 和 10MPa，材料密度为 1000kg/m3。在杆的底部边界的
30、两个水平方向设置静态边界条件，而杆的顶部为自由表面，在杆的底部施加水平方向的应力冲击荷载。根据材料参数，可得计算得到剪切波速为 100 m/s， 的乘积为 105。应力冲击荷载的幅值设置为 2105，根据公式（11-7）可以得到等效sCsC速度荷载幅值为 1m/s。计算中不考虑重力的作用，因此不需要进行初始应力的计算。分析杆的底部、中心点、顶部三个典型节点的速度响应，计算结果如 图 11-6 所示。可以发现模型底部产生的速度荷载幅值就是 1m/s，图中最后两个脉冲是从模型顶部自由面反射回来的波。在模型顶部节点的速度幅值是底部输入幅值的两倍，而且顶部反射回来的波传到模型底部以后，模型基本结束了振
31、动。这些现象说明采用静态边界的设置是合理有效的。计算文件如下：例 11.2：静态边界的例子newconfig dyngen zone brick size 1,1,50model elasprop shear 1e7 bulk 2e7ini dens 1000def setupomega = 2.0 * pi * freqpulse = 1.0 / freqendset freq=4.0setupdef waveif dytime pulsewave = 0.0elsewave = 0.5 * (1.0 - cos(omega * dytime)endifendrange name botto
32、m z=-.1 .1fix z range z=.5 55 ;将上部网格都施加数值向约束apply dquiet squiet range bottomapply sxz -2e5 hist wave syz 0.0 szz 0.0 range bottom ;-2e5 的系数来源于 的值sCapply nvel 0 plane norm 0,0,1 range bottomhist gp xvel 0,0,0hist gp xvel 0,0,25hist gp xvel 0,0,50hist dytimehist waveplot create hhhplot add hist 1 2 3
33、vs 4 plot showsolve age 250 m1 m图 11-5 静态边界条件示例FLAC3D 3.00Itasca Consulting Group, Inc.Mineapolis, MN USAStep 147915:13:54 Thu Apr 03 208History0.5 1.0 1.5 0.00.20.40.60.81.01.21.41.61.81 X-Velocity Gp 1 Linestyle -3.836e-02 1.01e+002 X-Velocity Gp 101 Linestyle -3.348e-02 1.00e+003 X-Velocity Gp 20
34、1 Linestyle -5.125e-02 1.95e+00 Vs.4 Dynamic Time 1.217e-02 1.789e+00图 11-6 在静态边界上输入应力荷载得到的响应2. 自由场边界对诸如大坝之类的地面结构进行动力反应分析时，在模型各侧面的边界条件须考虑为没有地面结构时的自由场运动。FLAC 3D 通过在模型四周生成二维和一维网格的方法来实现这种自由场边界条件（见图 11-7 所示） ，主体网格的侧边界通过阻尼器与自由场网格进行耦合，自由场网格的不平衡力施加到主体网格的边界上。由于自由场边界提供了与无限场地相同的效果，因此向上的面波在边界上不会产生扭曲。图 11-7 自由场
35、边界示意图对于 FLAC3D 而言，自由场边界模型包括 4 个平面网格和 4 个柱体网格，平面网格在模型边界上与主体网格是一一对应的，柱体网格相当于平面自由场网格的自由场边界。其中，平面自由场网格是二维计算，假设在面的法向无限延伸；柱体自由场网格是一维计算，假设在柱体两端无限延伸。下面给出一个简单的自由场边界条件的例子，见例 11.3。底部 中心顶部例 11.3：自由场边界条件的实例new;第一步：静力计算阶段config dynset dyn offgen zone brick size 6 3 2gen zone brick size 2 3 2 p0 0 0 2gen zone bric
36、k size 2 3 2 p0 4 0 2gen zone wedge size 1 3 2 p0 2 0 2gen zone wedge size 1 3 2 p0 4 3 2 p1 3 3 2 p2 4 0 2 p3 4 3 4 第二步：动力计算阶段set dyn ondef iniwaveper = 0.01 endiniwavedef wavewave = 0.5 * (1.0 - cos (2*pi*dytime/per)endfree x y z ran z -0.1 0.1 ;去掉模型底部原有的静力条件apply nquiet squiet dquiet ran z -0.1 0
37、.1 ;静态边界条件apply dstress 1.0 hist wave ran z -0.1 0.1 ;加动力荷载apply ff ;施加自由场边界条件group ff_cornergroup ff_side ran x 0 6group ff_side ran y 0 3group main_grid ran x 0 6 y 0 3set dyn time = 0 ;设置动力计算从 0s 开始hist reset ;清空已有的历史信息hist unbalhist dytime; 主体网格hist gp xvel 2 1 0hist gp xvel 2 1 5.0; 柱体网格hist gp
38、 xvel -1 -1 0hist gp xvel -1 -1 5.0; 平行于 y 方向的二维自由场网格hist gp xvel -1 0 0hist gp xvel -1 0 5.0; 平行于 x 方向的二维自由场网格hist gp xvel 2 -1 0hist gp xvel 2 -1 5.0solve age 0.015save 11-3_2.sav计算分两步，第一步进行重力作用下的初始应力计算，这里采用的是直接重力加载的方法生成初始应力，第二步进行动力计算。在设置动力边界条件时采用了 4 个步骤：去掉模型底部已有的静力约束条件； 步 骤 1施加静态边界条件； 步 骤 2施加动力荷载
39、； 步 骤 3施加自由场边界条件。 步 骤 4施加自由场边界条件以后，程序自动在主体网格周围生成一圈自由场网格。当计算模型的单元数较多时，这个生成过程需要花费较多的时间。可以对生成的自由场网格定义 group，如例 11.3 中就将二维自由场网格和一维自由场网格分别赋值了 group，这样便于模型的观察和主体网格结果的后处理。图 11-8 自由场边界施加前后的网格计算中对主体网格、柱体自由场网格和平面自由场网格对应位置的水平向速度进行了监测，计算结果见 图 11-9 所示。可以发现，主体网格与周围的自由场网格同步运动，达到了自由场网格的目的。FLAC3D 3.0Itasca Consultin
40、g Group, Inc.Mineapolis, MN USAStep 9217:04:15 Thu Apr 03 208History0.2 0.4 0.6 0.8 1.0 1.2 1.4x10-21.02.03.04.05.06.07.08.09.0x10-23 X-Velocity Gp 10 Linestyle 6.71e-06 9.06e-02 5 X-Velocity Gp 368 Linestyle6.71e-06 9.049e-02 7 X-Velocity Gp 145 Linestyle 6.71e-06 9.049e-02 9 X-Velocity Gp 23Linest
41、yle 6.71e-06 9.063e-02 Vs. 2 Dynamic Time5.735e-05 1.497e-02FLAC3D 3.0Itasca Consulting Group, Inc.Mineapolis, MN USAStep 9217:04:5 Thu Apr 03 208History0.2 0.4 0.6 0.8 1.0 1.2 1.4x10-2 0.0.20.40.60.81.0x10-14 X-Velocity Gp 102 Linestyle -6.480e-04 1.05e-01 6 X-Velocity Gp 380 Linestyle-3.849e-05 1.
42、01e-01 8 X-Velocity Gp 176 Linestyle -3.849e-05 1.01e-01 10 X-Velocity Gp 274Linestyle -6.91e-04 1.07e-01 Vs. 2 Dynamic Time5.735e-05 1.497e-02模型底部 模型顶部图 11-9 自由场边界条件的示例结果讨论：读者可以修改例 11.3 中边界条件设置的命令，比如去掉某个命令或者调换命令的前后位置，可以得到以下结论： free 命令在本例中可以删去，因为接下来设置静态边界本身就会将模型底部边界上已有的静力条件设置为 free。但是如果静力计算阶段模型底部的边界
43、条件设置为 fix x y z，而且施加的动力荷载形式为加速度或者速度，那么就需要首先将这些 fix 速度的静力条件设置为 free，否则会出现动力加载的错误。 本例中 app ff 的位置对结果没有影响，但是在动力计算中仍然要把 app ff 放在所有边界条件的后面，换句话说，静态边界条件、动力荷载施加等必须在自由场边界条件设置之前完成。另外，使用自由场边界条件还需要注意以下几点： 如果动力源只存在于模型的内部，那么可以不必设置自由场边界。 自由场边界设置对模型的形状有一定的要求：模型底部水平，重力方向为 z 向，四个侧面垂直，法向分别为 x、y 向。如果地震波的传播方向不是竖直方向的，则应
44、当进行坐标轴旋转，使得 z轴与地震波传播方向相同。这种情况下，重力方向将与 z 轴存在一定的夹角，而且模型边界也与水平面产生倾斜。 其他边界条件应该在 APPLY ff 之前进行设置。 APPLY ff 将边界上单元的属性、条件和变量全部转移至 ff 单元上，并且设置以后主体网格上的改动将不会被 FF 边界所响应，这一点与静态边界类似。 使用自由场边界对主体网格的本构模型没有要求，并可以进行竖向流体计算相耦合。 自由场边界进行的是小变形计算，主体网格可进行大变形计算，因此自由场边界上的变形要相对较小，如果自由场网格上的变形较大，那么需要扩大人工边界条件的选取。 存在 attach 的边界将不能
45、设置 FF 边界，而且主体网格边界上的 Interface 将不能连续到自由场网格。 可以对生成的自由场网格进行 group 赋值，这样在后处理的云图显示时去掉自由场网格，生成只有主体网格的图形。11.4.3 地震荷载的输入地震反应分析是 FLAC3D 动力计算的主要应用领域，所以有必要对地震荷载的输入做单独介绍。地震荷载常用表的文本文件格式读入到 FLAC3D，再施加到模型底部的边界上。对于地震荷载的动力边界条件选择可以参考以下标准：1. 刚性地基如果模型底部为岩石等模量较大的材料，可以在底部直接施加加速度或速度荷载，并采用自由场边界条件，模型底部无需施加静态边界条件。2. 柔性地基如果模型
46、底部的单元为土体，尤其是软土，则不能直接施加加速度和速度，而需要将加速度、速度转换成应力时程，再施加到模型底部。模型周围采用自由场边界条件，模型底部可采用静态边界条件（nquiet ，dquiet，squiet） 。11.5 力学阻尼阻尼的产生主要来源于材料的内部摩擦以及可能存在的接触表面的滑动。FLAC 3D 采用求解动力方程的方法解决两类力学问题：准静力问题和动力问题。在这两类问题中都要使用阻尼，但准静力问题需要更多的阻尼使得动力方程能够更快的收敛平衡。对于动力问题中的阻尼，需要在数值模拟中重现自然系统在动荷载作用下的阻尼大小。目前，FLAC 3D 动力计算提供了三种阻尼形式供用户选择，分
47、别是瑞利阻尼、局部阻尼和滞后阻尼。11.5.1 瑞利阻尼（Rayleigh damping）瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵 C 与刚度矩阵 K 和质量矩阵 M 有关：（11-9）其中， 为与质量成比例的阻尼常数， 为与刚度成比例的阻尼常数。瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。 图 11-10 为归一化的临界阻尼比与角频率之间的关系曲线，其中包含了仅设置刚度
48、分量、仅设置质量分量以及二者叠加的三种结果。可以看出，采用叠加的方法得到的阻尼比在较大的频率范围内保持定值，因此瑞利阻尼可以近似地反映岩土体具有的频率无关性。图 11-10 归一化的临界阻尼比与角频率之间的关系1. 瑞利阻尼的两个参数设置瑞利阻尼的命令为：SET dyn damp rayleigh 参数 1 参数 2根据 图 11-10 叠加结果的阻尼比曲线最小值可以确定瑞利阻尼的两个参数，分别是最小临界阻尼比（参数 1）和最小中心频率 （参数 2） ，可以按照式（11-10）进行计算：minmin（11-10）1/2min()注意：参数 2 中的单位是 Hz，而不是角频率的单位。2. 瑞利阻尼参数选择的原则很多读者最十分关心瑞利阻尼参数的选择。对于岩土材料而言，临界阻尼比的范围一般是 25%，而结构系统的临界阻尼比为 210%。在使用弹塑性模型进行动力计算时，相当多的能量消散于材料发生的塑性流动阶段，因此在进行大应变的动力分析时，只需要设置一个很小的阻尼比（比如 0.5%）就能满足要求，而且达到塑形以后，随着应力-应变滞回圈的扩大，能量的消散也逐渐明显。（1）中心频率的确定瑞利阻尼是与频率相关的，但是在一定的频率范围内，瑞利阻尼基本
