Как стать автором
Обновить

Задача о форме капли жидкости на потолке

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



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

image

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

Найдем уравнение для этих сил, приравняем их и из этого уравнения определим искомую поверхность. Начнем с сил поверхностного натяжения. Мысленно вырежем из капли цилиндр на расстоянии r и малой толщиной delta r. Жидкость, находящаяся в этом цилиндре удерживается силой поверхностного натяжения. Сила поверхностного натяжения T является касательной к поверхности, однако, нам необходимо знать силу по оси y, которая равна T умноженной на синус угла наклона поверхности к потолку т.е. оси r. Как известно, производная y по r равна тангенсу этого угла

image

Чтобы найти синус угла, выразим тангенс следующим образом

image

Откуда

image

Таким образом, интересующая нас сила равна

image

Как известно, сила поверхностного натяжения равна коэффициенту поверхностного натяжения альфа умноженного на длину края поверхности

image

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

image

Вынесем общую часть за скобки и разложим функцию по малому параметру

image

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

image

Преобразуем знаменатель

image

Сделаем два корня в знаменателе

image

За счет малости параметра, избавимся от корня во втором множителе в знаменателе

image

Также за счет малости параметра перетащим второй множитель из знаменателя в числитель

image

Раскроем скобки, приведем подобные и избавимся от слагаемых второго порядка малости

image

Приведем к общему знаменателю

image

Еще раскроем скобки, убрав слагаемые второго порядка малости

image

Окончательно получим

image

Полученная сила поверхностного натяжения должна уравновешивать вес жидкости в цилиндре, с силой тяжести равной

image

где ро — плотность жидкости, а g — ускорение свободного падения.
Здесь необходимо поставить знак минус, поскольку функция y отрицательна.

Введем параметр

image

Приравняв силу поверхностного натяжения к силе тяжести и сократив общие множители, получим

image

Выразим явно вторую производную

image

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

#include "stdio.h"
#include "math.h"

double old_t;
double old_y, old_dy;
double t;
double y, dy;
double delta;

const double A = 1000.0 * 9.8 / 0.073;
/*const double h = -0.0065;*/
const double h = -0.0045;


double funcY( double t, double yn, double dyn );
void run( void );

int main( void )
{
  int    i;
  int    j;
  FILE * OutputF;
  int    RetCode;

  RetCode = 0;

  OutputF = fopen( "result.txt", "wt" );

  if( OutputF != NULL ) {

    for( i = 0; i < 1; i++ ) {
      old_t = 0.0;
      t = 0.0;
      old_y = h + 0.001 * i;
      old_dy = 0.0;
      delta = 0.0001;



      while( old_y <= 0.0 ) {
        for( j = 0; j < 101; j++ ) {
          fprintf( OutputF, "%f %f %f\n", old_t * 1000.0 * cos( 0.0628 * j ),
                                          old_t * 1000.0 * sin( 0.0628 * j ),
                                          old_y * 1000.0 );
        }

        fprintf( OutputF, "\n");
        run();
      }
    }

    fclose( OutputF );
  }
  else {
    RetCode = 1;
  }

  return( RetCode );
}

void run( void )
{
  double tmp_tn,tmp_yn,tmp_dyn;
  double tk0y,tm0y,tk1y,tm1y,tk2y,tm2y,tk3y,tm3y;

  tk0y = delta * old_dy;

  tmp_tn = old_t;
  tmp_yn = old_y;
  tmp_dyn = old_dy;

  tm0y = delta * funcY(tmp_tn,tmp_yn,tmp_dyn);

  tmp_tn = old_t + (delta/2.0);
  tmp_yn = old_y + (tk0y/2.0);
  tmp_dyn = old_dy + (tm0y/2.0);

  tk1y = delta * (old_dy + (tm0y/2.0));
  tm1y = delta * funcY(tmp_tn,tmp_yn,tmp_dyn);

  tmp_yn = old_y + (tk1y/2.0);
  tmp_dyn = old_dy + (tm1y/2.0);

  tk2y = delta * (old_dy + (tm1y/2.0));
  tm2y = delta * funcY(tmp_tn,tmp_yn,tmp_dyn);

  tmp_tn = old_t + delta;
  tmp_yn = old_y + tk2y;
  tmp_dyn = old_dy + tm2y;

  tk3y = delta * (old_dy + tm2y);
  tm3y = delta * funcY(tmp_tn,tmp_yn,tmp_dyn);

  y = old_y + (tk0y + 2.0 * tk1y + 2.0 * tk2y + tk3y)/6.0;
  dy = old_dy + (tm0y+ 2.0 * tm1y + 2.0 * tm2y + tm3y)/6.0;

  t = t + delta;
  old_t = t;
  old_y = y;
  old_dy = dy;
}

double funcY( double t, double yn, double dyn )
{
  double tmp;

  tmp = 1 + dyn * dyn;

  if( t == 0.0 ) {
    return( -A * yn * tmp * sqrt( tmp ) );
  }
  else {
    return( -A * yn * tmp * sqrt( tmp ) - dyn * tmp / t );
  }
}

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

image

Эта программа кроме функции y выдает еще линии уровня по которым можно построить 3D модель капли. Эту задачку я оставляю читателям.

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

image

Интересно, что оказывается существует некоторый максимум объема в районе высоты в 5 мм, который поверхностное натяжение может удерживать. На этом мое исследование этой задачи заканчивается, пишите свои соображения в комментариях.
Теги: форма капли жидкостидифференциальные уравнения
Хабы: Программирование Физика Химия
Всего голосов 10: ↑10 и ↓0 +10
Комментарии 22
Комментарии Комментарии 22

Похожие публикации

Лучшие публикации за сутки