FreeFOAM The Cross-Platform CFD Toolkit
applyBoundaryLayer.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  applyBoundaryLayer
26 
27 Description
28  Apply a simplified boundary-layer model to the velocity and
29  turbulence fields based on the 1/7th power-law.
30 
31  The uniform boundary-layer thickness is either provided via the -ybl option
32  or calculated as the average of the distance to the wall scaled with
33  the thickness coefficient supplied via the option -Cbl. If both options
34  are provided -ybl is used.
35 
36 Usage
37 
38  - applyBoundaryLayer [OPTIONS]
39 
40  @param -case <dir>\n
41  Case directory.
42 
43  @param -parallel \n
44  Run in parallel.
45 
46  @param -help \n
47  Display help message.
48 
49  @param -doc \n
50  Display Doxygen API documentation page for this application.
51 
52  @param -srcDoc \n
53  Display Doxygen source documentation page for this application.
54 
55  @param -writenut \n
56  Output turbulent viscosity.
57 
58  @param -Cbl <scalar>\n
59  Coefficient for scaling of the average distance from the wall.
60 
61  @param -ybl <scalar>\n
62  Specify the uniform boundary layer thickness.
63 
64 \*---------------------------------------------------------------------------*/
65 
66 #include <finiteVolume/fvCFD.H>
67 #include <finiteVolume/wallDist.H>
68 
69 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70 
71 int main(int argc, char *argv[])
72 {
73  argList::validOptions.insert("ybl", "scalar");
74  argList::validOptions.insert("Cbl", "scalar");
75  argList::validOptions.insert("writenut", "");
76 
77 # include <OpenFOAM/setRootCase.H>
78 
79 # include <OpenFOAM/createTime.H>
80 # include <OpenFOAM/createMesh.H>
81 
82 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
83 
84  Info<< "Reading field U\n" << endl;
86  (
87  IOobject
88  (
89  "U",
90  runTime.timeName(),
91  mesh,
92  IOobject::MUST_READ,
93  IOobject::NO_WRITE
94  ),
95  mesh
96  );
97 
98 # include <finiteVolume/createPhi.H>
99 
100  Info<< "Calculating wall distance field" << endl;
102 
103  // Set the mean boundary-layer thickness
104  dimensionedScalar ybl("ybl", dimLength, 0);
105 
106  if (args.optionFound("ybl"))
107  {
108  // If the boundary-layer thickness is provided use it
109  ybl.value() = args.optionRead<scalar>("ybl");
110  }
111  else if (args.optionFound("Cbl"))
112  {
113  // Calculate boundary layer thickness as Cbl * mean distance to wall
114  ybl.value() = gAverage(y) * args.optionRead<scalar>("Cbl");
115  }
116  else
117  {
119  << "Neither option 'ybl' or 'Cbl' have been provided to calculate"
120  " the boundary-layer thickness"
121  << exit(FatalError);
122  }
123 
124  Info<< "\nCreating boundary-layer for U of thickness "
125  << ybl.value() << " m" << nl << endl;
126 
127  // Modify velocity by applying a 1/7th power law boundary-layer
128  // u/U0 = (y/ybl)^(1/7)
129  // assumes U0 is the same as the current cell velocity
130 
131  scalar yblv = ybl.value();
132  forAll(U, celli)
133  {
134  if (y[celli] <= yblv)
135  {
136  U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
137  }
138  }
139 
140  Info<< "Writing U" << endl;
141  U.write();
142 
143  // Update/re-write phi
144  phi = fvc::interpolate(U) & mesh.Sf();
145  phi.write();
146 
147  // Set turbulence constants
148  dimensionedScalar kappa("kappa", dimless, 0.41);
149  dimensionedScalar Cmu("Cmu", dimless, 0.09);
150 
151  // Read and modify turbulence fields if present
152 
153  IOobject epsilonHeader
154  (
155  "epsilon",
156  runTime.timeName(),
157  mesh,
158  IOobject::MUST_READ
159  );
160 
161  IOobject kHeader
162  (
163  "k",
164  runTime.timeName(),
165  mesh,
166  IOobject::MUST_READ
167  );
168 
169  IOobject nuTildaHeader
170  (
171  "nuTilda",
172  runTime.timeName(),
173  mesh,
174  IOobject::MUST_READ
175  );
176 
177  // First calculate nut
178  volScalarField nut
179  (
180  "nut",
181  sqr(kappa*min(y, ybl))*::sqrt(2)*mag(dev(symm(fvc::grad(U))))
182  );
183 
184  if (args.optionFound("writenut"))
185  {
186  Info<< "Writing nut" << endl;
187  nut.write();
188  }
189 
190 
191  // Read and modify turbulence fields if present
192 
193  if (nuTildaHeader.headerOk())
194  {
195  Info<< "Reading field nuTilda\n" << endl;
196  volScalarField nuTilda(nuTildaHeader, mesh);
197  nuTilda = nut;
198  nuTilda.correctBoundaryConditions();
199 
200  Info<< "Writing nuTilda\n" << endl;
201  nuTilda.write();
202  }
203 
204  if (kHeader.headerOk() && epsilonHeader.headerOk())
205  {
206  Info<< "Reading field k\n" << endl;
207  volScalarField k(kHeader, mesh);
208 
209  Info<< "Reading field epsilon\n" << endl;
210  volScalarField epsilon(epsilonHeader, mesh);
211 
212  scalar ck0 = ::pow(Cmu.value(), 0.25)*kappa.value();
213  k = sqr(nut/(ck0*min(y, ybl)));
214  k.correctBoundaryConditions();
215 
216  scalar ce0 = ::pow(Cmu.value(), 0.75)/kappa.value();
217  epsilon = ce0*k*sqrt(k)/min(y, ybl);
218  epsilon.correctBoundaryConditions();
219 
220  Info<< "Writing k\n" << endl;
221  k.write();
222 
223  Info<< "Writing epsilon\n" << endl;
224  epsilon.write();
225  }
226 
227  Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
228  << " ClockTime = " << runTime.elapsedClockTime() << " s"
229  << nl << endl;
230 
231  Info<< "End\n" << endl;
232 
233  return 0;
234 }
235 
236 
237 // ************************ vim: set sw=4 sts=4 et: ************************ //