В 1910 г. Л. Ричардсон представил в Королевское общество пятидесятистраничную статью, которая должна быть признана краеугольным камнем численного анализа дифференциальных уравнений в частных производных. До этого Шепперд выполнил некоторую фундаментальную работу по конечно-разностным операторам, однако вклад Ричардсона затмил все предыдущие исследования. Ричардсон разработал итерационные методы решения уравнения Лапласа, бигармонического уравнения и других уравнений. Он установил различие между стационарными задачами в зависимости от того, «можно или нельзя продолжить решение, отправляясь от некоторой части границы», т. е. в современной терминологии различал гиперболические и эллиптические задачи. Ричардсон тщательно изучил численное задание граничных условий, включая граничные условия в угловой точке и на бесконечности. Он получил оценки погрешности, дал метод экстраполяции полученных результатов при стремлении шага сетки к нулю, а также предложил проверять численные решения сравнением с точными решениями для тел простой формы, скажем для цилиндра. Наконец, он впервые фактически применил численные методы к такой практической задаче большого масштаба, как определение напряжений в каменной дамбе.
В итерационном методе Ричардсона для эллиптических уравнений на п-й итерации поочередно в каждом узле расчетной сетки удовлетворяется конечно-разностное уравнение, содержащее «старые» значения на (п-1)-й итерации в соседних узлах. В 1918 г. Либман показал, что можно значительно увеличить скорость сходимости просто за счет использования «новых» значений в узлах, как только они вычислены. В этой схеме «непрерывных замещений» на каждой n-й итерации используется некоторое число старых значений с (n-1)-й итерации и некоторое число новых значений с n-й итерации в соседних узлах. В каждом цикле итерационного метода Либмана наибольшие погрешности уменьшаются так же, как в двух циклах итерационного метода Ричардсона (Франкел [1950]).
Это сравнение служит примером специфики численного анализа уравнений в частных производных. Оказывается, что небольшое изменение конечно-разностных аппроксимаций, итерационных схем или трактовки граничных условий может дать большой выигрыш. Напротив, некоторые правдоподобные и на первый взгляд точные численные схемы могут приводить к полной катастрофе. Классическим историческим примером здесь является явная схема Ричардсона для параболического уравнения теплопроводности, в которой использовались конечно-разностные аппроксимации производных центральными разностями как по пространственным переменным, так и по времени. О'Брайен с соавторами [1950] показал, что эта схема безусловно неустойчива.
До появления ЭВМ основное внимание уделялось эллиптическим уравнениям. Первое строгое математическое доказательство сходимости и оценку погрешности итерационного метода Либмана для решения эллиптических уравнений дали Филлипс и Винер [1923]. В 1928 г. появилась классическая работа Куранта, Фридрихса и Леви. Эти авторы в основном интересовались использованием конечно-разностных методов как инструмента для исследований в чистой математике. Дискретизируя дифференциальные уравнения, доказывая сходимость дискретной системы к дифференциальной и, наконец, устанавливая существование решения дискретной системы алгебраическими методами, они доказывали теоремы существования и единственности для эллиптических, гиперболических и параболических систем дифференциальных уравнений. Эта работа определила направление практического получения конечно-разностных решений в последующие годы.
В СССР метод конечных разностей для уравнений в частных производных систематически начал разрабатываться в 30-х годах. К числу первых исследований в этой области относятся работы С.А.Гершгорина [1930], Л.А.Люстерника [1934], Д.Ю.Панова [1932-1933], И.Г.Петровского [1941] и др.
Первое численное решение уравнений в частных производных для задач гидродинамики вязкой жидкости было дано Томом в 1933 г. В 1938 г. Шортли и Уэллер разработали метод, являвшийся, по существу, более сложным вариантом метода Либмана. Они предложили блочную релаксацию, метод пробной функции, релаксацию погрешности, методы измельчения сетки и экстраполяцию погрешности. Они также впервые точно определили и исследовали скорость сходимости.
Саусвелл [1946] разработал более эффективный метод релаксации для численного решения эллиптических уравнений. В его методе релаксации невязки не проводятся вычисления последовательно в каждом узле сетки, а просматривается вся сетка для нахождения узлов с максимальными невязками, и именно в этих узлах вычисляются новые значения. (В случае стационарного уравнения теплопроводности невязка пропорциональна скорости накопления энергии в ячейке сетки; следовательно, стационарное состояние достигается, когда все невязки обращаются в нуль.) Фокс [1948] разработал усложненные варианты метода релаксации Саусвелла, введя схемы верхней и нижней релаксации (при которых невязки не полагаются точно равными нулю), способ выбора узлов сетки, в которых осуществляется релаксация, а также схему блочной релаксации.
В 1955 г. Аллен и Саусвелл применили метод релаксации Саусвелла для расчета вручную обтекания цилиндра вязкой несжимаемой жидкостью. В некоторых отношениях это была пионерская работа в численной гидродинамике. Для представления круговой границы на регулярной прямоугольной сетке использовалось конформное преобразование. Были получены численно устойчивые решения при числе Рейнольдса, равном 1000, что превышает физический предел устойчивости. При проведении вычислений авторы столкнулись с ясно выраженной тенденцией к неустойчивости при числе Рейнольдса, равном 100, и связали это с тенденцией к физической неустойчивости потока, предвосхитив тем самым современное понятие численного моделирования. Их работа может также считаться образцом финансирования научных исследований: на ее проведение Лондонскому имперскому колледжу в 1945 г. были выделены большие ассигнования фирмой по пошиву одежды!
Метод Саусвелла не так просто приспособить к использованию на ЭВМ. Вычислитель вручную просматривал матрицу в поисках максимальной невязки гораздо быстрее, чем производил арифметические операции. Для ЭВМ скорость просмотра матрицы не намного превышает скорость выполнения арифметических операций, и поэтому здесь становится более эффективным проведение релаксации последовательно во всех узлах сетки до сведения невязки к нулю, что идентично методу Либмана.
Таким образом, применение ЭВМ дало основание к дальнейшему развитию методов типа метода Либмана с использованием преимуществ идеи верхней релаксации Саусвелла. В 1950 г. Франкел (и в 1954 г. независимо от него Янг) разработал метод, который он назвал экстраполированным методом Либмана и который впоследствии стал называться методом последовательной верхней релаксации (Янг [1954]) или методом оптимальной верхней релаксации. Франкел подметил также аналогию между итеративным решением эллиптических уравнений и решением шагами по времени параболических уравнений, что имело важные последствия.
С развитием ЭВМ стали по-настоящему уделять внимание и уравнениям параболического типа, поскольку стало возможным рассчитывать нестационарные решения. В первой монографии Рихтмайера [1957], внесшей большой вклад в развитие, одномерной нестационарной гидродинамики, было приведено свыше десяти численных схем. В многомерном случае первым неявным методом был метод Кранка-Николсона, опубликованный в 1947 г. и требовавший итераций на каждом временном слое. Этот метод остается одним из самых популярных и лежит в основе широко используемого метода расчета неавтомодельных решений уравнений пограничного слоя (Блоттнер [1970]).
Невозможно точно определить, когда впервые была выдвинута идея асимптотического метода установления по времени, при которой для получения стационарного решения интегрируются уравнения нестационарного течения. Сомнительно, чтобы такая идея могла серьезно рассматриваться до появления ЭВМ.
Многие из пионерских работ в области вычислительной гидродинамики были выполнены в Лос-Аламосской лаборатории. Именно в Лос-Аламосе во время второй мировой войны фон Нейман разработал свой критерий устойчивости параболических конечно-разностных уравнений и дал метод исследования линеаризованной системы. Краткий отчет о его работах появился в открытой литературе лишь в 1950 г. (Чарни с соавторами; [1950]). В этой важной статье были впервые приведены расчеты метеорологических задач большого масштаба, в которых рассматривались нелинейные уравнения для вихря. Авторы выяснили, что в смысле устойчивости уравнения для вихря имеют преимущество над традиционными уравнениями для простейших физических переменных (скорость и давление), и привели эвристические обоснования своей трактовки нестационарной задачи как задачи с математически неполными условиями на входной и выходной границах.
Несмотря на то, что Нейманом в 1950 г. указывалось на необходимость введения понятия устойчивости конечно-разностных схем, четкие формулировки понятия устойчивости как аналога корректности для дифференциальных уравнений и теоремы о том, что из устойчивости и аппроксимации следует сходимость, впервые даны советскими математиками В. С.Рябеньким [1952] и А.Ф.Филипповым [1955].
В середине пятидесятых годов в работах Писмена и Ракфорда [1955], а также Дугласа и Ракфорда [1956] были предложены эффективные неявные методы для решения параболических уравнений, пригодные при произвольно больших шагах по времени. Под названием неявных схем метода чередующихся направлений они применялись и для решения эллиптических задач с использованием аналогии Франкела [1950] между продвижением решения по времени в параболических задачах и продвижением решения по итерациям в эллиптических задачах.
Неявные схемы чередующихся направлений, вероятно, наиболее популярны при расчетах задач о течениях несжимаемой жидкости, в которых используется уравнение переноса вихря.
В 1953 г. Дюфорт и Франкел опубликовали свою схему «чехарда» для параболических уравнений, которая, как и неявные схемы метода чередующихся направлений, пригодна для произвольно больших шагов по времени (при отсутствии конвективных членов), но сохраняет все преимущества чисто явных схем. Эта схема использована Харлоу и Фроммом [1963] при получении их широко известного численного решения, для нестационарной вихревой дорожки.
Требование единообразия и однородности вычислительного алгоритма для класса задач привело к понятию однородных разностных схем, которое было сформулировано А.Н.Тихоновым и А.А.Самарским [1950-1964]). С.К.Годунов [1958] использовал законы сохранения при выводе схем сквозного счета разрывных решений газовой динамики. Большое значение для решения задач гидродинамики и теплофизики имеет метод интегральных соотношений, предложенный А.А.Дородницыным и развитый О.М.Белоцерковским, П.И.Чушкиным и др., в котором использована частичная разностная аппроксимация уравнений, записанных в дивергентной форме на основе метода прямых. Эти методы сыграли существенную роль в формировании общего взгляда на конструкцию разностных схем для квазилинейных уравнений.
Статья Харлоу и Фромма [1965], опубликованная на страницах журнала Scientific American, была предназначена специально для того, чтобы привлечь внимание научной общественности Соединенных Штатов к возможностям вычислительной гидродинамики. Примерно в то же время во французском журнале La Houille Blanche появилась аналогичная статья Макано [1965]. В обеих этих статьях были впервые четко сформулированы понятия численного моделирования и численного эксперимента. Выходом этих статей можно датировать возникновение вычислительной теплофизики как отдельной дисциплины.
Вычислительная устойчивость всех упомянутых выше зависящих от времени решений была ограничена сверху по числу Рейнольдса (принципиально этот предел определяется сеточным числом Рейнольдса, т. е. числом, полученным по размеру шага ячейки конечно-разностной сетки). В 1966 г. Томан и Шевчик добились, по-видимому, неограниченной вычислительной устойчивости, используя для представления конвективных членов разности против потока и уделяя особое внимание граничным условиям. Их расчеты обтекания цилиндра простирались до чисел Рейнольдса, равных миллиону; они даже могли «вращать» цилиндр и получать магнусову подъемную силу, не сталкиваясь при этом с вычислительной неустойчивостью. Несмотря на то, что их схема имела лишь первый порядок точности, согласование полученных ими результатов с экспериментальными данными заставило переоценить важность формального порядка ошибок аппроксимации при разностном представлении дифференциальных уравнений в частных производных. В этой связи представляется важной работа Чена [1968], установившая существенное влияние численной постановки граничных условий.
Прямые (неитеративные) методы Фурье для численного решения эллиптического уравнения Пуассона были известны уже в течение некоторого времени (см., например, монографию Вазова и Форсайта [1960]), но не применялись к задачам гидродинамики. В 1965 году Хокни разработал родственный, но более быстродействующий метод, позволивший эффективно решать большие краевые задачи для уравнения Пуассона. После выхода этой работы прямые методы для уравнения Пуассона стали развиваться более интенсивно.
Описанные выше методы пригодны для расчета дозвуковых течений сжимаемой жидкости. Сверхзвуковые задачи отличаются от дозвуковых в нескольких важных аспектах, важнейшим из которых является возможность возникновения в сверхзвуковом течении ударных волн (т. е. разрывов в решениях).
Основополагающей работой для численного расчета гиперболических уравнений явилась статья Куранта, Фридрихса и Леви, опубликованная в 1928 г. Здесь обсуждались характеристические свойства уравнений и в общих чертах излагался известный метод характеристик. В этой работе было также получено и объяснено знаменитое необходимое условие устойчивости Куранта-Фридрихса-Леви, гласящее, что при расчетной сетке, не совпадающей с характеристической, область зависимости разностных уравнений должна по крайней мере включать в себя область зависимости дифференциальных уравнений. Это условие устойчивости КФЛ (которое в современной терминологии просто гласит, что число Куранта должно быть меньше единицы) справедливо для уравнений гидродинамики как в лагранжевых, так и в эйлеровых переменных.
Лагранжевы методы, в которых прослеживаются «частицы», были доведены до высокой степени совершенства в Лос-Аламосской лаборатории (Фромм [1961]). Вообще говоря, для двумерных задач эйлеровы методы предпочтительнее, однако при их использовании затрудняется расчет ударных волн. Если размер ячейки сетки не меньше, чем толщина ударной волны, то появляются осцилляции, снижающие точность. Эти осцилляции на дискретной сетке имеют физический смысл (Рихтмайер [1957]). Кинетическая энергия, высвобождающаяся из-за потери скорости при переходе через ударную волну, превращается во внутреннюю энергию случайных соударений молекул; при расчетах роль молекул играют ячейки конечно-разностной сетки.
Наиболее обычным подходом к расчету ударных волн на эйлеровой сетке является «размазывание» скачка на несколько ячеек сетки путем явного или неявного введения искусственной вязкости, не оказывающей влияния на решение на некотором расстоянии от ударных волн. В 1950 г. фон Нейман и Рихтмайер предложили схему искусственной вязкости, в которой «коэффициент вязкости» был пропорционален квадрату градиента скорости. Ладфорд, Полячек и Зегер [1953] просто брали большие значения физической вязкости в уравнениях течения вязкой жидкости на лагранжевой сетке, однако в их методе требовались нереально высокие значения вязкости.
Для размазывания скачка вместо явного введения искусственной вязкости можно использовать и неявную вязкость, имеющую место в конечно-разностных аппроксимациях. Это было осуществлено в широко известном методе частиц в ячейках (методе РIС), разработанном в Лос-Аламосе Эванс, Харлоу и др., а также в методе Лакса (Лакс [1954]) и в других методах.
В работе Лакса, опубликованной в 1954 г., сама численная схема гораздо менее важна, чем использованная форма дифференциальных уравнений — консервативная форма. Лакс показал, что преобразованием обычных уравнений гидродинамики, в которых зависимыми переменными являются скорость, плотность и температура, можно получить систему уравнений, в которой в качестве зависимых переменных служат количество движения, плотность и удельная внутренняя энергия торможения. Эта новая система уравнений отражает сущность физических законов сохранения и позволяет сохранять интегральные характеристики течения в конечно-разностной схеме. Такая система уравнений широко используется в настоящее время для расчета распространения ударных волн независимо от применяемых конечно-разностных схем, поскольку скорость плоской ударной волны точно рассчитывается любой устойчивой схемой (Лонгли [1960] и Гари [1964]).
Размазывание ударной волны при помощи неявной схемной вязкости осуществляется и в некоторых других методах. Так, в настоящее время широко применяется схема Лакса-Венд-роффа [1960] и ее двухшаговые варианты, например схема Рихтмайера (см. Рихтмайер [1963]). В методе РIС и в его модификации ЕIС (метод взрыва в ячейках), разработанной в 1964 г. Мадером, размазывание скачков достигается за счет введения конечного числа рассчитываемых частиц. Этот прием дает также возможность рассматривать поверхности раздела в жидкости (см. Харлоу и Уэлч [1965, 1966], а также Дали [1967]). В методе РIС, как и в более раннем методе Куранта-Изаксона-Риса [1952], используются односторонние разности для первых производных по пространству и таким образом вводится своего рода схемная вязкость, однако эти методы сохраняют истинные характеристические свойства дифференциальных уравнений. Хотя во всех этих методах неявно используются диссипативные члены, размазывающие ударные волны, для обеспечения устойчивости каждого из них в некоторых частных случаях требуется введение дополнительных членов с явной искусственной вязкостью.
Для уравнений и систем гиперболического типа в прямоугольных областях различные экономичные схемы построены и исследованы Е.Г.Дьяконовым [1963], А.Н.Коноваловым [1964], А.А.Самарским [1965].
Для численного решения системы уравнений реагирующих течений экономичные численные методы были построены и реализованы Л.А.Чудовым, Н.Н.Яненко, Б.Г.Кузнецовым, Д.Б.Сполдингом, Оран и др.