2017-08-21 07:26 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002655OpenFOAM[All Projects] Bugpublic2017-08-14 09:12
Reportertniemi 
Assigned Towill 
PrioritynormalSeveritymajorReproducibilityrandom
StatusresolvedResolutionfixed 
Product Versiondev 
Target VersionFixed in Version 
Summary0002655: ThermoParcel, numerical issue in calcHeatTransfer
DescriptionThe previous fix to Analytical.C (https://github.com/OpenFOAM/OpenFOAM-dev/commit/68f57d312655f45fcdb32cea3fdfb1df6a4c0b49) has exposed a problem in ThermoParcel.C calcHeatTransfer-function. Under certain conditions it can happen that the bp term becomes a large negative value. This produces a huge unphysical enthalpy source term which destroys the solution. The probability for this to happen is very small but it does happen occasionally. (Eq. once in a week of simulation in our test cases)

How the problem is triggered:

In line 299 in ThermoParcel.C, if (ap-T) is zero, bp will be divided with ROOTVSMALL. Mathematically from (ap-T)=0 it should follow that the numerator of bp is also zero. However, numerically the numerator can be nonzero leading to a large negative or positive bp. Large positive value is handled ok, but negative value is problematic.

Example values (no radiation):
Tc_=1238.109
T_=1238.79211
Sh=2.57176264e-5
htc=2960.9326
As=1.271486397570927e-08
->
ap-T = Tc_ + Sh/(As*htc)-T_ = 0
bp = 6.0*(Sh/As + htc*(Tc_ - T_))/ROOTVSMALL = -6.821210263296962e-12 / ROOTVSMALL = -6.821210263296962e+138

I have attached a patch where the bp term is written in mathematically equivalent form without the troublesome division. I have also attached a derivation for this patch.
TagsNo tags attached.
Attached Files
  • pdf file icon derivation.pdf (213,932 bytes) 2017-08-07 14:02
  • diff file icon patch.diff (1,052 bytes) 2017-08-07 14:02 -
    diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
    index 46454ec..0df5b2c 100644
    --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
    +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
    @@ -282,7 +282,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar As = this->areaS(d);
     
         scalar ap = Tc_ + Sh/(As*htc);
    -    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
    +    scalar bp = 6.0*htc/(rho*d*Cp_ + ROOTVSMALL);
         if (td.cloud().radiation())
         {
             tetIndices tetIs = this->currentTetIndices();
    @@ -294,9 +294,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
             scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_));
     
             ap += s/htc;
    -        bp += 6.0*s;
         }
    -    bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
     
         // Integrate to find the new parcel temperature
         IntegrationScheme<scalar>::integrationResult Tres =
    
    diff file icon patch.diff (1,052 bytes) 2017-08-07 14:02 +

-Relationships
+Relationships

-Notes

~0008536

will (manager)

Thanks for the report. Resolved by commit 631bd489c1d509462a581707b546b4b677fa98b0.

~0008569

tniemi (reporter)

We have been running with this patch for a while now with no issues. Could the patch be also added to 5.x?

~0008570

will (manager)

Yes, this seems sensible to put in 5.x, too. Commit 90ce63f6942d529d2aa6826b4729e78114ecdd4b.
+Notes

-Issue History
Date Modified Username Field Change
2017-08-07 14:01 tniemi New Issue
2017-08-07 14:02 tniemi File Added: derivation.pdf
2017-08-07 14:02 tniemi File Added: patch.diff
2017-08-07 17:41 will Assigned To => will
2017-08-07 17:41 will Status new => assigned
2017-08-07 18:09 will Status assigned => resolved
2017-08-07 18:09 will Resolution open => fixed
2017-08-07 18:09 will Note Added: 0008536
2017-08-14 08:14 tniemi Status resolved => feedback
2017-08-14 08:14 tniemi Resolution fixed => reopened
2017-08-14 08:14 tniemi Note Added: 0008569
2017-08-14 09:12 will Status feedback => resolved
2017-08-14 09:12 will Resolution reopened => fixed
2017-08-14 09:12 will Note Added: 0008570
+Issue History