{'fraction_non_alphanumeric': 0.035628588568000684, 'fraction_numerical': 0.012790299083040534, 'mean_word_length': 4.807041552625031, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Constitutive Relations using Symmetric Positive Definite Neural Networks

We present the Cholesky-factored symmetric positive definite neural network (SPD-NN) for modeling constitutive relations in dynamical equations. Instead of directly predicting the stress, the SPD-NN trains a neural network to predict the Cholesky factor of a tangent stiffness matrix, based on which the stress is calculated in the incremental form. As a result of the special structure, SPD-NN weakly imposes convexity on the strain energy function, satisfies time consistency for path-dependent materials, and therefore improves numerical stability, especially when the SPD-NN is used in finite element simulations. Depending on the types of available data, we propose two training methods, namely direct training for strain and stress pairs and indirect training for loads and displacement pairs. We demonstrate the effectiveness of SPD-NN on hyperelastic, elasto-plastic, and multiscale fiber-reinforced plate problems from solid mechanics. The generality and robustness of the SPD-NN make it a promising tool for a wide range of constitutive modeling applications.


Introduction
Material modeling aims to construct constitutive models to describe the relationship between strain and stress, in which the relationship may be hysteresis. The constitutive relations can be derived from microscopic interactions between multiscale structures or between atoms. However, the first-principles simulations, which resolve all these interactions, remain prohibitively expensive. Because of the computational difficulty, constitutive models are traditionally constructed with simplified assumptions and empirical relations. Parameters in the models are then calibrated on limited and coarse-scale tensile test data. These constitutive models lead to affordable simulations, and thus the models are very important for large-scale engineering and scientific applications. However, constructing these models is challenging because the models are usually mappings between high dimensional spaces and are different on a case by case basis.
Since the development of deep learning techniques, the deep neural network (NN) emerges as a promising alternative for constitutive modeling. For example, pioneering work has demonstrated the feasibility of constitutive modeling using neural networks in a wide variety of applications, such as Ghaboussi et al. [1] for modeling concrete, Ellis et al. [2] for modeling sands, Shen et al. [3] and Liang et al. [4] for modeling hyper-elastic materials, and Furukawa et. al. [5] for modeling viscoelastic materials. Recently, recurrent neural networks (RNNs), effective for history-dependent phenomena, have been applied to model multiscale multi-permeability poroplasticity [6], multiscaleplasticity [7], and multiscale one-dimenensional bars [8]. Additionally, neural networks are used to address more complex material behaviors, such as microcracking, brittle fracture, and crack propagation [9,10,11,12]. Moreover, attempts for new material designs are made by leveraging neural networks with microscopic structure parameters incorporated in the inputs [13,14,15,16]. The cited literature shows the potential capability of neural networks to represent complex constitutive relations. In this paper, we are more interested in using the NN-based constitutive model to predict material behavior, which involves embedding the neural network into finite element solvers.
When solving the conservation equations embedded with NN-based constitutive models that have no constraints on the neural networks, numerical instabilities are observed (see Section 4.1 for an example). We think one reason is the violation of the basic constraints related to these constitutive models. The idea of adding physics-related constraints to the neural network is not new. For example, Ling et al. [17] enforced isotropicity and cubic symmetry on the predicted strain-energy of the crystal elastic materials through building a basis of invariant inputs. Liu et al. [18,19] designed deep material networks, as composition of simple building blocks inspired from the two-phase linearly elastic model. Heider et al. [20] introduced coordinate-free invariant metrics for the objective function (the so-called loss function) to model anisotropic elastoplastic materials. However, despite the improved accuracy on the validation set by incorporating these constraints, the performances of the resulting conservation equations equipped with these physics constrained neural networks are not studied.
In the present work, we focus on the numerical stability aspect of the resulting hybrid modelthe conservation equations embedded with NN-based constitutive models. We propose a novel neural network architecture with customized output layers to fulfill two objectives.
• The first objective is related the convexity of the strain-energy, which plays significant role in the numerical stability. To this end, we propose to predict the tangent stiffness matrix and enforce it to be symmetric positive definite. The symmetric positive definiteness guarantees the weak convexity of the strain-energy.
• Another objective is the time consistency, related to path-dependent materials, such as elastoplastic material.
To achieve these objectives, instead of training a neural network to identify a nonlinear map directly between strain and stress, we train a neural network (with weights and biases θ) that maps the strain and other related quantities to a lower triangular matrix L θ (Cholesky factor). We construct the constitutive model in the following incremental form: where ∆ and ∆σ are the incremental strain and stress in Voigt notation 1 . In our formulation, the tangent stiffness matrix L θ L T θ is automatically symmetric positive semidefinite and the incremental form leads to time consistency.
In the present work, results based on neural networks trained on both direct input-output data and indirect full-field data [21,22,23] are presented. The robustness of our approach is demonstrated in different numerical applications, including a hyperelastic material, an elasto-plastic material, and a multiscale material. We have developed a software library that seamlessly combines traditional finite element methods and neural networks. The code is accessible online: https://github.com/kailaix/NNFEM.jl The remainder of this paper is organized as follows. We first introduce the background in Section 2, including the governing equations, different classical constitutive relations and their associated constraints. In Section 3, we present our constraint-embedded neural network architecture-SPD-NN, and the training procedures. Finally, we apply the learned NN-based constitutive relations for several solid mechanics problems, including a one-dimensional truss coupon, and several two-dimensional thin plate problems with hyperelastic material, elasto-plastic material and fiberreinforced multiscale material. Several issues related to neural networks is discussed in Section 5. We discuss a possible generalization of the approach in Section 6.

Governing Equations
The governing equation of a solid undergoing infinitesimal deformations can be written as where ρ is the mass density, u is the displacement vector, σ is the stress tensor, and ρb is the body force vector; Ω denotes the computational domain. The prescribed displacementū and the surface tractiont are imposed on the domain boundaries Γ u and Γ t with the outward unit normal n, where Γ u ∩ Γ t = ∅ and Γ u ∪ Γ t = ∂Ω.
To solve for the displacement u from Eq. (2), we also need the constitutive relations, which maps the deformation history of the structure to the stress: 1 For brevity, and σ represent both tensor form and vector form in Voigt notation depending on the context.
Here (t) is the strain tensor at time t related to the displacement vector and I(t) denotes all other quantities related to material states during time τ = [0, t), such as (τ ), σ(τ ), etc. The infinitesimal strain tensor is When we apply the finite element method to solve Eq. (2) numerically, we have the following semi-discrete equation at time t where we use the same notation u to denote the spatial discretization of the displacement vector u in Eq. (2), M is the discrete mass matrix, P and f are the discrete internal and external force vectors, x is the coordinate vector, and p is the parameter vector of external loads. We adopt the generalized α-method [24] with α m = −1 and α f = 0 for temporal discretization of Eq. (4). This generalized α-method allows for dissipating high frequency energy to damp high frequency modes, which is crucial for the robustness of the numerical solver when an approximate constitutive relation is used.
Remark 1. For structure undergoing finite or large deformations, the finite strain tensor reads σ in Eq.
(2) represents the first Piola-Kirchhoff stress tensor, the constitutive relation generally relates the finite strain tensor with the symmetric second Piola-Kirchhoff stress tensor The first Piola-Kirchhoff tensor and the second Piola-Kirchhoff tensor are related by

Constitutive Relations and Associated Constraints
Equation (3) represents one possible form of constitutive relations. Even when restricted to this specific form, the constitutive relations can describe a great variety of material properties. Nevertheless, most of the constitutive relations share extra constraints (i.e., objectivity, convex strain energy function, and second law of thermodynamics) that play significant roles in material stability and numerical stability. In what follows, several stability-and consistency-related constraints of constitutive relations are discussed. These constraints should be considered in any data-driven constitutive models.

Hyperelastic Materials
The constitutive relations of hyperelastic materials are path-independent and are related to the strain-energy density function ω( ), In general, the strain-energy density function is assumed to be convex, namely, the tangent stiffness matrix ∂ 2 ω( ) is symmetric positive definite (SPD). The assumption in one-dimensional is equivalent to that the strain-stress is monotonically increasing 2 . The convexity is crucial to both the stability of the material and also the numerical scheme.

Elasto-plastic Materials
The constitutive relations of elasto-plastic materials are rate-independent but path-dependent, and feature transition between elastic and plastic behaviors. Namely, under loading, the material behaves elastically until the initial yield stress σ Y is attained, and then undergoes permanent irreversible plastic deformations with further loading. The onset and continuance of plastic deformation is governed by a yield function For plastic deformation, the stress state must remain on the yield surface f = 0; and for elastic deformation, the yield function satisfies f < 0. The strain is assumed to be additively decomposed into elastic e and plastic p parts, as follows, The constitutive relation relates the stress and the elastic part of the strain: here C is the tangent stiffness tensor. The plastic strain rate is given by a flow rule, i.e., the associative flow rule as follows,˙ p =λ ∂f ∂σ whereλ is called the plastic rate parameter or the consistency parameter, which is non-zero only if f = 0. Naturally, the constitutive relation is written in the rate forṁ It is worth mentioning the tangent stiffness matrix H is symmetric positive semidefinite, and when strain hardening is considered, it becomes SPD. Moreover, the form in Eq. (11) implies time consistency, i.e., as ∆ → 0, ∆σ → 0, and also rate-independent, i.e., the tangent stiffness matrix is independent of the strain or stress rate.

Associated Constraints
Based on the discussion above, the following properties are crucial to be incorporated in datadriven constitutive models (1) Symmetry positive definiteness of the tangent stiffness matrix (i.e., strain energy convexity), which leads to non-singular stiff matrix. We also found that the SPD property is crucial for numerical stability.
(2) Time consistency, which is formulated as follows, It is crucial for the convergence of the numerical approximation when ∆t − → 0.

Methodology
In this section, we describe our method for learning the constitutive relations Eqs. (3) and (6). Our discussion will be divided into two parts: 1. the neural network architecture for approximating the mapping between the strain and the stress; 2. the direct and indirect training methods based on the types of available data.

Neural Network Architectures
We note that the neural network has already been used to approximate the constitutive relations. The neural network architectures in many literatures output stress directly (σ-NN) with the following form (recall that I stands for other relevant information up to the current time) or output stress increment directly (∆σ-NN): Those architectures are suitable for learning the constitutive relations when the strain and the stress data are available. However, these strain-stress relations expressed by Equations (13) and (14) do not satisfy certain physical constraints and thus may break numerical solvers (see Section 4.1), when we plug them into a numerical solver.
We propose an alternative architecture (SPD-NN) based on the incremental form, Here instead of outputting the stress or stress increment directly, the neural network outputs a lower triangular matrix L θ , which is the Cholesky factor of the tangent stiffness matrix H θ . The numerical approximation to Eq. (15) in the dynamic simulations has the following form here, the superscript n indicates the time step. An obvious advantage of Eq. (15) is the guarantee that the tangent stiffness matrix is a symmetric semidefinite matrix, which is true for commonly used constitutive relations. Additionally, when NN θ is bounded, we have lim which indicts the time-consistency for both path-dependent/independent constitutive relations. In contrast, both the σ-NN with σ n+1 = NN θ ( n+1 , n , σ n ) and the ∆σ-NN with σ n+1 = NN θ ( n+1 , n , σ n )+ σ n fail to satisfy the time-consistency condition Eq. (17) for path-dependent constitutive relations. Finally, when NN θ does not depend on strain rate, Eq. (15) implies that the constitutive law is rate-independent, which is a feature of the hyperelastic or elasto-plastic materials.
In the following, we discuss how to adapt SPD-NN for different materials.

• Linear Elasticity
For linear elastic material, the tangent stiffness matrix is independent of the strain or stress. Therefore, no neural network is needed, and the constitutive relation reads (18) where C θ is the parametric tangent stiffness tensor and the unknowns are simply entries in the tensor (no neural network).

• Nonlinear Elasticity
For nonlinear elastic material, the tangent stiffness matrix Eq. (7) depends only on the strain at the current time step. The constitutive relation Eq. (16) can be formulated as Note the input n+1 must be evaluated at the current time step n + 1 because it is impossible to determine the stress at time step n + 1 given only the information at time step n.

• Elasto-Plasticity
For elasto-plastic material, the material behavior features transition from elastic behavior to plastic behavior. To model this effect, we consider two types of constitutive relations: 1. In the linear elasticity region, the constitutive relation is approximated by 2. In the plasticity region, the constitutive relation is approximated by However, since we do not know when the transition occurs (the yield strength σ Y is not available and strain hardening could strength the material), we can relax the constitutive relation using a differentiable 3 transition function D(σ n ,σ Y ) (see Fig. 1), whose value is between 0 and 1, as follows, hereσ Y is the estimated yield strength, which does not need to be accurate. When the equivalent stress estimated from σ n (such as von Mises stress) 4 is smaller thanσ Y , the material behavior is assumed to be linear and described by Eq. (20); otherwise, it is assumed to be described by the plastic form Eq. (21). It is worth noting the plastic form Eq. (21) can degenerate to the linear elastic form Eq. (20), but not vice versa. Therefore, the estimated yield strengthσ Y should be smaller than the yield strength σ Y . In the present study, historical data for plastic form Eq. (21) span only one time step, more historical data might be needed for materials with strong hysteresis. Remark 2. Generally, we can write the strain-stress relations for a linear elastic material in Voigt notation, as follows: The corresponding Cholesky factor L θ of C θ is a full lower triangular matrix. For the materials considered in the present work, they are orthotropic, i.e., they have three mutually orthogonal planes of reflection symmetry. Thus the tangent stiffness matrix can be expressed by a block diagonal matrix (Orth-NN) Therefore, the associated Cholesky factor L θ has the following form Consequently, we can further simplify the neural network outputs to only the nonzero entries in Eq. (23). Equation (23) is the form of the Cholesky factor used in the present work.
Remark 3. Another restriction on plastically stable materials in small strains is given by the Drucker's postulate, which is related to the second law of thermodynamics. The Drucker's postulate can be stated in different ways, one of which is that the second-order plastic work is non-negative [26, p. 328] 0 ≤σ T˙ p (24) Bringing Eq. (10) and the rate form of Eq. (15) into Eq. (24) leads to The Drucker's postulate requires H θ C θ which brings an additional constraint to the tangent stiffness matrix besides SPD. We have attempted to design another neural network architecture imposing this constraint, However, this architecture does not improve the prediction and brings additional challenges in training the neural networks based on our observation. Therefore, this architecture Eq. (25) is not used in the present work.

Training Methods
Based on the types of available data, i.e., direct and indirect, we can employ different training methods. For the following discussion, we summarize the neural network based constitutive relations as follows plasticity Elasto-Plasticity

Direct Data
Direct data consist of the input and the output of the NN-based constitutive relations such as strain-stress pairs or strain-stress increments pairs. These data points come from experimental measurements and numerical simulation results. The comprehensive strain-stress data measurement relying on simple mechanical tests, such as tensile or bending tests, might be challenging. However, comprehensive strain-stress data generated from sub-scale simulations, such as representative volume element (RVE) simulations [27,28,29,30,31] or post-processed from direct numerical simulations, are widely used to train neural networks [13,14,17,6].
Mathematically, the direct data are given in terms of N sequences of strain-stress pairs at n time snapshots. ( . . , N Here the superscripts indicate time and the subscripts indicate the sequential number. We train the neural network by solving a minimization problem arg min

Indirect Data
Indirect data consist of deformation data from structure coupons under different load conditions. Deformations are measured by techniques such as digital image correlation or grid method [32]. These techniques can record complete heterogeneous fields, which are rich in the constitutive relations. The virtual fields method [21,33,34,35] has been designed to apply the finite element method (FEM) to bridge the full-field data with parametric constitutive relations. Recently, the method is generalized to an end-to-end training procedure to learn neural-network (or its counterparts) based constitutive relations from the full-field data [22].
The indirect data are given by N full-field deformation-load sequential data at n time snapshots Here the superscripts indicate time and the subscripts indicate the sequential number. We can compute the acceleration and the stress using the formulas Here ∆t is the time step, which is assumed to be constant. We train the neural network by solving a minimization problem arg min Here i is the index for time and j is the index for spacial degrees of freedom. It is worth mentioning, in Eq. (30), the predicted stress σ i j (θ) = M θ ( (u i j ), (u i−1 j ), σ i−1 j (θ)) depends on the predicted stress σ i−1 j (θ) from the last time step. The procedure of evaluating the residual at each time step in Eq. (30) is depicted in Fig. 2. The procedure resembles the recurrent neural network in deep learning, where the state variables are the stresses. Like the recurrent neural network, the stress predictions at different time steps must be sequentially computed and so does the back-propagated gradients, and thus the computation is hard to be parallelized. Additionally, the training of recurrent neural network suffers from exploding and vanishing gradients problem [36], which poses a challenge for indirect data training of SPD-NNs as well.

Initialization for Indirect Data Training
Another notable challenge in training the recurrent neural network based constitutive relation is the existence of many local minima, since the training data set is small in the present work, and different constitutive relations under the same external loads may produce similar displacements. A common strategy to find a good local minimum is to start from multiple random weights and biases. However, this strategy is expensive and does not guarantee a good initial guess. We thereby propose an initialization technique that yields a set of reasonably good initial weights and biases.
The key idea is that according to Eq. (4), we have the relation between the internal force and the stress field However, solving the stress field σ i j from Eq. (31) is generally an underdetermined problem, since the number of equations is fewer than the number of stress unknowns. To alleviate this, quadratic elements are used to represent the displacement field u, and in each quadratic element, a linear stress field is assumed and approximated (see Fig. 3). The number of stress unknowns is roughly 3(n x + 1)(n y + 1) for a 2D plate, and 6(n x + 1)(n y + 1)(n z + 1) for a 3D cube, the number of equations is roughly 2(2n x + 1)(2n y + 1) for a 2D plate, and 3(2n x + 1)(2n y + 1)(2n y + 1) for a 3D cube. Here n x , n y , and n z represent the number of quadratic elements in each direction. Therefore, the number of equations outnumbers the number of unknowns, and Eq. (31) becomes an overconstrained problem, the least-square fitting is applied for solving it. And then the stress field is approximated linearly at each Gaussian point in each quadratic element.
Once we solve for σ i j from Equation (31), we can use the technique in Section 3.2.1 to pre-train the neural network. Although the least square approximation of the stress field is poor (the error can be larger than 100%), but the approximation is qualitatively correct and sufficient for obtaining a good initial guess.

Applications
In this section, we present numerical results from solid mechanics for the proposed NN based constitutive relations: • Problem with a 1D truss coupon made of elasto-plastic materials under dynamic loading, which compares the proposed SPD-NN (16) and other neural network architectures, including σ-NN (13) and ∆σ-NN (14).
• Problems with 2D thin plates made of hyperelastic materials, elasto-plastic materials, and multiscale fiber-reinforced materials under dynamic loading, which demonstrate the effective-ness of the proposed SPD-NN (16) for learning path-independent, path dependent (hysteresis), and multiscale constitutive relations.
In all problems, the training data and test data of the strain/stress fields and the displacement fields are generated numerically.

1D Trusses with Elasto-plasticity
In this section, we consider 1D truss elements with a length L x = 1m, a cross section area A = 0.005m 2 and a density 8000kg/m 3 . The truss coupon, clamped on the left end, is tested under 5 loading conditions (see Fig. 4-left). The setting is corresponds to a uniaxial tensile test. The prescribed time-dependent load forcet consists of both loading and unloading parts and takes the formt The trusses are made of elasto-plastic materials. The young's modulus is The yield function with isotropic hardening has the form The yield strength is σ Y = 0.3 GPa, the plastic modulus is K = 200 9 GPa, the internal hardening variable follows the simplest evolutionary equatioṅ α =λ The truss coupon consists of 4 truss elements, which are modeled as geometric nonlinear truss elements [37, p. 63]. The time step size is ∆t = 0.001 s.
As for the SPD-NN Eq. (22), the estimated Young's modulus is assumed known, which can be separately calibrated using small external loads (e.g., using the method proposed in [22]), and the estimated yield strength isσ Y = 0.1 GPa, the transition function is here d = 0.1 is the nondimensional parameter.
Because the strain and stress pairs can be extracted from the uniaxial tensile tests, we can apply the direct input-output data training method Section 3.2.1 to train the neural network. Fig. 4-right shows the strain-stress curves exhibiting plastic loading and elastic unloading phenomena, which the neural networks are trained to learn. These curves start from the origin and rise with the slope E in the linear elastic region, until the stress reaches the yield stress. Then the curves enter the strain hardening region, where plastic deformations happen, until the elastic unloading. The slope of the unloading curve is typically equal to the slope in the elastic (initial) region of the stress-strain curve. The plasticity deformations result in permanent strains, which cannot be recovered by the elastic unloading. Since training neural networks involve highly non-convex optimization problems, we start from 10 different initial weights for all NN training.
To evaluate the quality of approximation, we propose two kinds of tests 1. NN test: extract the sequential strain-stress data ( i , σ i ) at each Gaussian quadrature point from the test data, and compare the predicted stress and the reference stress for each tuple ( i+1 , i , σ i ; σ i+1 ). 2. NN-FEM test: embed the learned constitutive relation M θ into the finite element framework, solve the governing equation of the truss coupon under the corresponding loading condition, and finally compare the predicted strain-stress paths and the truss deformations.
Because the optimization results depend on the initial guess for the neural network weights and biases, we pick the neural network with the minimal loss on the training set for the NN test and NN-FEM test. Other choices for selecting the candidate neural networks such as cross-validation [38] can also be used. It is worth noting that NN-FEM test is more challenging than the NN test, but is more relevant for predictive modeling. Indeed, the NN test does not take the numerical stability in the predictive modeling into account and only evaluates the NN's ability to fit the strain-stress curve. As we will see in the following examples, although SPD-NNs and other NN methods fit the strain-stress curves equally well in the NN test, SPD-NNs are significantly preferable due to their predicting power in new scenarios. For the NN test, the predicted strain-stress relations are reported in Fig. 6. All kinds of NN architectures give good results. For the NN-FEM test, the predicted displacement trajectories of the right end point and the predicted strain-stress curves for one Gaussian quadrature point on the right end element are reported in Fig. 7 and Fig. 8. Most of σ-NN and ∆σ-NN architectures are found to be numerically unstable. We attribute the instability to the deviation of the strain and stress pairs in the test loading condition from the training set, where the numerical errors are accumulated during the prediction process. More specifically, the violation of SPD constraints on the tangent stiffness matrix for σ-NN or ∆σ-NN makes them vulnerable to these errors and less robust to even slight extrapolations. In contrast, the proposed SPD-NN delivers numerically stable results and the predicted displacements and strain-stress curves overlap the exact ones.

Comparison of Different Neural Network Architectures for SPD-NN
In this experiment, we consider different neural network architectures. We vary widths, depths and activation functions of the neural network used in SPD-NN, while keeping other settings the same as Section 4.1.1. The errors in the NN-FEM test in terms of mean squared errors of the predicted displacements at the final time T are shown in Table 1. We see that the tanh activation function is in general more accurate than others if appropriate widths and depths of the neural network are chosen (For this case, the training set consists of 800 data, the number of neural network parameters is not supposed to outnumbers it too much). However, ReLU and leaky ReLU are more robust for deep and wide neural networks, despite being less accurate.

Choice ofσ Y and d
We consider the impact ofσ Y and d in the transition function D (see Eq. (32)) on the accuracy of the SPD-NN. As mentioned before,σ Y controls where the transition from the elastic form (Eq. (20)) to the plastic form (Eq. (21)) happens and d controls the sharpness of the transition (see Fig. 1). In this numerical experiment, we use a fully connected neural network with 3 hidden layers, 20 neurons in each layer, and the tanh activation function. We varyσ Y and d in Eq. (32). The error plot is shown in Fig. 9, where the error metric is the same as Section 4.1.2. The reported errors are the average errors of 10 simulations with different initial guesses. We can see that the accuracy of the SPD-NN is less sensitive to d. Additionally, as long as we choose a small enoughσ Y so that the corresponding plastic form (Eq. (21)) covers the plasticity regime, the SPD-NN is sufficiently expressive to approximate the constitutive relation. This justifies our choices d = 0.1 andσ Y = 0.1 GPa.

2D Thin Plate with Different Materials
In this section, we consider thin plate coupons of size L x = 10 cm by L y = 5 cm with the plane stress assumption (L z = 0.1 cm). These plates are made of different materials, including hyperelastic material (finite deformation), elasto-plastic material (infinitesimal deformation), and fiber-reinforced multiscale elasto-plastic material (infinitesimal deformation). For each case, the plate is tested under 13 loading conditions as depicted in Fig. 10 6.6×10 −5 NaN NaN NaN Table 1: Mean squared errors of the predicted displacements at the final time T for different neural network architectures, including tanh, ReLU (rectified linear unit), ELU (exponential linear unit), SELU (scaled exponential linear unit), and leaky ReLU (leaky rectified linear unit) with leakage 0.1. The width is the number of activation nodes in each hidden layer. The depth is the number of hidden layers. "NaN" denotes the numerical simulation fails (e.g., due to numerical instability) using the trained NN-based constitutive relations. • C) clamp on the left edge and impose force load on the bottom edge.
The total simulation time is T = 0.2s. We use A1-A5 and B1-B5 as training data and A6, B6, and C1 as test data. Both training procedures discussed in Sections 3.2.1 and 3.2.2 are applied. For the direct input-output data training, the strain-stress sequential data are extracted from all Gaussian points in the training sets. For the indirect data training, the full-filed displacement fields on the 21 by 11 grid (0.5 cm interval) from the training data are extracted. Therefore, this approach is potentially applicable to experimental data. The pre-training is required to obtain good initial guesses.t

Hyperelasticity
The plate is made of the incompressible Rivlin-Saunders material [39,40] with the density ρ = 800 kg/m 3 and the energy density function Here c 1 = 0.1863 MPa and c 2 = 0.00979 MPa, and I 1 , I 2 , I 3 are three scalar invariants of the right Cauchy-Green stretch tensor C = F F T = 2 + 1, where The incompressibility implies that J = 1. The plate is assumed to undergo finite deformations (see Remark 1), and the second Piola-Kirchhoff stress tensor reads here λ J is the Lagrangian multiplier, which can be calculated based on the plane stress assumption S 33 = 0. The plate domain is discretized by 20 × 10 quadratic quadrilateral elements. The time step size is ∆t = 0.001 s. The data sets are generated with load parameters (p 1 , p 2 , p 3 ) = (44800, 4480, 16800) N/m. Both the direct data training approach Section 3.2.1 and the indirect data training approach Section 3.2.2 are applied to train a SPD-NN: where the neural network consists of 4 hidden layers and 20 neurons in each layer. The predicted trajectories of displacements at the top-right and top-middle points as a function of time and the references for all test cases are depicted in Fig. 11. All predicted results are in good agreement with the references, and the direct input-output data training approach leads to slightly better results for case C1. The predicted von Mises stress fields at t = T 2 for all test cases and the references are depicted in Fig. 12.

Elasto-Plasticity
The plate is made of titanium, which is assumed to be elasto-plastic material with density ρ = 4200 kg/m 3 . The constitutive relation is σ = C here C denotes the isotropic plane stress stiffness tensor with Young's modulus E = 100 GPa and Poisson's ratio ν = 0.35. The von Mises yield function with isotropic hardening has the form The yield strength σ Y = 0.97 GPa and the plastic modulus K = 10 GPa, the internal hardening variable α follows the simplest evolutionary equatioṅ α =λ (36) This plate domain is discretized by 20 × 10 quadratic quadrilateral elements. And the time step size is ∆t = 0.001s. As for the SPD-NN Eq. (22), the estimated yield strength isσ Y = 0.32GPa, the transition function is here σ n vm is the computed von Mises stress at the previous time step and d = 0.1 denotes the nondimensional parameter. The tangent stiffness matrix C θ in the linear region is first estimated as following, and then used as constant in Eq. (22). The indirect data training approach in Section 3.2.2 is applied to extract the tangent stiffness matrix, and we obtain the following estimation For each components, the relative error is less than one percent. These errors are introduced by the discretization, including the interpolation of the displacement field on the observation grid and the estimation of the acceleration Eq. (27). It is worth mentioning, the direct input-output data training delivers exactly the same stiffness matrix as the reference.
Nonlinear region. The data sets are generated with load parameters (p 1 , p 2 , p 3 ) = (1.6, 0.16, 0.6) GN/m which are 10 times larger than these in the linear region. Both direct input-output data training (Section 3.2.1) and indirect data training (Section 3.2.2) are applied to train a SPD-NN with 5 hidden layers and 20 neurons in each layer. The predicted trajectories of the displacements at top-right and top-middle points for all test cases are depicted in Fig. 13, along with the references. SPD-NNs trained with both methods are able to predict the initial elastic behavior, the strainhardening region, and the unloading behavior. The SPD-NN obtained by the direct input-output data training performs better especially for the prediction of the yield strength the strain-hardening behavior. The predicted von Mises stress fields at t = T 2 and the references for all test cases are depicted in Section 4.2.2. Reasonable agreements are achieved.

Multiscale Fiber Reinforced Elasto-plasticity
The plate is made of the titanium-the same as Section 4.2.2-but reinforced by fibers made of SiC, which are assumed to be isotropic and elastic with ρ = 3200 kg/m 3 , E = 400 GPa, and ν = 0.35 These fibers are square shaped and uniformly distributed in the plate, with a diameter d = 0.25 cm and a fraction 25%. There are in total 800 fibers, as shown in Fig. 15. This plate domain is discretized by 200 × 400 quadratic quadrilateral elements, and 25 elements for each fiber. And the time step size is ∆t = 0.001s. which are 10 times larger than those in the linear region. Only the indirect data training approach is applied to train a SPD-NN with 5 hidden layers and 20 neurons in each layer. To enable direct input-output data training, homogenization is required to generate strain-stress data from RVE simulations [13,14,17,6], or extract strain-stress data from direct numerical simulations. This may be challenging in the context of experimental data, which necessitates the indirect training approach.
The predicted trajectories of displacements at top-right and top-middle points for all test cases are depicted in Fig. 18, along with the references. Comparing with the previous elaso-plasticity case, all displacements are smaller due to the SiC fiber reinforcement. The proposed SPD-NN gives a satisfactory approximation of the initial elastic behavior, the strain-hardening region, and the following unloading behavior. The predicted von Mises stress fields at t = T 2 and the references for all test cases are depicted in Fig. 19. The predicted and simulated homogenized stress fields are in reasonably good agreement. Although the solutions obtained by SPD-NNs do not capture local large stress concentrations near each fiber (at the level of the microstructure), the local recovery techniques [41] can be applied to estimate these local stress concentrations. For example, a RVE simulation can be conducted to estimate local stresses using the local strain from the coarse homogenized solution.
Accelerating Simulations with Neural Network Surrogates for Constitutive Modeling. It is also worth noting that each SPD-NN-based simulation is several order magnitude faster than the corresponding fiber resolved simulations. The CPU time for both the SPD-NN-based simuations and fiber resolved simulations are shown in Fig. 20. Note as we increase the external load, the CPU time increases because more elements undergo plastic deformations, which requires more expensive Raphson-Newton iterations in the numerical simulations. Still, the dramatic acceleration from around 24 hours to just a few minutes is impressive. However, the acceleration should be carefully interpreted in the context of surrogate models. First, the current benchmarks are based on serial execution, where the state-of-the-art FEM simulations usually involve parallelization. Nevertheless, these parallelization techniques for fiber resolved simulations are also directly applicable to SPD-NNs. Second, the speed of the fiber resolved simulations depend on the required resolution (e.g., the resolution on each fiber). If we use very coarse grids, the speed of the fiber resolved simulations may be comparable to or even faster than the SPD-NN-based simulations, though the coarse grids raise concerns on accuracy. Third, we have not compared the present model with other state-ofthe-art accelerating techniques for multiscale modeling, for example projection-based model order reduction [42] and self-consistent clustering [43]. Finally, as with any NN-based surrogate models, the SPD-NN only works on test data that does not deviate too much from the training data, where fiber resolved simulations are usually considered to be general and applicable in a much wider context.

Discussion on Neural Networks
The introduction and benchmarking of SPD-NNs raise many questions, some of which are partially answered in this paper while some others are worthwhile to be investigated further. We believe that the following discussion is important and central to the application of NN-based constitutive modeling and should be careful dealt for the development of SPD-NNs and other NN-based approaches.
Function approximators. In the present work, neural network is used as a basis function to approximate a complex constitutive relation, the strain-stress relation. Many other basis functions, including piecewise linear functions, radial basis functions, and radial basis function networks, can also be applied for the approximation of the Cholesky factor. The choice of neural networks is well-versed and justified by the following three reasons.
• First, the distribution of the strain data-the input to the neural network-is determined by experiment records, which in general does not evenly spread over the domain. The comparison study presented in [22] illustrates that neural network outperforms the other basis functions in terms of regularization and generalization properties when the data distribution is ill-behaved.
• Second, the strain-stress curves are non-smooth in the context of plastic deformations and elastic unloading. As a result of the non-smoothness, the approximation efficiency of typical basis functions is usually compromised. Nevertheless, neural networks exhibit the potential to capture the sharp transitions in the strain-stress relations and performs reasonably well for the non-smooth data.
• Third, the input and the output dimensions are relatively high (e.g., the input is 9D and the output is 4D in the elasto-plasticity and the multi-scale cases), which poses a chanllenge for traditional basis functions. For example, if we were to use the linear basis functions, to discretize the high dimensional input space, even if we only have 10 grid points per dimension, the total degrees of freedom is 10 9 , which is too costly regarding both the computation and the storage. Yet, neural networks are particularly convenient and useful for expressing mappings between high dimensional spaces [44] and requires no mesh in the input parameter domain.
All the aforementioned reasons motivate us to use neural networks in this work.
Optimization method. Most of neural networks in literature, especially in computer science communities, are trained with stochastic gradient methods (SGD). In our present work and the previous work [22], the optimization of the loss functions Eqs. (26) and (30) is done by the Limited-memory BFGS (L-BFGS-B) method [45] with the line search routine in [46], which attempts to enforce the Wolfe conditions [45] by a sequence of polynomial interpolations. Note BFGS is applicable in our case since the data sets are typically small and the neural network is reasonably deep and wide; otherwise, the memory requirement of L-BFGS-B is so high that SGD or other similar firstorder methods for training neural networks should be adopted. It is worth mentioning that the choice of BFGS optimizer is well-motivated and has long been adopted for scientific and engineering applications due to its fast convergence and robustness.
Data scaling. We observed that input data scaling significantly helps train SPD-NNs faster, reduces overfitting, and makes better predictions. In the present work, the inputs and outputs of the neural network are scaled to a similar magnitude. Specifically, we introduce a strain reference ref and a stress reference σ ref to scale strains and stresses. We also scale the tangent stiffness matrix by σ ref ref . For example, in the elasto-plasticity case, the Cholesky factor has the form , here E is the estimated Young's modulus. This guarantees the inputs for SPD-NN, especially for elasto-plasticitiy, ( i , i−1 , σ i−1 ) have similar magnitudes.
Local minima. Local minima are observed in Section 4.1, since training neural networks involves highly non-convex optimization problems. Although neural networks with minimal losses on the training set from 10 different initial weights perform well, neural networks with median losses are less satisfactory. This reveals the uncertainty with respect to the initial weights for most NN-based data-driven approaches. Wider neural networks will be considered in the future, since some theoretical and computational results [47,48] show quality of local minima tends to improve toward the global minimum value as depths and widths increase. As for the indirect data training approach, pre-training approach (see Section 3.2.2) produces acceptable initial weights and thus relieves the concern. For all these training processes, the optimization is terminated when the objective function is called 3000 times and 50000 times for pre-training and training, respectively. The direct input-output data training is about 4 times faster than the indirect data training.

Conclusion
Data-driven approaches continue to gain popularity for constructing constitutive models from high-fidelity simulations and high-resolution experiments. The incorporation of data-driven constitutive models into conservation equations leads to a hybrid model, namely a coupled system with differential equations to describe conservation laws and neural networks to describe the material properties.
To make these hybrid models numerically more robust, we introduce a novel neural network architecture, SPD-NN, where the neural network outputs the Cholesky factor of the tangent stiffness matrix instead of the stress or the stress increment. This neural network architecture weakly imposes convexity on the strain energy function (i.e., SPD tangent stiffness matrix). The incremental form of SPD-NN also preserves the time consistency. SPD-NN-based constitutive relations are tested on a 1D elasto-plastic truss problem and several 2D plate problems in which the plate is made of hyperelastic, elasto-plastic, and multiscale fiber-reinforced materials. When contrasting the SPD-NN with two other neural network architectures, we show that the SPD-NN exhibits better numerical stability in the resulting hybrid models. The general training approach and the improved numerical stability allow for the potential to extend SPD-NNs to other time-dependent physical systems, such as viscoelastic materials, where the constitutive relations are rate-dependent, and plastic materials with stronger hysteresis, where more history-dependent variables are required.
However, one limitation of the current approach is that the training process requires full field data, either for strain-stress pairs or displacement measurements. This may not be always possible; for example, the measurements may only be made on the surface of a 3D solid body. NN-based constitutive modeling with incomplete data remains to be investigated in the future.