Abstract
A numerical investigation is performed here using a NURBSbased finite element formulation applied to classical Computational Fluid Dynamics (CFD) and FluidStructure Interaction (FSI) problems. Model capabilities related to refinement techniques are analyzed using a finite element formulation with NURBS (non uniform rational Bsplines) basis functions, where Bsplines and loworder Lagrangian elements can be considered as particular cases. An explicit twostep TaylorGalerkin model is utilized for discretization of the fundamental flow equations and turbulence is considered using Large Eddy Simulation (LES) and the Smagorinsky’s subgrid 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 nonuniform parameterization lead to better predictions, compared with loworder Lagrangian models.
Keywords:
NURBSbased Finite Element Method; Computational Fluid Dynamics (CFD); FluidStructure Interaction; Large Eddy Simulation (LES); Rigidbody Dynamics; Arbitrary LagrangianEulerian (ALE) Formulation
Graphical Abstract
Keywords:
NURBSbased Finite Element Method; Computational Fluid Dynamics (CFD); FluidStructure Interaction; Large Eddy Simulation (LES); Rigidbody Dynamics; Arbitrary LagrangianEulerian (ALE) Formulation
1 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., 2013Zienkiewicz O. C., Taylor R. L., Zhu J. Z., (2013). The Finite Element Method for Fluid Dynamics, 7th Edition. ButterworthHeinemann, Oxford, UK.; Reddy and Gartling, 2010Reddy J. N., Gartling D. K., (2010). The Finite Element Method in Heat Transfer and Fluid Dynamics, 3rd Edition. CRC Press, Boca Raton, USA.). 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 smallscale 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 Bsplines) basis functions, which are extensively utilized in Computational Aided Design (CAD) (see Piegl and Tiller, 1997Piegl L., Tiller W., (1997). The NURBS Book, 2nd Edition. SpringerVerlag, Berlin Heidelberg.). 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 ^{p1}continuity order for pdegree basis functions, convex hull and variation diminishing properties (see, for instance, Piegl and Tiller, 1997Piegl L., Tiller W., (1997). The NURBS Book, 2nd Edition. SpringerVerlag, Berlin Heidelberg.). Recent advances have extended NURBS formulation by using Tsplines (Sederberg et al., 2003Sederberg T. N., Zhengs J. M., Bakenov A., Nasri A., (2003). Tsplines and TNURCCSs. ACM Transactions on Graphics 22: 477484.), which permitted local refinement and compatibility of adjacent patches efficiently. By using NURBS functions in a finite element formulation, preprocessing 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, 2014Townsend A., (2014). On the spline: a brief history of the computational curve. International Journal of Interior Architecture  Spatial Design: Applied Geometries, vol. 3. 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 Bspline functions conceived by Schoenberg (1946Schoenberg I. J., (1946). Contributions to the problem of approximation of equidistant data by analytic functions. Quarterly of Applied Mathematics 4: 4599.) were adopted.
Other concepts were introduced with the Bspline 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). Bsplines provide great flexibility for reproduction of geometries since the action of the control points is localized. In addition, Bsplines 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 Bsplines and knot vectors with nonuniform 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, 1997Piegl L., Tiller W., (1997). The NURBS Book, 2nd Edition. SpringerVerlag, Berlin Heidelberg.). NURBSbased finite element formulations have been widely applied to a variety of engineering problems, from solid mechanics (Espath et al., 2014Espath L. F. R., Braun A. L., Awruch A. M., Maghous S., (2014). NURBSbased threedimensional analysis of geometrically nonlinear elastic structures. European Journal of MechanicsA/Solids 47: 373390.; Espath et al., 2015) to material science (Gomez et al., 2008Gomez H., Calo V. M., Bazilevs Y., Hughes T. J. R., (2008). Isogeometric analysis of the CahnHilliard phasefield model. Computer methods in applied mechanics and engineering 197: 43334352.).
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 ^{pm} 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 NURBSbased 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 formulations with NURBS basis functions applied to CFD and FSI problems may be found in Bazilevs et al. (2006Bazilevs Y., Calo V. M., Zhang Y., Hughes T. J. R., (2006). Isogeometric fluidstructure interaction analysis with applications to arterial flow. Computational Mechanics 38: 310322.), Gomez et al. (2010Gomez H., Hughes T. J. R., Nogueira X., Calo V. M., (2010). Isogeometric analysis of the isothermal NavierStokesKorteweg equations. Computer Methods in Applied Mechanics and Engineering 199: 18281840.), Nielsen et al. (2011Nielsen P. N., Gersborg A. R., Gravesen J., Pedersen N. L., (2011). Discretizations in isogeometric analysis of NavierStokes flow. Computer methods in applied mechanics and engineering 200: 32423253.), Akkerman et al. (2011Akkerman I., Bazilevs Y., Kees C. E., Farthing M. W., (2011). Isogeometric analysis of freesurface flow. Journal of Computational Physics 230: 41374152.), and Kadapa et al. (2015Kadapa C., Dettmer W. G., Peric D., (2015). NURBS based leastsquares finite element methods for fluid and solid mechanics. International Journal for Numerical Methods in Engineering 101: 521539.).
Unlike the standard finite element refinement techniques, Bspline refinement can control element size as well as degree and continuity order of the basis functions. The basic refinement techniques associated with Bsplines 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 hrefinement 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 prefinement utilized in finite element modeling must be initially applied over basis functions with C ^{0}continuity, while prefinement for Bspline discretizations can be adopted using any continuity order. This aspect cannot be replicated with a standard finite element formulation. An exclusive technique for Bspline refinement is obtained by using degree elevation followed by knot insertion, which is called krefinement. One can notice that a NURBSbased 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. (2009Cottrell J. A., Hughes T. J. R., Bazilevs Y., (2009). Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009.).
In this work, a NURBSbased finite element formulation is proposed considering the explicit twostep TaylorGalerkin model, where spatial discretization is carried out taking into account Bsplines and NURBS basis functions. The fundamental flow equations are the NavierStokes equations and the mass conservation equation, which is described according to the pseudocompressibility hypothesis for incompressible flows and Newtonian fluids under isothermal conditions. In the present model, idealized 2D turbulent flows are simulated using a LEStype (Large Eddy Simulation) approach, where the Smagorinsky’s model is adopted for subgrid scale modeling. Fluidstructure 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.
2 FUNDAMENTAL EQUATIONS
2.1 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 twodimensional 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 pseudocompressibility hypothesis (see Chorin, 1967Chorin A. J., (1967). A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics, 2: 1226.). An Arbitrary LagrangianEulerian (ALE) formulation is adopted for the kinematical description of fluid motions when fluidstructure interaction is considered, and turbulence flows are approximately reproduced using Large Eddy Simulation (LES) and the Smagorinsky’s subgrid scale model (Smagorinsky, 1963Smagorinsky J., (1963). General circulation experiments with the primitive equations, I. the basic experiment. Monthly Weather Review 91: 99164.). Thus, the system of fundamental flow equations may be written as follows (see, for instance, White, 2005White F. M., (2005). Viscous Fluid Flow, 3rd Edition. Mc GrawHill, New York, USA.):
Momentum balance equation
Mass balance equation
where v_{i} and w_{i} are the components of the flow velocity vector v and mesh velocity vector w, respectively, which correspond to the idirection of a rectangular Cartesian coordinates system, where coordinates components are denoted by x_{i}, 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 (divergencefree velocity field). The fundamental flow equations are valid in the flow spatial domain Ω^{f}, which is bounded by
The fluid stress tensor σ is given according to the Newtonian constitutive model for viscous fluids, that is:
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 LEStype approach, where the fundamental flow equations are submitted to spatial filtering in which the flow field is decomposed into large and smallscale components. Large scales are solved directly utilizing the filtered equations, while scales smaller than the mesh resolution are defined by the subgrid stress tensor τ^{SGS}, which must be modeled considering a turbulence closure model. Turbulence models are utilized in order to represent the smallscale effects over the large scales. The subgrid stress tensor components
where
where C_{S} is the Smagorinsky constant, which must be specified according to the flow characteristics, with values usually ranging from 0.1 to 0.25, and
In order to solve the flow problem, initial conditions on the flow variables v_{i} and p must be specified. In addition, appropriate boundary conditions must also be defined on
where
2.2 Fluidstructure interaction
The fluidstructure coupling is accomplished by enforcing equilibrium and kinematical conditions on the fluidstructure interface
where n_{j} and
where x_{i} 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, 1979Thomas P. D., Lombard C. K., (1979). Geometric conservation law and its application to flow computations on moving grids. AIAA Journal 17: 10301037. for further information). According to Lesoinne and Farhat (1996Lesoinne M., Farhat C., (1996). Geometric conservation law for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations. Computer Methods in Applied Mechanics and Engineering 134: 7090.), 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 fluidstructure coupling, a conservative algorithm where the GCL is satisfied without violating the compatibility conditions can be obtained if the following equations are employed (seeLesoinne and Farhat, 1998Lesoinne M., Farhat C., (1998). Improved staggered algorithms for the serial and parallel solution of threedimensional nonlinear transient aeroelastic problems. In: Computational Mechanics. S Idelsohn, E Oñate and E Dvorakin (Eds.), Barcelona, CIMNE.):
where:
Figure 1 shows the coupling scheme adopted here for fluidstructure 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 (2009Braun A. L., Awruch A. M., (2009). Aerodynamic and aeroelastic analyses on the CAARC standard tall building model using numerical simulation. Computers and Structures 87: 564581.) 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 a_{ij} 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 w_{i} = 0. The influence coefficients are obtained with:
where d_{ij} is the Euclidian distance between the mesh points i and j and exponent n is a userdefined 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 twodimensional 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
The structural motion on the fluidstructure interface is obtained from the structural response evaluated at the center mass by considering the following kinematic relations (see Figure 2):
where
The kinematic relations above may also be expressed using matrix form, that is:
where:
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
3 NUMERICAL MODEL
3.1 NURBSbased finite element formulation
In a NURBSbased 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 (seeFigure 4). For twodimensional 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 nonzero knot spans.
The NURBS basis functions for twodimensional applications are given by:
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 w_{i,j} is related to the weight associated with the control point defined by the subindices i and j. The Bspline basis functions N_{i} are evaluated here using the Coxde Boor recursive formulation (Cox, 1972Cox M. G., (1972). Numerical evaluation of Bsplines. Journal of the Institute of Mathematics and its Applications, 10: 134149.; De Boor, 1972De Boor C., (1972). On calculation with Bsplines. Journal of Approximation Theory 6: 5062.), 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 η.
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 nonzero 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:
Finally, geometry and flow variables are discretized using the following NURBS approximations:
where R_{i} 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 v _{i}, δv _{i}, p, δp and x _{i}. 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.
3.2 The explicit twostep TaylorGalerkin model
In the present model, the fundamental flow equations are discretized using the explicit twostep TaylorGalerkin model, where secondorder Taylor series and the BubnovGalerkin method are adopted for time and space discretization considering a NURBSbased finite element framework.
Time discretization is initially applied to the fundamental flow equations considering a secondorder Taylor series expansion and a twostep 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:
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. R _{v}andR _{p} 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 BubnovGalerkin method, the weight functions are defined here using W _{v} = δv and W _{p} = δ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 GreenGauss theorem. An algebraic system of equations is obtained for the twostep TaylorGalerkin 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 twostep TaylorGalerkin scheme may be found in Kawahara and Hirano (1983Kawahara M., Hirano H., (1983). A finite element method for high Reynolds number viscous fluid flow using two step explicit scheme. International Journal for Numerical Methods in Fluids 3: 137163.).
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 = (x_{1}, x_{2}); the parametric space
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 Δt_{e} is the time increment referring to element e, α is the safety coefficient (0 ≤ α ≤ 1), Δx_{e} is the characteristic length of element e, c is the sound speed, and v_{e} is the flow characteristic speed. Although a multi time step formulation may be adopted for time integration, the smaller Δt_{e} is utilized in this work throughout the finite element mesh.
3.3 Fluidstructure coupling model
In the present work, a conservative partitioned model is adopted for fluidstructure coupling. Considering a fluid element Ω_{E} localized on the fluidstructure 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 fluidstructure interface (I) and mesh points located within the fluid flow domain (F). A matrix form of the rearranged system may be expressed as:
One can notice that only the first line of the system of equations above is relevant for the fluidstructure problem. In addition, the compatibility conditions are also expressed in matrix form considering an element E belonging to the fluidstructure interface as follows:
where:
Notice that N is the number of mesh points referring to element E and the matrices L _{i} and L’_{i} associated with mesh points out of the fluidstructure 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 F _{I} 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:
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:
while the equivalent load vector is given by:
Notice that all the fluid elements in contact with the fluidstructure interface (NELI) are included in the sum operation indicated by Eqs. (67), (68) and (69). Matrix
4 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 postprocessing stage, the commercial software Tecplot 9 was employed.
4.1 Walldriven cavity flow
The walldriven cavity flow problem is analyzed here in order to validate the numerical model proposed in this work for flow simulation with Bspline basis functions. In addition, the influence of mesh refinement over numerical results is also investigated. Incompressible flow regime is assumed and twodimensional flows are considered, which are characterized by the following Reynolds numbers: Re = 10^{3} and Re =10^{4}. Notice that a LEStype approach is employed for Re = 10^{4}, where the Smagorinsky’s constant is set to C_{S} = 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 noslip boundary conditions, while the top wall slides laterally with a constant velocity U = U_{0}.
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 Bspline and NURBS basis functions. The time increment utilized in the time integration of the flow equations is obtained from Eq. (57) with α = 0.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 corresponding predictions obtained by Ghia et al. (1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.), who adopted a secondorder 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. (1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.) 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. (1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.).
Influence of control mesh configuration: velocity profiles along the horizontal and vertical centerlines of the cavity for Re = 10^{3}; basis functions: p = q = 1; refinement parameter: (a) h = 0.001L, (b) h = 0.0025L.
Influence of mesh parameter h: velocity profiles along the horizontal and vertical centerlines of the cavity for Re = 10^{3}; basis functions: p = q = 1; control mesh configuration: (a) 80x80, (b) 100x100.
Influence of basis functions: velocity profiles along the horizontal and vertical centerlines of the cavity for Re = 10^{3}; control mesh configuration: 80x80; mesh parameter: (a) h = 0.0025L, (b) h = 0.005L.
Turbulence is also analyzed in the present application considering a LEStype approach with C_{S} = 0.15 (see Eq. 6) and flow conditions defined by Re = 10^{4}. An intermediate control mesh configuration (100x100) is utilized here with different basis functions (p = q = 1, 2 and 3) and two refinement parameters: h = 0.001L and h = 0.0025L. Results are presented in Figure 8, where comparisons are performed taking into account predictions obtained by Ghia et al. 1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411. for Re = 10^{4}.
Influence of turbulence: velocity profiles along the horizontal and vertical centerlines of the cavity for Re = 10^{4}; control mesh configuration: 100x100; mesh parameter: (a) h = 0.001L, (b) h = 0.0025L.
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. (1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.) 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 v _{1} and v _{2} are velocity vectors obtained from sample points located on the center lines of the cavity according to positions defined by Ghia et al. 1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411. and n is the number of sample points. Notice that each velocity vector corresponds to a specific refinement parameter h_{j} (j = 1,…,4), with h_{1} = 0.001L, h_{2} = 0.0025L, h_{3} = 0.005L, h_{4} = 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 v _{1} and v _{2} 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.
Euclidean norms:
Euclidean norms:
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 v _{1} and v _{2} 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 loworder basis functions, which are traditionally adopted by lagrangianbased 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. (1982Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.), with primary, secondary and tertiary recirculation zones reproduced correctly according to the Reynolds number adopted.
4.2 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, 1966Lienhard J. H., (1966). Synopsis of Lift, Drag and Vortex Frequency Data for Rigid Circular Cylinders. Washington State University, College of Engineering, Research Division Bulletin 300.), a twodimensional 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, 1997Piegl L., Tiller W., (1997). The NURBS Book, 2nd Edition. SpringerVerlag, Berlin Heidelberg.). 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. The time increment Δt is obtained using Eq. (57) with α = 0.3. Computational domain and boundary conditions utilized in the present analysis are shown in Figure 12.
Results referring to drag coefficient (C_{D}) and Strouhal number (St) are presented in Table 2, where predictions obtained in this work are compared with results obtained numerically by Henderson (1997Henderson R., (1997). Nonlinear dynamics and pattern formation in turbulent wake transition. Journal of Fluid Mechanics 352: 65112.). Drag (C_{D}) and lift (C_{L}) 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,
Flow on circular cylinder: drag coefficient (C_{D}) and Strouhal number (St) for 10 ≤ Re ≤ 1000.
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.
Flow around circular cylinder: pressure field and streamlines for different Reynolds numbers.
Distribution of pressure coefficient over the cylinder perimeter for different Reynolds numbers.
The influence of the basis functions on the present results is evaluated using a control mesh configuration with 100x70 control points and flow conditions defined by Re = 40. Two basis function configurations are compared: (a) quadratic functions in both parametric directions (p = q = 2) and (b) quadratic functions in the angular direction (p = 2) with linear functions in the radial direction (q = 1).
Drag coefficient (CD) and geometric characteristics of the recirculation zone behind the circular cylinder for Re = 40.
Results referring to drag coefficient (C_{D}) 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 (2002Wanderley L. B. V., Levi C. A., (2002). Validation of a finite difference method for the simulation of vortexinduced vibrations on a circular cylinder. Ocean Engineering 29: 445460.) 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 nonuniform 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.
4.3 Backwardfacing step flow
The backwardfacing step flow is investigated in the present example using multipatch refinement, where the computational domain is decomposed into parametrically independent subdomains. 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.
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, x_{e} = 1m, x_{t} = 30 m, L_{1} = 1 m, L_{2} = 17 m and L_{3} = 12 m. Noslip 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
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. (1983Armaly B. F., Durst F. J., Pereira C. F., Schönung B., (1983). Experimental and theoretical investigations of backwardfacing step flow. Journal of Fluid Mechanics 127: 473496.) and numerical results obtained by Williams and Baker (1997Williams P. T., Baker A. J., (1997). Numerical simulations of laminar flow over a 3D backwardfacing step. International Journal for Numerical Methods in Fluids 24: 11591183.) with a twodimensional 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 twodimensional 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 threedimensional 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 (1997Williams P. T., Baker A. J., (1997). Numerical simulations of laminar flow over a 3D backwardfacing step. International Journal for Numerical Methods in Fluids 24: 11591183.).
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) and pressure isolines are well defined near the patch interfaces, where oscillations are not identified.
4.4 Lockin analysis for elastically supported circular cylinder
In the present example, the lockin phenomenon observed in elasticallymounted circular cylinders submitted to fluid flow is numerically simulated. It is well known that the lockin phenomenon occurs for a specific range of flow speeds, where synchronization between the mechanical frequency of the structure and the vortexshedding frequency is obtained. During the lockin, the amplitude of oscillations is increased, although rarely exceeding half of the cylinder diameter (see, for instance, Simiu and Scanlan, 1996Simiu E., Scanlan R. H., (1996). Wind Effects on Structures: Fundamentals and Applications to Design, 3rd Edition. John Willey and Sons, New York, USA.). 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 nonuniform 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 lockin phenomenon is investigated using a fixed Reynolds number Re = 150 and different reduced velocities (V_{red} = V/f.D) within the range 3 ≤ V_{red} ≤ 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 u_{y} and C_{y} are the structure displacement and force coefficient transverse to the flow direction, ξ is the damping ratio and M_{red} = M/ρD^{2} is the reduced mass. In the present investigation, no damping is considered and a constant reduced mass M_{red} = 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 lockin response can be clearly identified for V_{red} = 4, while a structural response with relatively low displacement amplitudes are observed when V_{red} = 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, 2006Ahn H. T., Kallinderis Y., (2006). Strongly coupled flow/structure interactions with a geometrically conservative ALE scheme on general hybrid meshes. Journal of Computational Physics 219:671696.; Borazjani et al., 2008Borazjani I., Ge L., Sotiropoulos F., (2008). Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. Journal of Computational Physics 227: 75877620.), where one can see that the lockin interval can be identified within the reduced velocity range 4 ≤ V_{red} ≤ 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 V_{red} = 4 and V_{red} = 8 at t = 80 s. These results are similar to predictions obtained numerically by Borazjani et al. (2008Borazjani I., Ge L., Sotiropoulos F., (2008). Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. Journal of Computational Physics 227: 75877620.) using a FSI model and the immersed boundary method.
Lockin analysis. Instantaneous pressure and vorticity obtained at t = 80s: a) V_{red} = 4; b) V_{red} = 8.
4.5 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 crosssection 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, 1996Simiu E., Scanlan R. H., (1996). Wind Effects on Structures: Fundamentals and Applications to Design, 3rd Edition. John Willey and Sons, New York, USA.).
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, nonuniform 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^{4}. The structure properties are presented in dimensionless form as follows: m* = 20, c* = 0.0581195 and k* = 3.08425. The fluid flow is characterized by a Reynolds number Re = 250, where the dimensionless fluid parameters are defined as: v*∞ = 2.5, ρ* = 1.0 and μ* = 0.01. Notice that asterisks indicate dimensionless values.
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, 2006Dettmer W. G., Peric D., (2006). A computational framework for fluidrigid body interaction: Finite element formulation and applications. Computers Methods in Applied Mechanics and Engineering 195: 16331666.; Robertson et al., 2003Robertson I., Sherwin S. J., Bearman P. W., (2003). A numerical study of rotational and transverse galloping rectangular bodies. Journal of Fluids and Structures 17: 681699.). Results obtained here demonstrate a very good agreement with other predictions considering flow aspects and structural vibration characteristics.
4.6 Flow over elastically supported rectangular cylinder
Fluidrigid 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
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. (2001Sarrate J., Huerta A., Donea J., (2001). Arbitrary LagrangianEulerian formulation for fluidrigid body interaction. Computer Methods in Applied Mechanics and Engineering 190: 31713188.), 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: m_{t} = 195.57, c_{t} = 0.0325, k_{t} = 0.7864, I = 105.94, c_{r} = 0.0 and k_{r}= 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 nonuniform 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 nonuniform 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 nonuniform 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 fluidstructure 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.
Flow over rectangular cylinder: mesh configuration, pressure field and streamlines at t* = 439 and t* = 448.
5 CONCLUSIONS
In the present work a NURBSbased finite element formulation for incompressible fluid dynamics and fluidstructure interaction with rigidbody dynamics was presented. Model versatility with respect to spatial discretization procedures was demonstrated, where different techniques referring to NURBS discretization, such as multipatches, 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 nonuniform 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 NURBSbased finite element models may also be developed. For future works, the present formulation must be extended to include the threedimensional approach. Potential applications for the present model may be found in the field of Computational Wind Engineering, such as aerodynamic and aeroelastic analyses of longspan bridges, low and highrise buildings.
Acknowledgments
The authors would like to thank CNPq and CAPES for the financial support, and CESUPUFRGS and CENAPADUNICAMP for technical support and computational resources provided.
References
 Ahn H. T., Kallinderis Y., (2006). Strongly coupled flow/structure interactions with a geometrically conservative ALE scheme on general hybrid meshes. Journal of Computational Physics 219:671696.
 Akkerman I., Bazilevs Y., Kees C. E., Farthing M. W., (2011). Isogeometric analysis of freesurface flow. Journal of Computational Physics 230: 41374152.
 Armaly B. F., Durst F. J., Pereira C. F., Schönung B., (1983). Experimental and theoretical investigations of backwardfacing step flow. Journal of Fluid Mechanics 127: 473496.
 Bathe J., (1996). Finite Element Procedures, 1st Edition. Prentice Hall, Englewood CliffsNJ, USA.
 Bazilevs Y., Calo V. M., Zhang Y., Hughes T. J. R., (2006). Isogeometric fluidstructure interaction analysis with applications to arterial flow. Computational Mechanics 38: 310322.
 Bezier P., (1972). Numerical Control: Mathematics and Applications. Wiley, London.
 Borazjani I., Ge L., Sotiropoulos F., (2008). Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. Journal of Computational Physics 227: 75877620.
 Braun A. L., Awruch A. M., (2009). Aerodynamic and aeroelastic analyses on the CAARC standard tall building model using numerical simulation. Computers and Structures 87: 564581.
 Chorin A. J., (1967). A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics, 2: 1226.
 Cottrell J. A., Hughes T. J. R., Bazilevs Y., (2009). Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009.
 Cox M. G., (1972). Numerical evaluation of Bsplines. Journal of the Institute of Mathematics and its Applications, 10: 134149.
 De Boor C., (1972). On calculation with Bsplines. Journal of Approximation Theory 6: 5062.
 Dettmer W. G., Peric D., (2006). A computational framework for fluidrigid body interaction: Finite element formulation and applications. Computers Methods in Applied Mechanics and Engineering 195: 16331666.
 Espath L. F. R., Braun A. L., Awruch A. M., Dalcin L. D., (2015). A NURBSbased finite element model applied to geometrically nonlinear elastodynamics using a corotational approach. International Journal for Numerical Methods in Engineering 102: 18391868.
 Espath L. F. R., Braun A. L., Awruch A. M., Maghous S., (2014). NURBSbased threedimensional analysis of geometrically nonlinear elastic structures. European Journal of MechanicsA/Solids 47: 373390.
 Ghia U., Ghia K. N., Shin C. T., (1982). HighRe solutions for incompressible flow using the NavierStokes equations and multigrid method. Journal of Computational Physics 48: 378411.
 Gomez H., Calo V. M., Bazilevs Y., Hughes T. J. R., (2008). Isogeometric analysis of the CahnHilliard phasefield model. Computer methods in applied mechanics and engineering 197: 43334352.
 Gomez H., Hughes T. J. R., Nogueira X., Calo V. M., (2010). Isogeometric analysis of the isothermal NavierStokesKorteweg equations. Computer Methods in Applied Mechanics and Engineering 199: 18281840.
 Henderson R., (1997). Nonlinear dynamics and pattern formation in turbulent wake transition. Journal of Fluid Mechanics 352: 65112.
 Kadapa C., Dettmer W. G., Peric D., (2015). NURBS based leastsquares finite element methods for fluid and solid mechanics. International Journal for Numerical Methods in Engineering 101: 521539.
 Kawahara M., Hirano H., (1983). A finite element method for high Reynolds number viscous fluid flow using two step explicit scheme. International Journal for Numerical Methods in Fluids 3: 137163.
 Lesoinne M., Farhat C., (1996). Geometric conservation law for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations. Computer Methods in Applied Mechanics and Engineering 134: 7090.
 Lesoinne M., Farhat C., (1998). Improved staggered algorithms for the serial and parallel solution of threedimensional nonlinear transient aeroelastic problems. In: Computational Mechanics. S Idelsohn, E Oñate and E Dvorakin (Eds.), Barcelona, CIMNE.
 Lienhard J. H., (1966). Synopsis of Lift, Drag and Vortex Frequency Data for Rigid Circular Cylinders. Washington State University, College of Engineering, Research Division Bulletin 300.
 Nielsen P. N., Gersborg A. R., Gravesen J., Pedersen N. L., (2011). Discretizations in isogeometric analysis of NavierStokes flow. Computer methods in applied mechanics and engineering 200: 32423253.
 Piegl L., Tiller W., (1997). The NURBS Book, 2nd Edition. SpringerVerlag, Berlin Heidelberg.
 Reddy J. N., Gartling D. K., (2010). The Finite Element Method in Heat Transfer and Fluid Dynamics, 3rd Edition. CRC Press, Boca Raton, USA.
 Robertson I., Sherwin S. J., Bearman P. W., (2003). A numerical study of rotational and transverse galloping rectangular bodies. Journal of Fluids and Structures 17: 681699.
 Sarrate J., Huerta A., Donea J., (2001). Arbitrary LagrangianEulerian formulation for fluidrigid body interaction. Computer Methods in Applied Mechanics and Engineering 190: 31713188.
 Schoenberg I. J., (1946). Contributions to the problem of approximation of equidistant data by analytic functions. Quarterly of Applied Mathematics 4: 4599.
 Sederberg T. N., Zhengs J. M., Bakenov A., Nasri A., (2003). Tsplines and TNURCCSs. ACM Transactions on Graphics 22: 477484.
 Simiu E., Scanlan R. H., (1996). Wind Effects on Structures: Fundamentals and Applications to Design, 3rd Edition. John Willey and Sons, New York, USA.
 Smagorinsky J., (1963). General circulation experiments with the primitive equations, I. the basic experiment. Monthly Weather Review 91: 99164.
 Townsend A., (2014). On the spline: a brief history of the computational curve. International Journal of Interior Architecture  Spatial Design: Applied Geometries, vol. 3.
 Thomas P. D., Lombard C. K., (1979). Geometric conservation law and its application to flow computations on moving grids. AIAA Journal 17: 10301037.
 Wanderley L. B. V., Levi C. A., (2002). Validation of a finite difference method for the simulation of vortexinduced vibrations on a circular cylinder. Ocean Engineering 29: 445460.
 White F. M., (2005). Viscous Fluid Flow, 3rd Edition. Mc GrawHill, New York, USA.
 Williams P. T., Baker A. J., (1997). Numerical simulations of laminar flow over a 3D backwardfacing step. International Journal for Numerical Methods in Fluids 24: 11591183.
 Zienkiewicz O. C., Taylor R. L., Zhu J. Z., (2013). The Finite Element Method for Fluid Dynamics, 7th Edition. ButterworthHeinemann, Oxford, UK.
Publication Dates

Publication in this collection
03 Feb 2020 
Date of issue
2020
History

Received
20 Aug 2019 
Reviewed
25 Nov 2019 
Accepted
27 Nov 2019 
Published
05 Dec 2019