|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
首先得声明的是此文非我写（连翻译，带添油加醋自HaKan Nilson的 How to implement your own turbulence model），此code非我撰，原创性不是特别高，但是想来湍流是流体力学领域一个重要的研究方向。工程湍流也是一大热门。不少同仁一定花了大量时间在做这方面的工作，特别是不少从事工程湍流计算的朋友，一定苦于商业软件里湍流模型太少，想要实现一些修正模型又限制重重无从下手，想做些深入的工作比较难。再者本人最近一直在琢磨湍流模型－－真是个有趣却实在是烦人的问题。感觉总结一下希望能抛烂泥引出块金砖。希望拍砖者成群，让我收个盆满钵（那啥～～汗阿！语文不好，希望不要雷到大家）。
言规正传，有关OpenFOAM里C＋＋编程基础，请各位参看一下Junwei同学的博文。OpenFOAM 1.5的湍流模型位于$FOAM_INSTALL_DIR/src/turbulenceModels，与前面的版本不一样在OF1.5里大涡模拟（LES）和RANS（雷诺时均方法）都并到这个位置了。这里就讲讲RANS里的内容吧，RANS目录下有可压和不可压两块，以下内容均为基于不可压流的。
Part A 入门篇
对于初学者建议先将一个最接近你模型的湍流模型拷贝出来加以修改。在这里我以V2F模型为例子介绍一下简单的过程。V2F模型是近年来涌现的一种新的高级湍流模型（新就是比老的要年轻，高级就是复杂么。。。雷到了吧？），有较为精确的的处理自由剪切湍流的能力，尤其是近壁面湍流的计算上有着良好的表现。目前相关工作的主要研究人员就是那个大名鼎鼎的Durbin.。已经有不少商业软件如starCD，fluent都在最新的版本里添加了该功能，不过fluent里要使用V2F模型需要额外的license，有关V2F模型，可以仔细看看parneix的这篇参考文献。整体看来这个V2F模型里也用到湍动能k和耗散率Epsilon的方程，当然还有其他的方程。所以我们可以基本上拷贝kEpsilon模型里的两个文件到你的文件夹(~/V2F)下改名叫DurbinV2F.C 和DurbinV2F.H，并建立一个Make文件夹，在Make文件夹里创建两个文件一个是files 一个是optionn内容如下：
option:
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels/RAS/incompressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/transportModels
LIB_LIBS = \
-lfiniteVolume \
-lmeshTools
files
V2F/DurbinV2F.C //V2F是存放code的位置，DurbinV2F.C是该模型主程序
LIB = $(FOAM_USER_LIBBIN)/libaddedIncompressibleRASModels //生成的动态库文件存在$(FOAM_USER_LIBBIN)目录下 libaddedIncompressibleRASModels.so文件
接下来要做的是打开DurbinV2F.C 和DurbinV2F.H 把所有kEpsilon 改成DurbinV2F，这样就可以得到新的class V2F,当然里面的一切都还是kEpsilon
你可以简单的添加小的改变 比如在第89行加入：
Info << "Defining my own kEpsilon model" <<endl;
然后就可以编译了
在～/V2F文件夹下执行
wmake libso
这样一个新湍流模型的动态库 就“哦可”拉！
下面你可以试试看这个冒牌的V2F模型，该咋整呢？？ 接下去看...
请到$FOAM_RUN/tutorials/simpleFoam/pitzDaily/constant下
修改RASPorperties文件
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
RASModel DurbinV2F; //kEpsilon 把原来的kEpsilon换成DurbinV2F 下同 ;
turbulence on;
printCoeffs on;
laminarCoeffs
{
}
DurbinV2FCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
alphaEps 0.76923;
}
接着在$FOAM_RUN/tutorials/simpleFoam/system/文件夹下打开controlDict字典文件加如
libs (libaddedIncompressibleRASModels.so);
告诉OpenFoam调用你的新动态库，然后在$FOAM_RUN/tutorials/simpleFoam/pitzDaily 下先运行 BlockMesh 再运行 simpleFoam好了,我们的“冒牌V2F”可以计算了。
Part B 进阶篇
A note on new turbulence models--这一段 照猫画虎了，
1）OpenFoam里的湍流模型是虚基类RASModels的子类。
2）一般我们只能使用RASModels里定义过的成员函数，如果你需要其他成员函数的话，你需要到RASModels类里去添加，你需要修改src/turbulecenmodel/RAS/incompressible/RASModels的相关文件
3）你可以按照OpenFOAM里的文件结构来安排你的新模型。但也未必如此。
下面我们就开始v2f之旅吧。。
这里v2f模型是的源码来自于cfd－online－v2f-openfoam，当然是1.4.1的1.5的你可能下下来还需要小做调整。当然我已经弄好了。在最后么提供下载。当然这个模型也不是原始的V2F模型，是一个修正模型所以相关的修改需要参看一下这个文档，也可在文末附件里下载。
在这里先简单的介绍一下V2F模型，当我对这个模型也不是很熟悉，讲的不太准确，一切还是以上面这个参考文档为准。V2F模型是在标准kEpsilon模型基础上增加了一个虚拟雷诺主应力v2的控制方程（参考文献中的eqn12），以及f的方程 （参考文献中的eqn13）。在此基础上重新定义了涡黏性系数（参考文献中的eqn11）。并且利用Kolmogorov的湍流时间尺度T（参考文献中的eqn14）和空间尺度L（参考文献中的eqn15）对控制方程的部分进行了限制，以保证在在近壁面区（这里是指Y+<30的区域）不会出现奇异。当然相关的系数都会发生改变，这些也都在文献里具体给出。
那么首先我们在DurbinV2F.H文件里进行相关，成员变量和成员方法
/*---------------------------------------------------------------------------*\
Class DurbinV2F Declaration
\*---------------------------------------------------------------------------*/
class DurbinV2F
:
public RASModel
{
// Private data
dimensionedScalar Cmu;
dimensionedScalar CmuKE;
dimensionedScalar Ceps10;
dimensionedScalar Ceps11;
dimensionedScalar Ceps2;
dimensionedScalar C1;
dimensionedScalar C2;
dimensionedScalar CL;
dimensionedScalar CEta;
dimensionedScalar oneOnSigmaK;
dimensionedScalar oneOnSigmaEps;
dimensionedScalar yStarLim;
dimensionedScalar f0_;
wallDist yw_;
volScalarField k_;
volScalarField epsilon_;
volScalarField v2_; // 定义标量场v2_
volScalarField f_; //定义标量场f_
volScalarField nut_;
// Private member functions
tmp<volScalarField> T() const;//定义时间尺度的标量场模板T
tmp<volScalarField> L() const;//定义空间尺度的标量场模板L
public:
//- Runtime type information
TypeName("DurbinV2F");
// Constructors
.....
tmp<volScalarField> v2() const
{
return v2_;
}
tmp<volScalarField> f() const
{
return f_;
}
//在公共部分定义标量场模板v2 和L 用以返回这两个量
};
下面进入主程序DurbinV2F.C
首先我们需要增加对T和L的私有成员方法
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
tmp<volScalarField> DurbinV2F::T() const
{
volScalarField yStar_=pow(CmuKE,0.25)*sqrt(k_)*yw_/nu();
return max(k_/(epsilon_ + epsilonSmall_),
pos(yStarLim-yStar_)*6.0*sqrt(nu()/(epsilon_ + epsilonSmall_)));
}
tmp<volScalarField> DurbinV2F: () const
{
volScalarField yStar_=pow(CmuKE,0.25)*sqrt(k_)*yw_/nu();
return
CL*max(pow(k_,1.5)/(epsilon_ + epsilonSmall_),
pos(yStarLim-yStar_)
*CEta*pow(pow(nu(),3.0)/(epsilon_ + epsilonSmall_),0.25));
}
参考文献中的eqn14，15有定义
在构造函数里初始化时，这些系数都可以RASModelDict字典文件读取,也可以在这初始化
// from components
DurbinV2F: urbinV2F
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& lamTransportModel
)
:
RASModel(typeName, U, phi, lamTransportModel),
Cmu
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
coeffDict_,
0.09
)
),
.....
v2_
(
IOobject
(
"v2",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
f_
(
IOobject
(
"f",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
// Calculate viscosity (with Davidson correction - 2003)
nut_(min(CmuKE*sqr(k_)/(epsilon_ + epsilonSmall_),
Cmu*v2_*T()))
接着就要在void DurbinV2F::correct()这个方法里 进行相关修改
void DurbinV2F::correct()
{
......
volScalarField S2 = 2*magSqr(symm(fvc::grad(U_)));
volScalarField G = nut_*S2;
volScalarField T_ = T();
volScalarField Ceps1 = Ceps10*(scalar(1.0)+Ceps11*min(sqrt(k_/v2_),scalar(10.0)));//改变 C1定义
// Dissipation rate （Epsilon）equation
# include "epsilonV2FWallI.H"
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
Ceps1*G/T_
- fvm::Sp(Ceps2/T_, epsilon_) //用空间尺度T_对源项的限制
);
epsEqn().relax();
# include "wallDissipationV2FI.H"
solve(epsEqn);
bound(epsilon_, epsilon0_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp(1.0/T_, k_) //用空间尺度T_对源项的限制
// G - fvm::Sp(epsilon_/k_, k_)
);
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
接下来就是f方程 以及v2方程
// f equation
volScalarField L_ = L();
tmp<fvScalarMatrix> fEqn
(
// - fvm::laplacian(scalar(1.0),f_)
- fvm::laplacian(f_)
==
- fvm::Sp(1.0/sqr(L_),f_)
- ((C1-scalar(6.0))*v2_/k_ - 2.0/3.0*(C1-scalar(1.0)))/(sqr(L_)*T_)
+ C2*G/(k_*sqr(L_))
/* - fvm::Sp(1.0/sqr(L_),f_)
- ((C1-scalar(1.0))*(v2_/k_-scalar(2.0/3.0))/T_
- C2*G/k_
- 5.0*epsilon_*v2_/sqr(k_))/sqr(L_)*/
);
fEqn().relax();
solve(fEqn);
bound(f_, f0_);
// v2 equation
tmp<fvScalarMatrix> v2Eqn
(
fvm::ddt(v2_)
+ fvm::div(phi_, v2_)
- fvm::laplacian(DkEff(), v2_)
==
// k_*f_
// Davidson correction - 2003
min(k_*f_,
-((C1-scalar(6.0))*v2_ - 2.0/3.0*k_*(C1-scalar(1.0)))/T_ + C2*G)
- fvm::Sp(6.0*epsilon_/k_, v2_)
);
v2Eqn().relax();
solve(v2Eqn);
bound(v2_, k0_);
// Re-calculate viscosity (with Davidson correction - 2003)
nut_ = min(CmuKE*sqr(k_)/(epsilon_ + epsilonSmall_),
Cmu*v2_*T());
}
是不是虎头蛇尾了阿？呵呵 有那么点哦，希望大家不要介意。希望对大家有所帮助。
程序和例子在附件里。
文中所涉的参考文献 均可在本文末下载
［上海市应用数学与力学研究所 SIAMM］
[ 本帖最后由 OpenFOAM 于 2009-7-8 18:01 编辑 ] |
|
[复制链接]
IP卡
狗仔卡
发表于 2009-7-8 13:40:01
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
终于解决了
