# A geometrically exact model for thin magneto-elastic shells

Matteo Pezzulla111Current affiliation: Slender Structures Lab, Department of
Mechanical and Production Engineering, Århus University, Inge Lehmanns Gade
10, 8000 Århus C, Denmark<EMAIL_ADDRESS>Dong Yan<EMAIL_ADDRESS>Pedro M.
Reis<EMAIL_ADDRESS>Flexible Structures Laboratory, Institute of
Mechanical Engineering, École Polytechnique Fédérale de Lausanne, 1015
Lausanne, Switzerland

###### Abstract

We develop a reduced model for hard-magnetic, thin, linear-elastic shells that
can be actuated through an external magnetic field, with geometrically exact
strain measures. Assuming a reduced kinematics based on the Kirchhoff-Love
assumption, we derive a reduced two-dimensional magneto-elastic energy that
can be minimized through numerical analysis. In parallel, we simplify the
reduced energy by expanding it up to the second order in the displacement
field and provide a physical interpretation. Our theoretical analysis allows
us to identify and interpret the two primary mechanisms dictating the magneto-
elastic response: a combination of equivalent magnetic pressure and forces at
the first order, and distributed magnetic torques at the second order. We
contrast our reduced framework against a three-dimensional nonlinear model by
investigating three test cases involving the indentation and the pressure
buckling of shells under magnetic loading. We find excellent agreement between
the two approaches, thereby verifying our reduced model for shells undergoing
nonlinear and non-axisymmetric deformations. We believe that our model for
magneto-elastic shells will serve as a valuable tool for the rational design
of magnetic structures, enriching the set of reduced magnetic models.

###### keywords:

Magneto-rheological elastomer , Slender structures , Shells , Buckling

##  1 Introduction

Investigating the effect of non-mechanical stimuli on structures has long been
the subject of many research efforts, from the bending of bi-metallic
thermostats (Timoshenko, 1925) to the swelling-induced morphing of thin
structures (Kim et al., 2012). While these stimuli can induce extreme
deformations, the time scales of the actuation processes can be relatively
large, in the order of some seconds or minutes, a feature that might not be
ideal in tasks where fast actuation is required. By contrast, the coupling
between magnetism and elasticity represents a valuable avenue towards fast and
reversible actuation of soft structures (Lum et al., 2016; Kim et al., 2018a;
Zhao et al., 2019; Wang et al., 2020). A promising class of materials for such
actuations is that of magneto-rheological elastomers (MREs), which are
composites of magnetized (metallic) particles and a soft elastic polymeric
matrix that respond to an external applied magnetic field. Applications of
MREs range from minimal invasive procedures (Pancaldi et al., 2020) to soft
robotic actuators (Hu et al., 2018; Gu et al., 2020; Alapan et al., 2020) and
bio-medical devices (Kim et al., 2019).

Depending on the extent to which MREs can keep their residual magnetization
and how they respond to an external magnetic field, these materials can be
classified into different categories, such as superparamagnetic and soft-
ferromagnetic materials. For a more detailed discussion on these different
categories, we refer to (Bertotti, 1998; Sano et al., 2021). Past pioneering
studies have primarily focused on the deformation of structures made of
superparamagnetic (Moon and Pao, 1968; Cēbers, 2003; Cēbers and Javaitis,
2004; Cēbers and Cīrulis, 2007; Dreyfus et al., 2005; Roper et al., 2006;
Gerbal et al., 2015) or soft-ferromagnetic (Rigbi and Jilkén, 1983; Ginder et
al., 1999; Dorfmann and Ogden, 2003; Danas et al., 2012; Loukaides et al.,
2014; Seffen and Vidoli, 2016; Psarra et al., 2019) materials under an
external magnetic field. More recent efforts have turned to hard-magnetic soft
materials (Lum et al., 2016; Kim et al., 2018a; Zhao et al., 2019; Wang et
al., 2020; Yan et al., 2020). These materials are MREs made of hard-
ferromagnetic particles embedded in a soft elastomeric matrix. They retain a
permanent magnetization, have high coercivity (_i.e._ the necessary field
strength to erase the magnetization) upon saturation, while being mechanically
compliant due to the soft elastomeric matrix and their usual slender geometry
utilized in applications. Therefore, hard MREs present several advantages for
tasks where a fast and reversible actuation is required, as for example in
soft robotics, minimal invasive procedures, and bio-medical devices.

During the past decade, there have been several efforts on the modeling of
these materials, motivated by the numerous possible applications of hard MREs
(Danas et al., 2012). As an example, a nonlinear three-dimensional (3D) theory
for hard MREs was proposed by Zhao et al. (2019), where the Helmholtz free
energy of the system consists of elastic (neo-Hookean) and magneto-elastic
contributions. The basic assumption, following the physical observations in
(Bertotti, 1998), is that the magnetic flux density is linear with the applied
magnetic field. The model was then tested against experimental results, by
implementing the theoretical model in a finite element scheme in Abaqus,
resulting in a quantitative agreement between the two.

Based on the 3D model presented in (Zhao et al., 2019), reduced theories for
hard magnetic linear and nonlinear beams (Wang et al., 2020; Yan et al.,
2021), and rods (Sano et al., 2021) have been derived. In these studies, a
dimensional reduction procedure was performed on the 3D magneto-elastic
energy, assuming a reduced kinematics for the beams and the rods based on the
Kirchhoff-Love assumptions. Moreover, we have recently presented a study on
magneto-active axisymmetric shells made of hard MREs in (Yan et al., 2020),
where the coupling between mechanics and magnetism was leveraged to change the
stability onsets of shells undergoing pressure buckling (Hutchinson and
Thompson, 2018; Lee et al., 2016, 2019). In this study, experiments were
contrasted with a magnetic shell model for axisymmetric deformation and
geometrically exact strain measures (Yan et al., 2020). The results show that
the magnetic field can be used to tune the critical buckling pressure of
spherical shells, which are highly sensitive to imperfections (Hutchinson et
al., 1971; Hutchinson, 2016; Hutchinson and Thompson, 2018). Theoretical
models for non-axisymmetric deformations of magneto-active shells have
received less attention and have been derived for shallow shells only (Seffen
and Vidoli, 2016; Loukaides et al., 2014), with an understanding of the
coupling between magnetism and mechanics in the general case that remains
unclear. In our view, deriving reduced-order structural models should have the
twofold goal of providing alternatives to 3D models that can be solved
numerically with a reduced computational cost, and enabling a better
understanding of the mechanics of the problem given the lower complexity of
the reduced models.

In this paper, we derive a theory for thin, elastic, magnetic shells with
geometrically exact strain measures (Babcock, 1983; Koiter, 1969; Niordson,
1985), thereby generalizing the model that we presented in (Yan et al., 2020)
that was only limited to axisymmetric shells undergoing axisymmetric
deformations. A model with even simpler kinematics, such as with moderate
rotations (Sanders, 1963; Donnell, 1977), could also be derived to further
reduce the complexity of the model, but this derivation is beyond the scope of
the paper as we want to focus on the more general formulation with geometrical
exact strain measures. In particular, we perform a dimensional reduction of
the 3D magneto-elastic energy contribution presented in (Zhao et al., 2019),
by assuming a reduced kinematics according to the Kirchhoff-Love assumptions
(Niordson, 1985). The result is a reduced bidimensional (2D) energy that can
be minimized via a finite element scheme implemented in COMSOL Multiphysics, a
commercial software. Moreover, given that the reduced magnetic energy is
highly nonlinear and not amenable to physical interpretation, we expand it up
to the second order in the displacement field and untangle two different
contributions with a clear physical meaning. At the first order in the
displacement field, the magneto-elastic coupling can be interpreted as a
combination of equivalent magnetic pressure and in-plane forces. At the second
order, the magneto-elastic coupling is instead represented by distributed
torques. To verify our reduced model, we test it against the existing 3D model
derived by (Zhao et al., 2019), implemented in Abaqus, finding excellent
agreement in a set of three different test cases that we investigate in
detail. Specifically, we test the model in the cases of (i) non-axisymmetric
indentation under a magnetic field; (ii) pressure buckling under asymmetric
magnetic loading where the residual magnetization vector and the external
magnetic field are in the same plane; and (=iii) pressure buckling under
asymmetric magnetic loading, where the residual magnetization vector and the
external magnetic field are not in the same plane.

Our paper is organized as follows. In Sec. 2, we set the notation, describe
the geometry of the shell, and recall the reduced elastic energy of its mid-
surface. We perform the dimensional reduction on the magnetic energy of the
shell in Sec. 3. Then, in Sec. 4, we provide a mechanical interpretation to
the reduced magnetic energy. The numerical implementation of our model is
summarized in Sec. 5, whereas the validation with the three different problems
is presented in Sec. 6. In Sec. 7, we provide final discussions where we also
summarize our findings.

##  2 Preliminaries

For the sake of completeness and convenience to the reader, we start by
recalling some standard concepts of shell theory (Niordson, 1985) and
differential geometry (O’Neill, 1997; do Carmo, 2016) that we will need for
the analysis presented in this paper, since this material is often difficult
to find in a synthesized way in the vast and fragmented literature of shell
mechanics. We will also recall the reduced energy for linearly elastic thin
shells, with geometrically exact strain measures.

For the rest of the Section, we will follow mainly two references: the book on
shell theory by Niordson (1985), and the book on elementary differential
geometry by O’Neill (1997). When dealing with entities from 3D continuum
mechanics, we will follow the notation of Gurtin et al. (2010).

###  2.1 Geometry of the shell

A shell is a 3D body, $\mathcal{B}$, embedded in an Euclidean space,
$\mathcal{E}$, with Cartesian basis
$(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3})$. One of the characteristic
dimensions of a shell (the thickness) is much smaller than the other two
(along the surface), as depicted in the schematic diagram of Fig. 1. Shells
are typically described by their mid-surface $\mathcal{S}\subset\mathcal{E}$,
curved in its natural (stress-free) state. We denote the parametrization of
the mid-surface as
$\accentset{\circ}{\mathbf{r}}(\eta^{1},\eta^{2})\colon\mathcal{D}\rightarrow\mathcal{E}$.
Here, $(\eta^{1},\eta^{2})$ are curvilinear coordinates, $\mathcal{D}$ is the
domain of parametrization, and the accent $\accentset{\circ}{()}$ denotes
quantities in the undeformed configuration. As a 3D body, the shell is a stack
of surfaces (O’Neill, 1997), viewed as the Cartesian product
$\mathcal{B}=\mathcal{S}\times[-h/2,h/2]$, where $h(\eta^{1},\eta^{2})$ is the
thickness of the shell, which is, in general, a function of the curvilinear
coordinates.

Figure 1: Schematic representation of a non-axisymmetric shell, in its
reference configuration. The covariant base vectors $\mathbf{G}_{i}$ are drawn
at one point of the top surface, with the tangent plane. The thickness $h$ and
the mid-surface $\mathcal{S}$ are also shown.

Once the parametrization of the mid-surface is given, the covariant base
vectors associated to the surface can be computed as
$\accentset{\circ}{\mathbf{a}}_{\alpha}=\accentset{\circ}{\mathbf{r}},_{\alpha}$,
where $(),_{\alpha}$ denotes partial differentiation with respect to
$\eta^{\alpha}$, and Greek indices run from $1$ to $2$ (O’Neill, 1997). The
vectors $\accentset{\circ}{\mathbf{a}}_{\alpha}$ give rise to a covariant
basis for the mid-surface $\mathcal{S}$, whose covariant metric coefficients
can be determined as
$\accentset{\circ}{a}_{\alpha\beta}=\accentset{\circ}{\mathbf{a}}_{\alpha}\cdot\accentset{\circ}{\mathbf{a}}_{\beta}$,
also referred to as first fundamental form, where $(\cdot)$ denotes the inner
product in the Euclidean space $\mathcal{E}$. The contravariant metric
coefficients can be determined as the components of the inverse matrix
$\accentset{\circ}{a}^{\alpha\beta}=(\accentset{\circ}{a}_{\alpha\beta})^{-1}$,
so that the contravariant base vectors can be computed as
$\accentset{\circ}{\mathbf{a}}^{\alpha}=\accentset{\circ}{a}^{\alpha\beta}\accentset{\circ}{\mathbf{a}}_{\beta}$
(O’Neill, 1997). Finally, the normal unit vector to the mid-surface is
computed as
$\accentset{\circ}{\mathbf{n}}=(\accentset{\circ}{\mathbf{a}}_{1}\times\accentset{\circ}{\mathbf{a}}_{2})/|\accentset{\circ}{\mathbf{a}}_{1}\times\accentset{\circ}{\mathbf{a}}_{2}|$,
ensuring its perpendicularity to any tangent plane along the mid-surface.
Furthermore, the covariant components of the curvature tensor are defined as
$\accentset{\circ}{b}_{\alpha\beta}=\accentset{\circ}{\mathbf{n}}\cdot\accentset{\circ}{\mathbf{r}},_{\alpha\beta}$
(also referred to as second fundamental form).

With the definition of the normal vector $\accentset{\circ}{\mathbf{n}}$, the
parametrization of the shell body $\mathcal{B}$ can be written as
$\mathbf{X}=\accentset{\circ}{\mathbf{r}}+\eta^{3}\accentset{\circ}{\mathbf{n}}$,
where $\eta^{3}$ is the coordinate running along the normal vector. Then, with
Latin indices run from $1$ to $3$, the covariant base vectors associated with
$\mathbf{X}$ are $\mathbf{G}_{i}=\mathbf{X},_{i}$. By applying the definition
of the covariant base vectors $\mathbf{G}_{i}=\mathbf{X},_{i}$ to the
parametrization $\mathbf{X}$, one can derive

$\displaystyle\mathbf{G}_{\alpha}=\accentset{\circ}{\mathbf{a}}_{\alpha}+\eta^{3}\accentset{\circ}{\mathbf{n}},_{\alpha}\,,$
(1) $\displaystyle\mathbf{G}_{3}=\accentset{\circ}{\mathbf{n}}\,.$

Similarly to what was derived for the mid-surface (_i.e._ the first and second
fundamental forms $\accentset{\circ}{a}_{\alpha\beta}$ and
$\accentset{\circ}{b}_{\alpha\beta}$), the 3D covariant metric coefficients
can be computed as $G_{ij}=\mathbf{G}_{i}\cdot\mathbf{G}_{j}$, and the
contravariant metric coefficients as $G^{ij}=(G_{ij})^{-1}$. Then, the
contravariant base vectors can be computed as
$\mathbf{G}^{i}=G^{ij}\mathbf{G}_{j}$.

Finally, when dealing with a stack of surfaces, it is natural to ask how the
area measure varies along the thickness coordinate $\eta^{3}$. Not only is
this a remarkable result in differential geometry, but is also relevant to any
recipe for dimensional reduction involving shells. Defining
$\accentset{\circ}{G}=\sqrt{\det(G_{ij})}$ and
$\accentset{\circ}{a}=\sqrt{\det(\accentset{\circ}{a}_{\alpha\beta})}$, it can
be shown that (O’Neill, 1997)

$\accentset{\circ}{G}=(1-2\eta^{3}\accentset{\circ}{H}+(\eta^{3})^{2}\accentset{\circ}{K})\accentset{\circ}{a}\,,$
(2)

where $\accentset{\circ}{H}$ and $\accentset{\circ}{K}$ are the mean and
Gaussian curvatures of the undeformed mid-surface, respectively. This
geometric result explains how local areas vary within the stack of surfaces,
when moving along the normal. For thin shells, that is at first order in
$\eta^{3}$, changes in local areas are dictated by the mean curvature only,
which is why minimal surfaces have zero mean curvature (Deserno, 2004).

When the shell deforms in response to mechanical loads (e.g., distributed
forces and pressures) or non-mechanical loads (e.g., swelling or magnetic
fields), its material points are displaced, forming a new configuration,
termed deformed or actual configuration. The positions of the material points
in this new configuration are denoted as $\mathbf{x}$, as opposed to
$\mathbf{X}$ for the undeformed configuration. For thin shells, a reduced
kinematics is usually assumed following the Kirchhoff-Love assumptions,
according to which fibers normal to the mid-surface do not stretch nor shrink,
and remain normal to the mid-surface, upon deformation. Assuming such a
reduced kinematics means that the 3D deformed shape of the shell can be fully
described by the sole shape of its mid-surface.

All the quantities that have been computed in the undeformed state such as the
base vectors $\mathbf{G}_{i}$ and the fundamental forms
$\accentset{\circ}{a}_{\alpha\beta}$ and $\accentset{\circ}{b}_{\alpha\beta}$,
defined by the parametrization $\mathbf{X}$, can be now expressed in the
deformed state, $\mathbf{x}$. The parametrization of the deformed mid-surface
is denoted as $\mathbf{r}=\accentset{\circ}{\mathbf{r}}+\mathbf{u}$, where
$\mathbf{u}$ is the displacement of each material particle of the mid-surface,
from the undeformed to the deformed configurations. The displacement field is
usually expressed in the covariant basis as
$\mathbf{u}=u^{\alpha}\accentset{\circ}{\mathbf{a}}_{\alpha}+w\accentset{\circ}{\mathbf{n}}$,
where Einstein’s summation convention is implied and
$\check{\mathbf{u}}=u^{\alpha}\accentset{\circ}{\mathbf{a}}_{\alpha}$ is the
in-plane displacement vector field. The covariant base vectors of the deformed
mid-surface are $\mathbf{a}_{\alpha}=\mathbf{r},_{\alpha}$, whereas those
associated with $\mathbf{x}$ are $\mathbf{g}_{i}=\mathbf{x},_{i}$. The normal
vector is
$\mathbf{n}=(\mathbf{a}_{1}\times\mathbf{a}_{2})/|\mathbf{a}_{1}\times\mathbf{a}_{2}|$.
Finally, the first and second fundamental forms associated with the deformed
mid-surface are, respectively,
$a_{\alpha\beta}=\mathbf{a}_{\alpha}\cdot\mathbf{a}_{\beta}$ and
$b_{\alpha\beta}=\mathbf{n}\cdot\mathbf{r},_{\alpha\beta}$. Similarly to what
derived for the undeformed shell in Eq. (1), the base vectors $\mathbf{g}_{i}$
can be expressed as

$\displaystyle\mathbf{g}_{\alpha}=\mathbf{a}_{\alpha}+\eta^{3}\mathbf{n},_{\alpha}\,,$
(3) $\displaystyle\mathbf{g}_{3}=\mathbf{n}\,.$

With classic formalism describing the geometry of the shell presented above,
we can now move on to specify the deformation gradient and its expression for
thin shells.

###  2.2 The deformation gradient

As it will be made clear in Sec. 3, the deformation gradient $\mathbf{F}$ has
to be expressed in terms of the quantities describing the shell kinematics to
be able to perform the dimensional reduction.

The deformation gradient is a two-point tensor that maps the undeformed
tangent space to the deformed tangent space, defined as (Gurtin et al., 2010)

$\mathbf{F}=\mathbf{g}_{i}\otimes\mathbf{G}^{i}=\mathbf{g}_{\alpha}\otimes\mathbf{G}^{\alpha}+\mathbf{n}\otimes\accentset{\circ}{\mathbf{n}}\,.$
(4)

Before we can express $\mathbf{F}$ only in terms of surface quantities, we
first need to quantify how the contravariant components of the 3D metric
$G^{ij}$ vary as a function of the thickness coordinate $\eta^{3}$. Expressing
$G^{ij}$ as a function of $\eta^{3}$ is a common procedure for plates, where
the deformation gradient is a linear function of the thickness coordinate
(Pietraszkiewicz, 1980), whereas it is more intricate for shells because of
the existing rest curvature.

Before we can compute the contravariant metric coefficients of the undeformed
shell $G^{\alpha\beta}$, we first compute the covariant components as
(O’Neill, 1997)

$G_{\alpha\beta}=\mathbf{G}_{\alpha}\cdot\mathbf{G}_{\beta}=\accentset{\circ}{a}_{\alpha\beta}-2\eta^{3}\accentset{\circ}{b}_{\alpha\beta}+(\eta^{3})^{2}\accentset{\circ}{c}_{\alpha\beta}\,,$
(5)

where
$\accentset{\circ}{c}_{\alpha\beta}=\accentset{\circ}{\mathbf{n}},_{\alpha}\cdot\accentset{\circ}{\mathbf{n}},_{\beta}$
is the third fundamental form of the undeformed mid-surface. Recalling that,
by definition, $\mathbf{G}_{3}\cdot\mathbf{G}_{\alpha}=0$ and that
$\mathbf{G}_{3}\cdot\mathbf{G}_{3}=1$, the covariant 3D metric coefficients of
the undeformed shell are

$G_{ij}=\begin{pmatrix}G_{\alpha\beta}&\begin{matrix}0\\\ 0\end{matrix}\\\
\begin{matrix}0&0\end{matrix}&\\!\\!1\end{pmatrix}\,.$ (6)

The structure of the metric conveys that the normal direction, assumed to be a
coordinate line, is indeed normal to the mid-surface. Since $G_{\alpha 3}=0$,
we have that $G^{\alpha 3}=0$ and $G^{\alpha\beta}=(G_{\alpha\beta})^{-1}$. To
perform this inversion and simplify the algebra, we employ matrix notation,
where the covariant components of the first, second, and third fundamental
forms are denoted as $\mathbb{A}$, $\mathbb{B}$ and $\mathbb{C}$ so that the
covariant 3D metric may be expressed as

$\mathbb{G}=\mathbb{A}-2\eta^{3}\mathbb{B}+(\eta^{3})^{2}\mathbb{C}=\mathbb{A}(\mathbb{I}-2\eta^{3}\mathbb{A}^{-1}\mathbb{B}+(\eta^{3})^{2}\mathbb{A}^{-1}\mathbb{C})\,,$
(7)

where $\mathbb{I}$ is the identity matrix. We compute the inverse matrix of
$\mathbb{G}$ as

$\mathbb{G}^{-1}=(\mathbb{I}-2\eta^{3}\mathbb{A}^{-1}\mathbb{B}+(\eta^{3})^{2}\mathbb{A}^{-1}\mathbb{C})^{-1}\mathbb{A}^{-1}\,,$
(8)

and expand it up to second order in $\eta^{3}$ as

$\mathbb{G}^{-1}=\mathbb{A}^{-1}+2\eta^{3}\mathbb{A}^{-1}\mathbb{B}\mathbb{A}^{-1}-(\eta^{3})^{2}\mathbb{A}^{-1}\mathbb{C}\mathbb{A}^{-1}+4(\eta^{3})^{2}(\mathbb{A}^{-1}\mathbb{B})^{2}\mathbb{A}^{-1}+O((\eta^{3})^{3})\,.$
(9)

Returning to index notation, we notice that a pre-multiplication of a tensor
by $\mathbb{A}^{-1}$ corresponds to raising the first index of that tensor,
whereas a post-multiplication by $\mathbb{A}^{-1}$ corresponds to raising its
second index. Moreover, since
$\accentset{\circ}{c}_{\alpha\beta}=\accentset{\circ}{b}_{\alpha\gamma}\accentset{\circ}{b}^{\gamma}_{\,\,\beta}$,
we can write, in matrix notation,
$\mathbb{C}=\mathbb{B}\mathbb{A}^{-1}\mathbb{B}$. With these results, we can
rewrite Eq. (9) in index notation as

$G^{\alpha\beta}=\accentset{\circ}{a}^{\alpha\beta}+2\eta^{3}\accentset{\circ}{b}^{\alpha\beta}+3(\eta^{3})^{2}\accentset{\circ}{c}^{\alpha\beta}+O((\eta^{3})^{3})\,,$
(10)

which expresses the contravariant components of the 3D metric as a function of
the thickness coordinate and surface quantities only.

With the above results at hand, we can finally express the deformation
gradient, solely in terms of surface quantities, up to second order in the
thickness coordinate. By substituting Eqs.(1-3) and Eq. (10) into Eq. (4), we
obtain

$\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{a}_{\alpha}\otimes\accentset{\circ}{\mathbf{a}}^{\alpha}+\mathbf{n}\otimes\accentset{\circ}{\mathbf{n}}$
(11)
$\displaystyle+\eta^{3}(-b^{\alpha}_{\,\,\beta}+\accentset{\circ}{b}^{\alpha}_{\,\,\beta})\mathbf{a}_{\alpha}\otimes\accentset{\circ}{\mathbf{a}}^{\beta}$
$\displaystyle+(\eta^{3})^{2}(\accentset{\circ}{c}^{\alpha}_{\,\,\beta}-b_{\eta}^{\,\,\alpha}\accentset{\circ}{b}^{\eta}_{\,\,\beta})\mathbf{a}_{\alpha}\otimes\accentset{\circ}{\mathbf{a}}^{\beta}+O((\eta^{3})^{3})\,,$

where we used the Weingarten equations (do Carmo, 2016),
$\accentset{\circ}{\mathbf{n}},_{\alpha}=-\accentset{\circ}{b}_{\alpha}^{\,\,\beta}\accentset{\circ}{\mathbf{a}}_{\beta}$
and $\mathbf{n},_{\alpha}=-b_{\alpha}^{\,\,\beta}\mathbf{a}_{\beta}$. Note
that this expression of the deformation gradient reduces to that of a plate,
when $\accentset{\circ}{b}_{\alpha}^{\,\,\beta}=0$ (i.e., flat mid-surface).

We now seek to further simplify Eq. (11), in the limit of thin shells. To do
so, we define the principal curvatures in the undeformed and deformed
configuration, respectively, as
$(\accentset{\circ}{k}_{1},\accentset{\circ}{k}_{2})$ and $(k_{1},k_{2})$. In
the limit of thin linear elastic shells, the curvatures in both the undeformed
and deformed configurations have to be small compared to the average thickness
of the shell $\bar{h}$ (Niordson, 1985; Pietraszkiewicz, 1980), that is
$\bar{h}\cdot\max{(\accentset{\circ}{k}_{1},\accentset{\circ}{k}_{2})}\ll 1$
and $\bar{h}\cdot\max{(k_{1},k_{2})}\ll 1$. Therefore, in the limit of thin
linear elastic shells, the linear and quadratic terms of the deformation
gradient, Eq. (11), are negligible, with the deformation gradient that can
finally be expressed as

$\mathbf{F}=\mathbf{a}_{\alpha}\otimes\accentset{\circ}{\mathbf{a}}^{\alpha}+\mathbf{n}\otimes\accentset{\circ}{\mathbf{n}}\,.$
(12)

This reduced expression of $\mathbf{F}$, simplified for thin linear elastic
shells, will be instrumental to perform the dimensional reduction of the
magneto-elastic energy, as will be shown in Section 3.

###  2.3 The reduced elastic energy

Now that the expression of the deformation gradient for thin shells has been
derived, the last task needed before moving to the dimensional reduction
procedure is to recall the reduced elastic energy, as proposed by Koiter
(Koiter, 1945; Niordson, 1985) and used by many authors since then (Budiansky
and Hutchinson, 1972; Paulose and Nelson, 2013; Hutchinson, 2016). As commonly
adapted in the shell literature (Niordson, 1985), we define the following mid-
surface stretching and bending strain tensors

$E_{\alpha\beta}=\frac{1}{2}(a_{\alpha\beta}-\accentset{\circ}{a}_{\alpha\beta})\,$
(13)

$K_{\alpha\beta}=b_{\alpha\beta}-\accentset{\circ}{b}_{\alpha\beta}\,,$ (14)

the elastic energy of a linearly elastic shell can then be expressed as
(Niordson, 1985)

$\displaystyle\mathcal{U}_{\textup{e}}$
$\displaystyle=\int_{\mathcal{S}}\frac{Eh}{2(1-\nu^{2})}[(1-\nu)E_{\alpha}^{\beta}E_{\beta}^{\alpha}+\nu
E_{\alpha}^{\alpha}E_{\beta}^{\beta}]\,\mathrm{d}\accentset{\circ}{\omega}$
(15)
$\displaystyle+\int_{\mathcal{S}}\frac{Eh^{3}}{24(1-\nu^{2})}[(1-\nu)K_{\alpha}^{\beta}K_{\beta}^{\alpha}+\nu
K_{\alpha}^{\alpha}K_{\beta}^{\beta}]\,\mathrm{d}\accentset{\circ}{\omega}\,,$

where
$\mathrm{d}\accentset{\circ}{\omega}=\sqrt{|\mathrm{det}(\accentset{\circ}{a}_{\alpha\beta})|}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}=\accentset{\circ}{a}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}$
is the differential area, $E$ is the Young’s modulus of the material and $\nu$
its Poisson ratio. In the test cases that we will focus on in this work, the
material properties and the thickness will be assumed to be homogeneous. The
elastic energy in Eq. (15) can also be written in a more compact direct
notation, by introducing the trace operator in the undeformed metric “tr ”
such that for example $\mbox{\rm
tr\,}\mathbf{a}=\accentset{\circ}{a}^{\alpha\beta}a_{\alpha\beta}$, so that

$\displaystyle\mathcal{U}_{\textup{e}}$
$\displaystyle=\frac{Eh}{8(1-\nu^{2})}\int_{\mathcal{S}}[(1-\nu)\mbox{\rm
tr\,}(\mathbf{a}-\accentset{\circ}{\mathbf{a}})^{2}+\nu\mbox{\rm
tr\,}^{2}(\mathbf{a}-\accentset{\circ}{\mathbf{a}})]\,\mathrm{d}\accentset{\circ}{\omega}$
(16) $\displaystyle+\frac{D}{2}\int_{\mathcal{S}}[(1-\nu)\mbox{\rm
tr\,}(\mathbf{b}-\accentset{\circ}{\mathbf{b}})^{2}+\nu\mbox{\rm
tr\,}^{2}(\mathbf{b}-\accentset{\circ}{\mathbf{b}})]\,\mathrm{d}\accentset{\circ}{\omega}\,,$

where we have also defined the bending stiffness of the shell as
${D=1/12Eh^{3}/(1-\nu^{2})}$.

In this Section, we have described the geometry of the shell, writing the
deformation gradient for thin linear elastic shells, and recalling some basic
concepts in shell mechanics that will be useful for the remaining of the
paper, where we will introduce the magnetic Helmholtz free energy and study
the coupling between mechanics and magnetism. We are now ready to move on to
the next Section, where we will obtain the reduced energy for magneto-active
shells.

##  3 Reduced magneto-elastic energy via dimensional reduction

The goal of this Section is to reduce the dimension of the magnetic energy to
obtain a reduced order model for non-axisymmetric thin elastic magnetic
shells. This model will provide a much faster alternative to multiphysics 3D
simulations, in addition to being amenable to physical interpretation, as we
will show in Section 4.

The total energy $\mathcal{U}$ of the shell is the sum of its elastic energy
$\mathcal{U}_{\textup{e}}$, the potential energy of external loads such as
pressure, and the magnetic energy $\mathcal{U}_{\textup{m}}$, which we want to
reduce from 3D to 2D. The magnetic energy $\mathcal{U}_{\textup{m}}$ can be
expressed as the Helmholtz free energy for ideal hard-magnetic soft materials
(Kim et al., 2018b; Zhao et al., 2019), and is inherently 3D, meaning that it
is an energy per unit volume. The decomposition of the total energy into
elastic and non elastic terms is also found in other contexts, as in the case
of the swelling of hydrogels with the Flory-Rehner energy (Flory and Rehner,
1943a, b). This decomposition ensures that the elastic part can be linked to
any of the most common constitutive models, including the Kirchhoff Saint-
Venant strain energy used in the Koiter shell model or the neo-Hookean energy
commonly used for elastomeric materials.

For the magnetic Helmholtz free energy, we will consider the model proposed by
Zhao et al. (2019), which is valid for ideal hard-magnetic soft materials. The
basic underlying assumption, following the physical observations in (Bertotti,
1998), is that the magnetic flux density $\mathbf{B}$ of the hard-magnetic
soft material in the reference configuration is linearly related to the
applied magnetic field $\mathbf{H}$. This assumption is generally true for
hard-magnetic soft materials, where the field strength required for actuation
is much lower than the coercivity (Bertotti, 1998).

The magnetic energy is then simplified as the work required to align the
residual magnetic moment of the material along the external magnetic field for
3D scale-free materials, written as (Bertotti, 1998; Kim et al., 2018b; Zhao
et al., 2019)

$\mathcal{U}_{\textup{m}}=-\frac{1}{\mu_{0}}\int_{\mathcal{B}}\mathbf{F}\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}\
\textup{d}V\,,$ (17)

where $\mu_{0}$ is the vacuum permeability, $\mathbf{B}^{\textup{r}}$ is the
vector of residual magnetic flux density, $\mathbf{B}^{\textup{a}}$ is the
vector of externally applied magnetic flux density, and $V$ is the reference
volume of the 3D body $\mathcal{B}$.

In preparation for the dimensional reduction procedure, we now need to
represent the two magnetic fields, namely the residual magnetic flux density
and the applied magnetic flux density, via their Cartesian components, noting
that they can be functions of the spatial coordinates:
$\mathbf{B}^{\textup{r}}={{B}^{\textup{r}}}^{i}\mathbf{e}_{i}$ and
$\mathbf{B}^{\textup{a}}={{B}^{\textup{a}}}^{i}\mathbf{e}_{i}$. We then define
the normalized components as

$\displaystyle{\hat{B}^{\textup{r}}}^{i}$
$\displaystyle=\frac{{B^{\textup{r}}}^{i}}{|\mathbf{B}^{\textup{r}}|}\,,$ (18)
$\displaystyle{\hat{B}^{\textup{a}}}^{i}$
$\displaystyle=\frac{{B^{\textup{a}}}^{i}}{|\mathbf{B}^{\textup{a}}|}\,.$

Then, substituting the expressions for the deformation gradient (Eq. (12)) and
for the volume measure (Eq. (2)) in the 3D magnetic energy (Eq. (17)), and
integrating along the thickness, yields

$\mathcal{U}_{\textup{m}}=hE\lambda_{\textup{m}}\int_{\mathcal{S}}B^{ij}_{0}[(\mathbf{a}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+(\mathbf{n}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})]\,\mathrm{d}\accentset{\circ}{\omega}\,,$
(19)

where, given the thin shell assumption, we neglected $h\accentset{\circ}{H}$
and $h^{2}\accentset{\circ}{K}$ compared to unity. Indeed, the mean curvature
is of the order of the principal curvatures of the shell, while the Gaussian
curvature is equal to their product. Recall that, in Section 2, we highlighted
that the assumptions on the undeformed shape of the shell are
$h\cdot\max{(\accentset{\circ}{k}_{1},\accentset{\circ}{k}_{2})}\ll 1$, and,
hence, $h\accentset{\circ}{H}\ll 1$ and $h^{2}\accentset{\circ}{K}\ll 1$.
Moreover, according to Yan et al. (2020), we use the magneto-elastic
dimensionless parameter

$\lambda_{\textup{m}}=\frac{|\mathbf{B}^{\textup{r}}||\mathbf{B}^{\textup{a}}|}{\mu_{0}E}\,,$
(20)

which represents the ratio between an equivalent magnetic pressure and the
Young’s modulus. If $\lambda_{\textup{m}}\gg 1$, the shell will be
_magnetically compliant_ , whereas if $\lambda_{\textup{m}}\ll 1$, it will be
_magnetically stiff_. Furthermore, we have defined

$B_{0}^{ij}=-\frac{1}{h}\int_{-h/2}^{h/2}{\hat{B}^{\textup{r}}}^{i}{\hat{B}^{\textup{a}}}^{j}\,\textup{d}\eta^{3}\,,$
(21)

which can be seen as the Cartesian components of the
$\mathbf{B}_{0}=B_{0}^{ij}\mathbf{e}_{i}\otimes\mathbf{e}_{j}$ tensor. The
components $B_{0}^{ij}$ are the average of the product
${\bar{B}^{\textup{r}}}^{i}{\bar{B}^{\textup{a}}}^{j}$ along the thickness of
the shell, and quantify the interaction between the residual and the
externally applied magnetic flux densities.

The reduced energy in Eq. (19) is now suitable for the modeling of thin shells
made of hard MREs subject to external magnetic fields, which can be either
homogeneous or linear (constant gradient fields).

In the next Section, we will investigate the reduced energy to gain additional
insight and obtain a physical interpretation of how the external magnetic
field acts to deform the shell.

##  4 Reduced magneto-elastic energy as load potentials

The reduced magnetic energy in Eq. (19) is nonlinear and, therefore, opaque to
an intuitive understanding on how the magnetic field interacts with the
mechanics of the shell. Therefore, we seek to provide a physical
interpretation of the newly reduced magnetic energy, in the case where both
the residual and the applied magnetic fields are homogeneous. Previously, we
presented a similar approach, albeit limited to axisymmetric deformations, in
(Yan et al., 2020). The goal of the current Section is to expand the reduced
magnetic energy in Eq. (19) as a function of the displacement field and to
interpret it as the potential energy of _magnetic loads_.

We start by expanding the reduced magnetic energy up to second order in the
displacement field $\mathbf{u}$. The term that we need to expand is
$[(\mathbf{a}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+(\mathbf{n}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})]$.
We write
$\mathbf{a}_{\alpha}=\accentset{\circ}{\mathbf{a}}_{\alpha}+\delta\mathbf{a}_{\alpha}$
and $\mathbf{n}=\accentset{\circ}{\mathbf{n}}+\delta\mathbf{n}$, where
$\delta\mathbf{a}_{\alpha}$ and $\delta\mathbf{n}$ are the (finite) changes in
the base vectors corresponding to a (finite) displacement $\mathbf{u}$, and
notice that

$(\mathbf{a}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+(\mathbf{n}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})=\delta_{ij}+(\delta\mathbf{a}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+(\delta\mathbf{n}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})\,,$
(22)

where $\delta_{ij}$ is the Kronecker’s delta.

Then, since $\mathbf{r}=\accentset{\circ}{\mathbf{r}}+\mathbf{u}$, we can
derive the expression for $\delta\mathbf{a}_{\alpha}$ as follows

$\displaystyle\delta\mathbf{a}_{\alpha}$
$\displaystyle=(\mathbf{r}-\accentset{\circ}{\mathbf{r}}),_{\alpha}=(u^{\gamma}\accentset{\circ}{\mathbf{a}}_{\gamma}),_{\alpha}+(w\accentset{\circ}{\mathbf{n}}),_{\alpha}$
(23)
$\displaystyle=u^{\gamma},_{\alpha}\accentset{\circ}{\mathbf{a}}_{\gamma}+u^{\gamma}\accentset{\circ}{\mathbf{a}}_{\gamma},_{\alpha}+w,_{\alpha}\accentset{\circ}{\mathbf{n}}-w\accentset{\circ}{b}_{\alpha}^{\,\,\beta}\accentset{\circ}{\mathbf{a}}_{\beta}$
$\displaystyle=(u^{\eta},_{\alpha}+\Gamma^{\eta}_{\,\,\gamma\alpha}u^{\gamma})\accentset{\circ}{\mathbf{a}}_{\eta}+(\accentset{\circ}{b}_{\gamma\alpha}u^{\gamma}+w,_{\alpha})\accentset{\circ}{\mathbf{n}}-w\accentset{\circ}{b}_{\alpha}^{\,\,\beta}\accentset{\circ}{\mathbf{a}}_{\beta}$
$\displaystyle=(\nabla_{\alpha}u^{\eta}-w\accentset{\circ}{b}_{\alpha}^{\,\,\eta})\accentset{\circ}{\mathbf{a}}_{\eta}+(\accentset{\circ}{b}_{\gamma\alpha}u^{\gamma}+w,_{\alpha})\accentset{\circ}{\mathbf{n}}\,,$

where we used the Weingarten equation
$\accentset{\circ}{\mathbf{n}},_{\alpha}=-\accentset{\circ}{b}_{\alpha}^{\,\,\beta}\accentset{\circ}{\mathbf{a}}_{\beta}$,
$\Gamma^{\eta}_{\,\,\gamma\alpha}=\accentset{\circ}{\mathbf{a}}_{\gamma},_{\alpha}\cdot\accentset{\circ}{\mathbf{a}}^{\eta}$
is the Christoffel symbol of the second kind, and $\nabla_{\alpha}$ denotes
the covariant derivative along $\accentset{\circ}{\mathbf{a}}_{\alpha}$ (do
Carmo, 2016). We note that this expression is geometrically exact, given that
the covariant base vectors are linear functions of the displacement field and
its derivatives.

On the contrary, the normal vector is by definition a nonlinear function of
the displacement field. The expansion of $\delta\mathbf{n}$ up to second order
reads (Deserno, 2004)

$\delta\mathbf{n}=-\frac{1}{2}q_{\alpha}q^{\alpha}\accentset{\circ}{\mathbf{n}}-(q^{\alpha}-q^{\gamma}U^{\,\,\alpha}_{\gamma})\accentset{\circ}{\mathbf{a}}_{\alpha}+O(|\mathbf{u}|^{3})\,,$
(24)

where $q_{\alpha}=w,_{\alpha}+\accentset{\circ}{b}_{\gamma\alpha}u^{\gamma}$
are the covariant components of the rotation vector (Niordson, 1985), and
$U^{\,\,\alpha}_{\gamma}=\nabla_{\gamma}u^{\alpha}-w\accentset{\circ}{b}_{\gamma}^{\,\,\alpha}$
are the mixed components of the surface gradient of the displacement field,
that is $\nabla_{\\!\mathcal{S}}\mathbf{u}$.

Finally, substituting Eqs. (23), (24) and (22) into (19), we get

$\displaystyle\mathcal{U}_{\textup{m}}^{2}$
$\displaystyle=hE\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}\delta_{ij}\,\mathrm{d}\accentset{\circ}{\omega}$
(25) $\displaystyle+h$ $\displaystyle
E\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}U_{\alpha}^{\,\,\eta}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}-hE\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}\frac{1}{2}q_{\alpha}q^{\alpha}(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}$
$\displaystyle+h$ $\displaystyle
E\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}q^{\gamma}U^{\alpha}_{\,\,\gamma}(\accentset{\circ}{\mathbf{a}}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}\,,$

where $\mathcal{U}_{\textup{m}}^{2}$ denotes the reduced magnetic energy up to
second order in the displacement field. This expansion represents the first
step towards a mechanical interpretation of the reduced magnetic energy, where
the aim is to understand how the magnetic field loads the shell, in terms of
equivalent forces and torques. In the next two subsections, we will analyze
each term (apart from the constant term in the first line of Eq. (25)), and
aim at manipulating them to understand their physical meaning. For
convenience, we rewrite Eq. (25) as

$\mathcal{U}_{\textup{m}}^{2}=\mathcal{U}_{\textup{m}}^{1}+\mathcal{U}_{\textup{m}}^{q}+\mathcal{U}_{\textup{m}}^{\tau}\,,$
(26)

where we define

$\displaystyle\mathcal{U}_{\textup{m}}^{1}$
$\displaystyle=hE\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}U_{\alpha}^{\,\,\eta}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}\,,$
(27) $\displaystyle\mathcal{U}_{\textup{m}}^{q}$
$\displaystyle=-hE\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}\frac{1}{2}q_{\alpha}q^{\alpha}(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}\,,$
$\displaystyle\mathcal{U}_{\textup{m}}^{\tau}$
$\displaystyle=hE\lambda_{\textup{m}}\int_{\mathcal{S}}B_{0}^{ij}q^{\gamma}U^{\alpha}_{\,\,\gamma}(\accentset{\circ}{\mathbf{a}}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})\,\mathrm{d}\accentset{\circ}{\omega}\,.$

We will first tackle the linear term $\mathcal{U}_{\textup{m}}^{1}$ and then
move on to the two second order terms, namely $\mathcal{U}_{\textup{m}}^{q}$
and $\mathcal{U}_{\textup{m}}^{\tau}$.

###  4.1 Linear term of the reduced magnetic energy,
$\mathcal{U}_{\textup{m}}^{1}$

We expand the linear term, $\mathcal{U}_{\textup{m}}^{1}$, as

$\displaystyle\mathcal{U}_{\textup{m}}^{1}$
$\displaystyle=hE\lambda_{\textup{m}}\int_{\mathcal{S}}\Bigl{[}\nabla_{\alpha}u^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})-w\accentset{\circ}{b}_{\alpha}^{\,\,\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\Bigr{]}\,\mathrm{d}\accentset{\circ}{\omega}\,,$
(28)
$\displaystyle=hE\lambda_{\textup{m}}\int_{\mathcal{S}}\Bigl{[}\nabla_{\alpha}u^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})-\frac{1}{h}p_{\textup{mag}}w\Bigr{]}\,\mathrm{d}\accentset{\circ}{\omega}\,,$

where we defined an equivalent (dimensionless) magnetic pressure

$p_{\textup{mag}}=h\accentset{\circ}{b}_{\alpha}^{\,\,\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})=h\accentset{\circ}{\mathbf{b}}\colon\\!\mathbf{B}_{0}\,,$
(29)

which we can also express as the double inner product between the curvature
tensor
$\accentset{\circ}{\mathbf{b}}=\accentset{\circ}{b}_{\alpha\beta}\accentset{\circ}{\mathbf{a}}^{\alpha}\otimes\accentset{\circ}{\mathbf{a}}^{\beta}$
and the (dimensionless) magnetic tensor
$\mathbf{B}_{0}=B_{0}^{ij}\mathbf{e}_{i}\otimes\mathbf{e}_{j}$. To interpret
the first term of the r.h.s. in Eq. (28), we need to expand the covariant
derivative and integrate by parts, as follows

$\displaystyle\int_{\mathcal{S}}\nabla_{\alpha}u^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\accentset{\circ}{a}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}$
(30)
$\displaystyle=\int_{\mathcal{S}}[u^{\eta},_{\alpha}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+\Gamma_{\gamma\alpha}^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})u^{\gamma}]\accentset{\circ}{a}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}$
$\displaystyle=\int_{\mathcal{S}}(u^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\accentset{\circ}{a}),_{\alpha}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}+\frac{1}{h}\int_{\mathcal{S}}f_{\gamma}u^{\gamma}\accentset{\circ}{a}\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}\,.$

Here, we have defined an equivalent (dimensionless) membrane distributed force
$\mathbf{f}_{\textup{mag}}=~{}f_{\gamma}\accentset{\circ}{\mathbf{a}}^{\gamma}$
as

$f_{\gamma}=h\Gamma_{\gamma\alpha}^{\eta}B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\eta}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})-h\frac{(B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\gamma}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})\accentset{\circ}{a}),_{\alpha}}{\accentset{\circ}{a}}\,,$
(31)

where the last term is in the form of the covariant divergence of a vector
field. When the shells we consider have a fixed boundary or are boundary-free,
the divergence term in Eq. (30) integrates to zero and we can conclude that

$\mathcal{U}_{\textup{m}}^{1}=E\lambda_{\textup{m}}(f_{\gamma}u^{\gamma}-p_{\textup{mag}}w)=E\lambda_{\textup{m}}(\mathbf{f}_{\textup{mag}}\cdot\check{\mathbf{u}}-p_{\textup{mag}}w)\,.$
(32)

This equation suggests that, at the first order in the displacement field, the
reduced magneto-elastic energy can be interpreted as the energy potential of
distributed magnetic in-plane forces $\mathbf{f}_{\textup{mag}}$ and pressure
$p_{\textup{mag}}$. This interpretation of $\mathcal{U}_{\textup{m}}^{1}$
generalizes the one we derived for axisymmetric deformations presented in (Yan
et al., 2020), which also resulted in the definition of magnetic in-plane
forces and pressure.

###  4.2 Nonlinear terms of the reduced magnetic energy,
$\mathcal{U}_{\textup{m}}^{q}$ and $\mathcal{U}_{\textup{m}}^{\tau}$

We have just provided a mechanical interpretation to the linear term in the
reduced magneto-elastic energy, namely $\mathcal{U}_{\textup{m}}^{1}$, but we
still need to do the same for the second-order terms
$\mathcal{U}_{\textup{m}}^{q}$ and $\mathcal{U}_{\textup{m}}^{\tau}$.

We start by looking at $\mathcal{U}_{\textup{m}}^{q}$ which, as shown by Eq.
(27), is quadratic in the linear rotation vector $\mathbf{q}$. Indeed,
$\mathcal{U}_{\textup{m}}^{q}$ corresponds to the energy potential of a
distribution of linear torques throughout the shell, analogous to a
distribution of effective rotational springs of constant stiffness $k$,
proposed by Yan et al. (2020). This stiffness $k$, in our more general case,
can be defined as

$k=B_{0}^{ij}(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})=(\accentset{\circ}{\mathbf{n}}\otimes\accentset{\circ}{\mathbf{n}})\colon\\!\mathbf{B}_{0}\,,$
(33)

which, depending on its sign, can drastically change the nature of the
effective rotational springs. As an illustrative example, let us consider a
case where $\mathbf{B}^{\textup{r}}={{B}^{\textup{r}}}\mathbf{e}_{3}$ and
$\mathbf{B}^{\textup{a}}={{B}^{\textup{a}}}\mathbf{e}_{3}$, considered also in
(Yan et al., 2020). This specific case implies that
$k=-\hat{B}^{\textup{r}}\hat{B}^{\textup{a}}(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{3})^{2}$
or, stated differently, that the sign of the stiffness of the rotational
springs is opposite to that of
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}$. This simple example
conveys the importance of such a physical interpretation of the reduced
magneto-elastic energy, since a quick check of the sign of
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}$ provides a first
qualitative assessment of the case or problem at hand.

The last second order term in Eq. (26), namely
$\mathcal{U}_{\textup{m}}^{\tau}$, can be interpreted as the energy potential
of distributed torques, which depend on the surface gradient of the
displacement field. We define this dimensionless distributed torque as
$\boldsymbol{\tau}(\mathbf{u})~{}=\tau_{\gamma}\accentset{\circ}{\mathbf{a}}^{\gamma}$,
where

$\tau_{\gamma}=B_{0}^{ij}(\accentset{\circ}{\mathbf{a}}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})U^{\alpha}_{\,\,\gamma}\,$
(34)

is a linear function of the displacement field, since
$U^{\alpha}_{\,\,\gamma}$ is linear with the displacement field.

Using the above results on the mechanical interpretation, namely Eqs. (32),
(33) and (34), the reduced magnetic energy up to second order in the
displacement field reads

$\mathcal{U}_{\textup{m}}^{2}=E\lambda_{\textup{m}}\int(\mathbf{f}_{\textup{mag}}\cdot\check{\mathbf{u}}-p_{\textup{mag}}w)\,\mathrm{d}\accentset{\circ}{\omega}-\frac{1}{2}hE\lambda_{\textup{m}}\int
k\mathbf{q}\cdot\mathbf{q}\,\mathrm{d}\accentset{\circ}{\omega}+hE\lambda_{\textup{m}}\int\boldsymbol{\tau}(\mathbf{u})\cdot\mathbf{q}\,\mathrm{d}\accentset{\circ}{\omega}\,,$
(35)

for homogeneous magnetic fields and boundary-free (or clamped) shells.
Moreover, we recall that the above expression (Eq. (35)) of the reduced
magnetic energy is geometrically exact up to second order in the displacement
field, and that we assumed homogeneous Young’s modulus and thickness
throughout. In summary, we have shown that the reduced magneto-elastic energy
is equivalent to an energy potential of _magnetic loads_ , combining in-plane
forces, linear pressure, and torques, distributed across the mid-surface of
the shell.

In the next Section, we will contrast this second-order energy against the
fully nonlinear version, that is Eq. (19), to understand whether
$\mathcal{U}_{\textup{m}}^{2}$ represents a good approximation and whether
there are terms in the expansion that predominate over others.

##  5 Numerical implementation

###  5.1 Minimization of the reduced 2D energy

The reduced magneto-elastic energy in Eq. (19) is fully nonlinear and, as with
the majority of problems involving thin shells, its minimization (together
with the elastic energy) has to be tackled via numerical methods. The
mechanical interpretation provided in Section 4 can serve as a tool to
understand the basics of a specific problem and can guide its solution, but
quantitative solutions can only be obtained via the numerical minimization of
the reduced magneto-elastic energy (Eq. (19)). In particular cases, as we will
see in Section 6, the numerical minimization of $\mathcal{U}_{\textup{m}}^{2}$
might suffice to obtain a satisfactory solution of a problem, resulting in an
even faster numerical minimization scheme given that
$\mathcal{U}_{\textup{m}}^{2}$ is only quadratic in $\mathbf{u}$.

The total energy of the shell $\mathcal{U}$ is the sum of the elastic energy
$\mathcal{U}_{\textup{e}}$, the potential of the live pressure
$p\Delta\mathcal{V}$ (with $\Delta\mathcal{V}$ being the change in volume of
the shell), and the reduced magnetic energy $\mathcal{U}_{\textup{m}}$. For
our numerical scheme, we nondimensionalize the total energy by
$EhR^{2}/(8(1-\nu^{2}))$ to obtain

$\overline{\mathcal{U}}=\overline{\mathcal{U}}_{\textup{s}}+\frac{1}{3}\left(\frac{h}{R}\right)^{2}\overline{\mathcal{U}}_{\textup{b}}+\frac{p}{E}\Delta\overline{\mathcal{V}}+\lambda_{\textup{m}}\overline{\mathcal{U}}_{\textup{m}}\,,$
(36)

where $\overline{\mathcal{U}}_{\textup{s}}$,
$\overline{\mathcal{U}}_{\textup{b}}$, and $\Delta\overline{\mathcal{V}}$ are
the dimensionless stretching energy, bending energy, and change in volume. We
can write these terms as (Pezzulla and Reis, 2019)

$\displaystyle\overline{\mathcal{U}}_{\textup{s}}$
$\displaystyle=\frac{1}{R^{2}}\int[(1-\nu)\mbox{\rm
tr\,}(\mathbf{a}-\accentset{\circ}{\mathbf{a}})^{2}+\nu\mbox{\rm
tr\,}^{2}(\mathbf{a}-\accentset{\circ}{\mathbf{a}})]\,\mathrm{d}\accentset{\circ}{\omega}\,,$
(37) $\displaystyle\overline{\mathcal{U}}_{\textup{b}}$
$\displaystyle=\int[(1-\nu)\mbox{\rm
tr\,}(\mathbf{b}-\accentset{\circ}{\mathbf{b}})^{2}+\nu\mbox{\rm
tr\,}^{2}(\mathbf{b}-\accentset{\circ}{\mathbf{b}})]\,\mathrm{d}\accentset{\circ}{\omega}\,,$
$\displaystyle\Delta\overline{\mathcal{V}}$
$\displaystyle=\frac{8(1-\nu^{2})}{3hR^{2}}\biggl{[}\int\mathbf{r}\cdot\mathbf{n}\,\mathrm{d}\omega-\int\accentset{\circ}{\mathbf{r}}\cdot\accentset{\circ}{\mathbf{n}}\,\mathrm{d}\accentset{\circ}{\omega}\biggr{]}\,,$
$\displaystyle\overline{\mathcal{U}}_{\textup{m}}$
$\displaystyle=\frac{8(1-\nu^{2})}{R^{2}}\int
B^{ij}_{0}[(\mathbf{a}_{\alpha}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{a}}^{\alpha}\cdot\mathbf{e}_{i})+(\mathbf{n}\cdot\mathbf{e}_{j})(\accentset{\circ}{\mathbf{n}}\cdot\mathbf{e}_{i})]\,\mathrm{d}\accentset{\circ}{\omega}\,,$

where we added the (dimensionless) reduced magneto-elastic energy
$\overline{\mathcal{U}}_{\textup{m}}$ derived in Section 3. The energy written
in Eq. (36) is minimized numerically using the commercial software COMSOL
Multiphysics, similarly to how we described in (Pezzulla and Reis, 2019), even
though the procedure presented there was applied to a purely elastic energy,
without any magnetic contribution, and to a $1$D functional since the model
was axisymmetric. To discretize the domain of parametrization of the
undeformed mid-surface, we have made use of triangular elements using Argyris
shape functions (Argyris et al., 1968), a fifth-order polynomial approximation
that allows the derivatives of the state variables to be assigned at the
boundaries.

In figure 2(a), we present a representative example of a discretized domain of
parametrization $(\eta^{1},\eta^{2})$ for a hemispherical shell clamped at its
equator. In this example, the curvilinear coordinates $({\eta^{1},\eta^{2}})$
are spherical coordinates with $\eta^{1}\in[0,\pi/2]$ for the colatitude and
$\eta^{2}\in[0,2\pi)$ for the longitude. The color map indicates the
dimensionless normal displacement, $w/h$, induced by a poking force normal to
the surface at a point along the parallel at $45^{\circ}$. The corresponding
discretized deformed mid-surface in Euclidean space is shown in figure 2(b).

Figure 2: Domain of parametrization for a spherical shell using spherical
coordinates $\eta^{1}\in[0,\pi]$ and $\eta^{2}\in[0,2\pi)$ (a). The color
indicates the dimensionless normal displacement $w/h$ as a consequence of a
poking normal force at a point along the parallel at $45^{\circ}$. (b)
Deformed mid-surface of the shell in space.

In all the examples that we will present in Section 6, we will consider
hemispherical shells, clamped along the equator, in some cases with defects
precisely engineered on their surface (Lee et al., 2016). This spherical
geometry corresponds to the special case with $u^{\alpha}=0$, $w=0$, and the
clamped boundary translates to $w,_{1}=0$ at the equator ($\eta^{1}=\pi/2$)
(Niordson, 1985). Moreover, continuity is enforced between the boundary
segments at $\eta^{2}=0$ and $\eta^{2}=2\pi$.

###  5.2 FEM simulations of the 3D model

Our reduced magnetic shell model presented in Secs. 2 and 3 and its numerical
implementation presented above in Sec. 5.1 will be validated in Sec. 6 using
3D finite element modeling (3D FEM). The 3D FEM was performed in the
commercial software package Abaqus/Standard using the user-defined 8-node
brick element, proposed by Zhao et al. (2019) for modeling hard-magnetic
deformable solids under a uniform magnetic field. This FEM framework and has
been thoroughly validated by experiments (Zhao et al., 2019). The code was
also extended by some of the authors to include the case where the external
magnetic field is a constant-gradient field, and successfully applied to study
magnetic beams (Yan et al., 2021).

The user-defined element that we will use was developed based on the 3D
continuum theory of hard-magnetic soft materials (Zhao et al., 2019). The
elastic behavior of the material is described by the neo-Hookean model. The
magnetic interaction between the material and externally applied fields is
considered by introducing a magnetic part of the Cauchy stress, derived from
the 3D magnetic energy potential in Eq. (17).

In our FEM simulations, a 3D deformable hemispherical shell containing a
geometric defect was discretized by the user elements. The material of the
shell was assumed to be incompressible, with a bulk modulus 100 times larger
than its shear modulus. The shell was fixed at the equator and subjected to
combined magnetic and mechanical loading. Contact and distributed mechanical
loads (e.g., pressure) were imposed on the shell through a dummy mesh of C3D8R
solid elements, which shared the same nodes with the user elements. The dummy
material had a negligible elastic modulus ($10^{-20}\,$Pa) compared to the MRE
of the shell ($1.8\,$MPa). Geometric nonlinearities were taken into account
throughout the simulations. The geometric, material, and loading parameters,
provided in Sec. 6, are identical to those used in the shell model. We recall
that only problems where the external magnetic field is homogeneous will be
considered in Section 6, even though both the 3D FEM and the reduced magneto-
elastic energy can be used to study problems with constant-gradient magnetic
fields. We hope that future work will leverage the intrinsic ability of our
theory to consider cases of constant gradient fields.

##  6 Validation of the reduced magnetic shell model

We proceed by validating the reduced magnetic shell model against the 3D FEM
model for magnetorheological elastomers implemented in Abaqus via the study of
three different test cases:

  * (i)

Non-axisymmetric point indentation under magnetic loading (subsection 6.1),
which will test the model when the magnetic loading is axisymmetric while the
mechanical one is not;

  * (ii)

Pressure buckling under asymmetric magnetic loading where the residual
magnetization vector and the external magnetic field are in the same plane
(subsection 6.2), which will test the ability of the model to describe a sub-
critical instability under asymmetric magnetic conditions;

  * (iii)

Pressure buckling under asymmetric magnetic loading where the residual
magnetization vector and the external magnetic field are not in the same plane
(subsection 6.3), which will further test the model with the study of the
elastic instability when the residual and external magnetic fields are non-
coplanar.

Throughout the validation procedure, and for all test cases, (i), (ii), and
(iii), we will also investigate the roles of the reduced energy terms, namely
the linear term $\mathcal{U}_{\textup{m}}^{1}$, the second-order torque term
(quadratic in the rotation) $\mathcal{U}_{\textup{m}}^{q}$, and the full
second-order approximation $\mathcal{U}_{\textup{m}}^{2}$. We fix $R=25.2$ mm,
$h=0.32$ mm, $E=1.8$ MPa, and $\nu=0.5$, corresponding to the dimensionless
numbers $R/h=90$ and $\lambda_{\textup{m}}=0.00184$. Within all the three test
cases, we consider two loading sets identified by the sign of the scalar
product between the residual and the applied (homogeneous) magnetic fields,
namely $\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}>0$ and
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}<0$.

###  6.1 Point indentation under a magnetic field

The first validation test case considers the asymmetric point indentation of a
spherical shell under a uniform and vertical magnetic field, with
$\mathbf{B}^{\textup{a}}=B^{\textup{a}}\mathbf{e}_{3}$ and
$\mathbf{B}^{\textup{r}}=B^{\textup{r}}\mathbf{e}_{3}$. The indentation force
$F$ is set normal to the undeformed mid-surface of the shell and located at
$\eta^{1}=\pi/4$ and $\eta^{2}=0$. To include this external force into the
reduced shell model, we add its (dimensionless) potential to the total energy
in Eq. (36)

$\overline{\mathcal{P}}_{\textup{f}}=\frac{8(1-\nu^{2})}{EhR^{2}}\int_{\mathcal{S}}F\delta(\eta^{1}-\pi/4,\eta^{2})w\,\mathrm{d}\eta^{1}\mathrm{d}\eta^{2}\,,$
(38)

where the force is assumed to be positive if indenting into the shell and
$\delta(\eta^{1},\eta^{2})$ is the Dirac delta function.

In the 3D FEM simulations, we considered one half of the hemispherical shell
(i.e., one quarter of a full spherical shell) with symmetric conditions
imposed on the plane spanned by $\mathbf{e}_{1}$ and $\mathbf{e}_{3}$, in
which the indentation was exerted. The indenter was modeled as a rigid sphere
of radius $0.2\,$mm and Young’s modulus $1.8\,$GPa. The contact between the
indenter and the shell was assumed rigid with neither friction nor
penetration. We discretized the shell using a swept mesh with 10, 125, 100
seeds, respectively, in the thickness, equatorial, and meridional directions.
The mesh in the vicinity of the indentation point was further refined in order
to accurately describe the ensuing localized deformation, the characteristic
length of which scales as $\sqrt{Rh}$ (Abbasi et al., 2021). We have conducted
a convergence study to ensure that the results were independent of this
discretization. In each simulation, a uniform magnetic field at a given level
of flux density (at a fixed value of $\lambda_{\textup{m}}$) was first applied
on the shell, and then, under this fixed field, the shell was indented in a
second step.

Figure 3 shows the results of the simulations in terms of the dimensionless
force $FR/(2\pi D)$, versus the dimensionless normal displacement at the
indentation point, $w/h$. The agreement between our reduced model (solid
lines) and the 3D FEM model (symbols) is excellent. Moreover, the results also
show how the second-order energy $\mathcal{U}_{\textup{m}}^{2}$ reproduces the
results very well, as long as the displacement is not too large (lower than
approximately ten times the thickness). By contrast, the first-order energy
term $\mathcal{U}_{\textup{m}}^{1}$ and second-order torque energy term
$\mathcal{U}_{\textup{m}}^{q}$ fail at replicating the nonlinear behavior,
even for small displacements where, consequently, at least a second-order
reduced magnetic energy is required. Moreover, from the results in Fig. 3, we
notice that the magnetic field modifies the indentation response of the shell
similarly to pressure (Vella et al., 2012; Lazarus et al., 2012; Marthelot et
al., 2017; Hutchinson and Thompson, 2017), since the shell can be either
strengthened or weakened depending on whether the pressure is acting to
inflate or deflate the shell.

Figure 3: Asymmetric indentation of a spherical shell under a magnetic field
in terms of the dimensionless force $FR/(2\pi D)$ versus the dimensionless
normal displacement at the point of indentation $w/h$. The numerical results
from the reduced magnetic shell model (solid lines) are in excellent agreement
with 3D FEM simulations (symbols). Results in red refer to
$\lambda_{\textup{m}}=0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}>0$, whereas those in blue
refer to $\lambda_{\textup{m}}=-0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}<0$. The loading curve for
$\mathbf{B}^{\textup{a}}=\mathbf{0}$ is plotted in black. The results from the
second-order reduced energy $\mathcal{U}_{\textup{m}}^{2}$, the first-order
energy term $\mathcal{U}_{\textup{m}}^{1}$, and the torque energy term
$\mathcal{U}_{\textup{m}}^{q}$ are plotted as dotted lines, dashed lines, and
dotted-dashed lines, respectively.

###  6.2 Pressure buckling under asymmetric magnetic loading
with co-planar $\mathbf{B}^{\textup{r}}$ and $\mathbf{B}^{\textup{a}}$

The previous test case served to validate our reduced model in the case where
the mechanical loading is asymmetric, while the magnetic loading is symmetric.
Next, we want to further test our reduced model in the case where the magnetic
loading is asymmetric and the shell undergoes an elastic instability, namely
pressure-induced buckling. In this test, the residual magnetization vector
$\mathbf{B}^{\textup{r}}$ and the external applied magnetic field
$\mathbf{B}^{\textup{a}}$ are co-planar. If we denote by $\phi^{\textup{r}}$
and $\phi^{\textup{a}}$, the angles $\mathbf{B}^{\textup{r}}$ and
$\mathbf{B}^{\textup{a}}$, respectively, make with $\mathbf{e}_{3}$, we can
distinguish two sub-cases. In the first one, $\mathbf{B}^{\textup{r}}$ and
$\mathbf{B}^{\textup{a}}$ are parallel, that is
$\phi^{\textup{r}}=\phi^{\textup{a}}$ (Figure 4 (a)). In the second one,
$\mathbf{B}^{\textup{r}}$ will stay vertical, that is $\phi^{\textup{r}}=0$,
while $\mathbf{B}^{\textup{a}}$ will vary via its angle $\phi^{\textup{a}}$
(Figure 4 (b)).

Figure 4: Pressure buckling of a spherical shell under asymmetric magnetic
loading for $\lambda_{\textup{m}}=0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}>0$ (in red) and
$\lambda_{\textup{m}}=-0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}<0$ (in blue), with
$\delta/h=1.64$ and $\bar{\phi}_{o}=3.2$. (a) Change of the knockdown factor
$\Delta k_{\textup{d}}$ versus the angle of the shell magnetization vector
$\phi^{\textup{r}}$ and of the external applied field $\phi^{\textup{a}}$,
which are consistently aligned. (b) Change of the knockdown factor $\Delta
k_{\textup{d}}$ versus the angle of the external applied field
$\phi^{\textup{a}}$, with $\phi^{\textup{r}}=0$. Results from the reduced
magnetic energy are plotted as solid lines, whereas triangles represent
results from the three-dimensional model. Results from the second-order
reduced energy $\mathcal{U}_{\textup{m}}^{2}$, the first-order energy term
$\mathcal{U}_{\textup{m}}^{1}$, and the torque energy term
$\mathcal{U}_{\textup{m}}^{q}$ are plotted as circles, asterisks, and
diamonds, respectively.

Given that shells are highly sensitive to imperfections, their measured
critical buckling pressure, even in the absence of any external magnetic
fields, can be much lower than classical predictions for a perfect spherical
shell (Hutchinson et al., 1971; Lee et al., 2016). Consequently, it is common
to introduce the so-called knockdown factor, $k_{\textup{d}}$, defined as the
ratio between the buckling pressure of the shell and the theoretical buckling
pressure of the equivalent perfect shell,
$k_{\textup{d}}=p^{\textup{measured}}/p_{\textup{c}}$. Here, we consider a
spherical shell, clamped along the equator, with a precisely engineered
dimple-like defect at the pole. This defect can be modeled as a deviation from
the spherical mid-surface represented, in terms of the normal displacement, as
$w=-\delta(1-(\eta^{1}/\phi_{o})^{2})^{2}$. For the present test case and as a
representative example, we set the dimensionless defect amplitude to
$\delta/h=1.64$ and the dimensionless defect width to $\bar{\phi}_{o}=3.2$,
with $\bar{\phi}_{o}=\phi_{o}(\sqrt{12(1-\nu^{2})}R/h)^{1/2}$. We will present
the results in terms of the change of the knockdown factor $\Delta
k_{\textup{d}}$, which is defined as the difference between the knockdown
factor obtained when the external magnetic field is on, and that of the shell
without any magnetic fields.

This problem has also been studied in (Yan et al., 2020), where the Abaqus
user-defined solid element developed by Zhao et al. (2019) was employed to
predict the critical buckling load of hard-magnetic shells and tested
successfully against experiments. Next, we first verify our reduced shell
model by comparing with the 3D FEM results presented in (Yan et al., 2020).
Subsequently, using the proposed 2D model, we investigate the relevance of
magnetic energy terms at different orders, Eq. (25), in shell buckling.

Figure 4 (a) shows the results of the change in knockdown factor as a function
of the angle $\phi^{\textup{r}}$, for the case where
$\phi^{\textup{r}}=\phi^{\textup{a}}$. The reduced shell model (solid lines)
is in excellent agreement with the 3D FEM model (triangles). Moreover, the
second-order energy term $\mathcal{U}_{\textup{m}}^{2}$ (circles) replicates
the results obtained with the full energy almost exactly, being therefore a
good approximation of the reduced magnetic energy that can be used for the
study of this class of buckling problems. By contrast, neither the first-order
energy $\mathcal{U}_{\textup{m}}^{1}$ nor the (second-order) torque energy
term $\mathcal{U}_{\textup{m}}^{q}$ can accurately describe the changes in the
knockdown factor. We notice that the first-order energy is adequate to
describe the changes in the knockdown factor for
$\phi^{\textup{r}}\simeq\pi/2$, while the torque energy term can be used for
$\phi^{\textup{r}}\simeq 0$.

In figure 4 (b), we present results for the case when $\phi^{\textup{r}}=0$,
in terms of the change of the knockdown factor as a function of the angle
$\phi^{\textup{a}}$. The reduced shell model (solid lines) accurately
replicates the results from the 3D FEM model (triangles). Also in this case,
the results obtained with the second-order energy term (circles) are in very
good agreement with the full energy model. Moreover, we find that the (second-
order) torque energy term alone, is able to accurately describe the changes in
knockdown factor, while the first-order energy term gives an almost null
contribution.

###  6.3 Pressure buckling under asymmetric magnetic loading with non co-
planar $\mathbf{B}^{\textup{r}}$ and $\mathbf{B}^{\textup{a}}$

Until now, we have validated our reduced shell model in the cases where the
magnetic and mechanical loadings are either symmetric or asymmetric, with the
residual and external magnetic fields being co-planar. In this last test case,
we will evaluate our reduced shell model in the case where the shell undergoes
a pressure buckling instability, with $\mathbf{B}^{\textup{r}}$ and
$\mathbf{B}^{\textup{a}}$ being non co-planar. In particular, as shown in the
inset of Fig. 5, both the residual magnetization vector and the external
applied field, $\mathbf{B}^{\textup{r}}$ and $\mathbf{B}^{\textup{r}}$, are at
angle $\theta=\pi/4$ with respect to the equatorial plane, and are separated
by an angle $\phi$ along the longitudinal direction. Similarly to the previous
test cases, the shell contains a defect at the north pole with $\delta/h=1.64$
and $\bar{\phi}_{o}=3.2$.

Figure 5: Pressure buckling of a spherical shell under asymmetric magnetic
loading, for the more general case when $\mathbf{B}^{\textup{r}}$ and
$\mathbf{B}^{\textup{a}}$ are not in the same plane. $\mathbf{B}^{\textup{r}}$
and $\mathbf{B}^{\textup{a}}$ are at an angle $\theta=\pi/4$ with respect to
the equatorial plane, and separated by an angle $\phi$ along the longitudinal
direction. The change of the knockdown factor $\Delta k_{\textup{d}}$ versus
the angle $\phi$ is shown for $\lambda_{\textup{m}}=0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}>0$ (in red) and
$\lambda_{\textup{m}}=-0.00184$,
$\mathbf{B}^{\textup{r}}\cdot\mathbf{B}^{\textup{a}}<0$ (in blue). The
spherical shells contain a defect at the north pole with dimensionless
amplitude $\delta/h=1.64$ and dimensionless width $\bar{\phi}_{o}=3.2$.
Results from the second-order reduced energy $\mathcal{U}_{\textup{m}}^{2}$
(circles), the first-order energy term $\mathcal{U}_{\textup{m}}^{1}$
(asterisks), and the torque energy term $\mathcal{U}_{\textup{m}}^{q}$
(diamonds) are also shown.

In the 3D FEM of this test case, we modeled the full hemispherical shell with
no assumptions of symmetry. The Abaqus simulations were conducted by following
the procedure proposed in (Yan et al., 2020) for the pressure buckling of
hard-magnetic shells under a uniform field. A swept mesh was created by 8,
400, and 200 seeds in the thickness, equatorial, and meridional directions,
respectively, determined through a convergence study. In particular, the seeds
in the meridional direction were non-uniform with a bias ratio of 8, such that
the mesh was denser near the defect at the north pole, where buckling
occurred.

In Fig. 5, we plot the change of the knockdown factor versus the angle $\phi$.
We find that the results from the reduced shell model are in excellent
agreement with the 3D FEM model. Moreover, the second-order energy term
(circles) provides results that agree with the full energy, as in the previous
test case discussed in Sec. 6.2. Here, the (second-order) torque energy term
is not able to describe the changes in the knockdown factor, while the first-
order energy term is only successful for $\phi\simeq 0$.

##  7 Conclusion and discussion

We have developed a reduced model for magneto-active thin shells with
geometrically exact strain measures. Our model is based on a dimensional
reduction procedure of the Helmholtz free energy for ideal hard-magnetic soft
materials (Zhao et al., 2019). This dimensional reduction relies on the
slenderness of the shells ($h/R\ll 1$) and on the Kirchhoff-Love assumption on
the kinematics, according to which fibers normal to the mid-surface are
assumed not to stretch nor shrink and to remain normal to the mid-surface,
upon deformation. As a result, the reduced energy for thin shells comprises an
elastic energy, split into stretching and bending contributions, the energy
potential of the live pressure, and the reduced magnetic energy, which was
derived in the present work. This reduced energy is linear in the thickness of
the shell, but nonlinear in the displacement field. In its dimensionless form,
the reduced magnetic energy introduces a dimensionless parameter
$\lambda_{\textup{m}}$, which we refer to as the magneto-elastic parameter,
that can be interpreted as the ratio between the magnetic and the elastic
forces in the system.

Since the reduced magnetic energy is highly nonlinear, it is difficult to gain
physical insight without additional simplifications. Therefore, we proceeded
to expand the energy up to second order in the displacement field
$\mathbf{u}$, for the case where both the residual magnetization field and the
external applied field are homogeneous. This expansion yields a first-order
energy $\mathcal{U}_{\textup{m}}^{1}$, which is related to the work done by
equivalent magnetic forces and pressure, and two second-order torque energy
terms, with the most predominant term being proportional to the square of the
rotation vector, that is $\mathcal{U}_{\textup{m}}^{q}$.

We validated our reduced model against a 3D FEM model for MREs and studied
three different test cases, involving the indentation and the pressure
buckling of spherical shells under magnetic loading. In all the test cases,
the results from our reduced model were in excellent agreement with the 3D
FEM. Non-trivially, the second-order simplified model was able to replicate
the 3D FEM results for the buckling cases and for the indentation case, as
long as, in this last case, the normal displacement did not exceed around ten
times the thickness. The second-order energy is therefore a good approximation
of the reduced magnetic energy for the study of the onsets of the
instabilities.

Our model correctly predicts the nonlinear response of thin magnetic shells,
when subject to a combination of magnetic, force, and pressure loadings. In
all of the three test problems, the reduced model provided solutions that were
in excellent agreement with the 3D model, while being an order of magnitude
faster than its 3D counterpart. Indeed, we analyzed the computational
efficiency of the reduced magnetic shell model with respect to the three-
dimensional one, by evaluating the computational times for simulating the
indentation problem (figure 3) with $20$ points in the load-displacement
curves. For both the 3D and 2D models, we measured the computational times
when the minimum meshes that ensured convergence were employed. On a Dell
Precision 7820 workstation with a 12-core CPU (Intel Xeon Gold 6136) and
$192\,$GB RAM, the computational time for the 3D FEM model was $31$ min,
whereas the one for the 2D model was $3.5$ min. This tenfold increase in
efficiency could be useful for studies where an exploration of a vast
parametric space is desired to improve the understanding and design of
magnetic shells. Moreover, as we tested the model also in the presence of
elastic instabilities, we can conclude that the model is suitable to study
those problems where the (subcritical) instabilities, common in thin shells,
play a major role. Not only does this reduced model offer a more effective
alternative to 3D models, but it also provides a physical insight on the
coupling between the magnetic field and the nonlinear mechanics of thin
shells. For example, our analysis, applied to the pressure buckling case
studied by Yan et al. (2020), showed how effective rotational springs play a
major role in the phenomenon, and revealed the existence of a single governing
dimensionless parameter, $\Lambda_{\textup{m}}=\lambda_{\textup{m}}R/h$, which
summarizes all the geometrical and material properties of the system.

However, there are two main points that are worth of further consideration in
future research efforts. The first point regards the self long-range
interactions, which have been neglected in both our reduced model and in the
3D FEM model. While these interactions are not significant in the test cases
that we investigated for validation purposes, there might be applications were
they play a more prominent role. Consider for example a complete spherical
shell that is largely deforming under a magnetic field. In this situation,
parts of the shell in the northern hemisphere might become close to some parts
in the southern hemisphere, making long self-range interactions non-
negligible, similarly to magnetic helical rods (Sano et al., 2021) in near
self-contact. The second point regards a deeper understanding of the
simplified first-order and second-order energy terms. Indeed, while our
physical interpretation as forces, pressure, and torques simplifies the
understanding of the energy, it is not clear why some terms alone are able to
describe the deformation of the shell. For example, in the pressure buckling
represented in Fig. 4 (b), the second-order torque term
$\mathcal{U}_{\textup{m}}^{q}$ is sufficient to describe the entire
phenomenon, with no clear mechanical explanations.

Our proposed model augments the set of emerging reduced models for slender
structures, for planar beams (Yan et al., 2021; Wang et al., 2020; Ciambella
et al., 2018; Ciambella and Tomassetti, 2020) and Kirchhoff-like rods (Sano et
al., 2021). Moreover, our reduced model extends the shell model presented in
(Yan et al., 2020), which was restricted to axisymmetric deformations, now to
non-axisymmetric configurations. Such reduced-order models should act as
valuable tools for the predictive design of magnetic devices that are becoming
increasingly popular in soft robotics (Diller et al., 2014; Huang et al.,
2016; Pece et al., 2017; Seffen and Vidoli, 2016; Loukaides et al., 2014).
Indeed, these models are faster to simulate and more amenable to physical
interpretation. We believe that the derivation of a magnetic plate model will
be a natural next step towards the further and definitive enrichment of the
family of reduced magnetic models. We envision that the different reduced
magnetic models could eventually be combined to study complex magnetic
structures, comprising rods, plates, and shells elements.

## Acknowledgments

We are grateful to John W. Hutchinson for fruitful discussions. We thank
Alessandro Lucantonio for useful exchanges on the numerical implementation of
the model.

## References

  * Abbasi et al. (2021) Abbasi, A., Yan, D., Reis, P.M., 2021.  Probing the buckling of pressurized spherical shells.  J. Mech. Phys. Solids 155, 104545\.
  * Alapan et al. (2020) Alapan, Y., Karacakol, A.C., Guzelhan, S.N., Isik, I., Sitti, M., 2020.  Reprogrammable shape morphing of magnetic soft machines.  Sci. Adv. 6, eabc6414.
  * Argyris et al. (1968) Argyris, J.H., Fried, I., Scharpf, D.W., 1968.  The tuba family of plate elements for the matrix displacement method.  Aeronaut. J. 72, 701–709.
  * Babcock (1983) Babcock, C.D., 1983.  Shell stability.  J. Appl. Mech. 50, 935–940.
  * Bertotti (1998) Bertotti, G., 1998.  Hysteresis in Magnetism: for Physicists, Materials Scientists, and Engineers.  Academic press.
  * Budiansky and Hutchinson (1972) Budiansky, B., Hutchinson, J.W., 1972\.  Buckling of circular cylindrical shells under axial compression, in: Contributions to the Theory of Aircraft Structures. Delft University Press, The Netherlands, pp. 239–259.
  * do Carmo (2016) do Carmo, M., 2016.  Differential Geometry of Curves and Surfaces: Revised and Updated Second Edition.  Dover Books on Mathematics, Dover Publications.
  * Ciambella et al. (2018) Ciambella, J., Favata, A., Tomassetti, G., 2018.  A nonlinear theory for fibre-reinforced magneto-elastic rods.  P. Roy. Soc. Lond. A Mat. 474, 20170703\.
  * Ciambella and Tomassetti (2020) Ciambella, J., Tomassetti, G., 2020\.  A form-finding strategy for magneto-elastic actuators.  Int. J. Nonlin. Mech. 119, 103297\.
  * Cēbers (2003) Cēbers, A., 2003.  Dynamics of a chain of magnetic particles connected with elastic linkers.  J. Phys.: Condens. Matter 15, S1335–S1344.
  * Cēbers and Cīrulis (2007) Cēbers, A., Cīrulis, T., 2007\.  Magnetic elastica.  Phys. Rev. E 76, 031504\.
  * Cēbers and Javaitis (2004) Cēbers, A., Javaitis, I., 2004\.  Bending of flexible magnetic rods.  Phys. Rev. E 70, 021404\.
  * Danas et al. (2012) Danas, K., Kankanala, S.V., Triantafyllidis, N., 2012.  Experiments and modeling of iron-particle-filled magnetorheological elastomers.  J. Mech. Phys. Solids 60, 120–138.
  * Deserno (2004) Deserno, M., 2004.  Notes on Differential Geometry.
  * Diller et al. (2014) Diller, E., Zhuang, J., Lum, G.Z., Edwards, M.R., Sitti, M., 2014.  Continuously distributed magnetization profile for millimeter-scale elastomeric undulatory swimming.  Appl. Phys. Lett. 104, 174101\.
  * Donnell (1977) Donnell, L., 1977.  Beams, plates and shells.  New York: McGraw Hill.
  * Dorfmann and Ogden (2003) Dorfmann, A., Ogden, R.W., 2003\.  Magnetoelastic modelling of elastomers.  Eur. J. Mech. A Solids 22, 497–507.
  * Dreyfus et al. (2005) Dreyfus, R., Baudry, J., Roper, M.L., Fermigier, M., Stone, H.A., Bibette, J., 2005\.  Microscopic artificial swimmers.  Nature 437, 862–865.
  * Flory and Rehner (1943a) Flory, P.J., Rehner, J.J., 1943a.  Statistical mechanics of cross-linked polymer networks i. rubberlike elasticity.  J. Chem. Phys. 11, 512–520.
  * Flory and Rehner (1943b) Flory, P.J., Rehner, J.J., 1943b.  Statistical mechanics of cross-linked polymer networks ii. swelling.  J. Chem. Phys. 11, 521–526.
  * Gerbal et al. (2015) Gerbal, F., Wang, Y., Lyonnet, F., Bacri, J.C., Hocquet, T., Devaud, M., 2015\.  A refined theory of magnetoelastic buckling matches experiments with ferromagnetic and superparamagnetic rods.  Proc. Natl. Acad. Sci. U.S.A. 112, 7135–7140.
  * Ginder et al. (1999) Ginder, J.M., Nichols, M.E., Elie, L.D., Tardiff, J.L., 1999\.  Magnetorheological elastomers: properties and applications, in: Proc. SPIE 3675, Smart Structures and Materials 1999: Smart Materials Technologies, Society of Photo-Optical Instrumentation Engineers, Newport Beach, CA. pp. 131–138.
  * Gu et al. (2020) Gu, H., Boehler, Q., Cui, H., Secchi, E., Savorana, G., De Marco, C., Gervasoni, S., Peyron, Q., Huang, T.Y., Pane, S., Hirt, A.M., Ahmed, D., Nelson, B.J., 2020\.  Magnetic cilia carpets with programmable metachronal waves.  Nat. Commun. 11, 2637\.
  * Gurtin et al. (2010) Gurtin, M., Fried, E., Anand, L., 2010.  The Mechanics and Thermodynamics of Continua.  Cambridge University Press.
  * Hu et al. (2018) Hu, W., Lum, G.Z., Mastrangeli, M., Sitti, M., 2018\.  Small-scale soft-bodied robot with multimodal locomotion.  Nature 554, 81–85.
  * Huang et al. (2016) Huang, H.W., Sakar, M.S., Riederer, K., Shamsudhin, N., Petruska, A., Pané, S., Nelson, B.J., 2016.  Magnetic Microrobots with Addressable Shape Control .  IEEE International Conference on Robotics and Automation ICRA Stockholm, Sweden .
  * Hutchinson (2016) Hutchinson, J.W., 2016.  Buckling of spherical shells revisited.  P. Roy. Soc. Lond. A Mat. 472, 20160577\.
  * Hutchinson et al. (1971) Hutchinson, J.W., Muggeridge, D.B., Tennyson, R.C., 1971.  Effect of a local axisymmetric imperfection on the buckling behavior of a circular cylindrical shell under axial compression.  AIAA Journal 9, 48–52.
  * Hutchinson and Thompson (2017) Hutchinson, J.W., Thompson, J.M.T., 2017\.  Nonlinear buckling interaction for spherical shells subject to pressure and probing forces.  J. Appl. Mech. 84, 061001\.
  * Hutchinson and Thompson (2018) Hutchinson, J.W., Thompson, J.M.T., 2018\.  Imperfections and energy barriers in shell buckling.  Int. J. Solids Struct. 148-149, 157–168.
  * Kim et al. (2012) Kim, J., Hanna, J.A., Byun, M., Santangelo, C.D., Hayward, R.C., 2012.  Designing responsive buckled surfaces by halftone gel lithography.  Science 335, 1201–1205.
  * Kim et al. (2019) Kim, Y., Parada, G.A., Liu, S., Zhao, X., 2019.  Ferromagnetic soft continuum robots.  Sci. Robot. 4, eaax7329.
  * Kim et al. (2018a) Kim, Y., Yuk, H., Zhao, R., Chester, S.A., Zhao, X., 2018a.  Printing ferromagnetic domains for untethered fast-transforming soft materials.  Nature 558, 274\.
  * Kim et al. (2018b) Kim, Y., Yuk, H., Zhao, R., Chester, S.A., Zhao, X., 2018b.  Printing ferromagnetic domains for untethered fast-transforming soft materials.  Nature 558, 274–279.
  * Koiter (1945) Koiter, W.T., 1945.  Over de stabiliteit van het elastisch evenwicht.  Ph.D. thesis. Delft University of Technology. Delft, The Netherlands.
  * Koiter (1969) Koiter, W.T., 1969.  The nonlinear buckling behavior of a complete spherical shell under uniform external pressure, parts i, ii, iii & iv.  Proc. Kon. Ned. Ak. Wet. B72, 40–123.
  * Lazarus et al. (2012) Lazarus, A., Florijn, H.C.B., Reis, P.M., 2012.  Geometry-induced rigidity in nonspherical pressurized elastic shells.  Phys. Rev. Lett. 109, 144301\.
  * Lee et al. (2016) Lee, A., López Jiménez, F., Marthelot, J., Hutchinson, J.W., Reis, P.M., 2016.  The geometric role of precisely engineered imperfections on the critical buckling load of spherical elastic shells.  J. Appl. Mech. 83, 111005\.
  * Lee et al. (2019) Lee, A., Yan, D., Pezzulla, M., Holmes, D.P., Reis, P.M., 2019.  Evolution of critical buckling conditions in imperfect bilayer shells through residual swelling.  Soft Matter 15, 6134–6144.
  * Loukaides et al. (2014) Loukaides, E.G., Smoukov, S.K., Seffen, K.A., 2014.  Magnetic actuation and transition shapes of a bistable spherical cap.  Int. J. Smart Nano. Mater. 5, 270–282.
  * Lum et al. (2016) Lum, G.Z., Ye, Z., Dong, X., Marvi, H., Erin, O., Hu, W., Sitti, M., 2016.  Shape-programmable magnetic soft matter.  Proc. Natl. Acad. Sci. U.S.A. 113, E6007–E6015.
  * Marthelot et al. (2017) Marthelot, J., López Jiménez, F., Lee, A., Hutchinson, J.W., Reis, P.M., 2017.  Buckling of a pressurized hemispherical shell subjected to a probing force.  J. Appl. Mech. 84, 121005\.
  * Moon and Pao (1968) Moon, F.C., Pao, Y.H., 1968\.  Magnetoelastic buckling of a thin plate.  J. Appl. Mech 35, 53–58.
  * Niordson (1985) Niordson, F., 1985.  Shell Theory.  North-Holland Series in Applied Mathematics and Mechanics, Elsevier Science.
  * O’Neill (1997) O’Neill, B., 1997.  Elementary Differential Geometry.  Academic Press.
  * Pancaldi et al. (2020) Pancaldi, L., Dirix, P., Fanelli, A., Lima, A.M., Stergiopulos, N., Mosimann, P.J., Ghezzi, D., Sakar, M.S., 2020\.  Flow driven robotic navigation of microengineered endovascular probes.  Nat. Commun. 11, 1–14.
  * Paulose and Nelson (2013) Paulose, J., Nelson, D.R., 2013\.  Buckling pathways in spherical shells with soft spots.  Soft Matter 9, 8227–8245.
  * Pece et al. (2017) Pece, F., Zarate, J.J., Vechev, V., Besse, N., Gudozhnik, O., Shea, H., Hilliges, O., 2017.  MagTics, in: UIST ’17: The 30th Annual ACM Symposium on User Interface Software and Technology, ACM, New York, NY, USA. pp. 143–154.
  * Pezzulla and Reis (2019) Pezzulla, M., Reis, P.M., 2019\.  A Weak Form Implementation of Nonlinear Axisymmetric Shell Equations With Examples.  J. Appl. Mech. 86.
  * Pietraszkiewicz (1980) Pietraszkiewicz, W., 1980.  Finite Rotations in the Nonlinear Theory of Thin Shells. Springer Vienna, Vienna.  pp. 151–208.
  * Psarra et al. (2019) Psarra, E., Bodelot, L., Danas, K., 2019.  Wrinkling to crinkling transitions and curvature localization in a magnetoelastic film bonded to a non-magnetic substrate.  J. Mech. Phys. Solids 133, 103734\.
  * Rigbi and Jilkén (1983) Rigbi, Z., Jilkén, L., 1983\.  The response of an elastomer filled with soft ferrite to mechanical and magnetic influences.  J. Magn. Magn. Mater 37, 267–276.
  * Roper et al. (2006) Roper, M., Dreyfus, R., Baudry, J., Fermigier, M., Bibette, J., Stone, H.A., 2006\.  On the dynamics of magnetically driven elastic filaments.  J. Fluid Mech. 554, 167–190.
  * Sanders (1963) Sanders, J.L., 1963.  Nonlinear theories for thin shells.  Q. Appl. Math. 21, 21–36.
  * Sano et al. (2021) Sano, T., Pezzulla, M., Reis, P.M., 2021.  A Kirchhoff-like theory for hard magnetic rods undergeometrically nonlinear deformation in three dimensions.  arXiv:2106.15189[cond-mat.soft] .
  * Seffen and Vidoli (2016) Seffen, K.A., Vidoli, S., 2016\.  Eversion of bistable shells under magnetic actuation: a model of nonlinear shapes.  Smart Mater. Struct. 25, 065010\.
  * Timoshenko (1925) Timoshenko, S., 1925.  Analysis of bi-metal thermostats.  J. Opt. Soc. Am. 11, 233–255.
  * Vella et al. (2012) Vella, D., Ajdari, A., Vaziri, A., Boudaoud, A., 2012\.  Indentation of ellipsoidal and cylindrical elastic shells.  Phys. Rev. Lett. 109, 144302\.
  * Wang et al. (2020) Wang, L., Kim, Y., Guo, C.F., Zhao, X., 2020.  Hard-magnetic elastica.  J. Mech. Phys. Solids 142, 104045\.
  * Yan et al. (2021) Yan, D., Abbasi, A., Reis, P.M., 2021.  A comprehensive framework for hard-magnetic beams: reduced-order theory, 3D simulations, and experiments.  arXiv:2106.14878[cond-mat.soft] .
  * Yan et al. (2020) Yan, D., Pezzulla, M., Cruveiller, L., Abbasi, A., Reis, P.M., 2020.  Magneto-active elastic shells with tunable buckling strength.  Nat. Commun. 12, 2831\.
  * Zhao et al. (2019) Zhao, R., Kim, Y., Chester, S.A., Sharma, P., Zhao, X., 2019.  Mechanics of hard-magnetic soft materials.  J. Mech. Phys. Solids 124, 244 – 263.