diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C index 61c67a8726..379b109867 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2019 OpenFOAM Foundation - Copyright (C) 2017-2022 OpenCFD Ltd. + Copyright (C) 2017-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -189,105 +189,150 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const tmp tnutw = turbModel.nut(patchi); - const scalarField& nutw = tnutw(); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); const scalarField& y = turbModel.y()[patchi]; + const labelUList& faceCells = patch.faceCells(); + const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); const tmp tk = turbModel.k(); const volScalarField& k = tk(); - const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi]; - - const scalarField magGradUw(mag(Uw.snGrad())); - - const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); - const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75); - const scalar kappa = wallCoeffs_.kappa(); - const scalar yPlusLam = wallCoeffs_.yPlusLam(); - - // Set epsilon and G - forAll(nutw, facei) + // Calculate y-plus + const auto yPlus = [&](const label facei) -> scalar { - const label celli = patch.faceCells()[facei]; + return + ( + Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nuw[facei] + ); + }; - const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei]; + // Contribution from the viscous sublayer + const auto epsilonVis = [&](const label facei) -> scalar + { + return + ( + cornerWeights[facei]*2.0*k[faceCells[facei]]*nuw[facei] + / sqr(y[facei]) + ); + }; - const scalar w = cornerWeights[facei]; + // Contribution from the inertial sublayer + const auto epsilonLog = [&](const label facei) -> scalar + { + return + ( + cornerWeights[facei]*Cmu75*pow(k[faceCells[facei]], 1.5) + / (kappa*y[facei]) + ); + }; - // Contribution from the viscous sublayer - const scalar epsilonVis = w*2.0*k[celli]*nuw[facei]/sqr(y[facei]); - - // Contribution from the inertial sublayer - const scalar epsilonLog = w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]); - - switch (blender_) + switch (blender_) + { + case blenderType::STEPWISE: { - case blenderType::STEPWISE: + forAll(faceCells, facei) { - if (lowReCorrection_ && yPlus < yPlusLam) + if (lowReCorrection_ && yPlus(facei) < yPlusLam) { - epsilon0[celli] += epsilonVis; + epsilon0[faceCells[facei]] += epsilonVis(facei); } else { - epsilon0[celli] += epsilonLog; + epsilon0[faceCells[facei]] += epsilonLog(facei); } - break; - } - - case blenderType::BINOMIAL: - { - // (ME:Eqs. 15-16) - epsilon0[celli] += - pow - ( - pow(epsilonVis, n_) + pow(epsilonLog, n_), - scalar(1)/n_ - ); - break; - } - - case blenderType::MAX: - { - // (PH:Eq. 27) - epsilon0[celli] += max(epsilonVis, epsilonLog); - break; - } - - case blenderType::EXPONENTIAL: - { - // (PH:p. 193) - const scalar Gamma = 0.001*pow4(yPlus)/(scalar(1) + yPlus); - const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); - epsilon0[celli] += - epsilonVis*exp(-Gamma) + epsilonLog*exp(-invGamma); - break; - } - - case blenderType::TANH: - { - // (KAS:Eqs. 33-34) - const scalar phiTanh = tanh(pow4(0.1*yPlus)); - const scalar b1 = epsilonVis + epsilonLog; - const scalar b2 = - pow(pow(epsilonVis, 1.2) + pow(epsilonLog, 1.2), 1.0/1.2); - - epsilon0[celli] += phiTanh*b1 + (1 - phiTanh)*b2; - break; } + break; } - if (!lowReCorrection_ || (yPlus > yPlusLam)) + case blenderType::BINOMIAL: { - G0[celli] += - w + forAll(faceCells, facei) + { + // (ME:Eqs. 15-16) + epsilon0[faceCells[facei]] += + pow + ( + pow(epsilonVis(facei), n_) + pow(epsilonLog(facei), n_), + scalar(1)/n_ + ); + } + break; + } + + case blenderType::MAX: + { + forAll(faceCells, facei) + { + // (PH:Eq. 27) + epsilon0[faceCells[facei]] += + max(epsilonVis(facei), epsilonLog(facei)); + } + break; + } + + case blenderType::EXPONENTIAL: + { + forAll(faceCells, facei) + { + // (PH:p. 193) + const scalar yPlusFace = yPlus(facei); + const scalar Gamma = + 0.001*pow4(yPlusFace)/(scalar(1) + yPlusFace); + const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); + + epsilon0[faceCells[facei]] += + ( + epsilonVis(facei)*exp(-Gamma) + + epsilonLog(facei)*exp(-invGamma) + ); + } + break; + } + + case blenderType::TANH: + { + forAll(faceCells, facei) + { + // (KAS:Eqs. 33-34) + const scalar epsilonVisFace = epsilonVis(facei); + const scalar epsilonLogFace = epsilonLog(facei); + const scalar b1 = epsilonVisFace + epsilonLogFace; + const scalar b2 = + pow + ( + pow(epsilonVisFace, 1.2) + pow(epsilonLogFace, 1.2), + 1.0/1.2 + ); + const scalar phiTanh = tanh(pow4(0.1*yPlus(facei))); + + epsilon0[faceCells[facei]] += phiTanh*b1 + (1 - phiTanh)*b2; + } + break; + } + } + + const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi]; + const scalarField magGradUw(mag(Uw.snGrad())); + + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); + + forAll(faceCells, facei) + { + if (!lowReCorrection_ || (yPlus(facei) > yPlusLam)) + { + G0[faceCells[facei]] += + cornerWeights[facei] *(nutw[facei] + nuw[facei]) *magGradUw[facei] - *Cmu25*sqrt(k[celli]) + *Cmu25*sqrt(k[faceCells[facei]]) /(kappa*y[facei]); } } diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C index af072fe580..76feb57497 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2021 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -32,7 +32,6 @@ License #include "volFields.H" #include "addToRunTimeSelectionTable.H" - // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // Foam::tmp @@ -40,6 +39,10 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const { const label patchi = patch().index(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + const auto& turbModel = db().lookupObject ( IOobject::groupName @@ -52,88 +55,107 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi]; const scalarField magUp(mag(Uw.patchInternalField() - Uw)); - const tmp tnuw = turbModel.nu(patchi); - const scalarField& nuw = tnuw(); - - const scalar kappa = wallCoeffs_.kappa(); - const scalar E = wallCoeffs_.E(); - const scalar yPlusLam = wallCoeffs_.yPlusLam(); - tmp tyPlus = calcYPlus(magUp); const scalarField& yPlus = tyPlus(); + // Viscous sublayer contribution + const tmp tnutVis = turbModel.nu(patchi); + const scalarField& nutVis = tnutVis(); + + // Inertial sublayer contribution + const auto nutLog = [&](const label facei) -> scalar + { + const scalar yPlusFace = yPlus[facei]; + return + ( + nutVis[facei]*yPlusFace*kappa/log(max(E*yPlusFace, 1 + 1e-4)) + ); + }; + auto tnutw = tmp::New(patch().size(), Zero); auto& nutw = tnutw.ref(); - forAll(yPlus, facei) + switch (blender_) { - // Viscous sublayer contribution - const scalar nutVis = 0; - - // Inertial sublayer contribution - const scalar nutLog = - nuw[facei] - *(yPlus[facei]*kappa/log(max(E*yPlus[facei], 1 + 1e-4)) - 1.0); - - switch (blender_) + case blenderType::STEPWISE: { - case blenderType::STEPWISE: + forAll(nutw, facei) { if (yPlus[facei] > yPlusLam) { - nutw[facei] = nutLog; + nutw[facei] = nutLog(facei); } else { - nutw[facei] = nutVis; + nutw[facei] = nutVis[facei]; } - break; } + break; + } - case blenderType::MAX: + case blenderType::MAX: + { + forAll(nutw, facei) { // (PH:Eq. 27) - nutw[facei] = max(nutVis, nutLog); - break; + nutw[facei] = max(nutVis[facei], nutLog(facei)); } + break; + } - case blenderType::BINOMIAL: + case blenderType::BINOMIAL: + { + forAll(nutw, facei) { // (ME:Eqs. 15-16) nutw[facei] = pow ( - pow(nutVis, n_) + pow(nutLog, n_), + pow(nutVis[facei], n_) + pow(nutLog(facei), n_), scalar(1)/n_ ); - break; } + break; + } - case blenderType::EXPONENTIAL: + case blenderType::EXPONENTIAL: + { + forAll(nutw, facei) { // (PH:Eq. 31) const scalar Gamma = 0.01*pow4(yPlus[facei])/(1 + 5*yPlus[facei]); const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); - nutw[facei] = nutVis*exp(-Gamma) + nutLog*exp(-invGamma); - break; + nutw[facei] = + nutVis[facei]*exp(-Gamma) + nutLog(facei)*exp(-invGamma); } + break; + } - case blenderType::TANH: + case blenderType::TANH: + { + forAll(nutw, facei) { // (KAS:Eqs. 33-34) - const scalar phiTanh = tanh(pow4(0.1*yPlus[facei])); - const scalar b1 = nutVis + nutLog; + const scalar nutLogFace = nutLog(facei); + const scalar b1 = nutVis[facei] + nutLogFace; const scalar b2 = - pow(pow(nutVis, 1.2) + pow(nutLog, 1.2), 1.0/1.2); + pow + ( + pow(nutVis[facei], 1.2) + pow(nutLogFace, 1.2), + 1.0/1.2 + ); + const scalar phiTanh = tanh(pow4(0.1*yPlus[facei])); nutw[facei] = phiTanh*b1 + (1 - phiTanh)*b2; - break; } + break; } } + nutw -= nutVis; + return tnutw; } diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C index 8d9e40ea65..83ecf54a7b 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016, 2019 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -40,6 +40,11 @@ calcNut() const { const label patchi = patch().index(); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + const auto& turbModel = db().lookupObject ( IOobject::groupName @@ -49,93 +54,118 @@ calcNut() const ) ); + const labelUList& faceCells = patch().faceCells(); + const scalarField& y = turbModel.y()[patchi]; const tmp tk = turbModel.k(); const volScalarField& k = tk(); - const tmp tnuw = turbModel.nu(patchi); - const scalarField& nuw = tnuw(); + // Viscous sublayer contribution + const tmp tnutVis = turbModel.nu(patchi); + const scalarField& nutVis = tnutVis(); - const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); - const scalar kappa = wallCoeffs_.kappa(); - const scalar E = wallCoeffs_.E(); - const scalar yPlusLam = wallCoeffs_.yPlusLam(); + // Calculate y-plus + const auto yPlus = [&](const label facei) -> scalar + { + return (Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nutVis[facei]); + }; + + // Inertial sublayer contribution + const auto nutLog = [&](const label facei) -> scalar + { + const scalar yPlusFace = yPlus(facei); + return + ( + nutVis[facei]*yPlusFace*kappa + / log(max(E*yPlusFace, 1 + 1e-4)) + ); + }; auto tnutw = tmp::New(patch().size(), Zero); auto& nutw = tnutw.ref(); - forAll(nutw, facei) + switch (blender_) { - const label celli = patch().faceCells()[facei]; - - const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei]; - - // Viscous sublayer contribution - const scalar nutVis = 0; - - // Inertial sublayer contribution - const scalar nutLog = - nuw[facei]*(yPlus*kappa/log(max(E*yPlus, 1 + 1e-4)) - 1.0); - - switch (blender_) + case blenderType::STEPWISE: { - case blenderType::STEPWISE: + forAll(nutw, facei) { - if (yPlus > yPlusLam) + if (yPlus(facei) > yPlusLam) { - nutw[facei] = nutLog; + nutw[facei] = nutLog(facei); } else { - nutw[facei] = nutVis; + nutw[facei] = nutVis[facei]; } - break; } + break; + } - case blenderType::MAX: + case blenderType::MAX: + { + forAll(nutw, facei) { // (PH:Eq. 27) - nutw[facei] = max(nutVis, nutLog); - break; + nutw[facei] = max(nutVis[facei], nutLog(facei)); } + break; + } - case blenderType::BINOMIAL: + case blenderType::BINOMIAL: + { + forAll(nutw, facei) { // (ME:Eqs. 15-16) nutw[facei] = pow ( - pow(nutVis, n_) + pow(nutLog, n_), + pow(nutVis[facei], n_) + pow(nutLog(facei), n_), scalar(1)/n_ ); - break; } + break; + } - case blenderType::EXPONENTIAL: + case blenderType::EXPONENTIAL: + { + forAll(nutw, facei) { // (PH:Eq. 31) - const scalar Gamma = 0.01*pow4(yPlus)/(1 + 5*yPlus); + const scalar yPlusFace = yPlus(facei); + const scalar Gamma = 0.01*pow4(yPlusFace)/(1 + 5*yPlusFace); const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); - nutw[facei] = nutVis*exp(-Gamma) + nutLog*exp(-invGamma); - break; + nutw[facei] = + nutVis[facei]*exp(-Gamma) + nutLog(facei)*exp(-invGamma); } + break; + } - case blenderType::TANH: + case blenderType::TANH: + { + forAll(nutw, facei) { // (KAS:Eqs. 33-34) - const scalar phiTanh = tanh(pow4(0.1*yPlus)); - const scalar b1 = nutVis + nutLog; + const scalar nutLogFace = nutLog(facei); + const scalar b1 = nutVis[facei] + nutLogFace; const scalar b2 = - pow(pow(nutVis, 1.2) + pow(nutLog, 1.2), 1.0/1.2); + pow + ( + pow(nutVis[facei], 1.2) + pow(nutLogFace, 1.2), + 1.0/1.2 + ); + const scalar phiTanh = tanh(pow4(0.1*yPlus(facei))); nutw[facei] = phiTanh*b1 + (1 - phiTanh)*b2; - break; } + break; } } + nutw -= nutVis; + return tnutw; } @@ -231,8 +261,8 @@ yPlus() const tmp tkwc = k.boundaryField()[patchi].patchInternalField(); const scalarField& kwc = tkwc(); - tmp tnuw = turbModel.nu(patchi); - const scalarField& nuw = tnuw(); + tmp tnutVis = turbModel.nu(patchi); + const scalarField& nutVis = tnutVis(); tmp tnuEff = turbModel.nuEff(patchi); const scalarField& nuEff = tnuEff(); @@ -249,13 +279,13 @@ yPlus() const forAll(yPlus, facei) { // inertial sublayer - yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nuw[facei]; + yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nutVis[facei]; if (yPlusLam > yPlus[facei]) { // viscous sublayer yPlus[facei] = - y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nuw[facei]; + y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nutVis[facei]; } } diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 7339e1bd88..f7769ebc20 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016, 2019 OpenFOAM Foundation - Copyright (C) 2017-2022 OpenCFD Ltd. + Copyright (C) 2017-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -188,102 +188,146 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const tmp tnutw = turbModel.nut(patchi); - const scalarField& nutw = tnutw(); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); const scalarField& y = turbModel.y()[patchi]; + const labelUList& faceCells = patch.faceCells(); + + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); + const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); const tmp tk = turbModel.k(); const volScalarField& k = tk(); - const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi]; - - const scalarField magGradUw(mag(Uw.snGrad())); - - const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); - const scalar kappa = wallCoeffs_.kappa(); - const scalar yPlusLam = wallCoeffs_.yPlusLam(); - - // Set omega and G - forAll(nutw, facei) + // Calculate y-plus + const auto yPlus = [&](const label facei) -> scalar { - const label celli = patch.faceCells()[facei]; - const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei]; - const scalar w = cornerWeights[facei]; + return + ( + Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nuw[facei] + ); + }; - // Contribution from the viscous sublayer - const scalar omegaVis = 6.0*nuw[facei]/(beta1_*sqr(y[facei])); + // Contribution from the viscous sublayer + const auto omegaVis = [&](const label facei) -> scalar + { + return + ( + cornerWeights[facei]*6.0*nuw[facei]/(beta1_*sqr(y[facei])) + ); + }; - // Contribution from the inertial sublayer - const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa*y[facei]); + // Contribution from the inertial sublayer + const auto omegaLog = [&](const label facei) -> scalar + { + return + ( + cornerWeights[facei]*sqrt(k[faceCells[facei]]) + / (Cmu25*kappa*y[facei]) + ); + }; - switch (blender_) + switch (blender_) + { + case blenderType::STEPWISE: { - case blenderType::STEPWISE: + forAll(faceCells, facei) { - if (yPlus > yPlusLam) + if (yPlus(facei) > yPlusLam) { - omega0[celli] += w*omegaLog; + omega0[faceCells[facei]] += omegaLog(facei); } else { - omega0[celli] += w*omegaVis; + omega0[faceCells[facei]] += omegaVis(facei); } - break; - } - - case blenderType::BINOMIAL: - { - omega0[celli] += - w*pow - ( - pow(omegaVis, n_) + pow(omegaLog, n_), - scalar(1)/n_ - ); - break; - } - - case blenderType::MAX: - { - // (PH:Eq. 27) - omega0[celli] += max(omegaVis, omegaLog); - break; - } - - case blenderType::EXPONENTIAL: - { - // (PH:Eq. 31) - const scalar Gamma = 0.01*pow4(yPlus)/(1 + 5*yPlus); - const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); - - omega0[celli] += - w*(omegaVis*exp(-Gamma) + omegaLog*exp(-invGamma)); - break; - } - - case blenderType::TANH: - { - // (KAS:Eqs. 33-34) - const scalar phiTanh = tanh(pow4(0.1*yPlus)); - const scalar b1 = omegaVis + omegaLog; - const scalar b2 = - pow(pow(omegaVis, 1.2) + pow(omegaLog, 1.2), 1.0/1.2); - - omega0[celli] += phiTanh*b1 + (1 - phiTanh)*b2; - break; } + break; } - if (!(blender_ == blenderType::STEPWISE) || yPlus > yPlusLam) + case blenderType::BINOMIAL: { - G0[celli] += - w + forAll(faceCells, facei) + { + omega0[faceCells[facei]] += + pow + ( + pow(omegaVis(facei), n_) + pow(omegaLog(facei), n_), + scalar(1)/n_ + ); + } + break; + } + + case blenderType::MAX: + { + forAll(faceCells, facei) + { + // (PH:Eq. 27) + omega0[faceCells[facei]] += + max(omegaVis(facei), omegaLog(facei)); + } + break; + } + + case blenderType::EXPONENTIAL: + { + forAll(faceCells, facei) + { + // (PH:Eq. 31) + const scalar yPlusFace = yPlus(facei); + const scalar Gamma = 0.01*pow4(yPlusFace)/(1 + 5*yPlusFace); + const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL); + + omega0[faceCells[facei]] += + ( + omegaVis(facei)*exp(-Gamma) + + omegaLog(facei)*exp(-invGamma) + ); + } + break; + } + + case blenderType::TANH: + { + forAll(faceCells, facei) + { + // (KAS:Eqs. 33-34) + const scalar omegaVisFace = omegaVis(facei); + const scalar omegaLogFace = omegaLog(facei); + const scalar b1 = omegaVisFace + omegaLogFace; + const scalar b2 = + pow + ( + pow(omegaVisFace, 1.2) + pow(omegaLogFace, 1.2), + 1.0/1.2 + ); + const scalar phiTanh = tanh(pow4(0.1*yPlus(facei))); + + omega0[faceCells[facei]] += phiTanh*b1 + (1 - phiTanh)*b2; + } + break; + } + } + + const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi]; + const scalarField magGradUw(mag(Uw.snGrad())); + + forAll(faceCells, facei) + { + if (!(blender_ == blenderType::STEPWISE) || yPlus(facei) > yPlusLam) + { + G0[faceCells[facei]] += + cornerWeights[facei] *(nutw[facei] + nuw[facei]) *magGradUw[facei] - *Cmu25*sqrt(k[celli]) + *Cmu25*sqrt(k[faceCells[facei]]) /(kappa*y[facei]); } } diff --git a/src/thermoTools/derivedFvPatchFields/wallFunctions/sorptionWallFunction/sorptionWallFunctionFvPatchScalarField.C b/src/thermoTools/derivedFvPatchFields/wallFunctions/sorptionWallFunction/sorptionWallFunctionFvPatchScalarField.C index 8e4a513a37..1126f36c4c 100644 --- a/src/thermoTools/derivedFvPatchFields/wallFunctions/sorptionWallFunction/sorptionWallFunctionFvPatchScalarField.C +++ b/src/thermoTools/derivedFvPatchFields/wallFunctions/sorptionWallFunction/sorptionWallFunctionFvPatchScalarField.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2022 OpenCFD Ltd. + Copyright (C) 2022-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -64,9 +64,16 @@ static scalar calcYStarLam tmp sorptionWallFunctionFvPatchScalarField::yPlus() const { - // Calculate fields of interest const label patchi = patch().index(); + // Calculate the empirical constant given by (Jayatilleke, 1966) (FDC:Eq. 6) + const scalar Pc = + 9.24*(pow(Sc_/Sct_, 0.75) - 1)*(1 + 0.28*exp(-0.007*Sc_/Sct_)); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + + // Calculate fields of interest const auto& k = db().lookupObject(kName_); tmp tkwc = k.boundaryField()[patchi].patchInternalField(); const scalarField& kwc = tkwc.cref(); @@ -79,86 +86,116 @@ tmp sorptionWallFunctionFvPatchScalarField::yPlus() const tmp tywc = y.boundaryField()[patchi].patchInternalField(); const scalarField& ywc = tywc.cref(); + // (FDC:Eq. 3) + const auto yStar = [&](const label facei) -> scalar + { + return + ( + Cmu25*sqrt(kwc[facei])*ywc[facei]/nuwc[facei] + ); + }; - // Calculate the empirical constant given by (Jayatilleke, 1966) (FDC:Eq. 6) - const scalar Pc = - 9.24*(pow(Sc_/Sct_, 0.75) - 1)*(1 + 0.28*exp(-0.007*Sc_/Sct_)); - const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); - const scalar kappa = wallCoeffs_.kappa(); - const scalar E = wallCoeffs_.E(); + // (FDC:Eq. 4) + const auto yPlusVis = [&](const label facei) -> scalar + { + return (Sc_*yStar(facei)); + }; - auto tyPlus = tmp::New(patch().size(), Zero); + // (FDC:Eq. 5) + const auto yPlusLog = [&](const label facei) -> scalar + { + return + ( + Sct_*(log(max(E*yStar(facei), 1 + 1e-4))/kappa + Pc) + ); + }; + + auto tyPlus = tmp::New(patch().size()); auto& yPlus = tyPlus.ref(); - forAll(yPlus, facei) + switch (blender_) { - // (FDC:Eq. 3) - const scalar yStar = Cmu25*sqrt(kwc[facei])*ywc[facei]/nuwc[facei]; - - // (FDC:Eq. 4) - const scalar yPlusVis = Sc_*yStar; - - // (FDC:Eq. 5) - const scalar yPlusLog = Sct_*(log(max(E*yStar, 1 + 1e-4))/kappa + Pc); - - switch (blender_) + case blenderType::EXPONENTIAL: { - case blenderType::EXPONENTIAL: + forAll(yPlus, facei) { // (FDC:Eq. 2) + const scalar yStarFace = yStar(facei); const scalar Gamma = - 0.01*pow4(Sc_*yStar)/(1 + 5*pow3(Sc_)*yStar); + 0.01*pow4(Sc_*yStarFace)/(1 + 5*pow3(Sc_)*yStarFace); const scalar invGamma = scalar(1)/max(Gamma, ROOTVSMALL); // (FDC:Eq. 1) - yPlus[facei] = yPlusVis*exp(-Gamma) + yPlusLog*exp(-invGamma); - break; + yPlus[facei] = + ( + yPlusVis(facei)*exp(-Gamma) + + yPlusLog(facei)*exp(-invGamma) + ); } + break; + } - case blenderType::STEPWISE: + case blenderType::STEPWISE: + { + static const scalar yStarLam = + calcYStarLam(kappa, E, Sc_, Sct_, Pc); + + forAll(yPlus, facei) { - static const scalar yStarLam = - calcYStarLam(kappa, E, Sc_, Sct_, Pc); - // (F:Eq. 5.3) - if (yStar < yStarLam) + if (yStar(facei) < yStarLam) { - yPlus[facei] = yPlusVis; + yPlus[facei] = yPlusVis(facei); } else { - yPlus[facei] = yPlusLog; + yPlus[facei] = yPlusLog(facei); } - break; } + break; + } - case blenderType::BINOMIAL: + case blenderType::BINOMIAL: + { + forAll(yPlus, facei) { yPlus[facei] = pow ( - pow(yPlusVis, n_) + pow(yPlusLog, n_), + pow(yPlusVis(facei), n_) + pow(yPlusLog(facei), n_), scalar(1)/n_ ); - break; } + break; + } - case blenderType::MAX: + case blenderType::MAX: + { + forAll(yPlus, facei) { - yPlus[facei] = max(yPlusVis, yPlusLog); - break; + yPlus[facei] = max(yPlusVis(facei), yPlusLog(facei)); } + break; + } - case blenderType::TANH: + case blenderType::TANH: + { + forAll(yPlus, facei) { - const scalar phiTanh = tanh(pow4(0.1*yStar)); - const scalar b1 = yPlusVis + yPlusLog; + const scalar yPlusVisFace = yPlusVis(facei); + const scalar yPlusLogFace = yPlusLog(facei); + const scalar b1 = yPlusVisFace + yPlusLogFace; const scalar b2 = - pow(pow(yPlusVis, 1.2) + pow(yPlusLog, 1.2), 1.0/1.2); + pow + ( + pow(yPlusVisFace, 1.2) + pow(yPlusLogFace, 1.2), + 1.0/1.2 + ); + const scalar phiTanh = tanh(pow4(0.1*yStar(facei))); yPlus[facei] = phiTanh*b1 + (1 - phiTanh)*b2; - break; } + break; } }