当前位置： >>
岩石力学三维有限元分析的代数多重网格求解法
第 25 卷 第 11 期 2006 年 11 月岩石力学与工程学报 Chinese Journal of Rock Mechanics and EngineeringVol.25 No.11 Nov.，2006岩石力学三维有限元分析的代数多重网格求解法谢学斌 1，肖映雄 2，舒 适 3，潘长良 1(1. 中南大学 资源与安全工程学院，湖南 长沙 . 湘潭大学 基础力学与材料工程研究所，湖南 湘潭 411105； 3. 湘潭大学 计算与应用数学研究所，湖南 湘潭 411105)摘要：多重网格法是一种求解由偏微分方程边值问题所导出的代数方程组的快速算法，几何多重网格法存在某些 缺陷，影响它的推广应用。采用代数多重网格法求解岩石力学三维有限元离散线性方程组，简要介绍代数多重网 格三维粗网格形成方法与三维插值算子，利用研制的基于代数多重网格法的三维有限元程序进行一系列数值试验。 结果表明：代数多重网格法求解各种复杂计算条件下岩石力学三维有限元方程时具有良好的收敛特性和较强的适 应能力，计算效率远高于直接法求解器，为大规模岩土工程三维有限元分析提供一种快速有效的方法。 关键词：岩石力学；代数多重网格法；粗化技术；插值算子；三维有限元；数值方法 中图分类号：TU 45 文献标识码：A 文章编号：06)11C2358C06ALGEBRAIC MULTIGRID METHOD FOR THREE-DIMENSIONAL FINITE ELEMENT ANALYSIS OF ROCK MECHANICSXIE Xuebin1，XIAO Yingxiong2，SHU Shi3，PAN Changliang1(1. School of Resources and Safety Engineering，Central South University，Changsha，Hunan 410083，China； 2. Institute of Fundamental Mechanics and Material Engineering，Xiangtan University，Xiangtan，Hunan 411105，China； 3. Institute of Computational and Applied Mathematics，Xiangtan University，Xiangtan，Hunan 411105，China)Abstract：Multigrid method solver is of high numerical efficiency when used in solving linear equations derived from boundary-value problems of partial differential equations. There are some shortages in geometrical multigrid method which restricts its application area. The algebraic multigrid method is used to solve finite element linear equations which are derived from three-dimensional finite element analysis of rock mechanics and engineering. The three-dimensional coarse-grid selection method based on element agglomeration and the three-dimensional interpolation operator are briefly introduced. By using the newly developed three-dimensional finite element program based on algebraic multigrid method，four different numerical experiments are designed，and carried out to validate its convergence character，numerical efficiency and practical application to modeling excavation problem of rock engineering. The numerical experiments show that the algebraic multigrid method is of better stability，good convergence character and better adaptability，with much higher numerical efficiency with increasing of the number of linear equations and much less computer memory compared with direct method. Increasing. The algebraic multigrid method has much better numerical efficiency. The algebraic multigrid method is suitable and efficient for three-dimensional finite element modelling of large-scale geomechanical engineering. Key words：rock mechanics；algebraic multigrid method；coarsening technique；interpolation operator； three-dimensional finite element；numerical method收稿日期：；修回日期： 基金项目：国家自然科学基金资助项目()；高性能科学计算研究资助项目() 作者简介：谢学斌(1968C)，男，博士，1989 年毕业于中南工业大学采矿工程专业，现任副教授，主要从事岩土工程数值模拟研究和岩土灾害防治方 面的教学与研究工作。E-mail：xbxie@csu.edu.cn第 25 卷第 11 期谢学斌等. 岩石力学三维有限元分析的代数多重网格求解法T? 2359 ?1引言U 2 ，…， U N } ， U i ( i = 1，2，…，N)为 3×1 的块列阵； F 为综合考虑体力、边界力后的荷载列阵， 且 F = {F1 ， F2 ，…， FN }T ， Fi 为 3×1 的块列阵。 利用代数多重网格法求解上述由三维有限元分析导 出的线性方程组包括网格预处理和多重网格迭代 2 个过程，简述如下。 2.1 网格预处理方法 考虑 S 个网格层的情形，S 为最细网格层的层 数，nS := N 为最细网格层的节点数，3nS 阶矩阵 K S 为方程组式(1)的系数矩阵。 nk 为第 k 个网格层的 设 节点数，相应的离散解空间为 Vk ， K k 为该网格层 下的粗化线性代数方程组的系数矩阵， k = S ，研究适用于岩石力学问题有限元离散方程组 求解， 数值效率高， 需要计算机存储量少的新算法， 对于提高有限元法模拟三维超大规模岩土工程问题 的能力，具有较大的意义。多重网格法具有需要计 算机存储量少，计算速度快的优点，已有的理论研 究和数值试验[1～5]证明：多重网格法比直接法和其他常用单网格迭代法具有更高的数值效率，是现有 求解偏微分方程数值解的最有效的快速算法之一。 基于几何网格和物理意义的几何多重网格法在岩体 力学有限元分析中的应用存在较大的局限性，如当 求解区域为三维复杂边界条件或者网格非均匀剖分 形成复杂网格结构，以及材料为多介质或不连续情 形时，应用几何多重网格法都会遇到困难，因此目 前的研究仅局限于平面问题的有限元分析[5]。 代数多重网格法仅仅利用方程组的系数矩阵， 从刚度矩阵的代数结构出发，自动形成虚拟的粗细 网格，且不依赖于原问题的几何与物理性质，在计 算原理和方法上较几何多重网格法具有明显的优越 性[6]。作者曾将代数多重网格法用于岩体力学平面 问题有限元分析[7]，在此基础上，本文研究采用基 于聚集法的代数多重网格法求解岩石力学三维有限 元离散方程组，并采用一维稀疏存储技术，研制了 基于代数多重网格求解的岩石力学三维有限元软 件。数值试验结果表明，代数多重网格法求解三维 岩石力学有限元离散方程组是稳定的，具有良好的 收敛性能， 计算 CPU 时间远小于直接法， 具有较高 的数值效率。与平面问题有限元分析数值效率 对 比，本文发展的三维问题代数多重网格求解器，具 有更大的实用意义和更高的数值效率，为在微型机 上进行大规模岩土工程问题三维数值模拟提供了一 种较好的方法。[7]S ? 1 ，…，1。利用矩阵 K S 定义如下 nS 阶网格矩阵 ~ ~ K S = ( K ij ) s ，即~ ?0 ? K ij = ? ?1 ? ( K ij = 0) ( K ij ≠ 0) (2)~ 矩阵 K S 对应一个图， 记为 ?S = (VS ，E S ) ， S ， VE S 分别为顶点集合和边集合，称该图为最细层(第 S 层)网格剖分图。假设第 s 层粗网格 ?s ( s = S ， S ? 1 ，…，2)及相应的系数矩阵 K s = ( K ij ) s 已知，则更粗一层网格 ?s?1 及其系数矩阵 K s ?1 = ( K ij ) s ?1 可 由以下步骤[7]完成：(1) 利用下述各步，生成从第 s ? 1 层到第 s 层网格间的网格转移算子 Pss?1 ： ① 利用聚集法[8]，求第 s 层网格剖分图 ?s =(Vs ， E s ) 的极大不相关集，将它作为第 s ? 1 层节点集合，并记 ns ?1 为节点数。 ② 定义插值关系[7]， 由此确定从第 s ? 1 层网格 ~ 到第 s 层网格间的网格转移算子 Pss?1 。 ~ ③ 将 ns × ns ?1 阶转移算子 Pss?1 进行如下扩充： 对 i = 1 ，2，…， ns 及 j = 1，2，…， ns ?1 ，令 ~ ( Pss?1 ) 3i ? 2，3 j ? 2 = ( Pss?1 ) 3i ?1，3 j ?1 = ( Pss?1 ) 3i，3 j = ( Pss?1 ) i，j( Pss?1 ) 3i ?2，3 j ?1 = ( Pss?1 ) 3i ?2，3 j = ( Pss?1 ) 3i ?1，3 j ?2 =2代数多重网格法设在岩石力学三维有限元分析过程中，将研究( Pss?1 ) 3i ?1，3 j = ( Pss?1 ) 3i， j ? 2 = ( Pss?1 ) 3i， j ?1 = 0 3 3得到从第 s ? 1 层到第 s 层网格间阶为 (3ns ) × (3ns ?1 ) 的整体转移算子，记为 Pss?1 。对象离散化后，得到的有限元线性方程组为(2) 对上述得到的转移算子利用细网格系数矩阵 K s 进行磨光处理[9]， 得到能量极小意义下的转移 算子，仍记其为 Pss?1 。KU = F(1)式中： K 为总刚度矩阵，且 K = (K ij ) N × N ， K ij (i， j = 1，2，…，N)为 3×3 的块矩阵，N 为三维有限元 且 网格的节点数；U 为系统节点位移列阵， U = {U 1 ，(3) 作代数运算 K s ?1 = ( Pss?1 ) T K s ( Pss?1 ) ，得到第s ? 1 层的系数矩阵。重复上述各步骤，当粗点个数足够少时，网格? 2360 ?岩石力学与工程学报2006 年处理工作即告结束，并记所得的网格总层数为 S。 在网格预处理算法中，若 S = 1 ，则可用直接法或迭 则转入下述多重网格迭代。 代法直接求解； S ＞1， 若 2.2 代数多重网格迭代算法 这里仅介绍两重网格方法。一次引用 AMG(S，与平面应力条件下厚壁圆筒受均匀外压作用应力分 布完全相同，其应力可用解析法精确求解。 对上述模型首先进行有限元分析，计算材料参 数为：弹性模量 E = 20 000 MPa，泊松比 μ = 0.3，σ σ 计算结果中所有各单元的 σ Z 都约等于 0； r ， θ 沿径向分布情况分别见表 1 和 2。然后采用解析法求 得相应的解析解，结果详见表 1 和 2。表 1 径向应力σr 数值解和解析解比较表 Table 1 Comparison of radial stress σr between numerical and analytical results径向距离 /m 1.5 2.5 3.5 4.5 8.5 数值解 /MPa 5.592 5 8.514 8 9.261 8 9.596 8 9.920 1 解析解 /MPa 5.594 4 8.458 7 9.248 2 9.572 4 9.930 6 相对误差 /% 0.030 0.660 0.150 0.027 0.106p1 ， p2 ， p ) 迭代过程如下：(1) 选取初始值 U 0 ，进行 p1 次点块 GaussSeidel 迭代，结果记为 U1 。 (2) 求残量 r1 = F ? KU1 ，将残量 r1 进行限制得r2 = P12 r1 ，对 ?2 上的误差方程 K 2 e2 = r2 ，p 次引用 AMG ( S ? 1 ， p1 ， p2 ， p ) 过程，若 p = 1 为 VC循环方法， p = 2 为 WC循环方法，记相应结果为 e2 。(3) 粗网格校正，计算 U 2 = U 2 + P21e2 。 (4) 以校正值 U 2 为初始值，再进行 p2 次点块 Gauss-Seidel 迭代。在本文的有限元程序中，取 AMG(S，3，3，1)， 对应着代数多重网格 VC循环算法。3数值试验表2切向应力σθ数值解与解析解比较表Table 2 Comparison of tangential stress σθ between numerical3.1 厚壁圆筒受均匀外压作用应力分析 一均质各向同性线弹性材料的厚壁圆筒，圆筒 内径为 a = 1 m，外径 b = 12 m，在圆筒的外壁受径 向均布外压作用，压力大小为 q = 10 MPa。取 1/4 对称结构，采用如图 1 所示的有限元网格计算厚壁 圆筒受均布外压作用下的应力分布。模型中坐标 Z 方向与圆筒的轴线方向相同，X-Y 平面平行于厚壁 圆筒的横截面，坐标原点与圆心重合，共划分 25 个八节点等参单元。在圆筒的轴线方向只划分一层 单元，圆筒高度即为三维等参单元的厚度，将圆筒 高度取为 1 m，所有节点在 Z 方向的约束条件设置 为自由变形。上述特殊计算条件下的空间应力分布1.5 2.5 3.5 4.5 8.5and analytical results径向距离/m 数值解/MPa 14.520 5 11.663 8 10.816 8 10.493 9 10.216 5 解析解/MPa 14.545 5 11.681 1 10.891 6 10.567 2 10.208 8 相对误差/% 0.172 0.149 0.690 0.700 0.070分析表 1 和 2 的结果可知，径向应力 σ r 数值解 与解析解的相对误差为 0.027%～0.660%，切向应力σ θ 数值解与解析解相对误差为 0.070%～0.700%，说明本文的代数多重网格迭代求解方法在如算例 1 中的求解几何区域为复杂三维曲面且为非均匀网格 剖分的复杂条件下收敛于解析解，具有良好的数值 收敛特性。 3.2 节理化岩体受均布压力作用无侧限条件下应 力和位移分析 计算模型为一长方体岩体试件，模型尺寸为10.05 m×5.0 m×5.0 m(长×宽×高)。在试件的高度方向 X = 5.0 m 处有一贯通型水平节理面。整个模图 1 厚壁圆筒问题有限元分析模型和网格示意图 Fig.1 Diagram of finite element model and mesh for thick cylinder problem型共划分 165 个单元，288 个节点，其中 150 个为 八节点空间实体等参单元，用于模拟岩石材料，15 个为八节点等厚度节理单元，用于模拟贯通型水平第 25 卷第 11 期谢学斌等. 岩石力学三维有限元分析的代数多重网格求解法? 2361 ?节理面，相应单元编号分别为 26～30，81～85 和136～140，如图 2 所示。在岩体试件的顶面 X = 10.05 m 处作用有垂直于该平面的均布压力，其大小为 98 MPa。采用的计算参数为：岩石材料 E =表3若干单元应力数值计算结果Table 3 Numerical results of several elements in finite element model单元号 1 16 21 26 30 注：单元 26，30 为节理单元。2.0 × 104 MPa ， μ = 0.25 ， 节 理 单 元 K n = 2 000 MPa/m， K s = 200 MPa/m。将模型中 X = 0 m 处底面所有节点法向约束， 其他各面所有节点自由变形， 模拟无侧限单轴受压情形。在上述计算条件下，利 用本文的代数多重网格有限元程序进行三维线弹性 分析，并与解析结果进行对比。在线弹性且不考虑 体力的条件下，模型中各单元应力及节点位移都可 以用解析法求出，其中岩石材料单元应力 σ X = 98σX ，σn /MPa98.013 4 98.013 5 98.013 5 98.013 2 98.012 3MPa，其他各应力分量为 0；节理单元应力为：法向应力 σ n = 98 MPa，切向应力 τ s = 0，各节理单元 顶底面上对应点在法线方向(整体坐标 X 方向)的相 对位移为 0.049 m。有限元数值结果应力分布特点 为：所有各单元应力分量 σ X 或 σ n 均接近 98 MPa， 其他各应力分量都接近 0。为节省篇幅，本文仅列 出模型中不同位置的 3 个岩石单元的 σ X 和 2 个节 理单元的 σ n 的数值计算结果，详见表 3。不同位置 节点 X 方向的位移 U X 数值结果与解析解的对比详 见表 4。水平节理面表4若干节点 X 方向位移 UX 数值解与解析解对比 and analytical results for several nodesTable 4 Comparison of displacement UX between numerical节点号 13 25 31 37 43 55 67解析解/cm 0.98 1.96 2.45 7.35 7.84 8.82 9.80数值解/cm 0.980 1.960 2.450 7.351 7.841 8.823 9.8043.3 代数多重网格法数值效率比较 为便于说明代数多重网格求解器数值计算效 率，本文设计了均质各向同性线弹性条件下原岩自 重应力场计算的有限元算例。取有限元计算域为125 m×183 m×195 m(X ×Y×Z)，分别采用 3 种不同的非均匀网格剖分方案，以得到不同规模有限元方 程，详见表 5。在 Pentium III 微型机(内存 256 MB，CPU 为 550 MHz)的 Windows 环境下采用代数多重网格求解器完成了各方案的有限元分析，并同时测 试了完成各方案分析所需的 CPU 时间。所有各方案 的自重应力场分析结果均收敛于解析解，计算所耗CPU 时间详见表 5。为进行比较，对各方案同时采图2 节理岩体受均布压力作用有限元分析模型与网格图 Diagram of finite element model and mesh for rock sample with horizontal joint loaded by uniform external force 表5 Table 5 2 种不同求解方法消耗的 CPU 时间比较 Fig.2Comparison of used CPU time between direct method and algebraic multigrain method观察表 3， 中数值计算结果可知， 4 各单元应力 和节点位移均收敛于解析解，说明代数多重网格有 限元方法在求解含节理单元的岩土力学问题时具有 良好的数值收敛特性和较强的适应能力。求解方程数/个 10 080 12 600 25 280直接法用时/s 182.60 366.50 2 454.92代数多重网格法用时/s 6.07 8.22 27.66? 2362 ?岩石力学与工程学报2006 年用了直接法――改进的乔列斯基分解法(LU 分解)进 行了求解，所消耗的 CPU 时间见表 5。 分析表 5 可知：(1) 代数多重网格法所消耗的CPU 时间远小于直接法，直接法与代数多重网格法所耗时间之比分别为 30.08，44.58 和 88.75。(2) 随 着所求解方程组规模的扩大，代数多重网格法数值 效率优势更明显， 如当方程组个数从 10 080 个增加 至 25 280 个时，直接法求解消耗的 CPU 时间从182.60 s 急增到近 2 500 s，增加了 13 倍，而代数多重网格法求解 CPU 时间从 6.07 s 增加到 27.66 s， 只增加了 3.5 倍，而且代数多重网格求解消耗的单位：MPaCPU 绝对时间与直接法相差近 2 个数量级，直接法与代数多重网格法消耗的 CPU 时间之比则从约 30 倍增加到近 90 倍。 由于理论上多重网格法计算工作 量仅与未知数个数的一次方成正比，当求解方程组 规模增加到数十万或百万数量级时，将具有更大的 数值优势。因此，本文研制的三维代数多重网格求 解器用于分析大型、超大型岩石力学工程问题时具 有比较明显的优越性。 3.4 实际工程问题的模拟 为检验程序求解实际工程问题的能力，利用本 文的有限元程序对铜陵有色金属(集团)公司冬瓜山 铜矿床 I 号主矿体的开采过程进行了三维线弹性模 拟研究。 该矿体埋深 800～1 100 m， 走向长 1 800 m， 平均厚 30～50 m。工程地质条件良好，设计采用阶 段空场滞后充填采矿法。数值模拟时将模型简化成 顶板(大理岩)、矿体(含铜矽卡岩)和底板(粉砂岩)3 层多介质模型，采用非均匀网格剖分形成有限元模 型，采用的采场结构参数、岩体力学参数和原岩应 力场参见相关研究[10]。模拟研究中当两个空间采场 完全形成时，垂直采场长轴方向的典型剖面应力分 布如图 3～5 所示。 从应力分布图可知， 最大主应力 的分布规律为：在采场的四角和两个已采矿房之间 的矿柱为应力集中区，在采场的高度方向，中间部 位为应力降低区，在远离采场的区域应力大小与原 岩应力场一致。最小主应力的分布规律是：采场高 度方向中间部位为应力降低区，还出现了局部拉应 力区，远离采空区的位置应力大小与原岩应力场一 致。模拟研究得到的应力空间分布规律与采场地压 分布的一般规律完全一致，说明本文研制的基于代 数多重网格求解的三维有限元程序能够满足模拟岩 石力学工程问题的需要。图5 Fig.5 典型剖面最小主应力σ3 分布图 Contours of σ3 of a typical section单位：MPa 单位：MPa图3 Fig.3典型剖面最大主应力σ1 分布图 Contours of σ1 of a typical section图4 Fig.4典型剖面中间主应力σ2 分布图 Contours of σ2 of a typical section第 25 卷第 11 期谢学斌等. 岩石力学三维有限元分析的代数多重网格求解法? 2363 ?河海大学学报，)：9C17.(Shi Jinsong，Zai Junming.4结语Multigrid algorithm for finite element equations of elastic mechanical problems[J]. Journal of Hohai University，)：9C17.(in本文采用代数多重网格法求解岩石力学三维有 限元离散线性方程组。线弹性数值试验结果表明， 本文基于代数多重网格求解的有限元程序具有良好 的数值收敛特性和较强的适应能力。对于大规模的 离散方程组，所消耗的 CPU 时间远小于直接法，具 有很高的数值效率。同时由于采用了一维稀疏存储 技术，大量节省了计算机存储量，从而为在微机上 进行大规模和超大规模岩石力学工程问题的三维有 限元模拟提供了一种可行和快速的求解方法。 参考文献(References)：[1] 刘之行，封卫兵. 三维区域上的多重网格算法[J]. 西安交通大学学 报，)：94C97.(Liu Zhixing，Feng Weibing. Multigrid method to three-dimensional domains[J]. Journal of Xi′an Jiaotong [8] [7] [6] [5]Chinese)) 王建华， 殷宗泽， 赵维炳. 多重网格法在岩石力学与工程中的应用[J]. 岩石力学与工程学报，)：336C345.(Wang Jianhua，Yin Zongze，Zhao Weibing. The application of multigrid method to rock mechanics and engineering[J]. Chinese Journal of Rock Mechanics and Engineering，)：336C345.(in Chinese)) Stüben K. A review of algebraic multigrid[J]. Journal ofComputational and Applied Mathematics， 2001， 128(1/2)： 281C309. 谢学斌，肖映雄，潘长良，等. 代数多重网格法在岩体力学有限元 分析中的应用[J]. 工程力学，)：165C170.(Xie Xuebin， Xiao Yingxiong，Pan Changliang，et al. Application of algebraic multigrid method to finite element analysis of rock mechanics[J]. Engineering Mechanics，)：165C170.(in Chinese)) Chan T F，Xu J C，Zikatanov L. An agglomeration multigrid method for unstructured grids[J]. Contemporary Mathematics，： 67C81.University，)：94C97.(in Chinese)) [2] Parsons I D，Hall J F. The multigrid method in solid mechanics：part I―algorithm description and behavior[J]. International Journal of Numerical Methods in Engineering，)：719C737. [3] Parsons I D，Hall J F. The multigrid method in solid mechanics：part II―practical applications[J]. International Journal of Numerical Methods in Engineering，)：739C753. [4] 史金松，载俊明. 弹性力学问题中有限元方程的多重网格解法[J]. [9] Mandel J，Brezina M，Vanek P. Energy optimization of algebraic multigrid bases[J]. Computing，5C228. [10] 谢学斌. 硬岩矿床岩爆预测与控制的理论和技术及其应用研究[博 士学位论文][D]．长沙：中南工业大学，1999.(Xie Xuebin. Theory and technology of rockburst prediction and control for hard rock mines and its application[Ph. D. Thesis][D]. Changsha：Central South University of Technology，1999.(in Chinese))下期内容预告下期《岩石力学与工程学报》主要发表下列内容的文章： (1) 再论茅坪滑坡的复活机制与治理可行性； (2) 裂隙岩体应力渗流耦合模型在压力隧洞工程中的应用； (3) 基于坐标投影图解的结构面和块体的计算机描述； (4) 颗粒材料的离散颗粒及间隙水渗流模型与数值模拟； (5) 奉节白衣庵滑坡演化的工程地质与历史地质分析； (6) 节理岩石的应力波动与能量耗散； (7) 不同剪切速率下节理的强度特性研究； (8) 高地应力下地下工程稳定性分析与优化的局部能量释放率新指标研究； (9) 短文(研究进展与工程实录)。
赞助商链接
岩石力学有限差分方法,岩石力学离散元分析,岩石力学...{R} (10.8) 求解上式有限个线性代数方程组,得到...et 和 Rizze 分别提出的直接和间接两种数值解法。...5.3 三维静磁场的有限元分析 5.3.1 边值问题 ...分析磁场的旋度旋度方程,单元中每一代数方程的变量数...按不同的解法, A 的数值可能不同,即体现了 A ...分方法对三维声波方程进行离散,推导 了三维和 k 维...代数多重网格 1 引言 波动方程是一种主要描述声波,...有限元法则是电子计算机时代的产物,有多种专门的应用...三维时变、 高频电磁场分析 土壤和岩石中的 非定...一种在力学模型上进行近似数值计算的方法,所求得的...参量的代数方程组, 求解此离散方程组就得到有限元法...声场和热场分析于一体的大型大型 通用有限元分析软件...(波前、 稀疏矩阵) 、特征值求解法(分块 Lanczos...(分布式并行、 代数多重网格)等,用户可根据问题类型...固体力学的有限元分析中, 是最简单的三维有限元单元...量的代数方程组,进而借助矩 阵和计算机求解代数方程...有限元法的基本思想是 2、网格划分时,可以采用 成...通常,这些反问题的数值解法是计算能力的要 求,特别是当问题必须制定在三维上。...应用代数 多重网格算法的逆生物电场问题的有限元法的制定等。 制定了多重网格...有限元网格剖分采用四边形八节点等参单元,将边坡中...整治分析 中大多仅用它作为初步分析和估计,有限元...边坡稳定性三维有限元分析[J].岩石力学与工程学报,...(代数多重网格) 直接根据离散后的 Matrix 生成&...对 三维 27 点格式有限差分离散的 Poisson 方程,...有限元分析,也称 FEA,它把结构分解成离散的单元,...的代数方程组, 解这些方程组就可以求出物体上有限...成网格,这些网格称为单元;网格间互 1 《有限元法...维应力分析 黏弹性力学 二维与三维应力 分析;填筑和...
All rights reserved Powered by
www.tceic.com
copyright &copyright 。文档资料库内容来自网络，如有侵犯请联系客服。CFD中的多重网格求解器
CFD中的多重网格求解器
，Update：增加C++求解代码；
众所周知，在对控制方程进行离散后，会产生一个离散的数学方程组，这个方程组通常为稀疏的。其大小和复杂性跟维度、方程个数、离散格式有关。虽然我们可以使用任意一个方程组求解技术来求解这个方程组，但是目前的限制主要为计算机能力。目前的求解方法主要有直接求解法和迭代求解法。直接求解典型的方法为高斯消去法。在使用高斯消去法的时候，需要对求解系统进行$N*N*N$次操作（$N$为未知数数量），并且需要利用内存存储$N*N$个系数。迭代方法即为通过一系列迭代的方法获得最后的解。典型的迭代方法为Jacobi和Gauss-Seidel点迭代法。在迭代求解中，一般在每次迭代中操作$N$次即可，但是总的迭代数是不确定的。
目前在大部分CFD软件中，均采用迭代求解法来求解方程离散后的稀疏线性系统。然而，Jacobi和Gauss-Seidel点迭代法的思想虽然很简单，但是在处理大型方程组的时候收敛速率是非常慢的。其可以较快的消除短波误差分量，但是对长波误差分量的消除非常缓慢，因此他们并不适用于CFD求解。得益于多重网格求解技术，Gauss-Seidel点迭代法目前也在商业CFD软件中广泛应用。
多重网格方法的求解思想为首先在较细的网格上进行迭代来消除短波误差分量，然后再在粗网格上迭代消除长波误差分量。通过这种在不同粗细的网格上消除不同波长的误差分量的方法，可以快速的把各种波长的分量消除。在每个不同粗细量级的网格上，可以使用任何一种迭代求解技巧。
本文从一个包含三个未知数的方程组开始，简单介绍Gauss-Seidel点迭代法，然后介绍多重网格求解思想的基本步骤。最后使用基于几何的多重网格方法（Geometric multigrid, GMG）求解一个包含20个未知数的方程组并和传统的Gauss-Seidel点迭代法进行收敛速度对比。
一、Gauss-Seidel点迭代法
Gauss-Seidel点迭代法在大部分矩阵、数值分析教科书中均有提及。在这里用一个说明性的例子阐述Gauss-Seidel点迭代法的基本思想。假设我们有如下的方程组：
\begin{matrix}\tag{1}
2x_1+x_2+x_3=7 \\
-x_1+3x_2-x_3=2 \\
x_1-x_2+2x_3=5
\end{matrix}
对方程移项：
\begin{matrix}\tag{2}
x_1=(7-x_2-x_3)/2 \\
x_2=(2+x_1+x_3)/3 \\
x_3=(5-x_1+x_2)/2
\end{matrix}
假定初值为：
\begin{matrix}\tag{3}
x_1^0=0 \\
x_2^0=0 \\
\end{matrix}
将方程（3）中的$x_1^0$代入到方程（2）中有：
\begin{matrix}\tag{4}
x_1^1=(7-x_2^0-x_3^0)/2=(7-0-0)/2=3.5
\end{matrix}
将方程（4）求解的$x_1^1$和方程（3）中的$x_2^0$代入到方程（2）中有：
\begin{matrix}\tag{5}
x_2^1=(2+x_1^1+x_3^0)/3=(2+3.5+0)/3=1.8333
\end{matrix}
将方程（5）求解的$x_2^1$和方程（4）求解的$x_1^1$代入到方程（2）中有：
\begin{matrix}\tag{6}
x_3^1=(5-x_1^1+x_2^1)/2=(5-3.5+1.7
\end{matrix}
这样，我们有Gauss-Seidel点迭代法求解的第一次初始值：
\begin{matrix}\tag{7}
x_1^1=3.5 \\
x_2^1=1.8333 \\
x_3^1=1.6667
\end{matrix}
在接下来的迭代过程中，$x^1$的值即可当做初始值进行下一次迭代。针对方程(1)，10次迭代业已收敛。
二、多重网格求解思想
多重网格求解技术需要在不同粗细量级的网格上对同一个原始的PDE方程进行求解。考虑最细的原始网格均匀分布且节点距离为$h$，假设粗网格节点距离为$2h$，$4h$，$8h$...。首先我们需要在间距为$h$的网格上进行若干次迭代，然后跳跃至$2h$的网格上进行下一轮迭代。如图所示，最终在$h$网格上获得一个加速收敛的解。其具体步骤为：
图1. 多重网格V-cycle示意图，其中$h$为网格节点距离，横线表示不同粗细量级的网格，圆圈内的数字表示每层网格上的迭代数
1. Start multigrid iterations：
\begin{matrix}\tag{8}
\end{matrix}
其中$A_h$表示在间距为$h$的网格上对原始PDE离散后的矩阵系数。此处可以使用任意一种迭代技术求解，如使用Gauss-Seidel点迭代法迭代10步。在此步，可以快速的消除短波误差分量。迭代后我们有解$y_h$，以及残差$r_h$。
2. Restriction：
2.1:基于间距为$h$的网格上的残差$r_h$，需要通过Restriction对间距为$2h$的网格上的残差进行操作。一种简单的做法即为平均插值方法或带权插值方法（非均分网格）。需要注意的是，$r_h$共有$N$个分量，$r_{2h}$则有$2N$个分量。通过Restriction操作后，假设我们有3套粗网格，需要执行三次Restriction，最后有：$r_{2h}$，$r_{4h}$，$r_{8h}$。
2.2:在前一步，我们通过在间距为$h$的网格上对原始PDE离散后得到系数矩阵$A_h$，但是并不知道$A_{2h}$，$A_{4h}$，$A_{8h}$的具体值。多重网格可以使用基于几何（Geometric）或者基于代数（Algebraic）的方法来求出$A_{2h}$等矩阵系数。这俩种方法也即GMG和AMG。
2.3:在获得$r_{2h}$，$r_{4h}$，$r_{8h}$后，我们可以通过$A_{2h}$，$A_{4h}$，$A_{8h}$对误差$e_{2h}$，$e_{4h}$，$e_{8h}$进行迭代求解。通常并不需要很多步骤，简单几次即可。求解后我们有最终的$e_{2h}$，$e_{4h}$，$e_{8h}$。
2.4:通过误差$e_{2h}$，$e_{4h}$，$e_{8h}$对残差$r_{2h}$，$r_{4h}$，$r_{8h}$进行更新有最后的$r_{2h}'$，$r_{4h}'$，$r_{8h}'$。
3. Prolongation：
3.1:基于$e_{8h}$，插值计算$e_{4h}'$。然后对$e_{4h}$和$e_{4h}'$加和我们有新的$e_{4h}''$。以$e_{4h}''$作为初值，使用任意的迭代技术对$A_{4h} e_{4h}''=r_{4h}'$进行光顺（smooth）求解后有光顺的误差$e_{4h}'''$。在下文中用$e_{4h}$代替$e_{4h}'''$。
3.2:同理，基于$e_{4h}$（也即上文的$e_{4h}'''$），重复步骤3.1，我们有光顺后的误差$e_{2h}'''$。在下文中用$e_{2h}$代替$e_{2h}'''$。
3.3:同理，基于$e_{2h}$，最终我们有$e_{h}'''$。
4. Correction：
计算最终的解：$y=y+e_{h}'''$。至此我们完成了多重网格的一次sweep。重新进行第1步。
需要提及的是，在Prolongation的过程中，通常对误差进行光顺（smooth），Gauss-Seidel是一个非常好的光顺器。下面我们使用基于几何的多重网格（GMG）方法求解一个包含20个未知数的方程组来详细解释多重网格方法的具体实现步骤。
代码下载：使用C++对一个简单的算例进行多重网格求解，直接编译即可。
代码解析待续