diff options
| -rw-r--r-- | src/report/report.tex | 82 |
1 files changed, 44 insertions, 38 deletions
diff --git a/src/report/report.tex b/src/report/report.tex index ee3cb06..53c23d1 100644 --- a/src/report/report.tex +++ b/src/report/report.tex @@ -23,12 +23,18 @@ } \chapter{Введение} -Исследование отрывных течений имеет важное значение как с фундаментальной, так и прикладной точек зрения. Примером подобного течения является движение несжимаемой вязкой жидкости в плоском канале с обратным уступом, в котором присутствуют как отрыв потока на кромке уступа с его последующим присоединением вниз по течению к нижней стенке канала, так и зона рециркуляции жидкости сразу же за уступом. Причем в зависимости от соотношения инерции течения и вязких сил таких зон рециркуляции и сопровождающих их точек отрыва и присоединения потока к стенкам канала может быть несколько. При этом задача характеризуется простотой геометрии и зависимостью решения, вообще говоря, всего от двух параметров: числа Рейнольдса $Re$ и параметра расширения потока $E_R$, равного отношению полной высоты канала к высоте входного участка. Однако имеет место фактор, осложняющий задачу, — открытая выходная граница канала. Единая методология решения задач с открытыми границами, особенно для несжимаемых течений, до сих пор не создана, что и определяет общепринятую тактику локализации такой границы — чем она дальше от зоны основных возмущений потока, тем лучше. +Моделирование поведения несжимаемой вязкой жидкости - задача, имеющая большое значение в прикладной математике. В данной работе рассматривается движение несжимаемой вязкой жидкости в плоском канале с обратным уступом. Это течение интересно тем, что его можно визуально разбить на некоторые зоны - отрыв потока на кромке уступа, последующее его присоединение вниз по течению к нижней стенке канала, а также зона рециркуляции сразу за уступом. У этой задачи простая геометрия, и, как будет показано далее, решение её зависит лишь от двух параметров: числа Рейнольдса $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)$. +Рассматривается движение несжимаемой вязкой жидкости с постоянными характеристиками, такими как плотность $\rho$ и динамический коэффициент вязкости $\mu$. Считается, что поток двумерный и стационарный. + +Обозначим высоту канала за $H$ - это будет характерная длина в задаче. $l_c$ и $h_c$ — соответственно безразмерные длина и высота уступа, $L$ — безразмерная длина канала, $(1 − h_c)$ — безразмерная высота входного участка канала. По определению параметр расширения потока $E_R = (1 − h_c)$. + \begin{figure}[h] \centering @@ -36,6 +42,7 @@ \caption{Схема течения} \end{figure} + \section{Математическая постановка} Движение жидкости описывается системой нестационарных двумерных уравнений Навье-Стокса: @@ -49,26 +56,12 @@ \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$ и $v$ — горизонтальная и вертикальная компоненты скорости соответственно; $p$ — давление; $Re$ — число Рейнольдса $Re = \rho U_a H/\mu$, где $U_a$ — среднемассовая скорость потока на входе в канал. Время и давление обезразмерены на комплексы $H/U_a$ и $\rho U_a ^ 2$ соответственно. Нестационарная форма записи уравнений \ref{eq:ns2} и \ref{eq:ns3} обусловлена использованием метода установления при построении стационарного решения задачи. -\begin{equation} - u(t, y) = 6f(t)(y - h_c)(1 - y) / (1 - h_c)^2, -\end{equation} -\begin{equation} - v = 0, -\end{equation} +Будем считать, что в начальный момент времени жидкость покоится: $u = v = 0$. На левой границе $В4$ течение имеет постоянную горизонтальную скорость $u=u_0, v=0$. -где -\[ - 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 граничное условие имеет вид: +На стенках канала (границы В1, В2, В3, В5) используются граничные условия прилипания и непротекания. На выходе из канала $В6$ граничное условие имеет вид: \begin{equation} \label{eq:boundary_condition} \frac{\partial u}{\partial x} = \frac{\partial v}{\partial y} = 0. @@ -101,9 +94,10 @@ \clearpage \section{Смещённая шахматная сетка} - Описанные выше трудности могут быть разрешены, если учесть тот факт, что нам не нужно вычислять все переменные в одних и тех же узлах сетки. Если понадобится, можно использовать отдельную сетку хоть на каждую независимую переменную. Это особенно удобно в случае скорости - распределяя её на сетке, отличной от давления, мы избавляемся от проблемы. + Для реализации алгоритма будет используется смещённая шахматная сетка. Она решает вышеописанные проблемы следующим образом: значения разных переменных рассчитываются не обязательно в одних и тех же узловых точках. Дело в том, что на самом деле мы можем использовать отдельную сетку хоть на каждую независимую переменную. Именно это мы и делаем со скоростью - распределяя её в узловых точках, смещенных относительно узловых точек давления, мы избавляемся от проблем, описанных ранее. + При расположении в шахматном порядке сетке составляющие скорости рассчитываются для точек, лежащих на гранях контрольных объемов. Таким образом, составляющая скорости $u$ вдоль оси $x$ рассчитывается на гранях, перпендикулярных направлению оси $x$. - Для реализации алгоритма будет используется смещённая шахматная сетка. Её особенностью является то, что она не обязательно рассчитывает все переменные в одних и тех же узловых точках. При расположении в шахматном порядке сетке составляющие скорости рассчитываются для точек, лежащих на гранях контрольных объемов. Таким образом, составляющая скорости $u$ вдоль оси $x$ рассчитывается на гранях, перпендикулярных направлению оси $x$. +Точки, в которых рассчитывается горизонтальная компонента скорости $u$ смещены в направлении оси $x$. Аналогично, значения $v$ рассчитываются в точках, смещенных в направлении оси $y$. Ниже показана двумерная сетка, где узловые точки для $u$ и $v$ помещены на соответствующих гранях контрольных объёмов. \begin{figure}[h] \centering @@ -111,10 +105,12 @@ \caption{Смещённая шахматная сетка} \end{figure} -Следует отметить, что по отношению к узловым точкам основной сетки точки, в которых определяются $u$, смещены только в направлении оси $x$. - Легко выбрать способ размещения узловых точек для составляющих $v$ и $w$. Ниже показана двумерная сетка, где узловые точки для $u$ и $v$ помещены на соответствующих гранях контрольных объёмов. - -Прямым следствием введения шахматной сетки является то, что массовый расход через грани контрольного объема можно теперь определять без интреполяции соответствующей состовляющей скорости. Также имеются два важных преимущества: при использовании шахматной сетки только физические поля скорости могут удовлетворять уравнению неразрывности. Другое важное преимущество такой сетки заключается в том, что разность давлений между двумя соседними узловыми точками определяет составляющую скорости в точке, расположенной между этими узловыми точками. +Шахматная сетка обладает рядом дополнительных преимуществ: +\begin{itemize} + \item{Нам необязательно интерполировать значения скорости для вычисления других переменных} + \item{Разность давлений между соседними узлами определяет скорость в точке между этими узлами} + \item{Исключив неестественные зигзагообразные решения, мы гарантируем то, что только физические поля скорости будут удовлетворять уравнению неразрывности} +\end{itemize} \section{Уравнения количества движения} На рисунке \ref{fig:cvu} представлен контрольный объём для горизонтальной компоненты скорости $u$. Если не брать в расчёт наличие других переменных, он не представляет собой ничего необычного. Однако этот контрольный объём на самом деле \textbf{смещён} относительно контрольного объёма для узла основной сетки $P$. Смещение происходит только в направлении $x$, так что грани, перпендикулярные движению проходят через точки $P$ и $E$. Такое расположение реализует главное преимущество смещённой шахматной сетки - разница $P_P - P_E$ может использоваться для вычисления силы давления, действующей на контрольный объём для $u$. @@ -126,19 +122,19 @@ \label{fig:cvu} \end{figure} -Для расчёта коэффициента диффузии и массового расхода на гранях контрольного объема используется дискеретная формула, для $u_e$ она выглядит как: +Для расчёта коэффициента диффузии и массового расхода на гранях контрольного объема используется дискретная формула, для $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}) обозначает каждый из таких соседних элементов. +Здесь и в дальнейшем индексы $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^*$. +Уравнения количества движения можно решить только имея заданное поле давления $p$. Если при решении использовать неверное поле давления, найденное поле скорости не будет удовлетворять уравнению неразрывности. Такое поле скорости, рассчитанное с помощью угаданного $p^*$, обозначим как $u^*$ и $v^*$. Это поле скорости находится в результате решения следующих уравнений: @@ -150,7 +146,7 @@ \end{equation} \section{Поправки скорости и давления} -Найдём способ улучшения приближенного поля $p^*$ таким образом, чтобы результирующее поле скорости с каждой итерацией лучше удовлетворяло уравению неразрывности. Предположим, что истинное давление находится из соотношения: +Будем находить улучшенное приближение для $p^*$. Мы хотим, чтобы с каждой итерацией поле скорости всё лучше удовлетворяло уравнению неразрывности. Предположим, что истинное давление находится из соотношения: \begin{equation} p = p^* + p', @@ -190,12 +186,12 @@ \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} @@ -212,7 +208,7 @@ \end{split} \end{equation} -Если теперь вместо всех составляющих скорости подставить выражения из поправочных формул для скорости, то после группировки соответствующих членов получим следующее уравнение для сеточных значений $p’$: +Теперь, подставив значения составляющих скорости из поправочных формул и сгруппировав соответствующие члены, получим получим уравнение для $p’$: \begin{equation} a_P p_P' = a_E p_E ' + a_W p_W ' + a_N p_N' + a_S p_S ' + b, @@ -266,13 +262,13 @@ \section{Детальный разбор алгоритма} Метод называется полунеявным, т.к мы отбросили слагаемое $\sum a_{nb} u_{nb}'$. Оно представляет собой неявное влияние поправки давления на скорость; поправки давления в соседних узлах могут изменять соседние скорости, тем самым также создавая поправку к скорости к рассматриваемой точке. -Отбрасывание этого слагаемого было бы непозволительно, если бы оно приводило к получению неточных решений дискретных форм уравнений движения и неразрывности. Так получается, что решение, к которому сходится такой итерационный процесс, не содержит неточностей, вызванных отбрасыванием этого члена. В решении мы получаем поле давления такое, что поле скорости удовлетворяет уравнению неразрывности. Детали построения $p'$ не оказывают влияния на правильность решения. +Отбрасывание этого слагаемого было бы непозволительно, если бы оно приводило к получению неточных решений уравнение движения и неразрывности. Однако так получается, что конечное решение итерационного процесса не содержит неточностей, вызванных отбрасыванием этого члена. Поле скорости, полученное таким образом, всё ещё удовлетворяет уравнению неразрывности. Отсюда делаем вывод, что детали построения поправки давления $p'$ не влияют на корректность конечного решения. -Рассмотрим теперь последнюю итерацию, после которой мы объявим сходимость. К этому моменту у нас уже есть, в результате всех предыдущих итераций, некоторое поле давления. Используя его как $p^*$, мы решим уравнения движения чтобы получить $u^*, v^*$. Из этого поля скорости, мы вычислим массовый расход $b$ для уравнения поправки давления. Так как это последняя итерация, значение $b$ будет практически нулевым на всех контрольных объёмах. Тогда, $p' = 0$ по всей сетке будет решением уравнения поправки давления, а значения $u^*, v^*$ будут являться правильными значениями компонент скорости. Поэтому тот факт, что $b$ устремляется в ноль говорит нам о полученном решении. И конечно, это решение никак не содержит отклонений вызванных аппроксимациями, использованными для $p'$, т.к мы по факту не используем $p'$ в последней итерации. +Осталось понять, что будет индикатором сходимости такого итерационного процесса. Рассмотрим последнюю итерацию: используя предыдущие результаты, будем как обычно вычислять массовый расход $b$ для уравнения поправки давления. Однако он будет практически нулевым, чтобы выполнялось $p' = 0$ по всей сетке, что будет значить корректность $u^*, v^*$ как конечных значений компонент скорости. Поэтому тот факт, что $b$ устремляется в ноль говорит нам о полученном решении. И конечно, это решение никак не содержит отклонений вызванных аппроксимациями, использованными для $p'$, т.к мы по факту не используем $p'$ в последней итерации. Таким образом, массовый расход $b$ является полезным индикатором сходимости. Итерации следует продолжать, пока значение $b$ на всей сетке не станет достаточно малым. -Учитывая всё вышесказанное, уравнение поправки давления можно рассматривать как лишь промежуточный алгоритм, который приводит нас к правильному полю давления, и не отклоняет от точного решения. Если итерационный процесс сходится, все алгоритмы для $p'$ будут приводить к одинаковому решению. +Таким образом, можно сделать вывод, что алгоритм нахождения $p'$ - не более чем промежуточная формула. Выбор такой формулы не отклоняет нас от промежуточного решения (однако может влиять на сходимость!). Но если итерационный процесс всё-таки сходится, то любые формулы для $p'$ будут приводить нас к одинаковому правильному решению. Скорость сходимости, однако, зависит от этого алгоритма. Если мы отбросим слишком много слагаемых, то можем и вообще потерять сходимость. @@ -303,6 +299,8 @@ \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} @@ -311,7 +309,7 @@ \end{figure} \clearpage -Как мы можем видеть, увеличение числа $Re$ приводит к более выпуклому вихревому следу за ступенькой. Легко предположить, что его уменьшение приведёт к обратному эффекту - это и наблюдается: +Как мы можем видеть, увеличение числа $Re$ приводит к более выпуклому вихревому следу за ступенькой. Здесь горизонтальный размер вихря уже превышает $0.015$. Легко предположить, что его уменьшение приведёт к обратному эффекту - это и наблюдается: \begin{figure}[h] \centering @@ -334,6 +332,8 @@ \caption{Картина течения при $Re=10^{-4}$} \end{figure} +Как мы можем видеть, малые значения $Re$ дают меньший вихрь, однако зависимость явно не линейная. Размеры зоны рециркуляции при $Re = 10$ и $Re = 10^{-4}$ практически не отличаются. + \clearpage А вот примеры этого же канала с ещё большим значением $Re = 1000$: @@ -344,6 +344,8 @@ \caption{Картина течения при $Re=1000$} \end{figure} +Здесь зона присоединения начинается аж с $x = 0.045$. + \section{Сравнение результатов для разных размеров ступеньки} \begin{figure}[h] \centering @@ -359,18 +361,22 @@ \caption{Картина течения с размерами ступеньки $0.04\times0.01$} \end{figure} +Как и ожидалось, геометрия канала, безусловно влияет на картину течения. Ширина ступеньки здесь, на самом деле, не столь важна, как её высота. Самое главное, чтобы ширины ступеньки было достаточно для того, чтобы течение над ступенькой успело приобрести параболический профиль (в случае если мы задаём постоянную скорость на входе). При необходимости, такое течение можно сразу подать на вход, слегка изменив начальные условия. Отсюда можно сделать вывод, что параметр расширения потока $E_R$ является вторым определяющим фактором для течения, помимо числа Рейнольдса $Re$. + \chapter{Выводы} В результате анализа результатов были сделаны следующие выводы: \begin{enumerate} \item{На выходе из канала течение приобретает параболический профиль} - \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} |