|
[Sponsors] |
Some confusions about SST-DES (code and formulas) |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
April 21, 2018, 08:45 |
Some confusions about SST-DES (code and formulas)
|
#1 |
Member
王莹
Join Date: May 2017
Posts: 51
Rep Power: 9 |
Hello,
I'm trying to modify the SST-DES in OpenFOAM. I have found the main formulas of SST-DES and its code in OpenFOAM. However, when I compared the two parts so as to match the code with formulas one to one, it really confused me.....I can't even locate the key variables in the code. I wonder if I found the wrong formulas of SST-DES, so I enclose the formulas and code. Please, any suggestion is welcome. p.s. So sorry that I can't post the .doc file......Size of the file is limited. |
|
April 22, 2018, 02:27 |
|
#2 |
Member
Fatih Ertinaz
Join Date: Feb 2011
Location: Istanbul
Posts: 64
Rep Power: 15 |
Alisa
Short answer: Check Code:
$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C Find the constructor in the code you attached: Code:
template<class BasicTurbulenceModel> kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES ( someParametersHere ) kOmegaSST < LESeddyViscosity<BasicTurbulenceModel>, BasicTurbulenceModel > You can find the constructor implementation of the template class in $FOAM_SRC/TurbulenceModels/ turbulenceModels/* So now it is crucial to figure out how kOmegaSST is implemented. That information is in Code:
$FOAM_SRC/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C Code:
template<class BasicTurbulenceModel> kOmegaSST<BasicTurbulenceModel>::kOmegaSST ( some parameters ) : Foam::kOmegaSST < eddyViscosity<RASModel<BasicTurbulenceModel>>, BasicTurbulenceModel > This is the mathematical model for k-omega SST: So let's start with turbulence kinetic energy: Code:
// Turbulent kinetic energy equation tmp<fvScalarMatrix> kEqn ( fvm::ddt(alpha, rho, k_) + fvm::div(alphaRhoPhi, k_) - fvm::laplacian(alpha*rho*DkEff(F1), k_) == alpha()*rho()*Pk(G) - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_) - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_) + kSource() + fvOptions(alpha, rho, k_) ); F1 is the blending factor which is in charge of the transition between boundary layer (=1) and far field (=0). This is implemented as a protected member function and code below clearly corresponds to the F1 formula written above: Code:
tmp<volScalarField> CDkOmegaPlus = max ( CDkOmega, dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10) ); tmp<volScalarField> arg1 = min ( min ( max ( (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_), scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) ), (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_)) ), scalar(10) ); return tanh(pow4(arg1)); Code:
volScalarField CDkOmega ( (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_ ); Code:
//- Return the effective diffusivity for k tmp<volScalarField> DkEff(const volScalarField& F1) const { return tmp<volScalarField> ( new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu()) ); } Code:
tmp<volScalarField> alphaK(const volScalarField& F1) const { return blend(F1, alphaK1_, alphaK2_); } Code:
tmp<volScalarField::Internal> blend ( const volScalarField::Internal& F1, const dimensionedScalar& psi1, const dimensionedScalar& psi2 ) const { return F1*(psi1 - psi2) + psi2; } So now we should take a look at the definition of this->nut_. It is called kinematic eddy viscosity and is defined as: and the corresponding code is: Code:
this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2)); Code:
tmp<volScalarField> arg2 = min ( max ( (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_), scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) ), scalar(100) ); return tanh(sqr(arg2)); Code:
volScalarField S2(2*magSqr(symm(tgradU()))); I think this explains the third term. In general what you see in the code is hidden behind the DkEff function. When you dig deeper it becomes clearer. The first term on the RHS P_k is also explicitly implemented. You can find in the same file: Code:
$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C Code:
return min(G, (c1_*betaStar_)*this->k_()*this->omega_()); Code:
volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()()); volScalarField::Internal G(this->GName(), nut()*GbyNu); Now remember the definition of the tgradU that is written above. I don't have enough energy and time to look at the details of the GbyNu but it corresponds to the shear tensor. So the final term : Code:
- fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_) - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_) See epsilonByk function: Code:
return betaStar_*omega_(); References: + https://pingpong.chalmers.se/public/...o?item=3256524 + https://pingpong.chalmers.se/public/...OmegaSST-1.pdf Last edited by fertinaz; April 22, 2018 at 02:52. Reason: Added references |
|
April 24, 2018, 10:30 |
|
#3 | |
Member
王莹
Join Date: May 2017
Posts: 51
Rep Power: 9 |
Dear Fertinaz,
Thank you so much for your detailed explanation and I am also sorry for my late reply. After reading your reply, I have really known more about the frame of OF, which contains so much inheritance and deriving that I sometimes can not find the root of certain variable. Of course, more important reason is that my knowledge on CFD is very poor. Again, thank you so much. I greatly appreciate your generosity and patience. Quote:
|
||
Tags |
code, des, modify, sst, turbulence model |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Some confusions about SST-DES (code and formulas) | Alisa_W | OpenFOAM Programming & Development | 1 | April 12, 2018 23:29 |