Pull to refresh

OpenFOAM для чайников

Reading time 6 min
Views 22K
Некоторое время назад я искал описание операций, процессов, происходящих в библиотеке численного моделирования OpenFOAM. Нашел много абстрактных описаний работы метода конечных объёмов, классических разностных схем, различных физических уравнений. Мне же хотелось узнать более детально — откуда в таком-то выходном файле на такой-то итерации получились эти значения, какие выражения стоят за теми или иными параметрами в файлах настроек fvSchemes, fvSolution?
Для тех, кому это тоже интересно — эта статья. Те, кто хорошо знаком с OpenFOAM или с методами в нём реализованными — пишите о найденных ошибках и неточностях в личку.

На хабре уже была пара статей про OpenFOAM:
OpenFOAM на практике
OpenFOAM с точки зрения программиста-физика

Поэтому не буду останавливаться на том, что это «открытая (GPL) платформа для численнного моделирования, предназначенная для моделирования, связанного с решением уравнений в частных производных методом конечных объёмов, и широко используемая для решения задач механики сплошных сред».

Сегодня я на простом примере опишу операции, которые происходят в ходе расчёта в OpenFOAM.

Итак, дана геометрия — куб со стороной 1 метр:





Перед нами стоит задача смоделировать поток-распространение некоторого скалярного поля (температура, количество вещества), который задаётся следующим уравнением переноса (1) внутри объема тела.

(1)
image,

где image скалярная величина, например, выражает температуру [K] или концентрацию некоторого вещества, а image выражает перенос вещества, массовый поток [кг/с].

Это уравнение, например, используют для моделирования распространения тепла
,
где k — теплопроводность, а — температура [K].

Оператор дивергенции на самом деле это
оператор image.
Напомню, что существует оператор набла (оператор Гамильтона), который записывается следующим образом:
image,

где i, j, k — единичные векторы.
Если скалярно умножить оператор набла на векторную величину, то мы получим дивергенцию данного вектора:
image
«С точки зрения физики, дивергенция векторного поля является показателем того, в какой степени данная точка пространства является источником или стоком этого поля»

Если умножить оператор набла на скаляр, получается градиент этого скаляра:

image
Градиент показывает увеличение или уменьшение по какому-либо направлению величины скаляра image.

Граничные условия задачи следующие: есть грань-вход, грань-выход, остальные грани — гладкие стенки.



Разбиение объема куба на конечные объемы


Наша сетка будет очень простая — делим куб на 5 равных ячеек вдоль оси Z.



Много формул


Метод конечных объёмов предусматривает, что (1) в интегральной форме (2) будет выполняться для каждого конечного объёма image.

(2)
image,

где image — геометрический центр конечного объёма.

Центр конечного объема



Упростим, преобразуем первое слагаемое выражения (2) следующим образом:

(2.1) (HJ-3.12)*
image

image

image

Как видно — мы приняли, что скалярная величина изменяется внутри конечного объема линейно и значение величины image в некоторой точке image внутри конечного объёма можно вычислить как:

(2.2)
image,

где image

Для упрощения второго слагаемого выражения (2) используем обобщённую теорему Гаусса-Остроградского: интеграл от дивергенции векторного поля image по объёму image, равен потоку вектора через поверхность image, ограничивающую данный объём. На человеческом языке «сумма всех потоков в/из конечного объема равна сумме потоков через грани этого конечного объема»:

(2.3)
image,

где image замкнутая поверхность, ограничивающая объём image,
image — вектор, направленный по нормали от объёма image.
Вектор S


Учитывая то, что конечный объём ограничен набором плоских граней, можно выражение (2.3) преобразовать к сумме интегралов по поверхности:

(2.4) (HJ-3.13)
image,

где image выражает значение переменной в центре грани,
image — вектор площади, выходит из центра грани, направлен в сторону от ячейки (локально), в сторону от ячейки с меньшим индексом к ячейке с большим индексом (глобально).

Еще немного про вектор S


Чтобы не хранить одни и те же параметры вектора image два раза, т.к. очевидно, что у двух соседних ячеек вектор-нормали к грани между ячейками, направленный в сторону От центра ячейки будет различаться только направлением-знаком. Поэтому было создано owner-neighbour отношение между гранью и ячейкой. Если вектор площади image(глобальный, положительное направление от ячейки с меньшим индексом к ячейке с большим индексом) указывает ОТ центра image ячейки такое отношение между ячейкой и вектором image, а точнее между ячейкой и гранью, обозначается owner). В случае если этот вектор image указывает внутрь рассматриваемой ячейки, то отношение neighbour. Направление влияет на знак величины (+ для owner и — для neighbour) и это важно при суммировании см. далее.



(HJ-3.16)
image

Про разностные схемы


Значение image в центре грани вычисляется через значения image в центрах прилегающих ячеек — способ такого выражения носит название разностной схемы. В 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)
image

Согласно методу конечных объёмов проводится дискретизация по времени и выражение (3) записывается как:

(4)
image

Проинтегрируем (4):

(4.1)
image

Разделим левую и правую часть на image:

(5)
image

Данные для матрицы дискретизации


Теперь мы можем получить систему линейных уравнений для каждого конечного объёма image.

Ниже представлена нумерация узлов сетки, которую мы будем использовать.


Координаты узлов хранятся в /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), описывающее поведение величины

image

будет преобразовано в систему линейных алгебраических уравнений, каждое из которых вида:



или, согласно индексам точек на гранях



А еще все потоки в/из ячейки могут быть выражены в виде суммы



где, например, — коэффициент линеаризации потока в точке-центре ячейки E,
— коэффициент линеаризации потока в точке-центре грани,
— нелинейная часть (например, константа).

Согласно нумерации граней выражение примет вид:


С учетом граничных условий для элемента P_0 линейное алгебраическое уравнение может быть представлено в виде



… подставим ранее полученные коэффициенты…



Поток из inlet'a направлен в ячейку, поэтому имеет отрицательный знак.



Так как у нас в управляющем выражении присутствует кроме диффузионного еще и временной член, но конечное уравнение выглядит как



image

Для P = 1.

image

Для P = 4.

image

Система линейных алгебраических уравнений (СЛАУ) может быть представлена в матричном виде как

image,

где

=== 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.
И в итоге получается вектор значений image

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);


на основе которого получаются значения для вектора image

image

image

image

Затем вектор image подставляется в СЛАУ image и происходит новая итерация расчёта вектора image.

И так до тех пор, пока невязка не достигнет требуемых пределов.

Ссылки


* Некоторые уравнения в этой статье взяты из диссертации Ясака Хрвое (HJ — номер уравнения) и если кому-то захочется прочитать про них подробнее (http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdf)

Скачать файлы задачи можно здесь:
github.com/j-avdeev/case
Файлы решателя:
github.com/j-avdeev/matrHyper1Foam

В качестве бонуса — видео, как распространяется концентрация image.

Tags:
Hubs:
+6
Comments 2
Comments Comments 2

Articles