FreeFOAM The Cross-Platform CFD Toolkit
supersonicFreestreamFvPatchVectorField.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 \*---------------------------------------------------------------------------*/
25 
29 #include <finiteVolume/volFields.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
37 
39 (
40  const fvPatch& p,
42 )
43 :
44  mixedFvPatchVectorField(p, iF),
45  UInf_(vector::zero),
46  pInf_(0),
47  TInf_(0),
48  gamma_(0)
49 {
50  refValue() = patchInternalField();
51  refGrad() = vector::zero;
52  valueFraction() = 1;
53 }
54 
55 
57 (
59  const fvPatch& p,
61  const fvPatchFieldMapper& mapper
62 )
63 :
64  mixedFvPatchVectorField(ptf, p, iF, mapper),
65  UInf_(ptf.UInf_),
66  pInf_(ptf.pInf_),
67  TInf_(ptf.TInf_),
68  gamma_(ptf.gamma_)
69 {}
70 
71 
73 (
74  const fvPatch& p,
76  const dictionary& dict
77 )
78 :
79  mixedFvPatchVectorField(p, iF),
80  UInf_(dict.lookup("UInf")),
81  pInf_(readScalar(dict.lookup("pInf"))),
82  TInf_(readScalar(dict.lookup("TInf"))),
83  gamma_(readScalar(dict.lookup("gamma")))
84 {
85  if (dict.found("value"))
86  {
88  (
89  vectorField("value", dict, p.size())
90  );
91  }
92  else
93  {
94  fvPatchField<vector>::operator=(patchInternalField());
95  }
96 
97  refValue() = *this;
98  refGrad() = vector::zero;
99  valueFraction() = 1;
100 
101  if (pInf_ < SMALL)
102  {
104  (
105  "supersonicFreestreamFvPatchVectorField::"
106  "supersonicFreestreamFvPatchVectorField"
107  "(const fvPatch&, const vectorField&, const dictionary&)",
108  dict
109  ) << " unphysical pInf specified (pInf <= 0.0)"
110  << "\n on patch " << this->patch().name()
111  << " of field " << this->dimensionedInternalField().name()
112  << " in file " << this->dimensionedInternalField().objectPath()
113  << exit(FatalIOError);
114  }
115 }
116 
117 
119 (
121 )
122 :
123  mixedFvPatchVectorField(sfspvf),
124  UInf_(sfspvf.UInf_),
125  pInf_(sfspvf.pInf_),
126  TInf_(sfspvf.TInf_),
127  gamma_(sfspvf.gamma_)
128 {}
129 
130 
132 (
135 )
136 :
137  mixedFvPatchVectorField(sfspvf, iF),
138  UInf_(sfspvf.UInf_),
139  pInf_(sfspvf.pInf_),
140  TInf_(sfspvf.TInf_),
141  gamma_(sfspvf.gamma_)
142 {}
143 
144 
145 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
146 
148 {
149  if (!size() || updated())
150  {
151  return;
152  }
153 
154  const fvPatchField<scalar>& pT =
155  patch().lookupPatchField<volScalarField, scalar>("T");
156 
157  const fvPatchField<scalar>& pp =
158  patch().lookupPatchField<volScalarField, scalar>("p");
159 
160  const fvPatchField<scalar>& ppsi =
161  patch().lookupPatchField<volScalarField, scalar>("psi");
162 
163  // Need R of the free-stream flow. Assume R is independent of location
164  // along patch so use face 0
165  scalar R = 1.0/(ppsi[0]*pT[0]);
166 
167  scalar MachInf = mag(UInf_)/sqrt(gamma_*R*TInf_);
168 
169  if (MachInf < 1.0)
170  {
172  (
173  "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
174  ) << " MachInf < 1.0, free stream must be supersonic"
175  << "\n on patch " << this->patch().name()
176  << " of field " << this->dimensionedInternalField().name()
177  << " in file " << this->dimensionedInternalField().objectPath()
178  << exit(FatalError);
179  }
180 
181  vectorField& Up = refValue();
182  valueFraction() = 1;
183 
184  // get the near patch internal cell values
185  vectorField U = patchInternalField();
186 
187 
188  // Find the component of U normal to the free-stream flow and in the
189  // plane of the free-stream flow and the patch normal
190 
191  // Direction of the free-stream flow
192  vector UInfHat = UInf_/mag(UInf_);
193 
194  // Normal to the plane defined by the free-stream and the patch normal
195  vectorField nnInfHat = UInfHat ^ patch().nf();
196 
197  //vectorField UnInf = U - nnInfHat*(nnInfHat & U);
198  //vectorField Un = UnInf - UInfHat*(UInfHat & UnInf);
199  //vectorField nHatInf =
200  // (Un/(mag(Un) + SMALL)) * sign(patch().nf() & Un);
201 
202  // Normal to the free-stream in the plane defined by the free-stream
203  // and the patch normal
204  vectorField nHatInf = nnInfHat ^ UInfHat;
205 
206  // Component of U normal to the free-stream in the plane defined by the
207  // free-stream and the patch normal
208  vectorField Un = nHatInf*(nHatInf & U);
209 
210  // The tangential component is
211  vectorField Ut = U - Un;
212 
213  // Calculate the Prandtl-Meyer function of the free-stream
214  scalar nuMachInf =
215  sqrt((gamma_ + 1)/(gamma_ - 1))
216  *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(MachInf) - 1)))
217  - atan(sqr(MachInf) - 1);
218 
219 
220  // Set the patch boundary condition based on the Mach number and direction
221  // of the flow dictated by the boundary/free-stream pressure difference
222 
223  forAll(Up, facei)
224  {
225  if (pp[facei] >= pInf_) // If outflow
226  {
227  // Assume supersonic outflow and calculate the boundary velocity
228  // according to ???
229 
230  scalar fpp =
231  sqrt(sqr(MachInf) - 1)
232  /(gamma_*sqr(MachInf))*mag(Ut[facei])*log(pp[facei]/pInf_);
233 
234  Up[facei] = Ut[facei] + fpp*nHatInf[facei];
235 
236  // Calculate the Mach number of the boundary velocity
237  scalar Mach = mag(Up[facei])/sqrt(gamma_/ppsi[facei]);
238 
239  if (Mach <= 1) // If subsonic
240  {
241  // Zero-gradient subsonic outflow
242 
243  Up[facei] = U[facei];
244  valueFraction()[facei] = 0;
245  }
246  }
247  else // if inflow
248  {
249  // Calculate the Mach number of the boundary velocity
250  // from the boundary pressure assuming constant total pressure
251  // expansion from the free-stream
252  scalar Mach =
253  sqrt
254  (
255  (2/(gamma_ - 1))*(1 + ((gamma_ - 1)/2)*sqr(MachInf))
256  *pow(pp[facei]/pInf_, (1 - gamma_)/gamma_)
257  - 2/(gamma_ - 1)
258  );
259 
260  if (Mach > 1) // If supersonic
261  {
262  // Supersonic inflow is assumed to occur according to the
263  // Prandtl-Meyer expansion process
264 
265  scalar nuMachf =
266  sqrt((gamma_ + 1)/(gamma_ - 1))
267  *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(Mach) - 1)))
268  - atan(sqr(Mach) - 1);
269 
270  scalar fpp = (nuMachInf - nuMachf)*mag(Ut[facei]);
271 
272  Up[facei] = Ut[facei] + fpp*nHatInf[facei];
273  }
274  else // If subsonic
275  {
277  (
278  "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
279  ) << "unphysical subsonic inflow has been generated"
280  << "\n on patch " << this->patch().name()
281  << " of field " << this->dimensionedInternalField().name()
282  << " in file "
283  << this->dimensionedInternalField().objectPath()
284  << exit(FatalError);
285  }
286  }
287  }
288 
290 }
291 
292 
294 {
296  os.writeKeyword("UInf") << UInf_ << token::END_STATEMENT << nl;
297  os.writeKeyword("pInf") << pInf_ << token::END_STATEMENT << nl;
298  os.writeKeyword("TInf") << TInf_ << token::END_STATEMENT << nl;
299  os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
300  writeEntry("value", os);
301 }
302 
303 
304 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
305 
307 (
310 );
311 
312 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
313 
314 } // End namespace Foam
315 
316 // ************************ vim: set sw=4 sts=4 et: ************************ //