CFD中文网

@李东岳 李老师您好，D()是一个volVectorField。主要是发现建立两个fvMatrix时，前面一个fvMatrix会对后面的fvMatrix产生影响，主要体现在源项上。按理说这是两个独立的fvMatrix对象，不应该互相影响啊

各位好，
试了一下以下代码：

	fvVectorMatrix DEqn1
        (
            fvm::laplacian(impKf_, D(), "laplacian(DD,D)")
            - fvc::laplacian(impKf_, D(), "laplacian(DD,D)")
	);

        fvVectorMatrix DEqn2
       (
          fvm::laplacian(impKf_, D(), "laplacian(DD,D)")
          - fvc::laplacian(impKf_, D(), "laplacian(DD,D)")
       );

Info << "source of DEqn1: " << DEqn1.source() << endl;
Info << "source of DEqn2: " << DEqn2.source() << endl;

输出的结果为：

source of DEqn1: 4{(0 0 0)}
source of DEqn2: 4((0 0 0) (0 0 0) (0 0 0) (0 -27000 0))

为什么这两个fvMatrix的源项会有不同？初始化的参数都是相同的。

BTW，fvc::laplacian(impKf_, D(), "laplacian(DD,D)")的值为：

dimensions      [1 -2 -2 0 0 0 0];

internalField   nonuniform List<vector> 4((0 0 0) (0 0 0) (0 0 0) (0 -180000 0));

boundaryField
{
    top
    {
        type            zeroGradient;
    }
    left
    {
        type            zeroGradient;
    }
    right
    {
        type            zeroGradient;
    }
    bottom
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
}

请问大家有用过Cardiff的流固耦合代码solids4Foam吗？请问固体部分的阻尼要如何施加？solids4Foam中的固体力学部分是基于有限体积法来求解的。在传统的固体力学有限元法中，可以通过Rayleigh阻尼来给定结构阻尼，Rayleigh阻尼的阻尼矩阵是质量矩阵和刚度矩阵的线性叠加，solids4Foam中的一个固体模型unsNonLinGeoTotalLagSolid 中给定了一个质量比例阻尼，但没有刚度比例阻尼。
请问怎么在solids4Foam中对固体模型指定Rayleigh阻尼，或者对固体模型给定刚度比例阻尼？

@李东岳 感谢李老师，我去尝试一下

@oitocfd 也是在一个特定区域里面吗？有没有用过这种方法？https://www.cfd-online.com/Forums/openfoam-programming-development/158819-how-accelerate-findcell.html

想要找到一个区域内很多点对应的cell，用findCell()直接寻找要占据大部分的求解时间。请问有没有什么方法让findCell在指定的区域内查找。或者有没有什么加速findCell()查找网格的方法？

额，我应该早点看到这个帖子的，今天面试就问了这个问题，我有点懵，答得不太好。。。

@李东岳 感谢李老师！

原来是这样，感觉很神奇啊。请教李老师，为什么要对源项分开进行显式隐式的离散？

对对对，这样是对的，我之前把应变率和湍动能搞混了，谢谢李老师了！

李老师的“CFD中的大涡模拟”中的公式19，\begin{equation}
\nu_{\mathrm{SGS}}=\rho C_{\mathrm{SGS}} \Delta \sqrt{k_\mathrm{SGS}}
\label{Smag}
\end{equation}
CSGS和delta上面是不是少了平方？看了两本书都是有平方的。还有运动粘度的话是不是多乘了一个密度?

@Cr 原来是在这里啊，感谢感谢！

在各种源项类的声明和实现中经常能看到addSup函数，功能为将源项加入到控制方程中，如在actuationDiskSource.H中：

 //- Source term to momentum equation
virtual void addSup
 (
         fvMatrix<vector>& eqn,
         const label fieldi
 );

那么请问这个addSup函数是在什么时候被调用的？找了很久没有找到。在网上查了一些资料，基本可以确定这个函数是在求解器形成U方程的时候被调用了，如在simpleFoam的UEqn.H中：

     tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
       + turbulence->divDevSigma(U)
      ==
         fvOptions(U)
     );

大家都说fvOptions(U)就用到了addSup函数，可是这里addSup函数是如何被调用的？希望有大神能帮忙解答！

时隔九个月，今天重新看这个帖子，终于搞清楚了，记录一下。
确实是像李老师说的那样，Usource()就是源项，这里是手动加的源项，之前一直以为是系数矩阵的非对角元。
Udiag()代表U方程的对角线元素，来自IduMatrix，而Usource()是源项矩阵，来自fvMatrix。上面的代码之所以会出现Udiag()，其实就是像李老师说的那样，是为了增加对角占优。而多孔隙模型的源项就是从Usource加入的：

Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);

这里Cd就是公式（2）中括号中的部分。

@秋山结月
1.对的，最好另外再建立新的几何体，用这个新的几何体包裹要refine的局部区域。
2.这个问题倒是没有遇到过，是不是stl文件的问题？
3.snap结束后不应该直接没掉啊，即便是没有需要施加边界层的区域，也会有“No layers to generate ...”，这样的信息。看看你后面有没有这样的提示信息？

个人的一个经验是善用refinementRegions关键字，在需要加密的地方建立refinementRegions，将level提高些，应该可以优化一下网格。还有一个地方是可以试试调整surfaceFeatureDict中的includedAngle，可以优化角度变化比较剧烈处的网格，比如你第一张图那里。另外可以适当将maxLocalCells和maxGlobalCells提高些。这个还是熟能生巧，多试试改变那些参数。
SHM是分步来处理网格的，先是snap，最后才是加边界层。所以边界层网格在最后一个时间步内，不知道你有没有查看最后一个时间步的网格。

@东岳 在 请教porosityModel.C中diag()、source()和psi()的意思 中说：

这几行来看，应该就是人工手动加的源项。对角部分isoCd和源项isoCd部分可以抵消，但这样写可以增加对角占优

“对角部分isoCd和源项isoCd部分可以抵消”，这句话还是不太懂，我还得再消化消化。
但今天试着输出了一些网格的Usource()：
这是空间内存在多孔隙介质时的Usource()值：

这是空间内无多孔隙介质时的Usource()值：

发现确实没有多孔隙介质时Usource()的值要小很多，已经接近0了。

@东岳 在 请教porosityModel.C中diag()、source()和psi()的意思 中说：

http://dyfluid.com/icoFoam.html

感谢李老师的回复！
稍微有一点理解了，比如在icoFoam中这个源项是来自于上一个时间步的值：
\begin{equation}
S_{P}^{n}=\frac{V_{p}}{\Delta_{t}}U_{P}^{n}
\end{equation}
但如果是没有时间项的SIMPLE算法是不是这个源项UEqn.source()就没有值了呢？
主要是因为我的问题是在看porousSimpleFoam的代码时产生的。在porousSimpleFoam中的多孔隙模型DarcyForchheimer里面，源项的表达式为：
\begin{equation}
S=-(\mu d+\frac{\rho \left | U \right |}{2})U
\end{equation}
这一段代码应该是对应DarcyForchheimerTemplates.C中的apply()函数:

void Foam::porosityModels::DarcyForchheimer::apply
(
    scalarField& Udiag,
    vectorField& Usource,
    const scalarField& V,
    const RhoFieldType& rho,
    const scalarField& mu,
    const vectorField& U
) const
{
    forAll(cellZoneIDs_, zoneI)
    {
        const tensorField& dZones = D_[zoneI];
        const tensorField& fZones = F_[zoneI];

        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];

        forAll(cells, i)
        {
            const label celli = cells[i];
            const label j = this->fieldIndex(i);
            const tensor Cd =
                mu[celli]*dZones[j] + (rho[celli]*mag(U[celli]))*fZones[j];

            const scalar isoCd = tr(Cd);

            Udiag[celli] += V[celli]*isoCd;
            Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);
        }
    }
}

代码的最后

            Udiag[celli] += V[celli]*isoCd;
            Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);

给Udiag()和Usource()附了值，应该就是代表对U方程的修正，将源项的影响传入到U方程中，这样的理解是否正确？那在赋值之前，对于SIMPLE算法，Usoure()的值应该是为0的，是这样的吗？

啊，尴尬，题目写错了，应该是“请教DarcyForchheimer.C中diag()、source()和psi()的意思”。correct函数是在DarcyForchheimer.C中，不是在porosityModel.C中。正文中对应的地方也错了。

在OpenFoam的多孔介质模型中，代码porosityModel.C定义了correct函数，其中提取了U方程的以下三项内容：

const volVectorField& U=UEqn.psi();
scalarField& Udiag= UEqn.diag();
vectorField& Usource=UEqn.source();

diag()能够理解，应该是U方程系数矩阵的对角元素，那么source()和psi()呢？
李老师2012年在cfd-online里回复的一篇帖子里（The meanings of diag( ), source( ), psi( ) ）提到source()是代表系数矩阵的非对角元素
看了correct函数的代码也确实感觉source()是非对角元素，那为什么非对角元素要取一个“源项”的名字呢？而且对于psi()的意思也不太清楚。
关于多孔隙介质模型中源项的处理还有很多不理解的地方，希望能在这里向大神们多学习。