@SHUKK
你需要 libtorch,这个是C++库,提供了跟 pytorch 一样的功能。也就是说,用pytorch训练的模型,可以做到利用 libtorch 提供的 API 读入到任何 C++ 程序中来用。
针对Open FOAM,这个项目 https://gitlab.com/tmaric/openfoam-ml 可以作为入门参考。
xpqiu
帖子
-
-
@SHUKK
你的 g(i) 是用训练好的模型来计算对吧,那么其实你只需要用C++把你计算g(i)相关的代码写出来,写成湍流模型的成员函数,然后再在correctNonLinearStress里面调用这些函数就可以了。虽然g(i)不能写出显式的公式,但是也是代入一些流场变量,然后输出g(i)的数值吧,这样肯定也是可以写出可以调用的函数的。 -
@SHUKK
你可以参考一下这个项目:https://github.com/furstj/myTurbulenceModels里面有非线性湍流模型的例子,比如这个 EARSM 模型,https://github.com/furstj/myTurbulenceModels/tree/main/turbulenceModels/RAS/EARSM
你可以在这个 EARSM 模型的基础上来修改,基本上只需要其中的 correctNonlinearStress 函数就可以了。
-
@李东岳
只是知道用了这个函数就会引起计算错误,但是不确定是不是我其他代码有问题。
希望有大佬路过进来看看这个 bound 函数的实现到底为什么要这样。 -
不知道为什么这个 bound 函数里面要用到 average 这样的操作。
之所以关心这个问题,是因为这个 average 给我导致一个严重的bug -
很多有界变量在计算过程中为了防止数值问题,需要限定下限值,在OpenFOAM中一般都是使用bound 函数,代码见:src/finiteVolume/cfdTools/general/bound/bound.C。然而这个 bound 函数的实现却有些让人摸不着头脑:
Foam::volScalarField& Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound) { const scalar minVsf = min(vsf).value(); if (minVsf < lowerBound.value()) { Info<< "bounding " << vsf.name() << ", min: " << minVsf << " max: " << max(vsf).value() << " average: " << gAverage(vsf.primitiveField()) << endl; vsf.primitiveFieldRef() = max ( max ( vsf.primitiveField(), fvc::average(max(vsf, lowerBound))().primitiveField() * pos0(-vsf.primitiveField()) ), lowerBound.value() ); vsf.boundaryFieldRef() = max(vsf.boundaryField(), lowerBound.value()); } return vsf; }
一般的想法是,直接把小于 lowerBound 的重置为 lowerBound 就行了。但是 OpenFOAM 的实现是,对于某个变量v,当 0 < v < lowerBound的,令 v=lowBound (与预期一致);当 v < 0 时,
pos0(-vsf.primitiveField())
为 1,也就意味着用如下公式重新计算了下限值fvc::average(max(vsf, lowerBound))().primitiveField()
这种做法有什么数学上的解释吗?
-
@李东岳
参考文献:
Knacke, T. (2013).
Potential effects of Rhie & Chow type interpolations in airframe
noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed):
Accurate and efficient aeroacoustic prediction approaches for
airframe noise, VKI LS 2013-03. -
@Chen_hao
当这个模型用于多相流时,IOobject::groupName("omega", alphaRhoPhi.group()) 会变成类似 omega.air 这样的名字。 -
用 topoSet,把这些区域内的网格提取到 cellSet,提取过程会告诉你网格数量。
-
meshQualityControls 部分把 relaxed 里面的质量控制参数再放宽一些。然后 nRelaxedIter 设置为 0。
-
还可以用 topoSet 结合 createPatch 来实现,不过这样得到的面不一定是精确地 1/2,因为createPatch 只是从已经画好的网格上选一些面来生成一个新边界。
-
@soulx7 参考这个算例:
tutorials/incompressible/pisoFoam/LES/pitzDailyMapped -
@wg0632
Gauss LUST grad(U) 当中的 grad(U) 对应一个梯度计算格式,会从 gradSchemes 部分取查找。如果 gradSchemes 部分定义了 grad(U) ,则使用该部分的 grad(U) 对应的格式。如果gradSchemes 部分没有定义grad(U) ,则会使用 gradSchemes 部分的 default 对应的格式。如果 default 也没有则报错。所以,Gauss LUST grad(U) 和Gauss LUST unlimitedGrad(U) 的区别取决于 gradSchemes 里面如何定义的 grad(U) 和 unlimitedGrad(U)。
-
@BznW
direction 是截取的方向,比如 1 0 0 表示沿着 x 方向分块,nBin 是分块的数量,cumulative yes 表示输出累积结果。
因此,上述图片能达到的效果是,沿着 1 0 0 方向将部件分成 20 等份,然后会分别输出这20份的受力。因为 cumulative 是 yes,所以,第一个输出是第一份的受力,第二个输出是第一份和第二份累积的受力,因此类推。 -
A "trip" in the context of computational fluid dynamics (CFD) and the Spalart-Allmaras turbulence model refers to a turbulent trip. Specifically:
A turbulent trip is a device used to trigger a transition from laminar to turbulent flow in a fluid dynamics experiment or simulation. Fluids can flow in either a smooth, laminar state or a chaotic, turbulent state depending on conditions like velocity, viscosity, and surface roughness.
In CFD using the Spalart-Allmaras model, adding a "trip term" attempts to model the effects of an actual physical trip device that would be used in an experimental flow setup. This trip modeling triggers the transition to turbulence in the simulation.
So in the sentence you referenced, the "trip term" refers to an addition made to the Spalart-Allmaras equations that mimics the effect of a physical turbulent trip, forcing a transition from laminar to turbulent modeling. This trip modeling can be important to accurately replicating real experimental conditions in the CFD simulation.
-
@洱聿 集群如果使用作业调度软件来提交作业的话,应该生成日志文件(类似 xxx.out and/or xxx.err),里面可能有什么报错信息。
-
@Tens 在 关于并行中的reduce函数 中说:
List<scalar> np(20,0.0); for (label i=0;i<20;i++) { np[i] += xxxx; //每个时间步累加 reduce(np[i], sumOp<scalar>()); }
这段代码如果以比如 6 核运行,那么如果先不看 reduce 这句,则是这6个核各自都会创建一个长度为 20 的list,然后执行这个循环,最终得到的结果是每个核各自有一份自己的 np 值,如果xxxx是常数,则每个核的np值都一样,如果xxxx是一个不同核数取值不一样的数,则每个核中最终得到的 np 值不一样。
reduce(np[i], sumOp<scalar>()) 的作用是将每个核的 np[i] 全都累加起来(因为这里的第二个参数是sumOp,表示加和),然后将累加之后得到的值再分发给所有核,最终每个核中的 np 的值都一样。
所以加上reduce之后导致的结果就是,最终得到的每个核中的 np的值都一样,但是 np 的值显然会受核数影响。
举例说,假设 xxxx 是常数,等于10 ,则 单核运行时,np 的每个分量都是10,如果6核运行,则最终得到np 的每个分量都是 60。因为reduce之前,每个核的 np 都是10,reduce 的时候,将每个核的值累加起来,得到60。以此类推,核数越多,得到的 np 值越大。
-
可能是J-J 摩擦应力模型的问题。
上面图里面第一个是不考虑摩擦应力的结果,下面两个是考虑摩擦应力的结果。摩擦应力模型用的是 Laux,对应参考文献是:
A. Nikolopoulos, N. Nikolopoulos, N. Varveris, S. Karellas, P. Grammelis, and E.
Kakaras, Investigation of proper modeling of very dense granular flows in the
recirculation system of CFBs. Particuology, 2012, 10(6):699-709 -
@李东岳
对,摩擦应力(frictional stress)很重要。加上摩擦应力才能在双流体模型中模拟出来堆积角。 -
@hy1112006
哦,你的 refineMeshDict 里面需要一个 cellSet 来定义需要refine 的网格。这个 cellSet 也只是用来起这个作用吧。
所以,你需要在 refineMesh 之前,先 topoSet 把 cellSet 生成出来,然后 refineMesh
但是,在 decomposePar 的时候要排除对 cellSet 进行 decompose,因为我上一条回复说的原因。可以给 decomposePar 加一个选项,-noSets,这样在 decomposePar 的时候就不读取 cellSet 了,也就不会再触发你主楼遇到的错误了。 -
% 是整数求余,浮点数,不存在求余操作啊。
-
@hy1112006
先 refineMesh,后 topoSet 试试。
topoSet 生成的 cellZone 或者 cellSet 等,保存的网格ID是你refineMesh 之前的, refineMesh 之后,网格数量都变了,所以保存在 cellZone 或者 cellSet 里面的编号跟refine之后的网格都不匹配了。 -
https://www.open-mpi.org/faq/?category=openfabrics#ofa-device-error
这个链接里面的 53 应该是对应你的情况。 -
@vbcwl
我用的 openfoam.com 的版本,这个问题可能是 of9 的bug。 -
@xiaoyang-luan
RNG k-epsilon 预测不了转捩。
转捩相关的湍流模型,参考:https://turbmodels.larc.nasa.gov/
“Turbulence+Transition Models ” 部分 -
@vbcwl
哦,不好意思刚才看错了。-0.06788236085 这个数字很奇怪,没看出来是怎么来的。我试了一下也没复现出来。 -
@xiaoyang-luan
所以如果需要模拟 层流到湍流的转捩过程,需要特殊的湍流模型,一般的 k-epsilon, k-omega 模型确实是假定整个场都是湍流状态。 -
那是因为你打开来查看的那个核,恰好没有包含任何outer这个边界的面单元。 nonuniform 0() 表示一个空的list。
-
@chen_hao 在 关于修改湍流模型中非线性雷诺应力项的问题 中说:
您所指的加上雷诺应力非线性部分就是UEqn.H中turbulence->divDevRhoReff(rho, U)项的修改对吗?
不是,不需要修改这个函数,你仔细看看 ShihQuadraticKE 这个模型的继承关系,它继承的是 nonLinearEddyViscousity 类。这个类计算雷诺应力的方式不一样,除了线性项还有一个 nonLinearStress 项。你只要按照这个框架去写,把你的非线性项赋值给 nonLinearStress 就行了。
-
@chen_hao
你可以看看 src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE 这个湍流模型
这是一个非线性的湍流模型而且也只考虑了二阶非线性项,跟你要实现的那个类似。
这类模型除了修改 $\nu_t$ ,还需要把雷诺应力的非线性部分加上。 -
@vbcwl
哦,如果是这样的话,那只能通过你描述的
$$
y^+ = \frac{y u_\tau}{\nu}
$$
的方式来算 $y^+$ 了。postProcess 这种计算方式调用的是湍流模型的东西,适合用湍流模型的场景。 -
@vbcwl
请问你计算的时候,nut 的壁面边界条件设置的是什么?
postProcess 算出来的 yPlus 是调用 nut 的边界条件里面定义的 yPlus 函数算出来的。 -
@hitsc30
看你的离散格式,对流项 div(phi,U) 使用的是 linearUpwind,这个格式不适合用于 LES 仿真,数值耗散太大。建议尝试一下 LUST,或者 filteredLinear2 这类格式。
还有,如果用 PISO 求解器,p 不要用松弛。你现在用的是 0.3,这个是错的。
-
@星星星星晴
是,windows的文件名不能有冒号 -
@韬智tz
联网的目的无非是可以在线从更新源补充一些缺少的包,从而不需要手动去解决依赖问题。有一个办法是把你的系统的安装文件(一般是 ISO 文件)拿过来,做成一个本地更新源,然后就不联网也可以用类似 apt 或者 yum 这样的包管理器来补充依赖了。 -
@shrine
那还是没有安装好。 -
@shrine
好的,那应该是其他原因了。
可以先串行,再单机并行,再跨节点并行这样来测试,逐渐排查原因。
另外,可以在 controlDict 里面加上profiling { active true; cpuInfo false; memInfo false; sysInfo false; }
这样会输出来每一部分的耗时,可以对比看看差异。除了第一个 active 需要为 true,下面三个根据需要来选择true 或者 false。
-
@shrine
echo $FOAM_LIBBIN
看返回的是什么,比如我的返回的是
..../platforms/linux64GccDPInt32Opt/lib
注意看 linux64GccDPInt32Opt,后面如果是 Opt,那么就表示是优化模式,如果是 Debug,则是调试模式。
这个设置在 etc/bashrc 里面,
export WM_COMPILE_OPTION=Opt
通过这一句来指定 Opt 还是 Debug。 -
@shrine
一直在玩 ESI 的版本,从来没有遇到这种问题。一个猜想:是不是编译成 debug 模式了? -
@李东岳 在 OpenFOAM中matrix relax的bug 中说:
给定任意的松弛因子,不会影响计算结果。
这个要成立,前提条件是需要引入额外的迭代。不管是 matrix relax 还是 field relax,本质上都是把上个步骤的值跟当前步的值做一个混合之后来作为当前步的值,也就是上面提到的
$$
T^{n+1}=T^{n}+\alpha(T^{n+1}-T^n)
$$
matrix relax 虽然具体实现方式不同,但是本质不变。为避免混淆,这里使用 n+1 表示当前步,n表示上一步,n 更多的时候是内迭代步,而不是时间。
那么要想实现任意松弛因子对最终结果没有影响,唯一可能的就是迭代到一个稳定的值,也就是
$$
T^{n+1} = T^{n}
$$如果对 PISO 这样的非迭代式算法(指的是一个时间步内动量方程只求解一次)使用 matrix relax,而且设置很小的值,那么肯定会影响计算结果的啊,比如你设置U 的松弛为0.01,那么就相当于说每次都是使用99%的上一步的值加上1%的当前步算出来的值,来作为当前步的值,这样肯定就会表现为时间上基本无推进了。
OpenFOAM 里面的 PIMPLE 算法 的 matrix relax 是可以使用小于1的松弛因子的,因为 PIMPLE 算法可以理解为在时间推进步里面还有一层内迭代步,也是PIMPLE可以算比较大的 Co 数的主要原因。姑且认为从这一个时间步推进到下一个时间步都存在准确解,那么 PISO 相当于是通过一次求解动量方程,就把下一步的解算出来,PIMPLE 使用 relax,则相当于是把从上一个时间步到下一步时间步推进这个过程,使用一种类似稳态方法的方式来进行,使用小于1的松弛因子,同时增加迭代步数,这样计算会更稳定一些尤其是大时间步情况下。通过这样的迭代(多次求解动量方程和压力方程),可以逼近达到
$$
T^{n+1} = T^{n}
$$
这里的n就代表内迭代步了。从而实现虽然使用了松弛因子,但是松弛因子对结果基本无影响这个效果。另外,回到 OpenFOAM,或者说跟 OpenFOAM 类似的使用 collocated 网格的程序,都需要使用 Rhie-Chow 插值来解决棋盘压力问题。但是原始的 Rhie-Chow 方法都做不到仿真结果不受松弛因子影响。OpenFOAM 里面,可以尝试用最简单的二维 cavity 算例,用 simpleFoam 求解,一直算下去直到 initial residual 小于 1e-15,换不同的松弛因子来算,会发现虽然都能收敛,但是最终会收敛到不同的结果,虽然可能差异不会很大。
现在也有一些对 Rhie-Chow 的修正,比如
Cubero, A., & Fueyo, N. (2007). A compact momentum interpolation procedure for unsteady flows and relaxation. Numerical Heat Transfer, Part B: Fundamentals, 52(6), 507–529. https://doi.org/10.1080/10407790701563334
和
Cubero, A., Sánchez-Insa, A., & Fueyo, N. (2014). A consistent momentum interpolation method for steady and unsteady multiphase flows. Computers and Chemical Engineering, 62, 96–107. https://doi.org/10.1016/j.compchemeng.2013.12.002
-
@prometheus10
哦,明白了。你的前两种方式,等于是都在尝试修改 g 这个变量的值,但是查看 interIsoFoam求解器的代码,它其实是共用了 interFoam 里面的 UEqn.H ,在这个文件里面,重力是通过如下这段代码起作用的,fvc::reconstruct ( ( mixture.surfaceTensionForce() - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) ) * mesh.magSf() )
也就是说,真正进入到 UEqn 的是 ghf 这个变量。这个变量的定义在 gh.H 这个文件,如下:
dimensionedScalar ghRef ( mag(g.value()) > SMALL ? g & (cmptMag(g.value())/mag(g.value()))*hRef : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0) ); volScalarField gh("gh", (g & mesh.C()) - ghRef); surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
所以这就是为啥你修改 g 的值不起作用了。因为g 的值只是在开始阶段,读取进来,并据此计算了 ghf 的值,然后 ghf 进入 UEqn。要想让重力发生变化,需要 ghf进入 UEqn之前修改其值才行。或者,在计算 ghf 之前,修改 g 的值。
至于第三种体积力的方式为啥不起作用,我没看出来,主要是不清楚你说的改了 "改了interIsoFoam.C旁边的UEqn.H" 具体指的是什么,你是直接改了 interFoam 那边的 UEqn.H,然后重新编译的 interIsoFoam?
-
@prometheus10
constant 下面的 g 文件,只是一个参数设置文件,默认情况下求解器会从这个文件中读取g 的值。你在程序里面修改了g 的值,除非你把修改后的值重新写出到 g 文件,否则 g 文件就不会变了。也就是说,你修改的重力是否生效了,应该从流场结果来判断,而不是通过g文件是否发生了变化来判断。 -
@wwj
所以说这个坐标你要精心地选择,从你的图看,显然你的x方向坐标选择不合理。
假设有一层网格,x 坐标是 0.123456
那么你的 toposet 里面, x 最小值设置为 0.123455,最大值设置为0.123457 ,这样就能保证只选到一层了。 -
@李东岳
也不是抵触,infiniband 网络和 tcp 网络是共存的,他这样设置,应该是显式指定使用 tcp 网络,而没有使用 infiniband,所以速度就慢了。这个设置有的时候也是有用的,比如假设我有个工作站,没有 infiniband 也不需要,我只想单节点内的核之间通信。但是有的 mpi 它的 FI_PROVIDER 的默认值是 PSM2,这样的话如果不加参数,单节点并行也无法跑,加上 -genv FI_PROVIDER tcp 或者 -genv FI_PROVIDER shm 就可以正常跑了。
-
-genv FI_PROVIDER tcp
这一条表示你指定使用 tcp 网络通信,所以很可能你的节点间通信就没用到 infiniband。
建议先去掉 -genv FI_PROVIDER tcp ,这样mpi应该会默认选择一个可用且最快的选项。如果不行,那么参考
https://www.intel.com/content/www/us/en/develop/documentation/mpi-developer-guide-linux/top/running-applications/fabrics-control/ofi-providers-support.html
这里的说明选择一个跟你硬件匹配的 FI_PROVIDER。 -
给你一个格式参考一下
{ name samplingFace; type faceSet; action new; source boxToFace; box (2.6 0.75 0)(2.64 0.8 0.1); } { name samplingFace; type faceZoneSet; action new; source setToFaceZone; faceSet samplingFace; }
要创建 faceZone 得先创建 faceSet,上面第一部分是创建 faceSet,第二部分是根据 faceSet 创建 faceZone。
创建faceSet 的时候,你需要仔细设置 box 坐标,保证这个 box 只会框选到你需要的面(面心在这个box的face都会提取到faceSet中)。 -
@赵海盛
mesh.owner() 返回的是 内部面的 owner,跟 polyMesh 下面的owner 文件不一样,这个文件对应的是所有面的owner,对应的函数是 mesh.faceOwner() -
@wwj 在 OF如何创建内部面,该面不影响流场,但可以监测流量 中说:
改成这样:test2 { level (0 0); regions { trans {level (1 1); } } faceZone trans; }
-
你需要把 trans 单独拿出来作为一个单独的 stl文件,而不要让它是 test 文件的一部分。
-
首先,把你的监控面单独保存到一个 stl 里面
然后,snappyHexMesh 里面按照你上面图片那样,定义faceZone,比如定义为faceZone trans;
,这个trans 就是你的faceZone 名字了。functionObject里面,regionType 改为 faceZone,name 就是 trans,应该就可以了。
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
bound 函数的实现方法
bound 函数的实现方法
bound 函数的实现方法
用较小的时间步长,结果反而出问题了
IOobject::groupName("omega", alphaRhoPhi.group()) 这句话的意思是?
计算一定区域内网格数
SHM尖锐直角边界层添加
怎样才能取一个局部面作为一个新边界?
怎样把某个截面处的速度分布转移到边界上?
Gauss LUST grad(U)和Gauss LUST unlimitedGrad(U)
三维圆柱绕流升力沿管长分布
Wray-Agarwal湍流模型
mpirun并行显示运算中,但是实际log文件无输出
关于并行中的reduce函数
双流体方法颗粒堆积模拟
双流体方法颗粒堆积模拟
算例topoSet、refineMesh后,无法decomposePar
openfoam 如何求余
算例topoSet、refineMesh后,无法decomposePar
集群上OF不能跨节点并行
decomposePar后边界条件改变
RANS模型假设了整个域的湍流,那么它如何预测层流场和过渡场呢?
decomposePar后边界条件改变
RANS模型假设了整个域的湍流,那么它如何预测层流场和过渡场呢?
decomposePar后边界条件改变
关于修改湍流模型中非线性雷诺应力项的问题
关于修改湍流模型中非线性雷诺应力项的问题
yPlus在openfoam代码里面的实现
yPlus在openfoam代码里面的实现
钝体建筑扰流的大涡模拟
Field的命名
如何直接在linux下安装of
esi版本为何这么慢
esi版本为何这么慢
esi版本为何这么慢
esi版本为何这么慢
OpenFOAM中matrix relax的bug
interFoam修改随时间变化的重力
interFoam修改随时间变化的重力
OF如何创建内部面,该面不影响流场,但可以监测流量
集群上并行测试OpenFOAM,并行效率并没有比单节点提升
集群上并行测试OpenFOAM,并行效率并没有比单节点提升
OF如何创建内部面,该面不影响流场,但可以监测流量
求教forAll循环内部场问题
OF如何创建内部面,该面不影响流场,但可以监测流量
OF如何创建内部面,该面不影响流场,但可以监测流量
OF如何创建内部面,该面不影响流场,但可以监测流量