Deep material network with cohesive layers: Multi-stage training and interfacial failure analysis

TL;DR

This paper introduces an enhanced deep material network with cohesive layers for accurate, efficient multiscale modeling of heterogeneous materials, especially capturing interfacial failure, using a multi-stage training approach based on linear elastic simulations.

Contribution

The paper presents a novel data-driven model that incorporates cohesive layers into the deep material network, enabling better prediction of interfacial failure with reduced computational effort.

Findings

Accurate prediction of interfacial failure in composites.

Significant reduction in computational time.

Effective multi-stage training strategy using linear elastic data.

Abstract

A fundamental issue in multiscale materials modeling and design is the consideration of traction-separation behavior at the interface. By enriching the deep material network (DMN) with cohesive layers, the paper presents a novel data-driven material model which enables accurate and efficient prediction of multiscale responses 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 cohesive layer. It is shown that the enriched material network can be effectively optimized via a multi-stage training strategy, with training data generated only from linear elastic direct numerical simulation (DNS). The extrapolation capability of the method to unknown material and loading spaces is demonstrated through the debonding analysis of a unidirectional fiber-reinforced composite, where the interface behavior is governed by an irreversible softening mixed-mode cohesive law. Its predictive accuracy is validated against the nonlinear path-dependent DNS results, and the reduction in computational time is particularly significant.

Tables

Table 1. Table 1: Constitutive relations in the cohesive building block. The final expressions of the compliance matrix D and the residual strain δ ​ 𝜺 𝛿 𝜺 \delta\boldsymbol{\varepsilon} of are given in Eq. ( 2.27 ) and ( 2.28 )， respectively.

Inputs	Cohesive layer	$\Delta \text{d}=\text{G}\Delta \text{t}+\delta \text{d}$
	Material 0	$\Delta {𝜺}^0={\text{D}}^0\Delta {𝝈}^0+\delta {𝜺}^0$
Intermediate steps	(i) Rotation of material 0	$\Delta {𝜺}^r={\text{D}}^r\Delta {𝝈}^r+\delta {𝜺}^r$
	(ii) Homogenization	$\Delta \overline{𝜺}=\overline{\text{D}}\Delta \overline{𝝈}+\delta \overline{𝜺}$
Output step	(iii) Rotation of homogenized material	$\Delta 𝜺=\text{D}\Delta 𝝈+\delta 𝜺$

Table 2. Table 2: Material parameters of the fiber phase, the matrix phase and the interfaces used in the online stage.

Fiber	$E_1$ (GPa)	${\nu }_1$
	500	0.3
Matrix	$E_2$ (GPa)	${\nu }_2$	${\sigma }_2^y$ (GPa)	Hardening
	100	0.3	0.1	Eq. (5.2)
Interface	$K$ (GPa/m)	${\sigma }_c$ (GPa)	$G_c$ (GPa$\cdot$m)	$\beta$	$\zeta$ (GPa$\cdot$s)
	$1\times {10}^7$	0.5	$2.5\times {10}^{-6}$	0.5	$2\times {10}^{-5}$

Table 3. Table 3: Training results of UD RVE at stages I and II. Average training error e ¯ t ​ r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , average test error e ¯ t ​ r superscript ¯ 𝑒 𝑡 𝑟 \bar{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 subscript 𝑁 𝑐 4 N_{c}=4 .

		Stage I			Stage II
	Training ${\overline{e}}^{tr}$	Test ${\overline{e}}^{te}$	Maximum $e_s^{te}$	Training ${\overline{e}}^{tr}$	Test ${\overline{e}}^{te}$	Maximum $e_s^{te}$
$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%

Table 4. Table 4: Training results of the particle-reinforced RVE at stages I and II. Average training error e ¯ t ​ r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , average test error e ¯ t ​ r superscript ¯ 𝑒 𝑡 𝑟 \bar{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 subscript 𝑁 𝑐 4 N_{c}=4 .

		Stage I			Stage II
	Training ${\overline{e}}^{tr}$	Test ${\overline{e}}^{te}$	Maximum $e_s^{te}$	Training ${\overline{e}}^{tr}$	Test ${\overline{e}}^{te}$	Maximum $e_s^{te}$
$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%

Equations

$$z^j,{\alpha }_i^k,{\beta }_i^k,{\gamma }_i^k\text{with }i=1,2,...,N;j=1,2,...,2^{N-1};k=1,2,...,2^{i-1}.$$

$$w_i^k=\sum _{j=2^{N-i}(k-1)+1}^{2^{N-i}k}{w}_N^j=\sum _{j=2^{N-i}(k-1)+1}^{2^{N-i}k}a(z^j).$$

$$j=\mathcal{J}(q)\text{with }q=1,2,...,N_e.$$

$${\tilde{z}}_q^p,{\tilde{\alpha }}_q^p,{\tilde{\beta }}_q^p,{\tilde{\gamma }}_q^p\text{with }p=1,2,...,N_c.$$

$$N^{total}=\stackrel{\text{Material network}}{\overbrace{\left(7\times 2^{N-1}-3\right)}}+\stackrel{\text{Cohesive networks}}{\overbrace{N_e\times N_c}}.$$

$$\tilde{v}=\frac{1}{T_3}.$$

$${\tilde{v}}_q^p=\frac{a({\tilde{z}}_q^p)}{L}=\frac{\max ({\tilde{z}}_q^p,0)}{L},$$

$${\tilde{S}}_q^p={\tilde{v}}_q^p{\Omega }_q.$$

$$\Delta \text{d}=\text{G}\Delta \text{t}+\delta \text{d},$$

$$\text{G}=\left\{\begin{array}{ccc}{G}_{nn}& G_{ns}& G_{nt}\\ & G_{ss}& G_{st}\\ sym& & G_{tt}\end{array}\right\}.$$

$$\text{K}={\text{G}}^{-1}=\left\{\begin{array}{ccc}{K}_{nn}& K_{ns}& K_{nt}\\ & K_{ss}& K_{st}\\ sym& & K_{tt}\end{array}\right\}.$$

$$\mathbf{\sigma }=\{{\underset{˜}{\sigma }}_{11},{\underset{˜}{\sigma }}_{22},{\underset{˜}{\sigma }}_{33},\sqrt{2}{\underset{˜}{\sigma }}_{23},\sqrt{2}{\underset{˜}{\sigma }}_{13},\sqrt{2}{\underset{˜}{\sigma }}_{12}{\}}^T=\{{\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4,{\sigma }_5,{\sigma }_6{\}}^T$$

$$\mathbf{\epsilon }=\{{\underset{˜}{\epsilon }}_{11},{\underset{˜}{\epsilon }}_{22},{\underset{˜}{\epsilon }}_{33},\sqrt{2}{\underset{˜}{\epsilon }}_{23},\sqrt{2}{\underset{˜}{\epsilon }}_{13},\sqrt{2}{\underset{˜}{\epsilon }}_{12}{\}}^T=\{{\epsilon }_1,{\epsilon }_2,{\epsilon }_3,{\epsilon }_4,{\epsilon }_5,{\epsilon }_6{\}}^T,$$

$$\Delta {\mathbf{\epsilon }}^0={\text{D}}^0\Delta {\mathbf{\sigma }}^0+\delta {\mathbf{\epsilon }}^0,$$

$$\Delta {\mathbf{\epsilon }}^r={\text{D}}^r\Delta {\mathbf{\sigma }}^r+\delta {\mathbf{\epsilon }}^r.$$

$$\Delta \bar{\mathbf{\epsilon }}=\bar{\text{D}}\Delta \bar{\mathbf{\sigma }}+\delta \bar{\mathbf{\epsilon }}.$$

$$\Delta \mathbf{\epsilon }=\text{D}\Delta \mathbf{\sigma }+\delta \mathbf{\epsilon }.$$

$$\text{D}=\text{D}\left({\text{D}}^0,\text{G},\tilde{z},\tilde{\alpha },\tilde{\beta },\tilde{\gamma }\right)$$

$$\delta \mathbf{\epsilon }=\delta \mathbf{\epsilon }\left(\delta {\mathbf{\epsilon }}^0,\delta \text{d},\tilde{z},\tilde{\alpha },\tilde{\beta },\tilde{\gamma }\right).$$

$$\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })=\text{X}(\tilde{\alpha })\text{Y}(\tilde{\beta })\text{Z}(\tilde{\gamma }).$$

$${\text{D}}^r={\text{R}}^{-1}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma }){\text{D}}^0\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma }),\delta {\mathbf{\epsilon }}^r={\text{R}}^{-1}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })\delta {\mathbf{\epsilon }}^0.$$

$$\Delta t_n=\Delta {\bar{\sigma }}_3,\sqrt{2}\Delta t_s=\Delta {\bar{\sigma }}_4,\sqrt{2}\Delta t_t=\Delta {\bar{\sigma }}_5.$$

$$\Delta {\mathbf{\sigma }}^r=\Delta \bar{\mathbf{\sigma }},\Delta {\epsilon }_1^r=\Delta {\bar{\epsilon }}_1,\Delta {\epsilon }_2^r=\Delta {\bar{\epsilon }}_2,\Delta {\epsilon }_6^r=\Delta {\bar{\epsilon }}_6.$$

$$\Delta {\bar{\epsilon }}_3=\Delta {\epsilon }_3^r+\tilde{v}\Delta d_n,\Delta {\bar{\epsilon }}_4=\Delta {\epsilon }_4^r+\frac{\tilde{v}\Delta d_s}{\sqrt{2}},\Delta {\bar{\epsilon }}_5=\Delta {\epsilon }_5^r+\frac{\tilde{v}\Delta d_t}{\sqrt{2}}.$$

$$\bar{\text{D}}={\text{D}}^r+\tilde{v}\tilde{\text{G}},\delta \bar{\mathbf{\epsilon }}=\delta {\mathbf{\epsilon }}^r+\tilde{v}\delta \tilde{\text{d}},$$

$$\tilde{\text{G}}=\left\{\begin{array}{cccccc}0& & & & & \\ & 0& & & & \\ & & G_{nn}& G_{ns}/\sqrt{2}& G_{nt}/\sqrt{2}& \\ & & & G_{ss}/2& G_{st}/2& \\ & & sym& & G_{tt}/2& \\ & & & & & 0\end{array}\right\},\delta \tilde{\text{d}}=\left\{\begin{array}{c}0\\ 0\\ \delta d_n\\ \delta d_s/\sqrt{2}\\ \delta d_t/\sqrt{2}\\ 0\end{array}\right\}.$$

$$\text{D}=\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })\bar{\text{D}}{\text{R}}^{-1}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma }),\delta \mathbf{\epsilon }=\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })\delta \bar{\mathbf{\epsilon }}.$$

$$\text{D}=\text{D}\left({\text{D}}^0,\text{G},\tilde{v},\tilde{\alpha },\tilde{\beta },\tilde{\gamma }\right)={\text{D}}^0+\tilde{v}\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })\tilde{\text{G}}{\text{R}}^{-1}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })$$

$$\delta \mathbf{\epsilon }=\delta \mathbf{\epsilon }\left(\delta {\mathbf{\epsilon }}^0,\delta \text{d},\tilde{v},\tilde{\alpha },\tilde{\beta },\tilde{\gamma }\right)=\delta {\mathbf{\epsilon }}^0+\tilde{v}\text{R}(\tilde{\alpha },\tilde{\beta },\tilde{\gamma })\delta \tilde{\text{d}}.$$

$${\text{D}}_q^{N_c}={\text{D}}_q^0+\sum _{p=1}^{N_c}{\tilde{v}}_q^p\text{R}({\tilde{\alpha }}_q^p,{\tilde{\beta }}_q^p,{\tilde{\gamma }}_q^p){\tilde{\text{G}}}_q^p{\text{R}}^{-1}({\tilde{\alpha }}_q^p,{\tilde{\beta }}_q^p,{\tilde{\gamma }}_q^p),$$

Full text

1 Introduction

Materials have deformable interfaces across fine-scale phases, and an emerging need in multiscale materials modeling and design is to capture the effects of interfaces at various length and time scales. Some representative 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 ductile fracture of metals [3, 4], and grain-boundary effects in polycrystalline materials [5, 6]. One can see that the interface, in many cases, is the critical region involving damage initiation and evolution, thus, it plays an important role in the failure analysis of modern material systems. On the other hand, as the interface is a geometric object with one dimension less than the bulk material, it also results in various interesting phenomena 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. Analytical micromechanics methods based on the Eshelby’s solution [9] and mean-field theories [10, 11] can be enhanced to incorporate the interfacial effect, but they are usually limited to regular geometries and weakly nonlinear behaviors [12, 13, 8]. To resolve these issues, computational homogenization is emerging as an accurate way of deriving multiscale material responses from the representative volume element (RVE) [14], where the microstructures are explicitly constructed and discretized by numerical models.

In recent years, direct numerical simulation (DNS) of RVE has been widely used to understand the role of interface in heterogeneous materials. To describe the traction-separation relationships across the interfaces, 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 (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 inclusions to investigate the cause of void nucleation due to the inclusion debonding. Melro et al. [1] developed RVE models for polymer composite reinforced by unidirectional fibers and evaluated the influence of fiber-matrix interface on the strength properties of the composite. Zhao et al. [22] utilized the XFEM to simulate the multi-void nucleation process due to interface debonding under various loading conditions and particle size/shape distributions. However, limited by the large computational cost of DNS, these full-field RVE models 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 have been 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], and proper 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 responses of elastic and history-dependent materials [35, 36, 37, 38, 33]. Regarding interface modeling, Wang et al.[39] recently proposed a deep reinforcement learning framework based on directed graph for deriving the traction-separation law. As the training data is usually limited due to the expense of physical or numerical experiments, a key challenge faced by most data-driven methods is the credibility of extrapolation to 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 of the 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 building blocks. Essentially, all fitting parameters of DMN are interpretable with physical meanings related to the RVE morphology. Trained only by linear elastic data, DMN can be extrapolated to capture highly nonlinear material behaviors for various applications, including hyperelastic rubber composite under large deformation, polycrystalline materials with rate-dependent crystal plasticity, and various types of CFRPs. Recently, a transfer-learning framework of DMN is further introduced to generate a unified set of material databases covering the full-range structure-property relationship [24].

This paper extends the deep material network framework [41] to 3-dimensional heterogeneous materials with 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 describe the size effect of interface;

• •

Multi-stage training strategy to reduce the number of fitting parameters in the machine-learning process at each stage;

• •

Extrapolation of the enriched DMN from linear elasticity to interfacial failure analysis with irreversible softening cohesive laws.

The remainder of this paper is organized as follows. In Section 2, the new design of DMN architecture with cohesive layers is discussed in detail. Analytical solutions of the cohesive building block are derived to provide the basis of data propagation in the network. Section 3 presents the multi-stage training strategy based on linear elastic offline DNS data. Section 4 briefly introduces the mix-mode cohesive law and necessary viscous regularization for implicit analysis in the online stage. Section 5 shows the application of DMN to the interfacial debonding analysis of a unidirectional fiber-reinforced composite, and another example on a particle-reinforced composite is also provided in A. Concluding remarks and future work are given in Section 6.

2 New DMN architecture with enrichment

As a starting point, consider a material network based on a binary tree 111Current DMN architectures are based on the tree graph, in which the physical information of any node (e.g. the stiffness tensor 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. with a two-layer structure as the building block [40, 41], as shown on the right of Figure 1. The microscale constituents are input to the bottom layer, and the mechanical information such as the stiffness tensor and the residual stress is feed forward through the input-output functions defined by the physics-based building block. The macroscale material 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-1}$ 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-layer structures with unknown phase fractions and block orientations [41]. The fitting parameters of a material network are the activation $z$ of each bottom-layer node, and three rotation angles $(\alpha,\beta,\gamma)$ at each building block:

[TABLE]

The activation $z^{j}$ determines the weight of each node in the material network. For the $k$-th node in the $i$-th layer, its weight $w^{k}_{i}$ can be expressed as a summation of activations of its descendant nodes in the bottom layer:

[TABLE]

where $a(\cdot)$ is the the rectified linear unit (ReLU). Each building block takes the materials at two child nodes as inputs, and their weights determine the phase fraction of the two-layer structure. A homogenization operation is first performed to mix the two materials. The new homogenized material then undergoes a rotation defined by the angles $(\alpha,\beta,\gamma)$, and the rotated material is output to the mother node. All fitting parameters of DMN have physical meanings related to the microstructural topology, and they can be optimized through gradient-based methods by virtue of the existence of analytical solutions.

2.1 Design of cohesive networks

The original architecture of material network is created for heterogeneous materials with perfectly bonded phases. Its extension to a material with interfacial effect (e.g. debonding) is limited by the design of the building block, where the two input materials are perfectly bonded. In this regard, it is possible to consider the interface in the original two-layer building block to be deformable. However, this naive approach may suffer from two issues: 1) the interface orientation is predetermined by the building block, thus, loses the adaptivity; 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 the microscale constituents.

An alternative viewpoint regards the nodes at the bottom layer of the original network as degrees of freedom (DOF). In the online computation, each active bottom-layer node serves as an independent material point, which carries the pure information of a certain microscale constituent. This provides the physical basis to introduce the interfaces only at the bottom layer. In the paper, it is proposed to enrich each selected bottom-layer node with a so-called “cohesive network”, as shown in Figure 1. The index of the enriched node and its corresponding cohesive network is $q$, and the total number of enriched nodes is denoted by $N_{e}$. A mapping function $\mathcal{J}(q)$ finds the index of enriched node $q$ in the bottom layer of the original material network:

[TABLE]

For example, if the second cohesive network ($q=2$) enriches the third node in the bottom layer, we have $\mathcal{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 one single phase. For example, in the material network for unidirectional (UD) fiber-reinforced composite with deformable fiber-matrix interface to be tested in Section 5, only the active nodes of the fiber phase are linked to 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 microscale constituent 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 “cohesive building block” (see Section 2.3). For simplicity, all the cohesive networks are assumed to have the same number 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

[TABLE]

In this paper, the tilde symbol is used to indicate that the corresponding quantity is defined within the cohesive network. The activation of the $p$-th cohesive layer is $\tilde{z}_{q}^{p}$, and rotation angles of the $p$-th cohesive network are $\tilde{\alpha}_{q}^{p},\tilde{\beta}_{q}^{p},\tilde{\gamma}_{q}^{p}$. Physical meanings of all these fitting parameters will be explained in the following two subsections. Evidently, an important condition that the cohesive network should satisfy is that if all the cohesive layers are deactivated, the base material should be returned to the enriched node without any change.

The total number of parameters in the new material network with enriching cohesive networks is

[TABLE]

In the multi-stage training strategy to be discussed in Section 3, not all parameters are fitted at the same time. The basic idea is to train the original material network first and then the fitted parameters ($z,\alpha,\beta,\gamma$) are treated as constants in the subsequent training of the integrated network with cohesive layers. To ease the offline sampling and training processes, the fitting parameters are optimized based on linear elastic DNS data, which is sufficient for capturing essential physical interactions among the microscale material phases and interfaces. The learned model will be extrapolated to unknown material and loading spaces, such as plasticity, interfacial failure (debonding), and loading-unloading paths.

2.2 Activation of cohesive layers

In a generic cohesive building block shown in Figure 2, a non-negative reciprocal length parameter $\tilde{v}$ is introduced for the cohesive layer to represent its size effect. Physically, one can interpret $\tilde{v}$ as the inverse of effective thickness $T_{3}$ of the undeformed building block in the direction normal to the cohesive layer:

[TABLE]

Since the deformation of the cohesive layer is characterized by its separation displacement, this reciprocal length parameter $\tilde{v}$ will be further used to transform the displacement vector to the strain tensor contributed to the building block. It is noted that when $\tilde{v}=0$, the building block has infinite effective thickness and sees 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 $\tilde{v}_{q}^{p}$ is designed to be activated by the rectified linear unit (ReLU) $a(\cdot)$:

[TABLE]

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 help to keep the activations vary in a similar range as the rotation angles, and thus it eases the convergence of the gradient descent algorithm used for the training. The choice of characteristic length $L$ depends on the morphology of the microstructure. For example, in a unidirectional fiber-reinforced composite, $L$ can be chosen as the average diameter of the fibers. While for another polycrystalline material, $L$ can be chosen as the average grain size.

The activation of cohesive layer $\tilde{z}_{q}^{p}$ in Eq. (2.7) is a unitless fitting parameter to be learned. Similar to the activation of node $z$ in the original material network, once $\tilde{z}_{q}^{p}$ becomes negative during the training, the corresponding cohesive layer is deactivated, and it will not be activated again due to the vanished gradient of the ReLU function.

Assume the volume of the $q$-th enriched node is given by $\Omega_{q}$. Since the cohesive layers have zero volume, all the cohesive building blocks in the $q$-th cohesive network share the same volume $\Omega_{q}$. Following the definition of $\tilde{v}$ in Eq. (2.6, one can define the effective area of the $p$-th cohesive layer $\tilde{S}_{q}^{p}$ as

[TABLE]

Note that the RVE scaling effect due to the existence of deformable interfaces can be considered by scaling the characteristic length $L$ accordingly. For example, assume an RVE is trained to a DMN with a set of activations $\tilde{z}$ under $L$. If the RVE geometry is enlarged linearly by a scale factor of 2 and the cohesive law stays the same, the new DMN can be obtained by simply updating the characteristic length to $2L$ without changing the trained fitting parameters. According to Eq. (2.7), this is also equivalent to decreasing all the reciprocal length parameters $\tilde{v}$ by a factor of 0.5.

2.3 Cohesive building block

As depicted in Figure 2, a generic cohesive building block has two inputs: a planar cohesive layer and a bulk material 0. The cohesive layer has zero thickness, and its normal is in direction 3. In general, its constitutive relation for an arbitrary traction-separation law can be linearized as

[TABLE]

where $\Delta\textbf{d}$ is the increment of displacement vector, and $\Delta\textbf{t}$ is the increment of traction vector. The residual displacement vector $\delta\textbf{d}$ will appear in a nonlinear traction-separation law, including the bilinear one used in this work. For a 3-dimensional problem, the traction t is in the unit of force per unit area. It has three components: $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 independent components in 3D due to symmetry:

[TABLE]

The stiffness matrix K is defined as the inverse of G:

[TABLE]

Note that for history-dependent traction-separation laws, internal variables will also be carried at the cohesive layer, while the incremental constitutive relation 2.9 still applies at each loading step. For layers with perfectly plasticity, the stiffness in certain direction may become zero, and the corresponding elements in the compliance matrix G will be infinite. To avoid the potential numerical issues, a tiny stiffness number is recommended to be added to the diagonal terms in K. The same practice is also adopted in the paper for interfacial failure analysis (see Eq. (4.13), when the cohesive layer is fully separated or damaged.

In a small-strain setting, the Cauchy stress $\boldsymbol{\sigma}$ and the infinitesimal strain $\boldsymbol{\varepsilon}$ are used as stress and strain measures, respectively. Moreover, the stress and strain are vectorized under the Mandel notation:

[TABLE]

and

[TABLE]

where the subscripts 3, 4, and 5 denote the shear directions. Equipped with these conventions, the constitutive relation of Material 0 can be written as

[TABLE]

where $\Delta\boldsymbol{\varepsilon}^{0}$ and $\Delta\boldsymbol{\sigma}^{0}$ are the strain and stress increments, respectively. $\textbf{D}^{0}$ is its compliance matrix, which is symmetric and has 21 independent entries for a 3-dimensional problem. Additionally, $\delta\boldsymbol{\varepsilon}^{0}$ is the residual strain.

Three operations are introduced in the building block to propagate the mechanical information including the compliance matrix and residual strain/displacement. First, material 0 is rotated and the obtained constitutive relation is

[TABLE]

Then rotated material 0 is homogenized together with the cohesive layer, and the new constitutive relation is

[TABLE]

To recover the initial orientation of material 0, the homogenized material is rotated back, which gives the overall constitutive relation of the cohesive building block:

[TABLE]

The associated constitutive relations of the inputs, intermediate steps (i, ii), and output step (iii) are summarized in Table 1.

In the generic building block, the overall input-output function defined on the compliance matrix can be generally written as

[TABLE]

and the one for the residual strain/displacement is

[TABLE]

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 ($\tilde{\alpha},\tilde{\beta},\tilde{\gamma}$). The rotation matrix R can be decomposed as a product of three element rotation matrices,

[TABLE]

The expressions of these elementary rotation matrices can be found in B. The compliance matrix $\textbf{D}^{r}$ and residual strain $\delta\boldsymbol{\varepsilon}^{r}$ of the rotated material 0 are

[TABLE]

(ii) Homogenization: The equilibrium condition at the interface poses that the tractions on the cohesive layer should satisfy

[TABLE]

For the rotated material 0, the cohesive layer does not apply any displacement or stress constraints on its 1, 2, and 6 directions. Along with the equilibrium condition, we have

[TABLE]

The contribution of displacement vector at the cohesive layer $\Delta\textbf{d}$ to the homogenized strain $\Delta\bar{\boldsymbol{\varepsilon}}$ is scaled by the reciprocal length parameter $\tilde{v}$, which is interpreted as the inverse of effective thickness of the building block in 3 direction. The unknown components of the homogenized strain can be computed as

[TABLE]

Therefore, the compliance matrix and residual strain after the homogenization step are derived to be

[TABLE]

where $\tilde{\textbf{G}}$ and $\delta\tilde{\textbf{d}}$ are the modified compliance matrix and displacement vector of the cohesive layer:

[TABLE]

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 cohesive layer are acting on the stress components.

(iii) Rotation of homogenized material: The orientation of the material 0 is recovered by an inverse rotation $\textbf{R}^{-1}(\tilde{\alpha},\tilde{\beta},\tilde{\gamma})$ on the homogenized material, and the overall compliance matrix and residual strain of the cohesive building block are

[TABLE]

By combining Eq. 2.20, 2.24, and 2.26, the input-output functions of the cohesive building block are derived to be

[TABLE]

and

[TABLE]

Apparently, if $\tilde{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 are given in C. However, all examples presented in this paper are based on small-strain formulation.

2.4 Data transfer between cohesive and material networks

By applying the input-output functions of the cohesive building block iteratively as illustrated in Figure 1, the compliance matrix at the top of the cohesive network for the $q$-th enriched node with $N_{c}$ cohesive layers is obtained as

[TABLE]

and the residual strain is

[TABLE]

where $\textbf{D}_{q}^{0}$ and $\delta\boldsymbol{\varepsilon}_{q}^{0}$ are the compliance matrix and residual strain of the base material of the $q$-th cohesive network, respectively. Meanwhile, $\tilde{\textbf{G}}_{q}^{p}$ and $\delta\tilde{\textbf{d}}_{q}^{p}$ are modified from the compliance matrix and displacement vector 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 enriched node in the bottom layer ($i=N$) of the original material network. As the stiffness tensor C and residual stress $\delta\boldsymbol{\sigma}$ are primarily used in the material network, the following operations are performed:

[TABLE]

The definition of the mapping function $\mathcal{J}(q)$ is provided in Eq. 2.3.

3 Multi-stage offline training strategy

A two-stage training strategy is proposed in this work. The first stage learns the essential phase topology with no interfacial effect, while the second stage learns the length scale and orientation of the cohesive layers. In comparison to learning all the fitting parameters in one single step, the two-stage straining strategy decouples the interfacial effect from phase interactions, so that it has fewer parameters at each stage. Overall, it makes the landscape of the cost function less complicated for the learning algorithms to converge.

Figure 3 illustrated this strategy for a shallow material network with cohesive layers. Its number of layers $N$, number of enriched nodes $N_{e}$, number of cohesive layers in each cohesive network $N_{c}$, and mapping function $\mathcal{J}(q)$ between indices of the cohesive network and bottom-layer node in the original material network are listed below:

[TABLE]

Since the offline training processes are based on linear elastic data, the material network only propagates the information of compliance and stiffness matrices.

As shown in the figure, the material network without the enrichment of cohesive networks is first trained for the RVE with perfectly bonded interfaces. The overall stiffness matrix of the multi-layer network $\bar{\textbf{C}}^{rve}$ for a two-phase RVE can be written as a function of the stiffness matrix from each material phase ($\textbf{C}^{p1}$ and $\textbf{C}^{p2}$) and the fitting parameters ($z,\alpha,\beta,\gamma$):

[TABLE]

At training stage II, each active bottom-layer node of phase 1 is linked to a cohesive network to form a new network for the RVE with deformable interfaces. Therefore, the stiffness matrix of the cohesive layers $\textbf{K}^{c}$ enters as an extra input. Note that the fitting parameters of the two-layer building blocks obtained in the first training stage will be kept constant. With $(\tilde{z},\tilde{\alpha},\tilde{\beta},\tilde{\gamma})$ as the new fitting parameters, the overall stiffness matrix the multi-layer network $\tilde{\textbf{C}}^{rve}$ can be written as

[TABLE]

For both stages, the cost function $J$ is defined upon the mean square error (MSE):

[TABLE]

where $\bar{\textbf{C}}^{dns}_{s}$ is the overall stiffness matrix of sample $s$ computed from DNS. $N_{s}$ is the total number of training samples. The operator $||...||$ denotes the Frobenius matrix norm. While for training stage I, an extra regularization term is added to the cost function to constrain the magnitude of $\{z^{j}\}$ and makes the optimization problem well-posed [41].

Stochastic gradient descent (SGD) is used for minimizing the cost function, while the gradients are calculated 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.

3.1 Training stage I

The inputs $\textbf{C}^{p1}$ and $\textbf{C}^{p2}$ are assigned to the nodes in the bottom layer $N$ in a scheme described as below,

[TABLE]

This scheme is applicable to two-phase materials. The material network based on binary-tree structure can also handle materials with more than two phases, however, a different assigning scheme needs to be developed.

In the sampling of inputs, both material phases are assumed to be orthotropically elastic, whose compliance matrices can be written as

[TABLE]

where $pi$ is $p1$ or $p2$. To initiate material anisotropy, the tension moduli of each phase are first randomly sampled:

[TABLE]

where $U$ represents the uniform distribution. The moduli of phase 2 will be rescaled to create phase contrasts, and the rescaling factor $\bar{E}^{p2}$ is sampled from

[TABLE]

With the rescaling factor, the tension moduli of each phase are updated,

[TABLE]

After all the tension moduli are determined, the shear moduli are sampled from

[TABLE]

To guarantee the compliance matrices are positive definite, the Poisson’s ratios are sampled from

[TABLE]

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,\alpha,\beta,\gamma$) are initialized randomly following uniform distributions. 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].

3.2 Training stage II

In training stage II, the active nodes of phase 1 in the bottom layer are enriched by the cohesive network, while the input $\textbf{C}^{p2}$ is still assigned to the other active nodes of phase 2:

[TABLE]

Meanwhile, the input $\textbf{C}^{p1}$ is assigned to the base material of each cohesive network, and all the cohesive layers receive the same stiffness matrix $\textbf{K}^{c}$:

[TABLE]

The samples generated from training stage I are reused and augmented with new elastic cohesive properties. The cohesive layer is assumed to be isotropic in the two shear directions. In addition, the off-diagonal terms in $\textbf{K}^{c}$ are set to zero. As a result, $\textbf{K}^{c}$ has only two independent components $K_{nn}^{c}$ and $K_{ss}^{c}$:

[TABLE]

The elastic constant in the normal direction $K_{nn}^{c}$ is sampled as

[TABLE]

where $L$ is the characteristic length of the microstructure, and a discussion about the choice of $L$ can be found in Section 2.2. For the elastic constant in the shear direction $K_{ss}^{c}$, we have

[TABLE]

Although the new input space has 15 independent variables, only two of them need to be sampled at this stage and are further appended to the previous space. Similarly, Monte Carlo method is used for the sampling.

The fitting parameters $(\tilde{z},\tilde{\alpha},\tilde{\beta},\tilde{\gamma})$ are initialized randomly, while the parameters ($z,\alpha,\beta,\gamma$) trained from stage 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 at Eq. (2.29) and (2.30), one can see that the terms contributed by each cohesive layer are linearly added to the overall compliance matrix or residual strain. This nice property enables the merge of two cohesive layers regardless of their order in the network.

In each cohesive network, the similarity between any pair of cohesive layers will be evaluated. Let us assume that their fitting parameters are given by ($\tilde{z}^{p},\tilde{\alpha}^{p},\tilde{\beta}^{p},\tilde{\gamma}^{p}$) and ($\tilde{z}^{p^{\prime}},\tilde{\alpha}^{p^{\prime}},\tilde{\beta}^{p^{\prime}},\tilde{\gamma}^{p^{\prime}}$). Since the cohesive layers are set to behave isotropically in the shear directions, they are said to be similar if they share the same normal direction. Mathematically, it requires that the absolute value of the 33 component in a matrix $\textbf{R}^{p\text{-}p^{\prime}}$ is close to 1:

[TABLE]

with $\textbf{R}^{p\text{-}p^{\prime}}$ defined as

[TABLE]

In this case, the two cohesive layers will be merged, and their new activations are

[TABLE]

The $p^{\prime}$-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 extrapolated 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. Simplification first arises from the assumption that the cohesive surface is isotropic. This indicates that the resistance to sliding is independent of the direction of sliding. and only depends on the magnitude $d_{S}$:

[TABLE]

The total opening displacement vector d and traction t take the following forms,

[TABLE]

where n is the unit normal to the cohesive surface. An effective opening displacement $d_{m}$ is further introduced to simplify the mixed-mode cohesive law, and the free energy density per unit undeformed area $\phi$ is a function of $d_{m}$,

[TABLE]

where q are some internal variables which describe the irreversible processes. Accordingly, the effective traction is defined as

[TABLE]

The compressive stress state should not cause the cohesive layer to fail, so that it does not contribute to the free energy density. In addition, no friction effect is included in the cohesive model. The tensile and compressive cases are considered separately:

For the tensile case $d_{n}\geq 0$, the effective opening displacement $d_{m}$ is

[TABLE]

where the positive parameter $\beta$ defines the ratio of effects from normal and shear displacements. The cohesive law becomes

[TABLE] 2. 2.

For the compressive case $d_{n}<0$, the effective opening displacement $d_{m}$ is

[TABLE]

The cohesive law can be written as

[TABLE]

As shown in Figure 4, a bilinear cohesive law is adopted in this work. In the elastic loading regime, the modulus of cohesive layer is $K$. In the application to be shown in Section 5, this interfacial modulus is chosen to be a large value compared to the moduli of other material phases, thus, the interface is nearly perfectly bonded at the beginning. The effective opening displacement before the softening regime is denoted by $d_{c}$, and the corresponding effective traction is $\sigma_{c}$. The opening displacement at full failure $t_{m}=0$ is $d_{f}$. The critical energy releasing rate $G_{c}$ is equal to the area under the curve of $t_{m}$-${d_{m}}$,

[TABLE]

The only internal variable is the effective opening displacement at maximum traction $d_{0}$. Initially, $d_{0}$ is equal to $d_{c}$,

[TABLE]

and it should decrease along the softening path when $d_{m}$ exceeds the current $d_{0}$. Mathematically, the effective traction $t_{m}$ of the bilinear cohesive law can be calculated by

[TABLE]

and

[TABLE]

In practice, to avoid singularity at full separation of cohesive layer, a linear term $\kappa d_{m}$ is added to the effective traction,

[TABLE]

where $\kappa$ is a tiny stiffness parameter chosen as $\kappa/K=1\times 10^{-6}$ in this work.

When reloading, the curve will simply trace the unloading path and rejoin the cohesive law at $d_{m}=d_{0}$. Meanwhile, an elastic response is assumed for the compressive case, with the same modulus as the one in the bilinear law,

[TABLE]

More details on the derivation of the tangent stiffness matrix and viscous regularization for implicit analysis in the online stage are provided in D.

5 Application to unidirectional fiber-reinforced composite

5.1 Problem settings

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-reinforced composite is briefly discussed in A. The geometry of the UD RVE is shown in Figure 5 (a). The volume fraction of the fiber phase is 29.4%. The average diameter of the fibers is around 2.5 mm, so that the characteristic length $L$, a hyper-parameter for the cohesive networks, is set to be

[TABLE]

In the FE mesh for DNS, there are 69048 nodes, 63720 8-node hexahedron solid elements, and 4320 8-node zero-thickness cohesive elements. All the material parameters used in the online stage are given in Table 2.

The fiber phase is assumed to be linear elastic, and its Young’s modulus $E_{1}$ and Poisson’s ratio $\nu_{1}$ are 500 GPa and 0.3, respectively. Two cases are considered for the matrix phase. In the first case, the matrix phase is linear elastic, with $E_{2}=100$ GPa and $\nu_{2}=0.3$. In the second case, the matrix phase is an elasto-plastic material, which obeys von Mises plasticity with piecewise linear hardening. The material yields at 0.1 GPa, and the yield stress $\sigma^{Y}_{2}$ is a function of the equivalent plastic strain $\varepsilon^{p}_{2}$:

[TABLE]

The equivalent modulus ($K\times L$=25000 GPa) of the interface is much larger than the moduli of fiber and matrix phases, indicating that the interface is nearly perfectly bonded at the beginning. Following Eq. (4.9), the effective opening displacement at complete failure $d_{f}$ can be computed from the critical effective traction $\sigma_{c}$ and the critical energy releasing rate $G_{c}$:

[TABLE]

To restrain the rate-dependency from the numerical stabilization in Eq. D.7, the viscosity-like parameter $\zeta$ (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 set of linear elastic DNS of the UD RVE under 6 orthogonal loading conditions is performed, which typically takes 182 s using LS-DYNA® RVE package on 10 Intel® Xeon® CPU E5-2640 v4 2.40 GHz processors 222The same machine is used for all other computations in this work.. With 1000 samples in total, the offline data generation process costs 50 h. This time can be further reduced by 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 each epoch. 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 libraries for matrix manipulation and linear algebra.

5.2 Offline training results

In training stage I, the original material networks with three different depths, $N=5,7,$ and $9$, are considered. For $N=5$ and 7, the networks were trained for 20000 epochs. Since the network with $N=9$ has more fitting parameters, it was trained for 40000 epochs. The training and test errors are listed in Table 3, where the error of sample $s$ is defined as

[TABLE]

To show the fractions of active nodes in the bottom layer, treemap plots of all three DMNs are also provided in Figure 6. For the network with $N=5$, there remains only one active node in the bottom layer for the fiber phase. Evidently, this is not sufficient for capturing the anisotropic behaviors of the UD RVE, and the average test error $\bar{e}^{te}$ is more than 10%. One important feature of DMN is that the average training error and the average test error of unseen samples are almost the same as shown in Table 3. With a small amount of data, it requires no regularization scheme (e.g. L2 penalty or Dropout) to avoid the overfitting issue, mostly thanks to the embedded physics in the building block. More detailed results of training stage I, 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, where the new fitting parameters are also optimized using the SGD algorithm. For all the cases, the networks were trained for 5000 epochs. Figure 7 (a) shows the histories of test error for $N=5$, 7, and 9, with the same number of cohesive layers in each cohesive network $N_{c}=4$. Their training and test errors at the end of learning 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 been investigated for depth $N=9$. We can see from Figure 7 (b) that one cohesive layer $N_{c}=1$ is not able to capture the interfacial effect accurately. While for $N_{c}=2$, 4, and 6, the test errors are all reduced to be less than 0.5%. After 5000 epochs, the network with $N_{c}=4$ reaches the lowest test error, instead of the one with $N_{c}=6$. This can be explained by the fact that the network with $N_{c}=6$ has more redundant cohesive layers, 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 stage I, in any measure, is always less than the one in stage II. This is physically sound as the UD RVE with deformable fiber-matrix interfaces is a more complex material system. Moreover, train stage II inherits the parameters ($z,\alpha,\beta,\gamma$) from stage I, and they are fixed in the training. Following this strategy, the accuracy of 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 is the volume fraction of its enriched node in the bottom layer of original material network, and the y value is its reciprocal length parameter $\tilde{v}$. The remaining number of active cohesive layers for $N=5,7,$ and 9 are 2, 13, and 40. Since the network with $N=9$ has more enriched nodes, the mean volume fraction of the nodes is smaller, and more cohesive layers are deactivated. For all the cases, the mean reciprocal length parameter is around 0.05 $\text{mm}^{-1}$. This value is related to the characteristic length of the UD RVE, whose average diameter of the fibers is 25 mm.

Offline training time: All training processes are parallelized using 10 processors. Here the computational times in stage II for networks with different depths are listed and compared. For $N=5$, networks with $N_{c}=2,4,6$ took 0.4, 0.5, 0.6 hours. For $N=7$, the training times are 1.7, 2.0, 2.2 hours. For $N=9$, the training times are 6.7, 7.9, 8.2 hours. Therefore, $N_{c}$ has less effect on the training time than $N$, indicating that most of the computational cost still locates at the original material network.

5.3 Online extrapolation to interfacial failure analysis

For the UD composite, three different loading directions are evaluated in the online stage: transverse tension and compression, transverse shear, and longitudinal shear. Here the case under longitudinal tension along 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 Figure 5 (b), while the overall stress components in the other directions are zero. For example, in transverse tension and compression, the loading conditions can be written as

[TABLE]

Accordingly, the absolute strain rate $|\dot{\varepsilon}_{11}|$ is 1.0 s*-1*. The time step in both DNS and DMN is $\Delta t=1.0\times 10^{-4}$. For simplicity, $\varepsilon$ and $\sigma$ are used to denote the overall strain and stress in this section. Snapshots of DNS are taken at $t=0.004,0.012,0.024,$ and 0.030 s to illustrate the evolution of RVE deformation during the loading.

For all the online examples to appear in this section, the study shows that $N_{c}$ has little effect on the stress-strain predictions as long as $N_{c}\geq 2$, while the depth of the original material network $N$ is more prominent. Therefore, $N_{c}$ is set to be 4 by default in the online DMNs. Ideally, a 3D material block can be fully covered by 3 cohesive layers if they are oriented in three orthogonal directions. In this regard, $N_{c}=4$ tends to be sufficient for the solution space to capture the interfacial effects, and more layers would not help to improve the accuracy too much. This statement on the choice of $N_{c}$ also applies to the particle-reinforced RVE shown in 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) for plastic matrix, no interfacial failure happens at $t=0.004$ s. Since the fiber phase is stiffer than the matrix phase, its average von Mises stress is higher. The interfaces normal to the loading direction start to debond at $t=0.005$ s and reach the maximum opening at $t=0.012$ s. Due to the debonding, less load transfers from the matrix to the fibers, and more stress concentration appears in the matrix phase. Around $t=0.024$ s, closure of the interfaces happens. Interestingly, further compression results in debonding of the interfaces in the loading direction, which is related to the Poisson’s effect. The corresponding stress-strain curves from DNS and material networks with $N=5,7,9$ and $N_{c}=4$ are provided in 9 (c). Stiffness degradation due to the failure of interfaces is observed during unloading. Other softening behaviors are also captured by DMN with cohesive layers. For $N=7$ and 9, the DMN predictions agree with the corresponding DNS results very well.

The difference between transverse tension and compression is more evident on the plot for elastic matrix shown in Figure 9 (b). Since there is no plasticity in the matrix, during unloading, the stress-strain curve almost goes back to the origin with a degraded stiffness. However, in compression, the interfaces are in contact, so that the stiffness returns to a high value and the response is linear elastic without any failure. It appears that the DMN slightly underestimates the stiffness under compression. Nevertheless, these physical phenomena 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 shows that 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 can see from Figure 10 (b), the reverse-loading part of stress-strain curve still undergoes some softening effect during the plastic yielding, which is consistent with those findings in the DNS snapshots. For either case with 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 in predicting 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 interfacial failure is driven by the sliding deformation, which is different from the previous two loading cases. Moreover, after reversing the loading direction at $t=0.012$ s, the sliding locations still stay the same, indicating that the amount of newly failed interfaces is small. This is reflected in the corresponding stress-strain curve plotted 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 almost linear elastic, and passes the origin. For both cases, the DMN predictions with $N=7$ and 9 agree with the DNS results very well.

Online prediction time: The DMN with cohesive layers learns an efficient reduced-order model of the DNS RVE with deformable interfaces. Both DNS and DMN models have 300 loading steps. A typical DNS with plastic matrix took 5.6 h on 10 processors. Yet it is noteworthy that the original DNS with perfectly bonded interfaces (or no cohesive element) only took 1.6 h on 10 CPUs, implying that the softening interfaces consume more iterations in the implicit analysis. In comparison, DMNs with cohesive layers took 1.4, 5.8 and 33.3 s on one CPU for with $N=5,7$ and 9, respectively. In terms of the CPU time, the enriched material network with $N=9$ and $N_{c}=4$ is more than 6000 times faster than the DNS.

5.4 Predictions on local properties

Local properties can also be predicted by the DMN. In the online stage, each active node in the bottom layer of DMN is treated as an independent microscale material, representing either the fiber or the matrix phase in the UD composite. Meanwhile, each cohesive layer can be regarded as a part of the fiber-matrix interfaces, 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, one can 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 active bottom-layer nodes and 40 active cohesive layers after the training. While for the DNS, there are 63720 solid elements and 4320 cohesive elements.

Figure 12 shows the distributions of various local fields predicted by DMN and DNS with plastic matrix under uniaxial tension at $t=0.004$ s, while the strain in the loading direction is $\varepsilon_{11}=0.004$. The corresponding global stress-strain curves are given in Figure 9 (c). For the matrix and fiber phases, the distributions of von Mises stress are compared in Figure 12 (a). Since most of the interfaces are still bonded, the average von Mises stress of the fiber phase is more than two times larger than that of the matrix phase. For the interfaces, the normal traction ranges from -0.15 GPa to 0.45 GPa, which is less than the critical effective traction $\sigma_{c}=0.5$ GPa. Furthermore, the distributions of local fields from DMN and DNS at $t=0.012$ s and $\varepsilon_{11}=0.012$ are presented in Figure 13. Due to debonding, the profile of the normal traction in Figure 13 (b) has a large density at [math] GPa, which corresponds to those failed interfaces in DNS or cohesive layers in 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 phase becomes more dispersed, and the average stress is close to the one in the fiber phase. Although the DMN has much less DOF, it gives a good estimation on the distributions of local fields comparing to the DNS reference.

6 Conclusions and future work

A new architecture of DMN with cohesive layers is proposed to consider deformable interfaces inside a heterogeneous material. In the physics-based cohesive building block, the activation of a cohesive layer is governed by a reciprocal length parameter, which also defines the effective thickness of the building block in its normal direction. Due to the existence of analytical solutions, the new material network can still be trained using gradient-based optimization methods and the backpropagation algorithm. Moreover, a two-stage training strategy is utilized to learn the fitting parameters in the original material network and the cohesive networks separately.

The predictive capability of the method is demonstrated through the debonding analysis of the UD composite 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 irreversible cohesive law. Though the material network is trained with linear elastic data, it can accurately capture the nonlinear RVE responses both locally and globally with complexity coming from the plastic matrix and the softening interfaces. Involving less DOF, the enriched material network with $N=9$ and $N_{c}=4$ is more than 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 an efficient way of realizing concurrent multiscale simulations. With the addition of cohesive layers, it provides an effective data-driven material model for the interfacial failure analysis of a large-scale heterogeneous structure. Since measuring interface properties directly from the experiments can sometimes be challenging, DMN can also be used for inverse prediction of those unknown properties. Finally, with the transfer learning technique proposed in [24], the unified microstructural databases can potentially be created for the design of 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 Huang for helpful discussions. Z. Liu warmly thanks Dr. John O. Hallquist for his support to this research. This is 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 cohesive layers is applied to a particle-reinforced composite. Volume fraction of the particle phase is 22.6%. In terms of 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, while its particle phase shares the same properties with the fiber phase of the UD composite. The characteristic length $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, which are trained for 20000, 20000, and 40000 epochs, respectively. In training stage II, the number of cohesive layers in each cohesive network is $N_{c}=4$, and all three enriched DMNs are trained for 5000 epochs. The training and test errors of both stages are listed in Table 4. With depth $N=8$, the average training and test 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 14. The loading history of $\varepsilon_{11}$ is shown in Figure 5 (b), while the RVE is stress-free in the other directions. As shown in the DNS snapshot at $\varepsilon_{11}=0.012$, the particle-matrix interfaces are debonded in the loading direction. The matrix phase is made partially transparent to show the inner separated surfaces. Furthermore, the stress-strain curves predicted by DNS and DMNs are plotted in Figure 14 (b). Due to the contact of debonded interfaces, the particle-reinforced RVE is stiffer under compression. A good agreement is observed between DNS and DMNs with $N=6$ and 8.

Appendix B Elementary rotation matrices

In small-strain formulation, the elementary rotation matrices shown in Eq. (2.19) are given by

[TABLE]

[TABLE]

[TABLE]

The in-plane and output-plane rotation matrices $\textbf{r}^{p}$ and $\textbf{r}^{v}$ for an arbitrary angle $\theta$ are defined in Mandel notation as

[TABLE]

Appendix C Finite-strain solution of cohesive building block

For finite-strain problem, the deformation gradient $\utilde{\textbf{F}}$ and first Piola-Kirchhoff stress $\utilde{\textbf{P}}$ are chosen as the strain and stress measures, respectively. For 3-dimensional problems, $\utilde{\textbf{F}}$ and $\utilde{\textbf{P}}$ are vectorized to simplify the computation:

[TABLE]

[TABLE]

In the generic cohesive building block, the constitutive relation of material 0 in finite-strain formulation can be written as

[TABLE]

where $\textbf{B}^{0}$ and $\delta\textbf{F}^{0}$ are the finite-strain compliance matrices and residual deformation gradient, respectively. The overall constitutive relation of the building block is

[TABLE]

Based on the equilibrium condition, B and $\delta\textbf{F}$ are derived to be

[TABLE]

and

[TABLE]

The nonzero terms in $\tilde{\textbf{G}}^{f}$ and $\delta\tilde{\textbf{d}}^{f}$ are

[TABLE]

The compliance matrix G and residual displacement $\delta\textbf{d}$ of the cohesive layer are defined in Eq. (2.9). The definition of the finite-strain rotation matrix $\textbf{R}^{f}$ can be found in [41], which can also be decomposed into three 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 cohesive law needs to be derived, which is defined as

[TABLE]

Again, the tensile and compressive cases are treated differently:

For tensile case $d_{n}\geq 0$, the components of t are

[TABLE]

The independent components in K are

[TABLE]

[TABLE]

with

[TABLE] 2. 2.

For compression case $d_{n}<0$, the components of t are

[TABLE]

so that

[TABLE]

[TABLE]

The tangent stiffness matrix of the bilinear cohesive law may not be positive definite. At the point of failure, implicit calculations for both DNS and deep material network are unable to converge to an equilibrium solution, which makes it impossible to track the post-failure behavior. Therefore, an extra viscous term is introduced to overcome the convergence problem following the work by Gao et al.[43]:

[TABLE]

where $\zeta$ is a viscosity-like parameter that governs viscous energy dissipation. By selecting a sufficiently small time increment in the numerical scheme, one is able to make the tangent stiffness tensor of the softening material positive definite.