CN113761762A - 用于有限元数值模拟后验误差估计的平衡通量构造方法 - Google Patents

用于有限元数值模拟后验误差估计的平衡通量构造方法 Download PDF

Info

Publication number
CN113761762A
CN113761762A CN202110886391.1A CN202110886391A CN113761762A CN 113761762 A CN113761762 A CN 113761762A CN 202110886391 A CN202110886391 A CN 202110886391A CN 113761762 A CN113761762 A CN 113761762A
Authority
CN
China
Prior art keywords
kth
vector basis
construction method
finite element
unit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110886391.1A
Other languages
English (en)
Other versions
CN113761762B (zh
Inventor
韩峰
王建国
刘钰
付梅艳
闫凯
姚成宝
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwest Institute of Nuclear Technology
Original Assignee
Northwest Institute of Nuclear Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwest Institute of Nuclear Technology filed Critical Northwest Institute of Nuclear Technology
Priority to CN202110886391.1A priority Critical patent/CN113761762B/zh
Publication of CN113761762A publication Critical patent/CN113761762A/zh
Application granted granted Critical
Publication of CN113761762B publication Critical patent/CN113761762B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G06T17/205—Re-meshing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Computer Graphics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Complex Calculations (AREA)

Abstract

为克服平衡通量构造方法较为复杂、某些方法仅适用于二维问题的缺陷，本发明提出了一种用于有限元数值模拟后验误差估计的平衡通量构造方法，包括步骤：1)在集合Th中每个网格单元上定义矢量基函数；Th表示对计算区域Ω进行剖分所得到的网格单元的集合；2)在每个网格单元上计算所述矢量基函数的系数；3)在相邻网格单元公共边/面的两侧更新矢量基函数的系数；4)基于更新后的矢量基函数的系数和步骤1)定义的矢量基函数，计算平衡通量。本发明仅需要在Th中每个单元上求解一个三阶线性方程组和二阶线性方程组(d＝2)，或仅需要在Th中每个单元上求解一个四阶线性方程组和三阶线性方程组(d＝3)。

Description

用于有限元数值模拟后验误差估计的平衡通量构造方法
技术领域
本发明属于有限元数值模拟技术领域，具体涉及一种用于有限元数值模拟后验误差估计的平衡通量构造方法。
背景技术
在物理和工程上的偏微分方程的数值求解中，有限元方法发挥着重要的作用。它是当前数值模拟技术中一种重要的方法，被广泛应用于计算力学、计算流体力学、计算电磁学等科学与工程计算领域。
有限元方法的基本思想是对计算区域进行网格剖分，在每个网格上构造一个分片函数(多项式)，然后结合变分原理求解得到物理模型(物理问题的数学模型，通常表现为一组偏微分方程及初边值条件)的近似解，其实质是用有限维空间的离散解来逼近无穷维空间的连续解。由于有限元方法本质上是对物理模型的近似求解，因此近似解必然存在一定的误差(与物理模型的真解存在一定的差异)。误差的大小，决定了有限元数值解的可靠性和应用价值。
随着有限元方法数学理论的不断发展，逐渐建立了有限元方法的后验误差估计理论，并产生了多种后验误差估计方法，包括残量型估计、局部问题辅助型估计、分级型估计、重构型估计等。有限元后验误差估计方法能够对有限元数值解的误差给出量化估计，从而能够用于量化评价数值解的准确性。但是在上述方法给出的后验误差估计中，通常包含一个和网格正则性相关的未知常数，所以上述后验误差估计结果并没有直接用于数值解的精度估计，而是作为误差指示子更多的用于有限元网格自适应算法研究中。文献[1]提出了一种基于通量重构的完全可计算后验误差估计方法，解决了上述问题。
应用Vohralik提出的后验误差估计方法，必须构造一个平衡通量。Vohralik给出了一种平衡通量的构造方法。这种方法需要在每个节点的单元片(如图1所示)上求解一个局部Neumann混合有限元问题。从求解算法和程序实现两个方面来看，上述平衡通量构造都较为复杂。文献[2]给出了一种平衡通量的构造方法。这种通量构造方法也是在每个节点的单元片上构造平衡通量。且这种方法仅适用于二维问题。
文献[1]M.Vohralik,Guaranteed and fully robust a Posteriori errorestimates for conforming discretizations of Diffusion Problems withDiscontinuous coefficients,J.Sci.Comput.46(2011)397-438。
文献[2]R.Verfurth,Lecture notes:Adaptive finite element methods,URLhttp://www.ruhr-uni-bochum.de/num1/files/lectures/AdaptiveFEM.pdf,2015。
发明内容
为了克服平衡通量构造方法较为复杂、某些方法仅适用于二维问题的缺陷，本发明提出了一种用于有限元数值模拟后验误差估计的平衡通量构造方法，为计算平衡通量提供了一种新的技术途径，为进一步应用完全可计算后验误差估计方法估计工程应用问题的数值解的精度奠定了技术基础。
本发明的技术方案是：
用于有限元数值模拟后验误差估计的平衡通量构造方法，其特殊之处在于，包括以下步骤：
1)在集合Th中每个网格单元上定义矢量基函数；Th表示对计算区域Ω进行剖分所得到的网格单元的集合；
2)在每个网格单元上计算所述矢量基函数的系数；
3)在相邻网格单元公共边/面的两侧更新矢量基函数的系数；
4)基于更新后的矢量基函数的系数和步骤1)定义的矢量基函数，计算平衡通量。
若计算区域Ω被剖分为Ne个三角形网格单元，
步骤1)具体为：
1.1)令k＝1；
1.2)在第k个三角形网格单元上，定义和三边分别对应的矢量基函数γki(x)：
Figure BDA0003194352590000031
其中，aki表示第k个三角形网格单元K三个顶点的坐标，lki表示第k个三角形网格单元K三边的长度，S表示第k个三角形网格单元K的面积；
1.3)令k＝k+1，若k<Ne，则返回步骤1.2)；否则执行步骤2)。
进一步地，步骤2)具体为：
2.1)令k＝1；
2.2)根据第k个三角形网格单元K的信息，构造下述方程，求解方程得到在步骤1)中定义的第k个三角形网格单元K中与其三个边分别对应的矢量基函数的系数
Figure BDA0003194352590000032
i＝1,2,3；
Figure BDA0003194352590000033
2.3)令k＝k+1，若k<Ne，则返回步骤2.2)；否则执行步骤3)。
进一步地，步骤3)具体为：
3.1)令k＝1；
3.2)根据第k个公共边F两边三角形网格单元的信息，根据下式更新第k个公共边左边三角形网格单元KL,F中与边F对应的矢量基函数系数
Figure BDA0003194352590000034
和右边三角形网格单元KR,F中与边F对应的矢量基函数系数
Figure BDA0003194352590000035
的值，更新后的矢量基函数系数分别为
Figure BDA0003194352590000036
和
Figure BDA0003194352590000037
Figure BDA0003194352590000038
Figure BDA0003194352590000039
其中，SK,L和SK,R分别表示第k个公共边F左右两边三角形网格单元的面积，，
Figure BDA0003194352590000041
nF表示三角形网格单元KL,F在公共边F上的外法向单位矢量；
3.3)令k＝k+1，若k<Nf，则返回步骤3.2)；否则，执行步骤4)；Nf为计算区域Ω中公共边的总数。
进一步地，步骤4)具体为：
将步骤3)求得的更新后的矢量基函数系数带入下式计算得到平衡通量σh(x)：
Figure BDA0003194352590000042
其中，
Figure BDA0003194352590000043
Figure BDA0003194352590000044
为数值解在三角形网格单元K上的梯度。
若计算区域Ω被剖分为Ne个四面体网格单元，
步骤1)具体为：
1.1)令k＝1；
1.2)在第k个四面体网格单元上，定义和四个面分别对应的矢量函数
Figure BDA0003194352590000045
其中，aki表示第k个四面体网格单元K四个顶点的坐标，lki表示第k个四面体网格单元K四个面的面积，S表示第k个四个面网格单元K的体积；
1.3)令k＝k+1，若k<Ne，则返回步骤1.2)；否则执行步骤2)。
进一步地，步骤2)具体为：
2.1)令k＝1；
2.2)根据第k个四面体网格单元K的信息，构造下述方程，求解方程得到在步骤1)中定义的第k个四面体网格单元中与其四个面分别对应的矢量基函数的系数
Figure BDA0003194352590000046
Figure BDA0003194352590000051
2.3)令k＝k+1，若k<Ne，则返回步骤2.2)；否则执行步骤3)。
进一步地，步骤3)具体为：
3.1)令k＝1；
3.2)根据第k个公共面F两边四面体网格单元的信息，根据下式更新第k个公共面左边四面体网格单元KL,F中与面F对应的矢量基函数系数
Figure BDA0003194352590000052
和右边四面体网格单元KR,F中与面F对应的矢量基函数系数
Figure BDA0003194352590000053
的值，更新后的矢量基函数系数分别为
Figure BDA0003194352590000054
和
Figure BDA0003194352590000055
Figure BDA0003194352590000056
Figure BDA0003194352590000057
其中，SK,L和SK,R分别表示第k个公共面F左右两边四面体网格单元的体积，
Figure BDA0003194352590000058
nF表示四面体网格单元KL,F在公共面F上的外法向单位矢量；
3.3)令k＝k+1，若k<Nf，则返回步骤3.2)；否则，执行步骤4)；Nf为计算区域Ω中公共面的总数。
进一步地，步骤4)具体为：
将步骤3)求得的更新后的矢量基函数系数带入下式计算得到平衡通量σh(x)：
Figure BDA0003194352590000059
其中，
Figure BDA0003194352590000061
本发明的有益效果是：
1.本发明给出了一种简单易行的平衡通量构造方法，该方法仅需要在Th中每个单元上求解一个三阶线性方程组和二阶线性方程组(d＝2)，或仅需要在Th中每个单元上求解一个四阶线性方程组和三阶线性方程组(d＝3)。
2.本发明构造的平衡通量是一个全局平衡通量，可用于有限元数值模拟中数值解的误差估计。
3.本发明的平衡通量构造方法，基于对计算区域Ω进行三角形剖分或四面体剖分的基础上实现，能够更好地逼近曲线或曲面边界。
附图说明
图1是区域Ω的三角形网格剖分、节点a相邻的网格单元构成节点a的单元片ωa示意图。
图2是三角形单元K的示意图。
图3是由图2构成的两个相邻网格单元的示意图。
图4是四面体单元K的示意图。
图5是由图4构成的两个相邻网格单元的示意图。
图6是本发明的方法流程图。
图7是热传导问题计算区域。
图8是不同网格数下每个时间步的后验误差与真实误差。
图9是计算金属圆柱散射入射平面电磁波问题的示意图。
图10是横磁波照射圆柱导体解析解。
图11是初始网格。
图12是数值解二维云图(网格数3835)。
图13是数值解二维云图(网格数8272)。
图14是数值解二维云图(网格数28795)。
图15是数值解二维云图(网格数104806)。
图16是能量误差估计及后验误差估计。
具体实施方式
首先，以泊松方程边值问题为例，给出局部平衡通量的概念(对于其他类型的模型方程，可以类似的给出局部平衡通量定义，详细内容参见背景技术中的文献[1]和[2])。考虑泊松方程边值问题，对于f∈L2(Ω)，求函u：Ω→R满足方程
-Δu＝f inΩ
Figure BDA0003194352590000071
其中，
Figure BDA0003194352590000072
(d＝2，3)是问题的定义域。
根据有限元理论，方程(1)的变分形式为：求
Figure BDA0003194352590000073
满足
Figure BDA0003194352590000074
其中，
Figure BDA0003194352590000075
(·，·)Ω表示在区域Ω上的内积。
设uh是方程(2)的有限元近似解。由近似解构造一个矢量函数σh，如果σh满足
Figure BDA0003194352590000076
Figure BDA0003194352590000077
称σh为局部平衡通量函数。其中，
Figure BDA0003194352590000078
Th表示对计算区域Ω进行三角形(d＝2)或四面体(d＝3)剖分得到的网格单元的集合。
以下参照图6，分别对计算区域Ω进行三角形(d＝2)剖分和四面体(d＝3)剖分时，本发明提出的平衡通量构造方法进行详述。
一、对计算区域Ω进行三角形剖分
计算区域Ω被剖分为Ne个三角形网格单元，Ω中包含Nf个公共边，本发明提出的平衡通量构造方法，步骤如下：
第一步：在每个三角形网格单元上定义矢量基函数；
在集合Th中第k个三角形网格单元K(三角形定义如图2所示)上，定义与三角形网格单元K的每个边对应的矢量基函数γki(x)：
Figure BDA0003194352590000081
其中，aki(i＝1,2,3)表示第k个三角形网格单元K三个顶点的坐标，lki(i＝1,2,3)表示第k个三角形网格单元K三边的长度，S表示第k个三角形网格单元K的面积，k分别取1，2，…，Ne。
第二步：在每个三角形网格单元上求解矢量基函数的系数；
求解如下方程组
Figure BDA0003194352590000082
求得与第k个三角形网格单元K的三个边lki(i＝1,2,3)分别对应的矢量基函数的系数
Figure BDA0003194352590000083
其中
Figure BDA00031943525900000814
f是公式(1)中等式的右端项，f可以是任意可积函数，S表示第k个三角形网格单元K的面积，k分别取1，2，…，Ne。
第三步：在相邻三角形网格单元的公共边的两侧更新矢量基函数的系数；
设集合Th中任意两个相邻三角形网格单元的公共边为F，如图3所示。记根据第一步在公共边F左边的三角形网格单元KL,F中求得的与边F对应的矢量基函数系数值为
Figure BDA0003194352590000085
在公共边F右边的三角形网格单元KR,F中求得的F对应的矢量基函数系数
Figure BDA0003194352590000086
值为
Figure BDA0003194352590000087
计算
Figure BDA0003194352590000088
Figure BDA0003194352590000089
用
Figure BDA00031943525900000810
和
Figure BDA00031943525900000811
分别更新
Figure BDA00031943525900000812
和
Figure BDA00031943525900000813
的值。
其中，SK，L和SK，R分别表示公共边F左、右两边三角形网格单元的面积
Figure BDA0003194352590000091
nF表示三角形网格单元KL，F在公共边F上的外法向单位矢量，
Figure BDA0003194352590000092
为数值解在三角形网格单元KL，F上的梯度；
Figure BDA0003194352590000093
为数值解在三角形网格单元KR，F上的梯度。当SK，L＝SK，R时，不更新系数。
对集合Th中所有公共边完成第三步操作。
第四步：计算平衡通量；
在每个三角形网格单元中，令
Figure BDA0003194352590000094
其中，
Figure BDA0003194352590000095
为数值解在第k个三角形网格单元K上的梯度；
γk1(x)，γk2(x)，γk3(x)分别为第k个三角形网格单元K的三条边对应的矢量基函数；
Figure BDA0003194352590000096
分别为第k个三角形网格单元K的三条边对应的更新后矢量基函数的系数；
则平衡迪量为
Figure BDA0003194352590000097
σh(x)为本发明构造的平衡通量，该平衡通量可用于背景技术中文献[2]中后验误差估计方法。
二、对计算区域Ω进行四面体剖分
计算区域Ω被剖分为Ne个四面体网格单元，Ω中包含Nf个公共面，本发明提出的平衡通量构造方法，步骤如下：
第一步：在每个四面体网格单元上定义矢量基函数；
在集合Th中第k个四面体网格单元K(四面体定义如图4所示)上，定义与四面体网格单元K的每个面对应的矢量基函数γki(x)：
Figure BDA0003194352590000098
其中，aki(i＝1，2，3，4)表示第k个四面体网格单元K四个顶点的坐标，lki(i＝1，2，3，4)表示第k个四面体网格单元K四个面的面积，S表示第k个四个面体网格单元K的体积。
第二步：在每个四面体网格单元上求解矢量基函数的系数；
求解如下方程组
Figure BDA0003194352590000101
求得与第k个四面体网格单元K的四个面lki(i＝1，2，3，4)对应的矢量基函数的系数
Figure BDA0003194352590000102
其中
Figure BDA0003194352590000103
f是公式(1)中等式的右端项，f可以是任意可积函数，S表示第k个四面体网格单元K的体积，k分别取1，2，...，Ne。
第三步：在相邻四面体网格单元公共面的两侧更新矢量基函数的系数；
设Th中任意两个相邻四面体网格单元的公共面为F，如5所示。记根据第一步在公共面F左边四面体网格单元KL，F中求得的与面F对应的矢量基函数系数值为
Figure BDA0003194352590000104
在公共面F右边四面体网格单元KR，F中求得的与面F对应的矢量基函数系数值为
Figure BDA0003194352590000105
计算
Figure BDA0003194352590000106
Figure BDA0003194352590000107
用
Figure BDA0003194352590000108
和
Figure BDA0003194352590000109
分别更新
Figure BDA00031943525900001010
和
Figure BDA00031943525900001011
的值。
其中，SK，L和SK，R分别表示公共面F左右两边四面体网格单元的体积
Figure BDA00031943525900001012
nF表示四面体网格单元KL，F在公共面F上的外法向单位矢量，
Figure BDA00031943525900001013
为数值解在四面体网格单元KL，F上的梯度；
Figure BDA00031943525900001014
为数值解在四面体网格单元KR，F上的梯度。当SK，L＝SK，R时，不更新系数。
对集合Th中所有公共面完成第三步操作。
第四步：计算平衡通量；
在每个四面体网格单元中，令
Figure BDA0003194352590000111
其中，
Figure BDA0003194352590000112
为数值解在第k个四面体网格单元K上的梯度；
γk1(x)，γk2(x)，γk3(x)，γk4(x)分别为第k个四面体网格单元K的四个面对应的矢量基函数；
Figure BDA0003194352590000113
分别为第k个四面体网格单元K的四个面对应的更新后矢量基函数的系数；
则平衡通量为
Figure BDA0003194352590000114
σh(x)为本发明构造的平衡通量，该平衡通量可用于背景技术中文献[2]中后验误差估计方法。
下面再通过具体应用实例对本发明作进一步说明。
实例1：
计算热传导问题的后验误差估计。材料中的传热过程，可以由傅里叶传热方程描述：
Figure BDA0003194352590000115
其中，T表示材料温度，
Figure BDA0003194352590000116
表示热扩散系数，与材料的性质有关，其中k＝48W/m.k表示材料热传导系数，ρ＝7850kg/m3表示材料密度，cp＝461j/K表示材料比热容。计算区域如图7所示，其中L＝0.1m，H＝0.05m。
设初边值条件为
Figure BDA0003194352590000121
Figure BDA0003194352590000122
Tx＝0＝T1＝400K
Tx＝0.1＝T2＝300K
T(x，y，0)＝T0＝300K
上述问题有解析解，记为T(x，y，t)
Figure BDA0003194352590000123
式中，x，y表示空间坐标；t表示时间变量；T1和T2分别表示图7中计算区域左右两边的初始温度；T0表示材料的初始温度。
对于上述问题，应用经典有限元方法(单元基函数取一阶线性单元)计算得到数值解，记为第n个时间步下的解为
Figure BDA0003194352590000124
表示第n个时间步对应的时刻下，计算区域的所有网格节点处的温度值。令每个时间步的数值解与解析解差的能量范数为η*，
Figure BDA0003194352590000125
称η*为“真实误差”。
基于数值解，应用本发明构造平衡通量后，可以计算得到数值解的后验误差，记为η。计算结果如图8所示。图8中给出了在三个不同网格数目下，从第1到第30时间步，每个时间步的真实误差和后验误差。可以看到，基于本发明构造平衡通量计算得到的后验误差可以用于对数值解精度的估计。并且随着网格数的增加，后验误差逐渐逼近真实误差，这表明后验误差的估计精度随着网格数的增加也相应地增加。上述结果表明本发明给出的平衡通量构造方法，能够满足背景技术中文献[1]中后验误差估计方法对平衡通量的要求，从而能够正确计算出后验误差。
实例2：
计算金属圆柱散射入射平面电磁波问题的后验误差估计，问题如图9所示。在入射电磁波是平面波的情况下，上述问题归结为时谐场问题.在二维情况下，时谐电磁场波动方程可以分解为电场和磁场满足的标量亥姆霍兹方程，方程如下：
Figure BDA0003194352590000131
Figure BDA0003194352590000132
入射平面电磁波
Figure BDA0003194352590000133
为
Figure BDA0003194352590000134
式中，Ez表示沿图9中负z方向电场强度的值；Hy表示沿图9中y方向磁场强度的值；∈r表示介电常数；μr表示磁导率；k0表示波数；
Figure BDA0003194352590000135
表示散度算子；
Figure BDA0003194352590000136
表示梯度算子；φ表示平面电磁波入射角；j表示虚数单位。
为了将计算空间控制在有限空间范围内，在金属圆柱外增加截断边界。采用一阶吸收边界设置该截断边界。将一阶吸收边界条件加载到圆形边界上，散射场Esca需要满足
Figure BDA0003194352590000137
计算参数：电磁波波长λ＝3m，周期T＝10-8s，金属圆柱半径a＝λ/2，吸收边界半径为3λ/2，波数
Figure BDA0003194352590000138
磁导率μ0＝4π×10-7H/m，真空介电常数ε0＝8.854187817×10-12F/m。
上述问题电场的解析解为
Figure BDA0003194352590000139
式中，j表示虚数单位；Jn()表示第一类汉克尔函数；
Figure BDA00031943525900001310
表示汉克尔函数；ρ表示图9中，圆外任意一点到圆心的距离；φ表示平面电磁波入射角。
图10是横磁波照射圆柱导体电场散射场解析解。
采用经典有限元方法(单元基函数取一阶线性单元)求解上述问题，图11是初始网格，图12～15给出了在不同网格剖分数下的数值解。可以看到，随着网格数的增加，数值解逐渐逼近图10中的解析解。
用本发明构造平衡通量，计算数值解的后验误差。图16给出了不同网格数下真实误差与后验误差结果，表明随着网格数的增加，后验误差与真实误差均不断减小，后验误差逐渐接近真实误差，但后验误差始终是真实误差的一个上界。上述结果表明本发明给出的平衡通量构造方法，能够满足背景技术中文献[1]中后验误差估计方法对平衡通量的要求，从而能够正确计算出后验误差。这也表明本发明给出的平衡通量构造方法是正确的和有效的。

Claims (9)

1.用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于，包括以下步骤：
1)在集合Th中每个网格单元上定义矢量基函数；Th表示对计算区域Ω进行剖分所得到的网格单元的集合；
2)在每个网格单元上计算所述矢量基函数的系数；
3)在相邻网格单元公共边/面的两侧更新矢量基函数的系数；
4)基于更新后的矢量基函数的系数和步骤1)定义的矢量基函数，计算平衡通量。
2.根据权利要求1所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
计算区域Ω被剖分为Ne个三角形网格单元，
步骤1)具体为：
1.1)令k＝1；
1.2)在第k个三角形网格单元上，定义和三边分别对应的矢量基函数γki(x)：
Figure FDA0003194352580000011
其中，aki表示第k个三角形网格单元K三个顶点的坐标，lki表示第k个三角形网格单元K三边的长度，S表示第k个三角形网格单元K的面积；
1.3)令k＝k+1，若k<Ne，则返回步骤1.2)；否则执行步骤2)。
3.根据权利要求2所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤2)具体为：
2.1)令k＝1；
2.2)根据第k个三角形网格单元K的信息，构造下述方程，求解方程得到在步骤1)中定义的第k个三角形网格单元K中与其三个边分别对应的矢量基函数的系数
Figure FDA0003194352580000021
Figure FDA0003194352580000022
2.3)令k＝k+1，若k<Ne，则返回步骤2.2)；否则执行步骤3)。
4.根据权利要求3所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤3)具体为：
3.1)令k＝1；
3.2)根据第k个公共边F两边三角形网格单元的信息，根据下式更新第k个公共边左边三角形网格单元KL,F中与边F对应的矢量基函数系数
Figure FDA0003194352580000023
和右边三角形网格单元KR,F中与边F对应的矢量基函数系数
Figure FDA0003194352580000024
的值，更新后的矢量基函数系数分别为
Figure FDA0003194352580000025
和
Figure FDA0003194352580000026
Figure FDA0003194352580000027
Figure FDA0003194352580000028
其中，SK,L和SK,R分别表示第k个公共边F左右两边三角形网格单元的面积，
Figure FDA0003194352580000029
nF表示三角形网格单元KL,F在公共边F上的外法向单位矢量；
3.3)令k＝k+1，若k<Nf，则返回步骤3.2)；否则，执行步骤4)；Nf为计算区域Ω中公共边的总数。
5.根据权利要求4所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤4)具体为：
将步骤3)求得的更新后的矢量基函数系数带入下式计算得到平衡通量σh(x)：
Figure FDA0003194352580000031
其中，
Figure FDA0003194352580000032
Figure FDA0003194352580000033
为数值解在三角形网格单元K上的梯度。
6.根据权利要求1所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
计算区域Ω被剖分为Ne个四面体网格单元，
步骤1)具体为：
1.1)令k＝1；
1.2)在第k个四面体网格单元上，定义和四个面分别对应的矢量函数
Figure FDA0003194352580000034
其中，aki表示第k个四面体网格单元K四个顶点的坐标，lki表示第k个四面体网格单元K四个面的面积，S表示第k个四个面网格单元K的体积；
1.3)令k＝k+1，若k<Ne，则返回步骤1.2)；否则执行步骤2)。
7.根据权利要求6所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤2)具体为：
2.1)令k＝1；
2.2)根据第k个四面体网格单元K的信息，构造下述方程，求解方程得到在步骤1)中定义的第k个四面体网格单元中与其四个面分别对应的矢量基函数的系数
Figure FDA0003194352580000035
Figure FDA0003194352580000036
2.3)令k＝k+1，若k<Ne，则返回步骤2.2)；否则执行步骤3)。
8.根据权利要求7所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤3)具体为：
3.1)令k＝1；
3.2)根据第k个公共面F两边四面体网格单元的信息，根据下式更新第k个公共面左边四面体网格单元KL,F中与面F对应的矢量基函数系数
Figure FDA0003194352580000041
和右边四面体网格单元KR,F中与面F对应的矢量基函数系数
Figure FDA0003194352580000042
的值，更新后的矢量基函数系数分别为
Figure FDA0003194352580000043
和
Figure FDA0003194352580000044
Figure FDA0003194352580000045
Figure FDA0003194352580000046
其中，SK,L和SK,R分别表示第k个公共面F左右两边四面体网格单元的体积，
Figure FDA0003194352580000047
nF表示四面体网格单元KL,F在公共面F上的外法向单位矢量；
3.3)令k＝k+1，若k<Nf，则返回步骤3.2)；否则，执行步骤4)；Nf为计算区域Ω中公共面的总数。
9.根据权利要8所述的用于有限元数值模拟后验误差估计的平衡通量构造方法，其特征在于：
步骤4)具体为：
将步骤3)求得的更新后的矢量基函数系数带入下式计算得到平衡通量σh(x)：
Figure FDA0003194352580000048
其中，
Figure FDA0003194352580000051
CN202110886391.1A 2021-08-03 2021-08-03 用于电场/温度有限元数值解的后验误差估计方法 Active CN113761762B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110886391.1A CN113761762B (zh) 2021-08-03 2021-08-03 用于电场/温度有限元数值解的后验误差估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110886391.1A CN113761762B (zh) 2021-08-03 2021-08-03 用于电场/温度有限元数值解的后验误差估计方法

Publications (2)

Publication Number Publication Date
CN113761762A true CN113761762A (zh) 2021-12-07
CN113761762B CN113761762B (zh) 2023-10-20

Family

ID=78788406

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110886391.1A Active CN113761762B (zh) 2021-08-03 2021-08-03 用于电场/温度有限元数值解的后验误差估计方法

Country Status (1)

Country Link
CN (1) CN113761762B (zh)

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FI20030550A0 (fi) * 2003-04-11 2003-04-11 Pajunen Petri Menetelmä numeerisen ratkaisun tarkkuuden arvioimiseksi
JP2009301149A (ja) * 2008-06-10 2009-12-24 Kyushu Univ 非線形の有限要素法による構造解析方法、及びプログラム、記録媒体、シミュレーション装置
CN102446357A (zh) * 2011-11-23 2012-05-09 浙江工商大学 基于自适应有限元的水平集sar图像分割方法
CN110059327A (zh) * 2018-11-28 2019-07-26 电子科技大学 一种基于辐射换热的三维有限元模拟方法
CN110807289A (zh) * 2020-01-07 2020-02-18 北京唯智佳辰科技发展有限责任公司 基于后验误差估计的集成电路自适应有限元网格细分方法
CN111222241A (zh) * 2020-01-06 2020-06-02 中国人民解放军国防科技大学 热化学非平衡条件下流场数据的数值计算方法和装置
EP3671489A1 (de) * 2018-12-20 2020-06-24 Siemens Aktiengesellschaft Verfahren zum bestimmen von bauvorschriften für ein additives fertigungsverfahren, verfahren zum erstellen einer datenbank mit korrekturmassnahmen für die prozessführung eines additiven fertigungsverfahrens und computer-programmprodukt
CN111638556A (zh) * 2020-06-09 2020-09-08 东华理工大学 基于地空分解策略的大地电磁正演方法及装置、存储介质
CN112131762A (zh) * 2020-08-07 2020-12-25 上海大学 模拟马氏体相变的网格自适应有限元方法
CN112581624A (zh) * 2020-12-23 2021-03-30 中国计量大学 一种基于距离函数定义边界的二维有限元网格剖分算法
US11030365B1 (en) * 2016-05-11 2021-06-08 Comsol Ab Systems and methods for determining finite elements in physics simulation systems for modeling physical systems using common geometry shape function spaces

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FI20030550A0 (fi) * 2003-04-11 2003-04-11 Pajunen Petri Menetelmä numeerisen ratkaisun tarkkuuden arvioimiseksi
JP2009301149A (ja) * 2008-06-10 2009-12-24 Kyushu Univ 非線形の有限要素法による構造解析方法、及びプログラム、記録媒体、シミュレーション装置
CN102446357A (zh) * 2011-11-23 2012-05-09 浙江工商大学 基于自适应有限元的水平集sar图像分割方法
US11030365B1 (en) * 2016-05-11 2021-06-08 Comsol Ab Systems and methods for determining finite elements in physics simulation systems for modeling physical systems using common geometry shape function spaces
CN110059327A (zh) * 2018-11-28 2019-07-26 电子科技大学 一种基于辐射换热的三维有限元模拟方法
EP3671489A1 (de) * 2018-12-20 2020-06-24 Siemens Aktiengesellschaft Verfahren zum bestimmen von bauvorschriften für ein additives fertigungsverfahren, verfahren zum erstellen einer datenbank mit korrekturmassnahmen für die prozessführung eines additiven fertigungsverfahrens und computer-programmprodukt
CN111222241A (zh) * 2020-01-06 2020-06-02 中国人民解放军国防科技大学 热化学非平衡条件下流场数据的数值计算方法和装置
CN110807289A (zh) * 2020-01-07 2020-02-18 北京唯智佳辰科技发展有限责任公司 基于后验误差估计的集成电路自适应有限元网格细分方法
CN111638556A (zh) * 2020-06-09 2020-09-08 东华理工大学 基于地空分解策略的大地电磁正演方法及装置、存储介质
CN112131762A (zh) * 2020-08-07 2020-12-25 上海大学 模拟马氏体相变的网格自适应有限元方法
CN112581624A (zh) * 2020-12-23 2021-03-30 中国计量大学 一种基于距离函数定义边界的二维有限元网格剖分算法

Also Published As

Publication number Publication date
CN113761762B (zh) 2023-10-20

Similar Documents

Publication Publication Date Title
Neumann et al. Quantifying the influence of microstructure on effective conductivity and permeability: virtual materials testing
Andrieux et al. Solving Cauchy problems by minimizing an energy-like functional
Rosolen et al. An adaptive meshfree method for phase-field models of biomembranes. Part I: Approximation with maximum-entropy basis functions
Bui et al. A moving Kriging interpolation‐based meshless method for numerical simulation of Kirchhoff plate problems
Kong et al. A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity problems on domains with complex geometry
CN109684740B (zh) 一种基于混合网格及时间步长的电磁学多尺度计算方法
Masud et al. B-splines and NURBS based finite element methods for Kohn–Sham equations
CN111079278B (zh) 三维时域杂交间断伽辽金方法外加电磁源项的处理方法
CN113158527A (zh) 一种基于隐式fvfd计算频域电磁场的方法
Bourantas et al. Strong-form approach to elasticity: Hybrid finite difference-meshless collocation method (FDMCM)
Ho-Nguyen-Tan et al. Level set-based topology optimization for compliance and stress minimization of shell structures using trimmed quadrilateral shell meshes
Bhardwaj et al. Numerical solution of time fractional tricomi-type equation by an rbf based meshless method
CN116629079A (zh) 混合有限元空间构造及求解线弹性力学问题的方法及装置
Sladek et al. Local integro-differential equations with domain elements for the numerical solution of partial differential equations with variable coefficients
Bazile et al. Variational Multiscale error estimator for anisotropic adaptive fluid mechanic simulations: application to convection–diffusion problems
CN103984869B (zh) 一种预测复合材料热弹性有效属性和局部场的方法
Harvey et al. 3D topology optimization of sandwich structures with anisotropic shells
Zeng et al. Energy-conserved splitting spectral methods for two dimensional Maxwell’s equations
Doškář et al. Microstructure-informed reduced modes synthesized with Wang tiles and the Generalized Finite Element Method
Luo A POD-based reduced-order TSCFE extrapolation iterative format for two-dimensional heat equations
CN113761762A (zh) 用于有限元数值模拟后验误差估计的平衡通量构造方法
Saucan et al. Combinatorial ricci curvature and laplacians for image processing
Kalidindi et al. Spectral representation of higher-order localization relationships for elastic behavior of polycrystalline cubic materials
Rashidinia et al. A collocation method for the solution of nonlinear one-dimensional parabolic equations
Hu et al. Radial basis collocation method and quasi‐Newton iteration for nonlinear elliptic problems

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant