关于StochasticDispersionRAS原理疑问
-
最近在使用DPMFoam,尝试添加了随机游走模型,of版本比较老,只有StochasticDispersionRAS和GradientDispersionRAS,我使用的前者,看了看源代码,我想提出一下疑问:
通过阅读东岳老师以前的一篇文章:拉格朗日中的湍流分散力模型,我的理解是这个随机游走模型不应该是在粒子上添加一个随机游走速度吗?但是代码的最后有关Uc+UTurb,如果我理解没错,这不是连续相的速度加上一个湍流脉动速度么?难道是通过随机更改流体速度,间接再去更改颗粒受力,从而达到颗粒的随机?这个模型的原理是什么?(附上StochasticDispersionRAS.C代码)#include "StochasticDispersionRAS.H" #include "constants.H" using namespace Foam::constant::mathematical; // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::StochasticDispersionRAS ( const dictionary& dict, CloudType& owner ) : DispersionRASModel<CloudType>(dict, owner) {} template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::StochasticDispersionRAS ( const StochasticDispersionRAS<CloudType>& dm ) : DispersionRASModel<CloudType>(dm) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::~StochasticDispersionRAS() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> Foam::vector Foam::StochasticDispersionRAS<CloudType>::update ( const scalar dt, const label cellI, const vector& U, const vector& Uc, vector& UTurb, scalar& tTurb ) { cachedRandom& rnd = this->owner().rndGen(); const scalar cps = 0.16432; const scalar k = this->kPtr_->internalField()[cellI]; const scalar epsilon = this->epsilonPtr_->internalField()[cellI] + ROOTVSMALL; const scalar UrelMag = mag(U - Uc - UTurb); const scalar tTurbLoc = min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL)); // Parcel is perturbed by the turbulence if (dt < tTurbLoc) { tTurb += dt; if (tTurb > tTurbLoc) { tTurb = 0; const scalar sigma = sqrt(2*k/3.0); // Calculate a random direction dir distributed uniformly // in spherical coordinates const scalar theta = rnd.sample01<scalar>()*twoPi; const scalar u = 2*rnd.sample01<scalar>() - 1; const scalar a = sqrt(1 - sqr(u)); const vector dir(a*cos(theta), a*sin(theta), u); UTurb = sigma*mag(rnd.GaussNormal<scalar>())*dir; } } else { tTurb = GREAT; UTurb = vector::zero; } return Uc + UTurb; }