Accurate evaluation of the configurational forces in single-crystalline NiMnGa alloys under mechanical loading conditions

Acta Materialia 105 (2016) 306e316

Contents lists available at ScienceDirect

Acta Materialia journal homepage: www.elsevier.com/locate/actamat

Full length article

Accurate evaluation of the configurational forces in single-crystalline NiMnGa alloys under mechanical loading conditions Jiong Wang School of Civil Engineering and Transportation, South China University of Technology, 510640 Guangzhou, Guangdong, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 27 August 2015 Received in revised form 28 November 2015 Accepted 8 December 2015 Available online 8 January 2016

To properly predict the mechanical response of a single-crystalline NiMnGa sample, the configurational forces on the twin interfaces in the sample have to be calculated accurately. However, the conventional displacement-based finite element (DB-FE) method usually yields large numerical errors in calculating the contact stress values on the twin interfaces, which will further influence the evaluations of the configurational forces. In this paper, an interface element (IE) method is incorporated into the conventional FE framework. Through this FE-IE method, the contact stress is treated as one of the primary unknowns on the twin interfaces. By solving the governing system of the model numerically, the contact stresses on the twin interfaces can be directly obtained, based on which the configurational forces can be calculated. The obtained results are further used to predict the global mechanical response of the NiMnGa sample. For a simple example, the theoretical values of the resultant contact forces and the configurational forces on the twin interfaces can be determined. By comparing the numerical results obtained from the DB-FE method and the FE-IE method with the theoretical values, it is found that the FE-IE method can generate more accurate numerical solutions. Furthermore, the influence of the crosssectional shapes of NiMnGa samples on their mechanical response is investigated. © 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

Keywords: NiMnGa alloys Twin interfaces Configurational force Interface element method Mechanical properties

1. Introduction Single-crystalline NiMnGa alloys [1,2] are one of the most known and widely explored magnetic shape memory alloys (MSMAs). This group of materials can undergo significant strains (6e10%) under the applications of magnetic fields or mechanical loads. The response frequency of NiMnGa alloys can reach 1e2 kHz. These unusual properties make NiMnGa the appreciate candidates for various applications, e.g., actuators, sensors, energy harvesters and vibration dampers. During the past years, the magnetic and mechanical properties of NiMnGa alloys have been extensively investigated, through both experimental measurements [3e11] and theoretical modelings [12e18]. The underlying mechanism of the unusual properties of NiMnGa alloys has been found to be the fieldor stress-induced martensite variant reorientation. Generally, the NiMnGa samples adopted in the experiments are in the five-layer modulated (5M) martensite phase, which has a tetragonal lattice structure with one short-axis of length c ¼ 0:560 nm and two long-axes of length a ¼ 0:595 nm [4].

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.actamat.2015.12.018 1359-6454/© 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

Totally, three variants with different orientations can be identified in this phase (cf. Fig.1a). Under the applications of magnetic fields or mechanical loads, variant reorientations can take place, which leads to macroscopic deformations of the samples. However, if the three variants coexist in a sample, some textured twin structures will be formulated. In this case, the variant reorientation process in the sample is largely obstructed. Thus, in most of the experimental works [4,5,7,8,11], only the reorientations between two variants were carried out to obtain large reorientation strains. In this paper, the stress-induced variant reorientations in NiMnGa samples will be studied. One typical experimental procedure is illustrated in Fig.1b. Initially, the sample is in a single variant state (say variant two). A compressive force is then applied perpendicularly to the short-axis of the original variant, which induces the reorientation from the initial variant to a new variant (say variant one) with the short-axis parallel to the external loading. Accompanying the variant reorientation, a significant strain of the sample along the direction of the external force can be observed. Finally, the NiMnGa specimen attains another single variant state and the variant reorientation process is completed. One key feature that has been observed in the experiments is that the variant reorientation is realized through the nucleations and movements of

J. Wang / Acta Materialia 105 (2016) 306e316

307

Fig. 1. (a) The lattice structures of the three NiMnGa variants in the 5M martensite phase; (b) typical procedure of stress-induced variant reorientation in a NiMnGa sample.

twin interfaces in the NiMnGa sample (cf. Fig. 2). The twin interfaces are some flat planes with fixed orientations, which are arranged parallelly in the sample. The appearance of twin interfaces can also induce the deflection of the sample (cf. Fig. 6 in Ref. [14]). In our previous work [19,20], a variational constitutive model has been proposed to simulate the field- and stress-induced variant reorientations in single-crystalline MSMA samples. The governing system of the model is composed of the magnetic and mechanical field equations, the evolution laws for some internal variables and the criteria of twin interface movements. Especially, in the twin interface movement criteria, the expression of the configurational force on a twin interface is presented, which plays the role of ‘driving force’ for the movement of the twin interface. Only the mechanical part of the governing system derived in Ref. [19] will be considered in the current paper. To properly predict the mechanical behavior of a NiMnGa sample, the configurational forces on the twin interfaces need to be calculated accurately. However, it was found that the values of the configurational forces obtained from the conventional displacement-based finite element (DB-FE) method usually exhibit large numerical errors. In fact, for the purely mechanical loading case, the configurational forces

mainly depend on the contact stresses on the twin interfaces. By using the DB-FE method, the nodal displacements in the sample are first determined. Based on the displacement values and by virtue of the constitutive relations, the deformation gradients and the stress tensors can be evaluated in the elements, which are then mapped onto the twin interfaces to calculate the contact stresses and the configurational forces. Such a process may result in large errors in the numerical values of the contact stresses and the configurational forces. Especially, the compatibility condition of the stress tensors across the twin interface can be violated. The problem of contact stress between different media has been studied for a long time in the fields of civil engineering structures, composite structures and rock mechanics. In the literature, several approaches have been proposed to improve the numerical accuracy in evaluating the contact stresses, including the interface element method [21e25], the penalty function method [26,27], the boundary element method [28,29], etc. In current work, we shall adopt the interface element (IE) method proposed in Ref. [24]. In this method, the contact stress on a twin interface is treated as one of the primary unknowns. Additional equations can be proposed by considering the continuity of the displacements across the twin interfaces, which can be easily integrated in the FE system matrix.

Fig. 2. Development of the surface morphology of a NiMnGa prism-shaped specimen during the variant reorientation process induced by compression (taken from Fig. 2 in Ref. [6]).

308

J. Wang / Acta Materialia 105 (2016) 306e316

Through numerical computations, the contact stress values on the twin interfaces, together with the nodal displacements in the sample, can be directly obtained, based on which the configurational forces on the twin interfaces can be calculated. To verify the efficiency of the FE-IE method, a simple example is considered. In this example, the theoretical values of the configurational forces on the twin interfaces can be derived. By comparing the numerical results obtained from the DB-FE method and the FE-IE method with the theoretical results, it is found that the FE-IE method can generate more accurate numerical solutions. The structure of this paper is arranged as follows. In Section 2, the mechanical part of the governing system derived in Ref. [19] is introduced. In Section 3, the numerical scheme of the interface element method is presented, which is then integrated into an iterative numerical algorithm to solve the governing system. In Section 4, a simple example is considered to evaluate the numerical results obtained from the conventional DB-FE method and the FE-IE method. The influence of the cross-sectional shape of NiMnGa samples on the configurational forces will also be studied. In Section 5, it is discussed how to apply the IE method to evaluate the configurational forces under combined loading conditions. Besides that, some other open problems in current modeling approach are introduced. Finally, some conclusions are drawn.

2. The NiMnGa sample and the governing equations

(1)

where Fe is the elastic deformation gradient and Gi are the transformation deformation gradients in Uir ði ¼ 1; 2Þ. Furthermore, the following polar decompositions of Gi is considered

Gi ¼ ℚi ℙi ; i ¼ 1; 2;

(2)

where ℚi are orthogonal tensors representing rigid rotations and ℙi are symmetric tensors representing the principle stretches. For the two variants of 5M NiMnGa martensite, ℙi can be chosen as

ℙ1 ¼ I; ℙ2 ¼ I þ se1 5e1  s0 e2 5e2 :

(3) s0

RotGi ¼ 0; in Uir ; G2  G1 ¼ a5N; on vUT ;

where I is the identity tensor and s and are two small parameters representing the transformation stretches. As the transformation strain in variant reorientation is isochoric, s and s0 satisfy ð1 þ sÞð1  s0 Þ ¼ 1. Suppose that the transformation strain satisfies the following

(4)

where N is the unit normal vector on vUT . From the compatibility conditions (4) and by virtue of (2) and (3), it can be proved that ℚi is a constant tensor in each connected variant region [19]. Besides that, the normal vector N also has certain fixed orientations [19]. For simply twinned NiMnGa samples (i.e., all the twin interfaces in the sample are parallel planes), the rotation tensors ℚi ði ¼ 1; 2Þ and the normal vector N can be chosen as ℚ1 ¼ I and

ℚ2 ¼ þ

2½s þ 1 1 þ ½s þ 12

s½s þ 2 1 þ ½s þ 12 2

e1 5e1 þ

e1 5e2 

2½s þ 1

e2 5e2 þ e3 5e3

1 þ ½s þ 12

s½s þ 2 1 þ ½s þ 12

e2 5e1 ; 3

(5)

sþ1 1 6 7 N ¼ ±4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e1 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e2 5: 1 þ ½s þ 12 1 þ ½s þ 12 With the above preparations, the total energy functional of the NiMnGa sample can be proposed in the following form

G ðx; P; XÞ ¼

Suppose that the NiMnGa sample has the reference configuration Ur 3ℝ3 . Only two variants of the 5M martensite phase, say variant one and variant two (cf. Fig. 1), are considered in this model. For future convenience, a set of mutually orthogonal unit vectors fe1 ; e2 ; e3 g is constructed, where e1 and e2 direct along the c-axes of variant one and variant two, respectively. During the variant reorientation process, the two martensite variants coexist in the sample, which occupy the regions (closed sets) U1r and U2r . It is clear 1 2
represents the interior of a that U1r ∪U2r ¼ Ur and
Ur ∩
Ur ¼ ∅ (ð,Þ set). Denote vUT ¼ U1r ∩U2r , which just represents the twin interfaces in the sample. Notice that if multiple twin interfaces coexist in the sample (i.e., vUT ¼ vU1T ∪vU2T ∪/∪vUnT ), region Ur will be separated into several disjoined subregions. Denote X and x as the coordinates of a material points in the reference and current configurations, respectively. The deformation gradient tensor is given by F ¼ vx=vX. During the variant reorientation process, the material undergoes both the transformation (reorientation) strain and the elastic strain, thus the following multiplicative decomposition of F is adopted

F ¼ Fe Gi ; i ¼ 1; 2;

compatibility conditions in Ur [30].

Z 2 X i¼1

Z rr fðF; Gi ÞdV 

tA ,xdA;

(6)

vUr

Uir

where rr is the material density in the reference configuration, f is the elastic energy per unit mass and tA are the external traction applied on per unit referential area. It should be noted that besides the spatial placement field x, G also depends on the variant state distribution P :¼ fU1r ; U2r g in the sample. By calculating the variation of G with respect to x, one obtains the mechanical field equation together with the boundary and interface conditions, which are given by

DivS ¼ 0; in U1r ; U2r ; ST N ¼ tA ; on vU1r ; vU2r ; T ESF N ¼ 0; on vUT ;

(7)

where S ¼ rr vfðF; Gi Þ=vF is the nominal stress tensor, vUir are the sample surface areas associated with Uir ði ¼ 1; 2Þ and the square bracket E,F ¼ ð,Þ1  ð,Þ2 denotes the jump of the enclosed quantity across the twin interfaces. To derive the criteria of twin interface movements in the NiMnGa sample, one needs to calculate the variation of G with respect to the variant distribution P and consider the energy dissipation effect during the variant reorientation process. The detailed derivation procedure of the criteria can be found in Ref. [19], where the magnetic properties of NiMnGa alloys have also been taken into account. For the purely mechanical loading case, the twin interface movement criteria are given by

Z 

rr D x dAc < T

vUkT

8 > T > > > > < Vk ,Ns00 > > T > > > :

Z k

< Z

k

rr Dþ x dAc ;

¼ vUkT

k

rr Dþ x dAc 0 Vk ,N ¼ 0

vUkT

Z

¼

if Vk ,N < 0;

(8)

rr D x dAc ; if Vk ,N > 0:

vUkT

The above criteria are derived for each twinning plane vUkT in the

J. Wang / Acta Materialia 105 (2016) 306e316

sample (remember that vUT ¼ vU1T ∪vU2T ∪/∪vUnT ). In (8), T k is called the configurational force on vUkT (in the sense of Eshelby [31]), which has the following expression

Z T

k

Err fðF; Gi Þ  FN,ST NFdAc ;

¼

(9)

vUkT

D± x are the energy dissipation densities (per unit mass) during the forward and reverse variant reorientations, Vk is the configurational velocity of the twin interface vUkT and the normal vector N is supposed to be directing into the region of U2r . The quantities R D± ¼ rr D x dAc play the roles of blocking forces during the k vU variant reorientation processes. T The mechanical behavior of a single-crystalline NiMnGa sample can be described by using the field eqn (7) and the twin interface movement criteria (8). In fact, (7) and (8) just correspond to the mechanical part of the governing system derived in Ref. [19]. Remark: The expression of the configurational force (9) has been derived in the classical work of Abeyaratne [32]. In Ref. [19], we adopted the similar approach proposed in Ref. [32] and extended the previous results to model the response of MSMA samples under combined magnetic and mechanical loading conditions.

3. The numerical algorithm In Wang and Steinmann [20], we have proposed an iterative numerical algorithm to solve the governing system derived in Ref. [19]. For the mechanical part of the governing system (7) and (8), the algorithm proposed in Ref. [20] can be simplified. The flowchart of the algorithm is shown in Fig. 3. First, for each load step, the value of the external force tA applied on the sample surface is specified. Next, the initial variant distribution in the sample is assigned. Correspondingly, the initial positions of the twin

309

interfaces in the sample are determined. Suppose that the twin interfaces in the sample are fixed, by solving the mechanical field eqn (7), one can determine the displacement field, the deformation gradient, the stress tensor and other physical quantities in the sample. The obtained numerical results are then substituted into the criteria (8) to check the stability of each twin interface. If the inequalities (8)1 are satisfied, the twin interface will stay at its initial position and the calculation procedure is completed. Otherwise, (8)2 will be used to adjust the positions of the twin interfaces and the above procedure will be repeated again. It can be shown that this iterative numerical algorithm is efficient in solving the PDE system (7) and (8). However, there is one critical issue that deserves further investigations. To apply the twin interface movement criteria (8), the central task is to calculate the configurational force T k on the twin interfaces. For that purpose, the contact stress ST N on each twin interface needs to be calculated accurately. However, as mentioned before, the conventional DB-FE method usually results in large numerical errors in evaluating the contact stresses and the configurational forces. Especially, the compatibility condition (7)3 between the adjacent elements across the twin interface can be violated. In our previous work [20], we did not provide in-depth investigations on this problem. In current work, the interface element method proposed in Ref. [24] will be adopted to tackle this problem. The numerical scheme of this method is presented below, where the linear tetrahedral elements and the linear triangular elements are taken as examples. Suppose that one NiMnGa sample is meshed by using the linear tetrahedral elements. Inherited from the volume mesh, the twin interfaces in the sample are meshed with the linear triangular elements. Based on each triangular element Ge on the twin interface, the interface element can be constructed (cf. Fig. 4), which has the six nodes labeled as ni ði ¼ 1; /; 6Þ. The nodal displacement vector and external force vector of one interface element are denoted as

T  u ¼ u11 ; u12 ; u13 ; u21 ; u22 ; u23 ; /; u61 ; u62 ; u63 ; T 181  1 1 1 2 2 2 6 6 6 : f ¼ f1 ; f2 ; f3 ; f1 ; f2 ; f3 ; /; f1 ; f2 ; f3

(10)

181

During the deformation process, the nodes on the lower and upper faces of the interface element can undergo relative displacements, which is represented by

T  Du ¼ Du11 ; Du12 ; Du13 ; /; Du31 ; Du32 ; Du33

91

¼ ℂu:

(11)

The matrix ℂ is defined by

0

I3 ℂ ¼ @ 03 03

03 I3 03

03 03 I3

I3 03 03

03 I3 03

1 03 03 A ; I3 918

(12)

where I3 is the identity tensor and 03 is the zero tensor. With the relative nodal displacement Du, the relative displacement of the continuous element can be calculated through

Du ¼ ðDy1 ; Dy2 ; Dy3 ÞT ¼ ℕDu;

(13)

where ℕ contains the interpolation polynomials according to

0

N1 ℕ¼@ 0 0

Fig. 3. Flowchart of the iterative numerical algorithm for solving the PDE system (7) and (8).

0 N1 0

0 0 N1

N2 0 0

0 N2 0

0 0 N2

N3 0 0

0 N3 0

1 0 0 A : N3 39

(14)

and Ni ði ¼ 1; 2; 3Þ are the three shape functions of the triangular element. Due to the relative displacements, internal contact stress can

310

J. Wang / Acta Materialia 105 (2016) 306e316

Fig. 4. Illustrations of the mesh of the NiMnGa sample, the inherited mesh of a twin interface and the interface element on the twin interface.

also be induced in the interface element. By denoting the contact stress vectors on the nodal pairs as PðiÞ ¼ fP1i ; P2i ; P3i g ði ¼ 1; 2; 3Þ, the internal traction of the continuous element can be calculated by

T  P ¼ ℕP ¼ ℕ P11 ; P21 ; P31 ; P12 ; P22 ; P32 ; P13 ; P23 ; P33 :

(15)

Usually, P is further expressed in terms of the relative displacement Du through the constitutive equation [33,35]. However, in current interface element method [24], P will be treated as one of the primary unknowns. The linear element stiffness matrix K can be constructed by using the principle of virtual work. Induced by a virtual displacement, the amounts of internal and external works done in the interface element are given by

2

Z

6 dðDuÞ,PdA ¼ duT 4ℂT

dU ¼ Ge

Z

3 7 ℕT ℕdA5P;

(16)

Ge

dW ¼ du,f:

Then from the equation dU  dW ¼ 0 and the arbitrariness of du, one obtains

2 6 T 4ℂ

Z

3 7 ℕT ℕdA5P ¼ f :

(17)

Ge

Besides that, to ensure the continuity of the nodal displacement across the twin interface, the following constraint equation is proposed

Du ¼ ℂu ¼ 0:

(18)

Formally (17) and (18) can be written as

Z

0 B0 @ ℂ

ℂT Ge

1 ℕT ℕdA C u   f  ¼ : A P 0 0

(19)

The coefficient matrix on the left hand side of (19) is the element

stiffness matrix of the interface element, which can be integrated into the system matrix of FE calculations. Through conventional calculation procedure, both the nodal displacements in the sample and the contact stresses on the twin interfaces can be obtained. Because the contact stress vector P is chosen as the primary unknown, the accuracy of the numerical value of P can be improved. Some examples have been introduced in Ref. [24] to verify the efficiency of this finite element-interface element (FE-IE) method. In next section, this FE-IE method will be used to calculate the configurational forces on the twin interfaces in a NiMnGa sample.

4. Evaluation of the numerical results In this section, a simple example will be considered to evaluate the numerical results obtained from the DB-FE method and the FEIE method. The NiMnGa sample in this example has a rectangular parallelepiped shape with dimensions 20  4:25  4:25 mm3 in the reference configuration (cf. Fig. 5a). To conduct FE calculations, this sample is meshed into a number of tetrahedral elements. Through this mesh, the sample region is divided averagely into several subregions (cf. Fig. 5a). These subregions have the inclined band form and the interfaces between them are some parallel planes with the normal vector N given in (5)2. In fact, these subregion interfaces just provide the possible positions of the twin interfaces in the sample. In current example, only 17 interfaces will be taken into account, which are labeled as I8 ; /; I0 ; /; I8 in Fig. 5a. The region sandwiched between I8 and I8 is called the gauge region of the sample. The variant reorientation process is realized through the jump of the twin interface from its original position Ii to the neighboring interfaces Ii1 or Iiþ1. One benefit of this specific mesh is that the sample does not need to be re-meshed after each adjustment of the twin interfaces, thus the computational consumptions can be reduced. As the twin interface movement criteria (8) cannot describe the nucleations of new twin interfaces, the initially twin interface positions in the sample need to be assigned. By applying a uniformly distributed load tA on the YZ-faces along the X-axis (cf. Fig. 5b), the twin interfaces will move from the two ends to the middle portion of the NiMnGa sample, which corresponds to the variant

J. Wang / Acta Materialia 105 (2016) 306e316

311

be replaced (approximately) by the transformation deformation gradients G1 and G2 . By further using (7)3, the expression of T k can be reduced into

Z ST NdAc :

Τ k zðG1  G2 ÞN,

(21)

vUkT

Fig. 5. (a) Mesh of the NiMnGa sample, and the generated subregions and region interfaces; (b) the NiMnGa sample with the twin pair fI3 ; I3 g and illustration of the external loads.

reorientation from variant two to variant one. Due to the coexistence of the two variants, the configuration of the NiMnGa sample is deflected (cf. the experimental observations in Ref. [14]). Due to this deflection, the total moment of the external load t does not vanish. To ensure the equilibrium of the sample, another uniformly distributed tS is applied on the YZ-faces of the sample along the Yaxis (cf. Fig. 5b). The other faces of the sample are supposed to be traction free. It is further assumed that the two end sides of the sample can move freely in the YZ-planes, thus no other constraints will be generated due to the deflection of the sample. To obtain some concrete results, the following material parameters are adopted in the numerical calculations [3,8,9].

. . 4 2 rr ¼ 8:3  103 kg m3 ; D± x ¼ 8:063  10 N m ; . . k1 ¼ 1:70  1011 N m2 ; k2 ¼ 1:50  1011 N m2 ; . . k3 ¼ k4 ¼ 1:52  1011 N m2 ; k5 ¼ k6 ¼ 0:43  1011 N m2 ; s ¼ 0:0625;

(20) where k1  k6 are the elastic moduli related to the elastic energy density 4. Before the numerical calculations, some analyses will be given to derive the theoretical results for this example. As mentioned before, the central task in predicting the twin interface movements is to calculate the configurational force T k . For purely mechanical loading case, the expression of T k has been given in (9). It has been pointed out in Ref. [19] that the first term in the expression of T k is negligible compared with the second term (this fact has also been verified by the numerical results obtained in Ref. [20]). Besides that, under the applications of moderate external loads, e.g., sA ¼ jtA =AYZ j < 10MPa (AYZ denotes the area of the YZ-faces), the elastic strain in the NiMnGa sample is relatively small compared with the transformation strain. Thus, in the expression (9), the deformation gradients F on the two sides of the twin interface can

From (21), it can be seen that the value of T k directly depends on the resultant contact force on the twin interface vUkT . For this simple example, as the external loads on the sample surface have been specified (cf. Fig. 5b), it is easy to determine the resultant contact forces on the twin interfaces. For example, in the case shown in Fig. 5b, the resultant contact force on the twin interface I3 can be calculated by considering the equilibrium of the sample subregion on the right hand side of I3 , which equals fjtA j; jtS j; 0g. Next, we begin to introduce the numerical results obtained from the DB-FE method and the FE-IE method. First, the resultant contact forces and the configurational forces on fixed twin interfaces will be calculated. For simplicity, it is supposed that only two twin interfaces exist in the sample and they are arranged symmetrically at the positions fIi ; Ii g ði ¼ 1; /; 8Þ. The sample region between Ii and Ii is in the state of variant two and the other region is in the state of variant one. Due to symmetry, only one of the twin interfaces in the sample needs to be considered (e.g., the right twin interface Ii). For the purpose of illustration, the compressive stress values is chosen as sA ¼ 3 MPa, then the axial component (along the X-axis) and the lateral component (along the Y-axis) of the resultant contact forces on the different twin interfaces are shown in Fig. 6. The theoretical values of the resultant contact forces are also plotted as the dashed lines in Fig. 6 for comparison. From Fig. 6a and b, it can be seen that the results obtained from the DB-FE method exhibit large numerical errors. The contact force components on the two sides of the twin interfaces are not equal to each other, which does not satisfy the compatibility eqn (7)3. On the other hand, from Fig. 6c and d, it can be found that the resultant contact forces obtained from the FEIE method are exactly equal to the theoretical results. Besides that, the compatibility eqn (7)3 is also satisfied. Thus, the FE-IE method yields more accurate numerical results on the twin interfaces. Based on the numerical results of the resultant contact forces, the configurational forces on the twin interfaces can be calculated. In Fig. 7, we show the evolutions of the configurational forces on the different twin pairs accompanying the variation of sA. It can be seen that the results obtained from the DB-FE method and from the FEIE method exhibit obvious differences. In Fig. 7a, the evolution lines of the configurational forces highly depend on the positions of the twin interfaces, while this position effect becomes relatively small in the results shown in Fig. 7b. As the FE-IE method yields more accurate values of the resultant contact forces, the evolution lines shown in Fig. 7b should be more reliable. Thus, here I would like to correct one assertion we proposed in Ref. [20]. That is, under mechanical loading conditions and with the free end boundary conditions, the positions of the twin interfaces in the sample have no much influences on the configurational force values. In Fig. 7, the dashed horizontal line represents the blocking force D. By comparing the values of D and T k and through the similar analyses as those given in Ref. [20], one can determine the procedure of twin interface movements under the mechanical load tA. Correspondingly, the global mechanical response of the NiMnGa sample can be predicted. In Fig. 8a and b, we plot the stress-strain curves of NiMnGa sample predicted by using the DB-FE method and the FE-IE method, respectively. For the purpose of comparison, the experimental results reported in Fig. 2 of Chernenko et al. [5] are also plotted. It can be seen that response of the NiMnGa sample can

312

J. Wang / Acta Materialia 105 (2016) 306e316

Fig. 6. The resultant contact forces on the twin interfaces under the external compressive stress sA ¼ 3 MPa: (a) and (b) the axial and lateral force components obtained from the DB-FE method; (c) and (d) the axial and lateral force components obtained from the FE-IE method.

Fig. 7. Evolutions of the configurational forces on the different twin pairs accompanying the variation of the external compressive stress sA: (a) results obtained from the DB-FE method; (b) results obtained from the FE-IE method.

be predicted at a quantitative level by the current model. From the theoretical analyses, it can be further revealed that for the mechanical loading pattern considered in this section (i.e., axial compressive load and free end boundary conditions), the configurational force T k only depends on the magnitude of tA and the area of the sample's cross section, which does not depend on the geometrical shape of the cross section. To verify this issue, we consider four NiMnGa samples with different cross-sectional shapes: (a) square with dimensions 4.25  4.25 mm2; (b) rectangle with dimensions 6.01  3.005 mm2; (c) rectangle with dimensions 3.005  6.01 mm2; (d) circular area with radius 2.398 mm. The cross-sectional areas of these four samples are

equal. The lengths of these samples are set to be 20 mm. These four samples are also meshed by using tetrahedral elements such that several inclined subregions are generated in the samples (cf. Fig. 9). The interfaces between the subregions are some parallel flat planes with the normal vector N given in (5)2. Similar as before, only the 17 interfaces labeled as I8, …,I0, …,I8 are considered in each sample (cf. Fig. 9). These four samples are subjected to the same mechanical load pattern as that shown in Fig. 5b. For the purpose of comparison, it is also assumed that only two twin interfaces exist in each sample and they are arranged symmetrically at the positions {Ii,Ii} (i ¼ 1, …,8). By using the FE-IE method, the configurational forces on the twin

J. Wang / Acta Materialia 105 (2016) 306e316

313

Fig. 8. The predicted stress-strain response of the NiMnGa sample by using (a) the DB-FE method and (b) the FE-IE method (solid curves: the model predictions; dashed curves: the experimental results taken from Fig. 2 in Ref. [5]).

Fig. 9. The four NiMnGa samples with the different cross-sectional shapes: (a) square with dimensions 4.25  4.25 mm2; (b) rectangle with dimensions 6.01  3.005 mm2; (c) rectangle with dimensions 3.005  6.01 mm2; (d) circular area with radius 2.398 mm.

interfaces can be calculated and the mechanical responses of these four samples can be predicted. Accompanying the variations of the external loads, the evolution lines of the configurational forces are shown in Fig. 10. From Fig. 10, it can be seen that the configurational forces in the four samples have almost the same values. Thus, it is verified that the cross-sectional shapes of NiMnGa samples have no much influence on the configurational force values. Furthermore, it has been found that the four samples exhibit the same global stress-strain response as that shown in Fig. 8b. Remark: It has been noticed for a long time that the interface element methods may result in undesirable spurious oscillations of the stress field on the interfaces [33,34,24,35]. This problem has been investigated systematically during the past years. In Schellekens and Borst [33] and Vignollet and Borst [35], this oscillation phenomenon is attributed to the high dummy stiffness of the interface elements in combination with a Gauss numerical integration scheme. It has also been found that the performance of interface elements can be improved by using a nodal lumping

scheme or NewtoneCotes integration schemes. In the current work, the oscillations of the contact stress can also be found on the twin interfaces. However, to predict the mechanical behavior of a NiMnGa sample, only the resultant contact forces or the configurational forces on the twin interfaces need to be considered. Through systematic investigations, it has been found that the numerical values of the resultant contact forces and the configurational forces obtained through the interface element method are stable, which do not depend on the shapes of the sample or the mesh formats.

5. Discussions and future tasks In this paper, only the mechanical response of NiMnGa samples is considered. To further predict the magneto-mechanical behaviors of NiMnGa samples, it is necessary to calculate the accurate configurational forces on the twin interfaces under combined magnetic and mechanical loading conditions. In Ref. [19], the

314

J. Wang / Acta Materialia 105 (2016) 306e316

Fig. 10. Evolutions of the configurational forces on the different twin pairs accompanying the variation of the external load sA: (a) square with dimensions 4.25  4.25 mm2; (b) rectangle with dimensions 6.01  3.005 mm2; (c) rectangle with dimensions 3.005  6.01 mm2; (d) circular area with radius 2.398 mm.

expression of the configurational force under combined loading conditions has been derived (cf. Eqns. (39) and (49) in Ref. [19]). From this expression, it can be found that besides the deformation gradient and the contact stress, the configurational force also depends on some quantities related to the magnetic field, including the demagnetization field, the effective magnetization vector and some internal variables representing the magnetic domain structures. To obtain accurate configurational force values, these quantities also need to be evaluated precisely. In Ref. [20], an iterative numerical algorithm has been proposed to solve the coupled magnetic and mechanical governing system derived for MSMA samples. The procedure of the algorithm is mainly composed of the following three steps:  Step one. With the fixed twin interface positions and the initial displacement field values, the magnetic field equation is solved to determine the updated scalar magnetic potential in the whole space, based on which the demagnetization field and the internal variables related to the magnetic domain structures can be calculated. The effective magnetization in the sample can also be determined by using the updated internal variables.  Step two. With the fixed twin interface positions and the updated magnetic related quantities, the mechanical equilibrium equation is solved to calculate the updated displacement field in the sample, based on which the deformation gradient and the stress in the sample can be obtained.  Step three. The updated values of the magnetic and mechanical quantities are used to calculate the configurational forces on the twin interfaces. According to the twin interface movement criteria, the positions of the twin interfaces are adjusted. Through several iteration steps, the numerical solutions to the coupled governing system can be obtained. It is clear that the IE method introduced in this paper can be incorporated in the second

step to obtain accurate contact stresses on the twin interfaces. Similarly, in our future work, the IE method will be adopted in the first step to derive the accurate values of the magnetically related quantities on the twin interfaces. For example, through some comparisons, it can be found that the relationship between the scalar magnetic potential and the demagnetization field is analog to that of the displacement and the stress tensor. Besides that, the internal variables and the effective magnetization vector appearing in the configurational force directly depend on the demagnetization field. Thus, if the IE method is adopted in the first step, the demagnetization field can be chosen as the primary unknown for the nodes on the twin interfaces. The connection equations for the demagnetization field across the twin interfaces also need to be proposed. By using the FE-IE method, the magnetic field equation can be solved numerically. The demagnetization field on the twin interfaces, together with the scalar magnetic potential in the sample, can be obtained. Based on the numerical results obtained in the first two steps, the accurate configurational forces on the twin interfaces can be calculated in the third step. Then, the positions of the twin interfaces in the sample can be determined by using the twin interface movement criteria. It should be noted that in the model proposed in Ref. [19], we did not consider the detailed arrangements of the magnetic domain structures at the microscopic level, e.g., motion of each magnetic domain wall or rotation of each local magnetization vector. The effects of the magnetic domain structures on the global response of the NiMnGa sample are considered in an average sense. That is, some internal variables are adopted to represent the volume fractions of the different magnetic domains and the average rotation angle of the local magnetization vector. With these internal variables, the constitutive form of the effective magnetization vector in the sample is proposed, which is further used to evaluate the magnetic response of the sample. This treatment can simplify the works in theoretical modeling and numerical simulations. For

J. Wang / Acta Materialia 105 (2016) 306e316

example, if the IE method is incorporated in the first step of the above algorithm, once the accurate values of the demagnetization field on the twin interfaces are obtained, the internal variables and the effective magnetization vector on the twin interfaces can be calculated by using the constitutive equations. There is no need to further consider the evolutions of the magnetic domain structures at the microscopic level. Besides the accurate evaluations of the configurational forces under combined magnetic and mechanical loading conditions, the current model still has some other open problems that need to be addressed in our future work. For example, in aspect of theoretical modeling, the current model cannot predict the nucleations of new twin interfaces in NiMnGa samples. The twin interface movement criteria (8) can only be applied when the twin interfaces already exist in the sample. Thus, to simulate the variant reorientation process, we need to assign the initial positions of the twin interfaces in the sample. Due to this treatment, some key experimental features accompanying the nucleation process (e.g., the nucleation stress values, the preferred sites of the nucleated twin interfaces) cannot be captured. It is hoped that the twin interface nucleation criterion can be derived by extending the current modeling approach. Ideally, this criterion should be able to predict the number and sites of the nucleated twin interfaces as well as the corresponding nucleation stresses. In aspect of numerical simulations, currently we always adopt a specific mesh such that the NiMnGa sample is averagely divided into several subregions. The interfaces between the subregions provide the possible (discrete) positions for the twin interfaces. By adopting this scheme, the computational consumption can be significantly reduced. However, the response of the NiMnGa sample predicted by using this scheme is usually discontinuous, which will influence the accuracy of the simulation results. Besides that, in this specific mesh, the elements in the different sample regions (near or far from the twin interfaces) have almost the same size, which is not suitable to capture the sudden changes of the physical quantities across the twin interfaces. To overcome these shortcomings, it is necessary to explore the ‘re-meshing’ techniques in our future work. To improve the numerical efficiency, finer mesh will be generated in the neighboring regions of the twin interfaces and coarser mesh will be generated in the other regions. Accompanying the movements of the twin interfaces, the sample region will be re-meshed. By virtue of the re-meshing techniques, it is hoped that the accuracy of the model predictions can be improved and the existence of multiple twin interfaces in the sample can be simulated. Our current model is proposed at the continuum level and through the thermodynamic approach. In this model, the twin boundaries are viewed as some sharp interfaces with given orientations, where the lattice structures and the magnetic domain structures undergo abrupt changes. The twin interface movement criteria is automatically integrated in the governing system of the model. If the governing system is solved numerically, the distributions of the related physical quantities (e.g., displacement, stress, magnetization) as well as the arrangements of the twin interfaces in the sample can be determined. Then, the global magnetomechanical behavior of the sample can be simulated. In the literature, another approach has been proposed to simulate the variant reorientations in MSMAs, which is called phase field simulations [36e41]. In this approach, some variant and magnetization order parameters are adopted to represent the microstructure state. The twin boundaries as well as the magnetic domain walls are modeled as some diffusive interfaces where the order parameters undergo continuous but rapid changes. By using the GinzburgeLandau equations and the Landau-Lifshitz-Gilbert equations, the evolutions of the twin boundaries and the magnetic domain walls can be

315

simulated, based on which the (local) magneto-mechanical response of MSMAs can be predicted. In our opinion, the phase field models and the thermodynamic ‘sharp interface’ model have their own features. Thus, they should be properly chosen to fulfill different modeling objects. For example, our thermodynamic model is efficient in simulating the global magneto-mechanical behavior of a MSMA sample with practical dimensions and subject to practical loading conditions. To our knowledge, the overall sample size that can be simulated by phase field models is restricted to a range of microns due to computational complexities. Besides that, the phase field models usually adopt the periodic boundary conditions but not the practical boundary conditions on the sample edges. On the other hand, the detailed arrangements of the magnetic domain structures are not considered in our thermodynamic model, which are just represented by some internal variables. Thus, our thermodynamic model cannot provide detailed descriptions on the microscopic features of the material during the variant reorientation process (e.g., the interactions between the magnetic domain walls and the twin interfaces, the complex magnetic domain structures near the twin interfaces). In phase field models, the microstructure of the martensite twins and the magnetic domains can be simulated accurately. The texture effect of multiple variants and the interactions between the magnetic domain walls and twin boundaries can also be studied. In these aspects, the phase field approach demonstrates superiorities over our thermodynamic model. 6. Conclusions In this paper, an interface element method [24] was integrated into the conventional DB-FE method to evaluate the configurational forces on the twin interfaces in NiMnGa samples. One key feature of this interface element method is that the internal contact stress on the twin interfaces is taken as one of the primary unknowns. By solving the problem, the internal contact stresses on the twin interfaces, together with the nodal displacements in the sample, can be directly obtained. Based on the obtained contact stresses, the resultant contact forces and the configurational forces on the twin interfaces can be calculated. These results were further used to predict the global mechanical behavior of NiMnGa samples. To verify the effect of the FE-IE method, a NiMnGa sample with rectangular parallelepiped shape was considered, which was subjected to the compressive stresses applied on the two ends. Free end boundary conditions were proposed to remove other constraint stresses in the sample. Through some analyses, one can derive the theoretical values of the resultant contact forces and the configurational forces on the twin interfaces. On the other hand, this problem was solved numerically by using the conventional DBFE method and the FE-IE method. By comparing the numerical results and the theoretical results, it was found that the results obtained from the DB-FE method exhibit large numerical errors, while the FE-IE method yields accurate solutions. Furthermore, it was verified that the cross-sectional shape of NiMnGa samples has no much influence on the configurational force values. The current work proposed an efficient approach to improve the evaluations of the configurational forces in NiMnGa samples, based on which the global mechanical behaviors of NiMnGa samples can be predicted more accurately. In our future work, the interface element method will be further used to evaluate the configurational forces when the NiMnGa samples are subjected to combined magnetic and mechanical loading conditions. The obtained results will be important for the practical design and application of NiMnGa components. Besides that, our current model still have some other open problems (e.g., modeling of the twin interface nucleation processes, improvement of the numerical accuracy by

316

J. Wang / Acta Materialia 105 (2016) 306e316

adopting the re-meshing techniques). These open problems also need to be addressed in out future work. Acknowledgments The work described in this paper is supported by Fundamental Research Funds for the Central Universities from SCUT (Project No.: D215215w, D2155500) and the Guangdong Natural Science Funds for Distinguished Young Scholar (Project No.: 2015A030306009). References [1] P.J. Webster, K.R.A. Ziebeck, S.L. Town, M.S. Peak, Magnetic order and phase transformation in Ni2MnGa, Philos. Mag. B 49 (1984) 295e310. [2] K. Ullakko, J.K. Huang, C. Kantner, R.C. O'Handley, V.V. Kokorin, Large magnetic-field-induced strains in Ni2MnGa single crystals, Appl. Phys. Lett. 69 (1996) 1966e1968. [3] R. Tickle, R.D. James, Magnetic and magnetomechanical properties of Ni2MnGa, J. Magn. Magn. Mater. 195 (1999) 627e638. [4] O. Heczko, A. Sozinov, K. Ullakko, Giant field-induced reversible strain in magnetic shape memory NiMnGa alloy, IEEE Trans. Magn. 36 (2000) 3266e3268. [5] V.A. Chernenko, V.A. L'vov, P. Mullner, G. Kostorz, T. Takagi, Magnetic fieldinduced superelasticity of ferromagnetic thermoelastic martensites: experiments and modeling, Phys. Rev. B 69 (2004) 134410. [6] L. Straka, V. Nova0 k, M. Landa, O. Heczko, Acoustic emission of Ni-Mn-Ga magnetic shape memory alloy in different straining modes, Mater. Sci. Eng. A 374 (2004) 263e269. [7] H.E. Karaca, I. Karaman, B. Basaran, Y.I. Chumlyakov, H.J. Maier, Magnetic field and stress induced martensite reorientation in NiMnGa ferromagnetic shape memory alloy single crystals, Acta Mater. 54 (2006) 233e245. [8] L. Straka, Magnetic and Magneto-mechanical Properties of Ni-Mn-Ga Magnetic Shape Memory Alloys (Doctoral Dissertation), Helsinki University of Technology, 2007. [9] Y.W. Lai, Magnetic Microstructure and Actuation Dynamics of NiMnGa Magnetic Shape Memory Materials (Ph.D. thesis), Technical University of Dresden, 2009. [10] N. Scheerbaum, Y.W. Lai, T. Leisegang, M. Thomas, J. Liu, K. Khlopkov, J. McCord, S. Fahler, R. Trager, D.C. Meyer, L. Schultz, O. Gutfleisch, Constraint dependent twin variant distribution in Ni2MnGa single crystal, polycrystals and thin film: an EBSD study, Acta Mater. 58 (2010) 4629e4638. [11] X. Chen, Y.J. He, Z. Moumni, Twin boundary motion in NiMnGa single crystals under biaxial compression, Mater. Lett. 90 (2013) 72e75. [12] L. Hirsinger, C. Lexcellent, Modelling detwinning of martensite platelets under magnetic and (or) stress actions in Ni-Mn-Ga alloys, J. Magn. Magn. Mater. 254e255 (2003) 275e277. [13] B. Kiefer, D.C. Lagoudas, Magnetic field induced martensitic variant reorientation in magnetic shape memory alloys, in: Philosophical Magazine Special Issue: Recent Advances in Theoretical Mechanics, Honor of SES 2003 A.C. Eringen Medalist G.a. Maugin, vol. 85, 2005, pp. 4289e4329. [14] L. Suorsa, Performance and Modeling of Magnetic Shape Memory Actuators and Sensors, Doctoral Dissertation, Helsinki University of Technology, 2005. [15] N.N. Sarawate, M.J. Dapino, Magnetomechanical characterization and unified energy model for the quasistatic behavior of ferromagnetic shape memory NiMn-Ga, Smart Mater. Struct. 19 (2010) 035001. [16] J. Wang, P. Steinmann, A variational approach towards the modeling of magnetic field-induced strains in magnetic shape memory alloys,, J. Mech. Phys. Solids 60 (2012) 1179e1200. [17] X. Chen, Z. Moumni, Y.J. He, W.H. Zhang, A three-dimensional model of magneto-mechanical behaviors of martensite reorientation in ferromagnetic

shape memory alloys, J. Mech. Phys. Solids 64 (2014) 249e286. [18] K. Halder, D.C. Lagoudas, I. Karaman, Magnetic field-induced martensitic phase transformation in magnetic shape memory alloys: modeling and experiments, J. Mech. Phys. Solids 69 (2014) 33e66. [19] J. Wang, P. Steinmann, On the modeling of equilibrium twin interfaces in a single-crystalline magnetic shape memory alloy sample. I: theoretical formulation, Contin. Mech. Thermodyn. 26 (2014) 563e592. [20] J. Wang, P. Steinmann, On the modeling of equilibrium twin interfaces in a single-crystalline magnetic shape memory alloy sample. II: numerical algorithm, Contin. Mech. Thermodyn. (2015), http://dx.doi.org/10.1007/s00161014-0403-4. [21] A. Gens, I. Carol, E.E. Alonso, An interface element formulation for the analysis of soil-reinforcement interaction, Comput. Geotech. 7 (1988) 133e151. [22] J.C.J. Schellekens, R. de Borst, The application of interface elements and enriched or rate-dependent continua to micro-mechanical analyses of fracture in composites, Comput. Mech. 14 (1994) 68e83. [23] G. Alfano, M.A. Crisfield, Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, Int. J. Numer. Meth. Eng. 50 (2001) 1701e1736. [24] X.Y. Lei, Contact friction analysis with a simple interface element, Comput. Methods Appl. Mech. Eng. 190 (2001) 1955e1965. [25] A.L.G.A. Coutinho, M.A.D. Martins, R.M. Sydenstricker, J.L.D. Alves, L. Landau, Simple zero thickness kinematically consistent interface elements, Comput. Geotech. 30 (2003) 347e374. [26] A. Shirazi-Adl, An interface continuous stress penalty formulation for the finite element analysis of composite media, Comput. Struct. 33 (1989) 951e956. [27] K.N. Song, A procedure for the improved continuous stress field of composite media, KSME Int. J. 12 (1998) 871e880. [28] X.W. Gao, K. Yang, Interface integral BEM for solving multi-medium elasticity problems, Comput. Methods Appl. Mech. Eng. 198 (2009) 1429e1436. [29] K. Yang, W.Z. Feng, X.W. Gao, A new approach for computing hyper-singular interface stresses in IIBEM for solving multi-medium elasticity problems, Comput. Methods Appl. Mech. Eng. 287 (2015) 54e68. [30] P. Steinmann, Views on multiplicative elastoplasticity and the continuum theory of dislocations, Int. J. Eng. Sci. 34 (1996) 1717e1735. [31] J.D. Eshelby, The elastic energy-momentum tensor,, J. Elast. 5 (1975) 321e335. [32] R. Abeyaratne, An admissibility condition for equilibrium shocks in finite elasticity, J. Elast. 13 (1983) 175e184. [33] J.C.J. Schellekens, R. de Borst, On the numerical integration of interface elements, Int. J. Numer. Methods. Eng. 36 (1993) 43e66. [34] V.N. Kaliakin, J. Li, Insight into deficiencies associated with commonly used zero-thickness interface elements, Comput. Geotech. 17 (1995) 225e252. [35] J. Vignollet, S. May, R. de Borst, On the numerical integration of isogeometric interface elements, Int. J. Numer. Methods. Eng. 102 (2015) 1733e1749. [36] C.M. Landis, A continuum thermodynamics formulation for micro-magnetomechanics with applications to ferromagnetic shape memory alloys,, J. Mech. Phys. Solids 56 (2008) 3059e3076. [37] P.P. Wu, X.Q. Ma, J.X. Zhang, L.Q. Chen, Phase-field simulations of stress-strain behaviour in ferromagnetic shape memory alloy NiMnGa, J. App. Phys. 104 (2008) 073906. [38] Y.M. Jin, Effects of twin boundary mobility on domain microstructure evolution in magnetic shape memory alloys: Phase field simulation, Appl. Phys. Lett. 94 (2009) 062508. [39] L.J. Li, C.H. Lei, Y.C. Shu, J.Y. Li, Phase-field simulation of magnetoelastic couplings in ferromagnetic shape memory alloys, Acta. Mater. 59 (2011) 2648e2655. [40] C. Mennerich, F. Wendler, M. Jainta, B. Nestler, A phase-field model for the magnetic shape memory effect,, Arch. Mech. 63 (2011) 549e571. [41] Q. Peng, Y.J. He, Z. Moumni, A phase-field model on the hysteretic magnetomechanical behaviors of ferromagnetic shape memory alloy, Acta. Mater. 88 (2015) 13e24.