Numerical analysis of the magnetic shape memory effect based on the absolute nodal coordinate formulation

In this paper, the absolute nodal coordinate formulation (ANCF) is applied to simulate the magnetic shape memory effect. Using the absolute nodal coordinate formulation makes it possible to describe complicated or large deformation cases. The nonlinear bidirectional coupling terms between the mechanical and magnetic field are taken into account in the analysis of the single-crystalline Ni-Mn-Ga sample. A two-loop iteration procedure with variable steps is implemented to predict the magnetic-field-induced strain (MFIS) in the specimen under a changing external magnetic field and a constant auxiliary compression. In addition, the proposed approach is used to track the superelastic behavior of the magnetic shape memory alloy when subjected to a constant magnetic field. The approach effectively describes the hysteresis and superelastic phenomenon of the shape memory effect. The solution is compared here with solutions obtained using classical linear and quadratic quadrilateral elements. The deviation observed in the solution is discussed, and its cause is further clarified from a two-domain pure magnetostatic analysis of a permanent magnet. It is found that the accurate solution of such problems is associated with discontinuity of the normal component of the magnetic potential gradient across the domain interface. Special measures must be taken to make the absolute nodal coordinate formulation element compatible with the discontinuity. A mixed FEs strategy, which adopts ANCF FE in the displacement field solver and classical FE in the magnetic field solver, is proposed as an alternative option to rectify the problem, which is verified by predicting the MFIS and the dynamic mechanical response of a sample under cyclic compression.


Introduction
The magnetic shape memory effect (MSME) refers to the change in the physical shape of a magnetic material in the presence of a moderate magnetic field [28]. A magnetic shape memory alloy (MSMA) is a promising new smart material with certain prominent characteristics that include large magnetic-field-induced strain (MFIS) [23,25], high fatigue life, and ultra-fast actuation response [24]. Different aspects of the MSMA material have already been studied including transformation temperature [9], saturation magnetization, magnetic anisotropy energy [26], and the structure of magnetic domains [15]. Early experimental studies help reveal the underlying mechanisms behind the MSME phenomenon based on which different modeling methods have been proposed to predict the behavior of the MSMA. Existing modeling methods comprise the constrained theory [7], the statistical approach [9], the phase-field method [17], and others. Many models are constructed based on thermodynamics [11]. Among these, the models that get most attention use the variational principle [14,30] or the finite element method (FEM) because of their straightforward implementation procedures. In the finite element framework, an iterative numerical method has been proposed to obtain the displacement and magnetization of the specimen [31]. In Wang and Steinmann's work, several internal state variables are introduced to represent the micro-status of the specimen. Another approach that does not use numerical iteration was inspired by the Karush-Kuhn-Tucker conditions. It treats the internal state variables as an additional set of degrees of freedom and integrates the corresponding evolution laws into the finite-element-based equilibrium equations [1,5].
The above existing numerical analyses of the magnetic shape memory effect focus mainly on small deformation cases and even ignore the influence of the elastic deformation on the demagnetization field [5,31]. However, the MSMAs are known for their ability to develop large MFIS (up to 12% in single-crystalline MSMA with seven-layered orthorhombic martensite structure) [25]. Therefore, studying the behaviors of MSMA undergoing large deformation and with complicated deformation modes is practically important. How the specimen behavior is influenced by these conditions cannot be ignored. Additionally, research works on the numerical dynamic analysis of single-crystalline Ni-Mn-Ga are still rare [29]. Some works ignore the inertial effect to study the dynamic response [10,16] and do not integrate the equation of motion into the solving procedure. In a very recent paper [8], Fan and Wang pioneered the work on numerical analysis of dynamic mechanical response of a single-crystalline Ni-Mn-Ga sample by solving the equation of motion with the classical FE. Compared with the classical FE, the absolute nodal coordinate formulation (ANCF) offers a promising approach capable of accounting for large deformation and large reference motions and brings its inherent advantages to the future complicated dynamic analysis of MSMA.
Using the components of the deformation gradient as the nodal coordinates, the absolute nodal coordinate formulation imposes no restrictions on the amount of rotation or deformation within the finite element. Due to the simplicity and its consistency with strain and stress measures from general continuum mechanics [4], the absolute nodal coordinate formulation has been widely used in various applications especially for problems with geometrical and material nonlinearities [19,20,34,35]. Additionally, the ANCF has been verified as a powerful tool for multibody dynamic analysis. Since the kinematic description of ANCF is based on the global coordinate system, the ANCF method can act as a straightforward bridge to integrating the MSME simulation into a general multibody dynamic analysis, which is necessary for future different application scenarios of possible slender MSMA actuator/sensor. The objective of the paper is to develop an approach to simulating the behaviors of the magnetic shape memory alloy subjected to large deformation. To achieve this, a general large-rotation and large-deformation finite element formulation, i.e., the absolute nodal coordinate formulation is employed. Considering the bidirectional coupling effects between the mechanical and magnetic fields, a numerical procedure is presented comprising the mechanical FE solver, the magnetic FE solver, and the internal state variable evolution solver. The approach can qualitatively describe the full generation process of the magnetic-field-induced strain and capture the special features of the magnetic shape memory effect including hysteresis and the superelastic effect. The deviation between the result given by ANCF FE and classical FEs is further analyzed using a pure magnetostatic case. It is found that the gradient discontinuity issue on the interface need be overcome for the ANCF FE in the two-domain MSME simulation to obtain a more accurate solution. To fix the problem, a mixed FEs strategy, which adopts ANCF FE to solve the displacement field and classical FE for the magnetic field, is proposed.
The remainder of the paper is organized as follows: In Sect. 2, the microstate description of the MSMA specimen is briefly introduced. Section 3 constructs the ANCF finite element. The corresponding equilibrium equations based on the variational principle and the equation of motion is established. Subsequently, a numerical analysis framework is established to estimate the behavior of the MSMA based on adopting the newly constructed ANCF element as the approach to solve the displacement and the magnetic field. Section 4 presents several typical examples to verify the method. Meanwhile, the limitations of the ANCF method in MSME analysis are investigated, and a modification of the approach that is free of those limitations is proposed. With the modified approach, the dynamic mechanical response of a sample under cyclic compression is simulated. And finally, Sect. 5 offers the authors' conclusions. The magnetic shape memory effect (MSME) is caused by martensite variant reorientation and domain transformation at room temperature. Ever since Ullakko et al. [28] revealed the MFIS in 1996, the single-crystalline Ni-Mn-Ga alloy has been extensively studied. In the single-crystalline Ni-Mn-Ga alloys, three martensite structures are observed: five-layered modulated (5M) approximately tetragonal martensite with c/a< 1 and maximum theoretical strain of about 6%, seven-layered modulated (7M) orthorhombic martensite with maximum theoretical strain of about 12%, and non-modulated (NM) tetragonal martensite with c/a>1 with maximum theoretical strain of about 20%. No considerable magnetic shape memory effect has been reported in the last structure [27,36]. Compared with a 7M structure specimen, the samples in the 5M structure exist in a border range [36] and demonstrate best performance for practical application [22]. Therefore, a bulk MSMA specimen with 5M structure was studied in this work.
To get an evident MSME process, a compression assist is usually applied. A typical experimental procedure or the stress-assisted magnetic-field-induced strain generation procedure is depicted in Fig. 1a. After "training" under a cyclic and large mechanical stress, a single-crystalline Ni-Mn-Ga specimen in the 5M (5-layer modulated) structure exhibits a single variant state with the magnetic easy axis, i.e., the short c-axis, aligned with the external loading [15]. As shown in Fig. 1a, under the external magnetic field H a and the assisted compression t A , variant one begins to transform to variant two with the c-axis generally aligning with the magnetic field direction. With increasing the external magnetic field strength, variant two grows at the expense of the stress-favored variant one with simultaneous domain wall movement until the specimen evolves to the single variant and single domain state. Macroscopically, the specimen extends as a result of the internal variant reorientation, the rotation of the magnetization vector, and/or domain wall movement. Consequently, the final maximum induced strain is around 1 − c/a ≈ 0.06 [15].
In this study, a simple assumption based on the experiment [15] and adopted by Wang and Steinman [30,31] is used to describe the generation of MFIS. The corresponding micro-configuration of a two-variant 5M Ni-Mn-Ga sample is shown in Fig. 1b. The variants form distributed bands that are separated by twin interfaces. In each variant region, two kinds of magnetic domains are separated by 180 • domain walls. In each magnetic domain, the magnetization vector aligns with either the positive or negative direction of the corresponding easy axis with saturation magnitude M s . The external magnetic field will induce the movement of the magnetic domain walls in variant one and the rotation of the magnetization vectors in variant two as shown in Fig. 1b. Correspondingly, internal state variables α and θ are introduced to represent the magnetic domain volume fraction in variant one and magnetization vector rotation in variant two, respectively. In summary, only three internal state parameters, i.e., the variant two volume fraction ξ 2 , the magnetization vector rotation θ , and the magnetic domain volume fraction α are capable of fully describing the internal state of the specimen as shown in Fig. 1b. Accordingly, the elastic moduli tensor, the effective magnetization vector, and the total strain are dependent on these three parameters. The elastic moduli tensor can be written as [31] where C 1 and C 2 are the elastic moduli tensor of variant one and variant two, respectively. The tensor can be written in Voigt notation form as where κ i , i = 1, 2...6 are the moduli constants. Driven by the the magnetic domain wall motion, the rotation of the magnetization vector, and variant reorientation, the effective magnetization vectorM is expressed using the corresponding three internal state parameters as where e 2 is the unit direction vector of the external magnetic field measured in the current configuration. With the transformation from variant one to variant two, the transformation strain ξ 2 ε tr appears. As a result, the total strain is written as where ε el is the Green-Lagrange strain tensor. It is written as follows: where F = r ,X is the deformation gradient tensor, r is the position vector in the current deformed configuration, X is the material coordinate in the initial reference configuration, and the subscript following a comma denotes partial differentiation with respect to that variable. In Eq. (4), the maximum transformation strain tensor ε tr can be written in Voigt notation form as follows: where ε 1 and ε 2 represent the maximum transformation strains along the principle axes. According to the microstructure and the lattice dimensions of variant one and two, the induced strain in z direction does not appear during the variant transformation as shown in Fig. 1a. Therefore, it is reasonable to apply the plane strain assumption to simplify the simulation. In the plane strain case, the elastic moduli tensor in Voigt notation C 1 and C 2 degenerates to The elastic and maximum transformation strain in Voigt notation degenerate to

ANCF element construction and numerical analysis method
This Section briefly introduces the variational principle approach for the MSME problem, which leads to the FEM solution of the displacement and magnetic field. Then, the planar ANCF FE employing both the displacement and magnetic field degrees of freedom is constructed. Finally, a modified numerical simulation framework that integrates the ANCF FE solver and the internal state parameter evolution solver [30,31] is presented.

Variational principle approach
The governing partial differential equations and the evolution laws of the internal variables can be derived by calculating the variation of the total energy functional. The functional formulation presented in [31] is briefly introduced. The total energy functional of the whole magneto-mechanical system is where ψ is the scalar magnetic potential, U el is the elastic potential energy, U an is the magneto-crystalline anisotropy energy, U ze is the Zeeman energy, U sta is the magnetostatic energy, U mix is the mixture energy, U ml is the potential energy of the external mechanical load, and U dis is the total energy dissipation. All the terms are expressed under the Lagrangian description with respect to the reference configuration. In Eq. (9), the elastic potential energy U el can be written as: where r is the specimen region, ρ r is the mass density in the reference configuration, and φ is the elastic energy density. The magneto-crystalline anisotropy energy U an is written as: where K u is the anisotropy constant. The Zeeman energy, or the external field energy U ze is: where μ 0 is the magnetic permeability of free space and H a is the external magnetic field. The magnetostatic energy of the demagnetization field U sta can be written as: where R 3 is the analyzed whole space including the specimen domain r and the surrounding space domain r , H d = −F −T ∇ψ is the demagnetization field induced by the magnetization of the MSMA specimen measured in the current configuration, and J=det(F), i.e., the determinant of the deformation gradient tensor. According to the divergence theory and the Gauss's law for magnetism, the magnetostatic energy can also be written in the following form [30]: The mixture energy due to the interactions of different martensite variants and magnetic domains U mix is written as follows: where f ξ (ξ 2 ) and f α (α) are the energy densities arising from the mixture of the two variants and the interactions of the different magnetic domains in each variant region. The potential energy of mechanical loading U ml is: where t A is the external compression applied per unit referential area and ∂ r is the boundary, on which the compression is applied. Besides the above energy quantities, to capture the hysteresis phenomena observed in the experiments, the dissipation effect during the variant transformation process should be considered. The total energy dissipation U dis during the forward and reverse transformation processes can be written in the following form: where ξ 0 2 and ξ 2 are the initial and final volume fractions of variant two in the specimen. Accordingly, ξ 0 2 ≤ ξ 2 applies for the forward transformation and ξ 0 2 ≥ ξ 2 for the reverse transformation. D ± are positive functions called the dissipative resistances ("+" and "−" represent the forward and reverse transformation, respectively). By adding up the energy quantities together, the total energy functional for the magneto-mechanical system is: In the energy functional, the independent variables are the position vector r, the scalar magnetic potential ψ, and the three internal state variables ξ 2 , θ , and α. Specifically, in the FEM implementation, r and ψ are interpolated using the nodal vectors e u , e ψ of the corresponding element, and the internal state variables ξ 2 , θ , and α are directly calculated on the Gauss integration points of each element. The governing equations of the system are obtained by calculating the variations of G with respect to the independent variables r, ψ, and the internal state variables [30] as shown in the following: where H tot = H a + H dy , H a = H a · e 2 , and H dy = H d · e 2 are the components of the external and the demagnetization field in the external magnetic field direction. The inequalities as shown in Eqs. (21)(22)(23) are obtained based on thermodynamics [13]. Magnetostatic energy in the form of Eq. (14) is used in the variation of the total energy functional with respect to r. The governing equations indicate that the magnetic field and the displacement field are implicitly related via the internal state evolution of the MSMA. In addition, the influence of the deformation on the magnetic field explicitly goes into H d = −F −T ∇ψ via the deformation gradient tensor F, and the magnetic field influences the deformation via the variation of the magnetostatic energy as shown in Eq. (14). These two explicit coupling terms are ignored in [31] but considered in this work. The variational principle is a very general approach in formulating the constitutive models of the multiphysics problems such as the piezoelectric effect. Under the framework of variational principle, similar equations can be obtained in the piezoelectric effect [2,18]. However, there are essential differences between MSME and piezoelectric effect analysis. As shown in Sect. 2, the MSME involves with changes to the microcrystallographic structure. Correspondingly, several internal variables are introduced to describe the evolution of the micro-structure in the numerical simulation. The state of the MSMA specimen should be determined iteratively via equations obtained from Eqs. (21)(22)(23) since the internal variables depend on the current magnetic and displacement field. Thus, there must be an iterative state evolution solver in the MSME analysis procedure. As a result, most MSME simulations in the literature are quasi-static analyses [1,3,5,33]. By contrast, the numerical analysis of the piezoelectric effect is not associated with the micro-crystallographic structure. So there is no additional iterative solver besides the displacement and electric field solver in the piezoelectric effect analysis procedure [18].

ANCF element construction
As mentioned in Sect. 2, the classical MSME experiment as shown in Fig. 1a can be analyzed as a planar issue since there is no transformation strain along the z axis. For such planar rectangular domain, the planar quadrilateral element is the most appropriate element for the analysis [1,5]. Therefore, the planar quadrilateral element proposed by Olshevskiy et al. [21] is used in the analysis. For the magnetic field analysis, the scalar magnetic potential ψ is adopted as the basic variable, and the corresponding gradients ψ ,x , ψ ,y are also used as the nodal DOFs under the framework of the ANCF method, where x and y are the element coordinates in the straight configuration. By doing so, the ANCF method is checked in the MSME analysis by directly generalizing the original mechanical-element to a magneto-mechanical element. The schematic diagram of the element is shown in Fig. 2, As shown in Fig. 2, 12 variables, i.e., e 2,3,4, are used for the interpolation of the magnetic potential. Selecting the bases for the 12-term interpolation based on the Pascal triangle, a natural option is the same basis as that used in the displacement field interpolation, i.e., {1, x, y, x 2 , x y, y 2 , x 3 , x 2 y, x y 2 , y 3 }. Accordingly, the same shape functions s i , i = 1, 2, ..., 12 are obtained for the magnetic potential interpolation. The shape function matrix is shown as follows: where ξ = x/l and η = y/w are the local normalized coordinates in the straight configuration, and 0 ≤ ξ, η ≤ 1, and l and w are the element length and width in the straight configuration. The detailed expression of the shape functions s i , i = 1, 2, ..., 12 is shown in Appendix A. In the shape function matrices of the displacement field and magnetic field element, the identity matrix I ∈ R 3×3 and I ∈ R 1×1 in Eq. (24) where The deformation gradient with respect to the initial reference configuration is obtained as where x = [x, y] T are the element coordinates in the straight configuration, and the mapping matrix between the initial reference configuration and the straight configuration is denoted as J 0 . Denoting the initial nodal vector of an element as e u 0 , J 0 at the point with the local normalized coordinates ξ and η can be obtained as

Magnetic FEM equilibrium equation
The scalar potential ψ and its gradients ψ ,x , ψ ,y are used as the nodal degrees of freedom. The demagnetization field is H d = −F −T ∇ψ. Under the framework of FEM, the demagnetic field at an arbitrary point can be obtained as where S ,x and S ,y are the derivatives of the shape function matrix S with respect to the material coordinate. The dimension of the corresponding vectors and matrices is S, S ,x , S ,y ∈ R 1×n , and B ∈ R 2×n , where n is the number of the nodal coordinates in an element, and for the adopted ANCF element n=12. Substituting Eq. (24) into Eq. (19) results in the following equilibrium equation: where K ψ is the assembled system stiffness matrix consisting of two corresponding parts in the specimen region r , and the surrounding space region r , e ψ sys. is the assembled system nodal vector associated with the magnetic potential, and Q ψ r is the generalized force vector obtained in the specimen region. The corresponding expressions in an element domain elem are

Mechanical FEM equilibrium equation
Similar to the magnetic FE, the position field is interpolated as r = Se u . Recalling that ε = ε el + ξ 2 ε tr , the equilibrium equation can be written as where Q s sys. is the assembled system elastic force, Q tr sys. is the assembled system force associated with the martensite variant transformation, Q mag sys. is the assembled magneto-mechanical coupling force derived from the magnetostatic energy U sta as shown in Eq. (14), and Q ext sys. is the assembled system external generalized force. All the forces are defined in the specimen region, and the expressions of the internal forces on the element level are as follows: where the intermediate vectors S a , S b , S c , and S d are where (i, :), i = 1,2 means the i th row of the derivative of the shape function matrix in the reference configuration S ,X or S ,Y , which can be obtained by where J 0 (i, j) i, j = 1, 2 is the term at the i th row and j th column of the mapping matrix J 0 as shown in Eq. (27). The generalized nodal force of the external pressure on the element level is where n A is the unit outward normal vector of the end specimen surface, on which the compression with magnitude of t A is applied. Unlike the equilibrium equation of the magnetic FEM solver as shown in Eq. (29), the counterpart in the mechanical FEM solver Eq. (31) is nonlinear and solved via the Newton-Raphson iteration.

Equation of motion
To predict the dynamic mechanical response of a sample, the equation of motion should be established. Based on the principle of virtual work, the equation of motion on the element level can be written as where M elem is the element mass matrix, andë u is the second time derivative of the nodal vector e u associated with the position field. For the ANCF element, the element mass matrix is Clearly, the mass matrix remains constant in the integration loop of the equation of motion, which is calculated only once beforehand bringing computational merit. In this work, the equation of motion (37) is solved via an implicit integration generalized-α method [6].

Framework of numerical solution
Besides the displacement and magnetic field solution obtained from the FEM solver presented in Sect. 3.2, the internal state variables p = {ξ 2 , θ, α} need to be determined as the input variables for the FEM solver.
To some extent, the internal state variables act as the bridge that connects displacement and magnetic field. Therefore, the numerical solution procedure must integrate the sub-solver of the internal state variables p into the procedure. Considering the natural constraint −1 ≤ sin θ ≤ 1, the evolution law of sin θ can be determined by recalling Eq. (22), where F α = 1 + H a + H dy /H α cri /2, and H α cri is a critical magnetic field value. Equation (40) can be written based on Eq. (21) and the physical constraint 0 ≤ ξ 2 ≤ 1 as where ξ 2 * is the solution of F ± ξ 2 = ρ r (±D ± − π ξ 2 ) = 0. The dissipative resistances D ± (ξ 2 ) and intermediate function π ξ (r, ψ, ξ 2 , θ, α) are defined as where A i , B i , i = 1, 2 and B are constants in the dissipative resistances and mixture energy density. By substituting Eq. (41) into F ± ξ 2 = ρ r ±D ± − π ξ , it is possible to write F ± ξ 2 into a quadratic equation, where the coefficients are where C 21 = C 2 − C 1 , c ± 1 and c ± 2 are four constants associated with D ± , and c 3 (ψ, θ, α) is an intermediate variable independent of ξ 2 . The signs in front of c ± 1 and c ± 2 are "+" and "-", respectively, for the forward and reverse transformation process. Their expressions are Normally, there are two roots for the quadratic equation (40) evolution law. Here, due to the equality of ε 1 and ε 2 , the quadratic equation (42) degenerates to a linear equation by substituting Eqs. (7) and (8) into Eq. (43). As a result, ξ 2 can be uniquely determined by Eq. (40). Figure 3a shows the flow chart of the modified two-loop numerical simulation procedure based on [31,32]. The procedure comprises three solvers, i.e., the mechanical FE solver, the magnetic FE solver, and the evolution solver of the internal state variables. These solvers solve the corresponding partial differential equations to, respectively, obtain the displacement, the scalar magnetic potential, and the internal state variables of the specimen. Because the convergence speed and the computational cost for each numerical solver is different, the procedure is divided into two loops. Solving the displacement field is set as the outer loop, while the magnetic potential ψ and the internal state variables p = {ξ 2 , θ , α} are calculated in the internal loop. The iterative procedure keeps calculating and updating the state of the specimen until the residual is smaller than the specified tolerance. In the outer loop, the residual in the axial strain is adopted. To improve efficiency, a linear predictor is adopted to determine the initial value of each step based on the converged values obtained in the previous two steps. Specifically, the initial value of the variables in the next step is predicted as S (1) , where S is the collection of the unknown variables to be determined. The subscript "0" indicates the initial value, the subscript "1" indicates the value that has been updated after the solution in the corresponding solver, the superscripts "1", "0" and "-1" denote the next, current, and previous step, and λ is the ratio between the current and the previous step. Based on the static analysis, the dynamic analysis flow chart is shown in Fig. 3b. The compression is set as the time-dependent variable to introduce the algorithm. In each time step, the equation of motion is firstly solved, and the position DOFs r 1 , ψ (0) 1 , and p (0) 1 are passed to the next outer loop step if the strain difference between the current step and the last step in the inner loop is smaller than the tolerance.

Numerical examples
Numerical examples are presented in this Section to assess the performance of the proposed element with respect to the comparison with the classical FE. Firstly, two typical behaviors of the MSMA are simulated: the specimen under a changing magnetic field and the specimen undergoing various compressions. Accordingly, two key features of the MSME, i.e., the memory effect and the superelastic effect, are simulated. Then, a permanent magnet with the surrounding space is analyzed to further check the accuracy of the ANCF element in the pure magnetostatic analysis. Finally, the dynamic mechanical response of the MSMA specimen under a time-dependent compression is simulated. Figure 1a illustrates the typical experiment that examines magnetic-field-induced strain generation under auxiliary compression. This situation is simulated to verify the new element and the presented numerical procedure. The magnetic field problem is essentially unbounded. To simulate the problem under the framework of the FEM, a simple truncation technique is adopted with a ten-fold larger space surrounding the specimen. As shown in Fig 4a, the scalar magnetic potential on the boundary of the surrounding space is set to zero [31]. With regard to the mechanical constraints, nodes on the middle plane of the specimen are restricted along the x direction, and the position of the center node is fixed such that the rigid motion mode is removed to avoid a To induce the magnetic shape memory effect, auxiliary compression of magnitude t A is applied along the x axis, and an external magnetic field H a is applied along the positive y axis direction. Exploiting symmetry, only a quarter model is adopted. Correspondingly, the boundary conditions of the magnetic and displacement field change as shown in Fig. 4b. For the magnetic field, a zero scalar magnetic potential constraint is applied on the horizontal antisymmetric edge. The normal component of the potential gradient is set to zero on the symmetric edge nodes. For the displacement field, an additional constraint is applied at the y position on the horizontal antisymmetric edge nodes. The specimen is meshed using the FE with both e u and e ψ as nodal vectors. The surrounding space is modeled by the FE with only e ψ as nodal vectors. Internal state variables are directly calculated on the Gauss integration points of each element. Besides the constructed ANCF FE, classical FEs including the linear and 8-node quadrilateral FE are also adopted to simulate the MSME process.

Simulation of the magnetic shape memory effect
The single-crystalline Ni-Mn-Ga alloy in the 5M structure martensite state is analyzed, and the constitutive matrix is shown in Eq. (7). The material parameters adopted in the literature [31] are Other parameters associated with the magnetization, the anisotropy energy, the mixture energy, and the dissipative resistances are [31]  The meshing schematic diagram of the specimen and the surrounding space are shown in Fig. 5. The converged results obtained with n 1 = n 3 = 16 and n 2 = 8 are presented in the following Sections.

Magnetic-field-induced strain at a constant compression level
A specimen under a constant compression and various magnetic field strengths is analyzed quasi-statically by gradually increasing or decreasing the applied magnetic field. The external magnetic field H a = 1×10 6 e 2 A/m, and the initial step size is set as 2 × 10 4 A/m. Both the forward and the reverse process are simulated. Compression t A = 1.0 MPa is applied along the x axis on the specimen. Figure 6 shows the average magneticfield-induced strain and the effective magnetization magnitude. The numerical model effectively describes MSME hysteresis behavior. A full forward and reverse martensite variant reorientation process is obtained under the compression t A = 1.0 MPa. As shown in Fig. 6, different data points on the ANCF curve marked using arrows are numbered from 1 to 6 for the sake of description. Table 1 lists the internal state parameters and macro-variables in the corresponding configurations.
In the initial configuration 1, the specimen exhibits the single variant and domain state. Applying the external field, the specimen is magnetized via the rotation of the magnetization vector. At the beginning, the magnetization is not strong enough to initiate the variant reorientation. With the increase in the external field strength and the subsequent rise in magnetization, variant two appears in configuration 2. After configuration 2, as shown in Fig. 6b, the increase in the magnetization becomes faster due to both magnetization vector rotation and variant reorientation in this stage. Later, the variant two volume fraction reaches nearly 1 in configuration 3. Then, the rotation of the magnetization vector keeps increasing and nearly goes to 90 • in the final configuration 4. Consequently, the specimen is fully magnetized, and the induced strain nearly reaches its maximum value of 0.06. Then, at the beginning stage of the reverse process, e.g., in configuration 5, the magnetization vector in variant one rotates back towards the compression direction, but the induced strain remains mostly unchanged.  After the field decreases to a critical value, strain begins decreasing. In the next configuration 6, variant one recovers to 56% and the effective magnetization becomes 0.44M s due to the further decrease in the external field strength. With the removal of the external field, variant one almost totally recovers. As a result, the induced strain and the magnetization vanish. The different paths in the forward and reverse process clearly show the hysteresis response of the induced strain and magnetization.
To get a further interpretation on the evolution of the internal state of the specimen, the distribution of the magnitude of the effective magnetizationM and the variant two under two intermediate magnetic field strengths H a = 5.0 × 10 5 and 7.0 × 10 5 A/m is plotted in Figs. 7 and 8. As shown in Fig. 7, the distribution of magnetization is not uniform. Specifically, the maximum magnitude appears on the edges while the minimum magnetization is found in the center region. This is consistent with the illustration in [31]. Correspondingly, a nonuniform distribution of variant two is also observed in Fig. 8 before variant reorientation is complete. As a result, a non-homogeneous deformation mode appears regardless of the uniform uni-axial loading. In detail, little shrink appears on the vertical edge. This is because the elastic moduli tensor becomes different in various regions of the specimen due to the induced different internal state variables or the internal microstates in different local regions. However, this non-homogeneous phenomenon disappears with the increase in the magnetic field strength since the internal state variables reach their "saturation" value and are uniform throughout the specimen domain. This intermediate non-homogeneous phenomenon is not observed in the solution given by the linear FE due to its limited ability to describe non-uniform strain in an element.

The superelastic response of the MSMA subjected to a constant magnetic field
In this Section, the stress-strain response of the specimen under a constant external magnetic field is investigated. The superelastic effect under a constant magnetic field is another special feature of the magnetic shape memory alloy. The specimen is studied initially in the single variant two state. Under compression, variant two transforms to variant one. Consequently, the MSMA exhibits a superelastic response, where a small increase in compression induces a significant increase in strain. In contrast to the MFIS simulation in Sect. 4.1.1, the induced strain is negative here. Under the field of μ 0 H a = 1.0 T, compression is generally applied from 0 and 5 MPa in a quasi-static process. Figure 9 shows the stress-strain curve. The ANCF element effectively In summary, the introduced ANCF element is capable of simulating the special behavior of the MSMA and following the full process of the generation of MFIS. However, the ANCF element results do not fully agree with the results obtained from the two classical FEs as shown in Figs. 6 and 9. To explore the cause for this deviation, the new element is examined in a pure magnetostatic analysis in the following Section.

Magnetostatic analysis of a permanent magnet
A pure magnetostatic analysis as shown in Fig. 10 is performed using the ANCF element in this Section. The magnetic field inside a bulk permanent magnet domain r and the surrounding space r are analyzed. The As shown in Fig. 10a, the interface continuity conditions between two domains are where · indicates the jump between the variable value on the two interface sides, and n d is the normal vector of the interface directing from r to r . In other words, the normal component of B and the tangential component of H are continuous on the interface. In the framework of the variational principle, such interface continuity conditions enter into the FE formulation as the nodal loads on the interface r ∩ r as shown in Eq. (B.6) in the detailed derivation in Appendix B. Besides the interface continuity conditions, the magnetic insulation condition B · n = 0 is applied on the boundaries in this analysis, where n is the outward normal vector on the boundary ∂ r . As was the case for the analysis in Sect. 4.1, only a quarter model is analyzed as shown in Fig. 10b. All the conditions on the boundaries 3 ∼ 6 and the interface 1 ∼ 2 in the specific model are shown in Fig. 10b and Table 2.
Among the whole conditions, all the Neumann conditions automatically hold under the framework of the variational principle method or the finite element method as shown in Appendix B. The equality constraints of the scalar potential and the tangential gradient are automatically satisfied by assigning a unique value on the interface nodes using FEM as shown in Appendix C. The remaining continuity condition B · n = 0 leads to the nodal load on the interface Q = d (−M · n d )S T d d , where d is the interface. In summary, the only Fig. 11 Distribution of H y in the permanent magnet model boundary condition that needs to be explicitly applied is the Dirichlet condition on 3 . The model is analyzed using the classical linear quadrilateral element, the classical 8-node quadratic quadrilateral element, and the proposed ANCF magnetic quadrilateral element. The result provided by the commercial software COMSOL is used as the reference. In the element implementation, it is important to notice the difference between ANCF FE and the classical FE. Using magnetic potential gradients as additional nodal degrees of freedom, the ANCF element brings a disadvantage, i.e., extra explicit constraints on the nodal gradients should be applied to be consistent with the Neumann boundary conditions (BCs) on 4 ∼ 6 in the ANCF element implementation. The distribution of the component of the magnetic field H along the y axis is shown in Fig. 11. To get a clear look into the interface, the local region 0.02 × 0.01 m 2 is presented in Fig. 11. A unique range is employed in all the Figures for the sake of comparison. It is clear that the magnetic field calculated by the classical elements agrees well with that obtained by COMSOL, while the ANCF element cannot give a close result. This is because the ANCF element introduces gradients as nodal degrees of freedom, which implies that all the gradient components of ψ on the nodes are manually assigned to be equal on the two interface sides. These explicit equality constraints on the interface nodes largely change the local distribution of the normal gradient on the interface though the ANCF element cannot lead to a continuous normal gradient at the non-node side points (the proof is presented in Appendix C). However, based on the different constitutive relations as shown in Eq. (44) in two domains and the continuity conditions as shown in Eq. (45), in this specific analysis there is a jump in the y axis component of the magnetic field the H on the horizontal interface 1 . Looking into the value of H y on interface node A as shown in Fig. 10b, the classical linear and 8-node quadrilateral elements give a jump of 725 A/m and 750.3 A/m on the domain interface, respectively. These are close to the theoretical value 750 A/m. According to the above analysis, the jump in H y on interface 1 should be explicitly guaranteed on the nodes in the ANCF analysis. To handle the normal gradient gap issue, two equivalent options are available. For the first option, an extra set of nodes can be introduced on the interface. In the meantime, the constraints setting the magnetic potential and the tangential gradient of the two sets of nodes to be equal are applied. As a result, only the normal magnetic potential gradients of the second set of nodes are independent DOFs. Alternatively, rather than the extra set of nodes, another set of independent normal gradients on the interface nodes can be directly introduced. Either of them will bring a new set of independent gradients to the system. Due to the straightforward implementation, the first procedure is adopted in this paper to release the implicit full equality constraints on the magnetic potential gradients of the interface nodes in the original model. By doing so, the ANCF analysis is expected to be improved. Figure 12 shows the obtained distribution of H y . Apparently, the solution is improved, and the jump in H y on the interface is captured. However, the obtained gap value H y = 803 A/m on point A is still farther from the reference solution than values obtained using classical finite elements. This is because the jump constraint on H y is not exactly applied at the interface.
In summary, the cause for the deviation of the ANCF element solution is rooted in the direct adoption of magnetic potential gradients as FE nodal coordinates. This brings an added and even more complicated explicit constraint on the gradients to handle the natural BCs and interface gradient discontinuity. By contrast, in the case of the classical FE, such conditions are implicitly satisfied in the framework of variational principle.
In conclusion, both the classical linear and 8-node quadratic quadrilateral element can give a close result to that obtained by COMSOL. Nevertheless, the ANCF formulation is not recommended in the magnetic analysis with two or more different medium domains if no special operations are taken to handle the domain interface continuity conditions. The reasons are as follows. Firstly, an extra set of nodes or additional corresponding gradient DOFs need be introduced to capture the "jump" of the corresponding magnetic field component on the domain interface. Additionally, complicated constraints need to be explicitly applied on the interface nodes to guarantee the continuity conditions on the domain interface. To exactly describe the interface continuity condition, the gap value on the corresponding field component must be obtained in advance, which is virtually impossible in problems without analytical solutions. Therefore, the ANCF FE must be improved for the twodomain MSME analysis. As an alternative, using the ANCF FE for the deformation analysis but the classical FE for the magnetic field analysis can be a good option. From the perspective of interpolation, the displacement is described using the cubic 12-term interpolation while the magnetic potential is interpolated with lower order bases without gradients as interpolation conditions. Figure 13 shows the MSME simulation results obtained using mixed FEs, i.e., ANCF FE and linear FE.
Apparently, the mixed FEs strategy works well since it gives very close results to those obtained by the classical FE while keeping all potential benefits of the ANCF approach in the mechanical analysis involving large deformation. Additionally, the simulation time of the mixed FEs strategy is less than half (around 45%) that of the original ANCF simulation. The mixed FEs strategy provides an attractive way to integrate the ANCF method in the MSME analysis of the complicated specimen with large deformation. In this Section, the dynamic mechanical response of the MSMA specimen under a constant magnetic field and a high-frequency cyclic compression is simulated via the mixed FEs strategy. The same specimen as employed in Sect. 4.1 is adopted. Referring to [8,16], a linearly varying compression is used as the mechanical impact. As shown in Fig. 14a, the magnitude of the compression is 5 MPa, and the period is 0.02 s. The magnetic field is μ 0 H a = 1.0 T, which guarantees a full superelastic loop for the studied sample as shown in Fig. 9. As was the case in Sect. 4.1.2, the specimen is initially in the single variant two state. The generalized-α integrator with a numerical damping is used for the solution of the equation of motion. The mechanical response of the sample is shown in Fig. 14b.
As shown in Fig. 14b, the mechanical response of the MSMA sample consists of 6 stages. From A to B, the sample exhibits elastic axial contraction under the effect of compression t A . With the increase in compression, the transformation from variant 2 to variant 1 begins at point B, and the macro-axial strain jumps to the maximum induced strain 0.06 at point C. The increase in compression only leads to a tiny elastic growth in the axial strain. Corresponding to the peak value of 5 MPa at 0.01 s, the axial strain reaches the maximum at point D. Similarly, as compression drops, the sample demonstrates elastic behavior in segment DE and FG and decreases sharply from E to F as a result of the variant transformation. During the simulation, there is vibration on the strain at point C and F due to the inertial effect. This cannot be seen in the quasi-static analysis [10]. As was the case in [8], a switch-like response of the axial strain is obtained. The validity of the proposed approach is verified.
As a summary of the numerical examples Section, a clarification is made as follows. This work tries to check the possibility of applying the ANCF method to the MSME analysis not just to propose a new magnetomechanical element. As a first step, the ANCF method is implemented for both mechanical and magnetic field via a magneto-mechanical ANCF element. It is demonstrated that the approach can well capture the key behaviors of MSME. Meanwhile, the limitations of this approach are highlighted, which come from the magnetic part of the analysis. As a measure for overcoming those difficulties, a mixed FEs strategy is proposed, in which ANCF is only used in the mechanical part, and the magnetic equations are solved using classical FEs. This approach keeps inherent advantages of ANCF in the mechanical part but is free from the limitations mentioned above in the magnetic part of the analysis. This mixed approach has also an advantage in computational time with respect to the approach when ANCF is used in both mechanical and magnetic parts of the analysis.

Conclusions
This paper explores a possible way to simulate the magnetic shape memory effect with large deformation based on the absolute nodal coordinate formulation. The surrounding space is accounted for, and the nonlinear bidirectional coupling terms between the magnetic field and the mechanical field are considered. The newly constructed ANCF magneto-mechanical element is capable of describing the full generation process of the magnetic-field-induced strain and capturing the special features of the magnetic shape memory effect including field-strain/magnetization hysteresis and the superelastic effect. The difference between the results obtained using the ANCF element and the classical element is discussed, and the cause for the deviation is analyzed using a pure two-domain magnetostatic problem. The conclusion is drawn that the ANCF element should be used carefully with appropriate constraints for the two-domain MSME simulation such that the corresponding gradient discontinuity on the interface is suitably satisfied. A mixed FEs strategy, which adopts ANCF FE in the displacement field solver and classical FE in the magnetic field solver, is proposed as an alternative option to fix the problem. The mixed FEs approach shows high agreement with the classical FE on tracking the generation of MFIS in the sample under a varying magnetic field. Additionally, the method predicts the dynamic response of the MSMA sample under a cyclic compression. With relatively low computational cost, the mixed FEs approach can be an attractive way to integrate the ANCF method in the MSME analysis of a complicated specimen with large deformation. Future work will focus on the handling of the gradient jump issues in the two-domain magnetic shape memory effect simulation. To get a clear understanding on the boundary conditions and the interface continuity conditions, the variational principle of a general two-domain magnetostatic problem is presented here. The problem is shown in Fig. 15. One can assume that 1 is a permanent magnet domain with a constant magnetization M and 2 is a free space domain. Therefore, the constitutive relations are shown in Eq. (44). Recalling the illustration in Sect. 4.2, the governing differential equation is It should be noted that Eq. (B.1) is derived based on the constant magnetization assumption and holds on the whole domain. Without loss of generality, the boundary conditions are wherex andŷ are the unit direction vectors of the axes, and n is the unit outward normal vector on boundary 2 . The first and second term are classified into the Dirichlet BC and the natural BC, respectively. The natural BCs can be implicitly satisfied in the framework of variational principle as shown in the following formulation. Specifically, in the case of γ = 0 and q = 0, the second term degenerates to the magnetic insulation condition. Besides the BCs as shown in Eq. (B.2), the interface conditions as shown in Eq. (45) should be satisfied. Recalling the constitutive relations Eq. (44), the first term in Eq. (45) is where ψ + and ψ − are the potentials on the two sides of the interface, and n d is the interface normal vector directing from 1 (B.7) Based on the partial integration rules, and the divergence theorem, Eq. (B.7) can be written as where ∂ 1 = 1 + + d and ∂ 2 = 2 + − d are the full boundaries of domain 1 and 2 . Noting that ψ is constant on 1 , and therefore the spatial gradient vanishes, Eq. (B.8) can be simplified to  Clearly, the normal gradient of the magnetic potential on the element edge depends not only on the DOFs of the related two nodes but also the DOFs of the other two nodes. In other words, the normal gradient is not continuous between the two adjacent elements by default. Therefore, the element is not a C 1 -continuous element. As similar statement can also be found in [12]. Since the ANCF element cannot lead to a continuous normal gradient of the magnetic potential at the arbitrary non-node point on the element side, the derivation of the ANCF analysis result essentially comes from the use of gradients as DOFs at the nodes.