CC BY
66
Мезоскопическое моделирование электрогидродииамических течений
Д.А. Медведев, А.Л. Куперштох
Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, 630090, Россия
Представлена мезоскопическая модель, позволяющая описать широкий круг электрогидродинамических явлений. Гидродинамические течения, перенос электрического заряда и распространение тепла за счет конвекции и диффузии моделируются с помощью метода решеточных уравнений Больцмана, так же как и фазовые переходы и взаимодействие несмешивающихся жидкостей. Приведены результаты моделирования динамики пузырьков и капель в электрическом поле, которые показывают перспективность предложенного подхода для многофазных и многокомпонентных систем.
Mesoscopic simulation of electrohydrodynamic flows
D.A. Medvedev and A.L. Kupershtokh M.A. Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia
The paper proposes a mesoscopic model that allows describing a wide range of electrohydrodynamic phenomena. Hydrodynamic flows, electric charge transfer and heat propagation due to convection and diffusion are simulated with the lattice Boltzmann equation method, like phase transitions and interaction of immiscible fluids. We present simulation results for the dynamics of bubbles and drops in the electric field. The results demonstrate that the proposed approach can be applied to multiphase and multicomponent systems.
1. Введение
К электрогидродинамическим явлениям относят широкий класс процессов, происходящих в жидкости, помещенной в электрическое поле в контакте с электродами [1]. На поверхности электродов могут происходить инжекция заряда и различные электрохимические реакции (последние могут идти также в объеме вещества). На заряженное вещество действует сила со стороны электрического поля, одновременно объемный электрический заряд изменяет пространственное распределение потенциала. Действие силы может приводить к возникновению течений жидкости.
В случае многокомпонентной среды (например, эмульсия капель воды в нефти) возникают дополнительные силы, действующие на поляризационные заряды, образующиеся вблизи контакта веществ с различной диэлектрической проницаемостью [2]. Эти силы вызывают движение и деформацию капель и могут приводить как к фрагментации, так и к коагуляции капель. Ускорение слияния капель в электрическом поле используется,
в частности, для очистки нефти от воды [3]. Поэтому адекватное описание этого процесса представляет собой важную научно-техническую проблему.
Кроме того, при прохождении электрического тока через проводящее вещество происходит выделение тепла. Повышение температуры и действие растягивающих напряжений на жидкость могут приводить к фазовому переходу и образованию паровых пузырьков.
Появившиеся пузырьки затем растягиваются электрическими силами вдоль направления электрического поля. При достаточно высокой напряженности поля возможно дробление пузырьков, а также их рост за счет продолжающихся испарения и дегазации жидкости. Кроме того, поскольку электрическая прочность газа значительно ниже, чем у жидкости, в пузырьках могут происходить импульсные микропробои, во время которых внутри пузырьков возникает проводящая фаза. Это приводит к накоплению зарядов вблизи поверхности пузырьков под действием электрического поля. Поляризованные таким образом пузырьки еще больше вытяги-
© Медведев Д-A., Куперштск A-Л., 2GG6
ваются вдоль поля. Кроме того, происходит усиление поля в жидкости у их полюсов, увеличивающееся с ростом удлинения, что может приводить к инициированию электрического пробоя непосредственно жидкого диэлектрика. Динамика пузырьков в электрическом поле теоретически и экспериментально исследована в работах [2, 4, 5].
Электрический пробой жидкостей происходит в несколько стадий. Вначале возле одного из электродов возникает одно или несколько светящихся плазменных образований. Время между подачей напряжения и появлением таких образований называют статистическим временем запаздывания, его величина меняется от опыта к опыту. В зависимости от условий возможны различные механизмы зарождения пробоя: пузырьковый, тепловой, ионизационный и др. [6-13].
На следующей стадии из появившихся образований происходит рост проводящих каналов — стримеров — по направлению к противоположному электроду, которые могут также ветвиться. Картина стримерной структуры носит принципиально стохастический характер, то есть не воспроизводится в деталях от опыта к опыту [14-20]. Внутри стримеров электропроводность достаточно высока, поэтому в них происходит выделение энергии, что приводит к повышению температуры и расширению. Плотность вещества в стримерах, напротив, понижена [18]. Вокруг каждого расширяющегося канала возникают течения жидкости с расходящимися волнами сжатия. Происходит интерференция этих волн, что дополнительно осложняет картину течения. Возможно также обратное влияние волн сжатия на динамику стри-мерных каналов.
Моделирование электрогидродинамических течений представляет собой сложную задачу, в которой необходим учет многих явлений, происходящих одновременно. В данной работе представлена мезоскопическая модель, основанная на методе решеточного уравнения Больцмана (lattice Boltzmann equation, LBE), предложенная и развиваемая авторами [21-29].
2. Основные уравнения
Движение жидкости описывается уравнением неразрывности Эр
dt
+ divpu = 0
(1)
и уравнением Навье-Стокса
dpu
dt + diV П “p
= F + nV u +
z+-
graddiv u.
(2)
Здесь p — плотность жидкости; u — скорость течения;
п
$ = + puaue — невязкая часть тензора напря-
жений.
Уравнения для изменения плотности каждого вида носителей электрического заряда выглядят следующим образом Эп,-
dt
+ div niu =
= Di V ni - div
* " Abn E qi
+ si - ri .
(З)
Здесь Ь, — эффективная подвижность носителей заряда qi в электрическом поле; Di — их коэффициенты диффузии; si — скорости рождения носителей заряда за счет ионизации; г, — скорости их гибели за счет рекомбинации.
В случае равенства коэффициентов диффузии Д = = D можно умножить уравнения (3) на и просуммировать по i, что приводит к уравнению для суммарной плотности электрического заряда дq
dt
+ div qu = DV q - div aE.
(4)
Здесь q = £ д,п, — суммарная плотность свободных зарядов; ст = ХЬ,^,\п, ---локальная электропровод-
ность, зависящая от концентраций носителей заряда.
Распределение электрического потенциала находится из уравнения Пуассона
Шу(еУф) = -4я#, (5)
где е — диэлектрическая проницаемость, которая зависит от плотности вещества и может быть неоднородной по пространству.
Объемная сила, действующая на заряженную диэлектрическую жидкость в электрическом поле, выражается обычной формулой
F2 І
F = qE-—Vє +—V 8п 8п
2 * дє " F 2p — dp
(6)
Второе и третье слагаемые описывают действие электрического поля на поляризационные заряды в неоднородном диэлектрике и электрострикционные силы. В частности, для газов обычно е = 1 + р/р0. В этом слу-
= е -1 и формула (6) принимает вид (при
чае p
дє
dp
&т
отсутствии свободный зарядов)
F =— VF 2.
8п
(7)
Возможны и другие модели для зависимости е(р). Например, для неполярных жидкостей в первом приближении справедлива формула Клаузиуса-Моссотти [30]
є = І +
pp
І -pp/З’
где в = 4па/ т; а — поляризуемость; т — молекулярная масса. В этом случае имеем
p
дє
dp
v J
(є- І)(є + 2)
Для полярных жидкостей справедлива формула Он-загера-Кирквуда-Фрелиха [31]
Аер = (е-е^ )(2е+ е^),
где
2 4п * е^ + 2
А = gм4 ,
*ткТ) 3
Здесь — дипольный момент молекул в сильно разре-
женном газе; g — средний фактор корреляции, характеризующий ближний ориентационный порядок; е^ = = п2 — диэлектрическая проницаемость, связанная с деформационной (упругой) поляризацией; п — показатель преломления. При этом имеем
Р
. др
V у
где
В = -
= (е-е^ )В,
2у +1
4у-1- (У- 1)(2у +1)/У
, У = е/е^.
Интересно, что В ~ 1 (максимальное значение В ~ = 1.11 при у = 2.22). В работе [32] использовалось значение В =1.
3. Метод решения
Для решения системы уравнений (1)-(6) использовался метод расщепления по физическим процессам [33]. Полный шаг по времени разбивается на несколько подшагов, соответствующих различным физическим явлениям, подшаги выполняются один за другим. В задачах электрогидродинамики можно выделить следующие различные процессы:
1) гидродинамические течения,
2) конвективный перенос и диффузия носителей заряда,
3) вычисление электрического потенциала и переноса заряда токами проводимости,
4) действие объемных электростатических сил на жидкость в электрическом поле,
5) выделение тепла при прохождении электрического тока через проводящую жидкость,
6) фазовые переходы или взаимодействие несмеши-вающихся жидкостей.
3.1. Моделирование гидродинамических течений
Для моделирования гидродинамических течений применялся метод LBE [34-36], который является принципиально мезоскопическим. В отличие от классических методов расчета течений жидкости путем решения уравнений Навье-Стокса метод решеточных уравнений Больцмана рассматривает течение как движение ансамбля псевдочастиц, имеющих некоторую функцию распределения по дискретным скоростям. Метод LBE широко используется для моделирования течений жидкости, включая многофазные и многокомпонентные, благодаря его численной устойчивости и легкости реа-
лизации сложных граничных условий. Теоретическим обоснованием этого метода можно считать тот факт, что во втором порядке разложения Чепмена-Энскога из уравнений LBE получаются уравнения гидродинамики, то есть уравнения неразрывности и Навье-Стокса. Кроме того, показано, что решеточное уравнение Больцмана может быть выведено из непрерывного уравнения Больцмана с помощью подходящей дискретизации в импульсном и координатном пространстве [37, 38].
В методе LBE в качестве основных переменных используются одночастичные функции распределения Ык, определенные в узлах равномерной пространственной сетки для всего конечного набора возможных скоростей частиц с к. Уравнения эволюции имеют вид Ык (х + ск Д*, t + Д*) =
= Ык (х, 0 +
(р(х, t), и(х, t)) - Ык (х, t)
+ ДЫк, (8)
где Дг — шаг по времени (вектора сетки ек = с к Д*); ДЫк — изменение функций распределения за счет действия объемных сил. Плотность вещества в узле и скорость могут быть вычислены по формулам
р = ХЫк, и = ЕЫк ск/р•
кк
Для изотермических вариантов LBE-моделей жидкости обычно используется разложение равновесных функций распределения в ряд по скорости и до второго порядка, которое имеет вид
ЫГ (р, и) =
= рик
1++_(ски)
е
262
,(и)2
2е
(9)
Здесь е = кТ/т — нормированная температура. Набор векторов скорости ск и коэффициенты ч/к зависят от конкретной сетки. Вместе с тем, в случае действия объемных сил в уравнениях (8) и (9) следует использовать скорость и' = £ Ык ск/р. При этом реальная скок
рость течения, входящая в уравнения (1)-(4), выражается как и = и' + Ди/2, где Ди = Г Д */р — изменение скорости жидкости за счет действия силы Р в течение шага по времени Дг [39-41].
В двумерной модели D2Q9 используется квадратная
сетка, возможные значения скорости частиц ск
= 0,
к/Д*, Ли/Аt, где h — шаг сетки (рис. 1). Весовые коэффициенты равны w0 = 4/9, и>1-4 = 1/9, ^5-8 = 1/36. Безразмерное время релаксации т определяет коэффициент кинематической вязкости V = еД*(т-1/2). Для ряда изотермических вариантов LBE-моделей: D1Q3, D2Q9, D3Q15, D3Q19 нормированная температура равна е = (И/ Д )2/з.
Для учета действия объемных сил на жидкость в [4244] был разработан метод точной разности для LBE:
ДЫк = Ы?ф, и' + Ди) - Ы?ф, и'). (10)
т
в 2 5
3 V 1 1
уГ 1
? 4 Б
Рис. 1. Возможные скорости частиц для модели Б2р9
3.2. Конвективный перенос и диффузия носителей заряда
Уравнения для изменения плотности носителей заряда, учитывающие их перенос за счет течения жидкости, диффузию, а также процессы ионизации и рекомбинации:
дп
—L + Шу п, и = Д- V 2 п, + si - г,
д* г г г г г
решаются с помощью метода дополнительного LBE-компонента с нулевой массой (пассивный скаляр) [24], подобного методу, предложенному в [45].
Уравнения эволюции для функций распределения Qы, соответствующих каждому типу носителей заряда , записываются следующим образом:
Йй(х + скД*, * + Д*) = Йй(х> *) +
0^(п,(х *X u(x, *))- (х. *)
т, '
(11)
Равновесные функции распределения 0^ (п,, и) зависят от локальной концентрации носителей заряда данного типа п, = £ 0М
0м(п, ,и) = п^к
1 + Ск" + (ски) е 2е2
и
2е
Коэффициенты диффузии носителей заряда разных типов Д, = еД*(т, - 1/ 2) могут задаваться независимо друг от друга с помощью изменения времен релаксации т,.
Точные значения скоростей ионизации si и рекомбинации г, в жидкостях неизвестны, однако некоторые обсуждения и приближенные зависимости для слабо-проводящих жидкостей приведены в [46].
3.3. Вычисление электрического потенциала и переноса заряда токами проводимости
Уравнение Пуассона (5) для потенциала электрического поля решается совместно с уравнениями переноса
носителей электрического заряда, связанного с их подвижностью в электрическом поле
* "
дп,
д*
+ &у
-Г-Ащ Е
м
= 0.
Неявные по времени конечно-разностные уравнения для уравнений переноса носителей заряда и для уравнения Пуассона записываются в виде
* "
пп+1 = пп + Д* diy
\ч\
7 И+^^,_И+1
Ь,п VФ
(12)
Шу (еVфя+1) = -4п£ д,п”+1.
Полученную систему уравнений можно решать методом итераций относительно значений концентраций всех носителей заряда п"+1 и значений электрического потенциала фп+1 на новом шаге по времени для каждого узла сетки.
В простейшем случае, когда справедливо уравнение (4), для учета токов проводимости в методе расщепления достаточно использовать уравнение
— = -&у аЕ.
д*
При использовании неявной по времени конечноразностной схемы для этого уравнения значение полного заряда в узле на следующем временном слое может быть прямо подставлено в конечно-разностное уравнение Пуассона. В этом случае система уравнений (12) сводится к разностным уравнениям только для значений
п+1
потенциала на следующем временном слое ф, ^ , которые решаются итеративно методом релаксации [47] при соответствующих граничных условиях.
Объемная сила, действующая на жидкий диэлектрик в электрическом поле, вычисляется по формуле (6). Декартовы компоненты силы находятся с помощью аппроксимации производных центральными разностями. Например, для проекции электрического поля на ось х имеем (Ех ),,] =- (ф,+1,] -ф,_1;] )/2к. Для учета действия этих сил в LBE используется метод точной разности (10).
3.4. Фазовые переходы
Фазовые переходы жидкость - газ моделируются с помощью метода, предложенного в работах [48, 49]. Чтобы описать фазовые переходы, в этой модели вводятся силы притяжения, действующие между каждой парой соседних узлов решетки. Вектор результирующей силы имеет вид
Г(х) = ¥(р(х))£ Gk ¥(р(х + е к ))е к. (13)
к
Коэффициенты Gk > 0 различны для основных и диагональных направлений сетки; ^(р) — возрастающая функция плотности («эффективная масса»). Для двумерной модели D2Q9 используются коэффициенты G0 и G1 соответственно для основных и диагональных на-
правлений решетки. Для обеспечения изотропии взаимодействия необходимо соотношение G1 = G0/4• Уравнение состояния для этой изотермической модели D2Q9 имеет вид
Р = р0-3 Go^ 2(р).
(14)
Для некоторого вида функций ^(р) уравнение состояния (14) допускает фазовые переходы. В частности, нами использовалась функция, предложенная в работе [48],
¥(Р) = Р о (1 exp(p/po)). (15)
Критическая точка находится обычным образом (Эр/Эр = 0, д2р/Эр2 = 0) и для описанной модели соответствует G0* = 40/(3р0) и р*=р01п2. При р0 = 1 имеем для модели D2Q9: р* = ln2 = 0.693 и G0* = 4/9. При значениях G0 > G0* возможно сосуществование плотной (жидкость) и разреженной (пар) фаз вещества.
В данном методе граница между газом и жидкостью представляет собой тонкий переходный слой (шириной в несколько узлов сетки), в котором плотность монотонно меняется между равновесными значениями для разных фаз. Таким же образом описываются межфазные границы в модели фазового поля (phase-field method). Взаимодействие (13) обеспечивает поверхностное натяжение на границе жидкость - газ, которое зависит от величины G0 [49]. Величина коэффициента поверхностного натяжения была измерена с помощью расчета стационарного состояния круглой капли радиуса R, окруженной насыщенным паром. График разности давлений внутри и вне капли в зависимости от обратного радиуса 1/R представляет собой прямую линию (закон Лапласа pin - pout = X/R). Для значения G0 = 0.5, используемого в последующих расчетах, коэффициент поверхностного натяжения X = 0.019 в безразмерных единицах (при р0 = 1, т = 1). При уменьшении G0 коэффициент поверхностного натяжения уменьшается и обращается в ноль в критической точке, как и должно быть. При значениях G0 < G0* разделения фаз не происходит.
3.5. Взаимодействие несмешивающихся жидкостей
Взаимодействие двух несмешивающихся жидкостей также моделировалось с помощью метода, предложенного в работе [48]. В простейшем случае вводится взаимодействие между соседними узлами в форме
F(x) = V(Ps (x))£ £ GksX(x + ek ))ek •
X k
Индексы s и X соответствуют различным компонентам. В случае двух жидкостей возможные значения индексов 1 и 2. Величина р s = £ Nks — плотность ком-
k
понента s в данном узле. Здесь Nks — одночастичные функции распределения для данного компонента.
Суммарная плотность вещества в узле р = £р,,
суммарный импульс ри = £р,иs, где скорость каждого компонента иs = £ Nks ck/рs. Взаимодействие меж-
ду компонентами приводит к изменению их скоростей Ди* = ^ Д*/р*, которое используется в операторе столкновений для каждого компонента ^к? [48]. В наших вычислениях использовалось одно значение времени релаксации т для всех компонентов, т.е. вязкости жидкостей были одинаковыми.
Для моделирования динамики несмешивающихся жидкостей значения коэффициентов задавались как
= 0 GЪx = > °, причем 4 для
обеспечения изотропности взаимодействия. Функция эффективной массы выбиралась в простейшем виде ¥(р) = р.
4. Результаты
4.1. Динамика проводящих пузырьков в электрическом поле
Электрическая прочность газов гораздо меньше, чем у жидкостей, поэтому существующие и образующиеся в жидкости пузырьки газа и пара могут пробиваться при достижении определенного критического размера. После микроразрядов внутри пузырька возникает проводящая фаза (плазма). На рис. 2 представлены результаты моделирования динамики деформации и дробления проводящих пузырьков пара, изначально существующих в жидком диэлектрике, под действием электрического по-
ИГ
*
I
11
0 100 200 300 400 500 600 700
t
ИГ\
I
I
I
11
100 200 300 400 500 600 700 770 850
t
ИЛ
I
I
111
100 200 300 400 500 600 700 800 850
t
Рис. 2. Деформация и дробление парового пузырька в электрическом поле при различной электропроводности внутри пузырька а = 0.5 (a), 0.375 (б), 0.2 (в). Размер расчетной области 65x250
ля. Продемонстрировано несколько режимов их деформации, роста и дробления.
Наблюдается вытягивание пузырька вдоль направления приложенного электрического поля. При этом происходит его сжатие в поперечном направлении в районе экватора, так как здесь нет растягивающих электрических сил. При высокой проводимости газа в пузырьке в экваториальной части образуется перетяжка, и пузырек дробится на два меньших (рис. 2, а). При меньшей электропроводности газа внутри пузырька деформация протекает медленнее, и происходит отделение двух пузырьков с полюсов первоначального. В этом случае центральный пузырек не несет заряда и быстро исчезает путем конденсации пара (рис. 2, б). В случае низкой проводимости внешнее гидростатическое давление преобладает, и пузырек после первоначального удлинения исчезает совсем (рис. 2, в).
Похожие процессы наблюдались при пробое жидких диэлектриков [50]. При незавершенном пробое каналы стримеров распадались на цепочки пузырьков, которые затем быстро исчезали.
На рис. 3 представлены результаты моделирования коалесценции двух пузырьков. Диэлектрическая проницаемость жидкости принималась е = 1. Идеальная проводимость пузырьков моделировалась в приближении е = 100. После включения электрического поля пузырьки поляризуются и начинают притягиваться друг к другу.
\а]
• І • 1 ( в 1 &
• • • • V 1 1 9
0 600 920 960 1000 t 1100 1300 2600
б
И Iі і'
0 400 1020 1080 1120 1180 1400 1600
t
Рис. 3. Взаимодействие двух проводящих пузырьков в жидком диэлектрике из-за поляризации в электрическом поле. Пузы1рьки расположены вдоль электрического поля (а). Относительное расположение по отношению к электрическому полю — диагональное (б). Расчет на сетке 400x400. Показана центральная часть расчетной области 100x400
При их приближении на некоторое малое расстояние, тонкая перемычка жидкости, которая вытесняется из зазора между пузырьками, имеет радиальную скорость. Затем при инерционном движении плотность вещества в перемычке на оси уменьшается (т.е. происходит «кавитация») и перемычка исчезает. После этого жидкость, которая была вытеснена из зазора между пузырьками, продолжает по инерции двигаться, и это течение растягивает пузырек в «экваториальной плоскости». При слиянии пузырьков заряды в месте касания взаимно уничтожаются и далее происходят процессы, аналогичные показанным на рис. 2.
4.2. Капли в электрическом поле
На капли, диэлектрическая проницаемость которых е1 больше проницаемости окружающей жидкости е, действует во внешнем однородном поле пара сил, растягивающих капли вдоль поля. Эти силы пропорциональны Еа (1 - е/е!) и возрастают с увеличением отношения е^ е (Еа — средняя напряженность поля). В пределе е1/ е^го можно моделировать динамику идеально проводящих капель.
Поведение капель с различной проницаемостью е1, находящихся в жидкости с е = 1, в электрическом поле показано на рис. 4 и 5. В электрическом поле капля по-
| а є = 1 • є. = 100 б • \± 1 І£ I
д I І£ л ЧІ0 ^ 1* и О 13 О
Рис. 4. Динамика капли с диэлектрической проницаемостью, отличной от проницаемости основной жидкости, е! = 100, Еа = 0.035. Время t = 0 (а), 100 (б), 200 (в), 300 (г), 400 (д), 500 (е), 600 (ж) и 700 (з)
| а є = 1 • є. — 20 б • \± I к. т ш
І_д_ о О о іе Ч¥ 6*) и * 71 % $ 1_і * «
Рис. 5. Динамика капли с диэлектрической проницаемостью, отличной от проницаемости основной жидкости, е! = 20, Еа = 0.1. Время
I = 0 (а), 100 (б), 200 (в), 300 (г), 400 (д), 500 (е), 600 (ж) и 700 (з)
Рис. 6. Коалесценция двух проводящих капель в жидкости из-за их поляризации в электрическом поле. Среднее электрическое поле Еа = 0.3. Расчет на сетке 400x400. Показана центральная область 75x200
ляризуется и начинает растягиваться электрическими силами. Возникает течение, представляющее собой два вихревых диполя, двигающихся в противоположных направлениях (на рисунке — вверх и вниз). Такого типа течения наблюдаются в стратифицированной жидкости [51]. В случае, представленном на рис. 5, вихри более выражены, так как значение Е2 (1 - ^е1) больше.
На рис. 6 представлены результаты моделирования коалесценции двух «проводящих» капель (капель с диэлектрической проницаемостью е1 = 100) в жидкости с е = 1. Процесс электрокоалесценции интересен в связи с его использованием на практике, например, для ускорения коалесценции капель воды в нефти [3].
После включения электрического поля капли поляризуются и начинают притягиваться друг к другу. При их приближении на некоторое малое расстояние возникает тонкая перемычка жидкости, которая вытесняется из зазора между каплями. После касания жидкость в каплях продолжает по инерции двигаться навстречу в вертикальном направлении и в месте касания возникает течение, которое растягивает каплю в «экваториальной плоскости». Далее результирующая капля вытягивается, и возможны процессы, аналогичные показанным на рис. 4 и 5.
4.3. Электрострикция
Электропроводность вещества полагалась равной нулю, а диэлектрическая проницаемость — зависящей от плотности как е = 1 + р/р0. Рассмотрена эволюция парового пузырька в жидкости. В этом случае определяющими являются электрострикционные силы, т.е. действующая на вещество объемная сила выражается формулой (7). Результаты моделирования представлены на рис. 7, размер расчетной области 200x200. Электроды находятся сверху и снизу, граничные условия по горизонтали периодические. Темный цвет соответствует меньшей плотности. После подачи напряжения пузырек начинает вытягиваться и затем рвется на два. В целом, картина подобна случаю с проводящим пузырьком (рис. 2, а).
При уменьшении напряженности поля удлинение пузырька также уменьшается. Существует пороговое
поле, ниже которого фрагментации не происходит, что качественно согласуется с результатами [4, 5].
4.4. Тепловой пробой жидкости
Тепловой пробой жидкого диэлектрика возможен, когда его электропроводность сильно возрастает с ростом температуры, например, за счет ионизации. В этом случае возникшая возле одного из электродов (например, за счет полевой эмиссии или микроразряда возле острия) проводящая область нагревается электрическим током, вызывая переход окружающей жидкости в проводящее состояние. Нагретая область с повышенной концентрацией носителей заряда продолжает расти за счет продолжающегося энерговыделения и вытягивается к противоположному электроду под действием электрических сил. Если скорость нагрева превышает отвод тепла за счет теплопроводности, через некоторое время растущий канал достигает противоположного электрода, и наступает канальная стадия электрического пробоя жидкости.
В качестве первого приближения рассмотрим простейшую «toy model», в которой плотность вещества не зависит от температуры (а также нет фазовых переходов). Уравнение для изменения температуры за счет конвекции, теплопроводности и джоулева нагрева имеет вид
— + divTu = XV2T + П(Т)E .
dt cp
Здесь x — коэффициент температуропроводности, а слагаемое ст(Т)E2jcp представляет собой повышение температуры за счет энерговыделения при прохождении электрического тока через проводящее вещество. Электропроводность ст(Г) представляет собой быстрорастущую функцию от температуры. В наших расчетах при Т > Т0 + Tth использовалась зависимость ст = ст0 exp(- A (Т - Т0)), а при меньших температурах электропроводность считалась нулевой.
Изменение температуры за счет конвекции и диффузии моделировалось так же, как и перенос электричес-
а | • б • в 1 и
д 1 е ж 9 1 и § $
Рис. 7. Деформация пузырька в электрическом поле за счет электро-стрикции. Среднее электрическое поле Еа = 1.2, время г = 60 (а), 180 (б), 300 (в), 400 (г), 500 (д), 600 (е), 700 (ж) и 800 (з)
0 0 0 0
д 0 циж ▼
Рис. 8. Тепловой пробой жидкости. Время г = 100 (а), 200 (б), 300 (в), 400 (г), 500 (д), 600 (е), 700 (ж) и 800 (з)
кого заряда с помощью метода дополнительного LBE-компонента. Коэффициент температуропроводности X = 0Аг(тТ -1/2) регулировался с помощью изменения времени релаксации ТТ.
Результаты представлены на рис. 8, размер расчетной области 251x201. Граничные условия по горизонтали периодические, сверху и снизу расположены электроды. В середине нижнего электрода находится острие, с которого возможна инжекция заряда. Средняя напряженность поля Еа = 0.04, ст0 = 0.15, коэффициент А = = 0.1, теплоемкость ср = 0.001, начальная температура Т0 = 1, пороговая температура Г(Ь = 0.02. На рисунке темный цвет соответствует более высокой температуре. После подачи напряжения возле острия возникает нагретая область, которая затем растет из-за действия электрических сил и продолжающегося энерговыделения. Электрические силы приводят к возникновению течения в виде вихревого диполя. Нагретая область представляет собой проводящий канал, который растет в сторону противоположного электрода. После достижения каналом противоположного электрода через него начинает протекать ток, что приводит к сильному разогреву вещества и расширению канала.
При уменьшении напряженности поля развитие пробоя сильно замедляется. Существует пороговая напряженность поля, ниже которой токопроводящий канал не развивается, для использованных значений параметров она равна Еа* = 0.037.
Более сложная модель, учитывающая фазовые переходы для уравнений состояния общего вида, а также зависимость электропроводности и диэлектрической проницаемости от плотности вещества, разрабатывается в настоящее время.
5. Заключение
Представлена мезоскопическая модель, позволяющая описать широкий круг электрогидродинамических явлений — от возникновения электрогидродинамичес-ких течений до динамики образования и роста пузырьков, а также начальных стадий электрического пробоя диэлектриков. Гидродинамические течения, перенос электрического заряда и распространение тепла за счет конвекции и диффузии моделируются с помощью мето-
да LBE, так же как и фазовые переходы и взаимодействие несмешивающихся жидкостей. Изменение электрического потенциала и токи проводимости вычисляются с помощью конечно-разностных методов.
Приведенные примеры расчетов демонстрируют широкие возможности применения предложенного метода, который особенно полезен для моделирования задач со свободными границами (многофазные и многокомпонентные системы).
Работа выполнена при частичной финансовой поддержке РФФИ (гранты №№ 03-02-16474 и 05-0833715).
Литература
1. Остроумов Г.А. Взаимодействие электрических и гидродинамических полей. Физические основы электрогидродинамики. - М.: Наука, 1979. - 319 с.
2. Beroual A. Behaviour of charged and uncharged bubbles in dielectric liquids subjected to electric stress // J. Appl. Phys. - 1992. - V. 71. -No. 3. - P. 1142-1145.
3. Lundgaard L.E., Berg G., Pedersen A., Nilsen PJ. Electrocoalescence of water drop pairs in oil // Proc. 14th Int. Conf on Dielectric Liquids, IEEE No. 02CH37319. - Graz, Austria, 2002. - P. 215-219.
4. Garton C.G., Krasucki Z. Bubbles in insulating liquids: stability in an electric field // Proc. of Royal Society. - 1964. - V. A280. - No. 1381. -P. 211-226.
5. Ogata S., Tan K., Nishijima K., Chang J.-S. Development of improved bubble disruption and dispersion technique by an applied electric field method // American Institute of Chemical Engineers Journal. - 1985. -V. 31. - No. 1. - P. 62-69.
6. Арсентьев B.B. К теории импульсных разрядов в жидкой среде //
Прикладная механика и техническая физика. - 1965. - № 5. -С. 51-57.
7. Krasucki Z. Breakdown of liquid dielectrics // Proc. of Royal Society. -
1966. - V. A294. - No. 1438. - P. 393^04.
8. Halpern B., Gomer R. Field emission in liquids // J. Chem. Phys. -1969. - V. 51. - No. 3. - P. 1031-1047.
9. Watson P.K., Chadband W.G., Sadeghzadeh-Araghi M. The role of electrostatic and hydrodynamic forces in the negative-point breakdown of liquid dielectrics // IEEE Transactions on Electrical Insulation. -1991. - V. 26. - No. 4. - P. 543-559.
10. Sadeghzadeh-Araghi M., Qureshi M.I., Chadband W.G., Watson P.K. Measurement of the growth of cavities and of EHD instabilities during the negative-point breakdown of silicon fluids // IEEE Transactions on Electrical Insulation. - 1991. - V. 26. - No. 4. - P. 663-672.
11. Jones H.M., Kunhardt E.E. Development of pulsed dielectric breakdown in liquids // J. Phys. D: Appl. Phys. - 1995. - V. 28. - No. 1. -P. 178-188.
12. Sinkevich O.A., Smirnov P. V Heterogeneous mechanism of electrical breakdown of liquid dielectrics // J. of Moscow Physical Society. -1996. - V. 6. - P. 101.
13. Коробейников С.М. О роли пузырьков в электрическом пробое жидкостей. Сопоставление с экспериментом // Теплофизика высоких температур. - 1998. - Т. 36. - № 4. - С. 541-547.
14. Алхимов А.П., Воробьев B.B., Климкин В.Ф., Пономаренко А.Г., Солоухин РИ. О развитии электрического разряда в воде // Докл. АН СССР. - 1970. - Т. 194. - № 5. - С. 1052-1054.
15. Яншин Э.В., ОвчинниковИ.Т., ВершининЮ.Н. Оптические исследования предпробойных явлений в воде в наносекундном диапазоне // ЖТФ. - 1973. - Т. 43. - № 10. - С. 2067-2074.
16. Gournay P., Lesaint O. Study of the dynamics of gaseous positive streamer filaments in pentane // Proc. 11th Int. Conf. on Conduction and Breakdown in Dielectric Liquids, IEEE No. 93CH3204-5. - Ba-den-Dattwil, Switzerland, 1993. - P. 289-293.
17. Gavrilov I.M., Kukhta VR., Lopatin V.V., Petrov P.G. Dynamics of prebreakdown phenomena in water // IEEE Transactions on Dielectrics and Electrical Insulation. - 1994. - V. 1. - No. 3. - P. 496-502.
18. Gournay P., Lesaint O. On the gaseous nature of positive filamentary streamers in hydrocarbon liquids. II: Propagation, growth and collapse of gaseous filaments in pentane // J. Phys. D: Appl. Phys. - 1994. -V. 27. - No. 10. - P. 2117-2127.
19. Кухта В.Р., Лопатин В.В., Носков М.Д. Фрактальные характеристики приэлектродных образований при электрическом разряде в воде // Изв. вузов. Физика. - 1994. - № 7. - С. 89-92.
20. Кухта В.Р., Лопатин В.В., Носков М.Д. Применение фрактальной модели к описанию развития разряда в конденсированных диэлектриках // ЖТФ. - 1995. - Т. 65. - № 2. - С. 63-75.
21. Kupershtokh A.L., Medvedev D.A. Simulation of gas-dynamic flows during streamer propagation at liquid dielectrics breakdown // Conference record 1998 IEEE Int. Symposium on Electrical Insulation. -Arlington, VA, 1998. - P. 611-614.
22. Куперштох А.Л., Медведев Д.А. Моделирование газодинамических течений, сопровождающих распространение стримеров при пробое жидких диэлектриков // Современные проблемы электрофизики и электрогидродинамики жидкостей: Сборник докладов V Международной научной конференции. - Санкт-Петербург,
1998. - С. 126-130.
23. Kupershtokh A.L., Medvedev D.A. Simulation of hydrodynamic flows during streamer propagation in dielectric liquids // Proc. 13th Int. Conference on Dielectric Liquids, IEEE No. 99CH36213. - Nara, Japan,
1999. - P. 179-182.
24. Medvedev D.A., Kupershtokh A.L. Use of the lattice Boltzmann equation method to simulate charge transfer and electrohydrodynamic phenomena in dielectric liquids // Proc. 2nd Int. Workshop on Electrical Conduction, Convection and Breakdown in Fluids. - Grenoble, France,
2000. - P. 60-63.
25. Куперштох А.Л., Медведев Д.А. Структура и динамика «плазменных» каналов при пробое жидких диэлектриков // Динамика сплошной среды (Математические проблемы механики сплошных сред). - Вып. 116. - Новосибирск: Институт гидродинамики, 2000. - С. 137-141.
26. Куперштох А.Л., Медведев Д.А. Метод решеточного уравнения Больцмана в задачах электрогидродинамики // Динамика сплошной среды (Математические проблемы механики сплошных сред). - Вып. 118. - Новосибирск: Институт гидродинамики,
2001.- С. 117-121.
27. Medvedev D.A., Kupershtokh A.L. Modeling of electrohydrodynamic flows and micro-bubbles generation in dielectric liquid by lattice Boltzmann equation method // Proc. 14th Int. Conference on Dielectric Liquids, IEEE No. 02CH37319. - Graz, Austria, 2002. - P. 45-48.
28. Медведев Д.А., Ершов А.П., Куперштох А.Л. Численное исследование гидродинамических и электрогидродинамических неустойчивостей // Динамика сплошной среды (Математические проблемы механики сплошных сред). - Вып. 120. - Новосибирск: Институт гидродинамики, 2002. - С. 93-103.
29. Kupershtokh A.L., Medvedev D.A. Lattice Boltzmann equation method in electrohydrodynamic problems // Proc. 5th International EHD Workshop. - Poitiers, France, 2004. - P. 241-246.
30. Фейнман Р., Лейтон Р., Сэндс М. Фейнмановские лекции по физике. Т. 5. Электричество и магнетизм. - М.: Мир, 1977. - 300 с.
31. Шахпоронов М.И. Методы исследования теплового движения молекул и строения жидкостей. - М.: Изд-во Московского университета, 1963. - 281 с.
32. Коробейников С. М. Пузырьковая модель зажигания импульсного электрического разряда в жидкостях // Дисс. ... докт. физ.-мат. наук. - Новосибирск, 1997. - 310 с.
33. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. - Новосибирск: Наука, 1967. - 197 с.
34. McNamara G., Zanetti G. Use of the Boltzmann equation to simulate lattice-gas automata // Phys. Rev. Let. - 1988. - V. 61. - No. 20. -P. 2332-2335.
35. HigueraF.J., Jimunez J. Boltzmann approach to lattice gas simulation // Europhysics Letters. - 1989. - V. 9. - No. 7. - P. 663-668.
36. Qian Y.H., d’Hummres D., Lallemand P. Lattice BGK models for Navier-Stokes equation // Europhysics Letters. - 1992. - V. 17. -No. 6. - P. 479^84.
37. Abe T. Derivation of the lattice Boltzmann method by means of the discrete ordinate method for the Boltzmann equation // Journal of Computational Physics. - 1997. - V. 131. - No. 1. - P. 241-246.
38. He X., Luo L.-S. A priori derivation of the lattice Boltzmann equation // Physical Review E. - 1997. - V. 55. - No. 6. - P. 6333-6336.
39. Ginzburg I., Adler P.M. Boundary flow condition analysis for the three-dimensional lattice Boltzmann model // J. Phys. II France. -1994. - v. 4. - No. 2. - P. 191-214.
40. Shan X., Doolen G. Multicomponent lattice-Boltzmann model with interparticle interaction // J. of Statistical Physics. - 1995. - V. 81. -Nos. 1/2. - P. 379-393.
41. Ladd A.J.C., Verberg R. Lattice-Boltzmann simulations of particle-fluid suspension // J. of Statistical Physics. - 2001. - V. 104. - Nos. 5/6. -P. 1191-1251.
42. Kupershtokh A.L. Calculations of the action of electric forces in the lattice Boltzmann equation method using the difference of equilibrium distributin functions // Современные проблемы электрофизики и электрогидродинамики жидкостей: Сборник докладов VII Международной научной конференции. - Санкт-Петербург, 2003. -С. 152-155.
43. Куперштох А.Л. Учет действия объемных сил в решеточных уравнениях Больцмана // Вестник НГУ: Серия «Математика, механика и информатика». - 2004. - Т. 4. - № 2. - С. 75-96.
44. Kupershtokh A.L. New method of incorporating a body force term into the lattice Boltzmann equation // Proc. 5th International EHD Workshop. - Poitiers, France, 2004. - P. 241-246.
45. Shan X. Simulation of Rayleigh-Benard convection using a lattice Boltzmann method // Phys. Rev. E. - 1997. - V. 55. - No. 3. - P. 27802788.
46. Apfelbaum M.S. The prebreakdown EHD equations for liquid insulators // Современные проблемы электрофизики и электрогидродинамики жидкостей: Сборник докладов VII Международной научной конференции. - Санкт-Петербург, 2003. - С. 9-13.
47. Karpov D.I., Kupershtokh A.L. Models of streamer growth with «phy-sical» time and fractal characteristics of streamer structures // Conference record 1998 IEEE Int. Symposium on Electrical Insulation. -Arlington, VA, 1998. - P. 607-610.
48. Shan X., Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components // Phys. Rev. E. - 1993. - V. 47. -No. 3. - P. 1815-1819.
49. Shan X., Chen H. Simulation of nonideal gases and liquid-gas transitions by the lattice Boltzmann equation // Phys. Rev. E. - 1994. -V. 49. - No. 4. - P. 2941-2948.
50. Costeanu L., Lesaint O. On mechanisms involved in the propagation of subsonic positive streamers in cyclohexane // Proc. 14th Int. Conf. on Dielectric Liquids, IEEE No. 02CH37319. - Graz, Austria, 2002. -P. 143-146.
51. АфанасьевЯ.Д., Воропаев С.И. Горизонтальная затопленная струя в стратифицированной жидкости // Механика жидкости и газа. -1993. - № 6. - С. 10-16.
