用于多孔结构变形建模的混合有限元-深度学习方法

《Applied Mathematical Modelling》：Hybrid finite element-deep learning method for deformation modeling of porous structures

　　郝宇|陈嘉轩|徐汉西|段奎成|王权|孟思伟|刘贺|吴恒安•同时模拟多尺度下多孔材料的变形行为。•通过深度学习预测的插值方法体现孔结构效应。•该预测框架由多个并行DeepONets构成。•该方法能够在较低时间成本下实现精确的变形估算。1. 引言多孔材料不仅存在于骨骼、岩石等自然结构中，也见于混凝土、泡沫等人造结构中[[1], [2], [3]]。这类材料通常具有各向异性刚度、显著的尺寸效应、较高的冲击韧性以及优异的能量吸收能力，因此应用范围极为广泛。因此，准确评估多孔材料的力学性能对于其合理应用至关重要[[4], [5], [6]]。数值模拟为研究内部应力分布、变形机制以及微观结构特征的影响提供了有力工具，而这些内容通过实验很难获取[7,8]。然而，传统数值方法在处理多尺度多孔材料时存在固有局限，当材料中含有大量孔隙或孔隙尺度远小于整体材料尺度时，计算成本会极高[9]。为准确捕捉与微小孔隙相关的应力-应变分布，模拟需要极其精细的网格，这又会带来巨大的计算负担。因此，亟需一种新的数值方法，以高效且精确地评估多孔材料的力学性能，并揭示其详细的力学响应。为解决多尺度问题，计算均匀化方法——主要包括基于渐近分析的均匀化方法（AABH）[[10], [11], [12]]和代表体积元方法（RVE）[[13], [14], [15]]——为评估材料的有效性能提供了有力工具。AABH以渐近展开理论为基础，这是一种经典的数学方法，通过小特征参数的级数分解来近似求解复杂的多尺度物理和力学问题[10]。在该方法中，宏观单元节点处的应力场、位移场及其他变量会围绕某一小特征长度尺度进行渐近展开[11]。基于节点附近的这种渐近展开，可建立相应的控制方程，进而求解宏观与微观的应力-应变场[12]。李等人[16]在能量一致的框架下提出了双尺度渐近方法，用于分析含有周期性微裂纹的脆性材料的动态变形与断裂行为。近期，郭等人[17]成功拓展了AABH的适用范围，他们提出了一种将非周期性结构映射为周期性结构的映射方法，并运用AABH来描述非周期性结构的力学行为。尽管AABH能够实现高精度求解，并可获取宏观与微观变量，但其求解过程计算复杂且耗时较长，尤其是在处理具有复杂孔结构的多孔材料时[18]。除AABH外，RVE方法也是一种常用的均匀化方法，在多孔材料计算中应用更为频繁[13]。基于RVE，计算过程分为两个部分：一是对RVE的局部材料本构关系进行微观计算；二是将材料视为由大量周期性RVE构成的均匀介质，进而求解宏观边界值问题（BVP）[14]。这种处理方式有助于高效计算异质结构，但RVE的选取显然会影响宏观BVP的计算精度。因此，正确构建并合理处理RVE是计算均匀化方法的核心。最为通用的RVE方法是多层有限元方法（FE2）[[19], [20], [21], [22]]，该方法利用有限元方法（FEM）同时处理微观层面的RVE和宏观层面的BVP[22]。具体而言，会将本构关系和周期性边界条件应用于RVE，通过FEM计算有效应力σ?和有效应变??。随后利用RVE中得到的σ?与??之间的关系，计算求解宏观BVP所需的刚度矩阵。不过，尽管FE2方法能够实现多尺度计算，并且可以高精度、灵活地计算RVE的各类性能，但由于需要进行大量计算密集型的迭代运算，尤其是当RVE的结构复杂多样时，其计算成本通常相当高[23]。为降低FE2方法的巨大计算量，近期研究采用了模型降阶（MOR）技术[24,25]。该策略通过降低具有周期性边界条件的微观RVE问题的计算复杂度，从而加快RVE有效性能的计算速度，但会一定程度上牺牲精度。例如，阮等人[26]提出了一种基于交互作用的材料网络，用于降低RVE的计算成本。他们在RVE中设计了基于交互作用的材料网络架构，并利用FEM结果在Hill–Mandel条件约束下回归网络中的拟合参数。经过训练的网络能够高效处理具有微观孔隙的结构[14]。程等人[27]提出了一种基于FEM-Cluster的分析方法（FCA）来实现MOR，该方法近年来得到了逐步的研究与改进[13,28,29]。他们利用K均值聚类根据应变状态将RVE中的庞大有限元网格划分为多个簇，然后基于FEM的数据构建降阶模型的交互矩阵。尽管这些MOR方法确实降低了计算成本，但它们通常会将原始控制方程简化为低阶表达式，并在计算中舍弃局部自由度。这种简化会导致精度下降，限制其准确捕捉局部力学响应的能力。深度学习（DL）作为一种强大的非线性拟合工具，近年来被用于描述多孔材料等微观结构与其有效性能之间的关系[30,31]。DL常作为传统计算方法的补充手段[[32], [33], [34]]，它不仅能够大幅提升计算速度，还能保证精度。例如，韦伯等人[33]提出了一种基于均匀化理论的CNN，用于预测多孔介质的渗透率和扩散张量，这种方法借助均匀化理论生成数据集，避免了在宏观孔隙尺度上进行直接数值模拟的高昂计算成本。肖等人[34]开发了一种基于微观力学的神经网络，用于对含有圆柱正交各向异性纤维的热导复合材料进行周期性均匀化处理，同时运用迁移学习加速训练。里莫利和卡普阿诺[35]则应用不同的DL方法预测RVE中的受力情况，无需复杂计算位移场，也不需要迭代运算。这种方法能够提供精确的预测结果，表明DL可作为RVE计算的替代模型，实现多尺度模拟。他们还利用人工神经网络（ANN），根据位移的时间域数据预测二维RVE的受力时间域数据[36]。该方法能够以极低的时间成本根据有效应变计算RVE的有效应力，从而便于分析复杂均匀体内的波传播现象。尽管现有的基于DL的方法通过快速预测RVE的应力状态显著提升了计算速度，但它们仍存在与MOR类似的缺陷。大多数方法仅能提供RVE的有效应力与应变信息，无法通过DL模型获取内部的力学响应数据[37]，因此难以捕捉应力集中等现象，而这些现象对于准确评估材料失效、损伤演化以及其他局部力学行为至关重要[38]。因此，有必要优化DL的整合方式，提出一种高效、精确且适用于不同材料性能的多孔材料计算方法。幸运的是，近期的一些研究为此提供了思路，它们利用DL模型处理FEM中的插值函数。黄等人[39]提出了一种与问题无关的深度学习方法，用于拓扑优化。他们将二维计算域划分为低分辨率单元，利用DL根据材料密度预测插值函数。该方法使拓扑优化的所需时间大幅减少了两个数量级。在后续研究中[40,41]，他们将子结构的维度提升至三维，有效解决了三维拓扑优化问题。这一系列研究为快速、精确且通用性强的大规模多尺度计算提供了新思路。不过，目前所采用的深度学习方法主要基于传统的ANN，这在处理更大规模的问题时可能会限制计算效率。与此同时，张等人[42]提出了一种基于深度学习的智能插值函数，用于加速有限元计算。他们应用多个并行DL架构，预测具有不规则体素形状的有限元内的插值函数。该方法提高了不规则几何体的模拟效率，能够精准描述局部响应，但仅限于二维场景。对于三维多孔材料，需要一种新方法来准确捕捉复杂的空间拓扑结构和几何异质性，从而可靠地预测子结构内的插值函数，进而高效计算力学响应。在本研究中，我们提出了一种基于DL和FEM的混合计算方法，用于快速、精确地开展三维多孔材料的力学计算。该方法包括宏观尺度计算和RVE尺度（包含孔隙的子结构）计算，通过RVE中的隐式插值函数实现两种尺度之间的信息交换，而这些插值函数是由一种新型DL模型预测得到的。在该模型中，提出了三维多重并行DeepONets，用于预测三维多孔材料RVE中的插值函数分布。同时，还引入了ResNet架构[43]，用于约束与插值函数导数场相对应的应变场，从而减少由于神经网络固有特性导致的应变场噪声。基于三维插值函数，可以构建RVE的刚度矩阵并组合起来，求解多孔材料的全局平衡方程，同时还可以通过插值函数和节点位移准确获得每个子结构内的内部位移场。通过多种涉及不同多孔结构及力学性能的数值算例，进一步验证了所提出方法的有效性。结果表明，该方法为利用深度学习进行多孔结构的有限元分析提供了有效框架。首先，它在多尺度结构的变形建模中保持了较高精度，位移误差通常相对FEM参考解在3%以内。其次，该方法大幅降低了计算成本，在所测试的多尺度算例中，计算时间较传统的FEM和FE2方法减少了三到四个数量级。第三，它同样适用于那些因可用计算内存限制而无法使用传统FEM处理的超大尺寸多孔结构。在岩土工程等实际工程应用中，比如变形建模[44,45]，该方法展现出巨大潜力。对于由多种RVE构成的复杂地质体，该方法能够大幅提升效率，准确捕捉不同储层块体的应力分布特征，为实际应用中的变形预测提供有力支持。本文的其余部分安排如下：第2节介绍所提出的计算方法以及用于预测三维插值函数的神经网络架构。第3节展示该方法的验证结果，这部分包括四个方面：DL模型预测性能的验证、具有周期性RVE的三维多孔结构的力学响应计算、具有非周期性RVE的三维多孔结构的力学响应计算，以及具有多种RVE的超大尺寸多孔结构的力学响应计算。第4节总结了所提框架的创新之处、结果的准确性与可靠性，并对未来工作进行了展望。2. 方法论2.1. 带隐式插值的有限元框架我们所提出的用于分析多孔材料力学响应的计算方法将计算分为两个尺度：宏观尺度与RVE尺度。值得注意的是，在宏观尺度求解过程中，RVE被视为一个计算单元，而其内部的孔隙微观结构则会被完整保留。在求解出宏观位移场后，便可进一步计算每个RVE内的力学响应。详细的计算框架如图1所示。如图所示，两个尺度之间的信息交换是通过RVE中的隐式插值来实现的。在宏观计算方面，利用隐式插值通过Galerkin方法构建RVE的刚度矩阵。将各个RVE的单元刚度矩阵组合成整体刚度矩阵后，再运用有限元方法求解宏观力学响应问题。而在RVE计算中，则是根据宏观计算得到的节点位移以及插值函数来计算RVE内部的位移场。因此，通过插值函数的信息传递，能够实现多孔材料的快速计算。下一节将详细介绍用于获取RVE插值函数的深度学习方法。下载：高分辨率图像（395KB）下载：完整尺寸图像图1. 所提出的多孔材料计算方法，该方法通过深度学习模型预测的RVE中的隐式插值函数，将宏观尺度和RVE尺度上的力学响应计算相结合。考虑到多孔介质材料中的RVE，采用小位移和小应变假设，并假设计算域Ω内不存在体力分布，其控制方程和边界条件的加权残差形式可表示为[46]：(1)R=∫WΩ(?·σ)dΩ+∫ΓNWN(T?σ·n)dΓN+∫ΓDWD(u?u￣)dΓD，其中σ、ε和u分别表示RVE中点x处的应力、应变和位移，ΓD表示施加狄利克雷条件的边界部分，ΓN表示施加诺伊曼条件的边界部分，T为牵引向量，n为边界的向外法向向量，WΩ、WN和WD分别表示对应域和边界上与控制方程及边界条件相关的权重函数，它们具有相同的函数形式。通过加权残差公式，通过确保残差在域上的积分经过测试函数加权后为零，从而近似得到精确的控制方程，其形式可转化为：(2)R=∫ΩWΩ(?·DLu?)dΩ+∫ΓNWN(T?DLu?·n)dΓN，(3)uˇ=∑Niui=Nu?，其中L是将位移场与应变场联系起来的微分算子，D表示材料本构矩阵，ǔ表示基于插值函数的插值位移场，ui为节点位移。需要注意的是，方程(3)中的插值函数Ni是二阶张量而非标量，RVE内部结构的各向异性导致不同方向上的插值函数存在差异。由于位移场在不同空间方向上的插值特性不同，因此需要为每个方向构建不同的插值函数[42,47]，其表达式如下式所示：(4)Ni=[N?iu000N?iv000N?iw]，且N?iu≠N?iv≠N?iw。将Galerkin方法[48]应用于方程(2)并进行部分积分，即对加权残差和位移插值使用相同的插值函数，方程(2)可表示为：(5)R=?∫(?·N)TDLuˇdΩ+∫ΓNNTTdΓN。对于多孔材料，解域的边界包括内部边界和外部边界，前者代表固体基质与内部孔隙空间的界面，后者代表整体多孔结构与外部环境的界面。如前所述，本节的重点是高效处理具有复杂孔隙几何结构的材料的力学响应计算，因此暂时忽略内部孔隙中的气体压力对结构变形的影响。方程(5)可简化为：(6)R=?∫(LN)TDLNu=dΩ+∫ΓNNTTdΓN。当该残差等于0时，表明控制方程所描述的弹性力学问题在整个域内大致得到满足。刚度矩阵K与节点力向量T?之间的关系可如下获得：(7)Ku?=T?，(8)K=∫(LN)TDLNdΩ，(9)T==∫ΓNNTTdΓN。通过上述步骤，可以利用RVE上的插值函数建立相应的平衡方程。该方法与FE2方法（如附录A所述）的主要区别在于，FE2方法在RVE尺度上计算有效的应力-应变数据对以辅助宏观尺度计算，而该方法则利用插值函数在两个尺度之间传递计算信息。因此，FE2方法只能得到将局部孔隙视为均匀体的宏观力学响应结果，而该方法则能同时获得宏观力学响应以及RVE内部任意位置的力学响应。值得注意的是，上述推导对插值函数提出了特定要求。首先，这些函数应自动满足克罗内克δ函数的条件，即在相应节点上插值函数的值为1，而在其他所有节点上值为0。其次，为确保在刚体平移作用下不存在内部变形，插值函数需满足单位分解条件，即其在整个域内的总和为1。由于不同的RVE包含大小和空间分布各异的孔隙，从而导致复杂的结构，其插值函数没有显式表达式，而是以所提出的神经网络预测的空间标量场形式来表征，RVE中的每个点都被赋予相应的函数值。具体的预测方法将在下一节中详细阐述。此外，本节仅考虑所有角节点都位于固体骨架内的多孔结构，未来相关工作中可以实现更复杂形状的多孔结构。2.2 三维隐式插值函数预测模型在前一小节中，介绍了基于加权残差形式的插值函数在RVE尺度上的计算方法，旨在提高其计算效率。然而，由于孔隙结构的复杂性，要在RVE上建立插值函数颇具挑战性，难以给出准确合理的解析表达式。因此，借助深度学习这一强大的工具来提取变量映射关系，从RVE结构中导出以标量场形式的隐式插值函数。需要注意的是，RVE是具有孔隙的三维结构，因此神经网络需要合适的架构来快速准确地处理相关信息。同时，由于插值函数的导数场与应变场呈线性关系，神经网络在预测一阶导数时需要保持稳定性。但由于需要处理的 ??较为复杂，神经网络需要更多层才能形成准确的映射，但过深的网络会导致输出结果出现白噪声，使其一阶导数出现显著波动。为解决上述两个问题，本文的神经网络结构采用了三维多并行DeepONets架构来快速处理三维矩阵信息，同时采用ResNet架构来稳定预测一阶导数，并避免网络层过深的问题。因此，在接下来的内容中，我们将详细介绍三维多并行DeepONets架构和辅助的ResNet架构，然后全面阐述插值函数预测网络。2.2.1 三维多并行DeepONets架构如前所述，由于人工智能网络强大的拟合能力，已被广泛用于变量映射[49]。然而，对于三维预测而言，人工智能网络的预测效率非常低。相比之下，具有小核函数的标准卷积神经网络主要关注局部特征，若没有深度堆叠或特殊架构，要有效捕捉整个三维RVE的全局相关性仍面临挑战[50,51]。因此，受Zhang等人[42]工作的启发，本文提出了一种基于三维架构的多并行神经网络。当处理列向量输入I时，人工神经网络的操作可表示为如下形式：(10)Ii+1,j=act(∑l=1niwlIi,l+bj)=Λi,j°Ii，其中Ii,j表示第i层第j个神经元上的值，wl为权重，bj为偏置项，ni为该层的神经元数量。运算符Λi:Ii∈Rni→Ii+1∈Rni+1表示神经网络中相邻两层之间的可训练映射关系，○则表示ANN中的一层。正如Zhang等人[42]所讨论的，多并行（MP）架构是基于带有偏置层的卷积神经网络。在本研究中，该MP架构被扩展到三维，以便处理三维形式的RVE数据。具体而言，采用了带有偏置层的三维1×1卷积神经网络，其表达式为：(11)Ii+1,j,m,s,t=act(∑l=1ni,jCi,lhi,l,m,s,t+1×Cni,j+1)=C°Ii,m,s,t，其中I为五维特征图，m、s和t表示三维RVE数据的形状。索引i表示第i个特征图，j表示第j个通道，这里的通道指的是使用相同特征图但由不同卷积核计算出的不同输出。C为五维卷积核，其作用与ANN中的wl类似，1为来自偏置层的输入。该公式表明，带有偏置层的三维1×1卷积神经网络在数学形式上与ANN相同，但通过利用卷积神经网络的并行处理能力，显著提升了计算速度。为简化表述，这种带有偏置层的三维1×1卷积神经网络可定义为以下运算符：(12)Ii+1=Ci°ˉIi，其中Ii∈Rni,j，ni,m,ni,s,ni,t分别表示第i个特征图中j、m、s和t维度的大小。基于此结构，常见的基于ANN的模型也可以转换为多并行形式，如图2(a)所示。下载：高分辨率图像（614KB）下载：完整尺寸图像图2. 三维多并行DeepONets（3D MP-DeepONets）架构示意图，包括(a)多并行人工神经网络（MP-ANNs），(b)DeepONet模块，以及(c)集成网络架构。该网络的目的是预测RVE中每个体素上插值函数的三维分布，因此需要一种能够准确采样各体素上预测插值函数的神经网络架构。DeepONet[52]是一种基于ANN的预测算法，它在样本点上采样预测的隐式函数，已被广泛应用于解决各种问题，包括偏微分方程的求解[[53],[54],[55],[56],[57]]以及高精度数据的恢复[[58],[59],[60],[61]]。如图2(b)所示，该架构包含两个网络（分支网络和主干网络），其中分支网络将输入I映射到隐式函数F的过程可简写为：(14)N1:I(x)?RnI→F?RnF，其中I和F为一维向量，nI和nF分别表示它们的维度。此处x表示一组坐标点，每个点对应输入特征I的位置。主干网络则将输入的采样点y转换为隐式函数，如下所示：(15)N2:y?Rd→H?RnH，其中y和H为一维向量，d和nH分别表示它们的维度。对于三维问题，d的值为3。通过计算分支网络和主干网络输出之间的内积，即可得到对应采样点y处的函数值F，其表达式如下：(16)F(y)=|N1(I(x))·N2(y)|，其中|?|表示两个向量之间的内积。尽管DeepONet能够准确预测三维RVE内部插值函数的分布，但它一次只能预测一个体素的一个值，因此预测过程极为耗时。为此，采用了MP架构来优化基于ANN的DeepONet架构。具体而言，将所有体素上的输入信息I及其对应的空间坐标X=(x1,x2,…,xn)整合在一起，以三维矩阵的形式输入。该矩阵随后通过图2(c)所示的三维MP-DeepONet架构进行处理，其数学形式如下：(17)F(X^)=|(C1°ˉI)·(C2°ˉX)|，其中C1和C2分别表示不同分支中的网络。三维MP-DeepONets结构显著提高了RVE内插值函数分布的预测效率，缓解了传统ANN在解决三维问题时计算效率低的限制。这一架构 thus为基于深度学习实现三维多孔材料的快速精确力学计算提供了模型基础。2.2.2 辅助ResNet架构如前所述，本文建立的深度学习模型的目的是预测RVE内部位移插值函数场的分布。为预测三维随机有限元模型中的插值函数场，采用了卷积神经网络架构构建三维MP-DeepONets结构，以提高计算效率。然而，卷积神经网络存在梯度破碎问题[62]。梯度破碎是指输出中会出现白噪声，进而严重影响输出函数的导数精度。对于本研究而言，这意味着输出的插值函数场无法准确描述应变场，因此亟需引入新的网络架构来解决这一梯度破碎问题。具体而言，梯度破碎指的是在通过梯度下降优化的前馈神经网络中，随着网络深度的增加，梯度之间的相关性呈指数级下降的现象。其根本原因在于随着网络深度增加，数据中的空间相关性逐渐减弱，最终转变为白噪声。这样一来，由插值函数计算的位移场在随机初始化条件下会含有大量白噪声，从而导致预测的应变场不准确。为解决因网络深度增加而引发的梯度破碎问题，及其对预测插值函数一阶导数场的影响，本研究在深度学习模型中引入了图3所示的ResNet[43]结构。与传统神经网络架构相比，ResNet的核心特点在于增加了跳跃连接。传统神经网络直接学习目标映射关系，而ResNet则学习输入与期望映射之间的加权残差。图3展示了ResNet中输入层与输出层之间的数学关系。当α和β均为1且采用前述相同初始化方式时，其衰减速率明显低于传统神经网络。此外，将β参数设定在0.1到0.3之间，并在神经网络函数中加入正则化处理，可进一步降低衰减速率，有效缓解因网络层加深带来的梯度破碎问题。这一改进解决了深度神经网络预测结果中相邻区域导数出现大幅波动的问题，使得能够利用损失函数约束下的预测值一阶导数。借助ResNet，可以构建出能够精确预测RVE中位移插值函数且具备准确一阶导数场的神经网络架构，该架构被应用于后续将详细介绍的深度学习模型中。如图6所示，该网络的训练过程可分为两个阶段。在这些图中，一个“时代”代表对训练集进行一次遍历的训练过程，而一个“步骤”则代表在十次遍历后对训练集和测试集进行测试的过程。因此，一个步骤包含十个时代，在每个步骤中都会计算三个方向上的位移R2值：(19)R2=1?∑i=1n(y^i?yi)2∑i=1n(yˉi?yi)2，其中R2表示决定系数，y^i表示预测值，yi表示真实值。该模型在包含8,000个样本的带标签训练集上训练1,000个时代，每次运行的总训练时间约为6小时。在训练过程中，前200个时代采用不同的学习率策略，后续阶段的学习率约为前几个时代的一百分比，从而有助于神经网络的收敛。实验结果表明，神经网络预测的插值函数能够准确预测三个方向的位移场。此外，在训练过程中测试集的结果没有出现显著下降，这说明神经网络尚未达到过拟合阶段。这表明当前网络具有良好的泛化性能。

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

图6. 用于预测插值函数的深度学习模型的训练过程。(a) 损失函数的变化情况。(b-d) 训练过程中预测的U、V、W位移场的R2值变化情况。

为验证经过训练的深度学习模型所预测的插值函数的准确性，如图7所示，使用预测的插值函数在测试集中计算位移场，并与有限元分析的计算结果进行比较。图7(a-c)展示了各坐标轴方向上预测位移场与有限元分析结果之间的R2值分布情况。可以看出，大多数预测的位移场与有限元分析计算的位移场分布十分吻合。图7(d)随机展示了一个三维多孔材料在X方向上U位移场的等高线图。无论是在表面还是多孔材料的二维切片上，深度学习方法计算的位移场都与有限元分析结果高度一致。

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

图7. 使用预测的插值函数计算的位移场。(a-c) 分别为预测的U、V、W位移场的R2值分布。(d) 实际与预测的U位移场。(e) 沿两条选定线条的实际与预测位移值对比。

为进一步说明预测结果之间的数值误差，从二维切片中提取了两条穿过孔隙区域的路径，并在图7(e)中展示了沿这些路径的位移曲线。结果表明，神经网络模型能够准确识别孔隙结构的位置，这在孔隙界面处表现为明显的不连续性，同时在固体骨架上的结果与有限元分析得到的结果非常接近。尽管该模型能够准确计算位移场，但使用深度学习模型计算单个规则化多孔单元的位移场的时间已从有限元分析方法的60秒缩短至0.01秒，速度提升了三到四个数量级。这些结果表明，神经网络模型预测的插值函数能够高效且准确地处理三维多孔材料中的位移场计算问题。

基于训练好的神经网络，可以得到位移场预测中隐含的插值函数。这些插值函数的三维分布等高线图如图8(a)所示。该等高线图对应于位于(1,0,0)节点处的位移U。在该节点处，插值函数值为1，当接近其他节点附近的表面时，其值会呈近似线性衰减。此外，这些三维插值函数还满足刚体位移的约束，因为在固体区域内的任意一点，这些函数值的总和均为1。

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

图8. 对应于位于(1, 0, 0)节点的预测插值函数。(a) 3D规则化多孔单元内的空间分布。(b) 沿三条代表性线条的函数值，其中线路1和线路2对应图7中所示的内部路径，线路3则如图8(a)所示。

对于同一个节点，U和V的插值函数存在细微差异，这是由于内部孔隙分布的各向异性所致。为更详细地观察插值函数的空间分布特征，从3D多孔规则化多孔单元中提取了三条线条：图7中的线路1和线路2，以及图8(a)中的线路3。图8(b)展示了沿这三条提取出的线条的位移插值函数分布情况。可以看出，预测的插值函数能够准确区分固体区域和孔隙区域，孔隙区域的插值函数值被设定为0。此外，尽管不同位移分量的插值函数沿线条的分布总体相似，但在某些区域仍存在差异。因此，本工作中提出的将每个位移分量的预测分开处理的神经网络架构是获得更精确计算结果的必要手段。

为分析本文提出的神经网络架构如何利用3D规则化多孔单元的几何信息来推导方程(4)中所给出的插值函数，我们提取并可视化了图4中虚线框内的MP-ANNs模块的热图。这三个大小均为151的不同核的权重，通过与前面模块中提取的几何特征进行内积运算，用来计算插值函数Nu、Nv和Nw的值。如图9所示，在核的第七位置，Nu和Nw的权重接近0，而Nv的权重接近1。这意味着Nv的核利用了输入特征中第七位置的信息，而Nu和Nw的核则忽略了这一信息。在第六十八位和第一百三十位也观察到了类似现象。这表明，尽管Nu、Nv和Nw这些插值函数较为相似，但各个核会选择代表三个不同方向的相关信息来计算插值函数，而非让其余两个插值函数基于某个单一插值函数的微小变化。本质上，这说明网络结构是根据不同的信息来计算不同方向的插值函数，而不仅仅是拟合一个单一的插值函数。

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

图9. 卷积核中参数的热图。

为进一步验证所提出的用于三维多孔材料计算的算法（即深度学习方法）的性能，以下部分将通过一系列三维示例及其计算结果来进行说明。比较中采用了有限元分析方法和FE2方法。FE2方法先获取刚度矩阵，然后使用与深度学习方法相同的后处理程序进行计算，以此减少比较时的变量数量。

3.2 具有周期性规则化多孔单元的多孔结构计算

如图10所示，使用经典的有限元分析方法、FE2方法以及本文提出的深度学习方法，计算了在均匀压力作用下具有周期性规则化多孔单元的三维多孔结构的力学响应。该结构由120个几何形状相同的规则化多孔单元组成，每个单元的尺寸为1×1×1米，杨氏模量为1吉帕，泊松比为0.3。该结构在X方向上受到均匀的正压作用。需要指出的是，为了在保持大尺度结构和微观孔隙的同时获得有限元分析结果，该案例仅包含120个规则化多孔单元，更大的案例将在后面介绍。如表1所示，使用有限元分析方法分析该结构大约需要6分钟，而FE2方法和深度学习方法分别仅需约3分钟和0.1秒。

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

图10. 使用有限元分析方法、FE2方法和深度学习方法计算的具有周期性规则化多孔单元的3D多孔结构的计算结果。(a) 该分析案例的参数设置。(b) 有限元分析方法、FE2方法和深度学习方法的计算时间。(c-f) 三种方法得到的位移场结果。

表1. 使用有限元分析方法、FE2方法和深度学习方法计算具有周期性规则化多孔单元的多孔结构的计算时间

| 方法 | 计算时间 |
|------|----------|
| 有限元分析方法 | 6分钟 |
| FE2方法 | 3分钟 |
| 深度学习方法 | 0.1秒 |

在传统的有限元分析方法中，必须使用高度精细的网格对整个结构进行离散化，这就导致了极大的自由度数量。相比之下，FE2方法和本文提出的深度学习方法则是通过将对应于规则化多孔单元的单元刚度矩阵组合起来，来执行宏观尺度的有限元分析。这两种方法的主要区别在于刚度矩阵的构建方式。对于FE2方法，首先需要在COMSOL中通过对规则化多孔单元进行尺度归一化的数值处理，获得每个规则化多孔单元的等效本构矩阵。具体而言，要对规则化多孔单元施加六个独立的宏观应变加载情况，这对应于三维弹性理论中的六个应变分量。对于每一个加载情况，都要求解微观边界值问题，并计算规则化多孔单元的体积平均应力响应。然后，利用给定的宏观应变与平均应力之间的关系，来构建规则化多孔单元的归一化本构矩阵。基于这个等效本构矩阵，就可以为宏观尺度的有限元分析构建相应的单元刚度矩阵。这一预处理过程为一个规则化多孔单元大约需要3分钟。相比之下，深度学习方法则是直接利用训练好的神经网络来预测同一个规则化多孔单元的形函数，并在0.1秒内构建出相应的单元刚度矩阵。由于FE2方法和深度学习方法之后的全局刚度矩阵构建以及宏观尺度有限元分析的步骤是相同的，因此总计算时间的差异主要来自于规则化多孔单元级刚度矩阵的计算。对于由相同规则化多孔单元组成的结构，这种规则化多孔单元级的计算只需进行一次。因此，FE2方法完成规则化多孔单元刚度评估总共需要大约3分钟，而深度学习方法则在0.1秒左右即可完成这一步骤，从而大幅降低了整体计算成本。

计算结果的位移等高线图如图10(c-e)所示，三种方法得到的位移分布均较为相似，且呈近似线性分布，这与理论预期相符。此外，在压缩变形情况下，有限元分析方法得到的位移值为4.28×10^-4米，深度学习方法得到的类似数值为4.16×10^-4米，而FE2方法得到的位移值为3.69×10^-4米。为更直观地显示三种方法在位移场计算结果上的差异，如图10(b)所示，提取了与X方向变化相关的边缘处的位移变化情况。可以看出，三种方法在X方向上的位移分布都呈现线性特征。深度学习方法与有限元分析方法的结果非常接近，相对误差约为2.5%，而FE2方法的结果与有限元分析方法相比差异较大，相对误差约为13.8%。

此外，如图11所示，还对上述案例中规则化多孔单元的计算得到的应力场进行了比较。需要注意的是，FE2方法将规则化多孔单元视为一个归一化单元，因此无法得到其内部的应力场。所以，应力场的比较仅限于有限元分析方法和深度学习方法。图11(a)展示了与图10中所示的规则化多孔单元的冯·米塞斯应力场分布相对应的等高线图。可以看出，深度学习方法计算得到的应力场与有限元分析方法的结果在规则化多孔单元的表面以及各切片内部都非常吻合。有限元分析方法和深度学习算法分别计算得到的规则化多孔单元的有效应力分别为41.0千帕和42.1千帕。两个数值都与均匀介质的理论值40 kPa相符，这体现了DL算法在计算RVE应力时的准确性。下载：下载高分辨率图像（620KB）下载：下载全尺寸图像图11. 使用FEM和所提出的DL模型得到的(a)完整RVE以及(b)两个随机选取的内部切面的应力场。值得注意的是，RVE的孔隙附近会出现应力集中现象。在孔隙周边，冯·米塞斯应力在X方向上降低，而在Y和Z方向上上升。图11(b)中的FEM和DL结果都反映了这一现象。在切片1中，孔隙附近的X方向上的冯·米塞斯应力低于平均值；而在切片2中，可以观察到孔隙中心处的YZ平面上存在高于平均值的冯·米塞斯应力区域。这种应力集中现象的重现与FEM结果一致，进一步证明了DL算法在处理RVE计算时的准确性和有效性，同时也表明DL方法可用于分析多孔结构内的内部损伤。为了更直观地比较FEM和DL算法得到的应力场差异，并更详细地观察孔隙附近的应力变化，从图11(b)中的位置线提取了应力分布，如图12所示。其中，沿Line 1_1的分布显示出在Y=16附近应力急剧下降，随后穿过孔隙结构时迅速上升至高于平均水平的值。出现这种趋势是因为Line 1_1先穿过第一个孔隙的X方向区域，然后又穿过第二个孔隙的Y方向区域。值得注意的是，FEM和DL都捕捉到了这一现象。在FEM方法中，由于网格离散化的原因，在应力快速上升阶段会出现剧烈的波动。在该位置，某个FEM单元位于XY方向上，因此它捕捉到的是Y方向的应力集中特征，而非X方向的应力下降，从而导致应力呈现剧烈的波动趋势，这一点在图11的等值线图中也能观察到。下载：下载高分辨率图像（531KB）下载：下载全尺寸图像图12. 3D多孔结构中(a) Line 1_1、(b) Line 2_1、(c) Line 1_2、(d) Line 2_2的应力分布。从四条提取出的线上的应力分布可以看出，DL方法能够得到与FEM相近的力学响应，并且能够准确捕捉孔隙附近的应力变化。这一结果充分证明了DL方法在计算三维多孔材料的力学响应时，能够准确捕捉足够的局部细节。3.3. 具有非周期性RVE的多孔结构计算对于涉及多个具有不同形态的RVE的计算，DL方法具有显著优势，因为它基于不同RVE的预测插值函数。通过GPU处理实现并行化的预测过程，其计算时间约为0.1秒。而对于FE2方法来说，多个RVE意味着需要在RVE尺度上进行多次力学计算以获得各自的刚度矩阵，这大大增加了计算成本。为体现这一优势，对前一节中的示例进行了修改，加入了多个RVE。如图13(a)所示，在均匀压力作用下，对含有三种不同内部结构的RVE的多孔结构进行了力学响应计算。具体而言，在X从0到4米的范围内，使用图中最左侧RVE的孔隙结构，其杨氏模量为0.5 GPa，泊松比为0.3；在X从4到7米的范围内，使用中间RVE的孔隙结构，其杨氏模量为1 GPa，泊松比也为0.3；在X从7到10米的范围内，使用最右侧RVE的孔隙结构，其杨氏模量为1.5 GPa，泊松比仍为0.3。该示例中FEM和DL方法的计算过程与前节类似，没有明显差异。但对于FE2方法来说，由于存在三种不同的RVE，需要重复三次RVE刚度矩阵的计算。如表2所示，相应的计算时间大致与不同RVE的数量成线性关系。因此，与FEM和DL方法相比，FE2方法解决上述力学问题所需的时间更长。下载：下载高分辨率图像（405KB）下载：下载全尺寸图像图13. 通过FEM、FE2和DL方法得到的具有非周期性RVE的3D多孔结构的计算结果。(a)研究案例的参数设置。(b)FEM、FE2和DL的计算时间。(c-f)FEM、FE2和DL的位移场结果。表2. 使用FEM、FE2和DL方法计算具有非周期性RVE的多孔结构的计算时间。方法计算时间FEM6分钟FE29分钟DL0.1秒图13(c-e)中的计算结果表明，DL方法计算得到的变形与FEM算法相比误差约为2.4%，而FE2方法的误差则约为8.7%。因此，DL方法得到的结果更为精确。此外，图13(b)还展示了与X方向变化相关的边缘处的位移变化情况。可以看出，DL方法能更准确地再现FEM的计算结果。DL方法和FEM方法得到的结果不仅包括更大尺度下的几何力学响应，还包括RVE尺度下的力学响应。这些方法计算得到的应力场等值线图如图14所示。第一组、第三组和第四组以及第五组和第六组切片分别是从图13(a)中的矩形框从左到右提取的。在每一切片上，DL方法计算得到的应力场都与FEM的结果一致，能够捕捉到诸如孔隙位置处垂直于压缩应力面的应力集中现象，以及穿过孔隙结构时平行于压缩应力方向的应力下降现象。这进一步证明了DL在计算RVE内应力分布方面的准确性和高效性。下载：下载高分辨率图像（403KB）下载：下载全尺寸图像图14. 3D多孔结构中RVE的应力场，其中切片1和2、切片3和4以及切片5和6分别来自三种不同的RVE。同样值得注意的是，三种类型RVE内的平均应力相对接近。这是因为尽管它们的内部结构和杨氏模量不同，但都承受着相同大小的均匀压力。同时，上述示例也展示了DL算法在节省时间方面的优势：即便RVE的结构越来越复杂，DL算法的计算时间也基本保持不变。这是因为添加多个RVE结构只需通过神经网络预测它们的插值函数，而由于DL模型计算的并行性，这一步骤所需的计算时间极短（约0.001秒）。此外，该示例还表明DL方法适用于多孔材料的各种内部孔隙。对于不同的孔隙结构，只需进行简单的预测处理，无需重新训练网络。3.4. 具有多种RVE的超大尺度多孔结构计算最后，为了进一步凸显DL算法的高效性，将几何结构的规模扩大，以便包含更多的RVE。由于计算区域增大，FEM单元的数量超出了可直接计算的极限。因此，在这种情况下，仅比较DL算法和FE2算法得到的结果。具体示例如下。如图15(a)所示，该多尺度结构包含1200个形态相同但杨氏模量和泊松比不同的RVE，其规模比前几节中的示例大了一个数量级。此外，该结构还包含大量半径约为65毫米的小孔隙，而整个域的长度则为60米。宏观结构与其孔隙特征之间的尺度比超过了两个数量级。结构的杨氏模量和泊松比如图所示会发生变化。在这种情况下，如表3所示，DL算法仍然不需要大量的计算时间。所提出的方法只需通过组合预测的插值函数来反复构建刚度矩阵。这一过程相对简单且计算效率高，完成分析仅需大约10秒。相比之下，在FE2算法中，无法仅通过修改材料属性来更新刚度矩阵，尤其是当某个RVE的泊松比发生变化时。相反，必须进行多次RVE尺度的力学计算才能确定相应的刚度矩阵。因此，FE2方法在RVE尺度计算上需要大量的计算时间。在这个示例中，利用对称性仅分析了结构的一半，但仍涉及600个不同的RVE。对于FE2方法而言，每个RVE都需要进行独立的RVE尺度数值模拟，以获得其对应的本构矩阵和元素刚度矩阵。由于单个RVE的计算时间约为3分钟，因此总的预处理时间约为30小时。下载：下载高分辨率图像（347KB）下载：下载全尺寸图像图15. 通过DL和FE2方法得到的具有不同弹性常数的RVE的多尺度3D多孔结构的计算结果。(a)研究案例的参数设置。(b)DL和FE2的计算时间。(c-d)DL和FE2的位移场结果。表3. 使用FE2和DL方法计算具有多种RVE的超大尺度多孔结构的计算时间。方法计算时间FE230小时DL10秒如图15(b-d)所示，该示例的位移场计算结果表明，DL方法得到的数值略高于FE2方法的结果。这一结果与前面章节的结果一致，进一步证明了DL方法依然能够快速、准确地计算具有不同材料属性的复杂三维多孔结构的力学响应。此外，X方向边缘处的位移曲线呈单调递减趋势，且斜率逐渐变小。这是因为杨氏模量随X呈线性增加。不过，由于泊松比也会随X发生变化，这些位移曲线并不严格遵循由其二阶导数决定的固定下降速率。由于FE2方法无法获得局部RVE尺度的计算结果，因此仅提取了DL算法的应力场计算结果，如图16所示。在切片1和切片3中，可以观察到应力集中和应力水平下降的现象。在切片2中，由于该切片位于垂直于均匀压力方向的平面上，因此只能观察到孔隙附近的应力集中现象。此外，在切片2的固体区域中，还可以观察到局部最小应力值，这说明该点靠近该切片上不存在的孔隙，具体位于X方向上。在切片1的固体区域中，也可以观察到远离该切片上孔隙处的最大应力点，这表明该点靠近该切片上不存在的孔隙，具体位于Z方向上。除了孔隙附近区域外，其他位置的冯·米塞斯应力值约为39 kPa。该图表明，DL方法在处理具有不同材料属性的3D超大尺度多孔结构的力学响应计算时，依然能够快速得到合理的数值结果。下载：下载高分辨率图像（380KB）下载：下载全尺寸图像图16. 3D多尺度多孔结构中某个随机RVE的应力场，其中切片1-3是同一RVE中的三个正交切片。总之，验证过程分为四个部分，旨在说明DL方法在处理三维多孔结构时的高效性和准确性。在第一部分中，验证了DL方法的核心——即预测多孔RVE内插值函数的准确性。一个三维神经网络模型在包含10,000组具有1,000种不同孔隙结构的RVE位移场数据集上进行了训练和验证。在第二部分中，展示了深度学习方法在计算三维多孔结构力学响应方面的时间成本优势，以及其在计算结果方面相较于FE2算法更接近有限元分析的精度。在第三部分中，通过将几何结构改为包含三种类型的RVE，体现了深度学习方法在多RVE计算中的优势。在第四部分中，为了突出深度学习算法在计算时间成本和精度方面的优势，选取了一个由1,200个具有不同弹性常数的RVE组成的多尺度多孔结构，其中许多孔隙半径比整个结构的尺寸小两个数量级。通过这四个部分，本节逐步证明了深度学习算法的强大能力：它能够快速准确地计算三维多孔结构的力学响应，同时获得宏观和RVE尺度上的结果。该算法显示出在实际工程应用中大幅加速大规模有限元计算的能力，例如在岩土工程中的变形估算[44,45]。具体而言，地质构造中存在大量的孔隙-裂隙结构，其力学参数表现出明显的空间异质性[63]，如泊松比和杨氏模量。Hu等人[64]提出了一种基于CNN-LSTM的深度学习模型来预测井测数据，准确捕捉了岩石力学性质的变化规律。在此基础上，对于由不同RVE构成的复杂地质构造的有限元分析，所提出的方法能够在准确捕捉不同地下区域的应力分布模式的同时，显著提高计算效率，从而为地质开发实践提供可靠依据。

4. 结论
为了实现对多孔材料进行快速准确的有限元分析，本文提出了一种基于深度学习的新型方法，该方法包括宏观尺度和RVE尺度上的计算。与传统均质化方法不同，该方法的核心思想是通过隐式插值函数场来描述RVE的内部孔隙结构，并通过这些插值函数建立宏观尺度与RVE尺度之间的耦合关系。这样一来，不仅可以节省时间，还能获得每个RVE的内部应力分布。具体而言，利用插值函数构建刚度矩阵以计算宏观位移场，随后可从宏观节点位移出发计算每个RVE内的局部位移和应力场。这些隐式插值函数是由由多并行DeepONets和ResNets构成的新深度学习框架预测得到的。并行运算模块能够高效预测三维插值场，而ResNets则可稳定这些插值场的一阶导数，从而提高所得应力场的精度。

验证实验结果表明，与传统的有限元分析和FE2方法相比，深度学习方法的计算速度提升了三到四个数量级，其结果与有限元分析的吻合度更高，位移误差保持在3%以内。此外，深度学习方法还能捕捉到RVE内部的应力场和位移场分布，尤其是孔隙附近的应力变化，这体现了其相对于FE2方法的优越性。最后，所提出的深度学习算法能够在短短10秒内计算出含有大量孔隙的多孔材料的力学响应，这些孔隙的半径比整个结构尺寸小两个数量级以上。相比之下，由于所需元素数量极大，有限元分析方法受到硬件限制而无法完成计算。尽管仍有一些方面需要进一步研究以确保更广泛的应用，但该方法实现了对三维多孔结构的高效有限元分析。未来的工作将扩大数据集范围，纳入不同大小的RVE样本以及节点不一定完全位于固体区域内的更复杂案例。此外，还将改进神经网络架构，以处理具有不同大小和几何形状的更复杂多孔结构。经过进一步优化，该方法还可应用于先进材料的设计与优化[[65], [66], [67]]，使工程师能够以前所未有的效率对微观结构设计进行迭代优化。预计该方法将成为一种功能多样的预测工具，能够支持多种多尺度工程应用。

CRediT作者贡献说明
Hao Yu：概念构思、方法设计、实验研究、资源获取、初稿撰写、审稿与编辑、项目管理、结果验证。JiaXuan Chen：初稿撰写、审稿与编辑、结果可视化、软件应用、形式化分析。HanXi Xu：实验研究、方法设计、数据整理。KuiCheng Duan：结果可视化、审稿与编辑、结果验证。Quan Wang：概念构思、审稿与编辑、资金筹集、实验研究。SiWei Meng：方法设计、初稿撰写、资源获取、审稿与编辑、资金筹集。He Liu：监督指导、概念构思。HengAn Wu：资金筹集、监督指导、项目管理。