[电力/水利]FLAC动力分析

[电力/水利]FLAC动力分析 第11章 非线性动力反应分析 3D3DFLAC / FLAC可以进行非线性动力反应分析，而且具有强大的动力分析功能。本章以FLAC为 例，详细介绍了动力分析过程中的边界条件、阻尼形式、荷载要求等，并通过一些实例对个别问题做了 详细解答。 本章要点: , FLAC动力分析与等效线性方法的差别 , 动力分析时间步的确定方式及影响因素 , 动态多步的概念 , 动力荷载的形式及施加方法 , 动力边界条件的类型及适用条件 , 地震荷载输入的要点 , 三种阻尼形式的概念、参数确定及适用条件 , 网格尺寸的要求 , 输入荷载的校正 , 地震液化的模拟 , 完全非线性动力分析的步骤 Equation Section 11 11.1 概述 3D3DFLAC / FLAC可以进行二维或三维的完全动力分析，FLAC/FLAC中的动力分析功能是可选模 块，需要在程序中添加动力分析模块才可以进行。 3DFLAC中在动力分析前需要采用以下的命令: CONFIG dynamic 对于FLAC，在程序开始时的Model Options对话框中选择Dynamic复选框。 3D3DFLAC / FLAC中的动力分析并不是只能孤立进行的，还可以与其他FLAC/FLAC元素进行耦合。 (1)与结构单元相耦合，可以用来进行土与结构的动力相互作用。 (2)与流体计算相耦合，可以模拟动力作用下土体孔隙水压力的上升直至土体液化。 (3)与热力学计算相耦合，可以计算热力荷载和动力荷载的共同作用。 (3)采用大变形计算模式，可以分析岩土体在动力荷载作用下发生的大变形。 3DFLAC和FLAC可以模拟岩土体在外部(如地震)或内部(如风、爆炸、地铁振动)荷载作用下 的完全非线性响应，因此可以适用于土动力学、岩石动力学等学科的计算。 3D本章将以FLAC为例讨论动力计算的相关内容，FLAC的动力分析可以参照执行。 3D3D注意:FLAC和FLAC的动力计算十分复杂～读者在阅读本章内容之前要对FLAC的静力计算、流体计算十分熟悉～具体可以参阅本书的第7章和第12章的内容。 3D对于初次接触FLAC动力计算的读者，大多数都会提以下2个问题: 3D(1)FLAC动力分析与一般的等效线性方法有什么区别, 3D(2)FLAC动力分析怎么会采用静力本构模型，比如Mohr-Coulomb模型, 下面就这两个问题展开初步的讨论。 11.1.1 与等效线性方法的关系 在岩土地震工程中，等效线性方法广泛应用于计算地基土体中波的传播及土与结构的动力相互作用。 3D该方法已被工程师、科研人员广泛接受。而FLAC采用的完全非线性方法没有获得广泛使用，因此需要对这两种方法之间的差异做简要介绍。 1. 等效线性方法的特点 等效线性方法的基本原理是，假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。目前用于土动力分析的等效线性模型已有数种，根据骨干曲线的形状可以分为双直线模型、Ramberg-Osgood模型、Hardin-Drnevich模型等，其中又以Hardin模型使用最多。等效线性方法有如下的特点: , 使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小;在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。 , 不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。 , 塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量(而不是应变增量张量)是应力张量的函数。因此，塑性屈服的模拟不合理。 , 大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。 , 本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。 3D2. FLAC非线性方法的特点 3DFLAC采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。相对于等效线性方法而言，完全非线性分析方法主要有以下优点: , 可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部(local)阻尼，则在动力计算中阻尼参数将保持不变。 , 采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。 , 由于采用了弹塑性模型，因此程序可以自动计算永久变形。 , 采用合理的塑性方程，使得塑性应变增量与应力相联系。 , 可以方便地进行不同本构模型的比较。 , 可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。 3D另外，FLAC3.0已将等效线性方法中的模量衰减曲线以阻尼的形式嵌入到程序当中(见本章11.6.3 3D节)，使得FLAC的动力分析结果更易于被岩土地震工程师们所接受。 3D11.1.2 FLAC动力计算采用的本构模型 3DFLAC的动力计算可以采用任意的本构模型，比如弹性模型、Mohr-Coulomb模型。这一点很多读者都不能接受，他们普遍认为Mohr-Coulomb是静力本构模型，不适合用于动力分析，而应当采用更合适的Hardin模型。 3D3D其实这是对FLAC动力计算的误解。FLAC的原理是求解动力方程，所以从其算法上来说，不管是进行静力分析还是动力分析，其实质都是求解运动方程。只是对于静力分析而言，采用了特定的阻尼 3D方式以达到快速收敛的目的。所以，有的场合将FLAC的静力分析方法称为“拟动力方法”。相应的， 3DFLAC在进行动力分析时，通过求解动力方程理所当然地可以得到合适的动力问题解答。对于本构模型的选择，主要是描述单元的应力-应变关系，如果是弹塑性的，则考虑的是单元的屈服准则、流动法则等。 3D等效线性方法考虑土体的滞后性常常是通过将骨干曲线进行变换，比如Masing二倍法，而在FLAC的动力分析中，滞后性是通过阻尼来考虑，通过设置合适的阻尼形式和阻尼参数，同样可以描述土体在动力作用下的滞回曲线和滞回圈。 3D因此，FLAC动力分析中采用的本构模型可以选取任意模型，其参数也是对应静力本构模型的参数，关键是要设置合适的阻尼形式、阻尼参数、边界条件等，这些内容将在本章的后续内容中进行讲解。 11.2 动力时间步 动力计算中临界计算时间步的计算如下: ,,V,, (11-1) min,,t,,critfCA,,maxp,, 其中，为p波波速，与材料的体积模量K和剪切模量G有关，可以示为: Cp KG，4/3C, (11-2) p, f为四面体子单元(sub-zone)的体积，为与四面体子单元相关的最大表面积，表示遍历VAmin{}max 所有的单元，包括结构单元和接触面单元。 由于式(11-1)只是临界时间步的一个估计值，因此在使用中采用了一个安全系数，乘以0.5。因此，当采用无刚度比例的阻尼时，动力分析的时间步为: (11-3) ,,,tt2dcrit 如果采用了刚度比例的阻尼，那么为了保持数值稳定性，时间步必须减小。Belytschko(1983)提出了一个临界时间步的公式，其中考虑了刚度比例阻尼的影响。 ,t, ,,22 (11-4) ,,,,，,t1,,，,，,,,max 其中，为系统的最高特征频率，为该频率下的临界阻尼比。 ,,max 3D注意:FLAC在动力计算中～程序会根据数值计算的稳定性自动设臵动力计算时间步～一般不建议 读者对这个默认的时间步进行放大。甚至～在大应变计算过程中～如果出现很大的网格变形并导致 网格的几何错误时～还要对默认的时间步进行折减～降低动力时间步～以达到数值稳定的目的。 11.3 动态多步 3D由式(11-1)可知，FLAC动力计算中时间步需要遍历所有单元，取所有单元临界时间步中的最小值，因此时间步是由几何尺寸较小、模量较大的单元来确定的。因此，在计算中，尤其是在试算期间，要尽量避免较小的单元尺寸及刚度很大的材料，比如用实体单元来模拟较薄的混凝土墙，这种情况下必然会使动力时间步非常小，从而造成计算时间过长。可以通过采用结构单元或暂时不考虑混凝土墙的办法来进行试算，等到有关参数调试完成后再进行细化计算。 3D当计算模型中存在刚度差异较大、模型网格尺寸不均匀的情况时，FLAC可以采用“动态多步”(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 1 old_time = clock end 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函数和变量应具有相同的名称。因为设臵动力边 3D界条件命令中的hist关键词后面必须要跟随一个FISH函数名～FISH变量需要调用FLAC中的内臵标量dytime～该变量是动力计算的真实时间～通过调用可以给函数提供预订变化的数值。所以～一般FISH定义动荷载的方法如下,以定义函数xxx为例,: def xxx xxx = …dytime end app xvel = 1.0 hist xxx range … 3. 计算过程与输出结果 计算过程中命令窗口会提示动力计算的步数、动力时间和时间步，计算结束后可以将模型底部和墙 体顶部节点的水平速度时程输出，使用以下的命令: plot hist 1 2 v 4 输出结果见图11-2所示。读者可以采用设置动态多步和不设置动态多步两种情况分别进行计算，观 3D察图11-2中的速度时程曲线以及FLAC命令窗口中输出结果(图11-3和图11-4)。可以发现，时程曲 线、动力计算的时间步、迭代步数均一致，不同的是设置动态多步的情况下花费的计算时间较少。 注意:程序运行的时间将由于计算机配臵的不同而存在差别。 从图11-3中的输出信息可以看出，动态多步将模型中的375个单元分成了两类，并提供了不同的时 步乘子，其中只有小部分单元(墙体单元65个)的时步乘子为1，其他大部分单元(土体单元310个) 拥有较高的时步乘子，这样可以大大加快计算的进度。 1.2 FLAC3D 3.00 1.0Step 821616:46:10 Wed Apr 02 2008 0.8 History 0.6 1 X-Velocity Gp 55 Linestyle -1.000e+000 <-> 1.000e+000 0.4 2 X-Velocity Gp 688 Linestyle -1.192e+000 <-> 1.340e+000 0.2 Vs. 4 Dynamic Time 0.0 1.217e-003 <-> 9.993e-001 -0.2 -0.4 -0.6 -0.8 -1.0 2.0 4.0 6.0 8.0Itasca Consulting Group, Inc.x10^-1Minneapolis, MN USA 图11-2 动力计算结束时模型底部和顶部的水平速度时程曲线 图11-3 设置动态多步情况下的输出信息 图11-4 未设置动态多步情况下的输出信息 11.4 动力荷载和边界条件 3D利用FLAC进行动力计算时，有以下3个问题需要读者认真考虑: , 动力荷载和边界条件; , 力学阻尼; , 模型中波的传播。 本节及后续的两节将分别针对以上三个问题展开讨论。 11.4.1 动力荷载的类型与施加方法 3DFLAC可以在模型边界或内部节点施加动荷载来模拟材料受到外部或内部动力作用下的反应，程 序允许的动力荷载输入可以是:(1)加速度时程，(2)速度时程，(3)应力(压力)时程，(4)集中力 时程。 动力荷载的施加采用APPLY命令，另外，采用APPLY Interior命令可以将加速度、速度和力的时程 施加到模型内部的节点上。动力荷载的形式主要有2种: , FISH函数。FISH函数表达的动力荷载往往比较规则，也常用于试算阶段的动力输入，因为试 算时可以不用过多考虑荷载的频率、基线校正等问题。本章例11.1中已经给出了一个FISH函 数作为动力荷载的例子，这里不再赘述。 , TABLE命令定义的表。常用于离散的动力荷载输入，包括地震波、实测振动数据、不规则动力 输入等。下面简要介绍利用TABLE命令形成的表作为动力荷载的方法。 3D表是FLAC中的一种数据，表中的数据成对出现，相当于两列的表格。表建立的基本命令是: 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 3D在FLAC动力计算中，动力荷载往往数据点很多，用命令输入的方法显然不便，因此常用编辑文 本文件的方法进行表的读入与调用，编辑文本文件的表有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等格 式。 3D完成的文本文件需要进行读入操作才可以供FLAC调用。采用TABLE read命令进行读入引用: table 1(ID号) read 文件名 读入后，读者可以使用PLOT table和PRINT table命令来查看生成的表文件是否读入正确。在进行 3DFLAC动力边界条件设置时，使用APPLY命令调用已读入的表，下面的命令以施加水平向速度荷载为例。 app xvel 1.0 hist table 1 range … 命令中的1.0表示表格y向数据的乘子，可以方便地控制荷载幅值的大小。另外，动荷载的输入可以沿着x，y，z的任意方向施加或者模型边界的法向和切向施加。 注意:特定的边界条件不能在相同的边界上进行混合施加～否则程序会提示出错。 11.4.2 边界条件的设置 在动力问题中，模型周围边界条件的选取是一个主要内容，因为边界上会存在波的反射，对动力分析的结果产生影响。把分析模型的范围设置得越大，分析结果就越好，但较大的模型会导致巨大的计算 3D负担。FLAC中提供了静止(粘性)边界和自由场边界两种边界条件来减少模型边界上的波的反射。 1. 静态边界 3D3DFLAC中允许采用静态边界(也称粘性边界，吸收边界)条件来吸收边界上的入射波。FLAC中的静态边界是Lysmer和Kuhlemeyer(1969)提出的，具体做法是在模型的法向和切向分别设置自由的阻尼器从而实现吸收入射波的目的，阻尼器提供的法向和切向粘性力分别为式(11-5)和(11-6)。 (11-5) tCv,,,npn (11-6) tCv,,,sss 其中，，分别为模型边界上法向和切向的速度分量，为介质密度，，分别为p波和s,CvvCpnss波的波速。 这种静态边界对于入射角大于30度的入射波基本能够完全吸收，对于入射角较小的波，比如面波，虽然仍有一定的吸收能力，但吸收不完全。静态边界可以加在整体坐标系上，也可以加载在倾斜边界的法向和切向上。如果在倾斜边界的法向和切向施加静态边界，则需要同时使用nquiet，dquiet，squiet条件。 整体坐标系的静态边界条件设置使用命令为: APPLY xquiet (yquiet, zquiet) range … 使用倾斜边界上的静态边界条件命令为: APPLY nquiet dquiet squiet range … 注意: ,1,施加动力边界条件后～这些边界上原先的静力边界条件将被自动去掉,free,～在动力荷载施 加期间～程序始终自动计算边界上的作用力～用户不能将静态边界条件去掉,也不能在静态边界上 施加加速度、速度边界条件～因为静态边界上的作用力是根据边界上的速度分量计算得到的～如果 再施加速度荷载就会使静态边界失效。若需要在静态边界上输入动荷载～则只能输入应力时程。可 以将加速度、速度时程通过转换公式,11-7,和,11-8,形成应力时程施加到静态边界上。 (11-7) ,,=-2()Cvnpn (11-8) ,,,-2()Cvsss ,,ns式中的～分别为施加在静态边界上的法向应力和切向应力～公式中的系数2表示施加的能 量中只有一半是向上传播作为动力输入的～另一半向边界下部传播。公式中的负号是为了使应力施 加后节点的速度能与实际一致。 ,2,对于动力荷载来源于模型内部,如隧道中的列车振动问题,的情况～可以将动力荷载直接施加 在节点上～这种情况下使用静态边界可以有效减小人工边界上的反射～并且不需要再施加下面提到 的自由场边界。 ,3,动力计算过程中应避免静力荷载的变化。比如～在一个已经施加底部静态边界的模型上进行开 挖～会造成整个模型向上移动。因为施加静态边界时程序自动计算了施加在模型底部边界上的反力～ 这些边界反力不能与开挖后的模型相平衡～会引起模型的整体上移。 下面给出一个简单的静态边界的例子。 问题描述:如图11-5所示，一根竖直的弹性杆，高50m，宽1m，体积模量和剪切模量分别为20MPa3和10MPa，材料密度为1000kg/m。在杆的底部边界的两个水平方向设置静态边界条件，而杆的顶部为自由表面，在杆的底部施加水平方向的应力冲击荷载。根据材料参数，可得计算得到剪切波速为100 Cs55m/s，的乘积为10。应力冲击荷载的幅值设置为2×10，根据公式(11-7)可以得到等效速度荷载,Cs 幅值为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的系数来源于的值 ,Csapply 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 静态边界条件示例 1.8 FLAC3D 3.00 Step 1479 1.615:13:54 Thu Apr 03 2008 History 1.4顶部 1 X-Velocity Gp 1 Linestyle -3.836e-002 <-> 1.001e+000 2 X-Velocity Gp 101 1.2 Linestyle -3.348e-002 <-> 1.000e+000 3 X-Velocity Gp 201底部 中心 1.0 Linestyle -5.125e-002 <-> 1.995e+000 Vs. 0.8 4 Dynamic Time 1.217e-002 <-> 1.789e+000 0.6 0.4 0.2 0.0 0.5 1.0 1.5Itasca Consulting Group, Inc.Minneapolis, MN USA 图11-6 在静态边界上输入应力荷载得到的响应 2. 自由场边界 对诸如大坝之类的地面结构进行动力反应分析时，在模型各侧面的边界条件须考虑为没有地面结构 3D时的自由场运动。FLAC通过在模型四周生成二维和一维网格的方法来实现这种自由场边界条件(见图11-7所示)，主体网格的侧边界通过阻尼器与自由场网格进行耦合，自由场网格的不平衡力施加到主体网格的边界上。由于自由场边界提供了与无限场地相同的效果，因此向上的面波在边界上不会产生扭曲。 图11-7 自由场边界示意图 3D对于FLAC而言，自由场边界模型包括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个步骤: 去掉模型底部已有的静力约束条件; 施加静态边界条件; 施加动力荷载; 施加自由场边界条件。 施加自由场边界条件以后，程序自动在主体网格周围生成一圈自由场网格。当计算模型的单元数较 多时，这个生成过程需要花费较多的时间。可以对生成的自由场网格定义group，如例11.3中就将二维 自由场网格和一维自由场网格分别赋值了group，这样便于模型的观察和主体网格结果的后处理。 图11-8 自由场边界施加前后的网格 计算中对主体网格、柱体自由场网格和平面自由场网格对应位置的水平向速度进行了监测，计算结 果见图11-9所示。可以发现，主体网格与周围的自由场网格同步运动，达到了自由场网格的目的。 x10^-2x10^-1 9.0 1.0 8.0FLAC3D 3.00FLAC3D 3.00Step 992Step 99217:04:15 Thu Apr 03 200817:04:45 Thu Apr 03 2008 0.8 7.0HistoryHistory 3 X-Velocity Gp 10 4 X-Velocity Gp 102 Linestyle Linestyle 6.0 6.711e-006 <-> 9.060e-002 -6.480e-004 <-> 1.005e-001 5 X-Velocity Gp 368 6 X-Velocity Gp 380 0.6 Linestyle Linestyle 6.711e-006 <-> 9.049e-002 -3.849e-005 <-> 1.001e-001 5.0 7 X-Velocity Gp 145 8 X-Velocity Gp 176 Linestyle Linestyle 6.711e-006 <-> 9.049e-002 -3.849e-005 <-> 1.001e-001 9 X-Velocity Gp 232 10 X-Velocity Gp 274 4.0 Linestyle Linestyle 6.711e-006 <-> 9.063e-002 -6.919e-004 <-> 1.007e-001 0.4 Vs. Vs. 3.0 2 Dynamic Time 2 Dynamic Time 5.735e-005 <-> 1.497e-002 5.735e-005 <-> 1.497e-002 2.0 0.2 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.2 0.4 0.6 0.8 1.0 1.2 1.4Itasca Consulting Group, Inc.Itasca Consulting Group, Inc.x10^-2x10^-2Minneapolis, MN USAMinneapolis, MN USA 模型底部 模型顶部 图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 地震荷载的输入 3D地震反应分析是FLAC动力计算的主要应用领域，所以有必要对地震荷载的输入做单独介绍。地 3D震荷载常用表的文本文件格式读入到FLAC，再施加到模型底部的边界上。对于地震荷载的动力边界 条件选择可以参考以下: 1. 刚性地基 如果模型底部为岩石等模量较大的材料，可以在底部直接施加加速度或速度荷载，并采用自由场边界条件，模型底部无需施加静态边界条件。 2. 柔性地基 如果模型底部的单元为土体，尤其是软土，则不能直接施加加速度和速度，而需要将加速度、速度转换成应力时程，再施加到模型底部。模型周围采用自由场边界条件，模型底部可采用静态边界条件(nquiet，dquiet，squiet)。 11.5 力学阻尼 3D阻尼的产生主要来源于材料的内部摩擦以及可能存在的接触表面的滑动。FLAC采用求解动力方程的方法解决两类力学问题:准静力问题和动力问题。在这两类问题中都要使用阻尼，但准静力问题需要更多的阻尼使得动力方程能够更快的收敛平衡。对于动力问题中的阻尼，需要在数值模拟中重现自然 3D系统在动荷载作用下的阻尼大小。目前，FLAC动力计算提供了三种阻尼形式供用户选择，分别是瑞利阻尼、局部阻尼和滞后阻尼。 11.5.1 瑞利阻尼(Rayleigh damping) 瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关: CMK,，,, (11-9) ,其中，为与质量成比例的阻尼常数，为与刚度成比例的阻尼常数。 , 瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。图11-10为归一化的临界阻尼比与角频率之间的关系曲线，其中包含了仅设置刚度分量、仅设置质量分量以及二者叠加的三种结果。可以看出，采用叠加的方法得到的阻尼比在较大的频率范围内保持定值，因此瑞利阻尼可以近似地反映岩土体具有的频率无关性。 叠加 仅有刚度分量 ,=0 仅有质量分量 ,=0 图11-10 归一化的临界阻尼比与角频率之间的关系 1. 瑞利阻尼的两个参数 设置瑞利阻尼的命令为: SET dyn damp rayleigh 参数1 参数2 根据图11-10叠加结果的阻尼比曲线最小值可以确定瑞利阻尼的两个参数，分别是最小临界阻尼比 (参数1)和最小中心频率(参数2)，可以按照式(11-10)进行计算: ,,minmin 1/2,,,,,()min (11-10) 1/2,(/),,,min 注意:参数2中的单位是Hz～而不是角频率的单位。 2. 瑞利阻尼参数选择的原则 很多读者最十分关心瑞利阻尼参数的选择。对于岩土材料而言，临界阻尼比的范围一般是2,5%，而结构系统的临界阻尼比为2,10%。在使用弹塑性模型进行动力计算时，相当多的能量消散于材料发生的塑性流动阶段，因此在进行大应变的动力分析时，只需要设置一个很小的阻尼比(比如0.5%)就能满足要求，而且达到塑形以后，随着应力-应变滞回圈的扩大，能量的消散也逐渐明显。 (1)中心频率的确定 瑞利阻尼是与频率相关的，但是在一定的频率范围内，瑞利阻尼基本与频率无关。这个频率范围的最大频率常是最小频率的3倍。对于任何动力问题，可以对速度时程进行谱分析得到速度谱与频率之间的关系，见图11-11所示。可以逐渐调整，使得频率范围在,3×之间包含了动力能量的主要fffminminmin 部分，这时的就是瑞利阻尼的中心频率。在实际分析中，比如多种材料的土石坝动力分析，需要首fmin 先将分析模型假设为弹性材料进行动力计算，得到各种材料关键位置的功率谱，根据功率谱的分布来确定该区域的瑞利阻尼中心频率的大小。通过这种方法确定的中心频率，既不是输入频率，也不是系fmin 统的自振频率，而是二者的叠加。 主频区域 速度谱 ,3× ,频率 minmin 图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。读者在确定振荡周期时，可以 2 变量hist_rep的作用是将设置历史记录的间隔，默认为10，即每10步记录一次数据。由于本例中计算时间步较少，为了得到更平滑的记录曲线，将hist_rep取为1。 x10^-3将计算结果输出到文本文件，再通过Excel操作得到更精确的周期值。 Job Title: vertical displacement versus time (undamped)-0.1FLAC3D 3.00View Title: Step 15011:05:17 Fri Apr 04 2008-0.2 History 1 Z-Displacement Gp 56-0.3 Linestyle -1.091e-003 <-> -3.527e-006-0.4 Vs. 2 Dynamic Time 5.939e-004 <-> 8.909e-002-0.5 -0.6 -0.7 -0.8 -0.9 -1.0 2.0 4.0 6.0 8.0Itasca Consulting Group, Inc.x10^-2Minneapolis, MN USA 图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 3pause ;(2)只有质量分量 3 命令pause的作用是在命令文件执行中设置暂停，以便于用户观察暂停状态下的计算结果，在命令行中输入continue 可以继续执行，若输入命令中存在错误，就停止执行。pause后面可以跟keyword关键词，表示暂停状态下按任意键即可 继续。 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所示。可以看出，瑞利阻尼中刚度分量的设 3D置大大降低了FLAC中的动力时间步。 表11-1 瑞利阻尼三种参数格式下的计算时间步比较 (1)质量分量和刚度计算工况 (2)只有质量分量 (3)只有刚度分量 分量共同作用 时间步 7.8E-5 5.9E-4 3.9E-5 x10^-4 Job Title: vertical displacement versus time (stiffness damping only)-0.5FLAC3D 3.00View Title: Step 202515:29:37 Fri Apr 04 2008-1.0History 1 Z-Displacement Gp 56-1.5 Linestyle -5.371e-004 <-> -1.561e-008 Vs.-2.0 2 Dynamic Time 3.951e-005 <-> 8.001e-002-2.5 -3.0 -3.5 -4.0 -4.5 -5.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0Itasca Consulting Group, Inc.x10^-2Minneapolis, MN USA 图11-13 竖向位移时程曲线的临界阻尼状态 11.5.2 局部阻尼(Local damping) 3D局部阻尼是FLAC静力计算中采用的阻尼形式，但是它的一些特性可以用来进行动力计算。它在振动循环中通过在节点或结构单元节点上增加或减小质量的方法达到收敛，由于增加的单元质量和减小的相等，因此总体来说，系统保持质量守恒。当节点速度的符号改变时增加节点质量，当速度达到最大值(或最小值)时减小节点质量。因此损失的能量是最大瞬时应变能的一定比例，这个比值，,W/W,WW 是与频率无关和率无关的。因为是临界阻尼比D的函数: ,W/W (11-11) ,,,DL 其中，为局部阻尼系数，临界阻尼比的取值可以参考瑞利阻尼中的。 ,,Lmin 局部阻尼的设置命令为: SET dyn damp local 局部阻尼系数 局部阻尼系数不用求解系统的自振频率，而且相对于瑞利阻尼而言不会减小时间步，从这个意义上来说具有较大的优势。但局部阻尼只适合于简单问题的求解，实践证明设置局部阻尼不能有效地衰减复杂波形的高频部分，计算结果会产生一些高频“噪声”。因此，在使用时需慎重，最好将局部阻尼与瑞利阻尼的结果进行对比分析。 继续沿用例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输出动力计算信息了解瑞利阻尼中两种分量的比例。 x10^-4x10^-4 -1.0-1.0Job Title: vertical displacement versus time (5% Rayleigh damping)Job Title: vertical displacement versus time (5% Local damping)FLAC3D 3.00FLAC3D 3.00View Title: View Title: Step 1420Step 1179-2.016:02:52 Fri Apr 04 200816:02:35 Fri Apr 04 2008-2.0HistoryHistory 1 Z-Displacement Gp 56 1 Z-Displacement Gp 56-3.0-3.0 Linestyle Linestyle -9.881e-004 <-> -3.623e-005 -9.888e-004 <-> -4.737e-005 Vs. Vs.-4.0-4.0 2 Dynamic Time 2 Dynamic Time 2.466e-003 <-> 7.002e-001 2.970e-003 <-> 6.978e-001-5.0-5.0 -6.0-6.0 -7.0-7.0 -8.0-8.0 -9.0-9.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 1.0 2.0 3.0 4.0 5.0 6.0Itasca Consulting Group, Inc.Itasca Consulting Group, Inc.x10^-1x10^-1Minneapolis, MN USAMinneapolis, MN USA 瑞利阻尼 局部阻尼 图11-14 瑞利阻尼与局部阻尼的比较 11.5.3 滞后阻尼(Hysteretic damping) 3DFLAC将土动力学中岩土体的滞后特性用阻尼的形式加入到程序中。使用模量衰减系数来描述Ms土体的非线性特性。假设土体为理想粘弹性体，可以从模量衰减曲线上得到归一化的剪应力: , (11-12),,,Ms dMd,s ,,， (11-13) MM,tsdd,, 式中，为剪应变，为归一化的割线模量，为归一化的切线模量。则增量剪切模量可以表,GMMst 示为: (11-14) GGM,0t 其中，为小应变下的剪切模量。 G0 图11-15 模量衰减曲线(Seed & Idriss) 滞后阻尼是与材料无关的阻尼格式，在动力计算中，滞后阻尼可以满足Masing二倍法，从而构造土体在动力作用下的滞回圈。另外，滞后阻尼的优点是: , 可以直接采用动力试验中的模量衰减曲线; , 相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步; , 可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。 滞后阻尼与瑞利阻尼及局部阻尼的设置不同，采用的是初始条件INITIAL命令。 INITIAL damp hyst 3DFLAC中滞后阻尼提供了多种形式的割线模量衰减曲线模型，包括默认模型(default)、S型模型 (包括三参数模型Sig3和四参数模型Sig4)以及哈丁模型(Hardin)。下面对这几种模型做简要介绍。 1. 默认模型 默认模型中曲线可以用式(11-15)的三次方程来拟合。 Ms 2 (11-15) Mss=(32),s 其中 LL-2 (11-16) s,LL-21 (11-17) L,log(),10 和为默认模型的两个参数，表示曲线的循环应变范围。比如=-3，=1表示曲线中循LLMLLM12s12s-31环应变的最小值是0.001%(10)，最大值是10%(10)，因此可见默认模型的参数确定较简单，设置方 法是: INITIAL damp hyst default LL12 2. S型模型 S型模型包括三参数模型和四参数模型，采用如下的公式来拟合曲线。 MsSig3模型包含a，b，x三个参数: 0 a M, (11-18) s1exp(-(-)/)，Lxb0 设置命令为: INITIAL damp hyst sig3 a b x0 Sig4模型包含a，b，x和y四个参数: 00 ay+ M, (11-19) s01exp(-(-)/)，Lxb0 设置命令为: INITIAL damp hyst sig4 a b x0 y0 3. 哈丁模型 滞后阻尼中有一种Hardin / Drnevich模型，采用式(11-20)的双曲线公式来拟合曲线: Ms 1M, (11-20) s,,，1/ref 其中，为参考应变，一般取G/G=0.5时的对应的应变值。 ,maxref 哈丁模型只有一个参数，设置命令为: INITIAL damp hyst hardin ,ref 表11-2提供了图11-15中模量衰减曲线的拟合结果，表11-3提供了黏土的模量衰减曲线拟合结果， 读者在使用时可参考。 表11-2 Seed & Idriss模量衰减曲线的拟合结果 数据来源 默认模型 三参数模型 四参数模型 哈丁模型 a = 0.9762 a = 1.014 砂土 L = -3.325 b = -0.4393 1 = 0.06 ,b = -0.4792 ref图11-15 L = 0.823 = -1.285 x20x = -1.249 0y = 0.03154 0 表11-3 Seed & Idriss模量衰减曲线的拟合结果(2) 数据来源 默认模型 三参数模型 四参数模型 哈丁模型 a = 0.922 a = 1.017 L = -3.156 b = -0.481 1 = 0.234 ,b = -0.587 粘土 refL = 1.904 = -0.705 x20x = -0.633 0y = 0.0823 0 3D注意:滞后阻尼是FLAC新开发的一个技术～在使用过程中读者要反复调试～同时也要注意以下几 点: , 低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。 可以在中心频率上增加一个较小分量的Rayleigh阻尼(比如0.2%刚度比例)，这样也不会降低时步。 , 若初始剪应力不为0，剪应力-剪应变曲线可能不匹配，因此在生成初始应力时就要调用 Hyst阻尼，这一点至关重要。 , 滞后阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入 波的基频接近共振频率的时候，可能会导致动力响应的幅值偏大。 , 在设置滞后阻尼之前要做一次弹性无阻尼求解，获得各关键部位发生循环应变的最大水平， 若循环应变过大导致剪切模量过多的降低，那么使用滞后阻尼可能会存在问题。 , 即使应变较小，使用塑性模型也会增大应变，因此若模型存在广泛屈服的现象时，不能使 用滞后阻尼。 11.5.4 关于阻尼设置的一些讨论 3D力学阻尼的设置是FLAC动力计算中讨论最多的话题，在这些阻尼形式中，瑞利阻尼由于其理论与常规动力分析方法类似，而且实践证明，瑞利阻尼计算得到的加速度响应规律比较符合实际，因此最易为大家所接受，唯一也是最大的不足就是瑞利阻尼的计算时间步太小，导致动力计算时间过长，因此很多用户不得不使用局部阻尼来代替。 3D滞后阻尼是新版本FLAC的一个亮点，但在实际使用过程中读者可以发现使用时存在一定的困难，主要原因是滞后阻尼有过多的使用限制，且目前相关的参考资料极少。作者曾经尝试过一些滞后阻尼的算例，当模型较复杂时，很难得到满意的分析结果，尤其是在处理地震液化分析时更是如此。因此读者在使用滞后阻尼时需要慎重，从简单做起，逐步了解其功能后再应用于实践。 注意:不同的阻尼形式之间可以混合使用,对于不同的材料也可以按照初始条件的方式设臵不同的 阻尼形式和参数,动力分析中如果存在结构单元～需要采用SEL set damp命令～指定结构单元的阻 尼～否则计算会提示出错。 11.6 网格尺寸的要求 输入波形的频率成分和土体的波速特性会影响波传播的数值精度。Kuhlemeyer和Lysmer(1973)的研究表明，要想精确描述模型中波的传播，那么网格的尺寸必须要小于输入波形最高频率对应的波长,l 的1/8到1/10，也就是: 11,, (11-21) ,,,l~,,810,, 式中，是最高频率对应的波长。 , 可见在动力计算中，土体的模量越小，即土体越软，最大网格尺寸越小，划分的网格数量越多。注意到，任何离散化的介质都存在能量传播的上限频率，只有当输入荷载的频率小于这个上限频率时，计 3D算结果才有意义。因此上述公式不仅适用于FLAC，同样适用于其他基于时域的动力分析程序。 由于输入荷载的最大频率直接影响到单元的最大尺寸，所以对于脉冲荷载、爆炸荷载等这些频率范围很广的荷载形式，需要进行滤波处理，这部分内容在本章的11.7.1节中将会做介绍。 11.7 输入荷载的校正 这里的输入荷载一般指的是地震荷载，因为在地震反应分析中，常使用离散的荷载列表，因此在施加之前有必要进行滤波和基线校正。 11.7.1 滤波 滤波的目的是过滤掉原有波形中的高频分量，因为由式(11-21)可知，地震波的最大频率对网格尺寸的影响较大，最大频率越高，满足精度条件下的网格尺寸越小。采用滤波的方式，可以减小地震波的最大频率，从而增大计算所需的最小网格尺寸，减小单元数量，达到节约计算时间的目的。滤波可以通 3D过OriginPro，SeismoSignal等软件进行，也可以使用FLAC提供的FFT.fis函数进行。 11.7.2 基线校正 3D在FLAC地震动力分析中，输入波通常为加速度时程，若将输入的加速度进行积分得到的最终速度和最终位移不为0，则在动力计算结束时模型底部会出现继续的速度和残余的位移，此时需要对加速度时程进行基线校正。即通过在原始加速度时程上增加一个低频率的波形(多项式或周期函数)，使最终的速度和位移均为0。基线校正可以通过SeismoSignal软件进行，该软件提供了多种基线校正的方法，建议读者使用。 11.8 动孔压模型与土体的液化 3DFLAC可以进行动力与渗流的耦合分析，能够模拟砂土在动力作用下的孔压积累直至土体的液化，3DFLAC采用了Finn模型来描述这种孔压积累的效应。Finn模型的实质是在Mohr-Coulomb模型的基础上增加了动孔压的上升模式，并假定动孔压的上升与塑形体积应变增量相关。 ,设在有效应力为时砂土的一维回弹模量为，则对于不排水条件下孔隙水压力的增量与塑性,uE,r0 体积应变增量的关系为: ,,vd (11-22) ,,,uE,rvd 3DFLAC提供了2种不同的塑形体积应变增量公式，包括Finn模式和Byrne模式。 1. Finn模式 [143]Martin et al.(1975)的试验表明，塑性体积应变与循环剪应变幅值之间的关系与固结压力无关。为了实用目的，塑性体积应变增量仅是总的累积体积应变和剪应变的函数: ,,,,vdvd 2C,3vd (11-23) ,,,，()CC,,,vdvd12，C,,vd4 其中，，，和为模型常数。对于相对密度为45%的结晶二氧化硅砂: = 0.80， = 0.79，CCCCCC123412 = 0.45， = 0.73。 CC34 2. Byrne模式 Byrne(1991)提出了一种更简便的计算塑性体积应变增量的方法: ,,,,,vdvdCCexp ,, (11-24) 12,,,,,, 其中，和为两个参数，大多数情况下两者存在如下的关系: CC12 0.4C, (11-25) 2C1 参数与砂土的相对密度存在如下关系: C1 -2.5 (11-26) CD,7600)(r1 另外，相对密度与标准贯入击数存在一定的经验关系: 1/2 (11-27) DN,15(1)r60 因此，参数也可以通过标准贯入击数来得到: C1 -1.25 (11-28) CN,8.7()1160 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分别为孔压上升模式公式中的C，C，C和C。 1234, ff_latency表示两次应变反转之间的最小步数。 11.9 完全非线性动力耦合分析步骤 3DFLAC采用完全非线性的动力分析方法，可以考虑动力与渗流的耦合分析，模拟土体的液化。动 力分析是在静力分析的基础上进行的，因此在动力计算之前要进行静力的力学计算和渗流计算，得到正 确的应力场和渗流场。在动力计算前要考虑网格尺寸、边界条件、材料参数、阻尼类型、地震波调整等 问题，一般动力分析过程见图11-16所示。 开始 静力分析 生成网格 网格尺寸满足 否 要求, 是 静力边界条件和初始条件 力学平衡 流体平衡 震前初始状态 选择材料模型及参数 地震波的调整 动力分析 模型中刚度差 基线校正滤波 异较大, 是 设置动态多步计算 输入动力荷载 设置动力边界条件 选择阻尼类型 进行动力求解 结束 图11-16 完全非线性动力耦合问题求解流程图 11.10 应用实例——振动台液化试验模拟 下面用一个振动台液化试验的例子来说明动力分析的主要过程。在本节的分析中，比较了不同阻尼形式、流固耦合方式对计算结果的影响。读者在这个例子的基础上通过进一步细化改进，可以模拟其它条件下的振动台试验。 11.10.1 计算模型及参数 水平场地的液化模拟采用的计算模型见图11-17所示。为了减少计算时间，本例中没有采用真实的振动台尺寸，而采用了较大的模型(长30m，高11m)。其中，可液化砂土厚度为10m，其上部有一层厚度为1m的非液化层，垂直平面的方向取1m。计算中考虑了动力计算中的最大网格尺寸的要求，如表11-4所示，模型的网格划分见图11-18所示，共划分330个单元。为了便于不同计算工况的比较，在计算中对可液化层中的一个单元进行跟踪，分别记录该单元的孔压和平均有效应力随动荷载时间的变化曲线。 顶部非液化层 1 可液化层 10 监测单元 30单位:m 图11-17 水平场地液化计算模型 图11-18 水平场地液化计算的网格划分 表11-4 水平场地液化计算中最大单元尺寸的计算 , ,lG f cG,/,s (MPa) (Hz) (m) (m) (m/s) 30 70.71 5 14.14 1.4144 计算时，首先利用弹性模型使土体在重力作用下达到平衡，得到动力计算前的初始应力场，然后在 模型底部及两侧的水平方向 (x方向) 施加正弦的速度边界，见图11-19所示。速度时程考虑了振动台试 验激振的过程，从0 ~ 2 s速度逐渐增大达到最大值并稳定一定时间 (2 ~ 20s)，随后速度逐渐减小到0 (20 ~ 30s)。计算中采用的速度时程曲线见图11-20所示。 计算工况为: 3D, 本构模型:分别采用FLAC自带的Finn模型; , 阻尼形式:采用局部阻尼、瑞利阻尼和滞后阻尼，三种阻尼形式的参数见表111-5所示; , 耦合形式:分别进行考虑动力‐渗流的耦合及不考虑耦合的计算。 表111-5 阻尼形式及参数选择 阻尼形式 阻尼参数 局部阻尼 临界阻尼比D = 10% 瑞利阻尼 阻尼比D = 10%，f = 3.67 Hz min 滞后阻尼 参考应变 = 6% r 图11-19 水平场地液化计算采用的边界条件 0.10 -10.05 / ms0.00 速度-0.05 -0.10 0510152025303540时间 / s 图11-20 水平场地液化计算采用的速度荷载时程 11.10.2 计算过程 计算过程主要包括静力分析和动力分析，其中静力分析是动力分析的基础。 1. 静力分析 在进行动力分析之前要进行静力计算，获得振动施加前的初始应力状态。本例中初始应力状态的获 得分为两步: 设置弹性材料参数、流体模型参数、初始孔压场条件和边界条件，关闭流体模式和动力模 式，并设置流体模量为0，在重力作用下达到静力平衡。 定义了一个设置水平应力的FISH函数，_k0是可以自行设定的侧压力系数，函数运行结束 后再次达到平衡，完成初始应力的计算。 静力分析阶段的命令如例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 apply xvel = 1.0 hist sine_wave ran x 29.9 30.1 set dyn damp local .314 call ppr.dat4 set hist_rep 100 set large ;set dyn dt 3e-4 set mech rat 1e-20 def solve_ages loop n(1,39) save_file = '11-8_'+string(n)+'s.sav' command sol age n ;save save_file endcommand endloop end solve_ages save save_file hist write 20 21 22 23 24 30 31 32 33 34 vs 2 file 10-8_Outfile_pp.txt hist write 40 41 42 43 44 vs 2 file 10-8_Outfile_xdis.txt 11.10.3 计算结果与分析 将Finn模型应用不同的阻尼形式得到的结果进行了汇总，见图11-21所示。可以发现，利用不同阻 尼形式得到的计算存在一定的差异，但超孔压时程和平均应力时程曲线的变化趋势基本一致。 本书中仅提供了不耦合模式下的计算结果，读者可以自行开展耦合模式下的动力计算，可以发现， 在耦合模式下，在动荷载作用到20s后，随着荷载幅值的减小，超孔隙水压力会逐渐消散，而在不耦合 模式下，积累的超孔压将不会消散。因此在实际分析中，如果需要考虑动孔压在震动过程中的积累和消 散效应，则在流固耦合的同时也需要将渗流计算打开。 4ppr..dat文件的作用是监测计算过程中一些变量的值，由于监测的内容较多，为了使命令的可读性更强，所以将历史 变量的监测单独形成了一个文件。 200 150 局部阻尼 p / kPa 瑞利阻尼100 滞后阻尼 0 -50' / kPa 局部阻尼 m 瑞利阻尼,-100 滞后阻尼 0510152025303540 时间 / s 图11-21 阻尼形式对计算结果的影响 (不耦合) 11.11 本章小结 3D本章以FLAC为例，对动力分析中边界条件的设置、阻尼形式的选择、网格尺寸的要求、动荷载 3D的调整、土体的液化等内容进行了介绍。FLAC / FLAC软件在动力分析方面具有强大的功能，同时也具有一定的难度。读者在进行动力分析前要对静力计算非常熟悉，因为初始应力计算的正确才能保证后续动力计算的顺利开展。 3DFLAC的动力计算往往要耗费很多的时间，主要原因是由于有限差分法本身需要较小的时间步，因此建议读者在进行动力分析前首先进行简单模型的分析，熟悉边界条件、荷载施加、阻尼选择等方面的基本内容，虽然在这些“小例子”中会花费读者一定的时间，但是这样会大大提高分析的效率。