Численное моделирование деформации свободновисящих пленок под действием тепловой нагрузки Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 12, № 4, 2007

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДЕФОРМАЦИИ СВОБОДНОВИСЯЩИХ ПЛЕНОК ПОД ДЕЙСТВИЕМ ТЕПЛОВОЙ

НАГРУЗКИ*

А. С. Овчарова Институт гидродинамики им. М.А. Лаврентьева СО РАН,

Новосибирск, Россия e-mail: [email protected]

We consider a deformation of a thin film, which is hanging between two solid flat walls under the action of a thermal load. A two-dimensional model is applied to describe the motion of thin layers of viscous nonisothermal liquid under microgravity conditions. The model is based on the Navier—Stokes equations. We give a detailed description of the numerical algorithm for solution of this problem. Numerical analysis of the influence of thermal loads on the deformation and breaking of freely hanging thin films has been carried out. The mutual influence of capillary and thermo-capillary forces on thin film free surface position has been shown. Solutions to some model problems are presented.

Введение

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

Подавляющее число публикаций связано с течением тонких пленок по твердой подложке и лишь небольшая часть посвящена исследованиям процессов, возникающих в свободновисящих пленках. В работе [1] исследуется модель деформации свободной невесомой пленки жидкости, закрепленной по плоскому контуру и подверженной действию термокапиллярных сил, в длинноволновом приближении. Уравнение, описывающее эволюцию толщины пленки, имеет второй порядок по времени и четвертый порядок по продольным координатам. В настоящей работе модель деформации и разрыва свобод-новисящей пленки основана на уравнениях Навье—Стокса.

* Работа выполнена при финансовой поддержке Совета программы поддержки ведущих научных школ (проект НШ-5873.2006.1), Междисциплинарного интеграционного проекта СО РАН (№ 111).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

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

1. Математическая модель

Плоская пленка с плотностью р, кинематической вязкостью V и коэффициентом поверхностного натяжения а(Т) ограничена двумя твердыми плоскостями х = 0; х = Ь (рис. 1). Здесь у = 0 — плоскость симметрии (штрихпунктирная линия); f = f (Ь,х) — свободная поверхность.

Движение жидкости в пленке и теплообмен будем описывать системой уравнений Навье-Стокса, которая в переменных ф, ш, в имеет вид

дш д ( дф\ д ( дф\ .

— + — ш^- — ш^- = Аш; (1.1)

дЬ дх \ ду) ду \ дх

Аф + ш = 0; (1.2)

£+£ (вв1) —1 (вдф)=1Ав. (,з)

дЬ дх \ ду) ду \ дх) Рг

В качестве масштабов длины, скорости и давления выбраны величины к0 (половина толщины пленки), v/h0, соответственно, так что число Рейнольдса Ие = 1; Рг — число Прандтля; в = (Т — Т0)/АТ (Т0 — характерное значение температуры, АТ — перепад температур).

Функция тока введена соотношениями

дф дф и = —, V = — —. ду дх

Будем считать, что начальные условия заданы и соответствуют тому, что пленка плоская и находится в состоянии покоя.

При х = 0; х = Ь граничные условия соответствуют условиям прилипания:

ф = »• = 0

Для определения вихря ш на твердой стенке используется условие Тома [2]. На плоскости симметрии

ф = 0, ш = 0.

Рис. 1. Плоская пленка

Граничные условия для температуры всюду, где это необходимо, можно задать единообразным способом:

л пдО

аО + ¡3— = 7 (*).

Тогда, например, на плоскости симметрии а = 7 = 0; в = 1.

На свободной поверхности f (х,Ь) имеют место следующие условия. Определим векторы нормали и касательной к свободной поверхности f (х,Ь) в каждой ее точке:

— fx 1 I I 1 Ух

п

лЛ+Л' лЛ+ТХТ 1\Л+П' \А+7Х

и предположим, что коэффициент поверхностного натяжения ст(Т) зависит от температуры линейно:

1 $ст

ст(Т) = сто(1 - стт(Т - То)), сто = ст(То), стт =--™ |т=То, сто,стт > 0.

сто а±

Для определения свободной поверхности используем кинематическое условие

Л + >Л+£ ^ = 0, (1.4)

а граничные условия для ф и ш запишем в явном виде:

дф ди

(1.5)

ч Я дв ) дв

где у3 — решение уравнения

дуч „д 2у

ш = 21 | + + МП дв, (1.6)

+ у- ав = 2 + О (1-7)

где

■дШ + Са-1

ди

дУз

дЬ

д ' 1

дв Я

стт АТв сто

р. д / уп \ г Уп

- 2 дв ( Я )+ ^ + Л Я;

4-V-'

Нп = Л/лЛ + ТХ. (1.8)

Здесь у3, уп — касательная и нормальная составляющие скорости для точек, лежащих на свободной поверхности; Я — радиус кривизны поверхности Л(х,Ь); Са = роуои/сто — капиллярное число; Мп = сттДТ/(роуои) — число Марангони.

Замечание. Для реальных жидкостей величина сттДТв/сто ^ 1 и практически не влияет на положение свободной поверхности пленки. В дальнейшем изложении работы она опускается.

В выражении для О объединены вместе все главные и второстепенные члены правой части уравнения (1.7) относительно у8. Два слагаемых в первом объединении показывают, какие силы вызывают деформацию свободной поверхности пленки. Первые из них дш/ди представляют собой термокапиллярные силы, второе слагаемое Са-1д(1/Я)/дв — капиллярные силы. Что касается второго объединения (второстепенные члены), то его

можно было бы вообще опустить, например, в случае, если мы ищем стационарное решение задачи.

Для получения формулы (1.6) было использовано соотношение, выражающее непрерывность тангенциальных напряжений на свободной поверхности. Уравнение (1.7) было получено следующим образом. Если векторное уравнение Навье—Стокса, записанное в естественных переменных (V, Р), скалярно умножить на вектор в(Ь,х), то получим скалярное уравнение относительно vs. Производную дР/дв в правой части полученного уравнения можно исключить, если продифференцировать также по в соотношение, выражающее непрерывность нормальных напряжений на свободной поверхности. Уравнение (1.7) играет очень большую роль. Кроме того, что с его помощью можно получить граничные условия для искомых функций ф и ш в явном виде, оно, как будет показано ниже, — важный инструмент для исследования рассматриваемого процесса. Вообще говоря, именно уравнение (1.7) осуществляет обмен информацией между течением внутри области и ее свободной границей. Вывод уравнения (1.7) для стационарного случая приведен в [3].

2. Метод решения

При такой развязке граничных условий на свободной поверхности для решения задачи (1.1)—(1.8) можно использовать любые методы, применяемые в задачах тепло- и массопереноса в замкнутых областях в переменных ф, ш, в. Отметим,что в рамках рассматриваемой модели уравнения (1.1)—(1.3) одного типа и для их решения может быть использована одна и та же вычислительная процедура.

В настоящей работе применялся метод решения в фиксированной области. Область, занятая жидкостью, отображается на прямоугольник со сторонами 0 < £ < Ь, 0 < п < 1 с помощью преобразования

х = £ у = п^(х,ь).

Тогда все границы области, включая свободную поверхность, будут лежать на координатных линиях, а каждое из уравнений (1.1)—(1.3) можно записать в следующем виде:

дФ 1

дЬ Bf

А (В — В -- АФ дф д£\ 1 1 д£ 1 2 дп дп

+

д / дФ л дФ лждф +— I В12—т + В22— + АФ^

дп V д£

дп

д£

Здесь

для Ф = ш: для Ф = в: для Ф = ф:

В11 = f (£), В12 = п, В22

дФ

+ --+ Л2.

дп

1 + В122

f(£) '

(2.1) (2.2)

В =1, А =1, Я1 = ^п, Я2 = 0;

В = Рг, А = Рг, Я1 = у п, Я2 = 0;

В = !, А = 0, Я1 = 0, Я2 = Лш; Л

Л — итерационный параметр, используемый при решении уравнения Пуассона для ф.

Замечание. В данной работе уравнение Пуассона для функции тока ф решается итерационным методом. То есть вместо уравнения Дф = — ш методом установления решается уравнение дф/дЬ = А(Дф + ш). Физический смысл нестационарного решения здесь не играет роли, но когда решение такого уравнения диффузии приближается к стационарному, оно стремится к интересующему нас решению Пуассона (1.2) (см. [2, 4, 5]). Этот прием позволяет все три уравнения для определения искомых функций 9,ш и ф численно решать с помощью одной и той вычислительной процедуры. Основная трудность при решении уравнения Пуассона с переменными коэффициентами итерационным методом состоит в выборе оптимального параметра А. Задача усложняется тем, что при замене переменных и переходе от криволинейной области к прямоугольнику коэффициенты этого уравнения сами меняются на каждом шаге по времени. Если взять за основу соображения, изложенные в [5] относительно выбора оптимального параметра для прямоугольной области, то с их помощью можно получить соотношения, позволяющие выбрать оптимальный параметр для случая, рассматриваемого в настоящей работе. Оптимальный параметр А на равномерной сетке определяется из соотношения

0.5тА

де Дп

\/а*и а

— > 22

а

11

а.

22

Атта1Ъ Ат

Атаха22) Ап

/т

шт(Вп, В22), а11 = 4вт2

/т

шах(Вп, В22), а22 = 4 сое2

где В11,В22 — коэффициенты матрицы перехода В^. Введем обозначения:

дФ дФ дф и (Ф) = В11 + В12 дф — АФ дф, д£ дп дп

дФ дФ дф V (Ф) = В,2 _ + В22 - + АФ

Тогда уравнение (2.1) примет вид

пде

п Дп

дФ 1 дФ

Ж = В/ ^(Ф) + ^(Ф)] + * Щ + Л2-

(2.3)

а граничные условия на свободной поверхности (п = 1 ) можно записать следующим образом:

/ + Ц = 0- (2.4)

.)л/1 + В2.

(2.5)

(2.6)

дф = (уя + В12^„^1 + В2: дп

12

В

22

= ^ , 1 дУп 1 , мп ш ^ + ллтб2 де1 ^лд+в^ де-

где — решение уравнения

1

*

2

1

*

2

^ . дVs д ( 1 дVs . , .

«" +Ж = 2 ут+Ж её' + в- (2-7)

д^ -1 д / 1 \ д /vn\ „ V,

2

V к+1/2(Ф) = В12 Ф?к + В22 Ф^+1/2 + АФк+1/2

ик+1(Ф) = В11Фк+1 + В12Ф^+1/2 — АФк+1 (2.10)

В = — ^,2 + В22 Тп) + ОТ' — 2 ^ ) + П„ш — В.2 ^

Vn = Л/\/1 + В22. (2.8)

Заметим, что уравнение (2.3) представлено в дивергентной форме, а кинематическое условие (2.4) выражает закон сохранения массы.

Решение уравнения (2.3) для функций в (температура), ш (вихрь) и ф (функция тока) на каждом временном шаге будем искать по схеме стабилизирующей поправки [4, 6], взятой в форме

Ф&+1/2 _ 1

---= В] (Ф) + ^+1/2(Ф)] + Д1Фк+1/2 + Я2,

Ф^+1 _ Ф^+1/2 1

-т-= В] [Ц?+1(Ф) — ик(Ф)] , (2.9)

д£

{дф дп

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

Для реализации схемы (2.9), (2.10) в прямоугольнике, соответствующем преобразованной области, строится расчетная сетка стандартным образом:

£п = (п — 1)А£, А£ = Ь/ЫВ, п = 1,..,ЫЫ, NN = NB + 1,

пт = (т — 1)Ап, Ап =1/МВ, т =1,...,ММ, ММ = МВ + 1.

Дифференциальные выражения типа (а11Ф^, (а22Фп)п, (а1Ф)^, (2Ф)П, (а12Ф^)п, (а12Фпаппроксимируем со вторым порядком точности конечно-разностными аналогами Л11, Л22, Л1, Л2, Л12, Л21, которые имеют традиционное представление.

После замены в (2.9) производных соответствующими конечными разностями и подстановки вместо Уп(Ф) и Ц-(Ф) их разностных аналогов из (2.10) на каждом полушаге по времени для всех внутренних точек (п = 2,..., ЫВ; т = 2,..., МВ) получим систему линейных разностных уравнений относительно функции Ф(£п, пт). Система имеет трех-диагональную структуру с преобладанием диагональных элементов матрицы и может быть эффективно решена методом прогонки с учетом специфики граничных условий.

Уравнение (2.7) относится к уравнениям типа Бюргерса—Хопфа с правой частью, и решать его следует очень аккуратно, так как от этого зависит получение граничных условий для нахождения функций ф и ш. Для решения уравнения (2.7), следуя [7, 8],

была использована консервативная схема второго порядка по пространственной переменной, аппроксимирующая дивергентную форму записи уравнения (2.7):

дУ8 1 д / N2 о д ( 1 1 , п

аГ + 2де ^ =2де (тггЖ аё1 +

Левую гиперболическую часть этого уравнения аппроксимируем разностным оператором

МП+1 — МП , 1ЫП+1ЫПЙ — МП-^)П-1

т +2 2Де -

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

Первое слагаемое в правой части уравнения (2.7) аппроксимируем на верхнем временном шаге (к +1) разностным оператором Л11. Остальные дифференциальные операторы, входящие в Д аппроксимируются стандартным образом. Для производной дш/дп на свободной поверхности используется односторонняя аппроксимация второго порядка, причем значение дш/дп берется с предыдущего временного шага. В результате разностная схема, аппроксимирующая уравнение (2.7), имеет второй порядок по пространственной переменной и первый порядок по времени. Получаемая при этом система разностных уравнений имеет трехдиагональную структуру и решается методом прогонки. Граничные условия для решения (2.7) можно определить из граничных условий для ф, заданных на боковых стенках рассматриваемой области, а также из физической постановки задачи.

Общий алгоритм решения задачи осуществляется следующим образом. Пусть для момента времени ¿к = кт известно положение свободной поверхности /, а также распределение температуры 0, вихря ш и функции тока ф для всей области. По формуле (2.4) найдем новое положение свободной поверхности /, а по формулам (2.2) определим матрицу коэффициентов В11, В12, В22 уравнения (2.3). Наконец, решив это уравнение с новыми коэффициентами последовательно для температуры 0, вихря ш, функции тока ф, получим распределение этих функций для момента времени ¿к+1 = (к + 1)т на итерации в = 0. Итерации на каждом шаге по времени проводятся до тех пор, пока не будет выполнено условие

шах 1(Ф(ега))-+1 | — |(Ф(еп))'| <е шпах |(Ф(еп ))*+ч <е

для всех функций ш, 0, ф, где в — номер итерации; е — заданная точность.

На этом заканчивается расчет одного шага по времени. Все остальные шаги рассчитываются аналогично. Если мы ищем стационарное решение задачи, то оно считается полученным при выполнении условия

шах 1(/(Ы)к+К| — |(/(Ы)к| <е ш |(/(еп))к+к | <е-

где к — номер временного шага, К — некоторое заданное число шагов, е — заданная точность.

Из описания следует, что, во-первых, вычислительный процесс, организованный таким образом, не нарушает законов сохранения, а во-вторых, на каждом временном шаге и на каждой итерации сохраняется эллиптический вид аппроксимирующего оператора. Эти свойства и обеспечивают сходимость итерационного процесса в целом [5, 9].

3. Численные расчеты и обсуждение результатов

Рассмотрим деформацию и разрыв свободновисящей пленки под действием тепловой нагрузки. Наиболее интересным, на наш взгляд, является случай, когда капиллярное число Са и число Марангони Мп достаточно малы, а тепловая нагрузка недостаточно велика, чтобы привести к разрыву пленки. На рис. 2 представлены результаты расчета, проведенного для параметров Са-1 = 1500, Мп= 0.5, Рг= 1; отношение длины пленки к половине ее толщины равно 100. Безразмерная половина толщины пленки к0 равна единице. Температура в(х, у, в начальный момент времени равна нулю. На свободную поверхность пленки, находящейся в состоянии покоя, наложим тепловую нагрузку

в(х,г) = Аяп(пх/Ь), 0 < х < Ь, А = 1, (3.1)

где Ь — длина пленки. Для х = 0 и х = Ь в = 0;на плоскости симметрии дв/ди = 0.

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

С помощью уравнения (1.7) это обстоятельство можно прокомментировать следующим образом. Когда пленка находится в состоянии покоя, то правая часть уравнения (1.7) равна нулю. Но как только накладывается тепловая нагрузка, то правая часть уравнения (1.7) становится отлична от нуля за счет слагаемого дш/ди, где ш определяется формулой (1.6). Производная дш/ди имеет разные знаки по обе стороны от середины пленки, так как градиент температуры в этих областях направлен в разные стороны.

Рис. 2. Положение свободной поверхности пленки для разных моментов времени t: ti < t2 < t3 < t4 < tend; t2 — ti = P/2, t4 — t3 = P/2, P — период осцилляций; f (x,tend) — результат стационарного решения задачи

В силу этих причин касательная скорость также имеет разные знаки и направлена от середины пленки в разные стороны. Тоньше всего пленка в середине (рис. 2, /(х,Ь1)). Термокапиллярные силы как будто стремятся растянуть пленку в противоположных направлениях. Однако с дальнейшей деформацией пленки и уменьшением ее толщины растет кривизна свободной поверхности и слагаемое с производной от кривизны свободной поверхности Са-1д(1/Я)/дв в правой части уравнения (1.7) начинает играть существенную роль. Имея также разные знаки с дш/дп, оно растет вместе с уменьшением толщины пленки. Заметим, что в это же время термокапиллярные силы убывают, так как пленка прогревается. Как только слагаемое с производной от кривизны превысит значение дш/дп, утончение пленки прекращается, касательная скорость меняет знаки на противоположные и свободная поверхность стремится занять первоначальное положение (рис. 2, /(х,Ь2)). Далее процесс повторяется до тех пор, пока между капиллярными и термокапиллярными силами не будет достигнут некоторый компромисс — он приведет в дальнейшем к получению стационарного решения. На рис. 2 положение свободной поверхности, соответствующее стационарному решению, обозначено сплошной линией. Характерно, что на свободной поверхности имеются две узловые точки, где пленка сохраняет постоянную толщину, равную единице. Колебание свободной поверхности похоже на колебание тонкой мембраны, закрепленной на расстоянии Ь/4 от концов пленки.

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

Случай а:

0(ж,г) = А81п(пж/Ь), 0 < X < Ь, А =15.

1 «-1 1 о о

1.61.20.80.40-

б

1.61.2-

в

0.80.40-

Рис. 3. Положение свободной поверхности в случаях а —

в

Случай б:

Рис. 4. Пример расчета случая в

(х,г) = А, Ь/2 - йо < х < Ь/2 + ^0,

в остальных точках на свободной поверхности дв/ди = 0; й0 — половина толщины пленки.

Случай в:

в(х, г) = А, Ь/2 - й0 < х < Ь/2 + й0,

в остальных точках свободной поверхности в = 0.

Положение свободной поверхности во всех трех случаях а — в для одного и того же момента времени ¿* показано на рис. 3. Как видно из этого рисунка, пленка имеет наименьшую толщину в случае а. Однако если продолжить расчет, то в случае б толщина пленки будет и дальше уменьшаться до тех пор, пока это утончение не приведет к разрыву в ее середине. В случаях а и в свободная поверхность пленки будет совершать колебательные движения относительно некоторого положения равновесия. На рис. 3 сплошной линией обозначено стационарное решение для этих случаев.

И хотя небольшое увеличение коэффициента А приведет к разрыву пленки в случае а, в случае в для параметров, описанных выше, это увеличение приведет лишь к очень незначительному утончению пленки. Тем не менее для других параметров Са и Мп даже при меньшей тепловой нагрузке в этом случае может наступить разрыв пленки. На рис. 4 приведен пример расчета случая в, выполненного для параметров Са-1 = 50; Мп = 1.5; Рг = 1; А = 5. Разрыв пленки происходит так быстро, что возмущения, вызванные тепловой нагрузкой, не успевают добежать до концов пленки.

Заключение

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

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

[1] Пухначев В.В., ДувинкинА С.Б. Модель деформации и разрыва пленки под действием термокапиллярных сил // Изв. РАН. МЖГ. 2006. № 5. C. 89-107.

[2] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980.

[3] ОВЧАРОВА А.С. Метод расчета стационарных течений вязкой жидкости со свободной границей в переменных вихрь — функция тока // ПМТФ. 1998. Т. 39, № 2. С. 59-68.

[4] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[5] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

[6] Douglas J., jr., Rachford H.H., jr. On the numerical solution of heat conduction problems in two and three space variables // Transactions of the Amer. Math. Soc. 1956. Vol. 82, N 2. P. 421-439.

[7] Остапенко В.В. Метод теоретической оценки дисбалансов разностных схем на ударной волне // Докл. АН СССР. 1987. Т. 295, № 2. С. 292-297.

[8] Остапенко В.В. О построении разностных схем повышенной точности для сквозного расчета нестационарных ударных волн // Журн. вычисл. математики и мат. геофизики. 2000. Т. 40, № 12. С. 1857-1874.

[9] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1977.

Поступила в редакцию 15 февраля 2007 г., в переработанном виде —16 апреля 2007 г.