A Self-Learning Data-Driven Development of Failure Criteria of Unknown Anisotropic Ductile Materials with Deep Learning Neural Network

: This paper ﬁ rst proposes a new self-learning data-driven methodology that can develop the failure criteria of unknown anisotropic ductile materials from the minimal number of experimental tests. Establishing failure criteria of anisotropic ductile materials requires time-consuming tests and manual data evaluation. The proposed method can overcome such practical challenges. The methodology is formalized by combining four ideas: 1) The deep learning neural network (DLNN)-based material constitutive model, 2) Self-learning inverse ﬁ nite element (SELIFE) simulation, 3) Algorithmic identi ﬁ cation of failure points from the self-learned stress-strain curves and 4) Derivation of the failure criteria through symbolic regression of the genetic programming. Stress update and the algorithmic tangent operator were formulated in terms of DLNN parameters for nonlinear ﬁ nite element analysis. Then, the SELIFE simulation algorithm gradually makes the DLNN model learn highly complex multi-axial stress and strain relationships, being guided by the experimental boundary measurements. Following the failure point identi ﬁ cation, a self-learning data-driven failure criteria are eventually developed with the help of a reliable symbolic regression algorithm. The methodology and the self-learning data-driven failure criteria were veri ﬁ ed by comparing with a reference failure criteria and simulating with different materials orientations, respectively. generating suf ﬁ cient stress-strain training data from global experimental measurements. Symbolic regression technique such as genetic programming is useful for generating mathematical equations from


Introduction
Data-driven computational mechanics is one of the branches where the underlying laws such as boundary constraints, material constitutive law, or energy conservation law are replaced or collaborated with the experimental data in non-conventional schemes. Among them, the material constitutive law that has long been based on the traditional empiricism is relatively more influenced by experimental noises/ errors or uncertainties than other physics-based law associated with boundary value problems. Owing to advances in measurement science and technologies at multiple scales, traditional empiricism in the material constitutive modeling has started transferring to a data-driven paradigm. Not only does the inundation of data truly make it possible to advance computational mechanics but also the discovery of new physics-based material laws is enabled by emerging applications of data science, material informatics, and artificial intelligence (AI). The true data-driven (i.e., knowledge-based approach) approach was originally proposed by Ghaboussi and his co-workers [1]. They replaced the material constitutive model in the finite element model with an artificial neural network (ANN) form. The resurrection of artificial intelligence through recent deep learning advances has exploded in numerous applications to computational mechanics and material constitutive modeling [2][3][4][5]. In particular, Oishi et al. [3] proposed a new method of numerical quadrature enhanced by deep learning for the FEM stiffness matrices. Many researchers have applied ANN to various research fields such as cyclic plasticity [6,7], fiber-reinforced polymeric composites [8][9][10], rubber materials [11,12], rate-dependent materials [13,14], other composite materials [15,16], finite deformation with hyperelastic material [17] and traction-separation laws by reinforcement learning technique [18]. Recently, Liu et al. [19] also studied a sophisticated 3D network architecture of the deep material network for data-driven multiscale mechanics. Computational discretization by numerical methods such as the finite element method is required to apply ANNs in their problems [4,5].
Beyond the material constitutive modeling, more intrusive approaches of applying deep neural networks in diverse computational physics problems have also been researched recently [20][21][22][23][24][25]. Based on the idea of solving partial differential equations by ANN [26], Raisssi et al. [20] proposed physics-informed deep learning which models physics laws at random collocation points on the boundary and the domain. This method was improved by an adaptive collocation method for second-order boundary value problems (BVP) [24]. However, the current constraint in AI-based data-driven computational mechanics lies in a scarcity of data. Therefore, there is a pressing need for data-efficient training and modeling approach for solving such mathematically formulated physical problems with the help of deep learning neural networks.
The data-driven approach also encompasses the recent stress-strain data-driven computational mechanics approach originally proposed by Kirchodoerfer et al. [27]. This approach is also a model-free numerical method in the regime of computational mechanics. Their novelty lies in the fact that this approach does not require an explicit form of constitutive equations in their formulation. This approach has been further extended to diverse problems [27][28][29][30][31]. Eggersmann et al. [29] extended Kirchdoerfer and Oritz's data-driven method to inelasticity. Data from experiment tests have been mainly used for parameter identification [32] or model updating within the empiricism regime rather than replacing those laws or constraints in the boundary value problems. Leygue et al. [33] proposed a methodology called data-driven identification (DDI) that can build materials response data from full-field digital image correlation (DIC) data. They utilized a framework of data-driven computational mechanics (DDCM) by Kirchdoerfer et al. [27]. Stainier et al. [31] combined DDI with DDCM for bypassing the empiricism of material modeling. Further active studies on data-driven mechanics are undergoing by numerous researchers [34][35][36][37].
In this paper, we focus more on the aforementioned AI-based data-driven approach where the ANN model entirely replaced material constitutive law in the computational mechanics form [38,39]. ANN material models can predict the nonlinear multi-axial stress-strain relationships both under monotonic and cyclic loading [39]. Intrusive implementation technique of the ANN material constitutive model within finite element analysis codes is available [40,41]. However, one of the challenges with the ANN material constitutive model is the availability of comprehensive stress-strain training data from experiments, which is formidable in usual material tests. For tackling such challenges of ANN models, Ghaboussi et al. [1] proposed an online training methodology called an auto-progressive training whereby ANN material constitutive models are automatically trained in the course of nonlinear finite element analyses subjected to experimental boundary reaction forces and displacements. This powerful method has the advantage of generating sufficient stress-strain training data from global experimental measurements. Symbolic regression technique such as genetic programming is useful for generating mathematical equations from material response extracted by the auto-progressive training algorithm. The online auto-progressive training of ANN models is distinct from straightforward training of ANN models with stress-strain data. On the other hand, the deep learning neural network (DLNN) has a deeper hidden network of the perceptron. The backpropagation training algorithm and conventional activation functions such as hyperbolic tangent or sigmoid functions are not suitable for training the DLNN. Although the DLNN is advantageous in learning massive data, few research on its application to the material constitutive model, nonlinear finite element analysis, and the data-driven modeling have been researched to the best knowledge of authors.
In this paper, we propose a new self-learning data-driven modeling methodology that can discover the failure criteria of uncharacterized materials. Its novelty lies in the first DLNN-based material constitutive model formulated for nonlinear finite element analysis, on-line hybrid numerical/experimental training, data processing, subsequent symbolic regression, and their novel integration. The proposed method is verified through developing a new anisotropic initial failure criterion with a known reference failure criteria and synthetic simulated test data. Significance of the proposed self-learning data-driven modeling idea lies in its new unified and efficient ability to develop failure criteria of new materials such as process-dependent 3D printed materials from minimum experimental test data. This paper is organized as follows. Section 2 explains data-driven finite element analysis with DLNN material constitutive model. Section 3 describes the on-line hybrid numerical/experimental training called the self-learning inverse finite element simulation (SELIFE). It is followed by data processing ideas for identifying failure points and subsequent symbolic regression in Section 4. Section 5 presents a demonstration of the proposed methodology and verification of results with simulated tests. Finally, the conclusions are made in Section 6.

Data-Driven Finite Element Analysis with DLNN Material Constitutive Model
The data-driven FE simulation means the nonlinear FE simulations with the DLNN-based material constitutive model, which can be self-organized by experimental measurements. The DLNN version material constitutive model was first attempted in this paper with the derivation of the algorithmic tangent operator to be explained in the following.

DLNN Material Constitutive Model
Data features for input & output layers and DLNN material constitutive model in the data-driven nonlinear finite element model are shown in Fig. 1. Data features for the input and output layers are defined in terms of stresses and strains. Input features include stress-strain pairs and history-dependent internal variables. The internal variable means the sum of the total energy density at previous increment and incremental strain energy as Figure 1: A flow of the data-driven FEA with DLNN material constitutive model where, nÀ1 r i and nÀ1 e i are the stress and strain vector at (n − 1)-th increment, respectively; D n e i is the incremental strain, and nÀ1 & i is the internal variable. Especially, these variables transform one-to-many strain-to-stress mapping to one-to-one mapping between strain and stress values, resulting in a robust learning capability of hysteretic behavior of materials [39]. The sizes of the input and output features depend on the dimensionality of the finite elements. For example, in the case of the plane stress condition, twelve input nodes (i.e., n e 11 ; n e 22 ; n e 12 ; nÀ1 e 11 ; nÀ1 e 22 ; nÀ1 e 12 ; nÀ1 r 11 ; nÀ1 r 22 ; nÀ1 r 12 ; nÀ1 f 11 ; nÀ1 f 22 ; nÀ1 f 12 ) and three output nodes (i.e., n r 11 ; n r 22 ; n r 12 ) are needed while twenty four input nodes and six output nodes are required in the case of 3D brick linear elements accordingly.
The DLNN model is trained by a supervised learning scheme in the second step of Fig. 1. In particular, stress data from the output nodes are partially fed back into the input nodes of the DLNN recursively since the trained model is used in nonlinear FEA in that mode [39]. Before the training process, stresses and strains should be scaled by each of scale parameters since the stress and strain have different units. Those scaled values should be less than one for effective training.
where n r DLNN i is the scaled stress value from the DLNN model; n e DLNN j is the scaled strain value for the training of the DLNN model; and S r i and S e j are the constant scale values. Although DLNN can have a deeply hidden network maintaining sufficient learning capability, the formulation in the following assumes two hidden layers. Even if the DLNN can use more than five hidden layers, two hidden layers are sufficient to learn the complexity of the data in this study. Weight factors exist between each of the layers ( w rC ik ; w CB kl ; and w Be& lj ). The detailed notations of DLNN parameters are shown in Fig. 2. The Leaky-ReLU function [42] is utilized as an activation function for the DLNN, which is a key to learning with a deep network. Moreover, the adaptive moment estimation (Adam) optimizer [43] was utilized that is specially designed for training DLNN by minimizing the objective function. The objective function is expressed as Figure 2: The architecture of DLNN material constitutive model and notations where n r DLNN i; k ð Þ is the predicted stress component from DLNN; t i; k ð Þ is the corresponding stress component in the training database; i is the index for stress tensor; and k is the index for data pattern. The main advantage of the Adam method is that the step size is not affected by the gradient's rescaling. Besides, the step size can be updated by referring to the past gradient size. Even if the gradient increases, the step size is bound so that stable optimization is possible no matter what objective function is used. From the gradients g t of the objective function J W ð Þ, the Adam optimizer calculates two moment vectors; one is the average of g t called the 1st-moment vector m t and the other is the variance of g t called the 2nd-moment vector v t . Both moment vectors aids to easily escape out of local minima and to reliably reach a global minimum. However, Adam optimizer is suitable to use in online-training because it requires less effort to tune hyper-parameter (e.g., learning rate g) and save the total simulation time.
The gradient clipping method was applied to prevent the DLNN from exploding in this paper. If the gradient becomes larger during the training, weight parameters are updated with a large gradient and the DLNN model can diverge. It is the exploding gradient problem, which mostly happens in the recurrent neural networks. To deal with this problem, we made the gradient clipped with a limit value c. This method is called a gradient clipping [44].
TensorFlow [45,46] was utilized for training the DLNN material constitutive model. The DLNN-based material constitutive model was implemented in ABAQUS UMAT where weight factors of the trained DLNN are accessed through file opening-reading schemes and the Leaky-ReLU activation function is implemented. The trained DLNN will directly predict stresses at the next increment. Finally, the data-driven nonlinear FEA is conducted with the trained DLNN material constitutive model. We call the 3rd step of Fig. 1 as the forward DLNN-based FE simulation.

Algorithmic Tangent Operator for DLNN Material Constitutive Model
The algorithmic tangent operator (i.e., Jacobian matrix) is required for predicting nonlinear material constitutive behavior. The Jacobian matrix for the DLNN material constitutive model is derived as an explicit function of inputs, outputs and other DLNN parameters such as weight factors ( w rC ik ; w CB kl ; and w Be& lj ), scale factors (S r i and S e j ), the derivative form of the activation function and the activation function values from each of the hidden layers in the given DLNN ( n r DLNN i ; n C k ; and n B l ). The Jacobian matrix can be derived in terms of Leaky-ReLU activation function and other DLNN parameters as follows: where and : The convergence of this form in the case of multi-degrees of freedom FE models was confirmed as long as the DLNN correctly predicts stress values. The predicted stress and formulated material tangent stiffness matrix are updated at all Gauss points. The global stiffness matrix and external force vector for nonlinear FEA simulation can be updated from stress values and the Jacobian constitutive matrix.

Comparisons of ANN and DLNN Material Constitutive Models
To compare performances of the ANN and DLNN material constitutive models, a reference 2D plane stress FE model subjected to a tensile displacement loading was developed. The plasticity model with linear isotropic hardening was assumed for the reference model. The mesh and boundary conditions are depicted in Fig. 3 and material properties are summarized in Tab. 1.
From the reference simulation, stress and strain data were extracted from the whole gauss points within the FE model except for the right and left edges. The total number of training data was 15,040. Neural network structures were set to [12-20-20-3] for both ANN and DLNN models and they were trained until 50,000 epochs. The noticeable differences between ANN and DLNN models are summarized in Tab. 2.  The convergence of the cost function was compared for the ANN and DLNN models as shown in Fig. 4. The convergence of ANN was faster than DLNN until 10,000 epochs. However, the cost value convergence rate of DLNN outperformed ANN after 10,000 epoch, resulting in a much smaller cost value for DLNN than ANN.
Trained weight factors from the ANN and DLNN model at 50,000 epoch were used in the data-driven FEA. The global force-displacement responses were compared in Fig. 5 and the local responses were compared by enumerating von Mises stress contours as shown in Fig. 6. As expected, the DLNN-based FEA results were closer to the reference FEA results than ANN-based FEA.
For quantitative comparisons, the coefficient of determination (R 2 ) and root mean square (RMS) values were calculated by using the predicted stresses and strains from all FE integration points. They are defined as where n is the total number of data; t i is the true value in the training database; t is the mean of t i , and y i is the predicted value from DLNN. The averaged R 2 and averaged RMS values for each stress and strain component were calculated and summarized in Tabs. 3 and 4, respectively. Overall DLNN performance was better than ANN's. The mean R 2 values for DLNN are closed to one and much higher than the ones for ANN. Moreover, the mean RMS values of DLNN are relatively lower than the ones of ANN. According to the comparative analysis results, we conclude that DLNN is more suitable for material constitutive models than ANN in terms of training speed and accuracy of predictions.  The SELIFE can obtain complex 3D internal stress-strain states from a set of limited global reaction forces and displacements measured from material and structural tests. Applying inner surface deformation fields from the digital image correlation (DIC) sensor could be advantageous in identifying the 3D complex nonlinear stress-strain paths. Details of the methodology are described in the following.
The SELIFE simulation algorithm consists of repeatedly executing the forward DLNN-based FE simulations within an additional auto-adaptive training algorithm at each load (or time) increment. The SELIFE requires two additional iteration loops, which are DLNN pass, and DLNN auto-adaptive training cycle. Sweeping all load (or time) increments is called one DLNN Pass. Multiple DLNN passes may be necessary since the DLNN-based material model may not be trained with only one DLNN pass. In each of the load incremental steps, the SELIFE runs two independent forward DLNN-based FEA, which are the force-controlled analysis (FEM-A) and displacement-controlled analysis (FEM-B) performed within the auto-adaptive training cycles. The stresses and strains at the specified Gauss points from FEM-A and FEM-B, respectively, are appended to the training database. Stresses from FEM-A and strains from FEM-B are better for training the DLNN material constitutive model. Multiple auto-adaptive training cycles are performed until the predetermined number is reached or a convergence criterion in Eq. (11) is satisfied. The criterion is given as: where U FEM ÀA is the boundary displacement computed from FEM-A; U FEMÀB is the boundary displacement imposed to FEM-B; and Tol disp is the user-defined tolerance for auto-adaptive cycles. Thus, the criterion is checked as the criteria in each of the auto-adaptive training cycles. The DLNN model is gradually selforganized by the updated training dataset. Computational pseudo-codes, which include the two additional unique iteration loops for the SELIFE, are summarized in Tab. 7.
In Fig. 7, the iter1 and iter2 indicate gradual training of the DLNN model toward true stress-strain response through the auto-adaptive training cycles. After the SELIFE simulation, the DLNN model can be used in the forward DLNN-based nonlinear FE analysis. The SELIFE simulation is used to generate material "Big Data" in terms of stress-strain history data, which are subsequently used for establishing the failure criteria of new materials. All notations are summarized in Tab. 5. The pseudo-codes of the onlinetraining in the SELIFE simulation are illustrated in Tab. 6.      Mini-batch size for pre-training l 0 Initial mini-batch size l inc Increment of mini-batch size m i Current mini-batch size depending on j i b 1 and b 2 Exponential decays for the moments e Small constant value for preventing from dividing zero To develop failure criteria, the initial yield or failure stress data need to be identified based on the obtained multiaxial stress-strain paths from SELIFE analyses. We devised an identification algorithm that accurately pinpoint the failure initiation distinguishing different failure modes. The identification algorithm monitors the tangent stiffness along stress-strain curves. The algorithm for identifying failure initiation consists of two steps. The first step divides the whole multiaxial material response from the  Break // Escape out of auto-adaptive training cycle end end end SELIFE into nonlinear and linear response paths. The second step of the algorithm determines the stress response point corresponding to the initiation of failure within the whole range of the stress history.
The computational procedures of the first step are shown in Fig. 8a. At the beginning of the first step, stress and strain data are interpolated with a specific number of data points to reduce the effect of noise and error caused by trained DLNN on the calculation of tangent stiffness. It is named as the interpolated tangent stiffness expressed in Eq. (12).
where N indicates the number of interpolations. 50,000 data points were used to interpolate the stress-strain response in this paper. Then, the initial tangent stiffness was evaluated by using the interpolated tangent stiffness within the first two percent of the whole interpolation intervals as follows.
where n 0:02 describes the number of the first two percent interpolation interval.
Then the ratio of the interpolated tangent stiffness to the initial tangent stiffness (D) is evaluated in Eq. (14). It is followed by computation of the averaged ratio as in Eq. (15). The averaged ratio of the tangent stiffness d can classify the stress-strain path into nonlinear and linear regions. According to our experience, the averaged ratio of the tangent stiffness is close to one in the linear elastic range because the tangent stiffness in Eqs. (12) and (13). becomes close to each other within the linear elastic region. However, the averaged ratio of the tangent stiffness deviates from one with the nonlinear stress response regions because the interpolated tangent stiffness values in Eq. (12) become different from the initial tangent stiffness value in Eq. (13). Following this process, we still need an additional step to identify the initial failure point.
In the second step, initial failure points are identified from the nonlinear stress path. In the case of multiaxial stress state, the initial failure point corresponds to the point where the dominant stress component reaches the point earlier than other components Unlike the first step, d i À1 was used to capture the first deviating point from the initial tangent stiffness. According to our tests, the tolerance value in the range of 0.7~0.9 to be compared with d i À1 was appropriate to find the initial failure stresses. More detailed procedures for the second step of the proposed algorithm is depicted in Fig. 8b. Once the failure point is determined that corresponds to a specific stress component, then other stress components corresponding to the point are saved for subsequent development of failure criteria.

Symbolic Regression by Genetic Programming
Genetic programming (GP) is one of the machine learning variants, which is based on a bio-inspired technique. It can derive symbolic equations relating input data with output data out of a training database. It randomly produces tree architectures (e.g., a population consisting of individuals) and operates genetic mutation and cross-over process over the population. Throughout generations, the population converges to the best individual that gives the best performance. Then, the GP can provide a nonlinear equation by assembling weighted linear combinations from the final GP tree structure. An example of the GP tree model is shown in Fig. 9.
We utilized one of the genetic programming, called GPTIPS [48], implemented in MATLAB language, and generated an anisotropic yield surface out of the identified dataset. For the effective training of the GP, the data should be prepared. A more detailed process is described in the next section with an example simulation.

Verification of the SELIFE and Self-Learning Data-Driven Failure Criterion
For the demonstration of the proposed methodology, we developed a reference biaxial FE model with Hill's 48 anisotropic failure criterion, which is a reference solution but unknown in actual application to new materials. Therefore, Hill's 48 anisotropic failure criterion will be compared to the data-driven anisotropic failure criterion. Boundary force-displacement data from the reference simulated tests are deemed to be measurements from experimental tests. A large volume of multiaxial stress-strain data was self-learned and the identification algorithm was applied to discover initial yield stress points. Then, data-driven failure criteria were established by the GP. Finally, the data-driven failure criteria were verified by conducting simulations with different material orientations.

Simulated Reference Tests with Hill's Anisotropic Failure Criterion
The reference simulated test model is assumed to have anisotropic material with zero degree material orientation and subjected to four-types of displacement boundary conditions (DBC) summarized in Tab. 8. The specimen model is divide into five specific regions as shown in Fig. 10. Displacement of 0.01 mm magnitude was applied in Region A and B while Region C and D were fixed. The DBCs are summarized in Tab. 8. Stress-strain data were extracted from the Region E, where has 256 elements. The 2D plane stress condition was assumed and then the four-node bilinear plane stress (CPS4) element was utilized. Material properties for the simulation are summarized in Tab. 9. Hill's 48 anisotropic failure criterion [50] is defined as follows: where r ij are the stress components and F; G; H; L; M ; and N are constant parameters that express the current state of anisotropic behavior. The constant parameters are defined as functions of yield stress ratios R ij , the measured yields stress r y0 , and the material initial yield stress r y0 .   =12, 13, and 23) where s y0 ¼ r y0 ffiffi ffi 3 p For the plane stress condition, it is convenient to assume R ij as follows [51,52].
where r 0 and r 90 are the Lankford ratios. Moreover, the yield function is reduced to the plane stress condition as follows.
Applying Lankford values in Tab. 9, R ij in Eq. (18) are obtained. Eq. (21) provides the elliptic shape for anisotropic yield surface in the stress space (r 11 ; r 22 ; and, r 12 Þ. The intersections between the yield surface and each of axes define the measured yield stresses. They will be used for the verification of the new datadriven anisotropic yield surface. Applying the uniaxial tensile condition and the pure shear stress condition to Eq. (21), the measured yield stresses ( r 11 ; r 22 ; r 12 ) can be obtained. Both yield stress ratios and measured yield stresses are summarized in Tab. 10. The intersection points between the yield surface and axes were numerically calculated that are as shown in Fig. 11. Their values were identical to the measured yield stresses.
To confirm the anisotropic yield response, the reference simulation with boundary type 1 was investigated. The resultant force-displacement curves measured in Region A and B are different as shown in Fig. 11b.

SELIFE Simulation Results
SELIFE simulations were performed with all four boundary types separately. Fifteen nodes in each hidden layer of the DLNN were placed. The DLNN architecture in Eq. (22) was used for all SELIFE simulations. Online training was conducted with the same hyper-parameters for DLNN regardless of boundary types. All values for the SELIFE simulations are summarized in Tab. 11. r n ¼ r DLNN ½e n ; e nÀ1 ; r nÀ1 ; & n ð : 12 À 18 À 18 À 3 ½ Þ Effective plastic strain contour under tension-tension DBC (the boundary type 1) at the last analysis step is shown in Fig. 12. To show the local stress-strain path evolution from the SELIFE analysis, two elements were picked: 138th and 132nd elements representing elastic and plastic behavior, respectively.     Fig. 13 shows evidence of the gradual SELIFE learning of global responses. The global reaction forcedisplacement curve in x-direction, the global reaction force history in the y-direction, and reaction moment history from Region A are depicted in Figs. 13a, 13c, and 13e, respectively. Because of the material anisotropy and the effects of deformation in the other direction, the reaction moment was developed. The global reaction force history in the x-direction, the global reaction force-displacement in the y-direction, and reaction moment history from Region B are illustrated in Figs. 13b, 13d, and 13f, respectively. As evidenced, the global responses gradually approach their targeted curves as the SELIFE DLNN pass increases.
Figs. 14-16 compare local stresses trained by SELIFE with the results from the reference simulation. These results indicate that SELIFE is capable of learning true (i.e., reference) stress-strain relationships. The SELIFE learning progression of the elastoplastic behavior at the 132nd element is illustrated in Fig. 17. As the SELIFE DLNN pass increases, the stress-strain curves approach the reference curves. These results show promising evidence of the SELIFE performance that can learn nonlinear constitutive relationship from the experimentally measured global responses. These inverse learning results are followed by the identification algorithm to pinpoint failure initiation.

Identification of Failure Stress
By the proposed identification algorithm in Section 4.1, initial failure points were identified. As shown in Fig. 18, the identified failure points are indicated on the stress-strain curves with the failure points by Hill's anisotropic failure criterion. The initial yield points satisfying Hill's 48 anisotropic yield surface are exactly matching with yield points obtained by the proposed identification algorithm. The blue circle dots on each of the stress-strain curves are the initial yield points satisfying Hill's anisotropic failure criterion while the red asterisk dots are obtained from the identification algorithm. Although only three points are presented, the majority of the identified initial yield stresses were close to the initial yield stresses satisfying Hill's anisotropic failure criterion, according to our tests.
For example, in the case of Fig. 18a, the 11 stress component was the dominant one having maximum value and the initial yield stress was identified only with the 11 stress path. Next, the specific point within the path was saved and shared for the paired other stress component values (e.g., 22 and 12 components).
Following the feasibility study of the failure identification algorithm, two groups of the initial yield stresses were identified from all reference simulations and the self-learned stress-strain curves, respectively. Furthermore, all identified initial yield stresses, and their paired stress components were substituted into the Hill's anisotropic failure criterion Eq. (23).
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F r y 22 2 þ G r y 11 2 þ H r y 11 À r y where r y ij indicates the identified initial yield stresses from the self-learned stress-strain curves. Multiaxial stress states corresponding to the initial failure points from the self-learned stress-strain data are well-plotted on the anisotropic yield surface as shown in Fig. 19a. Moreover, the histogram of the difference defined in Eq. (23) in Fig. 19b indicates that most of the initial yield stresses are reasonably located on Hill's anisotropic yield surface. The average percent of the error distribution is 4.0654%. These results indicate fairly good agreement owing to the failure identification algorithm and autoadaptive algorithm.

Development of Data-Driven Failure Criterion by Genetic Programming
The failure dataset is prepared from all SELIFE simulations. We name the yield surface as a "selflearning data-driven yield surface". Furthermore, the data-driven yield surfaces are compared with Hill's anisotropic yield surface to verify the feasibility of the data-driven failure modeling methodology. Anisotropic constant parameters (F; G; H; and N) from the Hill's 48 anisotropic failure criterion should be determined by multiple experimental specimen tests.  However, the proposed self-learning data-driven methodology has a unique advantage of formulating anisotropic failure criterion without F; G; H; and N constants from a minimal number of experimental tests. We explain how to prepare training data to be used in the data-driven failure criteria modeling. Assuming a quadratic strength relationship, a general anisotropic failure criterion can be expressed as where r 2 yo indicates the initial yield stress of the given material and r 2 ij indicates the identified initial yield stresses. It is worth noting that F; G; H; and N constants are not included in Eq. (24). The inputs to the GP training are combinations of quadratic terms of the normal stresses (r 11 ; r 22 ), shear stress (r 11 À r 22 ) and initial yield stress (r yo ). r 2 12 was chosen as output of the general failure criterion. The input and output patterns for the GP training are summarized in Tab. 12.
Several parameters associated with the GP training are shown in Tab. 13 and the operational functions are shown in Tab. 14. The initial yield stress dataset was prepared by the proposed method and trained under the same conditions as in Tabs. 13 and 14.
Input data corresponding to the minimum r 2 12 value Minimum … ↓ ( Ascending order ) ↓ Input data corresponding to maximum r 2 12 value Maximum Finally, the GP symbolic regression was conducted with the database from all the SELIFE simulations. The data-driven anisotropic failure criterion is obtained as 0:2632 r 2 22 þ 0:2065 r 2 11 þ 0:7699 r 11 À r 22 Eq. (25) can produce a reasonable elliptic shape. The self-learning data-driven anisotropic yield surface Eq. (25) was compared with the reference Hill's 48 yield surface in Figs. 20a and 20b.
The present results are promising and innovative in the fact that the anisotropic yield surface of unknown material can be developed by using the minimal number of experimental tests.    25), measured yield stresses are back-calculated. Then, yield stress ratios (R 11 , R 22 , R 33 and R 12 ) were calculated by Eq. (18) and are listed in Tab. 15. Comparing with the reference yield stress ratios, yield stress ratios from the self-learning data-driven yield surface seem reasonably identical.
To verify the data-driven yield equations, additional simulations were conducted with different material orientations. In this study, we used a 22.5 degree of material orientation. The back-calculated anisotropic yield stress ratios in Tab. 15 were used for the verification simulations. Fig. 21 shows comparisons of global response under type 1 DBC for the two different anisotropic yield surfaces. Then, forcedisplacement data were extracted from Region A and B.
Moreover, stress-strain data were obtained from the specific elements and those are shown in Fig. 22. Comparing stress-strain curves, initial yield points from the self-learning data-driven yield surface are very close to those from the reference anisotropic yield surface. These results indicate that the proposed methodology can develop the data-driven failure criteria of unknown materials from experimental measurements.  6 Conclusions This paper presented a novel self-learning data-driven modeling methodology that can develop failure criteria of an unknown anisotropic ductile material. The methodology consists of the DLNN-based material constitutive model, the self-learning inverse finite element (SELIFE) simulation, failure identification algorithm through data-processing, and derivation of failure criteria through symbolic regression genetic programming (GP). The SEFIFE can gradually train the DLNN material constitutive model to follow true nonlinear stress-strain relationship being guided by the experimental measurements. In particular, stress update and the algorithmic tangent operator were derived in terms of DLNN parameters to be used in nonlinear finite element codes.
1. The DLNN material constitutive model was superior to the artificial neural network (ANN) based model in terms of computational efficiency and prediction accuracy. The DLNN material constitutive model could extend learning capabilities with massively hidden networks owing to the Leaky-ReLU activation function. 2. The SELIFE simulation renders the DLNN material constitutive model to self-learn multiaxial nonlinear stress-strain relationship with inputs of experimentally measured boundary reaction forces and displacements. Local stress distributions and the global responses of the self-learning model were fairly close to those of the reference model. Computational procedures for the SELIFE simulation were also presented in pseudocodes. 3. Identification techniques of the failure points based on tangent stiffnesses could reasonably pinpoint the failure points. 4. Finally, the GP application technique were presented to develop anisotropic failure criteria of unknown anisotropic ductile materials. 5. The proposed idea was verified by direct comparisons with the reference simulation results with Hill's 48 anisotropic failure criteria and further tests with material having different material orientations.
The proposed methodology has significances in capabilities of developing failure criteria of uncharacterized materials such as 3D printed materials having process-dependent properties. Beyond the Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.