066 岩土力学 0 年对比两种方法计算过程中单元刚度矩阵变化规律位移变化情况验证过程中发现, 现有数值流形法在计算过程中改变惯性矩阵重力场矩阵时, 忽略了质量守恒

Size: px

Start display at page:

Download "066 岩土力学 0 年对比两种方法计算过程中单元刚度矩阵变化规律位移变化情况验证过程中发现, 现有数值流形法在计算过程中改变惯性矩阵重力场矩阵时, 忽略了质量守恒"

学匹朱
7 years ago
Views:

1 第卷第 0 期岩土力学 Vol. No.0 0 年 0 月 Rock and Soil Mechanics Oct. 0 文章编号 : (0) 数值流形法中质量守恒的探讨林兴超, 汪小刚, 王玉杰, 李旭, 韩鑫晔 (. 中国水利水电科学研究院岩土工程研究所, 北京 00048;. 北京交通大学土木建筑工程学院, 北京 00044;. 北京航空航天大学交通科学与工程学院, 北京 009) 摘要 : 数值流形法是至少包含流形法 (NMM) 有限元法 (FEM) 和非连续变形分析 (DDA) 的数值方法体系将数值流形法中物理单元与数学单元完全重合, 去掉接触理论, 流形元能够回归到有限元, 将通过简单的板压缩数值试验验证这一点在以前的数值流形法法中, 质量守恒问题一直被忽视, 物理单元的质量会随着单元体积改变, 计算结果存在一定的误差通过改变计算过程中单元密度实现计算过程中的质量守恒, 完善了现有数值流形法的理论基础关键词 : 数值流形法 ; 有限元 ; 非连续变形分析 ; 刚度矩阵 ; 质量守恒中图分类号 :O 4 文献标识码 :A Discussion on mass-conservation in numerical manifold method LIN Xing-chao, WANG Xiao-gang, WANG Yu-jie, LI Xu, HAN Xin-e (. Department of Geotechnical Engineering, China Institute of Water Resources and Hdropower Research, Beijing 00048, China;. School of Civil Engineering &Architecture, Beijing Jiaotong Universit, Beijing 00044, China;. School of Transportation Science and Engineering, Beihang Universit, Beijing 009, China) Abstract: Numerical manifold method (NMM) is a numerical sstem which contains manifold method, finite elements and discontinuous deformation analsis (DDA). NMM can return to finite elements when phsical elements and mathematical elements in NMM are superposition completel; and contact theor is removed. It will be verified b plate compression numerical eperiment in the paper. In previous NNM, the problem of mass-conservation is neglected, which the weight of the phsical elements will change with changes in its volume; so there are certain error in the calculation results. Mass-conservation can be realized b changing the elements densit in the process of calculation; that will perfect the theoretical basis of numerical manifold methods. Ke words: numerical manifold method(nmm); finite elements(fe); discontinuous deformation analsic(dda); stiffness matri; mass- conservation 引言数值流形法 (numerical manifold method,nmm) 是石根华博士在非连续变形分析 ( discontinued deformation analsis,dda) 基础上, 融入了有限元和解析法, 于 99 年提出的一种至少包含有限元和 DDA 的数值方法体系 [-] 在数值流形法诞生的近 0 年来, 国内外学者对其进行了大量研究 : 王水林 [4] 张国新 [5] [6] MA 等通过数值流形法对裂纹扩展和节理岩体破坏过程进行模拟 ;KENJIRO 等 [7],JIANG 等 [8] [9] CHENG 等等对数值流形法在三维空间中的应用进行研究 ; 韩 [0] 有民等探讨了基于三角形有限元网格的裂纹扩 [] 展时物理覆盖与流形单元的生成算法 TIAN 等 [-] 对 4 种不同的四面体单元对比分析 ; 朱爱军研究数值流形法的固定边界约束处理方法进行, 通过改变广义节点上的覆盖函数达到边界固定的效果, 并对岩土工程中开挖卸荷过程进行模拟从已有文献看, 对数值流形法发展和应用的研究较多, 而在数值流形法与有限元对比验证的研究较少, 关于质量守恒问题未见文献提及为了便于对比分析, 将数值流形法中物理单元与数学单元完全重合, 去掉接触理论, 建立与有限元具有相同计算条件的数值流形法计算模型通过收稿日期 : 基金项目 : 国家自然科学基金重点项目资助 (No ); 十一五国家科技支撑项目资助 (No. 008BAB9B0,No. 008BAB9B0); 中国水利水电科学研究院博士生学位论文创新研究资助 (009 年 ) 第一作者简介 : 林兴超, 男,984 年生, 博士研究生, 主要从事岩土力学实验及数值计算方面工作 linc@iwhr.com

北京航空航天大学交通科学与工程学院, 北京 009) 摘要 : 数值流形法是至少包含流形法 (NMM) 有限元法 (FEM) 和非连续变形分析 (DDA) 的数值方法体系将数值流形法中物理单元与数学单元完全重合, 去掉接触理论, 流形元能够回归

2 066 岩土力学 0 年对比两种方法计算过程中单元刚度矩阵变化规律位移变化情况验证过程中发现, 现有数值流形法在计算过程中改变惯性矩阵重力场矩阵时, 忽略了质量守恒流形元采用单纯形积分方法, 在整个物理单元上进行质量和重力精确积分, 进而分步计算中改变惯性矩阵重力场矩阵现有流形元方法, 未考虑单元密度随着单元体积的改变, 因此, 单元的质量和重力会随着单元的体积改变而发生改变计算过程中随着材料的压缩或膨胀, 物体质量减小或增大这违背了质量守恒, 在计算大变形问题时, 会导致误差增大为实现计算过程中的质量守恒, 在流形元的每步计算中都需要对单元密度进行修正这一问题的讨论未见于其他文献本文的数值模拟结果表明, 在流形元的计算中改变单元密度以保证质量守恒是重要的, 并且是有效的在计算中考虑质量守恒将使流形元理论基础更为严密单元刚度矩阵本文对数值流形法中单元刚度矩阵变化规律进行了研究, 这里简要介绍流形元中的三角形单元刚度矩阵流形元刚度矩阵的积分域是物理单元, 它可以充满整个数学单元, 也可以只是数学单元的一部分当流形元的数学单元的广义节点和有限元的单元节点重合时, 流形元的数学单元和物理单元重合, 且与有限元中的有限单元重合此时, 数值流形法中的积分域和计算域一致, 则数值流形法与有限单元法就具有相同的理论框架考虑广义虎克定律下的弹性平面应力问题, 数值流形法与有限元一样, 应力 - 应变关系由式 () 给出 0 E 0 E z 0 0 () 0 E 式中 :E 0, E 分别为 0 0 弹性模量和泊松比节点位移函数表示为 u(, ) (, ) r r v D e() Te( ) De( ) ( Te() Te() Te() ) De() r D e() () 并有 u(, ) v(, ) v(, ) u(, ) () 其中 : we( i) (, ) 0 Te( i) (, ) 0 we( i) (, ) (4) ue( i) De( i) ve( i) (5) 及 we() (, ) f f f we() (, ) f f f we() (, ) f f f (6) ue( i) De( i) ve( i) (7) ue() v e() f 0 f 0 f 0 ue() 0 f 0 f 0 f (8) ve() f f f f f f u e() ve() 令 fi 0 B e( i) 0 fi, ( i,, ) fi f i (9) 则 D e() [ B B B ] D B D D e() e() e() e() e( i) e( i) e() 由单元 e 的弹性应力所产生的应变能 e d d A T T De Be EBe Dedd A T T De Be EBe ddde A T e T De S Be EBe De e (0) ()

这违背了质量守恒, 在计算大变形问题时, 会导致误差增大为实现计算过程中的质量守恒, 在流形元的每步计算中都需要对单元密度进行修正这一问题的讨论未见于其他文献本文的数值模拟结果表明, 在流形元的计算中改变单元密度

3 第 0 期林兴超等 : 数值流形法中质量守恒的探讨 067 式中 : S e 为单位的面积, 则单元刚度矩阵为 T e S Be E B e 算例分析为了对比相同计算条件下数值流形法与有限元的计算精度, 将数值流形法中物理单元与数学单元完全重合, 去掉接触理论, 建立与有限元具有相同计算条件的数值流形法计算模型, 对数值流形法计算结果与有限元程序 Abaqus 计算结果及理论解进行对比分析. 算例及计算模型算例采用平板重力压缩试验 ( 平面应变问题 ), 模型示意图见图, 计算边界及材料参数见表, 重力加速度取 0 m/s 流形元程序及 abaqus 程序计算模型如图所示由图可见, 流形元程序数学单元与物理单元完全重合, 数学单元网格与有限元程序计算模型网格划分相同, 计算模型由 55 个节点和 80 个单元组成 (a) 流形元 (b) abaqus 图计算网格 I Fig. Computational grid I 密度 / (kg/m ) 图平板压缩试验示意图 Fig. Sketch of flat conpress test 表算例计算边界及材料参数表 Table Calculation boundaries and materials 高 / mm 计算边界宽 / mm 弹性模量 E / MPa 变形参数泊松比计算条件 : 计算过程中只考虑重力作用, 流形元程序计算时间步长设为 0., 时步设为 400, 有限元分为 0 个增量步. 计算结果分析 () 数值流形法单元刚度矩阵分析由单元刚度矩阵计算公式可知, 三角形单元刚度矩阵为 6 6 对称矩阵, 因此, 结果中只统计刚度矩阵的下三角形矩阵, 计算结果如表所列在计算过程中, 数值流形法单元刚度矩阵根据变化后的数学单元广义节点坐标计算单元刚度矩阵, 单元刚度矩阵不断更新 () 监测点位移计算分析为了解不同条件下数值流形法计算结果有限元法计算结果和理论解间的对比情况, 布置了个监测点, 具体位置如图所示, 图中 (a, b) 表示监测点坐标表有限元数值流形法单元刚度矩阵对比 (80 号单元 ) Table Stiffness of NMM and FEM (Element 80) 比较流形元计算结果 ( 时步 ) 有限元计算结果 ( 增量步 ) 时步增量步时步增量步

算例及计算模型算例采用平板重力压缩试验 ( 平面应变问题 ), 模型示意图见图, 计算边界及材料参数见表, 重力加速度取 0 m/s 流形元程序及 abaqus 程序计算模型如图所示由图可见, 流形元程序数学单元与物理单元完全重合, 数学

4 068 岩土力学 0 年监测点 (, 9) 监测点 (, 8) 增大或减小本文将在第 4 节中重点讨论计算过程中质量守恒问题监测点 (, 5) 表 4 不同弹性模量条件下监测点计算结果比对 Table 4 Comparison between numerical solutions and theoretical solution (different elasticit modulus) 测点方向计算方法不同弹性模量 (MPa) 下位移 / mm 图监测点示意图 ( 单位 :mm) Fig. Location of monitoring points(unit: mm) 0, 板自重工况下理论值数值流形法计算和有限元计算结果如表所列有限元数值流形法理论值有限元数值流形法在自重条件下, 竖向位移 ( 方向 ) 有限元结果与理论解一致, 数值流形法计算结果与理论解有 0%~0.06% 的误差, 误差非常微小, 可以认为该工况下数值流形法计算程序能够得到精确解 ; 水平向有限元数值流形法理论值位移误差较大, 这是数值计算中由于浮点数的影响, 有限元在 0 附近产生较大误差, 有限元和数值流形法计数值流形法算结果基本一致, 两者间的差值不超过 0.0 m 以上计算结果表明, 小变形情况下数值流形法有限有限元元均能得到与理论解值一致的结果数值流形法理论值表监测点计算结果比对 Table Comparison between numerical solutions and theoretical solution 有限元数值流形法监测点坐标量值不同计算方法的位移 / m / mm 理论解数值流形法有限元不同弹性模量条件下理论值数值流形法计算结果有限元计算结果如表 4 所列在弹性模量为 MPa 条件下流形元计算结果与有限元计算结果的差值随着弹性模量的增大而减小监测点的方向位移 ( 为了便于对比, 将位移分别乘以 ) 随与弹性模量变化的示意图见图 4 由表 4 中数据和图 4 可以看出, 在弹性模量较大时, 数值流形法计算结果与理论解一致 ; 随着弹性模量的减小, 数值流形法计算结果逐渐偏离理论解结合刚度矩阵结果发现, 现有数值流形法在计算过程中忽略了质量守恒, 导致计算过程中质量图 4 监测点的方向位移随弹模变化曲线 Fig.4 displacement along curves of monitor point 4 数值流形法中质量守恒的探讨数值流形法通过动态方程求解静态问题, 计算中使用时步积分方案, 分步计算, 最终系统达到稳定状态在分步计算过程中, 刚度矩阵不断发生变化, 当产生大变形时, 计算结果与理论值相差较大 ( 表 4) 误差来源是求解静态方程过程中由于板压缩, 通过板重度在单元面积上积分得到的板重力减小所致

Location of monitoring points(unit: mm) 0, 板自重工况下理论值数值流形法计算和有限元计算结果如表所列有限元 0.09 0.09 0.09 0.09 0.09 数值流形法 0.568 0.05 0.0 0.09 0.09 理论值 4.950 4.950 4.950 4.950 4.950 有限元 4.

5 第 0 期林兴超等 : 数值流形法中质量守恒的探讨 069 数值流形法中, 对于二维问题, 单元质量可通过式 () 计算得到 G Ag () 式中 : 为单元密度, 为常量 ; A 为单元当前计算步下面积, 随节点位移变化而变化 ;g 为重力加速度, 也为常量因此, 通过式 () 计算得到的单元质量 G 为变量, 这与事实是相违背的, 没有遵守质量守恒定律在新的计算模拟中, 每一时步都对单元密度按照式 () 进行重新求解, 使计算过程中遵循质量守恒 M A () A A 将式 () 中代入式 () 中得到 G Ag C( 常量 ) (4) 在每一时步中修正单元密度后, 计算得到的结果如表 5 所列, 其结果表明, 考虑质量守恒后, 流形元计算得到结果的精度大大提高, 方向位移列出前 4 位均与理论解相同, 说明数值计算结果与理论解误差小于 0.% 密度无关的体积力, 如浮力, 现有流形元的计算方法是精确的, 即在整个单元体积上进行精确积分求解体积力数学单元物理单元两套网格的使用是数值流形法的突出优点之一, 它极大简化了大变形或裂纹扩展条件下模型变化过程中网格重构问题 [5] 本文算例为了方便对比分析将数学单元与物理单元完全重合, 削弱了数学单元物理单元两套网格使用的优势因此, 重新建立计算模型, 使数学单元不受物理网格的影响, 计算网格如图 5 所示, 关键测点的方向位移如表 6 所列从表中数据可看出, 数学单元与物理单元是否一致对计算结果影响不大监测点 (, 9) 监测点 (, 8) 监测点 (, 5) 测点表 5 修正单元密度后监测点计算结果比对 Table 5 Comparison between numerical solutions and theoretical solution (revise cell densit) 方向计算方法不同弹性模量 (MPa) 下位移 / mm 有限元数值流形法理论值有限元数值流形法有限元数值流形法理论值有限元数值流形法有限元数值流形法理论值有限元数值流形法需要指出的是 : 对于动力问题, 惯性矩阵起着非常重要的作用在流形元的计算中, 也必须采用修正的密度或者使用固定的质量对于和单元图 5 计算网格 II ( 单位 :mm) Fig.5 Computational grid II (unit: mm) 表 6 物理单元与数学单元是否一致的情况下监测点计算结果比对 Table 6 Comparison between grid I and grid II solution 测点方向计算网格不同弹性模量 (MPa) 下位移 / mm I II 为进一步验证本文对质量守恒考虑的必要性和可靠性, 对多层土基自重压密过程进行模拟, 计算参数及边界条件如图 6(a) 所示, 计算结果见图 6(b) 计算结果表明, 本算例中考虑质量守恒后的数值流形法计算结果与 abaqus 计算结果一致, 未考虑质量守恒的数值流形法计算误差在 0% 左右因此, 数值流形法在计算过程中对单元密度进行修正, 使其遵循质量守恒是必要的, 本文提出对单元密度的修正方法是可靠的图 6 中, 图中 (a, b) 表示模型控制点坐标 ; 土层参数中, E 为弹性模量 ; 为泊松比 ;d 为密度 ; 模型左右两侧为法向约束, 底部固定约束, 上部为自由边界

G Ag C( 常量 ) (4) 在每一时步中修正单元密度后, 计算得到的结果如表 5 所列, 其结果表明, 考虑质量守恒后, 流形元计算得到结果的精度大大提高, 方向位移列出前 4 位均与理论解相同, 说明数值计算结果与理论解误差小于 0.

6 070 岩土力学 0 年译. 北京 : 清华大学出版社, 997. [] SHI G H. Manifold method of material analsis[c]// Transactions of the 9th Arm Conference on Applied Mathematics and Computing. Minneapolis, Minnesota: [s. n.], 99: [] SHI G H. Modeling rock joints and blocks b manifold method[c]//proceedings of the rd US Rock Mechanics (a) 计算参数及边界条件 ( 单位 : m) (b) 位移计算结果图 6 多层土基模型及计算结果 Fig.6 Multilaer soil foundation model and results 5 结论 () 数值流形法是涵盖了有限元法 (FEM) 非连续变形分析法 (DDA) 的新的数值方法体系, 能够统一处理连续和非连续问题当流形元的数学单元广义节点和物理单元节点完全一致时, 数值流形法能够回归到有限元 () 用流形元分时步进行计算, 在每一时步中, 流形元根据变化后的节点坐标更新单元刚度矩阵, 单元刚度矩阵不断变化流形元在计算过程中不断更新单元坐标因此, 通过多个时间步长内的小变形累加, 可以实现大变形问题的模拟 () 在流形元的分步计算中, 可以通过修正单元密度以实现质量守恒计算结果表明, 考虑质量守恒的流形元模拟提高了计算的精度 (4) 通过算例分析, 物理单元与数学单元是否一致对计算结果影响不大数学单元物理单元两套网格的使用大大简化了网格剖分, 特别是大变形和裂纹扩展条件下网格重构问题, 体现了数值流形法在数值计算中的优势致谢 : 本文采用的数值流形法计算程序由石根华博士提供, 并得到了他的支持与指导, 在此深表感谢参考文献 [] 石根华. 数值流形方法与非连续变形分析 [M]. 裴觉民 Smposium. New Meico: [s. n.], 99: [4] 王水林, 冯夏庭, 葛修润. 高阶流形方法模拟裂纹扩展研究 [J]. 岩土力学, 00, 4(4): WANG shui-lin, FENG Xia-ting, GE Xiu-run. Stud on crack propagation modeling b high order manifold method[j]. Rock and Soil Mechanics, 00, 4(4): [5] 张国新, 赵研, 彭校初. 考虑岩桥断裂的岩质边坡倾倒破坏的流形元模拟 [J]. 岩石力学与工程学报, 007, 6(9): ZHANG Guo-in, ZHAO Yan, PENG Jiao-chu. Simulation of toppling failure of rock slope b numerical manifold considering fracture of rock bridges[j]. Chinese Journal of Rock Mechanics and Engineering, 007, 6(9): [6] MA G W, AN X M, ZHANG H H, et al. Modeling comple crack problems using the numerical manifold method[j]. International Journal of Fracture, 009, 56: [7] KENJIRO Terada, MAO Kurumatani. An integrated procedure for three-dimensional structural analsis with the finite cover method[j]. International Journal for Numerical Methods in Engineering, 005, 6: 0-. [8] JIANG Qing-hui, ZHOU Chuang-bing, LI Dian-qing. A three-dimensional numerical manifold method based on tetrahedral meshes[j]. Computers and Structures, 009, 87: [9] CHENG Y M, ZHANG Y H. Formulation of a threedimensional numerical manifold method with tetrahedron and heahedron elements[j]. Rock Mec- hanics and Rock Engineering, 008, 4(4): [0] 韩有民, 罗先启, 王水林, 等. 裂纹扩展时物理覆盖与流形单元的生成算法 [J]. 岩土工程学报, 005, 7(6): 下转第 074 页

Modeling rock joints and blocks b manifold method[c]//proceedings of the rd US Rock Mechanics (a) 计算参数及边界条件 ( 单位 : m) (b) 位移计算结果图 6 多层土基模型及计算结果 Fig.

7 第 0 期林兴超等 : 数值流形法中质量守恒的探讨 HAN You-min, LUO Xian-qi, WANG Shui-lin, et al. Formation algorithm of covers and manifold elements in NMM during propagation of cracks[j]. Chinese Journal of Geotechnical Engineering, 005, 7(6): [] TIAN Rong, MATSUBARA Hitoshi, YAGAWA Genki. Advanced 4-node tetrahedrons[j]. International Journal for Numerical Methods in Engineering, 006, 68: [] 朱爱军, 邓安福, 曾祥勇, 等. 岩体工程数值流形方法的固定边界约束处理方法 [J]. 岩石力学与工程学报, 005, 4(9): ZHU Ai-jun, DENG An-fu, ZENG Xiang-ong. Constraints treatment of fied boundar in numerical manifold method[j]. Chinese Journal of Rock Mechanics and Engineering, 005, 4(9): [] 朱爱军, 邓安福, 曾祥勇. 数值流形法对岩体工程开挖卸荷问题的模拟 [J]. 岩土力学, 006, 7(): ZHU Ai-jun, DENG An-fu, ZENG Xiang-ong. Numerical manifold method for simulation of ecavation unloading in geotechnical engineering[j]. Rock and Soil Mechanics, 006, 7(): 79-8.

International Journal for Numerical Methods in Engineering, 006, 68: 09 -. [] 朱爱军, 邓安福, 曾祥勇, 等. 岩体工程数值流形方法的固定边界约束处理方法 [J]. 岩石力学与工程学报, 005, 4(9): 58-587.

龚亚夫在重新思考基础教育英语教学的理念一文中援引的观点认为当跳出本族语主义的思维定式后需要重新思考许多相连带的问题比如许多发音的细微区别并不影响理解和

龚亚夫在重新思考基础教育英语教学的理念一文中援引的观点认为当跳出本族语主义的思维定式后需要重新思考许多相连带的问题比如许多发音的细微区别并不影响理解和语音语篇语感语域林大津毛浩然改革开放以来的英语热引发了大中小学英语教育整体规划问题在充分考虑地区学校和个体差异以及各家观点的基础上遵循实事求是逐级定位逐层分流因材施教的原则本研究所倡导的语音语篇语感语域

066 岩 土 力 学 0 年 对 比 两 种 方 法 计 算 过 程 中 单 元 刚 度 矩 阵 变 化 规 律 位 移 变 化 情 况 验 证 过 程 中 发 现, 现 有 数 值 流 形 法 在 计 算 过 程 中 改 变 惯 性 矩 阵 重 力 场 矩 阵 时, 忽 略 了 质 量 守 恒

066 岩土力学 0 年对比两种方法计算过程中单元刚度矩阵变化规律位移变化情况验证过程中发现, 现有数值流形法在计算过程中改变惯性矩阵重力场矩阵时, 忽略了质量守恒