#### View Issue Details

ID Project Category View Status Date Submitted Last Update 0002655 OpenFOAM Bug public 2017-08-07 14:01 2017-08-14 09:12 tniemi will normal major random resolved fixed dev 0002655: ThermoParcel, numerical issue in calcHeatTransfer The 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. No tags attached.

#### Activities

 2017-08-07 14:02 reporter derivation.pdf (213,932 bytes) 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::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::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::integrationResult Tres = ``` patch.diff (1,052 bytes) 2017-08-07 18:09 manager   ~0008536 Thanks for the report. Resolved by commit 631bd489c1d509462a581707b546b4b677fa98b0. 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? 2017-08-14 09:12 manager   ~0008570 Yes, this seems sensible to put in 5.x, too. Commit 90ce63f6942d529d2aa6826b4729e78114ecdd4b.

#### Issue History

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