diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H index a8c430061d..e4f9190cc5 100644 --- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H +++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H @@ -1,4 +1,4 @@ -volScalarField rAUrel = 1.0/UrelEqn().A(); +volScalarField rAUrel(1.0/UrelEqn().A()); Urel = rAUrel*UrelEqn().H(); if (pimple.nCorr() <= 1) diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C index e8dd45fa11..0c4177ea42 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -73,34 +73,23 @@ Foam::tmp Foam::GidaspowErgunWenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); // Wen and Yu (1966) - tmp tKWenYu = 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d(); - volScalarField& KWenYu = tKWenYu(); - - // Ergun - forAll (beta, cellj) - { - if (beta[cellj] <= 0.8) - { - KWenYu[cellj] = - 150.0*alpha_[cellj]*phaseb_.nu().value()*phaseb_.rho().value() - /sqr(beta[cellj]*phasea_.d().value()) - + 1.75*phaseb_.rho().value()*Ur[cellj] - /(beta[cellj]*phasea_.d().value()); - } - } - - return tKWenYu; + return + ( + pos(beta - 0.8) + *(0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d()) + + neg(beta - 0.8) + *( + 150.0*alpha_*phaseb_.nu()*phaseb_.rho()/(sqr(beta*phasea_.d())) + + 1.75*phaseb_.rho()*Ur/(beta*phasea_.d()) + ) + ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C index b9a7032759..8b8a284068 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,15 +72,12 @@ Foam::tmp Foam::GidaspowSchillerNaumann::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d(); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C index 0e69c8b58e..abbfcec43c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,15 +69,12 @@ Foam::tmp Foam::SchillerNaumann::K ) const { volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phaseb_.rho()*Ur/phasea_.d(); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C index d852bc0a8e..1b7780f8de 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,15 +70,11 @@ Foam::tmp Foam::SyamlalOBrien::K { volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField A(pow(beta, 4.14)); - volScalarField B(0.8*pow(beta, 1.28)); - - forAll (beta, celli) - { - if (beta[celli] > 0.85) - { - B[celli] = pow(beta[celli], 2.65); - } - } + volScalarField B + ( + neg(beta - 0.85)*(0.8*pow(beta, 1.28)) + + pos(beta - 0.85)*(pow(beta, 2.65)) + ); volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3))); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C index 00979db1e2..f34437f41e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,15 +72,11 @@ Foam::tmp Foam::WenYu::K volScalarField bp(pow(beta, -2.65)); volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3))); - volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re); - - forAll(Re, celli) - { - if (Re[celli] > 1000.0) - { - Cds[celli] = 0.44; - } - } + volScalarField Cds + ( + neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) + + pos(Re - 1000)*0.44 + ); return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d(); } diff --git a/applications/utilities/miscellaneous/foamDebugSwitches/foamDebugSwitches.C b/applications/utilities/miscellaneous/foamDebugSwitches/foamDebugSwitches.C index 8f0e53f4e2..7f84d96bc2 100644 --- a/applications/utilities/miscellaneous/foamDebugSwitches/foamDebugSwitches.C +++ b/applications/utilities/miscellaneous/foamDebugSwitches/foamDebugSwitches.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -61,7 +61,12 @@ int main(int argc, char *argv[]) if (args.optionFound("old") || args.optionFound("new")) { - dictionary controlDict(IFstream(findEtcFile("controlDict", true))()); + fileNameList controlDictFiles = findEtcFile("controlDict", true); + dictionary controlDict; + forAllReverse(controlDictFiles, cdfi) + { + controlDict.merge(dictionary(IFstream(controlDictFiles[cdfi])())); + } wordHashSet oldDebug ( diff --git a/etc/controlDict b/etc/controlDict index 9edaf1025e..352541fa75 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -24,7 +24,6 @@ Documentation "$WM_PROJECT_USER_DIR/html" "~OpenFOAM/html" "$WM_PROJECT_DIR/doc/Doxygen/html" - "$WM_PROJECT_DIR/doc/doxygen/html" ); doxySourceFileExts ( diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C index 48e4b3a4df..920d720cf5 100644 --- a/src/OSspecific/POSIX/POSIX.C +++ b/src/OSspecific/POSIX/POSIX.C @@ -57,6 +57,14 @@ Description #include +#ifdef USE_RANDOM +# include +# if INT_MAX != 2147483647 +# error "INT_MAX != 2147483647" +# error "The random number generator may not work!" +# endif +#endif + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::POSIX, 0); @@ -68,16 +76,19 @@ pid_t Foam::pid() return ::getpid(); } + pid_t Foam::ppid() { return ::getppid(); } + pid_t Foam::pgid() { return ::getpgrp(); } + bool Foam::env(const word& envName) { return ::getenv(envName.c_str()) != NULL; @@ -253,10 +264,11 @@ bool Foam::chDir(const fileName& dir) } -Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) +Foam::fileNameList Foam::findEtcFiles(const fileName& name, bool mandatory) { - // - // search for user files in + fileNameList results; + + // Search for user files in // * ~/.OpenFOAM/VERSION // * ~/.OpenFOAM // @@ -266,19 +278,17 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) fileName fullName = searchDir/FOAMversion/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } fullName = searchDir/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } } - - // - // search for group (site) files in + // Search for group (site) files in // * $WM_PROJECT_SITE/VERSION // * $WM_PROJECT_SITE // @@ -290,19 +300,18 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) fileName fullName = searchDir/FOAMversion/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } fullName = searchDir/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } } } else { - // // OR search for group (site) files in // * $WM_PROJECT_INST_DIR/site/VERSION // * $WM_PROJECT_INST_DIR/site @@ -313,20 +322,18 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) fileName fullName = searchDir/"site"/FOAMversion/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } fullName = searchDir/"site"/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } } } - - // - // search for other (shipped) files in + // Search for other (shipped) files in // * $WM_PROJECT_DIR/etc // searchDir = getEnv("WM_PROJECT_DIR"); @@ -335,24 +342,41 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) fileName fullName = searchDir/"etc"/name; if (isFile(fullName)) { - return fullName; + results.append(fullName); } } // Not found - // abort if the file is mandatory, otherwise return null - if (mandatory) + if (results.empty()) { - std::cerr - << "--> FOAM FATAL ERROR in Foam::findEtcFile() :" - " could not find mandatory file\n '" - << name.c_str() << "'\n\n" << std::endl; - ::exit(1); + // Abort if the file is mandatory, otherwise return null + if (mandatory) + { + std::cerr + << "--> FOAM FATAL ERROR in Foam::findEtcFiles() :" + " could not find mandatory file\n '" + << name.c_str() << "'\n\n" << std::endl; + ::exit(1); + } } - // Return null-constructed fileName rather than fileName::null - // to avoid cyclic dependencies in the construction of globals - return fileName(); + // Return list of matching paths or empty list if none found + return results; +} + + +Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory) +{ + fileNameList results(findEtcFiles(name, mandatory)); + + if (results.size()) + { + return results[0]; + } + else + { + return fileName(); + } } @@ -890,7 +914,6 @@ bool Foam::mvBak(const fileName& src, const std::string& ext) } - // Remove a file, returning true if successful otherwise false bool Foam::rm(const fileName& file) { @@ -1221,4 +1244,34 @@ Foam::fileNameList Foam::dlLoaded() } +void Foam::osRandomSeed(const label seed) +{ +#ifdef USE_RANDOM + srandom((unsigned int)seed); +#else + srand48(seed); +#endif +} + + +Foam::label Foam::osRandomInteger() +{ +#ifdef USE_RANDOM + return random(); +#else + return lrand48(); +#endif +} + + +Foam::scalar Foam::osRandomDouble() +{ +#ifdef USE_RANDOM + return (scalar)random()/INT_MAX; +#else + return drand48(); +#endif +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 1f1e25341e..3e52499e25 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -1,6 +1,6 @@ global/global.Cver -global/constants/constants.C -global/constants/dimensionedConstants.C +/* global/constants/constants.C in global.Cver */ +/* global/constants/dimensionedConstants.C in global.Cver */ global/argList/argList.C global/clock/clock.C @@ -136,7 +136,7 @@ $(StringStreams)/StringStreamsPrint.C Pstreams = $(Streams)/Pstreams $(Pstreams)/UIPstream.C $(Pstreams)/IPstream.C -$(Pstreams)/UPstream.C +/* $(Pstreams)/UPstream.C in global.Cver */ $(Pstreams)/UPstreamCommsStruct.C $(Pstreams)/Pstream.C $(Pstreams)/UOPstream.C @@ -181,7 +181,7 @@ $(IOobject)/IOobjectReadHeader.C $(IOobject)/IOobjectWriteHeader.C regIOobject = db/regIOobject -$(regIOobject)/regIOobject.C +/* $(regIOobject)/regIOobject.C in global.Cver */ $(regIOobject)/regIOobjectRead.C $(regIOobject)/regIOobjectWrite.C diff --git a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C index bab7ac3901..34dd4e80d1 100644 --- a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C +++ b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C @@ -197,8 +197,25 @@ Foam::functionEntries::codeStream::getFunction } } - // all processes must wait for compile to finish - reduce(create, orOp()); + //- We don't know whether this code was from IOdictionary + // (possibly read on master only) or from e.g. Field so cannot + // decide here. + //// all processes must wait for compile to finish - except if this + //// file is only read on the master + //bool masterOnly = + // ( + // regIOobject::fileModificationChecking + // == regIOobject::timeStampMaster + // ) + // || ( + // regIOobject::fileModificationChecking + // == regIOobject::inotifyMaster + // ); + // + //if (!masterOnly) + //{ + reduce(create, orOp()); + //} if (isA(topDict(parentDict))) { diff --git a/src/OpenFOAM/global/debug/debug.C b/src/OpenFOAM/global/debug/debug.C index b0c696987a..2cd9c861fe 100644 --- a/src/OpenFOAM/global/debug/debug.C +++ b/src/OpenFOAM/global/debug/debug.C @@ -75,10 +75,15 @@ Foam::dictionary& Foam::debug::controlDict() { if (!controlDictPtr_) { - controlDictPtr_ = new dictionary - ( - IFstream(findEtcFile("controlDict", true))() - ); + fileNameList controlDictFiles = findEtcFiles("controlDict", true); + controlDictPtr_ = new dictionary(); + forAllReverse(controlDictFiles, cdfi) + { + controlDictPtr_->merge + ( + dictionary(IFstream(controlDictFiles[cdfi])()) + ); + } } return *controlDictPtr_; diff --git a/src/OpenFOAM/global/global.Cver b/src/OpenFOAM/global/global.Cver index ef9a2ca9d2..3dc04c8b98 100644 --- a/src/OpenFOAM/global/global.Cver +++ b/src/OpenFOAM/global/global.Cver @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -65,6 +65,22 @@ bool Foam::JobInfo::constructed(false); #include "debug.C" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Read file modification checking switches + +#include "regIOobject.C" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Read parallel communication switches + +#include "UPstream.C" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Read constants + +#include "constants.C" +#include "dimensionedConstants.C" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Read and set cell models diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H index f01ad29647..119c669b5f 100644 --- a/src/OpenFOAM/include/OSspecific.H +++ b/src/OpenFOAM/include/OSspecific.H @@ -93,7 +93,7 @@ fileName cwd(); // else return false bool chDir(const fileName& dir); -//- Search for a file from user/group/shipped directories. +//- Search for files from user/group/shipped directories. // The search scheme allows for version-specific and // version-independent files using the following hierarchy: // - \b user settings: @@ -108,7 +108,14 @@ bool chDir(const fileName& dir); // - \b other (shipped) settings: // - $WM_PROJECT_DIR/etc/ // -// \return The full path name or fileName() if the name cannot be found +// \return The list of full paths of all the matching files or +// an empty list if the name cannot be found. +// Optionally abort if the file cannot be found +fileNameList findEtcFiles(const fileName&, bool mandatory=false); + +//- Search for a file using findEtcFiles. +// \return The full path name of the first file found which in the +// search hierarchy or an empty fileName if the name cannot be found. // Optionally abort if the file cannot be found fileName findEtcFile(const fileName&, bool mandatory=false); @@ -200,6 +207,18 @@ bool dlSymFound(void* handle, const std::string& symbol); fileNameList dlLoaded(); +// Low level random numbers. Use Random class instead. + +//- Seed random number generator. +void osRandomSeed(const label seed); + +//- Return random integer (uniform distribution between 0 and 2^31) +label osRandomInteger(); + +//- Return random double precision (uniform distribution between 0 and 1) +scalar osRandomDouble(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C index 704e5ef118..e2496ee648 100644 --- a/src/OpenFOAM/primitives/random/Random/Random.C +++ b/src/OpenFOAM/primitives/random/Random/Random.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "Random.H" +#include "OSspecific.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -37,11 +38,6 @@ namespace Foam # error "The random number generator may not work!" #endif -#ifdef USE_RANDOM -# include -#else -# include -#endif // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -57,22 +53,13 @@ Random::Random(const label seed) Seed = 1; } -# ifdef USE_RANDOM - srandom((unsigned int)Seed); -# else - srand48(Seed); -# endif - + osRandomSeed(Seed); } int Random::bit() { -# ifdef USE_RANDOM - if (random() > INT_MAX/2) -# else - if (lrand48() > INT_MAX/2) -# endif + if (osRandomInteger() > INT_MAX/2) { return 1; } @@ -85,11 +72,7 @@ int Random::bit() scalar Random::scalar01() { -# ifdef USE_RANDOM - return (scalar)random()/INT_MAX; -# else - return drand48(); -# endif + return osRandomDouble(); } @@ -140,11 +123,7 @@ tensor Random::tensor01() label Random::integer(const label lower, const label upper) { -# ifdef USE_RANDOM - return lower + (random() % (upper+1-lower)); -# else - return lower + (lrand48() % (upper+1-lower)); -# endif + return lower + (osRandomInteger() % (upper+1-lower)); } diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C index b48c1f4e63..ac0c270503 100644 --- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C +++ b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,7 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "cachedRandom.H" -#include +#include "OSspecific.H" #if INT_MAX != 2147483647 # error "INT_MAX != 2147483647" @@ -37,7 +37,7 @@ Foam::scalar Foam::cachedRandom::scalar01() { if (sampleI_ < 0) { - return drand48(); + return osRandomDouble(); } if (sampleI_ == samples_.size() - 1) @@ -76,7 +76,7 @@ Foam::cachedRandom::cachedRandom(const label seed, const label count) } // Initialise samples - srand48(seed_); + osRandomSeed(seed_); forAll(samples_, i) { samples_[i] = drand48(); @@ -98,7 +98,7 @@ Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset) ) << "Copy constructor called, but samples not being cached. " << "This may lead to non-repeatable behaviour" << endl; - srand48(seed_); + osRandomSeed(seed_); } else if (reset) { diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C index db585345f1..260fb328dc 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C @@ -59,10 +59,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField const DimensionedField& iF ) : - fixedValueFvPatchScalarField(p, iF), - phiName_("phi"), - rhoName_("rho"), - p0_(p.size(), 0.0), + totalPressureFvPatchScalarField(p, iF), fanCurve_(), direction_(ffdOut) {} @@ -76,10 +73,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField const fvPatchFieldMapper& mapper ) : - fixedValueFvPatchScalarField(ptf, p, iF, mapper), - phiName_(ptf.phiName_), - rhoName_(ptf.rhoName_), - p0_(ptf.p0_, mapper), + totalPressureFvPatchScalarField(ptf, p, iF, mapper), fanCurve_(ptf.fanCurve_), direction_(ptf.direction_) {} @@ -92,10 +86,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField const dictionary& dict ) : - fixedValueFvPatchScalarField(p, iF), - phiName_(dict.lookupOrDefault("phi", "phi")), - rhoName_(dict.lookupOrDefault("rho", "rho")), - p0_("p0", dict, p.size()), + totalPressureFvPatchScalarField(p, iF), fanCurve_(dict), direction_(fanFlowDirectionNames_.read(dict.lookup("direction"))) { @@ -109,10 +100,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField const fanPressureFvPatchScalarField& pfopsf ) : - fixedValueFvPatchScalarField(pfopsf), - phiName_(pfopsf.phiName_), - rhoName_(pfopsf.rhoName_), - p0_(pfopsf.p0_), + totalPressureFvPatchScalarField(pfopsf), fanCurve_(pfopsf.fanCurve_), direction_(pfopsf.direction_) {} @@ -124,10 +112,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField const DimensionedField& iF ) : - fixedValueFvPatchScalarField(pfopsf, iF), - phiName_(pfopsf.phiName_), - rhoName_(pfopsf.rhoName_), - p0_(pfopsf.p0_), + totalPressureFvPatchScalarField(pfopsf, iF), fanCurve_(pfopsf.fanCurve_), direction_(pfopsf.direction_) {} @@ -144,7 +129,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs() // Retrieve flux field const surfaceScalarField& phi = - db().lookupObject(phiName_); + db().lookupObject(phiName()); const fvsPatchField& phip = patch().patchField(phi); @@ -161,7 +146,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs() else if (phi.dimensions() == dimVelocity*dimArea*dimDensity) { const scalarField& rhop = - patch().lookupPatchField(rhoName_); + patch().lookupPatchField(rhoName()); aveFlowRate = dir*gSum(phip/rhop)/gSum(patch().magSf()); } else @@ -174,51 +159,23 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs() << exit(FatalError); } - // Normal flow through fan - if (aveFlowRate >= 0.0) - { - // Pressure drop for this flow rate - const scalar pdFan = fanCurve_(aveFlowRate); + // Pressure drop for this flow rate + const scalar pdFan = fanCurve_(max(aveFlowRate, 0.0)); - operator==(p0_ - dir*pdFan); - } - // Reverse flow - else - { - // Assume that fan has stalled if flow reversed - // i.e. apply dp for zero flow rate - const scalar pdFan = fanCurve_(0); - - // Flow speed across patch - scalarField Up = phip/(patch().magSf()); - - // Pressure drop associated withback flow = dynamic pressure - scalarField pdBackFlow = 0.5*magSqr(Up); - - if (phi.dimensions() == dimVelocity*dimArea*dimDensity) - { - const scalarField& rhop = - patch().lookupPatchField(rhoName_); - pdBackFlow /= rhop; - } - - operator==(p0_ - dir*(pdBackFlow + pdFan)); - } - - fixedValueFvPatchScalarField::updateCoeffs(); + totalPressureFvPatchScalarField::updateCoeffs + ( + p0() - dir*pdFan, + patch().lookupPatchField(UName()) + ); } void Foam::fanPressureFvPatchScalarField::write(Ostream& os) const { - fvPatchScalarField::write(os); - os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl; - os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl; + totalPressureFvPatchScalarField::write(os); fanCurve_.write(os); os.writeKeyword("direction") << fanFlowDirectionNames_[direction_] << token::END_STATEMENT << nl; - p0_.writeEntry("p0", os); - writeEntry("value", os); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H index 7c3d9ff411..6e99ca705a 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H @@ -25,7 +25,7 @@ Class Foam::fanPressureFvPatchScalarField Description - Assigns pressure inlet or outlet condition for a fan. + Assigns pressure inlet or outlet total pressure condition for a fan. User specifies: - pressure drop vs volumetric flow rate table (fan curve) file name; @@ -56,8 +56,8 @@ Description \endverbatim See Also - Foam::interpolationTable and - Foam::timeVaryingFlowRateInletVelocityFvPatchVectorField + Foam::totalPressureFvPatchScalarField and + Foam::interpolationTable SourceFiles fanPressureFvPatchScalarField.C @@ -67,8 +67,7 @@ SourceFiles #ifndef fanPressureFvPatchScalarField_H #define fanPressureFvPatchScalarField_H -#include "fvPatchFields.H" -#include "fixedValueFvPatchFields.H" +#include "totalPressureFvPatchScalarField.H" #include "interpolationTable.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -82,21 +81,10 @@ namespace Foam class fanPressureFvPatchScalarField : - public fixedValueFvPatchScalarField + public totalPressureFvPatchScalarField { - // Private data - //- Name of the flux transporting the field - word phiName_; - - //- Name of the density field - word rhoName_; - - //- Total pressure - scalarField p0_; - - //- Tabulated fan curve - interpolationTable fanCurve_; +public: //- Fan flow direction enum fanFlowDirection @@ -107,6 +95,13 @@ class fanPressureFvPatchScalarField static const NamedEnum fanFlowDirectionNames_; +private: + + // Private data + + //- Tabulated fan curve + interpolationTable fanCurve_; + //- Direction of flow through the fan relative to patch fanFlowDirection direction_; diff --git a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C index 9f3e49576b..10b284fa64 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C @@ -112,7 +112,7 @@ void Foam::rotatingTotalPressureFvPatchScalarField::updateCoeffs() + rotationVelocity ); - totalPressureFvPatchScalarField::updateCoeffs(Up); + totalPressureFvPatchScalarField::updateCoeffs(p0(), Up); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H index 087bf70b4a..01d9a73ad4 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C index 7defeec931..255a42d8c5 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C @@ -153,7 +153,11 @@ void Foam::totalPressureFvPatchScalarField::rmap } -void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) +void Foam::totalPressureFvPatchScalarField::updateCoeffs +( + const scalarField& p0p, + const vectorField& Up +) { if (updated()) { @@ -165,7 +169,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) if (psiName_ == "none" && rhoName_ == "none") { - operator==(p0_ - 0.5*(1.0 - pos(phip))*magSqr(Up)); + operator==(p0p - 0.5*(1.0 - pos(phip))*magSqr(Up)); } else if (rhoName_ == "none") { @@ -178,7 +182,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) operator== ( - p0_ + p0p /pow ( (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)), @@ -188,7 +192,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) } else { - operator==(p0_/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up))); + operator==(p0p/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up))); } } else if (psiName_ == "none") @@ -196,7 +200,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) const fvPatchField& rho = patch().lookupPatchField(rhoName_); - operator==(p0_ - 0.5*rho*(1.0 - pos(phip))*magSqr(Up)); + operator==(p0p - 0.5*rho*(1.0 - pos(phip))*magSqr(Up)); } else { @@ -220,7 +224,11 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up) void Foam::totalPressureFvPatchScalarField::updateCoeffs() { - updateCoeffs(patch().lookupPatchField(UName_)); + updateCoeffs + ( + p0(), + patch().lookupPatchField(UName()) + ); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H index 9aa386fe83..1eceb87faa 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -157,6 +157,45 @@ public: return UName_; } + //- Return the name of the flux field + const word& phiName() const + { + return phiName_; + } + + //- Return reference to the name of the flux field + // to allow adjustment + word& phiName() + { + return phiName_; + } + + //- Return the name of the density field + const word& rhoName() const + { + return rhoName_; + } + + //- Return reference to the name of the density field + // to allow adjustment + word& rhoName() + { + return rhoName_; + } + + //- Return the name of the compressibility field + const word& psiName() const + { + return psiName_; + } + + //- Return reference to the name of the compressibility field + // to allow adjustment + word& psiName() + { + return psiName_; + } + //- Return the heat capacity ratio scalar gamma() const { @@ -201,8 +240,12 @@ public: // Evaluation functions //- Update the coefficients associated with the patch field - // using the given patch velocity field - virtual void updateCoeffs(const vectorField& Up); + // using the given patch total pressure and velocity fields + virtual void updateCoeffs + ( + const scalarField& p0p, + const vectorField& Up + ); //- Update the coefficients associated with the patch field virtual void updateCoeffs(); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C index 812296e6f5..5580d3cb05 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C @@ -214,7 +214,8 @@ void Foam::PairSpringSliderDashpot::evaluatePair rHat_AB *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB)); - // Cohesion force + // Cohesion force, energy density multiplied by the area of + // particle-particle overlap if (cohesion_) { fN_AB += diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C index 9ae09dd98f..946ec3adb8 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -76,7 +76,8 @@ void Foam::WallLocalSpringSliderDashpot::evaluateWall typename CloudType::parcelType& p, const point& site, const WallSiteData& data, - scalar pREff + scalar pREff, + bool cohesion ) const { // wall patch index @@ -88,14 +89,18 @@ void Foam::WallLocalSpringSliderDashpot::evaluateWall scalar alpha = alpha_[wPI]; scalar b = b_[wPI]; scalar mu = mu_[wPI]; + scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI]; + cohesion = cohesion && cohesion_[wPI]; vector r_PW = p.position() - site; vector U_PW = p.U() - data.wallData(); - scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0); + scalar r_PW_mag = mag(r_PW); - vector rHat_PW = r_PW/(mag(r_PW) + VSMALL); + scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0); + + vector rHat_PW = r_PW/(r_PW_mag + VSMALL); scalar kN = (4.0/3.0)*sqrt(pREff)*Estar; @@ -105,6 +110,16 @@ void Foam::WallLocalSpringSliderDashpot::evaluateWall rHat_PW *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW)); + // Cohesion force, energy density multiplied by the area of wall/particle + // overlap + if (cohesion) + { + fN_PW += + -cohesionEnergyDensity + *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag)) + *rHat_PW; + } + p.f() += fN_PW; vector USlip_PW = @@ -168,6 +183,8 @@ Foam::WallLocalSpringSliderDashpot::WallLocalSpringSliderDashpot alpha_(), b_(), mu_(), + cohesionEnergyDensity_(), + cohesion_(), patchMap_(), maxEstarIndex_(-1), collisionResolutionSteps_ @@ -212,6 +229,8 @@ Foam::WallLocalSpringSliderDashpot::WallLocalSpringSliderDashpot alpha_.setSize(nWallPatches); b_.setSize(nWallPatches); mu_.setSize(nWallPatches); + cohesionEnergyDensity_.setSize(nWallPatches); + cohesion_.setSize(nWallPatches); scalar maxEstar = -GREAT; @@ -238,6 +257,13 @@ Foam::WallLocalSpringSliderDashpot::WallLocalSpringSliderDashpot mu_[wPI] = readScalar(patchCoeffDict.lookup("mu")); + cohesionEnergyDensity_[wPI] = readScalar + ( + patchCoeffDict.lookup("cohesionEnergyDensity") + ); + + cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > VSMALL); + if (Estar_[wPI] > maxEstar) { maxEstarIndex_ = wPI; @@ -325,20 +351,22 @@ void Foam::WallLocalSpringSliderDashpot::evaluateWall p, flatSitePoints[siteI], flatSiteData[siteI], - pREff + pREff, + true ); } forAll(sharpSitePoints, siteI) { - // Treating sharp sites like flat sites + // Treating sharp sites like flat sites, except suppress cohesion evaluateWall ( p, sharpSitePoints[siteI], sharpSiteData[siteI], - pREff + pREff, + false ); } } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H index 700ce321ab..c4ad801b13 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H @@ -65,6 +65,12 @@ class WallLocalSpringSliderDashpot //- Coefficient of friction in for tangential sliding scalarList mu_; + //- Cohesion energy density [J/m^3] + scalarList cohesionEnergyDensity_; + + // Switch cohesion on and off + boolList cohesion_; + //- Mapping the patch index to the model data labelList patchMap_; @@ -115,7 +121,8 @@ class WallLocalSpringSliderDashpot typename CloudType::parcelType& p, const point& site, const WallSiteData& data, - scalar pREff + scalar pREff, + bool cohesion ) const; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C index 5f14b8e274..1f547f0d5d 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,16 +77,19 @@ void Foam::WallSpringSliderDashpot::evaluateWall const point& site, const WallSiteData& data, scalar pREff, - scalar kN + scalar kN, + bool cohesion ) const { vector r_PW = p.position() - site; vector U_PW = p.U() - data.wallData(); - scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0); + scalar r_PW_mag = mag(r_PW); - vector rHat_PW = r_PW/(mag(r_PW) + VSMALL); + scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0); + + vector rHat_PW = r_PW/(r_PW_mag + VSMALL); scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag); @@ -94,6 +97,16 @@ void Foam::WallSpringSliderDashpot::evaluateWall rHat_PW *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW)); + // Cohesion force, energy density multiplied by the area of wall/particle + // overlap + if (cohesion) + { + fN_PW += + -cohesionEnergyDensity_ + *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag)) + *rHat_PW; + } + p.f() += fN_PW; vector USlip_PW = @@ -157,6 +170,11 @@ Foam::WallSpringSliderDashpot::WallSpringSliderDashpot alpha_(readScalar(this->coeffDict().lookup("alpha"))), b_(readScalar(this->coeffDict().lookup("b"))), mu_(readScalar(this->coeffDict().lookup("mu"))), + cohesionEnergyDensity_ + ( + readScalar(this->coeffDict().lookup("cohesionEnergyDensity")) + ), + cohesion_(false), collisionResolutionSteps_ ( readScalar @@ -183,6 +201,8 @@ Foam::WallSpringSliderDashpot::WallSpringSliderDashpot Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E); Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E)); + + cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL); } @@ -266,13 +286,14 @@ void Foam::WallSpringSliderDashpot::evaluateWall flatSitePoints[siteI], flatSiteData[siteI], pREff, - kN + kN, + cohesion_ ); } forAll(sharpSitePoints, siteI) { - // Treating sharp sites like flat sites + // Treating sharp sites like flat sites, except suppress cohesion evaluateWall ( @@ -280,7 +301,8 @@ void Foam::WallSpringSliderDashpot::evaluateWall sharpSitePoints[siteI], sharpSiteData[siteI], pREff, - kN + kN, + false ); } } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H index 10010c60c3..e6e93ddb24 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H @@ -65,6 +65,12 @@ class WallSpringSliderDashpot //- Coefficient of friction in for tangential sliding scalar mu_; + //- Cohesion energy density [J/m^3] + scalar cohesionEnergyDensity_; + + // Switch cohesion on and off + bool cohesion_; + //- The number of steps over which to resolve the minimum // harmonic approximation of the collision period scalar collisionResolutionSteps_; @@ -110,7 +116,8 @@ class WallSpringSliderDashpot const point& site, const WallSiteData& data, scalar pREff, - scalar kN + scalar kN, + bool cohesion ) const; diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C index 03c29aefd3..a7ec879794 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C @@ -744,7 +744,8 @@ void Foam::autoRefineDriver::doRefine 100 // maxIter ); - // Introduce baffles at surface intersections + // Introduce baffles at surface intersections. Remove sections unreachable + // from keepPoint. baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict); // Mesh is at its finest. Do optional zoning. diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index e9b137abc7..cac3aa685c 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -62,6 +62,7 @@ $(searchableSurface)/searchableSurfaces.C $(searchableSurface)/searchableSurfacesQueries.C $(searchableSurface)/searchableSurfaceWithGaps.C $(searchableSurface)/triSurfaceMesh.C +$(searchableSurface)/closedTriSurfaceMesh.C topoSets = sets/topoSets $(topoSets)/cellSet.C diff --git a/src/meshTools/searchableSurface/closedTriSurfaceMesh.C b/src/meshTools/searchableSurface/closedTriSurfaceMesh.C new file mode 100644 index 0000000000..14a505e356 --- /dev/null +++ b/src/meshTools/searchableSurface/closedTriSurfaceMesh.C @@ -0,0 +1,74 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "closedTriSurfaceMesh.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + +defineTypeNameAndDebug(closedTriSurfaceMesh, 0); +addToRunTimeSelectionTable(searchableSurface, closedTriSurfaceMesh, dict); + +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::closedTriSurfaceMesh::closedTriSurfaceMesh +( + const IOobject& io, + const triSurface& s +) +: + triSurfaceMesh(io, s) +{} + + +Foam::closedTriSurfaceMesh::closedTriSurfaceMesh(const IOobject& io) +: + triSurfaceMesh(io) +{} + + +Foam::closedTriSurfaceMesh::closedTriSurfaceMesh +( + const IOobject& io, + const dictionary& dict +) +: + triSurfaceMesh(io, dict) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::closedTriSurfaceMesh::~closedTriSurfaceMesh() +{} + + +// ************************************************************************* // diff --git a/src/meshTools/searchableSurface/closedTriSurfaceMesh.H b/src/meshTools/searchableSurface/closedTriSurfaceMesh.H new file mode 100644 index 0000000000..9290782719 --- /dev/null +++ b/src/meshTools/searchableSurface/closedTriSurfaceMesh.H @@ -0,0 +1,111 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. + \\/ 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 . + +Class + Foam::closedTriSurfaceMesh + +Description + A triSurfaceMesh where it is forced to check volumeTypes, used for surfaces + that are topologically non-manifold (small holes or multiple parts) but are + geometrically essentially closed + +SourceFiles + closedTriSurfaceMesh.C + +\*---------------------------------------------------------------------------*/ + +#ifndef closedTriSurfaceMesh_H +#define closedTriSurfaceMesh_H + +#include "triSurfaceMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class closedTriSurfaceMesh Declaration +\*---------------------------------------------------------------------------*/ + +class closedTriSurfaceMesh +: + public triSurfaceMesh +{ +private: + + //- Disallow default bitwise copy construct + closedTriSurfaceMesh(const closedTriSurfaceMesh&); + + //- Disallow default bitwise assignment + void operator=(const closedTriSurfaceMesh&); + + +public: + + //- Runtime type information + TypeName("closedTriSurfaceMesh"); + + + // Constructors + + //- Construct from triSurface + closedTriSurfaceMesh(const IOobject&, const triSurface&); + + //- Construct read. + closedTriSurfaceMesh(const IOobject& io); + + //- Construct from IO and dictionary (used by searchableSurface). + // Dictionary may contain a 'scale' entry (eg, 0.001: mm -> m) + closedTriSurfaceMesh + ( + const IOobject& io, + const dictionary& dict + ); + + + // Destructor + + virtual ~closedTriSurfaceMesh(); + + // Member Functions + + //- Whether supports volume type, forcing to true to force getVolumeType + // queries for this type + virtual bool hasVolumeType() const + { + return true; + } + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C index aa117c2794..ab26033eed 100644 --- a/src/meshTools/searchableSurface/triSurfaceMesh.C +++ b/src/meshTools/searchableSurface/triSurfaceMesh.C @@ -331,6 +331,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s) ), triSurface(s), tolerance_(indexedOctree::perturbTol()), + minQuality_(-1), maxTreeDepth_(10), surfaceClosed_(-1) {} @@ -375,6 +376,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io) ) ), tolerance_(indexedOctree::perturbTol()), + minQuality_(-1), maxTreeDepth_(10), surfaceClosed_(-1) {} @@ -422,6 +424,7 @@ Foam::triSurfaceMesh::triSurfaceMesh ) ), tolerance_(indexedOctree::perturbTol()), + minQuality_(-1), maxTreeDepth_(10), surfaceClosed_(-1) { @@ -444,6 +447,13 @@ Foam::triSurfaceMesh::triSurfaceMesh << tolerance_ << endl; } + // Have optional minimum quality for normal calculation + if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0) + { + Info<< searchableSurface::name() + << " : ignoring triangles with quality < " + << minQuality_ << " for normals calculation." << endl; + } // Have optional non-standard tree-depth to limit storage. if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0) @@ -804,22 +814,70 @@ void Foam::triSurfaceMesh::getNormal { normal.setSize(info.size()); - forAll(info, i) + if (minQuality_ >= 0) { - if (info[i].hit()) - { - label faceI = info[i].index(); - //- Cached: - //normal[i] = faceNormals()[faceI]; + // Make sure we don't use triangles with low quality since + // normal is not reliable. - //- Uncached - normal[i] = triSurface::operator[](faceI).normal(points()); - normal[i] /= mag(normal[i]) + VSMALL; - } - else + const triSurface& s = static_cast(*this); + const labelListList& faceFaces = s.faceFaces(); + + forAll(info, i) { - // Set to what? - normal[i] = vector::zero; + if (info[i].hit()) + { + label faceI = info[i].index(); + + scalar qual = s[faceI].tri(points()).quality(); + + if (qual < minQuality_) + { + // Search neighbouring triangles + const labelList& fFaces = faceFaces[faceI]; + + forAll(fFaces, j) + { + label nbrI = fFaces[j]; + scalar nbrQual = s[nbrI].tri(points()).quality(); + if (nbrQual > qual) + { + qual = nbrQual; + normal[i] = s[nbrI].normal(points()); + } + } + } + else + { + normal[i] = s[faceI].normal(points()); + } + normal[i] /= mag(normal[i]); + } + else + { + // Set to what? + normal[i] = vector::zero; + } + } + } + else + { + forAll(info, i) + { + if (info[i].hit()) + { + label faceI = info[i].index(); + //- Cached: + //normal[i] = faceNormals()[faceI]; + + //- Uncached + normal[i] = triSurface::operator[](faceI).normal(points()); + normal[i] /= mag(normal[i]) + VSMALL; + } + else + { + // Set to what? + normal[i] = vector::zero; + } } } } diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H index e1615ccc9b..19d2e78c88 100644 --- a/src/meshTools/searchableSurface/triSurfaceMesh.H +++ b/src/meshTools/searchableSurface/triSurfaceMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -31,6 +31,7 @@ Description - scale : scaling factor. - tolerance : relative tolerance for doing intersections (see triangle::intersection) + - minQuality: discard triangles with low quality when getting normal SourceFiles triSurfaceMesh.C @@ -70,6 +71,10 @@ private: //- Optional tolerance to use in searches scalar tolerance_; + //- Optional min triangle quality. Triangles below this get + // ignored for normal calculation + scalar minQuality_; + //- Optional max tree depth of octree label maxTreeDepth_; diff --git a/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C b/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C index 327c7b804a..1ea7ae253d 100644 --- a/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C +++ b/src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C @@ -27,6 +27,7 @@ License #include "dictionary.H" #include "Time.H" #include "IOobjectList.H" +#include "polyMesh.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -123,14 +124,26 @@ void Foam::partialWrite::write() else { // Delete all but marked objects + fileName dbDir; + if (isA(obr_)) + { + dbDir = dynamic_cast(obr_).dbDir(); + } + IOobjectList objects(obr_, obr_.time().timeName()); forAllConstIter(HashPtrTable, objects, iter) { if (!objectNames_.found(iter()->name())) { - const fileName f = obr_.time().timePath()/iter()->name(); - //Pout<< " rm " << f << endl; + const fileName f = + obr_.time().timePath() + /dbDir + /iter()->name(); + if (debug) + { + Pout<< " rm " << f << endl; + } rm(f); } } diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C index c001027d46..2f52c1b744 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -203,7 +203,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() // calculate the forces on the motion object from this data, then // update the positions - motion_.updatePosition(t.deltaTValue()); + motion_.updatePosition(t.deltaTValue(), t.deltaT0Value()); dictionary forcesDict; diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C index 1ab6cada87..889a5504d7 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -174,6 +174,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion() initialQ_(I), momentOfInertia_(diagTensor::one*VSMALL), mass_(VSMALL), + cDamp_(0.0), + aLim_(VGREAT), report_(false) {} @@ -190,6 +192,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion const point& initialCentreOfMass, const tensor& initialQ, const diagTensor& momentOfInertia, + scalar cDamp, + scalar aLim, bool report ) : @@ -211,6 +215,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion initialQ_(initialQ), momentOfInertia_(momentOfInertia), mass_(mass), + cDamp_(cDamp), + aLim_(aLim), report_(report) {} @@ -233,6 +239,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict) ), momentOfInertia_(dict.lookup("momentOfInertia")), mass_(readScalar(dict.lookup("mass"))), + cDamp_(dict.lookupOrDefault("accelerationDampingCoeff", 0.0)), + aLim_(dict.lookupOrDefault("accelerationLimit", VGREAT)), report_(dict.lookupOrDefault("report", false)) { addRestraints(dict); @@ -246,17 +254,19 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion const sixDoFRigidBodyMotion& sDoFRBM ) : - motionState_(sDoFRBM.motionState()), - restraints_(sDoFRBM.restraints()), - restraintNames_(sDoFRBM.restraintNames()), - constraints_(sDoFRBM.constraints()), - constraintNames_(sDoFRBM.constraintNames()), - maxConstraintIterations_(sDoFRBM.maxConstraintIterations()), - initialCentreOfMass_(sDoFRBM.initialCentreOfMass()), - initialQ_(sDoFRBM.initialQ()), - momentOfInertia_(sDoFRBM.momentOfInertia()), - mass_(sDoFRBM.mass()), - report_(sDoFRBM.report()) + motionState_(sDoFRBM.motionState_), + restraints_(sDoFRBM.restraints_), + restraintNames_(sDoFRBM.restraintNames_), + constraints_(sDoFRBM.constraints_), + constraintNames_(sDoFRBM.constraintNames_), + maxConstraintIterations_(sDoFRBM.maxConstraintIterations_), + initialCentreOfMass_(sDoFRBM.initialCentreOfMass_), + initialQ_(sDoFRBM.initialQ_), + momentOfInertia_(sDoFRBM.momentOfInertia_), + mass_(sDoFRBM.mass_), + cDamp_(sDoFRBM.cDamp_), + aLim_(sDoFRBM.aLim_), + report_(sDoFRBM.report_) {} @@ -358,7 +368,8 @@ void Foam::sixDoFRigidBodyMotion::addConstraints void Foam::sixDoFRigidBodyMotion::updatePosition ( - scalar deltaT + scalar deltaT, + scalar deltaT0 ) { // First leapfrog velocity adjust and motion part, required before @@ -366,9 +377,32 @@ void Foam::sixDoFRigidBodyMotion::updatePosition if (Pstream::master()) { - v() += 0.5*deltaT*a(); + vector aClip = a(); + scalar aMag = mag(aClip); - pi() += 0.5*deltaT*tau(); + if (aMag > SMALL) + { + aClip /= aMag; + } + + if (aMag > aLim_) + { + WarningIn + ( + "void Foam::sixDoFRigidBodyMotion::updatePosition" + "(" + "scalar deltaT, " + "scalar deltaT0" + ")" + ) + << "Limited acceleration " << a() + << " to " << aClip*aLim_ + << endl; + } + + v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_); + + pi() += 0.5*(1 - cDamp_)*deltaT0*tau(); // Leapfrog move part centreOfMass() += deltaT*v(); @@ -401,9 +435,33 @@ void Foam::sixDoFRigidBodyMotion::updateForce applyConstraints(deltaT); - v() += 0.5*deltaT*a(); + vector aClip = a(); + scalar aMag = mag(aClip); - pi() += 0.5*deltaT*tau(); + if (aMag > SMALL) + { + aClip /= aMag; + } + + if (aMag > aLim_) + { + WarningIn + ( + "void Foam::sixDoFRigidBodyMotion::updateForce" + "(" + "const vector& fGlobal, " + "const vector& tauGlobal, " + "scalar deltaT" + ")" + ) + << "Limited acceleration " << a() + << " to " << aClip*aLim_ + << endl; + } + + v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_); + + pi() += 0.5*(1 - cDamp_)*deltaT*tau(); if(report_) { diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H index a26c2378d1..ef4cbeb05a 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -116,6 +116,15 @@ class sixDoFRigidBodyMotion //- Mass of the body scalar mass_; + //- Acceleration damping coefficient. Modify applied acceleration: + // v1 = v0 + a*dt - cDamp*a*dt + // = v0 + dt*f*(1 - cDamp)/m + // Increases effective mass by 1/(1 - cDamp). + scalar cDamp_; + + //- Acceleration magnitude limit - clips large accelerations + scalar aLim_; + //- Switch to turn reporting of motion data on and off Switch report_; @@ -235,6 +244,8 @@ public: const point& initialCentreOfMass, const tensor& initialQ, const diagTensor& momentOfInertia, + scalar cDamp = 0.0, + scalar aLim = VGREAT, bool report = false ); @@ -260,10 +271,12 @@ public: void addConstraints(const dictionary& dict); //- First leapfrog velocity adjust and motion part, required - // before force calculation + // before force calculation. Takes old timestep for variable + // timestep cases. void updatePosition ( - scalar deltaT + scalar deltaT, + scalar deltaT0 ); //- Second leapfrog velocity adjust part, required after motion and diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C index 0312f396fe..aa2565e87c 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -40,6 +40,10 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const << momentOfInertia_ << token::END_STATEMENT << nl; os.writeKeyword("mass") << mass_ << token::END_STATEMENT << nl; + os.writeKeyword("accelerationDampingCoeff") + << cDamp_ << token::END_STATEMENT << nl; + os.writeKeyword("accelerationLimit") + << aLim_ << token::END_STATEMENT << nl; os.writeKeyword("report") << report_ << token::END_STATEMENT << nl; diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C index 7a045c04d4..fb168c8813 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C @@ -83,15 +83,20 @@ Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain ) const { vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0); + vector oldDir = refQ_ & refDir; + vector newDir = motion.orientation() & refDir; if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95) { // Directions getting close to the axis, change reference + refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1); - oldDir = refQ_ & refDir; - newDir = motion.orientation() & refDir; + + vector oldDir = refQ_ & refDir; + + vector newDir = motion.orientation() & refDir; } // Removing any axis component from oldDir and newDir and normalising diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C index 8581b3dd50..50cdb1460b 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C @@ -87,6 +87,7 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0); vector oldDir = refQ_ & refDir; + vector newDir = motion.orientation() & refDir; if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95) @@ -95,8 +96,9 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1); - oldDir = refQ_ & refDir; - newDir = motion.orientation() & refDir; + vector oldDir = refQ_ & refDir; + + vector newDir = motion.orientation() & refDir; } // Removing any axis component from oldDir and newDir and normalising diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C index 75771c7d84..5eae1b6bb8 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -147,7 +147,7 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() const polyMesh& mesh = this->dimensionedInternalField().mesh()(); const Time& t = mesh.time(); - motion_.updatePosition(t.deltaTValue()); + motion_.updatePosition(t.deltaTValue(), t.deltaT0Value()); vector gravity = vector::zero; diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict index 3e5b2dfbd0..fc5893222c 100644 --- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict +++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict @@ -16,19 +16,19 @@ FoamFile application simpleFoam; -startFrom latestTime; +startFrom latestTime; startTime 0; -stopAt nextWrite; +stopAt nextWrite; endTime 500; deltaT 1; -writeControl timeStep; +writeControl timeStep; -writeInterval 1; +writeInterval 1; purgeWrite 0; @@ -48,7 +48,7 @@ libs ( "libOpenFOAM.so" "libcompressibleTurbulenceModels.so" - "libincompressibleRASModels.so" + "libcompressibleRASModels.so" ); functions diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties index 7422e40c44..ccd5d89ef8 100644 --- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties +++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties @@ -112,6 +112,7 @@ subModels alpha 0.12; b 1.5; mu 0.43; + cohesionEnergyDensity 0; } frontAndBack { @@ -120,6 +121,7 @@ subModels alpha 0.12; b 1.5; mu 0.1; + cohesionEnergyDensity 0; } }; } diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties index 36448f7b30..3adb6ab771 100644 --- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties +++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties @@ -121,6 +121,7 @@ subModels alpha 0.12; b 1.5; mu 0.43; + cohesionEnergyDensity 0; } frontAndBack { @@ -129,6 +130,7 @@ subModels alpha 0.12; b 1.5; mu 0.1; + cohesionEnergyDensity 0; } }; }