Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0033 Acta Polytechnica 61(SI):33–48, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague NUMERICAL SIMULATION OF LID-DRIVEN CAVITY FLOW BY ISOGEOMETRIC ANALYSIS Bohumír Bastl∗, Marek Brandner, Jiří Egermaier, Hana Horníková, Kristýna Michálková, Eva Turnerová University of West Bohemia, Faculty of Applied Sciences, Department of Mathematics and European Centres of Excellence New Technologies for the Information Society, Univerzitní 8, 301 00, Plzeň, Czech Republic ∗ corresponding author: bastl@kma.zcu.cz Abstract. In this paper, we present numerical results obtained by an in-house incompressible fluid flow solver based on isogeometric analysis (IgA) for the standard benchmark problem for incompressible fluid flow simulation – lid-driven cavity flow. The steady Navier-Stokes equations are solved in their velocity-pressure formulation and we consider only inf-sup stable pairs of B-spline discretization spaces. The main aim of the paper is to compare the results from our IgA-based flow solver with the results obtained by a standard package based on finite element method with respect to degrees of freedom and stability of the solution. Further, the effectiveness of the recently introduced rIgA method for the steady Navier-Stokes equations is studied. The authors dedicate the paper to Professor K. Kozel on the occasion of his 80th birthday. Keywords: Isogeometric analysis, fluid flow simulation, B-splines, Navier-Stokes equations. 1. Introduction The fluid flow simulation is one of the fundamen- tal problems solved in engineering practice. During the history, many different numerical methods were proposed to solve this extremely challenging prob- lem, ranging from Finite Difference Methods (FDM), through Finite Volume Methods (FVM) to Finite El- ement Methods (FEM) with many variants of these methods (see e.g. [1–3]). Relatively recently, Hughes et al. (see [4, 5]) introduced a powerful numerical method for solving partial differential equations called isogeometric analysis (IgA). This method has many features common with FEM and the main motiva- tion is to bridge the gap between geometric modelling (Computer-Aided Design – CAD) and numerical sim- ulation (Finite Element Analysis – FEA). Optimally, one would like to use directly the description of an object in CAD system also for numerical simulations. However, this is not possible with the standard ap- proaches (FDM, FVM, FEM) since they require the construction of some computational mesh (typically, triangular/quadrilateral meshes for 2D problems and tetrahedral/hexahedral for 3D problems) which is, moreover, only an approximation of the original com- putational domain. The isogeometric analysis is based on an isopara- metric framework, i.e., the same basis functions are used for representation of the geometry and for the solution space. Currently, the de facto industrial stan- dard for representation of complex objects in CAD systems are Non-Uniform Rational B-spline (NURBS) objects (see [6]). This is the reason why isogeomet- ric analysis was formulated specifically for this class of objects. They allow for the exact representation of a wide variety of complex shapes in a compact form, i.e., by a small number of control points. Even though the NURBS description of an object in CAD system can be used directly for numerical simula- tions only in case when IgA-BEM (the isogeometric formulation of Boundary Element Method, see e.g. [7]) is used, the exact NURBS representation is still beneficial. Many methods providing geometrically exact bivariate/trivariate NURBS representation of planar/spatial domains, respectively, from its bound- ary were proposed in the last years (see e.g. [8–13]). Such NURBS representations then serve as computa- tional meshes used in solvers based on isogeometric analysis. Other advantages of the isogeometric analy- sis include e.g. higher accuracy of the solution with respect to a degree of freedom ([14]) or favourableness of higher continuity of NURBS basis functions for solving higher order PDEs ([15]). From a practical point of view, a typical implementation of IgA-based solver automatically allows for basis functions of the geometry and the solution space of any degree and continuity, without any additional effort needed for implementing the solver. This is typically not the case for codes based on FEM, where in p-version of FEM each degree of the basis functions of the solution space needs to be added/implemented individually. Because of many advantages of the IgA, the method has been successfully applied in various practical problems, like linear elasticity, structural vibrations, phase transi- tion phenomena, fluid flow simulation, plate and shell analysis, heat transfer analysis, shape optimization, etc., see e.g. [15–33]. On the other hand, several disadvantages of the isogeometric analysis were identified during the ex- tensive research in the last years. For example, the 33 https://doi.org/10.14311/AP.2021.61.0033 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica properties of NURBS objects induce that true local refinement is not possible for this class of objects. This drawback can be overcome by various generaliza- tions of NURBS objects, like T-splines, THB-splines or LR-splines (see e.g. [34]), which were recently stud- ied in the context of isogeometric analysis. Among other drawbacks belong e.g. time-consuming matrix assembly in the construction of the matrix of linear systems arising in Galerkin isogeometric methods (see e.g. [35–37]), degraded performance of direct solvers because of wide support of B-spline/NURBS basis functions (see [38]) or the necessity to modify exist- ing well-known stabilization techniques to suppress spurious oscillations. In this paper, we focus on one of the standard bench- mark problems for incompressible fluid flow simulation – lid-driven cavity flow. This problem is frequently studied in the literature. A well-done overview of the results for different Reynolds numbers obtained by different methods can be found e.g. in the recent paper [39], but even this extensive list is not complete. As it can be seen from the summary of results for this problem in [39], FDM is quite frequently used for ob- taining very accurate results for this problem. To the best of our knowledge, there exist also several papers which solve the lid-driven cavity problem by isogeo- metric analysis. The paper [15] applies isogeometric analysis to solve the fourth order stream function for- mulation of the Navier-Stokes equations and presents the results for Reynolds numbers up to 5000. Further, the paper [40] uses divergence-conforming B-splines which may be interpreted as smooth generalizations of Raviart-Thomas elements and presents the results for Reynolds numbers up to 1000. Finally, the paper [33] is based on velocity-pressure formulation of steady and unsteady Navier-Stokes equations, uses LBB compati- ble B-spline spaces which may be regarded as smooth generalizations of Taylor-Hood pairs of finite element spaces and presents the results for Reynolds numbers up to 1000. The main goal of this paper is to extend the study done in [33] for the lid-driven cavity problem to higher Reynolds number, but in a slightly different manner. We want to compare the solutions obtained by our in- house fluid flow solver based on isogeometric analysis (see [30, 31]) for different B-spline discretization spaces with solutions obtained by finite element method with respect to degrees of freedom and stability of the solu- tion. Further, we want to study the potential “stabi- lization effect” of high continuity B-spline discretiza- tion spaces (see [41]), i.e., if high continuity B-spline discretization spaces provide solutions with less oscilla- tions than standard finite element spaces and B-spline discretization spaces of low continuity (FEM-like) for the similar number of degrees of freedom. Moreover, we want to investigate if this potential stabilization effect of the high continuity B-spline discretizations dominates the effect of the SUPG stabilization for the low continuity B-spline discretizations. Finally, inspired by the recently introduced rIgA method (see [42, 43]), we want to explore the effectiveness of rIgA method for solving the Navier-Stokes equations, which has not yet been studied. The paper is organized as follows. Section 2 in- troduces the steady Navier-Stokes equations, their weak formulation and Galerkin discretization. In Section 3, we overview the fundamentals of B-spline objects and present the main ideas of isogeometric analysis. Section 4 is devoted to a very brief overview of Streamline Upwind/Petrov-Galerkin method used for stabilization of a numerical solution of convection- dominated problems typically for lower degree dis- cretization spaces. In Section 5, we formulate the benchmark problem studied in the paper – the lid- driven cavity flow – and present the reference solutions used for comparison with our results. Section 6 then demonstrates the main results of the paper, i.e., the results obtained by our in-house fluid flow solver based on isogeometric analysis and their comparison with reference solutions and other solutions obtained by a standard FEM package. In Section 7, we conclude the paper. 2. Navier-Stokes equations The mathematical model of incompressible viscous Newtonian flow is based on the incompressible Navier– Stokes equations (NSE). Let Ω ⊂ Rd be a bounded domain, d being the number of spatial dimensions, with the boundary ∂Ω consisting of two complemen- tary parts, Dirichlet ∂ΩD and Neumann ∂ΩN . The steady-state incompressible Navier–Stokes problem is given as a system of d + 1 differential equations together with boundary conditions −ν∆u + (u · ∇)u +∇p = 0 in Ω, ∇ · u = 0 in Ω, u = gD on ∂ΩD, ν ∂u ∂n − np = 0 on ∂ΩN , (1) where u is the flow velocity, p is the kinematic pressure, ν is the kinematic viscosity and gD is a given function. If the velocity is specified everywhere on the boundary, the pressure solution is only unique up to a hydrostatic constant. The laminar/turbulent behaviour of the fluid can be predicted by the dimensionless Reynolds number generally defined as the ratio of inertial forces to viscous forces Re = ûL ν , (2) where û is the mean velocity of the fluid and L is the characteristic linear dimension. 2.1. Galerkin discretization In order to derive the weak formulation of the prob- lem (1), we define velocity solution space V and test function space V0 as follows V = {u ∈ H1(Ω)d |u = gD on ∂ΩD}, V0 = {v ∈ H1(Ω)d |v = 0 on ∂ΩD}, (3) 34 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis where the boundary equalities are understood in the sense of traces and we assume that gD ∈ H1/2(∂Ω)d. By multiplying the first line of (1) by a test func- tion v ∈ V0 and the second line by a test function q ∈ L2(Ω), integrating over Ω and using Green’s the- orem, we obtain the weak formulation: find u = (u1, . . . , ud) ∈ V and p ∈ L2(Ω) such that ν ∫ Ω ∇u : ∇v + ∫ Ω (u · ∇u) · v− ∫ Ω p∇ · v = 0,∫ Ω q∇ · u = 0 (4) for all v = (v1, . . . , vd) ∈ V0 and q ∈ L2(Ω), where ∇u : ∇v = d∑ i=1 ∇ui · ∇vi = d∑ i=1 d∑ j=1 ∂ui ∂xj ∂vi ∂xj , (u · ∇u) · v = d∑ i=1 d∑ j=1 uj ∂ui ∂xj vi. The idea of Galerkin discretization method is to define finite dimensional subspaces V h ⊂ V , V h0 ⊂ V0, Wh ⊂ L2(Ω) together with their basis functions and solve the problem projected into these subspaces. Then, the discrete weak solutions uh ∈ V h and ph ∈ Wh are searched such that ν ∫ Ω ∇uh : ∇vh + ∫ Ω (uh · ∇uh) · vh− − ∫ Ω ph∇ · vh = 0,∫ Ω qh∇ · uh = 0 (5) for all vh ∈ V h0 and qh ∈Wh. Because of the non-linearity in the convective term, the problem needs to be linearized and solved iter- atively. Using Picard’s iteration method, we search un+1 h ∈ V h and pn+1 h ∈ Wh as a solution of the fol- lowing Oseen problem ν ∫ Ω ∇un+1 h : ∇vh + ∫ Ω (unh · ∇un+1 h ) · vh − (6) − ∫ Ω pn+1 h ∇ · vh = 0, ∫ Ω qh∇ · un+1 h = 0 for all vh ∈ V h0 and qh ∈ Wh. The initial guess u0 is often chosen to be zero, which means solving the Stokes equations in the first Picard iteration. Particular Galerkin-based discretization methods are then defined by the choice of the subspaces and their basis functions. Considering the basis of V h0 formed by vector functions (Qui (x), 0, . . . , 0︸ ︷︷ ︸ d ), . . . , (0, . . . , 0, Qui (x)︸ ︷︷ ︸ d ), i = 1, . . . , nu and the basis of Wh formed by functions Qpi (x), i = 1, . . . , np, we obtain after discretization a non- symmetric saddle-point linear system of the form[ F −BT B 0 ] [ u p ] = [ f g ] , (7) where F ∈ Rd·nu×d·nu is block diagonal with the di- agonal blocks consisting of the discretization of the viscous term and the linearized convective term, i.e., F = diag ( A + N(unh), . . . ,A + N(unh)︸ ︷︷ ︸ d ) and A = [Aij ] = ν ∫ Ω ∇Qui (x) · ∇Quj (x)dΩ  , N(u) = [Nij(u)] = ∫ Ω Qui (x) ( u · ∇Quj (x) ) dΩ  , BT ∈ Rd·nu×np and B ∈ Rnp×d·nu consist of d blocks, i.e., B = [ B1, . . . ,Bd ] and Bm = [Bmij ] = ∫ Ω Qpi (x) ( ∇Quj (x) · em ) dΩ  . As the right-hand side of (4) is equal to zero, f and g only reflect the given boundary conditions. 2.2. LBB condition To get a stable Galerkin discretization of a saddle- point problem, the so-called LBB (or inf-sup) condi- tion (see [44]) inf q∈Wh\{0} sup v∈V h\{0} ∫ Ω q∇ · v ||v||H1(Ω)d ||q||L2(Ω) ≥ β > 0 (8) has to be satisfied, where β is a constant indepen- dent of the mesh. Here, we consider only the case of so-called inf-sup stable combinations of velocity and pressure spaces ensuring that the inf-sup condition is satisfied. In standard finite elements, there are several known types of stable elements. One of the most popular is a Taylor–Hood element, where basis functions of degree k are used to approximate the pressure and basis functions of degree k + 1 are used for the velocities. 3. Isogeometric analysis In this section, we present fundamentals of isogeomet- ric analysis with a special focus on the problem solved in this paper – fluid flow simulation in a lid-driven cavity. 35 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica 0.5 1.0 1.5 2.0 2.5 3.0 0.2 0.4 0.6 0.8 1.0 (a). 0.5 1.0 1.5 2.0 2.5 3.0 0.2 0.4 0.6 0.8 1.0 (b). 0.5 1.0 1.5 2.0 2.5 3.0 0.2 0.4 0.6 0.8 1.0 (c). Figure 1. Comparison of 1D quadratic discretization bases on a “mesh” with three elements: (A) finite element basis, (B) isogeometric basis with C1 interele- ment continuity, (C) isogeometric basis with C0 in- terelement continuity. B-spline and NURBS objects are described using a B-spline basis. Its construction is based on a so- called knot vector Ξ = (ξ1, ξ2, . . . , ξn+k+1), where ξi ≤ ξi+1 are the knots, k is the basis degree and n is the number of basis functions. Note that the knots can be multiple. The i-th B-spline basis function of degree k, Ni,k(ξ), is defined recursively by Cox-de Boor formula as follows: for k = 0 : Ni,0(ξ) = { 1, ξi ≤ ξ < ξi+1, 0, otherwise, for k > 0 : Ni,k(ξ) = si,k(ξ)Ni,k−1(ξ)+ + ( 1− si+1,k(ξ) ) Ni+1,k−1(ξ), where si,k(ξ) :=  ξ − ξi ξi+k − ξi , ξi < ξi+k, arbitrary, otherwise. (9) Note that Ni,k−1(ξ) = 0 if ξi = ξi+k. One of the fundamental properties of B-spline basis functions is that they are Ck−mi continuous at the knot ξi, where mi is the multiplicity of the knot ξi. Other properties include non-negativity, partition of unity and local support. Generally, IgA-based computations are based on B- spline/NURBS surfaces (in 2D) or B-spline/NURBS volumes (in 3D). Moreover, as it is well-known that the tensor-product structure of B-spline/NURBS ob- jects prevents the possibility to describe an arbitrary computational domain by one B-spline/NURBS ob- ject. This is the reason why so-called multi-patch B-spline/NURBS meshes are typically used for the description of computational domains in isogeometric analysis. Our in-house solver used for obtaining results in Section 6 works with general, free-form multi-patch 2D and 3D B-spline domains. However, in this paper we focus on a geometrically very simple domain – cavity – which allows (under some other conditions mentioned bellow) to use only a so-called B-spline function. Let us consider a two-dimensional domain Ω = 〈0, 1〉 × 〈0, 1〉 ⊂ R2 representing the cavity (see Sec- tions 5 and 6) and let us consider a B-spline function f(ξ, ψ) of degree k and order of continuity r in both parameters ξ and ψ in the form f(ξ, ψ) = n∑ i=1 m∑ j=1 ci,jNi,k(ξ)Nj,k(ψ) = n∑ i=1 m∑ j=1 ci,jQi,j(ξ, ψ), (10) where ci,j are coefficients of the B-spline function, Ni,k(ξ) and Nj,k(ψ) are B-spline basis functions corre- sponding to knot vectors Ξ = (ξ1, ξ2, . . . , ξn+k+1),Ψ = (ψ1, ψ2, . . . , ψm+k+1), respectively, and the product of two univariate B-spline basis functions Ni,k(ξ)Nj,k(ψ) is denoted as one bivariate function Qi,j(ξ, ψ). As in our experiments we use open knot vectors in the form (0, . . . , 0︸ ︷︷ ︸ k+1 , tk+2, . . . , tn, 1, . . . , 1︸ ︷︷ ︸ k+1 ), where k is the degree of basis functions, the parametric domain of the B-spline function f and the computa- tional domain Ω coincide (both are represented by the unit square). Moreover, various degrees and con- tinuities of the B-spline computational meshes used in experiments in Section 6 are obtained by degree elevation and knot insertion (see [6]) starting from bilinear B-spline description of Ω. This means that for all used B-spline computational meshes a point (ξ, ψ) in the parametric domain coincides with the point (x, y) in the physical domain. Altogether, the above mentioned facts allow us to limit to B-spline functions, instead of more general B-spline surfaces. Let us denote the tensor product B-spline space spanned by the basis functions Qi,j(ξ, ψ) as Skr = span{Qi,j}n,mi=1,j=1 (11) 36 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis and the Cartesian product of these tensor product B-spline spaces as Sk r = Skr × Skr . For simplicity of notation, let us re-index the coeffi- cients and the basis functions with one index ` such that f(ξ, ψ) = N∑ `=1 c`Q`(ξ, ψ), N = n ·m. (12) A B-spline mesh used in isogeometric analysis is al- ready defined by the knot vectors Ξ,Ψ. If the subse- quent knots are different from each other, i.e. ξi 6= ξi+1 and ψj 6= ψj+1, then [ξi, ξi+1]× [ψj , ψj+1] defines an element of the computational mesh. Let us now define the discrete spaces V h, V h0 ,Wh as subspaces of V, V0,W generated by the basis functions defining the geometry of Ω. The velocity and pressure spaces must satisfy the condition (8). The inf-sup stability of some combinations of the discrete solution spaces in isogeometric analysis is addressed e.g. in [45, 46]. As examples of stable pairs of B-spline dis- cretization spaces, let us mention Taylor–Hood (TH) pairs of B-spline spaces (an analogue of the Taylor– Hood pairs of finite element spaces) which can be defined by taking the pressure basis functions equal to the geometry basis functions (i.e., Qp` (ξ, ψ) = Q`, ∀`) and the velocity basis functions {Quj (ξ, ψ)} as the B-spline basis functions obtained by degree elevation of the pressure basis functions {Qp` (ξ, ψ)}. Hence, the velocity basis functions are of bi-degree (k + 1, k + 1), both velocity and pressure bases are defined on the same mesh and have the same order of continuity. This yields the following definition V h = {u : u ∈ Sk+1 r ∧ u = gD on ∂ΩD}, V h0 = {v : v ∈ Sk+1 r ∧ v = 0 on ∂ΩD}, (13) Wh = {q : q ∈ Skr }. The discrete velocity uh ∈ V h and pressure ph ∈ Wh can be written as linear combinations of the basis functions, i.e., uh = Nu f∑ i=1 uiQui (x, y) + Nu∑ i=Nu f +1 u∗iQui (x, y),(14) ph = Np∑ j=1 pjQ p j (x, y), (15) where Qui (x, y) = Qui (ξ, ψ),∀i, Qpj (x, y) = Qpj (ξ, ψ),∀j (specifically in this case, see (10) and the discussion below) and ui ∈ R2 and pj ∈ R are unknown coeffi- cients, Nu, Np are number of velocity and pressure basis functions, respectively, and Nu f is the number of "free" velocity basis functions whose coefficients are not fixed due to Dirichlet boundary conditions. The coefficients u∗i come from a B-spline representation of the Dirichlet data at the boundary, which is assumed as known. After substituting (14) and (15) into the discrete linearized weak formulation, we obtain a linear system of the form (7) for the unknown coefficients ui and pj . Similarly to the finite element method, the sys- tem matrix is sparse thanks to local supports of the B-spline basis functions. For degree k equal to 1 the B-spline basis is the same as the corresponding finite element basis. For higher degrees, however, they differ from each other. For illustration, see the comparison of quadratic FEM and IgA bases on a 1D “mesh” consisting of three equal elements in Figure 1. Fig- ure 1(a) shows the standard quadratic FEM basis and Figure 1(b) shows a C1 continuous quadratic B-spline basis for the open knot vector Ξ = (0, 0, 0, 1, 2, 3, 3, 3). Obviously, we have less basis functions for the same mesh in isogeometric analysis, they are all point- wise non-negative and generally non-interpolatory, i.e. their coefficients do not correspond to nodal values of the IgA solution. Moreover, their values form partition of unity for arbitrary ξ ∈ (ξ1, ξn+k+1) and the support of most of the B-spline basis functions (except the ones close to the boundary) is k+ 1 knot spans. Both FEM and IgA bases result in the same bandwidth, since the maximum number of basis functions that a given basis function shares support with (including itself) is 2k + 1 in both cases. However, the FEM matrix is sparser since there are basis functions with supports containing only one element, thus sharing support with only k+ 1 basis functions. From the ma- trix structure point of view, a B-spline analogue of the FEM basis is the B-spline basis with C0 interelement continuity displayed in Figure 1(c), which is obtained for the knot vector Ξ = (0, 0, 0, 1, 1, 2, 2, 3, 3, 3). 4. Streamline Upwind/Petrov-Galerkin (SUPG) stabilization Since the isogeometric analysis is continuous Galerkin- based method, the numerical solution of convection dominated problems (for high Reynolds numbers) is usually polluted by spurious (unphysical) oscillations, which cause a loss of accuracy and stability. Stabi- lization techniques are investigated which improve the stability, without degrading accuracy. There are many linear and nonlinear stabilization techniques for convection-dominated problems. Their overview can be found, e.g., in [47]. In our case, we focused only on the SUPG method. The SUPG method proposed in [48] has become the most popular stabilization method. The weighted element residual is added to the Galerkin problem (5), i.e., nel∑ ie=1 ∫ Qie τSR(uh)uh · ∇vh, (16) 37 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica where R(uh) = uh · ∇uh − ν∆uh +∇ph (17) is the element residual, nel is the number of elements Qie , and τS is a non-negative SUPG stabilization parameter (cf. e.g. [49, 50]) τS = (( 2 deg(uh) ‖uh‖ h )2 + 9 ( 4ν h2 )2 )− 1 2 , (18) where deg(uh) is the degree of the velocity basis func- tions and h is the characteristic size of the element. Alternatives to the stabilization parameter τS are dis- cussed e.g. in [51]. Note that the second derivatives are not zero in the interior of the elements if higher order B-spline basis functions are used. 5. Lid-Driven Cavity Flow In this section, we introduce the standard benchmark problem for incompressible fluid flow simulation – lid- driven cavity flow – for which we present the results of our numerical experiments in Section 6. One of the reasons for this choice is a wide range of comparative data available in the literature. A thorough overview of results for lid-driven cavity flow is done in [39]. The authors summarize the papers dealing with this problem starting from 1966 up to the most recent results in 2014. For comparisons, the results presented in [52] are frequently used, where the authors used FDM and solved the problem for Reynolds numbers up to 10000. Other papers, which are usually used for comparisons are e.g. [53, 54]. For comparison with our results we chose two ref- erences. The first one is the above mentioned recent paper [39], where the authors solve lid-driven cavity flow by FVM for Reynolds numbers up to 5000. The results can be considered as very accurate because a very fine mesh with 1301×1301 grid points for the case of Reynolds number equal to 5000 was used (compare with [52], where FDM with 257× 257 grid points was used). The results presented in [55] is another source to which we compare our data. To the best of our knowledge, these results are unpublished in a journal. Here, the author uses FDM to solve the problem on several very fine graded grids ranging from 1201×1201 to 1661× 1661 grid points. Then, Richardson extrap- olation is used for improving the solutions’ accuracy and the extrapolated result is projected onto the grid of 1201× 1201 points. Comparing these results with [39], they are also considered as very accurate. Fig- ure 2 shows the comparison of the results presented in [39] and [55] where ux- and uy-velocity profiles along vertical/horizontal lines through the geometric center of the cavity are compared. The reason why we use also reference [55] is that other available sources usually present velocity pro- files only in the above mentioned two cross-sections through the computational domain – ux-velocity pro- files along a vertical line through the geometric center of the cavity and uy-velocity profiles along a horizon- tal line through the geometric center of the cavity. However, [55] presents velocity profiles also for some other horizontal and vertical cross-sections through the computational domain, which we will also use for the presentation of our results and for comparison in Section 6. The problem description is shown in Figure 3. Lid- driven cavity flow represents a flow inside a square container which is covered by a lid that slides sideways. From the existing literature, it is commonly accepted that there exists a steady solution of the problem for Reynolds numbers up to 8000 (in 2D, see e.g. [56, 57]). In our study presented in this paper, we use the kinematic viscosity ν = 0.0002 which implies Reynolds number equal to 5000. Generally, homogeneous Dirichlet boundary condi- tion for velocity is set on the left, right and bottom boundaries of the computational domain and constant velocity u = (1, 0) is prescribed on the lid part (top) of the boundary. However, since gD is a discontinuous function in this case, it does not lie in H1/2(∂Ω)d and thus the velocity solution cannot be searched in the space V defined in (3), see e.g. [58] for more details on this issue. For discretizations of such singular problem, slower convergence can be expected when refining the mesh. In order to ensure well-posedness of the problem, a so-called regularization is often used, which means that the discontinuous function gD is replaced by some suitable function from H1/2(∂Ω)d. For instance, we can define a function gεD such that its first component changes linearly from 0 to 1 on a part of the boundary of width ε near the upper corners as in [59]. However, in this paper, we discretize the original singular problem, i.e., no regularization is used here. Thus, we can compare our IgA solutions with the results available in the references mentioned above, where the singular lid-driven cavity is also considered. In fact, after IgA discretization, the finite dimensional problem is solved with a continuous boundary condi- tion gDh defined as a B-spline approximation of gD on the given mesh. In principle, there are two ways to approximate the discontinuous boundary condition, called “leaky” and “non-leaky” cavity. The leaky cavity approach considers the upper corners as part of the lid and the first component of gDh goes from 0 to 1 on some part of the vertical walls (dependent on the mesh) near the upper corners. We consider the non-leaky cavity approach, where, on the contrary, the upper corners are considered as part of the vertical walls. As Dirichlet boundary conditions for velocity are imposed on the whole boundary of the computational domain, the pressure is determined up to an arbitrary constant. Thus, to have the discrete pressure field uniquely defined, we fix the value of pressure to zero at the bottom left corner. 38 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis (a). (b). Figure 2. Comparison of reference solutions (Y. Papadopoulos stands for [55], T. A. AbdelMigid et al. for [39]): ux-velocity profile along a vertical line at x = 0.5 (A), uy-velocity profile along a horizontal line at y = 0.5 (B). Figure 3. Lid-driven cavity flow – the problem de- scription. 6. Results In this section, we present the numerical results ob- tained by our in-house isogeometric fluid flow solver for the model problem, lid-driven cavity flow, and compare them with reference solutions and with the numerical results obtained by standard package FEn- iCS (see [60]), which is based on the finite element method. Let us note that the simplicity of this bench- mark example is in favor of FEM because it is possible to represent the computational domain exactly also by finite element meshes. On more complicated domains, where finite element meshes only approximate the domain boundary, isogeometric analysis have an ad- vantage of exact representation of the computational domain by a B-spline/NURBS mesh. S2 0×S1 0 S4 0×S3 0 S4 2×S3 2 32× 32 9026 41666 9674 64× 64 36482 167298 37770 128× 128 146690 670466 149258 Table 1. Degrees of freedom (DOFs) for simulations on different meshes in the in-house IgA-based flow solver. Our in-house isogeometric fluid flow solver is imple- mented in C++ in the framework of G+Smo (Geom- etry and Simulation modules) library (see [61]). For all the simulations performed with our solver we used the kinematic viscosity ν = 0.0002 yielding Re = 5000 and we compare them with the results obtained by FEniCS and with the results of the reference solutions presented in [39] or [55] for the same Reynolds num- ber. Specifically, the results from [39] are used for the comparison of ux- and uy-velocity profiles along cross-sections of the computational domain at x = 0.5 and y = 0.5, respectively, and the results from [55] for comparison of velocity profiles along all other cross- sections of the computational domain. We tested several different inf-sup stable pairs of B-spline discretization spaces, i.e., • S2 0×S1 0 with linear basis functions for pressure and quadratic for velocity, which is close to the standard FEM, 39 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica Figure 4. Stream function (top) and vorticity (bot- tom) profiles obtained for S4 2 × S3 2 on the mesh com- posed of 128× 128 elements. The contour ranges for stream function and vorticity are ψiso ∈ [−0.12, 0.003] and ωiso ∈ [−10, 20], respectively. • S4 0×S3 0 with cubic basis functions for pressure and quartic for velocity, with C0 continuity across the elements, and • S4 2×S3 2 with cubic basis functions for pressure and quartic for velocity, with C2 continuity across the elements. For all these pairs of B-spline discretization spaces, we run the simulations on three different uniform B- spline meshes obtained in consecutive uniform mesh refinement levels, i.e., on B-spline meshes composed of 32× 32, 64× 64 and 128× 128 B-spline elements. The corresponding numbers of degrees of freedom for all these simulations are demonstrated in Table 1. Similarly, we tested different finite element dis- cretization spaces in FEniCS, i.e., • P2/P1 with linear basis functions for pressure and quadratic for velocity, and • P4/P3 with cubic basis functions for pressure and quartic for velocity. P2/P1 P4/P3 32× 32 9539 42691 64× 64 37507 169347 128× 128 148739 674563 256× 256 592387 2692611 Table 2. Degrees of freedom (DOFs) for simulations on different meshes in FEniCS. Let us emphasize that FEniCS is based on triangular meshes. The advantage of FEniCS is that it allows to choose higher degree finite element discretization spaces which can be compared with higher degree B-spline discretization spaces used in our isogeometric fluid flow solver. Further, for these pairs of finite elements discretization spaces, we run the simulations on four different meshes obtained in consecutive re- finement levels. The triangular meshes are composed of 32 × 32, 64 × 64, 128 × 128 and 256 × 256 quads divided into two triangles by a diagonal. The corre- sponding numbers of degrees of freedom for all these simulations are demonstrated in Table 2. In Figure 4, we present stream function and vorticity profiles which were obtained for B-spline discretization space pair S4 2×S3 2 (see [33] for details on stream func- tion and vorticity computation). We do not present stream function and vorticity profiles for other cases in this study because the figures are very similar and presenting these figures would not add value. Let us only note that visually both profiles match the corre- sponding profiles available in the literature very well (see e.g. [54]). Further, we present the comparison of ux-velocity profiles along vertical cross-sections of the computa- tional domain at x = 0.1, x = 0.5, x = 0.95 and uy- velocity profiles along horizontal cross-sections of the computational domain at y = 0.1, y = 0.5, y = 0.99 obtained for selected discretization spaces in study in Figures 5 and 6. Let us emphasize that except the velocity profiles along cross-sections at x = 0.5, y = 0.5, which are quite commonly presented in the literature, we demonstrate the results also for other cross-sections. The reason is that we want to com- pare the solutions obtained for different discretization spaces more thoroughly in the whole domain and the cross-sections near the boundary of the computational domain intersect the problematic parts of the solution. From Figures 5 and 6, we can conclude that most of the solutions obtained for selected discretization space pairs and selected meshes match the reference solution well, or even very well. The best solutions (the closest to the reference solution) are provided by discretization space pairs and meshes S4 2×S3 2 (64×64), P4/P3 (64× 64) and P2/P1 (128× 128). However, let us emphasize that the latter two cases have about four times more degrees of freedom (DOFs) than the former one which indicates the advantage of B-spline discretization spaces and isogeometric analysis in this 40 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis (a). (b). (c). Figure 5. Comparison ux-velocity profiles for vertical cross-sections of the computational domain at x = 0.1 (A), x = 0.5 (B), x = 0.95 (C) obtained for different discretization spaces, where S20S10 stands for S2 0 × S1 0 , S40S30 for S4 0 × S3 0 and S42S32 for S4 2 × S3 2 . 41 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica (a). (b). (c). Figure 6. Comparison uy-velocity profiles for horizontal cross-sections of the computational domain at y = 0.1 (A), y = 0.5 (B), y = 0.99 (C) obtained for different discretization spaces, where S20S10 stands for S2 0 × S1 0 , S40S30 for S4 0 × S3 0 and S42S32 for S4 2 × S3 2 . 42 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis ux at x = 0.1 ux at x = 0.5 uy at y = 0.1 uy at y = 0.5 S2 0×S1 0 (64× 64) 2.1 · 10−2 2.3 · 10−4 4.1 · 10−5 2.9 · 10−4 S2 0×S1 0 SUPG (64× 64) 6.8 · 10−3 3.9 · 10−5 1.2 · 10−5 6.5 · 10−5 S4 0 × S3 0 (32× 32) 2.4 · 10−3 3.4 · 10−5 6.9 · 10−7 1.2 · 10−5 S4 2 × S3 2 (64× 64) 5.1 · 10−6 3.8 · 10−5 6.5 · 10−8 2.4 · 10−6 P2/P1 (64× 64) 2.0 · 10−4 3.1 · 10−4 8.3 · 10−4 4.9 · 10−4 P2/P1 (128× 128) 2.5 · 10−5 1.2 · 10−6 9.0 · 10−8 1.3 · 10−6 P4/P3 (32× 32) 1.2 · 10−4 4.1 · 10−4 2.7 · 10−5 6.4 · 10−4 P4/P3 (64× 64) 3.1 · 10−6 1.7 · 10−6 8.8 · 10−8 1.9 · 10−6 Table 3. Residual sum of squares obtained for studied discretization space pairs and the corresponding reference solution at selected cross-sections through the computational domain. studied case. To have a fair comparison of higher degree finite element discretization space pairs with B-spline discretization space pairs, S4 2 × S3 2 (64× 64) has to be compared with P4/P3 (32×32) because they have the similar number of DOFs (see Table 1 and 2). Looking at these two cases, P4/P3 (32× 32) gives a significantly worse solution (more distant from the reference solution) than S4 2 × S3 2 (64× 64). Further, S4 0 × S3 0 (32× 32), which is similar to the case P4/P3 (32×32) regarding the number of DOFs, degree of dis- cretization spaces and continuity of the basis functions, behaves better than P4/P3 (32× 32) but worse than S4 2 ×S3 2 (64× 64), which has also the similar number of DOFs. To quantify these conclusions, we also com- puted residual sum of squares (RSS) of the obtained solutions and the reference solutions. For selected cross-sections, the results of RSS are demonstrated in Table 3 and they confirm the above mentioned conclusions. The only exception of the discussion presented in the previous paragraph is the uy-velocity profile along a horizontal cross-section at y = 0.99, see Figure 6(C). Near the upper right corner of the computational domain (right part of Figure 6(C)), the solutions ob- tained for B-spline discretization space pairs strongly oscillate. As the above mentioned cross-section is very close to the lid part of the boundary, these oscilla- tions are caused by a sudden change in the boundary condition near the upper corner of the computational domain and also by the fluid flow striking the wall of the cavity. Because of that, oscillations can be expected in this part of the computational domain. One can see that oscillations of the solution for S4 2×S3 2 are smaller in amplitude than for S4 0×S3 0 , but on a wider domain (which is caused by a wider support of B-spline basis functions of high continuity). Note that this behaviour can be also improved and the os- cillations reduced by considering finer computational meshes (at least 128×128 for S4 2×S3 2 ), by local refine- ment near the walls of the cavity or by stabilization of the numerical solution by SUPG method. It is also in- teresting that the FEM solutions obtained by FEniCS are less oscillative for comparable discretization spaces. On the other hand, near the upper left corner of the computational domain (left part of Figure 6(C)) the Initial mesh of S4 2×S3 2 Macro-elements size 64× 64 128× 128 128× 128 – 149258 64× 64 37770 151850 32× 32 39082 157106 16× 16 41778 167906 8× 8 47458 190658 4× 4 59970 240770 2× 2 89602 359426 1× 1 167298 670466 Table 4. Degrees of freedom (DOFs) for simulations in an in-house IgA-based flow solver with rIgA dis- cretization space pairs arising from S4 2×S3 2 . behaviour of the solutions significantly differs from the behaviour near the upper right corner. The solutions obtained by FEniCS oscillate. The solutions obtained for low continuity B-spline discretization space pairs (S2 0×S1 0 , S4 0×S3 0 ) also oscillate but less than FEniCS solutions. Even stabilization of the numerical solu- tion by SUPG method does not help to completely eliminate the oscillations in the case of B-spline dis- cretization space pair S2 0×S1 0 . The best solution in this part is the solution obtained for B-spline dis- cretization space pair S4 2×S3 2 which is without any oscillations and very close to the reference solution. Note that oscillations in the solution in this part of the computational domain are mainly caused by a sudden change in the boundary condition. The above mentioned observations indicate a small stabilization effect of B-spline discretization space pairs of higher continuity. Further, motivated by the recent results on the so- called rIgA (refined Isogeometric Analysis, see [42, 43]) we performed several experiments whether this ap- proach can speed up the time-consuming computa- tions also in the case of the Navier-Stokes equations. The main idea of rIgA approach lies in reduction of the continuity of highly continuous isogeometric discretiza- tion spaces over specific inter-element boundaries (C0 separators). The continuity reduction weakens the interconnection between the subdomains since it re- 43 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica duces the number of basis functions with support over the corresponding inter-element boundary, which im- proves the performance of the factorization using a multifrontal direct solver (see [62]). In our experiments, we start with the high conti- nuity isogeometric discretization space pair S4 2 × S3 2 and successively increase the number of C0 separators (and concurrently decrease the macro-element size). Insertion of C0 separator means inserting a multiple knot into the knot vector determining isogeometric discretization space which leads to the increase of the number of basis functions, and also degrees of freedom. Table 4 shows the numbers of degrees of freedom for rIgA discretization space pairs with decreasing macro- elements size arising from S4 2 × S3 2 considered on two different meshes 64 × 64 and 128 × 128. Note that macro-elements size 1×1 corresponds to discretization space pair S4 0×S3 0 used above. Figure 7 displays speed-up in solving the linear sys- tem (by a direct solver in PARDISO library) obtained from discretization of the steady Navier-Stokes equa- tions for lid-driven cavity problem in one Picard’s iteration. This speed-up is plotted with respect to the macro-elements size and the values on the y-axis correspond to the ratio of time obtained for standard high continuity isogeometric discretization space pair S4 2 × S3 2 and time for rIgA discretization for the par- ticular macro-elements size. All times are averages of 100 Picard’s iterations. For example, in the case of the mesh of 128× 128 elements solving the linear system obtained for macro-elements size 4×4 is 3.76× faster than solving the linear system without any C0 separators. Similarly, Figure 8 shows the reduction in the numbers of floating point operations needed for factorization during solving the linear system obtained by rIgA method for various macro-elements sizes by a direct solver (in PARDISO library) with respect to high continuity isogeometric discretization space pair S4 2 × S3 2 . It can be seen that our results for the steady Navier- Stokes equations are in a good agreement with results presented in [43] for the Stokes equations and they confirm the positive effect of rIgA for solving linear systems also in the case of the steady Navier-Stokes equations. Speed-up indicated by the reduction of floating point operations is not fully reached, as one can see by comparing Figure 7 and Figure 8, but this is also in an agreement with the results presented in [43]. Further, Figures 7 and 8 indicate that one can expect higher speed-up for finer meshes. The reason for this lies in the fact that we are not in a fully asymptotic regime with our computational grids where theoretical speed-up in FLOPs reaches O(p2) for the Stokes problem (see [43]). However, larger grid sizes are too time-consuming because we solve non-linear problem with many iterations required for convergence of the steady Navier-Stokes equations. As an assembly of matrices is considered as one of the major drawbacks of IgA frequently studied in the literature (see e.g. [35–37]), one could ask if increasing the number of degrees of freedom by insertion of C0 separators leads to more time-consuming assembly of matrices for the resulting linear system. If matrices are assembled element-wise, which is still frequently done in IgA-based solvers and which is also done in our case, and if an assembly of matrices only requires evaluation of B-spline basis functions and their deriva- tives for computation of local element matrices and assembling of local matrices to a global matrix (as in the case of mass matrix or stiffness matrix), the cost of an assembly of matrices in high continuity IgA and rIgA is identical. The reason lies in the fact that the number of basis functions per element is the same for high continuity IgA and rIgA. Moreover, during assembling of local matrices to the global matrix, the local matrices highly overlap in high continuity IgA and overlapping of the local matrices reduces in rIgA with decreasing macro-element size. This means that assembly time for rIGA with various macro-elements sizes is at most the same as for high continuity IgA. Figure 9 shows speed-up in assembly of matrices for the linear system in one Picard’s iteration of the steady Navier-Stokes equations discretized by rIgA method with various macro-elements sizes. It is seen that for small macro-elements sizes the time needed for assembly of matrices slightly increases. The rea- sons lies in the fact that we construct the velocity solution in each Picard’s iteration which is used in the construction of the matrix N(u) (see (7)). As the construction of the velocity solution depends on the number of degrees of freedom, which rapidly grows for smaller macro-element sizes (see Table 4), this leads to a small increase in the computational cost of an assembly of matrices. However, this behaviour does not significantly affect the overall effectiveness of the rIgA method. Finally, our experience with solving the Navier- Stokes and RANS (Reynolds-Averaged Navier-Stokes) equations by IgA have shown that decreased conti- nuity to only C0 inside the computational domain can lead to spurious oscillations along this C0 sepa- rator, especially for high Reynolds number fluid flow simulations in complex computational domains. This can be considered as a disadvantage of rIgA method, specifically in its usage for fluid flow simulation. 7. Conclusions In this paper, we have presented numerical results ob- tained by our in-house incompressible fluid flow solver based on isogeometric analysis and we have compared our results with the results obtained by a standard package based on finite element method FEniCS and with the selected reference solutions. For this compari- son, we have chosen the standard benchmark problem for incompressible fluid flow simulation – lid-driven cavity flow. Our experiments on uniform meshes have shown that even for coarse grids and low numbers of degrees 44 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis Figure 7. Time to solve the linear system arising from rIgA discretizations of the Navier-Stokes equation with a direct solver (PARDISO library) presented as speed-up gained by decreasing the macro-elements size. Figure 8. Floating point operations needed for factorization during solving the linear system arising from rIgA discretizations of the Navier-Stokes equation with a direct solver (PARDISO library) presented as speed-up (reduction of FLOPs) gained by decreasing the macro-elements size. Figure 9. Time needed for assembly of matrices in one Picard’s iteration in rIgA discretizations of the Navier-Stokes equation presented as speed-up gained by decreasing the macro-elements size. 45 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica of freedom, isogeometric analysis based on inf-sup stable B-spline discretization space pairs of higher de- gree and continuity can provide results which match reference solutions available in the literature very well. Further, comparing these results with results obtained by FEniCS (based on inf-sup stable finite element dis- cretization space pairs), isogeometric analysis has pro- duced the solutions with higher accuracy for meshes with the approximately same number of degrees of freedom. As the best case of the tested cases, we have identified B-spline discretization space pair S4 2 × S3 2 which produces an accurate solution even for a low number of degrees of freedom (about 40000 DOFs) and which corresponds to the solutions in FEniCS on the meshes with four times more degrees of free- dom. Moreover, the solutions obtained for S4 2 × S3 2 are superior to the solutions obtained for S4 0 × S3 0 for the similar number of degrees of freedom, which con- firms the superiority of B-spline discretization spaces of high continuity. Further, our experiments have indicated a small stabilization effect of high continuity B-spline discretization space pairs. However, this ef- fect is much smaller than the effect of stabilization by SUPG method. Finally, we have shown numerically that the idea of rIgA method is also valid for the steady Navier-Stokes equations, i.e., inserting suitable C0 separators helps to decrease the time needed for solving linear system using a direct solver. Acknowledgements This work was supported by the Czech Science Foundation (GA ČR) grant No. 19-04006S. We thank to all referees for their valuable comments which helped us to improve the paper. References [1] H. Versteeg, W. Malalasekera. An Introduction to Computational Fluid Dynamics - Second Edition. Prentice Hall, 2007. [2] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia, 2003. doi:10.1137/1.9780898718003. [3] H. C. Elman, D. Silvester, A. J. Wathen. Finite Elements and Fast Iterative Solvers With Applications in Incompressible Fluid Dynamics. Oxford University Press, 2014. doi:10.1093/acprof:oso/9780199678792.001.0001. [4] T. Hughes, J. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194(39–41):4135–4195, 2005. doi:10.1016/j.cma.2004.10.008. [5] J. Cottrell, T. Hughes, Y. Bazilevs. Isogeometric analysis: Toward integration of CAD and FEA. John Wiley & Sons, Ltd 2009. [6] L. Piegl, W. Tiller. The NURBS book. Springer Verlag, 1997. [7] A. Falini, C. Giannelli, T. Kanduč, et al. An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes. International Journal for Numerical Methods in Engineering 117(10):1038–1058, 2018. doi:10.1002/nme.5990. [8] G. Farin, D. Hansford. Discrete Coons patches. Computer Aided Geometric Design 16:691–700, 1999. doi:10.1016/s0167-8396(99)00031-x. [9] J. Gravesen, A. Evgrafov, D. Nguyen, P. Nørtoft. Planar parametrization in isogeometric analysis. In International Conference on Mathematical Methods for Curves and Surfaces, pp. 189–212. Springer, 2012. doi:10.1007/978-3-642-54382-1_11. [10] X. Nian, F. Chen. Planar domain parameterization for isogeometric analysis based on Teichmüller mapping. Computer Methods in Applied Mechanics and Engineering 311:41–55, 2016. doi:10.1016/j.cma.2016.07.035. [11] A. Mantzaflaris, B. Jüttler, B. Khoromskij, U. Langer. Low rank tensor methods in Galerkin-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 316:1062–1085, 2017. doi:10.1016/j.cma.2016.11.013. [12] M. Aigner, C. Heinrich, B. Jüttler, et al. Swept volume parameterization for isogeometric analysis. In IMA International Conference on Mathematics of Surfaces, pp. 19–44. Springer, 2009. doi:10.1007/978-3-642-03596-8_2. [13] P. Anděl, B. Bastl, K. Slabá. Parameterizations of generalized NURBS volumes of revolution. Engineering Mechanics 19:293–306, 2012. [14] T. Hughes, A. Reali, G. Sangalli. Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS. Computer Methods in Applied Mechanics and Engineering 197:4104–4124, 2008. doi:10.1016/j.cma.2008.04.006. [15] A. Tagliabue, L. Dedè, A. Quarteroni. Isogeometric Analysis and error estimates for high order partial differential equations in fluid dynamics. Computers & Fluids 102:277–303, 2014. doi:10.1016/j.compfluid.2014.07.002. [16] D. Benson, Y. Bazilevs, M. Hsu, T. Hughes. Isogeometric shell analysis: The Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering 199(5):276 – 289, 2010. Computational Geometry and Analysis, doi:10.1016/j.cma.2009.05.011. [17] F. Auricchio, L. B. da Veiga, A. Buffa, et al. A fully “locking-free” isogeometric approach for plane linear elasticity problems: A stream function formulation. Computer Methods in Applied Mechanics and Engineering 197(1):160 – 172, 2007. doi:10.1016/j.cma.2007.07.005. [18] A. Buffa, G. Sangalli, R. Vázquez. Isogeometric analysis in electromagnetics: B-splines approximation. Computer Methods in Applied Mechanics and Engineering 199(17):1143 – 1152, 2010. doi:10.1016/j.cma.2009.12.002. [19] H. Gómez, V. M. Calo, Y. Bazilevs, T. J. Hughes. Isogeometric analysis of the Cahn-Hilliard phase-field model. Computer Methods in Applied Mechanics and Engineering 197(49):4333 – 4352, 2008. doi:10.1016/j.cma.2008.05.003. 46 http://dx.doi.org/10.1137/1.9780898718003 http://dx.doi.org/10.1093/acprof:oso/9780199678792.001.0001 http://dx.doi.org/10.1016/j.cma.2004.10.008 http://dx.doi.org/10.1002/nme.5990 http://dx.doi.org/10.1016/s0167-8396(99)00031-x http://dx.doi.org/10.1007/978-3-642-54382-1_11 http://dx.doi.org/10.1016/j.cma.2016.07.035 http://dx.doi.org/10.1016/j.cma.2016.11.013 http://dx.doi.org/10.1007/978-3-642-03596-8_2 http://dx.doi.org/10.1016/j.cma.2008.04.006 http://dx.doi.org/10.1016/j.compfluid.2014.07.002 http://dx.doi.org/10.1016/j.cma.2009.05.011 http://dx.doi.org/10.1016/j.cma.2007.07.005 http://dx.doi.org/10.1016/j.cma.2009.12.002 http://dx.doi.org/10.1016/j.cma.2008.05.003 vol. 61 Special Issue/2021 Numerical simulation of lid-driven cavity flow by isogeometric analysis [20] X. Qian. Full analytical sensitivities in NURBS based isogeometric shape optimization. Computer Methods in Applied Mechanics and Engineering 199(29):2059 – 2071, 2010. doi:10.1016/j.cma.2010.03.005. [21] N. D. Manh, A. Evgrafov, A. R. Gersborg, J. Gravesen. Isogeometric shape optimization of vibrating membranes. Computer Methods in Applied Mechanics and Engineering 200(13):1343 – 1353, 2011. doi:10.1016/j.cma.2010.12.015. [22] W. A. Wall, M. A. Frenzel, C. Cyron. Isogeometric structural shape optimization. Computer Methods in Applied Mechanics and Engineering 197(33):2976 – 2988, 2008. doi:10.1016/j.cma.2008.01.025. [23] Y.-D. Seo, H.-J. Kim, S.-K. Youn. Shape optimization and its extension to topological design based on isogeometric analysis. International Journal of Solids and Structures 47(11):1618 – 1640, 2010. doi:10.1016/j.ijsolstr.2010.03.004. [24] N. Liu, A. E. Jeffers. Isogeometric analysis of laminated composite and functionally graded sandwich plates based on a layerwise displacement theory. Composite Structures 176:143 – 153, 2017. doi:10.1016/j.compstruct.2017.05.037. [25] N. Liu, A. E. Jeffers. A geometrically exact isogeometric Kirchhoff plate: Feature-preserving automatic meshing and C1 rational triangular Bézier spline discretizations. International Journal for Numerical Methods in Engineering 115(3):395–409, 2018. https://onlinelibrary.wiley.com/doi/pdf/ 10.1002/nme.5809 doi:10.1002/nme.5809. [26] N. Liu, X. Ren, J. Lua. An isogeometric continuum shell element for modeling the nonlinear response of functionally graded material structures. Composite Structures 237:111893, 2020. doi:10.1016/j.compstruct.2020.111893. [27] N. Liu, A. E. Jeffers. Adaptive isogeometric analysis in structural frames using a layer-based discretization to model spread of plasticity. Computers & Structures 196:1 – 11, 2018. doi:10.1016/j.compstruc.2017.10.016. [28] N. Liu, P. A. Beata, A. E. Jeffers. A mixed isogeometric analysis and control volume approach for heat transfer analysis of nonuniformly heated plates. Numerical Heat Transfer, Part B: Fundamentals 75(6):347–362, 2019. doi:10.1080/10407790.2019.1627801. [29] N. Liu, A. E. Jeffers. Feature-preserving rational Bézier triangles for isogeometric analysis of higher-order gradient damage models. Computer Methods in Applied Mechanics and Engineering 357:112585, 2019. doi:10.1016/j.cma.2019.112585. [30] B. Bastl, M. Brandner, J. Egermaier, et al. IgA-based solver for turbulence modelling on multipatch geometries. Advances in Engineering Software 113:7–18, 2017. doi:10.1016/j.advengsoft.2017.06.012. [31] B. Bastl, M. Brandner, J. Egermaier, et al. Isogeometric analysis for turbulent flow. Mathematics and Computers in Simulation 145:3–17, 2018. doi:10.1016/j.matcom.2016.05.010. [32] Y. Bazilevs, V. Calo, J. Cottrell, et al. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Computer methods in applied mechanics and engineering 197(1-4):173–201, 2007. doi:10.1016/j.cma.2007.07.016. [33] S. Hosseini, M. Möller, S. Turek. Isogeometric Analysis of the Navier-Stokes equations with Taylor-Hood B-spline elements. Appl Math Comput 267:264 – 281, 2015. doi:10.1016/j.amc.2015.03.104. [34] A. Falini, J. Špeh, B. Jüttler. Planar domain parameterization with THB-splines. Computer Aided Geometric Design 35–36:95–108, 2015. doi:10.1016/j.cagd.2015.03.014. [35] A. Mantzaflaris, B. Jüttler, B. N. Khoromskij, U. Langer. Matrix generation in isogeometric analysis by low rank tensor approximation. In International Conference on Curves and Surfaces, pp. 321–340. 2014. doi:10.1007/978-3-319-22804-4_24. [36] M. Bartoň, V. M. Calo. Gauss-Galerkin quadrature rules for quadratic and cubic spline spaces and their application to isogeometric analysis. Computer-Aided Design 82:57–67, 2017. doi:10.1016/j.cad.2016.07.003. [37] R. R. Hiemstra, G. Sangalli, M. Tani, et al. Fast formation and assembly of finite element matrices with application to isogeometric linear elasticity. Computer Methods in Applied Mechanics and Engineering 355:234–260, 2019. doi:10.1016/j.cma.2019.06.020. [38] N. Collier, D. Pardo, L. Dalcin, et al. The cost of continuity: A study of the performance of isogeometric finite elements using direct solvers. Comput Methods Appl Mech Engrg 213-216:353–361, 2012. doi:10.1016/j.cma.2011.11.002. [39] T. A. AbdelMigid, K. M. Saqr, M. A. Kotb, A. A. Aboelfarag. Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using gpu accelerated code. Alexandria Engineering Journal 56:123–135, 2017. doi:10.1016/j.aej.2016.09.013. [40] J. A. Evans, T. J. R. Hughes. Isogeometric divergence-conforming B-splines for the steady Navier–Stokes equations. Math Models Methods Appl Sci 23(8):1421–1478, 2013. doi:10.21236/ada560496. [41] V. John, L. Schumacher. A study of isogeometric analysis for scalar convection–diffusion equations. Applied Mathematics Letters 27:43–48, 2014. doi:10.1016/j.aml.2013.08.004. [42] D. Garcia, D. Pardo, L. Dalcin, et al. The value of continuity: Refined isogeometric analysis and fast direct solvers. Computer Methods in Applied Mechanics and Engineering 316:586–605, 2017. doi:10.1016/j.cma.2016.08.017. [43] D. Garcia, D. Pardo, V. M. Calo. Refined isogeometric analysis for fluid mechanics and electromagnetics. Computer Methods in Applied Mechanics and Engineering 356:598–628, 2019. doi:10.1016/j.cma.2019.06.011. [44] F. Brezzi, M. Fortin. Mixed and Hybrid Finite Element Methods. Springer-Verlag, 1991. 47 http://dx.doi.org/10.1016/j.cma.2010.03.005 http://dx.doi.org/10.1016/j.cma.2010.12.015 http://dx.doi.org/10.1016/j.cma.2008.01.025 http://dx.doi.org/10.1016/j.ijsolstr.2010.03.004 http://dx.doi.org/10.1016/j.compstruct.2017.05.037 https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5809 https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5809 http://dx.doi.org/10.1002/nme.5809 http://dx.doi.org/10.1016/j.compstruct.2020.111893 http://dx.doi.org/10.1016/j.compstruc.2017.10.016 http://dx.doi.org/10.1080/10407790.2019.1627801 http://dx.doi.org/10.1016/j.cma.2019.112585 http://dx.doi.org/10.1016/j.advengsoft.2017.06.012 http://dx.doi.org/10.1016/j.matcom.2016.05.010 http://dx.doi.org/10.1016/j.cma.2007.07.016 http://dx.doi.org/10.1016/j.amc.2015.03.104 http://dx.doi.org/10.1016/j.cagd.2015.03.014 http://dx.doi.org/10.1007/978-3-319-22804-4_24 http://dx.doi.org/10.1016/j.cad.2016.07.003 http://dx.doi.org/10.1016/j.cma.2019.06.020 http://dx.doi.org/10.1016/j.cma.2011.11.002 http://dx.doi.org/10.1016/j.aej.2016.09.013 http://dx.doi.org/10.21236/ada560496 http://dx.doi.org/10.1016/j.aml.2013.08.004 http://dx.doi.org/10.1016/j.cma.2016.08.017 http://dx.doi.org/10.1016/j.cma.2019.06.011 B. Bastl, M. Brandner, J. Egermaier et al. Acta Polytechnica [45] A. Bressan, G. Sangalli. Isogeometric discretizations of the Stokes problem: stability analysis by the macroelement technique. IMA Journal of Numerical Analysis 33(2):629–651, 2012. doi:10.1093/imanum/drr056. [46] P. N. Nielsen, A. R. Gersborg, J. Gravesen, N. L. Pedersen. Discretizations in isogeometric analysis of Navier–Stokes flow. Computer Methods in Applied Mechanics and Engineering 200(45):3242–3253, 2011. doi:10.1016/j.cma.2011.06.007. [47] V. John, P. Knobloch. On spurious oscillations at layers diminishing (SOLD) methods for convection- diffusion equations: Part I - A review. Computer Methods in Applied Mechanics and Engineering 196(17):2197 – 2215, 2007. doi:10.1016/j.cma.2006.11.013. [48] A. N. Brooks, T. J. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Computer methods in applied mechanics and engineering 32(1- 3):199–259, 1982. doi:10.1016/0045-7825(82)90071-8. [49] K. Takizawa, T. E. Tezduyar, Y. Otoguro. Stabilization and discontinuity-capturing parameters for space–time flow computations with finite element and isogeometric discretizations. Computational Mechanics 62(5):1169–1186, 2018. doi:10.1007/s00466-018-1557-x. [50] Y. Chai, J. Ouyang. Appropriate stabilized Galerkin approaches for solving two-dimensional coupled Burgers’ equations at high Reynolds numbers. Computers & Mathematics with Applications 79(5):1287–1301, 2019. doi:10.1016/j.camwa.2019.08.036. [51] V. John, E. Schmeyer. Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion. Computer methods in applied mechanics and engineering 198(3-4):475–494, 2008. doi:10.1016/j.cma.2008.08.016. [52] U. Ghia, K. Ghia, C. Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics 48(3):387–411, 1982. doi:10.1016/0021-9991(82)90058-4. [53] O. Botella, R. Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids 27(4):421– 433, 1998. doi:10.1016/s0045-7930(98)00002-4. [54] E. Erturk, T. C. Corke, C. Gökçöl. Numerical solutions of 2-d steady incompressible driven cavity flow at high reynolds numbers. International Journal for Numerical Methods in Fluids 48(7):747–774, 2005. doi:10.1002/fld.953. [55] Y. Papadopoulos. A driven cavity exploration. [Online; accessed 30-Jan-2020]. [56] C.-H. Bruneau, M. Saad. The 2D lid-driven cavity problem revisited. Computers & Fluids 35:326–348, 2006. doi:10.1016/j.compfluid.2004.12.004. [57] E. Erturk. Discussions on driven cavity flow. International Journal for Numerical Methods in Fluids 60(3):275–294, 2009. doi:10.1002/fld.1887. [58] Z. Cai, Y. Wang. An error estimate for two-dimensional stokes driven cavity flow. Mathematics of computation 78(266):771–787, 2009. [59] P. Matuszyk, M. Paszyński. Fully automatic hp adaptive finite element method for the Stokes problem in two dimensions. Computer Methods in Applied Mechanics and Engineering 197(51-52):4549–4558, 2008. doi:10.1016/j.cma.2008.05.027. [60] M. S. Alnaes, J. Blechta, J. Hake, et al. The FEniCS Project Version 1.5. Archive of Numerical Software 3:9–23, 2015. [61] B. Jüttler, U. Langer, A. Mantzaflaris, et al. Geometry + simulation modules: Implementing isogeometric analysis. In PAMM – Proceedings of Applied Mathematics and Mechanics, vol. 14, pp. 961–962. 2014. doi:10.1002/pamm.201410461. [62] I. S. Duff, J. K. Reid. The multifrontal solution of indefinite sparse symmetric linear. ACM Trans Math Softw 9(3):302–325, 1983. doi:10.1145/356044.356047. 48 http://dx.doi.org/10.1093/imanum/drr056 http://dx.doi.org/10.1016/j.cma.2011.06.007 http://dx.doi.org/10.1016/j.cma.2006.11.013 http://dx.doi.org/10.1016/0045-7825(82)90071-8 http://dx.doi.org/10.1007/s00466-018-1557-x http://dx.doi.org/10.1016/j.camwa.2019.08.036 http://dx.doi.org/10.1016/j.cma.2008.08.016 http://dx.doi.org/10.1016/0021-9991(82)90058-4 http://dx.doi.org/10.1016/s0045-7930(98)00002-4 http://dx.doi.org/10.1002/fld.953 http://dx.doi.org/10.1016/j.compfluid.2004.12.004 http://dx.doi.org/10.1002/fld.1887 http://dx.doi.org/10.1016/j.cma.2008.05.027 http://dx.doi.org/10.1002/pamm.201410461 http://dx.doi.org/10.1145/356044.356047 Acta Polytechnica 61(SI):33–48, 2021 1 Introduction 2 Navier-Stokes equations 2.1 Galerkin discretization 2.2 LBB condition 3 Isogeometric analysis 4 Streamline Upwind/Petrov-Galerkin (SUPG) stabilization 5 Lid-Driven Cavity Flow 6 Results 7 Conclusions Acknowledgements References