Решение краевых задач методом конечных элементов

Постановка задачи. Выберем в качестве ИТО одномерный стержень с коэффициентом теплопроводности l, показанный на рисунке 11.1-а. Стержень имеет теплоизолированную боковую поверхность. К левому концу стержня подводится тепловой поток заданной интенсивности q (Вт/см²). На правом конце стержня происходит конвективный обмен тепла с коэффициентом теплообмена – h (Вт/см^{2 о}С). Температура окружающей среды – Т_ос (^оС). Поскольку стержень теплоизолирован, потерь тепла через боковую поверхность не происходит. Требуется определить температурное поле вдоль стержня в установившемся режиме.

Известно, что для данной модели распределение температуры внутри стержня описывает следующее дифференциальное уравнение:

	l	д 2T	= 0	(11.1)
	дx 2
a)	б)
Рис. 10.1

При этом, поскольку в установившемся режиме в точках приложения (при х =0) и отвода (х =L) тепла тепловая энергия не должна «задерживаться», должны быть соблюдены следующие граничные условия:

– на левом конце стержня (х=0):

– на правом конце стержня (х=L):

Если тепло отводится от стержня, тепловой поток q должен быть положителен, в противном случае – отрицателен.

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

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

С физической точки зрения функционал (11.5) моделирует непрерывность теплового потока в установившемся тепловом режиме. Это означает, что в любой момент времени сумма подводимой (через поверхность S₁) к стержню и рассеиваемой им (через поверхность S₂) тепловой энергии равна энергии, сосредоточенной в объеме (V) стержня. В противном случае, не отводимый от стержня избыток тепловой энергии будет продолжать нагревать стержень, что противоречит условию установившегося режима.

Поскольку, с одной стороны, установившийся режим описывается дифференциальным уравнением (11.1) с граничными условиями (11.2 и 11.3), а, с другой стороны, функционал (11.4) достигает минимума именно в установившемся режиме, то минимум функционала (11.4) и является решением ДУ (11.1) с граничными условиями (11.3).

Температура стержня во всех точках сечения S₁ (S₂) одинакова и равна неизвестной пока (но постоянной в стационарном режиме) величине – Т₁ (Т₃). Учитывая, что в данном случае S₁ = S₂ = A и в силу сказанного, выражение (11.5) принимает вид:

Таким образом, исходное уравнение для определения температуры в каждой точке стержня методом МКЭ примет вид:

Реализация метода МКЭ включает этапы:

1. Определение подобластей (конечных элементов) и их узловых точек. В данном случае, стержень может быть разбит на два одномерных симплекс – элемента, как это показано на рисунке (10.1-б) с узловыми значениями Т₁, Т₂ и Т₃. Температура внутри элементов находится из формул:

ФФ здесь согласно (9.5) равны:

2. Вычисление частных производных, входящих в выражение (11.7):

3. Разделение интеграла в выражении (11.7) на два (по числу подобластей – конечных элементов, выделенных в пункте 1). Необходимость разбиения интеграла продиктована тем, что производная температуры по переменной х (градиент температуры по оси ОХ), входящая под знак интеграла, не является непрерывной в точке Т₃. Учитывая, что dV=Adx, где А – площадь сечения стержня (А₁ = А₂ = А₃ =А), после разделения и подстановки пределов интегрирования получаем выражение:

4., Проведение подстановки (11.9) в (11.10) и интегрирование:

5. Выражаем функционал через узловые значения температуры, для чего объединяем выражения (11.7) и (11.11):

Здесь приняты следующие обозначения:

С⁽¹⁾ = (А⁽¹⁾l⁽¹⁾/L⁽¹⁾); С⁽²⁾ = (А⁽²⁾l⁽²⁾/L⁽²⁾)

6. Получение системы алгебраических уравнений. Правильными значениями Т₁, Т₂ и Т₃ являются те, при которых величина функционала c достигает минимума. Приравнивая нулю первую производную функционала (11.12) по Т₁, получаем первое уравнение системы:

Аналогично получаем еще два уравнения:

Запишем полученную систему в матричной форме:

В более общей матричной форме система примет вид:

Матрица C в формуле (11.16) называется «глобальной матрицей жесткости». В контексте задачи переноса тепла –это – «глобальная матрица теплопроводности». Вектор-столбец F называется «глобальным вектором нагрузки». Искомый вектор [T] будем называть вектором решения.

Пример 11.1. Рассчитать температурное поле в круглом стержне с площадью поперечного сечения A=1 см² и длиной L=7, 5 см с теплоизолированными стенками. К левому концу стержня подводится тепловой поток q = 150 Вт/см². Коэффициент теплопроводности материала стержня и коэффициент конвективного теплообмена на правом конце стержня соответственно равны: l=75 Вт/(см · ^ОС), h = 10 Вт/(см² · ^ОС). Температура окружающей среды равна Т_ОС=40 ^ОС.

Решение.

1. Тепло подводится к стержню, поэтому тепловой поток q следует записывать со знаком «минус»: q = - 150 Вт/см².

2. Рассчитываем значение термов, входящих в коэффициенты матриц C и F:

С⁽¹⁾ =(А⁽¹⁾l⁽¹⁾/L⁽¹⁾)=(1·75/3, 75)=20Вт/(см·^ОС),

С⁽²⁾ =(А⁽²⁾l⁽²⁾/L⁽²⁾)=(1·75/3, 75)=20Вт/(см·^ОС),

hA₃=10Вт/(см·^ОС), -qA₁= -(-150)·1 = 150Вт/см,

hA₃T_OC=10·1·40 = 400Вт/см.

3. Окончательная система уравнений примет вид:

4. Решением полученной системы являются следующие узловые значения температуры: Т₁=70 ^оС, Т₂=62, 5 ^оС; Т₃=55 ^оС.