Рассмотренные в предыдущих лекциях разностные методы решения уравнений в частных производных демонстрировались на примере либо линейных задач, либо достаточно простых нелинейных уравнений с хорошо изученными свойствами. Такие задачи в реальной вычислительной практике обычно служат тестами для отработки методов решения более сложных нелинейных систем. Традиционным объектом приложения численных методов служат уравнения механики сплошной среды (МСС). Основных причин тому три. Первая — все математические модели МСС известны достаточно давно, хорошо исследованы и являются частью научной классики. Вторая — математические модели МСС нелинейны. Как правило, у линеаризованных уравнений очень узкая область применения. Третья — в результатах решения задач МСС на протяжении всего XX века была практическая заинтересованность, вызванная бурным развитием авиации, осуществлением наукоемких ядерных и космических программ в разных странах.
В данной книге ограничимся самой простой моделью МСС —
В основе построения математических моделей, описывающих поведение жидкостей и газов, лежит понятие о сплошной среде. Из молекулярной физики известно, что среда состоит из отдельных частиц (молекул, ионов, электронов, атомов), расстояние между которым существенно больше их собственных размеров. Длина свободного пробега частицы l (расстояние, пройденное частицей между двумя столкновениями) тем меньше, чем больше частиц заключено в единице объема, чем больше плотность среды. В механике жидкостей и газов рассматриваются среды, содержащие в единице объема большое количество частиц (много больше, чем число Авогадро — число частиц в одной грамм - молекуле вещества, $$N_A \approx 6 \cdot 10^{23} \quad (\mbox{г} \cdot \mbox{моль})^{- 1}$$ ).
В таких средах можно рассматривать лишь некоторые усредненные характеристики, не занимаясь изучением поведения каждой частицы в отдельности. В этом предположении заключается идея модели сплошной среды, непрерывно заполняющей пространство.
Количественным критерием применимости приближения сплошной среды может
служить неравенство $$l/L \ll 1$$, где L — характерный пространственный размер задачи (например, размер тела при внешнем обтекании потоком газа). В газах при нормальных условиях $$l \approx (10^{- 5} \div 10^{- 6})$$
см, поэтому приведенное условие для тел с размером более 1 см выполняется с достаточной точностью. Предложение о сплошности среды, по - видимому, берет свое начало от Эйлера, впервые рассмотревшего газ как непрерывную деформируемую субстанцию. Не вдаваясь в особенности получения
Эйлерова (недивергентная) форма одномерной системы уравнений газодинамики имеет вид
$$\begin{gather*} \frac{{\partial}{\rho}}{{\partial}t} + u \frac{{\partial}{\rho}}{{\partial x}} + {\rho}\frac{{\partial}u}{{\partial}x} = 0 \quad \mbox{(уравнение неразрывности), }\\ \frac{{\partial}u}{{\partial}t} + u \frac{{\partial u}}{{\partial}x} = - \frac{1}{\rho} \cdot \frac{{\partial}p}{{\partial}x} \quad \mbox{(уравнение движения), } \\ \frac{{\partial}e}{{\partial}t} + u \frac{\partial e}{{\partial}x} = - \frac{1}{\rho} \frac{{\partial}(pu)}{{\partial}x} \quad \mbox{(уравнение энергии), } \end{gather*}$$где e — удельная энергия, равная $$$ e = \varepsilon + \frac{u^2}{2}, \varepsilon $$$ удельная внутренняя энергия, u — скорость газа, $$\rho$$ — плотность среды, p — давление, t — время, x — декартова координата. Эту же систему уравнений в частных производных можно представить в матричной (характеристической) форме
где $$\mathbf{U}$$ — вектор - столбец, $$\mathbf{A}$$ — квадратная матрица 3 x 3.
Система замыкается уравнением состояния
$$f(\rho , \varepsilon , p) = 0,$$которое, например, для идеального газа имеет вид
$$\frac{p}{{\rho}(\gamma - 1)} - \varepsilon = 0.$$где $$\gamma$$ — безразмерная постоянная, равная отношению теплоемкости газа при постоянном давлении и теплоемкости при постоянном объеме — постоянная адиабаты.
В системе, записанной в
Не занимаясь выводом формул (это делается простыми алгебраическими преобразованиями), представим другие виды записи уравнения энергии, справедливые для приведенного выше уравнения состояния:
$$\begin{gather*} \frac{{\partial}p}{{\partial}t} + u \frac{{\partial p}}{{\partial}x} + \gamma p \frac{{\partial}u}{{\partial}x} = 0, \\ \frac{{\partial}p}{{\partial}t} + u \frac{{\partial p}}{{\partial}x} + c^2 {\rho}\frac{{\partial}u}{{\partial}x} = 0, \end{gather*}$$где $$c = \sqrt{\gamma {p/{\rho}}}$$ — адиабатическая скорость звука,
$$\begin{gather*} \frac{{{\partial}S}}{{\partial}t} + u \frac{{{\partial}S}}{{\partial}x} = 0, \\ \mbox{или} \quad \frac{dS}{dt} \left|_{tr}\right. = 0, \quad \mbox{где} \quad \frac{d}{dt} \left|_{tr}\right. = \frac{\partial }{{\partial}t} + u \frac{{\partial}}{{\partial}x}, \end{gather*}$$$$S = m(\varepsilon \rho )^{1 - \gamma }$$ — энтропия. Она, как следует из последней формулы, сохраняется вдоль траектории частицы идеального газа, т.е. на траектории уравнения
$$$ \frac{dX}{dt} = u(t, X), X(0) = X_0. $$$или в
Интегральная форма этих уравнений получается при использовании теоремы Гаусса - Остроградского
$$$ \int\int\limits_{\Omega}(\frac{{{\partial}{\mathbf{R}}}}{{\partial}t} + \frac{{{\partial}{\mathbf{Q}}}}{{\partial}x}){dtdx} = \oint\limits_{\Gamma}(\mathbf{R}dx - \mathbf{Q}dt) = 0, $$$где $$\Gamma$$ — граница замкнутой области интегрирования $$\Omega$$ в плоскости t, x. Здесь $$u(t, x), \rho (t, x), p(t, x)$$ — скорость, плотность и давление газа соответственно. В механике сплошных сред вводится эйлерово и лагранжево описание поведения среды. В первом случае наблюдатель полагается неподвижным, например, стоящим на берегу реки. Соответственно расчетная сетка будет неподвижной (фиксированная эйлерова сетка). Во втором случае полагаем, что наблюдатель движется вместе со средой, например, находится на лодке, плывущей по течению реки. В этом случае лагранжева расчетная сетка будет двигаться вместе с частицами среды.
Лагранжева форма одномерных
где $$$ I = \left({\frac{{\partial}x}{{\partial}\xi }}\right)^{-
1}, $$$ x и $$\xi$$ соответственно эйлерова и лагранжева
координаты, связь между которыми дается последним уравнением. Эта система в одномерном случае может быть записана в другом виде. Если ввести лагранжеву массовую координату $$\eta (\xi )$$, связанную с лагранжевой координатой $$\xi$$ дифференциальным уравнением $$$ \frac{{d \eta }}{{d \xi }} = {\rho}(0, \xi ), $$$ в массовых координатах последняя система запишется как
Эта система дополняется уравнениями, связывающими лагранжевы и эйлеровы координаты
$$$ \frac{dx}{dt} = u(t, \eta ), \quad \frac{{dx(0, \eta )}} {{d \eta }} = {v}(0, \eta ), $$$где $$v = \rho ^{ - 1}$$ — удельный объем.
Если система уравнений газодинамики записана в
то запись разностных схем, соответствующих методам Лакса - Вендроффа и Мак - Кормака, аналогична их записи для численного решения уравнения переноса (лекция 3).
Так, схема Лакса - Вендроффа может быть представлена в следующем виде:
$$\begin{gather*} \frac{\mathbf{\tilde{R}}_{m + 1/2} - 0.5(\mathbf{R}_{m + 1}^n + \mathbf{R}_m^n)}{\tau/2} + \frac{\mathbf{Q}_{m + 1}^n - \mathbf{Q}_m^n}{h} = 0, \\ \frac{\mathbf{\tilde{R}}_{m - 1/2} - 0.5(\mathbf{R}_m^n + \mathbf{R}_{m - 1}^n)} {\tau /2} + \frac{\mathbf{Q}_m^n - \mathbf{Q}_{m - 1}^n}{h} = 0. \end{gather*}$$(первый этап);
$$$ \frac{{{\mathbf{R}}_m^{n + 1} - {\mathbf{R}}_m^{n}}}{\tau} + \frac{{{\mathbf{\tilde{Q}}}_{m + 1/2} - {\mathbf{\tilde{Q}}}_{m - 1/2}}}{h} = 0 $$$(второй этап).
Схему МакКормака представим следующим образом:
$$\begin{gather*} \frac{\mathbf{\tilde{R}}_m - \mathbf{R}_m^n}{\tau} + \frac{\mathbf{Q}_{m + 1} - \mathbf{Q}_m}{h} = 0, \\ \frac{\mathbf{\tilde{R}}_{m - 1} - \mathbf{R}_{m - 1}^n}{\tau} + \frac{\mathbf{Q}_m - \mathbf{Q}_{m - 1}}{h} = 0 \end{gather*}$$(первый этап);
$$$ \frac{{{\mathbf{R}}_m^{n + 1} - 0.5({\mathbf{R}}_m^{n} + {\mathbf{\tilde{R}}}_m )}}{\tau} + \frac{{{\mathbf{\tilde{Q}}}_m - {\mathbf{\tilde{Q}}}_{m - 1}}}{{2h}} = 0 $$$(второй этап).
Так как использована
Этот метод (точнее, семейство методов) наиболее эффективен при численном решении задач, для которых существенным являются учет волновых процессов в сплошной среде. Подробнее о методе в книге [14.6], детальное описание — в монографии [14.7].
Матрица $$\mathbf{A}$$ системы
имеет вещественные собственные числа $$\lambda _{1} = u + c$$, $$\lambda _{2} = u$$, $$\lambda _{3} = u - c$$ и соответствующие им собственные векторы, $$\Omega _{i}$$ (например, левые, которые находятся при решении
Матрица из левых
Матрица $$\mathbf{A}$$ представима в виде $$\mathbf{A} = \mathbf{\Omega}^{- 1}\mathbf{{\Lambda}\Omega}$$, где $$\Lambda$$ — диагональная матрица, состоящая из собственных чисел матрицы $$\mathbf{A}$$:
$$\Lambda = diag(\lambda _{1}, \lambda _{2}, \lambda _{3}).$$Для получения разностной схемы введем обозначения
$$\begin{gather*} {\Lambda}^{+} = \frac{1}{2}({{\Lambda}} + \left|{{\Lambda}}\right|), {\Lambda}^{-} = \frac{1}{2}({{\Lambda}} - \left| {{\Lambda}}\right|), \\ {\mathbf{A}}^{+}= \frac{1}{2}({\mathbf{A}} + \left| {\mathbf{A}}\right|), {\mathbf{A}}^{-} = \frac{1}{2}({\mathbf{A}} - \left| {\mathbf{A}}\right|), {\sigma}= \tau/h. \end{gather*}$$
(рис 4.1) Умножим каждое из уравнений исходной системы газодинамики, записанных в
или, учитывая, что $$\omega _{i}$$ — левый
Построим разностную аппроксимацию полученной системы уравнений в частных производных с учетом знака собственных чисел (или направления характеристик). Для облегчения восприятия ограничимся простейшим методом первого порядка аппроксимации
$$$ {\omega }_i \frac{{{\mathbf{U}}_m^{n + 1} - {\mathbf{U}}_m^{n}}}{\tau} \mp \lambda_i{\omega }_i \frac{{{\mathbf{U}}_{m \mp 1}^{n} - {\mathbf{U}}_m^{n}}}{h} + {\omega }_i {\mathbf{f}} = 0. $$$Приведенное выше разностное уравнение может быть записано в матричной форме
$$\begin{gather*} {\Omega }({\mathbf{U}}_m^{n + 1} - {\mathbf{U}}_m^{n} ) - {\sigma}{{\Lambda}}^{+} {\Omega }({\mathbf{U}}_{m - 1}^{n} - {\mathbf{U}}_m^{n} ) + \\ + {\sigma}{{\Lambda}}^{-}{\Omega }({\mathbf{U}}_{m + 1}^{{n}} - {\mathbf{U}}_m^{n} ) + {\tau}{\Omega \mathbf{f}}_m^{n} = 0, \end{gather*}$$или в виде, разрешенном относительно верхнего слоя:
$$\begin{gather*} {\mathbf{U}}_m^{n + 1} = {\mathbf{U}}_m^{n} - \sigma \left[({\Omega }^{- 1}{{\Lambda}}^{+} {\Omega })_m^{n}({\mathbf{U}}_{m - 1}^{n} - {\mathbf{U}}_m^{n} ) -\right. \\ \left. - ({\Omega }^{- 1}{{\Lambda}}^{-} {\Omega })_m^{n} ({\mathbf{U}}_{m + 1}^{n} - {\mathbf{U}}_m^{n} )\right] + {\tau}{\mathbf{f}}_m^{n}, \end{gather*} $$или в компактной форме:
$${\mathbf{U}}_m^{n + 1} = {\mathbf{U}}_m^{n} + {\sigma}(- {\mathbf{B}}^{+} {\Delta}^{-}{\mathbf{U}} + {\mathbf{B}}^{-}{\Delta}^{+}{\mathbf{U}})_m^{n} + {\tau}{\mathbf{f}}_m^{n},$$где
$${\mathbf{B}}^{+} = {{\Omega}}^{- 1}{{\Lambda}}^{+}{{\Omega }}, {\mathbf{B}}^{-} = {\Omega }^{- 1}{{\Lambda}}^{-} {{\Omega}}.$$Учет направления характеристик позволяет получать устойчивые разностные
схемы для системы
Система уравнений газодинамики решается в области $$t \ge 0$$, $$\eta \in \left[{0, X}\right]$$, отрезок интегрирования разбивается на интервалы узлами $$\left\{{\eta_i }\right\}_0^{N}$$. Все интервалы заполнены газом, что соответствует приближению механики сплошной среды. Величины на интервалах считаются кусочно - постоянными. При численном решении определяются шесть функций: $$u, p, \varepsilon , \{ T, v, x\}$$ — скорость, давление, удельная внутренняя энергия, температура, удельный объем, эйлерова координата.
Для решения задачи используется система одномерных нестационарных уравнений в частных производных, описывающих поведение газа в лагранжевых переменных ( $$\eta$$ — лагранжева координата):
$$\begin{gather*} \frac{{dv}}{dt} - \frac{du}{{{\partial}\eta }} = 0, \\ \frac{du}{dt} + \frac{{{\partial}(p + Q)}}{{{\partial}\eta }} = 0, \\ \frac{d}{dt} \left({\varepsilon + \frac{{u^2}}{2}}\right) + \frac{{{\partial}\left[{(p + Q)u}\right]}}{{{\partial}\eta }} = \frac{{\partial}}{{{\partial}\eta }} \left[{a({T, v}) \frac{{\partial}t}{{{\partial}\eta }}}\right], \\ \frac{dx}{dt} = u. \end{gather*}$$Здесь a(T, v) — заданный коэффициент теплопроводности,
- Q = 0,
Коэффициент теплопроводности зависит от температуры и задается как
$$a(T, v) = T^{\gamma} a(T, v), \quad \gamma > 1, \quad a_1 \le a \le a_2 , \quad a_1 > 0.$$Начальные данные для рассматриваемой задачи будут
$$u(0, \eta ) = u_{0},\ x(0, \eta ) = x_{0},\ v(0, \eta ) = v_{0},\ T(0, \eta ) = T_{0} .$$Условия на границах выбираем следующие:
$$$ u(t, 0) = u_1 (t), \quad a \frac{{\partial}T}{{\partial}\eta }(t, 0) = P_1 (t), \quad T(t, X) = T_2 (t), $$$При u, T,
, v, x }известны на n слое и задача состоит в вычислении этих же функций на n + 1 слое. Сеточные функции { $$u_{m^{n + 1}}, T_{m + 1/2}^{n + 1}, v_{m + 1/2}^{n + 1}, x_{m^{n + 1}}$$ } вычисляются с помощью
В результате разностная схема записывается как
$$\begin{gather*} \frac{{{{v}}_{m + 1/2}^{n + 1} - {{v}}_{m + 1/2}^{n}}}{\tau} - \frac{{u_{m + 1}^{n + 1} - u_m^{n + 1}}}{2} = 0, m = 1, \ldots , M, \\ \frac{u_m^{n + 1} - u_m^{n}}{\tau} + \frac{{(p + Q)}_{m + 1/2}^{n + 1} - {(p + Q)}_{m - 1/2}^{n + 1}}{h_m} = 0 \\ {h_m = \eta_{m + 1/2} - \eta_{m - 1/2}, m = 1, \ldots , M - 1, } \\ \frac{1}{\tau} \left[{\varepsilon_{m + 1/2}^{n + 1} + \frac{{(u_m^{n + 1})^2 + (u_{m + 1}^{n + 1} )^2}}{4} - \varepsilon_{m + 1/2}^{n} - \frac{{(u_m^{n} )^2 + (u_{m + 1}^{n} )^2}}{4}}\right] + \\ + \frac{1}{{h_{m + 1/2}}} \left[{(p + Q)_{m + 1}^{n + 1}u_{m + 1}^{n + 1} - (p + Q)_m^{n + 1}u_m^{n + 1}}\right] = \\ = \frac{1}{{h_{m + 1/2}}} \left[{a_{m + 1} \frac{{T_{m + 3/2}^{n + 1} - T_{m + 1/2}^{n + 1}}}{{h_{m + 1}}} - a_m \frac{{T_{m + 1/2}^{n + 1} - T_{m - 1/2}^{n + 1}}}{{h_m }}}\right], \\ h_{m + 1/2} = \eta_{m + 1} - \eta_m , m = 1, \ldots , M - 1, \\ \frac{x_m^{n + 1} - x_m^{n}}{\tau} - \frac{u_m^{n + 1} - u_m^{n}}{2} = 0, \quad m = 0, \ldots , M. \end{gather*}$$К этим уравнениям системы добавим разностное выражение для вычисления
с фиктивным краевым условием QM + 1/2 = 0 и интерполяционное выражение для pm
Схема имеет второй порядок аппроксимации по координате, а при весовом коэффициенте, равном 0, 5, и второй порядок по времени, причем все точки спектра лежат на единичной окружности. Подробнее о данной схеме в [14.10].
Метод
Область интегрирования покрывается фиксированной в пространстве расчетной
сеткой, шаг которой h постоянен по обеим координатам x, y, ячейки занумерованы двумя индексами k, l.
В центре ячейки вычисляются величины $$u_{1kl}^{n}, \quad
u_{2kl}^{n}$$ (компоненты скорости газа), $$\varepsilon_{ikl}^{n}, m_{ikl}^{n}$$ где i — номер вещества. $$\varepsilon_{ikl}^{n}$$ — удельная внутренняя энергия газа с номером i, $$m_{ikl}^{n}$$ — масса этого вещества. Если
этого вещества в ячейке нет, то в ней и энергия, и масса полагаются равными нулю.
Предположим, что в каждой ячейке содержится несколько частиц (5 - 10), каждая из которых характеризуется координатами $$X_j^{n}, Y_j^{n}$$ массой $$\mu _{j}, i_{j}$$ — номер вещества, из которого состоит частица с номером j.
Шаг tn + 1 по вычисленным величинам $$\left\{u_1, u_2, \varepsilon_i, m_i\right\}_{kl}^{n}, \left\{{X, Y}\right\}_j^{n}$$ на нижнем слое tn.
На первом этапе расчета учитываются изменения искомых функций только за счет сил давления. При этом предположении разностные соотношения аппроксимируют уравнения
$$\begin{gather*} \frac{{{\partial}{\rho}}}{{\partial}t} = 0, \\ \frac{{{\partial}({\rho}u_1 )}}{{\partial}t} = - \frac{{\partial}p}{{\partial}x}, \\ \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} = - \frac{{\partial}p}{{\partial}y}, \\ \frac{{{\partial}E}}{{\partial}t} + \frac{{{\partial}(pu_1 )}}{{\partial}x} + \frac{{{\partial}(pu_2 )}}{{\partial}y} = 0, \end{gather*}$$где
$$$ E = {\rho}e = {\rho}\left[{\varepsilon + \frac{1}{2}(u_1^2 + u_2^2 )}\right]. $$$В расчетах участвуют также уравнения состояния для каждого газа
$$p_{i} = F_{i}(\varepsilon _{i}, \rho _{i}).$$На втором этапе аппроксимируются конвективные члены
$$\begin{gather*} \frac{{\partial}{\rho}}{{\partial}t} + \frac{{\partial}({\rho}u_1)}{{\partial}x} + \frac{{\partial}({\rho}u_2 )}{{\partial}y} = 0, \\ \frac{{\partial}{\rho}}{{\partial}t} + \frac{{\partial}({\rho}u_1)}{{\partial}x} + \frac{{\partial}({\rho}u_2 )}{{\partial}y} = 0, \\ \frac{{\partial}({\rho}u_1 )}{{\partial}t} + \frac{{\partial}({\rho}u_1^2 )}{{\partial}x} + \frac{{\partial}({\rho}u_1 u_2)}{{\partial}y} = 0, \\ \frac{{\partial}({\rho}u_2)}{{\partial}t} + \frac{{\partial}({\rho}u_1u_2)}{{\partial}x} + \frac{{\partial}({\rho}u_2)}{{\partial}y} = 0, \\ \frac{{\partial}E}{{\partial}t} + \frac{{\partial}(Eu_1)}{{\partial}x} + \frac{{\partial}(Eu_2)}{{\partial}y} = 0. \end{gather*}$$Опишем вычислительную процедуру на первом этапе. Известны $$u_1^{n}$$, $$u_2^{n}$$, mn, $$\varepsilon ^{n}$$, Xn, Yn (остальные индексы для простоты изложения опускаются). Сначала рассчитывается давление $$p_{ij}^{n}$$, исходя из предложения равенства давлений на границе двух сред p1 = p2 = ..., или
К этим уравнениям добавляется условие $$\sum {\gamma_i = h^2, }$$ поскольку $$\gamma _{i}$$ — часть объема h2 ячейки, занимаемого газом с номером i. По известной массе i газа находим его плотность: $$\rho_i = m_i^{n}/{\gamma}_i$$, а по известной удельной энергии $$\varepsilon_i^{n}$$ - давление $$p_i = F_i (\varepsilon_i^{n}, \gamma_i^{- 1} m_i^{n})$$.
Затем по закону Дальтона находится давление, $$p_{ikl}^{n}$$, которое приписывается к центру ячейки k, l. Система нелинейных алгебраических уравнений решается, вообще говоря, итерационным методом. В случае $$p = \rho f_{i}(\varepsilon )$$ выписывается ее явное решение. Далее находим предварительные значения расчитываемых величин, которые обозначим как $$\bar{u}_1 , \bar{u}_2 , \bar{m}, \bar{\varepsilon}$$. Первое из уравнений $$\frac{{{\partial}\rho }}{{\partial}t} = 0$$ закон сохранения массы — на сетке приобретает вид $$\bar{m}_{ikl} = m_{ikl}^{n}$$. Поскольку на первом этапе $$\frac{{\partial}{\rho}}{{\partial}t} = 0$$, то два следующих разностных уравнения — уравнения движения — записываются как
Здесь
$$\begin{gather*} p_{k + 1/2, l}^{n} = \frac{1}{2}(p_{kl}^{n} + p_{k + 1, l}^{n} ), \\ m_{kl}^{n} = \sum\limits_i {m_{ikl}^{n}}, \rho_{kl}^{n} = m_{kl}^{n}/h^2 . \end{gather*}$$Последняя из рассчитываемых величин — энергия. Дискретный аналог уравнения энергии в
Здесь
$$\begin{gather*} u_{1, k + 1/2, l}^{{n} + 1/2} = \frac{{\bar{u}_{1, {kl}} + u_{1, {kl}}^{n} + \bar{u}_{1, k + 1, l} + u_{1 , k + 1, l}^{n}}}{4}, \\ u_{2 , {k, l} + 1/2}^{{n} + 1/2} = \frac{{\bar{u}_{2 {kl}} + u_{2 {kl}}^{n} + \bar{u}_{2, {k, l} + 1} + u_{2 {k, l} - 1}}}{4}. \end{gather*}$$Вычислим величину $$e_{kl}^{n}$$ — энергию. Напомним,
что $$$ e = {\rho}(\varepsilon + \frac{{u_1^2 + u_{2}^{2}}}{2}). $$$ Тогда eh2 есть энергия в ячейке h x h:
где $$\rho h^{2}$$ — масса ячейки, равная $$m_{kl}^{n} = \sum\limits_i {m_{ikl}^{n} }$$. В таком случае, с учетом закона сохранения массы,
$$$ \left[{(h^2 {\rho}) \frac{{u_1^2 + u_2^2}}{2}}\right]_{kl}^{n} = \frac{1}{2}m_{kl}^{n} \left[{(u_{1, {kl}}^{n} )^2 + (u_{2 , {kl}}^{n} )^2}\right]. $$$Вычислим величину $$\left[{(h^2 {\rho}) \varepsilon }\right]_{kl}^{n} $$, имеющую смысл полной внутренней энергии в ячейке, зная массу $$m_{ikl}^{n}$$ и удельную внутреннюю энергию $$\varepsilon_{ikl}^{n}$$ вещества.
$$\left[{(h^2 {\rho}) \varepsilon }\right]_{kl}^{n} = \sum\limits_i {m_{ikl}^{n}} E_{ikl}^{n},$$теперь имеется алгоритм вычисления $$E_{kl}^{n}$$ и, следовательно, $$\tilde {E}_{kl}^{n}$$.
Из соотношения
$$$ \bar{e}_{kl} \cdot h^2 = (h^2 \rho_{kl} ) \bar{\varepsilon}_{kl} + (h^2 \rho_{kl} ) \cdot \frac{{(\bar{u}_{1, {kl}} )^2 + (\bar{u}_{2, {kl}} )^2}}{2} $$$находим величину полной удельной внутренней энергии $$\bar{\varepsilon}_{kl}$$. Однако искомыми являются значения удельной внутренней энергии для каждого вещества. Пусть $$\Delta \varepsilon _{i}$$ — изменение удельной внутренней энергии i вещества за первый этап шага по времени по i веществу. Зная mi (масса i вещества), запишем полное приращение полной удельной внутренней энергии в ячейке ( $$\Delta _{\varepsilon }$$ ) и приравняем его к уже полученному
Для определения изменения количества каждого газа нужно сделать некое правдоподобное предположение, например, считать, что все $$\Delta \varepsilon _{ikl}$$ одинаковы. Тогда $${N} \Delta \varepsilon_{ikl} = \bar{\varepsilon}_{kl}^{n} - \varepsilon_{kl}^{n}$$ и, соответственно, $$\bar{\varepsilon}_{ikl} = \varepsilon_{ikl}^{n} + \Delta \varepsilon_{ikl}$$. На этом первый этап расчета (предиктор) закончен.
Рассмотрим второй этап расчета.
Движение частиц описывается
которые могут быть приближены, например, с использованием явного
где скорости частиц $$$ \tilde{u}_k, \tilde{v}_k $$$ определяются интерполяцией величин $$$ \bar{u}_{1}, \bar{u}_{2} $$$ в ячейках, окружающих p частицу.
После этого рассчитывается перенос массы и вычисляется новая масса каждой ячейки. Для этого выделяются три группы частиц:
n + 1 слой в пределах ячейки, которые, очевидно, не вносят изменений в массу, импульс, энергию новой ячейки, т.е.$$(X_p^{n}, Y_p^{n}) \in \omega_{ij}, (X_p^{n + 1}, Y_p^{n + 1}) \in \omega_{ij},$$
где $$\omega _{ij}$$ — обозначение старой ячейки,На шаг по времени накладывается ограничение
$$$ {\tau} < \frac{h}{{\sqrt {u_1^2 + u_2^2}}}, $$$что означает запрет на перемещение частицы за один шаг больше, чем на одну
ячейку. Перемещение в данном методе возможно только в соседнюю ячейку. Предположим, что каждая p частица, перешедшая на n + 1 шаге по времени в соседнюю ячейку, переносит в нее массу mp. Это означает, что значение массы mikl считается путем сложения масс всех частиц типа i, для которых $$(X_p^{n + 1}, Y_p^{n + 1}) \in \omega_{ij}$$
Процедура вычисления импульса выполняется следующим образом. Компоненты
полного импульса частиц в ячейке ( k, l ) могут быть вычислены, как $$m_{kl}^{n} \overline{u}_{1, kl}, m_{kl}^{n} \overline{u}_{2, kl}$$, при этом p частица, покинувшая ячейку $$\omega _{ij}$$, уносит импульс $$m_p \overline{u}_{1, {kl}}, m_p u_{2, kl}$$. Изменение импульса в Skl за один шаг по времени будет
здесь символы суммирования означают суммирование по частицам ( p ),
покинувшим данную ячейку и пришедшим в нее, соответственно. После вычисления компонентов импульса каждой ячейки вычисляются компоненты скорости $$u_{1kl}^{n + 1}, u_{2kl}^{n + 1}$$.
Частица типа i, переходящая из Ske в другую
ячейку, переносит полную энергию
Тогда можно вычислить энергию i вещества в ячейке $$\omega _{kl}$$ на промежуточном шаге:
При t = tn + 1 полная удельная энергия изменится на величину
где знаки суммирования снова означают суммы по всем частицам, покинувшим ячейку $$\omega _{kl}$$ и пришедшим в нее, соответственно. Далее получим
$$$ \varepsilon_{ikl}^{n + 1} = \frac{{h^2}}{{m_{ikl}^{n + 1}}}E_{ikl}^{n + 1} - \frac{1}{2} \left[{(u_{1 {kl}}^{n + 1} )^2 + (u_{2 {kl}}^{n + 1} )^2}\right]. $$$Отметим недостатки этого метода. Во - первых, это дискретность плотности, что приводит при небольшом количестве частиц в ячейке к скачкообразным изменениям плотности. Во - вторых, аппроксимация исходных уравнений достигается при количестве частиц, стремящемся к бесконечности. Кроме того, этот метод требует значительно большего количества памяти, чем конечно - разностные методы, так как наряду с физическими характеристиками узлов необходимо хранить и свойства частиц в ячейках.
Подробное описание
Выпишем нелинейную систему уравнений одномерных движений идеальной сжимаемой жидкости в случае баротропных процессов. Она состоит из уравнения Эйлера
$$$ \frac{{{\partial} u}}{{{\partial} t}} + u \frac{{{\partial} u}}{{{\partial} x}} + \frac{1}{\rho} \frac{{{\partial} p}}{{{\partial} x}} = 0, $$$уравнения неразрывности
$$$ \frac{{{\partial}{\rho}}}{{{\partial} t}} + \rho \frac{{\partial}u} {{\partial}x} + u \frac{{{\partial}{\rho}}}{{\partial}x} = 0 $$$и условия баротропности
$$p = f(\rho ).$$Уравнения позволяют определить плотность $$\rho$$ и скорость u в зависимости от координаты x и времени t. Система не имеет решений, зависящих только от $$x \pm a_{0}t$$, но оказывается возможным найти решение этой системы, представляющее собой плоскую волну и являющееся обобщением решений вида $$f(x \pm a_{0}t)$$. Будем искать такие решения системы, для которых скорость u является функцией только плотности $$\rho.$$
$$$ u = \pm \int {\sqrt {\frac{{{dp}}}{{d {\rho}}}} \frac{{d \rho }}{\rho}}. $$$
Обозначим $$$ \frac{dp}{d{\rho}} = a^2 ({\rho}) $$$ и введем величину c = u + a. Какой физический смысл имеет величина c?
$$$ \frac{{\partial}\rho}{{\partial}T} + c(\rho ) \frac{{\partial}\rho}{{\partial}x} = 0, $$$
$$$ c \left({\rho}\right) = \sqrt {A \gamma } \left[{1 + \frac{2}{{\gamma - 1}}}\right] {\rho}^{{{1 \over 2}}(\gamma - 1)}, $$$
$$u = u_0 f({x/t}), {\rho} = \rho_0{\varphi}({x/t})$$.
Течения подобного типа — частный случай автомодельных течений, когда решение зависит от некоторой комбинации независимых переменных. Проиллюстрировать численными расчетами особенности распространения центрированных волн Римана.
Рассмотренные в предыдущих лекциях разностные методы решения уравнений в частных производных демонстрировались на примере либо линейных задач, либо достаточно простых нелинейных уравнений с хорошо изученными свойствами. Такие задачи в реальной вычислительной практике обычно служат тестами для отработки методов решения более сложных нелинейных систем. Традиционным объектом приложения численных методов служат уравнения механики сплошной среды (МСС). Основных причин тому три. Первая — все математические модели МСС известны достаточно давно, хорошо исследованы и являются частью научной классики. Вторая — математические модели МСС нелинейны. Как правило, у линеаризованных уравнений очень узкая область применения. Третья — в результатах решения задач МСС на протяжении всего XX века была практическая заинтересованность, вызванная бурным развитием авиации, осуществлением наукоемких ядерных и космических программ в разных странах.
В данной книге ограничимся самой простой моделью МСС —
В основе построения математических моделей, описывающих поведение жидкостей и газов, лежит понятие о сплошной среде. Из молекулярной физики известно, что среда состоит из отдельных частиц (молекул, ионов, электронов, атомов), расстояние между которым существенно больше их собственных размеров. Длина свободного пробега частицы l (расстояние, пройденное частицей между двумя столкновениями) тем меньше, чем больше частиц заключено в единице объема, чем больше плотность среды. В механике жидкостей и газов рассматриваются среды, содержащие в единице объема большое количество частиц (много больше, чем число Авогадро — число частиц в одной грамм - молекуле вещества, $$N_A \approx 6 \cdot 10^{23} \quad (\mbox{г} \cdot \mbox{моль})^{- 1}$$ ).
В таких средах можно рассматривать лишь некоторые усредненные характеристики, не занимаясь изучением поведения каждой частицы в отдельности. В этом предположении заключается идея модели сплошной среды, непрерывно заполняющей пространство.
Количественным критерием применимости приближения сплошной среды может
служить неравенство $$l/L \ll 1$$, где L — характерный пространственный размер задачи (например, размер тела при внешнем обтекании потоком газа). В газах при нормальных условиях $$l \approx (10^{- 5} \div 10^{- 6})$$
см, поэтому приведенное условие для тел с размером более 1 см выполняется с достаточной точностью. Предложение о сплошности среды, по - видимому, берет свое начало от Эйлера, впервые рассмотревшего газ как непрерывную деформируемую субстанцию. Не вдаваясь в особенности получения
Эйлерова (недивергентная) форма одномерной системы уравнений газодинамики имеет вид
$$\begin{gather*} \frac{{\partial}{\rho}}{{\partial}t} + u \frac{{\partial}{\rho}}{{\partial x}} + {\rho}\frac{{\partial}u}{{\partial}x} = 0 \quad \mbox{(уравнение неразрывности), }\\ \frac{{\partial}u}{{\partial}t} + u \frac{{\partial u}}{{\partial}x} = - \frac{1}{\rho} \cdot \frac{{\partial}p}{{\partial}x} \quad \mbox{(уравнение движения), } \\ \frac{{\partial}e}{{\partial}t} + u \frac{\partial e}{{\partial}x} = - \frac{1}{\rho} \frac{{\partial}(pu)}{{\partial}x} \quad \mbox{(уравнение энергии), } \end{gather*}$$где e — удельная энергия, равная $$$ e = \varepsilon + \frac{u^2}{2}, \varepsilon $$$ удельная внутренняя энергия, u — скорость газа, $$\rho$$ — плотность среды, p — давление, t — время, x — декартова координата. Эту же систему уравнений в частных производных можно представить в матричной (характеристической) форме
где $$\mathbf{U}$$ — вектор - столбец, $$\mathbf{A}$$ — квадратная матрица 3 x 3.
Система замыкается уравнением состояния
$$f(\rho , \varepsilon , p) = 0,$$которое, например, для идеального газа имеет вид
$$\frac{p}{{\rho}(\gamma - 1)} - \varepsilon = 0.$$где $$\gamma$$ — безразмерная постоянная, равная отношению теплоемкости газа при постоянном давлении и теплоемкости при постоянном объеме — постоянная адиабаты.
В системе, записанной в
Не занимаясь выводом формул (это делается простыми алгебраическими преобразованиями), представим другие виды записи уравнения энергии, справедливые для приведенного выше уравнения состояния:
$$\begin{gather*} \frac{{\partial}p}{{\partial}t} + u \frac{{\partial p}}{{\partial}x} + \gamma p \frac{{\partial}u}{{\partial}x} = 0, \\ \frac{{\partial}p}{{\partial}t} + u \frac{{\partial p}}{{\partial}x} + c^2 {\rho}\frac{{\partial}u}{{\partial}x} = 0, \end{gather*}$$где $$c = \sqrt{\gamma {p/{\rho}}}$$ — адиабатическая скорость звука,
$$\begin{gather*} \frac{{{\partial}S}}{{\partial}t} + u \frac{{{\partial}S}}{{\partial}x} = 0, \\ \mbox{или} \quad \frac{dS}{dt} \left|_{tr}\right. = 0, \quad \mbox{где} \quad \frac{d}{dt} \left|_{tr}\right. = \frac{\partial }{{\partial}t} + u \frac{{\partial}}{{\partial}x}, \end{gather*}$$$$S = m(\varepsilon \rho )^{1 - \gamma }$$ — энтропия. Она, как следует из последней формулы, сохраняется вдоль траектории частицы идеального газа, т.е. на траектории уравнения
$$$ \frac{dX}{dt} = u(t, X), X(0) = X_0. $$$или в
Интегральная форма этих уравнений получается при использовании теоремы Гаусса - Остроградского
$$$ \int\int\limits_{\Omega}(\frac{{{\partial}{\mathbf{R}}}}{{\partial}t} + \frac{{{\partial}{\mathbf{Q}}}}{{\partial}x}){dtdx} = \oint\limits_{\Gamma}(\mathbf{R}dx - \mathbf{Q}dt) = 0, $$$где $$\Gamma$$ — граница замкнутой области интегрирования $$\Omega$$ в плоскости t, x. Здесь $$u(t, x), \rho (t, x), p(t, x)$$ — скорость, плотность и давление газа соответственно. В механике сплошных сред вводится эйлерово и лагранжево описание поведения среды. В первом случае наблюдатель полагается неподвижным, например, стоящим на берегу реки. Соответственно расчетная сетка будет неподвижной (фиксированная эйлерова сетка). Во втором случае полагаем, что наблюдатель движется вместе со средой, например, находится на лодке, плывущей по течению реки. В этом случае лагранжева расчетная сетка будет двигаться вместе с частицами среды.
Лагранжева форма одномерных
где $$$ I = \left({\frac{{\partial}x}{{\partial}\xi }}\right)^{-
1}, $$$ x и $$\xi$$ соответственно эйлерова и лагранжева
координаты, связь между которыми дается последним уравнением. Эта система в одномерном случае может быть записана в другом виде. Если ввести лагранжеву массовую координату $$\eta (\xi )$$, связанную с лагранжевой координатой $$\xi$$ дифференциальным уравнением $$$ \frac{{d \eta }}{{d \xi }} = {\rho}(0, \xi ), $$$ в массовых координатах последняя система запишется как
Эта система дополняется уравнениями, связывающими лагранжевы и эйлеровы координаты
$$$ \frac{dx}{dt} = u(t, \eta ), \quad \frac{{dx(0, \eta )}} {{d \eta }} = {v}(0, \eta ), $$$где $$v = \rho ^{ - 1}$$ — удельный объем.
Если система уравнений газодинамики записана в
то запись разностных схем, соответствующих методам Лакса - Вендроффа и Мак - Кормака, аналогична их записи для численного решения уравнения переноса (лекция 3).
Так, схема Лакса - Вендроффа может быть представлена в следующем виде:
$$\begin{gather*} \frac{\mathbf{\tilde{R}}_{m + 1/2} - 0.5(\mathbf{R}_{m + 1}^n + \mathbf{R}_m^n)}{\tau/2} + \frac{\mathbf{Q}_{m + 1}^n - \mathbf{Q}_m^n}{h} = 0, \\ \frac{\mathbf{\tilde{R}}_{m - 1/2} - 0.5(\mathbf{R}_m^n + \mathbf{R}_{m - 1}^n)} {\tau /2} + \frac{\mathbf{Q}_m^n - \mathbf{Q}_{m - 1}^n}{h} = 0. \end{gather*}$$(первый этап);
$$$ \frac{{{\mathbf{R}}_m^{n + 1} - {\mathbf{R}}_m^{n}}}{\tau} + \frac{{{\mathbf{\tilde{Q}}}_{m + 1/2} - {\mathbf{\tilde{Q}}}_{m - 1/2}}}{h} = 0 $$$(второй этап).
Схему МакКормака представим следующим образом:
$$\begin{gather*} \frac{\mathbf{\tilde{R}}_m - \mathbf{R}_m^n}{\tau} + \frac{\mathbf{Q}_{m + 1} - \mathbf{Q}_m}{h} = 0, \\ \frac{\mathbf{\tilde{R}}_{m - 1} - \mathbf{R}_{m - 1}^n}{\tau} + \frac{\mathbf{Q}_m - \mathbf{Q}_{m - 1}}{h} = 0 \end{gather*}$$(первый этап);
$$$ \frac{{{\mathbf{R}}_m^{n + 1} - 0.5({\mathbf{R}}_m^{n} + {\mathbf{\tilde{R}}}_m )}}{\tau} + \frac{{{\mathbf{\tilde{Q}}}_m - {\mathbf{\tilde{Q}}}_{m - 1}}}{{2h}} = 0 $$$(второй этап).
Так как использована
Этот метод (точнее, семейство методов) наиболее эффективен при численном решении задач, для которых существенным являются учет волновых процессов в сплошной среде. Подробнее о методе в книге [14.6], детальное описание — в монографии [14.7].
Матрица $$\mathbf{A}$$ системы
имеет вещественные собственные числа $$\lambda _{1} = u + c$$, $$\lambda _{2} = u$$, $$\lambda _{3} = u - c$$ и соответствующие им собственные векторы, $$\Omega _{i}$$ (например, левые, которые находятся при решении
Матрица из левых
Матрица $$\mathbf{A}$$ представима в виде $$\mathbf{A} = \mathbf{\Omega}^{- 1}\mathbf{{\Lambda}\Omega}$$, где $$\Lambda$$ — диагональная матрица, состоящая из собственных чисел матрицы $$\mathbf{A}$$:
$$\Lambda = diag(\lambda _{1}, \lambda _{2}, \lambda _{3}).$$Для получения разностной схемы введем обозначения
$$\begin{gather*} {\Lambda}^{+} = \frac{1}{2}({{\Lambda}} + \left|{{\Lambda}}\right|), {\Lambda}^{-} = \frac{1}{2}({{\Lambda}} - \left| {{\Lambda}}\right|), \\ {\mathbf{A}}^{+}= \frac{1}{2}({\mathbf{A}} + \left| {\mathbf{A}}\right|), {\mathbf{A}}^{-} = \frac{1}{2}({\mathbf{A}} - \left| {\mathbf{A}}\right|), {\sigma}= \tau/h. \end{gather*}$$(рис 4.1) Умножим каждое из уравнений исходной системы газодинамики, записанных в
или, учитывая, что $$\omega _{i}$$ — левый
Построим разностную аппроксимацию полученной системы уравнений в частных производных с учетом знака собственных чисел (или направления характеристик). Для облегчения восприятия ограничимся простейшим методом первого порядка аппроксимации
$$$ {\omega }_i \frac{{{\mathbf{U}}_m^{n + 1} - {\mathbf{U}}_m^{n}}}{\tau} \mp \lambda_i{\omega }_i \frac{{{\mathbf{U}}_{m \mp 1}^{n} - {\mathbf{U}}_m^{n}}}{h} + {\omega }_i {\mathbf{f}} = 0. $$$Приведенное выше разностное уравнение может быть записано в матричной форме
$$\begin{gather*} {\Omega }({\mathbf{U}}_m^{n + 1} - {\mathbf{U}}_m^{n} ) - {\sigma}{{\Lambda}}^{+} {\Omega }({\mathbf{U}}_{m - 1}^{n} - {\mathbf{U}}_m^{n} ) + \\ + {\sigma}{{\Lambda}}^{-}{\Omega }({\mathbf{U}}_{m + 1}^{{n}} - {\mathbf{U}}_m^{n} ) + {\tau}{\Omega \mathbf{f}}_m^{n} = 0, \end{gather*}$$или в виде, разрешенном относительно верхнего слоя:
$$\begin{gather*} {\mathbf{U}}_m^{n + 1} = {\mathbf{U}}_m^{n} - \sigma \left[({\Omega }^{- 1}{{\Lambda}}^{+} {\Omega })_m^{n}({\mathbf{U}}_{m - 1}^{n} - {\mathbf{U}}_m^{n} ) -\right. \\ \left. - ({\Omega }^{- 1}{{\Lambda}}^{-} {\Omega })_m^{n} ({\mathbf{U}}_{m + 1}^{n} - {\mathbf{U}}_m^{n} )\right] + {\tau}{\mathbf{f}}_m^{n}, \end{gather*} $$или в компактной форме:
$${\mathbf{U}}_m^{n + 1} = {\mathbf{U}}_m^{n} + {\sigma}(- {\mathbf{B}}^{+} {\Delta}^{-}{\mathbf{U}} + {\mathbf{B}}^{-}{\Delta}^{+}{\mathbf{U}})_m^{n} + {\tau}{\mathbf{f}}_m^{n},$$где
$${\mathbf{B}}^{+} = {{\Omega}}^{- 1}{{\Lambda}}^{+}{{\Omega }}, {\mathbf{B}}^{-} = {\Omega }^{- 1}{{\Lambda}}^{-} {{\Omega}}.$$Учет направления характеристик позволяет получать устойчивые разностные
схемы для системы
Система уравнений газодинамики решается в области $$t \ge 0$$, $$\eta \in \left[{0, X}\right]$$, отрезок интегрирования разбивается на интервалы узлами $$\left\{{\eta_i }\right\}_0^{N}$$. Все интервалы заполнены газом, что соответствует приближению механики сплошной среды. Величины на интервалах считаются кусочно - постоянными. При численном решении определяются шесть функций: $$u, p, \varepsilon , \{ T, v, x\}$$ — скорость, давление, удельная внутренняя энергия, температура, удельный объем, эйлерова координата.
Для решения задачи используется система одномерных нестационарных уравнений в частных производных, описывающих поведение газа в лагранжевых переменных ( $$\eta$$ — лагранжева координата):
$$\begin{gather*} \frac{{dv}}{dt} - \frac{du}{{{\partial}\eta }} = 0, \\ \frac{du}{dt} + \frac{{{\partial}(p + Q)}}{{{\partial}\eta }} = 0, \\ \frac{d}{dt} \left({\varepsilon + \frac{{u^2}}{2}}\right) + \frac{{{\partial}\left[{(p + Q)u}\right]}}{{{\partial}\eta }} = \frac{{\partial}}{{{\partial}\eta }} \left[{a({T, v}) \frac{{\partial}t}{{{\partial}\eta }}}\right], \\ \frac{dx}{dt} = u. \end{gather*}$$Здесь a(T, v) — заданный коэффициент теплопроводности,
- Q = 0,
Коэффициент теплопроводности зависит от температуры и задается как
$$a(T, v) = T^{\gamma} a(T, v), \quad \gamma > 1, \quad a_1 \le a \le a_2 , \quad a_1 > 0.$$Начальные данные для рассматриваемой задачи будут
$$u(0, \eta ) = u_{0},\ x(0, \eta ) = x_{0},\ v(0, \eta ) = v_{0},\ T(0, \eta ) = T_{0} .$$Условия на границах выбираем следующие:
$$$ u(t, 0) = u_1 (t), \quad a \frac{{\partial}T}{{\partial}\eta }(t, 0) = P_1 (t), \quad T(t, X) = T_2 (t), $$$При u, T,
, v, x }известны на n слое и задача состоит в вычислении этих же функций на n + 1 слое. Сеточные функции { $$u_{m^{n + 1}}, T_{m + 1/2}^{n + 1}, v_{m + 1/2}^{n + 1}, x_{m^{n + 1}}$$ } вычисляются с помощью
В результате разностная схема записывается как
$$\begin{gather*} \frac{{{{v}}_{m + 1/2}^{n + 1} - {{v}}_{m + 1/2}^{n}}}{\tau} - \frac{{u_{m + 1}^{n + 1} - u_m^{n + 1}}}{2} = 0, m = 1, \ldots , M, \\ \frac{u_m^{n + 1} - u_m^{n}}{\tau} + \frac{{(p + Q)}_{m + 1/2}^{n + 1} - {(p + Q)}_{m - 1/2}^{n + 1}}{h_m} = 0 \\ {h_m = \eta_{m + 1/2} - \eta_{m - 1/2}, m = 1, \ldots , M - 1, } \\ \frac{1}{\tau} \left[{\varepsilon_{m + 1/2}^{n + 1} + \frac{{(u_m^{n + 1})^2 + (u_{m + 1}^{n + 1} )^2}}{4} - \varepsilon_{m + 1/2}^{n} - \frac{{(u_m^{n} )^2 + (u_{m + 1}^{n} )^2}}{4}}\right] + \\ + \frac{1}{{h_{m + 1/2}}} \left[{(p + Q)_{m + 1}^{n + 1}u_{m + 1}^{n + 1} - (p + Q)_m^{n + 1}u_m^{n + 1}}\right] = \\ = \frac{1}{{h_{m + 1/2}}} \left[{a_{m + 1} \frac{{T_{m + 3/2}^{n + 1} - T_{m + 1/2}^{n + 1}}}{{h_{m + 1}}} - a_m \frac{{T_{m + 1/2}^{n + 1} - T_{m - 1/2}^{n + 1}}}{{h_m }}}\right], \\ h_{m + 1/2} = \eta_{m + 1} - \eta_m , m = 1, \ldots , M - 1, \\ \frac{x_m^{n + 1} - x_m^{n}}{\tau} - \frac{u_m^{n + 1} - u_m^{n}}{2} = 0, \quad m = 0, \ldots , M. \end{gather*}$$К этим уравнениям системы добавим разностное выражение для вычисления
с фиктивным краевым условием QM + 1/2 = 0 и интерполяционное выражение для pm
Схема имеет второй порядок аппроксимации по координате, а при весовом коэффициенте, равном 0, 5, и второй порядок по времени, причем все точки спектра лежат на единичной окружности. Подробнее о данной схеме в [14.10].
Метод
Область интегрирования покрывается фиксированной в пространстве расчетной
сеткой, шаг которой h постоянен по обеим координатам x, y, ячейки занумерованы двумя индексами k, l.
В центре ячейки вычисляются величины $$u_{1kl}^{n}, \quad
u_{2kl}^{n}$$ (компоненты скорости газа), $$\varepsilon_{ikl}^{n}, m_{ikl}^{n}$$ где i — номер вещества. $$\varepsilon_{ikl}^{n}$$ — удельная внутренняя энергия газа с номером i, $$m_{ikl}^{n}$$ — масса этого вещества. Если
этого вещества в ячейке нет, то в ней и энергия, и масса полагаются равными нулю.
Предположим, что в каждой ячейке содержится несколько частиц (5 - 10), каждая из которых характеризуется координатами $$X_j^{n}, Y_j^{n}$$ массой $$\mu _{j}, i_{j}$$ — номер вещества, из которого состоит частица с номером j.
Шаг tn + 1 по вычисленным величинам $$\left\{u_1, u_2, \varepsilon_i, m_i\right\}_{kl}^{n}, \left\{{X, Y}\right\}_j^{n}$$ на нижнем слое tn.
На первом этапе расчета учитываются изменения искомых функций только за счет сил давления. При этом предположении разностные соотношения аппроксимируют уравнения
$$\begin{gather*} \frac{{{\partial}{\rho}}}{{\partial}t} = 0, \\ \frac{{{\partial}({\rho}u_1 )}}{{\partial}t} = - \frac{{\partial}p}{{\partial}x}, \\ \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} = - \frac{{\partial}p}{{\partial}y}, \\ \frac{{{\partial}E}}{{\partial}t} + \frac{{{\partial}(pu_1 )}}{{\partial}x} + \frac{{{\partial}(pu_2 )}}{{\partial}y} = 0, \end{gather*}$$где
$$$ E = {\rho}e = {\rho}\left[{\varepsilon + \frac{1}{2}(u_1^2 + u_2^2 )}\right]. $$$В расчетах участвуют также уравнения состояния для каждого газа
$$p_{i} = F_{i}(\varepsilon _{i}, \rho _{i}).$$На втором этапе аппроксимируются конвективные члены
$$\begin{gather*} \frac{{\partial}{\rho}}{{\partial}t} + \frac{{\partial}({\rho}u_1)}{{\partial}x} + \frac{{\partial}({\rho}u_2 )}{{\partial}y} = 0, \\ \frac{{\partial}{\rho}}{{\partial}t} + \frac{{\partial}({\rho}u_1)}{{\partial}x} + \frac{{\partial}({\rho}u_2 )}{{\partial}y} = 0, \\ \frac{{\partial}({\rho}u_1 )}{{\partial}t} + \frac{{\partial}({\rho}u_1^2 )}{{\partial}x} + \frac{{\partial}({\rho}u_1 u_2)}{{\partial}y} = 0, \\ \frac{{\partial}({\rho}u_2)}{{\partial}t} + \frac{{\partial}({\rho}u_1u_2)}{{\partial}x} + \frac{{\partial}({\rho}u_2)}{{\partial}y} = 0, \\ \frac{{\partial}E}{{\partial}t} + \frac{{\partial}(Eu_1)}{{\partial}x} + \frac{{\partial}(Eu_2)}{{\partial}y} = 0. \end{gather*}$$Опишем вычислительную процедуру на первом этапе. Известны $$u_1^{n}$$, $$u_2^{n}$$, mn, $$\varepsilon ^{n}$$, Xn, Yn (остальные индексы для простоты изложения опускаются). Сначала рассчитывается давление $$p_{ij}^{n}$$, исходя из предложения равенства давлений на границе двух сред p1 = p2 = ..., или
К этим уравнениям добавляется условие $$\sum {\gamma_i = h^2, }$$ поскольку $$\gamma _{i}$$ — часть объема h2 ячейки, занимаемого газом с номером i. По известной массе i газа находим его плотность: $$\rho_i = m_i^{n}/{\gamma}_i$$, а по известной удельной энергии $$\varepsilon_i^{n}$$ - давление $$p_i = F_i (\varepsilon_i^{n}, \gamma_i^{- 1} m_i^{n})$$.
Затем по закону Дальтона находится давление, $$p_{ikl}^{n}$$, которое приписывается к центру ячейки k, l. Система нелинейных алгебраических уравнений решается, вообще говоря, итерационным методом. В случае $$p = \rho f_{i}(\varepsilon )$$ выписывается ее явное решение. Далее находим предварительные значения расчитываемых величин, которые обозначим как $$\bar{u}_1 , \bar{u}_2 , \bar{m}, \bar{\varepsilon}$$. Первое из уравнений $$\frac{{{\partial}\rho }}{{\partial}t} = 0$$ закон сохранения массы — на сетке приобретает вид $$\bar{m}_{ikl} = m_{ikl}^{n}$$. Поскольку на первом этапе $$\frac{{\partial}{\rho}}{{\partial}t} = 0$$, то два следующих разностных уравнения — уравнения движения — записываются как
Здесь
$$\begin{gather*} p_{k + 1/2, l}^{n} = \frac{1}{2}(p_{kl}^{n} + p_{k + 1, l}^{n} ), \\ m_{kl}^{n} = \sum\limits_i {m_{ikl}^{n}}, \rho_{kl}^{n} = m_{kl}^{n}/h^2 . \end{gather*}$$Последняя из рассчитываемых величин — энергия. Дискретный аналог уравнения энергии в
Здесь
$$\begin{gather*} u_{1, k + 1/2, l}^{{n} + 1/2} = \frac{{\bar{u}_{1, {kl}} + u_{1, {kl}}^{n} + \bar{u}_{1, k + 1, l} + u_{1 , k + 1, l}^{n}}}{4}, \\ u_{2 , {k, l} + 1/2}^{{n} + 1/2} = \frac{{\bar{u}_{2 {kl}} + u_{2 {kl}}^{n} + \bar{u}_{2, {k, l} + 1} + u_{2 {k, l} - 1}}}{4}. \end{gather*}$$Вычислим величину $$e_{kl}^{n}$$ — энергию. Напомним,
что $$$ e = {\rho}(\varepsilon + \frac{{u_1^2 + u_{2}^{2}}}{2}). $$$ Тогда eh2 есть энергия в ячейке h x h:
где $$\rho h^{2}$$ — масса ячейки, равная $$m_{kl}^{n} = \sum\limits_i {m_{ikl}^{n} }$$. В таком случае, с учетом закона сохранения массы,
$$$ \left[{(h^2 {\rho}) \frac{{u_1^2 + u_2^2}}{2}}\right]_{kl}^{n} = \frac{1}{2}m_{kl}^{n} \left[{(u_{1, {kl}}^{n} )^2 + (u_{2 , {kl}}^{n} )^2}\right]. $$$Вычислим величину $$\left[{(h^2 {\rho}) \varepsilon }\right]_{kl}^{n} $$, имеющую смысл полной внутренней энергии в ячейке, зная массу $$m_{ikl}^{n}$$ и удельную внутреннюю энергию $$\varepsilon_{ikl}^{n}$$ вещества.
$$\left[{(h^2 {\rho}) \varepsilon }\right]_{kl}^{n} = \sum\limits_i {m_{ikl}^{n}} E_{ikl}^{n},$$теперь имеется алгоритм вычисления $$E_{kl}^{n}$$ и, следовательно, $$\tilde {E}_{kl}^{n}$$.
Из соотношения
$$$ \bar{e}_{kl} \cdot h^2 = (h^2 \rho_{kl} ) \bar{\varepsilon}_{kl} + (h^2 \rho_{kl} ) \cdot \frac{{(\bar{u}_{1, {kl}} )^2 + (\bar{u}_{2, {kl}} )^2}}{2} $$$находим величину полной удельной внутренней энергии $$\bar{\varepsilon}_{kl}$$. Однако искомыми являются значения удельной внутренней энергии для каждого вещества. Пусть $$\Delta \varepsilon _{i}$$ — изменение удельной внутренней энергии i вещества за первый этап шага по времени по i веществу. Зная mi (масса i вещества), запишем полное приращение полной удельной внутренней энергии в ячейке ( $$\Delta _{\varepsilon }$$ ) и приравняем его к уже полученному
Для определения изменения количества каждого газа нужно сделать некое правдоподобное предположение, например, считать, что все $$\Delta \varepsilon _{ikl}$$ одинаковы. Тогда $${N} \Delta \varepsilon_{ikl} = \bar{\varepsilon}_{kl}^{n} - \varepsilon_{kl}^{n}$$ и, соответственно, $$\bar{\varepsilon}_{ikl} = \varepsilon_{ikl}^{n} + \Delta \varepsilon_{ikl}$$. На этом первый этап расчета (предиктор) закончен.
Рассмотрим второй этап расчета.
Движение частиц описывается
которые могут быть приближены, например, с использованием явного
где скорости частиц $$$ \tilde{u}_k, \tilde{v}_k $$$ определяются интерполяцией величин $$$ \bar{u}_{1}, \bar{u}_{2} $$$ в ячейках, окружающих p частицу.
После этого рассчитывается перенос массы и вычисляется новая масса каждой ячейки. Для этого выделяются три группы частиц:
n + 1 слой в пределах ячейки, которые, очевидно, не вносят изменений в массу, импульс, энергию новой ячейки, т.е.$$(X_p^{n}, Y_p^{n}) \in \omega_{ij}, (X_p^{n + 1}, Y_p^{n + 1}) \in \omega_{ij},$$
где $$\omega _{ij}$$ — обозначение старой ячейки,На шаг по времени накладывается ограничение
$$$ {\tau} < \frac{h}{{\sqrt {u_1^2 + u_2^2}}}, $$$что означает запрет на перемещение частицы за один шаг больше, чем на одну
ячейку. Перемещение в данном методе возможно только в соседнюю ячейку. Предположим, что каждая p частица, перешедшая на n + 1 шаге по времени в соседнюю ячейку, переносит в нее массу mp. Это означает, что значение массы mikl считается путем сложения масс всех частиц типа i, для которых $$(X_p^{n + 1}, Y_p^{n + 1}) \in \omega_{ij}$$
Процедура вычисления импульса выполняется следующим образом. Компоненты
полного импульса частиц в ячейке ( k, l ) могут быть вычислены, как $$m_{kl}^{n} \overline{u}_{1, kl}, m_{kl}^{n} \overline{u}_{2, kl}$$, при этом p частица, покинувшая ячейку $$\omega _{ij}$$, уносит импульс $$m_p \overline{u}_{1, {kl}}, m_p u_{2, kl}$$. Изменение импульса в Skl за один шаг по времени будет
здесь символы суммирования означают суммирование по частицам ( p ),
покинувшим данную ячейку и пришедшим в нее, соответственно. После вычисления компонентов импульса каждой ячейки вычисляются компоненты скорости $$u_{1kl}^{n + 1}, u_{2kl}^{n + 1}$$.
Частица типа i, переходящая из Ske в другую
ячейку, переносит полную энергию
Тогда можно вычислить энергию i вещества в ячейке $$\omega _{kl}$$ на промежуточном шаге:
При t = tn + 1 полная удельная энергия изменится на величину
где знаки суммирования снова означают суммы по всем частицам, покинувшим ячейку $$\omega _{kl}$$ и пришедшим в нее, соответственно. Далее получим
$$$ \varepsilon_{ikl}^{n + 1} = \frac{{h^2}}{{m_{ikl}^{n + 1}}}E_{ikl}^{n + 1} - \frac{1}{2} \left[{(u_{1 {kl}}^{n + 1} )^2 + (u_{2 {kl}}^{n + 1} )^2}\right]. $$$Отметим недостатки этого метода. Во - первых, это дискретность плотности, что приводит при небольшом количестве частиц в ячейке к скачкообразным изменениям плотности. Во - вторых, аппроксимация исходных уравнений достигается при количестве частиц, стремящемся к бесконечности. Кроме того, этот метод требует значительно большего количества памяти, чем конечно - разностные методы, так как наряду с физическими характеристиками узлов необходимо хранить и свойства частиц в ячейках.
Подробное описание
Выпишем нелинейную систему уравнений одномерных движений идеальной сжимаемой жидкости в случае баротропных процессов. Она состоит из уравнения Эйлера
$$$ \frac{{{\partial} u}}{{{\partial} t}} + u \frac{{{\partial} u}}{{{\partial} x}} + \frac{1}{\rho} \frac{{{\partial} p}}{{{\partial} x}} = 0, $$$уравнения неразрывности
$$$ \frac{{{\partial}{\rho}}}{{{\partial} t}} + \rho \frac{{\partial}u} {{\partial}x} + u \frac{{{\partial}{\rho}}}{{\partial}x} = 0 $$$и условия баротропности
$$p = f(\rho ).$$Уравнения позволяют определить плотность $$\rho$$ и скорость u в зависимости от координаты x и времени t. Система не имеет решений, зависящих только от $$x \pm a_{0}t$$, но оказывается возможным найти решение этой системы, представляющее собой плоскую волну и являющееся обобщением решений вида $$f(x \pm a_{0}t)$$. Будем искать такие решения системы, для которых скорость u является функцией только плотности $$\rho.$$
$$$ u = \pm \int {\sqrt {\frac{{{dp}}}{{d {\rho}}}} \frac{{d \rho }}{\rho}}. $$$
Обозначим $$$ \frac{dp}{d{\rho}} = a^2 ({\rho}) $$$ и введем величину c = u + a. Какой физический смысл имеет величина c?
$$$ \frac{{\partial}\rho}{{\partial}T} + c(\rho ) \frac{{\partial}\rho}{{\partial}x} = 0, $$$
$$$ c \left({\rho}\right) = \sqrt {A \gamma } \left[{1 + \frac{2}{{\gamma - 1}}}\right] {\rho}^{{{1 \over 2}}(\gamma - 1)}, $$$
$$u = u_0 f({x/t}), {\rho} = \rho_0{\varphi}({x/t})$$.
Течения подобного типа — частный случай автомодельных течений, когда решение зависит от некоторой комбинации независимых переменных. Проиллюстрировать численными расчетами особенности распространения центрированных волн Римана.
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.