Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена.djvu), страница 7

Неявную схему Эйлера можно получить, если использовать разложение в ряд Тейлора в точке т~+,. Т! ж Т[+' — Т' (т»+,) Лт. Тогда придем к схеме и!+' =из+Лт)(т~э» и/+'), (!.34) прн реализации которой на каждом шаге необходимо решать в общем случае нелинейное уравнение относительно иг+'. Проведем сопоставление явной н неявной схем Эйлера. С точки зрения объема вычислений для одного шага явная схема имеет преимущество.

Только в случае, когда функция 1 (т, Т) линейна относительно Т, т. е. 1 (т) = а (т)Т + 6 (т). вычисления по неявной схеме не сложней, чем по явной, поскольку тогда уравнение (1.34) разрешается относительно и~+'. (1.35) и!+ ' = [из+ Ь (тзэ,) Лт)! [! — а (тэ,) Лт!. Если же ! (т, Т) — нелинейная функция Т, то приходится решать на каждом шаге уравнение (!.34) каким-либо итерационным методом. Рассмотрим теперь вопрос о погрешностях численных решений, получаемых по явной и неявной схемам Эйлера. Для этого введем понятия аппроксимации н устойчивости. Вернемся к разностному уравнению (!.33), переписав его в виде =)(т>, и!), (1.36) Ьт Поскольку это уравнение получается из соотношения (1.32) для Т~ путем отбрасывания малых членов, точное сеточное решение Т! не удовлетворяет уравнению (!.36).

Величину ф1, характеризующую чстепень неудовлетворения» функции Т! разностному аналогу (1.36) дифференциального уравнения (1.29), очевидно, можно определить как рассогласование (невязку) между левой и правой частями (1.36), которое возникает при замене и! на Т'.

т/+ — т/ о(я = — /(тп Т!). Лт Эта величина ф/ называется погрешностью аппроксимации ис- ходного дифференциального уравнения разностным уравнением. Из разложения в ряд Тейлора (1.32) нетрудно найти, что /Лот) дт ( а'т ~ Лт т. е. при достаточно малых Лт погрешность аппроксимации убывает пропорционально первой степени Ьт: )ф~) ВЛт, В= ! фг=о(Л ). (!.3Т) В этом случае говорят, что разностная схема (1.33) имеет первый порядок аппроксимации.

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

е. при увеличении номера временного слоя /. Для анализа этого вопроса исследуем поведение погрешности численного решения простейшего линейного уравнения теплового баланса йТ вЂ” + по, Т = О, Тм о = То, дт где то — темп охлаждения 151, а Т имеет смысл перегрева над температурой окружающей среды. Точное решение задачи имеет вид Т(т) = То г — "'. (1.39) Разностное решение, найденное по явной схеме Эйлера, можно записать в виде иг+' =иг(1 — тоЬт) =и'(1 — гложат)', а по неявной схеме Эйлера — в виде и/+' = их/(1+тойт) =ио/(1+гло Ма)/.

(1.41) Точное решение ограничено: Т (т) «„То и Т (т- оо) = О. Для явной схемы, как следует из (1.4О), разностное решение будет ограниченным только при ! 1 — гпо Лт ) < 1 или Лт < 2/поо. При Лт ) 2/т«разностное решение ие будет представлять собой осциллирующую функцию, амплитуда колебаний которой неограниченно возрастает с увеличением) при движении вдоль оси т. Таким обРазом, погРешность РешениЯ по Явной схеме также неогРаниченно возрастает.

Такое явление называется неустойчивостью розностной схемы. Разностное решение ис, найденное по неявной схеме (1.41), монотонно убывает с увеличением 1, и погрешность ее=- Т~' — иг остается ограниченной при любом шаге Лт. Условие устойчивости разностной схемы для решения задачи (1.29), (!.30) может быть сформулировано в одном из вариантов как требование ограниченности погрешности решения ее для любых тз при / — » оо: [ег~ - С, С= сопз1.

(1.43) Более подробно с проблемой устойчивости разностных схем для решения обыкновенных дифференциальных уравнений можно ознакомиться по книге 129). Таким образом, неявная схема Эйлера устойчива при любых значениях Лт, нли безусловно устойчива. Явная схема устойчива лишь при выполнении ограничения на значение шага (!.42), или условно устойчива. При попытках проводить расчеты с шагами Лт, превышающими предельно допустимые из условия устойчивости значения, происходит «раскачка» («разболтка») разностного решения, приводящая к абсурдным числовым результатам или даже к машинному оста- нову из-за переполнения разрядной сетки.

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

Если время переходного процесса всей системы определяется телами с большой тепловой инерцией, то может появиться необходимость проводить расчет с шагом Лт, который превышает постоянные времени тел с малой тепловой инерцией. Действительно, если выбрать шаг из условия Л«( 2!т,„, где т,„— максимальный из темпов охлаждения отдельных тел, то может потребоваться чрезвычайно большое число шагов для расчета всего нестационарного процесса. Рассмотренные нами схемы Эйлера имеют первый порядок аппроксимации. Для построения схем с более высоким порядком в разложении (1.32) нужно оставить члены более высокого порядка малос- 31 ти по Лт.

Например, удерживая члены второго порядка и учитывая, что ФТ д I 6Т'! д д(, д( дТ вЂ” = — ~ — 1=- — (7(т, Т)1= — + — —, дт1 дт ~ дт / дт дт дТ дт получаем иг и ' ==- иг -!- ЛтГ' (тго иг) + — ' [ ~ — ') + ~ — ) 1(тп и') 1. (1. 44) лР Г д(~г lд( У 2 (Г дт (ГдТ ) Использование схем, подобных (1 44), затруднено тем, что требуется вычислять производные от функции ) (т, Т), и в настоящее время они применяются не часто. Поэтому возникает необходимость в построении схем с высоким порядком точности, которые не содержали бы производных от Г (т, Т). К таким схемам относятся схемы Рунге — Кутта.

Схемы Рунге — Кутта. Лля получения схем Рунге — Кутта запишем приращение сеточной функции точного решения на промежутке !тп гз.ы1 в виде Тг~,— Т;= ~ )(т, Т)бт=!д;„Г, (1. 45) Схемы Эйлера можно трактовать как полученные путем замены интеграла в правой части (1.45) простейшими формулами 1дзэ, ж Лт~(тго иг) или 7пг+, Лт~(тзГ„иг+Г). Порядок точности разностной схемы можно повысить, если использовать более сложные квадратурные формулы, в которых производная 7 (т, Т) вычисляется не в одной, а в нескольких точках отрезка !тгп тз+,1.

При таком подходе возникает задача определения приближений Г и соответственно решения и в этих промежуточных точках. Эти приближения вычисляются последовательно по мере про- ДВИЖЕНИЯ ПО ОТРЕЗКУ (тго т;+,! От ТОЧКИ т; К ТОЧКЕ т;+,. Тая КаК ПРИ этом согласно (1.29) функция г(т, Т) равна производной от решения Т (т), то приближения для решения и строятся на основе оценок значений его производной — значений функции ) (т, и).

Поэтому в окончательных формулах приближение для 1 (т, Т) в определенной точке выражается через приближения Г'(т, Т) в предыдущих точках, см. ниже формулы (1.47), (1.48). Рассмотрим для примера простейшую схему, в которой для вычисления интеграла (1.45) используют две точки: тз и т;+,. Для ее реализации необходимо сначала построить первое приближение значения Т (т) в точке тз+, — и. Его находят, используя значения 7 (тп иГ), по явной формуле Эйлера, и=ит+Лтг(т;, иг). (1.46) Подчеркнем, что это первое приближение носит промежуточный характер.

Далее для определения окончательного значения игч х вычисляют интеграл (1.45) по формуле трапеций с учетом (1.46): и!+' = пг+ Лт ~ — ~(тл и() + — ~(т;+„и)~. (1.47) [2 2 Нетрудно доказать, что погрешность аппроксимации построенной схемы (1.46), (1.47) равна [ф[ = О (Лт'). Рассмотренная схема является простейшей двухэтапной схемой Рунге — Кутта. В ней интеграл определяется по двум точкам интервала [т;, тьы) и использУютсЯ два вычислениЯ фУнкции 7 (т, и) на одном шаге по времени. В общем случае при использовании для определения интеграла 7~ ~+, т точек на интервале [тл тг ы) получается льэтапная схема Рунге — Кутта. Первая точка для вычисления оценки производной ) (т, и) всегда совпадает с тв а остальные располагаются оптимальным образом сточки зрения обеспечения наивысшего прн данном т порядка аппроксимации.