View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0004193 | OpenFOAM | Patch | public | 2024-12-23 15:28 | 2024-12-23 16:09 |
| Reporter | wyldckat | Assigned To | |||
| Priority | normal | Severity | minor | Reproducibility | have not tried |
| Status | new | Resolution | open | ||
| Platform | GNU/Linux | OS | Other | OS Version | (please specify) |
| Product Version | 12 | ||||
| Summary | 0004193: Herschel-Bulkley implementation does not allow initializing with 'k' or 'tau0' values larger than ~1 | ||||
| Description | In the HerschelBulkley class, there is a division protection which relies on 'vSmall'. However, when dividing by 1e-308 (am using DP build), the numerator is multiplied by 1e308, which if large enough (e.g. 10 or 20), will result in overflow for double precision float. This overflow occurs at initialization, for example, with the following model settings: k 15.0; n 0.25; tau0 20.5; Attached is the proposed fix (this time with the correct line breaking bytes), to use 'small' instead of 'vSmall': - HerschelBulkley.C - the modified file to update 'src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C' in OpenFOAM 12. - HerschelBulkley_small_v1.patch - shows the modification. Attached is also a test case: - pitzDailySteadyHB.tar.gz It's based in 'incompressibleFluid/pitzDailySteady', adapted to use this rheology model, in which also changed the inlet vector just to try and stress test the laminar fluid flow profile. - In the file 'constant/momentumTransport', it's just a matter of changing the comment prefix in lines 41 and 42 to use coefficients which work without problems ('$HerschelBulkleyLow'), vs those which crash ('$HerschelBulkleyHigh'). | ||||
| Steps To Reproduce | tar -xf pitzDailySteadyHB.tar.gz cd pitzDailySteadyHB ./Allrun Results in SigFPE, using unmodified OpenFOAM 12 (and dev), as in attached file 'log'. | ||||
| Tags | No tags attached. | ||||
|
|
HerschelBulkley.C (3,255 bytes)
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2018-2024 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 "HerschelBulkley.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace laminarModels
{
namespace generalisedNewtonianViscosityModels
{
defineTypeNameAndDebug(HerschelBulkley, 0);
addToRunTimeSelectionTable
(
generalisedNewtonianViscosityModel,
HerschelBulkley,
dictionary
);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::
HerschelBulkley
(
const dictionary& viscosityProperties,
const Foam::viscosity& viscosity,
const volVectorField& U
)
:
strainRateViscosityModel(viscosityProperties, viscosity, U),
k_("k", dimKinematicViscosity, 0),
n_("n", dimless, 0),
tau0_("tau0", dimKinematicViscosity/dimTime, 0)
{
read(viscosityProperties);
correct();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::laminarModels::generalisedNewtonianViscosityModels::
HerschelBulkley::read
(
const dictionary& viscosityProperties
)
{
strainRateViscosityModel::read(viscosityProperties);
const dictionary& coeffs =
viscosityProperties.optionalSubDict(typeName + "Coeffs");
k_.read(coeffs);
n_.read(coeffs);
tau0_.read(coeffs);
return true;
}
Foam::tmp<Foam::volScalarField>
Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::
nu
(
const volScalarField& nu0,
const volScalarField& strainRate
) const
{
dimensionedScalar tone("tone", dimTime, 1.0);
dimensionedScalar rtone("rtone", dimless/dimTime, 1.0);
return
(
min
(
nu0,
(tau0_ + k_*rtone*pow(tone*strainRate, n_))
/max
(
strainRate,
dimensionedScalar ("small", dimless/dimTime, small)
)
)
);
}
// ************************************************************************* //
HerschelBulkley_small_v1.patch (1,029 bytes)
diff --git a/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C b/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
index c0511be812..4c267a03c3 100644
--- a/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
+++ b/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
@@ -108,7 +108,7 @@ nu
/max
(
strainRate,
- dimensionedScalar ("vSmall", dimless/dimTime, vSmall)
+ dimensionedScalar ("small", dimless/dimTime, small)
)
)
);
log (4,591 bytes)
Selecting generalised Newtonian model HerschelBulkley #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 ? in "/lib64/libc.so.6" #3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:? #4 void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:? #5 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:? #6 Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::nu(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:? #7 Foam::laminarModels::generalisedNewtonianViscosityModels::strainRateViscosityModel::correct() at ??:? #8 Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::HerschelBulkley(Foam::dictionary const&, Foam::viscosity const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:? #9 Foam::laminarModels::generalisedNewtonianViscosityModel::adddictionaryConstructorToTable<Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley>::New(Foam::dictionary const&, Foam::viscosity const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:? #10 Foam::laminarModels::generalisedNewtonianViscosityModel::New(Foam::dictionary const&, Foam::viscosity const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:? #11 Foam::laminarModel<Foam::incompressibleMomentumTransportModel>::adddictionaryConstructorToTable<Foam::laminarModels::generalisedNewtonian<Foam::incompressibleMomentumTransportModel> >::New(Foam::geometricOneField const&, Foam::geometricOneField const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::viscosity const&) at ??:? #12 Foam::laminarModel<Foam::incompressibleMomentumTransportModel>::New(Foam::geometricOneField const&, Foam::geometricOneField const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::viscosity const&) at ??:? #13 Foam::incompressibleMomentumTransportModel::adddictionaryConstructorToTable<Foam::laminarModel<Foam::incompressibleMomentumTransportModel> >::NewincompressibleMomentumTransportModel(Foam::geometricOneField const&, Foam::geometricOneField const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::viscosity const&) at ??:? #14 Foam::autoPtr<Foam::incompressibleMomentumTransportModel> Foam::momentumTransportModel::New<Foam::incompressibleMomentumTransportModel>(Foam::incompressibleMomentumTransportModel::alphaField const&, Foam::incompressibleMomentumTransportModel::rhoField const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::viscosity const&) at ??:? #15 Foam::incompressibleMomentumTransportModel::New(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::viscosity const&) at ??:? #16 Foam::solvers::incompressibleFluid::incompressibleFluid(Foam::fvMesh&) at ??:? #17 Foam::solver::addfvMeshConstructorToTable<Foam::solvers::incompressibleFluid>::New(Foam::fvMesh&) at ??:? #18 Foam::solver::New(Foam::word const&, Foam::fvMesh&) at ??:? #19 ? at ??:? #20 __libc_start_main in "/lib64/libc.so.6" #21 ? at ??:? Floating point exception (core dumped) |
|
|
Why choose small rather than rootVSmall? What is the smallest small which ensures no overflow from the division? |
|
|
Simple because 'small' seemed small enough and it was the one I knew without searching for all available options. I'll test 'rootVSmall' and let you know once the build update is complete (didn't check if it needs to build several template dependencies, simply ran 'Allwmake' at the root folder). |
|
|
Have confirmed that 'rootVSmall' is good enough. Please find in attachment the updated files: - HerschelBulkley.C - the modified file to update 'src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C' in OpenFOAM 12. - Assuming that MantisBT doesn't rename the file. - HerschelBulkley_small_v2.patch - shows the modification. HerschelBulkley-2.C (3,265 bytes)
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2018-2024 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 "HerschelBulkley.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace laminarModels
{
namespace generalisedNewtonianViscosityModels
{
defineTypeNameAndDebug(HerschelBulkley, 0);
addToRunTimeSelectionTable
(
generalisedNewtonianViscosityModel,
HerschelBulkley,
dictionary
);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::
HerschelBulkley
(
const dictionary& viscosityProperties,
const Foam::viscosity& viscosity,
const volVectorField& U
)
:
strainRateViscosityModel(viscosityProperties, viscosity, U),
k_("k", dimKinematicViscosity, 0),
n_("n", dimless, 0),
tau0_("tau0", dimKinematicViscosity/dimTime, 0)
{
read(viscosityProperties);
correct();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::laminarModels::generalisedNewtonianViscosityModels::
HerschelBulkley::read
(
const dictionary& viscosityProperties
)
{
strainRateViscosityModel::read(viscosityProperties);
const dictionary& coeffs =
viscosityProperties.optionalSubDict(typeName + "Coeffs");
k_.read(coeffs);
n_.read(coeffs);
tau0_.read(coeffs);
return true;
}
Foam::tmp<Foam::volScalarField>
Foam::laminarModels::generalisedNewtonianViscosityModels::HerschelBulkley::
nu
(
const volScalarField& nu0,
const volScalarField& strainRate
) const
{
dimensionedScalar tone("tone", dimTime, 1.0);
dimensionedScalar rtone("rtone", dimless/dimTime, 1.0);
return
(
min
(
nu0,
(tau0_ + k_*rtone*pow(tone*strainRate, n_))
/max
(
strainRate,
dimensionedScalar ("rootVSmall", dimless/dimTime, rootVSmall)
)
)
);
}
// ************************************************************************* //
HerschelBulkley_small_v2.patch (1,039 bytes)
diff --git a/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C b/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
index c0511be812..db94ef9f1f 100644
--- a/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
+++ b/src/MomentumTransportModels/momentumTransportModels/laminar/generalisedNewtonian/generalisedNewtonianViscosityModels/strainRateViscosityModels/HerschelBulkley/HerschelBulkley.C
@@ -108,7 +108,7 @@ nu
/max
(
strainRate,
- dimensionedScalar ("vSmall", dimless/dimTime, vSmall)
+ dimensionedScalar ("rootVSmall", dimless/dimTime, rootVSmall)
)
)
);
|
| Date Modified | Username | Field | Change |
|---|---|---|---|
| 2024-12-23 15:28 | wyldckat | New Issue | |
| 2024-12-23 15:28 | wyldckat | File Added: HerschelBulkley.C | |
| 2024-12-23 15:28 | wyldckat | File Added: HerschelBulkley_small_v1.patch | |
| 2024-12-23 15:28 | wyldckat | File Added: pitzDailySteadyHB.tar.gz | |
| 2024-12-23 15:28 | wyldckat | File Added: log | |
| 2024-12-23 15:42 | henry | Note Added: 0013492 | |
| 2024-12-23 15:56 | wyldckat | Note Added: 0013493 | |
| 2024-12-23 16:09 | wyldckat | Note Added: 0013494 | |
| 2024-12-23 16:09 | wyldckat | File Added: HerschelBulkley-2.C | |
| 2024-12-23 16:09 | wyldckat | File Added: HerschelBulkley_small_v2.patch |
