View Issue Details

IDProjectCategoryView StatusLast Update
0000153OpenFOAMBugpublic2011-03-02 10:48
Reporteruser131Assigned Tohenry  
Status resolvedResolutionfixed 
PlatformLinuxOSUbuntuOS Version10.10
Summary0000153: Herschel-Bulkley non-Newtonian viscosity model has term with sign error
DescriptionThe Herschel-Bulkley (HB) non-Newtonian viscosity model in OpenFOAM contains a term with a sign error. Outside of CFD, the HB model usually consists of a yield stress tao_y below which the shear rate sr is exactly zero, sr=0, and above which the stress tao is equal to tao_y plus a simple power law function of the shear rate, tao = tao_y + k*sr^n. The model in OpenFOAM allows for a transition to Newtonian behavior when the shear rate is very low, thereby eliminating the undesirable discontinuity in the function at zero stress. The Newtonian viscosity is nu0, and the transition occurs when the stress exceeds the threshold value given by the parameter tao0. Assuming that nu0 is not infinity, tao_y and tao0 are not the same but one can be calculated from the other by matching the Newtonian and non-Newtonian forms at the transition point resulting in: tao_y = tao0 - k*(tao0/nu0)^n. This correction to tao0 can be found in the source code for the HS model but INCORRECTLY with a plus sign rather than the minus sign in front of the k. In the code from the file HerschelBulkley.C, the model gives the effective viscosity nu by dividing the above expression for tao by a factor of sr:

    return (min(nu0_,(tau0_ + k_* rtone *( pow(tone * sr(), n_)
        + pow(tone*tau0_/nu0_,n_))) / (max(sr(), dimensionedScalar
        ("VSMALL", dimless/dimTime, VSMALL)))));

Note that the stress and viscosity are both given in kinematic form (i.e. divided by the density). The factors tone and rtone are both equal to unity but have dimensions of time and inverse time and are included to get the units correct for the result.

The CORRECTED VERSION with the sign change is:

    return (min(nu0_,(tau0_ + k_* rtone *( pow(tone * sr(), n_)
        - pow(tone*tau0_/nu0_,n_))) / (max(sr(), dimensionedScalar
        ("VSMALL", dimless/dimTime, VSMALL)))));

Note that when sr=tao0/nu0, the two pow() functions will now cancel out and give the expected result nu=tao0/sr. Without the indicated sign correction, tao0 has no obvious physical interpretation, as it is neither the stress at the transition point between the two forms nor is it the stress that would occur if the power law form was used as the shear rate approaches zero. The correction will have little effect in cases when the value of nu0 is set to a very large value to approximate infinity.

Steps To ReproduceN/A
TagsHerschel-Bulkley, non-Newtonian, rheology, viscosity



2011-03-02 10:48

manager   ~0000272

Thanks for the detailed bug-report.
Resolved by commit c4d7f891ae88c79ad22b6556a60701c90a30b698

Issue History

Date Modified Username Field Change
2011-03-01 23:41 user131 New Issue
2011-03-02 00:42 user131 Tag Attached: Herschel-Bulkley
2011-03-02 00:42 user131 Tag Attached: non-Newtonian
2011-03-02 00:42 user131 Tag Attached: rheology
2011-03-02 00:42 user131 Tag Attached: viscosity
2011-03-02 10:48 henry Note Added: 0000272
2011-03-02 10:48 henry Status new => resolved
2011-03-02 10:48 henry Resolution open => fixed
2011-03-02 10:48 henry Assigned To => henry