Go to the documentation of this file.00001 Info<< "Reading field U\n" << endl;
00002 volVectorField U
00003 (
00004 IOobject
00005 (
00006 "U",
00007 runTime.timeName(),
00008 mesh,
00009 IOobject::MUST_READ,
00010 IOobject::AUTO_WRITE
00011 ),
00012 mesh
00013 );
00014
00015
00016 Info<< "Creating face flux\n" << endl;
00017 surfaceScalarField phi
00018 (
00019 IOobject
00020 (
00021 "phi",
00022 runTime.timeName(),
00023 mesh,
00024 IOobject::NO_READ,
00025 IOobject::NO_WRITE
00026 ),
00027 mesh,
00028 dimensionedScalar("zero", mesh.Sf().dimensions()*U.dimensions(), 0.0)
00029 );
00030
00031
00032 singlePhaseTransportModel laminarTransport(U, phi);
00033
00034 autoPtr<incompressible::RASModel> turbulence
00035 (
00036 incompressible::RASModel::New(U, phi, laminarTransport)
00037 );
00038
00039
00040 IOdictionary transportProperties
00041 (
00042 IOobject
00043 (
00044 "transportProperties",
00045 runTime.constant(),
00046 mesh,
00047 IOobject::MUST_READ,
00048 IOobject::NO_WRITE
00049 )
00050 );
00051
00052 dimensionedVector Ubar
00053 (
00054 transportProperties.lookup("Ubar")
00055 );
00056
00057 vector flowDirection = (Ubar/mag(Ubar)).value();
00058 tensor flowMask = sqr(flowDirection);
00059
00060
00061
00062
00063 scalar nWallFaces(0);
00064 vector wallNormal(vector::zero);
00065
00066 const fvPatchList& patches = mesh.boundary();
00067
00068 forAll(patches, patchi)
00069 {
00070 const fvPatch& currPatch = patches[patchi];
00071
00072 if (isA<wallFvPatch>(currPatch))
00073 {
00074 forAll(currPatch, facei)
00075 {
00076 nWallFaces++;
00077
00078 if (nWallFaces == 1)
00079 {
00080 wallNormal =
00081 - mesh.Sf().boundaryField()[patchi][facei]
00082 /mesh.magSf().boundaryField()[patchi][facei];
00083 }
00084 else if (nWallFaces == 2)
00085 {
00086 vector wallNormal2 =
00087 mesh.Sf().boundaryField()[patchi][facei]
00088 /mesh.magSf().boundaryField()[patchi][facei];
00089
00090
00091 if
00092 (
00093 mag(wallNormal & wallNormal2) > 1.01
00094 ||mag(wallNormal & wallNormal2) < 0.99
00095 )
00096 {
00097 Info<< "boundaryFoam: wall faces are not parallel"
00098 << endl
00099 << abort(FatalError);
00100 }
00101 }
00102 else
00103 {
00104 Info<< "boundaryFoam: number of wall faces > 2"
00105 << endl
00106 << abort(FatalError);
00107 }
00108 }
00109 }
00110 }
00111
00112
00113
00114 scalarField y = wallNormal & mesh.C().internalField();
00115
00116
00117 dimensionedVector gradP
00118 (
00119 "gradP",
00120 dimensionSet(0, 1, -2, 0, 0),
00121 vector(0, 0, 0)
00122 );
00123
00124