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



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



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

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



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



Откуда



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



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



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



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



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



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



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



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



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



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



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



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



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



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



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

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



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



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



решить это уравнение в явном виде не представляется возможным, однако, можно решить задачу Коши для этого уравнения численными методами. Заметим, что для дифференциального уравнения второго порядка необходимо знать два начальных условия. В нашем случае можно задать начальный размер капли (то есть расстояние от потолка до ее нижней части, а также учесть, что при 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 );
  }
}

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



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

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



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