Этот пример показывает простой подход к двойному электромеханическому анализу конечных элементов электростатическим образом приводимого в движение микроэлектромеханического устройства (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. Вычислите загрузку и граничные условия для механического решения при помощи расчетных значений плотности заряда вдоль подвижного электрода. Электростатическим давлением на подвижный электрод дают
,
где величина электрической плотности потока и электрическая проницаемость рядом с подвижным электродом.
3. Вычислите деформацию подвижного электрода путем решения механической задачи FEA.
4. Обновите плотность заряда вдоль подвижного электрода при помощи расчетного смещения подвижного электрода,
где величина электрической плотности потока в деформированном электроде, величина электрической плотности потока в недеформированном электроде, расстояние между подвижными и фиксированными электродами в отсутствие приведения в действие, и смещение подвижного электрода в положении 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);
УЧП, управляющий этой проблемой, является уравнением Пуассона,
,
где коэффициент проницаемости и плотность заряда. Коэффициент проницаемости не влияет на результат в этом примере, пока коэффициент является постоянным. Предположение, что существует бесплатно в области, можно упростить уравнение Пуассона до уравнения Лапласа,
.
Задайте коэффициенты.
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
Задайте структурные свойства: модуль Молодежи 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.
Электростатическим давлением на подвижный электрод дают
,
где величина электрической плотности потока, электрическая проницаемость рядом с подвижным электродом, и величина электрического поля. Электрическое поле градиент электрического потенциала :
.
Решите механический FEA, чтобы вычислить деформацию подвижного электрода. Используя расчетное смещение подвижного электрода, обновите плотность заряда вдоль подвижного электрода.
где величина электрической плотности потока в деформированном электроде, величина электрической плотности потока в недеформированном электроде, расстояние между подвижными и фиксированными электродами в отсутствие приведения в действие, и смещение подвижного электрода в положении 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