# New Viscosity Model - Best Practice for Temporary Values

 Register Blogs Members List Search Today's Posts Mark Forums Read

 August 23, 2016, 10:52 New Viscosity Model - Best Practice for Temporary Values #1 New Member   Join Date: Aug 2016 Posts: 4 Rep Power: 9 Hello FOAMers, currently I am implementing a temperature-dependant viscosity model, to be concise, the Carreau-WLF-Model, for polymer melts. It is very similiar to a Cross-WLF-model. (See equations 2.29 and 2.32 in this preview for more details: http://files.hanser.de/hanser/docs/2..._Leseprobe.pdf) I have referred to the following materials: http://www.cfd-online.com/Forums/ope...time-loop.html http://www.tfd.chalmers.se/~hani/kur...nFoam%20v2.pdf http://www.cfd-online.com/Forums/ope...mperature.html Given the fact that the equation is - while not complex - rather long, I tried to split it up across several lines. See this snippet from GermanCarreau.C. Code: ```// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // Foam::tmp Foam::viscosityModels::GermanCarreau::calcNu() const { const volScalarField& Tcurrent=U_.mesh().lookupObject("T"); Foam::dimensionedScalar Tsp1 = Te_+scalar(50); Foam::dimensionedScalar termA = (scalar(8.86)*(Tmeasure_ - Tsp1)) / (scalar(101.6)+(Tmeasure_ -Tsp1)); Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1)); Foam::dimensionedScalar alpha_shift = pow(scalar(10), (termA - termB)); Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) )); return (eta_shifted/rho_.value() ); }``` This fails to compile, resulting in the following errors: Code: ```viscosityModels/GermanCarreau/GermanCarreau.C:58:82: error: conversion from ‘Foam::tmp >’ to non-scalar type ‘Foam::dimensionedScalar {aka Foam::dimensioned}’ requested Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent.internalField() - Tsp1)) / (scalar(101.6)+(Tcurrent.internalField() -Tsp1)); ^ viscosityModels/GermanCarreau/GermanCarreau.C:60:58: error: conversion from ‘Foam::tmp >’ to non-scalar type ‘Foam::dimensionedScalar {aka Foam::dimensioned}’ requested Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) )); ^ viscosityModels/GermanCarreau/GermanCarreau.C:62:24: error: could not convert ‘Foam::operator/(const Foam::dimensioned&, const Foam::dimensioned&) [with Type = double](Foam::dimensioned((* &((const Foam::viscosityModels::GermanCarreau*)this)->Foam::viscosityModels::GermanCarreau::rho_.Foam::dimensioned::value())))’ from ‘Foam::dimensioned’ to ‘Foam::tmp >’ return (eta_shifted/rho_.value() ^ viscosityModels/GermanCarreau/GermanCarreau.C: In constructor ‘Foam::viscosityModels::GermanCarreau::GermanCarreau(const Foam::word&, const Foam::dictionary&, const volVectorField&, const surfaceScalarField&)’: viscosityModels/GermanCarreau/GermanCarreau.C:81:5: error: class ‘Foam::viscosityModels::GermanCarreau’ does not have any field named ‘GermanCarreauCoeffs_’ GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")), ^ viscosityModels/GermanCarreau/GermanCarreau.C:82:22: error: ‘GermanCarreauCoeffs_’ was not declared in this scope A_("A", dimless, GermanCarreauCoeffs_), ^ viscosityModels/GermanCarreau/GermanCarreau.C:83:21: error: expected ‘)’ before ‘GermanCarreauCoeffs_’ B_("B", dimless GermanCarreauCoeffs_),``` As you can see, Scalars and volScalarFields cannot be used together. As far as I understand, the transport equation is evaluated for each point in the flow field - so how would I go about describing my temporary variables in terms of OpenFOAM data types? How would you approach such a problem? Is my approach of compartmentalizing the equation fundamentally sound? Can you point to a implented transport model that utilizes a very long formula in a similar manner? Thank you very much! Attached, you'll find GermanCarreau.C and GermanCarreau.H: GermanCarreau.C Code: ```/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 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 . \*---------------------------------------------------------------------------*/ #include "GermanCarreau.H" #include "addToRunTimeSelectionTable.H" #include "surfaceFields.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { namespace viscosityModels { defineTypeNameAndDebug(GermanCarreau, 0); addToRunTimeSelectionTable ( viscosityModel, GermanCarreau, dictionary ); } } // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // Foam::tmp Foam::viscosityModels::GermanCarreau::calcNu() const { const volScalarField& Tcurrent=U_.mesh().lookupObject("T"); Foam::dimensionedScalar Tsp1 = Te_+scalar(50); Foam::dimensionedScalar termA = (scalar(8.86)*(Tmeasure_ - Tsp1)) / (scalar(101.6)+(Tmeasure_ -Tsp1)); Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1)); Foam::dimensionedScalar alpha_shift = pow(scalar(10), (termA - termB)); Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) )); return (eta_shifted/rho_.value() ); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::viscosityModels::GermanCarreau::GermanCarreau ( const word& name, const dictionary& viscosityProperties, const volVectorField& U, const surfaceScalarField& phi ) : viscosityModel(name, viscosityProperties, U, phi), GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")), A_("A", dimless, GermanCarreauCoeffs_), B_("B", dimless GermanCarreauCoeffs_), C_("C", dimless, GermanCarreauCoeffs_), Te_(GermanCarreauCoeffs_.lookup("Te"), rho_("rho", dimless, GermanCarreauCoeffs_), Tmeasure_("Tmeasure", dimless, GermanCarreauCoeffs_), pFactor_("pFactor", dimless, GermanCarreauCoeffs_), nu_ ( IOobject ( name, U_.time().timeName(), U_.db(), IOobject::NO_READ, IOobject::AUTO_WRITE ), calcNu() ) {} // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // bool Foam::viscosityModels::GermanCarreau::read ( const dictionary& viscosityProperties ) { viscosityModel::read(viscosityProperties); GermanCarreau_ = viscosityProperties.subDict(typeName + "Coeffs"); GermanCarreauCoeffs_.lookup("A") >> A_; GermanCarreauCoeffs_.lookup("B") >> B_; GermanCarreauCoeffs_.lookup("C") >> C_; GermanCarreauCoeffs_.lookup("Te") >> Te_; GermanCarreauCoeffs_.lookup("rho") >> rho_; GermanCarreauCoeffs_.lookup("Tmeasure") >> Tmeasure_; GermanCarreauCoeffs_.lookup("pFactor") >> pFactor_; return true; } // ************************************************************************* //``` GermanCarreau.H Code: ```/*---------------------------------------------------------------------------*\ ========= | \\ / 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 . Class Foam::viscosityModels::GermanCarreau Description An incompressible isothermal Carreau Three-factor viscosity model. SourceFiles GermanCarreau.C \*---------------------------------------------------------------------------*/ #ifndef GermanCarreau_H #define GermanCarreau_H #include "viscosityModel.H" #include "dimensionedScalar.H" #include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { namespace viscosityModels { /*---------------------------------------------------------------------------*\ Class GermanCarreau Declaration \*---------------------------------------------------------------------------*/ class GermanCarreau : public viscosityModel { // Private data dictionary GermanCarreau_; dimensionedScalar A_; dimensionedScalar B_; dimensionedScalar C_; dimensionedScalar rho_; dimensionedScalar Te_; dimensionedScalar Tmeasure_; dimensionedScalar pFactor_; dimensionedScalar Tsp2_; dimensionedScalar termA_; dimensionedScalar termB_; dimensionedScalar eta_shifted_; dimensionedScalar alpha_shift_; volScalarField nu_; // Private Member Functions //- Calculate and return the laminar viscosity tmp calcNu() const; public: //- Runtime type information TypeName("GermanCarreau"); // Constructors //- Construct from components GermanCarreau ( const word& name, const dictionary& viscosityProperties, const volVectorField& U, const surfaceScalarField& phi ); //- Destructor ~GermanCarreau() {} // Member Functions //- Return the laminar viscosity tmp nu() const { return nu_; } //- Return the laminar viscosity for patch tmp nu(const label patchi) const { return nu_.boundaryField()[patchi]; } //- Correct the laminar viscosity void correct() { nu_ = calcNu(); } //- Read transportProperties dictionary bool read(const dictionary& viscosityProperties); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace viscosityModels } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* //```

 August 24, 2016, 03:36 #2 Senior Member   Gerhard Holzinger Join Date: Feb 2012 Location: Austria Posts: 339 Rep Power: 28 Hello, in Line 58 you try to compute a dimensionedScalar with a formula containing a volScalarField. This is equivalent to compute a scalar directly from a vector. You can compute the magnitude of the vector or the inner product to gain a scalar, but otherwise you won't get a scalar out of a vector. The same logic applies to scalar fields. Code: `Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1));` Thus, in your model the variables termB, alpha_shift and eta_shifted need to be volScalarFields, unless you do something about Tcurrent, e.g. apply spatial averaging to make a scalar from the field.

 August 24, 2016, 03:54 #3 New Member   Join Date: Aug 2016 Posts: 4 Rep Power: 9 First off: Thank you a lot for your very useful feedback, Gerhard. There is - however - a simple question remaining: As far as I understand, the viscosity model is evaluated for each point the flow channel. Otherwise, shear thinning [that is, the reduction of the kinematic viscosity nu in areas of high shear rate] could not be observed. How would I modify Code: `const volScalarField& Tcurrent=U_.mesh().lookupObject("T");{` so that Tcurrent is not the entire volScalarField, but rather the temperature of the node for which nu is currently being calculated? Thanks again!

 August 24, 2016, 04:36 #4 Senior Member   Gerhard Holzinger Join Date: Feb 2012 Location: Austria Posts: 339 Rep Power: 28 If the variables termB, alpha_shift and eta_shifted are volScalarFields, then your formula is automatically evaluated for all cells of the mesh. There is no need to change the expression for Tcurrent.

 August 25, 2016, 09:18 #5 New Member   Join Date: Aug 2016 Posts: 4 Rep Power: 9 Hello Foamers, in the end I did manage to make it work. I used a hack described here: http://www.cfd-online.com/Forums/ope...tml#post595370 so I could avoid having to manually create new volScalarFields. Also, I left out the calculation of Tstandard. Now, Ts is defined in the transportProperties file. A short test revealed that the results are accurate to within ca. 10%, which I blame on shear heating and my godawful mesh. Attached you'll find the compileable files GermanCarreau.C and GermanCarreau.H. However, the material model only runs with a non-isothermal solver such as the one found here: http://www.dhcae-tools.com/docs/visousHeatingSolver.pdf I hope you'll find my results useful! Code: ```/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 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 . \*---------------------------------------------------------------------------*/ #include "GermanCarreau.H" #include "addToRunTimeSelectionTable.H" #include "surfaceFields.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { namespace viscosityModels { defineTypeNameAndDebug(GermanCarreau, 0); addToRunTimeSelectionTable ( viscosityModel, GermanCarreau, dictionary ); } } // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // Foam::tmp Foam::viscosityModels::GermanCarreau::calcNu() const { const volScalarField& Tcurrent=U_.mesh().lookupObject("T"); volScalarField alpha_shift =Tcurrent; //instead of constructing a new volScalarField, we "copy" the old one alpha_shift *= scalar(0.0); //we then set all values of the volScalarField to zero alpha_shift.dimensions().reset(dimless); //also, since T is a Value in Kelvin, we make the new Field dimensionsless volScalarField termA =Tcurrent; termA *= scalar(0.0); termA.dimensions().reset(dimless); volScalarField termB =Tcurrent; termB *= scalar(0.0); termB.dimensions().reset(dimless); volScalarField Tnew =Tcurrent; Tnew.dimensions().reset(dimless); //alpha_shift += scalar(1.2); //Tsp1.dimensions().reset(dimLength*dimLength/dimTime); //turns into a kinematic viscosity termB = (scalar(8.86)*(Tnew - Ts_)) / (scalar(101.6)+(Tnew -Ts_)); termA = (scalar(8.86)*(Tmeasure_ - Ts_)) / (scalar(101.6)+(Tmeasure_ -Ts_)); alpha_shift = pow(scalar(10), (termA - termB)); return (nuA_*alpha_shift)/(scalar(1) + pow(alpha_shift*B_*strainRate(), C_)); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::viscosityModels::GermanCarreau::GermanCarreau ( const word& name, const dictionary& viscosityProperties, const volVectorField& U, const surfaceScalarField& phi ) : viscosityModel(name, viscosityProperties, U, phi), GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")), nuA_("nuA", dimViscosity, GermanCarreauCoeffs_), B_("B", dimTime, GermanCarreauCoeffs_), C_("C", dimless, GermanCarreauCoeffs_), Ts_("Ts", dimless, GermanCarreauCoeffs_), Tmeasure_("Tmeasure", dimless, GermanCarreauCoeffs_), nu_ ( IOobject ( name, U_.time().timeName(), U_.db(), IOobject::NO_READ, IOobject::AUTO_WRITE ), calcNu() ) {} // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // bool Foam::viscosityModels::GermanCarreau::read ( const dictionary& viscosityProperties ) { viscosityModel::read(viscosityProperties); GermanCarreauCoeffs_ = viscosityProperties.subDict(typeName + "Coeffs"); GermanCarreauCoeffs_.lookup("nuA") >> nuA_; GermanCarreauCoeffs_.lookup("B") >> B_; GermanCarreauCoeffs_.lookup("C") >> C_; GermanCarreauCoeffs_.lookup("Ts") >> Ts_; GermanCarreauCoeffs_.lookup("Tmeasure") >> Tmeasure_; return true; } // ************************************************************************* //``` Code: ```/*---------------------------------------------------------------------------*\ ========= | \\ / 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 . Class Foam::viscosityModels::GermanCarreau Description An incompressible isothermal Carreau Three-factor viscosity model with WLF correction for temperatures SourceFiles GermanCarreau.C \*---------------------------------------------------------------------------*/ #ifndef GermanCarreau_H #define GermanCarreau_H #include "viscosityModel.H" #include "dimensionedScalar.H" #include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { namespace viscosityModels { /*---------------------------------------------------------------------------*\ Class GermanCarreau Declaration \*---------------------------------------------------------------------------*/ class GermanCarreau : public viscosityModel { // Private data dictionary GermanCarreauCoeffs_; dimensionedScalar nuA_; dimensionedScalar B_; dimensionedScalar C_; dimensionedScalar Ts_; dimensionedScalar Tmeasure_; volScalarField nu_; // Private Member Functions //- Calculate and return the laminar viscosity tmp calcNu() const; public: //- Runtime type information TypeName("GermanCarreau"); // Constructors //- Construct from components GermanCarreau ( const word& name, const dictionary& viscosityProperties, const volVectorField& U, const surfaceScalarField& phi ); //- Destructor ~GermanCarreau() {} // Member Functions //- Return the laminar viscosity tmp nu() const { return nu_; } //- Return the laminar viscosity for patch tmp nu(const label patchi) const { return nu_.boundaryField()[patchi]; } //- Correct the laminar viscosity void correct() { nu_ = calcNu(); } //- Read transportProperties dictionary bool read(const dictionary& viscosityProperties); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace viscosityModels } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* //``` vanq likes this.

 March 6, 2018, 02:41 Compressible thermophysical transport model #6 Member   Alberto Join Date: Sep 2013 Posts: 37 Rep Power: 12 Hello everyone, I have seen some tutorials on how to modify the (incompressible) viscosity models (Powerlaw, Casson, etc.) to make them temperature dependent. I am trying to modify a compressible thermo physical model, such as Sutherland model (for viscosity) which depends on Pressure and Temperature, but not on the strainRate. The strainRate() is directly available for incompressible viscosityModels, but it is not available for compressible thermo physical model, such as Sutherland. Does any body know how to make the Sutherland model strainRate dependent? in addition to the Pressure and Temperature dependency?. I am using OpenFoam 5.x in Windows by BlueCFD. Thanks in advanced

 March 16, 2018, 03:51 #7 New Member   Join Date: Aug 2016 Posts: 4 Rep Power: 9 Dear malv, while I'm not certain that it's going to work in the compressible flow model, you might try to calculate the strain rate directly from the flow field. See this useful thread: StrainRate definition 141 vs 15 which eventually arrives at: Code: `strainrate= sqrt(2) * mag(symm(fvc::grad(U)))`