Продолжу рассматривать численное моделирование аэродинамики циклонных устройств. Начало здесь. Как я отмечал ранее, использование “классических” моделей турбулентности kOmegaSST и SA в их первоначальном виде приводит к завышению турбулентной вязкости и значительному смещению максимума тангенциальной скорости к боковой поверхности камеры. В данной статье введем в эти модели поправки на кривизну линий тока и сравним полученные профили безразмерной тангенциальной скорости с экспериментальными.
Математическую формулировку поправки на кривизну линий тока (поправка Спалларта-Шура) можно посмотреть на сайте NASA, или в [1, с.58]
Изменим исходный код оригинальной модели Спалларта-Алмараса (src/turbulenceModels/incompressible/RAS/SpalartAllmaras). Для этого вычислим поправку fr1 и в самой модели умножим на нее генерационный член
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
tmp< // находим тензор скоростей деформации volSymmTensorField> strainRateTensor (symm(fvc::grad(U_))); // находим тензор завихренности tmp<volTensorField> vorticityTensor (skew(fvc::grad(U_))); // находим модули (??) тензоров скоростей деформации и завихренности volScalarField S (sqrt(2.0)*mag(strainRateTensor())); // т.к. дальше модуль тензора завихренности используется в знаменателе он должен отличаться от нуля volScalarField omegaCaps (max(sqrt(2.0)*mag(vorticityTensor()),dimensionedScalar("minOmegaCaps",S.dimensions(),SMALL))); volScalarField D (sqrt(0.5*(sqr(S)+sqr(omegaCaps)))); // находим параметр r со звездой volScalarField rStar (S/omegaCaps); // находим полную производную тензора скорости деформации tmp<volSymmTensorField> DSDt (fvc::ddt(strainRateTensor())+fvc::div(phi_, strainRateTensor())); // находим параметр r с тильдой volScalarField rTilda (2.0*((strainRateTensor() & vorticityTensor()) && DSDt())/(pow(D,4))); // находим сам параметр, характеризующий закрутку потока volScalarField fr1 ((1+cr1_)/(1+rStar)*(2*rStar*(1-cr3_*atan(cr2_*rTilda)))-cr1_); //обнуляем ненужные переменные strainRateTensor.clear(); vorticityTensor.clear(); DSDt.clear(); S.clear(); omegaCaps.clear(); D.clear(); rStar.clear(); rTilda.clear(); |
Добавляем константы cr1, cr2 и cr3:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
cr1_ ( dimensioned<scalar>::lookupOrAddToDict ( "cr1", coeffDict_, 1 ) ), cr2_ ( dimensioned<scalar>::lookupOrAddToDict ( "cr2", coeffDict_, 12 ) ), cr3_ ( dimensioned<scalar>::lookupOrAddToDict ( "cr3", coeffDict_, 1 ) ), |
И умножаем генерационный член в правой части на найденный коэффициент fr1
|
1 2 3 4 5 6 7 8 9 10 |
tmp<fvScalarMatrix> nuTildaEqn ( fvm::ddt(nuTilda_) + fvm::div(phi_, nuTilda_) - fvm::laplacian(DnuTildaEff(), nuTilda_) - Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_)) == fr1*Cb1_*Stilda*nuTilda_ - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_/sqr(d_), nuTilda_) ); |
Для удобства меняем название модели турбулентности с SpalartAllmaras на любое другое, например SARC.
Аналогично поступаем и с моделью kOmegaSST. Математическую формулировку также можно найти на сайте NASA здесь. Эта поправка очень похожа на приведенную выше, поэтому без подробных комментариев.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
tmp<volSymmTensorField> strainRateTensor (symm(fvc::grad(U_))); tmp<volTensorField> vorticityTensor (skew(fvc::grad(U_))); volScalarField S (sqrt(2.0)*mag(strainRateTensor())); volScalarField omegaCaps (max((sqrt(2.0)*mag(vorticityTensor())),dimensionedScalar("minOmegaCaps", S.dimensions(),SMALL))); volScalarField D (sqrt(max(sqr(S),0.09*sqr(omega_)))); volScalarField rStar (S/omegaCaps); tmp<volSymmTensorField> DSDt (fvc::ddt(strainRateTensor())+fvc::div(phi_, strainRateTensor())); volScalarField rTilda (2.0*((strainRateTensor() & vorticityTensor()) && DSDt())/(omegaCaps*pow(D,3))); volScalarField fRotation ((1+cr1_)/(1+rStar)*(2*rStar*(1-cr3_*atan(cr2_*rTilda)))-cr1_); // для визуализации можно сделать запись коэффициента fr1 вместе с другими параметрами volScalarField fr1( IOobject( "fr1", runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::AUTO_WRITE ), (max(min(fRotation,scalar(1.25)),scalar(0.0))) ); if(runTime_.outputTime()) { fr1.write(); } strainRateTensor.clear(); vorticityTensor.clear(); DSDt.clear(); S.clear(); omegaCaps.clear(); D.clear(); rStar.clear(); rTilda.clear(); fRotation.clear(); |
Все исходные файл см. в приложенном архиве RCModels.tar. Распаковываем в любом удобном месте, заходим в директорию ../incompressible, пишем:
|
1 |
wmake libso |
Если при компиляции ошибок не возникло в /platforms/linux64GccDPOpt/lib должна появиться библиотека libincompressibleRASModels_RCModel.so. Чтобы ее использовать в controlDict кейса необходимо ее указать (см. приложенные файлы kOmegaSST.tar.gz kOmegaSSTRC.tar.gz SA.tar.gz SARC.tar.gz:
|
1 2 3 |
libs ( "libincompressibleRASModels_RCModel.so" ); |
Теперь можно использовать наши модифицированные модели турбулентности:
|
1 |
RASModel SARC; |
Сопоставление радиальных профилей тангенциальной скорости при dвых=0.3 выполнено на графиках ниже, видно что введение поправок на кривизну линий тока существенно увеличивается точность расчета: появляется ярко выраженный максимум скорости, положение максимума смещается к оси камеры.
Литература:
1. Смирнов, Е.М. Конспект лекций дисциплины “Течения вязкой жидкости и модели турбулентности: методы расчета турбулентных течений” / Е.М. Смирнов, А.В. Гарбарук- 2010. – 127 с.