Численные методы теплофизики

Численные методы теплофизики

Введение

Предмет вычислительной теплофизики

В прошлом в теплофизике, как и в других физических науках, существовало два основных метода исследований: теоретические и экспериментальные. Зададимся вопросом: в каком отношении к этим методам находится вычислительные методы? Можно ответить, что они являются отдельными новыми методами исследования, хотя и обладают некоторыми чертами теоретических и экспериментальных методов и скорее дополняют, чем заменяют их.

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

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

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

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

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

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

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

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

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

Принципы построения разностных схем

Основные понятия и обозначения теории разностных схем

Существо метода конечных разностей состоит в следующем.

Пусть состояние среды описывается дифференциальным уравнением изменения функции f непрерывного аргумента. Область непрерывного изменения аргумента заменяется некоторым конечным множеством точек (узлов). Расстояния между узлами называют шагами. Совокупность узлов составляет сетку, которая может быть равномерной или неравномерной. Если функция f зависит от одной переменной, то сетка является одномерной, если f является функцией нескольких переменных, то сетка является многомерной. Конечно-разностная сетка может быть равномерной по одной переменной и неравномерной – по другой. Например, на рисунке 1 изображена сетка, равномерная по х и неравномерная по у.

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

Таким образом, дифференциальное уравнение для функции f непрерывного аргумента заменяется алгебраическим конечно-разностным уравнением для сеточной функции f_D. Конечно разностная схема представляет собой систему дискретных алгебраических уравнений, аппроксимирующих дифференциальное уравнение с соответствующими граничными условиями. В качестве приближенного решения дифференциального уравнения принимается решение соответствующего разностного уравнения – сеточная функция – в виде одно- или многомерной таблицы.

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

В дальнейшем мы будем использовать следующее модельное уравнение:

(2.1)

Это уравнение выбрано потому, что оно учитывает локальное изменение искомой функции во времени (первое слагаемое в левой части), конвективный перенос (второе слагаемое в левой части) и молекулярный перенос (член в правой части уравнения). Коэффициент а равен кинематической вязкости n, если f представляет собой скорость u, коэффициенту диффузии D, если f является концентрацией с, и, наконец, это будет коэффициент температуропроводности, если f есть температура Т. В последнем случае можно сказать, что при f º Т уравнение (2.1) описывает нестационарный процесс изменения температуры в одномерном канале при течении в нем жидкости с постоянной скоростью u.

Описание неустойчивости

Рассмотрим модельное уравнение (2.1):

Конечно-разностная схема с разностями вперед по времени и центральными разностями по пространственной переменной имеет вид:

Перепишем это уравнение следующим образом:

. (3.4)

Введем следующие обозначения:

— число Куранта,

— диффузионное число.

После этого уравнение (3.4) запишется так:

. (3.5)

Пусть на n-ном слое по времени возникло малое возмущение e_iⁿ в виде, представленном на рисунке 5, а.

Рисунок 5.

e_iⁿ возрастает от точки к точке и имеет разные знаки в соседних узлах. Такие возмущения могут порождаться либо ошибками округления, либо поперечными движениями в реальной двумерной задаче. Возмущение, возникшее на n-ном слое, непременно проявится на (n+1)-ом слое. Если возмущение на следующих временных слоях уменьшается по абсолютному значению, то конечно-разностная схема является устойчивой. Таким образом, условие устойчивости в данном случае имеет следующий вид:

или:

. (3.6)

Проследим за развитием возмущения. Для этого запишем уравнение (3.5) с учетом наложенного на n-ном слое возмущения:

(3.7)

Вычтем из уравнения с возмущениями (3.7) «невозмущенное» уравнение (3.5):

Найдем отсюда приращение возмущения на (n+1)-ом слое:

. (3.8)

Рассмотрим это уравнение только с одним диффузионным членом, т.е. будем считать, что С=0 :

. (3.9)

Проанализируем это уравнение в точке i:

т.к. eⁿ_i₊₁> 0, eⁿ_i_-1> 0, e_iⁿ< 0. Следовательно, De_i > 0 и стремится скорректировать отрицательное возмущение в точке i. Аналогично, рассматривая De_i₊₁, имеем:

т.к. eⁿ_i₊₂< 0, eⁿ_i< 0, eⁿ_i₊₁> 0, т.е. положительное возмущение eⁿ_i₊₁ корректируется отрицательным приращением De_i₊₁. Таким образом, можно изобразить, каким будет приращение возмущения на (n+1)-ом слое (пунктирная линия на рисунке 5, б). Складывая графически e_iⁿ и De_i, получим вид возмущения на n-ном слое (жирная линия на рисунке 5, б):

e_iⁿ⁺¹=e_iⁿ+De_i.

Из рисунка 5, б видно, что в данном случае условие устойчивости (3.6) выполняется, и, следовательно, конечно-разностное уравнение (3.9) устойчиво . Однако если шаг Dt слишком велик, то поправка за счет приращения De_i₊₁ окажется чрезмерной. Для таких слишком больших Dt величина нового e_iⁿ⁺¹ будет больше начального возмущения e_iⁿ, как это показано на рисунке 5, в.

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

Динамическую неустойчивость можно устранить, наложив ограничения на шаг по времени Dt.

Рассмотрим теперь уравнение (3.8) с одним конвективным членом, то есть предположив, что d=0: Уравнение (3.8) в этом случае примет вид:

. (3.10)

Пусть u> 0, то есть С> 0. Приращение возмущения в точке i:

т.к. eⁿ_i₊₁> 0, eⁿ_i_-1> 0, но амплитуда e растет с ростом i, то есть eⁿ_i₊₁> eⁿ_i_-1, следовательно выражение в скобках положительно. Таким образом, приращение De_i, обусловленное конвекцией, усиливает возмущение e_iⁿ. Это означает, что ошибка растет монотонно (рисунок 5, г).

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

Таким образом, уравнение (3.5) при С=0 является условно устойчивым, а при d=0 — абсолютно неустойчивым. Если же одновременно С≠ 0 и d≠ 0, то конвективный и диффузионный члены будут взаимодействовать друг с другом, и в целом уравнение (3.5) может оказаться как устойчивым, так и неустойчивым. Чтобы ответить на вопрос об устойчивости уравнения (3.5), необходимо проанализировать его на устойчивость одним из методов, которые будут рассмотрены ниже.

Метод дискретных возмущений

Идея метода дискретных возмущений состоит в том, что в уравнение в каждую точку поочередно вводится дискретное возмущение e и прослеживается влияние этого возмущения на следующих временных слоях. Конечно-разностная схема будет устойчивой, если возмущение не возрастает, то есть выполняется условие (3.6).

Рассмотрим конечно-разностное уравнение (3.5). Введем возмущение e в точку (i, n):

Вычтя из этого уравнения с возмущением «невозмущенное» уравнение (3.5), получим уравнение для возмущений:

Отсюда:

Необходимо разрешить полученное неравенство относительно d:

-1≤ 1-2d≤ 1.

Рассмотрим отдельно правую и левую части этого неравенства:

а) 1-2d≤ 1 d≥ 0 — это условие выполняется всегда.

б) -1≤ 1-2d d≤ 1 — конечно-разностная схема (3.5) устойчива при выполнении этого условия.

Необходимо еще потребовать, чтобы знаки eⁿ⁺¹ и eⁿ были одинаковыми во избежание осцилляций, то есть должно выполняться условие:

В данном случае это условие означает выполнение следующего неравенства:

1-2d≥ 0 Þ d≤ 1/2.

Мы получили более жесткое ограничение, чем в б), которое включает в себя условие . При фиксированных Dх и a это условие накладывает ограничение на шаг по времени:

Это ограничение является жестким в смысле затрат времени для расчета на компьютере. Например, предположим, что расчет ведется с некоторым пространственным шагом Dх₁. Максимально возможный шаг по времени . Если надо провести расчет с вдвое меньшим шагом Dх₂=Dх₁/2, то шаг по времени , т.е. шаг по времени придется уменьшить в 4 раза, а затраты компьютерного времени в общей сложности увеличатся в 8 раз. В двумерной задаче уменьшение вдвое шагов Dх и Dy увеличивает затраты времени в 16 раз, а в трехмерной задаче - в 32 раза!

Введем теперь возмущение e в точку (i+1, n):

Вычтя из этого уравнения с возмущением «невозмущенное» уравнение (3.5), получим новое уравнение для возмущений:

Отсюда условие устойчивости имеет вид:

а) ,

(3.11)

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

Последнее неравенство можно переписать следующим образом:

Выражение слева представляет собой сеточное число Пекле: . Тогда последнее неравенство можно переписать так:

Pe_c< 2.

Это условие можно воспринимать как ограничение на шаг по x:

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

б)

или

Поскольку , то это неравенство выполняется всегда.

Проведя аналогичный анализ для точки (i-1, n), получим:

Отсюда условие устойчивости:

Левая часть этого неравенства выполняется всегда, и никаких ограничений на шаги сетки не накладывает. Правая часть неравенства дает следующее условие:

Это условие более ограничительное, чем (3.11), поэтому условие (3.11) можно не принимать во внимание.

Таким образом, метод дискретных возмущений для уравнения (3.5) дает три следующих условия:

а)

б)

в)

Из условия б: uDх< 2a, из условия a: . Следовательно, получаем еще одно условие:

г) .

Метод фон Неймана

Рассмотрим этот метод на примере того же конечно-разностного уравнения (3.5). Сеточная функция f_iⁿ представляется в виде разложения Фурье:

f_iⁿ = V_nе^Ii^q

где V_n – амплитуда отдельной компоненты с волновым числом на n-ном временном слое, q=kDx – фазовый угол, – мнимая единица. Аналогично

= V_n+1 е^-I(i±1)^q

Подставим эти выражения в уравнение (3.5):

Разделив на е^Ii^q, получим:

Учтем, что е^I^q+ е^-^I^q=2cosq, е^I^q- е^-^I^q=2Isinq:

V_n+1=V_n [1-2d(1-cosq)-ICsinq]

Определим множитель перехода G следующим образом:

V_n+1=GV_n

Для того, чтобы решение оставалось устойчивым, необходимо потребовать, чтобы

|G|≤ 1. (3.12)

В данном случае множитель перехода равен:

G=1-2d(1-cosq)-ICsinq

Он является комплексным выражением, поэтому условие устойчивости |G|≤ 1 приводится к неравенству:

[1-2d(1-cosq)]²-C²sin²q ≤ 1

Сделав соответствующие преобразования, получим:

C²(1+cosq) ≤ 4d[1-d(1-cosq)]

Рассмотрим два предельных случая.

а) cosq=-1: 0 ≤ 4d(1-2d).

Отсюда получаем первое условие устойчивости:

Такое же условие было получено методом дискретных возмущений.

б) cosq=1: .

Учитывая предыдущее неравенство, получим второе условие:

C ≤ 1.

Комбинируя эти два условия, можно получить условие устойчивости в следующем виде:

или Pe_c ≤ 2.

Все эти полученные условия справедливы только в случае линейного уравнения при u=const.

3.3.3 Метод практической устойчивости

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

Перепишем уравнение (3.5) в виде:

.

Условия устойчивости запишутся следующим образом:

;

;

– выполняется всегда;

– выполняется всегда.

Мы получили тремя разными методами одинаковые критерии устойчивости. Такой результат получается только для линейных дифференциальных уравнений с постоянными коэффициентами.

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

В 1952 г. В.С. Рябеньким была сформулирована теорема эквивалентности, которая устанавливает эквивалентность устойчивости и сходимости при выполнении следующих условий:

- решение дифференциального уравнения в частных производных должно непрерывным образом зависеть от начальных условий;

- конечно-разностное уравнение должно аппроксимировать дифференциальное уравнение в частных производных;

- устойчивость должна быть определена в форме фон Неймана.

При выполнении этих требований необходимое условие устойчивости становится и достаточным для сходимости.

Для нелинейных уравнений эта теорема не доказана, но на практике ее используют применительно и к нелинейным уравнениям.

Вообще говоря, исследование строгих определений аппроксимации, устойчивости, сходимости при Dх→ 0 и Dt→ 0 занятие зачастую бесплодное, т.к. реальные расчеты проводятся при конечных Dх и Dt.

Численные методы теплофизики

Введение

12 3 4 5 6 Следующая ⇒