1155 {
1156
1158
1160
1162
1163
1164 static constexpr MFloat minMFloat =
1165 1e-16;
1171
1176
1181
1188
1191
1192
1197
1200
1201
1202
1203 const MFloat dTKEdxi = F1B2 * (TKE[IPJP] + TKE[IPJ] - TKE[IJP] - TKE[IJ]);
1204 const MFloat dTKEdet = F1B2 * (TKE[IPJP] + TKE[IJP] - TKE[IPJ] - TKE[IJ]);
1205
1206 const MFloat depsdxi = F1B2 * (eps[IPJP] + eps[IPJ] - eps[IJP] - eps[IJ]);
1207 const MFloat depsdet = F1B2 * (eps[IPJP] + eps[IJP] - eps[IPJ] - eps[IJ]);
1208
1209 const MFloat metricTerms = cornerMetrics[
xsd * 2 +
xsd] * cornerMetrics[
ysd * 2 +
xsd]
1210 + cornerMetrics[
xsd * 2 +
ysd] * cornerMetrics[
ysd * 2 +
ysd];
1211
1213
1214 const MFloat sax1 = invCornerJac * dTKEdet * metricTerms;
1215 const MFloat sax2 = invCornerJac * depsdet * metricTerms;
1216 const MFloat say1 = invCornerJac * dTKEdxi * metricTerms;
1217 const MFloat say2 = invCornerJac * depsdxi * metricTerms;
1218
1219 eflux[0][IJ] = sax1;
1220 eflux[1][IJ] = sax2;
1221
1222 fflux[0][IJ] = say1;
1223 fflux[1][IJ] = say2;
1224
1229 const MFloat dudxi = F1B2 * (u[IPJP] + u[IPJ] - u[IJP] - u[IJ]);
1230 const MFloat dudet = F1B2 * (u[IPJP] + u[IJP] - u[IPJ] - u[IJ]);
1231
1232 const MFloat dvdxi = F1B2 * (v[IPJP] + v[IPJ] - v[IJP] - v[IJ]);
1233 const MFloat dvdet = F1B2 * (v[IPJP] + v[IJP] - v[IPJ] - v[IJ]);
1234
1235 const MFloat S11 = dudxi * cornerMetrics[
xsd * 2 +
xsd] + dudet * cornerMetrics[
ysd * 2 +
xsd];
1237 * (dudxi * cornerMetrics[
xsd * 2 +
ysd] + dudet * cornerMetrics[
ysd * 2 +
ysd]
1238 + dvdxi * cornerMetrics[
xsd * 2 +
xsd] + dvdet * cornerMetrics[
ysd * 2 +
xsd]);
1239 const MFloat S22 = dvdxi * cornerMetrics[
xsd * 2 +
ysd] + dvdet * cornerMetrics[
ysd * 2 +
ysd];
1240 const MFloat S = S11 * S11 + 2 * S12 * S12 + S22 * S22;
1241
1242 const MFloat muTurbAvg = F1B4 * (muTurb[IJP] + muTurb[IPJP] + muTurb[IJ] + muTurb[IPJ]);
1243
1244
1245
1246
1247 const MFloat mueOverRe = rRe0 * invCornerJac * muTurbAvg /
POW2(fac_nonDim);
1248
1250 }
1251 }
1252 }
1253
1254
1255
1256
1259 }
1260
1271
1272 const MFloat dTKEdxi = TKE[IPJ] - TKE[IJ];
1273 const MFloat depsdxi = eps[IPJ] - eps[IJ];
1274
1275 const MFloat muLamAvg_xi = F1B2 * (muLam[IJ] + muLam[IPJ]);
1276 const MFloat muTurbAvg_xi = F1B2 * (muTurb[IJ] + muTurb[IPJ]);
1277
1278
1280 0.5
1282 const MFloat Frj = rRe0 * surfJac;
1285
1288
1289 MFloat limiterVisc1 = 1.0;
1290 MFloat limiterVisc2 = 1.0;
1292
1294 limiterVisc1 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_k_xi));
1295 limiterVisc2 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_eps_xi));
1296 }
1297
1298 sa_1flux[0][IJ] = Frj * mu_k_xi * (temp_TKE + F1B2 * (eflux[0][IJ] + eflux[0][IJM])) * limiterVisc1;
1299 sa_2flux[0][IJ] = Frj * mu_eps_xi * (temp_eps + F1B2 * (eflux[1][IJ] + eflux[1][IJM])) * limiterVisc2;
1300 }
1301 }
1302
1310
1311 const MFloat dTKEdet = TKE[IJP] - TKE[IJ];
1312 const MFloat depsdet = eps[IJP] - eps[IJ];
1313
1314 const MFloat muLamAvg_eta = F1B2 * (muLam[IJ] + muLam[IJP]);
1315 const MFloat muTurbAvg_eta = F1B2 * (muTurb[IJ] + muTurb[IJP]);
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1329 0.5
1331 const MFloat Frj = rRe0 * surfJac;
1334
1337
1338 MFloat limiterVisc1 = 1.0;
1339 MFloat limiterVisc2 = 1.0;
1341
1343 limiterVisc1 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_k_eta));
1344 limiterVisc2 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_eps_eta));
1345 }
1346
1347 sa_1flux[1][IJ] = Frj * mu_k_eta * (temp_TKE + F1B2 * (fflux[0][IJ] + fflux[0][IMJ])) * limiterVisc1;
1348 sa_2flux[1][IJ] = Frj * mu_eps_eta * (temp_eps + F1B2 * (fflux[1][IJ] + fflux[1][IMJ])) * limiterVisc2;
1349 }
1350 }
1351
1352
1360
1361
1362
1364
1368 MFloat limiterVisc1 = 1.0;
1369 MFloat limiterVisc2 = 1.0;
1370
1371
1372
1373
1374
1375
1376
1377 const MFloat diffusion_TKE =
1378 (sa_1flux[0][IJ] - sa_1flux[0][IMJ]) + (sa_1flux[1][IJ] - sa_1flux[1][IJM]) * limiterVisc1;
1379
1380 const MFloat diffusion_eps =
1381 (sa_2flux[0][IJ] - sa_2flux[0][IMJ]) + (sa_2flux[1][IJ] - sa_2flux[1][IJM]) * limiterVisc2;
1382
1387
1392 } else {
1393
1394 const MFloat dudxi = 0.5 * (u[IPJ] - u[IMJ]);
1395 const MFloat dudet = 0.5 * (u[IJP] - u[IJM]);
1396 const MFloat dvdxi = 0.5 * (v[IPJ] - v[IMJ]);
1397 const MFloat dvdet = 0.5 * (v[IJP] - v[IJM]);
1398
1407
1408
1409
1410 P = ((rRe0 /
POW2(fac_nonDim)) * invCellJac * muTurb[IJ]
1411 * (2.0 *
POW2(dudx) + 2.0 *
POW2(dvdy) +
POW2(dudy + dvdx) - fac * F2B3 *
POW2(dudx + dvdy))
1412 - fac * F2B3 * rho[IJ] * TKE[IJ] * (dudx + dvdy))
1413 * invCellJac;
1414 }
1415
1416
1418 fac_nonDim *
m_solver->
m_Re0 * rho[IJ] *
POW2(TKE[IJ]) / muLam[IJ] / std::max(eps[IJ], minMFloat);
1419 const MFloat f2 = 1 - 0.222222222222222 * exp(-
POW2(Ret) * 0.027777777777777777 );
1420
1421 const MFloat Mt2 = 1.5 * TKE[IJ] * rho[IJ] / (gamma *
p[IJ]);
1422 const MFloat Fdist2 = 1.0 /
POW2(wallDist[IJ]);
1423 const MFloat k_diss1 = fac_nonDim * rho[IJ] * eps[IJ] * (1 + Mt2);
1424 const MFloat k_diss2 = 2.0 * rRe0 * muLam[IJ] * TKE[IJ] * Fdist2;
1425
1428 const MFloat yp = utau[IJ] * wallDist[IJ];
1430
1431
1435 const MFloat invRhoEps = 1.0 / std::max(rho[IJ] * eps[IJ], minMFloat);
1436
1439 + fac_nonDim * (1 + Mt2);
1440 const MFloat f2_limit = (RHS1 + k_diss2 * invRhoEps * (1 - exp_wall)) / (fac_nonDim *
RM_KEPS::C2);
1441 f2_ = std::min(1.0, std::max(f2, f2_limit));
1442 const MFloat diff = std::max(0.0, f2_limit - 1.0);
1443 f3 = std::max(0.0, 1.0 - fac_nonDim *
RM_KEPS::C2 * diff / (k_diss2 * invRhoEps * (1 - exp_wall)));
1444 }
1445
1446 const MFloat k_diss = k_diss1 + f3 * k_diss2;
1447 const MFloat eps_diss = tau * (fac_nonDim *
RM_KEPS::C2 * f2_ * rho[IJ] * eps[IJ] + f3 * k_diss2 * exp_wall);
1448
1451
1455
1456
1457
1462
1463
1464
1465 }
1466 }
1467
1468
1471 else
1473}
void diffusiveFluxCorrection()