applyBoundaryLayer: Updated to support LES cases more conveniently

applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C

@@ -63,7 +63,8 @@ int main(int argc, char *argv[])
     );
     argList::addOption
     (
-        "Cbl", "scalar",
+        "Cbl",
+        "scalar",
         "boundary-layer thickness as Cbl * mean distance to wall"
     );
     argList::addBoolOption
@@ -73,6 +74,24 @@ int main(int argc, char *argv[])
     );
 
     #include "setRootCase.H"
+
+    if (!args.optionFound("ybl") && !args.optionFound("Cbl"))
+    {
+        FatalErrorIn(args.executable())
+            << "Neither option 'ybl' or 'Cbl' have been provided to calculate "
+            << "the boundary-layer thickness.\n"
+            << "Please choose either 'ybl' OR 'Cbl'."
+            << exit(FatalError);
+    }
+    else if (args.optionFound("ybl") && args.optionFound("Cbl"))
+    {
+        FatalErrorIn(args.executable())
+            << "Both 'ybl' and 'Cbl' have been provided to calculate "
+            << "the boundary-layer thickness.\n"
+            << "Please choose either 'ybl' OR 'Cbl'."
+            << exit(FatalError);
+    }
+
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
@@ -93,7 +112,6 @@ int main(int argc, char *argv[])
             U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
         }
     }
-    U.correctBoundaryConditions();
 
     Info<< "Writing U\n" << endl;
     U.write();
@@ -102,90 +120,96 @@ int main(int argc, char *argv[])
     phi = fvc::interpolate(U) & mesh.Sf();
     phi.write();
 
-    // Calculate nut
-    tmp<volScalarField> tnut = turbulence->nut();
-    volScalarField& nut = tnut();
-    volScalarField S(mag(dev(symm(fvc::grad(U)))));
-    nut = sqr(kappa*min(y, ybl))*::sqrt(2)*S;
 
-    if (args.optionFound("writenut"))
+    if (isA<incompressible::RASModel>(turbulence()))
     {
-        Info<< "Writing " << nut.name() << nl << endl;
-        nut.write();
-    }
-
-    // Create G field - used by RAS wall functions
-    volScalarField G("RASModel::G", nut*2*sqr(S));
+        // Calculate nut
+        tmp<volScalarField> tnut = turbulence->nut();
+        volScalarField& nut = tnut();
+        volScalarField S(mag(dev(symm(fvc::grad(U)))));
+        nut = sqr(kappa*min(y, ybl))*::sqrt(2)*S;
 
+        if (args.optionFound("writenut"))
+        {
+            Info<< "Writing nut" << endl;
+            nut.write();
+        }
 
-    //--- Read and modify turbulence fields
-
-    // Turbulence k
-    tmp<volScalarField> tk = turbulence->k();
-    volScalarField& k = tk();
-    scalar ck0 = pow025(Cmu)*kappa;
-    k = sqr(nut/(ck0*min(y, ybl)));
-    k.correctBoundaryConditions();
+        // Create G field - used by RAS wall functions
+        volScalarField G("RASModel::G", nut*2*sqr(S));
 
-    Info<< "Writing " << k.name() << nl << endl;
-    k.write();
 
+        //--- Read and modify turbulence fields
 
-    // Turbulence epsilon
-    tmp<volScalarField> tepsilon = turbulence->epsilon();
-    volScalarField& epsilon = tepsilon();
-    scalar ce0 = ::pow(Cmu, 0.75)/kappa;
-    epsilon = ce0*k*sqrt(k)/min(y, ybl);
-    epsilon.correctBoundaryConditions();
+        // Turbulence k
+        tmp<volScalarField> tk = turbulence->k();
+        volScalarField& k = tk();
+        scalar ck0 = pow025(Cmu)*kappa;
+        k = sqr(nut/(ck0*min(y, ybl)));
+        k.correctBoundaryConditions();
 
-    Info<< "Writing " << epsilon.name() << nl << endl;
-    epsilon.write();
+        Info<< "Writing k\n" << endl;
+        k.write();
 
 
-    // Turbulence omega
-    IOobject omegaHeader
-    (
-        "omega",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false
-    );
+        // Turbulence epsilon
 
-    if (omegaHeader.headerOk())
-    {
-        volScalarField omega(omegaHeader, mesh);
-        omega = epsilon/(Cmu*k);
-        omega.correctBoundaryConditions();
+        tmp<volScalarField> tepsilon = turbulence->epsilon();
+        volScalarField& epsilon = tepsilon();
+        scalar ce0 = ::pow(Cmu, 0.75)/kappa;
+        epsilon = ce0*k*sqrt(k)/min(y, ybl);
+        epsilon.correctBoundaryConditions();
 
-        Info<< "Writing omega\n" << endl;
-        omega.write();
-    }
+        Info<< "Writing epsilon\n" << endl;
+        epsilon.write();
 
+        // Turbulence omega
+        IOobject omegaHeader
+        (
+            "omega",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
 
-    // Turbulence nuTilda
-    IOobject nuTildaHeader
-    (
-        "nuTilda",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false
-    );
+        if (omegaHeader.headerOk())
+        {
+            volScalarField omega(omegaHeader, mesh);
+            omega =
+                epsilon
+               /(
+                   Cmu*k+dimensionedScalar("VSMALL", k.dimensions(), VSMALL)
+                );
+            omega.correctBoundaryConditions();
+
+            Info<< "Writing omega\n" << endl;
+            omega.write();
+        }
 
-    if (nuTildaHeader.headerOk())
-    {
-        volScalarField nuTilda(nuTildaHeader, mesh);
-        nuTilda = nut;
-        nuTilda.correctBoundaryConditions();
+        // Turbulence nuTilda
+        IOobject nuTildaHeader
+        (
+            "nuTilda",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (nuTildaHeader.headerOk())
+        {
+            volScalarField nuTilda(nuTildaHeader, mesh);
+            nuTilda = nut;
+            nuTilda.correctBoundaryConditions();
 
-        Info<< "Writing nuTilda\n" << endl;
-        nuTilda.write();
+            Info<< "Writing nuTilda\n" << endl;
+            nuTilda.write();
+        }
     }
 
-
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
         << "  ClockTime = " << runTime.elapsedClockTime() << " s"
         << nl << endl;

applications/utilities/preProcessing/applyBoundaryLayer/createFields.H

@@ -62,13 +62,6 @@ License
         // Calculate boundary layer thickness as Cbl * mean distance to wall
         ybl.value() = gAverage(y) * args.optionRead<scalar>("Cbl");
     }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "Neither option 'ybl' or 'Cbl' have been provided to calculate "
-            << "the boundary-layer thickness"
-            << exit(FatalError);
-    }
 
     Info<< "\nCreating boundary-layer for U of thickness "
         << ybl.value() << " m" << nl << endl;