CC BY
51
УДК 516.6.531.33
В. В. Кузенов
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОСНОВНЫХ ПЛАЗМОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК В ЛАЗЕРНОМ ФАКЕЛЕ ВБЛИЗИ АЛЮМИНИЕВОЙ МИШЕНИ
Приведена постановка задачи и выполнено численное моделирование плазмогазодинамических процессов в приповерхностном лазерном факеле. Численная реализация плазмодинамической модели основана на многоблочной многосеточной технологии расчетов на неортогональных структурированных сетках с использованием схем расщепления по физическим процессам и направлениям. Решение расщепленных уравнений Рейнольдса найдено с помощью разработанного варианта нелинейной квазимонотонной компактной дифференциально-разностной схемы повышенного порядка точности, которая в пространственно гладкой части численного решения позволяет получить шестой порядок точности. E-mail: [email protected]
Ключевые слова: лазерное излучение, конденсированные среды, плазмо-динамическая модель, численные методы, структурированные сетки.
Вопросам взаимодействия лазерного излучения с конденсированными средами посвящено значительное число работ как теоретических, так и экспериментальных [1-8].
Плазменное образование, возникающее над поверхностью конденсированной преграды при действии на нее мощного лазерного излучения, может использоваться (в случае его осаждения на подложку) для получения и обработки многослойных металлических зеркал, создания различного вида пленок в полупроводниковых технологиях. Данный механизм осаждения может быть использован для образования слоя нанотрубок на подложке и т.д.
Определенный интерес проявляется в последнее время к применению мощных лазеров для прецизионной обработки металлических поверхностей ультракороткими лазерными импульсами. В случае использования импульсных лазеров пико- и фемтосекундного диапазонов качество обработки поверхностей значительно выше, чем при использовании относительно длинных (наносекундных) лазеров.
Одним из путей создания компактной ракетной техники является использование лазерных ракетных двигателей. В настоящее время известны два типа таких устройств:
первый — начальный импульс лазерного излучения фокусируется на поверхности металлической преграды и испаряет ее, а повторный импульс осуществляет оптический пробой, нагрев и ускорение [3];
второй тип основан на механизме, который может быть использован в лазерных воздушно-реактивных двигателях, — в этом случае пе-
редача импульса стенкам летательного аппарата происходит при действии на них ударной волны, возникающей (благодаря оптическому пробою газа сфокусированным лазерным излучением) в газовой среде, заполняющей рабочую камеру [4].
Отметим также, что лазерное излучение и связанное с ним объемное энерговыделение в определенной области сверхзвукового потока (например, в воздухозаборнике ГПВРД) позволяет изменять в нужную сторону структуру течения (систему взаимодействующих ударных волн) потока.
Метод построения адаптивных сеток. При решении задачи лазерного воздействия на металлические преграды необходимы простые и эффективные способы дискретизации расчетных областей. Такого рода сетки можно получить решая эллиптические уравнения, которые позволяют найти отображения расчетной криволинейной области в параметрический квадрат (в двумерном случае) и в параллелепипед (в трехмерном случае). Введем в декартовой системе координат ХУ2 прямоугольный параллелепипед АВГЕБСОН, который непрерывно дифференцируемым образом необходимо отобразить в криволинейный параллелепипед (гексаэдр) А' В' Г' Е' Б' С' О'Н'. При этом прямоугольная сетка, нанесенная на область АВГЕБСОН, образует гладкую криволинейную сетку в области А' В' Г'Е'Б' С' О'Н'.
Обозначим через г радиус-вектор в системе координат ХУ2 и введем вектор и = г * — г, характеризующий смещение точек. Здесь г иг * — радиус-векторы точек областей до (г € АВСБСОН) и после (г * € А'В'Г'Е'Б'С'О'Н') преобразования. Тогда уравнения, определяющие величины смещений их ,иу ,иг и применяемые для построения регулярных адаптивных сеток, близких к ортогональным, в трехмерной декартовой системе координат ХУ2 имеют вид [9]:
Ax-TT- ( F
ßUx
+
+
dx \ dx
(i - д 2 dz
(1 - *) д
2 dz
+ kopt
+
(1 - *) д
dUx dz dUx dz d_ dx
F
F
+
+
dy '
(1 + *) 2
(1 + 2
dUx dy d_ dx д dx
+
922933 dUx gii dx
— F
FT*
y
öU„
y
d_ x d_
x
d_ dy
F
F
dUz
z
dUz
z
d
z
9339h dUx 922 dy
+
9ii922 dUx
933 dz
+
+
= 0;
A d
Xy dy
'Fd3L
ч dy
+
(1 - *) д
2 dx
pdUy
x
+
+
(i - д 2 dz
FOUy
z
+ M д x
+
(1+ q)\ д
2 {Oy
922933 dUy
911
dx
+ д. Ш
y z
д_ dy
д_
z
+
= 0;
К-
z dz
F
ßU7,
+
z
(1 - *) д
+
(1 - *) д
x
F
Uz
x
+
F
Uz
y y
+
(1 + *) 2
(F0Ux
z x
+
z y
+
+ kopt ^ dx
д ( /922933 dUz
911
x
d_ dy
d_
z
9339h 9UZ 922 dy 9ii922 dU:
+
= 0.
9зз дг
Здесь принято, что все коэффициенты, входящие в данную систему уравнений, равны между собой Е = (V/, V/)а/2 + е, г € (ж, у, г}. Граничные условия, необходимые для решения данной системы уравнений, задаются следующим образом:
иъ\г = Гг \д(Л'Б'Е'Е'В'С'0'Н') — Гг\д(ЛБЕЕВССН) , г € (х,у,г} •
Компоненты ковариантного метрического тензора, входящие в систему уравнений, определяются с помощью соотношений
9
kj = si =
1 i = j 0,i = j
12 3
q = x,q = y,q = z.
« mm (к.) выбирается как величина,
hi&A'B'F'E' D'C'G'H'
9ък ^ д<? дqk
а=1
При этом / является функцией, управляющей адаптацией расчетных узлов, коэффициент а задает необходимую степень сгущения узлов расчетной сетки, коэффициенты Х^ > 0 и а € [-1,1]. Параметр е
имеющая ненулевое значение, и порядка шага сетки к, чтобы избежать особенности в тех узлах, где V/ = 0. Для численного решения уравнений, применяемых для построения регулярных адаптивных сеток, представленных в операторной форме Аи = 0, используется модифицированный вариант итерационного метода вариационного типа — метод минимальных невязок [10]. Здесь отметим также, что, при применении многоблочной технологии расчетов, рассматриваемая численная методика позволяет построить квазиортогональные структурированные сетки в областях сложной формы.
Математическая модель лазерного факела. Решение уравнений Рейнольдса конечно-разностным методом возможно, когда границы расчетной области совпадают с координатными линиями в некой обобщенной системе координат £пС • Эта система координат получается путем перехода от цилиндрической г,г,р к произвольной криволинейной системе координат £пС• Отметим также, что на равномерной сетке в обобщенной системе координат £пС сохраняется порядок конечно-разностной аппроксимации уравнений при сгущении узлов сетки в физическом пространстве с координатами г, г, р.
Введем преобразование координат вида
г = г(С,п,С), г = г(С,п,С), р = р(С,п,С).
Отметим, что при известных в физическом пространстве координатах узлов сетки в расчетной области £, п, С метрические коэффициенты в общем случае могут быть найдены путем численного дифференцирования [11].
При численном поиске решения систему уравнений, описывающих процессы в плазме лазерного факела, следует привести к безразмерному виду. Тогда исходные уравнения в безразмерных переменных позволят описывать целые классы течений. Для получения безразмерного вида отнесем все газодинамические переменные, входящие в систему уравнений, к их характерным значениям, а пространственные (£, п) и временную (£) переменные — к характерному размеру Ь* и характерному времени £ *. Введем следующие обозначения безразмерных переменных:
U* t L* i = _L L n = L L V« = V« V*' Vn = in " V*'
u = u V*' v = v V*' T = T T ' J- * e = L e
и PI L ßs = _ ßs ß *' As = As A Y = Cp * Cv * q = L ш q* ш*
Р
Р = —, Р
Причем, используя дополнительные алгебраические связи между
характерными величинами £ * = —*, е* = С,и*Т *, е* = V*2, Р* = р*У*2,
V*
Р
д* = Р^* получаем для числа Эйлера Ей = ——— и температурного
р *к2
V2
фактора в = — упрощенную форму записи: Ей = 1, в = 1.
е
Плазмодинамические процессы, протекающие в лазерной плазме, могут быть определены с помощью системы уравнений вязкой одно-температурной радиационной плазмодинамики. В безразмерных переменных эта система уравнений примет следующий вид:
dp + iö(jpVi) + iö(jpVv) = _; dt J dj J dn r
dpu | 1 d (JpuVg) | 1 d (JpuVv) _ dP dP a— + —■ dt J dj J dn r dj П dn a r Re'
dpv | i д (JpvVj:) | i д (JpvVn) _ dp dp apuv + S±. dt J dj J dn z dj dn a r Re'
dpe + 1 d(JpeVi + J ^ qt4) + 1 d(JpeVv + J E ) = dt + J dj + J dn =
_ P fd (JVg) d (JVn) 1 Pu peu S^ J \ dj dn j r r Re'
Y Re t
Se = ^D + -f div (As gradT) +-*ql.
Pr p*e*
Здесь Sr ,Sz математически описывают силы, возникающие в лазерном факеле за счет наличия в нем сил вязкого трения; Se — объемное энерговыделение, появляющееся из-за работы сил трения psD (где D
— диссипативная функция), переноса теплоты процессами теплопроводности div (As grad T) и энерговыделения QL, обусловленного действием лазерного излучения на плазму окружающей среды и паров материала преграды [12]; Re = *p* * — число Рейнольдса; Pr = р*
А*
— число Прандтля.
Величины Ql, Sr,Sz, div(As gradT), D определяются с помощью следующих выражений:
ql = Х^>(z,r = 0)рь exp( - — )exp
rl
- J Х^(x,r = 0)dx 0
nRL
_ 1 d (J {jr Orr + jz Orz }) , 1 d (J {nr Orr + nz Orz }) Sr — — — + ~ г
J
dj
J
n
+a-
u u u
2ps jr dj n dn — 2aps — r
1 d ( J {jr Ozr + jz Ozz}) , 1 d ( J {nr Ozr + nz Ozz}) , Orz
Sz = ~r-^--h i-я--h a—;
J dj J dn r
div (As grad T) = 1 d{AsJ (j2 + Т + As J (jr nr + ^) Tn} +
J
j
z
n
1
r
1 д {AsJ (nir + nziz) T« + As J (n2 + n2) Tn}
J
dn
As Г dT dT\
+aT\irdi + nr ;
D 2 |^err + ezz + J +
- V divV' 2
3
divV = -J
д (JV«) д (JVn)
di
dn
u
+ a-; r
du du dv dv u
err Sr^T + ' ezz + ' еи» a j
di dn di dn r
du du dv dv
erz = di + nz d^ + ^ di + nr dn
Orr = ßs
Ozz = ßs
4 ( du du\ 2 / dv d^ 2u
3 ^rdi + nrdn) - 3 [b di + nzdn) - a37.
4 dv dv 2 du du 2u
3 (,izdi + nzdn) - 3 (,irdi + nr dn) - a3r
Orz = ßs
du du
iz di + nz dn
+
dv dv di +"r dj
где u(r, Z' t)'V(r-' г, t) — проекции вектора скорости V(r, z, t) на оси R и Z; e — удельная внутренняя энергии плазмы; J = d (r z)/d (i' n)
— якобиан перехода от цилиндрической системы координат rz к криволинейной системе координат in; V« = iru + izv, Vn = nru + nzv
— контравариантные компоненты вектора скорости V в криволинейной системе координат i' n; P'P — плотность и давление плазмы;
Y2qi« 'Y2qin — проекции вектора плотности потока лучистой энергии
i i
Vна оси криволинейной системы координат i и n; a = 0 соответствует плоскому, a = 1 — осесимметричному случаям течения.
Для определения плотности плазмы аблирующих паров материала
металлической преграды в систему приведенных уравнений вводится
dPg V
дополнительное уравнение неразрывности + VVpg = 0, позволяющее определить пространственное положение контактной границы, которая отделяет границу аблирующих паров материала металлической преграды от газа окружающей среды.
Генерации магнитного поля, возникающего при взаимодействии лазерного излучения с различными мишенями, посвящено много работ (см., в частности, [13-15]). В зависимости от условий облучения магнитное поле может быть обусловлено различными механизмами: магнитотепловой неустойчивостью, резонансным поглощением
e
при наклонном облучении, термоЭДС в плазме и т.д. При действии на мишень излучения лазера с плотностью потока 108... 109 Вт/см2 первые два указанных механизма отсутствуют и спонтанное магнитное поле будет обусловлено термоЭДС. В этом случае основным является градиентный механизм генерации магнитного поля, когда вихревые токи в плазме образуются в результате неколлинеарности градиента электронной концентрации Vne и градиента температуры VT.
Уравнение генерации магнитного поля B c учетом источника магнитного поля в лазерной плазме и принципа суперпозиции (в линейном приближении) имеет вид
dB
-Ж = rot
V х B
c2 t* / rotB \ t* ck rv7 ^
" T~ 72 rOM - - ^Ti- [УПе X VT] ,
4n L2 \ о BL enP
<-e
где пе — концентрация электронов; Т — температура плазмы лазерного факела.
Уравнение переноса излучения представлено в виде системы уравнений диффузионного многогруппового приближения [16]:
+ + -U =4х,о,Т«
c dU; c dUi
+ xm = 0 + XiQi-n = o,
(1)
3 ' 3 дп
где иг (у, г,Ь) — плотность лучистой энергии в г-й спектральной группе; Хг — спектральный коэффициент поглощения.
Система уравнений, описывающих процессы нагрева и испарения материала поверхности металлической преграды под действием лазерного и теплового излучения из объема плазмы с плотностью (¡г без учета гидродинамических процессов в конденсированной среде, состоит из квазиодномерного уравнения теплопроводности в подвижной (связанной с фронтом волны испарения) системе координат с осью 02, перпендикулярной поверхности, и осью 0У, параллельной поверхности,
дтя д2т3 ТЛ дтя
= ам^^т + Уо-
дЬ дг2 дг
с граничными и начальными условиями
дТ
км(0, г, Ь) = (0, г, Ь) - р (0, г, Ь) V (0, г, Ь),
Т3 (г ^ ж, г, Ь) = То, Та(0, г, Ь = 0) = То
и системы уравнений, определяющих кинетику испарения поверхности конденсированного вещества в рамках модели с кнудсеновским слоем [17]
T (0,t) Ts (0,t)
1 + n
(y — 1) m (Y +1)2
— v n
(Y — 1) m (Y +1)"2
P (0,t) = Ts (0,t) Ps (0, t) V T (0, t)
m2 + 1 ) em2 erfc (m) —m 2 / \/n
m =
Ps (t) = pi exp V (0,t)
+ 1 Ts (0,t) 2 T (0,t)
1
RTi
+
1 — /nmem2 erfc (m) Ti - "
Ts (0, t)
(2)
'S
; p (0, y, t) v (0, y, t) = pmv),
У^ЩмУЛ
где Т (г, г, ¿) — температура конденсированной среды в точке г,г в момент времени ам,км,рм — соответственно коэффициенты температуропроводности, теплопроводности и плотность материала; У — скорость волны испарения; р3, Р.в — давление и плотность насыщенного пара конденсированного вещества при температуре поверхности Тв (г = 0, г, ¿); Я — универсальная газовая постоянная; Т1 — значение температуры поверхности конденсированной среды, соответствующее давлению насыщенного пара рь — скрытая теплота испарения; 3 — молярная масса пара; Т (0, ¿) ,р (0, ¿) ,У (0, ¿) — температура, плотность и скорость плазмы на внешней границе кнудсеновского слоя в точке (г = 0, г) на момент времени ¿; 7 — показатель адиабаты паров конденсированного вещества.
Коэффициент поглощения лазерного излучения хш, см-1, для СО2-лазера (Хьаг = 10,6 мкм) определяется с использованием механизма континуального поглощения, обратного механизму тормозного излучения электронов в условиях ЛТР [18]:
Хш = 2,82 • 10-29пе(п+ + 4п++)Т-3/2 ^ (2,17 • 103Тп-1/3) ,
где пе, п+, п++ — числовые концентрации электронов, одно- и двукратно заряженных ионов; Т — температура, К.
Коэффициент поглощения излучения неодимового ^-лазера (Хьаг = 1,06 мкм) в воздушной плазме находится с использованием формулы [16]
=
= ( а.)
nY J
1+
10,2
Y
2
_ 1,1 • 105PnP-1T-1/6
^Y) г 'm 4 п
1,55 + 2,9^3^ - 1
1,38 • 105
x exp '
. 1,36 • 104
1 — exp
T
где а — интегральный коэффициент поглощения в атомных линиях; Y — ширина атомных линий, определяемая штарковским уширением; T — температура, K; PN — парциальное давление атомов азота, атм; Pe — парциальное давление электронов, атм.
Расчет входящих в данную систему уравнений термодинамических (e (T, р), P (T, р)) и оптических (х (T, р)) параметров рабочих сред проводился в рамках приближения локального термодинамического равновесия с использованием компьютерной системы ASTEROID [19].
При расчете оптических характеристик разбиение спектра на группы осуществлялось с учетом особенностей электронной структуры атомов, входящих в состав излучающей плазмы и газа окружающей среды. Весь спектр был разделен на 7 групп с границами интервалов (эВ): 0,1; 3,14; 5,98; 6,52; 7,95; 9,96; 18,6; 200 (для воздуха).
Турбулентные коэффициенты вязкости (рт) и теплопроводности (As) рассчитываются с привлечением гипотезы Буссинеска, в соответствии с которой эффективная вязкость рт газового потока определяется по формуле рт = рт + Pt, где рт — динамический коэффициент вязкости, учитывающий атомно-молекулярные столкновительные процессы; pt — коэффициент турбулентной вязкости, для определения которого используется q—ш модель Кокли [20].
Численный метод решения. Численное решение разработанной нестационарной двумерной радиационно-газодинамической модели базируется на методе расщепления по физическим процессам и пространственным направлениям.
При этом решение двумерной нестационарной системы уравнений вязкой однотемпературной радиационной плазмодинамики, записанных в векторной полудивергентной форме в криволинейной системе координат
д£ dp dG dt дп Re
dFV + dGv + s
дп v
строится с использованием метода расщепления по физическим процессам и направлениям. Векторы, входящие в данную систему уравнений, имеют вид
p pVt pVn
pu puVt puVn
pv ; F — J pvVt ; G — J pvVn
pe peVt peVr,
pq pqVt pqVn
pш [шУ^ pdVr
а-
pu r
dP
dP ~ör
2
pu
Cr^rr + Vr—+ а-
дп P
r puv
r
дп
P f д(Щ) + d(JVn) |
J l dC
pqu pq
а
а
r
puu r
2ш*ш
drj
CfD
- ^IC
Ш*
Pu peu + а--+ а-
p e
-Ql
V.
L*
2
2 V
C,D
Vi
L*
- -шi-1"-ш divV - ш^ш2 3 L"
2 2 V - -ш*—"-ш divV ) - С2ш2ш2
з и 1 1
t
Fv — J
G v — J
0
Cr Orr + Cz Orz Cr Ozr + iz zz
^ (Ae (Cr2 + Cz2) T + Ae (Crnr + CzVz) Tn} PEq (Cr2 + C2) qt + ^E q (CrVr + CzVz) qV
PEw (C2 + C2) Ш + (CrVr + Cznz) шп 0
nr Orr + Vz Orz Vr zr + nz Ozz
^ (Ae (nr Cr + Vz Cz) T + Ae (v2 + n2) Tn}
PEq (ПгCr + VzCz) qt + PEq № + П2) %
PEw (Vr Cr + Vz Cz) <Ш + PEw (П2 + П22) Шп
Sy - J
2ps
a-
du du £r д£ + Vr
u
— 2aps — r
a
a
a
a
Orz
r
Pr r
PSq
r r
дТ дТ
£r ~е£ + Vr д дq + дя.
£r д£ + Vr дп . дш дш
д£ дп
+ Ps D
На первом временном дробном шаге г Е [г, г + Дг/3] с использованием явного интегроинтерполяционного метода и соответствующих краевых условий решается гиперболическая (невязкая) часть системы уравнений.
Для первого дробного шага г € [£,£ + Дг/3] использована нелинейная квазимонотонная компактная разностная схема повышенного порядка точности, которая в пространственно гладкой части численного решения позволяет достигнуть шестого порядка точности [21]:
ди„
i,j + Fi+1/2,j — Ft-1/2>3 +
G
ij+1/2
/2 — G i,j-1/2
+ Sij - 0.
дг Д£ Дп
При пространственной аппроксимации гиперболической (невязкой) части системы уравнений принято, что газодинамические параметры иП+1,иПз относятся к центрам расчетных ячеек, в то время как потоки Р1±-1/2 з ,^Пз±1/2 определяются на поверхностях, ограничивающих эти ячейки. При этом для повышения порядка аппроксимации разностной схемы следует восстановить газодинамические
т/й , Ъ л/'К, Ъ
параметры Уг±1/2 з ,Уг з±1/2 справа и слева от границ расчетных ячеек. Дискретные потоки ,3,Сг,з±1/2, необходимые при расчете по дифференциально-разностной схеме, находятся с помощью ве-
т г К, Ъ т г К, Ъ ^
личин ^¿±1/2з,Угз±1/2 на основе решения автомодельной задачи Римана о распаде произвольного разрыва газодинамических величин {р,У^,Уп,Р,е}. В результате чего получается хорошо физически и численно обусловленная вычислительная схема [22]. Здесь У+у2, и У+1/2. - левые и правые значения вектора переменных
У = {р,п,у,У^,Уп,е,Р} на границе ячеек с номерами г и г + 1, а У^з±1/2 (х = Ъ,Я) — левые и правые значения вектора переменных У на границе между ]-й и (] + 1)-й ячейками. При этом любая реконструируемая функция У (х) (х = {£,п}) представляется кусочно-линейными распределениями следующего вида: У (£) = У,3 +
r
+ R (Ind) ([6 - 6j], Y(n) = Yj + R (Ind) (дп) [П - Пj];
j У \
здесь R (Ind) — ограничитель значений производных I —- I и
V 6 / i,j
дП, %>3
Входящие в кусочно-линейные распределения У (С), У (п) про-
(дУ\ (дУ\
странственные производные ——- , —- вычисляются следу-
\dUii \д-п /ц
ющим образом. и и
1. Сначала для дискретной функции У определим приближенное значение ^ первой частной производной по пространственной переменной С или п с шестым (четвертым) порядком точности.
Для этого найдем первую производную . по переменной С или ц по обычной аппроксимационной формуле второго порядка точности
. _ Уг+1 — Уг-1
П _ 2А ,
где А — шаг пространственной сетки в направлении С или г/. Тогда приближенное значение первой производной ^ получается из решения системы уравнений с трехдиагональной матрицей. Эта линейная система уравнений представляет собой компактную схему для нахождения первой производной с четвертым порядком точности:
1 (Я+1 + + ¿и) _ .П. 6
Важно отметить, что данная компактная схема имеет меньший коэффициент при ошибке аппроксимации порядка О (А4), чем обычная схема четвертого порядка точности.
Для нахождения приближенного значения ^ первой частной производной по пространственным переменным С или ц с ошибкой аппроксимации на уровне О (А6) вычисления следует проводить по формуле
Fi = 1^ (' + ^ -1
д2 ,
6E + Д ) До
fi,
где До fi = fi+1 - fi-1, a2fi = fi+1 - 2fi + fi_i. Отметим, что приведенная формула является симметричной конечной разностью шестого порядка точности [23].
Монотонизированные приближенные значения F* первой частной производной определяются на основе следующих соотношений: F* = min mod (FR,FriR), где FR = min mod (Fi+1, Fi) и FR =
= min mod(Fi, Fi-1), min mod (a,b) = 1 (sign (a) + sign (b)) x x min (|a|, |b|);
2. Окончательные значения < —— , —— > определяются
М) и \дПЛ
i,3 \ 'S 1,3
из соотношения
dY\ _ I FENO, при (Y+i - Y) (Y - Yi-i) < 0;
dg) г,з 1 F*, при (Y+i - Yi) (Yi - Yi-i) > 0.
При этом для нахождения F^NO использована поправка из работы [24], имеющая следующий вид: F^NO (A, B) = min mod (A, B). В качестве аргументов A, B поправки F^NO (A, B) использованы соотношения A = (Y+i - Yi) - 2¿i+i/2 и B = (Yi - Yi-1) + 2^i—i/2, взятые
из работы [25], которые делают ее поправкой типа ENO (неосцилли-рующая поправка). Здесь величины ii+i/2, Si-i/2 представляют собой ограниченные значения вторых производных функции Yi.
Поясним, как находятся величины ii+i/2, Si-i/2, используемые для нахождения поправки типа ENO. Для этого на n-м временном слое со вторым порядком аппроксимации рассчитывается вторая производная , Yi+i - 2Yi + Yi-i
функции Yi: si = -—-. Далее по найденным значениям
Д2
si с помощью компактной конечно-разностной схемы определяются ограниченные вторые производные ii+i/2, Si-i/2 с четвертым порядком точности: i2 (¿i+i + 10Si + Si-i) = Si, Si+i/2 = min mod (Si+i, Si),
Si-i/2 = min mod (Si, Si-i).
Отметим, что для решения сложных задач аэротермодинамики широко используются TVD-схемы высокого порядка аппроксимации. Однако они имеют существенный недостаток: понижение порядка аппроксимации до первого в окрестности локальных экстремумов. В настоящей работе в области локальных экстремумов используется поправка FENO, которая избавлена от этого недостатка.
Однако если непосредственно воспользоваться значениями производных < ( д— ) , (д— ) >, то в области больших градиентов га-
{\ди i,3 \дП J ij)
зодинамические величины испытывают мелкомасштабные (3 %) колебания. Поэтому значения производных следует ограничить с помощью специальной поправки R (g, Ind), R (п, Ind).
Для этого проводится расчет ограничителей R (g, Ind) ,R (п, Ind) на основе следующей процедуры, которая повторяется несколько раз.
В начале процедуры исходный массив величин Yi,3 расширяется на величины
и
' Уя,г(0 при У^г + У^г+1 > 0, УеЖ) + УЬ>г+1(0
—,-— при У^ + У^+1 = 0,
, Уь,г+1(0 при У{,г + У^г+1 < 0
Уе,] (п) при Уп,3 + Уп,з+1 > 0, ¥я,3 (п) + Уь,]+1(п) при У + У =0
-2- при Уп,] + УП,]+1 = 0,
„ Уь,3+1(р) при Уп,] + Уп]+1 < 0.
Далее по тексту будем обозначать расширенный массив так же, как и раньше — У,].
В каждой ячейке с номером г для каждой восстанавливаемой величины У1,] осуществляется расчет индекса немонотонности 1ид,(У):
~ 1 —Уг+2,] + 16Уг+1,] — 30Уг,] + 16^-1,] — Y-2j \ 1и4 (У) = у-^2-1-т
\—Уг+2,] + 4Уг+1,] — 3Уг,] \ +1 \3Уг,] — 4Уг-1,] +У-2,] \
где в — малый параметр.
Далее для центра и границ г-й ячейки уточняется индекс немонотонности:
(!ий (Уг+1,] ) + 1иА (Уг] ) + 1иА (Уг-1] ))
Inde =
3
IndL = (Ind (Y-1,J) + Ind (Yj)), IndR = (Ind (Y+1j) + Ind (Yj)).
Регулируется нижняя граница ограничителей:
Rmin = (1 - Inde ) Wmax + Inde Wmin.
Определяются значения ограничителей R (g, Ind) ,R (n, Ind) на левой и правой границах ячейки с номером i:
R(g = gi—i/2,j, Ind) = (1 - IndL) + IndLRmin,
R{C = gi+1/2j,Ind) = (1 - Indu) + IndRRmin.
Далее проводится минимизация значений ограничителей R (Ind):
RLim (Ind) = min (R (g = gi—1/2j, Ind) , R (g = gi+i/2j, Ind)) .
Рассчитывается среднее по i-му распределению значение ограничителя RLim (Ind)
N
RLim,i
R = _i_
Rcp N .
На основе следующего соотношения определяются окончательные значения ограничителей Я (£, 1пв) , Я (п, 1пд):
Я% 1п= Ш1П {Rcp, ЯЬ%т,%-2,, Яыт^-Ъ ЯЫт,г, RLгm,г+1, ЯЬггп,г+2} .
Отметим, что при проведении расчетов были использованы следующие значения параметров: Сгтах = 0,5, Сгтт = 0,01, Wmax = 0,8, ^"т1п = 0,5, п = 2.
Исходная дифференциальная система уравнений относительно временной переменной £ — это система обыкновенных дифференциальных уравнений первого порядка, которая может быть разрешена с помощью многошагового метода Рунге-Кутты (в данной работе использован трехшаговый вариант метода [24]). Для этого приведем исходную систему уравнений к нормальной форме с выделенной в
левой части временной производной
dUi
13
= Ь I иц I, где Ь — правая
31 " \ гз,
часть исходной системы, не содержащая производных по времени. Тогда трехшаговый вариант метода Рунге-Кутты реализуется в виде следующей последовательности шагов:
U(1) =
13
U(2) = 3 U(0)
U0 + ^tL(Ui3
(0)
13
4
13
U (3) = 1 U(0)
13 з 13
1
+ 4
2
+ 3
U3 + AtL (U1
1j
(1)
1j
U3 + AtL i U
(2)
1j
1j
Известно, что такой способ поиска решения и^ относительно времени £ позволяет обеспечить положительность искомых функций (если в момент времени Ьп решение является положительным, то оно остается положительным и в момент времени £п+1).
На втором временном дробном шаге £ € [£ + Д£/3,£ + 2Д£/3] явным методом с использованием краевых условий прилипания решается параболическая (вязкая) часть системы уравнений. При этом первые и вторые производные по пространственным переменным, входящие в рассматриваемую систему уравнений второго временного дробного шага, определялись с помощью изложенной компактной схемы шестого порядка точности.
На третьем временном дробном шаге £ € [£ + 2Д£/3, £ + Д£] с использованием неявного метода Розенброкка рассчитывается жесткая часть системы уравнений д—ш модели Кокли.
Шаг по времени Д£, необходимый для интегрирования приведенной выше дифференциально-разностной схемы, выбирается из условия выполнения критерия устойчивости Куранта-Фридрихса-Леви.
При решении уравнений переноса излучения применен модифицированный попеременно-треугольный метод с использованием
Рис. 1. Схематичное изображение расчетной области и структуры лазерного факела, возникающего при взаимодействии лазерного излучения с газовыми средами вблизи металлической преграды:
1 — металлическая мишень; 2 — область эрозионных металлических паров; 3 — контактная граница; 4 — область ударно-сжатого газа окружающей металлическую преграду среды; 5 — фронт ударной волны; 6 — поток лазерного излучения
трехслойной итерационной схемы, в которой итерационный временной шаг находят с помощью метода сопряженных направлений. В этом случае эллиптическое уравнение, соответствующее многогрупповому диффузионному приближению уравнения переноса излучения, решали на основе одного из наиболее быстрых итерационных методов — модифицированным попеременно-треугольным методом [26].
Результаты проведенных расчетов. Конкретные расчеты проведены для случаев воздействия импульсов лазерного излучения прямоугольной формы с длительностями 5... 500 нс. Значение полной энергии лазерного (С02-лазер) излучения составляло 0,1... 3 Дж, размер пятна фокусировки — около 0,12 см. В расчетах в качестве материала металлической преграды был принят алюминий; окружа-
ющая среда — воздух.
Расчетная область (рис. 1) при проведении двумерных расчетов в системе координат тг и ^п представляла собой прямоугольник. В нижней части рис. 1 показана плоская металлическая мишень 1, на поверхность которой перпендикулярно подводится лазерный поток и через плоскую поверхность которой эрозионный поток паров материала мишени втекает в расчетную зону. Сверху она ограничена прямой линией, на которой заданы невозмущающие условия для выходящего
д2 /
из расчетной области потока:
8x1
= 0, где f = {p,u,v,e}
и xn —
координата, нормальная к граничной поверхности. Пространственное положение данной прямой определяется из условия, чтобы возмущения численного решения, возникающие на верхней границе расчетной области от выходящего потока, не искажали течение вблизи металлической облучаемой преграды. С правой стороны область интегрирования ограничивается осью симметрии, на которой задаются соответствующие условия симметрии течения плазмы лазерного факела. С левой стороны располагается поверхность, находящаяся на достаточном удалении от оси симметрии, чтобы на ней можно было задавать граничные условия, соответствующие условиям на бесконечности в невозмущенной газовой среде.
Важно отметить, что физические вопросы, обсуждаемые в настоящей работе, рассматривались экспериментально и теоретически многократно. Обзор работ (приблизительно за 20 лет), посвященный затрагиваемой тематике, выполнен в работе [27].
В проведенных расчетах (по разработанной численной методике в одномерном и двумерном вариантах реализации) наблюдается ряд особенностей динамики образования лазерной плазмы. Одной из важных особенностей исследуемой задачи является то, что для отмеченного случая воздействия лазерного излучения наряду с газодинамическим механизмом переноса энергии потоком лазерной плазмы существенную роль играет механизм переноса внутренней энергии широкополосным излучением.
Далее рассматриваются [28] результаты расчетов динамики разлета и структуры течения эрозионной Al-лазерной плазмы в воздушной окружающей среде для значений давлений в ней: P = 1 атм, P = 10-1 атм и P = 10-2 атм.
В связи с разработкой эксимерных лазеров (с активной средой HgBr2 или XeCl) с относительно высоким КПД, мощностью qLaz = 108 Вт/см2 и длительностью импульса ~1мкс определенный интерес представляет исследование действия лазерного излучения УФ-диапазона с энергией кванта ELaz«0,1... 6 эВ на различные металлические преграды. В этом случае возникает вопрос об изменении характера поглощения (значит, и энерговыделения) лазерного излучения и плазмодинамических характеристик лазерной плазмы при переходе от лазеров с длинами волн Л =1 мкм к лазерам с длинами волн Л = 0,1 мкм, поскольку для последнего случая взаимодействия энергия квантов лазерного излучения становится сравнима с характерным значением первого потенциала ионизации материала мишени.
Для моделирования взаимодействия с мишенью лазерного излучения УФ-диапазона были проведены дополнительные одномерные расчеты, в которых истинное монохроматическое лазерное излучение заменялось широкополосным излучением постоянной (по времени и спектру) интенсивности в одном из спектральных диапазонов: hv1 = 0,1... 3,14 эВ (Л = 12,4... 0,395 мкм), hv2 = 3,14... 5,9 эВ (Л = 0,395... 0,207 мкм) (что соответствует видимому и ближнему УФ-диапазонам спектра).
На рис. 2 (см. 4-ю полосу обложки) и рис. 3-8 приведены одномерные и двумерные пространственные распределения температуры T, K, P, атм, числа M, продольной скорости -у,км/с, магнитного давления Pmag, атм, создаваемого спонтанными электромагнитными полями, соответствующие двум вариантам расчетов. В тексте указаны значения яркостных температур эрозионной лазерной плазмы, функциональные зависимости от времени t импульса отдачи Imp (отношения механического импульса, передаваемого его приемнику с учетом
Индекс 72781
Иллюстрации к статье В.В. Кузенова, *
«Математическое моделирование основных плазмодинамических характеристик в лазерном факеле вблизи алюминиевой мишени»
Число
Маха
2.6
2.52857 2.4571 Д 2.38571 2.31429 2.24286 2.17143 2.1
2.02857 1.95714 1.88571 1.81429 1.74206 1.67143 1.6
1.52857 1.45714 1.38571 1.31429 1.24286 1.17143 1.1
1.02857 0.957143 0.885714 0.814286 0.742857 0.671429 0.6
0.528571 0.457143
Число Маха
2.В
2 52857 2 4 5714 2 38571 2.31429 2 2428В 2.17143 2 1
2.02057 1 95714 1.88571 1 81429 1 7420В 1 В7143 1 .В
1 52857 1 45714 1 38571 1.31429 1,2420В 1 17143 1.1
1 02857 0.957143 0 885714 0.81428В 0 742857 0.В71429 0 Б
0.520571 0 457143
Рис.
2. Расчетная адаптивная сетка (а) и поле течения недорасширенной струи плазмы (б)
а
Рис.3. Распределения температуры Т (1), продольной скорости V (2), числа М (3) в лазерном факеле вдоль оси X в момент времени £ = 1,85 мкс при воздействии CO2-лазера, 4 — фронт контактной границы, 5 — фронт ударной волны:
a — Р^ = 1 атм, £ = 1,85 мкс; б — Р^ = 0,1 атм, £ = 1,75 мкс; в — Р^ = 1 атм, £ = 1,77 мкс; г — Р^ = 0,1 атм, £ = 0,92 мкс
противодавления окружающей среды, к энергии лазерного импульса) и силы Е, действующей на поверхность облучаемой мишени со стороны лазерного факела.
Первая группа одномерных расчетов, показанная на рис.3, 4, отвечает следующим параметрам лазерного воздействия: дЬаг = 2 х х 107 Вт/см2, ЕЬаг = 0,3 Дж, = 500 нс. Эта группа расчетов позволяет оценить влияние давления (Р^ = 1 атм и Р^ = 0,1 атм) в окружающей металлическую мишень газовой среде на плазмодина-мические характеристики в лазерном факеле.
Следует отметить, что одномерные расчетные исследования, проведенные на основе разработанных математической и численной моделей, имеют количественный смысл лишь тогда, когда характерный размер области лазерного факела, охваченного движением, значительно
Рис. 4. Распределение плотности потока лазерного излучения (7) и
широкополосного собственного излучения лазерного факела в спектральных группах Ни 1 (2) и Ни3 (3) в момент времени £ = 0,92 мкс при воздействии Ни2
меньше радиуса пятна облучения мишени. Такого рода расчеты также относительно правильно моделируют плазмодинамические процессы вблизи оси лазерного факела или в случае равномерного по поверхности мишени воздействия лазерного излучения. При этом последний случай воздействия лазерного излучения соответствует плоскому варианту разлета лазерного факела. Таким образом, одномерные расчеты дают возможность оценить влияние пространственного фактора на характеристики лазерной плазмы.
Результаты проведенных расчетов указывают на сложную внутреннюю структуру (зависящую от интенсивности лазерного излучения, термодинамических параметров окружающей газовой среды, тепло-физических свойств материала мишени) плазменного образования лазерного факела, формирующегося вблизи металлической мишени.
В случае рассмотрения первой группы одномерных расчетов (яьаг = 2 • 107 Вт/см2, ЕЬаг = 0,3 Дж, = 500 нс), отдельные результаты которых приведены на рис. 3, 4, структура лазерного факела может быть описана следующим образом:
• внешней границей, отделяющей плазменное образование от невозмущенной газовой среды, является ударная волна, которая на рисунках обозначена цифрой 5;
• вблизи поверхности мишени формируется волна разрежения, состоящая из слоя паров и плазмы материала преграды с температурой более низкой, чем температура последующего плазменного слоя;
• внутри плазменного образования располагается контактная граница, на рисунках обозначенная цифрой 4. Она отделяет плазму материала мишени от плазмы окружающей газовой среды. Очевидно, что плазмодинамические и спектрально яркостные параметры существенно изменяются при переходе через контактную границу.
Такой вид структуры лазерного факела объясняется одновременным протеканием в плазме лазерного факела нескольких противоположно направленных процессов:
• первый процесс связан с экранированием оптически плотными парами материала преграды лазерного излучения от поверхности мишени, что приводит к периодическому во времени характеру процесса испарения мишени;
• второй процесс связан с газодинамическим растеканием плазмы из области интенсивного поглощения лазерного излучения.
Так, из графических зависимостей рис.3,г следует, что экранирование лазерного излучения (т.е. частичный захват излучения лазера материалом металлической преграды) совместно с интенсивным газодинамическим течением и сбросом внутренней энергии в окружающую среду через контактную границу и на поверхность металлической преграды (в виде широкополосного излучения) приводит к колоко-лообразному распределению температуры приповерхностной плазмы внутри лазерного факела.
Результаты первой группы одномерных расчетов, приведенные на рис.3, 4, позволяют отметить, что явление лазерного пробоя наблюдается только для лазерного излучения с энергиями квантов в диапазоне Нр2 = 3,14 ... 5,98 эВ (Л = 0,395 ... 0,207 мкм — ближний УФ-диапазон), хотя испарение преграды было и в остальных вариантах расчетов. В этом случае внутри лазерного факела после оптического пробоя достигаются температуры до 20 кК (рис. 3, в, г). При этом скорость распространения внешней границы лазерного факела, вызванного облучением мишени С02-лазером, составляет щв = 1,35 км/с (для Рж = 1 атм) и -уув = 2 км/с (для Рж = 0,1 атм), что приблизительно в 2 раза меньше этой же скорости в случае воздействия на мишень Нр2-лазером. В этой ситуации относительно высокие значения температур в плазменном образовании, которое появилось в результате облучения Нр2-лазером мишени, могут быть объяснены сильным поглощением лазерного излучения внутри лазерного факела. При этом объемное энерговыделение на порядок превышает энергию, выделившуюся в результате воздействия С02-лазером. Так, если обратить внимание на рис.4 (кривая 1 — плотность потока лазерного излучения ), то можно увидеть практически 4-кратное уменьшение лазерного потока в пространственной области, располагающейся между поверхностью мишени и контактной границей и занятой плазмой материала мишени. Аналогичный расчет, выполненный для С02-лазера, дает более слабое поглощение лазерного излучения в той же пространственной области.
Из проведенных двумерных расчетов следует, что процесс расширения лазерного факела во времени условно можно представить в виде двух стадий:
• первая стадия соответствует случаю, когда весь объем лазерного факела движется от конденсированной преграды в сторону невозмущенной газовой среды (давление за ударной волной сопоставимо с давлением внутри объема плазменного образования). При этом первая стадия может быть представлена в виде начальной фазы — фазы формирования лазерного факела — и фазы затухания, когда ударно-волновые характеристики лазерного факела деградируют с течением времени;
• вторая стадия соответствует случаю, когда лишь ударная волна (отделяющая лазерный факел от окружающей невозмущенной среды) движется в сторону невозмущенной газовой среды, а большая часть объема плазменного образования схлопывается на оси системы координат гг. В этом случае давление за ударной волной существенно (приблизительно на порядок) превышает давление внутри объема лазерного факела, а динамический напор внутри лазерного факела невелик. Такая особенность в пространственном распределении давления приводит к появлению возвратного движения внутри плазменного образования (часть плазмы начинает двигаться в сторону облучаемой мишени и оси симметрии плазменного образования).
На рис.5,б-8,а показаны пространственные распределения газодинамических параметров, полученные в первой группе двумерных расчетов = 2 • 107 Вт/см2, = 0,3 Дж, = 500 нс, Рж = 0,1 атм, Рж = 1 атм). Сразу отметим, что в данной группе расчетов (в отличие от одномерных расчетов) на ранних фазах образования и расширения лазерного факела наблюдается оптический пробой.
На рис. 5, а, в показаны двумерные распределения, соответствующие начальной фазе первой стадии формирования лазерного факела и демонстрируют волновую структуру сильно недорасширенной импульсной струи паров материала преграды, истекающей в затопленное пространство, состоящее из плазмы воздуха и эрозионных паров.
В этой фазе расширения (см. рис.5,а, в) вблизи поверхности мишени на границе лазерного факела и окружающей среды формируется ударно-волновой комплекс, состоящий из двух ударных волн (УВ) — внешней УВ, распространяющейся в невозмущенной газовой среде, и внутренней УВ, двигающейся внутри плазменного образования лазерного факела), разделенных между собой контактной границей.
Из пространственных распределений, приведенных на рис. 5, а, в, следует, что при достаточно больших значениях степени нерасчетно-Ра
сти п = —— ^ 100 (Ра — давление на поверхности мишени в облаке»
сти испарения) вблизи края пятна облучения мишени в окружающем воздухе на начальной стадии взаимодействия возникает система бочкообразных УВ (см. рис.5,а, в), распространяющихся в направлении
см 0,6 -—
0,5 -
0,4 -
0,3 7
0,2 -
0,1 -
0
-0,4
Рис. 6. Пространственное распределение скорости V в лазерном факеле в момент времени £ = 2,14 мкс
от оси симметрии плазменного образования. Вследствие того, что поперечные размеры первой бочки растут как К/т^ах ~ у/п/(1а — 1), осуществляется только первая бочка, размеры которой существенно превышают радиус пятна облучения. На более поздних стадиях расширения лазерный факел приобретает эллиптическую пространственную форму (рис. 1, рис. 5, г). Из расчетов также следует, что скорости V продольного расширения (расширение лазерного факела вдоль координаты ^), полученные при проведении двумерных и одномерных расчетов для одинаковых условий в окружающей среде и параметрах лазерного воздействия, имеют приблизительно близкие значения.
Таким образом, можно констатировать, что на общую динамику процессов плазмообразования существенное влияние оказывают эффекты двумерности, которые, в частности, проявляются в сильном радиальном расширении лазерного факела.
Сильное радиальное расширение лазерного факела в центрированной волне разрежения, которая располагается на краю пятна облучения, сопровождается более интенсивным (по сравнению одномерными вариантами расчетов) газодинамическим выбросом легко ионизуемых (что облегчает их оптический пробой) паров материала преграды. Эти пары, интенсивно поглощая падающее на поверхность мишени лазерное излучение, увеличивают свою скорость в волне разрежения, двигаются в сторону от мишени. При этом они быстро нагреваются и ионизируются — происходит переход пара из состояния с температурой Т = 2... 5 кК в сильно поглощающую плазму (формируется волна ионизации) с температурой Т = 15... 30 кК (см. рис.5, а, в). По мере увеличения интенсивности падающего лазерного излучения фронт волны ионизации продвигается по направлению к поверхности (в область относительно холодных паров вещества мишени). При до-
стижении пороговых значений лазерного излучения, равных для большинства металлов 107 ... 108 Вт/см2 (проведенные двумерные расчеты подтверждают данное утверждение), ионизация начинается практически от поверхности мишени.
Отметим, что в этой группе двумерных расчетов светодетонацион-ный режим движения переднего фронта лазерного факела не реализуется. Это связано с тем, что плотность потока лазерного излучения относительно не велика (qLaz = 2• 107 Вт/см2). Газ из окружающей среды, располагающийся за фронтом УВ, нагрет до сравнительно низких температур (T~8 кК) и прозрачен для падающего на мишень лазерного излучения. В этой ситуации движение переднего фронта лазерного факела может быть описано волной дефлаграционного типа.
На рис. 1, рис. 5, б, г, рис. 8, a показаны двумерные распределения температуры T(r, z) и продольной скорости v(r, z), которые отвечают второй временной стадии формирования лазерного факела (стадии схлопывания).
На второй временной стадии расширения лазерного факела подвод энергии лазерным излучением уже отсутствует (t > tLaz), т.е. отсутствует источник энерговыделения внутри плазменного образования, а значит, уменьшается подкачка энергии в ударно-волновой комплекс, в результате чего он деградирует:
• внутренняя УВ превращается в волну сжатия;
• интенсивность внешней УВ падает;
• температура внутри плазменного образования уменьшается от T = 15... 30 кК на начальной стадии до T = 7... 9 кК на второй временной стадии (см. рис. 5, б, г);
• у поверхности металлической преграды возникает возвратное движение внутри лазерного факела (рис. 6, 5, в), т.е. наблюдается схло-пывание части плазменной области к оси лазерного факела и поверхности металлической мишени.
На рис. 7, а, б приведены двумерные распределения магнитного давления Рм. Сразу отметим, что давление Рм не оказывает никакого влияния на процесс расширения лазерного факела. Для цилиндрически симметричной формы лазерного факела напряженность спонтанного магнитного поля имеет одну компоненту B = Bизолинии которой совпадают с изолиниями магнитного давления Рм и токов проводимости. Из приведенных рисунков видно, что на начальной фазе расширения вблизи мишени на границе плазма — окружающая среда возникает тороидальная система токов, которая с течением времени отодвигается от поверхности конденсированной среды со скоростью приблизительно равной 3 км/с. Динамический процесс движения тороидальной системы токов определяется выносом плазменного образования лазерного факела и генерацией и диффузией спонтанного магнитного поля в проводящей среде факела (см. рис. 7, б).
Рис. 7. Пространственное распределение магнитного давления Рмаг, атм, в лазерном факеле в момент времени £:
а — Яьаг = 2 • 107 Вт/см2, ЬЬаг = 500 нс, Р^ = 1 атм, Ь = 0,46 мкс; б — дьаг = 2 х х 107 Вт/см2, ЬЬаг = 500 нс, Р^ = 1 атм, Ь = 2,14 мкс; в — дЬаг = 2 • 109 Вт/см2, = 50 нс, Р^о = 0,1 атм, Ь = 0,31 мкс
Отметим также, что процесс схлопывания внутри плазменного образования сопровождается появлением второго временного максимума в пространственном распределении магнитного давления Рм (рис. 7, а, б), которое создается увеличивающимся во времени (градиентом концентрации электронов и температуры) в лазерном факеле спонтанным электромагнитным полем.
На процесс деградации УВ, который наблюдается по мере увеличения времени ¿, негативным образом накладывается пространственный характер разлета лазерного факела. Поскольку к этому моменту времени вся выделившаяся энергия заключена внутри объема лазерного факела V, то увеличение объема V ~ г3 будет сопровождаться сильным падением давления Р ~ — (ре ~ Р, р ~ 1^). Это, в свою очередь, приведет к значительному уменьшению скорости расширения лазерного факела V ~ УР/р ~ 1/г по мере удаления границы, разделяющей плазму факела между окружающей средой и металлической преградой.
Из приведенных на рис.5,а, б, 8,а распределений видно, что на структуру течения и значения всех газодинамических параметров в лазерном факеле большое влияние оказывает давление Р^ в невозмущенной газовой среде.
Благодаря снижению внешнего давления от Р^ = 1 атм до Р^ = = 0,1 атм приблизительно в 2 раза увеличивается продольная скорость движения (вдоль оси 2) границы плазменного образования (от ■уУВ~1 км/с, для Р^ = 1 атм до ^УВ = 2 км/с, для Р^ = 0,1 атм).
В отличие от варианта (Р^ = 1 атм) для случая Р^ = 0,1 атм интенсивность внешней УВ близка к интенсивности сильной удар-'рув (7 + 1)1
ной волны I -~ --- . Наблюдается сильный вынос испарен-
(7 - 1)У
ного и ионизованного материала мишени в расчетную область (см. рис. 5, г). Причем наиболее интенсивное движение наблюдается вдоль оси 2, контактная граница приобретает характерную форму клина, острая часть которого направлена навстречу потоку лазерного излучения. Температура в плазменной области возрастает по сравнению с вариантом Р^ = 1 атм, дЬаг = 2 • 107 Вт/см2, ЕЬаг = 0,3 Дж, ЬЬаг = 500 нс и составляет Т « 30 кК для начальной фазы разлета. В варианте Р^ = 0,1 атм давление сильно неравномерно распределено по объему лазерного факела. Его максимум (Р « 20 атм) находится вблизи фронта УВ. При этом внутри плазменного образования давление существенно ниже (Р ~ 2 атм). Такая форма распределения давления по объему и относительно небольшой уровень динамического напора приводят к возникновению внутри лазерного факела возвратного движения (рис. 8, а), направленного от поверхности мишени к оси сим-
Рис. 8. Пространственное распределение скорости V, м/^ в лазерном факеле в момент времени £:
а — Цьах = 2 • 107 Вт/см2, ЬЬаг = 500 нс, Р^ = 0,1 атм, Ь = 1,95 мкс; б — = 2 х х 109 Вт/см2, гЬаг = 50 нс, Р^ = 1 атм, Ь = 1,86 мкс; в — дЬаг = 2 • 109 Вт/см2, Ььах = 50 нс, Роо = 0,1 атм, Ь = 0,31 мкс
метрик плазменного образования (начинается вторая стадия процесса расширения).
Яркостная температура Тя^ неравномерно нагретой приповерхностной плазмы является одним из важнейших интегральных оптических параметров, позволяющих оценить спектрально-яркостные характеристики лазерного факела. В проведенных расчетах в качестве исходной информации для расчета Тя^ служили плотности спектральных потоков излучения ^ (г = 1, 2,3) в трех спектральных группах.
При известном значении ^ (г = 1, 2, 3) яркостная температура Тя^ в соответствии с определением [19] находится из решения уравнения
9г = "ТаТя4г (^ (Ж2,г) -F (хМ)) , F (х^) = -—--, Х^. ^
П
exp (x) — 1' ' кТя,
о
где а — постоянная Стефана-Больцмана; i — номер частотной группы энергий фотонов, для которых проводится вычисление Тя^; x2ii и x\ti — соответственно верхняя и нижняя частотные границы i-й спектральной группы.
В первой группе расчетов (qLaz = 2 • 107Вт/см2, ELaz = 0,3 Дж, tLaz = 500 нс, Р» = 0,1 атм, Р» = 1 атм) излучение (яркостные температуры) лазерного факела максимально на первой стадии процесса расширения и заметно отличается от излучения черного тела (Тя1 = 22 кК, Тя2 = 14 кК, Тя3 = 12 кК, t и 0,2 мкс). При этом наибольшие значения яркостных температур (Тя1 = 22 кК) эрозионной лазерной плазмы достигаются в первой спектральной группе (hv1 = 0,1... 3,14эВ) в момент времени t и 0,2 ... 0,4 мкс.
Уровни значений импульса отдачи Imp и силы F максимальны на начальной стадии расширения и в этой группе расчетов составляют: Imp и 16 H/МВт, F и 10 H.
Вторая группа расчетов показана на рис. 7, в; 8, б, в; 9.
Эти пространственные распределения соответствуют лазерному воздействию с параметрами qLaz = 2 • 109 Вт/см2, ELaz = 3 Дж, tLaz = 50 нс и дают возможность оценить влияние энергии, мощности лазерного воздействия, длительности лазерного импульса. В этой группе расчетов в сравнении с предыдущим вариантом на порядок увеличивается энергия в лазерном импульсе (ELaz = 3 Дж). Увеличение энергии в лазерном импульсе сопровождается сильным (на порядок) уменьшением длительности лазерного воздействия (tLaz = 50 нс).
Изменение лазерных параметров, сделанное в данной группе расчетов, приводит на начальной фазе расширения (показана на рис. 9, б, г) к увеличению значений и степени равномерности температурного поля (Т и 30 кК для Р» = 1 атм; Т и 50 кК для Р» = 0,1 атм). Большая равномерность по объему поля температуры
Уровень Г, К 15 30000
14 13 12 11
23000 26000 24000 22000 10 20000 9 13000 16000 14000 12000 10000 3000 3000 4000 2000
г, см
Уровень Т,К 10 46000
41111.1
36222.2
31333.3
26444.4
21555.6
16666.7
11777.8 6888.89 2000
Г, см
Уровень V, м/с
14 13 12 11 10 9 8 7 6 5 4 3 2 1
800.0 600.0
400.0 200.0
108.1 36.4
0.8 0.0 -200.0 -400.0 -600.0 -800.0 -1000.0 -1200.0
г, см
0,4 г, см
Уровень V, м/с 20 36000 19 34000 18 32000 17 30000 16 28000 15 26000 14 24000 13 22000 20000 18000 16000 14000 12000 10896.5 10163.5 9884.19 8000 6000 4000 2000
12 11 10 9
Рис. 9. Пространственное распределение температуры и скорости в лазерном факеле в момент времени
а — I = 0,62 мкс, Рос = 1 атм; б — I = 1,86 мкс, Роо = 1 атм; в — I = 0,06 мкс, Роо =0,1 атм; г — I = 0,31 мкс, Роо = 0,1 атм
в основном связана с выравнивающим действием собственного широкополосного излучения лазерного факела. При этом внутри лазерного факела в области максимального энерговыделения (вблизи оси симметрии исследуемого плазменного образования) располагается пространственный максимум температуры Т, давления Р и продольной скорости v. Такой вид пространственных распределений основных газодинамических параметров заметно отличается от распределений, характерных для сильного сферически симметричного взрыва.
В этой группе расчетов продольная скорость v движения границы плазменного образования, направленная от металлической преграды, существенно возрастает (v и 2,5 км/с для Р» = 1 атм; v и 30 км/с для Р» = 0,1 атм) в сравнении с предыдущим вариантом расчетов (это явление, по-видимому, связано с увеличение КПД преобразования лазерного излучения во внутреннюю энергию плазменного образования).
Увеличение продольной скорости v, сопровождающееся ее слабым радиальным спадом, приводит к пространственной асимметрии лазерного факела:
• продольная координата границы лазерного факела существенно больше радиальной (особенно при Р» =0,1 атм);
• продольная скорость v границы лазерного факела больше массовой скорости плазмы, что связано с возникновением тепловой волны (прекурсорной области прогрева), распространяющейся перед контактной границей в лазерном факеле.
Вторая стадия процесса расширения лазерного факела отражена на рис. 9, б, рис. 8, б. Отметим, что для варианта лазерного воздействия qLaz = 2 • 109 Вт/см2, ELaz = 3 Дж, tLaz = 50 нс, Р» = 0,1 атм вторая временная стадия наступает значительно позже.
Как видно из приведенного на рис. 8, б пространственного распределения, большая часть объема лазерного факела (приблизительно от координаты Z и 0,3 см продольная скорость v имеет отрицательные значения — v < 0) движется к облучаемой металлической преграде. Максимальные значения температуры (T и 7 кК; z £ [0, 0,3] см, r £ [-0,2, 0,2] см) достигаются в зоне схлопывания лазерного факела.
Общий уровень яркостных температур, характерных для данной группы расчетов, выше значений, полученных в предыдущем варианте расчетов (qLaz = 2 • 107 Вт/см2, ELaz = 0,3 Дж, tLaz = 500 нс, Р» = 0,1 атм, Р» = 1 атм), отличается от спектра излучения черного тела — максимален на первой стадии расширения лазерного факела:
• Тя1 = 47 кК, Тя2 = 30 кК, Тяз = 18 кК, t и 0,2 мкс для Р» = 1;
• Тя1 = 70 кК, Тя2 = 70 кК, Тяз = 30 кК, t и 0,06 мкс для Р» = 0,1 атм.
Как видно из распределений, показанных на рис. 7, в, в этом варианте расчетов большое влияние оказывает давление Рх в невозмущенной газовой среде.
Снижение давления в окружающей среде до уровня Р^ = 0,1 атм приводит к значительному увеличению скорости движения контактной границы в направлении оси Z (v ~ 30 км/с для Р^ = 0,1 атм), увеличению значений термодинамических параметров (ТУВ & 40 кК) за УВ (т.е. к усилению УВ) и формированию эллиптической формы плазменного образования.
Как показали выполненные расчеты, в импульсных струях плазмы, создаваемой импульсным воздействием лазерного излучения на плоскую металлическую преграду в газовой среде, на границе лазерной струи и невозмущенного газа (в отличие от плазменных струй, создаваемых стационарным источником плазмы) отсутствует тороидальная по форме вихревая структура.
Однако в целом характер течения в приосевой области соответствует натеканию струи плазмы на деформируемую газовую преграду: в этой области наблюдается структура, состоящая из двух ударных волн, разделенных контактной границей.
На рис.7 приведены линии уровня магнитного давления Рм. Из распределений, показанных на этом рисунке, следует, что максимальное значение спонтанного магнитного поля наблюдается между контактной границей и ударной волной в области значительных градиентов температуры и плотности. В заключение отметим, что значения импульса отдачи Imp максимальны на начальной стадии расширения, больше аналогичных значений, полученных в предыдущей группе расчетов, и для данных вариантов составляют: Imp & 68 Н/МВт для Р^ = 1 атм; Imp & 113 Н/МВт для Р^ = 0,1 атм.
Заключение. Разработана математическая модель приповерхностного лазерного факела, основанная на уравнениях радиационной плаз-модинамики, записанных в произвольных криволинейных координатах. Численно исследованы радиационные и газодинамические процессы, возникающие в приповерхностной лазерной плазме при воздействии на металлическую мишень излучения СО2-лазера. Проведены расчеты всех основных газодинамических и излучательных параметров лазерного факела и металлической преграды. Выполнен предварительный анализ закономерностей образования и разлета плазменного образования. Отмечено, что функциональная зависимость от времени магнитного давления Рм (t), создаваемого спонтанными магнитными полями, имеет несколько максимумов.
Работа выполнена в рамках программы фундаментальных исследований ОЭММПУ РАН и гранта РФФИ № 07-01-00133.
СПИСОК ЛИТЕРАТУРЫ
1. А н и с и м о в С. И., И м а с Я. А., Романов Г. С., Ходыко Ю. В. Действие излучения большой мощности на металлы. - М.: Наука, 1970. - 272 с.
2. В е д е н о в А. А., Г л а д у ш Г. Г. Физические процессы при лазерной обработке материалов. - М.: Энергоатомиздат, 1985. - 208 с.
3. Kantrovitz A. Propulsion to orbit by ground-based lasers // Astronautics and Aeronautics. - 1972. - Vol. 10. No 5. - P. 74-75.
4. L a s e r s Initiated blast wave for launch vehicle propulsion // AIAA Paper 98-3735, 1998.
5. H e n d e r s o n D. B. Preprint LA-UR-&&-1442, LLNL., USA. 1977.
6. P i r r i A. N., S c h l i e r R., N o r t n a m D. Momentum transfer and plasma formation above a surface with a high-power CO2 laser // Applied Physics Letter. -1972. - Vol. 21. No 3. - P. 79-81.
7. P i r r i A. N., M o n s l e r M. J., N e b d s i n e P. E. Propulsion by absorption of laser radiation // AIAA Journal. - 1974. - Vol. 12, No 9. - P. 1254-1261.
8. ВоробьевВ.С. Пробой эрозионного факела при квазистационарном облучении металлов лазером // Квантовая электроника. - 1986. - № 9. - С. 1909-1911.
9. К у з е н о в В. В. Об одном алгоритме построения регулярных адаптивных сеток / Тез. докл. XV Междунар. конф. по вычислительной механике и современным прикладным программным системам (ВМСППС-2007): - Алушта. Крым, 25-31 мая 2007. - М.: Вузовская книга, 2007. - 544 с.
10. Альшина Е. А., Болтнев А. А., К а ч е р О. А. Эмпирическое улучшение простейших градиентных методов // Математическое моделирование. -2005. - Т. 17, № 6. -C. 43-57.
11. Головачев Ю. П. Численное моделирование течений вязкого газа в ударном слое. - М.: Наука, 1996. - 374 с.
12. С у р ж и к о в С. Т. Физическая механика газовых разрядов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006. - 639 с.
13. Захаров Н. С., Коробейников В. П., Ш айнога И. С. Численное моделирование процессов разлета и генерации магнитных полей плазменных факелов // ДАН. - 1988. - Т. 299, № 3. - С. 624-626.
14. Г о р б у н о в В. А., Никольская Л. С., П е т р у х и н А. И., П у ш т а р и к В. А., Р ы б а к о в В. А. Структура магнитного поля плазмы лазерного факела при низких плотностях потока излучения // Квантовая электроника. - 1984. - Т. 11, № 2. - С. 349-354.
15. Р е п и н А. Ю., С т у п и ц к и й Е. Л., Ш а п р а н о в А. В. Динамика и взаимодействие с преградой тороидального плазменного сгустка, ионизационно-динамические характеристики и электромагнитное излучение // ТВТ. - 2004. -Т. 42, № 4. - С. 523-537.
16. Суржиков С. Т. Тепловое излучение газов и плазмы. - М.:Изд-во МГТУ им. Н.Э. Баумана, 2004. - 543 с.
17. Knight C. J. Theoretical modeling of rapid surface vaporization with back pressure // AIAA Journal, 1979. - Vol. 17. No. 5. - P. 519-523.
18. Р а й з е р Ю. П. Физика газового разряда. - М.: Наука, 1992. - 535 с.
19. S u r z h i k o v S. T. Computing system for solving radiative gasdynamic problems of entry and re-entry space vehicles // Proc. of the 1st International Workshop on Radiation of High Temperature Gases in Atmospheric Entry; 8-10 October 2003, Lisbon, Portugal. ESA-533, December 2003. - P. 111-118.
20. C o a c l e y T. J. Turbulence modeling methods for the compessible Navier-Stokes equations // AIAA. - 1983. - No. 830-1693.
21. Толстых А. И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. - М.: Наука, 1990. - 230 с.
22. Численное решение многомерных задач газовой динамики / Под ред. Го-дуноваС.К. - М.: Наука, 1976. -400 с.
23. С а в е л ь е в А. Д. Численное моделирование течения в окрестности хвостовой части осесимметричного тела с истекающей струей // ЖВМ и МФ. - 2007. -Т. 47, № 2. - C. 310-320.
24. Shu G. W., O s h e r S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. II // J. Comput. Phys., 1984. - Vol. 83. - P. 32-78.
25. Yang J. Y.,LombardC.K. Uniform second order accurate ENO schemes for the Euler equations of gas dynamics // AIAA Pap. - 1987. - No 1166. - 9 p.
26. Самарский А. А., Николаев Е. С. Методы решений сеточных уравнений. - М.: Энергоатомиздат, 1978. - 589 с.
27. S u r z h i k o v S. T., К u z e n o v V. V., C a p i t e 11 i M., C o l o n n a G. Numerical analysis on near-surface laser plasma in gases and vacuum // 44th AIAA Aerospace Sciences Meeting and Exhibit. 9-12 January 2006, Reno, Nevada, AIAA 2006-1174.
28. S u r z h i k o v S. T., К u z e n o v V. V., P e t r u s e v A. S. Radiation gas dynamics of aluminium laser plume in air // AIAA. - 2008-1108. - P. 1-8.
Статья поступила в редакцию 12.03.2009
Виктор Витальевич Кузенов родился в 1956 г., окончил в 1983 г. Московский государственный университет им. М.В. Ломоносова. Канд. техн. наук, старший научный сотрудник лаборатории "Радиационная газовая динамика" Института проблем механики им. А.Ю. Ишлинского РАН, доцент кафедры "Теплофизика" МГТУ им. Н.Э. Баумана. Автор более 70 научных работ в области теплофизики и радиационной газовой динамики.
V.V. Kuzenov (b. 1956) graduated from the Lomonosov Moscow State University in 1983. Ph. D. (Eng.), senior researcher of "Radiation Gas Dynamics" laboratory of A.Yu. Ishlinskii Institute for PGroblems of Mechanics, RAS, assoc. professor of "Thermal Physics" department of the Bauman Moscow State Technical University. Author of more than 70 publications in the field of thermal physics and radiation gas dynamics.
