FLCFLC3D基础与应用-陈育民

1、FLAC/FLAC3D基础与应用,陈育民 博士，副教授 河海大学土木与交通学院 Email: ,1,河海大学，2011年11月24日,几个问题,什么是FLAC？ 为什么要用FLAC？ FLAC能做什么？ FLAC为何这么流行？ 怎么学FLAC？,2,什么是FLAC？,Fast Lagrangian Analysis of Continua,3,为什么要用FLAC？,4,Abaqus：745条 Plaxis：80条 ADINA：310条,FLAC能做什么？,岩土工程中的绝大多数问题 土力学、岩石力学、防灾减灾、隧道、地下空间等 采矿工程中的大部分问题 水工结构中的部分问题 结构工程貌似无能为力,

2、国际通用的岩土工程专业分析程序,5,FLAC为何这么流行？,Charles Fairhurst 美国工程院、瑞典皇家工程院院士，国际岩石力学学科和岩石力学学会创始人之一，历任国际岩石力学学会主席和副主席，国际岩石力学学会Muller奖、美国岩石力学学会终生成就奖获得者。 Peter Cundall 美国工程院、英国皇家工程院院士，国际资深计算岩石力学学家。,6,怎么学习FLAC？,7,问：好学不？ 答：好！,问：难不？,答：,入门容易 深入较难,报告安排,第一讲：基本介绍、静力分析、前后处理 第二讲：接触面、FISH语言、流固耦合分析 第三讲：动力分析、自定义本构、结构单元 第四讲：FLAC（

3、2D）基本介绍与应用实例 讨论,8,第一讲,FLAC3D基本介绍、静力分析、前后处理,9,FLAC,大应变、小应变计算模式. 丰富的本构模型、提供自定义的本构模型功能 接触面可以模拟不同材料的接触 流固耦合实现土体的固结与渗流 拥有各种功能的结构单元类型，模拟土与结构的相互作用 强大的动力分析功能. 流变分析，拥有粘弹性模型和粘塑性模型 热力学分析.,10,Shear strainrate contours,FLAC3D,11,与FLAC类似，是FLAC的三维版本 与 FLAC拥有相同的优点,upstream,downstream,可选模块,可选模块包括: 热力学,热-力学耦合,热-流体-力学

4、耦合包括热传导和对流; 粘弹,粘朔性(蠕变)材料模型; 动力学分析,并可以模拟静边界和自由域 使用C+定义自己的模型,12,核废料储存中的热力学研究问题,New Features in FLAC3D Version 4.0,模拟颗粒状材料的硬化模型 自动网格重画功能，解决 bad-geometry 问题. 改进的interface 更快的渗流计算 更新的动力计算功能,13,Lagrangian格式动量平衡方程,14,F(t),m,牛顿运动定律,对于连续体,在静力平衡条件下，加速度项为0，方程变为平衡方程,自由落体的模拟,15,G = mg,S = 1/2gt2 = 20m,命令流： confi

5、g dyn gen zon bri size 1 1 1 ini x mul 0.1 y m 0.1 z m 0.1 model elas prop bulk 3e8 shear 1e8 ini dens 1000 set grav 0 0 -10 solve age 2,自由落体的模拟(movie),16,FLAC3D的求解过程,17,FLAC3D中的本构模型,开挖模型null 3个弹性模型 各向同性弹性 横观各向同性弹性 正交各向同性弹性 8个弹塑性模型 Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型、修正

6、剑桥模型和胡克布朗模型,18,分析问题的过程,19,建立网格,初始条件,边界条件,初始应力平衡,外荷载,求解,前处理,后处理,初始应力的生成,为什么要单独列出？ 分析过程中出现的很多问题都与初始应力是否合理有关 手册中的例子五花八门 是所有后续分析的基础！ 生成方法 弹性求解 更改强度参数的弹塑性求解 设置初始应力的弹塑性求解 存在水压力的初始应力生成 水下建筑的初始应力生成,20,弹性求解,21,gen zon bri size 1 1 2 m elas prop bulk 3e7 shear 1e7 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix

7、 y ran y 0 fix y ran y 1 ini dens 2000 set grav 0 0 -10 solve,Step = 162 sz = -40e3 sx = -21.54e3,更改强度参数的弹塑性求解,22,gen zon bri size 1 1 2 model mohr prop bulk 3e7 shear 1e7 c 1e10 f 15 ten 1e10 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 set grav 0 0 -10 solve

8、prop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0 solve,Step = 163 sz = -40e3 sx = -21.54e3,Or: solve elastic,设置初始应力的弹塑性求解,23,gen zon bri size 1 1 2 model mohr prop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 ini szz -40e3 grad 0 0

9、 20e3 ran z 0 2 ini syy -20e3 grad 0 0 10e3 ran z 0 2 ini sxx -20e3 grad 0 0 10e3 ran z 0 2 set grav 0 0 -10 solve,Step = 0 sz = -40e3 sx = -20e3,存在水压力的初始应力生成(1),24,gen zone brick size 1 1 2 model mohr prop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran

10、y 0 fix y ran y 1 ini dens 2000 ran z 0 1 ini dens 1500 ran z 1 2 ini szz -35e3 grad 0 0 20e3 ran z 0 1 ini syy -22.5e3 grad 0 0 15e3 ran z 0 1 ini sxx -22.5e3 grad 0 0 15e3 ran z 0 1 ini szz -30e3 grad 0 0 15e3 ran z 1 2 ini syy -15e3 grad 0 0 7.5e3 ran z 1 2 ini sxx -15e3 grad 0 0 7.5e3 ran z 1 2

11、ini pp 10e3 grad 0 0 -10e3 ran z 0 1 set grav 0 0 -10 solve,Step = 0 sz = -35e3 sx = -22.5e3,rd = rs nsrf,存在水压力的初始应力生成(2),25,config fluid gen zon bri size 1 1 2 model elas prop bu 3e7 sh 1e7 ini dens 1500 model fl_iso ini fdens=1000 fmod 0 prop por 0.5 set grav 0 0 -10 water dens 1500 water table fa

12、ce 0 0 1, 0 1 1, 1 1 1, 1 0 1 ini pp 10e3 grad 0 0 -10e3 ran z 0 1 ini szz -30e3 grad 0 0 15e3 ran z 1 2 . set fluid off solve,Step = 142 sz = -40e3 sx = -24e3,单元数较少产生的误差,水下建筑的初始应力生成,26,gen zon bri size 1 1 2 model m prop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10 fix z ran z 0 fix x ran x 0 fix x ran

13、 x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 ran z 0 2 ini szz -50e3 grad 0 0 20e3 ran z 0 2 ini syy -40e3 grad 0 0 15e3 ran z 0 2 ini sxx -40e3 grad 0 0 15e3 ran z 0 2 ini pp 30e3 grad 0 0 -10e3 ran z 0 2 app nstress -10e3 ran z 2 set grav 0 0 -10 solve,Step = 0 sz = -50e3 sx = -40e3,前后处理,27,FLAC

14、3D的前处理,28,FLAC3D网格生成的关键特征： FLAC3D 是命令驱动. 使用FLAC3D内置基元进行形状组合可形成复杂网格. 用户自定义FISH 函数可以用来修改基元网格以创建更为复杂的网格. 第三方软件导入.,FLAC3D网格基元,29,Tecplot的后处理,30,Tecplot的云图和等值线,31,Tecplot中云图的效果与flac3d的shade on的效果差不多,但是其出三维等值线的功能是卓越的.,Tecplot的切片功能,32,Tecplot的slice切片功能.和CAD一样,可以任意切剖面出图.最大优点是可以几个剖面同时出图.,33,整体和截面的数据以单元形式存在列表

15、中,可自由选择和组合多个单元出图.,外界模型的导入,复杂网格的生成难度大 接口编写 不同软件之间的精度差异 Group的定义 采用.flac3d的文件格式 .flac3d文件的格式 G 1 1.0e+00 1.0e+00 1.0e+00 Z B8 1 2 3 4 5 6 7 8 ZGROUP Soil 1 2 3 impgrid 变量和函数名不能与FLAC3D、FISH的保留字相冲突，不要采用过于简单的单词，比如a，hist等，这些都与保留字相冲突。 即使程序中存在与保留字相冲突的变量，FLAC3D也不会提供任何提示，所以提醒读者在编制FISH程序时尽量使用较长的、复杂的变量和函数名。,函数与

16、变量,对变量进行赋值时，不能使用当前函数的函数名放在“=”的右边，比如采用下面的定义abcd = abc + 1.0在FISH程序执行时会提出错误，因为这样会形成递归调用，这种调用方式在FISH程序中是不允许的。 变量和函数的作用是全局的，在命令中的任何地方修改变量的值都会立即生效，因此在实际应用中尽量避免不同的函数中含有相同的变量，因为这样可能会造成赋值错误，并难以检查。 在FLAC3D中可以用如下的命令来引用FISH函数和变量 PRINT用于查看函数和变量的数值； HISTORY命令可以对函数和变量的数值进行记录； SET命令用于变量的赋值。,主要语句,选择语句 CASEOF 表达式 默认

17、语句 CASE n1 表达式的值为n1时的语句 CASE n2 表达式的值为n2时的语句 ENDCASE,主要语句,条件语句 IF 条件表达式 THEN ELSE ENDIF,FISH中条件运算符没有“并”、“或”、“否”这样的符号,if aa 1.0 if aa 2.0 执行语句 endif endif,表达“1aa2”的条件,主要语句,循环语句 LOOP var (exp1, exp2) ENDLOOP 或者 LOOP WHILE 条件表达式 ENDLOOP,p_z = zone_head loop while p_z # null ;语句 p_z = z_next(p_z) endloo

18、p,p_gp = gp_head loop while p_ gp # null ;语句 p_gp = gp_next(p_z) endloop,单元遍历语句,节点遍历语句,主要语句,命令语句 COMMAND （FLAC3D命令） ENDCOMMAMD,应用实例：,土体的模量随小主应力变化,def E_modify p_z = zone_head d_k = 704 d_n = 0.38 d_pa = 101325.0 ;/标准大气压 loop while p_z # null sigma_3 = -1.0 * z_sig1(p_z) E_new = d_k * d_pa * (sigma_3

19、 / d_pa) d_n z_prop(p_z,young) = E_new p_z = z_next(p_z) endloop end E_modify,应用实例,获得最大位移的大小及发生位置,def find_max_disp p_gp = gp_head maxdisp_value = 0.0 maxdisp_gpid = 0 loop while p_gp # null disp_gp = sqrt(gp_xdisp(p_gp) 2 + gp_ydisp(p_gp) 2 + gp_zdisp(p_gp) 2) if disp_gp maxdisp_value maxdisp_value

20、 = disp_gp maxdisp_gpid = gp_id(p_gp) endif p_gp = gp_next(p_gp) endloop end find_max_disp print maxdisp_value maxdisp_gpid,FISH的编写习惯,FISH与建模,两个圆形隧道的连接部分,变直径的隧道部分,FISH检查,主要采用PRINT fish命令 查看变量的赋值是否合理，主要检查值为0的函数和变量，因为FISH程序中一般定义的变量都有实际的意义，输出为0的变量很可能是与保留字相冲突的变量 （如a就是apply的保留字） 或者由于编写笔误产生的变量 （如数字0与大写字母O

21、）。,流固耦合分析,很难！,流-固耦合分析(单相流),基本功能 理论框架 计算模式 渗流边界条件，初始条件 单渗流计算及渗流耦合计算,基本功能,渗流各向同性、各向异性 不同的渗流模型和属性 流体压力，涌入量，渗漏量和不渗水边界 抽水井、点源、体积源 饱和渗流可采用显式差分法、隐式差分法 非饱和渗流采用显式差分法 渗流-固体-热的耦合 流体和固体的耦合程度依赖于土体颗粒(骨架)的压缩程度，用Biot系数表示颗粒的可压缩程度。 循环荷载引起的动水压力变化和土体液化,理论框架,准静态Biot理论 多孔介质中遵循Darcy定律的单相渗流 描述多孔介质中流体渗流的变量 孔隙水压力，饱和度，特定排水向量的

22、三个分量 质量守恒定律 达西定律 本构定律 考虑流体响应孔隙水压力改变，饱和度改变，体积应变改变和温度改变,有效应力计算,无渗流模式 不设置CONFIG Fluid 孔压不改变 设置孔压分布 INITIAL pp WATER table WATER density WATER table face SET gravity 手动设置干湿密度,渗流模式 设置CONFIG fluid 设置土体干密度 渗流模型 MODEL fl_isotropic MODEL fl_anisotropic MODEL fl_null SET fluid off set WATER bulk = 0,渗流边界条件，初始

23、条件,默认的边界条件是不透水边界 孔隙压力自由(不透水边界 ) 固定孔隙水压力(透水边界) 如：水井 孔隙压力，孔隙率，饱和度和流体属性的初始分布可以用INITIAL命令或者PROPERTY命令定义。,单渗流计算及渗流耦合计算,时间比例 完全耦合分析方法 孔压固定分析(有效应力分析) 单渗流得到孔压分布 无渗流计算孔压的力学响应 流-固耦合计算,单渗流得到孔压分布,用途：排水沟；抽水井；耦合计算 计算步骤 CONFIG fluid SET mech off SET fluid implicit on/off MODEL fl_; PROP STEP; SOLVE age; SET fluid

24、ratio SET fluid off mech on PROP biot_c 0 (or INI fmod 0),无渗流计算孔压的力学响应,不排水短期响应 两种分析方法：干法和湿法 干法：Ku=K+a2M 两种破坏形式 WATER或INI获得常孔压，不排水的c, (孔压改变较小) =0,c=cu (MK+4/3G) 湿法：耦合体系的短期行为 使用排水的K, c, 若SET fluid off, Biot_mod(fmod)真实,力学过程的特征时间 流体扩散过程的特征时间,完全耦合分析方法,时间比例 短期行为 (不排水) ts(分析时间)tc 施加扰动的属性 流体扰动：渗流可不与力学过程耦合

25、力学扰动：耦合等级取决于流固刚度比 流固刚度比,流-固耦合计算,CONFIG fluid; M(Kf); K(渗透系数) 真实，则FLAC3D默认耦合计算 pevs evp 预估流/力特征时间 耦合计算前先达到一个平衡状态 SET fluid on mech off; SET fluid off mech on; STEP SET mech force; SET mech substep n auto; SET fluid substep m (=1) STEP:渗流步足够小,流固耦合的计算方法,手动调整的STEP求解 SET fluid on mech off STEP SET fluid

26、off mech on STEP 主从进程的SOLVE求解 SET mech force SET mech substep n auto (从进程) SET mech substep m(主进程) SOLVE age 自动STEP求解 STEP,渗流问题(CONFIG fluid)分析步骤,时间比例(ts, tc) 稳态 不排水状态 相当 扰动类型 力学扰动 孔压扰动 流固刚度比Rk 是否1 完全耦合模式 时间比例相当；力学扰动,心墙土坝的渗流 (1),new config fluid set fluid off gen zon brick p0 0 0 -10 size 20 1 10 ge

27、n zon brick p0 5 0 0 p1 15 0 0 p2 5 1 0 p3 9 0 5 p4 15 1 0 p5 9 1 5 p6 11 0 5 p7 11 1 5 size 10 1 5 group soil group dam ran x 5 7 z -5 0 group dam ran id 201 a id 211 a id 221 a id 231 a id 241 a group dam ran id 202 a id 212 a id 222 a id 232 a id 242 a m e prop bu 3e7 sh 1e7 ini pp 0 grad 0 0 -10

28、e3 ran z 0 -10 ini dens 2000 model fl_iso prop por 0.5 perm 1e-10 ini fden 1000 ften -1e10 ini sat 0.0 ran z 0 5 model fl_null ran gro dam ;ini pp 0 ran gro dam fix z ran z -10 fix x ran x 0 fix x ran x 20 fix y set grav 10 solve save elastic.sav,网格模型,初始孔压,心墙土坝的渗流 (2),rest elastic.sav ini xd 0 yd 0

29、zd 0 xv 0 yv 0 zv 0 app nstress -40e3 grad 0 0 10e3 ran z 0 4 x 0 9 solve save pressure.sav,竖向应力,沉降,心墙土坝的渗流 (3),rest pressure.sav set fluid on mech off ini fmod 2e3 ften 0.0 ran gro soil ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0 app pp 40e3 grad 0 0 -10e3 ran z 0 4 x 0 9 app pp 0 ran z 0 x 15 20 hist id=10

30、zone pp id 215 solve,Ratio = 1,Ratio = 1E-5,荷载引起的地基土体的超孔隙水压力,计算文件,; - apply load slowly - def ramp ramp = min(1.0,float(step)/200.0) end apply nstress = -40e3 hist ramp range x -.1 3.1 z 9.9 10.1 ; - fluid flow model - model fl_iso ini fmod 2e9 ; - pore pressure fixed at zero at the surface - fix pp

31、 0 range z 9.9 10.1 ; - settings - set fl off ; - test - step 7501 1因为本例中没有设置初始应力，这里只进行了750步的求解。,计算结果,对主从进程法的讨论,SET mech force 设置一个不平衡力的大小，达到这个不平衡力系统认为暂时达到平衡状态； SET mech substep n auto 设置力学进程为从进程，在主进程每执行一步中必须执行n步，当系统达到平衡时也可以少于n步； SET fluid substep m 设置流体进程为主进程。,讨论1：对收敛准则进行对比分析 不平衡力（force）：1E3、5E3、1E

32、4、5E4 不平衡力比（ratio）：1E-4、1E-3、1E-2、1E-1,采用1E-3的收敛准则既可以满足计算流固耦合过程中的计算精度要求，同时又具有较高的计算效率,对主从进程法的讨论(2),讨论2：子步数的影响,set mech sub 100 fluid sub 10 set mech sub 10 fluid sub 10 set mech sub 1 fluid sub 1,设置合理的子步数也很重要。设置过大，则会导致计算时间大大增加，过小又会造成计算结果的误差。,真空预压的简单模拟,孔压边界条件 tstc 长期分析(排水) Rk1 骨架很软 孔压扰动 进行biot_mod调整,砂

33、层,软土层,粘土层,PVD,2m,8m,10m,Data file:,数值分析过程(movie),FLAC3D非线性动力分析,76,非常复杂！,Said by Prof. Peter Cundall,为什么要用FLAC做动力分析？,FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。 可以计算塑性引起的永久变形以及孔隙水压力的消散。 土动力学中常用的等效线性方法无法直接处理上述问题。,77,动力模拟的3个重要问题,动力荷载与边界条件 材料响应与阻尼 土体液化,78,动力荷载,动力输入的类型 加速度时程 速度时程 应力(压力)时程 力

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

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

36、attach的边界将不能设置FF边界 边界上的Interface将不能连续 动力边界设置需在FF边界设置之前,83,Free-field边界与动力荷载,84,模型底部边界 fix施加速度或加速度荷载刚性边界 Free施加应力时程荷载柔性边界 对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载,Note that there is a factor of 2 because the input energy divides into a downward- & upward-propagating wave.,2. 材料响应与阻尼,85,连续的非线性，表观模量随着应变的增大而降低 对

37、所有循环应变等级均存在滞回特性，因此导致随着循环应变的增加阻尼比增大。阻尼是率相关的。 对于复杂波形的各个成分都产生阻尼。 剪切应变会产生的体积应变，相应的，随着剪应变循环次数的增加体积应变逐渐积累。,材料响应,86,土体在循环荷载作用下呈现出模量衰减和能量消散的特点，那么如何用非线性数值方法对其进行模拟呢？,Nonlinear characteristics of soils (Martin and Seed, 1979),试验得到的阻尼比、割线模量随循环剪应变的曲线,87,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-Ran

38、ge 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,等效线性方法,88,等效线性方法是岩土地震工程中模拟波的传播的最常用的方法。 假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。,Iteration toward strain-compatible shear modulus and damping ra

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

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

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

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

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

44、 curve reaches minimum at:,中心频率的选择,97,假设弹性模型计算土石坝不同材料的功率谱曲线 (assuming elastic material),frequency,frequency,frequency,frequency,frequency,阻尼比的选择,98,根据65%最大应变来选择阻尼比和模量衰减参参数,瑞利阻尼的缺点,必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率） 刚度比例项会导致时间步减小,99,滞后阻尼,100,FLAC Version 5.0 and FLAC3D Version 3.0 提供了滞后阻尼功能，该项

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

46、,102,可以直接采用动力试验中的模量衰减曲线； 相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步； 可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。,一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别,“Good” fit to Seed & Idris data for G/Gmax (sigmoidal 3-parameter function) note inconsistent damping result.,G/Gmax,D - % of critical,滞后阻尼的比较,滞后的默认模型,滞后的三参数模型,滞后的四参数模型,滞后的哈丁模型,滞后阻尼,低循环应变下得到的阻尼比

47、要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个小量的Rayleigh阻尼(0.2%刚度比例)，这样也不会降低时步； 若初始应力不为0，剪应力-剪应变曲线可能不匹配。因此在生成初始应力时就要调用Hyst阻尼； Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大； Hyst阻尼之前要做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的； 即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服

48、模型，不用Hyst阻尼,104,动孔压的生成液化,干沙剪应变循环加载试验 初始加载阶段，沙土通常先压实再膨胀。卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在。取决于初始孔隙率，这可能代表纯粹的压实 假定孔隙中充满水 对于常体积测试，有效应力降低，孔隙水压保持不变 对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小 有效应力为零时发生液化,105,动孔压的生成液化,因此孔隙水压增加不是液化的基本原因 由于颗粒间 (重组以后) 的低接触力导致有效应力的减小 描述液化的模型 高级模型：BSHP (边界面低塑性本构模型, Wang et al. 1990)

49、简单模型：MC + 体积应变增量模型 Finn模型： Byrne模型：,106,UBCTOT: Post-liquefaction,107,Assign post-liquefaction properties,t,g,t,PL-Finn,108,0：零有效应力状态 1：非零有效应力状态,0有效应力,非0有效应力,地震波的调整,基线校正 对于地震分析的加速度时程，其积分得到的速度和位移应归0 美国地质调查研究所 Basic Strong-Motion Accelerogram Processing Software (BAP) 对网格施加一个固定速度从而使残余的位移变为0 动力荷载的频率与单元

50、尺寸的双向调整 高频的输入要求单元尺寸很小 一定的单元尺寸对应输入的最大频率 一般进行滤波处理 滤掉低能量的高频 FFT.FIS Origin SeismoSignal,109,地震波的调整,110,El-Centro波,FFT,修正后的时程,修正后FFT,5Hz,PL-Finn模型的应用,阪神地震的码头分析,111,沉箱顶部水平残余位移最大达5m，平均为3.5m，残余沉降为1 2m，海侧倾斜角3 5 o。,sea Silt replaced sand stone land Sand caisson clay,FLAC3D grid 8050 zones 10386 gps,可液化的砂土,分析

51、参数,112,基本力学参数,流体参数,动力参数,PL-Finn模型参数,输入加速度,113,南北向地震波 最大值 = 0.6 g 竖直向地震波 最大值 = 0.2 g,接触面与自由边界条件,114,模型周围自由场 边界条件,计算结果,位移,115,震后水平位移等值线,沉箱顶点处的水平位移、沉降 时程曲线,3.44 m,1.99 m,计算结果（续）,超孔压比,116,震后的超孔压比云图,超孔压比时程曲线,A: 置换砂 B: 海砂（底部） C: 海砂（中部）,小结,FLAC的动力反应分析非常复杂，用户在进行动力分析之前要具有熟练的静力分析基础 FLAC的动力计算需要花费很多时间，主要原因是由于有限

52、差分法本身需要较小的时间步，因此建议在进行动力分析前首先进行简单模型的分析 “小例子”虽然会额外花费一定时间，但对于熟悉边界条件、荷载施加、阻尼选择等内容来说非常必要,自定义本构模型,118,自定义本构模型的基本方法,必要性 试验总结的本构模型 特定条件下的本构模型 交叉学科的本构模型 二次开发环境 自定义本构模型的功能 自定义本构模型的基本方法,119,二次开发环境,FLAC3D采用面向对象的语言标准C+编写 本构模型都是以动态连接库文件(.DLL文件)的形式提供 VC+2005或更高版本的开发环境 优点 自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平 自定义本构模型(.DLL

53、文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中,120,自定义本构模型的功能,主要功能：对给出的应变增量得到新的应力 辅助功能： 模型名称、版本 读写操作 模型文件的编写 基类(class Constitutive Model)的描述 成员函数的描述 模型的注册 模型与FLAC3D之间的信息交换 模型状态指示器的描述,121,自定义本构模型的基本方法,头文件(usermodel.h)中进行新的本构模型派生类的声明 修改模型的ID(100)、名称和版本 修改派生类的私有成员 C+文件(usermodel.cpp)中修改模型结构 (UserModel:UserMo

54、del(bool bRegister): ConstitutiveModel) const char *UserModel:Properties()函数 模型的参数名称字符串 const char *UserModel:States()函数 计算过程中的状态指示器,122,自定义本构模型的基本方法,double UserModel:GetProperty()和void UserModel: SetProperty()函数 const char * UserModel:Initialize()函数 参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值 const char * Use

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

56、例子(Duncan-Chang),125,其他成功的例子,南京水科院双屈服面模型,126,双屈服面,中心截面沉降云图,水平位移云图,小主应力云图,FLAC基本介绍与应用实例,127,FLAC,GIIC基本操作 Interface的建模 流体计算,128,GIIC基本操作,差分网格的建模思路 成层土建模 养成Mark的好习惯 修剪模型 撤销的最好方法,129,模型草图,计算网格数量,确定模型生成 方法,不规则形状生成,Attach, interface,参数化,运行FLAC,Build选项卡,130,Alter,131,Material,132,In-situ,133,Structure,134,Beam，梁单元可以用于各种类型的支护模拟，包括开挖支护、隧道中的支架等。可以在其两侧连接Interface单元来模拟岩土介质中的挡土墙，还可以通过Interface单元与FLAC网格相连以模拟土工格栅。 Liner，衬砌单元主要用于隧道衬砌，包括混凝土初衬或喷射混凝土初衬。 Cable，锚索单元不能承受弯矩，可以施加预应力，常用于模拟受拉构件，包括岩石中的锚杆。 Pile，桩单元常用于模拟地基中的桩。 Rockbolt，岩石锚杆单元常用于模拟岩石中的锚杆。 Strip，条形锚单元用来模拟加固提防或土坝中的多层条带型结构。 Support，支撑单元用来模