Устройства МЭМС

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

При включении, переключении и других функциях обработки сигналов и информации может использоваться деформация электрода, вызванная приложением напряжения между подвижным и неподвижным электродами. КЭМ предоставляет удобный инструмент для характеристики внутренних рабочих характеристик МЭМС-устройств и может прогнозировать температуры, напряжения, динамические характеристики отклика и возможные механизмы отказов. Одним из наиболее распространенных MEMS-переключателей является ряд консольных пучков, подвешенных на заземлении.

В этом примере для моделирования MEMS-переключателя используется следующая геометрия. Верхний электрод имеет длину 150 мкм и толщину 2 мкм. Модуль Юнга E равен 170 ГПа, а коэффициент Пуассона («Пуассона») - 0,34. Нижний электрод имеет длину 50 мкм и толщину 2 мкм и расположен на расстоянии 100 мкм от крайнего левого конца верхнего электрода. Зазор между верхним и нижним электродами составляет 2 мкм.

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

Подход к сопряженному электромеханическому анализу

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

1. Решите проблему электростатической АМКЭ в недеформированной геометрии с постоянным потенциалом, V0 на подвижном электроде.

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

$$\frac{\mathrm{}}{\mathrm{}}{}^{\mathrm{P=12ϵ\text{'}D|2}}$$,

где $$|D$$ | - величина плотности электрического потока, $$а$$ ϵ - диэлектрическая проницаемость рядом с подвижным электродом.

3. Вычислите деформацию подвижного электрода, решив проблему механической АМКЭ.

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

где $$|{}_{}Ddef(x$$) | - величина плотности электрического потока в деформированном $$электроде{,}_{\mathrm{}}|D0$$ (x) | - величина плотности электрического потока в недеформированном $$$$электроде, G - расстояние между подвижным и неподвижным электродами при отсутствии срабатывания$$,$$v (x) - смещение подвижного электрода в положении x вдоль его оси.

5. Повторяйте шаги 2-4 до тех пор, пока значения деформации электрода в последних двух итерациях не сойдутся.

Электростатический анализ

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

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

rect_domain = [3,4,1.75e-4,1.75e-4,-1.75e-4,-1.75e-4,-1.7e-5,1.3e-5,1.3e-5,-1.7e-5]';
rect_movable = [3,4,7.5e-5,7.5e-5,-7.5e-5,-7.5e-5,2.0e-6,4.0e-6,4.0e-6,2.0e-6]';
rect_fixed = [3,4,7.5e-5,7.5e-5,2.5e-5,2.5e-5,-2.0e-6,0,0,-2.0e-6]';
gd = [rect_domain,rect_movable,rect_fixed];

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

ns = char('rect_domain','rect_movable','rect_fixed');
ns = ns';

Создайте формулу, описывающую соединения и пересечения основных форм.

sf = 'rect_domain-(rect_movable+rect_fixed)';

Создайте геометрию с помощью decsg функция.

dl = decsg(gd,sf,ns);

Создайте модель PDE и включите геометрию в модель.

model = createpde;
geometryFromEdges(model,dl);

Постройте график геометрии.

pdegplot(model,'EdgeLabels','on','FaceLabels','on')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-2e-4, 2e-4,-4e-5, 4e-5])
axis square

Номера кромок в этой геометрии следующие:

Установите постоянные значения потенциалов 20 В для подвижного электрода и 0 В для неподвижного электрода и границы области.

V0 = 0;
V1 = 20;
applyBoundaryCondition(model,'dirichlet','Edge',[4,8,9,10],'u',V0);
applyBoundaryCondition(model,'dirichlet','Edge',[1,2,5,6],'u',V0);
applyBoundaryCondition(model,'dirichlet','Edge',[3,7,11,12],'u',V1);

PDE, регулирующий эту задачу, является уравнением Пуассона,

$$$$-∇⋅(ϵ∇V)=ρ,

где $ϵ$ - это коэффициент диэлектрической проницаемости, а $$\mathrm{start}$$- плотность заряда. Коэффициент диэлектрической проницаемости не влияет на результат в этом примере, пока коэффициент является постоянным. Предполагая, что в области нет заряда, можно упростить уравнение Пуассона до уравнения Лапласа,

$$\mathrm{\Delta V}\mathrm{=}$$0.

Укажите коэффициенты.

specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0);

Создайте относительно тонкую сетку.

hmax = 5e-6;
generateMesh(model,'Hmax',hmax);
pdeplot(model)
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-2e-4, 2e-4,-4e-5, 4e-5])
axis square

Решите модель.

results = solvepde(model);

Постройте график электрического потенциала во внешней области.

u = results.NodalSolution;
figure
pdeplot(model,'XYData',results.NodalSolution,...
              'ColorMap','jet');

title('Electric Potential');
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-2e-4, 2e-4,-4e-5, 4e-5])
axis square

Механический анализ

В части механического анализа этого примера вычисляется деформация подвижного электрода.

Создайте модель несущей конструкции.

structuralmodel = createpde('structural','static-planestress');

Создайте геометрию подвижного электрода и включите ее в модель. Постройте график геометрии.

dl = decsg(rect_movable);
geometryFromEdges(structuralmodel,dl);
pdegplot(structuralmodel,'EdgeLabels','on')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-1e-4, 1e-4,-1e-5, 1e-5])
axis square

Задайте структурные свойства: модуль Юнга $$E$$ равен 170 ГПа, а коэффициент Пуассона («$$\mathrm{Пуассона}$$») равен 0,34.

structuralProperties(structuralmodel,'YoungsModulus',170e9, ...
                                     'PoissonsRatio',0.34);

Задайте давление как граничную нагрузку на кромки. Давление имеет тенденцию втягивать проводник в поле независимо от знака поверхностного заряда. Для определения CalculateElectrostaticPressure см. раздел Функция электростатического давления.

pressureFcn = @(location,state) - ...
                CalculateElectrostaticPressure(results,[],location);
structuralBoundaryLoad(structuralmodel,'Edge',[1,2,4], ...
                                       'Pressure',pressureFcn, ...
                                       'Vectorized','on');

Укажите, что подвижный электрод закреплен на кромке 3.

structuralBC(structuralmodel,'Edge',3,'Constraint','fixed');

Создайте сетку.

hmax = 1e-6;
generateMesh(structuralmodel,'Hmax',hmax);
pdeplot(structuralmodel);
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-1e-4, 1e-4,-1e-5, 1e-5])
axis square

Решите уравнения.

R = solve(structuralmodel);

Постройте график перемещения подвижного электрода.

pdeplot(structuralmodel,'XYData',R.VonMisesStress, ...
                        'Deformation',R.Displacement, ...
                        'DeformationScaleFactor',10, ...
                        'ColorMap','jet');

title('von Mises Stress in Deflected Electrode')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-1e-4, 1e-4,-1e-5, 1e-5])
axis square

Найдите максимальное смещение.

maxdisp = max(abs(R.Displacement.uy));
fprintf('Finite element maximal tip deflection is: %12.4e\n', maxdisp);

Finite element maximal tip deflection is:   1.5630e-07

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

olddisp = 0;
while abs((maxdisp-olddisp)/maxdisp) > 1e-10
% Impose boundary conditions
    pressureFcn = @(location,state) - ...
                    CalculateElectrostaticPressure(results, R, location);
    bl = structuralBoundaryLoad(structuralmodel,'Edge',[1,2,4], ...
                                                'Pressure',pressureFcn, ...
                                                'Vectorized','on');
% Solve the equations
    R = solve(structuralmodel);
    olddisp = maxdisp;
    maxdisp = max(abs(R.Displacement.uy));
    delete(bl)
end

Постройте график смещения.

pdeplot(structuralmodel,'XYData',R.VonMisesStress, ...
                        'Deformation',R.Displacement, ...
                        'DeformationScaleFactor',10, ...
                        'ColorMap','jet');

title('von Mises Stress in Deflected Electrode')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([-1e-4, 1e-4,-1e-5, 1e-5])
axis square

Найдите максимальное смещение.

maxdisp = max(abs(R.Displacement.uy));
fprintf('Finite element maximal tip deflection is: %12.4e\n', maxdisp);

Finite element maximal tip deflection is:   1.8162e-07

Ссылки

[1] Сумант, П. С., Н. Р. Алуру и А. К. Кангелларис. «Методология быстрого конечно-элементного моделирования электростатически активируемых МЭМС». Международный журнал численных методов в инженерии. Том 77, номер 13, 2009, 1789-1808.

Функция электростатического давления

Электростатическое давление на подвижный электрод задается

$$\frac{\mathrm{}}{\mathrm{}}{}^{\mathrm{P=12ϵ\text{'}D|2}}$$,

где $$|D|=ϵ\text{'}E|$$ - величина плотности электрического потока, $$ϵ$$ - электрическая диэлектрическая проницаемость рядом с подвижным электродом, и $|\mathit{}E$ | - величина электрического поля. Электрическое $\mathit{поле}$ E представляет собой градиент электрического $\mathit{потенциала}$V:

$\mathit{}\mathit{E=-\nabla V}$.

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

где $$|{}_{}Ddef(x$$) | - величина плотности электрического потока в деформированном $$электроде{,}_{\mathrm{}}|D0$$ (x) | - величина плотности электрического потока в недеформированном $$$$электроде, G - расстояние между подвижным и неподвижным электродами при отсутствии срабатывания$$,$$v (x) - смещение подвижного электрода в положении x вдоль его оси. Первоначально подвижный электрод $$\mathrm{является\; недеформированным}$$, v (x) = 0, и ${\mathit{}}_{\mathrm{}}\mathit{}поэтому{\mathit{|}}_{\mathrm{}}\mathit{}Ddef$(x) |≈|D0 (x) |.

function ePressure = ...
    CalculateElectrostaticPressure(elecResults,structResults,location)
% Function to compute electrostatic pressure.
% structuralBoundaryLoad is used to specify the pressure load
% on the movable electrode.
% Inputs:
% elecResults: Electrostatic FEA results
% structResults: Mechanical FEA results (optional)
% location: The x,y coordinate where pressure is obtained
%
% Output:
% ePressure: Electrostatic pressure at location
%
% location.x : The x-coordinate of the points
% location.y : The y-coordinate of the points
xq = location.x;
yq = location.y;

% Compute the magnitude of the electric field from the potential
% difference.
[gradx,grady] = evaluateGradient(elecResults,xq,yq);
absE = sqrt(gradx.^2 + grady.^2);

% The permittivity of vacuum is 8.854*10^-12 farad/meter.
epsilon0 = 8.854e-12;

% Compute the magnitude of the electric flux density.
absD0 = epsilon0*absE;
absD = absD0;

% If structResults (deformation) is available,
% update the charge density along the movable electrode.
if ~isempty(structResults)
    % Displacement of the movable electrode at position x
    intrpDisp = interpolateDisplacement(structResults,xq,yq);
    vdisp = abs(intrpDisp.uy);
    G = 2e-6; % Gap 2 micron
    absD = absD0.*G./(G-vdisp);
end

% Compute the electrostatic pressure.
ePressure = absD.^2/(2*epsilon0);

end