Implementation of k-Exact Finite Volume Reconstruction in DLR’s Next-Generation CFD Solver: Flucs and its Comparison to Other Methods

(1)

,

STOCKHOLM SWEDEN 2017

Implementation of k-Exact Finite

Volume Reconstruction in DLR’s

Next-Generation CFD Solver: Flucs

and its Comparison to Other

Methods

SIDDHANT AGARWAL

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ENGINEERING SCIENCES

(2)

(3)

Acknowledgements

I would like to thank Dr. Philipp Schlatter for agreeing to supervise this thesis at KTH Royal Institute of Technology. I am also thankful to Dr. Norbert Kroll at DLR’s C2A2S2E (Center for Computer Applications in Aerospace Science) for receiving my application and forwarding it to relevant department heads.

Thanks to Dr. Ralf Hartmann - team leader of CFD algorithms - for supervising my thesis at DLR. He, along with Ms. Anne-Kathrin Kaiser at the International Office of Technische Universität Braunschweig (host ERASMUS university) also helped me wade through a near insurmountable pile of immigration paperwork and bureaucracy before and during my thesis in Germany. I am deeply indebted to Tobias Leicht and Axel Schwöppe of the CFD algorithms group for their relentless support throughout the thesis. With their never-ending patience, they helped me whenever I had C++ questions or needed to double-check a mathematical concept in the algorithms. Their feedback on the drafts of the thesis has also been very helpful.

Finally, I am grateful to my parents Neeru and Sameer Agarwal for their support in just about every way imaginable.

(4)

Reconstruction in DLR’s Next-Generation

CFD Solver: Flucs and its Comparison to

Other Methods

Siddhant Agarwal

∗

Master’s Thesis at DLR for KTH

December 14, 2017

Abstract

This thesis extended the order of the reconstruction of state for convective fluxes used by Finite Volume (FV) algorithm in DLR’s next-generation CFD solver: Flucs, from constant and linear to quadratic and cubic. Two approaches for calculating derivatives were implemented in Flucs and some test cases were tried. To allow for integration of moments within each cell and a higher-order integration of fluxes, the mesh used by Discontinuous Galerkin (DG) was fed to the FV algorithm. Insufficient geometric treatment of the boundary cells and the dummy cells in FV is believed to be detrimental to the order of error reduction in NACA0012 case and the smooth bump case. In the smooth bump case, the FV algorithms failed to show higher than second order error reduction because of this reason. The order of the schemes away from the boundaries was verified with the Ehrenfried Vortex test case. For at least structured meshes and unstructured meshes with quads, schemes of order k approached k+1 order accuracy on sufficiently fine meshes. The original goal of this thesis was partly accomplished and some further work in the code is expected. Detta examensarbete vidareutvecklade ordningen av rekonstruktion av tillståndet för konvektiva flöden som används av FV-algoritmen i DLR’s nästa generations CFD-lösare: Flucs, från konstant och linjär till kvadratisk och kubisk. Två metoder för beräkning av derivatorna genomfördes i Flucs och några testfall försökades. För att möjliggöra integrationen av momenten inom varje cell och integrationen av flöden på en högre ordning, matades det nät som används av DG till FV-algoritmen. Otillräcklig geometrisk behandling av gränscellerna och dummycellerna antas vara skadlig för ordning av felreduceringen i NACA0012-fallet och det Smmoth Bump-fallet. I det Smooth Bump-Fallet fick FV-algoritmerna inte visa högre än andra ordning av felreducering på grund av denna anledning. Ordningen av schemana bort från gränserna verifierades med Ehrenfried Vortex-testfallet. För åtminstone strukturerade nät och ostrukturerade nät med fyrkantiga celler nådde algorithmer med order k, k+1 ordernoggrannhet på tillräckliga bra nät. Det ursprungliga målet för denna avhandling uppnåddes delvis och ytterligare arbete i koden förväntas.

I. Introduction

The German Aerospace Center (DLR - Deutsches Zentrum für Luft- und Raumfahrt) is developing its next-generation solver called Flucs: Flexible Unstructured CFD Software. It is a part of the wider effort to move towards "virtual aircraft design and flight testing based on high-fidelity methods". The goal is to provide a platform for carrying out multidisciplinary simulations such as maneuver load simulations [15].

Accurate maneuver load simulations require the use of CFD simulations coupled with Computational Solid Mechanics (CSM) to account for structural defor-mations of flexible structures and their impact on the aerodynamic forces. The CFD-CSM loops, in addition,

∗_{agsiddhant@gmail.com}

Figure 1: Overall architecture of the FlowSimulator

also need a flight dynamics module for trim computations. Such a framework is provided through FSDM, which stands for Flow Simulator Data Manager [20]. Figure 1 shows some of the modules included.

(6)

of CFD solvers: ONERA’s elsA [5] or DLR’s TAU [8]. While TAU is currently used for many industry and research applications, the need for a new CFD solver arose due to three deciding factors: [16]

1. Full exploitation of new improvements in HPC (High Performance Computing) hardware

2. Need for modular code development for a highly flex-ible and maintainable software that can incorporate additions necessary for multidisciplinary simulations 3. Difficulty of integrating implicit solver techniques and

higher order discretizations into legacy codes Thus, began the development of Flucs. Just as TAU, Flucs will also eventually be available in the FSDM environment. The key focus of Flucs is the ability to fully exploit the existing and future parallel HPC infrastructures [14]. The parallelization occurs at two different levels - a domain decomposition where there is an overlap of communication and a parallel processing of a domain via shared-memory [16].

Flucs is being developed with state-of-the-art cod-ing practices. It has a highly modular structure with high level of abstraction, making future additions to the code convenient [16]. The version control of the code is made possible via git, enabling multiple developers to work and collaborate on the code at the same time [9].

Finally, Flucs will also boast some advanced nu-merical methods such as higher-order space and time discretization methods, that are often difficult to im-plement in legacy codes and an implicit solver which utilizes Automatic Differentiation (AD) for local resid-ual constribution when calculating a compact-stencil derivative-approximation for the Jacobian. The robust stability of implicit methods compared to explicit methods is well-known in the CFD community. There is also an increasing interest in higher-order methods for their accuracy and efficiency. It is even crucial, for example in LES, where higher than second order spatial and temporal accuracy is needed [21].

The purpose of this thesis is to implement a k-exact Finite-Volume reconstruction in Flucs. This will be an addition to the already existing spatial discretization schemes of Discontinuous Galerkin (DG) for an arbitrary degree and Finite Volume (FV) schemes of constant and linear reconstruction. The linear reconstruction scheme has two variants for calculating gradients, namely Green-Gauss and least-squares.

The layout of this paper is as follows. First, some

theory of the flow equations, the turbulence model and k-exact finite volume reconstruction are provided. Then, the implementation of the k-exact reconstruction algo-rithms in Flucs is documented. Finally, some results on the newly implemented algorithms and their comparison to already existing discretization schemes are presented. Some possible future improvements and additions are also explained.

II. Theory

i.

Governing Equations

The system of equations to be solved here is the Reynolds Averaged Navier-Stokes (RANS) equations coupled with the negative variant of Spalart Allmaras (SA-neg) turbu-lence model [1]. RANS models such as this one are an approximated, time-averaged version of the Navier-Stokes equations. ∂U ∂t + ∇ · ( −→ Fc(U) − −→ Fv(U,∇U)) =S(U,∇U) (1) In equation 1, the conservative variables U and the inviscid flux vectors are

U=          ρ ρu ρv ρw ρE ρ ˜ν          , Fcx=          ρu ρu2+p ρuv ρuw ρHu ρu ˜ν          , Fcy=          ρv ρuv ρv2+p ρvw ρHv ρv ˜ν          , (2) Fcz=          ρw ρuw ρvw ρw2+p ρHw ρw ˜ν          , (3)

where, H is the enthalpy and is equal to E+_ρp. The viscous flux vectors are

Fvx=             0 τxx τyx τzx uτxx+vτyx+wτzx+cp(Prµ + µT PrT )∂T ∂x 1 σ(µ+fnρ ˜ν) ∂ ˜ν ∂x             , (4) Fvy=             0 τxy τyy τzy uτxy+vτyy+wτzy+cp( µ Pr+ µT PrT )∂T ∂y 1 σ(µ+fnρ ˜ν) ∂ ˜ν ∂y             , (5)

(7)

Fvz =             0 τxz τyz τzz uτxz+vτyz+wτzz+ +cp( µ Pr+ µT PrT )∂T ∂z 1 σ(µ+ fnρ ˜ν) ∂ ˜ν ∂z             , (6)

The total viscous stress tensor τijis given by

τij=2(µ+µT)˙γij (7)

the shear rate tensor by

˙γij= 1 2 ∂Vi ∂xj +∂Vj ∂xi ! −1 3 ∂V_k ∂xk δij (8)

and the turbulent eddy viscosity by

µT=

ρ ˜ν fv1 ˜ν≥0

0 ˜ν<0 (9)

The source term S represents the turbulence model. The negative variant of the Spalart-Allmaras model is used due to its appropriateness for higher-order discretizations [1]. The source term, as follows,

S=           0 0 0 0 0 ρ(P−D) + 1 σcb2ρ∇˜ν· ∇˜ν− 1 σ(ν+˜ν)∇ρ· ∇˜ν           (10) comes from the equation for the evolution of turbulence working variable ˜ν: ∂ρ ˜ν ∂t + ∇ · (ρ ˜ν ~ V) − ∇ · (µ+fnρ ˜ν)∇˜ν σ =ρ(P−D) + 1 σcb2ρ∇˜ν· ∇˜ν− 1 σ(ν+ ˜ν)∇ρ· ∇˜ν (11)

In equation 11,V is the vector of Cartesian components~ of velocity, µ is the dynamic viscosity and ν the kinematic viscosity. P and D are production and destruction, respec-tively: P= ( cb1(1− ft2)S˜ ˜ν ˜ν≥0, cb1(1−ct3)S˜˜ν ˜ν<0 D=        µ0 cw1fw−cb1 κ2 ft2  ˜ν d 2 ˜ν≥0, −µ0cw1 ˜ν d 2 ˜ν<0 (12)

If S is the magnitude of vorticity, then modified vorticity is given by: ˜ S=      S+S¯ S¯ ≥ −cv2S S+S(c 2 v2S+cv3S¯) (cv3−2cv2)S−S¯ ¯ S< −cv2S ¯ S= ν 0_{˜ν f} ν2 κ2d2 (13)

Here d is he distance to the closest wall and f functions are defined as:

fv1 = χ 3 χ3+c3_v1, fv2 =1− χ 1+χ fv1 fn = cn1 +χ3 cn1−χ3, ft2=ct3 exp(−ct4χ 2₎ χ= µ 0_˜ν ν (14)

Note that the trip term for transition is not modeled in Flucs (ft1 = 0), while the ft2term is included. The de-struction term coefficient is:

fw=g 1+c6 w3 g6₊_c6 w3 , fv2=1− χ 1+χ fv1 g=r+cw2(r6−r) r=min µ0˜ν ˜ Sκ2_d2, 10 (15)

The values of constants used are: σ = 0.66, κ = 0.41, cv1 = 7.1, cv2 = 0.7, cv3= 0.9, cb1 = 0.1355, cb2 = 0.622, cw1 = cb1 κ2 + 1+cb2 σ , cw2 = 0.3, cw3 = 2.0, ct3 = 1.2, ct4=0.5, cn1=16.

The ideal gas equation is further used to introduce the relation between energy E and pressure p

E= p γ−1+ 1 2ρ(u 2₊_v2₊_w2₎ ₍₁₆₎ Here, γ=1.4, Pr=0.72 and PrT=0.9.

ii.

Spatial Discretization

The Finite Volume method numerically solves the con-servative form of equations. Hence, equation 1 can be integrated to get d dtU¯i+ 1 Vi Z Ωi ∇ · (~Fc− ~Fv)dV = Z Ωi S dV (17) where Ωi is the control volume. The definition of con-trol volume depends on whether the PDEs are discretized

(8)

Figure 2: A 2D control colume to integrate over

using cell-centered or node-centered formulations. Flucs uses the cell-centered method, i.e. the values of state vari-ables are stored at a discrete point located at the barycenter of each control volume. That means that a single control volumeΩilooks like, in the 2D case, figure 2. The control volumes are then, the mesh elements themselves which are bounded by a certain number of faces. A tetrahe-dron, for example consists of 4 faces. The average of the conservative variables is ¯ Ui = 1 Vi Z Ωi U(x, y, z)dV (18) The divergence of the volume integral of flux in equa-tion 17 is replaced with a surface integral through the divergence theorem: d dtU¯i+ 1 Vi I ∂Ωi (~Fc− ~Fv) · ~n dS= Z Ωi S· dV (19)

where~n is the unit vector (see figure 2) acting outward and normal to the face ∂Ωi and

H

∂Ωis the surface integral

over the volumeΩi.

The surface integral in equation 19 can be approxi-mated as the summation of integrals on all control volume faces: 1 Vi H ∂Ωi(~Fc− ~Fv) · ~n dS≈ NF ∑ j=1 R Ωi,j(~Fc− ~Fv) · ~n dS = NF ∑ j=1 N_G ∑ φ=1 ωg,φ~FcUL(~xj,φ), UR(~xj,φ) · ~n dS ! −N∑F j=1 NG ∑ φ=1 ωg,φ~FvU(~xj,φ),∇U(~xj,φ) · ~n dS ! (20)

where, NFis the number of faces for cell i and NG is the number of Gaussian quadrature points on the jthface of cell i, used to numerically integrate flux and ωg,φ are the integration weights at each quadrature point φ.

The number of quadrature points required depends on the desired order of accuracy and shape of the element.

Table 1: Number of Gaussian quadrature points needed. Triangle

Order k = 1 k = 2 k=3

Surface Integral 1 2 2

Volume Integral 1 3 4

Figure 3: First order accurate reconstruction; k = 0

These are summarized in table 1 [17].

Calculation of the convective flux ~Fc in equation 20 requires the state variables: ULon the left and URon the right side of each face j at each quadrature point φ.

ii.1 Reconstruction

The easiest way to obtain UL and URwould be to assume that are the same as Ui and Uj, respectively as shown in figure 3. This method would be first-order accurate. The next-order accurate method is to assume that the solution varies linearly within a control volume and hence, reconstruct it as a linear function, as shown in figure 4 [4]. Then, the left and right states of the solution can be written as:

UL =Ui+ ∇Ui·~rL

UR=Uj+ ∇Uj·~rR (21)

(9)

where,∇Uiand∇Uj are the gradients computed at cell centers i and j, respectively using either Green Gaus or the Least Squares method in the case of FV discretization. A modified version is the Barth and Jespersen method [3]:

UL=Ui+ψI(∇Ui·~rL)

UR=Uj+ψI ∇Uj·~rR (22) where the limiter functions ψI and ψJ can be used to prevent oscillations or spurious solutions in areas where high derivatives, such as shocks, exist [4].

ii.2 k-Exact Reconstruction

Gooch has presented a step by step approach for carrying out a k-exact FV reconstruction [19], and is briefly sum-marized here for completeness. A reconstruction can be made k-exact or (k+1)-order accurate if, the control vol-ume average ¯uiof any state variable u at reference location i is replaced by a Taylor expansion, as shown in equation 23, which is carried out through the kth derivative.

u_iR(~x) = k ∑ l=0 { 1 l![ l ∑ m=0 l! m!(l−m)! (x−xi)_∂x∂ l−m _m ∑ n=0 m! n!(m−n)! (y−yi)_∂y∂ m−n (z−zi)_∂z∂ n ]u|i } (24) where, k is the parameter specifying k-exact reconstruc-tion. k=1 is second-order accurate, k=2 is third-order accurate and so on.

To reconstruct ui in a control volume using equa-tion 24, the derivatives are needed. These derivatives can be computed by setting up a least squares problem, where an over-determined system of equations is solved. To build this system of equations, the conservation

Figure 5: Reconstruction stencils for an arbitrary (k+1)-order

accu-racy.

of the mean property is used, which requires that:

¯ ui= 1 Vi Z Ωi u_iR(~x− ~xi)dV (25)

where ¯ui is the exact average of cellΩ with a volume Vi and the reconstruction polynomial is u_iR. The right hand side can be expanded:

Here, the moments xn_ym_zl

ican be defined as follows:

xn_ym_zl i ≡ 1 Vi Z Ωi (x−xi)n(y−yi)m(z−zi)ldV (27)

This method of calculating the derivatives using least-squares requires a reconstruction stencil {Vj}i that can extend to neighboring cells and beyond, depending on the order of accuracy desired. The stencil can be obtained by recursively adding the neighboring cells, i.e. in the von Neumann neighborhood. For a 2D case, the stencil would look like figure 5. For example, for k=3, the fourth-order accurate stencil extends from cell i all the way up to cells marked with 4.

(10)

over a single control volume Vjis: 1 Vj R Ωju R i (~x)dV = u|i+ ∂u ∂x|i ¯xij+ ∂u ∂y|i ¯yij+ ∂u ∂z|i ¯zij +∂ 2_u ∂x2|i ¯x2_ij 2 + ∂2u ∂x∂y|i xyij + ∂ 2_u ∂x∂z|i xzij+ ∂2u ∂y2|i y2_ij 2 + ∂ 2_u ∂y∂z|i yzij+ ∂2u ∂z2|i z2_ij 2 +... (28) xn_ym_zl ij≡ 1 Vj Z Ωj (x−xi)n(y−yi)m(z−zi)ldV (29) In order to avoid computing moments of each control volume in the reconstruction stencil{Vj}iabout~xi,(~x−~xi) is replaced with(~x− ~xj) + (~xj− ~xi)modifying equation 28 into: 1 Vj R Ωju R i(~x)dV= u|i+ ∂u ∂x|i ¯xj+ xj−xi +∂u ∂y|i ¯yj+ yj−yi +∂u ∂z|i ¯zj+ zj−zi +∂ 2_u ∂x2|i x2 j +2xj xj−xi − xj−xi 2 2 + ∂ 2_u ∂x∂y|i(xyj+xj yj−yi + xj−xi yj + xj−xi yj−yi )+ .. . +∂ 2_u ∂z2|i z2 j +2zj zj−zi − zj−zi 2 2 +... (30)

1 V_i R Ωiu R i (~x− ~xi)dV = k ∑ l=0 { 1 l![ l ∑ m=0 l! m!(l−m)! b xij_∂x∂ l−m _m ∑ n=0 m! n!(m−n)! b yij_∂y∂ m−n b zij_∂z∂ n ]u|_i } (32) where, \ xmx_ymy_zmz ij ≡ 1 Vj R Ωj x−xj + xj−xi mx · y−yj+ yj−yimy · z−zj + zj−zi mz_dV = m_∑z kz=0 my ∑ ky=0 mx ∑ kx=0 mz! kz!(mz−kz)! · my! ky! my−ky! mx! kx!(mx−kx)! xj −xikx · yj−yiky · zj−zikz ·x(mx−kx)_y(my−ky)_z(mz−kz) j (33) Hence, the weighted least squares problem can be written as A·                       u ∂u ∂x ∂u ∂y ∂u ∂z 1 2 ∂2u ∂x2 ∂2u ∂x∂y .. .                       i =           ui wi1u1 wi2u2 .. . wiNuN           (34)

where, the matrix A is equal to            1 xi yi zi x2i . . . wi1 wi1bxi1 wi1byi1 wi1bzi1 wi1xb2i1 . . . wi2 wi2bxi2 wi2ybi2 wi2bzi2 wi2xb 2 i2 . . . .. . ... ... ... ... . .. wiN wiNbxiN wiNbyiN wiNbziN wiNxb2iN . . .            (35)

The system of equations presented in 34 and 35 can be solved to obtain the derivatives. N is the number of neigh-boring control volumes in the reconstruction stencil and geometric weights wij specify the relative importance of contribution from the nearby control volumes based on their distance from the reference point~xi:

wij= 1 ~x_j− ~x_i (36) This weighting is based on Mavriplis’ findings [18]. In his paper, Gooch also provides a way to incor-porate boundary conditions in the least-squares system

(11)

of 35. The boundary conditions are enforced at Gaussian quadrature points. A Dirichlet boundary condition u(~x) = f1(~x), for example, is enforced at each integration pointx~gon the boundary:

∂n = f2(~x),

the enforcing equation becomes: f2( ~xg) = ∂uR_i (~xg) ∂n = 5uR_i (~xg) ·ˆn =nx ∂u ∂x|i+ ∂2u ∂x2|i xg−xi +. . . +ny ∂u ∂y|i + ∂2u ∂y2|i xg−xi +. . . +nz ∂u ∂z|i+ ∂2u ∂z2|i xg−xi +. . . (38)

Equations 37 and 38 can then be used as additional con-straints to the least-squares system of 35.

ii.3 Boundary Conditions in Flucs

In Flucs, the concept of one layer of ’dummy cells’ added to the computational domain for dealing with boundary conditions is utilized. This is done to simplify the computation of fluxes, gradients, etc., as most quantities can now be computed in a single loop. Treatment of some boundary conditions is documented here.

Solid Wall

For inviscid flow, slip condition is applied at the boundary. That means, the velocity is tangential to the surface and for a velocity vector~v and a surface with normal~n,

~v· ~n=0 (39)

Hence, the velocity~ve in figure 6 would become

~ve= ~vi−2V2~n (40) where, V2=uinx+viny+winz. For a dummy cell, such as the one in figure 6, the pressure is approximated as simply

pe= pi (41)

Figure 6: Boundary treatment using dummy cells (marked with dotted

line)

The scalar quantities such as ρe and Ee are set equal to those in the inner cell at the boundary.

In case of viscous flow, the no-slip condition is im-posed and the fluid is assumed to have zero velocity relative to the surface. The pressure can once again be approximated using equation 41. Similarly,

ρe =ρi (42)

and

Ee =Ei (43)

The velocity components are simply reversed.~ve = −~vi. Inflow/Outflow

For a provided state variable ub at boundary, the state in the dummy cell can be calculated using linear interpolation

Ub =

Ui+Ue

2 (44)

It is also ensured that the average of the reconstructed states on left and right side of the face, UL and UR, respectively, is compatible with the boundary value Ub. Symmetry

There is no flux across a symmetry boundary, meaning that the velocity normal to the boundary is zero. Further-more, the gradient normal to the boundary of a scalar quantity and of a tangential velocity and gradient along the boundary of normal velocity are also all zero:

~n· ~5U=0

~n· ~5 ~v·~t

=0

~_t_{· ~}_{5 (~}_v_{· ~}_n_{) =}₀

(45)

where,~_{t is a vector tangential to the symmetry boundary.}

Scalar quantities like density and pressure in the dummy cell are set equal to those in the interior cell (equations 41 - 43)

(12)

Figure 7: Two cellsΩαandΩβ.

iii.

Efficient Parallelization Algorithm

As mentioned earlier, reconstruction of state variables at control volume faces requires a reconstruction stencil, whose size depends on the accuracy desired - figure 5. For k=1, or second order accuracy, only the cells in the Von Neumann neighborhood are needed, i.e. the cells that share a face with the cell in question. However, for an increasing k, the stencil grows bigger and the code must now somehow make sure that for each cell, the data of not only its neighbors, but also its neighbors’ neighbors (and so on..) is available. This would come at an immense computational cost in terms of communication on ditributed memory parallel systems. There is also a penalty in memory requirement because the least-squares system to be solved (equation 34) is bigger, thereby requiring more storage. To overcome this problem, an efficient k-exact reconstruction algorithm suggested by F. Haider et al is also implemented and tested in Flucs [10].

iii.1 Notation

Before the algorithm can be presented, first the notation used by Haider et al is described here. It is more succinct than the one used so far in this thesis and it would not be worth reformulating it. As shown in figure 7,Ωαand

Ωβare two adjacent control volumes. Equation 24 can be

written as: u(x) =u(xα) + _∑k j=1 1 j! d ∑ i1 ij=1 ∂u ∂xi1· · ·∂xij xα xi₁−xα;i1 · · ·xi_j−xα;ij (46) Next, a k-exact reconstruction operator is defined. For a space Sm(Rd)of symmetric tensors of rank m and a grid neighborhoodWαof cellΩα and k∈N with m≤k, a

k-exact mthderivative onWαis a linear map w

(m k) α :RN→ Sm(Rd)given by:

∑

β∈Wα w_αβ(m k)_;i 1 imu¯β= ∂u ∂xi1· · ·∂xim xα (47) The cell averages uγon a cellΩγ of the function u(x)in

equation 46 can be expressed as:

uγ=u(xα) + k

∑

j=1 1 j! d

∑

i1 ij=1 z(j)_αγ_;i 1 ij ∂u ∂xi1· · ·∂xij _x α (48) Here, z(j)_αγ_;i 1 ij , 1 |Ωγ| Z Ωγ xi1−xα;i1 · · ·xij −xα;ij dV (49) are the moments and|Ωγ|is the volume of the element γ.

Using equation 48, the reconstruction error of a k-exact mthderivative w(m k)α on a grid neighborhood Wα

at cellΩαapplied to a(k+1)-exact function u(x)can be

expressed as [12]: ∑ β∈Wα w(m k)_αβ_;i 1 imu¯β− ∂u ∂xi1· · ·∂xim xα = 1 (k+1)!β∈∑Wα w(m k)_αβ_;i 1 im d ∑ i₁ i_j=1 z(k+1)_αβ_;i 1 ij ∂u ∂xi1· · ·∂xij xα (50) Furthermore, a way of calculating (k+1)-exact(k+1) -order derivatives at cellΩαfrom(k)-exact(k)-order

deriva-tives in cellsΩβfor all β is needed. An operator=

(k+1) Wα is introduced: =(k+1)_W α : S (k+1)₍_Rd_{) → (}_S(k)₍_Rd₎₎m ₍₅₁₎ where, m is the number of neighbors inWα labeled β∈

{β1,· · ·, βm}. The βth component of a symmetric tensor b∈S(k+1)is defined as: =(k+1)_W α (b)β;i1 ik , d ∑ j=1 hαβ;jbi₁ i_kj + 1 (k+1)! ∑ γ∈W(_βk) w(k k)_βγ ∑d i₁ i_k+1=1 z(k+1)_βγ_;i 1 ik+1 bi1 ik+1 − 1 (k+1)! ∑ γ∈W(αk) w(k k)αγ d ∑ i1 ik+1=1 z(k+1)_αγ_;i 1 ik+1 bi1 ik+1 (52) Given a family of k-exact kthderivatives at cellΩβ: w

(k k)

β

with respective reconstructions stencilsW(k)_β and a func-tion u that is(k+1)-exact, the following identity holds:

=(k+1)_W α D(k+1) u|_x α =      w(k k)_β 1 [u] −w (k k) α [u] .. . w(k k)_β_m [u] −w(k k)α [u]      (53)

(13)

The construction of LHS of equation 53 depends on the vector of derivatives. If, say for second order deriva-tives, the symmetry is utilized then D(2) _u_|

xα consists of 6 unique entries:                        ∂2u ∂x2 ∂2u ∂x∂y ∂2u ∂x∂z ∂2u ∂y2 ∂2u ∂y∂z ∂2u ∂z2                        xα (54)

Then, the LHS for solving second order derivative can be constructed using equation 52 [11].

                          =(2)_W αβ1;xx = (2) Wαβ1;xy = (2) Wαβ1;xz . . . 0 0 =(2)_W αβ1;yx 0 . . . 0 0 0 =(2)_W αβ1;zx . . . = (2) Wαβ1;zz .. . . .. =(2)_W αβm;xx = (2) Wαβm;xy = (2) Wαβm;xz . . . 0 0 =(2)_W αβm;yx 0 . . . 0 0 0 =_W(2) αβm;zx . . . = (2) Wαβm;zz                           (55) If there exists a left inverse D_W_α of =(k+1)_W

α , then the

family of(k+1)-exact(k+1)th derivatives at cellΩα is

given by w(k+1 k+1)α [u],D (k+1) Wα      w(k k)_β 1 [u] −w (k k) α [u] .. . w(k k)_β_m [u] −wα(k k)[u]      (56)

Now, the algorithm for k-exact reconstruction using a compact stencil, i.e. only the cells in the immediate neigh-borhood, can be presented.

iii.2 Algorithm for Compact k-Exact Reconstruction

The algorithm described by Haider et al consists of three main steps:

1. Calculate w(1 1)_β [u]from the cell averages of the imme-diate neighbors (Von Neumann neighborhoodV(1)α )

for each cellΩαusing least-squares.

Figure 8: Flucs as a FSDM plugin [16].

2. Iterate from m=1 to m=k−1: compute in each cell Ωα, w (m+1 m+1) α [u]from w (1 1) β [u]for β∈V (1) α using equation 56.

3. Use equation 50 recursively for 1 ≤ m ≤ k−1 to calculate w(m k)α [u]from w

(m+1 m+1)

α [u].

Since, the notation used in this section is not always in-tuitive enough for understanding or implementing, the algorithm is written out in the simplest, expanded form for example cases of k=2 and k=3 in appendix VI and VII, respectively.

III. Implementation in Flucs

i.

Structure of Flucs

Flucs is written in C++11 to allow for high level of ab-straction while maintaining the code’s efficiency. Use of templates to write the code helps minimize the run-time overhead of these abstractions as the needed parameters can be passed on as inputs to a template [16]. The abstrac-tions also help utilize synergies, for example - between the two spatial discretization schemes: FV and DG. As seen from figure 8, most of the C++ framework needed for Flucs is provided to it via a base layer. This namespace is called Flis.

ii.

Example of Linear Reconstruction in Flucs

For FV linear reconstruction, for example, the inheritance structure in Flucs is shown in 9. The code is abstracted to the last possible level to avoid copying parts as much as possible. Linear reconstruction and constant recon-struction, for example, both use the same finite volume interface, which itself is inherited by a reconstruction interface common to DG and FV schemes.

(14)

Figure 9: Inheritance chain of linear reconstruction in Flucs.

Figure 10: Structure of Flucs.

iii.

Specifics of k-Exact Reconstruction

The already existing infrastructure for linear reconstruc-tion was used as a guideline to decide where to insert the relevant header files for k-exact reconstruction. From fig-ure 10, it is clear where the relevant k-exact reconstruction files are created. With this brief orientation to Flucs, now the specifics can be discussed.

iii.1 Integration of Moments

Recalling equation 49 for calculation of moments, it is implemented in Flucs as follows. A cell-centered Finite Volume mesh consists of just one point per cell where all the state variables are stored. That makes it difficult to calculate the moments inside a cell, because the required integration can not be performed with just one point. Hence, the code was modified to allow the Finite Volume algorithm to use a Discontinuous Galerkin mesh, which divides each cell into smaller volumes, each with its own center. The number of divisions depends on the degree of integration. This mesh with each cell (figure 11) consisting of multiple points can be used to calculate

Figure 11: Calculation of Moments

moments as follows: z(j)_αβ_;i 1 ij , 1 |Ωβ| R Ωβ xi1−xα;i1 · · ·xij−xα;ij dV = n ∑ p=1 n xi1;p−xα;i1 · · ·xij;p−xα;ij Vp o n ∑ p=1 Vp (57) An easier way from an implementation point of view is to calculate the moments of each cell zα, zβ separately, first

and then derive the necessary relations such as zαβusing

the following set of identities from [11]:

z(0)α =1, z (0) αβ =1, z(1)α =0 z (1) αβ =hαβ z(2)_αβ =h2_αβ+z(2)_β z(3)_αβ =h3_αβ+3hαβz (2) β +z (3) β (58)

iii.2 High-Order Flux Integration

Using the Discontinuous Galerkin mesh for Finite Volume calculations also enabled the use of higher-order integration. The idea is to use an appropriate number of quadrature points per cell edge to integrate flux to an appropriate degree.

Furthermore, while one normal is sufficient for a linear cell edge, at least two are required for a quadratic cell edge - such as in figure 12. Thus, higher-order mesh elements also require more than one integration point per edge.

iii.3 Derivatives Calculation: Haider

The code for calculation of derivatives for k > 1 using the Haider algorithm outlined in section II.iii works in Flucs in much the same way as the k = 1, i.e. linear

(15)

Figure 12: Higher-order integration for higher-order mesh elements

Figure 13: Implementation of Haider derivatives calculation algorithm

in Flucs

reconstruction explained in section III.ii works.

However, there is one key difference. The deriva-tives calculation for Haider algorithm receives linear reconstruction instead of constant reconstruction in case of k = 2. This is because, when at the boundaries, the algorithm for k=2 derivatives calculation, requires the first order derivatives of the state on the RHS of the least-squares system, instead of the constant values. In Flucs, these derivatives can be accessed via a member function of the ReconstructionInterface, GetStateAndGra-dients(), provided that the reconstruction is at least of linear order. For higher-order derivatives of degree m, each derivatives calculation must then receive an m−1 order reconstruction, as shown in figure 13.

iii.4 Derivatives Calculation: Gooch

The simpler derivatives calculation algorithm as explained earlier in section II was also implemtented in Flucs ii.2 to serve as a comparison to the Haider algorithm. Here, the idea is to solve the full least-squares system shown in equations 34-35, with the constraint of never accessing beyond the immediate neighboring cells.

This is done in a way analogous to the Haider al-gorithm. Consider the reconstruction stencil shown in figure 14 for calculating derivatives for k = 2 in cell α. First, a loop is carried out for each cell, where the information of other cells in their von Neumann neighborhood is stored, as shown in figure 15. During this step, the following information is stored in a matrix,

Figure 14: 2D representation of a reconstruction stencil at cell α

Figure 15: Example of the first loop at cell β1∈W(k)α for construction

of LHS0[β1]and RHS0[β1] say LHS0:            0 0 0 z(2)_β 1;xx . . . z (2) β1;zz h_β₁_γ₁;x hβ1γ1;y hβ1γ1;z z (2) γ1;xx · · · z (2) γ1;zz .. . h_β₁_γ₁;x hβ1γN;y hβ1γN;z z (2) γN;xx · · · z (2) γN;zz            (59) and RHS0: RHS0[β1] =        uβ1 uγ1 .. . uγN        (60)

Similarly, LHS0 and RHS0 can be calculated at cells β2,

β3and β4. Second, a loop is carried out through all cells once again, only this time, building LHS and RHS using LHS0[γ∈W(k)α ]and RHS0[γ∈W

(k)

α ], as shown in figure

16.

(16)

Figure 16: The Second loop at cell α for construction of LHS[α]and RHS[α]                   1 0+hαβ1;x 0+hαβ1;y . . . z (2) αβ1;zz 1 hαβ1;x+hβ1γ1;x hαβ1;y+hβ1γ1;y · · · z (2) αγ1;zz .. . 1 0+h_αβ_N;x 0+hαβN;y . . . z (2) αβN;zz .. . 1 h_αβ_N;x+hβNγ1;x hαβN;y+hβNγN;y · · · z (2) αγN;zz                   (61) and RHS[α]by: RHS[α] =                   uβ1 uγ1 .. . uβN .. . uγN                   (62)

It is worth noting that identities presented in equation 58 are used to calculate z(k)αγN from z

(k)

γN using h

(k)

αγN.

The advantage of implementing the algorithm in this way, where, at all times, only the information at the immediate information is accessed, is that it could offer similar communication cost advantages as in the Haider Algorithm.

Furthermore, unlike the Haider algorithm which requires the derivatives of order k−1 for calculating derivatives of order m, given m≤k, this algorithm only requires the state variables. Hence, instead of calculating a reconstruction of order k−1 and then passing it to the derivatives calculation, one can simply pass on the

Figure 17: Implementation of Gooch derivatives calculation algorithm

in Flucs

constant (k=0) construction to the derivatives calculation for any degree of k≥1, as shown in figure 17. There was enough time to implement it for k=2 and k=3 instead of just k=2 using Haider.

It must, however, be mentioned that the memory requirements of the Gooch algorithm are expected to be greater than those of Haider, because the RHS for each cell is now bigger. The memory savings of Haider could potentially be offset by the fact that since, say for k=2, a linear reconstruction must first be carried out and then passed on to the derivatives calculation. Nevertheless, the exact computational costs and relative advantages and disadvantages of each approach are ultimately subject to future investigation.

iii.5 Boundary Treatment

As documented in equations Equations 37 and 38, Gooch enforces a boundary condition within a cell at each flux integration point, when calculating the derivatives from the least-squares system. That makes the implementation in Flucs a bit complicated, as on the boundary, there is one layer of dummy cells and each dummy cell only has one point at which the state is stored. In the interest of time, the derivatives calculation is attempted with just one point in the dummy cell, which contributes one extra equation to the least-squares system.

Hence, in each loop over the mesh elements, only one dummy stencil is considered at the boundary, as shown in figures 18 and 19. Furthermore, the moments of these dummy cells are assumed to be the same as that of the inner cell they share the face with. In figure 19, for example, that means

z(k)_β₄ ,z(k)α z(k)γ8 ,z (k) β1 z(k)γ6 ,z (k) β3 (63)

(17)

Figure 18: Dummy cells at boundary cells in loop 1

Figure 19: Dummy cells at boundary cells in loop 2

IV. Results

Within the stipulated time for this thesis, this is what could be implemented in Flucs:

1. k=2 FV reconstruction using Haider algorithm for derivatives calculations

2. k=2 FV reconstruction using Gooch algorithm for derivatives calculations

3. k=3 FV reconstruction using Gooch algorithm for derivatives calculations

i.

Function Reconstruction

In order to debug the code and see if the calculated gradi-ents were correct, the computational domain was initial-ized with a polynomial of degree k such as

u(x) = u(xα) +ax+by+cz+dx2+2exy

+2 f xz+gy2+2hyz+iz2 (64) or,

u(x) = u(xα) +ax+by+cz+dx2+2exy+2 f xz

+gy2₊_2hyz₊_iz2₊_jx3₊_3kx2_y₊_3lx2_z₊_3mxy2

+3nxyz+3oxz2+py3+3qy2z+3ryz2+sz3 (65)

and so on and so forth. Then, the calculated gradients should ideally be equal to the coefficients of the polyno-mial up to machine accuracy.

It was observed that for k = 2, using both Haider and Gooch algorithms, k-exactness was observed in the inner cells. However, in cells near the boundary (up to two layers), the k-exact gradients could not be seen. This was also the case for k=3 using Gooch algorithm, where k-exactness is observed for inner cells up to 3 layers away from the boundary. The implications of this will become clear in the following subsection.

ii.

Euler Cases

There was only enough time to implement the code for reconstruction of convective fluxes and hence only Euler calculations were carried out (in 2D).

ii.1 NACA0012

A grid convergence study on the NACA0012 was carried out for a series of 6 meshes, with the element edge in a mesh halving with each level. The coarsest and finest meshes are shown in figure 20. The idea was to see if the error drops with(k+1)thorder for a reconstruction scheme of order k. The calculations were carried out for a Mach number of 0.5 and angle of attack of 2◦ using implicit Backward Euler temporal discretization scheme. It was ensured that the residuals dropped to machine accuracy.

The reference value of CL used for calculating the error is based on the highest CL value that was obtained for all the schemes and algorithms: CLre f = 0.2371,

whereas CDre f = 0.0 is used. First, the error in the

lift coefficient and drag coefficient for already existing schemes are plotted for comparison in figures 21 and 22.

It is seen that all schemes show second order ac-curacy. DG k = 2 scheme should theoretically show third-order accuracy; however, that would require the mesh to have higher order elements at solid walls in order to capture the geometry’s curvature exactly. The DG scheme’s order will be demonstrated in the next test case. The CL curves of all the schemes also show similar traits. It is also worth noting that least-squares (LS) and Green-Gauss (GG) gradient calculations both give very similar results. The slight discrepancy in the absolute error, or the vertical position of the LS curve can be matched to that of GG, if the inverse distance-based weighting of equation 36 is used.

(18)

Figure 20: Coarsest and finest mesh used in NACA0012 error conver-gence study. Cells(1/2) 101 ₁₀2 Error(C L ) 10-4 10-3 10-2 10-1 NACA0012

k=1 Green-Gauss, Integration Degree=1

k=1 least-squares, Integration Degree=1, 0 Weighting DG k=1

DG k=2

Figure 21: Error in CLfor NACA0012 using already existing schemes

Cells(1/2) 101 ₁₀2 Error(C D ) 10-4 10-3 10-2 10-1 100 NACA0012

k=1 Green-Gauss, Integration Degree=1

k=1 least-squares, Integration Degree=1, 0 Weighting DG k=1

DG k=2 2nd order 3rd order

Figure 22: Error in CDfor NACA0012 using already existing schemes

Now, the schemes implemented using Gooch in this thesis are plotted against k = 1 LS reconstruction for comparison. The Haider algorithm for one reason or another does not converge in this test case. Before going any further, it must be immediately noted that the point used for level 3 grid (the middle point) for k=3 is a fictional interpolatory data point. In practice, the code output negative drag at level 3 mesh. Furthermore, for level 5 and level 6 meshes (not plotted for k=3), the drag increases compared to level 4.

For k = 2, one can see a slight improvement in the order of the error between level 2 and level 3. However, with finer meshes, the error approaches the same value as that of k = 1. These results seem to be a strong indication that the earlier ’bandage solution’ attempted to deal with the calculation of gradients at boundary cells and its failure to return k-exact derivatives in cells near boundaries is detrimental to the solution.

ii.2 Smooth Bump

The next case tried was inviscid flow through a channel with a smooth bump [7]. The entropy in the flow field is expected to be constant and hence, the L2 norm of the error is computed as follows:

Error(Entropy) = v u u u t R Ω p p∞ · ρ_∞ ρ γ −1 2 dV V (66)

The coarsest and finest meshes used in this study are shown in figure 25. Furthermore, there are two variants of

(19)

Cells(1/2) 101 ₁₀2 Error(C L ) 10-4 10-3 10-2 10-1 NACA0012

k=1 LS, Integration Degree=2, 0 Weighting k=2 Gooch, Integration Degree=2 k=3 Gooch, Integration Degree=2

Figure 23: Error in CLfor NACA0012 using high order FV schemes

Cells(1/2) 101 ₁₀2 Error(C D ) 10-4 10-3 10-2 10-1 NACA0012

k=1 LS, Integration Degree=2, 0 Weighting k=2 Gooch, Integration Degree=2 k=3 Gooch, Integration Degree=2 2nd order

3rd order 4th order

Figure 24: Error in CDfor NACA0012 using high order FV schemes

Figure 25: Coarsest and finest Q1 meshes used in the smooth bump

study.

Figure 26: Upper: Q1 mesh with linear edges, lower: Q2 mesh with

(20)

structured meshes used. The Q1 mesh used linear edges for all cell elements, whereas Q2 meshes have quadratic edges, that would help them conform to the surface geometry of the bump better, as illustrated in figure 26. As before, the error in entropy using already exist-ing spatial discretization schemes is plotted in figure 27. For the DG scheme, it is observed that using Q2 meshes gives a lower error. For a polynomial of order 2 for example, one can even see the third order when Q2 meshes are used. The FV algorithm on the contrary, shows a slight deterioration in the results when a higher integration degree is used on a mesh with higher order elements. This is presumably a result of improper reconstruction at the boundary. In other words, one probably needs to modify the code to take the curvature of the elements into account, when adding equations to the least-squares system using dummy cells.

Then, the higher order FV schemes are plotted in figure 28. Not so surprisingly, all the higher-order schemes fail to demonstrate higher than second order. The absolute position of the curves changes vertically with respect to the k=1 curve, but the slope of the error shows no promise. This is another indication that the reconstruction at boundaries was not properly handled and that one would need to introduce improvements to the current implementation in Flucs, before the expected results can be obtained.

ii.3 Ehrenfried Vortex

The improper reconstruction at boundaries motivated for a test case, where the boundaries had minimal influence on the flow. While it might not be of physical interest in some cases, it would still be numerically relevant to see if the order of error reduction for higher order FV schemes can be observed.

The field is initialized with an Ehrenfried vortex super-posed on to a parallel flow [6]:

ρ ρ_∞ =1− 4ρv 2ρ∞ " 1+sin π 1− r R0 3 −π 2 !# , r≤R0 (67) where,4ρvis the prescribed density difference between the center of the vortex and the flow outside. Index∞ indicates reference value. R0 is the radius of the vortex. With the assumption of constant isentropy throughout the vortex, pressure can be given by

p p0 = ρ ρ0 γ (68) Cells(1/2) 101 ₁₀2 Error(Entropy) 10-4 10-3 10-2 Bump k=1 LS, Integration Degree=1, Q1 k=1 LS, Integration Degree=2, Q2 DG k=1 Q1 DG k=1 Q2 DG k=2 Q1 DG k=2 Q2 2nd Order 3rd Order

Figure 27: Error in Entropy for a smooth bump using already existing

schemes Cells(1/2) 101 ₁₀2 Error(Entropy) 10-4 10-3 Bump k=1 LS, Q2 k=2, Haider Q2 k=2, Gooch Q2 k=3, Gooch Q1 k=3, Gooch Q2 2nd Order 3rd Order 4th Order

Figure 28: Error in entropy for a smooth bump using high order FV

schemes; Integration degree used for all calculations here is equal to 2

(21)

Figure 29: Flow field initialized with the Ehrenfried vortex.

The tangential velocity is calculated by applying momen-tum equation in radial direction:

dp dr = ρv2_θ r (69) or, v_θ= s rdp_dr ρ (70) where, dp dr = ρ_∞γ ρ ρ_∞ γ−1₄ ρv 2ρ∞sin π 1− r R0 3! 3π R0 1− r R0 2 (71) After initializing the field with this vortex (figure 29), the vortex is transported through the domain for a given distance using explicit, 4thorder Runge-Kutta time stepping scheme. The obtained flow field can then be compared with the analytical solution at the final position of the vortex. It was ensured that at all times, the time step was small enough to not affect the solution. First, the test was attempted with the bump mesh from the previous case. However, only level 4 and level 5 meshes were fine enough for the vortex to be far enough from the boundary by more than 3 layers of cells.

k = 1 and k = 3 FV schemes were tested on the two Q1 meshes and the L2-norm of the error in average cell density of the analytical and the computed flow field is plotted in figure 30. There is a clear difference in the slope of the two schemes. However, it is also observed that k=1 does slightly better than expected while, k=3 a bit worse.

This further motivated some simulations on finer meshes. However, this time, unstructured meshes were used. The first variant consists of triangles and the second one consists of quads, as shown in figure 31. Both the meshes have some deliberately skewed cells that look

Cells(1/2) 60 70 80 90 100 110 120 130 140 150 Error( ; ) 10-5 10-4

Ehrenfried Vortex, Structured Mesh

k=1 LS, Integration Degree=2 k=3 Gooch, Integration Degree=2 2nd order

4th order

Figure 30: Error in density of k=3 FV scheme compared to k=1 LS scheme.

Figure 31: Coarsest unstructured meshes used for calculating error in

(22)

Cells(1/2) 102 103 Error( ; ) 10-6 10-5 10-4

Ehrenfried Vortex, Unstructured Mesh, Triangles

k=1 LS k=2 Gooch k=3 Gooch 2nd order 3rd order 4th order

Figure 32: Error in density of k-exact FV schemes with an

integra-tion degree of 2 compared to k = 1 LS scheme on the unstructured mesh with triangular elements.

like waves. The results for both variants of the mesh using Gooch derivatives calculation and an integration degree of 2 for moments and flux are plotted in figures 32 and 33. In case of the mesh made up of triangles, it is observed that the linear reconstruction scheme shows second order behavior. However, the order is not observed for k =2, even if the absolute error is smaller than that for k = 1. For k = 3, the scheme seems to approach 3rd order accuracy, instead of 4th. A possible reason for these results could be that the skewed cells in the mesh throw the least-squares calculations off. A distance-based weighting function could help with that. Nevertheless, it is interesting to see a difference in the absolute error of the different schemes. The error of k=3, for example, is one order of magnitude lower than that of k=1, even on the coarsest mesh.

For the mesh with quads, the order seems to get better as the mesh is refined. Both k=1 and k=2 seem to have the same slope - second order between level 1 and 2 and higher than second order between level 2 and 3. The reason for this behavior is not clear at the moment. With the Haider algorithm for k= 2, it can be seen that the curve approaches third order eventually. k=3 Gooch algorithm seems to approach 4thorder with finer meshes.

Cells(1/2) 100 150 200 250 300 350 400 450 500 Error( ; ) 10-6 10-5 10-4

Ehrenfried Vortex, Unstructured Mesh, Quads

k=1 LS k=2 Gooch k=2 Haider k=3 Gooch 2nd order 3rd order 4th order

Figure 33: Error in density of k-exact FV schemes with an

integra-tion degree of 2 compared to k = 1 LS scheme on the unstructured mesh with quad elements.

V. Conclusion

The scope of the topic was narrowed throughout the course of the thesis to produce meaningful results within the stipulated time of five to six months. This thesis extends the order of the reconstruction of state for convective fluxes in FV algorithm from constant and linear to quadratic and cubic. Two approaches for calculating derivatives were implemented in Flucs and some test cases were tried.

Minimal infrastructural changes had to be made in Flucs for k-exact reconstruction. To allow for integration of moments within each cell and a higher-order integra-tion of fluxes, the mesh used by DG was fed to the FV algorithm. The level of abstraction in Flucs made it a straightforward task.

Insufficient geometric treatment of the boundary cells and the dummy cells is believed to be detrimental to the NACA0012 case and the smooth bump case. In the smooth bump case, the FV algorithms failed to show higher than second order because of this reason. Further-more, the fact that, for example the linear reconstruction deteriorates a little when meshes with quadratic edges are used instead of meshes with linear edges, is an indication that the algorithm currently does not properly handle boundary curvature.

(23)

boundaries, a test case of Ehrenfried Vortex was used. Grid convergence studies were carried out on three to four different meshes and it was observed that different reconstruction schemes showed different gradients of the slope of error. For at least structured meshes and unstructured mesh with quads, the schemes of order k approached k+1 order accuracy on sufficiently fine meshes. Finer meshes might be required for more definitive results in case of triangular elements.

With the current progress in mind, possible future extensions can now be presented along with some helpful references.

i.

Future Work

1. RANS-negSA cases: With some minor modifications, higher-order integration of the viscous fluxes can be enabled in Flucs. Furthermore, the reconstructed state and k-exact first order gradients w(1 k)α [u] can

be forwarded to the viscous flux and source term calculation. With appropriate boundary-resolved meshes, the k-exact algorithm can be tested [2] [13]. 2. Weighting for LS system: It is foreseeable that

with highly stretched cells, such as the ones used for boundary layer resolution, the calculation of derivatives might suffer from accuracy issues. An appropriate distance-based weighting can be introduced for k≥2. [18] [19]

3. Arbitrary k: If deemed worth the time and effort, the reconstruction can be programmed for an arbitrary k. Currently, the polynomial order is hard coded. 4. Appropriate boundary treatment: The attempted

bandage solution at boundaries suffers from some possible issues. A bigger stencil might be needed for k-exactness in boundary cells [19].

Gooch also recommends "over-integration" for higher-order mesh elements.

Finally, the use of dummy cells in Flucs will need to be made compatible with a way of introducing as many equations into the least-squares system as there are gaussian quadrature points in a given cell. The calculation of moments of these dummy cells remains an open question for now.

5. Decouple moment and flux integration order: Right now, the order of integration for calculation of moments and flux is the same. It can be decoupled in the future for a closer investigation of their exclusive effect on the results.

6. Code optimization: The code can be further optimized with a bit of time investment. There are some redundancies that could not be addressed due to lack of time.

For example, all the derivatives calculation algorithms require geometric information from the mesh such as moments. These calculations could be carried out at the beginning and then made accessible to different order reconstructions and derivative calculations instead of calculating them again.

For the Gooch variant implemented in Flucs, the LS system built in 61 will have repetitions of certain cells and hence the related equations. These extra lines can be taken out.

7. Comparison of Computational Costs: With the optimized k-exact FV reconstruction versions, one can easily compare their computational costs with each other and the DG scheme, in order to determine the order up till which it might be worth going for different flow problems.

8. Introduction of limiters: The cases in this thesis were run for subsonic flows. For transonic and supersonic flows with high gradients due to shock, some work to introdce limiters is expected [19] [13].

(24)

References

[1] Allmaras, Steven R., et al. "Seventh International Conference on Computational Fluid Dynamics (IC-CFD7)." Modifications and Clarifications for the Im-plementation of the Spalart-Allmaras Turbulence Model, 13 July 2012, www.iccfd.org/iccfd7/assets/pdf/ papers/ICCFD7-1902_paper.pdf.

[2] Antoniadis, Antonis F., et al. "Assessment of High-Order Finite Volume Methods on Unstructured Meshes for RANS Solutions of Aeronautical Con-figurations." Computers & Fluids, vol. 146, 2017, pp. 86-104., doi:10.1016/j.compfluid.2017.01.002. [3] Barth, Timothy J., and Dennis C. Jespersen. "The

De-sign and Application of Upwind Schemes on Unstruc-tured Meshes." NASA Technical Reports Server, 1 Jan. 1989, ntrs.nasa.gov/search.jsp?R=19890037939. [4] Blazek, Jiri. Computational Fluid Dynamics: Principles

and Applications. Elsevier Science, 2001.

[5] Cambier, L., et al. "An Overview of the Multi-Purpose ElsA Flow Solver." Aerospace Lab, www.aerospacelab-journal.org/sites/www. aerospacelab-journal.org/files/AL2-10.pdf. [6] Ehrenfried, K., and G.E.A. Meier. "Ein

Finite-Volumen-Verfahren Zur Berechnung Von Insta-tionären, Transsonischen Strömungen Mit Wirbeln." DLR-Beauftragter, Elib, 1995, elib.dlr.de/37305/. [7] Galbraith , Marshall, and Carl Ollivier-Gooch. "BI2

- Inviscid Flow through a Channel with a Smooth Gaussian Bump." Cenaero.

[8] Gerhold, T. "Tau-Overview." 1st-day-theory. 2008. [9] "Git." Git-Scm, git-scm.com/.

[10] Haider, F., et al. "Parallel Implementa-tion of k-Exact Finite Volume Reconstruc-tion on Unstructured Grids." 12 June 2013, www.math.univ-metz.fr/~croisil/Postscript/ paper-honom-Jun-12-2013.pdf.

[11] Haider, Florian. "Discrétisation en maillage non struc-turé général et applications LES." l’Université Pierre et Marie Curie, 2009. HAL,tel.archives-ouvertes.fr/ tel-00813512/document.

[12] Haider, Florian, et al. "Finite Volumes for Complex Applications : Problems and Perspectives." HAL Archives-Ouvertes, Efficient Implementation of High Order Reconstruction in Finite Volume Methods, 27 Jan. 2011, hal.archives-ouvertes.fr/hal-01283666.

[13] Jalali, Alireza, and Carl Ollivier-Gooch. "Higher-Order Unstructured Finite Volume RANS Solution of Turbulent Compressible Flows." Computers & Fluids, vol. 143, 2017, pp. 32-47., doi:10.1016/j.compfluid. 2016.11.004.

[14] Jägersküpper, Jens. "The next-Generation CFD Solver Flucs - HPC Aspects." T-Systems-SfR HPCN workshop. 31 May 2017, Braunschweig. http://www.t-systems-sfr.com/e/downloads/ 2015/vortraege/1jaegerskuepper.pdf

[15] Kroll, N., et al. "DLR Project Digital-X: towards Vir-tual Aircraft Design and Flight Testing Based on High-Fidelity Methods." CEAS Aeronautical Jour-nal, vol. 7, no. 1, 2015, pp. 3-27., doi:10.1007/ s13272-015-0179-7.

[16] Leicht, Tobias, et al. "Deutscher Luft- Und Raum-fahrtkongress 2016." DLR-PROJECT DIGITAL-X NEXT GENERATION CFD SOLVER ’FLUCS’ . [17] Li, Wana. "High-Order Finite Volume Method for

the Compressible Flows." Efficient Implementation of High-Order Accurate Numerical Methods on Unstruc-tured Grids, Springer Theses, 2014, pp. 15-36.

[18] Mavriplis, Dimitri. "16th AIAA Computational Fluid Dynamics Conference." Revisiting the Least-Squares Procedure for Gradient Reconstruction on Unstructured Meshes. , 26 June 2003, arc.aiaa.org/doi/10.2514/ 6.2003-3986.

[19] Ollivier-Gooch, Carl. "High-Order Unstructured Mesh Finite Volume Methods." 37th Advanced CFD Lecture Series . Recent Developments in Higher Or-der Methods and Industrial Application in Aeronau-tics, 12 Dec. 2013, Sint-Genesius-Rode, Von Karman Institute for Fluid Dynamics.

[20] "PARA 2010 : Conference on State of the Art in Sci-entific and Parallel Computing." The FlowSimulator Framework for Massively Parallel CFD Applications, 9 June 2010, elib.dlr.de/67768/1/para10-paper-44. pdf.

[21] Spiegel, Seth C., et al. "A Survey of the Isentropic Euler Vortex Problem Using High-Order Methods." NASA Technical Reports Server, 22 June 2015, ntrs. nasa.gov/search.jsp?R=20150018403.

(25)

Figure 34: Reconstruction stencil for calculating derivatives in cellΩα.

VI. Appendix A

An example of k=2 is used here to demonstrate how the algorithm for a compact k-exact reconstruction presented in section iii.2 works. Consider a cellΩαin figure 34, where first and second order derivatives need to be calculated from

just its immediate neighborsΩβ1,Ωβ2,Ωβ3. The stencil is visuaized in 2D, but the calculations are carried out in 3D.

Step1: Calculate the following 1-exact derivatives using the least-squares method as shown in section ii.2 : ∂u ∂x, ∂u ∂y, ∂u ∂z xα , ∂u ∂x, ∂u ∂y, ∂u ∂z x_β1 , ∂u ∂x, ∂u ∂y, ∂u ∂z x_β2 , ∂u ∂x, ∂u ∂y, ∂u ∂z x_β3 (72)

Step2: Calculate w(2 2)α [u]using equation 56. The system of equations to be solved is set up using equation 52 :

=(2)_W α(b)β1;i1 , 3

∑

j=1 hαβ;jbi1j+ 1 2!

∑

γ∈W(1) β1 w(1 1)_β 1γ 3

∑

i₁ i2=1 z(2)_β 1γ;i1 i2 bi1 i2− 1 2!

∑

γ∈W(α1) w(1 1)αγ d

∑

i₁ i2=1 z(2)_αγ_;i 1 i2 bi1 i2 (73) or, =_W(2) α(b)β1;x, (xβ1−xα)bxx+ (yβ1−yα)bxy+ (zβ1−zα)bxz +1 2 ∂u ∂x _x β1 1 |Ωγ| R Ωγ x−xβ1 2 dV bxx+_|_Ω1 γ| R Ωγ x−xβ1 y−yβ1 dV bxy+ 1 |Ωγ| R Ωγ x−xβ1 z−zβ1 dV bxz −1 2 ∂u ∂x _x α 1 |Ωγ| R Ωγ(x−xα) 2_{dV b} xx+ _|_Ω1 γ| R Ωγ(x−xα) (y−yα)dV bxy+ 1 |Ωγ| R Ωγ(x−xα) (z−zα)dV bxz , (74) =_W(2) α(b)β1;y, (xβ1−xα)byx+ (yβ1−yα)byy+ (zβ1−zα)byz +1₂ ∂u ∂y _x β1 1 |Ωγ| R Ωγ y−yβ1 x−xβ1 dV byx+ 1 |Ωγ| R Ωγ y−yβ1 2 dV byy+_|_Ω1 γ| R Ωγ y−yβ1 z−zβ1 dV byz −1 2 ∂u ∂y _x α 1 |Ωγ| R Ωγ(y−yα) (x−xα)dV byx+ 1 |Ωγ| R Ωγ(y−yα) 2 dV byy+_|_Ω1 γ| R Ωγ(y−yα) (z−zα)dV byz , (75) =_W(2) α(b)β1;y, (xβ1−xα)bzx+ (yβ1 −yα)bzy+ (zβ1−zα)bzz +1₂ ∂u ∂z _x β1 1 |Ωγ| R Ωγ z−zβ1 x−x_β₁ dV bzx+_|_Ω1 γ| R Ωγ z−zβ1 y−y_β₁ dV bzy+ _|_Ω1 γ| R Ωγ z−zβ1 2 dV bzz −1 2 ∂u ∂z _x α 1 |Ωγ| R Ωγ(z−zα) (x−xα)dV bzx+ 1 |Ωγ| R Ωγ(z−zα) (y−yα)dV bzy+ 1 |Ωγ| R Ωγ(z−zα) 2 dV bzz , (76)

Implementation of k-Exact Finite Volume Reconstruction in DLR’s Next-Generation CFD Solver: Flucs and its Comparison to Other Methods

,

STOCKHOLM SWEDEN 2017