diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H index 4bfc58978e0d1db1f81a660fb2b7ad64e87de8d9..a85e136eaa69ddf2d16747c8f80f007263a848ca 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H @@ -39,6 +39,15 @@ if (mesh.nInternalFaces()) fvc::surfaceSum(mag(phi))().internalField() ); + forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter) + { + sumPhi = max + ( + sumPhi, + fvc::surfaceSum(mag(iter().phi()))().internalField() + ); + } + CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); meanCoNum = diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H deleted file mode 100644 index 9bd18a9fcc9d429690ea142360915ace22a6246d..0000000000000000000000000000000000000000 --- a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H +++ /dev/null @@ -1,12 +0,0 @@ -# include "CourantNo.H" - -// { -// scalar UrCoNum = 0.5*gMax -// ( -// fvc::surfaceSum(mag(phi1 - phi2))().internalField()/mesh.V().field() -// )*runTime.deltaTValue(); - -// Info<< "Max Ur Courant Number = " << UrCoNum << endl; - -// CoNum = max(CoNum, UrCoNum); -// } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H index 6b859f068f975db90db81ab1b163bd013df474cb..484ceae4ac904d5bcd7dc2555f3448e1630b1d75 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H @@ -42,7 +42,7 @@ dimensionedScalar("phi", dimArea*dimVelocity, 0) ); - multiphaseSystem fluid(mesh, phi); + multiphaseSystem fluid(U, phi); forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter) { @@ -77,9 +77,8 @@ scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue); - singlePhaseTransportModel laminarTransport(U, phi); autoPtr<incompressible::LESModel> sgsModel ( - incompressible::LESModel::New(U, phi, laminarTransport) + incompressible::LESModel::New(U, phi, fluid) ); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files index 43869fde400a882d8730fc01875fab9a05fcd7a9..ff81afc6346c3148a276bdf56a1f09101db51f3b 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files @@ -8,6 +8,7 @@ dragModels/Gibilaro/Gibilaro.C dragModels/WenYu/WenYu.C dragModels/SyamlalOBrien/SyamlalOBrien.C dragModels/blended/blended.C +dragModels/interface/interface.C heatTransferModels/heatTransferModel/heatTransferModel.C heatTransferModels/heatTransferModel/newHeatTransferModel.C diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C new file mode 100644 index 0000000000000000000000000000000000000000..0ef07da220005b3993d5e76fff5f8703c8d8db3d --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "interface.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace dragModels +{ + defineTypeNameAndDebug(interface, 0); + + addToRunTimeSelectionTable + ( + dragModel, + interface, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::dragModels::interface::interface +( + const dictionary& interfaceDict, + const phaseModel& phase1, + const phaseModel& phase2 +) +: + dragModel(interfaceDict, phase1, phase2) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::dragModels::interface::~interface() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp<Foam::volScalarField> Foam::dragModels::interface::K +( + const volScalarField& Ur +) const +{ + return tmp<volScalarField> + ( + new volScalarField + ( + IOobject + ( + "K", + Ur.mesh().time().timeName(), + Ur.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + Ur.mesh(), + dimensionedScalar("K", dimDensity/dimTime, 0) + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H new file mode 100644 index 0000000000000000000000000000000000000000..bc7c8223cf6f9e9e28727acc6d2ffbf353b8198c --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H @@ -0,0 +1,92 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::dragModels::interface + +Description + Drag between phase separated by a VoF resolved interface. + +SourceFiles + interface.C + +\*---------------------------------------------------------------------------*/ + +#ifndef interface_H +#define interface_H + +#include "dragModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace dragModels +{ + +/*---------------------------------------------------------------------------*\ + Class interface Declaration +\*---------------------------------------------------------------------------*/ + +class interface +: + public dragModel +{ + +public: + + //- Runtime type information + TypeName("interface"); + + + // Constructors + + //- Construct from components + interface + ( + const dictionary& interfaceDict, + const phaseModel& phase1, + const phaseModel& phase2 + ); + + + //- Destructor + virtual ~interface(); + + + // Member Functions + + tmp<volScalarField> K(const volScalarField& Ur) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace dragModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C index 91a71ce60b6906a21f1c964b25d31cc9f1a81f47..4ca01f1a3a89f3bc9e75697057a8fd42e0fd8920 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) #include "initContinuityErrs.H" #include "readTimeControls.H" #include "correctPhi.H" - #include "CourantNos.H" + #include "CourantNo.H" #include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -68,7 +68,7 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readTimeControls.H" - #include "CourantNos.H" + #include "CourantNo.H" #include "setDeltaT.H" runTime++; @@ -80,6 +80,7 @@ int main(int argc, char *argv[]) sgsModel->correct(); fluid.solve(); rho = fluid.rho(); + #include "zonePhaseVolumes.H" //#include "interfacialCoeffs.H" //#include "TEqns.H" diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C new file mode 100644 index 0000000000000000000000000000000000000000..635b3d21f84e05f6b3c5a2138f84e4b36847ce3e --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C @@ -0,0 +1,181 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "multiphaseFixedFluxPressureFvPatchScalarField.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::multiphaseFixedFluxPressureFvPatchScalarField:: +multiphaseFixedFluxPressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF +) +: + fixedGradientFvPatchScalarField(p, iF), + phi0Name_("phi0"), + phiName_("phi"), + rhoName_("rho") +{} + + +Foam::multiphaseFixedFluxPressureFvPatchScalarField:: +multiphaseFixedFluxPressureFvPatchScalarField +( + const multiphaseFixedFluxPressureFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const fvPatchFieldMapper& mapper +) +: + fixedGradientFvPatchScalarField(ptf, p, iF, mapper), + phi0Name_(ptf.phi0Name_), + phiName_(ptf.phiName_), + rhoName_(ptf.rhoName_) +{} + + +Foam::multiphaseFixedFluxPressureFvPatchScalarField:: +multiphaseFixedFluxPressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const dictionary& dict +) +: + fixedGradientFvPatchScalarField(p, iF), + phi0Name_(dict.lookupOrDefault<word>("phi0", "phi0")), + phiName_(dict.lookupOrDefault<word>("phi", "phi")), + rhoName_(dict.lookupOrDefault<word>("rho", "rho")) +{ + if (dict.found("gradient")) + { + gradient() = scalarField("gradient", dict, p.size()); + fixedGradientFvPatchScalarField::updateCoeffs(); + fixedGradientFvPatchScalarField::evaluate(); + } + else + { + fvPatchField<scalar>::operator=(patchInternalField()); + gradient() = 0.0; + } +} + + +Foam::multiphaseFixedFluxPressureFvPatchScalarField:: +multiphaseFixedFluxPressureFvPatchScalarField +( + const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf +) +: + fixedGradientFvPatchScalarField(wbppsf), + phi0Name_(wbppsf.phi0Name_), + phiName_(wbppsf.phiName_), + rhoName_(wbppsf.rhoName_) +{} + + +Foam::multiphaseFixedFluxPressureFvPatchScalarField:: +multiphaseFixedFluxPressureFvPatchScalarField +( + const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf, + const DimensionedField<scalar, volMesh>& iF +) +: + fixedGradientFvPatchScalarField(wbppsf, iF), + phi0Name_(wbppsf.phi0Name_), + phiName_(wbppsf.phiName_), + rhoName_(wbppsf.rhoName_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const surfaceScalarField& phi0 = + db().lookupObject<surfaceScalarField>(phi0Name_); + + const surfaceScalarField& phi = + db().lookupObject<surfaceScalarField>(phiName_); + + fvsPatchField<scalar> phi0p = + patch().patchField<surfaceScalarField, scalar>(phi0); + + fvsPatchField<scalar> phip = + patch().patchField<surfaceScalarField, scalar>(phi); + + if (phi.dimensions() == dimDensity*dimVelocity*dimArea) + { + const fvPatchField<scalar>& rhop = + patch().lookupPatchField<volScalarField, scalar>(rhoName_); + + phip /= rhop; + } + + const fvsPatchField<scalar>& Dpp = + patch().lookupPatchField<surfaceScalarField, scalar>("Dp"); + + gradient() = (phi0p - phip)/patch().magSf()/Dpp; + + fixedGradientFvPatchScalarField::updateCoeffs(); +} + + +void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write +( + Ostream& os +) const +{ + fvPatchScalarField::write(os); + writeEntryIfDifferent<word>(os, "phi0", "phi0", phi0Name_); + writeEntryIfDifferent<word>(os, "phi", "phi", phiName_); + writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_); + gradient().writeEntry("gradient", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchScalarField, + multiphaseFixedFluxPressureFvPatchScalarField + ); +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H new file mode 100644 index 0000000000000000000000000000000000000000..f1699309b1f22672f7b11d7feaca8484684f0e27 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::multiphaseFixedFluxPressureFvPatchScalarField + +Description + Foam::multiphaseFixedFluxPressureFvPatchScalarField + +SourceFiles + multiphaseFixedFluxPressureFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H +#define multiphaseFixedFluxPressureFvPatchScalarFields_H + +#include "fvPatchFields.H" +#include "fixedGradientFvPatchFields.H" +#include "Switch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class multiphaseFixedFluxPressureFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class multiphaseFixedFluxPressureFvPatchScalarField +: + public fixedGradientFvPatchScalarField +{ + // Private data + + //- Name of the predicted flux transporting the field + word phi0Name_; + + //- Name of the flux transporting the field + word phiName_; + + //- Name of the density field used to normalise the mass flux + // if neccessary + word rhoName_; + + +public: + + //- Runtime type information + TypeName("multiphaseFixedFluxPressure"); + + + // Constructors + + //- Construct from patch and internal field + multiphaseFixedFluxPressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField<scalar, volMesh>& + ); + + //- Construct from patch, internal field and dictionary + multiphaseFixedFluxPressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const dictionary& + ); + + //- Construct by mapping given + // multiphaseFixedFluxPressureFvPatchScalarField onto a new patch + multiphaseFixedFluxPressureFvPatchScalarField + ( + const multiphaseFixedFluxPressureFvPatchScalarField&, + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + multiphaseFixedFluxPressureFvPatchScalarField + ( + const multiphaseFixedFluxPressureFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp<fvPatchScalarField> clone() const + { + return tmp<fvPatchScalarField> + ( + new multiphaseFixedFluxPressureFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + multiphaseFixedFluxPressureFvPatchScalarField + ( + const multiphaseFixedFluxPressureFvPatchScalarField&, + const DimensionedField<scalar, volMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp<fvPatchScalarField> clone + ( + const DimensionedField<scalar, volMesh>& iF + ) const + { + return tmp<fvPatchScalarField> + ( + new multiphaseFixedFluxPressureFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 338250415367f7093bd9ad6cf51b1c98dc9b9f29..2403597d9cf5837bc5d33f3111f429673fc33b0c 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -441,24 +441,15 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K Foam::multiphaseSystem::multiphaseSystem ( - const fvMesh& mesh, + const volVectorField& U, const surfaceScalarField& phi ) : - IOdictionary - ( - IOobject - ( - "transportProperties", - mesh.time().constant(), - mesh, - IOobject::MUST_READ_IF_MODIFIED, - IOobject::NO_WRITE - ) - ), - phases_(lookup("phases"), phaseModel::iNew(mesh)), + transportModel(U, phi), + + phases_(lookup("phases"), phaseModel::iNew(U.mesh())), - mesh_(mesh), + mesh_(U.mesh()), phi_(phi), alphas_ @@ -514,7 +505,6 @@ Foam::multiphaseSystem::multiphaseSystem // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // - Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const { PtrDictionary<phaseModel>::const_iterator iter = phases_.begin(); @@ -530,6 +520,21 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const } +Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const +{ + PtrDictionary<phaseModel>::const_iterator iter = phases_.begin(); + + tmp<volScalarField> tnu = iter()*iter().nu(); + + for (++iter; iter != phases_.end(); ++iter) + { + tnu() += iter()*iter().nu(); + } + + return tnu; +} + + Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm ( const phaseModel& phase @@ -770,19 +775,81 @@ void Foam::multiphaseSystem::solve() label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); - volScalarField& alpha = phases_.first(); - if (nAlphaSubCycles > 1) { dimensionedScalar totalDeltaT = runTime.deltaT(); + PtrList<volScalarField> alpha0s(phases_.size()); + PtrList<surfaceScalarField> phiSums(phases_.size()); + + int phasei = 0; + forAllIter(PtrDictionary<phaseModel>, phases_, iter) + { + phaseModel& phase = iter(); + volScalarField& alpha = phase; + + alpha0s.set + ( + phasei, + new volScalarField(alpha.oldTime()) + ); + + phiSums.set + ( + phasei, + new surfaceScalarField + ( + IOobject + ( + "phiSum" + alpha.name(), + runTime.timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0) + ) + ); + + phasei++; + } + for ( - subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles); + subCycleTime alphaSubCycle + ( + const_cast<Time&>(runTime), + nAlphaSubCycles + ); !(++alphaSubCycle).end(); ) { solveAlphas(); + + int phasei = 0; + forAllIter(PtrDictionary<phaseModel>, phases_, iter) + { + phiSums[phasei] += (runTime.deltaT()/totalDeltaT)*iter().phi(); + phasei++; + } + } + + phasei = 0; + forAllIter(PtrDictionary<phaseModel>, phases_, iter) + { + phaseModel& phase = iter(); + volScalarField& alpha = phase; + + phase.phi() = phiSums[phasei]; + + // Correct the time index of the field + // to correspond to the global time + alpha.timeIndex() = runTime.timeIndex(); + + // Reset the old-time field value + alpha.oldTime() = alpha0s[phasei]; + alpha.oldTime().timeIndex() = runTime.timeIndex(); + + phasei++; } } else diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H index 408fdcdafa84d248d215abdf99c91f7705ff2dd3..4c3e52e0ea548db8647d3a625e5a143133b7d638 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H @@ -61,7 +61,7 @@ namespace Foam class multiphaseSystem : - public IOdictionary + public transportModel { public: @@ -247,7 +247,7 @@ public: //- Construct from components multiphaseSystem ( - const fvMesh& mesh, + const volVectorField& U, const surfaceScalarField& phi ); @@ -274,6 +274,9 @@ public: //- Return the mixture density tmp<volScalarField> rho() const; + //- Return the mixture laminar viscosity + tmp<volScalarField> nu() const; + //- Return the virtual-mass coefficient for the given phase tmp<volScalarField> Cvm(const phaseModel& phase) const; @@ -305,6 +308,10 @@ public: //- Solve for the mixture phase-fractions void solve(); + //- Dummy correct + void correct() + {} + //- Read base transportProperties dictionary bool read(); }; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index 3bd26add308b724631a9d11776a862cfecabe893..1d4a9b1daaf69c319cff22efeea936059ab32ae0 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -41,6 +41,20 @@ phasei++; } + surfaceScalarField phi0 + ( + IOobject + ( + "phi0", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("phi0", dimArea*dimVelocity, 0) + ); + phasei = 0; forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter) { @@ -60,7 +74,8 @@ + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi()) ); mrfZones.relativeFlux(phase.phi()); - phase.phi() += rAlphaAUfs[phasei]*(g & mesh.Sf()); + surfaceScalarField pphi0("pphi0", phase.phi()); + pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf()); multiphaseSystem::dragModelTable::const_iterator dmIter = fluid.dragModels().begin(); @@ -99,13 +114,16 @@ phasej++; } - phase.phi() += + pphi0 += fvc::interpolate ((1.0/phase.rho())*rAUs[phasei]*(*dcIter())) *phi0s[phasej]; } } + phi0 += alphafs[phasei]*pphi0; + phase.phi() = pphi0; + phasei++; } @@ -132,13 +150,13 @@ phasei++; } Dp = mag(Dp); - adjustPhi(phi, U, p); + adjustPhi(phi0, U, p); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqnIncomp ( - fvc::div(phi) + fvc::div(phi0) - fvm::laplacian(Dp, p) ); @@ -169,6 +187,7 @@ } // dgdt = + // ( // pos(alpha2)*(pEqnComp2 & p)/rho2 // - pos(alpha1)*(pEqnComp1 & p)/rho1 diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H new file mode 100644 index 0000000000000000000000000000000000000000..5f8f3fb19a97401c51a89a032b116c3fc7c51754 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H @@ -0,0 +1,26 @@ +{ + const scalarField& V = mesh.V(); + + forAll(mesh.cellZones(), czi) + { + const labelList& cellLabels = mesh.cellZones()[czi]; + + forAllConstIter(PtrDictionary<phaseModel>, fluid.phases(), iter) + { + const volScalarField& alpha = iter(); + scalar phaseVolume = 0; + + forAll(cellLabels, cli) + { + label celli = cellLabels[cli]; + phaseVolume += alpha[celli]*V[celli]; + } + + reduce(phaseVolume, sumOp<scalar>()); + + Info<< alpha.name() + << " phase volume in zone " << mesh.cellZones()[czi].name() + << " = " << phaseVolume*1e6 << " ml " << endl; + } + } +} diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H index 9f256bd6891346da0be3844abc2882fb2079bb46..c8e3ccba5611cd6eced4965bad4ffaf3e79e8041 100644 --- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H +++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H @@ -12,18 +12,80 @@ ) ); - dimensionedScalar rho(mechanicalProperties.lookup("rho")); - dimensionedScalar rhoE(mechanicalProperties.lookup("E")); - dimensionedScalar nu(mechanicalProperties.lookup("nu")); + const dictionary& rhoDict(mechanicalProperties.subDict("rho")); + word rhoType(rhoDict.lookup("rho")); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimMass/dimVolume, 0.0) + ); + + if (rhoType == "rhoInf") + { + rho = rhoDict.lookup("rhoInf"); + } + + volScalarField rhoE + ( + IOobject + ( + "E", + runTime.timeName(0), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimMass/dimLength/sqr(dimTime), 0.0) + ); + + const dictionary& EDict(mechanicalProperties.subDict("E")); + word EType(EDict.lookup("E")); + if (EType == "EInf") + { + rhoE = EDict.lookup("EInf"); + } + + + volScalarField nu + ( + IOobject + ( + "nu", + runTime.timeName(0), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimless, 0.0) + ); + + const dictionary& nuDict(mechanicalProperties.subDict("nu")); + word nuType(nuDict.lookup("nu")); + + if (nuType == "nuInf") + { + nu = nuDict.lookup("nuInf"); + } Info<< "Normalising E : E/rho\n" << endl; - dimensionedScalar E = rhoE/rho; + volScalarField E = rhoE/rho; Info<< "Calculating Lame's coefficients\n" << endl; - dimensionedScalar mu = E/(2.0*(1.0 + nu)); - dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu)); - dimensionedScalar threeK = E/(1.0 - 2.0*nu); + volScalarField mu = E/(2.0*(1.0 + nu)); + volScalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu)); + volScalarField threeK = E/(1.0 - 2.0*nu); Switch planeStress(mechanicalProperties.lookup("planeStress")); @@ -38,7 +100,3 @@ { Info<< "Plane Strain\n" << endl; } - - Info<< "mu = " << mu.value() << " Pa/rho\n"; - Info<< "lambda = " << lambda.value() << " Pa/rho\n"; - Info<< "threeK = " << threeK.value() << " Pa/rho\n"; diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H index 12e10a607b4dc4effa5e04a0f1b59697df48e7d7..4d7bb43fc95905006ce40a9f6350601784f6b7ef 100644 --- a/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H +++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H @@ -1,46 +1,122 @@ - Info<< "Reading thermal properties\n" << endl; +Info<< "Reading thermal properties\n" << endl; - IOdictionary thermalProperties +IOdictionary thermalProperties +( + IOobject + ( + "thermalProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) +); + +Switch thermalStress(thermalProperties.lookup("thermalStress")); + +volScalarField threeKalpha +( + IOobject + ( + "threeKalpha", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0) +); + + +volScalarField DT +( + IOobject + ( + "DT", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(0, 2, -1 , 0, 0), 0.0) +); + + +if (thermalStress) +{ + volScalarField C ( IOobject ( - "thermalProperties", - runTime.constant(), + "C", + runTime.timeName(0), mesh, - IOobject::MUST_READ_IF_MODIFIED, + IOobject::READ_IF_PRESENT, IOobject::NO_WRITE - ) + ), + mesh, + dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0) ); - Switch thermalStress(thermalProperties.lookup("thermalStress")); + const dictionary& CDict(thermalProperties.subDict("C")); + word CType(CDict.lookup("C")); + if (CType == "CInf") + { + C = CDict.lookup("CInf"); + } + - dimensionedScalar threeKalpha + volScalarField rhoK ( - "threeKalpha", - dimensionSet(0, 2, -2 , -1, 0), - 0 + IOobject + ( + "k", + runTime.timeName(0), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(1, 1, -3 , -1, 0), 0.0) ); - dimensionedScalar DT + const dictionary& kDict(thermalProperties.subDict("k")); + word kType(kDict.lookup("k")); + if (kType == "kInf") + { + rhoK = kDict.lookup("kInf"); + } + + volScalarField alpha ( - "DT", - dimensionSet(0, 2, -1 , 0, 0), - 0 + IOobject + ( + "alpha", + runTime.timeName(0), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("0", dimensionSet(0, 0, 0 , -1, 0), 0.0) ); - if (thermalStress) + const dictionary& alphaDict(thermalProperties.subDict("alpha")); + word alphaType(alphaDict.lookup("alpha")); + + if (alphaType == "alphaInf") { - dimensionedScalar C(thermalProperties.lookup("C")); - dimensionedScalar rhoK(thermalProperties.lookup("k")); - dimensionedScalar alpha(thermalProperties.lookup("alpha")); + alpha = alphaDict.lookup("alphaInf"); + } - Info<< "Normalising k : k/rho\n" << endl; - dimensionedScalar k = rhoK/rho; + Info<< "Normalising k : k/rho\n" << endl; + volScalarField k = rhoK/rho; - Info<< "Calculating thermal coefficients\n" << endl; + Info<< "Calculating thermal coefficients\n" << endl; - threeKalpha = threeK*alpha; - DT.value() = (k/C).value(); + threeKalpha = threeK*alpha; + DT = k/C; - Info<< "threeKalpha = " << threeKalpha.value() << " Pa/rho\n"; - } +} diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C b/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C index 85f3c67b38e49e6a702a61cee1f4ac1747dba722..b24239c0f8b212615f752f0220dacbe16ae97773 100644 --- a/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C +++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/tractionDisplacement/tractionDisplacementFvPatchVectorField.C @@ -149,14 +149,20 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs() const dictionary& thermalProperties = db().lookupObject<IOdictionary>("thermalProperties"); - dimensionedScalar rho(mechanicalProperties.lookup("rho")); - dimensionedScalar rhoE(mechanicalProperties.lookup("E")); - dimensionedScalar nu(mechanicalProperties.lookup("nu")); - dimensionedScalar E = rhoE/rho; - dimensionedScalar mu = E/(2.0*(1.0 + nu)); - dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu)); - dimensionedScalar threeK = E/(1.0 - 2.0*nu); + const fvPatchField<scalar>& rho = + patch().lookupPatchField<volScalarField, scalar>("rho"); + + const fvPatchField<scalar>& rhoE = + patch().lookupPatchField<volScalarField, scalar>("E"); + + const fvPatchField<scalar>& nu = + patch().lookupPatchField<volScalarField, scalar>("nu"); + + scalarField E = rhoE/rho; + scalarField mu = E/(2.0*(1.0 + nu)); + scalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu)); + scalarField threeK = E/(1.0 - 2.0*nu); Switch planeStress(mechanicalProperties.lookup("planeStress")); @@ -166,7 +172,7 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs() threeK = E/(1.0 - nu); } - scalar twoMuLambda = (2*mu + lambda).value(); + scalarField twoMuLambda = (2*mu + lambda); vectorField n(patch().nf()); @@ -175,7 +181,7 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs() gradient() = ( - (traction_ + pressure_*n)/rho.value() + (traction_ + pressure_*n)/rho + twoMuLambda*fvPatchField<vector>::snGrad() - (n & sigmaD) )/twoMuLambda; @@ -183,13 +189,13 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs() if (thermalStress) { - dimensionedScalar alpha(thermalProperties.lookup("alpha")); - dimensionedScalar threeKalpha = threeK*alpha; + const fvPatchField<scalar>& threeKalpha= + patch().lookupPatchField<volScalarField, scalar>("threeKalpha"); const fvPatchField<scalar>& T = patch().lookupPatchField<volScalarField, scalar>("T"); - gradient() += n*threeKalpha.value()*T/twoMuLambda; + gradient() += n*threeKalpha*T/twoMuLambda; } fixedGradientFvPatchVectorField::updateCoeffs(); diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C index e56267b0b15c7dfddedb9b5ef8dbedae58528b00..69dc8d66ba4d22a2e9afad83cad279405b5602b7 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C @@ -89,7 +89,7 @@ void Foam::extrude2DMesh::setRefinement forAll(mesh_.points(), pointI) { point newPoint(mesh_.points()[pointI]); - newPoint[extrudeDir] = thickness; + newPoint[extrudeDir] += thickness; meshMod.addPoint ( diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties index 9ae921391c681c7f0f0087232e9095ad295e62a6..c92d02a0f7d43d03d8a7a03776cb74b317baaa3a 100644 --- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties +++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties @@ -15,11 +15,23 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -rho rho [ 1 -3 0 0 0 0 0 ] 7854; +rho +{ + rho rhoInf; + rhoInf rhoInf [ 1 -3 0 0 0 0 0 ] 7854; +} -nu nu [ 0 0 0 0 0 0 0 ] 0.3; +nu +{ + nu nuInf; + nuInf nuInf [ 0 0 0 0 0 0 0 ] 0.3; +} -E E [ 1 -1 -2 0 0 0 0 ] 2e+11; +E +{ + E EInf; + EInf EInf [ 1 -1 -2 0 0 0 0 ] 2e+11; +} planeStress yes; diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties index 9d09fa79015d53d74ddd4f78ce3837f80eb55bc6..f68549dbba9a3c552ce22a71ee8fb3336c4fedbe 100644 --- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties +++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties @@ -15,11 +15,23 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -C C [ 0 2 -2 -1 0 0 0 ] 434; +C +{ + C CInf; + CInf CInf [ 0 2 -2 -1 0 0 0 ] 434; +} -k k [ 1 1 -3 -1 0 0 0 ] 60.5; +k +{ + k kInf; + kInf kInf [ 1 1 -3 -1 0 0 0 ] 60.5; +} -alpha alpha [ 0 0 0 -1 0 0 0 ] 1.1e-05; +alpha +{ + alpha alphaInf; + alphaInf alphaInf [ 0 0 0 -1 0 0 0 ] 1.1e-05; +} thermalStress no; diff --git a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes index 0f667a4132194c9ab934c08a69d6082be2d2832a..a7c9b4c8d280067e885263f954fb6f113acc9ea1 100644 --- a/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes +++ b/tutorials/stressAnalysis/solidDisplacementFoam/plateHole/system/fvSchemes @@ -20,6 +20,11 @@ d2dt2Schemes default steadyState; } +ddtSchemes +{ + default Euler; +} + gradSchemes { default leastSquares;