Docstoc

第五章_1_

Document Sample
第五章_1_ Powered By Docstoc
					第五章 有限差分法(1)FLAC2D

     制作人:
            王 金 安



     北 京 科 技 大 学
1. FLAC简介
 1.1 FLAC特点
 FLAC——Fast Lagrangian Analysis of Continua
 FLAC建立在拉格朗日算法基础上,采用有限差分显式
  算法来获得模型全部运动方程(包括内变量)的时间步长
  解,从而可以追踪材料的渐进破坏和垮落,这对研究采
  矿工程设计是非常重要的。
 FLAC适用模拟计算岩土材料力学行为,特别适合模拟
  大变形和扭曲,包括材料的高度非线性(应变硬化/软化)、
  不可逆剪切破坏和压密、粘弹(蠕变)、孔隙介质的应
  力—渗流耦合、热—力耦合以及动力学问题等。
1.2 岩体力学本构模型

(1) 各向同性弹性材料模型;
(2) 横观各向同性弹性材料模型;
(3) 莫尔-库仑弹塑材料模型;
(4) 应变软化/硬化塑性材料模型;
(5) 双屈服塑性材料模型;
(6) 遍布节理材料模型;
(7) 空单元模型,可用来模拟岩体开挖和开采。
         岩石应变软化/硬化模型

        岩石                    充填体
                       压力 P

1-3

             莫尔-库仑准则
              与应变软化
                              加载路径

                                     卸载路径




                                    体积应变 eV
                断层模型

                                  A ²à
                         N
                             T
断层
                  S ks       kn

            M                            P
     µ¥Ôª           LN                       B ²à
1.3 支护构件与本构模型
(1) 杆单元模型;
(2) 梁模型;
(3) 桩模型;
(4) 群支架模型。
                   预应力砂浆锚杆的模拟体系

         自由段            锚固段




         锚杆                                       砂浆体




          Kb            Kb                   Kb
               m                    m                 m
                   Sb                   Sb
                                                          锚固结点
砂浆剪切刚度                       cs,             cs, 


   砂浆粘结刚度                           砂浆粘聚力,摩擦力
     FLAC举例
 P




(a)巴西试验   (b) 边坡滑移
 (间接拉伸)
2. FLAC理论基础
 2.1 有限差分法
                                           y


             f  h
                    2
                         2 f   
f1  f 0  h          2
                         x     
                                 
             x 0 2           0                                12
                                                                                h

                                                               8   4    5


              f   h  2
                             2 f                   11       3   0    1   9
f 3  f 0  h             2
                             x     
                                     
              x 0 2              0                        7   2    6


                                                                   10

                                                                                    x
                                                           h



  f    f  f3                            f       f  f4
       1                                
                                            y      2
                                                   
  x  0   2h                                    0    2h

 2 f      f  f3  2 f0                 2 f      f  f4  2 f0

 x 2     1
                                          2
                                           y       2
                                                   
        0      h2                               0      h2
2.2 简单力学问题分析

                                    m
                      u
                  u, u, 
                             F(t)
•牛顿运动定律
                             
                            du
               F  ma  m
                            dt
•对于一个连续体,可写成广义形式:

                 du i ij
                  
                           gi
                 dt    x j
式中,  —— 体密度,
    xj —— 坐标矢量 (x,y)
    ij —— 应力张量分量
    gi —— 重力加速度
2.3 应力—应变关系
除了运动定律,连续介质必须服从本构关系——
应力与应变之间的关系
• 对于弹性材料:
                                 2   .          .
        ij :  ij  { ij ( K  G) e kk  2G eij }t
                                 3
式中,δij——Kronecker记号;
    t——时间步;
   G,K——分别是剪切模量和体积模量。

• 一般地,本构关系可表述为:
                                   .
                 ij : M ( ij , eij , k )

 式中                           .        .
                 .    1 u uj
                 eij  [ i      ]
                      2 x j xi
2.4 有限差分基本格式
•在有限差分法中, 每个在前一个公式中导出的量(运动量或应力应变)
由所在网格特定位置处的相关变量代数式所代替。


•代数式是完全显式的;在代数式右侧的所有量都是已知的。在一个计
算时步,FLAC网格每个元素(域或结点)与相邻单元表现出在物理上
是隔绝的。
                        应力

•计算循环的基本要求
选取的时步非常小,乃至在此时步间隔内    输出应力          非线性定律
实际信息不能从一个单元传递到另一个单
元(事实上,所有材料都有传播信息的某
种最大速度) 。
                                      应变
                             输入应变
  2.5 基本显式计算循环
                         对于所有结点
     速度                  平衡方程                                 结点力
                        (运动方程)
                          dui  ij
                           
                                   gi
                          dt   x j
    高斯定理                                                     Fi   ijn j L

                          对于所有单元
                         应力-应变关系
                          (本构方程)                             新的应力
    应变速率

                                     2   .          .
            ij :  ij  { ij ( K  G) e kk  2G eij }t
例如, 弹性                               3
FLAC的网格内部由三角形构成,三角形组合成四
边形单元。推导差分方程的方式如下:

                      ui(b)                             Fi
 a

         b        b                 ni(2)                    ni(1)
                         ΔS
                                            s(2) s(1)
             b
     a                        a
                                  ui(a)

 单元叠加            三角形单元                    结点力向量
                 速度向量
                          f
高斯定理:
           S ni fdS  A xi dA
用于推导任意形状单元的有限差分公式。                              u (i b ) 结点速度
                                                
                                           b
                                                  S

对于多边形,公式为:
                                                  a
                                                        u (i a ) 结点速度
                                                        
            f   1
                  f ni S
            xi  A S

下式用于计算应变增量, eij :

           ui
                    ui( a )  ui(b) n j S
                1
                              
           x j 2 A S
                 1  ui u j 
                          
           eij             t
                 2  x j xi 
                             
一旦计算出全部应力,可以从作用每个三角形边界上
产生的牵引力计算得到结点力。例如:
                                                                                         Fi
      1
  Fi   ij (n j S (1)  n j S ( 2) )
                (1)         ( 2)
                                                                       ni(2)
                                                                                              ni(1)
      2
                                                                                      s(1)
                                                                               s(2)
然后,用“传统”的中间差分
公式获得新的速度和位移:

  .                   .                                          t
                                              Fi
       ( t  t )          ( t  t / 2 )
                     ui
                                                         (t )
  ui
                                                                 m

  .                       .                 .
       ( t  t )                               ( t  t / 2 )
  xi                 x i  ui
                              (t )
                                                                 t   (…大变形模式)
在时间域上的解算方法
                                             数值结点
                       F
           应力                                                           力F
                                
                       u
                                                                        位移 u
                  显式                          x            隐式

 所有单元:
                                                             单元
  F  f u, 
       (非线性定律)
                            假设 (u)固定
                                                         F  Ku
 所有结点:                      假设 (F)固定                        结点
   u   F t
       
          m
                                                     m Ku  F
                                                         u

 重复n 个时步                     改正,若
                                    x min          每个时步内解算整套方程
                             t 
                                     Cp
 时步内不作迭代
                           p-波速度                    如果存在非线性,在时步
                                                    内迭
  信息在每个时步内
  不在单元之间传播
                  方法比较

    显式, 时间-进程                隐式, 静态

1. 可以遵循非线性定律而无需内部迭      1. 整个迭代过程需要遵循非线性定律。
代, 因为本构计算中位移被“冻结”
。

                        2. 解算时间增加 N2 甚至 N3。
2. 对同样问题,计算时间增加N3/2 。

                        3. 模拟物理不稳定性困难。
3. 物理不稳定不会引起数值不稳定。

                        4. 需要大内存,或大容量硬盘存储。
4. 因为无需存储矩阵,用少量内存
   可以模拟大型问题。

                        5. 模拟大变形明显需要更多的机时。
5. 大应变、大位移和转动模拟无需
   额外机时。
          动态衰减
•   在动态衰减过程中,结点按照牛顿定律运动。结点的
    加速度与非平衡力成正比。该算法确定这组位移将导
    致系统趋于平衡,或表征失稳状态。


•   关于动态衰减的两个重点:
① 选择时步
② 阻尼作用
                  时步
为满足数值稳定,时步必须满足::

                    xmin
               t 
                     Cp
式中, Cp 与 1 /mgp成正比。对于静态分析,结点质量按比例度
量 ,致使局部临界时步等于( t=1),这样提供了收敛的最优
时步。然后,对结点惯性质量进行调整来实现稳定性条件:


                 ( K  4G / 3) x m ax
                                  2
        m pg   
                        3A

注意:重力质量不受影响。
                     阻尼
 •   速度-比例阻尼导致影响计算的体积力。
 •   FLAC中采用了局部阻尼——结点处的阻尼力与不平衡力大小成正
     比,正负号的规定确保振动状态被阻滞。
 局部阻尼
• 在运动方程中引进阻尼力
                                                    t
                 ui  Fi  | Fi | sgn (ui )
                                           
 式中,Fi是不平衡力                                         m

• 阻尼力 Fd 为:
                  Fd   Fi sgn( ui )
                                   

• 在 FLAC 中,通过监测不平衡力比率(不平衡力 Fi 与施加力 Fm之比)
来确定静力状态

•当 Fi / Fm < 0.001时,则(计算机缺省)认为模型处于平衡状态。
1.6   FLAC计算流程举例

 如图,一维杆件用数个等尺寸的有限差分网格划分,杆件的
  密度为,杨氏模量为E。
                            单元                  速度、位移                           应力


           ①       ②        ③       ④                       i-1       i         i+1

       1       2        3       4       5             i-1         i       i+1         x

                   结点

对于固体材料,微分形式的本构方程为:
                                               u x
                                 xx        E                                            (1)
                                               x
 运动方程(或平衡方程)为:
                       2u x  xx
                      2                                  (2)
                       t     x

 假设杆件无侧向约束。对于单元i,与(1)式对应的中间有限
  差分公式为:
                                u x1 (t )  u x (t )
                                  i            i
                   xx (t )  E
                    i

                                         x                 (3)
  式中,(t)表示其变量是在时刻t确定的,上标i表示单元或结点编号。

 对结点i的有限差分形式运动方程为:
       . i      t   .       t    1
        {u x (t  )  u x (t  )}     { xx (t )   xx1 (t )} (4)
                                                      i
                         i                i

     t          2            2     x
          .   t     .       t    t
 或写成 u x (t     )  u x (t  )      { xx (t )   xx1 (t )} (5)
                                                     i
              i         i                i

              2              2    x
 积分后得出位移:
             t              .       t
       u (t  )  u x (t )  u x (t  )t
                                i
        i
        x
                    i
                                            (6)
             2                       2

 在显式算法中,所有有限差分方程右端的值均是已知的。因
  此,必须先用(3)式算出所有单元的应力,然后再由(5)式和
  (6)式计算所有结点的速度和位移。
 在概念上,这个过程等价于对变量值“同时”更新,而不是
  像其它方法那样,方程右端混存有 “新”、“旧”值,对变
  量值“依次”更新。
3. FLAC软件应用
3.1 建模

(1)设计计算模型的尺寸
(2)规划计算网格数目和分布
(3)安排工程对象(开挖、支护等)
(4)给出材料的力学参数
(5)确定边界条件
(6)计算模拟
FLAC 基本称谓




            单元编号




            结点编号
3.2 网格生成 Grid i,j
    例如:grid 30,20

     y




  j=20


                    x

             i=30
3.3 网格规划:Gen x1,y1 x2,y2 x3,y3 x4,y4
    例如:Gen 0,0 0,10 10,20 20,0
   2(x2,y2)                   3(x3,y3)




  1(x1,y1)                      4(x4,y4)
3.4 分区规划网格。例如:
Gen xI1,yI1 xI2,yI2 xI3,yI3 xI4,yI4 i=1,10 j=1,21 (I区)
Gen xII1,yII1 xII2,yII2 xII3,yII3 xII4,yII4 i=10,20 j=1,21 (II区)

                       i=10, j=21               i=21, j=21



       i=1, j=21
                                       II
                   I




        i=1, j=1        i=10                    i=21, j=1
3.5 特殊形状的网格
(1)圆形 gen circle xc,yc rad



              rad

             xc,yc
3.5 特殊形状的网格
(2)弧线 gen arc xc,yc xb,yb theta


                       xb,yb
                

               xc,yc
3.5 特殊形状的网格
(3)直线 gen line x1,y1 x2,y2


                  x2,y2




         x1,y1
(4)任意形状
tab 1 x1,y1, x2,y2, ,xn,yn, x1,y1
gen tab 1

                                (x3,y3)

                 (x2,y2)                     (x4,y4)




       (x1,y1)                            (x5,y5)
                           (x6,y6)
 3.6 赋给单元材料性质

mod e (弹性)
prop d 1800e-6 bu 12.5 sh 5.77 i=1,20 j=1,10
prop d 2400e-6 bu 1250 sh 577 i=1,20 j=11,20

mod m (弹塑性Mohr-Coulumb准则)
prop d 1800e-6 bu 12.5 sh 5.77 c 0 fri 20 ten 0.015
  reg i,j
                           *特别提示
(1)量纲统一

                      SI                             Imperial
Length       m         m          m         cm          ft           In
Density    kg/m3    103kg/m3   106kg/m3   106g/cm3   slugs/ft3   snails/in3
Force        N        kN         MN       Mdynes        Ibf         Ibf
Stress      Pa        kPa        MPa        bar       Ibf/ft2       psi
Gravity    m/sec2   m/sec2     m/sec2      cm/s2     ft/sec2      in/sec2
Stiffnes
           Pa/m      kPa/m     MPa/m      Bar/cm      Ibf/ft3      Ib/in3
s

                              E
(2)参数换算                K                    (bulk mod ulus)
                          3(1  2  )
                             E
                       G                     ( shear mod ulus)
                          2(1   )
  3.7 赋给模型边界条件


(1)固定边界 (结点)
Fix x i=1, j=1,21
Fix y i=1,21 j=1


               Fix x


                i=1

                       Fix y   j=1
3.7 赋给模型边界条件
(2)施加边界力 (结点)
  apply yf -10 i=1,21 j=21
或 apply syy -10 i=1,21 j=21
  apply xf -5 i=21, j=1,21
或 apply sxx -5 i=21, j=1,21
                              j=21




                              i=21
  3.7 赋给模型边界条件

(3)赋单元内应力 (单元)
ini sxx -10 i=1,20 j=1,20
ini syy -5 var 0 4 i=1,21 j=1,21

                                   j=21




                                     i=21
3.8   计算

  Set grav 9.81
  Set large
  Step 1000
  Save test.sav
 3.9 结果显示

Plot grid 显示网格
Plot bo 显示边界
Plot plas 显示塑性区
Plot sig1 fi 显示最大主应力1
Plot sig2 fi 显示最小主应力2
Plot sdif fi 显示主应力差(1- 2)
Plot str 显示主应力矢量场
Plot xdis fi 显示X方向位移
Plot ydis fi 显示Y方向位移
Plot disp 显示位移矢量场
 3.10 保存与调用结果命令

Call test.txt (或ca test.dat) 调用数据

Save test.sav 保存结果

New 重新开始

Rest test.sav 调用结果

Quit 退出程序
4. FLAC解题技巧
4.1 模型尺寸
4.2 模拟开挖
Mod nu i=6,15 j=5,12 (或 region i,j)



       i=6,15
       j=5,12
 4. FLAC解题技巧
4.3 模拟锚杆支护(端锚)
struct cable begin grid i1,j1 end grid i2,j2 seg n prop 1
stru prop 1 e 2e5 yield 0.5 a 0.235e-3 sbond 0 kbond
   0.0001
stru prop 1 sfri 30 peri 0.2723 den 7.5e-3


         i1,j1
                                端锚
                        i2,j2
4.3 模拟锚杆支护(全长锚固)
 struct cable begin grid i,j end x,y seg n prop 2
 stru prop 1 e 2e5 yield 0.5 a 0.235e-3 sbond 0.42
    kbond 5.37
 stru prop 1 sfri 30 perimeter 0.2723 den 7.5e-3



          i,j
                          x,y
                        全长锚
4.3 模拟锚杆支护(预应力锚固)
struct cable begin grid i1,j1 end i2,j2 seg n ten 10 prop 1
struct cable begin node i2,j2 end x,y seg n prop 2
stru prop 1 e 2e5 yield 0.5 a 0.235e-3 sbond 0. kbond 0.1
stru prop 1 sfri 0 perimeter 0.2723 den 7.5e-3
stru prop 2 e 2e5 yield 0.5 a 0.235e-3 sb 0.42 kb 5.37
stru prop 2 sfri 30 perimeter 0.2723 den 7.5e-3

                 自由段
         i1,j1      1        锚固段

                        i2,j2    2
                        node n       x,y
4.4 模拟断层

Mod nu j=37
int 1 as from i1,j1 to ,i1,j2 bs from i2,j1 to i2,j2
int 1 kn 1000 ks 400 fri 25 c 0.01 t 0.0001


                    i1,j2       i2,j2


                         as    bs



                 i1,j1        i2,j1
4. 5 其它技巧

(1)网格优化
(2)先弹性、后塑性
(3)分步开挖
(4)动力学问题模拟
(5)固流耦合问题模拟
(6)合理解释结果
(7)多种形式输出结果
(8)报告格式
返回目录

				
DOCUMENT INFO
Categories:
Tags:
Stats:
views:26
posted:12/11/2011
language:
pages:52