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
{
  \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{Шахматная сетка}
	Для реализации алгоритма используется шахматная сетка. Её особенностью является то, что она не обязательно рассчитывает все переменные в одних и тех же узловых точках (при желании можно использовать для каждой зависимой переменной свою сетку).
	При расположении в шахматном порядке сетке составляющие скорости рассчитываются для точек, лежащих на гранях контрольных объемов. Таким образом, составляющая скорости $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}

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

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

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

Уравнения количества движения можно решить только в том случае, если поле давления задано или каким-то образом найдено. Если при решении использовать неверное поле давления, найденное поле скорости не будет удовлетворять уравнению неразрывности. Выразим такое поле скорости, полученное с использованием приближенного поля давления $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' - y_N')
\end{equation}


\section{Уравнение для поправки давления}
Преобразуем уравнение неразрывности в уравнение для поправки давления. Препдположим, что плотность непосредственно зависит от давления. Уравнение неразрывности имеет вид:
\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}
\begin{split}
  \frac{(p_P -p_P^0)\Delta x \Delta y \Delta z}{\Delta t} +
  [(pu)_e - (pu)_w]\Delta y \Delta z +
  [(pv)_n - (pv)_s]\Delta z \Delta x + \\
  [(pw)_t - (pw)_b]\Delta x \Delta y = 0
\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 ' + a_T p_T ' + a_B p_B ' + b,
\end{equation}

где

\[
\begin{split}
  \begin{cases}
  a_E = p_e d_e \Delta y \Delta z; a_W = p_w d_w \Delta y \Delta z; a_N = p_N d_n \Delta z \Delta x; \\
  a_S = p_s d_s \Delta z \Delta x; a_T = p_t d_t \Delta x \Delta y; a_B = p_b d_b \Delta x \Delta y; \\
  a_P = a_E + a_W + a_N + a_S + a_T + a_B; \\
  \end{cases}
\end{split}
\]

\[
\begin{split}
  b = \frac{(p_P -p_P^0)\Delta x \Delta y \Delta z}{\Delta t} +
  [(pu^*)_e - (pu^*)_w]\Delta y \Delta z + \\
  + [(pv^*)_n - (pv^*)_s]\Delta z \Delta x +
  [(pw^*)_t - (pw^*)_b]\Delta x \Delta y.
\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

\chapter{Результаты}
Все результаты получены при $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}

\chapter{Код программы}
\lstinputlisting[language=C, caption=Исходный код, label=lst:code]{./src/source.cpp}

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