Введение в методы численного решения уравнений газовой динамики
4.5. Метод частиц в ячейках Харлоу (PIC method:Particle - In - Cell)
Метод PIC разработан Харлоу в Лос - Аламосской лаборатории (США) в 60 - х годах прошлого века для расчета процессов с большими деформациями исходной области интегрирования (расплескивание, разрушение).
Область интегрирования покрывается фиксированной в пространстве расчетной сеткой, шаг которой h постоянен по обеим координатам x, y, ячейки занумерованы двумя индексами k, l.
В центре ячейки вычисляются величины (компоненты скорости газа),
где i — номер вещества.
— удельная внутренняя энергия газа с номером i,
— масса этого вещества. Если
этого вещества в ячейке нет, то в ней и энергия, и масса полагаются равными нулю.
Предположим, что в каждой ячейке содержится несколько частиц (5 - 10), каждая из которых характеризуется координатами массой
— номер вещества, из которого состоит частица с номером j.
Шаг численного интегрирования состоит в расчете величин и
на верхнем временном слое tn + 1 по вычисленным величинам
на нижнем слое 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*}](/sites/default/files/tex_cache/0ad899d6fe0108cb990096c1af890310.png)
где
![$ E = {\rho}e = {\rho}\left[{\varepsilon + \frac{1}{2}(u_1^2 +
u_2^2 )}\right]. $](/sites/default/files/tex_cache/73519f2a9a0df158c1bfeaf2bb83bfce.png)
В расчетах участвуют также уравнения состояния для каждого газа
![p_{i} = F_{i}(\varepsilon _{i}, \rho _{i}).](/sites/default/files/tex_cache/d09662bb315136dbc72ab1c287830175.png)
На втором этапе аппроксимируются конвективные члены
![\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*}](/sites/default/files/tex_cache/2210c322ac915af70843bae61d5531ae.png)
Опишем вычислительную процедуру на первом этапе. Известны ,
, mn,
, Xn, Yn (остальные индексы для простоты изложения опускаются). Сначала рассчитывается давление
, исходя из предложения равенства давлений на границе двух сред p1 = p2 = ..., или
![F_1 (\varepsilon_1^{n}, \gamma_1^{- 1} m_1^{n} ) = F_2 (\varepsilon_2^{n}, \gamma_2^{- 1}m_2^{n} ) = \ldots](/sites/default/files/tex_cache/54e21e74bd25877e636ba31c181eef92.png)
К этим уравнениям добавляется условие поскольку
— часть объема h2 ячейки, занимаемого газом с номером i. По известной массе i газа находим его плотность:
, а по известной удельной энергии
- давление
.
Затем по закону Дальтона находится давление, , которое приписывается к центру ячейки k, l. Система нелинейных алгебраических уравнений решается, вообще говоря, итерационным методом. В случае
выписывается ее явное решение. Далее находим предварительные значения расчитываемых величин, которые обозначим как
. Первое из уравнений
![\frac{{{\partial}\rho }}{{\partial}t} = 0](/sites/default/files/tex_cache/fab95aaffbd591dcf01d6967827fdb89.png)
![\bar{m}_{ikl} = m_{ikl}^{n}](/sites/default/files/tex_cache/3ef73bae38e036d13d65cba54fe09099.png)
![\frac{{\partial}{\rho}}{{\partial}t} = 0](/sites/default/files/tex_cache/86ff281c6f964a50e0354b55b6f1e12c.png)
![\begin{gather*}
\rho_{kl}^{n} \frac{{\bar{u}_{1{kl}} - u_{1{kl}}^{n}}}{\tau} + \frac{{p_{k + 1/2, l}^{n} - p_{k - 1/2, l}^{n}}}{h} = 0, \\
\rho_{kl}^{n} \frac{{\bar{u}_{1{kl}} - u_{2{kl}}^{n}}}{\tau} + \frac{{p_{{k, l} + 1/2}^{n} - p_{{k, l} - 1/2}^{n}}}{h} = 0. \end{gather*}](/sites/default/files/tex_cache/6c54f1b4e9ccb0ef7e915a1e8a7e1557.png)
Здесь
![\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*}](/sites/default/files/tex_cache/76746c2b092ff5f9656cd7b0c3ac22fc.png)
Последняя из рассчитываемых величин — энергия. Дискретный аналог уравнения энергии в методе частиц в ячейках будет
![\begin{gather*}
\frac{{\bar{E}_{kl} - E_{kl}^{n}}}{\tau} + \frac{{p_{k + 1/2, l}^{n} \cdot u_{1 , k + 1/2, l}^{{n} + 1/2} - p_{k - 1/2, l}^{n} \cdot u_{1 , k - 1/2, l}^{{n} + 1/2}}}{h} + \\
+ \frac{{p_{{k, l} + 1/2}^{n} \cdot u_{2 , {k, l} + 1/2}^{{n} + 1/2} - p_{{k, l} - 1/2}^{n} \cdot u_{2{k, l} - 1/2}^{{n} + 1/2}}}{h} = 0. \end{gather*}](/sites/default/files/tex_cache/8042f7c2b340f3353ca0ec15c28b7376.png)
Здесь
![\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*}](/sites/default/files/tex_cache/d758e57f7128584af8ef66beede5b1a2.png)
Вычислим величину — энергию. Напомним,
что
![$ e = {\rho}(\varepsilon + \frac{{u_1^2 + u_{2}^{2}}}{2}). $](/sites/default/files/tex_cache/f6a52b6747d0c5d51ea3185125b808cd.png)
![$ {\rm E}_{kl}^{n} \cdot h^2 = \left[{(h^2 \cdot {\rho}) \cdot
\varepsilon + (h^2 {\rho}) \frac{{u_1^2 + u_2^2}}{2}}\right]_{kl}^{n}, $](/sites/default/files/tex_cache/f319c28c40c03c5bfa69f840014d5db5.png)
где — масса ячейки, равная
. В таком случае, с учетом закона сохранения массы,
![$ \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]. $](/sites/default/files/tex_cache/6f17cf54b2dea37942088e37d4dcd763.png)
Вычислим величину , имеющую смысл полной внутренней энергии в ячейке, зная массу
и удельную внутреннюю энергию
вещества.
![\left[{(h^2 {\rho}) \varepsilon }\right]_{kl}^{n} = \sum\limits_i {m_{ikl}^{n}} E_{ikl}^{n},](/sites/default/files/tex_cache/6ca52868d3898f593e397981086d549d.png)