Исследование термоупругих эффектов в металлах в рамках модифицированной модели термического пика Текст научной статьи по специальности «Физика»

УДК 519.624.3

Исследование термоупругих эффектов в металлах в рамках модифицированной модели термического пика

И. В. Амирханов, И. В. Пузынин, Т. П. Пузынина, И. Сархадов

Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д.6, Дубна, Московская область, Россия, 141980

На основе модифицированной путём добавления уравнений движения модели термического пика численно исследованы термоупругие волны, возникающие в образце никеля под воздействием тяжёлого урана высокой энергии 700 МэВ. Показано, что возникающие при этом максимальные напряжения термоупругих волн на два порядка превышают предел прочности материала на оси трека.

Определены временные зависимости температур электронного газа, кристаллической решётки и термоупругих напряжений для разных расстояний от оси трека. Описана также динамика изменения максимальной области, внутри которой термоупругие напряжения превышают предел прочности материала. По мере удаления от оси трека по радиусу амплитуды термоупругих волн уменьшаются обратно пропорционально радиусу.

Ключевые слова: термоупругость, численное моделирование, деформация, напряжение, модель термического пика.

1. Введение

При прохождении тяжёлой заряженной частицы через вещество вблизи траектории возникают локальные перегревы, приводящие к генерации термоупругих напряжений [1], которые могут превышать прочность материала. Изменение температуры материала происходит не только за счёт подвода тепла от внешних источников, но также и за счёт самого процесса деформации. При деформации возникают упругие и тепловые волны. Эти явления обычно исследуются в рамках теории термоупругости [2]. В модели термического пика тепловые эффекты описываются системой двух связанных уравнений для температур электронного газа и кристаллической решётки Те, Т.

ВТ

Се= У(АеУТе) - 9(Те - Т.) + Ае(г, 1), (1)

В Т

С^ = У(АгУТг) + д(Те - Т.) + Аг(г, 1), (2)

где С7, А7, А7; 7 = е,% — соответственно теплоёмкости, теплопроводности и функции источников для электронного газа и кристаллической решётки, д — константа взаимодействия электронной подсистемы с решёткой.

2. Постановка задачи

В работе [3] численно исследована модель эволюции термоупругих волн, возникающих под воздействием импульсного пучка ионов. В настоящей работе приведём исходные уравнения термоупругости в рамках модели термического пика.

Статья поступила в редакцию 18 декабря 2012 г. Работа выполнена при финансовой поддержке грантов РФФИ, № 12-01-00396-а, 11-01-00278-а.

Функция свободной энергии деформированного тела имеет вид [2]:

1 \ ^ КК 1 (ди ди \

Р = Ро(Т) - Ка(Т - То)иц + ц\щк - 3 ё^щА + —и2Ь щк = + -д^) '

Е 1 1

К = —--, ц = —--, — 1 < а < —,

3(1 - 2а) 2(1 + а) 2'

где а — число Пуассона, К, ц — называются соответственно модулями всестороннего сжатия и сдвига, Е — модуль Юнга, щ — 1-я компонента вектора перемещения и, а = 3аь — коэффициент объёмного теплового расширения (а^ — коэффициент линейного теплового расширения), 5iк — символ Кронекера, Ро, То — соответственно свободная энергия и начальная температура тела в неде-формированном состоянии. По заданной свободной энергии определяются тензор напряжений и энтропия

др / 1 \

а^ = дик = -Ка(Т - То)бгк + 2ц (щк - 3ккщА + КщАк,

5 = - дР = во(Т ) + Каии, д Т

где Я0(Т) — энтропия тела в недеформированном состоянии. По агк и 5 составляются уравнения движения и теплопроводности:

д2иг дац Т/Г дТ 1 ^ 1 ц \ д2щ 1 д2иг

д Uj = oai£ = _Ka<±L + (k + I) 1 + dt2 dxj dxj V 3 J dxjdxi

3/ dxjdxi dx2 ' rp°S fdSodT дщЛ dSo Cv

Tm =T\dT~dt + КаЖ) = V(XVT) + Q, дт = T-

Здесь Cv — теплоёмкость тела при постоянном объёме, Q — функция источника. Повторяя эту процедуру для второго уравнения системы (1)—(2) и учитывая, что Q = A(r, t) + g(Te — Tj), приходим к следующей системе уравнений:

д T

= V(\eVTe) — g(Te — Tj) + Ae(f, t), д T д

Cj-^ + KTja—div U = V(XjVTj) + g(Te — Tj) + Aj(r, t), pdU = —KaVTj + (К + |) Vdiv U + iV2u.

При решении этой системы уравнений необходимо наложить на функции Te, Tj, U начальные и граничные условия (исходя из физической постановки задачи).

Исходную систему можно решать, вводя скалярный и векторный потенциалы Ф и Ф, посредством U = grad Ф + rot Ф. При этом исходная система принимает вид

д T

Ce^f- = V(XeVTe) — 9(Te — Tj) + Ae(f, t), д T д

Cj~d^ + KTj^di АФ = V(^jVTj) + g(Te — Tj) + Aj(r, t), p-дОф = —Ka(Tj — To) +(k + АФ,

в2ф Л =, р~д¥ = аДФ'

Последние уравнения являются волновыми уравнениями продольных и поперечных волн в твёрдом теле. Как видно, в первом волновом уравнении присутствует источник продольных волн, который определяется тепловым расширением тела, а во втором уравнении источника нет. Поскольку в дальнейшем мы рассматриваем случай, когда тело в начальном состоянии находится в недеформирован-ном состоянии, и второе уравнение для этого случая имеет тривиальное нулевое решение, последнее волновое уравнение исключим из рассмотрения. Кроме этого, теплофизические параметры электронного газа и ионной кристаллической решётки полагаем константами, равными их значениям при комнатной температуре. Тогда исходная система в цилиндрической системе координат с учётом аксиальной симметрии, пренебрежением теплового потока по оси z и с добавлением начальных и граничных условий принимает вид:

ВТ

Се(То)-^ = Ле(То)ДТе - д(То)(Те - Тг) + Ле(г, г), ВТ В

Сг(То)ВТ + КТга-ДФ = Лi(То)ДТг + д(То)(Те - Тг) + Аг(г, В В

В 2Ф ( 4а \ В2 1 В

= -КЫЪ - То) + (л- + Дф; Д = В-2 + 1 - .

Ае(г, ¿) = В ехр(-г /г о) ехр

(I - Ыо)

2а2

Яе(0), Яе(0) = ( -

' х=о

-9 „ _ _ + _ 1П-15

( ВЕЛ

V Вг А

А^г, г) = 0, Го = 10-9 м, а1 = ¿о = 10-15 с

/ВФ \

Ф(Г, 0) = 0, - = 0, Те,г(Г, 0) = То,

<=о

С^Ь = 0-Т.. (Г- , = То, (£) = =( ^

агг = -Ла(Т - То) +(л - ^ ДФ + 2^?Ф

3 у ВГ 2 '

а,, = -Ла(Т -То)+ (л - |) ДФ + ^ВГ.

Константа В определяется из условия нормировки:

СЮ сю

12пгёг I А(г, <)^ = 5е(0); В

Здесь агг, а^ — ненулевые компоненты тензора термоупругих напряжений, возникающие под действием падающего тяжёлого заряженного иона на материал, ве (0) — ионизационные потери энергии иона на единицу длины при входе в металл.

Вводя безразмерные величины

t _ г Ф ^ ТР - Т агГ _ а

ДГ ■ Дг' * фо' ^ То'" То,аг"' ао

с = —, г = —, Ф = —, Т = —, Т = —, агг = —, а^

2

0

тах

1

исходную систему запишем в безразмерной форме (Дг = 100го, Д£ = 100£о).

При этом черту безразмерных величин, кроме безразмерной температуры кристаллической решётки, опускаем. Тогда исходная система принимает следующий вид:

дТ = кеАТ - де(Т -Т) + А(г, г), (3)

д Т д

— + РТ-АФ = кгАТ + дг(Т -Т), (4)

дг2 У ( ); дг2 + г дг' (5)

А(г, г) = Ао ехр(-100г) ехр[-5000(г - 0, 05)2], (6)

/дФ \

Ф(г, 0) = 0, — =0, Т(г, 0)=Т(г, 0) = 1, (7)

) 4=0

г=0

/дФ\

\дг

=1^1 =0,

д

(8)

агг = -С(Т - 1) + РАФ + дф, (9)

1 дФ

а^ = -С(Т - 1) + РАФ + - — . (10)

Здесь константа Ф0 = а0Аг2/2ц. При этом константы к1, д1 (7 = е, г), Р, С, Б, Ао, 2, Р определяются формулами

(ЩАг = д(Т)АЬ = . = Каа0 = КаП С7Аг2 , 91 = С1 , 7 = ^ Р = 2цСг , С = а0 ,

Б = 2цС Аг2 = ВАгБе(0) 2 = (3К + 4ц)Аг2 = 1

= рАг2 , 0 = Се(Т0)Т0 , у= ЗрАг2 , =2ц - 3'

Параметры модели Се, С», д, Хе, Лi,Sе(0), р, Е, а взяты при комнатной температуре Т0 = 300 К для никеля и имеют следующие значения (а = 0, 3; а0 = 2,2 ■ 106):

Се = 2, 73340-2 Дж/(см3■К), С = 3, 91 Дж/(см3■К), д = 4, 05-1012 Вт/(см3■К),

Ле = 4,1 Вт/(см ■ К), А» = 0,91 Вт/(см ■ К), 5е(0) = 1,157 ■ 10-7 Дж/см,

р = 8, 9 г/см3, Е = 20,2 ■ 1010 Па, а = 54 ■ 10-6 К-1.

При приведённых величинах безразмерные параметры к7, д7, (7 = е, г), Р, С, Б, А0, V2, Р, гШах принимают значения:

ке ~ 0,15; кг ~ 2,33 ■ 10-4; де ~ 14, 82; д1 ~ 0,1036; р ~ 3, 291 ■ 10-5; С ~ 1239, 54; Б ~ 2,164 ■ 10-2; А0 ~ 8,97 ■ 106; V2 ~ 30, 55 ■ 10-6; Р = 0, 75; гшах = 2.

Численно решая систему (3)—(6) с начальными и граничными условиями (7)—

(8), находим температуры Т(г, г), Т(г, г) и функцию Ф(г, г), а затем по формулам

(9), (10) компоненты напряжения агг, аПоскольку при а = 0 (число Пуассона) Р = 0 из последней формулы видно, что для этого случая с учётом граничного

условия, компонента напряжения а^ в точке г = 0 изменяется прямо пропорционально температуре. В работе [4] нами вычислена временная зависимость температуры в разных точках без учёта влияния упругих деформаций на температуру. Если считаем данное влияние слабым, то используя искомую зависимость температуры в точке = 0, можем оценить порядок термоупругих напряжений в оси трека. Поскольку С ~ 1000, Ттах(0, €) ~ 30, а^ ~ 30000 х ао ~ 6 ■ 1010, т.е. на оси трека термоупругое напряжение на два порядка превышает предел прочности материала апр ~ 108 и материал разрушается под действием термоупругого напряжения.

В настоящей работе проводим численное исследование временной динамики образования упругих деформаций и их взаимное влияние с температурой кристаллической решётки.

3. Метод численного решения

Введём равномерную сетку по переменным г, £ в уравнениях (3)—(10), то есть положим: { гк = кТг (к = 0,1, ••• ,т), tj = (,] = 0,1, ••• ,п)}; кг и — соответственно шаги по переменным г и ^ А- = А(г-, Р). Для решения системы уравнений (3)—(10) будем использовать следующую явную конечно-разностную схему порядка аппроксимации 0(Т + Л^) [5-7]:

Р + 1 ОЛj I 1

= у2ЛкФj -"(Т- - 1); к = 0,1, ••• ,т; (11)

Ф-+1 - 2Фк + Фк 2 —-Т-- = Лкф - "(Ч

Т1+1- Т

к Т-- = кеЛкТ - де(Ц -Т-) + А—; (12)

т

с j + 1 с

2Т

^к+1 - 2^к + ^к—1 . 1 ^к+1 - ^к-1,

с Г - ■ В - ■

-- + Ч&Лк(Ф,+1 - Ф,_-1) = кЛТ + дг(Т- - Т,); (13)

к2 + г- 2ТГ '

Ф—1 = Ф- = 0, ТГ-? = Т- = 1; к = 0,1, ••• , т; (14)

4ф1 - ф2 - 3фо =0 4Т, -Т2 - ЗТо =0 2Т 2Т

= 0; ^1,••• ,»;

4Ф^_1 - Ф^_2 - 3Фт 4Т^-1 - Т^-2 - ЗТ^ =0 2Т

2Т

г

4Т? _ Т7 — ЗТ

т— 1 2„г_2 О^т = 0

(15)

(16)

При решении дискретной системы уравнений (11)—(13) с начальными условиями (14) и граничными условиями (15)—(16) сначала решим первое уравнение относительно Ф-+1, а потом решим остальные уравнения относительно Т-+1 и

Тк.

При вычислении по данной схеме шаг по времени Т выбирается из условия устойчивости схемы [5—7]

ке^/И2г < 1/2. (17)

4. Обсуждение численных результатов

На рис. 1 представлены динамика расширения области до времени г = 40 ( = 10-12 с), внутри которой термоупругие напряжения превышает предел прочности материала, а также профили температуры кристаллической решётки и термоупругих напряжений за время = 40.

Рис. 1. Временная зависимость радиуса максимальной области: а — внутри которой напряжения термоупругих волн превышают предел прочности материала; б — профиль температуры кристаллической решётки; в, г — профиль термоупругих напряжений в момент времени £ = 40

На рис. 2 приведены временные зависимости температуры кристаллической решётки и термоупругих напряжений для разных расстояний от оси трека.

Рис. 2. Временная зависимость на разных расстояний от оси трека (1 — г = 0; 2 — г = 0, 05; 3 — г = 0,1; 4 — г = 0,15) для образца никеля облучаемого ионом урана: а — температуры кристаллической решётки; б — термоупругих

напряжений агг (г, ¿)

Температура кристаллической решётки превышает температуру плавления и испарения материала, т.е. вокруг трека иона может происходить плавление и испарение вещества. Кроме того, образующиеся термоупругие напряжения вокруг трека иона намного превосходят предел прочности материала, т.е. ион высокой энергии разрушает материал образца.

Максимальная область плавления, вычисленная в работе [4], намного меньше, чем области разрушения материала, вычисленные в настоящей работе, т.е. вклад термоупругих напряжений в разрушение материала больше вклада фазовых превращений.

5. Заключение

На основе модифицированной модели термического пика численно исследованы термоупругие волны, возникающие в образце никеля под воздействием тяжёлого урана высокой энергии 700 МэВ. Показано, что возникающие при этом максимальные напряжения термоупругих волн на два порядка превышают предел прочности материала на оси трека.

Определены временные зависимости температур электронного газа, кристаллической решётки и термоупругих напряжений для разных расстояний от оси трека. Описана также динамика изменения максимальной области, внутри которой термоупругие напряжения превышают предел прочности материала. По мере удаления от оси трека по радиусу амплитуды термоупругих волн уменьшаются обратно пропорционально радиусу.

Следует подчеркнуть, что при численном исследовании термоупругих эффектов с помощью упрощающих предположений задача приводилась к одномерной по пространству. Это позволяло получить качественные и количественные выводы.

Однако для более точной оценки термоупругих эффектов искомую систему лучше решать в двумерной по пространстве цилиндрической системе координат с учетом аксиальной симметрии. При этом источники в каждом из уравнений должны быть написаны с учетом потери энергии иона в обеих подсистем. Кроме этого, ион, двигаясь внутри мишени, потеряет свою энергию не мгновенно, а в течение времени своего движения до самой остановки. Учёт этих эффектов будет предметом нашего дальнейшего исследования.

Литература

1. Калиниченко А. И., Лазурик В. Т. Возбуждение акустических колебаний пучками заряженных частиц малой плотности // ЖЭТФ. — 1973. — Т. 65. — С. 2364-2368. [Kalinichenko A. I, Lazurik V. T. Vozbuzhdenie akusticheskikh kolebaniyj puchkami zaryazhennihkh chastic maloyj plotnosti // ZhEhTF. — 1973. — T. 65. — S. 2364-2368. ]

2. Ландау Л. Д., Лифшиц Е. М. Теория упругости. — 1970. — С. 248. [Landau L. D, Lifshic E. M. Teoriya uprugosti. — 1970. — С. 248. ]

3. Численное моделирование влияния вязкости на распространение термоупругих волн в образце металла, облучаемого импульсными пучками ионов / И. В. Амирханов, Е. В. Земляная, И. В. Пузынин и др. // Поверхность. — 2004. — № 10. — С. 92-99. [Chislennoe modelirovanie vliyaniya vyazkosti na rasprostranenie termouprugikh voln v obrazce metalla, obluchaemogo impuljsnihmi puchkami ionov / I. V. Amirkhanov, E. V. Zemlyanaya, I. V. Puzihnin и др. // Poverkhnostj. — 2004. — No 10. — S. 92-99. ]

4. Численное исследование температурных эффектов в материалах при облучении их тяжелыми ионами высоких энергий в рамках уравнений теплопроводности для электронов и решетки / И. В. Амирханов, А. Ю. Дидык, Е. В. Земляная и др. // Письма в ЭЧАЯ. — 2006. — Т. 3, № 1(130). — С. 63-75. [Chislennoe issledovanie temperaturnihkh ehffektov v materialakh pri obluchenii

ikh tyazhelihmi ionami vihsokikh ehnergiyj v ramkakh uravneniyj teploprovodnosti dlya ehlektronov i reshetki / I. V. Amirkhanov, A. Yu. Didihk, E. V. Zemlyanaya и др. // Pisjma v EhChAYa. — 2006. — T. 3, No 1(130). — S. 63-75. ] 5. Самарский А. А. Теория разностных схем. — М.: Наука, 1983. — 616 с. [Samarskiyj A. A. Teoriya raznostnihkh skhem. — M.: Nauka, 1983. — 616 s. ]

UDC 519.624.3

Investigation of Thermoelastic Effects in Metals in the Frame of the Modified Thermal Spike Model

I. V. Amirkhanov, I. V. Puzynin, T. P. Puzynina, I. Sarhadov

Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russia, 141980

A numerical investigation of a thermoelastic waves arising in a nickel sample under the action of the uranium heavy ion high energy 700 MeV is made on the basis of a thermal spike model which is modified by adding to it move equations. It is shown that the arising maximal strength of the thermoelastic waves exceeds hundred times ultimate strength of a samples material on the axis of a track.

The time dependence of the electron gas and ion lattice temperatures, of the thermoelastic waves strength are determined at a different distance from center of track. The time dynamic change is also determined of the maximal region within which thermoelastic strength exceeds ultimate strength of a samples material. By moving away from the center of a track along the radius, the amplitude of the thermoelastic waves will decrease inversely proportional to the radius.

Key words and phrases: thermoelasticity, thermoelastic waves, numerical simulation, deformation, stress, thermal spike model.