长江科学院院报

长江科学院院报 >

2025 , Vol. 42 >Issue 12: 160 - 169

DOI: https://doi.org/10.11988/ckyyb.20241090

水工结构与材料

基于任意网格的二维自适应分析及其优化

  • 宋文硕 , 1 ,
  • 苏海东 , 1, 2 ,
  • 颉志强 1, 2
展开
  • 1 长江科学院 材料与结构研究所, 武汉 430010
  • 2 水利部水工程安全与病害防治工程技术研究中心, 武汉 430010
苏海东(1968-),男,湖北武汉人,正高级工程师,博士,主要从事水工结构数值分析和计算力学方法研究。E-mail:

宋文硕(1996-),男,湖北恩施人,硕士研究生,主要从事计算力学方法研究。E-mail:

Copy editor: 刘运飞

收稿日期: 2024-10-25

  修回日期: 2024-12-13

  网络出版日期: 2025-12-11

收起

Two-dimensional Adaptive Analysis Based on Arbitrary Meshes and Its Optimization

  • SONG Wen-shuo , 1 ,
  • SU Hai-dong , 1, 2 ,
  • XIE Zhi-qiang 1, 2
Expand
  • 1 Material and Engineering Structure Research Department, Changjiang River Scientific Research Institute, Wuhan 430010, China
  • 2 Research Center of Water Engineering Safety and Disaster Prevention, Ministry of Water Resources, Wuhan 430010, China

Received date: 2024-10-25

  Revised date: 2024-12-13

  Online published: 2025-12-11

Fold

摘要

独立覆盖流形法采取“分区级数解”的思想,将计算域划分为若干分区,在各个分区内部采用多项式级数等完备级数直接逼近真实物理场函数;各分区之间采用窄条带连接,通过单位分解函数将各分区级数连接成整体近似函数,在收敛意义上具有C1连续性;计算网格具有任意形状、任意连接和任意加密的特性。基于计算网格的任意性,采用“一分为二”的网格分裂算法作为网格划分和细分的策略;基于收敛意义上的整体C1连续性,采用物理场导数的连续程度作为误差控制的指标。以上组成了基于任意网格的自适应分析策略。通过数值试验验证了这一策略的可行性。在此基础上提出优化方案:采用绝对误差指标控制近似函数导数的连续性,以简化误差判断;采用网格预划分的方式优化凹角局部区域的网格分布。最后,通过方孔、重力坝等算例进行验证,表明优化方案可以大幅减少网格数量,从而节省算力。

关键词： 网格自动划分; 自适应分析; 误差分析; 独立覆盖; 数值流形方法

本文引用格式

宋文硕 , 苏海东 , 颉志强 . 基于任意网格的二维自适应分析及其优化[J]. 长江科学院院报, 2025 , 42(12) : 160 -169 . DOI: 10.11988/ckyyb.20241090

Abstract

[Objective] This study aims to optimize the two-dimensional adaptive analysis strategy of the independent cover-based manifold method, focusing on addressing its deficiencies in error control and mesh distribution, thereby significantly enhancing computational efficiency and engineering practicality. [Methods] Based on the arbitrarily shaped and connected cover meshes of the independent cover-based manifold method, a “split-one-into-two” mesh splitting algorithm was employed for arbitrary refinement, and the degree of continuity of physical field derivatives was adopted as the error control indicator, forming an adaptive analysis strategy. An optimization scheme was proposed. 1) Adopting an absolute error indicator to replace the relative error indicator: the original relative error indicator tended to cause over-refinement in regions of minor stress and was overly sensitive in concave corner singularity regions. Using the absolute error indicator not only simplified the error judgment logic but also permitted larger error thresholds to be set near singular points such as concave corners, thereby effectively avoiding over-refinement. 2) Introducing a local mesh pre-partitioning and short strip elimination strategy: to address the issue of excessively high mesh density and irregular distribution in concave corner regions, a local pre-partitioning strategy was proposed, which pre-set the initial mesh in these regions by inwardly offsetting and reversely extending the edges of the concave corner. Simultaneously, an adjacent point merging algorithm was introduced during the mesh splitting process, which avoided the generation of extremely short connection strips and improved the conditioning of the system equations. [Results] Verification through two typical hydraulic structure examples, the square-hole and the gravity-dam model, demonstrated that the optimized scheme achieved a breakthrough improvement in computational efficiency. For the square-hole example, the original adaptive strategy generated 310 covers, corresponding to 6 520 degrees of freedom (DOFs). Under the same accuracy objective, the optimized scheme required only 59 covers and 933 DOFs. This represented a reduction of approximately 81% in the number of meshes and approximately 86% in DOFs. For the gravity-dam example, the original strategy generated 228 covers and 4 810 DOFs, whereas the optimized scheme required only 106 covers and 2 354 DOFs, achieving significant results of over 53% reduction in the number of meshes and 51% reduction in DOFs. The most notable achievement of the optimized scheme was in the effective suppression of mesh over-refinement near concave corner singularity regions. The calculation results demonstrated that the new strategy could generate more reasonable meshes, while ensuring computational accuracy, it substantially reduced the computational scale, and greatly enhanced the computational efficiency. [Conclusion] The proposed optimization strategy significantly enhances the efficiency of adaptive analysis while maintaining high accuracy. Through absolute error control and local mesh pre-partitioning, it effectively solves the problems of mesh over-refinement and unreasonable distribution near concave corner singularities, laying a foundation for subsequent three-dimensional adaptive analysis and engineering applications. Future research includes: criteria for selecting error thresholds and the highest order of cover series; further automating the local mesh pre-partitioning process to enable it to handle more complex geometries, ultimately achieving the goal of efficient and fully automatic simulation analysis for hydraulic structures.

Key words： automatic mesh generation; adaptive analysis; error analysis; independent covers; numerical manifold method

开放科学(资源服务)标识码(OSID):

0 引言

工程中的数值计算通常为求解各类偏微分方程,其中最基本、应用最为广泛的是基于连续介质理论(如连续介质力学)的偏微分方程[1]。在许多情况下,微分方程的解仅在受限区域内迅速变化,例如很多偏微分方程在凹角处具有奇异性[2];对于这样的问题,均匀网格建模的计算效率不高,而根据计算误差自动调整网格以匹配物理场的局部变化(即自适应分析)具有重要意义[3]。
有限元法[4-6]是目前最常用的数值方法。20世纪80年代,Babuška等[7-9]提出了基于有限元的自适应分析方法,通过误差估计指标来确定计算精度不满足要求的区域,并通过局部加密调整网格,以尽量少的自由度提高数值解的精度。共分为4步[10]:求解、误差估计、标记、网格调整。具体的,求解是指在当前的剖分网格上得到数值解;误差估计是指根据模型问题选出的误差估计指标,使用数值解来估计局部和全局误差;标记是指根据单元误差估计指标和标记策略,将需要加密的单元标记出来;加密是指根据网格加密算法,加密满足标记策略的单元,得到新网格。重复上述过程,就可以得到一系列逐步细化的网格,最后得到适应特定问题的自适应网格。
但由于有限元法的自身局限,在自适应分析方面仍存在一些困难[4-6,11]。首先是网格问题,网格有形状和连接方面的要求:网格形状要尽可能规则以获得理想精度;网格之间必须通过结点相连以保证连续性。以上问题带来了复杂求解域的网格剖分难题及自适应分析中的网格加密难题。其次是计算精度控制问题:数值计算常用完全多项式的阶次n来度量精度,但是有限元法中能实现n阶完全多项式的单元类型不多;有限元法常用的低阶插值方式,在凹角等奇异区由于物理场分布的复杂性,必须布置非常细密的网格;物理场的导数(如结构分析中的应力)精度不高等。
1991年,石根华、裴觉民将现代数学的流形思想引入数值计算[12-13],提出了数值流形法。在此基础上,2011年起,苏海东等提出了独立覆盖流形法[14-15]与分区级数解:采用任意形状和任意连接的分区网格(即所谓的“任意网格”),有望彻底解决网格剖分和网格加密难题;其中,多项式级数最常用,各分区的n阶多项式级数可以保证n阶精度,物理场导数随级数阶次升高而自然收敛。基于以上特点,采用一分为二的网格分裂算法进行任意的网格加密,以各分区之间的物理场导数的连续性程度作为误差控制指标,开展了自适应分析研究,甚至尝试了自动计算[15]。目前的计算程序可以求解常见的2阶偏微分方程,如弹性力学平衡微分方程、拉普拉斯方程、泊松方程等[16]。但自适应分析的计算效率还有待提高,主要是网格划分较多,特别是在凹角奇异区。
本文首先简要介绍独立覆盖流形法及其自适应分析的基本原理和步骤,然后给出自适应分析算例,重点针对凹角奇异区的问题,优化误差控制及网格布置方案,以提高自适应分析的效率,最后通过方孔和重力坝算例验证优化方案的有效性。

1 独立覆盖流形法的计算原理简介

如图1所示,独立覆盖流形法将计算域划分为若干分区,在分区之间采用窄条带连接,其中,图中的灰色区域称为独立覆盖,每个独立覆盖及其周边白色的条带合起来称为一个覆盖,独立覆盖即覆盖的独立区域。在各个覆盖上定义多项式级数等完备级数(称为覆盖级数):在独立覆盖i的内部,采用覆盖级数Vi作为近似函数V,直接逼近真实物理场函数;在各条带处采用单位分解函数 ${\varphi }_{i}$将相邻分区内的级数连接成整体近似函数,即近似函数 $V=\sum {\varphi }_{i}{V}_{i}$,从而实现近似函数的线性过渡。通过伽辽金法求得整体近似函数中的级数系数,形成所谓的“分区级数解”[17]。
图1 任意形状的覆盖系统

Fig.1 Cover system with arbitrary shapes

以图1中截取任意形状的分区(或覆盖)为例,如图2所示,讨论其收敛性[18]。设覆盖1上的真解为 ${V}_{1}^{\mathrm{*}}$,定义覆盖级数V1;覆盖2上的真解为 ${V}_{2}^{\mathrm{*}}$,定义覆盖级数V2;覆盖重叠部分的窄条带上的真解为V*。见图3[18]。在条带处沿其厚度方向建立一维局部坐标 $\xi$,单位分解函数 ${\varphi }_{i}$为一维有限元的线性形函数 ${\varphi }_{1}=\frac{1-\xi }{2},{\varphi }_{2}=\frac{1+\xi }{2}$。每个覆盖由于采用完备级数,根据经典的函数局部逼近理论,任意的连续函数在闭集的分区内都可以用完备级数逼近[19],因而随着级数项数的增加,显然有 ${V}_{1}➝{V}_{1}^{\mathrm{*}},{V}_{2}➝{V}_{2}^{\mathrm{*}}$,其中,在二者重叠的窄条部分,容易验证 $V={\varphi }_{1}{V}_{1}+{\varphi }_{2}{V}_{2}$满足C0连续性,考虑到真解的连续性有 ${V}^{\mathrm{*}}={V}_{1}^{\mathrm{*}}={V}_{2}^{\mathrm{*}}$,于是有: $\mathrm{V}={\varphi }_{1}{V}_{1}+{\varphi }_{2}{V}_{2}➝{\varphi }_{1}{{V}_{1}}^{\mathrm{*}}+{\varphi }_{2}{{V}_{2}}^{\mathrm{*}}={V}^{\mathrm{*}}$,即验证了物理场函数本身的收敛性。
图2 独立覆盖及其条带连接方式

Fig.2 Two independent covers and strip-connec-tion between them

图3 窄条区域内沿厚度方向的一维局部坐标[18]

Fig.3 One-dimensional local coordinate along thickness direction of strip[18]

进一步考查其导数的收敛性及连续性。由于采用了完备级数,在每个覆盖同样有 $\frac{\partial {V}_{1}}{\partial x}➝\frac{\partial {V}_{1}^{\mathrm{*}}}{\partial x},\frac{\partial {V}_{2}}{\partial x}➝\frac{\partial {V}_{2}^{\mathrm{*}}}{\partial x}$(以x方向的导数为例);其中,在覆盖重叠的窄条部分有 $\frac{\partial V}{\partial x}=\frac{\partial \left({\varphi }_{1}{V}_{1}+{\varphi }_{2}{V}_{2}\right)}{\partial x}={\varphi }_{1}\frac{\partial {V}_{1}}{\partial x}+{\varphi }_{2}\frac{\partial {V}_{2}}{\partial x}+\frac{{V}_{2}-{V}_{1}}{l}$,由于已经验证的物理场本身的收敛性和真解的连续性,有 $\frac{{V}_{2}-{V}_{1}}{l}➝\frac{{V}_{2}^{\mathrm{*}}-{V}_{1}^{\mathrm{*}}}{l}=0$,考虑到真解导数的连续性有 $\frac{\partial {{V}_{1}}^{\mathrm{*}}}{\partial \mathrm{x}}=\frac{\partial {{V}_{2}}^{\mathrm{*}}}{\partial \mathrm{x}}=\frac{\partial {V}^{\mathrm{*}}}{\partial x}$,于是有 $\frac{\partial {\mathrm{V}}_{}}{\partial \mathrm{x}}➝{\varphi }_{1}\frac{\partial {V}_{1}}{\partial x}+{\varphi }_{2}\frac{\partial {V}_{2}}{\partial x}➝{\varphi }_{1}\frac{\partial {{V}_{1}}^{\mathrm{*}}}{\partial \mathrm{x}}+{\varphi }_{2}\frac{\partial {{V}_{2}}^{\mathrm{*}}}{\partial \mathrm{x}}=\left({\varphi }_{1}+{\varphi }_{2}\right)\frac{\partial {V}^{\mathrm{*}}}{\partial x}=\frac{\partial {V}^{\mathrm{*}}}{\partial x}$,即验证了物理场导数的收敛性,且近似函数具有收敛意义上的C1连续性。
将上述讨论推广到图1的多个覆盖连接情况,其中,多个覆盖的连接区域中单位分解函数可以采用多边形有限元的形函数,或者分解成多个有限单元,可以得到同样的结论。但一般情况下,由于窄条带在多个覆盖连接区域的面积很小,单位分解函数没有严格定义的情况对计算结果影响很小,因此一般仅考虑覆盖之间的窄条带,且在计算程序中,条带是自动添加的[20]。
综上所述,在独立覆盖流形法中,随着覆盖级数的升阶,物理场本身及其导数均是收敛的。其中,整体近似函数的C1连续性作为收敛的一个必要条件,其严格程度反映了计算的精度,可以作为误差估计的一个指标。
同时,各分区或覆盖(以下两者等同)的形状和尺寸并不影响最终的逼近结果,而只可能影响逼近的快慢,即各覆盖可以是任意形状。另外如图1所示,各覆盖之间可以任意错位连接,展示了一种任意形状、任意连接的计算网格。由此得到覆盖网格的加密也具有任意性,即无论采用何种加密方式、独立覆盖细化成何种形状,以及细化后的覆盖之间如何连接,只需将各覆盖加密至适当大小,总可以做到在各覆盖内用适当阶次的完备级数逼近真实解[18]。这是该方法最基本的收敛原理。上述计算网格具有的任意形状、任意连接和任意加密的特性为自适应分析的网格加密提供了便利,采用一分为二的网格分裂算法并一直细分下去,这种任意的网格加密成为了一种可行的策略(由于条带为自动添加,下文关于计算网格的讨论均不考虑条带)。

2 基于任意网格的二维自适应分析方法

2.1 基于任意网格的自适应分析策略

目前,基于任意网格的h-p型混合自适应分析[21],其步骤为:
(1)对整个计算域进行凸分解,将凹角划分为凸体,如图4(a)的原计算域有一个凹角,可采用凹角面的双向切割方式,划分为图4(b)的形式(另一种方式为过凹角作最短垂线)。以此为计算网格进行一次计算,初始的多项式阶次为2阶。
图4 计算域网格划分及细分过程

Fig.4 Process of mesh partitioning and refinement in computational domain

(2)对计算结果进行误差估计,对不满足精度要求的分区网格进行标记。
(3)若标记分区的覆盖级数尚未达到预先设置的最高阶次(根据目前的经验,最高阶次常设置为4阶),对此分区进行升阶并再次计算。
(4)若标记分区的覆盖级数已达到预设的最高阶次,则需要对此分区进行网格加密,采用网格分裂算法[22],过形心作边界的最短垂线将该分区划分为2个分区(如图4(c)的网格是对图4(b)的网格1和网格2各加密一次的结果)。再从2阶开始计算。
(5)重复步骤(2)至步骤(4)步,直到满足精度要求为止。
以上计算流程中需要引入误差估计方法。考虑整体近似函数的C1连续性作为收敛的一个必要条件,包含对计算域内部相邻分区之间的连续性、计算域边界处与第二类边界条件的相符性的要求。对于结构分析而言,待评估的导数为分区边界面和计算域边界面上的法向正应力和切向剪应力。目前采用相对误差估计指标η[2,23],计算式为
$\mathrm{\eta }=\frac{{\left({\int }_{l}{(\mathrm{\epsilon }-\stackrel{~}{\epsilon })}^{2}\mathrm{d}l\right)}^{1/2}}{{\left({\int }_{l}{\stackrel{~}{\epsilon }}^{2}\mathrm{d}l\right)}^{1/2}}　。$
式中:ε为导数计算值; $\stackrel{~}{\epsilon }$为相对准确值;l为分区界面或边界面的尺寸。误差估计指标可以分为以下2种具体的情况:
(1)针对计算区域内部分区的边界面误差估计指标η1,如图5[21]所示,ε1和ε2分别为覆盖1和覆盖2的导数计算值(图中的黑线),ε取ε1和ε2中任意一个; $\stackrel{~}{\epsilon }=\frac{{\epsilon }_{1}+{\epsilon }_{2}}{2}$为两者在条形区的算术平均值。一旦η1超标,则对两侧网格中的较大面积分区进行标记,可以使网格细化趋向于均匀。
图5 相邻覆盖的ε分布[21]

(2)针对整个计算区域边界面为η2,其中,ε为边界面上的导数计算值; $\stackrel{~}{\epsilon }$为第二类边界条件的数值,若第二类边界条件的数值为0,则式(1)分母中的 $\stackrel{~}{\epsilon }$取为面内正应力[22]。一旦η2超标,则对相应分区进行标记。

Fig.5 ε distribution of two adjacent covers[21]

2.2 自适应分析算例

为了评估上述计算策略的有效性、精度与适用范围,以下依次进行3个数值算例的验证分析。第1个算例是无奇异点的流场计算,旨在验证该方法在无奇异工况下的基本收敛性与计算精度。然而,水工结构中广泛存在的计算域凹角点会引致奇异现象,对数值计算构成严峻挑战。为此,给出后续2个算例考察该上述策略在面对凹角奇异点时的表现。

2.2.1 不可压缩无旋流场

如图6(a)所示的两块平板之间的圆柱绕流算例[24],流场满足拉普拉斯方程,基于对称性,以圆柱中点为原点,建立如图6(b)所示的对称计算模型(其中c-g边为严格的圆弧)。
图6 圆柱绕流算例[24]

Fig.6 Example of flow around a cylinder[24]

采用前述的自适应分析流程,得到的最终计算网格如图7[16]所示,3个计算网格分别为采用最高阶次3、4、5的多项式覆盖级数在相同精度控制指标设置时自动生成的计算网格。计算得到的速度Vx的整体分布以及c-d边的分布与有限元参考解(划分2 000个8结点单元的细密有限元网格)的对比如图8和图9所示。为进一步对前述的自适应分析策略对误差阈值的敏感性进行初步探索,图10展示了不同误差阈值下c-d边上速度Vx与参考解的误差的变化规律(图中,v为自适应计算解,vh为有限元参考解)。
图7 流场分析的流形法网格[16]

Fig.7 Manifold method meshes for flow field analysis[16]

图8 流场分析的速度Vx分布

Fig.8 Distribution of velocity Vx in flow field

图9 c-d线上的速度Vx分布对比

Fig.9 Comparison of velocity Vx distributions along line c-d

图10 计算误差随误差阈值的变化规律

Fig.10 Variation of calculation error with error threshold

由图7可见,阶次越高,所需的网格越少;由图10可见,误差阈值越小,计算结果的误差越小。因此,覆盖级数最高阶次和误差控制阈值的设置需要综合考虑计算量的问题。根据研究经验,宜将级数最高阶次设置为4阶,下文算例以4阶为最高阶次;误差控制阈值的选取在不同类型的方程中差异较大,其取值规律有待进一步研究。本算例没有奇异区,在计算精度控制下可以得到稳定的计算结果。以下再讨论带有凹角奇异区的自适应分析算例。

2.2.2 方孔

中部开方孔的矩形板在上、下两端承受均布拉力p=10.00 kN/m2,采用如图11所示的1/4计算模型,左边和底部施加法向约束[21];按照前述的自适应策略,采用η1=η2=5%的误差控制指标进行h-p混合自适应分析。如图12所示,网格加密趋向于孔口的角点奇异部位,而在应力变化平缓的区域则是自动地采用大网格。选取孔口附近的σy、左侧边界(x=0 m处)和底部边界(y=0 m处)的法向应力σn,细密网格有限元分析与本文解的对比结果如图13、图14所示。
图11 方孔算例的1/4计算模型

Fig.11 Quarter computational model of square-hole example

图12 方孔算例自适应分析网格

Fig.12 Meshes for adaptive analysis of square-hole example

图13 孔口附近σy的分布

Fig.13 Distribution of σy near square hole

图14 方孔算例法向应力与细密网格有限元结果对比

Fig.14 Comparison of normal stresses with dense-mesh FEM results of square-hole example

图15 重力坝的计算模型

Fig.15 Computational model of gravity dam

2.2.3 重力坝

如图15所示的重力坝,坝高100 m,上游表面受水深为80 m的静水压力。坝体和基岩初始划分为2个分区,坝体弹性模量为26 GPa,泊松比为0.167,基岩弹性模量为30 GPa,泊松比为0.2。坝基底部全约束,上、下游侧法向约束。按照前述的自适应策略,采用η1=η2=5%的误差控制指标进行h-p混合自适应分析。如图16所示,网格加密趋向于坝踵和坝趾的角点奇异部位。选取全计算域的垂直向应力σy、上游面和坝基面的法向应力与细密网格的有限元分析结果对比,如图17、图18所示。
图16 重力坝算例自适应分析网格

Fig.16 Meshes for adaptive analysis of gravity dam example

图17 重力坝σy分布

Fig.17 Distribution of σy in gravity dam

图18 重力坝算例法向应力与细密网格有限元结果对比

Fig.18 Comparison of normal stresses with dense-mesh FEM results of gravity dam example

2.3 存在的不足

上节算例验证了基于任意网格的二维自适应分析策略的可行性,计算结果与参考解有很好的可比性。但仍存在误差控制和网格分布方面的问题,造成自适应划分的网格较多,影响了计算效率。
首先是误差控制问题。2.1节提出的相对误差指标,当物理场导数在某些局部区域的值很小,比如接近于0时,如果绝对误差稍大,则会得到较大的相对误差,比如0.02相对于0.01就有100%的相对误差,表现在某些应力较小的局部“过度加密网格”,而这些区域常常不是工程关注的重点,造成算力的浪费。另外在凹角区域,误差指标常常显得“过于严格”,原因是应力的角点奇异性会使得误差数值过大,造成网格不断加密而不能停止,如图19自适应生成的计算网格细节示意所示,2.2节给出的2个带有凹角的算例,自适应网格在凹角奇异区附近分布极细密且不规律,而局部应力奇异性的分布通常并不受关注,这也造成了算力浪费,将来开展三维分析更要考虑计算的效率,因此,不同部位应该分别选用合适的误差阈值。
图19 自适应生成的计算网格细节示意

Fig.19 Schematic details of adaptively generated meshes

其次是网格分布问题。2.1节提出的网格加密策略仅从待加密覆盖本身出发,并不考虑其与周围覆盖的连接关系,在加密过程中容易出现覆盖之间只通过很短的条带相连接的情况(例如图4(c)中2号覆盖和4号覆盖之间),使条带数目增多而增大计算量,且不利于线性方程组的性态。其次如图11的方孔算例,在对凹角进行凸分解时,由于凹角面过小,过早的凹角切割也会造成网格的过度加密。

3 自适应分析的优化方案

3.1 针对误差控制问题的优化方案

首先,采用绝对误差指标替代相对误差指标:
$\mathrm{\eta }=\frac{{\int }_{l}\left|\left.\epsilon -\stackrel{~}{\epsilon }\right|\mathrm{d}l\right.}{l}　。$
其次,如图20所示,对计算域不同区域设置不同的误差控制指标,特别针对凹角奇异区的很小范围内,可以设置较大的误差阈值,避免在凹角奇异区过度加密网格。
图20 误差指标的分区

Fig.20 Partition of error indicators

3.2 针对网格分布问题的优化方案

针对加密过程中容易出现覆盖之间通过极短条带相连接的问题,提出将相近点进行合并以消除短条带的方案。具体的,按照原策略加密某一覆盖时,过形心的最短垂线会与边界相交形成2个交点;然后检索交点所在的边上是否存在邻近点;若存在,就以此点代替交点;连接2个交点将待加密的覆盖一分为二,如图21所示的形成覆盖2和覆盖5的过程。
图21 合并相近点以消除短条带

Fig.21 Merging of adjacent points to eliminate short connection strips

其次,针对自适应网格在凹角奇异区域网格密度高、分布无规律的问题,提出局部网格预划分方案,如图22所示,分为两步:第一步,将凹角的两条边向计算域内部偏移一定的距离;第二步,反向延长凹角的边以消除凹角。
图22 局部网格预划分

Fig.22 Pre-partitioning of local meshes

3.3 优化方案的算例验证

以优化方案再次对第2节给出的方孔和重力坝算例进行自适应分析。

3.3.1 方孔自适应分析

优化方案的自适应网格(59个覆盖)相比于原策略的自适应网格(310个覆盖)大幅减少,如图23所示;优化方案的自由度(n=933)相比于原策略的自由度(n=6 520)大幅减少;计算效率显著提高。选取优化方案、原策略以及细密网格的有限元分析结果在左侧边界(x=0处)和底部边界(y=0处)的法向应力对比结果如图24所示,可见,网格和自由度虽然减少,但仍能保持很高的精度。
图23 方孔算例的优化网格

Fig.23 Optimized meshes for square-hole example

图24 优化前后的法向应力结果对比

Fig.24 Comparison of normal stresses before and after optimization

3.3.2 重力坝自适应分析

优化方案的自适应网格(106个覆盖)相比于原策略的自适应网格(228个覆盖)大幅减少,如图25所示;优化方案的自由度(n=2 354)相比于原策略的自由度(n=4 810)大幅减少;计算效率显著提高。全计算域的σy如图26所示,上游面和坝基面的法向应力与原策略以及细密网格的有限元分析结果对比如图27所示。优化方案在减少网格和自由度的同时仍能保持很高的精度。以上算例验证了优化方案的有效性。
图25 重力坝算例的优化网格

Fig.25 Optimized meshes for gravity-dam example

图26 优化方案的σy解

Fig.26 Solution of σy obtained from optimized scheme

图27 优化前后的法向应力以及与细密网格有限元结果对比

Fig.27 Comparison of normal stresses before and after optimization and with dense-mesh FEM results

4 结论

本文采取独立覆盖流形法这一新的数值方法,基于任意网格进行自适应分析,针对其自适应策略提出优化方案,可以得出以下结论:
(1)采取“分区级数解”思想的独立覆盖流形法在收敛意义上具有整体C1连续性,其计算网格具有任意形状、任意连接和任意加密的特性,为自适应分析提供了比有限元法更为有利的条件。
(2)采用“一分为二”的网格分裂算法作为网格划分和细分方法,采用物理场导数的连续性的程度作为误差控制指标,形成的自适应策略是可行的。
(3)采用绝对误差指标代替相对误差指标,误差控制问题就简化为不同区域的误差阈值的选取问题;针对凹角奇异区,采用局部网格预划分的方式可以优化网格分布,配合适当的误差阈值,从而使网格数量大幅减少而节省了算力。
这些工作为基于独立覆盖流形法的三维自适应分析打下了良好的基础。为了推动此方法走向实用性,有待进一步研究的内容包括:①选定合适的误差控制阈值以实现计算量和计算精度之间的平衡;②选定合适的覆盖级数最高阶次以实现覆盖级数阶次和覆盖数量之间的平衡;③针对各种凹角复杂情况,继续研究网格预划分的方式,最终实现全自动的网格优化;④研究h-p自适应的优化路径(即网格加密或级数升阶时机的合理选择),进一步提高计算效率。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
SU H, LIN S, XIE Z, et al. Fundamentals and Progress of the Manifold Method Based on Independent Covers[J]. Science China Technological Sciences, 2024, 67(4): 966-991.

DOI

[2]
郎美玲. 针对非线性椭圆型方程的自适应有限体积元法研究[D]. 济南: 济南大学, 2022.

(LANG Mei-ling. Adaptive Finite Volume Element Method for Nonlinear Elliptic Equations[D]. Jinan: University of Jinan, 2022. (in Chinese))

[3]
BRENNER S C, SCOTT L R. The Mathematical Theory of Finite Element Methods[M]. New York: Springer New York, 2002.

[4]
王勖成, 邵敏. 有限单元法基本原理和数值方法[M]. 2版. 北京: 清华大学出版社, 1997: 35-36.

(WANG Xu-cheng, SHAO Min. Basic Principle and Numerical Method of Finite Element Method[M]. 2nd Ed. Beijing: Tsinghua University Press, 1997: 35-36. (in Chinese))

[5]
ZIENKIEWICZ O C, TAYLOR R L, ZHU J Z. The Finite Element Method: Its Basis and Fundamentals[M]. Oxford,Tokyo: Butterworth-Heinemann, 2013: 43-79.

[6]
BELYTSCHKO T. 连续体和结构的非线性有限元[M]. 庄茁,译. 北京: 清华大学出版社, 2002: 4-11.

(BELYTSCHKO T. Nonlinear Finite Elements for Continua and Structures[M]. Translated by ZHUANG Zhuo. Beijing: Tsinghua University Press, 2002: 4-11. (in Chinese))

[7]
BABUŠKA I. Error Estimates for Adaptive Finite Element Computations[J]. SIAM Journal on Numerical Analysis, 1978, 15(4): 736-754.

DOI

[8]
BABUŠKA I, MILLER A. A Feedback Finite Element Method with a Posteriori Error Estimation: Part I. the Finite Element Method and Some Basic Properties of the a Posteriori Error Estimator[J]. Computer Methods in Applied Mechanics and Engineering, 1987, 61(1): 1-40.

DOI

[9]
BABUŠKA I, RHEINBOLDT W C. A Posteriori Error Estimates for the Finite Element Method[J]. International Journal for Numerical Methods in Engineering, 1978, 12(10): 1597-1615.

DOI

[10]
许雯. 二维椭圆偏微分方程基于高精度技术的自适应有限元方法[D]. 湘潭: 湘潭大学, 2021.

(XU Wen. High Accuracy Techniques Based Adaptive Finite Element Method for Two-dimensional Elliptic PDEs[D]. Xiangtan: Xiangtan University, 2021. (in Chinese))

[11]
WESSELING P. Principles of Computational Fluid Dynamics[M]. Heidelberg: Springer, 2001.

[12]
SHI G H. Manifold Method of Material Analysis[C]// Proceedings of the Ninth Army Conference on Applied Mathematics and Computing. Minneapolis, Minnesota, USA, 1991:51-76.

[13]
裴觉民. 数值流形方法与非连续变形分析[J]. 岩石力学与工程学报, 1997, 16(3):279-292.

(PEI Jue-min. Numerical Manifold Method and Discontinuous Deformation Analysis[J]. Chinese Journal of Rock Mechanics and Engineering, 1997, 16(3):279-292. (in Chinese))

[14]
祁勇峰, 苏海东, 崔建华. 部分重叠覆盖的数值流形方法初步研究[J]. 长江科学院院报, 2013, 30(1): 65-70.

(QI Yong-feng, SU Hai-dong, CUI Jian-hua. Preliminary Study on Numerical Manifold Method with Partially Overlapping Covers[J]. Journal of Yangtze River Scientific Research Institute, 2013, 30(1): 65-70. (in Chinese))

DOI

[15]
SU H D, QI Y F, GONG Y Q, et al. Preliminary Research of Numerical Manifold Method Based on Partly Overlapping Rectangular Covers[C]//DDA Commission of International Society for Rock Mechanics. Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD11), Fukuoka, Japan, August 27-29, 2013, London: Taylor & Francis Group, 2013: 341-347.

[16]
苏海东, 宋文硕, 龚亚琦, 等. 独立覆盖流形法的通用计算公式和通用程序设计: (三)算例验证[J]. 长江科学院院报, 2025, 42(4): 211-218.

DOI

(SU Hai-dong, SONG Wen-shuo, GONG Ya-qi, et al. General Formulas and Program Design for Manifold Method Based on Independent Covers Ⅲ: Example Verification[J]. Journal of Changjiang River Scientific Research Institute, 2025, 42(4): 211-218. (in Chinese))

[17]
苏海东, 韩陆超, 颉志强. 精确几何薄曲梁曲壳分析的分区级数解[J]. 长江科学院院报, 2022, 39(9):144-151.

DOI

(SU Hai-dong, HAN Lu-chao, XIE Zhi-qiang. Analysis of Thin Curved Beam and Curved Shell with Exact Geometry Using Piecewise-defined Series Solutions[J]. Journal of Yangtze River Scientific Research Institute, 2022, 39(9): 144-151. (in Chinese))

DOI

[18]
苏海东, 颉志强, 龚亚琦, 等. 基于独立覆盖的流形法的收敛性及覆盖网格特性[J]. 长江科学院院报, 2016, 33(2):131-136.

DOI

(SU Hai-dong, XIE Zhi-qiang, GONG Ya-qi, et al. Characteristics of Convergence and Cover Mesh in Numerical Manifold Method Based on Independent Covers[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(2): 131-136. (in Chinese))

DOI

[19]
吴顺唐. 数学分析[M]. 南京: 南京大学出版社,2022:29-46.

(WU Shun-tang. Mathematical Analysis[M]. Nanjing: Nanjing University Press, 2022: 29-46. (in Chinese))

[20]
苏海东, 杨震, 颉志强, 等. 独立覆盖流形法的通用计算公式和通用程序设计:(二)通用程序设计[J]. 长江科学院院报, 2025, 42(4): 202-210.

DOI

(SU Hai-dong, YANG Zhen, XIE Zhi-qiang, et al. General Formulas and Program Design for Manifold Method Based on Independent Covers Ⅱ. General Program Design[J]. Journal of Changjiang River Scientific Research Institute, 2025, 42(4): 202-210. (in Chinese))

[21]
苏海东, 龚亚琦, 颉志强, 等. 基于矩形独立覆盖初步实现结构静力分析的自动计算[J]. 长江科学院院报, 2016, 33(2): 144-150.

DOI

(SU Hai-dong, GONG Ya-qi, XIE Zhi-qiang, et al. Preliminary Implementation of Automatic Computation for Static Analysis of Structures Using NMM Based on Independent Rectangular Covers[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(2): 144-150. (in Chinese))

DOI

[22]
苏海东, 付志, 颉志强. 基于任意网格划分的二维自动计算[J]. 长江科学院院报, 2020, 37(7): 160-166.

DOI

(SU Hai-dong, FU Zhi, XIE Zhi-qiang. Automatic Two-dimensional Computation Based on Arbitrary Mesh Division[J]. Journal of Yangtze River Scientific Research Institute, 2020, 37(7): 160-166. (in Chinese))

DOI

[23]
刘亚军. 基于独立覆盖流形法的流体计算研究[D]. 武汉: 长江科学院, 2020.

(LIU Ya-jun. Research on Fluid Calculation Using Numerical Manifold Method Based on Independent Covers[D]. Wuhan: Changjiang River Scientiffic Research Institute, 2020. (in Chinese))

[24]
CHUNG T J. 流体动力学的有限元分析[M]. 张二骏, 龚崇准, 王木兰, 等, 译. 北京: 水利水电出版社, 1980: 203-213.

(CHUNG T J. Finite Element Analysis in Fluid Dynamics[M]. Translated by ZHANG Er-jun, GONG Chong-zhun, WANG Mu-lan, et al. Beijing: Water Resources and Electric Power Press, 1980: 203-213. (in Chinese))

Options
文章导航

/

〈 〉

*严禁使用本系统存储、处理、传输涉密信息*
版权所有 © 《长江科学院院报》编辑部
技术支持：北京玛格泰克科技发展有限公司