CC BY
32
УДК 539.3, 538.913
О динамике материала с изменяющейся микроструктурой
Н.Ф. Морозов1,2, Д.А. Индейцев1,2,3, Б.Н. Семенов1,2,3, С.А. Вакуленко1, Д.Ю. Скубов13, А.В. Лукин3, И.А. Попов3, Д.С. Вавилов1,4
1 Институт проблем машиноведения РАН, Санкт-Петербург, 199178, Россия 2 Санкт-Петербургский государственный университет, Санкт-Петербург, 199034, Россия 3 Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, 195251, Россия 4 Военно-космическая академия им. А.Ф. Можайского, Санкт-Петербург, 197198, Россия
В многочисленных экспериментальных работах по ударно-волновому нагружению металлов с помощью электронной микроскопии было обнаружено, что в определенном диапазоне скоростей ударника в материале могут происходить преобразования его кристаллической структуры. На макроуровне эти изменения находят отражение в виде энергетических потерь, связанных с образованием новой структуры, которые проявляются на временном профиле скорости тыльной стороны поверхности образца, содержащего ключевую информацию о свойствах материала. В данной статье с целью описания структурных превращений предложена двухкомпонентная модель материала с нелинейной силой внутреннего взаимодействия, учитывающая его периодическую структуру. Динамические уравнения записаны относительно перемещения центра масс компонентов, выступающего в роли измеряемого макропараметра, и их относительного смещения, являющегося внутренней степенью свободы, отвечающей за структурные преобразования. В рамках предложенной модели решена квазистатическая задача о кинематическом растяжении двухкомпонентного стержня, которая позволила определить параметры, обеспечивающие немонотонную зависимость напряжения от деформации, часто используемой при описании материалов, подверженных фазовым превращениям. В результате решения динамической задачи о нестационарном воздействии на материал, задающемся в виде короткого прямоугольного импульса, показан эффект гашения нестационарной волны, связанный с диссипацией ее энергии в изменение его структуры. На основе континуально-дискретной аналогии было получено аналитическое выражение, позволяющее оценить длительность процесса структурных преобразований и параметр, характеризующий силу внутреннего взаимодействия между компонентами. Сделанные выводы подтверждаются численным решением нелинейной задачи Коши, выполненным методом конечных разностей.
Ключевые слова: двухкомпонентная модель, динамика структурных преобразований, нелинейная сила взаимодействия, метод конечных разностей
On the dynamics of the material with changing microstructure
N.F. Morozov12, D.A. Indeitsev123, B.N. Semenov123, S.A. Vakulenko1, D.Yu. Skubov13, A.V. Lukin3, I.A. Popov3, and D.S. Vavilov14
1 Institute for Problems in Mechanical Engineering RAS, St. Petersburg, 199178, Russia
2 Saint Petersburg State University, St. Petersburg, 199034, Russia 3 Peter the Great St. Petersburg Polytechnic University, St. Petersburg, 195251, Russia 4 A.F. Mozhaysky Military-Space Academy, St. Petersburg, 197198, Russia
Numerous experimental studies on shock wave loading of metals have shown by electron microscopy that the crystal structure of the material can undergo transformation in a certain range of striker velocities. At the macroscale, these changes are observed as energy losses associated with the formation of a new structure. The losses are manifested on the temporal velocity profile of the back side of the sample surface which contains key information about the material properties. In this paper, a two-component model of a material with a nonlinear internal interaction force is proposed for the description of structural transformations, taking into account the periodic structure of the material. Dynamic equations are written with respect to the motion of the center of mass of the components acting as a measured macroparameter, as well as with respect to their relative displacement which is the internal degree of freedom responsible for structural transformations. The proposed model is applied to solve a quasi-static problem of the kinematic extension of a two-component rod in order to determine the parameters of a nonmonotonic stress-strain relationship that is often used in describing materials subjected to phase transformations. By solving a dynamic problem of nonstationary influence on the material by a short rectangular pulse, the effect of nonstationary wave damping is demonstrated which is associated with the wave energy dissipation in structural changes of the material. An analytical expression is obtained on the basis of a continuous-discrete analogy for estimating the duration of structural transformations and the parameter characterizing the force of internal interaction between the components. The conclusions are confirmed by a numerical solution of the nonlinear Cauchy problem within the finite difference framework.
Keywords: two-component model, structural transformation dynamics, nonlinear interaction force, finite difference method
© Морозов Н.Ф., Индейцев Д.А., Семенов Б.Н., Вакуленко С.А., Скубов Д.Ю., Лукин А.В., Попов И.А., Вавилов Д.С., 2017
1. Введение
На сегодняшний день одной нз важных проблем механики деформируемого твердого тела является исследование материалов с изменяющейся внутренней структурой при ударно-волновом нагружении, когда сложная перестройка кристаллической решетки происходит за чрезвычайно малые времена, не превышающие несколько микросекунд [1-3]. Проведенные эксперименты по высокоскоростному деформированию показали, что в определенном диапазоне скоростей ударника в результате прохождения волны по материалу в нем формируются устойчивые сетчатые образования размером 0.1-0.3 мкм, не исчезающие после снятия внешнего напряжения [4] и обладающие более высокой микротвердостью по сравнению с «чистым» материалом, что способствует увеличению его откольной прочности. На макроуровне процесс структурообразования проявляется в виде энергетических потерь, которые могут быть зафиксированы на временном профиле скорости тыльной стороны поверхности образца, содержащем ключевую информацию о динамических свойствах материала. Типичный временной профиль скорости металлического образца изображен на рис. 1 [5]. Он состоит из четырех основных участков: кривая ОА представляет собой упругий предвестник, за которым начинает развиваться пластический фронт АВ, далее зависимость выходит на плато ВС, где скорость остается почти постоянной, и, наконец, за ним следует фронт разгрузки (участок CD). При этом волновые сопротивления мишени и ударника выбираются так, что максимальное значение скорости поверхности и^, достигаемое в точке В, равно скорости ударника и1тр. В том случае когда в материале формируются новые структуры, данное соотношение не выполняется и образуется так называемый дефект скорости [6] Аи = и1тр - и& (рис. 2), показывающий, что часть энергии уходит на структурные преобразования. Процесс перестройки обычно сопровождается быстрыми осцилляциями на пластическом фронте (рис. 3), возникающими сразу после прохождения упругого предвестника, длительность которых не превышает 0.15-0.30 мкс. Они свидетельствуют о том, что
§ 2 зоо-§ £ ю а о ад 8 | 200Л О
0......
0 О 400 800 1200
Время, не
Рис. 1. Типичный временной профиль скорости
перестройка структуры представляет собой сложный динамический процесс, затрагивающий несколько масштабных уровней [7]. Отметим, что данный результат был получен для широкого класса металлических материалов, что указывает на необходимость создания моделей, учитывающих возможность перехода материала в новое состояние. Часто встречающаяся в литературе модель [8-11], в которой используется предположение о немонотонной связи между напряжением и деформацией, не подходит для описания динамики быстропротекаю-щих процессов при ударно-волновом нагружении, поэтому в данной работе рассматривается иной подход, основанный на механике многокомпонентных сред.
2. Двухкомпонентная модель
В этом случае твердое тело представляет собой совокупность N взаимодействующих континуумов, для каждого из которых в любой точке объема определены плотности рг- и скорости VI (i = 1, 2, ..., N [12-14]. Для простоты ограничимся рассмотрением только двух компонентов, представляющих собой кристаллические решетки, связанные нелинейной силой взаимодействия. Основная цель работы состоит в том, чтобы, используя данное представление о материале, показать возможность передачи энергии на микроуровень, которая в эксперименте выражается в виде энергетических потерь, приводящих к образованию дефекта скорости (рис. 2). Таким образом, не ставя перед собой задачу описания формы временного профиля скорости тыльной поверхности образца (рис. 1), мы не будем принимать во внимание пластические свойства материала и предположим, что каждый из компонентов подчиняется закону Гука. Тогда для одномерного случая их динамические уравнения принимают вид [15]
Рис. 2. Дефект скорости
Рис. 3. Осцилляции на пластическом фронте
где ui (i = 1, 2) — перемещение каждого из компонентов; Et — модуль Юнга; pi 0 — плотность в равновесном ненагруженном состоянии. Через R обозначена сила взаимодействия, которая должна учитывать возможность перехода материала из одного состояния в другое. Предположим, что ее аналитическое выражение состоит из двух слагаемых, первое из которых обеспечивает наличие нетривиального положения равновесия и отражает периодичность структуры, а второе представляет собой диссипативную составляющую
R = K sin (Az) + v z, (2)
где z = u1 - u2 обозначает относительное перемещение; параметр K определяет максимальное значение силы взаимодействия; коэффициент v характеризует диссипацию. Величина А = 2я/d обратно пропорциональна периоду кристаллической структуры.
Уравнения (1) записаны относительно двух континуумов, но экспериментально измеряемым параметром является некоторая средняя величина, вклад в которую вносится обоими компонентами. Если считать, что этой величиной, содержащей информацию о физическом состоянии системы, является перемещение центра масс U = (P10u1 + p20u2)/(p10 + р20), то уравнения двухком-понентной среды удобно переписать относительно центра масс U и разности смещений z = u1 - u 2 [16]: Э 2U 1 Э 2U Э2 z
Эх2 c2 dt2 дх2
д2 z 1 Э2 z
Э 2U
(3)
r-2 - "2 -dtT = вR(z, z) + Y^2
Эх c„ dt
где
2 = E1 + E2
С2 =
dt2
E1 E2(P10 +P20)
P10 +P ю' * (E1 + E2)P10P20* Параметры а, в и у определяются физическими свойствами материала:
а = -
Y =
E2p10 - E1p20 (E1 + E2)(P10 +P 20): E2p10 - E1p20 E1E2
в= E1 + E2 E1E2
По сути, относительное перемещение играет роль внутренней степени свободы, на которую при определенных условиях может происходить передача энергии от центра масс вследствие структурного преобразования. Таким образом, мы приходим к двухуровневой модели материала, первый из которых соответствует движению центра масс, а второй — относительному смещению, переход которого в новое устойчивое состояние равновесия в результате нестационарного воздействия на материал трактуется как изменение его микроструктуры. Любопытно отметить, что первое из уравнений системы (3) может быть получено из закона баланса импульса при условии, что определяющее уравнение материала задается следующим выражением:
'ди Эz
G = Gi +a2 =( El +E2)
-а Эх Эх
(4)
представляющим собой аналог уравнения соотношения Дюамеля-Неймана, обобщающего закон Гука для термоупругости, в котором роль дополнительной степени свободы выполняет температура.
3. Немонотонная определяющая диаграмма
Перед тем как обратиться к исследованию динамической проблемы, рассмотрим квазистатическую задачу о кинематическом растяжении двухкомпонентного стержня, которая представляет самостоятельный интерес с точки зрения возможности не постулировать существование немонотонной определяющей диаграммы, а получить ее в качестве результата решения и определить область параметров, обеспечивающих наличие неустойчивого участка на ней. Рассмотрим закрепленный с одного торца стержень конечной длины I, на противоположном конце которого задана зависимость смещения от времени и0^). Опуская производные по времени в системе уравнений (4), получим
Э2U Э2z Э2z в„ . п ч ту = , ТТ = PK sin(Xz). Эх Эх Эх
(5)
Первое из уравнений системы (5) интегрируется элементарно, и, следовательно, задача сводится к построению решения второго уравнения, которое нужно дополнить соответствующими граничными условиями
4=о = о, и|х=1 = и^). (6)
Нетрудно заметить, что граничных условий (6) не хватает, чтобы установить зависимость напряжения от деформации. Для корректной постановки задачи необходимо еще одно дополнительное граничное условие. Допустим, что напряжение а в сечении х = 0 распределено пропорционально плотности компонентов:
Эи1 Эх
_ Pi 0^
i = 1, 2.
(7)
х=0 р10 + р20
Во всяком случае это предположение не противоречит правилу суммирования напряжений, в соответствии с которым общее напряженное состояния образца равно
сумме напряжений в каждом из компонентов а = а1 + + а2. Если воспользоваться им, то краевая задача (5)-(7) после введения безразмерных переменных
5 = —, w = zX9 P = ъ l
принимает вид
- ö2 sin w = 0,
Ei + E2
u = -
и l''
£c(t)
Uo(t)
_ U 0
l
w
W
w5=0 = 0
9w Э5
Pr (E1 + E2)Xl
wL , =
5=1 E1 E2p0 Pc^^ ( - P),
(8)
15=1
где ö2 = KßXl2, r = ^2Pio - Elp 20, Po=Pio +P 20- После исключения величины P приходим к задаче со смешанными краевыми условиями
W55
w5=0 = 0
9w
э5
ö2sin w = 0,
(9)
w 5=1 +
p0 E1E2 _ P0(Ei + E2)Xl80
2=
■ 5=1 ^ r
Очевидно, что при равенстве скоростей продольной звуковой волны ci = д/Ejp¿о (i = 1, 2) относительного смещения не возникнет, поэтому предположим, что между скоростями имеется небольшое отличие, которое может заключаться как в модулях Юнга Et, так и в плотностях рi. Для определенности положим, что Р10 =Р20, а относительную разницу жесткостей обозначим через 5 = (E2 - E1)/E2. С использованием введенного обозначения краевая задача (9) может быть записана в виде
- Ф2 sin w = 0,
= 0, (10)
55'
w5=0
w 5=1 +
9w
э5
4(1 -5) 2(2 -S)Xl8r
5=1
Для исследования вопроса о существовании параметров, при которых определяющая диаграмма имеет немонотонный характер, предлагается искать приближенное решение задачи (10) в виде ^ = в^, где через в обозначен неизвестный коэффициент. Уравнение, устанавливающее связь между ним и смещением на торце, может быть получено на основе процедуры Галер-кина с использованием метода Гринберга [17]
0
cos 0
sin 0 0
0
1+-
4Х l (2 -5 )(1-5)
+
805_ = 0
- (11) 2(1 -5) ^ ;
Требование наличия экстремумов на определяющей диаграмме приводит к дополнительному уравнению относительно в
Рис. 4. Определяющая диаграмма. ö2 = 4.5 (1), 9.0 (2)
ö2 f 2 I 2ö2
—I -1 I sin 0—Vcos 0-1 = 0,
0 i 02 J 02
(12)
корни которого должны удовлетворять неравенству
Л2
ö2
1 + 5/(4(1 -5))
<0<ö2.
(13)
Численное решение граничной задачи (10) выполняется в программе AUTO 07P, широко используемой при анализе нелинейных проблем [18]. Выражения (12) и (13) позволяют найти значения параметров, обеспечивающих существование петлеобразной зависимости напряжения P от деформации 8 0(t), которая определяется через относительное смещение на торце следующим соотношением:
P (80) = 80
w(1) 5
(14)
2(2 -5)Я/'
устанавливающим связь между микропараметром, в роли которого выступает относительное смещение, и макрохарактеристикой материала. Результаты численного расчета приведены на рис. 4, где представлены две определяющие кривые. Модули Юнга, плотности и размер зерна выбирались в соответствии с характерными значениями для металлов, а значение К, характеризующее связь между компонентами, служит управляемым параметром, обеспечивающим наличие немонотонного участка на определяющей диаграмме. Его оценка будет произведена позднее при анализе динамической задачи. Кривая 1 была построена при следующих параметрах:
E1 = 1-1011 Н • м-2, E2 = 2 • 1011 Н • P1 =P 2 = 4 -103 кг - м-3,
м
К = 4-1010 Н • м-3, I = 5 • 10-3 м, при которых выполняется неравенство (13) и соблюдается единственность решения краевой задачи (10), а следовательно, и требование однозначного определения напряженного состояния по заданной деформации. Противоположная ситуация изображена с помощью кривой 2, которой соответствуют те же самые значения физических характеристик за исключением параметра К, значе-
Рис. 5. Точки возврата. е0 = 3 • 10-3 (1), 4 • 103 (2), 5 • 103 (3)
ние которого было увеличено в 2 раза. В этом случае неоднозначная зависимость относительного смещения на торце приводит к обращению в бесконечность касательного модуля, что показывает невозможность перехода материала в новое состояние при данном виде на-гружения. Полученные результаты приводят к вопросу о нахождении критического значения деформации е0сг, соответствующего точке экстремума на определяющей диаграмме, при достижении которого происходит структурное преобразование. Кривая 1 на рис. 4 была получена на основе численного решения, и поэтому было бы целесообразно попытаться найти чисто аналитический подход для оценки данной величины. При этом предполагается, что параметры материала выбраны таким образом, чтобы обеспечить существование неустойчивого участка на кривой Р(е0). Чтобы установить зависимость от параметров двухкомпонентной модели, будем искать решение нелинейного уравнения (10) в виде суммы
к( $ = м^пЙ) + ), (15)
где определяется решением линеаризованной задачи при близких к нулю относительных смещениях, а второе слагаемое представляет собой малое (w1 << 1) возмущение решения линейного оператора. После подстановки (15) в (10) приходим к уравнению второго порядка с переменными коэффициентами относительно функции W(£) и однородными граничными условиями
^ - 02 к ^ = 0 ^П - ^^. (16)
При структурной перестройке материала происходит смена формы решения, что математически выражается в наличии точек возврата, в которых переменный коэффициент обращается в 0. Для уравнения (16) данное условие записывается в виде
КИп = (17)
До критической точки на определяющей диаграмме (рис. 5), соответствующей началу неустойчивого участка, все возвратные точки лежат справа от 1. На рис. 5 показано, как с ростом деформации точка возврата постепенно оказывается внутри промежутка 0 1, на котором рассматривается нелинейная краевая задача. Уравнение (17) приводит к достаточно простому соотношению для определения критической деформации:
Е0„ . (18)
0сг 25Х1Ш 0 ^ ;
Отметим, что найденное по формуле (18) значение 8 0сг = 4 • 10-3 довольно хорошо согласуется с точкой экстремума на определяющей диаграмме.
Завершая исследование задачи о кинематическом растяжении двухкомпонентного стержня, оценим влияние дополнительной степени свободы на поле деформаций, которое в отсутствие относительного смещения является однородным. Из первого уравнения (5) следует, что в данном случае бессмысленно искать решение как линейную функцию от координаты. Предлагается изменить его форму наиболее простым способом, введя малую поправку и разыскивая его в виде суммы
к© = Д80)£ + ), (19)
где предполагается, что функция к2 мала по сравнению с Л(80). После подстановки (19) в (10) имеем
А(80) = 2^*, (20)
при этом функция к2(£) должна удовлетворять уравнению
к2-028ш( Д + к2) = 0 (21)
с однородными краевыми условиями. Используя предположение о малости к2, уравнение (21) приближенно можно представить как
(22)
к2-02эт( Л£) = 0,
после чего оно элементарно интегрируется:
к?
02 • /^ч
Л£) +
02 4(1 -5)cos Л
(23)
ЛА "" Л (5-2)2 На рис. 6 показано сравнение выражения (23) с численным решением исходной задачи. Из него видно, что асимптотика за исключением области приложения нагрузки качественно верно отражает поведение функции к2(£). Чтобы оценить ее влияние на центр масс, вернемся к первому из уравнений системы (5), которое в ранее введенных безразмерных переменных
И X л и /ч им)
£ = -, к = zX, и = у, 80(t) = принимает вид
Э!и =_5_^. (24)
Э^2 2X1 (2 -5) Э^2 Интегрируя данное уравнение и подставляя в него (23), получим
Рис. 6. Распределение w2(%). е„ = 6 • 103 (1), 7 • 103 (2)
ди = в К cos (А%) Э%=8() 4Ае 0 '
(25)
откуда следует, что отклонение от однородного деформированного состояния имеет периодический характер с частотой, зависящей от отношения размера образца к размеру микроструктуры, параметра 5 и заданной деформации в0. Из (25) нетрудно найти точку %^^ в которой это отклонение максимально. Дифференцируя (25) по % и приравнивая производную нулю, находим, что
% =пк (2-5) % ^ 2Х/5еп '
(26)
т.е. с увеличением деформации максимум смещается вглубь материала. Этот же результат виден из сопоставления с численным решением (рис. 8) при тех же значениях деформации, что и на рис. 7.
Подведем краткие итоги. В данной главе в рамках двухкомпонентной модели изучено влияние микроструктуры на свойства материала при квазистатическом растяжении стержня. Предложенный подход не позволил показать сам процесс перестройки формы относительного смещения решеток, выступающего в роли внутренней степени свободы, но несмотря на это были найдены параметры, при которых оно приводит к наличию петлеобразной определяющей диаграммы. Также был предложен способ оценки критического значения деформации, при котором начинается ниспадающая ветвь, и было проанализировано влияние микроструктуры на макродеформацию, выражающееся в том, что ее распределение в образце после перестройки приобретает гармонический характер со средней составляющей, равной значению приложенной нагрузки. Обращает на себя внимание узкий диапазон значений деформации, приходящийся на ниспадающую ветвь определяющей диаграммы даже при отличающихся в два раза модулях Юнга Е1 и Е2, т.е. при 5 = 0.5, что вряд ли допустимо с физической точки зрения. Более логичным кажется
Рис. 7. Отклонение от средней деформации. в0 = 6 • 10 3 (1), 7 • 10-3 (2)
предположение, что эти характеристики слабо отличаются друг от друга, а 5 является малым параметром (5 << 1). Тогда не выполняется неравенство (13), а значение критической деформации в соответствии с формулой (18) стремится к бесконечности. Таким образом, структурное преобразование материала в данной ситуации не реализуется. Этот вывод остается справедливым для квазистатического типа нагружения. В следующей главе будет показано, что при нестационарном воздействии ситуация меняется и влияние относительного смещения становится достаточно существенным.
4. Динамика структурного преобразования
Вернемся к динамическим уравнениям (4) и рассмотрим их в полубесконечной области 0 < х < Начальные условия считаются нулевыми. Для постановки граничного условия при х = 0, где образец подвергается
Рис. 8. Дисперсионные кривые
ударному воздействию, задаваемому в виде прямоугольного импульса конечной длительности, воспользуемся тем же предположением, которое было использовано при решении квазистатической задачи, именно положим, что напряжение на торце распределено пропорционально плотности компонентов. При переходе к новым переменным и = (р10щ + р20и2)/(р10 + Р20) и г = = щ - и2 граничные условия принимают вид
дх dz_ дх
(
= СТ
х=0
imp
(t)
Pi2o
р2о
Л
= СТ
,(t)
х=0
E1 (P10 +P20)2 Е2р10 - Е1р20 (P10 + р20) Е1Е2
Е2(р10 + р20)
(27)
где ) = ст0(Н^) - Н^ -10)). Через ^ здесь обозначена длительность импульса, а через Н (^ — функция Хевисайда. На противоположной границе будем использовать условия стандартного вида, при которых перемещение при х ^^ полагается равным нулю, что отвечает волнам, уходящим на бесконечность:
ии = o, 4^ = 0 (28)
Уравнения (4) совместно с граничными условиями (27), (28) представляют динамическую задачу о поведении материала с внутренней степенью свободы в виде относительного смещения компонентов, отвечающего за структурные преобразования. После перехода к безразмерным переменным
хю - * ~
X =-, % =fflt, и = иХ, г = гХ,
Ст = -
qXa
v = -
vro
ю* (Е1 + Е2)' КХ
где ю = с1Л/2К Х/ Е1 — частота отсечки, пренебрегая членами О(5) в коэффициентах, начально-краевую задачу можно представить в виде
Э 2U d2Ü
5 Э 2 i
дХ2 dt2 4 ЭХ2
д 2
дХ
dU
дХ
d2i . „ & йэ2U
■—r- = sin z + vz + 5—т-, dt2 dt2
(29)
= СТ
imp
(t),
di дх
= СТ
imp
(t )5.
График дисперсионных кривых для системы (29), описываемых выражением
ю4 -ю2(2к2 +1) + к4 + к2 = 0, (30)
показан на рис. 8. Он состоит из двух ветвей, нижняя из которых отвечает движению центра масс, а верхняя — относительному движению. Динамика последнего описывается вторым из уравнений (29), которое представляет собой нелинейное уравнение Клейна-Гордона, осложненное дополнительным слагаемым, выражающим влияние инерционных сил на относительное смещение. Данное уравнение занимает важное место в
теории нелинейных волн, т.к. оно часто встречается в разнообразных приложениях [19, 20], к которым относятся физика дислокаций, моделирование сейсмических явлений, описание джозефсоновского перехода и т.д. Точное решение нелинейного уравнения нетрудно получить для стационарных волн солитонного типа, здесь же речь идет о задаче Коши. Тогда кажется разумно с самого начала размышлять об использовании того или иного численного подхода, однако в большинстве случаев они не дают возможности выделить множество параметров, обеспечивающих переход системы в искомое состояние, и произвести оценку времени данного перехода. Следовательно, становится сложно заранее указать размеры области численного интегрирования. Также в нелинейных задачах существенную сложность представляет нахождение условий сходимости и устойчивости метода, связанных с выбором величины сеточного шага. В связи с вышеперечисленными вопросами появляется настоятельная потребность в приближенных подходах, позволяющих получить аналитические формулы, с помощью которых можно хотя бы качественно предсказать динамику модели и отделить физику процесса от эффектов, вносимых численной аппроксимацией. Для многих задач в качестве такого подхода может использоваться метод разложения на переменном интервале [21]. Он оказывается плодотворным при построении приближенного аналитического решения нестационарной задачи и может применяться для тестирования численного расчета. Суть этого метода состоит в том, что, аналогично процедуре Галеркина, решение дифференциального уравнения ищется в виде разложения по формам, однако длина промежутка, в отличие от классического подхода, не является постоянной величиной, а представляет собой неизвестную функцию времени. Главная особенность метода переменного интервала заключается в том, что он позволяет свести задачу в частных производных к обыкновенному дифференциальному уравнению, описывающему поведение дискретного элемента сплошной среды. Для обычного упругого стержня таким элементом служит пружинный маятник, тогда как для двухкомпонентной среды он представляет собой систему двух связанных нелинейной силой взаимодействия Q = Q0 sin X (х2 - х1) + ц(х2 - х1) осцилляторов с близкими собственными частотами (рис. 9), где через х1 и х2 обозначены перемещения масс m1 и m2.
В случае их равенства масс динамика системы описывается системой уравнений 5m
х + х = —,
4 (31)
ф + пф + 9 + Ksin ф = 5х.
Величины х = X(х1 +х2)/2 и Ф = Х(х 1 -х2) введены по аналогии с величинами U и z для континуальной
"1 m1
^ww-o
m2
Рис. 9. Связанные осцилляторы
модели. Дифференцирование в (1) производится по безразмерному времени т = ю^, где юг- = (г =
= 1, 2) — парциальные частоты. Малый параметр 5 = = (ю^ - ю2 Vю2 представляет собой их относительную разницу. Коэффициент к = Q0A/при нелинейном слагаемом пропорционален отношению жесткости пружины, характеризуемой модулем Е1, к жесткости связи между компонентами. В дальнейшем принимается, что он равен минимальному значению (к = 1), обеспечивающему невыпуклость потенциальной энергии относительного смещения. Для вывода системы из положения равновесия будем считать начальную скорость центра масс ¥0 отличной от нуля. При малых скоростях взаимодействия между х и ф не происходит, однако начиная с некоторого критического значения картина качественно меняется. Развитие колебаний х во времени, полученное численным интегрированием системы (1) при п = 0.022, показано на рис. 10 пунктирной линией. Характерной особенностью динамического процесса при достаточно большом значении ¥0 является наличие момента времени т*, разделяющего два режима колебаний. На первом участке дискретный элемент ведет себя как линейный осциллятор с близкими собственными частотами. После завершения переходного процесса его динамика также мало отличается от динамики линейной системы, в которой произошла расстройка частот. Получается, что нелинейная сила имеет существенное значение только в узком промежутке времени. Предположим, что работа совершается этой силой непосредственно в точке т*. Тогда нелинейное слагаемое приобретает вид к0$т ф50(т-т*), где через 50(т) обозначена дельта-функция. Параметр к0 подбирается таким образом, чтобы обеспечить переключение центра масс с режима биений на колебания с одной частотой. Такой подход, с одной стороны, сохраняет основные закономерности исходной задачи, а с другой — дает возможность получить аналитическое выражение для перемещения центра масс, которое будет использовано для анализа конти-
нуальной модели. Данную задачу можно трактовать как поиск управления, позволяющего изменить отклик системы на импульсное воздействие. Параметры к0 и т* находятся с помощью преобразования Лапласа fL = = Jo f (т)е-pTdT. При сделанных предположениях без учета вязкого трения решение системы уравнений (31) в пространстве изображений имеет вид
(p2 + 1)V0 - S/ 4 к0 sin Ф(т*)е"рт*
хЪ( p) =
ФЧ p) =
(p2 +1)2 -S2/4 V0S-( p2 +1) к0 sin ф(т,)е"pT*
(32)
(р2 +1)2 -574 '
Для переключения центра масс с режима биений на гармонический процесс, необходимо оставить в системе только одну собственную частоту. Этого можно добиться, потребовав, чтобы один из полюсов в выражении для xL (р) являлся корнем числителя. В результате данное условие приводит к следующим выражениям
для искомых параметров: nk
т* =-
1 -е/ 4
(33)
Kоsin ф(т*) = (-1)"+12К0, (34)
где через k обозначено неизвестное целое положительное число. Оно определяется из энергетического баланса, если приравнять работу сил трения на промежутке от 0 до т* скачку энергии в момент времени т* в задаче с дельта-функцией при выбранных значениях параметра k = 27. Аналитическое решение, полученное при нахождении обратного преобразования Лапласа от хъ(р), показано на рис. 10 сплошной линией. Видно, что оно достаточно хорошо согласуется с численным расчетом.
После того как процесс передачи энергии на внутреннюю степень свободы показан на примере дискретного элемента, вернемся к континуальным уравнениям (29) и воспользуемся методом переменного интервала. Будем разыскивать их решение в виде
и = ÍQn п(2п+1)Х Н (1 - X),
я=1 2
Численное решение Аналитическое решение
100 200 300 400 Рис. 10. Колебания центра масс
_ N П(2п + 1)Х TTÍ1 „. z = Е qn (%)cos—--— h (i - x).
п=1 2/
Умножая на форму и интегрируя в промежутке от 0 до / = где / — длина интервала, получаем следующую систему уравнений относительно функций Q(f) и q(t):
e + 02Q = 028q 2&imp(0
e+0 e =—
q + vq + 02q + 2 J1 (q) = -8Q - -
2&lmp8
(35)
l
где ^=л(2и + 1)/(2/) — частота, отвечающая своей форме. Через J1(q) здесь обозначена функция Бесселя I рода. После исключения из второго уравнения слагаемого б и введения новой переменной т = ^% система (35) приобретает вид
+ =8q - ИЦЙЮ 4 О lm
q+0q+q+ОтГ Ji( q)=8Q. 00
(36)
Эти уравнения практически не отличаются от уравнений дискретной модели (31), и, следовательно, к ним применимы все те же рассуждения, которые были сделаны для того, чтобы получить аналитическую формулу, показывающую уменьшение амплитуды колебаний центра масс за счет динамики внутренней переменной. Из аналогии со связанными осцилляторами сразу же делается вывод, что при слабом импульсе, величина которого ниже критического значения и которому соответствуют малые относительные перемещения, передача энергии от центра масс на другую степень свободы не происходит и волна движется по материалу, не замечая внутренней структуры. В дискретной модели момент переключения т* с режима биений на колебания с одной частотой определялся после замены исходной задачи на линейную, в которой сила взаимодействия действует не в течение малого промежутка времени, а непосредственно в точке т*. При этом мы пренебрегали силой вязкого трения и вместо импульсного воздействия на центр масс задавали начальную скорость. Проделаем эти же преобразования с уравнениями (36). Если внешнее воздействие ст^ (т) в правой части первого из них заменить дельта-функцией с050(т), то, в соответствии с формулой замены переменной 50(ах) = 50(х)/а, они приводятся к виду
e+e=8i _ 2ао8о(т) ^ ^ 4 OI ' 2
q+q+—r Ji(q)8o(T-T*) = 8q. О2
(37)
Параметр т* в дискретной модели определялся из условия гашения одной из собственных частот, которое имеет вид (33). Данное условие здесь сохраняется, тогда как аналог соотношения (34) записывается следующим образом:
12-
- t = 20
— t = 80
- - t = 140
- - t = 200
4-
0- ^
4lt-
V 1
-4-
0 50 100 150 200 x
Рис. 11. Распределение деформации, c0 = 11, t0 = 3
,2oo = (-1)4q) k=, 2
l
0
(38)
где qz = qг (т*). Он позволяет произвести оценку длительности процесса структурных преобразований, которую мы обозначим через %*. Для этого умножим последнее равенство на 1/%*, предполагая существование параметра, аналогичного т* в дискретной модели. С учетом введения новой переменной т = Ол и длины интервала / = % из (38) следует
(-1)
k+1
%* =т*ч-^-. (39)
| sm(oт*/4)cos т*
Подставляя в (39) числовые значения, взятые из дискретной модели 5 = 0.1, к = 27, находим ?* ~ 127. Сделанные выводы подтверждаются численным интегрированием уравнений, выполненным методом конечных разностей при V = 0.022. На рис. 11 показано распределение деформации в = Эи/Эх по длине, демонстрирующее эффект гашения прямоугольного импульса в результате структурного преобразования в материале, которому в данной модели соответствует переход относительного смещения в новое равновесное состояние (рис. 12). Чтобы понять, насколько полученный результат соответствует действительности, необходимо знать частоту отсечки. Примем, что ю ~ 2.5 •109 Гц, исходя из характерного периода колебаний на рис. 4. В таком случае и = I*/ю* ~ 5 • 10-8 с, что можно считать приемлемым результатом при сравнении с экспериментальными данными (рис. 3), учитывая грубость модели и приближенный характер вычислений. Зная частоту отсечки, становится возможным оценить неизвестный параметр К, характеризующий нелинейную силу взаимодействия. Используя соотношение ю ; = с1Л/2КХ/ Е1, получим, что К ~ 4 -1016 н/м3, т.е. данная величина на несколько порядков превышает то значение, которое было использовано при рассмотрении квазистатической
£
_________J t - 20 t - 80 ---t - 140 t - 200 l A
1' * 1 1' ' 1 1' 1 1 1' 1 г 1' ■ 1 |< 1 1 1' 1 1 t 1 -N.fr---ч \ ' 1
0 50 100 150 200
Рис. 12. Относительное смещение, о0 = 11, ?0 = 3
задачи, т.е. для того чтобы получить петлеобразную определяющую диаграмму, мы были вынуждены принять такие значения параметров, при которых материал вряд ли может существовать.
5. Заключение
Данная работа в очередной раз убеждает, что анализ колебательных систем, выполненный на основе динамической и статической постановки задачи, приводит к различным результатам. Более того, переход к динамике может привести к возникновению новых явлений, не происходящих в статическом приближении. Здесь это различие показано на примере материала с изменяющейся структурой при ударно-волновом нагружении. Модель основана на механике двухкомпонентной среды, состоящей из двух кристаллических решеток, связанных нелинейной силой взаимодействия, учитывающей периодическую структуру материала, где в качестве параметра, отвечающего за структурные преобразования, было выбрано их относительное смещение. Сначала была рассмотрена квазистатическая задача о кинематическом растяжении двухкомпонентного стержня с целью определения параметров, обеспечивающих немонотонную зависимость напряжения от деформации. В рамках данной задачи был предложен критерий определения критической деформации, соответствующей началу структурных преобразований, а также выявлен характер влияния относительного смещения на макропараметр, в роли которого выступает перемещение центра масс кристаллических решеток. При этом было отмечено, что параметры, отвечающие наличию неустойчивого участка на определяющей диаграмме, не соответствуют действительности, т.е. переход материала в новое состояние при данном типе нагружения невозможен. Чтобы осуществить структурное превращение, требуется изменить воздействие на материал, сделав его
динамическим. Оно задается в виде прямоугольного импульса, длительность которого много меньше времени пробега волны по образцу, что обуславливает проведение исследования динамических уравнений в полубесконечной области. Основным его результатом стала демонстрация эффекта гашения нестационарной волны вследствие процесса преобразования микроструктуры, длительность которого удалось определить, используя аналогию между уравнениями континуальной среды и ее дискретным аналогом, установленной на основе метода переменного интервала. Несмотря на грубость предложенной модели, в которой, к примеру, не учитываются пластические свойства материала, она позволила приближенно оценить параметр, отвечающий за нелинейно-упругую связь между компонентами. Вопрос об оценке параметра, характеризующего внутреннюю диссипацию материала, остается открытым и, по-видимому, требует более глубокого физического понимания процессов, происходящих на микроуровне.
Работа выполнена в рамках проекта РНФ № 15-1900182.
Литература
1. Duvall G.E., Graham R.A. Phase transitions under shock-wave loading // Rev. Modern Phys. - 1977. - V. 49. - No. 3. - P. 523.
2. Hixson R.S., Boness D.A., Shaner J.W., Moriarty J.A. Acoustic velocities and phase transitions in molybdenum under strong shock compression // Phys. Rev. Lett. - 1989. - V. 62. - No. 6. - P. 637640.
3. Канель Г.И., Фортов B.E., Разоренов С.В. Ударные волны в физике конденсированного состояния // УФН. - 2007. - Т. 177. - № 8. -С. 809-830.
4. Meshcheryakov Yu.I., Divakov A.K., Zhigacheva N.I., Makarevich I.P., Barakhtin B.K. Dynamic structures in shock-loaded copper // Phys. Rev. B. - 2008. - V. 78. - No. 6. - P. 064301.
5. Мещеряков Ю.И., Диваков А.К., Жигачева Н.И., Барахтин Б.К. Пороговые режимы и микромеханизмы динамического деформирования // Mater. Phys. Mech. - 2011. - Т. 11. - С. 23-59.
6. Барахтин Б.К., Мещеряков Ю.И., Савенков Г.Г. Динамические и фрактальные свойства стали СП-28 в условиях высокоскоростного нагружения // ЖТФ. - 1998. - Т. 68. - № 10. - С. 43-49.
7. Мещеряков Ю.И., Савенков Г.Г. Осцилляции фронта пластической волны в условиях высокоскоростного нагружения // ПМТФ. -2001. - Т. 42. - № 6. - С. 117-123.
8. Truskinovsky L., Zanzotto G. Ericksen's bar revisited: Energy wiggles // J. Mech. Phys. Solids. - 1996. - V. 44. - No. 8. - P. 1371-1408.
9. Rosakis P., Knowles J.K. Unstable kinetic relations and the dynamics of solid-solid phase transitions // J. Mech. Phys. Solids. - 1997. -V. 45. - No. 11. - P. 2055-2081.
10. Eremeev V.A., Freidin A.B., Sharipova L.L. Nonuniqueness and stability in problems of equilibrium of elastic two-phase bodies // Dokl. Phys. - 2003. - V. 48. - No. 7. - С. 359-363.
11. Gavrilov S.N., Shishkina E.V. On stretching of a bar capable of undergoing phase transitions // Continuum Mech. Thermodyn. - 2010. -V. 22. - No. 4. - P. 299-316.
12. Нигматулин Р.И. Динамика многофазных сред. Ч. 1. - М.: Наука. Физматлит, 1987. - 464 с.
13. Куропатенко В.Ф. Обмен импульсом и энергией в неравновесных многокомпонентных средах // ПМТФ. - 2005. - Т. 46. - № 1. -С. 7-15.
14. Vilchevskaya E.N., Ivanova E.A., Altenbach H. Description of liquid-gas phase transition in the frame of continuum mechanics // Continuum Mech. Thermodyn. - 2014. - V. 26. - No. 2. - P. 221245.
15. Индейцев Д.А., Наумов В.Н., Семенов Б.Н. Динамические эффекты в материалах со сложной структурой // Изв. РАН. МТТ. - 2007. -№ 5. - С. 17-39.
16. Aero E.L., Bulygin A.N., Pavlov Y.V. The nonlinear theory ofreorgani-zation of structure of superthin crystal layers at intensive loadings // Mater. Phys. Mech. - 2012. - V. 15. - P. 126-134.
17. Свешников А.Г., Боголюбов А.Н., КравцовВ.В. Лекции по математической физике. - М.: Изд-во МГУ, 2004.
18. Doedel E.J., Oldeman B.E. AUT0-07P: Continuation and Bifurcation Software for Ordinary Differential Equations. - Montreal, Canada: Concordia University, 2008.
19. Браун О.М., Кившарь Ю.С. Модель Френкеля-Конторовой. Концепции, методы, приложения. - М.: Физматлит, 2008. - 533 с.
20. Будылина Е.А., Мурачев Е.Г., Гарькина И.А., Данилов М.М. Решения уравнения Клейна-Гордона типа бегущей волны, сглаживающейся на бесконечности // Фундаментальные исследования. -2014. - № 5-5. - С. 1000-1005.
21. Слепян Л.И. Исследование нестационарных деформаций с помощью рядов, определенных на переменном интервале // Изв. АН СССР. Механика. - 1965. - № 4.
Поступила в редакцию 22.07.2017 г.
Сведения об авторах
Морозов Никита Федорович, д.ф.-м.н., акад. РАН, проф., зав. каф. СПбГУ, [email protected]
Индейцев Дмитрий Анатольевич, д.ф.-м.н., чл.-корр. РАН, проф., научн. рук. ИПМаш РАН, [email protected]
Семенов Борис Николаевич, к.ф.-м.н., доц. СПбГУ, [email protected]
Вакуленко Сергей Августович, д.ф.-м.н., проф., внс ИПМаш РАН, [email protected]
Скубов Дмитрий Юльевич, д.ф.-м.н., проф., внс ИПМаш РАН, [email protected]
Лукин Алексей Вячеславович, ассист. СПбПУ, [email protected]
Попов Иван Алексеевич, ассист. СПбПУ, [email protected]
Вавилов Дмитрий Сергеевич, препод. ВКА им. А.Ф. Можайского, [email protected]
