Numerical study on butterfly wings around inclusion based on damage evolution and semi-analytical method

Butterfly wings are closely related to the premature failure of rolling element bearings. In this study, butterfly formation is investigated using the developed semi-analytical three-dimensional (3D) contact model incorporating inclusion and material property degradation. The 3D elastic field introduced by inhomogeneous inclusion is solved by using numerical approaches, which include the equivalent inclusion method (EIM) and the conjugate gradient method (CGM). The accumulation of fatigue damage surrounding inclusions is described using continuum damage mechanics. The coupling between the development of the damaged zone and the stress field is considered. The effects of the inclusion properties on the contact status and butterfly formation are discussed in detail. The model provides a potential method for quantifying material defects and fatigue behavior in terms of the deterioration of material properties.


Introduction
The occurrence of premature failures due to white etching cracks (WECs) in rolling element bearings is a primary limitation to improving the service performance and reliability of advanced equipment, such as wind turbines, high-speed trains, and tunnel boring machines. It has been reported that failures can occur in as little as 5%-20% of the desired L10 life of bearings in wind turbine gearboxes (WTG) [1,2]. WECs are closely related to butterflies, which are symbiotic microstructures that contain white etching areas (WEAs), cracks, and inclusions. As one of the most important initiators of subsurface cracks, butterflies are characterized by a central stress concentrator (material defects such as inclusions) and surrounding microstructural alterations [3], as shown in Fig. 1. WEAs appear in pairs as wings on either side of the inclusion at angles of    30 50 to the rolling direction [4,5]. The formation and morphology of butterflies are strongly linked to the local stress state. The shape, size, rigidity, and distribution of inclusions may change the stress/strain behavior, thus affecting the initiation and evolution of butterflies [6,7]. Following the pioneering observations of butterflies by Jones [8] and Styri [9], failure modes related to WECs and butterflies have been continuously investigated over the past decades. Evans et al. [10,11] successfully reproduced the failure caused by WECs on a roller tribometer and a FAG-FE8 test rig. The evidence they obtained through the mapping and reconstruction of entire independent WECs using a combination of focused ion beam (FIB) and serial sectioning tomography techniques strongly supports the hypothesis that butterflies can be introduced by inclusions. In field bearing tests, bearings removed from a WTG were examined using X-ray microtomography, atomic force microscopy (AFM), and ultrasonic mapping. The crack morphology, microstructure alterations, and inclusions provide further support for the correlation between inclusions and butterflies [5,6,[12][13][14]. The cracks are preferentially initiated at the central inclusion and propagated in the direction of over-rolling. The critical depth of the crack initiation has a strong correlation with the positions of the maximum shear stress. The factors influencing butterflies, such as kinematic sliding, lubricant composition, and material properties have been extensively investigated in benchtop tests using ball-on-rod or triple-disc rolling contact fatigue rigs [15][16][17][18]. The direction and magnitude of the sliding are highly related to the formation of butterflies. It has been reported that WECs and butterflies can be formed with synthetic commercial oils as well as base oils under both positive and negative sliding in various types of bearing steels [19,20]. Accordingly, butterflies, as a common type of rolling contact fatigue, should be further investigated numerically owing to their complicated stress/strain behavior. A contact model that incorporates inclusion and damage evolution is necessary to understand the fundamental damage mechanisms of butterflies. This is the primary requirement for a numerical simulation of butterflies that solves for the subsurface stress field surrounding the inclusion and describes the microstructural alteration accompanying damage evolution. In early attempts made by Salehizadeh and Saka [21], Melander [22], and Yamakawa et al. [23], the local stress/strain fields around the inclusion were calculated using the finite element method (FEM). Several different types of inclusions (manganese sulfide, alumina, and titanium nitride) under rolling contact conditions were investigated. Raje et al. [24] proposed a discrete element approach to compute the elastic stress in an inhomogeneous medium. Ueda et al. [25] analytically predicted the location and direction of the initial crack emanating from an inclusion. Following McDowell's pioneering work [26], Alley and Neu [27] developed a numerical model of crystal plasticity using FEM. This influential simulation framework laid the foundation necessary to determine the effect of the microstructure in rolling contact fatigue (RCF). Based on this framework, the microstructural properties, including anisotropy, inclusion, and phase state, have been extensively studied [28][29][30]. Cerullo and Tvergaard [31] obtained the stress field surrounding the inclusion and evaluated the contact fatigue using the DangVan multiaxial criterion. Moghaddam et al. [32,33] investigated the formation of butterflies based on a FEM model and continuum damage mechanics (CDM). Guan et al. [34] studied the crack nucleation and propagation induced by metallic carbides using the Voronoi finite element model. To capture the stress variation around the inclusion, a sufficiently fine meshing scheme is necessary for the numerical model. This would result in a long computation time for the FEM calculation, especially for 3D contact under cyclic loading.
The recent development of semi-analytical methods (SAMs) provides an alternative way to simulate butterfly formation. In semi-analytical models, the computation domain is discretized into cuboidal elements with uniform stress and strain. The total stress is given by the superposition of the influence coefficient (IC) response of each element solved by an analytical method. Jacq et al. [35] and Nelias et al. [36,37] developed a 3D elastoplastic contact model using SAM. The equivalent inclusion method (EIM) assumes that an inhomogeneous area can be simulated as an inclusion with an equivalent eigenstrain. Based on this idea, Zhou et al. [38] proposed a semi-analytical solution that describes the disturbance stress field caused by multiple inclusions. Liu et al. [39] provided explicit solutions of the elastic field caused by an arbitrarily-shaped inclusion based on Galerkin vectors. In recent years, this SAM framework incorporating EIM has been extensively applied to the contact problems of inhomogeneous materials, including lubrication [40], partial slip [41], distributed inhomogeneities [42], and viscoelastic contact [43]. The combination of SAM and discrete convolution and www.Springer.com/journal/40544 | Friction fast Fourier transformation (DC-FFT) algorithm laid the foundation for efficiently and accurately simulating the damage evolution around inclusions.
In this work, a coupled-damage model incorporating inhomogeneous inclusions and material property degradation is proposed to investigate butterfly formation. The elastic field introduced by the inclusion and deterioration of material properties are obtained using SAM and EIM. The fatigue damage accumulation caused by the concentrated stress is described by CDM. Furthermore, the effects of inclusion depth and orientation angle on butterfly formation are discussed in detail. The findings from this study provide further insights into the mechanisms underlying butterfly formation. Although this numerical model can be further improved, for example by including the thermal effect and crack propagation, it provides a new approach for quantifying the damage evolution around the inclusion.

Contact problem description
The numerical model consists of four major components, as illustrated in Fig. 2. (i) The pressure and traction are determined using the compatibility conditions of the contact surfaces. (ii) The inhomogeneous inclusions (for both non-metallic inclusions and microstructural alteration) are replaced by an imaginary equivalent homogeneous inclusion with eigenstrain. (iii) The elastic field incorporating the disturbance stress and displacement due to the eigenstrain is obtained. (iv) The damage accumulation and elastic modulus variation around the inclusions are updated to represent the microstructural alteration. The influences of lubrication and roughness are neglected because only a limited amount of stress variation is introduced by the lubrication condition on or near the surface [44,45]. The gradient of the hardness is not considered because the simulated material is through hardening bearing steel. The four components are described in the following sections.
The general contact problem is governed by the gap and load balance equations as Eqs. (1) Here, ( , ) h x y is the gap between the two contact surfaces; 0 h and m h denote the rigid gap and initial macro-geometry, respectively; and 3 u is the surface normal displacement, which contains the elastic deformation e 3 , ( ) u x y and eigen-displacement r 3 , ( ) u x y in this work. Equation (2) represents the equilibrium of the surface pressure , ( ) p x y and the normal load n W in the contact area . The boundary constraints are given as The surface pressure can be solved based on the widely applied conjugate gradient method (CGM) [46]. The surface shear traction is simplified as the Coulomb friction. The coefficient of friction is obtained from the results of rings-to-roller tests [44].

Equivalent inclusion method for inhomogeneous inclusion
The inclusions and surrounding microstructural alterations can be treated as inhomogeneous inclusions because of their different elastic moduli [33]. As shown in Fig. 3, a deformable body with elastic modulus ijkl C contains several subdomains , which are assumed to be perfectly bonded with the matrix and have different elastic moduli  ijkl C . The subdomains are, in general, subjected to the initial eigenstrain  p kl due to self-equilibration and the applied stress  0 ij caused by an external load. The core idea in EIM is that an imaginary inclusion with the same elastic properties as the matrix can be enriched to simulate the influence of an inhomogeneous inclusion.
This inclusion contains the initial eigenstrain  p kl and the equivalent eigenstrain  * kl , which represents the disturbance caused by the inhomogeneity. Therefore, according to Hooke's law, the stress inside the subdomains can be expressed as Equations (4a) and (4b) indicate that the key step in EIM is to solve the equivalent eigenstrain  * kl . Rewriting the total strain in Eq. (4a) as can also be written as The total stress can be determined as where  0 ij is the stress induced by the external load.
 p ij and  * ij are the eigenstresses caused by the initial eigenstrain and equivalent eigenstrain, respectively. Equation (6) can be substituted into Eq. (5) to give the following: Equation (7) can be written in matrix form as follows: In Eq. (8), the three unknowns are * σ , p σ , and * ε .
The explicit expressions for the elastic constants, C and  1 C , are listed in Appendix A. In general, the mapping relationship from the eigenstrain to the eigenstress can be assumed to take the form of Equation (8) can then be written as Once the operation Y is determined, the governing equation of the equivalent eigenstrain becomes solvable, because * ε is the only unknown. In short, Eq. (10) can be simplified to a linear equation, that is,  * ε F a, for which the conjugate gradient method (CGM) can be used to obtain the solution. More details on the CGM and the solution of Eq. (10) can be found in Appendix B. A key part of this procedure, the operation Y associated with the eigenstrain problem, is explained in the next section.

Stress field and surface displacement
As summarized by Mura [47], incompatibility of the eigenstrains results in the creation of eigenstress fields. The general solution of elastic fields related to eigenstrain has been investigated by Chiu [48], Jacq et al. [35], Nelias et al. [36], and Liu et al. [39] over the decades. A special form of this problem considering elastoplastic contact was also discussed by the author in Ref. [49]. In general, the eigenstrain can be treated as an excitation in the elastic field and the eigenstress as its response. Once the point-to-point response function has been derived, the total stress can be determined via the summation of every element. The elastic deformation due to eigenstrain can be expressed in terms of the Galerkin vectors, F [50]: where  is the Poisson's ratio, and  is the Lamé constant. The commas in the subscripts indicate the derivatives with respect to the indices following the commas, e.g., Once the deformations are determined using Eq. (11), the eigenstress inside and outside the inclusion area Ω can be derived using Hooke's law: The entire calculation domain can be meshed into   M N L cuboid elements, each of which has a uniform eigenstrain and stress. After discretization, the stresses are given as the summation of cuboidal contributions in terms of the influence coefficients: Here,    , , denote the site of the element of interest, and    , , are the element indices along the x-, y-, T are the influence coefficients denoting the mapping relation from the unit eigenstrain to the eigenstress at a specific element. The explicit forms of the influence coefficients can be found in the appendices of Refs. [39,51]. Thus, the discrete form of operation Y for obtaining the eigenstress is determined by Eq. (13).
The calculation of the eigen-displacement can also be arranged into the form of an influence coefficient: The explicit forms of the terms in Eq. (14) can be found in Ref. [52]. The calculation of the eigenstress and deformation involves three-dimensional discrete convolution, for which 3D-DC-FFT [53,54] can be applied to accelerate the computation. It should be noted that the eigenstress calculation may be performed many times depending on the iterations of the CGM required to solve the equivalent eigenstrain.

Damage evolution
Butterfly formation is closely related to microstructural alterations, which can be described by the degradation of material properties [33]. Continuum damage mechanics (CDM) provides a potential method for evaluating the damage evolution and material deterioration. The phase transition from martensite to ferrite during butterfly formation is assumed to be represented by the damage variable D. In CDM theory, the stress-strain behavior is modified by the damage as follows [55]: The damage is reduced to a scalar variable under the assumption of isotropy. The damage-coupled constitutive equation then leads to the equivalent elastic modulus Based on the results of nanoindentation tests, Moghaddam et al. [33] determined the elastic modulus of butterfly wings to be  * 0.9 ijkl ijkl C C . Herein, the same value of the equivalent elastic modulus is applied, and the critical damage value is set to  * 0.1 D accordingly. It is assumed that the element of interest is transformed and that no more damage occurs after the damage accumulation reaches the critical value. In the elastic | https://mc03.manuscriptcentral.com/friction contact state, the evolution of the damage variable D can be described by [55]: where  r is the reference shear stress and m denotes the stress component related to the Weibull slope.
The parameters for the widely used AISI 52100 bearing steel taken from published works [44,45] are  10.1 m ,   r 6,113 MPa. The butterfly is a typical shear stress-dominated RCF failure mode. Alternating shear stress is the primary driving factor for butterfly formation [2]. The effect of the mean stress cannot be neglected because the stress disturbance around the inclusion is asymmetrical. To incorporate the effects of both the mean stress and the stress amplitude during the rolling contact process, the equivalent stress is represented by [32]: where  alt is the amplitude of the alternating shear stress, and  m is the mean shear stress. The jump-incycle method is introduced to capture the cyclic loading state intrinsic to high-cycle fatigue. In this procedure, it is assumed that the stress range of each element changes very little over a specific contact cycle of N [56][57][58]. Within the given finite cycles, the damage increment is given by The contact solver is executed at every contact location. At the end of the rolling passage, the equivalent stress is extracted from the stress history, and the damage rate d /d D N is calculated. Then, the total damage accumulation of each element is determined as The local elastic modulus of the matrix is updated, and the heterogeneity element is embedded in the damage site. The next rolling passage is repeated using the updated elastic modulus. To reduce the computing cost, the damage threshold is set as  t 0.001 D . Once the total damage accumulation of the element reaches the critical value, the transformed elements are added to the damage location of the butterfly map.

Numerical simulation 3.1 Parameter determination
The simulation was performed based on the results of the contact fatigue tests reported in a previous study [44]. The tests were conducted on a three-rings-on-roller test rig, as shown in Fig. 4(a). The diameters of the roller and rings were 12 mm and 54 mm, respectively. The material of the rings and roller was AISI 52100 through hardening tempered martensitic steel. The hardness of the roller and rings were 60 and 63 HRC , respectively. The Poisson's ratio and elastic modulus of the material matrix were 210 GPa and 0.3, respectively. The applied load in the test was approximately 500 N , for which a maximum Hertzian pressure of  h 1.9 GPa p and contact half width of  0.164 mm b were maintained. The tests were performed under different slide-to-roll ratios and rolling speeds, and the roller was sectioned and inspected to identify the sub-surface-initiated cracks. Figures 4(b) and 4(c) show the results of one test in which the butterfly was identified. The coefficient of friction is shown in Fig. 4(b). The mean value throughout the test was   0.081 , and the fatigue life was 42 million cycles. Herein, the coefficient of friction in the simulation was set to 0.08.
To simulate the test, the contact between the roller and rings was simplified into a three-dimensional line contact in which a rigid roller contacted with a deformable plane containing an inhomogeneous inclusion. The load and geometry conditions were the same as those of the RCF test. The Poisson's ratio and elastic modulus of the inclusion were set as 389 GPa and 0.25, respectively, to be consistent with those of 2 3 Al O , the most common type of oxide inclusion in 52100 bearing steel [59]. A spherical inclusion with the radius of 0.2b was set at the depth of  d 0.5 z b in the baseline case. The computation domain of   6 2 b b b was discretized into a   193 65 33 grid. The movement of the roller from  2 x b to  2 x b was decomposed into 17 contact steps. The cyclic loading was applied in the contact procedure via the jump-in-cycle method. The material properties and damage accumulation were updated for every loading block. The number of loading blocks was determined based on N and the tested fatigue life.

Numerical verification
To verify the computing accuracy and efficiency of the developed model, a comparison between the SAM solution and FEM results was conducted using the commercial software ABAQUS version 6.14. To remove the impact of the free boundary, a two-dimensional (2D) FE model and a three-dimensional (3D) quarter FE model were compared with the SAM solutions to evaluate the computing accuracy and efficiency, respectively. The geometry conditions of the 2D model matched those of the actual contact. A 12 mm-diameter roller containing a rigid inclusion was loaded on the 54 mm-diameter ring, as shown in Fig. 5(a). The fine mesh area in the region around the inclusion was partitioned with a local seed size of 0.0025 mm, and the element type was the eight-node plane strain reduced integration element (CPE8R). The 3D model was simplified as the contact pair of an analytical rigid roller with an equivalent radius of 4.9 mm and a quarter deformable cuboid based on Hertzian contact theory, as shown in Fig. 5(b). The size of the deformable body was fixed as   10 6 10 b b b along the x-, y-, and z-directions. The local seed size around the inclusion was set to 0.005 mm, which is close to the SAM element size (b/32). The total number of elements in the 3D model was 401,680, which is approximately half the number of SAM elements (413,985). The simulation was performed on a platform with an Intel(R) Core™ i5-8500 processor @ 3.00 GHz (6 cores, 6 threads) and 16 GB memory.
The von Mises stress for a square and a circular inclusion are compared in Fig. 5(d) and Fig. 5(e), respectively. The stress distributions of the SAM solutions agree well with the FEM results. The maximum stress in Fig. 5(d) obtained by SAM and FEM are 1,637.5 MPa and 1,642 MPa , respectively. In Fig. 5(e), several stress peaks can be found on the edge of the inclusion. These peaks resulted from the non-smooth interface of the circle for the discrete grids. The peaks would not affect the accuracy of the current model because the stress and damage within the inclusion are not considered in the work. The computation time for the 3D model is shown in Fig. 5(c)

Stress distribution and damage evolution
The stress history (i.e., time sequence of stress) is considered to be a major driver of butterfly formation and premature contact fatigue [20,60]. To capture the variation in stress history, a rolling contact procedure was simulated, and the variation of the stress components was documented. Figure 6 shows the variations in the alternating shear stress  xz during a rolling contact cycle. The sign of the shear stress depends on the direction of the z-axis. In this work, the z-axis points towards the depth direction. The minimum and maximum shear stresses are marked in black across the contour plots. It can be readily seen that the extremum of shear stress appears near the inclusion as the contact point moves close to the center. The maximum shear stress in Figs. 6(b) and 6(c) is 432.28 and 432.18 MPa, respectively, and increases to   Fig 6(d). An obvious stress concentration near the inclusion can be observed when the highstress area passes over the inclusion. The stress below the center point, which is denoted by the black solid line, is plotted in the second line of Fig. 6. The red portion of the line represents the stress within the inclusion, which will not be presented in the remainder of this work because the damage accumulation in the inclusion was neglected. It can be readily seen that the site of the maximum alternating stress changes as the roller moves. When the contact point reaches the center (Fig. 6(c)), a small divergence of the shear stress below the center point can be observed. This value should be zero in the contact condition without the inclusion. The variation of the alternating shear stress during the rolling contact process is asymmetric owing to the inclusion. This indicates that a nonzero mean shear stress exists around the inclusion.
A further analysis of the shear stress history around the inclusions is presented in Fig. 7. It should be noted that the inclusion and surrounding areas present a non-smooth interface owing to the intrinsic discrete characteristics of the SAM. Figure 7(a) illustrates the composite shear stress in a rolling contact cycle. Two distinguishable regions can be found in the contours of the equivalent shear stress. These high-stress regions are diagonally distributed at the top left and bottom right of the inclusion. This distribution is consistent with the experimental observation that the growth of butterfly wings occurs at a specific angle to the rolling direction. Figure 7 58 MPa, respectively. The variation of the equivalent stress components, that is, the alternating amplitude and mean shear stress, with the rotational angle  are shown in Fig. 7(c). It can be seen that the differences in the alternating amplitudes at different positions are limited, while the mean shear stress at the regions of the two wing regions is higher than that at the other regions. This implies that the mean shear stress is highly significant in butterfly initiation. It is hypothesized that the mean shear stress is related to surface traction, which introduces the disturbance of the shear stress around the inclusion.
The formation of butterfly wings during cyclic loading is shown in Fig. 8. The damaged zone, where microstructural alteration occurs, is represented by the small blue circles around the inclusion. The chronological order of butterfly formation illustrates the trend that the damage accumulates along a specific direction at an angle of approximately 45°. At the early stages, the left upper wing (which is closer to the surface) is more developed, as shown in Figs. 8(b) and 8(c). The two sides of the wings gradually propagate and become equally large as the damage accumulates. The emergence of the secondary wing is shown in Fig. 8(f). This process is consistent with the results reported by Moghaddam et al. [32].
The damaged elements are treated as inhomogeneous inclusions because of the degradation of the material properties (i.e., the reduction of elastic modulus). This deterioration process is accompanied by a variation in the stress field, as shown in Fig. 9. The deteriorated elements introduce stress disturbances in the surrounding area and cause further damage. The area of the lower wing is wider than that of the upper wing. One possible  reason for this is the larger area that experiences higher shear stress at the depth of the lower wing. Figure 10 shows the distribution of the damage around the inclusion on different XZ planes at  0 y , 0.0625b, 0.125b, 0.1825b and 0.25b. These figures present the same view as the serial sectioning of the tested samples. It can be readily seen that the damaged area deceases obviously with increasing distance from the mid-plane (  0 y ). However, damage still occurs outside of the planes containing the inclusion.  This suggests that butterfly wings without inclusions might be observed in the sectioned samples. Especially for spherical inclusions, the wings away from the middle plane may grow beyond the extent of the inclusion if the loading cycle is high enough.

Effect of inclusion depth
The inclusion depth is a significant factor influencing butterfly formation. The stress concentration site is highly related to the inclusion depth. In this section, the depth is set as   d 0.3 0.7 z b b to examine its effect on butterfly formation. Figure 11 shows the composite equivalent stresses at different depths at the first loading cycle. For the depth of

Effect of inclusion orientation angle
Irregular globular inclusions are a common type of material defects in engineering practice. In this section, an ellipsoid inclusion with the axial ratio of 2:1:1 at the depth of  0.5 z b is assumed to investigate the influence of the inclusion orientation angle on butterfly formation. As shown in Fig. 14, the inclusion orientation angle  is defined as the angle between the major axis of the ellipsoid and the rolling direction. Figure 15 illustrates butterfly formation under different  . It can be observed that the growth of the butterfly traces the same diagonal line for different direction angles. The upper wing is more developed at     0 60 , at which the long axis is closer to being parallel to the rolling direction. The initiation of the secondary pair of wings can be observed in these cases. At    120 and  150 , butterfly wings grow regularly around both sides of the long axis. Figure 16 further illustrates the number of damaged elements representing the propagation of butterfly wings. It can be readily seen that the damaged area is larger at acute direction angles compared to obtuse direction angles. One possible reason is the correspondence between the depth of the damaged regions and the depths at which the shear stress is higher. Another hypothesis is that the overall rolling direction might affect the stress history at different direction angles, which in turn changes the damage evolution process. This might explain the experimental observation of Bruce et al. [5] that longer cracks tend to be initiated from inclusions parallel to the raceway surface. Figure 17 shows the stress fields for different orientation angles after butterfly formation. The areas of stress concentration area correspond well with the butterfly wings. At    30 and  60 , discontinuous high stress distributions can be found in the area of lower wing. This is related to the scattered damaged elements observed in Figs. 15(b) and 15(c). At    120 and  150 , the stress distribution and damaged area are more extensive and form a new ellipsoid-like shape. The newly formed ellipsoid zone has a similar stress distribution to the original inclusion. It can be presumed that this similar shape would promote the further propagation of the damaged area, as per the butterfly wings (shown in Fig. 18) documented by Ganti et al. [61].

Conclusions
A semi-analytical 3D contact model incorporating inclusion and material property degradation was proposed to simulate butterfly formation. In the model, the elastic field caused by inhomogeneous and material deterioration is solved for by using a series of numerical methods. The accumulation of fatigue damage surrounding the inclusions is described by continuum damage mechanics. The coupling between the development of the damaged zone and the stress field is considered. The effects of the inclusion properties on the contact status and butterfly formation were discussed. The major conclusions are summarized as follows: (1) The distribution of alternating shear stress is asymmetric around the inclusion. The stress analysis shows that the mean shear stress plays a significant role in butterfly formation.
(2) The depth of the inclusion directly affects the initiation of the butterfly by changing the site of high stress concentration. The coupling between the damage evolution and the stress variation around the inclusion is a major manifestation of butterfly formation.
(3) For ellipsoidal inclusions, the orientation angle  significantly affects the formation of butterflies.
The damaged area is larger at     0 60 than at      90 150 . This work provides a potential method for quantifying material defects and fatigue behavior in terms of the deterioration of material properties. In addition, multiple inclusions can be considered in the proposed model. However, some influencing factors, such as the thermal effect and imperfect bonding at the inclusion-matrix interface, still require further study.