Анализ конечных элементов электростатическим образом приводимого в движение устройства MEMS

Этот пример показывает простой подход к двойному электромеханическому анализу конечных элементов электростатическим образом приводимого в движение микроэлектромеханического устройства (MEMS). Для простоты этот пример использует основанный на релаксации алгоритм, а не метод Ньютона, чтобы связать электростатическое и механические области.

Устройства MEMS

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

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

Этот пример использует следующую геометрию, чтобы смоделировать переключатель MEMS. Лучший электрод является 150 μm в длине и 2 μm в толщине. Модуль Молодежи E составляет 170 Гпа, и отношение Пуассона υ 0.34. Подовый электрод в электропечи является 50 μm в длине и 2 μm в толщине, и расположен 100 μm от крайнего левого конца лучшего электрода. Разрыв между верхними и нижними электродами является 2 μm.

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

Приблизьтесь для двойного электромеханического анализа

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

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

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

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

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

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

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

$| D_{d e f} (x) | \approx | D_{0} (x) | \frac{G}{G - v (x)},$

где $| D_{d e f} (x) |$ величина электрической плотности потока в деформированном электроде, $| D_{0} (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

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

Подвижный электрод: E3, E7, E11, E12
Фиксированный электрод: E4, E8, E9, E10
Доменный контур: E1, E2, E5, E6

Установите постоянные потенциальные значения 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) = ρ$ ,

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

$Δ V = 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 Гпа и отношение Пуассона $ν$ 0.34.

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

Задайте давление как граничную нагрузку на ребра. Давление имеет тенденцию вовлекать проводник в поле независимо от знака поверхностного заряда. Для определения CalculateElectrostaticPressure функционируйте, смотрите раздел Electrostatic Pressure Function в нижней части этой страницы.

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.5632e-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.8164e-07

Ссылки

[1] Sumant, P. S. Н. Р. Алуру и А. К. Кэнджеллэрис. “Методология для быстрого моделирования конечного элемента электростатическим образом приводимого в движение MEMS”. Международный журнал для численных методов в разработке. Vol 77, номер 13, 2009, 1789-1808.

Электростатическая функция давления

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

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

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

$E = - \nabla V$ .

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

$| D_{d e f} (x) | \approx | D_{0} (x) | \frac{G}{G - v (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

Документация