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

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

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

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

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

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

- способ численного решения дифференциального уравнения процесса;

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

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

Кроме того, программа выводит на печать или через графические средства результаты расчета.

Описание метода численных решений систем дифференциальных уравнений, подходов к созданию типовых программ приведены в [10, 11, 12]

Постановка задачи: определить изменение температур в теле во времени при задании по аналогии с аналитическими методами исходных данных:

- геометрия тела;

- теплофизические свойства материала тела;

- теплофизические свойства окружающей среды;

- граничные условия, определяющие характер теплообмена (I, II, III, или IV рода);

- начальное температурное состояние тела;

- продолжительность во времени процесса.

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

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

Задача называется устойчивой по входным данным, если решение непрерывно зависит от входных данных, т.е. при небольшой погрешности δx входных данных погрешность в решении также будет невелика. В противном случае решение неустойчиво.

При замене исходной краевой задачи конечно-разностной необходимо рассмотреть вопросы: 1) аппроксимации, т.е. метода приближения; 2) устойчивости конечно-разностного метода (алгоритма) решения задачи; 3) сходимости, т.е. приближения решения конечно-разностной задачи к решению исходной аналитической задачи при изменении так называемых шагов сеточной области.

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

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

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

Метод конечных элементов (МКЭ) оказался весьма эффективным численным методом решения большого круга физических и инженерных задач.

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

Идея метода Эйлера состоит в том, что значения, например, функционала

(4.53)

рассматриваются не на произвольных, непрерывных и допустимых в данной вариационной задаче кривых y = y(x), а на всевозможных ломаных, составленных из n прямолинейных звеньев с заданными абсциссами вершин, которые называются узлами:

а + Δх, а + 2Δх, ...., а + nΔх, где

, а ≤ х ≤ b

Если ординаты вершин ломаной линии обозначить через y(a), y₁, y₂ ... y_n, y(b), то функционал I(у) на таких ломаных становится функцией ординат y₁, y₂ ... y_n вершин этих ломаных, т.е

I = I(у₁, у₂, ..., у_n ), (4.54)

а так как ломаная вполне определяется заданием узловых ординат, то производные у^/ также могут быть выражены через конечные разности значений функции в узлах, а сам интеграл можно приближенно заменить суммой отдельных значений F, умноженных на Δx.

Таким образом, задача о минимуме функционала (4.53) сведена к хорошо разработанной задаче о минимуме функций многих переменных – узловых значений y₁, y₂ ... y_n. Для решения этой задачи, как известно, нужно составить следующую систему:

 (4.55)

.............

В результате решения систем уравнений (4.55) относительно узловых значений y₁, y₂ ... y_n находятся эти значения, а тем самым и ломаная, доставляющая минимум функционалу (4.53).

В методе конечных элементов (МКЭ), как и конечно-разностном методе Эйлера, составляются уравнения типа (4.55), из которых определяются неизвестные узловые значения y₁, y₂ ... y_n разыскиваемой функции y = y(x). Эти уравнения могут составляться либо на основе вариационных методов, либо на основе методов взвешенных невязок.

Принципиальное отличие МКЭ от конечно-разностного метода Эйлера для одномерной задачи заключается, пожалуй в том, что на выбранном конечном элементе (участке Δx) функция у(х) не обязательно заменяется линейной аппроксимацией, а может быть принята в виде полинома

(4.56)

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

Решение задачи методом элементарных балансов (академика А.П. Ваничева).

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

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

- выделенный элементарный объем находится со всех сторон внутри тела, и баланс тепла для всех граней определяется уравнением Фурье.

Для составления в качестве примера расчетных уравнений примем, что выделенный участок лежит внутри тела. Простейшая форма расчетного объема – параллелепипед с гранями Δx,Δy,Δz. Расчетными точками являются места пересечения плоскостей разбивки, т.е. угловые точки параллелепипеда (см. рис.4.6).

Рис. 4.6. Схема разбивки тела на элементы.

Температуры в расчетных точках снабжают индексами, характеризующими время и место. Температуру конкретной рассматриваемой точки называют просто t. Температуры в данный момент времени в соседних точках, находящихся на расстоянии Δx,Δy,Δz, обозначают соответственно . Температура расчетной точки в последующий момент времени, т.е. через промежуток времени Δτ обозначают .

Процесс распространения тепла со всех сторон описывается законом Фурье и определяется следующими теплофизическими свойствами вещества: коэффициентом теплопроводности, удельной теплоемкостью и плотностью. Изменением плотности пренебрегаем. Изменение коэффициента теплопроводности и удельной теплоемкости принимается по линейному закону от температуры: λ = А + Вt и с_р=С + Dt. Учитывая малые размеры элементарного объема, делаются допущения:

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

- величина среднего за время Δτ теплового потока ΔQ через какую-либо поверхность пропорциональна начальному в пределах отрезка времени значению температурного градиента;

- увеличение энтальпии элемента пропорционально приращению температуры в средней точке его объема.

Для получения расчетной формулы составим тепловой баланс для элемента, температура в центральной точке которого является расчетной t и t_τ+Δτ . Количество тепла, вошедшее в элемент за время Δτ через левую грань, параллельную плоскости YOZ, т.е. грань, лежащую в плоскости, отстоящей на расстоянии x = -Δx/2, на основании закона Фурье равно:

А через противоположную грань за тот же отрезок времени поступит количество тепла:

Количество теплоты, вошедшее в элемент через четыре другие грани, параллельные плоскости XOY и XOZ , определяется аналогично:

;

;

;

;

В силу линейного характера изменения температур в пределах расчетного элемента справедливы равенства:

;

;

и т.д.

С учетом этих равенств выражения для ΔQ₁, ΔQ_{2, ....} ΔQ₆ могут быть записаны:

; 

;

;

;  (*)

;

. 

Алгебраическая сумма количества теплоты, вошедшей за время Δτ через все грани в элемент, равна увеличению его энтальпии. Это может быть выражено в виде равенства:

Подставляя в это выражение вместо ΔQ₁, ΔQ_{2, ....} ΔQ₆ выражения (*) и решая полученное уравнение относительно интересующего нас значения температуры в следующий момент времени t_τ+Δτ получаем:

, (4.57)

где

;

;

;

;

.

Пользуясь формулой (4.13), можно по известному начальному распределению температур последовательно найти значения температур во всех расчетных точках в моменты времени τ+Δτ, τ+2Δτ, τ+3Δτ и т.д. до достижения требуемого времени расчета. Расчет будет успешным, если выполнится еще одно условие – отрезок времени Δτ должен быть достаточно мал, иначе допущение о линейном распределении температур будет слишком грубым приближением, и ход расчета пойдет «вразнос».

Расчетная формула (4.13) может быть представлена в виде:

, (**)

где

; 

;

......................................................................  (***)

......................................................................

. 

Допустимую величину Δτ можно определить из рассмотрения системы (**) и (***).

Расчетные зависимости (**) определяются геометрическими параметрами, температурами в теле. Среди температур можно выбрать максимальное и минимальное на данный момент времени значение. Для обеспечения устойчивости решения необходимо, чтобы результат лежал между ними. Можно видеть, что коэффициенты А_2, А₃ ...А₇ меняются монотонно от Δτ и имеют положительное значение. А коэффициент А₁ при произвольном выборе Δτ может изменяться от +1 до -∞. Максимально допустимой величиной расчетного отрезка времени является значение Δτ_макс является то, когда А₁ обращается в 0. Для запаса его надо уменьшить. Учитывая, что из-за зависимости от температур свойств вещества, и А₁ зависит от распределения температур, рассчитываются значения Δτ_макс для двух предельных значений из массива температур:



 (****)



Из двух найденных значений надо выбрать расчетный отрезок времени меньше меньшего.

В данном варианте составления уравнений тепловых потоков принята так называемая я в н а я схема расчета. Опорными в этом методе являются нулевые (или предшествующие) значения температур в точках тела и определяются явным образом температуры через расчетный интервал времени. Существенным недостатком этого подхода является строгое выполнение требования по расчетному отрезку времени Δτ.

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

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

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

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

- метод термокрасок;

- метод выплавления материалов с разными температурами плавления;

- метод измерения термопарами;

- метод измерения облученными алмазами и другие

имеет точность измерения на уровне не выше нескольких процентов, т.е. не выше ± (5-10⁰).

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

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

Самое простое решение получается тогда, когда градиентом температуры внутри тела можно пренебречь, т.е. когда мы имеем дело с телом, имеющим коэффициент теплопроводности,стремящийся к бесконечности и температуры тела во всем сечении равны друг другу в любой рассматриваемый отрезок времени (говорят, что критерий Био мал – Bi = α δ/λ→ 0).

В этом случае вместо дифференциального уравнения Фурье можно рассматривать просто уравнение баланса тепла:

(*)

Тепло , подводимое при помощи конвекции к телу с боковой поверхностью А и объемом V, идет на нагревание тела.

Интегрируя уравнение *, получаем:

Постоянная В определяется из начального условия: при τ =0 Т=Т₀; отношение объема к поверхности тела V/А = R_V = kR. Здесь R – половина толщины пластины, радиус шара или цилиндра. Постоянная k зависит от геометрии тела:

- для неограниченной пластины k = 1;

- для неограниченного цилиндра и квадратного стержня k= 1/2 ;

- для куба и шара k = 1/3.

С учетом отмеченого выше решение задачи запишется так:

(**)

Следовательно, процесс нагревания (Т_f >Т) или охлаждения (Т_f <Т) описывается простой экспонентой; он полностью определяется коэффициентом теплообмена и удельной теплоемкостью тела.

Мгновенный поток тепла

(***)

Дифференцируя решение (**) и подставляя результат в (***), получаем:

Суммарный поток тепла (общий расход тепла) за время от 0 до τ составляет:

(****)

Следовательно, значения мгновенного и суммарного потоков тепла могутбыть получены из графиков рис.4.7.

Рис.4.7. Охлаждение или остывание тел.

1- неограниченная пластина; 2- неограниченный стержень

квадратного сечения; 3- неограниченный цилиндр; 4- куб;

5- конечный цилиндр ℓ⁄d =1; 6- шар.

В общем случае решение задачи нагревания тела простейшей геометрической формы в среде с постоянной температурой представляется в виде суммы решений:

(*****)

Здесь μ_i –корни характеристического уравнения. Свойство этих корней таково, что с повышением порядка i величина самого корня стремится к нулю. С увеличением временного фактора Fо , т.е с ростом продолжительности процесса члены уравнения с порядком i>1 быстро стремятся к нулю. Определять процесс изменения температуры будет первый член этой суммы. Общая картина изменения температур при остывании (нагреве) будет происходить в соответствии с картиной на рис. 4.8.

Рис. 4.8. Логарифм температурной разности

как функция времени в процессе нагревания.

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

На стадии регулярного режима имеет место постоянство темпа остывания в полулогарифмических координатах, т.е.:

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