FLAC3D动力分析陈育民课件.ppt

1,FLAC/FLAC3D基本原理与应用,陈育民博士,河海大学科学研究院,2,第三讲,FLAC3D动力分析、自定义本构、结构单元,3,FLAC3D非线性动力分析,非常复杂！,SaidbyProf.PeterCundall,4,为什么要用FLAC做动力分析？,FLAC可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。可以计算塑性引起的永久变形以及孔隙水压力的消散。土动力学中常用的等效线性方法无法直接处理上述问题。,5,动力模拟的3个重要问题,动力荷载与边界条件材料响应与阻尼土体液化,6,动力荷载,动力输入的类型加速度时程速度时程应力(压力)时程力时程APPLYINTERIOR(内部)TABLEFISH,7,Quiet边界,静态(quiet,粘性)边界LysmerandKuhlemeyer(1969)模型边界法向和切向设置独立的阻尼器性能对于法向p波和s波能很好的吸收对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射人工边界仍应当足够远,8,Quiet边界应用,内部振动(如隧道中的列车振动问题)动力荷载直接施加在节点上使用Quiet边界减小人工边界上的反射不需要FF边界外部荷载的底部边界软土地基上的地震荷载不适合用加速度或速度边界条件使用应力条件t=-2Csrvs地震底部输入的侧向边界扭曲了入射波,quiet,quiet,quiet,9,Free-field边界,Cundalletal.(1980)自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上设置条件底部水平，重力方向为z向侧面垂直，法向分别为x,y向其他边界条件在APPLYff之前,相当于一个阻尼器,10,Free-field边界,APPLYff将边界上单元的属性、条件和变量全部转移ff单元上；设置以后主体网格上的改动将不会被FF边界所响应可存在任意的本构模型以及流体耦合(仅竖向)FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小存在attach的边界将不能设置FF边界边界上的Interface将不能连续动力边界设置需在FF边界设置之前,11,Free-field边界与动力荷载,模型底部边界fix施加速度或加速度荷载刚性边界Free施加应力时程荷载柔性边界对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载,Notethatthereisafactorof2becausetheinputenergydividesintoadownward-&upward-propagatingwave.,12,连续的非线性，表观模量随着应变的增大而降低对所有循环应变等级均存在滞回特性，因此导致随着循环应变的增加阻尼比增大。阻尼是率相关的。对于复杂波形的各个成分都产生阻尼。剪切应变会产生的体积应变，相应的，随着剪应变循环次数的增加体积应变逐渐积累。,2.材料响应与阻尼,13,土体在循环荷载作用下呈现出模量衰减和能量消散的特点，那么如何用非线性数值方法对其进行模拟呢？,Nonlinearcharacteristicsofsoils(MartinandSeed,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-RangeSandCurve,(Seed&Idriss,1970),SandFillInland:,Friction=32,hr=0.47,Go=440,SandFillunderRockDike:,Friction=30,hr=0.43,Go=440,15,等效线性方法是岩土地震工程中模拟波的传播的最常用的方法。假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。,等效线性方法,Iterationtowardstrain-compatibleshearmodulusanddampingratio(afterKramer,1996),16,等效线性方法的特点,使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移。塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。,17,完全非线性分析方法,FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。,18,完全非线性分析方法的特点,可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。由于采用了弹塑性模型，因此程序可以自动计算永久变形。采用合理的塑性方程，使得塑性应变增量与应力相联系。可以方便地进行不同本构模型的比较。可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。,19,使用弹塑性模型,附加考虑的因素：阻尼，对于屈服面以下应力的循环体积应变积累，是循环周数与幅值的函数模量衰减，基于平均应变水平的表格,最简单的弹塑性模型往往在描述累计塑性应变方面具有很好的效果，但是对于加速度放大系数的估算上效果不好。,20,简单的理想弹塑性本构模型仅仅在发生屈服时才会出现滞回特性,注意：即使这样粗糙的模型也能够作出连续的阻尼比和模量衰减曲线。,在屈服条件下会产生体积改变，但通常都是剪胀。,弹塑性模型,21,FLAC3D中的阻尼比,弹塑性本构模型使用瑞利（粘性）阻尼弹塑性本构模型使用滞后（HD）阻尼.复杂本构模型拥有连续的屈服应力应变关系和对应的加卸载响应（Wang，UBCSand）,22,瑞利阻尼,瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关：,瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应。,23,frequency,ratioofdampingtocritical,混合,仅有刚度分量,仅有质量分量,采用叠加的方法得到的阻尼比在较大的频率范围内保持定值（3:1）,Combinedcurvereachesminimumat:,瑞利阻尼,24,假设弹性模型计算土石坝不同材料的功率谱曲线(assumingelasticmaterial),frequency,frequency,frequency,frequency,frequency,中心频率的选择,25,根据65%最大应变来选择阻尼比和模量衰减参参数,阻尼比的选择,26,必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率）刚度比例项会导致时间步减小,瑞利阻尼的缺点,27,FLACVersion5.0andFLAC3DVersion3.0提供了滞后阻尼功能，该项阻尼独立于材料模型之外，实质是包含了一个与应变幅值相关的切线剪切模量乘子,如果割线模量是由衰减曲线来确定的，那么切线模量可以表示为：,FromSeed&Idriss(1970),Go=small-strainshearmodulus,滞后阻尼,28,表观应变是偏应变，根据先前反转的点进行累计得到的。使用堆栈的数据结构可以描述应变反转现象。,因此，在“迷你滞回圈”中的能量损失与主滞回圈的计算是一样的。,elasticmodelwithhysteresisdamping,滞后阻尼,特定的曲线可以使用滞后阻尼来描述，并在本构模型计算中根据应变的大小获得不同的切线模量。,29,滞后阻尼的特点,可以直接采用动力试验中的模量衰减曲线；相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步；可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。,一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别,“Good”fittoSeed&IdrisdataforG/Gmax(sigmoidal3-parameterfunction)noteinconsistentdampingresult.,G/Gmax,D-%ofcritical,30,滞后阻尼,低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个小量的Rayleigh阻尼(0.2%刚度比例)，这样也不会降低时步；若初始应力不为0，剪应力-剪应变曲线可能不匹配。因此在生成初始应力时就要调用Hyst阻尼；Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大；Hyst阻尼之前要做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的；即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服模型，不用Hyst阻尼,31,动孔压的生成液化,干沙剪应变循环加载试验初始加载阶段，沙土通常先压实再膨胀。卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在。取决于初始孔隙率，这可能代表纯粹的压实假定孔隙中充满水对于常体积测试，有效应力降低，孔隙水压保持不变对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小有效应力为零时发生液化,32,动孔压的生成液化,因此孔隙水压增加不是液化的基本原因由于颗粒间(重组以后)的低接触力导致有效应力的减小描述液化的模型高级模型：BSHP(边界面低塑性本构模型,Wangetal.1990)简单模型：MC+体积应变增量模型Finn模型：Byrne模型：,33,UBCTOT:Post-liquefaction,Assignpost-liquefactionproperties,t,g,t,34,PL-Finn,0：零有效应力状态1：非零有效应力状态,0有效应力,非0有效应力,35,地震波的调整,基线校正对于地震分析的加速度时程，其积分得到的速度和位移应归0美国地质调查研究所BasicStrong-MotionAccelerogramProcessingSoftware(BAP)对网格施加一个固定速度从而使残余的位移变为0动力荷载的频率与单元尺寸的双向调整高频的输入要求单元尺寸很小一定的单元尺寸对应输入的最大频率一般进行滤波处理滤掉低能量的高频FFT.FISOriginSeismoSignal,36,地震波的调整,El-Centro波,FFT,修正后的时程,修正后FFT,5Hz,2019/12/13,37,可编辑,38,PL-Finn模型的应用,阪神地震的码头分析,沉箱顶部水平残余位移最大达5m，平均为3.5m，残余沉降为12m，海侧倾斜角35o。,seaSiltreplacedsandstonelandSandcaissonclay,FLAC3Dgrid8050zones10386gps,可液化的砂土,39,分析参数,基本力学参数,流体参数,动力参数,PL-Finn模型参数,40,输入加速度,南北向地震波最大值=0.6g竖直向地震波最大值=0.2g,41,接触面与自由边界条件,沉箱与填石之间的接触面,模型周围自由场边界条件,42,计算结果,位移,震后水平位移等值线,沉箱顶点处的水平位移、沉降时程曲线,3.44m,1.99m,43,计算结果（续）,超孔压比,震后的超孔压比云图,超孔压比时程曲线,A:置换砂B:海砂（底部）C:海砂（中部）,44,自定义本构模型,45,自定义本构模型的基本方法,必要性试验总结的本构模型特定条件下的本构模型交叉学科的本构模型二次开发环境自定义本构模型的功能自定义本构模型的基本方法,46,二次开发环境,FLAC3D采用面向对象的语言标准C+编写本构模型都是以动态连接库文件(.DLL文件)的形式提供VC+2005或更高版本的开发环境优点自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中,47,自定义本构模型的功能,主要功能：对给出的应变增量得到新的应力辅助功能：模型名称、版本读写操作模型文件的编写基类(classConstitutiveModel)的描述成员函数的描述模型的注册模型与FLAC3D之间的信息交换模型状态指示器的描述,48,自定义本构模型的基本方法,头文件(usermodel.h)中进行新的本构模型派生类的声明修改模型的ID(100)、名称和版本修改派生类的私有成员C+文件(usermodel.cpp)中修改模型结构(UserModel:UserModel(boolbRegister):ConstitutiveModel)constchar*UserModel:Properties()函数模型的参数名称字符串constchar*UserModel:States()函数计算过程中的状态指示器,49,自定义本构模型的基本方法,doubleUserModel:GetProperty()和voidUserModel:SetProperty()函数constchar*UserModel:Initialize()函数参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值constchar*UserModel:Run()函数由应变增量计算得到应力增量，从而获得新的应力constchar*UserModel:SaveRestore()函数对计算结果进行保存。程序的调试在VC+的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(AdditionalDLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试；在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。,50,相关文件,mohr.h,mohr.cpp,51,一个例子(Duncan-Chang),52,其他成功的例子,南京水科院双屈服面模型,双屈服面,中心截面沉降云图,水平位移云图,小主应力云图,53,结构单元,54,FLAC3D中的结构单元,有限单元梁(beam)单元锚索(cable)单元桩(pile)单元锚杆:rockbolt壳(shell)单元格栅(geogrid)单元土工织物；土工格栅初衬(liner)单元,beam,cable,pile,shell,geogrid,liner,55,结构单元的应用,土与结构的相互作用桩基；基坑；边坡锚固地下硐室的支撑结构；采矿；盾构土工织物；土工合成材料结构不宜复杂岩土工程软件，不宜单纯的结构分析复杂结构的模拟很困难plot显示双向接触结构(挡土墙)FLAC3D3.1结构单元的厚度,56,实际问题与FLAC模型,selpileid=1beg000end0010nseg4,57,术语,58,Zone可以模拟任何物体Zones对于弯矩的实现十分困难使用zone模拟结构时会产生数量众多的节点，同时建模复杂如果不考虑结构小尺寸的细节响应（比如结构界面上的应力分布），则可以使用SEL来模拟结构,为什么使用SEL？,59,SELLinks,SELnodes通过link与其他node或zone进行联系两种类型的连接：node-zonenode-nodeLink有相应的属性(e.g.springs),Link可以与任何位置的grid进行联系，而不一定要与grid的坐标一致。,群桩=插秧,60,默认的连接属性,建模SEL结构模型时，程序自动建立结构node与zone的连接(node-zonelinks),61,Node-NodeLinks,SELnodes之间不会自动生成联系.必须手动设置node之间的联系(e.g.,beamandcable)这样才能使node产生相互作用.,62,连接SELs,2.建立新node-nodelinkfromSELAtoSELB,1.删除node-zonelinkonSELA,SELnodesconnecttogridorotherSELnodeswithlinks,63,错误！,连接SELs,64,连接SE