给定速度压力场,推算湍流粘度场
之前干过一个事,是“给定湍流粘度场,反推速度压力场”:
现在想做一个相反的事,就是“给定速度压力场,推算湍流粘度场”。
为什么要做这个事情?
比如我们现在通过测量,得知一个流体的速度压力场分布,这个分布是测量值,我们认为测量足够精确,那么这个测量值就是真值。既然得到了真值,我们当然希望通过真值来修正我们现有的湍流模型,那么最直接的一个问题,就是此时的雷诺应力是如何分布的?或者说,湍流粘度场是如何分布的?
假如测量值非常细致,DNS那种的,当然就直接求解雷诺应力了。但当测量值网络比较粗,不能直接求解时,该如何推测呢?这就引出了本次的工作。
如何做这个事情:
大致设计流程如下:
- 选定flow configuration,也就是跑什么case,几何形状如何,边界条件如何;
- 生成一个假想的测量值A,作为真值(倾向于用LES生成);
- 将这个真值场A normalize到即将计算湍流粘度的网络上,生成新场B;
- 将新场B带入守恒方程中,求解湍流粘度场。(重点)
- 验算,将湍流粘度场带入守恒方程,反推速度压力场,看看结果是否与原本一致。
我们这篇文章重点记录第四步的代码实现过程,因为其他三步都很简单,代码计算都在第四步。
现在开始做吧:
- 选定flow configuration
选定openfoam里自带的lid-driven cavity case, 几何形状之后根据参考实验配置调整。 - 生成一个假想的测量值A,作为真值
这个测量值A需要满足: - 比对实验结果
- mesh convergence study
参考了几篇文章中的实验结果:
这篇数据不充分,lid velocity和viscosity都没写,气死我了,浪费时间。
又找了一篇,唯一有价值的结果,参考的是上面那人的工作。
顺着找过去,原来那人还写了篇文章,里面数据清楚些。
结果还是没找到lid velocity和viscosity,只能靠自己猜了。工作流体是水,那么根据 这个表:
我感觉就是1e-6了。
文中作者说Re=3200
那么速度就是3200×1e-6/0.15=0.0213m/s,每秒钟两厘米,比较合理。
搞不懂这点东西为什么要藏着掖着,为啥不直接写出来。
我使用20×20×20,40×40×40,80×80×80三种网格运行LES模拟,并对比实验结果:
绘制区域来自下图的左下角四分之一区域:
可以看出结果不是很好,无论哪种网格都没有将试验结果很好match上,可能原因有:
- 模拟使用的几何形状不完全相同,工质粘度也不相同;
- 可能网格还不够细。
但我暂且使用80×80×80网格的结果A作为作为真值,原因如下:
- 80×80×80网格的运行速度就已经很慢了,我的电脑四核要跑一晚上才能达到稳态;
- 改变几何形状和工质也需要时间,我周五就要做个报告,来不及了。
3. 将这个真值场A normalize到即将计算湍流粘度的网格上,生成新场B;
即将计算湍流粘度的网格,我这里设定为20×20×20,也就是说,我们需要将一个80×80×80网格上的结果,投影到20×20×20的网格上。
想了想,80×80×80网格也来不及了,还是用80×80的二维吧。
如图所示的瞬态问题,先用lid velocity为1m/s运行到稳态,然后提高lid velocity到1.1,再次运行到稳态的这么一个过程。需要从80×80的网格投影到20×20的网格上。
先考虑openfoam里自带的mapfield功能,参考我之前的研究:
2.3和2.4版本的openfoam有三种map方法:
direct
: the meshes are assumed to be of identical topology, with one-to-one correspondence between cells, but perhaps with different addressing, e.g. order of points, faces, cells; mapping is volume conservativemapNearest
: cell values on the new mesh are sampled based on the nearest cell to the source mesh; not volume conservativecellVolumeWeight
(default): volume weighted for internal cells, using face area weight AMI for boundaries.
仅仅按照字面意思理解,direct肯定是不行的,direct只能用于两个网格一模一样的情况; mapNearest也是不行的,程序会寻找最近的点上的值而不做插值/平均处理。cellVolumeWeight看上去可以。我希望的平均方式是下图中的cell to cell方式,也叫做ensemble average。
我们运行一下试试吧:
比对20×20的左上角第一个网格的值,是否等于80×80的左上角4×4网格内16个值的ensemble average。
该网格的值为:
(0.000452898 0.0490641 0)
待map的16个网格的值分别为:
(0.102966 0.186011 0)
(0.119124 0.101878 0)
(0.174129 0.059435 0)
(0.207459 0.0416075 0)
(-0.00165368 0.056845 0)
(-0.0324625 0.0327002 0)
(-0.0167796 0.0270086 0)
(-0.00976475 0.0264007 0)
(-0.0256858 0.0709823 0)
(-0.0619708 0.0515074 0)
(-0.0754804 0.0393454 0)
(-0.0895898 0.0328625 0)
(-0.0239223 0.00882332 0)
(-0.067267 0.0133076 0)
(-0.087264 0.0186364 0)
(-0.104591 0.0176744 0)
用excel做平均,结果正是粗网格上的(0.000452898 0.0490641 0):
可见mapFields默认的cellVolumeWeight方法可行,是我们需要的ensemble average方法。
4. 将新场B带入守恒方程中,求解湍流粘度场。
思考了一下,需要重写solver,因为这个case我用pisoFoam跑,所以就修改pisoFoam solver。
重写的思路:
- 去掉time loop,让算法只读取某个特定时间的速度场压力场;
- 带入守恒方程中,创建一个新变量nut_new,让solve函数自己求解;
- 算完带回守恒方程验算,令输出守恒方程的residual,看看是不是都为0;
想了想,solve函数有没有这个能力啊?去看看东岳老师的博客。
……找不到了,东岳老师删得真快。
看了下代码,应该solve是有这能力的,反正先试试,这里mark一下,不对再回来排查。
修改了一下,核心语句应该是下面这几句:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UEqn.relax();
solve(UEqn == -fvc::grad(p));
看了一下,UEqn.relax()似乎没必要了,因为只对一个时间点读取数据,所以 Under-Relaxation等于啥也没做。
fvm::ddt(U)也可以不要了,也是因为只有一个时间点的数据,所以这项肯定也等于0.
fvm::div(phi, U)比较好改,改成fvc::div(phi, U)就行了,从矩阵改成常数列。
fvc::grad(p)不用改。
turbulence->divDevReff(U)要改,而且要大改。
首先需要知道divDevReff调用了什么函数,展开是什么样的:
之前的工作已经研究过了,结论是:
divDevReff(U)等价于:
- fvm::laplacian(nu_t_fix+turbulence->nu(), U)
- fvc::div((nu_t_fix+turbulence->nu())*dev(T(fvc::grad(U))))
这里把nu_t_fix换成我们需要的nut_new,然后稍加修改:
- fvm::laplacian(nut_new+turbulence->nu(), U)
- fvm::div((nut_new+turbulence->nu())*dev(T(fvc::grad(U))))
等等,不对劲啊?momentum equation一共是3个方程,解速度3个未知数正好,现在不解速度改解nut,不就变成3个方程解一个未知数了么?这样能解么?解出来对么?
(待解决)想不到特别好的办法,目前的想法是先令nut_new为一个1×3的vector,强行扩维,相当于对于每个维度的方程解一个nut_new。
按照这个思路,我们再看看原方程要怎么改。
- fvm::laplacian(nu_t_fix+turbulence->nu(), U)
拉普拉斯项好改。
- fvc::div((nu_t_fix+turbulence->nu())*dev(T(fvc::grad(U))))
后面这项不好改。
印象中incompressible flow似乎没有后面这项,思考后面这项能不能不要了。
参考CFDonline里对后面这项的说法:
再看看wiki里对同样步骤的推导:
非常奇怪,wiki里直接把后面这项去掉了。我看了好久CFDonline的解释,还是不得要领,感觉文中的意思似乎是说这种处理方法是为了OpenFOAM里的代码方便?
又查了查书 《turbulence flows (pope)》page 16-17:
按照这里的分析,也是可以去掉的。
我手推了一下:
理论上,确实是可以去掉的。
不清楚为啥OpenFOAM里要保留这一项,可能的原因:
- 代码习惯(能使用预先编译好的function,节省力气)
- 更好的收敛性
总而言之,为了代码方便,我决定大着胆子删了,以后出了事以后再说。
现在代码为:
fvVectorMatrix UEqn
(
+ fvc::div(phi, U)
//+ turbulence->divDevReff(U)
- fvm::laplacian(nut_new+turbulence->nu(), U)
);
solve(UEqn == -fvc::grad(p));
还是不太对。
- fvm::laplacian(nut_new+turbulence->nu(), U)
这项需要改改,现在这个写法变量还是U不是nut_new。需要改成类似下面这样的东西:
volVectorField laplacianU = fvc::laplacian(U)
- fvm::multiple(laplacianU, nut_new+turbulence->nu())
就是不知道openFOAM里有没有fvm::multiple这样的函数。
(想过能不能直接让Ax=b的b直接除以A,但想想不行,边界条件上不行,OpenFOAM里当作变量储存的matrix A和vector b都是不考虑边界条件的,最后solve的时候再把边界条件加上去一起求解,所以最好还是使用OpenFOAM里自带的function)
没找到。
怎么办。
想想能不能直接给fvVectorMatrix赋值,比如这样:
fvVectorMatrix UEqn
(
+ fvc::div(phi, U)
//+ turbulence->divDevReff(U)
//- fvm::laplacian(nut_new+turbulence->nu(), U)
- fvc::laplacian(turbulence->nu(), U)
);
UEqn.A() = - fvc::laplacian(U)
solve(UEqn == -fvc::grad(p));
A代表对角线diagonal element,H代表所有非对角线元素( 参考),我觉得行。(可能量纲上要再想想,不过问题不大,对照其他项总能看出来)
想得差不多了,试试吧:
结果问题很多,主要是变量类型不匹配之类的,非常复杂。还是回到之前的想法,再找找有没有fvm::multiple这样的函数。
我们需要找的其实是一个很简单的函数,就是需要在OpenFOAM里添加一个a x项,a是常数场,x是未知数变量场,为什么会这么难呢?我想起来以前有人问过在OpenFOAM里加一个关于速度变化的source term,找找看他们是怎么解决的吧。
想了想,两个问题不是一个问题,速度变化source term的问题很好解决,就直接加一个fvc格式的a ×U就行了,因为速度会随时间迭代直到收敛,而这里的nut_new必须一次算完。
烦死了,我实在没有办法,就老老实实自己编code算吧。
volVectorField b = - fvc::grad(p)
+ fvc::laplacian(turbulence->nu(), U)
- fvc::div(phi, U);
volVectorField a = - fvc::laplacian(U);
volVectorField x = a;
x.rename("nut_new");
forAll(mesh.cells(), cellI)
{
for (int i=0;i<2;i=i+1)
{
x[cellI][i] = b[cellI][i]/a[cellI][i];
}
}
Info << "b =" << b;
Info << "a =" << a;
Info << "x =" << x;
x.write();
运行之后的结果:
配合速度图:
运行结果可以看出在最靠近壁面的位置需要做出的修正是最多的。
5. 验算,将湍流粘度场带入守恒方程,反推速度压力场,看看结果是否与原本一致。
带入守恒方程,得到新速度场,与原场相比变化不大。
现在看看只用粗网格跑出来的结果:
对比一下x=0.5的位置上,速度magnitude在y轴上的分布:
可以看出,原本粗网格模拟结果(红线)相对真实值(绿线和蓝线)是十分不准的;经过我们添加湍流粘度场之后,预测结果(紫线)就比较准确了。
但预测结果还是跟真实值有一定距离,我理解误差来源可能是:
- 边界条件没做好,还有些问题;
- 算出来的nut _new自带边界条件,可能这里有问题;
- 离散化的scheme有些问题;
- 细网格没有用去掉第二项的方程跑。
理论上预测结果(紫线)应当与真实值(蓝线)完全重合,现在没重合就是代码还有问题,需要继续改。先上传一下,以后再改。