applications/solvers/surfaceTracking/bubbleInterTrackFoam/createBubble.H

   1     // Detect the space patch
   2     label spacePatchID = -1;
   3
   4     forAll (mesh.boundary(), patchI)
   5     {
   6         if (mesh.boundary()[patchI].name() == "space")
   7         {
   8             spacePatchID = patchI;
   9
  10             Info<< "Found space patch. ID: "
  11                 << spacePatchID << endl;
  12         }
  13     }
  14
  15     if (spacePatchID < 0)
  16     {
  17         FatalErrorIn(args.executable())
  18             << "Space patch not defined. Please make sure that "
  19             << "the space patch is named as space."
  20             << abort(FatalError);
  21     }
  22
  23
  24     IOdictionary mrfProperties
  25     (
  26         IOobject
  27         (
  28             "mrfProperties",
  29             runTime.constant(),
  30             mesh,
  31             IOobject::MUST_READ,
  32             IOobject::NO_WRITE
  33         )
  34     );
  35
  36     scalar lambdaFf
  37     (
  38         readScalar(mrfProperties.lookup("lambdaFf"))
  39     );
  40
  41     scalar lambdaF0
  42     (
  43         readScalar(mrfProperties.lookup("lambdaF0"))
  44     );
  45
  46
  47     // Parameters of moving reference frame
  48     IOdictionary movingReferenceFrame
  49     (
  50         IOobject
  51         (
  52             "movingReferenceFrame",
  53             runTime.timeName(),
  54             mesh,
  55             IOobject::READ_IF_PRESENT,
  56             IOobject::AUTO_WRITE
  57         )
  58     );
  59
  60
  61     // - Bubble center
  62     vector Cb;
  63     if (movingReferenceFrame.found("Cb"))
  64     {
  65         Cb = vector
  66         (
  67             movingReferenceFrame.lookup("Cb")
  68         );
  69     }
  70     else
  71     {
  72         if (interface.twoFluids())
  73         {
  74             Cb = gSum((1.0 - fluidIndicator.internalField())
  75                *mesh.C().internalField()*mesh.V())/
  76                 gSum((1.0 - fluidIndicator.internalField())*mesh.V());
  77         }
  78         else
  79         {
  80             scalar Va = gSum(mesh.V().field());
  81             vector Ra =
  82                 gSum(mesh.C().internalField()*mesh.V().field())
  83                /gSum(mesh.V().field());
  84
  85             Info << "Ra = " << Ra << endl;
  86
  87             scalar Vb =
  88               - gSum
  89                 (
  90                     mesh.Cf().boundaryField()[interface.aPatchID()]
  91                   & mesh.Sf().boundaryField()[interface.aPatchID()]
  92                 )/3;
  93
  94             scalar V = Va + Vb;
  95
  96             Cb = (V*Ra - Va*Ra)/Vb;
  97
  98             Info << "Current bubble centre: " << Cb << endl;
  99
 100             if (mesh.nGeometricD() != 3)
 101             {
 102                 FatalErrorIn(args.executable())
 103                     << "One-fluid bubble centroid calc "
 104                         << "is not implemented for 2d mesh"
 105                         << abort(FatalError);
 106             }
 107         }
 108
 109         movingReferenceFrame.add("Cb", Cb);
 110     }
 111
 112     // - Bubble centre at the start of calculation
 113     vector Cbf;
 114     if (movingReferenceFrame.found("Cbf"))
 115     {
 116         Cbf = vector
 117         (
 118             movingReferenceFrame.lookup("Cbf")
 119         );
 120     }
 121     else
 122     {
 123         Cbf = Cb;
 124
 125         movingReferenceFrame.add("Cbf", Cbf);
 126     }
 127
 128     // - Bubble position vector
 129     dimensionedVector XF ("XF", dimLength, vector::zero);
 130     if (movingReferenceFrame.found("XF"))
 131     {
 132         XF = dimensionedVector
 133         (
 134             movingReferenceFrame.lookup("XF")
 135         );
 136     }
 137     else
 138     {
 139         XF = dimensionedVector
 140         (
 141             "XF",
 142             dimLength,
 143             Cbf
 144         );
 145
 146         movingReferenceFrame.add("XF", XF);
 147     }
 148
 149
 150     // - Bubble velocity
 151     dimensionedVector UF ("UF", dimVelocity, vector::zero);
 152     if (movingReferenceFrame.found("UF"))
 153     {
 154         UF = dimensionedVector
 155         (
 156             movingReferenceFrame.lookup("UF")
 157         );
 158     }
 159     else
 160     {
 161         movingReferenceFrame.add("UF", UF);
 162     }
 163
 164
 165     // - Bubble acceleration
 166     dimensionedVector aF ("aF", dimAcceleration, vector::zero);
 167     if (movingReferenceFrame.found("aF"))
 168     {
 169         aF = dimensionedVector
 170         (
 171             movingReferenceFrame.lookup("aF")
 172         );
 173     }
 174     else
 175     {
 176         movingReferenceFrame.add("aF", aF);
 177     }
 178
 179
 180     IOdictionary bubbleProperties
 181     (
 182         IOobject
 183         (
 184             "bubbleProperties",
 185             runTime.timeName(),
 186             mesh,
 187             IOobject::READ_IF_PRESENT,
 188             IOobject::AUTO_WRITE
 189         )
 190     );
 191
 192     scalar Vb;
 193     if (bubbleProperties.found("Vb"))
 194     {
 195         Vb = readScalar
 196         (
 197             bubbleProperties.lookup("Vb")
 198         );
 199     }
 200     else
 201     {
 202         Vb = gSum((1.0 - fluidIndicator.internalField())*mesh.V());
 203
 204         if (!interface.twoFluids())
 205         {
 206             scalar Va = gSum(mesh.V().field());
 207
 208             boundBox box(mesh.C().boundaryField()[spacePatchID]);
 209             scalar r = (box.max().x() - box.min().x())/2;
 210
 211             scalar V = 4*M_PI*pow(r,3)/3;
 212
 213             Vb = V - Va;
 214
 215             if (mesh.nGeometricD() != 3)
 216             {
 217                 FatalErrorIn(args.executable())
 218                     << "One-fluid bubble centroid calc "
 219                         << "is not implemented for 2d mesh"
 220                         << abort(FatalError);
 221             }
 222         }
 223
 224         bubbleProperties.add("Vb", Vb);
 225     }
 226
 227     scalar Vbf;
 228     if (bubbleProperties.found("Vbf"))
 229     {
 230         Vbf = readScalar
 231         (
 232             bubbleProperties.lookup("Vbf")
 233         );
 234     }
 235     else
 236     {
 237         Vbf = Vb;
 238         bubbleProperties.add("Vbf", Vbf);
 239     }