1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "contactAngleForce.H"
27 #include "addToRunTimeSelectionTable.H"
29 #include "unitConversion.H"
30 #include "fvPatchField.H"
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 namespace regionModels
38 namespace surfaceFilmModels
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 defineTypeNameAndDebug(contactAngleForce, 0);
44 addToRunTimeSelectionTable(force, contactAngleForce, dictionary);
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 contactAngleForce::contactAngleForce
50 const surfaceFilmModel& owner,
51 const dictionary& dict
54 force(typeName, owner, dict),
55 deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
56 Ccf_(readScalar(coeffs_.lookup("Ccf"))),
57 rndGen_(label(0), -1),
60 distributionModels::distributionModel::New
62 coeffs_.subDict("contactAngleDistribution"),
69 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
71 contactAngleForce::~contactAngleForce()
75 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
77 tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
79 tmp<volVectorField> tForce
86 owner_.time().timeName(),
92 dimensionedVector("zero", dimForce/dimArea, vector::zero)
96 vectorField& force = tForce().internalField();
98 const labelUList& own = owner_.regionMesh().owner();
99 const labelUList& nbr = owner_.regionMesh().neighbour();
101 const scalarField& magSf = owner_.magSf();
103 const volScalarField& delta = owner_.delta();
104 const volScalarField& sigma = owner_.sigma();
109 pos(delta - dimensionedScalar("deltaWet", dimLength, deltaWet_))
111 volVectorField gradAlpha(fvc::grad(alpha));
114 scalarField nHits(owner_.regionMesh().nCells(), 0.0);
118 const label cellO = own[faceI];
119 const label cellN = nbr[faceI];
122 if ((delta[cellO] > deltaWet_) && (delta[cellN] < deltaWet_))
126 else if ((delta[cellO] < deltaWet_) && (delta[cellN] > deltaWet_))
133 // const scalar dx = Foam::sqrt(magSf[cellI]);
134 // bit of a cheat, but ok for regular meshes
135 const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
137 gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
138 scalar theta = cos(degToRad(distribution_->sample()));
139 force[cellI] += Ccf_*n*sigma[cellI]*(1.0 - theta)/dx;
144 forAll(delta.boundaryField(), patchI)
146 const fvPatchField<scalar>& df = delta.boundaryField()[patchI];
147 const scalarField& dx = df.patch().deltaCoeffs();
148 const labelUList& faceCells = df.patch().faceCells();
152 label cellO = faceCells[faceI];
154 if ((delta[cellO] > deltaWet_) && (df[faceI] < deltaWet_))
157 gradAlpha[cellO]/(mag(gradAlpha[cellO]) + ROOTVSMALL);
158 scalar theta = cos(degToRad(distribution_->sample()));
159 force[cellO] += Ccf_*n*sigma[cellO]*(1.0 - theta)/dx[faceI];
165 force /= (max(nHits, scalar(1.0))*magSf);
166 tForce().correctBoundaryConditions();
168 if (owner_.regionMesh().time().outputTime())
174 tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
182 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
184 } // End namespace surfaceFilmModels
185 } // End namespace regionModels
186 } // End namespace Foam
188 // ************************************************************************* //