A novel smoothed particle hydrodynamics and finite element coupling scheme for fluid-structure interaction: the sliding boundary particle approach

A novel numerical formulation for solving fluid-structure interaction (FSI) problems is proposed where the fluid field is spatially discretized using smoothed particle hydrodynamics (SPH) and the structural field using the finite element method (FEM). As compared to fully mesh- or grid-based FSI frameworks, due to the Lagrangian nature of SPH this framework can be easily extended to account for more complex fluids consisting of multiple phases and dynamic phase transitions. Moreover, this approach facilitates the handling of large deformations of the fluid domain respectively the fluid-structure interface without additional methodological and computational efforts. In particular, to achieve an accurate representation of interaction forces between fluid particles and structural elements also for strongly curved interface geometries, the novel sliding boundary particle approach is proposed to ensure full support of SPH particles close to the interface. The coupling of the fluid and the structural field is based on a Dirichlet-Neumann partitioned approach, where the fluid field is the Dirichlet partition with prescribed interface displacements and the structural field is the Neumann partition subject to interface forces. To overcome instabilities inherent to weakly coupled schemes an iterative fixed-point coupling scheme is employed. Several numerical examples in form of well-known benchmark tests are considered to validate the accuracy, stability, and robustness of the proposed formulation. Finally, the filling process of a highly flexible thin-walled balloon-like container is studied, representing a model problem close to potential application scenarios of the proposed scheme in the field of biomechanics.


Introduction
In many applications in science and engineering fluid-structure interaction (FSI) phenomena play an essential role in modeling and simulation, in particular, in some areas of biomechancis, e.g., digestion of food in the human stomach [1,2], referring to the authors target application. Besides the challenge to deal with large deformations of both fluid and structural domain, accurate modeling of fluid flow in biomechanics is even more demanding in the case of complex fluids including, e.g., multiple fluid phases and dynamic phase transitions (e.g. due to chemical reactions). Most current FSI frameworks utilize mesh-or gridbased methods, e.g., the finite element method (FEM), finite difference method (FDM), or finite volume method (FVM), which often require additional methodological and computational effort to capture the aforementioned phenomena. A promising approach to model complex fluids, e.g., the content of gastric lumen in the human stomach [1,2], is the method of smoothed particle hydrodynamics (SPH). SPH is a mesh-free discretization scheme that was originally and independently of one another introduced by Gingold and Monaghan [3] and Lucy [4] in 1977. While initially intended to study astrophysical problems, SPH gained increasing importance in other fields of computational fluid dynamics (CFD) since then. Due to its Lagrangian nature, SPH is very well suited for flow problems involving multiple phases, dynamic phase transitions, as well as complex interface topologies. Especially for many fluid-structure interaction scenarios in biomechanics it would therefore be desirable to discretize the fluid field with SPH whereas the solid field is often easier to handle with finite elements. To this end, a robust and efficient algorithm coupling SPH and FEM for the simulation of fluid-structure interactions is required.
On these grounds, this contribution proposes a novel numerical formulation for solving FSI problems where the fluid field is modeled using SPH and the structural field using FEM. Both sub-fields are coupled following a Dirichlet-Neumann partitioned approach. The fluid field is the Dirichlet partition with prescribed interface displacements or interface velocities, respectively, and the structural field is the Neumann partition subject to interface forces. That means, the interface forces are evaluated by the fluid solver utilizing the current interface displacements and interface velocities that are directly extracted from the structural field. Afterwards, the interface forces are applied to the structural solver enforcing conservation of linear momentum. An iterative fixed-point coupling scheme [5] is employed to satisfy dynamic equilibrium at the fluid-structure interface with respect to a predefinded convergence criterion. This so-called strong coupling of both sub-fields is crucial to overcome instabilities, e.g., due to the artificial added mass effect, that are known to occur for weakly coupled schemes in FSI [6,7].
One focus of this work lies on the crucial aspect of the treatment of deformable and strongly curved boundaries of the SPH domain as especially required for many FSI applications. In the literature several different formulations for modeling (rigid) boundaries in SPH are proposed. Among them are penalty-like repulsive force formulations [8,9,10], ghost particle formulations [11], boundary particle methods based on fixed layers of particles resembling rigid walls [12,13], or semi-analytical methods considering nonvanishing surface integrals due to missing kernel support [14,15,16]. For an overview on the advantages and disadvantages of the aforementioned methods the interested reader is refered to the literature, e.g., in [17,18]. In principle, all those methods modeling rigid boundaries in SPH naturally have the potential to serve as a basis also for the treatment of flexible structural boundaries in the context of FSI problems [19,20,21,22,23]. However, FSI applications, especially in biomechanics, are characterized by large deformations at the fluid-structure interface including strong curvature and large stretch. This requires a special treatment of boundaries in order to prevent loss of accuracy at the fluid-structure interface. To the best of the authors' knowledge, the existing methods are either missing the required accuracy, computationally expensive, or not capable of modeling deforming interfaces undergoing strong curvature and large stretch. To address this shortcoming of existing approaches, the novel sliding boundary particle approach is proposed. It is based on a transient set of virtual boundary particles regulary arranged around the current projection point of a fluid particle onto the fluid domain boundary. Moreover, a generalized formulation for the extrapolation of field variables from fluid to virtual boundary particles is proposed, which is inspired by the procedure of [13].
The present publication is organized as follows: To begin with, the governing equations for FSI problems are briefly introduced in Section 2, followed by a detailed presentation of the numerical methods and the computational framework being utilized with a focus on the evaluation of the interface forces and the coupling scheme, cf. Section 3. Finally, numerical results obtained with the proposed novel numerical formulation for solving FSI problems are shown in Section 4. For validation purposes, well-known CFD respectively FSI benchmark tests are studied confirming the accuracy and robustness of the proposed formulation. This is followed by an application-motivated academic example examining the filling process of a highly flexible thin-walled container.

Governing equations
At all times t ∈ [0, T ] the domain Ω of a fluid-structure interaction problem consists of a non-overlapping fluid domain Ω f and a structural domain Ω s that share a common interface Γ f s , with Ω = Ω f ∪ Ω s and Ω f ∩ Ω s = Γ f s , refer to Figure 1. This leads to the so-called geometric coupling condition that restricts the fluid and structural domains to perfectly match without any holes or gaps at the fluid-structure interface Γ f s . In the following, the (standard) governing equations of the fluid and structural field as well as the respective coupling condition for FSI are briefly given.
PSfrag replacements Figure 1: Domain Ω of a fluid-structure interaction problem consisting of two disjunct sub-domains for the fluid field Ω f and the structural field Ω s with shared common interface Γ f s .

Remark 1.
In the equations (1)-(2) governing the fluid field and (5) governing the structural field all time derivatives follow the motion of material points, i.e., the material derivative reads d(·) dt = ∂(·) ∂t + u · ∇(·). Furthermore, ∇(·) denotes within the setting of nonlinear continuum mechanics derivatives with respect to spatial coordinates while ∇ 0 (·) denotes derivatives with respect to material coordinates.

Fluid field
The fluid field is governed by the instationary Navier-Stokes equations in the domain Ω f in convective form consisting of the mass continuity equation and the momentum equation with viscous force f ν and body force b f each per unit mass. For a Newtonian fluid the viscous force is f ν = ν f ∇ 2 u f with kinematic viscosity ν f . The mass continuity equation (1) and the momentum equation (2) represent a system of four equations with the five unknowns, velocity u f , density ρ f , and pressure p f . The system of equations is closed with an equation of state p f = p f ρ f relating fluid density ρ f and pressure p f , cf. Section 3.1.5. The Navier-Stokes equations (1) and (2) are subject to the following initial conditions with initial density ρ f 0 and initial velocity u f 0 . In addition, Dirichlet and Neumann boundary conditions are applied on the fluid boundary with prescribed boundary velocityû f and boundary

Structural field
Considering the regime of finite deformations, the structural field is governed by the balance of linear momentum in the following local material form 3 with the material forms of density ρ s 0 and body force b s 0 , and the structural displacement d s as primary unknowns. The deformation of the structure is described by the deformation gradient F = ∇ 0 d s defining the Green-Lagrange strains E = 1 2 F T F − I . For simplicity, and as applicable and most often used in biomechanical problems, the second Piola-Kirchhoff stresses S are chosen to follow from a constitutive relation of the form S = ∂Ψ/∂E based on a hyperelastic strain energy function Ψ = Ψ(E). The partial differential equation (5) is subject to initial conditions for the structural displacement and velocity On the structural boundary Γ s = ∂Ω s \ Γ f s , Dirichlet and Neumann boundary conditions are prescribed with prescribed boundary displacementd s , boundary tractiont s 0 , and outward pointing unit normal vector N on Γ s in material description, where Γ s = Γ s D ∪ Γ s N and Γ s D ∩ Γ s N = ∅.

Coupling conditions
A geometric coupling condition results from restricting both the fluid and structural domain to match at the fluid-structure interface Γ f s as already described in the beginning of this section. In addition, the socalled kinematic coupling condition (or no-slip boundary condition) enforces a continuous fluid and structural velocity at the interface Γ f s . Consequently, these two conditions can be expressed as with the current position r f respectively r s of the fluid and structural field. Finally, the dynamic coupling condition ensures equilibrium of fluid and structural traction across the interface Γ f s

Numerical methods and computational framework
The purpose of this section is to present the methods for discretization and numerical solution of the fluid-structure interaction problem as described in Section 2. The discretization of the fluid field is based on smoothed particle hydrodynamics while the discretization of the structural field is based on the finite element method, as illustrated in Figure 2 (left). While in Sections 3.1 and 3.2 the basics of these two methods are recapitulated, the focus of this publication is set on the specific evaluation of interaction forces, cf. Section 3.3, introducing the sliding boundary particle approach, and the employed coupling algorithm, cf. Section 3.4, in terms of underlying methods. The presented computational framework is implemented in the in-house parallel multiphysics research code BACI (Bavarian Advanced Computational Initiative) [24].

Discretization of fluid field via smoothed particle hydrodynamics
The fluid field governed by the instationary Navier-Stokes equations (1) and (2) is solved using smoothed particle hydrodynamics following a weakly compressible approach [9,17,25]. For modeling fluid flow using SPH, several different formulations each with its own characteristics and benefits can be derived as reflected by the vast amount of literature. The aim of this section is to give a brief introduction into the basics of SPH and an overview of the formulation applied throughout this work. Note that the contribution resulting from the coupling condition of the fluid and structural field at the interface Γ f s , refer to Section 2.3, is omitted in this section and described in detail in Section 3.3. For ease of notation, in the following the index (·) f denoting fluid quantities, as used in Section 2, is dropped.

Approximation of field quantities via smoothing kernel
The fundamental concept of SPH is based on the approximation of a field quantity f via a smoothing operation and on the discretization of the domain Ω with discretization points, so-called particles. To begin with, a field quantity f on a domain Ω can be expressed exactly in integral form as making use of the Dirac delta function δ(r). Replacing the latter by a so-called smoothing kernel W (r, h), that fulfills certain required properties, cf. Remark 2 and [17], leads to an approximation of the field quantity f in smoothed integral form while committing a smoothing error.
Remark 2. The smoothing kernel W (r, h) is a monotonically decreasing, smooth function that depends on a distance r and a smoothing length h. The smoothing length h together with a scaling factor κ define the support radius of the smoothing kernel r c = κh. Compact support, i.e., W (r, h) = 0 for r > r c , as well as positivity, i.e., W (r, h) ≥ 0 for r ≤ r c , are typical properties of standard smoothing kernels W (r, h). In addition, the normalization property requires that Ω W (|r − r ′ |, h) dr ′ = 1. The Dirac delta function property lim h→0 W (r, h) = δ(r) ensures an exact representation of a field quantity f in the limit h → 0.
In a next step, the computational domain Ω is filled with discretization points or so-called particles j, each occupying a volume V j . Thus, the smoothed integral form of quantity f reduces in discretized form to a summation of contributions from all particles j in the domain Ω, cf. Remark 3, adding a discretization error [26]. A straightforward approach in SPH to determine the gradient of a quantity f follows directly by differentiation of equation (12) resulting in Note that this (simple) variant for an approximation of the gradient shows some particular disadvantages, hence, more advanced approximations for gradients are given in the literature [9] and also applied in this work [13,27], cf. Section 3.1.4.

Remark 3.
In general, contributions from all particles in the domain Ω are considered in the SPH approximation of a field quantity f , cf. equation (12). However, note that in practice due to the compact support of the smoothing kernel W only neighboring particles within the support radius r c need to be considered. This property is very beneficial as it reduces the computational effort of the method.
Applying the concept of SPH reduces the partial differential equations (1) and (2) to ordinary differential equations that are solved, i.e., evaluated and integrated in time, for all particles in the domain Ω (cf. Sections 3.1.4 and 3.1.7). The transient positions of particles are advected with the fluid velocity resembling the Lagrangian nature of the method. As a result, all fluid quantities are evaluated at and associated with particle positions, meaning each particle carries its corresponding fluid quantities.
Finally, in a post-processing step the continuous field quantity f is recovered from the discrete fluid quantities carried by each particle in the domain based on approximation (12) and the commonly known Shepard filterf Note that the denominator typically takes on values close to one inside the fluid domain and is mainly relevant for boundary regions with reduced support due to a lack of neighboring particles.
Remark 4. In the following, a quantity f evaluated for particle i at position r i is written as f i = f (r i ).
In addition, the short notation W ij = W (r ij , h) denotes the smoothing kernel W evaluated for particle i at position r i with neighboring particle j at position r j , where r ij = |r ij | = |r i − r j | is the absolute distance between particles i and j. Similarly, the derivative of the smoothing kernel W with respect to the absolute distance r ij is denoted by ∂W /∂r ij = ∂W (r ij , h)/∂r ij .
Remark 5. Herein, the smoothing of fluid quantities is carried out using a quintic spline smoothing kernel W (r, h) as defined in [12] with smoothing length h and compact support of the smoothing kernel with support radius r c = κh and scaling factor κ = 3.

Initial particle spacing
Within this contribution, the fluid domain is initially filled with particles located on a regular grid with particle spacing ∆x, thus in a d-dimensional space each particle initially occupies an effective volume of (∆x) d . The mass of a particle i is then set using the reference density according to m i = ρ 0 (∆x) d and remains constant throughout the simulation. In general, the initial particle spacing ∆x can be freely chosen, however, within this work the initial particle spacing ∆x is set equal to the smoothing length h = r c /κ .

Density summation
The density of a particle i is determined via summation of the respective smoothing kernel contributions of all neighboring particles j within the support radius r c This approach is typically denoted as density summation and results in an exact conservation of mass in the fluid domain, which can be shown in a straightforward manner considering the commonly applied normalization of the smoothing kernel to unity. It shall be noted that the density field may alternatively be obtained by discretization and integration of the mass continuity equation (1) [17].

Momentum equation
The momentum equation (2) is discretized following [13,27] including a transport velocity formulation to suppress the problem of tensile instability. It will be briefly recapitulated in the following. The transport velocity formulation relies on a constant background pressure p b that is applied to all particles and results in a contribution to the particle accelerations for in general disordered particle distributions. However, these additional acceleration contributions vanish for particle distributions fulfilling the partition of unity, thus fostering these desirable configurations. For the sake of brevity, the definition of the modified advection velocity and the additional terms in the momentum equation from the aforementioned transport velocity formulation are not discussed in the following and the reader is kindly referred to the original publication [27]. Altogether, the acceleration a i = du i /dt of a particle i results from summation of all acceleration contributions due to interaction with neighboring particles j and a body force as with volume V i = m i /ρ i of particle i, unit vector e ij = r i − r j /|r i − r j | = r ij /r ij pointing from particle j to particle i, relative velocity u ij = u i − u j , density-weighted inter-particle averaged pressurẽ and inter-particle averaged dynamic viscosityη In the following the acceleration contribution of a neighboring particle j to particle i is, for ease of notation, denoted as a ij , where a i = j a ij +b i . Note that the above given momentum formulation, cf. equation (16), exactly conserves linear momentum due to pairwise anti-symmetric particle forces which can easily be verified by using the property ∂W /∂r ij = ∂W /∂r ji of the smoothing kernel.

Equation of state
Following a weakly compressible approach, density ρ i and pressure p i of a particle i are linked via the equation of state with reference density ρ 0 , reference pressure p 0 = ρ 0 c 2 and artificial speed of sound c. Note that this commonly applied approach only represents deviations from the reference pressure, i.e., p i (ρ 0 ) = 0, and not the total pressure. Thus, free boundaries can be modeled by setting p = 0 (see also Section 3.1.6 below). To limit density fluctuations to an acceptable level, while still avoiding too severe time step restrictions, strategies are discussed in [12] on how to determine an appropriate value of the artificial speed of sound.

Boundary conditions
Rigid wall boundary conditions. Following the approach of [13], rigid wall boundary conditions are modeled using fixed boundary particles with quantities extrapolated from the fluid field based on a local force balance.
For more details the interested reader is referred to the aforementioned literature. In the numerical examples in Section 4 the channel walls are modeled using rigid wall boundary conditions.
Inflow and outflow boundary conditions. Open boundaries are modeled similar to [28] via defined inflow and outflow zones occupying so-called inflow respectively outflow particles. Thereby, full support of the interior fluid particles is maintained for density summation (15) and evaluation of the momentum equation (16) when considering contributions from neighboring inflow and outflow particles. At the inflow, i.e., the Dirichlet boundary, the desired inflow velocity is prescribed directly to all inflow particles, while the pressure field is extrapolated from the interior fluid particles i to the inflow particles k following At the outflow, i.e., the Neumann boundary, a zero pressure field is prescribed to all outflow particles. The density field of both inflow and outflow particles is determined from the pressure field with the equation of state (20). Finally, to determine consistent velocities of the outflow particles, the momentum equation (16) is evaluated for outflow particles considering interactions with neighboring fluid particles, boundary particles, and outflow particles.
Periodic boundary conditions. Imposing a periodic boundary condition in a specific spatial direction allows for particle interaction evaluation across opposite domain borders. Moreover, particles leaving the domain on one side are re-injecting on the opposite side. Periodic boundary conditions are commonly applied in SPH modeling of channel or shear flow.

Time integration scheme
The momentum equation (16) is integrated in time applying an explicit velocity-Verlet time integration scheme in kick-drift-kick form, also denoted as leapfrog scheme, as proposed by Monaghan [9]. In the absence of dissipative effects, the velocity-Verlet scheme is of second order accuracy and reversible in time [9].
In a first kick-step the particle accelerations a n i = ( du i /dt ) n determined in the previous time step n are used to compute intermediate particle velocities at n + 1/2 where ∆t is the time step size, before the particle positions at n + 1 are updated in a drift-step Using the particle positions r n+1 i and intermediate velocities u n+1/2 i , the particle densities ρ n+1 i and accelerations a n+1 i are updated following equations (15) and (16). In a final kick-step the particle velocities at n + 1 are determined To maintain stability of the time integration scheme, the time step size ∆t is restricted by the Courant-Friedrichs-Lewy (CFL) condition, the viscous condition, and the body force condition, refer to [12,27] for more details, with maximum fluid velocity u max and maximum body force b max .

Discretization of structural field via the finite element method
The discretization of the structural field, governed by the strong form of the balance of linear momentum (5), is based on the finite element method. Since it is not the focus of this work, the basics of the FEM are presented here only very briefly. For further informations the reader is referred to, e.g., [29,30].
Applying the method of weighted residuals, in the following interpreted as principle of virtual work, the weak form of the initial boundary value problem for the structural field is obtained as with the variation δd s of the primary unknown structural displacement d s . Herein, the contribution to the weak form resulting from the coupling condition of the fluid and structural field at the interface Γ f s (cf. Section 2.3) is omitted and instead treated in Section 3.3.
where H 1 denotes the Sobolev space of functions with square-integrable first derivatives, the weak form (26) is equivalent to the strong form of the balance of linear momentum (5).
The computational domain of the structural field Ω s is sub-divided into non-overlapping finite elements with nodes i. Hence, the structural displacement field d s is discretized introducing nodal displacements d s i of nodes i. The displacement field is approximated via using the Lagrange polynomials N e j with compact support inside element e. Within a Bubnov-Galerkin approach, the same Lagrange polynomials for trial and test functions are employed. Following the isoparametric concept, the parameter coordinates ξ used for the definition of the shape functions within a standard element geometry are mapped onto the physical coordinates applying the same shape functions also used for the displacement interpolation. Specifically, in the numerical examples in Section 4 finite elements based on first-order interpolation are employed.
Subsequently, the semi-discrete form is discretized in time applying a generalized-alpha time integration scheme. The resulting system of nonlinear equations in residual form is finally solved for the nodal structural displacements using a Newton-Raphson method.

SPH-FE interaction: a novel sliding boundary particle approach
In this section, a novel sliding boundary particle approach for the application in a fluid-structure interaction framework coupling SPH and FEM is proposed. In contrast to existing methods modeling boundaries in SPH, e.g., boundary particle methods, cf. Figure 3, the proposed method can handle also deforming interfaces undergoing strong curvature and large stretch, as typical for some FSI applications especially in biomechanics, while keeping the computational costs at a reasonable level. The following is mainly concerned with the evaluation of the interface force f f s at the fluid-structure interface Γ f s . The coupling of fluid and structural field following a Dirichlet-Neumann partitioned approach is subsequently described in Section 3.4.

Conforming interface mesh
Introducing an interface mesh on the fluid-structure interface Γ f s allows for exchange of interface displacement d f s and interface force f f s between the fluid and the structural field, cf. Figure 2, while keeping the fluid and structural solvers separated. For convenience, the interface mesh, which is purely introduced as one possibility to facilitate the displacement and load transfer between the solvers, can be chosen as an extraction or clone of the structural mesh at the fluid-structure interface Γ f s . But the proposed approach also works for non-matching meshes. For the interface mesh, again the iso-parametric concept is employed to describe the standard element geometry of interface elements e via parameter coordinates ξ and Lagrange polynomials N e j of corresponding nodes j. Note that for interface elements the parameter coordinates ξ are of one dimension lower compared to structural elements. Deduced from the geometric coupling condition (8) the current interface position is in the following depicted by r f s . 9 PSfrag replacements fluid particle boundary particle In case of a conforming mesh, the transfer of quantities between interface and structure is straightforward and is, hence, just briefly sketched here. Both, interface position r f s and interface displacement d f s can be extracted directly from the respective structural position r s and structural displacement d s . Similarly, the interface force f f s can be added directly to the respective structural force f s . In case non-matching interfaces are preferred or needed, e.g., because of special resolution demands of the two involved physical fields, the transfer of quantities between interface and structure could simply be done via a Mortar technique [31]. In comparison to the interface structure transfer, the transfer of quantities to the fluid field is more elaborate and will be covered in the following subsections.

Detection of closest projection point
The interaction evaluation is performed between fluid particles and interface elements. Consider a fluid particle i with support radius r c of the smoothing kernel W that is close to the fluid-structure interface Γ f s , cf. Figure 4.
In general, the closest projection point c e i of fluid particle i to the interface Γ f s is located on interface element e and lies within the support radius, i.e., r c e i − r i < r c . The position r c e i of point c e i can be described in iso-parametric coordinates ξ c e i on interface element e. As a result, the shape functions N e j ξ c e i of all nodes j of interface element e evaluated at the closest projection point c e i can be utilized to interpolate kinematic quantities, e.g., positions, velocities, and accelerations, at the closest projection point using nodal quantities and to distribute kinetic quantities, e.g., interaction forces, from the closest projection point to adjacent nodes. The closest projection point c e i of a fluid particle i to a neighboring interface element e is detected solving the following minimization problem with position r i of fluid particle i and positions r f s j of nodes j of interface element e. The solution of the minimization problem gives the iso-parametric coordinates ξ c e i of the closest projection point c e i on interface element e. Hence, the position of the closest projection point c e i results in As stated above only closest projection points c e i located within the support radius of fluid particle i are considered in the interaction evaluation, meaning in addition r c e i − r i < r c must be fulfilled. By definition, when evolving the position of a fluid particle i over time, also the position of the closest projection point c e i is changing, i.e., is sliding on the interface Γ f s . Remark 6. In the general case, the closest projection point of a particle is located on the surface of an interface element, as illustrated for instance in Figure 4. In addition, the two special cases of a convex and a concave angle between two neighboring interface elements are worth being discussed here. In the case a particle is located within the perpendicular straight lines of neighboring interface elements at a convex angle, cf. case 1 in Figure 5, a single closest projection point is considered that is located on the node respectively the edge being shared by those interface elements. For a particle located at a concave angle, cf. case 2 in Figure 5, multiple closest projection points on the surface of each of the interface elements are considered.
Remark 7. Note that very similar to typical contact problems in finite element analysis an extension to a C 1 -continuous representation of the structural geometry, e.g., by employing Hermite polynomials [32] or B-Splines [33,34] as shape functions, could be beneficial within the proposed sliding boundary particle approach in terms of a smoother interaction force evolution and help to abstain from the aforementioned case distinctions, cf. Remark 6.

Virtual boundary particles
The support of the smoothing kernel of a fluid particle i close to the fluid-structure interface Γ f s is truncated, i.e., fluid particle i experiences reduced contributions from neighboring fluid particles, cf. Figure 4. To overcome this issue, full support of the smoothing kernel of fluid particle i is retained by considering a set of virtual boundary particles k e i that contribute to the interaction evaluation of fluid particle i and are regularly and equidistantly arranged behind the closest projection point c e i as illustrated in Figure 4. This is achieved by a certain number of layers of virtual boundary particles with spacing ∆x among them. Accordingly, together with the closest projection point c e i , the set of virtual boundary particles k e i are sliding along the fluid-structure interface Γ f s following the movement of a fluid particle i, giving rise to the name of the proposed method: sliding boundary particle approach.
Remark 8. Within this work, as stated in Sections 3.1.1 and 3.1.2, a quintic spline smoothing kernel with support radius r c = 3h is applied with initial particle spacing ∆x equal to the smoothing length h. As a  All layers of virtual boundary particles are positioned perpendicular to the connection vector r c e i − r i of fluid particle i and its closest projection point c e i , where the first layer is at a distance of ∆x/2 behind the closest projection point c e i on interface element e. An orthonormal basis (e r , e s , e t ) with first base vector e r = r c e i − r i r c e i − r i is constructed [35]. Consequently, the position of all virtual boundary particles k e i can be given in terms of the particle spacing ∆x and the constructed orthonormal basis (e r , e s , e t ) as with integers m r ∈ {0, 1, . . . , (q − 1)} and m s , m t ∈ {−(q − 1), . . . , (q − 1)} where q = floor( r c /∆x ) defines the number of particles necessary to maintain full support of the smoothing kernel. Finally, the vector from fluid particle i to virtual boundary particle k e i is r k e i − r i , cf. Figure 4. Remark 9. The floor operator used herein is defined by floor(x) := max {k ∈ Z | k ≤ x} and returns the largest integer that is less than or equal to its argument x.

Interaction forces on fluid particles
A fluid particle i close to the fluid-structure interface Γ f s , i.e., for which the closest projection point c e i on interface element e is within the support radius r c e i − r i < r c of fluid particle i, additionally experiences contributions to the density summation (15) and the momentum evaluation (16) from all virtual boundary particles k e i for which r k e i − r i < r c holds, cf. Figure 4. As described in Section 3.1.3 the density field is computed via summation of the respective smoothing kernel contributions of neighboring fluid particles j, refer to equation (15). Hence, considering the additional contributions of virtual boundary particles k e i , the density summation for a fluid particle i reads ensuring full support of the smoothing kernel. Inspired by the treatment of boundary particles for rigid walls [13] the properties of virtual boundary particles k e i , i.e., density, pressure, and velocity, are extrapolated based on the corresponding quantities from 12 neighboring fluid particles j of closest projection point c e i on interface element e. The goal is to achieve an undisturbed pressure field of fluid particles close to the interface. Satisfying the kinematic coupling condition on the fluid-structure interface Γ f s , cf. equation (8), also called no-slip boundary condition, viscous forces are considered in the momentum equation. It shall be noted, that some boundary particle formulations in SPH are based on the assumption of zero normal pressure gradients close to the interface. However, in [13] it is shown, that including the pressure gradient obtained from a local force balance is beneficial to accurately model the pressure field of fluid particles close to the boundary. Therefore, a similar strategy is pursued in the following.
In a first step, the pressure p k e i of virtual boundary particles k e i is approximated based on a first order Taylor series expansion with center of expansion at The position r f can be interpreted as smoothed or averaged centroid position of the domain covered by the neighboring fluid particles j as illustrated in Figure 6. Hence, the pressure of virtual boundary particles k e i is determined following with smoothed pressure p f = j p j W c e i j j W c e i j and smoothed pressure gradient The latter is approximated based on a local force balance neglecting viscous forces as proposed in [13], cf. equation (2), with acceleration a c e i of the closest projection point c e i , cf. Remark 11. Applying the equation of state (20) of the respective interacting fluid particle i together with pressure p k e i , the density ρ k e i of virtual boundary particles k e i follows as Remark 10. Note that the approximation of the smoothed pressure gradient ∇p f , cf. equation (34), could be improved considering viscous forces in the local force balance, e.g., [36], however, at the cost of additional computational and algorithmic effort. In a next step, the velocity u k e i of virtual boundary particles k e i is approximated considering the kinematic coupling condition on the fluid-structure interface Γ f s , cf. equation (8), prescribing the velocity u c e i of closest projection point c e i , cf. Remark 11. Applying a first order Taylor series expansion with center of expansion at r f according to (32) gives the relation with smoothed velocity u f = j u j W c e i j j W c e i j and unit vector e r pointing from particle i to closest projection point c e i thus representing the wall normal vector as defined in Section 3.3.3. The quantity ∇ er u f denotes the smoothed directional derivative of the velocity in direction of e r and follows from equation (36) as 13 Sfrag replacements In addition to the acceleration contributions a ij of neighboring fluid particles j, the momentum equation (16) for a fluid particle i is extended by the acceleration contributions a ik e i of virtual boundary particles k e i related to the closest projection points c e i on interface elements e with and density-weighted inter-particle averaged pressurep ik e i as defined in (17). Remark 12. The extrapolation of pressure and velocity for virtual boundary particles k e i by the Taylor series expansions (33) and (38) requires the quantities · f to be evaluated only once for each closest projection point c e i , which is the main advantage of this procedure regarding computational costs.
Remark 13. Note that the contributions of virtual boundary particles k e i resulting from a background pressure p b as part of the transport velocity formulation [27] are also considered for fluid particles i, however, similar to Section 3.1.4 for ease of notation not pointed out here.

Nodal interface forces on interface elements
The coupling of the fluid and the structural field, cf. Figure 2, following a Dirichlet-Neumann partitioned approach (as discussed in more detail in Section 3.4) requires the evaluation of interface forces f f s . To enforce conservation of linear momentum at the fluid-structure interface Γ f s , in accordance with (19), the interface with mass m i of fluid particle i and acceleration contribution a ik e i of virtual boundary particle k e i on fluid particle i. Note that the above given formulation of the resulting force f e

Partitioned coupling approach
The fluid and the structural field are coupled following a Dirichlet-Neumann partitioned approach, where the fluid field is the Dirichlet partition with prescribed interface displacements d f s and the structural field is the Neumann partition subject to interface forces f f s , as illustrated in Figure 2 (right).
Introducing the field operators F and S for the fluid and the structural problem [5] both mapping the interface displacements d f s to interface forces equilibrium at the interface Γ f s is satisfied in case the condition holds. The inverse fluid and structural field operators mapping interface forces f f s to interface displacements are consequently defined as In [7] it is shown that weakly coupled schemes exhibit instabilities in FSI problems with incompressible flows due to the artificial added mass effect. To overcome those instabilities, a fixed-point coupling algorithm is employed to iteratively reach dynamic equilibrium of the fluid and the structural field at the interface with respect to a predefined convergence criterion, i.e., fluid and structural field are strongly coupled. Following a synchronous time stepping scheme the same time step size ∆t is set for both fluid and structural solver and is based on the in general more severe restrictions of the SPH time integration scheme, cf. equation (25).

Remark 14.
Note that the applied generalized alpha time integration scheme for the structural field being an implicit method in general allows for a larger time step size ∆t than possible for the fluid field solved using SPH. Thus, future research may focus on asynchronous time stepping and sub-stepping schemes in order to reduce computational costs.
The coupling algorithm applied herein is described in detail below as Algorithm 1. Convergence of the iterative coupling loop in Algorithm 1 is achieved in case the following criterion based on the increment of interface displacements ∆d f s n+1,i+1 is fulfilled with time step size ∆t, number of interface degrees of freedom n f s dof , and predefined tolerance for convergence ǫ.
Remark 15. In general, applying dynamic relaxation of the interface displacements d f s in each iteration of the coupling algorithm [5] can have a stabilizing effect and accelerate the convergence of the partitioned coupling. However, it shall be noted, that due to the restrictions of the time step size ∆t resulting from the SPH time integration scheme, an accelerating effect is not required with the proposed formulation, cf. examples 4.2.2 and 4.2.3.

Algorithm 1 Time loop of a Dirichlet-Neumann partitioned fixed-point fluid-structure interaction algorithm
∆t n f s dof < ǫ then ⊲ check convergence criterion, cf. equation (46) break end if i ← i + 1 ⊲ increment iteration counter end while n ← n + 1 ⊲ increment step counter end while

Numerical examples
The purpose of this section is to validate the novel sliding boundary particle approach and the proposed numerical formulation for solving fluid-structure interaction problems examining several numerical examples in two and three dimensions. The obtained results are assessed on the basis of analytical solutions and reference solutions given in the literature.

Validation of the sliding boundary particle approach
At first, the capabilities of the proposed method considering fluid flow in the presence of rigid and undeformable structures with a focus on the validation of the novel sliding boundary particle approach as presented in Section 3.3 are shown. The obtained results are compared to analytical solutions and reference solutions given in the literature both in a quantitative and qualitative manner. Additionally, as rigid and undeformable structures are considered, these examples can also be examined utilizing an implementation of the rigid wall boundary condition proposed in [13]. As a result, this allows for validation of the proposed sliding boundary particle approach against the established rigid wall boundary condition [13] within the context of rigid and undeformable structures. Finally, an example demonstrates the advantages of the proposed sliding boundary particle approach in the regime of large structural deformations. In all examples discussed in this section, the structural field is not solved, though, the fluid-structure interface is explicitly described either via an interface mesh or analytically by parameterization.

Hydrostatic pressure in a fluid between two parallel plates
The gap between two spatially fixed and undeforming parallel plates being a distance of L = 0.2 apart is filled with a Newtonian fluid of density ρ f = 1.0 and kinematic viscosity ν f = 1.0 × 10 −2 . A coordinate axis e q is introduced pointing in the direction perpendicular to the parallel plates with origin centered between the latter, cf. Figure 7(a). Finally, a body force of magnitude b q = 0.1 acting in direction e q is applied on the fluid. For this simple example the analytical solution for the pressure profile in the static equilibrium state is given to p(q) = ρ f b q q showing linear behavior.
The fluid domain between the two parallel plates is discretized by 40 layers of fluid particles, i.e., with an initial particle spacing ∆x = 5.0 × 10 −3 . The smoothing length h of the smoothing kernel is set equal to the initial particle spacing ∆x resulting in a support radius r c = 1.5 × 10 −2 . For the fluid phase, an artificial speed of sound c = 1.0 is chosen, leading to a reference pressure p 0 = 1.0. The background pressure p b is set equal to the reference pressure p 0 . The two parallel plates are modeled by a surface element each. The problem is solved with time step size ∆t = 3.125 × 10 −4 , cf. equation (25), until a static equilibrium state is reached.   Figure 8: Hydrostatic pressure in a fluid between two parallel plates: detailed view of boundary region with pressure values of fluid particles and (virtual) boundary particles using the proposed sliding boundary particle approach with (black circles) and without (blue circles) considering the pressure gradient (34) in equation (33), and the rigid wall boundary condition [13] also considering the pressure gradient (red circles) compared to analytical solution (black solid line).  (14). The results are in very good agreement with the analytical solution showing the capability of the proposed sliding boundary particle approach to capture linear pressure profiles near the boundary, cf. equation (33). In addition, the example is computed with an implementation of the rigid wall boundary condition [13] modeling the two parallel plates via fixed boundary particles. Comparing the result to those obtained with the proposed sliding boundary particle approach, cf. Figure 7(b), delivers apart from roundoff errors equivalent results for this example. Finally, a detailed view of the boundary region at q = 0.1 is given in Figure 8 showing the pressure values of fluid particles and (virtual) boundary particles obtained with the proposed sliding boundary particle approach and the rigid wall boundary condition [13]. In addition, a modified variant of the sliding boundary particle approach without considering the pressure gradient (34) in equation (33) is examined. With this modified variant, the pressure value of the fluid particle closest to the boundary, cf. Figure 8, clearly deviates from the expected linear pressure profile. Thus, considering the improved accuracy and the fact that the computational costs required for the extrapolation of pressure (and velocity) for virtual boundary particles are negligible, cf. Remark 12, in the following, the standard variant as proposed in Section 3.3.4 is applied.

Planar Taylor-Couette flow
In this example, a laminar, planar Taylor-Couette flow is considered. The gap between two coaxial cylinders with radii r 1 = 1.0 and r 2 = 2.0 is filled with a Newtonian fluid of density ρ f = 1.0 and kinematic viscosity ν f = 1.0. The inner cylinder is fixed, i.e., its angular velocity is ω 1 = 0.0, while the outer cylinder rotates with angular velocity ω 2 = 2.0 around its axis of symmetry. No-slip boundary conditions are applied between the fluid and the surfaces of the cylinders. The geometry and boundary conditions of the problem are illustrated in Figure 9(a). The Reynolds number of the problem is Re = ω 2 r 2 (r 2 − r 1 ) ν f = 4.0 with maximum velocity ω 2 r 2 and gap (r 2 − r 1 ) between the coaxial cylinders.   The fluid domain is discretized by fluid particles with initial particle spacing ∆x = 5.0 × 10 −2 . The smoothing length h is equal to the initial particle spacing ∆x resulting in a support radius r c = 1.5 × 10 −1 . The artificial speed of sound is set to c = 40.0, hence the reference pressure is p 0 = 1600.0. The background pressure p b is set equal to the reference pressure p 0 . In this example, the structural surfaces are described analytically by parameterization of the cylindrical surfaces in order to show the capabilities and flexibility of the proposed sliding boundary particle approach. However, it shall be noted that the geometry naturally could have been discretized by a finite element mesh. The problem is solved with time step size ∆t = 3.125 × 10 −4 , cf. equation (25), until a nearly stationary state is reached at time t = 2.0.
In Figure 9(b) the stationary velocity of the fluid in the gap between the cylinders, post-processed applying SPH approximation (14), is plotted over the radius r in angular direction θ at time t = 2.0. The result obtained with the proposed sliding boundary particle approach is compared to the result obtained with an implementation of the rigid wall boundary condition [13] and to the analytical solution of the problem [37]. Both methods show a high degree of conformity with the analytical solution. The deviation of the velocity in Figure 9(b) close to the cylindrical surfaces, i.e., at r = r 1 and r = r 2 , results from missing kernel support during post-processing. This phenomenon likewise occurs for both methods, but at a varying degree. Finally, the particle distribution at time t = 2.0 is shown in Figure 10 comparing the results of the sliding boundary particle approach with the rigid wall boundary condition [13]. In contrast to the rigid wall boundary condition with fixed boundary particles approximating the cylindrical shape, the sliding boundary particle approach does not suffer from geometry discretization errors thus resulting in an improved preservation of the solution symmetry (cf. Figure 10) and a decreased deviation from the analytical velocity profile (cf. Figure 9(b)). On the other hand, approaches were the cylindrical shape is discretized by boundary particles in ring-shaped arrangement suffer from a disturbed support of the smoothing kernel of fluid particles close to the cylindrical surface, similar than in Figure 3.
(a) Sliding boundary particle approach with parameterization of the cylindrical shape.
(b) Rigid wall boundary condition [13] with fixed boundary particles approximating the cylindrical shape.

Laminar flow around a rigid cylinder
A prominent CFD benchmark problem was proposed by Schäfer and Turek et al. [38] in the year 1996 and since then was considered in a huge variety of publications. The benchmark is concerned with the laminar flow around a rigid cylinder in a channel. Within this publication, the problem is utilized to validate the momentum exchange at the fluid-structure interface, i.e., at the surface of the cylinder, examining characteristic quantities such as the drag and the lift coefficient or the cycle duration of the time-periodic solution. In the following the focus is set on the two-dimensional, unsteady test case 2D-2 [38].   Figure 11. The channel is filled by a Newtonian fluid initially at rest with density ρ f = 1.0 and kinematic viscosity ν f = 1.0 × 10 −3 . It shall be noted that the problem setup is designed intentionally non-symmetric in order to initiate unsteady vortex shedding behind the cylinder. No-slip boundary conditions are applied at the bottom and top channel wall and on the surface of the cylinder. At the channel inflow, a parabolic, time dependent velocity profile u in = u(x = 0, y, t) is prescribed with components u x (x = 0, y, t) = u max 4y(H − y) H 2 τ (t) and u y (x = 0, y, t) = 0.0 where The maximum inflow velocity is set to u max = 1.5 resulting in a Reynolds number Re = u mean D ν f = 100 with mean velocity u mean = 2/3 u max = 1.0 for all times t ≥ 2.0. At the channel outflow a zero pressure condition p out = p(x = L, y, t) = 0.0 is applied. The fluid domain is discretized by fluid particles with initial particle spacing ∆x = 2.0 × 10 −3 . The smoothing length h is set equal to the initial particle spacing ∆x resulting in a support radius r c = 6.0×10 −3 of the smoothing kernel. An artificial speed of sound c = 12.5 is chosen for the fluid phase leading to a reference pressure p 0 = 156.25. The background pressure is set to p b = 312.5 and is on the order of the reference pressure as proposed by [27]. The bottom and top channel walls are modeled utilizing boundary particles according to [13] with spacing equal to the initial particle spacing ∆x. On account of the fact that the cylinder is fixed and undeformable only the surface of the cylinder is regularly discretized by 48 surface elements of same size that are considered in the computation of the fluid field, i.e., the structural field is not solved. The unsteady flow simulation is solved for times t ∈ [0, 8.0] with a time step size of ∆t = 4.0 × 10 −5 based on the time step size conditions defined in equation (25).
To allow for a quantitative comparison of the obtained results with existing reference solutions, the drag and the lift coefficient are defined as where f drag and f lif t denote the forces in x-respectively y-direction acting on the cylinder obtained from the sum of all force contributions of fluid particles acting on interface elements of the discretized surface of the cylinder, cf. equation (41). Figure 12 shows the drag coefficient c drag and the lift coefficient c lif t obtained for the fully developed time-periodic solution after approximately t = 5.0. Both drag and lift coefficient show typical fluctuations as common in SPH-based simulations (similar to an example in [27]), that result from disturbances of the density field [12] due to relative particle movement. Besides that, the obtained results are in good agreement to the lower bound (c drag = 3.2200, c lif t = 0.9900) and upper bound (c drag = 3.2400, c lif t = 1.0100) of the maximum drag and lift coefficient given by [38]. The shape of the curve of lift coefficient c lif t allows identifying periodic cycles of the solution with approximate cycle duration t cycle = 0.33 in close agreement to the result of [39]. In addition, the example is computed discretizing the cylinder with fixed boundary particles based on an implementation of the rigid wall boundary condition [13]. The results in form of the drag coefficient c drag and the lift coefficient c lif t are compared to those obtained with the proposed sliding boundary particle approach, cf. Figure 12, and likewise show the observed typical fluctuations. Note that the visible phase shift in the time-periodic solution of the lift coefficient c lif t is stemming from roundoff errors that influence the initiation of vortex shedding. Finally, Figure 13 shows the magnitude of the fluid velocity field for a periodic cycle from t 0 = 6.90 to t 1 = 7.23 at four equidistant points in time. At time t = 6.98 the present results of the velocity field visualized in Figure 13 are qualitatively in good agreement to the results of [39]. Altogether, the results of the CFD benchmark problem obtained with the sliding boundary particle approach represent the given reference solutions [38,39] both quantitatively and qualitatively in good approximation and further showcase the capabilities of the novel formulation to accurately model the momentum exchange at the fluid-structure interface.  Figure 12: Laminar flow around a rigid cylinder: drag coefficient c drag and lift coefficient c lif t using the proposed sliding boundary particle approach (black solid line) and the rigid wall boundary condition [13] (red dashed line) compared to the upper bounds given in reference solution [38] (blue dashed line).

Isochoric deformation of a box filled with a fluid
This example aims to demonstrate the advantages of the proposed sliding boundary particle approach over fixed (material) boundary particle methods in the case of large deformations at the fluid-structure interface. To this end, an academic example is examined utilizing both the proposed sliding boundary particle approach and an implementation of the rigid wall boundary condition [13] with boundary particles fixed to material points of the structure.
An initially quadratic structural box with inner edge length b = 0.1 and wall thickness d = 0.015 is filled by a Newtonian fluid initially at rest with density ρ f = 1.0 and kinematic viscosity ν f = 1.0 × 10 −2 . An isochoric deformation of the structural box to obtain a rectangular shape (with final edge lengths b x = 0.2 and b y = 0.05 starting from t = 2.5) is prescribed, defined by the deformation gradient and accordingly with det F = 1.0. It follows, that also the volume of the fluid within the structural box remains constant at all times. Consequently, in the final static equilibrium state the fluid density is expected to be constant throughout the entire fluid domain. The fluid domain within the structural box is discretized by fluid particles with initial particle spacing ∆x = 5.0 × 10 −3 . The smoothing length h of the smoothing kernel is set equal to the initial particle spacing ∆x resulting in a support radius r c = 1.5 × 10 −2 . For the fluid phase, an artificial speed of sound c = 1.0 is chosen, leading to a reference pressure p 0 = 1.0. The background pressure p b is set equal to the reference pressure p 0 . The walls of the structural box are either modeled by surface elements when using the proposed sliding boundary particle approach or by boundary particles fixed to material points of the structure when using the rigid wall boundary condition [13]. The problem is solved for times t ∈ [0, 10.0] with time step size ∆t = 3.125 × 10 −4 , cf. equation (25). Figure 14 shows the particle distribution obtained using the rigid wall boundary condition [13] with fixed (material) boundary particles in the initial configuration and at time t = 10.0. Presribing the deformation of the structural box naturally also distorts the initially regular arrangement of boundary particles fixed to material points of the structure, i.e., the boundary particle spacing is streched in horizontal direction and compressed in vertical direction, which is clearly visible at time t = 10.0. As a consequence, the support of the smoothing kernel of a fluid particle close to the interface is disturbed, also influencing the density (and pressure) field in the deformed fluid domain. Eventually, leakage of fluid particles through the fluidstructure interface occurs when the boundary particle spacing becomes too large, and accordingly, the fluid density within the structural box is significantly reduced with an average density error of approximately 7.5 %. The results obtained with the proposed sliding boundary particle approach are shown in Figure 15.
For the purposes of illustration, at time t = 10.0 the virtual boundary particles belonging to a fluid particle close to the upper edge and to a fluid particle close to the right edge are shown. Here the full benefits of the proposed sliding boundary particle approach become obvious: full support of the smoothing kernel of a fluid particle close to the interface is retained by a transient set of regularly arranged virtual boundary particles. As a result, an undisturbed density (and pressure) field is achieved in the deformed fluid domain, and consequently, no leakage of fluid particles through the fluid-structure interface occurs. Altogether, this example illustrates the advantages of the proposed sliding boundary particle approach over fixed (material) boundary particle methods when considering large deformations of the fluid-structure interface.
(a) Initial configuration with fluid particles (gray) and boundary particles (black).
(b) Final configuration with fluid particles colored with fluid density ranging from 0.9 (blue) to 1.1 (red) and boundary particles (black). Figure 14: Isochoric deformation of a box filled with a fluid: particle distribution obtained using the rigid wall boundary condition [13] with fixed (material) boundary particles at points in time t 0 = 0.0 and t 1 = 10.0.
(a) Initial configuration with fluid particles (gray) and surface elements (blue).
(b) Final configuration with fluid particles colored with fluid density ranging from 0.9 (blue) to 1.1 (red) and an illustration of the virtual boundary particles (black) belonging to two characteristic fluid particles. Figure 15: Isochoric deformation of a box filled with a fluid: particle distribution obtained using the proposed sliding boundary particle approach at points in time t 0 = 0.0 and t 1 = 10.0.

Validation of the fluid-structure interaction framework
Additional complexity is added by considering freely moving and deformable structures stressing the coupling of fluid and structural field following a Dirichlet-Neumann partitioned approach. Consequently, in the following examples also the structural field is solved. Analytical solutions and reference solutions given in the literature are used to validate the obtained results in quantitative and qualitative manner.

A rigid cylinder floating in a shear flow
The following example is based on the studies [40,41] stating that a rigid cylinder floating in a shear flow in a channel always migrates to the center of the channel independent of its initial position and velocity.
Here, this example serves as a further validation of the proposed method considering rigid body motion of the structural field. For validation, the obtained results are compared to [36] where both the fluid and the solid field are discretized using SPH.
A  In this example, the fluid domain is discretized by fluid particles with initial particle spacing ∆x = 1.0×10 −4 . The smoothing length h is equal to the initial particle spacing ∆x leading to a support radius r c = 3.0 × 10 −4 of the smoothing kernel. An artificial speed of sound c = 0.25 is chosen, resulting in a reference pressure p 0 = 0.0625 for the fluid phase. The background pressure p b is set equal to the reference pressure p 0 . The motion of the bottom and top channel walls is modeled using moving boundary particles according to [13]. A Saint Venant-Kirchhoff model with relatively high Young's modulus E s = 1.0 × 10 6 and Poisson's ratio ν s = 0.4 is applied for the structure in order to penalize deformation of the cylinder and allow primarily rigid body motions. The cylinder is regularly discretized by 144 first-order elements with 48 surface elements on the surface of the cylinder. Convergence of the iterative coupling algorithm is checked based on the tolerance ǫ = 1.0 × 10 −8 in equation (46). The problem is solved for times t ∈ [0, 60.0] with a time step size of ∆t = 1.0 × 10 −4 .
The vertical position of the center of the cylinder in the channel over time t is displayed in Figure 17. The cylinder migrates to the center line of the channel as expected. In addition, a quantitative comparison to the results given in [36] shows good agreement for the dynamcis of the solution.

Flow-induced oscillations of a flexible beam attached to a rigid cylinder
Based on the benchmark problem of a laminar flow around a rigid cylinder in a channel [38], cf. Section 4.1.3, a FSI benchmark was proposed by Turek and Hron [42] as modification of an example first described in [43]. The purpose of the example is to study flow-induced oscillations of a flexible beam attached to a rigid cylinder in a channel flow. In the following the two-dimensional test case FSI2 [42] is considered that is characterized by large structural displacements.
The problem setup (rectangular channel of length L = 2.2 and height H = 0.41, rigid cylinder of diameter D = 0.1 with center fixed at position (0.2, 0.2)) is very equal to the example of a laminar flow around a rigid cylinder discussed in Section 4.1.3. In addition, in this example a flexible beam of length l = 0.35 and height h = 0.02 is attached at the downstream end of the rigid cylinder, i.e., at position (0.25, 0.2) as illustrated in Figure 18. Note that also the length of the channel remains equal to the example in Section 4.1.3 and is thus slightly shorter than originally proposed for the benchmark [42]. A control point needed for evaluation of the results, e.g., in Figure 19  The fluid domain is discretized with fluid particles similar than in the example in Section 4.1.3 (initial particle spacing ∆x = 2.0 × 10 −3 , support radius r c = 6.0 × 10 −3 of the smoothing kernel, artificial speed of sound c = 12.5, reference pressure p 0 = 156.25). In this example, the background pressure is set to p b = 1250.0. Boundary particles according to [13] are utilized to model the bottom and top channel walls. The flexible beam as part of the structural domain is discretized by 35 × 3 first-order elements. The surface of the cylinder exposed to the fluid field, i.e., without considering the part where the flexible beam is attached, is discretized by 20 surface elements. Convergence of the iterative coupling of fluid and structural field is based on the tolerance ǫ = 1.0 × 10 −8 , cf. equation (46). The FSI problem is solved for times t ∈ [0, 12.0] with a time step size of ∆t = 4.0 × 10 −5 . In this example, convergence of the partitioned coupling loop, cf. Algorithm 1, is reached after an average number of approximately 5.45 iterations per time step, when averaging over all time steps of the given problem.
The vertical displacement d y of the control point at the tip of the flexible beam is displayed in Figure 19. In the present results the minimum and maximum displacement of the control point at the tip of the flexible beam in y-direction are approximately −0.08269 and 0.08309. This is in good agreement with the results given in the literature: [42] and [44] report a minimum and maximum displacement of −0.07937 and 0.08183 respectively −0.0803 and 0.0829. The solution of the FSI problem shows time-periodic cycles of the beam deflection after approximately t = 8.0 with a cycle duration t cycle ≈ 0.525 and a frequency f = 1/t cycle ≈ 1.905 (averaged over all time-periodic cylces), which is in good agreement with [42,44] (f = 1.90). In Figure 20 the magnitude of the fluid velocity field and the deformation of the structure for a periodic cycle from t 0 = 10.32 to t 1 = 10.84 at four equidistant points in time are shown. Especially, at times t = 10.45 and t = 10.71 the flexible beam experiences strong curvature. This is were approaches discretizing the structural domain by boundary particles fixed to structural material points suffer from a disturbed support of the smoothing kernel of neighboring fluid particles, cf. Figure 3. The novel sliding boundary particle approach by definition is not prone to that issue. In summary, the results of the FSI benchmark problem obtained with the sliding boundary particle approach are both quantitatively and qualitatively in good agreement with the given reference solutions [42,44].  Figure 19: Flow-induced oscillations of a flexible beam attached to a rigid cylinder: vertical displacement dy of the control point at the tip of the flexible beam using the proposed sliding boundary particle approach (black solid line) compared to the minimum and maximum displacements given in reference solution [44] (blue dashed line).

Inflation of an academic balloon-like problem
The filling process of a highly flexible thin-walled balloon-like container undergoing large deformations is studied in this example, representing a model problem close to potential application scenarios of the proposed scheme in the field of biomechanics.
An initially cubical structural geometry with inner edge length B = 3.0 and wall thickness d = 0.2 is inflated via a quadratic inlet of width and length b = 1.0 by a Newtonian fluid that is initially at rest with density ρ f = 1.0 and kinematic viscosity ν f = 5.0 × 10 −1 , cf. Figure 21(a). The constitutive behavior of the structure with density ρ s 0 = 1.0 is described by a Saint Venant-Kirchhoff model with Young's modulus E s = 1.0 × 10 2 and Poisson's ratio ν s = 0.45. A similar problem was first proposed in [45] with the purpose to study and solve the incompressibility dilemma in partitioned fluid-structure interaction with pure dirichlet fluid domains. This dilemma does not exist in our approach, given that SPH uses a weakly compressible approach. Here, the example is recapitulated on a three-dimensional geometry with modified fluid and   is prescribed. Note that the origin of the coordinate system (x, y, z) is located at the bottom left corner of the inflow area. Accordingly, the volume inside the academic balloon (without considering the volume of the inlet) can be determined analytically via for each time t with initial volume V 0 = B 3 and inflow area A in = b 2 . The fluid domain is discretized by fluid particles with initial particle spacing ∆x = 4.0 × 10 −2 . The smoothing length h is set equal to the initial particle spacing ∆x resulting in a support radius r c = 1.2×10 −1 of the smoothing kernel. The artificial speed of sound is set to c = 40.0, hence the reference pressure is p 0 = 1600.0. The background pressure p b is equal to the reference pressure p 0 . The walls of the fixed inlet are modeled utilizing boundary particles according to [13] with spacing equal to the initial particle spacing ∆x. The balloon-like structural domain is discretized by first-order elements with a cubic shape in the initial configuration and a characteristic element length of d resulting in one element over the wall thickness. This discretization is justified since the focus of this example is set on the coupling of fluid and structural field at the interface rather than a precise prediction of structural quantities such as the deformation field of the tank. The tolerance ǫ = 1.0×10 −8 in equation (46) is applied for the iterative coupling of fluid and structural field. The FSI Problem is solved for times t ∈ [0, 14.0] with a time step size of ∆t = 2.5 × 10 −4 . In this example, convergence of the partitioned coupling loop, cf. Algorithm 1, is reached after an average number of approximately 5.88 iterations per time step, when averaging over all time steps of the given problem.
The volume inside the academic balloon is determined summing up the effective volumes of respective particles j following V = j m j /ρ j and compared to the analytical solution (53), cf. Figure 21(b). The resulting volume is slightly below the analytically determined volume. This can be explained by the weakly compressible approach, cf. Section 3.1.5, applied in this SPH formulation leading in this example to a minor compression of the fluid phase with an average density error of approximately 1 %. Note that conservation of mass and accordingly (within the limits of a weakly compressible approach) conservation of volume is a characteristic property of SPH. By definition, this means that also the number of fluid particles is conserved. Therefore, the obtained results, among others, demonstrate that no leakage of fluid particles through the fluid-structure interface occurs. The rear half of the (inflated) structural geometry and a quarter section of the fluid domain are displayed in Figure 22 in the initial state and at time t = 12.65. The fluid velocity is post-processed applying SPH approximation (14). Note that at time t = 12.65 the volume inside the academic balloon has doubled, cf. analytical solution (53). In conclusion, this example is characterized by large structural deformations in form of strong curvature and stretch reaping the full benefits of the proposed sliding boundary particle approach in contrast to fixed (material) boundary particle methods. In the presence of large structural deformations, the latter class of boundary particle methods is characterized by insufficient kernel support, and eventually such methods are prone to leakage of fluid particles.

Conclusion and outlook
A novel smoothed particle hydrodynamics (SPH) and finite element (FE) coupling scheme for fluidstructure interaction, the sliding boundary particle approach, is presented in this publication. The coupled problem is solved via a Dirichlet-Neumann partitioned approach, with the fluid field (discretized via SPH) being the Dirichlet partition and the structural field (discretized via FE) being the Neumann partition. SPH is a mesh-free computational method that simplifies the treatment of both large deformations in the fluid domain as well as complex flow while avoiding additional methodological and computational effort compared to fully mesh-based methods. Introducing the sliding boundary particle approach for the treatment of deformable and strongly curved boundaries of the SPH domain in an accurate, robust, and computationally cheap manner, constitutes an important aspect of the proposed numerical formulation for solving FSI problems.
Several numerical examples showcase the capabilities of the novel numerical formulation. To begin with, the sliding boundary particle approach is validated examining two-dimensional examples driving certain characteristics of the proposed formulation. The numerical results obtained for the examples of a hydrostatic pressure in a fluid between two parallel plates, cf. Section 4.1.1, and a planar Taylor-Couette flow, cf. Section 4.1.2, are in very good agreement with the respective analytical solutions confirming the capability of the proposed method to model linear pressure profiles near the boundary and to account for no-slip boundary conditions at the boundary as required for high accuracy of the fluid velocity field. In a next step, numerical examples involving dynamic effects and large structural deformations are studied confirming the accuracy and robustness of the proposed formulation. This is, among others, demonstrated showing the results of well-known CFD respectively FSI benchmark problems, cf. Sections 4.1.3 and 4.2.2, as proposed in [38,42]. Altogether, the obtained numerical results are in very good agreement with the results given in the literature. Finally, a three-dimensional, application-focused example is considered examining the filling process of a highly flexible thin-walled container (cf. Section 4.2.3).
Future work may focus on an asynchronous time stepping scheme, e.g., a sub-cycling scheme of fluid and structural field, cf. Remark 14. Such an approach would allow to evolve the solution of the sub-fields with different time step sizes each best suitable for the underlying method respectively solver while reducing the overall computational effort. Besides that, the FSI framework may be extended to multiphase flow including the motion of rigid bodies. The framework developed herein will be a valuable tool for detailed studies of biomechanical problems involving complex flow, e.g., the human stomach during digestion [1,2].