针对土-结构相互作用问题专门设计的二维高阶边界条件

《Soil Dynamics and Earthquake Engineering》：Special development of 2D high-order boundary conditions for soil-structure interaction problems

　　摘要：本研究提出了一种新的二维高阶吸收边界条件，用于基于有限元方法的土-结构相互作用系统直接分析。该条件基于有限元框架中的应力-位移场，具有较高的计算效率和精度。所提出的边界条件包含了源自Lysmer–Kuhlemeyer算子与Higdon算子的终止项，以及基于沿吸收边界定义的辅助变量法向变分的两个补充方程。这一公式能够确保动态模拟中的数值稳定性，使得在截断边界处可无反射地吸收各类入射波和衰减波。此外，在正交吸收边界的交叉处采用了近似角条件，从而保证波的兼容性，避免角节点处的反射，同时无需引入额外的自由度。该边界条件可在有限元框架的单元面层面轻松实现，无需专门设计的边界单元或吸收层，因此能够以大幅缩减的计算域实现无界介质中波传播的精确建模。通过对接坝-基岩系统的频域和时域分析，与采用传统Lysmer–Kuhlemeyer吸收边界的更大规模模型相比，验证了该吸收边界条件的性能。在水平和垂直地震激励下的分析结果均显示出该公式的高精度、稳定性和高效性。

引言：考虑土-结构相互作用能够显著提升地震响应预测的准确性。在地震工程中，结构响应在很大程度上取决于基础土壤的建模方式以及地震荷载的施加方式。土壤建模方法多种多样，从简单的弹簧-阻尼器模型到解析阻抗函数，再到有限元法、边界元法以及缩放边界元法等全数值方法。这些方法旨在满足无限介质中的辐射条件，并考虑波的辐射阻尼效应。本研究重点关注的基于有限元方法的直接分析方法中，地震荷载通常分为三种输入机制：（1）自由场输入，（2）反卷积输入，（3）波注入输入。在第一种机制中，直接将自由场记录作用于土-结构界面，或等效地为结构施加惯性力[1]。第二种机制则在截断边界处施加反卷积后的自由场运动[2]。第三种机制在具备相关数据时，将近断层力或位移记录传递到包围局部计算域的围层中，通过对该层进行地震激励来进行直接分析[3,4]。此类方法包含多种形式，其生成激励层内地震荷载的策略各不相同。其中，域缩减法被认为是最具影响力且应用最广泛的波注入技术[[5], [6], [7]]。在这三种输入机制中均可以对结构的非线性行为进行建模。不过，自由场法和反卷积法依赖于土壤行为线性的假设[[8], [9], [10], [11]]。而在波注入方法中，非线性仅存在于激励层的有限区域内，周围的土壤区域则被假设为线性介质[[5], [6], [7]]。此外，当缺乏准确的近断层记录或地震源距离研究点较远时，波注入法在缩减计算域方面的效果会大打折扣[10]。无论地震源位置如何或非线性建模能力如何，要对土壤或岩石进行准确建模，都需要合适的辐射边界条件来描述波的传播过程[8,9,12]。换言之，缩减计算域的关键在于一种能够以较低成本准确再现无界波传播的稳健边界条件。若没有这样的条件，截断边界就会引入人为反射，从而污染内部解。而依赖传统的粘性边界则需使用较大的周围域来吸收出射波[10,11]。因此，在非线性分析中往往忽略土-结构相互作用现象[13]；即便加以考虑（例如通过域缩减法），所需计算域依然很大[10,11]。为解决这些问题，人们开发出了多种吸收边界技术，其中吸收边界条件和完美匹配层是最为突出且研究最为广泛的两种策略。随着波问题复杂性的增加，有效的边界处理方法也得到了显著发展。这类方法通常可按时间和空间分为全局方法和局部方法。全局方法（如边界元法、缩放边界元法、DtN方法、连续分数边界等）具有高精度，能够让截断边界置于靠近研究区域的位置，从而最大限度地减少系统的自由度[[14], [15], [16], [17], [18], [19], [20], [21], [22], [23]]。然而，由于它们具有时间上的非局部性以及需要复杂的卷积积分，计算成本较高。这些方法不在本文讨论范围之内，故不作详细阐述。相比之下，局部方法将时间和空间特性分开处理，虽然计算效率较高，但往往难以达到所需的精度。早期的局部方法，如经典的Lysmer–Kuhlemeyer吸收边界条件[24]，是为平面波的垂直入射而设计的；它们简单且稳健，但仅能提供一阶近似，主要适用于固体地球物理学领域[25]。后续的发展带来了更为复杂的技术，尤其是高阶吸收边界条件[26,27]以及完美匹配层[28,29]。这些现代方法提高了精度和灵活性，可应用于更广泛的波问题，包括标量场和矢量场[[30], [31], [32], [33]]、各向异性介质[[34], [35], [36]]，以及色散或对流环境[35,37]。本引言概述了这些吸收边界策略，分析了它们的优点与局限性，并简要介绍了旨在结合二者优势的最新混合公式。在波模拟中，如何在不产生非物理反射的情况下截断开放域是一个核心挑战，尤其是在长时间或高频问题中，微小的边界误差会不断累积。为此人们开发了多种边界处理方法，旨在模拟无界外部环境，使出射波能够无缝离开。理想的边界条件应具备以下特点：无反射、长时间稳定、精度高、计算效率高且易于实现。目前已有不同的方法被提出，每种方法都有其各自的优缺点：经典吸收边界条件与高阶吸收边界条件：1969年提出的Lysmer–Kuhlemeyer吸收边界条件[24,25]是最早也是最实用的弹性波模拟截断方法之一。该方法在边界上设置与边界垂直的阻尼器，对于垂直入射的P波和S波较为有效。尽管因其简单性和稳定性仍被广泛应用，但其精度仅为一阶，对于斜入射或复杂波形而言可能会产生显著反射[26,27]。为克服传统低阶吸收边界条件的局限性（如Engquist–Majda[38]和Bayliss–Turkel[39,40]提出的方法），人们引入了高阶吸收边界条件。从理论上讲，通过提高阶数，这些边界条件可以实现任意精度，进而降低出射波的反射系数[[41], [42], [43]]。高阶吸收边界条件通常是递归构建的：从低反射基算子开始（如Higdon型算子[27]），将这些基算子相乘得到更高阶的算子，引入辅助变量以消除高阶导数，构建与波动方程相对应的方程式，并确保边界上仅出现切向导数和时间导数。这些步骤使得高阶吸收边界条件能与标准的有限元和有限差分格式相兼容[44,45]。其中的典型代表包括Collino[46,47]、Givoli和Neta[41,42]以及Hagstrom和Warburton[43,48]的研究工作，后者还为标量波动方程提出了完全辐射边界条件[44,49]。其主要挑战在于消除辅助变量的法向导数的代数复杂性，尤其是在弹性问题中[45,50]，还有角点处理问题——在三维情况下，由于存在边缘和角点，这一问题会更加严重。这也解释了为何尽管在二维中已成功实现了高阶吸收边界条件，但在三维弹性动力学中仍不存在稳定且完全局部的此类条件[43,46,51]。完美匹配层：1990年代中期，Bérenger[28,29]提出了完美匹配层方法，其概念与吸收边界条件有所不同。该方法并非在边界上直接指定条件，而是在计算域周围添加一层虚拟的吸收层，在该层中对波动方程进行修改，从而使出射波呈指数衰减。这层吸收层被设计为与内部域“完美匹配”，这样任何进入吸收层的波都不会发生反射。理论分析表明，完美匹配层在连续形式下非常有效，但在离散化后其性能可能会下降[8,9,12]。这种性能下降源于对数值参数的敏感性、需要特定的阻尼或拉伸剖面，以及在离散层面确保匹配条件的困难[31,52]。此外，与高阶吸收边界条件不同，除非吸收层足够宽且离散化程度足够高，否则完美匹配层缺乏明确的收敛性概念——而这往往会导致较高的计算成本[44,50]。尽管如此，由于其易于实现，尤其是在处理角点及复杂几何形状时，完美匹配层仍被广泛使用。它也更适合处理各向异性和异质介质，不过也有报道称在各向异性弹性介质中会出现不稳定现象[34,36,53,54]。高阶吸收边界条件和完美匹配层均为解决波截断问题提供了可靠的方案，但两者各有特点，适用于不同的应用场景（见表1）。最近的研究揭示了这两种方法之间有趣的理论联系。例如，某种带有虚数拉伸特性的非标准完美匹配层设计，可以被证明等同于高阶完全辐射边界条件吸收边界条件[43,48]。这些发现表明，尽管两者在概念上存在差异，但在离散化后，它们或许可以被视为同一数学结构的不同表现形式。混合方法：鉴于局部吸收边界条件和完美匹配层各自的优缺点，研究人员提出了结合二者优势的混合技术。其中一种方法是双吸收边界法，它在由薄计算层分隔的两条平行人工边界上分别应用局部高阶吸收边界条件[55,56]。在双吸收边界法中，不仅在两条边界上，还在中间层中引入辅助变量，从而简化了公式，无需消除法向导数。该方法既保留了高阶吸收边界条件的收敛性保证，又具备了完美匹配层的几何灵活性。此外，双吸收边界法无需特殊的角点处理，在有限元格式中也很容易实现[57]。将其应用于弹性动力学时，已取得了令人满意的结果。一项重要创新在于，在递归公式中加入了经典的Lysmer–Kuhlemeyer条件作为终止步骤，从而保证了长时间稳定性——这是以往的高阶弹性吸收边界条件所不具备的特性[55,58]。通过结合标量Higdon算子以提高精度，以及矢量LK算子以确保稳定性，形成了一个既稳健又高效的方案。受早期关于瑞利波吸收研究的启发[59]，这种方法已被扩展用于处理均匀介质和分层弹性介质[58]。除了本质上为局部型的混合吸收边界条件（如Higdon型和完美匹配层）之外，还出现了其他将全局阻抗型方法转化为局部型方法的混合技术。这类方法主要分为两大类。第一类是从全局方法如缩放边界元法和DtN法出发，通过连续分数展开[[60], [61], [62], [63], [64], [65]]、有理逼近[[66], [67], [68]]、算子分解[[69], [70], [71]]或双渐近形式[[61,64,72]]等数值算法，衍生出高阶局部吸收边界条件。第二类则是缩放边界完美匹配层，它将缩放边界元法的几何结构与完美匹配层的吸收功能相结合，形成一种局部的三阶吸收层[[73], [74], [75], [76]]。这两类混合方法都能产生精确且局部的边界，但通常需要频域预处理、缩放中心或大量的辅助变量。尽管从理论上看具有有效性，但迄今为止，针对弹性动力学波传播的高阶吸收边界条件大多仅在简化的基准问题上进行了评估，例如一维问题或波导问题。在真实的土-结构相互作用问题中，尤其是在瞬态载荷作用和复杂几何形状条件下，这类边界条件的适用性、数值稳定性以及计算性能仍需进一步研究。而当传统公式被应用于瞬态问题时，其存在的已知稳定性问题更是加剧了这一局限。在实际应用中，当面对具有多个截断边界的域时，高阶边界条件可能会遇到困难，而土-结构相互作用问题中就经常出现这种情况[77]。特别是在由人工边界交汇形成的角节点处，波场的不一致性可能导致虚假反射或数值稳定性下降[78]。因此，对于通过基于有限元法的直接方法分析的SSI问题，所需的边界条件应具备以下特点：收敛性明确，在时域和频域分析中具有适当的稳定性，参数敏感性低，可在单元面层面应用（不依赖于单元类型或网格重生成），能简化角点处理，并且在保持高精度的同时大幅缩减计算域。相较于PML边界条件，此类边界条件更具优势，因为PML缺乏严格的收敛性，对离散化敏感，且无法用于处理来自外部源的入射波。本研究解决了这些问题，从而提出了所用的局部吸波边界条件。在本研究中，一种新的二维高阶吸波边界条件在有限元应力-位移框架下被开发出来，用以精确截断SSI分析中的无界域。该边界条件引入了由Lysmer–Kuhlemeyer算子和Higdon算子衍生的终止算子，同时还加入了基于辅助边界变量法向变分的额外控制方程。这种构造确保了时域和频域模拟中的数值稳定性，能够有效吸收出射的传播波和衰减波而不产生反射。此外，在正交吸波边界的交点处还引入了近似角点条件，以保障波的兼容性并消除角点反射，而无需增加额外的自由度。所提出的边界条件可直接在单元面层面实现，无需使用特殊的边界单元或吸波层，从而大幅缩减计算域。总之，与全局到局部的混合方法（SBFEM-CF、DtN-CF、SBPML）不同，所提出的方法无需频域预处理、坐标变换或特征值求解；其辅助变量完全为局部变量。据我们所知，这是首个应用于SSI问题的二维高阶局部ABC方法，已在频域和时域得到验证，同时具备以下优点：(i) 完全的局部性，无需频域预处理；(ii) 显式的无自由度角点条件；(iii) 直接在单元面层面实现，无需特殊边界单元或吸波层，因此计算效率高，更易于在实际SSI分析中应用，尤其是参数研究方面。需要强调的是，与文献中的常规做法一致，吸波边界条件是为线性弹性动力学的远场设计的，而潜在的非线性行为则通常局限于结构本身或近场土壤或岩石区域，具体取决于所采用的地震输入机制。通过水平及垂直地震激励下真实大坝-基础系统的频域和时域分析，证明了所提出方法的有效性。与采用传统Lysmer–Kuhlemeyer边界的更大规模模型相比，结果进一步验证了该方法的准确性、稳定性和计算效率。

问题描述
如图1所示，典型的二维土-结构相互作用（SSI）系统将半无限弹性土体视为由人工截断边界围成的有界子域。为保持通用性，假设在土体域的侧面有两个垂直截断边界（即ΓI和ΓII），其法向量分别沿x轴的正负方向指向外部。此外，在域的底部还设置了一个水平截断边界。

数值实验
在对带有假设吸波截断边界的无界域进行数值分析时，要获得足够的精度和效率，关键在于正确定义问题的自然条件以及所采用的截断边界条件类型。这些自然条件包括几何形状、材料属性与行为、加载类型，以及最终用于问题控制方程的离散化方案。

结论
通过将Lysmer–Kuhlemeyer终止项与Higdon算子相结合，开发出了一种用于SSI分析的高阶吸波边界条件。用零算子替换LK项可消除主要的稳定性问题。通过代数运算推导出了应力-位移公式，另外还通过变量变换引入了两个方程，用于处理辅助变量的法向导数，从而完成了动态分析的公式体系。这些边界条件……

作者声明
作者声明没有竞争性利益。

CRediT作者贡献说明
H. Abdorazaghi：形式分析、研究、方法论、项目管理、软件、验证、写作——初稿撰写、写作——审阅与编辑。V. Lotfi：方法论、项目管理、监督、验证、可视化。

利益冲突声明
作者声明没有已知的可能影响本文所述工作的财务利益或个人关系。

H. Abdorazaghi | V. Lotfi
伊朗德黑兰阿米尔卡比尔理工大学土木与环境工程系