但是RC插值没有离散，因为它是显式计算源项 -- “速度插值+高阶压力导数项”。

基于2011年发表在Journal of Computational Physics的文章Adaptive mesh refinement based on high order finite difference WENO scheme for multi-scale simulations，并根据您文章的网格自适应算法编写了一套计算程序，在计算双马赫反射算例过程中，遇到了一些问题，受困许久，希望大神指点迷津~ 在此诚挚感谢
遇到的问题具体问题如下：
在计算双马赫反射问题时，在Kelvin-Helmholtz rolls两侧流场出现较大的震荡，局部密度云图及网格自适应加密情况如图1所示，震荡区域存在粗细网格过渡,故怀疑在粗细网格过渡中存在算法错误或者代码错误（我已多次检查代码，未发现代码错误，但不能完全排除仍有代码错误存在）。
0_1536922532180_9e8d3f5c-b1c2-4828-bec3-8e472986a4d2-image.png
（a）局部密度云图（红框处为存在流场震荡处）
0_1536922552890_10e26979-d84c-4d89-956b-e81584a6adb5-image.png
（b）局部网格生成
0_1536922574620_04fd03c7-eca8-4851-a7e4-3c342cae3939-image.png
（c）总体密度云图
0_1536922591916_530a5e07-8d83-41c6-b1e3-338c5acb0b3a-image.png
（d）总体网格生成（黑色区域表示网格加密）
图1 局部密度云图及网格自适应加密情况
该程序采用C++语言编写，流场计算使用有限差分法，五阶WENO格式，Steger-Warming通量分裂，2D欧拉方程，时间推进采用3阶Runge-Kutta方法；网格自适应加密方法和加密准则与您的论文有所不同，采用基于原始根网格进行加密，每个网格单元生成3×3个子网格单元，加密准则使用密度梯度(∂ρ⁄(∂x+∂ρ⁄∂y))。一共使用三层网格，在计算域[0,3]×[0,1]中，原始根网格单元数量为300×100。Layer 0表示原始根网格，Layer 2表示最密的网格。Layer 2比Layer 1略大，覆盖Layer 1，两层网格的间隙为3个根网格单元的尺度，如图2所示。以此希望精准捕捉激波，并适当减少网格生成数量。子网格层的时间步长为父网格层的1/3。
父网格层每计算完一步，对父网格节守恒量（ρ，ρu，ρv，E）和守恒量对时间的导数[ρ’，（ρu）’，(ρv)’，E’]进行五阶加权(类似WENO加权)空间插值，求解出子网格层边界的虚网格节点的相应守恒量及其导数。并将前一时间步和本时间步插值出的通量及通量导数使用Hermite插值，以获得三个子时间步的相应变量，作为加密网格层的边界条件，并使用式（1）作为Runge-Kutta时间推进的中间步的边界条件。
0_1536922724975_5a02e4bb-a1eb-4fce-884b-9d73584cc27c-image.png
在子网格层三个时间步计算完成后，父网格层被子网格层覆盖的节点数值被子网格对应同位置节点更新。每一或两个根网格时间步重新生成一次加密网格。
0_1536922743488_5f7e9dc2-9e8c-44f9-88fb-73083e652ec1-image.png
（a）粗网格和细网格
0_1536922779864_ecb19b47-59e2-4210-97fb-4d0f9de2d195-image.png
b）不同层次网格
图2 网格加密示意图
其他相关的情况介绍：
在一次测试中使网格只加密，不减少，结果显示流场基本没有前述震荡，结果如图3所示。此算例应该可以证明该程序计算格式是正确的。然而，由于生成网格实在太多，不具实际应用意义。
0_1536922795688_29716dff-e9d7-4afb-a475-b7a750cebe66-image.png
（a）密度云图
0_1536922809856_e6bd76e2-af3b-4966-9f15-c3e59d18ea3b-image.png
（b）网格生成情况（黑色为加密区域）
图3 网格只加密不减少的计算结果

不知大家在调试程序过程中，是否遇到过类似的情况。
希望大家能够指出
（1）我对算法的理解是否存在错误，特别是粗细网格过渡部分的算法；
（2）解决该问题有何建议；
（3）能够告知粗细网格过渡区域，造成的流场震荡的原因可能有哪些;

再次感谢大家！

@东岳 这个边界条件已经整合到exteranlHeatFluxTemperature里面了，算的也是温度梯度，也就是$q=-\kappa \frac{\partial T}{\partial n}$，这种边界按照FDS技术手册里面应该是DNS的计算方法。
LES计算时给定的方法是$q=h \Delta T$,$h$是对流换热系数。openfoam里面的边界条件要么是固定值要么是梯度，这种差值不知道在of里面能不能实现

补充一下，不知道这些图中的这些数据点是如何提取的，谢谢大家

@shifc 不知道你现在解决了没有，我最近做格栅湍流也碰到了同样的问题，我的解决方法是谱空间的复速度对三个坐标轴必须都是共轭对称的。

你好👋

额，自问自答？

不好意思，那个算例里面求解器我这面编译不了（openfoam-4.x and 5.x）.. 估计后来问题已经解决了吧

这个是理论大坑啊。

如果压力基的求解器，求解的是速度和压力场，那么得用连续性方程来求压力。
如果是密度基的，连续性可以求密度，然后用状态方程求压力。

可压缩PBE方程形式的比较

咦、之前写的这个研究方向已经发出去了

SuperBee好像就会锐化界面。

Yes, SUPERBEE加了下风格式，是叫下风么？Down-wind

哈哈，东岳大神，感谢你做了这个平台 ... 感觉学习CFD就是要不断地重复才能领悟其中的真谛，当然有通俗易懂的资料会更为便捷。

@mohui 多谢支持，CFD任重道远比较耐玩！只要有兴趣，是能研究一辈子的，不像王者荣耀，玩几年到头了 哈哈 :mihu: 资料以后会持续更新，精力有限，稍微慢点

此问题同类于湍流动能耗散率的定义。湍流动能耗散率的定义为：
\begin{equation}
\varepsilon=2\nu_t\overline{S_{ij}S_{ij}}
\end{equation}
经过各向同性假定后有：
\begin{equation}
\varepsilon=\nu_t\overline{\frac{\p u_i'}{\p x_k}\frac{\p u_i'}{\p x_k}}
\end{equation}

Measurement of Lift Force of Air Bubble in Linear Shear Flow of Water Using L-Channel

在这篇文章中，作者研究了单一的不同粒径的气泡在水中的移动方向，主要研究升力系数问题。感兴趣可以看看。

下图就是实验，高速相机监控的不同直径气泡不同的切向移动。

0_1530002399465_捕获.JPG

页数按照pdf的页数，不是书上的页数

上册：

https://www.jianguoyun.com/p/DVcobNoQ9s3ZBhjejV8

73页【直线运动的速度】：导数的直观描述，CFD方程这面全是导数，大部分为偏导数，后面会介绍 75页、76页【公式1-4】 138页【公式3-3】：泰勒公式，CFD很多方程都是基于泰勒公式推导出来的

下册

https://www.jianguoyun.com/p/DQnanmcQ9s3ZBhjijV8

23页：流量计算公式 71-75页：偏导数，CFD这面都是偏导，重点看一下，理解一下 77页提到了个拉普拉斯方程，随便扫一眼 78页全微分，跟CFD的物质导数关联 79页3-3 81页有个例子 112页介绍了梯度，但是感觉描述的不是很简单，可以看看 141-144页 重积分，为之后理解有限体积法做准备 224-225页 曲面积分的定义，为之后理解有限体积法做准备。其中225页有个图，比较形象 230-231页，看高数中对通量的定义 243-245页，高数中明确的对速度散度的定义

为什么我这个推导的反了？应该是大于0 。

0_1529644803725_微信图片_20180622131807.jpg

@东岳 谢谢:baobao:

张量标识预习测验

下面的方程取自各个不同的CFD教材。方程并没有采用代码的形式而是直接截图，这样更加直观的可以感受到文献里面方程这种多样性的书写方式。各位老师和同学需要能够拆分为分量形式。

0_1529193840362_1.JPG

0_1529193880752_2.JPG

0_1529194030583_3.JPG

0_1529194061311_4.JPG

0_1529194127887_5.JPG

0_1529195669196_6.JPG

0_1529195734307_7.JPG

0_1529195783721_8.JPG

0_1529196217099_9.JPG

0_1529196245498_10.JPG

0_1529196323392_11.JPG

0_1529196403565_12.JPG

0_1529196500098_13.JPG

0_1529197325295_14.JPG

大都忽略重力源项，哪位大神给解释下为什么忽略重力源项

重力对流体的作用就是引起水动力学压力，如果是密闭的流体，其不会引起流动变化。因此可以忽略重力项，这样压力边界条件也可以好给出。

一般考虑重力的主要是重力为一种重要的体积力，如密度不均引起的浮力，其中密度不均进一步由温度导致。

还有一种情况要考虑重力是模拟流体的自由落体。

其他情况可能跟具体的应用有关，如旋转系、海洋热风那种普遍都添加重力。

如果不忽略，怎样恰当的处理？

如果考虑重力，通常求解的压力不是先前的压力p，而是压力p减去水力学压力$\rho g h$。从算法层面理解，重力引致下落，压力对这种效果抵消，采用这种压力计算的时候，可以省去这种来回下落-抵消的求解过程。求解出$p-\rho gh$之后，真实的压力即为$p-\rho gh+\rho gh$

@东岳 好的，麻烦您啦！我在一点一点仔细研究下。如果能研究明白，万一以后也有人问，我再做个帖子出来♪(^∇^*)

今天上到了4.4G，也还行

0_1528170382541_捕获.JPG

是否是真实的流场中湍流长度尺度呢？

Kolmogorov是一种假说：所有雷诺数足够高的湍流中，最小的尺度由一个普适性方程来决定，其参数为粘度和耗散率。这个假说建立在若干的假定基础之上，如各向同性。随后，有了Kolmogorov尺度。所以针对你的问题，至少各大教材都认为Kolmogorov尺度是最小尺度。

这个泰勒湍流尺度能否代表RANS求解中的尺度信息呢？

是的

在大涡模拟中的湍流长度如何让计算？

LES里面的湍流长度通过$\Delta$控制。大于$\Delta$被解析，小于$\Delta$被模化。不确定是否可以通过$\Delta$定义LES的湍流长度。

@lsprxd

$h=\frac{q}{T_w - T_b}$

其中：$q$为热流密度，$T_w$为壁面温度，$T_b$为参考温度

这样可以吗？

@东岳 是的，我看的文献大多数是关注速度，温度和涡量场变化，几乎没有涉及压力场变化，于是我编程计算的时候也没有考虑，结果老师问我为什么，我都不知道如何回答了:xinlei:

@winsway_zero 有空整理一下文献，感觉发在CFD界比较合适 :chigua: 类似还有这个帖子 http://www.cfd-china.com/topic/1056/3d-or-2d-les

@cfd中文网 东岳老师您好，我对第三点不是很明白。不太理解这个源项是什么意思，是数值耗散项吗？（如果是数值耗散项的话，这是什么原因产生的呢）可以有具体一点的方程可以参考吗？

我是按照《普朗特流体力学基础》第208页进行翻译的，感兴趣你也可以看看这一节，不过这里没有将太细 :quwan:

@东岳 是我搞错了，向量混合积跟微分算子不一样

在1954年，IBM创立了Fortran语言。

在1955年，多相流中的粒子法就开始进行了。在当时，其被取名为Particle In Cell (PIC)。

PIC由Los Alamos国家实验室的F.H. Harlow提出。在PIC中，拉格朗日粒子通过其坐标（x, y, z）进行定位并且具有质量。有些情况下，拉格朗日粒子还具有多个物理特征。例如在PIC中，拉格朗日粒子可能还具有大小的区别等。

欧拉欧拉那面由于已经在统计学方面已经处理过了，没有噪音。欧拉拉格朗日模拟这面目前有什么文章在处理么？找了找，还没看到。

As an alternative, Lagrangian methods ‘‘discretize” the density function into ‘‘parcels” that are simulated using Monte-Carlo methods. While quite accurate, as in any statistical approach, Lagrangian methods require a relatively large number of parcels to control statistical noise, and thus are computationally expensive.

@李东岳 http://dx.doi.org/10.1016/j.compchemeng.2015.09.008
这篇文章对三种重构方法做了对比。

谢谢各位，解决了:mihu:

之前有点问题，最终如下
\par $C_p$采用4次多项式分段拟合
\begin{equation}\label{equ:NS_Cp}
C_p(T)=a_1+a_2T+a_3T^2+a_4T^3+a_5T^4
\end{equation}
\par静焓
\begin{equation}\label{equ:NS_H}
H(T)=\int_{T_{0}}^{T_{x}}C_p(T)dT=a_1T+\frac{a_2}{2}T^2+\frac{a_3}{3}T^3+\frac{a_4}{4}T^4+\frac{a_5}{5}T^5+a_6
\end{equation}
\par熵
\begin{equation}\label{equ:NS_S}
S(T)=\int_{T_{0}}^{T_{x}}C_p(T)\frac{dT}{T}=a_1\ln{T}+a_2T+\frac{a_3}{2}T^2+\frac{a_4}{3}T^3+\frac{a_5}{4}T^4+a_7
\end{equation}
\par$\bullet$求解静温(已知总温和马赫数)
\begin{equation}\label{equNSUs}
U_s^2=2\left(H(T_{tot})-H(T_{sta})\right)
\end{equation}
\begin{equation}\label{equNSmach}
Mach^2=\frac{U^2}{\gamma(T)R_gT}=\frac{U^2}{\frac{C_p(T)}{C_p(T)-R_g}R_gT}
\end{equation}
\par由(\ref{equ:NSmach})和(\ref{equ:NSUs})得
\begin{equation}\label{equNSTsta}
T_{sta}=\frac{2\left(H(T_{tot})-H(T_{sta})\right)}{Mach^2\frac{C_p(T_{sta})R_g}{C_p(T_{sta})-R_g}}
\end{equation}
$\bullet$求解静压(已知总温、总压和静温)
\par由p等熵过程
\begin{equation}
ds = C_p(T)\frac{dT}{T} -R_g\frac{d p}{p}=0
\end{equation}
\par两边同时积分有
\begin{equation}
\int_{T_{tot}}^{T_{sta}}C_p(T)\frac{dT}{T} =\int_{p_{tot}}^{p_{sta}} R_g\frac{d p}{p}
\end{equation}
\par记
\begin{equation}
S(T_{x})=\int_{T_{0}}^{T_{x}}C_p(T)\frac{dT}{T}
\end{equation}
\par则
\begin{equation}
S(T_{sta}) - S(T_{tot}) = R_g\ln\frac{p_{sta}}{p_{tot}}
\end{equation}
\par那么
\begin{equation}\label{equNS_psta}
p_{sta}=p_{tot}e^{\left(\frac{S(T_{sta})-S(T_{tot})}{R_g}\right)}
\end{equation}

0_1519692297906_捕获.JPG
是否是这样？应力张量我可以理解。我的疑问是速度梯度在笛卡尔坐标系下变换之后直接就变成在极坐标系中了速度梯度了？

@李东岳 大神，这个我真没注意。好像不多。我是土木专业的，不是搞固体力学的。您可以再找找。

@cloud1990 在 VOF求解自由水面问题，求解空气区域吗？ 中说：

Image

Left: Correct jet shape predicted by TruVOF technique used in FLOW-3D. Right: Incorrect jet shape predicted by pseudo-VOF technique used by other CFD codes.

Image
Left: FLOW-3D‘s TruVOF technique predicts jet impingement on wall and some outflow. Right: Pseudo-VOF methods don’t predict realistic jetting of fluid on side walls.

从结果上看，flow3D 所谓的TruVOF是有问题的，air velocity 和　air resistance = 0?　这个可能比较适合粗网格计算，作为一个tricky的存在．但理论上，这个TruVOF是无法解决wind wave的问题．

0_1516841984913_捕获.PNG
动图1：河床沙石冲刷运动

河流的流动会改变河床的形状。底部的沙子以及大石头通常会分为两个区：

顶层的bed-load zone（有研究河海的么，这个怎么翻译？），其中的颗粒会随之移动； 底层：沙子以及石头会被困住，但是会缓慢的爬行；

在碎石较多的河床中，较大的沙石倾向于堆积在顶层，这种现象被称之为armoring。个人感觉之所以称之为armoring是因为较大的沙石趋向于在表面，类似一种武装。

实验表明，在这个区域，大沙石有一个净得向上的运动速度，小颗粒具有一个小的向下的移动速度。 目前对于这种现象的原因尚不清楚。不过现存了一些潜在的理论。例如支持程度比较高的巴西果效应（Brazil Nut Effect）。

0_1516841984913_捕获.PNG
动图2：有人晃荡大米里面的密度差不多的小石头

替代文字
动图3：进行的计算模拟

巴西果效应是指如果把两种颗粒的混合物置于容器中，然后施加外加的振荡，体积比较大的颗粒会上升到表层，而较小的颗粒会沉降到底部。关于这种古老的效应的动力学机制至今仍众说纷纭。

DEM模拟：文中作者采用LIGGGHTS对这个现象进行了模拟，下图可以看出，随着时间的推移，大颗粒漂浮出来。

0_1516841984913_捕获.PNG
图4：Nature Communication中的LIGGGHTS模拟

看起来很有意思 :cheeky:

本文来自Nature Communication

@李东岳 @lsprxd @赵一铭 对于标准壁面函数我是这样理解的：在求解流场的时候还是用原来的k方程和epsilon方程来求解，只不过是将在壁面第一层网格（有可能是网格面上的值，有可能是网格中心的值）上的信息用壁面函数求解的结果来修正，而且壁面函数其实只求解壁面第一层网格的信息，第二层网格之后就是输运方程求解的结果了。不知我这样理解是否正确？

@wwzhao

简单的感受了一下两种赋值方式，效率确实很大提高。

#include <stdio.h> #include <string> #include <chrono> #include <iostream> using namespace std ; class Timer { private: // Type aliases to make accessing nested type easier using clock_t = std::chrono::high_resolution_clock; using second_t = std::chrono::duration<double, std::ratio<1> >; std::chrono::time_point<clock_t> m_beg; public: Timer() : m_beg(clock_t::now()) { } void reset() { m_beg = clock_t::now(); } double elapsed() const { return std::chrono::duration_cast<second_t>(clock_t::now() - m_beg).count(); } }; class SomethingBefore11 { private: int m_array[5]; public: SomethingBefore11() // zero the member array { // If we want the array to have values, we'll have to use assignment here m_array[0] = 1; m_array[1] = 2; m_array[2] = 3; m_array[3] = 4; m_array[4] = 5; } }; class SomethingAfter11 { private: int m_array[5]; public: SomethingAfter11(): m_array { 1, 2, 3, 4, 5 } //zero the member array { } }; int main(){ Timer tBefore11; SomethingBefore11 m_array_before11; std::cout << "Time elapsed: " << tBefore11.elapsed() << ‘n’; Timer tAfter11; SomethingAfter11 m_array_after11; std::cout << "Time elapsed: " << tAfter11.elapsed() << ‘n’; return 0; } [xx OFtutorial0_helloWorld]$ whatAboutThisGuy Time elapsed: 7.506e-06 Time elapsed: 1.47e-07 [xx OFtutorial0_helloWorld]$ whatAboutThisGuy Time elapsed: 8.664e-06 Time elapsed: 1.9e-07 [xx OFtutorial0_helloWorld]$ whatAboutThisGuy Time elapsed: 7.646e-06 Time elapsed: 1.89e-07

最近看到了一个新的并行分解算法，KaHIP，依然存在疑问。不清楚这种不同的挂载方式的数学本质是什么。

http://algo2.iti.kit.edu/kahip/

替代文字

@李东岳 最后还是用了后处理的方法进行计算，又要做实验又要做模拟有点赶，用的是老板之前研究的方法几何平均粒径、几何标准差与浓度的关系，算了浓度再来得到颗粒数关系。不过布朗运动要加进去我感觉得用群体平衡方法。

@李东岳 谢谢李老师

@youmengtian
你用的什么显示器 :cheeky:

基类fvPatchField<Type>的snGrad()定义就是返回face normal gradient... 也就是$\nabla_{\perp}\phi=\nabla\phi\cdot\vec n $

我试试，不过感觉你的算法没毛病，M3守恒。
你考虑过用Wheeler求节点么？比PD要强健。

@anubis 集合重构不求解方程，而是通过通量变化计算体积分数，最后再进行界面重构。得到的体积分数必然是间断（指界面处发生跃变）的一个场。

下图的变量越界了。红的过了上界。蓝色过了下界。
解得纯对流方程，我这几年...一直在和有界和守恒作斗争 :big_mouth: :big_mouth:

......

0_1512380644852_捕获.JPG

各位大佬们，本人现在研究方向是slamming，能推荐一些书籍看看么，并且最近才接触openfoam，感觉要学的好多！

http://www.sciencedirect.com/science/article/pii/S0021999117301493

2011年Vikas提出来一种稳定的高阶矩对流格式之后，本文作者重新提出来一个新的数值格式。

新技术刚追完，马上又出新的。。

@wang123sd 这样啊，那就可能是另外一个问题了，用所谓显式代数应力模型会好一点，计算量没有高很多。用des或les的计算量会比rans至少高一个量级。

@李东岳 参考温度如何确定呢，老师，是环境的温度吗，如果我研究的是室内通风，参考温度是室外空气温度吗，谢谢

@李东岳 666:sunglasses: