Pull to refresh

Метод погруженной границы для чайников

Reading time6 min
Views15K
Отношение между «чистыми» и «прикладными математиками» основаны на доверии и понимании. «Чистые математики» не доверяют «прикладным математикам», а «прикладные математики» не понимают чистых математиков.


Некоторое время назад я столкнулся с тем, что не смог найти доступных материалов, в которых на великом и могучем описывался бы метод погруженной границы. Если в кратце, то это метод вычислительной гидродинамики, который позволяет расчитывать обтекание достаточно сложных по форме и динамике объектов. Так вот русскоязычных публикаций на эту тему было крайне недостаточно. «Не беда, будем читать работы зарубежных коллег» — подумалось мне. Но и тут ждал небольшой подвох — все имеющиеся материалы и публикации по этому методу были очень теоретичны, а мне (не уверен, возможно это не только моя особенность) обычно сложно сделать переход от теоретических выкладок к более-менее рабочему воплощению в коде. Поэтому для таких же несчастных, что и я (и с некоторой долей надежды на советы от умудренных опытом личностей) я решил сделать краткое описание этого метода и предложить самый простой способ его реализации.

Я буду рассматривать только оригинальный метод Пескина, т.к. на мой взгляд он проще, чем метод фиктивных ячеек (ghost-cells), метод усеченных ячеек (cut-cells) и все прочие их модификации.

Итак, википедия говорит нам, что «метод погруженной границы» (immersed boundary method) — это подход к моделированию взаимодействия между жидкостью и обтекаемым объектом (т.к. часто применяется для моделирования тонких волокон в биологических системах, то иногда используют термин «fiber»). И это чертовски верно =)
Основная особенность данного подхода заключается в том, что вводятся две раздельные сетки для расчета течения жидкости (в эйлеровых координатах) и для расчета параметров погруженной границы (тех самых волокон-фибр в лагранжевых координатах). Это позволяет использовать для расчета жидкости простые сетки и быстрые вычислительные методы, перенося всю сложность моделирования на этап взаимодействия с погруженной границей.

С чисто теоретической точки зрения, такой фокус проделывается с помощью добавления в уравнение Навье-Стокса силы воздействия погруженной границы на жидкость (точнее, плотности силы):



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






с краевыми условиями на погруженной границе, например, такого вида (условие прилипания aka «no-slip condition»):



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

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

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

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

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

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





Теперь вместо тысячи слов можно описать получившийся алгоритм в виде кода.

Расчет напряжения на узлах может выглядеть примерно так:

    for(int n = 0; n < boundary->nodes_count; ++n) {
        Node *node = &boundary->nodes[n];

        // get_area возвращает площадь поверхности, которая относится к данной точке
        node->x_force = -boundary->stiffness * (node->x - node->x_ref) * boundary->get_area();
        node->y_force = -boundary->stiffness * (node->y - node->y_ref) * boundary->get_area();
        node->z_force = -boundary->stiffness * (node->z - node->z_ref) * boundary->get_area();
    }

Распределение силы на точки жидкости:

    for(int n = 0; n < boundary->nodes_count; ++n) {

        // для узла границы находим ближайшую расчетную точку жидкости
        int x_int = index(boundary->nodes[n].x, COORD_X);
        int y_int = index(boundary->nodes[n].y, COORD_Y);
        int z_int = index(boundary->nodes[n].z, COORD_Z);

        // т.к. в дискретном варианте дельта-функции ненулевыми являются те узлы,
        // которые находятся в двух шагах по сетке, то необходимо распределить силу
        // на несколько (по две в каждую сторону) точек жидкости
        for(int i = x_int-1; i <= x_int + 2; ++i) {
            for(int j = y_int-1; j <= y_int + 2; ++j) {
                for(int k = z_int-1; k <= z_int + 2; ++k) {

                    // находим расстояние от узла границы до текущей точки жидкости
                    long double dist_x = fabs(boundary->nodes[n].x - coord(i, COORD_X));
                    long double dist_y = fabs(boundary->nodes[n].y - coord(j, COORD_Y));
                    long double dist_z = fabs(boundary->nodes[n].z - coord(k, COORD_Z));

                    long double weight_x = delta(dist_x, COORD_X);
                    long double weight_y = delta(dist_y, COORD_Y);
                    long double weight_z = delta(dist_z, COORD_Z);

                    // прибавляем полученную силу к данной точке жидкости
                    force_X[i][j][k] += (boundary->nodes[n].x_force * weight_x * weight_y * weight_z) * boundary->get_area();
                    force_Y[i][j][k] += (boundary->nodes[n].y_force * weight_x * weight_y * weight_z) * boundary->get_area();
                    force_Z[i][j][k] += (boundary->nodes[n].z_force * weight_x * weight_y * weight_z) * boundary->get_area();
                }
            }
        }

Любым угодным методом производим расчет течения. Затем интерполируем скорость на узлы жидкости:

    for(int n = 0; n < boundary->nodes_count; ++n)
    {
        // по аналогии с силой, находим ближайщую точку
        int x_int = index(boundary->nodes[n].x, COORD_X);
        int y_int = index(boundary->nodes[n].y, COORD_Y);
        int z_int = index(boundary->nodes[n].z, COORD_Z);

        // и собираем скорость с ближайших точек
        for(int i = x_int-1; i <= x_int + 2; ++i)
        {
            for(int j = y_int-1; j <= y_int + 2; ++j)
            {
                for(int k = z_int-1; k <= z_int + 2; ++k)
                {
                    long double dist_x = fabs(boundary->nodes[n].x - coord(i, COORD_X));
                    long double dist_y = fabs(boundary->nodes[n].y - coord(j, COORD_Y));
                    long double dist_z = fabs(boundary->nodes[n].z - coord(k, COORD_Z));

                    long double weight_x = delta(dist_x, COORD_X);
                    long double weight_y = delta(dist_y, COORD_Y);
                    long double weight_z = delta(dist_z, COORD_Z);

                    boundary->nodes[n].x_vel += (velocity_U[i][j][k] * weight_x * weight_y * weight_z) * CB(Hx[i]);
                    boundary->nodes[n].y_vel += (velocity_V[i][j][k] * weight_x * weight_y * weight_z) * CB(Hy[j]);
                    boundary->nodes[n].z_vel += (velocity_W[i][j][k] * weight_x * weight_y * weight_z) * CB(Hz[k]);
                }
            }
        }
    }

Обновляем координаты узлов погруженной границы:

    for(int n = 0; n < boundary->nodes_count; ++n)
    {
        boundary->nodes[n].x += boundary->nodes[n].x_vel;
        boundary->nodes[n].y += boundary->nodes[n].y_vel;
        boundary->nodes[n].z += boundary->nodes[n].z_vel;
    }


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

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


Tags:
Hubs:
Total votes 44: ↑41 and ↓3+38
Comments9

Articles