CFD中文网

史

@李东岳可以把日常的流体都看成稀薄流体呀。从各种benchmark来看的话，精度还不错，处于同一个水平。
现在主要的问题是，显式LBM并行效率特别的高，然后我就寻思着有限体积也可以显式呀，谁怕谁。结果性能卡在了压力泊松方程的求解上

史

有没有可以不通过压力修正方法来求解CFD的算法（LBM除外）？
事情是这样的：最近用PISO方法和一个做LBM的同事PK 求解CFD的性能，遇到了瓶颈：传统的投影法在求解压力泊松方程的时候会花费大量的时间。想知道有没有方法可以绕过求解压力-泊松方程？
我的一个思路就是：在不可压缩流体中引入较小的压缩性，然后写出压力和密度的状态方程，从而壁面压力方程的求解，不知道可不可行。
希望老铁们出谋划策，不想被LBM锤

史

@hongjiewang 在界面相变中说：

afvm::Sp(vDotvmcAlphal, alpha1) + avDotcAlphal
从你给的代码里面，没看出什么问题，可能是其他方面可以出问题了
相含量的变化不一定只和源项有关，还可能和宏观流动、求解误差等相关
建议跟踪一下你给的源项这个东西的值

a*fvm::Sp(vDotvmcAlphal, alpha1) + a*vDotcAlphal

如果这个值没有问题，那看一下你写的代码，看看你代码计算的过程是不是和现实中的物理过程有偏差。比如说这个地方相分数为1，但是你的源项可能是一个正值。因为你用他的散度去表示界面处的相变过程，这个表示本身就是不严谨的。
话说回来，你为啥不用相分数，来判断是否为相界面？代码如下：

forAll(a,celli)
{
    if (alpha1[celli] > 0.001 && alpha1[celli] < 0.999)
    {
        a[celli]=1;
    }
    else
    {
        a[celli]=0;
    }
}

史

@李东岳李老师，好久不见！最近重温你写的interFoam求解器的解析，里面有一句话

VOF作为网格依赖类求解器，加密网格研究网格无关性至关重要

这里所说的网格依赖具体的标准是什么？因为所有的数值模拟结果都和网格相关，但不是所有的求解器都叫网格依赖类求解器。所以我觉得应该要满足一定的标准才能叫网格依赖类求解器，那么有没有什么文献定义过赭红标准？还只是大家经验之谈，只是觉得VOF模型受网格影响大一些？

史

@hongjiewang 能帮到你我也很开心，继续加油~

史

@hongjiewang 你的a是一个volScalarField，而test是一个surfaceScalarField，两个变量不匹配。想办法把surfaeScalarField转化成一个volScalarField试一下

史

@Yu_Tian 都不需要，你只需要添加Ue的相应插值格式和求解器参数即可。也不用添加一个初始场，因为一开始你设置的是NO_READ属性，程序不会读初始场

史

@Yu_Tian 感谢分享，非常棒！学习了！这个失真的问题肯定会存在，只要不影响你想讨论的结果，都是可以接受的

史

@Yu_Tian 我这个图就是局部网格加密（AMR）结果的后处理图。对于局部网格加密算例，流线图没有影响，矢量图确实会在加密网格区域比较密集。矢量图的这个问题我也没找到比较好的处理方法，如果你有比较好的方法可以分享一下

史

@Yu_Tian 在createFields.H文件中，添加变量声明

volVectorField js
(
    IOobject
    (
        "js",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    -sigma*fvc::grad(Ue)
);

在求解完Ue后，更新一下js

js=-sigma*fvc::grad(Ue);

我一般这样通过求解器实现的。
也可以用postProcess实现，但是需要具体看一下他的源码，这个我不是很熟悉
放一个我之前写过的一个求边界“axis”上的最大速度，这个是参考之前版本的calc方法写的，后来calc全部整合到postProcess中了，但是在OpenFOAM 6里面，按照calc的思虑一样可以用

#include "timeSelector.H"
#include "calcType.H"

void calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
    IOobject UHeader
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ
    );

    if (phiHeader.headerOk())
    {
		volVectorField U(UHeader, mesh);
		forAll(U.boundaryField(), patchi)
		{
			const fvPatchVectorField & pU = U.boundaryField()[patchi];
			const fvPatch & pU_B = pU.patch();
			if(pU_B.name() == "axis")
			{
				Foam::Info<< "mag(U) max : " << max(mag(pU)).value() << Foam::endl;
				break;
			}
		}
	}
	else
	{
		Foam::Info<< "No U exists!" << Foam::endl;
	}

    Info<< "\nEnd\n" << endl;
}

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"


    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
   
        Foam::Info<< "Time = " << runTime.timeName() << Foam::endl;
   
        mesh.readUpdate();
   
        calc(args, runTime, mesh);
   
        Foam::Info<< Foam::endl;
    }
	
    Foam::Info<< "End\n" << Foam::endl;
   
    return 0;
}

史

@xiaolin leftEdge没有意义，leftObstacleEdge（rightObstacleEdge）是指运动的左（右）侧边界的位置
这个类其实是一个不太好用的类，你自己要改一下。这个建议你先去看一下源码，了解里面的操作流程，然后修改一下就可以自己用了，亲测可以

史

@Yu_Tian 抱歉，前面的公式少打了一个负号-，不过你已经找到这个错误了。
postProcess这个东西我还没用过，有空研究研究，你看看下面这种书写格式行不行

postProcess -func "grad(-1.0*Ue)"

当时我直接写在程序里面了，把他当成一个结果输出了

史

@xiaolin 我刚刚看了一下，我用的是OpenFOAM 6，前面写错了抱歉。
下面是movingConeTopoFvMesh中对应的构造函数（初始化函数）

 Foam::movingConeTopoFvMesh::movingConeTopoFvMesh(const IOobject& io)
 :
     topoChangerFvMesh(io),
     motionDict_
     (
         IOdictionary
         (
             IOobject
             (
                 "dynamicMeshDict",
                 time().constant(),
                 *this,
                 IOobject::MUST_READ_IF_MODIFIED,
                 IOobject::NO_WRITE,
                 false
             )
         ).optionalSubDict(typeName + "Coeffs")  //在movingConeTopoFvMeshCoeffs中寻找对应参数
     ),
     motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
     motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
     curMotionVel_
     (
         motionVelAmplitude_*sin(time().value()*pi/motionVelPeriod_)
     ),
     leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
     curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
     curRight_(readScalar(motionDict_.lookup("rightObstacleEdge")))
{
    ...
}

从构造函数可以看出来，所有的动网格参数要从一个movingConeTopoFvMeshCoeffs的字典中查找获得，所以，在设置动网格参数的时候，这个字典是不能少的，其对应结构如下：

movingConeTopoFvMeshCoeffs
{
    ...
}

你好像缺少了这个结构

史

@DY大世界久闻snappyHexMesh大名，遗憾从来没用过，平时计算的季和模型都比较简单抱歉，估计忙不上你什么忙
你在文中说网格和计算结果都会略有不同，有没有图可以展示一下他们的不同？

史

@xiaolin 你用的是哪个版本？我用OpenFOAM 7，这个动网格设置大约是这个风格

dynamicFvMeshLibs ("libmovingConeTopoFvMesh.so");

dynamicFvMesh      movingConeTopoFvMesh;

movingConeTopoXXXFvMeshCoeffs
{
    motionVelAmplitude      (0 0.00006 0);
    motionVelPeriod         1;

    //leftEdge                0;
    leftObstacleEdge        0.2;
    rightObstacleEdge       1.5;

    left
    {
        minThickness        2.5e-3;
        maxThickness        7.6e-3;
    }

    right
    {
        minThickness        1.2e-3;
        maxThickness        2.8e-3;
    }
}

史

@Yu_Tian 电场强度是电势（电压）的梯度

E = fvc::grad(Ue)  // Electric field intensity
j = sigma*E        // Current density, sigma--electric conductivity

史

@Yu_Tian 你说的是电场的streamline是么？是下图这种效果么？

这个直接用电场强度矢量做流线图就可以得到，不用再求导了

史

@东岳你这一提醒，我发现我上面的好像错了。这里说的加速度应该指的是物质导数，物质导数和局部导数有如下关系

我上面写的fvc::ddt(U)只是其中的局部导数。
那问题来了，在OpenFOAM中怎么求物质导数？下面的公式求物质导数看起来有点土

DUDt=fvc::ddt(U)+fvc::div(phi, U)

史

@DY大世界加速度是速度对时间的导数，应该是ddt(U)。grad(U)是空间的导数，相当于速度的梯度

史

@tangcheng 这种时候就不能忽略粗糙度和震动了，把这些因素都考虑进去再进行计算。这种时候两个理想平板之间的雷诺数计算方式已经不适用于你说的这种请款的计算了

CFD中文网

史浩

帖子