Isogeometric suitable coupling methods for partitioned multiphysics simulation with application to fluid–structure interaction

TL;DR

This paper introduces new coupling methods for isogeometric analysis in multiphysics simulations, improving accuracy and efficiency compared to traditional approaches.

Contribution

The paper introduces two spline-based coupling strategies tailored for isogeometric analysis, reducing communication overhead and preserving geometric accuracy.

Findings

Spline-based coupling methods reduce communication overhead while maintaining high-order continuity and exact geometry.

Experimental results align with theoretical predictions, confirming improved efficiency and accuracy.

The methods enable accurate interface communication between isogeometric and conventional solvers.

Abstract

This paper presents spline-based coupling methods for partitioned multiphysics simulations, specifically designed for isogeometric analysis (IGA) based solvers. Traditional vertex-based coupling approaches face significant challenges when applied to IGA solvers, including geometric accuracy issues, interpolation errors, and substantial communication overhead. The methodology draws on the IGA mathematical framework to deliver coupling solutions that preserve the high-order continuity and exact geometric representation of splines. We develop two complementary strategies: (1) a spline-vertex coupling method that enables efficient interaction between IGA and conventional solvers, and (2) a fully isogeometric coupling approach that maximizes accuracy for IGA-to-IGA communication. Both theoretical analysis and extensive numerical experiments demonstrate that our spline-based methods significantly reduce communication overhead compared to traditional approaches while simultaneously enhancing geometric accuracy through exact boundary representation and maintaining higher-order solution continuity across the coupled interfaces. We quantitatively confirm the communication efficiency benefits through systematic measurements of both transfer times and data volumes across various mesh refinement levels, with experimental results closely aligning with our theoretical predictions. Our benchmark studies further demonstrate the geometric fidelity advantages through exact boundary representation, while also highlighting how the inherent mathematical structure of splines naturally preserves solution derivatives across interfaces without requiring additional computation or specialized transfer algorithms. This work not only provides efficient coupling strategies tailored to IGA-based solvers but also establishes a practical bridge between IGA and traditional discretization methods in partitioned multiphysics simulations. By offering viable options for coupling conventional solvers with IGA-based components, our approach enables broader adoption of IGA in established simulation workflows while ensuring accurate and high-performance interface communications.

Figures

Partitioned heat conduction problem setup

Temperature field of a partitioned spline-mesh domain with varying $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _N$$\end{document}$ refinement levels $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{N}$$\end{document}$

Validation of derivative transfer via spline representation, relative $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}$ -error of the heat flux remains at machine precision level

Setup and parameters of the perpendicular flap test case. The original image is taken from the preCICE paper [[8](#CR8)]

Coupled fluid–structure interaction at selected time instances: fluid velocity magnitude (background) and structural pressure distribution (foreground) at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0.3\,\text {s}$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgr

Displacement of the perpendicular flap over time under fluid loading

Multi-field simulation of a system with 20 perpendicular flaps under fluid loading

Visual simulation results of IGA-IGA coupling: fluid velocity magnitude and isogeometric mesh at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.3\,\text {s}$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0

Perpendicular flap benchmark in pure IGA setting, information exchanged through the spline based communication method in a partitioned way coupled through preCICE

Turek Hron FSI benchmark setup

Illustration of vertex-based solvers coupled with IGA-based solvers, with information exchanged at quadrature points/vertices. Shaded region (in grey) depicts the FSI interface. Arrows indicate the direction of data exchange across the interface

Simulation results for the Turek-Hron FSI3 benchmark using G+Smo coupled with OpenFOAM. (a) Vertical displacement of the beam tip over time. (b) Beam displacement and fluid pressure distribution. (c) Beam displacement and fluid velocity magnitude for mesh $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=4$$\end{document}$ with degree $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepac

Spline-based coupling between IGA-based structural solver and vertex-based fluid solver, with interface data exchanged via continuous spline representation

IGA-IGA communication with spline representations on both fluid and structural domains, maximizing the benefits of isogeometric analysis

Constant load test setup for the vertical beam

Left: Comparison of tip displacement between reference solution and different coupling methods. Right: Relative $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}$ error of different coupling methods compared to the reference solution

Nonlinear load test setup for the vertical beam. The nonlinear load $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}(\xi )$$\end{document}$ varies linearly with parameter $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\e

Tip displacement versus time under nonlinear loading

Comparison of the communication overhead between vertex-based (FEM) and spline-based (IGA) coupling. The horizontal axis represents the number of parameter subdivisions per direction, serving as an equivalent indicator of the unique knots in IGA while reflecting mesh resolution in FEM. The vertical axis shows the total communication overhead in kilobytes

Experimental (Exp.) and theoretical (An.) comparison of communication overhead with respect to the number of parameter subdivisions *r* per direction, with for $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\in \{1,2,...,7\}$$\end{document}$

Full text

Introduction

Fluid–-structure interaction (FSI) remains a fundamental challenge in multiphysics simulations, governing a wide range of engineering applications, such as the aerodynamic response of aircraft wings, the pulsatile behavior of cardiovascular systems, and the structural resilience of bridges and offshore structures under hydrodynamic loading. These problems exhibit inherently bidirectional coupling, where the fluid exerts pressure on the structure, causing deformations that, in turn, modify the fluid domain, creating a nonlinear and dynamic feedback loop. Ensuring numerical stability, accuracy, and efficiency in such simulations is especially challenging when heterogeneous discretizations are utilized for the fluid and structural solvers, requiring robust interface coupling approaches [1–3].

Partitioned methods for multiphysics problems

FSI problems are typically tackled via one of the two following solution strategies: monolithic and partitioned formulations. The monolithic approach integrates fluid and structural solvers into a single, coupled system, yielding superior numerical stability and convergence behavior [4] but at the cost of flexibility, often necessitating tightly integrated code bases and customized solvers. In contrast, partitioned methods decouple the problem into distinct subdomains, allowing fluid and structural equations to be solved independently and with different solvers [5]. These subproblems communicate through the iterative exchange of boundary data across the interface. The modularity provided by partitioned schemes promotes solver reuse, algorithmic specialization, and compatibility with legacy software [6, 7]. Moreover, partitioned formulations are often favored over monolithic ones because they offer the flexibility to couple pre-existing solvers, facilitating the integration of specialized codes developed for each physical domain [8]. However, this same modularity introduces significant challenges: ensuring numerical stability, enforcing kinematic and dynamic consistency, and preserving accuracy across potentially non-conforming meshes are nontrivial tasks [9–11]. At the same time, a number of mature approaches successfully address these issues on non-matching interfaces, including (dual) mortar methods [12, 13], common-refinement based conservative/consistent transfer [7, 10], Nitsche-type formulations, and black-box mesh morphing with RBFs [14, 15]. We position our contribution relative to these by leveraging spline spaces at the interface to achieve smoothness, accuracy, and compact data exchange.

Partitioned methods involve two main coupling aspects: a temporal and spatial one. Temporal coupling controls how and when data are exchanged between the fluid and structural solvers. This can range from simple staggered approaches to more advanced methods, such as fixed-point or quasi-Newton iterations that improve stability and accuracy [16]. Spatial coupling, on the other hand, deals with how data is transferred across the interface, especially when the fluid and structural meshes do not match. This usually involves interpolation or projection between different meshes [7, 9, 14]. It is difficult to achieve both high accuracy and good computational performance at the same time. The conventional approach of increasing mesh resolution to improve accuracy directly conflicts with the need to minimize communication overhead, leading to diminishing returns as the size of the problem increases. This challenge is further compounded by the need to balance efficient data exchange rates with parallel execution capabilities [17], particularly in demanding applications such as marine and offshore engineering where large-scale flexible FSI simulations are increasingly critical [18].

Spline-based coupling: a paradigm shift

Isogeometric Analysis (IGA), pioneered by Hughes et al. [19], presents a powerful alternative to traditional finite element approaches by unifying geometry representation and solution approximation. In IGA, the same mathematical basis functions (typically Non-Uniform Rational B-Splines (NURBS) or B-splines) are used to define both the computational domain and the solution field [20]. This seamless integration enables exact geometric representation, even on coarse meshes, and supports higher-order continuity across elements: features especially advantageous for problems dominated by interface behavior.

Partitioned FSI methods have gained widespread adoption over the past decades due to their flexibility in leveraging specialized solvers for different physical domains. However, a persistent challenge in this context lies in the accurate and efficient transfer of data across the fluid–structure interface, particularly when nonmatching spatial discretizations are involved [21]. Spline-based coupling techniques, especially those employing NURBS or B-splines, have shown significant potential to address these issues. For example, Hosters et al. [1] first developed a NURBS-enhanced FSI framework using spline descriptions on both fluid and structure sides, enabling accurate and direct data exchange. Kamensky et al. [22] introduced a stabilized immersogeometric formulation, showing improved stability properties when using collocated constraints across thin, deforming interfaces.

Further developments in spline technology have expanded the scope of isogeometric coupling. Floating IGA (FLIGA), introduced by Hille et al. [23], enables analysis under extreme deformation by adapting basis functions along deformation-dependent directions. Similarly, Rosa et al. [24] proposed a blended IGA-FEM approach for modeling fracture propagation without remeshing, further illustrating the adaptability of spline-based representations to evolving geometries. Guarino et al. [25] presented a novel penalty-based coupling scheme for trimmed non-conformal shell patches, relevant for industrial applications where surface continuity and physical correctness must be preserved.

As the complexity of FSI problems increases, particularly in scenarios involving large structural deformations or compressible flows, robust coupling strategies become even more critical. While spline-based methods offer smooth and accurate interface representations, recent advances in finite element-based coupling have tackled these challenges from a different perspective. Rajanna et al. [26] proposed a compressible FSI framework using nonmatching finite element meshes, suitable for high-speed flow applications. Klöppel et al. [13] developed a dual mortar formulation to improve force and displacement consistency across nonconforming FEM interfaces. These contributions are important within the FEM paradigm, though they differ fundamentally in how geometry and continuity are handled compared to spline-based methods.

Despite advances in interface accuracy and stability, the issue of computational efficiency remains unresolved, particularly in large-scale, strongly coupled simulations. Data interpolation between nonmatching meshes can become a performance bottleneck. Radial Basis Function (RBF) methods have gained traction as black-box interpolants for multiphysics applications, as seen in the work of Lindner et al. [27] and its parallelized extension by Schneider et al. [15]. However, while RBFs offer flexibility, they may lack the conservation and variational structure desired in high-fidelity IGA settings. De Boer et al. [28] highlighted that consistent interpolation often outperforms conservative approaches in terms of accuracy and stability, particularly for quasi-1D FSI problems. In contrast, spline-based coupling inherently supports both consistency and high-order smoothness through function-space projection, offering a more principled route for interface treatment.

In contrast to existing point-based and spline-based coupling methods that mainly interpolate interface data at discrete points, this paper proposes a fundamentally new approach for partitioned FSI simulations by directly expressing interface fields as smooth, continuous spline functions. The key innovation lies in directly utilizing the inherent advantages of IGA to represent the interface as a continuous spline geometry, rather than as isolated interpolation points. In the proposed framework, both the structural and fluid participants maintain their own spline-defined interface fields, allowing for a consistent and symmetric field-to-field coupling. This shift in paradigm offers several compelling advantages: it preserves geometric fidelity for complex curvilinear boundaries, reduces communication overhead via compact control point exchange, improves numerical stability due to high-order continuity, and inherently handles nonmatching meshes through direct value evaluation, eliminating the need for ad hoc interpolation schemes.

The proposed methodology treats the interface as a mathematically structured entity rather than a numerical artifact. This approach maintains physical consistency while achieving high accuracy with reduced computational cost, effectively addressing the longstanding accuracy-efficiency trade-off in spatial coupling for partitioned FSI systems.

Research contributions and paper structure

This work introduces key contributions to multiphysics simulation. While the methods are broadly applicable to general multiphysics problems, we focus without loss of generality on fluid–structure interaction as a representative application throughout this study.

• A theoretical framework for interface coupling in partitioned FSI based on IGA, specifically addressing challenges with non-matching meshes while maintaining accuracy and computational efficiency.
• Two coupling strategies: (i) a hybrid approach that combines IGA-based structural solvers with traditional fluid solvers, and (ii) a fully isogeometric method for high-accuracy IGA-to-IGA coupling.
• An analysis of communication overhead, showing that our spline-based method significantly reduces communication cost, especially with finer meshes.
• Verification and validation using standard FSI benchmarks, demonstrating a strong balance between accuracy and efficiency.
• Implementation as an open-source software package integrated into the preCICE multiphysics framework, featuring a modular plug-and-play design compatible. The remainder of the paper is organized as follows. Section 2 introduces the governing equations used in this study. Section 3 provides an overview of IGA, which serves as the spatial discretization method for the structural domain in this work. Section 4 explains the numerical methods for both fluid and structural solvers and presents our spline-based coupling approach. Section 5 provides verification and validation results to demonstrate the accuracy of our method. The concluding Sect. 8 ends this study, summing up the main results and opening up for further research.

The methodology presented in this work has been implemented in gsPreCICE, developed as part of this research. As the first dedicated IGA solver adapter for the preCICE (Precise Code Interaction Coupling Environment) [8] framework, gsPreCICE is open-source and officially integrated as a module within the G+Smo (Geometry + Simulation Modules) library [29–31], which is a C++ library dedicated to IGA.

Governing equations

Partitioned algorithms tackle the distinct challenges posed by fluid dynamics and structural mechanics by separately solving each domain and exchanging data at a shared interface. While such approaches are broadly applicable to a wide range of multiphysics problems, the scope here is limited to FSI, involving the coupling of incompressible fluid flow with nonlinear elastic structural deformation in multidimensional spaces.

In this section, we briefly review the governing equations of incompressible fluid flow and nonlinear elasticity formulated in a total Lagrangian framework used throughout the paper. Clearly, restating these equations here serves two purposes: (a) it establishes a consistent notation that is essential for accurately describing our numerical methods and coupling strategies, and (b) it provides context for subsequent discussions on interface conditions and numerical discretizations.

Incompressible flow

The fluid dynamics within the domain are governed by the Navier–Stokes equations for incompressible flow, outlined as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\rho ^f \left( \frac{\partial \textbf{u}^f}{\partial t} + (\textbf{u}^f \cdot \nabla ) \textbf{u}^f \right) = - \nabla p^f + \nabla \cdot (\mu ^f \nabla \textbf{u}^f)\nonumber \\&\quad + \textbf{f}^f \quad \text {in } \Omega ^f\times [0,T], \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\nabla \cdot \textbf{u}^f = 0 \quad \text {in } \Omega ^f\times [0,T], \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\textbf{u}^f = \textbf{g}^f \quad \text {on } \Gamma _{D}^f \times [0, T], \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&(\mu ^f \nabla \textbf{u}^f - p^f \textbf{I}) \cdot \textbf{n}^f = \textbf{h}^f \quad \text {on } \Gamma _N^f\times [0,T]. \end{aligned}$$\end{document}

Within these equations, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho ^f$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^f$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{f}^f$$\end{document} represent the fluid’s density, dynamic viscosity, and external body force vector, respectively. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{D}^f$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{N}^f$$\end{document} denote the Dirichlet and Neumann boundaries of the fluid domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega ^f$$\end{document} . This set of equations, (1a) through (1d), is utilized to solve for the fluid’s velocity field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}^f(\textbf{x},t)$$\end{document} and pressure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^f(\textbf{x},t)$$\end{document} at any given time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \in [0,T]$$\end{document} .

Structure deformation

The solid displacement field, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^s(\textbf{x},t)$$\end{document} , which responds to external forces, is governed by Newton’s second law of motion, as expressed below:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d^2 \textbf{d}^s}{d t^2}-\frac{1}{\rho _0^s} \boldsymbol{\nabla }_0 \cdot \left( \textbf{S}^s \textbf{F}^T\right) = \textbf{b}^s \quad \text {in } \Omega _0^s \times [0,T] \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0^s$$\end{document} represents the solid’s density and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{b}^s$$\end{document} denotes the external body force acting on it. This equation is derived from a total Lagrangian framework and is hence formulated with respect to the unstrained reference configuration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _0^s$$\end{document} . All quantities and operators in this formulation are referenced with a subscript ’0’ to indicate this baseline state. Unlike the Cauchy stress tensor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}^s$$\end{document} , the internal stress here is expressed using the second Piola-Kirchhoff stress tensor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}^s = \operatorname {det}(\textbf{F}) \textbf{F}^{-1} \textbf{T}^s \textbf{F}^{-T}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{F}$$\end{document} is the deformation gradient. In this work, we adopt the Saint-Venant-Kirchhoff (SVK) material model with a linear stress–strain relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S} = \mathbb {C}: \textbf{E}$$\end{document} . The overall nonlinear structural response arises from geometric nonlinearity in the Green-Lagrange strain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{E} = \frac{1}{2}(\textbf{F}^T\textbf{F} - \textbf{I})$$\end{document} , which is essential for large deformations. The problem formulation is completed with initial zero displacement conditions and appropriate boundary constraints on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \Omega ^s$$\end{document} . The boundary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \Omega _0^s$$\end{document} is partitioned into two complementary regions: a Dirichlet boundary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _D^s$$\end{document} where displacement is prescribed, and a Neumann boundary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _N^s$$\end{document} where traction is prescribed. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}_0$$\end{document} denotes the unit outward normal vector to the reference configuration, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{g}^s$$\end{document} is the prescribed displacement boundary vector, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{h}^s$$\end{document} is the prescribed traction boundary vector. The boundary conditions are expressed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{d}^s&= \textbf{g}^s \quad \text {on } \Gamma _{D}^s \times [0, T] \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{F}\textbf{S}^s\textbf{n}_0&= \textbf{h}^s \quad \text {on } \Gamma _{N}^s \times [0, T] \end{aligned}$$\end{document}

Coupling conditions

In the partitioned approach, the fluid problem and the structure problem are solved independently, potentially by different solvers. For consistent coupling across the shared interface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma ^{fs} = \Gamma ^f\cap \Gamma ^s$$\end{document} , the interface solutions of both subproblems must fulfill kinematic continuity requirements:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{d}^f = \textbf{d}^s \quad \text {on } \Gamma ^{fs}\times [0,T] \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{u}^f = \textbf{u}^s \quad \text {on } \Gamma ^{fs}\times [0,T] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^f$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}^f$$\end{document} represent the displacement and velocity of the moving fluid boundary, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^s$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}^s$$\end{document} are those of the structure.

The implied accelerations are equivalent too. And the dynamic continuity for the equilibrium of stress:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{T}^f\cdot \textbf{n}^f = -\textbf{T}^s\cdot \textbf{n}^s \quad \text {on } \Gamma ^{fs}\times [0,T] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}^f$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}^s$$\end{document} are the interface unit normal vectors pointing outwards from the corresponding domains. Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}^f = -\textbf{n}^s$$\end{document} on the interface. While the kinematic constraints guarantee displacement and velocity continuity at the interface, the dynamic equilibrium condition maintains stress balance according to Newton’s third law (action-reaction principle). Together, these interface conditions preserve the physical consistency required for accurate fluid–structure interaction modeling [3].

In our partitioned framework, the fluid solver computes tractions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{t}^f = \textbf{T}^f \cdot \textbf{n}^f$$\end{document} on the current (deformed) interface. Since the structural formulation employs a total Lagrangian approach, these tractions must be pulled back to the reference configuration. Using Nanson’s formula ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}\,da = J\textbf{F}^{-T}\textbf{n}_0\,dA$$\end{document} ), the equivalent traction on the reference configuration becomes:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf{t}^s_0 = J^s (\textbf{F}^s)^{-1} \cdot \textbf{t}^f \quad \text {on } \Gamma ^{fs}_0} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J^s = \det (\textbf{F}^s)$$\end{document} is the Jacobian determinant of the structural deformation gradient. This transformation ensures that the resultant force is preserved between configurations, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _{\Gamma ^{fs}} \textbf{t}^f \, da = \int _{\Gamma ^{fs}_0} \textbf{t}^s_0 \, dA$$\end{document} , thereby maintaining physical consistency for large deformations.

Isogeometric analysis overview

In this work, we adopt a partitioned coupling approach, treating fluid and structure as separate subproblems that are solved independently. For the structural subproblem, we employ IGA, which offers both geometric exactness and higher continuity in the solution space.

IGA, introduced by Hughes et al. [19], enhances geometric accuracy and promotes closer integration between CAD models and numerical analysis by using the same basis functions for geometry representation and solution approximation. Typically, CAD geometries are represented by NURBS or b-splines, defined by a set of control points and associated rational basis functions. A NURBS surface can be mathematically expressed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \overline{\textbf{x}}(\boldsymbol{\xi })=\sum _{i=1}^{n_c} R_i^p(\boldsymbol{\xi }) \textbf{x}_i,\quad \text {with }\boldsymbol{\xi }=(\xi ^1,\xi ^2), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_c$$\end{document} is the number of control points, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{x}_i \in \mathbb {R}^{d}$$\end{document} represents the control points, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_i^p(\boldsymbol{\xi })$$\end{document} are rational basis functions of polynomial degree p defined on the parameter domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\Xi }\subset \mathbb {R}^d$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\in \{2,3\}$$\end{document} is the spatial dimension.

IGA provides an ideal framework for our coupling methods because it offers smooth and continuous representations of solution fields across the interface, which are critical for stable and accurate multiphysics coupling. In addition, it enables precise geometric descriptions, which help avoid approximation errors in curved or evolving interfaces.

Numerical methods

In partitioned FSI simulations, the fluid and structure subproblems are solved using separate solvers, each tailored to its respective physical field. This modular strategy offers considerable flexibility: established solvers can be reused or further developed independently, without requiring changes to the overall framework [8].

At the interface between fluid and structure, the two solvers often use different meshes and sometimes even different types of numerical methods. Because of this, we need a reliable way to transfer data (like forces and displacements) between them in a way that is accurate and physically consistent. One of the main contributions of this paper is a spline-based coupling method, which is especially effective when combined with IGA. This method helps maintain smooth and consistent data transfer across nonmatching meshes.

Building on this motivation, the remainder of this section focuses on spatial coupling strategies that bridge the fluid and structural solvers, with particular emphasis on the proposed spline-based methods that support the consistency and efficiency of the partitioned framework.

Evolution of spatial coupling schemes: from vertex-based to spline-based methods

Spatial coupling schemes are critical components in partitioned FSI simulations, determining the manner and accuracy of data exchange at the interface. We consider three representative approaches: vertex-vertex, spline-vertex, and spline-spline coupling, that range from conventional point-based methods to the proposed spline-based formulations.

Vertex-based coupling method: traditional approach and its limitations

The vertex-based coupling method facilitates communication between an IGA-based structural solver and traditional solvers based on discretization techniques such as finite volume (FV) or finite element methods (FEM). In this approach, the data exchanged between solvers is defined at discrete points along the interface, typically mesh vertices or quadrature points (Fig. 1).Fig. 1. Illustration of vertex-based solvers coupled with IGA-based solvers, with information exchanged at quadrature points/vertices. Shaded region (in grey) depicts the FSI interface. Arrows indicate the direction of data exchange across the interface

The vertex-based coupling procedure is performed in three main steps:

• Data mapping for the discrete forces: Vertex-based coupling requires discrete stress components at each fluid node along the FSI interface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma ^{fs}$$\end{document} . To transfer these stresses between nonmatching meshes, a range of mesh-to-mesh interpolation techniques are available, including nearest-neighbor, linear projection, and mortar methods. In this work, we use RBF interpolation, implemented via the preCICE library [15, 27], due to its flexibility, smoothness, and robustness for scattered, non-conforming interface meshes. Compared to nearest-neighbor or linear methods, RBF interpolation offers higher accuracy and smoother field transfer without requiring mesh alignment or auxiliary mesh constructions [15]. This enables accurate interpolation of stress values from the fluid mesh \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma ^f$$\end{document} to the structural mesh \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma ^s$$\end{document} . Each interface mesh is defined by a set of vertices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{ \textbf{x}_i \in \mathbb {R}^d \, | \, i = 1, \dots , N \}$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \in \{2,3\}$$\end{document} spatial dimensions. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\boldsymbol{\sigma }_i \in \mathbb {R}^m \, | \, i = 1, \dots , N\}$$\end{document} denote the stress tensor or vector values (with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m = d$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m = d(d+1)/2$$\end{document} depending on tensor representation) at the fluid mesh vertices. For RBF interpolation, each stress component is treated as a scalar field. Given the k-th stress component \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\sigma _{i,k} \in \mathbb {R} \, | \, i = 1, \dots , N\}$$\end{document} at the input vertices, the interpolant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {I}_{\sigma _k}: \mathbb {R}^d \rightarrow \mathbb {R}$$\end{document} for the k-th component is defined as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {I}_{\sigma _k}(\textbf{x}) = \sum _{i=1}^N \lambda _{i,k} \, \phi \left( \left\| \textbf{x} - \textbf{x}_i \right\| \right) + \beta _{0,k} + \boldsymbol{\beta }_{l,k}^T \textbf{x}, \end{aligned}$$\end{document}

where N is the number of input vertices, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{i,k} \in \mathbb {R}$$\end{document} are the RBF coefficients, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi : \mathbb {R} \rightarrow \mathbb {R}$$\end{document} is the radial basis function, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{0,k} \in \mathbb {R}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\beta }_{l,k} \in \mathbb {R}^d$$\end{document} are coefficients of a linear polynomial ensuring consistency for linear stress fields. Specifically, we employ Compact Polynomial Radial Basis Functions (Compact Poly-RBFs), which provide excellent smoothness properties with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^6$$\end{document} continuity. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^6$$\end{document} continuity of these basis functions ensures that higher-order derivatives are continuous, thereby avoiding spurious oscillations that can occur with globally supported RBF kernels. This property is particularly important for transferring sensitive quantities like stress components, where ill-conditioned interpolation matrices can produce non-physical oscillations. The interpolation conditions are enforced as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {I}_{\sigma _k}(\textbf{x}_i) = \sigma _{i,k}, \quad i = 1, \dots , N, \end{aligned}$$\end{document}

with regularization constraints to guarantee uniqueness:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{i=1}^N \lambda _{i,k} = 0 \quad \text {and} \quad \sum _{i=1}^N \lambda _{i,k} \textbf{x}_i = \textbf{0}. \end{aligned}$$\end{document}

• Transfer of forces onto the structure: The process is repeated for all stress components via preCICE. The interpolant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {I}$$\end{document} can yield a matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}^{f\rightarrow s}$$\end{document} defining a linear mapping of the input forces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{f}^f$$\end{document} from the fluid mesh to forces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{f}^s$$\end{document} in the structural mesh:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{f}^s = \mathcal {M}^{f\rightarrow s}\textbf{f}^f \quad \text {with} \quad \textbf{f}^f= \begin{pmatrix} \textbf{f}^f(\textbf{x}_1) \\ \textbf{f}^f(\textbf{x}_2) \\ \vdots \\ \textbf{f}^f(\textbf{x}_N) \end{pmatrix} \end{aligned}$$\end{document}

When the row-sum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}^{f\rightarrow s}$$\end{document} is equal to one, the interpolation is consistent, which means that constant values are interpolated exactly. On the other hand, when the column-sum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}^{f\rightarrow s}$$\end{document} equals one the interpolation is conservative.

• Transfer of deformation onto the fluid: The computed structural deformation has to be transferred back to the fluid in order to account for the mesh deformation. Once again, this procedure is straightforward due to the use of the RBF method in preCICE. The deformation is computed at every quadrature point of the structural element:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{d}^f = \mathcal {M}^{s\rightarrow f}\textbf{d}^s \quad \text {with} \quad \textbf{d}^s= \begin{pmatrix} \textbf{d}^s(\textbf{x}_1) \\ \textbf{d}^s(\textbf{x}_2) \\ \vdots \\ \textbf{d}^s(\textbf{x}_n) \end{pmatrix} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^s$$\end{document} represents the structural displacements at discrete structural points \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{x}_i$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^f$$\end{document} represents the fluid mesh displacements at discrete fluid mesh points. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {M}^{s\rightarrow f}$$\end{document} represents the interpolation operator, which maps the structural displacement onto the fluid mesh on the IGA mesh cell level. Although vertex-based methods are widely used due to their simplicity and compatibility with standard finite element frameworks, they face several significant challenges when dealing with complex interface dynamics [10, 11]:

1. Geometric accuracy issues, especially for complex or curved interfaces [32].
2. Interpolation errors that can lead to numerical inaccuracies in force and displacement transfer processes.
3. Significant increase in communication overhead as interface meshes are refined. These limitations motivate the development of more accurate and efficient spline-based coupling methods, which are introduced in this paper.

Spline-based coupling method

The spline-based coupling method leverages analytical spline functions to represent the interface, providing a more efficient and accurate approach for data exchange between solvers. Unlike vertex-based methods that operate on discrete points, spline-based coupling utilizes a continuous mathematical description of both interface geometry and field variables.

We propose two types of spline-based coupling methods:

• Spline-vertex coupling method: one solver is IGA-based.
• IGA-IGA coupling method: both solvers are IGA-based.

Spline-vertex coupling method

We first consider the spline-vertex coupling method, which is designed for scenarios in which only one side of the interface (typically the structure) is modeled using IGA, while the opposing side (e.g., the fluid) employs a standard finite element or finite volume discretization. In this case, the interface geometry is represented using a spline on the IGA side, and data is transferred between the spline surface and the vertex points of the non-IGA mesh using projection or interpolation techniques (Fig. 2).Fig. 2. Spline-based coupling between IGA-based structural solver and vertex-based fluid solver, with interface data exchanged via continuous spline representation

The coupling procedure consists of two key steps:

• (a) Force transfer (fluid to structure): Force data from the fluid solver is transferred to the structure using a RBF interpolation method:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{f}^s = \mathcal {M}^{f\rightarrow s}\textbf{f}^f \quad \text {with} \quad \textbf{f}^f= \begin{pmatrix} \textbf{f}^f(\textbf{x}_1) \\ \textbf{f}^f(\textbf{x}_2) \\ \vdots \\ \textbf{f}^f(\textbf{x}_N) \end{pmatrix} \end{aligned}$$\end{document}

This step maps the discrete fluid force data to the structural solver’s representation.

• (b) Displacement transfer (structure to fluid): The structural displacement field, computed in terms of control point displacements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}_i^s$$\end{document} , is evaluated at any interface point using:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{d}^s(\boldsymbol{\xi }) = \sum _{i=1}^{n_c} R_i^p(\boldsymbol{\xi }) \textbf{d}_i^s \end{aligned}$$\end{document}

The fluid solver samples this continuous field at required interface points to update its mesh.

IGA-IGA communication: optimal spline-based coupling

When both fluid and structural domains use isogeometric analysis or support spline representation (of solution/geometry), the IGA-IGA communication approach fully leverages spline representations on both sides of the interface, representing the most efficient coupling strategy within the spline-based family (Fig. 3).Fig. 3IGA-IGA communication with spline representations on both fluid and structural domains, maximizing the benefits of isogeometric analysis

The communication procedure involves:

• • Interface parameterization:* Both domains ideally share a common parametric representation for direct control point data exchange. When parameterizations differ, reparameterization through knot vector exchange and interpolation creates conforming interface splines. This mapping is computed once during initialization for stationary interfaces.
• • Stress transfer via control points:* Forces are directly expressed in terms of control point values. When reparameterization in step (a) introduces new nodes, force and stress information is evaluated at these newly inserted nodes to enable accurate data transfer:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{f}^s_c = \mathcal {T}_{f \rightarrow s} \textbf{f}^f_c \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{f}^s_c$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{f}^f_c$$\end{document} are the control point force vectors for the structural and fluid domains, respectively, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_{f \rightarrow s}$$\end{document} is the transformation operator established during reparameterization.

• • Displacement transfer via control points:* Similarly, structural displacements computed at control points transfer directly or through transformation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{d}^f_c = \mathcal {T}_{s \rightarrow f} \textbf{d}^s_c \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^s_c$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{d}^f_c$$\end{document} are the control point displacement vectors for the structural and fluid domains, respectively, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}_{s \rightarrow f}$$\end{document} is the transformation operator mapping between the two spline spaces. This approach provides minimal communication overhead, exact geometry representation, preservation of higher-order solution properties, enhanced numerical stability, and efficient handling of complex geometries. We will further discuss the advantage of reduced communication overhead in Sect. 5.

Coupling methods verification and validation

In this section, we present a comprehensive verification and validation of the spatial coupling methods described previously. The verification procedure employs both the constant load and nonlinear load test cases to assess the accuracy, efficiency, and stability of the proposed coupling approaches.

We first consider a vertical beam that is clamped at its base and subjected to a constant distributed load along its free edge. The setup ensures that the force transmission is uniform and can be analytically validated.

The problem setup, including geometry and loading conditions, is shown in Fig. 4. We compare the tip displacement computed by the standalone structural solver (reference solution) to the coupled solver solutions using vertex-based and spline-based coupling methods.Fig. 4. Constant load test setup for the vertical beam

The structural solver computes the initial beam deformation under the applied load. The force distribution is mapped to the fluid mesh using the vertex-based and spline-based methods described earlier. The reconstructed forces are transferred back to the structure to verify accuracy in communication. Finally, the resulting tip displacement is compared with the theoretical solution, which is the monolithic application of the force along the beam within one solver.Fig. 5. Left: Comparison of tip displacement between reference solution and different coupling methods. Right: Relative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} error of different coupling methods compared to the reference solution

Figure 5 compares the vertical beam tip displacements obtained via different coupling strategies against the monolithic reference solution (left), and the corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} -norm errors in the y- and z-displacement directions (right). All spline-based methods (spline-vertex and spline-spline) yield errors below machine precision ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{-15}$$\end{document} ), indicating perfect agreement with the reference. The vertex-vertex method achieves the same level of accuracy in this constant load scenario. This demonstrates the consistency of all tested coupling approaches.

To further validate the accuracy and robustness of the coupling methods, we designed a nonlinear load test case. In this scenario, the vertical beam is subjected to a nonlinear load \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}(\xi ) = -\xi \textbf{e}_z$$\end{document} that varies with parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} .Fig. 6. Nonlinear load test setup for the vertical beam. The nonlinear load \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}(\xi )$$\end{document} varies linearly with parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} is the normalized coordinate along the beam height, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi =0$$\end{document} at the fixed bottom end and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi =1$$\end{document} at the free top end

Figure 6 illustrates the nonlinear load test setup. In this test, we compare the reference solution (direct calculation by a monolithic solver) with solutions obtained using different coupling methods (vertex-based and spline-based) under identical load conditions.Fig. 7. Tip displacement versus time under nonlinear loading

Figure 7 shows the tip displacement over time under a nonlinear load. The spline-spline based, spline-vertex based, and vertex-vertex based coupling methods yield consistent results under complex loading conditions.

Communication overhead comparison and analysis

In partitioned multiphysics simulations, communication overhead is a critical metric for evaluating coupling method efficiency. We conducted both theoretical analysis and experimental validation of communication overhead for different coupling methods to assess the advantages of spline-based coupling methods compared to traditional vertex-based approaches.

Theoretical Communication Overhead

As described in Sect. 4.2, the communication overhead for spline-based coupling methods primarily depends on the number of control points and parametric representation, while vertex-based methods see overhead increase exponentially with mesh refinement and quadrature point quantities. Figure 8 illustrates this theoretical difference.

For vertex-based coupling using traditional Gauss quadrature with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p+1$$\end{document} quadrature points per parametric direction per element, the communication overhead O is proportional to the number of interface elements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{\text {element}}$$\end{document} , the polynomial degree p, the number of timesteps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t$$\end{document} , and the parametric dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{d}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} O_{\text {quad}} = N_t \cdot d \cdot n_{\text {element}} \cdot (p+1)^{\hat{d}} \end{aligned}$$\end{document}

For spline-based coupling, the total communication overhead is primarily governed by the number of control points, the physical field dimension d, and the polynomial degrees \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_i$$\end{document} along each parametric direction. Unlike mesh-based coupling (which requires per-node communication), spline-based coupling grows at a much slower rate:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} O_{\text {spline}} = N_t \cdot d \cdot \left( \prod _{i=1}^{\hat{d}} n_i \right) + \sum _{i=1}^{\hat{d}} (n_i + p_i + 1) \end{aligned}$$\end{document}

This expression consists of two main components:

• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t \cdot d \cdot \left( \prod _{i=1}^{\hat{d}} n_i \right) $$\end{document} : the total number of scalar values exchanged during time integration. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} denotes the number of control points along the i-th parametric direction (with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{d}$$\end{document} such directions in total), and d represents the number of coupled physical fields.
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{i=1}^{\hat{d}} (n_i + p_i + 1)$$\end{document} : the cumulative size of the knot vectors in each parametric direction. This geometric data is exchanged once during initialization and remains unchanged throughout the simulation. We use preCICE [8] to mediate the exchange of field data (e.g., displacement, velocity, pressure) between solvers. On the interface, both the geometry and field values are represented using B-splines or NURBS. The spline representation allows us to transmit compact and structured data, where:
• The knot vector (of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{i=1}^d (n_{i} + p_i + 1)$$\end{document} ) is transferred at the beginning of the simulation.
• The field values are exchanged at every coupling step through the control points. This comparison highlights a major advantage of spline-based coupling. As shown in Fig. 8, vertex-based methods lead to an exponential increase in communication overhead as mesh refinement and polynomial degree grow, making them computationally expensive for high-fidelity simulations. In contrast, spline-based methods maintain significantly lower communication costs, as their overhead scales more efficiently with refinement.Fig. 8. Comparison of the communication overhead between vertex-based (FEM) and spline-based (IGA) coupling. The horizontal axis represents the number of parameter subdivisions per direction, serving as an equivalent indicator of the unique knots in IGA while reflecting mesh resolution in FEM. The vertical axis shows the total communication overhead in kilobytes

The analysis shows that as the mesh becomes more refined, the difference in communication overhead between vertex-based and spline-based methods will grow significantly.

Experimental validation results

The experimental results confirm our theoretical predictions about communication overhead in different coupling methods. We measured both the actual data transfer volume during the initial timestep ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t=1$$\end{document} ) and total communication time per dimension over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t=10000$$\end{document} timestep in a vertical beam test case across various mesh refinement levels.

Our instrumentation specifically targeted only the coupling-related data transfer functions, ensuring that the measured values reflect purely the communication overhead of different coupling methods without contamination from other computational processes.

In the case of spline-based coupling, additional overhead arises from transmitting knot vectors alongside control point data. To support tensor-product B-spline representations in multiple parametric directions, knot vectors are embedded into a unified matrix structure with NaN padding. Since the knot vectors are defined separately for each parametric direction, we embed them into a unified matrix where each row corresponds to one direction. NaN values are used to fill the remaining entries. This ensures consistency in data formatting and facilitates reliable transfer and storage.

Here we consider a parametric dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{d} = 2$$\end{document} and a physical dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d = 3$$\end{document} . The tensor-product B-spline space is defined by two univariate knot vectors associated with each parametric direction:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Xi&= \{\xi _1, \xi _2, \ldots , \xi _{n_\xi + p_\xi + 1}\} \subset \mathbb {R}, \\ \mathscr {H}&= \{\eta _1, \eta _2, \ldots , \eta _{n_\eta + p_\eta + 1}\} \subset \mathbb {R}, \end{aligned}$$\end{document}

where each knot vector is a non-decreasing sequence of real numbers. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_\xi $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_\eta $$\end{document} denote the number of basis functions in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} - and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} -directions, respectively, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\xi $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\eta $$\end{document} are the corresponding polynomial degrees.

In the implementation, while arranging knot vectors into a unified matrix format to ensure consistent data transmission, we inevitably introduce some overhead.

Since the lengths of knot vectors can differ between parametric directions, NaN values are used as padding to maintain consistent matrix dimensions:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \texttt {knotMatrix} = \begin{bmatrix} \xi _1 & \xi _2 & \dots & \xi _{n_\xi + p_\xi +1} & \texttt {NaN} & \dots & \texttt {NaN} \\ \texttt {NaN} & \dots & \texttt {NaN} & \eta _1 & \eta _2 & \dots & \eta _{n_\eta + p_\eta + 1} \end{bmatrix} $$\end{document}

In the numerical experiment, we consider the symmetric case where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_\xi = n_\eta = n$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\xi = p_\eta = p$$\end{document} , so that both knot vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\Xi }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {H}$$\end{document} have the same length, namely \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n + p + 1$$\end{document} . To embed them into a unified matrix for transfer, we allocate a matrix of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \times (2n + 2p + 2)$$\end{document} , where each row stores one knot vector padded with NaN values in the complementary half. Consequently, the total number of NaN entries equals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2(n + p + 1)$$\end{document} .

All timing data were collected on a single compute node of the DelftBlue supercomputer [33]. Given that runtime performance can be affected by external processes and background system noise, all measurements were repeated ten times in an isolated environment with minimal external interference. The values reported represent the averaged results, ensuring reliability and minimizing the influence of transient fluctuations.Fig. 9. Experimental (Exp.) and theoretical (An.) comparison of communication overhead with respect to the number of parameter subdivisions r per direction, with for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\in \{1,2,...,7\}$$\end{document}

Figure 9 presents a comprehensive analysis of communication overhead for different coupling methods through four complementary perspectives. In Fig. 9a

Figure 9b provides a similar data volume comparison for spline-based coupling, where both theoretical and experimental results confirm the significantly reduced growth rate in communication overhead as mesh resolution increases. Finally, Fig. 9c translates these data volume differences into actual communication time measurements, validating that the theoretical efficiency gains of spline-based coupling translate to meaningful performance improvements in practice, with timing differences becoming increasingly pronounced at finer mesh resolutions.

As shown in Fig. 9b

Table 1 supports this observation by listing the absolute differences between experimental and analytical spline-based communication overheads at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=7$$\end{document} for various polynomial degrees p, corresponding to the results shown in Fig. 9b

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varepsilon _{\text {absolute}} = \left| O_{\text {exp}} - O_{\text {an}} \right| , \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O_{\text {exp}}$$\end{document} denotes the measured (experimental) communication overhead in kilobytes, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O_{\text {an}}$$\end{document} is the corresponding analytical estimate based on Eq. 18.Table 1. Absolute discrepancy (KB) between measured and analytical communication sizes at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=7$$\end{document} , with corresponding analytical number of NaN-padded entriespAbsolute discrepancy (KB)# Measured NaN entries# Analytical NaN entries^^22.078125266.0026632.1096875270.0427042.140375273.9727452.1721875278.04278Computed as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2(n + p + 1)$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2^r$$\end{document} is the number of basis functions per parametric direction and p is the spline degree

Key observations from comparing the experimental results (Fig. 9) with theoretical results (Fig. 8) include:

1. The communication overhead for Vertex-based coupling methods increases rapidly with mesh refinement in practice, which is consistent with the theoretical prediction.
2. Spline-based coupling methods (including spline-vertex and IGA-IGA approaches) demonstrate significantly lower growth rates in communication overhead, confirming the theoretical advantages.
3. The experimental curves closely match the theoretical predictions, although the absolute values differ due to implementation details and overhead from the coupling library. These results confirm the communication efficiency advantages of spline-based coupling methods, particularly when dealing with high-resolution interface meshes. The close agreement between theoretical predictions and experimental measurements validates our communication overhead model and provides a solid foundation for estimating performance in larger, more complex simulations.

Benchmark problems

So far we have demonstrated the superior communication efficiency of spline-based coupling methods in solvers that support spline representation, while also illustrating how vertex-based coupling methods can effectively interface with solvers utilizing alternative discretization techniques through the preCICE library. In this section, we extend our analysis to validate the accuracy and robustness of the proposed methods through a series of carefully selected benchmark problems that progressively increase in complexity and physical relevance. The solid domain is modeled using the G+Smo library [31], which provides native support for IGA, while the fluid domain is solved using OpenFOAM[34], a widely used finite volume-based CFD framework.

Partitioned heat conduction

In this example, we employed the spline-based communication method to simulate heat transfer across partitioned domains. Beyond the previously demonstrated advantages in communication efficiency and information preservation, this application reveals two additional benefits particularly relevant to FSI. Let’s first start with the benchmark setup.

The partitioned heat conduction setup, shown in Fig. 10, consists of two square regions with different thermal properties. The governing equation is the dimensionless, transient heat equation in two spatial dimensions:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \frac{\partial u}{\partial t}&= \nabla ^2 u + f & \text {in } \Omega = [0,2] \times [0,1],\quad t \in [0,T], \\ u&= g & \text {on } \partial \Omega ,\quad t \in [0,T], \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ u(x, y, t) $$\end{document} is the temperature field, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f $$\end{document} is a source term, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ g $$\end{document} denotes prescribed boundary data.

To validate the numerical method, we adopt a manufactured solution of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u_{\text {exact}}(x, y, t) = 1 + x^2 + \alpha y^2 + \beta t, \end{aligned}$$\end{document}

with constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \alpha , \beta \in \mathbb {R} $$\end{document} . Substituting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ u_{\text {exact}} $$\end{document} into (20) yields the corresponding source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f $$\end{document} and Dirichlet boundary condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ g $$\end{document} .

Following the formulation in [35], the domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Omega $$\end{document} is partitioned into two non-overlapping subdomains:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Omega _D = [0,1] \times [0,1], \qquad \Omega _N = [1,2] \times [0,1], $$\end{document}

such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Omega = \Omega _D \cup \Omega _N $$\end{document} . These subdomains are coupled along the shared interface

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Gamma = \Omega _D \cap \Omega _N = \{1\} \times [0,1], $$\end{document}

as shown in Fig. 10. This configuration is referred to as the partitioned heat equation.

During the coupling procedure, either the temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ u $$\end{document} or the normal component of the heat flux

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q_n:= \nabla u \cdot \textbf{n} $$\end{document}

is exchanged across \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Gamma $$\end{document} . Specifically, one subdomain applies a Dirichlet condition

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ u = u_\Gamma \quad \text {on } \Gamma , $$\end{document}

while the other enforces a Neumann condition

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \nabla u \cdot \textbf{n} = q_\Gamma \quad \text {on } \Gamma , $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ u_\Gamma $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q_\Gamma $$\end{document} represent the interface temperature and heat flux, respectively.Fig. 10. Partitioned heat conduction problem setup

Naturally non-conforming support: One of the main advantages of the proposed coupling strategy is its natural ability to handle non-conforming discretizations at the coupling interface. Even when the mesh is geometrically discontinuous across the interface, the use of spline-based representations allows us to preserve functional smoothness. In the current implementation, non-conforming meshes are assumed to be nested, meaning that the control points or knot vectors of one mesh are fully contained within those of the other.Fig. 11. Temperature field of a partitioned spline-mesh domain with varying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _N$$\end{document} refinement levels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{N}$$\end{document}

As illustrated in Fig. 11, the Dirichlet subdomain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _D$$\end{document} is discretized with a coarser mesh refined once in Fig. 11a and twice in Fig. 11c and Fig. 11e

Inherited derivatives via spline representation: One of the distinctive advantages of using spline-based representations for coupling is that the derivatives of the solution field are inherently preserved and transferred across the interface. Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _D$$\end{document} is represented as a spline function, its smoothness and continuity naturally extend to its derivatives. As a result, when exchanging data across non-conforming interfaces, both the field values and their gradients are consistently evaluated from the same spline representation, without the need for additional reconstruction or projection.Fig. 12. Validation of derivative transfer via spline representation, relative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} -error of the heat flux remains at machine precision level

As illustrated in Fig. 12, the spline-based temperature field received from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _D$$\end{document} participant on the interface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _D$$\end{document} is differentiated to evaluate the heat flux. These computed flux values are then compared against the exact analytical solution. The red curve presents the relative error at each time step on a logarithmic scale, while the blue dashed line denotes the mean error across all time steps, which remains consistently around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6.7 \times 10^{-15}$$\end{document} . This result demonstrates the high-fidelity derivative that can be accurately evaluated from the spline field without needing to transfer the derivative separately.

Perpendicular flap

The perpendicular flap benchmark is a widely used FSI test case [8] provided by preCICE to evaluate the performance and accuracy of partitioned coupling strategies between fluid and solid solvers. This benchmark is particularly useful for assessing solver interoperability and conducting a mesh convergence study.

We model a three-dimensional flow that interacts with an elastic flap fixed perpendicularly to the channel floor. The physical parameters for both the fluid and solid domains are detailed in Fig. 13. We employ vertex-vertex coupling with data exchange at quadrature points. The structure uses IGA discretization (single spline patch) while the fluid mesh is a structured grid with local refinement near the flap interface. Convergence results for structural refinement levels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\in \{1,2,3\}$$\end{document} and polynomial orders \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=3,4$$\end{document} are shown in Fig. 14.Fig. 13. Setup and parameters of the perpendicular flap test case. The original image is taken from the preCICE paper [8]

Figure 15 illustrates the time evolution of the flap tip displacement under fluid loading with time step size. The flap exhibits characteristic oscillatory behavior due to the interplay between fluid pressure and structural elasticity.Fig. 14. Coupled fluid–structure interaction at selected time instances: fluid velocity magnitude (background) and structural pressure distribution (foreground) at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0.3\,\text {s}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0.7\,\text {s}$$\end{document}

In Fig. 15, the dark blue curve corresponds to G+Smo, which is the only solver in this comparison based on IGA. It shows strong agreement with other well-established solvers in both the amplitude and phase of the oscillatory motion of the flap tip. This comparison is based on the tip displacement in the x-direction, using a time step size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t = 0.01\,\text {s}$$\end{document} .Fig. 15. Displacement of the perpendicular flap over time under fluid loading

To demonstrate applicability to large-scale systems, we extended the benchmark to a fluid domain of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20 \text {m}$$\end{document} width containing 20 perpendicular flaps, each spaced at a center-to-center distance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.05 \text {m}$$\end{document} . All other parameters are consistent with the perpendicular flap benchmark 13. This configuration employs vertex-based coupling. Such multi-flap systems represent practical applications, including arrays of flow control systems or flexible marine structures [18], where decomposition of the system into subdomains for efficient parallel computing across numerous interface components becomes critical for overall performance (Fig. 16).Fig. 16. Multi-field simulation of a system with 20 perpendicular flaps under fluid loading

The benchmark has also been extended to a pure IGA setting, where both the fluid and structural solvers are based on G+Smo (specifically, the gsIncompressibleFlow [36] and gsElasticity modules) and the information is exchanged through the spline-based communication method in a partitioned way coupled through preCICE. In this fully IGA framework, the fluid domain motion is handled using an Arbitrary Lagrangian–Eulerian (ALE) formulation. Specifically, we employ the Tangential Incremental Nonlinear Elasticity (TINE) method [37] to compute the mesh deformation, ensuring robustness under large interface displacements. For detailed comparison and advantages of this method, the authors kindly ask the readers refer to the paper by A. Shamanskiy [37] for more details. Figure 17 presents the simulation results at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.3\,\text {s}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.7\,\text {s}$$\end{document} , depicting both the fluid velocity field and the isogeometric mesh preserved by the TINE algorithm.Fig. 17. Visual simulation results of IGA-IGA coupling: fluid velocity magnitude and isogeometric mesh at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.3\,\text {s}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.7\,\text {s}$$\end{document}

To verify the proposed pure IGA coupling method, we benchmark the G+Smo-G+Smo framework against the hybrid OpenFOAM(v2404)-G+Smo configuration. Figure 18a presents the tip displacement comparison, showing that the pure IGA coupling yields results consistent with the vertex-based hybrid approach. Furthermore, Fig. 18b highlights the efficiency of the IGA-based method, which demonstrates significantly lower communication overhead compared to the vertex-based strategy.Fig. 18. Perpendicular flap benchmark in pure IGA setting, information exchanged through the spline based communication method in a partitioned way coupled through preCICE

Turek Hron FSI benchmark

The third benchmark involves two-dimensional flow-induced oscillations of an elastic structure mounted behind a rigid cylinder, commonly referred to as the Turek-Hron FSI3 benchmark. This configuration is widely used to evaluate the performance of FSI solvers, particularly under challenging conditions characterized by strong coupling effects and significant added-mass phenomena.

The geometry of the Turek-Hron FSI3 benchmark is illustrated in Fig. 19. The configuration consists of a flexible beam attached to the downstream side of a fixed circular cylinder, located inside a rectangular channel of height \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.41\,\textrm{m}$$\end{document} and length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.5\,\textrm{m}$$\end{document} . The diameter of the cylinder is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.1\,\textrm{m}$$\end{document} , and the beam has a length of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.35\,\textrm{m}$$\end{document} and a thickness of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.02\,\textrm{m}$$\end{document} . The inflow boundary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\text {in}}$$\end{document} is prescribed with a parabolic velocity profile, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\text {in}}(y) = 1.5 \bar{v} \frac{4y(0.41 - y)}{0.41^2}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{v}$$\end{document} is the average inlet velocity. No-slip boundary conditions are applied along the top and bottom channel walls, denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\text {no-slip}}$$\end{document} . The average inlet velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{v} = 2\,m/s$$\end{document} .Fig. 19. Turek Hron FSI benchmark setup

The fluid is assumed to be Newtonian and incompressible, with a density of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _f = 1000\,\mathrm {kg/m^3}$$\end{document} and a kinematic viscosity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _f = 10^{-4}\,\mathrm {m^2/s}$$\end{document} , resulting in a Reynolds number of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{Re} = 200$$\end{document} . The fluid domain is discretized using an OpenFOAM structured mesh generated with blockMesh, composed of hexahedral cells organized in 18 blocks, consisting of approximately 46k cells. The structure is modelled using the Saint-Venant–Kirchhoff material with a Young’s modulus of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E = 5.6\,\textrm{MPa}$$\end{document} and Poisson’s ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _s = 0.4$$\end{document} . The solid density is taken as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _s = 1000\,\mathrm {kg/m^3}$$\end{document} , which matches the fluid density to ensure neutral buoyancy. Coupling employs the vertex-vertex strategy with data exchange at interface quadrature points. On the structural side, we use IGA discretization with spline bases of degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=2$$\end{document} and uniform refinements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\in \{3,4\}$$\end{document} at the interface; on the fluid side, the structured hexahedral mesh is refined near the cylinder and beam.Fig. 20. Simulation results for the Turek-Hron FSI3 benchmark using G+Smo coupled with OpenFOAM. (a) Vertical displacement of the beam tip over time. (b) Beam displacement and fluid pressure distribution. (c) Beam displacement and fluid velocity magnitude for mesh \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=4$$\end{document} with degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=2$$\end{document} at time step \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=5.9 \text { s}$$\end{document}

The left plot in Fig. 20 shows the time evolution of the vertical displacement at the tip of the elastic beam. The results obtained using G+Smo with two refinement levels ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 3$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 4$$\end{document} , both with degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 2$$\end{document} ) are compared against reference data. The numerical solution from the finer mesh ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 4$$\end{document} ) closely follows the reference in both amplitude and phase, while the coarser mesh ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 3$$\end{document} ) slightly underestimates the peak displacements but still captures the oscillation frequency well.

This strong agreement demonstrates that our coupling framework for the IGA-based structure solver is capable of accurately resolving complex FSI dynamics. It highlights the method’s ability to perform high-fidelity and efficient coupling within libraries that support splines, while also remaining compatible with solvers that work without spline-based boundary representations.

Conclusion

In this work, we have put together a seamless, start-to-finish spline coupling workflow for partitioned multiphysics simulations, with a particular focus on fluid–structure interaction problems. Our approach demonstrates that splines not only serve as mathematical representations for geometry and solution fields, but also function effectively as an efficient data exchange format between different solvers.

1. Development of spline-based coupling methods for enhanced representation of interface fields.
2. Comprehensive analysis of communication overhead, demonstrating the efficiency of spline-based methods over traditional vertex-based approaches. Our theoretical model and experimental validation show that communication overhead for vertex-based coupling grows rapidly with mesh refinement, while spline-based methods demonstrate significantly lower growth rates.
3. Implementation of robust coupling strategies that naturally handle non-conforming but nested meshes and preserve solution derivatives across interfaces. The spline-based approach eliminates the need for complex interpolation schemes when meshes do not match at the interface, and inherently preserves both solution values and their gradients, as demonstrated in our heat conduction benchmark where derivative information was accurately maintained without requiring separate transfer operations. The proposed methods are verified using vertical beam test cases under both constant and nonlinear loading, and validated through established FSI benchmarks, including the perpendicular flap and Turek-Hron configurations. Although all benchmarks are two-dimensional, this restriction does not limit the generality of our findings. The primary goal is to assess the coupling strategy and its communication behavior, both of which extend naturally to three-dimensional settings.

Future work will focus on extending the proposed methods to more complex multiphysics problems, including those involving large deformations, topology changes, multiple coupled fields, and higher Reynolds number flows. While the current benchmarks demonstrate the effectiveness of spline-based coupling for established FSI test cases, future work will address higher Reynolds number regimes where steeper interface gradients and more complex flow physics pose additional challenges for interface data transfer. The smooth, high-order nature of spline representations should provide advantages in such demanding applications. In addition, we plan to implement adaptive local refinement using Truncated Hierarchical B-splines (THB-splines) [39], which enable targeted mesh refinement in regions with high solution gradients while preserving the advantages of spline-based coupling in 3D geometries. Additionally, we intend to adapt our coupling methodology to Kirchhoff–Love shell structures [40, 41], which are prevalent in many FSI applications such as aerospace and biomedical engineering. Another promising direction is the development of IGA-based mesh transfer techniques specifically optimized for FSI applications, which could further improve the handling of moving and deforming interfaces [42]. Finally, we will continue optimizing our implementation to support large-scale, industrial-grade simulations.