COMSOL 6.4 - Vibrating Micromirror with Viscous and Thermal Damping

Vibrating Micromirror with Viscous and Thermal Damping

Micromirror systems are part of small MEMS (microelectromechanical systems) systems for integrated scanning micromirror systems in telecommunications and consumer applications. The mirrors have low power consumption and the manufacturing cost is low. Their applications range from microscopy applications, medical imaging, scanners, and heads-up displays, to some fiber optics applications.

Figure 1: Geometry of the micromirror system including the fixed constraint, the torquing force components F, and the symmetry plane.

The dynamic vibrating behavior of the micromirror (or other MEMS resonators) and especially the damping is an important factor affecting the design of their actuator systems. In the development process of new micromirror systems, modeling the damping correctly, and thus predicting the performance of the system, can save both time and money.

This model presents a highly simplified micromirror geometry modeled with the Shell physics interface coupled to Thermoviscous Acoustics, Frequency Domain. The model compares the simulated damped vibration results to a simple oscillator model tuned using the results of an eigenfrequency study. For a more advanced micromirror model see the Prestressed Micromirror Vibrations: Thermoviscous–Thermoelasticity Coupling tutorial model also found in the Application Library.

Model Definition

A sketch of the idealized micromirror system modeled here is shown in Figure 1. The mirror is made of silicon and is surrounded by air. It is 0.5 mm by 0.5 mm (only half is shown in the figure because of symmetry) and has a thickness of 1 μm.

The solid mirror is modeled using shell elements and the surrounding air is modeled using thermoviscous acoustics. A pressure acoustics layer is used to truncate the computational domain. The model combines the Shell interface, the Thermoviscous Acoustics, Frequency Domain interface, and the Pressure Acoustics, Frequency Domain interface. The three physics interfaces are set up and combined using multiphysics couplings: the Acoustic–Thermoviscous Acoustics Boundary and the Thermoviscous Acoustic–Structure Boundary.

One goal with the model is to get a correct assessment of the damping the mirror experiences. Therefore, the detailed thermoviscous acoustics is used to model the air domain surrounding the micromirror. The interface includes thermal and viscous damping explicitly as it solves the full linearized Navier–Stokes, continuity, and energy equations, implemented in the Thermoviscous Acoustics, Frequency Domain interface.

An important parameter in such models is the thickness of the viscous and thermal boundary layers (also called the penetration depths). It is in these layers that the energy is dissipated (viscous drag and thermal conduction). For air, the layers are of approximately the same size. The thickness of the viscous boundary layer δv is given by the expression

(1)

since δv = 0.22 mm at 100 Hz. The boundary layer thickness can also be evaluated in postprocessing through the variables tash.d_visc and tash.d_therm.

The model is both solved using a frequency-domain sweep and using an eigenfrequency study. The frequency sweep results in a frequency response where the displacement (or velocity) is evaluated at the tip of the mirror (on the symmetry plane).

The response at the resonance yields the Q-factor Qr and the resonance frequency fr. Their expected values are computed based on the eigenfrequency analysis results. The Q-factor is defined as

where f0 is the resonance frequency and Δf is the peak width at half power (which corresponds to the 3 dB down width). The value of the band width is evaluated in the results using the functionality, see Figure 7 (bottom). The Q-factor is related to the attenuation rate α through

The result of the eigenfrequency study is a complex-valued eigenfrequency fc which relates to the other parameters through

(2)

Finally, the response of an underdamped harmonic oscillator can be written as

(3)

where F is the actuating force and v is the velocity of the oscillator. This expression is used as a model equation for the frequency response in the vicinity of the resonance frequency, based on parameters identified in the eigenfrequency analysis (using the relations in Equation 2).

Results and Discussion

The first two modes from the eigenfrequency study are shown in Figure 2. The torquing mode studied in the subsequent frequency domain analysis is the top figure where fc is at about 10.5 kHz, while a symmetric vibrating mode is found at about 13 kHz. The next mode (not shown) is at nearly 40 kHz. The modal von Mises stress for the first mode is shown in Figure 3. It is important to note the regions of high stress, these are also the regions where a finer mesh needs to be used to capture the full dynamics.

Figure 2: Displacement plot of the torquing (top) and symmetric (bottom) vibrating modes of the micromirror.

Figure 3: Von Mises stress distribution in the shell for first (torquing) mode.

The displacement of the micromirror with the torquing actuation is shown in Figure 4 at a frequency of 10 kHz. The acoustic temperature variations and the acoustic pressure distribution also at 10 kHz is shown in Figure 5. The instantaneous local velocity (acoustic variations in the velocity) is shown in Figure 6. A high velocity region is seen near the edge of the micromirror. The extent of this region into the surrounding air is given by the scale of the viscous boundary layer (also known as the viscous penetration depth, Equation 1). The thermal boundary layer (thermal penetration depth) can in the same way be identified in Figure 5 (top).

The displacement response of the system is plotted in Figure 7 (top). This plot shows the absolute value of the z-component of the displacement field |wshell| evaluated at the tip of the micromirror. The bottom part of this figure shows the power response, on the dB scale (normalized by the maximum), by plotting

and comparing to the model expression given in Equation 3. The parameters used in the analytical model expression stem from the eigenfrequency derived in the eigenfrequency study. The value of the eigenfrequency fc and the other parameters used in the analytical model are defined as variables (see the > > ). The values of the variables evaluated for the first eigenfrequency, are computed in the Evaluation Group in the model. Using the computed value leads to good agreement between the two power response curves, as seen in Figure 7 (bottom). The figure also shows the bandwidth of the frequency response, this value agrees well with the value extracted from the eigenfrequency analysis and evaluated in the Evaluation Group.

Figure 4: Displacement of the micromirror at 10 kHz subject to the torquing force.

Figure 5: (top) Temperature field in the thermoviscous acoustics domain around the micromirror, and (bottom) pressure isosurfaces showing the antisymmetric pressure distribution.

Figure 6: The instantaneous local velocity, with clear high velocity regions near the mirror edges. The extent of the high velocity region is given by the scale of the acoustic viscous boundary layer.

Figure 7: (top) Displacement response, and (bottom) power response with half-power bandwidth and analytical model curve.

Acoustics_Module/Vibrations_and_FSI/vibrating_micromirror

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file vibrating_micromirror_parameters.txt.

Geometry 1

The geometry sequence for the model (see Figure 1) is available in a file. If you want to create it from scratch yourself, you can follow the instructions in the Appendix — Geometry Modeling Instructions section. Otherwise, insert the geometry sequence as follows:

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file vibrating_micromirror_geom_sequence.mph.

In the toolbar, click

Click the

button in the toolbar.

Click the

button in the toolbar.

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file vibrating_micromirror_variables.txt.

Shell

In the toolbar, click

In the window for , type Shell in the text field.

Locate the section. From the list, choose .

Select Boundary 9 only.

Symmetry

In the toolbar, click

In the window for , type Symmetry in the text field.

Locate the section. From the list, choose .

Select Boundaries 1, 2, 6, 11, and 15 only.

Symmetry Edge

In the toolbar, click

In the window for , type Symmetry Edge in the text field.

Locate the section. From the list, choose .

Select Edge 9 only.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Silicon (mat2)

In the window for , locate the section.

From the list, choose .

Shell (shell)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Thickness and Offset 1

In the window, under > click .

In the window for , locate the section.

In the d0 text field, type h_mirror.

Fixed Constraint 1

In the toolbar, click

and choose .

Select Edge 13 only.

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Body Load 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the fV vector as


0	x
0	y
x/(0.25[mm])*1e5[N/m^3]	z

This corresponds to a torquing force proportional to the x-coordinate acting in the z direction. A load of 0 N/m3 act in the middle (at x = 0) rising to 1e5 N/m3 at the mirror edge.

Thermoviscous Acoustics, Frequency Domain (ta)

In the window, under click .

In the window for , locate the section.

Click

Select Domain 3 only.

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Pressure Acoustics, Frequency Domain (acpr)

In the window, under click .

Select Domains 1, 2, 4, and 5 only.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Spherical Wave Radiation 1

In the toolbar, click

and choose .

Select Boundaries 4, 5, 12, and 17 only.

Proceed to set up the Multiphysics couplings that couple Thermoviscous Acoustics to Pressure Acoustics and couple the interior shell to the surrounding thermoviscous acoustics domain.

Multiphysics

Acoustic–Thermoviscous Acoustic Boundary 1 (atb1)

In the toolbar, click

and choose > .

In the window for , locate the section.

From the list, choose .

When you select COMSOL automatically detects where a thermoviscous acoustics domain and a pressure acoustics domain meet and applies the coupling there.

Thermoviscous Acoustic–Structure Boundary 1 (tsb1)

In the toolbar, click

and choose > .

In the window for , locate the section.

From the list, choose .

Selecting here will result in the same setup as using the Shell selection.

Mesh 1

Free Triangular 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.8*dvisc.

To get a fully resolved acoustic boundary layer, it is also necessary to add a boundary layer mesh. Do this as the last mesh step.

Select the checkbox. In the associated text field, type 0.1*dvisc.

Select the checkbox. In the associated text field, type 1.4.

Select the checkbox. In the associated text field, type 0.4.

Size 2

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Points 4, 7, 15, and 18 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.2*dvisc.

Size 3

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Edges 8, 10–13, and 27–30 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.5*dvisc.

Click

The surface mesh of the mirror should look like this. Notice the higher mesh resolution regions. As seen in results, these are the regions with higher stress.

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 343[m/s]/f0/6.

In the text field, type 2.

Free Tetrahedral 1

In the toolbar, click

Boundary Layers 1

In the toolbar, click

Boundary Layer Properties

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 3.

From the list, choose .

In the text field, type 0.3*dvisc.

Click

The geometry with mesh should look like the figure below.

Proceed with the eigenfrequency study. The determined eigenfrequency including the complex part yields information about the attenuation in the system.

Study 1

Step 1: Eigenfrequency

In the window, under click .

In the window for , locate the section.

In the text field, type 10000.

From the list, choose .

Select the checkbox. In the associated text field, type 3.

These selections will search for the first three eigenmodes above 10 kHz. It is known that the first eigenfrequency is at about 10500 Hz (this can also be determined using a frequency study looking for the response peaks). The eigenfrequency search condition will make the solution converge faster. Selecting the will search for modes with both smaller and larger eigenfrequencies.

Again, turn off the generation of the default plots and generate and show the default solver to take a look at the solver suggestions.

In the window, click .

In the window for , locate the section.

Clear the checkbox.

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node.

Select and enable the first suggested iterative solver. In the eigenfrequency study for this size of problem, the direct solver is faster but it is more memory consuming.

In the window, expand the > > > node.

Right-click > > > > and choose .

In the toolbar, click

Result Templates

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Results

Mode Shape (shell)

Click the

button in the toolbar.

The default (first) selected eigenvalue is the vibrating torque mode at about 10,500 Hz. This is the vibration mode that will be studied next in the frequency domain study.

The first mode should look like the one in Figure 2 top.

The next eigenmode represents a symmetric vibration in the transverse direction happening at about 13100 Hz.

In the window for , locate the section.

From the list, choose .

In the toolbar, click

The second mode should look like the one in Figure 2 bottom.

The next mode happens at nearly 40 kHz and can be seen by selecting the last of the three computed eigenfrequencies.

Result Templates

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Results

Modal Stress (shell)

This second plot shows the von Mises stress distribution in the shell. To further improve the precision of the computed eigenfrequency, the mesh can be refined in the high-stress regions. Note that refining the mesh will rapidly lead to an increased computational time and memory consumption when solving.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Frequency Domain

In the window for , locate the section.

In the text field, type range(f0-deltaf,50,f0+deltaf).

Turn off the generation of default plots, only create the ones you need or use the . Plot the acoustic variations in temperature, pressure, and velocity as well as response curves. If turned on, the default plots for each physics interface will be generated.

In the window, click .

In the window for , locate the section.

Clear the checkbox.

Now, generate and show the default solver to take a look at the solver suggestions automatically generated.

Solution 2 (sol2)

In the toolbar, click

In the window, expand the node.

In this model, which couples pressure acoustics, thermoviscous acoustics, and shell, the default solver is a direct solver. Select and enable the first iterative suggestion. It is faster and more memory efficient than the direct solver.

In the window, expand the > > > node.

Right-click > > > > and choose .

In the toolbar, click

Results

In the window, expand the node.

Mirror 3D 1

In the window, expand the > node.

Right-click > and choose > .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Displacement

In the toolbar, click

In the window for , type Displacement in the text field.

Locate the section. From the list, choose .

Surface 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Deformation 1

Right-click and choose .

Displacement

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

In the toolbar, click

The plot should look like the one in Figure 4.

Temperature Variation

In the toolbar, click

In the window for , type Temperature Variation in the text field.

Locate the section. From the list, choose .

Locate the section. Select the checkbox.

Slice 1

Right-click and choose .

In the window for , locate the section.

In the text field, type ta.T_t.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

In the toolbar, click

The plot should look like the one in Figure 5 top.

Acoustic Pressure

In the toolbar, click

In the window for , type Acoustic Pressure in the text field.

Locate the section. From the list, choose .

Locate the section. Select the checkbox.

Isosurface 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > .

Locate the section. In the text field, type 20.

Locate the section. From the list, choose .

From the list, choose .

In the toolbar, click

The plot should look like the one in Figure 5 bottom.

Acoustic Velocity

In the window, right-click and choose .

In the window for , type Acoustic Velocity in the text field.

Locate the section. From the list, choose .

Locate the section. Select the checkbox.

Slice 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type ta.v_inst.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

In the toolbar, click

The plot should look like the one in Figure 6.

Proceed to setting up two 1D curves that represent the frequency response of the system by probing the displacement at the symmetry plane of the micromirror. The first plot gives the absolute value of the displacement. The second plot depicts the power (proportional to the velocity squared) on the dB scale and compared to a simple harmonic oscillator model tuned using parameters extracted from an eigenfrequency study of the system (instructions follow below). The parameters are given in the parameters list. Both curves are normalized by the maximum value.

Displacement Response

In the toolbar, click

In the window for , type Displacement Response in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

Point Graph 1

Right-click and choose .

Select Point 16 only.

In the window for , locate the section.

In the text field, type abs(w).

In the toolbar, click

The plot should look like the one in Figure 7 top. Proceed and create a plot to compare the simple resonator expression with the actual response of the system.

Power Response

In the toolbar, click

In the window for , type Power Response in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

Locate the section. From the list, choose .

Point Graph 1

Right-click and choose .

Select Point 16 only.

In the window for , locate the section.

In the text field, type 20*log10(abs(w*ta.iomega))-with(12,20*log10(abs(w*ta.iomega))).

Click to expand the section. Select the checkbox.

From the list, choose .

In the table, enter the following settings:


Legends
Simulated response

Global 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Expression	Description
P_simple	Simple harmonic oscillator model fit

Graph Marker 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 5.

Select the checkbox.

In the toolbar, click

The plot should look like the one in Figure 7 bottom. Finally, add an evaluation group to evaluate the harmonic oscillator parameters used in the simple analytical model.

Evaluation Group 1 - Harmonic Oscillator Parameters

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In the text field, type Evaluation Group 1 - Harmonic Oscillator Parameters.

Locate the section. Select the checkbox.

Global Evaluation 1

Right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Expression	Unit	Description
f_r	Hz	Resonance frequency
omega_r	Hz	Resonance angular frequency
alpha_r	Hz	Attenuation coefficient
Q_r	1	Q factor
df_r	Hz	Resonance half power width