Некоторое время назад я искал описание операций, процессов, происходящих в библиотеке численного моделирования OpenFOAM. Нашел много абстрактных описаний работы метода конечных объёмов, классических разностных схем, различных физических уравнений. Мне же хотелось узнать более детально — откуда в таком-то выходном файле на такой-то итерации получились эти значения, какие выражения стоят за теми или иными параметрами в файлах настроек fvSchemes, fvSolution?
Для тех, кому это тоже интересно — эта статья. Те, кто хорошо знаком с OpenFOAM или с методами в нём реализованными — пишите о найденных ошибках и неточностях в личку.
На хабре уже была пара статей про OpenFOAM:
OpenFOAM на практике
OpenFOAM с точки зрения программиста-физика
Поэтому не буду останавливаться на том, что это «открытая (GPL) платформа для численнного моделирования, предназначенная для моделирования, связанного с решением уравнений в частных производных методом конечных объёмов, и широко используемая для решения задач механики сплошных сред».
Сегодня я на простом примере опишу операции, которые происходят в ходе расчёта в OpenFOAM.
Итак, дана геометрия — куб со стороной 1 метр:
Перед нами стоит задача смоделировать поток-распространение некоторого скалярного поля (температура, количество вещества), который задаётся следующим уравнением переноса (1) внутри объема тела.
(1)
,
где
скалярная величина, например, выражает температуру [K] или концентрацию некоторого вещества, а 
выражает перенос вещества, массовый поток [кг/с].
Это уравнение, например, используют для моделирования распространения тепла
,
где k — теплопроводность, а
— температура [K].
Граничные условия задачи следующие: есть грань-вход, грань-выход, остальные грани — гладкие стенки.
Наша сетка будет очень простая — делим куб на 5 равных ячеек вдоль оси Z.
Метод конечных объёмов предусматривает, что (1) в интегральной форме (2) будет выполняться для каждого конечного объёма
.
(2)
,
где
— геометрический центр конечного объёма.
Упростим, преобразуем первое слагаемое выражения (2) следующим образом:
(2.1) (HJ-3.12)*
Как видно — мы приняли, что скалярная величина изменяется внутри конечного объема линейно и значение величины
в некоторой точке 
внутри конечного объёма можно вычислить как:
(2.2)
,
где
Для упрощения второго слагаемого выражения (2) используем обобщённую теорему Гаусса-Остроградского: интеграл от дивергенции векторного поля
по объёму 
, равен потоку вектора через поверхность 
, ограничивающую данный объём. На человеческом языке «сумма всех потоков в/из конечного объема равна сумме потоков через грани этого конечного объема»:
(2.3)
,
где
замкнутая поверхность, ограничивающая объём ,
— вектор, направленный по нормали от объёма 
.
Учитывая то, что конечный объём ограничен набором плоских граней, можно выражение (2.3) преобразовать к сумме интегралов по поверхности:
(2.4) (HJ-3.13)
,
где
выражает значение переменной в центре грани,
— вектор площади, выходит из центра грани, направлен в сторону от ячейки (локально), в сторону от ячейки с меньшим индексом к ячейке с большим индексом (глобально).
Чтобы не хранить одни и те же параметры вектора два раза, т.к. очевидно, что у двух соседних ячеек вектор-нормали к грани между ячейками, направленный в сторону От центра ячейки будет различаться только направлением-знаком. Поэтому было создано owner-neighbour отношение между гранью и ячейкой. Если вектор площади (глобальный, положительное направление от ячейки с меньшим индексом к ячейке с большим индексом) указывает ОТ центра 
ячейки такое отношение между ячейкой и вектором , а точнее между ячейкой и гранью, обозначается owner). В случае если этот вектор 
указывает внутрь рассматриваемой ячейки, то отношение neighbour. Направление влияет на знак величины (+ для owner и — для neighbour) и это важно при суммировании см. далее.
(HJ-3.16)
Значение в центре грани вычисляется через значения в центрах прилегающих ячеек — способ такого выражения носит название разностной схемы. В OpenFOAM тип разностной схемы задается в файле /system/fvSchemes:
Gauss — означает, что выбрана центральная разностная схема;
linear — означает, что интерполяция с центров ячеек на центры граней будет происходить линейно.
Допустим, что наша скалярная величина изменяется внутри конечного объема от центра к граням линейно. Тогда аппроксимированное в центре грани значение будет вычисляться согласно формуле:
Где веса
и 
рассчитываются как
Где
— объемы ячеек.
Для случаев скошенных ячеек существуют более сложные формулы расчета весов аппроксимации.
Таким образом, значения phi_f в центрах граней ячеек вычисляются на основе значений в центрах ячеек. Значения градиентов grad(phi) вычисляются на основе значений phi_f.
И весь этот алгоритм может быть представлен в виде следующего псевдокода.
Учитывая (2.1) и (2.4) выражение (2) принимает вид:
(3)
Согласно методу конечных объёмов проводится дискретизация по времени и выражение (3) записывается как:
(4)
Проинтегрируем (4):
(4.1)
Разделим левую и правую часть на
:
(5)
Теперь мы можем получить систему линейных уравнений для каждого конечного объёма
.
Ниже представлена нумерация узлов сетки, которую мы будем использовать.
Нумерация узлов-центров ячеек (50, 51 — центры граничных граней):
Нумерация узлов-центров граней:
Объемы элементов:
Коэффициенты интерполяции, необходимые для вычисления значений на гранях ячеек. Индекс «e» обозначает «правая грань ячейки». Правая относительно вида, как на рисунке «Нумерация узлов-центров ячеек»:
Для P = 0.
Выражение (5), описывающее поведение величины
будет преобразовано в систему линейных алгебраических уравнений, каждое из которых вида:
или, согласно индексам точек на гранях
А еще все потоки в/из ячейки могут быть выражены в виде суммы
где, например,
— коэффициент линеаризации потока в точке-центре ячейки E,
— коэффициент линеаризации потока в точке-центре грани,
— нелинейная часть (например, константа).
Согласно нумерации граней выражение примет вид:
С учетом граничных условий для элемента P_0 линейное алгебраическое уравнение может быть представлено в виде
… подставим ранее полученные коэффициенты…
Поток из inlet'a направлен в ячейку, поэтому имеет отрицательный знак.
Так как у нас в управляющем выражении присутствует кроме диффузионного еще и временной член, но конечное уравнение выглядит как
Для P = 1.
Для P = 4.
Система линейных алгебраических уравнений (СЛАУ) может быть представлена в матричном виде как
,
где
Далее полученная СЛАУ решается решателем, указанным в fvSchemes.
И в итоге получается вектор значений
на основе которого получаются значения для вектора
Затем вектор подставляется в СЛАУ и происходит новая итерация расчёта вектора .
И так до тех пор, пока невязка не достигнет требуемых пределов.
* Некоторые уравнения в этой статье взяты из диссертации Ясака Хрвое (HJ — номер уравнения) и если кому-то захочется прочитать про них подробнее (http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdf)
Скачать файлы задачи можно здесь:
github.com/j-avdeev/case
Файлы решателя:
github.com/j-avdeev/matrHyper1Foam
В качестве бонуса — видео, как распространяется концентрация.
Для тех, кому это тоже интересно — эта статья. Те, кто хорошо знаком с OpenFOAM или с методами в нём реализованными — пишите о найденных ошибках и неточностях в личку.
На хабре уже была пара статей про OpenFOAM:
OpenFOAM на практике
OpenFOAM с точки зрения программиста-физика
Поэтому не буду останавливаться на том, что это «открытая (GPL) платформа для численнного моделирования, предназначенная для моделирования, связанного с решением уравнений в частных производных методом конечных объёмов, и широко используемая для решения задач механики сплошных сред».
Сегодня я на простом примере опишу операции, которые происходят в ходе расчёта в OpenFOAM.
Итак, дана геометрия — куб со стороной 1 метр:
Перед нами стоит задача смоделировать поток-распространение некоторого скалярного поля (температура, количество вещества), который задаётся следующим уравнением переноса (1) внутри объема тела.
(1)
,где
скалярная величина, например, выражает температуру [K] или концентрацию некоторого вещества, а выражает перенос вещества, массовый поток [кг/с].Это уравнение, например, используют для моделирования распространения тепла
,где k — теплопроводность, а
— температура [K].Оператор дивергенции на самом деле это
оператор 
.
Напомню, что существует оператор набла (оператор Гамильтона), который записывается следующим образом:
,
где i, j, k — единичные векторы.
Если скалярно умножить оператор набла на векторную величину, то мы получим дивергенцию данного вектора:
«С точки зрения физики, дивергенция векторного поля является показателем того, в какой степени данная точка пространства является источником или стоком этого поля»
Если умножить оператор набла на скаляр, получается градиент этого скаляра:
Градиент показывает увеличение или уменьшение по какому-либо направлению величины скаляра
.
.Напомню, что существует оператор набла (оператор Гамильтона), который записывается следующим образом:
,где i, j, k — единичные векторы.
Если скалярно умножить оператор набла на векторную величину, то мы получим дивергенцию данного вектора:
«С точки зрения физики, дивергенция векторного поля является показателем того, в какой степени данная точка пространства является источником или стоком этого поля»
Если умножить оператор набла на скаляр, получается градиент этого скаляра:
Градиент показывает увеличение или уменьшение по какому-либо направлению величины скаляра
.Граничные условия задачи следующие: есть грань-вход, грань-выход, остальные грани — гладкие стенки.
Разбиение объема куба на конечные объемы
Наша сетка будет очень простая — делим куб на 5 равных ячеек вдоль оси Z.
Много формул
Метод конечных объёмов предусматривает, что (1) в интегральной форме (2) будет выполняться для каждого конечного объёма
.(2)
,где
— геометрический центр конечного объёма.Центр конечного объема
Упростим, преобразуем первое слагаемое выражения (2) следующим образом:
(2.1) (HJ-3.12)*
Как видно — мы приняли, что скалярная величина изменяется внутри конечного объема линейно и значение величины
в некоторой точке внутри конечного объёма можно вычислить как:(2.2)
,где
Для упрощения второго слагаемого выражения (2) используем обобщённую теорему Гаусса-Остроградского: интеграл от дивергенции векторного поля
по объёму , равен потоку вектора через поверхность , ограничивающую данный объём. На человеческом языке «сумма всех потоков в/из конечного объема равна сумме потоков через грани этого конечного объема»:(2.3)
,где
замкнутая поверхность, ограничивающая объём , — вектор, направленный по нормали от объёма .Вектор S
Учитывая то, что конечный объём ограничен набором плоских граней, можно выражение (2.3) преобразовать к сумме интегралов по поверхности:
(2.4) (HJ-3.13)
,где
выражает значение переменной в центре грани, — вектор площади, выходит из центра грани, направлен в сторону от ячейки (локально), в сторону от ячейки с меньшим индексом к ячейке с большим индексом (глобально).Еще немного про вектор S
Чтобы не хранить одни и те же параметры вектора
два раза, т.к. очевидно, что у двух соседних ячеек вектор-нормали к грани между ячейками, направленный в сторону От центра ячейки будет различаться только направлением-знаком. Поэтому было создано owner-neighbour отношение между гранью и ячейкой. Если вектор площади (глобальный, положительное направление от ячейки с меньшим индексом к ячейке с большим индексом) указывает ОТ центра ячейки такое отношение между ячейкой и вектором , а точнее между ячейкой и гранью, обозначается owner). В случае если этот вектор указывает внутрь рассматриваемой ячейки, то отношение neighbour. Направление влияет на знак величины (+ для owner и — для neighbour) и это важно при суммировании см. далее.(HJ-3.16)
Про разностные схемы
Значение
в центре грани вычисляется через значения в центрах прилегающих ячеек — способ такого выражения носит название разностной схемы. В OpenFOAM тип разностной схемы задается в файле /system/fvSchemes:divSchemes
{
default none;
div(phi,psi) Gauss linear;
}
Gauss — означает, что выбрана центральная разностная схема;
linear — означает, что интерполяция с центров ячеек на центры граней будет происходить линейно.
Допустим, что наша скалярная величина изменяется внутри конечного объема от центра к граням линейно. Тогда аппроксимированное в центре грани значение будет вычисляться согласно формуле:
Где веса
и рассчитываются какГде
— объемы ячеек.Для случаев скошенных ячеек существуют более сложные формулы расчета весов аппроксимации.
Таким образом, значения phi_f в центрах граней ячеек вычисляются на основе значений в центрах ячеек. Значения градиентов grad(phi) вычисляются на основе значений phi_f.
И весь этот алгоритм может быть представлен в виде следующего псевдокода.
1. Объявляем массив градиентов конечных объемов, инициализируем его нулями
2. Пробегаемся по всем внутренним граням (которые не граничные)
> Вычисляем flux_f = phi_f*S_f. Значения phi_f вычисляем на основе значений phi в центах ячеек
> Добавляем flux_f к градиенту элемента-owner и -flux_f к градиенту элемента-neighbour
3. Пробегаемся по всем граничным граням
> Вычисляем flux_f = phi_f*S_f
> Добавляем flux_f к градиенту элементу-owner (neighbour-элементов у граничных граней нет)
4. Пробегаемся по всем элементам
> Делим получившуюся сумму-градиент на объем элемента
Дискретизация по времени
Учитывая (2.1) и (2.4) выражение (2) принимает вид:
(3)
Согласно методу конечных объёмов проводится дискретизация по времени и выражение (3) записывается как:
(4)
Проинтегрируем (4):
(4.1)
Разделим левую и правую часть на
:(5)
Данные для матрицы дискретизации
Теперь мы можем получить систему линейных уравнений для каждого конечного объёма
.Ниже представлена нумерация узлов сетки, которую мы будем использовать.
Координаты узлов хранятся в /constant/polyMesh/points
24
(
(0 0 0)
(1 0 0)
(0 1 0)
(1 1 0)
(0 0 0.2)
(1 0 0.2)
(0 1 0.2)
(1 1 0.2)
(0 0 0.4)
(1 0 0.4)
(0 1 0.4)
(1 1 0.4)
(0 0 0.6)
(1 0 0.6)
(0 1 0.6)
(1 1 0.6)
(0 0 0.8)
(1 0 0.8)
(0 1 0.8)
(1 1 0.8)
(0 0 1)
(1 0 1)
(0 1 1)
(1 1 1)
)
Нумерация узлов-центров ячеек (50, 51 — центры граничных граней):
Нумерация узлов-центров граней:
Объемы элементов:
Коэффициенты интерполяции, необходимые для вычисления значений на гранях ячеек. Индекс «e» обозначает «правая грань ячейки». Правая относительно вида, как на рисунке «Нумерация узлов-центров ячеек»:
Формирование матрицы дискретизации
Для P = 0.
Выражение (5), описывающее поведение величины
будет преобразовано в систему линейных алгебраических уравнений, каждое из которых вида:
или, согласно индексам точек на гранях
А еще все потоки в/из ячейки могут быть выражены в виде суммы
где, например,
— коэффициент линеаризации потока в точке-центре ячейки E, — коэффициент линеаризации потока в точке-центре грани, — нелинейная часть (например, константа).Согласно нумерации граней выражение примет вид:
С учетом граничных условий для элемента P_0 линейное алгебраическое уравнение может быть представлено в виде
… подставим ранее полученные коэффициенты…
Поток из inlet'a направлен в ячейку, поэтому имеет отрицательный знак.
Так как у нас в управляющем выражении присутствует кроме диффузионного еще и временной член, но конечное уравнение выглядит как
Для P = 1.
Для P = 4.
Система линейных алгебраических уравнений (СЛАУ) может быть представлена в матричном виде как
,где
=== A(i,j) ===
40.5 0.5 0 0 0
-0.5 40 0.5 0 0
0 -0.5 40 0.5 0
0 0 -0.5 40 0.5
0 0 0 -0.5 40.5
=== b(i,j) ===
1
0
0
0
0
Далее полученная СЛАУ решается решателем, указанным в fvSchemes.
И в итоге получается вектор значений
psi = dimensions [0 0 0 0 0 0 0];
internalField nonuniform List<scalar> 5(0.0246875 0.000308546 3.85622e-06 4.81954e-08 5.95005e-10);
на основе которого получаются значения для вектора
Затем вектор
подставляется в СЛАУ и происходит новая итерация расчёта вектора .И так до тех пор, пока невязка не достигнет требуемых пределов.
Ссылки
* Некоторые уравнения в этой статье взяты из диссертации Ясака Хрвое (HJ — номер уравнения) и если кому-то захочется прочитать про них подробнее (http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdf)
Скачать файлы задачи можно здесь:
github.com/j-avdeev/case
Файлы решателя:
github.com/j-avdeev/matrHyper1Foam
В качестве бонуса — видео, как распространяется концентрация
.