钝体建筑扰流的大涡模拟
-
非常感谢 @gtian 提供的模型，最近和一个师弟在模拟TPU建筑时，迎风面风压总是计算整体偏小。我按照提供的模型及边界条件设置再计算了下，有几个疑问，再请教一下：
1. 速度入口目标值选择问题
描述：提供的边界条件选用了
turbulentInlet,参数设定为referenceField uniform (9 0 0); fluctuationScale (0.02 0.01 0.01);，这是一个平均速度为 9m/s 的均匀湍流场。但是TPU实验用的是指数率平均风速剖面，如下图所示。根据TPU数据库，对于建筑模型1:1:4，风速剖面指数为1/4的实验，参考风速为11m/s，即：$U(z)=U_{ref}\left ( \frac{z}{H_{ref}} \right )^{\alpha},U_{ref}=11m/s,H_{ref}=0.4m,\alpha =1/4$
问题：入口速度目标值不一致，计算结果跟TPU的风洞实验对比，是否会不太合适？
2. 平均风压系数的计算问题
描述：平均风压系数计算公式为：$ \bar{C}_{p}= \frac{\bar{p}}{\frac{1}{2} \rho U_H^2} $，其中 $\bar{p}$ 为平均压力，$\rho$为空气密度，$U_H$ 为建筑物高度处风速。对于TPU指数率剖面，本工况取11m/s。
问题：对于 9m/s 的均匀湍流场，$U_H$是如何确定？以及由OpenFOAM导出的压力值，计算风压系数时的分母是否需要除以密度$\rho$？3. 计算网格边界层设置问题
在计算钝体绕流问题时，计算域的底面、建筑物表面的边界层，是否需设置到符合壁面函数的
y+大小？如果没设置边界层，用了nutUSpaldingWallFunction壁面函数，对结果影响有多大？4. 根据 @gtian 提供的边界条件及模型，我们的计算结果
我们采用
Smagorinsky模型计算，其他的边界条件及求解格式都与提供算例一致，也用的是 9m/s 的均匀湍流场。从1s后开始统计平均场，总共计算24s。计算结果如下。可以看出最大平均风压为41Pa，计算最大风压系数约为：$41/(0.5\times 1.225\times 9\times 9)=0.826$。
5. 前期我们模拟TPU114模型的工况设置及结果
前期我们模拟TPU114模型设置，见下载地址：
链接：https://pan.baidu.com/s/161fNMd9disZoZtW-CZgaRg
提取码：rznq用fluentmeshing画的网格，采用groovyBC加载平均风剖面计算，$U(z)=U_{ref}\left ( \frac{z}{H_{ref}} \right )^{\alpha},U_{ref}=11m/s,H_{ref}=0.4m,\alpha =1/4$。
（1）建筑物平均风压结果：最大平均风压系数为 $54/(1/2\times 1.225\times 11\times 11)=0.728$，而实验值的最大平均风压系数为0.874，差距0.15。
（2）平均速度剖面：
（3）OpenFOAM与Fluent结果对比
①迎风面平均风压云图
②2/3H高度处的平均风压结果。OpenFOAM迎风面平均压力系数结果偏小0.1几，但fluent结果与实验值吻合程度非常好。其中两个侧面和背风面，有旋涡脱落，没有加脉动风，因此模拟结果较差是正常的。由于迎风面结果与平均风剖面大小是比较相关的，正常只加平均风，应该模拟结果是最好的。但目前OpenFOAM结果较差。
③除提供的算例设置外，我们尝试过修改网格、边界条件、湍流模型、fvSchemes、fvSolution，但都是这样的结果。可能一直认为加载平均风，OpenFOAM的迎风面平均风压系数应该和Fluent一样很吻合才比较合理，但一直没找到原因。
本来应该是一个很成熟的计算算例，但实在不知道什么因素影响了结果，请各位大佬指点迷津，实在感谢
-
@coolhhh 我之前没有算过风压这类问题，还挺有兴趣的。我对这个u_h，也是非常疑惑。我用的u_h是柱体中心横截面上（u_h）_max。之所以有疑惑，是关于这个u_h的选择位置的，因为我们一般算cp用到是远场处的速度（我之前做过3900圆柱绕流，https://doi.org/10.1063/5.0012358）。事实上，柱体上方会有一个强剪切区，流动分离后会形成回流区。柱体侧边同理。因此，我很疑惑的是u_h 是如何定义的？ 是 h高度处的最大速度吗？如果有流动分离，即使在初始边界上设置了u_h，那么这个u_h也会发生变化啊。反倒是都用，2/3高度处的流场数据（pmean-0）/ （0.5* u_0.667 ^2），我觉得更合理。因为不可压流p就是p/rho，可以看一下文件的量纲，所以就不需要在除密度了。
-
@coolhhh 还有，尽量不要用smagorinsky做亚格子模型，因为这个亚格子模型中的系数C是个定值。他会有两个问题，1.边界上，会overestimate边界附近的粘度。2. 是全场是一个固定系数的话，会忽视掉能量的反级串backscatter现象。one-equation或者dynamic的亚格子模型都比smagorinsky要好一些。对于静压力，我试一下你给的那个边界条件跑一下看看。还有就是关于总用时的。我一般分成两部分，一部分是从初识态到流场稳定（200-300个t*）。第二阶段，统计过程，100-300个t*。t*=特征长度/特征速度。比如H/u*，摩擦速度可以用sqrt（压差阻力/横面积）得到。
-
@gtian 在 钝体建筑扰流的大涡模拟 中说：
@coolhhh 什么是脉动风？怎么加的？
师兄这个我回答一下，平均风指的是随着时间风速不发生改变。脉动风指的是随着时间的改变，风速是会发生改变的。脉动风当前有两种主流的生成方法，一种是域前模拟法类似于风洞试验，前面有一个模拟带，好像该种方法是满足纳维斯托克斯方程。再者是合成法，主流的包括傅里叶合成法，常见的如中科大黄教授的DSRFG，西安大略的CDRFG和华南理工的NSRFG，我师兄也做得这个方面。涡合成法，该方向只知道SEM和罗银师兄的MSEM。我们当前的算例是基于CDRFG生成的脉动风速。
-
@gtian 几个问题个人理解如下：
不可压下得到的p是p/rho
这是我们一直犯的初级错误，也导致了找问题花了一个多月，非常感谢指正。推导的量纲如下：
0文件夹的
p量纲为 $m^2s^{-2}$
$\frac{p}{\rho}$的单位为 $\frac{Pa}{kg/m^3}=\frac{kg/(s^2m)}{kg/m^3}=m^2s^{-2}$因此不可压缩湍流输出的
p本质是$\frac{p}{\rho}$，希望大家都能避开这个初级错误。我之前没有算过风压这类问题，还挺有兴趣的。我对这个u_h，也是非常疑惑。我用的u_h是柱体中心横截面上（u_h）_max。之所以有疑惑，是关于这个u_h的选择位置的，因为我们一般算cp用到是远场处的速度（我之前做过3900圆柱绕流，https://doi.org/10.1063/5.0012358）。事实上，柱体上方会有一个强剪切区，流动分离后会形成回流区。柱体侧边同理。因此，我很疑惑的是u_h 是如何定义的？ 是 h高度处的最大速度吗？如果有流动分离，即使在初始边界上设置了u_h，那么这个u_h也会发生变化啊。反倒是都用，2/3高度处的流场数据（pmean-0）/ （0.5* u_0.667 ^2），我觉得更合理。因为不可压流p就是p/rho，可以看一下文件的量纲，所以就不需要在除密度了。
我们的研究方向为结构风工程，一般计算风压系数用的是无干扰状态下建筑物顶部高度处来流的速度$U_{H}$，在建筑物高度范围内是最大的风速，因此用$U_{H}$计算得到的平均风压系数最大值一般都不超过1。
还有，尽量不要用smagorinsky做亚格子模型，因为这个亚格子模型中的系数C是个定值。他会有两个问题，1.边界上，会overestimate边界附近的粘度。2. 是全场是一个固定系数的话，会忽视掉能量的反级串backscatter现象。one-equation或者dynamic的亚格子模型都比smagorinsky要好一些。对于静压力，我试一下你给的那个边界条件跑一下看看。还有就是关于总用时的。我一般分成两部分，一部分是从初识态到流场稳定（200-300个t*）。第二阶段，统计过程，100-300个t*。t*=特征长度/特征速度。比如H/u*，摩擦速度可以用sqrt（压差阻力/横面积）得到。
（1）目前对湍流模型还只停留在使用阶段，研究不深，非常感谢指导。
（2）现在这个静压力问题找到了，非常感谢。根据我之前的大量尝试，在计算平均风剖面结果时候，基本无论用什么湍流模型、离散格式、矩阵求解器、有无壁面函数，结果基本都是一致的。
（3）计算用时我也没有考虑到那么详细，你的研究非常专业。我之前是这么设置的：前面计算不稳定的时间，设置时间比计算域长度除以参考风速得到的时间再长点，大概让风场能跑完整个计算域；统计时间的长度，用的是模型缩尺比和速度缩尺比，折算到原型大概10min的长度，满足结构风工程一般采样时长要求来定的，当然计算资源足够的话，模拟越长越好。根据之前的计算经验，这么设置计算结果与实验值的吻合程度还算比较好。什么是脉动风？怎么加的？
（1）根据《风工程与结构抗风设计》(武岳, 2014)，脉动风定义如下。
（2）脉动风生成方法可参考综述文章：Inflow Turbulence Generation Methods。
我主要研究的是
Synthetic Random Fourier Method，加载过程例如：先用CDRFG方法生成脉动风速时程数据，用OpenFOAM自带的边界条件timeVaryingMappedFixedValue，这个边界条件是根据入口坐标插值得到每个时刻的瞬时风速实现脉动风的加载，可参考 LES定义入口速度的问题（DSRFG方法）。还有一个问题是关于inlet边界条件的。一般认为入口处是一条大气边界层的速度型曲线。比如：https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-inlet-atm-atmBoundaryLayer.html。 为什么入口的速度是一个指数型呢？
我们的研究领域主要是结构风工程。风工程中对平均风剖面主要分为对数率和指数率。《风工程与结构抗风设计》(武岳, 2014)中的相关描述如下。目前大部分结构风工程风洞实验室都采用的是指数率的剖面风场来做风洞实验，TPU数据库也是用的指数率剖面。
最后，再次感谢大家的讨论和帮助，感谢李老师，感谢CFD中文网
-
- 我们目前模拟的是高层建筑大涡模拟，不同地貌的粗糙度，对应不同的平均风剖面和湍流强度剖面形状。目前直接在入口施加平均风叠加脉动风，使用no-slip壁面，不采取其他操作，近壁面的平均风随着发展可能会有增大趋势，如下图所示。总体上这对于高层建筑的风压模拟结果影响不大。
(引自：Chen, L., et al., Consistency improved random flow generation method for large eddy simulation of atmospheric boundary layer. Journal of Wind Engineering and Industrial Aerodynamics, 2022. 229: p. 105147.)
2.有一种保持粗糙度的方式是地表布置粗糙元，如下图所示。对应的另外种方式就是施加浸没边界方式（也就是体积力源项）。我对这块没有研究，个人的疑惑是如何根据不同的粗糙地表直接选择合适的粗糙元或体积力源项来实现平均风的自保持，还是说需要试错的方式来确定？
（引自：Bervida, M., et al., Synthetic generation of the atmospheric boundary layer for wind loading assessment using spectral methods. Journal of Wind Engineering and Industrial Aerodynamics, 2020. 196: p. 104040.）3.也有文章说合适的壁面函数可以实现平均风的自保持。个人理解壁面函数的作用是用半经验公式近似模拟近壁面较粗糙网格的流动情况。OpenFOAM提供了几种壁面函数，觉得是没有考虑到风工程中不同地貌粗糙度问题。同样的疑惑是，当模拟不同粗糙度的风场，壁面函数如果能起到作用，那壁面函数是否得根据模拟的不同粗糙场地对应修改，然后还需要自定义植入？
- 我们目前模拟的是高层建筑大涡模拟，不同地貌的粗糙度，对应不同的平均风剖面和湍流强度剖面形状。目前直接在入口施加平均风叠加脉动风，使用no-slip壁面，不采取其他操作，近壁面的平均风随着发展可能会有增大趋势，如下图所示。总体上这对于高层建筑的风压模拟结果影响不大。
-
-
