这个VectorSpace是OpenFOAM中的primitive type,所有的primitive type都需要能够使用C++的萃取(traits)技术提炼基本信息。
话说OpenFOAM把C++的高级特性都用上了,这架构设计的水平实在是高!
这个VectorSpace是OpenFOAM中的primitive type,所有的primitive type都需要能够使用C++的萃取(traits)技术提炼基本信息。
话说OpenFOAM把C++的高级特性都用上了,这架构设计的水平实在是高!
总结的不错。
Su Sp 是灵活处理源项的方法,正确使用方法如下:
1、若源项为负,如你举的例子中的 -30
,一般写作 $\frac{-30}{k^{old}}k$,$k^{old}$ 代表上一时刻的值,采用 fvm::Sp
离散后,源项就进入 [A] 中对角元素,使对角元素更占优。
2、若源项为正,如你举的例子中的 50
,则直接用 fvm::Su
离散,源项进入右端项 [b] 中。
3、若源项是一个有正有负的 volScalarField,则利用 fvm::SuSp
离散,它会判断源项每个网格单元的正负号,然后选择用 Su 还是 Sp 的方法进行处理。
以上方法的最终目的都是增加矩阵对角占优,使线性方程组的迭代更稳定,更容易收敛。
OUTLET的value这行后面少了个分号。
没安装make
如果在声明变量的时候指定名称:
volScalarField rAU(1.0/UEqn.A(), "rAU");
那么可以在字典文件中用以下语句指定离散格式
laplacian(rAU,p) Gauss linear corrected
在声明rAU
时不指定名称:
volScalarField rAU(1.0/UEqn.A());
那么rAU
的名称通过量纲以及符号运算得到,1.0/UEqn.A()
的名称为(1|A(U))
,因此相应的离散格式应该写作:
laplacian((1|A(U)),p) Gauss linear corrected
思路没问题,建议你用OpenFOAM的mpirunDebug来调试并行程序。
@mohui 感谢分享!
试了一下cavity算例,可以正常跑完。
将system/fvSolution里面的压力求解器的relTol设为0就可以使残差降下来,不过相应的矩阵求解迭代步要增加。
编译PVReader需要paraview的头文件,也就是说你的paraview需要是自己编译的才行。
如果不想自己编译paraview,paraFoam -builtin
先凑合用吧。
更正一下问题:一个类的成员函数怎样在另一个类里使用?
可以通过对象实例+成员访问运算符.
或->
调用。
如:compressibleTwoPhaseMixture
类中定义了一个变量twoPhaseMixtureThermo& mixture
,那么可以通过mixture.rho1()
及mixture.rho2()
访问mixture
这个对象的相应成员函数。
调解后的时间步进的长度delta_t都要小于设置的delta_T
调节后的时间步大小取决于controlDict
中的maxCo
值,maxCo
的默认值为1。代码见src/finiteVolume/cfdTools/general/include/createTimeControls.H
以及
src/finiteVolume/cfdTools/general/include/readTimeControls.H
。
在每个delta_t中都要跑一遍nOuterCorrectors(50)和nCorrectors (2)对不?
没错。
干吗不用piso算法自己手动调节时间步长
可以通过增大maxCo
来增大自动调节的时间步长,据说PIMPLE算法的Co数可以达到10以上,不过对于用DES或LES计算的湍流问题另当别论。
如何给定nOuterCorrctors 的值呢?
这个只能靠经验了,不过一般肯定用不着50。
这句话怎么理解?什么叫high turbulent flow??
当simulationType
设置为RAS或LES时,求解的是考虑湍流粘度($\nu_t$)影响的RANS/LES方程,而湍流粘度取决于求解湍流输运方程后得到的各湍流系数。
turbOnFinalIterOnly
为false代表每个PIMPLE循环内求解多次湍流输运方程,每次PIMPLE循环求解RANS/LES方程时的湍流系数都用的是上个PIMPLE循环内得到的湍流系数。
turbOnFinalIterOnly
为true则代表每个PIMPLE循环内只求解一次湍流输运方程,每次PIMPLE循环求解RANS/LES方程时的湍流系数都用的是上个时间步得到的湍流系数。
高雷诺数流动为了精确得到速度压力场,需要每个PIMPLE循环内都求解多次湍流输运方程。
fvm和fvc是OpenFOAM中的两个命名空间,fvm中的函数(或称操作符)将场量离散,返回的是fvMatrix,而fvc中的函数则是显式调用,返回仍然是场量。
这个编译出现的warning是由于当前系统时间与NFS分区服务器时间不一致所导致的,不影响使用。解决这个warning的办法就是将两台机器的时间同步。
实际上用RANS来求解非定常问题均可以称之为URANS。
@supersoldier 通过求解器里面的runTime.write()
触发,runTime.write()
将调用Foam::Time::write()
向磁盘写入所有注册对象。
@cfd-china label
和 WM_LABEL_SIZE
的设置有关。WM_LABEL_SIZE=32
时,label
为 int32
;WM_LABEL_SIZE=64
时,label
为 int64
。
@CFDngu 回答一下你的第二和第三个问题。
第二个问题:建议使用 fvPatchField<Type>::operator==
,而不是 Field<Type>::operator=
。 fvPatchField<Type>::operator==
这个操作符被重载,有多个实现。
557 template<class Type>
558 void Foam::fvPatchField<Type>::operator==
559 (
560 const fvPatchField<Type>& ptf
561 )
562 {
563 Field<Type>::operator=(ptf);
564 }
565
566
567 template<class Type>
568 void Foam::fvPatchField<Type>::operator==
569 (
570 const Field<Type>& tf
571 )
572 {
573 Field<Type>::operator=(tf);
574 }
575
576
577 template<class Type>
578 void Foam::fvPatchField<Type>::operator==
579 (
580 const Type& t
581 )
582 {
583 Field<Type>::operator=(t);
584 }
第三个实现的参数类型是const Type&
,赋值后整个patchField都是相同的值。
第二个实现的参数类型是const Field<Type>& tf
,赋值后整个patchField的值可以不同。
第三个问题:不指定作用域调用的operator=(Uabsfinal);
默认是this->operator=(Uabsfinal);
@LXJCFD 审稿人是不是以为你用的是单相自由液面求解方法?
@cfd-china 离散的具体实现可参阅相关源码。
以fvm::ddt(T)
为例,调用的是
template<class Type>
tmp<fvMatrix<Type> >
ddt
(
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return fv::ddtScheme<Type>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
)().fvmDdt(vf);
}
这里调用了fv::ddtScheme<Type>::New
,这个函数是一个runtime construction的实现,即在运行时读取字典文件中的离散格式并构造相应的对象:
template<class Type>
tmp<ddtScheme<Type> > ddtScheme<Type>::New
(
const fvMesh& mesh,
Istream& schemeData
)
{
if (fv::debug)
{
Info<< "ddtScheme<Type>::New(const fvMesh&, Istream&) : "
"constructing ddtScheme<Type>"
<< endl;
}
if (schemeData.eof())
{
FatalIOErrorIn
(
"ddtScheme<Type>::New(const fvMesh&, Istream&)",
schemeData
) << "Ddt scheme not specified" << endl << endl
<< "Valid ddt schemes are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
const word schemeName(schemeData);
typename IstreamConstructorTable::iterator cstrIter =
IstreamConstructorTablePtr_->find(schemeName);
if (cstrIter == IstreamConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"ddtScheme<Type>::New(const fvMesh&, Istream&)",
schemeData
) << "Unknown ddt scheme " << schemeName << nl << nl
<< "Valid ddt schemes are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
return cstrIter()(mesh, schemeData);
}
以常见的Euler
格式为例,将会构造EulerDdtScheme
类型的对象,并用其fvmDdt()
方法返回的fvMatrix,看fvmDdt()的实现
:
template<class Type>
tmp<fvMatrix<Type> >
EulerDdtScheme<Type>::fvmDdt
(
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type> > tfvm
(
new fvMatrix<Type>
(
vf,
vf.dimensions()*dimVol/dimTime
)
);
fvMatrix<Type>& fvm = tfvm();
scalar rDeltaT = 1.0/mesh().time().deltaTValue();
fvm.diag() = rDeltaT*mesh().Vsc();
if (mesh().moving())
{
fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
}
else
{
fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
}
return tfvm;
}
可以看到fvm::ddt(T)
产生的fvm.diag()
为时间倒数乘以网格单元体积,其产生的源项fvm.source()
为时间倒数乘以上一时刻的T乘以网格单元体积。
对于fvm::div及fvm::laplacian可以做类似分析。
插一句,fvm::div离散不会产生source(),fvm::laplacian会产生source()。
@金石为开 想要显示自由面高度首先需要得到自由面形状,自由面形状通过做alpha.water=0.5的contour得到。得到自由面形状之后再用高度着色即可。
参考 $FOAM_ETC/caseDicts/postProcessing/flowRate/flowRatePatch
。
OpenFOAM-4.0 可以直接在controlDict中使用 #includeFunc
语句 [1],用法参考user-guide [2]。
[1] https://github.com/OpenFOAM/OpenFOAM-dev/commit/dd20d9086889ccbf993e8648d8e8f070f1b16b79
[2] http://cfd.direct/openfoam/user-guide/post-processing-cli/
@子仲无未 OpenFOAM 没有造波功能,waves2Foam是第三方实现,如果要在OpenFOAM的solver中使用waves2Foam则需要修改源码重新编译。
@李东岳 我认为截止尺度和widthCoeff不是一个东西。
试一试在 decomposeParDict 里面使用 preservePatches [1] 。
@史浩 可以试试mpirun的--bind-to-core参数
首先,我不认为coupled solver可以允许更大的时间步长。
其次,从理论上说SIMPLE算法不适用于瞬态问题,Fluent里面使用了Transient SIMPLE算法来计算瞬态问题,时间步长取很大也可以算下去,但是很显然时间项的精度会很低,算出来的东西不可信。
你这个问题很有意思,从你的推导分析来看,SuSp 确实是削弱了对角占优。
但其中有一处错误:V30 应该为 V-30。不过这改变的是源项,而与对角元素无关。我们一起分析下是怎么回事。
如果将方程按照这种形式写:$\frac{\partial k}{\partial t}=-(-A)-(-B)$
那么代码可以写成:== -fvm::SuSp(-A/k, k) - fvm::SuSp(-B/k, k)
。
实际计算时,假设 A=-30,B=50。
那么 $-A/k>0$,$-B/k<0$,fvm前面的负号与 ==
抵消,代码实际执行的是
fvm.diag() += V*30
fvm.source() -= V*(-50)*k
这又的确增加了对角占优,显然 fvm::SuSp
不能随心所欲地用!
-fvm::SuSp(-A/k, k)
和 fvm::SuSp(A/k, k)
的结果不一样,而前一种写法才是我们想要的结果。
我翻了 OpenFOAM 的代码后发现,凡是在 ==
右边使用的 Sp 和 SuSp 前面必有负号,凡是在 ==
左边使用的 Sp 前面必为正号 (代码中基本没有出现 fvm::Su
)。
看来 SuSp 是个深坑,用起来要非常小心才行。
因为源项通常写在控制方程等号的右边,为了保证源代码与控制方程形式统一, fvm::Su
fvm::Sp
以及 fvm::SuSp
必须写在 ==
右边。
你可以看一下 OpenFOAM 中湍流模型中的用法。
k_()
调用的是 GeometricField
的 operator()
,定义如下:
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
//- Return a const-reference to the dimensioned internal field
// Useful in the formulation of source-terms for FV equations
inline const Internal& operator()() const;
@random_ran 我猜可能是 int 类型长度限制的原因。你看一下 echo $WM_LABEL_SIZE
,默认是 32,如果要用超大网格,需要在编译前将 WM_LABEL_SIZE
设为 64。
@东岳 推荐李博一个小软件:mathpix,可以将图片转化成 LaTeX 代码,准确率很高。
@zwl 这些都是类模板,你需要找到对这些类模板进行特化的源文件。
参考:
thermophysicalModels/basic/psiThermo/psiThermos.C
thermophysicalModels/basic/rhoThermo/rhoThermos.C
@一二 可以说出来给大家分享一下哦。
可以使用 $ $ 来插入latex代码,例如:$\alpha$;
也可以使用:
\begin{equation}
\alpha + \nabla \cdot p = 0
\end{equation}
来插入居中带索引号的公式。目前公式需要刷新一下才能显示。
@simthere 你的问题其实可以简化为:为什么 $\int \rho dV = \rho_P \int dV$。(其中 $\rho_P$ 为格心处的流体密度)
下面给出简略证明:
我们假设P为单元格心,$\textbf x_P$ 为P点坐标,那么有如下定义
$\int (\textbf x - \textbf x_P) dV = 0$
于是密度对控制单元的积分为:
$\int \rho dV = \int [\rho_P + (\textbf x - \textbf x_P) \cdot (\nabla \rho_P)] dV = \rho_P \int dV + (\nabla \rho_P) \cdot \int (\textbf x - \textbf x_P) dV = \rho_P \int dV$
同理可以推广到压力、速度等物理量。
x方向流入质量:$(\rho u)|_{x} \Delta y \Delta z$
x方向流出质量:$(\rho u)|_{x+\Delta x} \Delta y \Delta z$
y方向流入质量:$(\rho v)|_{y} \Delta x \Delta z$
y方向流出质量:$(\rho v)|_{y+\Delta y} \Delta x \Delta z$
z方向流入质量:$(\rho w)|_{z} \Delta x \Delta y$
z方向流出质量:$(\rho w)|_{z+\Delta z} \Delta x \Delta y$
于是,
$\frac{\partial \rho}{\partial t}\Delta x \Delta y \Delta z = ((\rho u)|_{x} - (\rho u)|_{x+\Delta x}) \Delta y \Delta z + ((\rho v)|_{y} - (\rho v)|_{y+\Delta y}) \Delta x \Delta z + ((\rho w)|_{z} - (\rho w)|_{z+\Delta z}) \Delta x \Delta y$
两边除以 $\Delta x \Delta y \Delta z$
$\frac{\partial \rho}{\partial t} = \frac{(\rho u)|_{x} - (\rho u)|_{x+\Delta x}}{\Delta x} + \frac{(\rho v)|_{y} - (\rho v)|_{y+\Delta y}} {\Delta y} + \frac{(\rho w)|_{z} - (\rho w)|_{z+\Delta z}} {\Delta z}$
所以
$\frac{\partial \rho}{\partial t} = - \left ( \frac{\partial}{\partial x} (\rho u) + \frac{\partial}{\partial y} (\rho v) + \frac{\partial}{\partial z} (\rho w) \right ) = - \nabla \cdot \rho \textbf U$
对控制体积分,运用高斯定理,可得连续性方程的积分形式
$\frac{\partial}{\partial t} \int \rho dV = - \oint \rho \textbf U \cdot \textbf n dS$
我的理解是
See "15.5.5. Advection Scheme for Turbulence Equations"
Irrespective of the Advection Scheme setting, the advection scheme for the turbulence model equations depends on the Turbulence Numerics option: ...
Note that Turbulence Numerics option also specifies the transient scheme for the turbulence model equations (see Turbulence Numerics in the CFX-Pre User's Guide). This can be overridden in CFX-Pre in the Solver Control details view on the Equation Class Settings tab.
@东岳 硕博论文可以到中国知网的海外版下载pdf,网址:http://gb.oversea.cnki.net/