diff --git a/ReleaseNotes b/ReleaseNotes index 2293b6209..e0c886326 100644 --- a/ReleaseNotes +++ b/ReleaseNotes @@ -1,6 +1,7 @@ # -*- mode: org; -*- # #+TITLE: *Release notes for foam-extend-4.0* +#+TITLE: *Version 4.0 - Guimaraes* #+AUTHOR: foam-extend administrators: #+AUTHOR: Hrvoje Jasak #+AUTHOR: Håkan Nilsson diff --git a/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C b/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C index d68eb5680..44945b07b 100644 --- a/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C +++ b/applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C @@ -72,6 +72,11 @@ int main(int argc, char *argv[]) U.internalField() = vector::zero; } + p.correctBoundaryConditions(); + U.correctBoundaryConditions(); + + phi == (fvc::interpolate(U) & mesh.Sf()); + # include "volContinuity.H" # include "meshCourantNo.H" diff --git a/applications/solvers/compressible/sonicDyMFoam/UEqn.H b/applications/solvers/compressible/sonicDyMFoam/UEqn.H index 54edd4e23..a6908756c 100644 --- a/applications/solvers/compressible/sonicDyMFoam/UEqn.H +++ b/applications/solvers/compressible/sonicDyMFoam/UEqn.H @@ -7,7 +7,10 @@ UEqn.relax ( - mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter())) + mesh.solutionDict().equationRelaxationFactor + ( + U.select(pimple.finalIter()) + ) ); solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index 686acebaa..bd9dab41a 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -27,8 +27,6 @@ fvm::ddt(psis, p) + fvm::div(phid, p) // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); diff --git a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H index 10f11a421..3d8b68d58 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H index c93a3b888..d9378a743 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H @@ -1,13 +1,14 @@ { // Solve the enthalpy equation in total enthalpy formulation (see K) + T.storePrevIter(); + K = 0.5*(magSqr(U)); fvScalarMatrix hEqn ( - fvm::ddt(rho, h) - + fvc::ddt(rho, K) - + fvm::div(phi, h) + fvm::div(phi, h) + + fvm::SuSp(-fvc::div(phi), h) + fvc::div(phi, K) - fvm::laplacian(turbulence->alphaEff(), h) == @@ -16,7 +17,6 @@ ); hEqn.relax(); - hEqn.solve(); // Bounding of enthalpy taken out diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H index e954c7b13..3262ebc32 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H @@ -21,21 +21,11 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -49,7 +39,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/Make/options b/applications/solvers/compressible/steadyCompressibleMRFFoam/Make/options index b61f7ce7d..efe307e09 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/Make/options +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/Make/options @@ -2,7 +2,7 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ - -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \ + -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel EXE_LIBS = \ -lfiniteVolume \ diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H index 7f769de78..f8408bf59 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/UEqn.H @@ -1,10 +1,6 @@ - // Solve the momentum equation - U.storePrevIter(); - fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H index 3aebca26b..6cf655b45 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H @@ -24,6 +24,7 @@ ); hEqn.relax(); + hEqn.solve(); // Bounding of enthalpy taken out thermo.correct(); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H index 2fa37cfcc..720d3928b 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H @@ -5,13 +5,30 @@ Urel == U; mrfZones.relativeVelocity(Urel); + // Bound the relative velocity to preserve the i to h conversion bound + // HJ, 22/Mar/2017 + volScalarField magUrel = mag(Urel); + + if (max(magUrel) > UMax) + { + volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU) + + neg(magUrel - UMax); + UrelLimiter.max(scalar(0)); + UrelLimiter.min(scalar(1)); + + Urel *= UrelLimiter; + Urel.correctBoundaryConditions(); + } + // Create rotational velocity (= omega x r) Urot == U - Urel; + T.storePrevIter(); + fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) @@ -19,7 +36,6 @@ ); iEqn.relax(); - iEqn.solve(); // Calculate enthalpy out of rothalpy diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index 61fed64e8..21d7d1df3 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -16,7 +16,7 @@ mrfZones.correctBoundaryVelocity(U); // Calculate phi for boundary conditions - phi = rhof*fvc::interpolate(U) & mesh.Sf(); + phi = rhof*(fvc::interpolate(U) & mesh.Sf()); surfaceScalarField phid2 = rhoReff/rhof*phi; @@ -29,21 +29,11 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -57,7 +47,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H index 727548f3b..6f9413765 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UrelEqn ( - fvm::ddt(rho, Urel) - + fvm::div(phi, Urel) + fvm::div(phi, Urel) + turbulence->divDevRhoReff() + rho*SRF->Su() ); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H index c37796f06..36f61d359 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H @@ -1,10 +1,11 @@ { // Solve the rothalpy equation + T.storePrevIter(); fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) @@ -12,7 +13,6 @@ ); iEqn.relax(); - iEqn.solve(); // Calculate enthalpy out of rothalpy diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H index 59e03e2c4..eed8e6966 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H @@ -20,21 +20,11 @@ p.storePrevIter(); - volScalarField divPhid - ( - "divPhid", - fvc::div(phid) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) - // Convective flux relaxation terms - + fvm::SuSp(-divPhid, p) - + divPhid*p + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUrelA, p) ); @@ -48,7 +38,8 @@ } } -# include "compressibleContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalFoam/UEqn.H b/applications/solvers/compressible/steadyUniversalFoam/UEqn.H index 10f11a421..3d8b68d58 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/UEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyUniversalFoam/hEqn.H b/applications/solvers/compressible/steadyUniversalFoam/hEqn.H index 0a5c1e1db..d024e7df1 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/hEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/hEqn.H @@ -10,8 +10,8 @@ fvScalarMatrix hEqn ( - fvm::ddt(rho, h) - + fvm::div(phi, h) + fvm::div(phi, h) + + fvm::SuSp(-fvc::div(phi), h) - fvm::laplacian(turbulence->alphaEff(), h) == fvc::div(faceU, p, "div(U,p)") @@ -22,7 +22,6 @@ ); hEqn.relax(); - hEqn.solve(); // Bounding of enthalpy taken out diff --git a/applications/solvers/compressible/steadyUniversalFoam/pEqn.H b/applications/solvers/compressible/steadyUniversalFoam/pEqn.H index 8eec3318e..667bc89bb 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/pEqn.H +++ b/applications/solvers/compressible/steadyUniversalFoam/pEqn.H @@ -24,8 +24,7 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -39,8 +38,8 @@ } } - // Use custom continuity error check -# include "universalContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H b/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H index 31278008e..5dcfe3bf7 100644 --- a/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H +++ b/applications/solvers/compressible/steadyUniversalFoam/rhoFromP.H @@ -12,4 +12,5 @@ rho = Foam::min(rho, rhoMax); rho = Foam::max(rho, rhoMin); rho.relax(); + rho.correctBoundaryConditions(); } diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H index 7f769de78..919631569 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H @@ -3,8 +3,7 @@ fvVectorMatrix UEqn ( - fvm::ddt(rho, U) - + fvm::div(phi, U) + fvm::div(phi, U) + turbulence->divDevRhoReff() ); diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H index 915787b4b..1fc70ba27 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H @@ -3,14 +3,29 @@ Urel == U; mrfZones.relativeVelocity(Urel); + // Bound the relative velocity to preserve the i to h conversion bound + // HJ, 22/Mar/2017 + volScalarField magUrel = mag(Urel); + + if (max(magUrel) > UMax) + { + volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU) + + neg(magUrel - UMax); + UrelLimiter.max(scalar(0)); + UrelLimiter.min(scalar(1)); + + Urel *= UrelLimiter; + Urel.correctBoundaryConditions(); + } + // Create rotational velocity (= omega x r) Urot == U - Urel; fvScalarMatrix iEqn ( - fvm::ddt(rho, i) - + fvm::div(phi, i) - - fvm::laplacian(turbulence->alphaEff(), i) + fvm::div(phi, i) + + fvm::SuSp(-fvc::div(phi), i) + - fvm::laplacian(turbulence->alphaEff(), i) == // Viscous heating: note sign (devRhoReff has a minus in it) - (turbulence->devRhoReff() && fvc::grad(U)) @@ -23,8 +38,8 @@ // From rothalpy, calculate enthalpy after solution of rothalpy equation h = i + 0.5*(magSqr(Urot) - magSqr(Urel)); h.correctBoundaryConditions(); + // Update thermo for new h thermo.correct(); psis = thermo.psi()/thermo.Cp()*thermo.Cv(); } - diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H index bba75ac7f..23973366c 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H @@ -33,8 +33,7 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psis, p) - + fvm::div(phid, p) + fvm::div(phid, p) + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -48,8 +47,8 @@ } } - // Use custom continuity error check -# include "universalContinuityErrs.H" + // Use incompressible continuity error check: div(rho U) = 0 +# include "continuityErrs.H" // Relax the pressure p.relax(); diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H b/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H index 31278008e..5dcfe3bf7 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/rhoFromP.H @@ -12,4 +12,5 @@ rho = Foam::min(rho, rhoMax); rho = Foam::max(rho, rhoMin); rho.relax(); + rho.correctBoundaryConditions(); } diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H b/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H deleted file mode 100644 index 723f773e1..000000000 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/universalContinuityErrs.H +++ /dev/null @@ -1,49 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Global - continuityErrs - -Description - Calculates and prints the continuity errors. - -\*---------------------------------------------------------------------------*/ - -{ - volScalarField contErr = fvc::ddt(rho) + fvc::div(phi); - - sumLocalContErr = runTime.deltaT().value()* - mag(contErr)().weightedAverage(mesh.V()).value(); - - globalContErr = runTime.deltaT().value()* - contErr.weightedAverage(mesh.V()).value(); - - cumulativeContErr += globalContErr; - - Info<< "time step continuity errors : sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} - -// ************************************************************************* // diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C index 8e96be5a4..c7f481fa4 100644 --- a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C +++ b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C @@ -94,7 +94,7 @@ int main(int argc, char *argv[]) solve ( - tpEqn + tpEqn() == - fvc::div(U) ); @@ -112,7 +112,9 @@ int main(int argc, char *argv[]) U.correctBoundaryConditions(); p.correctBoundaryConditions(); - phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource; + phi = (fvc::interpolate(U) & mesh.Sf()) + + tpEqn().flux() + + tpresSource; } // Make flux relative in rotating zones diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H index 189324c0b..a189057e0 100644 --- a/applications/solvers/coupled/MRFPorousFoam/UEqn.H +++ b/applications/solvers/coupled/MRFPorousFoam/UEqn.H @@ -1,10 +1,4 @@ { - volScalarField divPhi - ( - "divPhi", - fvc::div(phi) - ); - // Update boundary velocity for consistency with the flux mrfZones.correctBoundaryVelocity(U); @@ -15,13 +9,10 @@ + turbulence->divDevReff() ); - // Add MRF sources - mrfZones.addCoriolis(UEqn); - - // Add porous sources + // Add MRF and porous sources implicitly. HJ, 18/Nov/2017 tmp tTU; - if (addPorosity) + if (addMRF || addPorosity) { tTU = tmp ( @@ -41,6 +32,10 @@ ); volTensorField& TU = tTU(); + // Add implicit MRF source as a Hodge dual of the rotational velocity + TU += *mrfZones.omega(); + + // Add implicit resistance pZones.addResistance(UEqn, TU); trTU = inv(TU + tensor(I)*UEqn.A()); @@ -58,7 +53,7 @@ // Insert momentum equation UpEqn.insertEquation(0, UEqn); - if (addPorosity) + if (addMRF || addPorosity) { // Manually over-ride the 3x3 block to handle the off-diagonal // part of the Ap coefficient @@ -68,11 +63,43 @@ const scalarField& V = mesh.V().field(); - // Note: insertion should only happen in porous cell zones + // Note: insertion should only happen in MRF and porous cell zones // HJ, 14/Mar/2016 + // Warning. Possible problem with a zone which is both MRF and porous + // The solution is to do the loop below everywhere + // Currently only inserting zones for efficiency. HJ, 18/Nov/2017 register label cellI; + forAll (mrfZones, mrfZoneI) + { + const labelList& curZoneCells = mrfZones[mrfZoneI].zone(); + + // Loop over all cells in the zone + forAll (curZoneCells, zcI) + { + cellI = curZoneCells[zcI]; + + const scalar& cellV = V[cellI]; + + const tensor& cellTU = TUIn[cellI]; + + CoeffField::squareType& cellDD = DD[cellI]; + + cellDD(0, 0) += cellV*cellTU.xx(); + cellDD(0, 1) += cellV*cellTU.xy(); + cellDD(0, 2) += cellV*cellTU.xz(); + + cellDD(1, 0) += cellV*cellTU.yx(); + cellDD(1, 1) += cellV*cellTU.yy(); + cellDD(2, 2) += cellV*cellTU.yz(); + + cellDD(2, 0) += cellV*cellTU.zx(); + cellDD(2, 1) += cellV*cellTU.zy(); + cellDD(2, 2) += cellV*cellTU.zz(); + } + } + forAll (pZones, pZoneI) { const labelList& curZoneCells = pZones[pZoneI].zone(); diff --git a/applications/solvers/coupled/MRFPorousFoam/createMRF.H b/applications/solvers/coupled/MRFPorousFoam/createMRF.H index 161446a8e..3f1204663 100644 --- a/applications/solvers/coupled/MRFPorousFoam/createMRF.H +++ b/applications/solvers/coupled/MRFPorousFoam/createMRF.H @@ -1,2 +1,9 @@ MRFZones mrfZones(mesh); mrfZones.correctBoundaryVelocity(U); + + bool addMRF = false; + + if (!mrfZones.empty()) + { + addMRF = true; + } diff --git a/applications/solvers/coupled/MRFPorousFoam/pEqn.H b/applications/solvers/coupled/MRFPorousFoam/pEqn.H index 6b1823c06..c14b41b7e 100644 --- a/applications/solvers/coupled/MRFPorousFoam/pEqn.H +++ b/applications/solvers/coupled/MRFPorousFoam/pEqn.H @@ -2,7 +2,7 @@ tmp tpEqn; tmp tpresSource; -if (addPorosity) +if (addMRF || addPorosity) { // Collect pressure source with tensorial 1/Ap const volTensorField& rTU = trTU(); @@ -10,7 +10,7 @@ if (addPorosity) tpresSource = ( (mesh.Sf() & fvc::interpolate(rTU)) - & fvc::interpolate(fvc::grad(p)) + & fvc::interpolate(fvc::grad(p)) ); const surfaceScalarField& presSource = tpresSource(); diff --git a/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H index 93c2bda47..1f917aa99 100644 --- a/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H +++ b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H @@ -1,5 +1,6 @@ label pRefCell = 0; scalar pRefValue = 0; + setRefCell ( p, diff --git a/applications/solvers/coupled/pUCoupledFoam/UEqn.H b/applications/solvers/coupled/pUCoupledFoam/UEqn.H index 759e22b83..e1535831e 100644 --- a/applications/solvers/coupled/pUCoupledFoam/UEqn.H +++ b/applications/solvers/coupled/pUCoupledFoam/UEqn.H @@ -1,10 +1,4 @@ { - volScalarField divPhi - ( - "divPhi", - fvc::div(phi) - ); - // Momentum equation fvVectorMatrix UEqn ( diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H index 9c11a3332..32bd0e614 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H @@ -8,7 +8,10 @@ UEqn.relax ( - mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter())) + mesh.solutionDict().equationRelaxationFactor + ( + U.select(pimple.finalIter()) + ) ); solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H index 530b4ef30..261366d08 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H @@ -12,4 +12,8 @@ hEqn.solve(); thermo.correct(); + + // Recalculate density + rho = thermo.rho(); + rho.correctBoundaryConditions(); } diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H index 45d692a63..dd396282e 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H @@ -1,6 +1,4 @@ { - rho = thermo.rho(); - rUA = 1.0/UEqn.A(); U = rUA*UEqn.H(); @@ -56,12 +54,10 @@ p.relax(); } -# include "rhoEqn.H" # include "compressibleContinuityErrs.H" U -= rUA*fvc::grad(p); U.correctBoundaryConditions(); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); - dpdt = fvc::ddt(p); } diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C index 587c19dc6..dd8707474 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C @@ -79,14 +79,15 @@ int main(int argc, char *argv[]) Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl; - // Make flux absolute - phi += meshFlux; + // Make the fluxes absolute (using the ddt(rho, U) scheme) + phi += fvc::interpolate(rho)*fvc::meshPhi(rho, U); bool meshChanged = mesh.update(); # include "volContinuity.H" - mesh.setBoundaryVelocity(U); + // Make the fluxes relative (using the ddt(rho, U) scheme) + phi -= fvc::interpolate(rho)*fvc::meshPhi(rho, U); if (meshChanged) { @@ -95,14 +96,7 @@ int main(int argc, char *argv[]) rho.correctBoundaryConditions(); } - meshFlux = fvc::interpolate(rho)*fvc::meshPhi(rho, U); - - phi = fvc::interpolate(rho) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)); - - DpDt = dpdt + fvc::div(phi/fvc::interpolate(rho), p) - - fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(U))*p; - + if (meshChanged) { # include "compressibleCourantNo.H" } @@ -111,22 +105,20 @@ int main(int argc, char *argv[]) while (pimple.loop()) { # include "rhoEqn.H" +# include "hEqn.H" # include "UEqn.H" // --- PISO loop while (pimple.correct()) { # include "pEqn.H" -# include "hEqn.H" } + + turbulence->correct(); } - turbulence->correct(); - # include "logSummary.H" - rho = thermo.rho(); - runTime.write(); # include "infoDataOutput.H" diff --git a/applications/solvers/immersedBoundary/SOLVER_COOKBOOK b/applications/solvers/immersedBoundary/SOLVER_COOKBOOK deleted file mode 100644 index 762319239..000000000 --- a/applications/solvers/immersedBoundary/SOLVER_COOKBOOK +++ /dev/null @@ -1,49 +0,0 @@ -Make/options: add - - -I$(LIB_SRC)/triSurface/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I../immersedBoundary/lnInclude - - - -ltriSurface \ - -lmeshTools \ - -L$(FOAM_USER_LIBBIN) -limmersedBoundary \ - - -createFields.H: add - -# include "createIbMasks.H" - - -solver: - -1) add immersed boundary headers - -#include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" - -2) on calculation of face fluxes (or their components), mask with faceIbMask - -3) before to adjustPhi, add: - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - -4) on explicit updates, add correctBoundaryConditions(); - -eg. - - U = fvc::reconstruct(phi); - U.correctBoundaryConditions(); - -5) On reports of continuity error, add masking: - - Info<< "IB-masked continuity error = " - << mag(cellIbMask*fvc::div(phi))().weightedAverage(mesh.V()).value() - << endl; - - or use immersedBoundaryContinuityErrs.H - -5) Chenge Courant number check to be IB-sensitive, using - - immersedBoundaryCourantNo.H diff --git a/applications/solvers/immersedBoundary/icoDyMIbFoam/Make/files b/applications/solvers/immersedBoundary/icoDyMIbFoam/Make/files deleted file mode 100644 index 5785f9fea..000000000 --- a/applications/solvers/immersedBoundary/icoDyMIbFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -icoDyMIbFoam.C - -EXE = $(FOAM_APPBIN)/icoDyMIbFoam diff --git a/applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H b/applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H deleted file mode 100644 index 5af8f976c..000000000 --- a/applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H +++ /dev/null @@ -1,56 +0,0 @@ - Info<< "Reading transportProperties\n" << endl; - - IOdictionary transportProperties - ( - IOobject - ( - "transportProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ_IF_MODIFIED, - IOobject::NO_WRITE - ) - ); - - dimensionedScalar nu - ( - transportProperties.lookup("nu") - ); - - - Info<< "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - - Info<< "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - -# include "createPhi.H" - - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, pimple.dict(), pRefCell, pRefValue); - mesh.schemesDict().setFluxRequired(p.name()); diff --git a/applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C b/applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C deleted file mode 100644 index e83440a27..000000000 --- a/applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C +++ /dev/null @@ -1,152 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Application - icoDyMOversetFoam - -Description - Transient solver for incompressible, laminar flow of Newtonian fluids - with dynamic mesh and immersed boundary mesh support. - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "dynamicFvMesh.H" -#include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" -#include "pimpleControl.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ -# include "setRootCase.H" -# include "createTime.H" -# include "createDynamicFvMesh.H" - - pimpleControl pimple(mesh); - -# include "createFields.H" -# include "initContinuityErrs.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.loop()) - { - Info<< "Time = " << runTime.timeName() << nl << endl; - - // Make the fluxes absolute - fvc::makeAbsolute(phi, U); - - bool meshChanged = mesh.update(); - reduce(meshChanged, orOp()); - Info<< "Mesh update" << meshChanged << endl; -# include "createIbMasks.H" - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); - -# include "CourantNo.H" - - // Pressure-velocity corrector - while (pimple.loop()) - { - fvVectorMatrix UEqn - ( - fvm::ddt(U) - + fvm::div(phi, U) - - fvm::laplacian(nu, U) - ); - - if (pimple.momentumPredictor()) - { - solve(UEqn == -fvc::grad(p)); - } - - // --- PISO loop - while (pimple.correct()) - { - volScalarField rUA = 1.0/UEqn.A(); - - U = rUA*UEqn.H(); - // Immersed boundary update - U.correctBoundaryConditions(); - - phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf()); - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - adjustPhi(phi, U, p); - - // Non-orthogonal pressure corrector loop - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(rUA, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - pEqn.solve - ( - mesh.solutionDict().solver - ( - p.select(pimple.finalInnerIter()) - ) - ); - - if (pimple.finalNonOrthogonalIter()) - { - phi -= pEqn.flux(); - } - } - -# include "immersedBoundaryContinuityErrs.H" - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); - - U -= rUA*fvc::grad(p); - U.correctBoundaryConditions(); - } - } - - runTime.write(); - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return(0); -} - - -// ************************************************************************* // diff --git a/applications/solvers/immersedBoundary/icoIbFoam/Make/files b/applications/solvers/immersedBoundary/icoIbFoam/Make/files deleted file mode 100644 index bccc003ee..000000000 --- a/applications/solvers/immersedBoundary/icoIbFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -icoIbFoam.C - -EXE = $(FOAM_APPBIN)/icoIbFoam diff --git a/applications/solvers/immersedBoundary/icoIbFoam/Make/options b/applications/solvers/immersedBoundary/icoIbFoam/Make/options deleted file mode 100644 index b2ecef30a..000000000 --- a/applications/solvers/immersedBoundary/icoIbFoam/Make/options +++ /dev/null @@ -1,13 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude - -EXE_LIBS = \ - -lfiniteVolume \ - -lmeshTools \ - -lsurfMesh \ - -lsampling \ - -ldynamicMesh \ - -limmersedBoundary \ - -llduSolvers diff --git a/applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C b/applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C deleted file mode 100644 index 5447f7f6c..000000000 --- a/applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C +++ /dev/null @@ -1,139 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Application - icoIbFoam - -Description - Transient solver for incompressible, laminar flow of Newtonian fluids - with immersed boundary support. - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" -#include "pimpleControl.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ -# include "setRootCase.H" - -# include "createTime.H" -# include "createMesh.H" - - pimpleControl pimple(mesh); - -# include "createIbMasks.H" -# include "createFields.H" -# include "initContinuityErrs.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.loop()) - { - Info<< "Time = " << runTime.timeName() << nl << endl; - -# include "CourantNo.H" - - // Pressure-velocity corrector - while (pimple.loop()) - { - fvVectorMatrix UEqn - ( - fvm::ddt(U) - + fvm::div(phi, U) - - fvm::laplacian(nu, U) - ); - - if (pimple.momentumPredictor()) - { - solve(UEqn == -fvc::grad(p)); - } - - // --- PISO loop - while (pimple.correct()) - { - volScalarField rUA = 1.0/UEqn.A(); - - U = rUA*UEqn.H(); - // Immersed boundary update - U.correctBoundaryConditions(); - - phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf()); - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - adjustPhi(phi, U, p); - - // Non-orthogonal pressure corrector loop - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(rUA, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - pEqn.solve - ( - mesh.solutionDict().solver - ( - p.select(pimple.finalInnerIter()) - ) - ); - - if (pimple.finalNonOrthogonalIter()) - { - phi -= pEqn.flux(); - } - } - -# include "immersedBoundaryContinuityErrs.H" - - U -= rUA*fvc::grad(p); - U.correctBoundaryConditions(); - } - } - - runTime.write(); - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return 0; -} - - -// ************************************************************************* // diff --git a/applications/solvers/immersedBoundary/interIbFoam/Make/files b/applications/solvers/immersedBoundary/interIbFoam/Make/files deleted file mode 100644 index 5d3b09dcd..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -interIbFoam.C - -EXE = $(FOAM_APPBIN)/interIbFoam diff --git a/applications/solvers/immersedBoundary/interIbFoam/UEqn.H b/applications/solvers/immersedBoundary/interIbFoam/UEqn.H deleted file mode 100644 index bd641eeee..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/UEqn.H +++ /dev/null @@ -1,27 +0,0 @@ - surfaceScalarField muEff - ( - "muEff", - twoPhaseProperties.muf() - + fvc::interpolate(rho*turbulence->nut()) - ); - - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(rhoPhi, U) - - fvm::laplacian(muEff, U) - - (fvc::grad(U) & fvc::grad(muEff)) - ); - - UEqn.relax(); - - if (pimple.momentumPredictor()) - { - solve - ( - UEqn - == - - gh*fvc::grad(rho) - - fvc::grad(pd) - ); - } diff --git a/applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H b/applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H deleted file mode 100644 index 320ea795c..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H +++ /dev/null @@ -1,56 +0,0 @@ -{ - word alphaScheme("div(phi,alpha)"); - word alpharScheme("div(phirb,alpha)"); - - // New formulation of phir: Co condition - surfaceScalarField phir - ( - "phir", - faceIbMask*interface.cAlpha()*interface.nHatf()* - min - ( - scalar(0.5)/ // This is ref compression Co number for cAlpha = 1 - ( - mesh.surfaceInterpolation::deltaCoeffs()* - mesh.time().deltaT() - ), - mag(phi/mesh.magSf()) - + dimensionedScalar("small", dimVelocity, 1e-3) - ) - ); - - fvScalarMatrix alpha1Eqn - ( - fvm::ddt(alpha1) - + fvm::div(phi, alpha1, alphaScheme) - + fvm::div - ( - -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), - alpha1, - alpharScheme - ) - ); - - alpha1Eqn.relax(); - alpha1Eqn.solve(); - - Info<< "Liquid phase volume fraction = " - << alpha1.weightedAverage(mesh.V()).value() - << " Min(alpha1) = " << min(cellIbMask*alpha1).value() - << " Max(alpha1) = " << max(cellIbMask*alpha1).value() - << endl; - - rhoPhi = faceIbMask*(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2); - - // Limit alpha to handle IB cells -// alpha1.max(0); -// alpha1.min(1); - - // Update of interface and density - interface.correct(); - twoPhaseProperties.correct(); - - // Calculate density using limited alpha1 - rho = twoPhaseProperties.rho(); - rho.correctBoundaryConditions(); -} diff --git a/applications/solvers/immersedBoundary/interIbFoam/correctPhi.H b/applications/solvers/immersedBoundary/interIbFoam/correctPhi.H deleted file mode 100644 index 2f5bd5e2f..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/correctPhi.H +++ /dev/null @@ -1,58 +0,0 @@ -{ -# include "continuityErrs.H" - - wordList pcorrTypes - ( - pd.boundaryField().size(), - zeroGradientFvPatchScalarField::typeName - ); - - for (label i=0; i turbulence - ( - incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) - ); diff --git a/applications/solvers/immersedBoundary/interIbFoam/limitU.H b/applications/solvers/immersedBoundary/interIbFoam/limitU.H deleted file mode 100644 index c014ae74f..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/limitU.H +++ /dev/null @@ -1,24 +0,0 @@ -{ - scalar limitMagU = readScalar(pimple.dict().lookup("limitMagU")); - - volScalarField magU(mag(U)); - - scalar maxMagU = max(magU).value(); - - Info<< "mag(U): max: " << maxMagU - << " avg: " << magU.weightedAverage(mesh.V()).value(); - - if (maxMagU > limitMagU) - { - U.internalField() *= - neg(magU.internalField() - limitMagU) - + pos(magU.internalField() - limitMagU)* - limitMagU/(magU.internalField() + SMALL); - U.correctBoundaryConditions(); - Info << " ...limiting" << endl; - } - else - { - Info << endl; - } -} diff --git a/applications/solvers/immersedBoundary/interIbFoam/pEqn.H b/applications/solvers/immersedBoundary/interIbFoam/pEqn.H deleted file mode 100644 index 88238b82c..000000000 --- a/applications/solvers/immersedBoundary/interIbFoam/pEqn.H +++ /dev/null @@ -1,53 +0,0 @@ -{ - volScalarField rUA = 1.0/UEqn.A(); - surfaceScalarField rUAf = fvc::interpolate(rUA); - - U = rUA*UEqn.H(); - // Immersed boundary update - U.correctBoundaryConditions(); -# include "limitU.H" - - surfaceScalarField phiU - ( - "phiU", - faceIbMask*(fvc::interpolate(U) & mesh.Sf()) - ); - - phi = phiU - + faceIbMask* - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rUAf*mesh.magSf(); - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - adjustPhi(phi, U, pd); - - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rUAf, pd, "laplacian(rAU,pd)") == fvc::div(phi) - ); - - pdEqn.setReference(pdRefCell, pdRefValue); - - pdEqn.solve - ( - mesh.solutionDict().solver(pd.select(pimple.finalInnerIter())) - ); - - if (pimple.finalNonOrthogonalIter()) - { - phi -= pdEqn.flux(); - } - } - - // Explicitly relax pressure - pd.relax(); - - U += rUA*fvc::reconstruct((phi - phiU)/rUAf); - U.correctBoundaryConditions(); -} diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files new file mode 100644 index 000000000..46c7c1f2a --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files @@ -0,0 +1,3 @@ +pimpleDyMIbFoam.C + +EXE = $(FOAM_APPBIN)/pimpleDyMIbFoam diff --git a/applications/solvers/immersedBoundary/interIbFoam/Make/options b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options similarity index 55% rename from applications/solvers/immersedBoundary/interIbFoam/Make/options rename to applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options index 41ed9d472..a011e27ec 100644 --- a/applications/solvers/immersedBoundary/interIbFoam/Make/options +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options @@ -1,22 +1,23 @@ EXE_INC = \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ - -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel EXE_LIBS = \ - -linterfaceProperties \ + -lfiniteVolume \ + -limmersedBoundary \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools \ -lincompressibleTransportModels \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleLESModels \ - -lfiniteVolume \ - -lmeshTools \ - -lsurfMesh \ - -lsampling \ - -ldynamicMesh \ - -limmersedBoundary \ - -llduSolvers + -llduSolvers \ + -L$(MESQUITE_LIB_DIR) -lmesquite diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H new file mode 100644 index 000000000..44acc8af0 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H @@ -0,0 +1,26 @@ + // Convection-diffusion matrix + fvVectorMatrix HUEqn + ( + fvm::div(phi, U) + + turbulence->divDevReff() + ); + + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + // Get under-relaxation factor + scalar UUrf = + mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter())); + + if (pimple.momentumPredictor()) + { + // Solve momentum predictor + solve + ( + ddtUEqn + + relax(HUEqn, UUrf) + == + - fvc::grad(p), + mesh.solutionDict().solver((U.select(pimple.finalIter()))) + ); + } diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H new file mode 100644 index 000000000..eda2e9eda --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H @@ -0,0 +1,31 @@ +{ + volScalarField pcorr("pcorr", p); + + // Initialise flux with interpolated velocity + phi = fvc::interpolate(U) & mesh.Sf(); + + adjustPhi(phi, U, pcorr); + + mesh.schemesDict().setFluxRequired(pcorr.name()); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(1/aU, pcorr) == fvc::div(phi) + ); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + + // Fluxes are corrected to absolute velocity and further corrected + // later. HJ, 6/Feb/2009 + } +# include "continuityErrs.H" +} + diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H new file mode 100644 index 000000000..ca5e25906 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H @@ -0,0 +1,11 @@ +#include "createTimeControls.H" + +bool correctPhi +( + pimple.dict().lookupOrDefault("correctPhi", false) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); diff --git a/applications/solvers/immersedBoundary/icoIbFoam/createFields.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H similarity index 61% rename from applications/solvers/immersedBoundary/icoIbFoam/createFields.H rename to applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H index 5af8f976c..74a461355 100644 --- a/applications/solvers/immersedBoundary/icoIbFoam/createFields.H +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H @@ -1,23 +1,3 @@ - Info<< "Reading transportProperties\n" << endl; - - IOdictionary transportProperties - ( - IOobject - ( - "transportProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ_IF_MODIFIED, - IOobject::NO_WRITE - ) - ); - - dimensionedScalar nu - ( - transportProperties.lookup("nu") - ); - - Info<< "Reading field p\n" << endl; volScalarField p ( @@ -54,3 +34,26 @@ scalar pRefValue = 0.0; setRefCell(p, pimple.dict(), pRefCell, pRefValue); mesh.schemesDict().setFluxRequired(p.name()); + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + + Info<< "Reading field aU if present\n" << endl; + volScalarField aU + ( + IOobject + ( + "aU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + 1/runTime.deltaT(), + zeroGradientFvPatchScalarField::typeName + ); diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H new file mode 100644 index 000000000..1c02221ba --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H @@ -0,0 +1,58 @@ +{ + p.boundaryField().updateCoeffs(); + + // Prepare clean Ap without time derivative contribution and + // without contribution from under-relaxation + // HJ, 26/Oct/2015 + aU = HUEqn.A(); + + // Store velocity under-relaxation point before using U for the flux + // precursor + U.storePrevIter(); + + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); + + // Global flux balance + adjustPhi(phi, U, p); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve + ( + mesh.solutionDict().solver(p.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + } + + // Explicitly relax pressure for momentum corrector except for last corrector + if (!pimple.finalIter()) + { + p.relax(); + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "movingMeshContinuityErrs.H" + + U = UUrf* + ( + 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ) + ) + + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C similarity index 66% rename from applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C rename to applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C index 8f9b2738d..f50b32a1c 100644 --- a/applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C @@ -22,19 +22,17 @@ License along with foam-extend. If not, see . Application - interIbFoam + pimpleDyMFoam.C Description - Solver for 2 incompressible, isothermal immiscible fluids using a VOF - (volume of fluid) phase-fraction based interface capturing approach, - with immersed boundary support - - The momentum and other fluid properties are of the "mixture" and a single - momentum equation is solved. + Transient solver for incompressible, flow of Newtonian fluids + with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm. Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. - For a two-fluid approach see twoPhaseEulerFoam. + Consistent formulation without time-step and relaxation dependence by Jasak + + Support for immersed boundary Author Hrvoje Jasak, Wikki Ltd. All rights reserved @@ -42,32 +40,30 @@ Author \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "interfaceProperties.H" -#include "twoPhaseMixture.H" +#include "singlePhaseTransportModel.H" #include "turbulenceModel.H" +#include "dynamicFvMesh.H" #include "pimpleControl.H" +#include "immersedBoundaryPolyPatch.H" #include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" +#include "emptyFvPatch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { # include "setRootCase.H" + # include "createTime.H" -# include "createMesh.H" +# include "createDynamicFvMesh.H" pimpleControl pimple(mesh); -# include "readGravitationalAcceleration.H" # include "initContinuityErrs.H" -# include "createFields.H" # include "createIbMasks.H" -# include "createTimeControls.H" -# include "correctPhi.H" -# include "CourantNo.H" -# include "setInitialDeltaT.H" +# include "createFields.H" +# include "createControls.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -75,19 +71,46 @@ int main(int argc, char *argv[]) while (runTime.run()) { -# include "readTimeControls.H" -# include "immersedBoundaryCourantNo.H" +# include "readControls.H" +# include "CourantNo.H" # include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; - // Pressure-velocity corrector + bool meshChanged = mesh.update(); + + U.correctBoundaryConditions(); + p.correctBoundaryConditions(); + aU.correctBoundaryConditions(); + phi.correctBoundaryConditions(); + turbulence->correct(); + +# include "updateIbMasks.H" +# include "volContinuity.H" + + if (runTime.outputTime()) + { + volScalarField divMeshPhi("divMeshPhi", mag(fvc::surfaceIntegrate(mesh.phi()))); + divMeshPhi.write(); + } + + if (checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + + // Fluxes will be corrected to absolute velocity + // HJ, 6/Feb/2009 +# include "correctPhi.H" + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + // --- PIMPLE loop while (pimple.loop()) { - twoPhaseProperties.correct(); - # include "UEqn.H" // --- PISO loop @@ -96,15 +119,6 @@ int main(int argc, char *argv[]) # include "pEqn.H" } -# include "immersedBoundaryContinuityErrs.H" - -# include "limitU.H" - - p = pd + rho*gh; - p.correctBoundaryConditions(); - -# include "alphaEqn.H" - turbulence->correct(); } diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H new file mode 100644 index 000000000..9f982e260 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H @@ -0,0 +1,5 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", false); + +checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false); diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/files b/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/files deleted file mode 100644 index 1a063f7e8..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -porousSimpleIbFoam.C - -EXE = $(FOAM_APPBIN)/porousSimpleIbFoam diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/options b/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/options deleted file mode 100644 index 212eee945..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/options +++ /dev/null @@ -1,21 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/turbulenceModels \ - -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude - -EXE_LIBS = \ - -lincompressibleRASModels \ - -lincompressibleTransportModels \ - -lincompressibleTurbulenceModel \ - -lfiniteVolume \ - -lmeshTools \ - -lsurfMesh \ - -lsampling \ - -ldynamicMesh \ - -limmersedBoundary \ - -limmersedBoundaryTurbulence \ - -llduSolvers diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H b/applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H deleted file mode 100644 index c6ab47056..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H +++ /dev/null @@ -1,12 +0,0 @@ - // Solve the Momentum equation - tmp UEqn - ( - fvm::div(phi, U) - + turbulence->divDevReff() - ); - - UEqn().relax(); - - pZones.addResistance(UEqn()); - - solve(UEqn() == -fvc::grad(p)); diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H b/applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H deleted file mode 100644 index a6ba38265..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H +++ /dev/null @@ -1,44 +0,0 @@ - Info << "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info << "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - -# include "createPhi.H" - - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, simple.dict(), pRefCell, pRefValue); - mesh.schemesDict().setFluxRequired(p.name()); - - singlePhaseTransportModel laminarTransport(U, phi); - - autoPtr turbulence - ( - incompressible::RASModel::New(U, phi, laminarTransport) - ); - - porousZones pZones(mesh); diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H b/applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H deleted file mode 100644 index 595165c48..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H +++ /dev/null @@ -1,41 +0,0 @@ -{ - volScalarField AU = UEqn().A(); - - U = UEqn().H()/AU; - // Immersed boundary update - U.correctBoundaryConditions(); - - UEqn.clear(); - phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf()); - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - adjustPhi(phi, U, p); - - // Non-orthogonal pressure corrector loop - while (simple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(1.0/AU, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi -= pEqn.flux(); - } - } - -# include "immersedBoundaryContinuityErrs.H" - - // Explicitly relax pressure for momentum corrector - p.relax(); - - // Momentum corrector - U -= fvc::grad(p)/AU; - U.correctBoundaryConditions(); -} diff --git a/applications/solvers/immersedBoundary/porousSimpleIbFoam/porousSimpleIbFoam.C b/applications/solvers/immersedBoundary/porousSimpleIbFoam/porousSimpleIbFoam.C deleted file mode 100644 index 3bf3a4fc3..000000000 --- a/applications/solvers/immersedBoundary/porousSimpleIbFoam/porousSimpleIbFoam.C +++ /dev/null @@ -1,88 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Application - simpleIbFoam - -Description - Steady-state solver for incompressible, turbulent flow - with porous zones and immersed boundary support. - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "singlePhaseTransportModel.H" -#include "RASModel.H" -#include "porousZones.H" -#include "simpleControl.H" - -#include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" - - simpleControl simple(mesh); - -# include "createIbMasks.H" -# include "createFields.H" -# include "initContinuityErrs.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (simple.loop()) - { - Info<< "Time = " << runTime.timeName() << nl << endl; - - // Pressure-velocity SIMPLE corrector - { -# include "UEqn.H" -# include "pEqn.H" - } - - turbulence->correct(); - - runTime.write(); - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return 0; -} - - -// ************************************************************************* // diff --git a/applications/solvers/immersedBoundary/potentialIbFoam/Make/files b/applications/solvers/immersedBoundary/potentialIbFoam/Make/files deleted file mode 100644 index a234c5801..000000000 --- a/applications/solvers/immersedBoundary/potentialIbFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -potentialIbFoam.C - -EXE = $(FOAM_APPBIN)/potentialIbFoam diff --git a/applications/solvers/immersedBoundary/potentialIbFoam/Make/options b/applications/solvers/immersedBoundary/potentialIbFoam/Make/options deleted file mode 100644 index 448c4b608..000000000 --- a/applications/solvers/immersedBoundary/potentialIbFoam/Make/options +++ /dev/null @@ -1,16 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ - -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude - -EXE_LIBS = \ - -lfiniteVolume \ - -lmeshTools \ - -lsurfMesh \ - -lsampling \ - -ldynamicMesh \ - -ldynamicFvMesh \ - -limmersedBoundary \ - -llduSolvers diff --git a/applications/solvers/immersedBoundary/potentialIbFoam/createFields.H b/applications/solvers/immersedBoundary/potentialIbFoam/createFields.H deleted file mode 100644 index 09ed08451..000000000 --- a/applications/solvers/immersedBoundary/potentialIbFoam/createFields.H +++ /dev/null @@ -1,50 +0,0 @@ - Info<< "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ), - mesh - ); - - // Do not reset p - - Info<< "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - // Do not reset U - - surfaceScalarField phi - ( - IOobject - ( - "phi", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - fvc::interpolate(U) & mesh.Sf() - ); - - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, simple.dict(), pRefCell, pRefValue); - mesh.schemesDict().setFluxRequired(p.name()); diff --git a/applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C b/applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C deleted file mode 100644 index bfd7273a9..000000000 --- a/applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C +++ /dev/null @@ -1,231 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Application - potentialIbFoam - -Description - Potential flow solver with immersed boundary support. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "immersedBoundaryFvPatch.H" -#include "immersedBoundaryAdjustPhi.H" -#include "simpleControl.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - argList::validOptions.insert("writep", ""); - -# include "setRootCase.H" - -# include "createTime.H" -# include "createMesh.H" - - simpleControl simple(mesh); - -# include "createIbMasks.H" -# include "createFields.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< nl << "Calculating potential flow" << endl; - - // Do correctors over the complete set - while (simple.correctNonOrthogonal()) - { - phi = faceIbMask*(linearInterpolate(U) & mesh.Sf()); - - // Adjust immersed boundary fluxes - immersedBoundaryAdjustPhi(phi, U); - - // Adjust fluxes - adjustPhi(phi, U, p); - - p.storePrevIter(); - - fvScalarMatrix pEqn - ( - fvm::laplacian - ( - dimensionedScalar - ( - "1", - dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0), - 1 - ), - p - ) - == - fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - pEqn.solve(); - - // Correct the flux - phi -= pEqn.flux(); - - if (!simple.finalNonOrthogonalIter()) - { - p.relax(); - } - - Info<< "p min " << gMin(p.internalField()) - << " max " << gMax(p.internalField()) - << " masked min " - << gMin(cellIbMask.internalField()*p.internalField()) - << " max " - << gMax(cellIbMask.internalField()*p.internalField()) - << endl; - - Info<< "continuity error = " - << mag - ( - fvc::div(faceIbMask*phi) - )().weightedAverage(mesh.V()).value() - << endl; - - Info<< "Contour continuity error = " - << mag(sum(phi.boundaryField())) - << endl; - - U = fvc::reconstruct(phi); - U.correctBoundaryConditions(); - - Info<< "Interpolated U error = " - << ( - sqrt - ( - sum - ( - sqr - ( - faceIbMask* - ( - fvc::interpolate(U) & mesh.Sf() - ) - - phi - ) - ) - )/sum(mesh.magSf()) - ).value() - << endl; - } - - // Calculate velocity magnitude - { - volScalarField magU = cellIbMask*mag(U); - - Info << "IB-masked mag(U): max: " << gMax(magU.internalField()) - << " min: " << gMin(magU.internalField()) << endl; - } - - // Force the write - U.write(); - phi.write(); - - cellIbMask.write(); - cellIbMaskExt.write(); - - if (args.optionFound("writep")) - { - // Find reference patch - label refPatch = -1; - scalar maxMagU = 0; - - // Go through all velocity patches and find the one that fixes - // velocity to the largest value - - forAll (U.boundaryField(), patchI) - { - const fvPatchVectorField& Upatch = U.boundaryField()[patchI]; - - if (Upatch.fixesValue()) - { - // Calculate mean velocity - scalar u = sum(mag(Upatch)); - label patchSize = Upatch.size(); - - reduce(u, sumOp()); - reduce(patchSize, sumOp