FLAC3D基础与应用(全部)

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

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

3、与应用实例 讨论,9,第一讲,FLAC3D基本介绍、静力分析、前后处理,10,软件介绍,Fast Lagrangian Analysis of Continua 美国Itasca咨询公司开发2D程序(1986) 1990年代初引入中国 有限差分法(FDM) 3D版本：DOS版2.0 2.1 3.0 3.14.0 2D版本：DOS版4.05.06.07.0,11,FLAC,大应变、小应变计算模式. 丰富的本构模型、提供自定义的本构模型功能 接触面可以模拟不同材料的接触 流固耦合实现土体的固结与渗流 拥有各种功能的结构单元类型，模拟土与结构的相互作用 强大的动力分析功能. 流变分析，拥有粘弹性模型

4、和粘塑性模型 热力学分析.,Shear strainrate contours,12,FLAC3D,与FLAC类似，是FLAC的三维版本 与 FLAC拥有相同的优点,upstream,downstream,13,基本特点,内置材料模型 连续介质非线性，大应变模拟 显式解题方案，为不稳定物理过程提供稳定解 界面或滑动面用来模拟可产生滑动或分离的离散面，从而模拟断层，节理或摩擦边界 内置材料模型丰富： 零模型, 三个弹性模型 (各向同性，横观各向同性和正交各向异性), 八个朔性模型 (德鲁克-布拉格, 摩尔-库伦, 应变硬化/软化，单一节理，双线性应变硬化/软化单一节理, 双屈服，修正剑桥粘土，霍

5、克-布朗),14,可选模块,可选模块包括: 热力学,热-力学耦合,热-流体-力学耦合包括热传导和对流; 粘弹,粘朔性(蠕变)材料模型; 动力学分析,并可以模拟静边界和自由域 使用C+定义自己的模型,核废料储存中的热力学研究问题,15,FLAC/FLAC3D基本原理,FLAC/FLAC3D利用有限差分，显示方案，动态松弛方法模拟连续体的非线性力学行为： 即使对准静态问题，程序仍然求解完整的动力学方程。这种方法的好处在于可以为物理非稳定过程例如塌方提供稳定解； 在 “松弛”方法中，使用阻尼来吸收动能以模拟系统的“静态”反应。 这种方法可以用比其它方案如解矩阵法更为真实有效地模拟塌方问题。,16,L

6、agrangian法,源自流体力学中的拉格朗日法 跟踪流体质点的运动状态 跟踪固体力学中结点，按时步用Lagrangian法研究网格节点的运动 节点和单元随材料移动，边界和接触面与单元的边缘一致 固体力学大变形理论,法国数学家、物理学家拉格朗日,17,混和离散技术,FLAC混和离散,+,/2,=,每个,为常应力/应变:,体积应变由整个四边形算出,. 应变偏量则有两个三角形,和,分别算出,(混合离散 过程),解题过程中网格坐标按照“拉格朗日方式更新” (网格随材料移动), 且为显式 (一个时步内局部变化不会影响邻域),18,FLAC3D混和离散,+,/2,=,19,FLAC3D混和离散,结构域离

7、散为可由四面体单元组合形成的五面体或六面体等单元; 以 为基本单元(常应力、常应变); 体应变的计算： ; 偏应变的计算： .,20,动态松弛,动态松弛法 在动态松弛法中，网格点根据牛顿运动定律运动. 网格点的速度与该点的不平衡力呈正比. 这种求解方法所决定的一系列位移将把系统带入平衡状态，或表明破坏模式. 在动态松弛法中有两个因素很重要: 时步的选择 阻尼效应,21,显式算法,显式解与隐式解的比较,显式,逐时推进,隐式, 静态,1.无需进行反复迭代来实现非线性本构关系 . 2. 类似问题求解时间呈 N3/2 规律增长 3. 物理非稳定性不会引起数值不稳定性. 4. 因为无需储存矩阵，用较小内

8、存即可模拟大尺度问题. 5.对大位移、大应变问题同样适合，无需额外的计算 .,1.需进行反复迭代来实现非线性本构关系 2.类似问题求解时间呈 N2 甚至 N3规律增长. 3.难以模拟物理非稳定性问题. 4.需存储刚度矩阵，需克服相关的带宽问题，需要的内存较大 . 5.对大位移、大应变问题需进行大量的计算 .,22,New Features in FLAC Version 6.0,使用Intel Fortran compiler拥有更快的计算速度 自动网格重画功能，解决 bad-geometry 问题. 新的模拟颗粒土材料的硬化模型 更新的通用网格生成工具,23,New Features in

9、FLAC3D Version 3.1,多处理器的并行计算功能 新结构单元类型 “Embedded Liner” 提供两个方向的接触作用，可以很好地模拟挡土墙 对四面体单元采用新的混合离散方法 “Nodal Mixed Discretization” 提供塑性问题更精确的解答 64位程序 包含命令手册、FISH手册和应用实例的帮助,24,New Features in FLAC3D Version 4.0,模拟颗粒状材料的硬化模型 自动网格重画功能，解决 bad-geometry 问题. 改进的interface 更快的渗流计算 更新的动力计算功能,25,Lagrangian格式动量平衡方程,F

10、(t),m,牛顿运动定律,对于连续体,在静力平衡条件下，加速度项为0，方程变为平衡方程,26,自由落体的模拟,G = mg,S = 1/2gt2 = 20m,命令流： config 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,27,自由落体的模拟(movie),28,FLAC3D中模型术语,节点,gridpoint：节点 zone：单元 boundary：边界,29,FL

11、AC3D的求解过程,30,FLAC3D中的本构模型,开挖模型null 3个弹性模型 各向同性弹性 横观各向同性弹性 正交各向同性弹性 8个弹塑性模型 Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型、修正剑桥模型和胡克布朗模型,31,FLAC3D中的本构模型,32,一个最简单的例子,gen zon bri size 3 3 3 ;建立网格(前处理) model elas ;材料参数 prop bulk 3e6 shear 1e6 ini dens 2000 ;初始条件 fix z ran z -.1 .1 ;边界

12、条件 fix x ran x -.1 .1 fix x ran x 2.9 3.1 fix y ran y -.1 .1 fix y ran y 2.9 3.1 set grav 0 0 -10 solve ;求解 app nstr -10e4 ran z 3 x 1 2 y 1 2 solve plo con zd ;后处理 切片功能,RUN FLAC3D,33,分析问题的过程,建立网格,初始条件,边界条件,初始应力平衡,外荷载,求解,前处理,后处理,34,FLAC3D的文件格式,保存文件 (*.sav) 含有所有状态变量和用户定义条件的二进制文件 数据文件 (*.dat) 数据文件由用户创

13、建的一种ASC格式的文件，它包括一系列的用于描述所分析问题的FLAC3D命令 FISH文件(*.fis) FISH程序文件 FLAC3D文件(*.flac3d) FLAC3D的网格信息文件 历史记录文件 (*.his) 记录输入输出历史值的文件 图形文件 图形文件(各种标准格式) 电影文件 (*.dcx) AVI或PCX图像文件，这些图像文件可以当作电影放映,35,初始应力的生成,为什么要单独列出？ 分析过程中出现的很多问题都与初始应力是否合理有关 手册中的例子五花八门 是所有后续分析的基础！ 生成方法 弹性求解 更改强度参数的弹塑性求解 设置初始应力的弹塑性求解 存在水压力的初始应力生成 水

14、下建筑的初始应力生成,36,弹性求解,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 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,37,更改强度参数的弹塑性求解,gen zon bri size 1 1 2 model mohr prop bulk 3e7 shear 1e7 c 1e10 f 1

15、5 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 prop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0 solve,Step = 163 sz = -40e3 sx = -21.54e3,Or: solve elastic,38,设置初始应力的弹塑性求解,gen zon bri size 1 1 2 model mohr prop bulk 3e7 shear 1e7 c 10e3

16、 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 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,39,存在水压力的初始应力生成(1),gen zone brick si

17、ze 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 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

18、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 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,40,存在水压力的初始应力生成(2),config fluid gen zon bri size 1 1 2 model elas prop bu 3e

19、7 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 face 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,单元数较少产生的误差,41,水下建筑的初始应

20、力生成,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 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

21、 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,42,前后处理,43,前后处理,基本前后处理 命令操作 菜单操作 dd 100GPa(岩石) 反分析方法：通过断层中岩石的变形与原岩的变形,单桩承载力分析,软土地基 bulk 1.6878E6 shear 3.6167E5 coh 15E3 fric 12 dens 1.73E3 桩体 bulk 5e9 shear 3.75e9 dens 2.5e3,0.5m,8m,10m,20m,计

22、算过程,施加桩顶荷载,计算结果,影响因素,7.50E+03,fric,g,k,3,coh/0.7,0.7fric,10g,10k,2,3.00E+04,20,100e,100e,1,coh,fric,ks,kn,水平因素,计算工况设计,1,2,3,3,(9),3,1,2,3,(8),2,3,1,3,(7),2,1,3,2,(6),1,3,2,2,(5),3,2,1,2,(4),3,3,3,1,(3),2,2,2,1,(2),1,1,1,1,(1),coh,fric,ks,kn,水平因素 工况,计算结果,Ks取1,Ks取2,Ks取3,96.9%,1,2,3,3,(9),22.4%,3,1,2,3

23、,(8),54.1%,2,3,1,3,(7),96.9%,2,1,3,2,(6),7.1%,1,3,2,2,(5),37.8%,3,2,1,2,(4),96.9%,3,3,3,1,(3),42.9%,2,2,2,1,(2),37.8%,1,1,1,1,(1),Ra差异度,coh,fric,ks,kn,水平因素,最优方案,2,2,2,2,11,1,1,2,2,10,1,1,2,2,最优方案,17.3%,6.8%,72.8%,11.9%,极差,52.4%,52.7%,96.9%,57.8%,k3,64.6%,59.2%,24.1%,47.3%,k2,47.3%,52.4%,43.2%,59.2%,

24、k1,coh,fric,ks,kn,水平因素,合理步骤,单桩分析,简单网格,接触面参数,多次试算,理论、实测,加密网格,接触参数,理想结果,群桩分析,Pile结构单元,单元参数,理想结果,Pile结构单元,单元参数,理想结果,挡土墙的接触面设置,对于未知问题的分析思路 3个独立的接触面 3个同ID的接触面 2个独立的接触面 2个独立的接触面并进行底部merge,挡墙,土体,1,2,3,Wall,Soil,挡土墙的接触面设置(2),Z = 5.16cm,Z = 1.14cm,Z = 不能初始平衡,Z = 0.35cm,？,挡土墙的接触面设置(3),X = 1.95cm,X = 1.91cm,X

25、= 1.47cm,？,挡土墙的接触面设置(3),前两种方法的差别实质 3个ID的独立接触面在相同位置产生互不影响的两个节点 共同ID的接触面在相同位置自动设置为1个节点 最终的结论需要您自己去判断！,3 interfaces, 2 IDs,1 interface, 1 ID,FISH语言,FISH语言简介,软件自带的编程语言 是否一定要学？ 视情况而定，需要时查询FISH变量即可 语法简单 xxx end_xxx 注意事项 与FLAC本身的关键字冲突 保留字不可缩写 变量可不定义，因此注意检查程序 print fish table, extra等命令使用,一个最简单的FISH程序,def ab

26、c abc = 1 + 2 * 3 abcd = 1.0 / 2.0 end abc print fish,数据格式,函数与变量,都可以在FISH函数中进行赋值，赋值操作与常规的编程语言类似，按照运算符的优先级先后顺序来执行。 函数和变量的赋值遵守数据类型的规则，即整型的计算结果为整型，浮点型的计算结果为浮点型，因此读者在进行除法运算、开方运算时都需要将数据类型设置为浮点型，数字尽量使用小数点以保证运算正确。 变量和函数名的命名规则是不能以数字开头，不能含有中文，并且不能包含如下的字符。 . , * / + - = # ( ) ; 变量和函数名不能与FLAC3D、FISH的保留字相冲突，不要采

27、用过于简单的单词，比如a，hist等，这些都与保留字相冲突。 即使程序中存在与保留字相冲突的变量，FLAC3D也不会提供任何提示，所以提醒读者在编制FISH程序时尽量使用较长的、复杂的变量和函数名。,函数与变量,对变量进行赋值时，不能使用当前函数的函数名放在“=”的右边，比如采用下面的定义 abcd = abc + 1.0 在FISH程序执行时会提出错误，因为这样会形成递归调用，这种调用方式在FISH程序中是不允许的。 变量和函数的作用是全局的，在命令中的任何地方修改变量的值都会立即生效，因此在实际应用中尽量避免不同的函数中含有相同的变量，因为这样可能会造成赋值错误，并难以检查。 在FLAC3

28、D中可以用如下的命令来引用FISH函数和变量 PRINT用于查看函数和变量的数值； HISTORY命令可以对函数和变量的数值进行记录； SET命令用于变量的赋值。,主要语句,选择语句 CASEOF 表达式 默认语句 CASE n1 表达式的值为n1时的语句 CASE n2 表达式的值为n2时的语句 ENDCASE,主要语句,条件语句 IF 条件表达式 THEN ELSE ENDIF,FISH中条件运算符没有“并”、“或”、“否”这样的符号,if aa 1.0 if aa 2.0 执行语句 endif endif,表达“1K+4/3G) 湿法：耦合体系的短期行为 使用排水的K, c, 若SET

29、fluid off, Biot_mod(fmod)真实,力学过程的特征时间 流体扩散过程的特征时间,完全耦合分析方法,时间比例 短期行为 (不排水) ts(分析时间)tc 施加扰动的属性 流体扰动：渗流可不与力学过程耦合 力学扰动：耦合等级取决于流固刚度比 流固刚度比,流-固耦合计算,CONFIG fluid; M(Kf); K(渗透系数) 真实，则FLAC3D默认耦合计算 pevs evp 预估流/力特征时间 耦合计算前先达到一个平衡状态 SET fluid on mech off; SET fluid off mech on; STEP SET mech force; SET mech s

30、ubstep n auto; SET fluid substep m (=1) STEP:渗流步足够小,流固耦合的计算方法,手动调整的STEP求解 SET fluid on mech off STEP SET fluid 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) 稳态 不排水状态 相当 扰动类型 力学扰动 孔压扰动 流固

31、刚度比Rk 是否1 完全耦合模式 时间比例相当；力学扰动,心墙土坝的渗流 (1),new config fluid set fluid off gen zon brick p0 0 0 -10 size 20 1 10 gen 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

32、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 -10e3 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

33、 0 fix x ran x 20 fix y set grav 10 solve save elastic.sav,网格模型,初始孔压,心墙土坝的渗流 (2),rest elastic.sav ini xd 0 yd 0 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 i

34、ni 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 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

35、.9 10.1 ; - fluid flow model - model fl_iso ini fmod 2e9 ; - pore pressure fixed at zero at the surface - fix pp 0 range z 9.9 10.1 ; - settings - set fl off ; - test - step 7501 1因为本例中没有设置初始应力，这里只进行了750步的求解。,计算结果,对主从进程法的讨论,SET mech force 设置一个不平衡力的大小，达到这个不平衡力系统认为暂时达到平衡状态； SET mech substep n auto 设置力

36、学进程为从进程，在主进程每执行一步中必须执行n步，当系统达到平衡时也可以少于n步； SET fluid substep m 设置流体进程为主进程。,讨论1：对收敛准则进行对比分析 不平衡力（force）：1E3、5E3、1E4、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

37、 sub 1,设置合理的子步数也很重要。设置过大，则会导致计算时间大大增加，过小又会造成计算结果的误差。,真空预压的简单模拟,孔压边界条件 tstc 长期分析(排水) Rk1 骨架很软 孔压扰动 进行biot_mod调整,砂层,软土层,粘土层,PVD,2m,8m,10m,Data file:,数值分析过程(movie),124,FLAC3D非线性动力分析,非常复杂！,Said by Prof. Peter Cundall,125,为什么要用FLAC做动力分析？,FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。 可以计算塑性引起的

38、永久变形以及孔隙水压力的消散。 土动力学中常用的等效线性方法无法直接处理上述问题。,126,动力模拟的3个重要问题,动力荷载与边界条件 材料响应与阻尼 土体液化,127,动力荷载,动力输入的类型 加速度时程 速度时程 应力(压力)时程 力时程 APPLY INTERIOR (内部) TABLE FISH,128,Quiet边界,静态(quiet,粘性)边界 Lysmer and Kuhlemeyer(1969) 模型边界法向和切向设置独立的阻尼器 性能 对于法向p波和s波能很好的吸收 对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射 人工边界仍应当足够远,129,Quiet边界应用,

39、内部振动(如隧道中的列车振动问题) 动力荷载直接施加在节点上 使用Quiet边界减小人工边界上的反射 不需要FF边界 外部荷载的底部边界 软土地基上的地震荷载不适合用加速度或速度边界条件 使用应力条件t = -2Csrvs 地震底部输入的侧向边界 扭曲了入射波,quiet,quiet,quiet,130,Free-field边界,Cundall et al. (1980) 自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上 设置条件 底部水平，重力方向为z向 侧面垂直，法向分别为x, y向 其他边界条件在APPLY ff之前,相当于一个阻尼器,131,Free-fi

40、eld边界,APPLY ff将边界上单元的属性、条件和变量全部转移ff单元上； 设置以后主体网格上的改动将不会被FF边界所响应 可存在任意的本构模型以及流体耦合(仅竖向) FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小 存在attach的边界将不能设置FF边界 边界上的Interface将不能连续 动力边界设置需在FF边界设置之前,132,Free-field边界与动力荷载,模型底部边界 fix施加速度或加速度荷载刚性边界 Free施加应力时程荷载柔性边界 对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载,Note that there is a facto

41、r of 2 because the input energy divides into a downward- & upward-propagating wave.,133,2. 材料响应与阻尼,连续的非线性，表观模量随着应变的增大而降低 对所有循环应变等级均存在滞回特性，因此导致随着循环应变的增加阻尼比增大。阻尼是率相关的。 对于复杂波形的各个成分都产生阻尼。 剪切应变会产生的体积应变，相应的，随着剪应变循环次数的增加体积应变逐渐积累。,134,材料响应,土体在循环荷载作用下呈现出模量衰减和能量消散的特点，那么如何用非线性数值方法对其进行模拟呢？,Nonlinear characteris

42、tics of soils (Martin and Seed, 1979),135,试验得到的阻尼比、割线模量随循环剪应变的曲线,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,136,等效线性方法,等效线性方法是岩土地震工程中模拟波的传播的最常用的

43、方法。 假定土体是粘弹性体，参照实验室得到的切线模量及阻尼比与剪应变幅值的关系曲线，对地震中每一单元的阻尼和模量重新赋值。,Iteration toward strain-compatible shear modulus and damping ratio (after Kramer, 1996),137,等效线性方法的特点,使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变。在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大。对于不同部位不同运动水平的特性存在空间变异性。 不能计算永久变形。等效线性方法模型在加荷与卸荷时模量相同，不能计算

44、土体在周期荷载作用下发生的剩余应变或位移。 塑形屈服模拟不合理。在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”。然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数。因此，塑性屈服的模拟不合理。 大应变时误差大。等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。 本构模型单一。等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性。方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能

45、反映与频率无关的滞回圈。另外，模形是率无关的，因此不能考虑率相关性。,138,完全非线性分析方法,FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程。,139,完全非线性分析方法的特点,可以遵循任何指定的非线性本构模型。如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数。如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变。 采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点。 由于采用了弹塑性模型，因此程序可以自动计算永久变形。

46、 采用合理的塑性方程，使得塑性应变增量与应力相联系。 可以方便地进行不同本构模型的比较。 可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。,140,使用弹塑性模型,附加考虑的因素： 阻尼，对于屈服面以下应力的循环 体积应变积累，是循环周数与幅值的函数 模量衰减，基于平均应变水平的表格,最简单的弹塑性模型往往在描述累计塑性应变方面具有很好的效果，但是对于加速度放大系数的估算上效果不好。,141,弹塑性模型,简单的理想弹塑性本构模型仅仅在发生屈服时才会出现滞回特性,注意： 即使

47、这样粗糙的模型也能够作出连续的阻尼比和模量衰减曲线。,在屈服条件下会产生体积改变，但通常都是剪胀。,142,FLAC3D中的阻尼比,弹塑性本构模型使用瑞利（粘性）阻尼 弹塑性本构模型使用滞后（HD）阻尼. 复杂本构模型拥有连续的屈服应力应变关系和对应的加卸载响应（Wang，UBCSand）,143,瑞利阻尼,瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅。在计算时，假设动力方程中的阻尼矩阵C与刚度矩阵K和质量矩阵M有关：,瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。虽然两个阻尼器本身是与频率有关的，但是通过选取合适的

48、系数，可以在有限的频率范围内近似获得频率无关的响应。,144,瑞利阻尼,frequency,ratio of damping to critical,混合,仅有刚度分量,仅有质量分量,采用叠加的方法得到的阻尼比在较大的频率范围内保持定值 （3:1）,Combined curve reaches minimum at:,145,中心频率的选择,假设弹性模型计算土石坝不同材料的功率谱曲线 (assuming elastic material),frequency,frequency,frequency,frequency,frequency,146,阻尼比的选择,根据65%最大应变来选择阻尼比和模

49、量衰减参参数,147,瑞利阻尼的缺点,必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率） 刚度比例项会导致时间步减小,148,滞后阻尼,FLAC Version 5.0 and FLAC3D Version 3.0 提供了滞后阻尼功能，该项阻尼独立于材料模型之外，实质是包含了一个与应变幅值相关的切线剪切模量乘子,如果割线模量是由衰减曲线来确定的，那么切线模量可以表示为：,From Seed & Idriss (1970),Go = small-strain shear modulus,149,滞后阻尼,表观应变是偏应变，根据先前反转的点进行累计得到的。使用堆栈的

50、数据结构可以描述应变反转现象。,因此，在“迷你滞回圈”中的能量损失与主滞回圈的计算是一样的。,elastic model with hysteresis damping,特定的曲线可以使用滞后阻尼来描述，并在本构模型计算中根据应变的大小获得不同的切线模量。,150,滞后阻尼的特点,可以直接采用动力试验中的模量衰减曲线； 相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步； 可以应用于任意的材料模型，且可以与其它阻尼格式同时使用。,一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别,“Good” fit to Seed & Idris data for G/Gmax (sigmoidal 3-pa

51、rameter function) note inconsistent damping result.,G/Gmax,D - % of critical,151,滞后阻尼,低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下。可以在中心频率上增加一个小量的Rayleigh阻尼(0.2%刚度比例)，这样也不会降低时步； 若初始应力不为0，剪应力-剪应变曲线可能不匹配。因此在生成初始应力时就要调用Hyst阻尼； Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大； Hyst阻尼之前要

52、做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的； 即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服模型，不用Hyst阻尼,152,动孔压的生成液化,干沙剪应变循环加载试验 初始加载阶段，沙土通常先压实再膨胀。卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在。取决于初始孔隙率，这可能代表纯粹的压实 假定孔隙中充满水 对于常体积测试，有效应力降低，孔隙水压保持不变 对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小 有效应力为零时发生液化,153,动孔压的生成

53、液化,因此孔隙水压增加不是液化的基本原因 由于颗粒间 (重组以后) 的低接触力导致有效应力的减小 描述液化的模型 高级模型：BSHP (边界面低塑性本构模型, Wang et al. 1990) 简单模型：MC + 体积应变增量模型 Finn模型： Byrne模型：,154,UBCTOT: Post-liquefaction,Assign post-liquefaction properties,t,g,t,155,PL-Finn,0：零有效应力状态 1：非零有效应力状态,0有效应力,非0有效应力,156,地震波的调整,基线校正 对于地震分析的加速度时程，其积分得到的速度和位移应归0 美国地质

54、调查研究所 Basic Strong-Motion Accelerogram Processing Software (BAP) 对网格施加一个固定速度从而使残余的位移变为0 动力荷载的频率与单元尺寸的双向调整 高频的输入要求单元尺寸很小 一定的单元尺寸对应输入的最大频率 一般进行滤波处理 滤掉低能量的高频 FFT.FIS Origin SeismoSignal,157,地震波的调整,El-Centro波,FFT,修正后的时程,修正后FFT,5Hz,158,PL-Finn模型的应用,阪神地震的码头分析,沉箱顶部水平残余位移最大达5m，平均为3.5m，残余沉降为1 2m，海侧倾斜角3 5 o。,sea Silt replaced sand stone land Sand caisson clay,FLAC3D grid