@pengdi 地面粗糙度的问题,我没有理了,这个不是我主要的研究方向,按标准的系数是没有问题了
SHUKK
帖子
-
-
@李东岳 目前跑出来的结果,有个小范围的浮动。也有可能是我在雷诺应力修正时候做了类似于DES那种,RANS跟LES的混合。
-
@李东岳 并行16线程,CFD+pytorch一共花了1个半小时左右,LES的并行36线程大概要10小时左右。CFD+pytorch的本质还是用RANS去跑,pytorch就是提供一个隐式方程,与原RANS湍流方程的时间差不多。当然硬件越好跑的越快。我们组的服务器还是从李老师你那里买的
-
libtorch跟OF耦合的问题已解决,纯粹是一个.so文件的名字大小写打错了 。目前1w步跑完结果与LES相似
-
@李东岳 上面的步骤是在编译MomentumTransportModels文件夹时候的。同时,我在simpleNNFoam中修改了我的需要连接的linearViscousStrees中的divDevTau(U, S, R),S跟R都已经在 createFields.H中创建了。
// Momentum predictor MRF.correctBoundaryVelocity(U); tmp<fvVectorMatrix> tUEqn ( fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevTau(U, S, R) // Data driven R-S + viscous component == fvModels.source(U) ); fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); fvConstraints.constrain(UEqn); if (simple.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); fvConstraints.constrain(U); }
但是我在wmake我的simpleNNFoam中出现了,数据类型不对的问题。
In file included from /home/user3/libtorch/include/c10/util/Exception.h:7, from /home/user3/libtorch/include/ATen/core/Generator.h:11, from /home/user3/libtorch/include/ATen/CPUGeneratorImpl.h:3, from /home/user3/libtorch/include/ATen/Context.h:3, from /home/user3/libtorch/include/ATen/ATen.h:7, from /home/user3/libtorch/include/torch/csrc/api/include/torch/types.h:3, from /home/user3/libtorch/include/torch/script.h:3, from /home/user3/OpenFOAM/user3-9/reynoldsNet/lnInclude/reynoldsNet.H:42, from /home/user3/OpenFOAM/user3-9/src/MomentumTransportModels/momentumTransportModels/lnInclude/linearViscousStress.H:38, from /home/user3/OpenFOAM/user3-9/src/MomentumTransportModels/momentumTransportModels/lnInclude/Stokes.H:39, from /home/user3/OpenFOAM/user3-9/src/MomentumTransportModels/momentumTransportModels/lnInclude/laminarModel.C:27, from /home/user3/OpenFOAM/user3-9/src/MomentumTransportModels/momentumTransportModels/lnInclude/laminarModel.H:193, from /home/user3/OpenFOAM/user3-9/src/MomentumTransportModels/incompressible/lnInclude/kinematicMomentumTransportModel.H:47, from simpleNNFoam.C:66: /home/user3/libtorch/include/c10/util/variant.h:1204:11: error: no matching function for call to 'Foam::data::data(const c10::detail_::impl<c10::SmallVector<c10::SymInt, 5>, at::Tensor>&)' 1204 | data(lib::forward<V>(v)), | ^~~~~~~~~~~~~~~~~~~~~~~~
我查看了 Foam::data::data(const objectRegistry& obr),要的是objectRegistry类,但是在我的reynoldsNet库中已经将libtorch输出的torch::Tensor转换成std::vector<float>然后带入输出,没有将是torch的数据直接带入OF中。这是有数据没按要求来转换吗?
-
@李东岳 李老师,我遇到一个wmake的问题。我把我训练完成的神经网络写成一个reynoldsNet库,打算在linearViscousStrees中调用,但是显示了无法识别的问题这是什么导致的?
wmakeLnIncludeAll: running wmakeLnInclude on dependent libraries: unknown option: '-I/home/user3/OpenFOAM/user3-9/reynoldsNet/lnInclude' Usage: wmakeLnInclude [OPTION] [dir] options: -update | -u update -silent | -s use 'silent' mode (do not echo command) -help | -h print the usage Link all the source files in the <dir> into <dir>/lnInclude Note The '-u' option forces an update when the lnInclude directory already exists and changes the default linking from 'ln -s' to 'ln -sf'.
-
@疏影横斜水清浅 应该是壁面粗糙度的问题,这不是我主要的研究方向
-
@疏影横斜水清浅 解决了,就是在地面时候有点保持不了
-
@李东岳 多谢李老师
-
@李东岳 李老师,我查了一下RANS的相关资料,nu是层流粘度,nut是湍流粘度,
nuEff=nu+nut,nut=Cmu*k^2/epsilon
,那是说在湍流时候,层流粘度nu是趋近于0,nuEff=nut=Cmu*k^2/epsilon
吗?不知道理解是否有问题?还请李老师指教一下。 -
@李东岳 李老师,那nu跟nut的表达式是什么,我看OF的源码中没有显示 ,nuEff=Cmu*k^2/epsilon吗?
-
各位老师,我问个简单的问题。 ,在RANS中nu_t就是一个简单的公式,但是我看在OF的incompressibleTurbulenceModel.C文件中,定义了好几个nu(),nut(),nuEff(),而在linearViscousStress.C中,只使用了nuEff(),这三个nu_t有什么区别吗?
-
@李东岳 是要像李老师学习才对,对OF的代码和耦合还在学习中
-
@李东岳 李老师,差不多是个意思,但是我们不打算找出g(i)=f(u)这个公式,我们是打算用神经网络结构来代替这个公式,也就是所谓的隐式公式,使用的神经网络是已经训练完成的,OF只需使用,不用再训练一次。现有的OF湍流模型的线性假设的公式是显式公式,也就是有明确的公式。目前的工作是要让OF程序能识别并使用神经网络来更新g(i),只用调用OF这一个程序。
-
@李东岳 李老师理解的没问题。为什么要将pytorch中的神经网络能被OF读取识别?主要是为了解决计算速度过慢的问题。可以理解为通过机器学习,获得训练完成的神经网络,也就是一个有关g(i)隐式的公式,若能被OF识别并有效输出,那我就不用每次迭代计算调用pytorch来进行计算,在两个程序中来回计算,可以直接在OF求解中同步计算,相当于直接提供显式公式的效果。对于速度慢的问题,其实不用担心,我只需调用训练完成的神经网络结构即可,相当于用一个已知的公式获得我所需要的g(i)数据。
-
@xpqiu 好的,多谢老师提供的资料。我之前打算是用caffe2来做为耦合的库,但是发现过于久远,可能有一些问题。我先了解一下libtorch的具体情况。
-
@李东岳 李老师,现在是可以在pytorch中每一步更新得到g(i),目前的问题是如何将pytorch中的神经网络能被OF读取识别,实现与OF耦合的过程。发散的是之前只在0文件下提供一个不变的标量场的方法。目前打算是将pytorch转变成c++的能识别神经网络库,g(i)的值可以通过这个库直接灌输到OF的每一个迭代步中。
-
@xpqiu 是的,应该是这个样的。但是有个问题,就是C++神经网络需要使用类似pytorch这样的环境,单加一个网络结构没办法识别计算,有那些类似pytorch这种的c++的神经网络编程语言吗?或者有什么神经网络语言是跟OF耦合的较好的呢?
-
@xpqiu 非常感谢,您的建议。我之前也是修改kepsilon湍流模型的,跟EARSM 模型一样,加了correctNonlinearStress函数。之前算发散的时候,我从g1到g10一个一个添加进去计算,发现只带入g1时候才能计算,并且得出来的效果与kOmegaSST的结果相似的,但是只要带入g2就开始发散了。后来,看之前做的过程时候发现,g(i),也就是g1到g10是一个有关gardU,k,epsilon的函数,目前无法写出其显式表达,应该用外接神经网络结构来计算更新。那为什么只带g1能算呢,因为g1*T1是等于OF的原有线性项假设,相当于没有改变其假设,只是修改了Cmu的值而已。
EARSM 模型的非线性项中,我看它的非线性项写法与我的类似,但是它的bate(i)是有明确的显式公式的,是能随着求解器迭代来计算更新的,而我的是隐式公式,目前没有明确公式。volScalarField beta1 = -N*(2.0*sqr(N) - 7.0*IIW) / Q; volScalarField beta3 = -12.0 * IV / (N * Q); volScalarField beta4 = -2.0 * (sqr(N) - 2.0*IIW) / Q; volScalarField beta6 = -6.0 * N / Q; volScalarField beta9 = 6.0 / Q; volScalarField Cmu = - 0.5 * (beta1 + IIW * beta6); this->nut_ = Cmu * this->k_ * tau; this->nut_.correctBoundaryConditions(); this->nonlinearStress_ = this->k_ * symm( beta3 * ( (W & W) - (1.0/3.0) * IIW * I ) + beta4 * ( (S & W) - (W & S) ) + beta6 * ( (S & W & W) + (W & W & S) - IIW * S - (2.0/3.0) * IV * I) + beta9 * ( (W & S & W & W) - (W & W & S & W) )
现在我要做的内容就让我的g(i)跟它bate(i)一样跟随迭代,为了达到我这个效果,不能单修改湍流模型,还要外加上我的c++神经网络来更新g(i)。目前就是编译的一些内容和步骤有点疑惑,需要请教。
-
@李东岳 好的,李老师。我先按您的方法试试。
-
@李东岳 李老师,是说可以模仿在constant/turbulenceProperties中修改湍流模型的Cmu这类常数数值的步骤一样吗?把g(i)的c++神经网络写到turbulenceProperties进去,让求解器读取,不用修改linearViscousStrees和IncompressibleTurbulenceModel。
-
@李东岳 多谢李老师对我目前所做的项目有个极高的评价,由于我已经是研三了,这是我课题的最后一部分,“联合培养”项目开展时间有点晚了,多谢李老师的好意。对于我目前所做的东西,感觉李老师理解稍有偏差。g(i),也就是g1到g10是一个有关gardU,k,epsilon的函数,目前无法写出其显式表达,因此打算用神经网络结构来表示,所以按之前步骤只放到0文件夹下面就是一个不变量,迭代一次就应该更新了,每一次OF求解迭代出来的g(i)应该是不一样的。之前我是想跟python耦合的,即求出一步就用学习完成的神经网络更新一次g(i),但是对于高网络数量的模型计算就十分缓慢,因此我才打算将在python的神经网络结构转换为c++能识别的神经网络结构,并重新编写linearViscousStrees和IncompressibleTurbulenceModel来进行OF耦合。目前,修改的代码写的差不多了,就是不知linearViscousStrees和IncompressibleTurbulenceModel的有效编译步骤该如何进行,还需请教一下李老师。
-
目前,通过这些步骤我有几个问题。
1、修改的linearViscousStrees和IncompressibleTurbulenceModel.C和.H文件,我也需要像simpleFoam求解器一样修改名称和USER位置来编译,避免与原有的OF内容冲突吗?或是有什么特殊的编译方法?
2、我看之前有关的帖子中,说只要修改linearViscousStrees.C和.H文件来达到修正效果,那这修正IncompressibleTurbulenceModel.C和.H文件的意义是什么呢?
有做过相关内容的老师们,解答一下吗? -
目前,是重新编写了simpleFoam,可以是在计算迭代中同时生成Sij跟Rij数据,为后面的linearViscousStrees提供相应数据,同时修改了动量方程的有关雷诺应力项。
tmp<fvVectorMatrix> tUEqn ( fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevRhoReff(U, S, R) == fvOptions(U)
我参照了之前的有关文献,同时修改了incompressible文件下的,incompressibleTurbulenceModel.H。
//- Return the source term for the momentum equation using NN surrogate virtual tmp<fvVectorMatrix> divDevReff ( volVectorField& U, volTensorField& S, volTensorField& R ) { }
incompressible/IncompressibleTurbulenceModel文件夹下的,IncompressibleTurbulenceModel.C和.H文件,添加了以下内容
template<class TransportModel> Foam::tmp<Foam::fvVectorMatrix> Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff ( volVectorField& U, volTensorField& S, volTensorField& R ) const { return divDevRhoReff(U, S, R); }
//- Return the source term for the momentum equation virtual tmp<fvVectorMatrix> divDevReff ( volVectorField& U, volTensorField& S, volTensorField& R ) const;
并修改了linearViscousStrees.C和.H文件,添加了以下内容
template<class BasicTurbulenceModel> Foam::tmp<Foam::fvVectorMatrix> Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff ( volVectorField& U, volTensorField& S, volTensorField& R ) const { NN结构:a_dd } return ( - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U) - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))) - this->a_dd) );
a_dd就是我要加入的非线性项,里面的内容主要是链接了reynoldsNet神经网络。
-
各位有经验的老师们,我想问一下加入修正雷诺应力项的修正湍流模型的问题。目前,我已经尝试了模仿ShihQuadraticKE湍流模型写法,重新写了一个新的湍流模型,具体修改的地方如下:
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // void kEpsilonTBNN::correctNut() { correctNonlinearStress(fvc::grad(U_)); } void kEpsilonTBNN::correctNonlinearStress(const volTensorField& gradU) { timeScale_=k_/epsilon_; // Linear (nut) nut_ = g1_*k_*timeScale_; nut_.correctBoundaryConditions(); // nolinear(tau_NL) volSymmTensorField S(timeScale_*symm(gradU)); volTensorField W(timeScale_*skew(gradU)); nonlinearStress_ = 2*k_ *( g2_ * twoSymm(S&W) + g3_ * dev(innerSqr(S)) + g4_ * dev(symm(W&W)) + g5_ * twoSymm(W&innerSqr(S)) + g6_ * dev(twoSymm(W&W&S)) + g7_ * twoSymm(W&S&W&W) + g8_ * twoSymm(S&W&innerSqr(S)) + g9_ * dev(twoSymm(innerSqr(S)&W&W)) + g10_ * twoSymm(W&innerSqr(S)&W&W) ); }
这是我主要的修正项,即Aij
其中,g1到g10是我通过机器学习出来的数据,
g(i)跟T(i)都是有关Sij跟Ri的公式。目前,按照这种方法计算是发散的,其主要的原因是g(i)跟T(i)是随每次迭代求解出的U,k,epsilon发生变化的,但是只通过读取0文件下的g(i)场只是一个不变的标量,从而导致求解发散,因此需要将学习出来的神经网络结构耦合到OpenFOAM的环境中。因此我打算修改simpleFoam,linearViscousStrees和新加的reynoldsNet神经网络结构。
-
@李东岳 谢谢李老师,我再回去看看是不是我公式的问题
-
@李东岳 李老师,我现在检查了OF有关非线性雷诺应力,也就是各向异性的代码和公式。
看您在NS无痛苦笔记中写的,是在线性项后面加上非线性项的,但是在OF代码中,为啥是减上非线性项呢?
-
void kEpsilonNN::correctNonlinearStress(const volTensorField& gradU) { timeScale_=k_/epsilon_; // Linear (nut) nut_ = -g1_*k_*timeScale_; nut_.correctBoundaryConditions(); // nolinear(tau_NL) volSymmTensorField S(timeScale_*symm(gradU)); volTensorField W(timeScale_*skew(gradU)); nonlinearStress_ =2*k_*(g2_ * twoSymm(S&W)+ g3_ * dev(innerSqr(S))+ g4_ * dev(symm(W&W)) + g5_ * twoSymm(W&innerSqr(S))+ g6_ * dev(twoSymm(W&W&S))+ g7_ * twoSymm(W&S&W&W) + g8_ * twoSymm(S&W&innerSqr(S)) + g9_ * dev(twoSymm(innerSqr(S)&W&W)) + g10_ * twoSymm(W&innerSqr(S)&W&W)); }
目前是在雷诺应力中加了各向异性项nonlinearStress,是按照LienCubicKE.C来修改的,其中g1到g10是一个标量场,在0文件夹下面读取。现在是湍流模型编译能通过,但是在调用simple求解器时候计算出现问题。
Time = 2 smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 0.0279014, No Iterations 2 smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 0.0209915, No Iterations 2 smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 0.0274975, No Iterations 2 GAMG: Solving for p, Initial residual = 1, Final residual = 0.0920741, No Iterations 4 time step continuity errors : sum local = 7.40969e+91, global = -1.79495e+89, cumulative = -1.79495e+89 smoothSolver: Solving for epsilon, Initial residual = 1, Final residual = 0.0321531, No Iterations 3 bounding epsilon, min: -1.33631e+104 max: 1.64016e+103 average: -1.01414e+99 smoothSolver: Solving for k, Initial residual = 1, Final residual = 0.072621, No Iterations 2 bounding k, min: -1.34214e+100 max: 1.0091e+99 average: -9.64939e+94 [1] #0 Foam::error::printStack(Foam::Ostream&)[6] #0 Foam::error::printStack(Foam::Ostream&)[7] #0 Foam::error::printStack(Foam::Ostream&)[8] #0 Foam::error::printStack(Foam::Ostream&)[9] #0 Foam::error::printStack(Foam::Ostream&)[12] #0 Foam::error::printStack(Foam::Ostream&)[13] #0 Foam::error::printStack(Foam::Ostream&)[15] #0 Foam::error::printStack(Foam::Ostream&)[2] #0 Foam::error::printStack(Foam::Ostream&)[5] #0 Foam::error::printStack(Foam::Ostream&)[3] #0 Foam::error::printStack(Foam::Ostream&)[11] #0 Foam::error::printStack(Foam::Ostream&)[10] [14] #0#0 [4] #0 Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)[0] #0 Foam::error::printStack(Foam::Ostream&) at ??:? [1] at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? [2] #1 Foam::sigFpe::sigHandler(int)[10] #1 Foam::sigFpe::sigHandler(int) at ??:? [15] #1 Foam::sigFpe::sigHandler(int) at ??:? at ??:? at ??:? at ??:? [5] #1 Foam::sigFpe::sigHandler(int)[13] #1 Foam::sigFpe::sigHandler(int)[14] #1 Foam::sigFpe::sigHandler(int) at at [6] #1 at ??:? ??:? Foam::sigFpe::sigHandler(int) at ??:? at ??:? ??:? at ??:? at ??:? at ??:? [7] #1 [8] #1 Foam::sigFpe::sigHandler(int)Foam::sigFpe::sigHandler(int)[12] #1 Foam::sigFpe::sigHandler(int)[0] #1 Foam::sigFpe::sigHandler(int)[4] #1 Foam::sigFpe::sigHandler(int)[9] #1 Foam::sigFpe::sigHandler(int)[3] #1 Foam::sigFpe::sigHandler(int)[11] #1 Foam::sigFpe::sigHandler(int) at ??:? [5] #2 ? at ??:? [10] #2 ? at ??:?
对此,我为了验证是不是我代码写错的问题,我用了一个简单的二维方管来验证,发现可以计算,效果还可以。我也把我修改的地方一一注释找问题,也没找到问题所在。有人知道是什么问题吗?
-
@李东岳 李老师,我看了应该是一样的,第一个是v2206,第二个是of9
const volScalarField::Internal divU ( fvc::div(fvc::absolute(this->phi(), U))().v() ); tmp<volTensorField> tgradU = fvc::grad(U); const volScalarField::Internal GbyNu ( this->type() + ":GbyNu", tgradU().v() && dev(twoSymm(tgradU().v())) ); const volScalarField::Internal G(this->GName(), nut()*GbyNu); tgradU.clear()
eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct(); volScalarField::Internal divU ( fvc::div(fvc::absolute(this->phi(), U))() ); tmp<volTensorField> tgradU = fvc::grad(U); volScalarField::Internal G ( this->GName(), nut()*(dev(twoSymm(tgradU().v())) && tgradU().v()) ); tgradU.clear();
-
@李东岳 李老师,我试了按您那个CFD张量公式的介绍来写,能编译成功,但是按LienCubicKE.C里面的epsilon方程和k方程来写,运算时候直接发散了。于是,我按标准的kEpsilon的公式来写时,发现v2206和OF9的方程不一样,这是有什么说法吗?第一个是v2206的,第二个是of9的,两者差异在GbyNu的位置,一个跟Cmu,一个跟epsilon。我看蛮多文章里面的公式都跟of9一样的。
tmp<fvScalarMatrix> epsEqn ( fvm::ddt(alpha, rho, epsilon_) + fvm::div(alphaRhoPhi, epsilon_) - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_) == C1_*alpha()*rho()*GbyNu*Cmu_*k_() //Cmu? - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_) - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_) + epsilonSource() + fvOptions(alpha, rho, epsilon_) );
tmp<fvScalarMatrix> epsEqn ( fvm::ddt(alpha, rho, epsilon_) + fvm::div(alphaRhoPhi, epsilon_) - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_) == C1_*alpha()*rho()*G*epsilon_()/k_() //epsilon? - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_) - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_) + epsilonSource() + fvModels.source(alpha, rho, epsilon_) );
-
@李东岳 李老师,我是照着LienCubicKE.C的模板来加雷诺应力非线性项的,我看它的代码里面也是一样的写法,但是为什么就它的编译成功了呢?
void LienCubicKE::correctNonlinearStress(const volTensorField& gradU) { volSymmTensorField S(symm(gradU)); volTensorField W(skew(gradU)); volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(S)); volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(W)); volScalarField Cmu((2.0/3.0)/(Cmu1_ + sBar + Cmu2_*wBar)); volScalarField fMu(this->fMu()); nut_ = Cmu*fMu*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); nonlinearStress_ = fMu*k_ *( // Quadratic terms sqr(k_/epsilon_)/(Cbeta_ + pow3(sBar)) *( Cbeta1_*dev(innerSqr(S)) + Cbeta2_*twoSymm(S&W) + Cbeta3_*dev(symm(W&W)) ) // Cubic terms - pow3(Cmu*k_/epsilon_) *( (Cgamma1_*magSqr(S) - Cgamma2_*magSqr(W))*S + Cgamma4_*twoSymm((innerSqr(S)&W)) ) ); }
目前参考其他贴子[链接文本](https://www.cfd-china.com/topic/4128/修改湍流模型-定义表达式的问题?_=1705911573991) 里面的改,
volTensorField T2((S&W)-(W&S)); volSymmTensorField T3(dev(innerSqr(S))); volTensorField T4(dev(W & W)); volTensorField T5((W & innerSqr(S)) - (innerSqr(S) & W)); volTensorField T6((W & W & S) - (dev2(S & W & W))); volTensorField T7((W & S & W & W) - (W & W & S & W)); volTensorField T8((S & W & innerSqr(S)) - (innerSqr(S) & W & S)); volTensorField T9((W & W & innerSqr(S)) + (dev2(innerSqr(S) & W & W))); volTensorField T10((W & innerSqr(S) & W & W) - (W & W & innerSqr(S) & W)); volTensorField::Internal nonlinearStress_(2*k_*(g2_ * T2+ g3_ * T3+ g4_ * T4+ g5_ * T5+ g6_ * T6+ g7_ * T7+ g8_ * T8+ g9_ * T9+ g10_ * T10 ));
是编译成功了,但是有个问题,这个内部场的写法是有问题吗?我求的是全流场的
-
-
这是我在雷诺应力中加入非线性的部分 ,目前是用复制的kEpsilon.C文件中进行修改,这里面的g1到g10都是标量,是一个系数。S跟W都是向量
void kEpsilonNN::correctNonlinearStress(const volTensorField& gradU) { timeScale_=k_/epsilon_; // Linear (nut) nut_ = -g1_*k_*timeScale_; nut_.correctBoundaryConditions(); // Quadratic (tau_NL) volSymmTensorField S(timeScale_*symm(gradU)); volTensorField W(timeScale_*skew(gradU)); nonlinearStress_ = 2 * k_ *( g2_ * ((S&W)-(W&S)) + g3_ * dev(innerSqr(S)) + g4_ * dev(W&W) + g5_ * ((W&innerSqr(S))-(innerSqr(S)&W)) + g6_ * ((W&W&S)-(dev2(S&W&W))) + g7_ * ((W&S&W&W)-(W&W&S&W)) + g8_ * ((S&W&innerSqr(S))-(innerSqr(S)&W&S)) + g9_ * ((W&W&innerSqr(S))+(dev2(innerSqr(S)&W&W))) + g10_ * ((W&innerSqr(S)&W&W)-(W&W&innerSqr(S)&W)) ); }
但是在编译的时候出现了这个问题。
error: no match for 'operator=' (operand types are 'Foam::volSymmTensorField' {aka 'Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh>'} and 'Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >') 81 | ); | ^
这是在我nonlinearStress_最后一个括号出现的。这是说我的nonlinearStress_没申明成向量还是啥呢?有人知道的吗?
还有人知道
, 是表示为 还是其他的呢? -
目前打算从稍微简单的kEpsilon或kOmage湍流模型进行修改
-
各位OF大佬,你们好!目前,已经通过神经网络对雷诺应力进行了学习,这是我用的雷诺应力公式,
其中,g(n)是我通过神经网络学习出来的,现在我想把雷诺应力的结果与OpenFOAM原有的湍流模型结合,利用RANS求解器来进一步变成U和p。由于是小白,就想问一下假如要这样做的实现湍流模型代码修改方向和编译的步骤。
1、在原有的雷诺应力公式加上各向异性项bij?
或是2、重新编写新的雷诺应力公式?
-
这是我的.C文件```
/---------------------------------------------------------------------------\\ / F ield OpenFOAM: The Open Source CFD Toolbox \ / O peration \ / A nd www.openfoam.com \/ M anipulation
Copyright (C) 2011-2017 OpenFOAM Foundation Copyright (C) 2019 OpenCFD Ltd.
License
This file is part of OpenFOAM.OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
*---------------------------------------------------------------------------*/
#include "kEpsilonNNQuadraticTrain.H"
#include "bound.H"
#include "wallFvPatch.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kEpsilonNNQuadraticTrain, 0);
addToRunTimeSelectionTable(RASModel, kEpsilonNNQuadraticTrain, dictionary);// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void kEpsilonNNQuadraticTrain::correctNut()
{
correctNonlinearStress(fvc::grad(U_));
}void kEpsilonNNQuadraticTrain::correctNonlinearStress(const volTensorField& gradU)
{
timeScale_=k_/epsilon_;// Linear (nut) nut_ = -g1_*k_*timeScale_; nut_.correctBoundaryConditions(); // Quadratic (tau_NL) volSymmTensorField S(timeScale_*symm(gradU)); volTensorField W(timeScale_*skew(gradU)); nonlinearStress_ = 2*k_ *( g2_ * twoSymm(S&W) + g3_ * dev(innerSqr(S)) + g4_ * dev(symm(W&W)) );
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
kEpsilonNNQuadraticTrain::kEpsilonNNQuadraticTrain
(
const geometricOneField& alpha,
const geometricOneField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
nonlinearEddyViscosityincompressible::RASModel
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
),Ceps1_ ( dimensioned<scalar>::lookupOrAddToDict ( "Ceps1", coeffDict_, 1.44 ) ), Ceps2_ ( dimensioned<scalar>::lookupOrAddToDict ( "Ceps2", coeffDict_, 1.92 ) ), sigmak_ ( dimensioned<scalar>::lookupOrAddToDict ( "sigmak", coeffDict_, 1.0 ) ), sigmaEps_ ( dimensioned<scalar>::lookupOrAddToDict ( "sigmaEps", coeffDict_, 1.3 ) ), k_ ( IOobject ( IOobject::groupName("k", alphaRhoPhi.group()), runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), epsilon_ ( IOobject ( IOobject::groupName("epsilon", alphaRhoPhi.group()), runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), g1_ ( IOobject ( "g1", runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), g2_ ( IOobject ( "g2", runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), g3_ ( IOobject ( "g3", runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), g4_ ( IOobject ( "g4", runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ ), timeScale_ ( IOobject ( "timeScale", runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh_, dimensionedScalar("timeScale", dimTime, scalar(0.0)) )
{
bound(k_, kMin_);
bound(epsilon_, epsilonMin_);if (type == typeName) { printCoeffs(type); }
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool kEpsilonNNQuadraticTrain::read()
{
if (nonlinearEddyViscosityincompressible::RASModel::read())
{
Ceps1_.readIfPresent(coeffDict());
Ceps2_.readIfPresent(coeffDict());
sigmak_.readIfPresent(coeffDict());
sigmaEps_.readIfPresent(coeffDict());return true; } return false;
}
void kEpsilonNNQuadraticTrain::correct()
{
if (!turbulence_)
{
return;
}nonlinearEddyViscosity<incompressible::RASModel>::correct(); tmp<volTensorField> tgradU = fvc::grad(U_); const volTensorField& gradU = tgradU(); volScalarField G ( GName(), (nut_*twoSymm(gradU) - nonlinearStress_) && gradU ); // Update epsilon and G at the wall epsilon_.boundaryFieldRef().updateCoeffs(); // Dissipation equation tmp<fvScalarMatrix> epsEqn ( fvm::ddt(epsilon_) + fvm::div(phi_, epsilon_) - fvm::laplacian(DepsilonEff(), epsilon_) == Ceps1_*G*epsilon_/k_ - fvm::Sp(Ceps2_*epsilon_/k_, epsilon_) ); epsEqn.ref().relax(); epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef()); solve(epsEqn); bound(epsilon_, epsilonMin_); // Turbulent kinetic energy equation tmp<fvScalarMatrix> kEqn ( fvm::ddt(k_) + fvm::div(phi_, k_) - fvm::laplacian(DkEff(), k_) == G - fvm::Sp(epsilon_/k_, k_) ); kEqn.ref().relax(); solve(kEqn); bound(k_, kMin_); // Re-calculate viscosity and non-linear stress correctNonlinearStress(gradU);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam// ************************************************************************* //
-
我主要编译的是一个基于kEpsilon模型修改的
-
请教各位老师,我试了我能找到的几乎所有编译新湍流模型的流程,到最后都是这个问题,有人知道怎么处理的吗?
-
@李东岳 李老师,按照你的提示,我尝试先simple -postProcess -func R,先在o中生成了R,然后再运行simple -postProcess -func probes,然后就出现这个问题了。
-
@李东岳 李老师,还是在system中写一个sampleDict然后指定点可以输出所需要的吗?
-
@李东岳 有点没明白,是要跟LES一样在control写出场的平均代码才能在probes输出吗?还是先postProcess -func R,再probes?
-
@李东岳 多谢李老师!刚刚去试了试,确实可以。那雷诺应力R呢?我只知道postPrcoess -func R 来获取全部的,LES的也是。
-
@李东岳 李老师,probes不是只能提取基本物理量吗?
-
大家好!对于OpenFOAM的基本物理量,例如U,P,都知道可以通过探针probes来进行提取,在postProcessing的probes中找到。那么非衍生物理量,例如K,Epsilon,Omega,雷诺应力等等,无法利用probes来直接提取。虽然在RANS的时间步文件中,我们直接得到计算得到的k,epsilon和omega中,但是这是全部网格的,从中提取出来的我们要的,需要知道特定点对应的网格位置。对于提取特定点的非基本物理量的方法,要去问“伟大”的chatGPT了吗?
-
@李东岳 好的,谢谢李老师!y2017kOmegaSST.zip
-
@李东岳 李老师,我最近尝试复现Yang2017的kOmegaSST,按照文献的设置和参数设置,k都无法保持稳定。不知道老师您是否复现出来了。
-
@李东岳 李老师,我使用的是指数率的风速,参考风速11m/s,高度0.4m。之前用LES算过yPlus是符合要求的,第一次用RANS算。
-
@李东岳 李老师,一开始设置y+=1,计算过程中不是应该是在1的小范围内来回波动的吗?我这个已经冲1一直到几百了,应该是有问题了吧
-
大家好!我使用的是kOmegaSST来对钝体绕流进行计算。对于kOmegaSST的yPlus一般要求在1左右,通过计算yPlus会随着时间不断地扩大,有人遇到过这个问题吗?
# y+ () # Time patch min max average 100 lowerwall 2.281880e-01 6.964871e+03 2.434004e+02 100 building 1.573129e-02 5.619403e+00 9.878271e-01 200 lowerwall 8.873161e-01 2.276240e+04 1.195805e+03 200 building 1.231834e-01 5.265857e+00 1.497548e+00 300 lowerwall 9.941788e+00 1.198720e+05 6.866907e+03 300 building 3.843790e-01 2.573107e+01 9.721960e+00 400 lowerwall 1.628298e+01 1.685941e+06 5.717353e+04 400 building 1.956629e+00 2.150759e+02 9.664974e+01 500 lowerwall 2.190860e+02 2.622048e+07 7.852843e+05 500 building 1.544316e+02 4.432274e+03 1.716193e+03
-
@李东岳李东岳老师,我看编写应该是没有问题的
ktype codedFixedValue; value uniform 0.375; //default value name kinlet; //name of new BC type code #{ const fvPatch& boundaryPatch = this->patch(); scalarField& vf = *this; forAll(vf, i) { scalar z = boundaryPatch.Cf()[i].y(); scalar D1 = -1.061; scalar D2 = 5.744; scalar alpha = 0.25; vf[i] = sqrt(D1*pow(z*400, alpha) + D2);
omega
type codedFixedValue; value uniform 0.1; //default value name omegainlet; //name of new BC type code #{ const fvPatch& boundaryPatch = this->patch(); scalarField& vf = *this; forAll(vf, i) { scalar z = boundaryPatch.Cf()[i].y(); scalar ur = 11.0; scalar alpha = 0.25; scalar Cmu = 0.06; scalar zr = 160; scalar u = ur*pow(z*400/zr, alpha); vf[i] = (alpha/sqrt(Cmu))*u/z;
u
type codedFixedValue; value uniform (0 0 0); //default value name uinlet; //name of new BC type code #{ const fvPatch& boundaryPatch = this->patch(); vectorField& vf = *this; forAll(vf, i) { scalar z = boundaryPatch.Cf()[i].y(); //scalar uStar = 0.511; //scalar z0 = 2.25e-4; //scalar kappa = 0.42; //vf[i].x() = uStar/kappa*log((z + z0)/z0); //vf[i].y() = 0.0; //vf[i].z() = 0.0; scalar ur = 11.0; scalar alpha = 0.25; scalar zr = 160; vf[i].x() = ur*pow(z*400/zr, alpha); vf[i].y() = 0.0; vf[i].z() = 0.0;
中性大气环境湍流动能的自保持 | 附有算例下载
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
中性大气环境湍流动能的自保持 | 附有算例下载
中性大气环境湍流动能的自保持 | 附有算例下载
OF中有效粘度nu_t的代码问题
OF中有效粘度nu_t的代码问题
OF中有效粘度nu_t的代码问题
OF中有效粘度nu_t的代码问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
原有的湍流模型加上非线性项雷诺应力的问题
场计算问题
场计算问题
场计算问题
湍流模型编译的问题
湍流模型编译的问题
湍流模型编译的问题
湍流模型编译的问题
湍流模型编译的问题
神经网络和OpenFOAM组合的新湍流模型问题讨论
神经网络和OpenFOAM组合的新湍流模型问题讨论
OpenFOAM中编译新的湍流模型报错
OpenFOAM中编译新的湍流模型报错
OpenFOAM中编译新的湍流模型报错
如何提取特定点的非基本物理量
如何提取特定点的非基本物理量
如何提取特定点的非基本物理量
如何提取特定点的非基本物理量
如何提取特定点的非基本物理量
如何提取特定点的非基本物理量
中性大气环境湍流动能的自保持 | 附有算例下载
中性大气环境湍流动能的自保持 | 附有算例下载
关于yPlus不随计算过程保持稳定的问题
关于yPlus不随计算过程保持稳定的问题
关于yPlus不随计算过程保持稳定的问题
openfoam运行报错