47 alphaK1_(dimensioned<scalar>::lookupOrAddToDict(
51 alphaK2_(dimensioned<scalar>::lookupOrAddToDict(
55 alphaOmega1_(dimensioned<scalar>::lookupOrAddToDict(
59 alphaOmega2_(dimensioned<scalar>::lookupOrAddToDict(
63 gamma1_(dimensioned<scalar>::lookupOrAddToDict(
67 gamma2_(dimensioned<scalar>::lookupOrAddToDict(
71 beta1_(dimensioned<scalar>::lookupOrAddToDict(
75 beta2_(dimensioned<scalar>::lookupOrAddToDict(
79 betaStar_(dimensioned<scalar>::lookupOrAddToDict(
83 a1_(dimensioned<scalar>::lookupOrAddToDict(
87 b1_(dimensioned<scalar>::lookupOrAddToDict(
91 c1_(dimensioned<scalar>::lookupOrAddToDict(
95 F3_(Switch::lookupOrAddToDict(
100 omega_(const_cast<volScalarField&>(
101 mesh_.thisDb().lookupObject<volScalarField>(
"omega"))),
105 mesh.time().timeName(),
110 dimensionedScalar(
"omegaRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
111 zeroGradientFvPatchField<scalar>::typeName),
113 const_cast<volScalarField&>(
114 mesh_.thisDb().lookupObject<volScalarField>(
"k"))),
118 mesh.time().timeName(),
123 dimensionedScalar(
"kRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
124 zeroGradientFvPatchField<scalar>::typeName),
125 y_(mesh_.thisDb().lookupObject<volScalarField>(
"yWall")),
130 mesh.time().timeName(),
132 IOobject::READ_IF_PRESENT,
133 IOobject::AUTO_WRITE),
135 dimensionedScalar(
"betaFIK", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
140 mesh.time().timeName(),
142 IOobject::READ_IF_PRESENT,
143 IOobject::AUTO_WRITE),
145 dimensionedScalar(
"betaFIOmega", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
151 omegaRes_.dimensions().reset(dimensionSet(0, 0, -2, 0, 0, 0, 0));
152 kRes_.dimensions().reset(dimensionSet(0, 2, -3, 0, 0, 0, 0));
156 omegaRes_.dimensions().reset(dimensionSet(1, -3, -2, 0, 0, 0, 0));
157 kRes_.dimensions().reset(dimensionSet(1, -1, -3, 0, 0, 0, 0));
161 label nWallFaces = 0;
164 if (
omega_.boundaryField()[patchI].type() ==
"omegaWallFunction"
165 &&
omega_.boundaryField()[patchI].size() > 0)
178 tmp<volTensorField> tgradU = fvc::grad(
U_);
179 volScalarField S2(2 * magSqr(symm(tgradU())));
180 volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
181 GPtr_.reset(
new volScalarField::Internal(
"kOmegaSST:G",
nut_ * GbyNu0));
188 const volScalarField& CDkOmega)
const
191 tmp<volScalarField> CDkOmegaPlus = max(
193 dimensionedScalar(
"1.0e-10", dimless / sqr(dimTime), 1.0e-10));
195 tmp<volScalarField> arg1 = min(
199 scalar(500) * (this->
nu()) / (sqr(
y_) *
omega_)),
203 return tanh(pow4(arg1));
209 tmp<volScalarField> arg2 = min(
212 scalar(500) * (this->
nu()) / (sqr(
y_) *
omega_)),
215 return tanh(sqr(arg2));
221 tmp<volScalarField> arg3 = min(
225 return 1 - tanh(pow4(arg3));
230 tmp<volScalarField> f23(
F2());
241 const volScalarField::Internal& GbyNu0,
242 const volScalarField::Internal& F2,
243 const volScalarField::Internal& S2)
const
252 const volScalarField::Internal& G)
const
258 const volScalarField& F1,
259 const volTensorField& gradU)
const
266 return tmp<fvScalarMatrix>(
274 return tmp<fvScalarMatrix>(
281 const volScalarField::Internal& S2,
282 const volScalarField::Internal& gamma,
283 const volScalarField::Internal& beta)
const
285 return tmp<fvScalarMatrix>(
319 word stateName = modelStates[idxI];
320 if (stateName ==
"nut")
322 modelStates[idxI] =
"omega";
323 modelStates.append(
"k");
336 const volVectorField
U =
mesh_.thisDb().lookupObject<volVectorField>(
"U");
337 tmp<volTensorField> tgradU = fvc::grad(
U);
338 volScalarField S2(2 * magSqr(symm(tgradU())));
342 nut_.correctBoundaryConditions();
359 k_.correctBoundaryConditions();
384 omega_.correctBoundaryConditions();
398 if (
omega_.boundaryField()[patchI].type() ==
"omegaWallFunction"
399 and
omega_.boundaryField()[patchI].size() > 0)
401 const UList<label>& faceCells =
mesh_.boundaryMesh()[patchI].faceCells();
423 if (
omega_.boundaryField()[patchI].type() ==
"omegaWallFunction"
424 &&
omega_.boundaryField()[patchI].size() > 0)
426 const UList<label>& faceCells =
mesh_.boundaryMesh()[patchI].faceCells();
431 omega_.boundaryFieldRef()[patchI][faceI] =
omega_[faceCells[faceI]];
487 label stateConSize = stateCon.size();
491 forAll(stateCon[idxI], idxJ)
493 word conStateName = stateCon[idxI][idxJ];
494 if (conStateName ==
"nut")
496 stateCon[idxI][idxJ] =
"omega";
497 stateCon[idxI].append(
"k");
507 forAll(stateCon[idxI], idxJ)
509 word conStateName = stateCon[idxI][idxJ];
510 if (conStateName ==
"U")
517 stateCon[idxI].append(
"U");
522 if (idxI != stateConSize - 1)
525 forAll(stateCon[idxI + 1], idxJ)
527 word conStateName = stateCon[idxI + 1][idxJ];
528 if (conStateName ==
"U")
535 stateCon[idxI + 1].append(
"U");
541 "In DAStateInfo, nut shows in the largest connectivity level! "
542 "This is not supported!")
589 if (
mesh_.thisDb().foundObject<volScalarField>(
"p"))
593 else if (
mesh_.thisDb().foundObject<volScalarField>(
"p_rgh"))
600 "Neither p nor p_rgh was found in mesh.thisDb()!"
601 "addModelResidualCon failed to setup turbulence residuals!")
611 {
"U",
"omega",
"k",
"phi"},
618 {
"U",
"omega",
"k",
"phi"},
628 {
"U",
"T", pName,
"omega",
"k",
"phi"},
629 {
"U",
"T", pName,
"omega",
"k"},
630 {
"U",
"T", pName,
"omega",
"k"}
635 {
"U",
"T", pName,
"omega",
"k",
"phi"},
636 {
"U",
"T", pName,
"omega",
"k"},
637 {
"U",
"T", pName,
"omega",
"k"}
655 dictionary dummyOptions;
656 dummyOptions.set(
"printToScreen", printToScreen);
685 word divKScheme =
"div(phi,k)";
686 word divOmegaScheme =
"div(phi,omega)";
689 label printToScreen = options.lookupOrDefault(
"printToScreen", 0);
693 isPC = options.getLabel(
"isPC");
697 divKScheme =
"div(pc)";
698 divOmegaScheme =
"div(pc)";
702 volScalarField
rho = this->
rho();
707 volScalarField::Internal divU(fvc::div(fvc::absolute(
phi_ / fvc::interpolate(
rho),
U_)));
709 tmp<volTensorField> tgradU = fvc::grad(
U_);
710 volScalarField S2(2 * magSqr(symm(tgradU())));
711 volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
712 volScalarField::Internal& G =
const_cast<volScalarField::Internal&
>(
GPtr_());
718 omega_.boundaryFieldRef().updateCoeffs();
727 volScalarField CDkOmega(
730 volScalarField
F1(this->
F1(CDkOmega));
731 volScalarField
F23(this->
F23());
735 volScalarField::Internal
gamma(this->
gamma(F1));
736 volScalarField::Internal
beta(this->
beta(F1));
739 tmp<fvScalarMatrix> omegaEqn(
754 omegaEqn.ref().relax();
755 omegaEqn.ref().boundaryManipulate(
omega_.boundaryFieldRef());
762 SolverPerformance<scalar> solverOmega =
solve(omegaEqn);
780 tmp<fvScalarMatrix> kEqn(
785 - fvm::SuSp((2.0 / 3.0) *
phase_() *
rho() * divU,
k_)
798 SolverPerformance<scalar> solverK =
solve(kEqn);
828 if (varName !=
"k" && varName !=
"omega")
831 "varName not valid. It has to be k or omega")
835 volScalarField
rho = this->
rho();
840 volScalarField::Internal divU(fvc::div(fvc::absolute(
phi_ / fvc::interpolate(
rho),
U_)));
842 tmp<volTensorField> tgradU = fvc::grad(
U_);
843 volScalarField S2(2 * magSqr(symm(tgradU())));
844 volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
845 volScalarField::Internal& G =
const_cast<volScalarField::Internal&
>(
GPtr_());
852 volScalarField CDkOmega(
855 volScalarField
F1(this->
F1(CDkOmega));
856 volScalarField
F23(this->
F23());
858 if (varName ==
"omega")
860 volScalarField::Internal
gamma(this->
gamma(F1));
861 volScalarField::Internal
beta(this->
beta(F1));
864 fvScalarMatrix omegaEqn(
885 upper = omegaEqn.upper();
886 lower = omegaEqn.lower();
888 else if (varName ==
"k")
896 - fvm::SuSp((2.0 / 3.0) *
phase_() *
rho() * divU,
k_)
903 upper = kEqn.upper();
904 lower = kEqn.lower();
914 tmp<volTensorField> tgradU = fvc::grad(
U_);
915 volScalarField S2(2 * magSqr(symm(tgradU())));
916 volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
917 volScalarField::Internal& G =
const_cast<volScalarField::Internal&
>(
GPtr_());
920 volScalarField
rho = this->
rho();
922 volScalarField CDkOmega(
925 volScalarField
F1(this->
F1(CDkOmega));
927 volScalarField::Internal P =
phase_() *
rho() *
Pk(G);
932 PoD[cellI] = P[cellI] / (
D[cellI] + P[cellI] + 1e-16);
943 tmp<volTensorField> tgradU = fvc::grad(
U_);
944 volScalarField S2(2 * magSqr(symm(tgradU())));
945 volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
946 volScalarField::Internal& G =
const_cast<volScalarField::Internal&
>(
GPtr_());
949 volScalarField
rho = this->
rho();
951 volScalarField CDkOmega(
954 volScalarField
F1(this->
F1(CDkOmega));
956 volScalarField::Internal P =
phase_() *
rho() *
Pk(G);
961 CoP[cellI] = C[cellI] / (P[cellI] + C[cellI] + 1e-16);