# B-Spline Impulse Response Functions of Rigid Bodies for Fluid-Structure Interaction Analysis

Authors: S. Zhou-Bowers; D. C. Rizos
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Reduced 3D dynamic fluid-structure interaction (FSI) models are proposed in this paper based on a direct time-domain B-spline boundary element method (BEM). These models are used to simulate the motion of rigid bodies in infinite or semi-infinite fluid media in real, or near real, time. B-spline impulse response function (BIRF) techniques are used within the BEM framework to compute the response of the hydrodynamic system to transient forces. Higher-order spatial and temporal discretization is used in developing the kinematic FSI model of rigid bodies and computing its BIRFs. Hydrodynamic effects on the massless rigid body generated by an arbitrary transient acceleration of the body are computed by a mere superposition of BIRFs. Finally, the dynamic models of rigid bodies including inertia effects are generated by introducing the kinematic interaction model to the governing equation of motion and solve for the response in a time-marching scheme. Verification examples are presented and demonstrate the stability, accuracy, and efficiency of the proposed technique.


---

## Body

## 1. Introduction

Fluid-structure interaction (FSI) is a highly complicated phenomenon which involves wave propagation in fluid domains, structural dynamics, and the dynamic interaction between these two coupled media. The response of the FSI system pertains to computing the dynamic response of the structure along with scattering and radiation of waves in the fluid. The latter is caused by the simultaneous action of the motion of the structure and propagating waves in the fluid domain. FSI problems are inherently two-phase problems, i.e., the fluid phase and the solid phase. Analytical solutions are generally restricted to extremely simple geometries and adopt oversimplifying assumptions for FSI analysis. Model-based simulations of the coupled multiphase problem use computational methods, such as the finite element method (FEM), the boundary element method (BEM), and coupled techniques [1–11]. Typically, separate models are developed for the fluid and solid phases. The FEM is used in general to model the structure and/or the fluid. However, modeling of infinite media using the FEM requires the use of special elements at the boundary of the truncated domain. The BEM is typically used to model the wave propagation in infinite fluid media because of its established advantages in dealing with systems of infinite or semi-infinite extents [11]. Direct time-domain BEM formulations of the external acoustic problem, like the one employed in this work, have an additional advantage since they eliminate any spurious modes appearing in frequency-domain formulation that require special techniques, like the one reported in [12]. However, frequency-domain formulations may employ the mirror-image technique to eliminate the need for discretization of the free field of the half space if present [13]. Solutions are obtained in a direct, “monolithic” fashion, where the algebraic equations are solved simultaneously [14, 15] or in a staggered manner [16, 17], by satisfying, explicitly or implicitly, equilibrium and compatibility conditions along the interface between the fluid domain and structure to account for FSI. Challenges arise from the typically large differences in physical characteristics of the two phases and the usually large deformations in one phase that commonly result to spatial and temporal discretization mismatch, as well as degradation of the solution over time [18]. These challenges are typically addressed at a cost of computational efficiency. To increase computational efficiency, reduced-order models (ROMs) have been developed within the FSI framework [19] using either eigenvalue analysis [20, 21] or transfer functions [22]. Detailed accounts and critical discussion of the various FEM, BEM, and FEM-BEM coupled formulations and solution techniques along with an extensive literature review of original works can be found in the state-of-the-art review chapter on FSI reported in [1].This work presents the development of a 3D hydrodynamic fluid-structure interaction mathematical model between the linearized fluid and rigid structure by coupling the BEM with the equation of motion of the rigid body that is expressed in an FEM sense. In general, the motion of the structure depends, among other forces, on the hydrodynamic forces developed by the interaction between the structure and the surrounding fluid due to the body’s motion. In this work, a B-spline impulse response function (BIRF) technique is used to compute the associated hydrodynamic forces. To this end, the BIRF of the fluid boundary, in the absence of any structures, is calculated first based on the time-domain boundary element method for wave propagation in fluid media published by the authors in [23]. In this work, a rigid body boundary element compatible to the aforementioned BEM is developed and used to compute the BIRFs of the rigid structure in contact with the fluid. The development of the rigid boundary element formulation for fluid-rigid body interaction and the calculation of the corresponding BIRFs is based on the equivalent work of the authors for soil-rigid body interaction analysis reported in [24–26]. The proposed kinematic model computes the hydrodynamic forces and moments at a reference point of the rigid body due to arbitrary acceleration applied at the same point based on a mere superposition of the BIRFs of the rigid structure. Finally, the dynamic equations of motion of the rigid body, as modeled by the FEM, are coupled in the direct time domain with the proposed kinematic interaction model to develop a complete FSI model for real-time simulation of the motion of massive rigid bodies floating on or submerged in the fluid. The following sections present a brief discussion of the BEM, detailed accounts of the rigid body boundary element, the proposed kinematic interaction models, their coupling with the general equation of motion, and the time-marching scheme of the solution. Verification of the proposed methods is presented through comparisons with analytical and FEM solutions of a representative system as obtained by commercial software. The stability, efficiency, and accuracy of the proposed methodology are demonstrated.


## 2. B-Spline BEM for 3D Scalar Wave Propagation

A time-domain BEM formulation for wave propagation through an inviscid, compressible acoustic medium was proposed by the authors [23] and briefly presented here for completeness. The governing equation of 3D scalar waves traveling in a linear, isotropic, inviscid, and homogeneous fluid is Lamb’s equation expressed as follows:(1)∇2px,t−1c2p¨x,t+fx,t=0,where p is the pressure field, x is the position vector, t represents time, f is the body source, ∇2 is the Laplacian operator, and c is the wave propagation velocity in the medium. Following well-established procedures [5], the governing equation can be expressed in a boundary integral equation (BIE) form as follows:(2)cξpξ,t=∫SGx,t,ξqx,t−Tx,t,ξpx,tdS,where x and ξ are the points in the domain or its boundary, the integration is over the boundary, S, of the domain, and G and T are the associated B-spline fundamental solutions, as reported in [23]. The field variables px,t and qx,t=∂px,t/∂n correspond to the pressure and its derivative in the direction of the outward normal, n, to the boundary, respectively. The term cξ is known as the “jump term” and depends on the location of the point ξ and the smoothness of the boundary in its neighborhood.Following the procedures in [23], the BIE (Equation (2)) can be written in a discrete form in space and time and expressed in a matrix form at time step N as follows:(3)12pN=∑n=1N+1GNqN−n+2−TNpN−n+2.It is assumed that the bounding surfaceS is discretized into NE number of boundary elements with NN number of boundary nodes based on a higher-order spatial discretization. Each node is assigned one degree of freedom (DOF). Vectors p and q are of size NN and contain nodal values of the pressure and its normal derivative, respectively. The coefficient matrices GN and TN are of size NN × NN and represent the influence of a boundary node on another boundary node, a typical element of which is a scalar value that is computed based on surface integrations of the fundamental solutions. Superscripts indicate the time step at which quantities are evaluated. The time step, Δt, depends on the duration (support), Dt, of the B-spline polynomial used in the derivation of the fundamental solutions and is equal to Δt = Dt/4 for the cubic B-splines adopted in this work. The BIRF of the boundary of the fluid domain due to a single load case pertains to the pressure at all NN nodes due to a unit B-spline impulse flux excitation of duration Dt applied at a single node, j, and is computed by Equation (3). Subsequently, the B-spline impulse response vectors can be computed for all excited nodes j = 1,2,…, NN and collected in a matrix form as follows:(4)BN=b1N,b2N,…,bjN,…,bNNN.Equation (4) is the complete BIRF of the boundary of the fluid domain expressed in a discrete form in space and time. Each element bij in the matrix represents the impulse response pressure of node i at time step N due to a “unit” B-spline excitation flux, applied at node j. This matrix is a characteristic of the system and needs to be computed only once for the specific geometry of the free surface of the solution domain. In the general case, the BIRF is computed for all nodes of the problem, and the impulse response matrix, BN, is square. In most practical cases, however, the BIRF of the system needs to be computed only for those nodes that are expected to carry a forced excitation at any time during the solution, yielding a rectangular matrix or a square matrix of significantly reduced size. Because of the small, finite duration of the B-spline impulse excitation and the wave attenuation, the number of time steps N that the BIRF matrices are obtained for is small compared to the response of the system to general excitations and limited to only a few time steps [23]. Once this characteristic response of the system is known, the transient scattered pressure field, p, at step N due to an arbitrary incident flux, q, can be computed as a mere superposition of the BIRFs as follows:(5)pN=∑n=1N+1BnqN−n+2=FqN+HN,where(6)F=2B1+B2and(7)HN=∑n=3N+1BnqN−n+2−B1qN−1.The proposed approach is very efficient especially when multiple load cases of prolonged duration are considered, since the BIRFs are independent of the external excitation and are typically of much shorter duration than the external excitation.


## 3. Kinematic Fluid-Structure Interaction Analysis of Rigid Bodies

This section discusses the development of the proposed kinematic FSI model for rigid bodies. This model is developed on the basis of BIRFs as they pertain to rigid bodies in fluid media. The computed BIRFs provide the relationship between applied acceleration and resulting hydrodynamic forces and account for kinematic fluid-structure interaction effects only. They will be incorporated in the dynamic models presented in Section4 at which time, inertia interaction effects will be accounted for. The motion of a rigid body can be represented by the motion of a reference point (RP). The RP can be anywhere in the solution domain and does not have to be physically located in the rigid body or on its surface boundary. For simplicity and without loss of generality, the RP is located at the center of gravity of the body, as depicted in Figure 1. The motion of the reference point is described by three translational and three rotational DOFs, identified in Figure 1. In this work, the driving excitation is applied at the reference point and corresponds to the translational and rotational accelerations exerted on the massless rigid body. Provided that the BIRFs of the wet surface of the rigid body (in the absence of the latter) have been calculated based on the procedures discussed in Section 2, the derivation of the BIRFs of the rigid body is discussed in the following.Figure 1
The 6 degrees of freedom of the motion of a rigid body and the definition of the reference point RP.


### 3.1. Treatment of Rigid Bodies

#### 3.1.1. Rigid Surface Boundary Element

The rigid surface boundary element is derived within the framework of the B-spline BEM for fluids and follows an approach similar to the one developed by the author for treatment of rigid bodies interacting with soils and solids [24]. This element is located at the wet surface of the rigid structure and is defined geometrically by 2 nodes. The first node, WS, is located on the wet surface of the rigid body, and the second node, RP, is located at the reference point of the rigid structure as indicated in Figure 2.Figure 2
Definition of the rigid surface element.The node on the wet surface is assigned three DOFs,u¨WS, that pertain to the translational components of the acceleration at this point. The reference point node has six acceleration DOFs, aRP=u¨RPθ¨RPT, pertaining to the three translational, u¨RP=u¨RPxu¨RPyu¨RPzT, and the three rotational, θ¨RP=θ¨RPxθ¨RPyθ¨RPzT, components. Since the body is assumed to be rigid, its deformation is ignored. Therefore, the acceleration of the reference point can be directly related to the acceleration of the wet surface node of the element through geometric rigid body conditions as follows:(8)u¨WS3×1=I3×3u¨RP3×1+H3×3θ¨RP3×1,where I is a 3 × 3 identity matrix and H is expressed as follows:(9)H=0Δz−Δy−Δz0ΔxΔy−Δx0,in which Δx = xRP − xWS, Δy = yRP − yWS, and Δz = zRP − zWS are the distance components between the reference point, RP, and the point on the wet surface, WS, in the x, y, and z directions, respectively. For every point WS on the wet surface of the rigid body, there is a corresponding point F on the fluid surface, as shown in Figure 2, whose motion is described by the acceleration in the direction of the outward normal to the fluid boundary u¨Fn. Under the assumption that the rigid body remains always in contact with the fluid in the normal direction and there is no friction between the rigid body and the fluid, compatibility of the normal component of the acceleration requires that(10)u¨Fn=−u¨WSn=−u¨WS⋅n,while the tangential components are not coupled. Therefore, introducing Equation (8) into Equation (10), the normal acceleration, u¨Fin, of each fluid point Fi on the wet surface is related to the acceleration, aRP=u¨RPTθ¨RPTT, applied at the reference point of the rigid body as follows:(11)u¨Fin=−I3×3u¨RP3×1+H3×3θ¨RP3×1⋅ni,where ni represents the outward normal vector at point i on the wet surface of the rigid body. Subsequently, the normal acceleration of all MM points on the wet surface is expressed in the vector form by collecting the nodal accelerations as follows:(12)aF=u¨F1nu¨F2n⋯u¨FMMnT.At the interface between the rigid body and the fluid media, flux field at fluid nodes on the wet surface is related to the acceleration as follows [27]:(13)∂p∂nF=qF=−ρaF,where qF=∂p/∂nF is the vector of the nodal flux normal to the wet surface and ρ is the mass density of the fluid medium. In view of Equations (7) and (13), the pressure distribution on the wet surface nodes can be computed by superposition of appropriate BIRFs as follows:(14)pFN=∑n=1N+1BN−n+2qFN=−ρ∑n=1N+1BN−n+2aFN,where superscripts indicate the time step, pF is the pressure at each node on the fluid surface caused by the applied acceleration, ρ is the density of the acoustic medium, and B is the B-spline impulse response matrix of the fluid surface, as determined in Section 2. Due to the previous assumption that the rigid body remains always in contact with the fluid, the pressure distribution along the wet surface of the rigid body can be computed as follows:(15)pWSN=−pFN=ρ∑n=1N+1BN−n+2aFN.In view of Equations (8)–(15), a rigid body boundary element is established that relates the hydrodynamic pressure that develops on the wet surface of the rigid body to the acceleration applied at the reference point of the rigid body.


#### 3.1.2. B-Spline Impulse Response of the Rigid Body

The BIRF of a rigid body pertains to the hydrodynamic forces and moments that develop at the reference point due to a B-spline acceleration applied at the reference point and can be computed as follows. A unit B-spline acceleration is applied in the direction of each degree of freedom of the reference point, and Equation (15) computes the hydrodynamic pressure on all nodes of the wet surface. However, the differential hydrodynamic force dFn in the direction normal to the wet surface at any point is computed based on the continuous pressure field pWS at this point as follows:(16)dFWSn=pWSdA,where dA is the differential area. Based on the field variable interpolation, the element discretization, and the directional cosines n=lmn of the outward normal vector at any point on an element, the Fx, Fy, and Fz components of the differential force dF vector are computed as follows:(17)dFWS=nTNspWSsdA,s=1,2,…,8.Subsequently, the differential force,dFRP, and corresponding moment, dMRP, at the reference point of the rigid body are computed as follows:(18a)dFRP=dFWS=nTNspWSsdA,(18b)dMRP=HdFWS=HnTNspWSsdA.Therefore, the contribution of each element, el, to the forces,FRPel, and moments, MRPel, at the reference point is determined by integrating Equation (18) over the area of the element, i.e.,(19a)FRPel=∫SeldFRP=∫SelnTNspWSsdA,(19b)MRPel=∫SeldMRP=∫SelHnTNspWSsdA.Finally, the resultant hydrodynamic forces,FRP, and moments, MRP, at the reference point of the rigid body are computed as a mere superposition of the element contributions computed in Equation (19) as follows:(20a)FRP=∑el=1NEFRPel=∑el=1NE∫SeldFEP=∑el=1NE∫SelnTNspWSsdA,(20b)MRP=∑el=1NEMRPel=∑el=1NE∫SeldMRP=∑el=1NE∫SelHnTNspWSsdA.To compute BIRFs of the rigid body, all six degrees of freedom of a rigid body need to be considered. The BIRF,BRP, of the rigid body is a 6 × 6 matrix. In order to obtain this matrix, unit acceleration with B-spline time modulation is applied in the direction of each of the six degrees of freedom, one at a time, and the corresponding forces and moments at the reference point of the rigid body are computed. This response represents one column in the BIRF matrix. Collecting all responses in a matrix form yields the BIRF matrix of the rigid body as follows:(21)BRP=FRPu¨xFRPu¨yFRPu¨zFRPθ¨xFRPθ¨yFRPθ¨zMRPu¨xMRPu¨yMRPu¨zMRPθ¨xMRPθ¨yMRPθ¨z,where the superscript indicates the acceleration component that generated the hydrodynamic force at the reference point. It should be noted that the reference point can be located in a way such that the coupling effects between the different modes may be ignored. Consequently, the BIRF matrix is a diagonal matrix, and the six modes of motion are completely uncoupled.


### 3.2. Kinematic Fluid-Structure Interaction Model

The last step in the development of the kinematic FSI model computes hydrodynamic forces and moments due to arbitrary acceleration. Once the BIRF of the reference point of the rigid body is evaluated as discussed in Section3.1.2, the equivalent hydrodynamic forces and moments, Fhyd, at the reference point of the system to an arbitrary acceleration applied at the reference point are computed at step N as a mere superposition of the BIRF responses as follows:(22)FhydN=∑n=1N+1BRPN−n+2aRPN.Therefore, kinematic FSI models can be obtained corresponding to the six primary modes of motion of the rigid body, as well as the associated coupled modes, if required.


## 3.1. Treatment of Rigid Bodies

### 3.1.1. Rigid Surface Boundary Element

The rigid surface boundary element is derived within the framework of the B-spline BEM for fluids and follows an approach similar to the one developed by the author for treatment of rigid bodies interacting with soils and solids [24]. This element is located at the wet surface of the rigid structure and is defined geometrically by 2 nodes. The first node, WS, is located on the wet surface of the rigid body, and the second node, RP, is located at the reference point of the rigid structure as indicated in Figure 2.Figure 2
Definition of the rigid surface element.The node on the wet surface is assigned three DOFs,u¨WS, that pertain to the translational components of the acceleration at this point. The reference point node has six acceleration DOFs, aRP=u¨RPθ¨RPT, pertaining to the three translational, u¨RP=u¨RPxu¨RPyu¨RPzT, and the three rotational, θ¨RP=θ¨RPxθ¨RPyθ¨RPzT, components. Since the body is assumed to be rigid, its deformation is ignored. Therefore, the acceleration of the reference point can be directly related to the acceleration of the wet surface node of the element through geometric rigid body conditions as follows:(8)u¨WS3×1=I3×3u¨RP3×1+H3×3θ¨RP3×1,where I is a 3 × 3 identity matrix and H is expressed as follows:(9)H=0Δz−Δy−Δz0ΔxΔy−Δx0,in which Δx = xRP − xWS, Δy = yRP − yWS, and Δz = zRP − zWS are the distance components between the reference point, RP, and the point on the wet surface, WS, in the x, y, and z directions, respectively. For every point WS on the wet surface of the rigid body, there is a corresponding point F on the fluid surface, as shown in Figure 2, whose motion is described by the acceleration in the direction of the outward normal to the fluid boundary u¨Fn. Under the assumption that the rigid body remains always in contact with the fluid in the normal direction and there is no friction between the rigid body and the fluid, compatibility of the normal component of the acceleration requires that(10)u¨Fn=−u¨WSn=−u¨WS⋅n,while the tangential components are not coupled. Therefore, introducing Equation (8) into Equation (10), the normal acceleration, u¨Fin, of each fluid point Fi on the wet surface is related to the acceleration, aRP=u¨RPTθ¨RPTT, applied at the reference point of the rigid body as follows:(11)u¨Fin=−I3×3u¨RP3×1+H3×3θ¨RP3×1⋅ni,where ni represents the outward normal vector at point i on the wet surface of the rigid body. Subsequently, the normal acceleration of all MM points on the wet surface is expressed in the vector form by collecting the nodal accelerations as follows:(12)aF=u¨F1nu¨F2n⋯u¨FMMnT.At the interface between the rigid body and the fluid media, flux field at fluid nodes on the wet surface is related to the acceleration as follows [27]:(13)∂p∂nF=qF=−ρaF,where qF=∂p/∂nF is the vector of the nodal flux normal to the wet surface and ρ is the mass density of the fluid medium. In view of Equations (7) and (13), the pressure distribution on the wet surface nodes can be computed by superposition of appropriate BIRFs as follows:(14)pFN=∑n=1N+1BN−n+2qFN=−ρ∑n=1N+1BN−n+2aFN,where superscripts indicate the time step, pF is the pressure at each node on the fluid surface caused by the applied acceleration, ρ is the density of the acoustic medium, and B is the B-spline impulse response matrix of the fluid surface, as determined in Section 2. Due to the previous assumption that the rigid body remains always in contact with the fluid, the pressure distribution along the wet surface of the rigid body can be computed as follows:(15)pWSN=−pFN=ρ∑n=1N+1BN−n+2aFN.In view of Equations (8)–(15), a rigid body boundary element is established that relates the hydrodynamic pressure that develops on the wet surface of the rigid body to the acceleration applied at the reference point of the rigid body.


### 3.1.2. B-Spline Impulse Response of the Rigid Body

The BIRF of a rigid body pertains to the hydrodynamic forces and moments that develop at the reference point due to a B-spline acceleration applied at the reference point and can be computed as follows. A unit B-spline acceleration is applied in the direction of each degree of freedom of the reference point, and Equation (15) computes the hydrodynamic pressure on all nodes of the wet surface. However, the differential hydrodynamic force dFn in the direction normal to the wet surface at any point is computed based on the continuous pressure field pWS at this point as follows:(16)dFWSn=pWSdA,where dA is the differential area. Based on the field variable interpolation, the element discretization, and the directional cosines n=lmn of the outward normal vector at any point on an element, the Fx, Fy, and Fz components of the differential force dF vector are computed as follows:(17)dFWS=nTNspWSsdA,s=1,2,…,8.Subsequently, the differential force,dFRP, and corresponding moment, dMRP, at the reference point of the rigid body are computed as follows:(18a)dFRP=dFWS=nTNspWSsdA,(18b)dMRP=HdFWS=HnTNspWSsdA.Therefore, the contribution of each element, el, to the forces,FRPel, and moments, MRPel, at the reference point is determined by integrating Equation (18) over the area of the element, i.e.,(19a)FRPel=∫SeldFRP=∫SelnTNspWSsdA,(19b)MRPel=∫SeldMRP=∫SelHnTNspWSsdA.Finally, the resultant hydrodynamic forces,FRP, and moments, MRP, at the reference point of the rigid body are computed as a mere superposition of the element contributions computed in Equation (19) as follows:(20a)FRP=∑el=1NEFRPel=∑el=1NE∫SeldFEP=∑el=1NE∫SelnTNspWSsdA,(20b)MRP=∑el=1NEMRPel=∑el=1NE∫SeldMRP=∑el=1NE∫SelHnTNspWSsdA.To compute BIRFs of the rigid body, all six degrees of freedom of a rigid body need to be considered. The BIRF,BRP, of the rigid body is a 6 × 6 matrix. In order to obtain this matrix, unit acceleration with B-spline time modulation is applied in the direction of each of the six degrees of freedom, one at a time, and the corresponding forces and moments at the reference point of the rigid body are computed. This response represents one column in the BIRF matrix. Collecting all responses in a matrix form yields the BIRF matrix of the rigid body as follows:(21)BRP=FRPu¨xFRPu¨yFRPu¨zFRPθ¨xFRPθ¨yFRPθ¨zMRPu¨xMRPu¨yMRPu¨zMRPθ¨xMRPθ¨yMRPθ¨z,where the superscript indicates the acceleration component that generated the hydrodynamic force at the reference point. It should be noted that the reference point can be located in a way such that the coupling effects between the different modes may be ignored. Consequently, the BIRF matrix is a diagonal matrix, and the six modes of motion are completely uncoupled.


## 3.1.1. Rigid Surface Boundary Element

The rigid surface boundary element is derived within the framework of the B-spline BEM for fluids and follows an approach similar to the one developed by the author for treatment of rigid bodies interacting with soils and solids [24]. This element is located at the wet surface of the rigid structure and is defined geometrically by 2 nodes. The first node, WS, is located on the wet surface of the rigid body, and the second node, RP, is located at the reference point of the rigid structure as indicated in Figure 2.Figure 2
Definition of the rigid surface element.The node on the wet surface is assigned three DOFs,u¨WS, that pertain to the translational components of the acceleration at this point. The reference point node has six acceleration DOFs, aRP=u¨RPθ¨RPT, pertaining to the three translational, u¨RP=u¨RPxu¨RPyu¨RPzT, and the three rotational, θ¨RP=θ¨RPxθ¨RPyθ¨RPzT, components. Since the body is assumed to be rigid, its deformation is ignored. Therefore, the acceleration of the reference point can be directly related to the acceleration of the wet surface node of the element through geometric rigid body conditions as follows:(8)u¨WS3×1=I3×3u¨RP3×1+H3×3θ¨RP3×1,where I is a 3 × 3 identity matrix and H is expressed as follows:(9)H=0Δz−Δy−Δz0ΔxΔy−Δx0,in which Δx = xRP − xWS, Δy = yRP − yWS, and Δz = zRP − zWS are the distance components between the reference point, RP, and the point on the wet surface, WS, in the x, y, and z directions, respectively. For every point WS on the wet surface of the rigid body, there is a corresponding point F on the fluid surface, as shown in Figure 2, whose motion is described by the acceleration in the direction of the outward normal to the fluid boundary u¨Fn. Under the assumption that the rigid body remains always in contact with the fluid in the normal direction and there is no friction between the rigid body and the fluid, compatibility of the normal component of the acceleration requires that(10)u¨Fn=−u¨WSn=−u¨WS⋅n,while the tangential components are not coupled. Therefore, introducing Equation (8) into Equation (10), the normal acceleration, u¨Fin, of each fluid point Fi on the wet surface is related to the acceleration, aRP=u¨RPTθ¨RPTT, applied at the reference point of the rigid body as follows:(11)u¨Fin=−I3×3u¨RP3×1+H3×3θ¨RP3×1⋅ni,where ni represents the outward normal vector at point i on the wet surface of the rigid body. Subsequently, the normal acceleration of all MM points on the wet surface is expressed in the vector form by collecting the nodal accelerations as follows:(12)aF=u¨F1nu¨F2n⋯u¨FMMnT.At the interface between the rigid body and the fluid media, flux field at fluid nodes on the wet surface is related to the acceleration as follows [27]:(13)∂p∂nF=qF=−ρaF,where qF=∂p/∂nF is the vector of the nodal flux normal to the wet surface and ρ is the mass density of the fluid medium. In view of Equations (7) and (13), the pressure distribution on the wet surface nodes can be computed by superposition of appropriate BIRFs as follows:(14)pFN=∑n=1N+1BN−n+2qFN=−ρ∑n=1N+1BN−n+2aFN,where superscripts indicate the time step, pF is the pressure at each node on the fluid surface caused by the applied acceleration, ρ is the density of the acoustic medium, and B is the B-spline impulse response matrix of the fluid surface, as determined in Section 2. Due to the previous assumption that the rigid body remains always in contact with the fluid, the pressure distribution along the wet surface of the rigid body can be computed as follows:(15)pWSN=−pFN=ρ∑n=1N+1BN−n+2aFN.In view of Equations (8)–(15), a rigid body boundary element is established that relates the hydrodynamic pressure that develops on the wet surface of the rigid body to the acceleration applied at the reference point of the rigid body.


## 3.1.2. B-Spline Impulse Response of the Rigid Body

The BIRF of a rigid body pertains to the hydrodynamic forces and moments that develop at the reference point due to a B-spline acceleration applied at the reference point and can be computed as follows. A unit B-spline acceleration is applied in the direction of each degree of freedom of the reference point, and Equation (15) computes the hydrodynamic pressure on all nodes of the wet surface. However, the differential hydrodynamic force dFn in the direction normal to the wet surface at any point is computed based on the continuous pressure field pWS at this point as follows:(16)dFWSn=pWSdA,where dA is the differential area. Based on the field variable interpolation, the element discretization, and the directional cosines n=lmn of the outward normal vector at any point on an element, the Fx, Fy, and Fz components of the differential force dF vector are computed as follows:(17)dFWS=nTNspWSsdA,s=1,2,…,8.Subsequently, the differential force,dFRP, and corresponding moment, dMRP, at the reference point of the rigid body are computed as follows:(18a)dFRP=dFWS=nTNspWSsdA,(18b)dMRP=HdFWS=HnTNspWSsdA.Therefore, the contribution of each element, el, to the forces,FRPel, and moments, MRPel, at the reference point is determined by integrating Equation (18) over the area of the element, i.e.,(19a)FRPel=∫SeldFRP=∫SelnTNspWSsdA,(19b)MRPel=∫SeldMRP=∫SelHnTNspWSsdA.Finally, the resultant hydrodynamic forces,FRP, and moments, MRP, at the reference point of the rigid body are computed as a mere superposition of the element contributions computed in Equation (19) as follows:(20a)FRP=∑el=1NEFRPel=∑el=1NE∫SeldFEP=∑el=1NE∫SelnTNspWSsdA,(20b)MRP=∑el=1NEMRPel=∑el=1NE∫SeldMRP=∑el=1NE∫SelHnTNspWSsdA.To compute BIRFs of the rigid body, all six degrees of freedom of a rigid body need to be considered. The BIRF,BRP, of the rigid body is a 6 × 6 matrix. In order to obtain this matrix, unit acceleration with B-spline time modulation is applied in the direction of each of the six degrees of freedom, one at a time, and the corresponding forces and moments at the reference point of the rigid body are computed. This response represents one column in the BIRF matrix. Collecting all responses in a matrix form yields the BIRF matrix of the rigid body as follows:(21)BRP=FRPu¨xFRPu¨yFRPu¨zFRPθ¨xFRPθ¨yFRPθ¨zMRPu¨xMRPu¨yMRPu¨zMRPθ¨xMRPθ¨yMRPθ¨z,where the superscript indicates the acceleration component that generated the hydrodynamic force at the reference point. It should be noted that the reference point can be located in a way such that the coupling effects between the different modes may be ignored. Consequently, the BIRF matrix is a diagonal matrix, and the six modes of motion are completely uncoupled.


## 3.2. Kinematic Fluid-Structure Interaction Model

The last step in the development of the kinematic FSI model computes hydrodynamic forces and moments due to arbitrary acceleration. Once the BIRF of the reference point of the rigid body is evaluated as discussed in Section3.1.2, the equivalent hydrodynamic forces and moments, Fhyd, at the reference point of the system to an arbitrary acceleration applied at the reference point are computed at step N as a mere superposition of the BIRF responses as follows:(22)FhydN=∑n=1N+1BRPN−n+2aRPN.Therefore, kinematic FSI models can be obtained corresponding to the six primary modes of motion of the rigid body, as well as the associated coupled modes, if required.


## 4. Dynamic FSI Models for Motion of Rigid Bodies

In this section, dynamic FSI models of the rigid body will be developed and discussed. The proposed dynamic models are 3D hydrodynamic fluid-structure interaction mathematical models that account for both inertia and kinematic interaction effects between a rigid structure and the surrounding fluid domain. The model is based on dynamic equilibrium and uses the characteristic BIRF of rigid bodies in fluid media, as developed in Section3. All six degrees of freedom discussed in the preceding section can be considered as required. The dynamic model is expressed in general by the equation of dynamic equilibrium of the system as follows:(23)Fnet=maRP,where m is the mass matrix in an FEM sense associated with the rigid body, and(24)aRP=u¨RPθ¨RPT=u¨RPxu¨RPyu¨RPzθ¨RPxθ¨RPyθ¨RPzTis the acceleration vector applied at the center of gravity of the rigid body. The force vector Fnet represents the net forces applied at the center of gravity of the rigid body and accounts for generalized driving forces and moments, Fexc, and hydrodynamic effects, Fhyd, and is expressed as follows:(25)Fnet=Fexc−Fhyd−Fother,where Fother represents other forces acting on the rigid body, such as drag forces, waves, and currents, among others. In this work, Fother is assumed to be zero, without loss of generality. The excitation forces, Fexc, are considered as the driving forces required accelerating the structure and are known. However, both the hydrodynamic forces, Fhyd, and the acceleration of the body, aRP, in Equation (23) are not known a priori. Consequently, Equation (23) cannot be solved for the body accelerations. In view of Equation (23) and provided that the BIRF of the rigid body is already computed, the following systems of equations need to be solved:(26a)Fexc−Fhyd=maRP,(26b)FhydN=∑n=1N+1BRPNaRPN−n+2.Different approaches can be used for solving Equation (26) such as simultaneous solution, staggered solution, and staggered solution with predictor-corrector schemes [1, 3, 28]. The simultaneous solution is the most straightforward and efficient approach and is adopted in this work, while the other two approaches are currently under investigation. Simultaneous solutions of Equation (26) combine the equations into a single system that can be solved in a time-marching scheme, implying that the time axes in both Equations (26a) and (26b) are discretized in a series of equal time steps i = 1, 2,…, N,…, M, where M is the total number of time steps that the solution is sought after. Subsequently, Equations (26a) and (26b) are both written for time step N. By separating known from unknown acceleration terms in Equation (26b), the latter is written as follows:(27)FhydN=BRP1aRPN+1+BRP2aRPN+∑n=3N+1BRPNaRPN−n+2,where aRPN and aRPN+1 are the unknown accelerations at steps N and N + 1, respectively. Following a procedure similar to the procedure for calculating the BIRF of the wet surface [23], the acceleration at the forward step is related to the acceleration at previous steps as follows:(28)aRPN+1=2aRPN−aRPN−1.Introducing Equation (28) into Equation (27) and factoring the unknown acceleration aRPN, the latter is expressed as follows:(29)FhydN=2BRP1+BRP2aRPN+∑n=3N+1BRPNaRPN−n+2−BRP1aRPN−1=mhydaRPN+HFhydN.Matrix,mhyd=2BRP1+BRP2, can be viewed as an equivalent hydrodynamic mass matrix, and vector HFhydN=∑n=3N+1BRPNaRPN−n+2−BRP1aRPN−1 is a force vector representing the effects of the history of the hydrodynamic forces on the current step N and is always known. Substituting Equation (29) into Equation (26a), the latter is expressed as follows:(30)FexcN=maRPN+mhydaRPN+HFhydN.Therefore, Equations (26a) and (26b) are combined into a single system of equations as follows:(31)F^N=m^aRPN,where F^N=FexcN−HFhydN and m^=m+mhyd. The system of equations shown in Equation (31) can be solved for the acceleration of the body in a time-marching scheme using conventional solvers for time integration. If the coefficient matrix, m^, is diagonal, indicating that the degrees of freedom are uncoupled, the proposed approach becomes very efficient. For example, such a case is met when only the surge, sway, and yaw modes are considered for a body with two vertical planes of symmetry and uniformly distributed mass lumped at the body’s center of gravity. However, if a coupled system is considered, the system of Equation (31) is still of a small size (6 degrees of freedom), the solution of which can be performed in a very efficient manner. The direct coupling of Equation (26) as proposed herein implies that the time step of the time-marching procedure for solving Equation (31) should be the one used for the calculation of the BIRF of the rigid body and presents a limitation of the adopted approach. The effects of this limitation will be assessed in the future work.


## 5. Verification Study

### 5.1. Problem Statement and Solution Process

This section presents the verification of the proposed model. The idealized body considered in this verification study is a rigid massive hollow sphere submerged in a fictitious fluid medium. A time history of the excitation force is applied at the center of gravity of the rigid spherical body, and the acceleration time history is computed by the proposed model and verified with results from the FEM models using the commercial software ANSYS. The rigid body is generated by attaching a rigid lining on a spherical cavity of radiusro = 212 m. The total mass of the rigid sphere is ms = 400,000 kg. The properties of the surrounding fictitious fluid medium are assumed as follows: (i) density ρf  = 1,000 kg/m3 and (ii) wave velocity c = 1,000 m/sec. The wet surface is discretized into 24 eight-node elements for a total of 74 boundary nodes, as shown in Figure 3. It is well established that the BEM does not require the dense discretization usually needed in FEM models for mesh convergence. A study on the mesh convergence of a sphere’s surface as it pertains to wave propagation in solids has been reported in [29]. These studies have shown the 24 quadratic elements used in the surface discretization yield solutions with a 0.1% mean error as compared to analytical solutions. The time history of the applied excitation force is also shown in Figure 3. Three time steps are used in the proposed study, i.e., Δt = 1/24 sec, Δt = 1/50 sec, and Δt = 1/80 sec. These steps correspond to three different 4th order B-spline excitations, each one of which is of duration 4Δt. The solution process consists of the following stages: (1) the BIRF of the fluid surface (wet surface) in the absence of the rigid lining is computed as described in Section 2 for each of the three B-spline excitations; (2) the BIRF of the rigid body is computed as discussed in Section 3. This BIRF represents the forces and moments at the center of gravity of the body when a B-spline acceleration is applied at the same point; and (3) the acceleration of the rigid body due to an applied excitation force is computed as detailed in Section 4. The three solution stages are presented in detail in Sections 5.2–5.4.Figure 3
Geometry (a), BEM discretization mesh of the wet surface (b), and time history of excitation force (c).
(a)
(b)
(c)


### 5.2. BIRF of Wet Surface and Verification

The BEM and the BIRF of the wet surface have been validated in [20]. It was demonstrated that the proposed BEM is highly accurate. For completeness, a summary of the verification study on the BIRF of the wet surface as reported in [23] is included here. First the BIRF of the wet surface is computed as discussed in Section 2, Equations (3) and (4), and three entries are shown in Figure 4(a). In particular, the time history of the BIRF of the pressure field at the North Pole, Equator, and South Pole is shown when a B-spline flux is applied at the North Pole.Figure 4
Wet surface solutions: (a) the BIRF of points on the “North Pole,” “Equator,” and “South Pole” due to B-spline flux at North Pole; (b) time history of the pressure field due to constant internal flux at specific points on the wet surface and comparison with analytical and other solutions in a nondimensional form where c is the pressure wave velocity,r is the radius of the sphere, t is time, p is the pressure field and q is the influx [23].
(a)
(b)Subsequently, the wet surface is subjected to a suddenly applied uniformly distributed internal flux, and the pressure field is computed as indicated in Equations (5)–(7) and is shown in Figure 4(b).The maximum error recorded is less than 1%, and the mean error is less than 0.1% compared to the analytical solution of the problem, as well as other solutions [30] reported in the literature. Extensive validation studies are published in [23].


### 5.3. BIRF of the Rigid Body

The next stage in the modeling process calculates the BIRF of the rigid body. In view of the excitation force being applied in thex direction at the reference point (center of sphere in this case), the BIRF of the rigid body pertains to the calculation of the time history of the force in the x direction that develops when acceleration in the x direction is applied at the reference point. This BIRF represents the first element in the BIRF shown in Equation (21). This is attributed to the symmetries of the rigid body and the particular location of its reference point that uncouples completely all vibration modes. It is noted that, in general cases, the dimension of the BIRF may be up to 6 × 6. To calculate the BIRF, unit acceleration with B-spline time modulation is applied in the x direction at the reference point of the rigid body, and the corresponding forces and moments at the reference point are computed, as discussed in Section 3.1.2. The assumption that the reference point coincides with the center of the sphere is not limiting, although it does not explicitly demonstrate the rotational vibration modes. It is noted that the resulting BIRF moments as computed in the present stage are zero for all practical purposes providing, thus, a strong verification that the formulation is valid. The BIRFs of the rigid sphere are shown in Figure 5 for B-spline excitation of the three different durations considered in this study. As expected, the longer duration B-spline excitations cause the higher amplitude responses and longer forced vibration response phases.Figure 5
B-spline impulse response functions of the rigid sphere for B-spline excitations of three different durations.


### 5.4. Response to Arbitrary Excitation Force

The last step in the analysis is the calculation of the acceleration of the body in response to a transient excitation force. It is noted that the tasks of computing the BIRF of the wet surface and the BIRF of the rigid body are independent of the actual excitation force and, thus, need be calculated only once for a given geometry of the rigid body. Given the time history of the excitation force shown in Figure3, Equations (26a) and (26b) are solved simultaneously as a single system of equations shown in Equation (31), for the acceleration of the body. Figures 6(a) and 6(b) show the calculated velocity and acceleration response time history as calculated based on two of the BIRFs of the rigid body shown in Figure 5, namely, the ones that correspond to time steps of Δt = 1/50 sec and 1/80 sec. The two BEM solutions are indistinguishable, the solution is stable even at later times, and the convergence of the solution is evident. The accuracy of the BEM solution is verified through comparisons with an FEM solution that is also shown in Figures 6(a) and 6(b). The details of the FEM model are discussed in the next section.Figure 6
Response of the rigid body to a transient excitation force: (a) velocity; (b) acceleration.
(a)
(b)


#### 5.4.1. Verification with ANSYS FEM Solution

The solution of the proposed FSI technique shown in Figure6 is verified with FEM solutions. To this end, 3D finite element models of the rigid sphere submerged in the fictitious fluid media are generated using the finite element commercial software ANSYS. The model consists of (i) the rigid sphere, (ii) the truncated surrounding fluid, and (iii) the solid-fluid interface. The rigid sphere is discretized using 8-node SOLID45 solid elements, and the fluid is discretized using 8-node FLUID30 fluid elements. The mesh ratio in the radial direction is 0.1 in order to create a denser mesh near the rigid surface to capture the higher variation of the field variables. The fluid domain needs to be truncated. To this end, the surface element FLUID130 is defined along the truncation surface of the fluid domain to simulate the infinite extents of the fluid. The fluid-structure interface boundary condition option is activated in ANSYS to represent the transfer of forces between the solid and the fluid at the wet surface. The rigid sphere is assumed as a highly stiff body in this study by using a high Young’s modulus value of 2.068 × 1011 k/m2 and zero Poisson’s ratio. The material properties of the fictitious fluid medium are defined by the density of the fluid ρf  = 1,000 kg/m3 and the wave velocity c = 1,000 m/sec. A series of studies on mesh and time step convergence have been conducted, and an acceptable solution with respect to convergence and computational effort was obtained for a time step dt = 1/50 sec and a relatively fine mesh with the infinite boundary located at 600 m away from the wet surface. The location of the infinite boundary at 600 m from the sphere’s center is determined based on a parametric study. This location is approximately three times the radius of the sphere and produced results that showed negligible, if any, influence from waves reflected at the infinite boundary and is consistent with the recommendations in the software manual. The transient excitation force is applied at the node located at the center of the rigid sphere in the FEM model, and the time history of the velocity and acceleration of the body is obtained and plotted in Figures 6(a) and 6(b). It is observed that the two methods are in good agreement. It is noted that FEM models require much higher mesh density than the proposed BEM-based technique to achieve a similar level of accuracy and that improved solutions may be obtained by further refinement of the mesh and the use of smaller time steps; however, such tasks become prohibitively expensive.


## 5.1. Problem Statement and Solution Process

This section presents the verification of the proposed model. The idealized body considered in this verification study is a rigid massive hollow sphere submerged in a fictitious fluid medium. A time history of the excitation force is applied at the center of gravity of the rigid spherical body, and the acceleration time history is computed by the proposed model and verified with results from the FEM models using the commercial software ANSYS. The rigid body is generated by attaching a rigid lining on a spherical cavity of radiusro = 212 m. The total mass of the rigid sphere is ms = 400,000 kg. The properties of the surrounding fictitious fluid medium are assumed as follows: (i) density ρf  = 1,000 kg/m3 and (ii) wave velocity c = 1,000 m/sec. The wet surface is discretized into 24 eight-node elements for a total of 74 boundary nodes, as shown in Figure 3. It is well established that the BEM does not require the dense discretization usually needed in FEM models for mesh convergence. A study on the mesh convergence of a sphere’s surface as it pertains to wave propagation in solids has been reported in [29]. These studies have shown the 24 quadratic elements used in the surface discretization yield solutions with a 0.1% mean error as compared to analytical solutions. The time history of the applied excitation force is also shown in Figure 3. Three time steps are used in the proposed study, i.e., Δt = 1/24 sec, Δt = 1/50 sec, and Δt = 1/80 sec. These steps correspond to three different 4th order B-spline excitations, each one of which is of duration 4Δt. The solution process consists of the following stages: (1) the BIRF of the fluid surface (wet surface) in the absence of the rigid lining is computed as described in Section 2 for each of the three B-spline excitations; (2) the BIRF of the rigid body is computed as discussed in Section 3. This BIRF represents the forces and moments at the center of gravity of the body when a B-spline acceleration is applied at the same point; and (3) the acceleration of the rigid body due to an applied excitation force is computed as detailed in Section 4. The three solution stages are presented in detail in Sections 5.2–5.4.Figure 3
Geometry (a), BEM discretization mesh of the wet surface (b), and time history of excitation force (c).
(a)
(b)
(c)


## 5.2. BIRF of Wet Surface and Verification

The BEM and the BIRF of the wet surface have been validated in [20]. It was demonstrated that the proposed BEM is highly accurate. For completeness, a summary of the verification study on the BIRF of the wet surface as reported in [23] is included here. First the BIRF of the wet surface is computed as discussed in Section 2, Equations (3) and (4), and three entries are shown in Figure 4(a). In particular, the time history of the BIRF of the pressure field at the North Pole, Equator, and South Pole is shown when a B-spline flux is applied at the North Pole.Figure 4
Wet surface solutions: (a) the BIRF of points on the “North Pole,” “Equator,” and “South Pole” due to B-spline flux at North Pole; (b) time history of the pressure field due to constant internal flux at specific points on the wet surface and comparison with analytical and other solutions in a nondimensional form where c is the pressure wave velocity,r is the radius of the sphere, t is time, p is the pressure field and q is the influx [23].
(a)
(b)Subsequently, the wet surface is subjected to a suddenly applied uniformly distributed internal flux, and the pressure field is computed as indicated in Equations (5)–(7) and is shown in Figure 4(b).The maximum error recorded is less than 1%, and the mean error is less than 0.1% compared to the analytical solution of the problem, as well as other solutions [30] reported in the literature. Extensive validation studies are published in [23].


## 5.3. BIRF of the Rigid Body

The next stage in the modeling process calculates the BIRF of the rigid body. In view of the excitation force being applied in thex direction at the reference point (center of sphere in this case), the BIRF of the rigid body pertains to the calculation of the time history of the force in the x direction that develops when acceleration in the x direction is applied at the reference point. This BIRF represents the first element in the BIRF shown in Equation (21). This is attributed to the symmetries of the rigid body and the particular location of its reference point that uncouples completely all vibration modes. It is noted that, in general cases, the dimension of the BIRF may be up to 6 × 6. To calculate the BIRF, unit acceleration with B-spline time modulation is applied in the x direction at the reference point of the rigid body, and the corresponding forces and moments at the reference point are computed, as discussed in Section 3.1.2. The assumption that the reference point coincides with the center of the sphere is not limiting, although it does not explicitly demonstrate the rotational vibration modes. It is noted that the resulting BIRF moments as computed in the present stage are zero for all practical purposes providing, thus, a strong verification that the formulation is valid. The BIRFs of the rigid sphere are shown in Figure 5 for B-spline excitation of the three different durations considered in this study. As expected, the longer duration B-spline excitations cause the higher amplitude responses and longer forced vibration response phases.Figure 5
B-spline impulse response functions of the rigid sphere for B-spline excitations of three different durations.


## 5.4. Response to Arbitrary Excitation Force

The last step in the analysis is the calculation of the acceleration of the body in response to a transient excitation force. It is noted that the tasks of computing the BIRF of the wet surface and the BIRF of the rigid body are independent of the actual excitation force and, thus, need be calculated only once for a given geometry of the rigid body. Given the time history of the excitation force shown in Figure3, Equations (26a) and (26b) are solved simultaneously as a single system of equations shown in Equation (31), for the acceleration of the body. Figures 6(a) and 6(b) show the calculated velocity and acceleration response time history as calculated based on two of the BIRFs of the rigid body shown in Figure 5, namely, the ones that correspond to time steps of Δt = 1/50 sec and 1/80 sec. The two BEM solutions are indistinguishable, the solution is stable even at later times, and the convergence of the solution is evident. The accuracy of the BEM solution is verified through comparisons with an FEM solution that is also shown in Figures 6(a) and 6(b). The details of the FEM model are discussed in the next section.Figure 6
Response of the rigid body to a transient excitation force: (a) velocity; (b) acceleration.
(a)
(b)


### 5.4.1. Verification with ANSYS FEM Solution

The solution of the proposed FSI technique shown in Figure6 is verified with FEM solutions. To this end, 3D finite element models of the rigid sphere submerged in the fictitious fluid media are generated using the finite element commercial software ANSYS. The model consists of (i) the rigid sphere, (ii) the truncated surrounding fluid, and (iii) the solid-fluid interface. The rigid sphere is discretized using 8-node SOLID45 solid elements, and the fluid is discretized using 8-node FLUID30 fluid elements. The mesh ratio in the radial direction is 0.1 in order to create a denser mesh near the rigid surface to capture the higher variation of the field variables. The fluid domain needs to be truncated. To this end, the surface element FLUID130 is defined along the truncation surface of the fluid domain to simulate the infinite extents of the fluid. The fluid-structure interface boundary condition option is activated in ANSYS to represent the transfer of forces between the solid and the fluid at the wet surface. The rigid sphere is assumed as a highly stiff body in this study by using a high Young’s modulus value of 2.068 × 1011 k/m2 and zero Poisson’s ratio. The material properties of the fictitious fluid medium are defined by the density of the fluid ρf  = 1,000 kg/m3 and the wave velocity c = 1,000 m/sec. A series of studies on mesh and time step convergence have been conducted, and an acceptable solution with respect to convergence and computational effort was obtained for a time step dt = 1/50 sec and a relatively fine mesh with the infinite boundary located at 600 m away from the wet surface. The location of the infinite boundary at 600 m from the sphere’s center is determined based on a parametric study. This location is approximately three times the radius of the sphere and produced results that showed negligible, if any, influence from waves reflected at the infinite boundary and is consistent with the recommendations in the software manual. The transient excitation force is applied at the node located at the center of the rigid sphere in the FEM model, and the time history of the velocity and acceleration of the body is obtained and plotted in Figures 6(a) and 6(b). It is observed that the two methods are in good agreement. It is noted that FEM models require much higher mesh density than the proposed BEM-based technique to achieve a similar level of accuracy and that improved solutions may be obtained by further refinement of the mesh and the use of smaller time steps; however, such tasks become prohibitively expensive.


## 5.4.1. Verification with ANSYS FEM Solution

The solution of the proposed FSI technique shown in Figure6 is verified with FEM solutions. To this end, 3D finite element models of the rigid sphere submerged in the fictitious fluid media are generated using the finite element commercial software ANSYS. The model consists of (i) the rigid sphere, (ii) the truncated surrounding fluid, and (iii) the solid-fluid interface. The rigid sphere is discretized using 8-node SOLID45 solid elements, and the fluid is discretized using 8-node FLUID30 fluid elements. The mesh ratio in the radial direction is 0.1 in order to create a denser mesh near the rigid surface to capture the higher variation of the field variables. The fluid domain needs to be truncated. To this end, the surface element FLUID130 is defined along the truncation surface of the fluid domain to simulate the infinite extents of the fluid. The fluid-structure interface boundary condition option is activated in ANSYS to represent the transfer of forces between the solid and the fluid at the wet surface. The rigid sphere is assumed as a highly stiff body in this study by using a high Young’s modulus value of 2.068 × 1011 k/m2 and zero Poisson’s ratio. The material properties of the fictitious fluid medium are defined by the density of the fluid ρf  = 1,000 kg/m3 and the wave velocity c = 1,000 m/sec. A series of studies on mesh and time step convergence have been conducted, and an acceptable solution with respect to convergence and computational effort was obtained for a time step dt = 1/50 sec and a relatively fine mesh with the infinite boundary located at 600 m away from the wet surface. The location of the infinite boundary at 600 m from the sphere’s center is determined based on a parametric study. This location is approximately three times the radius of the sphere and produced results that showed negligible, if any, influence from waves reflected at the infinite boundary and is consistent with the recommendations in the software manual. The transient excitation force is applied at the node located at the center of the rigid sphere in the FEM model, and the time history of the velocity and acceleration of the body is obtained and plotted in Figures 6(a) and 6(b). It is observed that the two methods are in good agreement. It is noted that FEM models require much higher mesh density than the proposed BEM-based technique to achieve a similar level of accuracy and that improved solutions may be obtained by further refinement of the mesh and the use of smaller time steps; however, such tasks become prohibitively expensive.


## 6. Conclusions

This work presents the development of an efficient model for dynamic fluid-structure interaction and the associated rigid body dynamic model. The proposed models consider the structure and the fluid combined as one system, while the excitation forces required to accelerate the body are treated as external forces to the system. The fluid medium is assumed to be a homogeneous, inviscid, compressible fluid with small displacements, and the structure is assumed to be rigid. The following tasks have been discussed: (i) a brief description of development of an advanced direct time-domain boundary element method for wave propagation problems in acoustic media; (ii) development of simplified models for the fluid-rigid structure kinematic interaction based on impulse response techniques and the BEM; and (iii) the development of models for the motion of the rigid body that integrate the kinematic interaction models into dynamic models in an efficient manner for real-time simulations. The following conclusions can be drawn from the development of these models:(i)
The B-spline direct time-domain boundary element formulation for wave propagation in fluid media is proved to be highly accurate, efficient, and stable.(ii)
The kinematic models are represented by the six degrees of freedom of a reference point. The models relate hydrodynamic forces and moments at the reference point of a rigid body to the acceleration applied at the same point. The reference point can be chosen so that the degrees of freedom are in general uncoupled, yielding an extremely efficient model for real-time simulations.(iii)
The proposed dynamic models are accurate as demonstrated through validations with FEM models. It should be noted that the proposed BEM is much more efficient than FEM both in terms of data preprocessing tasks and the reduced size of the problem.(iv)
The versatility and efficiency of the proposed models are demonstrated through an application example where the acceleration response of the structure is computed faster than real time.



---

*Source: 9760361-2018-10-11.xml*