一种基于转动惯量的多孔聚合物梁表面均匀化方法 向仁强、 李硕、 张明、 李莉

《Polymers》：A Moment-of-Inertia-Dependent Surface Homogenization Method for Porous Polymer Beams Renqiang Xiang, Shuo Li, Ming Zhang and Li Li

　　摘要：在多孔聚合物梁中观察到了明显的尺寸依赖性弯曲响应，尤其是在其厚度接近晶格常数尺度时。然而，尺寸依赖性与微观结构之间的关系仍不清楚。由于微观结构的复杂性，直接数值模拟在计算上非常昂贵，而忽略表面效应的经典多尺度方法所得结果与实际行为存在显著偏差。本研究开发了一种考虑表面驱动惯性矩的多孔聚合物梁等效模型，通过引入表面强度因子和表面厚度来捕捉尺寸依赖的杨氏模量。然后，基于由惯性矩依赖的表面均匀化方法生成的离线数据集，建立了一个在线预测框架，用于尺寸依赖的弯曲响应。所提出的框架在准确性和计算效率方面都得到了评估。结果表明，经典多尺度均匀化方法可能产生高达1108%的相对误差，而表面均匀化方法则将相对误差控制在4%以下。此外，与直接数值模拟相比，计算成本大幅降低。这项工作不仅揭示了超材料梁尺寸依赖行为的潜在惯性矩依赖表面机制，还为它们的结构设计和性能预测提供了一个准确高效的工具。

1. 引言
超材料是人工设计的材料，表现出非常规且常常非凡的特性[1,2,3]，包括优异的机械、声学和热学性能[4,5,6,7]，近年来由于制造技术的快速进步，特别是3D打印的广泛应用，引起了越来越多的关注。这些发展为超材料在航空航天、土木工程和生物医学系统等众多领域的应用和深入研究提供了可能[8,9,10,11]。在各种结构形式中，超材料梁具有可定制的刚度和可编程的机械响应，在智能系统中的振动隔离和能量收集方面展示了显著优势[12,13]，并且在工程设计中承重组件的设计中也具有广泛的适用性[14,15,16]。现有研究表明超材料中存在尺寸依赖效应[17,18,19]，这对于评估超材料组件的变形行为至关重要。当超材料结构受到剪切或弯曲载荷时，尺寸效应变得更为明显[20,21,22]，且梁结构沿厚度方向的晶格单元数量越少，在弯曲变形过程中观察到的尺寸依赖性越显著[23,24]。这种行为与拉伸变形下的尺寸依赖特性不同[25,26]，因为弯曲刚度对结构尺寸的敏感度明显高于拉伸刚度[27]。然而，在分析具有复杂微观结构的超材料梁的变形时，使用精细网格离散化并用高保真有限元方法求解需要大量计算资源[28,29]，甚至可能导致收敛问题或模型无法求解。因此，一种高效准确的均匀化方法对于预测超材料梁的变形行为将非常有益。
均匀化是一种从微观代表性体积元素提取机械性能的数值过程，该元素表征了超材料的整体行为。均匀化后，可以利用等效的机械性能在宏观尺度上模拟整个晶格，从而生成更简单的网格，显著减少计算时间。最初，均匀化理论是为了加速复合材料的模拟而引入的，现在它被广泛用于预测超材料结构的有效机械行为[30,31,32]。迄今为止已经开发了许多基于经典力学的均匀化方法。其中，一阶渐近展开方法是一个例子[33]，它基于多尺度展开的概念，将位移场展开为慢变量和快变量的级数，并通过解决一组微观问题获得宏观本构关系。另一种重要的方法是变分均匀化方法[34,35]，它基于最小势能原理，系统地将变分公式应用于多尺度问题。这两种方法都是基于代表性体积元素（RVE）框架[36]进行的，当处理具有周期性微观结构的介质时，RVE必须远小于特征宏观长度尺度，以便在整个结构中将微观和宏观行为视为两个独立的但耦合的系统，在微观尺度上将宏观变量近似为常数，从而实现变量的分离[37,38]。此外，经常使用半解析方法（如自洽方法和平均场理论）来确定非周期性多相材料的有效性能。这些方法引入了等效介质或嵌入式微观结构模型，通过统计平均原理建立了微观尺度和宏观尺度之间的联系[39,40]。在实际的数值计算中，基于RVE的有限元方法（FEM）已成为工程应用中最常用的均匀化工具。这种方法的关键在于正确定义RVE的边界值问题[41,42,43]。RVE方法利用有限元技术模拟单个周期性单元格，应用周期性边界条件（PBC）来模拟无限结构的响应行为[44,45]，从而提取宏观等效刚度张量或本构关系[46]。这些经典均匀化方法在处理周期性结构和介观问题时表现良好，但它们都假设材料响应是局部连续的，忽略了由相似的微观和宏观长度尺度引起的非局部效应和边界层影响。因此，当结构尺度接近或小于微观尺度，或者材料在弯曲过程中表现出显著的尺寸依赖性行为（如软化或硬化）时，经典方法往往无法准确预测其响应[47,48]。
随着实验和模拟研究的进步，许多研究表明，当微观结构的特征尺寸接近结构的整体几何尺度时，传统的基于尺度分离的均匀化方法无法准确预测超材料的机械行为[49,50,51]。为了解决这个问题，研究人员逐渐开发了一系列基于非经典连续介质理论的均匀化方法，旨在克服经典模型的局限性，更真实地捕捉超材料中的非局部响应、边界层效应和微宏观耦合特性。Gurtin-Murdoch表面弹性理论是第一个描述微观尺度材料机械行为的数学框架，通过在体材料表面引入独立的“表面能”或“表面应力”来扩展经典弹性理论，以表征微观和纳米尺度上的机械行为[52,53]。另一方面，应变梯度理论引入了应变的更高阶导数，建立了嵌入内在长度尺度参数的非局部本构关系，从而克服了经典柯西模型“零特征长度”的限制[54,55]。此外，Mindlin和Tiersten建立的微极理论[56]以及Fleck提出的应变梯度塑性理论[57]为这类模型提供了重要基础。基于这些理论的均匀化方法可以更准确地描述当结构厚度或特征长度减小到微观结构尺度时出现的尺寸效应[58]。另一类重要的方法基于Cosserat或微极连续介质理论，认为材料点不仅具有平移自由度，还具有独立的旋转自由度，引入了额外的耦合力矩和旋转刚度参数。这些理论特别适合描述在单元格内发生显著局部旋转、空洞、铰链或剪切滞后现象的结构，例如蜂窝状或以剪切为主导的超材料[59,60,61]。尽管基于非经典力学的均匀化理论可以在一定程度上捕捉尺寸效应和非局部响应，但当应用于具有复杂微观结构的超材料（如研究晶格微观结构和几何形状如何影响3D打印聚合物材料的机械行为）时，填充几何形状（立方体、螺旋形、线性等）和填充密度会对结构机械性能产生显著影响[62]。在这种情况下，准确确定模型参数（如内部长度尺度和高阶弹性模量）通常具有挑战性，这些参数的选择会对模型的预测精度产生显著影响。
在微观和纳米尺度上，已经广泛研究了可能由表面应力、表面弹性和残余表面应力引起的表面效应[63,64]，而在宏观尺度上，超材料中的表面效应是由微观特征与整体结构的变形相互作用引起的[65]。在复合材料弯曲性能的研究中，许多理论和实验工作提供了有价值的见解。例如，对功能梯度碳纳米管增强复合夹层梁的弯曲行为的研究表明，材料梯度和结构配置可以显著影响梁的弯曲刚度和变形特性[66]。同时，对碳纳米管增强纺织品基复合材料的实验和分析研究也表明，纳米级增强可以有效地提高结构的弯曲性能和承载能力[67]。此外，以前关于超材料结构中表面效应的研究主要集中在宏观拉伸变形上[68,69]，相对较少的研究关注表面效应对宏观弯曲变形下超材料梁的影响，这大大限制了超材料梁在工程应用中的结构设计和性能预测。显然，开发一种考虑表面效应的均匀化方法对于准确预测超材料梁的弯曲变形行为具有重要意义。
基于上述文献，我们的研究主要解决以下两个问题：首先，超材料梁在弯曲载荷下抵抗变形的能力受到表面效应的影响，但在不成立尺度分离假设的情况下，这种机制仍然不清楚。因此，有必要研究多孔聚合物梁在弯曲变形过程中表面效应的潜在机制，这对它们的性能预测和结构设计具有重大意义。其次，由于多尺度分析方法固有的无穷小假设不适用于表现出显著表面效应的多孔聚合物梁，这会导致计算上的重大不准确性。此外，构成多孔聚合物梁的晶格具有复杂的微观结构特征，因此进行直接数值模拟会导致计算效率低和成本高。为了解决这些挑战，我们旨在开发一种高效准确的均匀化方法，该方法在弯曲载荷下考虑表面效应，从而促进超材料在航空航天、土木工程和生物医学领域的广泛应用。
本文的结构如下：第2节介绍了多孔聚合物梁中观察到的显著尺寸依赖性弯曲变形，并使用高保真有限元分析研究了其潜在机制。第3节开发了一种考虑表面驱动惯性矩的多孔聚合物梁等效模型，并基于表面驱动均匀化方法生成的离线数据集建立了用于尺寸依赖弯曲响应的在线预测框架。第4节详细讨论了表面效应与微观结构特性之间的关系，并评估了所提出均匀化方法的准确性和计算效率。最后，第5节总结了研究的主要结论。图1展示了研究的整体示意图，说明了考虑表面效应的均匀化方法相对于忽略表面效应的高保真有限元方法和多尺度分析方法的优势。通过使用这种均匀化方法构建与晶格相关的数据库，可以准确高效地预测多孔聚合物梁的机械变形行为，这对于具有复杂微观结构的多孔聚合物梁的结构设计和性能预测具有重要意义。

2. 问题陈述
本研究的对象是由重复晶格组成的细长多孔聚合物梁，这些晶格具有复杂的微观结构，包括四种类型的晶格：Diamond、Primitive、Fischer-Koch S和Split P，它们可以通过偏移三周期最小表面（TPMS）转换为固体结构。通过调整孔隙率、晶格类型和晶格常数，可以调节它们的性能，如有效弹性模量和比刚度。晶格常数表示为a，孔隙率表示为p。在由TPMS晶格组成的多孔聚合物梁中，表面效应表现为：图1a展示了有效模量随梁厚变化的趋势，该趋势是通过将完整的微观结构放入高保真度有限元模拟中获得的。本研究中使用的材料是Photosensitive Resin 9400，其杨氏模量为E_s，密度为ρ，泊松比为μ。本研究用于探讨表面效应对超材料梁整体刚度影响的模型均基于线性弹性材料行为的假设。随着宏观厚度H的增加，有效弹性模量也增加。当H=1μm时，有效模量是E_s的7.74倍，且这一趋势随着厚度的进一步增加而持续。同时，体弹性模量相对于H=0μm的情况增加了12倍。这一现象对受弯曲载荷作用的多孔聚合物梁的变形响应产生了显著影响。然而，其背后的机制尚不清楚，这限制了其设计和预测的准确性。经典力学理论认为，在抵抗弯曲变形时，远离中性轴的区域贡献更大。然而，我们的研究发现，当受到载荷作用时，多孔聚合物梁的表面区域对抵抗弯曲变形的贡献小于核心区域，这由超材料梁的微观结构决定。图2展示了多孔聚合物梁中表面效应的本质机制。图2b显示了不同厚度H的多孔聚合物梁的局部视图，其中n表示不同的层，同一层内的晶格经历相同的变形。如图2c所示，当H=2μm时，在变形的超材料梁中，第n层的应力主要集中在中心区域，表明表面区域对抵抗弯曲变形的贡献很小。而对于H=1μm的梁，变形后核心层的应力主要集中在表面附近。同样，如图2d所示，当H=3μm时，该区域的应力分布与H=2μm的情况类似，而当H=4μm时，应力主要集中在表面附近，这由晶格的微观结构决定。图2. 多孔聚合物梁中的表面效应机制。（a）不同厚度多孔聚合物梁的表面效应。（b）不同厚度多孔聚合物梁的局部视图。（c）不同梁厚H下区域n的归一化应力分布。（d）不同梁厚H下区域n的归一化应力分布。图中的颜色条从0到1表示归一化的von Mises应力。对于每种结构，颜色越接近1表示该位置的应力越大，橙色和绿色框表示在不同厚度条件下的相同空间坐标区域。图3展示了基于表面效应的细长多孔聚合物梁的建模，由于超材料梁的厚度方向上有两个自由表面，位于这些自由表面的晶格在面对自由表面的边界上缺乏相互作用的结构，这与中间核心区域的晶格形成对比。核心区域的晶格被四周的微观结构包围，在弯曲变形过程中可以将载荷传递给相邻晶格，而表面区域的晶格受到的约束较少，这种约束不足导致表面区域抵抗弯曲变形的能力较弱。这种现象与纳米尺度的表面效应不同，后者源于表面原子与体原子的键合邻居较少或配位数较低[70]。在具有复杂微观结构的多孔聚合物梁中，表面效应主要因为表面区域在受载荷时受到约束不足而产生。对于TPMS结构，几何形状内的曲率变化也可能影响局部机械响应；然而，在目前的模型中，这种效应被隐含在所采用的几何表示中，没有作为独立的变量进行分析。此外，与纯拉伸载荷下的变形不同，当受到弯曲载荷时，位于远离中性轴的超材料梁的表面区域起着更重要的作用，因此表面效应的影响变得更加明显。特别是当超材料梁相对较薄时，晶格尺寸与宏观厚度相当，自由表面引起的约束不足不能被忽略，因为它显著影响梁抵抗弯曲变形的能力。因此，在计算受弯曲载荷作用的多孔聚合物梁的弯曲阻力时，需要考虑微观结构的表面弹性。然而，由于微观结构表面受到约束不足，无穷小应变假设被违反，导致一些基于经典力学的多尺度分析方法不再适用。因此，有必要开发一种基于经典力学的均质化方法，该方法考虑表面效应，以提高预测的准确性。在弯曲载荷下，较弱的表面区域可能与局部剪切或胞元变形模式相关联。然而，本研究的重点是整体结构刚度和宏观表面效应，由于模型的限制，未进行局部变形模式的系统分析。还应注意，用于增材制造的聚合物通常表现出粘弹性行为，这意味着它们的机械性能取决于时间、加载率和温度，这可能会进一步影响自由表面的变形机制。先前的研究表明，3D打印的聚合物样品在持续载荷下可以表现出明显的粘弹性蠕变和非线性响应，由于逐层制造过程和材料链的方向性沉积，它们还表现出与打印方向相关的各向异性机械行为[71,72]。这些现象表明，粘弹性和打印引起的各向异性可能会改变自由表面附近的局部刚度和应变分布，从而以纯线性弹性均质化无法捕捉的方式调节表面效应。然而，本研究基于理想的线性弹性材料，没有考虑聚合物粘弹性和打印引起的各向异性对表面效应的影响。未来的工作可能结合粘弹性本构模型和各向异性机械响应，以更准确地描述3D打印聚合物超材料的时间和方向依赖性机械行为。图3. 细长超材料结构在弯曲变形下的表面效应。根据上述描述，我们的研究重点解决以下两个问题：首先，从科学角度来看，多孔聚合物梁在非尺度分离条件下的表面效应的本质机制尚不清楚，这限制了性能预测和结构设计的准确性。因此，揭示多孔聚合物梁中表面效应的本质机制非常重要。其次，从技术角度来看，使用DNS计算多孔聚合物梁的机械性能需要大量的计算资源，而基于尺度分离假设的经典均质化方法的预测误差相对较大，特别是当构成超材料梁的晶格尺寸接近梁的宏观厚度时。因此，需要开发一种高效且准确的均质化方法，该方法考虑超材料梁的表面效应。通过解决上述两个问题，可以获得表面效应对多孔聚合物梁弯曲阻力影响的特点，并通过有限元分析推导出梁的尺寸依赖性参数，从而揭示由梁的微观结构引起的表面效应对其机械性能的影响。3. 基于表面驱动均质化的在线预测框架在本节中，建立了一个基于表面驱动均质化方法生成的离线数据集的在线预测框架，用于尺寸依赖性弯曲响应的计算。3.1. 考虑表面驱动惯性矩的多孔聚合物梁等效模型基于经典力学的多尺度分析方法在计算多孔聚合物梁的弯曲变形时经常忽略表面效应，这通常会导致较大的计算相对误差。当考虑表面效应时，由于表面区域远离中性轴，并且根据弯曲刚度与截面惯性矩之间的关系，表面区域对梁抵抗弯曲变形的贡献更大。因此，在梁的弯曲行为中不能忽略表面效应。在本节中，为了准确表征受弯曲载荷作用的超材料梁结构的有效刚性，引入了惯性矩，并开发了一种考虑多孔聚合物梁表面效应的均质化方法。此外，建立了一个与结构尺度相关的归一化弯曲刚度表达式，反映了表面效应对多孔聚合物梁机械响应的影响。为了描述表面效应对梁弯曲刚度的影响，建立了一个尺寸依赖的有效弹性模量模型，如图4所示。由于多孔超材料梁的表面区域缺乏一些相互作用的结构，因此受弯曲载荷作用的表面区域和核心区域的弹性张量不同。在图中，将梁的多孔内部结构转换为固体结构，同时考虑了表面效应。因此，沿宏观厚度方向，梁可以分为两个区域：由深蓝色表示的核心区域和由浅蓝色表示的表面区域。受表面影响的区域可以表示为σ_s(z)，其中z是笛卡尔坐标系中的厚度方向，核心区域对应于σ_c。假设表面区域的弹性模量为E_s，而核心区域（即体模量）的弹性模量为E_c。本研究中引入的参数并不严格代表表面区域的真实机械属性。相反，它是一个用于表征表面效应影响的有效参数，其值通过后续的均质化程序确定。它反映了表面区域在微观结构中对整体机械响应的贡献，可以视为表面效应的有效表示，从而有效地捕捉它们对宏观机械行为的影响。图4. 相当于超材料梁的固体梁的体模量和表面模量。此时，梁不同区域的有效弹性模量可以表示为方程（1）：(1) σ_s(z) = E_s * (1 - p) / (1 + p * μ)，其中E_s是表面区域的弹性模量，E_c是体模量，H是梁的厚度，l是表面区域的总厚度，它取决于晶格类型、孔隙率和晶格常数，但与材料类型无关。根据图3和方程（1），整个梁的有效弹性模量表示为Eeff = E_s * (1 - p) / (1 + p * μ)。然后，梁的有效弯曲刚度可以表示为σ_eff = Eeff * H / (1 + p * μ)。σ_s、Eeff和H之间的关系可以通过方程（2）表示，其中I_y是关于y轴的等效固体梁的总惯性矩，由方程（3）给出：(2) I_y = E_c * H^3 / (1 + p * μ)，(3) I_s = E_s * (1 - p) ^3 / (1 + p * μ^2)。表面区域的弹性模量与体弹性模量之间的关系假设为方程（6）：(6) σ_s / E_c = k * (1 - p)。然后，有效弹性模量与体弹性模量作为厚度H的函数可以表示为：(7) σ_eff = E_s * (1 - p) * (1 + p * μ) ^ (1 - 2p)，其中l是表面厚度。请注意，k和l受到晶格微观结构特征的影响，仅依赖于晶格的结构特性。弯曲刚度可以用来评估梁抵抗弯曲变形的能力。在考虑了表面效应并进行均匀化处理后，有效弯曲刚度为...（此处数值缺失），而通过基于经典力学的多尺度分析得到的均匀化弯曲刚度为...（此处数值缺失）。这两个数值的比值为：（8）这个无量纲比率描述了表面效应对梁有效弯曲刚度的影响。为了理解弯曲过程中表面效应的机制，需要校准两个内在参数k和l。然而，将完全解析的微观结构纳入直接高保真模拟（DHS）需要大量的计算资源和时间。因此，基于代表性体积元素（RVEs）的表面驱动均匀化方法是必要的。

3.2 构建在线预测框架
表面驱动的在线预测框架需要构建一个与晶格尺寸相关的离散数据集，从而能够在线预测具有不同厚度的多孔聚合物梁的机械性能。均匀化方法的一个关键方面是选择具有微观尺度周期性结构特征的代表性体积元素（RVE），该元素来自宏观上不均匀的结构。通过对RVE进行适当的计算和分析，可以获得其在微观尺度上的机械行为，从而准确预测整个结构的宏观性能。在本研究中，宏观不均匀结构是一个复杂的梁，需要从中提取一个代表性体积元素（RVE）来计算非均匀材料的本构参数。基于经典力学和多尺度分析的均匀化方法没有考虑表面效应，用于获得体弹性模量。这种方法依赖于无限小应变假设，即结构的微观尺度必须远小于宏观尺度[30]。因此，为了满足这一要求，必须在感兴趣的晶格边界施加周期性边界条件。由于本工作关注的是梁在-x平面内的弯曲行为，只需确保晶格在x和z方向上的边界满足周期性边界条件即可，如图5a所示。在弯曲过程中，晶格在x方向上主要受到拉力和压缩力的共同作用，而在z方向上施加周期性边界条件只是为了消除厚度方向上的表面效应影响。

图5. 两种不同均匀化方法的边界值问题。(a) 基于晶格的经典均匀化方法。(b) 基于全厚度RVE的表面效应均匀化方法。(c) 均匀化结构的变形。(d) RVE大小（即其所包含的单元格数量）对结果的影响：分析了由1个、2个和3个单元格组成的RVE的效果，x方向施加周期性边界条件，图中的色条表示von Mises应力，单位为Pa。因此，在计算x方向的体模量时，施加一个具有预定位移扰动的动力均匀边界条件，同时允许z方向无牵引力。同样，可以使用周期性边界条件来满足Hill–Mandel条件[69,73]。设RVE体积为V，其中A表示基质（固体）子域，Vp表示孔隙（空洞）子域。宏观（体积平均）应变和应力定义如下：（9）（10）其中δA是基质相的指示函数（在A中为1，在Vp中为0），孔隙相中的应力消失。有效本构关系为：（11）其中Eeff是从多尺度分析中得到的有效弹性张量。对于沿x方向的单轴均匀化，我们提取Eeff。然而，在微观尺度上，材料的微观结构本质上是非均匀的，当受到外部载荷作用时，内部的载荷传递路径变得非常复杂。此外，由于三周期最小表面（TPMS）独特的曲面几何形状，表面区域之间的连接倾向于表现出某些结构缺陷，这导致有效弹性模量降低。另外，随着梁的总厚度变薄，这些表面区域引起的表面效应变得越来越显著。在这种情况下，结构配置不再满足经典多尺度理论所需的无限小假设。因此，理论多尺度分析不再适合预测具有明显表面效应的微观结构的机械行为。为了提高预测材料宏观机械性能的准确性，提出了一种考虑表面效应的均匀化方法，以确保与实际微观结构的一致性。

如图5b所示，为了研究厚度对表面效应的影响，首先需要确保z方向的厚度保持不确定。随后，选择一个具有全厚度的代表性体积元素（RVE），并在x方向的两个边界施加周期性边界条件，以模拟一个长度尺度远大于其厚度的梁模型。这种配置还确保了边界满足位移连续性和牵引连续性条件。同样，由于分析仅关注梁在-x平面内的弯曲变形，y方向的边界被视为自由表面。梁的厚度沿z方向变化；因此，自由表面以这种方式定义，同时在梁长度方向上的RVE施加周期性边界条件。这使RVE能够有效捕获厚度变化对局部刚度和表面效应的影响。尽管RVE代表一个局部单元，但实际梁结构的全局边界条件不会改变自由表面沿z方向的物理意义或局部机械响应。基于这些条件，可以获得受到弯曲变形的全厚度RVE的有效弹性模量，用于构建与不同晶格配置相关的离线数据库。由于沿x方向施加周期性边界条件以模拟无限长的超材料梁，并且考虑到厚度尺度远小于长度尺度，因此采用Euler–Bernoulli梁理论进行分析。晶格在x方向边界受到方程（12）中描述的应变场的作用，从而引起弯曲变形，如图5c所示。此外，为了验证本研究中选择的RVE的合理性，分别构建了由1个、2个和3个单元格组成的RVE，并在x方向施加周期性边界条件进行比较，如图5d所示，孔隙率为0.9。结果表明，在相同的应变条件下，不同RVE大小的单元格的变形模式和应力分布基本一致。进一步计算表明，从三个代表性体积元素（RVEs）得到的有效弹性模量是相同的，这表明RVE中包含的单元格数量对研究结果没有影响。这证实了本研究中选择的RVE大小具有足够的代表性及收敛性，可以可靠地预测结构的宏观刚度以及表面效应的趋势。

（12）其中R是曲率半径。然后，根据晶格关于y轴的惯性矩公式（13），可以如下确定惯性矩I：（13）然而，当考虑晶格微观结构时，惯性矩I并不容易获得。根据图5b，具有多孔微观结构的TPMS相当于一个均匀模型，其截面近似为矩形截面，如图5c所示。此时，可以容易地得到关于y轴的均匀化惯性矩，表示为方程（14）：（14）其中Iy是关于y轴的均匀化惯性矩，W和H分别是截面的宽度和高度。在这个阶段，可以根据方程（15）得到x边界处的弯曲力矩M与均匀化有效弹性模量之间的关系：（15）通过使用表面驱动的均匀化方法，可以获得不同厚度多孔聚合物梁的有效弹性模量。值得注意的是，与传统多尺度分析不同，这种考虑表面效应的均匀化方法考虑了TPMS微观结构表面在厚度方向上的介质不连续性。因此，在计算Iy时，需要确保z方向上的表面是自由的，这与计算I时的条件正好相反。计算Iy时不考虑表面效应，并且仍然满足无限小假设，因此也必须在z方向上施加周期性边界条件。对于不同类型的晶格，可以使用上述均匀化方法获得表面厚度l和表面强度因子k。在本研究中，参数k和l共同影响宏观梁的整体刚度和表面效应；然而，它们在物理上是独立的。这两个参数通过数值校准来确定，以确保模型能够准确预测宏观结构响应。在当前研究的范围内，不假设k和l之间存在直接的数学关系。然后，可以建立由k和l组成的晶格相关数据库，如图6所示。通过构建这个离线数据库，可以准确地高效地实现具有不同晶格的多孔聚合物梁在受弯曲载荷作用下的机械变形行为的在线预测。应当注意的是，在本研究中，k和l是根据误差最小化标准进行校准的。具体来说，一旦确定了晶格的几何参数，所提出的均匀化模型就会显示出一组最佳的代表k和l值，这些值最能代表梁的真实机械响应。

4. 结果与讨论
在本节中，基于表面效应的有效弯曲刚度表达式（8）的有效性得到了验证。讨论了TPMS晶格类型、孔隙率和晶格常数对表面强度因子k和表面厚度l的影响，以及它们对超材料梁弯曲刚度的耦合效应。最后，通过构建具有复杂微观结构的多孔聚合物梁，评估了所提出均匀化方法的计算精度和效率。为了方便起见，使用归一化弯曲刚度来表征不同厚度下梁的有效弯曲刚度随表面效应的变化趋势。在本研究中，所采用的模型假设材料是均匀且无缺陷的，没有考虑3D打印过程中的潜在制造误差，如材料缺陷、孔隙率或残余应力。

4.1 表面驱动的尺寸依赖性弯曲刚度
以金刚石晶格为例，研究了内在参数k和l对不同厚度梁的有效弯曲刚度的耦合效应。在计算中，晶格常数设置为...（此处数值缺失），孔隙率为...（此处数值缺失）。通过增加z方向的晶格单元数量，改变了超材料梁的厚度，并获得了归一化弯曲刚度和厚度的变化趋势，从而量化了表面效应对不同厚度梁有效弯曲刚度的影响。结果表明，梁的有效弯曲刚度受到表面效应的显著影响，当厚度为1毫米时影响最大。如图7a所示，...（此处数值缺失）仅为16.7%，这意味着不考虑表面效应时，通过均匀化方法得到的弯曲刚度是考虑表面效应时的六倍。因此，研究表面效应对梁弯曲刚度的影响是极其必要的。然而，也可以观察到，随着厚度的不断增加，某个比率也随之增加，这表明表面效应对梁的弯曲刚度的影响随着厚度的增加而减小。接下来，将使用Diamond晶格作为具体例子来解释为什么归一化弯曲刚度随厚度变化的现象。归一化弯曲刚度随着厚度的增加而持续增加，这表明受表面效应影响的区域与核心区域相比表现出整体刚度的软化现象。这也可以从方程（8）中的两个内在参数k和l来理解。计算结果显示，表面强度因子k小于1，表明表面区域的弹性模量低于核心区域。因此，包含表面效应的梁的有效弯曲刚度表现出软化行为。此外，随着厚度的增加，表面效应对弯曲刚度的影响逐渐减小。这是因为当厚度较小时，表面厚度l与宏观厚度之比相对较大，对有效弯曲刚度的影响更为显著。然而，随着厚度的继续增加，这个比率减小，导致表面效应对有效弯曲刚度的贡献减弱。因此，观察到的刚度软化现象逐渐减弱。

上述分析表明，表面效应的强度取决于表面区域的微观结构，这决定了弯曲过程中两个内在参数k和l的值。k和l之间的耦合相互作用控制了表面效应对有效弯曲刚度的大小。值得注意的是，当梁的宏观厚度接近晶格尺寸时，表面效应变得特别显著。需要指出的是，当前研究是基于材料的线性弹性响应进行的。当考虑材料的非线性行为时，在较高应变或较大载荷下，非线性效应将不可避免地影响表面区域的力学行为。在本研究中，模型没有包括非线性因素。

4.2 微观结构依赖的表面效应

理论和数值模拟都表明，表面效应的大小与晶格的微观结构有关，在TPMS中，这取决于晶格类型、孔隙率和晶格常数。本节通过改变晶格的微观特征来研究表面效应的变化，并阐明了不同类型晶格、孔隙率和晶格常数对表面效应的影响，同时揭示了这些参数影响梁的有效弯曲刚度的根本机制。首先，研究了四种不同晶格类型在弯曲作用下的表面效应，其晶格常数和孔隙率已知。不同晶格类型的归一化弯曲刚度随厚度的变化趋势如图7a–d所示。在研究的四种TPMS晶格中，Diamond晶格表现出最显著的表面效应：当厚度为1毫米时，归一化弯曲刚度为16.7%，而在10毫米厚度时增加到85%。其次是Fischer–Koch S和Primitive晶格，其归一化弯曲刚度分别为51%和65%，而Split P晶格达到70%。值得注意的是，除了Diamond晶格外，其他三种晶格在10毫米厚度时的归一化弯曲刚度值都超过90%，表明在这一尺度下表面效应对有效弯曲刚度的影响非常小。与其他三种晶格相比，Diamond晶格表现出更强的表面效应。从微观角度来看，可以这样解释：与其核心区域相比，表面微观结构更为分散，导致在弯曲过程中缺乏有效的载荷传递路径。因此，表面区域更容易发生局部变形，从而产生更明显的表面刚度软化效应。这种描述旨在为宏观梁刚度的降低提供物理解释。

相比之下，其他三种晶格的表面微观结构连接相对紧密，限制了在弯曲载荷作用下的变形。因此，尽管仍然存在表面刚度软化现象，但其程度远小于Diamond晶格。特别是Split P晶格几乎没有孤立的表面区域，这解释了为什么在四种晶格类型中其刚度软化现象最弱。

其次，以Diamond晶格为例，研究了内在参数在不同孔隙率条件下对有效弯曲刚度的影响机制。分别设置了孔隙率为0.6、0.7、0.8和0.9的四种Diamond晶格，以观察内在参数和表面效应的变化模式。其对归一化弯曲刚度的影响如图8a所示。可以看出，随着孔隙率的降低，表面效应对归一化弯曲刚度的影响也越来越小。当厚度为1毫米时，随着孔隙率从0.9降低到0.6，归一化弯曲刚度从8.16%增加到24.5%，这说明对于同一类型的晶格，孔隙率越大，表面效应对归一化刚度的影响越显著。在分析孔隙率对表面效应的影响时，我们没有直接使用全局刚度的绝对变化来评估其大小，而是采用了相对测量方法。如图8a所示，对于相同厚度的梁，孔隙率越高，表面效应越明显。图8b还显示，内在参数k随着孔隙率的增加而减小，而l则表现出相反的趋势，表明表面的相对刚度软化变得更加显著。表面效应随孔隙率的增加是一个相对变化，并不与整体刚度的减小相矛盾。通过引入这一相对指标，我们可以有效地区分表面效应的贡献和全局刚度的降低。值得注意的是，当孔隙率为0.9且厚度为1毫米时，使用传统均匀化方法得到的弯曲刚度是使用考虑表面效应的均匀化方法得到的12倍。这一发现表明，在弯曲变形下，两种均匀化方法得到的梁响应之间的偏差可达12倍，这无疑与经典力学直觉相悖。

图8. 孔隙率对多孔聚合物梁归一化弯曲刚度的影响。(a) 不同孔隙率梁的归一化弯曲刚度，a代表晶格常数。(b) 不同孔隙率下的表面厚度l和表面强度因子k。模型指的是方程（8）中表达的模型，homo代表表面驱动的均匀化方法。图8b还示出了不同孔隙率下Diamond晶格的表面厚度l和表面强度因子k的变化趋势。可以观察到，随着孔隙率的降低，表面强度因子k逐渐增加，而表面厚度l持续减小。这表明随着孔隙率的降低，受表面效应影响的区域变得越来越小，表面区域的刚度逐渐增加，导致晶格的整体刚度软化效应减弱。这些发现表明，对于同一类型的晶格，降低孔隙率可以有效减弱表面效应对有效弯曲刚度的影响。此外，我们还对有限元模拟中使用的网格密度对内在参数k和l的敏感性进行了系统分析。在数值模拟中，考虑了三种不同的网格尺寸（粗、中等和细），以评估它们对校准参数的影响。如图8c所示，不同网格离散化下k和l的值保持高度一致，说明所提出的校准方法具有良好的鲁棒性和网格独立性。

此外，还研究了晶格常数对表面效应的影响，研究对象仍然是Diamond晶格。将晶格常数分别设置为1毫米、2毫米、3毫米和4毫米，得到的归一化弯曲刚度如图9所示。随着晶格常数的增加，表面厚度l也随之增加，表明晶格的尺寸依赖性行为变得更加明显。然而，尽管晶格常数持续增加，表面强度因子k保持不变，因为k用于表征表面区域的弹性模量与体材料弹性模量之间的比率。晶格常数的变化只改变了宏观几何尺度，而晶格的局部微观结构特征保持不变。因此，k不随晶格常数变化。这表明，在本研究的范围内，k可以被视为晶格微观结构特征的有效描述符，并且具有尺度独立性。

表1展示了不同晶格类型、晶格常数和孔隙率对应的体积模量，以及校准后的内在参数k和l。如前所述，k和l表现出耦合效应；当k较小而l较大时，表明晶格表面区域的刚度软化效应更为显著，从而对有效弯曲刚度的影响更大。这一结论在方程（8）中得到了体现，表1中Diamond晶格的内在参数k和l进一步验证了这一推断。此外，结果还表明k和l对表面行为的影响不是独立的，而是相互耦合的。因此，仅基于其中一个内在参数来评估表面效应对有效弯曲刚度的影响是不够的。必须同时考虑这两个参数，才能准确评估表面刚度软化的程度及其对超材料梁整体力学行为的影响。在本研究中，所提出的模型假设表面区域是均匀且无缺陷的，并且表面强度因子k是在理想化条件下通过数值校准确定的。在这个假设下，k是一个无量纲参数，其值不依赖于材料属性，反映了微观结构约束对表面刚度的贡献。然而，在实际的3D打印过程中，可能会出现诸如表面粗糙度、固化不完全、孔隙或残余应力等缺陷。这些因素可能会在微观尺度上局部改变表面区域的机械性能，从而影响表面硬度。尽管如此，如果这些缺陷在整个结构中相对均匀分布，它们对k的影响是有限的，并不会显著改变模型对宏观梁硬度的预测。表1. 不同晶格的固有参数k和l。4.3. 基于表面驱动的均质化方法的准确性和效率在前两部分中，我们详细分析了表面效应对多孔聚合物梁在弯曲变形中的作用机制，以及固有参数k和l的耦合影响。基于此，我们提出了一种考虑表面效应的均质化方法，可以预测多孔聚合物梁在受弯曲载荷作用下的变形行为。在本节中，如图10b所示，建立了一个包含表面效应的均质化欧拉-伯努利梁模型，以预测梁在弯曲过程中的机械响应。梁的宽度用W表示，厚度用H表示，根据欧拉-伯努利梁理论，长宽比应大于或等于10。通过应用考虑表面效应的均质化方法，可以使用公式(14)计算超材料梁的惯性矩I。图10. 不同厚度多孔聚合物梁的端位移比较。(a) 厚度为的超材料悬臂梁的端位移。(b) 受弯曲载荷的超材料及其均质化悬臂欧拉-伯努利梁模型。(c) 厚度为的超材料悬臂梁的端位移。(d) 不同均质化方法与高保真有限元分析的计算误差比较。(e) 在不同载荷条件下，孔隙率为0.9、梁长为20毫米、高度为1毫米的梁的载荷-位移曲线。根据欧拉-伯努利梁理论，超材料梁在x和z方向上的位移场可以表示为：其中和分别表示x和z方向的位移场，表示中性轴在z方向的位移。x方向的应变可以表示为：因此，超材料梁在弯曲变形后的曲率可以表示为：通过利用方程(8)中给出的有效弯曲刚度与经典弯曲刚度之间的关系，可以推导出有效弯曲刚度的表达式：根据欧拉-伯努利梁理论，平衡方程可以表示为：在悬臂梁的自由端施加固定弯曲力矩M，并使用本文提出的考虑表面效应的均质化方法来预测其弯曲响应，欧拉-伯努利悬臂梁的边界条件可以表示为：由此可以得到梁的位移解：其中梁在处的最大位移由下式给出：同时，以Diamond晶格作为一个典型例子，构建了一个具有复杂微观结构的超材料梁，如图10b所示，它仍然满足欧拉-伯努利梁理论。如果梁的厚度增加或剪切效应变得显著，欧拉-伯努利梁理论的假设可能不再成立，模型的适用性和准确性可能会受到影响。类似地，在梁的自由端施加了弯曲力矩，并采用高保真有限元方法研究了梁在弯曲载荷作用下的端位移。然后将高保真有限元分析的结果与所提出的考虑表面效应的均质化方法的预测结果进行比较，从而证明了该方法的计算准确性。为了提高研究的可靠性，设置了两组具有不同宏观厚度的超材料悬臂梁，分别为和，其中晶格常数分别为和，梁的长度为。对于每组，分别分配了0.6、0.7、0.8和0.9四种不同的孔隙率，以展示宏观厚度和晶格孔隙率对超材料悬臂梁表面效应的影响。在悬臂梁的自由端施加了弯曲力矩。使用经典多尺度均质化方法、考虑表面效应的均质化方法和高保真有限元方法获得了梁的端位移。如图10所示，结果表明，在多孔聚合物梁的弯曲变形中，表面效应非常显著。与高保真有限元结果相比，基于经典力学的多尺度均质化方法不再适用。从图10a,c可以看出，在弯曲力矩M的作用下，无论厚度H是1毫米还是2毫米，经典均质化方法预测的端位移与高保真有限元计算结果有显著偏差，而本文提出的考虑表面效应的均质化方法预测的结果与高保真有限元结果非常吻合。特别是当孔隙率较大时，如图10a所示，当和厚度为时，经典均质化方法预测的位移为0.348毫米，而这里提出的均质化方法预测的位移为4.202毫米，与高保真有限元结果4.233毫米非常吻合，是经典均质化预测的12倍。这清楚地表明，表面效应对超材料梁抗弯曲变形的能力的影响不可忽视。为了更直观地说明表面效应的影响，经典多尺度均质化方法、考虑表面效应的均质化方法和高保真有限元方法计算的弯曲刚度分别表示为、和。以作为参考真值，它们的相对误差通过以下公式计算：结果如图10d所示。从图中可以看出，不考虑表面效应的均质化方法计算的弯曲刚度的相对误差非常大。当时，最小误差达到233%，最大误差高达1108%，远超允许的误差范围。相比之下，考虑表面效应的均质化方法计算的弯曲刚度的相对误差始终低于4%，在允许的误差范围内。这表明该方法的预测是可靠的。应当注意的是，本研究提出的模型是基于欧拉-伯努利梁理论的，该理论假设梁截面在弯曲过程中保持平面且垂直于中性轴，从而忽略了横向剪切变形。当梁的厚度进一步增加时，剪切效应可能变得显著，可能导致模型预测与实际机械响应之间的较大差异。对于厚梁或剪切效应明显的情况，可能需要Timoshenko梁理论或其他考虑剪切变形的建模方法来保持预测准确性。然而，在相对误差低于4%的要求下，如图10d所示，当前模型对于以弯曲为主的情况非常适用，此时欧拉-伯努利梁假设仍然有效。此外，还讨论了晶格常数对超材料悬臂梁表面效应的影响，以进一步验证考虑表面效应的均质化方法的准确性。构建了两组晶格常数分别为a=1毫米和a=2毫米的梁。梁的宏观厚度和长度保持不变，每组设置了四种不同的孔隙率。的结果取自图10a。在梁的一端施加了固定弯曲力矩，并获得了端位移和弯曲刚度的相对误差，如图11b所示。结果表明，对于不同的晶格常数，考虑表面效应的均质化方法的预测与高保真有限元结果吻合得很好，误差都在4%以下，而经典均质化方法的误差明显较大。图10e展示了用于评估等效弯曲刚度的载荷-位移曲线。由于本研究假设材料具有线性弹性行为，因此计算出的位移与施加的载荷呈线性关系。此外，在不同载荷水平下，高保真有限元方法得到的结果与所提出的均质化方法的预测结果非常吻合，表明所提出的模型可以可靠地预测宏观梁的机械响应和表面效应。图11. 不同晶格常数多孔聚合物梁的端位移比较。(a) 晶格常数为的超材料梁的端位移。(b) 不同均质化方法计算的相对误差比较。为了证明本文提出的均质化方法在保持高准确性的同时可以显著提高计算效率，本文对高保真有限元方法和均质化方法的准确性和计算时间进行了比较。超材料梁由Diamond晶格单元组成，分别用L和H表示梁的宏观长度和厚度。晶格常数设定为，孔隙率为。边界条件和载荷与之前描述的相同。计算并比较了梁端的位移，结果见表2。表2. 表面驱动的均质化方法与高保真有限元方法之间的计算准确性和效率比较。所有模拟都在配备了两个AMD EPYC 7502处理器（每个处理器具有32个核心和2.50 GHz的时钟速度）和128 GB内存的高性能工作站上进行，模拟过程中使用了共享内存并行计算。如表2所示，表面驱动的均质化方法不仅确保了高计算准确性，还显著减少了计算时间。效率的提升归因于超材料梁中丰富的微观结构细节，模型识别和处理消耗了大部分数值模拟时间。特别是当超材料梁的尺寸较大时，计算资源的限制可能会阻止高保真有限元方法获得解决方案。然而，在进行均质化后，复杂的微观结构细节不再存在，从而简化并加速了整个解决过程。通过比较高保真有限元方法和考虑表面效应的均质化方法得到的结果，证明了所提出的均质化方法可以有效地预测受弯曲载荷作用的细长多孔聚合物梁的变形行为。在确保准确性的同时，与高保真有限元方法相比，在线计算时间显著减少。5. 结论对于具有复杂微观结构的多孔聚合物梁，当晶格尺寸接近梁的宏观厚度时，表面效应变得显著，在弯曲变形中不能忽略。此外，传统多尺度均质化所基于的经典连续介质假设在这种梁中不再成立，导致预测响应的相对误差较大。此外，多孔聚合物梁复杂的微观结构导致使用高保真数值方法时计算成本过高。因此，开发了一种基于表面驱动的均质化方法来更准确和高效地预测这些梁的弯曲响应。主要发现总结如下：(1) 薄型多孔聚合物梁中的表面效应机制源于顶部和底部的自由表面区域，在这些区域，晶格结构缺乏弯曲载荷作用下核心区域所具有的互锁作用，导致约束不足。因此，这些表面区域对变形的抵抗能力较弱。当超材料梁的宏观厚度接近晶格常数的尺度时，这种效应变得尤为显著，因为自由表面引入的约束缺陷大大降低了梁抵抗弯曲变形的能力。(2)基于所揭示的机制，建立了一个考虑表面驱动惯性矩的均质化模型。表面效应的大小通过两个内在参数来评估：表面厚度l和表面强度因子k。为了校准这些内在参数，开发了一种基于全厚度代表性体积元的表面驱动均质化方法，该方法能够更准确、更有效地预测这些梁的弯曲响应。随后，基于表面驱动均质化方法生成的离线数据集，建立了一个用于尺寸依赖性弯曲响应的在线预测框架。(3)对具有不同厚度、孔隙率和晶格常数的多孔聚合物梁进行了高保真有限元分析的数值模拟。结果表明，本文提出的考虑表面效应的均质化方法对梁在弯曲载荷作用下的变形预测与高保真模拟的结果非常吻合。相比之下，传统的多尺度均质化方法的预测误差最大可达1108%。此外，与高保真有限元方法相比，所提出的均质化方法不仅显著减少了计算时间，还使得复杂超材料梁模型的快速分析成为可能。总之，本研究不仅阐明了多孔聚合物梁在弯曲载荷作用下的变形机制（这些机制与表面驱动的尺寸效应有关），而且还提出了一种准确且高效的工具来预测其弯曲行为。这些发现为多孔聚合物梁在工程应用中的结构设计和性能预测提供了宝贵的指导。