On the wave dynamics in damaged shells interacting with the volume of the cavitating liquid

Full Text

Введение. Кавитационные процессы, происходящие в жидкости, превращают её в двухфазную среду с дисперсной фазой в виде пузырьков, наполненных газом и/или паром. Наличие в жидкости подобной фазы даже малой концентрации меняет характер распространения в ней ударных волн и волновые процессы в целом, приводит к возникновению дополнительных локальных механизмов нагружения и разрушения взаимодействующих с нею деформируемых сред (конструкций). Структура деформируемых сред также неоднородна, включает в себя как исходные, так и появляющиеся в процессе распространения ударных волн микродефекты, рост и объединение которых приводит к катастрофическому разрушению. Эти микродефекты, наблюдаемые микроскопией на мезоуровне структуры в виде пор и трещин, по существу также можно рассматривать в виде дисперсной фазы в нелинейно деформируемой среде подобно двухфазной жидкости с очень вязкой несущей фазой [1]. Анализ взаимодействия рассматриваемых сред оказывается достаточно сложной задачей из-за необходимости учёта реального состояния и взаимодействия фаз, фазовых превращений, истории волнового нагружения и особенностей течения. Разработка методов решения подобных задач до сих пор остаётся актуальной проблемой современной механики. В инженерных приложениях представляют интерес два основных аспекта этой проблемы. Традиционно это защита от кавитации сооружений и объектов техники, взаимодействующих с жидкостью, которая приводит к ухудшению их работы, генерации нежелательных шумов, вибрациям и разрушению. Другим сравнительно новым аспектом проблемы, обусловленным разработкой современных технологий, является не просто защита, а управление кавитацией для уменьшения сопротивления движущихся в жидкости объектов и ультразвуковой очистки поверхностей, для интенсификации процессов в химии и создания мощных акустических импульсных лазеров (сазеров) или, например, удаления камней в почках и др. [2-4]. Особенности распространения ударных волн во взаимодействующих средах определяются многими факторами - теплофизическими и механическими свойствами самих сред, содержанием и распределением пузырьков (кавитационного облака), вязкостью жидкости и наличием поверхностного натяжения, так же как гидродинамическими параметрами их течения, включая характеристики потока, уровни турбулентности и др. Происходящие при этом внутрифазовые и межфазовые явления и процессы являются чрезвычайно сложными для изучения, поскольку весьма локализованы, имеют различную физическую природу и различные пространственно-временные масштабы. Известные экспериментальные исследования, выполненные в этом направлении, отличаются невысокой достоверностью и ограниченной применимостью полученных результатов. Во многих случаях прямые измерения параметров физического поля оказываются пока недоступными [3, 5]. В математическом моделировании необходимо располагать соответствующими моделями сред, позволяющими описать локальные внутрифазовые и межфазовые процессы, возникающие в результате эволюции отдельной дисперсной частицы под действием ударной волны, и связать их с волновыми процессами в сжимаемых средах. Из-за отсутствия последовательного экспериментального изучения всех этих явлений и процессов, сопровождающих течения, в настоящее время нет надлежащего обоснования предлагаемых для их описания теоретических моделей. На практике к традиционно используемым в механике сплошных сред подходам обычно добавляется аппроксимация многофазной среды её составляющими в однофазном состоянии с известными физическими, термодинамическими и механическими параметрами. Поведение жидкости и пузырьков в ней, деформируемых сред и микроповреждений описывается далее раздельно на основе решения соответствующих краевых задач, а взаимодействие между ними - уравнениями связи, реально отражающими процессы, проис367 П е т у ш к о в В. А. ходящие на границах раздела фаз и сред, и обеспечивающих требуемую точность моделирования [3, 6, 7-10] и др. Возникающие здесь краевые задачи с математической точки зрения оказываются нелинейными, относятся к классу так называемых задач Стефана с подвижной границей и требуют для своего решения разработки соответствующих методов. В статье предлагается единый подход для моделирования динамики взаимодействующих между собой двухфазных жидкостей с парогазовыми пузырьками и деформируемой среды с микроповреждениями при распространении в них ударных волн. Разрушение деформируемой среды рассматривается при этом как эволюция микроповреждений - пор сферической формы, принимаемых по аналогии с кавитирующей жидкостью в виде пузырьков газа, объединение которых в процессе вязкопластического течения ведёт к образованию макротрещины. Решение сформулированной нелинейной краевой задачи динамики многофазной среды, включающей в себя уравнения взаимодействия фаз и фазовых превращений, строится на основе методов конечных разностей (МКР) и конечных элементов (МКЭ). Последний метод в силу своей дискретной природы, кроме того, хорошо сочетается с используемым в работе подходом к схематизации сред. Согласно этому подходу, в каждой макроскопической точке жидкой среды выделяется ячейка - конечный элемент, содержащая пробную дисперсную частицу и приходящую на неё часть несущей фазы. Распределение основных микропараметров внутри ячейки описывается уравнениями соответствующих микропроцессов с граничными условиями на их поверхностях. Применительно к деформируемой среде используются соотношения теории вязкопластического течения и континуальной механики повреждений. Выбор типа и размеров ячеек в МКЭ и/или разностных схем МКР определяется требуемой точностью аппроксимации геометрии задачи и её решения в рассматриваемом пространстве. Решение получаемых полудискретных уравнений МКЭ на временном слое осуществляется затем на основе оптимальных (рациональных) конечно-разностных схем и метода расщепления [10-13]. В качестве модельной в статье рассматривается задача об ударном нагружении оболочечной конструкции, составленной из соосных оболочек, пространство между которыми заполнено сплошной сжимаемой жидкостью. Для неё известны результаты, полученные численно в предположении упругого деформирования оболочек [14]. Эта же задача рассматривается далее и для более сложного случая ударного взаимодействия нелинейно деформируемых повреждаемых оболочек с кавитирующей жидкостью. 1. Краевая задача динамики гетерогенных сред. В R3 рассматривается движение многофазной среды, заполняющей открытую область Ω, ограниченную поверхностью S при распространении в ней ударных волн, где Ω ⊂ R3 , Ω = D ∪ S и S = R3 \Ω. Пусть D = DS ∪ DF , где DS и DF - подобласти, занимаемые нелинейно деформируемой средой с повреждениями и взаимодействующей с ней жидкостью с парогазовыми пузырьками соответственно. Положение частиц жидкости и деформируемой среды на их контактной по368 О волновой динамике повреждаемых оболочек . . . верхности Sc (xk , t) = DS ∩ DF в каждый момент времени t ∈ Dt = (0, τ ) определяется в декартовой системе координат xk , k = 1, 2, 3. В начальный момент времени t = 0 двухфазную среду в области 2 F DiF D = i=1 принимаем состоящей из сжимаемой жидкости - несущей фазы (i = 1) c объёмной концентрацией α1 с истинной плотностью ρ01 , и пузырьков газа или пара - дисперсной фазы (i = 2) с объёмной концентрацией α22 1, содержащих вещество с истинной плотностью ρ02 . Полагая, что в любом объёме жидкости находится n пузырьков сферической формы радиуса R, объёмная концентрация α2 , обычно не превышающая 5 %, определяется выражением1 4 α2 = πR3 n. 3 В дальнейшем истинную и приведённую плотности ρi вещества каждой i фазы будем считать связанными соотношением ρ0i = ρi /αi , следовательно, α1 + α2 = 1, и тем самым рассматривать лишь относительные доли объёма, занимаемого веществом каждой фазы, не отслеживая взаимное положение фаз. Необходимость учёта движения пузырьков относительно жидкости, приводящего к изменению их формы и разрушению, возникает для достаточно вязких жидкостей при резком изменении скорости движения частиц среды. Её можно оценить по значениям чисел Вебера и Бонда [6]. Нелинейно деформируемую среду с рассеянными в ней микроповреждениями, занимающую область 2 DS = DiS , i=1 будем рассматривать также в виде двухфазной сжимаемой жидкости с очень вязкой несущей фазой (i = 1) с истинной плотностью ρ03 и вязкостью γ3 . При этом повреждения или дисперсная фаза (i = 2) принимаются в виде достаточно мелких пузырьков газа сферической формы - микродефектов (пор, трещин, включений) с объёмной концентрацией α3 , обычно не превышающей 1-3 % в исходном до нагружения состоянии. Такие условия в жидкости реализуются при очень малых значениях чисел Рейнольдса p3 R Rep = 1, 0 γ3 /ρ3 ρ03 где p3 - давление несущей фазы среды. Принимая поведение газа в пузырьке изотермическим, можно пренебречь инерционными эффектами от движения стенок пузырьков и движением пузырьков относительно несущей фазы. 1 При большей концентрации α2 пузырьковая жидкость превращается сначала в пену, а затем - в газокапельную среду с изменением структуры течения. 369 П е т у ш к о в В. А. Разница между давлениями в несущей и дисперсной фазе уравновешивается в этом случае вязкими силами в несущей фазе. Математическая модель для описания гидродинамических течений рассматриваемых жидкости и деформируемой среды может быть получена как частный случай из уравнений сохранения для гетерогенных сред, учитывающих поведение каждой i фазы и взаимодействия между ними [6, 7, 12]: (αi ρi ),t + div(αi ρi vi ) = Gi , (αi ρi vi ),t + div(αi ρi vi vi + σi ) = mi - fi , (αi Ei ),t + div(αi Ei vi + αi σi vi + qi vi ) + αi,t σi + Qi,t = Li , (xk , t) ∈ Ω × Dt , Ω ∈ R3 , i = 1, 2, . . . (1) В дополнение к используемым обозначениям здесь vi = vik - вектор скорости макроскопического движения; σi = σikl - тензор напряжений, причём σikl = = σikl - ρi ∆vik ∆vil , где ∆vik - вектор пульсационной скорости; q - искусственная вязкость, используемая для сглаживания разрывов в профилях изменения рассматриваемых полевых функций; Qi - тепловая составляющая энергии Ei , обусловленной тепловыми потоками и внутренним тепловыделением; Gi , mi , fi , Li - члены межфазового взаимодействия, k, l = 1, 2, 3. Здесь и далее используются соглашения, принятые в тензорном исчислении. Уравнения (1) включают в себя осреднённые функции и их производные по координатам и времени в пространстве Ω × Dt . Осреднение мгновенных значений параметров процесса выполняется по микрообъему δVi , ограниченному поверхностью δSi с нормалью ni , занятому веществом i фазы. Для замыкания системы уравнений необходимо определить скорости передачи массы, импульса и энергии между фазами [6, 9]. Энергетическая связь взаимодействующих сред определяется членами Li = λij (θi - θj ) + µij (vik - vjk )2 + ηi Gi , где λij , µij и ηi - коэффициенты связи. Вследствие высоких скоростей деформирования взаимодействующих сред можно пренебречь теплообменом, полагая Qi = 0 и соответствующие термодинамические процессы адиабатическими. Влияние температуры θi на параметры сред будем учитывать только в рассматриваемой материальной точке области Ω. Воздействие внешних сил учитывается членом fi . Тензор напряжений σikl для вещества i фазы в гетерогенной среде принимается симметричным, а составляющие шарового тензора напряжений σikk принимаются пропорциональными её объёмной концентрации: σikl = -αi pi δ kl + Sikl . (2) Здесь pi - давление, Sikl - тензор сдвиговых напряжений, или девиатор напряжений. Обычно давления pi в фазах принимаются одинаковыми, что является условием их совместного деформирования, хотя вязкости фаз в выражении для тензора Sikl могут быть разными. При наличии поверхностного натяжения давления различных фаз не могут быть равны, и это учитывается величинами σi в приведённых уравнениях. 370 О волновой динамике повреждаемых оболочек . . . Условия на поверхности раздела дисперсной и жидкой фаз принимаются в виде равенства давлений с учётом вязкостей и поверхностного натяжения: (3) pi (ρi , θi ) = p, и аналогично для температур θi . Кроме того, вводится равенство интенсивностей фазовых преобразований, если они имеют место на поверхностях раздела фаз [9]. На границе контактного разрыва двухфазной жидкости с деформируемой средой, обозначим её Sc (xk , t) с индексом 3, эти условия дополняются следующими: 1 - σ 3 )n = 0, (σjk jk k (vj1 - vj3 ) = 0, x1k = x3k , (xk , t) ∈ Sc × Dt , (4) где nk и vj - соответственно вектор нормали и вектор скоростей на поверхности Sc (xk , t). Для однозначной разрешимости системы уравнений (1) необходимо ввести краевые условия, определяемые конкретной задачей, и уравнения состояния для рассматриваемых сред, обычно принимаемых в виде (5) pi = pi (ρi , ei ), где ei - внутренняя энергия. Состояние сжимаемой несущей фазы (i = 1) при распространении в ней ударных волн может быть описано, например, уравнением ударной адиабаты p1 = A ρ1 ρ1 -1 +B -1 ρ01 ρ01 2 + Cρ1 e1 , (6а) где ρ01 , ρ1 - начальная и текущая плотность жидкости; A, B, C - параметры, значения которых получены из наилучшего приближения к экспериментальной кривой. Для воды уравнение состояния может быть принято также в форме Тэта ρ1 κ p1 =A + B, (6б) p01 ρ01 справедливой до уровней давления порядка 105 МПа, где κ - экспериментально определяемый показатель политропы. При достижении давлением критического значения p∗1 происходит нарушение сплошности жидкости, которое учитывается дополнительными соотношениями: 1 A ρρ01 - 1 , p1 > -p∗1 , p1 = (6в) -p∗1 , p1 p∗1 . Для дисперсной фазы в жидкости (i = 2) ситуация с выбором уравнения состояния не столь однозначна. Как известно, основным фактором, определяющим эволюцию пузырька в жидкости, является скорость движения его стенки. В волне разрежения пузырёк может увеличиваться в размерах более чем на порядок от исходного, обычно составляющего 10-6 м. При сжатии пузырька в ударной волне в жидкости около его поверхности возникают 371 П е т у ш к о в В. А. большие градиенты давления и тонкие тепловые пограничные слои [9, 15]. От градиента давления зависит скорость сжатия пузырька, а от тепловых пограничных слоев - интенсивность образования в нём пара, масса которого так же сильно влияет на изменение радиуса. Гидродинамическое приближение (5), справедливое для описания большей части эволюции пузырька, оказывается неприменимым для моделирования заключительной стадии его разрушения вследствие экстремально высоких давлений и температур, возникающих внутри пузырька. Их величины могут достигать тысячи градусов и, соответственно, атмосфер, вызывая свечение пузырьков в жидкости, или сонолюминесценцию [3]. Поведение парогазовой смеси в этом случае определяется процессами, происходящими на межмолекулярных расстояниях с масштабами времени столкновения молекул. Поскольку описание заключительной стадии разрушения пузырьков выходит за рамки данной статьи, ограничимся следующим приближением: pb (x, t) = pg0 [V (t)/V0 ]k + pv - σC ◦ , (7) полагая содержимое пузырька состоящим из неконденсируемого газа и пара, которое следует политропному закону сжатия с константой k, и параметры его состояния постоянными всюду внутри пузырька. Эти предположения хорошо подтверждаются для широко диапазона эволюции пузырьков соответствием экспериментальным данным и численным результатам [15-17]. Здесь pg0 - начальное давление неконденсируемого газа внутри пузырька, pv - давление насыщенного пара, V0 - начальный объём пузырька, σ - коэффициент поверхностного натяжения, функция C ◦ = = C(x, t) - удвоенная локальная средняя кривизна в точке x - вводится для учёта возможного изменения формы пузырька. Константа k в данном случае определяется отношением удельных теплоёмкостей, связанных с составом смеси газа внутри пузырька; как известно, для идеального газа k = 1.4. Дополнительные трудности возникают при моделировании динамики дисперсной фазы в виде паровых пузырьков. Экспериментально установлено, что разрушение паровых пузырьков в результате обжатия внешним давлением может происходить уже при достижении ими размера не более 10 % от первоначального, и при этом начиная с некоторого критического размера они становится неустойчивыми [18]. Поскольку время разрушения tb ≈ 0.915R0 (ρ0 /|p0 |)1/2 пузырька радиуса R0 , содержащего пар с плотностью ρ0 и давлением на фронте волны p0 , оказывается гораздо меньшим характерного времени действия импульсного давления в жидкости, вполне обоснованным является использование ограничения на размер обжатого пузырька с последующим переходом к модели однофазной жидкости. Рэлей-Тейлоровская неустойчивость паровых пузырьков может контролироваться по текущим значениям чисел Рейнольдса, Бонда и Вебера [6, 19]. Для деформируемой поликристаллической среды с повреждениями (i = 3) введём в рассмотрение второй (мезо) уровень ее микроструктуры.2 Течение 2 Первые признаки разрушения поликристаллических сред (металлов) с размерами, наблюдаемыми оптической микроскопией на мезоуровне, как раз и проявляются в виде скопления пор и образования микротрещин [1, 21]. 372 О волновой динамике повреждаемых оболочек . . . среды во времени примем в виде жидкости с очень вязкой несущей фазой при скоростях Rep 1, а микроповреждения - в виде эквивалентных по размерам микропор сферической формы - тех же пузырьков, заполненных газом, однородно и изотропно с макроскопической статистической точки зрения распределённых по объему. Состояние газа в каждом пузырьке полагается изотермическим, а инерционные эффекты от движения стенок пузырьков и движения пузырьков относительно жидкости - пренебрежимо малыми. Разница между давлениями в несущей и дисперсной фазах уравновешивается вязкими силами в несущей фазе, которая в самом общем случае принимается сжимаемой. Тензор напряжений (2) будет иметь вид σjk = Sjk + p · δjk , где производные по времени от компонент напряжений и его девиатора Sjk принимаются в смысле Яумана-Нолла ∇ = S˙ jk - Sjr · ωrk - Skr · ωrj Sjk (8) в используемом лагранжевом описании среды с конвективными координатами X i . Здесь S˙ jk - полная производная напряжений Коши по времени, 1 ωjk = (vi,j - vj ,i ) 2 - тензор вихря скорости. Для сжимаемой при ударе деформируемой среды в качестве уравнения состояния (5) могут быть использованы известные соотношения Ренкина- Гюгонио на фронте распространяющейся ударной волны и уравнения Ми- Грюнейзена в форме [20] p = ρ0 Γ0 cv θm (1 + εd )Γ0 +1 ; θm = θm0 exp 2aεd /(1 + εd ) 1 + εd 2(Γ0 -a-1/3) . (9) Здесь Γ - коэффициент Грюнейзена; θm - температура плавления (θm0 - при нормальном давлении); a - постоянная материала; εd - мера сжимаемости повреждаемой среды, определяемая из условия εd = 1/J - 1; J = δVt /δV0 - якобиан преобразования между исходной C0 и текущей Ct конфигурациями, занимаемыми областью Ds × Dt . Для слабосжимаемых сред при относительно невысоких скоростях нагружения, обычно реализуемых на практике, можно ограничиться так называемым квазиакустическим приближением в форме, подобной (6а), когда ударная адиабата полагается совпадающей с изоэнтропой расширения, а давление является только функцией объёма. Рассматривая дисперсную фазу в виде пор - пузырьков газа, будем учитывать повреждаемость деформируемой среды как вследствие роста и слияния существующих пор вплоть до образования макротрещины, так и вследствие зарождения и эволюции новых. Этот процесс с учётом принятых допущений можно описать в виде непрерывной скалярно значимой функции пространства и времени ξ = ξ(xi , t). 373 П е т у ш к о в В. А. За меру повреждаемости ξ для всех (xi , t) ∈ Ds × Dt примем относительный объём микропор, содержащихся в элементарном объёме δV : ξ = δVd /δV, где δVd - часть его, заполненная микропорами. Уровни накапливаемых во времени повреждений определяются скоро˙ которая складывается из скоростей зарождения - стью повреждаемости ξ, ξ˙n и развития - ξ˙g микроповреждений, ξ˙ = ξ˙n + ξ˙g . (10) Для сжимаемой среды повреждения, накопленные за время ∆t, могут быть определены в предположении случайного характера зарождения микропор [10]: ξ(t + ∆t) = 8πN t Rn3 ∆t + ξ(t) exp 3∆t(p - pg )/4γ ; Nt Nt p>p0 p p0 = N0t exp (p - pn )/P1 ; (11) = 0, где p - растягивающее давление в среде; pn , pg - пороговое давление зарождения и роста микропор соответственно; Rn - параметр распределения размеров вновь образованных микропор; γ - вязкость материала; P1 и N0t - параметры материала; N t - скоростная функция числа зарождающихся микропор. Экспериментально установлено, что начальный уровень повреждённости поликристаллических конструкционных металлов обычно не превышает 3-4 %. Для слабосжимаемой среды изменение повреждённости за счёт эволюции уже существующих микроповреждений сферической формы может быть описано уравнением 1 ξ˙g = ψ(θ, ξ)(Σg - σg∗ (εpi , θ, ξ)). γ (12) Здесь при t = 0 ξ˙g = 0, ξg = ξ0 и ψ(θ, ξ) - функция роста микропор в рассматриваемой среде с учётом их взаимодействия. При сжатии пор ξ˙g < 0 и ψ(θ, ξ) = 1, σg∗ (εpi , θ, ξ) - пороговое значение напряжения для роста пор, зависящее от упрочнения, температуры и текущего уровня повреждённости, Σg = b1 I1 + b2 (I2 )1/2 + b3 (I3 )1/3 - интенсивность напряжений, возникающих при росте пор, основанная на инвариантах Ii напряжённого состояния в среде с порами, bi - параметры материала среды, i = 1, 2, 3. Образование новых микроповреждений (микропор, или «пузырьков») происходит в соответствии с нормальным законом распределения ξ˙n = Bp (σ˙ y + σ˙ m ) + Dp ε˙pi . (13) Здесь Bp и Dp - коэффициенты, значения которых зависят от выбора критерия появления микроповреждений: 374 О волновой динамике повреждаемых оболочек . . . - если используется пороговое напряжение σg∗ , то Bp = σ ∗2 FN √ exp - 2 , 2SN SN 2π Dp = 0, где σ ∗ = σy + σ m - σg∗ ; - если пороговая деформация ε∗g , то Dp = FN ε∗2 √ exp - 2 , 2SN SN 2π Bp = 0 и ε∗ = εpi - ε∗g ; FN , SN - параметры нормального закона распределения; ε˙pi - скорость изменения эффективной пластической деформации среды, которая определяется в данном случае из условия σ ij ε˙ij = (1 - ξ)σy ε˙pi ; компоненты тензора σ ij определяются уравнениями состояния для повреждаемой вязкопластической среды; 1 σ m = σ ij δij = p, 3 значения σg∗ или ε∗g , при которых происходит зарождение микропор, определяются экспериментально. Скорость изменения предела текучести материала σ˙ y может быть найдена из условия течения Мизеса для пористой среды. Для повреждаемой, зависящей от скорости нагружения деформируемой среды с анизотропным упрочнением это условие запишем в следующем виде: 1 f (σjk , θ, ξ) = Sˆjk Sˆkj - k2 (εpi , θ, ξ) = 0, 2 k(εpi , θ, ξ) = σys (εpi , θ)φ(ξ, ξF ). (14) Здесь k(εpi , θ, ξ) - размер поверхности текучести; σys (εpi , θ) - предел текучести; φ(ξ, ξF ) = (ξF - ξ)/(ξF - ξ0 ) - функция повреждённости среды; ξ0 - начальный, а ξF - разрушающий уровень повреждений соответственно; Sˆjk = Sjk - ρjk , ρjk - тензор микронапряжений, определяющий смещение поверхности текучести в соответствии с выражением ρ˙ jk = H ε˙pjk - 1 ∂k p ˆ ε˙ Sjk ; 2k ∂εplm lm H = H(εpi , θ) - угол наклона касательной к кривой деформирования в точке с накопленным уровнем деформации εpi . 375 П е т у ш к о в В. А. Соответствующие скорости вязкопластических деформаций ε˙pjk определяются выражением γ ∂ f˙ ε˙pjk = s Φ(F ) , F > 0, (15) σy ∂σjk где F = F (f, σys , εpi , ξ, θ); γ - вязкость материала деформируемой среды, зависящая от температуры θ. Связь между составляющими скоростей напряжений σ˙ jk и деформаций запишем в виде обобщённого закона Гука ˙ σ˙ = C : (ε˙ - ε˙p ) - αθI, где тензор упругости 2 C = 2µI + K - µ I ⊗ I = C(ξ, θ) 3 является функцией температуры и накопленных повреждений, ε˙ - скорость полной деформаций, α - коэффициент температурного расширения, I - единичный тензор. Фактом разрушения деформируемой среды является достижение параметром повреждаемости ξ предельного уровня ξF . В этом случае, как следует из (14), k = k∗ (εpi , θ, ξ) ξ=ξ = 0 F и материал конструкции теряет несущую способность. Для большинства поликристаллических металлов значения ξF находятся в диапазоне от 18 до 30 %. Вместе с тем вследствие накопления повреждений и возможного локального повышения температуры может происходить разупрочнение среды и её разрушению предшествовать образование локализованных полос сдвиговых деформаций. Учёт сдвигового механизма разрушения может быть принят в форме ηεpi /εF 1, где η - коэффициент объёмности напряжений, εpl - накопленная во времени пластическая деформация, а εF - предельная при разрушении деформация. Таким образом, краевая задача (1)-(5) распадается на систему уравнений, описывающих макроскопические процессы в деформируемой среде и связанную с этими процессами эволюцию микроповреждений. Для деформируемых сред в виде тонкостенных оболочек, широко используемых в конструкциях, взаимодействующих с жидкостями, может оказаться более эффективным с вычислительной точки зрения использование обобщённой структурной модели Афанасьева-Бессилинга, развитой на примере описания вязкопластического течения повреждаемых сред [10, 21]. В этом случае введённый ранее элементарный объём δV полагается состоящим из совокупности N связанных между собой первично изотропных структурных частиц - идеально упруговязкопластических субобъемов δv k , свойства которых определяют поведение среды в целом: N δv k . δV = k=1 376 О волновой динамике повреждаемых оболочек . . . Оптимальное число N таких объёмов может быть установлено из гистограммы статистического распределения пределов текучести, полученной измерениями микротвердости материала оболочки. Для большинства поликристаллических металлов оно не превышает 4-6 [21]. Это число отождествляется далее с числом слоёв оболочки по её толщине. Рассматривая исходную оболочку как многослойную, удаётся получить корректный характер изменения решения вдоль её вырожденного размера - толщины h. 2. Схемы численного решения задачи. Краевая задача (1)-(4), (8)-(13), описывающая высокоскоростные волновые процессы в жидкости с пузырьками газа или пара и в деформируемой среде с повреждениями, является сильно нелинейной. Ее решение возможно только численными методами, для выбора которых корректность постановки задачи играет важную роль. Если существование решения обычно следует из физической сущности задачи, то доказательство его единственности представляет самостоятельную и очень сложную для многомерного случая проблему. В качестве косвенного подтверждения корректности постановки задачи при ее численном решении в работе используется третье условие по Адамару - о непрерывной зависимости задачи от исходных данных, как это предложено А. Ф. Филипповым [13]. Задача решается численно по консервативным схемам, гарантирующим точность и вычислительную устойчивость, хорошо апробированным ранее и представленным, например, в [8-12, 20]. При этом используется известный подход, основанный на разложении решения по процессам, или метод «расщепления» [13], в силу различного характера, а также различных пространственных и временных масштабов явлений и процессов, происходящих в такой системе. Двухфазная жидкость аппроксимируется системой объёмных конечных элементов ∆k таким образом, что M h Ω≈Ω = ∆k , k=1 где M - общее число конечных элементов в жидкости и деформируемой среде и ∆l /∆m = ∅, если l = m. Деформируемая среда, взаимодействующая с жидкостью, аппроксимируется конечными элементами или разностной сеткой. Конечные элементы и/или разностная сетка, примыкающие по обе стороны к границе контактного разрыва Sc (xk , t), имеют на ней совпадающие узлы. Временной слой Dt = (0, τ ) протекания каждого исследуемого процесса аппроксимируется разностной сеткой Dt = {tn : tn = n∆t, n ∈ Z, 0 n∆t τ }. Интегрирование по времени осуществляется с использованием односторонней разностной схемы. Объёмные конечные элементы в жидкости рассматриваются в качестве ячеек, содержащих пробную дисперсную частицу (условный пузырек газа или пара), объём которой равен суммарному объёму находящихся в нем пузырьков. При этом полагается, что все пузырьки, содержащиеся в такой ячейке, 377 П е т у ш к о в В. А. колеблются одинаково. Изменение давления в жидкости связывается с изменением объёма находящейся в ней дисперсной фазы. При распространении ударных волн в кавитирующей жидкости парогазовые пузырьки могут изменять свою первоначальную форму, переходить в неустойчивое состояние и разрушаться. Учёт этих явлений оказывается чрезвычайно сложным для моделирования, поэтому в расчётах обычно вводится ограничение на минимальные размеры газового и парового пузырьков, по достижении которых они разрушаются. Значение шага интегрирования по времени ∆t для выбранных схем аппроксимации решения определяется с использованием известных критериев и уточняется в процессе численного эксперимента. 3. Численные результаты и обсуждение. В качестве примера рассматривается динамика канала бесконечной длины - пространства, образованного между двумя коаксиально расположенными деформируемыми оболочками и заполненного сплошной или двухфазной сжимаемой жидкостью - водой с пузырьками воздуха или пара (рис. 1, a). Рис. 1. [Figure 1] Внешняя оболочка подвергается действию ударной нагрузки (например, взрыва) с параметрами Pb = P b cos(φ), 0 φ π (рис. 1, c), где давление ∗ -3 P b действует в течение времени t = 0.26 · 10 с: P b = kp t при t < 0.5t∗ и P b = kp (t∗ - t) при 0.5t∗ < t < t∗ и достигающей максимального значения Pb∗ = 20 МПа, при t = 0.5t∗ (рис. 1, d). Внешняя и внутренняя оболочки имеют следующие размеры: R2 /R1 = = 0.5, h2 /h1 = 0.5 и h2 /R2 = 0.05, где R2 = 1 м, и изготовлены из стали со следующими физико-механическими свойствами: плотность ρ = 7810.0 кг/м3 , модуль упругости E = 180.0 ГПа, коэффициент Пуассона ν = 0.3, предел текучести σT = 200.0 МПа. Физико-механические характеристики материала оболочек и диаграммы его деформирования, зависящие от скорости нагружения и температуры, так же как параметры уравнений состояния и кинетики повреждений (13)- (15), приведены в [10]. Исходная повреждаемость ξ оболочек полагается равной 3 %. В начальный момент времени при нормальных условиях (P = 0.1 МПа) сплошная жидкая среда и двухфазная газожидкостная среда имеют комнатную температуру T0 = 293 K, а парожидкостная среда - температуру насы378 О волновой динамике повреждаемых оболочек . . . щения T0 = Ts (p0 ) = 373 K. Начальная плотность пара совпадает с плотностью насыщенного пара ρ020 = ρ02S . Дисперсная фаза в жидкости представляет собой парогазовые пузырьки сферической формы одинакового размера, распределённые в жидкой среде равномерно. Принимается, что в 1 см3 воды находится 1 000 пузырьков и объёмная концентрация дисперсной фазы составляет α2 = 0.01. Каждый отдельный пузырёк имеет начальный размер R0 = 0.134 · 10-3 м. При сжатии газового и парового пузырьков вводится ограничение на их минимальные g v размеры, соответственно Rmin = 5.0 · 10-5 м и Rmin = 5.0 · 10-6 м. Значения теплофизических параметров жидкости и дисперсной фазы заимствованы из [22] и считаются постоянными. Эффективные значения чисел Нуссельта Nui для вещества каждой из фаз приняты, как в [9]. Для воды при T0 = 373 K Nu1 = 195.0, для ее пара при тех же условиях Nu2 = 14.0, для воздуха при T0 = 293 K Nu2 = 16.5. В силу симметрии задачи рассматривается половина конструкции с очевидными в данном случае граничными условиями для жидкости и оболочек относительно плоскости симметрии и на торцах оболочек. Условия на поверхности контактного разрыва Sc (xk , t) между стенкой канала и жидкостью принимаются в форме (4). Начальное состояние для оболочек и жидкости в момент времени t = 0 полагается свободным от деформации. Интегрирование по времени осуществляется по явной схеме с шагом ∆t, значение которого принимается в соответствии с требованиями устойчивости вычислительного процесса [10, 11]. В качестве модельной вначале рассмотрим задачу, аналогичную представленной на рис. 1, с пространством между оболочками, заполненным водой при комнатной температуре. Известны результаты решения этой задачи в предположении упругого деформирования оболочек под действием того же импульса треугольной формы, как и на рис. 1, d, но меньшей продолжительности действия [14]. На рис. 2 показаны смещения во времени лобовых точек 1 и 5 контура наружной и внутренней оболочек (см. рис. 1, b): сплошные линии (1) - результаты, полученные по предложенной модели, точками (2) - результаты работы [14]. Из сравнения результатов следует, что представленная математическая модель позволяет адекватно и с приемлемой точностью описывать особенности импульсного взаимодействия деформируемых сред с жидкостью. Рис. 2. [Figure 2] 379 П е т у ш к о в В. А. На рис. 3 представлены результаты моделирования волновых процессов в сплошной (рис. 3, a) и кавитирующей с пузырьками газа (рис. 3, b) и пара (рис. 3, c) жидкостях под действием импульса давления (взрыва), заданного в форме, представленной на рис. 1, c, d. Кривые 1-5 соответствуют расчётным точкам конечно элементной сетки в поперечном сечении канала, показанном на рис. 1, b. Из рис. 3, a следует, что в сплошной среде волны давления мало отличаются по форме от внешнего воздействия. Волна давления с максимальной амплитудой Pmax = Pb∗ достигает стенки канала (точка 5 на рис. 1, b) в момент времени t = 0.35 · 10-3 c и отражается от его поверхности в момент времени t = 0.55·10-3 c. Скорость распространения волны давления в сплошной среде составляет порядка 1450 м/с. В момент времени 7.4·10-4 c отраженная волна разрежения догружает внешнюю оболочку канала (точка 1 на рис. 1, b). Из рис. 3, b следует, что наличие пузырьков газа малой концентрации в сжимаемой жидкости существенно меняет характер распространения в ней волн по сравнению со сплошной средой. Распространяющаяся от внешней оболочки волна в момент времени t = 0.37 · 10-3 c разделяется на две (кривая 4). Первая волна с максимальной амплитудой p = 35.0 МПа перемещается со средней скоростью v = 750 м/с и действует на внутреннюю оболочку в течение ∆t = 0.055 · 10-3 с, вторая с максимальной амплитудой p = 15.0 МПа перемещается со средней скоростью v = 647 м/с и действует на ту же оболочку в течение ∆t = 0.05 · 10-3 с. В этом случае форма возникающей в жидкости ударной волны отлича- Рис. 3. [Figure 3] 380 О волновой динамике повреждаемых оболочек . . . ется от формы внешнего воздействия. Амплитуда давления в волне меняется с частотой ωr = 1.089 · 105 с-1 , которая согласуется с частотой собственных колебаний изолированного пузырька воздуха в жидкости (воде) ωr = 1.533 · 105 с-1 , полученной в соответствии с формулой Миннаерта [6]. Отражения волны давления от поверхности внутренней оболочки не происходит, так как возникающая волна разрежения исчезает за счёт вызванного ею увеличения объема дисперсной фазы. В соответствии с этим внешняя оболочка не догружается отраженной волной, а внутренняя не разгружается ею. Максимальное воздействие среды на внешнюю оболочку (кривая 1 на рис. 3, b) имеет место в момент времени t = 1.35·104 c с амплитудой давления сжатия P = 18.9 МПа, а на внутреннюю (кривая 5 на рис. 3, b) - в момент времени t = 5.7 · 10-4 c с амплитудой давления сжатия P = 35.0 МПа и приходится на их лобовые части. Из рис. 3, c следует, что волновые процессы, происходящие в парожидкостной среде, также имеют ряд характерных особенностей. Дисперсные частицы за фронтом ударной волны обжимаются и исчезают, среда становится сплошной. В результате схлопывания отдельного парового пузырька возникает скачок давления разрежения, который накладывается на давление во фронте волны и уменьшает его амплитуду. Начиная с некоторого момента времени результирующая волна становится волной разрежения и исчезает за счёт увеличения объёма дисперсной фазы. До внутренней оболочки волна давления не доходит. Не происходит разделения распространяющейся ударной волны, а имеют место периодические колебания амплитуды давления на её фронте с частотой ωr = 1.095 · 105 c-1 , которая согласуется с частотой собственных колебаний изолированного пузырька пара в жидкости (воде) ωr = 1.24 · 105 c-1 , полученной в соответствии с формулой Миннаерта. Средняя скорость движения волны составляет порядка 920 м/с. Максимальное воздействие давление сжатия на внешнюю оболочку имеет место также в момент времени t = 1.35 · 10-4 c (кривая 1 на рис. 3, c), но с несколько большей амплитудой P = 20.8 МПа, и продолжается значительно дольше по сравнению с газожидкоcтной средой. С течением времени это воздействие превращается в давление разрежения с максимальным значением P = 7.7 МПа в момент времени t = 3.3 · 10-4 c и продолжается в течение ∆t ≈ 7.0 · 10-4 с. В случае сплошной среды догрузка внешней оболочки отраженной волной происходит с максимальной амплитудой P = 10.0 МПа и имеет существенно меньшую продолжительность ∆t ≈ 2.0 · 10-4 с. Выявленное поведение сжимаемой жидкости с дисперсной фазой в виде парогазовых пузырьков малой концентрации является характерным для подобных сред и подтверждается численно и экспериментально в ряде работ других авторов [3, 15, 17]. Различия в характере волновых процессов, обусловленные наличием дисперсной фазы и ее физической природой, существенно влияют на динамическую реакцию оболочечных конструкций (см. рис. 4). Здесь показаны распределения во времени давления вдоль контура внутренней оболочки канала в точках 1-11 для сплошной жидкости (рис. 4, a) и жидкости с пузырьками газа (рис. 4, b), а также соответствующие картины деформирования контура оболочки, построенные для одного и того же момента времени t = 9.75·10-4 с 381 П е т у ш к о в В. А. Рис. 4. [Figure 4] от начала нагружения. Поскольку для жидкости с пузырьками пара ударная волна не достигает внутренней оболочки, паровая фаза выступает для нее в роли защитного экрана. На рис. 4, с приведена картина деформирования внешней оболочки для того же момента времени. Максимальное смещение внутренней оболочки (точка 1 на рис. 4, b) под действием ударной волны сжатия происходит для жидкости с пузырьками газа и составляет порядка 6.3 · 10-3 м. Для сплошной жидкости (рис. 4, a) это смещение не превышает 3.0 · 10-3 м. С течением времени пиковые нагрузки перемещаются вдоль контура поперечного сечения оболочки и падают до нуля. Для жидкости с паровыми пузырьками (рис. 4, c) максимальные смещения внешней оболочки вдоль осей X1 и X2 (см. рис. 1, b) оказались сопоставимыми с аналогичными смещениями для сплошной жидкости и составили величины порядка 8.0 · 10-3 и 5.0 · 10-3 м соответственно. Они более чем в три раза превышают смещения внешней оболочки для случая жидкости с пузырьками газа. На рис. 5 красным цветом3 выделено деформированное состояние оболочек канала, заполненного сплошной жидкостью, для рассматриваемого момента времени. Здесь же показаны распределения пластических деформаций по толщине внешней оболочки в её наиболее нагруженном сечении A и зона поврежденности, соответствующая максимальному уровню ξ порядка 6% 3 382 Онлайн-версия статьи (doi: 10.14498/vsgtu1435) содержит цветные рисунки. О волновой динамике повреждаемых оболочек . . . Рис. 5. [Figure 5] (выделена голубым цветом). При этом максимальная повреждённость ξ внутренней оболочки составляет менее 4 %, а зона поврежденности располагается симметрично, аналогичной зоне на внешней оболочке. Выводы. Представлены математическая модель и результаты численного моделирования распространения ударных волн в нелинейно деформируемой повреждаемой среде и взаимодействующей с ней кавитирующей сжимаемой жидкости. Показано существенное влияние дисперсной фазы в виде парогазовых пузырьков малой концентрации в жидкости и микроповреждений в виде пор на уровни и кинетику повреждаемости в деформируемой среде, образованию в ней предельных состояний. Выявлена роль дисперсной фазы в жидкости в виде паровых пузырьков как защитного экрана для конструкций, подвергаемых воздействию интенсивных ударных волн (взрывам) в жидкости. Полученные в работе результаты создают основу для изучения условий формирования и распространения кавитационного облака в задачах управляемой кавитации и оптимального (рационального) проектирования объектов техники, подвергаемых действию кавитации.

About the authors

Vladimir A Petushkov

A. A. Blagonravov Mechanical Engineering Institute RAS

Email: pva_imash@bk.ru
(Dr. Phys. & Math. Sci.; pvaimash@bk.ru), Professor, Lab. of Mathematics Modeling 4, M. Khariton'evskii per., Moscow, 101990, Russian Federation

References

Supplementary files