contactForce.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
25 
26 #include "isothermalFilm.H"
27 #include "fvcGrad.H"
29 
30 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
31 
33 Foam::solvers::isothermalFilm::contactForce(const volScalarField& sigma) const
34 {
35  tmp<volVectorField::Internal> tForce
36  (
38  (
39  typedName("contactForce"),
40  mesh,
42  )
43  );
44 
45  vectorField& force = tForce.ref();
46 
47  const pointField& points = mesh.points();
48  const labelUList& own = mesh.owner();
49  const labelUList& nbr = mesh.neighbour();
50 
51  const scalarField& V = mesh.V();
52 
53  const volVectorField gradDelta(fvc::grad(constrainedField(delta)));
54 
55  const scalar deltaWet = this->deltaWet.value();
56 
57  boolList contactCells(delta.size(), false);
58 
59  forAll(nbr, facei)
60  {
61  const label cello = own[facei];
62  const label celln = nbr[facei];
63 
64  if ((delta[cello] > deltaWet) && (delta[celln] < deltaWet))
65  {
66  contactCells[cello] = true;
67  }
68  else if ((delta[cello] < deltaWet) && (delta[celln] > deltaWet))
69  {
70  contactCells[celln] = true;
71  }
72  }
73 
74  // Filter for film wall and surface patches
75  labelHashSet wallAndSurfacePatches(wallPatchIDs);
76  wallAndSurfacePatches.insert(surfacePatchID);
77 
78  const volScalarField::Boundary& deltaBf = delta.boundaryField();
79 
81  {
82  const fvPatch& p(mesh.boundary()[patchi]);
83 
84  // For internal coupled patches
85  if (p.coupled() && !wallAndSurfacePatches.found(patchi))
86  {
87  tmp<scalarField> tdeltan = deltaBf[patchi].patchNeighbourField();
88  const scalarField& deltan = tdeltan();
89 
90  forAll(deltan, facei)
91  {
92  const label celli = p.faceCells()[facei];
93 
94  if ((delta[celli] > deltaWet) && (deltan[facei] < deltaWet))
95  {
96  contactCells[celli] = true;
97  }
98  }
99  }
100  }
101 
102  forAll(wallPatchIDs, i)
103  {
104  const label patchi = wallPatchIDs[i];
105 
106  if
107  (
108  isA<filmContactAngleFvPatchScalarField>
109  (
111  )
112  )
113  {
114  const filmContactAngleFvPatchScalarField& contactAngle
115  (
116  refCast<const filmContactAngleFvPatchScalarField>
117  (
119  )
120  );
121 
122  const vectorField np
123  (
124  gradDelta.boundaryField()[patchi]
125  /(mag(gradDelta.boundaryField()[patchi]) + rootVSmall)
126  );
127 
128  const scalarField cosTheta
129  (
130  contactAngle.cosTheta(U.boundaryField()[patchi], np)
131  );
132 
133  const fvPatch& p(mesh.boundary()[patchi]);
134 
135  forAll(p, facei)
136  {
137  const label celli = p.faceCells()[facei];
138 
139  if (contactCells[celli])
140  {
141  const vector& n = np[facei];
142 
143  // Estimate the length of the contact line across the face
144  scalar clen = 0;
145  const face& f(p.patch()[facei]);
146  forAll(f, i)
147  {
148  clen += mag
149  (
150  n ^ (points[f[f.fcIndex(i)]] - points[f[i]])
151  );
152  }
153  clen /= 2;
154 
155  // Calculate the contact force for the film wall face
156  force[celli] =
157  n*sigma[celli]*(1 - cosTheta[facei])*clen/V[celli];
158  }
159  }
160  }
161  }
162 
163  return tForce;
164 }
165 
166 
167 // ************************************************************************* //
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &, const Field< Type > &)
Return a temporary field constructed from name, mesh,.
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const Type & value() const
Return const reference to value.
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:469
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:857
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:475
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const fvMesh & mesh
Region mesh.
Definition: solver.H:101
volScalarField & p
The thermodynamic pressure field.
const volVectorField & U
Reference to the film velocity field.
tmp< volScalarField > sigma() const
Return the film surface tension coefficient field.
const volScalarField & delta
Film thickness.
label surfacePatchID
Film surface patch ID.
dimensionedScalar deltaWet
Film thickness below which the surface is considered dry.
labelList wallPatchIDs
List of film wall patch IDs.
A class for managing temporary objects.
Definition: tmp.H:55
Calculate the gradient of the given field.
label patchi
const pointField & points
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
static const zero Zero
Definition: zero.H:97
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
const dimensionSet dimForce
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const dimensionSet dimVolume
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
dimensioned< scalar > mag(const dimensioned< Type > &)
word typedName(Name name)
Return the name of the object within the given type.
Definition: typeInfo.H:176
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Field< vector > vectorField
Specialisation of Field<T> for vector.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:213
UList< label > labelUList
Definition: UList.H:65
labelList f(nPoints)