关于 phi 和 fvVectorMatrix 的两个问题
-
刚开始学习openfoam，遇到两个问题， 前来请教。
我在cavity的算例里面建立了4*4的网格，如下
vertices ( (0 0 0) (1 0 0) (1 1 0) (0 1 0) (0 0 0.1) (1 0 0.1) (1 1 0.1) (0 1 0.1) ); blocks ( hex (0 1 2 3 4 5 6 7) (4 4 1) simpleGrading (1 1 1) );然后输出了flux phi的值
Info<< " phi " << phi << nl <<endl; Info<< " phi size " << phi.size() << nl <<endl;得到如下结果：
phi dimensions [0 3 -1 0 0 0 0]; internalField nonuniform List<scalar> 24 ( -2.28233e-06 2.28233e-06 -3.2402e-06 9.57866e-07 -2.28434e-06 -9.55866e-07 -2.28433e-06 -3.94416e-06 6.22649e-06 -5.34147e-06 2.35517e-06 -3.95016e-06 -2.34716e-06 -6.2345e-06 -6.72257e-06 1.29491e-05 -7.99247e-06 3.62507e-06 -6.71657e-06 -3.62306e-06 -1.29511e-05 1.29491e-05 1.65741e-05 1.29511e-05 ) ; boundaryField { movingWall { type calculated; value uniform 0; } fixedWalls { type calculated; value uniform 0; } frontAndBack { type empty; value nonuniform 0(); } } phi size 24我的第一个问题是， 为什么phi只在internal face上面赋值了， 边界上的face不用计算吗？在 fvm::div(phi, U) 也用到了phi值， 边界上的flux是怎么考虑的？
第二个问题是， 我声明了一个新的fvVectorMatrix
fvVectorMatrix UEqn_time_term ( fvm::ddt(U) );查看源项的值
Info<< " UEqn_time_term.s " << UEqn_time_term.source() << nl << endl;输出结果：
UEqn_time_term.s 16 ( (-5.83396e-06 5.74984e-06 0) (-1.42259e-05 2.43986e-06 0) (-1.42171e-05 -2.34961e-06 0) (-5.92331e-06 -5.6699e-06 0) (-1.00527e-05 2.19123e-05 0) (-2.40783e-05 8.37173e-06 0) (-2.41399e-05 -8.28764e-06 0) (-1.01346e-05 -2.19111e-05 0) (-1.39378e-05 4.90818e-05 0) (-3.80295e-05 1.57254e-05 0) (-3.80131e-05 -1.56865e-05 0) (-1.39231e-05 -4.90753e-05 0) (9.96865e-05 3.28844e-05 0) (7.57897e-05 9.67425e-06 0) (7.57862e-05 -9.65121e-06 0) (9.97123e-05 -3.2868e-05 0) )这个fvmatrix的源项不是应该为零吗？为什么会有值呢？
非常感谢！！
-
Hi,
我的第一个问题是， 为什么phi只在internal face上面赋值了， 边界上的face不用计算吗？在 fvm::div(phi, U) 也用到了phi值， 边界上的flux是怎么考虑的？
phi是通量值，只在面上有值，也就是说只有面有通量。看这个：http://www.dyfluid.com/flux.html 边界处的通量计算而来。要么给定，要么通过边界条件计算。这个fvmatrix的源项不是应该为零吗？为什么会有值呢？
不能是0. 是零的话：AU=0。你的速度都是零。具体的在离散过程中：
对时间的离散调用了显性的值。 -
东岳老师，非常感谢您的回复！
对于第一问题，我又进行了尝试，下面是我的理解，如果有不对的地方，还请您指正。
边界上的phi可以由边界条件得到，在cavity这个算例中，边界上没有法向的速度，因而phi全部为‘0’。而在internal face上的phi，可以由已知的速度场进行差值得到face上的速度，再乘以该face的面积和法向，就得到了phi。对于第二个问题，我又输出了U的值。
Info<< " U " << U << nl << endl;得到
U dimensions [0 1 -1 0 0 0 0]; internalField nonuniform List<vector> 16 ( (-0.00109257 0.00108567 0) (-0.00260853 0.000420926 0) (-0.00260827 -0.000416558 0) (-0.0010966 -0.00108186 0) (-0.00181643 0.00400068 0) (-0.00421015 0.00141225 0) (-0.00421276 -0.00140849 0) (-0.00182014 -0.00400052 0) (-0.0033937 0.00937929 0) (-0.00774594 0.0025887 0) (-0.00774527 -0.00258687 0) (-0.00339316 -0.00937899 0) (0.0210076 0.00646146 0) (0.0145416 0.00159495 0) (0.0145415 -0.00159377 0) (0.0210087 -0.00646063 0) ) ; boundaryField { movingWall { type fixedValue; value uniform (1 0 0); } fixedWalls { type noSlip; } frontAndBack { type empty; } }发现 UEqn_time_term.source() 里面的值对应为速度的值除以时间步长再乘以对应cell的体积。（时间步长为0.001，体积为0.01乘以0.025乘以0.025）
-
这一项吧？对应icoFoam解析公式6的第一项。
这个研究思路很好。小伙子我看好你。:sunglasses: 哈哈边界上的phi可以由边界条件得到，在cavity这个算例中，边界上没有法向的速度，因而phi全部为‘0’。而在internal face上的phi，可以由已知的速度场进行差值得到face上的速度，再乘以该face的面积和法向，就得到了phi。
是的
-
东岳老师，我在研究完icoFoam中的 ‘fvm::ddt(U)’ 和 ‘fvm::div(phi, U)’ 之后又进而研究了‘fvm::laplacian(nu, U)’，有如下几个问题没有弄明白。
为了看的清楚，附上网格编号的图。
我分别定义了不同的 fvVectorMatrix
fvVectorMatrix UEqn_time_term ( fvm::ddt(U) ); fvVectorMatrix UEqn_advection_term ( fvm::div(phi, U) ); fvVectorMatrix UEqn_diffusion ( fvm::laplacian(nu, U) );对他们进行输出之后发现 UEqn_time_term 和 UEqn_advection_term 的diag() 和A() 每个元素都是相差一个cell 的体积
UEqn.A() = UEqn.diag()/mesh.V()但是UEqn_diffusion的diag() 和A()只有在非边界的单元才满足上述关系。如下
UEqn_diffusion.diag 16 ( -0.0002 -0.0003 -0.0003 -0.0002 -0.0003 -0.0004 \\这个满足上面的关系 -0.0004 \\这个满足上面的关系 -0.0003 -0.0003 -0.0004 \\这个满足上面的关系 -0.0004 \\这个满足上面的关系 -0.0003 -0.0002 -0.0003 -0.0003 -0.0002 ) UEqn_diffusion.A dimensions [0 0 -1 0 0 0 0]; internalField nonuniform List<scalar> 16 ( -96 -80 -80 -96 -80 -64 -64 -80 -80 -64 -64 -80 -96 -80 -80 -96 ) ; boundaryField { movingWall { type extrapolatedCalculated; value nonuniform List<scalar> 4(-96 -80 -80 -96); } fixedWalls { type extrapolatedCalculated; value nonuniform List<scalar> 12 ( -96 -80 -80 -96 -96 -80 -80 -96 -96 -80 -80 -96 ) ; } frontAndBack { type empty; } }我的第一个问题是 在UEqn_diffusion中.A()是怎么计算的？
我在输出UEqn_diffusion的时候也发现它与其他的输出结果不同，在网上搜了半天也没有得到答案。
Info<< " UEqn_diffusion " << UEqn_diffusion << nl << endl;输出结果如下：
UEqn_diffusion false true true 16 \\这个应该是diag() ( -0.0002 -0.0003 -0.0003 -0.0002 -0.0003 -0.0004 -0.0004 -0.0003 -0.0003 -0.0004 -0.0004 -0.0003 -0.0002 -0.0003 -0.0003 -0.0002 ) 24 \\这个应该是upper ( 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 ) [0 4 -2 0 0 0 0] 16{(0 0 0)} \\我猜这个是source，这里面没有源项，4by4的网格，应该都是零。如果不对，请您指正。 3 \\从这个开始往下我不知道是什么了，其他的matrix这些项都是零，请问这些值都代表了什么？ ( 4((-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002)) 12 ( (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) (-0.0002 -0.0002 -0.0002) ) 0() ) 3 ( 4((-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0) (-0.0002 -0 -0)) 12{(-0 -0 -0)} 0() )第二个问题 写在了代码的注释里，主要是不明白这个输出的结果里面除了diag(), upper() 和source()之外的那些项是什么？
我在推倒UEqn_diffusion 的系数矩阵的时候发现，由边界条件得到边界上的速度值，那么最上面的滑动边界上的cell （编号为P）
\begin{equation}
S_{e}\nu_{e}\frac{u_{E}-u_{P}}{\Delta x} + S_{w}\nu_{w}\frac{u_{W}-u_{P}}{\Delta x} + S_{n}\nu_{n}\frac{u_{N}-u_{P}}{\Delta y} + S_{s}\nu_{s}\frac{u_{S}-u_{P}}{\Delta y}
\end{equation}里面的$u_{N}$可以根据边界条件差值
\begin{equation}
u_{N} = 2 u_{n} - u_{P}
\end{equation}其中$u_{n}$为上面滑动边界已知的速度。
但是这个边界条件并没有在系数矩阵的diag()中体现出来，这也带来了我的第三个问题，这个边界条件在 UEqn_diffusion 是怎么体现出来的？是否与第二个问题有关？
第四个问题就是也是又第三个问题引出来的，我试着求解了一下UEqn_diffusion，并输出了一下U的结果
solve(UEqn_diffusion); Info<< " U solution " << U << nl << endl;但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同， 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。
我的问题有点儿多，我也会继续查找相关的信息，期待您的解答！
非常感谢！
-
@cfd-china 谢谢您！我是刚开始学Openfoam，希望能多多向前辈们请教！
-
补充一下上面的帖子
solve(UEqn_diffusion); Info<< " U solution " << U << nl << endl;结果如下：
smoothSolver: Solving for Ux, Initial residual = 0.921181, Final residual = 9.32208e-06, No Iterations 11 smoothSolver: Solving for Uy, Initial residual = 0.999778, Final residual = 8.66271e-06, No Iterations 8 U solution dimensions [0 1 -1 0 0 0 0]; internalField nonuniform List<vector> 16 ( (0.0189067 -1.2469e-08 0) (0.0430646 -3.74215e-08 0) (0.0430639 -4.90001e-08 0) (0.0189057 -2.83642e-08 0) (0.0703756 -3.73926e-08 0) (0.153354 -1.08185e-07 0) (0.153353 -1.33266e-07 0) (0.070374 -6.9139e-08 0) (0.179619 -4.88555e-08 0) (0.34663 -1.32978e-07 0) (0.34663 -1.48122e-07 0) (0.179618 -6.64227e-08 0) (0.481091 -2.81724e-08 0) (0.706929 -6.88449e-08 0) (0.706929 -6.63391e-08 0) (0.481091 -2.54916e-08 0) ) ; boundaryField { movingWall { type fixedValue; value uniform (1 0 0); } fixedWalls { type noSlip; } frontAndBack { type empty; } }但这次与前面的 UEqn_time_term 和 UEqn_advection_term 不同， 没能找到其结果和UEqn.A() 与UEqn.H() 的关系。感觉openfoam求解了这个方程，从UEqn_diffusion的系数和source()项为零来看，求解结果应该是零，但是显然是我想错了。如果是$u_{P}$为未知，$u_{W}$等相邻的cell的速度用的是上一个时间步的值，那么线性方程组只存在对角线的值，不需要迭代，直接就能出结果。在这个问题上没有想明白，还希望您能指点一下这个方程是怎么求解的。
-
不好意思信息量太大。
如果对时间离散
fvm::ddt(T)，其离散后的为对角阵。因此如果源项为0，AT=0则T也为零。并且不需要迭代求解。对于对流项fvm::div(phi,T)，离散后非对角阵因此不一定T为0，且需要迭代求解。 -
@李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b，而是指加在网格单元中心的力源项。
-
@wwzhao
即使没有在方程等号右边添加源项的fvMatrix，比如：fvScalarMatrix TEqn ( fvm::ddt(T) );这个矩阵系统依然具有源项。显性离散的值全部进入了
fvMatrix的source()。P.S. 这个问题@mengweilm425 不是在一楼测试过了么？
用一个网格单元测试：#include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, 1 ); fvScalarMatrix TEqn ( fvm::ddt(T) ); Info<< "source" << TEqn.source() << endl; } -
针对@mengweilm425 提到的
A()和diag()。区别是A()=diag()/V；可测试：#include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, 1 ); fvScalarMatrix TEqn ( fvm::ddt(T) ); Info<< "source" << TEqn.source() << endl; Info<< "source" << TEqn.A() << endl; Info<< "source" << TEqn.diag() << endl;//diag() equals TEqn.A() multiply mesh.V() }源代码：
template<class Type> Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const { tmp<volScalarField> tAphi ( new volScalarField ( IOobject ( "A("+psi_.name()+')', psi_.instance(), psi_.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), psi_.mesh(), dimensions_/psi_.dimensions()/dimVol, extrapolatedCalculatedFvPatchScalarField::typeName ) ); tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V(); tAphi.ref().correctBoundaryConditions(); return tAphi; }template<class Type> Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const { tmp<scalarField> tdiag(new scalarField(diag())); addCmptAvBoundaryDiag(tdiag.ref()); return tdiag; } -
@李东岳 我是这么理解的，
fvm::ddt反映的是上一时刻（或迭代步）对当前时刻（或迭代步）的物理量的影响，OpenFOAM把这种影响当作源项处理，因此将离散后的系数放在source()里面。 -
lduMatrix::diag()和fvMatrix<Type>::D()并不等价。相比lduMatrix::diag()，fvMatrix<Type>::D()考虑了边界对系数矩阵的影响，具体代码对应这句：addCmptAvBoundaryDiag(tdiag.ref());。 -
有关diag()和D()严格来说是的。你可以测试一下然后给楼主提供更详细的信息。:sunglasses:
-
@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()。
-
非常感谢大家的回复！我现在的阶段还看不好代码，但是我又对fvMatrix做了一些测试。
还是如下的fvMatrix
fvVectorMatrix UEqn_diffusion ( fvm::laplacian(nu, U) );我发现在解这个方程的时候
\begin{equation}
[A][x]=[b]
\end{equation}矩阵[A] 实际上是 UEqn.D() + UEqn.upper() + UEqn.lower(), 而没有使用UEqn.diag()。UEqn.D()在UEqn.diag()的基础上加上了边界条件对当前网格的影响。上面方程中的[b]，也不是UEqn.source()，而是和边界条件中已知的速度相关的值。我直接求解得到的结果和openfoam算出来的结果是一样的。由于没有真正的去看代码，不知道Openfoam里面有没有相关的function可以直接输出上述$[A][x]=[b]$方程组中的[A]和[b]，还请各位大神指点。
再次感谢大家的帮助！
-
@mengweilm425 $[A][x]=[b]$ 中的 [A] 可以看成是 internal field 离散后的系数矩阵和 boundary conditions 对系数矩阵的影响两部分之和，其中 internal field 离散之后的系数矩阵表示为 upper()，diag() 和 lower()。boundary condition 对系数矩阵的影响则放在 internalCoeffs_ 里面。[b] 也一样，可以分为 internal field 产生的源项 source()，及 boundary conditions 对源项的影响 boundaryCoeffs_。
OpenFOAM中没有直接输出 [A] 和 [b] 的函数，关于如何将 boundary conditions 的影响加到 fvMatrix 里面，你可以参考 addBoundaryDiag，addCmptAvBoundaryDiag 及 addBoundarySource 这几个函数。
@wwzhao 在 关于 phi 和 fvVectorMatrix 的两个问题 中说：
@李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b，而是指加在网格单元中心的力源项。
上面的帖子描述不太准确，确切说 [b] 是 internal field 产生的源项 source() + boundary field 产生的源项 boundaryCoeffs_。
-
@wwzhao 谢谢您的回复，这就明朗多了。我再好好看看这部分代码。:happy:
-
@mengweilm425 您好 看了你的研究过程，觉得很好。我也想这样分步骤运行各个程序，以分清各个量含义。请问您是怎么做到的。我先编译了自己求解器，在文件中加入自己想输出的量，希望在输出中出现想看到的量。结果并没有看到想看的内容。请问您是如何。。
-
@vivian
我猜测你的求解器名字没有替换？如果编译后的是icoFoam2，记得输入icoFoam2而不是原来的icoFoam.. -
@cfd-china 您好，谢谢您的回复。是替换了的，是我给他速度赋值赋错啦。。
-
@vivian @cfd-china 看来我来晚了:happy:
-
厉害了我的哥哥姐姐们
论坛登录问题反馈可联系 li.dy@dyfluid.com