 |
Постановка задач расчета теплового процесса на дискретной модели(пример в 24)
|
|
|
|
В электронно-вычислительной аппаратуре могут происходить следующие процессы передачи тепла: конвекция, кондукция и лучеиспускание. Разностный метод не применим для расчета передачи тепла конвекцией и лучеиспусканием. Поэтому ограничимся рассмотрением конструкций, элементы которых помещены либо в твердый наполнитель, например, в пенопласт, компаунд или резину, либо в желеобразный наполнитель с малым коэффициентом текучести. То есть будем рассматривать конструкции, в которых происходит только передача тепла теплопроводностью (кондукция). Предположим, что блок ЭВА имеет прямоугольную форму, внутри которого находятся источники тепла – радиоэлементы, через которые протекает электрический ток. Блок залит наполнителем с коэффициентом теплопроводности К и удельной теплоемкостью С. Разобьем мысленно конструкцию на части прямоугольной формы, каждую из которых назовем элементом.
| . 
|
| Рис.2.4.1. | Рис.2.4.2. |
Для более высокой точности расчета выберем элементы одинаковых размеров, причем сами размеры элементов примем минимально возможными. В центре элемента выделим особую точку – узел сетки. Далее, попытаемся определить температуру в каждом узле сетки в каждый момент времени. Для простоты будем считать, что блок однороден, то есть входящие в него материалы имеют одинаковую теплоемкость и коэффициент теплопроводности. Температуру, определяемую в узле сетки с координатами x, y, z, в момент времени t, обозначим как QtX,Y,Z, а в следующий момент времени как Qt+1X,Y,Z. Размеры блока, координаты и мощность тепловыделяющих радиоэлементов будем считать заданными. Кроме того, для решения задачи должны быть заданы начальные и граничные условия.
|
|
|
В начальных условиях задачи необходимо указать температуру во всех узлах сетки блока в начальный момент времени. Обычно при рассмотрении переходных процессов за начальный момент времени выбирается момент включения электрических цепей под нагрузку. До этого момента температура во всех узлах считается одинаковой и равной наружной температуре, например, комнатной (20ОС или 293ОК).
В граничных узлах блока могут быть заданы различные граничные условия. Когда на границе задается значение самой функции, то есть температура – это граничное условие 1-го рода и решение получается наиболее простым. Однако, к сожалению, только при грубом упрощении нестационарной задачи (то есть задачи с изменением температуры во времени) можно считать температуру на поверхности бока заданной, например, равной наружной температуре. Наиболее близкими к реальным условиям являются граничные условия 2-го рода, когда задаются плотности теплового потока по всей наружной поверхности блока. Далее мы будем решать задачу с граничными условиями 2-го рода.
Метод конечных элементов
Метод конечных элементов представляет собой сеточный метод, который служит для решения задач микроуровня. Модель ТО задаётся здесь системой ДУ в частных производных с заданными краевыми условиями [8]. В отличие от метода конечных разностей метод конечных элементов основан на разбиении объекта на элементы конечных размеров. Поведение каждого элемента приближенно воспроизводит поведение той части объекта, которое он аппроксимирует. Точность решения зависит от размеров и типа конечных элементов.
В узлах конечно-элементной сетки на решение накладывается условие полной непрерывности искомой функции между элементами. На основе МКЭ построено большинство систем автоматизации инженерного анализа ТО, что говорит о его большой популярности.
|
|
|
Типы конечных элементов. При решении задач методом конечных элементов используются одномерные, двухмерные и трехмерные конечные элементы. В одномерной задаче конечные элементы – это отрезки прямой линии. Конечные элементы в двумерной задаче – это треугольные или четырёхугольные элементы. Последние элементы обычно выбираются со сторонами параллельными осям координат. Предпочтение обычно отдается трехмерным конечным элементам, что ведет к снижению размерности решающей системы уравнений. В качестве аппроксимирующей функции, действующей внутри элемента, обычно применяется полином, порядок которого и определяет тип конечного элемента. На практике используются три типа элементов: симплекс-элемент, комплекс-элемент и мультиплекс-элемент [17].
Симплекс - элементам соответствуют полиномы, содержащие константу и линейные члены. Одномерный симплекс-элемент аппроксимирует функцию линией: j(х) = a1 +a2 x. Полином, аппроксимирующий поведение функции в пределах треугольного симплекс-элемента, дается уравнением плоскости: j(х) = a1 +a2 x+ a2 y.
Комплекс–элементам соответствуют полиномиальные функции, содержащие константу, линейные члены, члены второго и более высоких порядков. Форма комплекс – элементов может быть такой же как у симплекс - элементов, но с дополнительными граничными или внутренними узлами. Интерполяционный полином для 2-мерного треугольного комплекс – элемента с шестью узлами имеет вид:
.
Для мультиплекс – элементов также используются полиномы, содержащие члены высокого порядка, Однако, для обеспечения непрерывности при переходе от одного мультиплекс – элемента к другому границы элементов должны быть параллельны осям координат.
Вывод уравнений, описывающих конечный элемент, продемонстрируем на числовом примере решения задачи определения изгиба консольной балки (рис.3.3).
рис.3.3
Условие задачи. Консоль длиной L0 =1,5 (м) жестко закреплена в точке х=0 и на свободном конце подвержена действию сосредоточенной силы величиной F =4500 (Н). Точки А, В, C, D делят балку на 5 равных частей. Конечно-элементная модель консоли приведена на рис.3.4. Требуется определить прогибы Y(x) консоли в узлах модели, если ее изгибная жесткость равна: EJ = 8,5×108 (Н×см2).
|
|
|
Рис.3.4.
Введем обозначение для кривизны балки: 
. По определению [21], изгибающий момент 
в произвольном поперечном сечении консольной балки численно равен алгебраической сумме моментов внешних сил, действующих по одну сторону от сечения, относительно той точки оси балки, через которую проходит сечение. Если момент, изгибает ось балки выпуклостью вверх, как в нашем случае, то ему приписывается знак минус. Числовые значения для кривизны балки в узловых точках приведены в табл.3.
Решение. Прогиб балки y(x) описывает ДУ вида:
.
Исследования данного ДУ методами вариационного анализа показывают, что в интересующем нас напряженном состоянии балки достигает экстремального значения функционал:
Таблица 3
Числовые значения кривизны балки в узловых точках
| Узел |  [см-1]
| Узел | [см-1]
| Узел | [см-1]
|
| –0,000794 | –0,000476 | –0,000159 | |||
| –0,000635 | –0,000318 | –0,000000 |
. (3.6)
Системы инженерного анализа конструкций, реализующие метод конечных элементов, не работают непосредственно с исходным ДУ. Целью их расчета является именно поиск экстремума функционала (3.6). Математически строго доказано, что найденная таким образом точка экстремума и является искомым решением соответствующему ему ДУ, в данном случае – уравнения (3.1). Однако, и с функционалом (3.6) системы инженерного анализа также непосредственно не работают. С чем же они имеют дело?
В основе большинства систем инженерного анализа лежит принцип ассемблирования локальных матриц жесткости всех конечных элементов, на которые предварительно разбит ТО, – в результирующую глобальную матрицу жесткости всего ТО. Физический смыслглобальной матрицы жесткости состоит в том, что ее элементы представляют коэффициенты при N неизвестных в полученной после указанного пересчета системе из N решающих линейных уравнений, где N – общее число узлов в конечно-элементной модели ТО. Указанные коэффициенты обычно представляют физические параметры моделируемого технического объекта.
Придерживаясь методики, изложенной в работе [8], покажем далее, каким образом можно на основе функционала (3.6) получить матрицы конечных линейных элементов. Обозначим общее количество конечных элементов, на которое разбивается ТО, символом Е (в текущем примере Е =5). Учитывая аддитивный характер функционала (3.6), представим его суммой Е функционалов вида:
|
|
|
, где: (3.7)
- е - номер конечного элемента;
- 
- элементарный функционал, представляющий собой сумму интегралов для произвольного (е -го) конечного элемента.
Для минимизации функционала (3.7) необходимо:
- задаться типом конечного элемента и получить полиномиальную функцию, описывающую поведение искомой функции в пределах конечного элемента;
- выразив полиномиальную функцию через узловые значения искомых перемещений Y, получить интерполяционный полином;
- вычислить производные от функционала (3.7) по всем узловым значениям: 
и, приравняв их к нулю, получить искомые матрицы элементов.
Отметим, что величины 
прогибов балки в узловых точках постоянны только в установившемся (напряженном) состоянии. В функционале (3.7) эти величины зависят от 
, поскольку сам функционал (3.7) моделирует переходный процесс нагружения консоли, в процессе которого величины меняются от нуля (в начале нагружения) до искомых постоянных величин (в конце нагружения).
Получение интерполяционного полинома проиллюстрируем, выбрав в качестве конечных элементов одномерные симплекс-элементы, один из которых показан на рис.3.5. Дополнительно введем несколько принятых в методе конечных элементов терминов.
Рис.3.5.
Одномерный симплекс-элемент представляет собой прямолинейный отрезок длины L с двумя узлами – по одному на каждом конце отрезка. Узлы обозначаются индексами i и j, значения функции в узлах – через 
и 
соответственно. Начало системы координат может не совпадать с началом координат. Полиномиальная функция Y для скалярной величины (в данном случае – линейного перемещения узла по оси OY) такова: Y = a1 +a2 x
Коэффициенты a1 и a2 определяются с помощью условий в узловых точках: j = 
при x = Xi и j = 
при x = Xj. Эти узловые условия приводят к системе двух уравнений:
= a1 +a2 Xi,
= a1 +a2 ,
решение которой дает: 
, 
. Подставляя найденные значения a1 и a2 в полиномиальную функцию, получим: 
, которое может быть переписано в виде: 
. Линейные функции от х в полученном таким образом интерполяционном полиноме называются функциями формы или интерполяционными функциями. Каждая функция формы должна быть снабжена нижним индексом для обозначения узла, к которому она относится. Произвольную функцию формы будем обозначать через N b. Запишем отдельно выражения для функций формы:
и 
, (3.8)
|
|
|
И перепишем интерполяционный полином в матричной форме:
. (3.9)
Функция Ni = 1 в узле с номером i и равна нулю в j -м узле. Аналогично функция Nj = 1 в узле с номером j и равна нулю в i- м узле. Эти значения характерны для функций формы. Они равны 1 в одном определенном узле и обращаются в 0 в остальных узлах.
Представление функционала 
в матричном виде. Введение в рассмотрение интерполяционного полинома позволяет существенно сократить последующие выкладки, представив функционал в матричном виде. С этой целью, учитывая соотношение (3.9), запишем интерполяционную формулу для произвольного симплекс – элемента в общем виде:
Y (е) = Ni(е)Yi + Nj(е) Yj = [ N(е) ] { Y }.(3.10)
Частная производная от функции Y(е) (3.10) по координате х, входящей в функционал (3.6) равна:
. (3.11)
Обозначив:
и 
, (3.12)
запишем (3.11) в матричной форме:
. (3.13)
Это позволяет получить выражение для функционала c(е) в матричной форме. Согласно выражениям (3.9), (3.12) и (3.13) имеем:
(3.14)
Перед вычислением производных от полученного функционала по всем узловым значениям покажем, что в матричном виде:
.
Левая часть приведенного тождества представляет собой искомую частную производную от квадрата частной производной Y(е) по Yi, представленную в матричной форме, которая по определению производной от сложной функции и с учетом (3.11) равна:
,
что и требовалось показать.
Вычисление производных от функционала (3.14) по всем узловым значениям перемещений. Для первого конечного элемента с учетом выражения (3.11) имеем:
.
Вектор { Y } не зависит от переменной интегрирования, поэтому, вынося его за скобки и вычисляя производную по второй узловой переменной, приходим к системе:
.
В данной системе длинными квадратными скобками выделены элементы, представляющие собой транспонированные матрицы [В(е)]T и [N(e)]T. Выделение позволяет кратко записать вклад отдельного конечного элемента в общую сумму (3.7) в матричном виде:
.
Данное выражение существенно сокращает запись решений задач методом конечных элементов и фактически заменяет следующую эквивалентную ей пару уравнений:
Введем обозначения:
, (3.15)
. (3.16)
Матрица 
в выражении (3.15) называется локальной матрицей жесткости конечного элемента. Матрица 
в выражении (3.16) называется вектором нагрузки конечного элемента. Матрицы и называются далее матрицами е-го конечного элемента.
Запишем в матричной форме соотношение, представляющее вклад отдельного конечного элемента в общую сумму (3.7):
. (3.17)
Приравнивая выражение (3.17) нулю, запишем локальную систему уравнений:
. (3.18)
Решающую систему уравнений получим, дифференцируя по узловым переменным функционал (3.6) и приравнивая нулю полученные производные. Учитывая выражения (3.17) и (3.18) имеем:
Приравнивая полученное выражение нулю, получим окончательную решающую систему линейных уравнений, которая известна под названием глобальной системы уравнений:
. (3.19)
В выражении (3.19) квадратная матрица 
является суммой:
(3.20)
и называется глобальной матрицей жесткости ТО, матрица-столбец: 
– глобальным вектором нагрузки, а матрица-столбец 
– вектором решения.
Таким образом, получены основные расчетные соотношения, достаточные для числового решения текущей задачи, которое включает этапы:
1. Задаем начальные условия. Учитывая факт закрепления левого конца балки, имеем: Y(0)=0.
2. Упрощаем матрицы элементов. Учитывая, что площадь поперечного сечения консоли (S) постоянна, переходим от интегралов по объему к интегралам по длине консоли, для чего: обозначим 
– длину конечного элемента, и после замены в формулах (70 и 71): 
, запишем матрицы элементов в упрощенном виде:
(3.21)
Кривизна балки в (3.21) осталась под знаком интеграла, так как она зависит от 
в пределах 
-го конечного элемента и определяется интерполяционным полиномом:
. (3.22)
3. Вычисляем матрицы функций формы и матрицы градиентов. Поместив узел 
в начало координат, вычислим функции формы для первого элемента по формуле (3.8):
и 
.
Таким образом, матрица функций формы первого конечного элемента примет вид:
.
Аналогично вычисляем матрицы функций формы остальных конечных элементов:
, 
,
.
Вычисляем матрицы градиентов всех конечных элементов по формулам (3.12): 
.
4. Вычисляем локальные матрицы жесткости. Учитывая, что элементы матриц градиентов не зависят от координаты x, локальная матрица жесткости для первого элемента примет вид:
.
Произведение матриц 
дает следующий результат:
Следовательно: 
и окончательно матрица жесткости 1-го конечного элемента примет вид:
Аналогично получаем матрицы жесткости для остальных конечных элементов:
Здесь в колонке справа стрелками указаны номера узлов конечно-элементной модели ТО, необходимые для правильного ассемблирования его глобальной матрицы жесткости, причем указанные номера строк совпадают с номерами столбцов соответствующей локальной матрицы. Так, отмеченному крестом элементу (–1) матрицы 
, соответствуют узлы 4 (по вертикали) и 3 (по горизонтали). Его вклад в глобальную матрицу жесткости (3.20) иллюстрируется в следующем пункте.
5. Вычисляем глобальную матрицу жесткости, которая ассемблируется по правилу прямой жесткости, которое гласит: элемент глобальной матрицы, расположенный на пересечении i-й строки и j-го столбца, равен сумме элементов всех локальных матриц, расположенных на пересечении строки и столбца, соответствующих узлам i и j. В нашем примере после сложения полученных локальных матриц результирующая глобальная матрица примет вид:
Здесь крестиком отмечен элемент (–1), о котором шла речь в предыдущем пункте.
6. Вычисляем локальные матрицы нагрузки. Учитывая формулы (3.21,3.22), имеем:
.
Подставляя сюда из пункта 3 выражения для функций формы: 
и 
и, вычисляя интегралы, получим:
.
Подставляя числовые значения 
и 
из таблицы 3:
[см],
окончательно получаем локальную матрицу нагрузки для первого конечного элемента:
Локальные матрицы нагрузки остальных конечных элементов примут вид:
7. Вычисляем глобальную матрицу нагрузки, складывая полученные локальные матрицы нагрузки по методу прямой жесткости:
8. Формируем решающую систему, по формуле (76):
Так как Y1=0, то из первого уравнения имеем: Y2 = –0,33345 [см].
Из второго уравнения: 
Аналогично далее: 
.
Сравнение полученных результатов с данными теоретического расчета и с данными, полученными в САПР AnSYS, приведено в табл.4. Расчеты совпадают с точностью до сотых долей миллиметра.
Таблица 4
Сравнение результатов расчета разными методами
| Node | UY(AnSys) | МКЭ | Теория |
| 0.0000 | 0.0000 | 0.0000 | |
| –0,33534E–2 | –0,33345Е–2 | –0,33350Е–2 | |
| –0,12456E–2 | –0,12384Е–2 | –0,12388Е–2 | |
| –0,25869E–2 | –0.25719Е–2 | –0.25729Е–2 | |
| –0,42157E–2 | –0,41929Е–2 | –0,41929Е–2 | |
| –0,59883E–2 | –0,59550Е–2 | –0,59559Е–2 |
|
|
|
