Exploring the 3D architectures of deep material network in data-driven   multiscale mechanics

Zeliang Liu; C.T. Wu

arXiv:1901.04832·cs.CE·July 7, 2020

Exploring the 3D architectures of deep material network in data-driven multiscale mechanics

Zeliang Liu, C.T. Wu

PDF

TL;DR

This paper advances the deep material network (DMN) to effectively model complex 3D multiscale heterogeneous materials with nonlinear behaviors, demonstrating high accuracy and efficiency through numerical validation on various composite systems.

Contribution

It introduces a new 3D multi-layer network structure and analytical solutions for building blocks, enhancing DMN's capability for multiscale nonlinear material modeling.

Findings

01

Effective training of 3D DMN using stochastic gradient descent.

02

Accurate modeling of complex multiscale materials with nonlinear behaviors.

03

Validation on diverse composite materials demonstrating robustness.

Abstract

This paper extends the deep material network (DMN) proposed by Liu et al. (2019) to tackle general 3-dimensional (3D) problems with arbitrary material and geometric nonlinearities. It discovers a new way of describing multiscale heterogeneous materials by a multi-layer network structure and mechanistic building blocks. The data-driven framework of DMN is discussed in detail about the offline training and online extrapolation stages. Analytical solutions of the 3D building block with a two-layer structure in both small- and finite-strain formulations are derived based on interfacial equilibrium conditions and kinematic constraints. With linear elastic data generated by direct numerical simulations on a representative volume element (RVE), the network can be effectively trained in the offline stage using stochastic gradient descent and advanced model compression algorithms. Efficiency and…

Tables7

Table 1. Table 1: Training results of the particle-reinforced composite. Average training error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , average test error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , maximum test error and predicted volume fraction v f 1 𝑣 subscript 𝑓 1 vf_{1} are provided for each DMN. Test errors of the linear FEM model are also shown.

	Epochs	Training ${\bar{e}}^{t r}$	Test ${\bar{e}}^{t e}$	Maximum $e_{s}^{t e}$	$v f_{1}$
$N = 4$	20000	7.61%	7.79%	17.9%	0.211 (-6.63%)
$N = 5$	20000	4.47%	4.49%	8.94%	0.220 (-2.65%)
$N = 6$	20000	1.34%	1.39%	4.46%	0.220 (-2.65%)
$N = 8$	40000	0.53%	0.59%	2.41%	0.224 (-0.88%)
Linear FEM	\	\	2.30%	12.5%	\

Table 2. Table 2: Material parameters for the hyperelastic particle-reinforced rubber composite with Mullins effect.

Matrix	$C_{10}$ (MPa)	$C_{01}$ (MPa)	$ν$	$η$	a (J $\cdot {cm}^{- 3}$ )	b
Particle	$C_{10}$ (MPa)	$C_{01}$ (MPa)	$ν$
Particle	100	0	0.3
Matrix	1.0	0.5	0.495	0.8	1.0	1.0

Table 3. Table 3: Training results of the polycrystalline RVE with random ODF. Average training error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , average test error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} and maximum test error are provided for each DMN.

	Epochs	Training ${\bar{e}}^{t r}$	Test ${\bar{e}}^{t e}$	Maximum $e_{s}^{t e}$
$N = 4$	20000	5.87%	6.00%	15.5%
$N = 6$	20000	1.16%	1.27%	3.64%
$N = 8$	20000	0.36%	0.43%	1.80%

Table 4. Table 4: Crystal plasticity parameters for a FCC metal with 12 slip planes [ 51 , 52 ] .

$C_{1111}$ (GPa)	$C_{1122}$ (GPa)	$C_{2323}$ (GPa)	${\dot{γ}}_{0}$ (s^-1)	$m$	initial $τ_{0}$ (MPa)
196.4	84.2	56.1	0.00242	58.8	171.85
$H$ (MPa)	$R$ (MPa)	$χ$	initial $a_{0}$ (MPa)	$h$ (MPa)	$r$ (MPa)
1.0	0.0	1.0	0.0	500.0	0.0

Table 5. Table 5: Training results of CFRP RVEs. Average training error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , average test error e ¯ t r superscript ¯ 𝑒 𝑡 𝑟 \bar{e}^{tr} , maximum test error and predicted volume fraction v f 1 𝑣 subscript 𝑓 1 vf_{1} are provided for each DMN. Test errors of the linear FEM model for woven RVE are also shown.

		Epochs	Training ${\bar{e}}^{t r}$	Test ${\bar{e}}^{t e}$	Maximum $e_{s}^{t e}$	$v f_{1}$
	$N = 5$	20000	11.1%	10.5%	47.2%	0.303 (+3.06%)
UD	$N = 7$	20000	0.79%	0.76%	3.46%	0.294 (+0.00%)
	$N = 9$	40000	0.23%	0.25%	1.34%	0.295 (+0.00%)
	$N = 4$	20000	34.6%	35.6%	83.9%	0.301 (-27.3%)
	$N = 6$	20000	11.6%	11.5%	38.1%	0.348 (-15.9%)
Woven	$N = 7$	20000	3.72%	3.48%	10.1%	0.424 (+2.42%)
	$N = 8$	20000	2.19%	2.15%	7.82%	0.420 (+1.45%)
	Linear FEM	\	\	6.12%	28.3%	\

Table 6. Table 6: Material parameters for the carbon fiber reinforced polymer (CFRP) composite [ 53 ] .

Carbon fiber	$E_{1}$ (GPa)	$E_{2}$ (GPa)	$E_{3}$ (GPa)	$G_{12}$ (GPa)	$G_{13}$ (GPa)
	19.8	19.8	245.0	5.9	29.2
	$G_{23}$ (GPa)	$ν_{12}$	$ν_{13}$	$ν_{23}$
	29.2	0.67	0.28	0.28
Epoxy	$E_{m}$ (GPa)	$ν_{m}$	$a_{1}$	$a_{2}$ (MPa)	$a_{3}$ (MPa)
Epoxy	3.8	0.387	200	10	20
Yarn (elastic)	$E_{1}$ (GPa)	$E_{2}$ (GPa)	$E_{3}$ (GPa)	$G_{12}$ (GPa)	$G_{13}$ (GPa)
	10.2	10.2	78.8	1.95	2.39
	$G_{23}$ (GPa)	$ν_{12}$	$ν_{13}$	$ν_{23}$
	2.39	0.60	0.35	0.35

Table 7. Table 7: Offline computational times for the particle-reinforced RVE.

	Training data generation	DMN training (20000 epochs)
	DNS (400 samples)	$N = 4$	$N = 6$	$N = 8$
$N_{c p u}$	10	10	10	10
Wall time (h)	39.5	5.4	16.7	43.0

Equations208

C_{N}^{o dd} = C^{p 1}, C_{N}^{e v e n} = C^{p 2} .

C_{N}^{o dd} = C^{p 1}, C_{N}^{e v e n} = C^{p 2} .

\overset{ˉ}{C}^{r v e} = \overset{ˉ}{C}_{1}^{1} .

\overset{ˉ}{C}^{r v e} = \overset{ˉ}{C}_{1}^{1} .

w_{N}^{j} = a (z^{j}) = max (z^{j}, 0),

w_{N}^{j} = a (z^{j}) = max (z^{j}, 0),

w_{i}^{k} = w_{i + 1}^{2 k - 1} + w_{i + 1}^{2 k} .

w_{i}^{k} = w_{i + 1}^{2 k - 1} + w_{i + 1}^{2 k} .

w_{i}^{k} = j = 2^{N - i} (k - 1) + 1 \sum 2^{N - i} k w_{N}^{j} = j = 2^{N - i} (k - 1) + 1 \sum 2^{N - i} k a (z^{j}) .

w_{i}^{k} = j = 2^{N - i} (k - 1) + 1 \sum 2^{N - i} k w_{N}^{j} = j = 2^{N - i} (k - 1) + 1 \sum 2^{N - i} k a (z^{j}) .

\frac{\partial w _{i}^{k}}{\partial z ^{j}} = {a^{'} (z^{j}) if k = ⌈ j / 2^{N - i} ⌉ 0 otherwise,

\frac{\partial w _{i}^{k}}{\partial z ^{j}} = {a^{'} (z^{j}) if k = ⌈ j / 2^{N - i} ⌉ 0 otherwise,

Output \overset{ˉ}{C}^{r v e} = f_{2} (Inputs C^{p 1}, C^{p 2}, z^{j = 1, 2, ..., 2^{N - 1}}, α_{i = 1, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, β_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, γ_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}} Fitting parameters) .

Output \overset{ˉ}{C}^{r v e} = f_{2} (Inputs C^{p 1}, C^{p 2}, z^{j = 1, 2, ..., 2^{N - 1}}, α_{i = 1, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, β_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, γ_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}} Fitting parameters) .

Output \overset{ˉ}{C}^{r v e} = f_{1} (Input C^{p 1}, z^{j = 1, 2, ..., 2^{N - 1}}, α_{i = 1, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, β_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, γ_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}} Fitting parameters) .

Output \overset{ˉ}{C}^{r v e} = f_{1} (Input C^{p 1}, z^{j = 1, 2, ..., 2^{N - 1}}, α_{i = 1, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, β_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}}, γ_{i = 1, 2, ..., N}^{k = 1, 2, ..., 2^{i - 1}} Fitting parameters) .

J (z, α, β, γ) = \frac{1}{2 N _{s}} s \sum J_{s} (z, α, β, γ) + λ (j \sum a (z^{j}) - 2^{N - 2})^{2} ，

J (z, α, β, γ) = \frac{1}{2 N _{s}} s \sum J_{s} (z, α, β, γ) + λ (j \sum a (z^{j}) - 2^{N - 2})^{2} ，

J_{s} (z, α, β, γ) = \frac{∣∣ C ˉ _{s}^{d n s} - f _{2} ( C _{s}^{p 1} , C _{s}^{p 2} , z , α , β , γ ) ∣ ∣ ^{2}}{∣∣ C ˉ _{s}^{d n s} ∣ ∣ ^{2}} for a two-phase RVE,

J_{s} (z, α, β, γ) = \frac{∣∣ C ˉ _{s}^{d n s} - f _{2} ( C _{s}^{p 1} , C _{s}^{p 2} , z , α , β , γ ) ∣ ∣ ^{2}}{∣∣ C ˉ _{s}^{d n s} ∣ ∣ ^{2}} for a two-phase RVE,

J_{s} (z, α, β, γ) = \frac{∣∣ C ˉ _{s}^{d n s} - f _{1} ( C _{s}^{p 1} , z , α , β , γ ) ∣ ∣ ^{2}}{∣∣ C ˉ _{s}^{d n s} ∣ ∣ ^{2}} for a single-phase RVE .

J_{s} (z, α, β, γ) = \frac{∣∣ C ˉ _{s}^{d n s} - f _{1} ( C _{s}^{p 1} , z , α , β , γ ) ∣ ∣ ^{2}}{∣∣ C ˉ _{s}^{d n s} ∣ ∣ ^{2}} for a single-phase RVE .

D^{p i} = ⎩ ⎨ ⎧ 1/ E_{11}^{p i} - ν_{12}^{p i} / E_{22}^{p i} 1/ E_{22}^{p i} - ν_{31}^{p i} / E_{11}^{p i} - ν_{23}^{p i} / E_{33}^{p i} 1/ E_{33}^{p i} 1/ (2 G_{23}^{p i}) 1/ (2 G_{31}^{p i}) 1/ (2 G_{12}^{p i}) ⎭ ⎬ ⎫ with i = 1, 2.

D^{p i} = ⎩ ⎨ ⎧ 1/ E_{11}^{p i} - ν_{12}^{p i} / E_{22}^{p i} 1/ E_{22}^{p i} - ν_{31}^{p i} / E_{11}^{p i} - ν_{23}^{p i} / E_{33}^{p i} 1/ E_{33}^{p i} 1/ (2 G_{23}^{p i}) 1/ (2 G_{31}^{p i}) 1/ (2 G_{12}^{p i}) ⎭ ⎬ ⎫ with i = 1, 2.

lo g_{10} (E_{11}^{p i}), lo g_{10} (E_{22}^{p i}), lo g_{10} (E_{33}^{p i}) \in U [- 1, 1],

lo g_{10} (E_{11}^{p i}), lo g_{10} (E_{22}^{p i}), lo g_{10} (E_{33}^{p i}) \in U [- 1, 1],

\overset{ˉ}{E}^{p 1} = 1, lo g_{10} (\overset{ˉ}{E}^{p 2}) \in U [- 3, 3] .

\overset{ˉ}{E}^{p 1} = 1, lo g_{10} (\overset{ˉ}{E}^{p 2}) \in U [- 3, 3] .

E_{k k}^{p i} \leftarrow \frac{E ˉ ^{p i}}{( E _{11}^{p i} E _{22}^{p i} E _{33}^{p i} ) ^{1/3}} E_{k k}^{p i} with k = 1, 2, 3.

E_{k k}^{p i} \leftarrow \frac{E ˉ ^{p i}}{( E _{11}^{p i} E _{22}^{p i} E _{33}^{p i} ) ^{1/3}} E_{k k}^{p i} with k = 1, 2, 3.

\frac{G _{12}^{p i}}{E _{11}^{p i} E _{22}^{p i}}, \frac{G _{23}^{p i}}{E _{22}^{p i} E _{33}^{p i}}, \frac{G _{31}^{p i}}{E _{33}^{p i} E _{11}^{p i}} \in U [0.25, 0.5] .

\frac{G _{12}^{p i}}{E _{11}^{p i} E _{22}^{p i}}, \frac{G _{23}^{p i}}{E _{22}^{p i} E _{33}^{p i}}, \frac{G _{31}^{p i}}{E _{33}^{p i} E _{11}^{p i}} \in U [0.25, 0.5] .

\frac{ν _{12}^{p i}}{E _{22}^{p i} / E _{11}^{p i}} \in U (0.0, 0.5), \frac{ν _{23}^{p i}}{E _{33}^{p i} / E _{22}^{p i}} \in U (0.0, 0.5), \frac{ν _{31}^{p i}}{E _{11}^{p i} / E _{33}^{p i}} \in U (0.0, 0.5) .

\frac{ν _{12}^{p i}}{E _{22}^{p i} / E _{11}^{p i}} \in U (0.0, 0.5), \frac{ν _{23}^{p i}}{E _{33}^{p i} / E _{22}^{p i}} \in U (0.0, 0.5), \frac{ν _{31}^{p i}}{E _{11}^{p i} / E _{33}^{p i}} \in U (0.0, 0.5) .

z^{j}^{(0)} \sim U (0.2, 0.8) and α_{i}^{k}^{(0)}, β_{i}^{k}^{(0)}, γ_{i}^{k}^{(0)} \sim U (- π /2, π /2) .

z^{j}^{(0)} \sim U (0.2, 0.8) and α_{i}^{k}^{(0)}, β_{i}^{k}^{(0)}, γ_{i}^{k}^{(0)} \sim U (- π /2, π /2) .

(R_{1})^{- 1} R_{2} = I .

(R_{1})^{- 1} R_{2} = I .

∣∣ R (λ_{i}) ∣ - 1∣ < t o l,

∣∣ R (λ_{i}) ∣ - 1∣ < t o l,

layer N forward (A, δ P) homogenization layer 1 (macroscale) backward (Δ F, Δ P) de-homogenization layer N .

layer N forward (A, δ P) homogenization layer 1 (macroscale) backward (Δ F, Δ P) de-homogenization layer N .

\overset{ˉ}{σ} = \overset{ˉ}{C} : \overset{ˉ}{ε} .

\overset{ˉ}{σ} = \overset{ˉ}{C} : \overset{ˉ}{ε} .

σ^{1} = \overset{ˉ}{C}^{1} : ε^{1}, σ^{2} = \overset{ˉ}{C}^{2} : ε^{2} .

σ^{1} = \overset{ˉ}{C}^{1} : ε^{1}, σ^{2} = \overset{ˉ}{C}^{2} : ε^{2} .

f_{1} = \frac{w ^{1}}{w ^{1} + w ^{2}}, f_{2} = 1 - f_{1} .

f_{1} = \frac{w ^{1}}{w ^{1} + w ^{2}}, f_{2} = 1 - f_{1} .

σ = {σ_{11}, σ_{22}, σ_{33}, 2 σ_{23}, 2 σ_{13}, 2 σ_{12}}^{T} = {σ_{1}, σ_{2}, σ_{3}, σ_{4}, σ_{5}, σ_{6}}^{T}

σ = {σ_{11}, σ_{22}, σ_{33}, 2 σ_{23}, 2 σ_{13}, 2 σ_{12}}^{T} = {σ_{1}, σ_{2}, σ_{3}, σ_{4}, σ_{5}, σ_{6}}^{T}

ε = {ε_{11}, ε_{22}, ε_{33}, 2 ε_{23}, 2 ε_{13}, 2 ε_{12}}^{T} = {ε_{1}, ε_{2}, ε_{3}, ε_{4}, ε_{5}, ε_{6}}^{T},

ε = {ε_{11}, ε_{22}, ε_{33}, 2 ε_{23}, 2 ε_{13}, 2 ε_{12}}^{T} = {ε_{1}, ε_{2}, ε_{3}, ε_{4}, ε_{5}, ε_{6}}^{T},

C = h_{C} (\overset{ˉ}{C}^{1}, \overset{ˉ}{C}^{2}, w^{1}, w^{2});

C = h_{C} (\overset{ˉ}{C}^{1}, \overset{ˉ}{C}^{2}, w^{1}, w^{2});

\overset{ˉ}{C} = r_{C} (C, α, β, γ) .

\overset{ˉ}{C} = r_{C} (C, α, β, γ) .

σ_{3}^{1} = σ_{3}^{2}, σ_{4}^{1} = σ_{4}^{2}, σ_{5}^{1} = σ_{5}^{2},

σ_{3}^{1} = σ_{3}^{2}, σ_{4}^{1} = σ_{4}^{2}, σ_{5}^{1} = σ_{5}^{2},

ε_{1}^{1} = ε_{1}^{2}, ε_{2}^{1} = ε_{2}^{2}, ε_{6}^{1} = ε_{6}^{2} .

ε_{1}^{1} = ε_{1}^{2}, ε_{2}^{1} = ε_{2}^{2}, ε_{6}^{1} = ε_{6}^{2} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Exploring the 3D architectures of deep material network in data-driven multiscale mechanics

Zeliang Liu

[email protected]

C.T. Wu

Abstract

This paper extends the deep material network (DMN) proposed by Liu et al. (2019) [1] to tackle general 3-dimensional (3D) problems with arbitrary material and geometric nonlinearities. It discovers a new way of describing multiscale heterogeneous materials by a multi-layer network structure and mechanistic building blocks. The data-driven framework of DMN is discussed in detail about the offline training and online extrapolation stages. Analytical solutions of the 3D building block with a two-layer structure in both small- and finite-strain formulations are derived based on interfacial equilibrium conditions and kinematic constraints. With linear elastic data generated by direct numerical simulations on a representative volume element (RVE), the network can be effectively trained in the offline stage using stochastic gradient descent and advanced model compression algorithms. Efficiency and accuracy of DMN on addressing the long-standing 3D RVE challenges with complex morphologies and material laws are validated through numerical experiments, including 1) hyperelastic particle-reinforced rubber composite with Mullins effect; 2) polycrystalline materials with rate-dependent crystal plasticity; 3) carbon fiber reinforced polymer (CFRP) composites with fiber anisotropic elasticity and matrix plasticity. In particular, we demonstrate a three-scale homogenization procedure of CFRP system by concatenating the microscale and mesoscale material networks. The complete learning and extrapolation procedures of DMN establish a reliable data-driven framework for multiscale material modeling and design.

keywords:

Machine learning, 3D building-block, hyperelasticity, crystal plasticity, CFRP composites, three-scale homogenization

1 Introduction

Modern material systems with properly designed microstructures offer new material engineering avenues for producing advantageous mechanical properties and functionalities in various applications. In aerospace and automotive industries, carbon fiber reinforced polymer (CFRP) composites have become attractive alternatives of metal materials due to their high strength/weight ratio induced by the interplay between carbon fibers and epoxy matrix [2, 3, 4]. Another typical example of heterogeneous materials is the particle-reinforced rubber composite where nanoparticles are added into the matrix to manipulate the overall mechanical properties, such as stiffness and viscoelastic behaviors [5, 6, 7, 8]. Moreover, in metal additive manufacturing, the performance of the final product is strongly affected by the microscopic polycrystalline microstructure, which could be controlled by the deposition and cooling processes [9, 10]. Efficient and physical descriptions of such multiscale materials will not only help to predict the performance of large-scale structures, but also accelerate the pace of discovery and design of new material systems.

The multiscale nature of heterogeneous materials poses limitations on single-scale empirical models, which lose the sight of microstructural interactions. They tend to fail in attempting to capture nonlinear or anisotropic responses, and may require burdensome calibration to find the model parameters. As a result, homogenization based on the concept of representative volume element (RVE) [11] has become an important approach to model multiscale materials [12]. Many analytical methods have been proposed, which adopt some micromechanics assumptions to simplify the full-field RVE problem, such as Hashin-Shtrikman bounds [13, 14], the Mori-Tanaka method [15] and self-consistent methods[16]. Since most analytical methods are derived based on the solution for regular geometries and simple material models (e.g. Eshelby’s solution for isotropic elastic materials [17]), it is usually difficult for them to consider complex microstructural morphologies, history-dependent materials and large deformations. Additionally, direct numerical simulation (DNS) tools, such as finite element [18], meshfree [19, 20] and fast Fourier transform (FFT)-based micromechanics methods [21, 22], are both flexible and accurate. However, a DNS model involves a detailed meshing of the RVE microstructures and requires tremendous computational cost, especially for 3D problems. Therefore, one of the fundamental issues in multiscale material modeling and design is how to find an accurate low-dimensional representation of the RVE for arbitrary morphologies and nonlinearities.

In the past decade, a plethora of data-driven material modeling methods have been proposed based on existing machine learning techniques. Depending on the type of training data in the offline stage, we classify these methods mainly into macroscopic and microscopic approaches. In macroscopic approaches, the stress-strain relations, or strain energy density functions, are directly fitted by regression methods, like deep neural network (DNN) [23, 24, 25, 26] and Kriging methods [26, 27]. Enabled by recent progresses in computer hardware systems, DNN becomes one of the most popular tools due to its large model generalities [28, 29], and has also stimulated applications across different engineering disciplines [30, 31, 32, 33]. However, the extrapolation capability of the macroscopic approaches to unknown material and loading spaces is usually limited by the lack of microscale physics. To overcome this problem, microscopic approaches take the full RVE stress or strain fields as the training data and use dimension-reduction techniques to find the reduced basis, such as non-uniform transformation field analysis[34, 35], proper orthogonal decomposition [36, 37], manifold learning methods [38, 39] and self-consistent clustering analysis (SCA)[40, 41, 42]. One potential limitation of the microscopic approaches is the difficulty of obtaining stress or strain fields from experiments, though advances in digital image correlation methods [43, 44] may provide a solution.

We recently proposed a novel data-driven material modeling method called “deep material network” (DMN) [1] and demonstrated its effectiveness for various challenging 2-dimensional RVE problems, including matrix-inclusion composite, amorphous material and anisotropic material with penetrating phase. It discovers a reduced representation of RVE based on a binary-tree network of mechanistic two-layer building blocks. Comparing to other data-driven methods, DMN has the following intriguing features: 1) Avoiding an extensive offline sampling stage; 2) eliminating the need for extra calibration and micromechanics assumption; 3) efficient online prediction without the danger of extrapolation [1]. With the physics embedded in the DMN model, we are able to extract essential microstructural information (e.g. volume fractions) from pure macroscopic mechanical data, which can also be measured from experiments. In the online extrapolation stage, efficient and accurate predictions are achieved for unknown materials, such as linear elasticity with high contrast of phase properties, nonlinear history-dependent plasticity and finite-strain hyperelasticity under large deformations. Different from methods based on solving the continuum equilibrium equation of the RVE, either in partial differential from (e.g. FEM) or integral forms (e.g. FFT-based methods and SCA), DMN uses a hierarchical network structure to propagate all physical quantities, such as the stress and strain. As a result, its computational time is proportional to the number of degrees of freedom (DOFs), which is advantageous over most existing homogenization methods.

This paper adds an important piece to the DMN framework to treat general 3-dimensional (3D) problems with both material and geometric nonlinearities. The global framework of DMN, including the offline and online stages, is introduced in Section 2. Theories of the 3D mechanistic building block are derived in Section 3. We apply the DMN to three representative multiscale material systems in Section 4, including the particle-reinforced rubber composite, polycrystalline materials and CFRP composites. Moreover, a three-scale homogenization procedure of CFRP material system is accomplished by simple concatenation of microscale and mesoscale material networks. We provide the computational cost of DMN in the offline and online stages in Section 5, and further discuss its efficiency comparing to DNS and reference FE models. Concluding remarks are given in Section 6.

2 The global framework of deep material network

2.1 Preliminaries

The global framework of deep material network is presented in Figure 1. As a common practice of data-driven methods, an offline stage is first performed to collect the training datasets and search for the optimum model that represents the input-output relation. We start from a given RVE mesh, and for each sample $s$ , we assign a set of random elastic stiffness matrices to microscale material phases as the inputs. For a two-phase RVE, these inputs are denoted by $\textbf{C}^{p1}_{s}$ and $\textbf{C}^{p2}_{s}$ . Design of experiments (DoE) is performed to effectively explore the high-dimensional sampling space of input variables. The full-field RVE models can be then analyzed by various DNS tools, such as finite element method and FFT-based method. In general, six orthogonal loading conditions need to be simulated for each sample in order to output all the components in an elastic stiffness matrix $\bar{\textbf{C}}^{dns}_{s}$ . Other than the training dataset, it is always important to prepare the test dataset with extra points sampled in the same design space, which can be used later on to check the quality of a fitted model for generalization.

Mechanistic machine learning. Once the training dataset is collected, it can be evaluated by machine learning algorithms to train the parameters in DMN. As shown in Figure 1, the fitting parameters in the model are activation $z$ and rotation angles $\alpha,\beta,\gamma$ , all of which have physical meanings related to microstructural geometry and mechanics, as will be discussed in Section 2.2. The training process is formulated as a supervised optimization problem, where the objective is to minimize the distance between the stiffness matrix predicted by the DMN and the one from DNS. Gradient-based methods are adopted to update the fitting parameters. Meanwhile, several network compression algorithms are also introduced to reduced the number of fitting parameters and improve the convergence speed during the training. Details on the data generation and machine learning of DMN can be found in Section 2.3.

Microstructural database. The trained DMN discovers a reduced topological representation of the original full-field RVE with fewer DOFs. Based on the mechanistic building block with embedded physics, it can be extrapolated to unknown materials and loading paths outside the original linear elastic sampling space, without violating any physical law in the homogenization. By virtue of its unique capability of extrapolation, DMN creates a so-called “microstructural database”, which has the potential to drive high-fidelity concurrent multiscale simulations. Its intrinsically parameterized structure also provides a mechanistic understanding of structure-property relations and enables material design across length scales.

2.2 Multi-layer network structure

We propose to use a collection of connected mechanistic building blocks to represent the RVE model. As shown in Figure 2, a binary-tree structure is chosen as the base of DMN, where each node has two child nodes. The top layer is denoted as Layer $1$ , and the bottom layer is denoted as Layer $N$ . Given a network with depth $N$ , there are $2^{N-1}$ nodes in the bottom layer and in total $(2^{N}-1)$ nodes within the whole network. Nodes can be deactivated during the training process, as a result, some nodes may only have one child node in the final network structure.

Physics of building blocks Other than the binary-tree network structure, another key ingredient of DMN is the use of mechanistic building blocks. The rule of a building block is similar to the one of “neurons” in artificial neural network, however, the transformation functions in the building block are derived based on the analytical solutions of mechanistic structures with physical meanings. Generally speaking, each building block contains one homogenization operation and one rotation operation. As will be shown in Section 3, the building block is proposed as a 3D two-layer structure, which is probably the simplest form of a heterogeneous composite. The material phases interact through the interface between the two layers. The effect of phase fraction is captured by the homogenization operation, which is related to the fitting parameter $z$ . To increase the flexibility of the DMN representation, the rotation operation. controlled by the fitting parameters ( $\alpha,\beta,\gamma$ ), is introduced to vary the global orientation of the block. Our hypothesis is that the original 3D RVE can be represented by patterns of the two-layer structures with unknown phase fractions and block orientations. At layer $i$ , the $k$ -th building block is denoted as $\mathfrak{B}_{i}^{k}$ . We use C to represent a stiffness matrix before the rotation operation, and the one after rotation is denoted by $\bar{\textbf{C}}$ .

Mechanical data flow. Inputs and outputs of current DMN models are all mechanical data, such as stiffness matrix and residual stress. In particular, a material network in the offline stage takes the stiffness matrices of microscale phases as the inputs and output the overall RVE stiffness matrix $\bar{\textbf{C}}^{rve}$ . Illustrations of how the stiffness matrices propagate from the input layer to the output layer are depicted in Figure 2 (a) for a two-phase material. The inputs $\textbf{C}^{p1}$ and $\textbf{C}^{p2}$ are first assigned to the nodes in the bottom layer $N$ in a scheme designed as below,

[TABLE]

At the end of forward propagation, the output of DMN is obtained as

[TABLE]

The framework can be generalized to single-phase or multi-phase materials. For a RVE with one single material phase, like the polycrystalline RVE that will be investigated in Section 4.2, the whole framework of DMN does not need to be altered, except that all the bottom-layer nodes in the network receive the same stiffness matrix $\textbf{C}^{p1}$ . Extensions to RVE problems with more than two material phases are also possible, such as redesigning the network structure to let each node take more child nodes, or introducing advanced schemes of assigning multiple material phases to the bottom layer without changing binary-tree structure. Nevertheless, we will focus on single-phase and two-phase RVEs in this paper. In addition, to consider both material and geometric nonlinearities in the online stage, we also need to add the residual stress to the mechanical data flow, as will be discussed in Section 2.4.

Geometric data flow. As shown in Figure 2 (b), the weights $w$ is added into the data flow to keep track of the physical portion of each node in the network. We use $w_{i}^{k}$ to denote the weight of the $k$ -th node at Layer $i$ . An important feature of DMN is that the weights in the bottom layer are activated through the rectified linear unit (ReLU) [45]. For $j\in[1,2^{N-1}]$ , we have

[TABLE]

where $a$ is the ReLU activation function and $z$ is called the “activation”. It should be noted that once a unit is deactivated ( $z^{j}<0$ ), the unit will never be activated again during the training process due to its vanishing gradient ( $a^{\prime}(z^{j})=0$ ). This feature is helpful for automatically simplifying the material network and increasing the training speed. Without losing the physical meaning of $w$ , weights of two child nodes are summed up and passed to their parent node. From Layer $i+1$ to $i$ , we have

[TABLE]

By performing the summation recursively, any weight $w_{i}^{k}$ can be expressed as a summation of weights of its descendant nodes at the bottom layer $N$ ,

[TABLE]

The activations $z^{j}$ determine all the weights in the network, hence, they are chosen as the independent fitting parameters in DMN. Additionally, the partial derivative of $w_{i}^{k}$ with respect to $z^{j}$ can be written as

[TABLE]

where the operator $\lceil\;\rceil$ rounds a number to the next larger integer.

2.3 Offline training based on linear elasticity

Let us zoom into the building block $\mathfrak{B}_{i}^{k}$ appeared in the offline stage. As illustrated in Figure 3, the building block has two inputs $\bar{\textbf{C}}_{i+1}^{2k-1}$ and $\bar{\textbf{C}}_{i+1}^{2k}$ , and there are two operations within each building block: Homogenization and rotation. The homogenization operation is defined by the transformation function $\boldsymbol{\mathfrak{h}}_{c}$ , which is determined by the underlying physical structure and takes geometric descriptors (e.g. the weights $w$ ) as parameters. The rotation operation, described by the function $\boldsymbol{\mathfrak{r}}_{c}$ , takes the stiffness matrix $\textbf{C}_{i}^{k}$ generated by the homogenization operation and rotates it based on the angles $\alpha_{i}^{k},\beta_{i}^{k}$ and $\gamma_{i}^{k}$ . The result after this 3D rotation is $\bar{\textbf{C}}_{i}^{k}$ . Explicit expressions of $\boldsymbol{\mathfrak{h}}_{c}$ and $\boldsymbol{\mathfrak{r}}_{c}$ for a generic two-layer structure can be found in Section 3.1. Note that a building block at the bottom layer is slightly different from ones in other layers, since it only takes one input and the homogenization operation is ignored as depicted in Figure 2 (a).

Each building block generates a more flexible material model than its inputs. In this way, a building block in the upper layer can capture material responses at a more complex level than one in the lower layer. Based on Eq. (2.1), (2.2) and (2.5), the overall stiffness matrix of the multi-layer network $\bar{\textbf{C}}^{rve}$ for a two-phase RVE can be written as a function of the stiffness matrix from each material phase ( $\textbf{C}^{p1}$ and $\textbf{C}^{p2}$ ) and the fitting parameters ( $z,\alpha,\beta,\gamma$ ),

[TABLE]

Similarly, the overall function of DMN for a single-phase RVE becomes

[TABLE]

For a network with depth $N$ , there are totally $(7\times 2^{N-1}-3)$ fitting parameters for a 3D problem. However, the actual number of active fitting parameters will be diminished on the fly during the training process due to nodal deactivation and model compression.

Cost function. A cost function based on the mean square error (MSE) and a regularization term [1] is formulated to characterize how close is the model prediction $\bar{\textbf{C}}^{rve}$ to the training reference $\bar{\textbf{C}}^{dns}$ :

[TABLE]

with

[TABLE]

or

[TABLE]

Here $s$ is the index of sample (or data point), and $N_{s}$ is the total number of training samples. In Eq. (2.10) and (2.11) , the error function is normalized by the squared norm of $\bar{\textbf{C}}^{dns}_{s}$ to remove the scaling effect. The operator $||...||$ denotes the Frobenius matrix norm. The regularization term, scaled by a positive hyper-parameter $\lambda$ , is introduced to make the optimization problem well-posed. Although $\lambda$ will not alter the optimum fitting parameters, in practice, it should be set appropriately to expedite the gradient descent algorithm in the training process.

Design of experiments As mentioned in Section 2.1, the training and test datasets can be generated through high-fidelity DNS of an RVE (e.g. FEM and FFT-based methods). Since the input and output stiffness matrices for DMN can also be measured from experiments, it is also promising to incorporate both numerical and experimental data in the training dataset following the same format. When generating the input data, the components in the compliance matrices $\textbf{D}^{p1}$ and $\textbf{D}^{p2}$ are first sampled through design of experiments, and then inverted to give the stiffness matrices.

In general, an anisotropic compliance matrix should be used in order to uniquely determine the material orientation for the bottom-layer nodes. However, we consider both material phases to be orthotropically elastic to reduce the dimension of the sampling space and ease the data generation process. Since the material models used in the online stage have orthotropic symmetry, the non-uniqueness of orientations will not affect the predicted results. In Mandel notation, the components in $\textbf{D}^{p1}$ and $\textbf{D}^{p2}$ are

[TABLE]

To initiate material anisotropy, the tension moduli of phase $i$ are first randomly assigned as

[TABLE]

where $U$ stands for uniform distribution. The moduli of each material phase will be rescaled to create overall phase contrasts, and the rescaling factor $\bar{E}^{pi}$ are sampled as

[TABLE]

Note that the rescaling factor of phase 1 is kept as constant to remove the scaling effect of RVE. We can then update the tension moduli of each phase

[TABLE]

After the tension moduli are determined, the shear moduli are sampled as

[TABLE]

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

[TABLE]

DoE based on the Monte Carlo sampling is performed to explore the input space. For each material system, we generate 500 DNS samples. The first 400 samples are selected as the training dataset, and the remaining 100 samples become the test dataset. The size of the test set is selected to be large enough so that the distribution of test error is sampled adequately. Ideally, It is always better to increase the size of the training set for covering more possible cases, however, this will put more loads on the data generation process using DNS. In our study, the number of training samples is determined when the average test error (of a trained model) asymptotes sufficiently close to the average training error [46].

As an example, we demonstrate the distributions of training and test samples for the particle-reinforced RVE in Figure 4. The DoE procedures for all the material systems are same in terms of the material models and sampling bounds. As a result, their distribution plots generated by random sampling are similar. It can be seen from Figure 4 (a) that the ratios between tension moduli along different directions for a single phase range from 0.01 to 100. Meanwhile, the contrasts of tension moduli between the two phases can go up to 10000 as shown in Figure 4 (b). Due to the strong material anisotropy and high phase contrasts in the training dataset, it will be challenging for methods relying on micromechanics assumptions (e.g. analytical mean-field methods [15, 16], SCA [40]) to capture the behaviors of all the samples. By preserving the physics in the building blocks, a DMN can cover the training dataset accurately with sufficient number of layers. In Section 4, we will validate its capability for various material systems.

SGD with model compression. Since the building block is designed to have analytical solutions for $\boldsymbol{\mathfrak{h}}_{C}$ and $\boldsymbol{\mathfrak{r}}_{C}$ , gradient-based optimization methods can be utilized to update the fitting parameters and minimize the cost function, and the gradient vector of the cost function $\nabla J$ is computed by the back-propagation algorithm. To accelerate the training speed, we use the stochastic gradient descent (SGD) to train the network. Instead of computing and averaging the gradients over all the samples at each step, the original dataset is divided into several mini-batches which will be used in a sequence of learning steps. Each time the algorithm has processed all the mini-batches and seen all the samples in the original dataset (or one “epoch”), the dataset will be randomly shuffled to minimize the sample bias. All fitting parameters are initialized randomly following a uniform distribution at the beginning of the training,

[TABLE]

Similar to 2D DMN [1], the 3D DMN training process can be accelerated through two model compression methods: 1) Deletion of the parent node with only one child node; 2) subtree merging based on the similarity search. In similarity search, both volume fractions and rotation matrices are compared layer-by-layer between two subtree structures. The representation of 3D rotation using three angles may encounter the loss of one degree of freedom, known as the “Gimbal lock”, in which case two different sets of angles may represent the same 3D rotation. Therefore, we propose to use the rotation matrix R in the similarity search algorithm. Mathematically, two rotations, denoted by $\textbf{R}_{1}$ and $\textbf{R}_{2}$ , are said to be the same if

[TABLE]

However, due to the use of orthotropic material in the offline training data, the similarity condition between two arbitrary rotation matrices $\textbf{R}_{1}$ and $\textbf{R}_{2}$ at the bottom layer is relaxed: All the eigenvalues $\lambda_{i}$ of the matrix $(\textbf{R}_{1})^{-1}\textbf{R}_{2}$ should satisfy

[TABLE]

where $tol$ is the tolerance and $\mathbb{R}(\lambda_{i})$ denotes the real part of $\lambda_{i}$ . If two subtrees are found similar to each other, they will be merged to one single branch. Each time before the similarity search and subtree merging, the whole material network will be reordered based on the nodal weights. To save the training time, the model compression operations are performed every 10 epochs in our study.

2.4 Online prediction

The online stage is intended to extrapolate the trained DMN to unknown material and loading spaces. Each active node in the bottom layer is treated as an independent material point, with its own loading path and internal variables (if any). For a general finite-strain problem, we will choose the deformation gradient F and the first Piola-Kirchhoff stress P as the strain and stress measures (see Section 3.2 for more details). With both material and geometric nonlinearities taken into account, the mechanical data now contains two parts: The tangent stiffness matrix A and the residual stress $\delta\textbf{P}$ . As presented in Figure 5, the inputs of the building block are $\bar{\textbf{A}}_{i+1}^{2k-1},\bar{\textbf{A}}_{i+1}^{2k},\delta\bar{\textbf{P}}_{i+1}^{2k-1}$ and $\delta\bar{\textbf{P}}_{i+1}^{2k}$ from lower layers. Two homogenization functions $\boldsymbol{\mathcal{H}}_{A}$ and $\boldsymbol{\mathcal{H}}_{P}$ are introduced separately to map the inputs to a new stiffness matrix ${\textbf{A}}_{i}^{k}$ and residual stress $\delta{\textbf{P}}_{i}^{k}$ , respectively. Then, they are rotated by functions $\boldsymbol{\mathfrak{R}}_{A}$ and $\boldsymbol{\mathfrak{R}}_{P}$ , and the outputs of building block $\mathfrak{B}_{i}^{k}$ are $\bar{\textbf{A}}_{i}^{k}$ and $\delta\bar{\textbf{P}}_{i}^{k}$ .

The nonlinear RVE problem is solved via Newton’s method. Each Newton iteration contains one forward homogenization process and one backward de-homogenization process, and the mechanical data in the forward and backward propagations are shown as below,

[TABLE]

At the start of each iteration, the increment of deformation gradient $\Delta\boldsymbol{F}_{N}^{j}$ is prescribed on each active node in the bottom layer $N$ , and the values of $\bar{\textbf{A}}_{N}^{j}$ and $\delta\bar{\textbf{P}}_{N}^{j}$ are obtained by evaluating the local constitutive law. Then the data of stiffness matrix and residual stress are feed forward to the top layer through the homogenization and rotation functions at each building block. The overall $\bar{\textbf{A}}^{rve}$ and $\delta\bar{\textbf{P}}^{rve}$ of the RVE are obtained as the results of one forward homogenization process. Given the macroscopic loading condition, we can then compute the overall increments of deformation gradient $\Delta\bar{\textbf{F}}^{rve}$ and stress $\Delta\bar{\textbf{P}}^{rve}$ . In the de-homogenization process, the increments of deformation gradients and stress are feed backward from the top layer to the bottom layer, using the stress-strain relations at each building block. Finally, we return to the bottom layer with the new increment of deformation $\Delta\boldsymbol{F}_{N}^{j(new)}$ at each active node. The convergence condition is evaluated by checking the relative difference between $\Delta\boldsymbol{F}_{N}^{j}$ and $\Delta\boldsymbol{F}_{N}^{j(new)}$ . If convergence is not yet achieved, $\Delta\boldsymbol{F}_{N}^{j}$ will be updated by $\Delta\boldsymbol{F}_{N}^{j(new)}$ and used for the next iteration. Upon convergence, the internal variables at each active node in the bottom layer will be updated, and the analysis moves on to the next loading step.

An important feature of DMN is that the total number of operations in one Newton’s iteration is proportional to the number of active nodes in the bottom layer, and we will study the computational cost of DMN for various 3D RVE problems in Section 5. Another interesting feature of DMN is that different networks can be easily concatenated to simulate a material system with more than two scales, and the new integrated network can be solved in the online stage without altering the algorithms. In Section 4.3, as an example, we will demonstrate a three-scale nonlinear homogenization of CFRP system by attaching networks of the microscale unidirectional fiber composite to the bottom-layer nodes in a network of the mesoscale woven composite.

3 Mechanistic building block in 3D

3.1 Fundamentals for linear elasticity

The theory of 3D building block is first discussed under small-strain assumption, which is mainly used in the offline training stage. The stress and strain measures are the Cauchy stress $\utilde{\boldsymbol{\sigma}}$ and the infinitesimal strain $\utilde{\boldsymbol{\varepsilon}}$ , and they are related by a fourth-order stiffness tensor $\utilde{{\textbf{C}}}$ . The overall stress-strain relation of the building block can be expressed as

[TABLE]

For a building block with two materials 1 and 2, we have

[TABLE]

Here a bar is placed above each individual stiffness tensor (e.g. $\utilde{\bar{\textbf{C}}}^{1}$ ) to denote that it is from the averaging of a previous building block in the network. With the weights $w^{1}$ and $w^{2}$ given, the volume fractions $f_{1}$ and $f_{2}$ can be computed by

[TABLE]

In Mandel notation, the stress and strain can be written as

[TABLE]

and

[TABLE]

where the subscripts 3, 4 and 5 denote the shear directions, and the corresponding stiffness matrix C has 21 independent components.

Similar to the building block in 2D [1], a two-layer structure is proposed in this paper, as shown in Figure 6. Due to its simplicity, it is possible to derive the analytical forms of its averaging functions, which are essential for calculating the gradients during the network training process. There are two operations evolved in the averaging of building block: 1) Homogenization given by the function $\boldsymbol{\mathfrak{h}}_{C}$ ,

[TABLE]

and 2) rotation of the whole two-layer structure given by $\boldsymbol{\mathfrak{r}}_{C}$ ,

[TABLE]

Homogenization function $\boldsymbol{\mathfrak{h}}_{C}$ . The analytical form of homogenization function $\boldsymbol{\mathfrak{h}}_{C}$ is derived based on the interfacial equilibrium conditions

[TABLE]

and the kinematic constraints

[TABLE]

First, assume an arbitrary overall strain after homogenization

[TABLE]

Rewrite the kinematic constraints in Eq. (3.8),

[TABLE]

and apply the averaging equation for strain $\boldsymbol{\varepsilon}=f_{1}\boldsymbol{\varepsilon}^{1}+f_{2}\boldsymbol{\varepsilon}^{2}$ ,

[TABLE]

Utilizing the constitutive law for each material and substituting Eq. (3.10) and (3.11) into the equilibrium condition in Eq. (3.7) yield the following equation for the strain components of material 1,

[TABLE]

where

[TABLE]

After solving the linear system, the unknown strain components $\varepsilon_{3}^{1}$ , $\varepsilon_{4}^{1}$ and $\varepsilon_{5}^{1}$ in material 1 can be expressed as a function of the overall strain,

[TABLE]

with

[TABLE]

Combining Eq. (3.10) and (3.14) gives the definition of strain concentration tensor of material 1 $\textbf{s}^{1}$ ,

[TABLE]

Rewrite the averaging equation of the stress based on the definition of $\textbf{s}^{1}$ :

[TABLE]

Finally, the homogenized stiffness matrix before the rotation operation is obtained:

[TABLE]

Rotation function $\boldsymbol{\mathfrak{r}}_{C}$ . Here we use the Tait-Bryan angles $(\alpha,\beta,\gamma)$ to parameterize the 3D rotation. Any rotation can be achieved by composing three element rotations, in other words, a rotation matrix R can be decomposed as a product of three element rotation matrices,

[TABLE]

The expressions of elementary rotation matrices are given in A. The rotation function $\boldsymbol{\mathfrak{r}}_{C}$ takes the following form:

[TABLE]

In practice, the 3D rotation is decomposed into 3 rotation steps, characterized by $\alpha$ , $\beta$ and $\gamma$ ,

[TABLE]

More details on computing the gradients for training are provided in B.

3.2 Extension to material and geometric nonlinearities

In this section, we extend the solutions of two-layer building block from linear elasticity to problems with both material and geometric nonlinearities, as to appear in the online extrapolation stage. Note that the same procedure can be degenerated to problems with only material nonlinearity (e.g. small-strain plasticity) or geometric nonlinearity (e.g. hyperelasticity), so that we will only focus on the general situation.

For finite-strain problem, the deformation gradient $\utilde{\textbf{F}}$ and first Piola-Kirchhoff stress $\utilde{\textbf{P}}$ are chosen as the strain and stress measures, respectively. Here, we defined the tangent elasticity tensor $\utilde{\textbf{A}}$ in the rate form as

[TABLE]

The first elasticity tensor has major symmetry, $\utilde{A}_{ijkl}=\utilde{A}_{klij}$ , but does not have minor symmetries. As a result, Voigt notation can not be applied in this case, instead, we vectorize the deformation gradient and first Piola-Kirchhoff stress in 3D based on all the 9 components in each tensor:

[TABLE]

Under the above notations, the stiffness matrix A is defined as

[TABLE]

Due to major symmetry, A is symmetric, and it has 45 independent components. For a general loading step, the stress-strain relation takes the following form,

[TABLE]

where $\Delta\textbf{P}$ and $\Delta\textbf{F}$ are the increments of stress and deformation gradient. Due to the existence of nonlinearity, the residual first Piola-Kirchoff stress $\delta\textbf{P}$ needs to be included in the relation.

In the two-layer building block, the constitutive laws of materials 1 and 2 are

[TABLE]

Again, we place the bar notations above the stiffness matrix and residual stress of each material (e.g. $\bar{\textbf{A}}^{1}$ , $\delta\bar{\textbf{P}}^{1}$ ) to keep the fact that these quantities should come from the averaging of previous building blocks in the network. As illustrated in Figure 5, two steps are involved in the averaging of the building block with nonlinearities: 1) Homogenization given by the functions $\boldsymbol{\mathcal{H}}_{A}$ and $\boldsymbol{\mathcal{H}}_{P}$

[TABLE]

and 2) rotation of the homogenized two-layer structure given by the functions $\boldsymbol{\mathfrak{R}}_{A}$ and $\boldsymbol{\mathfrak{R}}_{P}$ ,

[TABLE]

The intermediate constitutive law right after homogenization can be written as

[TABLE]

and the overall constitutive law after rotation is

[TABLE]

Homogenization functions $\boldsymbol{\mathcal{H}}_{A}$ and $\boldsymbol{\mathcal{H}}_{P}$ . Similar to the linear elastic case discussed in Section 3.1, analytical forms of the homogenization functions $\boldsymbol{\mathcal{H}}_{A}$ and $\boldsymbol{\mathcal{H}}_{P}$ can be derived based on the equilibrium conditions

[TABLE]

and the kinematic constraints

[TABLE]

In the derivation of the stiffness tensor, we first consider the situation with no residual stress and solve for the unknown components of the increment of deformation gradient $\Delta\textbf{F}^{1}$ of material 1 in the 3, 4, 6 directions. Then the strain concentration tensor $\textbf{S}^{1}$ of material 1, defined by $\Delta\textbf{F}^{1}=\textbf{S}^{1}\Delta\textbf{F}$ , is found to be

[TABLE]

with

[TABLE]

where definitions of $f_{1}$ and $f_{2}$ are given in Eq. (3.3). The homogenization function $\boldsymbol{\mathcal{H}}_{A}$ of the stiffness matrix A is derived to be

[TABLE]

To derive the homogenization function $\boldsymbol{\mathcal{H}}_{P}$ for the residual stress $\delta\textbf{P}$ , we first set the average strain after the homogenization $\Delta\textbf{F}$ in Eq. (3.29) to zero:

[TABLE]

With the kinematic constraints defined in Eq. (3.32) and the averaging conditions, we have

[TABLE]

Substituting the above equations to the equilibrium condition defined in Eq. (3.31) yields

[TABLE]

After solving for $\Delta{\textbf{F}}^{1}$ in material 1, the homogenized residual stress is given by the average stress increment $\Delta\textbf{P}$ in the structure,

[TABLE]

Therefore, the homogenization function $\boldsymbol{\mathcal{H}}_{P}$ is derived to be

[TABLE]

Rotation functions $\boldsymbol{\mathfrak{R}}_{A}$ and $\boldsymbol{\mathfrak{R}}_{P}$ . The rotation matrices are also required to be modified under finite-strain notation in Eq. (3.23). Based on the Tait-Bryan angles $\alpha,\beta,\gamma$ , the new rotation matrix $\textbf{R}^{f}$ can be decomposed as a product of three element rotation matrices,

[TABLE]

The expressions of finite-strain element rotation matrices $\textbf{X}^{f}$ , $\textbf{Y}^{f}$ and $\textbf{Z}^{f}$ are given in A. With the definition of the rotation matrix $\textbf{R}^{f}$ , the rotation functions can be expressed as

[TABLE]

and

[TABLE]

4 Applications

In this paper, we explore the 3D architectures of DMN for three typical multiscale material systems which are popular in both academic research and industrial applications: 1) Particle-reinforced rubber composite with Mullins effect in Section 4.1; 2) polycrystalline materials in Section 4.2; 3) carbon-fiber reinforced polymer (CFRP) composites in Section 4.3. In particular, three scales are considered in the CFRP example, where the corresponding two RVEs are the microscale unidirectional fiber (UD) composite and the mesoscale woven composite. We will first show the DMN training results based upon orthotropic linear elasticity. Then we will extrapolate the trained networks to nonlinear materials, including Mooney-Rivlin hyperelasticity with Mullins effect, isotropic von Mises plasticity and rate-dependent crystal plasticity. In the SGD algorithm, the mini-batch size is chosen to be 20, and there are 20 learning steps in each epoch (for 400 training samples). For all the cases, the coefficient $\lambda$ is set to 0.001, so that the regularization term will not affect the training at early stage with large MSE, but it will eventually limit the magnitudes of activations when the overall cost $J$ becomes smaller. Here we use $N_{a}$ to denote the number of active nodes in the bottom layer, so that $N_{a}\leq 2^{N-1}$ . Relative errors are computed at the end of training to evaluate the accuracy of DMN. For a sample $s$ , its error $e_{s}$ is defined as

[TABLE]

The average error of a whole dataset with S samples is

[TABLE]

The average errors of training and test datasets are denoted by $\bar{e}^{tr}$ and $\bar{e}^{te}$ , respectively. By default, finite-strain formulation with geometric nonlinearity is adopted in the online stages of all the RVEs.

4.1 Hyperelastic polymer composite with Mullins effect

The first RVE example is a two-phase particle-reinforced polymer composite. The geometry of the RVE is presented in Figure 7, together with the uniaxial responses of material models that will be studied in the online stage. Phase index of the particle material is 1, and there are four identical spherical particles embedded in the matrix. Volume fraction of the particle phase is 22.6%. The DNS RVE is discretized by 84693 nodes and 59628 10-node tetrahedron finite elements.

4.1.1 Offline evaluation

Histories of average training and test errors are plotted in Figure 8 for network depth $N=4,5,6,8$ , and the errors at the end of training are also listed in Table 1. Networks with $N=4,$ 5 and $6$ are trained for 20000 epochs. While for $N=8$ , we observed that the errors still decrease rapidly after 20000 epochs of training, so we choose to restart the learning process by doubling the learning rate, and train the network for another 20000 epochs. As we can see from the figure, a network with $N=4$ is not sufficient to capture the RVE responses, and its average training error stopped decreasing around 7.6%. Good accuracy can be achieved by the deeper networks with $N=6$ and $N=8$ . Due to the use of mechanistic building block, the difference between the training and test errors is almost negligible, indicating no over-fitting issue. For $N=8$ , the average training error $\bar{e}^{tr}$ was reduced to 0.53% after 40000 epochs of training, and the maximum error among all the samples in the test dataset is only 2.41%. Note that the contrast of stiffness between the two material phases can reach up to 1000 in the offline sampling space, indicating that the trained DMN discovers a good reduced representation of the original 3D RVE microstructure.

To evaluate the accuracy and efficiency of DMN objectively, we also introduce another FE model with linear elements as a reference and measure its performance on the test dataset with 100 samples. As shown in Table 1 and Figure 8 (b), the average test error $\bar{e}^{te}$ of the linear FE model is 2.30%, and the maximum test error is 12.5%. If we use $\bar{e}^{te}$ as the measure, the linear FE model performs comparably to the DMN model with $N=6$ and $N_{a}=13$ . Moreover, it has only 11236 nodes, hence, is more efficient than the original DNS model. The computational time of the linear FE model for the particle-reinforced RVE will be given in Section 5.

Treemaps of DMN at the end of training are provided in Figure 9 to show the fraction of each active bottom-layer node within the whole network. It should be noted that no rotation information is reflected in these treemaps. At the beginning of training, all the nodes in the bottom layer are active, so that we have $N_{a}=8,32,128$ for $N=4,6,8$ , respectively. At the end of training, these numbers are reduced to $N_{a}=4,13,28$ due to model compression. Volume fractions predicted by DMN based on the weights can be found in Table 1 and Figure 9. For $N=8$ , the relative difference between $vf_{1}$ predicted by DMN and the reference one of the full-field DNS model is only 0.88%. Although the training dataset only contains information of the overall RVE stiffness matrices, the trained DMN is capable of extracting essential geometric information of the RVE accurately from these mechanical data. The treemaps also provide some physical insights on the RVE microstructure: the particle block (dark) is always surrounded by multiple matrix blocks (light), indicating that the matrix phase dominates the overall RVE properties.

4.1.2 Online extrapolation

We start from a Mooney-Rivlin hyperelastic rubber model. To model the rubber as an unconstrained material, a hydrostatic work term scaled by the bulk modulus $K$ is included in the strain energy functional which is function of the relative volume $J$ [47]:

[TABLE]

In order to prevent volumetric work from contributing to the hydrostatic work, the first and second invariants are modified,

[TABLE]

where $I_{1},I_{2},I_{3}$ are the invariants of right Cauchy-Green tensor $C=F^{T}F$ . The independent parameters in the hyperelastic material model are chosen as $C_{10},C_{01}$ and Poisson’s ratio $\nu$ . The shear and bulk moduli can be computed by

[TABLE]

Mullins effect, a kind of stress softening phenomena, is considered in the matrix rubber material. As the rubber material is loaded, the breakage of links between rubber chains happens at different deformation levels, thereby leading to continuous damage in the macroscale. Note that the Mullins effect will cause dispersed softening in the RVE, rather than localized damage or even discrete cracks. For RVE problems with localizations, the current homogenization scheme based on the first-order Hill-Mandel condition is questionable and needs to be revised [12], so does the corresponding DMN framework. To model the Mullins effect, the resulting stress is reduced by a damage factor according to

[TABLE]

where $W_{d}$ is the current value of the deviatoric strain energy density, and $W_{d}^{max}$ is the peak deviatoric strain energy density reached up to this point in time. The damage factor is defined as

[TABLE]

where $\text{erf}()$ denotes the Gauss error function, and $\eta,a,b$ are the empirical material constants related to the Mullins effect. In addition, the particle phase is modeled as a Neo-Hookean material which is 100 times harder than the matrix phase in terms of $C_{10}$ . All the material parameters are provided in Table 2.

The stress-strain curves of the particle-reinforced rubber composite under uniaxial tension and shear loadings are shown in Figure 10. For both cases, the RVEs are unloaded to check the Mullins effect due to matrix damage. It can been seen from the figure that material responses for all the cases can be well captured by the DMN with $N=8$ , which has 28 active nodes in the bottom layer. As the DNS model contains much more DOFs (84693 nodes), the computation based on DMN is much faster than the DNS. Detailed investigation on the computational cost are provided in Section 5. Since the active nodes in the bottom layer serve as the microscale material points, the distributions of local quantities (e.g. strain, stress and internal variables) can be recovered by DMN. More results on this subject can be found in [1]. It is also possible to relate the bottom-layer nodes directly to parts in the real RVE microstructure using clustering techniques, which will be investigated in our future work.

4.2 Polycrystalline materials with rate-dependent crystal plasticity

Another application of the 3D DMN is to evaluate the effective properties of polycrystalline materials, such as Ni-based super-alloy. We used the software package DREAM.3D [48] to generate the polycrystalline RVEs of equiaxed grains with two different crystallographic orientation distribution functions (ODF): 1) Random ODF and 2) textured ODF with preferred grain orientations. For both cases, the nominal number of grains inside the RVE is set to 415. As shown in Figure 11 (a), the polycrystalline RVE is represented in a $45\times 45\times 45$ mesh with 91125 material points. Since the mesh is uniform, FFT-based micromechanics method [21] is chosen for the full-field DNS to accelerate the computation. In the online stage, we will extrapolate trained polycrystalline DMNs to finite-strain rate-dependent crystal plasticity, and the uniaxial responses of 25 randomly selected grains from the RVE with random ODF are presented in Figure 11 (b).

4.2.1 Offline evaluation

The polycrystalline RVE can be treated as a single-phase material. In other words, all the grains shares the same material model, but with different orientations. Hence, the DMN in the offline stage takes only one input $\textbf{C}^{p1}$ as shown in Eq. (2.8), while other parts of the training process remain the same as before.

Figure 12 presents the histories of the training and test errors for two RVEs with random and textured ODFs. In general, the training process of polycrystalline RVE is smoother than the one of a two-phase material (see Figure 8), because it has less deactivation or compression operations which may cause abrupt temporary changes in the error. DMNs with three network depths $N=4,6,8$ are investigated, and their initial number of active nodes in the bottom layers are $N_{a}=8,32,128$ , respectively. We can see from the figures that the errors of $N=8$ can be reduced to less than 0.5% for both random and textured ODFs after 20000 epochs. Since the slopes are not saturated yet, further reduction of the errors could be achieved by continuing the training for more epochs.

Errors of the RVE with random ODF at the end of training are also provided in Table 3. For $N=6$ , the mean training and test errors are close to 1%, and the maximum test error is only $3.64\%$ . In practice, the choice of the network depth $N$ is a trade-off between the desired accuracy and computational cost. For example, if 1.27% of average error is satisfactory for a certain application, DMN with $N=6$ would be a better candidate for representing the polycrystalline RVE, since it will be around four times faster than the one with $N=8$ in terms of both offline training and online extrapolation.

Treemaps of DMN at the end of training are plotted in Figure 13 for the RVE with random ODF. Since there is only one phase in the RVE, the volume fractions cannot be defined, thus, are not shown in the plots. In each polycrystalline DMN, the number of active bottom-layer nodes did not decrease during the training, except for $N=6$ where only one node is deactivated. After 20000 epochs, the weights of bottom-layer nodes are still uniformly distributed. In some sense, this is consistent with the DNS RVE model with equiaxed grains.

As the weights in the network do not change much, the main training effort is devoted to optimizing the angles in each building block. Therefore, it is also interesting to check whether our DMN can learn the distribution of angles or ODF purely from the mechanical data of stiffness matrices. In this work, we follow the Bunge convention of Euler angles ( $\phi_{1},\Phi,\phi_{2}$ ) [49] and use the MTEX package [50] to generate the pole figures. Figure 14 collects the (001),(011), (111) pole figures from DNS model and predicted by DMNs with network depth ranging from 4 to 9. Note that only a quarter of the pole figure is shown due to the symmetry of orthotropic materials used in the offline stage. By looking at these plots, we can conclude that DMN with sufficient number of layers can discover the nature of random ODF hidden in the mechanical dataset. Finding the ODF for the RVE with textured ODF can be more challenging, and the pole figures for this case are listed Figure 15. It shows that for $N>4$ , the DMN is able to recover the dense regions appearing in the DNS pole figures, and the best match with DNS is achieved by the network with $N=9$ . Due to the existence of dense regions in the textured ODF, the compression operations (e.g. subtree merging) are more likely to be triggered during the training process, and the resulting DMN usually has less number of active bottom-layer nodes, or a smaller $N_{a}$ , than the one for random ODF with the same network depth.

4.2.2 Online extrapolation

In the online stage, we evaluate the trained DMN on a finite-strain rate-dependent crystal plasticity model [51, 52]. To begin with, the deformation gradient $\mathbf{F}$ is decomposed multiplicatively as:

[TABLE]

where the plastic part $\mathbf{F}^{p}$ maps points in the reference configuration onto an intermediate one which is then mapped to the current one through the elastic part $\mathbf{F}^{e}$ . The effect of dislocation motion is modeled by relating the plastic velocity gradient $\tilde{\mathbf{L}}^{p}$ in the intermediate configuration to simple shear deformation $\gamma^{(\alpha)}$ :

[TABLE]

where $N_{\text{slip}}$ is the number of slip systems. For a crystal slip systems $(\alpha)$ in the intermediate configuration, $\dot{\gamma}^{(\alpha)}$ is a shear rate, $\tilde{\mathbf{s}}^{(\alpha)}$ is the slip direction, and $\tilde{\mathbf{n}}^{(\alpha)}$ is the slip plane normal. The relationship between $\tilde{\mathbf{L}}^{p}$ and $\mathbf{F}^{p}$ is given by

[TABLE]

We choose to formulate the constitutive laws of elasto-plasticity based on the Green strain $\mathbf{E}^{e}$ and Second Piola-Kirchhoff stress $\mathbf{S}^{e}$ , which are related by:

[TABLE]

where the elastic stiffness tensor $\tilde{\mathbf{C}}$ is defined in the intermediate configuration. The plastic shear rate in each slip system given by a phenomenological power law

[TABLE]

where ${\tau}^{(\alpha)}$ is the resolved shear stress, $a^{(\alpha)}$ is a back-stress that describes kinematic hardening, $\dot{\gamma}_{0}$ is a reference shear rate, $\tau^{(\alpha)}_{0}$ is a reference shear stress that accounts for isotropic hardening, and the exponent $m$ determines the material strain rate sensitivity. The shear stress ${\tau}^{(\alpha)}$ is resolved onto the slip directions with:

[TABLE]

where $\bm{\sigma}$ , $\mathbf{s}^{(\alpha)}$ and $\mathbf{n}^{(\alpha)}$ are the Cauchy stress, slip direction and slip plane normal in the current configuration, respectively. Moreover, the evolution equation of the reference shear stress $\tau^{(\alpha)}_{0}$ is defined as

[TABLE]

where $H$ is a direct hardening coefficient and $R$ is a dynamic recovery coefficient and $q^{\alpha\beta}$ is the latent hardening ratio given by:

[TABLE]

where $\chi$ is a latent hardening parameter. The back-stress $a^{(\alpha)}$ in Eq. (4.12) evolves based on the following expression,

[TABLE]

where $h$ and $r$ are direct and dynamic hardening factors respectively. The set of nonlinear equations from Eq. (4.8) to Eq. (4.16) need to be solved numerically in the crystal plasticity model. In this work, we defined the model for a FCC metal (e.g. Ni-based super-alloy) with 12 slip planes. In the elastic regime, the material behaves isotropically. All the crystal plasticity parameters are given in Table 4.

To show the rate-dependency in the crystal plastic model, two strain rates were considered: $\dot{\varepsilon}=1.0\times 10^{-4}$ and $\dot{\varepsilon}=1.0$ . The polycrystalline RVEs with random and textured ODFs are simulated by both DNS and DMN under uniaxial loading up to $F_{11}=1.2$ , and the results are compared in Figure 16. The stress-strain curves of $N=4$ are shown as the dotted lines, and we can see from the plots that it overestimates the yield stress for random ODF while underestimates the value for textured ODF. Thus, it appears to be too ambitious to use 8 active bottom-layer nodes to represent the a RVE with more than 400 grains subject to different orientations. For $N=6$ and $N=8$ , the networks can predict the hardening behavior very well under different strain rates.

4.3 Carbon fiber-reinforced polymer with three scales

The last material system that we have investigated using DMN is the carbon fiber reinforced polymer (CFRP) composite. A typical example of CFRP system is illustrated in Figure 17 (a). The microscale RVE is a unidirectional (UD) fiber composite with straight fibers penetrating through the matrix in one direction. The UD RVE is then homogenized to provide properties for the yarn phase within the mesoscale RVE of a woven composite, and the macroscale overall material properties are obtained by homogenizing the woven RVE. Compared to the previous two-scale examples, the CFRP system contains three scales with an additional mesoscale RVE, which poses more challenges on the efficiency and robustness of a reduced order model.

The UD RVE is discretized by 69048 nodes and 63720 8-node hexahedron finite elements, and the volume fraction of the fiber phase is $vf_{1}=29.5\%$ . Additionally, the woven RVE has 43226 nodes and 28475 10-node tetrahedron finite elements, and the volume fraction of its yarn phase is $vf_{1}=41.4\%$ . The uniaxial responses of online material models for the carbon fiber and epoxy matrix are presented in Figure 17 (b).

4.3.1 Offline evaluation

In the offline stage, the microscale UD RVE and the mesoscale woven RVE are trained separately to find their own DMN representations. Histories of the average training and test errors of UD and woven DMN are given in Figure 18. For all the cases, we doubled the learning rate of SGD algorithm after 10000 epochs to accelerate the training speed, which also causes more oscillations in the error histories. According to our study, the bottom-layer nodes in a DMN of UD composite is more often deactivated than the other RVE examples during the training, mainly due to its strongly anisotropic microstructure. A network with $N=4$ may unintentionally lose all its nodes for the fiber phase. Therefore, we choose to show the DMN results of UD composite from $N=5$ to 9. After 20000 epochs of training, errors for $N=7$ and 9 both reduced to be less than 1%.

Finding the optimum DMN that well represents the woven RVE appears to be a more challenging task. The lowest training error we achieved is 2.19% with $N=8$ , and the training process tends to be saturated in the end. This is somehow expected since the morphology of woven microstructure is more complex than the other cases. For example, it has penetrating phases along two different directions, as well as the entanglement of yarns. The current level of accuracy achieved by DMN is satisfactory in many practical applications, however, further reduction of errors may require more advanced learning techniques.

Table 5 summarizes the training and test errors for UD and woven RVEs, as well as the volume fractions of phase 1 predicted by the corresponding DMNs with various depths. A linear FE model with 5252 nodes is introduced for the woven composite as a reference, and its accuracy is evaluated on the test dataset. The average test error $\bar{e}^{te}$ and maximum test error of the linear FE model are 6.12% and 28.3%, respectively. Measured by $\bar{e}^{te}$ , its performance is close to the DMN model with $N=7$ and $N_{a}=16$ . The computational time of the linear FE model for woven composite will be given in Section 5. Treemaps are provided in Figure 19. Again, we can conclude that for both UD and woven composites, a DMN with sufficient number of layers is able to extract the volume fraction accurately from the training dataset with pure mechanical information.

4.3.2 Online extrapolation

Ideally, a full three-scale homogenization should be performed by linking every material point in the mesoscale woven RVE with a microscale UD RVE model. However, this is not feasible with the DNS models described above due to the high computational cost. On the other hand, solving the woven RVE implicitly also requires the full information of the overall stiffness tensor of the UD RVE, while the standard FEM solver only provides the values of stress and strain. Therefore, the DNS models of UD and woven RVEs are performed separately, and the obtained two-scale DNS homogenization results will be used to validate our DMN models.

For the UD composite, the fiber phase is considered as an orthotropic elastic material with high strength in the fiber direction. The matrix phase is modeled as an elasto-plastic epoxy material with an exponential hardening law. Its yield stress $\sigma_{Y}$ is determined by the hardening law as a function of the effective plastic strain $\bar{\varepsilon}^{pl}$ , which is a monotonically increasing internal state variable of the plastic material during the deformation. An exponential hardening law is assumed,

[TABLE]

where $a_{3}$ is the ultimate yield stress for large effective plastic strain, $(a_{3}-a_{2})$ represents the yield strength and $a_{1}$ is a dimensionless hardening constant. An isotropic J2 yield surface is further assumed,

[TABLE]

where $\bar{\sigma}$ is the von Mises equivalent stress. We choose to formulate the plastic law based on the Green-Lagrange strain and second Piola-Kirchhoff stress. Uniaxial responses of the fiber and epoxy materials are shown in Figure 17 (b). In the two-scale homogenization of woven RVE, the matrix phase shares the same epoxy material as the UD composite. The yarn phase is assumed to be an orthotropic elastic material, and its elastic constants are determined by homogenizing the UD composite. All the material parameters used in the online stage can be found in Table 6 111The hardening parameters in the published version were found to be incorrect. The parameters corresponding to the numerical examples are updated..

The stress-strain curves predicted by DNS and DMN for the UD composite are presented in Figure 20. Four loading-unloading conditions were considered: 11, 33 for uniaxial tension and 12, 23 for shear. We can see that the material responses in 33 direction are almost elastic, since the carbon fibers go straightly through the matrix and its Young’s modulus ( $E_{3}=245.0$ GPa) in this direction is 65 times larger than the one of epoxy matrix ( $E_{m}=3.8$ GPa). On the other hand, the overall stress-strain responses in other loading directions are dominated by the hardening of the epoxy matrix. For all these cases, the loading-unloading behaviors can be well captured by DMNs with $N\geq 7$ .

The stress-strain curves predicted by DNS and DMN for the woven composite are presented in Figure 21. Similarly, four loading-unloading conditions were considered. Figure 21 (a) shows the in-plane tension results in 11 direction. The overall stress-strain curve first slightly yields due to matrix hardening, and more interestingly, a stiffening effect is observed after the strain $\varepsilon_{11}$ reaches 1.5%. As illustrated in Figure 22 (a), the anisotropic yarn is straightened during the deformation, making it stiffer in the loading direction. Although the overall deformation is small ( $<2.0\%$ ), some segments within the yarn experience relatively large local rotations, which can be captured only if a model takes the geometric nonlinearity into account. To validate this hypothesis, we also loaded the DMN model for $N=8$ based on the small-strain formulation without geometric nonlinearity, and as expected, the stiffening effect was not reproduced as in the DNS or DMN based on finite-strain formulation (see Figure 22 (b)).

For the out-of-plane tension in 33 direction, the overall responses still show little hardening effect due to the yarn confinements in 11 and 22 directions. By looking at Figure 21 (d), we can also conclude that the out-of-plane shear in 23 direction is dominated by the epoxy matrix since the overall material is almost perfectly plastic in the end. On the other hand, the in-plane shear behavior in 12 direction is influenced more evenly by the yarn and matrix phases, as the DNS stress-strain curve shown in 21 (c) has a nontrivial hardening modulus. For all the cases, the predictions from DMN with $N=8$ match the DNS results very well, indicating that our reduced order model preserves the essential microstructural interactions in the original DNS model. One can also see that it will be very challenging to manually develop an empirical material model to accurately capture the complex responses of the woven composite. However, with the help of machine learning techniques, DMN can automatically find the optimum reduced description of the woven composite by learning the hidden physics in the training dataset.

Finally, we demonstrate a three-scale homogenization procedure for CFRP composites based on DMN. As shown in Figure 23, it can be realized through the concatenation of networks: Link every active bottom-layer node of the yarn phase in the woven composite with a copy of DMN for the UD composite. As a result, the online calculation of the new integrated network only requires constitutive models for the matrix and fiber materials. Here we take the concatenation of the UD DMN with $N^{ud}=7$ and woven DMN with $N^{woven}=8$ as an example, and we count each active node with an independent material input as a DOF. For the two-scale woven DMN, the number of DOFs is

[TABLE]

For the three-scale woven DMN, every yarn DOF is replaced by a UD DMN, and the new number of DOFs can be computed as

[TABLE]

In the online stage for the three-scale homogenization, we conduct the Newton’s method separately for the UD and woven parts in the integrated network. At each newton iteration of the woven DMN, it takes the converged stiffness matrices and residual stresses from the newton iterations of the copies of the UD DMN. In this way, all the algorithms developed for the two-scale homogenization do not need to be altered. However, for the best efficiency, the three-scale network shown in Figure 23 should be solved as a whole in one Newton iteration, so that no extra time is paid to the sub-cyclings in the UD structures.

The in-plane uniaxial tension and shear results predicted by the two-scale and three-scale woven DMNs are compared in Figure 24. Due to the consideration of matrix plasticity in the yarn phase in the three-scale model, more hardening effects are observed in the overall material responses for both cases. After the initial yielding, the three-scale in-plane tension test also shows a stiffening effect due to the straightening of fibers. For the three-scale in-plane shear test, the fully developed hardening behavior is chose to perfectly plasticity, while this phenomena is not captured by the two-scale model due to the lack of plasticity in the homogenized UD material. In summary, more physics within the CFRP system can be considered by the three-scale model based on DMN, while the one based on DNS is computationally infeasible.

5 A comment on computational cost

Last but not least, we shall discuss the computational cost of DMN in both online and offline stages. The computational time of generating the dataset and training the network with various depths for the particle-reinforced composite (see Section 4.1) are listed in Table 7. All the processes are parallelized using 10 Intel® Xeon® CPU E5-2640 v4 2.40 GHz processors. The generation of the training dataset with 400 samples took 39.5 h, or 356 s per sample. Note that to obtain all the components of the overall elastic stiffness matrix for each sample, we need to analyze the same RVE under 6 orthogonal loading conditions. Therefore, one RVE analysis of the particle-reinforced RVE takes around 59 s. Since the generation of each sample is independent, the wall time for generating the training data decreases linearly with the number of processors.

The DMN training times for $N=4$ , 6 and 8 are 5.4, 16.7 and 43.0 h, respectively. The time cost per epoch decreases as the number of active nodes in the bottom layer is reduced by the model compression algorithm. Here we take the network with $N=8$ as an example. The first 1000 epochs at the beginning of the training took 3.7 h, while the last 1000 epochs only cost 1.9 h. Our current parallelization algorithm is realized by distributing the training samples between various processors, or called “data level parallelism”. It has a good scalability as long as the number of processors is smaller than the mini-batch size in the SGD algorithm. Future improvements will include vectorizing the operations in the building block and parallelizing the computation at the model level, which will allow us to accelerate the training procedure of large-scale DMNs using more processing units (e.g. GPU).

The computational times of DNS and the online DMNs for the hyperelastic particle-reinforced RVE are given in Figure 25 (a). By using 10 cores, the DNS with 59628 10-node tetrahedron elements required 4860 s. All the tests based on DMN were performed on one core. A linear relationship (slope $\approx 1.14$ ) between the computational time and the number of DOFs $N_{a}$ is observed for network depth ranging from 4 to 9. When $N=8$ (with $N_{a}=28$ ), the online stage of DMN took 6.0 s, which is 810 times faster than the DNS in terms of the wall time, or 8100 times faster if measured by CPU time. Another plot on the computational times for the woven RVE is provided in Figure 25 (b). The DNS was limited to two-scale homogenization where the yarn phase is considered to be elastic, and it cost about 2440 s on one core. Both two-scale and three-scale homogenizations can be realized by DMN. As we can see from the figure, the two-scale DMN with $N^{\text{woven}}=8$ took 1.7 s, and the three-scale DMN with $N^{\text{woven}}=8$ and $N^{\text{woven}}=7$ took 39.4 s. Ideally, if the concatenated three-scale network is solved simultaneously without sub-cycling as discussed in Section 4.3, its computational time is estimated to be 11.0 s based on the numbers of DOFs for the two-scale and three-scale models given in Eq. (4.19) and (4.20).

It is more objective to evaluate the efficiency of DMN by comparing with a FEM model with similar accuracy. For the simulation of hyperelastic particle-reinforced composite, the linear FE model took 2420 s on one CPU. As shown in Section 4.1, its average test error on 100 random samples (2.30%) is close to that of the DMN with $N=6$ (1.38%), which only took 2.6 s on one CPU for the same number of loading steps. At the similar level of accuracy, the DMN is around 930 times faster than the linear FE model. For the two-scale homogenization of woven composite, the linear FE model is comparable to the DMN with $N=7$ in terms of the average test errors, 6.12% vs. 3.48%, while their computational times are 212 s vs. 0.65 s. Due to the complexity of woven microstructure, it is difficult to further coarsen the linear FE mesh with 5252 nodes. However, the material network with $N=7$ for woven composite has only 16 active nodes in the bottom layer, which may correspond to a $4\times 2\times 2$ FE mesh. Note that the current DMN framework is implemented as Python libraries, and it is possible to improve the efficiency through lower-level programming.

6 Conclusions

A new framework of deep material network (DMN) for modeling general 3D heterogeneous materials with both material and geometric nonlinearities is developed. Theories of the 3D two-layer building block are derived based on the homogenization equations and interfacial conditions. We show a complete machine learning procedure for training the DMN by first generating offline elastic datasets using DNS, and then optimizing the network fitting parameters with stochastic gradient descent and model compression algorithms. The 3D architectures of DMN are explored for three representative multiscale materials systems, including particle-reinforced hyperelastic rubber composite with Mullins effect, polycrystalline materials with finite-strain rate-dependent crystal plasticity, and elasto-plastic CFRP systems formed by UD and woven composites. For all the examples, we demonstrate the effectiveness of DMN on representing complex RVE microstructures and its accuracy of predicting highly nonlinear behaviors in the online stage. One unique feature of DMN is that it can discover the hidden geometric information, such as phase volume fraction and orientation distributions, from pure mechanical data. While DNS is limited to two-scale homogenization due to its high computational cost, we also show a three-scale DMN of the CFRP system through the concatenation of UD and woven networks, which captures more physics across multiple length scales. The proposed DMN and its machine learning procedure open new possibilities in the design and concurrent simulation of multiscale material systems. The idea of embedding physics into machine-learning models shows promise for the future of data-driven multiscale mechanics and physics.

Acknowledgment

Z. Liu and C.T. Wu warmly thank Dr. John O. Hallquist of LSTC for his support to this research. Z. Liu would like to thank Cheng Yu for his part in helpful discussions. This is a preprint version of an article published in Journal of the Mechanics and Physics of Solids. The final authenticated version is available online at: https://doi.org/10.1016/j.jmps.2019.03.004.

Appendix A Elementary rotation matrices

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

[TABLE]

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

[TABLE]

In finite-strain formulation, the elementary rotation matrices shown in Eq. (3.41 are given by

[TABLE]

where the in-plane and output-plane rotation matrices are

[TABLE]

for an arbitrary angle $\theta$ .

Appendix B Gradients for training

Derivatives need to be derived for computing the gradient vector of the cost function in the training process. Since the procedure of deriving the derivatives is standard based on the chain rule, we will only provide some essential steps here. All the derivatives of the building block that will appear during the training process are listed as below,

[TABLE]

The derivative of $(\hat{\textbf{C}}_{345})^{-1}$ with respect to an arbitrary variable $x$ (e.g. volume fraction or components in the stiffness matrix) will be necessary for computing the gradient of the strain concentration tensor $\textbf{s}^{1}$ in the homogenization function $\boldsymbol{\mathfrak{h}}_{C}$ :

[TABLE]

Moreover, the derivatives of volume fraction $f_{1}$ with respect to the weights $w^{1}$ and $w^{2}$ are

[TABLE]

Direct computation of derivatives of the rotation function $\boldsymbol{\mathfrak{r}}_{C}$ could be complicated. Instead, it is decomposed into three steps by utilizing Eq. (3.21), so that the derivatives of stiffness matrices with respect to the corresponding rotation angles are

[TABLE]

Specifically, $\textbf{X}^{\prime}(\alpha)$ denotes the derivative of the rotation matrix X to its rotation angle $\alpha$ . Furthermore, the derivatives of $\bar{\textbf{C}}$ with respect to the components in ${\textbf{C}}^{\alpha\beta}$ can be written as

[TABLE]

Similarly, we have

[TABLE]

Bibliography53

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Z. Liu, C. T. Wu, M. Koishi, A deep material network for multiscale topology learning and accelerated nonlinear modeling of heterogeneous materials, Computer Methods in Applied Mechanics and Engineering 345 (2019) 1138–1168.
2[2] I. M. Daniel, O. Ishai, I. M. Daniel, I. Daniel, Engineering mechanics of composite materials, Vol. 3, Oxford university press New York, 1994.
3[3] X. Deng, A. Korobenko, J. Yan, Y. Bazilevs, Isogeometric analysis of continuum damage in rotation-free composite shells, Computer Methods in Applied Mechanics and Engineering 284 (2015) 349–372.
4[4] B. Ren, C. T. Wu, P. Seleson, D. Zeng, D. Lyu, A peridynamic failure analysis of fiber-reinforced composite laminates using finite element discontinuous galerkin approximations, International Journal of Fracture 214 (1) (2018) 49–68.
5[5] A. C. Balazs, T. Emrick, T. P. Russell, Nanoparticle polymer composites: where two small worlds meet, Science 314 (5802) (2006) 1107–1110.
6[6] T. Ramanathan, A. Abdala, S. Stankovich, D. Dikin, M. Herrera-Alonso, R. Piner, D. Adamson, H. Schniepp, X. Chen, R. Ruoff, et al., Functionalized graphene sheets for polymer nanocomposites, Nature nanotechnology 3 (6) (2008) 327–331.
7[7] H. Xu, R. Liu, A. Choudhary, W. Chen, A machine learning-based design representation method for designing heterogeneous microstructures, Journal of Mechanical Design 137 (5) (2015) 051403.
8[8] Z. Liu, J. A. Moore, W. K. Liu, An extended micromechanics method for probing interphase properties in polymer nanocomposites, Journal of the Mechanics and Physics of Solids 95 (2016) 663–680.