Eversion of bistable shells under magnetic actuation: a model of nonlinear shapes

We model in closed form a proven bistable shell made from a magnetic rubber composite material. In particular, we incorporate a non-axisymmetrical displacement field, and we capture the nonlinear coupling between the actuated shape and the magnetic flux distribution around the shell. We are able to verify the bistable nature of the shell and we explore its eversion during magnetic actuation. We show that axisymmetrical eversion is natural for a perfect shell but that non-axisymmetrical eversion rapidly emerges under very small initial imperfections, as observed in experiments and in a computational analysis. We confirm the non-uniform shapes of shell and we study the stability of eversion by considering how the landscape of total potential and magnetic energies of the system changes during actuation.


Introduction
In an earlier paper [1], the first author, with others, considered the actuation capabilities of a novel magnetic rubber composite material for developing 'active' shape-changing structures. The material was made into a stress-free, shallow cap, approximately 36 mm across, which was then supported centrally in a convex manner. When a strong magnet was brought closer to the cap along its axis, the cap eventually 'everted' at a separation of 13 mm, and remained inside out after the magnet was removed, see figure 1. The cap shape is clearly uneven during eversion, only becoming axisymmetrical again at the end. A finite element analysis was conducted to confirm the proximity threshold for eversion, the shape of cap during eversion, and the final resting shape. A sequence of these shapes is repeated in figure 1 along with the variation of elastic strain energy over time, with some unusual features that are discussed later.
Finite element analysis reveals behaviour for a specific set-up, which may, or may not, be general. For broader insight, many analyses of different initial geometries are needed but this can be time-consuming, and critical information is not revealed in closed form-essential for further design and development. For this, a theoretical model of large-displacement eversion is required; however, relevant frameworks of the time were inadequate. For example, the assumption of uniform curvature used in studies on multistable shells [2,3] is clearly violated during eversion, whereas the emerging linearly varying curvature (LVC) model [4] considered non-uniform shapes of shell but actuation effects remained decoupled from the deformation they induced. The magnetic field applied to the shell is spatially nonlinear, and the level and distribution of actuation strains inherently depend on the shell displacements. Such strong coupling has not been modelled compactly, to our knowledge, and capturing this behaviour in finite element analysis was, in itself, non-trivial.
In this paper, we build a sophisticated theoretical model with non-uniform displacements and coupled magnetic actuation by augmenting the LVC model from [4], in order to compute the bistable shapes, to elucidate possible eversion 'paths', and to explore how axisymmetry may be 'lost' during eversion. The latter is not readily demonstrated by finite element analysis because the cap and magnetic field are nominally axisymmetrical; there has to be some kind of imperfection present for the simulation to follow an alternative equilibrium path of non-axisymmetrical shapes, see for example [5]. This was achieved w.r.t. figure 1 by imposing a small initial tilt of the cap about its pole relative to the magnetic field axis because of inevitable axial misalignments in practice. However, it is well known in perfectly symmetrical systems that symmetry-breaking behaviour can be the norm, for example in beam arches and spherical caps under axial point loads [6] and in vibrating shells [7]. In our model, we build this capability into the displacement field, which allows us to ascertain whether or not symmetry is lost during actuation under perfect conditions and to understand the influence that actuation has upon this assumption alongside the usual competing effects of bending and stretching in ordinary shells. The layout of study is as follows.
We reformulate the LVC model using a displacement field with axisymmetrical and antisymmetrical components. Their respective amplitudes are taken as control parameters and their expressions always ensure that the edge of the cap is free from tractions. We follow the recognised procedure developed originally in [2] and enhanced in [3] for finding equilibrium shapes: we find the in-plane stresses for the deformed cap by solving Gauss' compatibility equationidentical to one of the Föppl-von Kármán shallow plate equations; using these stresses and knowing the deformed curvatures, we can compute the total elastic strain energy stored in the deformed cap as functions of the control parameters. The effect of the applied magnetic field is incorporated as a potential energy function, which is added to the elastic components to yield the total energy of the system. By exploring the total energy 'landscape' for a given proximity of magnet, we can establish the nature of equilibria and how they change as the magnet moves. We do not include dynamic or inertial behaviour because that is beyond our scope; our focus is the symmetry character of eversion, which is the same for all time and thus, is defined by the quasi-static response of the shell. Correspondingly, we verify that the shell in figure 1 is bistable and that its everted shape is non-uniform. We show that the magnetic field, if strong enough, is able to render the natural configuration unstable i.e. to initiate eversion. We capture the everted shape of shell when the magnet is close enough so that contact takes place between the two. We confirm the general strain energy profile from finite element analysis shown in figure 1, and that the eversion threshold is nonlinearly dependent on the initial title angle; moreover, we predict the same threshold as per finite element analysis. We then make some general conclusions.

Kinematical assumptions
Consider a thin shallow shell with a circular planform  of radius R, planform area p = A R 2 , and constant thickness t e . Using cylindrical coordinates ( ) r f z , , , the natural stress-free configuration,  0 , is identified by the shell mid-surface coordinate, w 0 , in the transverse direction, ζ, with . The cap is 1.5 mm thick, approximately 36 mm wide, and is held at its centre. The magnet is moved closer to the cap and, at an axial distance of around 13 mm from the cap centre, the cap everts unevenly-without axisymmetry-before making contact with the magnet. When the magnet is retracted, the cap remains everted. Left bottom and right: finite element analysis shows a non-axisymmetrical transition for the cap when it is initially tilted by 5°r elative to the magnetic axis: the time taken is around 0.1 s. The plot indicates the elastic strain energy of deformation during the transition for two axial distances: d=13 mm, which results in eversion; d=14 mm, which does not. The feature highlighted in the red box is confirmed later.
This is a doubly curved, axisymmetrical shape with uniform and positive Gaussian curvature, R 1 c 2 , as per the stress-free cap in figure 1. If the shell is sufficiently curved and thin, it has one other stable, everted configuration, at least. A simple condition for bistability is derived in [4] for shallow isotropic shells: that the dimensionless ratio, ≔ l R t A c e , is less than a critical ratio, l  0.04 c . We also assume material isotropy and linear elastic behaviour.
Every shell configuration,  must satisfy the following boundary conditions on its edge at r = R: As the scalar parameter s varies,ˆ( ) r w s describes a deformed axisymmetrical shell with a free edge: clearly, the natural configuration, equation (2), is recovered when s = 0. For > s 0, eversion is deemed to be occurring when the curvature at the cap centre (r = 0), given bŷ Antisymmetrical distortions are addressed by complementing equation (4) with a second ansatz of circumferential variation. We find the simplest form to be: Hence, the general shape function satisfies the boundary conditions, equation (3), for every choice of s and t. Note that most shapes parameterised by equation (7) do not have a uniform curvature field.

Bending and stretching energies
For a deformed Föppl-von Kármán shell, the stored elastic strain energy is the sum of two quadratic components in bending,  b , and in stretching (or membranal),  s . For the bending energy, see for instance [8], we have where= w K is the standard curvature field of deformed shell, equation (2), = w K 0 0 is the curvature field of the natural stress-free shape, equation (1), and  is the bending stiffness matrix for an isotropic material: Y is the Young modulus and ν is the Poisson ratio. Using equation (7), the standard components of K arê( Note that ( ) vanishes for the natural configuration where = = s t 0. Denoting the membranal stress field via N, the stretching strain energy has the general expression: where  is the membranal stiffness defined in equation (9), see [8]. The membranal stress field can be found by solving Gauss' original compatibility equation for shallow shell displacements because we already know the Gaussian curvature properties of the shell; specifically, the difference in Gaussian curvature between the current and initial shapes constitutes a 'forcing' function for Gauss' equation. The midsurface strains are substituted by stresses reformulated in terms of an Airy function j for isotropic materials, and Gauss' equation, in its form defined in [4], now reads as with no in-plane tractions on the free edges demanding that Here, j DD indicates the standard biharmonic operator on j, and ¶ ¶n indicates derivation with respect to the direction n, the outward normal to  ¶ . The usual membranal stress components per unit length are recovered, see for instance [9], from: The forcing term in equation (15) can be computed in closed form from the generic shape variation in equation (7). This returns The generic Airy function is correspondingly written , and, using equations (14) and (17), the stretching strain energy can be computed to be

Total elastic energy: evidence of bistability
The total elastic energy is the sum of the two components and is a fourth order polynomial expression in s and t. We can solve ( ) for stationary values of s and t, in order to find the range of equilibria. We can also explore the variation of ( )  s t , Two energy minima corresponding to stable equilibria are clearly evident: the natural, stress-free state is located at ( ) ( ) = s t , 0,0, and the everted configuration at  Elsewhere on the landscape, the deformed shell is not stable and must 'move' towards either stable equilibria along a path of steepest descent-a longstanding principle embodied, for example, in the behaviour of chemical reactions [10]. When actuating forces are present, the landscape itself is not static but is modified by the work contribution from actuation-as described in, for example, [11,12]. We must therefore consider the total energy landscape including the magnetic potential energy, expressed in terms of the shape of shell, i.e. as a function of s and t.   By gradually increasing the level of actuation-by moving the magnet closer to the cap, we may observe how the stability character of the original elastic strain energy landscape changes. Eversion begins when the natural configuration loses its stability margin and becomes unstable in this landscape, leaving only the stable everted configuration. This is an important point, which is not usually noted in studies of actuator effectiveness, but which has been discussed thoughtfully in [13]. In becoming unstable, the landscape around the natural state changes from a well to a saddle, and eversion follows the path of steepest descent along a direction in which the local gradient first becomes negative. As noted before, we neglect inertial effects because this adds further modelling complexity; importantly, the properties of eversion, such as symmetry and threshold proximity, do not originate in the dynamics of response.
Before formally calculating the magnetic potential energy, it is worth pondering how eversion may proceed without axisymmetry. The magnetic field is axisymmetrical in nature and, therefore, exerts symmetrical forces upon the cap initially, provided the cap is itself perfectly axisymmetrical and isotropic. Although the actuator effect is absent from the landscape in figure 2, axisymmetrical deformation proceeds along t=0, along a path of steepest ascent and descent towards the everted shape, as highlighted. Either side of the ascending path, there are shape configurations for which  e is immediately lower and, intuitively, this part of the path not be globally stable. If, on the other hand, we move along a different path, it may be possible that local perturbations from  e are stable, giving a more secure route for eversion; one such route is also highlighted, which fulfils this property by rising least steeply away from the natural state, where the distance between successive contours is largest. Since t is now non-zero, the intermediate shapes of shell are non-axisymmetrical, as shown. If we now compute the following two energetic quantities: , signifying that axisymmetrical eversion is highly likely; a quantitative answer is only conveyed when we understand how actuation and imperfections distort the overall energy landscape, which is now performed.

Magnetic actuation
The magnetic field in [1] is generated by connecting together high-strength Neodymium N42 disk magnets along their axis into a cylindrical stack for a more intense field of flux. Calculating the resulting field is not trivial but we follow the scheme outlined in [1], originally based on the derivation in [14].
The magnetic field, B, is assumed to be axisymmetrical, and to vary with radial distance, r, and axial distance, z, from a reference point source. More formally is the magnetic potential originating from two parallel, oppositely charged monopolar disks of radius r=a separated by a perpendicular distance H, equivalent to the height of stack [14]. Here (·) E denotes the complete elliptic integral; the top of the stack lies within the plane z=0 and the bottom in the plane =z H; the vacuum permeability is m 0 and the magnetisation per unit volume of magnet material is M.
A material particle within the field, B, is subjected to a conservative force m m Figure 5. Spatial distribution of the potential energy ( ) U r z , m of the cylindrical magnetic field defined by the properties in table 2. r and z are orthogonal cylindrical coordinates with origin in the centre of the magnet. r and z are local shell coordinates with origin located a distance, d, from the magnet; q is the initial angle of tilt of the cap, whose initial cross-section is shown.
where ( ) m B is the magnetization induced by the magnetic field and U m is the potential energy associated with the magnetic field. It has been experimentally observed [15], that the magnetisation of a material particle of volume dV is well fitted by the following linear relation where dc is a dimensionless number indicating the magnetic susceptibility of the material composing the shell [16]. Using this approximation, the potential energy of the magnetic forces acting on every infinitesimal volume dV of shell is Equation (26) must be integrated over the volume of shell, in order to obtain the total magnetic energy associated with the shell. We may use the original shell cylindrical coordinates from section 2, resulting in However, since we want the shell to move inside the reference frame ( ) r z , of the magnet, we introduce the following change of coordinates where (·) h is the unit step function. Here, r * and z * indicate the coordinates of a generic point of the shell; d is the proximity distance between the centre of the shell and the upper surface of the magnet; and q is the angular rotation of the z axis (attached to the shell) with respect to the z (attached to the magnet), i.e. q is the tilt angle of the shell, see figure 5.
, t s, must be computed according to the ansatz equation (7), which will lead to an explicit dependance of the magnetic energy on the actual shell shape. The unit step function h is used when we want to avoid interpenetration of the magnet and shell by ensuring *  z 0.
Finally, we obtain the total potential energy of the shell as the following integral The energy functional,  m , depends on the shell configuration through the parameters (s, t) and the relative position and orientation through d and θ. The function has to be evaluated numerically: using a standard CPU Core at 2.8 GHz and WOLFRAM MATHEMATICA [17], 100 evaluations takes about 15 s. The magnetic properties and geometry of the magnet used in [1] are given in the following table.

Total energy
The total energy  of the actuated shell is the sum of elastic and magnetic energy components: Clearly, at large distances from the magnet (  ¥ d ), the magnetic energy is negligible and  reduces to the elastic strain energy landscape in figure 2 with two stable equilibria: this is the 'far-field' case. When d is reduced, the magnetic contribution increases, and  begins to differ from  e , which, in turn, affects the stability of the original equilibria. By way of example, we plot  in figure 6 as a function of s only, stipulating t=0 and q = 0, for three values of d above the experimental threshold of eversion: 1 m, 20 mm and 14 mm. For the largest d, the equilibria are unaltered see figure 3, whereas for = d 20 mm, the magnet is close enough to distort the initial state, observed as a slight movement of the original minimum along the s-axis. For the smallest d, the curve passing through the original minimum is now noticeably more rotated as this stationary point transforms into a point of inflexion. All everted configurations are local minima and always stable, and their locations along s increase beyond 0.995, indicating that they become more distorted than the far-field shape. Negative values of  are justified since we choose the initial configuration to have a vanishing energy in the far-field case.  If d is decreased further, the initial state becomes a point of inflexion and unstable in the direction > s 0. The system will naturally move along the path of steepest descent towards the stable everted configuration. However, we also need to ascertain the landscape elevation around this path when t is non-zero in view of the propensity for non-axisymmetrical eversion. This is now considered for caps with and without initial tilt angles.

Vanishing tilt angles: axisymmetrical eversion
For vanishing tilt angles, q = 0, it can be shown that the total energy,  , is symmetrical with respect to t=0, i.e.
By definition, the corresponding landscape is symmetrical about the t=0 axis, which heavily favours axisymmetrical eversion. We can test this hypothesis by setting d=7 mm so that the magnet is closer than required for eversion. The total energy is plotted in figure 7, where the original minimum at ( ) ( ) = s t , 0,0 is now a saddle and has lost its stability. The system will move away from this towards the stable everted minimum around s=1.25. Because of symmetry about t=0, the eversion path follows a decreasing ridge and the transition is axisymmetrical.
The spacing between points on the path has been calculated by solving the auxiliary equation where ∇ is the gradient function with respect to s and t and an overdot denotes differentiation with respect to time. This equation is a statement of damping in the direction of steepest descent and is analogous to the artificial 'numerical' damping      figure 1.
17.32 mm 18 mm 1.5 mm 750 kPa 0.5 Table 2. Properties of the magnet in figure 1.
p´-4 10 7 N A −11 .03 10 6 A m −1 0.20 25 mm 12 mm used in finite element analysis for achieving convergent solutions: when there are rapid changes in values, damping is higher and solution points are more finely resolved for the sake of accuracy. The equation above therefore provides a way of comparing to the time-oriented finite element data even though our framework does not include other dynamic effects; equation (31) is solved for fixed intervals of time so that the spacing between points is directly proportional to the magnitude of the local gradient. After calculating the (s, t) coordinates of each point, we extract the corresponding values of the elastic strain energy,  e , from the total landscape, which is rendered in figure 8(a). These values are then plotted for each time interval in figure 8(b). This profile, encouragingly, replicates the features present in the finite element output of figure 1, in particular, the highlighted second minimum. Strictly speaking, we need to consider nonaxisymmetrical eversion from our analysis, which is performed in the next section; however, these preliminary results give confidence to our solution methodology in view of more sophisticated multi-physics analysis software. The corresponding shapes of shell are also plotted in figure 8(c).

Non-vanishing tilt angles: non-axisymmetrical eversion
Symmetry conditions are not imposed upon the performance of the total energy,  , and the initial conditions specify d=13 mm and q =  5 , as in [1]. The landscape of  is plotted in figure 9 and, again, we solve equation (31) over fixed intervals of time. Again, the stability of the natural minimum near ( ) ( ) = s t , 0,0 is lost, which becomes a saddle; however, it is skewed in the direction of positive t, and the eversion path being directed away from the s-axis is therefore one of non-axisymmetrical shapes. When contact between the deformed shell and magnet is enforced by the condition of * > z 0 in equation (27), this path bifurcates as shown and follows a different eversion route. This diversion happens early on because the shell is more deformed on one side and makes contact well before full eversion. We see this also in the experiment in figure 1. Note that we are not modelling the contact forces, so the evaluation of the elastic strain energy after contact has been made may not be entirely correct. Both sets of path points are superposed onto the elastic strain energy landscape in figure 10(a), in order to construct the time-evolution profiles in figure 10(b). Again, the features compare well to the finite element output in figure 1. The deformed shapes for the contacting case are displayed in figure 10(c).

Critical values of the distance
We finally address the dependence of the eversion threshold, d, upon θ. We tackled this by performing a numerical continuation of the natural equilibrium branch for several tilt angles, see for instance [18]. For any θ, this branch originates from the point ( ) ( ) = s t , 0,0 as  ¥ d . The results are shown in figure 11 and can be summarised as follows.
• For a given angle of tilt, there exists a well defined critical value of distance, denoted as ( ) q d c , below which stability of the equilibrium branch is lost; • the function ( ) q d c for small tilt angles | |  q  15 monotonically increases from a minimum value of around 7 mm for q = 0; • a non-vanishing tilt angle 'encourages' eversion by enabling non-axisymmetrical shapes. This is evident from figure 11(a), where the solution branches are plotted in the three-dimensional space ( ) s t d , , : these are characterized by positive values of t well before their stability is lost; • a relevant sensitivity of the function ( ) q d c is evident near q = 0; one could expect the values  d 11 mm c to be rarely observed as any small imperfections in practice would lead to a non-axisymmetrical eversion.
Thus, in a perfect world, the proximity distance for eversion is around 50% of the value of 13 mm observed in experiments and in finite element analysis. This value quickly rises as the tilt increases, reaching 13 mm when θ is just three degrees. Thereafter, the rate of increase diminishes, so does the influence of the tilt angle, and eversion is nonaxisymmetrical.

Conclusions
We have developed a two degree-of-freedom theoretical model for predicting the general shapes of a linear elastic shell cap made of novel material, which ultimately everts under a magnetic field. We have advanced the general formulation for such problems by explicitly dividing the displacement field into axisymmetrical and antisymmetrical components, so that we may explore the behaviour directly under both assumptions. We have included magnetic actuation for the first time in such a formulation, and we have focussed on the total energy profile of the system, in order to understand its stability characteristics and how they change during actuation. Importantly, we have argued that because we have a closed system, actuation from the initial state will follow the path of maximum descent from this position. We have been able to interrogate this graphically via the landscapes we generate, for a more intuitive insight. We have obtained good agreement with results from a leading study, in particular, we have confirmed the proximity threshold for eversion, we have shown that an initial title angle is conducive for non-axisymmetrical eversion, and we have determined that the critical threshold for eversion is strongly controlled by the tilt angle when the angle is small. The goodness of these qualitative comparisons give us confidence about our model's overall robustness, which may then be used to predict detailed information such as transverse displacements, materials stresses etc throughout the cap, and furthermore, to describe the performance of other shapes of magnetic shell. We therefore hope our approach may inspire others to pursue modelling in this way in concert with computational studies.