Cell division in deep material networks applied to multiscale strain localization modeling
CCell division in deep material networks applied to multiscale strainlocalization modeling
Zeliang Liu a, ∗ a Ansys Inc., Livermore, CA, USA
Abstract
Despite the increasing importance of strain localization modeling (e.g., failure analysis) in computer-aidedengineering, there is a lack of effective approaches to consistently modeling related material behaviors acrossmultiple length scales. We aim to address this gap within the framework of deep material networks (DMN)- a physics-based machine learning model with embedded mechanics in the building blocks. A new celldivision scheme is proposed to track the scale transition through the network, and its consistency is ensuredby the physics of fitting parameters. Essentially, each microscale node in the bottom layer is described by anellipsoidal cell with its dimensions back-propagated from the macroscale material point. New crack surfacesin the cell are modeled by enriching cohesive layers, and failure algorithms are developed for crack initiationand evolution in the implicit DMN analysis. Besides single material point studies, we apply the multiscalemodel to concurrent multiscale simulations for the dynamic crush of a particle-reinforced composite tubeand various tests on carbon fiber reinforced polymer composites. For the latter, experimental validationson an off-axis tensile test specimen are also provided.
Keywords:
Deep learning, multiscale modeling, failure analysis, damage and fracture, composites
Contents1 Introduction 2 ∗ Corresponding author.
Email address: [email protected] (Zeliang Liu)
Preprint submitted to Elsevier January 19, 2021 a r X i v : . [ c s . C E ] J a n .1.2 Effects of relaxation time and network depth . . . . . . . . . . . . . . . . . . . . . . . 214.2 Dynamic crush of a composite tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224.3 Three-point bending tests with microstructure anisotropy . . . . . . . . . . . . . . . . . . . . 244.4 Composite coupon tests with experimental validation . . . . . . . . . . . . . . . . . . . . . . . 28 Over the past decades, multiscale modeling has become a key enabling technology for computer simulationand materials design in the aerospace, automotive, and consumer electronics industries. When dealingwith materials with expressive microstructures, engineers often find a single-scale phenomenological modelfails to capture the anisotropic nonlinear behaviors due to complex physical interactions across scales. Incontrast, the core of multiscale materials modeling stands on explicit microstructure representations andhomogenization schemes. For elastic and nonlinear plastic materials with scale separation, theories forhomogenization are well established. Most recent studies have been focusing on improving the efficiencyand accuracy of microstructure modeling. However, for materials modeling with strain localization, such asdamage and fracture, the homogenization condition is closely entangled with the microstructure’s length-scale variations, making it challenging to define a consistent multiscale model. In the paper, we will attemptto address this gap within the framework of Deep Material Network (DMN) [1, 2], which is a machinelearning model with physics-based building blocks.Meanwhile, the rise of Machine Learning (ML) has been continuously advancing the frontier of materialsmodeling and multiscale simulations. Although feedforward neural networks have been applied to materialsmodeling and characterization since the 1990s [3, 4], their popularity was limited mainly due to the lackof computational resources. Recent advancements of computer hardware systems and open-source MLplatforms stimulate the renaissance of deep learning. They allow researchers to effectively explore moresophisticated network architectures, such as convolutional neural networks [5] and recurrent neural networks[6, 7]. Importantly, the successes of these customized networks in computer vision and natural languageprocessing also motivate mechanicians to design new ML architectures specialized for materials modelingand multiscale simulations, which could potentially discover hidden physical relationships and outperformexisting knowledge-driven models in certain tasks.A category of knowledge-driven multiscale models is based on micromechanics theories, such as Mori-Tanaka method [8, 9] and self-consistent methods [10, 11]. Given the elegance of Eshelby’s solutions [12, 13]and mean-field schemes, they are widely adopted in current research and industrial applications. A majorbarrier to using micromechanics theories for predicting highly nonlinear material behaviors is that theanalytical solutions for idealized geometries cannot capture localized deformations, especially when it comesto failure analysis. Therefore, those predictions often lose the multiscale characteristics and require extracalibrations on the empirical model parameters. In contrast, the Representative Volume Element (RVE)-based methods simulate microstructures explicitly, and the discretized models are usually solved by FiniteElement (FE) methods [14, 15] or fast Fourier transform (FFT)-based methods [16]. Couplings betweenthe microscale RVE and the macroscale structure enable the so-called concurrent multiscale simulations,referred to as FE [17, 18] or FE-FFT [19] depending on the model type at each scale.2y formulating the microstructure as a boundary value problem, RVE-based methods offer more flex-ibility on the choices of material constitutive laws, including anisotropic elasticity, plasticity with historydependency, and nonlinear hyperelasticity under large deformations. However, there has been less consen-sus on using RVE for failure analysis and other straining localization problems, although they have becomeincreasingly important for multiscale modeling in computer-aided engineering. The uniqueness of failureanalysis originates from the loss of ellipticity when local damage or fracture appears in the RVE. As theequilibrium condition becomes unstable, the deformations tend to localize into a single layer of elements.In order to ensure the objectivity with respect to mesh size and shape in the RVE analysis, numericaltreatments are needed to regularize the ill-posed boundary value problem, and a few representatives aregradient-enhanced or non-local damage models [20, 21, 22, 23, 24], phase-field fracture [25, 26], crack bandtheories [27, 28], and cohesive zone models [29, 30, 31]. In particular, the cohesive zone models have beencoupled with extended/generalized finite element methods for crack growth modeling [32, 33]. The crackband model [27] regularizes the post-failure softening stiffness by the mesh size. Similar approaches havebeen widely adopted in commercial software due to their simplicity and effectiveness.Although the spurious mesh dependency could be resolved by the aforementioned numerical techniquesat a single scale, key challenges emerge when the multiscale coupling is considered. As a strain localizationband (e.g., a crack or a damage zone) cuts through the RVE and breaks the periodicity, prescribing propermacroscale boundary conditions becomes much more difficult. Besides, the homogenization results dependon the microscale RVE size, which does not necessarily match the macroscale length scale of localization.This mismatch will result in energy inconsistency of scale transition as the homogenized bulk stress-strainresponses are directly passed to the macroscale material point. Namely, if the RVE size is smaller thanthe macroscale length scale of localization (e.g., the mesh size), the macroscale energy dissipation would beover-predicted. In Section 1.2, we will discuss the problems for strain localization modeling in detail.Another challenge facing RVE-based methods is the efficiency issue since the direct numerical simulations(DNS) often contain thousands of degrees of freedom (DOF). It becomes the bottleneck for materials designand concurrent multiscale simulations, both of which require a recurrent assessment of the RVE model.In this regard, many ML methods and frameworks have been proposed to improve the materials models’efficiency by learning from experience or data. The training data can come from either physical experimentsor simulations, while the simulation data are usually easier to be acquired. One way of mining the data isto perform model reduction on the DNS model of RVEs, such as proper orthogonal decomposition (POD)[34, 35, 36], self-consistent clustering analysis (SCA) [37, 38, 39]. In comparison, methods have also beenformulated based on the stress-strain data – e.g., Gaussian process modeling [40, 41], model-free approaches[42, 43, 44, 45], feedforward neural networks [3, 46, 47, 48], recurrent neural networks [49, 50, 51, 52]. Weconsider DMN [1, 2] as a blend of both approaches: it finds a reduced-order representation of the DNSRVE model while the training only requires stress-strain data. Moreover, the transfer learning approach ofDMN for creating a unified microstructure database is proposed in [53] and applied to short-fiber reinforcedcomposites in [54]. The thermodynamic consistency of DMN has also been studied in [55].A few ML or data-driven methods have been applied to multiscale failure analysis or strain localizationmodeling in general. Oliver et al. [56] performed the model reduction via the POD method for multiscalefracture problems. A domain decomposition strategy was proposed to treat the regular domain and thelocalization band separately in the RVE. Liu et al. [38] applied the SCA method to multiscale damageanalysis, and the microscale damage parameters are calibrated to restore the energy consistency. Wang et al.[49] used recurrent neural networks to generate cohesive laws for modeling localized physical discontinuitiesat different length scales of multi-permeability porous media. Recently, DMNs with enriching cohesive layerswere developed for interfacial failure analysis [57], like debonding in composites. Although no localizationissue has been considered, the work establishes the basis of failure modeling within the DMN framework. When it comes to failure analysis or general strain localization problems, RVE-based methods encountercertain difficulties that have long been recognized by researchers in the field [15, 58, 59, 60]. To help explainthe issues, we introduce three length parameters involved in a typical multiscale simulation based on RVEmodeling: the macroscale length parameter h , the RVE size l RV E , and the microscale localization size3 … ℎ !" 𝑙 $ ℎ = 𝑙 $ = 𝑙 !" Microscale RVE ℎ Macroscale (a) Limitation on the RVE size … … Microscale RVE 𝑙 !" ℎ Macroscale ℎ = 𝑙 $ ≠ 𝑙 !" What is 𝑙 $ ? (b) Difficulties of applying boundary conditionsFigure 1: Issues of RVE modeling when localization is presented. The macroscale length parameter, the RVE size, and thelocalization size are denoted by h , l RV E , and l c , respectively. In (a), h is equal to the element size, while In (b), h is set bythe size of the nonlocal regularization scheme. Periodic boundary conditions are applied on the RVEs. l c . The element size usually determines the macroscale length parameter h because the softening regionlocalizes into a single layer of elements, as shown in Figure 1 (a). Alternatively, if one uses nonlocal orgradient-based regularization to reduce the mesh sensitivity, h is set by the averaging size in the numericalscheme. The RVE size l RV E is the edge length of the rectangular cell. The microscale localization size l c isthe homogenization length scale of a localization band in the microscale model. For example, if the periodicboundary conditions are applied on the RVE model, l c is equal to the vertical distance between the repeatingcracks. Therefore, depending on the modeling approach and the localization configuration, one can see that l c is not necessarily equal to the RVE size l RV E .Here, we summarize three main issues of RVE analysis that limits its applications to large-scale multiscalesimulations with localization:1. The optimum RVE size changes with macroscale length parameter h and localization configuration(e.g., crack orientation); otherwise, proper energy regularization must be introduced for a predefinedRVE size.2. The boundary conditions on a rectangular RVE model are not well defined due to the change ofperiodicity induced by localization.3. Solving the RVE with local material damage or softening can be time-consuming.For energy consistency, the localization size l c should be equal to the macroscale length parameter h , l c = h . For the crack configuration shown in Figure 1 (a), l c is equal to the RVE size. As a result, the optimumRVE size is set by the macroscale length parameter h . However, in a practical multiscale simulation, theelement size of the macroscale model is usually much larger than the microstructural characteristic length, h (cid:29) l RV E . Imposing l RV E = h would make the RVE model too expensive to solve, defeating the purposeof using multiscale simulations. Meanwhile, a single RVE with the predefined size is usually not sufficientbecause the optimum l RV E may change with the element size in the macroscale model or even the crackorientation (see Figure 1 (b)), which is not known before running the simulation.On the other hand, applying proper boundary conditions on a rectangular RVE model is questionable asthe localization surface (or band in 2-D space) breaks the periodicity of RVE in the normal direction to thesurface. For this reason, the existence of RVE under localization is often argued in the literature. As thescales are no longer separated, the homogenized stress-strain responses depend on the RVE size l RV E , whichis contradictory to the original definition of an RVE [10]. On the other hand, the boundary conditions wouldaffect the localization configuration in a non-physical way. For example, as shown in Figure 1 (b), the crackis distorted numerically to satisfy the periodic boundary conditions on the RVE model, while physically, itshould propagate straightly in 30 ◦ . To resolve this issue, one needs to adapt the aspect ratio of the RVE4ith respect to the crack orientation. Similar to the choice of RVE size, this becomes a “Chicken-and-Egg”problem since the crack orientation is unknown before analyzing the RVE.Last but not least, the appearance of material softening puts extra burdens on the numerical solver,while RVE models based DNS methods can be rather time-consuming, even for well-posed problems. Ifone intends to solve the RVE problem implicitly, extra numerical treatments like viscous regularization areoften required to make local tangent stiffness matrices positive-definite for the convergence of the nonlinearimplicit solver. Nevertheless, more iterations are needed for solving the RVE problem in the course oflocalization so that improving the model efficiency becomes more prominent. We will introduce a cell division scheme for consistent scale transition in the DMN framework to addressthe aforementioned issues in strain localization modeling. New algorithms of crack activation and evolutionwill be developed for the implicit failure analysis. By coupling DMNs with macroscale finite element models,we will demonstrate the applications of concurrent multiscale simulations to a broad range of materialsystems.The remainder of this paper is organized as follows. In Section 2, we propose the idea of cell divisionfor tracking the length scales inside a DMN. Analytical solutions of the two-layer building block are derivedunder the scenario of dividing an ellipsoidal (or elliptic in 2-D space) cell. The cell-division results on thebinary-tree networks are demonstrated for two representative microstructures. Section 3 focuses on DMN-related algorithms for failure analysis, including the cohesive law with viscous regularization, activationsof potential crack surfaces, and the overall implicit solution scheme. Section 4 provides examples rangingfrom single material point studies to concurrent multiscale simulations of nonlinear anisotropic composites.Advantages and limitations are discussed in Section 5. Conclusions and future work are summarized inSection 6.
2. Scale transition in DMN: Concept of cell division
Deep Material Network (DMN) is a mechanistic machine learning method for modeling materials acrossdifferent scales. Its key features include the physics-based building block with interpretable fitting parame-ters, extrapolation capability for material and geometric nonlinearities, and efficient inference with a smallnumber of DOF. The data-driven framework of DMN starts from the offline stage, where linear elastic DNSare performed to generate the training data. Gradient-based optimization algorithms are then used to findthe optimum fitting parameters. In the online prediction stage, the network with trained fitting parame-ters can be extrapolated to consider nonlinear material responses based on the mechanistic building block’sanalytical solutions. A more detailed description of the data-driven framework is provided in Appendix A.DMN has been applied to several representative multiscale material systems, including particle-reinforcedrubber composites, short or continuous fiber-reinforced composites, and polycrystalline materials. Formu-lated based on the network structure, DMN takes a different route of describing multiscale material behaviorsfrom RVE analysis, which solves a boundary value problem of partial differential equations. Importantly, itoffers the potentials of overcoming the aforementioned issues in RVE-based methods (see Section 1.2): theanalytical solutions of DMN building block avoid the difficulties of defining proper boundary conditions ona rectangular unit cell, and the reduced-order DMN model is more efficient than a full-field RVE analysis.However, to fully tackle the multiscale localization problem, a consistent scale transition scheme still needsto be developed within the framework.This motivates us to introduce the concept of “cell division”. Essentially, each node in the network isassociated with a so-called “cell” dedicated to encoding the corresponding node’s length scales. The cellof the top node of DMN represents macroscale material point, so that its dimensions are defined by themacroscale length parameters h . In a backward propagation process shown in Figure 2, cells are dividedhierarchically following the geometric structure of each building block, and the length scale information ofthe microscale DOF can be extracted from the cells at the bottom layer of the network.5 = 45°, 𝑓 ! = 0.25 𝜃 = 90°, 𝑓 ! = 0.2𝜃 = 0°, 𝑓 ! = 0.5 𝐀 " = 1.00 1.00ℎ = 2 Micro-cell 1
Cohesive layer 1Activation of a cohesive layer with normal 𝐧 $ 𝐧 $ 𝐀 ’! =8.50 −7.50−7.50 8.50 𝐀 ’’ =1.39 −0.39−0.39 1.39𝐀 (( =32.11 −0.39−0.39 1.39𝐀 (! =8.50 −7.50−7.50 14.15 𝐀 () =2.11 −0.39−0.39 1.39 Micro-cell 2 Micro-cell 3 Micro-cell 4
Cohesive layer 𝑝 … Cell 4 𝑣 = 12 𝓣𝓡 * 𝐧 $ with 𝓡(𝓣) +’ 𝓡 * = 𝐀 () Macro-cell ( ℎ = 2 ) Micro-cell 4: 𝐀 (’ =8.50 −7.50−7.50 14.15 Figure 2: Illustration of the cell division process for scale transition in 2-D space. The macroscale cell is represented by anellipse defined by the “scale tensor” A macro , which is related to the macro scale parameter h . In the back propagation process,each cell is divided based on the interface orientation θ and phase volume fraction f of the corresponding two-layer buildingblock. The nodes in the bottom layer are regarded as the DOF of DMN, whose length scales are defined by the micro-cells.Cohesive layers can be activated within a micro-cell to simulate the failure. To simplify the mathematical formulation of the cell division process, we assume each cell to be anellipsoid in 3-D space, or an ellipse in 2-D space. The implicit equation of an ellipsoidal cell can be expressedin a 3-D Cartesian coordinate system ( x, y, z ) as (cid:2) x y z (cid:3) A A A A A sym A xyz = 1 , (2.1)where we refer A as the “scale tensor” in the paper. The second-order tensor A is positively definite andsymmetric. Specifically, if a node possesses a spherical domain with diameter h , its scale tensor can bewritten as A = /h /h sym /h . (2.2)Accordingly, the cell division process for scale transition in DMN becomes the back-propagation modelingof the scale tensors. As illustrated in Figure 2 for a 2-D DMN with three network layers ( N = 3), thetop node, which represents the macroscale material point, has a scale tensor A macro , which is commonlydetermined by the element size or the nonlocal averaging size in the macroscale model. Based on theinterface’s orientation and the phase fraction in each building block, the mother node’s cell is divided intothe two child nodes’ cells. As information of the scale tensor propagates from the top layer to the bottomlayer, A macro → A , A → A , A , A , A , the scales of the micro-cells can be properly tracked. The celldivision’s mathematical formulation in a generic two-layer building block will be presented in Section 2.2.Crack planes can be activated in each micro-cell to model the microscale material failure. In thiswork, a crack plane is modeled by the cohesive-layer approach. The orientation of the cohesive layer, or thecorresponding crack plane, is determined by the base material’s current deformation state and the associatedfailure criteria. Once a cohesive layer is activated and attached to the DMN node, its constitutive behaviorswill be described by a traction-separation law. The scales of the micro-cell come into play when one needs6o determine the ”reciprocal length scale parameter” [57], v , which can be interpreted as the inverse ofthe effective thickness of the base material normal to the cohesive layer. Physically, it characterizes thecontribution of the separation displacement to the overall strain of the micro-cell. We will show how toderive the expression of v as a function of the micro-cell’s scale tensor A jN and the normal of crack plane n c in Section 3.2. Remark 1.
In our previous work for materials with predefined deformable interfaces [57], the orientationsand reciprocal length parameters of the cohesive layers are the fitting parameters optimized based on DNStraining data, and they will be kept the same in the prediction stage after training. By contrast, the cohesivelayers in this work are utilized to model the crack planes whose orientations and length scales depend on thebase material’s current deformation state during the online analysis. Therefore, these geometric parametersare not pre-determined from offline training.
Once the rotation angles and the reciprocal length parameters of the cohesive layers are determined, theenriched DMN can be solved through implicit analysis. At each iteration of Newton’s method, the stiffnesstensor C N and residual strain δ σ N of an enriched bottom-layer node/cell can be computed by C N = (cid:34) D base + N c (cid:88) p =1 v pc R pc ˜ G p ( R pc ) − (cid:35) − (2.3)and δ σ N = − C N (cid:34) δ ε base + N c (cid:88) p =1 v pc R pc δ ˜ d p (cid:35) , (2.4)where the compliance tensor D base and the residual strain δ ε base are obtained by evaluating the constitutivemodel of the base material. ˜ G p and ˜ d p are the compliance matrix and residual displacement vector comefrom the traction-separation law of the p -th cohesive layer. Note that these residual quantities appear onlywhen the system experiences geometric or material nonlinearities, such as plasticity in the base materialand softening in the cohesive layer. While this work will mainly focus on small-strain formulation, a shortdiscussion on the finite-strain formulation with geometric nonlinearity is provided in Appendix B. Moreover, N c denotes the total number of enriching cohesive layers in the base material. If the cohesive layer is assumedto behave isotropically in the crack plane (see Section 3.1), the rotation matrix R pc for the p -th cohesivelayer can be solely determined from its normal vector n pc .The information of stiffness tensors and residual stresses at the bottom layer are then propagated throughthe network all the way to the top node, which stores the macroscopic quantities: C macro and δ σ macro .Afterward, the macroscale boundary conditions are analyzed at the top node to fill any undefined macro-scopic stress or strain components. For example, in a 3-D concurrent multiscale simulation, the top nodereceives all six strain components of a material/integration point ∆ ε macro in the macroscale FE model, andthe missing stress components can be obtained as∆ σ macro = C macro ∆ ε macro + δ σ macro . (2.5)After the stress and strain information are propagated back to the base material and enriched cohesivelayers of each micro-cell at the bottom layer, the convergence is checked by comparing the updated microscaleincremental strains and displacement vectors with ones from the last iteration. Outside the loop for Newton’smethod, we also need algorithms to track the crack configuration and evaluate its convergence accordingly.The complete description of the solution scheme will be provided in Section 3.3. In a generic two-layer building block, we denote the scale tensor of the mother node as A , and onesof the child nodes as A and A . The unit normal to the two child materials’ interface is n . The volumefraction of the first child node is f , and the other child node has f = 1 − f . For a given ellipsoidal cell7 . . . − . . . n = [0 , T . . . . . . f − . . . − . . . n = [ √ / , √ / T . . . . . . f (a) A = (cid:20) . − . − . . (cid:21) . − . . . − . . . n = [1 , T . . . . . . f − . . . − . . . n = [ √ / , −√ / T . . . . . . f (b) A = (cid:20) . − . − .
39 1 . (cid:21) .Figure 3: Divisions of the mother node in a 2-D building block for different combinations of the scale tensors A and thenormal of the cutting plane n . Specifically, A s are picked from the two nodes at the second network layer in Figure 2. Thecell of the first child node with f = 0 . with a shape tensor A , we define the cutting surface, Λ( A , n ), as the intersection between the ellipsoid anda cutting plane, which has a normal n and passes through the center of the ellipsoid. Meanwhile, V ( A ) isthe volume of the ellipsoid.In our formulation, three conditions need to be satisfied during the cell division process:1. As centered at the origin, child cells are always within the mother cell.2. The mother and child cells share identical cutting surfaces with n normal to the building block’sinterface, Λ( A , n ) ∼ = Λ( A , n ) ∼ = Λ( A , n ) . (2.6)3. Consistency of volume fractions: f = V ( A ) /V ( A ) , with V ( A ) = V ( A ) + V ( A ) . (2.7)As we will show later in this section, the simplicity of the two-layer building block not only enables analyticalhomogenization of the mechanical quantities [1, 2] but also allows us to derive analytical functions of thecell-division process for scale transition under the above conditions.Figure 3 presents the cell-division results under four sets of A and n in 2-D space. For each case, theelliptic cells of the first child node are plotted for f varying from 0.0 to 1.0, with color ranging from red toblue. As f approaches 0, the cell of the first child node essentially falls onto the cutting surface and hasan infinitesimal thickness in the normal direction. One can also see that the child cells are always containedinside their mother cell, which is ensured by the first division condition. Remark 2.
The cell division for scale transition is analogical to the de-homogenization process for back-propagating physical quantities, such as stress and strain tensors. As the rotation of the two-layer buildingblock is performed at a separated step, the normal n of cutting plane in the cell-division model is always { , , } T in 3-D or { , } T in 2-D space. However, for generality, we develop the division model for anarbitrary normal direction. In Section 3.2, the solution will be reused to determine the reciprocal lengthparameter of a crack surface. Driving the formulations of A and A for A of an arbitrary ellipsoidal cell is not trivial. However,the process becomes much more straightforward for a unit sphere, because two principal axes of its childellipsoid are in the cutting plane, while the third axis is normal to the plane. Based on this observation, wepropose the following procedure for the derivation: 8. Rotate and scale the mother ellipsoid with A to a unit sphere with ˆ A = I . The transformation alsoapplies to the cutting plane, resulting in a new normal ˆ n .2. In the transformed configuration, cut the unit sphere based on ˆ n . Compute the scale tensors ˆ A andˆ A .3. Transform ˆ A and ˆ A back to A and A in the original configuration.We will derive the expressions of A and A for the 3-D building block, while the same procedure canalso be applied to the 2-D building block [1]. Through the eigenvalue analysis of A , we can find a rotationmatrix R which aligns the principal axes of the rotated ellipsoid A (cid:48) with the coordinate axes, A (cid:48) = λ λ λ = R T A R , (2.8)where λ , λ , and λ are the eigenvalues of A . For the rotated ellipsoid, the square root of λ i is equalto the inverse of the length of its i -th semi-axis. Therefore, another scaling operation is prescribed on therotated ellipsoid to transform it to a unit sphere, and the scaling matrix is T = / √ λ / √ λ / √ λ s.t. ˆ A = T R T A RT = I . (2.9)Let n s and n t denotes two orthogonal vectors in the original cutting plane, where we have n s × n t = n .After the transformation, the new vectors in the plane areˆ n s = T − R T n s , ˆ n t = T − R T n t . (2.10)The new normal to the transformed cutting plane isˆ n = ˆ n s × ˆ n t | ˆ n s × ˆ n t | = T R T n | T R T n | . (2.11)After the unit sphere is cut by the plane with normal ˆ n , the scale tensors of the two child nodes areobtained as ˆ A = I − (cid:18) − f ) (cid:19) ˆ n ⊗ ˆ n , ˆ A = I − (cid:18) − − f ) (cid:19) ˆ n ⊗ ˆ n , (2.12)which satisfy all three division conditions in the transformed configuration.The inverse transformations of ˆ A and ˆ A return the expressions for A and A in the original config-uration, A = RT − ˆ A T − R T , A = RT − ˆ A T − R T . (2.13)By substituting terms in Eq. (2.13) with Eq. (2.11) and (2.12), we arrive at A = A − (cid:18) − f ) (cid:19) n ⊗ n | T R T n | , A = A − (cid:18) − − f ) (cid:19) n ⊗ n | T R T n | . (2.14)Since both of the inverse rotation and scaling operations are affine transformations, A and A also satisfyall three division conditions in the original configuration.Meanwhile, the area of Λ( A , n ) can be calculated as S ( A , n ) = π | ˆ n s × ˆ n t | = π √ det A | T R T n | . (2.15)Implied by the second division condition in Eq 2.6, S ( A , n ) = S ( A , n ) = S ( A , n ) . (2.16)9 𝑓 = 22.6% 𝑁 = 4, 𝑁 !" = 4 𝑁 = 6, 𝑁 !" = 13 𝑁 = 8, 𝑁 !" = 28 (a) (b)(c)2 13 Figure 4: 3-D particle-reinforced composite: (a) RVE Geometry for finite element analysis, the volume fraction of the particlephase is 22.6%; (b) Tree-map plots of trained DMNs with depth N = 4, 6, and 8. The number of DOF, N dof , is shown on topof each plot; (c) Micro-cells arranged by the corresponding treemap plot with A macro = I . Micro-cells of the particle phaseare colored by red, and ones of the matrix phase blue. The macro-cell is plotted in gray at the center of each plot. This paper mainly focuses on two representative microstructures: a 3-D particle-reinforced composite anda 2-D composite with identical circular inclusions embedded in the matrix phase. For the 2-D composite,we will rely on the physical interpretations of DMN fitting parameters to directly transfer it to a 3-Dunidirectional-fiber composite, so that no extra training is required for generating the 3-D model. For bothcases, periodic boundary conditions are used in the RVE analysis based on finite element methods.The geometry of the 3-D particle-reinforced composite is shown in Figure 4 (a), and there are fouridentical spherical particles embedded in the matrix. The volume fraction of the particle phase is 22.6%. Toimprove the DNS training data’s accuracy, we mesh the RVE by 10-node tetrahedron finite elements. Theresulting DNS model has 84,693 nodes and 59,628 elements. Tree-map plots of the DMNs with differentdepth after training are given in Figure 4 (b). The number of DOF N dof is counted as the number of activenodes in the bottom layer. For N = 4, 6, and 8, we have N dof = 4, 13, and 28, respectively. More detailsabout the data generation and DMN training process can be found in Appendix A and our previous paper[2]. The micro-cell of each DOF can be obtained from the cell division process, and we place each micro-cellat the center of its corresponding block in the tree-map plot as shown in Figure 4 (c). For each case, themacro-cell has a scale tensor A macro = I , which appears as the gray sphere in the plot.The 2-D composite has identical circular inclusion embedded in the matrix with the volume fractionequal to 50%, as shown in Figure 5 (a). After discretization, the DNS model has 199,014 FE nodes and198,212 4-node 2-D plane strain elements, which guarantees that there are at least three layers of elementsbetween any two inclusions. The depth of DMN is set to 8, and there are 31 DOF after training. The 2-Dconfigurations for two A macro tensors are provided in Figure 5 (b). The 3-D configurations of micro-cellsafter the model transfer are shown in Figure 5 (c).Before diving into the details of model transfer from 2-D to 3-D composites, let us first look at the fittingparameters of the 2-D DMN with depth N : activations z j =1 , ,..., N − (2 D ) and rotation angles θ k =1 , ,..., i − i =1 , ,...,N .Physically, max(0 , z j D ) returns the weight of the j -th node at the bottom layer. Through forward summationof the weights, the phase fraction of every building block in the network can be determined. Meanwhile,the angle θ ki =1 controls the rotation or interface orientation of the k -th building block in the i -th layer. Asthe 3-D unidirectional fiber composite can be generated by stretching the 2-D model in the out-of-planedirection, there is no need to change the network topology in the plane. All the activations remain the same: z j (3 D ) = z j (2 D ) , for j = 1 , , ..., N − . (2.17)10 !" = 1 00 1 𝐀 !" = 2 00 1𝑣𝑓 = 50.0% 𝐀 !" = 1 0 00 1 00 0 1 𝐀 !" = 2 0 00 1 00 0 0.5 𝑁 = 8, 𝑁 &%’ = 31
Figure 5: 3-D unidirectional fiber composite from DMN model transfer. (a) The 2-D RVE with circular inclusions is modeledby finite elements, and the volume fraction of the inclusion phase is equal to 50.0%. After training, the 2-D DMN with 8layers has 31 DOF as shown in the treemap plot. The 3-D DMN models of the unidirectional fiber composite are createdby transferring the fitting parameters without extra training. (b) 2-D micro-cell configurations for two different macro scaletensors. (c) 3-D micro-cell configurations for two different macro scale tensors after model transfer.
For the 3-D rotation angles, we have α ki = θ ki , β ki = 0 , γ ki = 0 for i = 2 , , .., N and k = 1 , , ..., i − . (2.18)Moreover, to align the fibers along the third axis as shown in Figure 5(a) , we treat the rotation at the topnode ( i = 1) differently: α = θ , β = π , γ = 0 . (2.19) Remark 3.
It is possible to train the 3-D DMNs directly from a 3-D DNS model of the unidirectional fibercomposite. In fact, we have demonstrated this approach in [2]. However, in general, the 2-D DNS modeltakes less time for data generation, and the 2-D DMN also trains faster as it has fewer fitting parametersgiven the same depth. We realize it helpful to present this nice feature enabled by the physics-based buildingblock of DMN, while a complete discussion of model transfer is beyond the scope of this paper.
By the nature of model transfer in Eq. (2.17), the network configuration and weights have not beenchanged. Therefore, the 2-D and 3-D networks share the same treemap plots in Figure 5 (b). Meanwhile,since the cutting planes in the cell-division process are always parallel to the third axis, one can see fromFigure 5 (c) that all the micro-cells have the same dimension in the fiber direction as the macro-cell.
3. Failure Algorithms
As discussed in Section 2.2, crack planes are activated in the micro-cells to simulate the failure behaviorsin this work. The crack plane is modelled by a cohesive-layer model based on one-dimensional effectivetraction-displacement law [29, 30]. The total opening displacement vector d and the traction t can bewritten as d = d n n c + d S , t = t n n c + t S , (3.1)where n c is the unit normal to the crack plane. Different from n of the interface in the two-layer buildingblock, n c is not determined from offline training. Instead, as will be shown in Section 3.2, it depends11n the deformation state of the corresponding micro-cell during the online analysis. An effective openingdisplacement d m is further introduced to simplify the mixed-mode cohesive law, and the free energy densityper unit undeformed area φ becomes a function of d m , φ = φ ( d , q ) = φ ( d m , q ) , (3.2)where q are some internal variables for describing the irreversible processes and material states. Anothersimplification arises from the assumption that the cohesive law is isotropic in the crack plane. This indicatesthat the resistance to sliding is independent of the direction of sliding. and the effective displacement d m only depends on d n and the magnitude | d S | . Furthermore, the effective traction can be written as t m = ∂φ∂d m ( d m , q ) , (3.3)and it will be used to formulate the crack initiation criteria for activating the cohesive layers.The compressive normal stress should not cause the cohesive layer to fail, so that it does not contributeto the free energy density. In addition, no friction effect has been included in our cohesive model. Thetensile and compressive cases are considered separately as below,1. For the tensile case d n ≥
0, the effective opening displacement d m is d m = (cid:113) d n + β | d S | , (3.4)where the positive parameter β defines the ratio of effects from normal and shear displacements. Thecohesive law becomes t = ∂φ∂ d = ∂φ∂d m ∂d m ∂ d = t m d m (cid:0) d n n c + β d S (cid:1) . (3.5)Based on Eq. (3.1) and (3.4), the effective traction stress can be written as t m ( t , n c ) = (cid:113) t n + β − | t S | . (3.6)2. For the compressive case d n <
0, the effective opening displacement d m is d m = β | d S | . (3.7)The cohesive law becomes t = t n n + ∂φ∂ d S = t n n + ∂φ∂d m ∂d m ∂ d S = t n n + t m d m β d S . (3.8)Based on Eq. (3.1) and (3.8), the effective traction stress can be written as t m ( t , n c ) = β − | t S | . (3.9)As shown in Figure 6(a), a 1-D bilinear t m − d m effective cohesive law is adopted in this work. Tosimulate a perfectly bonded crack surface before the onset of failure, the initial stiffness of the elastic regime K should be set to a large value. The critical effective traction is denoted as t c , and the correspondingopening displacement as d c . The internal variables ( d , t ) store the feasible state with the maximumeffective traction. The cohesive layer fails completely at d m = d f with t m = 0. However, to avoid singularityin the analytical DMN calculation, we add a tiny stiffness κ to the entire t m − d m curve: t m ← t m + κd m . (3.10)For all the numerical examples to be studied in Section 4, we have K = 1 × GPa / mm , κ = 1 × − GPa / mm . (3.11)12 ! , 𝑡 ! 𝑑 " , 0 𝑑 ! 𝐾 ≫ 1 𝑑 , 𝜎 𝐺 ! = 12 𝑡 ! 𝑑 " 𝑡 ! (a) Effective cohesive law. 𝜎 ! 𝜏 ! 𝜎 ! 𝜎 " 𝜎 𝜎 ≤ 𝜎 " ≤ 𝜎 ! Potential crack initiation (b) Mohr’s circle in 3-D space.Figure 6: Cohesive layer for crack modeling. (a) The bilinear cohesive law for effective opening displacement d m and theeffective traction t m . The initial stiffness K is chosen to a large value, and the critical energy release rate G c is equal to thearea under the curve. An unloading-loading path during softening is highlighted as the red solid line. (b) Mohr’s circle fora 3-D stress state σ . The gray area covers the permissible normal and shear stress states, while the red solid line highlightspotential states for crack-plane initiation. Meanwhile, the normal component of traction t in Eq. (3.8) for d n < t n = Kd n , (3.12)which represents a nearly impenetrable crack surface. More details on the bilinear curve and the derivationof the tangent stiffness tensor for implicit analysis can be found in [57].Softening in the post-failure region is known to cause convergence difficulties in the implicit analysis.Meanwhile, when the total strain energy to release is higher than the fracture toughness of the homogenizedmaterial, the solution is hard to find during unstable crack propagation. To overcome these issues, viscousregularization is introduced in the DMN framework. In [57], we put the viscous effect directly on theopening displacement, which suits well in our debonding analysis with limited interfacial effects. However,for failure analysis interested in this paper, the separation displacement often dominates the overall straindue to localization, which causes nonphysical viscous stress even at the full failure of the material.Instead, we apply viscous regularization on the damage parameter of an inviscid backbone cohesivemodel, which is defined as D = σ σ c + K h ( d − d c ) , (3.13)where K h is a positive hardening stiffness of the reference undamaged material. Its magnitude is chosenbased on the base material’s stiffness and the cohesive layer’s characteristic length. Without viscosity, thishardening modulus does not affect the response of the cohesive layer. However, a properly selected K h willhelp improve the convergence of the implicit algorithm with viscous regularization. Specifically, we let theeffective stiffness of the cohesive layer ( K h /v c ) be equal to Young’s modulus of the base material E base , sothat K h = E base v c , (3.14)Here, v c is the reciprocal length parameter of the cohesive layer, and its expression will be provided in Eq.(3.30) when we discuss the activation of crack surfaces.The viscous damage parameter ˙ D v is controlled by the following evolution equation:˙ D v = 1 τ ( D − D v ) , (3.15)where τ is the viscosity coefficient representing the relaxation time of the viscous system. Here we use a13ackward Euler method to update D v . Finally, the effective traction stress after viscous regularization is t v = D v (cid:20) σ c + K h ( d − d c ) d d m (cid:21) , if d m < d or ˙ d m < D v [ σ c + K h ( d m − d c )] , otherwise. (3.16)In summary, the cohesive layers’ material behaviors are governed by four parameters: the critical effectivetraction t c , the critical energy release rate G c , the ratio β for defining the effective opening displacement,and the relaxation time τ for viscous regularization. Assume the stress state of the base material for an arbitrary micro-cell is σ . To decide whether a cracksurface should be activated, we will first search the plane(s) with the maximum effective traction, and itsnormal n c is n c = arg max n (cid:48) t m ( σ · n (cid:48) , n (cid:48) ) , (3.17)where the expressions of the effective traction t m are given in Eq. (3.6) and (3.9) for tension and compression,respectively.As shown in Figure 6(b), the Mohr circle of a three-dimensional stress state can be utilized to guide thesearching process. To arrive at the Mohr circle, we need to determine the values of principal stresses and theprincipal directions by the eigenvalue analysis of the stress tensor σ . Without loss of generality, we orderthe three eigenvalues as σ ≥ σ ≥ σ . (3.18)By the definition of t m , the optimum point should lie on the right half of the outer circle C , as highlightedby the red curve in the plot. For any point inside the permissible area, one can always find a point on C that has a larger effective traction. Searching on C also indicates a rotation along the second principal axis.The center and radius of C are ¯ σ = 12 ( σ + σ ) , ¯ τ = 12 ( σ − σ ) . (3.19)Here we use the angle θ to denote the rotation from the state with σ n = σ , so that the normal and shearstresses can be expressed as σ n = ¯ σ + ¯ τ cos 2 θ and τ n = ¯ τ sin 2 θ, (3.20)with θ ∈ [ − π/ , π/
4] corresponding to the right half of C Then, let us consider the effective traction in Eq. (3.6) for the tension loading case, or σ n > t m = (cid:112) (¯ σ + ¯ τ cos 2 θ ) + β − (¯ τ sin 2 θ ) . (3.21)The stationary points in terms of θ are θ ∗ = ±
12 arccos( ¯ σ ¯ τ ( β − −
1) ) with t (cid:48)(cid:48) m ( θ ∗ ) = (1 − β − )¯ τ sin θ ∗ . (3.22)Therefore, the solution θ ∗ exists as an potential global maximum point only if0 < ¯ σ ¯ τ ( β − − < σ > β < . (3.23)On the other hand, the effective traction in Eq. (3.9) for the compressive case, or σ n <
0, can be writtenas t m = β − | ¯ τ sin 2 θ | , (3.24)which reaches maximum when θ = ± π/
4. 14n summary, for any given σ and β , we only need to check the following candidates for the maximumeffective traction: θ = 0 , θ = ± π , (3.25) θ = ±
12 arccos( ¯ σ ¯ τ ( β − −
1) ) if 0 < ¯ σ ¯ τ ( β − − < σ > β < . Remark 4.
The failure algorithms are primarily designed for the specific cohesive law based on the effectivetraction-separation behavior in Section 3.1. Despite its simplicity, one limitation of the one-dimensionalmodel is that the fracture energy is constant regardless of the fracture mode. However, sometimes it becomesnecessary to consider the variances of fracture energy under different modes, e.g., Mode I and Mode IIfractures. In those situations, cohesive law models based on more general potentials can be used [31], whilethe viscous regularization scheme and the activations criteria may need to be revised accordingly.
When the maximum effective traction exceeds the critical value t c , t m ( σ · n c , n c ) > t c , (3.26)the crack surface with normal n c will be activated in a micro-cell. Depending on the solutions in Eq. (3.25),two crack surfaces might be activated at the same time. To differentiate the two crack surfaces numerically,a small perturbation is added to their t c values. Note that the displacement vector in the new cohesivelayers at loading step n is set to be d ( n − = σ ( n − · n c K , (3.27)where the superscript ( n −
1) labels a quantity from previous load step n −
1. This guarantees the equilibriumcondition after inserting the cohesive layers to the base material. Since the stiffness K is a large value (seeEq. (3.11)), the mismatch in the kinematic constraints does not affect the overall accuracy.Finally, we need to determine the reciprocal length parameter, as illustrated in Figure 2. Assume thescale tensor of the microcell is A jN . Similar to Eq. (2.15), the area of the crack surface is S ( A jN , n c ) = π (cid:113) det A jN | T R T n c | . (3.28)where T and R are defined in Eq. (2.8) and (2.9), respectively. Meanwhile, the volume of the micro-cellwith A jN is V = 4 π (cid:113) det A jN . (3.29)For an isotropic spherical micro-cell, the reciprocal length parameter v c can be naturally set equal to theinverse of the diameter of the sphere, with v c = 2 S/ V . Extending this formula to a general ellipsoidalmicro-cell yields v c ( A jN , n c ) = 2 S V = 12 | T R T n c | . (3.30)15 .3. Solution scheme for implicit DMN failure analysis Box 3.3.1 Solution scheme for implicit DMN analysis and crack activations
0. Initialization. M (0) = 0, L (0) = ∅ . Given A macro , perform cell division to get the scale tensor A jN for each micro-cell.1. For load step n , the macroscale strain increment is ∆ ε macro . Initialize the global crack configurationuse the converged one from the last load step n − M new ← M ( n − , L new ← L ( n −
2. Newton’s method:(a) If the number of iterations exceeds the limit, go to 5.(b) Evaluate the constitutive laws of the base material and enriched cohesive layers of each micro-cell to get the stiffness tensor and residual stress of the corresponding bottom-layer node(c) Forward propagation of stiffness tensors and residual stresses to the top node with C macro and δ σ macro (d) Compute the macroscale stress increment ∆ σ macro = C macro ∆ ε macro + δ σ macro (e) Backward propagation of stress and strain tensors from the top node to the bottom layer.Compute new incremental strain tensors of base materials and new incremental displacementvectors of cohesive layers(f) Check convergence. If not converged, go to 2(a)3. Crack activations:(a) Find the micro-cell index I c and surface normal n c , which satisfy the allowances (see Remark5) and maximize ∆ t = t m ( σ I c · n c , n c ) − t c (b) If ∆ t ≤
0, no more new crack surfaces, go to 4.(c) If ∆ t >
0, compute v c and initialize the internal variables q c of the new enriched cohesivelayer(s). Update the list M new ← M new + 1, L new ← L new + { I c , v c , n c , S } (d) Go to 2(a) to start a new round of iterations4. Global crack configuration reaches convergence. Update the internal variables, L ( n ) ← L new , M ( n ) ← M new
5. Invoke time step refinement, restore M ( n − , L ( n − , and all internal variables.To improve the overall convergence of the multiscale failure analysis, we propose a solution scheme whichdecouples the DMN implicit solver and the activation of new crack surfaces. First of all, let us define theglobal crack configuration L of DMN, L = {L , L , ..., L M } , (3.31)where M is the total number of crack surfaces. The list L i stores the information of the i -th crack surface(or cohesive layer in DMN formulation), L i = { I ic , v ic , n ic , S i } , (3.32)where I ic is the index of micro-cell where the crack surface is located, v ic is the reciprocal length parameter, n ic is the unit normal, and S i is the area of the crack surface. Note that once the cohesive layer for a crack16urface is activated, we also need to keep track of its internal variables q c in Newton’s method. At thebeginning of analysis, or load step 0, no crack surface exists in the network for an undamaged material, sowe have L (0) = ∅ . (3.33)In Box 4.3.1, we list the solution scheme for a material point (or integration point) in a concurrentmultiscale simulation, where the macro strain increment ∆ ε macro is applied at the top node of DMN.In all our simulations, the relative tolerances of the incremental strain tensors and the displacementvectors in Newton’s method (Step 2(f) in Box 4.1) are both set to 1 . × − for convergence, and themaximum iteration number is set to 40. Remark 5.
In Step 3(a), certain allowances can be specified to improve the scheme’s robustness for findingthe potential crack surfaces. In this paper, we introduce two allowances of crack activation in each micro-cell:(1)The total number of crack surfaces does not exceed 4; (2) The cosine similarity between the normals of anew crack surface and any existing one is less than √ / . The first allowance limits the memory cost of an enriched DMN, though adaptive memory allocationcould be implemented in the future to help relieve this constraint. The second allowance based on cosinesimilarity is introduced to suppress redundant crack activations. Due to viscous regularization, the stresssoftening is delayed, and the effective traction may exceed the critical value. As a result, new cracks maybe invoked at a similar orientation as the existing cracks, which is not desired for the analysis.If convergence cannot be reached for the current loading step, adaptive load step refinement is enabledto ease solving the nonlinear system. Essentially, the original load step is divided equally into two sub-steps. Note that the time increment should also be refined for the consistency of viscous regularization. Themaximum number of load step refinements is set to 10 in this paper. Detailed descriptions of the adaptiveload step refinement algorithm are provided in Box C.1.
4. Examples
In concurrent multiscale structural simulations, the macroscale FEA model can be either explicit orimplicit, since both the stress and the tangent stiffness tensor are available at the top node of DMN.However, in multiscale failure analyses with strain softening and localization, we will focus on the explicitdynamics in the macroscale FEA model, while the microscale DMN models are always solved implicitly.We begin with a single material point modeled by DMN in Section 4.1. The failure behaviors of theparticle-reinforced composite shown in Figure 4 will be evaluated under tension, compression, shear, andcyclic loadings. Meanwhile, we will study the effects of the macroscale length parameter h , the viscousrelaxation time τ , and the DMN depth N .Three representative examples of the concurrent multiscale simulation based on microscale DMNs willthen be presented in the following sections: 1) Section 4.2, the dynamic crush of a composite tube meshedby 3-D thin-shell elements; 2) Section 4.3, three-point bending test under 2-D plane strain condition; 3)Section 4.4, composite off-axis coupon test modeled by 3-D solid elements. Moreover, predictions from theoff-axis coupon model of a unidirectional carbon fiber reinforced polymer composite will be validated againstexperimental data from the literature [39]. In this section, the 3-D particle-reinforced composite is first trained by DMN and tested for its onlineextrapolation as a single material point. The particle phase is assumed to be linearly elastic, and the matrixphase is elasto-plastic with failure. All the material parameters are provided in Table 1. The default valueof the relaxation time τ is 1 . × − ms, while we will study its effect on the material responses in Section4.1.2. Additionally, the strain rate on the single material point is set to 1 . − .17 able 1: Material parameters of the particle reinforced composite with matrix failure. Particle E p (GPa) ν p E m (GPa) ν m σ Y (GPa) Hardening100.0 0.30 0.1 Eq. (4.1) t c (GPa) G c (GPa · mm) β τ (ms)0.15 6 × − . × − In terms of the matrix plasticity, von Mises plasticity with isotropic hardening surface is used, and theyield stress σ Y is described as a function of the effective plastic strain ε p , σ Y ( ε p ) = (cid:40) . · ε p , ε p ∈ [0 , . .
18 + 2 · ε p , ε p ∈ [0 . , ∞ ) GPa . (4.1)The macroscale material point is assumed to have the same length scale in all directions, so that thespherical macro-cell is represented by a isotropic scale tensor, A macro = 4 h I . (4.2)By default, the macro length scale h is set to 2 mm. The micro-cells after the cell division process basedon the trained network structure are shown in Figure 4 (c). Later in Section 4.1.1, we will also examine theeffects of h on the overall material responses.Figure 7 shows the stress-strain curves under six loading directions for the DMN with N = 8, whichhas 33 active DOF in the bottom layer. Because the particle-reinforced composite is nearly isotropic, it isexpected to behave similarly in three uniaxial tension loadings and in three shear loadings. For the uniaxialtension loadings, the composite reaches the maximum stress around 0.125 GPa, which is less than the criticaleffective traction t c of the matrix material due to stress concentration induced by the particle reinforcement.For the shear loadings, the plastic hardening effect is more prominent, and the softening starts at a lowerstress level of 0.10 GPa. In addition, treemap plots colored by the released energy per crack surface area G in each DOF at ε macro = 0 .
03 are presented above the stress-strain curves. A micro-cell with fully separatedcrack surfaces has G equal to G c ( G c = 0 . · mm, see Table 1), while G vanishes if there exists nocrack surface. Essentially, this contour treemap plot reflects the global crack configuration of DMN, whichis shown to vary with the applied macroscale boundary conditions.To demonstrate the physics introduced by the cohesive layers, we compare the material responses underuniaxial tension and compression loadings in Figure 8 (a). As the compressive stress does not contribute tothe effective traction t m of a potential crack surface in the matrix phase, the composite material reaches ahigher stress magnitude of 0 .
25 GPa under compression. For all the cases, the DMNs fail completely withnegligible residual stresses.
Remark 6.
Although not necessary for predicting the multiscale failure responses, it is useful to define aglobal damage indicator for tracking the failure process and informing discontinuous crack algorithms in themacroscale model. However, in the concurrent multiscale simulations to be shown in Section 4.3 and 4.4, weget around potential numerical issues like element distortions by using a volumetric strain-based criterion totrigger the element deletion in the macroscale model.
The cyclic loading behaviors with different strain magnitude ε max are shown in Figure 8 (b). When theDMN material is unloaded at ε = 0 .
005 in the middle of the softening process, it first goes through alinear elastic regime with a degraded stiffness. It then transits to a compressive state with almost the samestiffness as the undamaged material. The reloading curve closely follows the unloading curve but differs18 . . . . . . . G (GPa · mm) . . . . . . . G (GPa · mm) 𝜎 !! − 𝜀 !! 𝜎 "" − 𝜀 "" 𝜎 − 𝜀 . . . . . . . G (GPa · mm) .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . Normal strain . . . . . . . S t r e ss ( G P a ) æ ° " æ ° " æ ° " . . . . . . . G ( M P a · mm ) . . . . . . . G ( M P a · mm ) (a) Tension loadings. . . . . . . . G (GPa · mm) . . . . . . . G (GPa · mm) . . . . . . . G (GPa · mm) 𝜎 !" − 𝛾 !" 𝜎 " − 𝛾 " 𝜎 ! − 𝛾 ! .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . Shear strain . . . . . . . S t r e ss ( G P a ) æ ° ∞ æ ° ∞ æ ° ∞ . . . . . . . G ( M P a · mm ) . . . . . . . G ( M P a · mm ) (b) Shear loadings.Figure 7: DMN predictions for the particle-reinforced composite under three uniaxial tension loadings and three shear loadings.The number of layers N is 8, and the number of microscale DOF or micro-cells in the DMN is N dof = 33. For each loadingcase, the contour treemap plot of the released energy per crack area G at ε macro = 0 .
03 is placed above the stress-strain plot.In a micro-cell, G vanishes if no crack surface is activated. ° . ° . ° .
02 0 .
00 0 . Normal strain ° . ° . . . S t r e ss ( G P a ) æ ° " æ ° " æ ° " TensionCompression (a) Monotonic loadings in tension and compression. . . . . . . . " ° . ° . . . . æ ( G P a ) " max = 0 . " max = 0 . 𝑡 . . . . . . . " ° . ° . . . . æ ( G P a ) " max = 0 . " max = 0 . . . . . . . . " ° . ° . . . . æ ( G P a ) " max = 0 . " max = 0 . (b) Cyclic loadings with different strain magnitudes.Figure 8: Stress-strain curves of the particle reinforced composite under more general loading conditions. The number of layers N is 8. In (a), the uniaxial tension and compression loadings are evaluated independently. In (b), each uniaxial cyclic loadingpath contains two loading and one unloading processes between 0 and ε max , as shown in the embedded strain history plot. slightly due to plasticity induced at the end of compression. For the second case, the material is close tocomplete failure at ε = 0 . .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . ε . . . . . . . σ ( G P a ) h = 0 . mm h = 0 . mm h = 2 . mm h = 10 . mm h = 50 . mm (a) Uniaxial tension σ − ε . ° ° Macro length scale h (mm) . . . . . M a x i m u m s t r e ss ( G P a ) ° ° ° ° ° R e l ea s edene r g y ¶ ( k N · mm ) 𝑡𝑑 ! = 0.008 mm (b) Maximum stress and released energy.Figure 9: Effects of the macro length scale h . The macro scale tensor is A macro = 4 /h · I , and the number of layers N is 8.In (b), the maximum stress is extracted from the corresponding σ − ε curve in (a). The released energy Π is the total freeenergy φ of all the activated crack surfaces in the network. The dashed line highlights the full separation displacement of thecohesive layer d f = 0 .
008 mm. the scale transition of DMN enabled by the proposed cell division process overcomes some critical issuesof RVE-based models with strain localization, achieving a consistent and efficient framework for multiscalefailure analysis.
As discussed in Section 1.2, the macro length scale h is usually predetermined by the macroscale elementsize or the nonlocal regularization size. The cell division process consistently tracks the scale transition inDMN without retraining the network or introducing extra energy regularization for different macro lengthscales.Figure 9 (a) presents σ − ε curves of h varying from 0.08 mm to 50.0 mm. As we can see fromthe plots, the magnitude of the macroscale softening stiffness increases with h . When h is larger than 10.0mm, the stress-strain curve has a sharp drop and does not change much as h increases. In these cases, theactivated crack surfaces’ overall fracture energy is not sufficient to release the total strain energy in the bulkmaterial, and the excess energy is mainly dissipated via viscosity.Figure 9 (b) plots the maximum stress and released energy against the macro length scale h underuniaxial tension loadings. The maximum stress is directly extracted from the stress-strain curves in Figure9 (a), while the released energy Π is computed as the total free energy φ (see Eq. (3.2)) of all the cracksurfaces in L = {L , L , ..., L M } : Π = M (cid:88) i =1 φ i S i , (4.3)where S i is the surface area of the i -th crack as defined in Eq. (2.15). Note that this released energy isnot equal to the area under the stress-strain curve, due to the existence of matrix plasticity and viscousdissipation. The energy predictions for h < .
08 mm are not shown as the material point does not failcompletely at ε = 0 . h approaches the characteristic length of the cohesive layer d f = 0 .
008 mm, the material point’smaximum stress (or strength) reaches a plateau close to the critical effective traction of the cohesive layer t c = 0 .
15 GPa. When the material point has a small macro length scale, the reciprocal length parameters v c of the cohesive layers become so large that the global softening is delayed. Eventually, the overall strengthis governed by the cohesive layers in the matrix phase.20 .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . ε . . . . . . . . σ ( G P a ) τ = 1 × − ms τ = 1 × − ms τ = 1 × − ms τ = 1 × − ms τ = 1 × − ms (a) Relaxation time τ . .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . Strain . . . . . . . S t r e ss ( G P a ) Normal Shear N = 4 , N dof = 4 N = 6 , N dof = 13 N = 8 , N dof = 28 (b) Number of DMN layers N .Figure 10: Effects of model parameters to be selected by modelers. In (b), each stress-strain curve is averaged from all threeuniaxial tension (or shear) loadings. Another important observation is that the released energy in the material point is approximately pro-portional to h . This is physically sound since the material point has an isotropic scale tensor and theglobal (effective) crack surface is a two-dimensional manifold cutting through the material point. If oneuses volumetric damage mechanics in the matrix phase, the resulting energy would be proportional to h ,which requires post-regularizations based on h (e.g., element size) to restore the energy consistency [61, 38].However, in the DMN framework, the property Π ∝ h is reproduced naturally from the cell division processand the cohesive-layer enrichment. The purpose of putting the viscous regularization on the cohesive layer’s damage parameter is to overcomethe convergence difficulties in an implicit analysis. A larger relaxation time τ puts more viscosity or dampingin the system. On the one hand, we want to keep the added viscosity small enough to limit its effects onthe material responses. On the other hand, a more damped system tends to be more stable and requiresfewer time-step refinements to converge. Therefore, τ should be chosen properly to balance accuracy andefficiency.Figure 10 (a) shows σ − ε curves for τ ranging from 1 . × − to 1 . × − ms. For the prescribedstrain rate ˙ ε = 1 . − , we observe little viscous effect when τ is less than 1 . × − ms. However,smaller τ sacrifices the efficiency due to more time-step refinements in the implicit analysis. As τ increasesto 1 . × − ms, the onset of softening along the hardening curve is greatly delayed, and more strain energyis dissipated through the artificial viscosity. Although the macroscale strain rate we applied in this studyseems very high at first glance, it can often appear in the strain localization region, especially when themacro length scale h (e.g., element size) is tiny.Moreover, the number of DMN layers N , also referred to as the “network depth”, is an important hyper-parameter that controls the network’s complexity. Figure 10 (b) provides the stress-strain curves undertension and shear loadings for N = 4, 6, and 8. Since the network with N = 4 only has 4 DOF (seeFigure 4), we observe that its responses differ a lot in three orthogonal normal or shear loadings, while thenetwork with N = 8 shows much more isotropic failure behaviors as we can see from Figure 7. To betterdemonstrate the results, we only plot the average stress-strain curves. Overall, the maximum stress andsoftening stiffness predictions are consistent across different network depths, while deeper networks withmore DOF yield smoother responses.DNS results of RVE models are available for the nonlinear elastoplastic composite without matrix failure,so we will use them to validate the accuracy of DMN. Figure 11 shows the stress-strain curves from DMNand DNS under cyclic uniaxial tension and shear loadings. Good convergence of DMN to DNS results is21 .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . ε − . − . . . . σ ( G P a ) N = 4 , N dof = 4 N = 6 , N dof = 13 N = 8 , N dof = 28 DNS (a) Uniaxial tension σ − ε . .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . γ − . − . − . − . . . . . σ ( G P a ) N = 4 , N dof = 4 N = 6 , N dof = 13 N = 8 , N dof = 28 DNS (b) Shear σ − γ .Figure 11: Results of the particle-reinforced composite without matrix failure. The network depths are N = 4, 6, and 8. TheDNS results are denoted by circles. observed. For N = 6 and N = 8, the networks can predict the hardening behaviors very well for bothloading cases, although only linear elastic data are used for offline training. In this section, we apply the microscale DMNs of the particle reinforced composite to concurrent mul-tiscale simulations. In the macroscale, a symmetric crush tube is impacted by a moving wall on the topsurface. The tube dimensions are provided in Figure 12 (a), and the thickness of the tube is 2 mm. Due tothe model symmetry, only a quarter of the tube is simulated. The velocity of the rigid wall is v = 4 mm/ms.The macroscale crush tube is meshed by the Belytschko-Tsay shell elements with thickness stretch, and3 integration points are used across the thickness direction for each element. For the mesh refinement study,three element sizes are investigated: 8.32 mm, 4.16 mm, and 2.08 mm. Accordingly, the numbers of shellelements are 467, 1866, and 7464, respectively. As the tube is mainly under compressive loads, no elementdeletion will be triggered to sustain physical self-contacts in the model.Each integration/material point in the crush tube is coupled to a DMN with N = 8 and N dof = 28. Theoverall density of the DMN material is ρ = 5 × − kg / mm . Other material parameters of the particleand matrix phases can be found in Table 1. As no non-local or gradient-based regularization is used, we canregard the element size as the in-plane macro length scale. Moreover, the material axis 3 shown in Figure4 (a) is always normal to the shell plane, and the out-of-plane macro length scale is defined by the shellthickness. Given the thickness equal to 2 mm, the macro scale tensor A macro of each integration point is A macro = /h /h . mm − . (4.4)Figure 12 (b-d) show contour plots of the released energy on the deformed tube at T = 3, 9, and 15ms, The definition of the released energy in DMN is given in Eq. (4.3), while the element-wise value isaveraged over all three integration points. The mesh size of the tube model is 4.16 mm. As we can seefrom the snapshots, a hexagonal crack pattern is formed under the impact, and the failures mostly appearin the middle section. We also simulate an identical tube model except no matrix failure is considered inthe composite. The snapshots of the deformed tube at T = 3, 9, and 15 ms are provided in Figure 12 (e-g).As no crack surface is activated, the released energy is always zero in these plots. After the initial buckling,the plastic deformations are concentrated at the top section of the tube.Furthermore, we perform the mesh refinement study on the crush tube. First, the contour plots of thereleased energy are shown in Figure 13 (a-c) for mesh sizes 8.32 mm, 4.16 mm, and 2.08 mm, respectively.22 . mm m m mm (b) (c) (d)(e) (f) (g)(a) Π kN $ mmΠ kN $ mm
𝐹, 𝑣 = 4 mm/ms
Figure 12: The crush tube modeled by the particle-reinforced composite. (a) The geometry of the tube. The rigid wall ismoving downwards at 4 mm/ms. Only a quarter of the tube is simulated due to symmetry. (b-d) Snapshots of the tube withmatrix failure at T = 3, 9, 15 ms. (e-f) Snapshots of the tube without matrix failure at T = 3, 9, 15 ms. The mesh size of thetube is 4.16 mm. The deformed plots are colored by the released energy Π of each element. Π kN $ mm Π kN $ mm Π kN $ mm̅𝜀 ! ̅𝜀 ! ̅𝜀 ! (a) (b) (c)(d) (e) (f) Figure 13: Snapshots at T = 3 ms of failed tubes with different macroscale mesh sizes. (a-c) Contour plots of the releasedenergy Π for mesh sizes 8.32 mm, 4.16 mm, and 2.08 mm. (d-f) Contour plots of the average effective plastic strain at theinner surface for mesh sizes 8.32 mm, 4.16 mm, and 2.08 mm. The crack patterns predicted by the three models are similar in terms of the overall shape but differ slightly insome local regions. One may also observe that the magnitude of the released energy is nearly proportionalto the mesh size. In contrast, the released energy of the single material point shown in Figure 9 (b) is23
Time (ms) − − − − R i g i d - w a ll f o r c e ( k N ) Coarse, h = 8 . mmMedium, h = 4 . mmFine, h = 2 . mm (a) Composite tube with matrix failure. Time (ms) − − − − R i g i d - w a ll f o r c e ( k N ) Coarse, h = 8 . mmMedium, h = 4 . mmFine, h = 2 . mm (b) Composite tube without matrix failure.Figure 14: Histories of the rigid-wall force for mesh sizes 8.32 mm (Coarse), 4.16 mm (Medium), and 2.08 mm (Fine). proportional to the square of the macro length scale h . This is because the DMN macro-cell of the shellelement is anisotropic, with the macro length scales changes only in the shell plane.The average effective plastic strain ¯ ε p is another physical quantity to characterize the overall materialstate. It is defined as ¯ ε p = N dof (cid:88) i =1 ε ip f i , (4.5)where ε ip is the effective plastic strain of the i -th DOF, and f i is its volume fraction in the network. Typically, ε p vanishes if the phase is linear elastic. The contour plots of ¯ ε p on the inner surface for different elementsizes are provided in Figure 13 (d-f). More than the plastic deformation patterns, the magnitudes of theaverage effective plastic strain are around the same, indicating that the onset of failure is not sensitive tothe mesh size.Finally, we summarize the histories of rigid-wall force with and without matrix failure in Figure 14. Thefirst two spikes appear similarly in the two cases, which results from the initial elastic buckling of the tube.Afterward, the force on the composite tube with matrix failure drops to zero, while the one without matrixfailure stables at a finite value around -11 kN. Regarding the convergence on the mesh size, the medium andfine meshes predict very close results for the first case. While for the second case with only matrix plasticity,all three meshes output similar force history curves. From this section, we switch the focus to the unidirectional-fiber microstructure shown in Figure 5 (a)for modeling carbon fiber reinforced polymer composites. The volume fraction of the fiber phase is 50%.The number of layers is N = 8. After the offline training, 31 active DOF remains in the network, and 8of them belong to the fiber phase. The carbon fibers are modeled as an orthotropic linear elastic materialwith large stiffness in the fiber direction, and no fiber failure is considered in this paper. The epoxy matrixis modeled by isotropic von Mises plasticity with an exponential hardening law, σ Y ( ε p ) = ( σ y − σ u ) exp( − aε p ) + E h ε p + σ u , (4.6)where σ y represents the yielding strength, σ u is the ultimate yield stress for large effective plastic strain, E h is a linear hardening stiffness, and a is a dimensionless constant. The carbon fibers and epoxy are assumedto be perfectly bonded. The overall density of the UD composite is set to ρ = 1 . × − kg / mm . (4.7)24 .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . ε . . . . . . . . σ ( G P a ) h = 0 . mm h = 2 . mm h = 10 . mmExp. (a) Transverse tension with different h . .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . Strain . . . . . . . . S t r e ss ( G P a ) Transverse tensionLongitudinal tensionTransverse shearLongitudinal shearExp. (b) Various loading directions.Figure 15: Single material point tests on the 3-D unidirectional carbon fiber reinforced polymer composite. The macro scaletensor is A macro = 4 /h · I . In (b), the macro length scale h is 2.0 mm. Experimental results are plotted as circles ( ◦ ). All the microscale material parameters are provided in Table 2. Note that the elastic constants are measureddirectly from experiments on single-phase materials. The hardening parameters ( σ y , σ u , E h , a ) are fittedbased on the uniaxial tension and shear tests on the epoxy matrix. The critical energy release rate is set to3 × − GPa · mm. Specifically, the critical effective traction t c of the matrix phase is determined inverselyby matching the transverse tension curve of the UD composite. Table 2: Material parameters in the unidirectional-fiber microstructure for the CFRP composite. No fiber failure is consideredin this work.
Carbon fiber E (GPa) E (GPa) E (GPa) G (GPa) G (GPa) G (GPa)245.0 19.8 19.8 29.2 29.2 5.9 ν ν ν E m (GPa) ν m σ y (GPa) σ u (GPa) E h (GPa) a t c (GPa) G c (GPa · mm) β τ (ms)0 .
10 3 × − × − We first evaluate the responses of a single material point. Similarly, its macro-scale tensor is assumedto be A macro = (4 /h ) I , representing an isotropic sphere in space with diameter equal to 2 h . Figure 15 (a)shows the stress-strain curves under transverse (in-plane) tension for h ranging from 0.4 mm to 10.0 mm.As expected, the magnitude of the softening stiffness increases with h , while the strength predictions are allat the same level, around 0.065 GPa.The anisotropy of material responses induced by the microstructure is demonstrated in Figure 15 (b).Experimental data for longitudinal and transverse tension loadings are obtained from 0 ◦ and 90 ◦ coupontests, respectively. Without any parameter fitting, the DMN predictions of the elastic stiffness match theexperimental results very well. By adjusting the critical effective traction t c , the UD composite’s failurebehaviors under transverse tension can also be properly captured. In Section 4.4, the calibrated t c will befurther validated by the 10 ◦ coupon experiment. Another interesting finding is that the composite encountersnegligible failure under the longitudinal shear loading for the strain less than 0.05, because its deformationis shear-dominated and the matrix phase has less stress concentration comparing to other loading directions.We then apply the DMN to concurrent multiscale simulations for three-point bending tests of the UDcomposite. The geometry and boundary conditions are illustrated in Figure 16. No crack zone is predefined25 , 𝑣 = 0.2 mm/ms
300 mm400 mm mm ⨀ Figure 16: Illustration of the macroscale three-point bending model under 2-D plane-strain condition. Two cases with 1)out-of-plane and 2) in-plane 45 ◦ fiber orientations are considered. in the macroscale model so that damage initiations are handled inside the DMN. Although the microscaleDMN model is still 3-dimensional, 2-D plane-strain conditions are considered in the macroscale model. Theloading head moves downwards at a speed of v = 0 . ◦ angle to the horizontal axis. Mesh refinement study is performed for each case. The coarse, medium andfine mesh sizes are 5, 2.5, and 1.25 mm, respectively. The corresponding numbers of 2-D plane-strain finiteelements are 1280, 5120, and 20480, and each element has four integration points. No nonlocal regularizationis applied here, so the macroscale softening zone due to matrix failure will localize in one layer of elements.Therefore, we set the in-plane macro length scale h equal to the mesh size. Given the section thickness of 2mm, the macro scale tensor of 3-D DMN can be written as A macro = /h /h . mm − . (4.8) Remark 7.
Under the 2-D plane strain condition, the results should be independent of the section thickness.It requires that the crack surfaces activated in DMN are vertical to the 2-D plane in the global coordinatesystem, which has not been considered in this work. Under this extra constraint, the failure algorithms willneed to be modified as the crack surfaces are not necessarily activated in the direction with maximum effectivetraction of a 3-D stress state.
To avoid excessive element distortions, an element will be deleted if at least two of its integration pointshave the volumetric strain larger than 0.1, ε + ε + ε > . . (4.9)According to our study, this is sufficient to guarantee the DMN with h ≥ .
25 mm to fail completely or havenegligible residual stress.Figure 17 summarizes the snapshots of the three-point bending tests for various fiber orientations andmesh sizes, rendered by the average effective plastic strain ¯ ε p (see Eq. (4.5)). For all the samples, the26 = 12.8 ms 𝑇 = 13.5 ms Coarse: ℎ = 5 mm Medium: ℎ = 2.5 mm Fine: ℎ = 1.25 mm 0.004 ̅𝜀 ! (a) Case 1, out-of-plane fiber orientation. 𝑇 = 12.4 ms 𝑇 = 13.1 ms 0.0080.0000.004Coarse: ℎ = 5 mm Medium: ℎ = 2.5 mm Fine: ℎ = 1.25 mm ̅𝜀 ! (b) Case 2, in-plane 45 ◦ fiber orientation.Figure 17: Contour plots of the average effective plastic strain ¯ ε p for different fiber orientations and mesh sizes. Elementdeletions are enabled in both cases. macroscopic crack initializes from the bottom side, despite that the average effective plastic strain is largerat the loading area on the top side. This is physically sound as the composite is stronger under compressivestress states, which is captured naturally by the multiscale DMN with cohesive layers.In Case 1, the cracks propagate vertically for all the mesh sizes. At T = 13 . h = 5mm. In Case 2, the cracks first propagate in the vertical direction and then turn towards the fiber directionat 45 ◦ . Different from Case 1, where the cracks are mainly triggered by transverse tension loading, Case2 has a more complex stress state at the crack initiation point, with a combination of transverse tensionand longitudinal shear loadings. One can also see from the deformed coarse mesh in Figure 18 (b) thatshear deformations are prominent in the localization band. We did not consider the fiber failure and finitedeformations in the DMNs, which could influence the Case 2 results. Overall, the predicted crack paths areconsistent under mesh refinement.Figure 18 shows the histories of applied forces for different fiber orientations and mesh sizes. Due tothe strengthening effect from the carbon fibers, the models in Case 2 are stiffer than the ones in Case 1.Slight plastic yielding can be observed in Case 2, indicating more shear deformations in the models. Goodconvergence of the load curves concerning the mesh size is achieved in either case.27 . . . . . . . Time (ms) − . − . − . − . − . − . − . . Load ( k N ) Coarse, h = 5 mmMedium, h = 2 . mmFine, h = 1 . mm (a) Case 1, out-of-plane fiber direction. . . . . . . . Time (ms) − . − . − . − . − . − . − . . Load ( k N ) Coarse, h = 5 mmMedium, h = 2 . mmFine, h = 1 . mm (b) Case 2, in-plane 45 ◦ fiber direction.Figure 18: Histories of the applied forces in the three-point bending tests for mesh sizes 5 mm (Coarse), 2.5 mm (Medium),and 1.25 mm (Fine).
120 mm 45 mm45 mm . mm °𝑧 𝑥𝑥 𝑦 Figure 19: Illustration of 10 ◦ off-axis tensile coupon. Axes are shown on the right of the plots. The red dashed arrow denotesthe fiber direction of the composite. Only half of the model is simulated due to the symmetry of the middle x − y plane. In the last example, we apply the DMN of the unidirectional CFRP composite to a 10 ◦ off-axis tensilecoupon test with experimental validations [39]. All the material parameters remain the same as in Section4.3. The geometry of the tensile coupon is shown in Figure 19. The length, width, and thickness of thecoupon’s middle measuring section are 120 mm, 12.7 mm, and 2.42 mm, respectively. In the experiment,tab sections at the two ends are tightly clamped so that no movements in the x and z directions are allowed.The tensile loading velocity in the y direction is 0.0167 mm/s. However, to reduce the simulation time, weincrease this velocity to 0.0167 mm/ms, considering that the material plasticity and crack activation modelsare not rate-dependent, and it is still much less than the wave speed. Meanwhile, mass scaling is introducedso that the critical time step is no less than 1 . × − ms.Only the middle section is modeled by DMN, while the tab sections are described by an elastic materialwith the Young’s modulus equal to 8 GPa. The middle section is meshed uniformly by linear 8-node solidelements with sizes ∆ x = 0 .
808 mm, ∆ y = 0 .
808 mm, and ∆ z = 0 .
606 mm. The coupon is symmetricabout the central x − y plane, and the total number of elements in the symmetric model is 8320. Each solidelement has one integration point, which is linked to a microscale DMN. Based on the mesh sizes, the macroscale tensor of DMN is defined as A macro = / ∆ x / ∆ y / ∆ z = .
127 6 .
127 10 . mm − . (4.10)Like the three-point bending tests, element deletion will be triggered if the volumetric strain at an integrationpoint is larger than 0.1. 28 . . . . . . . Normal strain N o r m a l s t r e ss ( M P a ) ◦ model ◦ model (*)Exp. (a) Normal stress vs. normal strain. 𝑇: 60.00 → 70.00 ms;
Δ𝑇 = 0.25 ms; ̅𝜀 ! ∈ 0.0, 0.035 (b) Snapshots of crack propagation. Experiment:Multiscale simulation: (c) Coupon crack formations in experiment and multiscale simulation.Figure 20: Results of the 10 ◦ off-axis tensile coupon. In (a), experimental data of the 10 ◦ coupon is plotted as the circles. (*)Results of a 9 ◦ off-axis model is also presented to show the sensitivity of stress-strain curve to fiber orientation. Snapshots in(b) are colored by the average effective plastic strain ¯ ε p . The normal stress vs. normal strain curve predicted by the concurrent multiscale simulation is comparedwith the experimental results in Figure 20 (a). A good match of the curves before failure is observed.The multiscale simulation underestimates the maximum strain and stress, potentially due to higher stressconcentrations at the clamping areas in the numerical model. To demonstrate the effects of fiber orientation,we also simulated a 9 ◦ off-axis model. Although the angles are only differed by 1 ◦ , the stress-strain curveschange notably due to the strong anisotropy of the UD composite.Moreover, Figure 20 (b) shows the evolution of crack formations in the macroscale model, and Figure 20(c) compares the crack formations in the simulation and experiment. As we can see from the snapshots, thecrack starts from the side close to the gripping area and propagates along the 10 ◦ fiber direction, consistentwith the experiment. One discrepancy is that the crack initiates further away from the clamping area in theexperiment, which may explain the difference of failure-strain predictions in Figure 20 (a). Nevertheless,with basic microscale material models (e.g., elasticity, isotropic plasticity, 1-D effective cohesive law) andminimal material calibrations, the DMN-enabled multiscale predictions agree well with the experimentalvalidation data.
5. Advantages and limitations
At this point, we have demonstrated DMN’s capabilities of modeling multiscale systems with materialfailure. In terms of strain localization modeling, its advantages over traditional RVE-based methods are two-fold: 1) The macro length scales are naturally propagated to the microscale DOF so that the constraints on29 able 3: Wall times of various concurrent multiscale simulations on 10 CPUs. The simulated duration and the initial time stepare provided for each example. Note that we have applied mass scaling in the 10 ◦ off-axis tensile coupon test. Crush tube (4.2) Crush tube (4.2) 3-point bend (4.3) ◦ coupon (4.4) Damaged, 15 ms No damage, 15 ms Case 1, 16 ms Mass scaling, 100 msTimestep(ms) Cost(hour) Timestep(ms) Cost(hour) Timestep(ms) Cost(hour) Timestep(ms) Cost(hour)Coarse 2 . e − . e − . e − . e − . e − . e − . e − . e − . e − . e − N = 8and N dof = 28, implemented in FORTRAN, only took 0.1 s for the same amount of loading steps on oneCPU. Moreover, thanks to the nature of the hierarchical binary-tree structure, the computational cost ofDMN is proportional to the number of DOF in the network [1, 2].The wall times of various examples in Section 4 are listed in Table 3, including the crush tubes withand without matrix failure, the three-point bending test with out-of-plane fibers, and the 10 ◦ off-axis tensilecoupon test. For the same shell geometry and material properties, the critical time step in the macroscalesimulation is inversely proportional to the mesh size, while the number of elements is inversely proportionalto the square of the mesh size. If the mesh size is refined by half (e.g., Coarse → Medium), the computationaltime becomes approximately 8 times longer. Moreover, by comparing the times of crash tubes with andwithout damage (see the snapshots in Figure 12), we can conclude that the enrichment of cohesive layers formodeling new crack surfaces increases the computational cost by around 10% in our current implementation.The computational cost for the damaged crush tube with fine meshes is not provided in the table, asthe simulation was terminated around 5 ms due to the convergence difficulties. After the tube was buckledunder compression, some failed elements were highly distorted, and no convergence of the DMN implicitsolver was achieved even after 10 time-step refinements. For all the other cases, the simulations were finishedsuccessfully. Note that convergence issues due to severe element distortions could be avoided by the elementdeletion introduced in the three-point bending tests and the 10 ◦ off-axis tensile coupon test.Although not discussed in this paper, the network interpolation from transfer learning is also an appealingfeature of DMN. In [53, 54], databases of multiple microstructures are unified to cover the design space ofcomposites. The learned database can then be applied to concurrent multiscale simulations with localmicrostructure variations resulting from the manufacturing processes, such as the injection molding of shortfiber-reinforced composites and the metallic additive manufacturing.Several directions can be considered to further enhance the DMN framework in terms of both theoriesand implementations:1. The finite-strain formulation of DMN should be considered for more general problems. In Appendix B,we discuss several essential aspects of the finite-strain formulation, including its benefits for capturinglocal crack orientations, new algorithms for the crack initiation, and key steps of deriving the analyticalfunctions in the cohesive building bock.2. Nonlocal regularization methods can be used to reduce the mesh sensitivities in the macroscale model.In this work, the macro scale tensor of DMN is defined by the side lengths of a rectangular or cuboid30lement. However, for irregularly shaped elements, we need a more general method to relate the macroscale tensor to the mesh characteristics. Nonlocal regularization avoids this issue as the macro lengthscale is set by the nonlocal size.3. Element deletions cause losses of mass and momentum in the system. In this regard, more advancednumerical techniques should be used for crack modeling in the macroscale model, such as the extendedfinite element method [58, 62], and meshfree methods with bond-breaking mechanisms [63, 64]. Otherthan the volumetric strain adopted in this work, it is also helpful to define a handy global damageindicator for triggering these crack separation schemes.4. To further speed up the concurrent multiscale simulation, one can distribute the macroscale materialpoints to more CPUs. As the computational costs of a deforming DMN depend on the level ofmaterial nonlinearities and the number of enriching cohesive layers, the dynamic load balancing shouldbe equipped for better parallelization efficiency. GPU computing is also a promising direction toinvestigate.
6. Conclusions and future work
In this paper, we propose a new cell division scheme for consistent scale transitions within the DMNframework. In a multiscale failure analysis with strain localization effects, it overcomes the difficulties ofchoosing proper RVE sizes and applying boundary conditions. For a two-layer DMN building block, wederive the mathematical formulations of the cell division process, which enables the backward propagationof length scales from the macroscale material point to the microscale DOF. Algorithms for the activationof crack surfaces and the implicit failure analysis are proposed based on the cohesive-layer enrichment ofDMN. The cohesive layers are modeled by an effective traction-separation law with viscous regularizationon the damage parameter.We investigate two microstructures: the particle-reinforced composite and the unidirectional fiber com-posite. In particular, the network of the 3-D unidirectional-fiber composite is directly transferred from its2-D cross-section model based on the physical interpretations of fitting parameters. Parametric studies on asingle material point are performed to evaluate the effects of the macro length scale, the relaxation time, andthe number of network layers under various loading conditions. We provide three examples of concurrentmultiscale simulations with different element formulations in the macroscale, including 3-D thin shell, 2-Dplane-strain, and 3-D solid elements. Simulation results of the 10 ◦ off-axis CFRP tensile coupon test arefurther validated against the experimental data of the stress-strain responses and the crack formation.We believe the proposed cell division scheme for scale transition sets a strong basis of DMN for multiscalefailure analysis and general problems with strain localization across scales. Together with its accuracy,efficiency, and extrapolation capabilities, the DMN framework provides a feasible way of pushing machinelearning and data-driven multiscale materials modeling to computer-aided engineering at an industrial scale.It admits many straightforward extensions:1. The DMN framework for multiscale modeling and concurrent simulations can be applied to other typesof material systems, such as nano-particle reinforced rubber composites, polycrystalline materials, andvarious carbon fiber reinforced polymer composites.2. More complex material behaviors can be considered. For instance, in the unidirectional fiber compos-ite, one may model the fiber failure, matrix failure, and interfacial debonding simultaneously. Moregeneral cohesive laws could be adopted [31], while algorithms for the crack activation and tangentstiffness computation need to be modified accordingly. Multiphysics materials behaviors could also beconsidered in the framework.3. In the aspect of Integrated Computational Materials Engineering (ICME), it is important to derivethe process-structure-property relationship for a multiscale material system. In [54], we have demon-strated an integration of the manufacturing process simulation and microstructure-sensitive material31odels for short fiber reinforced composites based on DMN and transfer learning. With the new scaletransition scheme, this approach can be applied to a broader class of materials.4. Efficiency and physical interpretability of DMN make it suitable for materials design, optimization,and uncertainty quantification. Trained by linear elastic data, the DMN can be extrapolated topredict nonlinear material responses. Therefore, it reduces the sample-size requirement of offline datageneration, which tends to be time-consuming for deep learning models without embedded physics. Acknowledgments
Z. Liu would like to acknowledge Dr. Haoyan Wei, Dr. C.T. Wu, and Dr. Yong Guo for the helpfuldiscussions. Z. Liu would like to thank Dr. Cheng Yu for providing the adaptive time step refinementalgorithm and Dr. Jiaying Gao for sharing the coupon model and the NIST experimental data. Finally, Z.Liu would like to thank Xinying Yu for stimulating creativity.
Appendix A. Machine Learning of DMN: Data generation and optimization
Here, we summarize the key steps of training a DMN. For more detailed information, interested readersare referred to our previous papers for 2-D materials [1] and 3-D materials [2]. First, for the training basedon linear elastic data, the output function of a two-phase material in 3-D space can be written as¯ C rve (cid:124) (cid:123)(cid:122) (cid:125) Output = f ( C p , C p (cid:124) (cid:123)(cid:122) (cid:125) Inputs , Fitting parameters (cid:122) (cid:125)(cid:124) (cid:123) z j =1 , ,..., N − , α k =1 , ,..., i − i =1 ,...,N , β k =1 , ,..., i − i =1 , ,...,N , γ k =1 , ,..., i − i =1 , ,...,N ) , (A.1)where ¯ C rve is the composite material’s overall stiffness tensor predicted by the network, C p and C p arethe elastic stiffness tensors of microscale phases. In the design of experiments (or sampling) for C p and C p , we consider both phases to be orthotropically elastic.For a 3-D microstructure like the particle reinforced composite tested in this paper, we usually generate500 samples from finite element analyses on the same RVE geometry but with different material inputs.The first 400 samples are selected as the training dataset, and the remaining 100 samples go into the testdataset. Thanks to the embedded physics in its building block, DMN requires substantially less trainingdata than other machine learning problems relying on deep neural networks. According to our experiments,400 training samples are sufficient to achieve comparable training and test errors.A cost function based on the mean square error (MSE) is formulated to quantify the distance betweenthe DMN prediction and the training reference, J = 12 N s N s (cid:88) s =1 || ¯ C dnss − ¯ C rve || || ¯ C dnss || , (A.2)where ¯ C dnss is the overall stiffness matrix of s -th sample computed from DNS, and N s is the total numberof training samples. The operator || ... || denotes the Frobenius matrix norm.Since the homogenization and rotation operations in the two-layer building block both have analyticalfunctions, derivatives of the cost function with respect to all the fitting parameters in Eq. (A.1) can bederived following the chain rule in a backward propagation process. To accelerate the training speed, we usethe stochastic gradient descent (SGD) to optimize the cost function, and compression algorithms are alsointroduced to prune and merge the binary-tree network.For a given microstructure, the fitting parameters are initialized randomly following a uniform distribu-tion. However, if one has pre-trained networks for a similar microstructure, transfer learning can be used toease the training. As mentioned in our previous discussions, another application of transfer learning is tobuild a unified database of a set of microstructures within a material design space. As all the microstructures32 𝜃 𝜃 𝜋 − 𝜃𝑑 ! 𝑑 ! 𝜀 """ = 𝜀 "" 𝑑 ! −𝑑 ! Mat. 1Mat. 2 Mat. 1Mat. 2 𝐹 """ = 𝐹 "" , 𝐹 = 𝐹 Kinematic constraints: 𝜀 """ = 𝜀 "" 𝐹 """ = 𝐹 "" , 𝐹 ≠ 𝐹 Kinematic constraints:(a) (b)
Figure B.21: Small and finite formulations for two crack configurations of a 2-D building block. share the same topology structure through the initialization, DMNs for new intermediate microstructurescan be created by directly interpolating the fitting parameters.We implement the sampling, back-propagation, and optimization algorithms in Python, whereas theonline module for concurrent multiscale simulations is realized in FORTRAN. The off-shelf automatic dif-ferentiation functions in TensorFlow or PyTorch could also be utilized for model training more conveniently.
Appendix B. Extension to finite-strain formulation
Let us consider a 2-D building block with one crack surface in each child material, as shown in FigureB.21. Assume no normal displacement in the cohesive layers and no strain in the base material. Two cohesivelayers share the same reciprocal length parameter v c and the same magnitude of sliding displacement d t > ε = v c R c d = v cos θ sin θ √ θ cos θ sin θ cos θ −√ θ cos θ −√ θ cos θ √ θ cos θ cos θ − sin θ d t / √ = 12 vd t sin 2 θ − sin 2 θ √ θ . (B.1)with the notation ε = (cid:2) ε ε ε (cid:3) T . The kinematic constraint of the building block poses ε = ε (B.2)at the interface of the two child materials. As we can see Eq. (B.1), both configurations in Figure B.21 willsatisfy the kinematic constraint. Although the second case with two crack surfaces symmetric to the interfacerepresents a non-physical crack geometry, it cannot be distinguished in the small-strain formulation.This issue can be resolved by describing the system in the finite-strain formulation. Similarly, thedeformation gradient of a child material can be written as F = I + v c R fc ˜ d f = I + v c cos θ sin θ sin θ cos θ sin θ cos θ sin θ cos θ − sin θ cos θ − sin θ cos θ − sin θ cos θ sin θ cos θ cos θ − sin θ − sin θ cos θ sin θ cos θ − sin θ cos θ d t F = I + 12 v c d t sin 2 θ − sin 2 θ cos 2 θ + 1cos 2 θ − (B.3)with F = (cid:2) F F F F (cid:3) T and I = (cid:2) (cid:3) T . The kinematic constraints at the interface are F = F , F = F , (B.4)33hich only prefer the first crack configuration in Figure B.21 (a) for θ ∈ (0 , π ). Therefore, if one wants todetermine the exact crack orientations more than the overall material responses, the finite-strain formulationshould be used in general.In the remaining part of this section, we will briefly discuss the failure algorithms in DMN under thefinite-strain formulation. Assume the deformation gradient and the first Piola-Kirchhoff (PK1) stress of thebase material are F and P , respectively. The unit normal of a surface in the deformed configuration n isgiven by n = F − T N | F − T N | , (B.5)where N is the surface’s unit normal in the undeformed configuration. In contrast to Eq. (3.17), the planewith the maximum effective traction is obtained from N c = arg max N (cid:48) t m (cid:32) P · N (cid:48) , F − T N (cid:48) | F − T N (cid:48) | (cid:33) . (B.6)The definitions of the effective traction t m can be found in Eq. (3.6) and (3.9). Accordingly, the algorithmsproposed in Section 3.2 for finding the potential crack surfaces need to be modified.Once the crack surface is activated, the traction-separation law can be computed similar as before. Forexample, if | d · n c | ≥
0, the traction force per unit undeformed area can be written as t = t m d m (cid:2) β d + (1 − β )( d · n c ) n c (cid:3) , (B.7)where n c is the normal to the crack surface in the deformed configuration. Different from the small-strainformulation, t is not only related to the displacement vector d , but also affected by the deformation gradient F in the base material. Therefore, the incremental form of the traction-separation law becomes∆ t = ˜ K | F ∆ d + ˜ K geo | d ∆ F + δ t (B.8)with ˜ K = ∂ t ∂ d , ˜ K geo = ∂ t ∂ n c d n c d F , (B.9)where the extra geometric stiffness matrix ˜ K geo comes from the finite deformation of the base material.Note that if the deformation gradient is vectorized as 9 × K geo has a shape of 3 × − N c determined from Eq. (B.6), a rotation ofthe coordinate system is needed so that it becomes N c = (cid:8) (cid:9) T . The equilibrium conditions at theinterface are ∆ t = ∆ P , ∆ t = ∆ P , ∆ t = ∆ P , (B.10)and the cohesive layer does not impose any displacement or stress constraints on the base material in theother directions. Given the stress-strain relation in the base material, ∆ F = D ∆ P + δ F , the contributionof the enriching cohesive layer to the total deformation gradient of the micro-cell can be derived as ∆ F c ∆ F c ∆ F c = v c ∆ d = v c ˜ K − ∆ P ∆ P ∆ P − ˜ K − ˜ K geo D ∆ P − ˜ K − ˜ K geo δ F − ˜ K − δ t . (B.11) Appendix C. Adaptive time step refinement
For the implicit analysis of DMN with softening, the time step should be small enough so that themacroscale stiffness tensor at the top node is positive-definite. Meanwhile, the system with both nonlinear34lasticity in the base material and softening of enriching cohesive layers may need a smaller incrementfor convergence at certain load steps. Since a concurrent multiscale simulation commonly has thousandsof DMN instances with extensive loading paths, the adaptive load step refinement becomes necessary forbetter robustness and efficiency.
Box C.1 Algorithm for adaptive time step refinement
1. Initialize the macroscale strain increment ∆ ε macro and time increment ∆ t , and set N sub = 1 , i sub = 02. Compute ∆ ε macrosub = ∆ ε macro /N sub , and ∆ t sub = ∆ t/N sub
3. Use ∆ ε macrosub and ∆ t sub as the boundary conditions and time increments for the DMN failure analysislisted in Box 3.3.4. If the Newton’s method and the global crack configuration converge, i sub ← i sub + 1; else, N sub ← N sub , i sub = 2 i sub −
15. If i sub = N sub , the load step is completed; else, go to 2.The initial/critical time step is determined by the dilatational wave speed and the shortest size of anelement. For a general anisotropic material in either 2-D or 3-D spaces, the largest value in C macro , C macro ,and C macro of the macroscale stiffness tensor can be used to estimate the wave speed safely. Since the UDcomposite is stiffer in the fiber direction, the critical time step of Case 2 is smaller than the one of Case 1for the 2-D three-point bending tests in Section 4.3. References [1] 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.[2] 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.[3] J. Ghaboussi, J. Garrett Jr, X. Wu, Knowledge-based modeling of material behavior with neural networks, Journal ofengineering mechanics 117 (1) (1991) 132–153.[4] B. HKDH, Neural networks in materials science, ISIJ international 39 (10) (1999) 966–979.[5] Y. LeCun, Y. Bengio, et al., Convolutional networks for images, speech, and time series, The handbook of brain theoryand neural networks 3361 (10) (1995) 1995.[6] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural computation 9 (8) (1997) 1735–1780.[7] K. Cho, B. Van Merri¨enboer, D. Bahdanau, Y. Bengio, On the properties of neural machine translation: Encoder-decoderapproaches, arXiv preprint arXiv:1409.1259.[8] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Actametallurgica 21 (5) (1973) 571–574.[9] J. Qu, The effect of slightly weakened interfaces on the overall elastic properties of composite materials, Mechanics ofMaterials 14 (4) (1993) 269–281.[10] R. Hill, A self-consistent mechanics of composite materials, Journal of the Mechanics and Physics of Solids 13 (4) (1965)213–222.[11] R. Christensen, K. Lo, Solutions for effective shear properties in three phase sphere and cylinder models, Journal of theMechanics and Physics of Solids 27 (4) (1979) 315–330.[12] 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.[13] L. Liu, Solutions to the eshelby conjectures, Proceedings of the Royal Society A: Mathematical, Physical and EngineeringSciences 464 (2091) (2008) 573–594.[14] V. Kouznetsova, M. G. D. Geers, W. A. M. Brekelmans, Multi-scale constitutive modelling of heterogeneous materials witha gradient-enhanced computational homogenization scheme, International Journal for Numerical Methods in Engineering54 (8) (2002) 1235–1260.[15] T. Belytschko, S. Loehnert, J.-H. Song, Multiscale aggregating discontinuities: A method for circumventing loss of materialstability, International Journal for Numerical Methods in Engineering 73 (6) (2008) 869–894.[16] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complexmicrostructure, Computer Methods in Applied Mechanics and Engineering 157 (1–2) (1998) 69 – 94.
17] F. Feyel, J.-L. Chaboche, Fe 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ticomposite materials, Computer methods in applied mechanics and engineering 183 (3) (2000) 309–330.[18] F. Feyel, A multilevel finite element method (fe2) to describe the response of highly non-linear structures using generalizedcontinua, Computer Methods in Applied Mechanics and Engineering 192 (28–30) (2003) 3233 – 3244, ¡ce:title¿MultiscaleComputational Mechanics for Materials and Structures¡/ce:title¿.[19] J. Kochmann, S. Wulfinghoff, L. Ehle, J. Mayer, B. Svendsen, S. Reese, Efficient and accurate two-scale fe-fft-basedprediction of the effective material behavior of elasto-viscoplastic polycrystals, Computational Mechanics 61 (6) (2018)751–764.[20] R. De Borst, J. Pamin, R. Peerlings, L. Sluys, On gradient-enhanced damage and plasticity models for failure in quasi-brittle and frictional materials, Computational Mechanics 17 (1-2) (1995) 130–141.[21] R. H. Peerlings, R. de Borst, W. M. Brekelmans, J. De Vree, Gradient enhanced damage for quasi-brittle materials,International Journal for numerical methods in engineering 39 (19) (1996) 3391–3403.[22] M. Geers, R. De Borst, W. Brekelmans, R. Peerlings, Strain-based transient-gradient damage model for failure analyses,Computer methods in applied mechanics and engineering 160 (1-2) (1998) 133–153.[23] E. Kuhl, E. Ramm, R. de Borst, An anisotropic gradient damage model for quasi-brittle materials, Computer Methods inApplied Mechanics and Engineering 183 (1-2) (2000) 87–103.[24] Z. P. Bazant, M. Jir´asek, Nonlocal integral formulations of plasticity and damage: survey of progress, Journal of Engi-neering Mechanics 128 (11) (2002) 1119–1149.[25] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational prin-ciples and multi-field fe implementations, International Journal for Numerical Methods in Engineering 83 (10) (2010)1273–1311.[26] C. Miehe, M. Hofacker, L.-M. Schaenzel, F. Aldakheel, Phase field modeling of fracture in multi-physics problems. partii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic–plastic solids, Computer Methods inApplied Mechanics and Engineering 294 (2015) 486–522.[27] Z. P. Baˇzant, B. H. Oh, Crack band theory for fracture of concrete, Mat´eriaux et construction 16 (3) (1983) 155–177.[28] A. Gorgogianni, J. Eli´aˇs, J.-L. Le, Mechanism-based energy regularization in computational modeling of quasibrittlefracture, Journal of Applied Mechanics 87 (9).[29] G. T. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, International Journal of solidsand structures 33 (20-22) (1996) 2899–2938.[30] 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.[31] 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.[32] N. Mo¨es, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International journalfor numerical methods in engineering 46 (1) (1999) 131–150.[33] N. Mo¨es, T. Belytschko, Extended finite element method for cohesive crack growth, Engineering fracture mechanics 69 (7)(2002) 813–833.[34] J. Yvonnet, Q.-C. He, The reduced model multiscale method (r3m) for the non-linear homogenization of hyperelasticmedia at finite strains, Journal of Computational Physics 223 (1) (2007) 341–368.[35] F. Fritzen, O. Kunc, Two-stage data-driven homogenization for nonlinear solids using a reduced order model, EuropeanJournal of Mechanics-A/Solids 69 (2018) 201–220.[36] I. Rocha, F. van der Meer, L. Sluys, An adaptive domain-based pod/ecm hyper-reduced modeling framework withoutoffline training, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112650.[37] 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.[38] 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.[39] J. Gao, M. Shakoor, G. Domel, M. Merzkirch, G. Zhou, D. Zeng, X. Su, W. K. Liu, Predictive multiscale modeling forunidirectional carbon fiber reinforced polymers, Composites Science and Technology 186 (2020) 107922.[40] Z. Chen, T. Huang, Y. Shao, Y. Li, H. Xu, K. Avery, D. Zeng, W. Chen, X. Su, Multiscale finite element modeling of sheetmolding compound (smc) composite structure based on stochastic mesostructure reconstruction, Composite Structures188 (2018) 25–38.[41] R. Bostanabad, B. Liang, J. Gao, W. K. Liu, J. Cao, D. Zeng, X. Su, H. Xu, Y. Li, W. Chen, Uncertainty quantificationin multiscale simulation of woven fiber composites, Computer Methods in Applied Mechanics and Engineering 338 (2018)506–532.[42] T. Kirchdoerfer, M. Ortiz, Data-driven computational mechanics, Computer Methods in Applied Mechanics and Engi-neering 304 (2016) 81–101.[43] R. Ibanez, E. Abisset-Chavanne, J. V. Aguado, D. Gonzalez, E. Cueto, F. Chinesta, A manifold learning approach todata-driven computational elasticity and inelasticity, Archives of Computational Methods in Engineering 25 (1) (2018)47–57.[44] R. Eggersmann, T. Kirchdoerfer, S. Reese, L. Stainier, M. Ortiz, Model-free data-driven inelasticity, Computer Methodsin Applied Mechanics and Engineering 350 (2019) 81–99.[45] Q. He, J.-S. Chen, A physics-constrained data-driven approach based on locally convex reconstruction for noisy database,Computer Methods in Applied Mechanics and Engineering 363 (2020) 112791.[46] B. Le, J. Yvonnet, Q.-C. He, Computational homogenization of nonlinear elastic materials using neural networks, Inter- ational Journal for Numerical Methods in Engineering 104 (12) (2015) 1061–1084.[47] F. Fritzen, M. Fern´andez, F. Larsson, On-the-fly adaptivity for nonlinear twoscale simulations using artificial neuralnetworks and reduced order modeling, Frontiers in Materials 6 (2019) 75.[48] X. Lu, D. G. Giovanis, J. Yvonnet, V. Papadopoulos, F. Detrez, J. Bai, A data-driven computational homogenizationmethod based on neural networks for the nonlinear anisotropic electrical response of graphene/polymer nanocomposites,Computational Mechanics (2018) 1–15.[49] 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.[50] A. L. Frankel, R. E. Jones, C. Alleman, J. A. Templeton, Predicting the mechanical response of oligocrystals with deeplearning, Computational Materials Science 169 (2019) 109099.[51] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, M. Bessa, Deep learning predicts path-dependent plasticity,Proceedings of the National Academy of Sciences 116 (52) (2019) 26414–26420.[52] H. J. Logarzo, G. Capuano, J. J. Rimoli, Smart constitutive laws: Inelastic homogenization through machine learning,Computer Methods in Applied Mechanics and Engineering 373 (2020) 113482.[53] 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.[54] Z. Liu, H. Wei, T. Huang, C. Wu, Intelligent multiscale simulation based on process-guided composite database, arXivpreprint arXiv:2003.09491.[55] S. Gajek, M. Schneider, T. B¨ohlke, On the micromechanics of deep material networks, Journal of the Mechanics andPhysics of Solids (2020) 103984.[56] 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.[57] Z. Liu, Deep material network with cohesive layers: Multi-stage training and interfacial failure analysis, Computer Methodsin Applied Mechanics and Engineering 363 (2020) 112913.[58] T. Belytschko, J.-H. Song, Coarse-graining of multiscale crack propagation, International journal for numerical methodsin engineering 81 (5) (2010) 537–563.[59] Z. P. Bazant, Can multiscale-multiphysics methods predict softening damage and structural failure?, International Journalfor Multiscale Computational Engineering 8 (1).[60] 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.[61] Z. P. Baˇzant, B. H. Oh, Crack band theory for fracture of concrete, Mat´eriaux et construction 16 (3) (1983) 155–177.[62] Y. Wang, H. Waisman, From diffuse damage to sharp cohesive cracks: A coupled xfem framework for failure analysis ofquasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 299 (2016) 57–89.[63] C. Wu, Y. Wu, J. E. Crawford, J. M. Magallanes, Three-dimensional concrete impact and penetration simulations usingthe smoothed particle galerkin method, International Journal of Impact Engineering 106 (2017) 1–17.[64] B. Ren, C. Wu, E. Askari, A 3d discontinuous galerkin finite element method with the bond-based peridynamics modelfor dynamic brittle failure analysis, International Journal of Impact Engineering 99 (2017) 14–25.ational Journal for Numerical Methods in Engineering 104 (12) (2015) 1061–1084.[47] F. Fritzen, M. Fern´andez, F. Larsson, On-the-fly adaptivity for nonlinear twoscale simulations using artificial neuralnetworks and reduced order modeling, Frontiers in Materials 6 (2019) 75.[48] X. Lu, D. G. Giovanis, J. Yvonnet, V. Papadopoulos, F. Detrez, J. Bai, A data-driven computational homogenizationmethod based on neural networks for the nonlinear anisotropic electrical response of graphene/polymer nanocomposites,Computational Mechanics (2018) 1–15.[49] 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.[50] A. L. Frankel, R. E. Jones, C. Alleman, J. A. Templeton, Predicting the mechanical response of oligocrystals with deeplearning, Computational Materials Science 169 (2019) 109099.[51] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, M. Bessa, Deep learning predicts path-dependent plasticity,Proceedings of the National Academy of Sciences 116 (52) (2019) 26414–26420.[52] H. J. Logarzo, G. Capuano, J. J. Rimoli, Smart constitutive laws: Inelastic homogenization through machine learning,Computer Methods in Applied Mechanics and Engineering 373 (2020) 113482.[53] 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.[54] Z. Liu, H. Wei, T. Huang, C. Wu, Intelligent multiscale simulation based on process-guided composite database, arXivpreprint arXiv:2003.09491.[55] S. Gajek, M. Schneider, T. B¨ohlke, On the micromechanics of deep material networks, Journal of the Mechanics andPhysics of Solids (2020) 103984.[56] 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.[57] Z. Liu, Deep material network with cohesive layers: Multi-stage training and interfacial failure analysis, Computer Methodsin Applied Mechanics and Engineering 363 (2020) 112913.[58] T. Belytschko, J.-H. Song, Coarse-graining of multiscale crack propagation, International journal for numerical methodsin engineering 81 (5) (2010) 537–563.[59] Z. P. Bazant, Can multiscale-multiphysics methods predict softening damage and structural failure?, International Journalfor Multiscale Computational Engineering 8 (1).[60] 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.[61] Z. P. Baˇzant, B. H. Oh, Crack band theory for fracture of concrete, Mat´eriaux et construction 16 (3) (1983) 155–177.[62] Y. Wang, H. Waisman, From diffuse damage to sharp cohesive cracks: A coupled xfem framework for failure analysis ofquasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 299 (2016) 57–89.[63] C. Wu, Y. Wu, J. E. Crawford, J. M. Magallanes, Three-dimensional concrete impact and penetration simulations usingthe smoothed particle galerkin method, International Journal of Impact Engineering 106 (2017) 1–17.[64] B. Ren, C. Wu, E. Askari, A 3d discontinuous galerkin finite element method with the bond-based peridynamics modelfor dynamic brittle failure analysis, International Journal of Impact Engineering 99 (2017) 14–25.