MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
fvcartesiansyseqnrans.h
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#ifndef FVSYSEQNRANS_H_
8#define FVSYSEQNRANS_H_
9
10#include <algorithm>
11#include <array>
12#include <iterator>
13#include <limits>
14#include <numeric>
15#include "INCLUDE/maiaconstants.h"
16#include "INCLUDE/maiatypes.h"
17#include "MEMORY/alloc.h"
18#include "filter.h"
19#include "fvcartesiansyseqnns.h"
20#include "fvransmodelconstants.h"
21
22// Only needed (really: possible to use) in single-threaded applications
23#ifndef _OPENMP
24#include "MEMORY/scratch.h"
25#endif
26
27
28template <MInt nDim, class RANSModel>
29class FvSysEqnRANS : public FvSysEqnNS<nDim> {
30 // Declare parent a friend so that CRTP can access private methods/members
31 // friend class FvSysEqn<nDim>;
32
33 public:
34 FvSysEqnRANS(const MInt solverId, const MInt noSpecies);
35
36 template <MInt stencil = AUSM>
37 inline void Ausm(const MInt orientation, const MFloat upwindCoefficient, const MFloat A, const MFloat* const leftVars,
38 const MFloat* const rightVars, const MFloat* const NotUsed(srfcCoeff), MFloat* const flux);
39
40 inline void AusmBndryCorrection(const MInt orientation, const MFloat A, const MFloat* const leftVars,
41 const MFloat* const rightVars, MFloat* const flux);
42
43 // Dispatch function for classic five-point-style viscous flux schemes
44 template <MInt stencil>
45 inline void viscousFlux(const MInt orientation, const MFloat A, const MFloat* const vars0, const MFloat* const vars1,
46 const MFloat* const slope0, const MFloat* const slope1, const MFloat* const srfcCoeff,
47 const MFloat f0, const MFloat f1, MFloat* const flux) {
48 if(stencil == FIVE_POINT) {
49 viscousFluxFivePoint(orientation, A, vars0, vars1, slope0, slope1, srfcCoeff, f0, f1, flux);
50 } else {
51 TERMM(1, "Viscous Flux Scheme not implemented.");
52 }
53 }
54
55 // Dispatch function for compact viscous flux schemes
56 template <MInt stencil>
57 inline void viscousFlux(const MInt orientation, const MFloat A, const MBool isBndry,
58 const MFloat* const surfaceCoords, const MFloat* const coord0, const MFloat* const coord1,
59 const MFloat* const cellVars0, const MFloat* const cellVars1, const MFloat* const vars0,
60 const MFloat* const vars1, const MFloat* const slope0, const MFloat* const slope1,
61 const MFloat f0, const MFloat f1, MFloat* const flux) {
62 if(stencil == THREE_POINT) {
63 viscousFluxThreePoint(orientation, A, isBndry, surfaceCoords, coord0, coord1, cellVars0, cellVars1, vars0, vars1,
64 slope0, slope1, f0, f1, flux);
65 } else {
66 TERMM(1, "Viscous Flux Scheme not implemented.");
67 }
68 }
69
70 inline void viscousFluxFivePoint(const MInt orientation, const MFloat A, const MFloat* const vars0,
71 const MFloat* const vars1, const MFloat* const slope0, const MFloat* const slope1,
72 const MFloat* const NotUsed(srfcCoeff), const MFloat f0, const MFloat f1,
73 MFloat* const flux);
74
75 inline void viscousFluxThreePoint(const MInt orientation, const MFloat A, const MBool isBndry,
76 const MFloat* const surfaceCoords, const MFloat* const coord0,
77 const MFloat* const coord1, const MFloat* const cellVars0,
78 const MFloat* const cellVars1, const MFloat* const vars0, const MFloat* const vars1,
79 const MFloat* const slope0, const MFloat* const slope1, const MFloat f0,
80 const MFloat f1, MFloat* const flux);
81
82 inline void computePrimitiveVariables(const MFloat* const cvarsCell, MFloat* const pvarsCell,
83 const MFloat* const NotUsed(avarsCell));
84
85 inline void computeConservativeVariables(const MFloat* const pvarsCell, MFloat* const cvarsCell,
86 const MFloat* const NotUsed(avarsCell));
87
88 inline std::vector<std::vector<MFloat>> conservativeSlopes(const MFloat* const pvarsCell,
89 const MFloat* const cvarsCell,
90 const MFloat* const avarsCell,
91 const MFloat* const slopesCell);
92
93 void computeVolumeForces(const MInt cellId, MFloat* vars, MFloat* rhs, const MFloat cellVol,
94 const MFloat* const slopes, const MInt recData, const MInt* const recNghbr,
95 const MFloat* const recConst, const MInt noRecNghbr, MFloat dist) {
96 IF_CONSTEXPR(m_ransModel == RANS_SA_DV) {
97 computeVolumeForcesRANS_SA(cellId, vars, rhs, cellVol, slopes, recData, recNghbr, recConst, noRecNghbr, dist);
98 }
99 IF_CONSTEXPR(m_ransModel == RANS_FS) {
100 computeVolumeForcesRANS_FS(cellId, vars, rhs, cellVol, slopes, recData, recNghbr, recConst, noRecNghbr, dist);
101 }
102 IF_CONSTEXPR(m_ransModel == RANS_KOMEGA) {
103 computeVolumeForcesRANS_KOMEGA(cellId, vars, rhs, cellVol, slopes, recData, recNghbr, recConst, noRecNghbr, dist);
104 }
105 }
106
107 void computeVolumeForcesRANS_SA(const MInt cellId, MFloat* vars, MFloat* rhs, const MFloat cellVol,
108 const MFloat* const slopes, const MInt recData, const MInt* const /* recNghbr */,
109 const MFloat* const /* recConst */, const MInt noRecNghbr, MFloat dist);
110 void computeVolumeForcesRANS_FS(const MInt cellId, MFloat* vars, MFloat* rhs, MFloat cellVol,
111 const MFloat* const slopes, MInt recData, const MInt* const recNghbr,
112 const MFloat* const recConst, MInt noRecNghbr, MFloat /* dist */);
113 void computeVolumeForcesRANS_KOMEGA(const MInt cellId, MFloat* vars, MFloat* rhs, const MFloat cellVol,
114 const MFloat* const slopes, const MInt /* recData */,
115 const MInt* const /* recNghbr */, const MFloat* const /* recConst */,
116 const MInt /* noRecNghbr */, MFloat /* dist */);
117
118 public:
119 static constexpr MInt m_noRansEquations = RANSModel::noRansEquations;
120 static constexpr MInt m_ransModel = RANSModel::ransModel;
121
122 private:
123 using Base = FvSysEqnNS<nDim>;
124 using Base::index0;
125 using Base::index1;
126 using Base::m_noSpecies;
127 using Base::m_solverId;
128
129 // Equation specific values
130 using Base::m_F1BGammaMinusOne;
131 using Base::m_F1BPr;
132 using Base::m_gamma;
133 using Base::m_gammaMinusOne;
134 using Base::m_gFGMOrPr;
135 using Base::m_Pr;
136 using Base::m_referenceTemperature;
137 using Base::m_sutherlandConstant;
138 using Base::m_sutherlandConstantThermal;
139 using Base::m_sutherlandPlusOne;
140 using Base::m_sutherlandPlusOneThermal;
141
142 using Base::getArray012;
143
144 const MFloat m_Pr_t = 0.9;
145
146 MBool m_SACurvature = false;
147
148 public:
149 using Base::m_muInfinity;
150 using Base::m_Re0;
151
152 // Hold indices for primitive and conservative variables
153 struct ConservativeVariables;
154 struct FluxVariables;
155 struct PrimitiveVariables;
156 static constexpr MBool hasSC = false;
157 ConservativeVariables* CV = nullptr;
158 FluxVariables* FV = nullptr;
159 PrimitiveVariables* PV = nullptr;
160};
161
162
165template <MInt nDim, class RANSModel>
166struct FvSysEqnRANS<nDim, RANSModel>::ConservativeVariables {
167 static constexpr MInt m_noRansEquations = RANSModel::noRansEquations;
168
169 static constexpr MInt Segfault = std::numeric_limits<MInt>::min();
170
171 static constexpr MInt RHO_U = 0;
172 static constexpr MInt RHO_V = 1;
173 static constexpr MInt RHO_W = nDim == 3 ? 2 : Segfault;
174 static constexpr std::array<MInt, nDim> RHO_VV = getArray012();
175 static constexpr MInt RHO_E = nDim;
176 static constexpr MInt RHO = nDim + 1;
177 static constexpr MInt RHO_N = nDim + 2;
178 static constexpr MInt RHO_K = nDim + 2;
179 static constexpr MInt RHO_OMEGA = nDim + 3;
180 static constexpr MInt RHO_C = nDim + 2 + m_noRansEquations;
181
182 MUint m_noSpecies;
183 const MInt noVariables;
184
185 MInt* RHO_NN = nullptr;
186 MInt* RHO_Y = nullptr;
187
188 ConservativeVariables(const MInt noSpecies);
189 ~ConservativeVariables();
190};
191
192template <MInt nDim, class RANSModel>
193struct FvSysEqnRANS<nDim, RANSModel>::FluxVariables : ConservativeVariables {
194 FluxVariables(const MInt noSpecies);
195};
196
199template <MInt nDim, class RANSModel>
200struct FvSysEqnRANS<nDim, RANSModel>::PrimitiveVariables {
201 static constexpr MInt m_noRansEquations = RANSModel::noRansEquations;
202
203 static const MInt Segfault = std::numeric_limits<MInt>::min();
204
205 static constexpr MInt U = 0;
206 static constexpr MInt V = 1;
207 static constexpr MInt W = nDim == 3 ? 2 : Segfault;
208 static constexpr std::array<MInt, nDim> VV = getArray012();
209 static constexpr MInt RHO = nDim;
210 static constexpr MInt P = nDim + 1;
211 static constexpr MInt T = nDim + 1;
212
213 static constexpr MInt N = nDim + 2;
214 static constexpr MInt K = nDim + 2;
215 static constexpr MInt OMEGA = nDim + 3;
216 static constexpr MInt C = nDim + 2 + m_noRansEquations;
217
218 MUint m_noSpecies;
219 const MInt noVariables;
220
221 static std::vector<MString> varNames;
222
223 MInt* NN = nullptr;
224 MInt* Y = nullptr;
225
226 PrimitiveVariables(const MInt noSpecies);
227 void getPrimitiveVariableNames(MString* names);
228 ~PrimitiveVariables();
229};
230
231template <MInt nDim, class RANSModel>
232template <MInt stencil>
233inline void FvSysEqnRANS<nDim, RANSModel>::Ausm(const MInt orientation, const MFloat upwindCoefficient, const MFloat A,
234 const MFloat* const leftVars, const MFloat* const rightVars,
235 const MFloat* const NotUsed(srfcCoeff), MFloat* const flux) {
236 // catch the primitive variables rho and p,
237 // compute speed of sound, and interface mach number
238 const MFloat RHOL = leftVars[PV->RHO];
239 const MFloat PL = leftVars[PV->P];
240 const MFloat AL = sqrt(m_gamma * mMax(MFloatEps, PL / mMax(MFloatEps, RHOL)));
241 const MFloat ML = leftVars[orientation] / AL;
242
243 const MFloat RHOR = rightVars[PV->RHO];
244 const MFloat PR = rightVars[PV->P];
245 const MFloat AR = sqrt(m_gamma * mMax(MFloatEps, PR / mMax(MFloatEps, RHOR)));
246 const MFloat MR = rightVars[orientation] / AR;
247
248 // calculation of the resulting pressure and mach number on the surface
249 const MFloat MLR = 0.5 * (ML + MR);
250 const MFloat PLR = PL * (0.5 + upwindCoefficient * ML) + PR * (0.5 - upwindCoefficient * MR);
251
252 // calculation of the left and right rho*a
253 const MFloat RHO_AL = RHOL * AL;
254 const MFloat RHO_AR = RHOR * AR;
255
256 // calculation of the resulting mass flux through the surface
257 const MFloat RHO_U2 = 0.25 * (MLR * (RHO_AL + RHO_AR) + fabs(MLR) * (RHO_AL - RHO_AR));
258 const MFloat AbsRHO_U2 = fabs(RHO_U2);
259
260 // calculation of the energy:
261 const MFloat PLfRHOL = PL / RHOL;
262 const MFloat PRfRHOR = PR / RHOR;
263
264 // calculation of the velocity magnitude
265 const MFloat U2L = std::inner_product(&leftVars[PV->U], &leftVars[PV->U] + nDim, &leftVars[PV->U], 0.0);
266 const MFloat U2R = std::inner_product(&rightVars[PV->U], &rightVars[PV->U] + nDim, &rightVars[PV->U], 0.0);
267
268 const MFloat e0 = PLfRHOL * m_F1BGammaMinusOne + 0.5 * U2L + PLfRHOL;
269 const MFloat e1 = PRfRHOR * m_F1BGammaMinusOne + 0.5 * U2R + PRfRHOR;
270
271 std::array<MFloat, nDim> pFactor{};
272 pFactor[orientation] = 1.0;
273
274 for(MInt n = 0; n < nDim; n++) {
275 flux[FV->RHO_VV[n]] = (RHO_U2 * (leftVars[PV->VV[n]] + rightVars[PV->VV[n]])
276 + AbsRHO_U2 * (leftVars[PV->VV[n]] - rightVars[PV->VV[n]]) + PLR * pFactor[n])
277 * A;
278 }
279
280 flux[FV->RHO_E] = (RHO_U2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1)) * A;
281 flux[FV->RHO] = 2.0 * RHO_U2 * A;
282
283 // RANS specific flux calculations
284 for(MInt r = 0; r < PV->m_noRansEquations; ++r) {
285 const MFloat NL = leftVars[PV->NN[r]];
286 const MFloat NR = rightVars[PV->NN[r]];
287 flux[FV->RHO_NN[r]] = (RHO_U2 * (NL + NR) + AbsRHO_U2 * (NL - NR)) * A;
288 }
289
290 // Flux calculation for species transport
291 // TODO labels:FV Make noSpecies constexpr
292 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
293 const MFloat YsL = leftVars[PV->Y[s]];
294 const MFloat YsR = rightVars[PV->Y[s]];
295 flux[FV->RHO_Y[s]] = (RHO_U2 * (YsL + YsR) + AbsRHO_U2 * (YsL - YsR)) * A;
296 }
297}
298
299template <MInt nDim, class RANSModel>
300inline void FvSysEqnRANS<nDim, RANSModel>::AusmBndryCorrection(const MInt orientation, const MFloat A,
301 const MFloat* const leftVars,
302 const MFloat* const rightVars, MFloat* const flux) {
303 const MFloat PL = leftVars[PV->P];
304 const MFloat PR = rightVars[PV->P];
305 const MFloat PLR = 0.5 * (PR + PL);
306
307 std::array<MFloat, nDim> pFactor{};
308 pFactor[orientation] = 1.0;
309
310 for(MInt n = 0; n < nDim; n++) {
311 flux[FV->RHO_VV[n]] = PLR * pFactor[n] * A;
312 }
313
314 flux[FV->RHO_E] = 0.0;
315 flux[FV->RHO] = 0.0;
316
317 // RANS specific flux calculations
318 for(MInt r = 0; r < PV->m_noRansEquations; ++r) {
319 flux[FV->RHO_NN[r]] = 0.0;
320 }
321
322 // Flux calculation for species transport
323 // TODO labels:FV Make noSpecies constexpr
324 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
325 flux[FV->RHO_Y[s]] = 0.0;
326 }
327}
328
329template <MInt nDim, class RANSModel>
330inline void FvSysEqnRANS<nDim, RANSModel>::viscousFluxFivePoint(const MInt orientation, const MFloat A,
331 const MFloat* const vars0, const MFloat* const vars1,
332 const MFloat* const slope0, const MFloat* const slope1,
333 const MFloat* const NotUsed(srfcCoeff), const MFloat f0,
334 const MFloat f1, MFloat* const flux) {
335 static const MFloat rhoUDth = m_muInfinity * m_F1BPr;
336 static const MFloat F1BRe0 = F1 / m_Re0;
337
338 // calculate the primitve variables on the surface u,v,rho,p
339 const std::array<MFloat, nDim> velocity = [&] {
340 std::array<MFloat, nDim> tmp{};
341 for(MUint n = 0; n < nDim; n++) {
342 tmp[n] = F1B2 * (vars0[PV->VV[n]] + vars1[PV->VV[n]]);
343 }
344 return tmp;
345 }();
346
347 const MFloat rho = F1B2 * (vars0[PV->RHO] + vars1[PV->RHO]);
348 const MFloat Frho = F1 / rho;
349 const MFloat p = F1B2 * (vars0[PV->P] + vars1[PV->P]);
350
351 // Temperature on the surface T = gamma * p / rho
352 const MFloat T = m_gamma * p * Frho;
353
354 // Indices for the orientations
355 const MUint id0 = orientation;
356 const MUint id1 = index0[orientation];
357 const MUint id2 = index1[orientation];
358
359 // Compute A / Re
360 const MFloat dAOverRe = A * F1BRe0;
361
362 // calculate the viscosity with the sutherland law mue = T^3/2 * (1+S/T_0)(T + S/T_0)
363 const MFloat mue = (T * sqrt(T) * m_sutherlandPlusOne) / (T + m_sutherlandConstant);
364
365 // Turbulence modelling
366 MFloat mue_t = 0.0;
367
368 IF_CONSTEXPR(m_ransModel == RANS_SA_DV) {
369 const MFloat nuTilde = F1B2 * (vars0[PV->N] + vars1[PV->N]);
370 const MFloat nuLaminar = mue / rho;
371 const MFloat chi = nuTilde / (nuLaminar);
372 const MFloat fv1 = pow(chi, 3) / (pow(chi, 3) + RM_SA_DV::cv1to3);
373 const MFloat nuTurb = fv1 * nuTilde;
374 if(nuTilde > 0.0) {
375 mue_t = rho * nuTurb;
376 }
377 }
378 IF_CONSTEXPR(m_ransModel == RANS_FS) {
379 const MFloat nuTilde = F1B2 * (vars0[PV->N] + vars1[PV->N]);
380 const MFloat nuLaminar = mue / rho;
381 const MFloat chi = nuTilde / (nuLaminar);
382 const MFloat fv1 = pow(chi, 3) / (pow(chi, 3) + RM_FS::facv1to3);
383 const MFloat nuTurb = fv1 * nuTilde;
384 if(nuTilde > 0.0) {
385 mue_t = rho * nuTurb;
386 }
387 }
388 IF_CONSTEXPR(m_ransModel == RANS_KOMEGA) {
389 const MFloat k = F1B2 * (vars0[PV->K] + vars1[PV->K]);
390 const MFloat o = F1B2 * (vars0[PV->OMEGA] + vars1[PV->OMEGA]);
391 const MFloat nuTurb = k / o;
392 if(nuTurb > 0.0) {
393 mue_t = rho * nuTurb;
394 }
395 }
396
397
398 const MFloat lambda_t = mue_t / m_Pr_t;
399
400 // Wall-modeling with precomputed, approximated mue_wm
401 MFloat mue_wm = 0.0;
402 /*if(m_wmLES) {
403 if(*srfcs[srfcId].m_bndryCndId == 3399) {
404 MInt cellId = nghbrCellIds[2 * srfcId];
405 if(!a_isBndryCell(cellId)) cellId = nghbrCellIds[2 * srfcId + 1];
406 if(a_isBndryCell(cellId) && !a_isHalo(cellId)) {
407 mue_wm = computeWMViscositySpalding(cellId);
408 }
409 }
410 }*/
411
412 const MFloat lambda = (T * sqrt(T) * m_sutherlandPlusOneThermal) / (T + m_sutherlandConstantThermal);
413
414 const MUint sq0 = PV->P * nDim + id0;
415 const MUint sq1 = PV->RHO * nDim + id0;
416 const MFloat q = (lambda + lambda_t) * m_gFGMOrPr * Frho
417 * ((f0 * slope0[sq0] + f1 * slope1[sq0]) - p * Frho * (f0 * slope0[sq1] + f1 * slope1[sq1]));
418
419 // Compute the stress terms
420 const MUint s00 = id0 * nDim + id0;
421 const MUint s01 = id0 * nDim + id1;
422 const MUint s10 = id1 * nDim + id0;
423 const MUint s11 = id1 * nDim + id1;
424
425 std::array<MFloat, nDim> tau{};
426 IF_CONSTEXPR(nDim == 2) {
427 tau[id0] = (mue + mue_t + mue_wm)
428 * (f0 * (F4B3 * slope0[s00] - F2B3 * slope0[s11]) + f1 * (F4B3 * slope1[s00] - F2B3 * slope1[s11]));
429 tau[id1] = (mue + mue_t + mue_wm) * (f0 * (slope0[s01] + slope0[s10]) + f1 * (slope1[s01] + slope1[s10]));
430 }
431 else IF_CONSTEXPR(nDim == 3) {
432 const MUint s22 = id2 * nDim + id2;
433 const MUint s02 = id0 * nDim + id2;
434 const MUint s20 = id2 * nDim + id0;
435 tau[id0] = (mue + mue_t + mue_wm)
436 * (f0 * (F4B3 * slope0[s00] - F2B3 * (slope0[s11] + slope0[s22]))
437 + f1 * (F4B3 * slope1[s00] - F2B3 * (slope1[s11] + slope1[s22])));
438 tau[id1] = (mue + mue_t + mue_wm) * (f0 * (slope0[s01] + slope0[s10]) + f1 * (slope1[s01] + slope1[s10]));
439 tau[id2] = (mue + mue_t + mue_wm) * (f0 * (slope0[s02] + slope0[s20]) + f1 * (slope1[s02] + slope1[s20]));
440 }
441
442 // Compute the flux
443 for(MUint n = 0; n < nDim; n++) {
444 flux[FV->RHO_VV[n]] -= dAOverRe * tau[n];
445 }
446 flux[FV->RHO_E] -= dAOverRe * (std::inner_product(velocity.begin(), velocity.end(), tau.begin(), 0.0) + q);
447
448 IF_CONSTEXPR(m_ransModel == RANS_SA_DV) {
449 const MUint n0 = PV->N * nDim + id0;
450 const MUint rho0 = PV->RHO * nDim + id0;
451
452 const MFloat nuTilde = F1B2 * (vars0[PV->N] + vars1[PV->N]);
453 // const MFloat nuLaminar = mue/rho;
454 // source: Density Corrections for Turbulence Models," Aerospace Science and Technology, Vol. 4, 2000, pp. 1--11
455 MFloat viscflux_nu = RM_SA_DV::Fsigma
456 * (mue * (f0 * slope0[n0] + f1 * slope1[n0])
457 + sqrt(rho) * nuTilde
458 * (sqrt(rho) * (f0 * slope0[n0] + f1 * slope1[n0])
459 + nuTilde * F1 / (F2 * sqrt(rho)) * (f0 * slope0[rho0] + f1 * slope1[rho0])));
460
461 // no correction model
462 // MFloat viscflux_nu = RM_SA_DV::Fsigma * (nuLaminar + nuTilde) * (f0 * slope0[n0] + f1 * slope1[n0]);
463
464 flux[FV->RHO_N] -= dAOverRe * (viscflux_nu);
465 }
466
467 IF_CONSTEXPR(m_ransModel == RANS_FS) {
468 const MUint n0 = PV->N * nDim + id0;
469 const MFloat nuTilde = F1B2 * (vars0[PV->N] + vars1[PV->N]);
470 const MFloat nuLaminar = mue / rho;
471 MFloat viscflux_nu = rho * (nuLaminar + RM_FS::fasigma * nuTilde) * (f0 * slope0[n0] + f1 * slope1[n0]);
472 flux[FV->RHO_N] -= dAOverRe * (viscflux_nu);
473 }
474
475 IF_CONSTEXPR(m_ransModel == RANS_KOMEGA) {
476 const MUint n0 = PV->N * nDim + id0;
477 const MUint rho0 = PV->RHO * nDim + id0;
478
479 const MFloat nuTilde = F1B2 * (vars0[PV->N] + vars1[PV->N]);
480 MFloat viscflux_nu = RM_SA_DV::Fsigma
481 * (mue * (f0 * slope0[n0] + f1 * slope1[n0])
482 + sqrt(rho) * nuTilde
483 * (sqrt(rho) * (f0 * slope0[n0] + f1 * slope1[n0])
484 + nuTilde * F1 / (F2 * sqrt(rho)) * (f0 * slope0[rho0] + f1 * slope1[rho0])));
485
486 flux[FV->RHO_N] -= dAOverRe * (viscflux_nu);
487 const MUint nk = PV->K * nDim + id0;
488 const MUint no = PV->OMEGA * nDim + id0;
489
490 const MFloat k = F1B2 * (vars0[PV->K] + vars1[PV->K]);
491 const MFloat omega = F1B2 * (vars0[PV->OMEGA] + vars1[PV->OMEGA]);
492 const MFloat mueTurb = rho * k / omega;
493
494 MFloat viscflux_k = (mue + RM_KOMEGA::sigma_k * mueTurb) * (f0 * slope0[nk] + f1 * slope1[nk]);
495 MFloat viscflux_omega = (mue + RM_KOMEGA::sigma_o * mueTurb) * (f0 * slope0[no] + f1 * slope1[no]);
496
497 flux[FV->RHO_K] -= dAOverRe * (viscflux_k);
498 flux[FV->RHO_OMEGA] -= dAOverRe * (viscflux_omega);
499 }
500
501
502 // progress variable
503 const MFloat c = dAOverRe * rhoUDth;
504 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
505 const MUint i = PV->Y[s] * nDim + id0;
506 flux[FV->RHO_Y[s]] -= c * (f0 * slope0[i] + f1 * slope1[i]);
507 }
508}
509
510template <MInt nDim, class RANSModel>
511inline void FvSysEqnRANS<nDim, RANSModel>::viscousFluxThreePoint(
512 const MInt orientation, const MFloat A, const MBool isBndry, const MFloat* const surfaceCoords,
513 const MFloat* const coord0, const MFloat* const coord1, const MFloat* const cellVars0,
514 const MFloat* const cellVars1, const MFloat* const vars0, const MFloat* const vars1, const MFloat* const slope0,
515 const MFloat* const slope1, const MFloat f0, const MFloat f1, MFloat* const flux) {
516 static const MFloat rhoUDth = m_muInfinity * m_F1BPr;
517 static const MFloat F1BRe0 = F1 / m_Re0;
518
519 // calculate the primitve variables on the surface u,v,rho,p
520 const std::array<MFloat, nDim> velocity = [&] {
521 std::array<MFloat, nDim> tmp{};
522 for(MUint n = 0; n < nDim; n++) {
523 tmp[n] = F1B2 * (vars0[PV->VV[n]] + vars1[PV->VV[n]]);
524 }
525 return tmp;
526 }();
527
528 const MFloat rho = F1B2 * (vars0[PV->RHO] + vars1[PV->RHO]);
529 const MFloat Frho = F1 / rho;
530 const MFloat p = F1B2 * (vars0[PV->P] + vars1[PV->P]);
531
532 // Temperature on the surface T = gamma * p / rho
533 const MFloat T = m_gamma * p * Frho;
534
535 // Indices for the orientations
536 const MUint id0 = orientation;
537 const MUint id1 = index0[orientation];
538 const MUint id2 = index1[orientation];
539
540 // Compute A / Re
541 const MFloat dAOverRe = A * F1BRe0;
542
543 // calculate the viscosity with the sutherland law mue = T^3/2 * (1+S/T_0)(T + S/T_0)
544 const MFloat mue = (T * sqrt(T) * m_sutherlandPlusOne) / (T + m_sutherlandConstant);
545
546 const MFloat lambda = (T * sqrt(T) * m_sutherlandPlusOneThermal) / (T + m_sutherlandConstantThermal);
547
548 // TODO labels:FV @Julian, change to stack alloc
549 std::vector<MFloat> slopes(PV->noVariables * nDim);
550
551 /*#ifdef _VISCOUS_SLOPES_COMPACT_VARIANT_
552 const MFloat dx = F1B2 * (coord1[id0] - coord0[id0]);
553 #endif*/
554 for(MInt var = 0; var < PV->noVariables; var++) {
555 slopes[nDim * var + id0] =
556 cellVars1[var]
557 - cellVars0[var]
558 /*#ifdef _VISCOUS_SLOPES_COMPACT_VARIANT_
559 + (surfaceCoords[id0] - coord1[id0] + dx) * slope1[var * nDim + id0]
560 - (surfaceCoords[id0] - coord0[id0] - dx) * slope0[var * nDim + id0]
561 #endif*/
562 + (surfaceCoords[id1] - coord1[id1]) * slope1[var * nDim + id1]
563 - (surfaceCoords[id1] - coord0[id1]) * slope0[var * nDim + id1];
564 IF_CONSTEXPR(nDim == 3) {
565 slopes[nDim * var + id0] += (surfaceCoords[id2] - coord1[id2]) * slope1[var * nDim + id2]
566 - (surfaceCoords[id2] - coord0[id2]) * slope0[var * nDim + id2];
567 }
568 /*#ifdef _VISCOUS_SLOPES_COMPACT_VARIANT_
569 )
570 / (F2 * dx);
571 #else*/
572 // )
573 slopes[nDim * var + id0] /= (coord1[id0] - coord0[id0]);
574
575
576 //#endif
577 slopes[nDim * var + id1] = f0 * slope0[var * nDim + id1] + f1 * slope1[var * nDim + id1];
578 IF_CONSTEXPR(nDim == 3) {
579 slopes[nDim * var + id2] = f0 * slope0[var * nDim + id2] + f1 * slope1[var * nDim + id2];
580 }
581 if(isBndry) {
582 slopes[nDim * var + id0] = f0 * slope0[var * nDim + id0] + f1 * slope1[var * nDim + id0];
583 }
584 }
585
586 const MUint sq0 = PV->P * nDim + id0;
587 const MUint sq1 = PV->RHO * nDim + id0;
588
589 const MFloat q = lambda * m_gFGMOrPr * Frho * (slopes[sq0] - p * Frho * slopes[sq1]);
590
591 // Compute the stress terms
592 const MUint s00 = id0 * nDim + id0;
593 const MUint s01 = id0 * nDim + id1;
594 const MUint s10 = id1 * nDim + id0;
595 const MUint s11 = id1 * nDim + id1;
596
597 std::array<MFloat, nDim> tau{};
598 IF_CONSTEXPR(nDim == 2) {
599 tau[id0] = mue * (F4B3 * slopes[s00] - F2B3 * slopes[s11]);
600 tau[id1] = mue * (slopes[s01] + slopes[s10]);
601 }
602 else IF_CONSTEXPR(nDim == 3) {
603 const MUint s22 = id2 * nDim + id2;
604 const MUint s02 = id0 * nDim + id2;
605 const MUint s20 = id2 * nDim + id0;
606 tau[id0] = mue * (F4B3 * slopes[s00] - F2B3 * (slopes[s11] + slopes[s22]));
607 tau[id1] = mue * (slopes[s01] + slopes[s10]);
608 tau[id2] = mue * (slopes[s02] + slopes[s20]);
609 }
610
611 // Compute the flux
612 for(MUint n = 0; n < nDim; n++) {
613 flux[FV->RHO_VV[n]] -= dAOverRe * tau[n];
614 }
615 flux[FV->RHO_E] -= dAOverRe * (std::inner_product(velocity.begin(), velocity.end(), tau.begin(), 0.0) + q);
616
617 // progress variable
618 const MFloat c = dAOverRe * rhoUDth;
619 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
620 const MUint i = PV->Y[s] * nDim + id0;
621 flux[FV->RHO_Y[s]] -= c * (f0 * slope0[i] + f1 * slope1[i]);
622 }
623}
624
625template <MInt nDim, class RANSModel>
626inline void FvSysEqnRANS<nDim, RANSModel>::computePrimitiveVariables(const MFloat* const cvarsCell,
627 MFloat* const pvarsCell,
628 const MFloat* const NotUsed(avarsCell)) {
629 const MFloat fRho = F1 / cvarsCell[CV->RHO];
630 MFloat velPOW2 = F0;
631 for(MInt n = 0; n < nDim; ++n) { // compute velocity
632 pvarsCell[PV->VV[n]] = cvarsCell[CV->RHO_VV[n]] * fRho;
633 velPOW2 += POW2(pvarsCell[PV->VV[n]]);
634 }
635
636 // density and pressure:
637 pvarsCell[PV->RHO] = cvarsCell[CV->RHO]; // density
638 pvarsCell[PV->P] = m_gammaMinusOne * (cvarsCell[CV->RHO_E] - F1B2 * pvarsCell[PV->RHO] * velPOW2);
639
640 // RANS variables
641 for(MInt r = 0; r < PV->m_noRansEquations; r++) {
642 pvarsCell[PV->NN[r]] = cvarsCell[CV->RHO_NN[r]] * fRho;
643 }
644
645 // compute the species
646 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
647 pvarsCell[PV->Y[s]] = cvarsCell[CV->RHO_Y[s]] * fRho;
648 }
649}
650
651template <MInt nDim, class RANSModel>
652inline void FvSysEqnRANS<nDim, RANSModel>::computeConservativeVariables(const MFloat* const pvarsCell,
653 MFloat* const cvarsCell,
654 const MFloat* const NotUsed(avarsCell)) {
655 // compute rho v
656 MFloat velPOW2 = F0;
657 for(MInt vel = 0; vel < nDim; ++vel) {
658 const MFloat v = pvarsCell[vel];
659 cvarsCell[vel] = v * pvarsCell[PV->RHO];
660 velPOW2 += POW2(v);
661 }
662
663 // compute rho e
664 cvarsCell[CV->RHO_E] = pvarsCell[PV->P] * m_F1BGammaMinusOne + F1B2 * pvarsCell[PV->RHO] * velPOW2;
665
666 // copy the density
667 cvarsCell[CV->RHO] = pvarsCell[PV->RHO];
668
669 for(MInt r = 0; r < PV->m_noRansEquations; r++) {
670 cvarsCell[CV->RHO_NN[r]] = pvarsCell[PV->NN[r]] * pvarsCell[PV->RHO];
671 }
672
673 // compute rho Yi
674 for(MUint s = 0; s < PV->m_noSpecies; ++s) {
675 cvarsCell[CV->RHO_Y[s]] = pvarsCell[PV->Y[s]] * pvarsCell[PV->RHO];
676 }
677}
678
679template <MInt nDim, class RANSModel>
680inline std::vector<std::vector<MFloat>>
681FvSysEqnRANS<nDim, RANSModel>::conservativeSlopes(const MFloat* const pvarsCell, const MFloat* const NotUsed(cvarsCell),
682 const MFloat* const NotUsed(avarsCell),
683 const MFloat* const slopesCell) {
684 MFloat U2 = F0;
685 for(MInt i = 0; i < nDim; i++) {
686 U2 += POW2(pvarsCell[PV->VV[i]]);
687 }
688
689 std::vector<std::vector<MFloat>> dQ(CV->noVariables, std::vector<MFloat>(nDim));
690
691 for(MInt d = 0; d < nDim; d++) {
692 dQ[CV->RHO][d] = slopesCell[PV->RHO * nDim + d];
693 dQ[CV->RHO_E][d] = slopesCell[PV->P * nDim + d] / (m_gamma - F1) + F1B2 * U2 * slopesCell[PV->RHO * nDim + d];
694
695 for(MInt j = 0; j < nDim; j++) {
696 dQ[CV->RHO_VV[j]][d] =
697 pvarsCell[PV->VV[j]] * slopesCell[PV->RHO * nDim + d] + pvarsCell[PV->RHO] * slopesCell[PV->VV[j] * nDim + d];
698 dQ[CV->RHO_E][d] += pvarsCell[PV->RHO] * pvarsCell[PV->VV[j]] * slopesCell[PV->VV[j] * nDim + d];
699 }
700
701 for(MInt r = 0; r < m_noRansEquations; r++) {
702 dQ[CV->RHO_NN[r]][d] =
703 pvarsCell[PV->NN[r]] * slopesCell[PV->RHO * nDim + d] + pvarsCell[PV->RHO] * slopesCell[PV->NN[r] * nDim + d];
704 }
705
706 for(MUint s = 0; s < CV->m_noSpecies; ++s) {
707 dQ[CV->RHO_Y[s]][d] =
708 pvarsCell[PV->Y[s]] * slopesCell[PV->RHO * nDim + d] + pvarsCell[PV->RHO] * slopesCell[PV->Y[s] * nDim + d];
709 }
710 }
711
712 return dQ;
713}
714
715template <MInt nDim, class RANSModel>
716void FvSysEqnRANS<nDim, RANSModel>::computeVolumeForcesRANS_SA(const MInt cellId, MFloat* vars, MFloat* rhs,
717 const MFloat cellVol, const MFloat* const slopes,
718 const MInt recData, const MInt* const /* recNghbr */,
719 const MFloat* const /* recConst */,
720 const MInt noRecNghbr, MFloat dist) {
721 // return;
722
723 std::ignore = recData;
724 std::ignore = noRecNghbr;
725
726 const MFloat rRe = F1 / m_Re0;
727 const MFloat rho = vars[cellId * PV->noVariables + PV->RHO];
728 const MFloat p = vars[cellId * PV->noVariables + PV->P];
729 const MFloat nuTilde = vars[cellId * PV->noVariables + PV->N];
730 const MFloat T = m_gamma * p / rho;
731 const MFloat mue = (T * sqrt(T) * m_sutherlandPlusOne) / (T + m_sutherlandConstant);
732 const MFloat nuLaminar = mue / rho;
733
734 MFloat s = 0.0;
735 MFloat diff = 0.0;
736
737 // calculation of d(sqrt(rho) nu_tilde)/dx
738 // source: "Density Corrections for Turbulence Models", Aerospace Science and Technology, Vol. 4, 2000, pp. 1--11
739 for(MInt d = 0; d < nDim; d++) {
740 MInt indexRho = cellId * PV->noVariables * nDim + PV->RHO * nDim + d;
741 MInt indexNT = cellId * PV->noVariables * nDim + PV->N * nDim + d;
742 diff += rRe * RM_SA_DV::Fsigma * RM_SA_DV::cb2
743 * POW2((sqrt(rho) * slopes[indexNT] + nuTilde * slopes[indexRho] / (2 * sqrt(rho))));
744 // no correction model
745 // diff += rRe * RM_SA_DV::Fsigma * (rho * RM_SA_DV::cb2 * POW2(slopes[indexNT]) - (nuTilde + nuLaminar) *
746 // slopes[indexNT] * slopes[indexRho]);
747 }
748
749
750 for(MInt d1 = 0; d1 < nDim; d1++) {
751 for(MInt d2 = 0; d2 < nDim; d2++) {
752 MInt indexVV12 = cellId * PV->noVariables * nDim + PV->VV[d1] * nDim + d2;
753 MInt indexVV21 = cellId * PV->noVariables * nDim + PV->VV[d2] * nDim + d1;
754 MFloat wij = 0.5 * (slopes[indexVV12] - slopes[indexVV21]);
755 s += wij * wij;
756 }
757 }
758 s = sqrt(2 * s);
759
760 const MFloat chi = nuTilde / (nuLaminar);
761 const MFloat fv1 = pow(chi, 3) / (pow(chi, 3) + RM_SA_DV::cv1to3);
762 const MFloat Fv2 = F1 - (chi / (F1 + chi * fv1));
763 const MFloat Fdist2 = 1.0 / (dist * dist);
764 const MFloat term = nuTilde * Fdist2 * RM_SA_DV::Fkap2;
765 const MFloat stilde = s + term * Fv2 * rRe;
766 const MFloat r = mMin(10.0, rRe * term / stilde);
767 const MFloat g = r + RM_SA_DV::cw2 * (pow(r, 6) - r);
768 const MFloat Fwterm = (1 + RM_SA_DV::cw3to6) / (pow(g, 6) + RM_SA_DV::cw3to6);
769 const MFloat Fw = g * pow(Fwterm, (1.0 / 6.0));
770 const MFloat Ft2 = RM_SA_DV::Ft2; // RM_SA_DV::ct3 * exp(-RM_SA_DV::ct4*POW2(chi));
771
772 // RM_SA_DV::Ft2 = 0.0 (because simulations with no tripping term)
773 MFloat P = (F1 - Ft2) * RM_SA_DV::cb1 * /* fabs( */ nuTilde /* ) */ * mMax(stilde, 0.3 * s);
774 MFloat D = rRe * (RM_SA_DV::cw1 * Fw - RM_SA_DV::cb1 * RM_SA_DV::Fkap2 * Ft2) * pow(nuTilde, 2.0) * Fdist2;
775 MFloat prodDest = rho * (P - D);
776
777 rhs[CV->RHO_N] -= cellVol * (prodDest + diff);
778}
779
780template <MInt nDim, class RANSModel>
781void FvSysEqnRANS<nDim, RANSModel>::computeVolumeForcesRANS_FS(const MInt cellId, MFloat* vars, MFloat* rhs,
782 const MFloat cellVol, const MFloat* const slopes,
783 const MInt recData, const MInt* const recNghbr,
784 const MFloat* const recConst, const MInt noRecNghbr,
785 MFloat /* dist */) {
786 const MFloat rRe = F1 / m_Re0;
787
788 const MFloat rho = vars[cellId * PV->noVariables + PV->RHO];
789 const MFloat p = vars[cellId * PV->noVariables + PV->P];
790 const MFloat nuTilde = vars[cellId * PV->noVariables + PV->N];
791 const MFloat T = m_gamma * p / rho;
792 const MFloat mue = (T * sqrt(T) * m_sutherlandPlusOne) / (T + m_sutherlandConstant);
793 const MFloat nuLaminar = mue / rho;
794
795 MFloat SijSij = F0;
796 MFloat OijOjkSki = F0;
797 MFloat diff2 = F0;
798
799 for(MInt d1 = 0; d1 < nDim; d1++) {
800 MInt indexN1 = cellId * PV->noVariables * nDim + PV->N * nDim + d1;
801 MInt indexRHO1 = cellId * PV->noVariables * nDim + PV->RHO * nDim + d1;
802 diff2 += slopes[indexN1] * slopes[indexRHO1];
803 for(MInt d2 = 0; d2 < nDim; d2++) {
804 MInt indexVV12 = cellId * PV->noVariables * nDim + PV->VV[d1] * nDim + d2;
805 MInt indexVV21 = cellId * PV->noVariables * nDim + PV->VV[d2] * nDim + d1;
806 MFloat sij = 0.5 * (slopes[indexVV12] + slopes[indexVV21]);
807 MFloat wij = 0.5 * (slopes[indexVV12] - slopes[indexVV21]);
808 SijSij += sij * sij;
809 for(MInt d3 = 0; d3 < nDim; d3++) {
810 MInt indexVV13 = cellId * PV->noVariables * nDim + PV->VV[d1] * nDim + d3;
811 MInt indexVV31 = cellId * PV->noVariables * nDim + PV->VV[d3] * nDim + d1;
812 MInt indexVV23 = cellId * PV->noVariables * nDim + PV->VV[d2] * nDim + d3;
813 MInt indexVV32 = cellId * PV->noVariables * nDim + PV->VV[d3] * nDim + d2;
814 MFloat ski = 0.5 * (slopes[indexVV13] + slopes[indexVV31]);
815 MFloat wjk = 0.5 * (slopes[indexVV23] - slopes[indexVV32]);
816 OijOjkSki += wij * wjk * ski;
817 }
818 }
819 }
820
821 diff2 = rRe * (nuLaminar + RM_FS::fasigma * nuTilde) * diff2;
822
823 const MFloat omega = mMax(1e-16, F1 / sqrt(RM_FS::fabetcs) * sqrt(F2 * SijSij));
824 MFloat psi_o = fabs(OijOjkSki / (pow(RM_FS::fabetcs * omega, 3.0)));
825
826 MFloat term1 = 0.0;
827 MFloat term2 = 0.0;
828 MFloat domega[nDim]{};
829 MFloat dnutilde[nDim]{};
830
831 for(MInt nghbr = 0; nghbr < noRecNghbr; nghbr++) {
832 const MInt nghbrId = recNghbr[nghbr];
833
834 // calculation of omega for nghbr
835 MFloat SijSijNgbhr = 0;
836 for(MInt d1 = 0; d1 < nDim; d1++) {
837 for(MInt d2 = 0; d2 < nDim; d2++) {
838 MInt indexVV12 = nghbrId * PV->noVariables * nDim + PV->VV[d1] * nDim + d2;
839 MInt indexVV21 = nghbrId * PV->noVariables * nDim + PV->VV[d2] * nDim + d1;
840 MFloat sijNgbhr = 0.5 * (slopes[indexVV12] + slopes[indexVV21]);
841 SijSijNgbhr += sijNgbhr * sijNgbhr;
842 }
843 }
844 const MFloat omegaNgbhr = 1 / sqrt(RM_FS::fabetcs) * sqrt(2 * SijSijNgbhr);
845
846 const MFloat nutildeNgbhr = vars[nghbrId * PV->noVariables + PV->N];
847 for(MInt d = 0; d < nDim; d++) {
848 const MFloat recConst_ = recConst[nDim * (recData + nghbr) + d];
849 const MFloat deltaOmega = omegaNgbhr - omega;
850 const MFloat deltaNutilde = nutildeNgbhr - nuTilde;
851 domega[d] += recConst_ * deltaOmega;
852 dnutilde[d] += recConst_ * deltaNutilde;
853 }
854 }
855
856 for(MInt d = 0; d < nDim; ++d) {
857 term1 += domega[d] * domega[d];
858 term2 += dnutilde[d] * domega[d];
859 }
860
861 // psi_o = 0;
862 const MFloat psi_k = mMax(F0, F1 / pow(omega, F3) * (nuTilde * term1 + omega * term2));
863 const MFloat f_betc = (1 + RM_FS::fapsio1 * psi_o) / (1 + RM_FS::fapsio2 * psi_o);
864 const MFloat f_betcs = (1 + RM_FS::fapsik1 * psi_k) / (1 + RM_FS::fapsik2 * psi_k);
865 const MFloat beta = f_betc / f_betc * RM_FS::fabetc;
866 const MFloat beta_s = f_betcs / f_betcs * RM_FS::fabetcs;
867
868 MFloat Sijduidxj = F0;
869 MFloat duidxi = F0;
870 for(MInt d1 = 0; d1 < nDim; d1++) {
871 for(MInt d2 = 0; d2 < nDim; d2++) {
872 MInt indexVV12 = cellId * PV->noVariables * nDim + PV->VV[d1] * nDim + d2;
873 MInt indexVV21 = cellId * PV->noVariables * nDim + PV->VV[d2] * nDim + d1;
874 Sijduidxj += 0.5 * slopes[indexVV12] * (slopes[indexVV12] + slopes[indexVV21]);
875 }
876 MInt indexVV11 = cellId * PV->noVariables * nDim + PV->VV[d1] * nDim + d1;
877 duidxi += slopes[indexVV11];
878 }
879
880 // MFloat P_FS = mMin(0.0, 2 * (F1 - RM_FS::faalpha) * nuTilde * (Sijduidxj / omega - F1B3 * duidxi));
881 MFloat P_FS = 2 * (F1 - RM_FS::faalpha) * nuTilde * (Sijduidxj / omega - F1B3 * duidxi);
882 MFloat D_FS = (beta_s - beta) * nuTilde * omega;
883 MFloat diff1 = mMin(0.0, rRe * F2 * rho * (nuLaminar + RM_FS::fasigma * nuTilde) / omega * term2);
884
885 const MFloat prodDest_FS = rho * (P_FS - D_FS);
886 // const MFloat TU = m_turbulenceDegree;
887 // const MFloat T = m_TInfinity;
888 /* const MFloat mue_init = */
889 /* (m_TInfinity * sqrt(m_TInfinity) * m_sutherlandPlusOne) / (m_TInfinity + m_sutherlandConstant); */
890 /* const MFloat nu_init = mue_init / m_rhoInfinity; */
891 /* const MFloat k = Re * F3B2 * rho * pow(TU * m_UInfinity, F2) * RM_FS::faphi * pow(nuTilde / nu_init, F5); */
892
893 const MFloat k = 0.0;
894
895 const MFloat diff_FS = diff1 - diff2;
896
897 rhs[CV->RHO_N] -= cellVol * (prodDest_FS + diff_FS - k);
898}
899
900template <MInt nDim, class RANSModel>
901void FvSysEqnRANS<nDim, RANSModel>::computeVolumeForcesRANS_KOMEGA(const MInt cellId, MFloat* vars, MFloat* rhs,
902 const MFloat cellVol, const MFloat* const slopes,
903 const MInt /* recData */,
904 const MInt* const /* recNghbr */,
905 const MFloat* const /* recConst */,
906 const MInt /* noRecNghbr */, MFloat /* dist */) {
907 const MFloat Re = m_Re0;
908 const MFloat rRe = F1 / m_Re0;
909 const MFloat rho = vars[cellId * PV->noVariables + PV->RHO]; // vars(cellId,PV->RHO);//
910 const MFloat k = vars[cellId * PV->noVariables + PV->K]; // vars(cellId,PV->K);//
911 const MFloat omega = vars[cellId * PV->noVariables + PV->OMEGA]; // vars(cellId,PV->OMEGA);//
912
913 MFloat dukdxk = F0;
914 for(MInt d = 0; d < nDim; d++) {
915 dukdxk += slopes[PV->VV[d] * nDim + d];
916 }
917
918 MFloat omega_hat = omega; // max(omega,rRe*RM_KOMEGA::C_lim/sqrt(RM_KOMEGA::betas)*sqrt(2*Sij_Sij_));
919 // MFloat omega_hat = RM_KOMEGA::C_lim/sqrt(RM_KOMEGA::betas)*sqrt(2*SijSij));
920
921 const MFloat mue_Turb = /*m_Re0 */ rho * k / omega_hat;
922
923 MFloat tau = F0;
924 MFloat P = F0;
925 MFloat term = F0;
926 for(MInt d1 = 0; d1 < nDim; d1++) {
927 for(MInt d2 = 0; d2 < nDim; d2++) {
928 tau = rRe * mue_Turb * (slopes[PV->VV[d1] * nDim + d2] + slopes[PV->VV[d2] * nDim + d1]);
929 if(d1 == d2) {
930 tau -= (/*rRe * mue_Turb * F2B3 * dukdxk +*/ F2B3 * rho * k);
931 }
932 P += slopes[PV->VV[d1] * nDim + d2] * tau;
933 }
934 term += slopes[PV->K * nDim + d1] * slopes[PV->OMEGA * nDim + d1];
935 }
936
937 MFloat OijOjkSki = F0;
938 for(MInt d1 = 0; d1 < nDim; d1++) {
939 for(MInt d2 = 0; d2 < nDim; d2++) {
940 MFloat wij = 0.5 * (slopes[PV->VV[d1] * nDim + d2] - slopes[PV->VV[d2] * nDim + d1]);
941 for(MInt d3 = 0; d3 < nDim; d3++) {
942 MFloat ski = 0.5 * (slopes[PV->VV[d3] * nDim + d1] + slopes[PV->VV[d1] * nDim + d3]);
943 MFloat wjk = 0.5 * (slopes[PV->VV[d2] * nDim + d3] - slopes[PV->VV[d3] * nDim + d2]);
944 OijOjkSki += wij * wjk * ski;
945 }
946 }
947 }
948
949 MFloat X_omega = fabs(OijOjkSki / pow(RM_KOMEGA::betas * omega, F3));
950 MFloat f_beta = (F1 + RM_KOMEGA::fbeta1 * X_omega) / (F1 + RM_KOMEGA::fbeta2 * X_omega);
951 MFloat beta = RM_KOMEGA::beta * f_beta;
952
953 const MFloat P_k = P;
954 const MFloat P_o = RM_KOMEGA::alpha * omega / k * P;
955
956 const MFloat D_k = Re * RM_KOMEGA::betas * rho * k * omega;
957 const MFloat D_o = Re * beta * rho * POW2(omega);
958
959 // MFloat diff_o = F0;
960 // if(term > 0) {
961 // diff_o = rRe * rho * RM_KOMEGA::sigma_d / omega * term;
962 // }
963
964 const MFloat prodDest_K = P_k - D_k;
965 const MFloat prodDest_OMEGA = P_o - D_o;
966
967 rhs[CV->RHO_K] -= cellVol * (prodDest_K);
968 rhs[CV->RHO_OMEGA] -= cellVol * (prodDest_OMEGA);
969}
970
971#endif // FVSYSEQNRANS_H_
FvSysEqnNS::m_sutherlandConstant
MFloat m_sutherlandConstant
Definition: fvcartesiansyseqnns.h:303
FvSysEqnNS::getArray012
static constexpr std::array< MInt, nDim > getArray012()
Definition: fvcartesiansyseqnns.h:316
FvSysEqnNS::m_sutherlandPlusOneThermal
MFloat m_sutherlandPlusOneThermal
Definition: fvcartesiansyseqnns.h:304
FvSysEqnNS::m_sutherlandPlusOne
MFloat m_sutherlandPlusOne
Definition: fvcartesiansyseqnns.h:302
FvSysEqnNS::m_referenceTemperature
MFloat m_referenceTemperature
Definition: fvcartesiansyseqnns.h:301
FvSysEqnNS::m_F1BGammaMinusOne
MFloat m_F1BGammaMinusOne
Definition: fvcartesiansyseqnns.h:296
FvSysEqnNS::index1
static const MUint index1[nDim]
Definition: fvcartesiansyseqnns.h:353
FvSysEqnNS::index0
static const MUint index0[nDim]
Definition: fvcartesiansyseqnns.h:350
FvSysEqnNS::m_sutherlandConstantThermal
MFloat m_sutherlandConstantThermal
Definition: fvcartesiansyseqnns.h:305
FvSysEqnRANS::m_ransModel
static constexpr MInt m_ransModel
Definition: fvcartesiansyseqnrans.h:120
FvSysEqnRANS::viscousFlux
void viscousFlux(const MInt orientation, const MFloat A, const MBool isBndry, const MFloat *const surfaceCoords, const MFloat *const coord0, const MFloat *const coord1, const MFloat *const cellVars0, const MFloat *const cellVars1, const MFloat *const vars0, const MFloat *const vars1, const MFloat *const slope0, const MFloat *const slope1, const MFloat f0, const MFloat f1, MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:57
FvSysEqnRANS::Ausm
void Ausm(const MInt orientation, const MFloat upwindCoefficient, const MFloat A, const MFloat *const leftVars, const MFloat *const rightVars, const MFloat *const NotUsed(srfcCoeff), MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:233
FvSysEqnRANS::getArray012
static constexpr std::array< MInt, nDim > getArray012()
Definition: fvcartesiansyseqnns.h:316
FvSysEqnRANS::computeConservativeVariables
void computeConservativeVariables(const MFloat *const pvarsCell, MFloat *const cvarsCell, const MFloat *const NotUsed(avarsCell))
Definition: fvcartesiansyseqnrans.h:652
FvSysEqnRANS::viscousFluxFivePoint
void viscousFluxFivePoint(const MInt orientation, const MFloat A, const MFloat *const vars0, const MFloat *const vars1, const MFloat *const slope0, const MFloat *const slope1, const MFloat *const NotUsed(srfcCoeff), const MFloat f0, const MFloat f1, MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:330
FvSysEqnRANS::AusmBndryCorrection
void AusmBndryCorrection(const MInt orientation, const MFloat A, const MFloat *const leftVars, const MFloat *const rightVars, MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:300
FvSysEqnRANS::viscousFlux
void viscousFlux(const MInt orientation, const MFloat A, const MFloat *const vars0, const MFloat *const vars1, const MFloat *const slope0, const MFloat *const slope1, const MFloat *const srfcCoeff, const MFloat f0, const MFloat f1, MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:45
FvSysEqnRANS::computeVolumeForcesRANS_FS
void computeVolumeForcesRANS_FS(const MInt cellId, MFloat *vars, MFloat *rhs, MFloat cellVol, const MFloat *const slopes, MInt recData, const MInt *const recNghbr, const MFloat *const recConst, MInt noRecNghbr, MFloat)
Definition: fvcartesiansyseqnrans.h:781
FvSysEqnRANS::FV
FluxVariables * FV
Definition: fvcartesiansyseqnrans.h:158
FvSysEqnRANS::hasSC
static constexpr MBool hasSC
Definition: fvcartesiansyseqnrans.h:156
FvSysEqnRANS::computeVolumeForces
void computeVolumeForces(const MInt cellId, MFloat *vars, MFloat *rhs, const MFloat cellVol, const MFloat *const slopes, const MInt recData, const MInt *const recNghbr, const MFloat *const recConst, const MInt noRecNghbr, MFloat dist)
Definition: fvcartesiansyseqnrans.h:93
FvSysEqnRANS::CV
ConservativeVariables * CV
Definition: fvcartesiansyseqnrans.h:157
FvSysEqnRANS::computeVolumeForcesRANS_SA
void computeVolumeForcesRANS_SA(const MInt cellId, MFloat *vars, MFloat *rhs, const MFloat cellVol, const MFloat *const slopes, const MInt recData, const MInt *const, const MFloat *const, const MInt noRecNghbr, MFloat dist)
Definition: fvcartesiansyseqnrans.h:716
FvSysEqnRANS::viscousFluxThreePoint
void viscousFluxThreePoint(const MInt orientation, const MFloat A, const MBool isBndry, const MFloat *const surfaceCoords, const MFloat *const coord0, const MFloat *const coord1, const MFloat *const cellVars0, const MFloat *const cellVars1, const MFloat *const vars0, const MFloat *const vars1, const MFloat *const slope0, const MFloat *const slope1, const MFloat f0, const MFloat f1, MFloat *const flux)
Definition: fvcartesiansyseqnrans.h:511
FvSysEqnRANS::PV
PrimitiveVariables * PV
Definition: fvcartesiansyseqnrans.h:159
FvSysEqnRANS::conservativeSlopes
std::vector< std::vector< MFloat > > conservativeSlopes(const MFloat *const pvarsCell, const MFloat *const cvarsCell, const MFloat *const avarsCell, const MFloat *const slopesCell)
Definition: fvcartesiansyseqnrans.h:681
FvSysEqnRANS::m_noRansEquations
static constexpr MInt m_noRansEquations
Definition: fvcartesiansyseqnrans.h:119
FvSysEqnRANS::computePrimitiveVariables
void computePrimitiveVariables(const MFloat *const cvarsCell, MFloat *const pvarsCell, const MFloat *const NotUsed(avarsCell))
Definition: fvcartesiansyseqnrans.h:626
FvSysEqnRANS::computeVolumeForcesRANS_KOMEGA
void computeVolumeForcesRANS_KOMEGA(const MInt cellId, MFloat *vars, MFloat *rhs, const MFloat cellVol, const MFloat *const slopes, const MInt, const MInt *const, const MFloat *const, const MInt, MFloat)
Definition: fvcartesiansyseqnrans.h:901
no
@ no
Definition: enums.h:208
FIVE_POINT
@ FIVE_POINT
Definition: enums.h:184
THREE_POINT
@ THREE_POINT
Definition: enums.h:184
RANS_SA_DV
@ RANS_SA_DV
Definition: enums.h:54
RANS_FS
@ RANS_FS
Definition: enums.h:55
RANS_KOMEGA
@ RANS_KOMEGA
Definition: enums.h:56
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
MInt
int32_t MInt
Definition: maiatypes.h:62
MUint
uint32_t MUint
Definition: maiatypes.h:63
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54
FvSysEqnRANS::ConservativeVariables
Static indices for accessing conservative variables in nDim spatial dimensions.
Definition: fvcartesiansyseqnrans.h:166
FvSysEqnRANS::PrimitiveVariables
Static indices for accessing primitive variables in nDim spatial dimensions.
Definition: fvcartesiansyseqnrans.h:200
FvSysEqnRANS::PrimitiveVariables::varNames
static std::vector< MString > varNames
Definition: fvcartesiansyseqnrans.h:221
RANSModelConstants< RANS_FS >::faalpha
static constexpr MFloat faalpha
Definition: fvransmodelconstants.h:48
RANSModelConstants< RANS_FS >::fapsik2
static constexpr MFloat fapsik2
Definition: fvransmodelconstants.h:44
RANSModelConstants< RANS_FS >::fabetcs
static constexpr MFloat fabetcs
Definition: fvransmodelconstants.h:41
RANSModelConstants< RANS_FS >::fapsik1
static constexpr MFloat fapsik1
Definition: fvransmodelconstants.h:43
RANSModelConstants< RANS_FS >::fasigma
static constexpr MFloat fasigma
Definition: fvransmodelconstants.h:47
RANSModelConstants< RANS_FS >::facv1to3
static constexpr MFloat facv1to3
Definition: fvransmodelconstants.h:50
RANSModelConstants< RANS_FS >::fabetc
static constexpr MFloat fabetc
Definition: fvransmodelconstants.h:42
RANSModelConstants< RANS_FS >::fapsio2
static constexpr MFloat fapsio2
Definition: fvransmodelconstants.h:46
RANSModelConstants< RANS_FS >::fapsio1
static constexpr MFloat fapsio1
Definition: fvransmodelconstants.h:45
RANSModelConstants< RANS_KOMEGA >::sigma_o
static constexpr MFloat sigma_o
Definition: fvransmodelconstants.h:64
RANSModelConstants< RANS_KOMEGA >::sigma_k
static constexpr MFloat sigma_k
Definition: fvransmodelconstants.h:63
RANSModelConstants< RANS_KOMEGA >::beta
static constexpr MFloat beta
Definition: fvransmodelconstants.h:60
RANSModelConstants< RANS_KOMEGA >::betas
static constexpr MFloat betas
Definition: fvransmodelconstants.h:62
RANSModelConstants< RANS_KOMEGA >::alpha
static constexpr MFloat alpha
Definition: fvransmodelconstants.h:59
RANSModelConstants< RANS_KOMEGA >::fbeta1
static constexpr MFloat fbeta1
Definition: fvransmodelconstants.h:66
RANSModelConstants< RANS_KOMEGA >::fbeta2
static constexpr MFloat fbeta2
Definition: fvransmodelconstants.h:67
RANSModelConstants< RANS_SA_DV >::cw1
static constexpr MFloat cw1
Definition: fvransmodelconstants.h:28
RANSModelConstants< RANS_SA_DV >::Fsigma
static constexpr MFloat Fsigma
Definition: fvransmodelconstants.h:24
RANSModelConstants< RANS_SA_DV >::cv1to3
static constexpr MFloat cv1to3
Definition: fvransmodelconstants.h:31
RANSModelConstants< RANS_SA_DV >::cw3to6
static constexpr MFloat cw3to6
Definition: fvransmodelconstants.h:30
RANSModelConstants< RANS_SA_DV >::Ft2
static constexpr MFloat Ft2
Definition: fvransmodelconstants.h:32
RANSModelConstants< RANS_SA_DV >::cb2
static constexpr MFloat cb2
Definition: fvransmodelconstants.h:25
RANSModelConstants< RANS_SA_DV >::Fkap2
static constexpr MFloat Fkap2
Definition: fvransmodelconstants.h:26
RANSModelConstants< RANS_SA_DV >::cw2
static constexpr MFloat cw2
Definition: fvransmodelconstants.h:29
RANSModelConstants< RANS_SA_DV >::cb1
static constexpr MFloat cb1
Definition: fvransmodelconstants.h:27