Multiphysics field analysis and evolutionary optimization: Design of an electro-thermo-elastic microactuator

Abstract

In the paper, the optimization of an electro-thermo-elastic microactuator is proposed. In particular, the maximum temperature of the actuator is to be minimized, while the total displacement is to be maximized. For solving this problem, the Adaptive Gaussian Process-Assisted Differential Evolution AGDEMO method is applied.

Keywords

Electro-thermo-elastic microactuator multiphysics finite-element analysis efficient multiobjective design optimization

1. Introduction

The never-ending advance of technology at both the micrometer and the nanometer scales asks for the development of powerful and flexible modelling tools in order to help the design process of integrated devices and systems. In fact, design equations based on lumped-parameter models, like circuit models, are only able to size the device taking its main effect into account. However, a more accurate design cannot omit secondary effects like 3D effects and non-linear material properties; this generally implies the use of distributed-parameter models like field models. This is fortiori true when multiple physical domains coexist in the same device, giving rise to non-linear 3D coupled field models in three dimensions. The issue of computational cost, traditionally opposing lumped parameters against distributed parameters, is nowadays somehow obsolete, considering the broad availability of low-cost solid-state memories as well as multi-core processors. Moreover, real-life design models often incorporate various design criteria in mutual conflict which should be simultaneously minimized. Moving from this background, a contribution to multiobjective shape design based on coupled-field models is proposed here.

In the literature, interesting references can be found. For instance, in [1] the design optimization of an electromagnetic valve actuator is proposed, and a suitable combination of three design criteria is exploited. In turn, in [2] a shape optimization of comb drive electrostatic actuators, in order to achieve prescribed driving force profiles, is proposed. In [3] the magnetic field in a permanent-magnet spherical motor at no-load is recovered, after inverting the magnetic induction measured along an accessible surface; the final aim is to compute the on-load torque by means of the Lorentz’s law.

For solving the analysis problem, a multiphysics field analysis becomes more and more important for modelling the behaviour of an industrial device, like in [4, 5], or a MEMS device like [6]. Commercial software for solving this kind of problems are available, like Comsol Multiphysics [7].

In the last few years, automated optimization procedures for MEMS devices, based on stochastic optimization algorithms, have been proposed. In particular, the optimal design of a class of MEMS has been solved successfully with both evolutionary algorithms like in [8, 9] and cooperative algorithms like the biogeography-based optimization algorithm as shown in [10, 11] or the wind-driven optimization as shown in [12].

In the paper, evolutionary algorithms are linked with multiphysics finite-element analysis for approximating the Pareto front underlying the design problem: the case study is the optimal shape design of an electro-thermo-elastic microactuator.

The focused electro-thermo-elastic microactuator can be used in various systems, each of which may have different specifications on maximum temperature and displacement. Therefore, understanding the optimal trade-off relation of the microactuator design is important to support bespoke design for system integration. A natural way to achieve this is to approximate the Pareto front of the microactuator by multiobjective optimization methods. However, this is not trivial.

A major challenge is that the performance analysis of candidate designs has to be done by numerical simulations, which are computationally expensive. In particular, the FEA model built by COMSOL Multiphysics costs from 10 minutes to more than 1 hour per simulation for the targeted microactuator. Embedding the FEA model to existing multiobjective optimization methods may cost prohibitive optimization time. A routine method to address this issue is to employ surrogate modelling. A surrogate model is a computationally cheap mathematical approximation model aiming to estimate the output of the numerical simulation, which is often constructed by statistical learning techniques [13]. When collaborating surrogate modelling with the optimization, some computationally expensive numerical simulations can be replaced by cheap surrogate model predictions, the optimization time is therefore reduced.

There are two kinds of available surrogate-based optimization methods trying to solve this problem: (1) surrogate model-assisted multiobjective evolutionary algorithms (SA-MOEAs) [14, 15], (2) off-line surrogate model-based methods, which have been applied to MEMS design exploration [16].

Unfortunately, both of the above methods are difficult to solve the targeted problem. SA-MOEA is an emerging area in the computational intelligence field and to the best of our knowledge, there is no matured method although there is exciting progress. Most available SA-MOEAs can be classified into: (1) Methods simulating “Pareto-optimal” solutions predicted by the current surrogate model, although its quality may not be good enough [14, 15]; they are efficient but the generated solutions still have some distance to the true Pareto-front. (2) Methods using standard EAs for certain generations so as to build a good quality global surrogate model, and then starting to use surrogate models [17]; they can obtain near true Pareto-front solutions but the computing overhead is still not affordable for the targeted microactuator due to the generations using standard EAs. In addition, the targeted microactuator design exploration involves constraints, which are seldom involved by available SA-MOEAs. Off-line surrogate model-based methods first build a high-quality surrogate model using one-shot sampling. In the optimization, the numerical simulations will (mostly) be replaced by surrogate model prediction. Clearly, to obtain an accurate surrogate model, a sufficient number of samples are needed, which is determined by the volume of design space. In our microactuator, although there are 4 design variables (see Fig. 1), the ranges (L, dw) are quite large. A reasonable sampling (e.g., 2 $\mu$ m grid for all design variables, 5 $\mu$ m grid for L and 2 $\mu$ m grid for other design variables) needs around 5,000 to 10,000 simulations, which are not affordable.

Figure 1.

Geometry of the microactuator; the design variables of the inverse problem are also shown.

On the other hand, there are several high-performance surrogate model-assisted single objective evolutionary algorithms (SA-SOEAs) [18, 19]. A new framework, called surrogate model-aware evolutionary search framework is proposed in [19] and comparisons show up to 8 times speed improvement with comparable or better solution qualities with several popular SA-SOEAs based on more than 10 mathematical benchmark problems. The central idea of that framework is to improve the locations of training data points in order to construct high-quality surrogate models using fewer training data points (i.e. fewer simulations) [19, 20, 21]. Its application in MEMS design problems leads to more than 10 times speed improvement compared to standard EAs (wall clock time) [19]. Taking advantage of employing AGDEMO, a new framework is proposed to address the targeted microactuator.

2. The electro-thermo-elastic microactuator

In the early times of MEMS technology, the main actuation principle relied on electrostatic field [2], which is still moderately used. More recently, however, the techniques of actuation most performing in terms of position control fall in three categories: piezoelectric, thermal or magnetic actuation, respectively [3, 8]. In the paper, reference is made to a thermally actuated device because it is a clear example of a multiphysics domain which asks for a coupled-field model.

In particular, a electro-thermo-elastic microactuator is considered. This kind of microactuators can be applied in many different fields e.g. in civil engineering for monitoring the building deterioration and movements due to earthquakes, in electrical and electronics engineering for microrobotics, microrelays, microfluidics as well as in biomedical engineering for robotic surgery and miniature medical instrumentation.

Because these fields of application need very well performing devices, there is the need to optimize such kind of devices. In particular, the aspects to be improved have to be chosen, depending on the application. Usually, in the design of MEMS circuit models or lumped-parameter models of devices are preferred to field models. This is certainly convenient in view of a preliminary design, like in the prototyping stage; however, design optimization is recommended to be based on field or distributed-parameter models, in order to take into account secondary effects in the device behavior. In this area of industrial products, searching for a marginal improvement of the design of a class of devices is an effective way to innovate, which brings also the benefit of reducing the cost of prototyping.

Figure 2.

Current, temperature (left) and stress (right) distribution in the microactuator.

The device under study, shown in Fig. 1, has length L and width hh. An electric voltage is applied between two electrodes A and B; therefore, an electric current I flows in two out of three arms of the device (hot arms), while the third arm is current free (cold arm), as shown in Fig. 2.

The actuator is fixed to a substrate at the three arm ends, while three cylindrical bushings act on the cold arm, in order to make it rotate in the xy-plane (Fig. 2). The deformation due to the overheating of the hot arms with respect to the cold one is responsible for the rotation of the actuator.

Figure 3.

Mesh of the model, left, and a detail of it, right.

3. Analysis problem

A parametric finite-element model of the device has been developed; a typical mesh used, composed of about 8,000 three-dimensional elements, is shown in Fig. 3.

Table 1
Material properties of the finite element model

Material property	Value
Electric conductivity $\sigma$ [Sm ${}^{-1}$ ]	5 10 ${}^{4}$
Thermal conductivity k [Wm ${}^{-1}$ K ${}^{-1}$ ]	34
Thermal expansion coefficient $\alpha$ [K ${}^{-1}$ ]	2.6 10 ${}^{-6}$
Young modulus [Pa]	210 10 ${}^{9}$
Fracture toughness [Pa]	1.44–2.51 10 ${}^{9}$
Poisson coefficient	0.22

The material properties used in the model are listed in Table 1.

The following equations, which are coupled at the right-hand side level, are subsequently solved:

$\displaystyle{\rm{\bf J}}=\sigma{\rm{\bf E}}=-\sigma\nabla V$ (1) $\displaystyle Q={\bf J}\cdot{\bf E}=-\nabla\cdot(k\nabla T)$ (2) $\displaystyle\alpha(T-T_{\text{ref}})=\varepsilon^{Th}$ (3)

where ${{\bf J}}$ is the current density, ${{\bf E}}$ the electric field, Q the heat due to Joule-effect, T the temperature and h the convection coefficient which is different in the upper and in the lower surface, ${{\bm{\varepsilon}}}^{Th}$ the thermal part of the total strain.

The following boundary conditions are applied.

Conduction current problem

$\displaystyle V_{A}=\text{const},V_{B}=0$ (4) $\displaystyle{{\bf n}}\cdot{{\bf J}}=0$ (5)

elsewhere.

Thermal problem

$\displaystyle T=293.15∼{}K$ (6)

at the simple supports and at the hinges

$-{{\bm{n}}}\cdot(k\nabla T)=h(T_{ext}-T)$ (7)

elsewhere.

Mechanical problem

${{\bf u}}\cdot{\rm{\bf n}}=0$ (8)

at the simple supports,

${{\bf u}}=0$ (9)

at the hinges.

The three problems are coupled via the thermal heat Q (electric and thermal problems) and via the temperature T (thermal and elastic problems). Because non-linearities of material parameters are not taken into account, it is possible to solve the three problems subsequently. Therefore, a weakly-coupled analysis problem is dealt with; to solve it, a cascade algorithm like the successive substitution algorithm can be applied.

In general, in this kind of muliphysics problems the main sources of non-linearity depend on the electric and thermal properties of materials against temperature; in particular both electrical conductivity $\sigma$ and thermal conductivity k are temperature dependent [4, 5].

If the non-linearities are taken into account the problem to be solved is fully-coupled and this means a substantial increase of computational time. For the sake of an example, if the electrical conductivity is modelled as a temperature dependent function, the electrical and thermal problems have to be solved according to a nested loop: in fact, the electrical problem will give rise to a heating of the device, which, in turn, will modify the electrical conductivity, eventually changing the electrical solution and consequently the temperature and so on.

For the considered device, because the electrical power is low, approximately 10 mW, the variation range of the temperature is few hundreds K, so the electrical and thermal conductivities of the silicon have been assumed to be constant, referring to the average operating temperature. Also for the elastic problem, thanks to the assumption of small displacements, it is reasonable that the material properties are considered to be linear.

Under the hypothesis of neglecting self-weight and inertial forces, the continuum mechanics subproblem consists in finding the current second Piola-Kirchhoff stress tensor ${{\bf S}}$ , the Green-Lagrange strain tensor ${{\bm{\varepsilon}}}$ and the displacement vector ${{\bf u}}$ under the action of a given temperature field, say ${{\bf T}}$ , that translates into an assigned inelastic strain ${{\bm{\varepsilon}}}^{T}$ . An isotropic hyperelastic material of Mooney-Rivlin type is considered. If now ${{\bf X}}$ and ${{\bf x}}$ are respectively used to indicate the initial and current location of a given particle so that ${{\bf x}}={{\bf x}}({{\bf X}},t)$ gives the location at time $t$ of a particle with material coordinates ${{\bf X}}$ , upon introduction of the deformation gradient

${{\bm{F}}}=\frac{\partial{{\bf x}}}{\partial{{\bf X}}}=\nabla{{\bm{u}}}+{{\bf I}}$ (10)

and of the right Cauchy-Green strain tensor ${{\bf C}}$ that reads

${{\bm{C}}}={{\bm{F}}}^{T}{{\bm{F}}}$ (11)

one may define the Green-Lagrange strain tensor that is dual in the virtual work sense of the second Piola-Kirchhoff stress ${{\bf S}}$ and reads

${{\bm{\varepsilon}}}=\frac{1}{2}\left({{{\bf C}}-{{\bf I}}}\right)$ (12)

that plays a crucial role within the Comsol programming environment that has been used for the purpose of finding the optimal solution. The problem is in fact solved resorting to a Total Lagrangian Formulation that is briefly recalled next for completeness sake following the classical contributions [23, 24]. The basic equation to be solved writes

$\int\limits_{{}^{0}V}{{}^{t+}{}_{0}{}^{\Delta t}{{\bf S}}_{ij}{{\bm{\delta}}}^% {t+}{}_{0}{}^{\Delta t}\varepsilon_{ij}{{\bf d}}^{0}{{\bf V}}}={}^{t+\Delta t}% {{\bf R}}$ (13)

where the external virtual work ${}^{t+\Delta t}{{\bf R}}$ may be written as

${}^{t+\Delta t}{\bf R}=\int\limits_{{}^{t+\Delta t}{\bf V}}{{}^{t+\Delta t}}{% \bf f}_{i}^{B}{{\bm{\delta}}{\bf u}}_{i}{\bf d}^{t+\Delta t}{\bf V}+\int% \limits_{{}^{t+\Delta t}{\bf S}}{{}^{t+\Delta t}{\bf f}}_{i}^{S}{{\bm{\delta}}% {\bf u}}_{i}{\bf d}^{t+\Delta t}{\bf S}$ (14)

and

${}^{t+\Delta t}{{\bf f}}_{i}^{B}$ are the external volume forces at time $t+\Delta t$ ;

${}^{t+\Delta t}{{\bf f}}_{i}^{S}$ are the external surface tractions at time $t+\Delta t$ ;

${}^{t+\Delta t}{{\bf V}}_{f}$ is the volume at time $t+\Delta t$ ;

${}^{t+\Delta t}{{\bf S}}_{f}$ is the surface at time $t+\Delta t$ ;

$\delta{{\bf u}}_{i}$ are the virtual displacements at time $t+\Delta t$ .

One should notice that the time variable t is merely used to order events whereas inertial effects are not accounted for as already mentioned. After some algebra, the linearized governing equations may be shown to be

$\int\limits_{{}^{0}{\rm{\bf V}}}{}_{0}{\rm{\bf C}}_{ijrs}{\rm{\bf e}}_{rs}{\bm% {\delta}}{\bf e}_{ij}{{\bf d}}^{0}{{\bf V}}+\int\limits_{{}^{0}v}{}_{0}^{t}{{% \bf S}}_{ij}{{\bm{\delta}}}_{0}{\bm{\eta}}_{ij}{{\bf d}}^{0}{\bf V}={}^{t+% \Delta t}{\bf R}-\int\limits_{{}^{0}v}{}_{0}^{t}{{\bf S}}_{ij}{{\bm{\delta}}}_% {0}{{\bf e}}_{ij}{{\bf d}}^{0}{{\bf V}}$ (15)

where the linear and nonlinear incremental strains have been respectively defined as

$\displaystyle{}_{0}{\rm{\bf e}}_{ij}=\frac{1}{2}\left({{}_{0}{\rm{\bf u}}_{i,j% }+{}_{0}{\rm{\bf u}}_{j,i}+{}_{0}^{t}{\rm{\bf u}}_{k,i}{\rm{\bf u}}_{k,j}+{}_{% 0}{\rm{\bf u}}_{k,i}{}_{0}^{t}{\rm{\bf u}}_{k,j}}\right)$ (16) $\displaystyle{}_{0}{{\bm{\eta}}}_{ij}=\frac{1}{2}{}_{0}{\rm{\bf u}}_{k,i}{}_{0% }{\rm{\bf u}}_{k,j}$

It is important to note that Eq. (15) is linear in the incremental displacements since ${{\bm{\delta}}}_{0}{{\bf e}}_{ij}$ is independent of ${{\bf u}}_{i}$ .

For solving the whole forward problem it takes about 5 minutes on an Intel i7, 3.6 GHz, equipped with 16 GB of RAM.

4. Design exploration of the microactuator

The inverse problem reads as follows: acting on the design variables (see Fig. 1)

L length of the actuator

hh thickness of the actuator

dw width of the cold arm

d width of the hot arms

Find the minimum of

$f=\left[{\beta\frac{T_{\max}}{\gamma}-\left({1-\beta}\right)\frac{\left|u% \right|}{\delta}}\right]_{V=V_{\min}}+\left[{\beta\frac{T_{\max}}{\gamma}-% \left({1-\beta}\right)\frac{\left|u\right|}{\delta}}\right]_{V=V_{\max}}$ (17)

with $V_{\min}=$ 1 V, $V_{\max}=$ 5 V, $\gamma=$ 900 K, $\delta=$ 1.9 10 ${}^{-6}$ m, 0 $<\beta<$ 1, subject to the constraints:

Geometrical congruency (dw $>$ 2d)

$T_{\max}<$ 1500 K

Stress $<$ 1.44 GPa

The boundaries of the design variables are L [56–300] $\mu$ m, hh [2–5] $\mu$ m, dw [7–30] $\mu$ m, d [1–7] $\mu$ m.

The objective function Eq. (17) allows to minimize the maximum temperature T ${}_{\max}$ , normalized by $\gamma$ , and, at the same time, maximize the displacement u, normalized by $\delta$ . The minimization of the temperature is important for the device performance: the lower the overheating, the lower the material deterioration. On the other hand, a sizable displacement is needed, to guarantee a good operation of the device.

These two requirements are in conflict, because a high temperature allows for a high gradient of temperature and hence large displacement, on the contrary, if the temperature of the hot arm is low, a small displacement arises.

Figure 4.

Bespoke design exploration framework for the electro-thermo-elastic microactuator.

For solving the inverse problem, the new framework proposed, by applying the AGDEMO method, is described in Fig. 4.

A few clarifications are as follows: (1) Traditional multiobjective design exploration follows the idea of generating a complete Pareto front and directly selecting designs from it in system integration. This indicates a number of weight vectors, which largely increases the computational cost. In our framework, a few weight vectors are used to sparsely reflect the general shape of the Pareto front, providing an approximate understanding of the trade-off relation between Tmax and u. Based on this information, the designer is able to estimate suitable specifications (e.g. $\text{T}_{\max}\in\left[{400,405}\right]\text{K}$ , $\text{Stress}\in\left[{0.22,0.24}\right]\text{GPa}$ ) in system integration. AGDEMO will be run again to obtain the bespoke design. As above said, AGDEMO is efficient. (2) Because the number of weight vectors is small, there is ahigher risk of missing some parts of the Pareto front since neighbouring weight vectors are difficult to compensate each other. A normalization is thus carried out.

4.1 The AGDEMO method

This subsection provides a brief description of the AGDEMO method. More details are in [22]. The AGDEMO method is a Gaussian Process (GP) [24] surrogate model assisted evolutionary approach for MEMS optimization problems withfield-dependent 3D analysis. It has the following performances: (1) Achieves comparable results with MEMS optimization methods which directly embed numerical simulations to a standard EA; (2) More than an order of speed improvement compared to standard EAs and off-line surrogate model-based methods for MEMS design; (3) General enough for various kinds of MEMS. AGDEMO works as shown in Fig. 5.

Figure 5.

The AGDEMO Method [19].

In terms of constraint handling, AGDEMO uses the penalty function method [28]. The penalized function is often a piecewise function due to the constraints, but a reasonably continuous and smooth hypersurface is important for generating high-quality surrogate models. To address this problem, a separate GP model isconstructed for each performance, which is not piecewise. The penalized function value is then calculated by the predicted performance values. Experimental results verified this constraint handling method [19].

Following the parameter setting rules of [22], we use the following parameters for our implementation: $\alpha=$ 30, $\lambda=$ 30, $\tau=$ 8 $\times$ d (d is the number of design variables, 4). The penalty coefficient is set to 50. The high optimization quality and large efficiency improvement come from the surrogate model-aware evolutionary search framework [19], which has been verified by intensive empirical tests. We also verified it in our pilot experiments and such comparisons will not be repeated here.

5. Results

Different optimizations are carried out by varying the value of $\beta$ , i.e. $\beta=$ 0, $\beta=$ 0.25, $\beta=$ 0.5, $\beta=$ 0.75 and $\beta=$ 1. The results are shown in Table 2.

Table 2
Design variables of optimal devices and corresponding objective function for AGDEMO method

$\beta$	L [ $\mu m$ ]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]	o.f.	o.f. calls
0	299.73	4.99	29.24	2.33	$-$ 1.18	135
0.25	298.83	4.99	29.93	2.18	$-$ 0.67	123
0.5	299.82	4.83	29.45	1.94	$-$ 0.17	158
0.75	299.81	4.67	29.77	1.39	0.32	132
1	299.82	2.00	14.69	1.01	0.74	130

Figure 6.

History of the objective function for different values of $\beta$ .

In Fig. 6 the history of the objective function obtained for different values of $\beta$ re shown.

The constraints are fulfilled as shown in Table 4.

Table 3

Maximum value of the temperature and stress for each optimal point

$\beta$	${T}_{\max}$ [K]	Stress [GPa]
0	470.25	0.36
0.25	463.00	0.32
0.5	446.83	0.44
0.75	410.32	0.25
1	372.60	0.20

Table 4

Design variables of optimal devices and corresponding objective function, new objective function (20)

$\beta$	L	hh	dw	d	o.f.	o.f.
	[ $\mu$ m]	[ $\mu$ m]	[ $\mu$ m]	[ $\mu$ m]		calls
0	299.4	4.70	29.81	2.05	$-$ 4	159
0.25	298.8	4.97	29.53	1.99	$-$ 1.15	152
0.5	299.5	4.95	29.22	1.36	1.61	151
0.75	299.8	2.06	29.26	1.00	4.17	200
1	299.44	2.16	15.44	1.00	6.60	181

By defining the objective function space defined by functions $f_{1}$ and $f_{2}$ ,

$\displaystyle f_{1}=\left[{\left|u\right|}\right]_{V=V\min}+\left[{\left|u% \right|}\right]_{V=V\max}$ (18) $\displaystyle f_{2}=\left[{T_{\max}}\right]_{V=V\min}+\left[{T_{\max}}\right]_% {V=V\max}$ (19)

Figure 7 is obtained, in which the optimization results are shown; a random sampling of this objective space is added.

Table 5

Maximum value of the temperature and stress for each optimal point

$\beta$	T ${}_{\max}$ [K]	Stress [GPa]
0	453	0.37
0.25	451.6	0.21
0.5	409.2	0.22
0.75	372.5	0.26
1	373.59	0.19

Table 6

Objective function 1, s parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	$-$ 7.51 10 ${}^{-3}$	3.77 10 ${}^{-2}$	$-$ 6.05 10 ${}^{-4}$	5.32 10 ${}^{-3}$
0.25	3.91 10 ${}^{-3}$	3.48 10 ${}^{-2}$	5.04 10 ${}^{-3}$	3.95 10 ${}^{-2}$
0.5	1.61 10 ${}^{-3}$	4.09 10 ${}^{-2}$	4.99 10 ${}^{-3}$	1.06 10 ${}^{-1}$
0.75	2.11 10 ${}^{-3}$	2.60 10 ${}^{-2}$	1.84 10 ${}^{-3}$	4.06 10 ${}^{-1}$
1	6.12 10 ${}^{-4}$	1.20 10 ${}^{-2}$	9.81 10 ${}^{-3}$	5.57 10 ${}^{-1}$

Figure 7.

Objective space. Results of the five optimization runs are shown: $\beta=$ 0 (circle), $\beta=$ 0.25 (star), $\beta=$ 0.5 (diamond), $\beta=$ 0.75 (square), $\beta=$ 1 (triangle).

The temperature and strain field maps of the optimal devices found with $\beta=$ 0 and $\beta=$ 1 (see Fig. 7) are shown in Figs 8 and 9.

5.1 Results with normalized objective function

The results in Fig. 7 show that the AGDEMO method finds good solutions, seemingly, however, they tend to form a cluster. A new optimization is carried out considering the following function:

$f=\left[{\beta\frac{T_{\max}}{\left|{R_{2}-U_{2}}\right|}-\left({1\!-\!\beta}% \right)\frac{\left|u\right|}{\left|{R_{1}-U_{1}}\right|}}\right]_{V=V_{\min}}% \!\!+\!\!\left[{\beta\frac{T_{\max}}{\left|{R_{2}-U_{2}}\right|}\!-\!\left({1-% \beta}\right)\frac{\left|u\right|}{\left|{R_{1}-U_{1}}\right|}}\right]_{V=V_{% \max}}$ (20)

where ( $U_{1}$ , $U_{2}$ ) is the Utopia point, while ( $R_{1}$ , $R_{2}$ ) is the Nadir point calculated by considering the approximation of the Pareto front found by AGDEMO. In particular, $U_{1}=$ 2.24 10 ${}^{-6}$ m, $U_{2}=$ 669 K, $R_{1}=$ 1.68 10 ${}^{-6}$ m, $R_{2}=$ 770 K.

Table 7

Objective function 2, s parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	$-$ 1.22 10 ${}^{6}$	6.18 10 ${}^{6}$	3.21 10 ${}^{-4}$	5.74 10 ${}^{7}$
0.25	$-$ 1.15 10 ${}^{6}$	5.66 10 ${}^{6}$	2.97 10 ${}^{-4}$	5.99 10 ${}^{7}$
0.5	$-$ 1.06 10 ${}^{6}$	5.00 10 ${}^{6}$	2.68 10 ${}^{-4}$	6.28 10 ${}^{7}$
0.75	$-$ 8.32 10 ${}^{5}$	3.17 10 ${}^{6}$	1.79 10 ${}^{-4}$	7.23 10 ${}^{7}$
1	$-$ 5.87 10 ${}^{5}$	8.19 10 ${}^{6}$	4.90 10 ${}^{-4}$	6.34 10 ${}^{7}$

Table 8

Objective function weighted sum, s parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	3.95 10 ${}^{3}$	$-$ 1.99 10 ${}^{4}$	3.18 10 ${}^{2}$	$-$ 2.8 10 ${}^{3}$
0.25	$-$ 1.86 10 ${}^{3}$	$-$ 1.22 10 ${}^{4}$	$-$ 1.96 10 ${}^{3}$	1.04 10 ${}^{3}$
0.5	$-$ 1.01 10 ${}^{3}$	$-$ 7.98 10 ${}^{3}$	$-$ 1.30 10 ${}^{3}$	6.86 10 ${}^{3}$
0.75	$-$ 9.71 10 ${}^{2}$	$-$ 7.80 10 ${}^{2}$	$-$ 2.51 10 ${}^{2}$	6.74 10 ${}^{3}$
1	$-$ 6.52 10 ${}^{2}$	9.10 10 ${}^{3}$	2.59 10 ${}^{1}$	7.05 10 ${}^{4}$

Figure 8.

Temperature [K] (up) and stress [Nm ${}^{-2}$ ] (bottom) field maps for $\beta=$ 0.

The results are shown in Tables 4, 6 and Fig. 10.

The results obtained by the previous optimizations show that a long device is preferred. In fact, all the optimal microactuators are characterized by the maximum value of the length L (300 $\mu$ m).

The width of the hot arms and the thickness of the actuator are related to the electrical resistance between hinges A and B: decreasing hh and d means increasing the electrical resistance. Consequently, the current flowing in the hot arms decreases and the Joule losses decreases too. In fact, the points in Fig. 10 with lower temperatures (green and black triangles and green square) are characterized by the minimum values of hh and d (2 $\mu$ m and 1 $\mu$ m, respectively).

5.2 Sensitivity analysis

In view of a robust design, a sensitivity analysis has been carried out.

Table 9
Objective function 1, $\Delta$ f parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	$-$ 1.34 10 ${}^{-4}$	6.74 10 ${}^{-4}$	$-$ 1.08 10 ${}^{-5}$	9.512 10 ${}^{-5}$
0.25	7.01 10 ${}^{-5}$	6.23 10 ${}^{-4}$	9.02 10 ${}^{-5}$	7.07 10 ${}^{-4}$
0.5	2.90 10 ${}^{-5}$	7.39 10 ${}^{-4}$	9.01 10 ${}^{-5}$	1.92 10 ${}^{-3}$
0.75	4.05 10 ${}^{-5}$	5.00 10 ${}^{-4}$	3.54 10 ${}^{-5}$	7.81 10 ${}^{-3}$
1	1.46 10 ${}^{-5}$	2.86 10 ${}^{-3}$	1.17 10 ${}^{-4}$	1.32 10 ${}^{-2}$

Table 10

Objective function 2, $\Delta$ f parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	$-$ 6,34 10 ${}^{-5}$	3.21 10 ${}^{-4}$	2.22 10 ${}^{-6}$	2.98 10 ${}^{-3}$
0.25	$-$ 6,02 10 ${}^{-5}$	2.97 10 ${}^{-4}$	5.19 10 ${}^{-6}$	3.14 10 ${}^{-3}$
0.5	$-$ 5.68 10 ${}^{-5}$	2.68 10 ${}^{-4}$	9.00 10 ${}^{-7}$	3.36 10 ${}^{-3}$
0.75	$-$ 4.70 10 ${}^{-5}$	1.79 10 ${}^{-4}$	$-$ 6.09 10 ${}^{-7}$	4.08 10 ${}^{-3}$
1	$-$ 3.51 10 ${}^{-5}$	4.90 10 ${}^{-4}$	6.97 10 ${}^{-7}$	3.79 10 ${}^{-3}$

Table 11

Objective function weighted sum, $\Delta$ f parameter

$\beta$	L [ $\mu$ m]	hh [ $\mu$ m]	dw [ $\mu$ m]	d [ $\mu$ m]
0	$-$ 1.34 10 ${}^{-4}$	6.74 10 ${}^{-4}$	$-$ 1.08 10 ${}^{-5}$	9.51 10 ${}^{-5}$
0.25	$-$ 1.11 10 ${}^{-4}$	7.26 10 ${}^{-4}$	1.17 10 ${}^{-4}$	$-$ 6.22 10 ${}^{-5}$
0.5	2.40 10 ${}^{-4}$	1.90 10 ${}^{-3}$	3.10 10 ${}^{-4}$	$-$ 1.63 10 ${}^{-3}$
0.75	$-$ 1.23 10 ${}^{-4}$	$-$ 9.86 10 ${}^{-5}$	$-$ 3.17 10 ${}^{-5}$	8.53 10 ${}^{-4}$
1	$-$ 3.51 10 ${}^{-5}$	4.90 10 ${}^{-4}$	6.97 10 ${}^{-7}$	3.79 10 ${}^{-3}$

Figure 9.

Temperature [K] (up) and stress [Nm ${}^{-2}$ ] (bottom) field maps for $\beta=$ 1.

Figure 10.

Objective space. Results of the five optimization runs are shown: $\beta=$ 0 (circle), $\beta$ 0.25 (star), $\beta=$ 0.5 (diamond), $\beta=$ 0.75 (square), $\beta=$ 1 (triangle). In black the results obtained with the objective function (17), in green those obtained with the new objective function (20).

In particular, two sensitivity parameters have been considered for each objective function

$\displaystyle s=\frac{f^{+}-f^{-}}{2\delta}$ (21) $\displaystyle\Delta f=\frac{f_{i}^{+}-f_{i}^{-}}{f_{i}}$ (22)

where $\delta=$ 2 10 ${}^{-8}$ m is the perturbation step, $f^{+}$ and $f^{-}$ are the objective functions evaluated by varying each variable of the $\delta$ quantity ( $f^{+}$ is evaluated by adding $\delta$ and $f^{-}$ by subtracting it), and $f_{i}$ is the objective function value found by the optimization procedure.

The results about the s parameter are shown in Tables 6, 8 and 8, while those related to the $\Delta$ f parameter are shown in Tables 10, 10 and 11, for the three objective functions.

The Tables 6–11 could be very useful in view of a practical construction of this device. In particular, they show that, in view of a robust design, the most important design variables are hh and d, hence the thickness of the whole actuator and the width of the hot arms.

6. Conclusions and ongoing extensions

A multiphysics model coupled with an evolutionary algorithm for the multi-objective optimization of structures has been proposed with the aim of finding the optimal shape of an electro-thermo-elastic actuator. Coupling between electric, thermal and structural fields is assumed in such a way that a cascade of subproblems are set up and solved. The electric problem is solved first that originates a temperature field. The temperature field is in turn used to generate an imposed strain field that represents the forcing term of the continuum mechanics subproblem. A few multi-objective optimal problems are considered that approximate the Pareto optimal points with reference to two distinct objectives that depend on the displacements and on the temperature, the design variables being of geometric type and namely actuator thickness, length and width. Furthermore, a peculiar sensitivity analysis method has been proposed to ensure the robustness of the optimal design allowing to assess the design variables with respect to which the design happens to be most sensitive.

The application of more complex multiphysics models characterized by a stronger coupling as well as the inclusion of material nonlinearities are among the extensions that are currently under investigation and development. The achievement of such further goals would open the way to the last stage of the research that is expected to include the actual manufacture and test of a prototypical specimen.

Appendix

A.A Brief Description to Gaussian Process Surrogate Modeling.

To model an unknown function $y=f(x)$ , $x\in R^{d}$ , GP modeling assumes that $f(x)$ at any point $x$ is a Gaussian distributed stochastic variable with mean $\mu$ , and variance $\sigma$ , where $\mu$ and $\sigma$ are two constants independent of $x$ . For any $x$ , $f(x)$ is a sample of $\mu+\varepsilon(x)$ where $\varepsilon(x)\approx N\left({0,\sigma^{2}}\right)$ . By maximizing the likelihood function that $f(x)=y^{i}$ at $x=x^{i}$ , $i=1,\ldots,k$ (where $x^{1},\ldots,x^{k}\in R^{d}$ ) and their function values $y^{1},\ldots,y^{k}$ are $k$ training data points and best linear unbiased prediction:

$\hat{f}(x)=\hat{\mu}+r^{T}C^{-1}\left({y-1\hat{\mu}}\right)$ (23)

the mean squared error is:

$s^{2}(x)=\hat{\sigma}^{2}\left[{1-r^{T}C^{-1}r+\frac{\left({1-1^{T}C^{-1}r}% \right)^{2}}{1^{T}C^{-1}r}}\right]$ (24)

where $r=[{c({x,x^{1}}),\ldots,c({x,x^{k}})}]^{T}$ . C is a $K\times K$ matrix whose (i,j)-element is $c({x^{i},x^{j}})$ . $c({x^{i},x^{j}})$ is the correlation function between $x^{i}$ and $x^{j}$ . 1 is a K-dimensional column vector of ones.

A brief description to differential evolution algorithm

Differential Evolution is a popular population-based metaheuristic algorithm for continuous optimization. Suppose that P is a population and the best individual in Pis $x^{\textit{best}}$ . Let $x_{1},\ldots,x_{d}\in R^{d}$ be an individual solution in P. To generate a child solution $u_{1},\ldots,u_{d}$ for $x$ , DE works as follows.

A donor vector is first produced by mutation (DE/current-to-best/1):

$v_{i}=x^{i}+F\cdot\left({x^{best}-x^{i}}\right)+F\cdot\left({x^{r_{1}}-x^{r_{2% }}}\right)$ (25)

where $x^{r_{1}}$ and $x^{r_{2}}$ are two different solutions randomly selected from P and also different from $x^{\textit{best}}$ and $x^{i}$ . $F\in\left({0,2}\right]$ is the scaling factor. Then the following crossover operator is applied to produce the child $u$

Randomly select a variable index $j_{\textit{rand}}\in 1,\ldots,d$

For each $j=1,\ldots,d$ generate a uniformly distributed random number rand from (0,1) and set:

$\displaystyle u_{j}=\left\{{\begin{array}[]{ll}v_{j},&\textit{if}\ \left({% \textit{rand}\leqslant\textit{CR}}\right)|j=j_{\textit{rand}}\\ x_{j},&\text{otherwise}\\ \end{array}}\right.$ (26)

where $CR\in\left[{0,1}\right]$ is a constant called the crossover rate. More details can be found in [22].

References

Delinchant

Rakotoarison

H.L.

Ardon

Chabedec

and Cugat

, Gradient based optimization of semi-numerical models with symbolic sensitivity: Application to simple ferromagnetic MEMS switch device, IJAEM 30 (2009), 189–200.

Kotani

Yamada

Yamasaki

Ohkado

Izui

and Nishiwaki

, Driving force profile design in comb drive electrostatic actuators using a level set-based shape optimization method, Structural and Multidisciplinary Optimization 51(2) (2015), 369–383.

Grossard

Rotinat-Libersa

M.C.

Chaillet

and Boukallel

, Mechanical and control-oriented design of a monolithic piezoelectric microgripper using a new topological optimization method, IEEE Trans. Mechatronics 14(1) (2009), 32–45.

Di Barba

Dolezel

Karban

Kus

Mach

Mognaschi

M.E.

and Savini

, Multiphysics field analysis and multiobjective design optimization: a benchmark problem, Inverse Problems in Science and Engineering (IPSE) 22(7) (2014).

Di Barba

Dolezel

Mognaschi

M.E.

Savini

and Karban

, Non-linear multi-physics analysis and multi-objective optimization in electroheating applications, IEEE Transactions on Magnetics 50(2) (2014).

Aage

and Lazarov

B.S.

, Parallel framework for topology optimization using the method of moving asymptotes, Structural and Multidisciplinary Optimization 47(4) (2013), 493–505.

Comsol.com, Comsol’s Official Website. [online]. Available at https://www.comsol.com/ [Accessed 11st July 2016].

Di Barba

and Wiak

, Evolutionary Computing and Optimal Design of MEMS, IEEE Trans. on Mechatronics 20(4) (2015), 1660–1667.

Zhu

Wang

and Fan

, Evolutionary design optimization of MEMS: A brief review, In: IEEE International Conference on Industrial Technology (ICIT), (2010), 1683–1687.

10.

Di Barba

Dughiero

Mognaschi

M.E.

Savini

and Wiak

, Biogeography-Inspired Multiobjective Optimization and MEMS Design, IEEE Transactions on Magnetics 52(3) (2016).

11.

Di Barba

Mognaschi

M.E.

Savini

and Wiak

, Island biogeography as a paradigm for MEMS optimal design, IJAEM 51(s1) (2016), 97–105.

12.

Di Barba

, Multi-objective wind-driven optimisation and magnet design, Electronics Letters 52(14) (2016), 1216–1218.

13.

Koziel

and Leifsson

, Surrogate-based Modeling and Optimization: Applications in Engineering, Springer, 2013.

14.

Emmerich

Giannakoglou

and Naujoks

, Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels, IEEE Transactions on Evolutionary Computation 10(4) (2006), 421–439.

15.

Zhang

Liu

and Tsang

, Expensive multiobjective optimization by MOEA/D with Gaussian process model, IEEE Transactions on Evolutionary Computation 14(3) (2010), 456–474.

16.

Lee

Park

Niu

and Filipovic

, Design and optimisation of one-port rfmems resonators and related integrated circuits using ANN based macromodelling approach, IEE Proceedings-Circuits, Devices and Systems 153(5) (2006), 480–488.

17.

Lim

Jin

Ong

and Sendhoff

, Generalizing surrogate assistedevolutionary computation, IEEE Transactions on Evolutionary Computation 14(3) (2010), 329–355.

18.

Müller

Shoemaker

C.A.

and Piché

, SO-MI: A surrogate model algorithm for computationally expensive nonlinear mixed-integer black-box global optimization problems, Computers & Operations Research 40(5) (2013), 1383–1400.

19.

Liu

Zhang

and Gielen

, A gaussian process surrogate model assisted evolutionary algorithm for medium scale expensive optimization problems, IEEE Transactions on Evolutionary Computation 18(2) (2014), 180–192.

20.

Liu

Zhao

Reynaert

and Gielen

G.G.E.

, Gaspad: A general and efficient mm-wave integrated circuit synthesis method based on surrogate model assisted evolutionary algorithm, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 33 (2) (2014) 169–182.

21.

Liu

Chen

Zhang

Gielen

and Grout

, Behavioral study of the surrogate model-aware evolutionary search framework, in: IEEE Congress on Evolutionary Computation, 2014, pp. 715–722.

22.

Liu

and Nikolaeva

, Efficient global optimization of MEMS based on surrogate model assisted evolutionary algorithm. In 2016 Design, Automation & Test in Europe Conference & Exhibition (DATE), 555–558.

23.

Bathe

K.J.

Ramm

and Wilson

E.L.

, Finite Element Formulations for Large Deformation Dynamic Analysis, International Journal for Numerical Methods in Engineering 9 (1975), 353–386.

24.

Bathe

K.J.

, Finite Element Procedures, Prentice-Hall, 1996.

25.

Santner

Williams

and Notz

, The design and analysis of computer experiments. Springer, 2003.

26.

Stein

, Large sample properties of simulations using latin hypercube sampling, Technometrics (1987), 143–151.

27.

Storn

and Price

, Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces, Journal of Global Optimization 11(4) (1997), 341–359.

28.

Himmelblau

D.M.

, Applied nonlinear programming, McGraw-Hill Companies, 1972.

Multiphysics field analysis and evolutionary optimization: Design of an electro-thermo-elastic microactuator

Abstract

Keywords

1. Introduction

Table 1 Material properties of the finite element model

Table 2 Design variables of optimal devices and corresponding objective function for AGDEMO method

Table 9 Objective function 1, Δ f parameter

Appendix

A brief description to differential evolution algorithm

References

Table 1
Material properties of the finite element model

Table 2
Design variables of optimal devices and corresponding objective function for AGDEMO method

Table 9
Objective function 1, $\Delta$ f parameter