Deep material network with cohesive layers: Multi-stage training and interfacial failure analysis
DDeep material network with cohesive layers: Multi-stage training andinterfacial failure analysis
Zeliang Liu a, ∗ a Livermore Software Technology, an ANSYS company, Livermore, CA 94551, USA
Abstract
A fundamental issue in multiscale materials modeling and design is the consideration of traction-separationbehavior at the interface. By enriching the deep material network (DMN) with cohesive layers, the paperpresents a novel data-driven material model which enables accurate and efficient prediction of multiscaleresponses for heterogeneous materials with interfacial effect. In the newly invoked cohesive building block,the fitting parameters have physical meanings related to the length scale and orientation of the cohesivelayer. It is shown that the enriched material network can be effectively optimized via a multi-stage trainingstrategy, with training data generated only from linear elastic direct numerical simulation (DNS). Theextrapolation capability of the method to unknown material and loading spaces is demonstrated through thedebonding analysis of a unidirectional fiber-reinforced composite, where the interface behavior is governed byan irreversible softening mixed-mode cohesive law. Its predictive accuracy is validated against the nonlinearpath-dependent DNS results, and the reduction in computational time is particularly significant.
Keywords:
Machine learning; Model reduction; Path-dependency; Composites; Debonding analysis
1. Introduction
Materials have deformable interfaces across fine-scale phases, and an emerging need in multiscale ma-terials modeling and design is to capture the effects of interfaces at various length and time scales. Somerepresentative examples are fiber-matrix debonding in carbon fiber reinforced polymer (CFRP) composites[1], Mullins effect due to stress softening at the particle/matrix interface [2], void nucleation in ductilefracture of metals [3, 4], and grain-boundary effects in polycrystalline materials [5, 6]. One can see thatthe interface, in many cases, is the critical region involving damage initiation and evolution, thus, it playsan important role in the failure analysis of modern material systems. On the other hand, as the interfaceis a geometric object with one dimension less than the bulk material, it also results in various interestingphenomena related to the size effect of microstructures [7, 8].A plethora of multiscale material models have been developed to consider the interfacial effect. However,empirical models are problem-dependent and may lose the essential physics of the interface. Analyticalmicromechanics methods based on the Eshelby’s solution [9] and mean-field theories [10, 11] can be enhancedto incorporate the interfacial effect, but they are usually limited to regular geometries and weakly nonlinearbehaviors [12, 13, 8]. To resolve these issues, computational homogenization is emerging as an accurateway of deriving multiscale material responses from the representative volume element (RVE) [14], where themicrostructures are explicitly constructed and discretized by numerical models.In recent years, direct numerical simulation (DNS) of RVE has been widely used to understand therole of interface in heterogeneous materials. To describe the traction-separation relationships across theinterfaces, various cohesive laws have been developed [15] and further built into DNS tools with cohesive-element formulations, such as the finite element method (FEM) [16, 17], the extended finite element method ∗ Corresponding author.
Email address: [email protected] (Zeliang Liu) a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b XFEM) [18, 19], the generalized finite element method (GFEM) [5], and isogeometric analysis (IGA)[20, 21]. The pioneer work done by Needleman [7] simulated a periodic array of rigid spherical inclusionsto investigate the cause of void nucleation due to the inclusion debonding. Melro et al. [1] developed RVEmodels for polymer composite reinforced by unidirectional fibers and evaluated the influence of fiber-matrixinterface on the strength properties of the composite. Zhao et al. [22] utilized the XFEM to simulatethe multi-void nucleation process due to interface debonding under various loading conditions and particlesize/shape distributions. However, limited by the large computational cost of DNS, these full-field RVEmodels were mainly used to validate and calibrate macroscopic empirical material models.In the perspective of integrated multiscale frameworks like concurrent simulation and materials design[23, 24], the efficiency of RVE analysis arises as a fundamental issue. Many reduced-order methods havebeen proposed for fast RVE homogenization with interfacial failure, such as the Voronoi cell method [25,26], eigendeformation-based method[27, 28, 29], self-consistent clustering analysis [30, 23, 31, 32, 33], andproper orthogonal decomposition (POD) [34]. Specifically, the SCA and POD-based methods are data-driven since their reduced bases are extracted from mechanical data in the offline training stage. Feed-forward neural network and recurrent neural network have also been used to fit the stress-strain responsesof elastic and history-dependent materials [35, 36, 37, 38, 33]. Regarding interface modeling, Wang etal.[39] recently proposed a deep reinforcement learning framework based on directed graph for derivingthe traction-separation law. As the training data is usually limited due to the expense of physical ornumerical experiments, a key challenge faced by most data-driven methods is the credibility of extrapolationto unknown material and loading spaces, especially when nonlinear softening interfaces are involved.An alternative way of mining the mechanical data is to learn the hidden topology representation ofthe RVE. Following this idea, Liu et al. [40, 41] proposed the deep material network (DMN) for data-driven material modeling, which describes the RVE by a network structure based on physics-based buildingblocks. Essentially, all fitting parameters of DMN are interpretable with physical meanings related to theRVE morphology. Trained only by linear elastic data, DMN can be extrapolated to capture highly nonlinearmaterial behaviors for various applications, including hyperelastic rubber composite under large deformation,polycrystalline materials with rate-dependent crystal plasticity, and various types of CFRPs. Recently, atransfer-learning framework of DMN is further introduced to generate a unified set of material databasescovering the full-range structure-property relationship [24].This paper extends the deep material network framework [41] to 3-dimensional heterogeneous materialswith deformable interfaces. The main contributions of the work are: • New architecture of material network enriched by adaptable cohesive networks; • Physics-based cohesive building block where a reciprocal length parameter is introduced to describethe size effect of interface; • Multi-stage training strategy to reduce the number of fitting parameters in the machine-learningprocess at each stage; • Extrapolation of the enriched DMN from linear elasticity to interfacial failure analysis with irreversiblesoftening cohesive laws.The remainder of this paper is organized as follows. In Section 2, the new design of DMN architecturewith cohesive layers is discussed in detail. Analytical solutions of the cohesive building block are derived toprovide the basis of data propagation in the network. Section 3 presents the multi-stage training strategybased on linear elastic offline DNS data. Section 4 briefly introduces the mix-mode cohesive law and necessaryviscous regularization for implicit analysis in the online stage. Section 5 shows the application of DMN tothe interfacial debonding analysis of a unidirectional fiber-reinforced composite, and another example ona particle-reinforced composite is also provided in Appendix A. Concluding remarks and future work aregiven in Section 6. 2 wo-layer building block
Material network
CL 1 CL 𝑝 CL 𝑁 $𝑣 &’ = 𝑎 ̃𝑧 &’ /𝐿$𝑣 &. = 𝑎 ̃𝑧 &. /𝐿$𝑣 &/ = 𝑎 ̃𝑧 &/ /𝐿 Base material …… Cohesive network 𝒒 … Cohesive building block
Cohesive network Parameters: 𝑧, 𝛼, 𝛽, 𝛾
Parameters: ̃𝑧, $𝛼, 7𝛽, $𝛾 …… Enriched node 𝒒 Enriched node 1 … … Depth 𝑁 Outputs
Figure 1: Proposed architecture of deep material network enriched by cohesive networks. The depth of material network is N . The number of cohesive layers in each cohesive network is N c . The reciprocal length parameter of the p -th cohesive layerin cohesive network q is denoted by ˜ v pq , and L is the characteristic length of the microstructure predefined for the better oftraining (see Section 2.2).
2. New DMN architecture with enrichment
As a starting point, consider a material network based on a binary tree with a two-layer structure asthe building block [40, 41], as shown on the right of Figure 1. The microscale constituents are input tothe bottom layer, and the mechanical information such as the stiffness tensor and the residual stress is feedforward through the input-output functions defined by the physics-based building block. The macroscalematerial properties are extracted from the top node of the network. Let N be the depth of material network,so that there are in total 2 N − nodes at the bottom layer ( i = N ).The basic hypothesis of DMN is that the original 3D RVE can be represented by patterns of the two-layerstructures with unknown phase fractions and block orientations [41]. The fitting parameters of a materialnetwork are the activation z of each bottom-layer node, and three rotation angles ( α, β, γ ) at each buildingblock: z j , α ki , β ki , γ ki with i = 1 , , ..., N ; j = 1 , , ..., N − ; k = 1 , , ..., i − . (2.1)The activation z j determines the weight of each node in the material network. For the k -th node in the i -thlayer, its weight w ki can be expressed as a summation of activations of its descendant nodes in the bottomlayer: w ki = N − i k (cid:88) j =2 N − i ( k − w jN = N − i k (cid:88) j =2 N − i ( k − a ( z j ) . (2.2)where a ( · ) is the the rectified linear unit (ReLU). Each building block takes the materials at two child nodesas inputs, and their weights determine the phase fraction of the two-layer structure. A homogenizationoperation is first performed to mix the two materials. The new homogenized material then undergoesa rotation defined by the angles ( α, β, γ ), and the rotated material is output to the mother node. Allfitting parameters of DMN have physical meanings related to the microstructural topology, and they can beoptimized through gradient-based methods by virtue of the existence of analytical solutions. Current DMN architectures are based on the tree graph, in which the physical information of any node (e.g. the stiffnesstensor and weight) are collected by only one upper-layer node during forward propagation. For other types of graphs [39, 42],extra schemes are needed to take care how the information are distributed to multiple nodes in the upper layer. .1. Design of cohesive networks The original architecture of material network is created for heterogeneous materials with perfectly bondedphases. Its extension to a material with interfacial effect (e.g. debonding) is limited by the design of thebuilding block, where the two input materials are perfectly bonded. In this regard, it is possible to considerthe interface in the original two-layer building block to be deformable. However, this naive approach maysuffer from two issues: 1) the interface orientation is predetermined by the building block, thus, loses theadaptivity; and 2) the physical meaning of a deformable interface in the upper part of the network ( i < N )is not clear since the input materials of the corresponding building block usually represent a mixture of themicroscale constituents.An alternative viewpoint regards the nodes at the bottom layer of the original network as degrees offreedom (DOF). In the online computation, each active bottom-layer node serves as an independent materialpoint, which carries the pure information of a certain microscale constituent. This provides the physicalbasis to introduce the interfaces only at the bottom layer. In the paper, it is proposed to enrich each selectedbottom-layer node with a so-called “cohesive network”, as shown in Figure 1. The index of the enrichednode and its corresponding cohesive network is q , and the total number of enriched nodes is denoted by N e .A mapping function J ( q ) finds the index of enriched node q in the bottom layer of the original materialnetwork: j = J ( q ) with q = 1 , , ..., N e . (2.3)For example, if the second cohesive network ( q = 2) enriches the third node in the bottom layer, we have J (2) = 3.Note that only active nodes in the bottom layer (with z j >
0) are necessary to be enriched. Moreover,for a two-phase material with only one type of interface, it is usually sufficient to enrich the nodes for onesingle phase. For example, in the material network for unidirectional (UD) fiber-reinforced composite withdeformable fiber-matrix interface to be tested in Section 5, only the active nodes of the fiber phase are linkedto the cohesive networks. This design will help to reduce the redundancy in the machine learning model.In the cohesive network of the enriched node q , the base material should represent the same microscaleconstituent as the one previously at bottom layer of the original material network. As shown in Figure 1,the cohesive layers (CL), indexed by p , are added to the base material sequentially through the “cohesivebuilding block” (see Section 2.3). For simplicity, all the cohesive networks are assumed to have the samenumber of cohesive layers N c . Similar to the material network, the cohesive network has its own parameters,which needs to be optimized with respect to high-fidelity RVE data using gradient-based learning algorithms.The parameters of cohesive network q are˜ z pq , ˜ α pq , ˜ β pq , ˜ γ pq with p = 1 , , ..., N c . (2.4)In this paper, the tilde symbol is used to indicate that the corresponding quantity is defined within thecohesive network. The activation of the p -th cohesive layer is ˜ z pq , and rotation angles of the p -th cohesivenetwork are ˜ α pq , ˜ β pq , ˜ γ pq . Physical meanings of all these fitting parameters will be explained in the followingtwo subsections. Evidently, an important condition that the cohesive network should satisfy is that if allthe cohesive layers are deactivated, the base material should be returned to the enriched node without anychange.The total number of parameters in the new material network with enriching cohesive networks is N total = Material network (cid:122) (cid:125)(cid:124) (cid:123)(cid:0) × N − − (cid:1) + Cohesive networks (cid:122) (cid:125)(cid:124) (cid:123) N e × N c . (2.5)In the multi-stage training strategy to be discussed in Section 3, not all parameters are fitted at the sametime. The basic idea is to train the original material network first and then the fitted parameters ( z, α, β, γ )are treated as constants in the subsequent training of the integrated network with cohesive layers. To easethe offline sampling and training processes, the fitting parameters are optimized based on linear elastic DNSdata, which is sufficient for capturing essential physical interactions among the microscale material phasesand interfaces. The learned model will be extrapolated to unknown material and loading spaces, such asplasticity, interfacial failure (debonding), and loading-unloading paths.4 = 𝐃 𝐃 , 𝐆, &𝑣, &𝛼, )𝛽, &𝛾 Material 0: 𝐃 , 𝛿𝛆 Cohesive layer:
𝐆, 𝛿𝐝 𝑇 = 1/ &𝑣 𝐑( &𝛼 , )𝛽 , &𝛾)𝛿𝛆 = 𝛿𝛆 𝛿𝛆 , 𝛿𝐝, &𝑣, &𝛼, )𝛽, &𝛾𝐑 ( &𝛼 , )𝛽 , &𝛾) Figure 2: Illustration of the cohesive building block. ˜ v is the reciprocal length parameter, T is the effective thickness of thebuilding block, and (˜ α, ˜ β, ˜ γ ) are the rotation angles. In a generic cohesive building block shown in Figure 2, a non-negative reciprocal length parameter ˜ v isintroduced for the cohesive layer to represent its size effect. Physically, one can interpret ˜ v as the inverse ofeffective thickness T of the undeformed building block in the direction normal to the cohesive layer:˜ v = 1 T . (2.6)Since the deformation of the cohesive layer is characterized by its separation displacement, this reciprocallength parameter ˜ v will be further used to transform the displacement vector to the strain tensor contributedto the building block. It is noted that when ˜ v = 0, the building block has infinite effective thickness andsees no effect from the cohesive layer, indicating that it is perfectly bonded.For the p -th cohesive layer within the q -th cohesive network, its reciprocal length parameter ˜ v pq is designedto be activated by the rectified linear unit (ReLU) a ( · ):˜ v pq = a (˜ z pq ) L = max(˜ z pq , L , (2.7)where L is the characteristic length of the microstructure, which is a predefined constant.Theoretically, L does not affect the optimum solution. However, an appropriately selected L will helpto keep the activations vary in a similar range as the rotation angles, and thus it eases the convergence ofthe gradient descent algorithm used for the training. The choice of characteristic length L depends on themorphology of the microstructure. For example, in a unidirectional fiber-reinforced composite, L can bechosen as the average diameter of the fibers. While for another polycrystalline material, L can be chosen asthe average grain size.The activation of cohesive layer ˜ z pq in Eq. (2.7) is a unitless fitting parameter to be learned. Similar tothe activation of node z in the original material network, once ˜ z pq becomes negative during the training, thecorresponding cohesive layer is deactivated, and it will not be activated again due to the vanished gradientof the ReLU function.Assume the volume of the q -th enriched node is given by Ω q . Since the cohesive layers have zero volume,all the cohesive building blocks in the q -th cohesive network share the same volume Ω q . Following thedefinition of ˜ v in Eq. (2.6, one can define the effective area of the p -th cohesive layer ˜ S pq as˜ S pq = ˜ v pq Ω q . (2.8)Note that the RVE scaling effect due to the existence of deformable interfaces can be considered byscaling the characteristic length L accordingly. For example, assume an RVE is trained to a DMN with a setof activations ˜ z under L . If the RVE geometry is enlarged linearly by a scale factor of 2 and the cohesive lawstays the same, the new DMN can be obtained by simply updating the characteristic length to 2 L withoutchanging the trained fitting parameters. According to Eq. (2.7), this is also equivalent to decreasing all thereciprocal length parameters ˜ v by a factor of 0.5. 5 .3. Cohesive building block As depicted in Figure 2, a generic cohesive building block has two inputs: a planar cohesive layer anda bulk material 0. The cohesive layer has zero thickness, and its normal is in direction 3. In general, itsconstitutive relation for an arbitrary traction-separation law can be linearized as∆ d = G ∆ t + δ d , (2.9)where ∆ d is the increment of displacement vector, and ∆ t is the increment of traction vector. The residualdisplacement vector δ d will appear in a nonlinear traction-separation law, including the bilinear one usedin this work. For a 3-dimensional problem, the traction t is in the unit of force per unit area. It has threecomponents: t n , t s , and t t , which represent the normal and the two shear tractions, respectively. Similarly,the displacements in these directions are denoted by d n , d s , and d t . Moreover, the compliance matrix G relates the traction vector to the displacement vector under a linear elastic setting, which has 6 independentcomponents in 3D due to symmetry: G = G nn G ns G nt G ss G st sym G tt . (2.10)The stiffness matrix K is defined as the inverse of G : K = G − = K nn K ns K nt K ss K st sym K tt . (2.11)Note that for history-dependent traction-separation laws, internal variables will also be carried at thecohesive layer, while the incremental constitutive relation 2 . G will be infinite. To avoid the potential numerical issues, a tiny stiffness numberis recommended to be added to the diagonal terms in K . The same practice is also adopted in the paperfor interfacial failure analysis (see Eq. (4.13), when the cohesive layer is fully separated or damaged.In a small-strain setting, the Cauchy stress σ and the infinitesimal strain ε are used as stress and strainmeasures, respectively. Moreover, the stress and strain are vectorized under the Mandel notation: σ = { σ (cid:101) , σ (cid:101) , σ (cid:101) , √ σ (cid:101) , √ σ (cid:101) , √ σ (cid:101) } T = { σ , σ , σ , σ , σ , σ } T (2.12)and ε = { ε (cid:101) , ε (cid:101) , ε (cid:101) , √ ε (cid:101) , √ ε (cid:101) , √ ε (cid:101) } T = { ε , ε , ε , ε , ε , ε } T , where the subscripts 3, 4, and 5 denote the shear directions. Equipped with these conventions, the consti-tutive relation of Material 0 can be written as∆ ε = D ∆ σ + δ ε , (2.13)where ∆ ε and ∆ σ are the strain and stress increments, respectively. D is its compliance matrix, whichis symmetric and has 21 independent entries for a 3-dimensional problem. Additionally, δ ε is the residualstrain.Three operations are introduced in the building block to propagate the mechanical information includingthe compliance matrix and residual strain/displacement. First, material 0 is rotated and the obtainedconstitutive relation is ∆ ε r = D r ∆ σ r + δ ε r . (2.14)Then rotated material 0 is homogenized together with the cohesive layer, and the new constitutive relationis ∆¯ ε = ¯ D ∆ ¯ σ + δ ¯ ε . (2.15)6 able 1: Constitutive relations in the cohesive building block. The final expressions of the compliance matrix D and theresidual strain δ ε of are given in Eq. (2.27) and (2.28) respectively. Inputs Cohesive layer ∆ d = G ∆ t + δ d Material 0 ∆ ε = D ∆ σ + δ ε Intermediate steps (i) Rotation of material 0 ∆ ε r = D r ∆ σ r + δ ε r (ii) Homogenization ∆¯ ε = ¯ D ∆ ¯ σ + δ ¯ ε Output step (iii) Rotation of homogenized material ∆ ε = D ∆ σ + δ ε To recover the initial orientation of material 0, the homogenized material is rotated back, which gives theoverall constitutive relation of the cohesive building block:∆ ε = D ∆ σ + δ ε . (2.16)The associated constitutive relations of the inputs, intermediate steps (i, ii), and output step (iii) aresummarized in Table 1.In the generic building block, the overall input-output function defined on the compliance matrix can begenerally written as D = D (cid:16) D , G , ˜ z, ˜ α, ˜ β, ˜ γ (cid:17) (2.17)and the one for the residual strain/displacement is δ ε = δ ε (cid:16) δ ε , δ d , ˜ z, ˜ α, ˜ β, ˜ γ (cid:17) . (2.18)Analytical solutions of the input-output functions can be derived following the procedure described below. (i) Rotation of material 0: The rotation operation is defined by the rotation angles (˜ α, ˜ β, ˜ γ ). Therotation matrix R can be decomposed as a product of three element rotation matrices, R (˜ α, ˜ β, ˜ γ ) = X (˜ α ) Y ( ˜ β ) Z (˜ γ ) . (2.19)The expressions of these elementary rotation matrices can be found in Appendix B. The compliance matrix D r and residual strain δ ε r of the rotated material 0 are D r = R − (˜ α, ˜ β, ˜ γ ) D R (˜ α, ˜ β, ˜ γ ) , δ ε r = R − (˜ α, ˜ β, ˜ γ ) δ ε . (2.20) (ii) Homogenization: The equilibrium condition at the interface poses that the tractions on thecohesive layer should satisfy ∆ t n = ∆¯ σ , √ t s = ∆¯ σ , √ t t = ∆¯ σ . (2.21)For the rotated material 0, the cohesive layer does not apply any displacement or stress constraints on its1, 2, and 6 directions. Along with the equilibrium condition, we have∆ σ r = ∆ ¯ σ , ∆ ε r = ∆¯ ε , ∆ ε r = ∆¯ ε , ∆ ε r = ∆¯ ε . (2.22)The contribution of displacement vector at the cohesive layer ∆ d to the homogenized strain ∆¯ ε is scaled bythe reciprocal length parameter ˜ v , which is interpreted as the inverse of effective thickness of the buildingblock in 3 direction. The unknown components of the homogenized strain can be computed as∆¯ ε = ∆ ε r + ˜ v ∆ d n , ∆¯ ε = ∆ ε r + ˜ v ∆ d s √ , ∆¯ ε = ∆ ε r + ˜ v ∆ d t √ . (2.23)7herefore, the compliance matrix and residual strain after the homogenization step are derived to be¯ D = D r + ˜ v ˜ G , δ ¯ ε = δ ε r + ˜ vδ ˜ d , (2.24)where ˜ G and δ ˜ d are the modified compliance matrix and displacement vector of the cohesive layer:˜ G = G nn G ns / √ G nt / √ G ss / G st / sym G tt / , δ ˜ d = δd n δd s / √ δd t / √ . (2.25)Although it is equivalent to derive the homogenization functions based on the stiffness matrices, in practice,the expressions based on compliance matrices are more concise since all the constraints posed by the cohesivelayer are acting on the stress components. (iii) Rotation of homogenized material: The orientation of the material 0 is recovered by an inverserotation R − (˜ α, ˜ β, ˜ γ ) on the homogenized material, and the overall compliance matrix and residual strainof the cohesive building block are D = R (˜ α, ˜ β, ˜ γ ) ¯ DR − (˜ α, ˜ β, ˜ γ ) , δ ε = R (˜ α, ˜ β, ˜ γ ) δ ¯ ε . (2.26)By combining Eq. 2.20, 2.24, and 2.26, the input-output functions of the cohesive building block arederived to be D = D (cid:16) D , G , ˜ v, ˜ α, ˜ β, ˜ γ (cid:17) = D + ˜ v R (˜ α, ˜ β, ˜ γ ) ˜ GR − (˜ α, ˜ β, ˜ γ ) (2.27)and δ ε = δ ε (cid:16) δ ε , δ d , ˜ v, ˜ α, ˜ β, ˜ γ (cid:17) = δ ε + ˜ v R (˜ α, ˜ β, ˜ γ ) δ ˜ d . (2.28)Apparently, if ˜ v = 0, the cohesive building block sees no effect from the cohesive layer.Analytical solutions under finite-strain formulation also exist, and more details about the derivation aregiven in Appendix C. However, all examples presented in this paper are based on small-strain formulation. By applying the input-output functions of the cohesive building block iteratively as illustrated in Figure1, the compliance matrix at the top of the cohesive network for the q -th enriched node with N c cohesivelayers is obtained as D N c q = D q + N c (cid:88) p =1 ˜ v pq R (˜ α pq , ˜ β pq , ˜ γ pq ) ˜ G pq R − (˜ α pq , ˜ β pq , ˜ γ pq ) , (2.29)and the residual strain is δ ε N c q = δ ε q + N c (cid:88) p =1 ˜ v pq R (˜ α pq , ˜ β pq , ˜ γ pq ) δ ˜ d pq , (2.30)where D q and δ ε q are the compliance matrix and residual strain of the base material of the q -th cohesivenetwork, respectively. Meanwhile, ˜ G pq and δ ˜ d pq are modified from the compliance matrix and displacementvector of the p -th cohesive layer based on Eq 2.25.The information at the top node of the cohesive network is transferred to the corresponding enrichednode in the bottom layer ( i = N ) of the original material network. As the stiffness tensor C and residualstress δ σ are primarily used in the material network, the following operations are performed: C J ( q ) N = (cid:16) D N c q (cid:17) − , δ σ J ( q ) N = − (cid:16) D N c q (cid:17) − δ ε N c q . (2.31)The definition of the mapping function J ( q ) is provided in Eq. 2.3.8 ixed 𝑧, 𝛼, 𝛽, 𝛾 𝐂 ’( 𝐂 ’) 𝐂 ’’ 𝐂 ’* 𝐂 )( 𝐂 )) 𝐂 (( 𝐃 (, 𝐃 (( 𝐃 () 𝐃 (’ 𝐃 ), 𝐃 )( 𝐃 )) 𝐃 )’ 𝐆 (( 𝐆 () 𝐆 (’ 𝐆 )( 𝐆 )) 𝐆 )’ 𝐂 .( 𝐊 𝐂 .( 𝐊 𝐂 .) 𝐂 .) ̅𝐂 Input: 𝐂 .( , 𝐂 .) , 𝐊 Output: Fitting Parameters: ̃𝑧, 7𝛼, 8𝛽, 7𝛾
Training stage IITraining stage I
Input: 𝐂 .( , 𝐂 .) Output: Fitting Parameters: 𝑧, 𝛼, 𝛽, 𝛾
Varying 𝑧, 𝛼, 𝛽, 𝛾 𝐂 ’( 𝐂 ’) 𝐂 ’’ 𝐂 ’* 𝐂 )( 𝐂 )) 𝐂 (( ̅𝐂 𝐂 .) 𝐂 .) 𝐂 .( 𝐂 .( Trained 𝑧, 𝛼, 𝛽, 𝛾
Figure 3: Illustration of the multi-stage training strategy for deep material network with cohesive layers. The number of layersis N = 3. The number of enriched nodes is N e = 2, and the number of cohesive layers in each cohesive network is N c = 3. Themapping function is J (1) = 1 , J (2) = 3. The inputs of training stage I are the stiffness matrices of microscale phases C p and C p . The extra input of training stage II is the interfacial stiffness matrix K c .
3. Multi-stage offline training strategy
A two-stage training strategy is proposed in this work. The first stage learns the essential phase topologywith no interfacial effect, while the second stage learns the length scale and orientation of the cohesivelayers. In comparison to learning all the fitting parameters in one single step, the two-stage strainingstrategy decouples the interfacial effect from phase interactions, so that it has fewer parameters at eachstage. Overall, it makes the landscape of the cost function less complicated for the learning algorithms toconverge.Figure 3 illustrated this strategy for a shallow material network with cohesive layers. Its number oflayers N , number of enriched nodes N e , number of cohesive layers in each cohesive network N c , and mappingfunction J ( q ) between indices of the cohesive network and bottom-layer node in the original material networkare listed below: N = 3 , N e = 2 , N c = 3 , J (1) = 1 , J (2) = 3 . (3.1)Since the offline training processes are based on linear elastic data, the material network only propagatesthe information of compliance and stiffness matrices.As shown in the figure, the material network without the enrichment of cohesive networks is first trainedfor the RVE with perfectly bonded interfaces. The overall stiffness matrix of the multi-layer network ¯ C rve for a two-phase RVE can be written as a function of the stiffness matrix from each material phase ( C p and C p ) and the fitting parameters ( z, α, β, γ ):Training stage I: ¯ C rve (cid:124) (cid:123)(cid:122) (cid:125) Output = f ( Fitting parameters (cid:122) (cid:125)(cid:124) (cid:123) z j =1 , ,..., N − , α k =1 , ,..., i − i =1 , ,...,N , β ki , γ ki ; C p , C p (cid:124) (cid:123)(cid:122) (cid:125) Inputs ) . (3.2)At training stage II, each active bottom-layer node of phase 1 is linked to a cohesive network to form anew network for the RVE with deformable interfaces. Therefore, the stiffness matrix of the cohesive layers K c enters as an extra input. Note that the fitting parameters of the two-layer building blocks obtained in9he first training stage will be kept constant. With (˜ z, ˜ α, ˜ β, ˜ γ ) as the new fitting parameters, the overallstiffness matrix the multi-layer network ˜ C rve can be written asTraining stage II: ¯ C rve (cid:124) (cid:123)(cid:122) (cid:125) Output = ˜ f ( Fitting parameters (cid:122) (cid:125)(cid:124) (cid:123) ˜ z p =1 , ,...,N c q =1 , ,...,N e , ˜ α pq , ˜ β pa , ˜ γ pq ; C p , C p , K c (cid:124) (cid:123)(cid:122) (cid:125) Inputs , z j , α ki , β ki , γ ki ) (cid:124) (cid:123)(cid:122) (cid:125) Constants . (3.3)For both stages, the cost function J is defined upon the mean square error (MSE): J = 12 N s N s (cid:88) s =1 || ¯ C dnss − ¯ C rve || || ¯ C dnss || , (3.4)where ¯ C dnss is the overall stiffness matrix of sample s computed from DNS. N s is the total number oftraining samples. The operator || ... || denotes the Frobenius matrix norm. While for training stage I, anextra regularization term is added to the cost function to constrain the magnitude of { z j } and makes theoptimization problem well-posed [41].Stochastic gradient descent (SGD) is used for minimizing the cost function, while the gradients arecalculated through back-propagation algorithm enabled by the analytical solutions derived in Section 2.3.In the rest of this section, details on the sampling and training of each stage are provided. The inputs C p and C p are assigned to the nodes in the bottom layer N in a scheme described as below, C l − N = C p , C lN = C p with l = 1 , , ..., N − . (3.5)This scheme is applicable to two-phase materials. The material network based on binary-tree structurecan also handle materials with more than two phases, however, a different assigning scheme needs to bedeveloped.In the sampling of inputs, both material phases are assumed to be orthotropically elastic, whose compli-ance matrices can be written as D pi = (cid:0) C pi (cid:1) − = /E pi − ν pi /E pi − ν pi /E pi /E pi − ν pi /E pi sym 1 /E pi / (2 G pi ) 1 / (2 G pi ) 1 / (2 G pi ) , (3.6)where pi is p p
2. To initiate material anisotropy, the tension moduli of each phase are first randomlysampled: log ( E pi ) , log ( E pi ) , log ( E pi ) ∼ U [ − , , (3.7)where U represents the uniform distribution. The moduli of phase 2 will be rescaled to create phase contrasts,and the rescaling factor ¯ E p is sampled fromlog ( ¯ E p ) ∼ U [ − , . (3.8)With the rescaling factor, the tension moduli of each phase are updated, E p kk ← ¯ E p ( E p E p E p ) / E p kk with k = 1 , , . (3.9)10fter all the tension moduli are determined, the shear moduli are sampled from G pi (cid:113) E pi E pi , G pi (cid:113) E pi E pi , G pi (cid:113) E pi E pi ∼ U [0 . , . . (3.10)To guarantee the compliance matrices are positive definite, the Poisson’s ratios are sampled from ν pi (cid:113) E pi /E pi ∼ U (0 . , . , ν pi (cid:113) E pi /E pi ∼ U (0 . , . , ν pi (cid:113) E pi /E pi ∼ U (0 . , . . Monte Carlo method is used to sample the input space, which has in total 13 independent random variables.At the start of training, the fitting parameters ( z, α, β, γ ) are initialized randomly following uniformdistributions. During the training process, the network can be compressed through two operations: 1)Deletion of the parent node with only one child node; 2) subtree merging based on the similarity search [41].
In training stage II, the active nodes of phase 1 in the bottom layer are enriched by the cohesive network,while the input C p is still assigned to the other active nodes of phase 2: C lN = C p with l = 1 , , ..., N − . (3.11)Meanwhile, the input C p is assigned to the base material of each cohesive network, and all the cohesivelayers receive the same stiffness matrix K c : D q = (cid:0) C p (cid:1) − , G pq = ( K c ) − with q = 1 , , ..., N e ; p = 1 , , ..., N c . (3.12)The samples generated from training stage I are reused and augmented with new elastic cohesive proper-ties. The cohesive layer is assumed to be isotropic in the two shear directions. In addition, the off-diagonalterms in K c are set to zero. As a result, K c has only two independent components K cnn and K css : K c = ( G c ) − = K cnn K css
00 0 K css . (3.13)The elastic constant in the normal direction K cnn is sampled aslog ( K cnn × L ) ∈ U [ − , , (3.14)where L is the characteristic length of the microstructure, and a discussion about the choice of L can befound in Section 2.2. For the elastic constant in the shear direction K css , we havelog ( K css K cnn ) ∈ U [ − , . (3.15)Although the new input space has 15 independent variables, only two of them need to be sampled at thisstage and are further appended to the previous space. Similarly, Monte Carlo method is used for thesampling.The fitting parameters (˜ z, ˜ α, ˜ β, ˜ γ ) are initialized randomly, while the parameters ( z, α, β, γ ) trained fromstage I are fixed during stage II, as shown in Figure 3.Model compression can be introduced to reduce the redundancy of the cohesive layers. By looking atEq. (2.29) and (2.30), one can see that the terms contributed by each cohesive layer are linearly added tothe overall compliance matrix or residual strain. This nice property enables the merge of two cohesive layersregardless of their order in the network. 11n each cohesive network, the similarity between any pair of cohesive layers will be evaluated. Let usassume that their fitting parameters are given by (˜ z p , ˜ α p , ˜ β p , ˜ γ p ) and (˜ z p (cid:48) , ˜ α p (cid:48) , ˜ β p (cid:48) , ˜ γ p (cid:48) ). Since the cohesivelayers are set to behave isotropically in the shear directions, they are said to be similar if they share thesame normal direction. Mathematically, it requires that the absolute value of the 33 component in a matrix R p - p (cid:48) is close to 1: | R p - p (cid:48) | ≈ , with R p - p (cid:48) defined as R p - p (cid:48) = R (˜ α p , ˜ β p , ˜ γ p ) R − (˜ α p (cid:48) , ˜ β p (cid:48) , ˜ γ p (cid:48) ) . (3.16)In this case, the two cohesive layers will be merged, and their new activations are˜ z p ( new ) = ˜ z p + ˜ z p (cid:48) , ˜ z p (cid:48) ( new ) = 0 . (3.17)The p (cid:48) -th cohesive layer is deactivated, so that the size of the machine-learning model is reduced.
4. Cohesive law for online interfacial failure analysis
In the online stage, the DMN with cohesive layers, trained based on linear elastic data, can be extrap-olated to consider nonlinear interfacial behavior, such as the failure at the material interface. Specifically,the irreversible cohesive law proposed by Camacho and Ortiz [16, 17] is adopted in this work. Simplificationfirst arises from the assumption that the cohesive surface is isotropic. This indicates that the resistance tosliding is independent of the direction of sliding. and only depends on the magnitude d S : d S = (cid:113) d s + d t . (4.1)The total opening displacement vector d and traction t take the following forms, d = d n n + d S , t = t n n + t S , (4.2)where n is the unit normal to the cohesive surface. An effective opening displacement d m is further introducedto simplify the mixed-mode cohesive law, and the free energy density per unit undeformed area φ is a functionof d m , φ = φ ( d m , q ) , (4.3)where q are some internal variables which describe the irreversible processes. Accordingly, the effectivetraction is defined as t m = ∂φ∂d m ( d m , q ) . (4.4)The compressive stress state should not cause the cohesive layer to fail, so that it does not contributeto the free energy density. In addition, no friction effect is included in the cohesive model. The tensile andcompressive cases are considered separately:1. For the tensile case d n ≥
0, the effective opening displacement d m is d m = (cid:113) d n + β d S , (4.5)where the positive parameter β defines the ratio of effects from normal and shear displacements. Thecohesive law becomes t = ∂φ∂ d = t m d m (cid:0) d n n + β d S (cid:1) . (4.6)12 " , 𝜎 " 𝑑 & , 0 𝑑 ’ 𝑡 ’ 𝑡 ’) = 𝐾 𝑑 , , 𝜎 , 𝐺 " = 12 𝜎 " 𝑑 & Figure 4: The bilinear cohesive law expressed in terms of the effective opening displacement d m and the effective traction t m .The initial stiffness is K , and the critical energy releasing rate is denoted by G c .
2. For the compressive case d n <
0, the effective opening displacement d m is d m = β | d S | . (4.7)The cohesive law can be written as t = t n n + ∂φ∂ d S = t n n + t m d m β d S . (4.8)As shown in Figure 4, a bilinear cohesive law is adopted in this work. In the elastic loading regime, themodulus of cohesive layer is K . In the application to be shown in Section 5, this interfacial modulus is chosento be a large value compared to the moduli of other material phases, thus, the interface is nearly perfectlybonded at the beginning. The effective opening displacement before the softening regime is denoted by d c ,and the corresponding effective traction is σ c . The opening displacement at full failure t m = 0 is d f . Thecritical energy releasing rate G c is equal to the area under the curve of t m - d m , G c = 12 σ c d f . (4.9)The only internal variable is the effective opening displacement at maximum traction d . Initially, d isequal to d c , d = d c , (4.10)and it should decrease along the softening path when d m exceeds the current d . Mathematically, theeffective traction t m of the bilinear cohesive law can be calculated by t m = d m − d f d c − d f σ c if d m = d and ˙ d m ≥ , (4.11)and t m = d m d σ if d m < d or ˙ d m < . (4.12)In practice, to avoid singularity at full separation of cohesive layer, a linear term κd m is added to the effectivetraction, t m ← t m + κd m (4.13)where κ is a tiny stiffness parameter chosen as κ/K = 1 × − in this work.When reloading, the curve will simply trace the unloading path and rejoin the cohesive law at d m = d .Meanwhile, an elastic response is assumed for the compressive case, with the same modulus as the one inthe bilinear law, t n = Kd n if d n < . (4.14)More details on the derivation of the tangent stiffness matrix and viscous regularization for implicitanalysis in the online stage are provided in Appendix D.13 ohesive elements123 (a) Geometry and mesh. -
0. 600. 00 S t r a i n Time (s) (b) Loading path.Figure 5: Unidirectional fiber-reinforced RVE: (a) Volume fraction of the fiber phase is 29.4%, the FE model has 69048 nodes,63720 8-node hexahedron elements, and 4320 8-node zero-thickness cohesive elements; (b) History of strain component in theonline stage includes both loading and unloading paths. Snapshots of deformed RVEs are taken at t=0.04, 0.12, 0.24, 0.30 s.Table 2: Material parameters of the fiber phase, the matrix phase and the interfaces used in the online stage.
Fiber E (GPa) ν
500 0.3Matrix E (GPa) ν σ y (GPa) Hardening100 0.3 0.1 Eq. (5.2)Interface K (GPa/m) σ c (GPa) G c (GPa · m) β ζ (GPa · s)1 × . × − × −
5. Application to unidirectional fiber-reinforced composite
Here the deep material network with cohesive layers is applied to describe a two-phase unidirectional(UD) fiber-reinforced composite. The method is general, and another example of the particle-reinforcedcomposite is briefly discussed in Appendix A. The geometry of the UD RVE is shown in Figure 5 (a). Thevolume fraction of the fiber phase is 29.4%. The average diameter of the fibers is around 2.5 mm, so thatthe characteristic length L , a hyper-parameter for the cohesive networks, is set to be L = 2 . . (5.1)In the FE mesh for DNS, there are 69048 nodes, 63720 8-node hexahedron solid elements, and 43208-node zero-thickness cohesive elements. All the material parameters used in the online stage are given inTable 2.The fiber phase is assumed to be linear elastic, and its Young’s modulus E and Poisson’s ratio ν are 500GPa and 0.3, respectively. Two cases are considered for the matrix phase. In the first case, the matrix phaseis linear elastic, with E = 100 GPa and ν = 0 .
3. In the second case, the matrix phase is an elasto-plasticmaterial, which obeys von Mises plasticity with piecewise linear hardening. The material yields at 0.1 GPa,and the yield stress σ Y is a function of the equivalent plastic strain ε p : σ Y ( ε p ) = (cid:40) . ε p ε p ∈ [0 , . . ε p ε p ∈ [0 . , ∞ ) GPa . (5.2)The equivalent modulus ( K × L =25000 GPa) of the interface is much larger than the moduli of fiberand matrix phases, indicating that the interface is nearly perfectly bonded at the beginning. Following Eq.14 able 3: Training results of UD RVE at stages I and II. Average training error ¯ e tr , average test error ¯ e tr and maximum testerror are provided for each DMN. In stage II, the number of cohesive layers in each cohesive network is N c = 4. Stage I Stage II
Training ¯ e tr Test ¯ e te Maximum e tes Training ¯ e tr Test ¯ e te Maximum e tes N = 5 11.1% 10.5% 47.2% 11.9% 11.3% 47.5% N = 7 0.79% 0.76% 3.46% 0.88% 0.91% 3.82% N = 9 0.23% 0.25% 1.34% 0.31% 0.33% 2.45%(4.9), the effective opening displacement at complete failure d f can be computed from the critical effectivetraction σ c and the critical energy releasing rate G c : d f = 2 G c σ c = 0 .
01 mm . (5.3)To restrain the rate-dependency from the numerical stabilization in Eq. D.7, the viscosity-like parameter ζ (see the definition in Eq. (D.7)) is set to a small value.For both training stages, there are 400 training samples and 100 test samples. For each sample, a setof linear elastic DNS of the UD RVE under 6 orthogonal loading conditions is performed, which typicallytakes 182 s using LS-DYNA R (cid:13) RVE package on 10 Intel R (cid:13) Xeon R (cid:13) CPU E5-2640 v4 2.40 GHz processors .With 1000 samples in total, the offline data generation process costs 50 h. This time can be further reducedby distributing the jobs to more processors.In the SGD algorithm, the mini-batch size is chosen as 20, so that there are 20 learning steps in eachepoch. Additionally, a Bold driver is applied to dynamically adjust the learning rate during the training.The DMN framework integrated offline sampling, definition of the DMN data structure, offline training (e.g.backpropagation), and online prediction. Currently, it is implemented in Python, with necessary librariesfor matrix manipulation and linear algebra. In training stage I, the original material networks with three different depths, N = 5 , , and 9, areconsidered. For N = 5 and 7, the networks were trained for 20000 epochs. Since the network with N = 9has more fitting parameters, it was trained for 40000 epochs. The training and test errors are listed in Table3, where the error of sample s is defined as e s = || ¯ C dnss − ¯ C rves |||| ¯ C dnss || . (5.4)To show the fractions of active nodes in the bottom layer, treemap plots of all three DMNs are alsoprovided in Figure 6. For the network with N = 5, there remains only one active node in the bottom layerfor the fiber phase. Evidently, this is not sufficient for capturing the anisotropic behaviors of the UD RVE,and the average test error ¯ e te is more than 10%. One important feature of DMN is that the average trainingerror and the average test error of unseen samples are almost the same as shown in Table 3. With a smallamount of data, it requires no regularization scheme (e.g. L2 penalty or Dropout) to avoid the overfittingissue, mostly thanks to the embedded physics in the building block. More detailed results of training stageI, including the error histories, can be found in [41].In training stage II, the trained material networks from stage I are enriched with cohesive networks, wherethe new fitting parameters are also optimized using the SGD algorithm. For all the cases, the networks weretrained for 5000 epochs. Figure 7 (a) shows the histories of test error for N = 5, 7, and 9, with the same The same machine is used for all other computations in this work. a) N = 5 , N a = 4 (b) N = 7 , N a = 14 (c) N = 9 , N a = 60Figure 6: Treemaps of trained material networks from training stage I. The number of active nodes in the bottom layer N a isalso provided. The dark blocks belong to the fiber phase. Epoch − − T e s t e rr o r N = 5 N = 7 N = 9 (a) Different N , same N c = 4. Epoch − − T e s t e rr o r N c = 1 N c = 2 N c = 4 N c = 6 (b) Different N c , same N = 9.Figure 7: Histories of test error in training stage II. All the networks are trained for 5000 epochs. number of cohesive layers in each cohesive network N c = 4. Their training and test errors at the end oflearning are also provided in Table 3. Apparently, the network with N = 9 offers the most accurate results.The effect of the number of cohesive layers at each node N c on the training results has also beeninvestigated for depth N = 9. We can see from Figure 7 (b) that one cohesive layer N c = 1 is not able tocapture the interfacial effect accurately. While for N c = 2, 4, and 6, the test errors are all reduced to beless than 0.5%. After 5000 epochs, the network with N c = 4 reaches the lowest test error, instead of the onewith N c = 6. This can be explained by the fact that the network with N c = 6 has more redundant cohesivelayers, and the hyper-parameter L = 25 mm may not be optimal anymore.It is also interesting to compare the errors from training stage I and II in Table 3. The error at stageI, in any measure, is always less than the one in stage II. This is physically sound as the UD RVE withdeformable fiber-matrix interfaces is a more complex material system. Moreover, train stage II inherits theparameters ( z, α, β, γ ) from stage I, and they are fixed in the training. Following this strategy, the accuracyof enriched network is somehow limited by the original material network.Figure 8 shows the scatter plot of the cohesive layers for N c = 4. For each cohesive layer, the x value isthe volume fraction of its enriched node in the bottom layer of original material network, and the y valueis its reciprocal length parameter ˜ v . The remaining number of active cohesive layers for N = 5 , , and 9are 2, 13, and 40. Since the network with N = 9 has more enriched nodes, the mean volume fraction ofthe nodes is smaller, and more cohesive layers are deactivated. For all the cases, the mean reciprocal lengthparameter is around 0.05 mm − . This value is related to the characteristic length of the UD RVE, whoseaverage diameter of the fibers is 25 mm. 16 − − Volume fraction of enriched node . . . . . . . . R e c i p r o c a ll eng t h ˜ v ( mm − ) Deactivated cohesive layers N = 5 N = 7 N = 9 Figure 8: Scatter plots of the cohesive layers for N = 5, 7, 9 and N c = 4 after 5000 epochs of training. The y-axis shows thereciprocal length parameter ˜ v of a cohesive layer, and the x-axis shows the volume fraction of the corresponding bottom-layernode enriched by the cohesive layer. The deactivated cohesive layers have ˜ v = 0. Offline training time:
All training processes are parallelized using 10 processors. Here the computationaltimes in stage II for networks with different depths are listed and compared. For N = 5, networks with N c = 2 , , N = 7, the training times are 1.7, 2.0, 2.2 hours. For N = 9, thetraining times are 6.7, 7.9, 8.2 hours. Therefore, N c has less effect on the training time than N , indicatingthat most of the computational cost still locates at the original material network. For the UD composite, three different loading directions are evaluated in the online stage: transversetension and compression, transverse shear, and longitudinal shear. Here the case under longitudinal tensionalong the fiber direction is not presented, as there is almost no interfacial failure within the loading range.For each loading case, the overall strain component in the loading direction follows the path shown in Figure5 (b), while the overall stress components in the other directions are zero. For example, in transverse tensionand compression, the loading conditions can be written as ε = (cid:40) t t ∈ [0 , . . − t t ∈ [0 . , . , σ = σ = σ = σ = σ = 0 . (5.5)Accordingly, the absolute strain rate | ˙ ε | is 1.0 s − . The time step in both DNS and DMN is ∆ t = 1 . × − .For simplicity, ε and σ are used to denote the overall strain and stress in this section. Snapshots of DNSare taken at t = 0 . , . , . , and 0.030 s to illustrate the evolution of RVE deformation during theloading.For all the online examples to appear in this section, the study shows that N c has little effect on thestress-strain predictions as long as N c ≥
2, while the depth of the original material network N is moreprominent. Therefore, N c is set to be 4 by default in the online DMNs. Ideally, a 3D material block can befully covered by 3 cohesive layers if they are oriented in three orthogonal directions. In this regard, N c = 4tends to be sufficient for the solution space to capture the interfacial effects, and more layers would not helpto improve the accuracy too much. This statement on the choice of N c also applies to the particle-reinforcedRVE shown in Appendix A. However, it will be investigated more coherently in the future.Figure 9 summarizes the results of transverse tension and compression. As shown in Figure 9 (a) forplastic matrix, no interfacial failure happens at t = 0 .
004 s. Since the fiber phase is stiffer than the matrixphase, its average von Mises stress is higher. The interfaces normal to the loading direction start to debondat t = 0 .
005 s and reach the maximum opening at t = 0 .
012 s. Due to the debonding, less load transfersfrom the matrix to the fibers, and more stress concentration appears in the matrix phase. Around t = 0 . = 0.004 s 𝑡 = 0.012 s 𝑡 = 0.024 s 𝑡 = 0.030 s von M i s e s s t r e ss ( G P a ) (a) DNS snapshots of deformed RVE with plastic matrix. The displacements are magnified 10 times. − .
005 0 .
000 0 .
005 0 . ε − . − . − . − . . . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (b) σ vs. ε , elastic matrix. − .
005 0 .
000 0 .
005 0 . ε − . − . − . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (c) σ vs. ε , plastic matrix.Figure 9: Results of UD RVE with debonding interfaces under transverse tension and compression. Two cases are consideredfor the matrix phase: (b) Linear elasticity and (c) von Mises plasticity. Each cohesive network has 4 cohesive layers: N c = 4. s, closure of the interfaces happens. Interestingly, further compression results in debonding of the interfacesin the loading direction, which is related to the Poisson’s effect. The corresponding stress-strain curves fromDNS and material networks with N = 5 , , N c = 4 are provided in 9 (c). Stiffness degradation due tothe failure of interfaces is observed during unloading. Other softening behaviors are also captured by DMNwith cohesive layers. For N = 7 and 9, the DMN predictions agree with the corresponding DNS results verywell.The difference between transverse tension and compression is more evident on the plot for elastic matrixshown in Figure 9 (b). Since there is no plasticity in the matrix, during unloading, the stress-strain curvealmost goes back to the origin with a degraded stiffness. However, in compression, the interfaces are incontact, so that the stiffness returns to a high value and the response is linear elastic without any failure. Itappears that the DMN slightly underestimates the stiffness under compression. Nevertheless, these physicalphenomena related to the interface debonding and closure are all captured by DMNs.Figure 10 summarizes the results from transverse shear loading. The DNS snapshot at t = 0 .
012 s showsthat the debonding of interfaces is driven by the normal traction. After reversing the loading direction,previously debonded interfaces are closed, while new surfaces are created at different locations. As one cansee from Figure 10 (b), the reverse-loading part of stress-strain curve still undergoes some softening effectduring the plastic yielding, which is consistent with those findings in the DNS snapshots. For either casewith elastic or plastic matrix, the DMN with N = 9 is able to capture the overall stress-strain response,while the amount of interfacial failure is over-predicted. The network with N = 5 has a bad performance inpredicting the transverse shear behavior. Hence, representing the fiber phase with only one node in DMN(see Figure 6 (a)) is not sufficient for the debonding analysis of UD RVE.Figure 11 summarizes the results from longitudinal shear loading. The snapshots show that the interfacialfailure is driven by the sliding deformation, which is different from the previous two loading cases. Moreover,18 = 0.004 s 𝑡 = 0.012 s 𝑡 = 0.024 s 𝑡 = 0.030 s von M i s e s s t r e ss ( G P a ) (a) DNS snapshots of deformed RVE with plastic matrix. The displacements are magnified 10 times. − .
005 0 .
000 0 .
005 0 . ε − . . . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (b) σ vs. ε , elastic matrix. − .
005 0 .
000 0 .
005 0 . ε − . − . − . . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (c) σ vs. ε , plastic matrix.Figure 10: Results of UD RVE with debonding interfaces under transverse shear. Two cases are considered for the matrixphase: (b) Linear elasticity and (c) von Mises plasticity. Each cohesive network has 4 cohesive layers: N c = 4. after reversing the loading direction at t = 0 .
012 s, the sliding locations still stay the same, indicating thatthe amount of newly failed interfaces is small. This is reflected in the corresponding stress-strain curveplotted in Figure 11 (c), where the overall RVE behavior is elasto-plastic with negligible softening effect.Similarly, in Figure 11 (b) for elastic matrix, the stress-strain curve after reversing loading direction is almostlinear elastic, and passes the origin. For both cases, the DMN predictions with N = 7 and 9 agree with theDNS results very well. Online prediction time:
The DMN with cohesive layers learns an efficient reduced-order model of theDNS RVE with deformable interfaces. Both DNS and DMN models have 300 loading steps. A typical DNSwith plastic matrix took 5.6 h on 10 processors. Yet it is noteworthy that the original DNS with perfectlybonded interfaces (or no cohesive element) only took 1.6 h on 10 CPUs, implying that the softening interfacesconsume more iterations in the implicit analysis. In comparison, DMNs with cohesive layers took 1.4, 5.8and 33.3 s on one CPU for with N = 5 , N = 9 and N c = 4 is more than 6000 times faster than the DNS. Local properties can also be predicted by the DMN. In the online stage, each active node in the bottomlayer of DMN is treated as an independent microscale material, representing either the fiber or the matrixphase in the UD composite. Meanwhile, each cohesive layer can be regarded as a part of the fiber-matrixinterfaces, and its effective area is defined in Eq. (2.8), which is related to the reciprocal length parameter.With the volume fractions of the bottom-layer nodes and the effective areas of the cohesive layers, onecan reproduce the probability distribution of local fields in the 3D bulk material and the 2D interfaces,respectively. In this section, the DMN with N = 9 , N c = 4 is taken as an example, which has 60 active19 = 0.004 s 𝑡 = 0.012 s 𝑡 = 0.024 s 𝑡 = 0.030 s von M i s e s s t r e ss ( G P a ) (a) DNS snapshots of deformed RVE with plastic matrix. The displacements are magnified 10 times. − .
005 0 .
000 0 .
005 0 . ε − . . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (b) σ vs. ε , elastic matrix. − .
005 0 .
000 0 .
005 0 . ε − . − . − . . . . σ ( G P a ) N = 5 N = 7 N = 9 DNS (c) σ vs. ε , plastic matrix.Figure 11: Results of UD RVE with debonding interfaces under longitudinal shear. Two cases are considered for the matrixphase: (b) Linear elasticity and (c) von Mises plasticity. Each cohesive network has 4 cohesive layers: N c = 4. . . . . . . F i be r- ρ DMNDNS . . . . . . von Mises stress (GPa) M a t r i x - ρ DMNDNS (a) von Mises stress in fiber and matrix phases. . . . . . . . Normal traction (GPa) I n t e r f a c e - ⇢ DMNDNS (b) Normal traction at interface.Figure 12: Local fields predicted by DMN and DNS under uniaxial tension at t = 0 .
004 s and ε = 0 . bottom-layer nodes and 40 active cohesive layers after the training. While for the DNS, there are 63720solid elements and 4320 cohesive elements.Figure 12 shows the distributions of various local fields predicted by DMN and DNS with plastic matrixunder uniaxial tension at t = 0 .
004 s, while the strain in the loading direction is ε = 0 . . . . . . . F i be r- ρ DMNDNS . . . . . . von Mises stress (GPa) M a t r i x - ρ DMNDNS (a) von Mises stress in fiber and matrix phases. . . . . . . . Normal traction (GPa) I n t e r f a c e - ⇢ DMNDNS (b) Normal traction at interface.Figure 13: Local fields predicted by DMN and DNS under uniaxial tension at t = 0 .
012 s and ε = 0 . ing global stress-strain curves are given in Figure 9 (c). For the matrix and fiber phases, the distributionsof von Mises stress are compared in Figure 12 (a). Since most of the interfaces are still bonded, the averagevon Mises stress of the fiber phase is more than two times larger than that of the matrix phase. For theinterfaces, the normal traction ranges from -0.15 GPa to 0.45 GPa, which is less than the critical effectivetraction σ c = 0 . t = 0 .
012 sand ε = 0 .
012 are presented in Figure 13. Due to debonding, the profile of the normal traction in Figure13 (b) has a large density at 0 GPa, which corresponds to those failed interfaces in DNS or cohesive layersin DMN. As expected, the maximum normal traction is still less than the critical effective traction 0.5 GPa,while the compressive traction reaches -0.55 GPa. The distribution of von Mises stress in the matrix phasebecomes more dispersed, and the average stress is close to the one in the fiber phase. Although the DMNhas much less DOF, it gives a good estimation on the distributions of local fields comparing to the DNSreference.
6. Conclusions and future work
A new architecture of DMN with cohesive layers is proposed to consider deformable interfaces inside aheterogeneous material. In the physics-based cohesive building block, the activation of a cohesive layer isgoverned by a reciprocal length parameter, which also defines the effective thickness of the building blockin its normal direction. Due to the existence of analytical solutions, the new material network can stillbe trained using gradient-based optimization methods and the backpropagation algorithm. Moreover, atwo-stage training strategy is utilized to learn the fitting parameters in the original material network andthe cohesive networks separately.The predictive capability of the method is demonstrated through the debonding analysis of the UDcomposite under various loading conditions, including transverse tension and compression, transverse shear,and longitudinal shear. The behavior of the fiber-matrix interface is described by a mixed-mode irreversiblecohesive law. Though the material network is trained with linear elastic data, it can accurately capture thenonlinear RVE responses both locally and globally with complexity coming from the plastic matrix and thesoftening interfaces. Involving less DOF, the enriched material network with N = 9 and N c = 4 is morethan 6000 times faster than the DNS in terms of the CPU time.The proposed method opens new possibilities in materials engineering and design. First, DMN offers anefficient way of realizing concurrent multiscale simulations. With the addition of cohesive layers, it providesan effective data-driven material model for the interfacial failure analysis of a large-scale heterogeneousstructure. Since measuring interface properties directly from the experiments can sometimes be challenging,21 able A.4: Training results of the particle-reinforced RVE at stages I and II. Average training error ¯ e tr , average test error ¯ e tr and maximum test error are provided for each DMN. In stage II, the number of cohesive layers in each cohesive network is N c = 4. Stage I Stage II
Training ¯ e tr Test ¯ e te Maximum e tes Training ¯ e tr Test ¯ e te Maximum e tes N = 4 7.61% 7.79% 17.9% 8.12% 8.40% 17.9% N = 6 1.34% 1.39% 4.46% 1.47% 1.49% 2.81% N = 8 0.53% 0.59% 2.41% 0.64% 0.67% 2.47% von M i s e s s t r e ss (a) DNS snapshots at ε = 0 . − .
005 0 .
000 0 .
005 0 . ε − . − . − . . . . σ ( G P a ) N = 4 N = 6 N = 8 DNS (b) σ vs. ε , plastic matrix.Figure A.14: Particle-reinforced RVE: (a) The volume fraction of the particle phase is 22.6% while the FE model has 84693nodes, 59628 10-node tetrahedron elements, and 9647 6-node zero-thickness cohesive elements; (b) Stress-strain curves underuniaxial tension and compression. DMN can also be used for inverse prediction of those unknown properties. Finally, with the transfer learningtechnique proposed in [24], the unified microstructural databases can potentially be created for the designof various modern material systems with nontrivial interfacial effects.
Acknowledgment
Z. Liu would like to acknowledge Dr. C.T. Wu, Dr. Yong Guo, Dr. Tadashi Naito, and Tianyu Huangfor helpful discussions. Z. Liu warmly thanks Dr. John O. Hallquist for his support to this research. Thisis a preprint version of an article published in Computer Methods in Applied Mechanics and Engineering.The final authenticated version is available online at: https://doi.org/10.1016/j.cma.2020.112913 . Appendix A. Application to particle-reinforced composite
Other than the UD composite presented in the main content, the deep material network with cohesivelayers is applied to a particle-reinforced composite. Volume fraction of the particle phase is 22.6%. In termsof the DNS, the FE model has 84693 nodes, 59628 10-node tetrahedron elements, and 9647 6-node zero-thickness cohesive elements. The same material properties in Table 2 are assigned to the new RVE, while22ts particle phase shares the same properties with the fiber phase of the UD composite. The characteristiclength L is set to be 5.0 mm, which is same as the diameter of the particles.In training stage I, three material networks are evaluated with different depths N = 4, 6, and 8, whichare trained for 20000, 20000, and 40000 epochs, respectively. In training stage II, the number of cohesivelayers in each cohesive network is N c = 4, and all three enriched DMNs are trained for 5000 epochs. Thetraining and test errors of both stages are listed in Table A.4. With depth N = 8, the average training andtest errors of both stages I and II are reduced to be less than 1%.The results under uniaxial tension and compression are provided in Figure A.14. The loading historyof ε is shown in Figure 5 (b), while the RVE is stress-free in the other directions. As shown in the DNSsnapshot at ε = 0 . N = 6 and 8. Appendix B. Elementary rotation matrices
In small-strain formulation, the elementary rotation matrices shown in Eq. (2.19) are given byX (1 , = 1 , X ([2 , , , [2 , , ( α ) = r p ( α ) , X ([5 , , [5 , ( α ) = r v ( α ); (B.1)Y (2 , = 1 , Y ([1 , , , [1 , , ( β ) = r p ( − β ) , Y ([4 , , [4 , ( β ) = r v ( − β );Z (3 , = 1 , Z ([1 , , , [1 , , ( γ ) = r p ( γ ) , Z ([4 , , [4 , ( γ ) = r v ( γ ) . The in-plane and output-plane rotation matrices r p and r v for an arbitrary angle θ are defined in Mandelnotation as r p ( θ ) = cos θ sin θ √ θ cos θ sin θ cos θ −√ θ cos θ −√ θ cos θ √ θ cos θ cos θ − sin θ , r v ( θ ) = (cid:26) cos θ − sin θ sin θ cos θ (cid:27) . (B.2) Appendix C. Finite-strain solution of cohesive building block
For finite-strain problem, the deformation gradient F (cid:101) and first Piola-Kirchhoff stress P (cid:101) are chosen as thestrain and stress measures, respectively. For 3-dimensional problems, F (cid:101) and P (cid:101) are vectorized to simplifythe computation: F = { F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) , F (cid:101) } T = { F , F , F , ..., F } T , (C.1) P = { P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) , P (cid:101) } T = { P , P , P , ..., P } T . In the generic cohesive building block, the constitutive relation of material 0 in finite-strain formulation canbe written as ∆ F = B ∆ P + δ F , (C.2)where B and δ F are the finite-strain compliance matrices and residual deformation gradient, respectively.The overall constitutive relation of the building block is∆ F = B ∆ P + δ F . (C.3)Based on the equilibrium condition, B and δ F are derived to be B = B + ˜ v R f (˜ α, ˜ β, ˜ γ ) ˜ G f (cid:104) R f (˜ α, ˜ β, ˜ γ ) (cid:105) − (C.4)23nd δ F = δ F + ˜ v R f (˜ α, ˜ β, ˜ γ ) δ ˜ d f . (C.5)The nonzero terms in ˜ G f and δ ˜ d f are˜ G f ([3 , , , [3 , , = G nn G ns G nt G ss G st sym G tt , δ ˜ d f ([3 , , = δd n δd s δd t . (C.6)The compliance matrix G and residual displacement δ d of the cohesive layer are defined in Eq. (2.9). Thedefinition of the finite-strain rotation matrix R f can be found in [41], which can also be decomposed intothree elementary rotation matrices. Appendix D. Tangent stiffness matrix and viscous regularization
Since both DNS and DMN are based on implicit analysis, the tangent stiffness matrix of the cohesivelaw needs to be derived, which is defined as K = ∂ t ∂ d . (D.1)Again, the tensile and compressive cases are treated differently:1. For tensile case d n ≥
0, the components of t are t n = t m d m d n , t s = β t m d m d s , t t = β t m d m d t . (D.2)The independent components in K are K nn = t m d m + d n ∆ , K ns = β d n d s ∆ , K nt = β d n d t ∆ , (D.3) K ss = β t m d m + β d s ∆ , K ss = β t m d m + β d t ∆ , K st = β d s d t ∆with ∆ = (cid:18) t (cid:48) m − t m d m (cid:19) d m . (D.4)2. For compression case d n <
0, the components of t are t n = Kd n , t s = β t m d m d s , t t = β t m d m d t , (D.5)so that K nn = K, K ns = K nt = 0 , (D.6) K ss = β t m d m + β d s ∆ , K ss = β t m d m + β d t ∆ , K st = β d s d t ∆ . The tangent stiffness matrix of the bilinear cohesive law may not be positive definite. At the pointof failure, implicit calculations for both DNS and deep material network are unable to converge to anequilibrium solution, which makes it impossible to track the post-failure behavior. Therefore, an extraviscous term is introduced to overcome the convergence problem following the work by Gao et al.[43]: t ← t + ζ ddt (cid:18) d d f (cid:19) , (D.7)where ζ is a viscosity-like parameter that governs viscous energy dissipation. By selecting a sufficiently smalltime increment in the numerical scheme, one is able to make the tangent stiffness tensor of the softeningmaterial positive definite. 24 eferences [1] A. Melro, P. Camanho, F. A. Pires, S. Pinho, Micromechanical analysis of polymer composites reinforced by unidirectionalfibres: Part ii–micromechanical analyses, International Journal of Solids and Structures 50 (11-12) (2013) 1906–1915.[2] S. Cantournet, R. Desmorat, J. Besson, Mullins effect and cyclic stress softening of filled elastomers by internal slidingand friction thermodynamics model, International Journal of Solids and Structures 46 (11-12) (2009) 2255–2264.[3] F. J. Vernerey, C. McVeigh, W. K. Liu, B. Moran, D. Tewari, D. M. Parks, G. B. Olson, The 3-d computational modelingof shear-dominated ductile failure in steel, Jom 58 (12) (2006) 45–51.[4] C. McVeigh, F. Vernerey, W. K. Liu, B. Moran, G. Olson, An interactive micro-void shear localization mechanism in highstrength steels, Journal of the Mechanics and Physics of Solids 55 (2) (2007) 225–244.[5] A. Simone, C. A. Duarte, E. Van der Giessen, A generalized finite element method for polycrystals with discontinuousgrain boundaries, International Journal for Numerical Methods in Engineering 67 (8) (2006) 1122–1145.[6] H. D. Espinosa, P. D. Zavattieri, A grain level model for the study of failure initiation and evolution in polycrystallinebrittle materials. part i: Theory and numerical implementation, Mechanics of Materials 35 (3-6) (2003) 333–364.[7] A. Needleman, A continuum model for void nucleation by inclusion debonding, Journal of applied mechanics 54 (3) (1987)525–531.[8] H. Duan, J.-x. Wang, Z. Huang, B. L. Karihaloo, Size-dependent effective elastic constants of solids containing nano-inhomogeneities with interface stress, Journal of the Mechanics and Physics of Solids 53 (7) (2005) 1574–1596.[9] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. Lond.A 241 (1226) (1957) 376–396.[10] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Actametallurgica 21 (5) (1973) 571–574.[11] R. Hill, A self-consistent mechanics of composite materials, Journal of the Mechanics and Physics of Solids 13 (4) (1965)213–222.[12] J. Qu, The effect of slightly weakened interfaces on the overall elastic properties of composite materials, Mechanics ofMaterials 14 (4) (1993) 269–281.[13] H. Tan, Y. Huang, C. Liu, P. H. Geubelle, The mori–tanaka method for composite materials with nonlinear interfacedebonding, International Journal of Plasticity 21 (10) (2005) 1890–1918.[14] M. G. Geers, V. G. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges,Journal of computational and applied mathematics 234 (7) (2010) 2175–2182.[15] K. Park, G. H. Paulino, Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces,Applied Mechanics Reviews 64 (6) (2011) 060802.[16] G. T. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, International Journal of solidsand structures 33 (20-22) (1996) 2899–2938.[17] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis,International journal for numerical methods in engineering 44 (9) (1999) 1267–1282.[18] T. Hettich, E. Ramm, Interface material failure modeled by the extended finite-element method and level sets, ComputerMethods in Applied Mechanics and Engineering 195 (37-40) (2006) 4753–4767.[19] T. Belytschko, R. Gracie, G. Ventura, A review of extended/generalized finite element methods for material modeling,Modelling and Simulation in Materials Science and Engineering 17 (4) (2009) 043001.[20] X. Deng, A. Korobenko, J. Yan, Y. Bazilevs, Isogeometric analysis of continuum damage in rotation-free composite shells,Computer Methods in Applied Mechanics and Engineering 284 (2015) 349–372.[21] Y. Bazilevs, M. Pigazzini, A. Ellison, H. Kim, A new multi-layer approach for progressive damage simulation in compositelaminates based on isogeometric analysis and kirchhoff–love shells. part i: basic theory and modeling of delamination andtransverse shear, Computational Mechanics 62 (3) (2018) 563–585.[22] J. Zhao, O. Y. Kontsevoi, W. Xiong, J. Smith, Simulation-aided constitutive law development–assessment of low triaxialityvoid nucleation models via extended finite element method, Journal of the Mechanics and Physics of Solids 102 (2017)30–45.[23] Z. Liu, M. Fleming, W. K. Liu, Microstructural material database for self-consistent clustering analysis of elastoplasticstrain softening materials, Computer Methods in Applied Mechanics and Engineering 330 (2018) 547–577.[24] Z. Liu, C. Wu, M. Koishi, Transfer learning of deep material network for seamless structure–property predictions, Com-putational Mechanics 64 (2) (2019) 451–465.[25] S. Ghosh, Y. Ling, B. Majumdar, R. Kim, Interfacial debonding analysis in multiple fiber reinforced composites, Mechanicsof Materials 32 (10) (2000) 561–591.[26] S. Ghosh, J. Bai, P. Raghavan, Concurrent multi-level model for damage evolution in microstructurally debonding com-posites, Mechanics of Materials 39 (3) (2007) 241–266.[27] C. Oskay, J. Fish, Eigendeformation-based reduced order homogenization for failure analysis of heterogeneous materials,Computer Methods in Applied Mechanics and Engineering 196 (7) (2007) 1216 – 1243.[28] Z. Yuan, J. Fish, Multiple scale eigendeformation-based reduced order homogenization, Computer Methods in AppliedMechanics and Engineering 198 (21-26) (2009) 2016–2038.[29] X. Zhang, C. Oskay, Eigenstrain based reduced order homogenization for polycrystalline materials, Computer Methods inApplied Mechanics and Engineering 297 (2015) 408–436.[30] Z. Liu, M. Bessa, W. K. Liu, Self-consistent clustering analysis: An efficient multi-scale scheme for inelastic heterogeneousmaterials, Computer Methods in Applied Mechanics and Engineering 306 (2016) 319–341.
31] Z. Liu, O. L. Kafka, C. Yu, W. K. Liu, Data-driven self-consistent clustering analysis of heterogeneous materials withcrystal plasticity, in: Advances in Computational Plasticity, Springer, 2018, pp. 221–242.[32] M. Shakoor, J. Gao, Z. Liu, W. K. Liu, A data-driven multiscale theory for modeling damage and fracture of compositematerials, in: International Workshop on Meshfree Methods for Partial Differential Equations, Springer, 2017, pp. 135–148.[33] H. Li, O. L. Kafka, J. Gao, C. Yu, Y. Nie, L. Zhang, M. Tajdari, S. Tang, X. Guo, G. Li, et al., Clustering discretizationmethods for generation of material performance databases in machine learning and design optimization, ComputationalMechanics 64 (2) (2019) 281–305.[34] J. Oliver, M. Caicedo, A. Huespe, J. Hern´andez, E. Roubin, Reduced order modeling strategies for computational multiscalefracture, Computer Methods in Applied Mechanics and Engineering 313 (2017) 560–595.[35] J. Ghaboussi, J. Garrett Jr, X. Wu, Knowledge-based modeling of material behavior with neural networks, Journal ofengineering mechanics 117 (1) (1991) 132–153.[36] B. Le, J. Yvonnet, Q.-C. He, Computational homogenization of nonlinear elastic materials using neural networks, Inter-national Journal for Numerical Methods in Engineering 104 (12) (2015) 1061–1084.[37] M. Bessa, R. Bostanabad, Z. Liu, A. Hu, D. Apley, C. Brinson, W. Chen, W. Liu, A framework for data-driven analysisof materials under uncertainty: Countering the curse of dimensionality, Computer Methods in Applied Mechanics andEngineering 320 (2017) 633–667.[38] K. Wang, W. Sun, A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deeplearning, Computer Methods in Applied Mechanics and Engineering 334 (2018) 337–380.[39] K. Wang, W. Sun, Meta-modeling game for deriving theory-consistent, microstructure-based traction–separation laws viadeep reinforcement learning, Computer Methods in Applied Mechanics and Engineering 346 (2019) 216–241.[40] Z. Liu, C. T. Wu, M. Koishi, A deep material network for multiscale topology learning and accelerated nonlinear modelingof heterogeneous materials, Computer Methods in Applied Mechanics and Engineering 345 (2019) 1138–1168.[41] Z. Liu, C. Wu, Exploring the 3d architectures of deep material network in data-driven multiscale mechanics, Journal ofthe Mechanics and Physics of Solids 127 (2019) 20 – 46.[42] R. Banerjee, K. Sagiyama, G. Teichert, K. Garikipati, A graph theoretic framework for representation, exploration andanalysis on computed states of physical systems, Computer Methods in Applied Mechanics and Engineering 351 (2019)501–530.[43] Y. Gao, A. Bower, A simple technique for avoiding convergence problems in finite element simulations of crack nucleationand growth on cohesive interfaces, Modelling and Simulation in Materials Science and Engineering 12 (3) (2004) 453.31] Z. Liu, O. L. Kafka, C. Yu, W. K. Liu, Data-driven self-consistent clustering analysis of heterogeneous materials withcrystal plasticity, in: Advances in Computational Plasticity, Springer, 2018, pp. 221–242.[32] M. Shakoor, J. Gao, Z. Liu, W. K. Liu, A data-driven multiscale theory for modeling damage and fracture of compositematerials, in: International Workshop on Meshfree Methods for Partial Differential Equations, Springer, 2017, pp. 135–148.[33] H. Li, O. L. Kafka, J. Gao, C. Yu, Y. Nie, L. Zhang, M. Tajdari, S. Tang, X. Guo, G. Li, et al., Clustering discretizationmethods for generation of material performance databases in machine learning and design optimization, ComputationalMechanics 64 (2) (2019) 281–305.[34] J. Oliver, M. Caicedo, A. Huespe, J. Hern´andez, E. Roubin, Reduced order modeling strategies for computational multiscalefracture, Computer Methods in Applied Mechanics and Engineering 313 (2017) 560–595.[35] J. Ghaboussi, J. Garrett Jr, X. Wu, Knowledge-based modeling of material behavior with neural networks, Journal ofengineering mechanics 117 (1) (1991) 132–153.[36] B. Le, J. Yvonnet, Q.-C. He, Computational homogenization of nonlinear elastic materials using neural networks, Inter-national Journal for Numerical Methods in Engineering 104 (12) (2015) 1061–1084.[37] M. Bessa, R. Bostanabad, Z. Liu, A. Hu, D. Apley, C. Brinson, W. Chen, W. Liu, A framework for data-driven analysisof materials under uncertainty: Countering the curse of dimensionality, Computer Methods in Applied Mechanics andEngineering 320 (2017) 633–667.[38] K. Wang, W. Sun, A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deeplearning, Computer Methods in Applied Mechanics and Engineering 334 (2018) 337–380.[39] K. Wang, W. Sun, Meta-modeling game for deriving theory-consistent, microstructure-based traction–separation laws viadeep reinforcement learning, Computer Methods in Applied Mechanics and Engineering 346 (2019) 216–241.[40] Z. Liu, C. T. Wu, M. Koishi, A deep material network for multiscale topology learning and accelerated nonlinear modelingof heterogeneous materials, Computer Methods in Applied Mechanics and Engineering 345 (2019) 1138–1168.[41] Z. Liu, C. Wu, Exploring the 3d architectures of deep material network in data-driven multiscale mechanics, Journal ofthe Mechanics and Physics of Solids 127 (2019) 20 – 46.[42] R. Banerjee, K. Sagiyama, G. Teichert, K. Garikipati, A graph theoretic framework for representation, exploration andanalysis on computed states of physical systems, Computer Methods in Applied Mechanics and Engineering 351 (2019)501–530.[43] Y. Gao, A. Bower, A simple technique for avoiding convergence problems in finite element simulations of crack nucleationand growth on cohesive interfaces, Modelling and Simulation in Materials Science and Engineering 12 (3) (2004) 453.