FreeFOAM The Cross-Platform CFD Toolkit
multiphaseMixture.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 
26 #include "multiphaseMixture.H"
28 #include <OpenFOAM/Time.H>
29 #include <OpenFOAM/subCycle.H>
30 #include <finiteVolume/fvCFD.H>
31 
32 // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
33 
34 const scalar Foam::multiphaseMixture::convertToRad =
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 void Foam::multiphaseMixture::calcAlphas()
41 {
42  scalar level = 0.0;
43  alphas_ == 0.0;
44 
45  forAllIter(PtrDictionary<phase>, phases_, iter)
46  {
47  alphas_ += level*iter();
48  level += 1.0;
49  }
50 
51  alphas_.correctBoundaryConditions();
52 }
53 
54 
55 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
56 
58 (
59  const volVectorField& U,
60  const surfaceScalarField& phi
61 )
62 :
63  transportModel(U, phi),
64  phases_(lookup("phases"), phase::iNew(U, phi)),
65  refPhase_(*phases_.lookup(word(lookup("refPhase")))),
66 
67  mesh_(U.mesh()),
68  U_(U),
69  phi_(phi),
70 
71  rhoPhi_
72  (
73  IOobject
74  (
75  "rho*phi",
76  mesh_.time().timeName(),
77  mesh_,
78  IOobject::NO_READ,
79  IOobject::NO_WRITE
80  ),
81  mesh_,
82  dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
83  ),
84 
85  alphas_
86  (
87  IOobject
88  (
89  "alphas",
90  mesh_.time().timeName(),
91  mesh_,
92  IOobject::NO_READ,
93  IOobject::AUTO_WRITE
94  ),
95  mesh_,
96  dimensionedScalar("alphas", dimless, 0.0),
97  zeroGradientFvPatchScalarField::typeName
98  ),
99 
100  sigmas_(lookup("sigmas")),
101  dimSigma_(1, 0, -2, 0, 0),
102  deltaN_
103  (
104  "deltaN",
105  1e-8/pow(average(mesh_.V()), 1.0/3.0)
106  )
107 {
108  calcAlphas();
109  alphas_.write();
110 
111  forAllIter(PtrDictionary<phase>, phases_, iter)
112  {
113  alphaTable_.add(iter());
114  }
115 }
116 
117 
118 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
119 
121 {
123 
124  tmp<volScalarField> trho = iter()*iter().rho();
125 
126  for(++iter; iter != phases_.end(); ++iter)
127  {
128  trho() += iter()*iter().rho();
129  }
130 
131  return trho;
132 }
133 
134 
136 {
138 
139  tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
140 
141  for(++iter; iter != phases_.end(); ++iter)
142  {
143  tmu() += iter()*iter().rho()*iter().nu();
144  }
145 
146  return tmu;
147 }
148 
149 
151 {
153 
155  fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
156 
157  for(++iter; iter != phases_.end(); ++iter)
158  {
159  tmuf() +=
160  fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
161  }
162 
163  return tmuf;
164 }
165 
166 
168 {
169  return mu()/rho();
170 }
171 
172 
174 {
175  return muf()/fvc::interpolate(rho());
176 }
177 
178 
181 {
183  (
185  (
186  IOobject
187  (
188  "surfaceTensionForce",
189  mesh_.time().timeName(),
190  mesh_
191  ),
192  mesh_,
194  (
195  "surfaceTensionForce",
196  dimensionSet(1, -2, -2, 0, 0),
197  0.0
198  )
199  )
200  );
201 
202  surfaceScalarField& stf = tstf();
203 
204  forAllConstIter(PtrDictionary<phase>, phases_, iter1)
205  {
206  const phase& alpha1 = iter1();
207 
209  ++iter2;
210 
211  for(; iter2 != phases_.end(); ++iter2)
212  {
213  const phase& alpha2 = iter2();
214 
215  sigmaTable::const_iterator sigma =
216  sigmas_.find(interfacePair(alpha1, alpha2));
217 
218  if (sigma == sigmas_.end())
219  {
220  FatalErrorIn("multiphaseMixture::surfaceTensionForce() const")
221  << "Cannot find interface " << interfacePair(alpha1, alpha2)
222  << " in list of sigma values"
223  << exit(FatalError);
224  }
225 
226  stf += dimensionedScalar("sigma", dimSigma_, sigma())
227  *fvc::interpolate(K(alpha1, alpha2))*
228  (
229  fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
230  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
231  );
232  }
233  }
234 
235  return tstf;
236 }
237 
238 
240 {
241  forAllIter(PtrDictionary<phase>, phases_, iter)
242  {
243  iter().correct();
244  }
245 
246  const Time& runTime = mesh_.time();
247 
248  label nAlphaSubCycles
249  (
250  readLabel
251  (
252  mesh_.solutionDict().subDict("PISO").lookup("nAlphaSubCycles")
253  )
254  );
255 
256  label nAlphaCorr
257  (
258  readLabel(mesh_.solutionDict().subDict("PISO").lookup("nAlphaCorr"))
259  );
260 
261  bool cycleAlpha
262  (
263  Switch(mesh_.solutionDict().subDict("PISO").lookup("cycleAlpha"))
264  );
265 
266  scalar cAlpha
267  (
268  readScalar(mesh_.solutionDict().subDict("PISO").lookup("cAlpha"))
269  );
270 
271 
272  volScalarField& alpha = phases_.first();
273 
274  if (nAlphaSubCycles > 1)
275  {
276  surfaceScalarField rhoPhiSum = 0.0*rhoPhi_;
277  dimensionedScalar totalDeltaT = runTime.deltaT();
278 
279  for
280  (
281  subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
282  !(++alphaSubCycle).end();
283  )
284  {
285  solveAlphas(nAlphaCorr, cycleAlpha, cAlpha);
286  rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
287  }
288 
289  rhoPhi_ = rhoPhiSum;
290  }
291  else
292  {
293  solveAlphas(nAlphaCorr, cycleAlpha, cAlpha);
294  }
295 }
296 
297 
299 {}
300 
301 
302 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
303 (
304  const volScalarField& alpha1,
305  const volScalarField& alpha2
306 ) const
307 {
308  /*
309  // Cell gradient of alpha
310  volVectorField gradAlpha =
311  alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
312 
313  // Interpolated face-gradient of alpha
314  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
315  */
316 
317  surfaceVectorField gradAlphaf =
319  - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2));
320 
321  // Face unit interface normal
322  return gradAlphaf/(mag(gradAlphaf) + deltaN_);
323 }
324 
325 
326 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
327 (
328  const volScalarField& alpha1,
329  const volScalarField& alpha2
330 ) const
331 {
332  // Face unit interface normal flux
333  return nHatfv(alpha1, alpha2) & mesh_.Sf();
334 }
335 
336 
337 // Correction for the boundary condition on the unit normal nHat on
338 // walls to produce the correct contact angle.
339 
340 // The dynamic contact angle is calculated from the component of the
341 // velocity on the direction of the interface, parallel to the wall.
342 
343 void Foam::multiphaseMixture::correctContactAngle
344 (
345  const phase& alpha1,
346  const phase& alpha2,
348 ) const
349 {
350  const volScalarField::GeometricBoundaryField& gbf
351  = refPhase_.boundaryField();
352 
353  const fvBoundaryMesh& boundary = mesh_.boundary();
354 
355  forAll(boundary, patchi)
356  {
357  if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
358  {
360  refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
361 
362  vectorField& nHatPatch = nHatb[patchi];
363 
364  vectorField AfHatPatch =
365  mesh_.Sf().boundaryField()[patchi]
366  /mesh_.magSf().boundaryField()[patchi];
367 
369  const_iterator tp =
370  acap.thetaProps().find(interfacePair(alpha1, alpha2));
371 
372  if (tp == acap.thetaProps().end())
373  {
375  (
376  "multiphaseMixture::correctContactAngle"
377  "(const phase& alpha1, const phase& alpha2, "
378  "fvPatchVectorFieldField& nHatb) const"
379  ) << "Cannot find interface " << interfacePair(alpha1, alpha2)
380  << "\n in table of theta properties for patch "
381  << acap.patch().name()
382  << exit(FatalError);
383  }
384 
385  bool matched = (tp.key().first() == alpha1.name());
386 
387  scalar theta0 = convertToRad*tp().theta0(matched);
388  scalarField theta(boundary[patchi].size(), theta0);
389 
390  scalar uTheta = tp().uTheta();
391 
392  // Calculate the dynamic contact angle if required
393  if (uTheta > SMALL)
394  {
395  scalar thetaA = convertToRad*tp().thetaA(matched);
396  scalar thetaR = convertToRad*tp().thetaR(matched);
397 
398  // Calculated the component of the velocity parallel to the wall
399  vectorField Uwall =
400  U_.boundaryField()[patchi].patchInternalField()
401  - U_.boundaryField()[patchi];
402  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
403 
404  // Find the direction of the interface parallel to the wall
405  vectorField nWall =
406  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch;
407 
408  // Normalise nWall
409  nWall /= (mag(nWall) + SMALL);
410 
411  // Calculate Uwall resolved normal to the interface parallel to
412  // the interface
413  scalarField uwall = nWall & Uwall;
414 
415  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
416  }
417 
418 
419  // Reset nHatPatch to correspond to the contact angle
420 
421  scalarField a12 = nHatPatch & AfHatPatch;
422 
423  scalarField b1 = cos(theta);
424 
425  scalarField b2(nHatPatch.size());
426 
427  forAll(b2, facei)
428  {
429  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
430  }
431 
432  scalarField det = 1.0 - a12*a12;
433 
434  scalarField a = (b1 - a12*b2)/det;
435  scalarField b = (b2 - a12*b1)/det;
436 
437  nHatPatch = a*AfHatPatch + b*nHatPatch;
438 
439  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
440  }
441  }
442 }
443 
444 
446 (
447  const phase& alpha1,
448  const phase& alpha2
449 ) const
450 {
451  tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
452 
453  correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
454 
455  // Simple expression for curvature
456  return -fvc::div(tnHatfv & mesh_.Sf());
457 }
458 
459 
462 {
463  tmp<surfaceScalarField> tnearInt
464  (
466  (
467  IOobject
468  (
469  "nearInterface",
470  mesh_.time().timeName(),
471  mesh_
472  ),
473  mesh_,
474  dimensionedScalar("nearInterface", dimless, 0.0)
475  )
476  );
477 
478  forAllConstIter(PtrDictionary<phase>, phases_, iter)
479  {
480  surfaceScalarField alphaf = fvc::interpolate(iter());
481  tnearInt() = max(tnearInt(), pos(alphaf - 0.01)*pos(0.99 - alphaf));
482  }
483 
484  return tnearInt;
485 }
486 
487 
488 void Foam::multiphaseMixture::solveAlphas
489 (
490  const label nAlphaCorr,
491  const bool cycleAlpha,
492  const scalar cAlpha
493 )
494 {
495  static label nSolves=-1;
496  nSolves++;
497 
498  word alphaScheme("div(phi,alpha)");
499  word alphacScheme("div(phirb,alpha)");
500 
502  (
504  (
505  mesh_,
506  alphaTable_,
507  phi_,
508  mesh_.divScheme(alphaScheme)
509  )
510  );
511 
512  surfaceScalarField phic = mag(phi_/mesh_.magSf());
513  phic = min(cAlpha*phic, max(phic));
514 
515  for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
516  {
517  phase* refPhasePtr = &refPhase_;
518 
519  if (cycleAlpha)
520  {
521  PtrDictionary<phase>::iterator refPhaseIter = phases_.begin();
522  for(label i=0; i<nSolves%phases_.size(); i++)
523  {
524  ++refPhaseIter;
525  }
526  refPhasePtr = &refPhaseIter();
527  }
528 
529  phase& refPhase = *refPhasePtr;
530 
531  volScalarField refPhaseNew = refPhase;
532  refPhaseNew == 1.0;
533 
534  rhoPhi_ = phi_*refPhase.rho();
535 
536  forAllIter(PtrDictionary<phase>, phases_, iter)
537  {
538  phase& alpha = iter();
539 
540  if (&alpha == &refPhase) continue;
541 
542  fvScalarMatrix alphaEqn
543  (
544  fvm::ddt(alpha)
545  + mvConvection->fvmDiv(phi_, alpha)
546  );
547 
548  forAllIter(PtrDictionary<phase>, phases_, iter2)
549  {
550  phase& alpha2 = iter2();
551 
552  if (&alpha2 == &alpha) continue;
553 
554  surfaceScalarField phir = phic*nHatf(alpha, alpha2);
555  surfaceScalarField phirb12 =
556  -fvc::flux(-phir, alpha2, alphacScheme);
557 
558  alphaEqn += fvm::div(phirb12, alpha, alphacScheme);
559  }
560 
561  alphaEqn.solve(mesh_.solver("alpha"));
562 
563  rhoPhi_ += alphaEqn.flux()*(alpha.rho() - refPhase.rho());
564 
565  Info<< alpha.name() << " volume fraction, min, max = "
566  << alpha.weightedAverage(mesh_.V()).value()
567  << ' ' << min(alpha).value()
568  << ' ' << max(alpha).value()
569  << endl;
570 
571  refPhaseNew == refPhaseNew - alpha;
572  }
573 
574  refPhase == refPhaseNew;
575 
576  Info<< refPhase.name() << " volume fraction, min, max = "
577  << refPhase.weightedAverage(mesh_.V()).value()
578  << ' ' << min(refPhase).value()
579  << ' ' << max(refPhase).value()
580  << endl;
581  }
582 
583  calcAlphas();
584 }
585 
586 
588 {
589  if (transportModel::read())
590  {
591  bool readOK = true;
592 
593  PtrList<entry> phaseData(lookup("phases"));
594  label phasei = 0;
595 
596  forAllIter(PtrDictionary<phase>, phases_, iter)
597  {
598  readOK &= iter().read(phaseData[phasei++].dict());
599  }
600 
601  lookup("sigmas") >> sigmas_;
602 
603  return readOK;
604  }
605  else
606  {
607  return false;
608  }
609 }
610 
611 
612 // ************************ vim: set sw=4 sts=4 et: ************************ //