This paper first 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 finite element (SELIFE) simulation, 3) Algorithmic identification 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 finite 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 identification, 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 verified by comparing with a reference failure criteria and simulating with different materials orientations, respectively.
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–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–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–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–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–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
n−1ςi=n−1σin−1εi+n−1σiΔnεi
where, n−1σi and n−1εi are the stress and strain vector at (n − 1)-th increment, respectively; Δnε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ε11,nε22,nε12,n−1ε11,n−1ε22,n−1ε12,n−1σ11,n−1σ22,n−1σ12,n−1ζ11,n−1ζ22,n−1ζ12) and three output nodes (i.e., nσ11,nσ22,nσ12) are needed while twenty four input nodes and six output nodes are required in the case of 3D brick linear elements accordingly.
A flow of the data-driven FEA with DLNN material constitutive model
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.
nσiDLNN=nσiSiσ,where|nσiDLNN|<1
nεjDLNN=nεjSjε,where|nεjDLNN|<1
where nσiDLNN is the scaled stress value from the DLNN model; nεjDLNN is the scaled strain value for the training of the DLNN model; and Siσ and Sjε 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 (wikσC,wklCB, and wljBες). 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
where nσi,(k)DLNN is the predicted stress component from DLNN; ti,(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 gt of the objective function J(W), the Adam optimizer calculates two moment vectors; one is the average of gt called the 1st-moment vector mt and the other is the variance of gt called the 2nd-moment vector vt. 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 η) and save the total simulation time.
The architecture of DLNN material constitutive model and notations
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 (wikσC,wklCB, and wljBες), scale factors (Siσ and Sjε), the derivative form of the activation function and the activation function values from each of the hidden layers in the given DLNN (nσiDLNN,nCk, and nBl). The Jacobian matrix can be derived in terms of Leaky-ReLU activation function and other DLNN parameters as follows:
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.
Tensile simulation modelMaterial properties for the reference FE model
Properties
Values
Young’s modulus, E [MPa]
210,000
Poisson’s ratio, ν
0.3
Initial yield stress, σy0 [MPa]
182
Plastic modulus, K [MPa]
5500
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.
Different features between ANN and DLNN
ANN
DLNN
Implemented language
Fortran
TensorFlow
Optimizer
Resilient backpropagation
Adam
Activation function
Hyperbolic tangent
Leaky-ReLU
The computational time for training
59 min and 47 s
6 min
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.
Comparison of variations of cost functions for ANN and DLNN during training
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.
Comparison of global responses (force-displacement) at the right edgeComparison of von Mises stress contours. (a) Reference FE model with material properties of <xref ref-type="table" rid="table-1">Tab. 1</xref>. (b) ANN-based FEA results and (c) DLNN-based FEA results
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; ti is the true value in the training database; t¯ is the mean of ti, and yi 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.
Comparison of mean R<sup>2</sup> values for each component
Mean R^{2}
11 component
22 component
12 component
σijANN
0.92756
0.79557
0.87239
σijDLNN
0.99857
0.99744
0.99971
εijANN
0.80812
−0.32978
0.67306
εijDLNN
0.99932
0.98612
0.99963
Comparison of mean RMS values for each component
Mean RMS
11 component
22 component
12 component
σijANN
7.8012
4.0625
6.0421
σijDLNN
0.31521
0.35144
0.21614
εijANN
1.0965e-04
1.7392e-05
1.8988e-04
εijDLNN
4.4765e-06
2.0808e-06
8.5343e-06
Self-Learning Inverse Finite Element (SELIFE) Analysis
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:
edisp=UFEM−A−UFEM−B2UFEM−B2<Toldisp,
where UFEM−A is the boundary displacement computed from FEM-A; UFEM−B is the boundary displacement imposed to FEM-B; and Toldisp 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 self-organized 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 online-training in the SELIFE simulation are illustrated in Tab. 6.
Auto-adaptive learning of true material responseNotations and definitions
Notations
Descriptions
τi∈{τInitial,τ2,…,τFinal}
Current increment in SELIFE simulation
Nες,NB,NC,Nσ
Number of nodes in each layer
Epre
Pre-training Epoch
EAuto
Auto-adaptive training Epoch
α
The negative slope for Leaky-ReLU
η
Learning rate
c
Gradient clipping value
μpre
Mini-batch size for pre-training
μ0
Initial mini-batch size
μinc
Increment of mini-batch size
mi
Current mini-batch size depending on ξi
β1 and β2
Exponential decays for the moments
e
Small constant value for preventing from dividing zero
Pseudo-codes for DLNN in the SELIFE simulation
Online-training algorithm
Require:τi,Nες,NB,NC,Nσ,Epre,EAuto,α,η,c,μpre,μ0, and μinc
Output:W={wikσC,wklCB,wljBες} trained weight
1. Set hyper-parameters: α,η,c,μpre,μ0,μinc,Epre, and EAuto
2. Establish DLNN structure based on Nες,NB,NC, and Nσ
3. Setting type of training database, weights, mini-batch size, and Epoch
4. Scan train dataset {X(k),T(k)}k=1k=NDB to get the size of train dataset NDB
5. Obtain scaled train dataset {x(k),t(k)}k=1k=NDB using Eqs. (2) and (3)
6. Calculate the mini-batch learning iteration size : NMB=floor(NDBmi)
7. Data training
fori_epoch=1:Epoch
Data-preprocessing: Randomly shuffle and arrange scaled train dataset
forj=1:NMB
Sample scaled training data : {x(k),t(k)}k=1+mi⋅(j−1)k=mi⋅j
Obtain cost function (MSE) : J(W)=1mi[∑k=1+mi⋅(j−1)k=mi⋅j(nσi,(k)DLNN−ti,(k))2]
Calculate gradient: gW←∇WJ(W)
Gradient clipping [44]: gWclip←clip(gW,[−c,c])
Optimization [43]: W←W+η⋅Adam(gwclip,β1,β2,e) with default values
end
end
8. Save trained weight, W←W+η⋅Adam(gwclip,β1,β2,e)
Pseudo codes for SELIFE simulation
SELIFE algorithm
Require: ABAQUS FE input file with material information and boundary condition
Require: ABAQUS UMAT which implements DLNN feedforwarding and Jacobian matrix
Require: Text files which contain force and displacement boundary data
Require: Text file which controls SELIFE simulation (e.g. Epre,EAuto,α,η,c,μpre,μ0, the number of passes NpassDLNN, the number of cycles NcycleDLNN, and the size of DLNN)
Output: Stress-strain data and weight parameters at every point of time increment τi
1. Obtain the whole control variables of the SELIFE
2. Scan and save the ABAQUS FE input file and ABAQUS UMAT in string data type
3. Scan and save the force and displacement boundary data
4. Determine the simulation increment τi∈{τInitial,τ2,…,τFinal}
5. Define the number of gauss points, tensor, and degree of freedom which correspond to the ABAQUS FE model
6. Pre-training
7. SELIFE simulation start
foripass=1:NpassDLNN // Loop for DLNN pass
forτi=τInitial:τFinal // Loop for incremental FE simulation
foricycle=1:NcycleDLNN // Cycle for auto-adaptive training
Rewrite ABAQUS UMAT with trained weight parameters
FEM-A: DLNN-based FEA with force boundary data at τi
FEM-B: DLNN-based FEA with displacement boundary data at τi
Data extraction: stress data from FEM-A and strain data from FEM-B
Append the extracted stress-strain data to the training database
Auto-adaptive training (Step 6. Pre-training and Fig. 7)
If Check convergence criterion Eq. (11)
Break // Escape out of auto-adaptive training cycle
end
end
end
Development of Failure Criteria from Self-Learning Data for Unknown Ductile Anisotropic MaterialsIdentification Algorithm of Failure Initiation
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 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).
Δi=σi+1−σiεi+1−εii∈(1,2,…,N−1)
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.
Δ_=∑j=1nΔjn0.02
where n0.02 describes the number of the first two percent interpolation interval.
δi=ΔiΔ_i∈(1,2…,N−1)
⟨δ⟩=∑i=1N−1δiN−1
(a) Flowchart for the first step for stress history classification. (b) Flowchart for the second step of the identification algorithm
Then the ratio of the interpolated tangent stiffness to the initial tangent stiffness (Δ_) 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 δ 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, |δ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 |δ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.
An example of the symbolic regression model by multigene
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.
The biaxial specimen with boundary conditions and detail dimension in [mm] scale for the reference simulationFour types of displacement boundary conditions
Boundary types
Region A
Region B
Type 1
Tension
Tension
Type 2
Compression
Compression
Type 3
Compression
Tension
Type 4
Tension
Compression
Material properties of the anisotropic metal (DDQ mild steel [<xref ref-type="bibr" rid="ref-49">49</xref>])
Properties
Values
Young’s modulus, E [GPa]
206
Poisson’s ratio, ν
0.3
Initial yield stress, σy0 [MPa]
152.0
Lankford ratios : r0,r45,r90(state of anisotropy of the material)
2.64, 1.57, 2.17
Hill’s 48 anisotropic failure criterion [50] is defined as follows:
where σ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 Rij, the measured yields stress σ¯y0, and the material initial yield stress σy0.
Rij=σ¯ijτy0 (for ij =12, 13, and 23) where τy0=σy03
For the plane stress condition, it is convenient to assume Rij as follows [51,52].
R11=R13=R23=1
R22=r90(r0+1)r0(r90+1), R33=r90(r0+1)r0+r90, R12=3r90(r0+1)(2r0+1)(r0+r90) where r0 and r90 are the Lankford ratios. Moreover, the yield function is reduced to the plane stress condition as follows.
f(σ)=Fσ222+Gσ112+H(σ11−σ22)2+2Nσ122
When plastic deformation occurs, the yield function satisfies Eq. (21).
Fσ222+Gσ112+H(σ11−σ22)2+2Nσ122=σy02
Applying Lankford values in Tab. 9, Rij in Eq. (18) are obtained. Eq. (21) provides the elliptic shape for anisotropic yield surface in the stress space (σ11,σ22, and, σ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 data-driven anisotropic yield surface. Applying the uniaxial tensile condition and the pure shear stress condition to Eq. (21), the measured yield stresses (σ¯11,σ¯22,σ¯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.
(a) Hill’s 48 anisotropic yield surface with intersections and (b) force-displacement curves from Hill’s anisotropic failure criterionCalculated yield stress ratios and measured yield stress
Components
Yield stress ratios
Measured yield stresses [MPa]
11 component
R11=1.00
σ¯11=±152.00
22 component
R22=0.97
σ¯22=±147.67
12 component
R12=1.09
σ¯12=±95.73
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.
σn=σDLNN([εn,εn−1,σn−1,ςn]:[12−18−18−3])
Empirically set hyper-parameters for DLNN and parameters of the SELIFE simulation
Notations
# ofNN pass
Toldisp
Epre
EAuto
α
η
c
μpre
μ0
μinc
Values
10
2E-02
1000
5
0.1
5e-04
1e-04
49
200
300
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.
Effective plastic strain contour from the reference simulation under tension-tension displacement boundary condition (boundary type 1)
Fig. 13 shows evidence of the gradual SELIFE learning of global responses. The global reaction force-displacement 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.
Learning progression of global responses by SELIFE DLNN training; SELIFE simulation under tension-tension displacement boundary condition (boundary type 1). (a) Region A: F-D (x-direction). (b) Region B: F-time (x-direction). (c) Region A: F-time (y-direction). (d) Region B: F-D (y-direction). (e) Region A: F-time (moment). (f) Region B: F-time (moment)
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.
Stress (S11) contour comparison. (a) Result of the reference simulation and (b) result of SELIFE simulation at the last pass (the 10th pass)Stress (S22) contour comparison. (a) Result of the reference simulation and (b) result of SELIFE simulation at the last pass (the 10th pass)Stress (S12) contour comparison. (a) Result of the reference simulation and (b) result of SELIFE simulation at the last pass (the 10th pass)
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.
Stress-strain curves at the 3rd gauss point in the 132nd element by SELIFE DLNN passesIdentification 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.
Comparison of the initial yield stress positions. (a) At the 2nd gauss point in the 131st element. (b) At the 2nd gauss point in the 184th element, and (c) at the 2nd gauss point in the 240th element
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).
Fσ22y2+Gσ11y2+H(σ11y−σ22y)2+2Nσ12y2−σy0=error
where σijy 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 auto-adaptive algorithm.
The identified initial yield stresses from all SELIFE stress-strain data. (a) Identified initial yield stresses on the Hill’s yield surface. (b) Distribution of errors of the identified initial yield stresses in <xref ref-type="disp-formula" rid="eqn-23">Eq. (23)</xref>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 “self-learning 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
σ122=f(σyo2,σ222,σ112,(σ11−σ22)2)
where σyo2 indicates the initial yield stress of the given material and σij2 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 (σ11,σ22), shear stress (σ11−σ22) and initial yield stress (σyo). σ122 was chosen as output of the general failure criterion. The input and output patterns for the GP training are summarized in Tab. 12.
Input and output data pattern for the genetic program
Input data
Output data
X1=σyo2
X2=σ222
X3=σ112
X4=(σ11−σ22)2
Y1=σ122
Input data corresponding to the minimum σ122 value
Minimum
…
↓( Ascending order )↓
Input data corresponding to maximum σ122 value
Maximum
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.
Parameters for genetic program
Parameter description
Value
The number of population size
100
The number of generation
100
Maximum number of genes per individual
2
Maximum depth for gene
2
Activated functions for the genetic programming
Activated function name
Symbols
Times
×
Minus
−
Plus
+
Finally, the GP symbolic regression was conducted with the database from all the SELIFE simulations. The data-driven anisotropic failure criterion is obtained as
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.
Comparison of the self-learning data-driven anisotropic yield surfaces with Hill’s yield surface. (a) 3D stress space. (b) 2D stress plane
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.
Verification of the Self-Learning Data-Driven Failure criterion
Based on the anisotropic parameters (F,G,H, and N) of the data-driven yield surface i.e., Eq. (25), measured yield stresses are back-calculated. Then, yield stress ratios (R11, R22, R33 and R12) 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.
Comparisons of yield surface ratios for the reference yield surface and the self-learning data-driven yield surface (<inline-formula id="ieqn-183">
<alternatives><inline-graphic xlink:href="ieqn-183.png"/><tex-math id="tex-ieqn-183"><![CDATA[${\bi{R}_{11}}$]]></tex-math><mml:math id="mml-ieqn-183"><mml:mrow><mml:msub><mml:mi mathvariant="bold-italic">R</mml:mi><mml:mrow><mml:mn>11</mml:mn></mml:mrow></mml:msub></mml:mrow></mml:math>
</alternatives></inline-formula>, <inline-formula id="ieqn-184">
<alternatives><inline-graphic xlink:href="ieqn-184.png"/><tex-math id="tex-ieqn-184"><![CDATA[${\bi{R}_{22}}$]]></tex-math><mml:math id="mml-ieqn-184"><mml:mrow><mml:msub><mml:mi mathvariant="bold-italic">R</mml:mi><mml:mrow><mml:mn>22</mml:mn></mml:mrow></mml:msub></mml:mrow></mml:math>
</alternatives></inline-formula>,<inline-formula id="ieqn-185">
<alternatives><inline-graphic xlink:href="ieqn-185.png"/><tex-math id="tex-ieqn-185"><![CDATA[$\; {\bi{R}_{33}}$]]></tex-math><mml:math id="mml-ieqn-185"><mml:mspace width="thickmathspace"></mml:mspace><mml:mrow><mml:msub><mml:mi mathvariant="bold-italic">R</mml:mi><mml:mrow><mml:mn>33</mml:mn></mml:mrow></mml:msub></mml:mrow></mml:math>
</alternatives></inline-formula>, and <inline-formula id="ieqn-186">
<alternatives><inline-graphic xlink:href="ieqn-186.png"/><tex-math id="tex-ieqn-186"><![CDATA[${\bi{R}_{12}}$]]></tex-math><mml:math id="mml-ieqn-186"><mml:mrow><mml:msub><mml:mi mathvariant="bold-italic">R</mml:mi><mml:mrow><mml:mn>12</mml:mn></mml:mrow></mml:msub></mml:mrow></mml:math>
</alternatives></inline-formula>)
R11
R22
R33
R12
Reference material from Tab. 9 and Eq. (21)
1.0000
0.9715
1.2815
1.0909
Self-learning data-driven yield surface Eq. (25)
1.0120
0.9838
1.4591
0.9820
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, force-displacement data were extracted from Region A and B.
Force-displacement comparison from (a) the Region A and (b) the Region 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.
Comparison of stress-strain curves at the 3rd gauss point in the 186th element. (a) 11 component. (b) 22 componentConclusions
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.
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.
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.
Identification techniques of the failure points based on tangent stiffnesses could reasonably pinpoint the failure points.
Finally, the GP application technique were presented to develop anisotropic failure criteria of unknown anisotropic ductile materials.
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 generalization of the self-learned DLNN model, it can also construct material-related “Big Data.” The proposed methodology can also be extended to discover damage evolution laws in the future.
Funding Statement: This work was supported by the National Research Foundation of Korea (NRF) grant of the Korea government (MSIP) (2020R1A2B5B01001899) (Grantee: GJY, http://www.nrf.re.kr) and Institute of Engineering Research at Seoul National University (Grantee: GJY, http://www.snu.ac.kr). The authors are grateful for their supports.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
ReferencesJ.Ghaboussi, D. A.Pecknold, M.Zhang and R. M.Haj-Ali, “
Autoprogressive training of neural network constitutive models,”
C.Settgast, G.Hütter, M.Kuna and M.Abendroth, “
A hybrid approach to simulate the homogenized irreversible elastic-plastic deformations and damage of foams by neural networks,”
A.Oishi and G.Yagawa, “
Computational mechanics enhanced by deep learning,”
M. A.Bessa, R.Bostanabad, Z.Liu, A.Hu, D. W.Apleyet al., “
A framework for data-driven analysis of materials under uncertainty: Countering the curse of dimensionality,”
B. A.Le, J.Yvonnet and Q. C.He, “
Computational homogenization of nonlinear elastic materials using neural networks,”
G. J.Yun, J.Ghaboussi and A. S.Elnashai, “
A new neural network-based model for hysteretic behavior of materials,”
T.Furukawa and G.Yagawa, “
Implicit constitutive modeling for viscoplasticity using neural networks,”
R.Haj-Ali and H. K.Kim, “
Nonlinear constitutive models for FRP composites using artificial neural networks,”
Z.Zhang and K.Friedrich, “
Artificial neural networks applied to polymer composites: A review,”
H. ElKadi, “
Modeling the mechanical behavior of fiber-reinforced polymeric composite materials using artificial neural networks—A review,”
Y.Shen, K.Chandrashekhara, W. F.Breig and L. R.Oliver, “
Neural network based constitutive model for rubber material,”
B.Wang, J. H.Ma and Y. P.Wu, “
Application of artificial neural network in prediction of abrasion of rubber composites,”
S.Jung and J.Ghaboussi, “
Characterizing rate-dependent material behaviors in self-learning simulation,”
T.Furukawa and G.Yagawa, “
Implicit constitutive modelling for viscoplasticity using neural networks,”
R. M.Diaconescu, M.Barbuta and M.Harja, “
Prediction of properties of polymer concrete composite with tire rubber using neural networks,”
A. M.Hassan, A.Alrashdan, M. T.Hayajneh and A. T.Mayyas, “
Prediction of density, porosity and hardness in aluminum–copper-based composite materials using artificial neural network,”
V. M.Nguyen-Thanh, X.Zhuang and T.Rabczuk, “
A deep energy method for finite deformation hyperelasticity,”
K.Wang and W. C.Sun, “
Meta-modeling game for deriving theory-consistent, microstructure-based traction-separation laws via deep reinforcement learning,”
Z.Liu and C. T.Wu, “
Exploring the 3D architectures of deep material network in data-driven multiscale mechanics,”
M.Raissi and G. E.Karniadakis, “
Hidden physics models: Machine learning of nonlinear partial differential equations,”
J.Sirignano and K.Spiliopoulos, “
DGM: A deep learning algorithm for solving partial differential equations,”
J.Berg and K.Nystrom, “
A unified deep artificial neural network approach to partial differential equations in complex geometries,”
H.Guo, X.Zhuang and T.Rabczuk, “
A deep collocation method for the bending analysis of kirchhoff plate,”
C.Anitescu, E.Atroshchenko, N.Alajlan and T.Rabczuk, “
Artificial neural network methods for the solution of second order boundary value problems,”
E.Weinan and B.Yu, “
The deep ritz method: A deep learning-based numerical algorithm for solving variational problems,”
I. E.Lagaris, A.Likas and D. I.Fotiadis, “
Artificial neural networks for solving ordinary and partial differential equations,”
T.Kirchdoerfer and M.Ortiz, “
Data-driven computational mechanics,”
T.Kirchdoerfer and M.Ortiz, “
Data-driven computing in dynamics,”
R.Eggersmann, T.Kirchdoerfer, S.Reese, L.Stainier and M.Ortiz, “
Model-free data-driven inelasticity,”
T.Kirchdoerfer and M.Ortiz, “
Data driven computing with noisy material data sets,”
L.Stainier, A.Leygue and M.Ortiz, “
Model-free data-driven methods in mechanics: Material data identification and solvers,”
F.Abbassi, T.Belhadj, S.Mistou and A.Zghal, “
Parameter identification of a mechanical ductile damage using artificial neural networks in sheet metal forming,”
A.Leygue, M.Coret, J.Réthoré, L.Stainier and E.Verron, “
Data-based derivation of material response,”
P.Ladevèze, D.Néron and P. W.Gerbaud, “
Data-driven computation for history-dependent materials,”
W. K.Liu, G.Karniadakis, S.Tang and J.Yvonnet, “
A computational mechanics special issue on: Data-driven modeling and simulation—Theory, methods, and applications,”
M.Dalemat, M.Coret, A.Leygue and E.Verron, “
Measuring stress field without constitutive equation,”
L. T. K.Nguyen and M. A.Keip, “
A data-driven approach to nonlinear elasticity,”
J.Ghaboussi, J. H.Garrett and X.Wu, “
Knowledge-based modeling of material behavior with neural networks,”
G. J.Yun, J.Ghaboussi and A. S.Elnashai, “
A new neural network-based model for hysteretic behavior of materials,”
Y. M. A.Hashash, S.Jung and J.Ghaboussi, “
Numerical implementation of a neural network based material model in finite element analysis,”
G. J.Yun, A.Saleeb, S.Shang, W.Binienda and C.Menzemer, “
Improved selfsim for inverse extraction of non-uniform, nonlinear and inelastic constitutive behavior under cyclic loadings,”
A. L.Maas, A. Y.Hannun and A. Y.Ng, “
Rectifier nonlinearities improve neural network acoustic models,” in Proc. of the 30th Int. Conf. on Machine Learning, Atlanta, Georgia, USA, pp.
3,
2013.D. P.Kingma and J.Ba, “
Adam: a method for stochastic optimization,”
2014. [Online]. Available: https://arxiv.org/abs/1412.6980.R.Pascanu, T.Mikolov and Y.Bengio, “
On the difficulty of training recurrent neural networks,” in Proc. of the 30th Int. Conf. on Machine Learning, Atlanta, Georgia, USA, PMLR, vol. 28, no. 3, pp.
1310–
1318,
2013.M.Abadi, A.Agarwal, P.Barham, E.Brevdo, Z.Chenet al., “
Tensorflow: Large-scale machine learning on heterogeneous distributed systems,”
2016. [Online]. Available: https://arxiv.org/abs/1603.04467.A.Géron,
K.He, X.Zhang, S.Ren and J.Sun, “
Delving deep into rectifiers: Surpassing human-level performance on lmageNet classification,” in 2015 IEEE Int. Conf. on Computer Vision, Santiago, pp.
1026–
1034,
2015.D. P.Searson, D. E.Leahy and M. J.Willis, “
GPTIPS: An open source genetic programming toolbox for multigene symbolic regression,” in Proc. of the Int. Multi Conf. of Engineers and Computer Scientists, 2010 Vol I, IMECS2010, March17–19, 2010, Hong Kong, pp.
77–
80,
2010.Y.Choi, C. S.Han, J. K.Lee and R. H.Wagoner, “
Modeling multi-axial deformation of planar anisotropic elasto-plastic materials, part I: Theory,”
R.Hill, “
A theory of the yielding and plastic flow of anisotropic metals,”
S.Bagherzadeh, M. J.Mirnia and B. M.Dariani, “
Numerical and experimental investigations of hydro-mechanical deep drawing process of laminated aluminum/steel sheets,”
D. V.Hai, S.Itoh, T.Sakai, S.Kamado and Y.Kojima, “
Experimentally and numerical study on deep drawing process for magnesium alloy sheet at elevated temperatures,”