Technical Articles and Newsletters

Accelerating Finite Element Analysis in MATLAB with Parallel Computing




MEMS Devices









Figure 4. Cantilever switch modeled geometry.


Applying FEM to Coupled Electro-Mechanical Analysis




The electrostatic pressure on the movable electrode is given by \[P=\frac{1}{2\epsilon}|D|^2\]

\(|D|\) = Magnitude of the electric flux density
\(\epsilon\) = Electric permittivity next to the movable electrode



\[\left|D_{\mathrm{def}}(x)\right| \approx \left|D_0(x)\right|\frac{G}{G - v(x)}\]

\(|D_{\mathrm{def}}(x)|\) = Magnitude of the electric flux density in the deformed electrode
\(|D_0(x)|\) = Magnitude of the electric flux density in the undeformed electrode
\(G\) = Distance between the movable and fixed electrodes in the absence of actuation
\(v(x)\) = Displacement of the movable electrode at position x along its axis


The electrostatic analysis (step 1) involves five steps:

  • Draw the cantilever switch geometry.
  • Mesh the exterior domain.
  • Solve for the unknown electric potential in the exterior domain.
  • Plot the solution (Figure 5).


Mechanical Analysis

The mechanical analysis involves five steps:

  • Meshing the domain
  • Deriving element-level equations
  • Assembly
  • Imposing boundary conditions
  • Solving the equations

Meshing the Domain


Deriving Element-Level Equations

In this step, we assume an approximate solution for the differential equations over an element at selected points, or nodes. In our example, the solution is determined in terms of discrete values of \(\phi\) displacements in the \(x\) and \(y\) directions (2D analysis). The number of unknown primary fields at a node is called the degrees of freedom (DOF) at that node.

The governing differential equation is now applied to the domain of a single element (Figure 6). At the element level, the solution to the governing equation is replaced by a continuous function approximating the distribution of \(\phi\) over the element domain \(D^e\), expressed in terms of the unknown nodal values \(\phi_1\), \(\phi_2\), and \(\phi_3\) of the solution \(\phi\). A system of equations in terms of \(\phi_1\), \(\phi_2\), and \(\phi_3\) can then be formulated for the element.

Figure 6. FEA - Discretized domain and a single element.


We obtain the solution equations for the system by combining the solution equations of each element to ensure continuity at each node. In our example, the element-level stiffness and force matrices (\(K_e\) and \(F_e\)) are assembled to create a global stiffness matrix (\(K\)) and a force matrix (\(F\)) over the entire domain

Imposing Boundary Conditions

We impose the necessary boundary conditions at the boundary nodes and solve the global system of equations. The solution \(\phi(x,y)\) to the problem becomes a piecewise approximation, expressed in terms of the nodal values of \(\phi\). A system of linear algebraic equations results from the assembly procedure: \(K \phi = F\).


Solving the Equations

We solve the global system of equations \(K \phi = F\) for the displacements in the \(x\) and \(y\) directions.

Parallelizing the Problem


This step involves three main operations for each element:
1. Compute the element stiffness matrix (\(K_e\)) from the element-level solution equations. \(K_e\) is of size \(N_e \times N_e\)

\(N_e\) = \(n \times D\)
\(N_e\) is the number of DOF per element
\(n\) is the number of nodes per element
\(D\) is the number of DOF per node

2. Map the local positions of the \(K_e\) matrix values to their position in the global stiffness matrix.
3. Populate the global stiffness matrix (\(K\)) using the map with the element stiffness matrix values.

Number of Elements Degrees of Freedom (DOF) Assembly Time (seconds) Total Execution Time (seconds) Assembly Time/ Total Execution Time (%)
528 806 4.09 4.82 84.5
6584 7690 45.073 46.371 97.2
53550 55882 960.069 1005.448 95.5
460752 469566 64573.472 64616.662 99.9


As the number of elements increases, so do the number of iterations and the size of \(K\). Figure 7 and Table 1 compare the time taken for assembly to the total execution time in serial mode for the system with different DOF. Clearly, assembly is the most time-consuming portion, taking up more than 99% of the total execution time when the system has high DOF.


Stiffness Matrix Assembly

Fortunately, the final assembled \(K\) matrix is independent of the order in which the elements are chosen within the loop. We can evaluate the contributions of several elements to global stiffness matrix (\(K\)) simultaneously by distributing the computations across multiple MATLAB workers. The assembly operations are normally executed in a serial for-loop, which steps through each element and determines its contribution to the global stiffness matrix. We simply convert the serial for-loop to a parallel for-loop using the parfor construct in Parallel Computing Toolbox™ (Figure 8).

Figure 8. Serial and parallel approaches to stiffness matrix assembly.

In our example, the global stiffness matrix is the sum of the contributions of the element stiffness matrices across all the iterations of the loop. The parfor construct enables us to handle such reduction assignments (usually of the form \(r = f(\text{expr}, r)\)) automatically.




Degrees of Freedom (DOF) Total Execution Time (seconds)
Serial Mode Parallel Mode
150 0.53 1.15
250 0.77 1.18
806 4.82 1.37
2200 12.93 2.8
7690 46.37 6.1
28546 355.04 20.92
103822 3406.61 129.23
218862 14871.7 496.47
469566 64616.66 1911.71
866122 218674.88 6237.91

Table 2. Representative data used in Figures 9a and 9b.







Published 2010 - 91826v00

View Articles for Related Capabilities

View Articles for Related Industries