## Abstract

The paper presents a novel method for numerically modelling fluid–structure interactions. The method consists of solving the fluid-dynamics equations on an extended domain, where the computational mesh covers both fluid and solid structures. The fluid and solid velocities are relaxed to one another through a penalty force. The latter acts on a thin shell surrounding the solid structures. Additionally, the shell is represented on the extended domain by a non-zero shell-concentration field, which is obtained by conservatively mapping the shell mesh onto the extended mesh. The paper outlines the theory underpinning this novel method, referred to as the immersed-shell approach. It also shows how the coupling between a fluid- and a structural-dynamics solver is achieved. At this stage, results are shown for cases of fundamental interest.

## 1. Introduction

The analysis of fluid–structure interactions (FSIs) is important in many areas, for example, to simulate the dynamics of turbine blades, the wave impacts on floating bodies, the blood flow through arteries or pollutant dispersion in street canyons [1–3]. Numerical models are attractive in this context, because they can tackle different configurations while limiting expensive laboratory testing. For flows past rigid and still structures, the Navier–Stokes equations are usually solved on a mesh excluding the solid structures (defined-body (DB) method). This technique can give accurate results, when using appropriate discretization schemes and sufficient spatial resolution, because the structure shape and boundary conditions are represented exactly. However, when the structure moves in the fluid domain, re-meshing is necessary. This is the case, for example, for wind turbine rotors or floating wind turbines. In this context, the whole structure dynamically interacts with the surrounding wind and waves. Re-meshing is computationally expensive and might also yield highly distorted grids. In order to avoid this problem, another approach was developed. It consists of meshing the entire domain (containing both fluid and structures) and modelling the effect of the structures through surface or body forces. Because the effect of the structures is modelled through an additional forcing term in the fluid's momentum equation, the methods can be further divided into two categories depending on whether the force is added to the momentum equation before discretization (*continuous forcing*) or after discretization (*direct forcing*). This methodology underpins the so-called immersed-boundary method, proposed by Peskin [2] in a continuous-forcing context for modelling cardiac mechanics. Reviews and applications of immersed-boundary type methods are available in the literature [4–6]. Among the variants to Peskin's original formulation, Goldstein *et al.* [7] proposed to attach virtual springs and dampers to the immersed boundary, so that the forcing term acts as a feedback mechanism and tends to oppose any deviations from the desired values of velocity and position of the structure. This method suffers from limitations in the time stepping size, owing to the need to resolve the characteristic time scales of the spring–damper system. In order to overcome this drawback, Fadlun *et al.* [8] directly forced the desired velocity at the boundary nodes by modifying the entries in the matrix of the discrete momentum equations. This direct-forcing method, however, yielded force oscillations due to the interpolation method between the fixed Eulerian mesh and the moving boundary nodes. Uhlmann [9] solved this issue by combining Peskin's original formulation with the direct and explicit formulation of the FSI force. In this method, a volume was associated with the Lagrangian points placed on the structure boundary, hence forming a volumetric shell of one mesh width around the structure. The Lagrangian points followed a rigid-body motion, and smooth variants of Peskin's Dirac function were used to transfer information between Cartesian and Lagrangian points. Kempe & Fröhlich [10] extended this direct-forcing method to further improve the imposition of the boundary condition at the interface. These methods are particularly attractive for modelling particulate flows, for which numerical efficiency, particle collisions and adequate representation of the boundary layers are important. Recently, Viré *et al.* [11] presented an immersed-body method with a continuous-forcing approach based on a penalty forcing. In this context, the fluid and structural velocities were weakly relaxed to one another in the regions covered by the structures. This method relied on two distinct meshes used by two separate finite-element models: one mesh covering the entire domain (fluid mesh), and another mesh discretizing solely the structures (solid mesh). This method is a versatile way of modelling interactions between fluids and structures, because the use of separate models allows specific numerical requirements to be met for each material. One of the difficulties, however, is to ensure that the laws of physics are verified at the discrete level when exchanging information between the models. The authors proposed a novel algorithm to ensure spatial conservation of the penalty force and solid concentration field when projected from one mesh to the other. Importantly, spatial conservation was verified at a discrete level, and therefore was independent of the mesh resolution and the types of mesh used in the two models. This approach was validated on flows past rigid and deformable structures, as well as fixed actuator discs representing turbines [11,12]. Although the immersed-body approach has proven to be successful in modelling FSIs for various applications, it applies the coupling force across the regions covered by the structures. Thus, the action–reaction force between fluids and structures—although very small—is not strictly equal to zero within the structures. This is a common issue in immersed-body type methods. However, this is not physically meaningful and can potentially affect the level of stresses computed by the structural-dynamics model. In addition, the discrete accuracy of the force acting on the structure depends on the resolution of the solid mesh. This limits the number of structures that can be modelled in one simulation, owing to the computational cost associated with accurately discretizing each individual structure.

The present immersed-shell (IS) method aims at improving the representation of the action–reaction force at the fluid–structure interface. The proposed methodology has some similarities with both the immersed-body approach [11] and the direct-forcing methods [9,10]. First, a penalty force is added to the momentum balances of each material in order to relax the fluid and structural velocities to one another. Second, the force acts in the regions of non-zero solid concentrations. However, as opposed to the immersed-body method, the exchange of forces between fluid- and structural-dynamics models is entirely done via a shell surrounding the structures. As a result, the penalty force exactly vanishes within the structures. A better discrete representation of the boundary layer can also be achieved. Finally, as opposed to the existing direct-forcing methods [9,10], the shell thickness can be chosen independently of the mesh size. The projection of fields between the Eulerian fluid mesh and the shell mesh is also done differently, as explained herein. The paper presents the general formalism of the method and validates it on flows of fundamental interest. Further work will focus on implementing different wall parametrizations in the shell region, in order to tackle high Reynolds number flows that are relevant for practical aerodynamics applications. The aim of this paper is twofold. First, it presents the theoretical concepts underpinning the IS method. Section 2 focuses on the set of equations solved by the coupled models, whereas §3 details the projection steps required to transfer fields between fluid, shell and solid meshes. The sequence of steps required in the coupling approach is summarized in §4. Second, the paper demonstrates the capabilities of the present method on flows of fundamental interest (§5). The method is validated on flow past a fixed two-dimensional cylinder, an oscillating cylinder and a free-falling sphere. Finally, conclusions are drawn in §6, and a strategy is proposed to further apply the method to engineering problems, such as innovative offshore wind energy concepts.

## 2. Mathematical formulation

The fluid/ocean dynamics model ‘Fluidity’ is a finite-element open-source (LPGL) numerical tool that solves the non-hydrostatic Navier–Stokes equations on unstructured meshes. It can capture small-scale structures and large-scale processes simultaneously by adapting the meshes dynamically in time [13–15]. In this work, the Navier–Stokes equations are modified to account for the presence of immersed solids in the fluid domain. The idea underpinning the IS method is to solve the Navier–Stokes equations over a fluid mesh covering the extended domain *V* =*V* _{f}∪*V* _{s} (where the computational domain *V* contains fluids in *V* _{f} and structures in *V* _{s}). The first step is to fill the regions covered by the structures with the surrounding fluid, so that the fluid velocity *u*_{f} exists across the computational domain *V* . In this context, the continuity equation in the extended domain *V* states that
2.1
The second step consists of representing the effect of the structures on the flow. This is achieved by relaxing the fluid and structural velocities to one another in a thin shell surrounding the structures. The shell is meshed separately from the extended domain. A shell concentration field *α*_{sh} identifies, on the fluid mesh, the region occupied by the shell. It is computed by conservatively mapping the shell mesh onto the fluid mesh, as explained in §3. The relaxation between fluid and structural velocities is performed through a penalty force, which is non-zero in the fluid elements where the shell concentration field is non-zero. The penalty force *F*_{f} is added to the right-hand side of the fluid’s momentum equation, which yields
2.2
where *ρ*_{f} is the fluid density, *σ*_{f} is the total stress tensor (including the pressure contribution) and *B*_{f} represents additional body forces (e.g. due to gravity). The penalty force is expressed as
2.3
where the fluid *u*_{f} and structural *u*_{s} velocities are expressed in the thin shell, and
2.4
is a relaxation factor that dictates how fast the fluid and structural velocities equal one another at the interface. In equation (2.4), Δ*t* is the time step, *L* is the local edge length and *γ*=*l*_{e}/Δ is the ratio between the minimum edge length *l*_{e} of the fluid mesh around the shell and the shell thickness Δ, when the fluid mesh does not resolve the shell. In the limiting case where the shell is resolved (*l*_{e}<Δ), *γ*=1 is considered. In addition, the magnitude of the relaxation factor *β* is driven by viscous effects at small Reynolds numbers and inertial effects at large Reynolds numbers. In this study, it is assumed to be driven solely by inertial effects.

The structural dynamics is solved on a separate model. The combined finite–discrete element model ‘Y3D’ is coupled to ‘Fluidity’ and solves the deformable body equations using a finite-strain formulation [16,17]. This can compute the stresses and vibration modes of structures of any shape and stiffness. The equations governing the dynamics of deformable structures on the solid mesh are given by
2.5
where *d*_{s} denotes the solid displacement, is an operator dependent on the velocity gradient, *F*_{int} stands for the internal forces, *F*_{c} is the contact force when multiple structures impact on one another, and *F*_{ext} is the external force including the surface traction force and all the body forces other than the action–reaction force exerted by the fluid. In order to account for the presence of the surrounding fluid, the penalty force is added to the right-hand side of the structure’s momentum balance. Similar to the force in the fluid-dynamics equations, the penalty force acting on the structures is expressed as
2.6
By virtue of the action–reaction principle, the forces are such that
2.7
at a discrete level, as explained in §3. The equations of motion are discretized using finite elements. Time is discretized using a Crank–Nicolson scheme in the fluid-dynamics model, and a backward Euler scheme in the structural-dynamics model. The specific discretization used in ‘Fluidity’ and ‘Y3D’ is explained in reference [11]. The expression for the penalty force is analogous to the discontinuous Galerkin (DG) form of the viscous stresses, as highlighted hereafter. In the momentum balance of the fluid outlined by equation (2.2), the viscous stress tensor *σ*_{f} is given by
2.8
where the Kronecker symbol has the usual meaning, i.e. *δ*_{ii}=1 and *δ*_{ij}=0 for *i*≠*j*. By denoting the velocity components by *u*_{f}=(*u*,*v*,*w*), the viscous stress tensor expands to
2.9
and similarly for the other components. Figure 1*a* illustrates two neighbouring finite elements, where the solution fields are discontinuous at the nodes sharing the boundary between *E*^{−} and *E*^{+}. The DG discretization of the viscous force on such elements is given by
2.10
2.11
2.12
where *N*_{i} is the basis function at node *i*. After replacing the components of the stress tensor in equation (2.10) by their values from equation (2.9) (and similarly for *v* and *w*), DG methods commonly simplify the aforementioned equations by expressing the surface integrals as
2.13
where *α*_{xx}, *α*_{yy} and *α*_{zz} are penalization factors. Consequently, the above expressions penalize the jump in each component of the velocity at the element boundaries. This method of coupling the stress terms between the elements of a DG mesh is similar to the present coupling term that penalizes the discontinuity between fluid and structural velocities in the shell region (figure 1*b*).

## 3. Spatial projections between fluid, solid and shell meshes

In the equations presented in §2, the fluid velocity and pressure are computed on the fluid mesh, whereas the structural velocity is computed on the solid mesh. In order to couple these quantities, and compute the penalty force on each mesh, the fields need to be projected from fluid to shell mesh (and vice versa), and from the volumetric shell mesh to the solid surface mesh where the fluid forces are applied. This section shows how these projections are performed. The sequence in which the projections are done, and the fields they apply to, is explained in §4.

In order to project a field from shell to fluid mesh, the intersections between the two meshes are first identified. This can be done using either an R-tree spatial indexing algorithm [18] or an advancing front algorithm [19]. For each element of the shell mesh, a so-called *supermesh* is formed from the set of all the intersections, following a local supermeshing method [20]. Importantly, the order of the basis functions on the supermesh is equal to the maximum of their orders on the fluid and shell meshes. As a result, the projected field can be represented accurately on the supermesh. The mesh-to-mesh projection is performed via a Galerkin projection. The latter is used to obtain the shell concentration field *α*_{sh} on the fluid mesh, and also to project the coupling force from shell to fluid mesh. Given a donor mesh *D* and a target mesh *T* on a domain *Ω*, a Galerkin projection on a field *q* from *D* to *T* ensures that
3.1
by minimizing the *L*_{2} norm of the interpolation error [20]. In a discrete form, equation (3.1) involves solving a linear system of equations containing a mass matrix based on the basis functions of the target mesh and a mixed mass matrix formed of the basis functions of the target and donor meshes. For example, the shell concentration field *α*_{sh} on the fluid mesh is computed as
3.2

where the shell concentration field on the shell mesh is by definition a unitary field , and *n*_{f} (resp. *n*_{sh}) denotes the number of nodes in the fluid (resp. shell) mesh. In addition, superscripts refer to the mesh on which the quantity is defined (‘f’ for the fluid mesh and ‘sh’ for the shell mesh). Such a Galerkin projection ensures spatial conservation, but modifies the bounds of the projected quantity. It can be shown that lumping the mass matrix of the target mesh bounds the resulting solution without affecting the conservation property. The detailed projection methods between fluid and shell meshes are identical to those used in the immersed-body method [11]. This supermeshing approach differs from the transfer method used in similar direct-forcing methods [9,10,21]. The latter use a smoothed Dirac delta function over one mesh element around the structure, and therefore the interpolation between Lagrangian points and Eulerian nodes is performed via a weighted sum of regularized Dirac functions.

Once the fluid forces are known on the shell mesh, they need to be applied on the solid surface. Similarly, the solid quantities at the solid surface need to be mapped to the shell, before being projected to the fluid mesh. These procedures are described hereafter. The shell mesh is directly constructed from the triangular mesh discretizing the solid surface. This is illustrated in figure 2, where the shaded triangular element 1–2–3 lies on the solid surface. The shell is formed by extruding this element in the normal direction over a length Δ, which is the shell thickness. Thus, for each surface element at the solid surface, three tetrahedral elements are added to form the shell volumetric mesh, i.e. element 1–2–3–1′, element 2–2′–3–1′ and element 2′–3′–3–1′. When projecting a field from the solid surface to the volumetric shell, the shell is assumed to be much thinner than the characteristic length of the solid. Consequently, it is a good approximation to consider that the projected field *q* at a node on the outer shell equals the field at the corresponding node on the inner shell (*q*_{1′}=*q*_{1}). For example, the structure velocity at the solid surface is constant across the shell before projecting it to the fluid model. This assumption is not made when projecting a field from a volumetric shell to a solid surface. For the latter projection, the condition to satisfy is
3.3
where *S*_{s} denotes the surface of the structure and superscript notation identifies the mesh on which the field lies (i.e. ‘sh’ for the shell mesh and ‘surf’ for the surface mesh of the structure). At each shell nodes, the left-hand side of equation (3.3) is given by . The discrete values at the nodes 1′–2′–3′, which are located on the outer shell boundary, need to be mapped on the corresponding nodes 1–2–3 situated on the inner shell boundary. The mapping is such that the projected field at node 1 equals *q*_{1}*V* _{sh,1}+*q*_{1′}*V* _{sh,1′}.

## 4. Coupling algorithm between fluid- and structural-dynamics models

This section explains in more detail the coupling mechanism between the fluid- and structural-dynamics equations. The coupling is achieved through the penalty force across the shell, which ensures that, at a discrete level, the penalty force on the fluid mesh *F*_{f} is given by
4.1
In equation (4.1), denotes the conservative Galerkin projection from shell to fluid mesh of a field *q*, as detailed in §3. The detailed sequence of steps of the IS method is summarized hereafter. Steps (1)–(3) are done before the simulation starts, whereas the other steps are repeated at every time step.

(1) The computational domain covering both fluids and structures is meshed by the so-called fluid mesh, which is used to spatially discretize equation (2.2).

(2) The structural domain is meshed by the so-called solid mesh, which is used to spatially discretize equation (2.5).

(3) A shell is constructed around the structures (§3). The shell mesh is used to exchange the coupling term between fluids and structures. Note that the shell mesh is superimposed on the fluid mesh, and, importantly, does not necessarily coincide with the elements of the fluid mesh.

(4) A unitary field is projected from the shell mesh to the fluid mesh using a Galerkin projection. This yields a shell concentration field

*α*_{sh}on the fluid mesh, such that 4.2 where*V*_{sh}is the shell volume.(5) Two fields are projected from fluid to shell mesh: the pressure

*p*, and the part of the penalty force that depends on the fluid velocity, i.e. .(6) The pressure

*p*and force*F*^{(1)}_{s}are mapped from the shell to the solid surface (§3).(7) The structural velocity

*u*_{s}is computed by solving the structural-dynamics equations on the solid mesh, see equation (2.5). The surface integral of*p*represents the buoyancy force and dynamic pressure forces exerted by the surrounding fluid. It is included in the external forces in equation (2.5). The penalty force is also applied as a surface force in the structure’s momentum balance. Subcycling is required if the time step used to discretize the structural-dynamics equations is smaller than that used to discretize the fluid-dynamics equations. The total penalty force is updated based on the structural velocity at the previous time iteration*i*, as .(8) At the end of the structural time stepping loop, a time-averaged structural velocity is computed from the values at each time iteration

*i*.(9) The penalty force is mapped from the solid to the shell mesh, as explained in §3. The force on the shell mesh is then projected to the fluid mesh using a conservative Galerkin projection (§3).

(10) The fluid velocity

*u*_{f}is computed by solving the fluid-dynamics equations on the fluid mesh, equation (2.2).

## 5. Results

The algorithm described in the previous section is applied to three cases: flow past a two-dimensional cylinder, forced oscillation of a two-dimensional cylinder in a fluid at rest, and a fully coupled problem of flow past a sphere falling under the effect of gravity.

### (a) Flow past a two-dimensional cylinder

Flow past a two-dimensional cylinder is chosen as a verification/validation case owing to the number of references available in the literature. The computational configuration is as follows. The two-dimensional cylinder of diameter *D* is centred at a distance of 5*D* away from the inlet, 10*D* away from the sides, and 20*D* away from the outlet. It is subjected to a uniform flow of streamwise velocity *u*_{0}, and four Reynolds numbers are considered: *Re*_{D}=*u*_{0}*D*/*ν*=40, 100, 200, 500. In all the simulations, the Courant number is fixed at *CFL*=0.1. The fluid mesh is unstructured and does not change in the course of time. The typical value of the element edge lengths at the outer boundaries is *L*_{e}=*D*, whereas the minimum edge length *l*_{e} near the cylinder is indicated in table 1. The domain is discretized using a mixed finite-element discretization method, where discontinuous linear polynomials are used for the velocity field (i.e. P1-DG discretization) and continuous piecewise quadratic polynomials are used for the pressure field (i.e. P2 discretization). A well-known difficulty in representing solid structures as immersed in an extended computation domain is the fact that sharp fluid–structure interfaces are smeared on the computational cells. This affects the accuracy with which forces, and hence boundary conditions, are represented on the structure surface. A related issue concerns the computation of diagnostic variables. In particular, the integral of the pressure and viscous forces on the bluff-body surface is not straightforward owing to the smeared representation of the structure. Despite extensive literature on the numerical modelling of FSIs, this post-processing difficulty has been scarcely addressed [22,23]. This study shows that one of the main differences between DB and immersed-body type approaches is the location of the pressure extrema. In the DB approach, the maximum pressure occurs at the cylinder boundary, which is in accordance with experimental observations. In the IS approach, the pressure lines close in the shell as shown by figure 3. Thus, the pressure extrema are located close to the first node of the fluid mesh where the shell concentration field equals zero. This occurs at a distance *R*^{⋆}≈*R*+Δ+*l*_{e} from the cylinder centre, where *R* is the cylinder radius, Δ is the thickness of the shell surrounding the cylinder and *l*_{e} is the edge length of the fluid-mesh elements around the shell. This is only an estimation of *R*^{⋆}, because it further depends on where the shell and fluid meshes intersect one another within one fluid element. As a result, the aerodynamic forces in the IS results are obtained by integrating the pressure and viscous forces on the surface of the fictitious cylinder of diameter *D*^{⋆}=2*R*^{⋆}. Table 1 shows the drag coefficient obtained using the DB and IS approaches. The results obtained by the two methods are very close, and also match the values reported in the literature [24,25]. This is also apparent from the time evolutions of the lift *C*_{L} and drag *C*_{D} coefficients at *Re*=200, as shown by figure 4. The continuous line shows the result for the IS method, circles show the solution of Rajani *et al.* [26], and triangles show the solution of Meneghini *et al.* [27]. The agreement is acceptable, also considering the dispersion of the results given by other studies, see for example [28]. Finally, based on the frequency of oscillations of the drag coefficient at *Re*=100, the Strouhal number was found to equal *St*=*fD*/*u*_{0}=0.166, which is also in agreement with the literature.

### (b) Oscillating two-dimensional cylinder in a fluid at rest

The in-line oscillation of a cylinder in a fluid at rest is considered in this section. This is a classical benchmark problem for studying FSIs with moving objects and is well documented in the literature [29–31]. This test case is also chosen to show convergence of the numerical results with increasing resolution and decreasing shell thickness. The cylinder of diameter *D* is centred in the computational domain of size 55*D*×35*D*. The domain is discretized with an unstructured mesh, and the maximum element edge length on the outer boundaries equals *L*_{e}=*D*. The mesh is refined with a minimum element edge length in a region of size 7*D*×4*D* around the cylinder. Convergence of the results is shown by varying the minimum element edge length *l*_{e}, as indicated in table 2, in order to assess the effect of the mesh resolution on the computational results. The translational motion of the cylinder is prescribed as , where *A* is the oscillation amplitude. The Reynolds and Keulegan–Carpenter numbers are fixed to *Re*=*U*_{max}*D*/*ν*=100 and *KC*=*U*_{max}/*fD*=5, respectively, where *U*_{max} is the maximum velocity at the cylinder and *f* is the characteristic frequency of oscillation. The computations are run with a constant value of the time step Δ*t*=0.0013*T*, where *T*=*f*^{−1} is the period of oscillation. The simulations were run on the cx1 cluster of Imperial College London^{1} and the associated computational time (in hours) per period of oscillation *T* is indicated in table 2. Figure 5 shows the time evolution of the non-dimensionalized in-line force acting on the oscillating cylinder. The symbols indicate the experimental data [29], whereas the lines are the computational results obtained with the present method. Figure 5*a* shows convergence of the results for decreasing values of the shell thickness at a fixed value of the minimum edge length *l*_{e}/*D*=0.025 around the cylinder. It is clear that a sufficiently small shell thickness is necessary in order to have sensible results. This is expected given that the shell thickness modifies the structure shape, and therefore also impacts on the flow conditions and dynamics. Figure 5*b* demonstrates convergence for decreasing values of *l*_{e}/*D* at a fixed value of the shell thickness Δ/*D*=0.0125. Figure 5*b* shows that, for a small value of the shell thickness, the in-line force is less dependent on the mesh resolution around the cylinder.

### (c) Free fall of a sphere

The last validation test case is flow past a sphere, freely falling under the effect of gravity in a fluid at rest. In this configuration, the time evolutions of the structural velocity and forces (both drag and Basset forces) can be computed semi-analytically by solving the equation of motion of the sphere [32], i.e.
5.1
where *m*_{s}=*ρ*_{s}*V* _{s} is the sphere mass and *m*_{a}=*C*_{m}*ρ*_{f}*V* _{s} is the added mass arising from the relative acceleration of the sphere with respect to the fluid. The added mass coefficient of the sphere is set at , which is the theoretical value in inviscid flows [33]. The first term on the right-hand side of equation (5.1) is the difference between gravity and buoyancy forces. The second term represents the drag force due to pressure and shear effects acting on the sphere, *A* being the sphere cross-sectional area and *C*_{D} the drag coefficient. For non-creeping flows, the latter is given by the following empirical formula [34]:
5.2
The asymptotic value of the drag force is such that it balances gravity and buoyancy forces. The last term in equation (5.1) is the Basset force, which accounts for the time history of the sphere velocity [35,36]. In this study, equation (5.1) is solved using a fourth-order Runge–Kutta method and yields the time evolution of the sphere velocity. This semi-analytical solution is used to validate the results obtained with the IS method. The sphere of diameter *D* is dropped in a computational domain of size 10*D*×75*D*×10*D*. The minimum and maximum edge lengths of the fluid mesh are, respectively, *l*_{e}=0.025*D* (at the fluid–solid interface) and *L*_{e}=5*D* (at the domain outer boundaries). The fluid mesh is regenerated at a period of six time steps, in order to ensure that the fine resolution follows the sphere as it falls. The Courant number is fixed at *CFL*=0.2. The time step is also kept smaller than Δ*t*=10^{−2}, in order to avoid large initial values.

Figure 6 shows the time evolution of the sphere velocity for three values of the Reynolds number, which is computed as *Re*=2*u*_{t}*R*/*ν* (*u*_{t} being the terminal velocity). Figure 6 compare the numerical results (symbols) with the semi-analytical results (continuous lines), for two values of the solid-to-fluid density ratio: *ρ*_{s}/*ρ*_{f}=1.5 and *ρ*_{s}/*ρ*_{f}=4. It is shown that the numerical results agree well with the semi-analytical results. The relative error between the terminal velocity computed by the coupled models, *u*^{n}_{t}, and the semi-analytical value *u*^{a}_{t} is defined as
5.3
and the values in per cent are as follows: for *ρ*_{s}/*ρ*_{f}=1.5, (*Re*=7), (*Re*=20), (*Re*=100); whereas for *ρ*_{s}/*ρ*_{f}=4, (*Re*=7), (*Re*=20), (*Re*=100). The same test case was simulated with the original immersed-body method [11] and yielded an overestimation of the terminal velocity (compared with the semi-analytical solution). The present IS method underestimates the theoretical value of the terminal velocity, especially at large Reynolds numbers. This difference could be explained by the way the forcing is imposed in the momentum equations. In the immersed-body method, the volumetric force is applied across the entire sphere. The numerical results show a uniform fluid velocity across the sphere (on the fluid mesh), with a slight velocity overshoot upstream of the sphere and a slight velocity undershoot downstream of the sphere. By contrast, in the IS method, the volumetric force is solely applied in the thin shell surrounding the sphere. The numerical results show a uniform fluid velocity in the thin shell and a slight underestimation of the velocity inside the sphere. In the coupling procedure, this will in turn affect the sphere velocity in the structural-dynamics model. Further investigation is needed in order to understand these differences. However, despite the underestimation in the structural velocity, the discrete forcing term of the IS method seems to better represent the flow physics around the sphere.

## 6. Conclusion

This paper presents a novel numerical method for computing FSIs for a wide range of applications. It consists of immersing the solid structures in an extended fluid domain, and exchanging the coupling forces through a thin shell surrounding the structures. The use of a shell is advantageous (i) to resolve or parametrize the physical phenomena taking place at the fluid–structure interfaces and (ii) to reduce the computational cost particularly when large structures (or a large number of small structures) are investigated. Importantly, the IS method works with both structured and unstructured meshes, non-coinciding fluid and solid meshes, and different discretization schemes for the fluid- and structural-dynamics equations. This versatility is attractive for many engineering disciplines, including offshore wind energy. For example, the study of wave–structure interactions requires high-order representations of the velocity and pressure fields, in order to accurately model wave propagation over several wavelengths. This study shows that the method provides accurate results on a series of benchmarked test cases. The present analysis is limited to relatively low Reynolds numbers because of the absence of a turbulence model and the high computational cost associated with the present simulations. The simulation cost is due to the fact that none of the structural-dynamics models and the coupling algorithm are currently parallelized. However, the use of a thin shell surrounding the structure is very attractive to tackle high Reynolds numbers flows, for which boundary-layer resolution is critical. Future research will focus on improving the representation of physical phenomena in the IS to tackle the aero- and hydrodynamics of offshore wind turbines. In this context, several challenges are foreseen. As this study showed, the shell has to be sufficiently thin to get accurate numerical results. For slender structures, such as wind turbine blades, or complex geometries in general, having a very thin shell around the structure can be numerically challenging. Additionally, the implementation of a wall-layer model to reproduce the dynamics of high Reynolds number flows near the fluid–structure interface is non-trivial. The case of wind turbines is particularly challenging, owing to their high rotational speeds, and the complex wakes and tip vortices induced by the rotor. Hybrid subgrid-scale models will be developed for this purpose. Finally, moving grids will be considered to further reduce the computational cost associated with rotating turbine blades. Despite these challenges, the IS method is promising to improve the representation of the flow physics near complex fluid–structure interfaces.

## Funding statement

A.V. is supported by the European Union Seventh Framework Programme (FP7/2007-2013) under a Marie Curie Career Integration Grant (grant agreement PCIG13-GA-2013-618159). A.V. also acknowledges support from Imperial College London and its High Performance Computing Service, as well as Université Libre de Bruxelles. C.C.P. acknowledges support from the UK Engineering and Physical Sciences Research Council (Programme grant EP/K003976/1). The content of this paper reflects only the authors’ views and not those of the European Commission.

## Footnotes

One contribution of 17 to a theme issue ‘New perspectives in offshore wind energy’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.