A NURBS-based finite element formulation for incompressible fluid dynamics and fluid-structure interaction with rigid bodies

A numerical investigation is performed here using a NURBS-based finite element formulation applied to classical Computational Fluid Dynamics (CFD) and Fluid-Structure Interaction (FSI) problems. Model capabilities related to refinement techniques are analyzed using a finite element formulation with NURBS (non uniform rational B-splines) basis functions, where B-splines and low-order Lagrangian elements can be considered as particular cases. An explicit two-step Taylor-Galerkin model is utilized for discretization of the fundamental flow equations and turbulence is considered using Large Eddy Simulation (LES) and the Smagorinsky’s sub-grid scale model. FSI is considered using an ALE kinematic formulation and a conservative partitioned coupling scheme with rigid body approach for large rotations is adopted. CFD and FSI applications are analyzed to evaluate the accuracy associated with the different refinement procedures utilized. Results show that high order basis functions with appropriate refinement and non-uniform parameterization lead to better predictions, compared with low-order Lagrangian models.


INTRODUCTION
Although the use of finite volume models is still a common practice in the field of Computational Fluid Dynamics (CFD), the Finite Element Method (FEM) has gained some popularity in the last decades with significant advances observed in computers technology (see, for instance, Zienkiewicz et al., 2013;Reddy and Gartling, 2010). In the FEM context, the flow domain and the fundamental flow equations are spatially discretized using isoparametric finite elements with Lagrangian basis functions, where linear basis with C 0 -continuity are usually adopted. However, the discretization procedure may lead to approximation errors, which may be significant depending on the geometrical characteristics of the physical model to be investigated. In addition, highly nonlinear and small-scale problems, especially turbulent flows, demand for high order discretization in both spatial and time domain. In order to overcome these drawbacks some improvements have been proposed to the finite element formulation, such as the use of NURBS (non uniform rational B-splines) basis functions, which are extensively utilized in Computational Aided Design (CAD) (see Piegl and Tiller, 1997). With this improvement, many possibilities with respect to refinement procedures may be conceived, in spite of some shortcomings that are also observed.
In the field of Engineering Design, computational reproduction of physical models has been traditionally performed using numerical tools based on CAD technologies, where NURBS are very popular. NURBS can exactly reproduce all conic sections and present convenient mathematical properties, such as C p-1 -continuity order for p-degree basis functions, convex hull and variation diminishing properties (see, for instance, Piegl and Tiller, 1997). Recent advances have extended NURBS formulation by using T-splines (Sederberg et al., 2003), which permitted local refinement and compatibility of adjacent patches efficiently. By using NURBS functions in a finite element formulation, pre-processing and analysis procedures become unified, considering that the same numerical tools are employed.
However, when NURBS basis functions are utilized in finite element modeling, the numerical scheme must be reformulated in order to take into account a new framework for spatial interpolation. Basic concepts associated with control points and splines are introduced considering the formulation proposed in the automotive industry by De Casteljau and Bezier (see Townsend, 2014 for additional details), who utilized Bernstein polynomials and control points for manipulating geometric forms mathematically. De Casteljau observed that a curve could be accurately represented by manipulating some points around the curve and not along its length. By moving these points of influence (i.e., the control points), he noticed that the curve was modified such as moving weights in a boat builder spline. Later, the Bezier method was developed using similar concepts (see Bezier, 1972), which was surmounted by a recursive algorithm proposed by de Boor, where B-spline functions conceived by Schoenberg (1946) were adopted.
Other concepts were introduced with the B-spline functions, such as piecewise polynomials and knot vectors, which define the local support where a function is not null within the parametric space, considering that the parametric space is decomposed into knot spans by using breakpoints (knots). B-splines provide great flexibility for reproduction of geometries since the action of the control points is localized. In addition, B-splines can be seen as a generalization of Casteljau's algorithm, including the Bezier method as a special case. Further improvements were obtained from the aerospace industry, where the existing methods were unable to handle with designing and assembling of aircraft components. NURBS were then conceived considering rational B-splines and knot vectors with non-uniform knot spans. By using control points with variable weighting, it was possible to draw conics exactly as well as reproduce complex curves and surfaces accurately (see Piegl and Tiller, 1997). NURBS-based finite element formulations have been widely applied to a variety of engineering problems, from solid mechanics (Espath et al., 2014;Espath et al., 2015) to material science (Gomez et al., 2008).
A NURBS finite element formulation for CFD applications presents significant improvements over the classical Lagrangian finite element formulation. Complex flow phenomena, such as boundary layer and turbulent flows, separation, and reattachment are better reproduced considering that different flow regions can be discretized using distinct combinations of degree and continuity order associated with the interpolation functions. These basis functions are always smooth (C ∞ ) within the knot span and are C p-m continuous at the knots, where p is the polynomial degree and m is the knot multiplicity. Nevertheless, NURBS basis functions are not interpolatory in general since the control mesh defined by the control points does not conform to the actual geometry of the physical model. Notice that flow variables and coordinates are defined at the control points. The flow spatial domain may be decomposed into patches, depending on its geometric complexity, where independent parametric spaces are adopted considering that they must be compatible on patch interfaces. Every patch is also decomposed into knot spans, which define the element concept in a NURBS-based finite element formulation. Most of the classical flow problems can be simulated using a single patch, although multiple patches can also be utilized in order to obtain a more efficient discretization. Finite element Latin American Journal of Solids and Structures, 2020, 17(1), e242 3/37 formulations with NURBS basis functions applied to CFD and FSI problems may be found in Bazilevs et al. (2006), Gomez et al. (2010), Nielsen et al. (2011), Akkerman et al. (2011), and Kadapa et al. (2015). Unlike the standard finite element refinement techniques, B-spline refinement can control element size as well as degree and continuity order of the basis functions. The basic refinement techniques associated with B-splines are knot insertion and degree elevation, which do not change the physical model geometrically or parametrically. By using knot insertion, the solution space is enriched with additional elements, control points, and basis functions, but without changing the geometry. Knot insertion is equivalent to h-refinement adopted in the finite element modeling only when knots are inserted with multiplicity m = p in order to maintain the basis functions with C 0 continuity. In general, knots can be inserted with multiplicity m = 1, maintaining the original continuity of the basis functions, independently of the polynomial degree. This aspect cannot be replicated with a standard finite element formulation.
When order elevation is performed, knot multiplicity of the existing knots is also increased in order to maintain the original continuity order of the basis functions along the element boundaries and no new knots are created. Notice that the classical p-refinement utilized in finite element modeling must be initially applied over basis functions with C 0 -continuity, while p-refinement for B-spline discretizations can be adopted using any continuity order. This aspect cannot be replicated with a standard finite element formulation. An exclusive technique for B-spline refinement is obtained by using degree elevation followed by knot insertion, which is called k-refinement. One can notice that a NURBS-based finite element formulation becomes highly flexible with respect to continuity order of the basis functions and other refinement methods can be devised. Additional information on refinement methods for finite element formulations with NURBS basis functions are found in Cottrell et al. (2009).
In this work, a NURBS-based finite element formulation is proposed considering the explicit two-step Taylor-Galerkin model, where spatial discretization is carried out taking into account B-splines and NURBS basis functions. The fundamental flow equations are the Navier-Stokes equations and the mass conservation equation, which is described according to the pseudo-compressibility hypothesis for incompressible flows and Newtonian fluids under isothermal conditions. In the present model, idealized 2D turbulent flows are simulated using a LES-type (Large Eddy Simulation) approach, where the Smagorinsky's model is adopted for sub-grid scale modeling. Fluid-structure interaction (FSI) problems are reproduced using a conservative partitioned coupling model and a rigid body approach for large rotations. Finally, classical CFD and FSI applications are analyzed in order to validate the present methodology, where different refinement procedures are adopted and investigated.

Flow analysis
In the present model, the flow analysis is performed considering the following assumptions: fluid particles are described constitutively according to the Newtonian model for viscous fluids, and the flow is two-dimensional and restricted to the incompressible regime under isothermal conditions. In addition, fluid body forces are neglected and a mixed approach is adopted for variables definition, where the pressure field is explicitly evaluated using the pseudo-compressibility hypothesis (see Chorin, 1967). An Arbitrary Lagrangian-Eulerian (ALE) formulation is adopted for the kinematical description of fluid motions when fluid-structure interaction is considered, and turbulence flows are approximately reproduced using Large Eddy Simulation (LES) and the Smagorinsky's sub-grid scale model (Smagorinsky, 1963). Thus, the system of fundamental flow equations may be written as follows (see, for instance, White, 2005): Mass balance equation where vi and wi are the components of the flow velocity vector v and mesh velocity vector w, respectively, which correspond to the i-direction of a rectangular Cartesian coordinates system, where coordinates components are Latin American Journal of Solids and Structures, 2020, 17(1), e242 4/37 denoted by xi, p is the thermodynamic pressure, ρ is the fluid specific mass and c is the sound speed in the flow field.
Although the sound speed c usually has a strict physical meaning, this parameter is considered here as a numerical parameter in order to satisfy the incompressible flow condition (divergence-free velocity field with: where μ and λ are the dynamic and volumetric viscosities of the fluid and δij are the components of the Kroenecker's delta (δij = 1 for i = j; δij = 0 for i ≠ j). Idealized turbulence flows are simulated here using a LES-type approach, where the fundamental flow equations are submitted to spatial filtering in which the flow field is decomposed into large and small-scale components. Large scales are solved directly utilizing the filtered equations, while scales smaller than the mesh resolution are defined by the sub-grid stress tensor τ SGS , which must be modeled considering a turbulence closure model. Turbulence models are utilized in order to represent the small-scale effects over the large scales. The sub-grid stress tensor components SGS ij  are usually approximated according to the Boussinesq assumption, that is: where ij S are components of the filtered strain rate tensor. In the present work, eddy viscosity μt is obtained considering the Smagorinsky's sub-grid scale model, which may be expressed as follows: where CS is the Smagorinsky constant, which must be specified according to the flow characteristics, with values usually ranging from 0.1 to 0.25, and  is the characteristic length of the spatial filter, which may be locally defined as E    for a box filter (see Smagorinsky, 1963 for additional details), where ΩE is the finite element area referring to element E. When laminar flows are investigated, the sub-grid stress tensor components SGS ij  are omitted in Eq. (1).
It is important to notice that the present approach is only an approximation of the actual LES, considering that a turbulent flow is inherently three-dimensional. In order to solve the flow problem, initial conditions on the flow variables vi and p must be specified. In addition, appropriate boundary conditions must also be defined on Tf t  , which may be expressed as: Latin American Journal of Solids and Structures, 2020, 17(1), e242 5/37 where fsi t  (boundary representing the fluid-structure interface), v  (boundary with prescribed velocity i v  ), p  (boundary with prescribed pressure p  ) and   (boundary with prescribed traction t i  ) are complementary subsets of In Eq. (10), nj are components of the unit normal vector n evaluated at a point on boundary   . Notice that wi = 0 for points outside the ALE domain or when motions of the immersed body are not considered in the flow analysis.

Fluid-structure interaction
The fluid-structure coupling is accomplished by enforcing equilibrium and kinematical conditions on the fluid-structure interface fsi t  . The no-slip condition is assumed for viscous fluids, such that the relative velocity of fluid particles on the fluid-structure interface is set to zero. In this case, the equilibrium and compatibility equations may be expressed as follows: where nj and s t i are components of the unit normal vector n and the structure traction vector t s evaluated at a point on  . In addition, continuity conditions are also imposed on fsi t  when moving grids are adopted, which may be given as follows: where xi are components of the mesh position vector x referring to the flow spatial domain. In order to solve the flow problem on moving grids adequately, the geometric conservation law (GCL) must be satisfied (see Thomas and Lombard, 1979 for further information). According to Lesoinne and Farhat (1996), the GCL is satisfied in ALE finite element formulations if the mesh velocity vector w is calculated as: Considering that the trapezoidal form of the Newmark's method is utilized for solving the structural equation of motion and a partitioned algorithm is adopted for fluid-structure coupling, a conservative algorithm where the GCL is satisfied without violating the compatibility conditions can be obtained if the following equations are employed (see Lesoinne and Farhat, 1998): where:  Figure 1 shows the coupling scheme adopted here for fluid-structure interaction problems, where one can observe that the flow and structure analyses are performed sequentially. In addition, it is also observed that fluid and structure are displaced in time by a half time step Δt/2. The mesh motion must be arbitrarily defined when ALE formulation is adopted for the kinematical description of the fluid flow. In this work, a mesh motion scheme utilized by Braun and Awruch (2009) is employed, where the mesh velocity field is obtained as follows: where NS is the number of mesh points located on the boundaries of the ALE domain and aij are influence coefficients defined with mesh points i and j, considering that i are inner mesh points and j are mesh points belonging to the boundaries of the ALE domain, where wi = 0. The influence coefficients are obtained with: where dij is the Euclidian distance between the mesh points i and j and exponent n is a user-defined parameter, which is chosen according to the amplitude of the immersed body displacements. In this work, all FSI applications were analyzed using n = 4. The structural response is obtained here considering the equation of motion for two-dimensional rigid bodies with elastic constraints and viscous dampers, which may be expressed as: where M s is the mass matrix, C s is the damping matrix, K s is the stiffness matrix and U s c  , U s c  and U s c are the acceleration, velocity and displacement vectors evaluated at the center of mass of the structure, which are written using three degrees of freedom, with two translational displacements and one rotational displacement. The load vector applied at the center of mass of the structure is denoted by Q s c , with two force components and one moment component. Notice that s and c indicate that quantity refers to structure and center of mass, respectively.
The structural motion on the fluid-structure interface is obtained from the structural response evaluated at the center mass by considering the following kinematic relations (see Figure 2): where int u s  and int u s  are the velocity and acceleration vectors evaluated at a point on the fluid-structure interface, rc,int is the position vector referring to a point on the fluid-structure interface with respect to the center of mass of the structure, ω s is the vector of angular velocity and α s is the vector of angular acceleration.
The kinematic relations above may also be expressed using matrix form, that is: where: Latin American Journal of Solids and Structures, 2020, 17(1), e242 8/37 The fluid forces acting on the structure are obtained considering the equilibrium condition imposed on the fluidstructure interface. This load must be transferred to the center of mass of the structure in order to solve the equation of motion for rigid bodies. This procedure may be carried out using the following expressions (see Figure 3): where f int t is the fluid traction vector, which is evaluated using Eq. (10) for a point located at the fluid-structure interface, where the translation matrix L given by Eq. (26) is also evaluated.

Figure 3
Equilibrium conditions on the fluid-structure interface.

NURBS-based finite element formulation
In a NURBS-based finite element formulation for fluid dynamics the spatial domain is initially decomposed into patches according to the geometric complexity of the problem investigated, where the basis functions are defined using individual parametric spaces. In order to define the element concept in the present formulation, every patch is divided into knots spans specified in the different directions of the parametric space by using the knot vectors (see Figure 4). For two-dimensional problems, the following knot vectors may be utilized: where p and q are the polynomial degrees of the basis functions defined over the parametric directions ξ and η, respectively. The number of basis functions associated with the parametric directions ξ and η is defined by n+1 and m+1, respectively, which also defines the corresponding number of control points, and the number of elements is determined by the number of non-zero knot spans. The NURBS basis functions for two-dimensional applications are given by: Latin American Journal of Solids and Structures, 2020, 17(1), e242 9/37 where the subscripts i and j indicate the position of the control point in the index space and the superscripts p and q denote the polynomial degree of the basis functions. The weight term wi,j is related to the weight associated with the control point defined by the subindices i and j. The B-spline basis functions Ni are evaluated here using the Cox-de Boor recursive formulation (Cox, 1972;De Boor, 1972), which may be expressed as: where p is the polynomial degree of the basis function N(ξ) and i is the knot index. Notice that Eq. (34) is straightforwardly extended to basis functions associated with the parametric direction η.

Figure 4
Definitions on spatial discretization using a NURBS-based finite element formulation.
Considering that n+1 and m+1 denote the number of basis functions related to the parametric directions ξ and η, respectively, and the respective polynomial degrees are defined by p and q, element e is identified by determining the indices at which the corresponding non-zero knot span begins in the index space, that is: where p+1 ≤ i ≤ n and q+1 ≤ j ≤ m. The total number of elements in which the spatial field is discretized in the parametric domain is defined as: Latin American Journal of Solids and Structures, 2020, 17(1), e242 10/37 Finally, geometry and flow variables are discretized using the following NURBS approximations: where Ri is the NURBS basis function related to control point i, which is defined as a function of the parametric coordinates (ξ,η), and ncp is the number of global control points (or basis functions). The control point variables and geometry are specified by vi, δvi, p, δp and xi. Notice that the sum operations indicated above are performed over the total number of basis functions available (ncp), considering that basis functions support is highly localized.

The explicit two-step Taylor-Galerkin model
In the present model, the fundamental flow equations are discretized using the explicit two-step Taylor-Galerkin model, where second-order Taylor series and the Bubnov-Galerkin method are adopted for time and space discretization considering a NURBS-based finite element framework.
Time discretization is initially applied to the fundamental flow equations considering a second-order Taylor series expansion and a two-step time increment scheme, which leads to the following formulation: where: where Δt is the time increment. Notice that the flow variables must be evaluated at n+1/2 in order to obtain the velocity and pressure increments defined by Eqs. (42) and (43). The velocity and pressure fields at n+1/2 are obtained from: Latin American Journal of Solids and Structures, 2020, 17(1), e242 11/37 1 2 2 2 n n j n The velocity field at n+1/2 must be corrected considering that a pressure increment term was omitted in Eq. (44). The velocity correction is performed using the following equation: The weighted residual principle is adopted here in order to minimize approximation errors associated with geometry and flow variables, which may be written here as follows: where Ωe is spatial domain referring to element e and nel is the number of elements in the finite element mesh. Rv and Rp are the residual vectors referring to the momentum and mass balance equations, respectively, which are obtained by approximating the flow variables and geometry with the finite element interpolations given by Eqs. (37), (38) and (39). By using the Bubnov-Galerkin method, the weight functions are defined here using Wv = δv and Wp = δp, which are also approximated using velocity and pressure variations given by Eqs. (38) and (39). In order to reduce continuity requirements over the basis functions, the weak form of the flow formulation is adopted using the Green-Gauss theorem. An algebraic system of equations is obtained for the two-step Taylor-Galerkin model utilized in this work, where the flow variables are first evaluated at n+1/2 with the following expressions: The velocity field is then corrected with: Finally, the flow variables at n+1 are obtained using: where M is the fluid mass matrix, G is the gradient matrix, A is the advection matrix, D is the diffusion matrix and B is the balance diffusivity tensor. Detailed information on the two-step Taylor-Galerkin scheme may be found in Kawahara and Hirano (1983). The element matrices are numerically evaluated using Gaussian quadrature, where three different spaces are considered: the physical space Ωe, which is defined by the vector of Cartesian coordinates x = (x1, x2); the parametric space ˆe  , which is defined by the vector of parametric coordinates ξ = (ξ,η) and the quadrature space e   , which is defined by the vector of quadrature space coordinates . The mapping from physical to quadrature space is performed using the Jacobian transformation matrix as follows: The parametric and quadrature spaces are related according to the following equations: where the parametric space referring to element e is defined by The time increment Δt is locally determined using the Courant's stability condition, that is: where Δte is the time increment referring to element e, α is the safety coefficient (0 ≤ α ≤ 1), Δxe is the characteristic length of element e, c is the sound speed, and ve is the flow characteristic speed. Although a multi time step formulation may be adopted for time integration, the smaller Δte is utilized in this work throughout the finite element mesh.

Fluid-structure coupling model
In the present work, a conservative partitioned model is adopted for fluid-structure coupling. Considering a fluid element ΩE localized on the fluid-structure interface, the discretized momentum equation referring to element E may be written as: This equation is rearranged in order to define mesh points of element E belonging to the fluid-structure interface (I) and mesh points located within the fluid flow domain (F). A matrix form of the rearranged system may be expressed as: II  IF  I  II  II  IF  IF  I  I  I  I  f f Latin American Journal of Solids and Structures, 2020, 17(1), e242 13/37 One can notice that only the first line of the system of equations above is relevant for the fluid-structure problem. In addition, the compatibility conditions are also expressed in matrix form considering an element E belonging to the fluid-structure interface as follows: where: Notice that N is the number of mesh points referring to element E and the matrices Li and L'i associated with mesh points out of the fluid-structure interface are set to zero.
The load vector at the center of mass of the structure is obtained considering the equilibrium equation, which may be defined as: where FI is the fluid force vector corresponding to the first row of the system of equations given by Eq. (59). By substituting Eq. (59) into Eq. (63), one obtains the following expression: Finally, the compatibility conditions given by Eqs. (60) and (61) are imposed onto Eq. (64), which leads to:

Q T M TU T M v T A T + T D T + T M T U T A + T D v T G p
The equivalent equation of motion referring to the structure domain is obtained substituting Eq. (65) into Eq. (21), that is: where the equivalent mass matrix and the equivalent damping matrix are given by:

Q = T M v + T A + T D v T G p
Notice that all the fluid elements in contact with the fluid-structure interface (NELI) are included in the sum operation indicated by Eqs. (67), (68) and (69). Matrix C is nonsymmetric and nonlinear due to the advection f II A and translation T matrices. Considering that the time step adopted for time integration is generally small, the nonlinear matrix is linearized using ω3 from the last time step. The structural equation of motion defined by Eq. (66) is solved here using the implicit Newmark's method (see Bathe, 1996 for detailed information).

NUMERICAL APPLICATIONS
In this chapter, typical CFD and FSI problems are solved using the algorithms presented in the previous chapters. It's noteworthy that the preprocessing and analysis tools employed in this work were developed by the authors and implemented using Fortran programming language. For the post-processing stage, the commercial software Tecplot 9 was employed.

Wall-driven cavity flow
The wall-driven cavity flow problem is analyzed here in order to validate the numerical model proposed in this work for flow simulation with B-spline basis functions. In addition, the influence of mesh refinement over numerical results is also investigated. Incompressible flow regime is assumed and two-dimensional flows are considered, which are characterized by the following Reynolds numbers: Re = 10 3 and Re =10 4 . Notice that a LES-type approach is employed for Re = 10 4 , where the Smagorinsky's constant is set to CS = 0.15 (see Eq. 6). This classical problem consists of a square cavity with unit length walls (L = 1 m), where the side and bottom walls are submitted to no-slip boundary conditions, while the top wall slides laterally with a constant velocity U = U0.
A special discretization procedure is proposed in the present analysis considering that uniform knot vectors are adopted in both parametric directions and the control meshes are defined arbitrarily according to specifications presented in Table 1. The distribution of control points over the physical space is determined considering geometric progression from the smallest distance between adjacent control points, which is always localized next to the walls. The distance between the first two control points next to the walls is defined here as h and the following values are used for all the meshes shown in Table 1: h = 0.001L, h = 0.0025L, h = 0.005L, and h = 0.01L. Notice that unlike the discretization procedure utilized in this work, a control mesh is usually defined considering an original coarse mesh, which is successively refined using some of the existing refinement methods for B-spline and NURBS basis functions. The time increment utilized in the time integration of the flow equations is obtained from Eq. (57) with α = 0.3. Table 1 Meshes utilized in the cavity flow analysis.

Function degree (p = q)
1000 60 x 60 1, 2, 3 80 x 80 1, 2, 3 100 x 100 1, 2, 3 10000 80 x 80 1, 2, 3 100 x 100 1, 2, 3 120 x 120 1, 2, 3 The influence of the refinement level over the present numerical predictions is initially analyzed using three control mesh configurations (60x60, 80x80 and 100x100) with basis functions degree defined by p = q = 1 and refinement parameters h = 0.001L and h = 0.0025L. The flow is characterized with a Reynolds number Re = 10 3 . Notice that these mesh configurations reproduce classical finite element grids with Lagrangian basis functions and C 0 -continuity order. Figure 5 shows velocity profiles along the horizontal and vertical lines of the cavity referring to results obtained here and the Latin American Journal of Solids and Structures, 2020, 17(1), e242 15/37 corresponding predictions obtained by Ghia et al. (1982), who adopted a second-order finite difference model with 257x257 grid points.
One can see that a very good agreement is obtained with respect to the reference predictions by using h = 0.0025L and any of the mesh configurations proposed in the present application. On the other hand, when a mesh parameter h = 0.001L is adopted, differences are observed in the central region of the cavity, where predictions are improved as the number of control points of the mesh configuration is increased. Nevertheless, a very good agreement is also obtained near the wall regions for all mesh configurations utilized here. One can notice that excessive refinement near the wall regions (h = 0.001L) may lead to insufficient refinement in the central region of the cavity when coarse meshes are utilized.
In order to analyze the influence of mesh parameter h in the present study, two control mesh configurations (80x80 and 100x100) are adopted using basis functions with p = q = 1. Results referring to velocity profiles along the horizontal and vertical centerlines of the cavity are presented in Figure 6, where results obtained here are compared with numerical predictions obtained by Ghia et al. (1982) for Re = 10 3 .
The effect of mesh parameter h over the numerical predictions is clearly demonstrated. It is observed that insufficient results are obtained when h = 0.01L is adopted, independent of the mesh configuration utilized. One can also notice that the effect of using a more refined mesh configuration over the numerical results is small in this case.
The influence of basis functions over the numerical results is investigated here using an intermediate control mesh configuration (80x80) and two intermediate refinement parameters, h = 0.0025L and h = 0.005L. The present results are shown in Figure 7, which are compared with predictions obtained by Ghia et al. (1982).  The present results demonstrate that the control mesh configuration with mesh parameter h = 0.001L is clearly insufficient to reproduce the reference results, independent of the basis function degree utilized. In this sense, one can observe that by increasing the basis function degree, accuracy is gradually reduced, especially in the central region of the cavity. On the other hand, when the mesh parameter h value is increased to h = 0.0025L, the present results are improved. In addition, the influence of the basis function degree is now inverted, leading to better predictions when the function degree is increased, although the improvements are not significant in this case. Excessive refinement near the wall regions led again to insufficient refinement in the central region of the cavity when a coarse mesh was utilized.
An overall evaluation of the present results can be performed considering Figures 9 and 10, where mesh quality is determined using the Euclidean norm of vectors of velocity variation. These vectors are defined taking into account sample points proposed by Ghia et al. (1982) to define the velocity profiles along the center lines of the cavity (see Ghia et al. (1982) for detailed information). The Euclidean norms referring to velocity vector components are defined here as follows: where v1 and v2 are velocity vectors obtained from sample points located on the center lines of the cavity according to positions defined by Ghia et al. 1982 and n is the number of sample points. Notice that each velocity vector corresponds to a specific refinement parameter hj (j = 1,…,4), with h1 = 0.001L, h2 = 0.0025L, h3 = 0.005L, h4 = 0.01L. It is observed that mesh configurations with an insufficient number of control points lead to inaccurate predictions, where the use of smaller mesh parameters h and basis functions with higher degree deteriorates the numerical results due to inadequate spatial discretization. This effect may also be associated with the refinement procedure proposed for the present application, where the control points were distributed arbitrarily. For a control mesh configuration with an appropriate number of control points, mesh parameter and basis function degree, one can notice that v1 and v2 norms behave asymptotically with respect to mesh parameter h, converging usually to a determined error level as h → 0, which may also be null if the spatial discretization is optimized. In addition, basis functions with low continuity order lead to slightly better results when the present refinement procedure is adopted. The same effects formerly reported with respect to numerical predictions obtained with insufficient mesh configurations are observed here and amplified. Considering that a flow with higher Reynolds number is simulated, higher levels of spatial discretization are required. The asymptotic behavior of v1 and v2 norms as functions of the mesh parameter h is only identified for the control mesh configuration 120x120. Results indicate that the mesh configurations utilized in this case are not sufficiently optimized. This aspect may be associated with the spatial discretization methodology proposed in this work, where the control points are arbitrarily distributed over the computational domain. For a high Reynolds number flow, basis functions with higher degrees seem to lead to better results when compared with predictions obtained with low-order basis functions, which are traditionally adopted by lagrangian-based finite element formulations.
The cavity flow field is shown in Figure 11, where pressure field and streamlines obtained with the present formulation are presented. One can notice that the present predictions agree very well with the classic results obtained by Ghia et al. (1982), with primary, secondary and tertiary recirculation zones reproduced correctly according to the Reynolds number adopted.

Figure 11
Pressure field and streamlines.

Flow around circular cylinder
The flow around circular cylinder is analyzed here using a wide range of Reynolds numbers and different refinement procedures, where a viscous fluid under incompressible flow regime is utilized. Although it is well known that the flow regime around circular cylinders is turbulent for Re > 300 (see, for instance, Lienhard, 1966), a two-dimensional approach is adopted here as a first approximation to the actual problem. NURBS basis functions are adopted for spatial discretization.
The numerical model is initially validated using a computational domain constituted by a single patch and quadratic basis functions in both parametric directions (p = q = 2), where uniform knot vectors are employed. It is important to notice that quadratic NURBS functions must be employed along the angular direction (p), at least, in order to obtain a circular curve exactly (see Piegl and Tiller, 1997). The flow field is characterized using inflow velocity v∞ = 10 m/s and the following Reynolds numbers (Re): 10, 20, 30, 40, 50, 100, 300, 500, 700 and 1000, where control mesh configurations of 100x70, 120x90 and 168x120 control points are employed for 10 ≤ Re ≤ 50, 100 ≤ Re ≤ 300 and 500 ≤ Re ≤ 1000, respectively. Control points are distributed along the radial direction using geometric progression, considering that the smaller distance between two consecutive control points corresponds to the first pair of control points next to the cylinder surface. This smaller value is defined here as the mean distance between consecutive control points along the angular direction, which are only related to control points associated with the cylinder surface. Results referring to drag coefficient (CD) and Strouhal number (St) are presented in Table 2, where predictions obtained in this work are compared with results obtained numerically by Henderson (1997). Drag (CD) and lift (CL) coefficients, as well as the Strouhal number (St), are obtained here using the following expressions: where V is the undisturbed flow speed, D is the cylinder diameter, f is the vortex shedding frequency obtained from time histories of lift coefficient, f 1 t and f 2 t are the horizontal and vertical components of the fluid traction vector t f , which is evaluated using Eq. (10) for a point located at the fluid-structure interface, Γe is the boundary of element e in contact with the cylinder surface and nel is the number of fluid elements on the cylinder surface. An excellent agreement with the reference results can be observed. Flow characteristics may be observed in Figure 13, where streamlines and pressure fields are presented. One can notice that a recirculation zone is formed behind the cylinder for Re = 10, 20 and 30, while the von Karman vortex street with alternate vortex shedding is obtained for Re = 100, 500 and 1000, which demonstrates that the formulation proposed in this work can reproduce complex flow phenomena accurately. The distribution of pressure coefficient over the cylinder perimeter is shown in Figure 14 for some of the Reynolds numbers investigated here.   Table 3 Drag coefficient (CD) and geometric characteristics of the recirculation zone behind the circular cylinder for Re = 40.

Parameter
Mesh 1 (p = q = 2) Mesh 2 (p = 2; q = 1) Reference CD Results referring to drag coefficient (CD) and geometric characteristics of the recirculation zone obtained behind the circular cylinder are summarized in Table 3, which are defined considering the geometric parameters presented in Figure 15. The present results are compared with numerical predictions obtained by Wanderley and Levi (2002) using a finite difference model, where a good agreement is obtained using both the basis functions proposed here, although results obtained with p = q = 2 are slightly better. It is important to notice that a mesh configuration with p = 2 and q = 1 lead to significant reductions in terms of computational efforts when compared with the processing time spent by the mesh configuration with p = q = 2, considering that full Gauss quadrature is employed here for numerical evaluation of finite element quantities, such as element matrices and vectors.
Additional simulations are performed in order to evaluate the influence of refinement parameter h (smallest distance between adjacent control points), uniformity of the knot vectors and basis degree over the numerical results. For Re = 40, three refinement levels are utilized, which are defined by h = 2.09x10 -2 m, h = 3.14x10 -2 m and h = 4.71x10 -2 m. For Re = 100, the following refinement parameters are considered: h = 1.75x10 -2 m, h = 2.62x10 -2 m and h = 3.93x10 -2 m. Finally, h = 1.25x10 -2 m, h = 1.87x10 -2 m and h = 2.80x10 -2 m are utilized for Re = 1000. Configurations with linear and quadratic basis functions in both parametric directions are also considered.
Results obtained with the present formulation are shown in Figure 16, where drag coefficient and the Strouhal number are considered. It is observed that spatial discretizations with non-uniform knot vectors generally lead to better convergence when compared with predictions obtained using uniform knot vectors. In addition, degree elevation also improved convergence for flow conditions defined by Re = 40 and Re = 100. However, one can notice that similar trends are not observed for Re = 1000.

Backward-facing step flow
The backward-facing step flow is investigated in the present example using multi-patch refinement, where the computational domain is decomposed into parametrically independent sub-domains. This procedure is useful for applications with complex geometry or complex flow fields. Six different flow conditions are analyzed, which are characterized by the following Reynolds numbers: Re = 100, 200, 400, 600, 800 and 1000.
Latin American Journal of Solids and Structures, 2020, 17(1), e242 25/37 Five patches are utilized here, where mesh control configurations were defined considering a uniform distribution of control points and different combinations of basis function degrees along the respective parametric direction (p = q = 1; p = q = 2; p = 1 and q = 2 for patches 1 and 4 only, with p = q = 1 for the remaining patches). The control points are distributed over the patches as follows: • for Re = 100: patches 1 and 4 -60x20, patches 2 and 3 -200x20 and patch 5 -20x20; • for Re = 200: patches 1 and 4 -100x20, patches 2 and 3 -180x20 and patch 5 -20x20; • for Re = 400: patches 1 and 4 -200x24, patches 2 and 3 -220x24 and patch 5 -24x24; • for Re = 600: patches 1 and 4 -274x24, patches 2 and 3 -200x24 and patch 5 -24x24; • for Re = 800 and 1000: patches 1 and 4 -400x30, patches 2 and 3 -224x30 and patch 5 -30x30; The computational domain utilized in the present investigation is shown in Figure 17, where the patch distribution is identified. The geometric characteristics of the channel are defined as follows: h = 1m, s = 0.94 m, xe = 1m, xt = 30 m, L1 = 1 m, L2 = 17 m and L3 = 12 m. No-slip boundary conditions are adopted on the channel walls, while a velocity profile with parabolic distribution is prescribed at the channel entrance. The outflow condition is imposed considering p = 0 at the channel exit. The parabolic inflow velocity is defined using: where max V 1 is the maximum velocity (10 m/s) over the velocity profile and x2 is the vertical coordinate defined with respect to the coordinate system shown in Figure 17. The Reynolds number is calculated according to the expression presented in Armaly et al. (1983): Results obtained with the numerical model proposed in this work are presented in Table 4, where predictions referring to reattachment length in the recirculation region after the facing step are compared with experimental results obtained by Armaly et al. (1983) and numerical results obtained by Williams and Baker (1997) with a two-dimensional finite element model. The reattachment length is normalized with respect to the geometric parameter s (see Figure 17). One can observe that the reattachment lengths obtained here for Reynolds < 400 are close to experimental results obtained by Armaly et al. (1983), while a better agreement is obtained with respect to numerical predictions obtained by Williams and Baker (1997) for Re ≥ 400. It is important to notice that two-dimensional numerical models underestimate the extent of the primary separation region for Reynolds numbers greater than 400 when compared with predictions obtained experimentally. It has been postulated that this disagreement between physical and computational experiments is due to the onset of three-dimensional flow near Re = 400 (Williams and Baker, 1997). Figure 18 shows streamlines obtained numerically for the flow region after the backward facing step. One can notice that the reattachment length increases as the Reynolds number is increased. In addition, it is observed that a secondary recirculation region is formed along the superior wall of the channel for Re ≥ 600. These observations are in agreement with results presented by Williams and Baker (1997  Flow conditions can also be evaluated using the pressure fields obtained here as functions of the Reynolds number, which are presented in Figure 19. For Re = 600, 800 and 1000, these flow fields are obtained considering a time average over the last half of the simulation period. Taking into account that patch interfaces are critical for the present formulation, where basis functions are only C 0 -continuous, it is observed that streamlines ( Figure 18

Lock-in analysis for elastically supported circular cylinder
In the present example, the lock-in phenomenon observed in elastically-mounted circular cylinders submitted to fluid flow is numerically simulated. It is well known that the lock-in phenomenon occurs for a specific range of flow speeds, where synchronization between the mechanical frequency of the structure and the vortex-shedding frequency is obtained. During the lock-in, the amplitude of oscillations is increased, although rarely exceeding half of the cylinder diameter (see, for instance, Simiu and Scanlan, 1996). One can notice that the synchronization frequency is not necessarily the natural frequency of the structure.
The computational domain utilized here is the same as that presented in Figure 12, which is discretized now using a control mesh with 180x120 control points and non-uniform knot vectors. Quadratic NURBS basis functions are adopted in both parametric directions (p = q = 2). This control mesh configuration is obtained from a coarse grid, which is improved by using knot refinement. The grid spacing is gradually reduced towards the cylinder surface, where the smallest distance is Δx = 8.73x10 -2 m. The time step employed in the time integration of the flow equations is Δt = 3.23x10 -2 s (see Eq. 57).
The lock-in phenomenon is investigated using a fixed Reynolds number Re = 150 and different reduced velocities (Vred = V/f.D) within the range 3 ≤ Vred ≤ 8. The cylinder is restricted in angular and horizontal directions and free to vibrate perpendicularly to the flow, which is mechanically described with the following dimensionless form of the structure equation of motion: where uy and Cy are the structure displacement and force coefficient transverse to the flow direction, ξ is the damping ratio and Mred = M/ρD 2 is the reduced mass. In the present investigation, no damping is considered and a constant reduced mass Mred = 2 is adopted. The structural response obtained with the present formulation is shown in Figure 20 for some of the reduced velocities analyzed here. The typical lock-in response can be clearly identified for Vred = 4, while a structural response with relatively low displacement amplitudes are observed when Vred = 8 is considered. In Figure 21, the normalized maximum displacement of the structure is plotted against reduced velocity and compared with other numerical predictions (Ahn and Kallinderis, 2006;Borazjani et al., 2008), where one can see that the lock-in interval can be identified within the reduced velocity range 4 ≤ Vred ≤ 7, considering that the vibration amplitude is significantly reduced outside this range.
The flow field may be evaluated using Figure 22, where instantaneous pressure and vorticity are presented for Vred = 4 and Vred = 8 at t = 80 s. These results are similar to predictions obtained numerically by Borazjani et al. (2008) using a FSI model and the immersed boundary method.

Galloping analysis of an elastically supported square cylinder
Galloping is numerically simulated here considering a square cylinder with an elastic support mounted transverse to the flow direction. The galloping phenomenon is typically found in structures with special cross-section shapes, such as rectangular ad D shapes, which may exhibit large amplitude oscillations in the direction normal to the flow at lower frequencies than those associated with vortex shedding (see Simiu and Scanlan, 1996).
Geometry and boundary conditions utilized in the present application are found in Figure 23, where dimensionless values are indicated. The computational domain is discretized using 180x90 control points, non-uniform knot vectors with linear basis functions in the angular direction and quadratic basis functions in the radial direction. The smallest grid distance is found next to the cylinder surface with Δx* = 2.22 x10 -2 and the time step for time integration is set to Δt* = 1.26 x10 The structural response obtained here is presented in Figure 24 and instantaneous flow conditions are shown in Figure 25, where the pressure and vorticity fields are presented at t* = 301.93 and t* = 310.61. One can observe that the displacement amplitude is relatively large and limited, which characterizes the typical conditions for dynamic instability by Galloping. Table 5 summarizes the present results referring to maximum displacement amplitude and frequency of oscillation, which are compared with predictions obtained from other authors (Dettmer and Peric, 2006;Robertson et al., 2003). Results obtained here demonstrate a very good agreement with other predictions considering flow aspects and structural vibration characteristics.

Flow over elastically supported rectangular cylinder
Fluid-rigid body interaction is analyzed in this example considering a rectangular cylinder immersed in a viscous fluid flow. Two conditions are investigated: (a) the structure is released from an initial angular displacement with the fluid at rest and no structural damping; (b) the structure is submitted to uniform flow, where vertical and angular displacements are induced simultaneously. The computational domains and boundary conditions adopted in the present analyses are shown in Figure 26. In the first investigation, the cylinder is released from an initial angular displacement of θ0 = 5°, considering that fluid and structure are at rest. The structural properties are defined taking into account a torsional frequency 0.266 rad/s k I     , where kθ is the torsional stiffness, and I is the mass moment of inertia. No structural damping is considered in the present analysis and three Reynolds numbers are utilized. The Reynolds number associated with the flow is calculated here using a characteristic length L = b/2 = 1.25 m, where b is the base of the rectangle, and a characteristic flow speed   x 2 0 sin 2 2.90245 10 m/s V T L      , where T is the period of vibration. A control mesh with 161x71 control points is utilized, which is obtained considering knot refinement from an initial coarse mesh. Uniform knot vector in the angular direction and non-uniform knot vector in the radial direction are adopted, considering that linear basis functions are employed in both parametric directions. Results referring to the structural response obtained with the present formulation are shown in Figure 27. One can observe that the angular displacements are gradually reduced over the time due to the action of fluid viscosity since no structural damping is considered here. In addition, notice that flow damping gets more significant as the Reynolds number is reduced. The present predictions are in agreement with results presented by Sarrate et al. (2001), where a finite element model was utilized. In the second investigation proposed here, the rectangular cylinder is analyzed considering the action of a uniform flow (see Figure 26), which leads to vertical and angular displacements in the structure. The dimensionless parameters referring to the structure properties are set as follows: mt = 195.57, ct = 0.0325, kt = 0.7864, I = 105.94, cr = 0.0 and kr = 17.05, where m, I, c and k denote mass, mass moment of inertia, damping and stiffness constants and the subscripts t and r are related to translational and rotational degrees of freedom. The flow field is characterized by a Reynolds number Re = 1000.
A control mesh configuration with 161x71 control points is employed, which is similar to that utilized in the previous study. Comparisons are performed considering mesh configurations with uniform and non-uniform knot vectors. In addition, predictions obtained from a mesh configuration with linear basis functions in both parametric directions are compared with results obtained from a mesh with linear basis functions in the angular direction and quadratic basis functions along the radial direction. The structural responses obtained here are presented in Figure 28, where predictions referring to different mesh conditions are compared. Notice that the structural vibration is anticipated when uniform knot vectors are adopted, although the increase rate and oscillation frequency are the same as those obtained with non-uniform knot vectors. The use of basis with a higher degree also anticipates the onset of structural vibration without modifying the increase rate and frequency of oscillation. When linear basis are utilized in conjunction with non-uniform knot vectors, vertical displacements with greater amplitude are obtained in the early stages of the structural response.
The mesh configuration and the corresponding pressure field and streamlines for dimensionless time instants t* = 439 and t* = 448 are shown in Figure 29, where one can see that the mesh motion scheme adopted in the present work can accommodate the structural motion in the region next to the immersed body without excessive distortion. It is important to notice that using n = 4 in Eq. 20, one obtains a mesh motion such that elements near the immersed body show motion characteristics similar to rigid bodies, maintaining its original geometric aspects. On the other hand, elements located in the intermediate region between the fluid-structure interface and the outer border of the ALE domain present deformation similar to elastic bodies. As exponent n is reduced, this elastic aspect of the element motion is observed in the mesh region near the immersed body.

CONCLUSIONS
In the present work a NURBS-based finite element formulation for incompressible fluid dynamics and fluid-structure interaction with rigid-body dynamics was presented. Model versatility with respect to spatial discretization procedures was demonstrated, where different techniques referring to NURBS discretization, such as multi-patches, knot insertion, and degree elevation, were adopted. In addition, a control mesh discretization scheme was proposed, where control point locations are specified arbitrarily. The present investigation indicates that results obtained with the present model are very sensitive to the mesh configuration utilized. In this sense, if the control mesh distribution is not adequately set, good predictions are not obtained even with high degree basis functions. In this case, the use of functions with higher degree leads to deterioration of the numerical predictions. On the other hand, one can observe that improvements can be obtained by using higher degree basis functions when flows with high Reynolds numbers are investigated, especially for turbulent flows. The use of non-uniform knot vectors usually leads to better results when compared with predictions obtained with uniform knot vectors. This aspect is important in order to capture complex flow phenomena, especially in the boundary layer region. Computational efforts are usually high when the present formulation is applied to turbulent flows, considering that full Gauss quadrature is employed in this work to evaluate element matrices and vector. This drawback may be circumvented by using specialized techniques for numerical integration of NURBS basis functions. In this sense, a reduced integration formulation for NURBS-based finite element models may also be developed. For future works, the present formulation must be extended to include the three-dimensional approach. Potential applications for the present model may be found in the field of Computational Wind Engineering, such as aerodynamic and aeroelastic analyses of long-span bridges, low and high-rise buildings.