FLAC3D动力分析陈育民.ppt

1,FLAC / FLAC3D基本原理与应用,陈育民 博士,河海大学科学研究院,2,第三讲,FLAC3D动力分析、自定义本构、结构单元,3,FLAC3D非线性动力分析,非常复杂！,Said by Prof. Peter Cundall,4,为什么要用FLAC做动力分析？,FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。 可以计算塑性引起的永久变形以及孔隙水压力的消散。 土动力学中常用的等效线性方法无法直接处理上述问题。,5,动力模拟的3个重要问题,动力荷载与边界条件 材料响应与阻尼 土体液化,6,动力荷载,动力输入的类型 加速度时程 速度时程 应力(压力)时程 力时程 APPLY INTERIOR (内部) TABLE FISH,7,Quiet边界,静态(quiet,粘性)边界 Lysmer and Kuhlemeyer(1969) 模型边界法向和切向设置独立的阻尼器 性能 对于法向p波和s波能很好的吸收 对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射 人工边界仍应当足够远,8,Quiet边界应用,内部振动(如隧道中的列车振动问题) 动力荷载直接施加在节点上 使用Quiet边界减小人工边界上的反射 不需要FF边界 外部荷载的底部边界 软土地基上的地震荷载不适合用加速度或速度边界条件 使用应力条件t = -2Csrvs 地震底部输入的侧向边界 扭曲了入射波,quiet,quiet,quiet,9,Free-field边界,Cundall et al. (1980) 自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上 设置条件 底部水平，重力方向为z向 侧面垂直，法向分别为x, y向 其他边界条件在APPLY ff之前,相当于一个阻尼器,10,Free-field边界,APPLY ff将边界上单元的属性、条件和变量全部转移ff单元上； 设置以后主体网格上的改动将不会被FF边界所响应 可存在任意的本构模型以及流体耦合(仅竖向) FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小 存在attach的边界将不能设置FF边界 边界上的Interface将不能连续 动力边界设置需在FF边界设置之前,11,Free-field边界与动力荷载,模型底部边界 fix施加速度或加速度荷载刚性边界 Free施加应力时程荷载柔性边界 对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载,Note that there is a factor of 2 because the input energy divides into a downward- & upward-propagating wave.,12,连续的非线性，表观模量随着应变的增大而降低 对所有循环应变等级均存在滞回特性，因此导致随着循环应变的增加阻尼比增大。阻尼是率相关的。 对于复杂波形的各个成分都产生阻尼。 剪切应变会产生的体积应变，相应的，随着剪应变循环次数的增加体积应变逐渐积累。,2. 材料响应与阻尼,13,土体在循环荷载作用下呈现出模量衰减和能量消散的特点，那么如何用非线性数值方法对其进行模拟呢？,Nonlinear characteristics of soils (Martin and Seed, 1979),材料响应,14,试验得到的阻尼比、割线模量随循环剪应变的曲线,r,m,a,l,i,z,e,d,S,h,e,a,r,M,o,d,u,l,u,s,G,/,G,m,a,x,Mid-Range Sand Curve,(Seed & Idriss, 1970),Sand Fill Inland:,Friction =32, hr=0.47, Go=440,Sand Fill under Rock Dike:,Friction=30, hr=0.43, Go=440,15,等效线性方法是岩土地震工程中模拟波的传播的最常用的方法。 假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。,等效线性方法,Iteration toward strain-compatible shear modulus and damping ratio (after Kramer, 1996),16,等效线性方法的特点,使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。 不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。 塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。 大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。 本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。,17,完全非线性分析方法,FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。,18,完全非线性分析方法的特点,可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。 采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。 由于采用了弹塑性模型，因此程序可以自动计算永久变形。 采用合理的塑性方程，使得塑性应变增量与应力相联系。 可以方便地进行不同本构模型的比较。 可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。,19,使用弹塑性模型,附加考虑的因素： 阻尼，对于屈服面以下应力的循环 体积应变积累，是循环周数与幅值的函数 模量衰减，基于平均应变水平的表格,最简单的弹塑性模型往往在描述累计塑性应变方面具有很好的效果，但是对于加速度放大系数的估算上效果不好。,20,简单的理想弹塑性本构模型仅仅在发生屈服时才会出现滞回特性,注意： 即使这样粗糙的模型也能够作出连续的阻尼比和模量衰减曲线。,在屈服条件下会产生体积改变，但通常都是剪胀。,弹塑性模型,21,FLAC3D中的阻尼比,弹塑性本构模型使用瑞利（粘性）阻尼 弹塑性本构模型使用滞后（HD）阻尼. 复杂本构模型拥有连续的屈服应力应变关系和对应的加卸载响应（Wang，UBCSand）,22,瑞利阻尼,瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关：,瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。,23,frequency,ratio of damping to critical,混合,仅有刚度分量,仅有质量分量,采用叠加的方法得到的阻尼比在较大的频率范围内保持定值 （3:1）,Combined curve reaches minimum at:,瑞利阻尼,24,假设弹性模型计算土石坝不同材料的功率谱曲线 (assuming elastic material),frequency,frequency,frequency,frequency,frequency,中心频率的选择,25,根据65%最大应变来选择阻尼比和模量衰减参参数,阻尼比的选择,26,必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率） 刚度比例项会导致时间步减小,瑞利阻尼的缺点,27,FLAC Version 5.0 and FLAC3D Version 3.0 提供了滞后阻尼功能，该项阻尼独立于材料模型之外，实质是包含了一个与应变幅值相关的切线剪切模量乘子,如果割线模量是由衰减曲线来确定的，那么切线模量可以表示为：,From Seed & Idriss (1970),Go = small-strain shear modulus,滞后阻尼,28,表观应变是偏应变，根据先前反转的点进行累计得到的。使用堆栈的数据结构可以描述应变反转现象。,因此，在“迷你滞回圈”中的能量损失与主滞回圈的计算是一样的。,elastic model with hysteresis damping,滞后阻尼,特定的曲线可以使用滞后阻尼来描述，并在本构模型计算中根据应变的大小获得不同的切线模量。,29,滞后阻尼的特点,可以直接采用动力试验中的模量衰减曲线； 相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步； 可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。,一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别,“Good” fit to Seed & Idris data for G/Gmax (sigmoidal 3-parameter function) note inconsistent damping result.,G/Gmax,D - % of critical,30,滞后阻尼,低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个小量的Rayleigh阻尼(0.2%刚度比例)，这样也不会降低时步； 若初始应力不为0，剪应力-剪应变曲线可能不匹配。因此在生成初始应力时就要调用Hyst阻尼； Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大； Hyst阻尼之前要做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的； 即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服模型，不用Hyst阻尼,31,动孔压的生成液化,干沙剪应变循环加载试验 初始加载阶段，沙土通常先压实再膨胀。卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在。取决于初始孔隙率，这可能代表纯粹的压实 假定孔隙中充满水 对于常体积测试，有效应力降低，孔隙水压保持不变 对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小 有效应力为零时发生液化,32,动孔压的生成液化,因此孔隙水压增加不是液化的基本原因 由于颗粒间 (重组以后) 的低接触力导致有效应力的减小 描述液化的模型 高级模型：BSHP (边界面低塑性本构模型, Wang et al. 1990) 简单模型：MC + 体积应变增量模型 Finn模型： Byrne模型：,33,UBCTOT: Post-liquefaction,Assign post-liquefaction properties,t,g,t,34,PL-Finn,0：零有效应力状态 1：非零有效应力状态,0有效应力,非0有效应力,35,地震波的调整,基线校正 对于地震分析的加速度时程，其积分得到的速度和位移应归0 美国地质调查研究所 Basic Strong-Motion Accelerogram Processing Software (BAP) 对网格施加一个固定速度从而使残余的位移变为0 动力荷载的频率与单元尺寸的双向调整 高频的输入要求单元尺寸很小 一定的单元尺寸对应输入的最大频率 一般进行滤波处理 滤掉低能量的高频 FFT.FIS Origin SeismoSignal,36,地震波的调整,El-Centro波,FFT,修正后的时程,修正后FFT,5Hz,37,PL-Finn模型的应用,阪神地震的码头分析,沉箱顶部水平残余位移最大达5m，平均为3.5m，残余沉降为1 2m，海侧倾斜角3 5 o。,sea Silt replaced sand stone land Sand caisson clay,FLAC3D grid 8050 zones 10386 gps,可液化的砂土,38,分析参数,基本力学参数,流体参数,动力参数,PL-Finn模型参数,39,输入加速度,南北向地震波 最大值 = 0.6 g 竖直向地震波 最大值 = 0.2 g,40,接触面与自由边界条件,沉箱与填石之间 的接触面,模型周围自由场 边界条件,41,计算结果,位移,震后水平位移等值线,沉箱顶点处的水平位移、沉降 时程曲线,3.44 m,1.99 m,42,计算结果（续）,超孔压比,震后的超孔压比云图,超孔压比时程曲线,A: 置换砂 B: 海砂（底部） C: 海砂（中部）,43,自定义本构模型,44,自定义本构模型的基本方法,必要性 试验总结的本构模型 特定条件下的本构模型 交叉学科的本构模型 二次开发环境 自定义本构模型的功能 自定义本构模型的基本方法,45,二次开发环境,FLAC3D采用面向对象的语言标准C+编写 本构模型都是以动态连接库文件(.DLL文件)的形式提供 VC+2005或更高版本的开发环境 优点 自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平 自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中,46,自定义本构模型的功能,主要功能：对给出的应变增量得到新的应力 辅助功能： 模型名称、版本 读写操作 模型文件的编写 基类(class Constitutive Model)的描述 成员函数的描述 模型的注册 模型与FLAC3D之间的信息交换 模型状态指示器的描述,47,自定义本构模型的基本方法,头文件(usermodel.h)中进行新的本构模型派生类的声明 修改模型的ID(100)、名称和版本 修改派生类的私有成员 C+文件(usermodel.cpp)中修改模型结构 (UserModel:UserModel(bool bRegister): ConstitutiveModel) const char *UserModel:Properties()函数 模型的参数名称字符串 const char *UserModel:States()函数 计算过程中的状态指示器,48,自定义本构模型的基本方法,double UserModel:GetProperty()和void UserModel: SetProperty()函数 const char * UserModel:Initialize()函数 参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值 const char * UserModel:Run() 函数 由应变增量计算得到应力增量，从而获得新的应力 const char * UserModel:SaveRestore()函数 对计算结果进行保存。 程序的调试 在VC+的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(Additional DLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试； 在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。,49,相关文件,mohr.h,mohr.cpp,50,一个例子(Duncan-Chang),51,其他成功的例子,南京水科院双屈服面模型,双屈服面,中心截面沉降云图,水平位移云图,小主应力云图,52,结构单元,53,FLAC3D中的结构单元,有限单元 梁(beam)单元 锚索(cable)单元 桩(pile)单元 锚杆: rockbolt 壳(shell)单元 格栅(geogrid)单元 土工织物；土工格栅 初衬(liner)单元,beam,cable,pile,shell,geogrid,liner,54,结构单元的应用,土与结构的相互作用 桩基；基坑；边坡锚固 地下硐室的支撑结构；采矿；盾构 土工织物；土工合成材料 结构不宜复杂 岩土工程软件，不宜单纯的结构分析 复杂结构的模拟很困难 plot显示 双向接触结构(挡土墙)FLAC3D 3.1 结构单元的厚度,55,实际问题与FLAC模型,sel pile id=1 beg 0 0 0 end 0 0 10 nseg 4,56,术语,57,Zone可以模拟任何物体 Zones 对于弯矩的实现十分困难 使用zone模拟结构时会产生数量众多的节点，同时建模复杂 如果不考虑结构小尺寸的细节响应（比如结构界面上的应力分布），则可以使用SEL来模拟结构,为什么使用SEL？,58,SEL Links,SEL nodes 通过link与其他node或zone进行联系 两种类型的连接： node-zone node-node Link有相应的属性 (e.g. springs),Link可以与任何位置的grid进行联系，而不一定要与grid的坐标一致。,群桩 = 插秧,59,默认的连接属性,建模SEL结构模型时，程序自动建立结构node与zone的连接 (node-zone links),60,Node-Node Links,SEL nodes 之间不会自动生成联系. 必须手动设置node之间的联系 (e.g., beam and cable) 这样才能使node产生相互作用.,61,连接 SELs,2. 建立新 node-node link from SEL A to SEL B,1. 删除node-zone link on SEL A,SEL nodes connect to grid or other SEL nodes with links,62,错误！,连接 SELs,63,连接 SELs,正确！,64,Liner结构单元,三节点扁平有限单元 每个节点