Table Of ContentМИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ
ЯДЕРНЫЙ УНИВЕРСИТЕТ «МИФИ»
Г.Е. Красников, О.В. Нагорнов, Н.В. Старостин
Моделирование физических процессов
с использованием пакета
Comsol Multiphysics
Рекомендовано УМО «Ядерные физика и технологии»
в качестве учебного пособия
для студентов высших учебных заведений
Москва 2012
УДК 519.85(075)
ББК 22.18я7
К78
Красников Г.Е., Нагорнов О.В., Старостин Н.В. Моделирование физических
процессов с использованием пакета Comsol Multiphysics: Учебное пособие. М.:
НИЯУ МИФИ, 2012. 184 с.
Предназначено для изучения среды математического моделирования Comsol
Multiphysics. В пособии подробно рассматриваются ключевые методы работы с
данной системой и разбираются конкретные типовые задачи. Также в книге со-
держится руководство по математическому программированию на Comsol Script и
особенности взаимодействия пакета Comsol Multiphysics с системой Matlab.
Данное пособие – первое руководство по Comsol Multiphysics на русском язы-
ке.
Полезно для студентов 3 и 4 курсов, изучающих курс математического моде-
лирования.
Подготовлено в рамках Программы создания и развития НИЯУ МИФИ.
Рецензент д-р физ.-мат. наук, проф. Н.А. Кудряшов (НИЯУ МИФИ)
ISBN 978-5-7262-1688-1 © Национальный исследовательский
ядерный университет «МИФИ», 2012
Редактор Т.В. Волвенкова
Подписано в печать 15.11.2011. Формат 60х84 1/16
Уч.-изд. л. 11,75. Печ. л. 11,5. Тираж 120 экз.
Изд. № 1/57. Заказ № 21.
Национальный исследовательский ядерный университет «МИФИ».
115409, Москва, Каширское ш., 31
ООО «Полиграфический комплекс «Курчатовский».
144000, Московская область, г. Электросталь, ул. Красная, 42
ВВЕДЕНИЕ
Программа Comsol Multiphysics 3.5a (FEMLAB) – математиче-
ский пакет, предназначенный для численного решения задач раз-
личных областей физики. Пакет основан на методе конечных эле-
ментов, с помощью которого производятся все вычисления. Во
время выхода версии 3.2, компания Comsol переименовала пакет
FEMLAB в Comsol Multiphysics. В настоящей работе для упроще-
ния будет использоваться как новое, так и старое названия. Про-
грамма Comsol Multiphysics – это мощное средство для решения
сложных задач, сопровождаемых большими объёмами вычислений.
Возможность решать тот или иной класс задач реализована в виде
специальных прикладных режимов, при загрузке которых автома-
тически выбирается нужная система уравнений, в которой необхо-
димо только задать коэффициенты и граничные условия. В пакете
FEMLAB доступны для решения классы задач электростатики,
электродинамики, электромагнетизма, акустики, теплопереноса,
теории упругости, гидродинамики, а также классические диффе-
ренциальные уравнения, такие как уравнение Шрёдингера, уравне-
ние Гельмгольца, уравнение теплопроводности, волновое уравне-
ние и другие. В вышеуказанных классах задач допускается решать,
стационарные, временные задачи, задачи на собственные значения,
а также параметрические задачи.
1. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ
1.1. Теоретическое введение
Рассмотрим следующую задачу Коши для обыкновенного диф-
ференциального уравнения:
⎧u''(x)=q(x);
⎪ (1.1)
⎪x∈[a;b];
⎨
u(a)=0;
⎪
⎪⎩u'(b)=0;
3
где q(x) – произвольная известная функция от x. В общем виде
уравнение не решается аналитически. Однако можно упростить
задачу, потребовав выполнения решения уравнения не на всей не-
прерывной области [а; b], а лишь на определённом наборе проме-
жуточных точек. Сделаем это следующим образом.
Разобьём область решения [а; b] на конечное число отрезков
[x;x ], на каждом из которых выберем базисную функцию ϕ, по
i i+1 i
которым бы можно было разложить исходную функцию N .
u= ∑uϕ
i i
i=1
Функции должны быть такого вида, чтобы по ним можно было раз-
ложить любую функцию u(x), а также чтобы они образовывали ба-
зис на области решения, т.е. удовлетворяли условию: .
(ϕ;ϕ)=0
i j
Наиболее простыми для расчётов являются кусочно-линейные
функции в виде треугольников:
⎧ x−x
i−1 ;x∈[x ;x ]
⎪⎪x −x i−1 i
ϕ =⎨ i i−1 . (1.2)
i x −x
⎪ i+1 ;x∈[x ;x ]
⎪⎩x −x i i+1
i+1 i
Для упрощения вычислений выберем все элементы одинаково-
го размера h.
R(x)=u''(x)−q(x)=0 – невязка исходного уравнения.
Для решения задачи, необходимо минимизировать её по методу
Галёркина, умножая скалярно R(x) на каждую из базисных функ-
ций:
b b
(R(x);ϕ(x))=∫R(x)ϕ(x)dx=∫(u''(x)−q(x))ϕ(x)dx=0; (1.3)
i i i
a a
i∈1..N.
Преобразуем каждый интеграл по отдельности.
4
b b
∫u''(x)ϕ (x)dx=∫[(u'(x)ϕ (x))'−u'(x)ϕ '(x)]dx=
i i i
a a
b
=u'(x)ϕ (x)b −∫u'(x)ϕ '(x)dx.
i a i
a
Первое слагаемое равно нулю, поскольку все базисные функ-
ции обращаются в 0 на концах отрезка. Второй интеграл можно
разложить на две части
b xi 1 xi+1 1
∫u'(x)ϕ '(x)dx= ∫ u'(x) dx+ ∫ u'(x)(− )dx=
i h h
a xi−1 xi
1 u −2u +u
= (u −u −(u −u ))=− i+1 i i−1 .
h i i−1 i+1 i h
Таким образом получаем систему линейных алгебраических
уравнений с матрицей следующего вида:
⎛ 1 0 0 0 0 ⎞⎛ u ⎞ ⎛ f ⎞
⎜ ⎟⎜ 1 ⎟ ⎜ 1 ⎟
⎜− 1h 2h − 1h .. .. ⎟⎜ u ⎟ ⎜ f ⎟
2 2
⎜ .. .. .. .. .. ⎟⎜ .. ⎟=⎜ .. ⎟, (1.4)
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎜ .. .. − 1h 2h − 1h⎟⎜ u ⎟ ⎜ f ⎟
⎜⎜⎝ 0 0 0 − 1h 1h ⎟⎟⎠⎜⎝uNN+1⎟⎠ ⎜⎝ fNN+1⎟⎠
b
где, f = f = 0; f = ∫q(x)ϕ (x)dx;i = 2..N.
N+1 1 i i
a
Решив систему методом прогонки, получим коэффициенты
u ...u , которые задают исходный вид функции u(x) на опреде-
1 N+1
лённом наборе точек. При N →∞численное решение будет стре-
миться к точному.
5
Основное положение метода конечных элементов – переход от
решения уравнения на непрерывной области, к требованию выпол-
нения системы уравнений на ограниченном числе точек.
Уравнение (1.3), означающее ра-
венство невязки нулю, является сла-
бой формой дифференциального
уравнения. Оно означает выполне-
ние уравнения на определённом на-
боре точек.
В случае задачи двух перемен-
ных, в качестве конечных элементов
стоит выбрать шестигранные пира-
миды, как показано на рис.1.1, 1.2.
Рис. 1.1. Прямоугольная Линейная базисная функция будет
область, разбитая на конечные иметь вид шестигранной пирамиды с
элементы
вершиной в центре элемента (x;y).
i i
Вне элемента функция равна нулю. Далее раскладываем искомую
функцию по базису:
⎧ y−y
ϕ(1) = j−1;(x,y)∈S
⎪
ij b 1
⎪
⎪ x −x y−y
ϕ(2) =( i+1/2 + j−1/2);(x,y)∈S
⎪ ij a b 2
⎪
⎪ϕ(3) = xi+1−x;(x,y)∈S
⎪ ij b 3
⎪ (1.5)
⎪ y −y
ϕ =⎨ϕ(4) = i+1 ;(x,y)∈S
ij ij b 4
⎪
⎪ x−x y −y
⎪ϕ(5) =( i−1/2 + j+1/2 );(x,y)∈S
ij a b 5
⎪
⎪ x−x
ϕ(6) = i−1;(x,y)∈S
⎪ ij b 6
⎪
ϕ =0;(x,y)∉S
⎪ ij
⎪
⎩
6
N M
u(x,y)=∑∑αϕ(x,y). (1.6)
ij ij
i=1 j=1
Метод конечных элементов можно
распространить на систему диффе-
ренциальных уравнений. В этом слу-
чае необходимо разложить каждую
искомую переменную.
Представленные выше элементы
являются элементами первого порядка
или линейными элементами. Значения
Рис. 1.2. Шестиугольный
функции вычисляются только на уз-
конечный элемент
лах элемента. Базисные функции
имеют линейный вид. Помимо этих элементов существуют элемен-
ты второго, третьего и более высокого порядков. Вместо линейных
функций на них используются параболические, кубические
функции и функции более высоких порядков. Значения вычисля-
ются не только в узлах элемента, но и на его ребрах. Разложение
искомой функции на элементах второго порядка имеет гораздо бо-
лее сложный вид, чем (1.2). Использование элементов первого по-
рядка возможно, только если оператор L[u] содержит производные
только первого порядка. При нали-
чии производных третьего и более
высоких порядков в уравнении не-
обходимо использовать элементы
соответствующих порядков.
Вместо шестиугольных элемен-
тов в МКЭ могут также использо-
ваться треугольные (рис. 1.3, 1.4).
При использовании таких эле-
Рис. 1.3. Сетка из треугольных
ментов следует перейти к новым
элементов, сгненерированая
переменным ξ ;ξ ;ξ . Переход осу- cистемой FEMLAB
1 2 3
ществляется следующим образом:
7
⎛ 1 1 1 ⎞⎛ξ⎞ ⎛1⎞
⎜ ⎟⎜ 1⎟ ⎜ ⎟
⎜x x x ⎟⎜ξ ⎟=⎜x⎟. (1.7)
1 2 3 2
⎜ ⎟⎜ ⎟ ⎜ ⎟
y y y ξ y
⎝ ⎠⎝ ⎠ ⎝ ⎠
1 2 3 3
Необходимое условие на новую систему координат: ξ +ξ +ξ =1
1 2 3
В новых координатах функция раскладывается следующим об-
разом:
u(x,y)=u(x ,y )ξ+u(x ,y )ξ +u(x ,y )ξ. (1.8)
1 1 1 2 2 2 3 3
Разбиение на треугольные элементы (триангуляция) использу-
ется в системе FEMLAB.
1.2. Виды конечных элементов
В некоторых случаях влияние на
решение задачи может оказывать
выбор базисных функций. Рассмот-
рим основные типы таких функций.
Лагранжевы элементы – Наи-
более часто используемые базисные
функции в FEMLAB. Лагранжев
элемент k-го порядка задаёт значение
Рис. 1.4. Треугольный конечный функции на границах элемента, а
элемент
также на k точках самого элемента.
Например, лагранжев элемент второго порядка задаёт значения
функции на всех узлах конечного элемента, а также на середине
каждой стороны элемента. Подобная постановка условия позволяет
считать производные величины u второго порядка по пространст-
венным координатам.
Лагранжевы элементы второго порядка позволяют решать по-
давляющее большинство дифференциальных уравнений. В случае,
если уравнение содержит производные более высоких порядков,
следует выбирать порядок элементов, соответствующий макси-
мальной степени производной в уравнении.
8
Аргирисовы элементы. При аппроксимации переменной ла-
гранжевыми элементами возможны ситуации, когда первые произ-
водные на границах между элементами будут бесконечны. В неко-
торых уравнениях это может быть серьёзной проблемой для реше-
ния. В этом случае можно использовать аргирисовы элементы, ко-
торые задают конечные производные на границах и вторые произ-
водные в углах элемента. На каждом элементе искомая функция u
аппроксимируется полиномом 5-й степени в локальных координа-
тах ξ.
Аргирисовы элементы доступны только на треугольных или
тетраэдных (3D-случай) сетках.
Эрмитовы элементы. Эрмитовы элементы сходны с элемен-
тами Лагранжа за тем исключением, что в углах элемента опреде-
ляются непосредственно не значения функции, а определяются
лишь значения первых производных. Во внутренних точках ребра
элемента, как и в случае лагранжева элемента, определяются зна-
чения функции. Минимальный возможный порядок эрмитова эле-
мента k=3.
Пузырьковые элементы. Пузырьковые элементы (bubble ele-
ments) называются так благодаря своей форме. На границе области,
тестовая функция обращается в ноль, а в центре элемента достигает
максимума. Элементы доступны на любых сетках.
Помимо вышеперечисленных, в некоторых прикладных режи-
мах в системе FEMLAB используются следующие элементы: век-
торные, роторные, дивергентные, дискретные.
2. НАЧАЛО РАБОТЫ С FEMLAB
2.1. Установка
Пакет Comsol Multiphysics устанавливается с четырёх дисков.
При автоматическом запуске первого диска появляется следующее
окно (рис. 2.1):
9
Рис. 2.1. Начало инсталляции Comsol Multiphysics 3.5a
Выберем пункт New Installation. Следующим этапом осущест-
вляется установка лицензии программы (рис. 2.2). Лицензия к про-
грамме может поставляться в виде специального лицензионного
кода, либо содержаться в файле.
Рис. 2.2. Выбор лицензии для программы
10