CC BY
85
____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА
Том 157, кн. 1 Физико-математические науки
2015
УДК 519.6:532:533:539.3
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ ВЫСОКОСКОРОСТНОЙ СТРУИ НА ТВЕРДОЕ ТЕЛО
А.А. Аганин, Т.С. Гусева, Н.А. Хисматуллина
Аннотация
Предложена методика численного исследования ударного воздействия струи жидкости на приповерхностный слой упругопластического тела в условиях, близких к кавитационному разрушению. Методика состоит из трех этапов. Сначала в предположении абсолютной твердости тела рассчитывается ударно-волновая динамика в струе и в окружающем газе, определяется закон изменения давления на поверхности тела. Затем этот закон аппроксимируется сплайнами и простыми аналитическими зависимостями. Полученная аппроксимация нагрузки используется для расчета полей напряжений в приповерхностном слое тела и деформаций его поверхности. Для иллюстрации предлагаемого подхода выполнен расчет ударного воздействия высокоскоростной струи воды на твердое тело из никелевого сплава.
Ключевые слова: ударное воздействие струи, ударные волны, области текучести в теле, метод CIP-CUP, адаптивные сетки, метод Годунова.
Введение
Кавитационное разрушение - результат схлопывания кавитационных пузырьков вблизи поверхностей тел в жидкости. Как правило, говорят о его негативных последствиях, например разрушении гидроустройств [1]. Однако кавитационное разрушение имеет и широкое практическое применение, например, в технологиях кавитационной очистки [2]. Если в какой-либо области жидкости давление опускается и удерживается ниже некоторого критического уровня, а затем вновь поднимается, то возникающие кавитационные пузырьки сначала многократно увеличиваются в размерах, а затем сильно сжимаются. Вблизи поверхностей твердых тел пузырьки сжимаются несимметрично. При этом возможно образование на поверхности пузырька высокоскоростной струи, направленной к телу (рис. 1, а). В этом случае кинетическая энергия сходящегося потока жидкости кумулируется в виде кинетической энергии струи. Ударное воздействие таких струй является одной из причин кавитационной эрозии тел. Интенсивность воздействия максимальна, когда пузырек касается стенки, так как в этом случае струя бьет непосредственно по ней.
Несферическому схлопыванию пузырьков около твердой стенки посвящено много работ. В частности, в работе [3] с использованием метода граничных элементов рассчитывалось схлопываниие сфероидального пузырька. Установлено, что наиболее скоростными являются струи, диаметр которых мал по сравнению с размерами пузырька, так что влиянием его боковых стенок в ходе удара можно пренебречь. Таким образом, на первых этапах исследования ударное воздействие струи можно моделировать как воздействие цилиндрического столбика жидкости с полусферическим концом (рис. 1, б).
75
76
А.А. АГАНИН И ДР.
Рис. 1. Коллапс пузырька вблизи тела: (а) — фрагмент в окрестности пузырька, (б) — фрагмент в окрестности конца струи. Белый цвет — пузырек, серый — жидкость
Рис. 2. (а) Схема удара струи по телу. Черными стрелками обозначена ударная волна в жидкости, белыми стрелками — волны в теле; (б) удар капли диаметром 10 мм со скоростью 110 м/с по твердому телу, эксперимент [4]
При ударе кумулятивной струи по поверхности тела в самой струе возникает ударная волна, распространяющаяся вверх, а в приповерхностном слое тела формируются волны продольных и сдвиговых напряжений, расходящиеся от места воздействия (рис. 2, а). Наличие ударных волн в струе подтверждается близкими по условиям экспериментами [4] со скоростным ударом капли по твердой стенке (рис. 2, б).
Несмотря на большое число исследований (обзоры в работах [5, 6]), изучение напряжений и деформаций в теле при высокоскоростном ударе кумулятивной струи остается актуальным. Трудности экспериментального и теоретического исследований обусловлены малостью пространственно-временных масштабов ударного воздействия. Оно реализуется в области радиусом в десятые доли радиуса струи и в промежутке времени, за который ударная волна в теле пробегает эту область.
В настоящей работе предлагается методика численного исследования ударного воздействия струи на твердое тело в условиях кавитационного разрушения. Жидкость в струе и окружающий газ (содержимое пузырька) описываются уравнениями газовой динамики. Применяется численная методика на основе предложенной в [7] модификации метода CIP-CUP на динамически адаптивных soroban-сетках. Напряженно-деформированное состояние тела описывается уравнениями динамики идеального упругопластического полупространства, которые решаются методом С.К. Годунова. Расчет ударного воздействия струи на стенку разбивается на три стадии. Сначала в предположении абсолютной твердости тела и без учета вязкости жидкости (это соответствует достаточно большим пузырькам и телам, по свойствам близким к металлическим) рассчитывается ударно-волновая динамика в струе и окружающем газе. Результатом первой стадии является дискретная динамическая нагрузка на поверхности тела в ходе удара струи по стенке. На второй стадии эта нагрузка аппроксимируется сплайнами и аналитическими зависимостями. На третьей стадии полученная аппроксимация нагрузки используется для расчета деформаций поверхности тела и полей напряжения в приповерхностном слое.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
77
1. Постановка задачи
Рассматривается ударное воздействие высокоскоростной струи жидкости на плоскую поверхность твердого тела. Струя представляет собой цилиндрический столбик с полусферическим концом, с осью симметрии z, направленной перпендикулярно поверхности тела. За начало отсчета времени t = 0 принимается момент, когда конец струи касается поверхности тела в точке с координатами r = 0, z = 0, где r - радиальная координата, отсчитываемая от оси струи. Струя окружена газом, давление в струе и газе распределено однородно. В начальный момент времени имеет место конфигурация, показанная на рис. 1, б.
2. Основные положения методики расчета
2.1. Этап 1. Рассчитывается динамика жидкости в струе и окружающего газа в процессе удара струи по телу. Тело предполагается абсолютно твердым, вязкость жидкости не учитывается. Применяется подход с неявным «вылавливанием» межфазной границы. Вводится функция-идентификатор р такая, что в области жидкости р = 1, в области газа р = 0, а в малой окрестности их контактной границы р непрерывна и монотонно меняется от 0 до 1. Контактная граница при этом трактуется как переходная область, занятая смешанной средой. Значения функции-идентификатора сохраняются вдоль траекторий частиц среды, ее изменение описывается уравнением переноса
д + u -Ур = °. (1)
Для сквозного расчета динамики жидкости и газа используется метод CIP-CUP (CIP-Combined UP) [7]. Он базируется на уравнениях движения сжимаемой среды относительно неконсервативных переменных (плотность р, скорость u, давление р), что и позволяет осуществлять сквозной расчет без явного выделения межфазной границы [8]. Динамика контактирующих жидкости и газа без учета эффектов вязкости и теплопроводности описывается уравнениями
др
dt
+ u ■Ур
р У ■ u,
ди
д
+ u ■ yu
Ур
р
др
dt
+ u ■ Ур
рС1 У^ u. (2)
Здесь Cs = рСв\ + (1-р)Cs2, Csi = W+ЩГр, Csi - скорость звука, Г1,Б1 -
константы уравнения состояния Тэта для жидкости, Г2 = y, B2 =0 - для газа (y - показатель адиабаты).
Систему (1), (2) можно записать в виде
Of
- + u У = G, (3)
где f = (р, р, u,р), G = (0, -рУ ■ u, -р-1Ур, -рС13У ■ u).
Численное решение системы (3) разбивается на конвективную и неконвективную стадии
f * f n
Atn
+ u У
fП+1 f*
---;---- = G.
Atn
0,
(4)
(5)
Для расчета группы уравнений переноса (4) применяется полулагранжев метод CIP (Constrained Interpolation Profile) [9]. Решение уравнения переноса
дf dt
+ u ■Уf
0
78
А.А. АГАНИН И ДР.
У
о
4
8 9 10
6_______________7
3 4 5
|------•----•---------
3
2
/=1
х
2
j=1
Рис. 3. Двумерная soroban-сетка, j - номер направляющей линии, i - номер узла
в узле x в момент tn+1 = tn + Atn можно приближенно записать в виде f (x — uAtn), где аргумент представляет собой отправную точку - положение «лагранжевой частицы», которая за время Atn переместилась в узел x. Значение переносимой величины f в отправной точке определяется CIP-интерполяцией с использованием известных узловых значений f и ее пространственных производных fx,fy ,fXy в момент tn .
В рамках метода UP (Unified Procedure) [9] неконвективная часть (5) сводится к следующему уравнению относительно давлени pn+1
V
Vpn+1
р*
pn+1 — p* у . u*
p*C*S2Atn2 + Atn ’
которое решается итерационным методом. Затем пересчитываются значения скорости un+1 и плотности pn+1. При решении задач с интенсивными ударными волнами метод CIP-CUP применяется в сочетании с искусственной вязкостью.
Расчеты проводятся с использованием адаптивной soroban-сетки [7]. В двумерном случае это набор узлов, расположенных на параллельных направляющих линиях (рис. 3). На каждом временном шаге строится новая сетка, адаптирующаяся к решению и никак не связанная со старой сеткой. При этом взаимное расположение и количество как направляющих, так и узлов на них может изменяться.
При построении сетки сначала рассчитывается монитор-функция
M (x,y,t)
^ 1 + а
дх
2
+
ду
цЧ + в
d2f д 2f
д2 х д2у
(6)
Здесь f - какой-либо параметр решения, а, в - положительные константы. Направляющие и узлы сетки сгущаются в зонах с большими значениями мониторфункции и разрежаются там, где ее значения малы. При построении конечно-разностных аппроксимаций производных на soroban-сетке требуемые значения расчетных параметров в узлах шаблона схемы определяются CIP-интерполяциией [7].
Разработанные на основе описанной методики алгоритмы и программы были протестированы на задачах в одномерной и двумерной постановке. Численные решения сравнивались с известными аналитическими решениями и численными решениями других авторов [10-13].
2.2. Этап 2. Непосредственное использование дискретной динамической нагрузки, рассчитанной на первом этапе, затруднено из-за различия расчетных сеток первого и третьего этапов. Кроме того, полученная на первом этапе динамическая нагрузка на стенке на фоне выраженного общего характера изменения имеет высокочастотные осцилляции, предположительно, численной природы. Поэтому на этом этапе выявляется общий характер изменения нагрузки, для чего применяется аппроксимация сплайнами и простыми аналитическими зависимостями.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
79
Рис. 4. Аппроксимация распределения давления на поверхности тела в момент времени t. Основные точки - Ai, i = 1, 2, 3,4, вспомогательные точки - Bj , j = 1,..., 8
Характерный радиальный профиль нагрузки на поверхности тела при ударе по нему струи жидкости с полусферическим концом согласно [14-16] имеет вид, представленный на рис. 4. Величина нагрузки в области контакта струи и тела монотонно возрастает от центра (r = 0) к периферии. На границе области контакта r = rmax (t) она достигает максимума и далее резко падает до уровня давления в окружающем газе. Скорость роста нагрузки в более широкой центральной части области контакта остается относительно небольшой, тогда как в более узкой части на периферии она значительно возрастает. Подобный характер распределения давления на стенке сохраняется до тех пор, пока величина давления по порядку близка к максимальной.
Аппроксимация профилей нагрузки, полученных на первом этапе, осуществляется следующим образом. Радиальный профиль давления на стенке разбивается на три отрезка, ограниченных четырьмя основными точками - на рис. 4 им соответствуют точки А1, Л2, Аз, A4. Абсциссами этих точек являются соответственно ri = 0, r2 = 1/3rmax (t), r3 = 2/3rmax (t) и r4 = rmax (t). Сначала строится аппроксимация временной зависимости rmax (t), а также временных зависимостей давления в точках Ai, i = 1, 2, 3, 4. Для построения аппроксимации этих зависимостей применяется сплайн-интерполяция с использованием 10-15 опорных точек, соответствующих дискретным временным зависимостям, рассчитанным на первом этапе. Опорные точки выбираются так, чтобы в результате сплайн-интерполяции получились монотонно изменяющиеся функции, повторяющие общий характер изменения аппроксимируемых кривых. При этом для зависимостей p(ri,t), i = 1, 2, 3,4, в опорной точке t = 0 давление полагается равным давлению гидроудара (“water hammer pressure”) [4]
= Vj pi Di ps Ds
^ pl Dl + ps Ds ’
где Vj - скорость струи, pi , ps и Di , Ds - плотность и скорость распространения ударных волн в жидкости и теле соответственно. Для rmax (t) первая опорная точка есть rmax (0) = 0 .
Аппроксимация зависимости p(r4, t) усложнена рядом особенностей. Среди расчетных кривых первого этапа именно эта кривая имеет наибольшие высокочастотные осцилляции. В то же время ее общий характер отличается от p(ri, t), i = 1, 2, 3, наличием выраженного максимума. Действительно, согласно [14, 15] такой максимум пика давления на периферии области контакта струи и стенки, равный примерно 3p0, достигается в момент Tj [4]:
Rv
j
2 D2’
Dl = CS1 +
Г1 + 1
Vn.
j
j
2
80
А.А. АГАНИН И ДР.
При построении аппроксимации p(r4,t) рассчитанная временная зависимость с осцилляциями заменяется на гладкую кривую с острым пиком в точке максимума с сохранением общего характера изменения вне этого максимума. Множество точек этой кривой, которые служат опорными точками для сплайн-интерполяции, дополняется точкой (Tj, 3po). Интерполяция строится таким образом, чтобы максимум аппроксимирующей зависимости находился в этой точке.
С использованием построенных аппроксимаций rmax(t) и p(rj,t), i = 1, 2, 3,4, в каждый момент t аппроксимируется радиальный профиль давления на стенке. Из-за наличия пика давления на периферии области контакта интерполяция сплайнами не подходит. С учетом этого на отрезках (Ai, A2) и (A2, A3) применяется линейная интерполяция, а сплайн-интерполяция используется лишь на отрезке (A3, A4). Для реализации сплайн-интерполяции вводится ряд вспомогательных точек Bi, B2,..., Bn (рис. 4, n = 8). Их число и положение в каждый момент времени t определяется в процессе построения ломаной линии, соединяющей точки A3 и A4 . Углы наклона составляющих этой ломаной растут в геометрической прогрессии, коэффициент которой зависит от соотношения значений давления в точках A3 и A4 .
2.3. Этап 3. Рассчитывается динамика приповерхностного слоя тела под действием аппроксимированной динамической нагрузки, полученной на втором этапе. Линейные размеры струи в условиях кавитационного разрушения обычно намного меньше размеров тела, поэтому тело считается упругопластическим изотропным полупространством. Деформации и перемещения предполагаются малыми. С учетом этого для описания динамики тела используются следующие уравнения:
ди д(Srr — P) + dSrz +
Pdt dr
Srr — S(
dz
dP к (du dv и dt V dr + dz + r
dSrr 2 f 2 du dv и
dt 3 P\ dr dz r
dSvv 2g / du dv 2u\
dt 3 1 dr + dz r )
dv = dSrz + d(Szz — P) + Srz P dt dr dz r ,
dSzz 2 (du 2 dv u A
dt 3P\dr dz + r) ,
dSrz (du dv A
~dT = p\d~z + dr),
Здесь u, v - компоненты скорости по осям r и z соответственно; Srr, Szz, Sфф, Srz - компоненты девиатора S тензора напряжений а, записанные в цилиндрической системе координат; P - всестороннее (гидростатическое) давление; S = а + + Pg; g - метрический тензор; А = р(с2 — 2c|) и g = pc\ - параметры Ламе; K = А + 2/3g - коэффициент объемного расширения; ci и С2 - скорости распространения продольных и сдвиговых возмущений.
В пластических зонах выполняется условие текучести Мизеса аj = Y0, аj -интенсивность напряжений,
1
72
aj = (az — ar )2 + (az — аф )2 + (ar — аф)2 + бт^
где az = Szz — P, ar = Srr — P, аф = Sфф — P, Tr-z = Srz - компоненты тензора
напряжений а в цилиндрических координатах.
Если условие текучести нарушается, то есть если aj в некоторой точке расчетной области превышает предел текучести, то девиатор тензора напряжений корректируется по формуле [17]
S ___*0 S
Scorr
а
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
81
Область нагружения / я,. гя
1 т 1
тело 1 1 1 1 1 1 1 1
1 1 1 1
Рис. 5. Расчетная область
На поверхности тела z = 0 полагается
= -p(r,t), Trz = 0, (8)
где p(r, t) - аппроксимация давления, полученная на втором этапе по результатам вычислений первого этапа. В начальный момент времени t = 0 полагается
Cz = -Pg, = Сф = -\(\ + 2p)-1pg, Trz =0, u = v = 0, (9)
где pg - давление в окружающей струю газовой среде.
На бесконечном удалении от места приложения нагрузки параметры тела остаются невозмущенными.
В расчетах бесконечная область, занимаемая телом, 0 < r < ж, —ж < z < 0 заменяется конечной цилиндрической областью 0 < r < rs, hs < z < 0 (рис. 5). Внешние границы r = rs и z = hs являются «искусственными». Они выбираются так, чтобы отраженные от них волны в интересующем промежутке времени не сильно изменяли напряжения в области влияния струи. Удаленность искусственных границ зависит от принимаемых условий. В настоящей работе на искусственных границах ставятся неотражающие условия [18], что позволяет размеры расчетной области принять равными rs = —hs = 2.2R, где R - радиус струи [19].
Для численного решения задачи (7)-(9) применяется классический метод С.К. Годунова [20]. Шаг по времени определяется из условия устойчивости [20]
S Дг • Az At = — -Т----
ci Ar + Az
при S = 0.95.
Более подробное описание методики расчета динамики тела и иллюстрацию ее применения можно найти в [21, 22].
3. Ударное воздействие высокоскоростной струи воды на твердое тело из никелевого сплава
Для иллюстрации работоспособности предлагаемого подхода рассматривается задача об ударе водяной струи по телу из никелевого сплава. Все параметры в струе и окружающем газе распределены однородно. Радиус струи R = 100 мкм, ее скорость Vj = 300 м/с, газ в окрестности конца струи движется с той же скоростью. В области струи имеем: р = 103 кг/м3, Г = 7.15, Bi = 3047 бар, р = 1. В области газа: р = 1 кг/м3, Г2 = 1.4, B2 = 0, р = 0. Давление p в струе и газе всюду равно 1 бар. Никелевый сплав имеет плотность р = 8000 кг/м3, модуль упругости E = 196000 МПа, коэффициент Пуассона v = 0.3, его предел текучести при одноосном растяжении Yo = 125 МПа. В этих условиях давление гидроудара в случае абсолютно жесткой стенки po = 6.3 кбар, Tj = 3.4 нс.
82
А.А. АГАНИН И ДР.
Рис. 6. Осевое сечение области взаимодействия струи со стенкой. Левая колонка - изолинии давления и межфазная граница (р = 0.5): 1 - ударная волна, 2 - граница струи, 3 - пристеночная струйка. Правая колонка - вид сетки: а, б - момент t/rj ~ 0.3, в, г -t/rj ~ 1, д, е - t/rj ~ 3
3.1. Динамика жидкости в струе и нагрузка на поверхности тела.
На рис. 6 приведены результаты расчетов в три последовательных момента ударного воздействия струи на поверхность тела (стенку). Процесс протекает очень быстро в небольшой окрестности конца струи. Приведенные поля получены при расчете половины осевого сечения размером [0.4R, 0.2R]. Минимальный шаг сетки 2.5 • 10-4R, максимальный шаг - 1.6 • 10-2R. При построении мониторфункции (6) использовались пространственные производные давления и функции-идентификатора, а = 1,^ = 0.1. Переменный шаг по времени определялся из соотношения
Atn
kcrt min
1
(|un | + cn)
11
~— +—~— Ax Ay
-1
i
Хотя условие Куранта не имеет отношения к устойчивости схемы CIP-CUP, удобно использовать число Куранта kcrt для определения соотношения временного и пространственного шагов. В задачах с ударными волнами оптимальным является
0.1 < kcrt < 1, при kcrt > 1 может возникать чрезмерное размывание фронта ударной волны [13]. В данном случае kcrt = 0.95. Коэффициент искусственной вязкости в жидкости - 0.6 с удвоенной линейной составляющей, в газе - 0.9.
Процесс воздействия струи на тело во многом определяется ударной волной, распространяющейся вверх по струе. При этом сначала (до момента t ~ Tj, рис. 6, в), боковая граница этой волны совпадает с границей расширяющейся области контакта струи и стенки. В области сжатой жидкости, отделенной криволинейным фронтом ударной волны от невозмущенной части струи, формируется умеренное движение жидкости к кромке контакта струи и стенки. Скорость расширения области контакта уменьшается из-за увеличения угла наклона между границей струи и стенкой. В момент Tj скорость расширения становится равной скорости распространения ударной волны. В результате этого ударная волна отрывается от стенки и перемещается вверх по границе струи с формированием волны
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
83
Рис. 7. Профили давления на стенке в моменты t/Tj ~ 0.1, 0.3, 1, 2, 3, 4, 6
разрежения (рис. 6, д). Это приводит к тому, что в пристеночной области на периферии области контакта возникает тонкая струйка цилиндрически растекающейся жидкости.
На рис. 7 показано изменение распределения давления на стенке в процессе ударного воздействия струи. Для экономии расчет проводился сначала на грубой сетке, а затем уменьшающаяся окрестность момента t = Tj последовательно уточнялась на ряде более мелких сеток. Кривые 1-4 получены на сетке с минимальным шагом 4.7 • 10_5R, кривые 5, 6 - с минимальным шагом 1.25 • 10_4R, кривая 7 - с минимальным шагом 3• 10_4R. В начале ударного воздействия давление во всей области контакта достигает давления гидроудара р0 = 6.3 кбар. Затем распределение давления приобретает характер, описанный в разд. 2.2. Давление в центре области контакта начинает снижаться из-за упомянутого выше оттока сжатой жидкости к кромке контакта. А в узкой кольцевой области на периферии области контакта появляется выраженный максимум. Наибольшие значения давления возникают здесь в течение небольшого отрезка времени в окрестности момента Tj (рис. 7, кривая
3). При t > Tj пик давления на периферии контакта медленно спадает (рис. 7, кривые 4-7).
На рис. 8 приведена зависимость от времени максимального давления на стенке (при r = rmax (t)), а также зависимости давления в центре области контакта r = 0 и на расстоянии r = 1/3 rmax (t), r = 2/3 rmax (t). Часть кривой до t/Tj = 3 получена на сетке с минимальным шагом 4.7 • 10_5R, далее до t/Tj =4 - 1.25 • 10_4R, до t/Tj = 7 - 3 • 10_4R, до t/Tj = 12 - 1.25 • 10_3R и далее - 2.5 • 10_3R. Расчет на более мелкой сетке проводился до момента совпадения со значением, полученным на более грубой сетке. В ходе ударного воздействия струи различие давлений в центре области контакта и на расстояниях вплоть до r = 2/3 rmax (t) относительно невелико. На периферии области контакта на начальной стадии воздействия при t < Tj максимальное давление быстро нарастает. Наибольшие значения достигаются в интервале 0.8 < t/Tj < 2 и примерно в 3 (и в 4 в малой окрестности t/Tj = 1) раза превышают начальное давление гидроудара ро. Затем максимальное давление начинает снижаться, причем до уровня давления гидроудара р0 оно снижается только при t & 6Tj. При этом расхождение зависимостей максимального давления на периферии области контакта и давления в центре уменьшается и почти исчезает в момент t/Tj & 20, то есть распределение давления на стенке в области воздействия становится близким к однородному. При t/Tj & 35 давление в этой области падает до 1 кбар.
Кривая 1 на рис. 7 и кривая максимального давления при r = rmax(t) на рис. 8 характеризуются наличием высокочастотных осцилляций. Они имеют численное происхождение, обусловленное, в частности, тем, что криволинейная граница струи
84
А.А. АГАНИН И ДР.
о-------*-----■-----■-----■-------------
О 10 20 t/xj 30
Рис. 8. Временные зависимости максимального давления на стенке, достигаемого на расстоянии r = rmax (t) от центра, давления на стенке в центре области контакта r = 0 и на расстоянии r =1/3 rmax (t) , r = 2/3 Г max (t)
в начальный момент задается на декартовой сетке отрезками прямых линий. Во втором случае, кроме того, разброс значений максимального давления в численных данных обусловлен тем, что в каждый момент времени узлы адаптивной сетки попадают не точно в текущий пик давления, а лишь в различные точки в его окрестности, в которой давление стремительно снижается по мере удаления от максимума. Расчеты показывают, что при измельчении сетки немонотонность уменьшается, а пики давления на периферии области контакта, наоборот, становятся более острыми. Однако сходимости по значениям давления на пике в окрестности t = Tj достичь не удалось. Можно предположить, что дальнейшее уточнение этих пиков в целом не изменит картину динамической нагрузки, поскольку они реализуются очень кратковременно в точечной области на стенке, непрерывно удаляющейся от центра области контакта вследствие натекания струи.
3.2. Аппроксимация динамической нагрузки. Рис. 9 иллюстрирует аппроксимацию расчетной временной зависимости максимального давления на поверхности тела p(rmax(t),t). Аппроксимация получена с применением сплайнинтерполяции. На рисунке представлен лишь начальный промежуток времени, где функция p(rmax(t),t) меняется наиболее быстро, а расчетная кривая имеет наибольшие по амплитуде осцилляции. На основе известных аналитических оценок [4, 14] в качестве экстремального давления жидкости на поверхности тела при построении аппроксимационной зависимости выбрано давление, в три раза превышающее давление гидроудара ро. Видно, что там, где осцилляции расчетной кривой относительно небольшие (при t/Tj < 0.5 и t/Tj > 3), аппроксимирующая функция довольно хорошо согласуется с численным решением. На отрезке [0.5, 3] приближение строится так, чтобы оно проходило через выбранные опорные точки и представляло собой плавную монотонную кривую как в области ее роста (t < Tj), так и в области ее падения (t > Tj).
Аппроксимации расчетной временной зависимости p(rmax(t),t) при t/Tj > 6.5, а также зависимостей р(0, t), p(1/3rmax(t),t), p(2/3rmax(t),t) и координаты rmax(t) получаются относительно хорошими. Этого удается достичь без дополнительных усилий благодаря относительной гладкости аппроксимируемых кривых.
На рис. 10 показаны аппроксимации пространственного распределения давления для тех же моментов времени, что и на рис. 7. Видно, что приближения вполне удовлетворительно соответствуют кривым численного решения.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
85
Рис. 9. Аппроксимация временной зависимости максимального давления на стенке в интервале 0 < t < 6т (штриховая линия). Сплошная линия - результаты расчетов. Опорные точки для сплайн-интерполяции отмечены кружками
Рис. 10. Аппроксимация профилей давления на стенке в те же моменты времени, что и на рис. 7
3.3. Динамика приповерхностного слоя тела. Рис. 11 характеризует изменение поля интенсивности напряжений <ii в теле в девять последовательных моментов времени процесса ударного воздействия струи. Видно, что сразу после начала воздействия интенсивность напряжений <ii достигает своего предельного значения Yq . К моменту ti возмущения затрагивают лишь небольшую часть приповерхностного слоя тела в окрестности оси симметрии струи, а зона, где <ii = Yq , представляет собой куполообразную область, радиальные размеры которой чуть превышают текущий радиус области нагружения, а осевые размеры - в 1.5 раза меньше его. Это обусловлено тем, что в начале удара скорость расширения области нагружения больше, чем скорость распространения возмущений в теле.
С расширением области нагружения расширяется и возмущенная область в теле, а вместе с ней и зона с предельным значением интенсивности напряжений. К моменту t2 различие размеров этой области в радиальном и осевом направлениях уменьшилось до 1.2, радиальные размеры стали заметно больше радиуса области нагружения. Это связано с уменьшением скорости расширения области нагружения, которая, будучи в начале удара бесконечно большой, с течением времени быстро уменьшается.
В момент t3 границы как возмущенной области, так и зоны предельных значений интенсивности напряжений в теле имеют форму, близкую к полусферической. При этом радиус зоны, где = Yq , уже в полтора раза больше радиуса области нагружения поверхности тела.
К моменту t6 в окрестности оси симметрии у поверхности тела интенсивность напряжений опустилась ниже предельного значения Yq . Зона же с предельными
86
А.А. АГАНИН И ДР.
Рис. 11. Динамика приповерхностного слоя тела после удара струи (закрашена область, где Gi = Yo). Моменты ti -tg соответствуют t/rj ~ 0.6, 1, 4, 8, 9, 10, 11, 16, 38
значениями интенсивности напряжений распалась на две части, большую и маленькую. Внешняя граница большой части по форме подобна той, что в моменты £3 -1§. Увеличились лишь ее линейные размеры. Маленькая часть зоны с ai = Yq располагается вблизи нагруженной области поверхности тела у края этой области.
К моменту t8 размеры зоны с предельными значениями интенсивности напряжений уменьшились, а сама она отошла от поверхности тела. Вместе с тем в этот момент в дополнение к этой области у самого края нагруженной области поверхности также имеется небольшая зона с предельными значениями ai, имеющая форму тонкого кольца. К моменту t9 зона с предельными значениями интенсивности напряжений полностью исчезает.
На рис. 12 изображены профили поверхности тела в окрестности оси симметрии в те же девять моментов времени, что и на рис. 11. Для наглядности масштаб по z увеличен в 100 раз. В момент времени t\ радиус области нагружения поверхности тела еще мал, давление на периферии этой области - около 10 кбар. Ямка
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
87
Рис. 12. Динамика поверхности тела после удара струи (область тела закрашена темным цветом). Реальный масштаб по оси z на два порядка меньше, чем по r
на поверхности едва заметна, ее радиус совпадает с радиусом области нагружения. В момент t2 величина нагрузки на периферии области ее приложения достигает своего максимума. Глубина ямки на поверхности тела к этому моменту увеличилась, ее радиус по-прежнему совпадает с радиусом области нагружения. На краю ямки начинает формироваться выступ. К моменту t% ямка углубилась, высота выступа увеличилась и достигла своего максимума. Изменилась форма ямки. Теперь ее дно стало близким к плоскому, а ее боковая поверхность - близкой к конусообразной, сужающейся вниз. Ширина дна ямки примерно равна размерам области нагружения поверхности тела в момент t2, когда величина нагрузки на краю области нагружения достигает своего экстремального значения. В промежутке t^ -1^ давление на поверхности тела падает, глубина и радиус ямки продолжают расти. При этом форма ямки и размер ее дна остаются прежними (боковая поверхность несколько выполаживается), основание выступа увеличивается, его высота уменьшается. В момент t8 форма ямки прежняя, а ее глубина достигает максимума. К этому моменту радиус ямки в ее верхней части оказывается меньше радиуса области нагружения. Поверхность тела под краем нагрузки немного деформирована (продавлена вниз), выступ несколько опустился. К моменту tg и позднее глубина ямки немного уменьшается, тогда как ее форма не изменяется.
88
А.А. АГАНИН И ДР.
Таким образом, в рассматриваемом случае ударного воздействия струи на поверхности тела возникает осесимметричная микроямка с небольшим выступом на внешней границе. Дно ямки близко к плоскому, боковая поверхность близка к конической с расширением вверх. Максимальная глубина ямки 0.00524Д, ее радиус примерно равен 0.525Д. Максимальная высота выступа - 0.00064Д, что на порядок меньше максимальной глубины ямки.
Заключение
Представлены математическая модель и основные положения методики расчета ударного воздействия струи на приповерхностный слой упругопластического тела. Условия воздействия близки к условиям кавитационного разрушения, когда струя образуется при схлопывании касающегося тела кавитационного пузырька. Методика расчета состоит из трех этапов. Сначала в предположении абсолютной твердости тела и без учета вязкости жидкости (это соответствует достаточно большим пузырькам и телам, по свойствам близким к металлическим) рассчитывается ударно-волновая динамика жидкости в струе и окружающего газа, а также закон изменения давления на поверхности тела. Затем этот закон аппроксимируется сплайнами и простыми аналитическими зависимостями. Далее эта аппроксимация используется для расчета полей напряжений в приповерхностном слое тела и деформаций его поверхности.
Для иллюстрация представленной методики выполнен расчет ударного воздействия высокоскоростной струи воды на тело из никелевого сплава. Струя цилиндрическая с полусферическим концом, направлена ортогонально к поверхности тела. Скорость струи V = 300 м/с, радиус Д = 100 мкм. Показано, что в результате удара такой струи на поверхности тела возникает осесимметричная микроямка (радиусом ~ 0.5Д и глубиной ~ 0.005Д) с небольшим выступом на внешней границе (высотой ~ 0.0006Д).
Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 14-01-97004_ р_поволжье_а) и программы президиума РАН
(«Динамика деформируемого тела и откольные разрушения при ударном воздействии кавитационных пузырьков», координаторы Н.Ф. Морозов, Н.Л. Горячева).
Summary
A.A. Aganin, T.S. Guseva, N.A. Khismatullina. Numerical Simulation of High-Speed Jet Impact on a Rigid Body.
A technique is presented for numerical studying of the liquid jet impact on the surface layer of an elastic-plastic body under conditions close to cavitation damage. The technique comprises three stages. First, assuming that the body is perfectly rigid, the shock-wave dynamics in the jet and surrounding gas is calculated and the law of pressure variation on the body surface is derived. Then, the derived law is approximated by splines and simple analytic expressions. The obtained approximation is used to determine the fields of stresses in the surface layer of the body and deformations of its surface. The proposed approach is illustrated by computing the impact of high-speed water jet on the rigid body out of nickel alloy.
Keywords: jet impact, shock waves, yield areas in body, CIP-CUP method, adaptive grids, Godunov method.
Литература
1. Brennen C.E. Hydrodynamics of pumps. - N. Y.: Oxford Univ. Press, 1994. - 293 p.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УДАРНОГО ВОЗДЕЙСТВИЯ...
89
2. Kieser B, Phillion R., Smith S., McCartney T. The application of industrial scale ultrasonic cleaning to heat exchangers // Proc. Int. Conf. on Heat Exchanger Fouling and Cleaning / Eds. M.R. Malayeri, H. Muller-Steinhagen, A.P. Watkinson. - 2011. -P. 336-338.
3. Воинов О.В., Воинов В.В. О схеме захлопывания кавитационного пузырька около стенки и образования кумулятивной струи // Докл. АН СССР. - 1976. - Т. 227, № 1. -С. 63-66.
4. Field J.E., Dear J.P., Ogren J.E. The effects of target compliance on liquid drop impact // J. Appl. Phys. - 1989. - V. 65, No 2 - P. 533-540.
5. Field J.E., Camus J.-J., Tinguely M., Obreschkow D, Farhat M. Cavitation in impacted drops and jets and the effect on erosion damage thresholds // Wear. - 2012. - V. 290-291. -P. 154-160.
6. Lauterborn W., Kurz T. Physics of bubble oscillations // Rep. Prog. Phys. - 2010. -V. 73, No 10. - Art. 106501, P. 1-88.
7. Takizawa K., Yabe T., Tsugawa Y., Tezduyar T.E., Mizoe H. Computation of free-surface flows and fluid-object interactions with the CIP method based on adaptive meshless Soroban grids// Comput. Mech. - 2007. - V. 40, No 1. - P. 167-183.
8. Abgrall R., Karni S. Computations of compressible multifluids // J. Comp. Phys. -2001. - V. 169, No 2. - P. 594-623.
9. Yabe T., Xiao F., Utsumi T. The constrained interpolation profile method for multiphase analysis // J. Comput. Phys. - 2001. - V. 169, No 2. - P. 556-593.
10. Аганин А.А., Гусева Т.С. Методика расчета волн в жидкости и газе методом CIP-CUP с применением динамически-адаптивных Soroban-сеток // Вестн. Башкир. ун-та. - 2014. - Т. 19, № 2. - С. 368-380.
11. Аганин А.А., Гусева Т.С. Численное моделирование динамики неоднородных сжимаемых сред на основе метода CIP-CUP на адаптивных Soroban-сетках // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2014. - Т. 156, кн. 2. - С. 55-72.
12. Аганин А.А., Гусева Т.С. Численное моделирование контактного взаимодействия сжимаемых сред на эйлеровых сетках // Учен. зап. Казан. ун-та. Сер. физ.-матем. науки. - 2012. - Т. 154, кн. 4. - С. 74-99.
13. Аганин А.А., Гусева Т.С. Расчет контактного взаимодействия сжимаемых сред без явного выделения межфазных границ // Вестн. Башкир. ун-та. - 2013. - Т. 18, № 3. -С. 646-661.
14. Heymann F.J. High-speed impact between a liquid drop and a solid surface // J. Appl. Phys. - 1969. - V. 40, No 13. - P. 5113-5122.
15. Lesser M.B. Analytic solutions of liquid-drop impact problems // Proc. R. Soc. Lond. A. -1981. - V. 377. - P. 289-308.
16. Haller K.K., Ventikos Y.,Poulikakos D. and Monkewitz P. Computational study of highspeed liquid droplet impact // J. Appl. Phys. - 2002. - V. 92, No 5. - P. 2821-2828.
17. Уилкинс М.Л. Расчет упруго-пластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.
18. Ильгамов М.А., Гильманов А.Н. Неотражающие условия на границах расчетной области. - М.: ФИЗМАТЛИТ, 2003. - 240 с.
19. Аганин А.А., Малахов В.Г., Халитова Т.Ф., Хисматуллина Н.А. Расчет силового воздействия кавитационного пузырька на упругое тело // Вестн. Тат. гос. гуманит.-пед. ун-та. - 2010. - № 4 (22). - С. 6-13.
90
А.А. АГАНИН И ДР.
20. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. - М.: Наука, 1976. - 400 с.
21. Аганин А.А., Ильгамов М.А., Хисматуллина Н.А. Упруго-пластические деформации в теле при ударном воздействии кавитационного пузырька // Учен. зап. Казан. ун-та. Сер. физ.-матем. науки. - 2013. - Т. 155, кн. 2. - С. 131-143.
22. Аганин А.А., Хисматуллина Н.А. Ударное воздействие струи жидкости на упругопластическое тело // Учен. зап. Казан. ун-та. Сер. физ.-матем. науки - 2014. - Т. 156, кн. 2. - С. 72-86.
Поступила в редакцию 14.12.14
Аганин Александр Алексеевич - доктор физико-математических наук, профессор, заведующий лабораторией, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.
E-mail: aganin@kfti.knc.ru
Гусева Татьяна Сергеевна - кандидат физико-математических наук, старший научный сотрудник, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.
E-mail: ts.guseva@mail.ru
Хисматуллина Наиля Абдулхаевна - кандидат физико-математических наук, доцент, старший научный сотрудник, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.
E-mail: nailya_hism@mail.ru
