Wrinkling of soft magneto-active plates

Coupled magneto-mechanical wrinkling has appeared in many scenarios of engineering and biology. Hence, soft magneto-active (SMA) plates buckle when subject to critical uniform magnetic field normal to their wide surface. Here, we provide a systematic analysis of the wrinkling of SMA plates subject to an in-plane mechanical load and a transverse magnetic field?. We consider two loading modes: plane-strain loading and uni-axial loading, and two models of magneto-sensitive plates: the neo-Hookean ideal magneto-elastic model and the neo-Hookean magnetization saturation Langevin model. Our analysis relies on the theory of nonlinear magneto-elasticity and the associated linearized theory for superimposed perturbations. We derive the Stroh formulation of the governing equations of wrinkling, and combine it with the surface impedance method to obtain explicitly the bifurcation equations identifying the onset of symmetric and antisymmetric wrinkles. We also obtain analytical expressions of instability in the thin- and thick-plate limits. For thin plates, we make the link with classical Euler buckling solutions. We also perform an exhaustive numerical analysis to elucidate the effects of loading mode, load amplitude, and saturation magnetization on the nonlinear static response and bifurcation diagrams. We find that antisymmetric wrinkling modes always occur before symmetric modes. Increasing the pre-compression or heightening the magnetic field has a destabilizing effect for SMA plates, while the saturation magnetization enhances their stability. We show that the Euler buckling solutions are a good approximation to the exact bifurcation curves for thin plates.


Introduction
Soft magneto-active (SMA) materials, such as magneto-active elastomers, are a particularly promising kind of smart materials that can respond to magnetic field excitation.In general, these SMA materials are prepared by mixing micron-sized magnetizable particles (such as carbonyl iron and neodymium-iron-boron) into a non-magnetic elastomeric matrix such as rubber or silicone (Rigbi and Jilken, 1983;Ginder et al., 2002;Kim et al., 2018).They can then deform significantly under the simple, remote, and reversible actuation of magnetic fields.Moreover, their overall magneto-mechanical properties can be altered actively by applying suitable magnetic fields, thus resulting in tunable vibration and wave characteristics.
A challenging problem that arises in the study of SMA materials is how to model the strong nonlinearity and the magneto-mechanical coupling.Thus, a lot of academic interest has been devoted over the years to establish a general theoretical framework of nonlinear magneto-(visco)elasticity in order to describe appropriately the magneto-mechanical response of SMA materials (Tiersten, 1964;Brown, 1966;Pao and Yeh, 1973;Maugin, 1988;Brigadnov and Dorfmann, 2003;Dorfmann and Ogden, 2004;Bustamante, 2010;Destrade and Ogden, 2011;Saxena et al., 2013).Comprehensive reviews regarding the theoretical development of nonlinear magneto-elasticity include those by Kankanala and Triantafyllidis (2004) and Dorfmann and Ogden (2014).Furthermore, homogenization techniques based on micromechanical methods have also been developed to understand the connection between the magneto-active microstructures and the macroscopic physical or mechanical properties of SMA materials (Castañeda and Galipeau, 2011;Galipeau et al., 2014).
In a large variety of practical applications, the mechanical and magnetic loads (prestretch, magnetic field, etc.) influence the working performance of smart systems made of SMA materials and may lead to instability and even failure.In fact, it has long been observed that a magneto-elastic beam or plate in a uniform transverse magnetic field will buckle once the field reaches a critical value (Moon and Pao, 1968;Miya et al., 1978;Gerbal et al., 2015).
On the other hand, local buckling and other instability phenomena can be exploited to realize active pattern switching devices and reconfigurable metamaterials.Hence Psarra et al. (2017Psarra et al. ( , 2019) ) made use of wrinkling and crinkling instabilities of a thin SMA film bonded on a soft non-magnetic substrate to achieve surface pattern control through a combined action of magnetic field and uni-axial pre-compression.Goshkoderia et al. (2020) investigated experimentally instability-induced pattern evolutions in SMA elastomer composites driven by an applied magnetic field.
Thus, it is vital to theoretically study the stability of SMA structures and composites subject to coupled magneto-mechanical loads, so that we can provide solid guidance for simulations and experiments.From a theoretical point of view, the classical buckling problem of a magneto-elastic beam-plate was first addressed by Moon and Pao (1968), based on the thin-plate theory and the assumption of a linear ferromagnetic material.Pao and Yeh (1973) used a general theory of magneto-elasticity to re-examine this problem and to yield an identical antisymmetric buckling equation for thin plates.Following those works, many investigations looked at the same problem, trying to improve mathematical models to explain the discrepancy between experimental results and theoretical predictions (Wallerstein and Peach, 1972;Popelar, 1972;Dalrymple et al., 1974;Miya et al., 1978;Gerbal et al., 2015;Singh and Onck, 2018).Most of the above-mentioned works employed classical structural models to elucidate the magneto-mechanical coupling problem.More recently, using the theory of nonlinear magneto-elasticity, some researchers have explored the onset of instabilities of different SMA structures, including surface instabilities of isotropic SMA half-spaces (Otténio et al., 2008), buckling modes of rectangular SMA blocks undergoing plane-strain loading (Kankanala and Triantafyllidis, 2008), macroscopic instabilities of anisotropic SMA multilayered structures (Rudykh and Bertoldi, 2013), instabilities of a thin SMA layer resting on a soft non-magnetic substrate (Danas and Triantafyllidis, 2014), and instabilities of a cylindrical membrane (Reddy and Saxena, 2018).
The present work revisits the stability problem of SMA plates subject to an in-plane mechanical load and a uniform transverse magnetic field, and evaluates explicitly the onset of wrinkling instabilities.This work differs from previous works (Pao and Yeh, 1973;Kankanala and Triantafyllidis, 2008) in the following respects.(i) We adopt the nonlinear magneto-elasticity theory and the associated linearized theory developed by Dorfmann and Ogden (2004) and Otténio et al. (2008) to derive the governing equations of nonlinear static response and the bifurcation equations of wrinkles.These theories introduce a total stress tensor and a modified (or total) energy density function to express constitutive relations in a simple and compact form.(ii) Two mechanical loading modes are considered: plane-strain loading (Fig. 1(b)) and uni-axial loading (Fig. 1(c)).(iii) To overcome the complexity of conventional displacement-based method, we employ the Stroh formulation and the surface impedance method (Su et al., 2018) to obtain explicit expressions of the bifurcation equations of antisymmetric and symmetric wrinkling modes (Fig. 1(d) and 1(e)).(iv) We also manage to derive the explicit bifurcation equations corresponding to the thin-and thick-plate limits, and to establish the thin-plate approximate formulas.The paper is organized as follows.The equations of nonlinear magneto-elasticity are derived in Sec. 2. The linearized stability analysis is conducted in Sec. 3 to obtain the exact bifurcation equations for symmetric and antisymmetric wrinkles.Section 4 specializes the analytical expressions to the neo-Hookean magnetization saturation material model, including the ideal magneto-elastic model.Numerical calculations are carried out in Sec. 5 to illustrate the effects of loading mode, load amplitude, and saturation magnetization on nonlinear static response and bifurcation diagrams.Section 6 concludes the work with a summary.Some mathematical expressions or derivations are provided in Appendices A-C.

Finite magneto-elasticity theory
We first present the basic equations governing the finite magneto-elastic deformations of an incompressible soft magneto-elastic body, as developed by Dorfmann and Ogden (2004).
In the undeformed stress-free state, the body occupies in the Euclidean space a region B r (the boundary being ∂B r , with the outward normal N), which is taken to be a reference configuration.Any material point inside the body is then identified with its position vector X in B r .The application of both mechanical load and magnetic field deforms quasistatically the body from B r to the current configuration B (the boundary being ∂B, with the outward normal n).The particle located at X in B r now occupies the position x = χ(X) in B, where the vector function χ is a one-to-one, orientation-preserving mapping with a sufficiently regular property.The associated deformation gradient tensor is F = Grad χ = ∂χ/∂X (F iα = ∂x i /∂X α in Cartesian components) with Grad being the gradient operator in B r .
Note that Greek indices are associated with B r and Roman indices with B.
The Jacobian of the deformation gradient is J = det F, and it is equal to 1 at all times for incompressible materials.The left and right Cauchy-Green deformation tensors b = FF T and C = F T F are used as the deformation measures, where the superscript T signifies the transpose.
In the Eulerian description, the equilibrium equations in the absence of mechanical body forces, and the Maxwell equations in the absence of time dependence, free charges and electric currents, are where τ is the total Cauchy stress tensor including the contribution of magnetic body forces, B is the Eulerian magnetic induction vector and and H is the Eulerian magnetic field vector; curl and div are the curl and divergence operators in B, respectively.The conservation of angular momentum implies that τ is symmetric.
The Eulerian magnetization vector M is defined by the standard relation where µ 0 = 4π × 10 −7 N • A −2 is the magnetic permeability in vacuum.
In vacuum, there is no magnetization and Eq. ( 2) reduces to B = µ 0 H , where a star superscript identifies a quantity exterior to the material.The Maxwell stress tensor τ in vacuum is where I is the identity tensor.The external fields B and H satisfy div B = 0 and curl H = 0, which leads to div τ = 0.
The traction and magnetic boundary conditions are written in Eulerian form as where t a is the applied mechanical traction vector per unit area of ∂B.Note that the magnetic traction vector t m induced by the external Maxwell stress τ is t m = τ n.Using Nanson's formula, we can transform the governing equations (1) and boundary conditions (4) into the Lagrangian form.
Following the theory of nonlinear magneto-elasticity (Dorfmann and Ogden, 2004), it is convenient to express the nonlinear constitutive relations for incompressible magneto-elastic materials in terms of a total energy function, or modified free energy function, Ω (F, B L ) per unit reference volume as where T = F −1 τ is the total nominal stress tensor, the nominal magnetic field vector B L = F −1 B and the nominal magnetic induction vector H L = F T H are the Lagrangian counterparts of B and H, respectively, and p is a Lagrange multiplier related to the incompressibility constraint, which can be determined from the equilibrium equations and boundary conditions.The Eulerian counterparts of Eq. ( 5) are For an incompressible isotropic material, the energy function Ω can be reduced to a function depending only on the following five invariants: According to Eqs. ( 6) and ( 7), the total Cauchy stress tensor τ and the Eulerian magnetic field vector H are then derived as where Ω n = ∂Ω/∂I n (n = 1, 2, 4, 5, 6).

Pure homogeneous deformation of a plate
We now specialize the above equations of nonlinear magneto-elasticity to the homogeneous deformation of an SMA plate subject to an in-plane mechanical load and a uniform magnetic field along its thickness direction.In this work, we consider two mechanical loading modes: plane-strain loading (see Fig. 1(b)) and uni-axial loading (see Fig. 1(c)).
We point out that to simplify the mathematical modelling and obtain explicit analytical solutions to the nonlinear response and the bifurcation criteria of SMA plates, we make the following assumptions.(i) The plate is entirely immersed in the external magnetic field and subjected to an applied uniform magnetic field in the thickness direction.(ii) The plate is of infinite lateral extent, i.e., the plate lateral dimensions are much larger than its thickness.
(iii) The magnetic poles are very close to the plate, see Fig. 1(b) and 1(c) for sketches and Fig. 1(f) for the experimental apparatus of Psarra et al. (2017).These assumptions ensure that the field heterogeneity is only confined to the plate edges, while far away from the edges, the induced stress and magnetic field distributions can be envisioned as uniform, which is a good approximation in the sense of Saint-Venant's principle.A similar processing method has been adopted by Su et al. (2020) to analyze the wrinkling of soft electro-active plates immersed in external electric fields.
Let (X 1 , X 2 , X 3 ) and (x 1 , x 2 , x 3 ) be the Cartesian coordinates in the reference and current configurations, respectively, along the (E 1 , E 2 , E 3 ) and (e 1 , e 2 , e 3 ) unit vectors.In the undeformed configuration, the plate has uniform thickness 2H along the X 2 direction and lateral lengths L 1 , L 3 in the (X 1 , X 3 ) plane.The pure homogeneous deformations corresponding to the two loading modes are defined by x 3 = λ 3 X 3 , where λ i (i = 1, 2, 3) is the principal stretch in the X i direction and λ 2 = (λ 1 λ 3 ) −1 by incompressibility.The resulting deformation gradient tensor is the diagonal matrix in the e i ⊗ E j basis.The current thickness and lateral lengths of the deformed plate are 2h, 1 and 3 , respectively.
We take the underlying Eulerian magnetic induction vector B to be in the The five independent invariants in Eq. ( 7) are written now as Thus, we can define a reduced energy function ω as ω (λ 1 , λ 3 , I 4 ) = Ω (I 1 , I 2 , I 4 , I 5 , I 6 ).
It follows from the magnetic boundary conditions (4) 2,3 that the non-zero components of B and H are B 2 = B 2 and H 2 = µ −1 0 B 2 , respectively.Thus, the non-zero components of the Maxwell stress (3) are These Maxwell stress components generate the magnetic traction vector t m = τ n.
For uni-axial mechanical loading in the x 1 direction (Fig. 1(c)), there are no mechanical tractions applied on the faces x 2 = ±h and x 3 = ± 3 /2, only magnetic tractions.The traction boundary conditions (4) 1 yield where s 1 is the nominal mechanical traction per unit area of ∂B r applied on the faces x 1 = ± 1 /2.Thus, we deduce from Eqs. ( 10) and ( 13) that the governing equations of the nonlinear response for uni-axial loading are Note that Eq. ( 14) 2 determines the induced principal stretch λ 3 in terms of the stretch λ 1 and magnetic field B 2 .Then s 1 and H 2 are calculated from Eqs. (14) 1,3 , respectively.

Linearized stability analysis
In this section we employ the linearized incremental theory of magneto-elasticity (Otténio et al., 2008;Destrade and Ogden, 2011) and the Stroh formalism (Su et al., 2018) to investigate the formation of small-amplitude wrinkles, signaling the onset of wrinkling instability of the SMA plate, for the two loading modes described in Sec.2.2.

Incremental governing equations
Consider an infinitesimal incremental mechanical displacement u = ẋ along with an updated incremental magnetic induction vector ḂL0 , superimposed on the finitely deformed configuration reached via x = χ(X).Here and henceforth, a superposed dot indicates an increment.
The incremental balance laws and the incremental incompressibility condition are formulated in Eulerian, or updated Lagrangian, form as where Ṫ0 = F Ṫ, ḂL0 = F ḂL and ḢL0 = F −T ḢL are the push-forward versions of the corresponding Lagrangian increments Ṫ, ḂL and ḢL , respectively, and L = grad u is the incremental displacement gradient.The resulting push-forward quantities are identified with a subscript 0.
The linearized incremental constitutive equations for incompressible SMA materials read where A 0 , Γ 0 and K 0 are, respectively, fourth-, third-and second-order tensors, which are referred to as instantaneous magneto-elastic moduli tensors (see Otténio et al. (2008) and Destrade and Ogden (2011) for their general expressions).In index notation, these magnetoelastic moduli tensors are given by where A, Γ and K are the relevant referential magneto-elastic moduli tensors, which are the instantaneous moduli tensors have the symmetries A 0piqj = A 0qjpi , Γ 0piq = Γ 0ipq , and Using the incremental form of the rotational balance condition FT = (FT) T , we have the following relation between A 0 and τ for an incompressible material The incremental fields exterior to the material read where Ḃ and Ḣ are to satisfy div Ḃ = 0 and curl Ḣ = 0, respectively, and hence τ * is divergence-free, i.e., div τ = 0.
The mechanical and magnetic boundary conditions for the incremental fields can be expressed, in updated Lagrangian form, as where ṫA 0 da = ṫA dA, with t A being the applied mechanical traction vector per unit area of ∂B r (i.e., t A dA = t a da).Note that dA and da are area elements of the reference and current configurations, respectively.
In this work we focus on incremental two-dimensional solutions independent of x 3 , such that u 3 = ḂL03 = ḢL03 = 0, and hence From Eq. ( 20) 3 , the incremental magnetic field vector ḢL0 is curl-free and thus an incremental magnetic scalar potential φ can be introduced such that ḢL0 = −grad φ, with components The incremental balance laws and incompressibility condition (20) 1,2,4 become For pure homogeneous deformation of the SMA plate subject to the transverse magnetic field, we have F ij = 0 for i = j and B 1 = B 3 = 0.The magneto-elastic moduli tensors A 0 , Γ 0 and K 0 satisfy (Otténio et al., 2008) Consequently, using Eqs.( 26) and ( 28), the incremental constitutive relations (21) are written in terms of u i (i = 1, 2) and φ as and where the effective material parameters c ij , e ij and µ ij are defined as

Stroh formulation and its resolution for plates
The two-dimensional incremental solutions, for the wrinkling instability with sinusoidal shape along the x 1 direction and amplitude variations along the x 2 direction (see Fig. 1(d) and 1(e)), are sought in the form where i = √ −1, U 1 , U 2 , Φ, Σ 21 , Σ 22 and ∆ are non-dimensional functions of kx 2 only, k = 2π/L is the wrinkling wavenumber with L being the wavelength of the wrinkles, and G is the initial shear modulus (in Pa) of the SMA plate in the absence of magnetic field.
Following a standard derivation procedure (Su et al., 2018(Su et al., , 2019)), we can rewrite the incremental governing equations ( 27), ( 29) and ( 30) in the following non-dimensional Stroh form: where For the constant Stroh matrix N considered here, we seek solutions to Eq. ( 33) in the form, Substituting Eq. ( 34) into Eq.( 33) yields an eigenvalue problem N − qI η 0 = 0, with eigenvalue q and eigenvector η 0 .The requirement of non-trivial solutions results in det N − qI = 0, which gives the following bi-cubic characteristic equation in q: It is clear that the characteristic equation and hence the eigenvalues do not depend on τ 22 (appearing in the Stroh matrix (A.1)) for any choice of energy density function, despite the presence of external magnetic field.
So far, there is no restriction on the specific form of energy density function Ω.To proceed further and demonstrate the possibility of obtaining concise and analytical eigen-equation and eigenvectors, we henceforth consider the magneto-elastic material to be characterized by the following Mooney-Rivlin magneto-elastic model, where F is an arbitrary function of I 5 only and β ∈ [0, 1] is a constant.Note that Eq. ( 37) with β = 0 reduces to the neo-Hookean magneto-elastic model, including, for example, the ideal model with no saturation and the magnetization saturation Langevin model, which we consider in Sec. 4.
Using Eqs. ( 37) and (A.5), the dimensionless parameters a − g in Eq. (A.2) are calculated as where B 2 = B 2 / √ Gµ 0 is the dimensionless applied magnetic induction; F 5 and F 55 are defined as Then we find that the characteristic equation ( 35) factorizes as which yields six pure imaginary eigenvalues as where the real numbers p 1 , p 2 , and p 3 are The associated eigenvectors are derived as where the asterisk * denotes the complex conjugate.We observe from Eq. ( 43) that where η (j) i is the i-th component of the eigenvector η (j) .

Incremental boundary conditions
We take the applied mechanical traction t A as a dead load (i.e., ṫA 0 = ṫA = 0), so that the general incremental boundary conditions ( 25) are specialized to the two-dimensional problem at x 2 = ±h, as From curl Ḣ = 0, we deduce the existence of an incremental magnetic scalar potential Substituting Eq. ( 46) 3,4 into div Ḃ = 0 results in the Laplace equation for φ , To satisfy the decay condition φ → 0 as x 2 → ±∞, we take the solutions to Eq. ( 47) which are localized near the interfaces x 2 = ±h, as where A + and A − are arbitrary constants.Thus, the associated incremental Maxwell stress tensor (24) 2 has non-zero components and for x 2 > h and x 2 < −h, respectively.
Similarly, at the plate bottom surface x 2 = −h, we obtain and where the surface impedance matrix in the lower half-space is
In the thick-plate or short-wave limit (kH → ∞), the tanh functions in Eq. ( 64) are replaced by 1 and the bifurcation criteria for both symmetric and antisymmetric modes reduce to Note that the bifurcation equation ( 65) identifies the surface wrinkling instability for the magneto-elastic half-space, which can also be derived based on the surface impedance method shown in Appendix B. For B 2 = 0, Eq. ( 65) reduces to λ 2 1 λ 3 3 + λ 2 1 λ 3 2 + 3λ 2 1 λ 3 − 1 = 0, corresponding to the surface instability of a purely elastic half-space (Flavin, 1963).
In the thin-plate or long-wave limit (kH → 0), we find that the antisymmetric bifurcation equation ( 64) can be rearranged as For B 2 = 0, Eq. ( 66) recovers the critical buckling condition λ 2 1 λ 3 = 1 for the purely elastic case, which results in λ cr 1 = 1.0 for both the uni-axial and plane-strain loading.This indicates that in the absence of magnetic field, the infinite-long plate or thin-plate buckles immediately when subject to an extremely small in-plane compression.

Specialization to the neo-Hookean magneto-elastic solid
For definiteness, we now specialize the previous results to the neo-Hookean (β = 0) ideal magneto-elastic model and the neo-Hookean magnetization saturation Langevin model, which are characterized by the following energy functions, respectively, and where m s is the saturation magnetization and χ is the magnetic susceptibility that is associated with the relative magnetic permeability µ r by µ r = 1/ (1 − χ).
Note that the magneto-elastic material models (67) and ( 68) take no account of particleparticle interactions and thus neglect magneto-mechanical coupling in terms of pure material magnetostriction.However, it is emphasized that the neo-Hookean magnetization saturation model ( 68) allows for a satisfactory quantitative and very good qualitative agreement with the experimental data presented in previous works (Psarra et al., 2017(Psarra et al., , 2019;;Goshkoderia et al., 2020), although it is anticipated to be less accurate in the post-bifurcation regime, especially when wrinkles develop substantially due to the large shear strains.As a result, more elaborate magneto-mechanical models such as the ones proposed recently by Mukherjee et al. ( 2020) are required to explore the effects of the strain-stiffening and the constituent phase properties (such as particle volume fraction, particle-particle interactions, etc) on the wrinkling instability of SMA structures, but such models are out of the scope of this paper.
Using Eqs. ( 2), (8) 2 and (67) or (68), we get the governing equations of the magnetization response for the ideal model, and for the saturation Langevin model.In the limit of small magnetic field (B → 0), one can verify that the saturation magnetization response ( 70) is compatible with the linear magnetization response (69).

Now introduce the following dimensionless quantities in terms of the initial shear modulus
G and the vacuum magnetic permeability µ 0 : Inserting Eqs. ( 67) and ( 68) into Eq.( 39) gives for the ideal model, and for the saturation Langevin model.
Thus, the magnetization responses ( 69) and ( 70) to the applied transverse magnetic field are written in non-dimensional form, as where for uni-axial loading, and for plane-strain loading.Practically, we determine the response by prescribing the prestretch λ 1 or λ and solving (75) for s 1 and λ 3 (uni-axial loading) or solving (76) for s 1 and s 3 (plane-strain loading).
For the neo-Hookean saturation Langevin model, the bifurcation or buckling equations ( 64)-( 66) are the same except that β = 0 and F 5 and F 55 are given in Eq. ( 73).

Results and discussion
We first conduct numerical calculations in Sec.5.1 to investigate quantitatively the nonlinear static response of incompressible SMA plates subject to mechanical and magnetic loads.Section 5.2 then focuses on the bifurcation analysis to calculate the critical mechanical and/or magnetic field generating the wrinkling instability.
We consider two different loading modes (plane-strain and uni-axial loading) and two neo-Hookean magneto-elastic models (( 67) and ( 68)).The material properties used in the numerical computations are taken as G = 10 kPa, χ = 0.4 and µ 0 m s = 0.5 T, which are obtained from experiments with a class of magnetorheological elastomers (Psarra et al., 2017(Psarra et al., , 2019) ) consisting of a soft silicone mixed with iron particles at a volume fraction of 20%.More details about the fabrication technique can be found in the paper by Psarra et al. (2017).

Nonlinear static response
The nonlinear static response of an SMA plate subject to mechanical and magnetic loads is calculated from Eqs. ( 74)-( 76).For plane-strain 2(a) shows the magnetization response of the dimensionless magnetization M 2 = M 2 /m s versus the dimensionless magnetic induction field B 2 = B 2 /m s (see Eq. ( 74)) for the two material models ( 67) and ( 68).We see from Fig. 2(a) that a neo-Hookean ideal magneto-elastic plate exhibits a linear magnetization response, whereas the response of a plate with magnetization saturation is nonlinear.The magnetization responses of the two material models are essentially identical for B 2 ≤ 1.0 (i.e., B 2 ≤ 4.46).The magnetization begins to saturate at B 2 6.0 (i.e., B 2 26.76).We note from Eqs. (73) 1 and (74) that the magnetization response of each material model is independent of the mechanical stretch ratio.
For plane-strain loading, the effect of the magnetic induction field B 2 on the nominal mechanical traction s 1 applied to the SMA plate is plotted in Fig. 2(b) for three different values of pre-stretch λ = 0.75, 1.0, 2.0.Clearly, when B 2 increases, the required mechanical traction increases monotonically, which indicates that the SMA plate has an in-plane con-traction trend because of the increasing external Maxwell stress.For the three pre-stretches λ = 0.75, 1.0, 2.0, the mechanical tractions corresponding to the two material models overlap up to B 2 3.8, 4.2, 5.0, respectively.However, the difference predicted by the two material models becomes more evident with subsequent increases in B 2 .Specifically, at the same level of B 2 , the plate with saturation magnetization effect requires a smaller mechanical traction as compared to the ideal magneto-elastic plate.This is because the presence of saturation magnetization will induce a smaller in-plane contraction trend.For uni-axial loading, the magnetization responses of the two material models are essentially the same as those for plane-strain loading, as shown in Fig. 2(a).For three different values of pre-stretch λ 1 = 0.75, 1.0, 2.0, Fig. 3(a) and 3(b) display the effect of the magnetic induction field B 2 on the stretch ratio λ 3 and the nominal mechanical traction s 1 , respectively, for the two material models.Again, we point out that the induced mechanical and magnetic field distributions are assumed to be uniform when solving the nonlinear static response, as explained previously, and that the material models ( 67) and ( 68) neglect pure material magnetostriction.We observe from Fig. 3 Furthermore, in the whole B 2 range of interest, the mechanical tractions based on the two material models are almost identical, because the saturation magnetization does not alter the contraction trend in the x 1 direction and just affects the compression amount in the x 3 direction.

Bifurcation analysis
We now examine the critical values of stretch in the x 1 direction and of transverse magnetic induction field for which antisymmetric (see Fig. 1(d)) and symmetric (see Fig. 1(e)) modes of wrinkling instability appear.
For plane-strain loading, antisymmetric modes are identified by bifurcation equations ( 64) and ( 77) with λ 1 = λ, λ 3 = 1 for the two material models.The critical buckling fields of antisymmetric modes for uni-axial loading are calculated by solving bifurcation equations ( 64) and ( 77) together with nonlinear static response (75) 2 .The corresponding critical fields of symmetric modes are calculated by using coth to replace tanh in Eqs. ( 64) and ( 77).The wrinkling criteria for thin-and thick-plate limits are obtained by evaluating Eqs. ( 65) and (66) or Eqs. ( 78) and ( 79) for the two material models.For plane-strain loading and uni-axial loading, Fig. 4(a) and 4(b) illustrate the critical combinations of stretch ratio λ or λ 1 and dimensionless magnetic induction field B 2 of antisymmetric wrinkling modes for the two material models.Specifically, Fig. 4 shows the results corresponding to the thin-and thick-plate limits and a representative value of kH = 1.
Since it is found that antisymmetric modes always occur first, Fig. 4 does not display the results for symmetric modes, which are addressed below.
We first focus on the results for the ideal magneto-elastic model in Fig. 4. In the absence of magnetic field (B 2 = 0, purely elastic case), a thin plate with kH → 0 is unstable for λ < 1 (λ 1 < 1) and is stable for any λ > 1 (λ 1 > 1) for plane-strain (uni-axial) loading.For both loading modes, a larger value of the parameter kH requires an increasing compression to induce instability.In the thick-plate limit kH → ∞, we recover the well-known critical compression stretches for surface instability in the purely elastic case, namely λ cr = 0.544 and λ cr 1 = 0.444 for plane-strain and uni-axial loading, respectively (Beatty and Pan, 1998).
For a fixed non-zero B 2 , the variation trends of the critical stretch λ cr or λ cr 1 with increasing kH are qualitatively the same as that for B 2 = 0. Besides, the critical magnetic field B cr 2 , for a given stretch λ or λ 1 , increases monotonically with kH, resulting in enhanced stability.
Furthermore, for a given kH, Fig. 4(a) shows that the critical stretch λ cr exhibits a monotonous increase when B 2 goes up, which means that the plate become more and more unstable and the magnetic field has a destabilizing effect.In particular, thin plates with kH → 0 are unstable in tension (λ cr > 1) for non-zero B 2 , while the plate with kH = 1 has a wrinkling instability in tension for B 2 2.8.Similar phenomena are observed in Fig. 4(b) for uni-axial loading.
We now evaluate the effect of saturation magnetization on the stability.Fig. 4(a) shows that for plane-strain loading, the critical stretch λ cr predicted by the saturation Langevin model coincides with that based on the ideal model for small and moderate values of B 2 , the range of which depends on kH.For example, the results based on the two material models are almost the same when B 2 3.0, 3.5, 4.0 for kH = 0, 1, ∞, respectively.However, the presence of saturation magnetization reduces remarkably the critical cr for a large value of B .These also found in Fig. 4(b) for uni-axial loading.Nevertheless, compared with plane-strain loading, the effect of saturation magnetization on the critical fields is weaker for uni-axial loading because there is no constraint in the x 3 direction.

Critical stretch for a prescribed magnetic load
For a prescribed magnetic induction field, we first determine the critical stretch of the underlying deformed configuration for which antisymmetric and symmetric wrinkling modes are induced.Specifically, for plane-strain (uni-axial) loading, Fig. 5(a) (Fig. 6(a)) displays the variation of the critical stretch λ cr (λ cr 1 ) with kH for the neo-Hookean magnetization saturation SMA plates subject to three representative values of B 2 = 0, 2.5, 5.0 (B 2 = 0, 2.5, 4.0), wherein the antisymmetric and symmetric solutions are represented by the solid and dashed-dotted lines, respectively.We find from Figs. 5(a) and 6(a) that antisymmetric modes always occur before symmetric modes become possible for any value of kH.Therefore, to realize a symmetric buckling mode we must in principle suppress the appearance of the antisymmetric mode.
For a fixed B 2 , the critical stretch λ cr (or λ cr 1 ) required to initiate the antisymmetric instability decreases monotonically with increasing kH, asymptotically approaching the surface instability of the thick-plate limit when kH → ∞.However, the symmetric bifurcation curves exhibit an opposite trend.Thus, the stable range of combinations of λ (or λ 1 ) and kH is determined by the region above the solid line for a fixed B 2 .Moreover, the critical stretch λ cr (or λ cr 1 ) for any value of kH is shifted upwards when raising B 2 , thus indicating that the SMA plate is destabilized by the application of an increasing magnetic field.
In particular, for plane-strain loading with B 2 = 0, 2.5, 5.0, the critical stretches λ cr of a plate with kH → 0 are 1.000, 1.426, 2.282, while those with kH → ∞ are 0.544, 0.608, 0.961, respectively.For uni-axial loading with B 2 = 0, 2.5, 4.0, the critical stretches λ cr 1 of the thin-plate limit are equal to 1.000, 2.016, 3.080, while those of the thick-plate limit are 0.444, 0.639, 0.918, respectively.The effect of small magnetic field on the stability is more significant for uni-axial loading than for plane-strain loading.
Figure 6: Uni-axial loading: (a) critical stretch λ cr 1 as a function of kH for anti-symmetric (solid lines) and symmetric (dashed-dotted lines) modes of the neo-Hookean saturation Langevin plates subject to three prescribed values of B2 = 0.0, 2.5, 4.0; (b) λ cr a function of kH for anti-symmetric modes of the neo-Hookean ideal (dashed lines) and saturation Langevin (solid lines) magneto-elastic plates subject to four fixed values of B2 = 0.0, 2.5, 4.0, 5.0.
For antisymmetric modes, Figs.5(b) and 6(b) illustrate how the saturation magnetization affects the bifurcation curves (λ cr versus kH and λ cr 1 versus kH) for plane-strain and uni-axial loading, respectively.We see that the bifurcation curves based on the two material overlap in the entire range of kH for B 2 ≤ 2.5.This means that the critical stretch is hardly affected by the saturation magnetization for small to moderate values of the magnetic field.
However, the saturation magnetization plays an important role in determining the critical stretch λ cr or λ cr 1 for large values of B 2 : the figures show that it then stabilizes the SMA plate, and that its effect is much stronger for plane-strain loading than for uni-axial loading.

Critical magnetic induction field for a fixed pre-stretch
We now evaluate the effect of pre-stretch on the bifurcation curves (B  symmetric modes, both asymptotically tending to the surface instability of the thick-plate limit.The stable region of B 2 and kH for a fixed pre-stretch is enclosed below the solid line for antisymmetric modes.Furthermore, an increase in the pre-stretch results in a larger value of B cr 2 for a given kH.This means that increasing the pre-stretch enhances the stability of SMA plates.Besides, we see from Figs. 7(a) and 8(a) that in the absence of pre-stretch (i.e., λ = 1 or λ 1 = 1), a thin plate (kH → 0) buckles immediately when subject to an extremely small transverse magnetic field.By contrast, for a thin plate with a pre-stretch (say λ = 1.25 or λ 1 = 1.25), a non-zero magnetic induction field B cr 2 2.0 (plane-strain loading) or B cr 2 1.25 (uni-axial loading) is required to trigger the buckling instability.Interestingly, if an SMA plate with small values of kH is subject to a pre-compression (i.e., the pre-stretch is less than 1), the underlying configuration is unstable for any applied B 2 .For example, a plate with λ = 0.8 or λ 1 = 0.8 is unstable for kH 0.97 under plane-strain loading and for kH 0.85 under uni-axial loading.Physically, this means that plates with small values of kH do not support pre-compression even if there is no magnetic field.
For antisymmetric modes, the effect of saturation magnetization on the bifurcation curves (B cr 2 versus kH) is highlighted in Figs.7(b) and 8(b) for plane-strain and uni-axial loadings, respectively.We observe that for a given pre-stretch, the two material models predict an identical critical magnetic field B cr 2 for small and moderate values of kH.For example, the bifurcation curves of a plate without pre-stretch (λ = 1 or λ 1 = 1) overlap for kH 1.0 and kH 1.4 under plane-strain and uni-axial loading, respectively.For a large value of kH with a fixed pre-stretch, the saturation Langevin model leads to a higher B cr 2 compared with the prediction of the ideal magneto-elastic model.On the other hand, for a large value of kH, the predicted difference based on the two material models become larger and larger when increasing the pre-stretch.

Euler's buckling approximations
For the magneto-elastic coupling case, it is useful to establish thin-plate approximate equations (i.e., the Euler buckling solutions) of the antisymmetric wrinkling modes, as they always occurs first.The derivation procedure is provided in Appendix C in detail.We specialize the analysis to the neo-Hookean ideal magneto-elastic model ( 67) because its predicted bifurcation curves coincide with those based on the magnetization saturation model for small and moderate values of kH, as shown in Figs.5-8.
For plane-strain loading, we find that the critical stretch is approximated as where of a slender plate under plane-strain loading (Beatty and Pan, 1998).In Appendix C we also establish the quadratic expansions of the critical magnetic field and the corresponding expansions for the uni-axial loading mode.We see from Fig. 9 that in the absence of magnetic field (B 2 = 0), the λ cr − kH curve and the λ cr 1 − kH curve for the thin plate should be approximated quadratically, as in the purely elastic case (Beatty and Pan, 1998).For B 2 = 0, the earliest correction for the stretch is of the first order in kH.For plane-strain loading, the first-order Euler solutions provide enough accuracy to approximate the exact solutions without having to resort to the second-order correction.However, for uni-axial loading, the linear approximations are not great and the quadratic corrections are required to approximate the exact bifurcation curves.
Fig. 10 shows that in the absence of pre-stretch (λ = 1 or λ 1 = 1), the first-order Euler buckling solutions can approximate well the B cr 2 − kH curve for the thin plate under the two loading modes, as described by Kankanala and Triantafyllidis (2008).But for a pre-stretch larger than 1, the B cr 2 − kH curve should be approximated quadratically for small values of kH to enlarge the effective range of Euler's buckling solutions.

Conclusions
In the framework of nonlinear magneto-elasticity theory and its associated incremental formulation, we presented a comprehensive theoretical analysis of the wrinkling instability of SMA plates under the combined action of transverse magnetic field and in-plane mechanical loading.We discussed two loading modes (plane-strain loading and uni-axial loading) and two types of neo-Hookean magneto-elastic material models (ideal model and magnetization saturation model).Employing the Stroh formulation and the surface impedance method, we derived explicit bifurcation equations of symmetric and antisymmetric modes and obtained their corresponding thin-and thick-plate limits analytically.Finally, we conducted detailed calculations to demonstrate the dependence of the nonlinear static response and bifurcation diagrams on the loading mode, load amplitude, and saturation magnetization.Our main observations are summarized below: (1) In contrast to the ideal model, introducting saturation magnetization results in a nonlinear magnetization response.Its effect on nonlinear mechanical response and critical buckling fields is more significant for plane-strain loading than for uni-axial loading.
(2) Antisymmetric wrinkling modes always appear before symmetric modes are expressed.
(3) Increasing the pre-compression and the magnetic field weakens the stability of SMA plates.However, the saturation magnetization effect strengthens their stability, especially for large pre-stretch or high magnetic field.
(4) The thin-plate approximate formulas agree well with the exact bifurcation curves for thin plates.
Note that we made the assumption of uniform magneto-mechanical biasing fields in this work to simplify the mathematical modelling of infinite SMA plates and to derive explicit analytical solutions.One factor that influences the uniformity of biasing fields is the shape and geometrical size of the SMA specimen.The mechanical and magnetic field distributions are usually non-uniform for ellipsoidal and cylindrical SMA specimens (see, for example, Martin et al. (2006); Rambausek and Keip (2018); Lefèvre et al. (2017Lefèvre et al. ( , 2020))).In addition, the SMA plate slenderness may affect the magneto-mechanical field distributions.If the biasing fields are not uniform (i.e., they vary with spatial coordinates), then the Stroh formulation for wrinkling instabilities becomes a set of first-order differential equations with variable coefficients, which are difficult to solve analytically.Nonetheless, as shown by the work of Shuvalov et al. (2004) for example, the Stroh formalism is still useful and forms the basis of a robust numerical resolution with the surface impedance matrix method.Nonuniform biasing fields are beyond the scope of this preliminary paper and are worthy of further research.

Appendix B. Bifurcation equation of surface instability based on the surface impedance method
Assume that the SMA half-space in the reference and current configurations occupies the region X 2 ≥ 0 and x 2 ≥ 0, respectively.To satisfy the decay condition at x 2 → +∞ in the half-space, we only keep the eigenvalues in Eq. ( 36) with positive imaginary parts.Thus, we take the three eigenvalues q 1 , q 2 , q 3 according to Eqs. ( 41) and ( 42).The general solution to the Stroh formulation (33) is written as A j η (j) e iq j kx 2 = 3 j=1 A j η (j) e −p j kx 2 , (B.1) where q j = ip j (j = 1, 2, 3) with p j > 0. In matrix form, Eq. (B.1) is expressed as where where Z = −iP 2 P −1 1 is the surface impedance matrix of the half-space, through which the generalized traction and displacement vectors at the face x 2 = 0 are connected.
In view of Eqs. ( 55) and ( 56), the surface impedance matrix exterior to the material is (C.8)

Figure 1 :
Figure 1: (a)-(e): Schematic diagram of a rectangular SMA plate with Cartesian coordinates and geometry.(a): Undeformed configuration before activation.(b): Deformed configuration subject to plane-strain loading (λ3 = 1) with a transverse magnetic field B2 and (c): Deformed configuration subject to uni-axial loading with a transverse magnetic field B2 (the magnetic poles are assumed to be very close to the plate, to ensure quasi-uniformity of the mechanical and magnetic fields).(d)-(e): Antisymmetric and symmetric modes of wrinkling instability.(f): Antisymmetric wrinkling due to the application of an external transverse magnetic field to a pre-compressed SMA plate glued onto an inert substrate (taken from Psarra et al. (2017)).

Figure 2 :
Figure 2: Plane-strain loading: (a) magnetization response of the dimensionless magnetization M 2 = M 2/m s versus the dimensionless magnetic induction field B 2 = B2/m s from Eq. (74); (b) response of the dimensionless nominal mechanical traction s1 versus the dimensionless magnetic induction field B2 from Eq. (76)1 for three different values of pre-stretch λ = 0.75, 1.0, 2.0.Solid lines correspond to the ideal magneto-elastic model while dashed lines represent the saturation Langevin model.

Figure 3 :
Figure 3: Uni-axial loading for three different values of pre-stretch λ1 = 0.75, 1.0, 2.0: (a) response of the stretch ratio λ3 versus the dimensionless magnetic induction field B2 from Eq. (75)2; (b) response of the dimensionless nominal mechanical traction s1 versus B2 from Eq. (75)1.Solid lines correspond to model while dashed lines represent the saturation Langevin model.
(a) that the stretch ratio λ 3 decreases notably with increasing B 2 only due to the magnetic traction induced by the external Maxwell stress.The stretch λ 3 for the ideal model is slightly lower than that predicted by the saturation Langevin model, because the saturation magnetization generates a smaller stretch λ 2 in the thickness direction.It is clear from Fig.3(b) that the mechanical traction s 1 increases gradually with B 2 .After B 2 reaches a certain value, s 1 remains unchanged.This is because the in-plane elongation trend due to the compression in the x 3 direction counteracts the in-plane contraction tendency due to the external Maxwell stress in the x 1 direction.

Figure 4 :
Figure 4: Bifurcation curves (or critical combinations of stretch ratio and dimensionless magnetic induction field) of antisymmetric wrinkling modes for the thin-and thick-plate limits and kH = 1.0:(a) plane-strain loading (λ versus B2); (b) uni-axial loading (λ1 versus B2).Solid to the model while dashed represent the saturation Langevin model.

Figure 5 :
Figure5: Plane-strain loading: (a) critical stretch λ cr as a function of kH for antisymmetric (solid lines) and symmetric (dashed-dotted lines) modes of the neo-Hookean saturation Langevin plates subject to three prescribed values of B2 = 0.0, 2.5, 5.0; (b) λ cr as a function of kH for antisymmetric modes of the neo-Hookean ideal (dashed lines) and saturation Langevin (solid lines) magneto-elastic plates subject to four fixed values of B2 = 0.0, 2.5, 4.0, 5.0.

cr2
versus kH) of antisymmetric and symmetric wrinkling modes.For plane-strain (uni-axial) loading, Fig. 7(a) (Fig. 8(a)) shows the critical magnetic induction field B cr 2 as a function of kH for the neo-Hookean magnetization saturation SMA plates for four different levels of pre-stretch λ or λ 1 = 0.8, 0.9, 1.0, 1.25.Note that the solid and dashed-dotted lines denote the antisymmetric and symmetric solutions, respectively.We observe from Figs. 7(a) and 8(a) that antisymmetric wrinkling modes are always triggered before the symmetric modes in the entire range of kH, which is analogous to what Figs.5(a) and 6(a) show.As kH increases, the critical magnetic field B cr 2 for a given pre-stretch increases gradually for antisymmetric modes and decreases monotonically for

Figure 7 :
Figure 7: Plane-strain loading for four fixed values of pre-stretch λ = 0.8, 0.9, 1.0, 1.25: (a) critical magnetic induction field B cr 2 as a function of kH for antisymmetric (solid lines) and symmetric (dashed-dotted lines) modes of the neo-Hookean saturation Langevin plates; (b) B cr 2 as a function of kH for antisymmetric modes of the neo-Hookean ideal (dashed lines) and saturation Langevin (solid lines) magneto-elastic plates.

Figure 8 :
Figure 8: Uni-axial loading for four fixed values of pre-stretch λ1 = 0.8, 0.9, 1.0, 1.25: (a) critical magnetic induction field B cr 2 as a function of kH for antisymmetric (solid lines) and symmetric (dashed-dotted lines) modes of the neo-Hookean saturation Langevin plates; (b) B cr 2 as a function of kH for antisymmetric modes of the neo-Hookean ideal (dashed lines) and saturation Langevin (solid lines) magneto-elastic plates.

Fig. 9
Fig. 9 compares the critical stretch λ cr or λ cr 1 versus kH based on the exact solutions to that calculated by the thin-plate buckling approximations.Fig. 10 illustrates the bifurcation curves of the critical magnetic induction field B cr 2 versus kH obtained by the exact solutions and the Euler buckling solutions.The results for plane-strain loading are shown in Figs.9(a) and 10(a) while those for uni-axial loading are depicted in Figs.9(b) and 10(b).The solid, dashed, and dashed-dotted lines represent, respectively, the exact solutions, the first-order and second-order Euler solutions.

Figure 10 :
Figure 10: Critical magnetic induction field B cr 2 as a function of kH for antisymmetric modes of a neo-Hookean ideal magneto-elastic plate subject to different fixed pre-stretch: (a) plane-strain loading with λ = 1.0, 1.1, 1.25, 1.4; (b) uni-axial loading with λ1 = 1.0, 1.25, 1.5.The solid lines, dashed lines, and dashed-dotted lines represent, respectively, the exact solutions, the first-order and second-order Euler buckling solutions.