Устройства MEMS

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

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

Этот пример использует следующую геометрию, чтобы смоделировать переключатель MEMS. Верхний электрод имеет длину 150 мкм и толщину 2 мкм. Модуль Е Юнга является 170 GPa, и отношение Пуассона 0,34. Нижний электрод имеет длину 50 мкм и толщину 2 мкм и расположен на расстоянии 100 мкм от крайнего левого конца верхнего электрода. Зазор между верхним и нижним электродами составляет 2 мкм.

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

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

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

1. Решите электростатическую задачу FEA в недеформированной геометрии с постоянным потенциалом, V0 на подвижном электроде.

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

$$P=\frac{1}{2ϵ}|D{|}^2$$,

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

3. Вычислите деформацию подвижного электрода, решив механическую задачу FEA.

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

где $$|D_{def}\left(x\right)|$$ - величина плотности электрического потока в деформированном электроде, $$|D_0\left(x\right)|$$ - величина плотности электрического потока в недеформированном электроде, $$G$$ - расстояние между подвижным и неподвижным электродами при отсутствии приведения в действие, и $$v\left(x\right)$$- перемещение подвижного электрода в положении 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);

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

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);

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

$$-\nabla \cdot (ϵ\nabla V)=\rho$$,

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

$$\Delta V=0$$.

Задайте коэффициенты.

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

Сгенерируйте относительно мелкий mesh.

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 ГПа и коэффициент Пуассона $$\nu$$ составляет 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');

Сгенерируйте mesh.

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] Sumant, P. S., N. R. Aluru, and A. C. Cangellaris. Методология быстрого конечного Элемента моделирования МЭМС с электростатическим приводом. International Journal for Numerical Methods in Engineering. Vol 77, Число 13, 2009, 1789-1808.

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

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

$$P=\frac{1}{2ϵ}|D{|}^2$$,

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

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

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

где $$|D_{def}\left(x\right)|$$ - величина плотности электрического потока в деформированном электроде, $$|D_0\left(x\right)|$$ - величина плотности электрического потока в недеформированном электроде, $$G$$ - расстояние между подвижным и неподвижным электродами при отсутствии приведения в действие, и $$v\left(x\right)$$- перемещение подвижного электрода в положении x вдоль его оси. Первоначально подвижный электрод недеформирован, $$v\left(x\right)=0$$и, следовательно, $\left|{\mathit{D}}_{\mathrm{def}}\left(\mathit{x}\right)\right|\approx \left|{\mathit{D}}_0\left(\mathit{x}\right)\right|$.

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