Диссертация (1145359), страница 3
Текст из файла (страница 3)
Суть его можно свестик следующему. Вместо точного значения какой-либо газодинамической величины f (r) используют ее сглаженное значение < f (r) >, которое определяют припомощи интегрального интерполированияZ< f (r) >= f (x) W (r − x; h) dx ,(1.1)По существу это есть процедура свертки функции f (r). Предполагается, что интерполяционное ядро W (r; h) нормировано на 1 и стремится к дельта-функциипри h → 0.Если функцию W (r) выбрать сферически-симметричной, то точностьпредставления величины f (r) ее сглаженным значением будет O(h2 ) (см., например, [32]).
Следующий шаг сводится к оценке многомерного интеграла ввыражении (1.1) методом Монте-Карло. Если значения f (r) известны для каких-либо N точек, то < f > (ri ) оценивается как< f (ri ) >=NXj=1f (rj )W (ri − rj ; h) ,< n(rj ) >(1.2)где n(r) — плотность распределения выбранных точек.Газодинамическое течение можно описывать как ансамбль движущихсяэлементов газа.
При численном моделировании мы можем выбрать только конечное число элементов (N ), но чем больше число элементов в ансамбле и чемменьше их размеры, тем ближе такое описание к описанию непрерывной среды. В методе SPH элементы газа представляются частицами конечного размера. Положение частиц изменяется согласно уравнению движения; каждой из15них соответствует значение температуры газа в данной точке; их скорость естьлокальная скорость течения; распределение же плотности частиц n(r) дает распределение плотности газа ρ(r).
При таком подходе с учетом свойств заданнойфункции W (r; h) оценка сглаженных значений гидродинамических величин методом Монте-Карло (1.2) есть не что иное, как процедура сглаживания по ансамблю элементов газа в объеме размером порядка h3 ; h называется длинойсглаживания и характеризует размер частиц в ансамбле. Если каждой частицеприписать массу mj так, чтоPmj — полная масса газа, то сглаженное значе-ние плотности ρi в данной точке ri определяется суммой масс частиц mj с весомW (ri − rj ; h) из окрестности размером порядка hXρi =mj W (ri − rj ; h) .(1.3)jЗначение производной от f (r) оценивают при помощи процедуры сглаживания (1.1) самой этой величины с ядром ∇W (r − x; h)Z< ∇f (r) >= f (x) ∇W (r − x; h) dx .(1.4)Описанный формализм, будучи примененным к системе гидродинамических уравнений, сводит их к обыкновенным дифференциальным, решать которые значительно легче, чем уравнения в частных производных.
Если ядро взятьдостаточно компактным, например, в виде сплайна, как это предложено в [36]1 − 1.5 q 2 + 0.75 q 3 , 0 6 q 6 1 ,1W =0.25 (2 − q)3 ,1 6 q 6 2,A 0,q > 2,(1.5)где q = r/h, то суммирование в (1.2) нужно производить лишь по небольшомучислу соседей вокруг данной частицы в окрестности размером 2 h. Для трехмерного случая A = π h3 .Строгое обоснование описанного метода и его улучшенных модификаций,а также обсуждение многих идейных и философских сторон SPH можно найти,например, в [32, 36–38].16В последние годы в развитии метода SPH достигнут большой прогресс, иэто дает ему возможность уверенно конкурировать с конечно-разностными алгоритмами.
В первую очередь, такой прогресс связан с использованием переменной длины сглаживания для каждой частицы — h(ri ; t) [32], что значительнорасширяет динамический диапазон пространственного разрешения и позволяеткорректно моделировать эволюцию объектов с быстро меняющейся структуройи плотностью. Во-вторых, благодаря введению дифференцированого шага повремени при интегрировании уравнений, преодолевается глобальное ограничение на временной шаг, которое накладывается условием Куранта [32].
При этомзатраты машинного времени существенно уменьшаются. Таким образом, длямоделирования объектов с произвольной и далекой от симметрии структурой,какими, например, являются формирующиеся в результате взаимодействия галактик приливные и кольцеобразные детали, SPH метод оказывается наиболееподходящим, если только решаемая задача не требует высокоточного расчетапрофилей возникающих ударных волн. В противном случае предпочтительными являются конечно-разностные методы.1.1.2.
Численная реализацияОсновные уравненияС точки зрения программирования SPH алгоритм довольно прост. При использовании машин типа рабочих станций эффективность программ, основанных на данном методе, достаточно велика (см., например, [38]). Для адаптацииметода к вычислительным возможностям персональных компьютеров я былавынуждена выбрать один из наиболее простых вариантов SPH с постояннойдлиной сглаживания, который и описан в начале раздела.При моделировании газовых течений во взаимодействующих галактикахвполне оправданным считается изотермическое приближение [39, 40]. Фактически, метод SPH из-за ограничения на разрешение по массе не может одинако-17во корректно описывать различные фазы межзвездной среды.
Поэтому, чтобыподавить образование плотных облаков, функцию охлаждения приходится обрезать для температур ниже ≈ 104 K. С другой стороны, время высвечиванияобычно меньше динамического, и температура газа остается постоянной и близкой к этому пределу [40, 41].В систему уравнений газодинамики, описывающих изотермическое течениегаза во внешнем гравитационном поле, входят:уравнение неразрывностиdρ= ρ ∇v ,dt(1.6)dv1= − ∇P − ∇Φ ,dtρ(1.7)уравнение движениягде ρ — плотность газа, P — давление, а Φ — гравитационный потенциал, изамыкает систему уравнение состояния идеального газа. Для изотермическогослучаяP = c2 ρ ,(1.8)где c = const — скорость звука.
При движении во внешнем поле Φ являетсязаданной функцией.Существуют различные варианты перехода от гидродинамических уравнений (1.6) – (1.8) к SPH уравнениям (подробнее об этом см. в [37]. Я остановиласьна следующей системе:ρi = mXW (ri − rj ; h) ,(1.9)jdri= vi ,dtdvi1= − ∇Pi − ∇Φi ,dtρi2Pi = c ρ i .Уравнение (1.9) записано для частиц одинаковой массы m.(1.10)(1.11)(1.12)18Выражение для гидродинамического ускорения можно представить в виде√√−2 Pi (∇ Pi )/ρi . Тогда процедура сглаживания (1.1) и (1.4) с учетом (1.2)и (1.12) приводит к−X 2 c2∇Pi= −m∇i W (ri − rj ; h) .√ρiρρijj(1.13)При больших числах Маха давление не способно воспрепятствовать пересечению орбит частиц. В такой ситуации в действие вступает молекулярнаявязкость.
При численном моделировании вводится ее аналог — искусственнаявязкость. Ускорение за счет сил искусственной вязкости можно записать какaivisc= −mXQij ∇i W (ri − rj ; h) ,(1.14)jгде Qij — вклад вязкости в градиент давления. Существуют различные формыпредставления Qij .
Использовалось два следующих выражения (о достоинствахи недостатках такого выбора см. [32]). Одно из них может быть таким [37] (−α c µ + β µ2 )/ρ , (v − v )(r − r ) 6 0 ,ijijijijijQij = 0,(vi − vj )(ri − rj ) > 0 ,(1.15)2+η 2 ), ρij = (ρi +ρj )/2, rij = |ri −rj |, η ' 0.1 h. Пагде µij = h (vi −vj )(ri −rj )/(rijраметры α и β — аналоги коэффициентов вязкости в уравнении Навье-Стокса.Обычно выбирают α ' 1, β ' 2.Можно также ввести искусственную вязкость, зависящую от дивергенцииполя скоростей [32], q /ρ + q /ρ , (v − v )(r − r ) 6 0 ,i ijjijijQij = 0,(vi − vj )(ri − rj ) > 0 ,(1.16)где α c h|∇ · v | + β h2 |∇ · v |2 , ∇ · v 6 0 ,iiiqi = 0,∇ · vi > 0 .PСглаженная оценка ∇ · vi = −(m/ρi ) j (vi − vj ) ∇i W (ri − rj ; h).(1.17)В свою численную схему я включила оба выражения для вязкости и предусмотрели возможность переключения с одного на другое.19Вычислительная схемаДля решения уравнений (1.9) – (1.12) использовалась явная схема с перешагиванием (leap frog), обеспечивающая второй порядок точности,(n+1/2)= ri(n−1/2)(n+1)= viri(n)vi(n),(1.18)(n+1/2).(1.19)+ δt vi+ δt aiТак как ускорения в (1.19) зависят от скорости через искусственную вязкость,то для сохранения второго порядка точности вычисление скорости осуществлялось в два этапа.
Сначала в процедуре, определяющей новые положения частиц,делалась предварительная оценка скорости(n+1/2)(n)vi(n+1/2)Затем новые координаты ri(n+1/2)этого через ri(n+1/2), ρi= vi+δt (n−1/2)a.2 i(n+1/2)использовались для вычисления ρi(n+1/2)и vi(n+1/2)находилось ускорение ai(n+1)в соответствии с (1.19) определялось значение vi(0)Если заданы начальные значения ri(0)= ri(0)+ δt/2 vi.(0)(1/2)(1/2). И наконец,и vi , то для первого шага уравне-ние (1.18) неприменимо. В этом случае оценка riследует из ri, послевторого порядка точности(0)+ 1/2 (δt/2)2 ai .Выбор шага интегрированияТак как для решения уравнений (1.9) – (1.12) использовалась явная схема,то выбор временного шага ограничен условием Куранта (см., например, [32])∆t =0.3 h,bi + c + 1.2 (α c + β di )(1.20)bi = h|∇ · vi |.
Для искусственной вязкости (1.15) di = maxj |µij |, а для (1.16)di = h|∇ · vi |, если ∇ · vi < 0 и di = 0, если ∇ · vi ≥ 0. Шаг интегрирования δt выбирался не больше, чем ∆t, но так, чтобы δt = ∆t0 /2n , где ∆t0 —максимально возможный шаг (задаваемый параметр), а n > 0 — целое. Если20шаг интегрирования менялся в процессе счета, то для сохранения устойчивостисхемы новые координаты определялись из(n+1/2)ri(n−1/2)= ri+δtold + δtnew (n) δt2old − δt2new (n−1/2)vi −ai.28Алгоритм суммированияДля ядра в виде сплайна (1.5) суммирование в (1.9) и (1.13) нужно производить лишь по небольшому числу соседей, расположенных в окрестностиразмером ∼ 2 h от данной частицы.
Если длина сглаживания постоянна, тоодин из быстрых путей поиска соседей — использование сетки с длиной ячейки2 h и алгоритма связанных списков [42]. При реализации этот путь не всегдабывает эффективным из-за больших требований к объему оперативной памяти.Тем не менее в описываемой схеме я остановилась на данном алгоритме, так какон наиболее прост в реализации.Критерием выбора h было условие достаточно большого среднего числачастиц (порядка 20–30), по которым нужно производить сглаживание гидродинамических величин.Начальные и граничные условияВ начальный момент времени частицы распределяются в пространстве всоответствии с задаваемым профилем плотности. Начальные значения плотности определяются согласно процедуре сглаживания (1.9).
Что касается скоростей, то для уменьшения влияния флуктуаций, вводимых в систему процедуройслучайного задания данных, применялся стандартный способ [32] сглаживанияполя скоростей < vi > (0) = (m/ρi )P(0)j(0)vj W (ri(0)− rj ; h).В описываемой схеме предполагаются свободные граничные условия, т.е.давление на границе распределения газа считается пренебрежимо малым. Такоеприближение оправдано, если моделируется течение холодного газа.211.1.3. ТестыБыло проведено два тестовых эксперимента. В первом — моделировалосьодномерное течение адиабатического газа. Для этого в вычислительную схемубыло добавлено уравнение энергии du/dt = −P/ρ∇v, где u — удельная тепловая энергия, а уравнение состояния записано в виде P = (γ − 1)ρu, γ —показатель адиабаты. Решалась так называемая задача Сода о распространении ударной волны в трубе.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
