View Issue Details

IDProjectCategoryView StatusLast Update
0002655OpenFOAMBugpublic2017-08-14 09:12
Reportertniemi Assigned Towill  
PrioritynormalSeveritymajorReproducibilityrandom
Status resolvedResolutionfixed 
Product Versiondev 
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.

Activities

tniemi

2017-08-07 14:02

reporter  

derivation.pdf (213,932 bytes)

tniemi

2017-08-07 14:02

reporter  

patch.diff (1,052 bytes)   
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 =
patch.diff (1,052 bytes)   

will

2017-08-07 18:09

manager   ~0008536

Thanks for the report. Resolved by commit 631bd489c1d509462a581707b546b4b677fa98b0.

tniemi

2017-08-14 08:14

reporter   ~0008569

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

will

2017-08-14 09:12

manager   ~0008570

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

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