@xuqiming 没有理解错。
因为OpenFOAM里默认的transport model是sutherland公式,输运特性参数里的粘度的计算方法为:
$$
\mu=A_s \frac {T^{3/2}}{T+T_s}
$$
可以看出需要 $A_s$ 和 $T_s$ 这两个参数。这是一个非常简单的transport model。如果要用原生的任何燃烧求解器包括reactingFoam,那就是默认用sutherland公式,就只需要用OpenFOAM里 $A_s$ 和 $T_s$ 两个参数的transport文件。
关注分子输运特性的组分,chemkin形式的transport文件里的六列分子特性参数是用来计算该组分粘度和比热容的,最终得到混合物的特性的,这种就是mixture-averaged/multicomponent transport model。OpenFOAM原生代码里没有植入这俩模型。
看你用的是OpenFOAM-7,多说一句,OpenFOAM-8及以上的org版本增加了基于Fickian和MaxwellStefan模型求扩散系数的,这对求组分扩散非常重要,应该值得去关注一下,虽然这些模型用着也不是很方便。
wangfei9088
帖子
-
chemkinToFoam -
关于在求解器中使用热物理库中的函数@huangyuhui723 我明白你的意思了。
@huangyuhui723 在 关于在求解器中使用热物理库中的函数 中说:
通过其它计算更新了T之后想对he赋值
如果是求解TEqn,可以是得到温度T后更新he。OpenFOAM里的求解器一般都是求解EEqn,得到he后通过thermo.correct()更新T。
如果已知边界温度T后要对he赋值,应该在求解EEqn之前对he边界赋值,然后求解EEqn。大致代码是这样:
volScalarField& he = thermo.he(); //****************************************//加这一段 forAll(p.boundaryField(), patchi) { forAll(p.boundaryField()[patchi], facei) { he.boundaryFieldRef()[patchi][facei] = function(p.boundaryField()[patchi][facei], T.boundaryField()[patchi][facei]); } } //****************************************// fvScalarMatrix EEqn ( ...... );
这个function函数就是T与he的关系,类似上面,混合物的he通过各组分的显焓的质量分数加权得到,显焓用janaf里的关系确定。
供参考。
-
关于在求解器中使用热物理库中的函数@huangyuhui723 以of10为例。
对的,he是通过heThermo.C里这个函数构造的。
he_ ( IOobject ( BasicThermo::phasePropertyName ( MixtureType::thermoType::heName(), phaseName ), mesh.time().timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), volScalarFieldProperty //这个函数确定cell和patch的值 ( "he", dimEnergy/dimMass, &MixtureType::cellThermoMixture, //cell &MixtureType::patchFaceThermoMixture, //patch &MixtureType::thermoMixtureType::HE, this->p_, this->T_ ), this->heBoundaryTypes(), this->heBoundaryBaseTypes() ),
能看到he的边界是通过这个函数确定的:
template<class BasicThermo, class MixtureType> template<class PatchFaceMixture, class Method, class ... Args> Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::patchFieldProperty ( PatchFaceMixture patchFaceMixture, Method psiMethod, const label patchi, const Args& ... args ) const { tmp<scalarField> tPsi ( new scalarField(this->T_.boundaryField()[patchi].size()) ); scalarField& psi = tPsi.ref(); forAll(this->T_.boundaryField()[patchi], facei) { psi[facei] = ((this->*patchFaceMixture)(patchi, facei).*psiMethod) ( args[facei] ... ); //计算he边界上的值 } return tPsi; }
关键是这个patchFaceMixture函数,调用的是coefficientMultiComponentMixture.C里的:
template<class ThermoType> const typename Foam::coefficientMultiComponentMixture<ThermoType>::thermoMixtureType& Foam::coefficientMultiComponentMixture<ThermoType>::patchFaceThermoMixture ( const label patchi, const label facei ) const { mixture_ = this->Y()[0].boundaryField()[patchi][facei] *this->specieThermos()[0]; for (label i=1; i<this->Y().size(); i++) { mixture_ += this->Y()[i].boundaryField()[patchi][facei] *this->specieThermos()[i]; } return mixture_; }
可以看出混合物的he的边界值是通过上面这个函数更新的,即各个组分的质量分数加权求得的。如果需要更改边界值,修改这个函数。但是注意,直接修改这个函数会改变所有特性的值包括he。
speciesThermos()[i]指的是组分i的特性,OpenFOAM源代码中每个组分的焓值调用的是thermoI.H里的:
template<class Thermo, template<class> class Type> inline Foam::scalar Foam::species::thermo<Thermo, Type>::HE(const scalar p, const scalar T) const { return Type<thermo<Thermo, Type>>::HE(*this, p, T); }
调用的是sensibleEnthalpy.H里的:
scalar HE ( const Thermo& thermo, const scalar p, const scalar T ) const { return thermo.Hs(p, T); }
调用的是janafThermoI.H里的:
template<class EquationOfState> inline Foam::scalar Foam::janafThermo<EquationOfState>::Hs //显焓 ( const scalar p, const scalar T ) const { return Ha(p, T) - Hf(); } template<class EquationOfState> inline Foam::scalar Foam::janafThermo<EquationOfState>::Ha //总焓 ( const scalar p, const scalar T ) const { const coeffArray& a = coeffs(T); return ( ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T + a[5] ) + EquationOfState::H(p, T); } template<class EquationOfState> inline Foam::scalar Foam::janafThermo<EquationOfState>::Hf() const //生成焓 { const coeffArray& a = lowCpCoeffs_; return ( ( (((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd + a[0] )*Tstd + a[5] ); }
也就是说,如果要更改每个组分的特性,需要改变不同的janaf系数,这一般是不可能这么做的。当然,可以使用不同的thermo模型。
上面仅以coefficientMultiComponentMixture和janaf为例说明这些函数的调用关系。OpenFOAM0-10中很多调用关系,比如&MixtureType::thermoMixtureType::HE,理解会有点难。
如果觉得OpenFOAM0-10中这些关系难以理解,建议看看ESI版本的OpenFOAM(比如OpenFOAM-v2012)和基金会以前的版本(比如OpenFOAM-6和OpenFOAM-7),这些版本的代码直接很多,但意思是一样的,会容易理解很多。
供参考。
-
关于openfoam中计算气体密度的疑问@脉动通风 是的。
$$\rho = \frac {1} {RT} p= \psi p$$
$$ \psi = \frac {1} {RT}=\frac {1} {\frac{RR}{M} T}$$
在perfectGasI.H文件中可以看到:
$$ \psi = \frac {1} {RT}$$template<class Specie> inline Foam::scalar Foam::perfectGas<Specie>::psi(scalar p, scalar T) const { return 1.0/(this->R()*T); }
this->R()可以在speciesI.H文件里看到:(从这里可以看出来RR除以摩尔质量)
$$ R = \frac{RR}{M} $$inline scalar specie::R() const { return RR/molWeight_; }
这个RR在thermodynamicConstants.C文件中:
// Note: the 1e3 converts from /mol to /kmol for consistency with the // SI choice of kg rather than g for mass. // This is not appropriate for USCS and will be changed to an entry in // the DimensionedConstants dictionary in etc/controlDict const scalar RR = 1e3*physicoChemical::R.value();
physicoChemical::R.value()就是通用气体常数8.314。因为molWeight单位是kg/kmol,RR要乘以1000.
-
reactingFoam甲烷空气预混燃烧后处理@风云5091 在 reactingFoam甲烷空气预混燃烧后处理 中说:
怎么计算层流火焰的燃烧速度
laminar flame speed可以用这个公式[1]:
$$
S_{L} =- \frac {1} {\rho_1 Y_{F,1}} \int_{-\infty}^{\infty} \dot w_F dx
$$@风云5091 在 reactingFoam甲烷空气预混燃烧后处理 中说:
火焰厚度
flame thickness可以用这个公式[2]:
$$
\delta_{F} =\frac {T_{max} - T_{min}} {\rm max |\nabla T|}
$$@风云5091 在 reactingFoam甲烷空气预混燃烧后处理 中说:
reactingFoam似乎是一个不可压缩求解器,因此结果文件中并没有密度文件
reactingFoam是不可压缩求解器,applications/solvers/combustion/reactingFoam/createFields.H里密度是默认不输出的。要输出密度,需要将rho改为:
volScalarField rho ( IOobject ( "rho", runTime.timeName(), //mesh //comment mesh, //add IOobject::NO_READ, //add IOobject::AUTO_WRITE //add ), thermo.rho() );
@风云5091 在 reactingFoam甲烷空气预混燃烧后处理 中说:
质量分数分布转换成摩尔浓度
摩尔浓度与质量分数的关系:
$$
[X_k] = \rho \frac {Y_k} {W_k}
$$可以在createFields.H里定义一组摩尔浓度的体积标量场,size为化学组分的数量,设置成AUTO_WRITE。最简单直接的就是通过质量分数后处理得到。
@风云5091 在 reactingFoam甲烷空气预混燃烧后处理 中说:
红圈那块感觉不太合理?
红圈位置的速度不合理可能有两个原因:计算时间太短了,需要增大end time。如果没有改善,就是因为Courant number太大了,办法是调小maxCo。
供参考。
References:
[1] https://doi.org/10.1016/j.cpc.2018.11.011
[2] https://doi.org/10.1016/j.ijheatmasstransfer.2020.120127 -
请教各位老师或大佬一个简单的离散问题@李东岳 哦哦,懂了,那我都试试吧。谢谢李老师及时回复。 感恩的心
-
请教各位老师或大佬一个简单的离散问题请教各位老师或大佬一个简单的离散问题,求C的方程其中一项是与另一个场$\alpha$的梯度有关,D是常数,是这样的:
$$ \nabla \cdot (D \frac{\nabla \alpha}{\alpha} C)$$下面两行的写法在求解器中有什么区别吗?我该用哪个?或者都不对。感谢回复。
phiV = fvc::interpolate(D/α*fvc::grad(α)) & mesh.Sf(); phiV = fvc::interpolate(D/α)*fvc::snGrad(α); fvm::div(phiV, C, "div(phiV,C)")
-
rhoCentralFoam发散,出现Maximum number of iterations exceeded@gooseEast 在 rhoCentralFoam发散,出现Maximum number of iterations exceeded 中说:
在thermoI.H文件第46行标量f是什么
OpenFOAM代码中,这个f是Cp,定压比热容,单位是[J/(kg K)]。这个T函数是通过定压比热容的值迭代求出温度。limit,F,dFdt分别对应的是hConstThermoI.H文件里limit,Cp,dCpdt三个函数。
-
paraview 不显示坐标轴数值@liujm 假设你用的版本是OpenFOAM-10,其他版本也一样。
路径:OpenFOAM-10/etc/config.sh/
打开paraview文件,将69行取消注释,注释掉72行。
./makeParaView重新编译就是5.0.1版本了。#export ParaView_VERSION=5.0.1 //69行 #export ParaView_VERSION=5.4.0 #export ParaView_VERSION=5.5.0 export ParaView_VERSION=5.6.3 //72行
-
chemkin在网上找到的C12H26化学反应机理文件,不管是直接运行还是通过chemkinToFoam转换都显示类似格式错误的问题
因为OpenFOAM里chemisrtyReader和chemkinToFoam只会读transportProperties这个文件里的数据,类似下面这种。用在sutherland公式里。
".*" { transport { As 1.512e-06; Ts 120.; } } "H2" { transport { As 6.362e-07; Ts 72.; } } "CO2" { transport { As 1.572e-06; Ts 240.; } }
我在网上找到的trans文件都长类似的样子
chemkin软件的三个文件之一就是这个六列分子特性数据的trans.dat,这才是正确的。只是OpenFOAM因为用sutherland公式就不需要这个文件了。sutherland公式需要的是上面的As和Ts。
使用OpenFOAM自带的chemkinToFoam命令:chemkinToFoam chem.inp therm.dat transportProperties reactions.dat thermo.dat
-
porousGasificationFoam运行ODESolver报错@xuanze 哈哈,可以呀,我邮箱wangfei9088@hotmail.com。
-
porousGasificationFoam运行ODESolver报错谢谢东岳老师。
@xuanze 在 porousGasificationFoam运行ODESolver报错 中说:
确实需要降低时间步长,但是不知道怎样确定降低多少
时间步长由controlDict文件里的deltaT,maxDeltaT和maxCo这些共同决定,取最小值,一般调小maxCo即可,比如0.2试试。
@xuanze 在 porousGasificationFoam运行ODESolver报错 中说:
还有其他适合化学反应的ODE求解器吗
如果chemistryProperties文件里chemistryType的solver是ode的话,odeCoeffs里的sover就是目前case的ODESolver。看图片用的是seulex。低温燃烧和点火等尽量不要使用OpenFOAM里自带的seulex,因为可能有低温燃烧特性捕捉不准确,点火延迟时间预测不准,不稳定等这些问题。用它要先验证下。有文献用SIBS,可以试试。
-
运行结束后输出密度rho@myheart 是的。
-
运行结束后输出密度rho@myheart rhoReactingFoam也不能算可压。
reactingFoam和rhoReactingFoam分别用的是psiReactionThermo和rhoReactionThermo,区别是基于可压缩性和密度的热物性(thermodynamic properties),不是流场。 -
运行结束后输出密度rho@myheart reactingFoam不能算可压。
把createFields.H里rho的构造函数加两行就可以输出密度了,如下:
volScalarField rho ( IOobject ( "rho", runTime.timeName(), mesh, //加逗号 IOobject::NO_READ, // IOobject::AUTO_WRITE // ), thermo.rho() );
-
重新编译热物理库@ann 我不是大佬,这是大佬的代码,供参考。可能要根据你的OpenFOAM版本稍微改改。
https://github.com/ZhangYanTJU/ZYsprayFoam/blob/master/logSummary.H -
Wray-Agarwal湍流模型没人试试SA么,那我还是把坑填上吧。
@xpqiu 在 Wray-Agarwal湍流模型 中说:
A "trip" in the context of
@李东岳 在 Wray-Agarwal湍流模型 中说:
只需要搞ft2就可以。要简单不少。那个ft1太长了。
对,我也这么理解的。
含$f_{t2}$两项,不含$f_{t1}$的SA模型,of10版本。mySpalartAllmaras里wmake,pitzDaily里simpleFoam。感兴趣的可以试试,已测试,没问题,供参考。mySpalartAllmarasTurbulenceModel.tar.xz
注:11楼和14楼的源项少了一个系数$C_{b1}$,应该是:
$$
S_{\tilde \nu} = -\alpha \rho C_{b1} f_{t2} \tilde \nu \tilde S + \alpha \rho \frac{C_{b1}}{\kappa ^2} f_{t2} \left[ \frac{\tilde \nu}{y} \right]^2
$$ -
Wray-Agarwal湍流模型 -
基于reactingFoam的管道内预混氢气点火爆炸@尚善若水 在 基于reactingFoam的管道内预混氢气点火爆炸 中说:
应该都没有特别考虑这个功能,这个也许可以 https://github.com/JieSun-pku/detonationFoam
对,这个应该可以。
-
Wray-Agarwal湍流模型@李东岳 哈哈,好。
根据文献和网页,写一下SA湍流模型的公式,方便在OpenFOAM中植入。
$$
\frac{\partial \tilde \nu}{\partial t}+\nabla \cdot (\mathbf{U} \tilde \nu) = c_{b1}(1-f_{t2}) \tilde S \tilde \nu +\frac{1}{\sigma} \left[ \nabla \cdot ((\nu + \tilde \nu)\nabla \tilde \nu) +c_{b2} (\nabla \tilde \nu)^2) \right] - \left[ c_{w1} f_w - \frac{c_{b1}}{\kappa ^2} f_{t2} \right] \left[ \frac{\tilde \nu}{d} \right]^2 + f_{t1} \Delta \mathbf{U} ^2
$$$$
\nu_t = \tilde \nu f_{v1} ,
f_{v1} = \frac{\chi ^3}{\chi ^3 + c_{v1} ^3} ,
\chi = \frac{\nu_t} {\nu}
$$$$
\tilde S = \Omega + \frac {\tilde \nu} {\kappa ^2 d^2} f_{v2} ,
\Omega = \sqrt {2} |\mathbf{W}|,
\mathbf{W} = \frac{1}{2} (\nabla \mathbf{U} - \nabla \mathbf{U} ^T) ,
f_{v2} = 1-\frac{\chi}{1+\chi f_{v1}}
$$
说明:原文用的是$S$表示涡量,这里和网页保持一致,用的$\Omega$。$$
f_w = g \left[\frac {1+c_{w3}^6}{g^6+c_{w3}^6} \right]^{1/6}
$$$$
g = r + c_{w2} (r^6 - r) ,
r = min \left[\frac{\tilde \nu}{\tilde S \kappa ^2 d^2}, 10 \right]
$$$$
f_{t2} = c_{t3} exp(- c_{t4} \chi ^2)
$$$$
f_{t1} = c_{t1} g_t exp \left(- c_{t2} \frac{\omega_t ^2}{\Delta \mathbf{U} ^2}[d^2+g_t ^2 d_t^2] \right) ,
g_t = min \left(0.1, \frac{\Delta \mathbf{U}}{\omega_t \Delta x_t} \right)
$$$$
c_{b1} = 0.1355 ,
\sigma = 2/3 ,
c_{b2} = 0.622 ,
\kappa = 0.41 ,
c_{w1} = \frac{c_{b1}}{\kappa ^2} + \frac{1+c_{b2}}{\sigma}
$$$$
c_{w2} = 0.3 ,
c_{w3} = 2 ,
c_{v1} = 7.1 ,
c_{t1} = 1 ,
c_{t2} = 2 ,
c_{t3} = 1.2 ,
c_{t4} = 0.5
$$壁面边界条件:
$$
\tilde \nu_{wall} = 0 ,
\nu_{t,wall} = 0
$$far field边界条件:
注:方程有无最后一项,far field边界条件是不一样的,如下:
$$
0 \leq \tilde \nu_{farfield} < \frac{1}{10} \nu_{\infty} ,
0 \leq \nu_{t,farfield} < 0.27940 \times 10^{-6} \nu_{\infty}
$$$$
3\nu_{\infty} \leq \tilde \nu_{farfield} \leq 5 \nu_{\infty} ,
0.210438\nu_{\infty} \leq \nu_{t,farfield} \leq 1.294234 \nu_{\infty}
$$