首先谢谢李老师的回复，我的计算模型长这个样子
最常见的报错信息是这个样子：
--> FOAM FATAL ERROR: (openfoam-2306)
Maximum number of iterations exceeded: 100 when starting from T0:300 old T:18296.6 new T:3283.45 f:1.71386e+07 p:100000 tol:0.03
简言之就是发散了。
上面这个是二维网格，计算结果长这个样子（都给把计算结果对称了）：速度云图
温度云图：
感觉计算结果就是错的，因为在Cathode是wall类型，从Inlet附近到Cathode附近流线应该有弯曲的表现。
阴极我设定的边界条件3500K，它最高温度才3800K，磁场产生的源项作用效果十分小。
谢谢李老师提的reconstruct方法，我研究试试。
另外还想请教大家，在OpenFOAM中可以引入柱坐标系吗，轴向，径向，切向。就像Fluent中计算Axisymmetric Swirl模型，将2D拓展成旋转3D计算域。
同样求推荐资料或者书籍，比如介绍OpenFOAM数据类型，场的类型，遍历边界，遍历内部场，将数学公式转化为OpenFOAM代码，等等。
大家早上好，我是想通过OpenFOAM研究磁流体的同学，我自己做的求解器（下面详说）收敛性极差，而且感觉算得不太对（本身代码就感觉写得很low）。主要想请教大家看看我将数学方程转化为代码的形式哪里有问题，给我些提高指导意见。
首先需要添加电磁场方程
$\nabla \cdot (\sigma \nabla \phi) = 0$（标量运输方程）
$\nabla \cdot (\nabla \mathbf{A}) = -\mu_0 \mathbf{j}$ （矢量运输方程）
其中$\phi$和$\mathbf{A}$是新添加的场电势和磁矢量势，$\sigma$是电导率，$\mu_0$是一个常数，$\mathbf{j} = - \sigma\nabla \phi$是电流密度。
需要在动量方程添加一个源项洛伦兹力$\mathbf{F}_{Lorentz} = \mathbf{j} \times \mathbf{B}$, 其中$\mathbf{B} = \nabla \times \mathbf{A}$。
需要在能量方程添加源项焦耳热$S=\sigma \mathbf{E}^2$和电子运输焓$\frac{5}{2}\frac{k_B}{e}\mathbf{j}\cdot \nabla T$，其中$\mathbf{E} = -\nabla \phi$，$k_B$是玻尔兹曼常数，$e$是元电荷量（常数），$T$是温度。
我以buoyantSimpleFoam为基础，添加了如下代码来求解磁场方程：
solve(fvm::laplacian(sigma, PotE));
solve
(
fvm::laplacian(A) == Foam::constant::electromagnetic::mu0*sigma*fvc::grad(PotE)
);
动量方程在buoyantSimpleFoam的UEqn.H上只添加了洛伦兹力源项代码如下：
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
//添加动量方程源项洛伦兹力
- sigma*(-fvc::grad(PotE) ^ fvc::curl(A))
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvOptions.constrain(UEqn);
if (simple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvOptions.correct(U);
}
能量方程就是在原有的方程基础上添加了两个源项代码如下：
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::div(phi, he)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
)
- fvm::laplacian(turbulence->alphaEff(), he)
// 添加能量方程源项焦耳热
- (fvc::grad(PotE) & fvc::grad(PotE))*sigma
// 添加能量方程源项电子运输焓
- (
(2.50*Foam::constant::physicoChemical::k / Foam::constant::electromagnetic::e)
* sigma
* (-fvc::grad(PotE) & fvc::grad(he))
/ thermo.Cp()
)
==
rho*(U&g)
+ radiation->Sh(thermo, he)
+ fvOptions(rho, he)
);
EEqn.relax();
fvOptions.constrain(EEqn);
EEqn.solve();
fvOptions.correct(he);
thermo.correct();
radiation->correct();
}
其中$\nabla T$用$\frac{\nabla he}{thermo.Cp()}$表示了。
谢谢大家的批评指正。
想请教大家标量的梯度不应该是个矢量吗，为什么设置边界条件的时候fixedGradient是个标量？
我在使用buoyantSimpleFoam的时候也出现了如下报错：
--> FOAM FATAL ERROR: (openfoam-2306) Maximum number of iterations exceeded: 100 when starting from T0:300 old T:20112.8 new T:-884.425 f:1.0349e+07 p:105000 tol:0.03
这段报错代码位于$FOAM_SRC/thermophysicalModels/specie/thermo/thermo/thermoI.H: 67~86行的do循环里，这段应该是热物理模型计算温度的代码。
这里有个项目（第24页）绕过了循环求解温度，直接通过temperature-enthalpy表格插值来计算温度：https://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2017/DanielHummel/hummelReport.pdf
想请教各位大佬在thermoI.H文件第46行标量f是什么，项目里说是current
enthalpy，是比焓吗？单位是 J/kg 吗？
上面的报错信息在尤其在设定为第二类边界条件的时候基本不能运行起来，所以我打算仿照上面这个项目的办法绕过原来的循环求解温度，通过比焓（J/kg）-温度（K）来实现温度场的计算。
附一个thermoI.H文件：
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "thermo.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Thermo, template<class> class Type>
inline Foam::species::thermo<Thermo, Type>::thermo
(
const Thermo& sp
)
:
Thermo(sp)
{}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::T
(
scalar f,
scalar p,
scalar T0,
scalar (thermo<Thermo, Type>::*F)(const scalar, const scalar) const,
scalar (thermo<Thermo, Type>::*dFdT)(const scalar, const scalar)
const,
scalar (thermo<Thermo, Type>::*limit)(const scalar) const
) const
{
if (T0 < 0)
{
FatalErrorInFunction
<< "Negative initial temperature T0: " << T0
<< abort(FatalError);
}
scalar Test = T0;
scalar Tnew = T0;
scalar Ttol = T0*tol_;
int iter = 0;
do
{
Test = Tnew;
Tnew =
(this->*limit)
(Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));
if (iter++ > maxIter_)
{
FatalErrorInFunction
<< "Maximum number of iterations exceeded: " << maxIter_
<< " when starting from T0:" << T0
<< " old T:" << Test << " new T:" << Tnew
<< " f:" << f
<< " p:" << p
<< " tol:" << Ttol
<< abort(FatalError);
}
} while (mag(Tnew - Test) > Ttol);
return Tnew;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Thermo, template<class> class Type>
inline Foam::species::thermo<Thermo, Type>::thermo
(
const word& name,
const thermo& st
)
:
Thermo(name, st)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Thermo, template<class> class Type>
inline Foam::word
Foam::species::thermo<Thermo, Type>::heName()
{
return Type<thermo<Thermo, Type>>::energyName();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Cpv(const scalar p, const scalar T) const
{
return Type<thermo<Thermo, Type>>::Cpv(*this, p, T);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::gamma(const scalar p, const scalar T) const
{
#ifdef __clang__
volatile const scalar Cp = this->Cp(p, T);
#else
const scalar Cp = this->Cp(p, T);
#endif
return Cp/(Cp - this->CpMCv(p, T));
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::CpByCpv
(
const scalar p,
const scalar T
) const
{
return Type<thermo<Thermo, Type>>::CpByCpv(*this, p, T);
}
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);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::G(const scalar p, const scalar T) const
{
return this->Ha(p, T) - T*this->S(p, T);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::A(const scalar p, const scalar T) const
{
return this->Ea(p, T) - T*this->S(p, T);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::cp(const scalar p, const scalar T) const
{
return this->Cp(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::ha(const scalar p, const scalar T) const
{
return this->Ha(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::hs(const scalar p, const scalar T) const
{
return this->Hs(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::hc() const
{
return this->Hc()*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::s(const scalar p, const scalar T) const
{
return this->S(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::he(const scalar p, const scalar T) const
{
return this->HE(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::cv(const scalar p, const scalar T) const
{
return this->Cv(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::es(const scalar p, const scalar T) const
{
return this->Es(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::ea(const scalar p, const scalar T) const
{
return this->Ea(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::g(const scalar p, const scalar T) const
{
return this->G(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::a(const scalar p, const scalar T) const
{
return this->A(p, T)*this->W();
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::K(const scalar p, const scalar T) const
{
scalar arg = -this->Y()*this->G(Pstd, T)/(RR*T);
if (arg < 600)
{
return exp(arg);
}
else
{
return VGREAT;
}
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kp(const scalar p, const scalar T) const
{
return K(p, T);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kc(const scalar p, const scalar T) const
{
const scalar nm = this->Y()/this->W();
if (equal(nm, SMALL))
{
return Kp(p, T);
}
else
{
return Kp(p, T)*pow(Pstd/(RR*T), nm);
}
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::Kx
(
const scalar p,
const scalar T
) const
{
const scalar nm = this->Y()/this->W();
if (equal(nm, SMALL))
{
return Kp(p, T);
}
else
{
return Kp(p, T)*pow(Pstd/p, nm);
}
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::Kn
(
const scalar p,
const scalar T,
const scalar n
) const
{
const scalar nm = this->Y()/this->W();
if (equal(nm, SMALL))
{
return Kp(p, T);
}
else
{
return Kp(p, T)*pow(n*Pstd/p, nm);
}
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THE
(
const scalar he,
const scalar p,
const scalar T0
) const
{
return Type<thermo<Thermo, Type>>::THE(*this, he, p, T0);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THs
(
const scalar hs,
const scalar p,
const scalar T0
) const
{
return T
(
hs,
p,
T0,
&thermo<Thermo, Type>::Hs,
&thermo<Thermo, Type>::Cp,
&thermo<Thermo, Type>::limit
);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THa
(
const scalar ha,
const scalar p,
const scalar T0
) const
{
return T
(
ha,
p,
T0,
&thermo<Thermo, Type>::Ha,
&thermo<Thermo, Type>::Cp,
&thermo<Thermo, Type>::limit
);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::TEs
(
const scalar es,
const scalar p,
const scalar T0
) const
{
return T
(
es,
p,
T0,
&thermo<Thermo, Type>::Es,
&thermo<Thermo, Type>::Cv,
&thermo<Thermo, Type>::limit
);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::TEa
(
const scalar ea,
const scalar p,
const scalar T0
) const
{
return T
(
ea,
p,
T0,
&thermo<Thermo, Type>::Ea,
&thermo<Thermo, Type>::Cv,
&thermo<Thermo, Type>::limit
);
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::dKcdTbyKc
(
const scalar p,
const scalar T
) const
{
const scalar dKcdTbyKc =
(this->S(Pstd, T) + this->Gstd(T)/T)*this->Y()/(RR*T);
const scalar nm = this->Y()/this->W();
if (equal(nm, SMALL))
{
return dKcdTbyKc;
}
else
{
return dKcdTbyKc - nm/T;
}
}
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::dcpdT(const scalar p, const scalar T) const
{
return this->dCpdT(p, T)*this->W();
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Thermo, template<class> class Type>
inline void Foam::species::thermo<Thermo, Type>::operator+=
(
const thermo<Thermo, Type>& st
)
{
Thermo::operator+=(st);
}
template<class Thermo, template<class> class Type>
inline void Foam::species::thermo<Thermo, Type>::operator*=(const scalar s)
{
Thermo::operator*=(s);
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Thermo, template<class> class Type>
inline Foam::species::thermo<Thermo, Type> Foam::species::operator+
(
const thermo<Thermo, Type>& st1,
const thermo<Thermo, Type>& st2
)
{
return thermo<Thermo, Type>
(
static_cast<const Thermo&>(st1) + static_cast<const Thermo&>(st2)
);
}
template<class Thermo, template<class> class Type>
inline Foam::species::thermo<Thermo, Type> Foam::species::operator*
(
const scalar s,
const thermo<Thermo, Type>& st
)
{
return thermo<Thermo, Type>
(
s*static_cast<const Thermo&>(st)
);
}
template<class Thermo, template<class> class Type>
inline Foam::species::thermo<Thermo, Type> Foam::species::operator==
(
const thermo<Thermo, Type>& st1,
const thermo<Thermo, Type>& st2
)
{
return thermo<Thermo, Type>
(
static_cast<const Thermo&>(st1) == static_cast<const Thermo&>(st2)
);
}
// ************************************************************************* //
我在使用buoyantSimpleFoam时需要将流体密度、黏度、热导率、定压比热容等设置成温度的插值函数，通过如下的thermophysicalProperties字典文件能够实现依照温度插值的效果吗？
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type heRhoThermo;
mixture pureMixture;
transport tabulated;
thermo hTabulated;
equationOfState icoTabulated;
specie specie;
energy sensibleEnthalpy;
}
mixture
{
specie
{
molWeight 1;
}
equationOfState
{
rho
(
(300 1)
(400 1)
...
(1000 1)
);
}
thermodynamics
{
Hf 0;
Sf 0;
Cp
(
(300 1)
(400 1)
...
(1000 1)
);
}
transport
{
mu
(
(300 1)
(400 1)
...
(1000 1)
);
kappa
(
(300 1)
(400 1)
...
(1000 1)
);
}
}
// ************************************************************************* //
最近也在折腾做旋转对称圆柱型计算域，找到个这：https://www.fangzhenxiu.com/post/1826059?code=031TTpFa1IW6xG0skbHa1OsJy23TTpFt&state=afterWeixin
根据这个操作画出来的网格虽然是楔形块，但是楔面并不参照哪一坐标平面对称，依然不能用在OpenFOAM上。
又找到一个这个：https://www.fangzhenxiu.com/post/1829359
作者好像通过两次旋转实现了对称的楔面，但我没有复现出来。
最后还是打算用blockMesh来画了。
删除掉0文件下CATHODE的gradientOld行代码，即可在paraview里查看0时刻场信息。
上面李老师的方法对于小白同学依然门槛太高，下边是我通过swak4Foam（OpenFOAM-v2306环境下）中的funkySetBoundaryField来实现的：
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object funkySetBoundaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
currentDensity
{
field PotE;
expressions
(
{
target gradient;
patchName CATHODE;
variables
(
"Jcath0=1;"
"Rc=1;"
"nc=1;"
"sigma0=1;"
"y=pos().y;"
"r=1;"
);
expression "-Jcath0*exp(-pow(mag(y) / Rc, nc)) / sigma0";
}
);
}
// ************************************************************************* //
想要实现速度边界条件为高斯分布的同学可以参考代码修改field和target，想要了解其他信息，可以参见：https://www.cfd-online.com/Forums/openfoam-pre-processing/96280-funkysetboundaryfields-manipulation-existing-field.html
哦哦好的好的，谢谢李老师
是的，李老师你讲的没错，J_cath = C1 * exp(- (r / C2)^C3)，假设是一个垂直X轴的面，C1、C2、C3均为已知参数，r是距X轴的距离。这样就是一个非均匀的fixedGradient。该怎么实现呢？
首先恭喜一下我终于顺利成为该论坛的用户！
我创建了一个新场电势，其中阴极边界条件如图：
Jcath是一个仅与边界相关的位置函数，sigma是电导率。
类似于codeFixedValue，想请教怎么对gradient边界条件编程。
如有案例，不胜感激。
另外有没有小伙伴做磁流体或者等离子体仿真，我们可以详细交流。