Efficient micromagnetic–elastic simulations based on a perturbed Lagrangian function

Ferromagnetic materials are gaining an increasingly important role in a wide range of industrial branches. Improving magnetic properties is essential to meet the resulting increased demand, which requires a deep understanding of magnetic as well as magneto-mechanical mechanisms. These mechanisms rely on different competing energy contributions¹ within the magnetic solid and result in the total energy of the magnet. This total energy can be minimized by adopting certain magnetization patterns. These magnetization patterns (i.e., magnetic domains) describe energetically stable states (energy minima of the total energy) and depend on external factors, such as mechanical stresses, external magnetic fields, or other influences, see Refs. 2–4. Their evolution can be described precisely by the micromagnetic theory and the underlying Landau–Lifshitz equation. The latter is the result of the findings by Landau and Lifshitz⁵ and Brown.⁶ Phenomenologically motivated modifications on the Landau–Lifshitz equation by Gilbert⁷ lead to the Landau–Lifshitz–Gilbert (LLG) equation, which will be applied within the context of this work in its Gilbert form. Today, the LLG equation is applied to numerically investigate the magnetic properties of microstructures and analyze the influence of varying chemical and geometrical compositions as done by, e.g., Refs. 8–12. This work is limited to isothermal processes far below the Curie temperature. Therefore, a constant length of the magnetization vectors $M= M s m$ is assumed, where $M s$ denotes the material specific saturation magnetization and $m$ is the so-called unit director. Formally, this can be expressed as $‖ M‖= M s$ and $‖ m‖=1$ (in the context of this work, $‖∙‖= ∙ 2$ hold true). The length or unit constraint allows the magnetization vectors to only move in a lateral direction and led in the past to different fulfillment approaches. A very straight forward way to enforce the unit constraints is to introduce the magnetization vectors in terms of Cartesian coordinates and add a penalty term to the formulation.^13–16 This penalty term penalizes the deviations of the magnetization vectors from the unit constraint, depending on an initially defined penalty parameter. A slightly different technique, also based on penalizing the deviation of the magnetization vectors from the unit length, is introduced in Refs. 17–19. Here, the constraint is controlled by a Lagrange multiplier that recalculates the intensity of the penalization in each iteration step. A solution procedure that relies on the staggered solution of the system of equations and the renormalization of the magnetization vectors was proposed in Ref. 20. Other methods that fulfill the unit constraint a priori involve a reformulation of the LLG equation, e.g., in terms of spherical coordinates as done in Refs. 21–24 or in terms of Rodrigues parameters in combination with an exponential map as presented in Refs. 25 and 26. Recently, some of the approaches mentioned above have been compared and their individual advantages and disadvantages have been highlighted.²⁷ The other approach was proposed by Alouges²⁸ where conformal finite elements are used to restrict the magnetization vectors onto the unit sphere. Since this work deals with finite element based micromagnetics, the above presented schemes are also based on the same numerical scheme. However, since other numerical discretization methods have also proven to be efficient in micromagnetism, a selection of works is presented below. Open source programs such as OOMMF²⁹ or MuMax3,³⁰ both based on the finite difference method, are widely used. A key advantage of the finite difference method is the ability to perform fast GPU-based simulations since the underlying discretization is a uniform grid. However, this also highlights the disadvantages of the method. In particular, complex microstructures often require local refinements, which can be realized much easier using the finite element method.³⁰ For a more comprehensive literature review, we also mention here the novel approach of neural networks in micromagnetism.^31,32 Neural networks offer a fast ability to predict magnetization dynamics of complex structures.

The present work deals with the well-known static condensation of the Lagrange multiplier on the element level and, thus, with the reduction of the degrees of freedom within each element node by one. The static condensation has been intensively applied within typical engineering problems as contact formulations or mixed elements for incompressibility.³³ However, in the context of micromagnetic simulations, to the best of our knowledge, the static condensation of the system of equations has only been briefly described in Ref. 34. The micromagnetic models utilized within this paper are implemented into the framework of AceGen/AceFEM.^35–37 

The paper is structured as follows: First, the required field equations are introduced and the corresponding boundary conditions are defined. Second, an overview of the most basic contributions to the enthalpy function is given. Third, the discrete element matrices (i.e., system vector and matrix) of the whole system are sketched, followed by the static condensation procedure of the perturbed Lagrange multiplier. Fourth, the workability of the proposed method is validated using suitable numerical examples for the static, dynamic, and micromagnetic–mechanically coupled verification. Additionally, the computational costs are compared against the non-condensed perturbed Lagrangian formulation. Finally, the results are evaluated and discussed.

The equations introduced above are solved numerically using the finite element method (FEM). The general approach of the FEM is to transfer the equations Eqs. (1), (3), (5), and (7) into their weak forms and build their corresponding linear increments. The resulting system of equations can be subsequently solved using Newton-like iterative methods. A comprehensive overview of FEM can be obtained by the study of, e.g., Refs. 33, 35, and 42. The following formulation results in a fully coupled system of equations which, in the context of this work, is solved in a monolithic way by applying Newton’s method. However, in order to account for mechanical effects in micromagnetism, there are other methods besides the fully coupled and monolithically solved method. In the work of Xiao et al.,⁴³ two such methods are compared. These are the so-called unidirectionally coupled (UC) method and the bidirectionally coupled (BC) method, also known as staggered approach. Both models are similar in that mechanical and magnetic quantities are calculated separately. However, the solution routines differ fundamentally from each other. The UC method first solves the mechanical quantities and then provides them for the calculation of the magnetic quantities so that the magnetization can be influenced by the previously calculated strains. The BC method extends this procedure by providing the magnetic quantities during the solution of the mechanical quantities and repeats this process until the solution satisfies a certain criterion, e.g., both residual norms are smaller than a predefined tolerance $ε tol$⁠. Since the UC method couples only in one direction, it is suitable only for weakly coupled phenomena, such as weakly magnetostrictive materials. The BC method considers both coupling directions, hence it is also suitable for the simulation of strongly coupled phenomena. However, the here presented approach is based on a fully coupled system of equations so that the strength of the coupling effects does not matter.

This section can be summarized as follows: The Lagrange multiplier is included as an additional degree of freedom within the coupled system of equations [Eq. (25)] that is solved incrementally for each time step using Newton’s method. Thus, the value of the Lagrange multiplier (as well as all other degrees of freedom) changes with each additional global iteration of Newton’s method until the time step converges with reaching $‖ I R _ e , it‖≤ 10 − 10$⁠. If the time step reached convergence, the values of the Lagrange multiplier correspond to the required penalization intensity of the magnetization vector. Since the Lagrange multiplier in the system of equations is an additional degree of freedom, it belongs to the coupled system of equations and needs to be solved. This increases the computational effort significantly. By perturbing the Lagrange multiplier, i.e., adding the term $λ 2 2 κ L$⁠, non-zero entries on the corresponding main diagonal of the system matrix are generated by differentiating twice. These entries are introduced in Eq. (25) as $I K _ λ λ it$⁠. This allows the inversion of this sub-matrix to $I K _ it - 1 λ λ$ and, finally, the static condensation of the Lagrange multiplier on element level in Eqs. (26) and (27). The update of the condensed Lagrange multiplier follows the rule introduced in Eq. (28). Thus, the intensity (value of the Lagrange multiplier) of the penalty is determined anew within each iteration.

In this section, the functionality of the presented method is demonstrated. This includes the preservation of the unit length during the simulation, as well as the physically correct representation of micromagnetic–mechanically coupled effects. First, the $μ$Mag standard problem $#$4 is used to verify the ability of the framework presented here to handle dynamic problems. Subsequently, in order to verify the mechanical coupling, a micromagnetic–mechanically coupled boundary value problem is investigated, which was taken from literature.⁴⁴ Finally, an Ni thin film is analyzed with respect to the evolution of its magnetic domains under the influence of a defect structure. All simulations included in this work follow the undimensionalization procedure described in Reichel et al..²⁷ 

A well-known example to validate dynamic micromagnetism is the standard problem $#4$⁠.⁴⁵ The considered boundary value problem corresponds to a cuboid shaped thin film of dimensions $500×125×3 nm 3$⁠. This thin film is considered to mimic the material behavior of Permalloy with the following material parameters assigned to the simulation: A $exc=1.3× 10 − 11$ J/m, M $s=8.0× 10 5$ A/m, K $ani=0$⁠, and $α=0.02$⁠. The initial magnetization configuration corresponds to the so-called S-state, compare Fig. 1(a). This S-state can be obtained by a full saturation along the [1,1,1] axis caused by a magnetic field and a subsequent relaxation by slowly removing this field. After the relaxation a constant magnetic field $μ 0 H ¯= [ − 24.6 , 4.3 , 0.0 ] T$ mT is applied to the thin film, causing an instantaneous reversal (The second field, defined in the original problem, is not considered in this work since it proceeds analogously to the presented one.). The thin film is discretized with $100×30×3$ quadratic tetrahedron elements along the corresponding axis. The time integration was accomplished by applying time increments $Δt= 10 − 4$ ns.

For the comparison of the time evolution of the averaged magnetization components, the results generated by the software MagTens⁴⁶ provided on the $μ$Mag homepage have been used. The small deviations of the here obtained solution to the standard problem $#$4 can be eliminated with further decreasing the time increments or applying a different and more sophisticated time integration scheme. The latter one will be accomplished within the next publications. Also, the magnetization distribution in Fig. 1(b), when m $1$ equals zero for the first time indicates very good agreement with the patterns provided on the $μ$Mag homepage.⁴⁵ In addition to the validation of the dynamical properties of the proposed model, a comparison with the implicit midpoint rule^47–49 proved the proposed method to be equal for the convergence behavior and the required simulation time. A direct comparison of the relative residual norms of the condensed perturbed Lagrange method and the Midpoint rule is shown in Table II for the time step no. 100 of 101 at time 0.1485 ns. The applied increment corresponds to $Δt=0.15× 10 − 2$ ns. The direct comparison of these residual norms after the respective iterations in the corresponding time step shows an equally good convergence behavior of both methods.

This section serves as a validation of the micromagnetic-mechanically coupling. Hence, the magnetostrictive material nickel (Ni) is used for the subsequent investigations. In this context, only the effect of mechanical strains on the effective hysteresis behavior is considered. This effect is generally known as Villari or inverse magnetostrictive effect.⁴⁰ 

While magnetostriction describes the change of the shape of the magnet under varying magnetization, inverse magnetostriction describes the change of the magnetic properties under the influence of shape-changing effects, such as external stresses or strains acting on the magnetic solid. This phenomenon has been observed in nanowires where both experimental and numerical results have been compared.⁵⁰ However, the numerical simulations performed are purely micromagnetic in nature. In contrast, a fully micromagnetic–mechanically coupled material model was used by Liang et al.⁴⁴ to correctly incorporate the influence of the Villari effect. They could show in a comparison that the experimental results (Bur et al.⁵¹) and numerical results agree well. Hence, the boundary value problem⁴⁴ is adopted within this work to verify the method presented above, i.e., the magnetic Ni nanostructure with the dimensions $300×100×35$ nm $3$ is fixed on top of a pure mechanical carrier substrate with dimensions $1000×1000×100$ nm $3$⁠. Since magnetic fields, such as stray fields, evolve around magnetic materials, the free space surrounding the magnetic matter must be considered. This can be done by various hybrid formulations coupling the FEM to surface discretization methods^38,39,52 or in a very straight forward way by discretizing some of the surrounding space with finite elements. A rule of thumb for estimating the size of the area to be discretized is given in Ref. 53. A sketch of the boundary value problem is given in Fig. 2, while the corresponding material parameters can be taken from Table III. The Gilbert damping parameter is set to be $α=0.5$⁠.⁴⁴ 

To analyze the influence of mechanical strains on the Ni structure, mechanical boundary conditions are applied on the surface of the substrate such that the relative strain $με=ε× 10 − 6=( ε 22− ε 11$⁠) within the substrate yield $−$1210, 0, and 1060, where $ε 11=−ν ε 22$ and $ν=0.2$⁠. This implies the homogeneous support of the substrate is prescribed as $u x( x=−500 nm)= u y( y=−500 nm)= u z( z=−100 nm)=0 nm$⁠, and the mechanical loads correspond to the values given in Table IV.

To obtain hysteresis loops, the Ni film is subjected to an alternating external magnetic field, following the load path shown in Fig. 3(a) with a maximal field intensity of $μ 0 H ¯= [ 0.1 , 0 , 0 ] T$T. Initially, the mechanically unstressed nanostructure is subjected to the alternating field, which generates a hysteresis, see Fig. 3(b). Subsequently, the carrier substrate is exposed to the mechanical displacements so that the Ni nanostructure also deforms. These deformed states are now treated with the alternating magnetic field. The corresponding hystereses are shown along with the hysteresis of the undeformed structure in Fig. 3(b) and indicate good agreement with the results from the literature.⁴⁴ The influence of the mechanical loads on the hysteresis behavior of the material can be clearly seen in the hystereses shown in Fig. 3(b). Negative relative strains produce a much more pronounced coercivity compared to an undistorted structure, while positive strains have the opposite effect. The origin of this behavior lies in the so-called Villari effect also known as inverse magnetostriction.⁴⁰ 

The Lagrange multipliers considered here do not necessarily require continuity across element boundaries, thus these degrees of freedom can be approximated both continuously and discontinuously. The discontinuity of the degrees of freedom ultimately allows their condensation. Classical Lagrange multipliers can lead to numerical difficulties due to zero entries on the main diagonal. Thus, in this work, the Lagrange multiplier has been extended by a quadratic perturbation term $− λ 2 2 κ$ so that the main diagonal contains non-zero entries. Thus, numerical difficulties in the solution process can be circumvented.^33,41 In addition, the number of constraints must be at least less than or equal to the number of quantities to be constrained, i.e., number of $d λ≤ d m$⁠. Hence, different interpolation orders for the Lagrange multiplier are examined in the following for their performance.

The considered discretization of the individual degrees of freedom $Ξ:={ u _,φ, m _}$ and $λ$ are restricted in this work to tetrahedral elements, denoted as $T$⁠. The interpolation order of the tetrahedral elements is indicated as a subscript index as $T i$ with i = 0,1,2 for constant, linear, and quadratic polynomials. The degrees of freedom $Ξ$ will be interpolated using linear and quadratic interpolation functions and, therefore, denoted as $T 1 Ξ$ or $T 2 Ξ$⁠. Since the here considered Lagrange multiplier can be interpolated continuously, discontinuously, or discontinuously with condensation, additional subscript indices indicate their continuity as $T i , c λ$⁠, $T i , d λ$⁠, or $T i , d , cond λ$⁠. This allows for six different interpolations of the considered micromagnetic formulation, if the interpolation order of the Lagrange multiplier is equal or an order below the other degrees of freedom, as shown in Table V. The time comparison is performed using the undeformed Ni structure from Sec. V B. All six variants of interpolation are used to calculate the demagnetization branch of the hysteresis loop of the Ni stripe. Compared to the non-condensed parameters, there is a significant speed up of the simulations, while the results remain the same.

In this example, the magnetization behavior of three different Ni thin film boundary value problems is analyzed. All considered thin films have the dimensions $1000×500×50 nm 3$⁠, and the material properties are already given in Table III. The first considered thin film (TF1) is free of mechanical imperfections and defects and thus represents a kind of ideal reference. Subsequently, a defect (rectangular on the surface) is placed within the thin film structure (TF2). This defect completely intersects the film and has dimensions $135×85×50 nm 3$⁠. For the last boundary value problem considered, mechanical displacements u $x=0.2$ nm are applied to the inner surfaces of the defect structure (TF3). An example of the latter is shown in Fig. 4. The differences to the previously described boundary value problems are thus mainly in the defect structure and the mechanical displacements within the defect.

First, all structures are initialized with a normally distributed random magnetization distribution and subsequently relaxed to their energetically most favorable and stable state. The initial as well as relaxed states of the first boundary value problem, i.e., TF1, is shown in Fig. 5. Its relaxed state corresponds to a double vortex formation with well balanced domain formations. Since the initial configurations of the defect containing thin films (TF2 and TF3) are very similar to Fig. 5(a), besides of the hole within the structure, they are omitted here. Their relaxed configurations also corresponds to a double vortex formation with slightly differing characteristics, caused by the different types of defects (hole and hole with displacement), as presented in Fig. 6. The distribution of the domains of TF2 and TF3 appears to be clearly affected by the defect. As can be seen in Figs. 6(a) and 6(b), the left positively aligned m $1$ domain is pinned to the defect so that a uniform size distribution of all domains, as for TF1 [Fig. 5(b)] cannot be achieved. However, the domains in TF2 each have an equally sized, antiparallel oriented domain within each vortex, while the mechanical influence in TF3 aligns them asymmetrically. The result can be seen clearly when the vortex cores are connected to each other by a line. This line is parallel to the x $1$-axis in Fig. 6(a) for TF2, while it is slightly tilted in Fig. 6 for TF3. Furthermore, the magnetization to the right of the defect points entirely in the x $2$-direction. Without the defect a deflection in the x $1$-direction would still be expected at this point. Hence, the exchange interactions from the left to the right of the defect can be treated as exchange decoupled.

To obtain the hysteresis loops of the three thin film structures, an alternating magnetic field is applied, similar to that described in Fig. 3(a). The hysteresis loops of all three considered boundary value problems are very similar, with little differences [Fig. 7(a)]. The maximal coercivity is reached by the nanostructure containing the defect without mechanical influences (TF2) even though it is very close to the value obtained for the displacement influenced defect nanostructure (TF3), as presented in Fig. 7(b). The reason for these, even if minor, different values in coercivity is the pinning of the magnetization vectors between the edge of the thin film and the defect structure. This delays the cascade like reversal of the domains compared to the defect-free structure.

This work addresses the simulation of micromagnetic–elastically coupled systems. As in micromagnetism, the coupled theory used here involves the numerically challenging task of length preservation of the magnetization vectors during the simulation. A very straightforward method to enforce this is to use Lagrange multipliers. Since Lagrange multipliers lead to zero entries on the main diagonal in the system matrix, which can lead to problems in solving the system of equations, it is convenient to use a perturbed Lagrange multiplier. Here, a quadratic penalty-like term is added to the energy potential which provides non-zero entries on the main diagonal of the system matrix. This characteristic allows a static condensation of the Lagrange multiplier at the element level and thus reduction of the local system of equations by one degree of freedom per node. This methodology is a well-known approach to engineering problems, such as contact problems and mixed formulations, which to our knowledge has not yet been applied to micromagnetism. It should be mentioned that the magnetization vectors in our simulations still corresponded to the unit length, even for the perturbation of the Lagrange multiplier. To demonstrate the general workability of the previously introduced method, several test cases have been considered and evaluated, including the well-known $μ$Mag standard problem $#4$ for dynamic simulations as well as micromagnetic elastic simulations of a nickel nanostructure. These comparisons provided good agreements and proved the introduced framework to be able to simulate dynamic, quasi-static, as well as magneto-elastically coupled boundary value problems. To quantify the performance of the method, comparative simulations with non-condensed perturbed Lagrange multipliers were performed. Here, the condensation of the Lagrange parameter has shown a reduction in computational costs of up to 29.8 $%$ for a comparison with the continuously discretized parameter. Based on this verification, an Ni thin film with defect structures was analyzed regarding the influence of the defect on the domain formation. The obtained results suggest that a correlation exists between geometric defects such as holes, strains, and the resulting domain arrangement. However, this influence will be studied much more intensively in later work, including how to target defects and generally distortions in microstructures to gain eventually advantages for effective hysteretic behavior.

We gratefully acknowledge the financial support of the German Research Foundation (DFG) in the framework of the CRC/TRR 270, project A07 “Scale-bridging of magneto-mechanical mesostructures of additive manufactured and severe plastically deformed materials,” and also Z-Inf therein, Project No. 405553726.

The authors have no conflicts to disclose.

Maximilian Reichel: Conceptualization (lead); Investigation (lead); Methodology (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Rainer Niekamp: Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). Jörg Schröder: Funding acquisition (lead); Resources (lead); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.