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 |
| Reporter | tniemi | Assigned To | will | ||
| Priority | normal | Severity | major | Reproducibility | random |
| Status | resolved | Resolution | fixed | ||
| Product Version | dev | ||||
| Summary | 0002655: ThermoParcel, numerical issue in calcHeatTransfer | ||||
| Description | 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. | ||||
| Tags | No tags attached. | ||||
|
|
|
|
|
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 =
|
|
|
Thanks for the report. Resolved by commit 631bd489c1d509462a581707b546b4b677fa98b0. |
|
|
We have been running with this patch for a while now with no issues. Could the patch be also added to 5.x? |
|
|
Yes, this seems sensible to put in 5.x, too. Commit 90ce63f6942d529d2aa6826b4729e78114ecdd4b. |
| 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 |
