FLAC3D基本原理和应用特点第三讲FLAC3D动力分析、自定义本构以及结构单元

1、FLAC / FLAC3D基本原理和应用特点1第三讲FLAC3D动力分析、自定义本构以及结构单元FLAC3D非线性动力分析3非常复杂！Said by Prof. Peter Cundall为什么要用FLAC做动力分析？FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。可以计算塑性引起的永久变形以及孔隙水压力的消散。土动力学中常用的等效线性方法无法直接处理上述问题。 4动力模拟的3个重要问题动力荷载与边界条件材料响应与阻尼土体液化5动力荷载动力输入的类型加速度时程速度时程应力(压力)时程力时程APPLY INTERIOR (内部)T

2、ABLEFISH6Quiet边界静态(quiet,粘性)边界Lysmer and Kuhlemeyer(1969) 模型边界法向和切向设置独立的阻尼器性能对于法向p波和s波能很好的吸收对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射人工边界仍应当足够远7Quiet边界应用内部振动(如隧道中的列车振动问题)动力荷载直接施加在节点上使用Quiet边界减小人工边界上的反射不需要FF边界外部荷载的底部边界软土地基上的地震荷载不适合用加速度或速度边界条件使用应力条件t = -2Csrvs地震底部输入的侧向边界扭曲了入射波8quietquietquietFree-field边界Cundall e

3、t al. (1980)自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上设置条件底部水平，重力方向为z向侧面垂直，法向分别为x, y向其他边界条件在APPLY ff之前9相当于一个阻尼器Free-field边界APPLY ff将边界上单元的属性、条件和变量全部转移ff单元上；设置以后主体网格上的改动将不会被FF边界所响应可存在任意的本构模型以及流体耦合(仅竖向)FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小存在attach的边界将不能设置FF边界边界上的Interface将不能连续动力边界设置需在FF边界设置之前10Free-field边界与

4、动力荷载模型底部边界fix施加速度或加速度荷载刚性边界Free施加应力时程荷载柔性边界对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载11Note that there is a factor of 2 because the input energy divides into a downward- & upward-propagating wave.2. 材料响应与阻尼12连续的非线性，表观模量随着应变的增大而降低对所有循环应变等级均存在滞回特性，因此导致随着循环应变的增加阻尼比增大。阻尼是率相关的。对于复杂波形的各个成分都产生阻尼。剪切应变会产生的体积应变，相应的，随着剪应

5、变循环次数的增加体积应变逐渐积累。材料响应13土体在循环荷载作用下呈现出模量衰减和能量消散的特点，那么如何用非线性数值方法对其进行模拟呢？Nonlinear characteristics of soils (Martin and Seed, 1979)试验得到的阻尼比、割线模量随循环剪应变的曲线14 0.00010.0010.010.11Shear Strain Amplitude (%)01020304050Damping Ratio (%)0.00010.0010.010.11Shear Strain Amplitude (%)0.00.10.20.30.40.50.60.70.80.9

6、1.0Normalized Shear Modulus, G/GmaxMid-Range Sand Curve(Seed & Idriss, 1970)Sand Fill Inland: Friction =32, hr=0.47, Go=440Sand Fill under Rock Dike: Friction=30, hr=0.43, Go=440等效线性方法15等效线性方法是岩土地震工程中模拟波的传播的最常用的方法。假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。Iteration toward strain-compa

7、tible shear modulus and damping ratio (after Kramer, 1996)等效线性方法的特点使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量

8、张量）是应力张量的函数。因此，塑性屈服的模拟不合理。大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。16完全非线性分析方法 FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运

9、动方程。 17完全非线性分析方法的特点可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。由于采用了弹塑性模型，因此程序可以自动计算永久变形。采用合理的塑性方程，使得塑性应变增量与应力相联系。可以方便地进行不同本构模型的比较。可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能

10、会动态地减小从而降低土体的抗剪强度。 18使用弹塑性模型附加考虑的因素：阻尼，对于屈服面以下应力的循环体积应变积累，是循环周数与幅值的函数模量衰减，基于平均应变水平的表格19最简单的弹塑性模型往往在描述累计塑性应变方面具有很好的效果，但是对于加速度放大系数的估算上效果不好。弹塑性模型20简单的理想弹塑性本构模型仅仅在发生屈服时才会出现滞回特性strainstress注意： 即使这样粗糙的模型也能够作出连续的阻尼比和模量衰减曲线。在屈服条件下会产生体积改变，但通常都是剪胀。FLAC3D中的阻尼比弹塑性本构模型使用瑞利（粘性）阻尼弹塑性本构模型使用滞后（HD）阻尼.复杂本构模型拥有连续的屈服应力应

11、变关系和对应的加卸载响应（Wang，UBCSand）21瑞利阻尼22瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关： 瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。 瑞利阻尼23frequencyratio of damping to critical混合仅有刚度分量仅有质量分量采用叠加的方法得到的阻尼比在较大的频率范围内保持定值 （3:1）Combi

12、ned curve reachesminimum at:中心频率的选择24假设弹性模型计算土石坝不同材料的功率谱曲线 (assuming elastic material)frequencyfrequencyfrequencyfrequencyfrequency阻尼比的选择25根据65%最大应变来选择阻尼比和模量衰减参参数瑞利阻尼的缺点必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率）刚度比例项会导致时间步减小26滞后阻尼27FLAC Version 5.0 and FLAC3D Version 3.0 提供了滞后阻尼功能，该项阻尼独立于材料模型之外，实质是包含

13、了一个与应变幅值相关的切线剪切模量乘子如果割线模量是由衰减曲线来确定的，那么切线模量可以表示为：From Seed & Idriss (1970)Go = small-strain shear modulus滞后阻尼28表观应变是偏应变，根据先前反转的点进行累计得到的。使用堆栈的数据结构可以描述应变反转现象。因此，在“迷你滞回圈”中的能量损失与主滞回圈的计算是一样的。elastic model with hysteresis damping特定的曲线可以使用滞后阻尼来描述，并在本构模型计算中根据应变的大小获得不同的切线模量。滞后阻尼的特点29可以直接采用动力试验中的模量衰减曲线；相对于瑞利阻尼

14、而言，滞后阻尼不影响动力计算的时间步；可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别“Good” fit to Seed & Idris data for G/Gmax (sigmoidal 3-parameter function) note inconsistent damping result.G/GmaxD - % of critical滞后阻尼低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个小量的Rayleigh阻尼(0.2%刚度比例)，这样也不会降低时步；若初始应力不为0

15、，剪应力-剪应变曲线可能不匹配。因此在生成初始应力时就要调用Hyst阻尼；Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大；Hyst阻尼之前要做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的；即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服模型，不用Hyst阻尼30动孔压的生成液化干沙剪应变循环加载试验初始加载阶段，沙土通常先压实再膨胀。卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在。

16、取决于初始孔隙率，这可能代表纯粹的压实假定孔隙中充满水对于常体积测试，有效应力降低，孔隙水压保持不变对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小有效应力为零时发生液化31动孔压的生成液化因此孔隙水压增加不是液化的基本原因由于颗粒间 (重组以后) 的低接触力导致有效应力的减小描述液化的模型高级模型：BSHP (边界面低塑性本构模型, Wang et al. 1990)简单模型：MC + 体积应变增量模型Finn模型：Byrne模型：32UBCTOT: Post-liquefaction33Assign post-liquefaction propertiestgtPL

17、-Finn34开始初始液化判断?零有效应力计算YES非零有效应力计算零有效应力判断?NO普通的Finn模型YESNO动力时间？结束NOYES0：零有效应力状态1：非零有效应力状态0有效应力非0有效应力地震波的调整基线校正对于地震分析的加速度时程，其积分得到的速度和位移应归0美国地质调查研究所Basic Strong-Motion Accelerogram Processing Software (BAP)对网格施加一个固定速度从而使残余的位移变为0动力荷载的频率与单元尺寸的双向调整高频的输入要求单元尺寸很小一定的单元尺寸对应输入的最大频率一般进行滤波处理滤掉低能量的高频OriginSeismo

18、Signal35地震波的调整36El-Centro波FFT修正后的时程修正后FFT5HzPL-Finn模型的应用阪神地震的码头分析37沉箱顶部水平残余位移最大达5m，平均为，残余沉降为1 2m，海侧倾斜角3 5 o。 sea Siltreplaced sandstoneland SandcaissonclayFLAC3D grid8050 zones10386 gps可液化的砂土分析参数Group nameConstitutive modelrd(kg/m3)E(MPa)mc(kPa)j(o)clayMC1350500.333020sea siltMC1250200.33030replaced

19、 sandMC1350150.33037land sandMC135013.70.33036stoneMC15501000.33040caissonElastic350020000.17Group nameFluid modelK (cm/s)PorositynDamping ratio Dliquefied parametersPL-Finn parametersclayfl_iso1.0E-60.450.05sea siltfl_iso1.0E-050.450.05replaced sandfl_iso1.0E-030.450.05Byrne ModelDr = 40 %C1 = 0.75

20、1C2 = 0.533C3 = 0pprc = 0.99k0 = 3105.4n0 = 0.3225k1 = 5503.1n1 = 0.1739land sandfl_iso1.0E-030.450.05Byrne ModelDr = 25 %C1 = 2.432C2 = 0.164C3 = 0pprc = 0.99k0 = 3105.4n0 = 0.3225k1 = 5503.1n1 = 0.1739stonefl_iso1.0E-010.450.05caissonfl_null0.0538基本力学参数 流体参数动力参数PL-Finn模型参数输入加速度39南北向地震波 最大值 = 0.6 g

21、竖直向地震波 最大值 = 0.2 g接触面与自由边界条件40沉箱与填石之间的接触面模型周围自由场边界条件计算结果位移41震后水平位移等值线沉箱顶点处的水平位移、沉降时程曲线3.44 m1.99 m计算结果（续）超孔压比42震后的超孔压比云图超孔压比时程曲线 A: 置换砂 B: 海砂（底部）C: 海砂（中部）自定义本构模型43自定义本构模型的基本方法必要性试验总结的本构模型特定条件下的本构模型交叉学科的本构模型二次开发环境自定义本构模型的功能自定义本构模型的基本方法44二次开发环境FLAC3D采用面向对象的语言标准C+编写本构模型都是以动态连接库文件(.DLL文件)的形式提供VC+2005或更高

22、版本的开发环境优点自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平 自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中45自定义本构模型的功能主要功能：对给出的应变增量得到新的应力辅助功能：模型名称、版本读写操作模型文件的编写基类(class Constitutive Model)的描述成员函数的描述模型的注册模型与FLAC3D之间的信息交换模型状态指示器的描述 46自定义本构模型的基本方法头文件(usermodel.h)中进行新的本构模型派生类的声明修改模型的ID(100)、名称和版本修改派生类的私有成员C+文件(usermo

23、del.cpp)中修改模型结构(UserModel:UserModel(bool bRegister): ConstitutiveModel)const char *UserModel:Properties()函数模型的参数名称字符串const char *UserModel:States()函数计算过程中的状态指示器47自定义本构模型的基本方法double UserModel:GetProperty()和void UserModel: SetProperty()函数const char * UserModel:Initialize()函数参数和状态指示器的初始化，并对派生类声明中定义的私有变

24、量进行赋值const char * UserModel:Run() 函数由应变增量计算得到应力增量，从而获得新的应力const char * UserModel:SaveRestore()函数对计算结果进行保存。程序的调试在VC+的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(Additional DLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试；在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。48相关文件49一个例子(Dunc

25、an-Chang)50其他成功的例子南京水科院双屈服面模型51双屈服面中心截面沉降云图水平位移云图小主应力云图结构单元52FLAC3D中的结构单元有限单元梁(beam)单元锚索(cable)单元桩(pile)单元锚杆: rockbolt壳(shell)单元格栅(geogrid)单元土工织物；土工格栅初衬(liner)单元53beamcablepileshellgeogridliner结构单元的应用土与结构的相互作用桩基；基坑；边坡锚固地下硐室的支撑结构；采矿；盾构土工织物；土工合成材料结构不宜复杂岩土工程软件，不宜单纯的结构分析复杂结构的模拟很困难plot显示双向接触结构(挡土墙)结构单元的厚

26、度54实际问题与FLAC模型55（a）实际问题中的桩（b）FLAC3D中的桩单元结构节点node结构构件SELs结构单元xyzsel pile id=1 beg 0 0 0 end 0 0 10 nseg 4术语56为什么使用SEL？Zone可以模拟任何物体Zones 对于弯矩的实现十分困难使用zone模拟结构时会产生数量众多的节点，同时建模复杂 如果不考虑结构小尺寸的细节响应（比如结构界面上的应力分布），则可以使用SEL来模拟结构57SEL LinksSEL nodes 通过link与其他node或zone进行联系两种类型的连接：node-zone node-nodeLink有相应的属性 (e.g. springs)Link可以与任何位置的grid进行联系，而不一定要与grid的坐标一致。58群桩 = 插秧默认的连接属性59建模SEL结构模型时，程序自动建立结构node与zone的连接 (node-zone links)Node-Node LinksSEL nodes 之间不会自动生成联系. 必须手动设置node之间的联系 (e.g., beam and cable) 这样才能使node产生相互作用.60连接 SELs612. 建立新 node-node link from SEL A to SEL B1. 删除node-zone link on SEL A