Метод конечных элементов в задачах теплопроводности (стр. 7 )

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

, , (6.2.4)

Применив метод Галеркина к дифференциальному уравнению (6.2.1), для -го временнόго элемента получим систему уравнений в интегральной форме:

. (6.2.5)

Рис. 6.1

(6.2.3)

где – матричная строка базисных временных функций; – вектор-столбец всей совокупности значений температуры в -ом и -oм узлах -го временного элемента.

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

, . (6.2.4)

Применив метод Галеркина к дифференциальному уравнению (6.2.1), для -го элемента получим систему уравнений в интегральной форме:

.

Интегрирование этого выражения с помощью -координат приводит к системе алгебраических линейных уравнений:

Разрешим ее относительно {}:

.

Вводя обозначения:

, , (6.2.64)

последнее выражение запишем компактно:

. (6.2.75)

Уравнение (6.2.64) содержит две неизвестные – совокупности значений температур в -ом и в -ом узлах -го элемента. Для обеспечения непрерывности интерполяционной функции (6.2.3) на границе между в общем для соседних элементов узлеами должно выполняться условие:

. (6.2.86)

Подставляя (6.2.86) в (6.2.75) и опуская одинаковый для всех членов уравнения индекс , получим рекуррентное уравнение, позволяющее выразить координатные узловые значения температуры на -ом временнόом шаге через совокупность их значений на предыдущем – -м шаге:

. (6.2.97)

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

На первом временнόом шаге в качестве будет фигурировать, очевидно, начальное – задаваемое – условие (6.2.2) в -ом временном узле, что и позволяет найти по уравнению (6.2.95) значения температур в -ом узле этого же – с – временного элемента:

. (6.2.108)

Определение по уравнению (6.2.108) и, тем самым, согласно (6.2.86) и значений в -ом узле второго временнόого элемента с , позволяет организовать последовательный процесс в соответствии с рекуррентным уравнением (6.2.97).

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

Другие способы решения системы динамических уравнений (6.2.1) описаны в [2, 6] (решение методом конечных разностей).

6.З Учет температурной зависимости теплофизических параметров

Полученное в п. 6.2 выражение (6.2.107) является решением линейного уравнения теплопроводности с линейными граничными условиями, так как все теплофизические параметры задачи полагались постоянными. Это ограничение общности полученного решения может быть снято программными методами [1].

При небольшой величине временного шага изменения температуры на -ом шаге будут невелики,; в силу чего параметры задачи можно считать постоянными на этом временном элементе. Это тем более справедливо, если учесть, что температурная зависимость параметров довольно слабая – (через а обозначен теплофизический параметр) невелико. Однако постепенно изменения температуры будут накапливаться и на каком-то временнόом шаге значения параметров следует перевычислить соответственно средней температуре элемента по заложенной в физическом каталоге процедуре:

(6.3.1)

Здесь под понимаются: коэффициент теплопроводности ; объемная теплоемкость , ; .

Для установления момента времени, в который нужно производить перерасчет параметров, на каждом шаге , начиная с , достаточно сравнить среднеобъемную температуру элемента (см.(5.4.1)) с ее значением на первом временном отрезке:

. (6.3.2)

По задаваемой величине , например, ≈ 100К, можно теперь определить номер к-го временнόого шага, с которого нужно перевычислить по (6.3.1) параметры тех элементов e′, для которых условие оказалось выполненным. Для этих элементов находятся разности:

,

,

и в стандартизованных матрицах элементов e′коэффициенты заменяются на , что позволяет найти поправочные числовые матрицы и .Все найденные указанным способом поправочные матрицы заносятся по известной процедуре (см. п. 5.2) в глобальные матрицы и , что приводит к поправленным на температурную зависимость параметров матрицам и . Относительно поправки объемной теплоемкости следует иметь в виду, что температурная зависимость ср и ρ носит противоположный характер: если ср растет с температурой, то ρ – падает, так что их произведение крайне слабо зависит от температуры. Поэтому поправлять матрицу теплоемкости на температурную зависимость едва ли целесообразно.

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

и .

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

,

который так же находится на каждом .

Таким образом, рекуррентное уравнение (6.2.107) с внесенными поправками на температурную зависимость параметров, по которому нужно проводить расчет температур на следующем – шаге, будет иметь вид [1]:

. (6.3.3)

На последующих временных шагах температура "поправленных"

На последующих временных шагах температура "поправленных" элементов сравнивается с их же температурой на шаге, а температура остальных элементов – с их температурой на шаге. Следовательно, неравенство (6.3.2) после -го шага расщепляется на два неравенства:

,

(6.3.4)

.

При выполнении любого из них процедура внесения поправок повторяется.

6.45 Радиационный компонент теплообмена

Конвективный компонент теплообмена характеризуется его интенсивностью – коэффициентом теплоотдачи α. Его зависимость от температуры довольно слабая, поскольку температура отнесения определяется как средняя температур теплоносителя и обтекаемого тела – , которые влияют на числа Рейнольдса и Нуссельта – параметры конвективного процесса теплоотдачи – лишь через слабую зависимость от температуры коэффициента теплопроводности и динамической вязкости газового теплоносителя. Поправка лежит в пределах 5÷10% и ею можно пренебречь.

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

Вторая отличительная особенность переноса излучения – не обязательное наличие среды для обмена энергией излучения между участками, в то время как для переноса энергии конвекцией или теплопроводностью присутствие физической среды обязательно. Если среда отсутствует, то излучение становится единственноым возможным способом переноса тепла. Так, например, для отвода отработанного тепла от космических энергетических установок может быть использовано только излучение [19, 27, 32, 47 – -52].

Радиационный теплообмен с окружающей средой внешних элементов в диффузно-сером приближении можно описать выражением [1, 16]:

, , (6.45.1)

где – постоянная Стефана-Больцмана ; – интегральная

65

степень черноты; – среднеповерхностная температура; – номер внешней "радиационной" поверхности элемента.

Радиационный теплообмен между и элементами – внутренний теплообмен, описывается более сложной формулой, основанной на зональном методе расчета лучистого теплообмена в замкнутой системе изотермических поверхностей [1, 16]:

, (6.45.2)

где – средний разрешающий угловой коэффициент излучения, величина которого находится решением системы уравнений:

, (6.45.3)

включающих известные средние геометрические угловые коэффициенты излучения . В частности, для камеры стационарного плазменного двигателя (СПД) они даны в работах [17, 18]. Под температурами Tm и Tn в (6.4.2) понимаются среднеповерхностные температуры этих поверхностей.

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

; .

Процедура учета температурной зависимости аналогична описанной в п. 6.3.

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

66

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

Сравнение температур (по столбцам) при разных показывает, что они

существенно падают при одновременном росте их градиента в радиальном направлении. Большие различия температур в первой и второй строках свидетель-

Таблица 6

ε r	Температура в узлах, К
1	2	3	4	5	6	7	8	9
0,0 0,3 0,9	710 592 556	708 590 553	705 587 551	646 525 488	643 521 483	640 517 477	600 470 432	598 467 428	595 464 422

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

Таблица 7

Величина температурного градиента

τ, мин.	Т1, К α ≠ 0, ε = 0	∂Т1/∂τ	Т2, К α ≠ 0, ε ≠ 0	∂Т2/∂τ
0,5 1,0 1,5 2,0 2,5	342,9 352,8 354/9 355,3 355,4	20,2 4,2 3,2 0,2	297,3 307,6 317,1 323,8 324,1	20,6 19,0 13,4 0,6

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

6.5 Сходимость, полнота и согласованность, точность

и пространственная изотропия

Остановимся кратко на некоторых вопросах математического характера.

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

Точность – это мера близости численного решения к точному, или истинному. Устойчивость определяется ростом ошибок при выполнении отдельных вычислительных операций. Неустойчивое вычисление является результатом аппроксимации, округления или других ошибок, которые неограниченно накапливаются, вследствие чего истинное решение тонет в ошибках. Сходимость – это постепенное приближение последовательно вычисляемых решений к предельному по мере того, как уточняются некоторые вычислительные параметры, такие, как размер элемента или ранг базисной функции в пробном решении. Термин “сходимость“ в этом же смысле применяется и к итерационной процедуре, в которой некоторые или все результаты одного вычисления становятся входной информацией для другого (повторного) вычисления. В сходящейся процедуре разница между последовательными результатами уменьшается, стремясь в пределе к нулю. Более подробные сведения по этим вопросам можно найти в курсах по численному анализу и методам вычислений.

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

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

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

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