面向纤维增强稳定岩土材料逐层 3D 打印的算法化有限元框架

《Advances in Engineering Software》：An algorithmic finite element framework for layer-by-layer 3D printing of fiber-reinforced stabilized geomaterials

　　

•自动化有限元框架模拟稳定化岩土材料的 3D 打印。
•几何感知算法自动将网格单元分组为打印层。
•刚度和强度随层龄、固化速率和纤维含量演变。
•预测结果与解析基准相符，误差低于百分之二。
•纤维使延展性提高十三倍；更快的固化使可打印墙高增加三倍。

引言
逐层 3D 打印已成为一项重要的数字制造技术 [1], [2], [3]，对建筑工程具有广泛的影响，为减少材料浪费、加快施工进度以及实现此前通过传统方法无法实现的几何复杂的结构设计提供了新机会 [3], [4], [5], [6], [7]。虽然基于挤出成型的 3D 打印已达到显著成熟度 [1], [13], [14]，但将 3D 打印更广泛地应用于基于水泥的材料以外的替代岩土材料和稳定化系统 [15], [16]，包括石灰处理土、生物聚合物改性粘土、微生物诱导碳酸钙沉淀 (MICP) 以及混合稳定剂组合，仍受到关键计算和算法差距的制约。逐层打印的顺序性质，其中每个新沉积的层必须在自身机械性能通过固化或凝胶化过程演变的同时，支撑逐渐增加的上覆荷载，提出了超越单个材料系统的独特建模挑战。这些挑战因添加离散纤维增强而加剧，引入了各向异性材料行为、依赖纤维取向的刚度以及在复合材料尺度上增强的裂缝桥接能力。迄今为止，文献中缺乏能够在通用、模块化算法结构内统一这些物理现象的鲁棒计算框架，特别是那些在标准岩土本构模型中集成自动化几何处理、智能单元管理、时间相关材料演变和异质稳定剂效应的框架 [3], [5], [17]。

推进岩土材料逐层 3D 打印计算建模的根本动力源于建筑部门减少碳排放的迫切需求，同时保持或提高结构性能。仅波特兰水泥生产就占全球人为二氧化碳排放量的 7%–8% [12], [17]，使建筑业成为最大的人为温室气体来源之一。这场碳危机催生了全球对低碳或零碳稳定化技术的研究，包括基于氢氧化钙（熟石灰）的系统、基于生物聚合物和酶的固结以及微生物诱导碳酸钙沉淀。每种替代稳定化方法都具有独特的可持续性优势：稳定土比水泥具有更低的隐含能耗；生物聚合物改性粘土使用可再生、可生物降解的结合剂，加工要求最小；经过 MICP 处理的土壤通过生物代谢实现固结，将隐含能耗降低至接近零，同时提供原位处理的潜力 [18], [19], [20], [21], [22], [23]。然而，每种稳定剂都表现出独特的固化动力学、强度发展轨迹以及与纤维增强的相互作用，创造了多样化的材料系统，这需要统一的、基于物理的计算框架，而不是临时性的、特定材料的模型。加速这些可持续材料在实践中的部署需要计算工具，使设计师和工程师能够快速可靠地预测结构性能、优化工艺参数并评估可行性，而无需广泛的实验活动。

离散纤维增强与稳定化岩土材料的集成已在基于水泥、基于石灰和基于生物聚合物稳定的系统中经过广泛的实验验证，持续显示出拉伸能力、延展性、裂缝控制和峰值后能量吸收的改善 [4], [24], [25], [26]。纤维增强通过众所周知的微力学原理运作：纤维跨越不连续处桥接传播的裂缝，跨纤维 - 基体界面的摩擦应力传递，以及渐进式纤维拔出和断裂，将脆性破坏转化为延性、吸能行为 [4], [24], [25]。在逐层 3D 打印的背景下，纤维增强通过其潜在能力桥接弱的层间界面并减少挤出过程中优先丝状取向引起的机械各向异性，具有额外的意义 [27], [28], [29]。然而，尽管有这些明显的机械效益和已验证的原理，纤维增强与 3D 打印模拟的计算集成仍然有限，特别是当结合时间相关材料演变和特定稳定剂的固化动力学时。

目前 3D 打印模拟的有限元方法主要开发于金属增材制造（粉末床熔融、定向能量沉积）和基于水泥的 3D 打印领域 [1], [13], [14]，但每种方法在应用于更广泛的具有纤维增强的 3D 打印稳定化岩土材料挑战时都表现出显著局限性。

金属增材制造模拟方法：金属增材制造的计算框架，实施于 Abaqus、ANSYS 和 Simufact 等商业软件中 [2], [13], [30]，采用耦合热力学分析，渐进式单元激活与激光扫描或电子束路径同步。这些方法将热分析（跟踪热源移动、相变、热膨胀）与力学分析（应力、塑性应变、变形、残余应力）耦合，单元激活与工具路径几何形状同步。然而，金属增材制造与基于挤出的岩土材料 3D 打印在空间尺度（激光光斑直径 0.1–1 mm 对比挤出喷嘴 5–50 mm）、热梯度强度（10?–10? K/m 对比 10–102 K/m）、时间尺度（每层毫秒对比每层秒到分钟）和固结机制（熔融 - 凝固对比物理沉积和固化）上存在根本差异。这些区别意味着金属增材制造模拟架构与岩土材料 3D 打印根本不匹配，使其直接转移不切实际。此外，金属增材制造框架强调热效应和残余应力，这些在室温下基于挤出的岩土材料 3D 打印中是次要或不存在的。

基于水泥的 3D 打印框架：对于基于水泥的 3D 混凝土打印 (3DCP) [2], [31], [32]，计算建模已取得显著进展，多个研究小组提出了数值策略。Vantyghem 等人 [11] 提出了结合单元生成和摩尔 - 库仑塑性的有限元策略，以预测打印期间的结构能力并识别失效模式。然而，他们的方法采用显式模型变更，为每层或层段创建单独的分析步骤，这是一种繁琐的策略，对于需要数百个单元集的复杂几何形状变得难以管理。研究人员开发了光滑粒子流体动力学 (SPH) 模型，捕捉挤出过程中的新混凝土流动、层几何演变和纤维取向，但 SPH 方法对于全尺寸结构分析计算成本高昂，且其重点仍然是新鲜状态行为而不是硬化状态承载能力 [33], [34]。Cervera 等人 [30] 实施了具有逐层网格粗化的自适应有限元策略以降低计算成本，但这种方法需要复杂的自适应网格划分算法，并且对时间相关材料属性演变或纤维增强效应提供的指导有限。最关键的是，已发布的 3DCP 计算模型通常分别专注于新鲜状态行为（流动、流变学、层间结合）[35], [36], [37] 或硬化状态结构分析，在整个打印和固化过程中对时间相关材料演变的集成有限，且对纤维增强效应的处理除轶事性纳入外微乎其微。

特定材料局限性：现有框架在很大程度上是为基于水泥的系统量身定制的，固化动力学、材料参数和本构假设被硬编码或显式参数化以适应波特兰水泥水化。当应用于替代稳定剂（石灰、生物聚合物、MICP）时，这些框架缺乏算法灵活性来适应不同的固化时间常数、强度发展曲线和特定稳定剂的材料响应。此外，扩展到新材料系统或打印几何形状通常需要大量的人工返工：重新创建几何模型、手动将网格划分为单元集、重写输入文件以及重新校准材料参数。这种缺乏自动化和通用性严重限制了在统一计算平台上进行参数研究以比较不同稳定剂或几何配置的能力。

单元激活策略挑战：实施渐进式单元激活的三种主要策略各自存在显著的计算和实际局限性：
•显式模型变更方法：Abaqus *MODEL CHANGE 关键字允许在分析步骤之间显式激活或停用单元。对于具有 nL 层和每层 nS 段的结构，这需要 nL×nS 个单独的分析步骤，对于现实几何形状可能需要数百个步骤。每个步骤都需要自己的边界条件规范、载荷增量和收敛容差，导致繁琐的输入文件管理、频繁的求解器初始化导致计算开销增加，以及严重受损的有效参数研究能力。具有 20 层和每层 8 段的墙几何形状需要 160 个分析步骤，这种模型设置和管理负担变得令人望而却步。
•通过事件系列进行渐进式激活：Abaqus *ELEMENT PROGRESSIVE ACTIVATION 配合 *EVENT SERIES 提供了更多的自动化，基于规定的激活场或工具路径 - 网格交点算法激活单元。然而，这种方法需要规定的工具路径/激活场与有限元网格之间精确的几何同步，对于几何修改或网格细化灵活性有限。此外，渐进式激活算法与 USDFLD 子程序（时间相关材料演变所需）之间的相互作用没有很好的记录，导致在激活瞬态期间关于单元龄初始化以及场变量更新的不确定性。
•准静态低刚度方法：给非活动单元分配极低的刚度（例如 E≈0.001 MPa）然后在激活时增加刚度，这避免了创建大量步骤，但引入了数值问题。接触刚度定义变得成问题，需要仔细调整罚参数以保持数值稳定性。从非活动到活动刚度的转变引入了人为的刚度变化，可能触发不需要的塑性变形或收敛困难。此外，低刚度的非活动单元为全球系统行为贡献了虚假柔度，可能掩盖局部应力集中或失稳。然而，这些缺点并非策略固有的：它们主要发生在将非活动刚度规定为比例不当的绝对值、激活引起应力跳变或接触罚函数针对非活动刚度进行调整时。本框架采用此策略，并对这些问题设有明确的保护措施（相对刚度缩放、无应力激活、可忽略的非活动自重以及事后验证），详见第 2.10 节。

这些方法中没有一种为几何感知的单元集管理、时间相关材料演变以及集成在单一、可扩展算法框架内的纤维增强效应提供令人满意的自动化。

缺乏通用材料演变实现：沉积后的时间相关材料属性演变，由特定稳定剂固化、水化或凝胶化过程驱动，在现有框架中得到不足的处理。虽然水泥水化动力学已经在经验和理论上被广泛表征 [33], [35], [38]，大多数 3D 打印模拟框架要么完全忽略时间依赖性（假设属性恒定），要么实施硬编码为水泥的特定稳定剂固化模型。生物聚合物固化动力学、石灰火山灰反应和 MICP 沉淀遵循不同的时间轨迹，具有不同的时间常数和最终属性值。结合这些不同的特征需要依赖场变量的本构实现，其中材料参数与单元“龄”，即单元激活后经过的时间，通过从实验测量校准的函数联系起来。然而，大多数已发布的框架没有提供透明的、通用的机制，在不求助于自定义 UMAT 子程序的情况下，在标准摩尔 - 库仑塑性中实现这种依赖龄的材料演变，这大大增加了实现复杂性并降低了在有限元平台之间的可移植性。

纤维增强集成缺陷：离散纤维增强引入了现有框架处理不足的异质性和各向异性。虽然纤维对有效属性的影响通过微力学均匀化在理论上被充分理解，但实际集成到 3D 打印模拟框架中仍然是临时的。大多数现有框架要么完全忽略纤维效应，要么采用经验强度增强因子，而不系统地依赖于纤维体积分数、取向或固化状态。微尺度纤维几何（长度、直径、长细比、界面粘结强度）、宏观尺度取向分布（随机性、挤出诱导的对齐）和有效属性（弹性模量、凝聚力、峰值后延展性）之间的关系本质上是非线性的且依赖于稳定剂 [27]。挤出过程引起沉积方向上的优先纤维取向，形成横向各向同性，其特征是平行于打印方向的刚度和强度高于垂直方向 [27], [29], [39], [40]。现有框架缺乏在标准本构模型内指定依赖纤维含量和取向的材料参数的系统机制，迫使从业者要么忽略纤维效应，要么实施自定义本构子程序 [41]。

除了现有方法的局限性之外，纤维增强稳定化岩土材料逐层 3D 打印固有的基本计算挑战是多方面的且相互关联的。首先，顺序材料沉积过程从根本上违反了标准有限元假设，即问题域是固定的，所有材料从分析开始时即存在 [11], [42]。每个单元必须在特定时间“生成”或“激活”，对应于其层位置和打印顺序，其随后的应力 - 应变历史由后续沉积层施加的演变载荷塑造。对于应用于具有复杂几何形状（墙、环、晶格结构）岩土材料的基于挤出的材料沉积，目前尚不存在普遍适用的、几何感知的激活方法（上述已回顾了各个激活策略的权衡）。

沉积后的材料属性演变，由特定稳定剂固化、水化或凝胶化过程驱动，引入了弹性模量、凝聚力和内摩擦角的显著时间依赖性，这违反了标准本构模型的时间不变参数假设 [19], [20], [21]。水泥水化动力学由考虑温度效应的等效龄指数或 S 形函数描述 [33], [35], [38]。生物聚合物固化动力学遵循由凝胶化、交联和水分重新分布驱动的类似时间依赖轨迹，具有不同的时间常数和最终属性值 [19]。

石灰火山灰反应发生在更长的时间尺度上（天到周），而 MICP 沉淀取决于细菌活性和营养可用性，引入了生物变异性 [18]。逐层构造由于层间界面较弱和每层内的平行纤维结构，本质上会产生各向异性的结构行为 [27], [29], [34], [43]。实验研究一致记录了各种系统的层间/层内强度比 [29], [39], [40], [43] 范围在 0.5 到 0.9 之间。界面弱点源于非矩形挤制截面导致的接触不完整、表面水分和刚度条件、底层固化不完整，以及几何缺陷如空洞和微裂纹 [43], [44]。数值建模需要显式内聚区单元及牵引 - 分离定律，或通过场变量依赖的材料软化在层边界处设置减弱强度区，或使用完全各向异性塑性模型。每种方法都在物理保真度、实现复杂度和计算成本之间涉及权衡。通用框架必须提供灵活性，在需要时纳入显式界面建模，同时支持用于快速参数探索的更简单模型。在本研究中，演示的案例采用了这些选项中最简单的一种，即通过连续位移场实现完美的层间粘结；然而，该框架仍然支持通过边界单元行的场变量依赖强度减弱来实现近似界面减弱，显式内聚区界面单元作为模块化扩展保留（讨论部分）。

离散纤维增强引入了空间异质性和各向异性，虽然理论上适合显式纤维级离散化，但在包含数百万根纤维的结构尺度上计算上是不被允许的 [27]。实用策略依赖于均匀化方法，通过微观力学关系或经验校准将纤维效应纳入有效连续介质特性 [4], [24], [25]。挤制引起的纤维对齐可在层尺度上产生横向各向同性。本框架有意采用均匀化表示，其中纤维效应通过体积分数依赖的有效（各向同性）特性进入，避免了会增加复杂度和降低可移植性的完全自定义各向异性本构实现；显式方向依赖各向异性作为模块化扩展保留（见讨论、局限性和未来研究）。

实际的 3D 打印结构表现出多样的几何形状，平面墙、圆柱环、弯曲壳体、晶格填充，每种都需要适当的离散化和网格划分策略，以便在模拟期间系统地激活层和层段。对于墙体等简单几何形状，基于轴对齐包围盒的笛卡尔坐标划分提供了一种自然的方法。对于圆柱几何形状，圆柱坐标变换允许按垂直和角度位置划分。对于复杂的自由形式几何形状，划分可能需要基于参数化工具路径的符号距离计算、投影到参数化表面或基于邻近度的聚类算法 [30]。手动单元集创建和激活调度是繁琐且容易出错的。通用框架必须包含感知几何的单元集创建算法，自动将任意有限元网格划分为层和层段组，并具有鲁棒的验证检查，确保完全覆盖、无重叠，以及跨网格细化水平的一致性。

稳定化地质材料的本构建模需要塑性模型来捕捉压力依赖的屈服、塑性剪胀以及应变硬化或应变软化行为。莫尔 - 库仑屈服准则仍然是岩土工程实践中最广泛采用的模型，因其简单性、物理透明性和在所有主要有限元代码中的可用性 [45], [46], [47]。对于 3D 打印背景下的纤维增强稳定化地质材料，通过非关联流动法则和场变量依赖的内聚力和摩擦角增强的莫尔 - 库仑塑性，在物理保真度和可处理性之间提供了务实的平衡 [45], [46]。

总之，目前没有可用的框架能同时提供感知几何的自动单元管理、跨多种稳定化系统的广义时间依赖材料演化，以及在统一、模块化、可扩展的算法平台内系统的纤维增强集成。

本文的目的是提出一种通用计算框架，专门为纤维增强稳定化地质材料的逐层 3D 打印而开发，通过集成算法创新克服这些局限性。框架被设计为：（i）跨多种稳定化系统（水泥、石灰、生物聚合物、MICP、混合组合）和纤维类型通用，只需材料特定校准而无需算法重新实现；（ii）结构模块化，具有可替换的几何类型、单元激活策略、材料演化函数和本构模型组件；（iii）几何处理、网格划分和输入文件生成自动化，大幅减少用户工作量和出错倾向；（iv）可扩展至高级功能，如显式内聚区界面单元、各向异性织物描述、多物理场耦合和替代本构模型，仅需对核心进行最小修改；（v）通过全面的算法文档实现实现透明，确保可重复性、验证和其他研究小组的采用。

该框架集成了基于 Python 的几何生成、网格划分和输入文件预处理自动化，以及采用标准莫尔 - 库仑塑性增强的 Abaqus/Standard 有限元分析，通过 USDFLD 实现材料参数演化，链接到纤维含量（通过有效、均匀化特性）和单元年龄。主要贡献包括：1.感知几何的单元集创建算法：适用于笛卡尔（墙）和圆柱（环）坐标系的通用算法，自动将有限元网格划分为层段组（附录 A 中的算法 A.3 和 A.4），具有覆盖一致性和无重叠的验证检查。讨论了自由形式几何形状的扩展路径。2.统一单元激活策略：基于金属增材制造和焊接模拟中成熟的单元激活方法，将鲁棒的单分析激活方法与 USDFLD 子程序和感知几何的单元集集成，避免了显式模型更改（过多的分析步骤）和渐进激活（几何约束）的局限性，同时保持计算效率。3.基于 USDFLD 的材料演化实现：使用标准、文档完善的 USDFLD 场变量机制，实现将弹性模量、内聚力和摩擦角同时与纤维体积分数和单元年龄耦合，通过用户提供的校准函数，适应不同的稳定化固化动力学，而无需自定义 UMAT 开发。4.自动输入文件处理工具：基于 Python 的处理器，解析、验证并增强 Abaqus .inp 文件，添加场依赖材料表、激活调度逻辑、全面输出规范和鲁棒性检查，自动化繁琐且容易出错的模型准备。5.全面的框架文档：通过伪代码和实现示例提供详细的算法描述，确保透明理解、可重复性和采用。6.广义适用性演示：参数研究探索框架在不同纤维含量、固化时间常数和本构参数下的灵活性，确立该方法在不同稳定化 - 纤维 - 几何组合下通用效用的信心。

**章节片段**

**数学公式与本构框架**
本节阐述了纤维增强稳定化地质材料逐层 3D 打印有限元建模的数学基础。该公式涵盖了准静态弹塑性控制方程、场变量依赖本构关系和渐进单元激活策略。始终强调物理清晰度和计算可处理性，以促进在标准有限元代码中的实现。

**算法实现与框架设计**
本节提出的计算框架代表了生成、自动化和管理逐层 3D 打印有限元模型的完整端到端解决方案。该框架集成了多种算法创新与实用软件工程，消除手动、容易出错的模型准备步骤，同时保持对多种几何形状、材料系统和打印参数的灵活性。本节提供了足够的实现细节，以便

**验证与参数研究**
本节通过多种验证方法和全面的参数分析展示了框架的能力。验证策略采用三种互补方法：（i）简单情况的分析基准，（ii）与理想化测试场景的比较，以及（iii）系统的参数敏感性研究。这些共同确立了框架数值实现的信心，同时探索其在不同场景下的预测效用。

**讨论**
本节将工作置于现有文献中进行背景化，检查其理论和实践意义，承认当前局限性，并确定有前景的未来研究方向。

**结论**
本文提出了一种用于纤维增强稳定化地质材料逐层 3D 打印计算建模的综合算法框架。通过集成基于 Python 的自动化、USDFLD 场变量实现和 Abaqus/Standard 有限元分析，该框架消除了 3D 打印结构系统计算设计的主要障碍。主要贡献和发现总结如下：1.该框架适应多种

**CRediT 作者贡献声明**
Nitin Tiwari: 写作 – 初稿，验证，监督，软件，项目管理，方法论，调查，正式分析，数据整理，概念化。Suvechha Dhakal: 写作 – 审查与编辑，软件，调查，正式分析，数据整理。Bijay Nepal: 写作 – 审查与编辑，验证，软件，数据整理。

**利益冲突声明**
作者声明他们没有已知的竞争财务利益或个人关系，这些可能会影响本文报告的工作。

Nitin Tiwari|Suvechha Dhakal|Bijay Nepal