(完整word版)FLAC动力分析

当计算模型中存在刚度差异较大、模型网格尺寸不均匀的情况时，FLAC3D可以采用“动态多步”（Dynamic Multi-stepping）的过程来有效减少计算所需要的时间。在此过程中，模型单元和节点按照相近最大时步进行分组和排序，然后每个组在特定的时步下运行，信息在适当的时候在单元之间进行交换。
动态多步的调用采用如下命令：
（1）FLAC3D动力分析与一般的等效线性方法有什么区别？
（2）FLAC3D动力分析怎么会采用静力本构模型，比如Mohr-Coulomb模型？
下面就这两个问题展开初步的讨论。
11.
在岩土地震工程中，等效线性方法广泛应用于计算地基土体中波的传播及土与结构的动力相互作用。该方法已被工程师、科研人员广泛接受。而FLAC3D采用的完全非线性方法没有获得广泛使用，因此需要对这两种方法之间的差异做简要介绍。
old_time = clock
end
setup ;执行变量赋值
def wave ;定义动荷载函数
wave = sin(omega * dytime);定义动荷载变量
end
apply xvel = 1 hist wave range z=-.1 .1;施加动荷载
apply zvel = 0 range z=-.1 .1
本章将以FLAC3D为例讨论动力计算的相关内容，FLAC的动力分析可以参照执行。
注意：FLAC和FLAC3D的动力计算十分复杂，读者在阅读本章内容之前要对FLAC3D的静力计算、流体计算十分熟悉，具体可以参阅本书的第7章和第12章的内容。
对于初次接触FLAC3D动力计算的读者，大多数都会提以下2个问题：
其实这是对FLAC3D动力计算的误解。FLAC3D的原理是求解动力方程，所以从其算法上来说，不管是进行静力分析还是动力分析，其实质都是求解运动方程。只是对于静力分析而言，采用了特定的阻尼方式以达到快速收敛的目的。所以，有的场合将FLAC3D的静力分析方法称为“拟动力方法”。相应的，FLAC3D在进行动力分析时，通过求解动力方程理所当然地可以得到合适的动力问题解答。对于本构模型的选择，主要是描述单元的应力-应变关系，如果是弹塑性的，则考虑的是单元的屈服准则、流动法则等。
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
使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。
不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。
1.等效线性方法的特点
等效线性方法的基本原理是，假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。目前用于土动力分析的等效线性模型已有数种，根据骨干曲线的形状可以分为双直线模型、Ramberg-Osgood模型、Hardin-Drnevich模型等，其中又以Hardin模型使用最多。等效线性方法有如下的特点：
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
（2）与流体计算相耦合，可以模拟动力作用下土体孔隙水压力的上升直至土体液化。
（3）与热力学计算相耦合，可以计算热力荷载和动力荷载的共同作用。
（3）采用大变形计算模式，可以分析岩土体在动力荷载作用下发生的大变形。
FLAC和FLAC3D可以模拟岩土体在外部（如地震）或内部（如风、爆炸、地铁振动）荷载作用下的完全非线性响应，因此可以适用于土动力学、岩石动力学等学科的计算。
（11-3）
如果采用了刚度比例的阻尼，那么为了保持数值稳定性，时间步必须减小。Belytschko（1983）提出了一个临界时间步的公式 ，其中考虑了刚度比例阻尼的影响。
（11-4）
其中， 为系统的最高特征频率， 为该频率下的临界阻尼比。
注意：FLAC3D在动力计算中，程序会根据数值计算的稳定性自动设置动力计算时间步，一般不建议读者对这个默认的时间步进行放大。甚至，在大应变计算过程中，如果出现很大的网格变形并导致网格的几何错误时，还要对默认的时间步进行折减，降低动力时间步，以达到数值稳定的目的。
塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。
大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。
等效线性方法考虑土体的滞后性常常是通过将骨干曲线进行变换，比如Masing二倍法，而在FLAC3D的动力分析中，滞后性是通过阻尼来考虑，通过设置合适的阻尼形式和阻尼参数，同样可以描述土体在动力作用下的滞回曲线和滞回圈。
因此，FLAC3D动力分析中采用的本构模型可以选取任意模型，其参数也是对应静力本构模型的参数，关键是要设置合适的阻尼形式、阻尼参数、边界条件等，这些内容将在本章的后续内容中进行讲解。
print tim;输出计算时间
print dyn ;输出动力计算相关信息
save mult1.sav
注意：动力计算中必须设置材料的密度，若模型中存在结构单元，也必须设置结构单元的密度，否则会出错。采用FISH函数定义动力荷载时，FISH函数和变量应具有相同的名称。因为设置动力边界条件命令中的hist关键词后面必须要跟随一个FISH函数名，FISH变量需要调用FLAC3D中的内置标量dytime，该变量是动力计算的真实时间，通过调用可以给函数提供预订变化的数值。所以，一般FISH定义动荷载的方法如下（以定义函数xxx为例）：
11.
动力计算中临界计算时间步的计算如下：
（11-1）
其中， 为p波波速，与材料的体积模量K和剪切模量G有关，可以表示为：
（11-2）
为四面体子单元（sub-zone）的体积， 为与四面体子单元相关的最大表面积， 表示遍历所有的单元，包括结构单元和接触面单元。
由于式（11-1）只是临界时间步的一个估计值，因此在使用中采用了一个安全系数，乘以0.5。因此，当采用无刚度比例的阻尼时，动力分析的时间步为：
本章要点：
FLAC动力分析与等效线性方法的差别
动力分析时间步的确定方式及影响因素
动态多步的概念
动力荷载的形式及施加方法
动力边界条件的类型及适用条件
地震荷载输入的要点
三种阻尼形式的概念、参数确定及适用条件
网格尺寸的要求
输入荷载的校正
地震液化的模拟
完全非线性动力分析的步骤
11.
FLAC/FLAC3D可以进行二维或三维的完全动力分析，FLAC/FLAC3D中的动力分析功能是可选模块，需要在程序中添加动力分析模块才可以进行。
本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。
采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。
由于采用了弹塑性模型，因此程序可以自动计算永久变形。
采用合理的塑性方程，使得塑性应变增量与应力相联系。
可以方便地进行不同本构模型的比较。
可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。
SET dyn multi on
下面用一个简单的例子来描述动态多步的应用效果。同时读者可以从例子中了解到利用FISH函数来编写简单的动力荷载的方法。
1.问题描述
如图111所示，土体的深度为10 m，挡土墙的高度为5 m，两者的模量差异为20倍。动力荷载从模型底部输入，主要分析目的是了解动态多步对计算时间的影响。
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
FLAC3D中在动力分析前需要采用以下的命令：
CONFIG dynamic
对于FLAC，在程序开始时的Model Options对话框中选择Dynamic复选框。
FLAC /FLAC3D中的动力分析并不是只能孤立进行的，还可以与其他FLAC/FLAC3D元素进行耦合。
（1）与结构单元相耦合，可以用来进行土与结构的动力相互作用。
def xxx
xxx =…dytime
end
app xvel = 1.0 hist xxx range…
3.计算过程与输出结果
计算过程中命令窗口会提示动力计算的步数、动力时间和时间步，计算结束后可以将模型底部和墙体顶部节点的水平速度时程输出，使用以下的命令：