Partitioned Co-Simulation for CAD-integrated Vibroacoustic Problems in Unbounded Domains Juan Ignacio Camarottia,1,∗, Philip Leb,d,1 , Yinshan Caic,d , Ricky Aristioa , Dionysios Panagiotopoulose,d , Elke Deckersb,d , Roland Wüchnera a Chair of Structural Analysis, Technical University of Munich, Arcisstr. 21, 80333 München, Germany b Department of Mechanical Engineering, Campus Diepenbeek, KU Leuven, Wetenschapspark 27, Diepenbeek, B-3590, Belgium c Department of Mechanical Engineering, KU Leuven, Celestijnenlaan 300, B-3001, Heverlee, Belgium d FlandersMake @KU Leuven, Leuven, Belgium
arXiv:2609.15384v1 [cs.CE] 14 Sep 2026
e Department of Mechanical Engineering, Campus De Nayer, KU Leuven, Jan Pieter de Nayerlaan 5, B-2680, St. Katelijne Waver, Belgium
Abstract Vibroacoustic analysis often requires coupling structural and acoustic solvers based on different numerical formulations and discretizations, making monolithic implementations intrusive and limiting software modularity and reuse. This work presents a partitioned co-simulation framework for exterior vibroacoustic analysis that couples an Isogeometric boundary representation analysis (IBRA) structural solver with an isogeometric boundary element method (IGA-BEM) acoustic solver. The methodology operates directly on the computer-aided design (CAD) boundary representation, preserving the exact geometry throughout the analysis and supporting both weak and strong coupling between non-conforming discretizations. A key contribution is the extension of the Aitken dynamic relaxation and Interface Quasi-Newton with Inverse Least-Squares (IQN-ILS) convergence accelerators to complex-valued interface quantities, allowing the coupling iterations to account directly for both amplitude and phase information. The approach is validated using one-way and two-way coupled vibroacoustic benchmark problems involving thin-shell structures and exterior acoustic domains. The results show excellent agreement with monolithic reference solutions, while the proposed complex-valued convergence accelerators improve the robustness and convergence behavior of the strongly coupled solution procedure without compromising solution accuracy. These results demonstrate that the proposed approach provides an accurate, robust, and modular approach for CAD-integrated frequency-domain vibroacoustic analysis. Keywords: Vibroacoustic, Co-simulation, Boundary Element Method (BEM), Isogeometric Analysis (IGA), Isogeometric B-Rep Analysis (IBRA), Complex-valued Convergence Accelerators
1. Introduction The control of structural vibration and noise radiation has become an increasingly important aspect of modern engineering design in sectors such as automotive, aerospace, transportation, and renewable energy. The widespread adoption of lightweight structures further increases the importance of accurately predicting structural–acoustic interactions, as these systems are generally more susceptible to vibroacoustic excitation. Consequently, reliable computational tools have become essential for evaluating and optimizing the dynamic and acoustic performance of engineering systems during the virtual design stage, reducing the reliance on costly experimental prototypes. In many industrial applications, frequency-domain analyses are of particular interest because they enable the characterization ∗ Corresponding author
Email addresses: [email protected] (Juan Ignacio Camarotti), [email protected] (Philip Le), [email protected] (Yinshan Cai), [email protected] (Ricky Aristio), [email protected] (Dionysios Panagiotopoulos), [email protected] (Elke Deckers), [email protected] (Roland Wüchner) 1 Juan Ignacio Camarotti and Philip Le share first authorship and contributed equally to this work. Preprint submitted to Engineering with Computers
September 15, 2026
of steady-state structural and acoustic responses over broad frequency ranges, providing valuable information for noise, vibration, and harshness (NVH) assessment and product optimization. For structural dynamics and bounded acoustic domains, the Finite Element Method (FEM) [1, 2, 3] provides a flexible and widely adopted numerical method. However, exterior acoustic problems require special treatment of the unbounded domain, typically through domain truncation techniques such as infinite elements [4, 5] or perfectly matched layers (PML) [6]. In contrast, the Boundary Element Method (BEM) [7, 8] inherently satisfies the Sommerfeld radiation condition while reducing the problem dimensionality to the boundary. Therefore, coupled FEM–BEM formulations have become an established approach for exterior structural-acoustic interaction problems. These formulations are, however, commonly implemented using monolithic solution strategies, where the structural and acoustic systems are assembled into a single coupled system. In parallel with these developments, significant advances have also been achieved in CAD-integrated numerical discretizations. Isogeometric Analysis (IGA) [9] and its extension to Isogeometric B-Rep Analysis (IBRA) [10, 11, 12] enable simulation directly on spline-based CAD boundary representations, avoiding geometry reconstruction while preserving exact geometric representation throughout the analysis process. Compared with conventional finite and boundary element discretizations, isogeometric formulations such as IGA-FEM and IGA-BEM [13, 14] provide highly smooth approximation spaces and reduce pollution effects [15], which are particularly advantageous for wave propagation problems. The combination of CAD-based structural discretizations with IGA-BEM is especially attractive for exterior acoustic problems, since the boundary discretization required by BEM is directly available from the CAD model. Several monolithic isogeometric structural acoustic formulations have therefore been developed for frequency-domain analyses [16, 17, 18]. Despite these advances, most frequency-domain vibroacoustic formulations remain based on monolithic coupling strategies [19, 20], in which the structural and acoustic governing equations are assembled into a single coupled system. Although such approaches provide accurate and mathematically consistent solutions, they require intrusive integration of the participating solvers and limit the reuse of independently developed simulation software. This restriction becomes particularly relevant in industrial workflows, where structural and acoustic analyses are often performed using dedicated commercial or proprietary solvers with different numerical formulations, meshes, and non-conforming interface discretizations. Since these solvers are often available only as executable software without access to their source code, intrusive modifications required by monolithic coupling strategies are generally infeasible. Integrating such heterogeneous software within a unified monolithic approach therefore requires substantial implementation effort and reduces software modularity, interoperability, and maintainability. Partitioned approaches address these limitations by decomposing the coupled problem into individual physical subproblems, each governed by its own discretization and solution algorithm, while enforcing the coupling conditions through the exchange of interface quantities across the coupling interface. When these subproblems are solved using independent simulation codes whose execution, synchronization, and communication are coordinated through a dedicated coupling tool, the partitioned methodology is commonly referred to as co-simulation [21, 22, 23]. Such a framework preserves the independence of the participating solvers while coordinating their execution and communication throughout the coupled analysis. Several partitioned structural–acoustic formulations have been investigated in recent years. Rodrı́guez-Tembleque et al. [24] developed a frequency-domain FEM-BEM partitioned formulation using mortar and localized Lagrange multiplier coupling strategies for non-matching interfaces. Bunting and Miller [25] introduced a staggered coupling approach for time-domain structural acoustics, while Kersschot et al. [26] demonstrated the applicability of partitioned co-simulation for time-domain vibroacoustic systems by coupling a linearized Euler flow-acoustic solver with a structural solver through the preCICE communication library [27]. Although these studies demonstrate the potential of partitioned methodologies, the integration of CAD-integrated discretizations within a co-simulation framework for frequency-domain vibroacoustic analysis has received no attention. Moreover, the iterative solution of strongly coupled partitioned problems poses additional challenges related to the robustness and efficiency of the coupling procedure, making effective convergence acceleration essential. A key ingredient of strongly coupled partitioned simulations is the use of convergence acceleration techniques to stabilize and accelerate the fixed-point iterations enforcing the coupling conditions. These methods can generally be classified into two groups [28, 29]. The first comprises relaxation-based methods, which improve convergence by applying a global scaling factor to the interface correction, the most widely used representative being the Aitken dynamic relaxation method [28, 30]. The second comprises quasi-Newton methods, which exploit information from previous coupling iterations to approximate the interface Jacobian and construct improved interface corrections. Among these 2
approaches, the Interface Quasi-Newton with Inverse Least-Squares (IQN-ILS) method [29, 31] has become one of the most widely adopted techniques in partitioned multiphysics simulations. In contrast to time-domain formulations, frequency-domain vibroacoustic analyses involve complex-valued interface quantities that simultaneously contain amplitude and phase information. However, existing implementations of Aitken, IQN-ILS, and related convergence accelerators have been developed primarily for real-valued interface vectors, requiring the real and imaginary parts to be treated separately. Such a treatment neglects the inherent complex structure of the coupled problem and may deteriorate convergence or increase the number of coupling iterations in strongly coupled vibroacoustic simulations, particularly in the vicinity of structural resonances. In that context, this work proposes a co-simulation framework for the partitioned solution of CAD-integrated frequency-domain exterior vibroacoustic problems, combining an IBRA structural formulation with an IGA-BEM acoustic solver. The proposed framework enables the non-intrusive coupling of independently developed structural and acoustic solvers while supporting non-conforming spline-based interface discretizations through nearest-neighbor, nearest-element, and mortar-based mapping strategies. Furthermore, complex-valued extensions of the Aitken dynamic relaxation and the IQN-ILS convergence accelerators are introduced to improve the robustness and efficiency of strongly coupled frequency-domain simulations. The proposed methodology is validated through representative one-way and two-way coupled problems and compared against corresponding monolithic and analytical reference solutions. The remainder of this paper is organized as follows. Section 2 presents the governing equations of the coupled frequency-domain vibroacoustic problem, including the structural and acoustic formulations, coupling conditions, and monolithic and partitioned solution strategies. Section 3 describes the proposed co-simulation framework, including the coupling algorithms, and the communication strategy between the participating solvers. Section 4 addresses the treatment of the coupling interface, presenting the considered mapping techniques for non-conforming discretizations and the convergence acceleration strategies employed for the strongly coupled partitioned solution. Section 5 assesses the proposed methodology through representative one-way and two-way coupled benchmark problems and compares the obtained results against corresponding monolithic reference solutions. Finally, Section 6 summarizes the main findings and concludes the paper. 2. Problem formulation This section introduces the mathematical formulation of the coupled vibroacoustic problem considered in this work. First, the FEM and indirect BEM discretized governing equations of the structural and acoustic subproblems are presented. Subsequently, the vibroacoustic coupling conditions are introduced and the corresponding monolithic and partitioned solution strategies are discussed, providing the foundation for the partitioned co-simulation framework developed in the following sections. 2.1. Structural model Assuming steady-state harmonic excitation at angular frequency ω, the structural response is governed by −ω2 M + iωC + K û(ω) = F̂m (ω) + F̂a (ω).
(1)
Here, M, C, and K ∈ RNs ×Ns denote the mass, damping, and stiffness matrices, respectively, while û(ω) ∈ CNs contains the complex-valued structural degrees of freedom. The right-hand side consists of a mechanical load contribution F̂m (ω) ∈ CNs and an acoustic load contribution F̂a (ω) ∈ CNs arising from the pressure acting on the vibroacoustic interface Γsa . The acoustic loading contribution is obtained from the pressure field computed by the acoustic solver and transferred to the structural interface degrees of freedom through an appropriate mapping procedure. In this paper, the structural domain is discretized using spline-based isogeometric shell elements operating directly on the CAD boundary representation, thereby preserving the exact geometry throughout the analysis. Both the Kirchhoff–Love shell formulation of Kiendl et al. [32] and the Reissner–Mindlin shell formulation proposed by Benson et al. [33] are considered. While the former requires at least C 1 -continuity of the displacement field both within patches and across patch interfaces, the latter requires C 0 across patch interfaces. Since the focus of the present work is the partitioned vibroacoustic coupling framework rather than the structural discretization itself, the interested reader is referred to [10, 32, 33] for details on the employed shell formulations. 3
Rayleigh damping is optionally considered in this paper, with C = αM + βK,
(2)
where α and β denote the Rayleigh damping coefficients. 2.2. Acoustic model Given the homogeneous Helmholtz equation, deploying the Green’s function G : R3 × R3 → C and considering two sides of the boundary surface Γ, the indirect boundary integral equation # Z " ∂G(r, rf ) − G(r, rf )σ̂(rf ) dΓ(rf ), r ∈ Ω \ Γ, (3) p̂a (r) = µ̂(rf ) ∂n Γa can be derived. In this equation, p̂a (r) denotes the acoustic pressure at the field point r in domain Ω and r f is the source point on Γ. In addition, the single σ(r f ) and double layer potential µ(r f ) are the difference of the normal pressure gradient and pressures between two sides of Γ, respectively. Within the coupled vibroacoustic problem, the structural and acoustic subdomains exchange information through the vibroacoustic interface Γsa . The normal structural velocity v̂n is communicated to the acoustic solver and imposed as a Neumann boundary condition, ∂ p̂a = −iωρa v̂n , r ∈ Γsa , (4) ∂n where ρa denotes the fluid density and ω ∈ Ψ is the angular frequency with Ψ := [ωmin , ωmax ]. The acoustic boundary is represented using spline-based isogeometric discretizations, ensuring an exact geometric description and a geometrically consistent interface with the structural model. Practical CAD models are typically composed of multiple NURBS patches with independent parameterizations and discretizations. Because neighboring patches are generally non-conforming, continuity across shared interfaces must be enforced through an appropriate coupling strategy. For acoustic analyses, at least C 0 -continuity is required throughout the computational domain. In this work, this requirement is fulfilled using the strong patch coupling technique for non-conforming NURBS surfaces proposed by Coox et al. [34]. Following a Galerkin discretization of the boundary integral equations, the acoustic problem is expressed as A(ω) x(ω) = b(ω),
ω ∈ Ψ,
(5)
where A : Ψ → CN×N , x : Ψ → CN contains the discrete unknowns associated with the layer potentials, and b : Ψ → CN . Further details on indirect IGA-BEM formulations for acoustic analysis can be found in [14]. After solving the acoustic problem, the interface loads resulting from the acoustic pressure field are transferred to the structural model through the interface mapping procedure and applied to the shell mid-surface. Further details on the transfer of the interface quantities and the resulting vibroacoustic coupling are provided in Section 2.3. 2.3. Vibroacoustic coupling conditions and solution strategies The vibroacoustic interaction between the structural and acoustic domains is enforced through coupling conditions prescribed on the common interface Γsa , ensuring consistency of kinematic and dynamic quantities in the frequency domain. These interface conditions can be enforced either within a monolithic formulation, where both fields are solved simultaneously, or within a partitioned framework, where independent structural and acoustic solvers exchange interface quantities. The kinematic coupling condition enforces continuity of the normal velocity at the interface. The normal particle velocity of the acoustic medium is prescribed by the structural motion according to v̂n (r, ω) = iω ûn (r, ω),
r ∈ Γsa ,
(6)
where ûn denotes the normal component of the structural displacement field and v̂n is imposed as a Neumann boundary condition on the acoustic problem. 4
The dynamic coupling condition enforces equilibrium of tractions at the vibroacoustic interface. For thin-walled structures, the structural loading arises from the acoustic pressure acting on the shell surface. The resulting acoustic traction is given by t̂a (r) = µ̂(r) n,
r ∈ Γsa ,
(7)
where n denotes the outward unit normal vector of the structural domain. The corresponding discrete load contribution reads Z Fa = µ̂(r) n Ni (r) dΓ, (8) Γsa
with Ni denoting the acoustic basis functions. Based on these coupling conditions, vibroacoustic systems can be solved using either monolithic or partitioned strategies. In a monolithic formulation, the structural and acoustic equations are assembled into a single coupled system and solved simultaneously. After spatial discretization, the frequency-domain problem can be written as # " #" # " As (ω) Csa (ω) F̂s (ω) û(ω) = . (9) Cas (ω) − ω21ρa Aa (ω) x̂(ω) − ω21ρa F̂a (ω) Here, As (ω) = −ω2 M + iωC + K denotes the structural dynamic stiffness operator introduced in Eq. (1), while Aa (ω) represents the acoustic operator. The matrices Csa and Cas are the discrete coupling operators enforcing the vibroacoustic interface conditions. Monolithic formulations provide a mathematically consistent reference solution in which compatibility and equilibrium are satisfied implicitly through the global system. Alternatively, partitioned strategies retain separate specialized solvers for the structural and acoustic subproblems and enforce the coupling through an iterative exchange of interface quantities. The coupled problem can then be written as As (ω)û(k) (ω) = F̂s (ω) + Csa x̂(k−1) (ω),
(10)
Aa (ω)x̂ (ω) = F̂a (ω) + Cas û (ω),
(11)
(k)
(k)
where (k) denotes the coupling iteration. A single exchange per frequency step yields a weakly coupled strategy, whereas iterative enforcement of the interface conditions results in a strongly coupled scheme. In practical vibroacoustic applications, the structural and acoustic interface discretizations are generally nonconforming. Consequently, for partitioned approaches, the exchanged interface quantities cannot be transferred directly and require dedicated mapping operators. Within the present framework, all interface data transfers are performed through the Kratos MappingApplication [35, 36]. Structural displacements are transferred from the structural to the acoustic interface, while acoustic forces are mapped back conservatively to the structural side. The employed mapping techniques are described in detail in Section 4.1. While monolithic approaches provide a valuable reference solution, partitioned formulations preserve solver modularity and enable the non-intrusive coupling of independent structural and acoustic solvers. This is particularly attractive for FE–BE vibroacoustic analyses, where specialized solvers can be reused without modification. The present work therefore adopts a partitioned co-simulation framework in which the structural and acoustic solvers remain fully independent and communicate exclusively through mapped interface quantities. 3. Co-simulation strategies for partitioned vibroacoustic systems This section presents the partitioned co-simulation framework adopted in the present work. First, the considered coupling algorithms and the associated vibroacoustic workflow are introduced. Subsequently, the remote-controlled co-simulation strategy and the communication mechanism between the structural and acoustic solvers are described. Particular attention is devoted to the transfer of interface quantities between non-conforming discretizations through dedicated mapping operators. Finally, the fixed-point formulation of the coupled problem and the employed complexvalued convergence acceleration techniques are presented. 5
3.1. Coupling algorithms and workflow In coupled multi-physics simulations, the participating subproblems interact through the exchange of interface quantities. The coupling algorithm defines how this interaction is resolved and, in particular, how interface conditions are satisfied across the coupled domains. Depending on the strength of the interaction and the desired accuracy, different coupling strategies can be employed. In the following, we distinguish between one-way (unidirectional) coupling and two-way (bidirectional) coupling. For two-way coupling, both weak (explicit) and strong (implicit) schemes are considered. In addition, the order in which interface data are exchanged and incorporated into each solver can be described by different communication patterns, most commonly Jacobi- and Gauss–Seidel-type coupling. These concepts form the basis of the partitioned vibroacoustic framework considered in this work and are discussed in detail in the following. 3.1.1. Coupling strategies In one-way coupling, the interaction is assumed to be unidirectional, meaning that one subproblem influences the other while the feedback effect is neglected. In this case, interface data are transferred from the driving subproblem to the driven subproblem, but no information is transferred back. One-way coupling is computationally efficient and can be sufficient when the neglected feedback remains small compared to the dominant interaction. In two-way coupled simulations, both subproblems influence each other through the bidirectional exchange of interface quantities. Weak coupling, also referred to as explicit, staggered, or loose coupling, performs a single exchange of interface data and solves each subproblem once per frequency step (see Figure 1(a)). The absence of an iterative enforcement of the interface conditions can reduce robustness and accuracy, particularly in strongly coupled configurations or when the coupled response is highly sensitive to the interface quantities. Strong coupling (Figure 1(b)), also referred to as implicit or iterative coupling, aims to enforce the interface conditions by iterating between the participating subproblems until convergence is achieved. In this setting, the coupled solution is obtained by minimizing an interface residual that measures the violation of interface equilibrium NΓ and/or compatibility. Let y(k) denote the interface quantities exchanged between the coupled subproblems at Γ ∈ C coupling iteration k, where NΓ is the number of interface degrees of freedom. The corresponding interface residual is NΓ denoted by r(k) Γ ∈ C , and convergence is achieved when its norm falls below the prescribed tolerance ε. Strong coupling generally improves stability and accuracy compared to weak coupling, particularly in problems with pronounced interaction between the participating fields, at the expense of increased computational cost due to multiple solver evaluations per coupling step. In practice, relaxation or acceleration techniques can be employed to reduce the number of coupling iterations required for convergence. fn
fn k ←k+1
P1
P1
q(1),n+1,k Γ
q(1),n+1 Γ
P2
P2
q(2),n+1,k Γ no
∥rmΓ ∥ ≤ ε yes
f n+1
f
(a) Weak (explicit) coupling.
(b) Strong (implicit) coupling.
n+1
Figure 1: Weak (explicit) and strong (implicit) coupling schemes using a Gauss–Seidel execution pattern.
Besides the choice of coupling strategy, the coupling procedure can also be characterized by the employed communication pattern [37, 38]. In a Jacobi-type communication pattern, each participant uses interface data from the 6
previous coupling iteration to advance its solution, such that the participants can be evaluated independently. In contrast, a Gauss–Seidel-type communication pattern uses the latest available interface information, meaning that the subproblems are executed sequentially and downstream solvers immediately incorporate the most recent updates provided by upstream solvers. Jacobi patterns can offer higher potential for parallel execution, whereas Gauss–Seidel patterns often exhibit improved convergence properties due to the use of up-to-date interface information [37, 38]. In the present work, one-way coupling as well as two-way weak and strong coupling strategies are investigated. For the two-way coupled simulations, a Gauss–Seidel communication pattern is employed. 3.1.2. Vibroacoustic coupling workflow The overall partitioned frequency-domain vibroacoustic workflow is illustrated in Figure 2. The previously introduced coupling strategies are applied through a staggered exchange of interface quantities between the structural and acoustic solvers. For a given excitation frequency, the isogeometric structural problem is solved in the open-source multiphysics framework Kratos Multiphysics [35, 36], while the acoustic Helmholtz problem is solved using a MATLAB-based isogeometric boundary element solver. The structural interface displacements are transferred to the acoustic solver and imposed as boundary conditions. The resulting acoustic pressures are converted into equivalent loads acting on the structure and transferred back to the structural solver. Since the structural and acoustic interfaces may be discretized independently, the exchanged quantities are transferred through the mapping operators described in Section 4.1. For weak coupling, the interface quantities are exchanged only once per frequency step. In contrast, strong coupling iteratively updates the exchanged quantities until the prescribed interface convergence criterion is satisfied.
Start Connect the solvers and initialize co-simulation Set frequency ωi and k = 0
Structural solve in Kratos Ks(ωi)uk+1 = fΓk
Map and export interface displacement uk+1 Γ
uk+1 Γ : Kratos → MATLAB
Map and export acoustic interface load fΓk+1
k ←k+1
i←i+1
Acoustic BEM solve in MATLAB H(ωi)pk+1 = Gvnk+1
fΓk+1: MATLAB → Kratos
Compute interface residual rk+1 Γ no
krk+1 Γ k≤ε no
yes
yes
last frequency?
End
Figure 2: Workflow of the partitioned frequency-domain vibroacoustic coupling strategy.
7
3.2. Integration of the acoustic solver through a remote-controlled approach Since the external acoustic solver is not natively accessible by the Kratos co-simulation framework, a dedicated integration mechanism is required to enable its participation in the coupled simulation. To this end, the acoustic solver is integrated through a remote-controlled solver wrapper following the concept introduced by Bucher [23]. In this approach, the external solver does not implement any coupling logic. Instead, it exposes a set of callback routines associated with different stages of the solution procedure, such as importing interface data, solving the acoustic problem, exporting interface quantities, and performing post-processing operations. These routines are registered during the initialization stage and subsequently invoked by the co-simulation controller according to the selected coupling algorithm. After the registration stage, the external acoustic solver relinquishes control to the co-simulation framework and remains idle until the execution of a registered operation is requested. The main advantage of this approach is the complete separation between the solver implementation and the coupling algorithm. As a result, the same acoustic solver wrapper can be used for both weakly and strongly coupled simulations without modifying the MATLAB code. Changes in the coupling strategy therefore only require modifications within the co-simulation framework rather than within the external solver implementation. Furthermore, the centralized management of the execution sequence reduces the risk of synchronization issues and deadlocks when compared to classical co-simulation approaches. A representative pseudo-code implementation of the remote-controlled solver wrapper used for the acoustic participant is provided in Appendix A. 3.3. Data exchange between the solvers Since the structural and acoustic solvers execute as independent processes, an interprocess communication (IPC) mechanism is required for the exchange of interface quantities. In the present work, this functionality is provided by the tool CoSimIO, which acts as a detached communication interface between the participating solvers. The resulting data exchange architecture is illustrated in Figure 3.
Structural solver Harmonic analysis Kratos StructuralMechanicsApplication
IPC (CoSimIO)
Coupling interface Kratos CoSimulationApplication
IPC (CoSimIO)
Acoustic solver BEM Helmholtz solver MATLAB
Figure 3: Data exchange architecture in the partitioned vibroacoustic co-simulation: interface quantities are exchanged between the structural harmonic solver in Kratos and the MATLAB-based acoustic BEM Helmholtz solver through CoSimIO.
In the present implementation, a file-based communication approach is employed. Interface quantities are exchanged through a dedicated communication directory, while synchronization between both participants is ensured through auxiliary control files. This guarantees that data are accessed only after being completely written by the corresponding solver and thereby prevents race conditions during the coupling procedure. The choice of a file-based IPC approach is motivated by its robustness, simplicity, and portability across different programming languages and execution environments. Although this strategy introduces a higher communication latency compared to shared-memory or socket-based approaches, the associated overhead remains acceptable for the frequency-domain analyses considered in this work. 4. Interface treatment In partitioned vibroacoustic simulations, the enforcement of the coupling conditions across non-conforming interfaces relies on two key ingredients. First, interface quantities must be transferred consistently between the generally different discretizations employed by the structural and acoustic solvers. Second, for strongly coupled problems, the resulting fixed-point iterations require robust convergence acceleration techniques to ensure stability and efficiency. 8
The present section therefore addresses both aspects of the interface treatment, namely the transfer of interface quantities between non-conforming discretizations and the convergence acceleration strategies employed in the partitioned solution procedure. 4.1. Data mapping between non-conforming discretizations In partitioned vibroacoustic co-simulation, the structural and acoustic solvers generally employ independent discretizations of the vibroacoustic interface. As a consequence, the corresponding interface meshes are in general non-conforming, differing in topology, resolution, and numerical representation. To enable a consistent exchange of interface quantities between the two solvers, a data transfer operator is required to map fields between the respective interface discretizations. Although the presented framework is not restricted to a specific discretization technology, the considered mapping procedures are particularly well suited for partitioned IGA–IGA coupling, where both structural and acoustic interfaces are represented by spline-based geometries. In such settings, independently parameterized interfaces and nonconforming discretizations naturally arise, making robust and geometrically consistent transfer operators especially important for the coupled analysis procedure. In contrast to standard FEM discretizations, where nodal coordinates are typically located on the physical interface geometry, the control points defining spline-based geometries generally do not lie on the geometry itself. As a consequence, direct control-point-to-control-point transfer strategies commonly employed in FEM–FEM coupling are, in general, not suitable for isogeometric interface coupling. From an algebraic perspective, the data transfer between the interface discretizations can be expressed as a linear mapping. Let uo ∈ CNo denote a vector collecting interface quantities defined on the origin discretization and ud ∈ CNd the corresponding quantities on the destination discretization, where No and Nd denote the dimensions of the origin and destination interface vectors, respectively. The mapping operation is written as ud = H uo ,
(12)
where H ∈ RNd ×No is the mapping operator whose entries depend on the geometric relationship between the two interface discretizations and on the chosen mapping strategy. For vector-valued interface quantities, such as structural displacements and forces, the mapping is applied independently to each component using the same mapping operator. To ensure consistency of the coupled problem, force quantities are transferred using a conservative mapping. Let Fd ∈ CNd denote a vector of forces defined on the destination interface. The corresponding forces on the origin discretization, Fo ∈ CNo , are obtained as Fo = HT Fd .
(13)
This choice guarantees that the virtual work associated with the interface quantities is preserved under the mapping, FTo uo = FTd ud ,
(14)
and is therefore essential for maintaining physical consistency in the partitioned vibroacoustic coupling. The overall partitioned vibroacoustic coupling framework together with the transfer of interface quantities between the structural and acoustic domains is illustrated in Figure 4. In the present framework, different mapping strategies are explored for the transfer of interface quantities between the acoustic and structural interfaces. The nearest-neighbor and nearest-element mappers perform pointwise transfers between interface integration points based on geometric proximity and projection operations on the underlying spline interface representations. In contrast, the mortar mapper formulates the transfer as a weak variational projection between the interface discretizations, leading to a controlpoint-to-control-point transfer operator assembled through numerical integration over the spline coupling interface. 4.1.1. Nearest-neighbor mapping In the nearest-neighbor mapping, each destination integration point is associated with the spatially closest integration point on the origin interface (Figure 5). Let xd , xo ∈ R3 denote points on the destination and origin interfaces, respectively, and let O represent the set of origin integration points. The mapped value ũ(xd ) is obtained from the value at the closest origin integration point, 9
Consistent mapping ua,R (ω) = Hus,R (ω) ua,I (ω) = Hus,I (ω) real and imaginary displacement transfer
Structural domain:
Acoustic domain:
Frequency-domain elastodynamics with IBRA
Helmholtz equation with IGA-BEM
• Kratos IgaApplication
• MATLAB IGA-BEM solver coupled through Python
−ω 2 Ms + jωCs + Ks us = Fs,R + jFs,I
A(ω) x(ω) = b(ω)
real and imaginary load transfer
Fs,R (ω) = HT Fa,R (ω) Fs,I (ω) = HT Fa,I (ω)
Conservative mapping
Figure 4: Overview of the proposed partitioned vibroacoustic co-simulation framework. The structural harmonic solver and the acoustic Helmholtz solver exchange real and imaginary interface quantities through consistent displacement transfer and conservative load mapping operators.
ũ(xd ) = u(xo ),
xo = arg min ∥xd − xi ∥ .
(15)
xi ∈O
From an algebraic perspective, the nearest-neighbor mapper defines a sparse binary transfer matrix H ∈ RNd ×No , where No and Nd denote the dimensions of the origin and destination interface vectors, respectively. Each destination degree of freedom is associated with exactly one origin integration point. As a result, each row of H contains a single nonzero entry equal to unity, corresponding to the nearest origin point, while all remaining entries vanish. The nearest-neighbor mapping is straightforward to implement, since the transfer is based solely on the spatial proximity of interface integration points. In practice, the closest-point search can be performed efficiently using spatial search structures such as k-d trees [39] or octrees [40]. However, the method does not account for the local element geometry or interpolation within the origin discretization. As a consequence, the mapping accuracy may deteriorate for coarse or highly non-conforming interfaces, particularly when transferring strongly varying interface fields. origin integration points destination integration points
1.0 xdi
Γd
0.7 0.4 0.1
xoj(i)
Γo −1.0 j(i) = arg minj xdi − xoj 2 ũdi = uoj(i)
z
−0.2 0.5 −0.5
0.0 x
0.0 0.5 1.0
−0.5
y
Figure 5: Nearest-neighbor mapping between non-matching origin and destination interfaces.
10
4.1.2. Nearest-element mapping In the nearest-element mapping (Figure 6), each destination integration point on the interface Γd is mapped by geometrically projecting it onto the origin interface Γo , which is discretized using isogeometric analysis. For a given destination integration point xd ∈ R3 , the closest origin surface patch is first identified and xd is projected onto this patch to obtain the corresponding parametric coordinates ξio , ηoi ∈ Ωo . The mapped quantity is then evaluated by interpolating the origin field at the projected location using the isogeometric shape functions of the associated origin patch. The mapped value can be expressed as ũ(xd ) =
ne X
N oj ξio , ηoi uoj ,
(16)
j=1
where N oj denote the shape functions of the origin discretization, uoj are the corresponding origin degrees of freedom, and ξio , ηoi are the parametric coordinates associated with the geometric projection of xd onto the origin surface. From an algebraic perspective, the nearest-element mapping defines a sparse transfer operator H ∈ RNd ×No , where No and Nd denote the dimensions of the origin and destination interface vectors, respectively. Each destination integration point is associated with the local basis functions of the projected origin element. In contrast to the nearestneighbor mapper, multiple nonzero entries may appear in each row of H, corresponding to the active shape functions of the origin patch at the projected parametric location. The values of these entries are given by the evaluated isogeometric basis functions N oj (ξio , ηoi ). Compared to the nearest-neighbor approach, the nearest-element mapping accounts for the interpolation properties of the origin discretization, resulting in a smoother and more accurate transfer of interface quantities. This is particularly advantageous for curved IGA interfaces and strongly non-conforming discretizations, where pointwise nearest-neighbor transfers may introduce interpolation errors and discontinuities. Origin parametric space
η
destination integration point projected point
(ξio , ηio )
1.000
ξ find (ξio , ηio )
xdi
Γd
0.663
projection
0.325 xoi
−0.012
Γo −1.0 xoi = P ΠΓo (xdi ) e ũdi = nj=1 Njo (ξio , ηio ) uoj
z
−0.350 0.5 −0.5
0.0
0.0 x
0.5
1.0
−0.5
y
Figure 6: Nearest Element (closest geometric projection) mapper for IGA: the destination integration point is projected onto the origin IGA patch to obtain (ξio , ηoi ), and the mapped value is evaluated as ũdi = No (ξio , ηoi )uo .
4.1.3. IGA-IGA Mortar Mapper In contrast to the nearest-element mapping, where the transfer is performed pointwise through geometric projection (closest-point projection), the mortar mapper formulates the interface transfer as a variational projection problem over the coupling interface. The objective is to determine a destination field that minimizes the mismatch with the origin field in an integral sense, thereby providing a consistent transfer between non-matching discretizations. The present formulation extends the mortar-based isogeometric mapping concepts proposed in [41] towards partitioned IGA–IGA vibroacoustic coupling involving non-conforming NURBS interfaces. 11
Let Γo and Γd denote the origin and destination interfaces, respectively. Furthermore, let uo : Γo ⊂ R3 → C and ud : Γd ⊂ R3 → C denote the complex-valued scalar interface fields on the origin and destination interfaces. The mapped destination field ud is obtained by minimizing the L2 -norm of the error between the origin and destination fields over the interface, Z 1 (ud (x) − uo (x))2 dΓ. (17) Ψ(ud ) = 2 Γ Enforcing stationarity of the functional leads to the weak form Z Nd (x) (ud (x) − uo (x)) dΓ = 0, Γ
(18)
where Nd denotes the shape functions of the destination interface. The origin and destination fields are approximated using their corresponding isogeometric basis functions, ud (x) ≈ uhd (x) = NTd (x)ud ,
(19)
uo (x) ≈ uho (x) = NTo (x)uo ,
(20)
where No (x) ∈ RNo and Nd (x) ∈ RNd collect the active isogeometric basis functions of the origin and destination interfaces, respectively, while uo ∈ CNo and ud ∈ CNd are the corresponding vectors of interface degrees of freedom. Substituting the discrete approximations into the weak form (Eq. 18) yields the mortar system Z Z Nd NTd dΓ ud = Nd NTo dΓ uo , (21) Γd
Γdo
which can be written in matrix form as MDD ud = MDO uo , Nd ×Nd
(22)
Nd ×No
where MDD ∈ R and MDO ∈ R are the mortar matrices resulting from the integration of the destination basis functions and the destination-origin basis function products, respectively. Assuming that the matrix MDD is invertible, the mapping operator is obtained as ud = HDO uo ,
HDO = M−1 DD MDO .
(23)
For the IGA–IGA mortar mapper, the evaluation of the coupling integrals requires the construction of a common integration domain between both interfaces (Γdo ). Since the origin and destination IGA patches generally possess different parametric discretizations, the knot spans of the destination interface are first mapped to the physical space and subsequently projected onto the parametric space of the origin interface. This procedure defines overlapping regions where the mortar integrals are evaluated consistently. Figure 7 illustrates this procedure. The destination knot spans are mapped from the destination parametric space to the destination physical surface and then projected onto the origin parametric space. The resulting projected regions are subdivided along the knot lines of the origin discretization, generating integration subdomains that are conforming with both parameterizations. Numerical quadrature points are then constructed within these subdomains in the origin parametric space. To evaluate the destination basis functions Nd , each quadrature point is mapped to the physical interface and subsequently projected back onto the destination parametric space. The corresponding origin and destination basis functions, No and Nd , are then evaluated at their respective parametric coordinates and used to assemble the mortar matrices MDD and MDO . 4.2. Fixed-point problem and convergence accelerators In strongly coupled partitioned vibroacoustic simulations, the interface equilibrium conditions are enforced iteratively through the exchange of interface quantities between the participating solvers. As a consequence, the coupled problem can be interpreted as a fixed-point problem defined on the vibroacoustic interface. Let Γ denote the coupling interface between the structural and acoustic domains. Furthermore, let uΓ : Γ ⊂ R3 → 3 C and FΓ : Γ ⊂ R3 → C3 denote the complex-valued interface displacement and load fields, respectively. For a 12
IGA-IGA Mortar Mapper - Definition of the Integration Domain
Origin Physical Space
Destination Physical Space
0.45
0.9
0.15 z
0.3
−0.15
−0.3
−0.45 1.0
−0.9 1.5 1.0 0.5 0.0 y −0.5 −1.0 −1.5
0.5 −1.0
−0.5
0.0 0.0 x
−1.5 −1.0 −0.5
y
−0.5
0.5 1.0
−1.0
0.0 0.5 x 1.0
projection to origin parameter space
1.5
z
projection to destination physical space
ηd
Origin parametric space
ηo
1.0
integration points for the mortar mapper
0.8
Destination parametric space
1.2
Ωo
Ωd
subdivided mapped rectangle
1.0
untrimmed knot span
0.6 0.8 trimming curves
0.4 0.6 0.2
mapped triangles
0.4
0.0 0.0
0.2
0.4
0.6
0.8
0.2 0.2
ξo
1.0
0.4
0.6
0.8
1.0
1.2
ξd
trimmed knot span (triangulated)
Figure 7: Construction of the integration domain for the IGA-IGA mortar mapper. Destination knot spans are projected onto the origin parameter space and subdivided along the origin knot lines, yielding conforming integration subdomains used for the evaluation of the mortar coupling integrals.
given interface displacement field uΓ , the acoustic solver computes the corresponding acoustic pressure field and the resulting interface load, FΓ = A(uΓ ),
(24)
where A(·) denotes the acoustic solution operator. Similarly, for a given interface load FΓ , the structural solver computes the corresponding structural response and the updated interface displacement, ũΓ = S(FΓ ),
(25)
where S(·) represents the structural solution operator. The coupled vibroacoustic problem can therefore be written as uΓ = S (A(uΓ )) ,
(26)
which corresponds to a fixed-point problem on the interface unknowns. Introducing the residual operator r(uΓ ) = S (A(uΓ )) − uΓ ,
(27)
r(uΓ ) = 0.
(28)
the coupled problem is equivalently written as
13
At the discrete level, the interface displacement, load, and residual fields are represented by the complex-valued vectors uΓ ∈ CNx , FΓ ∈ CNx , and rΓ ∈ CNx , respectively, where N x = No or N x = Nd depending on the selected interface representation used by the coupling algorithm. A standard Gauss–Seidel partitioned iteration corresponds to the fixed-point update u(k+1) = S A(u(k) (29) Γ Γ ) . where k denotes the coupling iteration index. Although this procedure preserves the modularity of the individual solvers, the fixed-point iterations may converge slowly or even diverge for strongly coupled vibroacoustic configurations, particularly in the vicinity of resonance frequencies where the interaction between the structural and acoustic fields becomes more pronounced. In such cases, convergence acceleration techniques become essential to stabilize and accelerate the iterative coupling procedure. Convergence accelerators generally operate on the residual of the fixed-point iteration. Let the interface residual at coupling iteration k be defined as r(k) = x̃(k+1) − x(k) ,
r(k) , x(k) ∈ CNx .
(30)
Here, x(k) denotes the current interface solution and x̃(k+1) the interface values obtained after one complete coupling step. The vector x is used as a generic notation for the accelerated interface quantity and may represent either interface displacements or interface loads, depending on the selected coupling formulation. The standard Aitken relaxation updates the interface solution according to [30] x(k+1) = x(k) + α(k) r(k) ,
(31)
where α(k) is a dynamically updated relaxation factor. In the classical real-valued formulation, the relaxation parameter is computed as T r(k−1) r(k) − r(k−1) (32) α(k) = −α(k−1) . r(k) − r(k−1) T r(k) − r(k−1) The IQN-ILS method follows a different strategy [29, 31]. Instead of applying a scalar relaxation to the interface residual, IQN-ILS constructs an approximation of the inverse interface Jacobian directly from previous coupling iterations. The method assumes that changes in the interface solution and the corresponding residual variations approximately satisfy ∂r ∆x ≈ J−1 ∆r, J= , (33) ∂x where J ∈ CNx ×Nx denotes the Jacobian of the discrete residual with respect to the interface solution vector. The differences between successive coupling iterations are defined as ∆ri = r(i+1) − r(i) ,
∆x̃i = x̃(i+1) − x̃(i) .
These vectors are assembled into the matrices h i V = ∆r1 ∆r2 · · · ∆rm ,
h W = ∆x̃1
∆x̃2
(34)
···
i ∆x̃m .
(35)
where V, W ∈ CNx ×m , ∆ri ∈ CNx , and ∆x̃i ∈ CNx . The parameter m denotes the number of residual and solution difference vectors retained from previous coupling iterations, i.e., the number of columns of the history matrices V and W. The current residual is then approximated as a linear combination of previously observed residual differences by solving the least-squares problem Vc ≈ −r(k) , (36) Finally, the interface correction is computed as x(k+1) = x(k) + Wc + r(k) . 14
(37)
Since most co-simulation infrastructures exchange interface quantities only as real-valued vectors, complex-valued interface fields arising in frequency-domain simulations are commonly transferred as separated real and imaginary components. A straightforward extension therefore consists of applying conventional real-valued convergence accelerators independently to both fields. However, this treatment neglects the intrinsic coupling between amplitude and phase information and may lead to deteriorated convergence behavior. To address this issue, the present work reconstructs the interface residual directly in the complex domain and performs the convergence acceleration using complex-valued interface quantities. For the Aitken method, two different complex-valued formulations are investigated. The first formulation computes the relaxation parameter directly from the complex-valued residual using a Hermitian inner product, while the second formulation evaluates the relaxation parameter from the component-wise modulus of the complex residual. In the following, these approaches are referred to as the complex Aitken and complex modulus Aitken methods, respectively. In contrast, the IQN-ILS method is extended by constructing the quasi-Newton approximation directly in the complex domain. For both Aitken and IQN-ILS methods, the interface residual is reconstructed as a single complex-valued quantity, (k) (k) r(k) c = rreal + i rimag ,
(k) N x (k) Nx r(k) c ∈ C , rreal , rimag ∈ R .
(38)
For the complex Aitken, the relaxation parameter is then computed using the Hermitian inner product ⟨a, b⟩ = aH b, where (·)H denotes the conjugate transpose, which yields the complex-valued Aitken update H r(k−1) r(k) − r(k−1) α(k) = −α(k−1) . r(k) − r(k−1) H r(k) − r(k−1)
(39)
Since the relaxation parameter is complex-valued, the interface correction (k) (k) ∆x(k) c = α rc .
(40)
can be interpreted geometrically as a simultaneous scaling and rotation of the complex residual vector in the complex plane. The magnitude of α(k) controls the amplitude of the interface correction, while its phase determines the rotation applied to the correction direction. In contrast to applying two independent real-valued convergence accelerators separately to the real and imaginary interface fields, the present formulation treats the interface quantities as a unified complex-valued residual field. Consequently, a single relaxation parameter is applied simultaneously to both components, thereby preserving the intrinsic coupling between amplitude and phase information. Although the relaxation parameter is computed from complex-valued residuals, the numerical experiments indicate that unrestricted phase rotations of the scalar relaxation coefficient may deteriorate the robustness of the fixedpoint iterations, particularly near resonance frequencies. Improved convergence behavior was obtained when the relaxation parameter remained predominantly aligned with the real axis. Therefore, the phase of the complex relaxation parameter is regularized by constraining it within a prescribed interval around the real axis. More specifically, the phase of α(k) is bounded around either 0 or π, depending on the sign of its real part. These observations also motivated the investigation of an alternative formulation in which the relaxation parameter remains purely real-valued while still being applied to the reconstructed complex-valued residual field. Figure 8 illustrates the admissible region of the complex relaxation parameter in the complex plane. The admissible phase interval is restricted to sectors of opening angle ∆ϕmax around the positive and negative real axes, thereby preventing excessively large rotations of the interface correction direction during the fixed-point iterations. This strategy preserves the coupled complex-valued treatment of the interface residual while preventing excessively large phase rotations of the scalar relaxation parameter. The resulting formulation can therefore be interpreted as a phase-regularized extension of the classical real-valued Aitken relaxation method. An alternative relaxation strategy is also investigated in which the Aitken relaxation parameter is computed from the component-wise magnitude of the complex residual vector. In this formulation, the complex residual is transformed into the real-valued residual (k) r(k) (41) abs = rc . 15
Im(α)
20◦ = ∆ϕmax
20◦ = ∆ϕmax α(k) Re(α)
A20 = α ∈ C : | arg(α)| ≤ 20 ∨ | arg(−α)| ≤ 20 Figure 8: Schematic representation of the phase regularization applied to the complex-valued Aitken relaxation parameter. The admissible phase variation is restricted to sectors of opening angle ∆ϕmax around the positive and negative real axes in order to avoid excessively large rotations of the interface correction direction.
where the absolute value is applied component-wise. The relaxation parameter is then computed using the classical real-valued Aitken expression, (k−1) T (k) rabs rabs − r(k−1) abs α(k) = −α(k−1) (k) . (k−1) T r(k) rabs − r(k−1) abs − rabs abs
(42)
The resulting scalar relaxation factor remains real-valued, but is applied to the full complex residual field, (k) (k) ∆x(k) c = α rc .
The interface correction is subsequently decomposed into its real and imaginary contributions, (k) (k) (k) , ∆x = ℑ ∆xc . ∆x(k) = ℜ ∆x c imag real
(43)
(44)
which are transferred back to the individual solvers as independent real-valued interface updates through the cosimulation communication routines. A similar extension is introduced for the IQN-ILS method. In this case, the matrices of residual and solution differences are assembled directly using complex-valued interface quantities, and the least-squares problem is solved in the complex domain, Vc ≈ −r(k) , V ∈ CNx ×m , r(k) ∈ CNx , c ∈ Cm . (45) In contrast to scalar relaxation methods, the complex-valued IQN-ILS formulation does not compute a single complex relaxation coefficient. Instead, the method constructs a multi-dimensional approximation of the inverse interface Jacobian directly in the complex domain. As a consequence, the interface correction is not restricted to a single scaling and rotation of the residual vector, but may involve a more general transformation of the complex-valued interface residual based on information gathered from previous coupling iterations. As a consequence, the coefficients defining the quasi-Newton update become complex-valued and naturally account for the phase relation between the structural and acoustic fields. Although the quasi-Newton operations are performed using complex-valued quantities, the final interface correction is again decomposed into real and imaginary components before being communicated to the individual solvers. 5. Numerical results This section presents three benchmark problems to assess the accuracy, robustness, and convergence behavior of the proposed partitioned vibroacoustic co-simulation framework. The examples are organized in a progressive 16
manner, increasing in complexity from one-way coupled configurations to fully bidirectional vibroacoustic interaction problems. Unless otherwise stated, the acoustic response is evaluated over frequency sweeps with a frequency step of ∆ f = 1 Hz and reported in terms of the magnitude of the complex acoustic pressure, |p|. All partitioned simulations employ a Gauss–Seidel communication pattern in which the structural and acoustic solvers are executed sequentially and the most recently available interface quantities are exchanged between both participants during the coupling procedure. Since the monolithic and partitioned formulations are implemented in different software environments, direct wall-clock time comparisons are not considered. The partitioned approach is instead evaluated in terms of solution accuracy, convergence behavior, and the number of coupling iterations. 5.1. One-way coupled rectangular plate The first numerical example considers a one-way coupled vibroacoustic benchmark to assess the influence of the interface mapping strategy on the predicted acoustic response. The structural displacements are transferred to the acoustic solver and imposed as boundary conditions for the Helmholtz problem, while acoustic feedback is neglected. Consequently, the observed differences between the monolithic and partitioned solutions can be primarily attributed to the transfer of interface quantities rather than to coupling-iteration errors. Since the coupling is unidirectional, no coupling iterations are required and, consequently, no convergence accelerator is employed in this example. The benchmark consists of a simply supported rectangular plate with dimensions L = 1 m and w = 0.5 m and, subjected to a unit harmonic surface load acting normal to the plate surface, as illustrated in Figure 9. Pressure Probe @ [0.5 m, 0.25 m, 1.0 m]
p = 1 Pa
w = 0.5 m
z
y x L=1m
Figure 9: Simply supported rectangular plate subjected to a uniform transverse surface load
The plate is modeled as a thin, homogeneous isotropic shell with thickness t = 5×10−3 m, density ρ = 7800 kg/m3 , Young’s modulus E = 200 × 109 Pa, and Poisson’s ratio ν = 0.3. The surrounding acoustic medium is assumed to be air with sound speed c = 340 m/s and density ρ f = 1.225 kg/m3 . The vibroacoustic response is evaluated over the frequency range 0–1000 Hz. Both undamped and damped structural configurations are considered. While undamped systems are not representative of practical engineering applications, they provide a particularly demanding benchmark for partitioned coupling. For the damped case, Rayleigh damping is introduced using α = 6.8 and β = 4.8 × 10−6 , corresponding to modal damping ratios of approximately 0.7–0.9% over the frequency range investigated. For the initial validation of the partitioned framework, interface quantities are transferred using the nearest-element mapping approach described in Section 4.1.2. The influence of alternative mapping strategies on the vibroacoustic response is investigated subsequently in Section 5.1.1. The structural interface is discretized using quadratic NURBS basis functions (p = q = 2) with 30 × 30 knot spans, while the acoustic interface employs a quadratic NURBS discretization with 18 × 9 knot spans. The acoustic discretization is selected according to [42]. The radiated acoustic field is evaluated in terms of the acoustic pressure magnitude at the observation point x = [0.5, 0.25, 1.0], shown in Figure 9. Figure 10 compares the acoustic pressure magnitude obtained using the monolithic reference solution and the proposed partitioned co-simulation framework for both the undamped and damped structural configurations. In both 17
cases, very good agreement is observed across the investigated frequency range, with the partitioned solution accurately reproducing both the resonance frequencies and the pressure amplitudes of the monolithic reference. Slightly larger discrepancies may occur in the vicinity of resonance peaks, where small errors introduced by the transfer of the interface quantities can be amplified by the increased sensitivity of the acoustic response. The damped configuration exhibits smoother resonance peaks and lower pressure amplitudes compared to the undamped case, while maintaining the same overall frequency response trend. 0.30
Monolithic Partitioned (Nearest Element Mapper)
Monolithic Partitioned (Nearest Element Mapper)
0.0200
0.25
Absolute pressure |p̂a | [Pa]
Absolute pressure |p̂a | [Pa]
0.0175 0.0150
0.20
0.0125
0.15
0.0100 0.0075
0.10
0.0050
0.05
0.0025 0.00
0.0000 0
200
400
600
800
1000
0
Frequency f [Hz]
200
400
600
800
1000
Frequency f [Hz]
(a) Without damping
(b) With damping
Figure 10: Comparison of the acoustic pressure magnitude obtained with the monolithic and partitioned one-way vibroacoustic coupling approaches: (a) undamped structural model and (b) damped structural model.
The undamped configuration is included primarily for completeness, since some level of damping is generally present in practical systems. Moreover, the dynamic stiffness matrix becomes singular at the exact eigenfrequencies and ill-conditioned in their vicinity, making the response highly sensitive to small discretization and mapping errors. Results close to the resonances should therefore be interpreted with particular care. 5.1.1. Influence of the interface mapper on the vibroacoustic response The previous results demonstrated that the partitioned framework is able to accurately reproduce the monolithic vibroacoustic response when using the nearest-element mapper. In the following, the sensitivity of the acoustic response to the selected interface transfer strategy is investigated. To emphasize the influence of the mapping procedure, the acoustic interface is intentionally discretized much more coarsely than the structural interface. The comparison is performed using a refined quadratic NURBS structural discretization with 30×30 knot spans and a coarse quadratic NURBS acoustic discretization with approximately 10×6 knot spans. To quantify the influence of the interface transfer procedure, the relative pressure amplitude error with respect to the monolithic reference solution is computed as erel ( f ) =
p̂a,part ( f ) − p̂a,mono ( f ) p̂a,mono ( f )
× 100%.
(46)
Figure 11 compares the relative acoustic pressure error at the probe location with respect to the monolithic reference solution, computed using Eq. (46), for the different interface transfer strategies and acoustic interface discretizations in the damped configuration. For the strongly non-conforming 10 × 6 acoustic discretization shown in Figure 11(a), the influence of the interface transfer strategy is clearly visible. All mapping approaches reproduce the overall acoustic response with good accuracy over most of the investigated frequency range, while larger deviations occur in the vicinity of structural resonances. In these regions, the vibroacoustic response becomes particularly sensitive to small perturbations in the transferred interface quantities. Among the considered approaches, the nearest-element mapper generally yields the 18
Relative pressure amplitude error [%]
Relative pressure amplitude error [%]
7
Nearest Element Mapper Nearest Neighbor Mapper Mortar Mapper
40
30
20
10
0
Nearest Element Mapper Nearest Neighbor Mapper Mortar Mapper
6 5 4 3 2 1 0
250
500
750
1000
1250
1500
1750
2000
200
Frequency f [Hz]
400
600
800
1000
Frequency f [Hz]
(a) Coarse acoustic discretization, 10 × 6
(b) Refined acoustic discretization, 18 × 9
Figure 11: Relative acoustic pressure error between the monolithic and partitioned vibroacoustic solutions for the nearest-element, nearest-neighbor, and mortar mapping strategies. For the strongly non-conforming 10 × 6 acoustic discretization in (a), the influence of the interface transfer strategy is clearly visible. Refining the acoustic interface to 18 × 9 knot spans in (b) significantly reduces the pressure errors and the differences between the considered mapping approaches.
smallest pressure errors, whereas the nearest-neighbor mapper exhibits larger error peaks. The mortar mapper shows a similar qualitative behavior, although larger deviations are observed for the intentionally coarse acoustic destination discretization. The influence of the interface transfer strategy is substantially reduced when the acoustic interface is refined to 18 × 9 knot spans, as shown in Figure 11(b). The relative pressure errors decrease significantly for all considered mapping approaches, with the nearest-element and mortar mappers remaining below 1% throughout the investigated frequency range. This reduction demonstrates that the larger discrepancies observed for the 10 × 6 case are strongly influenced by the severe mismatch between the interface discretizations. It should be emphasized that the coarse configuration is deliberately designed to amplify the influence of the interface transfer procedure and therefore represents a particularly challenging non-conforming coupling case. The reported errors should consequently not be interpreted as the accuracy expected for practical discretizations. Rather, the comparison demonstrates that the choice of interface mapper can noticeably influence the predicted acoustic response for strongly non-conforming interfaces, whereas these differences become considerably smaller as the acoustic interface discretization is refined. 5.2. Two-way coupled rectangular plate submerged in water The second example extends the previous configuration to a fully two-way coupled vibroacoustic problem. The same simply supported rectangular plate is subjected to a harmonic structural surface load and to the acoustic excitation generated by a monopole source located within the surrounding water domain (Figure 12). The structural and acoustic fields interact through the bidirectional exchange of interface quantities, which are transferred using the nearest-element mapping approach described in Section 4.1.2. In contrast to the one-way coupled configuration, the present problem introduces a strong vibroacoustic interaction and therefore constitutes a significantly more demanding partitioned solution procedure. Owing to the substantially higher density of water compared to air, the acoustic feedback acting on the structure is considerably stronger, making this benchmark particularly suitable for assessing the robustness of the proposed partitioned framework and the effectiveness of the complex convergence acceleration strategies introduced in Section 4.2. The structural model employs the same material properties as in Example 1. However, to obtain a representative submerged configuration, the plate thickness is increased to t = 0.01 m. Only the damped structural configuration is considered in this example. Structural damping is modeled using the same Rayleigh damping coefficients as in Example 1, namely α = 6.8 and β = 4.8 × 10−6 . For all strongly coupled simulations presented in this example, convergence is assessed using the absolute norm of the interface residual according to 19
3
Acoustic monopole @ [0.5 m, 0.25 m, 1.0 m], Qm = 1 ms
p = 1 Pa
w = 0.5 m
Pressure Probe @ [0.5 m, 0.25 m, 0.2 m]
z
y x L=1m
Figure 12: Two-way coupled vibroacoustic benchmark problem consisting of a simply supported rectangular plate subjected to a harmonic surface load and coupled with an acoustic monopole source located in the surrounding water domain. The structural and acoustic fields interact through the bidirectional exchange of interface quantities, while the sound pressure level (SPL) is evaluated at a probe position located above the plate surface.
∥r(k) ∥ ≤ 10−6 , (k)
(47)
where r denotes the interface residual at coupling iteration k. The selected tolerance was found to provide partitioned solutions in very close agreement with the corresponding monolithic reference while maintaining a reasonable computational cost. For the present strongly coupled configuration, the partitioned fixed-point iteration without convergence acceleration fails to converge. The use of a convergence accelerator is therefore required to obtain a converged partitioned solution. Figure 13 compares the acoustic response obtained with the monolithic reference solution and the strongly coupled partitioned formulation using the proposed complex-valued IQN-ILS, complex Aitken, and complex modulus Aitken convergence accelerators. The response is reported in terms of both the acoustic pressure magnitude and the sound pressure level (SPL). Although the pressure magnitude is used as the primary quantity throughout this work, the SPL representation is additionally included here because its logarithmic scale facilitates the visual comparison of the different solutions over the large dynamic range associated with the resonance peaks. The corresponding number of coupling iterations required for convergence at each excitation frequency is also reported. As shown in Figures 13a and 13b, the complex IQN-ILS accelerator remains in very close agreement with the monolithic reference solution throughout the investigated frequency range, demonstrating the capability of the proposed partitioned framework to accurately reproduce the strongly coupled vibroacoustic response. The SPL representation provides a complementary visualization of this agreement, particularly in frequency regions where the acoustic pressure varies over several orders of magnitude. This agreement is further quantified by the relative acoustic pressure error shown in Figure 14, computed with respect to the monolithic reference solution using Eq. (46). The error remains below 0.25% over most of the investigated frequency range, confirming the accuracy of the complex IQN-ILS coupling strategy. As shown in Figure 13c, the number of coupling iterations required by the complex IQN-ILS method increases gradually with frequency, reflecting the stronger acoustic–structure interaction at higher excitation frequencies. Nevertheless, fewer than ten iterations are required below approximately 400 Hz, and convergence is achieved well below the prescribed maximum of 50 iterations over the complete frequency sweep. In contrast, as shown also in Figures 13a and 13b, both Aitken-based accelerators exhibit significantly poorer convergence behavior. The complex Aitken produces solutions only up to approximately 427 Hz, while the complex modulus Aitken method does so only up to approximately 191 Hz. Beyond these frequencies, the fixed-point iterations become numerically unstable and eventually produce non-finite values, preventing the computation of a partitioned solution. Even within the frequency ranges where the methods remain operational, both Aitken variants generally require more coupling iterations than the complex IQN-ILS approach and frequently fail to satisfy the prescribed convergence criterion. In Figure 13a, the pressure range is limited to preserve the visibility of the monolithic and 20
Monolithic Complex IQN-ILS Complex Aitken Complex Modulus Aitken
14 Absolute pressure |p̂a | [Pa]
12
Monolithic Complex IQN-ILS Complex Aitken Complex Modulus Aitken
120 110 SPL [dB]
10 8 6
100 90 80
4 70 2 60 0 0
200
400 600 Frequency f [Hz]
800
1000
0
200
(a) Acoustic pressure magnitude.
400 600 Frequency f [Hz]
800
1000
(b) Sound pressure level.
55 Complex IQN-ILS Complex Aitken Complex Modulus Aitken
50 45
Iterations [-]
40 35 30 25 20 15 10 5 0
0
200
400 600 Frequency f [Hz]
800
1000
(c) Number of coupling iterations.
Figure 13: Comparison of the proposed complex-valued convergence accelerators for the submerged two-way coupled vibroacoustic benchmark: (a) acoustic pressure magnitude, (b) sound pressure level, and (c) number of coupling iterations required for convergence. The acoustic responses are compared with the monolithic reference solution, and the maximum number of coupling iterations is set to 50.
Relative pressure amplitude error [%]
Complex IQN-ILS
0.20
0.15
0.10
0.05
0.00 0
200
400 600 Frequency f [Hz]
800
1000
Figure 14: Relative acoustic pressure error between the monolithic and partitioned solutions obtained using the complex IQN-ILS convergence accelerator.
21
complex IQN-ILS solutions, since the very large pressure amplitudes obtained with the poorly converged Aitken solutions would otherwise dominate the linear scale. These larger deviations remain visible in the SPL representation in Figure 13b due to its logarithmic scale. The convergence behavior of the Aitken-based accelerators observed in the present study is consistent with previous findings in the partitioned FSI literature. For strongly coupled problems, particularly those exhibiting pronounced added-mass effects, scalar relaxation methods such as Aitken are known to exhibit reduced robustness and slower convergence than interface quasi-Newton methods [29, 43]. In contrast, IQN-ILS approximates the inverse interface Jacobian from previous coupling iterations, thereby providing a more robust and efficient coupling strategy for strongly coupled partitioned simulations. The behavior observed for the present submerged vibroacoustic problem is therefore in agreement with these findings from the partitioned FSI literature. These results indicate that the proposed complex IQN-ILS formulation provides a more robust and efficient coupling strategy than the investigated Aitken-based alternatives for the present benchmark problem. In particular, it remains in very close agreement with the monolithic solution while converging over the entire investigated frequency range. To further investigate the importance of the proposed complex-valued formulation, the results are compared against an alternative strategy in which the real and imaginary interface quantities are accelerated independently. Figure 15 compares the corresponding SPL with the monolithic reference solution and reports the associated number of coupling iterations. 250 225
50 45
200
40
175
35
Iterations [-]
SPL [dB]
55
Monolithic IQN-ILS separated fields Aitken separated fields
150 125
30 25 20 15
100
10
75
IQN-ILS separated fields Aitken separated fields
5
50 0
200
400 600 Frequency f [Hz]
800
0
1000
(a) Comparison of the monolithic solution with the strongly coupled partitioned solutions obtained when accelerating the real and imaginary interface quantities independently.
0
200
400 600 Frequency f [Hz]
800
1000
(b) Number of coupling iterations required by the separated-field acceleration strategies.
Figure 15: Performance of the separated-field acceleration strategies for the submerged two-way coupled vibroacoustic benchmark. The left figure compares the SPL response with the monolithic reference solution, while the right figure reports the corresponding number of coupling iterations required for convergence. The maximum number of coupling iterations is set to 50.
As shown in Figure 15a, the separated-field IQN-ILS formulation provides finite partitioned solutions over the entire investigated frequency range. However, the agreement with the monolithic reference deteriorates with increasing frequency, with substantial deviations occurring at higher frequencies. The separated-field Aitken formulation exhibits an even poorer behavior, with deviations from the monolithic reference developing at lower frequencies. For clarity, the SPL curves of both formulations are truncated once these deviations become sufficiently large to compromise the readability of the comparison. The convergence histories in Figure 15b provide further insight into this behavior. For the separated-field IQNILS formulation, the number of coupling iterations increases substantially with frequency and frequently reaches the prescribed maximum of 50 iterations. Thus, although finite solutions are obtained over the entire frequency range, the prescribed convergence criterion is not satisfied at many frequencies. For the separated-field Aitken formulation, the convergence behavior is even less robust. In addition to frequently reaching the maximum number of coupling iterations, finite interface values are obtained only up to approximately 570 Hz, beyond which the fixed-point iterations become numerically unstable and generate non-finite values. Compared with the genuinely complex-valued formulation, independently accelerating the real and imaginary 22
interface quantities leads to reduced convergence robustness and a noticeable deterioration of the predicted acoustic response. These results indicate that preserving the coupling between the amplitude and phase of the interface quantities during the convergence acceleration is important for accurately solving strongly coupled frequency-domain vibroacoustic problems. 5.3. Plane-wave excitation of a spherical shell coupled to interior and exterior acoustic domains The final example considers a thin spherical shell separating an enclosed acoustic fluid from an unbounded exterior acoustic domain. The shell is excited by a harmonic plane wave propagating through the exterior fluid and interacts dynamically with the acoustic fields on both sides of its midsurface. The incident wave excites the shell, whose vibration generates an outgoing scattered field in the exterior domain and an acoustic field within the enclosed cavity. This configuration is used to assess the accuracy and robustness of the proposed partitioned framework in the presence of strong bidirectional fluid–structure interaction and independently discretized fluid–structure interfaces. An analytical solution based on the spherical-harmonic formulation of Junger and Feit [20] is available for this problem and provides the exterior acoustic pressure used as the reference solution for the undamped configuration. The problem configuration is illustrated in Figure 16. The shell has midsurface radius a and thickness t, and is free in space. A harmonic plane wave of pressure amplitude A propagates in the positive z-direction. The exterior acoustic pressure is evaluated at the observation point xobs = (0, 0, 5) m, located on the positive z-axis.
xSPL = (0, 0, 5)m
y z x
O
+z Figure 16: Problem configuration of the spherical shell separating an enclosed acoustic fluid from an unbounded exterior fluid. The shell is excited by a harmonic plane wave propagating in the positive z-direction, and the exterior acoustic pressure is evaluated at xobs = (0, 0, 5) m.
The material, geometrical, and acoustic properties employed in this example are summarized in Table 1. For the numerical simulations, independent multipatch CAD representations are employed for the structural and acoustic domains, as shown in Figure 17. The structural shell is represented by two trimmed quadratic NURBS patches. Due to repeated knots in the underlying spline parameterization, the shell midsurface is only C 0 -continuous along the corresponding knot lines and patch interface. Consequently, the isogeometric Kirchhoff–Love shell formulation cannot be employed, as it requires a globally C 1 -continuous displacement field over the shell midsurface. The structural domain is therefore discretized using the isogeometric Reissner–Mindlin shell formulation of Benson et al. [33], in which the translational and rotational degrees of freedom are treated as independent unknowns. Since both the displacement and rotational fields require C 0 -continuity, their continuity across the patch interface is weakly enforced through a penalty formulation. The corresponding penalty energy is defined as Z 1 Πpen = αu ∥⟦u⟧∥2 + αθ ∥⟦θ⟧∥2 dΓ, (48) 2 Γint 23
Table 1: Material, geometrical, and acoustic properties of the spherical shell benchmark.
Property Exterior fluid density Exterior speed of sound Interior fluid density Interior speed of sound Structural density Young’s modulus Poisson’s ratio Shell midsurface radius Shell thickness Incident pressure amplitude
Patch 0 Patch 1
Symbol ρe ce ρi ci ρs E ν a t A
Value 1000 1500 1000 1500 7850 210 0.30 3.5 0.05 1
Unit kg/m3 m/s kg/m3 m/s kg/m3 GPa – m m Pa
Patch 0 Patch 1 Patch 2 Patch 3 Patch 4 Patch 5
4
4
2
2
0
z
0 −2
−2
−4
−4 −4
−2 x
0 2 4
−4
−2
2
0
−4
4
−2 x
0 2 4
y
−4
−2
2
0
4
y
Figure 17: Multipatch discretizations employed for the spherical shell benchmark. Left: structural shell represented by two trimmed NURBS patches. Right: acoustic boundary represented by six NURBS patches. The different patch layouts form a geometrically coincident but independently parameterized non-conforming fluid–structure interface.
where ⟦u⟧, ⟦θ⟧ ∈ R3 denote the displacement and rotational jumps across the patch interface, respectively, while αu and αθ are the corresponding displacement and rotational penalty parameters. In the present simulations, both penalty parameters are set to αu = αθ = 1014 . The resulting contribution is assembled into the structural stiffness matrix to weakly enforce continuity between adjacent patches. Further details on the formulation and implementation of the penalty coupling are given in [11]. The acoustic boundary is represented by six quartic NURBS patches. Since the acoustic formulation only requires C 0 -continuity of the pressure field, continuity across adjacent patches is enforced through the strong patch coupling formulation introduced in Section 2.2. Together with the two-patch structural representation, this results in a geometrically coincident but independently parameterized non-conforming fluid–structure interface. The proposed partitioned framework is evaluated for both undamped and Rayleigh-damped configurations using the same damping coefficients and convergence criteria introduced in Section 5.2. All partitioned simulations employ the complex IQN–ILS convergence accelerator together with the nearest-element mapping technique for the transfer of interface quantities between the structural and acoustic discretizations. The frequency response is computed over the range 20 ≤ f ≤ 100 Hz using a frequency increment of ∆ f = 0.2 Hz. The structural discretization is obtained through k-refinement of the initial quadratic NURBS representation. The 24
z
basis is first degree-elevated to cubic order and subsequently refined by knot insertion to obtain 60 knot spans in each parametric direction of every patch. The resulting structural model comprises 8978 control points, each possessing six degrees of freedom corresponding to three translational and three rotational components. The acoustic boundary is uniformly refined by knot insertion to obtain eight knot spans in each parametric direction of every patch. The final acoustic discretization comprises 728 control points associated with quartic NURBS basis functions, each carrying a single acoustic degree of freedom. 5.3.1. Results and discussion The resulting acoustic pressure magnitude, |p|, at the observation point xobs = (0, 0, 5) m is shown in Figure 18 for the undamped and Rayleigh-damped configurations. 14
Analytical [20] Partitioned co-simulation Monolithic in-house code
Monolithic in-house code Commercial software (COMSOL) Partitioned co-simulation
5 Absolute pressure |p̂a | [Pa]
Absolute pressure |p̂a | [Pa]
12
4
10 8
3
6
2
4
1
2 0
0 20
30
40
50 60 70 Frequency f [Hz]
80
90
100
20
(a) Undamped configuration.
30
40
50
60 70 Frequency f [Hz]
80
90
100
(b) Rayleigh-damped configuration.
Figure 18: Acoustic pressure magnitude at the observation point xobs = (0, 0, 5) m: (a) undamped configuration compared with the analytical solution of Junger and Feit [20] and the in-house monolithic solution, and (b) Rayleigh-damped configuration compared with the in-house monolithic and COMSOL Multiphysics solutions.
For the undamped configuration, Figure 18(a) compares the partitioned solution with both the analytical reference solution of Junger and Feit [20] and the monolithic solution obtained with the in-house research code. The monolithic model employs a six-patch representation of the spherical geometry for both the structural and acoustic domains, with each patch described by a NURBS surface with p = q = 4, where p and q denote the polynomial degrees in the two parametric directions. The structural and acoustic discretizations comprise 26 and 8 knot spans per patch, respectively. An isogeometric Kirchhoff–Love shell formulation is used for the structural domain, with continuity between the structural patches enforced using the multipatch coupling formulation proposed by Coox et al. [44]. In contrast, the partitioned model employs a two-patch structural representation together with an isogeometric Reissner– Mindlin shell formulation, for which continuity between the structural patches is weakly enforced using the penalty formulation described above, while retaining the six-patch representation for the acoustic domain. Very good agreement is obtained among the analytical, monolithic, and partitioned solutions over most of the investigated frequency range. The partitioned model accurately reproduces the overall acoustic pressure response and the resonance locations, while larger deviations occur primarily in the vicinity of the resonances. In the undamped case, the response is particularly sensitive to small perturbations in these frequency regions. Consequently, differences associated with the structural formulation and spatial discretization, as well as errors introduced by the interface mapping and iterative strong coupling procedure, can be significantly amplified near resonance. For the Rayleigh-damped configuration, Figure 18(b) compares the partitioned solution with the in-house monolithic solution and an independent reference obtained using COMSOL Multiphysics. In contrast to the proposed IBRA–IGA-BEM framework, the COMSOL model employs a monolithic FEM–FEM formulation with quadratic finite elements for both the structural and acoustic domains, resulting in a total of 109,531 degrees of freedom. Very good agreement is again obtained over the investigated frequency range. The partitioned solution closely follows the in-house monolithic reference, while somewhat larger differences with respect to the COMSOL solution occur 25
at higher frequencies and in the vicinity of resonance peaks. These differences should be interpreted in view of the different numerical formulations and spatial discretizations employed by the considered models. In addition to the pressure magnitude, the phase of the complex acoustic pressure provides a complementary measure of the agreement between the different solution strategies. Figure 19 compares the resulting phase angle obtained with the partitioned co-simulation formulation, the in-house monolithic code, and COMSOL Multiphysics. Monolithic in-house code Commercial software (COMSOL) Partitioned co-simulation
Pressure phase angle φpˆa [◦ ]
150 100 50 0 −50 −100 −150 20
30
40
50 60 70 Frequency f [Hz]
80
90
100
Figure 19: Phase angle of the complex acoustic pressure p̂a at the observation point xobs = (0, 0, 5) m for the Rayleigh-damped configuration. The partitioned solution is compared with the in-house monolithic and COMSOL Multiphysics reference solutions.
The phase comparison complements the pressure-magnitude results in Figure 18(b). The partitioned formulation combined with the proposed complex IQN-ILS accelerator closely reproduces the phase evolution predicted by the reference solutions, including the characteristic phase variations associated with the resonant response. This demonstrates that the proposed complex-valued coupling strategy reproduces not only the magnitude of the acoustic pressure but also its phase behavior, thereby accurately capturing the complex-valued acoustic response. The accuracy of the damped partitioned solution and the convergence behavior of the strongly coupled procedure are examined in more detail in Figure 20. Panel (a) shows the relative pressure amplitude error with respect to the in-house monolithic solution. Panel (b) reports the number of coupling iterations required for convergence for both the undamped and Rayleigh-damped configurations. 45
Damped Undamped
40 35
2.0
Iterations [-]
Relative pressure amplitude error [%]
Monolithic vs partitioned 2.5
1.5
30 25 20
1.0
15 10
0.5
5 0.0
20
30
40
50 60 70 Frequency f [Hz]
80
20
90
(a) Relative pressure amplitude error for the Rayleigh-damped case.
30
40
50 60 70 Frequency f [Hz]
80
90
100
(b) Number of coupling iterations.
Figure 20: Accuracy and convergence behavior of the partitioned solution using complex IQN-ILS for the spherical shell subjected to plane-wave excitation: (a) relative pressure amplitude error; and (b) number of coupling iterations for the undamped and Rayleigh-damped configurations.
As shown in Figure 20(a), the partitioned and monolithic solutions remain in close agreement in pressure ampli26
tude over most of the investigated frequency range. The three peaks in the relative pressure amplitude error correspond to the resonances, with the largest error occurring at the highest-frequency resonance. Since the monolithic and partitioned models differ in their structural representation and shell formulation, and the partitioned model additionally involves interface mapping and iterative strong coupling, the observed differences cannot be attributed exclusively to the coupling procedure. The concentration of error at the resonances is consistent with the increased sensitivity of the vibroacoustic response in these regions. The convergence histories in Figure 20(b) show a general increase in the number of coupling iterations with excitation frequency, together with pronounced peaks in the vicinity of the structural resonances. For the undamped structural problem, the dynamic stiffness matrix, K − ω2 M, becomes singular at the exact eigenfrequencies and increasingly ill-conditioned in their vicinity. Small perturbations in the transferred interface quantities can therefore produce comparatively large changes in the structural response, making the strongly coupled partitioned iterations more difficult to converge. The simultaneous increase in the pressure differences and coupling iterations toward the upper end of the frequency range further reflects the increased numerical sensitivity of the coupled problem in this region. Rayleigh damping regularizes the structural response in the vicinity of resonance and consequently reduces the pronounced peaks in the number of coupling iterations. Away from these frequency regions, the undamped and damped configurations exhibit comparable convergence behavior. Overall, the results demonstrate that the proposed partitioned framework remains robust over the investigated frequency range, while highlighting the increased numerical sensitivity of the strongly coupled problem near structural resonances and the beneficial effect of structural damping on the convergence behavior. 6. Conclusions This work presented a partitioned co-simulation strategy for frequency-domain vibroacoustic analysis combining isogeometric B-Rep structural models (IBRA) with an external IGA-BEM acoustic solver for unbounded domains. The proposed framework enables a fully modular and non-intrusive coupling procedure, while supporting fully CADintegrated workflows based on spline-native structural and acoustic interface representations. The numerical examples demonstrated very good agreement between the partitioned and monolithic vibroacoustic solutions for both one-way and two-way coupled problems. In particular, the proposed strategy accurately captured the acoustic pressure response and the resonance behavior of the investigated structures over the considered frequency range. The mapping studies highlighted the importance of the interface transfer operator for non-conforming spline-based interfaces. Among the investigated approaches, the nearest-element mapper provided the most accurate and reliable transfer behavior for strongly non-conforming discretizations and highly oscillatory interface fields. For strongly coupled vibroacoustic problems, the coupling iterations became highly sensitive near structural resonance frequencies. Since most existing co-simulation infrastructures exchange interface data as real-valued vectors, frequency-domain complex quantities are commonly transferred by separating their real and imaginary components. This naturally motivates applying standard convergence accelerators independently to both fields. However, the presented numerical results showed that this separated-field treatment leads to a deterioration of the convergence behavior and significantly reduced robustness, particularly for strongly coupled configurations and near resonance frequencies. To address this limitation, complex-valued extensions of the Aitken and IQN-ILS accelerators were introduced by reconstructing the complex interface residuals in the co-simulation tool and performing the convergence acceleration directly in the complex domain. The results demonstrated that consistently accounting for the coupled amplitudephase relation of the interface fields significantly improves the stability and convergence behavior of the partitioned solution procedure. Among the investigated strategies, the complex IQN-ILS formulation provided the best overall performance in terms of accuracy, robustness, and iteration count. Overall, the proposed partitioned IGA-BEM framework provides a flexible and accurate alternative to monolithic vibroacoustic formulations while preserving solver modularity and supporting fully CAD-integrated analysis of coupled structural-acoustic systems with non-conforming spline-based interfaces.
27
Appendix A. Implementation of the remote-controlled acoustic participant The remote-controlled execution strategy adopted for the acoustic participant is implemented through a Python interface layer that registers a set of callback functions in CoSimIO. This layer acts as a wrapper between the cosimulation infrastructure and the MATLAB-based acoustic BEM solver, handling data exchange, solver execution, and communication with the structural participant. Listing 1 summarizes the essential structure of the acoustic participant used in the present work. 1
# acoustic solver initializes ...
2 3
CoSimIO.Connect(...)
# establish connection
4 5 6 7 8
# define functions to be registered def ImportData(identifier): if identifier == "disp_real": CoSimIO.ImportData(Re(u_Gamma))
9 10 11
if identifier == "disp_imag": CoSimIO.ImportData(Im(u_Gamma))
12 13 14 15 16
def SolveSolutionStep(): # reconstruct complex interface displacement u_Gamma = Re(u_Gamma) + i * Im(u_Gamma)
17 18 19
# solve acoustic BEM problem f_Gamma = BEM(u_Gamma, omega)
20 21 22
# split complex acoustic load Re(f_Gamma), Im(f_Gamma)
23 24 25 26 27
def ExportData(identifier): if identifier == "load_real": CoSimIO.ExportData(Re(f_Gamma))
28 29 30
if identifier == "load_imaginary": CoSimIO.ExportData(Im(f_Gamma))
31 32 33 34 35 36
def ExportMesh(): # construct acoustic interface mesh from integration points X_Gamma = ComputeIntegrationPoints() CoSimIO.ExportMesh(X_Gamma)
37 38 39 40 41
def FinalizeSolutionStep(): # acoustic postprocessing SPL = Postprocess(f_Gamma, omega)
42 43 44 45 46 47 48 49
# after defining the functions they are registered in CoSimIO CoSimIO.Register(ImportData) CoSimIO.Register(SolveSolutionStep) CoSimIO.Register(ExportData) CoSimIO.Register(ExportMesh) CoSimIO.Register(FinalizeSolutionStep)
28
50 51 52
# after all functions are registered, execution control is handed over CoSimIO.Run(...)
53 54
CoSimIO.Disconnect(...)
# stop connection Listing 1: Remote-controlled acoustic BEM participant.
Declarations Conflict of interest. The authors have no competing interests to declare that are relevant to the content of this article. Data availability Data will be made available on request. Acknowledgments The authors gratefully acknowledge the Design for IGA-type discretization workflows (GECKO) project. The Design for IGA-type discretization workflows has received funding from the European Union’s Horizon Europe research and Innovation programme under grant agreement No. 101073106, Call: HORIZON-MSCA-2021-DN-01. Views and opinions expressed are however those of the authors and do not necessarily reflect those of the European Union. The European Union cannot be held responsible for them. References [1] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th Edition, Butterworth-Heinemann, Oxford, 2013. [2] K.-J. Bathe, Finite Element Procedures, Klaus-Jürgen Bathe, Watertown, MA, 2006. [3] J. Fish, T. Belytschko, A first course in finite elements, Vol. 1, Wiley New York, 2007. [4] D. S. Burnett, A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion, The Journal of the Acoustical Society of America 96 (5) (1994) 2798–2816. [5] R. Astley, J. Hamilton, Numerical studies of conjugated infinite elements for acoustical radiation, Journal of Computational Acoustics 8 (01) (2000) 1–24. [6] J.-P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of computational physics 114 (2) (1994) 185–200. [7] C. A. Brebbia, The boundary element method for engineers, Pentech Press (1978). [8] P. Banerjee, P. Banerjee, R. Butterfield, Boundary element methods in engineering science, McGraw-Hill Book Company (UK), 1981. [9] T. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39) (2005) 4135–4195. doi:https://doi.org/10.1016/j.cma.2004.10.008. [10] M. Breitenberger, A. Apostolatos, P. Bucher, R. Wüchner, K.-U. Bletzinger, Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures, Computational Methods in Applied Mechanics and Engineering 284 (2015) 401–457. [11] T. Teschemacher, A. M. Bauer, T. Oberbichler, M. Breitenberger, R. Rossi, R. Wüchner, K.-U. Bletzinger, Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis, Advanced Modeling and Simulation in Engineering Sciences 5 (2018) 19. [12] T. Teschemacher, A. M. Bauer, R. Aristio, M. Meßmer, R. Wüchner, K.-U. Bletzinger, Concepts of data collection for the CAD-integrated isogeometric analysis, Engineering with Computers 38 (6) (2022) 5675–5693. [13] R. Simpson, M. Scott, M. Taus, D. Thomas, H. Lian, Acoustic isogeometric boundary element analysis, Computer Methods in Applied Mechanics and Engineering 269 (2014) 265–290. doi:https://doi.org/10.1016/j.cma.2013.10.026. [14] L. Coox, O. Atak, D. Vandepitte, W. Desmet, An isogeometric indirect boundary element method for solving acoustic problems in openboundary domains, Computer Methods in Applied Mechanics and Engineering 316 (2017) 186–208, special Issue on Isogeometric Analysis: Progress and Challenges. doi:https://doi.org/10.1016/j.cma.2016.05.039. [15] G. C. Diwan, M. S. Mohamed, Pollution studies for high order isogeometric analysis and finite element for acoustic problems, Computer Methods in Applied Mechanics and Engineering 350 (2019) 701–718. doi:https://doi.org/10.1016/j.cma.2019.03.031. [16] Z. Liu, M. Majeed, F. Cirak, R. N. Simpson, Isogeometric fem-bem coupled structural-acoustic analysis of shells using subdivision surfaces, International Journal for Numerical Methods in Engineering 113 (9) (2018) 1507–1530. [17] Y. Wu, C. Dong, H. Yang, A 3d isogeometric FE-IBE coupling method for acoustic-structural interaction problems with complex coupling models, Ocean Engineering 218 (2020) 108183.
29
[18] Y. Wu, C. Dong, H. Yang, F. Sun, Isogeometric symmetric FE-BE coupling method for acoustic-structural interaction, Applied Mathematics and Computation 393 (2021) 125758. [19] G. C. Everstine, F. M. Henderson, Coupled finite element/boundary element approach for fluid–structure interaction, The Journal of the Acoustical Society of America 87 (5) (1990) 1938–1947. doi:10.1121/1.399186. [20] M. C. Junger, D. Feit, Sound, Structures, and Their Interaction, MIT Press, Cambridge, MA, 1986. [21] S. Sicklinger, C. Lerch, R. Wüchner, K.-U. Bletzinger, Fully coupled co-simulation of a wind turbine emergency brake maneuver, Journal of Wind Engineering and Industrial Aerodynamics 144 (2015) 134–145, selected papers from the 6th International Symposium on Computational Wind Engineering CWE 2014. doi:https://doi.org/10.1016/j.jweia.2015.03.021. URL https://www.sciencedirect.com/science/article/pii/S0167610515000811 [22] S. Sicklinger, V. Belsky, B. Engelmann, H. Elmqvist, H. Olsson, R. Wüchner, K.-U. Bletzinger, Interface jacobian-based co-simulation, International Journal for Numerical Methods in Engineering 98 (6) (2014) 418–444. arXiv:https://onlinelibrary.wiley.com/doi/ pdf/10.1002/nme.4637, doi:https://doi.org/10.1002/nme.4637. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.4637 [23] P. L. Bucher, Cosimulation and mapping for large scale engineering applications, Ph.D. thesis, Technische Universität München (2024). URL https://mediatum.ub.tum.de/1714023 [24] L. Rodrı́guez-Tembleque, J. A. González, A. Cerrato, Partitioned solution strategies for coupled bem–fem acoustic fluid–structure interaction problems, Computers & Structures 152 (2015) 45–58. doi:https://doi.org/10.1016/j.compstruc.2015.02.018. URL https://www.sciencedirect.com/science/article/pii/S0045794915000565 [25] G. Bunting, S. T. Miller, Partitioned coupling for structural acoustics, Journal of Vibration and Acoustics 142 (1) (2019) 011012. arXiv:https://asmedigitalcollection.asme.org/vibrationacoustics/article-pdf/142/1/011012/6448035/ vib_142_1_011012.pdf, doi:10.1115/1.4045215. URL https://doi.org/10.1115/1.4045215 [26] J. Kersschot, H. Denayer, W. De Roeck, W. Desmet, Simulation of strong vibro-acoustic coupling effects in ducts using a partitioned approach in the time domain, in: Proceedings of the International Conference on Noise and Vibration Engineering (ISMA 2020), Leuven, Belgium, 2020, pp. 285–294. [27] G. Chourdakis, K. Davis, B. Rodenberg, M. Schulte, F. Simonis, B. Uekermann, G. Abrams, H. Bungartz, L. Cheung Yau, I. Desai, K. Eder, R. Hertrich, F. Lindner, A. Rusch, D. Sashko, D. Schneider, A. Totounferoush, D. Volland, P. Vollmer, O. Koseomur, preCICE v2: A sustainable and user-friendly coupling library [version 2; peer review: 2 approved], Open Research Europe 2 (51) (2022). doi:10.12688/ openreseurope.14445.2. URL https://doi.org/10.12688/openreseurope.14445.2 [28] U. Küttler, W. A. Wall, Fixed-point fluid–structure interaction solvers with dynamic relaxation, Computational Mechanics 43 (1) (2008) 61–72. [29] J. Degroote, R. Haelterman, S. Annerel, P. Bruggeman, J. Vierendeels, Performance of partitioned procedures in fluid–structure interaction, Computers & Structures 88 (7) (2010) 446–457. doi:https://doi.org/10.1016/j.compstruc.2009.12.006. URL https://www.sciencedirect.com/science/article/pii/S0045794909003022 [30] B. M. Irons, R. C. Tuck, A version of the aitken accelerator for computer iteration, International Journal for Numerical Methods in Engineering 1 (3) (1969) 275–277. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620010306, doi:https://doi.org/ 10.1002/nme.1620010306. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620010306 [31] J. Degroote, K.-J. Bathe, J. Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid–structure interaction, Computers & Structures 87 (11) (2009) 793–801, fifth MIT Conference on Computational Fluid and Solid Mechanics. doi:https://doi.org/10.1016/j.compstruc.2008.11.013. URL https://www.sciencedirect.com/science/article/pii/S0045794908002605 [32] J. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner, Isogeometric shell analysis with kirchhoff-love elements, Computer Methods in Applied Mechanics and Engineering 198 (49–52) (2009) 3902–3914. doi:10.1016/j.cma.2009.08.013. [33] D. Benson, Y. Bazilevs, M. Hsu, T. Hughes, Isogeometric shell analysis: The reissner–mindlin shell, Computer Methods in Applied Mechanics and Engineering 199 (5) (2010) 276–289, computational Geometry and Analysis. doi:https://doi.org/10.1016/j.cma.2009.05. 011. URL https://www.sciencedirect.com/science/article/pii/S0045782509001820 [34] L. Coox, F. Greco, O. Atak, D. Vandepitte, W. Desmet, A robust patch coupling method for nurbs-based isogeometric analysis of nonconforming multipatch surfaces, Computer Methods in Applied Mechanics and Engineering 316 (2017) 235–260, special Issue on Isogeometric Analysis: Progress and Challenges. doi:https://doi.org/10.1016/j.cma.2016.06.022. URL https://www.sciencedirect.com/science/article/pii/S0045782516306120 [35] P. Dadvand, R. Rossi, E. Oñate, An object-oriented environment for developing finite element codes for multi-disciplinary applications, Archives of Computational Methods in Engineering 17 (3) (2010) 253–297. doi:https://doi.org/10.1007/s11831-010-9045-2. [36] P. Dadvand, R. Rossi, M. Gil, X. Martorell, J. Cotela, E. Juanpere, S. Idelsohn, E. Oñate, Migration of a generic multi-physics framework to hpc environments, Computers & Fluids 80 (2013) 301–309. doi:https://doi.org/10.1016/j.compfluid.2012.02.004. [37] H. G. Matthies, J. Steindorf, Partitioned strong coupling algorithms for fluid–structure interaction, Computers & Structures 81 (2003) 805– 812. doi:10.1016/S0045-7949(02)00409-1. [38] M. Mehl, B. Uekermann, H. Bijl, F. Blom, H.-J. Bungartz, A. H. van Zuijlen, Parallel coupling numerics for partitioned fluid–structure interaction simulations, Computers & Mathematics with Applications 71 (4) (2016) 869–891. doi:10.1016/j.camwa.2015.12.025. [39] J. L. Bentley, Multidimensional binary search trees used for associative searching, Communications of the ACM 18 (9) (1975) 509–517. doi:10.1145/361002.361007. [40] D. Meagher, Geometric modeling using octree encoding, Computer Graphics and Image Processing 19 (2) (1982) 129–147. doi:10.1016/ 0146-664X(82)90104-6.
30
[41] A. Apostolatos, A. Emiroğlu, S. Shayegan, et al., An isogeometric b-rep mortar-based mapping method for non-matching grids in fluidstructure interaction, Advanced Modeling and Simulation in Engineering Sciences 8 (9) (2021). doi:10.1186/s40323-021-00190-9. URL https://doi.org/10.1186/s40323-021-00190-9 [42] S. Marburg, Six boundary elements per wavelength: Is that enough?, Journal of Computational Acoustics 10 (01) (2002) 25–51. doi: 10.1142/S0218396X02001401. [43] A. Delaissé, J. Degroote, Quasi-newton methods for partitioned simulation of fluid-structure interaction, Archives of Computational Methods in Engineering 30 (2023) 2515–2558. [44] L. Coox, F. Maurin, F. Greco, E. Deckers, D. Vandepitte, W. Desmet, A flexible approach for coupling nurbs patches in rotationless isogeometric analysis of kirchhoff–love shells, Computer Methods in Applied Mechanics and Engineering 325 (2017) 505–531. doi:https://doi.org/10.1016/j.cma.2017.07.022. URL https://www.sciencedirect.com/science/article/pii/S0045782517305674
31