[Из песочницы] OpenFOAM для чайников

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

Для тех, кому это тоже интересно — эта статья. Те, кто хорошо знаком с OpenFOAM или с методами в нём реализованными — пишите о найденных ошибках, если таковые имеются — лучше в личку.

В данной статье я на простом примере-задаче опишу операции, которые происходят в ходе расчёта в OpenFOAM.

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

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

Построение сетки — делим куб на 5 равных частей вдоль оси Z.

Поведение потока задаётся уравнением переноса (1).

(1)

,

где выражает концентрацию некоторого вещества, а выражает перенос вещества, поток.

Если читатель не знаком или забыл оператор дивергенции

Разбиение по пространству

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

(2)

,

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

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

(2.1) (HJ-3.12)*

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

,

где

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

,

где замкнутая поверхность, ограничивающая объём ,

— вектор, по величине равный площади бесконечномалого участка на грани, направлен вектор по нормали от объёма .

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

(2.4) (HJ-3.13)

,

где выражает значение переменной в центре грани,

— вектор площади, выходит из центра грани, направлен от центра ячейки.

Вектор площади указывает от центра , если мы его рассматриваем относительно нашей ячейки (в OpenFOAM такое отношение между ячейкой и вектором , а точнее между ячейкой и гранью, обозначается owner). Для соседних ячеек этот вектор указывает внутрь (отношение neighbour). Направление влияет на знак величины и это важно при суммировании.

(HJ-3.16)

Численно это площадь одной грани. А вот значение выражается через значения в других ячейках — способ такого выражения носит название разностной схемы. В нашей задаче в файле fvSchemes тип разностной схемы задан следующим образом:

divSchemes
{
    default         none;
    div(phi,psi) Gauss linear;
}

Gauss — означает, что выбрана центральная разностная схема;

linear — означает, что интерполяция с центров ячеек на центры поверхностей будет происходить линейно.

Значение в центре грани будет вычисляться по следующей схеме:

(HJ-3.19)

,

Учитывая (2.1) и (2.4) выражение (2) принимает вид:

(3)

Дискретизация по времени

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

(4)

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

(4.1)

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

(5)

Формирование матрицы дискретизации

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

Для P = 0 выражение (5) примет вид:

Для P=1:

Для P=4:

Так как на выходе задано граничное условие

outlet
{
    type                zeroGradient;
}

то получаем

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

,

где

=== 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://ift.tt/PSJmBe)

Скачать файлы задачи можно здесь:

http://ift.tt/PSJmBj

Файлы решателя:

http://ift.tt/PSJmBl

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

This entry passed through the Full-Text RSS service — if this is your content and you're reading it on someone else's site, please read the FAQ at http://ift.tt/jcXqJW.