CN120409160B

CN120409160B - 一种基于多场耦合与动态网格处理的仿真方法及系统

Info

Publication number: CN120409160B
Application number: CN202510929370.1A
Authority: CN
Inventors: 程经理; 刘骁; 曾磊; 桂业伟; 徐航
Original assignee: Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Current assignee: Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date: 2025-07-07
Filing date: 2025-07-07
Publication date: 2025-09-02
Anticipated expiration: 2045-07-07
Also published as: CN120409160A

Abstract

本发明公开了一种基于多场耦合与动态网格处理的仿真方法及系统，属于计算模拟领域，包括步骤：建立并求解耦合物理场模型：对求解域进行有限元空间离散，建立并求解描述至少包含温度场、位移和应力场的耦合偏微分方程组，获得这些物理场在当前状态下的数值解；再基于物理准则进行渐进式单元移除；最后执行网格数据结构更新与求解器同步。本发明可以实现金属材料在物理量达到临界条件时的渐进式移除，并在此过程中动态、鲁棒地处理网格变化。

Description

一种基于多场耦合与动态网格处理的仿真方法及系统

技术领域

本发明涉及计算模拟领域，更为具体的，涉及一种基于多场耦合与动态网格处理的仿真方法及系统。

背景技术

许多工程应用，特别是航空航天领域（如热防护系统、发动机部件等），金属材料常处于高温、高应力、高速气流冲刷等耦合物理环境中，经历复杂的性能退化乃至质量损失过程（如烧蚀、剥蚀）。准确模拟这些过程对于部件设计、寿命预测和安全性评估至关重要。

现有数值模拟方法面临挑战：1)如何有效处理因材料损失导致的计算域边界形状和拓扑结构的动态演化，传统固定网格难以适用；2)如何建立既反映物理机制又能保证计算稳定性的材料移除准则和网格处理策略，避免例如数值不稳定或产生不符合物理实际的悬浮“孤岛”单元。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于多场耦合与动态网格处理的仿真方法及系统，实现金属材料在物理量（如应力）达到临界条件时的渐进式移除，并在此过程中动态、鲁棒地处理网格变化。

本发明的目的是通过以下方案实现的：

一种基于多场耦合与动态网格处理的仿真方法，包括以下步骤：

建立并求解耦合物理场模型：对求解域进行有限元空间离散，建立并求解描述至少包含温度场、位移和应力场的耦合偏微分方程组，获得这些物理场在当前状态下的数值解；

再基于物理准则进行渐进式单元移除；

最后执行网格数据结构更新与求解器同步。

进一步地，所述耦合物理场模型的材料本构关系包括温度场对力学行为的耦合效应。

进一步地，所述基于物理准则进行渐进式单元移除，具体包括如下子步骤：

S11，状态检查与单元标记：基于获得的物理场的数值解，遍历预定检查范围内的活动单元，并根据至少一个预设的物理量阈值，将满足移除条件的单元标记为“待移除单元”；

S12，新暴露边界的识别与更新：对于每一个被标记为“待移除单元”E_del，在其被物理移除之前，检查其所有活动且未被标记为待移除的邻居单元N，若E_del与N之间的共享内部面F_internal当前不属于“暴露边界”，则将F_internal更新为新的“暴露边界”，确认为新的被侵蚀面；

S13，单元的物理移除：从有限元网格数据结构中删除所有当前被标记为“待移除单元”的单元；

S14，迭代剥蚀：重复执行步骤S11~步骤S13，检查范围聚焦于所有当前为“暴露边界”的活动单元，直至满足迭代终止条件。

进一步地，在基于物理准则进行渐进式单元移除与执行网格数据结构更新与求解器同步之间，还包括步骤：

S21，在执行完渐进式单元移除操作后，从预定义的“稳定边界”出发，通过图搜索算法识别所有与“稳定边界”连通的活动单元；

S22，将所有活动但未被识别为连通的单元判定为“孤岛单元”；

S23，从有限元网格数据结构中删除所有识别出的“孤岛单元”。

进一步地，所述执行网格数据结构更新与求解器同步，具体包括如下子步骤：

S31，更新连接关系与状态：在执行完单元移除后，修改网格数据结构，以反映单元移除导致的节点与单元之间连接关系的改变，以及单元和节点活动状态的更新；

S32，重建邻接与边界信息：基于更新后的连接关系，重新计算或更新活动单元间的邻接信息，并维护或重建构成当前所有已定义边界的几何实体及其拓扑信息；

S33，对并行计算同步并行环境数据：通过处理器间通信，同步因单元移除导致的各处理器所持有的网格部分之间共享边界的几何和拓扑信息，确保全局网格数据的一致性；

S34，通知主求解框架：通过内部机制通知主有限元求解框架，表明网格数据结构已发生变更，以触发求解器内部状态的重新初始化。

进一步地，求解所述耦合偏微分方程组具体采用牛顿法进行非线性迭代求解，其中每次牛顿迭代中产生的线性方程组采用预条件迭代法求解。

进一步地，在步骤S11中，所述状态检查与单元标记包括子步骤：

如果是基于应力的状态检查与单元标记，则所述检查范围为所有当前邻接于至少一个已定义的“暴露边界”的活动单元；还要计算每个被检查单元内的应力指标。

进一步地，步骤S14中，所述迭代终止条件，包括以下至少一项：

在一轮迭代的状态检查与单元标记步骤中，没有新的单元被标记为“待移除单元”；或

在迭代开始时，“待移除单元集合”即为空；或

达到预设的最大迭代次数。

进一步地，步骤S21中，所述图搜索算法包括广度优先搜索 BFS、深度优先搜索DFS，或者连通性分析方法。

一种基于多场耦合与动态网格处理的仿真系统，包括计算机设备，所述计算机设备包括处理器和存储器，所述存储器中存储有计算机程序，当所述计算机程序被所述处理器加载时执行如上任一项所述的方法。

本发明的有益效果包括：

（1）模拟渐进过程：特别是基于应力的迭代移除机制，能够模拟材料逐层损失的物理过程，而非一次性移除。

（2）提升鲁棒性：孤岛单元检测与移除机制，显著提高了大规模单元移除仿真的数值稳定性和结果的物理合理性。

（3）灵活性与适用性：可适用于基于阈值的材料失效准则（温度、应力等），并可为更复杂的航空航天烧蚀/侵蚀模型提供动态网格处理的基础框架。对于烧蚀/侵蚀模型，需要进一步扩展边界的影响，考虑化学反应带来的烧蚀影响以及颗粒冲击带来的剥蚀影响。

（4）耦合效应考虑：通过在耦合求解框架内执行，保证了物理场演化与几何形状变化之间的相互影响得以考虑。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1a为基于温度的网格删除策略；

图1b为基于应力的网格删除策略；

图2a为基于温度的网格删除示意图；

图2b为基于应力的网格删除流程图；

图3为基于温度的网格删除的时变计算结果；

图4a为基于应力的网格删除的时变计算结果中的温度分布图；

图4b为基于应力的网格删除的时变计算结果中的应力分布图；其中，vonmises_stress表示应力。

具体实施方式

本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。

本发明旨在提供一种能够在耦合多物理场有限元模拟框架下，有效处理金属材料因超过物理阈值（如应力）而发生的渐进式移除，并在此过程中动态更新计算域边界、维持数值稳定性的数值方法，主要解决在模拟材料逐层损失时，如何自动识别并标记失效单元、如何根据移除过程更新暴露边界、如何避免因单元移除产生数值孤岛，以及如何将这些过程与主物理场求解器有效集成的问题。

具体而言，在一个耦合多物理场有限元计算框架内，提供一种基于多场耦合与动态网格处理的仿真方法，包括以下步骤：

步骤(a)，建立并求解耦合物理场模型。

对求解域进行有限元空间离散，建立并求解描述至少包含温度场和位移/应力场的耦合偏微分方程组，获得这些物理场在当前状态下的数值解；其中，模型的材料本构关系包含温度场对力学行为的耦合效应，例如通过热应变影响总应变和应力计算。

具体而言，采用有限元法 (Finite Element Method, FEM) 对求解域进行空间离散。这为后续基于单元的物理量计算和单元移除操作提供了基础。针对需要模拟物理计算域的物理现象，建立其数学描述，即偏微分方程组 (PDEs)，主要包括：

能量守恒(非稳态热传导)：描述温度场T的时空演化。其强形式（不考虑体积热源）为：

相应的弱形式为：对所有合适的测试函数，求解T使得：

其中，是材料密度，c_p是材料比热容，k是材料热导率，q_n是通过边界Γ_q的法向热流密度，n是边界Γ_q的单位外法线向量，V是微元体积，S是微元面积，是整个计算域。公式(2)从左到右分别表示瞬态项、扩散项以及边界热流。

动量守恒(准静态力平衡方程)：描述材料在载荷和温度变化下的位移场u和应力场的平衡。由于未考虑惯性项，该方程简化为力的平衡方程，其强形式为：

相应的弱形式为：对所有合适的测试函数，求解u使得：

其中，σ是柯西应力张量，是体积力，是施加在边界Γ_t上的面力。

材料本构与耦合：连接不同物理场的桥梁。

热应变：结构在外部热流以及压力作用下，材料产生变形，例如，线性热膨胀，其中α是热膨胀系数，Tref是参考温度，I是二阶单位张量。

应变分解：总应变由位移梯度计算得到，且被分解为机械应变和热应变，即。

弹性应力：机械应变通过弹性张量C产生应力：。

耦合关系：通过以上关系，温度场T(x, t)的变化会引起热应变，进而影响机械应变和应力场σ，构成了热对力的驱动。

数值求解策略：

时间离散：对于瞬态问题，采用时间积分方法（例如，隐式的后向欧拉法、Crank-Nicolson 法等）将时间导数离散化。

非线性求解：耦合的方程组系统在离散后通常形成一个大型的非线性代数方程组R(U,t)=0，其中U是包含所有节点未知自由度（如温度、位移分量）的全局向量。使用牛顿法进行迭代求解：

在第k次迭代，计算残差向量R(Uk)。

计算雅可比矩阵（切线刚度矩阵）。

求解线性系统得到解增量。

更新解，其中λ是线搜索步长。

重复直至残差满足收敛判据。

线性求解：上述的线性系统通常规模巨大且稀疏，采用预条件迭代法求解，如广义最小残差法(GMRES) 或双共轭梯度稳定法(BiCGStab)，并配合高效的预条件器，例如代数多重网格 (AMG) 或不完全 LU 分解。

步骤(b)，实现基于物理准则的渐进式单元移除。

1）状态检查与单元标记：基于步骤(a) 获得的物理场解，遍历预定检查范围内的活动单元，并根据至少一个预设的物理量阈值（例如，单元应力指标超过应力阈值）将满足移除条件的单元标记为“待移除单元”；

2）新暴露边界的自动识别与更新：对于每一个被标记为“待移除单元”E_del，在其被物理移除之前，检查其所有活动且未被标记为待移除的邻居单元N，若E_del与N之间的共享内部面F_internal当前不属于“暴露边界”，则将F_internal更新为新的“暴露边界”，确认为新的被侵蚀面；

3）单元的物理移除：从有限元网格数据结构中删除所有当前被标记为“待移除单元”的单元；

4）迭代剥蚀策略：重复执行步骤 1）中的基于应力流程、步骤2) 和步骤3)，检查范围聚焦于所有当前为“暴露边界”的活动单元，直至满足迭代终止条件。

具体而言，在求解过程的特定阶段（例如，每个时间步结束后），执行以下基于物理准则的单元移除流程，这里给出了两种方式：一种是基于温度阈值的网格删除方式，另一种是基于应力阈值的网格删除方式。基于温度的网格删除方式用于测试算法的功能实用性，不具备应用实际环境；而基于应力阈值的网格删除功能，可用于实际金属材料的渐进损失模拟过程。

1）状态检查与标记（图1a中的步骤B）：

基于温度：遍历所有受侵蚀边界的单元温度。计算这些单元的温度指标（如单元平均温度，积分点最大温度等等）。若温度指标超过预设的温度阈值，则将该单元标记为“待移除”。如图1a所示。

基于应力：遍历所有受侵蚀边界的单元应力。计算这些单元的应力指标（如单元平均应力，积分点的最大 von Mises 应力）。若应力指标超过预设的应力阈值，则将该单元标记为“待移除”。如图1b所示。

2）边界更新（图1a和图1b的步骤C）：对于每一个被标记为“待移除”的单元 E_del：遍历其所有邻居单元 N。如果邻居 N 是活动的，且未被标记为“待移除”，则将单元 E_del 与邻居 N 之间的内部共享面 F_internal，添加到“暴露边界”的定义中，由此确认新的被侵蚀面。此步骤在物理移除单元之前完成，确保模拟材料移除后新表面的暴露。

3）单元物理移除（图1a和图1b的步骤D）：从网格数据结构中，删除所有被标记为“待移除”的单元。

4）迭代剥蚀(针对应力准则，（图1b中的步骤E、步骤C、步骤D）)：重复步骤 1）中的基于应力流程至-步骤3）。即，在新暴露的边界上以及原暴露边界上对应的未移除单元中，再次检查应力是否超限，并进行标记、边界更新、移除。此迭代过程模拟材料的逐层剥蚀，直至在一个迭代循环中没有新的单元被标记，或达到预设的最大迭代次数。对于温度准则，则没有此项限制，当受侵蚀面温度高于阈值，则被直接全部删除（如图2a所示），而应力准则由于材料受热/力载荷作用，材料内部也可能出现较大应力，材料的删除条件往往内部单元同样满足，而材料的损失往往是从被侵蚀面发展的。因此，限制了检索方式为逐层检查，递归删除网格（如图2b所示）。

步骤(c)，孤岛单元检测与处理（可选择）。

在执行完步骤 (b) 的移除操作后，从预定义的“稳定边界”出发，通过图搜索算法识别所有与“稳定边界”连通的活动单元；

将所有活动但未被识别为连通的单元判定为“孤岛单元”；

从有限元网格数据结构中删除所有识别出的“孤岛单元”。

具体而言，建议在应力准则下增加此项的连通性分析，（图1b中的步骤F）：在确认没有单元超限被标记，或者达到预设的最大迭代次数之后，开始进行“孤岛”单元检测分析，主要是进行连通性分析，包括：

定义一个或多个“稳定边界”标识符，实际过程中代表固定或不会被完全移除的区域边界。

从所有邻接于“稳定边界”的活动单元出发。

采用图搜索算法（如广度优先搜索 BFS 或深度优先搜索 DFS）遍历所有通过共享面或节点相互连接的活动单元，标记所有可达单元。

识别网格中所有处于活动状态但未被标记为可达的单元，这些即为“孤岛”单元。

从有限元网格数据结构中移除所有识别出的“孤岛”单元。

步骤(d)，网格数据结构更新与求解器同步。具体而言，执行如下步骤：

1）更新连接关系与状态：在执行完步骤 (b) 和/或 (c) 中的单元移除后，修改所述网格数据结构，以反映单元移除导致的节点与单元之间连接关系的改变，以及单元和节点活动状态的更新；

2）重建邻接与边界信息：基于更新后的连接关系，重新计算或更新活动单元间的邻接信息，并维护或重建构成当前所有已定义边界（包括新生成的“暴露边界”）的几何实体及其拓扑信息；

3）(若为并行计算) 同步并行环境数据：通过处理器间通信，同步因单元移除导致的各处理器所持有的网格部分之间共享边界的几何和拓扑信息，确保全局网格数据的一致性；

4）通知主求解框架：通过内部机制通知主有限元求解框架，表明网格数据结构已发生变更，以触发必要的求解器内部状态（如自由度映射、稀疏矩阵模式）的重新初始化。

为了更清楚地说明本发明的技术方案，下面结合附图对本发明进行详细描述。图1a~图1b展示了本发明方法中两种基于不同物理准则的网格单元移除策略的流程图。其中，图1a示意了基于温度的网格删除策略。在一个时间步内，该策略主要执行一次性的检查：在耦合场计算后（时间步开始），检查初始暴露边界上的单元温度（步骤 B），若超过阈值，则更新相邻边界（步骤 C）并删除标记单元（步骤 D），然后进入时间步结束阶段进行网格更新。这代表了一种根据当前温度场快照进行单元移除的模式。这种网格删除策略较为浅显，实际工作环境中，因为温度导致材料的删除是不符合物理规律。在此基础上，开展模拟基于应力导致的材料删除策略，模拟真实的热力耦合环境下，金属材料的渐进损伤过程。图1b示意了基于应力的网格删除策略。该策略体现了更复杂的迭代式、逐层移除逻辑：在耦合场计算后检查初始暴露边界应力并标记（步骤 B）；若有单元被标记，则更新相邻边界（步骤 C），然后删除标记单元（步骤 D）；关键在于删除操作后，再次检查当前所有暴露边界（包括新暴露的边界）上的单元应力（步骤 E），若发现新的超限单元，则返回步骤 C，形成一个步骤C、步骤D、步骤E、步骤C的迭代循环，模拟材料的逐层剥蚀，直到没有新的超限单元被发现；在整个迭代移除过程结束后（无论是否有单元被移除），执行孤岛检测与删除（步骤 F），以确保网格的物理连通性；最后进入时间步结束阶段进行最终网格更新。该策略的核心在于通过迭代检查新暴露面，模拟了由应力驱动的材料渐进式物理剥蚀过程。

图2a~图2b示意性展示了基于温度的网格删除策略（如图1a流程）可能产生的结果。根据热流分布（这里展示了非均匀的边界影响），温度超限的区域可能一次性导致单层或多层单元被移除，反映的是达到热阈值的瞬时状态，而非逐层推进。对应的计算结果如图3所示（对应顶部是热流边界，左右两侧是绝热边界，底部恒温），这里给出的初始热流边界是基于坐标的函数，当新暴露的边界出现，依据新边界的坐标进行生成新的热流边界用于后续计算，同样的方式可扩展至压力边界。后续对于航空航天领域的内外流场向结构内部传递的初始时刻热流分布进行参数化，新生成的边界依据坐标进行参数化热流边界更新。依据坐标来实现边界的更新能够较符合新的边界面临的实际环境，从而模拟材料的损失后退。

图3示意性展示了基于应力的网格删除策略（如图1b流程）的核心逻辑和效果。它强调了在一个时间步内增加的迭代逻辑：通过逐层检查（对应图1b的 C ->D ->E 循环）、删除单元、进行网格拓扑更新、通过邻域搜索确定新暴露边界、并再次检查这些新边界上的单元是否超限，从而实现材料从初始暴露面向内部渐进式移除的效果，更贴近物理剥蚀过程。对应的计算结果如图4a~图4b所示（对应顶部是恒定热流及压力边界，左右两侧是绝热边界，底部300K恒温以及固支边界）。对于此种网格删除策略，主要是考虑到材料损失由被侵蚀面开始，且应力超限不仅仅存在于受侵蚀面，需要从顶部开始逐层检查以及递层删除。

需要说明的是，图2a~图2b和图3是两种不同的计算边界下进行的数值模拟过程。

需要说明的是，本发明其他实施例场景中，可在通用的有限元分析软件平台基础上，通过开发用户自定义模块（例如 User Objects 或子程序）来实现，例如包括：

数据结构：使用集合存储待移除、已移除、新暴露边界的单元 ID。

网格接口：需要访问和修改有限元网格库提供的接口，以获取单元邻接关系、边界信息，并执行单元删除、边界添加/删除操作。

模块集成：将单元移除和孤岛检测模块嵌入到主求解器的流程控制中，例如在每个时间步求解完成之后、更新网格之前调用。

信息传递：确保物理场求解器能访问到最新的网格信息，并且单元移除模块能获取到当前时间步计算得到的物理场变量值（温度、应力等）。

连通性算法：图搜索算法（BFS 或 DFS）用于孤岛检测。

以上实施方式，均在本发明权利要求书保护范围之内。

描述于本发明实施例中所涉及到的单元可以通过软件的方式实现，也可以通过硬件的方式来实现，所描述的单元也可以设置在处理器中。其中，这些单元的名称在某种情况下并不构成对该单元本身的限定。

根据本发明实施例的一个方面，提供了一种计算机程序产品或计算机程序，该计算机程序产品或计算机程序包括计算机指令，该计算机指令存储在计算机可读存储介质中。计算机设备的处理器从计算机可读存储介质读取该计算机指令，处理器执行该计算机指令，使得该计算机设备执行上述各种可选实现方式中提供的方法。

作为另一方面，本发明实施例还提供了一种计算机可读介质，该计算机可读介质可以是上述实施例中描述的电子设备中所包含的；也可以是单独存在，而未装配入该电子设备中。上述计算机可读介质承载有一个或者多个程序，当上述一个或者多个程序被一个该电子设备执行时，使得该电子设备实现上述实施例中所述的方法。

Claims

1.基于多场耦合与动态网格处理的仿真方法，其特征在于，包括以下步骤：

再基于物理准则进行渐进式单元移除；

最后执行网格数据结构更新与求解器同步；

所述基于物理准则进行渐进式单元移除，具体包括如下子步骤：

S14，迭代剥蚀：重复执行步骤S11~步骤S13，检查范围聚焦于所有当前为“暴露边界”的活动单元，直至满足迭代终止条件；

在基于物理准则进行渐进式单元移除与执行网格数据结构更新与求解器同步之间，还包括步骤：

S23，从有限元网格数据结构中删除所有识别出的“孤岛单元”；

在步骤S11中，所述状态检查与单元标记包括子步骤：

2.根据权利要求1所述的基于多场耦合与动态网格处理的仿真方法，其特征在于，所述耦合物理场模型的材料本构关系包括温度场对力学行为的耦合效应。

3.根据权利要求1~2任一项所述的基于多场耦合与动态网格处理的仿真方法，其特征在于，所述执行网格数据结构更新与求解器同步，具体包括如下子步骤：

4.根据权利要求1所述的基于多场耦合与动态网格处理的仿真方法，其特征在于，求解所述耦合偏微分方程组具体采用牛顿法进行非线性迭代求解，其中每次牛顿迭代中产生的线性方程组采用预条件迭代法求解。

5.根据权利要求1所述的基于多场耦合与动态网格处理的仿真方法，其特征在于，步骤S14中，所述迭代终止条件，包括以下至少一项：

在迭代开始时，“待移除单元集合”即为空；或

达到预设的最大迭代次数。

6.根据权利要求1所述的基于多场耦合与动态网格处理的仿真方法，其特征在于，步骤S21中，所述图搜索算法包括广度优先搜索 BFS、深度优先搜索 DFS，或者连通性分析方法。

7.一种基于多场耦合与动态网格处理的仿真系统，其特征在于，包括计算机设备，所述计算机设备包括处理器和存储器，所述存储器中存储有计算机程序，当所述计算机程序被所述处理器加载时执行如权利要求1~2任一项所述的方法。