IgA-Based Solver for turbulence modelling on multipatch geometries

ARTICLE IN PRESS

JID: ADES

[m5G;July 5, 2017;14:58]

Advances in Engineering Software 0 0 0 (2017) 1–12

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

Research paper

IgA-Based Solver for turbulence modelling on multipatch geometries Bohumír Bastl∗, Marek Brandner, Jirˇí Egermaier, Kristýna Michálková, Eva Turnerová ˇ Czech Republic NTIS – New Technologies for the Information Society, Faculty of Applied Sciences, University of West Bohemia, Univerzitní 8, 301 00 Plzen,

a r t i c l e

i n f o

Article history: Received 29 February 2016 Revised 22 May 2017 Accepted 18 June 2017 Available online xxx Keywords: NURBS Objects Isogeometric analysis Navier–Stokes equations RANS Equations k − ω Model

a b s t r a c t This paper is focused on numerical solving of RANS (Reynolds-Averaged Navier-Stokes) equation with k − ω model for simulation of turbulent flows in 3D. The solver which is based on a recently proposed approach called isogeometric analysis is presented. This numerical method is based on isoparametric approach, i.e., the same basis functions are used for the description of a geometry of a computational domain and also for the representation of a solution. As computational domains are described by NURBS objects in isogeometric analysis, any real application requires to handle the so-called multipatch domains, where the computational domain is composed of more parts and each part is represented by one NURBS object. In our solver, discontinuous Galerkin method is used to connect different NURBS patches into one computational domain. The results of the solver are demonstrated on a standard benchmark example – backward facing step.

1. Introduction The development of the numerical methods to simulate fluid flows is fundamental in practice. The goal is to find numerical method, which is efficient, and the numerical solution of the given problem, which is the most accurate. The classical numerical methods, which are used to find approximate solution of incompressible fluid flow, are finite difference methods, finite volume methods, finite element methods, spectral methods etc. Finite difference method (FDM) [1] is the oldest one and it is simple to implement, but it is rarely applied for fluid flow simulation because it is difficult to implement this method in the case of complex geometries. On the other hand, finite volume method (FVM) is a common approach and it is frequently used in CFD codes because it follows continuous conservation laws of physics, see e.g. [1,2]. FVM is also easy to implement, it requires only low memory, even in case of turbulent simulation where finer grids are used. Further, FVM is robust, flexible and provides a reasonable numerical solution even on coarser grid. One of drawbacks is that too much diffusion can be added using methods of lower order. In the case of high-order FVM the size of reconstruction stencils is increasing and this approach becomes computationally expensive [3]. Finally, the finite element method (FEM) [1,4,5] can be applied for solving fluid

∗

Corresponding author. . E-mail addresses: [email protected] (B. Bastl), [email protected] (M. Brandner), [email protected] (J. Egermaier), [email protected] (K. Michálková), [email protected] (E. Turnerová).

© 2017 Elsevier Ltd. All rights reserved.

flow simulation. However, this method is not very popular in CFD community because of higher memory requirements [6]. Nevertheless, one of its advantages is that this method is able to deal relatively easily with problems that are defined on complex geometries and provide higher accuracy even on coarse grid compared to FVM and FDM methods. FEM method is one of the Galerkin approaches which are based on weak formulation of the problem and polynomial approximation of the solution is usually used. To avoid oscillations in the convection-dominated case a suitable stabilization has to be used (SUPG, SOLD, GLS, etc.). Robust, efficient and high-order stabilized methods are not available now and there is still urgent need to construct better non-linear stabilization techniques [7]. Recently, modification of FEM based on B-spline/NURBS objects appeared (cf. [8,9]). This approach is known as isogeometric analysis and it is similar to FEM, only triangular/tetrahedral meshes typically used in FEM are replaced by meshes composed of parts of NURBS surfaces/volumes representing a computational domain. This is more suitable for consequent computations, because it allows to avoid the time-consuming step of generating triangular/tetrahedral meshes and directly perform computations. Moreover, since the discretization of a computational domain is always exact, this approach reduces errors in the computational analysis and is more suitable for formulation of automatic shape optimization algorithms. Turbulent flow is extensively studied in engineering practice in these days. Big effort is made to produce numerical methods which describe turbulent flow precisely as much as possible. Moreover,

http://dx.doi.org/10.1016/j.advengsoft.2017.06.012 0965-9978/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

JID: ADES 2

ARTICLE IN PRESS

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

the numerical methods have to be able to cover a wide range of applications for both, the external and internal flows. The most common approaches are DNS (Direct Numerical Simulation), RANS (Reynolds-Averaged Navier-Stokes) and LES (Large Eddy Simulation). Using DNS approach, all details are simulated till the smallest scales of the flow [2,10]. But the number of grid points in spatial discretization must be proportional to Re9/4 , where Re is the Reynolds number, otherwise the simulation becomes unstable. Moreover, in case of unsteady problem, the time step has to be sufficiently small to resolve the movement of the fastest fluctuations. Thus, the direct numerical simulation is applied only for low Reynolds number flows and simple geometries. In the case of a fully developed turbulent flow, certain simplifications need to be involved. In case of LES method, only the large scales (large energetic structures) are simulated and the small scales are modeled (cf. [2,11]). Thus, it is less expensive than DNS and can be applied to higher Reynolds number flows with sufficient accuracy. RANS approach represent another simplification and it is the most common method for simulation of turbulent flows. Because all turbulent scales are modeled, it has low memory requirements (cf. [2,12]). RANS equations are time-averaged equations of fluid flow motion. Similarly to LES, flow variables (velocity and pressure) are separated (the so-called Reynolds decomposition). But in case of RANS, the variables are decomposed into a mean (timeaveraged) component and a fluctuating component and the NavierStokes equations are integrated over a time interval. The attention is focused only on the mean flow and the effects of the turbulence on the mean flow, while the turbulent fluctuations are not resolved in details which is often sufficient in practice. Time-averaging (or ensemble/Reynolds averaging) causes that extra unknown term appears (Reynolds stress) which must be approximated using special turbulent models, where the Boussinesq assumptions are usually applied (cf. [2,12–15]). Turbulence models differ in the way of approximation of the Reynolds stress term, their applicability and computational requirements. In modern engineering applications, the commonly used turbulent models are Spalart-Allmaras, k −  , k − ω and SST models (cf. [2,16,17]). The main goal of this paper is to combine RANS approach to modeling turbulent flows with a novel method, called isogeometric analysis (IgA). Moreover, computational domains are considered to be composed of more NURBS patches, the so-called multipatch domains in isogeometric analysis terminology. As the main purpose of the solver, presented in this paper, is a fluid flow simulation in water turbines, the most general case of non-conforming NURBS meshes (non-conforming multipatch NURBS domains) is considered. Thus, the solver combines continuous IgA-Galerkin method on different NURBS patches with discontinuous IgA-Galerkin approach on interfaces between NURBS patches. This is still a challenging task in current research, as usually conforming or nested meshes (see Section 3.2) are considered and/or only linear problems are solved (typically, elliptic problems), cf. e.g. [18,19]. The paper is organized as follows. Section 2 very briefly recalls basic facts about NURBS objects and Navier-Stokes equations. In Section 3, theoretical aspects of the solution of NavierStokes equations and RANS equations with k − ω model, based on isogeometric analysis on multipatch NURBS domains, are presented. Section 4 then presents several implementation details mainly related to the usage of isogeometric analysis as the numerical method for solving these problems and also to discontinuous IgA-Galerkin approach for handling multipatch domains. Section 5 is devoted to the presentation of results of the solver. Finally, Section 6 concludes the paper.

Fig. 1. B-spline basis functions.

2. Preliminaries 2.1. NURBS Objects NURBS surface (cf. [8,9]) of degree r, s is determined by a control net P (of control points Pi, j , i = 0, . . . , n, j = 0, . . . , m), weights wi, j of these control points and two knot vectors U = (u0 , . . . , un+r+1 ), V = (v0 , . . . , vm+s+1 ) and it is given by a parametrization n  m 

S (u, v ) =

i=0 j=0

wi, j P i, j Ni,r (u )M j,s (v )

n  m  i=0 j=0

wi, j Ni,r (u )M j,s (v )

=

n  m 

P i, j Ri, j (u, v ).

(1)

i=0 j=0

B-spline basis functions Ni, r (u) and Mj, s (v) are determined by knot vectors U and V and degrees r and s, respectively, by a formula (Mj, s (v) is constructed in the similar way)



Ni,0 (u ) = Ni,r (u ) =

1 ui ≤ u < ui+1 0 otherwise

(2)

u − ui ui+r+1 − u Ni,r−1 (u ) + Ni+1,r−1 (u ). ui+r − ui ui+r+1 − ui+1

Knot vector is a non-decreasing sequence of real numbers which determines the distribution of a parameter on the corresponding curve/surface. B-spline basis functions (see Fig. 1) of degree r are C r−1 -continuous in general. Knot repeated k times in the knot vector decreases the continuity of B-spline basis functions by k − 1. Support of B-spline basis functions is local – it is nonzero only on the interval [ti , ti+r+1 ] in the parameter space and each B-spline basis function is non-negative, i.e., Ni, r (t) ≥ 0, ∀t. 2.2. Navier–Stokes equations The model of viscous flow of an incompressible Newtonian fluid can be described by the Navier-Stokes equations in the common form (cf. [4–6])

∂u + ∇ p + u · ∇ u − ν u = f , in  × (0, T ), ∂t ∇ · u = 0, in  × (0, T ),

(3)

where  ⊂ Rd (dimension d = 1, 2 or 3) is the computational domain, T > 0 is the final time, u = u(x, t ) is the vector function describing flow velocity, p = p(x, t ) is the kinematic pressure function, ν describes kinematic viscosity and f additional body forces acting on the fluid. The laminar/turbulent behaviour of the fluid can be predicted by the dimensionless Reynolds number generally defined as the ratio of momentum forces to viscous forces

Re =

uˆ · L

ν

,

(4)

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

ARTICLE IN PRESS

JID: ADES

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

where uˆ is the mean velocity of the fluid and L is the characteristic linear dimension. Generally, this model can be used for any Reynolds number, but it may require very fine meshes and very small time steps to provide a reasonable solution. The initial-boundary value problem is considered as the system (3) together with suitable initial conditions and the following boundary conditions

u=w

∂u ν − np = 0 ∂n

on

∂ D × 0, T (Dirichlet b. c.),

on

∂ N × 0, T (Neumann b. c.).

(5)

computational domain. In 3D, the solution has the form u

uh =

nd 

3. Flow analysis with isogeometric analysis This section is focused on solving Navier-Stokes equations on singlepatch and multipatch NURBS domains with the help of isogeometric analysis. Singlepatch or multipatch NURBS domain are meant to be computational domains described by one NURBS object (surface in 2D or volume in 3D) or composed of several NURBS objects joined along common interfaces, respectively.

3.1. Navier–Stokes equations on singlepatch domains

i=1

V = {u ∈ H 1 ()d |u = w on





−  



 

p∇ · vd =

∇ u : ∇ vd + 



q ∇ · ud  = 0 ,



vh ·

−



∂ uh d + ν ∂t



fh =

nv 







ph ∇ · vh d =

⎡

A+C ⎢ 0 ⎣ 0 B1 C ⎢0 +⎣ 0 0

f h · vh d,

0 A+C 0 B2

0 0 A+C B3

⎤

⎡

⎤⎡

 

A = Ai j

1≤i≤nud ,1≤ j≤nud



Bm = Bmi j



 

C = Ci j



f m = t C |C

  ∗

· fmi



⎡⎤

(11)

⎤

∗ 0 u1 0 ⎥ u∗2 ∗ ∗⎦ A +C u∗3 B∗3

0 A + C∗ 0 B∗2 ∗

 

A∗ = A i j

,

1≤i≤n p ,1≤ j≤nud

1≤i≤nud ,1≤ j≤nud

⎤

−BT1 un1+1 f1 −BT2 ⎥⎢un2+1 ⎥ ⎢f 2⎥ ⎦⎣un+1 ⎦ = ⎣f ⎦+ −BT3 3 3 0 pn+1 0

0 n A∗ + C ∗ u1 0⎥ n ⎢ 0 u −⎣ C⎦ 2n 0 u3 0 B∗1

0 C 0 0



(uh · ∇ uh )vh d  qh ∇ · uh d = 0 

Ai j = ν t

Ci j = (7)

∇ uh : ∇ vh d + 

(10)

For general function f, the approximation fh can be obtained with the help of L2 projection into a linear space spanned by basis functions {Rui }1≤i≤nuv . Then the discrete problem can be written (implicit Euler method is used for the time discretization, superscripts n and n + 1 denote the values at time layers tn and tn+1 ) in the matrix form as



f · vd, ∀v ∈ V0 ,



( f1i , f2i , f3i )T Rui .

i=1

Bmi j = t

(u · ∇ u )vd

∀q ∈ L2 ().



i=1

,

B∗m

1≤i≤nud ,nud +1≤ j≤nuv



= Bmi j

 

C ∗ = Ci j

,



,

,

1≤i≤n p ,nud +1≤ j≤nuv

1≤i≤nud ,nud +1≤ j≤nuv

,

,

1≤i≤nuv

(12)

Further, Galerkin approach using NURBS basis functions is performed. Let us define the finite dimensional subspaces V h ⊂ V, V0h ⊂ V0 , Wh ⊂ L2 () then the discrete weak solutions uh ∈ Vh and ph ∈ Wh are searched such that



pi Rip ,

(6)

Then a weak formulation of the initial-boundary value problem is find u ∈ V and p ∈ L2 () [4] such that

∂u d + ν ∂t

i=nud +1

where nud is the number of points where the Dirichlet boundary condition is not defined, nuv = nu + nud . Further, f is written as a linear combination of velocity basis functions, i.e.,

where 3.1.1. Galerkin approach For simplicity, only Stokes equations (convective term is removed in (3)) will be considered in this section. The convective term is discussed separately in the next section. Let V be a velocity solution space and V0 be the corresponding space of test functions, i.e.,

v·

p

n 

(u∗1i , u∗2i , u∗3i )T Rui , ph =

(9)

⎡

This section is devoted to the solution of the Navier–Stokes problem (3), (5) on a domain described by a single NURBS patch [8].



u

nv 

(u1i , u2i , u3i )T Rui +

u

If the velocity is specified everywhere on the boundary, the pressure solution is only unique up to a hydrostatic constant.

∂ D } , V0 = {v ∈ H 1 ()d |v = 0 on ∂ D }.

3

(8)

for all test functions vh ∈ V0h and qh ∈ Wh . Similarly to Galerkin approach for classical finite element method, the approximated solution uh , ph respectively, is given by linear combination of the nu , p p np basis functions Rui ∈ V h , Ri ∈ W h respectively, where Rui and Ri are NURBS basis function obtained from a NURBS description of a



 

 

(∇ Rui · J−1 ) · (∇ Ruj · J−1 )| det J|d,



Rip (∇ Ruj · J −1 ) · em



| det J|d,

(13)

Rui Ruj | det J|d,

m = {1, 2, 3} and J is Jacobi matrix of a mapping from the parametric domain to the computational domain. 3.1.2. Nonlinear iteration Because of non-linearity of Navier–Stokes equations, it is necessary to solve the problem iteratively with linear problem in every step. One of the possibilities is to use the so-called Picard’s method [7]. For the sake of simplicity, stationary problem is mentioned in the following paragraphs, as time derivative does not affect explanation of this method. First, a non-linear residuum of the weak formulation is computed by the values uk and pk . The residuum Rk and rk satisfy



Rk =





+

f · vd − ν 





pk ∇ · vd,

∇ uk : ∇ vd − 

rk = −







(uk · ∇ uk )vd +

q ∇ · u k d .

(14)

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

ARTICLE IN PRESS

JID: ADES 4

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

An exact solution u and p is considered to be the sum of the solution of the current iteration and fluctuation of the solution, i.e.,

u = uk + δ u k , p = pk + δ pk .

(15)

Substituting (15) into the weak formulation and using the equalities (14)





f · vd − ν



∇ δ uk : ∇ vd −



(uk · ∇δ uk )vd −      − (δ uk · ∇ uk )vd − (δ uk · ∇δ uk )vd + δ pk ∇ · vd,    k k r =− q ∇ · δ u d . (16)

Rk =



 The quadratic term  (δ uk · ∇δ uk )vd and also the linear term  k k  (δ u · ∇ u )vd are assumed to be sufficiently small and can be neglected. The linear problem (16) is obtained for fluctuations δ uk ∈ V0 and δ pk ∈ L2 (), which define the next step uk+1 = uk + δ uk , pk+1 = pk + δ pk . Therefore, uk+1 ∈ V and pk+1 ∈ L2 () are searched such that the following relation is valid

ν

 

−

∇ uk+1 : ∇ vd +







pk+1 ∇ · vd =





(uk · ∇ uk+1 )vd −



f · vd,





q∇ · uk+1 d = 0

(17)

for all functions v ∈ V0 and q ∈ L2 (). The solution of the Stokes problem is used to be the initial condition (u0 , p0 ) of the Picard’s iteration (17). According to the mentioned Galerkin approach above, the finite dimensional spaces V h ⊂ V, V0h ⊂ V0 , Wh ⊂ L2 () are defined. Let us assume that ukh ∈ V h and pkh ∈ W h are the numerical solutions of the previous iteration step, then the approximations ukh+1 and pkh+1 are searched such that all functions vh ∈

ν

 

−

∇ ukh+1 : ∇ vh d +







pkh+1 ∇ · vh d =





V0h

and qh ∈

Wh

satisfy

A + N ( uk ) ⎢ 0 ⎣ 0 B1



0 A + N ( uk ) 0 B2

⎡

f · vh d,





qh ∇ · ukh+1 d = 0. (18)

⎤⎡

⎤

uk1+1 −BT1 T −B2 ⎥⎢uk2+1 ⎥ ⎦⎣uk+1 ⎦ = −BT3 3 pk+1 0

0 0 A + N ( uk ) B3

⎤

(19)

where A, A∗ , Bm , B∗m are defined by (12) and





Rui



u

nv 

(

)

u1l , u2l , u3l Rul

· (∇

Ruj

·J

−1

 1≤i≤nud ,1≤ j≤nud

,



N ∗ (u ) = Ni j (u )

) | det J|d,

l=1

(21)

f 1 − (A∗ + N ∗ (uk )) · u∗1 ⎢ f 2 − (A∗ + N ∗ (uk )) · u∗2 ⎥ =⎣ ⎦, f 3 − (A∗ + N ∗ (uk )) · u∗3 ∗ ∗ ∗ ∗ ∗ ∗ B1 · u1 + B2 · u2 + B3 · u3

N (u ) = Ni j (u )

Ni j (u ) =





(ukh · ∇ ukh+1 )vh d −

The solution ukh , pkh is again written as a linear combination of NURBS basis functions (see (9)) and it is substituted into the system (18). The sequence of the solutions (ukh , pkh ) ∈ V h × W h converges to the weak solution. The system has the matrix form

⎡

Fig. 2. Examples of multipatch domains with different type of connection. Top: Conforming NURBS meshes; Middle: Nested NURBS meshes; Bottom: Nonconforming NURBS meshes.

3.2. Navier–Stokes equations on multipatch domains Theory of NURBS objects directly implies that it is not possible (or reasonable, in some cases) to describe an object of arbitrary topology by one NURBS object [9]. One of the typical examples is a rectangular plate with two holes. The reason for this lies in a regular control net describing NURBS object, which is composed of m × n (in 2D) and l × m × n (in 3D) control points. Thus, if the isogeometric analysis is used for numerical solving of partial differential equations, it is usually necessary to decompose a computational domain into subdomains, which are suitable for description by one NURBS object. Such a domain is then composed of multiple NURBS patches and will be called multipatch domain in the following. Then, any solver based on isogeometric analysis working on multipatch domains has to be able to join NURBS patches along their interfaces into one computational domain. Generally, multipatch domains are usually distinguished by the type of connection of patches into: •

 1≤i≤nud ,nud +1≤ j≤nuv

,

(20)

•

conforming NURBS meshes – patches have along a common interface the same elements (discretization) and the same control nets describing this interface (see Fig. 2 (top)), nested NURBS meshes – discretization of the common interface of left patch is obtained as a uniform refinement of a dis-

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

ARTICLE IN PRESS

JID: ADES

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

•

cretization of the interface of the right patch, or vice versa (see Fig. 2 (middle)), non-conforming NURBS meshes – description of the common interface between patches is completely independent on both patches (see Fig. 2 (bottom)).

Conforming meshes can be handled easily – it is enough to identify the corresponding control points determining the common interface on both patches. The case of nested meshes is more complex, but it is still possible to find a direct relation between control points determining a common interface on both patches. The most complex case is then represented by non-conforming meshes which is usually handled differently.

( u, v ) = (σ , τ ) =

¯ = 

nel 

¯ i, 

i ∩  j = ∅, i = j,

i=1

ˆ := 

nel 

i ,

=

i=1

nel 

(22)

¯ i∩ ¯j= 

i, j=1,i = j

nel 

∂ i \ ∂ .

(23)

 

u · vd,

σ : τ d ,

for vector functions,

(28)

for second-order tensors.

Analogously,

( p, q )ϒ =



ϒ

pqdϒ

(29)

denotes L2 scalar product of scalars p, q in any domain Y ⊂ ∪ ∂ . The standard Interior Penalty Method (IPM) [20] is used for the Navier-Stokes problem, i.e. find uh ∈ Vh and ph ∈ Wh such that ∀v ∈ Vh and ∀q ∈ Wh

 3.2.1. Discontinuos Galerkin method One of the possible approaches to multipatch domains described by non-conforming meshes is to use discontinuous Galerkin method. The Application of the discontinuous Galerkin methods for solving the Navier–Stokes equations (based on the classical FEM) is described in [20] in details. The principles are briefly described in the following. Let  ⊂ R3 be an open bounded domain with the boundary ∂ . Suppose that  is partitioned in nel disjoint subdomains i with boundaries ∂ i , that define an internal interface . The following notation will be used



5

 ∂ uh , v + aIP (uh , v ) + c (uh ; uh , v ) − ∂t −b(v, ph ) + ({ ph }, n · v ) = lIP (v ), b(uh , q ) − ({q}, n · uh  ) = 0,

(30)

where

aIP (u, v ) = (ν u, ∇ v ) +

γ h

(n  u, n  v ) − (ν{∇ u},

n  v ) − (n  u, ν{∇ v} ) , lIP (v ) = ( f , v ),

b( v , p ) =





q∇ · vd

(31)

and the trilinear form associated to the convective term (u · ∇ )u is defined as

c (u; w, v ) = ((u · ∇ )w, v ) + (

i=1

1 1 {u · n}w − |{u · n}|w, v ) . 2 2 (32)

Remark 1. Let us emphasize that in the context of FEM (cf. [20]), subdomains i are represented by elements in discretization, but in our context these subdomains represent NURBS patches composing a computational domain, i.e. discontinuous Galerkin method is used only to connect NURBS patches along common interfaces.

A positive scalar γ , the so-called penalty parameter, must be large enough to ensure the coercivity of the symmetric bilinear form aIP . The characteristic mesh size is denoted by h. Construction of the matrix form of this problem is described in the Section 4.

The incompressible Navier–Stokes equations are considered in the form (3) with boundary conditions (5), where

3.3. RANS equations on multipatch domains

∂  = ∂ D ∪ ∂ N ,

The Navier-Stokes equations can describe turbulent incompressible flow. The turbulent flow contains lots of eddies of different sizes which are changing in time. RANS equations are the most common model for simulation of the turbulent flows. It is based on decomposition of the solution into a time-averaged value and a fluctuation value. In 3D, the solution u = u(x1 , x2 , x3 , t ) and p = p(x1 , x2 , x3 , t ) is decomposed in the following way

∂ D ∩ ∂ N = ∅

(24)

and interface conditions  n  u

∂u ν − n p ∂n

=

0,

on ,

=

0,

on .

(25)

Let us define the jump · and mean {·} operators along the interface using values from the elements to the left and to the right of the interface, namely

  =

i +  j 

on , on ∂ ,

 {} =

1 2



i + 12  j

on , on ∂ . (26)

For instance, in the case of two neighboring subdomains i and j , whose exterior unit normals are denoted ni and nj (evidently ni = −n j ), and n denotes exterior unit normal along ∂ , then



p i ni + p j n j = ni ( p i − p j ) pn

 pn =

on , on ∂ 

(27)

for scalar function p. L2 scalar product in  is defined as

( p, q ) =





pqd,

for scalar functions,

u(x1 , x2 , x3 , t ) = u¯ (x1 , x2 , x3 ) + u (x1 , x2 , x3 , t ),

(33)

p(x1 , x2 , x3 , t ) = p¯ (x1 , x2 , x3 ) + p (x1 , x2 , x3 , t ),

(34)

where u¯ , p¯ are time-averaged values and u , p are fluctuation ones. Substituting (33) and (34) into the system (3) leads to (see [2] for details)

∂ u¯ + u¯ · ∇ u¯ = −∇ p¯ + ν u¯ − u · ∇ u , ∇ · u¯ = 0. ∂t

(35)

The attention is focused only on the mean flow u¯ , p¯ and the effects of the turbulence on the mean flow. While the turbulent fluctuations u , p are not resolved anymore. Obviously, an extra term appears in the RANS equations which is unknown, the so called Reynolds stress

u · ∇ u  .

(36)

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

ARTICLE IN PRESS

JID: ADES 6

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

This stress term requires modeling to close the RANS equations. Thus, the so called Boussinesq hypothesis is applied, i.e. the following assumption is considered, see more e.g. in [12] for details

−u · ∇ u = ∇ ·

   2  νT ∇ u¯ + (∇ u¯ )T − kI , 3

(37)

where k is the turbulent kinetic energy and I is the identity matrix and (in 3D case)

⎡ ∂ u¯ 1

∂ u¯ 1 ∂ x1 + ∂ x1 ∂ u¯ 1 T ⎣ ∇ u¯ + (∇ u¯ ) = ∂ x2 + ∂∂ ux¯12 ∂ u¯ 3 ∂ u¯ 1 ∂ x1 + ∂ x3

∂ u¯ 3 ⎤ ∂ u¯ 1 ∂ x3 + ∂ x1 ∂ u¯ 3 ⎦ ∂ u¯ 2 ∂ x3 + ∂ x2 . ∂ u¯ 3 ∂ u¯ 3 ∂ x3 + ∂ x3

∂ u¯ 1 ∂ u¯ 2 ∂ x2 + ∂ x1 ∂ u¯ 2 ∂ u¯ 2 ∂ x2 + ∂ x2 ∂ u¯ 3 ∂ u¯ 2 ∂ x2 + ∂ x3

νT =

ω

,

ν u¯ +∇ ·(νT ∇ u¯ )

(39)

Low Reynolds number model (LRN) One of the common problems in turbulence modelling is computing turbulent flow influenced by the adjacent wall. There is a boundary layer, where the velocity changes from the no-slip condition at the wall to its free stream value. The standard method of solving this problem is to apply a very fine mesh close to the wall. This is the so-called integration method, which necessitates an LRN type of turbulence model. On the other hand, this approach does not require any wall-function approximation which is common tool for obtaining near wall values in other methods. The region under wall influence diminishes in the case of higher Reynolds numbers. Low Reynolds number version of Wilcox (2006) k − ω two-equation model This model (which we use in our solver) is described in [21] in more detail. The equations for k and ω are defined as

∂k + u¯ · ∇ k = P − β ∗ kω + ∇ ·[(ν + σk νt )∇ k], ∂t ∂ω γω σ + u¯ · ∇ω = P − βω2 + ∇ ·[(ν + σω νt )∇ω] + d ∇ k · ∇ω, ∂t k ω (40)

ω

,

⎛

⎜ β = 0.09⎝ ∗

 σd =

ReT =

1 , 8

νω

β∗ α∗

ReT Rβ

ReT Rβ

where

aIP (u¯ , v ) = (ν∇ u¯ , ∇ v ) +

j

,

α∗ =

(41)

α0∗ + 1+

ReT Rk ReT Rk

,

α0∗ =

β0 3

,

(43)

γ h

(44)

(n  u¯ , n  v ) −

(45)

d (u¯ , v ) = −(νT {∇ u¯ + (∇ u¯ ) }, n  v ) + T

+ ( νT



 ∇ u¯ + (∇ u¯ )T , ∇ v ) − (n  u¯ , νT {∇ v} )

and

lIP (v ) =

2 2 (k, ∇ · v ) − ({k}, n · v ) . 3 3

(46)

Further, implicit Euler method is again used for time discretization and Picard’s method is used for linearization of the problem. Then the weak formulation of system (39) is to find u¯ ∈ V, p¯ ∈ L2 () such that ∀v ∈ V and ∀q ∈ W

1





n+1 n+1 , v + aIP (u¯ , v) + u¯



j

7 . 8

−(ν{∇ u¯ }, n  v ) − (n  u¯ , ν{∇ v} ) ,

n+1

, v ) − b(v, pn+1 ) +

1  n  n+1 + { pn+1 }, n · v − d (u¯ , v) = u¯ , v + lIP (v ) t   n+1 n+1 b(u¯ , q ) − {q}, n · u¯  = 0.

4 ⎞

4

Clim =

¯ , v ) = ((u¯ · ∇ )w ¯ , v) + c (u¯ ; w 1  1 + {u¯ · n}w¯  − |{u¯ · n}|w¯ , v , 2 2 b(v, p) = ( p, ∇ · v ),

n

 Re  13 α0 + RωT ⎟ (α ∗ )−1 , ⎠, γ = 25 T 1 + Re Rω

,

 ∂ u¯ , v + aIP (u¯ , v ) + c (u¯ ; u¯ , v ) − ∂t −b(v, p) + ({ p}, n · v ) + d (u¯ , v ) = lIP (v ), ∀v ∈ V, b(u¯ , q ) − ({q}, n · u¯  ) = 0, ∀q ∈ W,

+c (u¯ ; u¯

for ∂∂xk ∂∂ω ≤ 0, xj j for ∂∂xk ∂∂ω > 0, x

0,

k

1+





2Si j Si j



 1  P = νt ∇ u¯ + (∇ u¯ )T , 2 100 β270 +

(42)

3.3.1. Discontinuous Galerkin method for RANS This section is devoted to the application of discontinuous Galerkin method on the RANS equations. The goal is to use the model based on RANS equations on multipatch domains with discontinuous Galerkin technique on interfaces between the patches. If IPM is similarly applied as in the Section 3.2.1 on the system (39), the weak formulation has the form

t

where

νt = α

%

ωˆ = max ω, Clim

2 ∂ u¯ T + u¯ ·∇ u¯ = −∇ p¯ + ∇ ·[(ν + νT )∇ u¯ ] +∇ · (νT (∇ u¯ ) ) − ∇ k, ∂t 3   

k



(38)

∇ · u¯ = 0.

1 9

α0 = , β0 = 0.0708.

Rω = 2.61,

As mentioned above, solving the two equations model leads to the approximation of the Reynolds stresses (37). However, using LRN Wilcox k − ω model (2006), the turbulent viscosity ν T is not defined by the expression (38), but the modified expression is used instead

where

where ω is the specific dissipation. Our goal is to determine these quantities k and ω in order to approximate (37). Several approaches can be applied but the common way is k − ω model which is described in the following paragraphs. Substituting (37) into (35), RANS equations can be written in the form

∗

Rk = 6,

k ω νT = α ∗ , (thus νT = νt ), ωˆ ωˆ

Function ν T is the so-called turbulent viscosity defined by

k

Rβ = 8,



(47)



RANS equations are coupled with the k − ω model (40). Thus, it is necessary to apply discontinuous Galerkin method also on (40), similarly to RANS equations in this section. Therefore, the following weak formulation is obtained



 ∂k , v + ck (u¯ , k, v ) + ak (k, v ) + (β ∗ kω, v ) = (P, v ), ∂t   ∂ω , v + cω (u¯ , k, v ) + aω (k, v ) + ∂t

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

ARTICLE IN PRESS

JID: ADES

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

+(βω2 , v ) − dω (ω, k, v ) =

γ ω k



P, v ,

(48)

∀v ∈ W, where the bilinear and trilinear forms are defined by the following equalities

•

•

•

1 ({u¯ · n}k − |{u¯ · n}|k, v ), 2 1 cω (u¯ , ω, v ) = (u¯ · ∇ω, v ) + ({u¯ · n}ω − |{u¯ · n}|ω, v ), 2 ck (u¯ , k, v ) = (u¯ · ∇ k, v ) +

ak (k, v ) = ( (ν + σk νt )∇ k, ∇ v ) +

γ

h

( n  k, n  v )

−((ν + σk νt ){∇ k}, n  v ) − −(n  k, (ν + σk νt ){∇v} ) ,

aω (ω, v ) = ( (ν + σω νt )∇ ω, ∇ v ) +

τ

h

(49)

7

evaluation of basis functions and their derivatives at Gauss points on a given element, evaluation of inverse Jacobi matrix and determinant of Jacobi matrix at Gauss points on a given element, evaluation of integrals using Gauss quadrature on a given element to form local matrix.

From the computation time point of view, the fastest way to implement this type of computation in Wolfram Mathematica is to exploit automatic conversion of functions into C code, compilation of such a code with C compiler and linking the output dll file back to Wolfram Mathematica. Then, whenever this function is called, the computation is done by this compiled function. To be more specific, any function defined as

( n  ω , n  v )

−((ν + σω νt ){∇ω}, n  v )

−(n  ω, (ν + σω νt ){∇v} ) , dω (ω, k, v ) =

σ

 ∇ k · ∇ω, v ω & & σ ' '  1  σ  d + ∇ k ·n ω −  d ∇ k ·n ω, v . 2 ω ω d

The fully discretized problem (using implicit Euler method and Picard’s method) has the weak formulation

1  n+1  n k , v + ck (u¯ , kn+1 , v ) + t + ak (kn+1 , v ) + (β ∗ kn+1 ωn , v ) =

1

t

( kn , v ) + ( P n , v ),

1  n+1  (50) ω , v + cω (u¯ n , ωn+1 , v ) + aω (ωn+1 , v ) + t  n  1 ω + (βωn ωn+1 , v ) − dω (ωn+1 , kn , v ) = (ω n , v ) + γ n P n , v . t k

is converted to C code as mentioned above. For good performance of such compiled function, it is necessary to write as simple code as possible and to avoid using high-level Wolfram Mathematica functions which will require calling Wolfram Mathematica Kernel from generated C code. The computation needed for assembling of matrices mentioned above are realized with the help of this C conversion and compilation. For example, function which returns values of integrals in matrix A for all pairs of non-zero basis functions on a given element (see (13)) is

4. Implementation This section is devoted to some implementation details related to our IgA-based flow solver. Our IgA-based solver is implemented in Wolfram Mathematica. One of the reasons for this choice is faster development speed of a code in this software package than in low level programming languages. Nonetheless, we put an effort to exploit all available features provided by Wolfram Mathematica to make this solver as fast as possible. If a geometric definition of a computational domain is neglected, the crucial parts of the solver are: •

•

assembling matrices arising from different terms in weak formulation of the problem and setting up the resulting linear system of equations, solving linear system of equations.

On the contrary to the standard FEM, time spent on assembling matrices necessary for setting the final system of linear equations is not negligible with respect to time spent on solving this linear system, it can easily be comparable with respect to the computation time. This is a well-known problem of IgA and there is a huge effort to reduce time-complexity of this step of isogeometric analysis in the research community (see [22,23]). The reasons for that are obvious: basis functions are polynomials of a higher degree, or even rational functions which are more difficult to evaluate, and because of higher degree and/or rational basis functions, Gaussian quadrature with more Gauss points has to be used to obtain accurate enough result. Focusing on assembling matrices, the main parts which need to be computed are:

where gradsingp contains gradients of basis functions on the element, detingp determinant of Jacobi matrix evaluated at Gauss points and gweights are Gauss weights following from Gaussian quadrature. Concerning the second main part of the solver, solving linear systems of equations, standard function LinearSolve[] is used, which performs well in many cases, even in comparison with libraries like PETSc (see [24]). It allows to select from different solvers: • •

direct solvers, iterative solvers (preconditioned conjugate gradient, preconditioned stabilized biconjugate gradients, preconditioned GMRES).

As the structure of linear systems arising from IgA applied on flow problems on multipatch domains is very complicated (see e.g. Fig. 4), all tested preconditioning strategies recommended in the literature fail to provide accurate enough solutions of these linear systems. This is the reason why currently direct solver is used in our solver. We further investigate possible preconditioning strategies for this case because direct solvers are too memoryconsuming when large problems are solved. The rest of this section is devoted to some implementation details related to discontinuous Galerkin method used for coupling NURBS patches along common interfaces. Let us recall the weak formulation presented in (30). Terms which do not contain “jump” or “mean” operators follow from the

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

JID: ADES 8

ARTICLE IN PRESS

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

coupling matrix is

(

)

Mr,r Ms,r

Mr,s , Ms,s

where r, s are indices of NURBS patches defining a common interface. These indices also determine where this matrix will be put in the resulting block matrix. Now, a closer look is taken to all new terms following from application of discontinuous Galerkin method along common interfaces between NURBS patches. Term (n  u, n  v ) can be expressed in matrix form as follows

(

K Fig. 3. Structure of matrix A for computational domain composed of 4 NURBS patches without adding coupling matrices.

r,s

Rui,r Ruj,r = − r,s Rui,s Ruj,r

)  − r,s Rui,r Ruj,s  , u u r,s Ri,s R j,s

(51)

where {r, s} are indices of NURBS patches representing a common interface, r, s is the corresponding part of for this common interface, w stands for x, y or z, nr, x is then x-coordinate of outer unit normal of NURBS patch no. r, Rui,r is i-th NURBS basis function for velocity of r-th NURBS patch, i, j runs through global indices of nonzero NURBS basis functions on r, s (for index i, NURBS basis functions corresponding to control points with given Dirichlet boundary condition are excluded). For each common interface one matrix Kr,s w is obtained and the resulting matrix is then just a sum of all of these matrices

K= Fig. 4. Left: Structure of top-left block of matrix (57) after including all couples for computational domain composed of 4 NURBS patches; Right: Structure of (57) after including all couples for computational domain composed of 4 NURBS patches.

r,s



Kr,s .

(52)

{r,s}

Similarly, term ({ ph }, n · v ) can be included into the matrix formulation of the problem via matrices Lw obtained as a sum of matrices

(1 

Lr,s w =

2

nr,w Rui,r R pj,r u p r,s nr,w Ri,s R j,r

r,s

− 12

)



1 2

nr,w Rui,r R pj,s . u p r,s nr,w Ri,s R j,s

r,s

− 12

(53)

Transposed matrices then correspond to the term ({q}, n · uh  ) . The last term following from steady Stokes problem is (ν{∇ u}, n  v ) which can be analogously described by matrices



Fig. 5. Backward facing step – the computational domain is composed of four NURBS patches. Patches are distinguished by colors.

standard Stokes formulation on a single-patch domain, i.e., the corresponding matrix formulation is shown in (11) and represented by matrices A, B1 , B2 and B3 . First of all, description of the structure of matrices A, B1 , B2 , B3 without including any couples between NURBS patches is needed. Control points (and the corresponding basis functions) are indexed incrementally through all patches, i.e., if the last control point of NURBS patch no. 1 has an index n, then the first control point of NURBS patch no. 2 has index n + 1. For instance, the structure of matrix A for multipatch domain composed of four NURBS patches is shown in Fig. 3. One can see that the matrix is decomposed into four block matrices where each block matrix corresponds to one NURBS patch. Since indices of control points (basis functions) are global, as mentioned above, such a structure of the whole matrix directly follows. As each couple is defined between two NURBS patches, this block structure of matrices is used during assembling of coupling matrices. For computational domain composed of four NURBS patches, each matrix can be taken as composed of 4 × 4 block matrices (see Fig. 3). Blocks determined by volumetric integrals always correspond only to one NURBS patch and always lie on the diagonal in this block matrix. Further, typical structure of a

Mr,s u,v

=

∂ Ru

nr,u Rui,r ∂vj,r  ∂ Ru − r,s nr,u Rui,s ∂vj,r r,s



∂ Ru

nr,u Rui,r ∂vj,s  ∂ Ru , − r,s nr,u Rui,s ∂vj,s r,s

(54)

where u, v again stands for x, y or z. Term (n  u, ν{∇ v} ) again corresponds to transposed matrices Mr,s u,v . For the steady Navier-Stokes problem, the coupling term ( 12 {u · n}u − 12 |{u · n}|u, v ) following from convective term of this problem can be expressed as

(

P(uk )r,s =

r,s {u

( Q(uk )r,s =

)  − r,s {uk · n}Rui,r Ruj,s  , k u u r,s {u · n}Ri,s R j,s

· n}Rui,r Ruj,r − r,s {uk · n}Rui,s Ruj,r k

|{uk · n}|Rui,r Ruj,r − r,s |{uk · n}|Rui,s Ruj,r r,s

(55)

)  − r,s |{uk · n}|Rui,r Ruj,s  , (56) k u u r,s |{u · n}|Ri,s R j,s

where uk is the velocity from previous iteration step of Picard iteration. The final matrix of a linear system for steady Stokes problem on multipatch domain is then obtained by adding coupling matrices to the original block matrix, i.e.,

⎡

A + γh K + M ⎢ 0 ⎣ 0 B1 − L x

0 A + γh K + M 0 B2 − L y

0 0 A + γh K + M B3 − L z

⎤

−BT1 + LTx −BT2 + LTy ⎥ ⎦, (57) −BT3 + LTz 0

where

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

JID: ADES

ARTICLE IN PRESS

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

Fig. 6. Velocity of the flow in the whole domain. Upper figure - IgA based solver, lower figure - OpenFOAM solver.

Fig. 7. Detail of the velocity field. Upper figure - IgA based solver, lower figure - OpenFOAM solver.

Fig. 8. Velocity vectors. Upper figure - IgA based solver, lower figure - OpenFOAM solver.

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

9

ARTICLE IN PRESS

JID: ADES 10

[m5G;July 5, 2017;14:58]

B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

Fig. 9. Profiles of x-axis component of the velocity. X coordinates are describes in the brackets. Left figure - IgA based solver, right figure - OpenFOAM solver.

Fig. 10. Pressure contours. Upper figure - IgA based solver, lower figure - OpenFOAM solver.

M

=





−ν Mxx + Myy + Mzz + MTxx + MTyy + MTzz ,

(58)

The structure of the top-left block of this matrix and of the whole matrix is shown in Fig. 4. For steady Navier-Stokes problem, N(uk ) + P(uk ) − Q(uk ) is added to the main block diagonal in each step of Picard iteration.

Remark 2. Analogous approach can be used for the case of RANS equations with added terms with turbulent viscosity and turbulent kinetic energy (cf. (44)) and also for the equations representing k − ω model (cf. (48)).

5. Results Results obtained with using our solver for a standard test example – backward facing step in 3D are presented in this section. The computational domain is shown in Fig. 5. Inflow of the domain is set on the front surface of blue patch, where paraboloid-like input velocity profile with maximum umax = (umax , 0, 0 ) is prescribed. Outflow of the domain is set on back surface of yellow/green patches, i.e., homogeneous Neumann boundary condition is prescribed on these surfaces. Other boundary surfaces of the domain are considered to be solid walls, i.e., Dirichlet boundary condition with u = 0 is prescribed. For pressure, homo-

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

JID: ADES

ARTICLE IN PRESS B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

[m5G;July 5, 2017;14:58] 11

Fig. 11. X-axis components of velocities near the bottom. Based on these results, reattachment lengths are 29.63 (IgA solver) and 27.95 (OpenFOAM solver). The value, visually identified from the graph in [25], is 30.

geneous Neumann boundary condition is prescribed on the whole boundary. The following examples differ in the proportions of the computational domain and values of velocity umax and kinematic viscosity ν . Turbulence model is used only if is stated. 5.1. Comparative 3D steady flow The proportions of the computational domain and flow parameters are chosen in the agreement with the experiment published in [25]. The step is located at x = 0 and the height is 4.9. Inlet boundary is located at x = −17.15 and its height is 5.2. Outlet boundary is located at x = 220. The width of the domain is 90. The input velocity parabolic profile has the maximum umax = (30, 0, 0 ). Kinematic viscosity is set ν = 0.7509 such that the results can be visually compared with the presented ones in [25] for

Re =

2 umax · D = 277 3 ν

(59)

with hydraulic diameter D = 10.4. Free open source software OpenFOAM was used for the detail comparison. Computational mesh was refined in the step neighborhood and contains 650 0 0 0 volumes. The selected numerical scheme is based on bounded Gauss upwind divergence scheme, which is in detail described in [26]. Comparison of the results computed by our IgA solver (3672 elements mesh) and OpenFOAM solver are made in the middle plane of the domain and depicted at the Figs. 6, 7, 8, 9, 10, 11. It can be seen that velocity fields are very similar in the whole computational domain including eddy area behind the step. This is also in the concordance with the velocity profiles presented in [25]. Pressure contours are also very similar with small difference in the range of values. 5.2. 3D flow with RANS model For this case, the inlet and step heights and their distance are equal and set to 1. The maximum inlet velocity umax = 3 and the

Fig. 12. Solution of RANS equations with k − ω model after 1, 2 and 3 time units.

kinematic viscosity ν = 0.0 0 04 are used. This corresponds to the flow with Reynolds number

Re =

2 umax · D = 10 0 0 0 3 ν

(60)

with hydraulic diameter D = 2. Altogether, mesh contains 3360 NURBS elements. Timestep is chosen as t = 0.01. Let us mention also the setting of the turbulence model, especially of k − ω model. Boundary conditions for RANS equations (i.e. for time-averaged velocity and pressure) are set as mentioned above and of k − ω are set as follows: •

turbulent kinetic energy k – on inflow boundary k is set with the help of formula [21]

k=

3 (u¯ I )2 , 2

(61)

where u¯ is flow velocity and I is turbulence intensity (I = 0.01 here). Further, k = 0 on solid walls and homogeneous Neumann boundary condition is used on outflow boundary.

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012

JID: ADES 12 •

ARTICLE IN PRESS B. Bastl et al. / Advances in Engineering Software 000 (2017) 1–12

specific dissipation ω – on inflow boundary ω is set with the help of formula [21]

ω

√ k = , 0.07D

(62)

where D is inflow width (i.e., D = 1 here). Further, ω = 10 0 0 on solid walls and homogeneous Neumann boundary condition is used on outflow boundary. Initial conditions for all quantities had to be specified, i.e., • •

[m5G;July 5, 2017;14:58]

velocity, pressure – solution of the steady Stokes problem, turbulence kinetic energy k, specific dissipation ω – solution of the steady and linearized k − ω model presented in (48) for initial velocity field,

Fig. 12 shows streamlines of velocity obtained in the time evolution for 1, 2 and 3 time units. The growing eddy can be seen behind the step, which disintegrates with increasing in time, and also a new eddy near the upper wall. 6. Conclusion and future work The paper was devoted to the numerical simulation of turbulence flows based on RANS equations with k − ω model with the help of isogeometric analysis. This recently proposed method allows to eliminate discretization error even for very complex computational domainss. For example, water turbines for which the domain is exactly described by a composition of NURBS objects. The weak formulation of RANS equations with LRN k-omega was presented in the context of isogeometric analysis on, generally, non-conformal multipatch NURBS domains. The results of the implemented solver were presented on several examples and compared with previously published results and results obtained with standard CFD package. In the future, we want to focus on testing and comparison of other turbulence models (e.g. SST model, Spalart-Allmaras) and their applicability for solving different problems. Acknowledgments The authors of this paper were supported by Technology agency of the Czech Republic through the project TA03011157. References [1] Saad Y. Iterative methods for sparse linear systems. Philadelphia: SIAM; 2003. [2] Versteeg H, Malalasekera W. An introduction to computational fluid dynamics - second edition. Prentice Hall; 2007.

[3] Ekaterinaris JA. High-order accurate, low numerical diffusion methods for aerodynamics. Progress in Aersopace Sciences 2005;41:192–300. [4] Elman HC, Silvester D, Wathen AJ. Finite elements and fast iterative solvers with applications in incompressible fluid dynamics. Oxford University Press; 2014. [5] Gerbeau J, Farhat C. The finite element method for fluid mechanics. Stanford University; 2009. [6] Segal A. Finite element methods for the incompressible Navier–Stokes equations. Delft University of Technology: Burgerscentrum, Research School for Fluid Mechanics; 2012. [7] Augustin M, Caiazzo A, Fiebach A, Fuhrmann J, John V, Linke A, et al. An assessment of discretizations for convection–dominated convection–diffusion equations. Comput Methods Appl Mech Eng 2011;200:3395–409. [8] Hughes T, Cottrell J, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194(3941):4135–95. [9] Cottrell J, Hughes T, Bazilevs Y. Isogeometric analysis. John Wiley & Sons, Ltd; 2009. [10] Moin P, Mahesh K. Direct Numerical Simulation: A tool in turbulence research. Annu Rev Fluid Mech 1998;30:539–78. [11] Berselli L, Iliescu T, Layton W. Mathematics of large eddy simulation of turbulent flows. Springer; 2006. [12] Davidson L. Fluid mechanics, turbulent flow and turbulence modeling. Chalmers University of Technology: Lecture Notes in MSc courses; 2013. [13] Furbo E, Harju J, Nilsson H. Evaluation of turbulence models for prediction of flow separation at a smooth surface. Sweden: Tech. rep., Uppsala Universitet; 2009. [14] Furbo E. Evaluation of RANS turbulence models for ow problems with significant impact of boundary layers. Sweden: Uppsala Universitet; 2010. Master’s thesis [15] Nichols RH. Turbulence models and their application to complex flows. Tech. rep., University of Alabama at Birmingham; 2010. [16] Kuzmin D, Mierka O. On the implementation of the k −  turbulence model in incompressible flow solvers based on a finite element discretization. Germany: Tech. rep., University of Dortmund; 2006. [17] Petit O. Towards full predictions of the unsteady incompressible flow in rotating machines using openFOAM. Sweden: Chalmers University of Technology; 2012. Ph.D. thesis. [18] Nguyen VP, Kerfriden P, Brino M, Bordas SPA, Bonisoli E. Nitsches method for two and three dimensional NURBS patch coupling. Comput Mech 2014;53(6):1163–82. [19] Langer U, Mantzaflaris A, Moore SE, Toulopoulos I. Multipatch discontinuous Galerkin isogeometric analysis. Tech rep, NFN Technical Report 2014;18. [20] Montlaur A. High-order discontinuous Galerkin methods for incompressible flows. Spain: Escola Politécnica Superior de Castelldefels, Barcelona; 2009. Ph.D. thesis. [21] Wilcox D. Formulation of the k-omega turbulence model revisited. AIAA J 2008;46:2823–38. [22] Hughes T, Reali A, Sangalli G. Efficient quadrature for NURBS-based isogeometric analysis. Comput Methods Appl Mech Engrg 2010;199:301–13. [23] Mantzaflaris A, Jüttler B. Integration by interpolation and look-up for Galerkin-based isogeometric analysis. Comput Methods Appl Mech Eng 2015;284:373–400. [24] PETSc (Portable, Extensible Toolkit for Scientific computation) library(2015). URL http://www.mcs.anl.gov/petsc/. [25] Jiang BN, Hou LJ, Lin TL. Least-squares finite element solutions for three-dimensional backward-facing step flow. Tech rep, NASA Technical Memorandum 1993;106353. [26] Openfoam user guide. URL http://cfd.direct/openfoam/user-guide/.

Please cite this article as: B. Bastl et al., IgA-Based Solver for turbulence modelling on multipatch geometries, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.06.012