Стационарное ламинарное течение в трубе

В данной статье рассмотрю каким образом считается простейшее течение. В качестве примера возьмем течение на начальном участке трубы. Отсутствие необходимости определять турбулентные характеристики значительно упрощает задачу, т.к. решается обычное уравнение Навье-Стокса, без турбулентной вязкости или напряжений Рейнольдса:
  \frac {dV} {dt} = F - \frac {1} {\rho} \nabla p + \nu \nabla ^2 V
где V – скорость потока, p – давление, F – массовые силы, t – время, \rho – плотность, \nu – молекулярная вязкость.
В статье сделаем следующее:
1) Сопоставим численное и аналитическое решение увеличения скорости на оси трубы;
2) Сопоставим параболический профиль с профилем скорости после смыкания пограничных слоев;
3) Сравним скорость расчета при использовании различных методов решения матриц.
plot_XY_tube_reedit

Первая проблема, возникающая при моделировании любого течения, это сложность получения связанных полей скорости-давления. Поле давления является искомыми, при этом оно не выражено явно и находится из уравнения неразрывности (т.е если “правильное” поле давления подставить в уравнение Навье-Стокса, можно получить поле скоростей, удовлетворяющее уравнению неразрывности):
  \nabla V = 0
Решение следующее:
1) Предварительное задание поле давление p*;
2) Определение поля скоростей u*;
3) По данному полю скоростей вычисляется поправка к исходному полю давлений p';
4) Добавляем поправку давлений p’ к p*;
5) Расчет u по новому полю давлений;
6) Возвращаемся к п.2.
И так повторяем пока расхождение не будет меньше 1e-5. Сам метод SIMPLE хорошо изложен в переводе книги Патанкара [1].

Ламинарное течение в трубе характерно следующим:
1) На начальном участке происходит образование пограничных слоев, которые “вытесняют” поток газа и происходит увеличение скорости на оси;
2) В установившемся режиме профиль – параболический, скорость на оси в 2 раза выше средней.

См. приложенный файл laminarTube.tar.gz

Примем следующую геометрию сетки: диаметр – 15 мм, длина – 140 мм. Принципиальная схема следующая:
Схема

Ниже приведена расчетная сетка, как видно к стенкам толщина ячеек уменьшается, для турбулентных течений необходимо соблюдать y+<1. Для ламинарных я нигде такого требования не нашел, но очевидно, что именно у стенок происходит наибольшее изменение скорости и именно там необходимо уменьшать размеры ячеек.
plot_XYplot_XYZ
Для начальных граничных условий наверно пояснений не надо. На входе задаем фиксированное значение скорости и нулевой градиент давления, на выходе нулевой градиент скорости и фиксированное значение давление (в примере 0, но можно поставить например атмосферное, разницы нет), на стенках нулевая скорость и нулевой градиент давления.
Вообще simpleFoam используется для несжимаемых турбулентных течений, но его можно использовать и для ламинарных, для этого надо в constant/RASProperties указать

и в constant/turbulenceProperties указываем

.
Не смотря на то, что солвер стационарный указывать временной шаг все-равно надо, но он не несет никакой смысловой нагрузки, разве что названия папок, system/controlDict:

Решение сходится замечательно, поэтому можно не мелочиться и указывать самые точные схемы дискретизации дифференциальных уравнений в system/fvSchemes:

Основные особенности заключаются в файле system/fvSolvers:

Решение сходится после 524 итераций, сделал небольшую серию, чтобы определить наилучшие параметры решения матриц

№ п/п           Матрица p         Матрица U                     Время
Опыт1          PCGDIC                smoothSolver               982,89с
Опыт2          PCGDIC                DILUPBiCG                     947,25c
Опыт3          GAMG                  smoothSolver               626,46c //nCellsInCoarsestLevel 400;
Опыт4          GAMG                  DILUPBiCG                     614,41c //nCellsInCoarsestLevel 400;
Опыт5          GAMG                  GAMG                              702,29c //nCellsInCoarsestLevel 100;
Опыт6          GAMG                  GAMG                              666,56c //nCellsInCoarsestLevel 250;
Опыт7          GAMG                  GAMG                              617,04c //nCellsInCoarsestLevel 400;
Опыт8          GAMG                  GAMG                              614,29c //nCellsInCoarsestLevel 550;

Видно, что опыты 4,8 самые быстрые (в примере аналогично опыту 4), их и рекомендую использовать.
Ниже приведен график сходимости до 1e-6 полученный через foamLog, видно монотонная сходимость, как и нужно:

residual

Профиль скорости снимаем по оси трубы и по радиусу в конце трубы:

Сопоставим профиль скорости с параболическим:

u(r)

Как видно сходится хорошо, а вот сопоставление изменения скорости по оси не очень, не известно из-за чего, быть может из-за допущений в аналитическом расчете [2]:

u0(x)

plot_XY_tube_reedit
Литература:
1. Патанкар С. – Численные методы решения задач теплообмена и динамики жидкости, 1984 г. – М::Энергоатомиздат.
2. Слезкин Н.А, – Динамика вязкой жидкости. 1955 -М::Гос. изд-во.тех.-теор. лит-ры.