基于hPolynomial热物理模型实现Cp多段多项式拟合过程中遇到的问题
-
从这个图来看，Cp是一个存在极大值的曲线。如果楼主能确认T的函数也是存在极值的话。按照 @xpqiu 的解释，温度求解使用的为牛顿迭代法。牛顿迭代法如果某一解位于极值附近，必然发散不会收敛，则导致温度越界。你可以试试别的迭代法。
Test=-144277 Tnew=-288653 F=9.35162 f=-1.35014e+06 Cp=9.35162 Test=-288653 Tnew=-433029 F=9.35162 f=-1.35014e+06 Cp=9.35162 Test=-433029 Tnew=-577405 F=9.35162 f=-1.35014e+06 Cp=9.35162 Test=-577405 Tnew=-721781 F=9.35162 f=-1.35014e+06 Cp=9.35162 看起来这些温度越界了如果确认是这个问题的话，我建议采用Ridder求解。非常简单，并且对于下图这种存在极值的曲线也可以求解，或许OpenFOAM基金会也会对你这个工作感兴趣。
个人浅见，最好确定一下是不是我说的这个路子。
-
-
@random_ran 多谢。目前试验的拟合还在180K以内，所以还没涉及到Cp梯度大的情况。
-
我最初看错了范围，我以为你的曲线里包含了200K附近的那个间断点，所以才会说是牛顿法不适合。现在发现你的温度范围是100-180K，这一段就不包含间断点了，所以牛顿法应该没问题。
从你上面提供的调试信息，我认为主要问题在于，当T < 100K 和 T > 180K 时，你的cp函数将无定义，return 值是无定义的。你看出现开始发散的地方：Test=100.516 Tnew=99.569 F=-1.34695e+06 f=-1.35014e+06 Cp=3365.81 Test=99.569 Tnew=-144277 F=9.35162 f=-1.35014e+06 Cp=9.35162Tnew = 99.569 K ，小于100，此时，下一步得到的 Tnew就完全错了。我认为就是因为 99.569K对于的cp函数返回值无定义，程序不知道给你返回什么。
所以，我认为你程序里必须给出 T < 100 和 T > 180 时 cp 的定义。
另外，还有一个问题是，理论上如果能量值在正常范围内，应该牛顿法迭代得到的温度值不会小于100k，我相信你也是这么认为的。但是，据我的理解，那个迭代的函数里面，参数 f 对应的应该是能量（焓或者内能）。从你输出的信息来看，你的 f 的值似乎一直都是负数，这应该是不对的。所以，可能你的算例设置也有点问题。
以上仅供参考，欢迎讨论。
-
"Tnew = 99.569 K ，小于100，此时，下一步得到的 Tnew就完全错了。我认为就是因为 99.569K对于的cp函数返回值无定义，程序不知道给你返回什么。
所以，我认为你程序里必须给出 T < 100 和 T > 180 时 cp 的定义。"
——————————————————————————————
这个我有点疑问，我的初始值在130K或者更高，壁面加热在170K左右，算出来的温度怎么会低于100K，而且当我限定了cp 在 T < 100 和 T > 180 时 的情况，比如T < 100令T=100，T > 180 令T=180，结果用不了几步还是发散。“还有一个问题是，理论上如果能量值在正常范围内，应该牛顿法迭代得到的温度值不会小于100k，我相信你也是这么认为的。但是，据我的理解，那个迭代的函数里面，参数 f 对应的应该是能量（焓或者内能）。从你输出的信息来看，你的 f 的值似乎一直都是负数，这应该是不对的。”
———————————————————————————————
这个跟hPolynomial本身焓的设定有关系，他有一个参考值（偏移量）：hCoeffs_[0] += Hf_ - hCoeffs_.value(Tstd_);template<int PolySize> Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const { scalar val = this->v_[0]; // avoid costly pow() in calculation scalar powX = 1; for (label i=1; i<PolySize; ++i) { powX *= x; val += this->v_[i]*powX; } if (logActive_) { val += logCoeff_*log(x); } return val; }理论上应该没问题，如果有问题那么在100～150K或150～180K也应当有问题。
-
有了点新发现，当在hPolynomialThermoI.H（src/thermophysicalModels/specie/thermo/hPolynomial/）的ha()函数中添加如下内容，能正常求解：
template<class EquationOfState, int PolySize> inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::ha ( const scalar p, const scalar T ) const { if ( T >= 100 && T < 150){ return hCoeffs_.value(T); } else { return hCoeffs_.value(T); } }但是如果改为下面这种就求解不了（当然Cp也是跟着改的）：
template<class EquationOfState, int PolySize> inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::ha ( const scalar p, const scalar T ) const { if ( T >= 100 && T < 150){ return hCoeffs_.value(T); } else { if ( T >= 150 && T < 180 ){ return hCoeffs1_.value(T);（实际上这里改为其他的比如多项式，查表之类的都一样） } else { return hCoeffs_.value(T); } } }template<class EquationOfState, int PolySize> inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp ( const scalar p, const scalar T ) const { // return CpCoeffs_.value(T); if ( T >= 100 && T < 150){ return CpCoeffs_.value(T); } else { if ( T >= 150 && T < 180 ){ return CpCoeffs1_.value(T); } else { return CpCoeffs_.value(T); } } }感觉问题是由这里引起的。
-
@李东岳 @xpqiu
我想我可能找到原因了：因为单段拟合时，Cp系数组成的多项式计算的Cp连续（Cp连续），进而由Cp系数算出的焓he的系数组成的多项式算出的焓he连续（焓值连续）。而当多段拟合时，由于各段Cp系数不一致导致算出的焓值he与其它多项式的焓值不连续，导致出现焓值跳跃点（焓值间断）。在焓值间断点算出的（h1-h0）异常而导致温度异常，最终发散。
对于这种情况是否还有其他出路？求教！:crying:CpCoeffs1_ *= this->W(); hCoeffs_ = CpCoeffs_.integral(); hCoeffs_[0] += Hf_ - hCoeffs_.value(Tstd_); return hCoeffs_.value(T);template<int PolySize> Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const { scalar val = this->v_[0]; // avoid costly pow() in calculation scalar powX = 1; for (label i=1; i<PolySize; ++i) { powX *= x; val += this->v_[i]*powX; } if (logActive_) { val += logCoeff_*log(x); } return val; }template<int PolySize> typename Foam::Polynomial<PolySize>::intPolyType Foam::Polynomial<PolySize>::integral(const scalar intConstant) const { intPolyType newCoeffs; newCoeffs[0] = intConstant; forAll(*this, i) { newCoeffs[i+1] = this->v_[i]/(i + 1); } return newCoeffs; }
