diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C index 7b0a694e90..41845e6488 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C +++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -213,25 +213,17 @@ void reactingOneDim::updateFields() } -void reactingOneDim::updateMesh(const scalarField& mass0) +void reactingOneDim::updateMesh(const scalarField& deltaV) { - if (!moveMesh_) - { - return; - } - - const scalarField newV(mass0/rho_); - - Info<< "Initial/final volumes = " << gSum(regionMesh().V()) << ", " - << gSum(newV) << " [m3]" << endl; + Info<< "Initial/final volumes = " << gSum(deltaV) << endl; // move the mesh - const labelList moveMap = moveMesh(regionMesh().V() - newV, minimumDelta_); + const labelList moveMap = moveMesh(deltaV, minimumDelta_); // flag any cells that have not moved as non-reacting forAll(moveMap, i) { - if (moveMap[i] == 0) + if (moveMap[i] == 1) { solidChemistry_->setCellReacting(i, false); } @@ -246,28 +238,26 @@ void reactingOneDim::solveContinuity() Info<< "reactingOneDim::solveContinuity()" << endl; } - const scalarField mass0 = rho_*regionMesh().V(); - fvScalarMatrix rhoEqn - ( - fvm::ddt(rho_) - == - - solidChemistry_->RRg() - ); - - if (regionMesh().moving()) + if (!moveMesh_) { - surfaceScalarField phiRhoMesh + fvScalarMatrix rhoEqn ( - fvc::interpolate(rho_)*regionMesh().phi() + fvm::ddt(rho_) + == + - solidChemistry_->RRg() ); - rhoEqn += fvc::div(phiRhoMesh); + rhoEqn.solve(); } - rhoEqn.solve(); + if (moveMesh_) + { + const scalarField deltaV = + -solidChemistry_->RRg()*regionMesh().V()/rho_; - updateMesh(mass0); + updateMesh(deltaV); + } } @@ -298,7 +288,7 @@ void reactingOneDim::solveSpeciesMass() fvc::interpolate(Yi*rho_)*regionMesh().phi() ); - YiEqn += fvc::div(phiYiRhoMesh); + YiEqn -= fvc::div(phiYiRhoMesh); } @@ -329,7 +319,6 @@ void reactingOneDim::solveEnergy() - fvc::laplacian(kappa(), T()) == chemistrySh_ - - fvm::Sp(solidChemistry_->RRg(), h_) ); if (gasHSource_) @@ -351,7 +340,7 @@ void reactingOneDim::solveEnergy() fvc::interpolate(rho_*h_)*regionMesh().phi() ); - hEqn += fvc::div(phihMesh); + hEqn -= fvc::div(phihMesh); } hEqn.relax(); @@ -682,12 +671,6 @@ const surfaceScalarField& reactingOneDim::phiGas() const void reactingOneDim::preEvolveRegion() { pyrolysisModel::preEvolveRegion(); - - // Initialise all cells as able to react - forAll(h_, cellI) - { - solidChemistry_->setCellReacting(cellI, true); - } } diff --git a/src/regionModels/regionModel/regionModel/regionModel.C b/src/regionModels/regionModel/regionModel/regionModel.C index b1d516ce03..7bedbc5302 100644 --- a/src/regionModels/regionModel/regionModel/regionModel.C +++ b/src/regionModels/regionModel/regionModel/regionModel.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -201,7 +201,6 @@ bool Foam::regionModels::regionModel::read(const dictionary& dict) } infoOutput_.readIfPresent("infoOutput", dict); - return true; } else @@ -511,8 +510,6 @@ void Foam::regionModels::regionModel::evolve() Info<< "\nEvolving " << modelName_ << " for region " << regionMesh().name() << endl; - //read(); - preEvolveRegion(); evolveRegion(); diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.C b/src/regionModels/regionModel/regionModel1D/regionModel1D.C index d5e80beec8..072add55b1 100644 --- a/src/regionModels/regionModel/regionModel1D/regionModel1D.C +++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -75,6 +75,19 @@ void Foam::regionModels::regionModel1D::initialise() const polyBoundaryMesh& rbm = regionMesh().boundaryMesh(); + forAll(intCoupledPatchIDs_, i) + { + const label patchI = intCoupledPatchIDs_[i]; + const polyPatch& ppCoupled = rbm[patchI]; + localPyrolysisFaceI += ppCoupled.size(); + } + + boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI); + boundaryFaceFaces_.setSize(localPyrolysisFaceI); + boundaryFaceCells_.setSize(localPyrolysisFaceI); + + localPyrolysisFaceI = 0; + forAll(intCoupledPatchIDs_, i) { const label patchI = intCoupledPatchIDs_[i]; @@ -83,7 +96,6 @@ void Foam::regionModels::regionModel1D::initialise() { label faceI = ppCoupled.start() + localFaceI; label cellI = -1; - label nFaces = 0; label nCells = 0; do { @@ -99,14 +111,14 @@ void Foam::regionModels::regionModel1D::initialise() nCells++; cellIDs.append(cellI); const cell& cFaces = regionMesh().cells()[cellI]; - faceI = cFaces.opposingFaceLabel(faceI, regionMesh().faces()); faceIDs.append(faceI); - nFaces++; + label face0 = + cFaces.opposingFaceLabel(faceI, regionMesh().faces()); + faceI = face0; } while (regionMesh().isInternalFace(faceI)); boundaryFaceOppositeFace_[localPyrolysisFaceI] = faceI; - faceIDs.remove(); //remove boundary face. - nFaces--; + //faceIDs.remove(); //remove boundary face. boundaryFaceFaces_[localPyrolysisFaceI].transfer(faceIDs); boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs); @@ -115,10 +127,8 @@ void Foam::regionModels::regionModel1D::initialise() nLayers_ = nCells; } } - - boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI); - boundaryFaceFaces_.setSize(localPyrolysisFaceI); - boundaryFaceCells_.setSize(localPyrolysisFaceI); + faceIDs.clear(); + cellIDs.clear(); surfaceScalarField& nMagSf = nMagSfPtr_(); @@ -128,16 +138,22 @@ void Foam::regionModels::regionModel1D::initialise() const label patchI = intCoupledPatchIDs_[i]; const polyPatch& ppCoupled = rbm[patchI]; const vectorField& pNormals = ppCoupled.faceNormals(); + nMagSf.boundaryField()[patchI] = regionMesh().Sf().boundaryField()[patchI] & pNormals; + forAll(pNormals, localFaceI) { - const vector& n = pNormals[localFaceI]; + const vector n = pNormals[localFaceI]; const labelList& faces = boundaryFaceFaces_[localPyrolysisFaceI++]; - forAll (faces, faceI) + + forAll (faces, faceI) //faceI = 0 is on boundary { - const label faceID = faces[faceI]; - nMagSf[faceID] = regionMesh().Sf()[faceID] & n; + if (faceI > 0) + { + const label faceID = faces[faceI]; + nMagSf[faceID] = regionMesh().Sf()[faceID] & n; + } } } } @@ -150,8 +166,6 @@ bool Foam::regionModels::regionModel1D::read() { if (regionModel::read()) { - moveMesh_.readIfPresent("moveMesh", coeffs_); - return true; } else @@ -199,72 +213,65 @@ Foam::tmp Foam::regionModels::regionModel1D::moveMesh forAll(intCoupledPatchIDs_, localPatchI) { label patchI = intCoupledPatchIDs_[localPatchI]; - const polyPatch pp = bm[patchI]; - const vectorField& cf = regionMesh().Cf().boundaryField()[patchI]; + const polyPatch& pp = bm[patchI]; - forAll(pp, patchFaceI) + forAll (pp, patchFaceI) { const labelList& faces = boundaryFaceFaces_[totalFaceId]; const labelList& cells = boundaryFaceCells_[totalFaceId]; + const label oFace = boundaryFaceOppositeFace_[totalFaceId]; const vector n = pp.faceNormals()[patchFaceI]; const vector sf = pp.faceAreas()[patchFaceI]; - List oldCf(faces.size() + 1); - oldCf[0] = cf[patchFaceI]; - forAll(faces, i) + List oldCf(faces.size() + 1, vector::zero); + List frozen(faces.size(), false); + + forAll (faces, i) { - oldCf[i + 1] = regionMesh().faceCentres()[faces[i]]; + oldCf[i] = regionMesh().faceCentres()[faces[i]]; } - vector newDelta = vector::zero; - point nbrCf = oldCf[0]; + oldCf[faces.size()] = regionMesh().faceCentres()[oFace]; - forAll(faces, i) + forAll (faces, i) { - const label faceI = faces[i]; const label cellI = cells[i]; + if (mag(oldCf[i + 1] - oldCf[i]) < minDelta) + { + frozen[i] = true; + cellMoveMap[cellI] = 1; + } + } + + vectorField newDelta(cells.size() + 1, vector::zero); + + label j = 0; + forAllReverse (cells, i) + { + const label cellI = cells[i]; + newDelta[j+1] = (deltaV[cellI]/mag(sf))*n + newDelta[j]; + j++; + } + + forAll (faces, i) + { + const label faceI = faces[i]; const face f = regionMesh().faces()[faceI]; - newDelta += (deltaV[cellI]/mag(sf))*n; - - vector localDelta = vector::zero; forAll(f, pti) { const label pointI = f[pti]; - if - ( - mag((nbrCf - (oldPoints[pointI] + newDelta)) & n) - > minDelta - ) + if (!frozen[i]) { - newPoints[pointI] = oldPoints[pointI] + newDelta; - localDelta = newDelta; - cellMoveMap[cellI] = 1; + newPoints[pointI] = + oldPoints[pointI] + newDelta[newDelta.size() - 1 - i]; } } - nbrCf = oldCf[i + 1] + localDelta; - } - // Modify boundary - const label bFaceI = boundaryFaceOppositeFace_[totalFaceId]; - const face f = regionMesh().faces()[bFaceI]; - const label cellI = cells[cells.size() - 1]; - newDelta += (deltaV[cellI]/mag(sf))*n; - forAll(f, pti) - { - const label pointI = f[pti]; - if - ( - mag((nbrCf - (oldPoints[pointI] + newDelta)) & n) - > minDelta - ) - { - newPoints[pointI] = oldPoints[pointI] + newDelta; - cellMoveMap[cellI] = 1; - } } + totalFaceId ++; } } @@ -307,16 +314,15 @@ Foam::regionModels::regionModel1D::regionModel1D boundaryFaceOppositeFace_(regionMesh().nCells()), nLayers_(0), nMagSfPtr_(NULL), - moveMesh_(true) + moveMesh_(false) { if (active_) { constructMeshObjects(); initialise(); - if (readFields) { - read(); + moveMesh_.readIfPresent("moveMesh", coeffs_); } } } @@ -343,10 +349,9 @@ Foam::regionModels::regionModel1D::regionModel1D { constructMeshObjects(); initialise(); - if (readFields) { - read(dict); + moveMesh_.readIfPresent("moveMesh", coeffs_); } } }