CC BY
05. Telyakovskii S.A. On the behavior near the origin of the sine series with convex coefficients // Publ. Inst. Math. (Beograd). 1995. 58(72). 43-50.
6. Солодов А.П. Точные константы в двусторонней оценке С.А. Теляковского суммы ряда по синусам с выпуклой последовательностью коэффициентов // Матем. заметки. 2020. 107, № 6. 906-921.
7. Попов А.Ю., Солодов А.П. Оптимальные на отрезке [п/2,п] двусторонние оценки суммы синус-ряда с выпуклой последовательностью коэффициентов // Матем. заметки. 2022. 112, № 2. 317-320.
8. Попов А.Ю. Уточнение оценок сумм синус-рядов с монотонными и косинус-рядов с выпуклыми коэффициентами // Матем. заметки. 2021. 109, № 5. 768-780.
Поступила в редакцию 30.08.2023
УДК 519.633.6
О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ КОЛЕБАНИИ В ХОЛОДНОЙ, НО ВЯЗКОЙ плазме
Е. В. Чижонков1
Численно анализируется влияние вязкости на нерелятивистские колебания холодной плазмы. С этой целью построена неявная разностная схема типа Мак-Кормака, имеющая более слабое ограничение на устойчивость, чем явная схема, и реализуемая без итераций, что увеличивает ее вычислительную эффективность в десятки раз. Показано, что учет вязкости плазмы может быть причиной не только затухания амплитуды плазменных колебаний, но и изменения формы решения. При увеличении коэффициента вязкости у решения наблюдается седловая точка, которая сохраняется во времени.
Ключевые слова: численное моделирование, нерелятивистские колебания, холодная вязкая плазма, неявная схема Мак-Кормака, седловая точка.
The effect of viscosity on non-relativistic oscillations of cold plasma is numerically analyzed. For this purpose, an implicit difference scheme of the MacCormack type is constructed, which has a weaker restriction on stability in comparison with the explicit scheme. The scheme is implemented without iterations, which increases its computational efficiency tenfold. It is shown that taking into account the plasma viscosity can cause not only attenuation of the amplitude of plasma oscillations, but also a change in the shape of the solution. With an increase in the viscosity coefficient, a saddle point is observed in the solution, which is preserved in time.
Key words: numerical simulation, non-relativistic oscillations, cold viscous plasma, implicit MacCormack scheme, saddle point.
DOI: 10.55959/MSU0579-9368-1-65-4-5
Введение. Гидродинамическая модель холодной плазмы хорошо известна и достаточно подробно описана в учебниках и монографиях по физике плазмы [1-4]. В настоящее время внимание к этой модели обусловлено в первую очередь задачами, относящимися к распространению сверхмощных лазерных импульсов в плазме [5, 6]. Подобные постановки напрямую связаны с приложением результатов, удостоенных Нобелевской премии по физике 2018 г. Приведем следующие примеры практически важных задач этой тематики: лазерное ускорение электронов и ионов, быстрое зажигание термоядерного синтеза, ядерные реакции в луче лазера, синхротронное и субмиллиметровое излучение и пр. [7]. Численному моделированию колебаний в холодной плазме, а также кильватерных волн, возбуждаемых коротким мощным лазерным импульсом, посвящена монография [8].
Напомним, что гидродинамическая модель "холодной" плазмы, в которой температура электронов формально полагается равной нулю, является точным математическим следствием кинетической модели, основанной на системе уравнений Власова-Максвелла (см., например, [1, 3]). Однако такое гидродинамическое приближение является безусловной идеализацией физических процессов,
1 Чижонков Евгений Владимирович — доктор физ.-мат. наук, проф. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: [email protected].
Chizhonkov Evgenyi Vladimirovich — Doctor of Physical and Mathematical Sciences, Professor, Lomonosov Moscow State University, Faculty of Mechanics and Mathematics, Chair of Computational Mathematics.
© Чижонков Е. В., 2024 © Chizhonkov E. V., 2024
так как температура плазмы может быть необязательно большой, но необходимо положительной. По этой причине в работе [9] в модель холодной плазмы было внесено уточнение на случай малых температур, т.е. "теплой" плазмы. Главный смысл модификации модели — добавление в уравнение, выражающее закон сохранения импульса, новых слагаемых, описывающих электронное давление, вязкость и электрон-ионные соударения.
Влияние различных зависимостей давления от электронной плотности было рассмотрено в [10], причем был сделан вывод, что учет давления трансформирует ленгмюровские колебания в бегущие волны. В свою очередь влияние электрон-ионных соударений было теоретически исследовано в [11] и численно в работе [12], где было показано затухание амплитуды ленгмюровских колебаний при увеличении коэффициента соударений. Таким образом, вне области внимания оставалось только качественное влияние вязкого слагаемого на решение исходной дифференциальной задачи.
Однако здесь следует отметить важное обстоятельство: для численного решения систем гиперболического типа, к которым относятся уравнения холодной плазмы, традиционно используются явные разностные схемы (см., например, [8, 13] и цитируемую там литературу), а учет вязкости сразу же приводит в них к жестким ограничениям на устойчивость. Поэтому введение в модель слагаемого, описывающего вязкость, диктует необходимость разработки новых или соответствующих модификаций ранее известных численных алгоритмов (с более слабыми ограничениями на устойчивость вычислений). В настоящей работе с целью моделирования ленгмюровских колебаний при учете вязкости плазмы построена новая разностная схема сквозного счета типа неявной схемы Мак-Кормака [14]. Ее главной особенностью является экономичность: замена классической явной схемы Мак-Кормака на неявную в численном алгоритме в данном случае требует незначительного увеличения вычислительной работы (в пределах 10%), причем неявная процедура носит безытерационный характер, что существенно увеличивает ее привлекательность. Требования на устойчивость в неявной схеме качественно ослаблены, что достигается за счет ухудшения асимптотической погрешности аппроксимации. Однако при доминирующем переносе (т.е. достаточно малых значениях коэффициента вязкости, как это имеет место в "теплой" плазме) этим ухудшением можно пренебречь.
Настоящая работа имеет следующую структуру. В первом пункте приведена простейшая (плоская, пространственно одномерная) постановка задачи, описывающая нерелятивистские колебания холодной плазмы с учетом вязкости в эйлеровых переменных. Во втором пункте приведена оригинальная неявная схема Мак-Кормака для модельной линейной скалярной задачи типа уравнения конвекции-диффузии, а также отмечены наиболее важные вычислительные свойства этой схемы. В третьем пункте построена неявная схема Мак-Кормака для нелинейных уравнений с вязкостью из первого пункта. Наконец, в четвертом пункте описаны результаты вычислительных экспериментов, иллюстрирующих влияние вязкого слагаемого на форму ленгмюровских колебаний, порождаемых коротким мощным лазерным импульсом. В заключении систематизированы результаты проведенных исследований.
1. Постановка задачи для колебаний вязкой плазмы. Будем считать плазму нерелятивистской электронной жидкостью, пренебрегая рекомбинационными эффектами и движением ионов. Тогда в рамках модели холодной плазмы ее плоские одномерные колебания можно описать безразмерной системой уравнений
где V — скорость электронов, Е — электрическое поле, р и в — обезразмеренные координаты по пространству и времени соответственно. К системе (1), (2) обычно добавляют уравнение
ЩР,0) = 1-«, (3)
характеризующее безразмерную плотность электронов N. Формула (3) является частным случаем теоремы Гаусса [15], которая в дифференциальной размерной форме имеет вид ё1у Е = 4 п е(п — п0). Здесь е < 0 — заряд электрона, по — значение невозмущенной электронной плотности. Подробный вывод уравнений (1)-(3) можно найти в различных источниках (см., например, [8, с. 19]).
Система из двух первых уравнений относится к гиперболическому типу и хорошо изучена как аналитически (см., например, [16, 17]), так и численно [18]. Для таких систем существует локально
по времени единственное решение задачи Коши того же класса, что и начальные данные. Также известно, что для таких систем потеря решением гладкости происходит по одному из следующих сценариев: либо сами компоненты решения в течение конечного времени обращаются в бесконечность, либо они остаются ограниченными, но в бесконечность обращаются их производные [19]. Последняя возможность реализуется, например, для однородных законов сохранения, к которым относятся уравнения газовой динамики, где возникновение особенности соответствует образованию ударной волны.
Как уже говорилось выше, в работе [9] в модель холодной плазмы было внесено уточнение на случай малых температур. В используемых здесь обозначениях модификация модели может быть выражена как добавление в правую часть уравнения (1) новых слагаемых:
9У_ дУ_ }_дР_
V д^
--2 - чУ,
где Р — электронное давление; V ^ 0 и п ^ 0 — постоянные коэффициенты вязкости и соударений, зависящие от температуры. С учетом специфического интереса к влиянию вязкости на решение рассмотрим модификацию уравнения (1), содержащую только одно новое слагаемое:
дУ дУ_ у д2У
V > 0.
(4)
Формализуем постановку задачи: нас будет интересовать в полуплоскости {(р, в) : р £ М, в > 0} решение задачи Коши для уравнений (2)-(4) с начальными условиями
V(р, 0) = Vо(р), Е(р, 0) = Ео(р), р £ М.
(5)
Наиболее естественным выбором начальных условий (5) является имитация возмущений электрического поля, которые порождаются в разреженной плазме коротким мощным лазерным импульсом при его фокусировке в линию (этого можно добиться, используя цилиндрическую линзу, см. детали в [20]):
ЫР) = (^)2реХР{"2р|}' Мр) = °> (6)
где а*,р* — амплитуда и ширина импульса.
2. Неявная схема Мак-Кормака для модельной постановки. Следуя оригинальной идее [14], напомним основные формулы неявной схемы Мак-Кормака для простейшего уравнения конвекции-диффузии с постоянными коэффициентами а и V ^ 0:
ди ди д2и _
яГ +а7Г = 1У7Г2> жеК> д1 дх дх2
снабженного начальным условием и(х, 0) = и0(х), таким, что и0(х) = 0 при \х\ ^ ж.
Для построения и анализа численных методов решения задачи Коши (7) полезной является формула [21]:
u(x,t) = j uq (х — at + 2£\/vtj ехр {— £2}
Определим дискретизацию независимых переменных с помощью постоянных параметров т и h так, что
tn = пт, п ^ 0; Xi = ih, i = 0, ±1, ±2,...,
и будем обозначать зависимую переменную u(x,t) в узле сетки (xi,tn) через и'П-
Запишем неявную схему Мак-Кормака в традиционной форме "предиктор-корректор", считая известными величины и'П:
1) предиктор с результатом Up реализуется формулами
пт VT
Д< = -— К+1 - и?) + F К+1 - 2и? + uU) ,
\T \ \т (8) '^'I^Pa^ .'^'^P Р Tí Г- Р
1 + ^)6ирг=Аи? + —6ирг+1, ирг=и? + 5ир,
где вычисления проводятся в сторону уменьшения индекса: г = ...,к + 1,к, к — 1,... 2) корректор с результатом иС реализуется формулами
^р = -^(ир-ир_1)+^(ир+1-2ир+ир_1),
\т\ \т (9)
1 + — ) 6и1 = Аир + — 6иЧ_!, и1 = и™ + 5и1,
где вычисления проводятся в сторону увеличения индекса: г = ...,к — 1,к,к + 1,....
В формулах (8), (9) верхний индекс р (или с) обозначает шаг предиктор (или корректор), а п — временной слой Ьп; X — постоянный параметр схемы, который будет определен ниже.
Окончательная формула, формирующая решение на следующем временном слое с номером (п + 1), имеет вид
иП+ = « + иС) /2.
Следует отметить, что при X = 0 неявная схема трансформируется в обычную явную схему Мак-Кормака [22] с условием устойчивости
1
Т ^ \а\/Н + 2ту//г2 (10)
и первым дифференциальным приближением [23] вида
т2 3 аЬ2 д3и иЬ2 д4и , 2 , 2.
Щ = Ьои ~ ~6Ь°и + + °(Т + Л
где Ь0и = иихх — аих.
Если же выбирать X из условия [14]
1 Г, , 2и Ь
А ^ - тах < |а| + ----, 0
2 [ Ь т
то неявная схема (при X > 0) будет безусловно устойчивой, а в правую часть первого дифференциального приближения добавятся слагаемые [24]
т 2
, , , h \ д4и . Л h\ д3 и
Иными словами, в неявной схеме на гладких решениях качественное ослабление требований на устойчивость достигается за счет ухудшения асимптотической погрешности аппроксимации. Однако
при доминирующем переносе (т.е. достаточно малых V, как это имеет место в "теплой" плазме) этим ухудшением можно пренебречь.
3. Неявная схема Мак-Кормака для вязкой плазмы. Приведем систему (2), (4) к удобной в рассматриваемом случае векторной форме:
где оператор А^) является линейным и диагональным: А = V(р, в) I, I — единичная (2х2)-матрица; В — диагональная (2 х 2)-матрица с элементами Ьц = V, Ь22 = 0; и = ^,Е)Т, Б = (—Е^)т — вектор-функции, рассматриваемые в полуплоскости {(р, в) : в ^ 0, р £ М}.
Определим дискретизацию независимых переменных с помощью постоянных параметров т и Н так, что
вп = пт, п ^ 0; рг = гН, г = 0, ±1, ±2,...,
и будем обозначать зависимую переменную И(р, в) в узле сетки (рг,вп) через ИП.
Введем полезные обозначения для операторов разностей "вперед" В+ и "назад" , у которых аргумент может быть как векторным, так и скалярным:
В+¥г = ¥г+1 - ¥г, Б-¥г = ¥г - ¥г-ъ
и запишем неявную схему Мак-Кормака (при Л = 1) для системы (11), считая известными величины Vn Еп Nп:
'г > г ' 1 Ч '
1) предиктор с результатом И:
положим V+-y|2 = (^п + /2, определим матрицу
^1/2 = ^(|С+1/21 )+ШВ
и последовательно вычислим
Аи? = ~А (*£1/2) + Д+2Г Ц? + тв™
г
т
/-Л-СТ!н1/21?+]ЛЗ? = Ди?> (12)
В+ +
2) корректор с результатом И£:
положим V--l|2 = (У[ + /2, определим матрицу
С_1|2 = А(\^1|20 +
2 В
^Н
и последовательно вычислим
(1 + ХТ-С11/2П-)бЩ = А и?,
и = ип + ¿И.
(13)
В формулах (12), (13) верхний индекс р (или с) обозначает шаг предиктор (или корректор), а п — временной слой вп; Л — постоянный параметр схемы, который будет определен ниже.
Окончательные формулы, формирующие решение на следующем временном слое с номером (п + 1), имеют вид
и?- = Н1±Д =
4. Численные эксперименты. Для однозначного задания начальных условий и сохранения преемственности с предыдущими расчетами (см., например, [8, 25]) зафиксируем параметры в (6):
р* = 0.6, а* = 0.414 и заметим, что на больших расстояниях от прямой р = 0 плазменные колебания не возбуждаются:
V(р ±с», в) = 0, Е(р ±с», в) = 0.
Однако в целях численного моделирования плазменных колебаний расчетную область необходимо ограничить: определим ее по переменной р как отрезок [—й,й], на концах которого следует задать искусственные граничные условия. Обсуждению их построения посвящен раздел 3.6 в [8], здесь же мы ограничимся "обрезанием" бесконечной области с помощью однородных граничных условий первого рода:
V (±й,в) = Е (±й,в) = 0.
Отметим, что такой выбор одинаково удобен при использовании как явной, так и неявной разностной схемы типа Мак-Кормака. Конечно, параметр й следует выбирать достаточно большим. В силу экспоненциального затухания функции Ео(р) достаточно положить й = 4.5р*. В этом случае имеем ехр2{—й2/р2} ~ 2.5768 • 10_18. Это означает, что при вычислениях с двойной точностью величина скачка начальной функции Ео в точках р = ±й соизмерима с машинной точностью, т.е. с обычной погрешностью округления данных. Иными словами, при численном моделировании колебаний эффект обрезания начальных условий совершенно не будет заметен, что полностью соответствует понятию "искусственной границы". Учитывая сказанное, параметр й, характеризующий искусственную границу, положим следующим: с! = 4.5 р*. В качестве ограничения расчетной области по переменной в примем отрезок [0, 6п], считая, что трех периодов достаточно для наблюдения за решением.
Напомним, что исходные уравнения обезразмерены так, что плазменная частота шр =
(4пе2п0/т)1/2, где т — масса электрона (остальные обозначения определены ранее), становится равной единице, т.е. период ленгмюровских колебаний составляет 2п. Поэтому процесс колебаний в рамках нерелятивистской гидродинамической модели холодной плазмы при отсутствии вязкости (и = 0) кратко состоит в следующем. Начальное пространственное распределение электронной плотности N, имеющее один глобальный минимум как следствие формул (3) и (6), приводит к избытку положительного заряда в окрестности начала координат (т.е. при р = 0). По этой причине начинается движение электронов в направлении центра области, что через четверть периода вызывает равновесное распределение электронов N(р,п/2) = 1. Однако электроны движутся с ускорением и "проскакивают" положение равновесия, что в середине периода колебаний инициирует распределение плотности с глобальным максимумом также при р = 0. В свою очередь сформировавшийся избыток электронов в центре области приводит к их движению в обратном направлении и еще через половину периода пространственное распределение электронной плотности возвращается к начальному (как при в = 0). В результате наблюдается 2п-периодическое по времени решение в соответствии с аналитическими результатами [16] (подробные доказательства приведены в [17]).
Обратим внимание, что при взятых параметрах электронная плотность N при р = 0 колеблется в диапазоне [0.52,10.96] (границы указаны с двумя верными знаками после десятичной точки), т.е. изменяется на периоде примерно в 20 раз. Для модельных начальных условий (6) соответствующая аналитическая формула
Щр = 0,в) = --!~а а = Е'0(р =
1 — а(1 — еов в) \р*/
следует из [18]. В результате функция Е(р,в), от которой зависит такое значительное изменение плотности (3), является не очень гладкой, что требует выбора достаточно малых значений параметра дискретизации по пространству Ь (при расчетах — порядка 10_3).
Хорошей иллюстрацией к описанию невязких колебаний является рис. 1,а, на котором изображена динамика электронной плотности, и рис. 1,Ъ, где представлены пространственные распределения плотности в некоторые характерные моменты времени (в = 0, п/2, п). При отсутствии вязкости максимальное по области значение Nmax (изображается сплошной линией) синхронизировано со значением в центре области N¡^8 (изображается пунктирной линией): их максимумы одинаковы по величине и совпадают во времени. При этом локальные максимумы по времени одновременно являются локальными максимумами в пространстве. Следует обратить внимание на то, что в силу нелинейности исходных уравнений относительно малые начальные отклонения плотности от равновесного значения приводят к большим возмущениям в моменты времени в = п(2к + 1), к = 0,1,... (см. рис. 1,Ъ). При наличии законов сохранения заряда и энергии [26] это свидетельствует о значительном изменении гладкости решения на каждом периоде колебаний.
0 4 8 12 16 20 в -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 р
Рис. 1. Коэффициент вязкости V = 0: а — динамика плотности (сплошная линия — максимальное по области значение Nmax, пунктирная линия — значение в центре области Naxis); Ь — пространственные распределения плотности (пунктирная линия — график N0 при в = 0, штриховая линия — график Nп/2
при в = п/2, сплошная линия — график N. при в = п)
Если коэффициент вязкости V отличен от нуля, но достаточно мал, например V = 10_3, то качественное изменение решения состоит в затухании амплитуды колебаний, что наглядно представлено на рис. 2. На рис. 2,а , как и на рис. 1,а, изображена динамика максимального по области значения Nmax и значения в центре области N¡^3. Легко заметить, что первый во времени максимум плотности уменьшился примерно в 3.5 раза по сравнению с невязким случаям, причем далее величины этих максимумов монотонно убывают. Однако синхронизация в динамике по-прежнему присутствует, т.е. эффект вязкости пока не может изменить локаций в пространстве и времени локальных экстремумов: максимальное скопление электронов, как и минимальное, наблюдается строго в начале координат. На рис. 2,Ь следует обратить внимание на то, что в моменты времени в = п/2 + пк, к = 0,1,..., равновесное значение плотности N(р,п/2) = 1 не достигается: в окрестности начала координат наблюдается небольшой избыток положительного заряда.
О 4 8 12 16 20 0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 р
Рис. 2. Коэффициент вязкости V = 10~3: а — динамика плотности (сплошная линия — максимальное по области значение Nmax, пунктирная линия — значение в центре области Naxis); Ь — пространственные распределения плотности (пунктирная линия — график N0 при в = 0, штриховая линия — график N^/2
при в = п/2, сплошная линия — график N. при в = п)
Определенный интерес представляет качественное изменение решения, наблюдаемое при увеличении коэффициента вязкости. На рис.3 приведены результаты численных экспериментов при V = 10_2. Динамика электронной плотности сигнализирует, что вязкость может преобразовать локальный максимум по пространству и времени в седловую точку. Об этом свидетельствует несовпадение в моменты в = п(2к + 1), к = 0,1,..., представленных на рис. 3,а графиков (максимальное по области значения Nmax и значение в центре области Naxis). Локальные по времени максимумы плотности при р = 0 не являются локальными максимумами в пространстве. Это хорошо иллюстрируют графики на рис. 3,Ь. Кулоновские силы в центре области не могут изменить направление выпуклости графиков: изначальный "провал" плотности N0 не компенсируется даже приближенным равновесием (график N^2 значительно отличается от константы). В результате формируется
седловая точка (на графике N — "ямка" в окрестности прямой р = 0 ), которая никуда не исчезает с течением времени.
_I_I_I_I_I_I___I_I_I_I_I_I_
0 4 8 12 16 20 9 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 р
Рис. 3. Коэффициент вязкости V = 10 2: а — динамика плотности (сплошная линия — максимальное по области значение Жтах, пунктирная линия — значение в центре области Naxis); Ь — пространственные распределения плотности (пунктирная линия — график N0 при в = 0, штриховая линия — график N^/2
при в = п/2, сплошная линия — график Nп при в = п)
Таким образом, можно зафиксировать, что увеличение коэффициента вязкости не только усиливает монотонное затухание амплитуды колебаний, но и приводит к качественному изменению формы колебаний — возникновению седловой точки, которая сохраняется во времени. Подобный эффект наблюдали и другие авторы [9, 27], но в представленных в этих работах результатах учитывались дополнительно электрон-ионные соударения, хотя и с малыми значениями коэффициентов. В настоящей работе демонстрируется, что причиной как затухания амплитуды колебаний, так и формирования неисчезающей седловой точки может являться только учет вязкости плазмы.
Отметим при этом, что увеличение коэффициента вязкости имеет естественный (физический) предел, связанный с исходным ограничением математической модели — небольшой температурой плазмы [9]. Иными словами, влияние вязкости на колебания холодной плазмы не следует доводить до физического абсурда, когда механизм вязкости будет доминировать над процессом переноса (речь идет о значениях коэффициента вязкости порядка 10"1 и более).
Прокомментируем вычислительную сторону проведенных экспериментов. Как уже говорилось выше, характерной особенностью задачи является необходимость применения достаточно подробных пространственных сеток. Например, чтобы добиться в невязком расчете (V = 0) относительной погрешности для максимума электронной плотности не более 2%, следует выбирать пространственный шаг Н = 0.25 • 10"3 и менее. При этом явная схема (при А = 0) типа Мак-Кормака [25] допускала предельное значение шага по времени т = 10 Н, так как коэффициент |а| = \У(р, в)|тах в условии (10) не превышал 0.1. Для сравнения заметим, что в предложенной выше неявной схеме (при А = 1) шаг по времени ограничений сверху не имел, но при расчетах не превосходил т = 100 Н, чтобы не ухудшить реальную (не асимптотическую!) погрешность аппроксимации. Таким образом, в невязком случае объем вычислений можно сократить как минимум в 10 раз. Уточним, что все представленные выше графики были получены независимо, с использованием каждой из схем и визуально неотличимы.
Рассмотрим пример ограничения шага из-за вязкости. В расчете с V = 10"3 при том же значении Н явная схема типа Мак-Кормака допускала значение шага по времени т = Н/20 (увеличение т в 2 раза приводило к неустойчивости). В свою очередь расчет по неявной схеме с т = Н приводил к данным, близким к устойчивым явным вычислениям, т.е. наблюдалась экономия вычислений в 20 раз.
Еще более жесткое ограничение для устойчивости присутствует в явной схеме при V = 10"2: прежний шаг Н = 0.25 • 10"3 требовал значения т = Н/200 (увеличение т в 2 раза порождало неустойчивость). Но при увеличении коэффициента вязкости решение становится более гладким, поэтому необходимость в указанном ранее значении Н отсутствовала. Без потери точности можно использовать Н = 0.5 • 10"3 с ограничением т = Н/100, хотя расчет по неявной схеме с т = Н, как и прежде, приводил к визуально неотличимым от устойчивых результатам.
Суммируя вышесказанное о расчетах, можно утверждать, что учет вязкости при моделировании плазменных колебаний требует обязательного применения неявных разностных схем. Конечно, при решении нелинейных задач это обычно приводит к дополнительным внутренним итерациям, но
использование явных схем для подобных постановок может повлечь необоснованную утилизацию вычислительных ресурсов.
Заключение. В работе численно анализируется влияние вязкости на нерелятивистские колебания холодной плазмы, возбуждаемые коротким мощным лазерным импульсом. С этой целью построена неявная разностная схема типа Мак-Кормака, имеющая более слабое ограничение на устойчивость, чем явная схема. При этом неявная схема имеет безытерационный характер реализации, что увеличивает ее вычислительную эффективность в десятки раз. Показано, что учет вязкости плазмы может быть причиной не только затухания амплитуды плазменных колебаний, но и изменения формы решения. При увеличении коэффициента вязкости у решения наблюдается седловая точка, которая сохраняется во времени, что согласуется с результатами, полученными другими авторами.
Результаты работы могут быть обобщены на постановки задач большей размерности и учет дополнительных физических факторов в модели плазмы (электрон-ионные соударения, тепловое давление и пр.).
СПИСОК ЛИТЕРАТУРЫ
1. Александров А.Ф., Богданкевич Л.С., Рухадзе А.А. Основы электродинамики плазмы. М.: Высшая школа, 1978.
2. Гинзбург В.Л., Рухадзе А.А. Волны в магнитоактивной плазме. М.: Наука, 1975.
3. Силин В.П. Введение в кинетическую теорию газов. М.: Наука, 1971.
4. Силин В.П., Рухадзе А.А. Электромагнитные свойства плазмы и плазмоподобных сред. 2-е изд., испр. М.: Книжный дом ЛИБРОКОМ, 2012.
5. Esarey E., Schroeder C.B, Leemans W.P. Physics of laser-driven plasma-based electron accelerators // Rev. Mod. Phys. 2009. 81. 1229-1285.
6. Bulanov S.V., Esirkepov T.Zh., Hayashi Y. et al. On some theoretical problems of laser wake-field accelerators // J. Plasma Phys. 2016. 82, N 3. 905820308(1-55).
7. Горбунов Л.М. Зачем нужны сверхмощные лазерные импульсы? // Природа. 2007. 21, № 4. 11-20.
8. Чижонков Е.В. Математические аспекты моделирования колебаний и кильватерных волн в плазме. М.: Физматлит, 2018.
9. Skorupski A.A., Infeld E. Nonlinear electron oscillations in a viscous and resistive plasma // Phys. Rev. Lett. E. 2010. 81. 056406(1-11).
10. Chizhonkov E.V., Frolov A.A. Effect of electron temperature on formation of travelling waves in plasma: Kinetic and hydrodynamic models // Russ. J. Numer. Anal. Math. Modelling. 2023. 38, N 2. 63-74.
11. Rozanova O., Chizhonkov E., Delova M. Exact thresholds in the dynamics of cold plasma with electron-ion collisions // AIP Conf. Proc. 2020. 2302, N 1. 060012(1-12).
12. Chizhonkov E.V., Delova M.I., Rozanova O.S. High precision methods for solving a system of cold plasma equations taking into account electron-ion collisions // Russ. J. Numer. Anal. Math. Modelling. 2021. 36, N 3. 139-155.
13. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд., испр. и доп. М.: Физматлит, 2012.
14. MacCormack R. W. A numerical method for solving the equations of compressible viscous flow // AIAA J. 1982. 20, N 9. 1275-1281.
15. Davidson R.C. Methods in Nonlinear Plasma Theory. N.Y.: Academic Press, 1972.
16. Розанова О.С., Чижонков Е.В. О существовании глобального решения одной гиперболической задачи // Докл. РАН. Математика, информатика, процессы управления. 2020. 492, № 1. 97-100.
17. Rozanova O.S., Chizhonkov E.V. On the conditions for the breaking of oscillations in a cold plasma // Z. angew. Math. und Phys. 2021. 72, N 13. 1-15.
18. Розанова О. С., Чижонков Е.В. Об аналитическом и численном решении одномерных уравнений холодной плазмы // Журн. вычисл. матем. и матем. физ. 2021. 61, № 9. 1508-1527.
19. Dafermos С.M. Hyperbolic Conservation Laws in Continuum Physics. The 4th ed. Berlin; Heidelberg: Springer, 2016.
20. Sheppard C.J.R. Cylindrical lenses — focusing and imaging: a review // Appl. Opt. 2013. 52, N 4. 538-545.
21. Morton K.W., Sobeye I.J. Discretization of a convection-diffusion equation // IMA J. Numer. Anal. 1993. 13. 141-160.
22. MacCormack R. W. The effect of viscosity in hypervelocity impact cratering //J. Spacecr. Rockets. 2003. 40, N 5. 757-763.
23. Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985.
24. Furst J., Furmanek P. An implicit MacCormack scheme for unsteady flow calculations // Computers & Fluids. 2011. 46. 231-236.
25. Чижонков Е.В. О схемах второго порядка точности для моделирования плазменных колебаний // Вы-числ. методы и программир. 2020. 21. 115-128.
26. Фролов А.А., Чижонков Е.В. О применении закона сохранения энергии в модели холодной плазмы // Журн. вычисл. матем. и матем. физ. 2020. 60, № 3. 503-519.
27. Verma P.S., Soni J.K., Sengupta S., Kaw P.K. Nonlinear oscillations in a cold dissipative plasma // Phys. Plasma. 2010. 17. 044503(1-4).
Поступила в редакцию 08.09.2023
