new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
)
这个主要是初始化的作用。这个fvScalarMatrix
得各种系数都是0。所以就是个0
new fvScalarMatrix
(
k_,
dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
)
这个主要是初始化的作用。这个fvScalarMatrix
得各种系数都是0。所以就是个0
这二者都可以得到cs=0.065,Ce=1.048,Ck=0.065
是标准系数。然后$k_{sgs}$也跟这两个系数的比有关。不过取别的值也能达到这个效果。 @coolhhh 在这方面应该有更深的见解。
$$
\begin{array}{l}
\mathbf{F}_y(t)=-\rho\left(g \sin \theta(t)+y \omega(t)^2+z \varepsilon(t)+2 \omega(t) u_y(t)\right) \boldsymbol{j} \\
\mathbf{F}_z(t)=-\rho\left(g \cos \theta(t)+z \omega(t)^2-y \varepsilon(t)-2 \omega(t) u_z(t)\right) \boldsymbol{k}
\end{array}
$$
这个方程里面非粗体的$\omega,\varepsilon$以及$\theta_m,T$你知道是什么么
在循环前植入:
我忽然注意到你写的这个。你要在循环中更新一下这个量。你更新没有。
当气泡尺寸过小时,特别是对于水中自由浮升的微气泡,就下面这篇帖子曾说到的,interFoam算不出来。
大家主要是从壁面效应带来的阻力,解释了微气泡无法浮升。
在远离壁面的时候,也无法上升么?
这玩意挺有意思。不过你的方程里面epsilon、omega、Ur之类不知道是啥,有相关sci可以参考么
在大涡模拟中使用dynamicSmagorinsky模型,模拟水流与颗粒耦合
哪个求解器
我没有这两个求解器啊老铁 ..
@Prometheus10 老铁,你这个帖子我移动到这里了,https://cfd-china.com/topic/6754, 工作日我看看
厉害啊老铁
楼上说的正确。OpenFOAM没有加滤波的就是网格滤波。加滤波的simple是2倍网格,laplace滤波尺度可以变换。请参考无痛笔记
@zhangxy 是一个网格的尺度
A turbulent trip is a device...
厉害了各位老铁。我知道了。感觉这个trip应该不需要植入,如果没有trip的话。那对于没有trip的SA模型,就是
$$
\frac{\p\alpha \rho \tilde{\nu}}{\p t} +\nabla\cdot \left( \alpha \rho \bfU \tilde{\nu} \right)- \nabla\cdot\left(\alpha \rho D_{\tilde{\nu}}\nabla \tilde{\nu} \right)
= \alpha\rho \frac{C_{b2}}{\sigma_{\nu_t}} |\nabla \tilde{\nu}|^2
+\alpha\rho C_{b1} \tilde{\nu}\tilde{S}
-\alpha\rho C_{w1} f_w\frac{\tilde{\nu}^2}{{y}^2}+S_{\tilde{\nu}}
$$
$$
S_{\tilde{\nu}}=-\alpha\rho f_{t2}\tilde{\nu}\tilde{S}
+\alpha\rho f_{t2}\frac{C_{b1}}{\kappa^2}\frac{\tilde{\nu}^2}{y^2}
$$
只需要搞ft2就可以。要简单不少。那个ft1太长了。
$$
\frac{\p\alpha \rho \tilde{\nu}}{\p t} +\nabla\cdot \left( \alpha \rho \bfU \tilde{\nu} \right)- \nabla\cdot\left(\alpha \rho D_{\tilde{\nu}}\nabla \tilde{\nu} \right)
= \alpha\rho \frac{C_{b2}}{\sigma_{\nu_t}} |\nabla \tilde{\nu}|^2
+\alpha\rho C_{b1} \tilde{\nu}\tilde{S}
-\alpha\rho C_{w1} f_w\frac{\tilde{\nu}^2}{{y}^2}+S_{\tilde{\nu}}
$$
$$
S_{\tilde{\nu}}=-\alpha\rho f_{t2}\tilde{\nu}\tilde{S}+\alpha\rho f_{t1}(\Delta\bfU)^2+\alpha\rho f_{t2}\frac{C_{b1}}{\kappa^2}\frac{\tilde{\nu}^2}{y^2}
$$
我对比了一些,主要是多了个上面这个源项。所以OpenFOAM植入这个应该就是所谓的没有$f_{t2}$项与trip项的SA模型。但是这个$\Delta\bfU$还有后面的$\Delta x$怎么理解?
$\Delta\bfU$ is the difference between the velocity at the field point and that at the trip (on the wall)
https://cfd-china.com/topic/6644
你看看这个,我觉得可以用symmetry边界
我的服务器客户之前遇到过这个问题,也是anaconda跟openfoam并行计算冲突,后来他们吧anaconda卸载了就没事了。我没搞过anaconda,是个知识盲区。看看路过的大佬怎么说吧
@wangfei9088 厉害啊,方程12要不要植入搞一下。算是一个低雷诺数SpalartAllmaras湍流模型。
但是,实际上它只是Spalart Allmaras文献里的方程4,只适用于高雷诺数。对于低雷诺数或者更宽的适用范围,应该用方程(12)的。
另外我把你的这句话加到了笔记里面,说的太好了。
这方面属实属于我的知识盲区 老铁
卧槽,神了,老铁你太厉害了,服了!
我把Wray-Agarwal在OpenFOAM-10里面植入进去了,但是算出来的速度不太对。我检查了几遍也没发现哪里植入错了。感兴趣的可以在我的代码上检查检查。
直接去TurbulenceModeling文件夹下wmake,然后去testCase算simpleFoam就可以
下图一个是Wray-Agarwal算的,一个是SpalartAllmaras算的。
https://turbmodels.larc.nasa.gov/wray_agarwal.html
这里面把Wray-Agarwal湍流模型的公式重新写一下,方便在OpenFOAM中植入
$$
S=\sqrt{2}|\bfS|, \bfS=\frac{\nabla\bfU+\nabla\bfU^T}{2}, W=\sqrt{2}|{\bf{W}}|, {\bf{W}}=\frac{\nabla\bfU-\nabla\bfU^T}{2}
$$
$$
\eta=S \max\left(1, \left|\frac{W}{S}\right|\right) ,k=\frac{\nu_t S}{\sqrt{C_\mu}},\omega=\frac{S}{\sqrt{C_\mu}},
$$
$$
\mathrm{arg_1}=\frac{v+R}{2}\frac{\eta^2}{C_\mu k\omega},f_1=\mathrm{tanh}(\mathrm{arg_1}^4)
$$
$$
C_1=f_1(C_{1kw} - C_{1k\varepsilon}) + C_{1k\varepsilon},
\sigma_R=f_1(\sigma_{kw}-\sigma_{k\varepsilon}) + \sigma_{k\varepsilon}
$$
$$
C_{2kw}=\frac{C_{1kw}}{\kappa^2}+\sigma_{kw},
C_{2k\varepsilon}=\frac{C_{1k\varepsilon}}{\kappa^2}+\sigma_{k\varepsilon}
$$
\begin{split}
\frac{\p R}{\p t}+\nabla\cdot(\bfU R) &-\nabla\cdot((\sigma_R R+\nu)\nabla R)=
\\
& C_1 R S + f_1 C_{2kw}\frac{R}{S}(\nabla R)\cdot(\nabla S)-(1-f_1)\min\left(C_{2ke}R^2 \frac{|\nabla S|^2}{S^2}, C_m|\nabla R|^2\right)
\end{split}
$$
\chi=R/v,f_\mu=\frac{\chi^3}{\chi^2+C_w^3},\nu_t=\rho f_\mu R,
$$
$$
\sigma_{kw}=0.72,\sigma_{k\varepsilon}=1.0,C_{1kw}=0.0829, C_{1k\varepsilon}=0.1284
$$
$$
\kappa=0.41,C_w=8.54,C_\mu=0.09, C_m=8.0
$$
我刚才更新了下笔记,你看下4.2.5 RANS-LES 混合模型、DES 模型
还有这个:10.1007/s10494-005-6916-y
alpha的时间离散格式是什么
那就奇怪了,不知道咋回事。你这个方程式没问题的?
开始时间步循环计算alpha的开头,
alpha方程之前的alpha表示t时刻已知的alpha,alpha方程求解之后表示t+dt时刻的alpha
你写在了alpha.H的前面还是后面
相当小了,可以。
用blockMesh来画这个也是个狠人!
@xiezhuoyu 大神能debug出来更狠!
while (runTime.loop())
{
Info<< "Time = " << runTime.userTimeName() << nl << endl;
#include "CourantNo.H"
// Pressure-velocity PISO corrector
{
fvModels.correct();
#include "UEqn.H"
如果不进行动量预测 U都是当前时间步
// --- PISO loop
while (piso.correct())
{
#include "pEqn.H"
}
在这面之后U是下一个时间步的
}
viscosity->correct();
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
你那个我看不太懂,用上面的代码举例是这样
开始循环U是表示当前已知时间步U(n),现在求的是U(n+1),U.oldTime()表示U(n-1)
这个对
我按照$\nu_t=(C_s\Delta)^2\sqrt{2\bfS:\bfS}$的方式植入进入了。从云图来看没啥区别。希望有大佬来详细验证一下是否均衡。下面这个文件,可以在OpenFOAM-10下进行编译,步骤:
wmake
cd pitzDaily
./Allrun
注意算例中的$C_k$的值忘记改了,应该改为0.065
个人疑问是植入:$ \nu _{sgs} =\left ( C_s \Delta\right )^2 \sqrt{2\mathbf{D} :\mathbf{D} }$ 那么式 (1d) 中的k是怎么计算的?这会不会就是两个方法的主要差别?
@coolhhh 应该是从$\nu_t$反推出来的,$k=\nu_t^2/(C_k\Delta)^2$
这方面可以详细讨论下。大佬有没有相关的sci看看 @kcol
这个问题很难搞 玄学
twoSymm(fvc::grad(U))
正解
直接blockMesh,然后buoyantFoam跑就可以,我跑到了20多秒没啥问题,流场看起来也正确。
https://github.com/lhb8125/swOpenFOAM
之前见过他们这个东西。他们吧GAMG里面的smoother添加了PCG。没详细玩过。也不好评价。
blocks
(
//- block0
hex (0 1 2 3 4 5 6 7) (33 33 56)
simpleGrading
(
(
(0.39716 0.45455 3.7975) (0.20569 0.09091 1) (0.39716 0.45455 0.26333)
)
(
(0.39716 0.45455 3.7975) (0.20569 0.09091 1) (0.39716 0.45455 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block1
hex (3 2 9 8 7 6 13 12) (33 365 56)
simpleGrading
(
(
(0.39716 0.45455 3.7975) (0.20569 0.09091 1) (0.39716 0.45455 0.26333)
)
(
(0.01826 0.0411 3.7975) (0.96348 0.91781 1) (0.01826 0.0411 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block2
hex (8 9 10 11 12 13 14 15) (33 41 56)
simpleGrading
(
(
(0.39716 0.45455 3.7975) (0.20569 0.09091 1) (0.39716 0.45455 0.26333)
)
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block3
hex (1 36 37 2 5 39 38 6) (129 33 56)
simpleGrading
(
(
(0.05252 0.11628 3.7975) (0.94748 0.88372 1) (0 0 1)
)
(
(0.39716 0.45455 3.7975) (0.20569 0.09091 1) (0.39716 0.45455 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block4
hex (9 16 19 10 13 20 23 14) (15 41 56)
simpleGrading
(
1
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block5
hex (16 17 18 19 20 21 22 23) (41 41 56)
simpleGrading
(
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block6
hex (27 26 17 16 31 30 21 20) (41 338 56)
simpleGrading
(
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.0198 0.04438 3.7975) (0.96041 0.91124 1) (0.0198 0.04438 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block7
hex (24 25 26 27 28 29 30 31) (41 41 56)
simpleGrading
(
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
//- block8
hex (25 32 33 26 29 35 34 30) (102 41 56)
simpleGrading
(
(
(0.0676 0.14706 3.7975) (0.9324 0.85294 1) (0 0 1)
)
(
(0.26477 0.36585 3.7975) (0.47046 0.26829 1) (0.26477 0.36585 0.26333)
)
(
(0.17054 0.26786 3.7975) (0.65891 0.46429 1) (0.17054 0.26786 0.26333)
)
)
);
你把上面的blockMesh的block网格减少,处理成1万个网格一下,然后用跟我一样的形式发上来
算例网格少的发上来我给你看下
fixedFluxPressure不对么
Cd这个在forceCoeff.C
里面计算出来的,这个我看没有公开的接口,也不知道哪个是先声明的。另外一个方法是是你自己计算一下Cd。
$$
Cd=\frac{\left( \sum p_f\bfS_f + \sum \bfS_f\cdot \tau_f \right)\cdot\mathbf{direction}_{drag}}{(\sum|\bfS_f|) \frac{1}{2}\rho|\bfU|^2}
$$
里面的$\tau$定义在eddyViscosity里面的sigma()
,被称之为R
。不过我不确定能否在codedFixedValue里面去引用(不确定这个前后关系是哪个先声明的)
不太好整。
你把计算域画大,比较远的地方网格比较少,也没影响也不增加计算量。
为了减少网格数量,液池画的较小,模拟的时候,液面高度和曲率发生了很大的变化,
如果总体积固定的话,毛细管里面的液体上去了,外面的液体不可避免的要下降没办法啊
@2019201300 厉害厉害,LES么?
@李东岳 在 Smagorinsky模型系数问题 中说:
你的意思就是,OpenFOAM里面植入的应该是$\nu_{sgs} = (C_s^{of})^2 \Delta^2 \sqrt{ \mathbf{D} :\mathbf{D} } $,其中$(C_s^{of})^2=C_s^2\sqrt{2}$
而不是$\nu_{sgs} = C_s^2 \Delta^2 \sqrt{ \mathbf{D} :\mathbf{D} } $
是这样么
https://bugs.openfoam.org/view.php?id=816
这个确实是个bug,不过OpenFOAM新版已经把动态Smagorinsky删掉了,所以无关紧要了。
你说的没毛病,老铁
可以做到没问题 他那个有调节箭头大小的 好像是叫tip 你找找
那,只要用到phi的地方就会更新phi的值吗,还是有什么代码会暗中计算一下
pEqn.H里面更新了phi
我在程序里写phi=fvc::flux(U),也能正常计算是吗
能
tmp<fvVectorMatrix> tUEqn
你这个方程里面的phi2,应该需要显性的更新一下