маслов. Практикум предназначен для проведения лабораторных занятий и са мостоятельной работы студентов факультета Ф, обучающихся по спе циальности Ядерные реакторы и энергетические установки
 Скачать 1.6 Mb.
|
|
1.1.4. Выбор шага интегрирования и оценка погрешности численного решения Обычно при реализации численных схем решения ОДУ преду- сматривают автоматический выбор величины шага интегрирова- ния. Этот выбор основывается на оценке локальной погрешности численного решения на одном шаге интегрирования, т.е. оценка погрешности величины u i+1 в точке t i+1 выполняется в предположе- нии, что в начале интегрирования в t i значение искомой функции известно точно. Одношаговые разностные схемы. При использовании одношаго- вых разностных схем для вычисления локальной погрешности ε i+1 на отрезке [t i , t i+1 ] искомая величина u i+1 определяется по одной и той же схеме с шагом ∆t 1 = (t i+1 - t i ) и шагом ∆t 2 = ∆t 1 /2. Оценка для локальной погрешности при меньшем шаге может быть получена по формуле ) 1 2 /( 1 2 2 1 1 1 − − = ε ∆ + ∆ + ∆ + p t i t i t i u u . (1.24) 19 Здесь 2 1 t i u ∆ + , 1 1 t i u ∆ + – значения искомой функции, рассчитанные с ша- гом ∆t 2 и ∆t 1 , соответственно; р – порядок аппроксимации разност- ной схемы. Оценка вида (1.24) справедлива при использовании не только одношаговых, но любых других разностных схем. При автоматическом выборе шага задают максимальное и ми- нимальное допустимые значения погрешности ε max и ε min (ε max > ε min ). Если оценка локальной погрешности ε i+ 1 лежит в пре- делах заданных допустимых значений, то шаг интегрирования ∆t на следующем шаге не меняется. Если ε i+ 1 > ε max , то расчет повто- ряется от точки t i с половинным шагом. В случае ε i+ 1 < ε min сле- дующий шаг выполняется с удвоенным шагом. Многошаговые разностные схемы. В случае использования ме- тода прогноза-коррекции для неявных многошаговых разностных схем оценка локальной погрешности может быть получена с по- мощью вычисляемых при реализации схемы значений искомой функции на этапах прогноза и коррекции и не требует дополни- тельных расчетов. Это справедливо при использовании схем одно- го порядка аппроксимации на обоих этапах. Действительно, если порядок аппроксимации схем равен р, то можно выразить точное значение искомой функции T i+ 1 = T(t i+ 1 ) на этапах прогноза и коррекции по формулам (1.25) и (1.26) соответ- ственно: , ) ( 1 1 1 0 1 1 + + + + + ξ ∆ + = p p p i i dt T d t A u T (1.25) ) ( 1 1 1 1 1 + + + + + η ∆ + = p p p s i i dt T d t B u T (1.26) Здесь А и В – постоянные коэффициенты, известные из сопостав- ления разностной схемы и разложения в ряд Тейлора; ξ, η – неиз- вестные промежуточные точки отрезка [t i , t i+ 1 ], в которых выполня- ется расчет р+1 производной. Считая, что р+1 производная при- 20 мерно постоянна на рассматриваемом отрезке интегрирования, из системы (1.25) и (1.26) можно получить ). ( 1 0 1 1 1 s i i s i i u u A B B u T + + + + − − + ≈ (1.27) Второй член в правой части выражения (1.27) является оценкой локальной погрешности для метода прогноза-коррекции. Оценка полной погрешности. Локальная погрешность не учиты- вает накопления погрешности в ходе всего расчета. Фактическая погрешность решения неизвестна, но при устойчивости используе- мой разностной схемы является ограниченной и соизмеримой с суммой локальных погрешностей на отдельных шагах. Оценку полной погрешности искомой функции в точке t i+ 1 мож- но получить из сравнения двух решений 1 + ′ i u и 1 + ′′ i u , полученных во всей области [0, t i+ 1 ] с постоянными шагами t′ ∆ и t ′′ ∆ ( t′ ∆ > t ′′ ∆ ), соответственно ( ) ( ) 1 1 1 1 1 − ′′ ∆ ′ ∆ ′ − ′′ ≈ ′′ − = ε + + + + p i i i i t t u u u T . (1.28) 1.2. Описание программы ODY_lab численного решения задачи Коши для ОДУ Учебная программа ODY_lab предназначена для изучения раз- ностных схем, используемых при численном решении задачи Коши для ОДУ. Она позволяет пользователю: − сформулировать исходную задачу Коши, выбрав рассматри- ваемую функцию f(t, T), задав начальное значение искомой функции Т и конечное рассматриваемое значение независимой переменной t; 21 − выбрать метод решения, указав используемую разностную схе- му и определив используемый шаг изменения независимой пе- ременной; − проанализировать полученные результаты, представляемые в табличном и графическом виде. 1.2.1. Особенности программы В программе реализована возможность задания двух видов функции f(t, T): 1. AT T t f ≡ ) , ( , где А – изменяемый пользователем коэффициент. Функция данного вида часто встречается в задачах охлаждения и нагрева. 2. ) /( ) ( ) , ( t C BT A T t f − − ≡ . Значения коэффициентов A, B, C за- даются пользователем. Функция данного вида обеспечивает не- линейность рассматриваемой задачи и описывает, например процесс реактивного движения. Выбранные функции с одной стороны позволяют легко проана- лизировать особенности поведения рассматриваемых разностных схем, а с другой имеют простое аналитическое решение, необходи- мое для определения точных погрешностей численного расчета. Для численного решения задачи в программе можно выбрать один из восьми методов численного решения: 1) явная схема Эйлера; 2) модифицированная схема Эйлера; 3) исправленная схема Эйлера; 4) схема Рунге-Кутта 4-го порядка; 5) явная схема Адамса 2-го порядка; 6) явная схема Адамса 4-го порядка; 7) метод прогноза-коррекции со схемами Адамса 2-го порядка; 8) метод прогноза-коррекции со схемами Адамса 4-го порядка. Шаг изменения независимой переменной в программе явно за- дается пользователем и не меняется в ходе расчета. Это позволяет исследовать его влияние на сходимость рассматриваемых разност- ных схем, выполнять оценки локальной и полной погрешности по- лученных результатов. 22 1.2.2. Интерфейс программы Программа написана на языке C++ с помощью пакета Microsoft Visual C++ 6.0 в виде диалогового приложения Windows. На рис.1.1 показано диалоговое окно программы в момент ее за- пуска. В программе независимая переменная описывается как вре- мя процесса. После выбора параметров рассматриваемой задачи и параметров ее численного решения, задаваемых в левой части диалогового ок- на программы, нажатие на кнопку «Расчет» приводит к появлению в правой части окна результатов расчета. Они выводятся как в виде графика изменения искомой функции Т от независимой перемен- ной t, так и в виде таблицы для просмотра полученных значений Т во всех узлах расчетной сетки (рис. 1.2). В случае большого коли- чества узлов сетки просмотр всей таблицы осуществляется с по- мощью вертикальной полосы прокрутки, расположенной справа от таблицы. Рис. 1.1. Диалоговое окно программы в момент ее запуска 23 Рис. 1.2. Диалоговое окно программы после выполненного расчета При неполном задании всех необходимых входных параметров программы появляется предупреждающее сообщение, показанное на рис. 1.3, а проведение расчета блокируется. Рис. 1.3. Предупреждающее сообщение при неполном задании входных параметров программы 24 Изменение любого входного параметра программы приводит к очистке диалогового поля «Результаты расчета», поэтому пред- ставленные результаты всегда соответствуют отображаемым на экране входным параметрам. 1.3. Варианты лабораторных работ 1.3.1. Исследование сходимости методов разностного решения ОДУ Цель работы Определение скорости сходимости различных методов числен- ного решения ОДУ и влияния на нее порядка аппроксимации раз- ностной схемы. Задачи работы 1. В соответствии с индивидуальным заданием, выданным препо- давателем, получить разностное решение заданной задачи Ко- ши с помощью нескольких разностных схем при различных шагах изменения независимой переменной. 2. Для каждого из рассмотренных методов численного решения определить погрешность расчета для всех исследованных ша- гов изменения независимой переменной путем сравнения с точным аналитическим решением. 3. Построить зависимости погрешности расчета от величины шага и проанализировать полученные результаты. 4. Оформить отчет о выполнении работы, в котором описать по- ставленную задачу Коши и используемые для ее решения раз- ностные схемы, привести результаты численного решения, по- строенные зависимости погрешности от величины шага и свои выводы по работе. Варианты задания Индивидуальный вариант задания содержит следующую ин- формацию. 1. Описание функции f(t, T) – ее вид и значения коэффициентов. 2. Начальное значение искомой величины. 25 3. Значение независимой переменной, при котором анализируют- ся результаты (конечное время процесса в программе). 4. Список используемых методов численного решения. 5. Набор шагов изменения независимой переменной, используе- мый в расчетах. 6. Точное аналитическое решение рассматриваемой задачи при заданном значении независимой переменной. 1.3.2. Исследование устойчивости методов разностного решения ОДУ Цель работы Определение границы устойчивости различных методов чис- ленного решения ОДУ. Задачи работы 1. В соответствии с индивидуальным заданием, выданным препо- давателем, получить разностное решение заданной задачи Ко- ши с помощью нескольких разностных схем при изменении шага независимой переменной в указанном диапазоне. 2. Для каждого из рассмотренных методов численного решения определить значение шага независимой переменной, при кото- ром наблюдается потеря устойчивости численного решения. 3. Построить зависимости искомой функции T(t) полученные при шагах изменения независимой переменной 99%, 100% и 101% от порогового значения. 4. Оформить отчет о выполнении работы, в котором описать по- ставленную задачу Коши и используемые для ее решения раз- ностные схемы, привести результаты численного решения, по- лученные пороговые значения шага и свои выводы по работе. Варианты задания Индивидуальный вариант задания содержит следующую ин- формацию. 1. Описание функции f(t, T) – ее вид и значения коэффициентов. 2. Начальное значение искомой величины. 3. Конечное значение независимой переменной, до которого про- водится расчет (конечное время процесса в программе). 26 4. Список используемых методов численного решения. 5. Диапазон изменения шага независимой переменной, исполь- зуемый в расчетах. 1.3.3. Оценка погрешности методов разностного решения ОДУ Цель работы Оценка локальной и полной погрешности результатов числен- ного решения ОДУ и определения их связи с истинной погрешно- стью в зависимости от используемой разностной схемы. Задачи работы 1. В соответствии с индивидуальным заданием, выданным препо- давателем, получить разностное решение заданной задачи Ко- ши с помощью нескольких разностных схем при двух значени- ях шага изменения независимой переменной. 2. Оценить полную и локальную погрешность расчета на первом шаге изменения независимой переменной и в момент заверше- ния расчета для каждого из рассмотренных методов численного решения. 3. Для каждого из рассмотренных методов численного решения определить истинную погрешность расчета при меньшем из ис- следованных шагов путем сравнения с точным аналитическим решением при значении независимой переменной, соответст- вующем первому большему шагу и моменту завершения расче- та. 4. Сравнить полученные значения локальной, полной и истинной погрешности. 5. Оформить отчет о выполнении работы, в котором описать по- ставленную задачу Коши и используемые для ее решения раз- ностные схемы, привести результаты численного решения, сравнения погрешностей и свои выводы по работе. Варианты задания Индивидуальный вариант задания содержит следующую ин- формацию. 1. Описание функции f(t, T) – ее вид и значения коэффициентов. 27 2. Начальное значение искомой величины. 3. Конечное значение независимой переменной, до которого про- водится расчет (конечное время процесса в программе). 4. Список используемых методов численного решения. 5. Набор шагов изменения независимой переменной, используе- мый в расчетах. 6. Точное аналитическое решение рассматриваемой задачи при заданном значении независимой переменной. 1.4. Контрольные вопросы 1. Что такое сходимость, аппроксимация и устойчивость разност- ной схемы? 2. Каковы достоинства и недостатки одношаговых и многошаго- вых разностных схем решения ОДУ? 3. Что такое схема Эйлера? 4. Как получить расчетные формулы исправленного метода Эйле- ра? 5. Какой порядок аппроксимации имеет модифицированный ме- тод Эйлера? 6. Что такое метод Рунге-Кутта 4-го порядка аппроксимации? 7. Как получаются расчетные формулы явной схемы Адамса 2-го порядка аппроксимации? 8. Каков порядок аппроксимирующего полинома при построении неявной схемы Адамса 2-го порядка аппроксимации? 9. Что такое явная схема Адамса 4-го порядка аппроксимации? 10. Что такое неявная схема Адамса 4-го порядка аппроксимации? 11. Как выбирается число итераций в методе прогноза-коррекции? 12. Как получаются расчетные формулы метода прогноза- коррекции на основе схем Адамса 2-го порядка аппроксима- ции? 13. Какие расчетные формулы используются в методе прогноза- коррекции на основе схем Адамса 4-го порядка аппроксима- ции? 14. Как оценить погрешность численного решения? 15. Каковы основные возможности учебной программы ODY_lab? 28 2. КОНЕЧНО-РАЗНОСТНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ТЕПЛООБМЕНА Пространственно-временные поля температур, скоростей, дав- лений и других теплофизических параметров при моделировании процессов тепломассопереноса определяются из решения краевых задач для дифференциальных уравнений в частных производных в заданных областях изменения пространственных переменных и на заданных временных интервалах. В настоящее время применение вычислительной техники и численных методов позволяет получать приближенные решения многомерных, нелинейных, нестационар- ных задач теплообмена, для которых использование точных и при- ближенных аналитических методов не представляется возможным. При использовании численных методов на первом этапе реше- ния задачи выполняется дискретизация пространственной и вре- менной областей, в ходе которой в этих областях задаются узловые точки. На втором этапе составляется система алгебраических урав- нений относительно значений искомых функций в этих узловых точках. Полученная система уравнений решается на третьем этапе. Основными численными методами решения уравнений в част- ных производных являются: метод конечных разностей и метод конечных элементов. Они отличаются способами получения сис- темы уравнений для значений искомых функций в узловых точках. Метод конечных разностей базируется непосредственно на диффе- ренциальном уравнении и граничных условиях, а метод конечных элементов — на эквивалентной вариационной постановке задачи. В данном разделе на примере уравнения теплопроводности рас- смотрим основные понятия теории численных методов, свойства и методы создания конечно-разностных схем. 2.1. Основные понятия теории разностных схем 2.1.1. Разностная схема и разностное решение Основные понятия теории разностных схем разберем на приме- ре одномерного нестационарного уравнения теплопроводности для пластины с внутренним источником тепловыделения 29 0 , 0 , max 2 2 t t l x q x T t T c v ≤ < < < + ∂ ∂ λ = ∂ ∂ ρ (2.1) На границах пластины заданы граничные условия третьего рода , , 0 , 0 , 0 l l x l q T x T = ⎥⎦ ⎤ ⎢⎣ ⎡ α + ∂ ∂ λ = m (2.2) а начальное условие имеет вид ) ( ) , ( 0 0 | x T t x T t = = (2.3) Решением задачи (2.1) – (2.3) является функция ) , ( t x T , заданная в непрерывной области { } { } max 0 0 t t l x ≤ ≤ × ≤ ≤ = Ω При использовании численных методов в пространственной об- ласти выбирается некоторое конечное число значений координаты x 1 , x 2 , …, x N (узлы пространственной сетки), для временной пере- менной также выбирается конечное число значений t 0 , t 1 , …, t J (уз- лов временной сетки). Целью является определение значений тем- пературы j n T в узлах пространственной сетки х n в моменты време- ни t j : T n j = T(x n , t j ), n = 1, ... N; j = 0, ... J, (2.4) т.е. значения искомой функции находятся в дискретной области t h ∆ Ω , (рис. 2.1): { } { } , , , , 0 1 , J N t h t t x x K K × = Ω ∆ 30 Для упрощения будем считать пространственное и временное разбиения равномерными с шагами h по координате х и ∆t по вре- мени: х n = (п – 1)h, h = l/(N – 1), n = l, ..., N; t j = j ∆t, ∆t = t max /J, j = 0, ..., J. Рис. 2.1. Пространственно-временная сетка, используемая при численном решении Уравнения для определения j n T получим из основной задачи (2.1) – (2.3). Выразим производные дТ/дх и д 2 Т/дх 2 в точке (х п , t j )через зна- чения функции j n T в этой точке и в некоторых соседних узлах сет- ки. Из определения производной имеем , ) ( 1 t t T T t T j n j n j n ∆ δ + ∆ − = ∂ ∂ − (2.5) где ) ( t j n ∆ δ – величина, стремящаяся к нулю при 0 → ∆t . Конечную разность 1 − − j n j n T T называют «разностью назад» или левой разно- стью. 0 = x 1 x 2 …, x n-1 x 1 x n x N = l t 1 t 2 t j ∆t h (x n , t j ) t max = t J : : :: t x 31 Выражение для ) ( t j n ∆ δ можно получить, выразив в (2.5) 1 − j n T с помощью разложения в ряд Тейлора относительно точки (х п , t j ): 2 2 2 2 1 K − ∆ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∆ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ − = − t t T t t T T T j n j n j n j n После несложных преобразований получим: 6 2 ) ( 2 3 3 2 2 K + ∆ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ − ∆ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ = ∆ δ t t T t t T t j n j n j n (2.6) Из выражения (2.6) следует, что при достаточно малых ∆t вы- полняется неравенство ) ( 1 t А t j n ∆ ≤ ∆ δ Здесь А 1 – константа. Таким образом, ) ( t j n ∆ δ является величиной порядка ∆t, т.е. o(∆t). Аналогичным образом можно построить аппроксимацию для временной производной с помощью «разности вперед» (или правой разности): ) ( 1 t o t T T t T j n j n ∆ + ∆ − = ∂ ∂ + (2.7) Достоинства и недостатки этих двух способов аппроксимации производной дТ/дt рассмотрим позднее. При построении выражения для второй производной д 2 Т/дх 2 ис- пользуем значения искомой функции j n T в трех соседних узлах пространственной сетки: 32 ). ( / ) 2 ( ) ( / ) ( / ) ( 2 1 1 1 1 2 2 h h T T T h h h T T h T T x T x x T j n j n j n j n j n j n j n j n j n j n j n γ + + − = γ + − − − = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + − − + (2.8) Используя представление j n T 1 − и j n T 1 + с помощью рядов Тейло- ра относительно точки ( х п , t j ), можно показать, что в выражении (2.8) ) ( ) ( 2 h o h j n = γ Подставив выражения для производных (2.5) и (2.8) в уравнение теплопроводности (2.1), получим с q h T T T с t T T v j n j n j n j n j n j n j n ρ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ γ + + − ρ λ = δ + ∆ − − + − 2 1 1 1 2 . (2.9) Отбрасывая величины j n δ и j n γ , получим уравнение для опреде- ления приближенных значений j n u искомой величины j n T ( ) с q u u u h a t u u v j n j n j n j n j n ρ + + − = ∆ − − + − 1 1 2 1 2 . (2.10) Здесь a =λ/ρc – коэффициент температуропроводности. Уравнения (2.10) можно записать для всех внутренних про- странственных узлов (п = 2, ..., N – 1). Уравнения для j u 1 и j N u по- лучим из граничных условий (2.2). Заменяя производные их разно- стными аналогами и отбрасывая члены порядка о(h), получим l j N l j N j N j j j q u h u u q u h u u = α + − λ = α + − λ − −1 0 1 0 1 2 , . (2.11) 33 В начальный момент времени из начального условия (2.3) име- ем: N n x T u n n , , 1 , ) ( 0 0 K = = . (2.12) В численных методах дискретное множество { } N n n x 1 = называется пространственной сеткой, дискретное множество { } J j j t 0 = – времен- ной сеткой,дискретное множество (область) t h ∆ Ω , – пространст- венно-временной сеткой.Совокупность значений ) , ( j n j n t x T T = в узлах пространственно-временной сетки называется сеточной функцией точного решения.Совокупность приближенных значе- ний j n u называется сеточной функцией разностного решения или просто разностным решением. Разница между j n T и j n u называется погрешностью разностного (численного) решения. Обозначим ее через j n j n j n u T − = ε . Система алгебраических уравнений (2.10) – (2.12), соответствующая исходной дифференциальной задаче (2.1) – (2.3), называется разностной схемой. |