The In-plane Expansion of Fractured Thermally Pre-stressed Glass Panes

An equivalent temperature difference model for engineering glass design

Jens H. Nielsen Michael A. Kraus Jens Schneider

Assoc. Prof. Dr. J. H. Nielsen Technical University of Denmark

Prof. Dr. J. Schneider Technical University of Darmstadt

Dr. M. A. Kraus ETH Zürich

<EMAIL_ADDRESS>
The present paper is concerned with deriving simplified design equations and charts for modelling in-plane expansion of fractured thermally pre-stressed glass panes using the method of equivalent temperature differences (ETD) together with a thermal expansion analogy for strains. The starting point is a theoretical method based on linear elastic fracture mechanics merged with approaches from stochastic geometry to predict the 2D-macro-scale fragmentation of glass. The approach is based on two influencing parameters of glass: (i) fragment particle size, $\delta$, and (ii) fracture particle intensity, $\lambda$, which are related to the pre-stress induced strain energy density, $U_\text{D}$, before fracture. Further Finite Element (FE) analysis of single cylindrical glass particles allow for establishing functional relations of the glass fragment particle dimensions, the pre-stress level and the resulting maximum in-plane deformation. Combining the two parts of two-parameter fracture pattern modelling and FE results on fragment expansion, formulas and engineering design charts for quantifying the in-plane expansion of thermally pre-stressed glass panes due to fracturing via an ETD approach is derived and provided within this paper. Two examples from engineering practice serve as demonstrators on how to use our ETD approach to compute the equivalent temperature difference and resulting internal forces as well as deformations. This approach serves furthermore as a basis to estimate secondary effects (such as fracture-expansion-induced deformations or stresses) on support structures or remaining parts of glass laminates in form of handy ETD load cases within analytical as well as FE analysis.

§ INTRODUCTION AND MOTIVATION

The use of glass as a structural material in engineering requires the analysis and prediction of its behaviour in the intact state, during the fracture process and, after the fracture process is completed, in the post-fracture limit state. Thermally tempered glass is often used as a structural material due to its superior strength properties compared to annealed glass. Due to the large amount of residual stresses present, tempered glass typically fractures into small dices that are less harmful than shards from broken annealed glass. It is therefore also called safety glass if the number of fragments per unit area is large enough. Contrary, if used for laminated glass, the post-fracture behaviour of tempered glass is often considered unfavourable due to the small fragments that lead to a global membrane-like structural behaviour (where almost no bending stiffness is left) of the laminated glass plates in case all glass plies are fractured. Another interesting effect is the rapid in-plane expansion of the tempered glass fragments during the fracturing process. The fracture front travels at about 1500m/s <cit.> and a significant part of the strain energy is released during this fracturing process. Within the process, each fragment expands, eventually resulting in a significant global expansion of the fractured glass plate. This expansion can lead to (i) flying debris in monolithic (single layer) glass plates, (ii) to global bending in glass laminates, where only one or several plies of a multi-laminate set-up fracture, and (iii) to delamination of the interlayer or in-plane failure of the broken glass if the expansion is strongly impeded by the remaining plies. To account for this effect in engineering glass design and to provide a handy, yet realistic, estimation of this effect, a simplified model for computing the increase in size of a glass pane made of thermally pre-stressed glass after fracture was developed within this paper. A clear analogy to thermal expansion of continuum materials given an equivalent temperature difference (ETD) is hereby followed.

§ THEORETICAL BACKGROUND AND STATE OF THE ART

This section lays the foundation for the deduction of our engineering approximation model of in-plane expansion of fractured pre-stressed glass panes.

§.§ Background on Thermally Pre-Stressed Glass

The strength of standard float glass is governed by its tensile strength, which itself is significantly influenced by small flaws in the surface, reducing the typical engineering tensile strength of annealed float glass to somewhere in the interval 30 to 100MPa, see e.g. <cit.>. Inducing a residual stress state by thermal tempering, a greater resistance to external loads, together with a certain fracture pattern in case of failure can be obtained. This leads to a desired level of safety with respect to human injuries due to the smaller fragments. For these reasons, thermally pre-stressed glass is also known as tempered safety glass. The residual stress obtained during the tempering process is characterised by its approximately parabolic distribution through the thickness, where compressive stresses reside on the outer surfaces, which is balanced by internal tensile stresses in the mid-plane, cf. Figure <ref>. This parabolic stress distribution, $ \sigma(z) $, can be written in terms of the surface stress, $ \sigma_\text{s} $, and glass thickness, $h$, as:

\begin{equation} \centering
\sigma(z)=\frac{1}{2}\sigma_\text{s}(1-3\zeta^2),\qquad\zeta=\frac{2z}{h}
\label{eq:sig_z1}
\end{equation}

Stress distribution in the far field area of a tempered glass plate and a sketch showing a contour line at zero stress (dotted line)

using symbols defined in Figure <ref>. The parabolic stress distribution is in equilibrium and symmetric about the mid-plane. The magnitude of the surface stress is approximately twice the tensile stress, $ (2\sigma_\text{m}=-\sigma_\text{s}) $. The zero stress level is at a depth of approximately $ 21\% $ of the thickness, $ h $, from the surface, known as the compressive zone depth <cit.>. The surface flaws then are in a permanent state of compression by the compressive residual stress at the surface, which has to be exceeded by externally imposed stresses due to loadings before failure and fracture of the glass pane can occur <cit.>. The magnitude of residual stresses depends on processing- and material parameters and is not within the scope of this paper. For a deeper insight on this, the reader is referred to other literature such as <cit.>.

Fragmentation of glass plates with the same thickness of 12 mm for (left) low stored strain energy ($U = 78.1~J/m^2$), (right) high stored strain energy ($U = 354.7~J/m^2$), from <cit.>

§.§ Fracture of Thermally Pre-Stressed Glasses and Glass Fracture Pattern Modelling

For the case of fracture, thermally pre-stressed glass panes fragmentize completely into many pieces, if the equilibrated residual stress state within the glass plate is disturbed sufficiently and it holds an elastic strain energy large enough. A commonly used example demonstrating this are the so-called "Prince Rupert's drops", possessing bulbous heads and thin tails. These glass drops can withstand high impact or pressure applied to the head, but ”explode” immediately into small particles if the tail is broken, see e.g. <cit.>. The fragmentation is the direct consequence of the release of elastic strain energy stored inside the material due to the residual stress state. The fragment size depends on the amount of the released energy. Small fragments are caused by a high energy release such as the high residual stress state found in tempered glass originating from the quenching process. Lower residual stress states result in larger fragments due to lower stored strain energy (cf. Figure <ref>). Thus, not only the stress but also the thickness of the glass plate plays a role in determining the strain energy. The strain energy density, $U_\text{D}$, which is the strain energy, $U$, normalised with the thickness. The strain energy density then becomes a thickness independent quantity for characterising the energy state of a thermally pre-stressed glass pane, and can then be derived to yield <cit.>:

\begin{equation} \centering
U_\text{D} =\frac{U}{h}= \frac{1}{5} \frac{(1-\nu)}{E} \sigma_\text{s}^2 = \frac{4}{5} \frac{(1-\nu)}{E} \sigma_\text{m}^2
\label{eq:U_04}
\end{equation}

where $U_\text{D}$ is the amount of elastic strain energy stored in the system per unit volume and thus only depends on the residual stress and the material properties. The above equation is fully in line with other derivations of the strain energy in tempered glass as provided by e.g. <cit.>.

Hexagonal close packing (HCP) distribution of points and resulting honeycomb pattern due to the Voronoi tessellation of HCP distributed seed points, from <cit.>

The papers <cit.> discuss the properties of fragment size, $\delta$, fragmentation intensity, $\lambda$, and the tessellation pattern that result from the fragmentation process, where especially in <cit.> a statistical evaluation and Bayesian treatment of the fracture pattern model is presented. The fracture pattern model combines an energy criterion of linear elastic fracture mechanics and tessellations induced by random point patterns. Statistical analysis of the glass fracture patterns of a comprehensive experimental programme, consisting of thermally pre-stressed glass panes with different thicknesses and levels of thermal pre-stress, allowed for a sound and exhaustive investigation the fracture pattern of tempered glass in order to determine characteristics of the fragmentation pattern (e.g fragment size, $\delta$, fracture intensity, $\lambda$, etc.). The basic modelling approach consists of the idea, that the final fracture pattern is a Voronoi tessellation induced by a stochastic point process, whose parameters can be inferred by statistical evaluation of pictures of several fractured glass specimen. By calibration of a stochastic point process and consecutive tessellation of the region of interest, statistically identically distributed realisations of fracture patterns of a glass pane can be generated. The evaluations there quintessentially show, that the size, shape and number of fragments strongly and non-linearly depend on the strain energy density, $U_\text{D}$, (cf. Figure <ref>) and the fracture pattern may be approximated by a Voronoi-tessellation induced by a Matérn-Hardcore-Point-Process. The literature, <cit.> show, that the fracture pattern is on average a Hexagonal Close Packing (HCP) with the uniform distance $\delta_\text{HCP}$ between any two adjacent nuclei. This is caused by the dynamic fracture properties as derived by <cit.> and experimentally verified for tempered glass in <cit.>. Thus the mean fracture pattern of thermally pre-stressed glass panes is a regular honeycomb with hexagonal cells, cf. Figure <ref>.

Plot of the intensity $\log_{10}(\lambda)$ - strain energy density $\log_{10}(U_\text{D})$ relation, marking two scales for the fracture pattern

The spread in fragment sizes is often found to follow the power law size distributions <cit.>, having only dimensionless fit parameters and contain no characteristic length scale, i.e. they are scale invariant. However, <cit.> can show, that unstressed glass plates follow a hierarchical breakup process with power law size distribution while stressed glass plates follow a random (Poisson) process with fragments showing an exponential size distribution with a natural characteristic length as a fit parameter linked to the residual stress.  <cit.> found, that the characteristic length scale of the exponential size distribution is approximately the thickness of the plate, $h$, which is in agreement with the findings of <cit.>.

In cf. Figure <ref> the variables of the experiments from <cit.> "pre-stress" and "thickness" are graphically encoded by colour and symbols respectively. Our analysis delivers two patterns for the relationship between energy density $U_\text{D}$ and the fragment intensity $\lambda$ (both in $\log_{10}$-scale). It is especially interesting, that the patterns are associated with two scales:

log_10(U_D) =
0.381 ·log_10(λ) + 4.78     if  λ∈10^[-3;-2]

0.643 ·log_10(λ) + 5.30     if  λ∈10^[-2;-1]

Despite that novel finding of two separated functions for relating fracture pattern parameters and strain energy density, $U_D$, the expansion model derived in Section <ref> will be based on the more simple relation reported in <cit.> for the relation between the fragment size parameter $\delta = 2r_0$ and the strain energy density $U_\text{D}$ is used:

\begin{equation}\label{eq_r_0_simple}
U_\text{D} = \frac{\SI{122.1}{\joule\per\meter\squared}}{\delta} \quad \Leftrightarrow \quad r_0 = \frac{\SI{61.05}{\joule\per\meter\squared}}{U_\text{D}}.
\end{equation}

Thermally pre-stressed glasses used for building applications reside only in the scale from -2 to -1 for $\log_{10}(\lambda)$ in Figure <ref>, which correspond to a range of 3 to 12 mm for $\delta$ in Figure <ref>. Statistical evaluation of that relation  yields $R^2=0.93$ given the data in <cit.> and hence delivers a suitable and simple model for further analysis within the context of this paper.

Elastic strain energy density, $U_\text{D}$, versus Fragment size parameter, $\delta$. Adapted and enhanced using data from <cit.>.

Having established relations between the per-stress level of the glass pane and the characteristic size of the mean glass fracture particle, the next section is concerned with the computation of the expansion of thermally pre-stressed glass fragments.

§.§ Expansion of thermally Pre-Stressed Glass Fragments

The release of residual stresses in tempered glass leads to an overall in-plane expansion of a glass plate due to straining of the individual fragments. In <cit.> it is shown how a single fragment is deforming when the pre-stress is released.

The expansion of a tempered glass plate due to fragmentation is investigated by searching for the deformations of an average fragment and then integrate the contributions over a specific glass plate to obtain the total expansion as indicated in Figure <ref>.

An efficient (in terms of computational costs) axi-symmetric FE model as described in <cit.> is applied. This indicates that each fragment is, initially, considered as a cylinder with height equal to the glass thickness, $h$, and radius, $r$, representing the in-plane size of the fragment Figure <ref>.

Initially the cylinder is considered stress free, however, in the second step a parabolic stress distribution (over the height) is applied by means of a prescribed temperature field.
In the third step boundary conditions representing the neighbouring glass are removed which will represent the fragmentation of the glass. However, due to linearity (linear elastic material, small displacements and deformations) we can skip some of the intermediate steps and apply the stress state directly on the cylinder without boundary conditions (except for those needed to prevent rigid body motions). Furthermore, symmetry can be utilised in order to reduce the computational costs even further and only 1/4 of the cylindrical cross-section was meshed as indicated in Figure <ref>. Due to the high stress variation in the fragment, a dense finite element mesh is required. For these calculations second order displacement elements were used with at least[In the FE analysis, several different models were used as the geometry (fragments size) was one of the key parameters investigated.] 100 elements through the (full height) of the fragment. This model is extremely efficient and for the results presented in this paper more than 20000 computations with varying parameters were carried out. The principles of the FE-model was experimentally validated in the paper <cit.>. The output from a single simulation, as shown in Figure <ref>, provides both stresses and deformations of a fragment.

Deformations in a fragment ($\sigma_\text{s}=\si{-75~}{MPa}$, $h=\si{8~}{mm}$, $r_0=\si{5~}{mm}$, $E=\si{70~}{GPa}$ and $\nu=\si{0.23}$). (a) Deformed cylindrical fragment. (b) Radial displacements, $u_\text{r}$, and a typical mesh. (c) Magnitude of the displacements. Arrows indicates direction and magnitude. The red square indicates the initial shape. All plots are shown with displacements magnified 100 times.

The primary result of interest for this study was the in-plane expansion of the fragments due to their horizontal deformation at the top and the bottom of each fragment. For this, it is assumed, that the stress distribution is homogeneous throughout the glass plate, which is a fairly good approximation from a distance of already two times the thickness from the edges. As glass plates are typically very thin compared to their width and length, the influence of the stress distribution at the edges is neglected.

When a pre-stressed glass plate fractures, the neighbouring fragments will "push" each other, resulting in a net expansion of the plate, as sketched in Figure <ref>, which is used for calculating the total free expansion.

Expansion of a glass plate due to deformations in neighbouring fragments. The red lines indicates the undeformed fragment geometry.

From the FEM-model, as shown in Figure <ref>, a relative maximum in-plane expansion of a fragment can be found by simply dividing the maximum in-plane displacement of the fragment, $u_\text{r}$, by the size of the (undeformed) fragment, $r_0$:

\begin{equation}
\varepsilon_\text{r} = \frac{u_\text{r}}{r_0}.
\end{equation}

This quantity is comparable to a strain and is denoted, $\varepsilon_r$, and referred to as the maximum radial strain.

By assuming that all fragments in a glass plate can be represented (on average) by the fragment used in the model, it is then possible to estimate the total free in-plane expansion of a thermally tempered glass plate. However, since we assume all fragments to expand equally, the total free expansion in any in-plane direction, $\Delta u_i$, can be calculated by simply multiplying the given plate dimension with the maximum radial strain for a representative fragment:

\begin{equation}\label{eq:ui_epsr}
\Delta u_i = \ell_i\cdot \varepsilon_\text{r}
\end{equation}

where subscript $i$ indicates a direction and $\ell_i$ is the in-plane dimension in the $i$'th direction.

From the parametric study $u_\text{r}$ is recalculated to $\varepsilon_\text{r}$ using Eq. (<ref>). A plot showing this strain as a function of the surface residual stress, $\sigma_\text{s}$, for different thicknesses can be seen in Figure <ref>.

Radial strain of a fragment as function of the residual surface stress for different thicknesses of the glass and a fixed fragment size ($r_0=\SI{5}{\mm}$)

From the Figure it is seen that the response is linear and it is therefore reasonable to normalise the strain with $\sigma_\text{s}$. Doing this allows us to plot the variations with the fragment size, $r_0$. This is shown in Figure <ref>.

Nomalised radial strain versus fragment size for different thicknesses.

From the Figure <ref> it can be noticed that all curves seems to have the same overall shape and normalising the horizontal axis with the thickness, $h$, yields all curves to coincide. This is shown in Figure <ref>.

Normalised radial strain versus normalised fragment size.

§ EQUIVALENT TEMPERATURE DIFFERENCE (ETD) MODEL FOR IN-PLANE EXPANSION OF FRACTURED THERMALLY PRE-STRESSED GLASS PANES

Modelling concrete shrinkage effects and primary as well as secondary effects on the composite structure via equivalent temperature differences (ETD) and induced linear expansion and / or curvature is well known and established for steel-concrete composite structures in both academia and engineering practice <cit.>. This paper takes the foundations laid so far to elaborate a equivalent temperature differences (ETD) model upon a thermal expansion analogy for strains to provide a handy method of estimating the average amount of in-plane expansion of thermally pre-stressed glass panes.

§.§ Deriving the ETD model for the free expansion of tempered glass at failure

In analogy to the definition of thermal strains via a thermal expansion coefficient, $\alpha_\text{th}$, and a governing temperature difference, $\Delta T$,

\begin{equation}
\label{eq:ThermalStrain}
\varepsilon_\text{th}  =\alpha_\text{th} \Delta T
\end{equation}

a fracture expansion coefficient, $\alpha_\text{fr}$ can be defined in order to establish a relation between the free expansion strain caused by fragmentation of the tempered glass, $\varepsilon_\text{fr}$, and the residual stress state, quantified by $\sigma_\text{s}<0$;

\begin{equation}\label{eq:fracStrain}
\varepsilon_\text{fr} = \alpha_\text{fr}\frac{\nu-1}{E}\sigma_\text{s}.
\end{equation}

The term, $(\nu-1)/E$, is governing the in-plane behaviour of a plate and, $E$, and, $\nu$, represents Young's modulus and Poisson's ratio respectively. The free expansion strain caused by the fracture, $\varepsilon_\text{fr}$, can be interpreted as the total strain of a given tempered glass plate upon fragmentation.

Enforcing compatibility of the ”strain from fracture”, $\varepsilon_\text{fr}$, with the maximum radial strain, $\varepsilon_\text{r}$, delivers:

\begin{equation}
\varepsilon_\text{r} = \varepsilon_\text{fr}.
\end{equation}

The fracture expansion coefficient, $\alpha_\text{fr}$, from Eq. (<ref>) can be determined through the FE-analysis on the radial strain, $\varepsilon_\text{r}$, as carried out above and summarised in the plot shown in Figure <ref>. Multiplying the curve in Figure <ref> with $E/(\nu-1)$ as indicated in Eq. (<ref>) the fracture expansion coefficient, $\alpha_\text{fr}$, can be found as a function of the fragment size relative to the glass thickness, $r_0/h$, as shown in Figure <ref>. The upper limit for the fracture expansion coefficient is $\alpha_\text{fr}=1$, which corresponds to a zero fragment size. The physical interpretation of this is that if the glass is completely pulverised, all residual stresses are released and are all converted linear elastically into the fracture expansion strain. For (unrealistically) large fragment sizes the curve tends towards zero indicating that the relative amount of residual stresses converted into fracture expansion strain approaches zero.

In Figure <ref>, the first and, from a practical point of view, most relevant part of the curve is fitted using a hyperbolic secant function, $\sech{x}=\frac{1}{\cosh{x}}$, as this function have the right properties with horizontal asymptotes for $x=0$ and $x\rightarrow \infty$. A function on the form:

\begin{equation}\label{eq_alp0}
\alpha_\text{fr}\left(\frac{r_0}{h}\right) = a_1\sech{\left(b_1\frac{r_0}{h}\right)}+\left(1-a_1\right)
\end{equation}

was therefore fitted to the plot for $\frac{r_0}{h}\leq 3.7$ with relatively good agreement as shown in Figure <ref>. This value corresponds roughly to $\sigma_\text{s}=\SI{-50}{MPa}$ for a $\SI{3}{\mm}$ glass which we will consider maximum fragment size for a standard thickness thermally pre-stressed glass.

Relation between $\alpha_\text{fr}$ and the ratio between the fragment size (radius) and the glass thickness, $r/h$. Notice the zoom providing results for more typical values of $r/h$.

Within this paper in the sense of a statistical first order expectation approximation it is assumed, that the radial strain shown in Figure <ref> is fully contributing to the plate's expansion after fracture. Hence, this paper specifies the nomenclature "$\left\langle \cdot \right\rangle$" to formally emphasise, that the derived quantity is to be interpreted as a statistical first-order approximation of the expected value of the respective quantity, so that Eq. <ref> and  <ref> yield:

\begin{equation}
\label{eq:ThermalAnalogyStrainOfFracture}
\left\langle \varepsilon_\text{th} \right\rangle = \left\langle \varepsilon_\text{fr} \right\rangle \hspace{.2cm} \iff \hspace{.2cm} \alpha_\text{th} \left\langle \Delta T \right\rangle = \alpha_\text{fr}\frac{\nu-1}{E}\sigma_\text{s}
\end{equation}

which after rearranging yields:

\begin{equation}
\label{eq:ThermalAnalogyStrainOfFracture2}
\left\langle \Delta T \right\rangle = \frac{\alpha_\text{fr}}{ \alpha_\text{th} } \frac{\nu-1}{E}\sigma_\text{s}
\end{equation}

Using the proposed ETD method allows the use of both analytical and commercial software to compute an estimate of the effects of the stress state in adjacent structural elements, see e.g. the examples in Section <ref> .

§ SUMMARY AND APPLICATION EXAMPLES

In the previous section a model relating the free expansion of a tempered glass plate with the residual surface stresses, $\sigma_\text{s}$, was derived. In this section we will summarise and provide some examples of usage, repeating some of the key equations in the model for the convenience of the reader.

The fracture expansion strain can be calculated from Eq. (<ref>):

\begin{equation}\label{eq_exp}
\varepsilon_\text{fr} = \alpha_\text{fr}\frac{\nu-1}{E}\sigma_\text{s} \quad \Leftrightarrow \quad \alpha_\text{fr}=\frac{\varepsilon_\text{fr}}{\sigma_\text{s}}\frac{E}{\nu-1}
\end{equation}

In this equation Young's modulus, $E$, and Poisson's ratio, $\nu$, for glass can be found in Table <ref>. The fracture expansion coefficient, $\alpha_\text{fr}$, can be estimated from the FEM study reported in Figure <ref> in which the most relevant part for tempered glass is fitted by an expression in form of:

\begin{equation}\label{eq_alp}
\alpha_\text{fr}\left(\frac{r_0}{h}\right) = a_1\sech{\left(b_1\frac{r_0}{h}\right)}+\left(1-a_1\right)
\end{equation}

where the constants $a_1$ and $b_1$ can be found in Table <ref>. The thickness of the glass, $h$, is assumed known.

The mean fragment size, $r_0$, can be estimated from <cit.> (cf. Eq.  (<ref>)) and the relation $\delta=2r_0$ as:

\begin{equation}\label{eq_r0}
r_0 = \frac{a_2}{U_\text{D}}
\end{equation}

where, $a_2$, can be found in Table <ref> and, $U_\text{D}$, is the strain energy density for tempered glass given by Eq. (<ref>):

\begin{equation}\label{eq_UD}
U_\text{D} = \frac{1-\nu}{5E}\sigma_\text{s}^2
\end{equation}

Parameters to be used for the ETD model.

Parameter                   Reference

$a_1$=0.862            2*Figure <ref>

$b_1$ = 1.686
$a_2$ = 61.05        Eq. (<ref>)

Combining Eqs. (<ref>),(<ref>), (<ref>) and Eq. (<ref>) we obtain:

\begin{equation}
\begin{split}
\varepsilon_\text{fr} =&
\frac{\nu-1}{E}\left(a_1\sech{\left(\frac{5 E a_{2} b_{1}}{h \sigma_\text{s}^{2} \left(1 - \nu\right)} \right)}+1-a_1\right)\sigma_\text{s}
\end{split}
\end{equation}

Inserting all constants from Table <ref> and rearranging, the expression can be written as:

\begin{equation}\label{eq:eps_simple}
\begin{split}
\varepsilon_\text{fr} &=
\frac{-\sigma_\text{s}}{\SI{658.8e3}{\MPa}}\left[1\phantom{6.246\sech{\left(\frac{\SI{4.679e16}{\meter\MPa\squared}}{h\sigma_s^2}\right)}}\right.\\
\end{split}
\end{equation}

From Eq. (<ref>) a plot showing the fracture strain, $\varepsilon_\text{fr}$, as a function of the residual surface stress, $\sigma_\text{s}$, for different glass thicknesses, $h$, is generated and shown in Figure <ref>:

The fracture strain as a function of the residual surface compression for different glass thicknesses.

§.§ Example 1: ETD Model Prediction for the expansion due to fracture of a mono tempered glass pane

In this example we will consider a $\ell_x \times \ell_y = \SI{2.5x1}{\m}$, 10 thick monolithic tempered glass pane. The residual surface stress is measured to $\sigma_s = \SI{-85}{MPa}$. The total free expansion in case of failure can now be estimated using Eq. (<ref>) and Eq. (<ref>) or Figure <ref>.

According to Eq. (<ref>), the fracture strain in the glass pane is $\varepsilon_\text{fr}=\num{791e-6}$. The total expansion in the $x$ and $y$ directions, then becomes:

\begin{equation}
\begin{split}
u_x &= \varepsilon_\text{fr}\cdot\ell_x =\num{791e-6}\cdot\SI{2500}{\mm}= \SI{1.98}{\mm}\\
u_y &= \varepsilon_\text{fr}\cdot\ell_x =\num{791e-6}\cdot\SI{1000}{\mm}=\SI{0.791}{\mm} \\
\end{split}
\end{equation}

The equivalent temperature for use inside a Finite Element software can be computed using Eq. <ref>, which for this example yields:

\begin{equation}
\label{eq:ThermalAnalogyStrainOfFracture_example1}
\left\langle \Delta T \right\rangle = 88 K
\end{equation}

assuming $\alpha_\text{th}=\SI{9.1e-6}{\per\kelvin}$ which is commonly used for glass.

§.§ Example 2: ETD model prediction for partly fractured laminated glass

The model can also be used for estimating the extra load on intact panes in (partly) fractured laminated glass.

Considering a two ply laminated glass as shown in Figure <ref>a. In the initial configuration, no external stresses are present and the two plies have the thickness $h_1$ and $h_2$.

Laminated glass, ply 1 is tempered glass. a) shows the initial configuration, b) shows the free expansion of the tempered glass and c) shows the total elongation of the laminated glass.

In the intermediate pseudo configuration, see Figure <ref>b, the fractured tempered glass (ply 1) is shown with its free expansion and the intact glass (ply 2) is not affected by this. Obviously, the expansion will be transferred through the interlayer and a first approximation to the final state can be found by assuming a stiff interlayer and no rotations (bending), see Figure <ref>c. Now we can find the forces in the final configuration by first equal length for the two plies:

\begin{equation}
\begin{split}
L(1+\varepsilon_\text{fr})+L(1+\varepsilon_\text{fr})\cdot\varepsilon_1 &= L+L\cdot\varepsilon_2 \qquad \Leftrightarrow\\
\varepsilon_\text{fr}+(1+\varepsilon_\text{fr})\cdot\varepsilon_1 &= \varepsilon_2\\
\end{split}
\end{equation}

Now applying Hookes law and requiring equilibrium, $F_1=-F_2=F$ we find:

\begin{equation}
\varepsilon_\text{fr}+(1+\varepsilon_\text{fr})\cdot\frac{-F}{E_1h_1} - \frac{F}{E_2h_2}=0
\end{equation}

from which the force, $F$, and the stress in the intact ply can be found as:

\begin{equation}\label{eq:example02}
\end{equation}

As an example, one could consider the 10 glass with a surface residual stress, $\sigma_{\text{s}}=\SI{-85}{MPa}$ from Section <ref> for both layers in a two-ply laminated glass plate. If the one layer fails, the stresses in the intact layer can be estimated from Eq. (<ref>). Assuming same stiffness for both plies, $h_1=h_2=\SI{10}{\mm}$ and $E_1=E_2=\SI{70}{\GPa}$ we find, in this case, the stress in the intact pane to be $\sigma=\SI{28}{\MPa}$.

An unknown in the model is the compressive stiffness for fractured tempered glass. However, to the best knowledge of the authors, such investigations have not yet been published. It is expected that the stiffness of the broken glass is lower compared to the intact glass and is likely to be strain dependent.

The plot in Figure <ref> shows an example of the tensile stress in the intact layer, $\sigma_2$, as a function of the compressive stiffness of the broken layer, $E_1$, for varying thicknesses of the intact layer, $h_2$. From the results in the figure, it is seen that combining a 15 fully tempered glass with e.g. a 6 thick annealed glass may cause problems in case of failure of the tempered glass as the peak stress in the annealed glass may reach up to 65.

Peak stress in intact ply, $\sigma_2$, as a function of the stiffness in the broken ply, $E_1$, for different thicknesses of the intact ply, $h_2$ (in mm). The broken ply is having a thickness of $h=\SI{15}{\mm}$ and a residual surface stress of $\sigma_\text{s}=\SI{-120}{\MPa}$.

The proposed model does not take into account the stiffness of the inter-layer and bending of the plies due to asymmetric failure. However, the fracture velocity of tempered glass has been measured using high-speed cameras to be approximately 1466 <cit.>, indicating a very high loading rate of the inter-layer and thereby a dynamic problem. It is well known, that common interlayers, such as Poly Vinyl Butyral (PVB), Ethylene Vinyl Acetate (EVA) and ionomers e.g. SentryGlas®, shows an increase in stiffness with the loading rate. This supports the assumption of a full shear transfer between the plies, however, it also indicates that some of the load might actually be carried by the interlayer, which is not accounted for in the current model. Due to the dynamic nature of the problem, the strength of the intact glass may also be higher than what is often assumed for quasi static problems, see e.g. <cit.> where a review on available strength data for soda-lime-silica glass is given.

§ SUMMARY AND CONCLUSION

This paper first presented experimental and theoretical background on the glass fracture process for thermally pre-stressed glasses. The 2D-macro-scale fragmentation of glass can be basically described by the fragment particle size, $\delta$, and the fracture particle intensity, $\lambda$, which both are related to the pre-stress induced strain energy density, $U_\text{D}$, before fracture. Then further details on Finite Element (FE) simulations of single cylindrical glass particles are reported, which allowed to establish functional relations of the glass fragment particle dimensions, the pre-stress level and the resulting maximum in-plane deformation. These results are then combined with the two-parameter fracture pattern modelling to furnish an equivalent temperature differences (ETD) for describing the in-plane expansion of thermally pre-stressed glass panes due to fracturing. Finally, two examples from engineering practice demonstrated the application of the developed graphs and formulas for further use in analytical as well as FE analysis of fractured glass laminates.

Further analysis of the fracture particle statistics vs. strain energy density proved existence of two glass fracture statistic regions. Despite that novel finding within Section <ref>, for engineering practice it is sufficient to only concentrate on one of the two fracture domains (3mm to 12mm fragment particle size). The idea of ETD then was applicable straight forward for relating FE analysis results of the deformations of a single fragment upon failure of the glass with the total free expansion of a piece of tempered glass from its initial dimensions and residual surface compressive stress by using a thermal strain analogy. Here, the derived ETD model provides a tool for estimating the free (unconstrained) expansion of tempered glass at failure, which is not possible at all at the moment. To that end, our approach in form of handy ETD load cases within both, analytical as well as FE analysis, allows for the estimation of (i) secondary effects in the fractured laminate such as fracture-expansion-induced deformations or stresses, and (ii) effects due to second order influences on residual load bearing capacity of the fractured glass laminate as well as of support structures or remaining parts of glass laminates. Furthermore, a simple analytical model for estimating the peak stress in the intact ply in a partly broken laminated glass plate is provided. From the model it is found that the peak stress must be considered relevant and the model also suggest that care should be taken if laminating glass plies with too different thicknesses or mixing both tempered and annealed glass. Future research needs to address experimental validation of this ETD model on various glass laminates, where level of pre-stress as well as laminate size and glass thicknesses are varied. Furthermore, this ETD model needs to be enhances for the influence of different interlayer types as these possess pronounced differences in stiffness and hence shear coupling of the glass panes during fracture and speed of redistribution of internal forces in the post-fractured state.


# VLEIBot: A New $\boldsymbol{45}$-mg Swimming Microrobot Driven by a
Bioinspired Anguilliform Propulsor
††thanks: This work was partially funded by the Washington State University
(WSU) Foundation and the Palouse Club through a Cougar Cage Award to N. O.
Pérez-Arancibia. Additional funding was provided by the WSU Voiland College of
Engineering and Architecture through a start-up fund to N. O. Pérez-Arancibia.
††thanks: E. K. Blankenship and C. K. Trygstad contributed equally to this
work.††thanks: The authors are with the School of Mechanical and Materials
Engineering, Washington State University (WSU), Pullman, WA 99164, USA.
Corresponding authors’ e-mail<EMAIL_ADDRESS>(C. K. T.);
<EMAIL_ADDRESS>(N. O. P.-A.).

Elijah K. Blankenship, Conor K. Trygstad, Francisco M. F. R. Gonçalves, and
Néstor O. Pérez-Arancibia

###### Abstract

This paper presents the VLEIBot111We pronounce this word as vl-ahy-bot, thus
producing the same sound in the pronunciations of eigenvalue and Einstein.
(Very Little Eel-Inspired roBot), a
$\boldsymbol{45}$-mg/$\boldsymbol{23}$-$\text{mm}^{\boldsymbol{3}}$
microrobotic swimmer that is propelled by a bioinspired anguilliform
propulsor. The propulsor is excited by a single $\boldsymbol{6}$-mg high-work-
density (HWD) microactuator and undulates periodically due to wave propagation
phenomena generated by fluid-structure interaction (FSI) during swimming. The
microactuator is composed of a carbon-fiber beam, which functions as a leaf
spring, and shape-memory alloy (SMA) wires, which deform cyclically when
excited periodically using Joule heating. The VLEIBot can swim at speeds as
high as $\boldsymbol{15.1}\,\text{mm}\cdot\text{s}^{\boldsymbol{-1}}$
($\boldsymbol{0.33}\,\text{Bl}\cdot\text{s}^{\boldsymbol{-1}}$) when driven
with a heuristically-optimized propulsor. To improve maneuverability, we
evolved the VLEIBot design into the
$\boldsymbol{90}$-mg/$\boldsymbol{47}$-$\text{mm}^{\boldsymbol{3}}$ VLEIBot+,
which is driven by two propulsors and fully controllable in the two-
dimensional ($\boldsymbol{2}$D) space. The VLEIBot+ can swim at speeds as high
as $\boldsymbol{16.1}\,\text{mm}\cdot\text{s}^{\boldsymbol{-1}}$
($\boldsymbol{0.35}\,\text{Bl}\cdot\text{s}^{\boldsymbol{-1}}$), when driven
with heuristically-optimized propulsors, and achieves turning rates as high as
$\boldsymbol{0.28}\,\text{rad}\cdot\text{s}^{\boldsymbol{-1}}$, when tracking
path references. The measured root-mean-square (RMS) values of the tracking
errors are as low as $\boldsymbol{4}\,\text{mm}$.

##  I Introduction

The recent development of extremely light ($1$ to $10$ mg) and small ($0.45$
to $1.89$ mm3) fast (up to $40\,\text{Hz}$) actuators based on shape-memory
alloy (SMA) technology enabled the creation of new types of microrobots, which
can operate at low voltages and with simple excitation electronics [1, 2, 3].
These advancements predict that the vision of creating swarms and schools of
insect-scale robots capable of performing essential tasks for humans will
become a reality soon. Here, we present two new microrobotic surface swimmers,
whose propulsion mechanism is inspired by the anguilliform swimming mode at a
moderate ($\sim 10^{2}$ to $\sim 10^{4}$) Reynolds number ($Re$). These
microrobots, shown in Fig. 1, are: (i) the VLEIBot, which weighs $45$ mg,
operates at frequencies up to $20$ Hz, and achieves speeds as high as
$15.1\,\text{mm}\cdot\text{s}^{-1}$ ($0.33\,\text{Bl}\cdot\text{s}^{-1}$);
and, (ii) the VLEIBot+, which weighs $90\,$ mg, operates at frequencies up to
$10$ Hz, and achieves speeds as high as $16.1\,\text{mm}\cdot\text{s}^{-1}$
($0.35\,\text{Bl}\cdot\text{s}^{-1}$). Additionally, because of its dual
propulsion mechanism, the VLEIBot+ is fully controllable in the two-
dimensional ($2$D) space. We envision that in coordination with other flying
and crawling insect-scale robots—advanced evolutions of those in [1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]—teams of VLEIBot-type
microswimmers will assist humans in aquatic search and rescue missions,
inspection of vessels and marine infrastructure, hydroponic agriculture,
aquaculture tasks, research of shallow-water reefs, and minimally-invasive
continuous automatic monitoring of water quality in reservoirs, just to
mention a few.

Figure 1: Photograph of the VLEIBot and VLEIBot+. The VLEIBot (left) is a
$45$-mg swimming microrobot driven by a new bioinspired propulsor. The
VLEIBot+ (right) is a $90$-mg controllable microrobotic swimmer designed to
achieve high maneuverability in the two-dimensional ($2$D) space.

The VLEIBot design exploits the surface tension of water to stay afloat and
the fluid-structure interaction (FSI) of a flexible thin undulating tail with
water to propel itself forward. In this case, the propulsor is driven by the
$6$-mg high-work-density (HWD) SMA-based actuator introduced in [1] and its
design was loosely inspired by the anguilliform swimming mode exhibited by
European eels (Anguilla anguilla L.) [19], yellow-bellied sea snakes (Pelamis
platurus) [20], and young African clawed frog (ACF) tadpoles (Xenopus laevis)
[21]. Anguilliform swimming is believed to be generated by sub-ambient
pressure zones in the troughs of the swimmer’s body waves and, for this
reason, highly efficient regarding cost of transport. This is an FSI
phenomenon observed at moderate $Re$s ($\sim 10^{2}$ to $\sim 10^{4}$) [22,
23]. Given that the $Re$s associated with the motion of the VLEIBot’s
propulsor are in the moderate range ($\sim 10^{-1}$ to $\sim 10^{3}$), we
hypothesize that the same hydrodynamic effect occurs in this case and we
considered it during the design process. In the past, researchers have
distributively embedded SMA wires in semi-rigid bodies to generate undulation
[24, 25, 26, 27, 28, 29]; however, to our best knowledge, the VLEIBot’s design
is the first with onboard actuation that removes the need for complex
propulsor configurations and fully exploits FSI phenomena to generate a wave
that propagates along the swimmer’s tail, thus producing locomotion. This is
the key design element that enables anguilliform-inspired swimming at the mm-
to-cm scale for the first time.

The most common actuation methods used to drive swimming microrobots are based
on electromagnetic [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
44, 45, 46, 47, 48], piezoelectric [49, 50, 51, 52, 53, 54, 55], and ionic-
polymer [56, 57, 58] technologies. While electromagnetic actuation has enabled
the functionality of very small and fast swimmers, some of which can operate
inside the human body, this approach does not represent a path to full
autonomy at the mm-to-cm scale because it requires the generation of large
external magnetic fields in laboratory environments. The main advantage of
piezoelectric actuation is its wide frequency bandwidth; however, it requires
very high excitation voltages ($\sim 300\,\text{V}$) and, therefore, complex
excitation electronics. Wide actuation bandwidths are required for creating
insect-scale flapping-wing flyers, but the observed undulation frequencies in
anguilliform swimmers do not overpass $10\,\text{Hz}$ [22]. Ionic-polymer
actuation provides unique new capabilities for soft robotic systems; however,
its relatively low force outputs significantly limit its applicability to the
development of aquatic propulsion systems. In contrast with other actuators
for microrobotics, the actuator that drives the VLEIBot is not only light
($6\,\text{mg}$) and fast (up to $30\,\text{Hz}$), but also can generate large
forces ($\sim 84\,\text{mN}$) using very low voltages ($5$ to $25\,\text{V}$)
and simple excitation electronics. To heuristically optimize the hydrodynamic
design and operation parameters of the VLEIBot, we used a Vicon motion-capture
system to investigate the relationship between the robot’s average speed and
tail geometry through a series of characterization experiments. We discovered
that, among the cases tested, the fastest speed of
$15.1\,\text{mm}\cdot\text{s}^{-1}$
($0.33\,\text{Bl}\cdot\text{s}^{\boldsymbol{-1}}$) is attained using a
parabolic tail with a constant aspect ratio of $0.41$ and length of
$26\,\text{mm}$, undulating at a frequency of $5\,\text{Hz}$.

Simple swimming experiments showed that, due to asymmetrical force production
of the SMA-based exciting actuator and high sensitivity to external
disturbances, the VLEIBot exhibits low maneuverability. To address this
performance challenge, we developed the $90$-mg VLEIBot+, a modular design
that achieves full $2$D controllability by using two $6$-mg SMA actuators that
independently drive two anguilliform-inspired propulsors. From a control
perspective, this approach follows the methodology first introduced in [2]. In
the case presented here, we achieved high turning capabilities. Using a Vicon
motion-capture system, to evaluate and demonstrate the locomotion performance
that the VLEIBot+ can achieve, we performed several swimming open-loop and
closed-loop control experiments. In addition to swimming at speeds on the
order of $16\,\text{mm}\cdot\text{s}^{-1}$
($0.35\,\text{Bl}\cdot\text{s}^{\boldsymbol{-1}}$), the VLEIBot+ can turn at
rates as high as $0.28\,\text{rad}\cdot\text{s}^{-1}$. The measured root-mean-
square (RMS) tracking-error values are as low as $3.94\,\text{mm}$ when
operating at the optimal frequency of $5\,\text{Hz}$. These results are highly
significant because it is well known that as the $Re$ decreases, swimming
becomes increasingly challenging by the traditional methods used at the human
scale ($Re\approx 10^{4}$) [59]. The VLEIBot and VLEIBot+ highlight the
possibility of creating mm-to-cm–scale autonomous swimming microrobots that
can work together in schools (teams) to perform complex tasks useful to
humans.

The rest of the paper is organized as follows. Section II presents the design
and fabrication of the $45$-mg VLEIBot. Section III describes the experimental
setup used in the swimming experiments, and discusses the characterization and
performance evaluation of two sets of bioinspired soft propulsors. Section IV
discusses the development of the VLEIBot+ and presents closed-loop swimming
experiments. Finally, Section V states some conclusions regarding the
presented research.

Figure 2: Robotic design and functionality. (a) Components and assembly of the
VLEIBot. The VLEIBot is composed of four main components: (i) a rigid body
that supports the weight of the robot on water due to surface tension; (ii) a
$6$-mg HWD SMA-based actuator with passive hinges installed at its two distal
ends; (iii) a planar four-bar transmission mechanism that maps the output
displacement generated by the SMA-based actuator into the angular oscillation
that excites the undulating tail of the swimmer; and, (iv) an anguilliform-
inspired soft tail made of fluoropolymer film. The triplet
$\left\\{\boldsymbol{n}_{1},\boldsymbol{n}_{2},\boldsymbol{n}_{3}\right\\}$
denotes the inertial frame used to describe the kinematics of the system. (b)
Transmission functionality. The planar four-bar mechanism can be adjusted to
have a constant bias angle, $\phi_{\text{b}}$, by shifting the installation
point of the fixed end by a distance $d_{\text{b}}$. (c) Forward locomotion.
The VLEIBot is designed to function with a symmetric stroke envelope (sweeping
area). The SMA-based actuator produces the periodic displacement output that
is mapped by the transmission into the large stroke angles required to
undulate the robot’s tail. The triplet
$\left\\{\boldsymbol{b}_{1},\boldsymbol{b}_{2},\boldsymbol{b}_{3}\right\\}$
denotes the body-fixed frame used to describe the kinematics of the system.
(d) Right turn. Theoretically, a right-biased undulation produces a negative
torque with respect to the body frame and, as a consequence, a right turn. (e)
Left turn. Theoretically, a left-biased undulation produces a positive torque
with respect to the body frame and, as a consequence, a left turn.

##  II Design and Fabrication

As shown in Fig. 2(a), the VLEIBot has four main components: (i) a rigid body
that supports the weight of the entire robot on water due to surface tension;
(ii) a $6$-mg SMA-based actuator, which is installed in the robot with passive
hinges attached to its two distal ends; (iii) a planar four-bar transmission;
and, (iv) a soft tail heuristically designed to generate anguilliform-inspired
locomotion. All composing parts were designed and fabricated using the methods
described in [6, 9, 8, 7] and references therein. The body of the robot is
made of carbon fiber (CF); the actuator is made of CF, FR4, and SMA wires,
according to the configuration in [1, 2]; the transmission and hinges are made
of CF and Kapton; and, the soft tail is made of $25$-µm-thick fluoropolymer
film (AirTech A$4000$R$14417$). The basic configuration of the body is simply
the perpendicular intersection of two flat CF pieces with significant areas
relative to their volumes. One piece interacts horizontally with the water
and, because of a relatively high surface tension, maintains the robot afloat;
the other piece functions as a keel, thus providing stability and high lateral
drag compared to that acting on the tail. The $6$-mg HWD SMA-based actuator
excites the four-bar transmission; the transmission excites the flexible tail,
which periodically undulates during swimming to produce thrust. To select the
geometrical parameters that maximize the transmission ratio, we used the
nonlinear relationship discussed in [60]. Also, we incorporate
interchangeability and enforce alignment in the tail installation process by
using a keyway system integrated into the transmission.

As discussed in [1], actuators of the type employed to drive the VLEIBot
accumulate a steady-state bias during operation at pulse-width modulation
(PWM) frequencies higher than $1\,\text{Hz}$ because of the limited time
available for the composing SMA wires to cool down within an actuation cycle.
To address this issue, the transmission is tuned and installed with a fixed
bias angle, $\phi_{\text{b}}$, as depicted in Fig. 2(b). This bias angle is
created by simply displacing the fixed end of the transmission, as graphically
defined in Fig. 2(b), by a distance $d_{\text{b}}$; namely,
$\phi_{\text{b}}=T\cdot d_{\text{b}}$, where $T$ is the transmission ratio.
For small displacements, the transmission ratio can be estimated as
$T=\delta^{-1}$ [60], where $\delta$ is the offset between the fixed and free
ends of the mechanism, as graphically defined in Fig. 2(b). In the case of
PWM-based Joule heating, because of the $\phi_{\text{b}}$ adjustment, when the
actuator contracts during an actuation cycle, the direction of bending is
opposite to that of $\phi_{\text{b}}$. This effect allows us to generate a
symmetric stroke envelope (sweeping area) with respect to the body-fixed
$\boldsymbol{b}_{1}$-$\boldsymbol{b}_{3}$ plane, for any given operational
frequency, as depicted in Fig. 2(c). Unfortunately, different operational
frequencies require different $\phi_{\text{b}}$ adjustments. Nonetheless, test
results indicate that the VLEIBot can achieve straight forward locomotion in
open loop after its transmission is tuned for a predefined operational
condition. The tuning procedure is very simple. First, we install the
transmission with a bias angle of approximately $0.3$ rad; then, we excite the
VLEIBot with a signal that produces significant forward locomotion, selected
through the tests presented and discussed in Section III; last, we visually
determine the magnitude and direction of $\phi_{\text{b}}$ required to adjust
the transmission. We repeat this process, if necessary, until we observe a
symmetric stroke envelope.

By generating a symmetric stroke envelope, we ensure that the total
hydrodynamic force produced by the undulating tail has a lateral average of
zero in steady state and, therefore, its resulting average direction is
ideally aligned with the body-fixed $\boldsymbol{b}_{1}$ axis. However, due to
external disturbances, such as the pull of the tether wire and wind gusts, the
VLEIBot often deviates from a given straight path in open loop operation. This
observation highlights the need for feedback control, even after the stroke
envelope has been tuned. Theoretically, basic control strategies can be
implemented by simply modulating the PWM actuator-exciting signal and thus
create an asymmetrical tail undulation with respect to the body-fixed
$\boldsymbol{b}_{1}$ axis (see Fig. 2(c)). As illustrated in Figs. 2(d) and
(e), theoretically, an asymmetrical stroke envelope produces a torque that
rotates the robot during swimming, thus endowing it with turning capabilities
and controllability in $2$D. Unfortunately, empirical evidence indicates that
this turning ability of the robot is highly sensitive to slight miscalibration
of the stroke envelope and external disturbances. Fortunately, as discussed in
Section IV, this steerability issue can be solved with a modular design that
puts together two VLEIBot platforms to create the VLEIBot+.

Figure 3: Experimental setup used in the tail-characterization swimming tests.
(a) Experimental setup. A Simulink Real-Time host-target system and an AD/DA
board (National Instruments PCI-$6229$) are used to generate the PWM signal
that excites the robot. This signal is then filtered through a waveform
amplifier (Accel Instruments TS$250$-$02$) to provide the power necessary to
periodically Joule heat the SMA-based actuator that drives the VLEIBot. During
the swimming tests, the VLEIBot is placed in a pool filled with water and a
four-VK16-camera Vicon motion-capture system is used to measure the
instantaneous position and orientation of the robot at a rate of
$250\,\text{Hz}$. (b) Parameters of the parabolic tails with constant aspect
ratio. For the tail-characterization experiments, we kept the aspect ratio
constant at $\AR_{\text{p}}=0.41$ and varied the tail length, $l_{\text{p}}$,
in increments of $2\,\text{mm}$ over the range $\left[0:28\right]\,\text{mm}$.
The height $h_{\text{p}}$ directly depends on the length $l_{\text{p}}$. (c)
Parameters of the rectangular tails with constant height. For the tail-
characterization experiments, we kept the height, $h_{\text{r}}$, constant at
$4\,\text{mm}$ and varied the length, $l_{\text{r}}$, in increments of
$5\,\text{mm}$ over the range $[0:50]\,\text{mm}$.

##  III Tail Characterization

###  III-A Experimental Setup and Swimming Tests

An illustration of the experimental setup used in the swimming tests is shown
in Fig. 3(a). As seen here, to generate the PWM signal that excites the robot,
and for signal processing and control, we use a Simulink Real-Time host-target
system and a National Instruments PCI-$6229$ analog-digital/digital-analog
(AD/DA) board. After the PWM signal is generated, this is filtered through a
waveform amplifier (Accel Instruments TS$250$-$02$) to supply the power
necessary to cyclically Joule heat the SMA material of the robot’s actuator.
During the experiments, a VLEIBot prototype is placed in a pool filled with
water while a four-VK16-camera Vicon motion-capture system measures and
records its instantaneous position, $\boldsymbol{r}$, and attitude at a
frequency of $250\,\text{Hz}$. We estimate the velocity, $\boldsymbol{v}$, and
angular velocity, $\boldsymbol{\omega}$, of the robot using a discrete-time
derivative algorithm and, then, these data are further processed using a low-
pass finite impulse response (FIR) filter with order $10^{3}$ and cutoff
frequency of $40\,\text{Hz}$. For consistency with the Vicon system, the host-
target system samples, processes, generates, and records signals at
$250\,\text{Hz}$.

To characterize the effects of tail geometry and size on the swimming
performance of the VLEIBot, we experimentally investigated two basic tail
geometries: (i) a low-aspect-ratio parabolic-shaped tail, loosely inspired by
the tails of the European eel and young ACF tadpole; and, (ii) a thin
rectangular-shaped tail, loosely inspired by the tail of the yellow-bellied
sea snake. During the characterization process, we tested several aspect
ratios of the two basic geometries selected, over a range of excitation
frequencies. Then, using the measured data, we heuristically determined the
best-performing tail for the VLEIBot regarding speed of locomotion. In the
tail-characterization experiments and analyses, we estimate the Reynolds
number as

$\displaystyle Re=\frac{lv\rho}{\mu},\vspace{-0.7ex}$ (1)

in which $l$ and $v$ are the length (including the tail) and speed of the
swimmer, respectively; and, $\rho$ and $\mu$ are the density and dynamic
viscosity of water at $20\,^{\circ}\textrm{C}$, respectively. Considering all
achieved speeds and tested tails, the $Re$ is inside the range
$\left[10^{-1}:10^{3}\right]$. As a reference, the $Re$ for the swimming of
the young ACF tadpole is on the order of 250 [21]. Also, similarly to the
VLEIBot, the young ACF tadpole has a wide rigid head and a thin highly-
flexible tail. Furthermore, the lack of osseous or cartilaginous support in
their tails [23], prompt us to speculate that the young ACF tadpoles might
take advantage of FSI phenomena to locomote.

The similarities between the VLEIBot and young ACF tadpoles are the main
reasons to use the anatomy of these animals as a model for tail design.
Specifically, we selected a constant aspect ratio, specified as
$\AR_{\text{p}}=h_{\text{p}}^{2}\cdot A_{\text{p}}^{-1}$, where $h_{\text{p}}$
and $A_{\text{p}}$ are respectively the height and area of the tail, both
graphically defined in Fig. 3(b). Assuming a parabolic shape, we determined
that the aspect ratio is approximately $0.41$ for the tails of young ACF
tadpoles. By keeping the aspect ratio of each tested tail fixed, its height
becomes a function of its length, $l_{\text{p}}$, which is defined in Fig.
3(b). With this in mind, for the constant aspect ratio $\AR_{\text{p}}=0.41$,
we made and tested $14$ different tails in the set $\left[2:2:28\right]$ mm.
To design the tested rectangular tails, loosely inspired by the tail geometry
of yellow-bellied sea snakes, we simply selected the constant height
$h_{\text{r}}=4$ mm, as defined in Fig. 3(c), and varied the length,
$l_{\text{r}}$, as illustrated in Fig. 3(c), over the range
$\left[5:50\right]$ mm in increments of $5\,\text{mm}$. To evaluate the
performance associated with the tested tails, we measured the average forward
speed of the robot across all the considered tail lengths and shapes. For each
tested pair of exciting frequency and tail shape, we performed $10$ swimming
tests. The tested frequencies are $1$, $5$, $10$, $15$, and $20$ Hz with a
heuristically chosen duty cycle—defined as the percentage of the PWM period
for which the signal is at its on voltage—of $8\,\%$, for the excitation of
$1\,\text{Hz}$, and of $10\,\%$, for the other tested frequencies. An on
voltage of $22\,\text{V}$ was selected such that the actuator is provided with
$250\,\text{mA}$ of current; higher current values were observed to burn the
thin $53$ AWG tether wires that power the robot. For analysis, the average
speed of each swimming experiment was computed from $20\,\text{s}$ of data.
Then, we found the mean and empirical standard deviation (ESD) of the $10$
tests corresponding to each pair of PWM frequency and tail shape. These data
are summarized in Fig. 4 and discussed next.

Figure 4: Tail-characterization data. (a) Characterization of parabolic tails
with constant aspect ratio. During the swimming tests, the aspect ratio was
kept constant at $\AR_{\text{p}}=0.41$, while we investigated the relationship
between tail length and average forward speed across frequencies of $1$, $5$,
$10$, $15$, and $20$ Hz. We see an increasing trend in average speed as tail
length increases. At the tail length of $26\,\text{mm}$ and frequency of
$5\,\text{Hz}$, we measured the maximum average speed of
$15.1\,\text{mm}\cdot\text{s}^{-1}$ ($0.33\,\text{Bl}\cdot\text{s}^{-1}$).
When the tail length, $l_{\text{p}}$, was increased further, the weight of the
tail tended to drag the robot under the surface of the water. (b)
Characterization of rectangular tails with constant height. During the
swimming tests, the tail height was kept constant at
$h_{\text{r}}=4\,\text{mm}$, while we investigated the relationship between
tail length and speed. We varied the tail length in increments of
$5\,\text{mm}$ over the range $\left[5:50\right]\,\text{mm}$. It can be
observed that high-frequency actuation is not able to produce thrust with any
of the investigated tails, and $1\,\text{Hz}$ is the best actuation frequency
for all the tested profiles. At the length of $20\,\text{mm}$, we see the best
swimming performance of $4.7\,\text{mm}\cdot\text{s}^{-1}$
($0.12\,\text{Bl}\cdot\text{s}^{-1}$).  Figure 5: Locomotion experiments
performed using optimized parabolic $\boldsymbol{26}$-mm tail with constant
aspect ratio of $\boldsymbol{0.41}$. (a) Photographic composite from video
footage showing the VLEIBot swimming at $1\,\text{Hz}$. The duty cycle of the
driving PWM signal is $8\,\%$. (b) Photographic composite from video footage
showing the VLEIBot swimming at $5\,\text{Hz}$. The duty cycle of the driving
PWM signal is $12\,\%$. (c) Photographic composite from video footage showing
the VLEIBot+ swimming at $1$ Hz. The duty cycles of the driving PWM signals
are $8\,\%$. (d) Photographic composite from video footage showing the
VLEIBot+ swimming at $5\,\text{Hz}$. The duty cycles of the driving PWM
signals are $12\,\%$. Video footage of these experiments can be viewed in the
accompanying supplementary movie.

###  III-B Results and Discussion

The experimental data corresponding to the parabolic tails with a constant
aspect ratio $\AR_{\text{p}}=0.41$ are shown in Fig. 4(a). Here, we can
clearly see a trend; over the range $\left[2:26\right]\,\text{mm}$, as the
tail length, $l_{\text{p}}$, increases, the average speed, $\bar{v}$, tends to
increase. After passing a length of $26\,\text{mm}$, $\bar{v}$ tends to
decrease because of the relatively high weight of the tail compared to that of
the robot. With a frequency of $5\,\text{Hz}$ and
$l_{\text{p}}=26\,\text{mm}$, we obtained the best swimming performance in
terms of speed, $15.1\,\text{mm}\cdot\text{s}^{-1}$
($0.33\,\text{Bl}\cdot\text{s}^{-1}$). The experimental data corresponding to
the rectangular tails with constant height $h_{\text{r}}=4\,\text{mm}$ are
shown in Fig. 4(b). In these tests, for all tail lengths, the best frequency
is $1\,\text{Hz}$; for this frequency, the best average speed is
$4.7\,\text{mm}\cdot\text{s}^{-1}$ ($0.12\,\text{Bl}\cdot\text{s}^{-1}$),
corresponding to a tail length of $20\,\text{mm}$. At frequencies higher than
$5\,\text{Hz}$, essentially, no locomotion can be produced with the tested
rectangular tails. Overall, the swimming performances achieved with the
parabolic tails are significantly better than those obtained with the
rectangular tails. Through this heuristic process, we selected the parabolic
$26$-mm tail for the VLEIBot prototype. Photographic composites of frames
taken at intervals of $6\,\text{s}$ over $24\,\text{s}$, with the optimized
VLEIBot swimming at $1$ and $5\,\text{Hz}$ are shown in Figs. 5(a) and (b),
respectively. Video footage of these two experiments are shown in the
accompanying supplementary movie.

##  IV The VLEIBot+: A Control Approach

###  IV-A Controllable Design

As explained in previous sections, after multiple open-loop swimming
experiments, we determined that the VLEIBot design lacks the steerability
required to execute high-speed turning maneuvers. To address this issue, we
developed the $90$-mg dual-propulsor VLEIBot+, which is fast, highly
maneuverable, and fully controllable in the $2$D space. Photographic
composites of frames taken at intervals of $3\,\text{s}$ over $12\,\text{s}$,
with a VLEIBot+ prototype swimming at $1$ and $5\,\text{Hz}$, are shown in
Figs. 5(c) and (d), respectively. Video footage of these two experiments are
shown in the accompanying supplementary movie. A VLEIBot+ is built by
connecting two VLEIBot prototypes in the parallel and symmetrical
configuration depicted in Fig. 6(a), using two CF beams. The separation
between the two composing bodies is chosen to be $6\,\text{mm}$ in order to
enable the cancelation, due to symmetry, of undesired torques produced by the
undulating tails during forward locomotion; this notion is illustrated in Fig.
6(b). The swimming tests corresponding to Figs. 5(c) and (d) compellingly
demonstrate an almost perfect cancelation of the opposing torques generated by
the two sides of the tested prototype. Due to this effect, VLEIBot+ prototypes
are fully controllable in the $2$D space.

Figure 6: The VLEIBot+ and block diagram of the control scheme. (a) Assembly
of two VLEIBots in parallel to create the VLEIBot+. (b) Turning functionality
of the VLEIBot+. The separation between the two propulsors of the VLEIBot+
enables a dependable production of turning torque by introducing a net force
differential. (c) Architecture used to control the VLEIBot+ during swimming.
The lateral position controller receives as input the lateral position error,
$r_{\text{e},2}=r_{\text{d},2}-r_{2}$, and computes the desired yaw angle,
$\psi_{\text{d}}$. Next, the heading controller receives as input the yaw-
angle error, $\psi_{\text{e}}=\psi_{\text{d}}-\psi$, and computes the heading
control signal, $u_{\psi}$. Last, the actuator mapping receives as inputs
$u_{\psi}$ and the constant $u_{v}$, and maps them into the duty cycles
$u_{\text{l}}$ and $u_{\text{r}}$, which are used to define the PWM signals
that excite the left and right SMA-based actuators of the robot, respectively.

###  IV-B Control Strategy

To enable the VLEIBot+ to execute turning maneuvers and track references in
the $2$D space, we implemented the control scheme in Fig. 6(c), which is
almost identical to that in [2]. Here, the pair $\left\\{r_{1},r_{2}\right\\}$
(first two components of $\boldsymbol{r}$) is the instantaneous position of
the robot’s body-fixed frame in the $2$D space relative to the inertial frame,
and $\psi$ is the instantaneous yaw angle, measured with the Vicon system. In
this scheme, the lateral position controller computes the desired heading
according to

$\displaystyle\psi_{\text{d}}=k_{\text{p},2}r_{\text{e},2}+k_{\text{i},2}\int_{0}^{t}r_{\text{e},2}\,d\tau,\vspace{-0.7ex}$
(2)

where $k_{\text{p},2}$ and $k_{\text{i},2}$ are proportional and integral
controller gains; and, $r_{\text{e},2}=r_{\text{d},2}-r_{2}$ is the lateral-
position error computed as the difference between the desired and measured
lateral positions, $r_{\text{d},2}$ and $r_{2}$, respectively. The heading
controller computes the control signal
$u_{\psi}=k_{\text{p},\psi}\psi_{\text{e}}$, where $k_{\text{p},\psi}$ is a
proportional controller gain and $\psi_{\text{e}}=\psi_{\text{d}}-\psi$ is the
yaw-angle error computed as the difference between the desired and current yaw
angles, $\psi_{\text{d}}$ and $\psi$, respectively. Last, the actuator mapping
receives $u_{\psi}$ and a constant $u_{v}$, and maps them into the duty-cycle
values $u_{\text{l}}$ and $u_{\text{r}}$ used to define the PWM signals that
respectively drive the left and right actuators of the robot, according to

$\displaystyle\left[\begin{array}[]{c}u_{\text{l}}\\\ \vspace{-2ex}\hfil\\\
u_{\text{r}}\\\
\end{array}\right]=\left[\begin{array}[]{cc}k_{\text{l}}&-k_{\text{l}}\\\
\vspace{-2ex}\hfil\\\ k_{\text{r}}&k_{\text{r}}\\\
\end{array}\right]\left[\begin{array}[]{c}u_{v}\\\ \vspace{-2ex}\hfil\\\
u_{\psi}\\\ \end{array}\right],$ (12)

where $k_{\text{l}}$ and $k_{\text{r}}$ are tunable mapping constants.

Figure 7: VLEIBot+ locomotion experiments. (a) $20$-s open-loop swimming
experiment at $5$ Hz. (i) Photographic composite of frames taken at intervals
of $5\,\text{s}$ from video footage of the VLEIBot+ swimming forward in open
loop. Both actuators of the robot were excited using a $5$-Hz PWM signal with
a duty cycle of $12\,\%$. (ii) Measured position of the VLEIBot+ over the
course of the open-loop experiment. (iii) Lateral error of the VLEIBot+ over
the course of the experiment. We measured an RMS-error value of
$20.4\,\text{mm}$. (b) $20$-s closed-loop swimming experiment at
$5\,\text{Hz}$. (i) Photographic composite of frames taken at intervals of
$5\,\text{s}$ from video footage of the VLEIBot+ swimming forward in closed
loop. (ii) Measured position of the VLEIBot+ over the course of the closed-
loop experiment. (iii) Lateral error of the VLEIBot+ over the course of the
experiment. We measured an RMS-error value of $3.94\,\text{mm}$. (c) $18$-s
closed-loop left-turn experiment at $5\,\text{Hz}$. (i) Photographic composite
of frames taken at intervals of $6\,\text{s}$ from video footage of the
VLEIBot+ swimming in closed loop. (ii) Measured position of the VLEIBot+ over
the course of the closed-loop experiment. We measure a maximum angular turning
rate of $0.25\,\text{rad}\cdot\text{s}^{-1}$. (iii) Lateral error of the
VLEIBot+ over the course of the experiment. We measured an RMS-error value of
$20.8\,\text{mm}$. (d) $18$-s closed-loop right-turn experiment at
$5\,\text{Hz}$. (i) Photographic composite of frames taken at intervals of
$6\,\text{s}$ from video footage of the VLEIBot+ swimming in closed loop. (ii)
Measured position of the VLEIBot+ over the course of the closed-loop
experiment we measure a maximum angular turning rate of
$0.28\,\text{rad}\cdot\text{s}^{-1}$. (iii) Lateral error of the VLEIBot+ over
the course of the experiment. We measured an RMS-error value of
$5.45\,\text{mm}$. Video footage of the experiments can be viewed in the
accompanying supplementary movie.

###  IV-C Swimming Experiments

To evaluate and demonstrate the swimming capabilities of the VLEIBot+, we used
the same setup depicted in Fig. 3(a) and described in Section III-A. Fig. 7(a)
summarizes the experimental data of a VLEIBot+ prototype swimming for
$20\,\text{s}$ at $5\,\text{Hz}$ in open loop, after being tuned to nominally
follow a straight path. In this test, the two robot’s actuators were excited
using PWM signals with an on voltage of $22\,\text{V}$ and duty cycle of
$12\,\%$. Fig. 7(a.i) shows a photographic composite of frames taken at
intervals of $5\,\text{s}$, of the swimming test. And, Figs. 7(a.ii) and
(a.iii) show the robot’s $2$D position and lateral-position error during the
test. The maximum average forward speed, computed over a $1$ second window, is
$13.6\,\text{mm}\cdot\text{s}^{-1}$ ($0.30\,\text{Bl}\cdot\text{s}^{-1}$). For
this experiment, the RMS value of the lateral-position error is
$20.4\,\text{mm}$. Fig. 7(b) summarizes the experimental data of a VLEIBot+
prototype swimming for $20\,\text{s}$ at $5\,\text{Hz}$ in closed loop,
following a straight-line reference and driven by the controller described in
Section IV-B. In this case, $k_{\text{p,2}}=2\,\text{rad}\cdot\text{m}^{-1}$,
$k_{\text{i,2}}=0.1\,\text{rad}\cdot\text{m}^{-1}$,
$k_{\psi}=0.2\,\text{rad}^{-1}$, $u_{v}=0.11$ (per unit), $k_{\text{l}}=1$,
and $k_{\text{r}}=1$; the on PWM voltage was set at $22\,\text{V}$. Fig.
7(b.i) shows a photographic composite of frames taken at intervals of
$5\,\text{s}$, of the swimming test. And, Figs. 7(b.ii) and (b.iii) show the
robot’s position in the $2$D space and lateral-position error during the test.
The maximum average forward speed is $14.8\,\text{mm}\cdot\text{s}^{-1}$
($0.32\,\text{Bl}\cdot\text{s}^{-1}$), and the corresponding RMS value of the
lateral-position error is $3.94\,\text{mm}$; this RMS-error value represents a
decrease of $80\,\%$ with respect to that obtained in open loop.

Fig. 7(c) summarizes the experimental data of a VLEIBot+ prototype during the
execution of a $90$-degree left-turn maneuver. Fig. 7(c.i) shows a
photographic composite of frames taken at intervals of $6\,\text{s}$, of the
swimming test. And, Figs. 7(c.ii) and (c.iii) show the robot’s position in the
$2$D space and lateral-position error during the test. The measured maximum
average speed and turning rate are $16.1\,\text{mm}\cdot\text{s}^{-1}$
($0.35\,\text{Bl}\cdot\text{s}^{-1}$) and $0.25\,\text{rad}\cdot\text{s}^{-1}$
respectively, and the corresponding RMS value of the lateral-position error is
$20.8\,\text{mm}$. Fig. 7(d) summarizes the experimental data of a VLEIBot+
prototype during the execution of a $90$-degree right-turn maneuver. Fig.
7(d.i) shows a photographic composite of frames taken at intervals of
$6\,\text{s}$, of the swimming test. And, Figs. 7(d.ii) and (d.iii) show the
robot’s position in the 2D space and lateral-position error during the test.
The measured maximum average speed and turning rate are
$10.0\,\text{mm}\cdot\text{s}^{-1}$ ($0.22\,\text{Bl}\cdot\text{s}^{-1}$) and
$0.28\,\text{rad}\cdot\text{s}^{-1}$, respectively, and the corresponding RMS
value of the lateral-position error is $5.45\,\text{mm}$.

##  V Conclusions

We presented the design and development of two new microrobotic swimmers: the
$45$-mg/$23$-$\text{mm}^{3}$ VLEIBot and its high-performance evolution, the
$90$-mg/$47$-$\text{mm}^{3}$ VLEIBot+. We designed the propulsor for the
VLEIBot drawing inspiration from the anguilliform swimming mode exhibited by
some aquatic animals such as the European eel, yellow-bellied sea snake, and
young ACF tadpole. The final heuristically-optimized propulsor design, loosely
inspired by the tail of the young ACF tadpole, has a parabolic shape, a length
of $26\,\text{mm}$, and an aspect ratio of $0.41$, and produces the best
results in terms of speed at an actuation frequency of $5\,\text{Hz}$. We
arrived to these parameters after performing hundreds of tail-characterization
experiments. With the optimized propulsor, the VLEIBot can achieve forward
speeds as high as $15.1\,\text{mm}\cdot\text{s}^{-1}$
($0.33\,\text{Bl}\cdot\text{s}^{-1}$) at a frequency of $5\,\text{Hz}$;
however, its maneuverability is extremely limited. With optimized propulsors,
the VLEIBot+ can achieve forward speeds as high as
$16.1\,\text{mm}\cdot\text{s}^{-1}$ ($0.35\,\text{Bl}\cdot\text{s}^{-1}$) at a
frequency of $5\,\text{Hz}$; furthermore, because of its high maneuverability
in the $2$D space, it can perform maneuvers at angular speeds as high as
$0.28\,\text{rad}\cdot\text{s}^{-1}$. To our best knowledge, the VLEIBot and
VLEIBot+ are the first self-propelled anguilliform-inspired microswimmers.
Furthermore, these robots are two of the very few that fully exploit FSI
phenomena to generate the waves that propagate along their tails to produce
locomotion.

## References

  * [1] X.-T. Nguyen, A. A. Calderón, A. Rigo, J. Z. Ge, and N. O. Pérez-Arancibia, “SMALLBug: A 30-mg Crawling Robot Driven by a High-Frequency Flexible SMA Microactuator,” _IEEE Robot. Automat. Lett._ , vol. 5, no. 4, pp. 6796–6803, Oct. 2020.
  * [2] R. M. Bena, X.-T. Nguyen, A. A. Calderón, A. Rigo, and N. O. Pérez-Arancibia, “SMARTI: A 60-mg Steerable Robot Driven by High-Frequency Shape-Memory Alloy Actuation,” _IEEE Robot. Automat. Lett._ , vol. 6, no. 4, pp. 8173–8180, Oct. 2021.
  * [3] C. K. Trygstad, X.-T. Nguyen, and N. O. Pérez-Arancibia, “A New $1$-mg Fast Unimorph SMA-Based Actuator for Microrobotics,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst. (IROS)_ , Detroit, MI, USA, Oct. 2023, pp. 2693–2700.
  * [4] N. O. Pérez-Arancibia, K. Y. Ma, K. C. Galloway, J. D. Greenberg, and R. J. Wood, “First controlled vertical flight of a biologically inspired microrobot,” _Bioinspir. Biomim._ , vol. 6, no. 3, Sep. 2011, Art. no. 036009.
  * [5] N. O. Pérez-Arancibia, P.-E. J. Duhamel, K. Y. Ma, and R. J. Wood, “Model-Free Control of a Hovering Flapping-Wing Microrobot,” _J. Intell. Robot. Syst._ , vol. 77, no. 1, pp. 95–111, Jan. 2015.
  * [6] X. Yang, Y. Chen, L. Chang, A. A. Calderón, and N. O. Pérez-Arancibia, “Bee+: A 95-mg Four-Winged Insect-Scale Flying Robot Driven by Twinned Unimorph Actuators,” _IEEE Robot. Automat. Lett._ , vol. 4, no. 4, pp. 4270–4277, Oct. 2019.
  * [7] R. M. Bena, X.-T. Nguyen, X. Yang, A. A. Calderón, Y. Chen, and N. O. Pérez-Arancibia, “A Multiplatform Position Control Scheme for Flying Robotic Insects,” _J. Intell. Robot. Syst._ , vol. 105, no. 1, May 2022, Art. no. 19.
  * [8] R. M. Bena, X. Yang, A. A. Calderón, and N. O. Pérez-Arancibia, “High-Performance Six-DOF Flight Control of the Bee++: An Inclined-Stroke-Plane Approach,” _IEEE Trans. Robot._ , vol. 39, no. 2, pp. 1668–1684, Apr. 2023.
  * [9] X. Yang, L. Chang, and N. O. Pérez-Arancibia, “An 88-milligram insect-scale autonomous crawling robot driven by a catalytic artificial muscle,” _Sci. Robot._ , vol. 5, no. 45, Aug. 2020, Art. no. eaba0015.
  * [10] S. B. Fuller, Z. E. Teoh, P. Chirattananon, N. O. Pérez-Arancibia, J. Greenberg, and R. J. Wood, “Stabilizing air dampers for hovering aerial robotics: design, insect-scale flight tests, and scaling,” _Auton. Robot._ , vol. 41, no. 8, pp. 1555–1573, Dec. 2017.
  * [11] N. O. Pérez-Arancibia, P.-E. J. Duhamel, K. Y. Ma, and R. J. Wood, “Model-Free Control of a Hovering Flapping-Wing Microrobot,” _J. Intell. Robot. Syst._ , vol. 77, no. 1, pp. 95–111, Jan. 2015.
  * [12] ——, “Model-Free Control of a Flapping-Wing Flying Microrobot,” in _Proc. 16th Int. Conf. Adv. Robot. (ICAR)_ , Montevideo, Uruguay, Nov. 2013, pp. 1–8.
  * [13] P.-E. J. Duhamel, N. O. Pérez-Arancibia, G. L. Barrows, and R. J. Wood, “Biologically Inspired Optical-Flow Sensing for Altitude Control of Flapping-Wing Microrobots,” _IEEE/ASME Trans. Mechatron._ , vol. 18, no. 2, pp. 556–568, Apr. 2013.
  * [14] Z. E. Teoh, S. B. Fuller, P. Chirarttananon, N. O. Pérez-Arancibia, J. D. Greenberg, and R. J. Wood, “A Hovering Flapping-Wing Microrobot with Altitude Control and Passive Upright Stability,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst. (IROS)_ , Vilamoura, Algarve, Portugal, Oct. 2012, pp. 3209–3216.
  * [15] P.-E. J. Duhamel, N. O. Pérez-Arancibia, G. L. Barrows, and R. J. Wood, “Altitude Feedback Control of a Flapping-Wing Microrobot Using an On-Board Biologically Inspired Optical Flow Sensor,” in _Proc. IEEE Int. Conf. Robot. Autom. (ICRA)_ , Saint Paul, Minnesota, USA, May 2012, pp. 4228–4235.
  * [16] N. O. Pérez-Arancibia, J. P. Whitney, and R. J. Wood, “Lift Force Control of a Flapping-Wing Microrobot,” in _Proc. Amer. Control Conf. (ACC)_ , San Francisco, CA, USA, Jun. 2011, pp. 4761–4768.
  * [17] N. O. Pérez-Arancibia, P. Chirarattananon, B. M. Finio, and R. J. Wood, “Pitch-Angle Feedback Control of a Biologically Inspired Flapping-Wing Microrobot,” in _Proc. Int. Conf. Robot. Biomim. (ROBIO)_ , Karon Beach, Thailand, Dec. 2011, pp. 1495–1502.
  * [18] N. O. Pérez-Arancibia, J. P. Whitney, and R. J. Wood, “Lift Force Control of Flapping-Wing Microrobots Using Adaptive Feedforward Schemes,” _IEEE/ASME Trans. Mechatron._ , vol. 18, no. 1, pp. 155–168, Feb. 2013.
  * [19] A. P. Palstra and G. E. E. J. M. van den Thillart, “Swimming physiology of European silver eels (Anguilla anguilla L.): energetic costs and effects on sexual maturation and reproduction,” _Fish Physiol. Biochem._ , vol. 36, no. 3, pp. 297–322, Apr. 2010.
  * [20] J. B. Graham, W. R. Lowell, I. Rubinoff, and J. Motta, “Surface and Subsurface Swimming of the Sea Snake Pelamis Platurus,” _J. Exp. Biol._ , vol. 127, no. 1, pp. 27–44, Jan. 1986.
  * [21] A. Roberts, N. A. Hill, and R. Hicks, “Simple Mechanisms Organise Orientation of Escape Swimming in Embryos and Hatchling Tadpoles of Xenopus Laevis,” _J. Exp. Biol._ , vol. 203, no. 12, pp. 1869–1885, Jun. 2000.
  * [22] N. B. Tack, K. T. Du Clos, and B. J. Gemmell, “Anguilliform Locomotion across a Natural Range of Swimming Speeds,” _Fluids_ , vol. 6, no. 3, Mar. 2021, Art. no. 127.
  * [23] H. Liu, R. Wassersug, and K. Kawachi, “The Three-Dimensional Hydrodynamics of Tadpole Locomotion,” _J. Exp. Biol._ , vol. 200, no. 22, pp. 2807–2819, Nov. 1997.
  * [24] T. Fukuda, H. Hosokai, and I. Kikuchi, “Distributed Type of Actuators by Shape Memory Alloy and its Application to Underwater Mobile Robotic Mechanism,” in _Proc. IEEE Int. Conf. Robot. Autom. (ICRA)_ , Cincinnati, OH, USA, May 1990, pp. 1316–1321.
  * [25] L. J. Garner, L. N. Wilson, D. C. Lagoudas, and O. K. Rediniotis, “Development of a shape memory alloy actuated biomimetic vehicle,” _Smart Mater. Struct._ , vol. 9, no. 5, pp. 673–683, Oct. 2000.
  * [26] C. Rossi, J. Colorado, W. Coral, and A. Barrientos, “Bending continuous structures with SMAs: a novel robotic fish design,” _Bioinsp. Biomim._ , vol. 6, no. 4, Dec. 2011, Art. no. 045005.
  * [27] K.-J. Cho, E. Hawkes, C. Quinn, and R. J. Wood, “Design, fabrication and analysis of a body-caudal fin propulsion system for a microrobotic fish,” in _Proc. IEEE Int. Conf. Robot. Autom. (ICRA)_ , Pasadena, CA, USA, May 2008, pp. 706–711.
  * [28] Z. Wang, G. Hang, J. Li, Y. Wang, and K. Xiao, “A micro-robot fish with embedded SMA wire actuated flexible biomimetic fin,” _Sens. Actuators A: Phys._ , vol. 144, no. 2, pp. 354–360, Jun. 2008.
  * [29] L. Shi, S. Guo, and K. Asaka, “A Novel Jellyfish-like Biomimetic Microrobot,” in _Proc. IEEE/ICME Int. Conf. Comp. Med. Eng. (CME)_ , Gold Coast, QLD, Australia, Jul. 2010, pp. 277–281.
  * [30] D. Byun, J. Choi, K. Cha, J.-O. Park, and S. Park, “Swimming microrobot actuated by two pairs of Helmholtz coils system,” _Mechatronics_ , vol. 21, no. 1, pp. 357–364, Feb. 2011.
  * [31] Q. Ze, S. Wu, J. Dai, S. Leanza, G. Ikeda, P. C. Yang, G. Iaccarino, and R. R. Zhao, “Spinning-enabled wireless amphibious origami millirobot,” _Nat. Commun._ , vol. 13, Jun. 2022, Art. no. 3118.
  * [32] J. Chen, H. Hu, and Y. Wang, “Magnetic-driven 3D-printed biodegradable swimming microrobots,” _Smart Mater. Struct._ , vol. 32, no. 8, Jul. 2023, Art. no. 085014.
  * [33] L. Tan and D. J. Cappelleri, “Design, Fabrication, and Characterization of a Helical Adaptive Multi-Material MicroRobot (HAMMR),” _IEEE Robot. Automat. Lett._ , vol. 8, no. 3, pp. 1723–1730, Mar. 2023.
  * [34] S. N. Tabatabaei, J. Lapointe, and S. Martel, “Shrinkable Hydrogel-Based Magnetic Microrobots for Interventions in the Vascular Network,” _Adv. Robot._ , vol. 25, no. 8, pp. 1049–1067, Apr. 2011.
  * [35] S. Fusco, M. S. Sakar, S. Kennedy, C. Peters, R. Bottani, F. Starsich, A. Mao, G. A. Sotiriou, S. Pané, S. E. Pratsinis, D. Mooney, and B. J. Nelson, “An integrated microrobotic platform for on-demand, targeted therapeutic interventions,” _Adv. Mater._ , vol. 26, no. 6, pp. 952–957, Feb. 2014.
  * [36] H. Li, J. Tan, and M. Zhang, “Dynamics Modeling and Analysis of a Swimming Microrobot for Controlled Drug Delivery,” _IEEE Trans. Autom. Sci. Eng._ , vol. 6, no. 2, pp. 220–227, Apr. 2009.
  * [37] F. Z. Temel and S. Yesilyurt, “Confined swimming of bio-inspired microrobots in rectangular channels,” _Bioinspir. Biomim._ , vol. 10, no. 1, Feb. 2015, Art. no. 016015.
  * [38] S. Palagi, V. Pensabene, L. Beccai, B. Mazzolai, A. Menciassi, and P. Dario, “Design and development of a soft magnetically-propelled swimming microrobot,” in _Proc. IEEE Int. Conf. Robot. Autom. (ICRA)_ , Shanghai, China, May 2011, pp. 5109–5114.
  * [39] O. S. Pak, W. Gao, J. Wang, and E. Lauga, “High-speed propulsion of flexible nanowire motors: Theory and experiments,” _Soft Matter_ , vol. 7, no. 18, Sep. 2011, Art. no. 8169.
  * [40] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, “Microscopic artificial swimmers,” _Nature_ , vol. 437, no. 6, pp. 862–865, Oct. 2005.
  * [41] W. Gao, S. Sattayasamitsathit, K. M. Manesh, D. Weihs, and J. Wang, “Magnetically Powered Flexible Metal Nanowire Motors,” _J. Am. Chem. Soc._ , vol. 132, no. 41, pp. 14 403–14 405, Sep. 2010.
  * [42] T. Xu, G. Hwang, N. Andreff, and S. Régnier, “The rotational propulsion characteristics of scaled-up helical microswimmers with different heads and magnetic positioning,” in _Proc. IEEE/ASME Int. Conf. Adv. Intell. Mechatron. (AIM)_ , Wollongong, NSW, Australia, Jul. 2013, pp. 1114–1120.
  * [43] A. Ghosh and P. Fischer, “Controlled Propulsion of Artificial Magnetic Nanostructured Propellers,” _Nano Lett._ , vol. 9, no. 6, pp. 2243–2245, May 2009.
  * [44] C. E. Sing, L. Schmid, M. F. Schneider, and A. Alexander-Katz, “Controlled surface-induced flows from the motion of self-assembled colloidal walkers,” _Proc. Natl. Acad. Sci._ , vol. 107, no. 2, pp. 535–540, Dec. 2009.
  * [45] L. Zhang, T. Petit, Y. Lu, B. E. Kratochvil, K. E. Peyer, R. Pei, J. Lou, and B. J. Nelson, “Controlled Propulsion and Cargo Transport of Rotating Nickle Nanowires near a Patterned Solid Surface,” _ACS Nano_ , vol. 4, no. 10, pp. 6228–6234, Sep. 2010.
  * [46] W. Liu, X. Jia, F. Wang, and Z. Jia, “An in-pipe wireless swimming microrobot driven by giant magnetostrictive film,” _Sens. Actuators A: Phys._ , vol. 160, no. 2, pp. 101–108, May 2010.
  * [47] P. Tierno, R. Golestanian, I. Pagonabarraga, and F. Sagués, “Magnetically Actuated Colloidal Microswimmers,” _J. Phys. Chem. B_ , vol. 112, no. 51, pp. 16 525–16 528, Nov. 2008.
  * [48] S. Floyd, C. Pawashe, and M. Sitti, “An Untethered Magnetically Actuated Micro-Robot Capable of Motion on Arbitrary Surfaces,” in _Proc. IEEE/ASME Int. Conf. Robot. Autom. (ICRA)_ , Pasadena, CA. USA, May. 2008, pp. 419–424.
  * [49] Q. Zhao, S. Liu, J. Chen, G. He, J. Di, L. Zhao, T. Su, M. Zhang, and Z. Hou, “Fast-moving piezoelectric micro-robotic fish with double caudal fins,” _Robot. Auton. Syst._ , vol. 140, Jun. 2021, Art. no 103733.
  * [50] K. Li, X. Zhou, Y. Liu, J. Sun, X. Tian, H. Zheng, . Zhang, J. Deng, J. Liu, W. Chen, and J. Zhao, “A 5 cm-scale piezoelectric jetting agile underwater robot,” _Adv. Intell. Syst._ , vol. 5, no. 4, Apr. 2023, Art. no 2200262.
  * [51] Y. S. Song and M. Sitti, “Surface-Tension-Driven Biologically Inspired Water Strider Robots: Theory and Experiments,” _IEEE Trans. Robot._ , vol. 23, no. 3, pp. 578–589, Jun. 2007.
  * [52] Y. Chen, N. Doshi, B. Goldberg, H. Wang, and R. J. Wood, “Controllable water surface to underwater transition through electrowetting in a hybrid terrestrial-aquatic microrobot,” _Nature Commu._ , vol. 9, Jun. 2018, art. no. 2495.
  * [53] X. Deng and S. Avadhanula, “Biomimetic Micro Underwater Vehicle with Oscillating Fin Propulsion: System Design and Force Measurement,” in _Proc. IEEE/ASME Int. Conf. Robot. Autom. (ICRA)_ , Barcelona, Spain, Apr. 2005, pp. 3312–3317.
  * [54] A. Ming, S. Park, Y. Nagata, and M. Shimojo, “Development of Underwater Robots using Piezoelectric Fiber Composite,” in _Proc. IEEE/ASME Int. Conf. Robot. Autom. (ICRA)_ , Kobe, Japan, May 2009, pp. 3821–3826.
  * [55] T. Fukuda, A. Kawamoto, F. Arai, and H. Matsuura, “Mechanism and Swimming Experiment of Micro Mobile Robot in Water,” in _Proc. IEEE/ASME Int. Conf. Robot. Autom. (ICRA)_ , San Diego, CA, USA, May 1994, pp. 814–819.
  * [56] B. Kim, D.-H. Kim, J. Jung, and J.-O. Park, “A biomimetic undulatory tadpole robot using ionic polymer–metal composite actuators,” _Smart Mater. Struct._ , vol. 14, no. 6, pp. 1579–1585, Dec. 2005.
  * [57] N. Kamamichi, M. Yamakita, K. Asaka, and Z.-W. Luo, “A Snake-like Swimming Robot Using IPMC Actuator/Sensor,” in _Proc. IEEE/ASME Int. Conf. Robot. Autom. (ICRA)_ , Orlando, Florida, May 2006, pp. 1812–1817.
  * [58] S. Guo, T. Fukuda, and K. Asaka, “A New Type of Fish-Like Underwater Microrobot,” _IEEE/ASME Trans. Robot._ , vol. 8, no. 1, pp. 136–141, Mar. 2003.
  * [59] E. M. Purcell, “Life at low Reynolds number,” _Am. J. Phys._ , vol. 45, no. 1, pp. 3–11, Jan. 1977.
  * [60] B. M. Finio, N. O. Pérez-Arancibia, and R. J. Wood, “System identification and linear time-invariant modeling of an insect-sized flapping-wing micro air vehicle,” in _Proc. IEEE Int. Conf. Intell. Robots Syst. (IROS)_ , San Francisco, CA. USA, Sep. 2011, pp. 1107–1114.