path: root/src/report/report.tex

\documentclass[a4paper,14pt]{extreport}

\usepackage{./include/bsumain}
\usepackage{./include/bsudiplomatitle14}
\usepackage[utf8x]{inputenc}
\usepackage[hidelinks]{hyperref}
\usepackage{listings}
\usepackage{amsmath}
\graphicspath{ {./src/report/images/} }

\author{Соколов Евгений Викторович}
\title{Численное моделирование течения в плоском канале с обратным уступом}
% \supervisor{Репников Василий Иванович}
\mentor{Тетерев Александр Владимирович}
\subfaculty{Кафедра вычислительной математики}

\begin{document}

\maketitle
{
  \renewcommand{\contentsname}{Содержание}
  \tableofcontents
}

\chapter{Введение}
Моделирование поведения несжимаемой вязкой жидкости - задача, имеющая большое значение в прикладной математике. В данной работе рассматривается движение несжимаемой вязкой жидкости в плоском канале с обратным уступом. Это течение интересно тем, что его можно визуально разбить на некоторые зоны - отрыв потока на кромке уступа, последующее его присоединение вниз по течению к нижней стенке канала, а также зона рециркуляции сразу за уступом. У этой задачи простая геометрия, и, как будет показано далее, решение её зависит лишь от двух параметров: числа Рейнольдса $Re$ и параметра расширения потока $E_R$ (отношения полной высоты канала к высоте входного участка).

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

В силу вышеуказанных причин рассматриваемая задача является практически идеальной для тестирования и демонстрации возможностей новых вычислительных технологий, разрабатываемых для моделирования отрывных течений в полуоткрытых областях. Соответственно в литературе можно найти большое количество работ, посвященных изучению несжимаемых течений вязкой жидкости в плоском канале с обратным уступом в двумерном приближении. Большинство из них выполнено для чисел Рейнольдса, не превышающих $Re = 1000$ (здесь и далее число Рейнольдса строится по полной высоте канала и среднемассовой входной скорости потока). Однако существует несколько публикаций, в которых стационарное решение задачи найдено для более высоких значений $Re$ , вплоть до 3000. И совсем немного статей посвящены задачам при значениях числа Рейнольдса меньше единицы вплоть до $Re = 10^{−4}$. Основной их целью было получение вихрей Моффатта у основания уступа, поскольку при столь низких значениях $Re$ эффект инерции потока пренебрежимо мал по сравнению с вязкими силами и, следовательно, в этом месте движение жидкости имеет чисто стоксовский характер.

\chapter{Постановка задачи}
\section{Схема течения}
Рассматривается движение несжимаемой вязкой жидкости с постоянными характеристиками, такими как плотность $\rho$ и динамический коэффициент вязкости $\mu$. Считается, что поток двумерный и стационарный.

Обозначим высоту канала за $H$ - это будет характерная длина в задаче. $l_c$ и $h_c$ — соответственно безразмерные длина и высота уступа, $L$ — безразмерная длина канала, $(1 − h_c)$ — безразмерная высота входного участка канала. По определению параметр расширения потока $E_R = (1 − h_c)$.


\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{flow-scheme}
    \caption{Схема течения}
\end{figure}


\section{Математическая постановка}
Движение жидкости описывается системой нестационарных двумерных уравнений Навье-Стокса:

\begin{equation}
  \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0,
\end{equation}
\begin{equation} \label{eq:ns2}
  \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = - \frac{\partial p}{\partial x} + \frac{1}{Re}(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}),
\end{equation}
\begin{equation} \label{eq:ns3}
  \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = - \frac{\partial p}{\partial y} + \frac{1}{Re}(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}),
\end{equation}

Здесь: $u$ и $v$ — горизонтальная и вертикальная компоненты скорости соответственно; $p$ — давление; $Re$ — число Рейнольдса $Re = \rho U_a H/\mu$, где $U_a$ — среднемассовая скорость потока на входе в канал. Время и давление обезразмерены на комплексы $H/U_a$ и $\rho U_a ^ 2$ соответственно. Нестационарная форма записи уравнений \ref{eq:ns2} и \ref{eq:ns3} обусловлена использованием метода установления при построении стационарного решения задачи.

Будем считать, что в начальный момент времени жидкость покоится: $u = v = 0$. На левой границе $В4$ течение имеет постоянную горизонтальную скорость $u=u_0, v=0$.


На стенках канала (границы В1, В2, В3, В5) используются граничные условия прилипания и непротекания. На выходе из канала $В6$ граничное условие имеет вид:

\begin{equation} \label{eq:boundary_condition}
  \frac{\partial u}{\partial x} = \frac{\partial v}{\partial y} = 0.
\end{equation}

Здесь необходимо иметь в виду, что в силу отсутствия точной информации о поведении решения на открытой границе В6 условие \ref{eq:boundary_condition} следует рассматривать как приближенное. Это означает, что в процессе итерационного построения решения задачи, исходя из требования баланса массы во всей расчетной области, значения компонент u и v на каждой итерации подвергаются малой корректировке.

\chapter{Метод решения}
\section{Основная сложность}
Основной сложностью вычисления поля скорости является тот факт, что поле давления неизвестно. Уравнение движения будет содержать разницу между двумя \textbf{чередующимися} узлами сетки, а не \textbf{соседними}. Следствием этого является то что давление берётся на сетке, отличающейся от применённой, из-за чего точность решения может снижаться.

Но есть более серьезная проблема - рассмотрим распределение давления на рисунке. Такое поле давления сложно назвать реалистичным, однако в любой точке $P$ соответствующая разница $P_W - P_E$ обращается в нуль. Сокрушительным последствием этого будет то, что такое поле давления будет \textbf{казаться} однородным, хотя таковым и не является.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{zigzag-pressure-field}
    \caption{Зигзагообразное поле давления}
\end{figure}

Разница ещё более заметна в двумерном случае. Например, абсолютно случайно выбранное поле давления, распределенное определенным образом, не будет применять давление ни в направлении $x$, ни $y$ оси. И только вдруг такое поле появится в течение итерационного процесса, ничего это уже не остановит вплоть до того, как процесс сойдётся.

Если взять гладкое поле давления, являющееся решением, то можно построить любое количество таких же "зигзагообразных" полей, также являющихся решениями, т.к уравнения движения в этом случае не будут затронуты. Численный метод, допускающий такие абсурдные решения нам не подходит, поэтому нужно искать что-то другое.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{zigzag-pressure-field-2d}
    \caption{Двумерное зигзагообразное поле давления}
\end{figure}

\clearpage

\section{Смещённая шахматная сетка}
	Для реализации алгоритма будет используется смещённая шахматная сетка. Она решает вышеописанные проблемы следующим образом: значения разных переменных рассчитываются не обязательно в одних и тех же узловых точках. Дело в том, что на самом деле мы можем использовать отдельную сетку хоть на каждую независимую переменную. Именно это мы и делаем со скоростью - распределяя её в узловых точках, смещенных относительно узловых точек давления, мы избавляемся от проблем, описанных ранее.
  При расположении в шахматном порядке сетке составляющие скорости рассчитываются для точек, лежащих на гранях контрольных объемов. Таким образом, составляющая скорости $u$ вдоль оси $x$ рассчитывается на гранях, перпендикулярных направлению оси $x$.

Точки, в которых рассчитывается горизонтальная компонента скорости $u$ смещены в направлении оси $x$. Аналогично, значения $v$ рассчитываются в точках, смещенных в направлении оси $y$. Ниже показана двумерная сетка, где узловые точки для $u$ и $v$ помещены на соответствующих гранях контрольных объёмов.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{staggered-grid-2}
    \caption{Смещённая шахматная сетка}
\end{figure}

Шахматная сетка обладает рядом дополнительных преимуществ:
\begin{itemize}
  \item{Нам необязательно интерполировать значения скорости для вычисления других переменных}
  \item{Разность давлений между соседними узлами определяет скорость в точке между этими узлами}
  \item{Исключив неестественные зигзагообразные решения, мы гарантируем то, что только физические поля скорости будут удовлетворять уравнению неразрывности}
\end{itemize}

\section{Уравнения количества движения}
На рисунке \ref{fig:cvu} представлен контрольный объём для горизонтальной компоненты скорости $u$. Если не брать в расчёт наличие других переменных, он не представляет собой ничего необычного. Однако этот контрольный объём на самом деле \textbf{смещён} относительно контрольного объёма для узла основной сетки $P$. Смещение происходит только в направлении $x$, так что грани, перпендикулярные движению проходят через точки $P$ и $E$. Такое расположение реализует главное преимущество смещённой шахматной сетки - разница $P_P - P_E$ может использоваться для вычисления силы давления, действующей на контрольный объём для $u$.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{control-volume-u}
    \caption{Контрольный объём для $u$}
    \label{fig:cvu}
\end{figure}

Для расчёта коэффициента диффузии и массового расхода на гранях контрольного объема используется дискретная формула, для $u_e$ она выглядит как:

\begin{equation} \label{eq:momentum1}
  a_e u_e = \sum a_{nb} u_{nb} + b + (p_p - p_E)A_{e^*}.
\end{equation}

Здесь и в дальнейшем индексы $n, e, s, w$ (\emph{North, East, South, West}) обозначают элементы, смещённые в северном, восточном, южном и западном направлении соответственно. Пересчёт всех соседних элементов обозначен индексом $nb$ (\emph{neighbour}).

В уравнении \ref{eq:momentum1} коэффициенты $a_{nb}$ учитывают влияние конвекции и диффузии на гранях контрольного объёма. Слагаемое $b$ представляет собой массовый расход в контрольном объёме. Градиент давления порождает последнее слагаемое - это сила давления действующая на контрольный объём, где $A_e$ - площадь на которую действует это самое давление. В двумерном случае $A_e$ это просто $\Delta y \times 1$.

Здесь и в дальнейшем индексы $n, e, s, w$ (\emph{North, East, South, West}) обозначают элементы, смещённые в северном, восточном, южном и западном направлении соответственно. Индекс $nb$ (\emph{neighbour}) обозначает каждый из таких соседних элементов.

Уравнения количества движения можно решить только имея заданное поле давления $p$. Если при решении использовать неверное поле давления, найденное поле скорости не будет удовлетворять уравнению неразрывности. Такое поле скорости, рассчитанное с помощью угаданного $p^*$, обозначим как $u^*$ и $v^*$.

Это поле скорости находится в результате решения следующих уравнений:

\begin{equation} \label{eq:momentum2}
  a_e u_e^* = \sum a_{nb} u_{nb}^* + b + (p_p^* - p_E^*)A_{e^*};
\end{equation}
\begin{equation}
  a_n u_n^* = \sum a_{nb} u_{nb}^* + b + (p_p^* - p_N^*)A_{n^*};
\end{equation}

\section{Поправки скорости и давления}
Будем находить улучшенное приближение для $p^*$. Мы хотим, чтобы с каждой итерацией поле скорости всё лучше удовлетворяло уравнению неразрывности. Предположим, что истинное давление находится из соотношения:

\begin{equation}
  p = p^* + p',
\end{equation}

где $p’$ назовём поправкой давления. Аналогичным образом введём соотвествующие поправки составляющих скорости:

\begin{equation}
  u = u^* + u'; v = v^* + v';
\end{equation}

Вычитая \ref{eq:momentum2} из \ref{eq:momentum1} и отбрасывая добавочный член, получаем:
\begin{equation}
  a_e u_e' = (p_P' -p_E')A_e
\end{equation}

или

\begin{equation}
  u_e' = d_e(p_P' -p_E')
\end{equation}

где

\begin{equation}
  d_e \equiv A_e/a_{e^*}
\end{equation}

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

\begin{equation}
  u_e = u_e^* + d_e(p_P' - p_E')
\end{equation}
\begin{equation}
  v_n = v_n^* + d_n(p_P' - p_N')
\end{equation}


\section{Уравнение для поправки давления}
Запишем уравнение неразрывности:
\begin{equation} \label{eq:pcorrect}
  \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} = 0.
\end{equation}

Предположим, что плотность непосредственно зависит от давления. Тогда преобразуем его в уравнение для поправки давления, проинтегрировав по заштрихованному контрольному объему:
\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{control-volume}
    \caption{Контрольный объем}
\end{figure}

Предположим, что значение плотности во всем контрольном объеме равно $\rho_P$. Также будем считать, что значение массовой скорости на всей грани контрольного объема определяется значение составляющей скорости $u_e$ в точке $e$. Тогда проинтегрированная форма уравнения неразрывности \ref{eq:pcorrect} становится:

\begin{equation}
\begin{split}
  \frac{(\rho_P - \rho_P^0)\Delta x \Delta y}{\Delta t} +
  [(\rho u)_e - (\rho u)_w]\Delta y +
  [(\rho v)_n - (\rho v)_s]\Delta x
\end{split}
\end{equation}

Теперь, подставив значения составляющих скорости из поправочных формул и сгруппировав соответствующие члены, получим получим уравнение для $p’$:

\begin{equation}
  a_P p_P' = a_E p_E ' + a_W p_W ' + a_N p_N' + a_S p_S ' + b,
\end{equation}

где

\[
\begin{split}
  \begin{cases}
  a_E = {\rho}_e d_e \Delta y; \\
  a_W = {\rho}_w d_w \Delta y; \\
  a_N = {\rho}_n d_n \Delta x; \\
  a_S = {\rho}_s d_s \Delta x; \\
  a_P = a_E + a_W + a_N + a_S; \\
  \end{cases}
\end{split}
\]

\[
\begin{split}
  b = \frac{(\rho_P - \rho_P^0)\Delta x \Delta y}{\Delta t} +
  [(\rho u^*)_e - (\rho u^*)_w]\Delta y + [(\rho v^*)_n - (\rho v^*)_s]\Delta x
\end{split}
\]


\section{Алгоритм SIMPLE}
\textbf{SIMPLE} (Semi-Implicit Method for Pressure-Linked Equations) — полунеявный метод для связанных давлением уравнений.
Последовательность операций:

\begin{enumerate}
  \item{Задание поля давления $p^*$}
  \item{Решение уравнений движения для получения $u^*$, $v^*$}
  \item{Решение уравнения для $p'$}
  \item{Расчет $p$ путем добавления $p'$ к $p^*$}
  \item{Расчет $u, v$ с учетом соответствующих значений со звездочкой и с помощью формул для поправки скорости}
  \item{Решение дискретных аналогов для других Ф (таких как температура, концентрация и турбулентные характеристики), если они влияют на поле течение через физические свойства жидкости}
  \item{Представление скорректированного давления $p$ как нового $p^*$, возвращение к пункту 2 и повторение всей процедуры до тех пор, пока не будет получено сходящееся решение.}
\end{enumerate}

Блок-схема:
\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{scheme}
    \caption{Блок-схема алгоритма SIMPLE}
\end{figure}

\clearpage

\section{Детальный разбор алгоритма}
Метод называется полунеявным, т.к мы отбросили слагаемое $\sum a_{nb} u_{nb}'$. Оно представляет собой неявное влияние поправки давления на скорость; поправки давления в соседних узлах могут изменять соседние скорости, тем самым также создавая поправку к скорости к рассматриваемой точке.

Отбрасывание этого слагаемого было бы непозволительно, если бы оно приводило к получению неточных решений уравнение движения и неразрывности. Однако так получается, что конечное решение итерационного процесса не содержит неточностей, вызванных отбрасыванием этого члена. Поле скорости, полученное таким образом, всё ещё удовлетворяет уравнению неразрывности. Отсюда делаем вывод, что детали построения поправки давления $p'$ не влияют на корректность конечного решения.

Осталось понять, что будет индикатором сходимости такого итерационного процесса. Рассмотрим последнюю итерацию: используя предыдущие результаты, будем как обычно вычислять массовый расход $b$ для уравнения поправки давления. Однако он будет практически нулевым, чтобы выполнялось $p' = 0$ по всей сетке, что будет значить корректность $u^*, v^*$ как конечных значений компонент скорости. Поэтому тот факт, что $b$ устремляется в ноль говорит нам о полученном решении. И конечно, это решение никак не содержит отклонений вызванных аппроксимациями, использованными для $p'$, т.к мы по факту не используем $p'$ в последней итерации.

Таким образом, массовый расход $b$ является полезным индикатором сходимости. Итерации следует продолжать, пока значение $b$ на всей сетке не станет достаточно малым.

Таким образом, можно сделать вывод, что алгоритм нахождения $p'$ - не более чем промежуточная формула. Выбор такой формулы не отклоняет нас от промежуточного решения (однако может влиять на сходимость!). Но если итерационный процесс всё-таки сходится, то любые формулы для $p'$ будут приводить нас к одинаковому правильному решению.

Скорость сходимости, однако, зависит от этого алгоритма. Если мы отбросим слишком много слагаемых, то можем и вообще потерять сходимость.

Алгоритм нахождения $p'$, используемый нами, также склонен приводить к расхождению процесса, если не использовать релаксацию. Достаточно будет задать коэффициент релаксации для $u^*, v^*$ около $0.5$. Для давления будем использовать $\alpha_p = 0.8$:

\begin{equation}
  p = p^* + \alpha_p \cdot p'
\end{equation}


\chapter{Результаты}
\section{Течение в пустом канале}
В пустом канале течение с постоянной горизонтальной скоростью на входе приобретает параболический профиль:
\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/empty/0.001/raw}
    \includegraphics[width=\textwidth]{results/empty/0.001/streamplot}
    \caption{Картина течения в пустом канале}
\end{figure}


\section{Сравнение результатов для разных чисел $Re$}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re300/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re300/0.001/streamplot}
    \caption{Картина течения при $Re=300$}
\end{figure}

Здесь хорошо наблюдается зона рециркуляции за ступенькой - явно виден вихрь прямо за ступенькой. Примерно на $x = 0.35$ поток соединяется с нижней стенкой и вновь продолжает горизонтальное движение по каналу. Горизонтальные размеры вихря примерно $0.01$.

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re700/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re700/0.001/streamplot}
    \caption{Картина течения при $Re=700$}
\end{figure}

\clearpage
Как мы можем видеть, увеличение числа $Re$ приводит к более выпуклому вихревому следу за ступенькой. Здесь горизонтальный размер вихря уже превышает $0.015$. Легко предположить, что его уменьшение приведёт к обратному эффекту - это и наблюдается:

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re10/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re10/0.001/streamplot}
    \caption{Картина течения при $Re=10$}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re0.1/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re0.1/0.001/streamplot}
    \caption{Картина течения при $Re=0.1$}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re0.0001/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re0.0001/0.001/streamplot}
    \caption{Картина течения при $Re=10^{-4}$}
\end{figure}

Как мы можем видеть, малые значения $Re$ дают меньший вихрь, однако зависимость явно не линейная. Размеры зоны рециркуляции при $Re = 10$ и $Re = 10^{-4}$ практически не отличаются.

\clearpage
А вот примеры этого же канала с ещё большим значением $Re = 1000$:

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re1000/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.02x0.02_Re1000/0.001/streamplot}
    \caption{Картина течения при $Re=1000$}
\end{figure}

Здесь зона присоединения начинается аж с $x = 0.045$.

\section{Сравнение результатов для разных размеров ступеньки}
\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.01x0.03_Re300/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.01x0.03_Re300/0.001/streamplot}
    \caption{Картина течения с размерами ступеньки $0.01\times0.03$}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.04x0.01_Re300/0.001/plot}
    \includegraphics[width=\textwidth]{results/0.05x0.05_0.04x0.01_Re300/0.001/streamplot}
    \caption{Картина течения с размерами ступеньки $0.04\times0.01$}
\end{figure}

Как и ожидалось, геометрия канала, безусловно влияет на картину течения. Ширина ступеньки здесь, на самом деле, не столь важна, как её высота. Самое главное, чтобы ширины ступеньки было достаточно для того, чтобы течение над ступенькой успело приобрести параболический профиль (в случае если мы задаём постоянную скорость на входе). При необходимости, такое течение можно сразу подать на вход, слегка изменив начальные условия. Отсюда можно сделать вывод, что параметр расширения потока $E_R$ является вторым определяющим фактором для течения, помимо числа Рейнольдса $Re$.

\chapter{Выводы}
В результате анализа результатов были сделаны следующие выводы:
\begin{enumerate}
  \item{На выходе из канала течение приобретает параболический профиль}
  \item{Течение в канале зависит от конфигурации ступеньки}
  \item{Высота ступеньки (и соответственно параметр расширения потока $E_R$) играет гораздо большую роль, чем её ширина}
  \item{За ступенькой образуется вихрь}
  \item{Размер и выпуклость вихря зависят также от числа Рейнольдса}
  \item{Бoльшие числа Рейнольдса дают более выпуклый вихрь и наоборот}
  \item{Линейные размеры зоны рециркуляции находятся в нелинейной зависимости от числа Рейнольдса}
\end{enumerate}

\chapter{Заключение}
В результате дипломной работы было смоделировано течение жидкости в плоском канале с обратным уступом. Был проведён анализ полученных результатов и сделаны соответствующие выводы. Было подтверждено, что определяющими факторами такой задачи являются параметр расширения потока $E_R$ и число Рейнольдса $Re$.

\chapter{Код программы}
\lstinputlisting[language=Python, caption=Исходный код, label=lst:code]{./src/simple.py}
\lstinputlisting[language=Python, caption=Исходный код (функции для исследования), label=lst:code]{./src/research.py}
\lstinputlisting[language=Python, caption=Исходный код (функции отрисовки), label=lst:code]{./src/plotter.py}

\chapter{Список использованных источников}
\begin{enumerate}
  \item{Patankar SV. Numerical Heat Transfer and Fluid Flow, 1984}
  \item{Versteeg HK, Malalasekera W. An Introduction to Computational Fluid Dynamics, 1995}
  \item{Devic I. Fluid simulation with SIMPLE method using graphic processors, 2014}
  \item{Ferziger JH, Peric M. Computational Methods for Fluid Dynamics, 2002}
  \item{Андерсон Д, Таннехилл Дж, Плетчер Р. Вычислительная гидромеханика и теплообмен, 1990}
\end{enumerate}

\end{document}