Page 122 - 2024年第55卷第1期
P. 122

在不计体力和裂缝面荷载情况下,由虚功原理可导出多边形比例边界有限元控制方程                                         [15,19] 为:
                                                { {u( ξ )} }  =- [Z] { {u( ξ )} }
                                               ξ                                                        (1)
                                                 { q( ξ )}
                                                         , ξ      { q( ξ )}
              式中:{ u( ξ )}为节点位移函数;{q( ξ )}为节点力函数;[Z]为 Hamilton矩阵                    [19] 。
                  采用特征值分解或者 Schur分解            [18] 均可求解控制方程(1)。以特征值分解为例,多边形内任一点
              的位移{u( ξ ,η )}可表示为节点位移函数{u( ξ )}沿环向边界的插值,具体表达式为:
                                                                           (u)  〈 λ b 〉
                                    { u( ξ ,η )} = [N( η )]{u( ξ )} = [N( η )][ Φ b  ] ξ  {c}           (2)
                       (u)
              式中:[ Φ b    ]为多边形比例边界有限单元的位移模态;[N( η )]为边界上线性单元节点的插值函数;
                为矩阵[Z]的特征值;{c}为积分常数向量。
              λ b
                  由位移导出对应多边形内任一点的应变和应力,其表达式为:
                                              ( ε )   - 1      λ 1   λ 2        λ n - 2
                                                                             n - 2
                                                             1
                                                                   2
                               { ε ( ξ ,η )} = [ Φ b  ( η )] ξ diag(c ξ ,c ξ ,…,c ξ ,0,0)               (3)
                                                ( ε )   - 1      λ 1   λ 2         λ n - 2
                                                               1
                            { σ ( ξ ,η )} = [D][ Φ b  ( η )] ξ diag(c ξ ,c ξ ,…,c ξ ,0,0)               (4)
                                                                               n - 2
                                                                     2
                       ( ε )
              式中:[ Φ b   ( η )]为多 边形比 例边 界有 限单 元的应 变 模态,具体 推导 及表 达 式 详见 文 献 [18 - 19];
              [ D]为材料弹性矩阵。
              2.2 刚度矩阵和质量矩阵 对于任意的多边形比例边界有限单元,单元刚度矩阵[K]表达式为:
                                                                                           e
                                                           (q)    (u)  - 1
                                                     e
                                                  [K] = [ Φ b  ][ Φ b  ]                                (5)
                     (q)
              式中[ Φ b   ]为多边形比例边界有限单元的节点力特征向量集。
                  多边形比例边界有限单元质量矩阵由边界上线性单元根据单元节点连接组装得到,其表达式为:
                                                        (u)  - 1 T     (u)  - 1
                                              [ M ] = ([ Φ b  ] )[m][ Φ b  ]                            (6)
                                                 e
              式中[ m]为系数矩阵的缩略表达式,其仅取决于多边形单元的密度和单元形状                                   [15] 。
              2.3 多边形比例边界有限元静动力方程 对于由多边形单元离散的任意结构体系,结构整体刚度阵
              [ K]和整体质量阵[M]类似于有限元法通过单元组装形成,因此静力平衡方程和动力方程分别为:
                                                      [ K]{d} = {F}                                     (7)
                                                   ··       ·
                                             [M]{d} + [C]{d} + [K]{d} = {F}                             (8)
                                                 ·    ··
              式中:[ C]为结构阻尼矩阵;{d}、{d}、{d}分别为节点位移、速度和加速度;{F}为荷载向量。
                  实际工程动力分析需要考虑阻尼的影响,采用瑞利阻尼,阻尼矩阵通过结构整体刚度阵和质量阵
              的线性组合得到,其表达式为:
                                                                  [K]                                   (9)
                                                   [ C] = α 0 [M] + α 1
                         为组合系数,均与结构的频率和材料阻尼相关。
              式中 α 0 、α 1
                  采用 Newmark方法求解式(8)可获得任意时刻结构的位移场、速度场和加速度场,应变场和应力
              场可分别由式( 3)(4)通过位移场导出。
                  从总的流程来说,对于采用多边形比例边界有限单元离散的结构动态断裂模拟,类似于有限元
              法,每一时间步先进行结构动力分析,在此基础上求解结构断裂参数并进行动态断裂分析。

              3 广义动态应力强度因子及裂缝扩展准则


              3.1 广义裂缝动态应力强度因子 对于无限平面内以速度为 v扩展传播的稳态裂缝,其动态应力强度

              因子  [24] 计算公式为:
                                                                        eq
                                           dyna
                                          K ( θ ,t)                   K ( θ ,t)
                                           dyna {  Ⅰ  } =[k(v) k(v)]    eq {  Ⅰ  }                     (10)
                                          K ( θ ,t)      Ⅰ      Ⅱ     K ( θ ,t)
                                                                        Ⅱ
                                           Ⅱ
              式中:应力强度因子 K上标 “eq” 表示依据瞬时应力场计算的瞬态应力强度因子,上标 “dyna” 代表
              动态应力强度因子;k(v)、k(v)与裂缝传播速度和材料特性相关,其表达式为:
                                  Ⅰ       Ⅱ
                                                                                                —  1 1 7 —
   117   118   119   120   121   122   123   124   125   126   127