Сеточно-характеристический метод для задач лазерного термоядерного синтеза

1. Наиболее полные математические модели рассматриваемых задач основаны на использовании уравнений магнитной гидродинамики с учетом реальных свойств плазмы (группа ''уравнений состояния") и большого числа различных физических процессов: поглощение внешнего лазерного излучения, электронная теплопроводность, электрон-ионная столкновительная релаксация, собственное излучение плазмы, перенос энергии быстрыми электронами и ионами, спонтанные магнитные поля и рентгеновское излучение, термоядерное горение и т. д. Уточнение математического описания группы ''уравнений состояния", а также различных процессов в плазме составляет предмет интенсивно ведущихся физических исследований.

В настоящей главе используется упрощенная математическая модель, воспроизводящая основные черты рассматриваемых задач, и основное внимание уделяется исследованию многомерных эффектов, сопровождающих сжатие и нагрев микромишеней. Из разнообразных физических процессов учитываются первые три из перечисленных выше (т.е. поглощение внешнего лазерного излучения, электронная теплопроводность и электрон-ионная столкновительная релаксация). Электронная и ионная компоненты плазмы рассматриваются как идеальный газ.

В двумерном нестационарном случае с учетом сделанных предположений система уравнений газовой динамики в фиксированной в пространстве ортогональной системе координат s = х_х, г = х₂, = х₃ может быть записана в виде.

Здесь к = const — показатель адиабаты; c_et c_f = const — удельные теплоемкости электронной и ионной компонент при постоянном объеме; геометрические параметры, А = НН₂Н_Ъ> Si ₃ = Н_уН_Ъу S₂3 ⁼ Н₂Н₃ зависят от выбора системы координат и характеризуют элементарный объем Adsdrd^ и его поверхность; Н_у, Н_1у Н_ъ- коэффициенты Лямэ. Для используемой далее сферической системы координат с радиальной координатой г, азимутальным и меридиональным углами s, у (см. (1.12) из гл. V) имеем Hi = г, Н₂ = 1, #₃ = г sin s.

Для величин Q_s, Q_Ty Q_ei принято следующее описание (см.: [245]): 238.

Здесь р_с — критическая плотность. Плотности потока энергии падающего р⁺ и отраженного от поверхности р (/, s, г) = р_с лазерного излучения q~ определяются из одномерных уравнений переноса (вдоль направления г).

с граничными условиями.

В частном случае полного поглощения на слое с критической плотностью имеем

где r_r(r, s) — текущий радиус внешней границы мишени, q^t, s) — плотность потока падающего лазерного излучения на границе мишени.

Решение уравнений (1) отыскивается в области t > 0_У 0 <$?, 0.<

^^гг (^^s) ^ПР^И начальных условиях:

На внешней границе мишени (граница с вакуумом) требуется выполнение условий.

В центре мишени (г = 0) и на поверхностях s = 0, s* в качестве граничных условий используются соответствующие условия симметрии. В случае, если мишень является оболочкой с вакуумной полостью, внутренняя граница мишени до момента ''схлопывания" является границей с вакуумом и условия на ней аналогичны (8).

В некоторых из рассматриваемых ниже задач (''конические" мишени, разд. 3) система уравнений (1) дополнялась уравнением для массовой концентрации? одной из компонент двухкомпонентной смеси идеальных газов

при этом теплофизические свойства плазмы (удельные теплоемкости при постоянном объеме для электронной и ионной компонент с_е, С/, коэффициент поглощения К, коэффициент электронной теплопроводности к_е и т. д.) определялись обычным для многокомпонентной смеси образом, т. е. как функции плотности, температуры и массовой концентрации.

В некоторых расчетах вместо одного из уравнений энергии (для электронной компоненты) использовалась полная энергия Е = е_е + e_t + (и² + v²)/2.

и соответствующая дивергентная форма этого уравнения.

Соотношения (1)—(10) выписаны для фиксированной в пространстве системы координат, т. е. принят несколько нетрадиционный для данного класса задач эйлеров подход. Выбор переменных Эйлера предпочтительнее при значительных сдвиговых деформациях материала, возможных в многомерных задачах. Он позволяет в ряде случаев более детально описать поглощение внешнего лазерного излучения в малоплотной ''короне" разлетающейся части материала мишени, где в традиционном лагранжевом подходе довольно сложно обеспечить необходимый пространственный шаг разностной сетки (так же как и вблизи центра мишени) и др. Вместе с тем в переменных Эйлера в ряде случаев трудно обеспечить необходимую подробность разностной сетки в области с большой плотностью материала, в особенности в задачах для очень тонких оболочек. При использовании переменных Лагранжа в подобных задачах, когда в областях с высокой плотностью автоматически группируется большое число узлов разностной сетки, решение в этих областях воспроизводится более детально.

В связи с этим при численном решении некоторых одномерных задач, где преимущества лагранжева подхода более очевидны, и в основном с методическими целями проводились также расчеты с использованием переменных Лагранжа и консервативного варианта сеточно-характеристического метода (4.2.11). В этом случае при использовании массовой лагранжевой координаты т так, что dm = pr^vdr, имеем.

Qet, уравнения состояния, начальные и граничные условия определяются соотношениями (2), (3), (7), (8). Коэффициент v =0,1,2 соответственно для плоской, цилиндрической и сферической геометрии.

2. Системы уравнений (1), (11) в отсутствие члена, связанного с электронной теплопроводностью в правой части уравнения энергии для электронной компоненты (Gr), имеют гиперболический тип, и для них могут быть использованы те или иные разностные схемы из числа рассмотренных в гл. Ill, IV (например, явные консервативные схемы (4.2.11), (4.4.6)) с добавлением соответствующей аппроксимации Q_T. Присутствие указанного члена с электронной теплопроводностью, вообще говоря, делающего эти системы уравнений параболическими, вызывает необходимость их неявной аппроксимации, поскольку явная аппроксимация, в отличие от условий типа (4.2.12), (4.4.5), приводит к более жесткому ограничению на шаг интегрирования т типа г < h²1(2 шах к_е), характерному для уравнения теплопроводности, и делает расчеты затруднительными. По аналогичной причине неявной аппроксимации требует столкновительный член в правой части уравнений (1), (11), т. е. Q_e{, явная аппроксимация которого в некоторых случаях (например, на переднем фронте тепловой волны, распространяющейся по холодному газу) приводит к не менее жестким ограничениям на г.

В данной главе для решения описанных выше краевых задач используется именно такой подход к построению общего вычислительного алгоритма, связанный с расщеплением исходных уравнений по ''физическим процессам". При этом для аппроксимации гиперболической части систем уравнений (1), (11) использовался консервативный вариант сеточно-характеристического метода, т. е. явные разностные схемы, аналогичные (4.4.6), (4.2.11) соответственно для двумерных и одномерных нестационарных задач, а члены, связанные с электронной теплопроводностью и электронионной столкновительной релаксацией, аппроксимировались по неявной схеме.

При реализации на ЭВМ удобнее вначале аппроксимировать левую часть систем (1) по схеме (4.4.6) (или (11) по аналогичной схеме (4.2.11)) в одномерном нестационарном случае при использовании переменных Лагранжа), что позволяет вычислить окончательные значения плотности Р"У = Рт/ • компонент вектора скорости м"⁺/ = и_т1, и" ⁺/ = v_ml и пред;

варительные значения электронной Т_ет1 и ионной T_iml температур. В частности, в случае двух пространственных переменных (система (1)), выбирая в качестве искомого вектора u = | р, рм, pv, ре_е, ре Л, имеем.

В этих соотношениях о_;— = г /Л, — (/ =1,2), все функции от р, и, и, Т_е, T_t в узловых точках разностной сетки (tⁿ_f s_m, r_/), (/", s_m±j, г,),.

(/", ^r/±i) вычисляются в соответствии с (2), (3), (6), а в полуцелых точках с индексами m ± ½, / и m, I ±½ усредняются по их значениям в ближайших узловых точках. Геометрические параметры S_i3, S₂₃, Д, Д¹, как отмечалось, характеризуют поверхность и объем элементарных ячеек с центрами в узлах разностной сетки (S₁₃ = H_lH₃, S₂₃ = Я₂Я₃, Д = Д¹ = Я₂Я₂Я₃). Чтобы учесть «координатную» особенность (точнее, неопределенность), возникающую при расчетах в сферической системе координат узловых точек на оси симметрии (s = 0, s = л) ив центре мишени (г = 0), выражения для S₁₃, S₂₃, Д, Д¹ для этих точек соответствующим образом видоизменялись.

покрывающей фиксированную в пространстве область интегрирования -hi + h_x, 0 = const. Искомые параметры в граничных узловых точках /*/=(/- 1) Л₂, / = 1,. .. , L, s_x = = (Л/ - 2) Л, =

= s. + /ij, а также в центре мишени вычисляются с использованием условий симметрии течения относительно ''линий" s = 0, s = s_# и точки г = 0. Параметры в остальных граничных узловых точках гl = (L — 1) Л₂ -R, 5_m = (m — 2) Л, m = 2, … — 1 определяются линейной экстраполяцией.

(либо просто сносом) по их значениям в ближайших внутренних узловых точках. До того момента, как возмущения, возникающие на внешней границе мишени г = г_т (0, s) (расположенной внутри области интегрирования ^гг (0> ^s) < R)> достигнут внешней границы области интегрирования г — R, параметры в этих точках остаются невозмущенными, т. е. теми же, что были заданы в начальных условиях (в которых в области вакуума задавался некоторый, слабо влияющий на решение фон р = р*^р₀, Т_е-Т_{= Г* < 0i и = v = 0). Затем на этой границе довольно быстро устанавливается сверхзвуковое течение (разлет за область интегрирования), так что влияние покидающей область интегрирования массы с малой плотностью (и указанной процедуры экстраполяции) на остальное течение несущественно. Это влияние заключается в том, что при t < t_k в области г > R (с малой плотностью) происходит частичное поглощение падающего излучения, которое в данном подходе не учитывается.

На следующем этапе производится учет электронной теплопроводности и обмена энергией между электронной и ионной компонентами (вычисление окончательных значений ТГ^¹, по следующей неявной схеме:

получим типичную для нелинейного уравнения теплопроводности краевую разностную задачу относительно ГЦ**, которая решалась итерациями с использованием одномерных прогонок на лучах т = const. В некоторых расчетах на заключительном этапе вместо (14) использовались другие неявные экономичные схемы, с использованием прогонки как по лучам т = const, так и в поперечном направлении. Схема решения системы одномерных уравнений (11) при использовании переменных Лагранжа аналогична рассмотренной выше и отличается лишь выбором вектора и = /р, и, Е_у €/(и несколько иным видом матрицы с = (Л/2) Я" ¹!А112. В случае использования (10) вместо уравнения энергии для электронной компоненты в (1) также несколько видоизменяется вид матрице* в (13). Там, где использовалось уравнение (9) для массовой концентрации ?, его аппроксимация производилась по явной схеме, как и вся гиперболическая часть системы (1), с учетом характеристических направлений этого уравнения, т. е. с учетом знаков и, v.