При расчете теплообмена излучением часто возникает проблема с нахождением угловых коэффициентов, которые представляют собой геометрическую характеристику взаимного положения двух поверхностей в пространстве. Если случай какой-нибудь из простых, например, соосные квадрат и круг и пр., аналитически найти угловые коэффициенты никаких трудностей не вызывает. Однако, когда мы имеем сложные криволинейные поверхности в объеме, задача становится нетривиальной и требует численных методов. В openFoam есть для этого специальная утилита viewFactorsGen, которую мы в данной статье протестируем на простейшем случае: куб внутри цилиндра с длиной L=D.
Численный расчет угловый коэффициентов (они же views factors в англоязычной литературе) производится следующим образом:
1) Разбиваются поверхности, чем меньше площадь получаемых граней тем точнее расчет
2) Составляется матрица “видимости” граней
3) Для каждой пары “видимых” граней рассчитывается угловой коэффициент с учетом их площади, расстояния и углом между ними
4) Для каждой поверхности угловые коэффициенты суммируются.
Для двух поверхностей, с количество элементарных площадок ~10000, получаются матрицы 10000×10000. Чтобы сгруппировать грани и уменьшить время расчета используется утилита faceAgglomerate, однако в этом есть существенная проблема. В исходном коде данной утилиты допускается максимальное количество элементов (переменная maxDynListLength) в количестве 100000, это чаще всего не хватает, поэтому необходимо это исправить. Идем “/applications/utilities/preProcessing/viewFactorsGen” и изменяем shootRays.H. Можно просто задать большое число для этой константы, но вот здесь предложен более удобный вариант, когда переменная берется из constant/viewFactorsGenDict:
|
1 2 3 4 |
//const label maxDynListLength = 100000; const label maxDynListLength = viewFactorDict.lookupOrDefault<label>("maxDynListLength", 100000); if (Pstream::master()) { Info << "\n maxDynListLength: " << maxDynListLength << endl; } |
Все рассмотренные в данной статье случаи в приложенном архиве: radiationCubeInCylinder.tar.gz
В папке начальных значений в нашем случае используется только один файл Qr – тепловой поток излучением:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
// размерность Вт/м^2 = кг/с^3 dimensions [ 1 0 -3 0 0 0 0 ]; internalField uniform 0; boundaryField { "cube|cylinder" { // другие типы граничных условий можно посмотреть в src/thermophysicalModels/radiationModels/derivedFvPatchFields type greyDiffusiveRadiationViewFactor; value uniform 0; // режим излучения: либо solidRadiation (для сопряженных задач) либо lookup (только поверхности) emissivityMode lookup; // излучательная способность emissivity uniform 1; // внешний тепловой поток Qro uniform 0; } } |
Единственный файл настроек – constant/viewFactorsDict:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
// записывать в 0 матрицы facesAgglomeration и viewFactorField (если надо - ставим true) writeViewFactorMatrix false; writeFacesAgglomeration false; // та самая переменная, о которой говорили выше maxDynListLength 100000000; // создавать в папке с задачей файлы *.obj и *.vtk для визуализации dumpRays false; // настройки поверхностей cube { // количество граней после укрупнения faceAgglomerate-ом (чем больше тем точнее) nFacesInCoarsestLevel 1000; // минимальный угол "видимости" граней featureAngle 0.1; } cylinder { nFacesInCoarsestLevel 4000; featureAngle 0.1; } |
В результате работы утилиты viewFactorsGen в лог запишутся необходимые нам угловые коэффициенты:
|
1 2 3 4 |
F00: 1.65174e-19 F01: 0.998229 F10: 0.496542 F11: 0.506226 |
В нашем случае поверхность куба имеет номер 0, поверхность цилиндра 1. Таким образом:
F00 – угловой коэффициент куба самого на себя (т.к. он не имеет выпуклых поверхностей F00=0)
F01 – угловой коэффициент куба на цилиндр (т.к. цилиндр полностью окружает куб F01=1)
F10 – угловой коэффициент цилиндра на куб, аналитически рассчитывается по формуле
F11 – угловой коэффициент цилиндра сам на себя, аналитически рассчитывается по формуле
Сопоставление коэффициента F10 полученного численно и аналитически выполнено на графике ниже, расхождение менее 0.2%.
