СРАВНЕНИЕ ВЫЧИСЛИТЕЛЬНОЙ СЛОЖНОСТИ МЕТОДОВ ГРАНИЧНЫХ И КОНЕЧНЫХ ЭЛЕМЕНТОВ НА ПРИМЕРЕ ЗАДАЧИ ТЕПЛОВОГО МОДЕЛИРОВАНИЯ СТЕРЖНЯ ИГОЛЬЧАТОГО РАДИАТОРА Текст научной статьи по специальности «Физика»

Электронный журнал «Труды МАИ». Выпуск № 72

www.mai.ru/science/trudy/

УДК. 621.38

Сравнение вычислительной сложности методов граничных и конечных элементов на примере задачи теплового моделирования стержня игольчатого радиатора

Добряков В.А., Енгалычев А.Н.*, Назаров A.B.**

Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское ш., 4, Москва, А-80, ГСП-3, 125993, Россия

*e-mail: aengalychev@ inbox. ru; **e-mail: naza46@post.ru

Аннотация

Известно, что метод граничных элементов (МГЭ) уменьшает размерность исходной задачи на единицу, делая двумерные задачи одномерными, а трехмерные -двумерными. Следует ожидать также пропорционального снижения вычислительной сложности МГЭ по сравнению с методом конечных элементов (МКЭ). В настоящей работе проведена количественная оценка применения обоих методов к задаче теплового моделирования одного и того же технического объекта, которая показывает, что МГЭ позволяет находить более точное решение при меньших затратах машинного времени. Ключевые слова:

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

Введение

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

ханизма теплопроводности. Внедрение в авиационную и информационно-телекоммуникационную технику так называемых теплорассеивающих полимерных композитов (ТРПК) позволило многократно повысить теплопроводность радиаторов. Основной сферой применения ТРПК в настоящее время является изготовление так называемых термоинтерфейсов — изделий, обеспечивающих передачу тепла непосредственно от теплогенерирующего источника в окружающую среду. Типичным представителем термоинтерфейсов являются именно радиаторы охлаждения. В последнее время наблюдается переход от традиционной ребристой к игольчатой конструкции радиаторов, что дает увеличение эффективности охлаждения радиатора на 60-100%. при сохранении установочных габаритов. Термоинтерфейсы постепенно заменяют собой традиционный алюминий за счет использования наполнителей с высокой теплопроводностью (до 200 Вт/(м-К). В результате снижается полетный вес и стоимость устройств охлаждения. В этой связи стержень игольчатого радиатора выбран в настоящей работе в качестве технического объекта для теплового моделирования.

Информация о распределении температурного поля (фазовой переменной) в стержне игольчатого радиатора получается путем вычислительного эксперимента с использованием математических моделей, описываемых дифференциальными уравнениями в частных производных - в данном случае уравнением Лапласа. В современных САПР инженерного анализа решение подобных уравнений обычно находится методом конечных разностей или методом конечных элементов МКЭ. Наряду с ними в практике инженерного анализа все чаще применяется метод граничных элементов (МГЭ) [1, 2, 7]. В отличие от МКЭ, МГЭ позволяет определить значения фазовых переменных во внутренних точках дискретизированной области (©) на основании информации о поведении искомых физических величин и/или их производных на границе (в) области ©.

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

Постоянный граничный элемент аппроксимирует граничный сегмент прямой линией, которая соединяет крайние его точки. Узел помещается в середине прямой линии, а граничная величина принимается постоянной вдоль всего ГЭ и равной значению в узловой точке. Линейный граничный элемент представляется прямой линией, соединяющей две его конечные точки (узлы). Фазовая переменная линейно изменяется между двумя узловыми значениями. Геометрия квадратичного (параболического) граничного элемента представляет собой параболическую дугу с тремя узлами, расположенными по краям и в середине дуги.

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

На рис. 1 показан однородный стержень длиной Ь [см] и сечением Я [см ], ограниченный точками Q и Р, температуру в которых обозначим соответственно Т(О) и Т(Р) [0С]. В точке приложения нагрузки (В) с координатой 4 [см] находится точечный источник тепла интенсивности д=д(В) [Вт/см ]. Точка наблюдения или точка поля (Р) имеет абсциссу х[см]. В одномерном приближении температура Т(х) [0С] во всех внутренних точках удовлетворяет уравнению Лапласа:

д 2Т

дх2

а тепловой поток описывается уравнением:

= 0 (1)

д(х ) = -ЛдТ (2)

дх

При нулевых граничных условиях решение уравнений (1 и 2) в любой точке объекта, кроме х=4, можно записать в виде [1]:

Т(х) = д • х, Я(х) = -Я при 0<х<4 (3)

Ь л Ь

Т(х) = д4 ' д(х) = -д4 при 4< х<Ь (4)

Ь л Ь

При величина T(x) однозначно определяется из (3) и/или (4) , а при переходе точки наблюдения через точку приложения нагрузки, q(x) меняется скачком, равным q(B).

В

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

Пример 1

Построить график изменения температуры T(x) и интенсивности теплового потока q(x) вдоль цилиндрического стержня с теплоизолированными боковыми стенками, если его длина L= 150 см, площадь поперечного сечения S = 1 см и в сечении на расстоянии £ = 50 см от левого конца действует источник тепла интенсивностью q=150 Вт/см . Внешними силами температура в граничных точках поддерживается равной нулю. Коэффициент теплопроводности материала стержня принять равным Л =75 Вт/(см oC).

Решение

На участке стержня 0<x<^, функция T(x) в выражении (3) линейна относительно x, поэтому график T(x) представляет прямую линию. Так как T(Q) = T(0)=0 (по условию задачи), то для ее построения достаточно найти ординату T(x) в точке x=£:

T(x=50) = T(x) = q^^•x = = 150(150-50) •50 -66,7 (oC)

L Л

150 75

На этом же участке интенсивность теплового потока, согласно (3) постоянна и

равна:

q(x) = -q iL-H = = -(-150)

L

(150 - 50) 150

q(x) = 100 (Вт)

Здесь и далее: если тепловой поток подводится к объекту, то он берется со знаком минус.

Вычисления для участка x<L выполняются аналогично. На рис. 2 приведены искомые графики, на которых наглядно видно, что полученное решение является

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

Если в уравнения (3) и (4) подставить q =1, то функции Т(х) и д(х), которые в них фигурируют, могут выступать в качестве функций влияния источника на величины соответственно температуры и теплового потока в точке наблюдения.

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

50 100 150

Рис. 2. График температурного поля внутри стержня

В общем случае построение решения методом МГЭ состоит из следующих пяти шагов:

- получение фундаментальных решений для неограниченной области;

- вывод соотношений, связывающих неизвестные и заданные узловые значения в области © с границей Г;

- получение дополнительных энергетических соотношений;

- решение полученной в пунктах 2 и 3 системы уравнений;

- вычисление неизвестных узловых значений.

Рассмотрим более подробно шаги реализации каждого этапа МГЭ на примере вычисления температурного поля однородного стержня.

Получение фундаментальных решений для неограниченной области

Непосредственно в неограниченных системах можно определить не саму функцию Т(х), а лишь ее отклонение относительно некоторого фиксированного значения. Примем в качестве такого значения величину температуры Т(Ь)=0 в сечениях ±Ь

стержня. По аналогии с выражением (3), тепловой поток q(R) [Вт/см ], действующий в сечении R[см], поднимет температуру Т(Р') в сечении Р' стержня на величину:

Т (Р') = ШУ(2Ш( I - I г |) = ^ЮО(Р\ Ю (5)

где 0(Р',Ю) - выражение, в котором сгруппированы все члены равенства кроме q(R), - представляет, по сути, оператор, зависящий от упорядоченной пары (Р', Ю), причем такой что, будучи умноженный на q(R) , даст значение температуры в точке Р'. Переходя далее к абсолютным координатам, перепишем (5) в виде:

Т(х) = q(£)G(x& (6)

Рис. 3. Неограниченная одномерная область

Обозначим q(r) -тепловой поток в точке на расстоянии г от источника. Учитывая выражение (2) и дифференцируя (5):

q(r) =[^/2]^) = ^) F(P',R) (7)

где выражение F(P', R), в котором сгруппированы все члены равенства кроме q(R), -представляет оператор, зависящий от упорядоченной пары (Р', R), причем такой, что, будучи умноженный на q(R) , даст значение интенсивности потока в сечении Р'. Переходя далее к абсолютным координатам, перепишем (7) в виде:

ф) = q(Z)F(x& (8)

Функция sgn(r) принимает значение +1 (-1) в сечениях справа (слева) от источника теплового потока q (рис.2). В самой же точке приложения потока функция sgn(r) не определена, что позволяет автоматически учесть скачкообразное изменение q(г) при г —^ 0.

/П I Ф

а)

б)

Рис.4 Стержень (а), помещенный в неограниченную одномерную область ±оо (б)

Реальный объект в неограниченной системе

Поместим реальный технический объект (рис.4а) в неограниченную систему (рис.4б). На рис.4 реальные значения параметров в отличие от тех же параметров в неограниченной системе помечены звездочками. Поскольку границами одномерной области (2Р) являются просто ее две точки (2 и Р), достаточно ввести в этих точках по одному граничному элементу. Далее в обоих граничных элементах введем фиктивные источники неизвестной пока интенсивности q(Q) и q(P). Влияние введенных источников на каждую внутреннюю точку Р' может быть выражено в виде (6 и 8). Каждый из источников - q(P), q(Q) и q(B), показанных на рис.5, оказывает влияние на температуру сечения в точке х. Принцип суперпозиции позволяет оценить суммарное их действие в виде:

Т(х) = q(Q)G(x,Q) + q(P)G(x,L) + ^В)С(х£) (9)

ф) = тЩх^) + q(P)F(x,L) + (10)

Полученные равенства справедливы для любой точки наблюдения в пределах отрезка QP. Выберем в качестве таковой точку левее точки Р, т.е. точку с абсциссой (Ь- при 4^0. По формулам (5) и (9) находим температуры в точках P и Q: Ш) = ^ ШУ(2Х)] [^ - 0) + [q(P)/(2X)](L - L) + [ф)/(2Х)] ^ -После некоторых упрощений получим:

Ш) = [1/(2Х)] [(L ^) + 0 + (L - 4) ^В)] (11)

Аналогично находим: Т^) = [1/(2Х)](0 ■q(Q) + L q(P) + В)]. По формулам (7 и 10) находим интенсивность потока в точках P и Q:

^) = (1/(2) Ш) - ^) - ^В)]

^) = (1/(2) [т) - ^) + ^В)] (12)

Потребуем теперь, чтобы граничные условия в точках P и Q фиктивной системы в точности совпадали с условиями реальной задачи. С этой целью в уравнениях (11) и (12) сделаем замену: T(Q)=T*(Q) и q(P)= q*(P), после чего (11) и (12) примут вид:

Т*{0)^2-к= + 0 ) + ^ - &(В)] (13)

2q*(P) = ^) - ^) + ф) (14)

Решая систему (13) и (14), получим обе неизвестные величины: q(Q) и д(Р), которые после подстановки в выражения (9) и (10) позволят вычислить лишь относительные значения искомых величин (см. «Получение фундаментальных решений»). Для обеспечения единственности решения потребуем, чтобы обращался в нуль полный поток через бесконечно удаленные границы, что эквивалентно обращению в 0 суммы потоков q(P), q(Q) и q(B). С этой целью в выражении (5) заменим источник q(R) на источник ^(Я). Замена вызовет коррекцию температуры на заранее неизвестную константу С, после чего выражение (5) примет вид:

Т(Р') = [^(Я)/(2Х)] • (Ь - | г | ) + С и вытекающие из него уравнения (13 и 14) так же изменятся и дополнятся условием равенства нулю суммы всех потоков, включая модифицированные: 2ЛТ*(Q) = Ь q'(Q) + 0 • q'(P) + (Ь - 4) q(B)] + 2ЛС 2^*(Р) = q(Q) - q(P) + q(B)

Ш) - ^(Р) + Ф)=0

В матричной форме данная система уравнений примет вид:

(15)

Ь 0 2Л V ^ ) ' Т * ь-4

1 -1 0 < ^ = <2q * ^) > — 1 ' q(B)

1 1 0 С 0 1

Решая систему(15), получим выражения для вычисления q'(Q), q'(P) и С:

= д^)-д(В), q'{p) = -q{p) С = Т * (а)-* ^ )+&{в)]

(16)

Подставив эти равенства в уравнения (9, 10), можно вычислить Т(х) и q(x) в любой точке P'(x) внутри интервала (PQ).

Так, при 0<х<4из равенств (5 и 9) имеем:

Т(х) [(Ь- х) (-g(B)) + (Ь- 4 + хЫВ) -q* х]+ С 2Л

х х

откуда с учетом (16): (х) = T*(Q) - —q*(P) + — q(B) и

Л Л

q (х) = 0,5^^ - q(B) + q*(P) - ц^)] = q*(P) ^) Аналогично на интервале 4 <х< Ь имеем:

Т(х) = Т*^) - (х/X) + ($/X) ■ д(В) ф) = q*(P)

(18)

Пример 2

Рассчитать температурное поле в круглом стержне с теплоизолированными стенками, показанном на рис. 5. Площадь поперечного сечения стержня А=1 см , его длина L=7,5 см. Температура на левом конце стрежня поддерживается постоянной и равной 70оС. Тепловой поток, отводящий тепло от правого конца стержня, равен 150 Вт/(см2),

как это показано на рис. 6. Коэффициент теплопроводности материала стержня Х=75 Вт/(см оС).

Решение

Из формулы (18) получаем q(x)=q*(P)=150 Вт/(см ). Так как на интервале (PQ) внутри стержня источников тепла нет, то q(B)= 0, и формула (17) примет вид:

Т*(х) = T*(Q) - (х/X) ■ q*(P) =

= 70 - (х / 75) -150 = 70 -2х

Таким образом, температура линейно убывает и на правом конце стержня составит:

Т(х=7,5) = 70 - 2 7,5 = 55 оС

Рис. 5. Стержень с теплоизолированными боковыми стенками

Рис. 6. Конечно-элементная модель стержня

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

Получим расчетные формулы, которые позволяют решить данную задачу (пример 2) методом конечных элементов. Для этого необходимо минимизировать следующий функционал [3]:

Л

дТ (х) дх

2Л

(У +1 ^Т (х)) (18,

(19)

где: 8 -площадь поверхности стержня, по которой идет теплообмен с окружающей средой;

Т (х) -температура в точке стержня с абсциссой х;

У - объем стержня.

Если представить исходный стержень одним конечным элементом, то невозможно будет подтвердить линейный характер изменения температуры вдоль стержня, поэтому нужно взять 2 конечных симплекс-элемента для моделирования, как показано на рис. 6, и получим решающую систему уравнений этапами:

1. Записываем интерполяционный полином, выражающий температуру в точке с абсциссой х (для 1-го конечного элемента) через искомые узловые значения температур Т и Т.:

Т (х) = МТ + М.Т. (20)

где: и N. - функции формы, равные соответственно: ы, = ((х;. -х)/Ь) и N. =((х - X, ) / Ь) .

2. Интерполяционные полиномы для 1-го и 2-го конечного элемента примут

вид:

Т[1] (х) = М1[1]Т1 + М21]Т2 и Т[2] (х) = м22]Т2 + ы32]Т3, где: ы[1] =((х2 -х)/ь) , м1 =((х - х1)/ь) и м22] =((х3 - х)/ь), м32] =((х- х2)/ь) . 3. Вычисляем частные производные, входящие в функционал (19):

(21)

дТ [1]( х)' дх

Т1 + Т 2 Ь

дТ [2]( х)' дх

Т 2 + Т 3 Ь

(22)

4. Разбиваем объемный интеграл в выражении (19) на два линейных интеграла (по числу конечных элементов), учитывая постоянство сечения стержня по длине ((У=А(х):

2

8

Г а " дТ (х)" 2 Л ¿V = ААхг " дТ (х)"

2 V _ дх _ У 2 1 2 Х1 _ дх _

X 3

л АА Г Эх + — |

2

X 2

дТ (х)' дx

дх

(23)

5. Выполняем подстановку (22) в (23) и интегрируем:

А

дТ (x) дx

2\

dV = А- (- Тх + Т2 )2 + А- (- Т2 + Т3 )2

(24)

6. Учитывая, что:

- теплообмен стержня идет только по поверхностям S1 и S2 (по условиям задачи);

- во всех точках поверхности S1(S2) искомая температура и величина теплового потока в установившемся режиме постоянна;

- тепловой поток подводится к поверхности S1 и отводится от поверхности S2; определяем температуры:

- в узле 1: q = -д , Т (x = X1) = Т1,

- в узле 2: q = 0, Т(х = X2) = 72,

- в узле 3: q = + q, т(х = X3) = Т3.

7. Вычисляем поверхностный интеграл в выражении (19). С учетом сказанного в п.6:

| {дТ(х)) dS = -qT1 А + 0 + qT3А.

(25)

АА

8. Введём обозначение: С = —, а суммируя выражения (24) и (25), получаем

исходный функционал для минимизации:

Си 2

Х = С(Т]2 - 2Т]Т2 + 272 - 2Т2Т3 + Т32)2 - дАТ] + дАТ3.

(26)

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

2

I

V

S

26) на максимум-минимум: приравняем нулю производную от х по каждой узловой переменной, что приводит нас к системе:

СТ1 - СТ2 + 0 = дЛ

- СТ1 + 2СТ2 - СТ3 = 0 0 - СТ2 + СТ3 =-дЛ

10.Решаем полученную систему уравнений: -поскольку С=20 Вт/ °С; дЛ=150 Вт и (по условию) Т1 = 70 °С, то из 1-го уравнения получаем: 20x70 - 20хТ2 =150, откуда Т2 =62,5 0С.

2дЛ

-выполняя вычитание 2-го уравнения из 1-го, имеем: Т1 -Т3 ; или

70 - Т3 = 300; отсюда Т3 = 55 0С.

3 20 3

Полученная одинаковая разность температур между узлами 1 и 2, и соответственно, между узлами 2 и 3 (7,5 0С) подтверждает линейный характер изменения температурного поля вдоль стержня с изолированными стенками в установившемся режиме, что совпадает с решением этой задачи методом граничных элементов.

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

Приведенное выше решение задач подтверждает выводы, сделанные в работе [4] о том, что «... как только получена необходимая информация о границе объекта, так сразу могут быть вычислены значения переменных, описывающих решение, в любых последовательно выбираемых внутренних точках, и что решение полностью непрерывно всюду внутри исследуемого объекта». Данный вывод согласуется с расчетным выражением, полученным выше при решении второго примера методом граничных элементов, а также с результатами исследований Лаша и Бенерджи [1, 5], которые показывают, что время решения одних и тех же трехмерных задач методом

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

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

Библиографический список

1 .Бенерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках // Пер. с англ. М.: Мир, 1984. 494 с.

2.Кацикаделис Дж. Т. Граничные элементы. Теория и приложения. Пер. с англ. - М.: Изд-во Ассоциации строительных вузов. 2007. - 348 с. ил.

З.Зенкевич С, Старфилд А. Метод конечных элементов в технике/ Пер. с англ. М.: Мир, 1987. 328 с.

4.http://masters.donntu.edu.ua/2000/fvti/sugonyak/du/mke&mge.htm(10/06/13)

5.Lachat, J.C. Further developments of the boundary integral techniques for elasto-statics: Ph.D. thes / J.C. Lachat. Southampton Univ., 1975.

6.Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979.

7.Дембицкий Н.Л., Назаров А.В. Модели и методы в задачах автоматизированного конструирования радиоприемных устройств. - М.: МАИ, 2011.