When deriving models to predict path-dependent responses of materials that lack material symmetry or exhibit complex size-dependent behaviors, the smallest number of variables (e.g., stress measures, stress invariants, internal variables) required to replicate constitutive responses increases, and the dimension of the parametric space in which the model is formulated also becomes higher. This increase in variables often leads to the increase of material parameters due to the need for additional support to control the geometry in the high-dimensional space.
As such, modelers must decide a trade-off between simplicity and sophistication [1]. Increasing the number of material parameters is often undesirable or only used as the last measure to capture the phenomenology precisely [2]. This preference for simpler models is not limited to constitutive theories of solids and has a long history in science (p.398, Newton [3])and philosophy [4±6]. In the early ages of the development of plasticity theories, experimental data were relatively limited in quantities and lacked precision afforded by the state-of-the-art instruments [1,7,8]. Hence, assumptions on material symmetry and limitations on the number of variables used to describe the deformation mechanisms become common strategies that maximize what [1] described as the trade-off between simplicity and sophistication.
Increasing the precision of a material model may unavoidably increase the least number of material parameters necessary to describe the material behaviors. As such, many simpler models that involve more hypotheses often become the blueprint of more elaborated models of higher dimensions. This approach is often adopted by plasticity models in, for instance, the following scenarios.
1. There is a need to capture additional causal mechanisms that are not describable with the current existing variables [9]. For instance, the anisotropy and non-coaxiality induced by fabric evolution require incorporating the fabric tensors into the constitutive laws [10±12]. 2. There is a need to enrich the description of a model such that the model can be further generalized for broader applications. For instance, the Drucker±Prager model can be viewed as a generalization of the Von Mises model by introducing the dependence of the yielding behaviors with respect to the pressure (position of the hydrostatic axis in the principal stress space). de Borst [13] introduces micropolar constitutive models for geo-materials by introducing the additional couple stress terms and length scale parameter in a Drucker± Pragger model. Multi-physics constitutive models, such as [14], also adopt an extension strategy where a pure solid mechanics model is enhanced by introducing additional variables such as temperature, degree of saturation, and volume fraction of void that are necessary to capture the multiphysical coupling [14±18]). 3. Finally, there have also been cases in which a model is amended to circumvent known limitations. A classic example is the usage of micropolar theory to circumvent the pathological mesh dependence of plasticity models in the softening regimes. (e.g. [19±21]).
Nevertheless, this incremental strategy to amend models exhibits disadvantages. First of all, obtaining plasticity data from either representative elementary volume simulations or experiments could be costly. As such, the lack of foresight of the model dimensions may lead to a biased data acquisition strategy where data might not distribute well enough to calibrate and test the learned model. For models expressed in higher-dimensional space (e.g., anisotropic elasticity model that requires not just the strain invariants but also the principal directions, micropolar and micromorphic models that require higher-order terms to fully describe the kinematics), data may appear to be sparser and hence requires more data points to characterize the behaviors of comparable complexity [22].
Machine learning approaches with advanced architecture, such as the recurrent neural networks for sequential learning [23±26], 1D convolutional neural network [27], and transformer with attention mechanisms [28], can be, in theory, trained to replicate these high-dimensional constitutive responses. However, the vanishing/exploding gradients, the increased computational cost, and the increased demand for data can all become the bottleneck of these approaches. Furthermore, the higher dimensionality of the model also makes the learned model more vulnerable to overfitting (incapable of generalized prediction) [29] and the results more difficult to interpret properly [27,30,31].
On the other hand, there are multiple attempts to represent the plastic yield surfaces via parametrized surface or implicit functions. Vlassis and Sun [27], for instance, introduce the usage of the signed distance function to implicitly represent the yield surface in the stress-internal-variable space. The additional signed distance property enables the plastic flow to be a unit gradient which simplifies the calculation of the plastic multiplier, and allows one to incorporate the plastic flow direction into the Sobolev training of yield function. The resultant yield function is then parametrized via an MLP architecture. Meanwhile, Coombs et al. [32],Coombs and Motlagh [33,34] leverage the flexibility afforded by non-uniformed rational B-splines (NURBs) to parametrize yield function. Xiao and Sun [35], on the other hand, represents the yield surface as a manifold and train neural networks to paramterize an atlas of coordinate charts to form three-variable yield surface of complex shapes.
While the component/modular-based learning approach [36,37] enables one to reformulate the plasticity learning problem without enforcing the memory effect through the architecture, the challenges of training and interpreting those models remain. Another less explored route to model plasticity is to use the kernel method in which a nonlinear transformation maps the material data in a higher-dimensional feature space to make it easier to perform clustering (e.g., unsupervised learning) or insert hyperplane (e.g., classification). The advantage of the kernel method is that it is generally more robust (due to the non-parametric nature and the implicit feature mapping without any explicit feature engineering) and less sensitive to outliers [38]. More importantly, the construction of a higher-dimensional feature space also makes it possible to interpret the relations between the features and the learned plasticity models. However, a key technical barrier of the classical kernel method is that the performance of the method is highly sensitive to the specific kernel function (e.g., polynomials, radial basis function) used to generate the feature space. As different data sets may require different kernel functions to achieve good performance, selecting the kernel function becomes a time-consuming trial-and-error process.
This paper introduces a neural kernel approach to generate micropolar and micromorphic plasticity models. Here we want to leverage the robustness and interpretability of the kernel method afforded by the feature space while generating a data-dependent kernel tailored to our specific need to capture the high-dimensional constitutive responses of the materials with complex internal structures. This treatment provides us with a unified data-driven approach to recognize the pattern of the data (through the data-dependent kernel), regardless of the data dimensions. Combined with a simple kernel ridge regression, we may generate a yield function of arbitrary input space dimensions without explicitly handcrafting the feature space. This trait is input for us to automate the process of generating yield surface from data of dimensions higher than 3, especially when data is sparse.
The organization of the rest of the paper is as follows. For completeness, we first review the plasticity theory of micropolar and micromorphic materials and the related Hill±Mandel lemma necessary for generating the multiscale data for computational homogenization. This review explains the difficulties in formulating yield surfaces in a high-dimensional stress space (Section 2). We then introduce the neural kernel method formulated for the highdimensional space. In particular, we explain both the theory of the neural kernel method, as well as the strategy we adopted to train both the data-dependent neural network kernels and obtain the coefficients used to interpolate the yield function in the high-dimensional feature space. The return mapping algorithm that adopts the neural kernel yield function is also included (Section 3). This is followed by a collection of representative numerical examples in which we provide verifications and demonstrate how the proposed approach can be applied to complex data obtained from direct numerical simulations of layered pressure-sensitive materials (Section 4). We then summarize our major findings in Section 5. Additional numerical tests and the detailed procedure necessary for third-party inspection and validations are provided in the Appendix.
For completeness, we review the theory of micropolar and micromorphic continua (i.e., higher-order continua). Cosserat and Cosserat [39] is credited with introducing the first high-order continuum theory in which the concept of micro-rotation is introduced to describe the effects of internal structures on the constitutive responses. Different variations, for instance, the couple stress theory [40±43], which derives energy density as a function of both strain and the curl of strain. The generalized micromorphic continuum theory, which introduces the concepts of micro-deformation and the corresponding energy-conjugate stress measures as a generalization of the kinematics, have been introduced in the 1960s [44±47]. The extension of the higher-order continua theory to the finite deformation range has been formulated by Toupin [41] where an action density is derived such that it is invariant under the group of Euclidean displacements in a Hamiltonian mechanics framework. MÈ uhlhaus and Vardoulakis [48] and Peerlings et al. [49] further extended the higher-order theories for elastoplasticity problems and examined the regularization effect of higher-order continuum theories. Steinmann [50] extends the multiplicative kinematics theory to formulate micropolar elastoplasticity in the geometric nonlinear regime. Recently, the multiscale micropolar [51] and micromorphic [52] constitutive modeling has been derived in the finite deformation regime, while [53] have introduced a linear relaxed version of micromorphic models to capture the wave propagation in meta-materials. A comprehensive review of connections between the higher-order and non-local continuum theories can be found in [54].
For simplicity, we restrict our learned models to be within the infinitesimal deformation regime. We then introduce machine learning to generate three classes of plasticity models based on Cauchy, micropolar, and micromorphic continuum theories. The micro-deformation χ i j describes the configuration of the directors, i.e., the internal micro-structures (e.g., voids and inclusions) contained in the material point sampled at an arbitrary position x, as shown in Fig. 1. The material point of an effective medium is also called the representative volume element (RVE). For a Cauchy continuum, the directors are much smaller than the RVE, so the kinematic configuration of the directors χ i j mechanically affects the RVE much less than the strain tensor,
such that the χ i j can be neglected. In comparison, the micromorphic theory considers that the size effect of the deformable directors is not negligible. As such χ i j must be considered to describe the distortion of the RVE internal structure. Micropolar continua is a sub-class of micromorphic continua in which the directors can be assumed to be rigid (e.g., rigid inclusions) such that
where θ k describes the rigid rotation and ϵ i jk is the Levi-Civita permutation symbol. Based on the kinematic configuration of higher-order continua described by u i and χ i j , the boundary value problems can be solved given the constitutive laws of micropolar and micromorphic continua, which are reviewed in Section 2.1, followed by the RVE homogenization scheme that models the constitutive law by multiscale simulation (Section 2.2) and data-driven approaches to learn the constitutive law (Section 2.3).
To solve the boundary value problems (BVP) of higher-order continua for the kinematic configuration u i and χ i j , we need to consider the kinematic relation, balance of linear and angular momentum, and constitutive law as summarized in the upper half of Table 1 [55]. The kinematic relation defines the strain tensor as the difference between the displacement gradient u i, j and the micro-deformation χ i j or micro-rotation -ϵ i jk θ k ; the gradient of micro-deformation G i jk or the gradient micro-rotation κ i j (i.e., curvature tensor) are also introduced to describe the higher-order deformation. The balance law of linear and angular momentum states the relationship between the Table 1 The summary of BVP components of the micropolar and micromorphic materials.
Kinematic relation
Cauchy stress σ ji and higher-order couple stress m ji or generalized stress ζ i jk , which are computed given ε i j , κ i j , and G i jk based on the constitutive law.
To complete the solution of BVP, the elastoplastic constitutive law in the lower half of Table 1 needs to be found to model the relation between the kinematic modes and the higher-order stresses, i.e., between ε i j and σ ji , between κ i j and m ji , and between G i jk and ζ i jk . The constitutive law consists of the elasticity model, yield function, KKT condition, and plastic flow rule, as summarized in Table 1. The elasticity model finds the elastic energy functional W such that the elastic stresses are work conjugates of the kinematic modes. The core component of the plasticity model [56] is the yield function f that defines the onset of plastic yielding at f = 0. The elastic region in the stress space requires f < 0. When the stresses touch the yield surface, then f = 0. The permanent plastic deformation grows at the rate of Λ and in the direction of gradients of f according to the associative flow rule, where Λ is called the plastic multiplier. The KKT condition shown in Table 1 implies that if the yield surface is not touched, equivalent to f < 0, then Λ = 0 and no permanent plastic deformation is growing.
The micromorphic elastoplastic constitutive law can then be derived as the return mapping algorithm, i.e., given current elastic deformation (ε e i j (t), G e i jk (t)) and incremental deformation (∆ε i j , ∆G i jk ), finding the stresses at the next time step (σ ji (t + ∆t), ζ i jk (t + ∆t)).
The first step is to find the trial stresses (σ tr ji , ζ tr i jk ), assuming the incremental deformation is elastic, which is equivalent to f (σ tr ji , ζ tr i jk ) < 0 as shown in Eq. ( 3).
[ σ ji (t + ∆t)
If the incremental deformation is not elastic, i.e., f (σ tr ji , ζ tr i jk ) ≥ 0, then a correction of the trial stresses should be made to ensure f (σ ji (t + ∆t), ζ i jk (t + ∆t)) = 0, and the final stresses can be found based on the associative plastic flow rule as shown in Eq. ( 5). The elasticity tensors are defined as,
The incremental stress update expressed in Voigt notation reads,
The two equations above summarize the return mapping algorithm for micromorphic continua. The micropolar continua is a special case of the micromorphic continua where the micro-deformation is restricted to be rotational only. As such, the return mapping algorithm can be implemented by re-expressing the micro-rotation gradient κ i j and the couple stress m ji in terms of the micro-deformation gradient G i jk and generalized stress ζ i jk , i.e., The micromorphic return mapping algorithm then return the constitutive updates and the gradient of micro-rotation and couple stress can then be recovered via
ϵ imn G mn j , -ϵ imn ζ mn j ). ( 7)
The constitutive relation of the material point at position x, shown in Fig. 1, is obtained by studying the RVE domain Ω with the local coordinate Y as shown in Fig. 1. The kinematic modes denoted as εi j , κi j , and Ḡi jk are prescribed by deforming the RVE with the displacement boundary condition shown in Eq. ( 8) and Fig. 2 [51,52], where the boundary condition is the linear combination of the prescribed kinematic modes. This boundary condition is admissible, as proven by Hill's Lemma in Appendix A.
The local boundary value problem is then solved given the prescribed boundary condition to solve for the local displacement and stress field. The homogenized stress and generalized stress [52] can be computed as Eq. ( 9), which satisfies the Hill±Mandel's condition as proven in Appendix A.
for micropolar continua
The elastoplastic constitutive law can be sufficiently modeled by machine learning by training the elastic energy functional and the plastic yield function independently. The elastic energy functional, W (ε e i j , κ e i j ) or W (ε e i j , G e i jk ), can be simply modeled by a neural network consisting of multi-layer perceptrons (MLP), with the strain measures as the input and the elastic stored energy as output as the output. In our implementation, we adopt the Voigt vectorized notation used in [57] to train the neural network elasticity models. For brevity, the supervised learning of elasticity for micromorphic continua will not be discussed in great detail. Interested readers may refer to, for instance, Vlassis and Sun [27]. On the other hand, the supervised learning for the yielding function and the corresponding hardening laws are formulated in the next section.
Remark 1 (Neural Network Architecture). The architecture of the MLP, which we adopted in this study, is shown in Eq. ( 10). The 3-layer architecture is composed of neurons equipped with Exponential Linear Unit (ELU) activation function. The activation function of the output layer A 3 can be ELU or identity map, depending on how the MLP is used: A 3 should be an identity map mostly except when the MLP is used to construct the neural kernel (NK) architecture.
, where
This neural network design is used for both the elastic stored energy functional as well as the yield surface because the derivatives of ELU are sufficiently differentiable. This smoothness may improve the robustness of the optimization process and alleviates the vanishing and exploding gradient problems.
This section describes (1) how to use the neural kernel (NK) method to reconstruct the yield surface for micromorphic continua and (2) provides the implementation details necessary to incorporate the learned model into a return mapping algorithm. Here, we represent the yield surface in a multi-dimensional parametric space via an implicit scalar signed distance level set function f (x), such that the yield surface geometry is recovered at f (x) = 0 where x stores the higher-order stress components of the implicit yield function and the internal variables.
The general framework of our NK method follows the work of Williams et al. [38] on Neural Kernel Field (NKF), while we adopt the kernel function architecture as deep neural networks for the generalizability to arbitrary dimensions. A general workflow of this framework is presented in Fig. 3 containing three major steps:
1. Generate the labeled narrow band level set data given the yield stress point x and plastic flow direction n as discussed in Section 3.1, 2. Train the kernel coefficients for the kernel function associated with the NN-based feature map φ θ , after which a two-step training may be needed for the NN weight and bias θ to ensure that the sign of the yield function is correct, as shown in Algorithm 1, and 3. Predict the level set yield function f θ (x) as a linear combination of the basis kernel functions and locate the surface at f θ (x) = 0 (see Section 3.2).
The workflow shown in Fig. 3 presents the 2D and 3D views of the surface for the readability purpose, but the NK model is able to reconstruct a much higher-dimensional yield surface. After training the NK-based yield function, the elastoplastic constitutive law is reproduced by the return mapping algorithm presented in Section 3.3 and Algorithm 2.
The raw data set for surface reconstruction consists of point coordinates sampled from the ground-truth surface and the corresponding normal vectors, such that the data set is described in the form of
where x i are point coordinates and n i are the surface outward unit normal vector at x i ; in the context of the plastic yield function, x i are yield stress points and n i are the plastic flow direction. For supervision purposes, we create two labeled datasets following the concept of the narrow band level set [58]. The first data set D is generated by perturbating the spatial coordinate of surface points in the normal direction, and labeling them by the distance to the true surface as follows:
where ϵ is a small number indicating the distance of perturbation for the surface points. The first data set controls the trained level set function to be zeros at the surface, and the gradient of the function is equal to the unit normal as shown in Fig. 4. In our numerical examples, the data are centered and scaled such that the average norm of each data point becomes one, and ϵ is tuned as a hyperparameter between 0.01 and 0.1.
The second set of labeled data V is needed to ensure the occupancy condition [59], i.e., the function should be negative inside the surface and positive outside the surface, and there should not be additional holes in the region enclosed by the surface, as shown in Fig. 4. V is sampled in the d-dimensional Euclidean space excluding the points on the surface, where y vol = 1 if x vol is a point outside the surface, and y vol = -1 otherwise.
This subsection describes the proposed NK method for surface reconstruction in detail. Before we go through the general framework of NK, we present a brief overview of classical kernel methods to facilitate further demonstrations of the NK method. Classical kernel methods for regression problems adopt a pre-defined kernel function with a set of trainable kernel coefficients α i to approximate the function y(x) given a labeled dataset D = {(x i , y i )} as follows:
where the hat over ŷ indicates an approximation of function y. The kernel coefficients α j are trained by directly solving the following system of linear equations:
where λ is the tunable hyperparameter used for regularization and data denoising.
It is theoretically proven that the conventional kernel method does not accurately make predictions of unseen features if the spatial dimension d is much smaller than the number of data |D| [60]. In this sense, the neural network is adopted in kernel methods in order to increase the representation power of this machine learning model, where a trainable NN φ θ is introduced to create a map from a d-dimensional input space to an h-dimensional feature space. The kernel regression is then implemented in the feature space enforcing h and |D| in the same order of magnitude. The φ θ implemented in [38] follows an architecture similar to C-OccNet [59] and is not applicable in higher-dimensional space. For generalizability to higher-dimensional inputs, we establish φ θ as a multi-layer perceptron (MLP).
Instead of the kernel map of the original feature vector x shown in Eq. ( 13), the arguments of the kernel function are replaced by the feature map φ θ (x). As a result, we introduce the MLP weight θ in addition to coefficients α j as the trainable parameters, and the architecture of the neural kernel f θ : R d -→ R is described as follows:
where K N S (x, y) = ψ T (x)ψ( y) is the neural spline kernel [61] induced from a single-layer non-trainable neural network ψ. In this paper, K N S is replaced by K ′ (x, y) = x T y such that the neural spline kernel is integrated into φ θ as an additional layer; which results in the following simplified NK architecture:
We next present the loss function L(α j , θ) adopted as follows:
where γ is a tunable hyperparameter. By minimizing this loss function, we achieve two goals in recovering the correct surface geometry: (1) constrain the value of f θ to zero and the gradient of f θ to fit the actual surface normal direction, which is satisfied as the first term goes to zero; and (2) enforce the occupancy condition so that f θ predicts the correct sign inside and outside the yield surface, which is satisfied as the second term goes to zero. Notice that the second goal characterizes a binary classification problem, but the conventionally used binary cross entropy (BCE) loss is not included in the loss function, which is because the usage of logarithm becomes problematic with negative values in this case.
In the training process, α j and θ are updated in an asynchronous fashion: α j are trained by directly solving Eq. ( 14), but θ is updated using gradient-based optimization given fixed α j and ∂ L/∂θ . We summarize the training routine in Algorithm 1.
Require: Training data set D, occupancy data set V , learning rate η, hyperparameters λ and ϵ.
1. Setup the feature map neural network φ θ (x) with the trainable parameters θ . Define
Given pairs of feature vectors x i , x j from D, compute the kernel matrix K i j
Assemble K i j = φ θ T (x i )φ θ (x j ) for (x i , y i ), (x j , y j ) ∈ D 3. Solve the linear kernel equations for the kernel coefficients α j .
Solve (K i j + λδ i j )α j = y i for (x i , y i ) ∈ D. 4. Given the level set yield function f θ (x), compute the loss function L(α j , θ ) and its gradient ∂ L/∂θ.
Compute
Differentiate L(α j , θ ) for ∂ L/∂θ with loss.backward(). 5. Update θ given η and ∂ L/∂θ using gradient-based optimizer like ADAM.
Run torch.optim.Adam(θ ,η).step(). 6. Repeat Steps 2-5 until the loss function converges and output f θ (x).
Remark 2 (Hyperparameter Tuning).
1. γ < ϵ < r min can be a set of good hyperparameters in Eqs. (11) and (17), where r min is found by first centering the data, i.e. translating the data such that the centroid goes to the origin, and then computing the distance from the closest data point to the origin. 2. The output dimension h should be large enough, which is equivalent to increasing the dimension of the basis level set function shown in Fig. 3. Ideally, h should be comparable with |D|. 3. Increasing λ would also improve the convergence of the loss function when the data are noisy. In both examples of this paper, λ = 0.01 is used.
In the Sobolev training technique with deep learning, we consider the influence of the derivative of the neural network function with respect to the network input in the loss function, such that we enforce the prediction accuracy for the neural network derivative with respect to its input. For the global loss function, we directly supply the original loss in Eq. ( 17) with the MSE between the groundtruth and predicted surface normal directions:
where L Sob is the Sobolev loss function we adopt, γ ′ is a hyperparameter controlling the influence of derivative term in the global loss. We will compare NK with the MLP-based level set (MLP-LS) method, which is documented in [27]; the details of the MLP-LS method are provided in Appendix C.
The return mapping algorithm for micromorphic materials is presented in Algorithm 2, which follows the return mapping theory shown in Eqs. ( 3) and (5). We further assume that the yield surface is able to evolve, and such hardening process is governed by an internal variable Λ, i.e., the magnitude of the cumulative plastic strain (plastic multiplier), such that f θ (x, Λ) represents a family of yield functions evolving with Λ. For micropolar materials, convert the (κ i j , m ji ) into (G i jk , ζ i jk ) = (-ϵ i jl κ lk , -1 2 ϵ i jl m lk ) before Algorithm 2 and convert back by (κ i j , m ji ) = (- 1 2 ϵ imn G mn j , -ϵ imn ζ mn j ) after the return mapping. The notation 3 G and 3
ζ are used to represent the 3rd-order gradient of the micro-deformation tensor and micromorphic generalized stress tensor respectively.
The implementation of the return mapping algorithm requires the elastic energy functional W (ε, 3 G), the yield function f θ (x, Λ), where x is the vector consisting of the components of σ and 3 ζ ; given the pre-trained plastic yield function, the plastic flow direction ∇ x f θ (x, Λ) can be derived by pytorch automatic differentiation and ∇ Λ f θ (x, Λ) = 0 in the case of perfect plasticity. Unless otherwise stated, f θ (x, Λ) is derived from the pretrained machine learning model, either from NK or MLP-LS, and W (ε, 3 G) comes from the Sobolev training of elastic energy functional [62].
Given all the necessary ingredients, Algorithm 2 first computes the elastic trial Cauchy and higher-order stresses following Eq. ( 3). If yielding is detected by f θ (x tr , Λ) ≤ 0, a system of nonlinear equations following Eq. ( 5) would be solved to find the vectorized stresses x and the increment in the plastic multiplier ∆Λ. Otherwise, the trial states are directly output.
Algorithm 2 Return mapping algorithm for higher-order continuum.
Require: Current internal variable Λ, elastic strain ε e n , elastic gradient of micro-deformation .
Matrixize
⎦ at (ϵ e,tr n+1 ,
G e,tr n+1 ).
2. Check the yield condition and perform return mapping if yield is detected.
else Solve for x and ∆Λ, such that
Extract σ n+1 and 3
Remark 3 (Stress Integration for Non-Convex Yield Surfaces). As pointed out by Lin and Bažant [63], there do exist non-convex yield surfaces in the literature, such as the Argyris yield surface (cf. Argyris et al. [64]) and the Barcelona Basic Model (cf. Sheng et al. [65]), which purposely introduce non-convexity for the sake of matching the phenomenological responses observed from experiments. As we did not enforce the convexity of the yield function explicitly, the resultant yield functions (see Figs. 11 and14) are found to have concave regions. A robust implicit stress integration may require specific treatment to find the first intersection between the non-convex yield function and an elastic trial stress point (cf. Pedroso et al. [66] and Sheng et al. [65]). In our case, we follow the treatment in GlÈ uge and Bucci [67], in which an incremental trial step smaller than the radius of the curvature of the yield surface is used such that the return mapping algorithm may yield a unique corrected state.
In this section, two examples of yield surface reconstruction with MLP-LS and NK are presented, each followed by a performance evaluation. In the first example, both methods are verified by the analytical micropolar J2 yield surface [68], followed by a short case study comparing the performance of both methods given limited or missing data. In the second example, both methods are validated by a direct numerical simulation (DNS) data upscaled from finite element simulations. For brevity, we only present the comparisons of yield surfaces obtained from MLP-LS and MK methods here. The results obtained from other alternative approaches (e.g., Gaussian Kernel, and NURBS) are presented in the Appendix.
The first example applies the MLP-LS and NK methods to reconstruct the micropolar J2 yield surface [68] with a data set inferred from the yield function in Eq. ( 19). This micropolar J2 plasticity model can be viewed as a generalized version of the classical J2 plasticity, where the additional terms with respect to the couple stress tensor m and micropolar length scale l are introduced to capture the size effect. The spatial dimension is reduced to 2D assuming plane stress, i.e. σ 33 = 0, such that only the components σ 11 , σ 22 , σ 12 , m 13 , and m 23 are non-zero in the yield function, where the mean stress p = (σ 11 + σ 22 )/3 and the deviatoric stress s = σp I.
where Y indicates the yield stress. The numerical specimen used for verification follows linear micropolar elasticity [68,69],
with bulk modulus K = 100 3 MPa, shear modulus µ = 50 MPa, coupled shear modulus µ c = 0, and yield stress Y = √ 2 3 MPa. Due to the micropolar effect, the strain tensor is divided into the hydrostatic part tr(ε)I, deviatoric part ε dev and the skew-symmetric part ε skw .
The performance of the two methods is examined by visualizing cross sections of the yield surface and the results of the return mapping algorithm. According to the test results, both methods are able to detect the yield point on the yield surface accurately and produce correct higher-order stress curves when the trained yield function is integrated into a return mapping algorithm. An additional case study with limited and missing data is conducted to evaluate the performance of both methods given low-quality data.
The yield data points generated by Algorithm 3 in Appendix B are used to train both the MLP-LS and the NK models. The MLP-LS model, constructed by a 3-layer MLP shown in Eq. ( 10) function and 64 neurons in each hidden layer, is trained with the batch size of 1000 for 500 epochs. The NK model, whose feature map φ θ is constructed by a 3-layer MLP architecture with 64 neurons in each hidden and output layer, is trained for 500 epochs. The training and validation loss histories are shown in Fig. 5, where the validation loss can be less than the training loss because the training loss function has an additional term controlling the sign of the yield function, as shown in Eq. (17).
The smooth yield surface is reconstructed by the MLP-LS and NK. The cross sections on σ 11σ 22 , σ 11σ 12 , and m 13m 23 planes are shown in Fig. 6, where both methods accurately reconstruct a yield surface that respects the ground truth yield points. Therefore, we consider that both methods are generalizable to higher-dimensional yield surfaces given a set of sufficient and well-distributed data.
In addition to the yield surface reconstruction, the return mapping algorithm 2 is implemented to verify both MLP-LS and NK methods. The material is loaded elastoplastically in a single kinematic mode, i.e. one of the uniaxial tension, shear, and bending modes, and then unloaded elastically. The test results of the return mapping, as shown in Fig. 7, indicate that both methods can be integrated into the return mapping algorithm and are able to produce valid stress curves.
It has been shown that given a set of sufficient and well-distributed yield point data, both MLP-LS and NK methods are able to reconstruct the yield surface accurately, as shown in Fig. 6. However, when the sufficient and well-distributed data are not available, the performance of the two methods is of interest. Therefore, two sets of non-ideal data are generated from the training set to test the performance of the two methods. The first data set, called the set of limited data, is sampled from the training set with a limited size via numpy.random.choice(), such that the data set is as well-distributed as the training set but has a much smaller data size. The other data set, called the set of missing data, is generated by removing a cluster of yield point data from the training set, such that the data set has a missing data patch and is considered poorly distributed.
In this case study, the limited data set with 800 data is sampled from the training set, and the missing data is generated by removing a cluster of data as shown in Fig. 8 (MIDDLE), which contains 13 665 data after the data cluster is removed. It is observed from Fig. 8 that given limited data that follows the same distribution with the test data, the yield surface can still be reconstructed accurately by both methods, shown in Fig. 8 (LEFT); however, when the missing data set is used for training, it is observed in Fig. 8 (MIDDLE and RIGHT) that NK produces a significantly more accurate prediction of the test data and stress history than MLP-LS, which reflects the robustness against the missing data.
The second example presents the application of the neural kernel method to reconstruct micropolar and micromorphic yield surfaces upscaled from direct numerical simulations. We select an idealized microstructure commonly used to represent shale, i.e. a layered material consisting of hard and soft constituents [70±72]. Readers interested in previous mathematical and neural network modeling efforts of layered geomaterials may refer to Semnani et al. [73], Zhao et al. [74], Borja et al. [75] and Xiao and Sun [35]. To the best knowledge of the authors, there has not yet been any attempt to model the micropolar and micromorphic plasticity of shale via deep learning.
The RVE data is obtained from finite element simulations on a domain consisting of layered clay materials composed of a hard and a soft constituent with intact interfaces (see Fig. 9). These layer constituents are assumed to be Cauchy continua with elasto-plastic behaviors characterized by the classical Cam-Clay model but with different material parameters. Since the constituents are assumed to exhibit no higher-order effects, the higher-order effect of the effective medium is stemmed from the spatial heterogeneity of the micro-structures. The elastic response is captured by the following elasticity energy function
where p 0 is the initial pressure, ε e v0 is the initial volumetric deformation, µ 0 indicates a constant shear modulus, and C r is the elastic re-compression ratio. The yield function with the hardening law is captured by,
where p c the preconsolidation pressure, M is the slope of the critical state line. The hardening law that governs the evolution of p c is expressed as follows:
where C c is the plastic compression index and p c0 indicates the initial preconsolidation pressure. This data set is generated from the same finite element domain used in [35]. The boundary value problem that generates the data set is solved also by the same finite element solver (cf. Xiao and Sun [35]) that employs the deal.ii library [76]. To capture the high-order constitutive behaviors, data are collected by applying the admissible boundary conditions, solving local BVP, and homogenizing the higher-order stresses according to Eqs. ( 8) and ( 9) in Section 2.
The DNS yield data are first sampled by loading the RVE at evenly parametrized deformation rate by Algorithm 4 in Appendix B, and the yield points are recorded when permanent deformation is detected.
We present the yield surfaces in a higher-dimensional stress space. The results are evaluated by comparing the accuracy and robustness of the yielding and hardening behaviors from unseen stress paths predicted by the NK and MLP-LS (which serves as the benchmark model).
To present a yield surface in a higher-dimensional stress space, the cross sections of the yield surfaces are presented by projecting the surface onto the stress planes defined by the combinations of two stress components. Two cases are studied: micropolar and micromorphic yield surfaces in higher-dimensional stress spaces. The results show that both MLP-LS and NK are able to capture the complex features of the higher-dimensional yield surface where data is sufficient. However, when data is sparsely distributed at some parts of the yield surface, NK outperforms MLP-LS in terms of extrapolating the unseen data.
Different sets of hyperparameters are used to learn the two yield surfaces with different dimensions. The micropolar yield function consists of 5 stress components, learned by the deep architecture with an input layer of 5 neurons. The MLP architecture consists of two dense hidden layers of 256 neurons and an output layer of 1 neuron, while the NK architecture consists of three dense layers of 512 neurons. The MLP-LS is trained with the batch size of 100 and the learning rate of 0.01 for 1000 epochs, while NK is trained for 500 epochs with the batch size of 200 and the learning rate of 0.0001, with the other hyperparameters being γ = 0.1, ϵ = 0.02, λ = 0.01. In the micromorphic case, the yield function consists of 9 stress components, such that the input layer with 9 neurons is used. The MLP with the same hidden dimension is trained with the batch size of 100 and the learning rate of 0.005 for 1000 epochs, while the NK architecture that consists of three dense layers of 512, 512, and 1024 neurons is trained for 100 epochs with the batch size of 500, learning rate of 0.0002, with the other hyperparameters being γ = 0.1, ϵ = 0.05, λ = 0.01.
The training and validation loss histories of different combinations of yield surfaces and models are shown in Fig. 10. The training loss is generally higher than the validation loss due to the additional term in the training loss function Eq. ( 17) that controls the sign of the yield function.
The strikes in the loss history are probably because of the mini-batch gradient-descent optimization of the loss function, where the loss is not guaranteed to be consistently decreasing, and the gradient evaluated on some data batches may be very large to create the instability of the loss history. In general, the NK has a higher training loss but a lower validation loss than MLP-LS which is more likely to overfit the data. The validation loss is defined by | f θ (x i )| given the test data x i , which does not reflect the accuracy of prediction defined by ∥x i -xi ∥, where xi is the predicted yield point ( f θ ( xi ) = 0), but the convergence of the validation loss reflects a decent accuracy of prediction.
Since the yield surface for higher-order continua depends on more than 3 variables, it is not feasible to fully visualize the learned yield surface geometrically in a three-dimensional space. As such, we projected the highdimensional yield surface onto 2D stress planes (where the rest of the stress components and internal variables are fixed) to demonstrate the geometrical features of the yield surfaces in the high-dimensional models. We then further examine the results via 2D prediction vs. ground-truth plots, samples of stress paths, and stress±strain curves for individual stress/strain components.
Micropolar yield surfaces projected onto 2D stress planes. The micropolar yield surface is visualized by the cross sections projected on 10 stress planes of different combinations of stress components, as shown in Fig. 11; the ground truth yield data are projected to the 10 stress planes as well and compared with the yield points predicted by NK and MLP-LS. We observe that both NK and MLP-LS are able to capture the complex features of the higher-dimensional yield data. Furthermore, The yield surfaces parametrized by both approaches are geometrically similar in most regions, except for a few locations where the training data is sparse, e.g., the top left corner of the yield surface projected to the first stress plane in Fig. 11. This convergence of the learned function indicates that the data is sufficiently abundant to constrain the geometry of the yield surfaces. As such, the feature extraction enabled by the data-dependent kernel of NK does not lead to significant difference in the data-rich regions.
To have a better understanding of the difference between the two predicted surfaces, we inspect the yield surface projected on the first stress plane and generate two additional unseen data as shown in Fig. 12. It is observed that although MLP-LS fits the training data more accurately, NK outperforms MLP-LS when extrapolating the missing data. The same phenomenon is observed in Fig. 8 in the previous example. The reason why NK is better at extrapolating is probably that the function space of the NK model is spanned by a controlled number of basis kernel functions as shown in Fig. 3, such that the pattern of the yield surface is recognized by the linear Galerkin projections of the data to the finite-dimensional function space; the yield surface can then be extrapolated more reasonably based on the learned pattern. Accuracy on unseen data. The unseen test yield data (sampled within the identical distribution of the train data) are compared with the yield point predictions component-wisely in Fig. 13. In this Prediction vs. Ground Truth figure, a perfect result would be a straight line with a 45-degree inclination angle. In our case, the micropolar yield surface generated by MK and MLP-LS both demonstrate sufficient accuracy, with the MLP-LS performing slightly better. This performance difference could be attributed by the different ways NK and MLP-LS parametrize the yield surface. In the NK case, the inductive bias obtained from the training is generated via the finite-dimensional space spanned by the data-dependent basis kernels which is used to express the yield surface. In the MLP-LS case, the yield surface is directly parametrized by the neural network weights. This setting is less restrictive than the kernel approach where the learned yield function must be an element of the space spanned by the kernel basis. As such, assuming that the data is sufficiently populated in this micropolar experiment, the performance gain of the MLP-LS might have been attributed to the higher expressivity of the MLP neural network [77].
Compared to the micropolar case, the dimension of the micromorphic yield surface is higher. With 9 stress components and an internal variable as the input, the yield surface would require significantly more data to populate the parametric space. Theoretically, the number of data needed to produce similar performance increases exponentially with the dimension [78], which means if 10 3 is needed to reconstruct the micropolar yield surface accurately, then approximately 10 6 data is needed in the micromorphic case in order to reach the similar performance. This data sparsity manifested by the high dimensionality might contribute to the fact that there are so few attempts to hand-craft yield surfaces for micromorphic continua and most of them are simply an extension of the Cauchy continuum counterpart. As such, many of the micromorphic and micropolar simulations are often related to F E M 2 multiscale approach that is computationally expensive [79,80]. Since it is increasingly unrealistic to always expect that we will have sufficient data when we increase the dimensionality of the model to incorporate the micromorphic effect, the robustness of the learning algorithms in the sparse data regime becomes critical.
In this micromorphic example, 16 904 data points (split into 13 500 training data and 3404 test data) are used to train and test the model. In other words, we have only increased the training data by 2.25 times (cf. Appendix B.2) whereas the data demand is expected to be increased by 3 orders. This sparsity of data is intentional, as we would like to investigate whether the feature space generated from the data-dependent kernel may enable us to generate a more robust inductive bias for extrapolation when the data is limited.
Micromorphic yield surfaces projected onto 2D stress planes. We first examine the NK and MLP-LS yield surfaces by projecting them onto 36 different stress planes. By comparing Fig. 14 with Fig. 11, this micromorphic data is distributed in a much sparser manner due to the increased dimensionality. many of the planes shown in Fig. 14, there is not even one single data point on the same plane (the orange point(s)). The sparsity showcased in these 2D planes indicates that the inductive bias inferred from the rest of the data becomes the only dominant factor dictating the prediction accuracy in those data-missing regions.
At the location regions far away from the DNS data (see ), MLP-LS method tends to generate yield surfaces of more complex shapes that vary significantly among different 2D stress planes. On the other hand, the NK yield functions (see Fig. 14) have significantly less concave regions and generally maintain a convex shape. The yield surface projected on different stress planes also exhibits more consistent geometrical patterns than those obtained from MLP-LS. This difference in the resultant yield functions is attributed to the sparsity of the data, which makes the learned function depend more significantly on the hypothesis sets employed by the NK and MLP-LS models (cf. Mohri et al. [78]).
Remarkably, in this data-limited region, the NK model seems to be capable of exploiting the structure and similarity of the data in the feature space. This exploitation on the structure and similarity of data seems to be helpful in preventing overfitting (and hence the less complex yield surface) as well as enabling the learned model to be consistent with the underlying physical laws obeyed by the data. In particular, while both the NK and MLP-LS algorithms are subjected to the DNS data set compatible to the thermodynamics principles, only the NK model yields a convex yield surface compatible with the thermodynamics constraint where data are sparse. This result is particularly interesting, because the convexity of the yield function has not been explicitly enforced via loss function or specific neural network architecture design.
Accuracy on unseen data. As shown in Fig. 15, both NK and MLP-LS are able to predict the micromorphic test data with reasonable accuracy in the unseen data. However, due to the high dimensionality and sparsity of the data, the similar accuracy in predicting a limited set of test data does not necessarily imply the similarity in the geometry of the learned yield function.
As such, the test data sampled within the same distribution of the training may not be sufficient in painting a complete picture of the prediction performance. Instead, an adversarial sampling on regions distant from the training data may provide more useful insight on the robustness of the model, as shown in numerical examples demonstrated in Figs. 8 and12.
This subsection shows the hardening process by the evolution of the yield surface and validates the NK return mapping algorithm (Algorithm 2) by the benchmark DNS stress history. To realistically reproduce the stress curve, the hardening process is modeled by introducing the magnitude of cumulative plastic strain as the internal variable, equivalent to an additional dimension of the higher-order yield surface; the initial yield surface expands as the internal variable increases. To further validate the constitutive modeling with hardening, the benchmark DNS stress history as the ground truth is first generated by controlling the deformation rate and solving the local BVP with the prescribed time-dependent boundary conditions. Another stress history is then produced by the NK return mapping algorithm and validated by the DNS benchmark. Given that the micropolar and micromorphic data could be too sparse, the DNS simulations that generate the data are run sequentially with the machine learning step described in Section 3. This adaptive strategy enables us to sample additional stress data at the locations where local support of the NK-based yield function is needed to amend the neural kernels. Note that a more rigorous active sampling strategy could potentially be derived via deep reinforcement learning, as shown in [81,82]. The rational design of As shown in Algorithm 2, the implementation of the NK-enabled return mapping algorithm requires two ingredients: the elastic energy functional and the plastic yield function that evolves with an internal variable. The elastic energy functional with respect to strain and higher-order kinematic modes is first trained by Sobolev training with the loss function of Eq. ( 24), where the predicted elastic energy and its gradients are both constrained (cf. [27]). To simulate the hardening process upon yielding, the yield function is expressed as a function of both the individual stress components and internal variable Λ is trained via the NK Algorithm 1 where Λ is the magnitude of cumulative plastic strain.
The micropolar return mapping example is first presented to introduce a test case with the microstrain being neglected. The deformation rate is controlled such that ε11 = ε22 = -6e -7, κ32 = 1.9e -4 for t < 200 and the opposite deformation rate is applied for t > 200, where t is the pseudo time. The stress history is then simulated from DNS under the elastoplastic loading and elastic unloading with the prescribed boundary condition in Eq. ( 8); all the stress and couple stress components are homogenized. The stress history is projected onto the stress plane of the combination of each stress component and m 23 , which is the component conjugate to the main kinematic mode κ 32 . The yield surfaces with different internal variables are reconstructed by NK and projected to the same stress planes, as shown in Fig. 16. It is observed that as the stress propagates, the internal variable increases, and the yield surface evolves such that the elastic region expands.
In addition to the DNS stress history as the ground truth, another stress history from the return mapping algorithm of NK is also simulated by prescribing the same deformation history. The material is first loaded elastically such that the stress increments follow the elastic energy functional W (ε, κ); when the yield is detected, the return mapping algorithm is initiated and the internal variable is incremented to model the hardening process; finally, the elastic unloading is applied, and the stress decreases with the curve showing the permanent plastic deformation. Since some stress components remain zero as shown in Fig. 16, only the nonzero components are presented in Fig. 17, and the consistency in the component-wise comparison validates the micropolar NK return mapping.
One interesting finding observed from Fig. 17 is that the bending kinematic mode can affect the Cauchy stress. Most of the handwritten models are chiral, i.e. higher-order kinematic modes and stresses being independent of the Cauchy stress and strain [51] to avoid modeling the complex coupling effect of the non-chiral material. The geomaterial we are studying is a non-chiral material according to the stress pattern shown in Fig. 18, because under the bending mode, the tension part starts to yield such that the stress stops increasing, while the compression part remains elastic with the normal stresses decreasing. As a result, the homogenized σ 11 and σ 22 decrease under the higher-order modes, and our NK method is able to capture this non-chiral effect.
In the micromorphic case, the stress history involves more components due to the additional kinematic modes. The deformation modes are controlled by the rate of ε11 = ε22 = -6e -7, Ġ222 = 3.5e -4 for t < 200 and the opposite rate is applied for t > 200. The stress history is then simulated from DNS under the elastoplastic loading and elastic unloading with the prescribed boundary condition in Eq. ( 8); all the stress and couple stress components are homogenized. The stress history is projected onto the stress plane of the combination of each stress component and ζ 222 . The yield surfaces evolving with the internal variable are reconstructed by NK and projected to the same stress planes, as shown in Fig. 19.
Similar to the micropolar case, the stress history simulated by DNS and predicted by NK return mapping under elastoplastic loading and elastic unloading, with more kinematic modes involved. The non-zero components of the stress history are compared and found to be consistent as shown in Fig. 20, which validates the micromorphic NK return mapping algorithm and shows the robustness of the NK method in modeling higher-order plasticity problems.
This paper presents the NK and MLP-LS method to reconstruct the higher-order yield surface in order to model the path-dependent constitutive law of a material with a complex micro-structure, e.g., a layered geo-material. The DNS constitutive data are first collected by solving local BVPs over the RVE domain, and then used to train the MLP-based elastic energy functional and NK-based narrow band yield function, which reproduce the path-dependent constitutive relation via return mapping algorithm. Two examples are presented to evaluate the performance of NK compared to MLP-LS. In the first example, the NK and MLP-LS are both able to reproduce a simple analytical plasticity yield model, verified by the yield surface and path-dependent stress curves. However, in the following case study, NK significantly outperforms MLP-LS in the accuracy of the yield surface and the stress curves reproduction, given the poorly distributed data with missing patches. In the second example, the two methods are validated by the micropolar and micromorphic DNS constitutive data of a layered geo-material. It is observed from the micropolar yield surface that NK outperforms MLP-LS in extrapolating unseen data, which is also observed in the first example. We consider the reason is that NK learns the pattern of the higher-dimensional surface via projecting the data onto a finite-dimensional kernel function space and extrapolates the data from the learned pattern. The micropolar return mapping results show that the material is non-chiral (i.e. Cauchy and higher-order stresses are coupled), which is also observed from the distortion of the internal RVE structure; the NK method is able to model the non-chiral effect that is hard to reflect by the handwritten models. The micromorphic results of yield surface reconstruction and return mapping further show that the NK method is generalizable to a much higher-dimensional stress space and decently reproduces the path-dependent constitutive responses given limited and missing data. the authors are supported by the National Science Foundation, United States under grant contracts CMMI-1846875 and the Dynamic Materials and Interactions Program from the Air Force Office of Scientific Research under grant contracts FA9550-21-1-0391 with the major equipment supported by FA9550-21-1-0027, and the MURI, United States of America Grant No. FA9550-19-1-0318. These supports are gratefully acknowledged. The views and conclusions contained in this document are those of the authors, and should not be interpreted as representing the official policies, either expressed or implied, of the sponsors, including the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. The views and conclusions contained in this document are those of the authors, and should not be interpreted as representing the official policies, either expressed or implied, of the sponsors, including the Army Research Laboratory, United States or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.
The second term of Eq. ( 25), named ªpart 2º, is the surface integral involving the micromorphic generalized stress. The ªpart 2º surface integral can be expanded by the divergence theorem and simplified similarly to the previous derivation.
Adding the equations Eqs. ( 27) and ( 28), the Hill's lemma Eq. ( 25) can be proven. Given the balance equation σ ji + ζ i jk,k = 0, all additional terms can be canceled other than the final result consisting of 5 terms of volume integrals, as shown in Eq. ( 29). This concludes the proof of Eq. 1
In the case where the micropolar continuum is considered, the generalized stress is replaced by the couple stress m ji = -ϵ imn ζ mn j , and the micro-deformation is replaced by the micro-rotation θ i = -1 2 ϵ i jk χ jk , where ϵ i jk is the Levi-Civita permutation symbol. As a result, the curvature tensor becomes κ i j = θ i, j = -1 2 ϵ imn G mn j , and Hill's lemma can be reduced to Eq. (30).
For the Cauchy continuum, the Hill's Lemma is further reduced to Eq. ( 31) by neglecting the micro-deformation and the higher-order generalized stress.
The admissible boundary condition of the RVE can be derived from the Hill's Lemma, which should satisfy the Hill±Mandel's condition in the form of Eq. (32), such that the left-hand side of Eq. ( 25) vanishes.
The Hill±Mandel's condition is satisfied if the 3 terms of the right-hand side of Eq. ( 25), including two surface integrals and one volume integral, become zero, as shown in (33).
One admissible boundary condition that satisfies the Hill±Mandel's condition can be derived based on the macroscopic strain tensor εi j and the gradient of micro-deformation Ḡi jk [51,52]. In the case where the base materials of RVE are Cauchy continuum (i.e. Cauchy to higher-order upscaling), χ i j = 0 internally and χi j = 0, such that ūi, j = εi j , the boundary condition shown in Eq. ( 34) can be derived, which satisfies Hill±Mandel's condition by satisfying Eq. (33).
The 3rd-order tensor G i jk is further split into the symmetric part G sym i jk = (G i jk +G jik )/2 and the skew-symmetric part G skw i jk = (G i jk -G jik )/2, such that the higher order modes can be separated into micropolar bending modes and microstrain modes. The micropolar behavior of the material can be studied by prescribing the bending modes κ i j = θ i, j = -1 2 ϵ imn G mn j = -1 2 ϵ imn G skw mn j . As a results, all the characteristic kinematic modes in 2D case are presented in Fig. 2, such that the prescribed boundary is the linear combination of characteristic
In the special case when the micropolar or Cauchy continuum is considered, the boundary condition is reduced to the form of Eq. ( 35), where Ḡi jk becomes ϵ i jl κlk for micropolar continuum and becomes zero for Cauchy continuum.
The homogenization of stress and generalized stress should still satisfy Hill±Mandel's condition in the case of Cauchy to micromorphic homogenization. In the case of Cauchy to higher-order upscaling, the ζ i jk G i jk term in Eq. ( 32) vanishes because ζ i jk vanishes in Cauchy continua and G i jk is redefined as G i jk = u i, jk . The Hill±Mandel's condition can be rewritten in the form of Eq. (36).
σ ji ε i j = σ ji εi j + ζi jk Ḡi jk (36) The left hand side of Eq. ( 36) can be rewritten into the surface integral terms by divergence theorem as shown in Eq. (37).
Based on the admissible boundary condition shown in Eq. ( 34), Eq. ( 37) is rewritten into the form with respect to εi j and Ḡi jk as shown in Eq. (38).
Subtracting equation Eq. ( 38) from Eq. ( 36), the homogenized stress σ ji and generalized stress ζi jk can be derived in Eq. (39).
For micromorphic continuum, σi j is the homogenized Cauchy stress conjugate to εi j and ζi jk is the generalized stress conjugate to Ḡi jk as computed in Eq. (39). For Cauchy continuum, only σi j is effective as G i jk is neglected. For micropolar continuum, the coupled stress mi j is derived from the skew-symmetric part of ζi jk , such that mi j = -ϵ jmn ζmni is conjugate to the curvature tensor κi j = -1 2 ϵ imn Ḡmnj . In the special case where only the symmetric part of ζi jk is considered, [83] to represent the microstrain generalized stress. quantification, but in our application, we only use the GP kernel to predict a deterministic yield surface; please refer to Williams and Rasmussen [84] for more mathematical background of GP.
The other alternative method is the NURBS-based plasticity [32] that approximates the yield function with basis spline functions, where the isotropic hardening can be modeled by the evolving control points with the internal variable [33], and the non-associative flow rule can be obtained via the gradient of the plastic potential extrapolated by Non-Uniform Rational Basic Splines (NURBS) [34].
Basis spline (B-spline) is one popular method for fitting curves or surfaces [85]. NURBS, as one kind of Bspline method, has been used to reconstruct the yield surface [32], where the NURBS model, with the control points selected in a particular way, fits the 3D ellipsoidal surface accurately. In this section, we generalize the NURBS model to surface fitting in arbitrary dimensional Euclidean spaces. We will introduce the adapted NURBS implementation from the following perspectives: (1) the basis functions, (2) parametrization of higher-dimensional yield surfaces, and (3) surface fitting with NURBS.
In general, we adopt the multi-variate B-spline basis functions as compositions of B-spline basis functions with a single variable. To this end, we first introduce the single-variate B-spline basis functions with the independent variable denoted as ϑ, which form a family of polynomials. These functions are defined piecewisely given a knot vector t that specifies a partition of interval [0, 1]: 0 = t 1 ≤ t 2 ≤ • • • ≤ t i ≤ • • • ≤ t q = 1, where q is the length of the knot vector. The single-variate basis functions B(•) of zero degree are shown to be piecewise constant functions as follows:
where the superscript with bracket indicates the degree of polynomials and the subscript i is an index indicator within the basis family. We further express the single-variate basis function of higher degrees k recursively as follows:
where we adopt k = 3, and the knot vector t = [0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1] in our implementation, which indicates that the length of the knot vector q = 11 and the number of basis function is qk -1 = 7. We denote the number of basis function members as |B|.
The composition of B i (for simplicity the superscript (k) is omitted) that formulates multi-variate B-spline basis functions is then expressed as follows:
for single variable,
The yield surface living in d-dimensional space describes a d -1-dimensional manifold and hence could be parametrized with d -1 independent variables ϑ = (ϑ 1 , ϑ 2 , . . . , ϑ d-1 ). We adopt the generalized spherical parametrization for describing the yield surface mathematically in Eq. ( 45), where the independent variables are angle indicators within the range [0, 1] and there is a dependent variable ρ(ϑ) indicating the radial distance.
j=1 sin π ϑ j ) sin 2π ϑ i if i = d -1 for higher dimensions (45)
The point cloud data x (i) sampled from the yield surface can be mapped into the labeled data (ϑ (i) , ρ (i) ) through the parametrization shown in Eq. (45). The NURBS model is described by Eq. (46), where c ∈ R |B| d-1 is an array of control points.
To fit the NURBS model with the data sampled from the yield surface, the control points c are found such that the difference between the prediction ρ(i) = N(ϑ (i) ) • c and the label ρ (i) is minimized. Therefore, the control points are computed through the regularized least square method as shown in Eq. (47), where [ρ] i = ρ (i) and [N(ϑ)] i j = N j (ϑ (i) ).
To compare the performances of the approaches outlined in Appendices C.1±C.3, we obtained the learned yield functions via these approaches for three data sets we used in this paper, i.e., the classical J2 von Mises yield function, a generalized J2 von Mises yield function with the additional couple stress terms, and the micropolar DNS data inferred from finite element simulations of the layered materials. The results are shown in Fig. 22.
In the first case, the actual learned function only depends on the J2 stress and the data set is of low dimension. As such, the results in Fig. 22(a) indicate that all four approaches (neural kernel method (red), classical Gaussian kernel (black), level set MLP (green), and NURBS (yellow)) perform well in this simple case. In the second case, the Von Mises yield function is amended with a regularization term (see Eq. ( 19)). This dependence on coupled stress also breaks the symmetry of the force stress, and hence leads to the dimension of the parametric space to be increased from 3 to 5. This increase of dimensionality leads to the NURBS yield function failing, but the rest ). In the last case, the DNS data set is used where only a portion of data is purposely missing (see Fig. 22(c)), the dimensionality of the data is identical to the second case, but the data exhibits lower symmetry. Visual inspections reveal that the neural kernel generates yield function with geometric features consistent with those of data. On the other hand, the MLP level set is capable to generate non-oscillatory yield functions, whereas both the NURBS and classical GP kernel may lead to spurious oscillation in the learned function.
Note that, in the last case where data are missing in one portion of the parametric space, the robustness of the learned function in the extrapolated regime depends strongly on the types of inductive biases of the learning algorithm [86]. Since the neural kernel method employs a data-dependent adaptive kernel, it is more expressive than the classical GP and the NURBS approach where the parametric space is spanned by a pre-determined set of bases.