diff --git a/applications/solvers/coupled/Allwmake b/applications/solvers/coupled/Allwmake index 9ce03ddae..e6320de43 100755 --- a/applications/solvers/coupled/Allwmake +++ b/applications/solvers/coupled/Allwmake @@ -8,6 +8,6 @@ wmake blockCoupledScalarTransportFoam wmake conjugateHeatFoam wmake conjugateHeatSimpleFoam wmake pUCoupledFoam -wmake pUCoupledFullPicardFoam -wmake pUCoupledSemiPicardFoam +wmake MRFPorousFoam + # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C new file mode 100644 index 000000000..f9fb25c1d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C @@ -0,0 +1,118 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Application + MRFPorousFoam + +Description + Steady-state solver for incompressible, turbulent flow, with implicit + coupling between pressure and velocity achieved by fvBlockMatrix. + Turbulence is in this version solved using the existing turbulence + structure. + + Added support for MRF and porous zones + +Authors + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "fvBlockMatrix.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" +#include "MRFZones.H" +#include "porousZones.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "createFields.H" +# include "createMRF.H" +# include "createPorous.H" +# include "initContinuityErrs.H" +# include "initConvergenceCheck.H" + + Info<< "\nStarting time loop\n" << endl; + while (runTime.loop()) + { +# include "readBlockSolverControls.H" +# include "readFieldBounds.H" + + Info<< "Time = " << runTime.timeName() << nl << endl; + + p.storePrevIter(); + + // Initialize the Up block system (matrix, source and reference to Up) + fvBlockMatrix UpEqn(Up); + + // Assemble and insert momentum equation +# include "UEqn.H" + + // Assemble and insert pressure equation +# include "pEqn.H" + + // Assemble and insert coupling terms +# include "couplingTerms.H" + + // Solve the block matrix + maxResidual = cmptMax(UpEqn.solve().initialResidual()); + + // Retrieve solution + UpEqn.retrieveSolution(0, U.internalField()); + UpEqn.retrieveSolution(3, p.internalField()); + + U.correctBoundaryConditions(); + p.correctBoundaryConditions(); + + phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; + + // Make flux relative in rotating zones + mrfZones.relativeFlux(phi); + +# include "continuityErrs.H" + +# include "boundPU.H" + + p.relax(); + + turbulence->correct(); + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + +# include "convergenceCheck.H" + } + + Info<< "End\n" << endl; + + return 0; +} diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/files b/applications/solvers/coupled/MRFPorousFoam/Make/files new file mode 100644 index 000000000..2ae8a1e4f --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/Make/files @@ -0,0 +1,3 @@ +MRFPorousFoam.C + +EXE = $(FOAM_APPBIN)/MRFPorousFoam diff --git a/applications/solvers/coupled/MRFPorousFoam/Make/options b/applications/solvers/coupled/MRFPorousFoam/Make/options new file mode 100644 index 000000000..5d2b99e6d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/Make/options @@ -0,0 +1,12 @@ +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 + +EXE_LIBS = \ + -lincompressibleRASModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -llduSolvers diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H new file mode 100644 index 000000000..ed0e2bc44 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/UEqn.H @@ -0,0 +1,15 @@ + // Momentum equation + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + + turbulence->divDevReff(U) + ); + + // Add MRF and porous sources + mrfZones.addCoriolis(UEqn); + pZones.addResistance(UEqn); + + UEqn.relax(); + + UpEqn.insertEquation(0, UEqn); diff --git a/applications/solvers/coupled/MRFPorousFoam/boundPU.H b/applications/solvers/coupled/MRFPorousFoam/boundPU.H new file mode 100644 index 000000000..eca655e8c --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/boundPU.H @@ -0,0 +1,32 @@ +{ + // Bound the pressure + dimensionedScalar p1 = min(p); + dimensionedScalar p2 = max(p); + + if (p1 < pMin || p2 > pMax) + { + Info<< "p: " << p1.value() << " " << p2.value() + << ". Bounding." << endl; + + p.max(pMin); + p.min(pMax); + p.correctBoundaryConditions(); + } + + // Bound the velocity + volScalarField magU = mag(U); + dimensionedScalar U1 = max(magU); + + if (U1 > UMax) + { + Info<< "U: " << U1.value() << ". Bounding." << endl; + + volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU) + + neg(magU - UMax); + Ulimiter.max(scalar(0)); + Ulimiter.min(scalar(1)); + + U *= Ulimiter; + U.correctBoundaryConditions(); + } +} diff --git a/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H new file mode 100644 index 000000000..9c7734131 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/convergenceCheck.H @@ -0,0 +1,8 @@ +// Check convergence +if (maxResidual < convergenceCriterion) +{ + Info<< "reached convergence criterion: " << convergenceCriterion << endl; + runTime.writeAndEnd(); + Info<< "latestTime = " << runTime.timeName() << endl; +} + diff --git a/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H b/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/couplingTerms.H @@ -0,0 +1,15 @@ +{ + // Calculate grad p coupling matrix. Needs to be here if one uses + // gradient schemes with limiters. VV, 9/June/2014 + BlockLduSystem pInU(fvm::grad(p)); + + // Calculate div U coupling. Could be calculated only once since + // it is only geometry dependent. VV, 9/June/2014 + BlockLduSystem UInp(fvm::UDiv(U)); + + // Last argument in insertBlockCoupling says if the column direction + // should be incremented. This is needed for arbitrary positioning + // of U and p in the system. This could be better. VV, 30/April/2014 + UpEqn.insertBlockCoupling(0, 3, pInU, true); + UpEqn.insertBlockCoupling(3, 0, UInp, false); +} diff --git a/applications/solvers/coupled/MRFPorousFoam/createFields.H b/applications/solvers/coupled/MRFPorousFoam/createFields.H new file mode 100644 index 000000000..51403c4a7 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createFields.H @@ -0,0 +1,52 @@ +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" + +singlePhaseTransportModel laminarTransport(U, phi); +autoPtr turbulence +( + incompressible::RASModel::New(U, phi, laminarTransport) +); + +// Block vector field for velocity (first entry) and pressure (second +// entry). +Info << "Creating field Up\n" << endl; +volVector4Field Up +( + IOobject + ( + "Up", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector4("zero", dimless, vector4::zero) +); diff --git a/applications/solvers/coupled/MRFPorousFoam/createMRF.H b/applications/solvers/coupled/MRFPorousFoam/createMRF.H new file mode 100644 index 000000000..161446a8e --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createMRF.H @@ -0,0 +1,2 @@ + MRFZones mrfZones(mesh); + mrfZones.correctBoundaryVelocity(U); diff --git a/applications/solvers/coupled/MRFPorousFoam/createPorous.H b/applications/solvers/coupled/MRFPorousFoam/createPorous.H new file mode 100644 index 000000000..430b466aa --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/createPorous.H @@ -0,0 +1 @@ + porousZones pZones(mesh); diff --git a/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H new file mode 100644 index 000000000..482bc421d --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H @@ -0,0 +1,4 @@ +// initialize values for convergence checks + + scalar maxResidual = 0; + scalar convergenceCriterion = 0; diff --git a/applications/solvers/coupled/MRFPorousFoam/pEqn.H b/applications/solvers/coupled/MRFPorousFoam/pEqn.H new file mode 100644 index 000000000..9029b873f --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/pEqn.H @@ -0,0 +1,23 @@ +// Pressure parts of the continuity equation +surfaceScalarField rUAf +( + "rUAf", + fvc::interpolate(1.0/UEqn.A()) +); + +surfaceScalarField presSource +( + "presSource", + rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf()) +); + +fvScalarMatrix pEqn +( + - fvm::laplacian(rUAf, p) + == + - fvc::div(presSource) +); + +pEqn.setReference(pRefCell, pRefValue); + +UpEqn.insertEquation(3, pEqn); diff --git a/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H new file mode 100644 index 000000000..93c2bda47 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/readBlockSolverControls.H @@ -0,0 +1,15 @@ + label pRefCell = 0; + scalar pRefValue = 0; + setRefCell + ( + p, + mesh.solutionDict().subDict("blockSolver"), + pRefCell, + pRefValue + ); + + mesh.solutionDict().subDict("blockSolver").readIfPresent + ( + "convergence", + convergenceCriterion + ); diff --git a/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H b/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H new file mode 100644 index 000000000..4a6c6f787 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/readFieldBounds.H @@ -0,0 +1,14 @@ + // Read field bounds + dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds"); + + // Pressure bounds + dimensionedScalar pMin("pMin", p.dimensions(), 0); + dimensionedScalar pMax("pMax", p.dimensions(), 0); + + fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value(); + + // Velocity bound + dimensionedScalar UMax("UMax", U.dimensions(), 0); + + fieldBounds.lookup(U.name()) >> UMax.value(); + dimensionedScalar smallU("smallU", dimVelocity, 1e-10); diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C index cc2b1f7e0..25b389b87 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C @@ -38,13 +38,8 @@ Description #include "fieldTypes.H" #include "Time.H" #include "fvMesh.H" +#include "fvBlockMatrix.H" -#include "blockLduSolvers.H" -#include "VectorNFieldTypes.H" -#include "volVectorNFields.H" -#include "blockVectorNMatrices.H" - -#include "blockMatrixTools.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -92,53 +87,27 @@ int main(int argc, char *argv[]) TsEqn.relax(); // Prepare block system - BlockLduMatrix blockM(mesh); - - //- Transfer the coupled interface list for processor/cyclic/etc. - // boundaries - blockM.interfaces() = blockT.boundaryField().blockInterfaces(); - - // Grab block diagonal and set it to zero - Field& d = blockM.diag().asSquare(); - d = tensor2::zero; - - // Grab linear off-diagonal and set it to zero - Field& l = blockM.lower().asLinear(); - Field& u = blockM.upper().asLinear(); - u = vector2::zero; - l = vector2::zero; - - vector2Field& blockX = blockT.internalField(); - vector2Field blockB(mesh.nCells(), vector2::zero); + fvBlockMatrix blockM(blockT); //- Inset equations into block Matrix - blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB); - blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB); + blockM.insertEquation(0, TEqn); + blockM.insertEquation(1, TsEqn); - //- Add off-diagonal terms and remove from block source - forAll(d, i) - { - d[i](0, 1) = -alpha.value()*mesh.V()[i]; - d[i](1, 0) = -alpha.value()*mesh.V()[i]; + //- Add off-diagonal coupling terms + scalarField coupling(mesh.nCells(), -alpha.value()); - blockB[i][0] -= alpha.value()*blockX[i][1]*mesh.V()[i]; - blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i]; - } + blockM.insertEquationCoupling(0, 1, coupling); + blockM.insertEquationCoupling(1, 0, coupling); + + // Update source coupling: coupling terms eliminated from source + blockM.updateSourceCoupling(); //- Block coupled solver call - BlockSolverPerformance solverPerf = - BlockLduSolver::New - ( - blockT.name(), - blockM, - mesh.solutionDict().solver(blockT.name()) - )->solve(blockX, blockB); - - solverPerf.print(); + blockM.solve(); // Retrieve solution - blockMatrixTools::blockRetrieve(0, T.internalField(), blockX); - blockMatrixTools::blockRetrieve(1, Ts.internalField(), blockX); + blockM.retrieveSolution(0, T.internalField()); + blockM.retrieveSolution(1, Ts.internalField()); T.correctBoundaryConditions(); Ts.correctBoundaryConditions(); diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H index 2f4eb6aaa..d3d4b42a8 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H @@ -44,15 +44,6 @@ dimensionedVector2("zero", dimless, vector2::zero) ); - { - vector2Field& blockX = blockT.internalField(); - - blockMatrixTools::blockInsert(0, T.internalField(), blockX); - blockMatrixTools::blockInsert(1, Ts.internalField(), blockX); - } - - blockT.write(); - Info<< "Reading field U\n" << endl; volVectorField U diff --git a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C index f8b206c57..b66a430df 100644 --- a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C +++ b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C @@ -46,7 +46,7 @@ ersPointSource::ersPointSource alpha_(readScalar(dict.lookup("alpha"))), direction_(dict.lookup("direction")) { - q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), 0.0); + q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), scalar(0)); } diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C index 4804383a1..2af8a4675 100644 --- a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C +++ b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C @@ -30,18 +30,10 @@ Description #include "contactPatchPair.H" #include "volFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from components -contactPatchPair::contactPatchPair + Foam::contactPatchPair::contactPatchPair ( const fvMesh& m, const label master, @@ -75,6 +67,4 @@ contactPatchPair::contactPatchPair {} -} // End namespace Foam - // ************************************************************************* // diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C index 9e5cf7992..59a70f0e4 100644 --- a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C +++ b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C @@ -302,7 +302,7 @@ void Foam::contactPatchPair::correct ( slavePressure ), - 0.0 + scalar(0) ); // Calculate master traction, using the positive part of diff --git a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H index 50bb49ddd..a9f4c0b83 100644 --- a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H +++ b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H @@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0; // scalarField effTraction = // cohesiveZone.internalField() * // mag(traction.internalField()); - scalarField normalTraction = - cohesiveZone.internalField() * - ( n.internalField() & traction.internalField() ); - normalTraction = max(normalTraction, 0.0); // only consider tensile tractions - scalarField shearTraction = - cohesiveZone.internalField() * - mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); + scalarField normalTraction = + cohesiveZone.internalField()* + ( n.internalField() & traction.internalField() ); + + // only consider tensile tractions + normalTraction = max(normalTraction, scalar(0)); + scalarField shearTraction = + cohesiveZone.internalField() * + mag( (I - Foam::sqr(n.internalField())) & traction.internalField() ); // the traction fraction is monitored to decide which faces to break: // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face @@ -27,16 +29,18 @@ nCoupledFacesToBreak = 0; scalarField effTractionFraction(normalTraction.size(), 0.0); if(cohesivePatchUPtr) - { - effTractionFraction = - (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); - } + { + effTractionFraction = + (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + + (shearTraction/tauMaxI)*(shearTraction/tauMaxI); + } else - { - // solidCohesiveFixedModeMix only uses sigmaMax - effTractionFraction = - (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); - } + { + // solidCohesiveFixedModeMix only uses sigmaMax + effTractionFraction = + (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI); + } maxEffTractionFraction = gMax(effTractionFraction); SLList