Radial Basis Functions Vector Fields Interpolation for Complex Fluid Structure Interaction Problems

: Fluid structure interaction (FSI) is a complex phenomenon that in several applications cannot be neglected. Given its complexity and multi-disciplinarity the solution of FSI problems is difﬁcult and time consuming, requiring not only the solution of the structural and ﬂuid domains, but also the use of expensive numerical methods to couple the two physics and to properly update the numerical grid. Advanced mesh morphing can be used to embed into the ﬂuid grid the vector ﬁelds resulting from structural calculations. The main advantage is that such embedding and the related computational costs occur only at initialization of the computation. A proper combination of embedded vector ﬁelds can be used to tackle steady and transient FSI problems by structural modes superposition, for the case of linear structures, or to impose a full non-linear displacement time history. Radial basis functions interpolation, a powerful and precise meshless tool, is used in this work to combine the vector ﬁelds and propagate their effect to the full ﬂuid domain of interest. A review of industrial high ﬁdelity FSI problems tackled by means of the proposed method and RBF is given for steady, transient, and non-linear transient FSI problems.


Introduction
The interaction between a fluid and a structure (FSI) is a complex phenomenon involving several disciplines and physics. The governing laws of the structure and of the surrounding fluid are coupled by means of the boundaries, when deformable or moving. Pressures and forces are indeed exerted by the fluid to the wetted surfaces, generating stress and strain fields causing a shape change of the structure that, depending on the magnitude of the deformations, can influence back the fluid flow in a closed loop. When dealing with steady problems, the interaction between fluid and structure consists in finding the configuration in which the elastic forces and the loads exerted by the fluid are in equilibrium. In transient simulations by the other hand, the effects of inertia must be taken into account, and the interaction between dynamically moving structures and transient waves in the flow can bring to very complex couplings in which slight variations in the numerical methods adopted to solve the FSI problem can make the difference between faithful and diverging simulations. Given this, the scientific community has devoted a lot of attention over the years to this problem, proposing several approaches to tackle FSI at different levels of detail. Several books and reviews are available in literature, among the most interesting are the works by Bazilevs et al. [1], Morand and Ohayon [2], Mittal and Iaccarino [3], and others [4][5][6]. Numerically, available FSI methods can be broadly divided in two groups according to the degree of integration of the two physics involved. When the fluid and structural systems are solved in the same framework, with a mutually strong integration in terms of equations, the approach is called monolithic [7]. On the other hand, when the fluid and structural problems are managed by different independent codes the approach is called partitioned [8]. In both cases there are advantages and drawbacks. The monolithic approach can be the most stable and powerful, since both physics can be intertwined having direct access to the values of interest. Being able to solve at the same time the computational fluid dynamic (CFD) and computational structural mechanics (CSM) problems, allows to rely on more refined and sophisticated approaches, taking advantage of several possible numerical methods to couple the two physics. Stability can be greatly improved over partitioned approaches since data exchange is continuously performed. These positive characteristics are however achieved at the cost of a more difficult development and coding. Software maintenance can be problematic given the largest codebase and monolithic methods are often more computationally demanding. On the other hand the partitioned approach has the huge advantage of allowing the use of well-established and reliable software, solving, separately, fluid and structural problems. The drawback in this case is, however, the cost associated with the data exchange required by the coupling, that can be complex and expensive. Information, as a matter of fact, must be exchanged at each FSI iteration between CFD and CSM codes, transferring the loads calculated from the CFD to the structural model, and the resulting displacements back from CSM to the fluid dynamic numerical grid. This task is sensitive, especially because the requirements in terms of refinement, element typology, and dimension of the computational grid depend on the particular solver, and for this reason CSM and CFD meshes do not generally match in the sense that numerical grid elements, nodes, connectivity, and data locations are different. There is, for this reason, the need for numerical algorithms able to properly interpolate back and forth loads and displacements, with a proper equilibrium between accuracy and computational costs while respecting the requirements in terms of load balance when dealing with forces, precision and smoothness when dealing with displacements. Over the years, several mapping methods have been developed and are available in literature [9][10][11], such as the standard mortar method (SMM) [12] and the force reaction method [13]. Another drawback given by data transfer for partitioned approaches, when using commercial CSM and CFD solvers, is also given by the technical problem of exchanging information since, when using a closed source program offering limited degree of customization or scripting, transferring data by means of external files is often the only available option that can turn into a major bottleneck in terms of computational cost. Another important task, that can be a major challenge, is also given by the update of the CFD numerical grid according to imported displacements. To face this task several methods were proposed in literature and, among them, two proved to be suitable for FSI problems: the immersed boundary (IB) method and the arbitrary Lagrangian-Eulerian (ALE) method. IB was successfully employed for solving complex CFD problems involving flow interaction with non-trivial geometries [14,15]; IB method was also successfully applied to FSI simulation of prosthetic aorta valve [16][17][18][19][20] and successfully employed with variational transfer for coupling the fluid and solid subproblems [21]. The ALE method takes advantage of an hybrid description of finite elements, in which the grid is not fixed (as in Eulerian approach) or attached to the material (as in the Lagrangian approach), but has an intermediate velocity value [22,23]: during the simulation, internal nodes can move in order to optimize their shape as the simulation evolve, and boundary elements can move following material displacement. Both IB and ALE characteristics allow to easily update the fluid-solid interface configuration during FSI simulation, but one of the major drawbacks of these methods is that mesh elements can suffer from high distortion and quality degradation, which will affect both numerical solvers performances and stability. Another viable approach to CFD grid update can be given by domain remeshing, but its feasibility is limited by the computational cost required and by the remeshing noise introduced. Automatic remeshing, moreover, can introduce reliability problems that can influence a good outcome of the simulation. For this reason, a powerful approach to this problem can be given by mesh morphing. By propagating mesh deformations in the volume with morphing allows to reduce the computational cost, since it is faster than remeshing and-maintaining the mesh topology unaltered-it is possible to hot restart a CFD simulation on the deformed mesh from the previous CFD solution, further accelerating the calculation. Several mesh morphing methods exist and are employed in commercial software, among them the most common are the free-form deformation (FFD) [24,25] and the boundary displacement method or pseudo-solid [26,27]. Another class of methods, that is gaining traction over the years, is based on radial basis functions (RBF). RBF are a special class of functions able to interpolate everywhere in the space scalar values given at points. This ability is ideal for mesh morphing, since displacements must be retrieved in the volume when applied as boundary conditions on wetted surfaces. By prescribing scalar values at points, there is a higher degree of freedom and precision with respect to other methods, such as FFD, that do not operate directly on the surfaces. Using this technology, a linear problem must be solved to obtain the interpolation field, but the result can be stored and quickly employed when required. Furthermore, being RBF-based morphing a meshless method, meaning that a continuous deformation field in the space is generated, the same shape variation can be employed for different meshes, regardless of their different mesh refinement or element typology. If compared to IB or ALE methods, the setup of the shape modification in RBF based mesh morphing can request more user efforts. Nevertheless, these efforts will generate less distorted modified meshes and with higher element quality: this will result in a lower impact on CFD solver performances and stability. Meshless characteristics of RBF-based mesh morphing can be efficiently exploited for FSI simulations using the modal approach. According to this method, the modal shapes of the deformable surfaces are at first computed by means of the CSM solver and stored as shape variations. This approach moves the problem on the fluid solver side. Both for the linear case, where modal shapes and their frequency represent the basic nature of the dynamic and static behavior of a structure, and in the non-linear one, where keyframes of the deformation evolution are captured from the FEA run, the combination of applied vector fields can be properly over imposed to adapt the full CFD grid to the required deformation state. By importing the RBF vector fields into the CFD solver, their combination can be tweaked according to a physical law based on the loads retrieved during the calculation, turning the CFD mesh into an intrinsically aeroelastic model for the case on linear structures, or to smoothly fade across the known key configurations. In this way, the expensive interpolation task required at each iteration of the CFD run is avoided and the volume mesh adapted to the deformed shape. It is worthwhile to note that vector interpolation work to map the required vector fields onto the CFD mesh is performed once, as the shapes are all known in advance, and during the run the shape evolution is achieved acting only on the scalar weights that are imposed to apply each vector field with the required intensity. The effectiveness of such approach [28] has been demonstrated by controlling the movement of the walls for the case of a prosthetic aortic valve, driving the movement of the walls according to FEA computes time history [29] and for the case of an ascending aorta where the movement of the vessel wall is imposed by in vivo acquired evolving shape [30]. This paper is conceived as a "walkthrough" of the possible approaches that can be followed, using the modal superposition method based on RBF morphing, for several different classes of FSI problems. This work is not thought of as a general review of the FSI methods available in literature, but a recollection of FSI problems successfully solved using the proposed method by the authors. For each application a very brief introduction is given, redirecting the reader interested in the problem to published papers in which more information is available. From this point of view, this work is an update of a previous paper presented by the same authors [31] to which new cases and approaches were added.
Three classes of challenging industrial cases will be tackled: steady problems, transient FSI cases with structural movement prescribed in advance and fully transient simulations.
The paper is arranged as follows: at first the material and methods at the core of the approach are introduced to the reader. The RBF mathematical framework is presented together with the theory behind modal superposition and non-linear time history shape evolution, then the workflow and related information are given on how to properly couple the tools to tackle FSI cases. To close the work, several industrial applications and examples of use of the proposed method are presented.

RBF Background Theory
RBF are a mathematical tool able to interpolate and extrapolate everywhere in the space scalar values defined at points, we call centers or source points. The resulting field is retrieved on a distance basis, and its calculation entails the solution of a linear system of order equal to the number of centers [32]. A single solution can be performed for a set of source points, once, to compute the coefficients of the linear system that can be used later to directly retrieve interpolated or extrapolated values elsewhere in the space as the summation of their radial contribution. The RBF can be optionally enriched by using a polynomial term: in which the interpolant s is composed by a radial function ϕ and a polynomial h, that allows to retrieve analytically the polynomial fraction of the interpolation and rigid motions. h has order m − 1, where m is the order of the RBF function ϕ. N is the number of source points having, each point i, position x k i . As it can be seen in (1), the interpolant s(x) for any point in the space x depends on the radial contribution of each source point i, function of the Euclidean norm between the required point in the space and the coordinate x k i . The polynomial h can be written as: and so Equation (1) can be expanded as: The solution of the system is found by enforcing that the RBF gives the exact prescribed value g i at each source point i, and that polynomial term satisfies the orthogonality conditions: for all polynomials q with a degree less than or equal to that of polynomial h ( [33]). As a result of the interpolation, the RBF coefficients γ i and polynomial terms β i are obtained. Depending on the RBF kernel chosen, the order of h varies. As shown in [34], if the RBF is conditionally positive definite, a unique interpolant exists. The linear system required to compute γ i and β i can be written in matricial form where U is the distance matrix in which all the mutual distances are computed as and P is the matrix required by the polynomial interpolation in the form Since RBF interpolation allows to retrieve scalar values, to interpolate the vector fields required by mesh morphing in the three dimensions, a system can be employed to calculate for each point in the space the components along x, y, and z: The behavior of the interpolation between source points depends on the kind of basis function employed. A first categorization can be performed depending on the sphere of influence of each source point: with globally supported functions each source point influences the interpolation everywhere, also at significant distances, but with a law based on the Euclidean norm. With compactly supported functions, on the other hand, each source point has an influence only inside a given support radius R sup . Compactly supported functions, such as Wendland's functions, bring to a sparse distance matrix U, allowing the use of fast linear sparse solvers. In Table 1 some of the most used RBF kernels, globally and compactly supported, are shown.

RBF with Compact Support
As previously described, an ideal field of application for RBF is mesh morphing. Being able to interpolate in the space a scalar value defined at discrete points, allows to retrieve a vector field imposed on chosen boundaries. Shape variations can be then smoothly propagated on the computational mesh to surfaces or through volume elements. For transient and steady FSI of linear structures, the modal shapes obtained by means of a FEM simulation can be imported in terms of source points and their displacement for the wetted surfaces, and enriched adding fixed source points to deform the computational mesh only where required leaving fixed other boundaries. A similar approach can be used for imposing to the fluid mesh the evolution of a non-linear structure, with the difference that, in this case, the incremental vector fields to move from a keyframe to the subsequent one are applied. The resulting RBF problem can be linearly amplified by means of an amplification factor and several vector fields over imposed to obtain complex shapes. With this approach the CFD grid is made parametric at the beginning; all the vector fields are in fact interpolated by RBF onto the full volume mesh at the beginning and so, at run time, the deformed mesh is obtained by a simple span of the pre-computed information. In this work, the RBF-based morphing required by the modal superposition method was accomplished using the commercial morpher RBF Morph add-on in ANSYS Fluent. In Figure 1, an example of a possible setup for an RBF problem is shown on a HIRENASD model. Displacements, extracted from a FEM eigenvalue simulation and relative to the first modal shape, are applied in this case as a vector field to the wing wetted surface. The fuselage is left undeformed by imposing a set of source points on its surface with zero motion, and a buffer is left in between the two sets in order to allow a smooth blending between moving and fixed surfaces. The volume domain is finally wrapped by a set of points arranged on a bounding box with zero motion, in order to reduce the morphing field only to the mesh elements inside it. In Figure 1, it is also visible the geometry after morphing, with an amplification factor of 0.5.
RBF-based mesh morphing exhibits features that make the approach attractive with respect to other well-known algorithms, such as the already cited FFD and boundary displacement method. At the core of FFD methods there is a parametric representation of the volume wrapping the deformable area, using geometrical entities described by means of Bezier basis functions [25], but also b-splines [35], NURBS [36] and t-splines [37]. Mesh morphing is performed by moving the control points used to define the Bernstein polynomials, obtaining a cost efficient and smooth deformation of all the points wrapped by the parametric volume. This method, however, lacks of an accurate pointwise control of the surfaces, since there is not a direct action on the geometry but on control points scattered at a certain distance, making this method unfeasible to be used to transfer back and forth precise nodal displacements between CFD and CSM solvers. By the other hand, the pseudo-solid or boundary displacement method allows to modify the shape of the computational grid with an high pointwise level of detail, by solving a structural auxiliary model in which the modified geometry is obtained by means of prescribed displacements or applying ad-hoc loads [38]. Displacements are propagated in the volume as the result of the auxiliary structural deformation that can be finely tuned, controlling the behavior by using different material stiffnesses on an element by element basis [26]. To improve the deformation definition advanced approaches, such as the use of rigid elements [39] or hybrid techniques [40] are suggested in literature. Although this method, in contrast to FFD, exhibits a pointwise accurate control on surfaces, some downsides can be encountered with respect to RBF. Each time that a shape variation must be generated, the mathematical systems required for the solution of the elastic problem must be solved, building and decomposing the reduced stiffness matrix. This expensive task can be inefficient or unfeasible when dealing with the very large meshes that are often required by CFD. Fluid dynamic computation requires several mesh features, such as boundary layers and wake refinement, that increase cell and node counts up to hundreds of millions, dimensions that cannot be easily processed by means of a structural solver. As an example, for a CFD mesh composed of 6 millions of cells, the boundary displacement method requires 4 h for a single shape modification, while RBF-based morphing only 5 min. Moreover, some elements used commonly in CFD, such as polyhedral or mosaic elements, are not implemented generally in FEM structural solvers, requiring a complex conversion back and forth to hexahedral or tetrahedral elements. RBF-based morphing overcomes FFD and pseudo-solid limitations, by providing a pointwise accurate action where desired with a meshless flavor, allowing to prescribe displacements at node level if needed without the requirements of matching nodal positions. From a computational point of view, with RBF meshes of hundred of millions of elements can be easily tackled with thousands of source points, grouped near the area to be modified [32,41]. To understand how the quality of the mesh is preserved using the pseudo-solid and RBF methods, two test cases were explored and compared. At first a simple 2D plate with a circular hole was deformed [41] by translating the hole while maintaining fixed the boundaries, as shown in Figure 2, in which the nodal displacements achieved using the boundary displacement method are compared with those obtained using RBF. Mesh quality is inspected for several degrees of displacement of the central hole, in terms of the worst elemental quality and mean quality averaged on the whole mesh, as shown in Figure 3. As it can be seen, for this case the FEM solution propagates smoothly the deformation on the mesh, obtaining a better mean quality through the mesh when increasing the displacement. In the most deformed elements, however, the RBF method helps to preserve a better quality with respect to FEM, a crucial aspect when dealing with CFD computations. Another interesting test case comparing the boundary displacement method and RBF, was published in [42], in which A 3D problem dealing with the displacement of a cube inside a wind tunnel was investigated. In Figure 4 the comparison of the ground deformation for a prescribed displacement of the cube is shown for the pseudo-solid (left) and RBF methods (right). The elemental deformation introduced near the corners is clearly visible in the pseudo-solid case. Being RBF based on distance on the other hand, the elements closer to source points tend to maintain elemental ratio and shape, feature especially useful when dealing with boundary layers that must be kept as much undeformed as possible. In Figure 5 results in terms of mean skewness and element quality are shown for the cube test case. Differently to the holed plate, in this case RBF performs better for all the range of displacement, given the complexity introduced by the right angles.

Parametric Mesh Formulation
The RBF interpolant defined in Equation (9) can be used to update position for each node of the CFD mesh, since the only input needed are the original node position and the RBF interpolant itself By computing the appropriate amplification during CFD computation and including them into the CFD solver, it is possible to make the mesh parametric, since it can deform under CFD loads thanks to the modal superposition implementation. In the next sections it will be illustrated how the deformed shape of the fluid immersed bodies can be represented linearly combining structural modes, each weighted by the relative modal coordinate η m , see Equation (15), or by smoothly adding the contribution of intermediate incremental displacement A i , see Equation (24). In this way, the CFD mesh is parametric with respect to the shapes. Since these are evaluated accordingly to the modal linear theory, it is possible to evaluate them once and the use of the linear combination is preferred to the RBF Formula (10). In Equation (11) (15) and (24)) and ∆u m are the vectors of the m-th shape: these are precomputed once at the beginning thanks to Equation (10) and are stored in memory to be used in CFD mesh updating.

Background to Modal Analysis and Modal Superposition
In structural mechanics, modal analysis is used to evaluate undamped free vibrations modes of a system, each associated to a natural shape and a natural frequency. One of the use of modal analysis results in FEM modeling, is the evaluation of the system structural response, both static and dynamic and under linearity conditions. If no damping is acting on the system, the nodal amplitudes (i.e., modes of the structure) u and relative natural frequencies can be evaluated by solving the well-known eigenvalue problem of the system where K and M are, respectively, stiffness matrix and mass matrix, ω 2 is an eigenvalue, ω is a natural frequency. System in Equation (12) explains that a natural mode (or vibration mode) is a structure deformed configuration in which elastic resistance and inertial load are in a balance ( [43]). Mechanical systems behave as low-pass filters, modes associated with lowest frequencies have highest energy levels and are most important since from a physical point of view they take over the other ones. Therefore, to perform modal superposition analyses, a good approximation is achieved using a reduced set of low frequency modes, obtaining thus a reduction in DOF to be solved. It is worthwhile to remark that, since the eigenvalue problem solution is a subspace of the eigenvector problem, as reported in [44], the eigenvector amplitude and sign may depend on the solution algorithm adopted. Usually, in the solution process, eigenvectors are normalized with respect to mass, imposing for each m-th mode ∆u m that ∆u T m M∆u m = 1 (13) and then A principal characteristic of modal analysis is the spectral decomposition: the modes are orthogonal with respect to each other and they form an orthogonal basis in the modal coordinates (displacements) η. Thanks to this, system dynamic response can be approximated with the truncated summation of each mode response u = n ∑ m=1 ∆u m η m = ∆u η (15) in which η m is the modal coordinate relative to each retained mode, ∆u m is the modal shape and m is the number of retained modes. Due to the orthogonality of the modal basis, each mode solution is independent from the others and can be sees as a single DOF dynamic system: in this way, the mass and stiffness matrices become diagonal, and the system equations can be written as follows in whichη m and η m are, respectively, acceleration and displacement of the modal coordinate, and F m the load value for the m-th mode. In case the mass normalization is applied, M mm value is the unity.

Modal Superposition for Static Analysis
Modal superposition can be used in static structural analysis to approximate the solution by superposing natural modes, under the assumption of linear behavior system. According to this, Equation (16) can be written as The modal forces F m can be obtained as reported in next section, see Equation (21), considering loads related to pressure and shear stresses acting over the entire structure. By computing η m from Equation (17) it is possible to evaluate new structure configuration reverting from modal to natural coordinates thanks to Equation (15).

Unsteady Analysis Modal Approach
In case the FSI study cannot be represented using steady formulation reported in Section 2.4, it is possible to evaluate each modal coordinate time evolution using the following formula, as reported in [44] in which η 0 andη 0 are, respectively, the modal coordinate and the modal velocity at initial time. ζ is the damping factor and ω d are the circular frequencies for the analyzed system, which can be expressed as In Equation (18), the integral expression on the second sum term in the right part of the equation is generally referred to as the Duhamel integral. This integral states that, in a linear system, the reaction to a force f (t) can be evaluated as the sum of all differential responses generated during the loading history: so it is possible to evaluate the whole response to f (t) as the sum of responses to the single impulses that constitute the time evolution of f (t) itself. In numerical analysis, in order to apply this formulation, it is necessary to assume the force as constant during each simulation time step. In this case, Equation (19) is considered valid for each time step: η 0 andη 0 are used to evaluate modal coordinate at the end of a time step starting from their values at the beginning of the time step: Equation (20) can be successfully employed in the implementation of a time marching algorithm. Additionally, in the case of unsteady analyses, the structural deformation can be determined by summing contribution of a limited number of retained modes, weighted by time dependent modal coordinate, evaluated according to Equation (20). Natural coordinates can be obtained using Equation (15).

Setup of Modal FSI Analysis
By bringing together and integrating the RBF mesh morphing technique presented in Section 2.1 and the theory highlighted in Section 2.3, it is possible to generate an FSI approach suitable for solving a wide range of industrial problems. The starting point of this method is the natural modes extraction of the structure involved in the FSI problem, which will be used to parametrize the fluid domain. Modal shapes can be obtained not only from FEM calculation, as presented in this work, but also from experiments or numerical or theoretical approaches, as, for example, from the beam theory. In classical FSI approaches, the so called wetted surfaces of both structural and fluid-dynamics models, are accurately represented. Furthermore, since CFD and FEA grids generally have a different spacing due to different solver requirements, an interpolation of data exchanged between models is required. In the modal approach here proposed and presented, this level of detail needed for data interpolation between non-conformal meshes, is not required. Thanks to RBF interpolation, even a low level detail structural mesh can be used to extract natural modes that can be used to modify CFD mesh. This allows the use of simple beam FEA models even if the classical two-way approach cannot be used.
As anticipated, the modal approach allows to make parametric the CFD model on modal coordinates and the model shape is updated during the CFD calculation. Each retained natural mode shape is weighted using modal coordinates, which are related to the modal forces evaluated from pressure and friction forces on wetted surfaces. Since modal forces can be seen as an "interface" between structural and fluid-dynamics solvers, they play an important role, so the method used to evaluate them will be presented. In the proposed approach, modal loads are evaluated on the CFD model, integrating actions on nodes of the wetted surfaces and projecting them from the natural coordinates onto each modal shape. This projection of nodal forces from CFD onto the nodal displacement of the modal shape is summed for all nodes of wetted surfaces It is worthwhile to remark that, when integrating modal forces calculation in a parallel CFD solver, particular care to the way the are evaluated has to be paid. First of all, because pressure and shear loads are usually known at centroids, and because in a parallel calculation the cells on the boundaries of each partition could be summed twice.
In Figure 6, an overview of the methodology steps is given. Using the modal shapes, an RBF solution for each retained mode is generated and stored in the database of parametric mesh models. As stated before, for each structure it is necessary to store this database only once. Then, the mesh update is implemented within the CFD solver using RBF mesh morphing. With these tools, structure response to CFD loads can be evaluated in the modal space: mesh morpher is initialized with stored RBF solutions and the routines to evaluate modal forces, which are called during the CFD computation, provide the weights for each selected modal shape to evaluate the structure deformation. To ensure the coupling of the fluid and solid subproblem, a different convergence criterion for steady and unsteady FSI analyses is adopted. In steady FSI simulation the criterion is the deformed shape stability: after each FSI iteration the displacement of the most significant CFD mesh node is evaluated and compared with the previous iteration. Authors observed that after the first FSI iteration an over-elongation of about 20% is usually registered and that the steady FSI simulation converges with a tolerance lower than 10 −3 after 3-5 iterations. In case of unsteady FSI simulations, the coupling between solid and fluid problem is weak, being the CFD grid updated only after each time step and not at any inner CFD iteration, so the convergence criterion adopted is the time step bisection solution convergence: in the FSI setup phase an initial time step for the CFD solver is selected and the FSI solution evaluated. In case the time step is too big the simulation will not evolve because of the high mesh deformations, so the time step is divided in half: a convergent time step is reached when at the next time step bisection the FSI solution does not change. It must be noted that by using the Duhamel integral Equation (18) stability is assured, in contrast to the explicit coupling that requires a fine tuning of the time-step in order to provide stable fluid-structural coupling. With the mode embedding approach it is possible to solve different FSI problems, as, for example:

1.
Transient FSI in which structural displacements are known in advance and, therefore, imposed to the structure. It is the case of flapping devices experiencing complex motion, which was computed using FEA, multi-body, or observed in experimental tests. Acceleration of structural modes can be used to prepare ROM for non-linear flutter analysis; 2.
Steady FSI in which structural displacements are evaluated on the CFD mesh. It can be applied in aeronautical and motorsport applications, and can also be used to perform optimization considering coupled response. In some cases structural deflection can influence in a strong way the CFD mesh around wetted surfaces; 3.
Full coupled transient FSI in which the interaction between fluid flow and structure is captured by a time marching solution.
For the listed FSI problem, the update of modal coordinates and forces is performed with some differences. In the first case, structure time evolution is known in advance and imposed to it: with the mode superposition approach it is possible to obtain the desired deformation and to evaluate the modal forces relative to the structure motion. For the second and third FSI typology problems, the structural motion is the result of the CFD forces acting on the wetted surfaces. The implementation is more complex since the CFD solver has to handle the structural coordinates evolution; with the software tools adopted for cases described in this work, these functionalities are already developed and ready to be used. In case of steady FSI problem, the structural shape update is performed considering that, adopting mass normalization of structure natural modes, the modal coordinates can be expressed using Equation (17): Then, it is possible to obtain the parametric mesh formulation as follows The mesh update can be performed after a specific number of CFD solver iterations: a good practice in FSI runs involving wings is to update shapes less than ten times during the whole steady calculation. With transient FSI simulations it is necessary to store both modal coordinates and relative derivatives. The update of modal coordinates is performed according to Equation (20), using a modal load assumed to be constant within each time step, the shape update is performed at the end of each time step adopting Equation (11).
Regarding the generality of the proposed method, specially in application with steady FSI problems, some doubts can be expressed about the optimal number of natural modes to be retained (see [45]). It has been demonstrated that, for each new structure to be studied, a modal shapes qualification procedure can be put in place, if no consideration from previous experience can be of help [46]. The procedure requires, obviously, the generation of a FEA and CFD model of the structure. CFD model loads acting on wetted surfaces are extracted and used as input for the FEA model. The FEA model is used both for extracting a large enough number of modal shapes and to evaluate static deformation of the structure under CFD loads. The difference field between the FEM deformed and the modal amplified solution (composed by a reduced number of natural modes) represents the nodal truncation error that can be used to decide whether the selected modes can adequately approximate static structural deformation.
The approach presented here has, if compared with other ones reviewed in the introduction, both advantages and drawbacks. The modal approach requires only one structural calculation: the modes extraction to be performed once before CFD solution. Data extracted, namely the modal shapes with relative frequencies, can be used even if a change in the CFD boundary conditions occurs, being data which depend on the structure configuration. In this way the number of FEA calculations are minimized, if compared to the number of runs needed for the 2-way approach. Furthermore, integrating into the CFD model the capability to deform the structure via the parametric mesh paradigm, all the data transfers and interpolations needed in other FSI approaches are no longer necessary, reducing time to complete a FSI simulation. The adoption of a mesh morphing methodology apply an additional speed-up to the FSI simulation, since there is no need to perform remeshing of the CFD model and data mapping on FEA model surface mesh: this has a great impact in transient FSI analysis for which, if compared to the ANSYS system coupling method, a 10× speed-up can be achieved. The above described advantages come with a major limitation of the method: being based on linear modal theory, it can be applied only to linear problems. Efforts are currently being performed to overcome this limitation.

Setup of Non-Linear Imposed Wall Motion
In some applications there is the need to adapt the fluid dynamic domain to a prescribed non-linear motion of the wetted surfaces. In this field of application there is for example the case in which the non-linear behavior of a structural device is physically stimulated by a load known in advance, and that is partially influenced by the flow. In this case the interaction between fluid and structure can be significantly simplified, relying to a 1D coupling in which the structure is simulated in advance on CSM and then the obtained history of motion repeated on the CFD mesh by applying the displacement vector field to the deformable surfaces. To tackle this non-linear motion by means of a linear method, such as RBF, the complex field describing the displacements in time is partitioned in a given number of steps, that can be overimposed in series to replicate the non-linearity. Depending on how complex the displacement is, the number of steps can be increased. In Figure 7, an example of the described morphing strategy is shown, in which the non-linear opening and closing of an aortic valve is partitioned in a number of linear steps.
The nodal displacements inside each step are controlled by means of an amplification factor, that can vary from 0 to 1 according to a prescribed law. For a linear variation in time a law can be: (24) in which the step i, replicating the motion between times t i and t i+1 , has amplification A i (t).
Other laws of motion, such as sinusoidal or quadratic, can be also considered depending on the requirements.

Industrial Applications of FSI Modal Approach
In this section, a broad review of several applications tackled by making use of the methods previously described is given, pertaining to diverse fields of engineering. The problems, the strategies and the results are briefly presented for each application to give a picture of what it is possible to achieve by coupling an RBF-based method with the modal approach for three different classes of application: Transient FSI with prescribed displacement field, steady FSI, and unsteady FSI. Among the FSI cases here presented in the next section, two of them, the HIRENASD case [47] and the 3X4 tubes array in an inline cross-flow heat exchanger [48] are reference benchmarks in steady FSI simulations and in unsteady flow induced vibrations, respectively. Another interesting FSI-RBF benchmark not presented here, tackled by the authors and available in literature, is RIBES [49,50]. Being applications collected from literature, for each case the source is cited to allow readers to have a more detailed description of numerical setup (such as CFD solver, FEM solver, RBF setup) and CFD settings (such as flow conditions, turbulence models, FEM loads and constraints). The presented method was conceived to be general with respect to the solvers, meaning that any CFD and CSM softwares can be used. Thanks to the low level of distortion introduced in the computational grid, the same CFD settings prescribed by the state of the art best practices, and employed for a rigid simulation, can be adopted for the FSI problem. The same level of detail achieved by means of a traditional full detailed CFD simulation can be obtained through the FSI calculation.

Transient FSI with Prescribed Displacement Field Examples
In many cases, there are systems in which the displacement field of some boundaries are known in advance, and the analyst is interested in knowing how this prescribed time history can influence other structural elements and the flow. A notable example of this class of applications are, for instance, the flapping devices [51], in which a law of motion is applied to the wings but their deformation and interaction with the surrounding flow in unknown. In this case, the thrust is generated by means of the transient FSI interaction between the flexible wings and the flow. Another example of transient interaction with imposed displacement is given in [47], in which the authors excite in sequence each structural mode in a transient modal-based FSI problem, with the aim of obtaining a frequency response function (FRF) of the pressure with respect to displacements. A similar approach was pursued in [52] on an AGARD 445.6 wing using ANSYS Fluent as CFD solver and MSC. Nastran to solve the structural problem: a reduced order model (ROM) of the aeroelastic system around its steady flight condition was achieved by exciting by means of a quasi-Gaussian function each of the first four modes. These transient FSI problems were employed to generate the generalised aerodynamic force matrices (GAF). This method proved its reliability by reducing the errors in the transonic dip identification from 30%, typical values of low-fidelity computations based on doublet lattice methods, to 1%. In Figure 8, on the left, the flutter velocity results are shown for the AGARD 445.6 case, on the right the first four modes employed for the transient FSI simulation are displayed. Another interesting example of the transient FSI problem with prescribed displacements is shown in [29] for a biomedical application. With the aim of reducing the complexities related to a fully transient 2-way FSI problem, authors managed to reproduce a full systolic cycle of a polymeric prosthetic heart valve (P-PHV) by means of a 1-way coupling, taking into account structural non-linearities and inertia effects. At first the structural system was studied by applying a physiological time-varying pressure uniformly applied to the ventricular portion of the valve, in order to simulate the opening pressure. Nodal displacements were then extracted from the CSM solver ANSYS Mechanical and imported in ANSYS Fluent by means of the RBF-based morpher RBF Morph. In Figure 9, on the right, the blood velocity contour obtained during four systole instants of time are shown, together with the corresponding valvular positions and the flow lines. In Figure 9, on the left, the comparison with a workflow based on remeshing is shown, demonstrating the quality of the approach, with a speed-up of the procedure quantifiable in 16×.

Steady FSI Examples
As a demonstration of the static modal approach for FSI using RBF morphing, the steady aeroelasticity of the Piaggio Aerospace P1XX wind tunnel model was studied in [46] using NX Nastran to compute the eigenvalue problem and ANSYS Fluent to solve the CFD. The behavior of the business jet class airplane was explored in a transonic range, comparing obtained results with those achieved by means of a traditional 2-way approach and with experimental data. Only the wing was considered flexible, and a convergence study with respect to embedded modes was carried out, demonstrating that a good result can be obtained by using only the first two eigenvalues. In Figure 10, left, the steady aeroelastic configuration of the model is shown, with deformations amplified by a factor 10 to highlight wing flexibility. In the picture on the right, the pressure contours computed on the wing suction side obtained by means of the modal superposition approach and the 2-way method are compared to those observed experimentally. The maximum difference between the 2-way and the modal approach in terms of displacements is in the order of 0.03 mm, equal to 3.1%. Another interesting problem, tackled by means of a steady modal superposition approach was conducted on the HIRENASD test case [47] developed for the NASA Aeroelastic Prediction Workshop (AePW). In this case the focus of the study was to demonstrate that the modal superposition method is suitable to be used for problems in which the structural and fluid dynamic models have a discrepancy. In Figure 11, on the right, the structural model of the HIRENASD model is shown: the wing wetted surface is attached to a rig specially conceived to tune numerical and experimental data that does not match the CFD domain shown on the left. Wing root, on the structural problem, is also left free to vertically translate and so this movement must be taken into account also for the FSI problem. A special RBF setup was put in place, to apply the modal shapes extracted from the eigenvalue problem on the wetted surface, fixing all the fuselage but a buffer area supplying a deformable region to accommodate the vertical displacement of the wing root. On Figure 11, left, the baseline geometry is shown together with the deformed configuration pertaining to the second mode, the wing root vertical displacement is clearly visible. Results of the simulation, achieved using ANSYS Fluent and NX Nastran were compared to those achieved by means of a traditional 2-way method and experimental testing in terms of tip displacement, increasing the number of modes employed for the modal FSI approach. In Figure 12, results are shown demonstrating that no more than 4 modal shapes are required to achieve a satisfactory result. A similar field of application is given by cases in which the level of detail of the CFD and CSM models is very different. In this case the meshless property of the RBF method allows to tackle the problem, applying structural modes achieved for instance on a geometry modeled using plate elements on a solid 3D CFD case, as shown in Figure 13 left. With the same strategy an efficient shape optimization can be carried on FSI cases in which the design variation influences the flow but has a limited impact on the structural behavior. As an example of this class of problems, the optimization carried within the EU FORTISSIMO 2 project [54] on the Piaggio Aerospace Avanti EVO aircraft is here shown, taking advantage of the RBF4AERO optimization platform [55]. Thanks to the multi-solver feature of the method proposed, in this case the CFD simulation was carried on SU2, while the structural problem was solved in Nastran. Three shape parameters on the winglet modifying the cant, sweep, and twist angles were conceived to improve the efficiency of the aircraft, as shown in Figure 13, right. For each design point (DP) the shape parameters were at first applied according to the value prescribed by the design of experiment (DOE) table, and then the baseline modal shapes-extracted from CSM only once prior to the optimization-were used to carry the FSI simulation. Although the structural response of the wing was not changed from DP to DP, in this way the final results depend strictly from the forces applied to the winglet, which vary according to the shape variation. The results of the optimisation are shown in Figure 14 on the right. The chart shows that there are candidate points able to provide for improved efficiency, optimizing the winglet on cruise conditions. On the left side, the deformed configuration for baseline is shown together with the rigid original model. A very similar approach was shown in [56], using OpenFOAM and ANSYS APDL for the optimization of the WATTsUP electric aircraft propeller, with an increase in the weighted sum of efficiencies at two different flight regimes equal to 4.0%.

Transient FSI Examples
The RBF mesh morphing based modal superposition method main characteristics, i.e., the low computational requirements and the high flexibility, can be a big advantage in dealing with transient FSI studies. In [57], authors made a successful attempt in reproducing the test case illustrated in [58]: a NACA0009 hydrofoil immersed in water. The main scope of the study was to numerically reproduce the vortex induced vibrations and the lock-in/lock-off frequency according to data available in literature. The whole study was performed in ANSYS Fluent and Mechanical, adapting the boundary conditions reproduced in Figure 15 and varying the flow velocity. Numerical simulations were performed using the first six natural modes of the NACA profile, then the fluid velocities simulated were set at 16 m/s and 22 m/s. Displacements of one point near the trailing edge was monitored and studied in the frequency domain. Obtained results, reported in Figure 15, right, fit very well with available experimental data, reporting, respectively, 909.91 H z and 1202.0 Hz values for lock-in and lock-off frequencies. Analyzing FSI results, the hydrofoil encounter resonance effects at 909.91 Hz when the flow velocity is 16 m/s, which is lower than the first natural frequency obtained extracting structure eigenvalues (1133.8 Hz). To investigate the structure natural frequencies when immersed in water, an additional FSI simulation was performed imposing a zero velocity to the fluid and amplifying each modal shape at initial simulation time. With the same settings a simulation in air was also performed. To evaluate the structural natural modes in water and air, a spectral analysis was performed on displacement converted from the time domain into the frequency domain: results are shown in Figure 16. Data reported in Figure 16, left, are obtained in air and match very well with results from eigenvalue analysis. Simulation in water results (depicted in Figure 16, right) highlights that the first natural frequency is 891.9 Hz, a value that matches with experimental data.
The modal superposition method allows to easily impose a prescribed load or deformation, as for the immersed hydrofoil, by using the modal coordinates or using the modal forces of Equation (21), as for the study reported in [59], in which at each FSI cycle the modal force generated by the store weight was first added to the modal forces evaluated on wetted surfaces and then removed to simulate the store separation. Transient modal superposition can also be successfully applied when more than one structure has to be studied in an FSI analysis. In [60], a 3X4 tubes array in an inline crossflow heat exchanger was analyzed in ANSYS Workbench using Fluent and Mechanical, and results compared with the ones found in [48]. In this case, since the tubes were identical, only one eigenvalue analysis was performed to extract modal shapes and natural frequencies for all the tubes. To apply the RBF mesh morphing only a simple translation was needed. In Figure 17, a section of the computational domain and a view of the first three modal shapes is shown.

Advanced Approaches for Limiting Mesh Distortion
One of the advantages in applying shape modification via RBF based mesh morphing is the resultant high quality of the distorted mesh at the end of the process, if compared with different methods. Unfortunately, this advantage is highly reduced if between CFD surfaces on which displacement are imposed only few rows of volume elements are placed. This limits the 'buffer' of elements on which the displacement can be smoothed, rapidly decreasing the morphed volume mesh quality.
This can limit the application of the modal approach in all cases in which small clearances and gaps are present, especially if large deformations occur. An example of this problem can be found in [61], where a thermowell FSI study was tackled coupling ANSYS Fluent with ANSYS Mechanical by means of RBF Morph. In such applications, the end face of the thermal sensor protection is very near the tube surface and applying mesh morphing will result in highly distorted cell, as depicted in Figure 19. To overcome this limit, a combination of RBF shape modifications can be successfully applied. The procedure foresees that nodes in the shadow area of the thermowell (node on the tube surface identified by projecting the thermowell tip on it) have to follow tip by rotating around and moving along the tube axis. Due to the fact that the RBF solutions are built on the modal shapes of the thermowell only, an additional shape modification is needed to project back morphed nodes on the tube surface (see Figure 20). Thanks to this mesh morphing procedure, the results obtained allowed to successfully perform transient FSI analysis, successfully catching the vortex induced vibrations phenomenon.

Conclusions
In this work an overview of the possible problems that can be tackled by relying to a proper combination of embedded vector fields using RBF morphing was given. CFD and CSM models can be coupled by means of the modal embedding performed on the fluid dynamic environment, using RBF based mesh morphing to make the CFD model parametric with respect to the modal coordinates. Several application examples were illustrated, considering steady and transient FSI problems using structural modes superposition, for the case of linear structures, or imposing a full non-linear displacement time history. In all these applications, the proposed method allowed to obtain reliable results with limited computational results if compared with classical FSI approaches. Nevertheless, the proposed method has limitations: being based on mesh morphing, it is possible that after deformation, the mesh can result too distorted to be employed; furthermore, by relying on modes extracted by means of an eigenvalue problem, the field of application is limited to linear deformation problems. To overcome the first, a strategy to reduce mesh deformation was proposed in the paper and studied in literature. For the latter, a research based on the non-linear correction of modal shapes is under way, to broaden the field of application of the method also for cases outside the linear range of deformation.
An additional consideration can be performed for unsteady FSI analyses. The procedure described in this work employs a weak coupling between the solid and fluid subproblems, since the CFD grid update is performed at each time-step. An improvement of this procedure is currently being investigated: thanks to the low cost of the structural solution (which is evaluated once and integrated by means of modal superposition inside the CFD solver) the CFD mesh update can be computed at each internal iteration at any given time-step, generating a strong coupling.
Funding: Researches and activities presented in this paper received funding from the following EU research projects:

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Symbols
The following symbols are used in this manuscript: Vector of structural displacements U Interpolation matrix used for RBF fit X CFD Mesh nodes position X CFD 0 Mesh nodes in the undeformed position

Definitions, Acronyms and Abbreviations
The following abbreviations are used in this manuscript: