Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 17
Текст из файла (страница 17)
Действительно, из условия аппроксимации а из условия устойчивости ()е'„(! ( В))ф„'(), но тогда 1'т ))е„'()=11ш()Т„' и„''~1=(), т. е. разностное решение сходится к точному, причем порядок точ- ности схемы совпадает с порядком аппроксимации. $ ЗЛ. ЯВНАЯ И НЕЯВНАЯ СХЕМЫ Для решения нестационарных одномерных задач теплопроводности можно построить большое число «разумных» разностных схем.
Однако мы рассмотрим только те разностные схемы, важность которых подтверждена вычислительной практикой. В связи с наличием в нестацнояарном уравнении теплопроводности двух дифференциальных операторов — по временной и пространственной переменным — различают два вида схем: явные и неявные. Рассмотрим особенности этих схем на примере решения одномерной нестационарной задачи (3.1) — (3.3) на равномерных пространственной и временной сетках (см. рис. 3.1).
' Разностную схему для определения разностиого решения будем по-прежнему строить, заменяя в уравнении (3.1) и граничных условиях (3.2), (3.3) производные конечными разностями. Рассмотрим аппроксимацию производной по времени. В принципе для построения соотношений, аппроксимирующих временную производную, в 1-й момент времени можно использовать значения температур в различные моменты времени: Т~', Тт-,', Тт-„», .... Однако на практике в подавляющем большинстве случаев используются только значения температуры в 1-й и (1 — 1)-й моменты времени. Такие схемы называются двухслойными (повремени).
Значительно реже учитывают значение температуры в (1 — 2)-й момент времени и получают трехслойные схемы. Дальше мы будем рассматривать только двухслойные схемы. В этом случае производную по времени аппроксимируют разностью назад вт тг — тг — ' в лт При двухслойной схеме пространственный дифференциальный оператор целесообразно аппроксимировать также на основе значений температуры в )сй и (1 — 1)-й моменты времени. При этом наибольшее распространение получили два «крайиих» случая.
В первом случае при аппроксимации используются только значения температуры для искомого, «настоящего» /-го момента времени: д»У — ж (Т'„+ ~ — 2Т'+ Т'„.,)(й» и« вЂ” иу' л « Лт — (и'„'1 — 2и,', '+и„' 1)+д,. (3.23) а« Для аппроксимации граничных условий (3.2) пока воспользуемся простейшим способом замены производных по координате правой и тевой разностями.
Соответственно получаем; и« вЂ” и, «« — " а«"' =Ч« Ь (3.24) и,„— ид ! ! х та,и~=дь а (3. 25) Начальное условие для разностной схемы задается точным образом: и„' == Т„(х„), и = 1„, Ж. (3.26) В 4 3.1 мы уже доказали, что записанные разностные уравнения аппроксимируют уравнение теплопроводности с порядком О (Лт+ -й«) и граничные условия с порядком О (Ь). Уравнения (3.22) или (3.23) вместе с уравнениями (3.24)— — (3.26) образуют разностные схемы, позволяющие найти сеточную функцию и~. Рассмотрим, в чем заключается принципиальная разница между схемами, использующими уравнения (3.22) и (3.23).
Начнем со схемы (3.23) — (3.26). Уравнение (3.23) позволяет выразить в явном виде неизвестное значение иг' сеточной функции на «новом» временном слое / через известные значения сеточной функции на предыдущем (! — 1)-м слое: и!, ай« (и,', 1' — 2и'„' —,-и',...'1)!Ч+ д, А«7ср+ и'„'.
(3.27) Так как в начальный момент времени значения и'„заданы условием (3.26), то по формуле (3.27) можно найти сначала и„' для внутренних узлов и --= 2, ..., А! — 1, а затем из граничных уравнений !3.24), (3.25) определить и, 'и ич. Аналогичная процедура проводится для отыскания сеточной функции на втором временном слое и т. д.
80 ! способ аппроксимации второй производной д«Т)дх«был рассмотрен в. 6 3.1), а во втором — только значения температуры для «прошлого», предыдущего момента времени (! — 1): аг ж (Т„'~ — 2Т'„' + Т„''1 ) )й«. Соответственно получают два различных разностных уравнения, аппроксимирующих уравнение теплопроводности (3.1): г р " " — — '(и'„.г ~ — 2и'„+ и'„~) + д„ (3.22) Л« а« Разностная схема (3.23) — (3.26) называется явной, так как позволяет искомые значения сеточной функции и( в явном виде выразить через найденные ранее значения ие — '.
Алгоритм численного расчета по явной схеме очень прост и легко программируется. Перейдем к разностной схеме, заданной уравнениями (3.22), (3.24) — (3.26). Здесь ситуация сложнее, поскольку в каждое уравнение вида (3.22) кроме неизвестного значения и,', для и-й пространственной точки входят еще два искомых значения сеточной функции и,(+, н и„~', для соседних (и — !)-й и (п + !)-й точек. Поэтому рассмотренный выше для явной схемы прием получения явной формулы для неизвестного значения и„' в этой ситуации не проходит. Все искомые значения (иг„)„н, оказываются «завязанными» друг сдругом в общую нераспадающуюся систему уравнений.
Эта система состоит из (!У вЂ” 2) уравнений (3.22) для внутренних узлов и двух уравнений (3.24), (3.25), соответствующих граничным условиям. Всего имеем А( уравнений относительно А! неизвестных (иД и,. Таким образом, в данном случае на каждом временном слое значения сеточной функции и„' определяются не по явным формулам, а из решения системы Ф уравнений, поэтому рассмотренная разностная схема называется неявной. Эффективный алгоритм решения системы уравнений (3.22), (3.24), (3.25) рассмотрим ниже. На первый взгляд явная схема предпочтительнее, так как она имеет такой же порядок аппроксимации 0 (Лт + Ь), как и неявная, но не требует решения на каждом шаге по времени систем У уравнений.
Однако более подробный анализ показывает, что явная схема условно устойчивая,т. е. устойчивая при определенном ограничении на величину шага по времени Лт. Условие устойчивости для явной схемы (3.23) — (3.25) имеет вид Лт ~ Лт„„=й»/2а. (3.28) Из условия устойчивости следует, что измельчение пространственной сетки должно сопровождаться измельчением временной сетки. Например, при увеличении числа пространственных узлов Ж в 4 раза, требуется увеличить число шагов по времени в ! 6 раз.
Необходимость соблюдения условия (328) приводит к тому, что при определении шага по времени для решения реальной нестационарной задачи мы не можем исходить только из характера протекания во времени изучаемого физического процесса. Это в ряде случаев приводит к неприемлемым затратам машинного времени. Кроме того, при неоправданно большом числе временных шагов может начать проявляться погрешность округления, возникающая в ЭВМ при реализации арифметических операций. Неявная схема (3.22) безусловно устойчивая, т. е. явление неустойчивости не возникает при любых величинах Лт. Поэтому при решении задачи по неявной схеме величину шага по времени задают 8! только из соображений обеспечения требуемой погрешности численного 7 решения.
,)=1 Проиллюстрируем «физический механизм» возникновения неустойчиво- О 7 ! ! 7 сти при рас е по явной схеме на х-7 -! д ие! «+г примере плоской стенки без источнид ков теплоты. Положим, что начальная температура стенки равна нулю -3 во всех точках пространственной сетки, кроме одной точки с номером Рис. 3.4 п=й(рнс. 3.4): иа=О, л=1, .„, Л1, и „-й й; и$ ="1. Пусть й = 1, а = 1. Выберем Лт = 1, т. е. величина шага вдвое превышает условие устойчивости (3.28). Расчет по явной схеме ведется по формуле (3.2!), которая вданном случае записывается в виде: пгв 3 ! ( / — ! 1 1 ' — 1 1 — 1 7 — 1 ип = — (и„+ ! — 2и'„+ пи !)+ им ! Выполняя расчет, получим следующие данные: Номер точки по проетраиетаеиио«моор«ипате Номер шага по иргмеии )=о 1=1 1=2 ! — ! 3 0 1 — 2 Из рис.
3.4 видно, что возникает «разболтка» и получается разностное решение с возрастающими по величине от слоя к слою значениями !иД, не имеющее ничего общего с точным решением, представляющим функцию, которая принимает только положительные значения и убывает в точке с номером й. «Разболтка» началась из-за того, что температура и» «упала» в отрицательную область и стала меньше, чем температуры соседних точек, что противоречит физическому смыслу. Нетрудно убедиться, что при расчете на границе устойчивости с шагом Лт = О,5 получим иат ! = 0,5, и»! = О, и» ! =- О,5, т.
е. все температуры положительны и не превышают температур соответствующих им соседних точек на предыдущем временном слое. Дальнейший расчет с шагом Лт = 0,5 привел бы к получению колеблющегося решения иг' с убывающей нормой ~!иф~. Расчет же по неявной схеме при любом Лт дает решение„правильно отражающее качественный характер изменения температуры. Рассмотренному отличию в поведении решений, полученных по явной н неявной схемам, можно дать следующее физическое обьяс- нение. ~.ри расчете по явной н неявной схемам предполагается, что функция меняется линейно на интервале [т! „ т![, но значение производной по времени при явной схеме вычисляется по значениям искомой функции в начале временного интервала, поэтому прираще ние искомой функции (и! — и!„— ') не зависит от получаемых зна чений, а абсолютная величина этого приращения пропорциональна шагу, см. (3.27).
В результате при некотором критическом шаге Лт мы можем получить новые значения иг, противоречащие физическому смыслу задачи, как это и было в рассмотренном примере. В неявной же схеме приращение (и! — и! — „') зависит от всех значений и1' на новом временном слое, т. е. имеется как бы кобратная связь», не позволяющая получать абсурдные приращения сеточной функции.