VI Международная студенческая научная конференция Студенческий научный форум - 2014

Характерной особенностью современной науки и высокотехнологического производства является широкое использование математического моделирования, средств вычислительной техники и специализированного программного обеспечения для решения сложных научных и производственных задач. Поэтому базовый учебный план подготовки бакалавров и магистров по направлению 220700 Автоматизация технологических процессов и производств включает в себя изучение таких дисциплин как «Программно-математический инструментарий», «Математическое моделирование», «Основы компьютерного моделирования технологических процессов», «Моделирование систем и процессов».

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

В предложенных лабораторных работах излагается методика расчета тепловых полей в проводнике с помощью модуля PDE-Toolbox пакета MATLAB.

1. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ. КРАТКИЕ ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ

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

Tx=i=0nTiNix . (1)

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

1. дискретизировать область;

2. выбрать базисные функции;

3. перейти к дискретной постановке задачи;

4. вычислить коэффициенты в системе линейных алгебраических уравнений;

5. решить систему уравнений;

6. проверить решение на адекватность.

Рассмотрим подробнее каждый этап применения МКЭ. В процессе дискретизации область решения разбивается на множество элементов. Это могут быть как треугольники, так и тетраэдры (для трехмерного случая). На узлах сетки будет определяться искомое решение, а форма обрабатываемой области может быть любой. Здесь, также, проявляется одно из основных достоинств метода конечных элементов – можно увеличить частоту сетки там, где нужно более точное решение (рис. 1б), или наоборот, уменьшить число элементов там, где столь высокая точность не нужна (рис. 1а).

Далее используется метод Галёркина. Для его применения нужно подобрать базисные функции, которые будут удовлетворять краевым условиям и в пределах бесконечного количества элементов базиса образуют полную систему. То есть, если ω=x1,x2,…, xn– значения функции, то должно выполняться следующее условие (2)

Nix= 1, x=xi 0, x∈ω{xi} . (2)

Это значит, что подставляя какое-либо значение функции в узел сетки, мы получим значение температуры Ti, при условии, что Nix в этом узле равно единице. Во всех остальных узлах сетки базисная функция должна равняться нулю.

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

В общем случае можно найти решение краевой задачи для уравнения вида Lφx=0. Здесь оператор L[ ] может содержать частные или полные производные искомой функции. Решение представляется в виде разложения по базису (1), затем необходимо подставить приближенное решение в исходное дифференциальное уравнение и вычислить погрешность измерений

Li=0nTiNi(x)=Mx.

Для неоднородного уравнения Lu=f(x) погрешность (невязка) будет определяться из выражения Mx=Lu-f(x).

После этого выдвигается условие ортогональности невязки к базисным функциям, то есть

abMxφixdx=0.

После этого умножаем на каждую из базисных функций наше уравнение

ΩLi=0nTiNi(x)*Njxdx=Ωfx*Njxdx, j=1,N

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

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

Долгое время широкому распространению МКЭ мешало отсутствие алгоритмов автоматического разбиения области на «почти равносторонние» треугольники (погрешность, в зависимости от вариации метода, обратно пропорциональна синусу или самого острого, или самого тупого угла в разбиении). Впрочем, эту задачу удалось успешно решить (алгоритмы основаны на триангуляции Делоне), что дало возможность создавать полностью автоматические САПР, основанные на методе конечных элементов.

Лабораторная работа №1

ИЗУЧЕНИЕ ВОЗМОЖНОСТЕЙ ПАКЕТА MATLAB ДЛЯ РЕШЕНИЯ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ

Цель работы: Изучить возможности встроенного в MATLAB модуля PDE-Toolbox. Решить простейшую стационарную задачу теплопроводности (распределение тепла на теплоизолированной пластине с точечным источником тепла).

1. Пояснения к работе

В программном пакете MATLAB предлагается встроенный инструмент PDETOOLBOX(Partial Differential Equation Toolbox), работающий с уравнениями эллиптического/гиперболического типа и определяющий их решение с помощью метода конечных элементов. Геометрия описываемой области в данном случае может быть только двухмерной.

Запустим MATLAB и введем в рабочее поле команду pdeinit, затем Enter. Откроется рабочее окно инструмента PDE Toolbox (рис. 1).

Рис. 1.− Рабочее окно инструмента PDEToolbox

Решение задачи теплопроводности начинается с описания геометрии исследуемой области. Для этого предусмотрен специальный инструментарий. Сперва зададим границы рабочего поля axeslimits. Для этого перейдем в меню Options, далее Axes Limits… (рис.2).

Рис. 2. – Меню AxesLimits

Зададим нужные нам значения x∈-5;5, y∈(-1;1). Затем нужно описать геометрию исследуемой области.

Для этого выберем в меню Draw, Ellipse/circle (centered) и левой кнопкой мыши построим эллипс произвольного размера и центром в начале координат (рис. 3).

Рис. 3.−Построение исследуемой области

Теперь построим на данном эллипсе окружность, с меньшим радиусом. Она будет играть роль точечного источника тепла. Программа автоматически присваивает ей обозначение С1 (рис. 4).

Рис. 4.− Исследуемая область с точечным источником тепла

На этом описание геометрии исследуемой области завершается. Следующий этап – задание краевых условий. Перейдем в пункт меню Boundary, далее Boundary Mode. Так как пластина теплоизолированная, то на ее краях будет выполняться условие Неймана. Для выделения всей границы эллипса E1 необходимо зажать клавишу Shift и поочередно зажать на каждую часть границы (красная линия) левой кнопкой мыши (рис.5). Затем двойным щелчком открыть окно настройки краевых условий.

Рис. 5.− Выделение границы, для задания краевых условий

В случае термоизоляции пластины положим значения Heat flux (тепловой поток) и Heat transfer coefficient (коэффициент конвективной теплопередачи) равными нулю (рис.6).

Рис. 6.− Окно настройки краевых условий

Перейдем к следующему этапу решения задачи − спецификации коэффициентов дифференциального уравнения. Для этого перейдем в пункт меню PDE, далее PDE Mode. Отметим также галочкой Show subdomain labels (показать метки подобластей) (рис. 7).

Рис. 7.− Спецификация коэффициентов дифференциального уравнения

Зададим коэффициенты дифференциального уравнения для зоны 1. Известны коэффициенты теплопроводности каждого из материалов, мощность внутренних источников тепла, а также внешняя температура (рис.8а). Аналогично − для зоны 2. (рис.8б).

Рис. 8а – Спецификация коэффициентов для зоны 1

Рис. 8б – Спецификация коэффициентов для зоны 2

После спецификации коэффициентов диф. уравнения строим вычислительную сетку (см. Триангуляция Делоне). Количество узлов и ее частоту можно регулировать. Переходим в меню Mesh, далее Mesh Mode (рис. 9).

Рис. 9.− Меню Mesh Mode

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

После этого можно искать решение дифференциального уравнения. Переходим в меню Solve, далее Solve PDE. Решение будет представлено в следующем виде: (рис. 10).

Рис. 10.− Искомое распределение теплового поля

Теперь можно настроить визуализацию решения. Перейдем в меню Plot, далее Parameters… Поставим галочки на против Contour (контур) и Arrows (температурный градиент) (рис. 11).

Рис. 11.− Окно настройки визуализации решения

Тогда вид решения изменится следующим образом: (рис. 12).

Рис. 12.− Распределение теплового поля

Из рисунка можно сделать выводы о распределении тепла по поверхности пластины.

2. Содержание отчета

1. Наименование работы, ее цель и краткие теоретические сведения.

2. Пошаговое решение поставленной задачи, с необходимыми пояснениями (геометрия области и мощность тепловых источников задается преподавателем).

3. Выводы по работе.

3.Контрольные вопросы

1. Дифференциальные уравнения в частных производных. Общие сведения.

2. Метод конечных элементов, применение его в различных областях науки и техники.

3. Метод Галёркина и подбор базисной функции.

4. Триангуляция Делоне.

Лабораторная работа №2

СТАЦИОНАРНАЯ ЗАДАЧА РАСПРЕДЕЛЕНИЯ ТЕПЛА В ПРОВОДНИКЕ

Цель работы :Углубить знания о работе встроенного в MATLAB модуля PDE-Toolbox и методе конечных элементов. Решить задачу распределения тепла в проводнике.

1. Пояснения к работе

Рассмотрим поэтапно стационарную задачу распределения тепла в медном кабеле. Пользовательский интерфейс представлен на рис.1

Рис. 1.− Пользовательский интерфейс

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

Зайдем в пункт меню Options (настройки) и далее Axes Limits. Зададим там нужные границы рабочего поля (рис. 2). Выберем Draw, Ellipse/circle (centered) и правой кнопкой мыши построим круг радиусом 8 мм с центром в начале координат. (рис. 3)

Рис. 2.− Задание границ рабочего поля

Рис. 3.−Построение исследуемой области

Мы начертили внешнюю оболочку кабеля. Теперь нужно начертить контур проводника и воздушного зазора. Аналогично чертим центрированные окружности радиусом 4,5 и 5 мм соответственно. (рис. 4)

Рис. 4.−Построение контура проводника и воздушного зазора

Программа дает каждой области наименование, в данном случае их здесь три C1, C2 и C3. В окне Set Formula можно исключить какое-либо множество из расчета. В данном случае все остается по умолчанию: C1+C2+C3.

Мы обозначили геометрию области расчета. Теперь необходимо задать краевые условия. Переходим к следующему этапу решения данной задачи. Выбираем панель Boundary, далее активируем Boundary Mode (рис. 5).

Рис. 5.− Выделение границы, для задания краевых условий

В открывшемся диалоговом окне (рис. 6) выбирают тип граничного условия (в рассматриваемой задаче – обобщенное условие Нейманна (Neumann)) и задают коэффициенты g (Heat flux, плотность теплового потока) и q (Heat transfer coefficient), коэффициент теплоотдачи. При указанных числовых значениях мы специфицировали граничное условие 3-го рода, приняв, что температура объекта T отсчитывается от температуры окружающей среды. Коэффициент теплопроводности тела (обозначение k в уравнении Boundary condition equation) будет определен далее при спецификации коэффициентов дифференциального уравнения.

Рис. 6.−Задание граничного условия

Граничные условия определены, и мы можем переходить к следующему этапу решения задачи – спецификации коэффициентов дифференциального уравнения. Для этого активируем панель PDE, и далее PDE Mode. Для удобства также отметим галочкой параметр Show Subdomain Labels (Показать метки подобластей) (рис. 7).

Рис. 7.− Спецификации коэффициентов дифференциального уравнения

1 – изоляционный материал;

2 – воздушный зазор;

3 – медный проводник.

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

Начнем с медного проводника (область 3). Спецификация коэффициентов представлена на рис. 8:

Рис.8.−Задание коэффициентов дифференциального уравнения

Из рисунка видны наши заданные величины – коэффициент теплопроводности k = 385(Вт/(м·K)), мощность источника тепла Q = 0,1кДж, а также внешняя температура Text = 20^OC.

Аналогичным образом поступим и с остальными областями 2 и 1 соответственно: (рис. 9, рис. 10)

Рис.9.−Задание коэффициентов для области 1

Рис. 10.−Задание коэффициентов для области 2

Всё, мы задали все необходимые параметры и можно искать решение задачи. Перейдем в меню Mesh (сетка) и построим её (рис. 11) Можно увеличивать её частоту, если необходима высокая точность решения.

Рис. 11.−Построение сетки на исследуемой области

Теперь активируем меню Solve (решение) и нажмем на Solve PDE. Решение задачи будет представлено графическом виде (рис. 12).

Рис. 12.− Решение задачи в графическом виде

Визуализация решения весьма гибко настраивается с помощью меню Plots и Parameters (Plot Selection) (рис. 13). Сделаем визуализацию решения трехмерной.

Рис. 13.−Трехмерное представление решения

2. Содержание отчета

1. Наименование работы, ее цель и краткие теоретические сведения.

2. Пошаговое решение поставленной задачи, с необходимыми пояснениями и в соответствии с вариантом.

3. Выводы по работе.

3.Контрольные вопросы

1. Коэффициенты теплопроводности различных веществ

2. Стационарная задача распределения тепла

3. Зависимость температуры проводника от параметров протекающего в нем тока

4. Теплопроводность воздуха

Варианты заданий:

Вариант	1	2	3	4	5
d1, мм	8	9	6	5	11
d2, мм	4,5	5	4,3	3,2	6
d3, мм	5	5,5	4,8	3,8	6,5
Материал проводника	медь	алюминий	алюминий	медь	алюминий
Материал изоляции	резина	ПВХ	полиэтилен	полипропилен	ПВХ

Лабораторная работа №3

СТАЦИОНАРНАЯ НЕЛИНЕЙНАЯ ЗАДАЧА РАСПРЕДЕЛЕНИЯ ТЕПЛА В ПРОВОДНИКЕ

Цель работы :Углубить знания о работе встроенного в MATLAB модуля PDE-Toolbox и методе конечных элементов. Решить нелинейную задачу распределения тепла в проводнике.

1. Пояснения к работе

В качестве второго примера рассмотрим распределение тепла в четырехпроводном кабеле. Геометрия исследуемой области показана на рис.1.

Рис. 1.− Геометрия исследуемой области

В данном случае следует учесть зависимость коэффициента теплопроводности от температуры вещества. Например, для зоны 2 (воздушный зазор): пусть коэффициент теплопроводности воздуха будет увеличиваться на 0,01 при повышении температуры на один градус (рис. 2).

Рис. 2.− Зависимость коэффициента теплопроводности от температуры

В поле k (Coeff. of heat conduction) – коэффициент теплопроводности. Здесь зададим зависимость, с которой этот коэффициент будет меняться с изменением температуры. Уравнение задается согласно синтаксису Matlab, значение внешней температуры записывается теперь не в ячейку External temperature, а в уравнение K1^.(u+20). С этого момента будет начинаться отсчет. Мы можем аналогичным образом задать зависимость для изменения мощности источника тепла и т.п.

Зададим также изменение коэффициента теплопроводности для остальных зон. Например, для медных проводников уравнение будет иметь вид, представленное в следующем окне (рис.3).

Рис. 3.− Задание коэффициента теплопроводности для медных проводников

Для изоляционного материала (рис. 4):

Рис. 4.− Задание коэффициента теплопроводности для изоляционного материала

Далее строим вычислительную сетку на заданной области (аналогично предыдущей работе). После этого нужно решить данную задачу. Т.к. она стала нелинейной, то необходимо использовать нелинейный решатель. Для этого активируем меню Solve, далее Parameters…(рис.5).

Рис. 5.−Задание параметров решения

По умолчанию здесь выставлено значение нелинейной погрешности равное 1e ^{− 4}. В данном случае система способна решить задачу с такой степенью точности. Если задать более сложные условия, то значение нелинейной погрешности можно уменьшить. Зададим теперь нужные значения мощности источников тепла в проводниках и можем запускать решение задачи (рис. 6).

Рис. 6

2. Содержание отчета

1. Наименование работы, ее цель и краткие теоретические сведения.

2. Пошаговое решение поставленной задачи, с необходимыми пояснениями, в соответствии с вариантом

3. Визуализация решения задачи в двумерном и трехмерном виде

4. Выводы по работе.

3. Контрольные вопросы

1. Нелинейные дифференциальные уравнения в частных производных

2. Зависимость теплопроводности различных веществ от их температуры

Варианты заданий: d₁ – диаметр кабеля, d₂ – диаметр внутренней части кабеля (без изоляции), d₃ – диаметр проводника

Вариант	1	2	3	4	5
d1, мм	8	9	15	13	11
d2, мм	6	7	10	9	8
d3, мм	5	5,5	4,8	3,8	6,5
Материал проводника	медь	алюминий	алюминий	медь	алюминий
Материал изоляции	резина	ПВХ	полиэтилен	полипропилен	ПВХ
Переменные коэффициенты (зона)	проводник	воздушный зазор	проводник	изоляция	воздушный зазор
K1	0,07	0,026	0,05	0,04	0,023
Внешняя температура	20	18	15	32	24

ЛИТЕРАТУРА

1. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: В 2–х томах. М.: Мир, 1990. 728 с.

2. Базаров И.П. Термодинамика. М.:Высшая школа. 1991. 376 с.

3. Вукалович М.П., Новиков И.И. Техническая термодинамика. М.: Энергия, 1968. 496 с.

4. Галлагер Р. Метод конечных элементов. Основы: Пер. с англ. − М.: Мир, 1984

5. Зенкевич О., Морган К. Конечные элементы и аппроксимация: Пер. с англ. − М.: Мир, 1986.

Код для цитирования: Скопировать