在ARCHER2平台上，使用openFOAM对多亿体素的多孔材料微观CT图像进行可扩展的CFD（计算流体动力学）模拟

《Journal of Computational Science》：Scalable CFD simulations in multi-billion voxel micro-CT images of porous materials using openFOAM on ARCHER2

　　朱利安·梅斯|加文·J·普林格尔|汉娜·P·门克
赫瑞瓦特大学，英国

**摘要**
本研究探讨了使用高性能计算（HPC）通过计算流体动力学（CFD）来模拟多亿体素微CT图像中多孔材料的流动和传输过程。研究分析了两种不同的岩石样本——本特海默砂岩和埃斯塔伊亚德碳酸盐岩，这两种岩石分别代表了两种不同的地质构造。本特海默砂岩图像的尺寸为1,950 × 1,950 × 10,800体素，分辨率为6微米，共包含410亿个体素，具有相对均匀的结构；而埃斯塔伊亚德碳酸盐岩图像的尺寸为1,144 × 1,144 × 6,000体素，分辨率为3.9676微米，包含80亿个体素，表现出更大的异质性，其中包括微孔区域。这些图像通过我们基于OpenFOAM的数值求解器GeoChemFoam进行直接CFD模拟，利用了英国超级计算机ARCHER2的计算资源。本研究的一个关键方面是采用了Darcy-Brinkman-Stokes方法，该方法使用体积指示函数来表示固体表面，而不是复杂的网格。这使得可以高效且可扩展地生成结构化的、与体素对齐的笛卡尔网格。研究通过子体积分解技术评估了强 scaling 和弱 scaling 情况，证明了由于ARCHER2强大的计算能力，无需降低图像尺寸即可完成全分辨率CFD模拟。这项工作展示了HPC在处理大型高分辨率微CT数据时进行详细全尺度模拟的潜力。该方法依赖于一种网格策略，该策略利用基于体积指示函数的简单、可并行化的笛卡尔网格，无需复杂的表面贴合网格，从而能够实现对地质和工程应用中流动和传输过程的可扩展模拟。

**1. 引言**
在多孔材料中模拟流体流动和传输过程对于理解水文地质学[1]、化学工程[2]和材料科学[3]等领域的现象至关重要。成像技术的进步，特别是X射线微计算机断层扫描（micro-CT），使得能够获取高分辨率的3D数据集，揭示岩石和其他多孔材料的复杂孔隙结构[4][5]。然而，这些数据集日益增长的分辨率和体积带来了重大的计算挑战。

传统上，多孔材料中的流动和传输模拟是通过孔隙网络建模（PNM）来实现的，该方法将复杂的孔隙空间简化为一个由通道和节点组成的互连网络[6]。尽管PNM在计算效率上具有优势，并能捕捉关键的宏观特性，但它常常依赖于某些假设，这些假设限制了其在孔隙尺度上表现流动和传输细节的能力。相比之下，直接数值模拟（DNS）通过在孔隙几何结构上直接求解控制流体流动方程（无论是基于网格的Stokes/Navier–Stokes方程离散化方法[7]还是介观动力学方法如格子玻尔兹曼方法（LBM）[8]），由于计算基础设施的进步而越来越受到青睐。如今，大规模并行超级计算机的可用性、处理器速度和内存带宽的提高，以及像OpenFOAM这样的可扩展开源CFD框架，使得以前计算上不可行的高保真度大规模模拟成为可能。因此，DNS为研究复杂多孔介质中的流动和传输现象提供了一种强大而准确的方法。

尽管有这些优点，DNS仍面临重大的计算挑战。尽管微CT成像现在可以生成每个方向数千个体素的数据集，但由于内存和处理能力的限制，DNS通常仅限于每个维度几百个体素的较小子体积[10]。对于更大规模的模拟，研究人员经常返回使用PNM来克服这些限制。解决这一难题需要可扩展的计算框架和对高性能计算资源的访问。

将CFD应用于微CT图像的一个重大挑战在于计算网格的生成[11]。传统的CFD方法需要对孔隙空间进行网格划分，并在孔隙表面上施加边界条件来求解守恒方程。然而，在传统的孔隙尺度CFD工作流程中，从高分辨率微CT图像构建符合边界条件的计算域通常涉及表面重建（例如STL提取）、孔隙空间提取（单元移除）和细化/裁剪操作，这些操作在计算上非常复杂且耗费大量内存，在大问题规模下可能成为瓶颈[12]。对于非常大的（数十亿个体素）数据集，由于并行生成复杂网格的效率不高，这一问题变得更加突出。最近，王等人[13]使用LBM方法对包含超过300亿个单元的图像中的质子交换膜燃料电池进行了大规模、物理准确的建模，使用了超过20,000个CPU核心。然而，LBM使用的规则格子无法准确捕捉固体表面的复杂性[14]，并且LBM在扩展到反应传输建模时也存在一定困难[15]。

为了解决这些问题，浸没边界方法（IBM）作为一种有前景的替代方案应运而生[16]。IBM通过允许在结构化的笛卡尔网格上进行模拟，消除了对体适应网格的需求。已经开发出多种IBM技术，如连续驱动[17]、使用水平集的离散驱动[18]、幽灵流体[19]或裁剪单元[20]方法，以及基于粘度[21]或Darcy-Brinkman-Stokes（DBS）公式的惩罚方法[22]来模拟障碍物的存在。关于这些方法的全面回顾，请参考Mittal和Seo的工作[23]。

最近，Trebotich和Graves[20]使用带有裁剪单元方法的IBM来模拟大型二元多孔材料图像中的流体流动。他们展示了良好的可扩展性，使用大约80亿个单元的合成图像在130,000个CPU核心上进行了模拟，每处理器约30,000个单元时的可扩展性较弱。他们的方法也被应用于真实的微CT图像，包括贝德福德石灰岩（2048 × 2048 × 320体素）的CT扫描和页岩（1920 × 1600 × 640体素）的FIB-SEM扫描。尽管结果令人满意，但这些图像的规模仍小于典型的CT扫描，且模拟计算成本较高，每个时间步长需要20到40秒。这种高计算成本可能是由于通过裁剪单元方法显式实现了无滑移边界条件所致。

DBS方法为孔隙尺度建模提供了一个实用的微连续介质框架，因为它避免了在尖锐重建的流体-固体界面上施加显式无滑移条件的需要。这使得该方法特别适合于多孔结构仅部分解析或随时间演变的情况。除了二元孔隙/固体图像外，DBS还可以通过整合局部孔隙度和渗透率值来表示非二元数据集和欠解析区域，从而实现对亚体素孔隙度和双重孔隙行为的微连续介质描述[24]。DBS方法已成功应用于多孔材料的微CT图像中的流动和热传递[25]、碳酸盐岩中的矿物溶解和沉淀[26][27][28]，以及多尺度图像中的流动[29][30][31]。最近的研究进一步展示了DBS在生物膜生长和营养运输[32]以及生物矿化和MICP相关过程[33][34]中的应用，突显了微连续介质公式在广泛的应用领域中的独特建模能力。

在中等问题规模下，DBS模拟每个时间步长的成本可能高于仅基于孔隙无滑移条件的DNS，因为线性系统是在包含固体和界面单元的完整体素域上求解的。其主要优势在于工作流程避免了显式的界面重建和符合边界条件的孔隙空间网格划分（表面重建、孔隙空间提取、细化/裁剪和网格变形）。这些步骤破坏了结构化的网格索引，并引入了非局部/全局的网格操作，随着图像规模的增加而变得越来越昂贵和耗内存。通过保留一个结构化的、与体素索引（i,j,k）和网格实体（单元/面/点）之间存在确定性映射的笛卡尔网格，我们可以在分布式内存中直接生成完整的并行网格并初始化相应的场，从而实现了使用现代HPC系统对数十亿个体素微CT图像的可扩展CFD模拟，如下面的扩展研究和完整8-410亿个体素模拟所示。

我们的工作流程已集成到基于OpenFOAM的求解器GeoChemFoam中（可在https://github.com/GeoChemFoam获取），并通过模块加载（gcfoam）在英国超级计算机ARCHER2上使用。我们将该工作流程应用于两个微CT图像（本特海默砂岩和埃斯塔伊亚德碳酸盐岩），分别包含410亿和80亿个体素。我们的目标是在ARCHER2上使用多达1000个节点对完整图像进行流动和传输模拟，这代表了我们在ARCHER2上的资源分配限制。我们展示了强 scaling 和弱 scaling 结果，并证明了在不到一分钟的时间内为整个图像创建计算网格、在不到一小时内完成流动模拟以及在不到10小时内完成完整传输模拟的能力。与标准表面贴合网格工作流程不同，后者由于内存和预处理要求过高而无法在此规模上应用，我们的方法能够实现对数十亿个体素微CT数据集的高效、高保真度建模。

**2. 方法**

**2.1. 控制方程**
整个域Ω由流体域Ωf和固体域Ωs组成。不可压缩牛顿流体Ωf的稳态层流运动遵循Navier–Stokes方程：
(1) ??u = 0,
(2) ??u?u = ??p + ν?2u，
在流体-固体界面Γ处施加无滑移边界条件：
(3) u = 0 at Γ，
其中u（m/s）是速度，p（m2/s2）是运动压力，ν（m2/s）是运动粘度。我们的模拟采用恒定的运动压降ΔP，即在入口和出口处压力固定且速度梯度为零，固体壁面上也施加无滑移边界条件。样本的渗透率定义为：
(4) K = νUDLΔP，
其中UD是达西速度，L是样本长度。达西速度定义为：
(5) UD = QA，
其中Q是穿过样本的体积流量，A是样本的横截面积。如果雷诺数（Re）较小，即表征粘性力和惯性相对影响的参数较小（即Re = UDLpore?ν ? 1），则渗透率与压降无关。在我们的模拟中，我们使用水的粘度（ν = 10?? m2/s），并调整压降使得Re = 0.001，在这种情况下Navier–Stokes方程可以简化为Stokes方程：
(7) ??u = 0,
(8) 0 = ??p + ν?2u。

对于在Ωf中传输的物种，其守恒可以用对流-扩散方程表示：
(9) ?C?t + ??Cu ? ??D?C = 0，
其中C是物种的无量纲浓度，D（m2/s）是物种的分子扩散系数。Péclet数（Pe）用于表征对流和扩散的相对重要性：
(10) Pe = UDLpore?D。

在我们的模拟中，入口处的浓度恒定为1.0，出口和固体壁面上的浓度梯度为零。扩散系数设为D = 10?? m2/s，适用于水中的化学物种。通过调整压降使得Re = 0.001，此时Navier–Stokes方程简化为Stokes方程：
(7) ??u = 0,
(8) 0 = ??p + ν?2u。

**2.2. Darcy-Brinkman-Stokes公式**
为了捕捉复杂固体几何结构的影响，我们采用了体积平均方法。流体-固体界面用每个控制体积中的流体体积分数??f来描述。质量和动量守恒根据体积平均属性来解决：
(11) u? = 1V∫VfudV,
(12) p?f = 1Vf∫VfpdV,
(13) C?f = 1Vf∫VfCdV，
其中Vf和V分别是控制体积中的流体体积和总体积。注意，速度是相对于V（即达西速度）进行缩放的，而压力和浓度是相对于Vf进行缩放的。体积平均速度满足整个域Ω上的Darcy-Stokes-Brinkman方程：
(14) ??u? = 0,
(15) ?0 = ??p?f + ??ν?f?u? ? νKf?1u?，
其中Kf是控制体积的渗透率。νKf?1u?表示流体相与固体相之间的动量交换，即达西阻力。这一项在固体和微孔相中占主导地位，在流体相中则消失。需要注意的是，为了避免除以零的情况，假设在固体相中??f等于一个较小的值（例如0.0001）。

体积平均浓度满足体积平均对流-扩散方程：
(16) ????fC?f?t + ??C?fu? ? ???fD??C?f = 0，
其中D?是体素内的扩散张量。在完全包含在流体域内的控制体积中，扩散张量简单地等于分子扩散系数。然而，在包含固体或欠解析孔隙的单元中，扩散张量是速度场（流体动力扩散）和未解析结构的迂曲度的函数。为了简化，我们假设：
(17) ?D? = ?fD。

这种方法用于模拟多孔介质中的流动和传输已经得到了广泛的验证[24][26][29][37][38]。这里使用了一种简化的闭合方法来提供一致的可扩展性基准；更先进的弥散/迂曲模型可以直接整合进来[31]。下载：下载高清晰度图像（534KB）下载：下载全尺寸图像图1. Bentheimer的微CT图像，以及已发布的四标签分割结果（来自[39]）。黑色表示Viton套管，深灰色表示已解析的固体，浅灰色表示未解析的小类（在本研究中被视为固体），白色表示已解析的孔隙空间。固体表面经过平滑处理，以去除分割引起的人为粗糙度。

2.3. 样本1：Bentheimer砂岩
Bentheimer是一种浅海砂岩，主要由石英组成，含量约为95%，含有少量的长石和粘土，具有良好的颗粒大小分布，平均孔径为40μm[40]。微CT图像来自数字岩石门户[39]。该图像包含1950 × 1950 × 10,800个体素，分辨率为6μm，由Zeiss Versa 510 X射线显微镜拍摄。存储库[39]提供的分割图像包含四标签分割结果（Viton套管、已解析的固体、未解析区域和已解析的孔隙空间），这是通过手动灰度阈值划分得到的（三个阈值定义了四个类别；详见[41]）。

经过分割后，从图像中得到的孔隙率为0.18，使用Ruska气体渗透仪[40]测量的渗透率为1.64 D。尽管Bentheimer是一种相对均匀的介质，可以在较小的范围内分析其几何特性（如孔隙率和孔径分布），但如此大的范围对于表征流动特性（如绝对渗透率和相对渗透率）的代表性基本体积（REV）非常有用[41]。

在我们的研究中，我们假设未解析的部分为固体（即??f=0.0001）。这类占体素的约1%，主要是分割伪影，而不是目标微孔人口。然后使用拉格朗日平滑器对固体表面进行平滑处理，以消除分割引起的人为粗糙度（图1）。这种平滑步骤会产生一层具有部分体积孔隙率的界面单元（?0.000197%）组成，含有少量石英。它显示出双峰的孔径分布，颗粒内部有微孔隙。宏观孔隙的平均半径约为10μm，而微孔隙的平均半径约为0.4μm[37]。微CT图像来自BGS存储库[43]。该图像包含1202 × 1236 × 6000个体素，分辨率为3.9676μm，由Zeiss XRM-510微CT拍摄。使用Fiji中的Weka3D机器学习分割算法[44]对宏观孔隙进行分割。然后使用掺杂盐水的岩石图像来识别固体颗粒和未连接的微孔隙。接着对孔隙空间、未连接的微孔隙和固体颗粒进行标记，剩余的灰度值用于根据孔隙率将连接的微孔颗粒分割为12个微孔相（即标签[37]）。由于在固体/Viton界面存在分割问题，我们选择了一个尺寸为1144 × 1144 × 6000的子体积，以避免这些问题，同时仍保持足够大的问题尺寸以进行可扩展性和性能评估。图2显示了分割后的图像切片。从图像中得到的样本孔隙率为0.34，使用Keller PD-33X差压传感器[37]测量的渗透率为25 mD。分配给微孔标签（1–12）的孔隙率和渗透率值取自Menke等人的已发布的Esteillades数据集[37]。在那项工作中，使用Zeiss Ultra Nano-CT从代表性的微孔区域（标签4，?=0.42）钻取的65μm直径子样本中获取了纳米CT图像，并对其进行了分析以提取微结构描述符（例如，孔隙和颗粒大小统计）。这些描述符用于参数化基于对象的孔隙网络生成器，生成了具有不同孔隙率的合成微孔实现。对这些实现的渗透率模拟得出了孔隙率-渗透率关系，Menke等人[37]利用该关系根据每个标签的孔隙率为其分配渗透率值（表1）。对于已解析的孔隙空间和固体相，我们假设Kf?1=0和Kf=10?6 mD。

下载：下载高清晰度图像（1MB）下载：下载全尺寸图像图2. 分割为14个标签的Esteillades微CT图像的中间切片。

表1. Esteillades微CT图像每个分割相（即标签）的总体积分数、孔隙率和渗透率。
标签 分数 ? Kf (mD)
0 0.1 0.0
1 0.1 1.0
∞ 10.1 7.5
2 0.5 7.5
3 0.0 6.0
4 0.5 7.0
5 0.4 6.0
6 0.3 7.0
7 0.0 4.8
8 0.4 5.4
9 0.0 4.8
10 0.3 5.0
11 0.2 6.5
12 0.0 7.0
13 0.0 1.9
14 0.0 0.1

2.5. 并行网格化
CFD的第一步是构建计算网格。在OpenFOAM中，网格 définit dans le dossier constant/polyMesh d'un fichier de cas[9]。此目录包含多个文件，描述计算网格的几何形状、连通性和属性。首先，points fichier liste toutes les sommets (c'est-à-dire les coordonnées 3D) du réseau. Ces points forment les coins des blocs de grille de calcul. Ces points sont partagés par les cellules voisines et définissent la géométrie. Pour les blocs de grille hexahédriques, chaque cellule est composée de 8 points, et chaque point peut appartenir à jusqu'à 8 blocs de grille. Ensuite, face fichier décrit les faces du réseau en fonction des indices de sommets. Pour les blocs de grille hexahédriques, une face est composée de 4 points. Les faces peuvent être classées en faces internes (appartenant à deux cellules voisines) et faces de bordure (appartenant à un bloc de grille et à la frontière du domaine).

Dans OpenFOAM, les faces internes sont indexées en premier avant toutes les faces de bordure. Les fichiers owner et neighbour affectent les faces aux cellules : owner spécifie la cellule qui “possède” chaque face, tandis que neighbour spécifie la cellule voisine, uniquement pour les faces internes. Les faces de bordure, qui n'ont pas de cellules voisines, ne apparaissent que dans le fichier owner et sont détaillées davantage dans le fichier boundary, qui spécifie les conditions de bordure et les régions. L'ordre strict des faces internes en premier assure l'alignement avec le fichier neighbour, car il dépend de l'ordre successif pour définir la connectivité. Il est important de noter que la représentation d'un seul réseau n'est pas unique, car elle dépend de l'ordre des points, des faces et des cellules, qui est défini par l'algorithme.

Dans un réseau parallèle (MPI) d'OpenFOAM, les blocs de grille sont divisés en sous- domaines, chacun défini dans son propre dossier processeurJ, où J est le rang du processeur. Chaque dossier de processeur contient son propre sous-dossier constant/polyMesh avec des fichiers de grille locaux. Par conséquent, le réseau MPI peut être considéré comme une collection de sous-réseaux, un pour chaque processeur. Le fichier points définit les sommets du réseau du processeur de manière indépendante. Les fichiers faces, owner et neighbour décrivent les faces des sous-réseaux en fonction des faces internes et des faces de bordure. Les faces partagées entre les blocs de grille appartenant à différents processeurs sont des faces internes dans le réseau global, mais apparaissent comme des faces de bordure dans les sous-réseaux du processeur. Ces limites supplémentaires sont identifiées comme procBoundary dans le fichier boundary. De plus, quatre fichiers supplémentaires — cellProcAddressing, pointProcAddressing, faceProcAddressing et boundaryProcAddressing — sont utilisés pour mapper les indices dans les sous-réseaux du processeur à ceux du réseau global. Ces fichiers permettent une communication et un échange de données efficaces entre les sous-domaines lors des calculs parallèles.

Dans un réseau non structuré, il n'y a pas de relation directe entre les indices des points et leurs coordonnées, ce qui rend difficile la construction rapide et efficace du réseau, en particulier en parallèle. Cependant, notre approche numérique permet l'utilisation d'un réseau cartésien structuré, où les indices des points sont directement liés aux coordonnées. Cette capacité est particulièrement avantageuse pour le traitement parallèle.

Dans notre approche, les sous-réseaux du processeur sont définis de manière régulière, comme le montre la figure 3. Chaque processeur peut alors construire son propre réseau local indépendamment. De plus, il est simple de calculer l'indice global d'un point de réseau en fonction de son indice local dans son propre processeur et du rang du processeur. Ainsi, chaque processeur peut générer ses propres fichiers de grille locaux, y compris cellProcAddressing, pointProcAddressing, faceProcAddressing et boundaryProcAddressing, facilitant une haute scalability.

Pour déterminer la correspondance entre les indices globaux et locaux, nous considérons une image micro-CT de dimensions NX×NY×NZ, et une décomposition parallèle définie par NPX×NPY×NPZ processeurs dans les directions x, y et z respectivement. Pour chaque direction J∈{X,Y,Z}, si NJ est entièrement divisible par NPJ, alors le nombre de cellules par processeur dans cette direction, NJPP, est constant et égal à QJ=NJ/NPJ.

下载：下载高清晰度图像（357KB）下载：下载全尺寸图像图3. Cartésien régulier décomposé en sous-domaines. Le domaine complet (gauche) est divisé en sous-domaines de taille égale (droite), chacun attribué à un processeur. Les indices des processeurs augmentent d'abord dans la direction x, puis y, et enfin z, permettant une correspondance prédicable des sous-domaines.

Dans les cas où le nombre de processeurs NPJ ne divise pas entièrement le nombre de cellules NJ, nous définissons le reste de la division comme r, de sorte que NJ=QJ×NJP+r. Dans cette situation, les premiers r processeurs dans la direction J re?oivent des sous-blocs qui contiennent une ligne de cellules de plus (QJ+1) que les autres processeurs (QJ). Cela entra?ne une distribution non uniforme des données. Cependant, le schéma d'indexation permet à chaque processeur de calculer les indices globaux de ses cellules locales directement à partir de son rang et de ses indices locaux, sans nécessiter de communication.

La dernière étape du processus de création de réseau consiste à définir les champs de porosité et de perméabilité locaux. Dans notre approche, ces champs sont représentés par deux objets volScalarField : eps, correspondant à la porosité locale, et Kinv, représentant l'inverse de la perméabilité locale. Dans OpenFOAM, un volScalarField est un champ défini sur les centres des cellules (arrangement colocalisé) du réseau de calcul, utilisé pour stocker les variables principales telles que la pression, la vitesse et la concentration (arrangement colocalisé), ainsi que d'autres champs comme la porosité et la perméabilité. Les flux de surface sont définis sur les faces des cellules en tant qu'objets surfaceScalarField (par exemple, ?=U?S). La couplage pression-vitesse est géré par la procédure STANDARD SIMPLE pour éviter la séparation pression-vitesse sur les réseaux colocalisés[9].

Les objets volScalarField eps et Kinv sont initialisés en lisant leurs valeurs depuis le dossier 0 du cas d'essai au début de la simulation. Lors des exécutions parallèles, les deux champs sont décomposés de sorte que chaque processeur initialise ses cellules internes et ses conditions de bordure — y compris procBoundary — en lisant des données depuis le dossier correspondant procBoundaryJ/0. Chaque fichier d'entrée contient une liste de valeurs correspondant aux valeurs des cellules internes, ordonnées selon le numéronement interne des cellules, ainsi que les valeurs pour procBoundary et toutes les limites globales attribuées à ce processeur.

En utilisant un réseau cartésien structuré qui correspond précisément aux dimensions de l'image micro-CT, ainsi qu'une décomposition parallèle régulière, nous assurons une correspondance directe entre le numéronement interne des cellules locales et des cellules de chaque processeur et leurs indices correspondants dans l'image micro-CT. Par conséquent, chaque processeur peut générer ses propres fichiers eps et Kinv, ce qui soutient une implémentation parallèle hautement scalable.

Chaque processeur ouvre l'image micro-CT, qui est stockée dans un format .raw avec les dimensions NX, NY et NZ. Ce format correspond à NZ couches, où chaque couche est un tableau bidimensionnel de taille NX×NY. Afin de minimiser l'utilisation de la mémoire et les co?ts d'E/S, l'image est lue couche par couche, et un processeur accède uniquement à une couche donnée si elle contient des cellules locales internes ou des limites. Les indices globaux de ces cellules sont ensuite mappés sur les indices locaux correspondants dans le processeur et stockés dans un tableau tridimensionnel de taille NXPP×NYPP×NZPP, reflétant le nombre de cellules attribuées à chaque processeur. Chaque étiquette de voxel est ensuite traduite en valeurs correspondantes pour eps et Kinv. Pour l'ensemble de données Bentheimer, les étiquettes 1, 2 et 3 représentent des régions solides (??f=0.0001), tandis que l'étiquette 4 représente l'espace poreux (??f=1.0). L'image est smoothingée, et le champ de perméabilité est calculé à l'aide de l'équation (18).对于Esteillades样本，体素标签根据表1.2.6中提供的值映射到孔隙度??f和渗透率K。作为参考，我们在这里定义了通常用于具有无滑移边界条件的二元孔隙/固体图像的孔隙尺度DNS的标准OpenFOAM工作流程（图4）。从一个分割的图像开始，首先使用等值面提取方法将孔隙-固体界面重建为三角化表面（例如，STL）。然后生成一个笛卡尔背景网格（例如，使用blockMesh），并对其进行分区以进行并行执行（decomposePar）。在这个阶段，网格是一个结构化的体素对齐的格子，因此网格实体（单元/面/点）遵循可以直观映射到体素索引（i,j,k）的确定性索引。随后，使用snappyHexMesh将网格转换为孔隙空间中的边界一致网格，该网格移除了固体区域中的单元，并通过迭代细化、调整和变形网格来使单元面与重建的界面对齐，同时满足网格质量约束[9]。在这些操作之后，网格在索引意义上不再是一个结构化的格子：单元编号和连通性不再直接映射到体素索引，需要全局网格操作来保持连通性和质量。相比之下，本工作中使用的DBS工作流程保留了一个体素对齐的结构化网格，并避免显式的表面重建和边界一致的孔隙空间网格划分。

在实践中，这种标准工作流程在中等问题规模下就已经遇到了可扩展性限制，因为blockMesh是串行执行的，必须在单个进程上构建完整的背景网格。随着目标分辨率的提高，这个串行阶段的内存占用和运行时间可能成为瓶颈，甚至在应用边界一致步骤之前。因此，在标准工作流程中实现6003基准测试时，我们生成了一个更粗糙的背景网格（例如1503，即粗2-3倍），使用decomposePar对其进行分解，然后在snappyHexMesh中并行执行所需的细化。这个过程产生了有效的6003背景分辨率，同时避免了在单个串行步骤中构建完整精细背景网格的需要。最后，在具有无滑移条件的结果孔隙空间网格上执行流求解器。

下载：下载高分辨率图像（325KB）
下载：下载全尺寸图像

图4. 标准OpenFOAM边界一致工作流程的示意图（表面重建到STL，blockMesh，decomposePar，snappyHexMesh）与提出的DBS工作流程的比较，后者直接并行构建体素对齐的结构化网格并初始化eps。

3. 结果
在本节中，我们通过直接与标准OpenFOAM方法比较以及强适温和弱适庞性研究来评估我们工作流程在网格生成和求解器性能方面的效率。这里，标准OpenFOAM方法指的是第2.6节中总结的边界一致的孔隙空间工作流程；强适温和弱适庞性遵循下面给出的标准HPC定义。这些基准测试使用了之前描述的Bentheimer和Estaillades样本的裁剪子体积。由于目标是在不超过1000个节点的情况下模拟完整图像，因此选定的配置对于Bentheimer限制为每个核心最多约320,826个单元，对于Esteillades限制为每个核心最多61,347个单元。

虽然网格生成算法是确定性的，但其在ARCHER2上的运行时间受到文件系统I/O性能的显著影响。我们观察到I/O开销可能占总网格生成时间的10%到99%，主要是由于共享Lustre文件系统的竞争[45]。例如，在65,536个核心（512个节点）上运行Bentheimer样本的网格生成步骤五次，总运行时间分别为545秒、52秒、221秒、612秒和3720秒。然而，在所有情况下，实际计算时间都在25到26秒之间，其余时间主要由文件系统性能的变化所主导。为了确保公平和一致的比较，我们只报告计算网格生成时间，不包括I/O时间。

流模拟使用标准的半隐式压力相关方程方法（SIMPLE）和欠松弛进行[35]。每次非线性迭代涉及解决一系列线性系统。OpenFOAM中提供了两种线性求解器：几何-代数多网格（GAMG）求解器和预处理双共轭梯度稳定（PBiCGStab）求解器[9]。求解器（simpleDBSFOAM）、测试案例和参数设置的完整实现细节可以在GeoChemFoam项目仓库中找到，网址为https://github.com/GeoChemFoam。

为了与标准方法和强适庞性测试进行比较，模拟运行到稳态，定义为压力残差低于10^-4。在每次SIMPLE迭代中，线性系统解决到相对容差0.1。传输求解器在100个固定时间步长Δt=0.001秒内执行，相当于平均CFL数0.05。对于每个时间步长，使用PBiCGStab求解线化方程，直到残差降低到10^-9。

在强适庞性测试中，由于所需的计算量保持不变，强适庞性并行效率Effs定义为：
Effs = (20) Effs = TRef / (TNP / NP)，
其中T是模拟的计算时间，NP是使用的处理器（核心）数量，而TRef和NP,Ref是基线情况的值。理想情况下，基线情况在1个核心上运行，或者由于资源限制，在可能的最低核心数量上运行。对于以下强适温和弱适庞性研究，由于内存限制，使用了单个节点（128个核心）作为基线情况。这个选择反映了目标的大规模Bentheimer模拟的实际限制，该模拟至少需要100个节点的内存容量，同时也符合大规模OpenFOAM案例的现实部署场景。相应的记忆需求估计为每个单元大约300字节，包括每个单元约200字节的网格拓扑和几何数据（点、面、所有权者-邻居连接性和几何场）以及每个单元约100字节的场变量和求解器数据。这些数字与之前报告的OpenFOAM内存和可扩展性特性[46]、[47]、[48]、[49]一致。使用这个完整的节点作为并行效率基线，使我们能够在可用资源和时间的限制内评估并行性能。

在弱适庞性测试中，随着核心数量的增加，域的大小也会增加，求解器的收敛行为可能会因为几何复杂性的变化而受到影响。为了将并行效率与收敛变异性分开，我们固定了每个时间步长的迭代次数：GAMG为一次迭代（在相对容差0.1时足以近似收敛），PBiCGStab为30次迭代（以达到可比的精度）。传输求解器同样使用PBiCGStab在每个时间步长执行5次固定迭代。这种固定迭代策略确保了每个时间步长的计算工作量保持不变，使我们能够独立于算法收敛性来评估并行性的影响。弱适庞性并行效率Effw定义为：
Effw = (21) Effw = TRef / T。

3.1. 与标准方法的比较
为了评估我们提出工作流程的实际优势，我们将其与OpenFOAM中通常使用的基于网格的传统方法进行了比较，后者涉及生成一个符合形状的网格。DBS公式在完整的体素域上解决流动和传输问题，这在中等规模下相对于仅考虑孔隙的无滑移DNS会增加每个时间步的成本。这种方法的一个明显优势是它能够直接整合欠解析的孔隙度，如Esteillades样本中所做的那样。然而，这一小节的关键信息是，当图像大小超出边界一致预处理仍然可行的范围时，这种方法也对二元孔隙/固体图像变得有利。

因此，我们从一个中等大小的基准测试开始：Bentheimer的6003子体积（2.16亿个体素），从样本的中间提取。这个基准测试的规模比第4节中针对的数十亿体素域小两个数量级。随后，我们将域的大小增加到10403个体素，以展示所提出的工作流程如何扩展到更大的问题，并突出标准边界一致方法的局限性。

对于每种工作流程，首先生成网格，然后使用GAMG线性求解器解决流动直到收敛。表2总结了在不同节点和核心数量（从1个节点（128个核心）到5个节点（640个核心）下测试的两种工作流程的计时结果。

标准工作流程涉及多步网格生成过程：blockMesh定义一个粗糙的背景网格，decomposePar将域划分以便并行执行，snappyHexMesh将网格与几何形状对齐。最初的blockMesh步骤在单个核心上运行，需要大量内存，这阻止了直接生成6003背景网格。为了避免这一限制，首先创建了一个1503块结构的更粗糙的网格，然后使用decomposePar进行划分，最后在snappyHexMesh的并行执行期间将其细化四倍，有效地生成了一个6003背景网格。snappyHexMesh工具使用STL表面表示法来移除固体相内的网格单元，然后迭代地细化和变形网格，将单元面调整到界面以准确捕获固体-流体边界[9]。

表2. 标准工作流程和新工作流程之间的模拟时间比较。标准方法使用blockMesh、decomposePar和snappyHexMesh进行网格生成。新方法直接使用结构化笛卡尔网格并行构建网格。

核心/节点
标准方法
新方法
Empty Cell
raw2
stl（秒）
blockMesh（秒）
decomposePar（秒）
snappyHexMesh（秒）
createMesh（秒）
createEps（秒）
128/11
262
367
33
83
33
21
256/21
262
377
27
31
71
11
384/3
126
2396
257
107
51
12/4
126
2310
52
51
85
64
0/5
126
2311
62
376

尽管进行了这种细化，标准工作流程仍然存在几个可扩展性限制。decomposePar工具虽然旨在支持并行性，但实际上是串行操作的，并且随着分区数量的增加而变得缓慢。这个瓶颈在大规模模拟中尤其成问题。此外，snappyHexMesh包括几个仅部分并行化的阶段，它们并不随域的大小或处理器数量良好的扩展。因此，对于大型域来说，网格生成仍然耗时且效率低下，当域的网格需要超过单个共享节点上的可用内存时，在计算上变得不可行。

相比之下，我们的工作流程直接使用与体素化输入图像对齐的结构化笛卡尔网格并行构建网格。这消除了所有串行预处理步骤。每个核心独立构建其部分的网格，无需进程间通信。因此，网格生成随着核心数量的增加而高效扩展，从128个核心的54秒减少到640核心的仅10秒。这种扩展行为在3.2强适庞性分析和3.3弱适庞性分析中进一步分析。

求解器运行时间和收敛行为在表3中进行了比较。虽然我们方法中的流求解器时间一致更长，但这是在保留整个体素域（包括固体和流体相）的预期结果。在标准工作流程中，网格仅包含流体区域，在我们的案例中大约有4000万个单元（图5a）。相比之下，我们的方法解决了整个6003=2.16亿个单元，大约多5.4倍（图5b）。

求解器时间的增加反映了这种差异。在所有核心数量下，两种工作流程的流动时间比率在5到6之间，这与单元数量的增加一致。重要的是，尽管每个核心的求解器时间更高，但我们的方法仍然能够有效地扩展资源。相比之下，标准工作流程由于串行网格生成阶段而显示出有限的可扩展性。

下载：下载高分辨率图像（517KB）
下载：下载全尺寸图像

图5. 计算域：（a）标准方法（4000万个网格块）；（b）新方法（2.25亿个网格块）；（c）新方法中的流体和界面单元（5300万个网格块）。

对于两种方法，非线性迭代次数几乎保持不变，表明域分解对收敛性没有显著影响。有趣的是，新方法始终需要大约10%更少的非线性迭代次数，这可能是由于使用了结构化笛卡尔网格，避免了单元偏斜并确保了均匀的间距。然而，每次非线性步骤的线性迭代次数更多，这是由于系统中单元和未知数的数量更多。尽管有这些差异，所有情况下计算出的渗透率是相同的，验证了新方法的数值准确性。

为了更好地说明标准工作流程的可扩展性限制，我们评估了可以根据节点数量网格化和模拟的最大域大小。图6显示了每种方法成功处理的最大体素域大小。对于标准工作流程，我们观察到大约2.25亿个体素（6083）的硬性上限，超过这个限度后，由于内存限制，snappyHexMesh会失败。这个限制与节点数量无关，表明边界拟合预处理的关键阶段（如单元格移除、对齐等）并没有完全分布式处理，而是涉及非局部的网格操作，这些操作会对内存使用产生上限。表3. 标准工作流程和新工作流程下求解器收敛性的比较。流动时间反映了单相稳态求解器的运行时间。

**核心数** | **标准方法** | **新方法** |
| --- | --- | --- |
| 128 | 12109 | 16412.4×10^-12 |
| 256 | 259721558 | 28972.4×10^-12 |
| 302 | 18681987 | 38432.4×10^-12 |
| 384 | 197387 | 5122.4×10^-12 |
| 512 | 2329721572 | 26872.4×10^-12 |
| 640 | 175971163 | 113187 |
| 1131 | 11979 | 1987 |

相比之下，我们的工作流程在处理最大可处理域大小时几乎呈线性增长，与节点数量成正比。使用五个节点，我们成功处理了一个超过10亿体素（10^403）的域。这证实了网格创建是完全分布式的，仅需要本地内存，不需要全局网格组装。这种可扩展性使我们的方法非常适合大型岩石样本的模拟，直接实现了我们能够进行十亿体素图像计算的目标。

**图6.** 标准工作流程和新工作流程下，可处理的最大图像大小（以十亿体素为单位）随节点数量的变化。标准方法的处理能力大约限制在2.25亿体素。新工作流程则与节点数量成线性增长。

尽管由于网格规模较大，求解器仍然较慢，但可以采用几种策略来提高性能。首先，增加核心数可以显著加快速度，特别是在网格构建期间，因为具有很强的可扩展性。其次，可以使用subSetMesh等工具在网格生成后移除固体区域，只保留流体和界面单元格（见图5c）。然而，这一步骤需要大量内存并且扩展性较差，因此我们在当前工作流程中避免了它。相反，未来的工作将重点放在局部网格适应和多分辨率策略上，例如在主要是固体区域进行粗化处理，同时在孔隙空间和界面附近保持体素级分辨率，以减少自由度而不牺牲精度。

**总结来说，** 虽然标准方法产生的网格较小且求解器运行速度较快，但它存在可扩展性差和内存限制的问题。我们提出的工作流程通过完全并行的网格生成和分布式内存使用消除了这些限制。该方法为非常大的域的高分辨率模拟打开了大门，远远超出了当前网格工具的极限。

**表4.** Bentheimer子域（1950×1950×30单元格）的强扩展性结果。并行效率是相对于单节点配置的。**

| 核心数 | 单元格/核心 | 网格生成（GAMG） | 网格生成（PBiCGStab） | 运输（某算法） | 空单元格 |
| --- | --- | --- | --- | --- | --- |
| 128 | 11,485,352 | 49.8 | 1.00 | 7.50 | 1.00 |
| 256 | 2742,676 | 24.8 | 1.00 | 3.33 | 1.02 |
| 302 | 19222.4×10^-12 | 1558 | 2.4×10^-12 | 2597 | 2.4×10^-12 |
| 384 | 19872.4×10^-12 | 3873 | 1987 | 2.4×10^-12 | 512 | 2.4×10^-12 |
| 512 | 2329721572 | 2687 | 2.4×10^-12 | 1426 | 2.4×10^-12 |
| 640 | 175971163 | 1131 | 1985 | 2.4×10^-12 | 640 | 1.75 |
| 1131 | 11979 | 1985 | 2.4×10^-12 | 640 | 1.75 |

**图6.** 显示了标准工作流程和新工作流程下，可处理的最大图像大小（以十亿体素为单位）随节点数量的变化。标准方法的处理能力大约限制在2.25亿体素。新工作流程则与节点数量成线性增长。

**表5.** Estaillades子域（1144×1144×150单元格）的强扩展性结果。并行效率是相对于单节点配置的。**

尽管由于网格规模较大，求解器仍然较慢，但可以采用几种策略来提高性能。首先，增加核心数可以显著加快速度，特别是在网格构建期间，因为具有很强的可扩展性。其次，可以使用subSetMesh等工具在网格生成后移除固体区域，只保留流体和界面单元格（见图5c）。然而，这一步骤需要大量内存并且扩展性较差，因此我们在当前工作流程中避免了它。相反，未来的工作将重点放在局部网格适应和多分辨率策略上，例如在主要是固体区域进行粗化处理，同时在孔隙空间和界面附近保持体素级分辨率，以减少自由度而不牺牲精度。

**3.2. 强扩展性分析**
为了评估强扩展性，我们考虑了由于内存限制在单核上可以处理的z方向上的最大子体积（即切片数量）。对于Bentheimer域，这对应于1950×1950×50；对于Estaillades，则是1144×1144×150。强扩展性测试涉及在保持问题规模不变的情况下增加核心数（即Bentheimer为190百万单元格，Estaillades为196百万单元格）。随着核心数的增加，每个核心的单元格数量减少，从而可以评估每个求解器在固定工作负荷条件下的资源利用效率。

**表4和表5**展示了每种配置下以及模拟每个阶段的执行时间和相应的并行效率。结果包括网格生成（总秒数）以及流动和运输求解器（每时间步的秒数，SPT）及其各自的并行效率。图7展示了两个域的强扩展性趋势。在理想情况下，执行时间随核心数线性减少，从而实现并行效率为一。偏离这一理想情况反映了通信开销、负载不平衡和算法限制的影响。虚线表示架构无法扩展到超过1000节点分配限制的极限。

**表4和表5**还展示了每种配置下网格生成阶段的效率。效率在整个核心数量范围内保持稳定，甚至在超过10,000个核心时也是如此。这种行为突显了网格构建期间有限的通信需求，因为一旦域被分解，大部分通信都是局部的。尽管并行网格生成引入了额外的数据结构（如procAddressing和procBoundary），但这些结构只初始化一次，并且不会随着核心数的增加而增长，从而保持性能的一致性。

在低到中等核心数时，求解器的并行效率呈现超线性增长。这种行为主要归因于缓存效应和内存利用的提高。随着问题分布在更多处理器上，每个核心上的工作集变小，可以更好地适应缓存层次结构，减少内存访问延迟并提高算术吞吐量。然而，随着核心数的继续增加，每个处理器上的单元格数量最终变得太少，无法抵消并行化引入的通信和同步开销。超过这个点后，进程间通信的成本（特别是在压力校正、多网格粗化和Krylov求解器简化等集体操作期间）开始占主导地位[50]。因此，效率开始下降，因为每个处理器的计算量不足以掩盖通信延迟。这一转变标志着给定问题大小的强扩展性的实际限制。然而，在我们的案例中，这一限制发生在ARCHER2的最大分配限制之外。

**表6.** Bentheimer案例的弱扩展性配置。每一行对应于每个核心固定单元格数量的配置（371,338）。计算时间以每时间步（STP）的秒数报告。并行效率是相对于单节点配置计算的。

**表7.** Estaillades案例的弱扩展性配置。每一行对应于每个核心固定单元格数量的配置（95,855）。计算时间以每时间步（STP）的秒数报告。并行效率是相对于单节点配置计算的。

**4. 多亿体素模拟**
基于弱扩展性研究中确定的配置，我们现在对两个代表性岩石样本（Bentheimer和Estaillades）进行了全核模拟：Bentheimer（1950×1950×10,800体素）和Estaillades（1144×1144×6,000体素）。这两个域都在体素尺度上进行了离散化，分别产生了大约410亿和80亿个单元格。根据之前建立的最佳负载平衡，我们为Bentheimer分配了864个节点（即每个核心371,338个单元格），为Estaillades分配了640个节点（即每个核心95,855个单元格）。

**图8.** 弱扩展性结果表明，**Bentheimer（每个核心371,338单元格）和Estaillades（每个核心95,855单元格）**在大多数配置下都能保持高效。网格生成显示出接近理想的可扩展性，即使在最大的核心数量下，效率也始终高于0.95。这种行为反映了网格算法的完全分布式特性，每个处理器独立构建其局部域，进程间通信最少。

**5. 多亿体素模拟**
基于弱扩展性研究中的配置，我们现在对两个代表性岩石样本进行了全核模拟：Bentheimer（1950×1950×10,800体素）和Estaillades（1144×1144×6,000体素）。两个域都是在体素尺度上离散化的，分别产生了大约410亿和80亿个单元格。基于之前建立的最佳负载平衡，我们为Bentheimer分配了864个节点（即每个核心371,338个单元格），为Estaillades分配了640个节点（即每个核心95,855个单元格）。这些配置在ARCHER2的1000节点限制范围内。

**结论**
弱扩展性分析的目的是评估系统在保持每个核心单元格数量不变的情况下，如何处理不断增加的问题规模。我们从强扩展性研究中选出的最大配置开始，然后调整域规模和核心数量，以探索系统在不同扩展场景下的表现，同时确保我们在ARCHER2的1000节点限制范围内。图(a)显示了速度幅值分布，按达西速度进行了缩放，证实了Estaillades中的分布范围更广。图(b)绘制了沿核心长度的平均浓度，突出了Bentheimer中的明显前沿以及Estaillades中更为平缓的过渡，这表明其扩散作用更强。下载：下载高分辨率图像（516KB）下载：下载全尺寸图像

图9. 圆柱形岩芯样本内的压力分布，通过裁剪右上象限来显示内部场。（a）Bentheimer；（b）Estaillades。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像

图10. 圆柱形岩芯样本内的速度分布，通过裁剪右上象限来显示内部场。（a）Bentheimer；（b）Esteillades。图13将模拟得到的渗透率与实验值进行了比较——Bentheimer的实验数据来自[41]，Esteillades的实验数据来自Menke等人[37]。结果显示两者吻合良好，Bentheimer的误差为10.4%，Esteillades的误差为19.5%。直接从图像数据得出的孔隙率也与实验值相符。下载：下载高分辨率图像（624KB）下载：下载全尺寸图像

图11. 注入0.5孔隙体积后的溶质浓度分布，通过裁剪右上象限来显示内部场。（a）Bentheimer；（b）Estaillades。下载：下载高分辨率图像（328KB）下载：下载全尺寸图像

图12. （a）按达西速度标准化的孔隙空间内的速度分布。（b）注入0.5孔隙体积后沿核心轴的平均溶质浓度。

为了研究域大小的影响，我们还在较小的子体积上计算了渗透率。Bentheimer沿z轴被分为三个相等的部分（每个部分1950 × 1950 × 3600个体素），而Esteillades被对半划分（1144 × 1144 × 3000个体素）。这些子域的孔隙率和渗透率也显示在图13中。尽管子体积模拟与实验结果仍然基本一致，但误差有所增加——Bentheimer的孔隙率误差高达2.82%，渗透率误差高达22.6%，尽管样本本身是均匀的。对于Esteillades，孔隙率的误差为3.0%，渗透率的误差为52.5%。这些差异，尤其是在异质性更强的样本中，强调了使用大域（数十亿个体素）进行准确表征的重要性。

这些结果证明了我们的求解器能够高效地执行数十亿个体素、高度异质多孔介质中的流动和传输的高保真模拟，采用了针对弱尺度优化过的并行配置。与实验数据的良好吻合证明了该方法的稳健性及其在多孔系统的高分辨率建模中的适用性。下载：下载高分辨率图像（292KB）下载：下载全尺寸图像

图13. Bentheimer和Esteillades样本的渗透率与孔隙率关系。蓝色标记代表全芯模拟结果，红色标记表示子体积，黑色标记显示实验数据。

5. 结论

我们提出了一种可扩展且高效的工作流程，用于在具有隐式固体质表示的体素对齐结构化网格上模拟单相流动和溶质传输。虽然我们的实现是在Darcy-Brinkman-Stokes框架内展示的，但核心方法依赖于避免流体-固体界面的显式重建，可以扩展到任何支持这种表示的方法。这种灵活性使得该框架能够稳健且通用地用于各种物理公式的高分辨率孔隙尺度建模。

结构化网格架构允许直接进行域分解，每个子域可以独立处理，通信开销最小。这种设计显著简化了并行实现，并且对于大规模模拟非常有效。与标准非结构化网格化工作流程相比，我们的方法显著减少了预处理时间和内存消耗，使其更适合极端规模的应用。性能评估表明，包括网格生成、流动和传输在内的整个模拟流程能够高效地扩展到数万个核心。特别是网格生成展示了理想的可扩展性，如果我们忽略共享系统上变化的I/O时间的话。虽然流动和传输求解器在大约10万个核心的数量级上保持了一致的性能，保持了至少50%的并行效率。

我们将该方法应用于模拟两个大型、真实的岩石样本——Bentheimer砂岩和Esteillades石灰岩——使用了多达410亿个单元格。模拟结果显示，在均匀的Bentheimer样本中捕获了清晰的溶质前沿，在异质性更强的Esteillades样本中显示了广泛的扩散现象，反映了它们各自的孔隙结构。与实验测量的比较表明，我们的全芯模拟结果与实验数据吻合良好：Bentheimer的渗透率误差在10.4%以内，Esteillades的误差在19.5%以内。然而，子体积分析显示了显著的偏差——Bentheimer的渗透率误差高达22.6%，Esteillades的渗透率误差高达52.5%——这强调了在足够大的尺度上进行建模以捕捉异质介质中的代表性行为的重要性。

总体而言，该工作流程能够在无需预处理或几何简化的情况下直接进行复杂多孔介质中的流动和传输的高保真模拟。它为数字岩石物理学、不确定性量化以及实验岩芯尺度数据的解释提供了一个实用且可扩展的工具。

未来的工作将集中在将这一框架扩展到支持多相流动和反应传输。这些扩展将引入更多的复杂性，包括相间的非线性耦合、相界面动力学以及潜在的强烈化学反应。应对这些挑战需要继续优化求解器组件、制定管理局部复杂性的适应策略，并支持新兴的计算架构（包括GPU加速）。从长远来看，这些发展旨在扩大基于图像的建模范围，以解决科学和工业应用中越来越真实和多尺度的多孔介质系统问题。

CRediT作者贡献声明

Julien Maes：撰写——原始草稿、软件、方法论、概念化。Gavin J. Pringle：撰写——审稿与编辑、软件、方法论。Hannah P. Menke：撰写——审稿与编辑、资金获取、数据管理、概念化。

在撰写过程中使用生成式AI和AI辅助技术的声明

在准备这项工作时，作者使用了chatGPT进行语言校正。使用该工具/服务后，作者根据需要对内容进行了审阅和编辑，并对已发表文章的内容负全责。