Методы решения задач теплопроводности

Формула (11) позволяет вычислить значения в узлах горизонтального ряда (k+1) по значениям , находящимся только в одном предшествующем ряду (k). Поэтому с помощью формулы (11) можно найти значения узлах первого горизонтального ряда (при ) по известным из краевых условий значениям температуры в узлах самой оси 0x (при ). Получив таким образом значения в первом ряду, по той же формуле находим значения в узлах второго горизонтального ряда (т. е. при ). Этот процесс построения можно продолжать как угодно далеко, так как значения температуры в узлах прямых x=0 и x=L будут известны из граничных условий.

Формулу (11) можно вывести, применяя законы Фурье и Ньютона к составлению тепловых балансов элементов, на которые разбито тело. Этот способ вывода расчетных формул широко практикуется в настоящее время многими зарубежными исследователями.

Перепишем, следуя Д. Ю. Панову, формулу (11) в более удобном виде:

(12)

Выбирая соотношения между шагами l и h различным образом, из (12) можно получить ряд частных соотношений. Так, например, при

(13)

при

(14)

при

(15)

и вообще при

 (16)

Особенно простой вид формула (12) получает при p=2:

(17)

Эта последняя формула, называемая формулой Э. Шмидта, была впервые получена Д. Ю. Пановым в 1938 г. и имеет большое практическое преимущество по сравнению с формулой (14) и особенно (15), так как при том же h величина l самая большая, а следовательно, и объем вычислительной работы при применении (17) сокращается в 3 раза по сравнению с (14) и в 6 раз по сравнению с (15). Следует отметить, из-за своей простоты формула (17) широко используется при графическом решении нестационарных задач переноса.

Исследования показывают, что при p=1 получаем уже расходящуюся вычислительную схему. Вообще следует отметить, что при решении нестационарных уравнений в частных производных параболического типа вопросы соотношения между h и l, а также ошибка округления в численном решении играют первостепенную роль, ибо ими определяются сходимость и устойчивость получаемых решений. Из строгих теоретических рассуждений следует:

1) пользоваться формулой (16) можно только при ;

2) самый большой шаг l дает формула (17), т. е. при p=2;

3) формула (16) тем точнее, чем больше p.

Рассмотренная сетка для численного интегрирования уравнения (6) удобна, когда задача решается при граничных условиях первого рода: в этом случае граничные прямые x=0 и x=L принадлежат самой сетке. Если уравнение решается при граничных условиях третьего рода, практика вычислений и теоретические исследования показывают, что для повышения точности определения потенциала на границах следует вводить дополнительные узловые точки, лежащие вне изучаемой области. Например, решая уравнение (6) при граничных условиях

, (18)

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

(19)

а значение температуры на самой поверхности, т. е. T(L,kl), берем как среднее арифметическое значение температур в точках A и C:

(20)

Тогда условие (18) запишется так:

(21)

или переходя к приближенным значениям , получим после преобразований

(22)

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

Значение же температуры на самой граничной прямой x=L определяется по формуле

(23)

что после преобразований дает

(24)

где

(25)

Значение температуры в узлах вспомогательной прямой находится по формуле

 (26)

На левой поверхности получим

(27)

Методика решения дифференциального уравнения теплопроводности с источниками не отличается от изложенной выше. Метод конечных разностей позволяет успешно решать как одномерные, так и двух- и трехмерные задачи. Случай, когда на область изменения переменных x и y наносится квадратная сетка, полностью исследован Ш. Е. Микеладзе. Треугольные и полярные сетки рассмотрены П. П. Юшковым. Необходимо отметить, что полярные сетки особенно удобны для решения задач с осевой симметрией.

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

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

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

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

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

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

(28)

в то время как - производная по времени аппроксимируется по прежней формуле (7).

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

Известно, что уравнение

 (29)

устойчиво и его решение при , сводится к решению дифференциального уравнения теплопроводности.

Покажем только, как в случае одной пространственной переменной решать систему уравнений на каждом слое. Нетрудно убедится, что (29) для неизвестных можно переписать в виде

(30)

где образованы из и значений , которые известны.

Пусть краевые условия (первого, второго или третьего рода) имеют вид

(31)

при этом и . Будем искать зависимость вида

(32)

которой должно удовлетворять решение (30).

Определяя при помощи (32) и используя (30), имеем

(33)

Потребуем, чтобы (33) было следствием (31); так как коэффициенты при одинаковых и свободные члены должны быть пропорциональны, будем иметь

; (34)

По этим формулам можно найти (i=N-1, N-2, …, 1), начиная с . При этом будет

Зная , вычисляем q_i (i=N-1, …, 1; q_N=).

Взяв (32) при i=1 и первое из соотношений (31), решим систему уравнений для и ; неравенства и обеспечивают ее разрешимость. По соотношениям (32) находим последовательно все остальные . Ввиду полученных для оценок в этих операциях не будет происходить существенного накопления и роста погрешностей. Этот метод носит название метода прогонки.

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

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

3.1.3 Условие устойчивости

Наиболее естественной конечно-разностной задачей для аппроксимации дифференциальной задачи:

является следующая [9]:

(35)

Если положить

(36)

то разностное уравнение можно переписать так:

 (37)

Исходя из начальных значений , можно последовательно вычислять значения в узловых точках временных рядов k=1, 2, … с помощью формулы (37).

Если

(38)

то правая часть уравнения (37) является линейной комбинацией значений с неотрицательными коэффициентами, сумма которых равна единице. Следовательно, , как и , во всех узлах сетки удовлетворяет неравенствам

, если (39)

коль скоро .

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

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

Пусть

Тогда решения уравнения (35), соответствующее этим начальным условиям, представляет собой отклонение, вызванное отдельной погрешностью величины в точке x=0. Если , то мы докажем по индукции, что знак в узлах будет чередоваться, т. е. если отлично от нуля в двух соседних узловых точках, то значения в этих точках будут противоположного знака.

Это утверждение непосредственно проверяется при k=0 и k=l. Допустим теперь, что оно уже проведено для всех узлов сетки, лежащих под и на временном ряду k, и рассмотрим узел (i, k), в котором . Тогда правая часть уравнения (37) положительна, т. е. . В случае рассуждения аналогичны. Этим индукция завершается.

Ввиду чередования знака соотношение (37) дает

(40)

Пусть

где сумма распространена на все узлы ряда k. Так как самое большее этих членов отличны от нуля, то сумма S(t) конечна. Согласно (40), мы имеем

С другой стороны, должен существовать по крайней мере один узел (n,k) временного ряда k, в котором

Правая часть этого неравенства экспоненциально растет при , потому что . Иными словами, разностное уравнение (35) неустойчиво при .

3.1.4 Сходимость и погрешность метода

Исследуем сходимость к [6]. Допустим сначала, что при , функция T имеет непрерывные ограниченные частные производные до третьего порядка по и до шестого порядка по x. Так как , то эти два условия для данного простейшего дифференциального уравнения эквивалентны. Многие формулы, приводимые ниже, могут быть упрощены подстановкой вместо . Однако это уже невозможно для других простых дифференциальных уравнений, которые можно было исследовать аналогичным методом. Пусть обозначает левую часть равенства (35). По формуле Тейлора и в силу того, что мы получаем

 (41)

Вычитая это из , мы видим, что погрешность метода является решением разностной краевой задачи

(42)

где

(43)

Разностное уравнение в (42) можно переписать так:

(44)

Если обозначает верхнюю грань для при , , и если - аналогичная верхняя грань для , то из уравнения (44) при получается неравенство: , = 0, k, 2k, …. Так как m(0)=0 и так как - неубывающая функция ф, то по индукции получаем, что , и поэтому

(45)

Эту оценку погрешности метода можно использовать независимо от вида . Опираясь на представление (43), мы заключаем что погрешность метода при сделанных предположениях имеет точный порядок , за исключением случая , когда она имеет порядок . Это доказывает сходимость к .

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

4. Расчет температурного поля пластины

Рассмотрим практическую задачу, используемую в инженерных расчетах - распространение тепла в пластине с течением времени.

Условия задачи следующие: Определить нестационарные температурные поля в неограниченной пластине (алюминиевый сплав ) с теплофизическими свойствами: , , , если на одной поверхности, при x=L, задано изменение температуры со временем , , а другая поверхность теплоизолирована. Толщина пластины: L=10 см. Начальное распределение на всей пластине равно .

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

(1)

где a - коэффициент теплопроводности определятся теплофизическими свойствами тела .

Начальные условия:

(2)

Граничные условия:

(3)

 (4)

Формула (4) - описывает теплоизолированную поверхность, где нулевой слой взят за пределами пластины и на этом слое температура всегда равна .

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

А решение в интегральной форме приводилось выше, и имело менее громоздкий результат (п. 3.3, формулы (2.36)-(2.39)), но имело существенный минус, так как довольно сложно составить программу нахождения оригинала по изображению.

Для соблюдения поставленных критериев решения удобнее воспользоваться численными методами решения, в частности методом конечных разностей (МКР).

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

(1а)

(2а)

 (3а)

(4а)

при , где .

Данный метод обычно требует определенного количества однообразных вычислительных операций. Поэтому была разработана программа Teplopr (см. Приложение 1).

Программа Teplopr - написана на языке Pascal и работает по явной схеме. Программа имеет довольно простую структуру циклов (см. Приложение 2) и проста в применении.

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

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

Для примера можно ввести следующие данные: dx=1 (h=1), dt=100 (l=100). Тогда kF=0.3 (т. е. ), и мы получаем приемлемый результат (см. Приложение 3). Дальнейшие расчеты (после получасового промежутка), при постоянном граничном условии стороны x=L, показали, что в процессе нагревания пластины происходит ее равномерное нагревание и она выходит на стационарный режим.

Также можно попробовать ввести данные заведомо неудовлетворяющие критерию устойчивости: dx=1(h=0.7), dt=100(l=100). Тогда kF=0,612 (т. е. ), и мы получаем совершенно противоречивые результаты (см. Приложение 4). Также в приложение 5 дана сравнительная диаграмма двух полученных результатов.

Заключение

В дипломной работе была рассмотрена математическая постановка задач нестационарной теплопроводности.

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

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

Составлена программа Teplopr численного расчета температурного поля пластины. В таблице приведены значения температуры для любого сечения пластины. На графике наглядно видно как изменяется температура вдоль пластины в зависимости от времени прогрева поверхности.

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

Список использованной литературы

1. Беляев Н.М., Рядно А.А. «Метод нестационарной теплопроводности» М. Высшая школа. 1978. 328с.

2. Лыков А.В. «Теория теплопроводности». М. Высшая школа. 1967. 600с.

3. Патанкар С.В. «Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах.» М. Издательство МЭИ, 2003. 312с.

4. Пискунов Н.С. «Дифференциальное и интегральное исчисления» (том II). М. Интеграл-пресс. 2002. 410с.

5. Рихтмайер Р. Мортон К. «Разностные методы решения краевых задач.» М. Издательство Мир. 1972. 380с.

6. Сабитов К.Б. «Уравнения математической физики.» М. Высшая школа. 2003. 255с.

7. Форсайт Дж., Вазов В. «Разностные методы решения дифференциальных уравнений в частных производных.» М. Иностранная литература. 1963. 488с.

8. Юшков П.П. «Приближенное решение задач нестационарной теплопроводности методом конечных разностей.» Труды Института энергетики АН БССР (выпуск 6.). 1958. 203с.

9. Яненко Н.Н. «Метод дробных шагов многомерных задач математической физики.» Новосибирск. Наука. 1967. 195с.

Размещено на Allbest.ru