applications/solvers/surfaceTracking/freeSurface/freeSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "freeSurface.H"
  30
  31 #include "volFields.H"
  32 #include "transformField.H"
  33
  34 #include "emptyFaPatch.H"
  35 #include "wedgeFaPatch.H"
  36 #include "wallFvPatch.H"
  37
  38 #include "EulerDdtScheme.H"
  39 #include "CrankNicholsonDdtScheme.H"
  40 #include "backwardDdtScheme.H"
  41
  42 #include "tetFemMatrices.H"
  43 #include "tetPointFields.H"
  44 #include "faceTetPolyPatch.H"
  45 #include "tetPolyPatchInterpolation.H"
  46 #include "fixedValueTetPolyPatchFields.H"
  47 #include "fixedValuePointPatchFields.H"
  48 #include "twoDPointCorrector.H"
  49
  50 #include "slipFvPatchFields.H"
  51 #include "symmetryFvPatchFields.H"
  52 #include "fixedGradientFvPatchFields.H"
  53 #include "zeroGradientCorrectedFvPatchFields.H"
  54 #include "fixedGradientCorrectedFvPatchFields.H"
  55 #include "fixedValueCorrectedFvPatchFields.H"
  56
  57 #include "primitivePatchInterpolation.H"
  58
  59 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  60
  61 namespace Foam
  62 {
  63
  64 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  65
  66 defineTypeNameAndDebug(freeSurface, 0);
  67
  68
  69 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  70
  71 void freeSurface::clearOut()
  72 {
  73     deleteDemandDrivenData(interpolatorABPtr_);
  74     deleteDemandDrivenData(interpolatorBAPtr_);
  75     deleteDemandDrivenData(controlPointsPtr_);
  76     deleteDemandDrivenData(motionPointsMaskPtr_);
  77     deleteDemandDrivenData(pointsDisplacementDirPtr_);
  78     deleteDemandDrivenData(facesDisplacementDirPtr_);
  79     deleteDemandDrivenData(totalDisplacementPtr_);
  80     deleteDemandDrivenData(aMeshPtr_);
  81     deleteDemandDrivenData(UsPtr_);
  82     deleteDemandDrivenData(phisPtr_);
  83     deleteDemandDrivenData(surfactConcPtr_);
  84     deleteDemandDrivenData(surfaceTensionPtr_);
  85     deleteDemandDrivenData(surfactantPtr_);
  86     deleteDemandDrivenData(fluidIndicatorPtr_);
  87 }
  88
  89 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  90
  91 freeSurface::freeSurface
  92 (
  93     dynamicFvMesh& m,
  94     const volScalarField& rho,
  95     volVectorField& Ub,
  96     volScalarField& Pb,
  97     const surfaceScalarField& sfPhi
  98 )
  99 :
 100     IOdictionary
 101     (
 102         IOobject
 103         (
 104             "freeSurfaceProperties",
 105             Ub.mesh().time().constant(),
 106             Ub.mesh(),
 107             IOobject::MUST_READ,
 108             IOobject::NO_WRITE
 109         )
 110     ),
 111     mesh_(m),
 112     rho_(rho),
 113     U_(Ub),
 114     p_(Pb),
 115     phi_(sfPhi),
 116     curTimeIndex_(Ub.mesh().time().timeIndex()),
 117     twoFluids_
 118     (
 119         this->lookup("twoFluids")
 120     ),
 121     normalMotionDir_
 122     (
 123         this->lookup("normalMotionDir")
 124     ),
 125     motionDir_(0, 0, 0),
 126     cleanInterface_
 127     (
 128         this->lookup("cleanInterface")
 129     ),
 130     aPatchID_(-1),
 131     bPatchID_(-1),
 132     muFluidA_
 133     (
 134         this->lookup("muFluidA")
 135     ),
 136     muFluidB_
 137     (
 138         this->lookup("muFluidB")
 139     ),
 140     rhoFluidA_
 141     (
 142         this->lookup("rhoFluidA")
 143     ),
 144     rhoFluidB_
 145     (
 146         this->lookup("rhoFluidB")
 147     ),
 148     g_(this->lookup("g")),
 149     cleanInterfaceSurfTension_
 150     (
 151         this->lookup("surfaceTension")
 152     ),
 153     fixedFreeSurfacePatches_
 154     (
 155         this->lookup("fixedFreeSurfacePatches")
 156     ),
 157     pointNormalsCorrectionPatches_
 158     (
 159         this->lookup("pointNormalsCorrectionPatches")
 160     ),
 161     nFreeSurfCorr_
 162     (
 163         readInt(this->lookup("nFreeSurfaceCorrectors"))
 164     ),
 165     smoothing_(false),
 166     interpolatorABPtr_(NULL),
 167     interpolatorBAPtr_(NULL),
 168     controlPointsPtr_(NULL),
 169     motionPointsMaskPtr_(NULL),
 170     pointsDisplacementDirPtr_(NULL),
 171     facesDisplacementDirPtr_(NULL),
 172     totalDisplacementPtr_(NULL),
 173     aMeshPtr_(NULL),
 174     UsPtr_(NULL),
 175     phisPtr_(NULL),
 176     surfactConcPtr_(NULL),
 177     surfaceTensionPtr_(NULL),
 178     surfactantPtr_(NULL),
 179     fluidIndicatorPtr_(NULL)
 180 {
 181     //Read motion direction
 182     if (!normalMotionDir_)
 183     {
 184         motionDir_ = vector(this->lookup("motionDir"));
 185         motionDir_ /= mag(motionDir_) + SMALL;
 186     }
 187
 188     // Set point normal correction patches
 189     boolList& correction = aMesh().correctPatchPointNormals();
 190
 191     forAll(pointNormalsCorrectionPatches_, patchI)
 192     {
 193         word patchName = pointNormalsCorrectionPatches_[patchI];
 194
 195         label patchID = aMesh().boundary().findPatchID(patchName);
 196
 197         if(patchID == -1)
 198         {
 199             FatalErrorIn
 200             (
 201                 "freeSurface::freeSurface(...)"
 202             )   << "Patch name for point normals correction does not exist"
 203                 << abort(FatalError);
 204         }
 205
 206         correction[patchID] = true;
 207     }
 208
 209     // Clear geometry
 210     aMesh().movePoints();
 211
 212
 213     // Detect the free surface patch
 214     forAll (mesh().boundary(), patchI)
 215     {
 216         if(mesh().boundary()[patchI].name() == "freeSurface")
 217         {
 218             aPatchID_ = patchI;
 219
 220             Info<< "Found free surface patch. ID: " << aPatchID_
 221                 << endl;
 222         }
 223     }
 224
 225     if(aPatchID() == -1)
 226     {
 227         FatalErrorIn("freeSurface::freeSurface(...)")
 228             << "Free surface patch not defined.  Please make sure that "
 229                 << " the free surface patches is named as freeSurface"
 230                 << abort(FatalError);
 231     }
 232
 233
 234     // Detect the free surface shadow patch
 235     if (twoFluids())
 236     {
 237         forAll (mesh().boundary(), patchI)
 238         {
 239             if(mesh().boundary()[patchI].name() == "freeSurfaceShadow")
 240             {
 241                 bPatchID_ = patchI;
 242
 243                 Info<< "Found free surface shadow patch. ID: "
 244                     << bPatchID_ << endl;
 245             }
 246         }
 247
 248         if(bPatchID() == -1)
 249         {
 250             FatalErrorIn("freeSurface::freeSurface(...)")
 251                 << "Free surface shadow patch not defined. "
 252                     << "Please make sure that the free surface shadow patch "
 253                     << "is named as freeSurfaceShadow."
 254                     << abort(FatalError);
 255         }
 256     }
 257
 258
 259     // Mark free surface boundary points
 260     // which belonge to processor patches
 261     forAll(aMesh().boundary(), patchI)
 262     {
 263         if
 264         (
 265             aMesh().boundary()[patchI].type()
 266          == processorFaPatch::typeName
 267         )
 268         {
 269             const labelList& patchPoints =
 270                 aMesh().boundary()[patchI].pointLabels();
 271
 272             forAll(patchPoints, pointI)
 273             {
 274                 motionPointsMask()[patchPoints[pointI]] = -1;
 275             }
 276         }
 277     }
 278
 279
 280     // Mark fixed free surface boundary points
 281     forAll(fixedFreeSurfacePatches_, patchI)
 282     {
 283         label fixedPatchID =
 284             aMesh().boundary().findPatchID
 285             (
 286                 fixedFreeSurfacePatches_[patchI]
 287             );
 288
 289         if(fixedPatchID == -1)
 290         {
 291             FatalErrorIn("freeSurface::freeSurface(...)")
 292                 << "Wrong faPatch name in the fixedFreeSurfacePatches list"
 293                     << " defined in the freeSurfaceProperties dictionary"
 294                     << abort(FatalError);
 295         }
 296
 297         const labelList& patchPoints =
 298             aMesh().boundary()[fixedPatchID].pointLabels();
 299
 300         forAll(patchPoints, pointI)
 301         {
 302             motionPointsMask()[patchPoints[pointI]] = 0;
 303         }
 304     }
 305
 306
 307     // Mark free-surface boundary point
 308     // at the axis of 2-D axisymmetic cases
 309     forAll(aMesh().boundary(), patchI)
 310     {
 311         if
 312         (
 313             aMesh().boundary()[patchI].type()
 314          == wedgeFaPatch::typeName
 315         )
 316         {
 317             const wedgeFaPatch& wedgePatch =
 318                 refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 319
 320             if(wedgePatch.axisPoint() > -1)
 321             {
 322                 motionPointsMask()[wedgePatch.axisPoint()] = 0;
 323
 324                 Info << "Axis point: "
 325                     << wedgePatch.axisPoint()
 326                     << "vector: "
 327                     << aMesh().points()[wedgePatch.axisPoint()] << endl;
 328             }
 329         }
 330     }
 331
 332
 333     // Read free-surface points total displacement if present
 334     readTotalDisplacement();
 335
 336
 337     // Read control points positions if present
 338     controlPoints();
 339
 340
 341     // Check if smoothing switch is set
 342     if (this->found("smoothing"))
 343     {
 344         smoothing_ = Switch(this->lookup("smoothing"));
 345     }
 346 }
 347
 348
 349 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
 350
 351 freeSurface::~freeSurface()
 352 {
 353     clearOut();
 354 }
 355
 356
 357 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 358
 359 void freeSurface::updateDisplacementDirections()
 360 {
 361     if(normalMotionDir())
 362     {
 363         // Update point displacement correction
 364         pointsDisplacementDir() = aMesh().pointAreaNormals();
 365
 366         // Correcte point displacement direction
 367         // at the "centerline" symmetryPlane which represents the axis
 368         // of an axisymmetric case
 369         forAll(aMesh().boundary(), patchI)
 370         {
 371             if(aMesh().boundary()[patchI].type() == wedgeFaPatch::typeName)
 372             {
 373                 const wedgeFaPatch& wedgePatch =
 374                     refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 375
 376                 vector axis = wedgePatch.axis();
 377
 378                 label centerLinePatchID =
 379                     aMesh().boundary().findPatchID("centerline");
 380
 381                 if(centerLinePatchID != -1)
 382                 {
 383                     const labelList& pointLabels =
 384                         aMesh().boundary()[centerLinePatchID].pointLabels();
 385
 386                     forAll(pointLabels, pointI)
 387                     {
 388                         vector dir =
 389                             pointsDisplacementDir()[pointLabels[pointI]];
 390
 391                         dir = (dir&axis)*axis;
 392                         dir /= mag(dir);
 393
 394                         pointsDisplacementDir()[pointLabels[pointI]] = dir;
 395                     }
 396                 }
 397                 else
 398                 {
 399                     Info << "Warning: centerline polyPatch does not exist. "
 400                         << "Free surface points displacement directions "
 401                         << "will not be corrected at the axis (centerline)"
 402                         << endl;
 403                 }
 404
 405                 break;
 406             }
 407         }
 408
 409         // Update face displacement direction
 410         facesDisplacementDir() =
 411             aMesh().faceAreaNormals().internalField();
 412
 413         // Correction of control points postion
 414         const vectorField& Cf = aMesh().areaCentres().internalField();
 415
 416         controlPoints() =
 417             facesDisplacementDir()
 418            *(facesDisplacementDir()&(controlPoints() - Cf))
 419           + Cf;
 420     }
 421 }
 422
 423
 424 bool freeSurface::predictPoints()
 425 {
 426     // Smooth interface
 427
 428     if (smoothing_)
 429     {
 430         controlPoints() = aMesh().areaCentres().internalField();
 431         movePoints(scalarField(controlPoints().size(), 0));
 432         movePoints(-fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()]);
 433     }
 434
 435     for
 436     (
 437         int freeSurfCorr=0;
 438         freeSurfCorr<nFreeSurfCorr_;
 439         freeSurfCorr++
 440     )
 441     {
 442         movePoints(phi_.boundaryField()[aPatchID()]);
 443     }
 444
 445     return true;
 446 }
 447
 448
 449 bool freeSurface::correctPoints()
 450 {
 451     for
 452     (
 453         int freeSurfCorr=0;
 454         freeSurfCorr<nFreeSurfCorr_;
 455         freeSurfCorr++
 456     )
 457     {
 458         movePoints(phi_.boundaryField()[aPatchID()]);
 459     }
 460
 461     return true;
 462 }
 463
 464
 465 bool freeSurface::movePoints(const scalarField& interfacePhi)
 466 {
 467     pointField newMeshPoints = mesh().points();
 468
 469     scalarField sweptVolCorr =
 470         interfacePhi
 471       - fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()];
 472
 473     word ddtScheme
 474     (
 475         mesh().schemesDict().ddtScheme
 476         (
 477             "ddt(" + rho().name() + ',' + U().name()+')'
 478         )
 479     );
 480
 481     if
 482     (
 483         ddtScheme
 484      == fv::CrankNicholsonDdtScheme<vector>::typeName
 485     )
 486     {
 487         sweptVolCorr *= (1.0/2.0)*DB().deltaT().value();
 488     }
 489     else if
 490     (
 491         ddtScheme
 492      == fv::EulerDdtScheme<vector>::typeName
 493     )
 494     {
 495         sweptVolCorr *= DB().deltaT().value();
 496     }
 497     else if
 498     (
 499         ddtScheme
 500      == fv::backwardDdtScheme<vector>::typeName
 501     )
 502     {
 503         if (DB().timeIndex() == 1)
 504         {
 505             sweptVolCorr *= DB().deltaT().value();
 506         }
 507         else
 508         {
 509             sweptVolCorr *= (2.0/3.0)*DB().deltaT().value();
 510         }
 511     }
 512     else
 513     {
 514         FatalErrorIn("freeSurface::movePoints()")
 515             << "Unsupported temporal differencing scheme : "
 516                 << ddtScheme
 517                 << abort(FatalError);
 518     }
 519
 520     const scalarField& Sf = aMesh().S();
 521     const vectorField& Nf = aMesh().faceAreaNormals().internalField();
 522
 523     scalarField deltaH =
 524         sweptVolCorr/(Sf*(Nf & facesDisplacementDir()));
 525
 526     pointField displacement = pointDisplacement(deltaH);
 527
 528
 529     // Move only free-surface points
 530
 531     const labelList& meshPointsA =
 532         mesh().boundaryMesh()[aPatchID()].meshPoints();
 533
 534     forAll (displacement, pointI)
 535     {
 536         newMeshPoints[meshPointsA[pointI]] += displacement[pointI];
 537     }
 538
 539     if(twoFluids_)
 540     {
 541         const labelList& meshPointsB =
 542             mesh().boundaryMesh()[bPatchID_].meshPoints();
 543
 544         pointField displacementB =
 545             interpolatorAB().pointInterpolate
 546             (
 547                 displacement
 548             );
 549
 550         forAll (displacementB, pointI)
 551         {
 552             newMeshPoints[meshPointsB[pointI]] += displacementB[pointI];
 553         }
 554     }
 555
 556     // Update total displacement field
 557
 558     if(totalDisplacementPtr_ && (curTimeIndex_ < DB().timeIndex()))
 559     {
 560         FatalErrorIn("freeSurface::movePoints()")
 561             << "Total displacement of free surface points "
 562                 << "from previous time step is not absorbed by the mesh."
 563                 << abort(FatalError);
 564     }
 565     else if (curTimeIndex_ < DB().timeIndex())
 566     {
 567         totalDisplacement() = displacement;
 568
 569         curTimeIndex_ = DB().timeIndex();
 570     }
 571     else
 572     {
 573         totalDisplacement() += displacement;
 574     }
 575
 576     twoDPointCorrector twoDPointCorr(mesh());
 577
 578     twoDPointCorr.correctPoints(newMeshPoints);
 579
 580     mesh().movePoints(newMeshPoints);
 581
 582     // faMesh motion is done automatically, using meshObject
 583     // HJ, 8/Aug/2011
 584 //     aMesh().movePoints(mesh().points());
 585
 586
 587     // Move correctedFvPatchField fvSubMeshes
 588
 589     forAll(U().boundaryField(), patchI)
 590     {
 591         if
 592         (
 593             (
 594                 U().boundaryField()[patchI].type()
 595              == fixedGradientCorrectedFvPatchField<vector>::typeName
 596             )
 597             ||
 598             (
 599                 U().boundaryField()[patchI].type()
 600              == fixedValueCorrectedFvPatchField<vector>::typeName
 601             )
 602             ||
 603             (
 604                 U().boundaryField()[patchI].type()
 605              == zeroGradientCorrectedFvPatchField<vector>::typeName
 606             )
 607         )
 608         {
 609             correctedFvPatchField<vector>& pU =
 610                 refCast<correctedFvPatchField<vector> >
 611                 (
 612                     U().boundaryField()[patchI]
 613                 );
 614
 615             pU.movePatchSubMesh();
 616         }
 617     }
 618
 619     forAll(p().boundaryField(), patchI)
 620     {
 621         if
 622         (
 623             (
 624                 p().boundaryField()[patchI].type()
 625              == fixedGradientCorrectedFvPatchField<scalar>::typeName
 626             )
 627             ||
 628             (
 629                 p().boundaryField()[patchI].type()
 630              == fixedValueCorrectedFvPatchField<scalar>::typeName
 631             )
 632             ||
 633             (
 634                 p().boundaryField()[patchI].type()
 635              == zeroGradientCorrectedFvPatchField<scalar>::typeName
 636             )
 637         )
 638         {
 639             correctedFvPatchField<scalar>& pP =
 640                 refCast<correctedFvPatchField<scalar> >
 641                 (
 642                     p().boundaryField()[patchI]
 643                 );
 644
 645             pP.movePatchSubMesh();
 646         }
 647     }
 648
 649     return true;
 650 }
 651
 652
 653 bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
 654 {
 655     if(totalDisplacementPtr_)
 656     {
 657         pointField newPoints = mesh().points();
 658
 659         const labelList& meshPointsA =
 660             mesh().boundaryMesh()[aPatchID()].meshPoints();
 661
 662         forAll (totalDisplacement(), pointI)
 663         {
 664             newPoints[meshPointsA[pointI]] -= totalDisplacement()[pointI];
 665         }
 666
 667
 668         // Check mesh motion solver type
 669         bool feMotionSolver =
 670             mesh().objectRegistry::foundObject<tetPointVectorField>
 671             (
 672                 "motionU"
 673             );
 674         bool fvMotionSolver =
 675             mesh().objectRegistry::foundObject<pointVectorField>
 676             (
 677                 "pointMotionU"
 678             );
 679
 680         if (feMotionSolver)
 681         {
 682             tetPointVectorField& motionU =
 683                 const_cast<tetPointVectorField&>
 684                 (
 685                     mesh().objectRegistry::
 686                     lookupObject<tetPointVectorField>
 687                     (
 688                         "motionU"
 689                     )
 690                 );
 691
 692             fixedValueTetPolyPatchVectorField& motionUaPatch =
 693                 refCast<fixedValueTetPolyPatchVectorField>
 694                 (
 695                     motionU.boundaryField()[aPatchID()]
 696                 );
 697
 698             tetPolyPatchInterpolation tppiAPatch
 699             (
 700                 refCast<const faceTetPolyPatch>
 701                 (
 702                     motionUaPatch.patch()
 703                 )
 704             );
 705
 706             motionUaPatch ==
 707                 tppiAPatch.pointToPointInterpolate
 708                 (
 709                     totalDisplacement()/DB().deltaT().value()
 710                 );
 711
 712             if(twoFluids_)
 713             {
 714                 const labelList& meshPointsB =
 715                     mesh().boundaryMesh()[bPatchID()].meshPoints();
 716
 717                 pointField totDisplacementB =
 718                     interpolatorAB().pointInterpolate
 719                     (
 720                         totalDisplacement()
 721                     );
 722
 723                 forAll (totDisplacementB, pointI)
 724                 {
 725                     newPoints[meshPointsB[pointI]] -=
 726                         totDisplacementB[pointI];
 727                 }
 728
 729                 fixedValueTetPolyPatchVectorField& motionUbPatch =
 730                     refCast<fixedValueTetPolyPatchVectorField>
 731                     (
 732                         motionU.boundaryField()[bPatchID()]
 733                     );
 734
 735                 tetPolyPatchInterpolation tppiBPatch
 736                 (
 737                     refCast<const faceTetPolyPatch>(motionUbPatch.patch())
 738                 );
 739
 740                 motionUbPatch ==
 741                     tppiBPatch.pointToPointInterpolate
 742                     (
 743                         totDisplacementB/DB().deltaT().value()
 744                     );
 745             }
 746         }
 747         else if (fvMotionSolver)
 748         {
 749             pointVectorField& motionU =
 750                 const_cast<pointVectorField&>
 751                 (
 752                     mesh().objectRegistry::
 753                     lookupObject<pointVectorField>
 754                     (
 755                         "pointMotionU"
 756                     )
 757                 );
 758
 759             fixedValuePointPatchVectorField& motionUaPatch =
 760                 refCast<fixedValuePointPatchVectorField>
 761                 (
 762                     motionU.boundaryField()[aPatchID()]
 763                 );
 764
 765             motionUaPatch ==
 766                 totalDisplacement()/DB().deltaT().value();
 767
 768             if(twoFluids_)
 769             {
 770                 const labelList& meshPointsB =
 771                     mesh().boundaryMesh()[bPatchID()].meshPoints();
 772
 773                 pointField totDisplacementB =
 774                     interpolatorAB().pointInterpolate
 775                     (
 776                         totalDisplacement()
 777                     );
 778
 779                 forAll (totDisplacementB, pointI)
 780                 {
 781                     newPoints[meshPointsB[pointI]] -=
 782                         totDisplacementB[pointI];
 783                 }
 784
 785                 fixedValuePointPatchVectorField& motionUbPatch =
 786                     refCast<fixedValuePointPatchVectorField>
 787                     (
 788                         motionU.boundaryField()[bPatchID()]
 789                     );
 790
 791                 motionUbPatch ==
 792                     totDisplacementB/DB().deltaT().value();
 793             }
 794         }
 795
 796         twoDPointCorrector twoDPointCorr(mesh());
 797
 798         twoDPointCorr.correctPoints(newPoints);
 799
 800         mesh().movePoints(newPoints);
 801
 802         deleteDemandDrivenData(totalDisplacementPtr_);
 803
 804         mesh().update();
 805
 806         // faMesh motion is done automatically, using meshObject
 807         // HJ, 8/Aug/2011
 808 //         aMesh().movePoints(mesh().points());
 809
 810         // Move correctedFvPatchField fvSubMeshes
 811
 812         forAll(U().boundaryField(), patchI)
 813         {
 814             if
 815             (
 816                 (
 817                     U().boundaryField()[patchI].type()
 818                  == fixedGradientCorrectedFvPatchField<vector>::typeName
 819                 )
 820                 ||
 821                 (
 822                     U().boundaryField()[patchI].type()
 823                  == fixedValueCorrectedFvPatchField<vector>::typeName
 824                 )
 825                 ||
 826                 (
 827                     U().boundaryField()[patchI].type()
 828                  == zeroGradientCorrectedFvPatchField<vector>::typeName
 829                 )
 830             )
 831             {
 832                 correctedFvPatchField<vector>& aU =
 833                     refCast<correctedFvPatchField<vector> >
 834                     (
 835                         U().boundaryField()[patchI]
 836                     );
 837
 838                 aU.movePatchSubMesh();
 839             }
 840         }
 841
 842         forAll(p().boundaryField(), patchI)
 843         {
 844             if
 845             (
 846                 (
 847                     p().boundaryField()[patchI].type()
 848                  == fixedGradientCorrectedFvPatchField<scalar>::typeName
 849                 )
 850                 ||
 851                 (
 852                     p().boundaryField()[patchI].type()
 853                  == fixedValueCorrectedFvPatchField<scalar>::typeName
 854                 )
 855                 ||
 856                 (
 857                     p().boundaryField()[patchI].type()
 858                  == zeroGradientCorrectedFvPatchField<scalar>::typeName
 859                 )
 860             )
 861             {
 862                 correctedFvPatchField<scalar>& aP =
 863                     refCast<correctedFvPatchField<scalar> >
 864                     (
 865                         p().boundaryField()[patchI]
 866                     );
 867
 868                 aP.movePatchSubMesh();
 869             }
 870         }
 871     }
 872
 873     return true;
 874 }
 875
 876
 877 bool freeSurface::moveMeshPoints()
 878 {
 879         scalarField sweptVolCorr =
 880             phi_.boundaryField()[aPatchID()]
 881           - fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()];
 882
 883         word ddtScheme
 884         (
 885             mesh().schemesDict().ddtScheme
 886             (
 887                 "ddt(" + rho().name() + ',' + U().name()+')'
 888             )
 889         );
 890
 891         if
 892         (
 893             ddtScheme
 894          == fv::CrankNicholsonDdtScheme<vector>::typeName
 895         )
 896         {
 897             sweptVolCorr *= (1.0/2.0)*DB().deltaT().value();
 898         }
 899         else if
 900         (
 901             ddtScheme
 902          == fv::EulerDdtScheme<vector>::typeName
 903         )
 904         {
 905             sweptVolCorr *= DB().deltaT().value();
 906         }
 907         else if
 908         (
 909             ddtScheme
 910          == fv::backwardDdtScheme<vector>::typeName
 911         )
 912         {
 913             sweptVolCorr *= (2.0/3.0)*DB().deltaT().value();
 914         }
 915         else
 916         {
 917             FatalErrorIn("freeSurface::movePoints()")
 918                 << "Unsupported temporal differencing scheme : "
 919                 << ddtScheme
 920                 << abort(FatalError);
 921         }
 922
 923
 924         const scalarField& Sf = aMesh().S();
 925         const vectorField& Nf = aMesh().faceAreaNormals().internalField();
 926
 927         scalarField deltaH =
 928             sweptVolCorr/(Sf*(Nf & facesDisplacementDir()));
 929
 930
 931         pointField displacement = pointDisplacement(deltaH);
 932
 933
 934         //-- Set mesh motion boundary conditions
 935
 936         tetPointVectorField& motionU =
 937             const_cast<tetPointVectorField&>
 938             (
 939                 mesh().objectRegistry::
 940                 lookupObject<tetPointVectorField>
 941                 (
 942                     "motionU"
 943                 )
 944             );
 945
 946         fixedValueTetPolyPatchVectorField& motionUaPatch =
 947             refCast<fixedValueTetPolyPatchVectorField>
 948             (
 949                 motionU.boundaryField()[aPatchID()]
 950             );
 951
 952         tetPolyPatchInterpolation tppiAPatch
 953         (
 954             refCast<const faceTetPolyPatch>
 955             (
 956                 motionUaPatch.patch()
 957             )
 958         );
 959
 960         motionUaPatch ==
 961             tppiAPatch.pointToPointInterpolate
 962             (
 963                 displacement/DB().deltaT().value()
 964             );
 965
 966         if (twoFluids())
 967         {
 968             fixedValueTetPolyPatchVectorField& motionUbPatch =
 969                 refCast<fixedValueTetPolyPatchVectorField>
 970                 (
 971                     motionU.boundaryField()[bPatchID()]
 972                 );
 973
 974             tetPolyPatchInterpolation tppiBPatch
 975             (
 976                 refCast<const faceTetPolyPatch>(motionUbPatch.patch())
 977             );
 978
 979             motionUbPatch ==
 980                 tppiBPatch.pointToPointInterpolate
 981                 (
 982                     interpolatorAB().pointInterpolate
 983                     (
 984                         displacement/DB().deltaT().value()
 985                     )
 986                 );
 987         }
 988
 989         mesh().update();
 990
 991         // faMesh motion is done automatically, using meshObject
 992         // HJ, 8/Aug/2011
 993 //         aMesh().movePoints(mesh().points());
 994
 995
 996         // Move correctedFvPatchField fvSubMeshes
 997
 998         forAll(U().boundaryField(), patchI)
 999         {
1000             if
1001             (
1002                 (
1003                     U().boundaryField()[patchI].type()
1004                  == fixedGradientCorrectedFvPatchField<vector>::typeName
1005                 )
1006                 ||
1007                 (
1008                     U().boundaryField()[patchI].type()
1009                  == fixedValueCorrectedFvPatchField<vector>::typeName
1010                 )
1011                 ||
1012                 (
1013                     U().boundaryField()[patchI].type()
1014                  == zeroGradientCorrectedFvPatchField<vector>::typeName
1015                 )
1016             )
1017             {
1018                 correctedFvPatchField<vector>& aU =
1019                     refCast<correctedFvPatchField<vector> >
1020                     (
1021                         U().boundaryField()[patchI]
1022                     );
1023
1024                 aU.movePatchSubMesh();
1025             }
1026         }
1027
1028         forAll(p().boundaryField(), patchI)
1029         {
1030             if
1031             (
1032                 (
1033                     p().boundaryField()[patchI].type()
1034                  == fixedGradientCorrectedFvPatchField<scalar>::typeName
1035                 )
1036                 ||
1037                 (
1038                     p().boundaryField()[patchI].type()
1039                  == fixedValueCorrectedFvPatchField<scalar>::typeName
1040                 )
1041                 ||
1042                 (
1043                     p().boundaryField()[patchI].type()
1044                  == zeroGradientCorrectedFvPatchField<scalar>::typeName
1045                 )
1046             )
1047             {
1048                 correctedFvPatchField<scalar>& aP =
1049                     refCast<correctedFvPatchField<scalar> >
1050                     (
1051                         p().boundaryField()[patchI]
1052                     );
1053
1054                 aP.movePatchSubMesh();
1055             }
1056         }
1057
1058     return true;
1059 }
1060
1061
1062 void freeSurface::updateBoundaryConditions()
1063 {
1064     updateVelocity();
1065     updateSurfactantConcentration();
1066     updatePressure();
1067 }
1068
1069
1070 void freeSurface::updateVelocity()
1071 {
1072     if(twoFluids())
1073     {
1074         vectorField nA = mesh().boundary()[aPatchID()].nf();
1075
1076         vectorField nB = mesh().boundary()[bPatchID()].nf();
1077
1078         scalarField DnB = interpolatorBA().faceInterpolate
1079         (
1080             mesh().boundary()[bPatchID()].deltaCoeffs()
1081         );
1082
1083         scalarField DnA = mesh().boundary()[aPatchID()].deltaCoeffs();
1084
1085
1086         vectorField UtPA =
1087             U().boundaryField()[aPatchID()].patchInternalField();
1088
1089         if
1090         (
1091             U().boundaryField()[aPatchID()].type()
1092          == fixedGradientCorrectedFvPatchField<vector>::typeName
1093         )
1094         {
1095             fixedGradientCorrectedFvPatchField<vector>& aU =
1096                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1097                 (
1098                     U().boundaryField()[aPatchID()]
1099                 );
1100
1101             UtPA += aU.corrVecGrad();
1102         }
1103
1104         UtPA -= nA*(nA & UtPA);
1105
1106
1107         vectorField UtPB = interpolatorBA().faceInterpolate
1108         (
1109             U().boundaryField()[bPatchID()].patchInternalField()
1110         );
1111
1112         if
1113         (
1114             U().boundaryField()[bPatchID()].type()
1115          == fixedValueCorrectedFvPatchField<vector>::typeName
1116         )
1117         {
1118             fixedValueCorrectedFvPatchField<vector>& bU =
1119                 refCast<fixedValueCorrectedFvPatchField<vector> >
1120                 (
1121                     U().boundaryField()[bPatchID()]
1122                 );
1123
1124             UtPB += interpolatorBA().faceInterpolate(bU.corrVecGrad());
1125         }
1126
1127         UtPB -= nA*(nA & UtPB);
1128
1129         vectorField UtFs = muFluidA().value()*DnA*UtPA
1130           + muFluidB().value()*DnB*UtPB;
1131
1132         vectorField UnFs =
1133             nA*phi_.boundaryField()[aPatchID()]
1134            /mesh().boundary()[aPatchID()].magSf();
1135
1136         Us().internalField() += UnFs - nA*(nA&Us().internalField());
1137         correctUsBoundaryConditions();
1138
1139         UtFs -= (muFluidA().value() - muFluidB().value())*
1140             (fac::grad(Us())&aMesh().faceAreaNormals())().internalField();
1141
1142
1143         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
1144
1145         if(!cleanInterface())
1146         {
1147             tangentialSurfaceTensionForce =
1148                 surfaceTensionGrad()().internalField();
1149         }
1150         else
1151         {
1152             vectorField surfaceTensionForce =
1153                 cleanInterfaceSurfTension().value()
1154                *fac::edgeIntegrate
1155                 (
1156                     aMesh().Le()*aMesh().edgeLengthCorrection()
1157                 )().internalField();
1158
1159             tangentialSurfaceTensionForce =
1160                 surfaceTensionForce
1161               - cleanInterfaceSurfTension().value()
1162                *aMesh().faceCurvatures().internalField()*nA;
1163         }
1164
1165         UtFs += tangentialSurfaceTensionForce;
1166
1167         UtFs /= muFluidA().value()*DnA + muFluidB().value()*DnB + VSMALL;
1168
1169         Us().internalField() = UnFs + UtFs;
1170         correctUsBoundaryConditions();
1171
1172         // Store old-time velocity field U()
1173         U().oldTime();
1174
1175         U().boundaryField()[bPatchID()] ==
1176             interpolatorAB().faceInterpolate(UtFs)
1177           + nB*fvc::meshPhi(rho(),U())().boundaryField()[bPatchID()]/
1178             mesh().boundary()[bPatchID()].magSf();
1179
1180         if
1181         (
1182             p().boundaryField()[bPatchID()].type()
1183          == fixedGradientFvPatchField<scalar>::typeName
1184         )
1185         {
1186             fixedGradientFvPatchField<scalar>& pB =
1187                 refCast<fixedGradientFvPatchField<scalar> >
1188                 (
1189                     p().boundaryField()[bPatchID()]
1190                 );
1191
1192             pB.gradient() =
1193                - rhoFluidB().value()
1194                 *(
1195                      nB&fvc::ddt(U())().boundaryField()[bPatchID()]
1196                  );
1197         }
1198
1199
1200         // Update fixedGradient boundary condition on patch A
1201
1202         vectorField nGradU =
1203             muFluidB().value()*(UtPB - UtFs)*DnA
1204           + tangentialSurfaceTensionForce
1205           - muFluidA().value()*nA*fac::div(Us())().internalField()
1206           + (muFluidB().value() - muFluidA().value())
1207            *(fac::grad(Us())().internalField()&nA);
1208
1209         nGradU /= muFluidA().value() + VSMALL;
1210
1211         if
1212         (
1213             U().boundaryField()[aPatchID()].type()
1214          == fixedGradientCorrectedFvPatchField<vector>::typeName
1215         )
1216         {
1217             fixedGradientCorrectedFvPatchField<vector>& aU =
1218                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1219                 (
1220                     U().boundaryField()[aPatchID()]
1221                 );
1222
1223             aU.gradient() = nGradU;
1224         }
1225         else if
1226         (
1227             U().boundaryField()[aPatchID()].type()
1228          == fixedGradientFvPatchField<vector>::typeName
1229         )
1230         {
1231             fixedGradientFvPatchField<vector>& aU =
1232                 refCast<fixedGradientFvPatchField<vector> >
1233                 (
1234                     U().boundaryField()[aPatchID()]
1235                 );
1236
1237             aU.gradient() = nGradU;
1238         }
1239         else
1240         {
1241             FatalErrorIn("freeSurface::updateVelocity()")
1242                 << "Bounary condition on " << U().name()
1243                     <<  " for freeSurface patch is "
1244                     << U().boundaryField()[aPatchID()].type()
1245                     << ", instead "
1246                     << fixedGradientCorrectedFvPatchField<vector>::typeName
1247                     << " or "
1248                     << fixedGradientFvPatchField<vector>::typeName
1249                     << abort(FatalError);
1250         }
1251     }
1252     else
1253     {
1254         vectorField nA = aMesh().faceAreaNormals().internalField();
1255
1256         vectorField UnFs =
1257             nA*phi_.boundaryField()[aPatchID()]
1258            /mesh().boundary()[aPatchID()].magSf();
1259
1260         // Correct normal component of free-surface velocity
1261         Us().internalField() += UnFs - nA*(nA&Us().internalField());
1262         correctUsBoundaryConditions();
1263
1264         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
1265
1266         if(!cleanInterface())
1267         {
1268             tangentialSurfaceTensionForce =
1269                 surfaceTensionGrad()().internalField();
1270         }
1271         else
1272         {
1273             vectorField surfaceTensionForce =
1274                 cleanInterfaceSurfTension().value()
1275                *fac::edgeIntegrate
1276                 (
1277                     aMesh().Le()*aMesh().edgeLengthCorrection()
1278                 )().internalField();
1279
1280             tangentialSurfaceTensionForce =
1281                 surfaceTensionForce
1282               - cleanInterfaceSurfTension().value()
1283                *aMesh().faceCurvatures().internalField()*nA;
1284
1285             if (muFluidA().value() < SMALL)
1286             {
1287                 tangentialSurfaceTensionForce = vector::zero;
1288             }
1289         }
1290
1291         vectorField tnGradU =
1292             tangentialSurfaceTensionForce/(muFluidA().value() + VSMALL)
1293           - (fac::grad(Us())&aMesh().faceAreaNormals())().internalField();
1294
1295         vectorField UtPA =
1296             U().boundaryField()[aPatchID()].patchInternalField();
1297         UtPA -= nA*(nA & UtPA);
1298
1299         scalarField DnA = mesh().boundary()[aPatchID()].deltaCoeffs();
1300
1301         vectorField UtFs = UtPA + tnGradU/DnA;
1302
1303         Us().internalField() = UtFs + UnFs;
1304         correctUsBoundaryConditions();
1305
1306         vectorField nGradU =
1307             tangentialSurfaceTensionForce/(muFluidA().value() + VSMALL)
1308           - nA*fac::div(Us())().internalField()
1309           - (fac::grad(Us())().internalField()&nA);
1310
1311         if
1312         (
1313             U().boundaryField()[aPatchID()].type()
1314          == fixedGradientCorrectedFvPatchField<vector>::typeName
1315         )
1316         {
1317             fixedGradientCorrectedFvPatchField<vector>& aU =
1318                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1319                 (
1320                     U().boundaryField()[aPatchID()]
1321                 );
1322
1323             aU.gradient() = nGradU;
1324         }
1325         else if
1326         (
1327             U().boundaryField()[aPatchID()].type()
1328          == fixedGradientFvPatchField<vector>::typeName
1329         )
1330         {
1331             fixedGradientFvPatchField<vector>& aU =
1332                 refCast<fixedGradientFvPatchField<vector> >
1333                 (
1334                     U().boundaryField()[aPatchID()]
1335                 );
1336
1337             aU.gradient() = nGradU;
1338         }
1339         else
1340         {
1341             FatalErrorIn("freeSurface::updateVelocity()")
1342                 << "Bounary condition on " << U().name()
1343                     <<  " for freeSurface patch is "
1344                     << U().boundaryField()[aPatchID()].type()
1345                     << ", instead "
1346                     << fixedGradientCorrectedFvPatchField<vector>::typeName
1347                     << " or "
1348                     << fixedGradientFvPatchField<vector>::typeName
1349                     << abort(FatalError);
1350         }
1351     }
1352 }
1353
1354
1355 void freeSurface::updatePressure()
1356 {
1357     // Correct pressure boundary condition at the free-surface
1358
1359     vectorField nA = mesh().boundary()[aPatchID()].nf();
1360
1361     if(twoFluids())
1362     {
1363         scalarField pA =
1364             interpolatorBA().faceInterpolate
1365             (
1366                 p().boundaryField()[bPatchID()]
1367             );
1368
1369         const scalarField& K = aMesh().faceCurvatures().internalField();
1370
1371         Info << "Free surface curvature: min = " << gMin(K)
1372             << ", max = " << gMax(K)
1373             << ", average = " << gAverage(K) << endl << flush;
1374
1375         if(cleanInterface())
1376         {
1377 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
1378             pA -= cleanInterfaceSurfTension().value()*K;
1379         }
1380         else
1381         {
1382             scalarField surfTensionK =
1383                 surfaceTension().internalField()*K;
1384
1385             pA -= surfTensionK - gAverage(surfTensionK);
1386         }
1387
1388         pA -= 2.0*(muFluidA().value() - muFluidB().value())
1389             *fac::div(Us())().internalField();
1390
1391 //         vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1392         vector R0 = vector::zero;
1393
1394         pA -= (rhoFluidA().value() - rhoFluidB().value())*
1395             (
1396                 g_.value()
1397               & (
1398                     mesh().C().boundaryField()[aPatchID()]
1399                   - R0
1400                 )
1401             );
1402
1403         p().boundaryField()[aPatchID()] == pA;
1404     }
1405     else
1406     {
1407 //         vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1408         vector R0 = vector::zero;
1409
1410         scalarField pA =
1411           - rhoFluidA().value()*
1412             (
1413                 g_.value()
1414               & (
1415                     mesh().C().boundaryField()[aPatchID()]
1416                   - R0
1417                 )
1418             );
1419
1420         const scalarField& K = aMesh().faceCurvatures().internalField();
1421
1422         Info << "Free surface curvature: min = " << gMin(K)
1423             << ", max = " << gMax(K) << ", average = " << gAverage(K)
1424             << endl;
1425
1426         if(cleanInterface())
1427         {
1428 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
1429             pA -= cleanInterfaceSurfTension().value()*K;
1430         }
1431         else
1432         {
1433             scalarField surfTensionK =
1434                 surfaceTension().internalField()*K;
1435
1436             pA -= surfTensionK - gAverage(surfTensionK);
1437         }
1438
1439         pA -= 2.0*muFluidA().value()*fac::div(Us())().internalField();
1440
1441         p().boundaryField()[aPatchID()] == pA;
1442     }
1443
1444
1445     // Set modified pressure at patches with fixed apsolute
1446     // pressure
1447
1448 //     vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1449     vector R0 = vector::zero;
1450
1451     for (int patchI=0; patchI < p().boundaryField().size(); patchI++)
1452     {
1453         if
1454         (
1455             p().boundaryField()[patchI].type()
1456          == fixedValueFvPatchScalarField::typeName
1457         )
1458         {
1459             if (patchI != aPatchID())
1460             {
1461                 p().boundaryField()[patchI] ==
1462                   - rho().boundaryField()[patchI]
1463                    *(g_.value()&(mesh().C().boundaryField()[patchI] - R0));
1464             }
1465         }
1466     }
1467 }
1468
1469
1470 void freeSurface::updateSurfaceFlux()
1471 {
1472     Phis() = fac::interpolate(Us()) & aMesh().Le();
1473 }
1474
1475
1476 void freeSurface::updateSurfactantConcentration()
1477 {
1478     if(!cleanInterface())
1479     {
1480         Info << "Correct surfactant concentration" << endl << flush;
1481
1482         updateSurfaceFlux();
1483
1484         // Crate and solve the surfactanta transport equation
1485         faScalarMatrix CsEqn
1486         (
1487             fam::ddt(surfactantConcentration())
1488           + fam::div(Phis(), surfactantConcentration())
1489           - fam::laplacian
1490             (
1491                 surfactant().surfactDiffusion(),
1492                 surfactantConcentration()
1493             )
1494         );
1495
1496
1497         if(surfactant().soluble())
1498         {
1499             const scalarField& C =
1500                 mesh().boundary()[aPatchID()]
1501                .lookupPatchField<volScalarField, scalar>("C");
1502
1503             areaScalarField Cb
1504             (
1505                 IOobject
1506                 (
1507                     "Cb",
1508                     DB().timeName(),
1509                     mesh(),
1510                     IOobject::NO_READ,
1511                     IOobject::NO_WRITE
1512                 ),
1513                 aMesh(),
1514                 dimensioned<scalar>("Cb", dimMoles/dimVolume, 0),
1515                 zeroGradientFaPatchScalarField::typeName
1516             );
1517
1518             Cb.internalField() = C;
1519             Cb.correctBoundaryConditions();
1520
1521             CsEqn +=
1522                 fam::Sp
1523                 (
1524                     surfactant().surfactAdsorptionCoeff()*Cb
1525                   + surfactant().surfactAdsorptionCoeff()
1526                    *surfactant().surfactDesorptionCoeff(),
1527                     surfactantConcentration()
1528                 )
1529               - surfactant().surfactAdsorptionCoeff()
1530                *Cb*surfactant().surfactSaturatedConc();
1531         }
1532
1533         CsEqn.solve();
1534
1535         Info << "Correct surface tension" << endl;
1536
1537         surfaceTension() =
1538             cleanInterfaceSurfTension()
1539           + surfactant().surfactR()
1540            *surfactant().surfactT()
1541            *surfactant().surfactSaturatedConc()
1542            *log(1.0 - surfactantConcentration()
1543            /surfactant().surfactSaturatedConc());
1544
1545         if(neg(min(surfaceTension().internalField())))
1546         {
1547             FatalErrorIn
1548             (
1549                 "void freeSurface::correctSurfactantConcentration()"
1550             )   << "Surface tension is negative"
1551                     << abort(FatalError);
1552         }
1553     }
1554 }
1555
1556
1557 void freeSurface::correctUsBoundaryConditions()
1558 {
1559     forAll(Us().boundaryField(), patchI)
1560     {
1561         if
1562         (
1563             Us().boundaryField()[patchI].type()
1564          == calculatedFaPatchVectorField::typeName
1565         )
1566         {
1567             vectorField& pUs = Us().boundaryField()[patchI];
1568
1569             pUs = Us().boundaryField()[patchI].patchInternalField();
1570
1571             label ngbPolyPatchID =
1572                 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1573
1574             if(ngbPolyPatchID != -1)
1575             {
1576                 if
1577                 (
1578                     (
1579                         isA<slipFvPatchVectorField>
1580                         (
1581                             U().boundaryField()[ngbPolyPatchID]
1582                         )
1583                     )
1584                  ||
1585                     (
1586                         isA<symmetryFvPatchVectorField>
1587                         (
1588                             U().boundaryField()[ngbPolyPatchID]
1589                         )
1590                     )
1591                 )
1592                 {
1593                     vectorField N =
1594                         aMesh().boundary()[patchI].ngbPolyPatchFaceNormals();
1595
1596                     pUs -= N*(N&pUs);
1597                 }
1598             }
1599         }
1600     }
1601
1602     Us().correctBoundaryConditions();
1603 }
1604
1605
1606 vector freeSurface::totalPressureForce() const
1607 {
1608     const scalarField& S = aMesh().S();
1609
1610     const vectorField& n = aMesh().faceAreaNormals().internalField();
1611
1612     const scalarField& P = p().boundaryField()[aPatchID()];
1613
1614     vectorField pressureForces = S*P*n;
1615
1616     return gSum(pressureForces);
1617 }
1618
1619
1620 vector freeSurface::totalViscousForce() const
1621 {
1622     const scalarField& S = aMesh().S();
1623     const vectorField& n = aMesh().faceAreaNormals().internalField();
1624
1625     vectorField nGradU =
1626         U().boundaryField()[aPatchID()].snGrad();
1627
1628     vectorField viscousForces =
1629       - muFluidA().value()*S
1630        *(
1631             nGradU
1632           + (fac::grad(Us())().internalField()&n)
1633           - (n*fac::div(Us())().internalField())
1634         );
1635
1636     return gSum(viscousForces);
1637 }
1638
1639
1640 vector freeSurface::totalSurfaceTensionForce() const
1641 {
1642     const scalarField& S = aMesh().S();
1643
1644     const vectorField& n = aMesh().faceAreaNormals().internalField();
1645
1646     const scalarField& K = aMesh().faceCurvatures().internalField();
1647
1648     vectorField surfTensionForces(n.size(), vector::zero);
1649
1650     if(cleanInterface())
1651     {
1652         surfTensionForces =
1653             S*cleanInterfaceSurfTension().value()
1654            *fac::edgeIntegrate
1655             (
1656                 aMesh().Le()*aMesh().edgeLengthCorrection()
1657             )().internalField();
1658     }
1659     else
1660     {
1661         surfTensionForces *= surfaceTension().internalField()*K;
1662     }
1663
1664     return gSum(surfTensionForces);
1665 }
1666
1667
1668 void freeSurface::initializeControlPointsPosition()
1669 {
1670     scalarField deltaH = scalarField(aMesh().nFaces(), 0.0);
1671
1672     pointField displacement = pointDisplacement(deltaH);
1673
1674     const faceList& faces = aMesh().faces();
1675     const pointField& points = aMesh().points();
1676
1677
1678     pointField newPoints = points + displacement;
1679
1680     scalarField sweptVol(faces.size(), 0.0);
1681
1682     forAll(faces, faceI)
1683     {
1684         sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints);
1685     }
1686
1687     vectorField faceArea(faces.size(), vector::zero);
1688
1689     forAll (faceArea, faceI)
1690     {
1691         faceArea[faceI] = faces[faceI].normal(newPoints);
1692     }
1693
1694     forAll(deltaH, faceI)
1695     {
1696         deltaH[faceI] = sweptVol[faceI]/
1697             (faceArea[faceI] & facesDisplacementDir()[faceI]);
1698     }
1699
1700     displacement = pointDisplacement(deltaH);
1701 }
1702
1703
1704 scalar freeSurface::maxCourantNumber()
1705 {
1706     scalar CoNum = 0;
1707
1708     if(cleanInterface())
1709     {
1710         const scalarField& dE =aMesh().lPN();
1711
1712         CoNum = gMax
1713         (
1714             DB().deltaT().value()/
1715             sqrt
1716             (
1717                 rhoFluidA().value()*dE*dE*dE/
1718                 2.0/M_PI/(cleanInterfaceSurfTension().value() + SMALL)
1719             )
1720         );
1721     }
1722     else
1723     {
1724         scalarField sigmaE =
1725             linearEdgeInterpolate(surfaceTension())().internalField()
1726           + SMALL;
1727
1728         const scalarField& dE =aMesh().lPN();
1729
1730         CoNum = gMax
1731         (
1732             DB().deltaT().value()/
1733             sqrt
1734             (
1735                 rhoFluidA().value()*dE*dE*dE/
1736                 2.0/M_PI/sigmaE
1737             )
1738         );
1739     }
1740
1741     return CoNum;
1742 }
1743
1744
1745 void freeSurface::updateProperties()
1746 {
1747     muFluidA_ = dimensionedScalar(this->lookup("muFluidA"));
1748
1749     muFluidB_ = dimensionedScalar(this->lookup("muFluidB"));
1750
1751     rhoFluidA_ = dimensionedScalar(this->lookup("rhoFluidA"));
1752
1753     rhoFluidB_ = dimensionedScalar(this->lookup("rhoFluidB"));
1754
1755     g_ = dimensionedVector(this->lookup("g"));
1756
1757     cleanInterfaceSurfTension_ =
1758         dimensionedScalar(this->lookup("surfaceTension"));
1759 }
1760
1761
1762 void freeSurface::writeVTK() const
1763 {
1764     aMesh().patch().writeVTK
1765     (
1766         DB().timePath()/"freeSurface",
1767         aMesh().patch(),
1768         aMesh().patch().points()
1769     );
1770 }
1771
1772
1773 void freeSurface::writeVTKControlPoints()
1774 {
1775     // Write patch and points into VTK
1776     fileName name(DB().timePath()/"freeSurfaceControlPoints");
1777     OFstream mps(name + ".vtk");
1778
1779     mps << "# vtk DataFile Version 2.0" << nl
1780         << name << ".vtk" << nl
1781         << "ASCII" << nl
1782         << "DATASET POLYDATA" << nl
1783         << "POINTS " << controlPoints().size() << " float" << nl;
1784
1785     forAll(controlPoints(), pointI)
1786     {
1787         mps << controlPoints()[pointI].x() << ' '
1788             << controlPoints()[pointI].y() << ' '
1789             << controlPoints()[pointI].z() << nl;
1790     }
1791
1792     // Write vertices
1793     mps << "VERTICES " << controlPoints().size() << ' '
1794         << controlPoints().size()*2 << nl;
1795
1796     forAll(controlPoints(), pointI)
1797     {
1798         mps << 1 << ' ' << pointI << nl;
1799     }
1800 }
1801
1802
1803 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1804
1805 } // End namespace Foam
1806
1807 // ************************************************************************* //