Численный анализ одного нелинейного эллиптического уравнения, возникающего при моделировании микроэлектромеханических систем Текст научной статьи по специальности «Математика»

УДК517.95 : 517.988

ЧИСЛЕННЫЙ АНАЛИЗ ОДНОГО НЕЛИНЕЙНОГО ЭЛЛИПТИЧЕСКОГО УРАВНЕНИЯ, ВОЗНИКАЮЩЕГО ПРИ МОДЕЛИРОВАНИИ МИКРОЭЛЕКТРОМЕХАНИЧЕСКИХ СИСТЕМ

КОНЧАКОВСКАЯ О.С, СИДОРОВ М.В._

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

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

Введение

Актуальность исследования. Современные исследования микросистемной техники (МСТ) сконцентрированы на создании и применении разнообразных методов математического моделирования, целью которых является описание механических (деформация), электромагнитных (проводимость, диэлектрическая и магнитная проницаемость), оптических и других свойств [5]. Создаваемые микроразмер ные устройства широко применяются в различных областях науки и техники: авиационной технике, автомобилестроении, робототехнике, атомной энергетике, системах связи, медицине.

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

Таким образом, разработка и совершенствование существующих средств математического моделирования и численного анализа задач, возникающих при проектировании МСТ, является актуальной научной задачей.

В данной работе рассмотрена проблема математического моделирования и численного анализа первой краевой задачи для нелинейного эллиптического урав-

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

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

Для достижения поставленной цели необходимо:

- сформулировать задачу математического моделирования и численного анализа процессов в микросистемном устройстве;

- разработать метод решения задачи расчета микроси-темного устройства;

- провести вычислительные эксперименты для различных параметров модели;

- провести анализ адекватности полученного решения.

Настоящая работа продолжает исследования, начатые в [2, 3].

1. Постановка задачи

Применение МСТ для различных областей физики обусловило разработку микроэлектромеханических (МЭМС), микрооптоэлектромеханических (МОЭМС), микроакустоэлектромеханических (МАЭМС) и других систем.

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

Тонкая металлизированная пленка |

1 V 1

диэлектрическая мембрана

Неподвижный электрод

Рис. 1. Схема работы простейшей электростатической МЭМС

Рассмотрим механическую систему, состоящую из упругого тела и приложенных к нему внешних «мертвых» сил, т.е. сил, сохраняющих величину и направление при деформациях системы; тело считаем закрепленным таким образом, что его перемещения как жесткого целого исключены [1].

Полная потенциальная энергия такой консервативной силы в нагруженном состоянии определяется суммой

Е = R + W ,

где R - потенциальная энергия деформации тела; W -потенциал внешних сил.

Электростатические методы активации устройств основаны на законе Кулона [9]: сила взаимодействия неподвижных зарядов q1 иq2 прямо пропорциональна произведению зарядов и обратно пропорциональна квадрату расстояния между ними, который в приведенных единицах измерения имеет вид

Б =

^ • q2

В общем случае, когда электрические заряды распределены неравномерно по двум электродам, расстояние между электродами равно

г = L + v(x),

где L > 0 - расстояние между двумя электродами при отсутствии деформаций мембраны; v(x) - величина деформации мембраны.

Потенциал внешних сил определяется равенством

w =

Ь :+V« '

здесь а>0 - константа, которая определяет относительную диэлектрическую проницаемость

среды; Ь - ограниченная область в Я2 .

Если в системе присутствуют упругие деформации, то V Ф 0 .

Потенциальная энергия деформации - это энергия, которая накапливается в теле при его упругой деформации.

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

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

Т 2

Р = 1Т ^х

Ь 2

где Т > 0 - константа напряжения.

Потенциальная энергия изменения объема - это энергия, накапливаемая за счет изменения объема рассматриваемого элементарного объема (одинакового изменения всех его размеров без искажения формы), которая равна

Q = I"

2 , Ь 2

где Б = —у ; ^ - толщина мембраны; у - модуль

3(1 -V2)

Юнга; V - коэффициент Пуассона.

Таким о бразом, полная потенциальная энергия системы равна

Е = Р + Q + W = | |Т|Уу|2 + БД^2 -^к

Ы2' ' 2' ' : + V.

По теореме Лагранжа состояние равновесия консервативной механической системы устойчиво тогда и только тогда, когда её полная потенциальная энергия минимальна [7]. Необходимое условие минимума функционала е (уравнение Эйлера-Остроградского) имеет вид

т. ^.2 а • V2

TДv - БД V =-

2 в П,

(1)

(Ъ + V)

v|sQ = 0. (2)

Поскольку нижний электрод тонкий, то его толщиной можно пренебречь относительно его размеров (Б = 0), и задача (1), (2) примет вид

ТДу =

в Ь,

а • V2

(Ъ + V)2

ЧдЬ = 0.

Тогда после замены V = -и получим нелинейную краевую задачу

-Ди = 4 1 2 в Ь, (Ъ - и)2

и|дЬ = 0,

где Дх^) =

а • V2 Т

Пусть Ъ = 1, Дхьх2) = X . Тогда получим задачу

-Д и = -

X

(1 - и)2

в Ь,

(3)

г

0 < U < 1 В Q , u|5Q =

2. Метод численного анализа

Для решения задачи (3), (4) воспользуемся методом двусторонних приближений [4, 6].

На конусе К неотрицательных в C(Q) функций введем в рассмотрение нелинейное операторное уравнение

u = Tu,

где

Tu = ds,

Q (1 - u(s))2

G(x, s) - функция Грина первой краевой задачи для оператора -Д.

Рассмотрим некоторые свойства оператора Т.

Оператор Т является непрерывным и монотонным на конусе К, так как

VUj, U2 e K Uj < U2 < 1 ^ Tuj < TU2 .

Ввиду монотонности оператора Т построим инвариантный конусный отрезок

< Vo,Wo >= {u| Vo < U< wq} ,

T < vq , wq >c< vq , wq > .

Для этого рассмотрим схему последовательных приближений

, ч „г G(x, s)

u n+i (x) = XJ-- ds, n = 0, 1, 2, .... (5)

Q (1 - Un (s))2

Положим vq = 0, wq = P= const e (0, 1). Тогда полу-

I Ч ,[ пг \ ds Хи0 (х) Wl (х) = X | G(x, 8)-2 = —^г

Ь (1 -Р)2 (1 -Р)2 '

V! (х) = X I G(x, = Хи0 (х) Ь

где %(х) = | G(x, . Ь

Необходимо выбором Р достичь выполнения условия Хи0(х) <в

Wl < Wo, т.е. (1 р)^ _ . Значит, р определяется системой неравенств

Х-тах | G( х, s)ds < Р(1 -Р)2, 0 <р<1. (6)

хеЬ 0

Пусть Ь - круг радиуса R. Функция Грина для круга радиусаR имеет вид

G(x, s) = — ln-----ln R

где х = (х1,Х2), 8 = (Sl,S2); точки 8 ,81 - точки, симметричные относительно окружности радиуса R ;

гх8, г 1 - расстояние между точками х и 8, х и 81

соответственно; р = ^/s2 +s -

R2

Тогда тахи0(х) = ——. Значит, система неравенств

хеЬ 4

(6) для определения Р принимает вид >2

Л R— <Р(1 -Р)2, 0 <р<1.

(7)

В работе [8] было доказано, что задача (3), (4) в

4 9

2 4

единичном круге в R имеет при 0 < X < — един-

ственное положительное решение, при

4

— < X < 0,593...

9

существует несколько положительных решений, а при Х> 0,593... - задача решений не имеет. Таким образом, задача (3), (4) в круге радиусаR имеет

единственное решение при X e I 0.

9R2 )

Вследствие монотонности оператора Т справедливы неравенства 0 = v0 < v1 < и* < w1 < w0 = р . Продолжая этот процесс, получаем

0 = v0 < v1 <... < vn <... < и* <

< ... < wn < ... < Wl < Wo = Р . Итак, имеет место следующая теорема.

Теорема. Оператор Tu = X J

G(x, s)

2

ds, где Q

Q (1 - u(s))

круг радиуса R, имеет при X е I 0;-I единствен-

\ 9R2)

ную неподвижную точку на конусном отрезке <Vo,Wo >, Wo =Р, Vo = 0, причемР определяется неравенствами (7). Последовательные приближения, формируемые по правилу (5), двусторонне сходятся к неподвижной точке оператора т .

3. Результаты численного анализа

Вычислительный эксперимент для задачи (3), (4) был проведен в круге радиуса R = 1. Тогда рассматриваемая задача эквивалентна нелинейному интегральному уравнению

u(x) = xj^xs)yds

Q (1 - u(s))2

2n rxs 2n Prxs1

(8)

и имеет единственное решение при X е ^0; — = 0,(4)J .

Для вычислительного эксперимента возьмем X = 0,3 . Тогда итерационный процесс решения уравнения (8) строим по схеме

и(к+1)(х) = 0,3 •!

G(x,s)

, (к) 2<1я, к = 0,1,2,.... (9) Ь (1 - и(к)(Ю)2

Для определения в воспользуемся соотношениями (7). Получим 0,090710 < в < 0,663889 .

На рис. 2 представлены верхние (непрерывная линия) и нижние приближения (пунктирная линия)

при п = 0, 1,2,3,4,5 в сечении Х2 = 0 .

то с точностью

0,25 -10 имеем

и « и 5 =

W5 + V5 2

i = 0 2, ф=—, 1 = 0, 5. Поверхность приближенно' 3 10

го решения и5(х) и его линииуровня представлены на рис. 3 и 4 соответственно.

Значения нормы приближенного решения задачи (3),

(4) с точностью е = 0,25 -10 5 в зависимости от параметра X представлены в табл. 3 и на рис. 5.

Рис. 2. Графикиw п (х^ 0) и vn (х^ 0) при п = 0, 1,2,3,4,5

В табл. 1 приведены значения wп (х) и vn (х) при п = 0, 1,2,3,4,5 в точках области Ь с полярными

координатами (Р;, Ф3), где р; = 0,2i, Ф3 = ~",

i = 0, 4,3 = 0, 5 (значения в остальных четвертях симметричны).

Как видно из табл. 1, в каждой из выбранных точек vn (х) < w п (х) значения w п (х) убывают, а значения vn (х) - возрастают. Поскольку

V.; < и* < w5, тах^5(х)-v5(x)) = 0,5-10-5,

хеЬ

Рис. 3. Поверхность приближенного решения ^(х)

Полученные результаты были доложены на XX Международном молодежном форуме «Радиоэлектроника и молодежь в XXI веке» (Харьков, ХНУРЭ, 19 -21 апреля 2016) и на 19 Всеукраинской (14 Международной) студенческой научной конференции по прикладной математике и информатике „СНКПМ1-2016" (Львов, ЛНУ им. И.Франко, 14 - 15 апреля 2016) [2, 3].

В табл. 2 приведены значения приближенного решения из(х) с точностью е = 0,25-10-5 в точках области Ь с полярными координатами (р;, Ф3), где р; = 0,25;,

Рис. 4. Линии уровня приближенного решения из (х)

ф 0 п п 3п 2п п п

р 10 5 10 5 2

0,090710 0,090710 0,090710 0,090710 0,090710 0,090710 0

0,090710 0,090710 0,090710 0,090710 0,090710 0,090710 1

0 0,086424 0,086424 0,086424 0,086424 0,086424 0,086424 2

0,085720 0,085720 0,085720 0,085720 0,085720 0,085720 3

0,085626 0,085626 0,085626 0,085626 0,085626 0,085626 4

0,085614 0,085614 0,085614 0,085614 0,085614 0,085614 5

0,090710 0,090710 0,090710 0,090710 0,090710 0,090710 0

0,084816 0,084832 0,084849 0,084848 0,084847 0,084824 1

0,25 0,080543 0,080560 0,080578 0,080580 0,080582 0,080562 2

0,079895 0,079912 0,079929 0,079931 0,079932 0,079912 3

0,079809 0,079827 0,079844 0,079845 0,079847 0,079827 4

0,079798 0,079816 0,079833 0,079835 0,079836 0,0798162 5

0,090710 0,090710 0,090710 0,090710 0,090710 0,090710 0

0,067997 0,068044 0,068090 0,068098 0,068106 0,068034 1

0,5 0,063980 0,064026 0,064072 0,064086 0,064100 0,064037 2

0,063485 0,063530 0,063576 0,063589 0,063602 0,063538 3

0,063422 0,063467 0,063513 0,063526 0,063539 0,063476 4

0,0634139 0,063459 0,063505 0,063518 0,0635313 0,063468 5

0 0 0 0 0 0 0

0,075000 0,075000 0,075000 0,075000 0,075000 0,075000 1

0 0,084250 0,084250 0,084250 0,084250 0,084250 0,084250 2

0,085438 0,085438 0,085438 0,085438 0,085438 0,085438 3

0,085590 0,085590 0,085590 0,085590 0,085590 0,085590 4

0,085609 0,085609 0,085609 0,085609 0,085609 0,085609 5

0 0 0 0 0 0 0

0,070126 0,070140 0,070154 0,070153 0,070152 0,070133 1

0,25 0,078864 0,078581 0,078598 0,078599 0,078600 0,078581 2

0,079640 0,079657 0,079674 0,079676 0,079677 0,079657 3

0,079777 0,079794 0,079811 0,079813 0,079814 0,079795 4

0,079794 0,079812 0,079829 0,079830 0,079832 0,079812 5

0 0 0 0 0 0 0

0,056221 0,056259 0,056297 0,056304 0,056310 0,056251 1

0,5 0,062513 0,062557 0,062602 0,062615 0,062627 0,062565 2

0,063299 0,063344 0,063389 0,063340 0,063416 0,063352 3

0,063398 0,063444 0,063489 0,063502 0,063516 0,063452 4

0,063411 0,063456 0,063502 0,063515 0,063528 0,063465 5

Таблица 2. Значения приближенного решения ^(х) в точках области О

\ Ф Р \ 0 п 10 п 7 3п 10 2п Т п 2

0 0,085612 0,085612 0,085612 0,085612 0,085612 0,085612

0,25 0,079796 0,079814 0,079831 0,079832 0,079834 0,079814

0,5 0,063412 0,063458 0,063504 0,063517 0,063530 0,063467

Таблица 3. Значения нормы приближенного решения в

зависимости от параметра X

X HI C X IH C X HI C

0,02 0,005038 0, 22 0,060318 0,42 0,128525

0,04 0,06 0,010153 0,015348 0, 24 0, 26 0,066437 0,072687 0,44 0,46 0,136425 0,144587

0,08 0,10 0,020631 0,026001 0, 28 0,30 0,079076 0,085611 0,48 0,50 0,153037 0,161804

0,12 0,14 0,031462 0,037020 0,32 0,34 0,092304 0,099164 0,52 0,54 0,170920 0,180429

0,16 0,18 0,042678 0,048446 0,36 0,38 0,106202 0,113431 0,56 0,58 0,190377 0,200825

0,20 0,054323 0, 40 0,120867 0,592 0,207367

iu lie

Рис. 5. Значения нормы приближенного решения в зависимости от параметра X

Выводы

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

Литература: 1. Алфутов Н.А. Основы расчета на устойчивость упругих систем. М.: Машиностроение, 1978. 312 с. 2. Кончаковская О.С. Численный анализ одной нелинейной краевой задачи, моделирующей микроэлектромеханическую систему // Радиоэлектроника и молодёжь в XXI веке: Материалы ХХ Международного молодёжного форума, 19 - 21 апреля 2016 г. Харьков: ХНУРЭ, 2016. Т. 7. С. 88 - 89. 3. Кончаковська О.С. Чисельний аналiз одше! нелшшно! крайово! задач^ яка моделюе мжроелек-тромехашчну систему // Студентська наукова конферен-щя з прикладно! математики та шформатики: Науковi пращ XIX Всеукрашсько! (XIV Мiжнародноl) конференцп, 14 - 15 квгтня 2016 р. Львш: ЛНУ ш. I. Франка, 2016. С. 78 - 79. 4. Красносельский М.А. Положительные решения операторных уравнений. М.: ГИФМЛ, 1962. 394 с. 5. Му-

хуров Н.И., Ефремов Г.И. Электромеханические микроустройства. Минск: Беларус. навука, 2012. 257 с. 6. Опой-цев В.И., Хуродзе Т.А. Нелинейные операторы в пространствах с конусом. Тбилиси: Изд-во Тбилис. ун-та, 1984. 246 с. 7. Циглер Г. Основы теории устойчивости конструкций. М.: Мир, 1971. 192 с. 8. Y. Guo and N. Ghoussoub. On the partial differential equations of electrostatic MEMS devices: stationary case. Submitted (2005). 9. F.Lin, Y. Yisong.Nonlinear non-local elliptic equation modeling electrostatic actuation. A Proceedings of The Royal Society. 05/2007; 463(2081) : 1323 - 1337.

Transliterated bibliography: 1. AlfutovN.A. Osnovy rascheta na ustojchivost' uprugih sistem. M.: Mashinostroenie, 1978. 312 p. 2. Konchakovskaja O.S. Chislennyj analiz odnoj nelinejnoj kraevoj zadachi, modelirujushhej mikrojelektromehanicheskuju sistemu // Radiojelektronika i molodjozh' v XXI veke: Materialy XXI Mezhdunarodnogo molodjozhnogo foruma, 19 - 21 aprelja 2016 g. Har'kov: KhNURE, 2016. T. 7. P. 88 - 89. 3. Konchakovs'ka O.S. Chysel'nyj analiz odnijei' nelinijnoi' krajovoi' zadachi, jaka modeljuje mikroelektromehanichnu systemu // Students'ka naukova konferencija z prykladnoi' mate-matyky ta informatyky: Naukovi praci XIX Vseukrai'ns'koi' (XIV Mizhnarodnoi') konferencii', 14 - 15 kvitnja 2016 r. L'viv: LNU im. I. Franka, 2016. P. 78 - 79. 4. Krasnosel'skijM.A. Polozhitel'nye reshenija operatornyh uravnenij. M.: GIFML, 1962. 394 p. 5. Muhurov N.I., Efremov G.I. Jelektromehanicheskie mikroustrojstva. Minsk: Belarus. navuka, 2012. 257 p. 6. Opojcev V.I.,Hurodze T.A. Nelinejnye operatory v prostranstvah s konusom. Tbilisi: Izd-vo Tbilis. un-ta, 1984. 246 p. 7. Cigler G. Osnovy teorii ustojchivosti konstrukcij. M.: Mir, 1971. 192 p. 8. Y. Guo and N. Ghoussoub. On the partial differential equations of electrostatic MEMS devices: stationary case. Submitted (2005). 9. F. Lin, Y. Yisong. Nonlinear non-local elliptic equation modeling electrostatic actuation. A Proceedings of The Royal Society. 05/2007; 463(2081) : 1323- 1337.

Поступила в редколлегию 30.04.2016

Рецензент: д-р физ.-мат. наук, проф. Литвин О.Н.

Кончаковская Оксана Сергеевна, студентка гр. СА-12-1

факультета прикладной математики и менеджмента ХНУРЭ. Научные интересы: математическое моделирование и вычислительная математика, программирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.

Сидоров Максим Викторович, канд. физ.-мат. наук, доцент каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, численные методы, математическая физика, теория R-функций и её приложения, стохастический анализ и его приложения. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.

Konchakovskaya Oksana Sergeevna, student of Faculty of Applied Mathematics and Management KhNURE. Research interests: mathematical modeling and computational mathematics, programming. Address: Ukraine, 61166, Kharkov, Lenin Ave, 14, phone +38 (057) 7021436.

Sidorov Maxim Victorovich, Ph.D. in Physics and Maths, associate professor, associate professor of Department of Applied Mathematics KhNURE. Research interests: mathematical modeling, numerical methods, mathematical physics, R-functions theory and its applications, stochastic analysis and its applications. Address: Ukraine, 61166, Kharkov, Lenin Ave, 14, phone +38 (057) 7021436.