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 интенсивность течения плавно нарастает до максимального значения по закону

\begin{equation}
  u(t, y) = 6f(t)(y - h_c)(1 - y) / (1 - h_c)^2,
\end{equation}
\begin{equation}
  v = 0,
\end{equation}

где

\[
  f(t) =
  \begin{cases}
    0.5 (sin(0.5\pi(2t / t_m - 1)) + 1), & 0 \leq t \leq t_m, \\
    1, & t_m < t.
  \end{cases}
\]

На стенках канала (границы В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$.

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

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

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

\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^*$, обозначим как $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 p}{\partial t} + \frac{\partial (pu)}{\partial x} + \frac{\partial (pv)}{\partial y} + \frac{\partial (pw)}{\partial z} = 0.
\end{equation}

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

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

\begin{equation}
\begin{split}
  \frac{(p_P -p_P^0)\Delta x \Delta y}{\Delta t} +
  [(pu)_e - (pu)_w]\Delta y +
  [(pv)_n - (pv)_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 = p_e d_e \Delta y; \\
  a_W = p_w d_w \Delta y; \\
  a_N = p_n d_n \Delta x; \\
  a_S = p_s d_s \Delta x; \\
  a_P = a_E + a_W + a_N + a_S; \\
  \end{cases}
\end{split}
\]

\[
\begin{split}
  b = \frac{(p_P -p_P^0)\Delta x \Delta y}{\Delta t} +
  [(pu^*)_e - (pu^*)_w]\Delta y + [(pv^*)_n - (pv^*)_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'$ не оказывают влияния на правильность решения.

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

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

Учитывая всё вышесказанное, уравнение поправки давления можно рассматривать как лишь промежуточный алгоритм, который приводит нас к правильному полю давления, и не отклоняет от точного решения. Если итерационный процесс сходится, все алгоритмы для $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}

\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$ приводит к более выпуклому вихревому следу за ступенькой. Легко предположить, что его уменьшение приведёт к обратному эффекту - это и наблюдается:

\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}

\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}

\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}

\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}