From e8615a3ce9adfb476703b80cc6cbcc7d81138fda Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 22 May 2018 10:31:51 +0100 Subject: [PATCH] Temporary implicit block boundary conditions and porous/MRF --- .../solvers/coupled/MRFPorousFoam/UEqn.H | 78 +------------------ .../coupled/MRFPorousFoam/addBlockCoupledBC.H | 56 +++++++++++++ .../MRFPorousFoam/addImplicitMRFPorous.H | 77 ++++++++++++++++++ .../solvers/coupled/pUCoupledFoam/UEqn.H | 2 + .../coupled/pUCoupledFoam/addBlockCoupledBC.H | 56 +++++++++++++ 5 files changed, 193 insertions(+), 76 deletions(-) create mode 100644 applications/solvers/coupled/MRFPorousFoam/addBlockCoupledBC.H create mode 100644 applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/addBlockCoupledBC.H diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H index a189057e0..55c804f37 100644 --- a/applications/solvers/coupled/MRFPorousFoam/UEqn.H +++ b/applications/solvers/coupled/MRFPorousFoam/UEqn.H @@ -53,80 +53,6 @@ // Insert momentum equation UpEqn.insertEquation(0, UEqn); - if (addMRF || addPorosity) - { - // Manually over-ride the 3x3 block to handle the off-diagonal - // part of the Ap coefficient - const tensorField& TUIn = tTU().internalField(); - - CoeffField::squareTypeField& DD = UpEqn.diag().asSquare(); - - const scalarField& V = mesh.V().field(); - - // 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(); - - // 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(); - } - } - } +# include "addImplicitMRFPorous.H" +# include "addBlockCoupledBC.H" } diff --git a/applications/solvers/coupled/MRFPorousFoam/addBlockCoupledBC.H b/applications/solvers/coupled/MRFPorousFoam/addBlockCoupledBC.H new file mode 100644 index 000000000..b80aeb90e --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/addBlockCoupledBC.H @@ -0,0 +1,56 @@ + // Hack block-coupled boundary conditions: due for rewrite + const volScalarField nuEff = turbulence->nuEff(); + + forAll (U.boundaryField(), patchI) + { + if (U.boundaryField()[patchI].blockCoupled()) + { + // Insert correcting fully implicit coupling coefficient + + const labelList fc = mesh.boundary()[patchI].faceCells(); + + const fvPatchVectorField& Up = U.boundaryField()[patchI]; + + // Warning: hacked for nuEff in viscosity + const scalarField nutpMagSf = + nuEff.boundaryField()[patchI]* + mesh.magSf().boundaryField()[patchI]; + + // Get boundary condition contribution to matrix diagonal + tensorField patchDiag = + -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf; + + // Get matrix diagonal + CoeffField::squareTypeField& blockDiag = + UpEqn.diag().asSquare(); + + forAll (fc, faceI) + { + blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx(); + blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy(); + blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz(); + + blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx(); + blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy(); + blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz(); + + blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx(); + blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy(); + blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz(); + } + + // Get boundary condition contribution to matrix source + vectorField patchSource = + -Up.blockGradientBoundaryCoeffs()*nutpMagSf; + + // Get matrix source + Field& blockSource = UpEqn.source(); + + forAll (fc, faceI) + { + blockSource[fc[faceI]](0) -= patchSource[faceI](0); + blockSource[fc[faceI]](1) -= patchSource[faceI](1); + blockSource[fc[faceI]](2) -= patchSource[faceI](2); + } + } + } diff --git a/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H b/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H new file mode 100644 index 000000000..5ba812f51 --- /dev/null +++ b/applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H @@ -0,0 +1,77 @@ + + if (addMRF || addPorosity) + { + // Manually over-ride the 3x3 block to handle the off-diagonal + // part of the Ap coefficient + const tensorField& TUIn = tTU().internalField(); + + CoeffField::squareTypeField& DD = UpEqn.diag().asSquare(); + + const scalarField& V = mesh.V().field(); + + // 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(); + + // 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(); + } + } + } diff --git a/applications/solvers/coupled/pUCoupledFoam/UEqn.H b/applications/solvers/coupled/pUCoupledFoam/UEqn.H index e1535831e..3e9907023 100644 --- a/applications/solvers/coupled/pUCoupledFoam/UEqn.H +++ b/applications/solvers/coupled/pUCoupledFoam/UEqn.H @@ -12,4 +12,6 @@ UEqn.relax(); UpEqn.insertEquation(0, UEqn); + +# include "addBlockCoupledBC.H" } diff --git a/applications/solvers/coupled/pUCoupledFoam/addBlockCoupledBC.H b/applications/solvers/coupled/pUCoupledFoam/addBlockCoupledBC.H new file mode 100644 index 000000000..b80aeb90e --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/addBlockCoupledBC.H @@ -0,0 +1,56 @@ + // Hack block-coupled boundary conditions: due for rewrite + const volScalarField nuEff = turbulence->nuEff(); + + forAll (U.boundaryField(), patchI) + { + if (U.boundaryField()[patchI].blockCoupled()) + { + // Insert correcting fully implicit coupling coefficient + + const labelList fc = mesh.boundary()[patchI].faceCells(); + + const fvPatchVectorField& Up = U.boundaryField()[patchI]; + + // Warning: hacked for nuEff in viscosity + const scalarField nutpMagSf = + nuEff.boundaryField()[patchI]* + mesh.magSf().boundaryField()[patchI]; + + // Get boundary condition contribution to matrix diagonal + tensorField patchDiag = + -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf; + + // Get matrix diagonal + CoeffField::squareTypeField& blockDiag = + UpEqn.diag().asSquare(); + + forAll (fc, faceI) + { + blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx(); + blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy(); + blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz(); + + blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx(); + blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy(); + blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz(); + + blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx(); + blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy(); + blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz(); + } + + // Get boundary condition contribution to matrix source + vectorField patchSource = + -Up.blockGradientBoundaryCoeffs()*nutpMagSf; + + // Get matrix source + Field& blockSource = UpEqn.source(); + + forAll (fc, faceI) + { + blockSource[fc[faceI]](0) -= patchSource[faceI](0); + blockSource[fc[faceI]](1) -= patchSource[faceI](1); + blockSource[fc[faceI]](2) -= patchSource[faceI](2); + } + } + }