Расчёт течений в суспензиях

Вестник Челябинского государственного университета. 2010. № 24 (205).

Физика. Вып. 8. С. 23-30.

Н. Б. Волков, А. Е. Майер, В. В. Погорелко, А. П. Яловец

расчёт течений в суспензиях

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

ключевые слова: модель гетерогенной среды, скорость звука, смещения включений, суспензия.

В последние годы особый интерес представляет моделирование волновых процессов в многофазных средах . Одним из таких процессов, интересным с точки зрения практических приложений и моделирования, является распространение звука в двухфазной среде . В частности, в работах [1-3] изучалось движение частиц, взвешенных в жидкости или газе, в результате действия звуковых волн . В настоящее время большое внимание уделяется ультразвуковому исследованию суспензий [4] . В составе суспензии можно выделить жидкую несущую среду и включения различной объемной доли с размерами от 10 нм до 100 мкм . Взаимодействие компонент суспензии (механическое, тепловое, электрическое, химическое и пр .) в значительной мере определяет ее свойства. В настоящей работе нас в первую очередь будет интересовать механическое взаимодействие компонент и его влияние на поведение суспензии в условиях динамического нагружения

Для описания течения суспензий в условиях динамического нагружения будем опираться на получивший широкое распространение подход многоскоростных взаимопроникающих и взаимодействующих континуумов [57] . В рамках этого подхода поведение каждого компонента определяется законами сохранения Существенной проблемой при реализации данного подхода является определение удельного объема компонент суспензии В существующих моделях эта проблема обычно решается за счет упрощающих предположений, например, о несжимаемости включений [5-6] либо о равенстве давлений несущей фазы и включений [7] . При рассмотрении динамических процессов в суспензиях такие предположения часто являются неприемлемыми

В данной работе в рамках подхода многоскоростной гетерогенной среды предложена мате-

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

Модель гетерогенной среды. Рассмотрим многокомпонентную среду, состоящую из несущей среды (матрицы) и включений разного размера и термодинамического состояния . Будем считать, что размеры включений в смеси во много раз больше расстояний между молекулами . Данное допущение позволяет рассматривать гетерогенную среду как совокупность многоскоростных взаимодействующих континуумов [5-7] .

Пусть V — некоторый элемент объема гетерогенной смеси; V3) — объем j-компоненты в элементе V; а(/) — объемная доля j-компоненты .

Сформулируем основные положения модели .

1 . Усредненный по физически малому объему тензор деформации гетерогенной смеси ицс имеет вид

% =Е“

7

где ы*к) — тензор деформации j-компоненты .

Данное выражение получается путем стандартной процедуры усреднения тензора деформации по элементарному объему, включающего подобласти, содержащие лишь одну компоненту. Из (1) естественным образом вытекает аддитивность приращений объемов компонент среды

Ь¥ = ^Ь¥(3

]

Для приращений вектора смещения можно выполнить аналогичную процедуру усреднения и получить следующее выражение для тензора

скоростей деформации = ^а(3)у(£}, где ^—

3

тензор скоростей деформации j-компоненты .

2 . Элементарная работа А((], совершаемая j-компонентой над всей гетерогенной средой в единицу времени в единице объема, имеет вид

где — тензор напряжений j-компоненты .

Действительно, А} должна определяться напряженным состоянием данной компоненты и деформацией среды в целом и, кроме того, должна удовлетворять выражению для элементарной работы, совершаемой при деформации в единицу времени в единице объема гетерогенной сме-

си

3 . Все компоненты среды равноправны и движутся в самосогласованном поле напряжений , создаваемом всеми компонентами среды:

= 1

а

•

Для получения (3) в качестве термодинамического потенциала гетерогенной среды выберем внутреннюю энергию единицы ее объема . Поскольку данный термодинамический потенциал обладает свойством аддитивности, запишем:

а

№), *(])),

где Е(']'>, Б(-/) — внутренняя энергия и энтропия единицы объема j-компоненты . Так как компоненты тензора напряжений имеют вид

°,к =

Кдигк )

[8], то, записав с учетом (2) основ-

ное термодинамическое соотношение для деформируемых тел

аЕ (]) =т и) (]) +аи) аи,к,

можно из (4) получить (3) .

4 . Компоненты в гетерогенной среде стремятся к локальному термодинамическому равновесию . Время установления локального термодинамического равновесия между компонентами в некотором малом объеме конечно и определяется геометрией компонент и их объемной долей . При описании процесса установления равновесия необходимо учесть релаксацию напряжений компонент и обусловленное этим процессом изменение плотности каждой компоненты и ее внутренней энергии .

Сформулируем систему уравнений механики для компонент гетерогенной среды в интегральной форме для некоторого элемента объема V в подвижной системе координат, связанной с движением центра масс . Для учета взаимодействия компонент в уравнение движения и в уравнение для внутренней энергии необходимо добавить слагаемые, описывающие теплообмен и тре -ние между компонентами гетерогенной среды Полная производная по времени примет вид

йс/йг = д/дг + ( усУ) ,

где V- = ^ V(**^ ^ т( — скорость цент-

] / ] ра масс, V(1 \ ш(~1 1 — соответственно скорость и масса j-компоненты в объеме V В дальнейшем для обозначения производной по времени в системе центра масс будем использовать также

точку. Обозначим через ^^ = V(^^ относительную скорость j-компоненты .

Уравнение непрерывности для j-компоненты полностью эквивалентно уравнению непрерывности, приведенному в [5], и имеет вид

Жст

ж

где р(^^ ^ /V = а^^\ р(7'^ = т(7'^/V(7'^ —

соответственно приведенная и истинная плотности j-компоненты; правая часть уравнения описывает поток массы через поверхность объема V, обусловленный переходом в систему центра

масс . При записи (5) учтено, что — = (Уу ) •

Уравнение движения для j-компоненты в поле напряжений (3) имеет вид

Ж

= а

Ъгк^к

-V£ В

к

где правая часть уравнения описывает соответственно поверхностную силу, силы трения между компонентами, поток импульса через поверхность элемента объема V, связанный с переходом в подвижную систему координат, Б{ ]к) — коэффициент трения между j- и к-компонентами, для которого справедливо соотношение В(]к) = В(к]).

Уравнение для внутренней энергии j-компоненты можно записать в виде

—(т(])и(])) = |а(3^к^У1кёУ - а(3) ф $)^ +

л У ^

+{Е Q(]к ) (т (3) - т (к) )у+

V к

р

V к Яу

где и() — внутренняя энергия единицы массы j-компоненты; слагаемые в правой части уравнения соответственно учитывают работу компоненты при деформации смеси (2), теплопроводность, теплообмен между компонентами, работу сил трения, поток энергии через поверхность элемента объема V, связанный с переходом в подвижную систему координат. В уравнении

введены обозначения: = —к — по-

дхк

ток тепла, обусловленный теплопроводностью, где к(']> — коэффициент теплопроводности j-компоненты; ^  — коэффициенты теплообмена между j-компонентой и к-компонентой;

Р1

(Д) ----

доля тепла, передаваемая _)-компоненте

при трении с к-компонентой [5], причем Р(к) = = 1 - Р1^1 .

Удобно представить =- Р(] ')дік + ^, где

рО) — шаровая часть, БЦ] — девиатор напряжений для _)-компоненты . Шаровая часть тензора напряжений находится из уравнения состояния, а уравнение для девиатора тензора напряжений _)-компоненты имеет вид

&

1

где

уЦ) = уЯ) -- уУ )5

у1к

II иїк з

модуль

сдвига для j-компоненты . Для жидкой компоненты модуль сдвига будем полагать равным нулю Пластическое течение описывается путем сохранения девиатора напряжений на пределе текучести (условие текучести Мизеса) [9] .

Запишем уравнение для изменения объема j-компоненты . Изменение объемной доли компоненты за счет перетекания массы вещества через поверхность объема V определяется выражением

. (ч т ч (л т ч

а 7 = ——— = а

т

О')

Поскольку V(* ^ = а( * ^V + а(* ^V, уравнение для изменения объема j-компоненты (без учета релаксации) будет иметь вид

V

И У)

V т  — +

Важным физическим процессом в многокомпонентной среде является выравнивание (релаксация) напряжений в компонентах Поэтому система уравнений (5)-(9) дополняется уравнениями релаксации компонент, скорость которой зависит от геометрии включений

Рассмотрим в сферической системе координат отдельное неподвижное включение радиуса Я(/), находящееся в жидкости с давлением .

За счет разности напряжений в жидкости и во включении происходит изменение объемов компонент, стремящихся к равновесию по напряжениям . В течение времени & в результате деформации включения возмущается слой шириной с( 1)Ы, , где С(1) — продольная скорость звука в j-компоненте .

Уравнение движения для вещества отдельного включения можно записать в виде

р( 1) V1) = д°гг

дг

где г — координата в сферической системе коор-

і)

динат, Уг — радиальная составляющая скоро -сти вещества включения

Начальные условия имеют вид

(г < !„) = а<11,.((1„),

(г > Яи>, 1„) = Р<і> ().

Решая уравнение движения для вещества отдельного включения, можно получить приближенное выражение для изменения плотности компонент:

р( (+1) = р( ()([п)+ р{Ь) + ст( ] ) -^

+(1 -е Хг),

сеЛТ

где С2 = (с{] (с{Ь))2

где сеГГ (С1 ) + V (ь) _(]) (с ) >

и) 1 я( У} }

а тУ) =----------------—

3 с2

- характерное время релак-

е//

сации .

При этом необходимо учесть изменение объема и внутренней энергии компонент

Коэффициенты трения В(к') для включения сферической формы находились по формуле Стокса . Коэффициенты теплообмена между

включениями и жидкостью были найдены из решения уравнения теплопроводности в окружающей сферическую частицу жидкости

Из уравнений (5)-(10) легко получить уравнения для среды в целом, для чего следует просуммировать каждое из уравнений по всем компонентам Получившиеся уравнения имеют обычный вид, как и в односкоростном случае, но зависят от относительного движения составляющих . Интегрируя получившиеся уравнения по всему объему гетерогенной среды, учитывая, что в случае изолированной системы потоки массы, импульса и энергии через ее поверхность равны нулю, получим законы сохранения массы импульса и энергии

Дополняя систему уравнений (5)-(10) уравнениями состояния, мы получим замкнутую систему уравнений для суспензии, учитывающую взаимодействие фаз Данная система уравнений решалась методом, описанным в работе [10] .

Акустическое приближение. Пусть по суспензии вдоль оси z распространяется плоская продольная волна малой амплитуды . В этом случае уравнение движения (6) в акустическом приближении в эйлеровых координатах можно переписать в виде

дг

&

-XВ 1)

Рассмотрим случай, когда суспензия состоит из двух компонент: несущей среды (обозначим индексом «1») и включений радиуса Я^ (индекс «2») . Согласно (3), тензор напряжений в этом случае можно записать в виде

=(1 - а) Р

С

'

,

где с(1), с(2) — соответственно скорость звука в жидкости и продольная скорость звука в материале включения; е(1), е(2) определяют деформацию (относительное изменение объема) несущей среды и включений в системе центра масс

Изменение объемной доли компоненты а^ ) при фиксированном объеме V может быть обусловле -но перетоком ее массы через поверхность объема и процессом релаксации компонент к термоди-

намическому равновесию . При малой амплитуде волны перетоком массы компоненты через поверхность элемента объема V можно пренебречь и изменение объемной доли компонент обусловлено только релаксацией компонент к одинаковому давлению . Полная релаксация напряжений в рассматриваемой суспензии будет происходить, если выполняется условие Я / СI22 1/ ю,

где ю — частота волны .

Из сказанного следует, что в отсутствие релаксации деформации компонент имеют вид

(1) _ О (2) - о - к =1,2

1р

к=1,2

В случае мгновенной релаксации компонент по напряжению деформацию компонент можно записать в следующем виде:

(1) ЪУ 8 =80 +

В(2) =80 -

5У

V(2)'

В выражении (14) 5У — приращение объемов компонент за счет релаксации, которое можно найти, приравняв напряжения компонент:

р(1)с(1)2

ЪУ

V'

_ ~(2)(2)2

= Р

'I

ЪУ

--

V'

.

Учитывая оба предельных случая — (11) и (12), можно выражение для тензора напряжений в суспензии записать в виде

= О

к=1,2

ди(к) & '

Входящая в (16) константа В имеет вид

где Вд = 0 в отсутствие релаксации компонент, (1 -а)а (р(1)с(1)2 -р(2)сР)2}2

р(1) +р(2) (1 -а)р(2)с(2)2 + ар(1)с(1)2 в случае мгновенной релаксации

Подставляя выражения (16) с соответствующим коэффициентом (17) в уравнение движения (11), получим систему уравнений для векторов смещения и(1) и и(2):

Р

Я 21,

(1) ^ = (1 -а) Б

Яі2

Я 2« (2) (2) Я и

2„а)

Я 2 и

+ Р

2,(2)'

Я 2 и

Яг2

9 а Я 2 Д2 Яі

і

Яі2

= аБ

дУ !)

Яг 2

+ Р

Система уравнений (18) представляет собой волновые уравнения для суспензии, содержащей включения одного сорта. Будем искать решение (18) в виде плоской волны:

Подставляя (19) в (18), получим однородную систему уравнений, которая будет отличаться для двух рассматриваемых предельных случаев релаксации по напряжению только видом коэффициента В:

Л(1) В11 + Л(2) В„ = 0

12

где

Л"1 В21 + Л112 В22 = 0,

п 2/1 >ш2 9 1 а .

Вц — — ш + (1 — а)Бк + — —(1)- п—2іш

Р(2) 2 . . ~

В12 — (1 — а)^]уПк —Т “ЙТ /ш>

2 р(1) 'я2

9 1 а .

2 р(1) Я2

г, Р'*' гл »2 9 1 а .

В21 — <*^5) Ы — 2 р(2) ПЯ2 “

П 2 тлт2 9 1 а .

В99 — —ш +аУк +-----------тчтЛ—т ^ш.

22 2 р(2) V

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

Скорость звука в суспензии. Если частицы в суспензии малы так, что имеют место полное увлечение их жидкостью и мгновенная релаксация компонент по напряжению, то скорость распространения малых возмущений в такой среде можно вычислить аналитически Пусть плоская продольная волна распространяется вдоль оси z . Тензор напряжения о22 и продольную скорость звука в упругой среде можно записать в виде [8]

Я 2и

Яг 2

_9 л-аА

где К = К + 4ц/3, К, ц — соответственно модули всестороннего сжатия и модуль сдвига, игг — относительное изменение объема среды, р — массовая плотность . Для жидкости полагаем модуль сдвига равным нулю

Применяя, в соответствии с формулами (1) и (3), выражения (21) к суспензии и ее компонентам, можно найти К для суспензии:

К =а-К

'(2) и

и

и

где К = К(2 + 4|0,(2)/з, к(2\— моду-

ли всестороннего сжатия и сдвига включений, К1 — модуль всестороннего сжатия несущей среды, и(2, и<У2 — относительные изменения объема компонент, и22 = а• и^ + (1 - а) • и £? — относительное изменение объема всей среды, определенное по формуле (1), а — объемная доля частиц

Если при распространении волны с частотой ю в суспензии с размерами частиц Я выполняется условие Я / с(2) « 1/ ш, где с(2) — продольная скорость звука в материале частицы, то имеет место полная релаксация компонент среды по напряжению (при ю = 3 МГц для частиц железа имеем Я ^ 1 мм):

Учитывая, что массовая плотность суспензии равна р = а -р(2) + (1 - а)р(1), можно из (21)-(23) получить выражение для скорости звука в суспензии:

с = с

Рг

2

у

где

Р = р(2)/р(1), у = е(2)/с(1), р(2), р(1) -

истин-

.

ные плотности частиц и жидкости, с — скорость звука в жидкости, а — объемная доля включений . Если в формуле (24) положить у = с(2)/ с, где с(2) — объемная скорость

звука в материале включении, то она перейдет в формулу Вуда [11] .

Если частицы в суспензии такие крупные, что релаксация компонент по напряжению не происходит, то скорость распространения малых возмущений в такой среде тоже можно вычислить аналитически . В этом случае выполняется условие Я/с(2) ~ 1/га, релаксация компонент по напряжению отсутствует, и сжимаемость суспензии определяется в основном сжимаемостью жидкости, а не включений . В этом случае выражение для скорости звука в суспензии можно записать в следующем виде:

с = с

1

1 + (Ру-1)а

1 + (р-1)а '

Расчеты и обсуждение результатов. Справедливость формулы (24) проверялась ее сравнением с результатами, полученными в экспериментах [4], где была измерена зависимость скорости звука в воде, в которой присутствовали частицы SiO2 радиуса 30 нм, от их объемной доли . В данной работе отмечается, что массовая плотность частиц двуокиси кремния несколько меньше массовой плотности макроскопических объемов, которая обычно приводится в справочниках . Согласно данным [12], массовая плотность частиц 8Ю2, размеры которых меньше 75 мкм, достигает значений 2140 кг/м3, что существенно меньше плотности макроскопического объема (2650 кг/м3) . Здесь при сравнении с экспериментом использовалось значение плотности 2170 кг/м3 .

На рис . 1 проводится сравнение с экспериментальными данными [4] результатов, полученных по формуле (24), результатов решения системы уравнений (20), а также результатов расчетов по изложенной выше модели . Разброс экспериментальных значений в [4] объясняется сильной зависимостью скорости звука от температуры суспензии . Как видно из рис . 1, результаты расчета скорости звука по формуле (24) и на основе построенной здесь модели динамики суспензии хорошо согласуются с результатами эксперимента [4] . Отметим, что при определении скорости звука в суспензии на основе модели решалась задача о распространении волны малой амплитуды (102—103 Па) с частотой ю = 3 МГц, возбуждаемой на границе . Для расчетов были использованы двухчленное уравнение состояния воды и уравнение состояния Ми — Грюнайзена для включений

а

Рис. 1. Зависимость скорости звука в суспензии от объемной доли включений стекла ЗЮ2. Радиус включений К = 30 нм.

• — экспериментальные данные работы [4] — аналитическая кривая по формуле (24);

• — решение системы уравнений (20);

— численные расчеты

Формула (24) позволяет объяснить немонотонную зависимость скорости звука в суспензии от объемной доли включений . Действительно, скорость звука в суспензии определяется соотношением К /р. С увеличением объемной доли включений плотность суспензии р растет быстрее, чем уменьшается сжимаемость суспензии (К ) . Этим и объясняется начальное уменьшение скорости звука в суспензии при малых объемных долях включений . При больших объемных долях включений сжимаемость суспензии

(К ) уменьшается быстрее, чем растет плотность суспензии р, что приводит к увеличению скорости звука в суспензии

Приведем результаты численных исследований скорости звука в суспензии, состоящей из воды и включений железа различных размеров и объемной доли . Для определения скорости звука в суспензии моделировалось прохождение ультразвуковой волны с частотой ю = 3 МГц . Результаты расчетов приведены на рис 2

Отметим, что при малых радиусах частиц численное значение скорости звука, полученное на основе модели динамики суспензии, стремится к аналитическому решению (24) . Некоторое отличие численного и аналитического решения

при малых радиусах включений обусловлено использованием приближения полного увлече -ния включений жидкостью при выводе формулы (24) .

а

Рис. 2. Зависимость скорости звука в суспензии от объемной доли включений железа: 1 — полная релаксация компонент по напряжению (сплошная — расчет по формуле (24), маркеры о — решение системы уравнений (20) для R = 0,1 мкм); 2 — отсутствие релаксации компонент по напряжению (сплошная — расчет по формуле (25), маркеры • — решение системы уравнений (20) для К = 1 мм).

Численные расчеты: 3 — R = 0,1 мкм;

4 — R = 1 мкм; 5 — R = 10 мкм

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

На рис . 2 кривая 3 соответствует почти полному увлечению частиц жидкостью, кривая 5 — минимальному увлечению частиц . Кривая 4 при малых объемных долях, когда мала суммарная поверхность частиц и мала сила трения, соответствует случаю минимального увлечения частиц жидкостью и по мере увеличения объемной доли и, как следствие, увеличения силы трения приближается к режиму полного увлечения

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

их часть в • а объединилась в кластеры радиуса Яс = Я3 п (п — количество частиц в кластере), то количество таких кластеров в объеме суспензии V будет равно Лс = гЛ0 / п, где N — количество частиц радиуса Я в этом же объеме, соот-ветствующее объемной доле включений а Исследования зависимости скорости звука в суспензии от доли кластеров радиуса Яс в общем числе включений показали, что объединение включений в кластеры приводит к увеличению скорости звука в суспензии Если кластеры образуются из частиц малых размеров, то имеет место существенное увеличение скорости звука . Например, при фиксированной объемной доле включений а = 0,2 образование кластеров радиуса Яс = 1 мкм из частиц Я = 0,1 мкм приводит к увеличению скорости звука от 1140 до 1520 м/с при изменении доли в в интервале 0-1 . Если же кластеры образуются из крупных частиц, то ско -рость звука в такой суспензии меняется мало Так, в случае образования кластеров радиуса Яс = 10 мкм из частиц радиусом Я - 1 мкм скорость звука меняется в интервале 1540-1550 м/с при всех в .

На рис 3 показаны отношения амплитуды смещений включений к амплитуде смещения воды при различных объемных долях в зависимости от частоты волны, найденные из решения системы уравнений (20) для случая полной релаксации компонент по напряжениям

1 л 0.8 -0.6 -

<

0.40.2 -

о И—................ —. —.... —...... —«—г-1

0.1 1 10 100 1000

V, МГ ц

Рис. 3. Зависимость отношения амплитуд смещений включений и воды от частоты волны.

Радиус включений R = 0,1 мкм. Объемные доли включений: • — а(2 = 0,1; • — а(2) = 0,3;  — а(2) = 0,5;  — а(2 = 0,7;  — а(2) = 0,9

В области частот V < 0,3 МГц при любых объемных долях включений смещения включений и воды совпадают, т. е . имеет место полное увлечение частиц жидкостью при прохождении волны . В области частот V > 3000 МГц смещения включений при любых объемных долях становятся минимальными из-за их инерционности, поскольку массовая плотность (истинная) больше плотности жидкости .

В диапазоне частот 0,3 < V < 3000 МГц наблюдается переходная область, в которой происходит плавный переход от режима полного увлечения к режиму минимального увлечения частиц проходящей волной . В переходной области при фиксированной частоте волны увеличение объемной доли включений приводит к увеличению их смещения по отношению к смещению жидкости, что обусловлено увеличением силы трения между компонентами суспензии

При увеличении радиуса частиц характер рассматриваемой зависимости не изменяется, но происходит смещение картины в область меньших частот. Так, например, при увеличении радиуса частиц на порядок зависимости смещаются в область меньших на два порядка частот Такая закономерность объясняется тем, что сила трения между компонентами в единице объема суспензии, как видно из уравнений (18), обратно пропорциональна квадрату радиуса включений Случай отсутствия релаксации в двухкомпонентной суспензии реализуется для крупных частиц при частотах, когда Я / С>> 1/ ю. Например, при радиусе включений Я(/) = 1 мм данное неравенство выполняется при частотах, превышающих 5 МГц . В этом диапазоне частот движение включений соответствует режиму минимального увлечения

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

1 . Горьков, Л . П . О силах, действующих на малую частицу в акустическом поле в идеальной жидкости / Л . П . Горьков // Докл . Акад . наук СССР. 1961 . Т. 140, № 1 . С 88-91 .

2 . Каневский, И . Н . Постоянные силы, возникающие в звуковом поле / И . Н . Каневский // Акуст. журн. 1961 . Т. 7, № 1 . С. 3-16 .

3 . Шутилов, В. А . Основы физики ультразвука / В . А . Шутилов . Л. : ЛГУ, 1980. 280 с .

4 . Dukhin, A . S . Ultrasound for characterizing colloids / A . S . Dukhin, P. G. Goetz . Elsevier, 2002 .

5 . Нигматулин, Р. И . Динамика многофазных сред / Р. И . Нигматулин. М . : Наука, 1987

6 . Нигматулин, Р. И . Основы механики гетерогенных сред / Р И Нигматулин М : Наука, 1978

7. Куропатенко, В . Ф. Модель гетерогенной среды / В Ф Куропатенко // Докл Акад наук 2005. Т. 403, № 6 . С. 761-763 .

8 Ландау, Л Д Теоретическая физика / Л . Д . Ландау, Е. М . Лифшиц. М . : Физматлит, 1987 Т 7: Теория упругости

9 Уилкинс, М Л Расчет упругопластических течений / М . Л . Уилкинс // Вычисл . методы в гидродинамике М : Мир, 1967

10 . Яловец, А . П . Расчет течений среды при воздействии интенсивных потоков заряженных частиц / А П Яловец // Приклад механика и техн физика. 1997. Т. 38, № 1 . С 151-166 .

11 . Mavko, G. The rock physics handbook: tools for seismic analysis in porous media / G. Mavko, T. Mukerji, J. Dvorkin . Cambridge University Press, 1998

12 . Compendium of shock wave data / M . van Thiel (Ed .) // Lawrence Livermore Laboratory Report . 1977. UCRL-50108. P. 356-357.