Вычислительная макромеханика как обобщение идей С.К. Годунова

Вычислительный эксперимент стал неотъемлемым инструментом исследования во многих отраслях науки и техники. Ещё сравнительно недавно двухмерное моделирование стационарного обтекания охлаждаемой турбинной лопатки авиационного двигателя казалось волшебством, а сегодня уже никого не удивить цифровыми двойниками сложнейших технических устройств. Однако и у математического моделирования физико-химических процессов, основу которого составляет триада «модель–алгоритм–программа», есть недостатки: попытки создания всё более точных математических моделей приводят, как правило, к заметному усложнению вычислительных алгоритмов и компьютерных программ для численного решения основополагающих уравнений. К факторам, существенно затрудняющим моделирование, следует отнести неустранимую погрешность моделей (неточные начальные и/или граничные условия, приближённые уравнения состояния и другое), сложную геометрию (пористость тел, шероховатость омываемых поверхностей и прочее) и комплексный характер моделируемых процессов (термодинамическая неравновесность, турбулентность, многомасштабность (multiscale simulation), совокупность одновременно протекающих процессов (multiphysics simulation) и другое). По мере усложнения изучаемых систем стала очевидной необходимость снижения издержек, связанных с написанием и отладкой компьютерных программ (встала проблема автоматизации моделирования или формализации вычислений).

Первоначально разработка новых вычислительных технологий велась преимущественно в США и СССР в 50–60-х годах прошлого века и была нацелена на создание оружия массового уничтожения и средств его доставки. Секретностью исследований существенно ограничивалась возможность опубликования разработанных численных методов и проведения международных конференций. Только в 1978 году выдающийся английский учёный Брайн Сполдинг (Brian Spalding) задумал написание открытого кода вычислительной гидродинамики, способного обрабатывать все процессы, связанные с течением жидкости или газа. Следует заметить, что тогда доминировали вычислительные алгоритмы решения упрощённых уравнений Навье–Стокса (например, в приближении пограничного слоя). В то время стандартная статья по вычислительной гидродинамике выглядела следующим образом: 1) постановка задачи, 2) анализ отдельных членов уравнений Навье–Стокса, 3) упрощение уравнений Навье–Стокса, 4) численное решение упрощённых уравнений и сопоставление с экспериментом. Фактически результаты каждой статьи показывали, какими членами этих уравнений можно пренебречь в каждом конкретном случае. Использование полных уравнений Навье–Стокса для моделирования процессов гидродинамики и теплообмена при помощи кода PHOENICS стало революционным шагом (http: //www. cham. со. uk). Однако сразу выяснилось, что основой единого кода должна быть универсальная вычислительная технология, которая давала бы возможность эффективно решать широкий класс нелинейных (начально-)краевых задач механики сплошных сред. Разработка подобной технологии является отдельной и весьма трудоёмкой проблемой [1].

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

Статья опубликована в открытом доступе по лицензии CC BY 4.0

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

С другой стороны, идея универсальных алгоритмов заключается в выборе компонентов независимо от конкретной задачи, равномерно для как можно большего класса задач. Это означает, что универсальные алгоритмы строятся без задействования априорной информации о решаемых задачах, но при условии, что их алгоритмическая трудоёмкость сравнима с трудоёмкостью оптимизированных алгоритмов. На подобных алгоритмах базируются программы, устроенные по принципу «чёрного ящика» (black-box software), где не так важно минимальное время счёта для отдельной задачи [2] . Невозможность прибегнуть к априорной информации упрощает структуру универсальных алгоритмов, но затрудняет их разработку.

В настоящее время известен ряд комплексов программ, существенно упрощающих математическое моделирование физико-химических процессов в приближении сплошной среды: Fluent, Comsol Multiphysics, STAR-CD, отечественный Логос и другие. Однако данные пакеты программ ещё не обладают требуемым уровнем формализации вычислений, и успешность их применения определяется опытом моделирования конкретного пользователя, который вручную строит вычислительную сетку и выбирает проблемно-зависимые компоненты алгоритма. Использование современных программных средств для моделирования физико-химических процессов основано на ручном построении сетки и расстановке галочек в выпадающих меню в надежде, что собранный таким образом алгоритм будет вычислять достаточно быстро (то есть будет близок к оптимальному). Поскольку сложность современных научных и технических задач исключает возможность их теоретического анализа, то заранее невозможно предсказать продолжительность счёта в последовательном или параллельном исполнениях. Тем не менее, современные программы с дружественным пользовательским интерфейсом и хорошо проработанным постпроцессингом демонстрируют существенный прогресс в области автоматизации математического моделирования физико-химических процессов, хотя за прошедшие 30–45 лет можно было сделать гораздо больше. Тут нет ничего удивительного, поскольку требования к вычислительным алгоритмам, заложенным в программы, устроенные по принципу «чёрного ящика», и к самим программам сформулированы сравнительно недавно [1] .

Математическое моделирование является междисциплинарным направлением, развитием которого занимаются математики, физики, химики, инженеры, профессионалы в области компьютерных технологий, а также узко ориентированные специалисты, причём каждый из них использует привычную для себя терминологию и оперирует привычными для себя понятиями. Поэтому очень важно выработать единую терминологию и систему приоритетов. Например, термин «устроенный по принципу чёрного ящика» у математиков-вычислителей, связанных с математическим моделированием, вызывает ассоциацию с вычислительным алгоритмом, предназначенным для решения систем линейных алгебраических уравнений (СЛАУ), в котором входными данными служат матрица коэффициентов, вектор правых частей и некоторое начальное приближение к вектору неизвестных. Поэтому после продолжительных дискуссий в [1] было предложено следующее определение: «Программное обеспечение для математического моделирования физико-химических процессов в приближении сплошной среды определяют как устроенное по принципу «чёрного ящика» (black-box software), если от пользователя не требуется никаких дополнительных действий, кроме физической формулировки задачи, состоящей из геометрии области, граничных и начальных условий, перечня решаемых уравнений (уравнения теплопроводности, уравнений Навье–Стокса, уравнений Максвелла и так далее) и описания сред. Пользователю не требуется никаких знаний о численных методах, высокопроизводительных и параллельных вычислениях». Из данного определения очевидно, насколько современные программные средства далеки от необходимого уровня формализации вычислений.

Работы советских математиков (метод конечного объёма, А.А. Самарский, 1958; схема Годунова, С.К. Годунов, 1959; многосеточный метод, Р.П. Федоренко, 1961 и другие) на многие десятилетия вперёд определили развитие важнейших направлений вычислительной математики, связанных с численным решением (начально-)краевых задач. Достаточно сказать, что спустя 65 лет высказанные ими идеи не только востребованы, но и составляют вычислительную основу современных программных средств для математического моделирования физикохимических процессов. Более того, дальнейшее совершенствование численных методов, в частности создание универсальных алгоритмов, фактически является эволюционированием результатов, полученных в середине прошлого века советской математической школой.

Программы, устроенные по принципу «чёрного ящика», должны иметь в основе итерационные методы с минимальным количеством проблемно-зависимых компонентов, чтобы избежать затруднений при оптимизации. Другими словами, универсальные алгоритмы должны быть конкурентоспособными по отношению к оптимизированным алгоритмам, но без использования априорной информации о решаемых задачах. Основными требованиями к универсальному алгоритму для программного обеспечения, устроенного по принципу «чёрного ящика», являются: универсальность — минимальное количество дополнительных проблемно-зависимых компонентов, помимо упорядочения неизвестных и критерия останова; эффективность — минимально возможное время решения различных задач в последовательном исполнении; параллелизм реализации, близкий к полному. Кроме того, математическое моделирование физико-химических процессов желательно проводить при минимальном задействовании ресурсов компьютера. Данные требования к универсальному алгоритму образуют проблему «универсальность–эффективность–параллелизм». Идеальный вычислительный алгоритм для перспективного программного обеспечения, устроенного по принципу «чёрного ящика», выглядит следующим образом:

• а)    не содержит проблемно-зависимых компонентов (универсальность);

• б)    обладает оптимальной алгоритмической трудоёмкостью в последовательном исполнении, то есть вычислительная работа (количество арифметических операций) W пропорциональна числу неизвестных N : W ∼ N (эффективность);

• в)    обладает полным параллелизмом (то есть время счёта уменьшается в p раз при решении задач на параллельном компьютере, состоящем из p вычислительных узлов) (параллелизм).

Очевидно, что создать идеальный алгоритм невозможно, поскольку перечисленные требования противоречат друг другу. Например, прямой метод Гаусса с выбором главного элемента для решения СЛАУ является универсальным (не содержит проблемно-зависимые компоненты), но ему не свойственна эффективность: совершается W ∼ N ³ арифметических операций. Метод Якоби обладает полным параллелизмом, но не эффективен и не универсален. Классические многосеточные методы характеризуются оптимальной трудоёмкостью ( W ∼ N арифметических операций), но не универсальны и плохо распараллеливаются на грубых сетках. Аналогичные недостатки можно легко найти у каждого метода решения СЛАУ. Поэтому было предпринято огромное количество попыток развития отдельных свойств алгоритмов за счёт других (например, повышения эффективность за счёт универсальности и другого), но практической пользы от таких изысканий не последовало.

Долгое время не удавалось разработать эффективный алгоритм для численного решения дискретных краевых задач, которые могут быть сведены к СЛАУ с плохообусловленной разреженной матрицей коэффициентов высокого порядка. Первые результаты Р.П. Федоренко, Н.С. Бахвалова, Г.П. Астраханцева, называемые в западной литературе «советским периодом в многосеточных методах», выглядели настолько фантастическими, что научному сообществу понадобилось более 20 лет для осознания их перспективности. В основе многосеточных методов лежит использование быстрой (не зависящей от величины шага сетки) сходимости итерационных методов на первых итерациях. Однако в 80-е годы прошлого века количество исследований в области многосеточных методов резко возросло, и этот период в западной литературе получил название многосеточного бума. Доминирующим направлением развития являлось оптимизирование многосеточных методов: компоненты (сглаживающие процедуры, циклы, операторы пролонгации и сужения, и другое) выбирались исходя из априорной информации о решаемой краевой задаче. Фактически каждая новая краевая задача приводила к новым многосеточным алгоритмам.

В 80-х годах прошлого века началось бурное развитие вычислительной техники и методов решения научных и инженерных задач. В 1988 году Пауль О. Фредериксон (Paul O. Frederickson) и Оливер А. МакБрайан (Oliver A. McBryan) предложили параллельный суперсходящийся многосеточный метод (Parallel Superconvergent Multigrid Method — PSMG:) [3]. Основная идея PSMG заключается в построении двух грубых сеток, состоящих из чётных и нечётных узлов более мелкой сетки. Первоначально цель разработки PSMG заключалась в снижении алгоритмической трудоёмкости многосеточных методов посредством экстраполяции решений с грубых сеток на мелкие. Данный подход трудно использовать для индустриальных задач, но тем не менее Пауль О. Фредериксон и Оливер А. МакБрайан получили два важных теоретических результата:

• 1)    PSMG является односеточным алгоритмом, поскольку число неизвестных на мелкой сетке и двух грубых сетках одинаково. Данное обстоятельство не столь очевидно, однако основополагающую идею Р.П. Федоренко можно с успехом использовать и в односеточных итерационных методах;

• 2)    отсутствие общих узлов на грубых сетках позволяет эффективно распараллеливать сглаживающие итерации. Действительно, сглаживание на грубых сетках можно осуществлять независимо как друг от друга, причём без обмена данными, так и от применяемого сглаживателя. Этим PSMG отличается от классических многосеточных методов, которые трудно распараллеливать на грубых сетках [2] .

В 1990 году в Центральном институте авиационного моторостроения имени П.И. Баранова (Baranov Central Institute of Aviation Motors, Moscow, USSR) была предложена Универсальная Многосеточная Технология — УМТ (Robust Multigrid Technique — RMT) как алгоритм решения широкого класса нелинейных (начально-)краевых задач [4] . В RMT более мелкая сетка представляется в виде объединения трёх более грубых сеток, которые не имеют между собой общих узлов и граней конечных объёмов (то есть посредством агломерации конечных объёмов). В ходе дальнейших исследований установлено:

• 1)    RMT, по сравнению с соответствующим односеточным сглаживателем, содержит единственный дополнительный проблемно-зависимый компонент, а именно количество сглаживающих итераций на каждом сеточном уровне. Однако оптимальный сглаживатель, оптимальное количество сглаживающих итераций и оптимальное упорядочение неизвестных возможно определить в результате прямых вычислительных экспериментов на грубых уровнях и последующей экстраполяции на самую мелкую сетку. Оптимизация проблемно-зависимых компонентов в RMT осуществляется без контроля пользователем вычислительного процесса и приводит к увеличению времени счёта на несколько процентов;

• 2)    минимизация количества проблемно-зависимых компонентов приводит к увеличению алгоритмической трудоёмкости RMT до W ∼ N log N арифметических операций, где N есть число неизвестных, то есть RMT обладает близкой к оптимальной алгоритмической трудоёмкостью;

• 3)    RMT, как и PSMG, является односеточным алгоритмом, который удобнее описывать как многосеточный, и использует основной многосеточный принцип, который заключается в аппроксимации гладкой (длинноволновой) части погрешности на более грубых сетках, а негладкая или грубая часть уменьшается за небольшое количество итераций (независимо от шага сетки) с помощью базового итерационного метода на мелкой сетке. RMT обладает практически полным параллелизмом;

• 4)    самая мощная стратегия отыскания поправок в RMT позволяет применять простейшие сглаживатели для решения научно-технических задач при минимальном использовании ресурсов компьютера, если матрица коэффициентов результирующей СЛАУ не хранится, а вычисляется в процессе решения.

Последовательность реализации RMT такова: сначала строится самая грубая воксельная (трёхмерная, состоящая из 3D-пикселей — вокселей) сетка, которая не согласована с границами области; далее отыскивается решение прямым методом. Потом эта сетка измельчается, осуществляется оптимизация проблемно-зависимых компонентов RMT (сглаживателя, упорядочения неизвестных, количества сглаживающих итераций и другого); снова отыскивается решение и выполняется его апостериорный анализ. В случае необходимости вычислительная сетка локально измельчается. После достижения некоторого критерия останова (например, по количеству узлов или по минимальному шагу сетки), получется достаточно точное приближение к решению и выделены подобласти сгущения сетки. Далее строится вычислительная сетка, согласованная с границами области (body-fitted grid) путём перестроения воксельной сетки с сохранением подобластей структурированности [1] . Затем производится пролонгация решения на воксельной сетке на сетку, согласованную с границами области, и выполняется несколько сглаживающих итераций для удаления погрешностей пролонгации. В общем случае RMT выступает в качестве многосеточного предобуславливателя, а указанная выше процедура получила название метода вспомогательного пространства (Auxiliary Space Method), в котором вспомогательным пространством является воксельная сетка. После завершения сглаживания на сетке, согласованной с границами области, многосеточная итерация считается завершённой. RMT в достаточной мере соответствует всем требованиям, предъявляемым к вычислительным технологиям в программном обеспечении, устроенном по принципу «чёрного ящика».

Следующий вопрос заключается в согласовании математических моделей физико-химических процессов в приближении сплошной среды и вычислительных технологий численного решения представляющих их уравнений. Самым простым подходом стали алгебраические многосеточные методы, предназначенные для решения СЛАУ без использования информации о вычислительной сетке и математической модели. Безусловно, это не самый лучший подход, поскольку ещё Н.Е. Жуковский учил, что «механик должен составлять интегрируемые уравнения». Это означает, что математические модели следует формулировать с оглядкой на методы их решения. Кроме того, алгоритмы решения линейных или глобально линеаризованных систем алгебраических уравнений мало интересны для практических приложений.

Однако прежде, чем приступить к согласованию моделей и вычислительных технологий, необходимо кратко напомнить историю (вычислительной) механики жидкости и газа (МЖГ). В становлении МЖГ можно выделить три характерных периода — аналитический, переходный и вычислительный [5, 6]. Основной задачей механиков в аналитический период (XVII век – первая половина XVIII века) являлся вывод и анализ фундаментальных уравнений, верификация их решений на основе имеющихся экспериментальных данных и применение полученных результатов к решению инженерных задач (совершенствование речных и морских судов, различных гидротехнических сооружения и дугого). Основой МЖГ стала гипотеза сплошности, согласно которой жидкость или газ заменялись непрерывно распределенной по пространству фиктивной сплошной средой (континуумом), обладающей в макромасштабах физическими свойствами реальной жидкости, и фундаментальные законы сохранения массы, импульса и энергии. Казалось, проще всего разделить исходную область на малые подобласти и для каждой из них записать законы сохранения. Однако даже для малого количества подобластей результирующие уравнения могли быть столь громоздкими, что их было трудно записать гусиным пером в мерцающем свете восковых свечей. Поскольку молекулярное строение жидкостей и газов в то время практически оставалось не изученным, то континуум был доопределён свойством бесконечной делимости в микромасштабах. В результате к каждой малой подобласти стало возможным применить операцию предельного перехода (фактически стянуть подобласть в точку) и записать дифференциальные уравнения, уже не зависящие от формы подобласти. Сначала были сформулированы уравнения Эйлера, а затем, по мере изучения процессов теплопроводности и вязкого трения, и уравнения Навье–Стокса. Далее, в простейших случаях получены точные решения уравнений Навье–Стокса (УНС), сформулированы приближённые теории (теория пограничного слоя [7]) и новые фундаментальные проблемы (существование и единственность решений УНС — задача тысячелетия или Millennium Prize Problems).

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

p ( v )=MP

.

Дифференциальное описание физико-химических процессов подразумевает, что каждая функция определена и дифференцируема в каждой точке. Операция предельного перехода (стягивание подобласти V в точку x ∈ V ) лишает плотность физического смысла:

( У r M (V ) p(x)= lim----- ^P( ) V m o V

Аналогично теряют физический смысл давление, температура, концентрация и другие макропараметры.

Операция предельного перехода добавляет макропараметрам новые, не характерные для них свойства. Рассмотрим два соседних объёма V ₁ и V ₂ , которые имеют общую грань. В каждом из них можно рассчитать плотность традиционным образом как p 1 = M 1 (V 1 )/V 1 и p ₂ = M 2 (V 2 )/V 2 . Давление можно представить как полевую переменную, постоянную внутри каждого объёма и разрывную (а точнее, неопределённую) на его гранях, поскольку геометрическая размерность грани на единицу меньше геометрической размерности объёма, и плотность на грани нельзя выразить как отношение массы к «объёму» грани. В результате операции предельного перехода «плотность» получается в виде достаточно гладкой функции, что противоречит гипотезе сплошности.

Правильнее всего, было бы признать приближённый характер дифференциального описания физикохимических процессов в приближении сплошной среды. Однако недостатки дифференциальных УНС пытались маскировать при помощи бесконечно малого физического объёма (БМФО), содержащего достаточное количество частиц, чтобы считать среду сплошной: если известно значение некоторого макропараметра в некоторой точке БМФО, то значение того же макропараметра во всём БМФО будет приблизительно таким же. Иногда в литературе такой подход называют БМФО как «точка» (Рис. 1) .

Появление в середине прошлого века компьютеров ещё больше запутало ситуацию. Вычислительная техника не оперирует с дифференциальными уравнениями, поэтому стало необходимым заменить краевую задачу (КЗ) или начально-краевую задачу (НКЗ) приближённой алгебраической задачей. Эта замена называется аппроксимацией и является отдельной проблемой. Роль дифференциальных уравнений в математическом моделировании лучше всего представлена в [8]: «в реализуемой цепочке традиционного моделирования дифференциальное и интегральное исчисление

БМФО как «точка»

^Qz = ^hz

z

^Qx =0

модель:

КЗ или НКЗ

аппроксимация итерационные методы

модель: СНАУ

итерационные методы: Зейдель+Ньютон

численное решение

численное решение

Рис. 1. Дифференциальный и алгебраический подходы к математическому описанию физико-химических процессов в приближении сплошной среды

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

Первые же компьютерные расчёты показали существенную ограниченность использования дифференциальных уравнений для исследования отдельных классов течений. Актуальная проблема распространения ударных волн (разрыва параметров внутри области) была решена фактически без употребления дифференциальных уравнений [9, 10] . Помимо этого, ставшая легендарной схема С.К. Годунова позволяла производить расчёты без явного выделения особенностей решения в предположении постоянного значения параметров в каждой ячейке и убедительно демонстрировала преимущества тщательного учёта физических аспектов моделируемых процессов над математическим формализмом при выражении потоковых переменных через консервативные переменные. К сожалению, идеи С.К. Годунова вовремя не получили должного развития.

В вычислительной макромеханике БМФО является не точкой, а минимальным объёмом, размеры которого таковы, что заключённая в нём среда считается сплошной, но бесконечной делимости этой среды уже не требуется (Рис. 1) [5, 6] . Минимальные размеры конечных объёмов ограничены именно требованием сплошности среды и не позволяют выполнить операцию предельного перехода, которая приводит к производным. С точки зрения физики, каждый конечный объём является открытой термодинамической системой, находящейся в термодинамическом равновесии, поэтому все макропараметры (плотность, давление, температура, концентрации и другие) являются постоянными внутри каждого объёма. Гипотеза сплошности, которая ограничивает минимальные размеры конечных объёмов, и гипотеза локального термодинамического равновесия, которая обуславливает максимальные размеры конечных объёмов, дают возможность решать широкие классы научных и инженерных задач без привлечения дифференциального исчисления. Сходные идеи высказываются достаточно давно: «Другими словами, мы можем получить алгебраическую формулировку, избегая произвольного процесса дискретизации дифференциальных уравнений. Эта формулировка обладает большим достоинством, поскольку сохраняет тесную связь между математическим описанием и описываемым физическим явлением» [11] , но широкого распространения они пока не получили. Математические модели исследуемых физико-химических процессов, состоящие из законов сохранения макропараметров, феноменологических уравнений и уравнений состояния, строятся в каждом объёме, который по предположению находится в состоянии термодинамического равновесия. В этом случае понятия «математическая модель» и «разностная схема» совпадают.

Переходный период становления МЖГ (вторая первая половина XX века – настоящее время) характеризуется сменой математического аппарата: УНС решаются приближённо, численными методами, при помощи быстродействующих компьютеров, то есть произошла формальная смена решателя. (Начально-)краевая задача для УНС заменяется более простой — алгебраической — задачей, погрешность замены является контролируемым параметром. Численное решение УНС породило множество математических проблем: построение (адаптивных) сеток в областях со сложной геометрией; выбор метода и порядка аппроксимации (начально-)краевых задач; теоретический анализ существования и единственности решения дискретных аналогов УНС; параллельное решение седловых задач с оптимальными усилиями и так далее. Другими словами, к существующим проблемам при оперировании (интегро-)дифференциальными уравнениями добавились новые проблемы, связанные с переходом к их дискретным аналогам. Фактически современная вычислительная физика превратилась в раздел вычислительной математики, связанный с формальным решением (интегро-)дифференциальных уравнений всё более и более точными численными методами без учёта физической природы моделируемых процессов. В переходный период не было понято самое главное: использование принципиально новых методов решения неизбежно приводит к необходимости пересмотра фундаментальных уравнений.

По мере усложнения моделируемых систем стала очевидной потребность изменения существующих представлений об уравнениях вычислительной гидрогазодинамики. В 50-е годы прошлого века в связи с развитием оружия массового уничтожения возникла потребность в расчёте течений с ударными волнами, а формальная замена производных конечными разностями не позволяла достигнуть требуемой точности. Дж. фон Нейман одним из первых предложил использовать интегральные тождества и энтропийное неравенство для построения разностных схем для уравнений газовой динамики вместо традиционной аппроксимации дифференциальных уравнений в областях, где функции, описывающие течение, являются достаточно гладкими [10] . Таким образом возникла дилемма: либо корректное определение производных и некорректное определение макропараметров (БМФО как «точка»), либо корректное определение макропараметров и конечные разности вместо производных при явном ограничении минимальных шагов сетки (БМФО как объём) (Рис. 1) .

Третий (вычислительный) период связан с объединением фундаментальных уравнений МЖГ и методов их численного решения в единую физико-математическую методологию, в совокупность подходов, приёмов и методов теоретической физики и вычислительной математики, являющуюся целостностной технологией численного изучения физико-химических процессов. В [5] для данной методологии автором предложено название вычислительная макромеханика, где термин «вычислительная» означает решение физико-химических задач методами вычислительной математики без использования интегрального и дифференциального исчисления, а термин «макромеханика» подразумевает явное наличие макромасштабов при физико-математическом представлении течения жидкости или газа в приближении сплошной среды.

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

• 2.    Постулаты вычислительной макромеханики

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

Постулат 1 . Для каждого физико-химического процесса существует минимальный пространственный (ℏ _s ) и временной (ℏ _t ) масштабы, при которых фундаментальные уравнения, полученные в приближении сплошной среды, имеют физический смысл (гипотеза сплошности). Субпроцессы, протекающие при меньших масштабах, могут быть учтены только при помощи феноменологических законов. Именно наличие минимальных масштабов не позволяет стягивать конечные объёмы (КО) в точки и описывать физико-химические процессы при помощи частных производных.

Постулат 2 . Каждый физико-химический процесс протекает в пространственной области П за время [0,T] . Данная область Л х [0,T] представима в виде объединения КО П _т х (t (n⁺¹⁾ -t ⁽ⁿ⁾ ) так, что

MN

Л X [0,T]= [ [ П т X (t ⁽ n ⁺¹⁾ — t ⁽ n) ), m=1n=1

где t ⁽¹⁾ =0 и t N ⁺¹⁾ = T , причём min dim П _т >  h s и t ⁽ⁿ⁺¹⁾ — t ^{( n )} ^ h _t .

Замечание 1. В задачах в ^m ычислительной МЖГ минимальный пространственный размер каждого КО Ω _m должен быть таким, чтобы жидкость или газ в нём можно было считать сплошной средой. Например, согласно молекулярно-кинетической теории, критерием сплошности газа в каждом КО П _т служит число Кнудсена ( Kn s ), равное отношению средней длины свободного пробега молекул l к минимальному пространственному размеру ℏ s :

Kn _s =1/^.

При Kn _s <  10 ^{- 3} газ в каждом КО П _т можно считать сплошной средой, дифференциальным УНС соответствует Kn _s = +м .

Аналогично необходимо соблюсти условие малости числа Кнудсена, определяемого отношением среднего времени между столкновениями молекул газа t _m к минимальному шагу по времени: min (t ⁽ⁿ⁺¹⁾ — t ⁽ⁿ⁾ ) , то есть n

Kn t = t m /min (t ⁽ⁿ⁺¹⁾ — t ⁽ⁿ⁾ ) <  10 ^{- 3} .

З а м е ч а н и е 2. Граница сплошности среды весьма условна, и не существует чётких критериев сплошности для газов, жидкостей и твёрдых тел.

Постулат 3 . Каждый КО считается открытой термодинамической системой, находящейся в термодинамическом равновесии. Это позволяет полагать макропараметры (плотность, давление, температуру и другие) постоянными внутри КО в соответствии с законами равновесной термодинамики и разрывными на его границах. Все макропараметры рассматриваются как функции положения в пространстве и времени и могут быть закреплены в любой внутренней точке КО для вычислений. Любые равновесные соотношения в каждом КО могут использоваться для моделирования неравновесных физико-химических процессов в приближении сплошной среды (гипотеза локального термодинамического равновесия).

З а м е ч а н и е . Значения макропараметров на гранях КО строго определить нельзя, поскольку геометрическая размерность грани на единицу меньше геометрической размерности объёма. Пример: пусть кубический КО заполнен газом, гранью куба ( d _v = 3 ) является квадрат ( d _f = 2 ), поэтому нельзя найти температуру газа на грани тем же статистическим способом, как и температуру газа в кубе. Температуру газа на грани КО можно установить из условия равенства тепловых потоков через общую грань соседних КО.

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

Cопоставим основные положения вычислительной макромеханики со схемой С.К. Годунова [9] :

• 1)    Схема С.К. Годунова основана на аппроксимации методом конечного объёма и выражении потоковых переменных

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

• 2)    Использование дискретных уравнений в схеме С.К. Годунова снижает требования к гладкости решений и позволяет рассматривать задачи с разрывными и непрерывными решениями унифицированным образом. Аналогично в вычислительной макромеханике унифицированным образом можно решать задачи с непрерывными и разрывными решениями и/или граничными условиями.

• 3)    В исходной схеме С.К. Годунова значения функций считались постоянными в каждом конечном объёме. В вычислительной макромеханике постоянство макропараметров внутри каждого конечного объёма обусловлено гипотезой локального термодинамического равновесия.

• 4)    Традиционно построение схемы С.К. Годунова начинается с интегрирования дифференциального уравнения переноса по конечному объёму, хотя дискретные уравнения можно получить из законов сохранения макропараметров. В вычислительной макромеханике тоже необходимо прибегать к обоим подходам при выводе основополагающих уравнений, хотя первый подход позволяет использовать преимущества дифференциального исчисления для преобразования дифференциальных уравнений (например, сводить конвективно-диффузионное уравнение к чисто диффузионному уравнению), но при формальном интегрировании важно помнить, что все макропараметры постоянны в каждом конечном объёме, находящемся в состоянии локального термодинамического равновесия.

• 3.    Стационарное течение в каверне с движущейся крышкой при малых числах Рейнольдса

  Одна из самых распространённых тестовых задач стационарного течения в двухмерной каверне с движущейся со скоростью U _w крышкой может послужить для демонстрации преимуществ вычислительной макромеханики. Достаточно сказать, что первые результаты моделирования вихревого течения в каверне получены ещё в 1961 году [12] . Для наглядности течение будем считать настолько медленным, что конвективным переносом можно пренебречь, и математическая модель вихревого течения в каверне будет включать следующие безразмерные уравнения:

Таким образом, вычислительную макромеханику можно рассматривать как развитие идей С.К. Годунова и обобщение на задачи не только с конвективным, но и с диффузионным переносом.

• а) уравнение неразрывности

∂u ∂v

F + F = 0;

∂x ∂x

• б)    уравнение движения по x

1 ( Q Uu

∂ 2 u

Re \ dx² + dy ²

∂p dx ’

• в)    уравнение движения по y

1 d d Vv

∂ ² v

Re \ dx² + dy ²

∂p ∂y.

Геометрия области и граничные условия показаны на рисунке 2.

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

Рис. 2. Течение в каверне с движущейся крышкой ( а ) и разнесённая сетка ( б )

А и Б (Рис. 2) не позволяет получить классическое решение краевой задачи для дифференциального уравнения (1). Причиной данного эффекта являются повышенные требования к гладкости решений и граничных условий, обусловленные операцией предельного перехода. При стандартном построении разнесённой сетки [13] в случае, когда линии сетки параллельны основанию каверны и крышке, удаётся выколоть точки разрыва и численно решить дифференциальную задачу, которая не имеет аналитического решения. Однако отсутствует сходимость этого численного решения к аналитическому, если устремить шаг разнесённой сетки к нулю.

Рис. 3. Разнесённая сетка в каверне с движущейся крышкой

Следуя С.К. Годунову, перейдём к дискретному описанию течения в каверне, понизив требования к гладкости решений и граничных условий. Точная аппроксимация граничных условий Дирихле (прилипания и непротекания) подразумевает расположение граничных узлов сетки на границе области. Для этого следует повернуть сетку на 45 ° относительно каверны, как показано на рисунке 3.

Использование разнесённой сетки для численного решения дискретных краевых задач подразумевает, что для каждого уравнения математической модели используется свой КО (Рис. 4) . Вертикальные и горизонтальные линии сетки (как объединение граней КО) являются линиями разрыва соответствующих компонент вектора скорости. Кроме того, граничные условия в вычислительной макромеханике задаются не в точках, а на КО, пересекающихся с границей расчётной области. Методика вывода основополагающих уравнений МЖГ в дисконтинуальном приближении описана в [5] , поэтому здесь приведены только окончательные безразмерные уравнения:

• а)    закон сохранения массы

u i+ij u ij ₊ v ij+1 v ij = ₀. hh

• б)    проекция закона сохранения импульса на ось i

1 u i - 1j 2u ij ⁺u i+1j    1 u ij - 1 2u ij ⁺u ij + 1 _ pij p i - ij _ о

Re       h ^{2      +}Re      h           h =

• в)    проекция закона сохранения импульса на ось j

• ¹    v i - 1j ^2v ij ⁺ v i+1j . ¹ v ij - 1 ^2v ij ⁺ v ij + 1    p ij p ij - 1

Re         h ²         ' Re         h ²               h

Здесь u , v и p есть компоненты скорости и давление; h >  hs — шаг равномерной сетки; Re — число Рейнольдса.

Рис. 4. Граничные условия в вычислительной макромеханике

v = 0

v = 0

Численное решение системы (2) – (4) осложняется отсутствием уравнения для отыскания давления p _ij . В сегрегированных (раздельных) алгоритмах (например, в наиболее популярном методе SIMPLE) для отыскания давления, как правило, используется уравнение Пуассона, установленное некоторым образом. Но, если дифференциальное уравнение Пуассона ещё можно как-то получить, то корректно сформулировать граничные условия (Неймана) для него не представляется возможным [4] . Данные о дополнительных граничных условиях для давления в исходной постановке задачи отсутствуют. Хотя уравнения Стокса достаточно хорошо изучены [14] , совершенно непонятно различие между решениями исходной задачи и задачи с дополнительными граничными условиями для давления.

В настоящее время одна из самых значительных и популярных тем в вычислительной математике — это исследование сеточных седловых задач, имеющих выраженную блочную структуру. Это объясняется широтой приложений и новизной идей, необходимых для разработки эффективных методов их решения [14] . Наиболее перспективными (особенно для мультифизичных приложений) являются совместные итерационные алгоритмы на основе метода Ванки [15] . Перепишем систему (2) – (4) в виде:

u i+1j - u ij +v ij + 1 - v ij = ⁰,

- ⁴u ij + u i-ij +u i+ij +u ij — i +u ij+1 - p ij +P i - 1j =⁰, ^-4v ij ⁺ v i — 1 j +v i+1j ⁺v ij - 1 +v ij+1 ^-p ij ^+p ij - 1 = ^0.

Здесь p ij = Rehp ij . Следуя методологии метода Ванки, придадим полученным уравнениям матричную форму:

Г4

-4

\-1

-4    1

1 -4

-1    1

• -1      u ij

¹     u i+1j

• ^{- 1}       v ij

1     v ij +1

0      p ij

b 1

b 2

^b 3

b ₄

где компоненты вектора неизвестных есть b1 = -ui-1j -uij-1 -uij+1 -pi-1j , b2 = -ui+2j-ui+1j-1 — ui+1j+1 + pi+1j , b3 = -vi-1j -vi+1j -vij-1 -pij-1, b4 = -vij+2 — vi-1j+1 -vi+1j + 1 + pij + 1.

Основная идея метода Ванки состоит в совместном решении уравнений (3) (первые два уравнения системы (5) ),

(4) (третье и четвёртое системы (5) ) и (2) (последнее уравнение системы (5) ). Поскольку система (5) состоит

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

Рис. 5. Распределение давления в каверне с движущейся крышкой

Систему уравнений (5) можно решить точно при помощи метода Гаусса с выбором главного

элемента:

uij ui+1j vij vij+1 pij J

^{1 - 13}

-7

-3

\-15

^{- 3} \ 3

-7

-13

lb' \ b2

^b 3

b 4

Таким образом, с каждым КО связана система линейных алгебраических уравнений типа (5) , решение которой позволяет определить компоненты скорости u _ij , u _{i +1 j} , v _ij , v _{ij +1} на гранях КО и давление p ij в центре КО, причём компоненты скорости точно удовлетворяют закону сохранения массы (2) . Итерация метода Ванки состоит в последовательном решении систем типа (5) , каждая из которых связана с определённым КО. К достоинству метода Ванки следует отнести отсутствие проблемно-зависимых релаксационных параметров, поэтому метод Ванки часто используется в качестве сглаживателя в многосеточных методах.

Для наглядности результаты вычислений на сетке с h =1/9 и Re= 1 , показанной на рисунке 3, представлены на рисунке 5. Столь грубая сетка необходима для иллюстрации «растрового» характера изменения давления, обусловленного постоянством на КО. В итерациях Ванки давление не фиксируется в каком-либо узле сетки, но при визуализации полагается, что давление в левом углу рассматриваемой единичной каверны равно нулю.

Как видим, обобщение идей С.К. Годунова на дозвуковую МЖГ приводит к тому же результату — чисто алгебраическому (без аппроксимации дифференциальных уравнений) описанию физико-химических процессов физически значимыми макропараметрами. Снижение требований к гладкости решений позволяет моделировать процессы с разрывными решениями и разрывными граничными условиями.

• 4.    Дискуссия: вычислительная математика и вычислительная физика

С учётом накопленного опыта математического моделирования физико-химических процессов и сегодняшнего уровня развития численных методов представляется целесообразным классифицировать работы следующим образом:

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

– к вычислительной физике следует отнести исследования физических задач, решаемых численными методами при помощи вычислительной техники. Все константы, переменные, параметры и функции должны обладать чётким физическим смыслом, за исключением специально оговорённых случаев (например, турбулентной вязкости при моделировании турбулентности в приближении RANS). Основная цель вычислительной физики — как можно более точное физико-математическое описание исходного объекта или физико-химических процессов и последующее компьютерное изучение основополагающих уравнений для обобщения полученных результатов на объект или процесс. В вычислительной макромеханике дискретные уравнения являются точными, а дифференциальные — приближёнными.

Для наглядности рассмотрим на единичном отрезке [0,1] краевую задачу для одномерного уравнения теплопроводности:

(6a)

(6b) (6c)

q‘ = -f, с граничными условиями

T (0)=0, q(1)=0, где q = XT, есть удельная плотность теплового потока, A(x) — известный коэффициент теплопроводности, а T — неизвестная температура. Задача исследования состоит в сравнении точных решений дифференциальной и дискретной задач теплопроводности (то есть распределений температуры) и анализе различий между ними.

Повторное интегрирование (6a) в пределах [0,x] c учётом граничных условий (6b) и (6c) позволяет получить на отрезке [0,1] функцию температуры T (x) :

x1

T (x)= /¹ Л«)d^. о ^Л( ^

Температура традиционным образом известна в каждой точке области [0,1], что противоречит её определению как макропараметра.

Построим на единичном отрезке [0,1] равномерную сетку с шагом h и узлами xV. Значения температуры в узлах сетки xiv составят:

i-1 ^{x v n +1}

T ^{(x V )} =X / ж d n=1 _v

xn где функция F (£) имеет вид

F ® = Jf (С М.

ξ

Точное решение (7) краевой задачи (6) получено без учёта гипотезы сплошности.

Теперь рассмотрим на единичном отрезке [0,1] точную краевую задачу для дискретного 1D уравнения теплопроводности:

qh(xf )-qh(xf-i) =-hsfh(xv),                                       (8a) с граничными условиями T h(xV)=0,                                            (8b) qh(xfO=0,                                       (8c) где qh(xf ) = Ah(xf)T h(xv+i)-T h(xv)                                   (8d) i            i              ℏs есть удельная плотность теплового потока, λh — коэффициент теплопроводности, а ℏs — минимально возможный шаг равномерной сетки, при котором среду в КО сетки можно считать сплошной. Шаг вычислительной сетки не может быть меньше hs: h > hs, при этом ^ = 1/h есть параметр дискретизации. Ограничение на величину шага сетки (h ⩾ ℏs) обусловлено гипотезой сплошности.

В вычислительной макромеханике краевая задача, сформулированная при минимально возможном шаге сетки h = h s , называется точной, а её решение — точным решением. Ограничение h = h s есть следствие справедливости гипотезы сплошности.

Напомним разницу между функциями q , λ , T , f в задаче (6) и q ^h , λ ^h , T ^h , f ^h в задаче (8) : функции в задаче (6) являются достаточное количество раз дифференцируемыми, а функции в задаче (8) — кусочно-постоянными, то есть они постоянны в пределах КО и разрывны на их гранях. Для удобства будем использовать значение функции в точках сетки, полагая его постоянным для всего открытого КО, например:

f h ⁽x ^{v )} = f = ^{f (x) I x f-i f }.

Суммирование (как дискретный аналог интегрирования) дискретного уравнения теплопроводности (8a) приводит к следующему равенству:

ℵ+1                                            ℵ+1

£ ( q⁽x m ^{) -} q^(x m-i ⁾ ) =qCxUOqCx f-O ^{- h} s ^^{f (}x m ).

m=i                                          m=i

Отсюда с учётом граничного условия (8c) нетрудно получить следующее уравнение:

ℵ+1

qO-iO hs X f CO m=i или, согласно определению (8d), в виде:

ℵ+1

T (x n ) - T (x n — i )= h 2   f     X f (x m ).

Mxn -1 ) m=n

Суммирование (как дискретный аналог интегрирования) данного уравнения i                                 i

X ( T (xO TOi) ) = h 2 X T7M X f (x vm )

n=2                      n=2 xn-1 m=n с учётом граничного условия (8b) позволяет рассчитать значения температуры в узлах сетки:

i

T h(xV)= % X         X f h(xm), i = 2,3,..O1.

n=

Уравнение (9) получено в предположении, что гипотеза сплошности справедлива при условии h ⩾ ℏ _s .

Теперь остаётся сравнить точные решения задач (6) и (8) , чтобы оценить влияние границ применимости гипотезы сплошной среды на точность получаемых результатов. Разность температур (7) и (9) составляет:

xv n                 2    ℵ+1

d^ - w s a Ef ^h (x m ) | .                   (10)

j A(e)   ООО о / xn-1

Полагая, что F (х)/А(х) является достаточно гладкой функцией, и используя разложение в ряд Тейлора в окрестности точки х П _{- 1} , приходим к выражению:

F (X ₌ Х£ d k F (X         f _{}к =} ^F (х П- i ) ,X£ d k F (€)

А(Х k=₀ k! de А(е x n- 1     n ^{- 1} ) А(х П-1 ) ^k! d£ k А(е

x

f

n

-

(i-х П-1 ) ^к ,

при этом грань х П _{- 1} расположена посередине между узлами х П _-1 и х П :

x

x n-1

= x n—1 ^{- х} П-1 = ^ s /2.

Тогда первый член в левой части (10) можно представить следующим образом:

x _n

/ F1UЙ F(х П - 1 ) +v 1 d F^

J X(e ^ s A(4- 1 )+£l k! de О

^ v                                    k 1

x n- 1

x

f

n

¹ x ^v

I ^(e

n

^- x n-1 ^{) k} de

Интеграл

b

J = I (x - cy m dx.

a

b — c = c - a = e/^’

равен:

b

J ^-cy m d =

a

0, h m⁺¹

для чётных m,

2 ^m m+1 ’

для нечётных m,

то есть

v xn

J < - xn xv n-1

1 Й 2l + 1 ) 2l de=1      ’

^{11    4}   4 l 2 l + 1 ’

l = 0’1’2,3’... .

С учётом этого выражение (11) принимает окончательный вид:

^v x _n

Z

v xn-1

ь F (х П - 1) . X 1          d ² l F (i)

А(Х ^ s А(х П-1 )⁺   4 l (2l +1)! di ² A(e

f n-1

Воспользуемся разложением в ряд Тейлора в окрестности точки х П _{- 1}

1          К+1 x m          К+1 x m ■ ^

.ечхе l!

F (х П-1 )= / f (Z )dZ = £ / f (x)dx = £ / £ _f                m-n _f                 m-n _f l-0

(x — xm^dx

для преобразования отношения F (х П _{- 1} ) /А(х П _{- 1} ) . Вычисляя здесь интеграл по аналогии с предыдущим, получаем выражение

^{К + 1} /            ⁺ “            1 fc2l + 1 \

F (х П-1 )=Е м (х т )+Еf ^(2l) (х т );! X    ’

4 ^l 2l+1

m-n              l-1

из которого вытекает:

F ^(х П-1 )

^А(х П-1 )

f   (х т )+

^А(х п-1 ) m-n

^А(х П-1 )

К+1 +^

EE f ⁽²¹⁾ (х т ) m-n l-1

1 ^l⁺² 4 l 21 + 1.

И выражение (12) принимает вид:

x

/

x

v n

F^= аю ^e

- 1

Ь 2     ^{К + 1}               1     ⁺“ 1 t 2l+1

X       X s

А(х П-1 )^nf^{( m} ' А(х П-1 ) l— 4 l (2l +1)!

К + 1

е X f ⁽²¹⁾ (х т )+

d ² f (е di ²⁷ Ж

После подстановки этого результата в (10) разность температур, являющихся точными решениями задач (6) и (8) , принимает вид:

T(x ^v)- h^h(x^v} X__ ¹__ X J1X( X fC²¹^}     F ⁽^⁾

⁽ i ^-    ( i ^- n=₂M4-1 ) 4^l (2l + 1)!    s^   ⁽    ' d?l A(^

Положим, что все производные ограничены max I f (2l)(x) I < Cf, xe[0,i]'                1 f

max xe[0,i]

d²¹ F (x)

dx ^{2 1} A(x)

⩽ C F

и рассмотрим первый член разложения (l = 1 ):

3    h ³            1

max I T ⁽x ^{v )}-T ^{h (}x x ^{v )} | < .^ X A f )

(MK+1- n)c _f + CF) .

Поскольку

М^ + 1-п)С > <  eC C _} = C _f ,

i

X

1 ^-J

i - n         ℵ

min A(x) min A(x)

x£[0,1]            xG[0,1]

1       1

h min A(x), xe[0,i]

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

max | T (x V )- T ^h (x V ) | i

^ C f + Cf ^

min A(x) 24

xe[0,i]

В вычислительной математике дифференцируемая достаточное количество раз функция T (x V ) (7) считается точным решением, а дискретная функция T ^h (x V ) (9) — приближённым решением, причём

Th(xv) ^T(xv) при С ^ °, то есть при отсутствии ограничений на величину шага сетки, накладываемых гипотезой сплошности. В индустриальных приложениях справедлива оценка погрешности приближённого решения T h :

max 1T (xv)- T h(xv) 1 < Ch2, где h ≫ ℏs (Рис. 6).

вычислительная математика

(Н)КЗ ⁰                                       J

1                      h_s               СХЕМА ( h » h_s )

ФПФ h ® 0)                                 *

0       вычислительная макромеханика

Рис. 6. Два подхода к математическому описанию физико-химических процессов

В вычислительной макромеханике существует не одно, а два приближения к решению: первое — формальное, может быть получено в пределе h _s ^ 0 с нарушением гипотезы сплошности. Формальный предельный переход h _s ^ 0 приводит к формально-предельной форме (ФПФ) дискретной краевой задачи (8) , а именно к дифференциальной краевой задаче (6) (Рис. 6) . С учётом гипотезы сплошности погрешность аналитического решения составит:

max | T (x v ) -T h ^ ) | < C%,

i иначе говоря для достаточно гладких решений различие между решениями пренебрежимо мало на сетках с предельно малым шагом ℏs .

С другой стороны, решение на сетках с предельно малым шагом ℏs не всегда оправдано, поэтому для второго приближения при h ≫ ℏs следует проанализировать погрешность решения на грубых сетках. Рассмотрим изменение погрешности при аппроксимации на многосеточной структуре RMT (иерархии структурированных сеток, построенных путём утроения шага), порождаемой самой мелкой сеткой с шагом ℏs . Закон сохранения энергии на многосеточной структуре принимает вид:

γ qf )-q(xf{i-i})=-hs 12 f (xb+k), k = -Y где параметр γ определяется как

Y =(3l -1)/2, l есть номер сеточного уровня, а {i} — отображение индексов граней и узлов сеток. Суммирование, как аналог интегрирования, позволяет получить следующее выражение для удельной плотности теплового потока:

* l +1 Y

q(x{i-i})=^sXX f(x{n}+k), n=i k=-y при этом ty — параметр дискретизации сеток уровня l. Отсюда нетрудно придти к точному решению краевой задачи (8) на многосеточной структуре:

i                * i ^{+ 1} Y                            i                   * + 1

T h ( x {i} ) = ^l X    f XX f h ( x {n}+k ) = ^l X    f        X f h ( x n ) .

m=2 ' ^{{ m — 1}}/ n=mk= - Y                   m=2 ' ^{ m ^{- 1}}' n={m}- y

Для самой мелкой сетки ( l = 0 , {i} = i ) справедливо соотношение:

i       ,      N+1

T h (x {i} )= ^ 2 X      f X f h (x n ).

m=2    m ^{- 1} n=m

Тогда погрешность численного решения на сетках уровня l составит:

max|Th(xvJ-T h(xv JI < 3lxmf( max A(x) - min A(x)), i 1        {i}         v {i}           max A2(x) i£[c,1]          xe[0,1]

x€[0,1]

то есть погрешность пропорциональна 3 l , где l есть номер сеточного уровня.

• 5.    Заключение

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

Вычислительную макромеханику можно рассматривать как развитие и обобщение идей С.К. Годунова на задачи не только с конвективным, но и с диффузионным переносом.

Работа выполнена при поддержке Российского научного фонда, проект № 23-19-00734 .