Medical Journal of the Chinese People's Armed Police Forces

Journal of Changjiang River Scientific Research Institute >

2025 , Vol. 42 >Issue 4: 202 - 210

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

Numerical Manifold Method Based on Independent Covers

General Formulas and Program Design for Manifold Method Based on Independent Covers Ⅱ: General Program Design

  • SU Hai-dong , 1, 2 ,
  • YANG Zhen 1 ,
  • XIE Zhi-qiang 1, 2 ,
  • QI Yong-feng 1, 2 ,
  • GONG Ya-qi 1, 2
Expand
  • 1 Material and Engineering Structure Department, Changjiang River Scientific Research Institute,Wuhan 430010, China
  • 2 Research Center of Water Engineering Safety and Disaster Prevention of Ministry of Water Resources, Wuhan 430010, China

Received date: 2024-02-02

  Revised date: 2024-05-06

  Online published: 2024-12-27

Fold

Abstract

Based on the general calculation formula of the manifold method based on independent covers presented in the previous article, we provide the flowchart of the calculation program. First, we summarize the integration methods for various geometric shapes (such as partitions, stripes, and boundary faces) that may appear in one- to three-dimensional spaces. On this basis, we develop integration programs according to simplex geometric elements of points, lines, faces, and bodies. This approach ensures the universality for any mesh shape. Next, we propose a programming strategy that separates the integration module from the integrand function module. The arbitrary combination of these two modules endows the program with extensibility and the potential to achieve universality in solving partial differential equations. Moreover, the universality of series is realized through the determination of series formulas, corresponding coordinates, coordinate transformation matrices, and series matrices. In addition, all calculation parameters can be input via formulas using user subroutines, thus achieving universality of input parameters. Ultimately, with relatively less program code, we can conduct one- to three- dimensional steady-state and transient analyses of the differential equations of motion in elasticity, conduction equations, and wave equations, including one to three types of boundary conditions.

Key words： partial differential equations; series solutions; mesh division; exact geometry; independent covers; numerical manifold method

Cite this article

SU Hai-dong , YANG Zhen , XIE Zhi-qiang , QI Yong-feng , GONG Ya-qi . 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 . DOI: 10.11988/ckyyb.20240113

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

0 引言

在《独立覆盖流形法的通用计算公式和通用程序设计——(一)通用计算公式》[1]中,将第i个覆盖对应于第mi项级数的形函数矩阵 T m i表达为单位分解函数φi、坐标转换矩阵Li和级数矩阵 N m i的乘积,具体讨论了各种情况下的形函数矩阵及其导数。对于弹性力学运动微分方程、传导方程、波动方程的一维至三维稳态和瞬态分析及其各类边界条件,由形函数及其偏导数矩阵形成单元质量矩阵、刚度矩阵和各种荷载向量,组集成整体矩阵后(瞬态分析再通过时间步的离散)形成线性方程组进行求解。
本文是独立覆盖流形法的通用计算公式和通用程序设计的第二部分——通用程序设计。首先介绍程序流程和程序设计思路,对一维至三维可能出现的各种几何形体(包括分区、条带和边界面,其含义见文献[1])的积分方式进行总结,提出积分子程序设计方案;然后将积分坐标转换成被积函数中的各种坐标,以计算被积函数并得到各种矩阵;通过裂纹尖端解析级数的计算流程给予具体示例;最后总结程序在通用性设计上的特点。

1 程序流程和设计思路

1.1 坐标定义

独立覆盖流形法适用于任意形状的覆盖网格,这意味着几何描述有一定的复杂性,而且程序的底层运行是由几何数据推动的,因此有必要重申一下描述几何数据的各种坐标,包括输入坐标 x ·、整体坐标x、参数坐标 x ∧、积分坐标x'、覆盖坐标 x ˙、级数坐标 x ¯:
(1)输入坐标 x ·是几何模型输入数据采用的坐标,包括通常情况下的整体直角坐标x(简称整体坐标),或用于曲线或曲面描述的参数坐标 x ∧。
(2)积分坐标x'是数值积分时采用的坐标,比如:在直线或平面边界所围成的分区积分中,通常为整体坐标;由分区界面(或本质边界条件所在的边界面)向内部扩展形成条带(或边界条带),这些面及其法向形成条带的积分坐标系,条带处的单位分解函数φi(x')由x'的法向坐标确定;其他如荷载作用的边界面、分区中的曲线(曲面)和其下面的直线(平面)所围区域等也要确定积分坐标系(见下文)。
(3)在各覆盖中定义覆盖坐标系 x ˙,最常见的是采用整体坐标x,也可能根据材料的方向、裂纹尖端所在覆盖的裂纹走向来定义。对于边界条带中的边界覆盖有 x ˙=x'。由 x ˙定义的Li( x ˙)矩阵用于位移的坐标转换。
(4)级数坐标 x ¯即级数矩阵 N m i( x ¯)中的坐标,如二维裂尖附近的解析级数采用的(r,θ)极坐标,多项式级数中将 x ˙归一化的坐标 x ¯= x ˙ - x ˙ i l x ˙ i,其中 x ˙ i为覆盖的中心坐标, l x ˙ i表示该覆盖在各方向的尺度。

1.2 程序流程

整个程序的计算流程如图1和图2所示。
图1 几何坐标及φi、Li和 N m i计算流程

Fig.1 Geometric coordinates and calculation process of φi, Li and N m i

图2 计算流程

Fig.2 Computational flow chart

首先由图1可知,由输入坐标 x ·形成一维至三维的各种几何体,通过数值积分得到积分点处的积分坐标x',转换到覆盖坐标 x ˙,进而得到级数坐标 x ¯,分别用于计算单位分解函数φi、坐标转换矩阵Li和级数矩阵 N m i。再分别由x'、 x ˙和 x ¯与整体坐标x的关系,计算 ∂ φ i ∂ x、 ∂ L i ∂ x和 ∂ N m i ∂ x。
再进入图2,从左上角框开始,此框向下,计算 ∂ T m i ∂ x= ∂ ∂ x(φi Li N m i)(关于 ∂ ∂ x(Li N m i)各种情况的讨论见文献[1]),放入应变子矩阵 B m i相应位置(A为微分算子矩阵),然后计算刚度子矩阵 K m i m j,形成单元刚度矩阵Ke。图中,Ω、Γ分别表示区域和边界。
左上角框向右,得到形函数子矩阵 T m i,然后分别计算质量子矩阵 M m i m j,各种荷载子向量 F m i,或第三类边界条件的刚度子矩阵 K m i m j,形成单元质量矩阵Me、刚度矩阵Ke和荷载向量Fe。图中ρ为密度,D和Db为参数矩阵,f和q分别为体荷载(或点荷载)和面荷载,L为边界面上的坐标转换矩阵,u'0 为边界位移,具体表达式见文献[1]。
单元矩阵和向量组集成整体矩阵和向量K、M、F和u,进入右下角的几个框,形成线性方程组进行求解,其中,瞬态分析考虑时间步的离散,Δt为时间步长,阻尼矩阵C=αM+βK,α、β和θ为参数,u0、 ∂ u 0 ∂ t、F0为前一步的物理量、时间导数和荷载,F1为当前步的荷载,u1为待求的当前步的物理量,具体见文献[1]。

1.3 总体思路

为了程序设计的通用性和简便性,本文提出将积分模块与被积函数模块分开考虑的思路:针对各种几何形体(包括一维至三维各种形状的分区、条带和边界面)确定不同的数值积分方案;根据偏微分方程的各项组成(如2阶、1阶、0阶偏微分项和非齐次项)及不同的级数类型,计算各种矩阵向量中的被积函数;两者之间任意组合——在积分点处调用被积函数,就可用较短的程序达成各式各样的功能。
对于积分模块,为适应复杂几何形状,以实现一维至三维的各种分区、条带及边界面的积分,将积分对象以点、线、面、体的几何元素表达,并拆解成单纯形(0维点、1维线、2维三角形、3维四面体)或单纯形与法向线段的组合,采用相对应的数值积分方案,得到积分点处的积分坐标x',然后转换成覆盖坐标 x ˙和级数坐标 x ¯,供被积函数模块调用。
对于被积函数模块,根据不同的级数类型,按单位分解函数、坐标转换矩阵和级数矩阵进行分类,计算各种矩阵和向量,主要内容见文献[1]。
图2涵盖了文献[1]的主要公式,对整个计算流程都有清晰的表述。文献[1]详述了各公式的计算细节,本文主要讨论图1中的2个细节:①如何由输入坐标 x ·形成各种积分形体,根据不同的积分方式得到积分坐标x';②如何将积分坐标x'转换到覆盖坐标 x ˙和级数坐标 x ¯;再加上程序在单元矩阵组装、公用数据模块和计算参数输入等方面的安排。这些内容将在以下各章节依次讨论。

2 各种几何形体及其积分方式

1.2节中的单元矩阵和向量积分,分为区域积分和边界积分两类,其中区域积分又分为2个步骤:先不考虑条带,按φi进行“分区积分”;然后再自动添加条带,考虑不同的φi进行“条带积分”。以下就分区、条带和边界3种类型,结合计算模型的输入数据,讨论可能出现的各种积分形体及其积分方式。

2.1 分区积分

任意形状分区积分的主要手段是拆解成单纯形,分别为一维线段、二维三角形、三维四面体。对于常用的整体坐标系下定义的多项式级数,可采用单纯形精确积分法[2-3]。本文主要讨论数值积分法,采用统一的积分公式[4],即
∫ Ω F d Ω = ∑ j n j J H j F j 　 。
式中:F为被积函数;nj为积分点个数(具体的积分点数根据积分精度确定,下同);J为雅克比系数(为长度、面积或体积,均考虑有向性[2]);Hj为积分点的权值(查数值积分表[4]);Fj为积分点上的被积函数值。
计算模型的几何输入坐标为 x ·,通常是整体坐标x,也可能是曲线或曲面的参数坐标 x ∧。一般情况下,单纯形(包括在参数坐标 x ∧投影面上建立的单纯形)内某个积分点处的积分坐标,即
x ' = ∑ k n d i m + 1 ξ k x · k 　 。  
式中: x · k为单纯形第k个顶点(或端点)的输入坐标;ξk为积分点位置的局部坐标(查数值积分表[4]);ndim为维数。

2.1.1 一维分区

如图3(a)所示,一维分区为一条直线(一维单纯形),计算模型的输入数据为2个端点的一维坐标x=(x)。式(1)采用一维Gauss积分公式[4],J为线段长度的一半。
图3 一维分区及条带示意图(点为分区的界面)

Fig.3 One-dimensional partition and strips (points represent the interfaces of partitions)

2.1.2 二维分区

考虑如图4(a)所示的包含曲线边界的二维分区。计算模型的输入数据为依次连接的顶点坐标,其中,曲线段的端点采用曲线的一维参数坐标 x ∧=( x ∧),其他端点采用整体二维坐标x=(x,y)。
图4 二维分区

Fig.4 Two-dimensional partition

如图4(b)所示,先不考虑曲线及其下面的直线所围区域,进行多边形积分:按顶点顺序依次连接形成三角形(二维单纯形);在每个三角形中,采用式(1)的二维Hammer积分公式[4],J为三角形的面积,图中显示了某个三角形中的4个积分点。
如图4(c)所示,对于曲线和直线所围区域,建立积分坐标系x'=(x',y')(图中表示为t-r)。在此坐标系下,积分由 ∫dS和 ∫dh两部分合成(即一维单纯形+法向线段),即
∫ Ω F d Ω = ∫ ∫ Ω F d S d h = ∑ j s n s ∑ j h n h J s H j s J h H j h F j 　 。
直线段 ∫dS采用一维Gauss积分(图中可见积分点个数ns=6),每个积分点再沿法向进行 ∫dh一维Gauss积分(图中可见积分点个数nh=2),积分上限落在曲线上(采用二分法求得曲线与法向线段的交点[5]),从而实现精确几何边界的模拟。雅克比系数Js和Jh为各自线段长度的一半, H j s和 H j h为相应的一维各积分点的权值。按 ∫dS和 ∫dh 2个方向的一维线段,分别根据式(2)计算积分坐标x'。

2.1.3 三维分区

考虑如图5(a)所示的包含曲面边界的三维分区。计算模型依次输入描述分区轮廓的平面和曲面,每个面为依次连接的顶点坐标,其中,平面的各顶点采用整体三维坐标x=(x,y,z)输入,曲面的各顶点采用曲面的二维参数坐标 x ∧=( x ∧, y ∧)。
图5 三维分区

Fig.5 Three-dimensional partition

与二维分区类似,先不考虑曲面及其下面的平面所围区域,进行多面体积分:如图5(b)所示的三维多面体区域,将形心依次与各面连接形成锥体;每个面按图4(b)的方式划分为三角形,与形心连成四面体(三维单纯形),采用式(1)的三维Hammer积分公式[4],J为四面体的体积,图5(b)中显示了某个四面体中的4个积分点。
如图5(c)所示,类似的按式(3)进行曲面和平面所围区域的积分(即二维单纯形+法向线段):以平面及其法向建立积分坐标系x'=(x',y',z'),其中选取一条线作为x'轴,在平面上按图4(b)的方式划分为三角形,采用二维Hammer积分公式进行 ∫dS积分(图中可见4个积分点),每个积分点再沿法向进行 ∫dh积分,采用一维Gauss积分公式(图中可见2个积分点),积分上限落在曲面上(采用二分法求得曲面与法向线段的交点),从而实现精确几何边界面的模拟。积分坐标也分为 ∫dS和 ∫dh 2个方向按式(2)分别计算。

2.2 条带积分

如图1(b)、图6(a)和图7(a)所示的一维至三维的一般条带(仅考虑两个覆盖之间的窄条带),计算模型的输入数据为分区的交界面,一维为点,二维为线,三维为面,图中显示为直线或平面,但可以是曲线或曲面,输入方式与2.1节中的点、线、面相同。
图6 二维条带

Fig.6 Two-dimensional strips

图7 三维条带

Fig.7 Three-dimensional strips

从分区的交界面沿法向朝两侧扩展,以图7(a)的三维条带为例,每一侧沿交界面及其法向建立积分坐标系x'=(x',y',z'),图中表述为t-s-r,其中选取一条线作为t轴,分界面的法向为r轴,并定义原点。两侧厚度分别为h1和h2(称为半条带),由相应的分区面积或体积确定,以保证条带占有整个分区面积或体积的一定比例(如1%)。每部分按式(3)进行积分(即单纯形+法向线段),其中 ∫dS为交界面上的积分, ∫dh沿法向采用等厚度的一维Gauss积分。
对于 ∫dS:一维条带 ∫dS=1;二维条带为线积分,采用一维Gauss积分(曲线投影到参数坐标 x ∧下积分,Js还要考虑映射关系);三维条带为面积分,按图4(b)划分为三角形,采用二维Hammer积分(曲面投影到参数坐标 x ∧=( x ∧, y ∧)下积分,Js考虑映射关系),还可能存在如图7(a)所示的直线与曲线所围区域,仍然如图4(c)所示,考虑面内两个方向的一维Gauss积分。积分坐标分为 ∫dS和 ∫dh 2个方向按式(2)分别计算。
如图1(c)、图6(b)和图7(b)所示的一维至三维的边界条带,由边界面(一维为点,二维为线,三维为面,可以是直线或曲线、平面或曲面)向内侧扩展出h厚度的条带,按以上一般条带的半条带方式进行积分。
在条带的各积分点上,被积函数分为2部分:一部分是条带积分值,考虑单位分解函数φi(具体公式见文献[1])计算局部近似函数,另一部分是需要扣除的重复计算的分区积分值,按φi=1的独立覆盖方式计算。

2.3 边界积分

如图3(a)的灰点、图4(a)和图5(a)灰色线段代表的边界荷载面(不一定与分区的表面重合),输入方式与2.1节中的点、线、面相同。边界线(面)按式(4)进行边界积分。
∫ Γ F d Γ = ∑ j n j J H j F j 　 。
式中:Γ为线或面的积分域,分别映射到线或面的积分坐标x'=(x')或x'=(x',y')上,应用一维Gauss线积分公式,或按图4(b)的方式划分为三角形采用二维Hammer积分公式,曲线和曲面情况的J考虑映射关系。

2.4 小结

综上所述,关于各种积分方式,需要掌握以下几条原则:
(1)几何形体的输入方式为其边界描述,一般为边界顶点坐标 x ·,包括整体坐标x或参数坐标 x ∧ 2种输入方式,后者表示边界为曲线或曲面,还需输入方程。
(2)分区、边界面拆分为单纯形,按式(1)或式(4)进行数值积分。积分坐标x'按式(2)计算。
(3)条带,分区中的直线与曲线、平面与曲面所围区域,需要建立积分坐标系x',采用单纯形+法向线段(考虑不同的线段长度)的方式,按式(3)分为∫dS 和 ∫dh 2个方向合成积分,积分坐标也分为2个方向按式(2)分别计算。对于曲线或曲面,积分公式中的J还要考虑映射关系。

3 以单纯形几何元素为基础的积分子程序

各类几何形体,其基本的几何元素是点、线、面、体,拆分成单纯形(分区)或单纯形+法向线段(条带或曲线+直线、曲面+平面所围区域),通过点、线、面、体4个最基本的单纯形积分子程序来实现第2节的所有积分,并分多种情况输出积分坐标x'。图8列出了各个子程序的积分类型,简要说明如下:
图8 由基本几何元素扩展成各类积分

Fig.8 Expanding from basic geometric entities to various integrals

(1)点、线、面程序,由维数ndim确定顶点坐标,线、面还分为直线(平面)和曲线(曲面)两类,对应于不同的J。然后分为直接输出积分坐标、作为界面再加上法向一维线段扩展(包括等厚度或变厚度)后输出积分坐标的两类过程。
(2)直接输出积分坐标的过程,用于模拟一维分区(线,如图3(a)所示),二维分区(面,如图4(b)所示),或计算一维至三维的边界荷载(如图3(a)、图4(a)和图5(a)所示的灰色点或线段),或计算内部刚度或荷载(如内部钢筋或钢板,传导方程的汇或源)。
(3)点、线、面作为分区界面或本质边界条件作用面,分为一维(见图3(b)和图3(c))、二维(图6)、三维(图7)情况,在界面上得到积分点后,由法向扩展成等厚度的条带,再输出法向一维线积分上的积分坐标。
(4)对于二维分区中的曲线和直线所围区域,如图4(c)所示,三维条带所包含的曲线和直线所围区域,如图7(a)所示,确定直线上的积分点后,再沿法向采用变厚度的方式进行一维积分并输出积分坐标。对于三维分区中的曲面和平面所围区域,如图5(c)所示,确定平面上的积分点后,再沿法向采用变厚度的方式进行一维积分并输出积分坐标。
(5)体程序作为三维分区,如图5(b)所示,直接输出整体坐标系下的积分坐标。

4 坐标转换及示例

4.1 坐标转换

第2节和第3节的各种积分,只为得到积分坐标x',然后转换成覆盖坐标 x ˙、级数坐标 x ¯后,供被积函数调用。图1的坐标转换具体流程如图9所示,说明如下:
图9 坐标转换流程

Fig.9 Coordinate conversion flow chart

(1)最左侧的输入坐标框中,几何体的输入坐标用 x · d表示,可能是整体坐标xd,或曲线、曲面的参数坐标 x ∧ d。此框向右,在框①中分别取出积分域的第k个顶点坐标 x · k,按式(2)生成积分坐标x'。
(2)积分坐标x'按分区和条带分两类,生成用于过渡的整体坐标xt:分区框中,对于基本的单纯形有xt=x',或者如图4(c)或图5(c)情况的曲(曲线或曲面)+直(直线或平面)所围区域计算xt=L'x'+x0s,L'为图中的局部坐标系的坐标转换矩阵,x0s为局部坐标系的原点坐标;条带框中,对于某一侧为独立覆盖的情况,也按照xt=L'x'+x0s计算xt,L'为条带局部坐标系的坐标转换矩阵。
(3)过渡坐标xt通过框②生成覆盖坐标 x ˙, x ˙= L i T(xt-x0),其中Li为覆盖坐标系的转换矩阵,x0为其原点坐标,比如图10所示的裂纹尖端所在的覆盖,由裂纹走向及裂尖位置确定覆盖坐标系。
图10 裂纹尖端所在的覆盖(二维问题)

Fig.10 The cover containing crack tip (two-dimensional issue)

(4)覆盖坐标 x ˙通过框③分为两类,生成级数坐标 x ¯:一般的多项式级数情况,归一化生成 x ¯= x ˙ - x ˙ i l x ˙ i;如图10所示的裂纹情况,在覆盖坐标系下求得 x ¯=[r θ]T。
(5)条带框中,边界条带的边界覆盖,覆盖坐标 x ˙=x',在框③中归一化得到级数坐标 x ¯。其中,直边界条的xi为整体坐标输入,需要通过 x ˙ i=L'T(xi-x0s)转化到条带局部坐标。
(6)与条带框相关(向下)的是单位分解函数φi(x')及 ∂ φ i ∂ x ',再结合xt=L'x'+x0s得 ∂ x ' ∂ x(条带框向下),进一步得到 ∂ φ i ∂ x。
(7)与覆盖坐标框 x ˙相关(向下)的是覆盖坐标转换矩阵Li( x ˙),对于边界覆盖为Li(x'),只有在曲线或曲面情况下才计算 ∂ L i ∂ x ',结合条带框向下的 ∂ x ' ∂ x得到 ∂ L i ∂ x。
(8)与级数坐标框 x ˙相关(向下)的是级数矩阵 N m i( x ¯)及 ∂ N m i ∂ x ¯,结合框③向下得到的 ∂ x ¯ ∂ x ˙,继续向左结合由框②得到的 ∂ x ˙ ∂ x(对于边界覆盖,结合 ∂ x ' ∂ x)求得 ∂ N m i ∂ x。

4.2 计算流程示例

以二维裂纹分析的解析级数为例,具体说明计算流程。
如图10所示,在包含裂纹尖端的分区中,根据裂纹走向确定覆盖坐标系 x ˙=( x ˙, y ˙),与整体坐标之间的转换矩阵为Li,裂尖坐标x0=(x0,y0)作为原点。位移级数公式为[6]
$\tilde{\boldsymbol{u}}_{i}=\left\{\begin{array}{c}\tilde{u} \\\tilde{v}\end{array}\right\}_{i}=\sum_{m_{i}=1}^{n_{i}} \frac{1}{2} r^{\left(m_{i}-1\right) / 2}\left[\begin{array}{ll}f_{1}^{x}(\theta) & f_{2}^{x}(\theta) \\f_{1}^{y}(\theta) & f_{2}^{y}(\theta)\end{array}\right]_{m_{i}}\left\{\begin{array}{l}a \\b\end{array}\right\}_{m_{i}} 。$
式中:ni为级数项数,各f函数参考文献[6](类似的有三维裂纹分析公式[7]); N m i为非对角阵。级数坐标 x ¯=( x ¯, y ¯)在覆盖坐标系下定义为极坐标(r,θ)。独立覆盖单元的自由度个数为npde·ni,其中,npde=2。
由裂尖依次连接分区顶点形成单纯形(三角形,如虚线所示),采用二维Hammer积分,对于任一积分点,由三角形各顶点的坐标x结合式(2)得到积分点坐标x',为整体坐标xt表示。然后计算覆盖坐标 x ˙= L i T(xt-x0)。由极坐标(r,θ)与覆盖坐标 x ˙的关系求得级数坐标 x ¯及 ∂ x ¯ ∂ x ˙,代入式(5)得到级数矩阵 N m i( x ¯)并求得 ∂ N m i ∂ x ¯。再由 x ˙= L T i(xt-x0)得到 ∂ x ˙ ∂ x=Li,再结合 ∂ x ¯ ∂ x ˙得到 ∂ N m i ∂ x。考虑到Li为常量矩阵, ∂ L i ∂ x=0;φi=1, ∂ φ i ∂ x=0,则图2左上框的计算完成。图2中的剩余操作参考文献[1]进行。
图10中与其他覆盖相接的右侧边界,向内扩展成半条带,建立条带坐标系x'=(x',y'),x0s=(x0s,y0s)为坐标系的原点(可选为界面线段的中点)。条带内采用 ∫dS和 ∫dh双向一维Gauss积分,分别按式(2)得到积分点坐标x',其中,先确定 ∫dS的积分点位置x',再由此位置沿法向确定 ∫dh积分点位置y'。然后转换到整体坐标xt=L'x'+x0s。进一步转换到覆盖坐标 x ˙= L T i(xt-x0)之后同上操作。条带另一侧连接的其他覆盖,也是通过整体坐标xt转换到其覆盖坐标 x ˙。
以上步骤同时演示了在程序中加入新的级数类型的一般过程,需要确定以下内容:①以级数公式确定级数的自由度数、级数坐标 x ¯和级数矩阵 N m i( x ¯)(为对角或非对角矩阵);②覆盖坐标 x ˙和坐标转换矩阵Li( x ˙)(为常量或变量矩阵);③ x ˙、 x ¯及x之间的相互关系。

5 矩阵组装及数据模块

5.1 矩阵组装

首先讨论自由度编号。按照独立覆盖(包括边界覆盖)的编号,依次对自由度进行排序,形成整体自由度编号。其中,独立覆盖单元的自由度个数为npde· ni,npde为本文所讨论的几种偏微分方程组的方程个数。对于常用的多项式级数,其阶次pi与项数ni之间的关系是ni=[ ∏ j = 1 n d i m(pi+j)]/(ndim!)。条带单元的自由度为两侧覆盖的自由度之和。
在自由度编号过程中,本质边界条件通过边界条带的边界覆盖施加:在边界覆盖中去掉不含有法向坐标r的所有自由度,使之自动满足边界条件。
单元矩阵的计算:对于每个覆盖的每一项级数,首先计算子矩阵 T m i和 B m i,其维数分别为npde×npde和npde×nb,nb为2阶微分项相关参数矩阵(如弹性矩阵D)的维数;以 K m i m j= ∫Ω B m i TD B m jdΩ为例,是不同级数项的子矩阵 B m i之间的相互乘积,对于条带,还包括不同覆盖的每个级数项之间的乘积,然后根据mi级数项和mj级数项在单元自由度中的排序,放入单元矩阵的合适位置。
各单元矩阵的系数,根据在整体自由度排序的位置,放入整体方程组中。

5.2 公用数据模块与用户子程序

用户输入的计算模型数据,包括:①总体控制信息——独立覆盖(含边界覆盖)的个数,流形单元(含独立覆盖单元和条带单元)的个数、维数、偏微分方程个数等;②每个独立覆盖的形心坐标、半条带的厚度、级数阶次、材料参数号、边界覆盖的约束信息等;③每个流形单元相关的独立覆盖号,描述单元几何信息的边界顶点坐标及其坐标系号(整体坐标系号为0,特殊的覆盖坐标系、曲线或曲面坐标系需要用户依次输入)等;④定义了第二、第三类边界条件的边界几何信息及边界参数号;(5)材料参数和初边值条件参数。这些数据放入公用数据模块,贯穿于整个程序流程,供各个子程序(子过程)调用,一般不可修改。
为了程序的通用性和可扩展性,设立了用户子程序接口,用于输入以下内容:
(1)曲线、曲面方程。常用的几种类型,如圆弧或球面、柱面,多项式描述的曲线或曲面等,可以在输入文件中按指定类型输入其参数,如圆弧半径和圆心,多项式阶次及每项系数。也可以通过用户子程序,直接定义曲线或曲面方程。
(2)参数公式。程序的所有计算参数,如材料参数、初边值条件参数等,可以在输入文件中直接输入数值,或者表述为“fn”,通过用户子程序的第n个公式输入。用户子程序可调用公用数据模块,引入坐标、时间等信息输入这些公式。

6 结论

对一维至三维各种几何形体(包括分区、条带和边界面)的积分方式进行总结,提出基于点、线、面、体的单纯形几何元素实现上述各种积分的积分程序设计方案,并能通过用户子程序输入各种几何外形,实现了模拟任意几何边界的形状通用性。提出将积分模块与被积函数模块分开考虑的编程思路,积分模块考虑各种积分形体对象,而被积函数由偏微分方程的各项组成及级数类型确定,两者可任意组合,使程序具备了扩展性,有望实现偏微分方程求解的通用性。通过级数公式和相应的各种坐标以及坐标转换矩阵、级数矩阵的确定,实现了级数的通用性。可通过用户子程序进行所有计算参数的公式化输入,实现输入参数的通用性。
最终可用较少的程序代码(目前的主体程序不到3 000行Fortran代码),实现了弹性力学运动微分方程、传导方程、波动方程的一维至三维(维数的通用性)稳态和瞬态分析(含一类至三类边界条件),包含了高阶级数逼近、任意形状和任意连接的网格、精确几何边界的模拟及本质边界条件的准确施加、裂纹尖端附近解析级数的应用等独立覆盖流形法的特色功能。梁板壳的程序设计将另文介绍。
文献[8]将给出一维至三维的稳态和瞬态分析算例进行验证,涵盖了位移场、温度场、渗流场、声场、以势函数表示的静电场和不可压缩无旋流场。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
苏海东. 独立覆盖流形法的通用计算公式和通用程序设计:(一)通用计算公式[J]. 长江科学院院报, 2025, 42(4):193-201,210.

(SU Hai-dong. General Formulas and Program Design for Manifold Method Based on Independent Covers Ⅰ. General Formulas[J]. Journal of Changjiang River Scientific Research Institute, 2025, 42(4):193-201,210.(in Chinese )

[2]
SHI G H. Simplex Integration for Manifold Method, FEM, DDA and Analytical Analysis[C]//Proceedings of the First International Forum on Discontinuous Deformation Analysis and Simulations of Discontinuous Media. Berkeley, Califonia, USA, 1996:205-262.

[3]
LIN S, XIE Z. A New Recursive Formula for Integration of Polynomial over Simplex[J]. Applied Mathematics and Computation, 2020, 376: 125140.

[4]
ZIENKIEWICZ O C, TAYLOR R L. The Finite Element Method[M].5th Edition. Oxford, UK: Butterworth-Heinemann, 2000.

[5]
苏海东, 付志, 颉志强. 基于任意形状网格和精确几何边界的数值计算[J]. 长江科学院院报, 2020, 37(7):167-174.

DOI

(SU Hai-dong, FU Zhi, XIE Zhi-qiang. Numerical Computations Based on Cover Meshes with Arbitrary Shapes and on Exact Geometric Boundaries[J]. Journal of Yangtze River Scientific Research Institute, 2020, 37(7): 167-174.) (in Chinese)

DOI

[6]
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.

[7]
祁勇峰, 苏海东, 龚亚琦. 基于新型流形法的三维应力强度因子求解[J]. 土木与环境工程学报(中英文), 2021, 43(5):58-65.

(QI Yong-feng, SU Hai-dong, GONG Ya-qi. Computing 3D Stress Intensity Factors Based on New Manifold Method[J]. Journal of Civil and Environmental Engineering, 2021, 43(5):58-65.) (in Chinese)

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

(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)

[9]
SMITH I M, GRIFFITHS D V. 有限元方法编程[M]. 王崧, 周坚鑫, 王来, 等, 译. 北京: 电子工业出版社, 2003.

(SMITH I M, GRIFFITHS D V. Programming the Finite Element Method[M]. Translated by WANGSong, ZHOUJian-xin, WANGLai, et al. Beijing: Publishing House of Electronics Industry, 2003.) (in Chinese)

Options
Outlines

/

〈 〉

Copyright © Journal of Changjiang River Scientific Research Institute, All Rights Reserved.
Powered by Beijing Magtech Co. Ltd.