path: root/src/report/report.tex

\documentclass[a4paper,14pt]{extreport}

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

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

\begin{document}

\maketitle
\tableofcontents
\newpage

\section{Постановка задачи}
\subsection{Схема течения}
Рассматривается движение несжимаемой вязкой жидкости с постоянными характеристиками, такими как плотность $\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}

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

\begin{equation}
  \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0,
\end{equation}
\begin{equation}
  \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}
  \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$ соответственно. Нестационарная форма записи уравнений (2) и (3) обусловлена использованием метода установления при построении стационарного решения задачи.  Считается, что в начальный момент времени жидкость покоится: $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}
  \frac{\partial u}{\partial x} = \frac{\partial v}{\partial y} = 0.
\end{equation}

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

\section{Метод решения}
\subsection{Шахматная сетка}
	Для реализации алгоритма используется шахматная сетка. Её особенностью является то, что она не обязательно рассчитывает все переменные в одних и тех же узловых точках (при желании можно использовать для каждой зависимой переменной свою сетку).
	При расположении в шахматном порядке сетке составляющие скорости рассчитываются для точек, лежащих на гранях контрольных объемов. Таким образом, составляющая скорости $u$ вдоль оси $x$ рассчитывается на гранях, перпендикулярных направлению оси $x$.

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

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

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\textwidth]{staggered-grid-2}
    \caption{Сетка и узловые точки}
\end{figure}

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

\subsection{Уравнения количества движения}
Для расчёта коэффицента диффузии и массого расхода на гранях контрольного объема используется дискеретная формула:

\begin{equation}
  a_e u_e = \sum a_nb u_nb + b + (p_p - p_E)A_{e^*}.
\end{equation}

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

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

\begin{equation}
  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}

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

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

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

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

Вычитая (5.8) из (5.6) и отбрасывая добавочный член, получаем:
\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' - y_N')
\end{equation}


\subsection{Уравнение для поправки давления}
Преобразуем уравнение неразрывности в уравнение для поправки давления. Препдположим, что плотность непосредственно зависит от давления. Уравнение неразрывности имеет вид:
\begin{equation}
  \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$.

\begin{equation}
  TODO
\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 ' + a_T p_T ' + a_B p_B ' + b,
\end{equation}

где

\begin{equation}
 TODO
\end{equation}


\subsection{Алгоритм SIMPLE}
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{Результаты}
Все результаты получены при $Re = 400$:

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{u}
    \includegraphics[width=0.5\textwidth]{v}
    \caption{Компоненты скорости после распределения $(u, v)$}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{vector-field}
    \caption{Картина течения на входе в канал}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{whirl}
    \caption{Вихревой след за ступенькой}
\end{figure}

\end{document}