hybrid Dirichlet formulation for wall-bounded flow problems including turbulent flow

Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A mixed/hybrid Dirichlet formulation for wall-bounded flow problems including turbulent flow Peter Gamnitzer a,1, Volker Gravemeier a,b,⇑, Wolfgang A. Wall a a b

Institute for Computational Mechanics, Technische Universität München, Boltzmannstr. 15, 85747 Garching, Germany Emmy Noether Research Group ‘‘Computational Multiscale Methods for Turbulent Combustion’’, Technische Universität München, Boltzmannstr. 15, 85747 Garching, Germany

a r t i c l e

i n f o

Article history: Received 12 September 2011 Received in revised form 14 June 2012 Accepted 16 June 2012 Available online 7 July 2012 Keywords: Embedded Dirichlet formulation Weak boundary conditions Law of the wall Variational multiscale method Turbulent channel flow

a b s t r a c t In this study, a new method for the weak imposition of Dirichlet boundary conditions in convectiondominated convection–diffusion and wall-bounded turbulent flow problems is developed. It is based on an embedded Dirichlet formulation which utilises an additional discontinuous stress field for imposing boundary conditions. Values at outflow boundaries may not be sufficiently controllable on underresolved meshes, e.g., when using Nitsche’s method for weak enforcement of Dirichlet boundary conditions. Hence, a consistent modification of an earlier-proposed constant-parameter embedded Dirichlet formulation that allows to regain control if required is introduced. A separate treatment of boundary conditions in wall-normal and tangential direction is enabled. This split is physically motivated, and it turns out to be crucial for a successful application of the method. On the one hand, it allows for including a formulation based on the law of the wall in tangential direction. On the other hand, a consistent modification derived from theoretical considerations for convection–diffusion problems is applied in wall-normal direction, preventing unphysical penetration of walls. Numerical results obtained for turbulent channel flow at modest friction Reynolds number Res ¼ 395 employing meshes with and without adequate resolution of the boundary layer are presented. Superior performance of the proposed method compared to strongly-imposed boundary conditions is demonstrated by an error reduction of up to 50% for mean velocity results in the core of the channel on coarse discretisations. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Many applications of engineering interest involve turbulent flow at very high Reynolds numbers. A promising approach for the simulation of such flows is large eddy simulation (LES), where the larger, energetic scales are resolved while the more universal, smaller scales are accounted for by a model. Boundary layers in such flows are especially challenging for the LES approach. In the region next to the wall, the dynamically important motions become progressively smaller with increasing Reynolds number. Consequently, the resolution requirements for LES in those regions can become prohibitively high. A well established approach to overcome this problem is to model the so-called inner layer at the wall. As a result, if the viscous wall region is modelled and only the outer layer needs to be resolved, substantial computational savings are enabled. The

⇑ Corresponding author. Tel.: +49 (0)89 289 15245; fax: +49 (0)89 289 15301. E-mail addresses: [email protected] (P. Gamnitzer), vgravem@lnm. mw.tum.de (V. Gravemeier), [email protected] (W.A. Wall). 1 Present address: Institute of Basic Sciences in Civil Engineering, Unit of Strength of Materials and Structural Analysis, Universität Innsbruck, Technikerstr. 13, 6020 Innsbruck, Austria. 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.06.013

simplest modelling approach for the inner layer is a wall-stress model such as the early ones proposed by Deardorff [1] and Schumann [2], see, e.g., the reviews in [3–5]. In these models, the effect of the inner layer is replaced by a tangential stress acting on the wall. This stress can be related to the outer velocity, for instance, by a wall law. In wall-normal direction, an impermeability condition is usually assumed. Note that a significant flux across the boundary would lead to a potentially unphysical transport of momentum and thus cause incorrect predictions of the drag as well as instabilities. For a recent overview on other approaches to wall-layer modelling such as zonal approaches and hybrid RANS/LES methods, the reader is referred to the review in [5], and for more on the latter, hybrid RANS/LES methods, to another review in [6] published in the same year. The idea to incorporate wall laws in the finite element analysis of turbulent flow was recently brought up in [7]. Therein, a weak imposition of the Dirichlet boundary condition via Nitsche’s method [8] was proposed for convection–diffusion and fluid flow problems. In a particular version of the approach, an adjustment of the underlying Nitsche parameter according to a formulation of the law of the wall due to Spalding [9] was demonstrated in [10,11]. In those publications, it was shown that mean velocity profiles can be significantly improved by a weak imposition of

23

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

the boundary condition at the wall. In this study a new strategy for weakly imposing Dirichlet boundary conditions in flows at high Reynolds number is proposed. It is (i) based on the embedded Dirichlet method as proposed in [12], (ii) includes a wall-stress model to impose a tangential traction boundary condition, and (iii) guarantees the correct impermeability of the wall by a consistent modification of the original approach. The impermeability of the wall is ensured even though our approach enforces the Dirichlet boundary condition weakly in both directions, tangential and normal. This treatment allows for straightforward handling of complex, curved geometries, since local coordinate systems are not required to enforce a strong Dirichlet boundary condition in the wall-normal direction. The method we develop in this paper might be used, for instance, in conjunction with a quasi-static residualbased finite element method (see, e.g. [13–15]), a time-dependent residual-based finite element method (see, e.g. [16,17]), or an algebraic variational multiscale multigrid method (AVM3) (see, e.g. [18–20]). Those methods have been shown to be well suited for LES of turbulent flows. However, our approach can also be used in a rather standard stabilised finite element method, as it was proposed in, e.g., [21–23] (see also [24] for an overview). The imposition of the boundary condition presented here does not interfere with the extra terms added for the variational multiscale modelling. The remainder of this paper is organised as follows. In Section 2, the weak imposition of boundary conditions for the convection–diffusion equation is presented. A modified parameter-based formulation is proposed which better ‘controls’ boundary values on coarse meshes. The choice of this parameter is based on an analysis of a simplified model problem. The implementation of the method is briefly described, and results obtained from the modified approach are compared to results obtained from a constant-parameter formulation and Nitsche’s method. In Section 3, a framework for applying boundary conditions weakly for wallbounded turbulent flows using our method is presented. A separate consideration of boundary conditions in tangential and normal direction based on the modified approach shown in Section 2 in wall-normal direction is proposed. In tangential direction, the wall-law-based approach from [10] is adopted to the present setting. This combined approach and its implementation is described in detail, and results from turbulent channel flow computations, which allow for a quantitative comparison to other methods, are presented. In Section 4, we draw conclusions and give a short outlook to future research directions.

ða  $u; v ÞX þ jð$u; $v ÞX  ðf ; v ÞX  ðq  n; v Þ@ X ¼ 0:

ð4Þ

This statement is formulated using the usual L2 -inner products ð; ÞX on X and on @ X, respectively. The outward-pointing unit normal vector on the domain boundary is denoted by n. The flux term on the boundary arises from partial integration and will subsequently be used to weakly enforce the Dirichlet boundary condition. To obtain a finite element formulation, the domain is partitioned into non-overlapping elements Xe with a characteristic length h. Based on these element domains and associated basis functions, finite-dimensional subspaces of the infinite-dimensional spaces of test and solution functions are generated. In our approach, an approximate solution uh for u is computed by solving the following stabilised discrete problem:



       a  $uh ; v h X þ j $uh ; $v h X  f ; v h X  qh  n; v h @ X X   1 þ a  $uh  jDuh  f ; sa  $v h Xe þ qh ; r h X 

e

h

 $u ; r

h

 X



 

sB

  g  u n; rh @ X ¼ 0: h

j

ð5Þ

Problem (5) is a linear equation for the solution of the unknown discrete problem variable uh and the unknown discrete fluxes qh . Accordingly, the method has to be considered as a mixed method. The test function corresponding to the problem variable is denoted by v h , and the test function for the fluxes in equation (5) is denoted by r h . The first line in Eq. (5) is the Galerkin part including the consistency term arising from partial integration and Eq. (3). The second line corresponds to the residual-based SUPG stabilisation term proposed in [21]. Furthermore, following the approach of [12], the first two terms in the third line represent a weak formulation of the balance of fluxes according to (3). The last term in the third line enforces the boundary condition on @ X weakly including an additional parameter sB . For one-dimensional examples, the following stabilisation parameter definition is used:

s¼

  h kakh min 1; 2kak 6j

ð6Þ

Although not used here, various other definitions and generalisations of this stabilisation parameter are available (see, e.g., [25, 26]). For multi-dimensional problems, the element length h can be replaced, for instance, by a stream length as in [23], and kak can be interpreted as the Euclidian norm of a. In [12], it was proposed to use a constant parameter for the enforcement of the boundary condition, that is,

2. Convection–diffusion

sB ¼ 1:0:

2.1. Strong and weak form of convection–diffusion equation

However, we propose a modification of the parameter as follows:

The strong form of the (stationary) convection–diffusion problem on domain X subjected to Dirichlet boundary conditions on the domain boundary @ X is given by

a  $u  jDu  f ¼ 0 in X;

ð1Þ

u  g ¼ 0 on @ X:

ð2Þ

In this equation, u is the transported quantity, a is the convective velocity vector, which is assumed to be solenoidal, j the positive diffusion coefficient, f a source term, and g the prescribed boundary value. The diffusive flux q vector associated with this problem is given by

q  j$u ¼ 0 in X:

ð3Þ

Using appropriate function spaces, Eq. (1) can be restated in its weak form:

ð7Þ

sB ¼ sB ðPee Þ ¼ 1 þ 2Pee ;

ð8Þ

with the element Peclet number

Pee ¼

kakh

j

:

ð9Þ

Unlike for the SUPG stabilisation parameter, a generalisation of definition (8) to multidimensions will not be based on a definition of h as a stream length but rather as the size of the element normal to the boundary, as it was done in [10]. The rationale on which this modification is based will be addressed in Section 2.2. This modification allows to keep control of the boundary condition even for convection-dominated flows on coarse meshes. Eq. (5) can be interpreted as a consistent penalty formulation. The penalty parameter is either constant or bounded with respect to mesh refinement. In the limit h ! 0, parameter sB ðPee Þ convergences to the constant value of the original method from Eq. (7). Hence, we can expect that the convergence properties of the original method are not altered

24

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

for the modified approach when considering meshes which are sufficiently refined towards the boundary. We use the same basis functions (e.g., linear basis functions) both for the primary problem variable and fluxes in the following. For uh and v h , continuous basis functions are chosen. For the fluxes the basis functions are allowed to be discontinuous across element boundaries, allowing for local condensation of the fluxes on the element level. Since equation (5) depends only on flux values at the boundaries, this condensation is required only for elements adjacent to the boundary. The boundary elements are assumed to constitute a subdomain which is denoted by XH , such that Eq. (5) can be reformulated as



       a  $uh ; v h X þ j $uh ; $v h X  f ; v h X  qh  n; v h @ X X   1 þ a  $uh  jDuh  f ; sa  $v h Xe þ qh ; r h XH

Fig. 1. Domains and function definitions used in a one-dimensional model problem for linear basis functions.

j

e

       $uh ; r h XH  sB g  uh n; r h @ X ¼ 0:

ð10Þ

To sum up, the extra fluxes are only considered on a minor part of the domain next to the boundary, while in the interior of the domain a standard stabilised finite element method is used. All additional flux degrees of freedom are eliminated on the element level by static condensation. Thus, the term ‘mixed/hybrid’ appears appropriate for this method, as proposed in [12]. 2.2. Analysis for a simplified model problem In this subsection, we would like to provide evidence for the choice of the parameter as given in Eq. (8). For this purpose, a sufficient condition for the uniqueness of the solution for a simplified problem without SUPG stabilisation is given. We show how the parameter choice according to (8), even for coarse discretisations, allows to include the value at the boundary in the statement such that uniqueness of the solution is ensured. The present analysis is motivated by considerations in [27] on a symmetric method for weakly imposing boundary conditions in the context of embedded meshes. A bilinear form associated with equation (10) reads:

~ e  n > x for every unit normal n on @ Xe \ @ X. It poses the property n a restriction on the curvature of the boundary as well as a restriction on the minimum angle of kinks allowed on @ Xe \ @ X. ~e With this assumption, definition (14) could be made using n instead of n and the proof could be done in the same steps now yielding a definition of sB depending also on x. However, for the ease of clarity, we display the proof only for the initial simplifying assumption. In definitions (13) and (14), we only selected special members of the finite-dimensional function spaces which will turn out beneficial for the subsequent proof. The second summand in (14) will allow to prove that the bilinear form for the special test functions (13) and (14) is bounded from below by a positive expression including the norm of the unknown variable on the boundary. We obtain

  1 jb h u@ X n B uh ; qh ; uh ; qh þ h sB       1  h h q ; q XH ¼ a  $uh ; uh X þ j $uh ; $uh X  qh  n; uh @X þ

jsB

        B uh ; qh ; v h ; rh ¼ a  $uh ; v h X þ j $uh ; $v h X  qh  n; v h @X      1 þ qh ; r h XH  $uh ; rh XH þ sB uh n; r h @ X : ð11Þ



j

For simplicity, we assume that n is constant at each element boundary Xe \ @ X (i.e., straight edges and only one boundary edge per element are considered). Furthermore, it is assumed that all elements at the boundary are characterised by element length h. We introduce the function uh@ X as an element of the finite-dimensional function space for u which equals uh on the domain boundary,

  uh@ X @ X ¼ uh @ X ;

v h ¼ uh :

ð13Þ

Elementwise in XH , we define

1

sB

qh þ

jb h

$uh ; qh

h

$uh  n; uh@ X

uh@ X n;

ð14Þ

where b is a constant, positive parameter which will be specified later on. The above expression is well defined based on the uniqueness of the normal according to the simplifying assumption we made after equation (11). A more general result could be obtained based on the weaker assumption that there exists a global constant x > 12 such that for every element there exists a unit vector n~ e with

XH



XH

þ

sB jb  h

uh ; uh

 @X

:

ð15Þ

Taking into account that the first and last summand in the second line cancel out, using

Z  2 1 a  uh $uh dX ¼ a  $ uh dX 2 X X Z Z   2   2 1 1 h ¼ $ a u a  n uh dX dX ¼ 2 X 2 @X   1 P  kak uh ; uh @ X ; 2

  a  $uh ; uh X ¼

ð12Þ

but vanishes at all interior nodes of X. This function is sketched in Fig. 1. We choose the test functions such that

rh ¼





   b þ uh ; qh  n @X þ qh  n; uh@ X XH h

1

sB jb 

Z

ð16Þ

and the Cauchy–Schwarz inequality, we obtain

  1 jb h u@ X n B uh ; qh ; uh ; qh þ h sB 2 2 1 1 qh 2 H P  kak uh @X þ j $uh X þ X 2 jsB  þ

1 $uh

XH

sB sB jb h 2 h

u

h q H  b qh H uh H  jb $uh H uh H @X X @X X X X X h h

@X

:

Furthermore, Young’s inequality yields

ð17Þ

25

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

  1 jb h B uh ; qh ; uh ; qh þ u@ X n h sB h 2 2 1 1 qh 2 H P  kak  u @ X þ j $uh X þ X 2 jsB h 2 h 2 ! h 2 js  h 2 ! B q XH u@ X XH 1 j $u X H b q XH Ch    þ  jsB þ h 2 Ch sB 2 2j 2 h 2 h 2 ! jb h $u XH u@X XH sB jb h 2 þ u @X :  þ h h 2 2h

  B uh ; qh ; v h ; r h ¼ 0: Since (24) holds for all test functions ular that

 1 B uh ; qh ; uh ; qh þ

sB

ð18Þ

Here, C is a constant from the inequality

h 2 u H 6 Ch uh 2 ¼ Ch uh 2 ; @X X @X @X @X

ð19Þ

which provides an upper bound for the L2 -norm of the function uh@ X on the domain in terms of its L2 -norm on the boundary. Inequality (19) results from the shape of the support of uh@ X , which extends only for the size h of one element from the boundary to the interior of the domain. The functions uh@ X are zero on all interior nodes by definition. The nonzero part of uh@ X is defined exclusively by finite element basis functions associated with boundary nodes. Since these functions take their absolute maximum value on the boundary, we can expect C to be smaller than one. Applying equation (19) and using XH  X results in

  1 jb h u@ X n B uh ; qh ; uh ; qh þ h sB

2 2 2 2 j $uh X qh XH 1 1 qh 2 H   P  kak  uh @ X þ j $uh X þ X 2 jsB 2sB 2jsB Cb qh 2 H  bjsB uh 2  jb $uh 2  jbC uh 2  X X @X @X 2jsB 2h 2  2h  sB jb h 2 j jb h 2 1 Cb qh 2 H u @X ¼ j  þ  $u X þ  X h 2 2jsB 2jsB 2sB   h 2 sB jb jbC 1 þ ð20Þ   kak u @ X : 2h 2h 2

Assuming C < 1, we take b ¼ 12 and obtain

 1 B uh ; qh ; uh ; qh þ

sB

   3 j h 2 j $u X uh@ X n P 4 2h 2sB   h 2 1 C q H þ  X 2jsB 4jsB   2 j 1 þ ðsB  C Þ  kak uh @ X : 4h 2

j

ð21Þ The positive factors for all three norms hselected value for b ensures $u ; qh H , and uh if sB is sufficiently large. This is the case X X @X for the choice

sB :¼ 1 þ 2

kakh

j

¼ 1 þ 2Pee ;

ð22Þ

which implies

 1 B uh ; qh ; uh ; qh þ

sB

ð24Þ

 j 2 uh@ X n P $uh X 2h 4

j

  1 C qh 2 H 1 X 2ðj þ 2kakhÞ 2 j 2 þ ð1  C Þ uh @ X : 4h ð23Þ

j 2h

 uh@X n ¼ 0:



v h ; rh



, this implies in partic-

ð25Þ

From Eq. (23) we immediately obtain uh  0 on @ X and $uh  0 on X, hence uh  0 on X. Furthermore, qh  0 is implied on XH , and thus, the solution is unique. To highlight the consequences of parameter choice (8) once more, we again take a look at the structure of inequality (21). On the right-hand side of the inequality, we find three contributions. The first summand, which provides uniqueness of the solution gradients, as well as the second summand, which provides equality of the fluxes, are not sensitive with respect to the choice of sB . Any choice sB P 1, including the original approach, will provide positive coefficients and thus control over solution gradients and fluxes. The third term, which controls the magnitude of the solution, requires more attention. For the choice sB ¼ 1, the factor would be

ð1  C Þ

j 4h

1  kak: 2

ð26Þ

If h is sufficiently small, expression (26) will always be positive and thus the approach from [12] is sufficient to control the boundary value in all cases. Contrarily, if h is large and convection is strong, the expression may become negative, and our analysis is not sufficient to prove the uniqueness of the solution anymore. Nevertheless, in the numerical examples, we will observe that, in practice, singular problems do not occur, even on coarse meshes. The boundary value at the inflow will always be represented quite accurately, providing the required control on the magnitude of the solution. Only the value at the outflow boundary will not be controlled sufficiently on coarse meshes when using sB ¼ 1:0 instead of (8) if the element Peclet number is high. We would like to emphasise that our modified approach is consistent with the original approach proposed in [12], i.e., for both the zero-convection limit (kak ! 0) and for the zeromesh size limit (h ! 0), the original approach is recovered. The implementation of the method is shown in detail in Appendix A. Remark. For the proof of stability, it would have been sufficient to introduce the Peclet-number-based parameter only on the inflow part of the boundary. However, since we introduced the method not only for stability purposes but also to ensure the nopenetration condition on the wall more vigorously, we will not use this refined result in the next sections. The distinction between in-and outflow-parts can be obtained from an extended version of inequality (16):

1 2

Z

Z  2  2 1 a  n uh dX þ a  n uh dX P |ffl{zffl} 2 @ Xin @ Xout P0 Z    h h 1 1 h 2 a  n u dX P  kak u ; u @ X P in 2 @ Xin 2

 2 1 a  n uh dX ¼ 2 @X

Z

þ

This result represents a sufficient condition for the uniqueness     of the solution. Let uh1 ; qh1 and uh2 ; qh2 be two solutions of the simplified model problem. Then, according to the linearity of the problem, the functions uh :¼ uh1  uh2 and qh :¼ qh1  qh2 constitute a solution of the homogeneous variational problem

This extended estimate causes a modification of the last term in (21), which proves stability also if sB is modified only on @ Xin . From this theoretical consideration and the fact that the original approach [12] turned out to be unstable for the channel-flow computations in Section 3, we conclude that the modification of sB according to (8) is an analogon to the stabilising convective term present on inflow boundaries in Nitsche’s approach. For the latter approach, an instability as observed for the present method with a constant parameter does not occur.

26

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

Fig. 2. Solution profiles for

j ¼ 0:0001; a ¼ 1:0 (left) and j ¼ 0:01; a ¼ 1:0 (right).

2.3. Numerical examples Numerical tests are conducted for a one-dimensional outflow boundary-layer problem on the domain X ¼ ½0; 1. The chosen convective velocity is a ¼ 1, and a vanishing source term f ¼ 0 is assumed. The Dirichlet value at the inflow boundary is uðx ¼ 0Þ ¼ 1, and at the outflow boundary uðx ¼ 1Þ ¼ 0. Initially, a diffusivity j ¼ 0:0001 is chosen. The resulting problem is strongly convection-dominated for all discretisations chosen below, characterised by a Peclet number

Pe ¼

a1

j

¼ 10000:

ð27Þ

The results for both choices (7) and (8) of the parameter sB on five equally spaced meshes encompassing 8–512 linearly-interpolated elements are shown in the left half of Fig. 2. For comparison, results obtained with Nitsche’s method as described in [7] are also shown in Fig. 2. It can be observed that the impact of the outflow boundary condition is very small both for the constant-parameter and the Nitsche approach on these relatively coarse meshes. As expected, this is not the case at the inflow boundary. Although not explicitly observable, the boundary value on that side is represented almost exact. The variant parameter (8) puts more emphasis on a fulfilment of the outflow boundary condition. It is obvious that there is a trade-off for this improved representation in that the L2

27

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

convergence is sacrificed for coarse meshes, in order to get control over the boundary value; see Fig. 3. The best choice among these approaches will, to some extent, always be defined by the way the error of a solution is defined by the observer. The best approach will be different if the boundary value and the gradient on the boundary are considered important or if the L2 -error is the only measure for the quality of the solution. With this respect, the method using parameter (8) is located between the strong approach and the Nitsche approach. To gain more insight into the nature of the method, we also studied the problem for a smaller Peclet number Pe ¼ 100 (i.e., j ¼ 0:01). This time, the results are shown in the right half of Fig. 2, again displaying from top to bottom Nitsche’s approach, the present method using a constant parameter and the present method using the Peclet-number-based definition (8). Each time, the same five equally spaced meshes consisting of 8–512 elements are used. Again, we observe that the method with parameter (8) allows for a better representation of the boundary condition on coarse meshes. For the finer meshes, all three approaches show a similar performance. As predicted in the last section, the method with parameter (8) converges to the constant-parameter approach, if the mesh is sufficiently fine (i.e., h is small enough compared to j=kak). This behaviour can also be observed in Fig. 4. Note that also for the error in the H1 norm the method with parameter choice (8) turns out to be beneficial. The better approximation of the boundary value improves the gradient at the outflow boundary and hence decreases the coarse-mesh error in the H1 norm. For both norms, the H1 -norm and the L2 -norm, an almost optimal rate of convergence was obtained for sufficiently refined discretisations. For a specific application, it has to be decided whether it is more appropriate to choose the approach with the smaller L2 -error or the approach with improved control of the boundary value. Dependent on this decision, the proposed method with either a fixed sB according to [12] or a variant sB according to (8) should be chosen. As can be seen from Figs. 3 and 4, one of these two approaches is always superior to Nitsche’s method. In various applications, the improved control of the boundary value is essential for stability or physical reasons. One example is the imposition of a no-penetration boundary condition in wall-normal direction in a turbulent flow problem, as mentioned in the introduction, and as it will be discussed in the next section in more detail. In this case, parameter choice (8) is the key for the success of the method.

3.1. Governing equations and weak form The convective strong form of the incompressible Navier– Stokes equations is given as

@u þ ðu  rÞu  $  ðp1 þ 2meðuÞÞ  f ¼ 0 in X; @t

ð28Þ

$  u ¼ 0 in X;

ð29Þ

where u is the fluid velocity, eðuÞ the (symmetric) rate-of-strain tensor, p the pressure, m the kinematic viscosity, and f a body-force vector. The momentum equation is based on the assumption of an incompressible Newtonian fluid, i.e., on a stress–pressure/velocity relation of the form

ð30Þ

For the Navier–Stokes equations, the stresses play the role of the fluxes in the convection–diffusion equation studied in the previous section. Dirichlet boundary conditions as

u  g ¼ 0 on @ X

g ¼ 0:

ð32Þ

Nevertheless, to simplify the identification of terms enforcing the boundary condition, we will keep g in the equations as long as possible. Based on suitable function spaces, equations (28) and (29) can be restated in a weak form using the L2 -inner products ð; ÞX on X and @ X, respectively.

  @u ;v þ ððu  rÞu; v ÞX  ðp; $  v ÞX þ 2mðeðuÞ; eðv ÞÞX @t X  ðf ; v ÞX  ðr  n; v Þ@X þ ð$  u; qÞX ¼ 0

ð33Þ

A generalised-alpha approach is used for time integration as proposed in [28]. The weak form is discretised in space using a stabilised finite element method. Boundary conditions are imposed weakly, in analogy to what was described for the convection– diffusion equation. In the statement of the equations we will use the residual of the momentum and continuity equation,

   @uh  h þ u  r uh þ $ph  2m$  e uh  f ; @t rhC ¼ $  uh

rhM ¼

ð34Þ ð35Þ

and stabilisation parameters as given in [29,14]:

sC



4 þ uh  Guh þ C I m2 G : G Dt 2 0 !2 11 X X @nj A : ¼ @sM @xi i j

sM ¼

1=2 ð36Þ

ð37Þ

Here, C I is a positive constant and G the second-rank metric tensor,

G¼



T @n @n ; @x @x

ð38Þ

@n which is defined based on the inverse Jacobian @x of the element mapping between the reference (unit) element and the physical domain. The obtained discretised equation reads:

 h            @u þ uh  r uh ; v h X  ph ; $  v h X þ 2m e uh ; e v h X ; vh @t X Xh         f ; v h X  rh  n; v h @ X þ $  uh ; qh X þ sC rhC ; $  v h Xe

3. Wall-bounded turbulent flow

r ¼ p1 þ 2meðuÞ in X:

are assumed, except for periodic boundary conditions at part of the boundary. The pressure in this setting is defined only up to a constant, and has to be fixed for numerical computations. No-slip conditions are assumed at walls:

ð31Þ

þ

 

sM rhM ; rqh

Xe

þ







sM rhM ; uh  r v h

h M rM



 i

e

 Xe







sM rhM  r uh ; v h

 1  h h 1  h  s s  r v Xe þ r ;r X þ p 1; r h X 2m 2m          e uh ; rh X  sB;t g  uh  n; rh @X      ðsB;n  sB;t Þ g  uh  n n  n; rh @ X ¼ 0: h M rM ;





h

 Xe

ð39Þ

The first two lines form the Galerkin part of the equation. It is followed by the least-squares incompressibility, PSPG, and SUPG stabilisation terms in the third line, see [21–24]. The remaining extra terms are the so-called cross-and Reynolds-stress terms in the fourth line, which are taken into account in the residual-based variational multiscale modelling of turbulence described in [14]. In analogy to the convection–diffusion problem, a hybrid approach is used to weakly impose the Dirichlet boundary condition. Extra stresses are introduced on elements next to the boundary. They are determined by the last three lines of the equation. Again, these discontinuous discrete fluxes will be eliminated locally within each element. The introduction of two parameters, sB;n and sB;t , allows for

28

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

Fig. 3. Convergence plot for

j ¼ 0:0001; a ¼ 1:0. (A) Present method with parameter (8), (B) present method with sB ¼ 1 and (C) Nitsche’s method.

Fig. 4. Convergence plot for

j ¼ 0:01; a ¼ 1:0. (A) Present method with parameter (8), (B) present method with sB ¼ 1 and (C) Nitsche’s method.

a separate treatment of tangential and normal direction. We will discuss the corresponding split of the stresses below. For sB;n ¼ sB;t ¼ 1:0, the original form from [12] is recovered. The separate treatment is beneficial especially in the context of wall-bounded turbulent flow. It enables a weak, wall-law-type treatment of the velocities in tangential direction, while a strict enforcement of the condition in normal direction takes care of the required impermeability condition. The improved control of the boundary value in wall-normal direction avoids undesired transport of momentum across the (weakly enforced) boundary. It is also consistent with properties of turbulent channel flows as shown, e.g., in [30]. In the vicinity of the wall, that is, within the viscous sublayer, the root-mean-square values of the velocity fluctuations in tangential direction are known as Oðyþ Þ, while the wall-normal fluctu to scale  ations scale as O ðyþ Þ2 , and thus are much smaller close to the wall. 3.1.1. Decomposition of stresses and tractions An additive split of the additional stresses on elements next to the boundary is considered: h

h con

r ¼r

þr

h bc :

ð40Þ

The first part is a consistency part which accounts for terms arising from partial integration:



rhcon ; rh

 X

      ¼  ph 1; rh X þ 2m e uh ; r h X :

ð41Þ

The second part corresponds to the stresses enforcing the boundary condition:



rhbc ; rh

 X

    ¼ rhbc;t ; r h þ rhbc;n ; r h : X

X

ð42Þ

The latter contains contributions counteracting both the tangential deviations,



rhbc;t ; rh

 X

¼ 2m









sB;t ð1  n  nÞ g  uh  n; rh

 @X

;

ð43Þ

and normal deviations,



rhbc;n ; rh

 X

    ¼ 2m sB;n g  uh  n n  n; rh @X ;

ð44Þ

from the exact boundary condition g. The test functions r h as well as the unknown stresses rh are assumed to be symmetric tensors. The additional stresses and corresponding test functions are defined based on the same basis functions as velocity and pressure; however, they are assumed discontinuous across element boundaries. For velocity and pressure an equal-order interpolation is used, as usual within a stabilised finite element method. Applying the decomposition to the weak form of the momentum equation (33) yields

       rh  n; v h @ X ¼  rhcon  n; v h @ X  rhbc;n  n; v h

@X

   rhbc;t  n; v h : @X

ð45Þ Thus, the traction on the boundary is composed of a consistency part, a normal part enforcing the normal boundary condition and a tangential part enforcing the tangential boundary condition. Remark. This clear interpretation of the additional stresses is a remarkable advantage of the proposed method in comparison to ‘Nitsche-type’ methods. For those methods, a similar split can be done for the additional penalty-type term, see [11]. Nevertheless, those methods involve other additional terms, which also contribute to the overall forcing associated with the boundary

29

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

condition. In particular, the viscous adjoint-consistency term, arising due to the derivative of the test function involved, is responsible for discrete forces acting not only on boundary nodes but also on nodes next to the boundary.

the distance to the wall in wall units,

3.1.2. Determination of the tangential parameter using the law of the wall As indicated above, the part of the boundary condition enforcing the tangential velocities can be interpreted as a Neumann (traction) boundary condition:

uþ :¼



t bc;t ; v h

 @X

rhbc;t ; rh

 Xe

us

m

y;

ð53Þ

and the velocity in wall units, respectively:

u : us

ð54Þ

Using these dimensionless parameters, the law of the wall according to [9] is given by an implicit definition in residual form:

Rwall-law ðus ; y; m; u; j; BÞ

  ¼ rhbc;t  n; v h :

ð46Þ

@X

The discontinuous approximation of the stresses enables a decoupled, element-local computation of the stresses enforcing this condition:



yþ :¼

     ¼ sB;t;e 2m ð1  n  nÞ g  uh  n; r h @ Xe \@ X :

ð47Þ

In our computations, we estimate the ab initio unknown parameter sB;t;e using the formulation of the law of the wall as proposed in [9]. For this purpose, the traction t bc;t is computed using that empirical relation, and the result is compared to the traction resulting from the enforcement of the tangential boundary condition, which is computed based on the most recent (iterative) value of the velocity:

kf t k ¼ t bc;t : areað@ Xe \ @ XÞ

ð48Þ

Here, kf t k denotes the norm of the total tangential force acting on the boundary element. This procedure is in analogy to the one in [10]. The traction on the left-hand side is obtained by solving an element-local, linear system for the element stresses: h

M e rr rhbc;t;e ¼ sB;t;e be ;

þ

¼ yþ  uþ  ejB  eju  1:0  juþ 

ðjuþ Þ2 ðjuþ Þ3  2 6

! ¼0 ð55Þ

Model and material parameters are the von Kármán constant j ¼ 0:4, the log-law constant B ¼ 5:5 and the kinematic viscosity m. Furthermore, u ¼ uh is the local fluid velocity on the boundary and y is an assumed thickness of the unresolved layer. The local fluid velocity is available at each integration point on the surface from the most recent step of the nonlinear iteration that is used to solve the global problem in each time step. Following [10], the assumed thickness of the unresolved layer is taken to be hB =4, where hB is the size of the boundary element in wall-normal direction. It is computed based on the metric tensor G from equation (38) and the outward unit normal n:

2 hB ¼ pffiffiffiffiffiffiffiffiffiffiffiffi : nT Gn

ð56Þ

Thus, Eq. (55) is a scalar, nonlinear equation for us , which can be solved locally at each integration point using a Newton-type itera tion. Once us is determined, the resulting stress t bc;t can be computed using Eq. (52), and the tangential parameter can be obtained via:

ð49Þ

t bc;t  areað@ Xe \ @ XÞ : ~ f t

  h where matrix M e rr results from rhbc;t ; r h . The right hand side be

sB;t;e ¼

is a discrete representation of the tangential part of the stresses normalised by sB;t;e ,

Please note that the denominator may actually be computed independent of sB;t;e , because of the linear relationship between sB;t;e and rhbc;t .

Xe



    2m ð1  n  nÞ g  uh  n; rh @Xe \@X :

Here, uh is the most recent value from the Newton-type nonlinear iteration that is used to solve the nonlinear problem in each time step. System (49) can be used to compute the required part of the ~ hbc;t;e on each element adjacent to discrete normalised stresses r the boundary:

~ hbc;t;e

r

 1 h ¼ M e rr be :

ð50Þ

Using the element’s basis functions, these discrete quantities can be interpolated to the boundary surface and then integrated, yielding the normalised force on the surface:

Z ~ f t ¼

@ Xe \@ X

~h rbc;t  n da:

ð51Þ

Thus the left hand side of Eq. (48) is defined as a linear function of the unknown parameter sB;t;e . In the following, we describe how, first, the traction on the right-hand side of Eq. (48) can be obtained from the law of the wall, and second, the unknown parameter can be determined. For this purpose, we first introduce the friction velocity

s ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi t bc;t ; us :¼

q

ð52Þ

ð57Þ

3.1.3. Setting the wall-normal parameter The wall-normal parameter is chosen based on the same asymptotic behaviour as the one for the convection–diffusion problem in Section 2:

sB;n ¼ 1 þ 2Ree ;

ð58Þ

with the element Reynolds number

Ree ¼

kuC khB

m

:

ð59Þ

Here, kuC k is a characteristic velocity, and hB is the characteristic length of the boundary element in wall-normal direction subject to Eq. (56). For the channel computations in the next section, the characteristic velocity kuC k is determined locally in each element by simply setting it to the maximum velocity on the respective volume element. For practical use, it is still recommended to monitor the size of sB;n and the quality of the resulting approximation, respectively. This can be checked by measuring the flux across the boundary and the fluctuations in wall-normal direction on the boundary. If the fluctuations in wall-normal direction become too large, the calculations are expected to yield results of poor quality. In severe cases, for excessively underresolved boundary layers, it

30

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

might even be observed that they become unstable. We believe that this phenomenon is of the same nature as the well-known outflow instabilities in flows at high Reynolds numbers, occurring at Neumann-type outflow boundaries which have a local recirculation. 3.2. Application of the method to turbulent channel flow The method introduced before is now applied to the wellknown test case of turbulent flow in a plane channel. The problem setup is depicted in Fig. 5. The flow is computed in a volume of size 2p  2  23 p in streamwise, wall-normal, and spanwise direction, respectively (i.e., the channel half-height is d ¼ 1:0). Periodicity is assumed in streamwise and spanwise direction. In this setting, the pressure is determined only up to a constant mode. We solve the corresponding singular linear problems iteratively on a factor space using the discrete projector given in [31] embedded in a GMRES method [32]. With the (kinematic) viscosity m ¼ 0:0001472, and a constant pressure gradient $p ¼ 0:00337204, the friction Reynolds number is Res ¼ ums d ¼ 395. Reference data for this problem based on a direct numerical simulation (DNS) can be found in [33]. A generalised-alpha time-integration scheme with the second-order-accurate parameter setting aM ¼ 5=6; aF ¼ 2=3, and c ¼ 2=3, which corresponds to moderate damping (q1 ¼ 0:5), is chosen. The time-step size is Dt ¼ 0:03. Spatial discretisation of the problem is done using equidistant meshes of 32; 64, and 128 elements in all spatial directions. For comparison, solutions are also provided for meshes which are based on a hyperbolic mesh stretching subject to

h : ½1; 1 ! ½1; 1;

y # hðyÞ ¼

tanh ðC stretch  yÞ : tanh ðC stretch Þ

ð60Þ

The constants for the three resolutions are C stretch ¼ 2:75; 2:3, and 1.9, respectively, chosen such that the first node off the wall is loþ cated at a distance between yþ first ¼ 1:1 and yfirst ¼ 1:3. In contrast, for the unstretched meshes, the first node off the wall is located at yþ first 25; 12, and 6, respectively. This choice clearly does not fulfil usual resolution requirements for boundary layers (yþ 1), see, e.g., [34].

All computational results are given in form of statistics for mean streamwise velocity and velocity fluctuations in all spatial directions, as usual. The mean values are based on averages in space and time, obtained after the solution has reached a statistically steady state (see [17] for details on the sampling procedure). On the unstretched meshes, a comparison of the mixed/hybrid approach both to the ‘Nitsche-type’ formulation proposed in [7], which uses a strong imposition of the no-penetration condition, and to a strong imposition of the no-slip condition is enabled. Results obtained with the mixed/hybrid method on the unstretched mesh are always compared to the respective results on the corresponding stretched mesh. Furthermore, on the coarsest equidistant mesh, an additional result comparing the ‘Nitschetype’ formulation with weakly imposed no-penetration condition to the mixed-hybrid method will be presented. On the coarsest equidistant mesh, 323 , the solution using strong boundary conditions substantially overpredicts the mean velocity in the middle of the channel, see Fig. 6. Both weak formulations cannot be considered accurate either, but they reduce the overprediction by about 50% and get notably closer to the solution on the stretched mesh. For the stretched mesh, the wall-normal size hB of elements on the wall is very small, causing the boundary condition to be enforced rather strictly. The obtained solution improves, being almost equal to the solution using strong boundary conditions on the stretched, problem-suited, boundary-layer mesh. Such results for strong boundary conditions can also be found in [14]. For the stretched mesh, the solutions for weak and strong imposition of the boundary conditions were also reported to be almost equal in [10] for the ‘Nitsche-type’ approach used therein. The mean velocity computed with Nitsche’s method is more accurate in the core of the channel, but the boundary value is better represented by the mixed/hybrid approach. We do not attribute this difference to the way the boundary condition is imposed in wallnormal direction but rather to the way it is imposed in tangential direction. A notable benefit due to the weak imposition of the boundary condition for the predicted fluctuations is not observed. The numerical value obtained for Res was only 370 for Nitsche’s method, while it reached 392 for the mixed/hybrid approach. The severe underprediction observed for Nitsche’s method can be attributed to the way the forces on the boundary are computed.

Fig. 5. Setup for the channel flow problem.

31

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

For our approach, forces on the boundary are computed directly from the contributions of the terms enforcing the boundary condition to the residual, which acts as the right-hand side of the linear incremental system in the global solution process. Hence, non-zero values at boundary nodes are extracted and summed up. In Nitsche’s method, contributions to the residual have non-zero entries not only at boundary nodes but also at nodes next to the boundary. This is due to the adjoint viscous term, as already addressed in a remark in Section 3.1.1. The respective nodal forces are not included in the computation of the total force, causing the imbalance between forcing and computed friction on the wall. The fluctuations in wall-normal direction of the mixed/hybrid approach on the wall are very small, as intended by the design of our normal parameter. For a ‘Nitsche-type’ approach with weakly imposed no-penetration condition, the corresponding fluctuations are contained in Fig. 7. It can be observed that if Nitsche’s method is used also to define the no-penetration condition, a measurable amount of fluctuation of the wall-normal velocity persists on the

boundary. The Reynoldsstress on the wall is predicted about 30 times closer to zero for our newly proposed approach, indicating that for Nitsche’s approach although the method remains stable there is a significantly higher transport of momentum by velocity fluctuations across the ‘impermeable’ boundary. For the medium mesh, 643 , all solutions improve, as expected. As observable in Fig. 8, the weak approaches are still significantly better than the solution using strongly enforced boundary conditions on an equidistant mesh. Values obtained for Res on the equidistant mesh are 380 for Nitsche’s method and 393 for the mixed/hybrid method on this discretisation. A refinement to 1283 elements shows that all solutions further approach the DNS solution. The weak approaches on the equidistant mesh now provide a solution close to the one obtained by our weak method on the stretched mesh, which is within 4% of the DNS solution for the mean value at the core of the channel. The strong approach on the unstretched mesh deviates by about 10%. For such fine resolutions, differences between weak and strong imposition of

40 u+ 35 30 25 20 15 DNS 10

strong Dirichlet Nitsche’s method

5

Mixed/Hybrid Dirichlet Mixed/Hybrid Dirichlet (stretched mesh)

0 0

50

100

150

200

250

300

350

y+

5 DNS strong Dirichlet Nitsche’s method

4

Mixed/Hybrid Dirichlet Mixed/Hybrid Dirichlet (stretched mesh) 3 +

rms u

2 rms w+ 1 rms v+ 0 0

50

100

150

200

250

300

350

y+

Fig. 6. Mean streamwise velocity and rms of velocity fluctuations (u: streamwise, v: wall-normal, w: spanwise) in wall coordinates for 323 meshes at Res ¼ 395.

32

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

boundary conditions as well as between stretched and unstretched meshes become increasingly smaller, as expected.

0.9 +

rms v

0.8

4. Conclusions and outlook

0.7

0.6

0.5 Mixed/Hybrid Nitsche’s method 0.4

0.3

0.2

0.1

0 0

100

200

y

300

+

Fig. 7. Rms of velocity fluctuations in wall-normal direction in wall coordinates for an equidistant 323 mesh at Res ¼ 395. Both solutions displayed here impose the nopenetration condition weakly.

In this work, an embedded Dirichlet formulation has been proposed for weakly enforcing boundary conditions for convectiondominated convection–diffusion and turbulent flow problems. For turbulent flow, tangential and normal components of the traction enforcing the boundary condition are separately addressed, allowing for a physically-motivated distinction in the respective definition of the boundary condition. A wall-law based formulation was invoked for the tangential component, while the normal component was treated by a parameter derived from theoretical considerations for convection–diffusion problems at high Peclet number. Wall impermeability both at high Reynolds numbers and on coarse meshes is guaranteed, and thus, accurate and robust computations in practical applications are enabled. With this respect, the approach presented in this paper is clearly superior to the original approach proposed in [12] which was observed to fail for the turbulent channel flow benchmark on the coarsest mesh investigated. The impact of the modification was studied for convection–diffusion problems, showing an improved prediction of boundary values. Our stability analysis showed further an analogy between our modified approach and the way stability is enhanced in the ‘Nitsche-type’ approach by an extra advection-controlled term on the inflow. For refined meshes, the present approach consistently converges towards the constant-parameter embedded

30 u+ 25

20

15 DNS 10

strong Dirichlet

5

Mixed/Hybrid Dirichlet

Nitsche’s method Mixed/Hybrid Dirichlet (stretched mesh) 0 0

50

100

150

200

250

300

350

y+

4 DNS 3.5

strong Dirichlet Nitsche’s method

3 rms u

+

Mixed/Hybrid Dirichlet

2.5

Mixed/Hybrid Dirichlet (stretched mesh)

2 rms w+

1.5 1

rms v+

0.5 0 0

50

100

150

200

250

300

350

y+

Fig. 8. Mean streamwise velocity and rms of velocity fluctuations (u: streamwise, v: wall-normal, w: spanwise) in wall coordinates for 643 meshes at Res ¼ 395.

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

Dirichlet formulation initially proposed for interfaces in [12]. We would like to emphasise that the inability of that constant-parameter formulation to impose boundary values on outflow boundaries becomes problematic only if particularly sharp boundary layers are present. If it is used as an embedded boundary condition in the interior of the domain, where no sharp boundary layers are present, that approach provides satisfactory results. This has been demonstrated in [12]. For turbulent channel flow superior performance of the weak formulation compared to a strong imposition of wall boundary conditions has been observed. On coarse meshes, the overprediction of the mean streamwise velocity is reduced by almost 50%. For refined meshes, the weak solution approaches the strong solution. We have also compared the proposed mixed/hybrid method in which all components of the no-slip condition are weakly imposed to a weak imposition of the boundary condition using Nitsche’s method as proposed in [7], where a strong imposition of the boundary condition in wall-normal direction was used. We have been able to identify advantages and disadvantages for both approaches. On the one hand, Nitsche’s method predicted 5–10% more accurate values for the mean velocity in the core of the channel. On the other hand, boundary values were predicted more accurately by the mixed/hybrid method. It remains to be investigated whether the assumption of taking the assumed distance to the wall to be one quarter of the size of the first element, which was the choice for Nitsche’s method in [10], is also the best choice for the present mixed/hybrid method. Major advantages of the mixed/hybrid approach are, on the one hand, that the associated forces are localised on the boundary, and on the other hand, that the physical interpretation of the boundary condition as a traction is particularly evident. We consider these advantages beneficial for applications of the method in the context of turbulent fluid–structure interaction problems we intend to pursue in the future. Further improvement of the present method in the future might address the fact that, up to now, ‘convective’ Reynolds stresses at the wall are suppressed by the strict enforcement of the wall-normal velocity. For the resolutions and moderately high Reynolds numbers investigated in this study, this appears to be a valid assumption. Nevertheless, when going to even higher Reynolds numbers, additional input to the wall-layer modelling process might be generated from RANS-based models in the near-wall region as in other state-of-the-art wall-layer modelling approaches (see, e.g., [5] for an overview).

33

Appendix A. Implementation for the convection–diffusion problem An algorithmic realisation of the hybrid/mixed Dirichlet approach elaborated on in Section 2 is provided. To simplify a generalisation to nonlinear problems, an incremental framework is assumed for the finite-element solution of equation (5). In the following, the term ‘location vector’ will refer to the mapping connecting local element degrees of freedom to their global numbers. Table A.1 shows how the weak imposition of the boundary condition is included in a general stabilised finite element solution process for the convection–diffusion equation. It contains a loop over all boundary elements, in which the Dirichlet condition is actually imposed. This loop is described in more detail in Table A.2. The algorithm shows how the additional degrees of freedom for the fluxes are eliminated via an element-local static condensation procedure. This is ensured due to the fact that the local problems for the fluxes are decoupled as a result of the discontinuous basis functions. Consequently, the extra fluxes do not occur on the global level outside the element call. All data required for the computation of the fluxes, like for instance the shape functions, is already available on the parent (volume) element. Thus, the mixed/hybrid approach does not require more additional data structures than the primal Nitsche approach, even in parallel. However, for optimal performance, we recommend to include the boundary elements in the load balancing scheme in parallel computations. For an efficient implementation we furthermore recommend to exploit the special block-diagonal structure of the element ‘mass’-matrix of fluxes for its inversion in the static condensation procedure. For convenience, the boundary integrals are evaluated for all basis functions of the parent element. This way, several zeros are computed for basis functions associated with inner nodes, but the local assembly procedure is simplified, because the matrices K e unr and K e qnv as well as the corresponding vectors f e $ur and f e unr already have the correct size for the matrix–matrix and matrix–vector multiplications involved in the computation of Be and f e . Note that the basis functions for both problem variable and flux are based on the same basis functions. They have to be evaluated once at every integration point of the volume element and once at every integration point of the boundary element.

Appendix B. Implementation for flow problems Acknowledgement The support of the second author by the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged.

The implementation of the mixed-hybrid boundary condition for flow problems is very similar to what has been shown for the convection–diffusion problem in Appendix A. In every Newton step of the solution procedure of the nonlinear flow problem, an

Table A.1 Finite element algorithm for the solution of the convection–diffusion equation based on a weak mixed/hybrid imposition of Dirichlet boundary conditions. set initial values uh ¼ uh0 loop all elements e in X do Gauss-loop for integrals on element e:       f e; Ke $ a  $uh ; v h Xe þ j $uh ; $v h Xe  f ; v h Xe þ   þ a  $uh  jDuh  f ; sa  $v h Xe do assembly based on the element’s location vector obtain system matrix K obtain right-hand side f do boundary element loop to impose boundary condition (Table A.2) solve the linear system K  Duh ¼ f to compute increment Duh update solution by increment uh ¼ uh0 þ Duh do output

34

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35 Table A.2 Boundary element loop for imposing mixed/hybrid weak Dirichlet boundary conditions for convection–diffusion problems. loop all boundary elements b on @ X get parent (volume) element e owning boundary element b do Gauss-loop for volume integrals on element e: M e qr

$ þ j1 ðqh ; r h ÞXe $ ð$uh ; r h ÞXe

f e $ur ; K e $ur  1 do matrix inversion and compute M e qr

do loop for boundary integrals on boundary element b:   $  qh  n; v h @ Xe \@ X  h   $ þsB u  g n; rh @ Xe \@ X

K e qnv f e unr ; K e unr compute matrix

 1   Be ¼ K e qnv M e qr K e $ur þ K e unr

and right-hand side contribution f

e

 1  ¼ K e qnv M e qr f e $ur þ f

e

unr



do assembly based on location vector of parent element: add Be to system matrix K add f e to right-hand side f of linear system

Table B.3 Boundary element loop to impose mixed/hybrid weak Dirichlet boundary conditions based on the law of the wall according to [9] for fluid problems. loop all boundary elements b on @ X get parent (volume) element e owning boundary element b do Gauss-loop for volume integrals on element e: M e rr f e pr ; K e pr f e eðuÞr ; K e eðuÞr  do matrix inversion and compute M e rr

$ þ 21m ðrh ; r h ÞXe $ þ 21m ðph 1; rh ÞXe   $ ðe uh ; rh ÞXe 1

determine normalised tangential stresses by integration on b:     h $ ð1  n  nÞ g  uh  n; rh @ Xe \@ X be  1 h 1 compute sB;t;e  rhbc;t;e ¼ M e rr be and 1  rhbc;t  n by quadrature evaluate sB;t;e L2 ð@ Xe \@ XÞ

do loop for boundary integrals on boundary element b:   K e rnv $  rh  n; v h @ Xe \@ X compute t bc;t using Spalding’s law and set

sB;t;e ¼ k bc;t k 1 t

areað@ Xe \@ XÞ

sB;t;e kf t k

K e urn ; f e urn

sB;n;e ¼ sB;n;e ðRee Þ, see Eq. (58)     $  sB;t g  uh  n; r h @ Xe \@ X ,      sB;n  sB;t g  uh  n n  n; r h @ Xe \@ X

multiply matrices to obtain the linearisations with respect to u  1   Be u ¼ K e rnv M e rr K e eðuÞr þ K e urn multiply matrices to obtain the linearisations with respect to p  1   Be p ¼ K e rnv M e rr K e pr compute right-hand side contribution  1   f e ¼ K e rnv M e rr f e eðuÞr þ f e urn þ f e pr do assembly based on location vector of parent element: add Be u and Be p to system matrix K add f e to right-hand side f of linear system

algorithm analogous to Table A.1 is employed. The Gauss-loop for the integrals on the element e is now used to compute a right-hand side vector f and a tangent matrix K for the weak form of the Navier–Stokes equations. As usual, they include contributions from the SUPG, PSPG and least-squares incompressibility stabilisation terms, and if required, the additional cross-and Reynolds stress

terms from the residual-based variational multiscale modelling of turbulence. Major changes appear only on the boundary element loop. Here, the separate treatment of the tangential and the normal parameter, the pressure, and the usage of the law of the wall have to be accounted for. Table B.3 provides an overview on this boundary element loop.

P. Gamnitzer et al. / Comput. Methods Appl. Mech. Engrg. 245–246 (2012) 22–35

References [1] J.W. Deardorff, A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers, J. Fluid Mech. 41 (1970) 453–480. [2] U. Schumann, Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli, J. Comput. Phys. 18 (1975) 376–404. [3] W. Cabot, P. Moin, Approximate wall boundary conditions in the large-eddy simulation of high Reynolds number flow, Flow Turbul. Combust. 63 (1999) 269–291. [4] U. Piomelli, E. Balaras, Wall-layer models for large-eddy simulation, Annu. Rev. Fluid Mech. 34 (2002) 349–374. [5] U. Piomelli, Wall-layer models for large-eddy simulations, Prog. Aerosp. Sci. 44 (2008) 437–446. [6] J. Fröhlich, D. von Terzi, Hybrid LES/RANS methods for the simulation of turbulent flows, Prog. Aerosp. Sci. 44 (2008) 349–377. [7] Y. Bazilevs, T.J.R. Hughes, Weak imposition of Dirichlet boundary conditions in fluid mechanics, Comput. Fluids 36 (2007) 12–26. [8] J. Nitsche, Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind, Abh. Math. Semin. Univ. Hambg. 36 (1971) 9–15. [9] D.B. Spalding, A single formula for the law of the wall, J. Appl. Mech. 28 (1961) 444–458. [10] Y. Bazilevs, C. Michler, V.M. Calo, T.J.R. Hughes, Weak Dirichlet boundary conditions for wall-bounded turbulent flows, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4853–4862. [11] Y. Bazilevs, C. Michler, V.M. Calo, T.J.R. Hughes, Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes, Comput. Methods Appl. Mech. Engrg. 199 (2010) 780–790. [12] A. Gerstenberger, W.A. Wall, An embedded Dirichlet formulation for 3D continua, Int. J. Numer. Methods Engrg. 82 (2010) 537–563. [13] V.M. Calo, Residual-based Multiscale Turbulence Modeling: Finite Volume Simulations of Bypass Transition, Ph.D. Thesis, Stanford University, 2005. [14] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (2007) 173–201. [15] V. Gravemeier, W.A. Wall, Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number, Int. J. Numer. Methods Fluids 65 (2011) 1260–1278. [16] R. Codina, J. Principe, O. Guasch, S. Badia, Time dependent subscales in the stabilized finite element approximation of incompressible flow problems, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2413–2430. [17] P. Gamnitzer, V. Gravemeier, W.A. Wall, Time-dependent subgrid scales in residual-based large eddy simulation of turbulent channel flow, Comput. Methods Appl. Mech. Engrg. 199 (2010) 819–827.

35

[18] V. Gravemeier, W.A. Wall, An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number, J. Comput. Phys. 229 (2010) 6047–6070. [19] V. Gravemeier, M.W. Gee, M. Kronbichler, W.A. Wall, An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow, Comput. Methods Appl. Mech. Engrg. 199 (2010) 853–864. [20] V. Gravemeier, M. Kronbichler, M.W. Gee, W.A. Wall, An algebraic variational multiscale-multigrid method for large-eddy simulation: generalized-alpha time integration, Fourier analysis and application to turbulent flow past a square-section cylinder, Comput. Mech. 47 (2011) 217–233. [21] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equation, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [22] P. Hansbo, A. Szepessy, A velocity–pressure streamline diffusion finite element method for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 84 (1990) 175–192. [23] T. Tezduyar, S. Mittal, S.E. Ray, R. Shih, Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity–pressure elements, Comput. Methods Appl. Mech. Engrg. 95 (1992) 221–242. [24] M. Braack, E. Burman, V. John, G. Lube, Stabilized finite element methods for the generalized Oseen problem, Comput. Methods Appl. Mech. Engrg. 196 (2007) 853–866. [25] R. Codina, Stabilization of incompressibility and convection through orthogonal subscales in finite element methods, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1579–1599. [26] L.P. Franca, F. Valentin, On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1785–1800. [27] J. Baiges, R. Codina, F. Henke, S. Shahmiri, W.A. Wall, A symmetric method for weakly imposing Dirichlet boundary conditions in embedded finite element meshes, Int. J. Numer. Methods Engrg. 90 (2012) 636–658. [28] K.E. Jansen, C.H. Whiting, G.M. Hulbert, A generalized-method for integrating the filtered Navier–Stokes equations with a stabilized finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2000) 305–319. [29] C.A. Taylor, T.J.R. Hughes, C.K. Zarins, Finite element modeling of blood flow in arteries, Comput. Methods Appl. Mech. Engrg. 158 (1998) 155–196. [30] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [31] P. Bochev, R.B. Lehoucq, On the finite element solution of the pure Neumann problem, SIAM Rev. 47 (2005) 50–66. [32] Y. Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856– 869. [33] R.D. Moser, J. Kim, N.N. Mansour, Direct numerical simulation of turbulent channel flow up to Res = 590, Phys. Fluids 11 (1999) 943–945. [34] P. Sagaut, Large Eddy Simulation for Incompressible Flows, third ed., Springer, 2006.