Исследование реологии льда на основе численного моделирования медленного удара

Полный текст

Лед играет центральную роль в Арктическом регионе, который привлекает большой интерес исследователей. Каждый тип природного льда имеет отличные от других типов механические свойства [1], определяемые его кристаллической структурой [2, 3]. Модель изотропной линейной теории упругости выбрана базовой моделью для поликристаллического льда [4]. Несмотря на это, подбор подходящих параметров для данной модели осложняется влиянием на лед различных факторов, таких как температура [5, 6], соленость [7, 8] и пористость [9].

Другой особенностью льда является изменение его поведения при различных способах деформирования. При прохождении акустической волны малой амплитуды лед деформируется упруго. Затем при приложении существенных нагрузок либо осуществляется переход в пластический режим, либо происходит хрупкое разрушение [10]. Некоторые модели позволяют учесть оба эффекта [11, 12]. Вязкость льда рассматривается в [13], а в [14] изучается его ползучесть. Однако для медленного столкновения вопрос о выборе подходящей модели не до конца решен.

В связи с данным фактом в сообщении исследуется поведение льда при низкоскоростном ударе шаровым индентором. Целью ставится численное моделирование лабораторного эксперимента, проведенного в лаборатории Института проблем механики им. А.Ю. Ишлинского Российской академии наук. Изучаются реологические модели упругопластичности с критерием текучести фон Мизеса [15] и Мизеса–Шлейхера [16] и упругости с упругопластическим включением в зоне удара. Решение определяющей системы уравнений производится сеточно-характеристическим методом [17, 18] на структурированных сетках.

ЛАБОРАТОРНЫЙ ЭКСПЕРИМЕНТ

Лабораторный эксперимент проводился с ледяными дисками в холодильной камере, в которой поддерживалась постоянная температура –10 ˚С. В качестве индентора применен стальной шар, внутри которого жестко крепился пьезоэлектрический акселерометр. Скорость шара контролировалась высотой его поднятия, в начале соударения она равнялась 0.56 м/с. Ледяной диск опирался на массивную плиту с возможностью скольжения. Для регистрации фронта прошедшей волны на тыльной поверхности ледяного диска в точке на линии удара также крепился акселерометр. Геометрические характеристики эксперимента, а также экспериментальные графики показаны на рис. 1.

 

Рис. 1. Сверху – общий вид расчетной области и параметры сетки, 2D. Снизу слева – экспериментальные графики, синяя кривая – с приемника в шаре, фиолетовая – с приемника на нижней поверхности льда. Снизу справа – макет льда в модели упругости с упругопластическим включением в форме полукруга заданного радиуса.

 

РЕОЛОГИЧЕСКИЕ МОДЕЛИ

Для моделирования эксперимента рассматривались несколько реологических моделей. В качестве основной определяющей системы уравнений использовалась гиперболическая система уравнений изотропной линейной теории упругости [4], в которой неизвестными выступают напряжения σ и скорости v. Данная модель описывает продольную и поперечную волны со скоростями с_p, c_s. В итоге характеристики c_p, c_s и плотность ρ можно использовать как параметры упругой модели. В шаре и плите они задавались равными: c_p = 5700 м/с, с_s = 3100 м/с, ρ = 7800 кг/м³, во льду: с_p = 3600 м/с, с_s = 1942 м/с, ρ = 917 кг/м³.

Для решения системы определяющих уравнений использовался сеточно-характеристический метод [17, 18]. Согласно методу, производится расщепление по физическим процессам и переход к инвариантам Римана. Таким образом, изначальная система сводится к системе независимых уравнений переноса, для решения которых использовалась сеточно-характеристическая схема 3-го порядка, которая монотоноизировалась с помощью сеточно-характеристического критерия монотонности [19, 20].

Для учета пластического поведения льда модель упругости модифицировалась. Использовалась версия модели Прандтля–Рейсса [15], в которой производилась коррекция девиатора тензора напряжения, , при нарушении критериев текучести: . В качестве критериев были выбраны критерий фон Мизеса [15]: 0.5s_ijs_ij – k² > 0, k – предел упругости, а также критерий Мизеса–Шлейхера [16], в котором присутствует зависимость предела упругости от давления : . На иллюстрациях для данной модели нами используется сокращение УП (упругопластичность).

Последняя рассматриваемая модель льда заключается в выделении из упругой среды области пластичности в зоне удара. В данной работе эта область имеет форму полукруга заданного постоянного радиуса с центром в центре ледяного диска как на рис. 1 снизу справа. В ячейках, входящих в границы полукруга, использовалась вышеописанная модель упругопластичности, в остальных – модель упругости. Данная модель получила название “модель упругости с УП включением”. В расчетах варьировались предел упругости k (в дальнейшем единицы измерения, Па, будут опущены), параметры в критерии Мизеса–Шлейхера k₀ (Па), a, радиус упругопластического включения r (мм).

ПОСТАНОВКА ЗАДАЧИ

Для численных расчетов в двумерном случае расчетная область разбивалась на сегменты согласно рис. 1: 1, 2 соответствуют шару, 3 – льду, 4 – плите. В каждом сегменте строилась структурированная сетка, число ячеек вдоль горизонтальной N_x и вертикальной N_y осей указаны на рисунке, шаг по времени равнялся 5 ∙ 10^–8 c. Сетки под номерами 2 формировались вращением сетки 2*. Между областями 1–2 ставилось контактное условие полного слипания, между 2–3 и 3–4 – условие проскальзывания в области контакта. На боковых и нижней границах подставки использовалось условие поглощения, на поверхности шара, льда, выступающей верхней поверхности подставки была установлена свободная граница. В качестве начального условия в сетках шара задавалась скорость соударения 0.56 м/с, сетки сдвигались с помощью коррекции по Лагранжу. Узлы сеток шара и льда, находящиеся на расстоянии 0.05 мм друг от друга, считались контактирующими. Изначальное расстояние между ними также равнялось 0.05 мм.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

На рис. 2 представлены волновые картины в расчетах по модели упругости с упругопластическим включением (а) и модели упругопластичности с условием Мизеса–Шлейхера (б). В первом случае образуются структуры, которые, по-видимому, соответствуют образованию трещин [21] и наблюдаются в других численных экспериментах [22]. На рис. 2а также видна характерная лунка, образованная в этих структурах. Лед начинает приподниматься по краям зоны контакта с шаром в момент времени 0.375 мс после начала расчета, до тех пор, пока лунка полностью не формируется к 0.75 мс. Данные явления не наблюдаются на рис. 2б. Таким образом проявляется недостаток изображенной модели.

 

Рис. 2. Волновые картины в расчетах по модели упругости с упругопластическим включением r = 7.5, k = 3 ∙ 105 (а); по модели упругопластичности с условием Мизеса–Шлейхера в центральной области льда, k0 = 3 ∙ 105, a = 0.5 (б).

 

Сравнение с лабораторным экспериментом производилось на основе скорости шара, рассчитанной по силе F из экспериментального графика и массе шара m по формуле

 

Дополнительно применялось преобразование

 

для сведения конечной скорости к нулю. Таким образом, получаются графики, которые должны соответствовать модулю вертикальной проекции скорости. Для проведения оценок в качестве первого критерия для подбора параметров было выбрано время, когда ускорение шара меняет свой знак и скорость шара достигает минимума, t_v=0 = 0.7 мс.

Дополнительной величиной для сравнения была выбрана глубина осадки шара, координата его нижней точки. Затем получившаяся вертикальная проекция скорости интегрировалась для получения координаты:

 

Для интегрирования использовалась формула Симпсона. На основе полученной величины был выработан второй критерий отбора параметров – максимальная глубина осадки: x_max = 0.2253 мм.

В результате расчетов были построены кривые, изображенные на рис. 3: модуль вертикальной проекции тензора напряжений σ_yy, глубина осадки x и модуль скорости v_y ≈ v от времени в нижней точке шара. Можно заметить качественную разницу между кривой скорости из эксперимента и кривыми, полученными в результате расчетов, что может быть связано с несовершенством выбранных моделей. Рисунок 3 также позволяет определить основные особенности рассматриваемых моделей. Для модели упругости характерна наивысшая амплитуда σ_yy, наименьшее значение x_max, а также наименьшее время t_v=0. При использовании во льду модели упругопластичности с критерием текучести фон Мизеса с достаточно большим пределом упругости k решение будет совпадать с моделью упругости. При уменьшении k значительно падает амплитуда σ_yy, увеличивается x_max и время t_v=0, а кривые скорости сдвигаются вправо. Для данной модели удается подобрать параметр k = 7.5 ∙ 10⁵ так, чтобы получить значение критериев близко к экспериментальным (ярко-голубая кривая на рис. 3 во втором ряду: t_v=0 = 0.7065, x_max = 0.2281).

 

Рис. 3. Результаты расчетов. Ряды имеют одинаковую легенду. 1-й столбец – модуль вертикальной проекции тензора напряжений от времени в нижней точке шара, 2-й – координата, 3-й – модуль скорости.

 

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

Следующей рассмотренной моделью стала модель упругопластичности с критерием Мизеса–Шлейхера. В данном случае изменение параметра k₀ сходно с изменением k. Наибольший интерес представляет изменение a. Увеличение a приводит к увеличению суммарного предела упругости, однако при этом отсутствует стремление результатов к модели упругости, что можно наблюдать при a > 1 на верхнем ряду на рис. 3. К тому же на графиках модуля вертикальной проекции тензора напряжений начинают образовываться осцилляции, сходные с эффектами на деформационных кривых в других экспериментах [23]. Таким образом проявляется упругопластическое поведение мишени.

Благодаря введению зависимости предела упругости от давления открывается возможность для получения бесконечного множества параметров, позволяющих воспроизвести значения t_v=0, x_max близко к экспериментальным. При этом появляется ограничение, связанное с тем, что при k₀ ≥ 7.5 ∙ 10⁵ невозможно подобрать подходящий параметр a, так как увеличение a сдвигает кривые скорости влево, дальше от нужного значения t_v=0.

В последней рассматриваемой модели, модели упругости с упругопластическим включением с критерием фон Мизеса, центральным параметром является радиус включения r. Согласно среднему ряду на рис. 3, при нулевом значении r производится расчет по модели упругости. Затем при увеличении r понижается амплитуда σ_yy, причем это происходит неравномерно: сильнее в начале столкновения. Также увеличивается x_max и время t_v=0. Для k = 3 ∙ 10⁵ для того, чтобы подобрать значение t_v=0, близкое к экспериментальному, необходимо искать необходимое значение радиуса в диапазоне от 7.5 до 8 мм. Однако при этом совпадение по координате x не будет наблюдаться.

В случае варьирования k для разных r можно получить бесконечное множество подходящих параметров согласно выбранным критериям. Несмотря на это, согласно нижнему ряду на рис. 3, для r ≤ 7 не удается подобрать k. Это связано с тем, что для того чтобы подобрать значение t_v=0, близкое к экспериментальному, в данном случае необходимо уменьшать k и сдвигать кривые скорости вправо. При этом возникает предел, при котором предел упругости становится настолько малым, что критерий фон Мизеса начинает выполняться во всех точках включения. Таким образом, дальнейшее уменьшение k не будет оказывать значимого влияния на результаты.

Примерно при таких же радиусах r ≈ 7 меняется поведение кривых | σ_yy | (t). При больших радиусах изменение кривых при увеличении k гораздо более равномерное: амплитуда на всех этапах соударения растет примерно одинаково. При малых радиусах на поздних этапах столкновения напряжения постепенно увеличиваются. Также кривая координаты при k = 3 ∙ 10⁵ и r = 7 демонстрирует несовершенство выбранного критерия подбора параметров по x_max. В данном случае x_max = 0.2266 близко к нужному значению в 0.2253. Однако при этом из-за того, что t_v=0 оказывается значительно меньше экспериментального, график x располагается левее. При больших радиусах удается подобрать k так, чтобы график сдвигался вправо, что позволяет более точно воспроизвести вид экспериментальной кривой даже при x_max > 0.2253. Из-за этого одним из направлений дальнейшей работы можно указать выработку новых критериев подбора параметров, основанных, например, на мере отклонения результатов расчетов от эксперимента.

В дополнение для r ∊ [8, 10] удается воспроизвести значения t_v=0 и x_max близко к экспериментальным. Для r = 8 (пурпурная кривая на нижнем ряду на рис. 3) при k = 6 ∙ 10⁵ они равны 0.6805 и 0.2302. Для r = 9 (желтая кривая) при k = 6.5 ∙ 10⁵ – 0.7045 и 0.234 и при k = 7 ∙ 10⁵ (ярко-голубая кривая) – 0.6785 и 0.2259. Для r = 10 (серо-голубая кривая) при k = 7 ∙ 10⁵ – 0.7035, 0.2308.

Для количественного воспроизведения вышеописанных трендов на рис. 4 были построены графики t_v=0 и x_max в зависимости от параметров рассматриваемых моделей. Можно отметить высокую корреляцию между выбранными критериями. Для рассмотренных моделей подобраны аппроксимирующие кривые, найдены параметры, при которых выполняется: t_v=0 = 0.7 и x_max = 0.2253. Для начала зависимость критериев от k в модели упругопластичности с критерием фон Мизеса (темно-синяя кривая в третьем ряду) имеет ниспадающий тренд. Данная кривая располагается правее кривых для моделей с критерием Мизеса–Шлейхера и с упругопластическим включением.

 

Рис. 4. Максимальная глубина осадки (слева) и момент времени, когда модуль скорости минимален (справа), в зависимости от параметров рассматриваемых моделей. В каждом ряду сплошные линии соответствуют кривым с одинаковыми параметрами (легенда по центру).

 

Для модели упругопластичности с критерием Мизеса–Шлейхера графики критериев от a экспоненциально затухают, зависимость от k₀ имеет форму ниспадающей гиперболы. Однако при малых k₀ и больших a > 1 качество аппроксимации падает, появляется новое ограничение: уменьшение k₀ не приводит к росту значений критериев. С другой стороны, при k₀ ≥ 0.8 и k₀ = 0.3 визуально воспроизводится невозможность подбора подходящего a.

В случае модели упругости с упругопластическим включением графики t_v=0 и x_max от k имеют немного более пологий ниспадающий тренд. Для r ≤ 10 появляется возможность для линейной аппроксимации кривых. При уменьшении r наблюдается уменьшение коэффициента наклона прямых и их амплитуды. Снова наблюдается разделение поведения кривых по радиусу: оказывается, что для r ≤ 8 критерии становятся почти постоянными при всех значениях k. При устремлении r к 0 проявляется стремление к модели упругости. Данный факт подтверждается графиками t_v=0 и x_max от r. Аппроксимирующие кривые в данном случае хорошо описываются функцией арктангенса. Однако при малых k поведение кривой похоже на более простую степенную зависимость.

Результаты рассчитанных прогнозов по аппроксимирующим кривым представлены в табл. 1. Заметно, что полученные предсказанные значения довольно близки друг к другу. В результате построенные аппроксимации позволяют сужать область допустимых параметров моделей. Однако более полную информацию можно получить при рассмотрении многомерных зависимостей критериев от параметров, что можно указать в качестве цели для дальнейшей работы.

 

Таблица 1. Рассчитанные значения параметров по аппроксимирующим кривым на рис. 4

№ ряда	Фиксированный параметр	Прогноз по tv=0	Прогноз по xmax
1	k0 = 5 ∙ 105	a = 0.17	a = 0.175
2	a = 0.01	k0 = 7.4 ∙ 105	k0 = 7.44 ∙ 105
2	a = 0.1	k0 = 6.14 ∙ 105	k0 = 6.26 ∙ 105
2	a = 0.2	k0 = 4.75 ∙ 105	k0 = 4.8 ∙ 105
4	k = 3 ∙ 105	r = 7.79	r = 6.9
4	k = 5 ∙ 105	r = 10.02	r = 9.16
4	k = 8 ∙ 105	r = 24.9	r = 23.24
4	k = 8 ∙ 104	r = 7.79	r = 6.71

 

ЗАКЛЮЧЕНИЕ

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

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Исследование выполнено при финансовой поддержке гранта Российского научного фонда (проект № 23-21-00384).

КОНФЛИКТ ИНТЕРЕСОВ

Авторы данной работы заявляют, что у них нет конфликта интересов.

Об авторах

И. Б. Петров

Московский физико-технический институт (национальный исследовательский университет)

Автор, ответственный за переписку.
Email: petrov@mipt.ru

Член-корреспондент РАН

Россия, Долгопрудный, Московская обл.

Е. К. Гусева

Московский физико-технический институт (национальный исследовательский университет); Институт проблем механики им. А.Ю. Ишлинского Российской академии наук

Email: guseva.ek@phystech.edu
Россия, Долгопрудный, Московская обл.; Москва

В. И. Голубев

Московский физико-технический институт (национальный исследовательский университет)

Email: golubev.vi@mipt.ru
Россия, Долгопрудный, Московская обл.

В. П. Епифанов

Институт проблем механики им. А.Ю. Ишлинского Российской академии наук

Email: evp@ipmnet.ru
Россия, Москва

Список литературы

Дополнительные файлы