Finite Element Analysis of Electrostatically Actuated MEMS Device

This example shows a simple approach to the coupled electromechanical finite element analysis of an electrostatically actuated micro-electromechanical (MEMS) device. For simplicity, this example uses the relaxation-based algorithm rather than the Newton method to couple the electrostatic and the mechanical domains.

MEMS Devices

MEMS devices typically consist of movable thin beams or electrodes with a high aspect ratio that are suspended over a fixed electrode.

Actuation, switching, and other signal and information processing functions can use the electrode deformation caused by the application of voltage between the movable and fixed electrodes. FEM provides a convenient tool for characterizing the inner workings of MEMS devices, and can predict temperatures, stresses, dynamic response characteristics, and possible failure mechanisms. One of the most common MEMS switches is the series of cantilever beams suspended over a ground electrode.

This example uses the following geometry to model a MEMS switch. The top electrode is 150 μm in length and 2 μm in thickness. The Young’s modulus E is 170 GPa, and the Poisson ratio υ is 0.34. The bottom electrode is 50 μm in length and 2 μm in thickness, and is located 100 μm from the leftmost end of the top electrode. The gap between the top and bottom electrodes is 2 μm.

A voltage applied between the top electrode and the ground plane induces electrostatic charges on the surface of the conductors which, in turn, leads to electrostatic forces acting normal to the surface of the conductors. Because the ground plane is fixed, the electrostatic forces deform only the top electrode. When the beam deforms, the charge redistributes on the surface of the conductors. The resultant electrostatic forces and the deformation of the beam also change. This process continues until the system reaches a state of equilibrium.

Approach for Coupled Electromechanical Analysis

For simplicity, this example uses the relaxation-based algorithm rather than the Newton method to couple the electrostatic and the mechanical domains. The example follows these steps:

1. Solve the electrostatic FEA problem in the nondeformed geometry with the constant potential V0 on the movable electrode.

2. Compute the load and boundary conditions for the mechanical solution by using the calculated values of the charge density along the movable electrode. The electrostatic pressure on the movable electrode is given by

P=12ϵ|D|2,

where |D| is the magnitude of the electric flux density and ϵ is the electric permittivity next to the movable electrode.

3. Compute the deformation of the movable electrode by solving the mechanical FEA problem.

4. Update the charge density along the movable electrode by using the calculated displacement of the movable electrode,

|Ddef(x)||D0(x)|GG-v(x),

where |Ddef(x)| is the magnitude of the electric flux density in the deformed electrode, |D0(x)| is the magnitude of the electric flux density in the undeformed electrode, G is the distance between the movable and fixed electrodes in the absence of actuation, and v(x)is the displacement of the movable electrode at position x along its axis.

5. Repeat steps 2–4 until the electrode deformation values in the last two iterations converge.

Electrostatic Analysis

In the electrostatic analysis part of this example, you compute the electric potential around the electrodes.

First, create the cantilever switch geometry by using the constructive solid geometry (CSG) modeling approach. The geometry for electrostatic analysis consists of three rectangles represented by a matrix. Each column of the matrix describes a basic shape.

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

Create a name for each basic shape. Specify the names as a matrix whose columns contain the names of the corresponding columns in the basic shape matrix.

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

Create a formula describing the unions and intersections of basic shapes.

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

Create the geometry by using the decsg function.

dl = decsg(gd,sf,ns);

Create a PDE model and include the geometry in the model.

model = createpde;
geometryFromEdges(model,dl);

Plot the geometry.

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

The edge numbers in this geometry are as follows:

  • Movable electrode: E3, E7, E11, E12

  • Fixed electrode: E4, E8, E9, E10

  • Domain boundary: E1, E2, E5, E6

Set constant potential values of 20 V to the movable electrode and 0 V to the fixed electrode and domain boundary.

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

The PDE governing this problem is the Poisson equation,

-(ϵV)=ρ,

where ϵ is the coefficient of permittivity and ρ is the charge density. The coefficient of permittivity does not affect the result in this example as long as the coefficient is constant. Assuming that there is no charge in the domain, you can simplify the Poisson equation to the Laplace equation,

ΔV=0.

Specify the coefficients.

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

Generate a relatively fine 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

Solve the model.

results = solvepde(model);

Plot the electric potential in the exterior domain.

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

Mechanical Analysis

In the mechanical analysis part of this example, you compute the deformation of the movable electrode.

Create a structural model.

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

Create the movable electrode geometry and include it in the model. Plot the geometry.

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

Specify the structural properties: the Young's modulus E is 170 GPa and the Poisson ratio ν is 0.34.

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

Specify the pressure as a boundary load on the edges. The pressure tends to draw the conductor into the field regardless of the sign of the surface charge. For the definition of the CalculateElectrostaticPressure function, see the Electrostatic Pressure Function section at the bottom of this page.

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

Specify that the movable electrode is fixed at edge 3.

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

Generate a 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

Solve the equations.

R = solve(structuralmodel);

Plot the displacement for the movable electrode.

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

Find the maximal displacement.

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

Repeatedly update the charge density along the movable electrode and solve the model until the electrode deformation values converge.

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

Plot the displacement.

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

Find the maximal displacement.

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

References

[1] Sumant, P. S., N. R. Aluru, and A. C. Cangellaris. “A Methodology for Fast Finite Element Modeling of Electrostatically Actuated MEMS.” International Journal for Numerical Methods in Engineering. Vol 77, Number 13, 2009, 1789-1808.

Electrostatic Pressure Function

The electrostatic pressure on the movable electrode is given by

P=12ϵ|D|2,

where |D|=ϵ|E| is the magnitude of the electric flux density, ϵ is the electric permittivity next to the movable electrode, and |E| is the magnitude of the electric field. The electric field E is the gradient of the electric potential V:

E=-V.

Solve the mechanical FEA to compute the deformation of the movable electrode. Using the calculated displacement of the movable electrode, update the charge density along the movable electrode.

|Ddef(x)||D0(x)|GG-v(x),

where |Ddef(x)| is the magnitude of the electric flux density in the deformed electrode, |D0(x)| is the magnitude of the electric flux density in the undeformed electrode, G is the distance between the movable and fixed electrodes in the absence of actuation, and v(x)is the displacement of the movable electrode at position x along its axis. Initially, the movable electrode is undeformed, v(x)=0, and therefore, |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