Метод конечных объемов. Метод конечных объёмов Свойства дискретных схем

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

На хабре уже была пара статей про OpenFOAM:

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

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

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

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

(1)
,

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

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

Оператор дивергенции на самом деле это

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

Где i, j, k - единичные векторы.
Если скалярно умножить оператор набла на векторную величину, то мы получим дивергенцию данного вектора:

«С точки зрения физики, дивергенция векторного поля является показателем того, в какой степени данная точка пространства является источником или стоком этого поля»

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

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

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

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

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

Много формул

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

(2)
,

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

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

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

(2.1) (HJ-3.12)*

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

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

(2.3)
,

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

Вектор S

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

(2.4) (HJ-3.13)
,

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

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

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

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

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

Psi = dimensions ; internalField nonuniform List 5(0.0246875 0.000308546 3.85622e-06 4.81954e-08 5.95005e-10);

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

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

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

Ссылки

* Некоторые уравнения в этой статье взяты из диссертации Ясака Хрвое (HJ - номер уравнения) и если кому-то захочется прочитать про них подробнее (

Ранее поминался метод подобластей, послуживший отправной точкой для ряда численных методов. Одним из таких методов является метод конечных объемов. Этот же метод является представителем еще одного широко распространившегося класса – интегральных методов. От классической формы записи метода подобластей взято разбиение расчетной области на подобласти и интегрирование невязки по подобласти. Отличием является отсутствие явной записи аппроксимирующей (пробной) функции. Но, по-прежнему, пытаемся «точно» решить уравнение в каждой подобласти. Поэтому по подобласти интегрируется исходное уравнение. Интегральные методы характеризуются тем, что сначала берется интеграл от дифференциального уравнения, получается интегральная форма записи уравнения. Затем уравнение в этой форме применяют к отдельным ячейкам сетки. В данном случае ячейки и подобласти – это одно и то же.

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

Как правило, интегралы от уравнений имеют простой и понятный физический смысл. Для примера рассмотрим уравнение неразрывности. Исходное дифференциальное уравнение записывается

проинтегрируем его по объему V, имеющему поверхность S, и по времени в интервале от t 0 до t 1 . При интегрировании производных воспользуемся формулой Стокса (частные случаи ее носят названия формул Грина и Остроградского-Гаусса). В результате получаем

В этой записи разность между первыми двумя интегралами означает изменение массы в заданном объеме за рассматриваемый интервал времени. А двойной интеграл показывает массу, втекающую в данный объем через ограничивающую его поверхность за тот же промежуток времени. Естественно, раз речь идет о численных методах, то эти интегралы считаются приближенно. И здесь начинаются вопросы аппроксимации, аналогичные тем, что рассматривались в методе конечных разностей.

Рассмотрим один из простейших случаев – двумерная прямоугольная равномерная сетка. В методе конечных объемов обычно значения функций определяются не в узлах сетки, а в центрах ячеек. Индексируются, соответственно, тоже не линии сетки в каждом направлении, а слои ячеек (см. рис.).

j-1

j+1

k-1

k+1

Для данного случая интегральная форма уравнения запишется так

Как видим, в данном случае мы получили обычное уравнение, какое могли написать и с помощью метода конечных разностей. Значит, к нему можно применять и те же методы исследования устойчивости. (Вопрос «на засыпку»: а устойчива ли данная схема?)

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

j-1

j+1

k-1

k+1

Например, для случая показанного на рисунке интегральная форма уравнения будет иметь вид

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

Для закрепления распишем, как будет выглядеть применение метода конечных объемов к одному из уравнений сохранения импульса. Возьмем плоский стационарный случай для однозарядных ионов. Пренебрежем вязкостью и упругими столкновениями. Получаем уравнение

Для прямоугольной сетки (см рис. выше) получаем

Простейшая аппроксимация такого уравнения запишется так

после сокращений получаем формулу

алгоритм программа моделирование

Отправной точкой метода конечных объёмов (МКО) является интегральная формулировка законов сохранения массы, импульса, энергии и др. Балансовые соотношения записываются для небольшого контрольного объема; их дискретный аналог получается суммированием по всем граням выделенного объема потоков массы, импульса и т.д., вычисленных по каким - либо квадратурным формулам. Поскольку интегральная формулировка законов сохранения не накладывает ограничений на форму контрольного объема, МКО пригоден для дискретизации уравнений гидрогазодинамики как на структурированных, так и на неструктурированных сетках с различной формой ячеек, что, в принципе, полностью решает проблему сложной геометрии расчетной области.

Следует заметить, однако, что использование неструктурированных сеток является довольно сложным в алгоритмическом отношении, трудоемким при реализации и ресурсоемким при проведении расчетов, в особенности при решении трехмерных задач. Это связано как с многообразием возможных форм ячеек расчетной сетки, так и с необходимостью применения более сложных методов для решения системы алгебраических уравнений, не имеющей определенной структуры. Практика последних лет показывает, что развитые разработки вычислительных средств, базирующихся на использовании неструктурированных сеток, по силам лишь достаточно крупным компаниям, имеющим соответствующие людские и финансовые ресурсы. Гораздо более экономичным оказывается использование блочно-структурированных сеток, предполагающее разбиение области течения на несколько подобластей (блоков) относительно простой формы, в каждой из которых строится своя расчетная сетка. В целом такая составная сетка не является структурированной, однако внутри каждого блока сохраняется обычная индексная нумерация узлов, что позволяет использовать эффективные алгоритмы, разработанные для структурированных сеток. Фактически, для перехода от одноблочной сетки к многоблочной необходимо лишь организовать стыковку блоков, т.е. обмен данными между соприкасающимися подобластями для учета их взаимного влияния. Заметим также, что разбиение задачи на отдельные относительно независимые блоки естественным образом вписывается в концепцию параллельных вычислений на кластерных системах с обработкой отдельных блоков на разных процессорах (компьютерах). Все это делает использование блочно-структурированных сеток в сочетании с МКО сравнительно простым, но чрезвычайно эффективным средством расширения геометрии решаемых задач, что исключительно важно для небольших университетских групп, разрабатывающих собственные программы в области гидрогазодинамики.

Отмеченные выше достоинства МКО послужили основанием к тому, что в начале 1990-х гг. именно этот подход с ориентацией на использование блочно-структурированных сеток был выбран авторами в качестве основы для разработки собственного пакета программ широкого профиля для задач гидрогазодинамики и конвективного теплообмена.

Описание

Неформальное

Выбирается некоторая замкнутая область течения жидкости или газа, для которой производится поиск полей макроскопических величин (например, скорости, давления), описывающих состояние среды во времени и удовлетворяющих определенным законам, сформулированным математически. Наиболее используемыми являются законы сохранения в Эйлеровых переменных.

Для любой величины , в каждой точке пространства, окруженной некоторым замкнутым конечным объемом , в момент времени существует следующая зависимость: общее количество величины в объеме может изменяться за счет следующих факторов:

Другими словами, при формулировке МКО используется физическая интерпретация исследуемой величины. Например, при решении задач переноса тепла используется закон сохранения тепла в каждом контрольном объеме.

Математическое

Модификации

Литература

Патанкар С. В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах = Computation of conduction and Duct Flow Heat Transfer: Пер. с англ. - М.: Издательство МЭИ, 2003. - 312 с.

См. также

Wikimedia Foundation . 2010 .

Метод квадратичного решета
Метод конечных отношений

Смотреть что такое "Метод конечных объёмов" в других словарях:

Метод конечных элементов - Решение методом конечных элементов двухмерной магнитостатической задачи (линии и цвет означают направление и величину магнитной индукции) … Википедия

Computer-aided engineering - CAE (англ. Computer aided engineering) общее название для программ и программных пакетов, предназначенных для решения различных инженерных задач: расчётов, анализа и симуляции физических процессов. Расчётная часть пакетов чаще всего… … Википедия

Вычислительная гидродинамика - Вычислительная гидродинамика (англ. Computational fluid dynamics, CFD) подраздел механики сплошных сред, включающий совокупность физических, математических и численных методов, предназначенных для вычисления характеристик потоковых… … Википедия

Прямое численное моделирование - (англ. DNS (Direct Numerical Simulation)) один из методов численного моделирования течений жидкости или газа. Метод основан на численном решении системы уравнений Навье Стокса и позволяет моделировать в общем случае движение вязких… … Википедия

Matrix Template Library - Тип Математическое ПО Операционная система Linux, Unix, Mac OS X, Windows Языки интерфейса C++ Лицензия Boost Software License … Википедия

МКО - машинно котельное отделение Словарь: С. Фадеев. Словарь сокращений современного русского языка. С. Пб.: Политехника, 1997. 527 с. МКО Межамериканский комитет обороны воен. Словарь: Словарь сокращений и аббревиатур армии и спецслужб. Сост. А. А.… … Словарь сокращений и аббревиатур

Компьютерное моделирование - краш теста методом конечных элементов. Компьютерная модель (англ. computer model), или численная мод … Википедия

Численное моделирование - Компьютерное моделирование является одним из эффективных методов изучения сложных систем. Компьютерные модели проще и удобнее исследовать в силу их возможности проводить т.н. вычислительные эксперименты, в тех случаях когда реальные эксперименты… … Википедия

ГАЗОВАЯ ДИНАМИКА - раздел гидроаэромеханики, в к ром изучается движение сжимаемых сплошных сред (газа, плазмы) и их вз ствие с тв. телами. Как часть физики, Г. д. связана с термодинамикой и акустикой. Св во сжимаемости состоит в способности в ва изменять свой… … Физическая энциклопедия

Механика сплошных сред - изучает движение и равновесие газов, жидкостей и деформируемых твёрдых тел. Моделью реальных тел в М. с. с. является сплошная среда (СС); в такой среде все характеристики вещества являются непрерывными функциями пространственных координат и… … Энциклопедия техники

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

Рис. 13. Расчетная сетка, используемая для решения уравнения (31)

методом конечных объемов

Используя теорему о среднем можно записать

где Δх, Δу – длины граней ячейки, x W – абсцисса левой ("западной") границы ячейки А, x Е – абсцисса правой ("восточной") границы, у N – ордината верхней ("северной") границы, у S – ордината нижней ("южной") границы, S * – средняя по ячейке скорость тепловыделения. Индекс у производных (*), в левой части (32), указывает на то, что их следует рассматривать как средние значения, определенные таким образом, чтобы правильно представить тепловые потоки на каждой из границ. С учетом данного обстоятельства, дискретный аналог (32) может быть получен без затруднений [Патанкар].

Таким образом, уравнение (32) описывает баланс тепла (закон сохранения энергии) в пределах ячейки А. При условии правильного описания тепловых потоков между ячейками, система, составленная из уравнений вида (32), примененных к каждому контрольному объему, будет верно описывать баланс тепла во всей расчетной области.

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

5. Свойства дискретных схем

5.1 Точность

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

Когда говорят о точности дискретной схемы, обычно имеют в виду погрешность аппроксимации производных 27 . В частности, если погрешность аппроксимации сопоставима со второй степенью шага расчетной сетки, то говорят, что дискретная схема имеет второй порядок точности. Более подробно этот вопрос рассматривался в § 3.

5.2 Согласованность

Дискретная схема называется согласованной с исходным дифференциальным уравнением, если при измельчении расчетной сетки погрешность аппроксимации (см. § 3) стремится к нулю,

Известны расчетные схемы, у которых для достижения согласованности необходимо выполнение дополнительных условий, [Андерсон и К]. Поскольку проверка согласованности расчетных схем является задачей разработчиков (а не пользователей) программного обеспечения более подробно этот вопрос здесь обсуждаться не будет.