Расчет равновесной поверхности капли жидкости

Размещено на http://www.allbest.ru

Размещено на http://www.allbest.ru

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ

БЕЛОРУССКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

Факультет прикладной математики и информатики

Кафедра вычислительной математики

Расчет равновесной поверхности капли жидкости

Лаврентик Антон Леонидович

Отчет по с/л

(“Методы решения задач со свободными границами ”)

студента 4 курса 5 группы

Преподаватель Будник Анатолий Михайлович

доцент кафедры выч. мат., канд. физ.-мат. наук

Минск 2012

Постановка задачи

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

Необходимо:

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

Решить полученную задачу итерационно-разностным методом.

Найти решение при следующих значениях физических параметров:

капля жидкость равновесная поверхность

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

Решение

Построение безразмерной математической модели

Для полного описания состояния жидкости достаточно пяти функций:



Первые 3 функции - это компоненты вектора скорости элементарного объема жидкости с координатами в момент времени . Последние две функции - это давление и плотность.

Так как нашей целью является нахождение именно равновесной формы капли жидкости, то

Будем считать жидкость несжимаемой, в этом случае будем иметь:

Получим уравнения, связывающие все эти величины.

На основе уравнения Лапласа и используя соотношения Эйлера, мы можем окончательно получить уравнение равновесия:

где - кривизны главных нормальных сечений равновесной поверхности, - потенциал поля массовых сил.

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

 на линии

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

Введем систему координат, как показано на рисунке:

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

.

Тогда будет равно:

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

Таким образом, имеем дифференциальное уравнение:

где - некоторая константа, которая будет вычислена позже.

К этому уравнению добавятся дополнительные условия:

Условия в граничных точках, следующие из условий Дюпре-Юнга:

:

:

Естественное условие: .

Условие на объем: .

На основе этих данных построим безразмерную модель.

Введем безразмерные величины

Связь между производными будет иметь вид:

Тогда дифференциальное уравнение примет вид:

Введем безразмерное число Бонда , которое характеризует отношение гравитационных сил к капиллярным, получаем:

Выразим из условия на объем:

Получаем уравнение:

Теперь вычислим константу . Для этого сначала проинтегрируем это уравнение по отрезку :

Вычислим интегралы и воспользуемся безразмерным граничным условием . Выразив из полученного уравнения , получаем:

Окончательная безразмерная модель вместе с дополнительными условиями примет вид:

где .

Построение вычислительного алгоритма

В области определения введем сетку

где - число разбиений области. Количество узлов в такой сетке будет равно .

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

Теперь перепишем уравнение, описывающее поверхность капли в виде

где .

Во внутренних точках сетки запишем разностную схему, аппроксимирующую это уравнение со вторым порядком точности:

При получении формулы для использовалось граничное условие .

Осталось аппроксимировать граничные условия. Так как требуется аппроксимация не ниже второго порядка точности, то нужно будет выполнить процедуру повышения порядка точности. Получим уравнения:



Функцию выразим из основного дифференциального уравнения:

Используя граничные условия, получим следующие формулы для и :

Таким образом, разностная схема для решения задачи принимает вид:

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

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

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

Задаем первое приближение к решению, например можно взять прямую .

Запускаем итерационный процесс по . На каждой итерации выполняем следующие действия:

по известному на данный момент приближению вычисляем интеграл

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

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

используя известное на данный момент приближение , а также интеграл и константу , вычисляем коэффициенты и свободные члены СЛАУ.

решаем СЛАУ методом разностной прогонки

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

Реализация алгоритма в КТС Wolfram Mathematica 8.0.1 приведена в приложении А.

Чтобы исследовать зависимость формы поверхности от соотношения действующих сил проведем расчет для 3 различных значений числа Бонда.

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

Результат вычислений

Для заданных в п.3 постановки задачи физических параметров получаем следующее решение:

Рисунок 1 - равновесная форма капли жидкости в безразмерных координатах

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

Рисунок 2 - равновесная поверхность капли

Проиллюстрируем зависимость формы поверхности от величины числа Бонда:

Рисунок 3 - формы капли в зависимости от величины Бонда

;

;

Приложение А

\[Rho]=1; (*плотность жидкости, г/см^3*)

\[Sigma]=72.75; (*коэффициент поверхностного натяжения, дин/см*)

g=981; (*ускорение свободного падения, см/сек^2*)

V=1.07; (*объем капли, см^3*)

\[Alpha]=80; (*угол смачивания, градусы*)

\[Alpha]=\[Alpha]*\[Pi]/180;

Bo=(\[Rho]*g)/\[Sigma]*V^(2/3); (*число Бонда*)

n=1000; (*количество разбиений*)

h=1/n; (*величина шага*)

\[CurlyEpsilon]=h^3; (*допустимая погрешность*)

z=Table[0,{i,1,n+1}];

zprev=Table[0,{i,1,n+1}];

r=Table[(i-1)*h,{i,1,n+1}];

q=0;

(*построение начального приближения*)

For[i=1,i<=n+1,i++,

z[[i]]=1-r[[i]];

];

p={Null,Null,Null};

For[w=0,w<=2,w++,

g=w*981;

Bo=(\[Rho]*g)/\[Sigma]*V^(2/3); (*число Бонда*)

(*запускаем итерационный процесс*)

(*на каждой итерации получаем уточненный профиль капли*)

q=0;

zprev=Table[0,{i,1,n+1}];

While[Norm[zprev-z,\[Infinity]]>\[CurlyEpsilon]&&q<1000,

zprev=z;

k=Table[0,{i,1,n}];

k[[1]]=r[[1]]/Sqrt[1+((z[[2]]-z[[1]])/h)^2];

For[i=2,i<= n,i++,

k[[i]]=r[[i]]/Sqrt[1+((z[[i+1]]-z[[i-1]])/(2*h))^2];

];

I1=2*\[Pi]*h((r[[1]]*z[[1]]+r[[n+1]]*z[[n+1]])/2+\!\(

\*UnderoverscriptBox[\(\[Sum]\), \(i = 2\), \(n\)]\((r[\([\)\(i\)\(]\)]*z[\([\)\(i\)\(]\)])\)\));

const = -2*Sin[\[Alpha]]-Bo/\[Pi]*I1^(1/3);

z[[n+1]]=0;

z[[n]]=z[[n+1]]+h*Tan[\[Alpha]]+h^2/2*(const+Sin[\[Alpha]])/(Cos[\[Alpha]])^3;

a=Table[0,{i,1,n-1}];

For[i=2,i<=n-1,i++,

a[[i]]=(k[[i-1]]+k[[i]])/(2*h^2);

];

c=Table[0,{i,1,n-1}];

c[[1]]=-(1/h);

For[i=2,i<=n-1,i++,

c[[i]]=-((k[[i-1]]+2*k[[i]]+k[[i+1]])/(2*h^2))-Bo*r[[i]]*I1^(-(2/3));

];

b=Table[0,{i,1,n-1}];

b[[1]]=1/h;

For[i=2,i<=n-2,i++,

b[[i]]=(k[[i]]+k[[i+1]])/(2*h^2);

];

b[[n-1]]=0;

F=Table[0,{i,1,n-1}];

F[[1]]=h^2/2*1/2*Bo*z[[1]];

For[i=2,i<=n-2,i++,

F[[i]]=const*r[[i]];

];

F[[n-1]]=const*r[[n-1]]-(k[[n-1]]+k[[n]])/(2*h^2)*z[[n]];

(*для решения системы применяем метод разностной прогонки*)

\[Alpha]1=Table[0,{i,1,n-1}];

\[Beta]1=Table[0,{i,1,n-1}];

\[Alpha]1[[2]]=-b[[1]]/c[[1]];

\[Beta]1[[2]]=F[[1]]/c[[1]];

For[i=1,i<=n-2,i++,

\[Alpha]1[[i+1]]=-b[[i]]/(c[[i]]+\[Alpha]1[[i]]*a[[i]]);

\[Beta]1[[i+1]]=(F[[i]]-\[Beta]1[[i]]*a[[i]])/(c[[i]]+\[Alpha]1[[i]]*a[[i]]);

];

z[[n-1]]=(F[[n-1]]-a[[n-1]]*\[Beta]1[[n-1]])/(c[[n-1]]+a[[n-1]]*\[Alpha]1[[n-1]]);

For[i=n-2,i>=1,i--,z[[i]]=\[Alpha]1[[i+1]]*z[[i+1]]+\[Beta]1[[i+1]]];

q++;

];

Print["Число бонда: ", Bo];

Print["Число итераций: ", q];

p[[w+1]]=ListLinePlot[Table[{(i-1)*h/I1,z[[i]]/I1},{i,1,n+1}],AspectRatio->Automatic,PlotStyle->{Blue,{Red,Thick},Green}[[w+1]]];

];

Show[p[[1]],p[[2]],p[[3]],PlotRange->All]

Размещено на Allbest.ru