Как научиться CFD с нуля: путь, который учит проверять ответы

Красивая картинка и три вопроса

Студент показывает вам красивое поле течения вокруг автомобиля. Плавные цвета, правдоподобные линии тока, коэффициент сопротивления с тремя знаками после запятой. Вы задаёте три вопроса: какой был y+ на поверхности, сошёлся ли баланс массы и насколько изменилось сопротивление при измельчении сетки? Ни на один нет ответа.

Этот студент освоил программу для CFD. Сам CFD он ещё не освоил.

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

Что значит научиться CFD?

Это три вещи, а программа идёт четвёртой и значит меньше всего.

• Физика. Что говорят уравнения сохранения, что делает каждый член и какие члены доминируют в вашей задаче. Числа Рейнольдса, Пекле и Маха говорят о течении до построения сетки больше, чем большинство цветных картинок после расчёта.
• Численные методы. Как непрерывное уравнение превращается в алгебру и какую ошибку это вносит. Любой ответ метода конечных объёмов содержит погрешность дискретизации, и умение состоит в том, чтобы знать её величину и характер.
• Суждение. Решить, что моделировать, проверить, что числа сошлись, балансы замкнуты и результат не чувствителен к сетке, и сравнить с чем-то независимым. Это верификация и валидация, и именно за это платят работодатели и заказчики.
• Программа служит интерфейсом ко всему перечисленному. Её можно освоить за несколько недель, когда остальное уже на месте, и она всё равно меняется с каждым выпуском.

Какая математика нужна и зачем

Меньше, чем многие боятся, и больше, чем признаёт большинство учебных роликов.

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

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

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

Напишите маленький решатель и сломайте его

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

Хорошая бесплатная отправная точка CFD Python, «12 шагов к уравнениям Навье-Стокса» Лорены Барбы и Гилберта Форсайта: блокноты Jupyter, которые ведут от одномерной линейной конвекции через диффузию, уравнение Бюргерса, уравнения Лапласа и Пуассона к двумерному течению в каверне и в канале. Материалы распространяются по лицензии CC-BY 4.0, код по BSD-3, так что их можно свободно использовать и адаптировать.

Затем проведите тот единственный эксперимент, который каждый инженер CFD должен проделать своими руками.

Пример расчёта: сеточное число Пекле

Возьмём стационарную одномерную конвекцию и диффузию скаляра φ на отрезке 0 ≤ x ≤ L при φ(0) = 0 и φ(L) = 1. Уравнение и его точное решение:

При Pe = 50 решение почти равно нулю по всей области и резко растёт в тонком слое у правого конца. При x = 0,9 L точное значение около exp(−5) ≈ 0,007.

Разобьём отрезок на N одинаковых ячеек. Всё определяет сеточное число Пекле Pe_Δ = ρuΔx / Γ = Pe / N. На равномерной сетке каждая схема сводится к рекуррентному соотношению, решения которого растут от узла к узлу с постоянным множителем, и этот множитель показывает, что делает схема:

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

На 10 ячейках Pe_Δ = 5:

На 40 ячейках Pe_Δ = 1,25: центральная схема даёт отношение около 4,3 при точном 3,5, гладко и близко; схема против потока даёт 2,25 и всё ещё заметно размазывает решение.

Отсюда следуют два урока, и это самые важные уроки статьи.

Центральная схема осциллирует при Pe_Δ > 2. Отношение становится отрицательным, и в решении появляются колебания, к физике отношения не имеющие. Это классический результат из учебника Патанкара, и именно поэтому коммерческие решатели в большинстве течений не используют чистую центральную схему для конвекции.

Схема против потока никогда не осциллирует и ошибается так, что выглядит правильной. Ошибка усечения из разложения Тейлора выше, (Δx/2)φ″, действует в точности как дополнительная диффузия. Схема ведёт себя так, будто коэффициент диффузии равен:

На 10 ячейках это 3,5 Γ. Ваше гладкое решение по схеме против потока в главном порядке является решением задачи с числом Пекле около 50/3,5 ≈ 14, а вовсе не 50. На экране об этом ничего не говорит. Изолинии чистые, невязки красиво падают, а ответ относится к другому течению.

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

Любой современный решатель предлагает ограниченные схемы повышенного порядка (второй порядок против потока с ограничителями, QUICK, семейства TVD), которые находятся между этими крайностями. Понять, что делают эти настройки, вы сможете, только увидев, как ошибаются обе крайности.

Проверьте порядок точности

Теперь измельчайте сетку. Вычислите ошибку относительно точного решения на двух сетках с коэффициентом измельчения r. Если e₁ ошибка на грубой сетке, а e₂ на мелкой, наблюдаемый порядок равен:

Пример расчёта с условными числами. Схема второго порядка даёт ошибки 0,020 и 0,0051 на сетках с шагами Δx и Δx/2, значит p = ln(3,92)/ln(2) ≈ 1,97. Верно. Если схема, которая должна иметь второй порядок, даёт 1,0, где-то есть что-то первого порядка: граничное условие, градиент, интерполяция на грани. Найти это лучшая численная школа, какая бывает. Тот же расчёт лежит в основе индекса сеточной сходимости, пошагово разобранного в нашей статье о сеточной независимости.

Поднимайтесь по лестнице задач с опубликованными ответами

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

Ступень 3 заслуживает замечания. Длина ламинарного начального участка примерно 0,05 Re D, так что при Re = 1000 профилю нужно около 50 диаметров, чтобы установиться. Короткая труба даёт коэффициент трения заметно выше 64/Re, и это развивающееся течение делает то, что должно. Научиться отличать одно от другого и есть цель.

На ступени 6 большинство самоучек обнаруживает, что модель турбулентности это модель. Сравните профиль средней скорости в пристеночных единицах с DNS один раз при первой ячейке на y+ ≈ 1, а другой раз с пристеночными функциями при y+ ≈ 50. Статья о k-omega SST и k-epsilon в этом блоге объясняет, почему они ведут себя по-разному.

Ступень 5 к тому же хорошая проверка для любого, кто или что берётся ставить для вас расчёты CFD; статья об ИИ-агентах для CFD использует её именно так.

Потом осваивайте программу

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

Хорошо освойте один пакет, прежде чем пробовать несколько. Для учёбы достаточно двух бесплатных путей. OpenFOAM открыт, работает в Linux или в Windows через WSL. Ansys Student бесплатен для студентов и преподавателей, работает в Windows и допускает до 1 миллиона ячеек в задачах течения на четырёх ядрах, только для учебных целей. Миллиона ячеек хватает для каждой ступени описанной лестницы, и это не случайно.

Повторите ступени 2, 3 и 5 в пакете. Если ваша самописная каверна и коммерческая совпадают друг с другом и с Ghia, вы проверили обе и освоили настройки пакета через сравнение, а не через доверие.

Книги, которые стоит иметь

• Versteeg и Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Самый понятный первый учебник по методу конечных объёмов, которым пользуются коммерческие коды.
• Патанкар, «Численные методы решения задач теплообмена и динамики жидкости». Короткая, старая и по-прежнему лучшая книга о том, почему схемы ведут себя так, как ведут.
• Ferziger, Perić и Street, Computational Methods for Fluid Dynamics. Справочник по численным методам, когда основы уже пройдены.
• Moukalled, Mangani и Darwish, The Finite Volume Method in Computational Fluid Dynamics. Мост от теории к реализации в OpenFOAM.
• Pope, Turbulent Flows. Для момента, когда вопрос «какая модель турбулентности?» становится настоящим.

План на 16 недель

Мой вариант для человека с инженерным образованием и примерно восемью-десятью часами в неделю. Каждый блок заканчивается результатом, который можно показать.

Последний результат как раз тот, благодаря которому вас возьмут на работу. Отчёт со словами «перепад давления 1245 Па при численной неопределённости 1,6%, расхождение с корреляцией в пределах 4%» стоит больше любого количества цветных картинок.

Типичные ошибки и как они выглядят на экране

• Гладкое решение, которое ни с чем не сравнивали. Его дала схема против потока. Оно может отличаться от истины в разы, как показывает пример.
• Колебания, «вылеченные» переходом на первый порядок. Колебания говорили, что сетка слишком груба для схемы. Первый порядок прячет сообщение и сохраняет ошибку.
• Учебный пример, повторённый щелчок за щелчком. Он работает, а объяснить, почему каждая настройка именно такая, вы не можете. Меняйте по одной настройке и предсказывайте эффект до запуска.
• Коэффициент трения, «неверный» на 30% в короткой трубе. Это начальный участок, и решатель прав.
• Турбулентный расчёт без проверки y+. Пристеночная обработка ни разу не была в своей области применимости, и модель делала не то, что вы думаете.
• Невязки как единственное свидетельство сходимости. Следите за величиной, которую будете указывать, и проверяйте баланс массы; статья о сходимости в этом блоге объясняет почему.

Что бы сделал я

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

Если вы уже работаете в коммерческом пакете, но численной подготовки не получали, пройдите тот же путь в сжатом виде: эксперимент с числом Пекле, каверна против Ghia, исследование на трёх сетках для знакомой задачи. Три выходных. Вы начнёте видеть, под какими из своих старых результатов вы бы ещё подписались.

Если нужна помощь в системе, наш курс «CFD-моделирование течений жидкостей и газов» идёт в том же порядке: постановка, физика, моделирование, верификация и интерпретация, а курс «Теплообмен в инженерных приложениях» охватывает ступень 7 и дальше. В любом случае держитесь правила: ни один результат не покидает ваш стол, пока его не сравнили с чем-то, что вы не считали сами.